1 理论介绍
模板匹配是在一幅图像中寻找一个特定目标的方法之一，这种方法的原理非常简单，遍历图像中的每一个可能的位置，比较各处与模板是否“相似”，当相似度足够高时，就认为找到了我们的目标。OpenCV提供了6种模板匹配算法：
- 平方差匹配法CV_TM_SQDIFF
- 归一化平方差匹配法CV_TM_SQDIFF_NORMED
- 相关匹配法CV_TM_CCORR
- 归一化相关匹配法CV_TM_CCORR_NORMED
- 相关系数匹配法CV_TM_CCOEFF
- 归一化相关系数匹配法CV_TM_CCOEFF_NORMED
用T表示模板图像，I表示待匹配图像，切模板图像的宽为w高为h，用R表示匹配结果，匹配过程如下图所示：
上述6中匹配方法可用以下公式进行描述：
2 示例代码
下面给出方法6的python代码
1 import numpy as np 2 import cv2 3 4 def EM(pModel, width, height): 5 sum = np.double(0.0) 6 for i in range(0,height): 7 for j in range(0,width): 8 sum += pModel[i][j] 9 return sum 10 11 def EM2(pModel, width, height): 12 sum = np.double(0.0) 13 for i in range(0,height): 14 for j in range(0,width): 15 sum += pModel[i][j]*1.0*pModel[i][j] 16 return sum 17 18 def EI(pToSearch, l, h, u, v, pModel, width, height): 19 sum = np.double(0.0) 20 roi = pToSearch[v:v+height, u:u+width] 21 for i in range(0,height): 22 for j in range(0,width): 23 sum += roi[i][j] 24 return sum 25 26 def EI2(pToSearch, l, h, u, v, pModel, width, height): 27 sum = np.double(0.0) 28 roi = pToSearch[v:v+height, u:u+width] 29 for i in range(0,height): 30 for j in range(0,width): 31 sum += roi[i][j]*1.0*roi[i][j] 32 return sum 33 34 def EIM(pToSearch, l, h, u, v, pModel, width, height): 35 sum = np.double(0.0) 36 roi = pToSearch[v:v+height, u:u+width] 37 for i in range(0,height): 38 for j in range(0,width): 39 sum += pModel[i][j]*1.0*roi[i][j] 40 return sum 41 42 def Match(pToSearch, l, h, pModel, width, height): 43 uMax = l-width 44 vMax = h-height 45 N = width*height 46 len = (uMax+1)*(vMax+1) 47 MatchRec = [0.0 for x in range(0, len)] 48 k = 0 49 50 M = EM(pModel,width,height) 51 M2 = EM2(pModel,width,height) 52 for p in range(0, uMax+1): 53 for q in range(0, vMax+1): 54 I = EI(pToSearch,l,h,p,q,pModel,width,height) 55 I2 = EI2(pToSearch,l,h,p,q,pModel,width,height) 56 IM = EIM(pToSearch,l,h,p,q,pModel,width,height) 57 58 numerator=(N*IM-I*M)*(N*IM-I*M) 59 denominator=(N*I2-I*I)*(N*M2-M*M) 60 61 ret = numerator/denominator 62 MatchRec[k]=ret 63 k+=1 64 65 val = 0 66 k = 0 67 x = y = 0 68 for p in range(0, uMax+1): 69 for q in range(0, vMax+1): 70 if MatchRec[k] > val: 71 val = MatchRec[k] 72 x = p 73 y = q 74 k+=1 75 print "val: %f"%val 76 return (x, y) 77 78 def main(): 79 img = cv2.imread('niu.jpg', cv2.IMREAD_GRAYSCALE) 80 temp = cv2.imread('temp.png', cv2.IMREAD_GRAYSCALE) 81 82 print temp.shape 83 imgHt, imgWd = img.shape 84 tempHt, tempWd = temp.shape 85 #print EM(temp, tempWd, tempHt) 86 (x, y) = Match(img, imgWd, imgHt, temp, tempWd, tempHt) 87 cv2.rectangle(img, (x, y), (x+tempWd, y+tempHt), (0,0,0), 2) 88 cv2.imshow("temp", temp) 89 cv2.imshow("result", img) 90 cv2.waitKey(0) 91 cv2.destroyAllWindows() 92 93 if __name__ == '__main__': 94 main()
代码58行中的N就是公式(6)中的w*h，由于python代码运行速度比较慢，代码的58、59行相当于对公式(6)的分子分母都进行了平方操作，并且分子分母都乘以了N方，以减小计算量，所以代码61行的ret相当于公式(6)中的R(x,y)的平方，
为了更快的进行算法验证，用上述代码进行验证时请尽量选用较小的匹配图像及模板图像，下图显示了我的匹配结果（待匹配图像295x184模板69x46用了十几分钟）：
3 OpenCV源码
较新版本的OpenCV库中的模板匹配已经进行了较多的算法改进，直接看新版本中的算法需要了解很多相关理论知识，所以我们结合OpenCV0.9.5的源码进行讲解，该版本的源码基本上是C风格代码更容易进行理解（如果要对
OpenCV源码进行研究，建议用该版本进行入门），仍以归一化相关系数匹配法为例进行分析。
1 /* 2 * pImage: 待匹配图像 3 * image: 待匹配图像宽（width*depth并已4字节对齐） 4 * roiSize: 待匹配图像尺寸 5 * pTemplate: 模板图像 6 * templStep: 模板图像宽 7 * templSize: 模板图像尺寸 8 * pResult: 匹配结果 9 * resultStep: 匹配结果宽 10 * pBuffer: 中间结果数据缓存 11 */ 12 IPCVAPI_IMPL( CvStatus, icvMatchTemplate_CoeffNormed_32f_C1R, 13 (const float *pImage, int imageStep, CvSize roiSize, 14 const float *pTemplate, int templStep, CvSize templSize, 15 float *pResult, int resultStep, void *pBuffer) ) 16 { 17 float *imgBuf = 0; // 待匹配图像相关数据 18 float *templBuf = 0; // 模板图像数据 19 double *sumBuf = 0; // 待匹配图像遍历块单行和 20 double *sqsumBuf = 0; // 待匹配图像遍历块单行平方和 21 double *resNum = 0; // 模板图像和待匹配图像遍历块内积 22 double *resDenom = 0; // 待匹配图像遍历块累加和及待匹配图像遍历块平方累加和 23 double templCoeff = 0; // 模板图像均分差倒数 24 double templSum = 0; // 模板图像累加和 25 26 int winLen = templSize.width * templSize.height; 27 double winCoeff = 1. / (winLen + DBL_EPSILON); // + DBL_EPSILON 加一个小整数防止分母为零 28 29 CvSize resultSize = cvSize( roiSize.width - templSize.width + 1, 30 roiSize.height - templSize.height + 1 ); 31 int x, y; 32 33 // 计算并为imgBuf、templBuf、sumBuf、sqsumBuf、resNum、resDenom分配存储空间 34 CvStatus result = icvMatchTemplateEntry( pImage, imageStep, roiSize, 35 pTemplate, templStep, templSize, 36 pResult, resultStep, pBuffer, 37 cv32f, 1, 1, 38 (void **) &imgBuf, (void **) &templBuf, 39 (void **) &sumBuf, (void **) &sqsumBuf, 40 (void **) &resNum, (void **) &resDenom ); 41 42 if( result != CV_OK ) 43 return result; 44 45 imageStep /= sizeof_float; 46 templStep /= sizeof_float; 47 resultStep /= sizeof_float; 48 49 /* calc common statistics for template and image */ 50 { 51 const float *rowPtr = (const float *) imgBuf; 52 double templSqsum = icvCrossCorr_32f_C1( templBuf, templBuf, winLen ); // 模板图像平方累加和 53 54 templSum = icvSumPixels_32f_C1( templBuf, winLen ); // 模板图像累加和 55 templCoeff = (double) templSqsum - ((double) templSum) * templSum * winCoeff; // 模板图像均方差的平方 56 templCoeff = icvInvSqrt64d( fabs( templCoeff ) + FLT_EPSILON ); // 模板图像均方差倒数 57 58 for( y = 0; y < roiSize.height; y++, rowPtr += templSize.width ) 59 { 60 sumBuf[y] = icvSumPixels_32f_C1( rowPtr, templSize.width ); // 待匹配图像按模板图像宽度求每行之和（遍历位置第一列） 61 sqsumBuf[y] = icvCrossCorr_32f_C1( rowPtr, rowPtr, templSize.width ); // 待匹配图像按模板图像宽度求每行平方之和（遍历位置第一列） 62 } 63 } 64 65 /* main loop - through x coordinate of the result */ 66 for( x = 0; x < resultSize.width; x++ ) 67 { 68 double sum = 0; 69 double sqsum = 0; 70 float *imgPtr = imgBuf + x; // 待匹配图像起始位置 71 72 /* update sums and image band buffer */ // 如果不是第1列需重新更新sumBuf，更新后sumBuf为遍历位置第x列每行之和（行宽为模板图像宽） 73 if( x > 0 ) 74 { 75 const float *src = pImage + x + templSize.width - 1; 76 float *dst = imgPtr - 1; 77 float out_val = dst[0]; 78 79 dst += templSize.width; 80 81 for( y = 0; y < roiSize.height; y++, src += imageStep, dst += templSize.width ) 82 { 83 float in_val = src[0]; 84 85 sumBuf[y] += in_val - out_val; 86 sqsumBuf[y] += (in_val - out_val) * (in_val + out_val); 87 out_val = dst[0]; 88 dst[0] = (float) in_val; 89 } 90 } 91 92 for( y = 0; y < templSize.height; y++ ) // 求遍历位置第x列，第1行处遍历块累加和sum及平方累加和sqsum 93 { 94 sum += sumBuf[y]; 95 sqsum += sqsumBuf[y]; 96 } 97 98 for( y = 0; y < resultSize.height; y++, imgPtr += templSize.width ) 99 { 100 double res = icvCrossCorr_32f_C1( imgPtr, templBuf, winLen ); // 求模板图像和待匹配图像y行x列处遍历块的内积 101 102 if( y > 0 ) // 如果不是第1行需更新遍历块累加和sum及平方累加和sqsum 103 { 104 sum -= sumBuf[y - 1]; 105 sum += sumBuf[y + templSize.height - 1]; 106 sqsum -= sqsumBuf[y - 1]; 107 sqsum += sqsumBuf[y + templSize.height - 1]; 108 } 109 resNum[y] = res; 110 resDenom[y] = sum; 111 resDenom[y + resultSize.height] = sqsum; 112 } 113 114 for( y = 0; y < resultSize.height; y++ ) 115 { 116 double sum = ((double) resDenom[y]); 117 double wsum = winCoeff * sum; 118 double res = ((double) resNum[y]) - wsum * templSum; 119 double nrm_s = ((double) resDenom[y + resultSize.height]) - wsum * sum; 120 121 res *= templCoeff * icvInvSqrt64d( fabs( nrm_s ) + FLT_EPSILON ); 122 pResult[x + y * resultStep] = (float) res; 123 } 124 } 125 126 return CV_OK; 127 }
以上代码是归一化相关系数法核心函数icvMatchTemplate_CoeffNormed_32f_C1R的源码，我已经在源码中进行了详细的注释，读者需自己再进行理解，需要进一步说明的是：
代码118行res就是计算公式（6）的分子部分，代码56行templCoeff就是计算公式(6)分母的左半部分，代码121行icvInvSqrt64d函数就是在计算公式(6)分母的右半部分，该行res的最终结果正是公式(6)中的R(x,y)。
4 结束语
OpenCV0.9.5源码下载：http://download.csdn.net/detail/weiwei22844/9547820