基于模板匹配的目标检测（2）-opencv模板匹配的加速原理

搜了很多博文，几乎没有博主去分析opencv的matchtemplate接口的底层加速原理。因此，我们选择相对较为简单的TM_CCORR_NORMED归一化相关匹配，对其加速方案进行叙述。本文相关源码和素材均以上传，可免费下载：opencv模板匹配加速原理源码和测试图像

涉及知识点

• FFT加速相关计算
• 积分图

常规思路

在没有去阅读opencv的相关源码之前，我一直以为opencv的模板匹配是通过SIMD指令、OpenCL、GPU等并行计算工具进行深度优化加速的。通过阅读源码才发现，opencv对模板匹配的加速方案，远比我们想象要来得巧妙。

原理加速

卷积定理指出：时域的卷积对应着频域的点乘（区分深度学习中的卷积，严格意义上来说，深度学习中的卷积应该叫相关），时域的相关对应着频域的共轭点乘。
卷积操作：卷积核进行转置、滑窗求点积。
相关操作：直接使用核进行滑窗求点积。
有了这个理论基础，我们就来分析归一化相关匹配的公式：

对于分子部分，就是一个普通的相关操作，直接使用FFT进行加速计算。
对于分母部分，是需要计算模板对应区域的平方和进行开根号，可以理解为计算图像任意矩形区域的平方和。这部分，opencv采用的方案是积分图。

积分图

图像积分图是计算机视觉中一个重要的数据结构，主要用于快速计算图像中任意矩形区域的像素和。
图像积分图是与输入图像的大小相同，它的每个元素存储了从图像左上角原点到对应位置的矩形区域内所有像素的灰度值之和。换句话说，对于图像中的每个像素(i, j)，其积分图像值I(i, j)等于该像素及其左上角所有像素的灰度值之和。
从直观来说，一张图像就是一个矩形，这个矩形中每个像素点的积分值，就是以图像左上角像素点为左上角顶点，以该像素点为右下角顶点的矩形中包含的所有元素之和。如下图所示，点(x,y)处的积分值为矩形区域中所有像素之和（包括该点）。

在得到积分图之后，如果需要快速计算图像上任意矩形区域的和值，那就可以使用以下方法。
设图像为i，对应的积分图为I，计算矩形区域（A、B、C、D）的和值，计算方法如下：

可以简单的通过4个点的积分图的值进行线性运算，计算出任意矩形区域的和值，比使用两层循环进行遍历简洁、省时。

加速原理验证

为了验证加速原理是否正确，使用opencv模拟实现归一化相关匹配。代码如下：
首先，我们使用循环的滑窗的方式进行实现，作为比较基准：

/*
* 使用朴素的循环实现
*
*/
cv::Mat principleCORR(cv::Mat& input, cv::Mat t)
{
	//计算模板的平方和
	float tsqsum = 0.0f;
	for (int r = 0; r < t.rows; r++)
	{
		unsigned char* row = t.ptr(r);
		for (int c = 0; c < t.cols; c++)
		{
			tsqsum += row[c] * row[c];
		}
	}
	//开方
	tsqsum = std::sqrtf(tsqsum);
	//相似度矩阵缓存申请
	cv::Mat ncc = cv::Mat::zeros(input.rows - t.rows + 1, input.cols - t.cols + 1, CV_32FC1);
	//滑窗循环
	for (int i = 0; i < input.rows - t.rows + 1; i++)
		for (int j = 0; j < input.cols - t.cols + 1; j++)
		{
			//点积
			double dot = 0.0f;
			//平方和
			float sqsum = 0.0f;
			//模板循环
			for (int m = 0; m < t.rows; m++)
			{
				//待匹配图像的行首地址
				unsigned char* img_row = input.ptr(i + m);
				//模板图像的行首地址
				unsigned char* t_row = t.ptr(m);
				for (int n = 0; n < t.cols; n++)
				{
					float t_value = t_row[n];
					float img_value = img_row[n+j];
					//分子公式（相乘）
					dot += t_value * img_value;
					//分母公式（平方和）
					sqsum += img_value * img_value;
				}
			}
			//相似度值
			double nc = dot / (tsqsum * std::sqrtf(sqsum));
			ncc.at<float>(i, j) = nc;
		}

	return ncc;
}

根据上一章节的加速原理分析，使用opencv的进行仿真模拟实现。

/*
* 模拟opencv的加速方案模拟实现
*
*/
cv::Mat fastCORR(cv::Mat& input, cv::Mat t)
{
	/*
	* 1、利用时域相关频域复共轭相乘的定理，完成分子的计算
	*/
	//将模板图像进行补0填充
	cv::Mat fill_template;
	cv::copyMakeBorder(t, fill_template, 0, input.rows - t.rows, 
		0, input.cols - t.cols, cv::BORDER_CONSTANT, cv::Scalar::all(0));
	//傅里叶变换
	cv::Mat input_fft, template_fft;
	input.convertTo(input_fft, CV_32F);
	cv::dft(input_fft, input_fft, cv::DFT_COMPLEX_OUTPUT);
	fill_template.convertTo(template_fft, CV_32F);
	cv::dft(template_fft, template_fft, cv::DFT_COMPLEX_OUTPUT);
	//频域共轭相乘
	cv::Mat dot_fft;
	cv::mulSpectrums(input_fft, template_fft, dot_fft, 0, true);
	//逆傅里叶变换
	cv::idft(dot_fft, dot_fft, cv::DFT_REAL_OUTPUT | cv::DFT_SCALE);

	/*
	* 2、利用积分图完成分母的计算，同时完成输出
	*/
	//计算积分图
	cv::Mat sq, ssq;//和、平方和
	cv::integral(input, sq, ssq, CV_64F);
	//提前完成模板的平方和、开根号
	float tsqsum = 0.0f;
	for (int r = 0; r < t.rows; r++)
	{
		unsigned char* row = t.ptr(r);
		for (int c = 0; c < t.cols; c++)
		{
			tsqsum += row[c] * row[c];
		}
	}
	//开方
	tsqsum = std::sqrtf(tsqsum);

	//循环遍历完公式计算，并完成输出
	cv::Mat ncc = cv::Mat::zeros(input.rows - t.rows + 1, input.cols - t.cols + 1, CV_32FC1);
	for (int r = 0; r < input.rows - t.rows + 1; r++)
		for (int c = 0; c < input.cols - t.cols + 1; c++)
		{
			//取分子的值
			float dot = dot_fft.at<float>(r, c);
			//利用积分图完成分母计算
			double A = ssq.at<double>(r, c);
			double B = ssq.at<double>(r, c + t.cols);
			double C = ssq.at<double>(r + t.rows, c);
			double D = ssq.at<double>(r + t.rows, c + t.cols);
			float sqsum = std::sqrtf(D + A - B - C);
			ncc.at<float>(r, c) = dot / (tsqsum * sqsum);
		}
		
	return ncc;
}

同时，在为了验证加速原理是否正确，使用opencv的matchtemplate进行基准对比。main函数实现如下：

int main()
{
	//以灰度模式读取待匹配图像
	cv::Mat image = cv::imread("./待匹配图像.png", cv::IMREAD_GRAYSCALE);
	//以灰度模式读取模板图像
	cv::Mat template_img = cv::imread("./模板图像.png", cv::IMREAD_GRAYSCALE);

	//opencv API
	cv::Mat cv_template_ncc;
	auto start = cv::getTickCount();
	cv::matchTemplate(image, template_img, cv_template_ncc, cv::TM_CCORR_NORMED);
	auto end = cv::getTickCount();
	auto cost_time = (end - start) / cv::getTickFrequency() * 1000;
	printf("opencv matchTemplate 耗时：%f ms\n", cost_time);
	cv::normalize(cv_template_ncc, cv_template_ncc, 1.0, 0.0, cv::NORM_MINMAX);
	cv::imshow("opencv-matchTemplate实现", cv_template_ncc);
	
	//模拟opencv的加速方案实现
	cv::Mat fast_ncc;
	start = cv::getTickCount();
	fast_ncc = fastCORR(image, template_img);
	end = cv::getTickCount();
	cost_time = (end - start) / cv::getTickFrequency() * 1000;
	printf("模拟opencv的加速方案实现耗时：%f ms\n", cost_time);
	cv::normalize(fast_ncc, fast_ncc, 1.0, 0.0, cv::NORM_MINMAX);
	cv::imshow("模拟opencv加速方案实现", fast_ncc);


	//测试朴素实现的归一化相关匹配
	start = cv::getTickCount();
	cv::Mat princi_ncc = principleCORR(image, template_img);
	end = cv::getTickCount();
	cost_time = (end - start) / cv::getTickFrequency() * 1000;
	printf("朴素循环实现耗时：%f ms\n", cost_time);
	cv::normalize(princi_ncc, princi_ncc, 1.0, 0.0, cv::NORM_MINMAX);
	cv::imshow("朴素循环实现", princi_ncc);

	cv::imshow("模板图像", template_img);
	cv::imshow("待匹配图像", image);

	cv::Point pt_max, pt_min;
	double dmin, dmax;
	cv::minMaxLoc(princi_ncc, &dmin, &dmax, &pt_min, &pt_max);
	cv::Mat draw;
	cv::cvtColor(image, draw, cv::COLOR_GRAY2BGR);
	cv::rectangle(draw, cv::Rect(pt_max.x, pt_max.y, template_img.cols, template_img.rows),
		cv::Scalar(0, 255, 0), 2);
	cv::imshow("匹配结果", draw);
	cv::waitKey();
	return 0;
}

对比验证结果如下：

实验分析

• 从相似度矩阵分析，加速原理是正确无误的。
• 从耗时分析，朴素循环耗时最长，opencv的耗时最短，模拟加速原理居中。耗时的绝对值不具备参考意义，但是，相对值是可信的。加速原理确实比纯粹的循环实现节省时间，提速非常明显。但是，依然比不上opencv自带的匹配函数，主要是因为，opencv在实现中，综合运用了SIMD、分块DFT等手段。
• 积分图是常用到的图像处理算法加速手段，在日常算法开发中应合理应用。
• 图像处理算法的加速，应该结合算法原理优化和硬件加速。