图像相似度测量与模板匹配总结

本文主要总结了进行目标跟踪、检测中经常使用到的图像相似度测量和模板匹配方法，并给出了具体的基于OpenCV的代码实现。

模板匹配是一种在源图像中寻找与图像patch最相似的技术，常常用来进行目标的识别、跟踪与检测。其中最相似肯定是基于某种相似度准则来讲的，也就是需要进行相似度的测量。另外，寻找就需要在图像上进行逐行、逐列的patch窗口扫描，当然也不一定需要逐行逐列的扫描，当几个像素的误差比计算速度来的不重要时就可以设置扫描的行列步进值，以加快扫描和计算的时间消耗。下面就对相似度测量和模板匹配进行介绍（所有的图像都假定是灰度图）。

图像相似度测量

直方图方法

方法描述：有两幅图像patch(当然也可是整幅图像)，分别计算两幅图像的直方图，并将直方图进行归一化，然后按照某种距离度量的标准进行相似度的测量。

方法的思想：基于简单的向量相似度来对图像相似度进行度量。

优点：直方图能够很好的归一化，比如256个bin条，那么即使是不同分辨率的图像都可以直接通过其直方图来计算相似度，计算量适中。比较适合描述难以自动分割的图像。

缺点：直方图反应的是图像灰度值得概率分布，并没有图像的空间位置信息在里面，因此，常常出现误判；从信息论来讲，通过直方图转换，信息丢失量较大，因此单一的通过直方图进行匹配显得有点力不从心。

基于opencv的实现，我把它封装为函数（输入输出都是灰度图像，且设定的灰度级为8，及0-255），有需要的可以直接拿去使用

<span style="font-size:18px;"><span style="font-size:18px;">double getHistSimilarity(const Mat& I1, const Mat& I2)
{
	int histSize = 256;
	float range[] = {0,256};
	const float* histRange = {range};
	bool uniform = true;
	bool accumulate = false;

	Mat hist1,hist2;

	calcHist(&I1,1,0,Mat(),hist1,1,&histSize,&histRange,uniform,accumulate);
	normalize(hist1,hist1,0,1,NORM_MINMAX,-1,Mat());

	calcHist(&I2,1,0,Mat(),hist2,1,&histSize,&histRange,uniform,accumulate);
	normalize(hist2,hist2,0,1,NORM_MINMAX,-1,Mat());

	return compareHist(hist1, hist2, CV_COMP_CORREL);

}</span></span>

/* Histogram comparison methods */
enum
{
CV_COMP_CORREL =0,
CV_COMP_CHISQR =1,
CV_COMP_INTERSECT =2,
CV_COMP_BHATTACHARYYA =3,
CV_COMP_HELLINGER =CV_COMP_BHATTACHARYYA
};
分别是：相关度，卡方，相交系数和巴氏距离，其计算公式如下图所示：

直方图方法的一些改进的出发点大致是为了改善直方图无法利用空间位置信息的弱点，比如FragTrack算法，将图像再进行分割成横纵的patch，然后对于每一个patch搜索与之匹配的直方图，来计算两幅图像的相似度，从而融入了直方图对应的位置信息，但计算效率肯定不高。

矩阵分解的方法

方法描述：将图像patch做矩阵分解，比如SVD奇异值分解和NMF非负矩阵分解等，然后再做相似度的计算。

方法思想：因为图像本身来讲就是一个矩阵，可以依靠矩阵分解获取一些更加鲁棒的特征来对图像进行相似度的计算。

基于SVD分解的方法优点：奇异值的稳定性，比例不变性，旋转不变性和压缩性。即奇异值分解是基于整体的表示，不但具有正交变换、旋转、位移、镜像映射等代数和几何上的不变性,而且具有良好的稳定性和抗噪性，广泛应用于模式识别与图像分析中。对图像进行奇异值分解的目的是得到唯一、稳定的特征描述，降低特征空间的维度，提高抗干扰能力。

基于SVD分解的方法缺点是：奇异值分解得到的奇异矢量中有负数存在，不能很好的解释其物理意义。

基于特征点方法

基于峰值信噪比（PSNR）的方法

当我们想检查压缩视频带来的细微差异的时候，就需要构建一个能够逐帧比较差视频差异的系统。最
常用的比较算法是PSNR( Peak signal-to-noise ratio)。这是个使用“局部均值误差”来判断差异的最简单的方法，假设有这两幅图像：I1和I2，它们的行列数分别是i，j，有c个通道。每个像素的每个通道的值占用一个字节，值域[0,255]。注意当两幅图像的相同的话，MSE的值会变成0。这样会导致PSNR的公式会除以0而变得没有意义。所以我们需要单独的处理这样的特殊情况。此外由于像素的动态范围很广，在处理时会使用对数变换来缩小范围。

基于OpenCV的PSNR方法实现，同样将其封装为函数，可以直接拿来调用：

<span style="font-size:18px;">double getPSNR(const Mat& I1, const Mat& I2)
{
	Mat s1;
	absdiff(I1, I2, s1);       // |I1 - I2|
	s1.convertTo(s1, CV_32F);  // cannot make a square on 8 bits
	s1 = s1.mul(s1);           // |I1 - I2|^2

	Scalar s = sum(s1);         // sum elements per channel

	double sse = s.val[0] + s.val[1] + s.val[2]; // sum channels

	if( sse <= 1e-10) // for small values return zero
		return 0;
	else
	{
		double  mse =sse /(double)(I1.channels() * I1.total());
		double psnr = 10.0*log10((255*255)/mse);
		return psnr;
	}
}</span>

基于结构相似性（SSIM,structural similarity (SSIM) index measurement）的方法

基于OpenCV的SSIM方法实现，同样将其封装为函数，可以直接拿来调用：

<span style="font-size:18px;">Scalar getMSSIM( const Mat& i1, const Mat& i2)
{
	const double C1 = 6.5025, C2 = 58.5225;
	/***************************** INITS **********************************/
	int d     = CV_32F;

	Mat I1, I2;
	i1.convertTo(I1, d);           // cannot calculate on one byte large values
	i2.convertTo(I2, d);

	Mat I2_2   = I2.mul(I2);        // I2^2
	Mat I1_2   = I1.mul(I1);        // I1^2
	Mat I1_I2  = I1.mul(I2);        // I1 * I2

	/*************************** END INITS **********************************/

	Mat mu1, mu2;   // PRELIMINARY COMPUTING
	GaussianBlur(I1, mu1, Size(11, 11), 1.5);
	GaussianBlur(I2, mu2, Size(11, 11), 1.5);

	Mat mu1_2   =   mu1.mul(mu1);
	Mat mu2_2   =   mu2.mul(mu2);
	Mat mu1_mu2 =   mu1.mul(mu2);

	Mat sigma1_2, sigma2_2, sigma12;

	GaussianBlur(I1_2, sigma1_2, Size(11, 11), 1.5);
	sigma1_2 -= mu1_2;

	GaussianBlur(I2_2, sigma2_2, Size(11, 11), 1.5);
	sigma2_2 -= mu2_2;

	GaussianBlur(I1_I2, sigma12, Size(11, 11), 1.5);
	sigma12 -= mu1_mu2;

	///////////////////////////////// FORMULA ////////////////////////////////
	Mat t1, t2, t3;

	t1 = 2 * mu1_mu2 + C1;
	t2 = 2 * sigma12 + C2;
	t3 = t1.mul(t2);              // t3 = ((2*mu1_mu2 + C1).*(2*sigma12 + C2))

	t1 = mu1_2 + mu2_2 + C1;
	t2 = sigma1_2 + sigma2_2 + C2;
	t1 = t1.mul(t2);               // t1 =((mu1_2 + mu2_2 + C1).*(sigma1_2 + sigma2_2 + C2))

	Mat ssim_map;
	divide(t3, t1, ssim_map);      // ssim_map =  t3./t1;

	Scalar mssim = mean( ssim_map ); // mssim = average of ssim map
	return mssim;
}
</span>

图像模板匹配

一般而言，源图像与模板图像patch尺寸一样的话，可以直接使用上面介绍的图像相似度测量的方法；如果源图像与模板图像尺寸不一样，通常需要进行滑动匹配窗口，扫面个整幅图像获得最好的匹配patch。

在OpenCV中对应的函数为：matchTemplate()：函数功能是在输入图像中滑动窗口寻找各个位置与模板图像patch的相似度。相似度的评价标准（匹配方法）有：CV_TM_SQDIFF平方差匹配法（相似度越高，值越小），CV_TM_CCORR相关匹配法（采用乘法操作，相似度越高值越大），CV_TM_CCOEFF相关系数匹配法（1表示最好的匹配，-1表示最差的匹配）。

CV_TM_SQDIFF计算公式：

CV_TM_CCORR计算公式：

有一种新的用来计算相似度或者进行距离度量的方法：EMD，Earth Mover's Distances

EMD is defined as the minimal cost that must be paid to transform one histograminto the other, where there is a “ground distance” between the basic featuresthat are aggregated into the histogram。

光线变化能引起图像颜色值的漂移，尽管漂移没有改变颜色直方图的形状，但漂移引起了颜色值位置的变化，从而可能导致匹配策略失效。而EMD是一种度量准则，度量怎样将一个直方图转变为另一个直方图的形状，包括移动直方图的部分（或全部）到一个新的位置，可以在任意维度的直方图上进行这种度量。

2015-7-24