G2S: The GeoStatistical Server

An extensible geostatistical framework connecting programming interfaces to simulation algorithms through a shared server.

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

Direct Sampling DS

Native Direct Sampling is available as ds, DS, or DirectSampling. It is a CPU-only implementation integrated with the modern G2S simulation foundation used by QS and AS: random or explicit simulation paths, vector and full simulation, per-node control maps, interface upload plumbing, informed-neighbor filtering, kernel handling, and OpenMP path-level parallelism.

The old DS-like implementation remains available as ds-l, dsl, DirectSamplingLike, and DS-L. Those names continue to resolve to the legacy executable and are intentionally separate from native ds.

Parameters

Usage: result = g2s(flag1, value1, flag2, value2, ...)

Schema outputs include result["simulation"] / result.simulation, result["indexmap"] / result.indexmap, result["time"] / result.time, result["progress"] / result.progress, and result["job_id"] / result.job_id.

Flag Description Mandatory
-a Use ds, DS, or DirectSampling. ✔
-ti Training image. Multiple TIs are supported. ✔
-di Destination image. NaN values are simulated; finite values are conditioning data unless --forceSimulation is set. ✔
-dt Data type per variable: 0 continuous, 1 categorical. ✔
-th Acceptance threshold. DS accepts the first scanned candidate at or below this mismatch. ✔
-f / -mer Maximum exploration ratio. TI candidates are visited without replacement through a deterministic pseudo-random permutation. ✔
-n Number of informed neighbors to use, scalar or per-variable. ✔ unless -ni is provided
-ki Optional kernel image. Kernel weights are honored using the true kernel flat index after neighborhood sorting.  
-sp Optional simulation path. Equal path priorities are resolved by index for reproducible ordering. Default is a random path.  
-j Number of OpenMP path-level workers.  
-s Global seed. DS candidate permutations and stochastic transforms use stateless hash-based draws from this seed.  
-ii Optional per-node TI-id map. Invalid values (NaN, non-finite, negative, non-integer, or out-of-range) mean unrestricted TI selection.  
-ni Optional per-node neighbor-count map.  
-kii Optional per-node kernel-index map.  
-kvi Optional per-node exploration-ratio map.  
-cti Treat each TI as circular during TI-side lookup.  
-csim Treat the simulation domain as circular for neighbor lookup.  
-wPO Enable vector simulation path optimization. Equal dependency depths are resolved by path index for same-seed reproducibility. This mirrors the QS flag and is ignored by full simulation.  
-fs Full simulation: one path entry per variable.  
--forceSimulation Simulate finite destination values too.  
-cn / -cnorm Continuous rooted Lp mismatch norm. Provide one strictly positive value for all continuous variables or one per continuous variable. Default is 2.  

Transforms

Native DS can apply stochastic local transforms to the search-pattern offsets before scanning TI candidates. The TI itself is never transformed. Transformed offsets use nearest-neighbor rounding for CPU lookup.

Flag Description
-rmi Rotation-center map. In 2D it has one angle channel in radians. In 3D it has four quaternion channels (qx,qy,qz,qw).
-rti Rotation-tolerance map with the same channel count as -rmi. Values are finite and non-negative.
-smi Isotropic scale-center map, one strictly positive finite channel.
-sti Isotropic scale tolerance, one finite non-negative channel. It is a half-width in log-scale space.

For each simulated node, DS draws a local transform around the nominal maps. Zero tolerance is deterministic. Offsets are transformed in this order: scale, rotate, nearest-neighbor rounding. 3D quaternions are normalized before use; invalid or near-zero quaternions fall back to identity rotation.

Transform cache bins are deterministic: 1 degree for 2D angles and 0.05 for log-scale. Quantized values are used both for cache keys and transformed-offset construction.

Mismatch

Categorical mismatch is a normalized weighted mismatch count. Continuous mismatch is a kernel-weighted rooted Lp mismatch using -cn / -cnorm. The common continuous powers 1 and 2 use direct arithmetic instead of generic pow calls. Mixed-variable DS uses a normalized sum across active variables and ignores non-finite neighbor or TI values without consuming support.

For each simulated node, candidate order is a deterministic pseudo-random permutation over the flattened allowed TI cells. The permutation is keyed by the global seed, local per-node seed, simulation path order, and variable, so a fixed explicit path remains reproducible while changing -s changes the TI visit order. In parallel runs, DS waits for earlier-path strict-informed neighbors before scoring candidates, so thread timing does not change the data event.

Examples

The following examples assume the G2S server is running. They are fully unconditional, create all-NaN destination images with the same shape as the loaded training image, and use -j 1.00001 for path-level parallel execution.

The Python examples load remote TIFFs with urllib and Pillow; they do not require requests. The MATLAB examples use imread with a websave fallback for older MATLAB releases. Legacy positional-output versions remain under legacy_example/ and pass -legacy_output directly.

Continuous simulation

#This code requires the G2S server to be running

from io import BytesIO
from urllib.request import urlopen

import matplotlib.pyplot as plt
import numpy as np
from PIL import Image
from g2s import g2s


def load_tiff(url):
    with urlopen(url) as response:
        return np.array(Image.open(BytesIO(response.read())), dtype=np.float32)


ti = load_tiff("https://raw.githubusercontent.com/GAIA-UNIL/TrainingImagesTIFF/master/stone.tiff")
di = np.full(ti.shape, np.nan, dtype=np.float32)

_schema_result = g2s(
    "-a",
    "ds",
    "-ti",
    ti,
    "-di",
    di,
    "-dt",
    [0],
    "-th",
    0.08,
    "-f",
    0.35,
    "-n",
    40,
    "-j",
    1.00001,
    "-s",
    123)
simulation = _schema_result["simulation"]
index = _schema_result["indexmap"]

fig, axes = plt.subplots(1, 2, figsize=(7, 3))
for ax, image, title in zip(axes, [ti, simulation], ["Stone TI", "Native DS"]):
    ax.imshow(image, cmap="gray")
    ax.set_title(title)
    ax.axis("off")
plt.tight_layout()
plt.show()


%This code requires the G2S server to be running

% This code requires the G2S server to be running.
% Load example continuous training image.
url = 'https://raw.githubusercontent.com/GAIA-UNIL/TrainingImagesTIFF/master/stone.tiff';
try
    ti = imread(url);
catch
    localTi = websave(fullfile(tempdir, 'stone.tiff'), url);
    ti = imread(localTi);
end
ti = single(ti);
di = single(nan(size(ti)));

[schemaResult, ~] = g2s('-a', 'ds', ...
                          '-ti', ti, ...
                          '-di', di, ...
                          '-dt', [0], ...
                          '-th', 0.08, ...
                          '-f', 0.35, ...
                          '-n', 40, ...
                          '-j', 1.00001, ...
                          '-s', 123);
simulation = schemaResult.simulation;
index = schemaResult.indexmap;

figure;
sgtitle('Native DS continuous unconditional simulation');
subplot(1, 2, 1);
imshow(ti, []);
title('Stone TI');
subplot(1, 2, 2);
imshow(simulation, []);
title('Native DS');

Categorical simulation

#This code requires the G2S server to be running

from io import BytesIO
from urllib.request import urlopen

import matplotlib.pyplot as plt
import numpy as np
from PIL import Image
from g2s import g2s


def load_tiff(url):
    with urlopen(url) as response:
        return np.array(Image.open(BytesIO(response.read())), dtype=np.float32)


ti = load_tiff("https://raw.githubusercontent.com/GAIA-UNIL/TrainingImagesTIFF/master/strebelle.tiff")
di = np.full(ti.shape, np.nan, dtype=np.float32)

_schema_result = g2s(
    "-a",
    "DirectSampling",
    "-ti",
    ti,
    "-di",
    di,
    "-dt",
    [1],
    "-th",
    0.12,
    "-f",
    0.4,
    "-n",
    48,
    "-j",
    1.00001,
    "-s",
    321)
simulation = _schema_result["simulation"]
index = _schema_result["indexmap"]

fig, axes = plt.subplots(1, 2, figsize=(7, 3))
for ax, image, title in zip(axes, [ti, simulation], ["Strebelle TI", "Native DS"]):
    ax.imshow(image, cmap="tab20")
    ax.set_title(title)
    ax.axis("off")
plt.tight_layout()
plt.show()


%This code requires the G2S server to be running

% This code requires the G2S server to be running.
% Load example categorical training image.
url = 'https://raw.githubusercontent.com/GAIA-UNIL/TrainingImagesTIFF/master/strebelle.tiff';
try
    ti = imread(url);
catch
    localTi = websave(fullfile(tempdir, 'strebelle.tiff'), url);
    ti = imread(localTi);
end
ti = single(ti);
di = single(nan(size(ti)));

[schemaResult, ~] = g2s('-a', 'DirectSampling', ...
                          '-ti', ti, ...
                          '-di', di, ...
                          '-dt', [1], ...
                          '-th', 0.12, ...
                          '-f', 0.4, ...
                          '-n', 48, ...
                          '-j', 1.00001, ...
                          '-s', 321);
simulation = schemaResult.simulation;
index = schemaResult.indexmap;

figure;
sgtitle('Native DS categorical unconditional simulation');
subplot(1, 2, 1);
imagesc(ti);
axis image off;
title('Strebelle TI');
subplot(1, 2, 2);
imagesc(simulation);
axis image off;
title('Native DS');
colormap(parula);

Local rotation and scale transforms

#This code requires the G2S server to be running

from io import BytesIO
from urllib.request import urlopen

import matplotlib.pyplot as plt
import numpy as np
from PIL import Image
from g2s import g2s


def load_tiff(url):
    with urlopen(url) as response:
        return np.array(Image.open(BytesIO(response.read())), dtype=np.float32)


ti = load_tiff("https://raw.githubusercontent.com/GAIA-UNIL/TrainingImagesTIFF/master/stone.tiff")
di = np.full(ti.shape, np.nan, dtype=np.float32)

yy, xx = np.indices(di.shape, dtype=np.float32)
cy = (di.shape[0] - 1) / 2.0
cx = (di.shape[1] - 1) / 2.0
radius = np.sqrt((yy - cy) ** 2 + (xx - cx) ** 2)
radius = radius / np.nanmax(radius)

rotation_map = (radius * np.pi / 3.0).astype(np.float32)
rotation_tolerance = np.full(di.shape, np.deg2rad(5.0), dtype=np.float32)
scale_map = (0.9 + 0.25 * radius).astype(np.float32)
scale_tolerance = np.full(di.shape, 0.05, dtype=np.float32)

_schema_result = g2s(
    "-a",
    "DS",
    "-ti",
    ti,
    "-di",
    di,
    "-dt",
    [0],
    "-th",
    0.08,
    "-f",
    0.35,
    "-n",
    40,
    "-j",
    1.00001,
    "-s",
    456,
    "-rmi",
    rotation_map,
    "-rti",
    rotation_tolerance,
    "-smi",
    scale_map,
    "-sti",
    scale_tolerance)
simulation = _schema_result["simulation"]
index = _schema_result["indexmap"]

fig, axes = plt.subplots(1, 4, figsize=(12, 3))
for ax, image, title, cmap in zip(
    axes,
    [ti, rotation_map, scale_map, simulation],
    ["Stone TI", "Rotation center", "Scale center", "Native DS"],
    ["gray", "twilight", "viridis", "gray"],
):
    ax.imshow(image, cmap=cmap)
    ax.set_title(title)
    ax.axis("off")
plt.tight_layout()
plt.show()


%This code requires the G2S server to be running

% This code requires the G2S server to be running.
% Load example continuous training image.
url = 'https://raw.githubusercontent.com/GAIA-UNIL/TrainingImagesTIFF/master/stone.tiff';
try
    ti = imread(url);
catch
    localTi = websave(fullfile(tempdir, 'stone.tiff'), url);
    ti = imread(localTi);
end
ti = single(ti);
di = single(nan(size(ti)));

[xx, yy] = meshgrid(single(1:size(di, 2)), single(1:size(di, 1)));
cy = single(size(di, 1) + 1) / 2;
cx = single(size(di, 2) + 1) / 2;
radius = sqrt((yy - cy).^2 + (xx - cx).^2);
radius = radius ./ max(radius(:));

rotationMap = single(radius * pi / 3);
rotationTolerance = single(ones(size(di)) * (5 * pi / 180));
scaleMap = single(0.9 + 0.25 * radius);
scaleTolerance = single(ones(size(di)) * 0.05);

[schemaResult, ~] = g2s('-a', 'DS', ...
                          '-ti', ti, ...
                          '-di', di, ...
                          '-dt', [0], ...
                          '-th', 0.08, ...
                          '-f', 0.35, ...
                          '-n', 40, ...
                          '-j', 1.00001, ...
                          '-s', 456, ...
                          '-rmi', rotationMap, ...
                          '-rti', rotationTolerance, ...
                          '-smi', scaleMap, ...
                          '-sti', scaleTolerance);
simulation = schemaResult.simulation;
index = schemaResult.indexmap;

figure;
sgtitle('Native DS transform unconditional simulation');
subplot(1, 4, 1);
imshow(ti, []);
title('Stone TI');
subplot(1, 4, 2);
imagesc(rotationMap);
axis image off;
title('Rotation center');
subplot(1, 4, 3);
imagesc(scaleMap);
axis image off;
title('Scale center');
subplot(1, 4, 4);
imshow(simulation, []);
title('Native DS');

Mixed multivariate simulation

#This code requires the G2S server to be running

from io import BytesIO
from urllib.request import urlopen

import matplotlib.pyplot as plt
import numpy as np
from PIL import Image
from g2s import g2s


def load_tiff(url):
    with urlopen(url) as response:
        return np.array(Image.open(BytesIO(response.read())), dtype=np.float32)


categorical = load_tiff("https://raw.githubusercontent.com/GAIA-UNIL/TrainingImagesTIFF/master/strebelle.tiff")
continuous = (
    categorical
    + np.roll(categorical, 1, axis=0)
    + np.roll(categorical, -1, axis=0)
    + np.roll(categorical, 1, axis=1)
    + np.roll(categorical, -1, axis=1)
) / 5.0
continuous = (continuous - np.nanmin(continuous)) / (np.nanmax(continuous) - np.nanmin(continuous))

ti = np.stack([categorical, continuous.astype(np.float32)], axis=-1).astype(np.float32)
di = np.full(ti.shape, np.nan, dtype=np.float32)

_schema_result = g2s(
    "-a",
    "ds",
    "-ti",
    ti,
    "-di",
    di,
    "-dt",
    [1, 0],
    "-th",
    0.15,
    "-f",
    0.4,
    "-n",
    48,
    "-j",
    1.00001,
    "-cnorm",
    2.0,
    "-fs",
    "-s",
    789)
simulation = _schema_result["simulation"]
index = _schema_result["indexmap"]

fig, axes = plt.subplots(1, 4, figsize=(12, 3))
images = [ti[..., 0], ti[..., 1], simulation[..., 0], simulation[..., 1]]
titles = ["Categorical TI", "Derived continuous TI", "DS category", "DS continuous"]
cmaps = ["tab20", "viridis", "tab20", "viridis"]
for ax, image, title, cmap in zip(axes, images, titles, cmaps):
    ax.imshow(image, cmap=cmap)
    ax.set_title(title)
    ax.axis("off")
plt.tight_layout()
plt.show()


%This code requires the G2S server to be running

% This code requires the G2S server to be running.
% Build a two-variable TI from the categorical Strebelle image.
url = 'https://raw.githubusercontent.com/GAIA-UNIL/TrainingImagesTIFF/master/strebelle.tiff';
try
    categorical = imread(url);
catch
    localTi = websave(fullfile(tempdir, 'strebelle.tiff'), url);
    categorical = imread(localTi);
end
categorical = single(categorical);
continuous = (categorical + circshift(categorical, 1, 1) + circshift(categorical, -1, 1) + ...
              circshift(categorical, 1, 2) + circshift(categorical, -1, 2)) / 5;
continuous = (continuous - min(continuous(:))) / (max(continuous(:)) - min(continuous(:)));

ti = cat(3, categorical, single(continuous));
di = single(nan(size(ti)));

[schemaResult, ~] = g2s('-a', 'ds', ...
                          '-ti', ti, ...
                          '-di', di, ...
                          '-dt', [1, 0], ...
                          '-th', 0.15, ...
                          '-f', 0.4, ...
                          '-n', 48, ...
                          '-j', 1.00001, ...
                          '-cnorm', 2.0, ...
                          '-fs', ...
                          '-s', 789);
simulation = schemaResult.simulation;
index = schemaResult.indexmap;

figure;
sgtitle('Native DS full mixed unconditional simulation');
subplot(1, 4, 1);
imagesc(ti(:, :, 1));
axis image off;
title('Categorical TI');
subplot(1, 4, 2);
imagesc(ti(:, :, 2));
axis image off;
title('Derived continuous TI');
subplot(1, 4, 3);
imagesc(simulation(:, :, 1));
axis image off;
title('DS category');
subplot(1, 4, 4);
imagesc(simulation(:, :, 2));
axis image off;
title('DS continuous');
colormap(parula);