Weighted correlation#

Theory#

The weighted correlation characterises the average indicator \(\mathcal{I}\) around high weight factor \(\mathcal{W}\).

Mathematically the weighted correlation reads

\[\mathcal{P} (\Delta \vec{x}) = \frac{ \sum_i \mathcal{W}(\vec{x}_i) \; \mathcal{I}(\vec{x}_i + \Delta x) }{ \sum_i \mathcal{W}(\vec{x}_i) \; }\]

Additionally pixels can be masked, for instance to ignore \(\mathcal{I}\) everywhere where \(\mathcal{W}\) is non-zero. The masked correlation reads

\[\mathcal{P} (\Delta \vec{x}) = \frac{ \sum_{i}\; \mathcal{W} (\vec{x}_i) \; [ \mathcal{I} (1-\mathcal{M}) ] (\vec{x}_i + \Delta \vec{x}) \; }{ \sum_{i}\; \mathcal{W} (\vec{x}_i) \; (1-\mathcal{M})\, (\vec{x}_i + \Delta \vec{x}) \; }\]

where all pixels where \(\mathcal{M}(\vec{x}_i) = 1\) are ignored; all pixels for which \(\mathcal{M}(\vec{x}_i) = 0\) are considered as normal.

Example#

W2.py

W2.cpp

Note

Like for the 2-point correlation, a mask can be used. Similarly, the average can be extended to that of an ensemble of images.

Python#

import GooseEYE
import numpy as np
import prrng

# image + "damage" + correlation
# ------------------------------

# square grid of circles
N = 15
M = 500
row = np.linspace(0, M, N)
col = np.linspace(0, M, N)
row, col = np.meshgrid(row, col, indexing="ij")  # ('indexing' only for comparison with C++ code)
row = row.reshape(-1)
col = col.reshape(-1)
r = float(M) / float(N) / 4.0 * np.ones(N * N)

# random perturbation
rng = prrng.pcg32(0)
row += rng.normal(shape=[N * N], mu=0, sigma=M / N)
col += rng.normal(shape=[N * N], mu=0, sigma=M / N)
r *= rng.normal(shape=[N * N], mu=0, sigma=2) + 0.1

# generate image, extract 'volume-fraction' for plotting
img = GooseEYE.dummy_circles((M, M), np.round(row), np.round(col), np.round(r))
phi = np.mean(img)

# create 'damage' -> right of inclusion
col += 1.1 * r
r *= 0.4
W = GooseEYE.dummy_circles((M, M), np.round(row), np.round(col), np.round(r))
W[img == 1] = 0

# weighted correlation
WI = GooseEYE.W2((101, 101), W, img, fmask=W)

# gray-scale image + correlation
# ------------------------------

# convert to gray-scale image and introduce noise
Igr = np.array(img, copy=True).astype(float)
Igr += 0.1 * (2.0 * rng.random(Igr.shape) - 1.0) + 0.1
Igr /= 1.2

# mean intensity (for bounds)
Iav = np.mean(Igr)

# weighted correlation
WIgr = GooseEYE.W2((101, 101), W.astype(float), Igr, fmask=W)

C++#

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

#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()
{
    // image + "damage" + correlation
    // ------------------------------

    // square grid of circles
    size_t N = 15;
    size_t M = 500;
    auto row = xt::linspace<double>(0, M, N);
    auto col = xt::linspace<double>(0, M, N);
    xt::xarray<double> rowmat;
    xt::xarray<double> colmat;
    std::tie(rowmat, colmat) = xt::meshgrid(row, col);
    rowmat = xt::ravel(rowmat);
    colmat = xt::ravel(colmat);
    xt::xtensor<double, 1> r = (double)(M) / (double)(N) / 4.0 * xt::ones<double>({N * N});

    // random perturbation
    prrng::pcg32 rng(0);
    rowmat += rng.normal({N * N}, 0.0, (double)(M) / (double)(N));
    colmat += rng.normal({N * N}, 0.0, (double)(M) / (double)(N));
    auto dr = rng.normal({N * N}, 0.0, 2.0) + 0.1;
    r = r * dr;

    // generate image
    auto I = GooseEYE::dummy_circles({M, M}, xt::round(rowmat), xt::round(colmat), xt::round(r));

    // create 'damage' -> right of inclusion
    colmat += 1.1 * r;
    r *= 0.4;
    auto W = GooseEYE::dummy_circles({M, M}, xt::round(rowmat), xt::round(colmat), xt::round(r));
    W = xt::where(xt::equal(I, 1), 0, W);

    // weighted correlation
    auto WI = GooseEYE::W2({101, 101}, W, I, W);

    // gray-scale image + correlation
    // ------------------------------

    // convert to gray-scale image and introduce noise
    xt::xarray<double> Igr = I;
    Igr += 0.1 * (2.0 * rng.random(Igr.shape()) - 1.0) + 0.1;
    Igr /= 1.2;

    // weighted correlation
    auto WIgr = GooseEYE::W2({101, 101}, xt::xarray<double>(W), Igr, W);

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

    return 0;
}

Collapse to single point#

To calculate the probability of the inclusion directly next to a weight site (i.e. the red circles in the example above and below) the ‘collapsed correlation’ is calculated. The distance to the edge of the site, \(\vec{\delta}_i\) is therefore corrected for as follows:

\[\mathcal{P} (\Delta \vec{x}) = \frac{ \sum_{i}\; \mathcal{W} (\vec{x}_i) \; \mathcal{I} (\vec{x}_i + \Delta \vec{x} + \vec{\delta}_i) \; }{ \sum_{i}\; \mathcal{W} (\vec{x}_i) \; }\]

Similarly to the above, a mask may be introduced as follows:

\[\mathcal{P} (\Delta \vec{x}) = \frac{ \sum_{i}\; \mathcal{W} (\vec{x}_i) \; [ \mathcal{I} (1-\mathcal{M}) ] (\vec{x}_i + \Delta \vec{x} + \vec{\delta}_i) \; }{ \sum_{i}\; \mathcal{W} (\vec{x}_i) \; (1-\mathcal{M})\, (\vec{x}_i + \Delta \vec{x} + \vec{\delta}_i) \; }\]