模板匹配&&频域相位相关

        最近在做图像拼接项目，原打算用大众化的sift点匹配然后拼接。发现sift算法运算量巨大，不能满足拼接的实时要求，同时特征点法需要不同程度的干预（比如，需要随机抽样来删选优质的匹配点对），不满足自动匹准的要求。最后，只得尝试着用基于区域匹配的方法。    

      基于区域的方法其实就是模板匹配了，opencv已经为我们提供函数:

void cv::matchTemplate( InputArray _img, InputArray _templ, OutputArray _result, int method )；
        模板匹配就是在一幅图像中寻找和模板图像最相似的区域。该函数的功能为，在输入源图像中滑动框，寻找各个位置与模板图像的相似度，并将结果保存在结果矩阵中。该矩阵的每一个点的灰度表示与模板T的匹配程度。然后可以通过函数minMaxLoc（）定位矩阵中的最大值（or最小值）。

匹配的方法有：

CV_TM_SQDIFF 平方差匹配法，最好的匹配为0，值越大匹配越差

CV_TM_SQDIFF_NORMED 归一化平方差匹配法

CV_TM_CCORR 相关匹配法，采用乘法操作，数值越大表明匹配越好

CV_TM_CCORR_NORMED 归一化相关匹配法

CV_TM_CCOEFF 相关系数匹配法，最好的匹配为1，-1表示最差的匹配

CV_TM_CCOEFF_NORMED 归一化相关系数匹配法

前面两种方法为越小的值表示越匹配，后四种方法值越大越匹配。

      当然，纯粹调用opencv的模板匹配来拼接图像，那真是太naive了。参考模板匹配，提出利用频域相位相关法拼接，大致思路如下：

1，对图像做必要的预处理；

2，对图像进行傅里叶变换，并进行相关度计算，求取两幅图像的重叠区域；

3，在重叠区域选取合适的子块计算图像的旋转因子。

4，线性融合。

最后稍微了解下opencv提供的源码:

void matchTemplate( const Mat& _img, const Mat& _templ, Mat& result, int method )
{
CV_Assert( CV_TM_SQDIFF <= method && method <= CV_TM_CCOEFF_NORMED );
//numType用来表示模板匹配的方式，0表示相关匹配法，1表示相关系数匹配法，2表示平方差匹配法
//isNormed表示是否进行归一化处理，true表示进行归一化，false表示不进行归一化处理
int numType = method == CV_TM_CCORR || method == CV_TM_CCORR_NORMED ? 0 :
method == CV_TM_CCOEFF || method == CV_TM_CCOEFF_NORMED ? 1 : 2;
bool isNormed = method == CV_TM_CCORR_NORMED ||
method == CV_TM_SQDIFF_NORMED ||
method == CV_TM_CCOEFF_NORMED;
//判断两幅图像的大小关系，如果输入的原始图像比匹配图像要小，则将原始图像作为模板，原来的模板图像作为搜索图
Mat img = _img, templ = _templ;
if( img.rows < templ.rows || img.cols < templ.cols )
std::swap(img, templ);

CV_Assert( (img.depth() == CV_8U || img.depth() == CV_32F) &&
img.type() == templ.type() );

//crossCorr函数是将输入图像做了一次DFT变换（离散傅里叶变换），将空间域的图像转换到频率域中来进行处理，并将处理的结果存放在result中
int cn = img.channels();
crossCorr( img, templ, result,
Size(img.cols - templ.cols + 1, img.rows - templ.rows + 1),
CV_32F, Point(0,0), 0, 0);

//如果是相关匹配方法，此处已经计算完毕，返回
if( method == CV_TM_CCORR )
return;

//将模板看作单位1，计算每一个像元所占的百分比（也可以理解为整个模板面积为1，计算每个像元的面积）
double invArea = 1./((double)templ.rows * templ.cols);

Mat sum, sqsum;
Scalar templMean, templSdv;
double *q0 = 0, *q1 = 0, *q2 = 0, *q3 = 0;
double templNorm = 0, templSum2 = 0;

//相关系数匹配算法
if( method == CV_TM_CCOEFF )
{
integral(img, sum, CV_64F);//对原始图像进行求和
templMean = mean(templ);//计算模板图像的均值向量
}
else//其他匹配算法
{
integral(img, sum, sqsum, CV_64F);//计算原始图像的和以及平方和
meanStdDev( templ, templMean, templSdv );//计算模板图像的均值向量和方差向量

templNorm = CV_SQR(templSdv[0]) + CV_SQR(templSdv[1]) +
CV_SQR(templSdv[2]) + CV_SQR(templSdv[3]);//计算所有通道的方差和

if( templNorm < DBL_EPSILON && method == CV_TM_CCOEFF_NORMED )
{//如果所有通道的方差的和等于0，并且使用的方法是归一化相关系数匹配方法，则返回
result = Scalar::all(1);
return;
}

templSum2 = templNorm +
CV_SQR(templMean[0]) + CV_SQR(templMean[1]) +
CV_SQR(templMean[2]) + CV_SQR(templMean[3]);//计算所有通道的均值的平方和

if( numType != 1 )//匹配方式不是相关系数，对模板均值向量和templNorm重新赋值
{
templMean = Scalar::all(0);
templNorm = templSum2;
}

templSum2 /= invArea;
templNorm = sqrt(templNorm);
templNorm /= sqrt(invArea); // care of accuracy here

q0 = (double*)sqsum.data;
q1 = q0 + templ.cols*cn;
q2 = (double*)(sqsum.data + templ.rows*sqsum.step);
q3 = q2 + templ.cols*cn;
}

int sumstep = sum.data ? (int)(sum.step / sizeof(double)) : 0;
int sqstep = sqsum.data ? (int)(sqsum.step / sizeof(double)) : 0;

int i, j, k;

for( i = 0; i < result.rows; i++ )
{
float* rrow = (float*)(result.data + i*result.step);
int idx = i * sumstep;
int idx2 = i * sqstep;

for( j = 0; j < result.cols; j++, idx += cn, idx2 += cn )
{
double num = rrow[j], t;
double wndMean2 = 0, wndSum2 = 0;

if( numType == 1 )
{
for( k = 0; k < cn; k++ )
{
t = p0[idx+k] - p1[idx+k] - p2[idx+k] + p3[idx+k];
wndMean2 += CV_SQR(t);
num -= t*templMean[k];
}

wndMean2 *= invArea;
}

if( isNormed || numType == 2 )
{
for( k = 0; k < cn; k++ )
{
t = q0[idx2+k] - q1[idx2+k] - q2[idx2+k] + q3[idx2+k];
wndSum2 += t;
}

if( numType == 2 )
{
num = wndSum2 - 2*num + templSum2;
num = MAX(num, 0.);
}
}

if( isNormed )
{
t = sqrt(MAX(wndSum2 - wndMean2,0))*templNorm;
if( fabs(num) < t )
num /= t;
else if( fabs(num) < t*1.125 )
num = num > 0 ? 1 : -1;
else
num = method != CV_TM_SQDIFF_NORMED ? 0 : 1;
}

rrow[j] = (float)num;
}
}
}