1""" 2===================================== 3Image Registration 4===================================== 5 6In this example, we use phase cross-correlation to identify the 7relative shift between two similar-sized images. 8 9The ``phase_cross_correlation`` function uses cross-correlation in 10Fourier space, optionally employing an upsampled matrix-multiplication 11DFT to achieve arbitrary subpixel precision [1]_. 12 13.. [1] Manuel Guizar-Sicairos, Samuel T. Thurman, and James R. Fienup, 14 "Efficient subpixel image registration algorithms," Optics Letters 33, 15 156-158 (2008). :DOI:`10.1364/OL.33.000156` 16 17""" 18import numpy as np 19import matplotlib.pyplot as plt 20 21from skimage import data 22from skimage.registration import phase_cross_correlation 23from skimage.registration._phase_cross_correlation import _upsampled_dft 24from scipy.ndimage import fourier_shift 25 26image = data.camera() 27shift = (-22.4, 13.32) 28# The shift corresponds to the pixel offset relative to the reference image 29offset_image = fourier_shift(np.fft.fftn(image), shift) 30offset_image = np.fft.ifftn(offset_image) 31print(f'Known offset (y, x): {shift}') 32 33# pixel precision first 34shift, error, diffphase = phase_cross_correlation(image, offset_image) 35 36fig = plt.figure(figsize=(8, 3)) 37ax1 = plt.subplot(1, 3, 1) 38ax2 = plt.subplot(1, 3, 2, sharex=ax1, sharey=ax1) 39ax3 = plt.subplot(1, 3, 3) 40 41ax1.imshow(image, cmap='gray') 42ax1.set_axis_off() 43ax1.set_title('Reference image') 44 45ax2.imshow(offset_image.real, cmap='gray') 46ax2.set_axis_off() 47ax2.set_title('Offset image') 48 49# Show the output of a cross-correlation to show what the algorithm is 50# doing behind the scenes 51image_product = np.fft.fft2(image) * np.fft.fft2(offset_image).conj() 52cc_image = np.fft.fftshift(np.fft.ifft2(image_product)) 53ax3.imshow(cc_image.real) 54ax3.set_axis_off() 55ax3.set_title("Cross-correlation") 56 57plt.show() 58 59print(f'Detected pixel offset (y, x): {shift}') 60 61# subpixel precision 62shift, error, diffphase = phase_cross_correlation(image, offset_image, 63 upsample_factor=100) 64 65fig = plt.figure(figsize=(8, 3)) 66ax1 = plt.subplot(1, 3, 1) 67ax2 = plt.subplot(1, 3, 2, sharex=ax1, sharey=ax1) 68ax3 = plt.subplot(1, 3, 3) 69 70ax1.imshow(image, cmap='gray') 71ax1.set_axis_off() 72ax1.set_title('Reference image') 73 74ax2.imshow(offset_image.real, cmap='gray') 75ax2.set_axis_off() 76ax2.set_title('Offset image') 77 78# Calculate the upsampled DFT, again to show what the algorithm is doing 79# behind the scenes. Constants correspond to calculated values in routine. 80# See source code for details. 81cc_image = _upsampled_dft(image_product, 150, 100, (shift*100)+75).conj() 82ax3.imshow(cc_image.real) 83ax3.set_axis_off() 84ax3.set_title("Supersampled XC sub-area") 85 86 87plt.show() 88 89print(f'Detected subpixel offset (y, x): {shift}') 90