OpenCV 模板匹配 matchTemplate 源码解析

本文最后更新于:2022年11月9日 下午

OpenCV 实现了图像平移模板匹配的功能,封装在函数接口 matchTemplate 中,本文解析该功能的实现源码。

简介

例程选取

  • 之前我们记录过模板匹配函数用法,损失函数分为 差值平方和相关度去均值相关度 三种,并且每种损失可以选择是否归一化
  • 其中归一化的去均值相关度 (对应 CV_TM_CCOEFF_NORMED)运算过程最为复杂,此处以该损失函数为例解读源码,其余函数可以以此类推
  • 为了简单且不失一般性,我们假设输入的图是单通道的数据,多通道以此类推
  • 用 $I$ 表示待匹配图像(大图),$T$ 表示模板图像(小图),$w,h$ 表示模板宽高,计算公式:
$$ R(x, y)=\frac{\sum_{x^{\prime}, y^{\prime}}\left(T^{\prime}\left(x^{\prime}, y^{\prime}\right) \cdot I^{\prime}\left(x+x^{\prime}, y+y^{\prime}\right)\right)}{\sqrt{\sum_{x^{\prime}, y^{\prime}} T^{\prime}\left(x^{\prime}, y^{\prime}\right)^{2} \cdot \sum_{x^{\prime}, y^{\prime}} I^{\prime}\left(x+x^{\prime}, y+y^{\prime}\right)^{2}}} $$
  • 其中:
$$ \begin{array}{l} T^{\prime}\left(x^{\prime}, y^{\prime}\right)=T\left(x^{\prime}, y^{\prime}\right)-1 /(w \cdot h) \cdot \sum_{x^{\prime \prime}, y^{\prime \prime}} T\left(x^{\prime \prime}, y^{\prime \prime}\right) \\ I^{\prime}\left(x+x^{\prime}, y+y^{\prime}\right)=I\left(x+x^{\prime}, y+y^{\prime}\right)-1 /(w \cdot h) \cdot \sum_{x^{\prime \prime}, y^{\prime \prime}} I\left(x+x^{\prime \prime}, y+y^{\prime \prime}\right) \end{array} $$

源码解析

生成内积图

  • 几种损失函数最核心的计算都离不开模板在原图中的卷积运算,因此所有模板匹配都预先计算好了卷积图
  • 这部分运算在matchTemplate 函数中实现,源码
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
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 变量中

计算 CCOEFF_NORMED 损失

  • 不考虑 mask 的情况下,OpenCV 模板匹配核心用的是 common_matchTemplate 函数
  • 我们定义待匹配的单通道图像(大图)为 $I$,模板单通道图像(小图)为 $T$,宽度$W$,高度$H$,均值 $Mean$,标准差 $Std$
  • 变量会带下标,例如: $W_T$ 表示模板图像的宽度
  • 源码:
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
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;
}
}
}
}
1 行
  • 函数定义
  • img 表示我们定义的 $I$,templ 表示我们定义的 $T$,result 表示已经算好的卷积图
  • method 是算法标识, cn 我理解为通道数
6-10 行
  • 确定损失函数和是否归一化,假设我们选择的算法是 CV_TM_CCOEFF_NORMED ,此处会是 numType=1isNormed=true
12 行
  • $invArea$ 为 $\frac{1}{W_TH_T}$
14-17 行
  • 定义 $I$ 的部分和矩阵 $sum$ 和部分平方和矩阵 $sqsum$

    含义为 $sum_{p,q}$ 中元素值为 $\sum_{i=0}{p-1}\sum_{j=0}{q-1}I_{i,j}$ ,这样可以以 $O(1)$ 复杂度计算 $I$ 中任意矩形包围的元素的和

  • 定义 $T$ 的均值 $templMean$ ,标准差 $templSdv$

  • 定义 $sqsum$ 的矩形四角指针 $q0,q1,q2,q3$

  • 定义临时变量 $templNorm$ , $templSum2$

26 行
  • 因为 method=CV_TM_CCOEFF_NORMED ,我们进到 else 代码段
  • 该行计算 $sum$ 和 $sqsum$
27 行
  • 该行计算 $templMean = Mean_T$,$ templSdv = Std_T$
29-43 行
  • $templNorm = Std_T^2$
  • $templSum2 = templNorm + Mean_T^2 = Std_T^2 + Mean_T^2$
45-47 行
  • $templSum2 = (Std_T^2 + Mean_T^2)W_TH_T$
  • $templNorm = Std_T \sqrt{W_TH_T}$
50-53 行
  • $q0,q1,q2,q3$ 分别指向 $sqsum$ 中 $T$ 同样大小的最左上角矩形四个点,分别是左上角、右上角、左下角、右下角
57-60 行
  • $p0,p1,p2,p3$ 分别指向 $sum$ 中 $T$ 同样大小的最左上角矩形四个点,分别是左上角、右上角、左下角、右下角
68-74 行
  • 整理图像的行列下标,对于理解上不用理,记住 $i$ 表示行,$j$ 表示列就可以了,初始值均为 $0$
75-76 行
  • 定义 $num$ 为 当前 $i,j$ 的 $ T$ 在 $I$ 上卷积的结果,即

$$
num = \sum_{k=0}^{H_T} \sum_{l=0}^{W_T}T_{k,l}I_{k+i,l+j}
$$

  • $t$ 为临时变量
  • 定义窗口临时变量 $wndMean2 = 0$, $wndSum2 = 0 $
80-85 行
  • 由于 $numType = 1$,因此进入了 $80$ 行的 if 代码段
  • $t$ 为 $I$ 在当前位置($i,j$)下和 $T$ 相同大小的矩形范围内的元素和,即:

$$
t =\sum_{k=0}^{H_T} \sum_{l=0}^{W_T} I_{k+i,l+j}
$$

  • $wndMean2 = t^2$

$$
wndMean2=(\sum_{k=0}^{H_T} \sum_{l=0}^{W_T} I_{k+i,l+j})^2
$$

  • $num$ 更新:
$$ num =\sum_{k=0}^{H_T} \sum_{l=0}^{W_T}T_{k,l}I_{k+i,l+j}-\sum_{k=0}^{H_T} \sum_{l=0}^{W_T} I_{k+i,l+j}Mean_T $$
  • 这个 $num$ 表示的就是原始损失函数 $R$ 中的分子,我们此处使用官方文档的表示方法:
$$ R(x, y)=\frac{\sum_{x^{\prime}, y^{\prime}}\left(T^{\prime}\left(x^{\prime}, y^{\prime}\right) \cdot I^{\prime}\left(x+x^{\prime}, y+y^{\prime}\right)\right)}{\sqrt{\sum_{x^{\prime}, y^{\prime}} T^{\prime}\left(x^{\prime}, y^{\prime}\right)^{2} \cdot \sum_{x^{\prime}, y^{\prime}} I^{\prime}\left(x+x^{\prime}, y+y^{\prime}\right)^{2}}} $$
  • 分子部分:
$$ \begin{array}{l} &\sum_{x^{\prime}, y^{\prime}}\left(T^{\prime}\left(x^{\prime}, y^{\prime}\right) \cdot I^{\prime}\left(x+x^{\prime}, y+y^{\prime}\right)\right)\\ =&\sum_{x^{\prime}, y^{\prime}}(T\left(x^{\prime}, y^{\prime}\right)-1 /(w \cdot h) \cdot \sum_{x^{\prime \prime}, y^{\prime \prime}} T\left(x^{\prime \prime}, y^{\prime \prime}\right)) (I\left(x+x^{\prime}, y+y^{\prime}\right)-1 /(w \cdot h) \cdot \sum_{x^{\prime \prime}, y^{\prime \prime}} I\left(x+x^{\prime \prime}, y+y^{\prime \prime}\right) )\\ =&\sum_{x^{\prime}, y^{\prime}}(T\left(x^{\prime}, y^{\prime}\right)-Mean_T)(I\left(x+x^{\prime}, y+y^{\prime}\right)-Mean_{I,T})\\ =&\sum_{x^{\prime}, y^{\prime}}T\left(x^{\prime}, y^{\prime}\right)I\left(x+x^{\prime}, y+y^{\prime}\right)-\sum_{x^{\prime}, y^{\prime}}T\left(x^{\prime}, y^{\prime}\right)Mean_{I,T}-\sum_{x^{\prime}, y^{\prime}}I\left(x+x^{\prime}, y+y^{\prime}\right)Mean_T+W_TH_TMean_TMean_{I,T}\\ =&\sum_{x^{\prime}, y^{\prime}}T\left(x^{\prime}, y^{\prime}\right)I\left(x+x^{\prime}, y+y^{\prime}\right)-W_TH_TMean_TMean_{I,T}-W_TH_TMean_TMean_{I,T}+W_TH_TMean_TMean_{I,T}\\ =&\sum_{x^{\prime}, y^{\prime}}T\left(x^{\prime}, y^{\prime}\right)I\left(x+x^{\prime}, y+y^{\prime}\right)-W_TH_TMean_TMean_{I,T}\\ =&\sum_{x^{\prime}, y^{\prime}}T\left(x^{\prime}, y^{\prime}\right)I\left(x+x^{\prime}, y+y^{\prime}\right)-\sum_{x^{\prime}, y^{\prime}}I\left(x+x^{\prime}, y+y^{\prime}\right)Mean_T \end{array} $$
  • 其中 $Mean_{I,T}$ 表示 $I$ 在当前位置和 $T$ 同样大小区域中的均值
  • 而这个形式和 $sum$ 表示的含义是完全相同的

至此 $R$ 中的分子已经得到了,之后算分母

87 行

$$
wndMean2 = \frac {t^2}{W_TH_T}= \frac {(\sum_{k=0}^{H_T} \sum_{l=0}^{W_T} I_{k+i,l+j})^2}{W_TH_T}
$$

94-95 行
  • $t$ 为 $I$ 在当前位置($i,j$)下和 $T$ 相同大小的矩形范围内的元素平方和,即:

$$
t =\sum_{k=0}^{H_T} \sum_{l=0}^{W_T} I_{k+i,l+j}^2
$$

  • $wndSum2 = t = \sum_{k=0}^{H_T} \sum_{l=0}^{W_T} I_{k+i,l+j}^2$
107-111 行
$$ 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 $$
  • $diff2$ 表示归一化项中 $\sum_{x^{\prime}, y^{\prime}} I^{\prime}\left(x+x^{\prime}, y+y^{\prime}\right)^{2}$ 的结果,其计算过程与方差类似,结果为数据的平方的均值减均值的平方乘以数据数量。
  • 之前计算得到 $templNorm = Std_T \sqrt{W_TH_T}$,这其实表示归一化项中的 $ \sqrt{\sum_{x^{\prime}, y^{\prime}} T^{\prime}\left(x^{\prime}, y^{\prime}\right)^{2}}$,因为标准差的平方是方差,方差乘以数据数量为 $\sum_{x^{\prime}, y^{\prime}} T^{\prime}\left(x^{\prime}, y^{\prime}\right)^{2}$
  • 因此分母 $t = \sqrt{diff2} \cdot templNorm $,就表示 $R$ 式计算中的分母
114 行
  • 之前的 $num$ 表示了 $R$ 中的分子,分吗母为 $t$,那么除一下就好了,得到了目标 $R$ 式计算的结果
121 行
  • 将按照 $R$ 计算的结果填入之前 $num$ 的位置
循环得到匹配结果矩阵
  • 对每个可选的位置重复上述流程,计算得到最终的匹配结果矩阵
  • 将结果返回给用户

参考资料


OpenCV 模板匹配 matchTemplate 源码解析
https://www.zywvvd.com/notes/study/image-processing/opencv/opencv-matchTemplate/opencv-matchtemplate-src/opencv-matchtemplate-src/
作者
Yiwei Zhang
发布于
2022年11月9日
许可协议