Subpixel interpolation errors

Image correlation can easily be more accurate than pixel in displacements. However for measuring very small displacements, systematic interpolation errors can appear that have a small bias towards integer or 0.5 pixel displacements.

This example illustrates and studies this systematic error.

This error is typically reduced using higher order interpolation (order 3 instead of order 1).

We will also show an elegant mitigation from [Wantz2025] called “Shift-DVC” that we have implemented in 2D and 3D.

We will apply synthetic displacements to an existing pattern, and it is essential to use a different interpolator than what will be used for the correlation, for this we will Fourier shifts.

Import modules

import numpy
import matplotlib.pyplot as plt
import spam.deformation
import spam.DIC
import spam.datasets

Define function to shift images in Fourier space

def fourier_shift(im, shift_z=0, shift_y=0, shift_x=0):
    nz, ny, nx = im.shape

    # FFT freq
    kz = numpy.fft.fftfreq(nz)
    ky = numpy.fft.fftfreq(ny)
    kx = numpy.fft.fftfreq(nx)

    # put them in the freq grid
    kz, ky, kx = numpy.meshgrid(kz, ky, kx, indexing="ij")

    # phase shift
    phase = numpy.exp(-2j * numpy.pi * (shift_z * kz + shift_y * ky + shift_x * kx))

    # returb shifted
    return numpy.real(numpy.fft.ifftn(numpy.fft.fftn(im) * phase))

Here we will load the data and synthetically apply subpixel displacements from [0, 1] px in N steps

# Load data
im = spam.datasets.loadSnow()[0:50, 0:50, 0:50]

N = 20
steps = numpy.linspace(0, 1, N + 1)

# N x im series of displaced images
ims = numpy.zeros((len(steps), *im.shape))

for n, step in enumerate(steps):
    ims[n] = fourier_shift(im, shift_x=step)

Here we will do image correlations, varying the shift option and the interpolation order

for shift in [True, False]:
    for order in [1, 3]:
        displacements = numpy.zeros((N + 1, 3))
        for i in range(N + 1):
            reg = spam.DIC.register(
                ims[0],
                ims[i],
                # margin=4,
                interpolationOrder=order,
                # deltaPhiMin=0.001,
                shift=shift,
                # PhiInit = spam.deformation.computePhi({'t': [0.,0.,-1.]})
                # verbose=True
            )
            displacements[i] = reg["Phi"][0:3, -1]
        print(displacements[:, 2])
        plt.plot(steps, displacements[:, 2] - steps, ".-", label=f"{order = }, {shift = }")
plt.xlabel("Applied displacement (px)")
plt.ylabel("Error (px)\n Measured displacement - Applied Displacement")
plt.legend()
plt.show()

[0.    0.05  0.1   0.151 0.202 0.254 0.305 0.354 0.403 0.452 0.5   0.548
 0.597 0.646 0.696 0.746 0.798 0.849 0.9   0.95  1.   ]
[0.    0.05  0.101 0.151 0.201 0.251 0.301 0.351 0.401 0.451 0.5   0.55
 0.6   0.649 0.699 0.749 0.799 0.849 0.899 0.95  1.   ]
[0.    0.059 0.115 0.168 0.219 0.268 0.316 0.363 0.408 0.454 0.499 0.544
 0.589 0.635 0.682 0.73  0.779 0.831 0.884 0.941 1.   ]
[0.    0.051 0.103 0.153 0.204 0.254 0.304 0.353 0.402 0.451 0.5   0.548
 0.597 0.646 0.696 0.746 0.796 0.846 0.897 0.949 1.   ]