图像比较之模板匹配

1.模板匹配基本原理概述

      当我们比较两幅图像的时候，首先面对的基本问题是：什么时候两幅图像才是一样或比较相似的，这两幅图像的相似程度如何衡量？当然，比较一般的方法是，当两幅图像的所有像素灰度值一样的时候，我们认为这样幅图是一样的。这种比较方法在某些特定的应用领域是可行的，比如在恒定光照环境和相机内部环境下，检测连续两帧图像的变化。简单的比较像素之间的差值在大多数应用场合下是不太合适的。噪声、量化误差、微小的光照变化、极小的平移或旋转在两幅图像像素简单差值的时候将产生较大的像素差值，但是人眼观察看来，两幅图像仍然是一样的。显然，人类知觉能够感应更为广泛的相似性内容，即使在单个像素差值较大的情况下，也能够利用诸如结构和图像内容来识别两幅图像的相似性。在结构或语义层面比较图像是一个比较困难的问题，同时也是一个比较有趣的研究领域。

      在这我们介绍的是一种在像素层面相对简单的图像比较方法。即在一幅较大图像定位一幅给定的子图像-----模板图像，即通常所说的模板匹配。这种情况经常发生，比如在一个图像场景中定位一个特定的物体，或者是在图像序列中追踪最某些特定模式。模板匹配的基本原理很简单：在待搜寻的图像中，移动模板图像，在每一个位置测量待搜寻图像的子图像和模板图像的差值，当相似度达到最大时，记录其相应的位置。但实际情况却不是这么简单，如何选取合适的距离测量方法？当亮度或对比度改变时应该如何处理？匹配中总的距离差值为多少才可以被认为是相似度比较高？这些问题都需要根据实用情况加以考虑。

     以上图像都“相同”吗？在原始图像（a）和其他五幅图像(b-f)中,简单的比较图像像素间的差值将会得到较大的距离值。

       模板匹配的基本原理。以两幅图像的原点为参考点，参考图像R在待搜寻的图像I中平移（r,s）个单位，待搜寻图像的大小和模板图像的大小确定的最大的搜寻区域。

2.灰度图像中的模板匹配

     灰度图像中的模板匹配，主要是找到模板图像R和搜寻图像I中的子图像相同或最相似的位置。以下式子表示模板在原始图像偏移（r,s）单位，

    参考图像R分别在水平方法和竖直方向平移r和s个单位，那么模板匹配的问题可以概述为：

   给定搜寻图像I和模板图像R。找到平移后的参考图像以及搜寻图像中所对应的子图像相似度最大的位置。

2.1图像匹配中的距离函数

   模板匹配中最重要的是找到一种对灰度和对比度变化具有较好鲁棒性的相似度测量函数。

    为了测量图像间的相似度程度，我们计算每一次平移（r,s）后的参考图像和搜寻图像中所对应的子图像的“距离”d(r,s)（如下图所示）。二维灰度图像中的测量函数有如下基本的几种：

二维图像中距离测量函数的表示

Sum ofabsolute differences:

Maximumdifference:

                                                                

Sum ofsquared differences:

                                                                                    

由于SSD函数的特性，在统计和最优化领域被经常用到。为了找到参考图像在搜寻图像中的最佳匹配位置，只要使得SSD函数取得最小值即可。即    

         上式中的B是参考图像中所有像素灰度值的平方和，是一个常量（独立于r,s）因此在计算其最小值的时候可以忽略。A(r,s)表示的是搜寻图像在（r,s）的子图像中所有像素灰度值的平方和，C(r,s)即称为搜寻图像和参考图像的线性互相关函数，通常可以表示为：

                                                             

     当R 和I超出边界时，其值为零，故上式亦可表示为：

      如果我们假设A(r,s)在搜寻图像中是一个常量，那么SSD中计算其最佳匹配位置的时候，其值是可以忽略的，且当C(r,s)取得最大值的时候，参考图像和搜寻图像中的子图像最相似。这种情况下，SSD的最小值可以通过计算C(r,s)的最大值获得。

2.2归一化的互相关

    实际应用当中，上述假定A(r,s)为一常量并不总是成立。因此，上述的互相关的结果将随着搜寻图像中像素灰度值的改变而产生较大变化。归一化的互相关把参考图像和当前的子图像的能量考虑进去：

     当参考图像和搜寻图像子图像灰度值都是正数的时候，Cn(r,s)的值总是为[0,1]范围内，与图像其他像素灰度值无关。当Cn(r,s)等于1时，表明在平移位置（r,s），参考图像和子图像达到最大相似度；相反，当Cn(r,s)等于0时，表明在平移位置（r,s），参考图像和子图像完全不匹配。当子图像中的所有灰度值改变时，归一化的互相关Cn(r,s)也会极大的改变。

2.2相关系数

    为了解决上述问题，我们引入子图像和模板图像中的灰度平均值。上述归一化函数可以改写为：

   子图像和参考图像平均值分别定义为：

    其中，K表示参考图像中像素的数量。在统计学中，上述表达式被称为相关系数。然而，与统计学中作为全局测量方法不同的是，CL（r,s）被定义是一个局部的测量函数。CL（r,s）的值域范围为[-1,1]。当CL(r,s)等于1时，表明在平移位置（r,s），参考图像和子图像达到最大相似度；相反，当CL(r,s)等于0时，表明在平移位置（r,s），参考图像和子图像完全不匹配。表达式：

表示的是K倍模板图像R中像素灰度值方差，这是一个常量，只需要计算一次即可。由于方差存在如下关系：（具体可参考http://imagej.net/Integral_Image_Filters）

        具体推到如下所示：

     因此，可以SR改写为

  带入到CL(r,s)，可得：

  由此，可得到一个高效的计算相关系数的方法。

3.代码实现

int Corre_Coefficient(IMAGE_TYPE *OriginalBMP_img,IMAGE_TYPE *TempBMP_img){DWORD Source_width,Source_height,Template_width,Template_height;WORD  Source_biBitCount,Template_biBitCount;T_U8 *Source_dst,*Template_dst,*source_img,*template_img,*result,*result_dst;T_U32 Source_linebyte,Template_linebyte,K,SumR = 0,SumR2 = 0,SumI,SumI2,SumIR,AI,AR;T_U32 i,j,r,s,x,y;double MeanR = 0,MeanSR = 0,Crs,MeanI = 0,Max = -999,tempValue;BITMAPFILEHEADER Source_bf,Template_bf,Result_bf;BITMAPINFOHEADER Source_bi,Template_bi,Result_bi;FILE *CorrectionBMP_fp = fopen("Corre_Coefft.bmp","wb");RGBQUAD colortable[256];RGBQUAD* pColorTable1 = colortable;if(NULL == CorrectionBMP_fp){printf("Can't open CorrelationCoefficient.bmp.\n");return -1;} for (i = 0; i<256; ++i)      {          colortable[i].rgbBlue = i;          colortable[i].rgbGreen = i;          colortable[i].rgbRed = i;          colortable[i].rgbReserved = 0;      }  memset(&Source_bf, 0, sizeof(Source_bf));memset(&Source_bi, 0, sizeof(Source_bi));source_img = OriginalBMP_img;memcpy(&Source_bf,source_img,14);memcpy(&Source_bi,&source_img[14],40);Source_height = Source_bi.biHeight;Source_width  = Source_bi.biWidth;Source_biBitCount = Source_bi.biBitCount;//每一个像素由24 bits表示，即RGB分量每一个分量用8 bits表示Source_linebyte = (Source_width * Source_bi.biBitCount / 8 + 3) / 4 * 4; //Source_linebyte = WIDTHBYTES(Source_width*Source_bi.biBitCount);Source_dst = source_img+54+1024;memset(&Template_bf, 0, sizeof(Template_bf));memset(&Template_bi, 0, sizeof(Template_bi));template_img = TempBMP_img;memcpy(&Template_bf,TempBMP_img,14);memcpy(&Template_bi,&TempBMP_img[14],40);Template_height = Template_bi.biHeight;Template_width  = Template_bi.biWidth;Template_biBitCount = Template_bi.biBitCount;//每一个像素由24 bits表示，即RGB分量每一个分量用8 bits表示Template_linebyte = (Template_width * Template_biBitCount / 8 + 3) / 4 * 4; //Template_linebyte = WIDTHBYTES(Template_width*Template_bi.biBitCount);Template_dst = template_img+54+1024;result = (T_U8 *)malloc(Source_height*Source_linebyte+BMPHEADSIZE+1024);memset(&Result_bf, 0, sizeof(Source_bf));memset(&Result_bi, 0, sizeof(Source_bi));memcpy(&Result_bf,source_img,14);memcpy(&Result_bi,&source_img[14],40);memcpy(pColorTable1,&source_img[54],sizeof(RGBQUAD)*256);memcpy(result,source_img,Source_height*Source_linebyte+BMPHEADSIZE+1024);memset(result+54+1024,0,Source_height*Source_linebyte);result_dst = result + 54+1024;K = Template_height*Template_width;for (i = 0;i < Template_height;i++){for (j = 0;j < Template_width;j++){SumR = SumR+Template_dst[i*Template_linebyte+j];SumR2 = SumR2+Template_dst[i*Template_linebyte+j]*Template_dst[i*Template_linebyte+j];}}MeanR = (double)SumR/K;MeanSR = sqrt(SumR2-(double)(SumR*SumR)/K);for (r = 0; r < Source_height-Template_height;r++){for (s = 0;s < Source_width - Template_width;s++){SumI = 0;SumI2 = 0;SumIR = 0;for (i = 0;i < Template_height;i++){for (j = 0;j < Template_width;j++){AI = Source_dst[(r+i)*Source_linebyte+s+j];AR = Template_dst[i*Template_linebyte+j];SumI = SumI+AI;SumI2 = SumI2+AI*AI;SumIR = SumIR+AI*AR;}}tempValue = sqrt(SumI2-(double)(SumI*SumI)/K);Crs = (SumIR-MeanR*SumI)/(tempValue*MeanSR);if(Crs > Max ){Max = Crs;x = s;y = r;}}}for (i = y; i < y+Template_height;i++){for (j = x;j < x+Template_width;j++){result_dst[i*Source_linebyte+j] = Source_dst[i*Source_linebyte+j];}}fwrite(result,Source_height*Source_linebyte+BMPHEADSIZE+1024,1,CorrectionBMP_fp);fclose(CorrectionBMP_fp);CorrectionBMP_fp = NULL;free(result);return 0;}

4.匹配效果

左上角为模板图像，左下角为在原始图像中匹配到的图像，右侧图像为原始图像