标准的基于欧式距离的模板匹配算法优源码化和实现(附源代码)。

     很久没有出去溜达了,今天天气好,就放松放松去,晚上在办公室没啥事,把以前写的一个基于标准的欧式距离的模板匹配代码共享吧。

     opencv有模板匹配的代码,我没看他是如何优化的,所以不管他吧,我只描述我自己实现。

     基于欧式距离的模板匹配就是遍历被匹配图的每一个像素,然后计算以该像素为中心,和模板图重叠部分的像素的欧式距离,当模板图越大时,计算就急剧增加,因此做优化才能有真正的实用价值。

     两个标量的欧式距离表达式为 (a - b) * (a - b),展开后为 a^2 + b^ 2 - 2ab,我们每一个像素点的计算就是WM * HM个像素色阶值的距离的累加和(WM和HM分别为模板图的宽度和高度),模板匹配中,模板图所有像素的平方和是固定的,可以提前计算,而被匹配图中每个像素点周边WM * HM的像素的平方和可以使用类似BoxBlur中懒惰算法快速的得到,而只有两者的成绩项是必须每个点重新计算,这也是整个计算过程中最为耗时的部分,如果直接用C的代码写出来,恐怕等到花儿都谢了。

     我在图像处理中任意核卷积(matlab中conv2函数)的快速实现一文中曾经给出过一种基于SSE的的快速卷积的算法,他可以一次性计算出16个字节的乘法,速度因此也得到了大的提升,因此,完全可以用在上述的计算a * b的过程中,这样我们的模板匹配速度就能有质的提高。

    计算模板图的像素自乘平方和代码非常简单,也没啥耗时,简单代码如下:

int GetPowerSum(TMatrix *Src)            //    无需注释
{
    if (Src == NULL || Src->Data == NULL) return 0;
    if (Src->Depth != IS_DEPTH_8U) return 0;

    int X, Y, Sum, Width = Src->Width, Height = Src->Height;
    unsigned char *LinePS;
    
    if (Src->Channel == 1)
    {
        for (Y = 0, Sum = 0; Y < Height; Y++)
        {
            LinePS = Src->Data + Y * Src->WidthStep;
            for (X = 0; X < Width; X++)
            {
                Sum += LinePS[X] * LinePS[X];
            }
        }
    }
    else
    {
        for (Y = 0, Sum = 0; Y < Height; Y++)
        {
            LinePS = Src->Data + Y * Src->WidthStep;
            for (X = 0; X < Width; X++)
            {
                Sum += LinePS[0] * LinePS[0] + LinePS[1] * LinePS[1] + LinePS[2] * LinePS[2];
                LinePS += 3;
            }
        }
    }
    return Sum;
}  而计算被匹配图中每个像素为中心,WH*WM范围内像素的自乘平方和的O(1)算法也比较简单:
/// <summary>
/// 计算图像的局部平方和,速度已经优化,支持1和3通道图像。(2015.10.5日)
/// </summary>
/// <param name="Src">待求平方和的图像。</param>
/// <param name="Dest">平方和数据,需要使用int类型矩阵保存,大小为Src->Width - SizeX + 1, Src->Height - SizeY + 1,程序内部分配数据。</param>
/// <param name="SizeX">在水平方向使用的模板大小,如果是半径模式,对应的量为2 * Radius + 1。</param>
/// <param name="SizeY">在垂直方向使用的模板大小,如果是半径模式,对应的量为2 * Radius + 1。</param>
/// <remarks> 1:使用了类似BoxBlur里的优化算法,耗时和参数基本无关。</remarks>
/// <remarks> 2:也可以使用积分图实现。</remarks>

IS_RET GetLocalSquareSum(TMatrix *Src, TMatrix **Dest, int SizeX, int SizeY)
{
    if (Src == NULL || Src->Data == NULL) return IS_RET_ERR_NULLREFERENCE;
    if (Src->Depth != IS_DEPTH_8U || Src->Channel == 4) return IS_RET_ERR_NOTSUPPORTED;
    if (SizeX < 0 || SizeY < 0) return IS_RET_ERR_ARGUMENTOUTOFRANGE;
    
    int X, Y, Z, SrcW, SrcH, DestW, DestH, LastIndex, NextIndex, Sum;
    int *ColSum, *LinePD;    
    unsigned char *SamplePS, *LastAddress, *NextAddress;
    IS_RET Ret = IS_RET_OK;

    SrcW = Src->Width, SrcH = Src->Height;
    DestW = SrcW - SizeX + 1, DestH = SrcH - SizeY + 1;

    Ret = IS_CreateMatrix(DestW, DestH, IS_DEPTH_32S, 1, Dest);                                
    if (Ret != IS_RET_OK) goto Done;
    ColSum = (int*)IS_AllocMemory(SrcW * sizeof(int), true);
    if (ColSum == NULL) {Ret = IS_RET_ERR_OUTOFMEMORY; goto Done;}

    if (Src->Channel == 1)
    {
        for (Y = 0; Y < DestH; Y++)
        {
            LinePD = (int *)((*Dest)->Data + Y * (*Dest)->WidthStep);
            if (Y == 0)
            {
                for (X = 0; X < SrcW; X++)
                {
                    Sum = 0;
                    for (Z = 0; Z < SizeY; Z++)
                    {
                        SamplePS = Src->Data + Z * Src->WidthStep + X;
                        Sum += SamplePS[0] * SamplePS[0] ;
                    }
                    ColSum[X] = Sum;
                }
            }
            else
            {
                LastAddress = Src->Data + (Y - 1) * Src->WidthStep;                    
                NextAddress = Src->Data + (Y + SizeY - 1) * Src->WidthStep;            
                for (X = 0; X < SrcW; X++)
                {
                    ColSum[X] -= LastAddress[X] * LastAddress[X] - NextAddress[X] * NextAddress[X];
                }
            }
            for (X = 0; X < DestW; X++)
            {
                if (X == 0)
                {
                    Sum = 0;
                    for (Z = 0; Z < SizeX; Z++)
                    {
                        Sum += ColSum[Z];
                    }
                }
                else
                {
                    Sum -= ColSum[X - 1] - ColSum[X + SizeX - 1];
                }
                LinePD[X] = Sum;
            }
        }
    }
    else if (Src->Channel == 3)
    {
        for (Y = 0; Y < DestH; Y++)
        {
            LinePD = (int *)((*Dest)->Data + Y * (*Dest)->WidthStep);
            if (Y == 0)
            {
                for (X = 0; X < SrcW; X++)
                {
                    Sum = 0;
                    for (Z = 0; Z < SizeY; Z++)
                    {
                        SamplePS = Src->Data + Z * Src->WidthStep + X * 3;            //    三通道累加到一起
                        Sum += SamplePS[0] * SamplePS[0] + SamplePS[1] * SamplePS[1] + SamplePS[2] * SamplePS[2];
                    }
                    ColSum[X] = Sum;
                }
            }
            else
            {
                LastAddress = Src->Data + (Y - 1) * Src->WidthStep;                    
                NextAddress = Src->Data + (Y + SizeY - 1) * Src->WidthStep;    
                for (X = 0; X < SrcW; X++)
                {
                    ColSum[X] += NextAddress[0] * NextAddress[0] + NextAddress[1] * NextAddress[1] + NextAddress[2] * NextAddress[2] - LastAddress[0] * LastAddress[0] - LastAddress[1] * LastAddress[1] - LastAddress[2] * LastAddress[2];
                    LastAddress += 3;
                    NextAddress += 3;
                }
            }
            for (X = 0; X < DestW; X++)
            {
                if (X == 0)
                {
                    Sum = 0;        
                    for (Z = 0; Z < SizeX; Z++)
                    {
                        Sum += ColSum[Z];
                    }
                }
                else
                {
                    Sum -= ColSum[X - 1] - ColSum[X + SizeX - 1];
                }
                LinePD[X] = Sum;
            }
        }
    }
Done:
    IS_FreeMemory(ColSum);
    return Ret;
}
  上述代码思路类似于BoxBlur的实现方式,如果还想更快点,可以参考解析opencv中Box Filter的实现并提出进一步加速的方案(源码共享)一文的基于SSE的处理方式,有兴趣的朋友可以自研。

       其实速度也不快,但是有些应用场合模板图很小(比如16*16的),被匹配图也不大,比如640 * 480的,这个时候大概也就30ms左右吧,如果是灰度的匹配那就能更快了。

       其实代码如果想优化,还是可以用线程并行的。

      代码下载:http://files.cnblogs.com/files/Imageshop/MatchTemplate.rar(解压密码: Buy me a beer)

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

发表于

我来说两句

0 条评论
登录 后参与评论

相关文章

来自专栏程序员叨叨叨

8.4 CG 标准函数库

和 C 的标准函数库类似,Cg 提供了一系列内建的标准函数。这些函数用于执行数学上的通用计算或通用算法(纹理映射等),例如,需要求取入射光线的反射光线方向向量可...

1225
来自专栏aCloudDeveloper

算法导论第四章分治策略剖根问底(二)

   在上一篇中,通过一个求连续子数组的最大和的例子讲解,想必我们已经大概了然了分治策略和递归式的含义,可能会比较模糊,知道但不能用语言清晰地描述出来。但没关系...

2076
来自专栏张俊红

机器学习中的参数调整

总第102篇 前言 我们知道每个模型都有很多参数是可以调节的,比如SVM中使用什么样的核函数以及C值的大小,决策树中树的深度等。在特征选好、基础模型选好以后我们...

3987
来自专栏数据小魔方

Stata特别篇(下)——多变量图表汇总!

今天跟大家分享Stata特别篇的下篇——多变量图表汇总! 在多变量图表中,增加的变量仅仅限于定距变量,也可以是定类变量。 打开数据集: use "D:\Sta...

4565
来自专栏SeanCheney的专栏

Numpy和MatplotlibPython科学计算——Numpy线性代数模块(linalg)随机模块(random)Python的可视化包 – Matplotlib2D图表3D图表图像显示

Python科学计算——Numpy Numpy(Numerical Python extensions)是一个第三方的Python包,用于科学计算。这个库的前身...

4374
来自专栏大数据文摘

改变计算技术的 9 个伟大算法

1943
来自专栏机器学习之旅

GolVe向量化做文本分类向量化文本分类

第一种是常规方法的one-hot-encoding的方法,常见的比如tf-idf生成的0-1的稀疏矩阵来代表原文本:

1624
来自专栏人工智能LeadAI

GolVe向量化做文本分类

第一种是常规方法的one-hot-encoding的方法,常见的比如tf-idf生成的0-1的稀疏矩阵来代表原文本:

1113
来自专栏aCloudDeveloper

算法导论第十五章 动态规划

      写在前面:从本章开始,算法导论章节进入第四部分:高级设计和分析技术。在读的过程中,可以明显感觉到本章内容跟之前章节的内容要复杂得多。这么来说,之前章...

1905
来自专栏zhisheng

学习算法之路

一个搞ACM的需要掌握的算法的sheet。 第一阶段:练经典常用算法,下面的每个算法给我打上十到二十遍,同时自己精简代码,因为太常用,所以要练到写时不用想,10...

3685

扫码关注云+社区