使用极坐标变换和对数极坐标变换进行配准

相位相关 (registration.phase_cross_correlation )是确定相似图像对之间的平移偏移的有效方法。然而，这种方法依赖于图像之间几乎没有旋转/缩放差异，这在真实世界的示例中是典型的。

利用对数极坐标变换和频域平移不变性这两个几何性质，可以恢复两幅图像之间的旋转和缩放差异。首先，笛卡尔空间中的旋转变成了沿角坐标的平移 (\(\theta\) )对数极空间的轴。其次，笛卡尔空间中的缩放转换为沿径向坐标的平移 (\(\rho = \ln\sqrt{{x^2 + y^2}}\) )的对数极空间。最后，空间域的平移差异不会影响频域的震级谱。

在这一系列示例中，我们以这些概念为基础来展示对数极坐标的转换 transform.warp_polar 可以与相位相关一起使用，以恢复也具有平移偏移的两个图像之间的旋转和缩放差异。

使用极坐标变换恢复旋转差

在第一个示例中，我们考虑两个图像的简单情况，这两个图像仅在围绕共同中心点的旋转方面不同。通过将这些图像重新映射到极空间，旋转差变成简单的平移差，可以通过相位相关来恢复。

import numpy as np
import matplotlib.pyplot as plt

from skimage import data
from skimage.registration import phase_cross_correlation
from skimage.transform import warp_polar, rotate, rescale
from skimage.util import img_as_float

radius = 705
angle = 35
image = data.retina()
image = img_as_float(image)
rotated = rotate(image, angle)
image_polar = warp_polar(image, radius=radius, multichannel=True)
rotated_polar = warp_polar(rotated, radius=radius, multichannel=True)

fig, axes = plt.subplots(2, 2, figsize=(8, 8))
ax = axes.ravel()
ax[0].set_title("Original")
ax[0].imshow(image)
ax[1].set_title("Rotated")
ax[1].imshow(rotated)
ax[2].set_title("Polar-Transformed Original")
ax[2].imshow(image_polar)
ax[3].set_title("Polar-Transformed Rotated")
ax[3].imshow(rotated_polar)
plt.show()

shifts, error, phasediff = phase_cross_correlation(image_polar, rotated_polar)
print("Expected value for counterclockwise rotation in degrees: "
      f"{angle}")
print("Recovered value for counterclockwise rotation: "
      f"{shifts[0]}")

输出：

/scikit-image/doc/examples/registration/plot_register_rotation.py:48: FutureWarning:

`multichannel` is a deprecated argument name for `warp_polar`. It will be removed in version 1.0. Please use `channel_axis` instead.

/scikit-image/doc/examples/registration/plot_register_rotation.py:49: FutureWarning:

`multichannel` is a deprecated argument name for `warp_polar`. It will be removed in version 1.0. Please use `channel_axis` instead.

Expected value for counterclockwise rotation in degrees: 35
Recovered value for counterclockwise rotation: 35.0

使用对数极坐标变换恢复旋转和缩放差异

在第二个示例中，图像在旋转和缩放方面都不同(请注意轴刻度值)。通过将这些图像重新映射到对数极空间，我们可以像以前一样恢复旋转，现在还可以通过相位相关来进行缩放。

# radius must be large enough to capture useful info in larger image
radius = 1500
angle = 53.7
scale = 2.2
image = data.retina()
image = img_as_float(image)
rotated = rotate(image, angle)
rescaled = rescale(rotated, scale, multichannel=True)
image_polar = warp_polar(image, radius=radius,
                         scaling='log', multichannel=True)
rescaled_polar = warp_polar(rescaled, radius=radius,
                            scaling='log', multichannel=True)

fig, axes = plt.subplots(2, 2, figsize=(8, 8))
ax = axes.ravel()
ax[0].set_title("Original")
ax[0].imshow(image)
ax[1].set_title("Rotated and Rescaled")
ax[1].imshow(rescaled)
ax[2].set_title("Log-Polar-Transformed Original")
ax[2].imshow(image_polar)
ax[3].set_title("Log-Polar-Transformed Rotated and Rescaled")
ax[3].imshow(rescaled_polar)
plt.show()

# setting `upsample_factor` can increase precision
shifts, error, phasediff = phase_cross_correlation(image_polar, rescaled_polar,
                                                   upsample_factor=20)
shiftr, shiftc = shifts[:2]

# Calculate scale factor from translation
klog = radius / np.log(radius)
shift_scale = 1 / (np.exp(shiftc / klog))

print(f"Expected value for cc rotation in degrees: {angle}")
print(f"Recovered value for cc rotation: {shiftr}")
print()
print(f"Expected value for scaling difference: {scale}")
print(f"Recovered value for scaling difference: {shift_scale}")

输出：

/scikit-image/doc/examples/registration/plot_register_rotation.py:85: FutureWarning:

`multichannel` is a deprecated argument name for `rescale`. It will be removed in version 1.0. Please use `channel_axis` instead.

/scikit-image/doc/examples/registration/plot_register_rotation.py:86: FutureWarning:

`multichannel` is a deprecated argument name for `warp_polar`. It will be removed in version 1.0. Please use `channel_axis` instead.

/scikit-image/doc/examples/registration/plot_register_rotation.py:88: FutureWarning:

`multichannel` is a deprecated argument name for `warp_polar`. It will be removed in version 1.0. Please use `channel_axis` instead.

Expected value for cc rotation in degrees: 53.7
Recovered value for cc rotation: 53.85

Expected value for scaling difference: 2.2
Recovered value for scaling difference: 2.199797163552113

转译图像上的寄存器旋转和缩放.第1部分

仅当要配准的图像共享一个中心时，上述示例才有效。然而，更常见的情况是，还存在对要配准的两个图像之间的差异的转换组件。记录旋转、缩放和平移的一种方法是首先校正旋转和缩放，然后求解平移。通过对傅立叶变换图像的幅度谱进行处理，可以解决平移图像的旋转和比例差异。

在下一个示例中，我们首先展示当两个图像在旋转、缩放和平移方面不同时，上述方法是如何失败的。

from skimage.color import rgb2gray
from skimage.filters import window, difference_of_gaussians
from scipy.fftpack import fft2, fftshift

angle = 24
scale = 1.4
shiftr = 30
shiftc = 15

image = rgb2gray(data.retina())
translated = image[shiftr:, shiftc:]
rotated = rotate(translated, angle)
rescaled = rescale(rotated, scale)
sizer, sizec = image.shape
rts_image = rescaled[:sizer, :sizec]

# When center is not shared, log-polar transform is not helpful!
radius = 705
warped_image = warp_polar(image, radius=radius, scaling="log")
warped_rts = warp_polar(rts_image, radius=radius, scaling="log")
shifts, error, phasediff = phase_cross_correlation(warped_image, warped_rts,
                                                   upsample_factor=20)
shiftr, shiftc = shifts[:2]
klog = radius / np.log(radius)
shift_scale = 1 / (np.exp(shiftc / klog))

fig, axes = plt.subplots(2, 2, figsize=(8, 8))
ax = axes.ravel()
ax[0].set_title("Original Image")
ax[0].imshow(image, cmap='gray')
ax[1].set_title("Modified Image")
ax[1].imshow(rts_image, cmap='gray')
ax[2].set_title("Log-Polar-Transformed Original")
ax[2].imshow(warped_image)
ax[3].set_title("Log-Polar-Transformed Modified")
ax[3].imshow(warped_rts)
fig.suptitle('log-polar-based registration fails when no shared center')
plt.show()

print(f"Expected value for cc rotation in degrees: {angle}")
print(f"Recovered value for cc rotation: {shiftr}")
print()
print(f"Expected value for scaling difference: {scale}")
print(f"Recovered value for scaling difference: {shift_scale}")

输出：

Expected value for cc rotation in degrees: 24
Recovered value for cc rotation: 9.05

Expected value for scaling difference: 1.4
Recovered value for scaling difference: 1.0476106934350098

平移图像上的寄存器旋转和缩放.第2部分

接下来，我们将展示旋转和缩放差异，而不是平移差异，在图像的频率幅度谱中是如何明显的。这些差异可以通过将震级谱视为图像本身，并应用上面采用的相同的对数-极+相位相关方法来恢复。

# First, band-pass filter both images
image = difference_of_gaussians(image, 5, 20)
rts_image = difference_of_gaussians(rts_image, 5, 20)

# window images
wimage = image * window('hann', image.shape)
rts_wimage = rts_image * window('hann', image.shape)

# work with shifted FFT magnitudes
image_fs = np.abs(fftshift(fft2(wimage)))
rts_fs = np.abs(fftshift(fft2(rts_wimage)))

# Create log-polar transformed FFT mag images and register
shape = image_fs.shape
radius = shape[0] // 8  # only take lower frequencies
warped_image_fs = warp_polar(image_fs, radius=radius, output_shape=shape,
                             scaling='log', order=0)
warped_rts_fs = warp_polar(rts_fs, radius=radius, output_shape=shape,
                           scaling='log', order=0)

warped_image_fs = warped_image_fs[:shape[0] // 2, :]  # only use half of FFT
warped_rts_fs = warped_rts_fs[:shape[0] // 2, :]
shifts, error, phasediff = phase_cross_correlation(warped_image_fs,
                                                   warped_rts_fs,
                                                   upsample_factor=10)

# Use translation parameters to calculate rotation and scaling parameters
shiftr, shiftc = shifts[:2]
recovered_angle = (360 / shape[0]) * shiftr
klog = shape[1] / np.log(radius)
shift_scale = np.exp(shiftc / klog)

fig, axes = plt.subplots(2, 2, figsize=(8, 8))
ax = axes.ravel()
ax[0].set_title("Original Image FFT\n(magnitude; zoomed)")
center = np.array(shape) // 2
ax[0].imshow(image_fs[center[0] - radius:center[0] + radius,
                      center[1] - radius:center[1] + radius],
             cmap='magma')
ax[1].set_title("Modified Image FFT\n(magnitude; zoomed)")
ax[1].imshow(rts_fs[center[0] - radius:center[0] + radius,
                    center[1] - radius:center[1] + radius],
             cmap='magma')
ax[2].set_title("Log-Polar-Transformed\nOriginal FFT")
ax[2].imshow(warped_image_fs, cmap='magma')
ax[3].set_title("Log-Polar-Transformed\nModified FFT")
ax[3].imshow(warped_rts_fs, cmap='magma')
fig.suptitle('Working in frequency domain can recover rotation and scaling')
plt.show()

print(f"Expected value for cc rotation in degrees: {angle}")
print(f"Recovered value for cc rotation: {recovered_angle}")
print()
print(f"Expected value for scaling difference: {scale}")
print(f"Recovered value for scaling difference: {shift_scale}")

输出：

Expected value for cc rotation in degrees: 24
Recovered value for cc rotation: 0.05102763997165131

Expected value for scaling difference: 1.4
Recovered value for scaling difference: 0.46663985889861054

关于这种方法的几点注解

应该指出的是，这种方法依赖于必须提前选择的几个参数，对于这些参数，没有明显的最佳选择：

1.图像应该应用一定程度的带通滤波，特别是为了去除高频，这里的不同选择可能会影响结果。带通滤光器还使问题复杂化，因为由于要配准的图像将在比例上不同，并且这些比例差异是未知的，所以任何带通滤光器都将必然衰减图像之间的不同特征。例如，在这里的最后一个例子中，对数极坐标变换的星等谱看起来并不“相似”。然而，如果你仔细观察，这些光谱中有一些常见的模式，它们确实通过相位相关得到了很好的对准，正如所展示的那样。

2.必须使用圆形对称的窗口对图像进行加窗，以消除图像边界带来的光谱泄漏。没有明确的最佳窗口选择。

最后，我们注意到，尺度的巨大变化将极大地改变震级谱，特别是因为尺度的巨大变化通常会伴随着一些裁剪和信息内容的损失。文献建议保持在1.8-2倍的比例变化范围内 1 2 。这对大多数生物成像应用来说都是很好的。

参考文献

1

B.S.Reddy和B.N.Chatterji。一种基于FFT的平移、旋转和比例不变图像配准技术。IEEE传输图像处理，5(8)：1266-1271,1996。 DOI:10.1109/83.506761

2

Tzimiropoulos，Georgios和Tania Stathaki。“基于FFT的稳健尺度不变图像配准。”第四届海洋直接转矩控制技术会议。2009年。 DOI:10.1109/TPAMI.2010.107