標準的基於歐式距離的模板匹配算法優源碼化和實現(附源代碼)。

     好久沒有出去溜達了，今每天氣好，就放鬆放鬆去，晚上在辦公室沒啥事，把之前寫的一個基於標準的歐式距離的模板匹配代碼共享吧。

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

       其實代碼若是想優化，仍是能夠用線程並行的。

a

　　    代碼下載：http://files.cnblogs.com/files/Imageshop/MatchTemplate.rar（解壓密碼： Buy me a beer）