Reconstruct microstructure from pymks two point correlation data

Please excuse the delayed response, I did not receive a notification of the opened issue.

The easiest approach would be to use the TwoPointCorrelation descriptor of MCRpy, which is a differentiable approximatio to the TwoPointCorrelation in pyMKS. Or even better, you could use the three-point correlations, which are called Correlations in MCRpy.

Alternatively, if you really need to use the version from pymks, then you can convert the two. The conversion formula is given in this work https://www.researchgate.net/publication/352036973_Reconstructing_random_heterogeneous_media_through_differentiable_optimization in Equation 27. Note that pyMKS provides all symmetrically equivalent entries of S2, which MCRpy does not. Also, the entries are sorted in a different way. An example implementation for this conversion is given here:

import numpy as np

from pymks import TwoPointCorrelation
import mcrpy

cutoff = 25
limit_to = cutoff + 1
ms = np.load('mcrpy/example_microstructures/pymks_ms_64x64_2.npy') # this example is included in the MCRpy GitHub repo

pymks_correlations = TwoPointCorrelation(correlations=[(1, 1)], cutoff=cutoff, periodic_boundary=True).transform(
        np.expand_dims(np.stack([1-ms, ms], axis=-1), axis=0)) # note that we compute only the autocorrelation of phase 1

settings = mcrpy.CharacterizationSettings(
    descriptor_types=['TwoPointCorrelations'], limit_to=limit_to, use_multigrid_descriptor=False, use_multiphase=False)
mcrpy_correlations = mcrpy.characterize(ms, settings)['TwoPointCorrelations']

def compute_factor(settings: mcrpy.CharacterizationSettings) -> float:
    zl = 1 / (1 + np.exp(0.75 * settings.threshold_steepness))
    zu = 1 / (1 + np.exp(-(0.25 * settings.threshold_steepness) ))
    a = 1 / (zu - zl)
    b = -a * zl
    sigmoid = lambda z: 1/(1 + np.exp(-z))
    return sigmoid(-0.25 * settings.threshold_steepness) * a + b

def correct_thresholding(pcs: np.array, settings: mcrpy.CharacterizationSettings) -> np.array:
    ft05 = compute_factor(settings) # Equation 27
    factor = 2 * ft05 ** 2
    pcs[1:] = (1 - factor) * pcs[1:] + factor * pcs[0]
    return pcs

def sort_pymks_to_mcrpy(pymks_s2: np.ndarray, limit_to: int) -> np.ndarray:
    assert pymks_s2.shape[0] == 1
    assert pymks_s2.shape[-1] == 1
    assert len(pymks_s2.shape) == 4
    pymks_s2 = pymks_s2[0, ..., 0]

    out = []
    for i in range(limit_to):
        for j in range(limit_to):
            out.append(pymks_s2[cutoff + i, cutoff + j])
    for i in range(1, limit_to):
        for j in range(1, limit_to):
            out.append(pymks_s2[cutoff + i, cutoff - j])
    return np.array(out)

def convert_pymks_to_mcrpy(pymks_correlations: np.ndarray, settings: mcrpy.CharacterizationSettings) -> np.ndarray:
    return correct_thresholding(sort_pymks_to_mcrpy(pymks_correlations, settings.limit_to), settings)

def error(a, b):
    return np.max(np.abs(a - b))

print(error(convert_pymks_to_mcrpy(pymks_correlations, settings), mcrpy_correlations)) # should be around 1E-9

If you have any further questions, please send an e-mail to paul.seibert@tu-dresden.de

NEFM-TUDresden / MCRpy

Reconstruct microstructure from pymks two point correlation data #1