频域相位相关模板匹配

• • 相位相关基本概念
  • 频域相位相关模板匹配
  • 相关链接

相位相关基本概念

相位相关 (Phase Correlation) 是一种可用于估计两个相似图像之间的相对平移偏移的方法。它通常用于图像匹配，并依赖于图像的频域表示，通常通过快速傅里叶变换计算。以下公式说明了使用相位相关在目标图像内匹配模板图像的步骤，本节中，我们使用频域相位解决相关模板匹配：

FB=F{fB},FT=F{fT}R=FB∘FT∗∣FB∘FT∗∣r=F−1{R}(Δx,Δy)=argmaxx,y{r}F_B=\mathcal F\{f_B\},F_T=\mathcal F\{f_T\} \\ R=\frac {F_B\circ F_T^*} {|F_B\circ F_T^*|} \\ r=\mathcal F^{-1}\{R\} \\ (\Delta x,\Delta y)=arg\underset{x,y}{max}\{r\} FB​=F{fB​},FT​=F{fT​}R=∣FB​∘FT∗​∣FB​∘FT∗​​r=F−1{R}(Δx,Δy)=argx,ymax​{r}

其中，fBf_BfB​ 是目标图像，fTf_TfT​ 是模板图像，F\mathcal FF 是傅里叶变换函数，∘\circ∘ 是逐元素操作，RRR 是交叉功率谱，rrr 是归一化相位相关。

频域相位相关模板匹配

接下来，我们使用 Python 执行频域相位相关模板匹配。

(1) 首先，导入所需库，并读取目标和模板图像：

import scipy.fftpack as fp
from skimage.io import imread
from skimage.color import rgb2gray, gray2rgb
from skimage.draw import rectangle_perimeter
import numpy as np
import matplotlib.pylab as plt
from mpl_toolkits.mplot3d import Axes3D  # noqa: F401 unused import
from mpl_toolkits.axes_grid1 import make_axes_locatable
from matplotlib.ticker import LinearLocator, FormatStrFormatterdef plot_3d(X, Y, Z, cmap=plt.cm.seismic):fig = plt.figure(figsize=(20,20))ax = fig.gca(projection='3d')# Plot the surface.surf = ax.plot_surface(X, Y, Z, cmap=cmap, linewidth=5, antialiased=False)#ax.plot_wireframe(X, Y, Z, rstride=10, cstride=10)#ax.set_zscale("log", nonposx='clip')#ax.zaxis.set_scale('log')ax.zaxis.set_major_locator(LinearLocator(10))ax.zaxis.set_major_formatter(FormatStrFormatter('%.02f'))ax.set_xlabel('F1', size=20)ax.set_ylabel('F2', size=20)ax.set_zlabel('Freq Response', size=20)#ax.set_zlim((-40,10))# Add a color bar which maps values to colors.fig.colorbar(surf) #, shrink=0.15, aspect=10)#plt.title('Frequency Response of the Gaussian Kernel')plt.show()    im = 255*rgb2gray(imread('1.png'))
im_tm = 255*rgb2gray(imread('1.png'))[100:250,200:350]

(2) 计算两个图像的离散 2D 傅立叶变换：

F = fp.fftn(im)
F_tm = fp.fftn(im_tm, shape=im.shape)

(3) 通过获取第二个结果的复共轭，按元素将傅里叶变换相乘，并按元素归一化乘积结果，计算交叉功率谱；通过应用逆傅立叶变换来获得归一化的互相关值；计算最佳的匹配位置，并使用 np.argmax() 函数计算确定归一化互相关矩阵中峰值的位置：

F_cc = F * np.conj(F_tm)
c = (fp.ifftn(F_cc/np.abs(F_cc))).real
i, j = np.unravel_index(c.argmax(), c.shape)
print(i, j)

(4) 围绕最佳匹配位置绘制矩形：

im2 = (gray2rgb(im)).astype(np.uint8)
rr, cc = rectangle_perimeter((i,j), end=(i + im_tm.shape[0], j + im_tm.shape[1]), shape=im.shape)
for x in range(-2,2):for y in range(-2,2):im2[rr + x, cc + y] = (255,0,0)

(5) 绘制模板和目标图像以及最佳匹配的位置。另外，绘制交叉功率谱以可视化峰值：

plt.figure(figsize=(2,3))
plt.gray()
plt.imshow(im_tm), plt.title('template', size=10), plt.axis('off')
plt.show()
fig, ax = plt.subplots(1, 2, sharey=True, figsize=(12,7))
ax[0].imshow(im), ax[0].set_title('target', size=10)
ax[1].imshow(im2), ax[1].set_title('matched template', size=10)
for a in ax.ravel():a.set_axis_off()
plt.tight_layout()
plt.show()
Y = np.arange(F_cc.shape[0])
X = np.arange(F_cc.shape[1])
X, Y = np.meshgrid(X, Y)
Z = c
plot_3d(X,Y,Z, cmap='YlOrRd')

相关链接

Python图像处理【1】图像与视频处理基础
Python图像处理【2】探索Python图像处理库
Python图像处理【3】Python图像处理库应用
Python图像处理【4】图像线性变换
Python图像处理【5】图像扭曲/逆扭曲
Python图像处理【6】通过哈希查找重复和类似的图像
Python图像处理【7】采样、卷积与离散傅里叶变换
Python图像处理【8】使用低通滤波器模糊图像
Python图像处理【9】使用高通滤波器执行边缘检测