NCC算法实现及其优化

本文将集中探讨一种实现相对简单，效果较好的模板匹配算法（NCC）

\[R(x,y)= \frac{ \sum_{x',y'} (T'(x',y') \cdot I'(x+x',y+y')) }{ \sqrt{\sum_{x',y'}T'(x',y')^2 \cdot \sum_{x',y'} I'(x+x',y+y')^2} } \\ \\ \begin{array}{l} T'(x',y')=T(x',y') - 1/(w \cdot h) \cdot \sum _{ x'',y''} T(x'',y'') \\ I'(x+x',y+y')=I(x+x',y+y') - 1/(w \cdot h) \cdot \sum _{x'',y''} I(x+x'',y+y'') \end{array} \]

具体原理是高中知识，也就是Spearma相关系数，这里就不赘述了。

基本实现代码

float MatchTemplateCPU(int8_t* pSearchImage,
    int8_t* pTemplate,
    int nSearchImageWidth,
    int nSearchImageHeight,
    int nTemplateWidth,
    int nTemplateHeight,
    int& nMatchPointX,
    int& nMatchPointY
)
{
    //R(x,y)= \frac{ \sum_{x',y'} (T'(x',y') \cdot I'(x+x',y+y')) }{ \sqrt{\sum_{x',y'}T'(x',y')^2 \cdot \sum_{x',y'} I'(x+x',y+y')^2} }
    //T'(x', y')=T(x', y') - 1/(w \cdot h) \cdot \sum _{ x'',y''} T(x'',y'')
    //I'(x + x',y+y') = I(x + x',y+y') - 1 / (w \cdot h) \cdot \sum _{ x'',y'' } I(x + x'', y + y'')

//预处理T和S
    float fAvgT = 0;
    for (int i = 0; i < nTemplateHeight; i++)
    {
        for (int j = 0; j < nTemplateWidth; j++)
		{
            fAvgT += pTemplate[i * nTemplateWidth + j];
		}
    }
    fAvgT = fAvgT / (nTemplateHeight * nTemplateWidth);

    float fSigmaT = 0;

    float *fT=new float[nTemplateHeight*nTemplateWidth];

    for (int i = 0; i < nTemplateHeight; i++)
	{
		for (int j = 0; j < nTemplateWidth; j++)
		{
            fT[i * nTemplateWidth + j] = pTemplate[i * nTemplateWidth + j] - fAvgT;
			fSigmaT += fT[i * nTemplateWidth + j] * fT[i * nTemplateWidth + j];
		}
	}
//结束预处理
    int nRetHeight=nSearchImageHeight-nTemplateHeight;
    int nRetWidth=nSearchImageWidth-nTemplateWidth;
    float* m_pRet=new float[nRetHeight*nRetWidth];

    float fMaxScore = -1;

    int sze_Template = nTemplateHeight * nTemplateWidth;

    for (int i = 0; i < nRetHeight; i++)
    {
        for (int j = 0; j < nRetWidth; j++)
        {
            float fSigmaS=0;
            for (int m = 0; m < nTemplateHeight; m++)
            {
                for (int n = 0; n < nTemplateWidth; n++)
                {
                    fSigmaS += pSearchImage[(i + m) * nSearchImageWidth + j + n];
                }
            }


            float fSigmaST = 0;
            for (int m = 0; m < nTemplateHeight; m++)
            {
                for (int n = 0; n < nTemplateWidth; n++)
                {
					fSigmaST += fT[m * nTemplateWidth + n] * (pSearchImage[(i + m) * nSearchImageWidth + j + n] - fSigmaS/ sze_Template);
                }
            }

            m_pRet[i * nRetWidth + j] = fSigmaST / sqrt(fSigmaS * fSigmaT);
            if (m_pRet[i * nRetWidth + j] > fMaxScore)
			{
				fMaxScore = m_pRet[i * nRetWidth + j];
				nMatchPointX = j;
				nMatchPointY = i;
			}
        }
    }
    return fMaxScore;
}

这里的这段代码实际运行速度非常慢。下一章我将讲解如何使用FFT算法加速计算。
优化思路
根据公式不难得出NCC算法的复杂度为$$\O(n^4)$$,是一种计算复杂度很大的算法。在笔者 i9-13900HX 处理器单线程上仍需要将近12s才能完成（1000 * 1000上匹配100 * 100）。
因此需要先把复杂度降低再考虑其他方向的优化。