注意
前往末尾以下載完整的範例程式碼。或透過 Binder 在您的瀏覽器中執行此範例
使用極座標和對數極座標轉換進行配準#
相位相關 (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, channel_axis=-1)
rotated_polar = warp_polar(rotated, radius=radius, channel_axis=-1)
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, normalization=None
)
print(f'Expected value for counterclockwise rotation in degrees: ' f'{angle}')
print(f'Recovered value for counterclockwise rotation: ' f'{shifts[0]}')

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, channel_axis=-1)
image_polar = warp_polar(image, radius=radius, scaling='log', channel_axis=-1)
rescaled_polar = warp_polar(rescaled, radius=radius, scaling='log', channel_axis=-1)
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, normalization=None
)
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}')

Expected value for cc rotation in degrees: 53.7
Recovered value for cc rotation: 53.75
Expected value for scaling difference: 2.2
Recovered value for scaling difference: 2.1981889915232165
在平移影像上註冊旋轉和縮放 - 第 1 部分#
只有當要註冊的影像共享一個中心時,以上範例才有效。然而,更常見的情況是,兩個要註冊的影像之間的差異也包含平移分量。註冊旋轉、縮放和平移的一種方法是先校正旋轉和縮放,然後求解平移。可以透過處理傅立葉變換影像的幅度頻譜來解決平移影像的旋轉和縮放差異。
在下一個範例中,我們首先展示當兩個影像在旋轉、縮放和平移方面存在差異時,上述方法如何失敗。
from skimage.color import rgb2gray
from skimage.filters import window, difference_of_gaussians
from scipy.fft 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, normalization=None
)
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: -167.55
Expected value for scaling difference: 1.4
Recovered value for scaling difference: 25.110458986143573
在平移影像上註冊旋轉和縮放 - 第 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, normalization=None
)
# 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: 23.753366406803682
Expected value for scaling difference: 1.4
Recovered value for scaling difference: 1.3901762721757436
關於這種方法的一些注意事項#
應該注意的是,這種方法依賴於一些必須預先選擇的參數,並且沒有明確的最佳選擇
1. 影像應該應用一定程度的帶通濾波,特別是去除高頻,並且此處的不同選擇可能會影響結果。帶通濾波器也會使情況複雜化,因為由於要註冊的影像在比例上會有所不同,並且這些比例差異未知,因此任何帶通濾波器都必然會衰減影像之間的不同特徵。例如,對數極座標轉換後的幅度頻譜在最後一個範例中看起來並不「相似」。然而,如果您仔細觀察,這些頻譜中確實存在一些共同的模式,並且它們最終會如示範那樣通過相位相關很好地對齊。
2. 必須使用具有圓對稱性的視窗對影像進行視窗化,以消除來自影像邊界的頻譜洩漏。沒有明確的最佳視窗選擇。
最後,我們注意到,比例的巨大變化會顯著改變幅度頻譜,特別是因為比例的巨大變化通常會伴隨著一些裁剪和資訊內容的損失。文獻建議保持在 1.8-2 倍的比例變化範圍內 [1] [2]。這對於大多數生物影像應用來說是足夠的。
參考文獻#
腳本的總執行時間:(0 分鐘 7.940 秒)