OpenCV 模板匹配 matchTemplate 源码解析

2022-11-24 17:26:34 浏览数 (1)

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

简介

  • OpenCV 实现的 matchTemplate 速度很快，核心提速在于使用了卷积加速和累加和技巧
  • 参考 OpenCV 版本 4.4.0，源码位于: opencvsourcesmodulesimgprocsrctemplmatch.cpp
  • 官方文档：https://docs.opencv.org/2.4.13.7/modules/imgproc/doc/object_detection.html?highlight=matchtemplate#cv2.matchTemplate

例程选取

  • 之前我们记录过模板匹配函数用法，损失函数分为 差值平方和，相关度，去均值相关度 三种，并且每种损失可以选择是否归一化
  • 其中归一化的去均值相关度 （对应 CV_TM_CCOEFF_NORMED）运算过程最为复杂，此处以该损失函数为例解读源码，其余函数可以以此类推
  • 为了简单且不失一般性，我们假设输入的图是单通道的数据，多通道以此类推
  • 用 I 表示待匹配图像（大图），T 表示模板图像（小图），w,h 表示模板宽高，计算公式：

源码解析

生成内积图
  • 几种损失函数最核心的计算都离不开模板在原图中的卷积运算，因此所有模板匹配都预先计算好了卷积图
  • 这部分运算在matchTemplate 函数中实现，源码
代码语言：javascript复制
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 表示模板图像的宽度
  • 源码：
代码语言：javascript复制
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=1 和 isNormed=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 更新：

其中 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 的位置
循环得到匹配结果矩阵
  • 对每个可选的位置重复上述流程，计算得到最终的匹配结果矩阵
  • 将结果返回给用户

参考资料

  • https://cloud.tencent.com/developer/article/2068694
  • https://blog.csdn.net/github_38148039/article/details/109469238
  • https://www.yumpu.com/en/document/view/47330650/intelr-ipp-quick-reference
c++ opencv https 网络安全

0 人点赞