OpenCV 实现了图像平移模板匹配的功能,封装在函数接口
matchTemplate
中,本文解析该功能的实现源码。
matchTemplate
速度很快,核心提速在于使用了卷积加速和累加和技巧opencv\sources\modules\imgproc\src\templmatch.cpp
CV_TM_CCOEFF_NORMED
)运算过程最为复杂,此处以该损失函数为例解读源码,其余函数可以以此类推matchTemplate
函数中实现,源码void cv::matchTemplate( InputArray _img, InputArray _templ, OutputArray _result, int method, InputArray _mask )
{
CV_INSTRUMENT_REGION();
int type = _img.type(), depth = CV_MAT_DEPTH(type), cn = CV_MAT_CN(type);
CV_Assert( CV_TM_SQDIFF <= method && method <= CV_TM_CCOEFF_NORMED );
CV_Assert( (depth == CV_8U || depth == CV_32F) && type == _templ.type() && _img.dims() <= 2 );
if (!_mask.empty())
{
cv::matchTemplateMask(_img, _templ, _result, method, _mask);
return;
}
bool needswap = _img.size().height < _templ.size().height || _img.size().width < _templ.size().width;
if (needswap)
{
CV_Assert(_img.size().height <= _templ.size().height && _img.size().width <= _templ.size().width);
}
CV_OCL_RUN(_img.dims() <= 2 && _result.isUMat(),
(!needswap ? ocl_matchTemplate(_img, _templ, _result, method) : ocl_matchTemplate(_templ, _img, _result, method)))
Mat img = _img.getMat(), templ = _templ.getMat();
if (needswap)
std::swap(img, templ);
Size corrSize(img.cols - templ.cols + 1, img.rows - templ.rows + 1);
_result.create(corrSize, CV_32F);
Mat result = _result.getMat();
CV_IPP_RUN_FAST(ipp_matchTemplate(img, templ, result, method))
crossCorr( img, templ, result, Point(0,0), 0, 0);
common_matchTemplate(img, templ, result, method, cn);
}
result
为卷积图结果,计算时用了 CV_IPP_RUN_FAST
加速ipp全称:Intel® Integrated Performance Primitives英特尔®集成性能原语(Intel®IPP)是一种软件库,提供了广泛的功能,包括通用信号和图像处理、计算机视觉、数据压缩和字符串操作。 如果在英特尔的处理器上使用,OpenCV就会自动使用一种免费的英特尔集成性能原语库(IPP)的子集,IPP 8.x(IPPICV)。IPPICV可以在编译阶段链接到OpenCV,这样一来,会替代相应的低级优化的C语言代码(在cmake中设置WITH_IPP=ON/OFF开启或者关闭这一功能,默认情况为开启)。
result
变量中common_matchTemplate
函数static void common_matchTemplate( Mat& img, Mat& templ, Mat& result, int method, int cn )
{
if( method == CV_TM_CCORR )
return;
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;
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 = templSdv[0]*templSdv[0] + templSdv[1]*templSdv[1] + templSdv[2]*templSdv[2] + templSdv[3]*templSdv[3];
if( templNorm < DBL_EPSILON && method == CV_TM_CCOEFF_NORMED )
{
result = Scalar::all(1);
return;
}
templSum2 = templNorm + templMean[0]*templMean[0] + templMean[1]*templMean[1] + templMean[2]*templMean[2] + templMean[3]*templMean[3];
if( numType != 1 )
{
templMean = Scalar::all(0);
templNorm = templSum2;
}
templSum2 /= invArea;
templNorm = std::sqrt(templNorm);
templNorm /= std::sqrt(invArea); // care of accuracy here
CV_Assert(sqsum.data != NULL);
q0 = (double*)sqsum.data;
q1 = q0 + templ.cols*cn;
q2 = (double*)(sqsum.data + templ.rows*sqsum.step);
q3 = q2 + templ.cols*cn;
}
CV_Assert(sum.data != NULL);
double* p0 = (double*)sum.data;
double* p1 = p0 + templ.cols*cn;
double* p2 = (double*)(sum.data + templ.rows*sum.step);
double* p3 = p2 + 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 = result.ptr<float>(i);
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 += t*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 )
{
double diff2 = MAX(wndSum2 - wndMean2, 0);
if (diff2 <= std::min(0.5, 10 * FLT_EPSILON * wndSum2))
t = 0; // avoid rounding errors
else
t = std::sqrt(diff2)*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;
}
}
}
}
method
是算法标识, cn
我理解为通道数CV_TM_CCOEFF_NORMED
,此处会是 numType=1
和 isNormed=true
method=CV_TM_CCOEFF_NORMED
,我们进到 else
代码段其中 Mean_{I,T} 表示 I 在当前位置和 T 同样大小区域中的均值
至此 R 中的分子已经得到了,之后算分母
$$ diff2 = wndSum2 - wndMean2\ = \sum_{k=0}^{H_T} \sum_{l=0}^{W_T} I_{k+i,l+j}^2-(\sum_{k=0}^{H_T} \sum_{l=0}^{W_T} I_{k+i,l+j})^2 $$