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