G2S: The GeoStatistical Server

A free and flexible multiple point (geo)statistics framework including state-of-the-art algorithms: QuickSampling and Narrow Distribution Selection

View project on GitHub Download .zip file Download .tar.gz file

QuickSampling QS

Parameters for QS

Usage: [sim,index,time,finalprogress,jobid] = g2s(flag1,value1, flag2,value2, ...) Outputs: sim = simulation, index = index of the simulated values in the flattened TI, time = simulation time, finalprogress = final progression of the simulation (normally 100), jobid = job ID.

Flag Description Mandatory
-ti Training images (one or more images). If multivariate, the last dimension should be the same size as the number of variables, and should also match the size of the array given for the parameter dt. NaN values in the training image will be ignored. Unlike other MPS-algorithms, if there are multiple variables they will not be automatically normalized to be in the same range. ✔
-di Destination image (one image, aka simulation grid). The size of di will be the size of the simulation grid. di can be identical as ti for gap-filing. NaN values will be simulated. Non-NaN values will be considered as conditioning data. ✔
-dt Data type. 0 → continuous and 1 → categorical. This is where the number of variables is defined. Provide a list containing zeros or ones representing the data type of each variable. ✔
-k Number of best candidates to consider ∈ [1 ∞]. ✔
-n N closest neighbors to consider. If multiple variables: use a single N value if identical for all variables or use an array of N values if each variable has a different N. ✔
-ki Image of the weighting kernel. Can be used to normalize the variables. If multiple variables: use a single kernel value if identical for all variables or if each variable has a different kernel, stack all kernels along the last dimension. The number of kernels should then match the size of the array given for the parameter dt. Possibility to use different kernels for different parts of the simulation with the “-kii” parameter explained below.  
-f f=1/k equivalent to f of DS with a threshold to 0.  
-j To run in parallel (default is single core). To use as follows ‘-j’, N1, N2, N3 (all three are optional but N3 needs N2, which in turn needs N1). Use integer values to specify a number of threads (or logical cores). Use decimal values ∈]0,1[ to indicate fraction of the maximum number of available logical cores (e.g., 0.5=50% of all available logical cores). N1 threads used to parallelize the path (path-level) Default: the maximum number of threads available. N2 threads used to parallelize over training images (if many TIs are available, each is scanned on a different core). Default: 1. N3 threads used to parallelize FFTs (node-level). Default: 1. N1 and N2 are recommended over N3. N1 is usually more efficient than N2, but requires more memory.  
-sp Simulation path, array of the size of the simulation grid containing values that specify the simulation path (from low to high). Default is a random path. Equal values are accepted but will always be ordered in the same way (i.e. not random). -∞ values are not simulated. In case of multiple variables, a vector simulation is default (same path for all variables) and the simulation path should be one dimension less than the number of variables. If you prefer a full simulation, provide an array containing the path for each variable and use the “-fs” flag below.  
-s Random seed value.  
-W_GPU Use integrated GPU if available.  
-W_CUDA Use Nvidia Cuda compatible GPU: specify the device id. (e.g., '-W_CUDA',0,1 to use the first two GPUs).  

Advanced parameters

Flag Description Mandatory
-ii Array that specifies for each pixel which training image to sample from. Default: all training images are searched for the best match.  
-kii Array that specifies for each pixel which kernel to use. Useful if you want different weight fields when simulating different parts of the di, or if you want certain covariables to only influence parts of your di (e.g. you want a surface variable to only influence the simulation of the top subsurface layer). Usage: provide a list of full kernels in -ki, and only a single image for -kii to specify which full kernel to use at which pixel.  
-ni Array that specifies for each pixel which the number of neighbors.  
-kvi Array that specifies for each pixel which the number of best candidates to consider.  
-far Fast and risky 😄, like -ii but with a random input (experimental).  
-cti With this flag QS will treat the training image(s) as periodic (aka circular or cyclic) over each dimension.  
-csim With this flag QS will make sure to create a periodic (aka circular or cyclic) simulation over each dimension.  
-adsim Augmented dimentionality simulation: allows for 3D simulation using 2D training image, only for categories (Coming maybe some day!).  
-fs Full simulation: follows a different simulation path for each variable (as opposed to vector simulation, where the same simulation path is used for all variables).  
-nV No Verbatim, i.e. prohibits neighbors in the training image to be neighbors in the simulation. (experimental).  
--forceSimulation Restimulate a value even if already existing.  

Examples

Below are several examples showcasing different applications of QS. For these examples the G2S server should be installed and running, either on your own machine or remotely. A Google Colab notebook with more examples and an automatic installation of G2S can be found here.

Unconditional simulation

#This code requires the G2S server to be running

import numpy
from PIL import Image
import requests
from io import BytesIO
from g2s import g2s
import matplotlib.pyplot as plt

# load example training image ('stone')
ti = numpy.array(Image.open(BytesIO(requests.get(
    'https://raw.githubusercontent.com/GAIA-UNIL/TrainingImagesTIFF/master/stone.tiff').content)));

# QS call using G2S
simulation,*_=g2s('-a','qs',
                 '-ti',ti,
                 '-di',numpy.zeros((200,200))*numpy.nan,
                 '-dt',[0], #Zero for continuous variables
                 '-k',1.2,
                 '-n',50,
                 '-j',0.5);

#Display results 
fig, (ax1, ax2) = plt.subplots(1, 2,figsize=(7,4))
fig.suptitle('QS Unconditional simulation',size='xx-large')
ax1.imshow(ti)
ax1.set_title('Training image');
ax1.axis('off');
ax2.imshow(simulation)
ax2.set_title('Simulation');
ax2.axis('off');
plt.show()


%This code requires the G2S server to be running

%load data
ti=imread('https://raw.githubusercontent.com/GAIA-UNIL/TrainingImagesTIFF/master/stone.tiff');

% QS call using G2S
simulation=g2s('-a','qs','-ti',ti,'-di',nan(200,200),'-dt',[0],'-k',1.2,'-n',50,'-j',0.5);

%Display results 
sgtitle('Unconditional simulation');
subplot(1,2,1);
imshow(ti);
title('Training image');
subplot(1,2,2);
imshow(simulation);
title('Simulation');

Conditional simulation

#This code requires the G2S server to be running

import numpy
from PIL import Image
import requests
from io import BytesIO
from g2s import g2s
import matplotlib.pyplot as plt

# load example training image ('stone')
ti = numpy.array(Image.open(BytesIO(requests.get(
    'https://raw.githubusercontent.com/GAIA-UNIL/TrainingImagesTIFF/master/stone.tiff').content)));

# create empty grid and add conditioning data
conditioning = numpy.zeros((200,200))*numpy.nan; 
# fill the grid with 50 random points
conditioning.flat[numpy.random.permutation(conditioning.size)[:50]]=ti.flat[numpy.random.permutation(ti.size)[:50]];

# QS call using G2S
simulation,*_=g2s('-a','qs', 
                 '-ti',ti,
                 '-di',conditioning,
                 '-dt',[0],
                 '-k',1.2,
                 '-n',30,
                 '-j',0.5);

# Display results 
fig, (ax1, ax2, ax3) = plt.subplots(1, 3,figsize=(7,4),subplot_kw={'aspect':'equal'})
fig.suptitle('QS Conditional simulation',size='xx-large',y=0.9)
ax1.imshow(ti)
ax1.set_title('Training image');
ax1.axis('off');
ax2.scatter(*numpy.meshgrid(numpy.arange(200),numpy.arange(200,0,-1)),s=5,c=conditioning,marker='.')
ax2.set_title('Conditioning');
ax2.axis('off');
ax3.imshow(simulation)
ax3.set_title('Simulation');
ax3.axis('off');
plt.show()


%This code requires the G2S server to be running

%load data
ti=imread('https://raw.githubusercontent.com/GAIA-UNIL/TrainingImagesTIFF/master/stone.tiff');
% empty grid 
conditioning=nan(200,200); 
% fill the grid with 50 random points
conditioning(randperm(numel(conditioning),50))=ti(randperm(numel(ti),50));

% QS call using G2S
simulation=g2s('-a','qs','-ti',ti,'-di',conditioning,'-dt',[0],'-k',1.2,'-n',50,'-j',0.5);

%Display results 
sgtitle('Conditional simulation');
subplot(1,3,1);
imshow(ti);
title('Training image');
subplot(1,3,2);
imshow(conditioning);
title('Conditioning');
subplot(1,3,3);
imshow(simulation);
title('Simulation');

Simulation with multiple Training Images

#This code requires the G2S server to be running

import numpy
from PIL import Image
import requests
from io import BytesIO
from g2s import g2s
import matplotlib.pyplot as plt

#ti1 contains horizontal lines and ti2 vertical lines
ti1 = numpy.tile(numpy.sin(range(150)),(150,1))
ti2 = ti1.transpose()

# QS call using only the horizontal lines as TI
simulation1,index1,*_ = g2s('-a','qs',
                         '-ti',ti1,
                         '-di',numpy.zeros((150,150))*numpy.nan,
                         '-dt',[1], #1 for categorical variable
                         '-k',1.2,
                         '-n',25,
                         '-j',0.5);

# QS call using both horizontal and vertical lines as TI's
simulation2,index2,*_ = g2s('-a','qs',
                         '-ti',[ti1,ti2],
                         '-di',numpy.zeros((150,150))*numpy.nan,
                         '-dt',[1],
                         '-k',1.2,
                         '-n',25,
                         '-j',0.5);

#Display results
fig, ([[ax1,ax2],[ax3,ax4]]) = plt.subplots(2,2,figsize=(10,10),sharex = True,sharey = True)
fig.suptitle('QS Multiple TI simulation',size='xx-large')
ax1.imshow(ti1)
ax1.set_title('Training image 1');
ax1.axis('off');
ax2.imshow(ti2)
ax2.set_title('Training image 2');
ax2.axis('off');
ax3.imshow(simulation1)
ax3.set_title('Simulation 1');
ax3.axis('off');
ax4.imshow(simulation2)
ax4.set_title('Simulation 2');
ax4.axis('off');


%This code requires the G2S server to be running

%ti1 contains horizontal lines and ti2 vertical lines
ti=repmat((sin(1:150)>0)',1,150);
ti1=ti;
ti2=rot90(ti,1);

%QS call using only the horizontal lines as TI
simulation1=g2s('-a','qs', ...
            '-di',nan(150,150,1), ...
            '-ti',{ti1}, ...
            '-dt',[1], ...
            '-k',1.2, ...
            '-n',25, ...
            '-j',0.5);

% QS call using both horizontal and vertical lines as TI's
simulation2=g2s('-a','qs', ...
            '-di',nan(150,150,1), ...
            '-ti',{ti1 ti2}, ...
            '-dt',[1], ...
            '-k',1.2, ...
            '-n',25, ...
            '-j',0.5);

%Display results 
figure(4);clf
sgtitle('multi-TI');
subplot(2,2,1);
imshow(ti1);
title('Training image 1');
subplot(2,2,2);
imshow(ti2);
title('Training image 2');
subplot(2,2,3);
imshow(simulation1);
title('Simulation with training image 1');
subplot(2,2,4);
imshow(simulation2);
title('Simulation with training images 1 and 2');

Multivariate simulation

#This code requires the G2S server to be running

import numpy
import requests
from io import BytesIO
from g2s import g2s
import matplotlib.pyplot as plt
import tifffile as tf

ti = tf.imread(BytesIO(requests.get(
    'https://raw.githubusercontent.com/GAIA-UNIL/TrainingImagesTIFF/master/ti_3_variables.tiff').content))

###Example 1 - Complete multivariate simulation
#Here all three variables are entirely unkwown and will be simulated simultaneously

# QS call using G2S (with dt set to two continuous and one categorical variable)
simulation1,index1,*_=g2s('-a','qs',
                 '-ti',ti,
                 '-di',numpy.zeros([400,400,3])*numpy.nan,
                 '-dt',[0,0,1],
                 '-k',1.2,
                 '-n',50,
                 '-j',0.5);

#Display results 
fig, ([[ax1,ax2,ax3],[ax4,ax5,ax6],[ax7,ax8,ax9]]) = plt.subplots(3,3,figsize=(15,10),sharex = True,sharey = True)
fig.suptitle('QS Multivariate simulation - example 2',size='xx-large')
ax1.imshow(ti[:,:,0])
ax1.set_title('Training image dim 1');
ax1.axis('off');
ax2.imshow(ti[:,:,1])
ax2.set_title('Training image dim 2');
ax2.axis('off');
ax3.imshow(ti[:,:,2])
ax3.set_title('Training image dim 3');
ax3.axis('off');
ax4.imshow(numpy.zeros([400,400])*numpy.nan)
ax4.set_title('Destination image dim 1')
ax4.axis('off')
ax5.imshow(numpy.zeros([400,400])*numpy.nan)
ax5.set_title('Destination image dim 2')
ax5.set_aspect('equal')
ax5.axis('off')
ax6.imshow(numpy.zeros([400,400])*numpy.nan)
ax6.set_title('Destination image dim 3')
ax6.axis('off')
ax7.imshow(simulation1[:,:,0])
ax7.set_title('Simulation dim 1');
ax7.axis('off');
ax8.imshow(simulation1[:,:,1])
ax8.set_title('Simulation dim 2');
ax8.axis('off');
ax9.imshow(simulation1[:,:,2])
ax9.set_title('Simulation dim 3');
ax9.axis('off');



### Example 2 - Multivariate simulation with partially informed covariables
#In many situations we will have fully and/or partially informed covariables (e.g. topography + point measurements)
#To simulate a different environment than the training images, we will use the results from example 1 above
#Variable 2 is partially informed and variable 3 is fully informed

#Take random values from the previous simulation of variable 2 and add them to the di
di_var2 = numpy.zeros([400,400])*numpy.nan
di_var2.flat[numpy.random.permutation(di_var2.size)[:500]] = ti[:,:,1].flat[numpy.random.permutation(ti[:,:,1].size)[:500]]
di_example2 = numpy.stack([numpy.zeros([400,400])*numpy.nan,
                                 di_var2,
                                 simulation1[:,:,2]],axis = 2)

# QS call using G2S (with dt set to two continuous and one categorical variable)
simulation2,index,*_=g2s('-a','qs',
                 '-ti',ti,
                 '-di',di_example2,
                 '-dt',[0,0,1],
                 '-k',1.2,
                 '-n',50,
                 '-j',0.5);

#Display results 
fig, ([[ax1,ax2,ax3],[ax4,ax5,ax6],[ax7,ax8,ax9]]) = plt.subplots(3,3,figsize=(15,10),sharex = True,sharey = True)
fig.suptitle('QS Multivariate simulation - example 2',size='xx-large')
ax1.imshow(ti[:,:,0])
ax1.set_title('Training image dim 1');
ax1.axis('off');
ax2.imshow(ti[:,:,1])
ax2.set_title('Training image dim 2');
ax2.axis('off');
ax3.imshow(ti[:,:,2])
ax3.set_title('Training image dim 3');
ax3.axis('off');
ax4.imshow(di_example2[:,:,0])
ax4.set_title('Destination image dim 1')
ax4.axis('off')
ax5.scatter(*numpy.meshgrid(numpy.arange(400),numpy.arange(400,0,-1)),s=5,c=di_example2[:,:,1],marker='.')
ax5.set_title('Destination image dim 2')
ax5.set_aspect('equal')
ax5.axis('off')
ax6.imshow(di_example2[:,:,2])
ax6.set_title('Destination image dim 3')
ax6.axis('off')
ax7.imshow(simulation2[:,:,0])
ax7.set_title('Simulation dim 1');
ax7.axis('off');
ax8.imshow(simulation2[:,:,1])
ax8.set_title('Simulation dim 2');
ax8.axis('off');
ax9.imshow(simulation2[:,:,2])
ax9.set_title('Simulation dim 3');
ax9.axis('off');



%This code requires the G2S server to be running

%load data
ti=imread('https://raw.githubusercontent.com/GAIA-UNIL/TrainingImagesTIFF/master/ti_3_variables.tiff');

% QS call using G2S (with dt set to two continuous and one categorical
% variable)
simulation=g2s('-a','qs','-ti',ti,'-di',nan(400,400,3),'-dt',[0,0,1],'-k',1.2,'-n',30,'-j',0.5);

%Display results 
sgtitle('Unconditional simulation');
subplot(2,3,1);
imshow(ti(:,:,1));
title('Training image dim 1');
subplot(2,3,2);
imshow(ti(:,:,2));
title('Training image dim 2');
subplot(2,3,3);
imshow(ti(:,:,3));
title('Training image dim 3');
subplot(2,3,4);
imshow(simulation(:,:,1));
title('Simulation dim 1');
subplot(2,3,5);
imshow(simulation(:,:,2));
title('Simulation dim 2');
subplot(2,3,6);
imshow(simulation(:,:,3));
title('Simulation dim 3');
colormap winter;

Gap filling

#This code requires the G2S server to be running

import numpy
from PIL import Image
import requests
from io import BytesIO
from g2s import g2s
import matplotlib.pyplot as plt

# load example training image ('stone') and cut out a part of it
tiWithGap = numpy.array(Image.open(BytesIO(requests.get(
    'https://raw.githubusercontent.com/GAIA-UNIL/TrainingImagesTIFF/master/stone.tiff').content)));
tiWithGap[60:140,60:140]=numpy.nan;

# QS call using G2S
simulation,*_=g2s('-a','qs',
                 '-ti',tiWithGap,
                 '-di',tiWithGap,
                 '-dt',[0],
                 '-k',1.2,
                 '-n',25,
                 '-j',0.5);

#Display results 
fig, (ax1, ax2) = plt.subplots(1, 2,figsize=(7,4))
fig.suptitle('QS Gap filling',size='xx-large')
ax1.imshow(tiWithGap)
ax1.set_title('Training image');
ax1.axis('off');
ax2.imshow(simulation)
ax2.set_title('Simulation');
ax2.axis('off');
plt.show()


%This code requires the G2S server to be running

%load data
tiWithGap=imread('https://raw.githubusercontent.com/GAIA-UNIL/TrainingImagesTIFF/master/stone.tiff');
tiWithGap(60:140,60:140)=nan;

% QS call using G2S
simulation=g2s('-a','qs','-ti',tiWithGap,'-di',tiWithGap,'-dt',[0],'-k',1.2,'-n',25,'-j',0.5);

%Display results 
sgtitle('Gap filling');
subplot(1,2,1);
imshow(tiWithGap);
title('Training image');
subplot(1,2,2);
imshow(simulation);
title('Simulation');

Downscaling

#This code requires the G2S server to be running

import numpy as np
from PIL import Image
from g2s import g2s
import matplotlib.pyplot as plt
from matplotlib import colors
from io import BytesIO
import requests

#Load example training image "Dunes gobi"
ti_full = Image.open(BytesIO(requests.get(
    'https://raw.githubusercontent.com/GAIA-UNIL/TrainingImagesTIFF/master/dunes_gobi.tiff').content))
#Keep one version at fine resolution
ti_fine = np.array(ti_full)
#Make a synthetic coarse resolution image by upscaling and reprojecting to the original grid 
ti_coarse = np.array(ti_full.resize(
    (int(ti_full.width/5),int(ti_full.height/5))).resize(
    (ti_full.width,ti_full.height),resample=Image.NEAREST))

#Display the full training image at both resolutions
norm =colors.Normalize(vmin=ti_fine.min(),vmax=ti_fine.max())
f,(ax1,ax2) = plt.subplots(2,1,figsize=(12,12))
ax1.imshow(ti_fine,norm=norm)
ax1.set_title('Fine resolution training image')
ax2.imshow(ti_coarse, norm =norm)
ax2.set_title('Coarse resolution training image')

#Crop half of the image to be used as ti
ti = np.stack((ti_fine[:500,:500],ti_coarse[:500,:500]),axis=2)
#Crop upper right corner to be used as di
di_coarse = ti_coarse[:200,-200:]
di_fine   = np.zeros((200,200))*np.nan
di=np.stack((di_fine,di_coarse),axis=2)
#dt consists of two zeros representing two continuous variables
dt = [0]*ti.shape[-1]

# QS call using G2S
simulation,index,*_=g2s('-a','qs',
                         '-ti',ti,
                         '-di',di,
                         '-dt',dt,
                         '-k',1.0,
                         '-n',30,
                         '-j',0.5)

#Display results 
norm =colors.Normalize(vmin=ti_fine.min(),vmax=ti_fine.max())
f,((ax1,ax2),(ax3,ax4)) = plt.subplots(2,2,figsize=(12,12),sharex=True,sharey=True)
plt.subplots_adjust(wspace=0.1,hspace=0.1)
plt.suptitle('QS Downscaling',size='xx-large')
ax1.imshow(di_coarse,norm=norm)
ax1.set_title('Coarse res di')
ax2.imshow(simulation[:,:,0], norm =norm)
ax2.set_title('Simulation')
ax3.imshow(index)
ax3.set_title('Index')
ax4.imshow(ti_fine[:200,-200:],norm=norm)
ax4.set_title('True image')
ax4.set_xticks(np.arange(0,200,25),np.arange(ti_full.width-200,ti_full.width,25));


%This code requires the G2S server to be running

%load example training image 'Dunes gobi'
ti_fine=imread('https://raw.githubusercontent.com/GAIA-UNIL/TrainingImagesTIFF/master/dunes_gobi.tiff');

%create an artificial coarse resolution image on the same grid 
ti_coarse = imresize(imresize(ti_fine,0.2,'nearest'),5,'nearest');

%display the full training image at both resolutions
figure(1);
tiledlayout(2,1);
nexttile;
image(ti_fine);
title('Fine TI');
nexttile;
image(ti_coarse);
title('Coarse TI');

ti_size = 500;
di_size = 200;

%crop half of the image, to be used as ti
ti = cat(3,ti_fine(1:ti_size,1:ti_size),ti_coarse(1:ti_size,1:ti_size));
%crop upper right corner to be used as di
di_coarse = ti_coarse(1:di_size,(end-di_size+1):end);
di_fine  = nan(di_size,di_size);
di = cat(3,di_fine,di_coarse);


% QS call using G2S
[simulation,index,time]=g2s('-a','qs',...
                '-ti',ti,...
                '-di',di,...
                '-dt',[0 0],... #Zero for continuous variables
                '-k',1.0,...
                '-n',30,...
                '-j',0.5);

figure(2);
tiledlayout(2,2)
nexttile;
image(di_coarse);
title('Coarse DI');
nexttile;
image(simulation(:,:,1));
title('Simulation');
nexttile;
image(index(:,:,1),'CDataMapping','scaled');
title('Index');
nexttile;
image(ti_fine(1:di_size,(end-di_size+1):end));
title('True image');

3D simulation

#This code requires the G2S server to be running

import numpy as np
import requests
from io import BytesIO
from g2s import g2s
import matplotlib.pyplot as plt
import tifffile as tf
from matplotlib import cm

#Load example boolean training image "Damier 3D" and crop it to reduce computational requirements
dim = 30
ti = tf.imread(BytesIO(requests.get(
    'https://raw.githubusercontent.com/GAIA-UNIL/TrainingImagesTIFF/master/Damier3D.tiff').content))[:dim,:dim,:dim]

# QS call using G2S
simulation,index,time,*_=g2s('-a','qs',
                             '-ti',ti,
                             '-di',np.zeros_like(ti)*np.nan,
                             '-dt',[1], #1 for categorical variables
                             '-k',1.2,
                             '-n',30)

#Display results
fig = plt.figure(figsize=(20,15))
ax1=fig.add_subplot(121,projection='3d')
ax1.voxels(ti,alpha=0.9,facecolor='tab:blue',edgecolor='black')
ax1.view_init(azim=45)
ax1.set_title('3D Training image')
ax2=fig.add_subplot(122,projection='3d')
ax2.voxels(simulation,alpha=0.9,facecolor='tab:blue',edgecolor='black')
ax2.view_init(azim=45)
ax2.set_title('QS 3D simulation')

#Display the indices of the ti-values found in the di-image
viridis = cm.get_cmap('viridis',index.max())
colormap = viridis(index)
fig = plt.figure(figsize=(10,10))
ax1=fig.add_subplot(111,projection='3d')
ax1.voxels(index,facecolors=colormap)
ax1.view_init(azim=45);


%This code requires the G2S server to be running

%% 3D multivariate

clear;home;clf
%download TI (first variable, continuous)
ti1=imread('https://raw.githubusercontent.com/GAIA-UNIL/TrainingImagesTIFF/master/fold_continuous.tiff');
%create second variable (categorical) by thresholding the first
ti2=ti1>0.5;
ti=cat(4,ti1,ti2);
ti=double(ti(1:80,1:50,101:120,:));

%display TI
subplot(2,2,1)
[X,Y,Z] = meshgrid(1:size(ti,2),1:size(ti,1),1:size(ti,3));
slice(X,Y,Z,ti(:,:,:,1),[1 size(ti,2)],[1 size(ti,1)],[1 size(ti,3)])
shading flat
camproj('perspective')
axis equal
title('TI variable 1');

subplot(2,2,2)
[X,Y,Z] = meshgrid(1:size(ti,2),1:size(ti,1),1:size(ti,3));
slice(X,Y,Z,ti(:,:,:,2),[1 size(ti,2)],[1 size(ti,1)],[1 size(ti,3)])
shading flat
camproj('perspective')
axis equal
title('TI variable 2');

% QS call using G2S
tic
simulation=g2s('-a','qs','-ti',ti,'-di',nan(80,50,10,2),'-dt',[0,1],'-k',1.2,'-n',{20 10},'-j',0.5);
toc

%display simulation
subplot(2,2,3)
[X,Y,Z] = meshgrid(1:size(simulation,2),1:size(simulation,1),1:size(simulation,3));
slice(X,Y,Z,simulation(:,:,:,1),[1 size(simulation,2)],[1 size(simulation,1)],[1 size(simulation,3)])
shading flat
camproj('perspective')
axis equal
title('Simulation variable 1');

subplot(2,2,4)
[X,Y,Z] = meshgrid(1:size(simulation,2),1:size(simulation,1),1:size(simulation,3));
slice(X,Y,Z,simulation(:,:,:,2),[1 size(simulation,2)],[1 size(simulation,1)],[1 size(simulation,3)])
shading flat
camproj('perspective')
axis equal
title('Simulation variable 2');
sgtitle('3D Simulation with two variables');

Asynchronous mode

#This code requires the G2S server to be running

import numpy
from PIL import Image
import requests
from io import BytesIO
from g2s import g2s
import matplotlib.pyplot as plt
import time

# load example training image ('stone')
ti = numpy.array(Image.open(BytesIO(requests.get(
    'https://raw.githubusercontent.com/GAIA-UNIL/TrainingImagesTIFF/master/stone.tiff').content)));


# asynchronous QS call using G2S with the "-submitOnly" flag
jobid_1=g2s('-a','qs', 
                 '-submitOnly',
                 '-ti',ti,
                 '-di',numpy.zeros((200,200))*numpy.nan,
                 '-dt',[0],
                 '-k',1.2,
                 '-n',50,
                 '-j',0.5);

# second asynchronous call that waits for job 1 to finish using the "-after" flag
jobid_2 = g2s('-a','qs', 
                  '-after',jobid_1,
                  '-submitOnly',
                  '-ti',ti,
                  '-di',numpy.zeros((200,200))*numpy.nan,
                  '-dt',[0],
                  '-k',1.2,
                  '-n',50,
                  '-j',0.5);

# check the status of both jobs in 2-second intervals using the "-statusOnly" flag
status_1 = 0
status_2 = 0
while status_2 < 95:
    time.sleep(2)
    status_1,*_ = g2s('-statusOnly',jobid_1)
    status_2,*_ = g2s('-statusOnly',jobid_2)
    print('Status jobs 1 & 2:   ', status_1,status_2)

# retrieve the simulation results from the server using the "-waitAndDownload" flag
# if the simulation would not be ready yet this call would wait for it to be ready
sim1,*_ = g2s('-waitAndDownload',jobid_1)
sim2,*_ = g2s('-waitAndDownload',jobid_2)

# display results 
fig, (ax1, ax2,ax3) = plt.subplots(1, 3,figsize=(7,4))
fig.suptitle('QS Unconditional simulation',size='xx-large')
ax1.imshow(ti)
ax1.set_title('Training image');
ax1.axis('off');
ax2.imshow(sim1)
ax2.set_title('Simulation 1');
ax2.axis('off');
ax3.imshow(sim2)
ax3.set_title('Simulation 2');
ax3.axis('off');
plt.show()


%This code requires the G2S server to be running

% load example traning image ('stone')
ti=imread('https://raw.githubusercontent.com/GAIA-UNIL/TrainingImagesTIFF/master/stone.tiff');

% asynchronous QS call using G2S with the "-submitOnly" flag
jobid_1=g2s('-a','qs','-submitOnly','-ti',ti,'-di',nan(200,200),'-dt',[0],'-k',1.2,'-n',50,'-j',0.5);

% 2nd async call that waits for job 1 to finish using the "-after" flag
jobid_2=g2s('-a','qs','-after',jobid_1,'-submitOnly','-ti',ti,'-di',nan(200,200),'-dt',[0],'-k',1.2,'-n',50,'-j',0.5);

% check the status of both jobs in 2-second intervals using the
% "-statusOnly" flag
status_1 = 0;
status_2 = 0;
while status_2 < 95
    pause(2);
    status_1=g2s('-statusOnly',jobid_1);
    status_2=g2s('-statusOnly',jobid_2);
    fprintf('Status jobs 1 & 2:   %s %s\n', status_1, status_2);
end

% retrieve the simulation results from the server using the
%"-waitAndDownload" flag. If the simulation would not be ready yet this
%call would wait for it to be ready
sim1=g2s('-waitAndDownload',jobid_1);
sim2=g2s('-waitAndDownload',jobid_2);

% display results
sgtitle('Unconditional asynchronous simulation');
subplot(1,3,1);
imshow(ti);
title('Training image');
subplot(1,3,2);
imshow(sim1);
title('Simulation 1');
subplot(1,3,3);
imshow(sim2);
title('Simulation 2');

Publication

Gravey, M., & Mariethoz, G. (2020). QuickSampling v1.0: a robust and simplified pixel-based multiple-point simulation approach. Geoscientific Model Development, 13(6), 2611–2630. https://doi.org/10.5194/gmd-13-2611-2020

Benchmarking

This code CAN be used for benchmarking (and I invite you to do so 😉), the code needs to run natively on macOS or on Linux using the Intel Compiler with MKL library. The version needs to be reported, and the time needs to be the time reported by the algorithm (that is the true execution time without taking into account interfaces overhead).

When benchmarking, the code should NOT be used inside a Virtual Machine or through WSL on Windows 10+.

Advanced

This page provides some additional information on how to do common tasks. It also serves as personal note to myself 😁

Multiple realization

A question that is frequently asked is how to do multiple realizations. Currently there are three different ways to do it, each one of them has pros and cons.

The lazy solution

The simplest solution is to do a for-loop over the realizations. However, in each step the algorithm needs to wait for the data to load, and this creates an overhead

sims1=numpy.empty((250,250,numberOfSimulation));
for i in range(numberOfSimulation):
  sims1[:,:,i],*_=g2s('-sa',computationServer,'-a','qs','-ti',ti,'-di',numpy.zeros((250,250))*numpy.nan,'-dt',numpy.ones((1,)),'-k',1.2,'-n',50,'-j',0.5,'-s',100+i);

For-loop without overhead

By using the G2S submission queue, it’s possible to remove most of the overhead. This is the most versatile solution and it is recommended over the first solution, especially in case the computations are run on another machine. Furthermore, in this solution the parameters can be changed for each realization.

ds=numpy.empty((numberOfSimulation,), dtype='long');
for i in range(numberOfSimulation):
  ids[i]=g2s('-sa',computationServer,'-a','qs','-ti',ti,'-di',numpy.zeros((250,250))*numpy.nan,'-dt',numpy.ones((1,)),'-k',1.2,'-n',50,'-j',0.5,'-s',100+i,'-submitOnly');

sims2=numpy.empty((250,250,numberOfSimulation));
for i in range(numberOfSimulation):
  sims2[:,:,i],*_=g2s('-sa',computationServer,'-waitAndDownload',ids[i]);

Overdimensioning of the destination image

This is the third solution to get multiple simulations at once. Although this solution looks easier, it has more limitations (e.g., being limited same size, same parameters, etc.), and therefore it is ⚠️ not guaranteed to stay functional in the future. This last solution should only be considered in case of extremely parallelized simulations (number of threads >50) and/or extremely small simulations (less than 1e4 pixels per simulation)

sims3,*_=g2s('-sa',computationServer,'-a','qs','-ti',ti,'-di',numpy.zeros((numberOfSimulation,250,250))*numpy.nan,'-dt',numpy.ones((1,)),'-k',1.2,'-n',50,'-j',0.5);

How to do a simulation by segment

Although QS tends to reduce issues related to non-stationarity, it won’t remove them completely. Therefore in cases with high non-stationarity, the challenge is to respect each zone separately. Here we assume that we have different training images for each specific zone.

Using a secondary variable

The first solution for doing a simulation by segment is to add a secondary variable with the information of non-stationarity. In this case the setting of the algorithm has to be identical for both training images.

ti1=imread('https://raw.githubusercontent.com/GAIA-UNIL/TrainingImagesTIFF/master/bangladesh.tiff');
ti2=imread('https://raw.githubusercontent.com/GAIA-UNIL/TrainingImagesTIFF/master/strebelle.tiff');

dt=[0];
N=300;

ti1_bis=cat(3,ti1,zeros(size(ti1)));
ti2_bis=cat(3,ti2,ones(size(ti2)));

localKernel=cat(3,kernel,padarray(15,25*[1,1]));
% localKernel=cat(3,kernel,kernel);

sim=g2s('-a','qs','-ti',ti1_bis,ti2_bis,'-di',cat(3,nan(size(tiSelection)),tiSelection),'-dt',[0,1],'-n',50,'-k',1.5,'-j',0.5,'-ki',localKernel);
imshow(medfilt2(sim(:,:,1)))
drawnow;
pause(5)

Using training image indices

The second solution for doing a simulation by segment is to use a training image index map. This requires that each training image represents a stationary subdomain. Also in this case the setting of the algorithm needs to be identical between training images.


ti1=imread('https://raw.githubusercontent.com/GAIA-UNIL/TrainingImagesTIFF/master/bangladesh.tiff');
ti2=imread('https://raw.githubusercontent.com/GAIA-UNIL/TrainingImagesTIFF/master/strebelle.tiff');

dt=[0];
N=300;

sim=g2s('-a','qs','-ti',ti1,ti2,'-di',nan(size(tiSelection)),'-dt',dt,'-n',50,'-k',1.5,'-j',0.5,'-ki',kernel,'-ii',tiSelection);
imshow(medfilt2(sim))
pause(5)

The optimal solution

Using the -ni and -kii … TODO!!

Simulating each subdomain sequentially

The last solution for doing a simulation by segment is to do simulations per subdomain. This allows for changing the settings between each subdomain. However each subdomain is simulated sequentially and therefore can be at the origin of some artifacts, regarding the order of simulations.

ti1=imread('https://raw.githubusercontent.com/GAIA-UNIL/TrainingImagesTIFF/master/bangladesh.tiff');
ti2=imread('https://raw.githubusercontent.com/GAIA-UNIL/TrainingImagesTIFF/master/strebelle.tiff');

dt=[0];
N=300;

path=reshape(randperm(numel(tiSelection)),size(tiSelection));
p1=path;
p1(tiSelection==0)=-inf;
p2=path;
p2(tiSelection==1)=-inf;
sim=g2s('-a','qs','-ti',ti1,'-di',nan(size(tiSelection)),'-dt',dt,'-n',50,'-k',1.5,'-j',0.5,'-ki',kernel,'-sp',p1);
sim=g2s('-a','qs','-ti',ti2,'-di',sim,'-dt',dt,'-n',50,'-k',1.5,'-j',0.5,'-ki',kernel,'-sp',p2);
imshow(medfilt2(sim))
pause(5)