Phase cross correlation is not a phase cross correlation

[JulesScholler]

scikit-image/skimage/registration/_phase_cross_correlation.py

Lines 118 to 288 in 2bbfce4

	def phase_cross_correlation(reference_image, moving_image, *,
	upsample_factor=1, space="real",
	return_error=True, reference_mask=None,
	moving_mask=None, overlap_ratio=0.3):
	"""Efficient subpixel image translation registration by cross-correlation.
	This code gives the same precision as the FFT upsampled cross-correlation
	in a fraction of the computation time and with reduced memory requirements.
	It obtains an initial estimate of the cross-correlation peak by an FFT and
	then refines the shift estimation by upsampling the DFT only in a small
	neighborhood of that estimate by means of a matrix-multiply DFT.
	Parameters
	----------
	reference_image : array
	Reference image.
	moving_image : array
	Image to register. Must be same dimensionality as
	``reference_image``.
	upsample_factor : int, optional
	Upsampling factor. Images will be registered to within
	``1 / upsample_factor`` of a pixel. For example
	``upsample_factor == 20`` means the images will be registered
	within 1/20th of a pixel. Default is 1 (no upsampling).
	Not used if any of ``reference_mask`` or ``moving_mask`` is not None.
	space : string, one of "real" or "fourier", optional
	Defines how the algorithm interprets input data. "real" means
	data will be FFT'd to compute the correlation, while "fourier"
	data will bypass FFT of input data. Case insensitive. Not
	used if any of ``reference_mask`` or ``moving_mask`` is not
	None.
	return_error : bool, optional
	Returns error and phase difference if on, otherwise only
	shifts are returned. Has noeffect if any of ``reference_mask`` or
	``moving_mask`` is not None. In this case only shifts is returned.
	reference_mask : ndarray
	Boolean mask for ``reference_image``. The mask should evaluate
	to ``True`` (or 1) on valid pixels. ``reference_mask`` should
	have the same shape as ``reference_image``.
	moving_mask : ndarray or None, optional
	Boolean mask for ``moving_image``. The mask should evaluate to ``True``
	(or 1) on valid pixels. ``moving_mask`` should have the same shape
	as ``moving_image``. If ``None``, ``reference_mask`` will be used.
	overlap_ratio : float, optional
	Minimum allowed overlap ratio between images. The correlation for
	translations corresponding with an overlap ratio lower than this
	threshold will be ignored. A lower `overlap_ratio` leads to smaller
	maximum translation, while a higher `overlap_ratio` leads to greater
	robustness against spurious matches due to small overlap between
	masked images. Used only if one of ``reference_mask`` or
	``moving_mask`` is None.
	Returns
	-------
	shifts : ndarray
	Shift vector (in pixels) required to register ``moving_image``
	with ``reference_image``. Axis ordering is consistent with
	numpy (e.g. Z, Y, X)
	error : float
	Translation invariant normalized RMS error between
	``reference_image`` and ``moving_image``.
	phasediff : float
	Global phase difference between the two images (should be
	zero if images are non-negative).
	References
	----------
	.. [1] Manuel Guizar-Sicairos, Samuel T. Thurman, and James R. Fienup,
	"Efficient subpixel image registration algorithms,"
	Optics Letters 33, 156-158 (2008). :DOI:`10.1364/OL.33.000156`
	.. [2] James R. Fienup, "Invariant error metrics for image reconstruction"
	Optics Letters 36, 8352-8357 (1997). :DOI:`10.1364/AO.36.008352`
	.. [3] Dirk Padfield. Masked Object Registration in the Fourier Domain.
	IEEE Transactions on Image Processing, vol. 21(5),
	pp. 2706-2718 (2012). :DOI:`10.1109/TIP.2011.2181402`
	.. [4] D. Padfield. "Masked FFT registration". In Proc. Computer Vision and
	Pattern Recognition, pp. 2918-2925 (2010).
	:DOI:`10.1109/CVPR.2010.5540032`
	"""
	if (reference_mask is not None) or (moving_mask is not None):
	return _masked_phase_cross_correlation(reference_image, moving_image,
	reference_mask, moving_mask,
	overlap_ratio)
	# images must be the same shape
	if reference_image.shape != moving_image.shape:
	raise ValueError("images must be same shape")
	# assume complex data is already in Fourier space
	if space.lower() == 'fourier':
	src_freq = reference_image
	target_freq = moving_image
	# real data needs to be fft'd.
	elif space.lower() == 'real':
	src_freq = fftn(reference_image)
	target_freq = fftn(moving_image)
	else:
	raise ValueError('space argument must be "real" of "fourier"')
	# Whole-pixel shift - Compute cross-correlation by an IFFT
	shape = src_freq.shape
	image_product = src_freq * target_freq.conj()
	cross_correlation = ifftn(image_product)
	# Locate maximum
	maxima = np.unravel_index(np.argmax(np.abs(cross_correlation)),
	cross_correlation.shape)
	midpoints = np.array([np.fix(axis_size / 2) for axis_size in shape])
	float_dtype = image_product.real.dtype
	shifts = np.stack(maxima).astype(float_dtype, copy=False)
	shifts[shifts > midpoints] -= np.array(shape)[shifts > midpoints]
	if upsample_factor == 1:
	if return_error:
	src_amp = np.sum(np.real(src_freq * src_freq.conj()))
	src_amp /= src_freq.size
	target_amp = np.sum(np.real(target_freq * target_freq.conj()))
	target_amp /= target_freq.size
	CCmax = cross_correlation[maxima]
	# If upsampling > 1, then refine estimate with matrix multiply DFT
	else:
	# Initial shift estimate in upsampled grid
	upsample_factor = np.array(upsample_factor, dtype=float_dtype)
	shifts = np.round(shifts * upsample_factor) / upsample_factor
	upsampled_region_size = np.ceil(upsample_factor * 1.5)
	# Center of output array at dftshift + 1
	dftshift = np.fix(upsampled_region_size / 2.0)
	# Matrix multiply DFT around the current shift estimate
	sample_region_offset = dftshift - shifts*upsample_factor
	cross_correlation = _upsampled_dft(image_product.conj(),
	upsampled_region_size,
	upsample_factor,
	sample_region_offset).conj()
	# Locate maximum and map back to original pixel grid
	maxima = np.unravel_index(np.argmax(np.abs(cross_correlation)),
	cross_correlation.shape)
	CCmax = cross_correlation[maxima]
	maxima = np.stack(maxima).astype(float_dtype, copy=False)
	maxima -= dftshift
	shifts += maxima / upsample_factor
	if return_error:
	src_amp = np.sum(np.real(src_freq * src_freq.conj()))
	target_amp = np.sum(np.real(target_freq * target_freq.conj()))
	# If its only one row or column the shift along that dimension has no
	# effect. We set to zero.
	for dim in range(src_freq.ndim):
	if shape[dim] == 1:
	shifts[dim] = 0
	if return_error:
	# Redirect user to masked_phase_cross_correlation if NaNs are observed
	if np.isnan(CCmax) or np.isnan(src_amp) or np.isnan(target_amp):
	raise ValueError(
	"NaN values found, please remove NaNs from your "
	"input data or use the `reference_mask`/`moving_mask` "
	"keywords, eg: "
	"phase_cross_correlation(reference_image, moving_image, "
	"reference_mask=~np.isnan(reference_image), "
	"moving_mask=~np.isnan(moving_image))")
	return shifts, _compute_error(CCmax, src_amp, target_amp),\
	_compute_phasediff(CCmax)
	else:
	return shifts

For this function to compute the phase cross correlation the image_product should be divided by np.abs(image_product) line 220. Also the overlap_ratio is not used (so it doesn't actually work as intended) and thus very confusing.

The phase cross correlation is explained here in detail: https://en.wikipedia.org/wiki/Phase_correlation

[grlee77]

Thanks for reporting this @JulesScholler! Indeed it seems this is a long-standing bug (that appears to originate in the Matlab-central implementation mentioned in the comment at the top of the file).

Here is a simple illustration of the effect of the proposed change on a small example:

import matplotlib.pyplot as plt
import numpy as np
from scipy import fft
from skimage import data, img_as_float

src = img_as_float(data.camera()[128:256, 128:256])
src = src + 0.05 * np.random.standard_normal(src.shape)
target = np.roll(src, (15, -10), axis=(0, 1))
target = target + 0.05 * np.random.standard_normal(target.shape)
src_freq = fft.fftn(src)
target_freq = fft.fftn(target)

# current implementation
image_product1 = src_freq * target_freq.conj()
cross_correlation1 = fft.ifftn(image_product1)

# fixed implementation
image_product = image_product1 / np.abs(image_product1)
cross_correlation = fft.ifftn(image_product)

fig, axes = plt.subplots(1, 2)
axes[0].imshow(np.abs(cross_correlation1), cmap=plt.cm.gray)
axes[0].set_title('Existing Implementation')
axes[1].imshow(np.abs(cross_correlation), cmap=plt.cm.gray)
axes[1].set_title('Proposed Implementation')
for ax in axes:
    ax.set_axis_off()
plt.tight_layout()
plt.show()

In the figure below, the left image is the magnitude of the cross_correlation variable from the current implementation, whereas the right image is with the proposed change. It can be seen that the peak value is in the same location in this case, but the peak is much sharper and more isolated with the fix:

[grlee77]

Do you happen to have a test case that would succeed with this change, but currently fails without it? (Our existing test suite passes both for the existing implementation as well as with this proposed modification)

[grlee77]

With slight modification to avoid potential divide by zero:

    eps = np.finfo(image_product.real.dtype).eps
    image_product /= (np.abs(image_product) + eps)

this fix resolves the case brought up in #5460 as well

[JulesScholler]

Hello @grlee77!

Great that you added a regularization to avoid dividing by zero.

Did you manage to find a case where it would fail with the normal cross correlation and not with the phase cross correlation? Otherwise I suggest to occult part of the moving image by adding a black/white disc or square.

Let me know if you want me to test some things!

[petsuter]

Could it be this broke the RMS error return value? #7078