2-point probability#

Theory#

This measure determines the typical distance over which two points (pixels, voxels, …) are related to each other. This is best understood by considering a binary 2D image, wherein each pixel is either black or white. This is described by the following indicator function, indicating the ‘greyscale’ of a pixel at position \(\vec{x}_i\):

\[\begin{split}\mathcal{I}(\vec{x}_i) = \begin{cases} 1 \quad & \mathrm{if}\;\; \vec{x}_i\; \in \mathrm{white} \\ 0 \quad & \mathrm{if}\;\; \vec{x}_i\; \in \mathrm{black} \\ \end{cases}\end{split}\]

The 2-point probability, \(S_2\), is the probability that two points, at a certain distance \(\Delta \vec{x}\) are both white. I.e.

\[S_2 (\Delta \vec{x}) = P \big\{ \mathcal{I}(\vec{x}) = 1 , \mathcal{I}(\vec{x} + \Delta\vec{x}) = 1 \big\}\]

Two limits can directly be identified. If \(\Delta\vec{x} = \vec{0}\), \(S_2\) is simply the probability that a point is white: the (volume) fraction of white, \(\varphi\). I.e.

\[S_2 ( || \Delta \vec{x} || = 0) = \varphi\]

The two points are completely uncorrelated if \(|| \Delta\vec{x} ||\) is sufficiently large (i.e. larger than the correlation length). In this case, both a point at \(\vec{x}\) and at \(\vec{x} + \Delta \vec{x}\) have a probability \(\varphi\) to be white, and thus

\[S_2 ( || \Delta \vec{x} || \rightarrow \infty) = \varphi^2\]

In between these extremes, \(S_2\) decays from \(\varphi\) towards the asymptotic value of \(\varphi^2\).

Example#

This example is based on a simple, periodic, image comprising circular white inclusions embedded in a black background. The figure shows from left to right: the image, the 2-point probability \(S_2\) in two dimensions, and a cross-section of this result in the middle of the region-of-interest along the horizontal axis.

Note

The Python-code can be used for the plotting: The complete code is included in the download. Note that to obtain the same plot one should download and install the matplotlib-styles available in GooseMPL.

Note

All functions make the assumption of the images being periodic. If this assumption is not reasonable be sure to specify the periodic option (that defaults True).

S2.py

import GooseEYE
import numpy as np

# generate image, extract 'volume-fraction' for plotting
img = GooseEYE.dummy_circles((500, 500))
phi = np.mean(img)

# 2-point probability
S2 = GooseEYE.S2((101, 101), img, img)

S2.cpp

#include <GooseEYE/GooseEYE.h>
#include <highfive/H5Easy.hpp>

#define MYASSERT(expr) MYASSERT_IMPL(expr, __FILE__, __LINE__)
#define MYASSERT_IMPL(expr, file, line) \
    if (!(expr)) { \
        throw std::runtime_error( \
            std::string(file) + ':' + std::to_string(line) + \
            ": assertion failed (" #expr ") \n\t"); \
    }

int main()
{
    // generate image
    auto I = GooseEYE::dummy_circles({500, 500});

    // 2-point probability
    auto S2 = GooseEYE::S2({101, 101}, I, I);

    // check against previous versions
    H5Easy::File data("S2.h5", H5Easy::File::ReadOnly);
    MYASSERT(xt::all(xt::equal(I, H5Easy::load<decltype(I)>(data, "I"))));
    MYASSERT(xt::allclose(S2, H5Easy::load<decltype(S2)>(data, "S2")));

    return 0;
}

Masked correlation#

This function (as most of GooseEYE’s functions) also has the possibility to mask certain pixels (which can be used for example to exclude acquisition artefacts from the measurement). The image’s mask is a binary matrix of exactly the same shape as the image. For each pixel in the mask with value 1, the corresponding pixel in the image is ignored. The normalisation is corrected for the reduced amount of data points, whereby the number of data points is no longer constant over the region-of-interest.

S2_mask.py,

S2_mask.cpp

Ensemble average#

To compute the ensemble average of a statistic, one constructs an Ensemble with a certain shape for the region-of-interest, and then adds the result per image to it. Consider the following example.

S2_ensemble.py

S2_ensemble.cpp

Note

An ensemble is used to compute the mean using a selection (ensemble) of relative small measurements. See Wikipedia.

import GooseEYE
import numpy as np

ensemble = GooseEYE.Ensemble((101, 101))

for i in range(5):
    img = GooseEYE.dummy_circles((200, 200))
    ensemble.S2(img, img)

S2 = ensemble.result()

#include <GooseEYE/GooseEYE.h>
#include <highfive/H5Easy.hpp>

#define MYASSERT(expr) MYASSERT_IMPL(expr, __FILE__, __LINE__)
#define MYASSERT_IMPL(expr, file, line) \
    if (!(expr)) { \
        throw std::runtime_error( \
            std::string(file) + ':' + std::to_string(line) + \
            ": assertion failed (" #expr ") \n\t"); \
    }

int main()
{
    GooseEYE::Ensemble ensemble({101, 101});

    for (size_t i = 0; i < 5; ++i) {
        auto I = GooseEYE::dummy_circles({200, 200});
        ensemble.S2(I, I);
    }

    auto S2 = ensemble.result();

    // check against previous versions
    H5Easy::File data("S2_ensemble.h5", H5Easy::File::ReadOnly);
    MYASSERT(xt::allclose(S2, H5Easy::load<decltype(S2)>(data, "S2")));

    return 0;
}

Auto-correlation#

The the greyscale generalisation of the 2-point probability (for floating-point images (with \(0 \leq \mathcal{I}(\vec{x}_i) \leq 1)\)) corresponds to:

\[S_2 (\Delta \vec{x}) = \frac{1}{N} \sum\limits_{i=1}^N \mathcal{I} (\vec{x}_i) \, \mathcal{I} (\vec{x}_i + \Delta \vec{x}) \equiv \mathcal{I} (\vec{x}) \star \mathcal{I} (\vec{x})\]

where the \(\star\) represent the convolution, in this case of \(\mathcal{I}\) with itself. Along the same arguments as for the 2-point probability, limit values can be obtained. In this case:

\[\begin{split}S_2(\Delta \vec{x} = 0) &= \langle \mathcal{I}^2 \rangle \\ S_2(\Delta \vec{x} \rightarrow \infty) &= \langle \mathcal{I} \rangle^2\end{split}\]

where the brackets \(\langle \ldots \rangle\) denotes the spatial average.

S2_autocorrelation.py

S2_autocorrelation.cpp