GEONE - DEESSE - Incomplete images

Main points addressed

• reconstruction of an incomplete image (i.e. having uninformed cells)

Import what is required

import numpy as np
import matplotlib.pyplot as plt
import time
import os

# import package 'geone'
import geone as gn

# Show version of python and version of geone
import sys
print(sys.version_info)
print('geone version: ' + gn.__version__)

sys.version_info(major=3, minor=13, micro=7, releaselevel='final', serial=0)
geone version: 1.3.0

Incomplete (training) image

Read the incomplete image from the file ti_inc.txt, precising the value used in the file as missing data (uninformed cell).

# Read file
data_dir = 'data'
filename = os.path.join(data_dir, 'ti_inc.txt')
im = gn.img.readImageTxt(filename, missing_value=-9999999)

# Values
im.get_unique()

array([0., 1., 2.])

# Setting for categories / colors
categ_val = [0, 1, 2]
categ_col = ['lightblue', 'darkgreen', 'orange']

plt.figure(figsize=(5,5))
gn.imgplot.drawImage2D(im, categ=True, categVal=categ_val, categCol=categ_col, title='Incomplete image')
plt.show()

# Set lines delineating the missing part
ind = np.where(np.isnan(im.val))
iy0, iy1 = ind[2].min(), ind[2].max()
ix0, ix1 = ind[3].min(), ind[3].max()
x0 =  im.x()[ix0] - 0.5*im.sx
x1 =  im.x()[ix1] + 0.5*im.sx
y0 =  im.y()[iy0] - 0.5*im.sy
y1 =  im.y()[iy1] + 0.5*im.sy
line_x = [x0, x1, x1, x0, x0]
line_y = [y0, y0, y1, y1, y0]

plt.figure(figsize=(5,5))
gn.imgplot.drawImage2D(im, categ=True, categVal=categ_val, categCol=categ_col, title='Incomplete image')
plt.plot(line_x, line_y, c='red', ls='dashed')
plt.show()

Reconstruction of the image

Deesse can be used to reconstruct the missing part of the image, i.e. to simulate the variable in uninformed cells. For that, the same incomplete image is specified as the training image and as hard data set. In this situation, the image is called a training data set.

Note: of course, an incomplete training image can be used for simulation on a different simulation grid with or without hard data.

Fill the input structure for deesse and launch deesse

nreal = 20
deesse_input = gn.deesseinterface.DeesseInput(
    nx=im.nx, ny=im.ny, nz=im.nz, # set same dimensions grid as for the incomplete image
    sx=im.sx, sy=im.sy, sz=im.sz, # set same cell units grid as for the incomplete image
    ox=im.ox, oy=im.oy, oz=im.oz, # set same origin grid as for the incomplete image
    nv=1, varname=im.varname,     # number of variable(s), name of the variable(s),
                                  #    as for the incomplete image
    TI=im,                        # set the incomplete image as TI
    dataImage=im,                 # set the incomplete image as hard data
    distanceType='categorical',
    nneighboringNode=24,
    distanceThreshold=0.02,
    maxScanFraction=0.25,
    npostProcessingPathMax=1,
    seed=444,
    nrealization=nreal)

# Run deesse
t1 = time.time() # start time
deesse_output = gn.deesseinterface.deesseRun(deesse_input)
t2 = time.time() # end time
print(f'Elapsed time: {t2-t1:.2g} sec')

deesseRun: DeeSse running... [VERSION 3.2 / BUILD NUMBER 20230914 / OpenMP 19 thread(s)]
deesseRun: DeeSse run complete
deesseRun: warnings encountered (1 times in all):
#  1: WARNING 00010: a variable in a training image is not exhaustively informed
Elapsed time: 9.4 sec

Retrieve the reconstructed images, do some statistics (and display)

# Retrieve the results
sim = deesse_output['sim']

# Do some statistics on the realizations
# compute the pixel-wise proportion for the given categories
all_sim_stats = gn.img.imageListCategProp(sim, categ_val)

# Color map for proportion of each category
prop_col=['blue', 'darkgreen', 'orange'] # colors for the proportion maps
prop_cmap = [gn.customcolors.custom_cmap(['white', c]) for c in prop_col]

# Display
plt.subplots(2, 3, figsize=(17,8), sharex=True, sharey=True) # 1 x 3 sub-plots
for i in range(3):
    plt.subplot(2, 3, i+1) # select next sub-plot
    gn.imgplot.drawImage2D(sim[i], categ=True, categVal=categ_val, categCol=categ_col)
    plt.plot(line_x, line_y, c='red', ls='dashed') # add line delineating the reconstructed part
    plt.title(f'Real. #{i}')

for i in range(3):
    plt.subplot(2, 3, i+4) # select next sub-plot
    gn.imgplot.drawImage2D(all_sim_stats, iv=i, cmap=prop_cmap[i],
                           title=f'Prop. of categ. {i} (over {nreal} real.)')
    plt.plot(line_x, line_y, c='black', ls='dashed') # add line delineating the reconstructed part
plt.show()

Reconstruction using pyramids

Enabling pyramids implies multi-resolution simulations, which can help to better reproduce the spatial structures. Below, the pyramid type is set to categorical_auto: one pyramid for the indicator variable of each category except one is built and used.

See jupyter notebook ex_deesse_16_advanced_use_of_pyramids.ipynb for more details about pyramids.

# Deesse input
#   - with 2 additional levels (npyramidLevel=2)
#     and a reduction factor of 2 along x and y axes between the original image
#     and the 1st pyramid level and between the 1st pyramid level to the second one
#     (kx=[2, 2], ky=[2, 2], kz=[0, 0]: do not apply reduction along z axis)
pyrGenParams = gn.deesseinterface.PyramidGeneralParameters(npyramidLevel=2, kx=[2, 2], ky=[2, 2], kz=[0, 0])
pyrParams = gn.deesseinterface.PyramidParameters(nlevel=2, pyramidType='categorical_auto')

deesse_input = gn.deesseinterface.DeesseInput(
    nx=im.nx, ny=im.ny, nz=im.nz,
    sx=im.sx, sy=im.sy, sz=im.sz,
    ox=im.ox, oy=im.oy, oz=im.oz,
    nv=1, varname=im.varname,
    TI=im,
    dataImage=im,
    distanceType='categorical',
    nneighboringNode=24,
    distanceThreshold=0.02,
    maxScanFraction=0.25,
    pyramidGeneralParameters=pyrGenParams, # pyramid general parameters
    pyramidParameters=pyrParams,           # pyramid parameters for each variable
    npostProcessingPathMax=1,
    seed=444,
    nrealization=nreal)

# Run deesse
t1 = time.time() # start time
deesse_output = gn.deesseinterface.deesseRun(deesse_input)
t2 = time.time() # end time
print(f'Elapsed time: {t2-t1:.2g} sec')

deesseRun: DeeSse running... [VERSION 3.2 / BUILD NUMBER 20230914 / OpenMP 19 thread(s)]
deesseRun: DeeSse run complete
deesseRun: warnings encountered (3 times in all):
#  1: WARNING 00010: a variable in a training image is not exhaustively informed
#  2: WARNING 00210: pryramid: a variable in a training image is not exhaustively informed (in an additional level)
Elapsed time: 16 sec

# Retrieve the results
sim = deesse_output['sim']

# Do some statistics on the realizations
# compute the pixel-wise proportion for the given categories
all_sim_stats = gn.img.imageListCategProp(sim, categ_val)

# Color map for proportion of each category
prop_col=['blue', 'darkgreen', 'orange'] # colors for the proportion maps
prop_cmap = [gn.customcolors.custom_cmap(['white', c]) for c in prop_col]

# Display
plt.subplots(2, 3, figsize=(17,8), sharex=True, sharey=True) # 1 x 3 sub-plots
for i in range(3):
    plt.subplot(2, 3, i+1) # select next sub-plot
    gn.imgplot.drawImage2D(sim[i], categ=True, categVal=categ_val, categCol=categ_col)
    plt.plot(line_x, line_y, c='red', ls='dashed') # add line delineating the reconstructed part
    plt.title(f'Real. #{i}')

for i in range(3):
    plt.subplot(2, 3, i+4) # select next sub-plot
    gn.imgplot.drawImage2D(all_sim_stats, iv=i, cmap=prop_cmap[i],
                           title=f'Prop. of categ. {i} (over {nreal} real.)')
    plt.plot(line_x, line_y, c='black', ls='dashed') # add line delineating the reconstructed part
plt.show()

Reconstruction using pyramids (2)

Here the pyramid type is set to categorical_to_continuous: the pyramid is built for one continuous variable derived from the categorical variable (accounting for the connections between adjacent nodes).

# Deesse input
pyrParams = gn.deesseinterface.PyramidParameters(nlevel=2, pyramidType='categorical_to_continuous')

deesse_input = gn.deesseinterface.DeesseInput(
    nx=im.nx, ny=im.ny, nz=im.nz,
    sx=im.sx, sy=im.sy, sz=im.sz,
    ox=im.ox, oy=im.oy, oz=im.oz,
    nv=1, varname=im.varname,
    TI=im,
    dataImage=im,
    distanceType='categorical',
    nneighboringNode=24,
    distanceThreshold=0.02,
    maxScanFraction=0.25,
    pyramidGeneralParameters=pyrGenParams, # pyramid general parameters
    pyramidParameters=pyrParams,           # pyramid parameters for each variable
    npostProcessingPathMax=1,
    seed=444,
    nrealization=nreal)

# Run deesse
t1 = time.time() # start time
deesse_output = gn.deesseinterface.deesseRun(deesse_input)
t2 = time.time() # end time
print(f'Elapsed time: {t2-t1:.2g} sec')

deesseRun: DeeSse running... [VERSION 3.2 / BUILD NUMBER 20230914 / OpenMP 19 thread(s)]
deesseRun: DeeSse run complete
deesseRun: warnings encountered (3 times in all):
#  1: WARNING 00010: a variable in a training image is not exhaustively informed
#  2: WARNING 00210: pryramid: a variable in a training image is not exhaustively informed (in an additional level)
Elapsed time: 8.1 sec

# Retrieve the results
sim = deesse_output['sim']

# Do some statistics on the realizations
# compute the pixel-wise proportion for the given categories
all_sim_stats = gn.img.imageListCategProp(sim, categ_val)

# Color map for proportion of each category
prop_col=['blue', 'darkgreen', 'orange'] # colors for the proportion maps
prop_cmap = [gn.customcolors.custom_cmap(['white', c]) for c in prop_col]

# Display
plt.subplots(2, 3, figsize=(17,8), sharex=True, sharey=True) # 1 x 3 sub-plots
for i in range(3):
    plt.subplot(2, 3, i+1) # select next sub-plot
    gn.imgplot.drawImage2D(sim[i], categ=True, categVal=categ_val, categCol=categ_col)
    plt.plot(line_x, line_y, c='red', ls='dashed') # add line delineating the reconstructed part
    plt.title(f'Real. #{i}')

for i in range(3):
    plt.subplot(2, 3, i+4) # select next sub-plot
    gn.imgplot.drawImage2D(all_sim_stats, iv=i, cmap=prop_cmap[i],
                           title=f'Prop. of categ. {i} (over {nreal} real.)')
    plt.plot(line_x, line_y, c='black', ls='dashed') # add line delineating the reconstructed part
plt.show()