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

来源：互联网发布：电气画图软件编辑：程序博客网时间：2024/05/07 05:06

一下文章转载自：http://www.cnblogs.com/Imageshop/p/6108764.html，正文如下：

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

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

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左右吧，如果是灰度的匹配那就能更快了。

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

0 0