Very fast template matching（非常快的模板匹配）

Integral image
代数矩的快速计算
快速模板匹配
算法步骤
结果

主要讲述文章“Very Fast Template Matching” 中用到的特征和算法的基本思路。（用尽量少的公式）

Integral image

本文作者从2001 年人脸识别神作 “Robust Real-time Object Detection” 中提到的图片特征 integral image中得到启发，发现这个特征可以作为模板匹配的特征，并且通过快速计算图片某个位置对应的代数矩达到快速匹配的目的。假定 I ( x , y ) I(x,y) I(x,y)代表一张图片，则该图片中任意一个点 ( x , y ) (x,y) (x,y)对应的integral image的值是该点所在图片中左方和上方所有像素的和，如图：

Integral image的好处是，当一张图片的integral image已知，这张图片的任何一个区域的像素的总和可以以极快的速度计算出来。

例如，如果我们对图片 I ( x , y ) I(x,y) I(x,y)中的区域 D D D感兴趣，那么在integral image中我们通过计算点1,2，3,4对应的值就能得到图片 I ( x , y ) I(x,y) I(x,y)中的区域 D D D的像素的总和，也就是integral image中的： 4 + 1 − 3 − 2 4+1-3-2 4+1−3−2。

代数矩的快速计算

上一部分简要介绍了integral image的定义和在图片像素计算中的优点，这里以integral image为基础，介绍图片特定的点 ( x , y ) (x,y) (x,y)的代数矩的计算方式。这也是本文作者的创新点。
一个矩形区域 I I I 的代数矩定义为：
m p q = ∑ x = x a x b ∑ y = y a y b x p y q I ( x , y ) (1) m_{pq} = \sum_{x=x_a}^{x_b}\sum_{y = y_a}^{y_b}x^py^qI(x,y)\tag{1} mpq=x=xa∑xby=ya∑ybxpyqI(x,y)(1)
假设 x 0 , y 0 x_0,y_0 x0,y0代表该矩形的中心，定义这个矩形的中心距为：
u p q = ∑ x = x a x b ∑ y = y a y b ( x − x 0 ) p ( y − y 0 ) q I ( x , y ) (2) u_{pq} = \sum_{x=x_a}^{x_b}\sum_{y = y_a}^{y_b}{(x-x_0)}^p{(y-y_0)}^qI(x,y)\tag{2} upq=x=xa∑xby=ya∑yb(x−x0)p(y−y0)qI(x,y)(2)
中心距和代数矩满足展开式：
u 00 = m 00 u_{00} = m_{00} u00=m00
u 10 = m 10 − x 0 m 00 u_{10} = m_{10}-x_0m_{00} u10=m10−x0m00
u 01 = m 01 − y 0 m 00 u_{01} = m_{01}-y_0m_{00} u01=m01−y0m00
u 20 = m 20 − 2 x 0 m 10 + x 0 2 m 00 u_{20} = m_{20}-2x_0m_{10}+x_0^2m_{00} u20=m20−2x0m10+x02m00
u 11 = m 11 − x 0 m 01 − y 0 m 10 + x 0 y 0 m 00 u_{11} = m_{11}-x_0m_{01}-y_0m_{10}+x_0y_0m_{00} u11=m11−x0m01−y0m10+x0y0m00
u 02 = m 02 − 2 y 0 m 01 + y 0 2 m 00 u_{02} = m_{02}-2y_0m_{01}+y_0^2m_{00} u02=m02−2y0m01+y02m00
…
计算代数矩和中心距的原因是，作者希望通过以中心距的一组值作为特征用于模板匹配。计算代数矩的速度类似于计算integral image，在计算时能够对图像矩阵进行直接操作。

快速模板匹配

设模板匹配任务中的模板为 T T T图片为 I I I，像素值均表示为 f ( x , y ) f(x,y) f(x,y)，一般选 p = 2 , q = 2 p=2,q=2 p=2,q=2
通过上述介绍，我们不难得到模板 T T T的中心距 U T = ( u 00 T , . . . , u p q T ) U_T = (u_{00}^T,...,u_{pq}^T) UT=(u00T,...,upqT)，定义
U I ( u , v ) = ( u 00 I , . . . , u p q I ) U_I(u,v) = (u_{00}^I,...,u_{pq}^I) UI(u,v)=(u00I,...,upqI)为中心在 ( u , v ) (u,v) (u,v)的矩形的中心距。通过计算模板与图片中每一个点对应的中心距的“距离”，来判定模板匹配的位置。例如
h ( u , v ) = ∑ p , q u p q I ( u , v ) u p q t ( ∑ p , q ( u p q I ( u , v ) ) 2 ) (3) h(u,v) = \frac{\sum_{p,q} u_{pq}^I(u,v)u_{pq}^t}{\sqrt{(\sum_{p,q}(u_{pq}^I(u,v))^2)}}\tag{3} h(u,v)=(∑p,q(upqI(u,v))2) ∑p,qupqI(u,v)upqt(3)
h ( u , v ) h(u,v) h(u,v)的值最大的点对应最准确的匹配点。

但是这样计算出的 h ( u , v ) h(u,v) h(u,v)对于噪声比较敏感，所以文章提供了另外一个思路:
在以归一化互相关为特征的模板匹配计算方式为：
h ( u , v ) = ∑ x , y I ( u + x , v + y ) T ( x , y ) ( ∑ x , y I 2 ( u + x , v + y ) (4) h(u,v) = \frac{\sum_{x,y} I(u+x,v+y)T(x,y)}{\sqrt{(\sum_{x,y} I^2(u+x,v+y)}}\tag{4} h(u,v)=(∑x,yI2(u+x,v+y) ∑x,yI(u+x,v+y)T(x,y)(4)
通过多项式 p p p对像素 I I I进行近似，将上式转变为：
h ( u , v ) = U I ′ ( u , v ) B − 1 U t U I ′ ( u , v ) B − 1 U I ( u , v ) (6) h(u,v) = \frac{U_{I}'(u,v)B^{-1}U_{t}}{\sqrt{U_{I}'(u,v)B^{-1}U_{I}(u,v)}}\tag{6} h(u,v)=UI′(u,v)B−1UI(u,v) UI′(u,v)B−1Ut(6)
推导请移步论文

算法步骤

1.计算出模板的中心距
2. 计算出图片的每一阶integral image
3. 计算出图片中每一个可能成为矩形中心的点中心距
4. 根据计算出的图片和模板的中心距计算特征
5. 得到特征最大的点，也就是模板匹配的结果。

结果

模板和图片如下：

我的matlab仿真结果：

代码移步：Matlab代码