Matlab模板匹配实现图像运动估计（频域实现）

引言：在利用显微镜观测量物体的长度时，其精度非常高，然而，其行程却有限。一个测量较长物体的长度时，一种方法是在物体的左侧开始拍摄图像，然后移动相机（或者物体），再拍图像，直到拍摄到物体的右侧为止。对于其中拍摄的相邻的图像，如下图所示（图片由余师兄提供），需要计算它们的运动距离。

模板匹配：对于上面的图像，可以明显看出其往右移了，但对于精密测量来讲，求其移动的距离并精确到像素级非常重要。方法是在右图像上移动左图像，并比较其重叠区域的相似度程度，当其足够相似时，则认为找到了左右图像的关系。
一般的，相似度程度可以由SSD（误差平方和）或NCC（归一化积相关）来表征，当然，还有很多其它的方式可以表征，如SAD，但是个人认为SSD和NCC的性能更好，且其为线性函数，可以进行频域实现。
对于SSD，其评价函数为：

H(p,q)=∑U∑U(f(x,y)−g(i,j))2N(p,q)

其中的f和g为需要匹配的图像，U为f和g在位置(p,q)时的重叠区域，f和g的下标(x,y)和(i,j)不同，是因为重叠区域U在f和g的位置不同引起的，N(p,q)表示重叠区域U的点的个数。上式可以化简为：

H(p,q)=∑U∑Uf(x,y)2+∑U∑Ug(i,j)2−2∑U∑Uf(x,y)g(i,j)N(p,q)

对于∑U∑Uf(x,y)2和 ∑U∑Ug(x,y)2，可以用如下Matlab程序实现：

function result = ImageAllCumSum(curr)[m, n] = size(curr); result = zeros(m*2-1, n*2-1); result(1:m, 1:n) = cumsum(cumsum(curr, 1), 2); result(1:m, n:2*n-1) = flip(cumsum(cumsum(flip(curr, 2), 1), 2), 2); result(m:2*m-1, 1:n) =  flip(cumsum(cumsum(flip(curr, 1), 1), 2), 1); result(m:2*m-1, n:2*n-1) = rot90(cumsum(cumsum(rot90(curr, 2), 1), 2), 2);

计算f和g的平方和只需输入

result = ImageAllCumSum(f.*f);

而对于N(p,q)可以以如下方式实现：

function result = ImageAllSize(curr) [m, n] = size(curr);resultTemp = ones(m, n);result = ImageAllCumSum(resultTemp);

对于∑U∑Uf(x,y)g(i,j)，利用卷积进行计算将十分便利。程序如下：

function result = ImageAllConvolution(prev, curr) prevFFT = ImageAllFFT(rot90(prev, 2));currFFT = ImageAllFFT(curr);result = ifft2(currFFT.*prevFFT);function result = ImageAllFFT(curr)[m, n] = size(curr);result = zeros(m*2-1, n*2-1); result(1:m,1:n) = curr; result = fft2(result);

对于计算上述f和g的模板匹配为：

function result = ImageAllFftSsd(prev, curr) convRes = ImageAllConvolution(prev, curr);squareCumSumPrev = ImageAllCumSum(rot90(prev.*prev, 2));squareCumSumCurr = ImageAllCumSum(curr.*curr);sizeImage = ImageAllSize(prev);result = (squareCumSumPrev + squareCumSumCurr - 2*real(convRes))./sizeImage;

对于NCC，其评价函数为：

H(p,q)=∑U∑U(f(x,y)−f−(x,y))(g(x,y)−g−(x,y))∑U∑U(f(x,y)−f−(x,y))2∑U∑U(g(x,y)−g−(x,y))2−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−√

其中的f和g为需要匹配的图像，U为f和g在位置(p,q)时的重叠区域，f和g的下标(x,y)和(i,j)不同，是因为重叠区域U在f和g的位置不同引起的，f−(x,y)和g−(x,y)表示重叠区域U内f和g的均值。上式可以化简为：

H(p,q)=∑U∑Uf(x,y)g(x,y)−∑U∑Uf−(x,y)g−(i,j)∑U∑U(f(x,y)2−f−(x,y)2)∑U∑U(g(x,y)2−g−(x,y)2)−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−√

NCC可以用如下Matlab程序实现：

function result = ImageAllFftNcc(prev, curr)convRes = ImageAllConvolution(prev, curr); sizeRes = ImageAllSize(curr); CumSumPrev = ImageAllCumSum(rot90(prev, 2));CumSumCurr = ImageAllCumSum(curr); squareCumSumPrev = ImageAllCumSum(rot90(prev.*prev, 2)); squareCumSumCurr = ImageAllCumSum(curr.*curr); result = (convRes - CumSumPrev.*CumSumCurr./sizeRes);resultDonominator = sqrt((squareCumSumPrev - CumSumPrev.*CumSumPrev./sizeRes) .* ...(squareCumSumCurr - CumSumCurr.*CumSumCurr./sizeRes));result = result./resultDonominator;

以下为程序的Demo：

im1 = rgb2gray(im2double(imread('1.jpg')));im2 = rgb2gray(im2double(imread('2.jpg')));SSD = ImageAllFftSsd(im1, im2);SSD = SSD(11:end-10, 11:end-10);[x, y] = meshgrid(1:size(SSD, 2), 1:size(SSD, 1));surf(x, y, SSD);NCC = ImageAllFftNnc(im1, im2);NCC = NCC(11:end-10, 11:end-10);[x, y] = meshgrid(1:size(NCC, 2), 1:size(NCC, 1));surf(x, y, NCC);