模板匹配加速——opencv

介绍

模板匹配是一个图像处理问题，当对象的姿势（X、Y、+）未知时，它使用模板图像在另一个搜索图像中查找其位置。在这篇文章中，我们实现一个算法，该算法使用对象的边缘信息来识别搜索图像中的对象。

背景

由于模板匹配的速度和可靠性问题，模板匹配本质上是一个棘手的问题。当对象部分可见或与其他对象混合时，该解决方案应针对亮度变化保持稳健，最重要的是，该算法的计算效率应高。解决这个问题主要有两种方法，基于灰值的匹配（或基于区域的匹配）和基于特征的匹配（非基于区域的匹配）。

基于灰值的方法：在基于灰值的匹配中，规范化交叉关联（NCC）算法是从过去开始认识的。这通常通过减去均值和除以标准差来在每个步骤中完成。模板 t（x， y） 与子图像 f（x， y） 的交叉相关性是：

其中 n 是 t（x、y） 和 f（x， y） 中的像素数。•维基]

尽管此方法针对线性照明变化是健壮的，但当对象部分可见或对象与其他对象混合时，算法将失败。此外，该算法的计算成本很高，因为它需要计算模板图像中所有像素与搜索图像之间的相关性。

基于特征的方法：在图像处理领域采用几种基于特征的模板匹配方法。与基于边的对象识别一样，对象边缘是用于匹配的要素，在通用霍夫变换中，对象几何特征将用于匹配。

在这篇文章中，我们实现一个算法，该算法使用对象的边缘信息来识别搜索图像中的对象。此实现使用开源计算机视觉库作为平台。

编译示例代码

我们使用 OpenCV 2.0 和 Visual studio 2008 来开发此代码。要编译示例代码，我们需要安装 OpenCV。

OpenCV可以从这里免费下载。OpenCV（开源Computer Vision）是一个用于实时计算机视觉的编程函数库。下载 OpenCV 并将其安装到您的系统中。安装信息可从这里读取。

我们需要配置我们的可视化工作室环境。此信息可在此处阅读。

算法

在这里，我们解释一种基于边缘的模板匹配技术。边缘可以定义为数字图像中的点，其中图像亮度会急剧变化或具有不连续性。从技术上讲，它是一种离散的分化操作，计算图像强度函数梯度的近似值。

边缘检测的方法有很多，但大多数方法可以分为两类：基于搜索和基于零交叉。基于搜索的方法首先计算边强度的度量，通常是一阶导数表达式（如梯度幅度）来检测边缘，然后使用对边缘的局部方向（通常是梯度方向）的计算估计来搜索梯度幅度的局部方向最大值。在这里，我们使用这样的方法实现由索贝尔称为索贝尔运算符。操作员计算每个点的图像强度的渐变，给出从浅到暗的最大可能增加方向以及该方向的变化速率。

我们使用这些梯度或导数在X方向和Y方向进行匹配。

此算法涉及两个步骤。首先，我们需要创建模板图像的基于边缘的模型，然后使用此模型在搜索图像中搜索。

创建基于边的模板模型

我们首先从模板图像的边缘创建一个数据集或模板模型，用于在搜索图像中查找该对象的姿势。

在这里，我们使用 Canny 边缘检测方法的变体来查找边缘。你可以在这里阅读更多关于坎尼的边缘检测。对于边缘提取，Canny 使用以下步骤：

第 1 步：查找图像的强度渐变

在模板图像上使用 Sobel 筛选器，该筛选器返回 X （Gx） 和 Y （Gy） 方向的渐变。在此梯度中，我们将使用以下公式计算边缘幅度和方向：

我们使用 OpenCV 函数来查找这些值。

cvSobel( src, gx, 1,0, 3 ); //gradient in X direction
cvSobel( src, gy, 0, 1, 3 ); //gradient in Y direction

for( i = 1; i < Ssize.height-1; i++ )
{
    for( j = 1; j < Ssize.width-1; j++ )
    {          
        _sdx = (short*)(gx->data.ptr + gx->step*i);
        _sdy = (short*)(gy->data.ptr + gy->step*i);
        fdx = _sdx[j]; fdy = _sdy[j];
        // read x, y derivatives

        //Magnitude = Sqrt(gx^2 +gy^2)
        MagG = sqrt((float)(fdx*fdx) + (float)(fdy*fdy));
        //Direction = invtan (Gy / Gx)
        direction =cvFastArctan((float)fdy,(float)fdx);
        magMat[i][j] = MagG;
            
        if(MagG>MaxGradient)
            MaxGradient=MagG;
            // get maximum gradient value for normalizing.

            
        // get closest angle from 0, 45, 90, 135 set
        if ( (direction>0 && direction < 22.5) || 
              (direction >157.5 && direction < 202.5) || 
              (direction>337.5 && direction<360)  )
            direction = 0;
        else if ( (direction>22.5 && direction < 67.5) || 
                  (direction >202.5 && direction <247.5)  )
            direction = 45;
        else if ( (direction >67.5 && direction < 112.5)||
                  (direction>247.5 && direction<292.5) )
            direction = 90;
        else if ( (direction >112.5 && direction < 157.5)||
                  (direction>292.5 && direction<337.5) )
            direction = 135;
        else
            direction = 0;
            
        orients[count] = (int)direction;
        count++;
    }
}

找到边缘方向后，下一步是关联图像中可跟踪的边缘方向。描述周围像素有四种可能的方向：0 度、45 度、90 度和 135 度。我们指定所有方向到这些角度。

第 2 步：应用非最大抑制

找到边缘方向后，我们将进行非最大抑制算法。非最大抑制沿边缘方向跟踪左右像素，如果当前像素幅度小于左右像素幅度，则禁止抑制。这将导致图像变薄。

for( i = 1; i < Ssize.height-1; i++ )
{
    for( j = 1; j < Ssize.width-1; j++ )
    {
        switch ( orients[count] )
        {
            case 0:
                leftPixel  = magMat[i][j-1];
                rightPixel = magMat[i][j+1];
                break;
            case 45:
                leftPixel  = magMat[i-1][j+1];
                rightPixel = magMat[i+1][j-1];
                break;
            case 90:
                leftPixel  = magMat[i-1][j];
                rightPixel = magMat[i+1][j];
                break;
            case 135:
                leftPixel  = magMat[i-1][j-1];
                rightPixel = magMat[i+1][j+1];
                break;
        }
        // compare current pixels value with adjacent pixels
        if (( magMat[i][j] < leftPixel ) || (magMat[i][j] < rightPixel ) )
            (nmsEdges->data.ptr + nmsEdges->step*i)[j]=0;
        Else
            (nmsEdges->data.ptr + nmsEdges->step*i)[j]=
                                 (uchar)(magMat[i][j]/MaxGradient*255);
        count++;
    }
}

第 3 步：执行滞后阈值

使用滞后设置阈值需要两个阈值：高和低。我们应用一个高阈值来标出那些边缘， 我们可以相当肯定是真实的。从这些开始，使用之前派生的方向信息，其他边缘可以通过图像进行跟踪。在跟踪边时，我们应用较低的阈值，允许我们跟踪边缘的微弱部分，只要我们找到一个起点。

_sdx = (short*)(gx->data.ptr + gx->step*i);
_sdy = (short*)(gy->data.ptr + gy->step*i);
fdx = _sdx[j]; fdy = _sdy[j];
    
MagG = sqrt(fdx*fdx + fdy*fdy); //Magnitude = Sqrt(gx^2 +gy^2)
DirG =cvFastArctan((float)fdy,(float)fdx);     //Direction = tan(y/x)

////((uchar*)(imgGDir->imageData + imgGDir->widthStep*i))[j]= MagG;
flag=1;
if(((double)((nmsEdges->data.ptr + nmsEdges->step*i))[j]) < maxContrast)
{
    if(((double)((nmsEdges->data.ptr + nmsEdges->step*i))[j])< minContrast)
    {
        (nmsEdges->data.ptr + nmsEdges->step*i)[j]=0;
        flag=0; // remove from edge
        ////((uchar*)(imgGDir->imageData + imgGDir->widthStep*i))[j]=0;
    }
    else
    {   // if any of 8 neighboring pixel is not greater than max contraxt remove from edge
        if( (((double)((nmsEdges->data.ptr + nmsEdges->step*(i-1)))[j-1]) < maxContrast) &&
            (((double)((nmsEdges->data.ptr + nmsEdges->step*(i-1)))[j]) < maxContrast)   &&
            (((double)((nmsEdges->data.ptr + nmsEdges->step*(i-1)))[j+1]) < maxContrast) &&
            (((double)((nmsEdges->data.ptr + nmsEdges->step*i))[j-1]) < maxContrast)     &&
            (((double)((nmsEdges->data.ptr + nmsEdges->step*i))[j+1]) < maxContrast)     &&
            (((double)((nmsEdges->data.ptr + nmsEdges->step*(i+1)))[j-1]) < maxContrast) &&
            (((double)((nmsEdges->data.ptr + nmsEdges->step*(i+1)))[j]) < maxContrast)   &&
            (((double)((nmsEdges->data.ptr + nmsEdges->step*(i+1)))[j+1]) < maxContrast))
        {
            (nmsEdges->data.ptr + nmsEdges->step*i)[j]=0;
            flag=0;
            ////((uchar*)(imgGDir->imageData + imgGDir->widthStep*i))[j]=0;
        }
    }
}

第 4 步：保存数据集

提取边后，我们将所选边的 X 和 Y 导数以及坐标信息作为模板模型。这些坐标将被重新排列，以反映起点作为重心。

查找基于边的模板模型

算法中的下一个任务是使用模板模型在搜索图像中查找对象。我们可以看到我们从包含一组点的模板图像创建的模型：，及其在 X 和 Y 方向的渐变，其中 i = 1 ...n，n是模板 （T） 数据集中的元素数。

我们还可以在搜索图像 （S） 中找到渐变，其中 u = 1...搜索图像中的列数。

在匹配过程中，应使用相似性度量度将模板模型与所有位置的搜索图像进行比较。相似性度量背后的理念是采取模板图像的梯度矢量的所有规范化点乘量的总和，并在模型数据集的所有点上搜索图像。这将导致搜索图像中每个点的分数。这可表述如下：

如果模板模型和搜索图像之间完全匹配，则此函数将返回分数 1。分数对应于搜索图像中可见的对象部分。如果搜索图像中不存在对象，则分数将为 0。

cvSobel( src, Sdx, 1, 0, 3 );  // find X derivatives
cvSobel( src, Sdy, 0, 1, 3 ); // find Y derivatives
for( i = 0; i < Ssize.height; i++ )
{
    for( j = 0; j < Ssize.width; j++ )
    { 
        partialSum = 0; // initilize partialSum measure
        for(m=0;m<noOfCordinates;m++)
        {
            curX    = i + cordinates[m].x ;    // template X coordinate
            curY    = j + cordinates[m].y ; // template Y coordinate
            iTx    = edgeDerivativeX[m];    // template X derivative
            iTy    = edgeDerivativeY[m];    // template Y derivative

            if(curX<0 ||curY<0||curX>Ssize.height-1 ||curY>Ssize.width-1)
                continue;
                     
            _Sdx = (short*)(Sdx->data.ptr + Sdx->step*(curX));
            _Sdy = (short*)(Sdy->data.ptr + Sdy->step*(curX));
                        
            iSx=_Sdx[curY]; // get curresponding  X derivative from source image
            iSy=_Sdy[curY];// get curresponding  Y derivative from source image
                        
            if((iSx!=0 || iSy!=0) && (iTx!=0 || iTy!=0))
            {
                //partial Sum  = Sum of(((Source X derivative* Template X drivative)
                //+ Source Y derivative * Template Y derivative)) / Edge
                //magnitude of(Template)* edge magnitude of(Source))
                partialSum = partialSum + ((iSx*iTx)+(iSy*iTy))*
                            (edgeMagnitude[m] * matGradMag[curX][curY]);
                                    
            }

在实际情况下，我们需要加快搜索过程。这可以通过各种方法实现。第一个方法是使用平均的属性。找到相似性度量时，如果我们可以为相似性度量值设置最低分数 （S^min），则无需评估模板模型中的所有点。要检查特定点 J 的部分分数 S_u，v，我们必须找到点 m. Sm 的部分总和，可以定义如下：

显然，该金额的剩余条款较小或等于 1。因此，如果 ，我们可以停止评估。

另一个标准可能是，任何时间点的部分分数应大于最低分数。即.使用此条件时，匹配速度将非常快。但问题是，如果首先检查对象的缺失部分，则部分总和将很低。在这种情况下，该对象的实例将不被视为匹配项。我们可以用另一个条件修改这一点，其中我们检查模板模型的第一部分与安全停止标准，其余与硬条件 ， 。用户可以指定贪婪参数 （g），其中使用硬条件检查模板模型的分数。因此，如果 g=1，模板模型中的所有点都使用硬条件进行检查，如果 g=0，则所有点将仅使用安全条件进行检查。我们可以制定如下程序。

部分分数的评估可以停止在以下：

// stoping criterias to search for model
double normMinScore = minScore /noOfCordinates; // precompute minumum score 
double normGreediness = ((1- greediness * minScore)/(1-greediness)) /noOfCordinates;
// precompute greedniness

sumOfCoords = m + 1;
partialScore = partialSum /sumOfCoords ;
// check termination criteria
// if partial score score is less than the score than
// needed to make the required score at that position
// break serching at that coordinate.
if( partialScore < (MIN((minScore -1) + 
        normGreediness*sumOfCoords,normMinScore*  sumOfCoords)))
    break;

这种相似性度量具有以下几个优点：相似性度量与非线性照明变化是不变的，因为所有梯度矢量都已规范化。由于边缘滤波上没有分段，因此将显示对照明任意变化的真不变性。更重要的是，当对象部分可见或与其他对象混合时，这种相似性度量是健壮的。

增强

此算法可能提供各种增强功能。为了进一步加快搜索过程，可以使用金字塔式方法。在这种情况下，搜索以小图像大小的低分辨率开始。这对应于金字塔的顶部。如果搜索在此阶段成功，则搜索将继续在金字塔的下一个级别，该级别表示更高分辨率的图像。以这种方式，搜索继续，因此，结果被优化，直到原始图像大小，即，到达金字塔的底部。

另一个增强是可能的，通过扩展旋转和缩放算法。这可以通过创建用于旋转和缩放的模板模型以及使用所有这些模板模型执行搜索来完成。

引用

1. 机器视觉算法和应用 [卡斯滕·斯特格、马库斯·乌尔里希、克里斯蒂安·维德曼]
2. 数字图像处理 [拉斐尔·冈萨雷斯，理查德·尤金·伍兹]
3. http://dasl.mem.drexel.edu/alumni/bGreen/www.pages.drexel.edu/_weg22/can_tut.html
4. https://www.codeproject.com/articles/99457/edge-based-template-matching