文章目录
  1. 1. 模板匹配简介
  2. 2. NCC算法实现
    1. 2.0.0.1. 快速傅里叶变换计算卷积
    2. 2.0.0.2. 图像金字塔
    3. 2.0.0.3. 旋转匹配
    4. 2.0.0.4. End

模板匹配是机器视觉算法中非常重要的一项技术,文本介绍的是一种基于灰度值的匹配方法,参考了《机器视觉算法与应用》相关章节。重点讲解了其中匹配分数计算的原理及代码实现。

模板匹配简介

模板匹配(Template Matching)是在目标图像中快速定位模板图像的一种技术,这种技术非常广泛的应用在实际的机器视觉应用当中。比如在各种的质量检测或者识别应用中定位参考目标的位置,如下面左图中在快递单中定位logo,右图中清点工件的个数:




模板匹配实际应用

本文将介绍一种最简单的模板匹配技术,基于NCC模板匹配。
不凡假设模板的大小是mxm, 目标图像的大小是nxn,n>m。我们这里先不考虑旋转和缩放的情况,一种最简单的想法就是在nxn的目标图像中用mxm的窗口遍历图像,比较每个位置的窗口图像和模板图像的差异。自然而然可以想到计算每个像素差的平方和来表示目标和模板的差异:

$$d_{f,t}^2(u,v)=\sum_{x,y}[f(x,y)-t(x-u,y-u)]^2$$

将上述式子展开:

$$d_{f,t}^2(u,v)=\sum_{x,y}(f^2(x,y) - 2f(x,y)t(x-u,y-u) + t^2(x-u,y-v))$$

其中 $\sum t^2(x-u,y-v)$ 是常量。如果式子$\sum f^2(x,y)$近似常量的话,剩下的cross correlation:
$$c(u,v) = \sum_{x,y}f(x,y)t(x-u,y-u)$$
就可以作为目标和模板的匹配度量。

实际上$\sum f^2(x,y)$不会是常量,如果仅使用cross correlation作为度量会有问题。比如全白的图像的匹配值会比实际匹配目标的匹配值高,而且$c(u,v)$的值会受模板大小的影响,同时会对光线非常敏感。因此,需要对cross correlation进行归一化:
$$\gamma(u,v) = \frac{\sum_{x,y}[f(x,y)-\overline{f}_{u,v}][t(x-u,y-v)-\overline{t}]}{\sqrt{\sum_{x,y}[f(x,y)-\overline{f}_{u,v}]^2[t(x-u,y-v)-\overline{t}]^2}}$$
其中${ \overline{t} }$是模板图像的平均值,${ \overline{f}_{u,v} }$是搜索窗口图像的平均值。得到的这个匹配度量称为归一化相关系数(Normalized cross correlation),它能够有效的避免上述问题,是一个非常实用的匹配度量。同时我们注意NCC度量的范围在[-1,1]。

NCC算法实现

如果直接实现NCC算法,无疑会产生巨大的运算量。因此下面将会提到两种方法来提高匹配的性能。

快速傅里叶变换计算卷积

快速傅里叶变换在图像处理中的应用(目前我了解的)有两个:用来过滤图像,如过滤桌面的纹理,使得目标能够更好的被分割。另一个就是用来快速计算卷积。(关于傅里叶变换的相关实现以后细讲)

上述NCC度量公式可以展开后可以得到:
记$f’(x,y) \equiv f(x,y)-\overline{f}_{u,v}$

$$ 分子 = \sum_{x,y}f’(x,y)t(x-u,y-v)$$

$$ 分母 = \sqrt{\sum_{x,y}[t(x-u,y-v)-\overline{t}]^2[f^2(x,y)-\overline{f}_{u,v}^2]}$$
(推倒过程并不复杂,略去)

分子注意到是一个卷积的形式,可以用快速傅里叶变换来进行计算。分母中$\sum_{x,y}[t(x-u,y-v)-\overline{t}]^2 (1)$和$\overline{f}_{u,v}$可以事先计算好。对于$f^2(x,y)$可以用积分图来进行快速计算。下面给出OpenCV关于NCC的实现,下面关键地方给出注释:

//快速傅里叶变换计算卷积,结果存放在result中
crossCorr( img, templ, result, result.size(), result.type(), Point(0,0), 0, 0);

if( method == CV_TM_CCORR )
    return;

double invArea = 1./((double)templ.rows * templ.cols);
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//以下是NCC实现,支持彩色图像
{
    integral(img, sum, sqsum, CV_64F);//计算积分图像,sum为f(x,y)积分图,sqsum为f^2(x,y)积分图
    meanStdDev( templ, templMean, templSdv );

    //模板的方差,用来计算分母中的(1)
    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

    q0 = (double*)sqsum.data;
    q1 = q0 + templ.cols*cn;  // The first row and column of the integral image are all zeros
    q2 = (double*)(sqsum.data + templ.rows*sqsum.step);
    q3 = q2 + templ.cols*cn;  
}

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 )
        {
            t = std::sqrt(MAX(wndSum2 - wndMean2,0))*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;
    }
}
图像金字塔

即使采用了上述的方法优化性能,所需的计算量也是极大的,因为要在所有的图像位置全部计算出匹配度量。为了节约开销,可以构建图像金字塔。




金字塔图像


最底层(第0层)为原始图像,每往上一层,图像宽度高度缩小一倍,图像大小缩小4倍。
我们根据模板大小选择金字塔的层数,一般构造2~3层金字塔。
在金字塔的顶层,对所有的图像位置计算匹配度量。我们选择合适的阈值过滤匹配度量,对过滤后的图像选取局部最大值作为潜在的匹配点。然后根据潜在匹配点的位置,逐层向下搜索,在下面一层我们只需要在潜在匹配点一定范围进行搜索,一般来说10x10的范围已经十分足够了。
因此我们可以看到,假设金字塔构建了3层,那么大约计算量缩小了4x4x4=64倍。这是非常巨大的性能优化,可以说如果没有这一步的优化,NCC难以在实际中应用。当然这里的图像金字塔的方法也可以被其他方法替代。
同时,这个图像金字塔在其他的匹配中也可以应用到,比如基于边缘的模板匹配,基于彩色信息的模板匹配等等。这是一个非常使用的优化技术。

旋转匹配

事实上,以上的介绍的模板匹配并没有考虑旋转目标的匹配。这实现起来会更加复杂点,感兴趣的读者可以尝试着自己实现,需要注意的是如果对模板进行360全方位进行旋转的话要考虑旋转后的模板会有“黑边”(这个就要考虑带掩膜的模板了)。同时和只考虑平移是类似的,只需要在金字塔最顶层进行全方位的比对(此时角度搜索的步长可以不用十分精细),然后对潜在匹配点在小角度范围搜索。

End

NCC是基于图像灰度的一种匹配技术,在大部分情况下都可以适用。但是一定程度下还是会受光照不均匀影响,另外对于目标物体存在混乱,遮挡等情况匹配效果十分不理想。下面有机会的话会介绍基于边缘的模板匹配。

PS:本文适合有一定图像处理基础的人阅读,由于公式编辑实在麻烦,省去部分推倒过程。

文章目录
  1. 1. 模板匹配简介
  2. 2. NCC算法实现
    1. 2.0.0.1. 快速傅里叶变换计算卷积
    2. 2.0.0.2. 图像金字塔
    3. 2.0.0.3. 旋转匹配
    4. 2.0.0.4. End