[Python] Why doesn’t the barycenter method detect subpixel displacements where correlation...

I’m working with X-ray imaging data. I have a reference image containing a structured pattern, and a sample image where this pattern is slightly distorted due to the presence of a physical sample.

My goal is to track local distortions of the pattern with subpixel accuracy.

Using local cross-correlation methods, I’m able to detect the expected subpixel displacements. However, when I use a method based on the local barycenter (center of mass) of intensity values in a sliding window, it only seems to detect shifts near the edges of the sample. Inside the sample — despite the fact that the pattern is visibly deformed — the computed displacement is nearly zero.

The barycenter formula I’m using is:

bx(x, y) = \frac{\sum I(i, j) \cdot (i - c)}{\sum I(i, j)} by(x, y) = \frac{\sum I(i, j) \cdot (j - c)}{\sum I(i, j)}

where (i, j) iterate over a window of size w × w centered at (x, y), and c = (w - 1) / 2.

The displacement is then computed as:

dx = bx_{sample} - bx_{reference}
dy = by_{sample} - by_{reference}

I also tried applying a Gaussian weight in the barycenter computation to give more importance to the central pixels of each window, but this did not change the outcome.

Additionally, I applied the same method to the gradients of the image (∇x and ∇y) instead of the raw intensity values, hoping to enhance sensitivity to local shifts — but again, the result remains the same: no displacement is detected inside the sample region.

I’m trying to understand why this barycenter-based approach fails to detect motion inside the deformed region, and whether there’s a way to adapt or improve it to make it a viable alternative to correlation-based methods.

Any insight or suggestions would be greatly appreciated!

import numpy as np
from numpy.lib.stride_tricks import sliding_window_view

def barycentre(img, window_size=9, iterations=100):
# Padding to handle borders and perform calculation on central values of the image
pad_size = window_size // 2
padded_array = np.pad(img, pad_size, mode='edge')

# Hanscomb transform
img = 2 * np.sqrt(img + (3 / 8))

windows = sliding_window_view(padded_array, (window_size, window_size)).astype(np.float32)

X, Y = np.mgrid[0:window_size, 0:window_size] - pad_size

# Optionally apply a Gaussian weight (commented out)
# G = np.exp(-(X**2 + Y**2) / (2 * (pad_size/2)**2))
# G = G / G.sum()
# windows *= G[None, None, ...]

Ix = np.sum(windows * X[None, None, ...], axis=(-2, -1))
Iy = np.sum(windows * Y[None, None, ...], axis=(-2, -1))
I = np.sum(windows, axis=(-2, -1)) + 1e-8

barycentre_x = Ix / I
barycentre_y = Iy / I

return barycentre_x, barycentre_y

# Example usage:
barycentre_x, barycentre_y = barycentre(img, window_size=3, iterations=100)
barycentre_x_ref, barycentre_y_ref = barycentre(ref, window_size=3, iterations=100)
dx = barycentre_x - barycentre_x_ref
dy = barycentre_y - barycentre_y_ref

Continue reading...