【工程應用六】 繼續聊一聊高效率的模板匹配演算法（分水嶺助威+蒙版提速）。
總是寫很長的複雜的文章，目前發現真的有點無法靜心去弄了，感覺寫程式碼的動力要比寫文章強大的多，所以，往後的文章還是寫的剪短一點吧。
繼續聊一聊模板匹配。 最近這方面也出了一些新的資料，說明還是有人關注他的。
我最近一個月的研究成果主要有以下幾個方面。
一、頂層金字塔的候選點選擇改由分水嶺相關演算法實現（用時10天）。
頂層的金字塔，我們是全圖計算相關得分值的。當計算完所有的頂層金字塔得分後，我們得到了不同角度不同位置的一個全方位的候選點資訊，接下來我們的目標就是從這些點中選擇合適的候選點。
這裡有幾個指標可以作為初步篩選的依據：
1、最小的得分值。
2、重疊的區域。
在【工程應用一】 多目標多角度的快速模板匹配演算法（基於NCC，效果無限接近Halcon中……..) 一文中，我曾分享過如下的程式碼：
Point getNextMinLoc(Mat &result, Point minLoc, int maxValue, int templatW, int templatH) { int startX = minLoc.x - templatW / 3; int startY = minLoc.y - templatH / 3; int endX = minLoc.x + templatW / 3; int endY = minLoc.y + templatH / 3; if (startX < 0 || startY < 0) { startX = 0; startY = 0; } if (endX > result.cols - 1 || endY > result.rows - 1) { endX = result.cols - 1; endY = result.rows - 1; } int y, x; for (y = startY; y < endY; y++) { for (x = startX; x < endX; x++) { float *data = result.ptr<float>(y); data[x] = maxValue; } } double new_minValue, new_maxValue; Point new_minLoc, new_maxLoc; minMaxLoc(result, &new_minValue, &new_maxValue, &new_minLoc, &new_maxLoc); return new_minLoc; }
他通過不斷的迭代，每次以剩餘數據中最大值為候選點，並且逐步去除部分領域的方法來獲取候選點。 我所能看到的開源項目里基本已這個程式碼為藍圖來實現他。
這個方法是可行的，我一直在用，但是他也是有缺陷的，當模板比較小的時候，我們的金字塔層數不夠多，這個時候這個函數本身的計算的耗時就較為明顯了，而且還有一個問題，就是他會返回相對來說較多的候選點。造成後續的進一步篩選的計算量加大。
我一直在尋找更為科學的辦法，直到最近偶爾在一個地方看到一個這樣的演算法效果。
這不正是我們頂層金字塔需要的演算法嗎？
我嘗試把幾個測試圖的頂層金字塔的得分數轉換為影像，分別如下所示：
可以看到，他們都是類似的這種有局部最亮點的影像，那如何用演算法實現呢，後來我在ImageJ里發現一個功能（如上圖所示介面的Process菜單下的FindMaxma），基本就是這個功能的翻版：
於是我就去ImageJ里找這個演算法的程式碼，在MaximumFinder.java里找到了相關的資料，程式碼有1300多行，說垛也不多，不過我去描了一下，還是過於複雜了，關鍵是沒有相關參考文章，無法理解其程式碼的意義，不過一個核心的意思就是利用了分水嶺演算法，並且ImageJ里的一些二值分割演算法里也用到這個。
知道了他是用的分水嶺演算法，那就好辦了，我同樣在ImageJ的網站了找到了這個：//imagej.nih.gov/ij/plugins/watershed.html，他提供了最原始的分水嶺實現程式碼，對應文章為：The Watershed Transform: Definitions, Algorithms and Parallelization Strategies”。裡面了用到了一些特殊的結構。在github上還可以找到一個對應的C版本的程式碼，不過那個裡面有很多delete *p,建議刪除，很影響速度。
我憑著我的聰明才智，把哪些什麼Queue, List等等複雜的數據結構體都拋棄不用，即提高了程式碼速度，也減少了記憶體佔用量，也基本實現了這個演算法，而用在外面這裡，需要先把得分的數據整形化，我測試把他劃分為為1024個整形，應該就足夠了。
分水嶺的計算過程把影像分成一個一個的分開的塊，外面有了塊的標記後，選取每個塊的最大值作為候選點的位置和得分值即可。
說起來簡單，但是做起來難啊， 前前後後我這折騰這個過程，也用了10來天時間，結果就是對於大部分測試圖，整體速度有40%的提高（因為頂層找到的候選點少了，同時找候選點的時間也少一些），對於那些本身目標就非常明確的圖，區別不大，對於模板很小的圖，頂層計算佔據了很大耗時的影像，速度有200%的提高。
另外，單獨說明一點，在我的測試中，僅僅依據最大得分選擇候選點有可能會丟失一些目標，核心原因是頂層金字塔的角度量化方面可能到底局部得分偏向於某一側，解決辦法是檢測通過分水嶺獲得的最大值周邊3*3領域的點的角度和最大值處的角度的差異，如果差異明顯，則周邊的點最好也納入候選點系列。
二、增加形狀匹配的蒙版功能（4天搞定）
形狀匹配的準確性和提取到的形狀邊緣運算元的精確度有著很大的關係，在有些應用中，我們選擇的模板可能有部分區域的邊緣特徵是不需要的，或者模板有部分噪音過於嚴重，會對檢測結果有很大的影響，這時候帶有蒙版功能的形狀匹配就非常有必要了。
要實現這個功能，理論上來說是不複雜的，只要把哪些處於非有效區域的邊緣特徵點剔除掉即可。但是在實際的編碼過程中，還是有幾點要注意：
1、我們需要為蒙版影像也創建金字塔，那麼客戶提供的蒙版一般為二值圖，在創建金字塔的過程中，因為是2*2插值縮放，必然會產生非二值的像素結果，處理辦法是放鬆這個結果，只要處理的結果大於0，則賦值為白色，否則為黑色，如果不做這個處理，或用普通的127作為二值的閾值，調試發現會丟失目標。
2、即使是這樣，可能還不夠，在上述二值化的基礎上，最好還對邊緣進行一次半徑為1的Dilate。
其實有了蒙版不是壞事，雖然在創建模型的時候速度會慢一些，但是後續因為特徵點的減少，這個查找目標的速度反而會快一點，比如下面這個莫版圖，我們最關係的其實周邊的橢圓的形狀，而橢圓內部有什麼我們不在乎，所以增加了一個蒙版，在沒有增加蒙版前共有13510 個邊緣點（其實大部分邊緣點都是無意義的點，設置蒙版後，只強調了有用的邊緣資訊，只有2828個特徵點了。這也就意味著只需要匹配更少的特徵數據了，速度和精確度都有所提高。
模板圖 蒙版圖 原始模板的特徵點 帶蒙版的特徵點
三、蒙版可為空版本+對比度自動設置（2天）
有的時候可能還是不需要蒙版的，所以這個函數還是要考慮這個功能，即傳入空指針就調用沒有蒙版的函數。
另外，基於形狀的匹配有個對比度和最小對比度的參數，一般客戶還是希望自動化，這裡取個簡單的演算法，直接用模板影像的OSTU二值化的那個參數作為對比度的值，最小對比度取其1/2或者1/4吧。
四、基於NCC版本的蒙版功能（廢了我7天）
原來一直搞不定這個演算法，主要是因為不曉得如何快速的計算中間的有些函數了，後來還是想起來那個行程編碼，還是搞定了，具體的實現其實還是在上我上一篇博文之間就已經實現了，當然現在也可以借用那個博文來描述本演算法的過程，詳見：超越OpenCV速度的MorphologyEx函數實現（特別是對於二值圖，速度是CV的4倍左右），只是那裡是求最大值，這裡是求累加值或者平方累加值。
用類似如下程式碼更改即可：
StdT = sqrtf(IM_Max((PowSumT - (float)SumT * SumT / ValidPixel) / ValidPixel, 0.000001)); // 模板的均方差,IM_Max還是主要為了防止計算精度有誤差，導致小於0的情況出現 float *SumST = (float *)malloc(ValidW * ValidH * sizeof(float)); // 模板和原圖相乘的卷積和 int *Sum = (int *)calloc(ValidCol * ValidW, sizeof(int)); // 保存一行像素每個有效列的累計值 int *SumP = (int *)calloc(ValidCol * ValidW, sizeof(int)); // 保存一行像素每個有效列的累計平方值 int *SumS_X = (int *)calloc(ValidW, sizeof(int)); float *PowSumS_X = (float *)calloc(ValidW, sizeof(float)); if ((SumST == NULL) || (Sum == NULL) || (SumP == NULL) || (SumS_X == NULL) || (PowSumS_X == NULL)) { Status = IM_STATUS_OUTOFMEMORY; goto FreeMemory; } // 注意這裡沒有使用Mask的資訊，是因為Template之中Mask為0的地方像素值也為0，因此相乘還是0，所以可以不用帶入Mask。 Status = IM_FastConv2(Src, Template.Bmp, SumST); // 整體計算出乘積卷積，缺點是耗用了記憶體，優點是能夠能提高一定的速度 if (Status != IM_STATUS_OK) goto FreeMemory; int BlockSize = 8, Block = ValidW / BlockSize; for (int Y = 0; Y < ValidH; Y++) { float *LinePN = NCC + Y * ValidW; float *LineST = SumST + Y * ValidW; if (Y == 0) // 第一行第一個點，完整的進行計算 { for (int XX = 0, Index = 0; XX < TemplateW; XX++) { if (RL_V[XX].Amount != 0) { int *LinePS = Sum + Index * ValidW; int *LinePP = SumP + Index * ValidW; for (int Z = 0; Z < RL_V[XX].Amount; Z++) { unsigned char *LinePT = Src.Data + RL_V[XX].SE[Z].Start * StrideS + XX; for (int YY = RL_V[XX].SE[Z].Start; YY <= RL_V[XX].SE[Z].End; YY++) { for (int X = 0; X < Block * BlockSize; X += BlockSize) // 對速度基本沒有啥影響 { __m128i SrcV = _mm_cvtepu8_epi16(_mm_loadl_epi64((__m128i *)(LinePT + X))); // 8個16位數據 __m128i PowerV = _mm_mullo_epi16(SrcV, SrcV); // 8個16位數據的平方 _mm_storeu_si128((__m128i *)(LinePS + X + 0), _mm_add_epi32(_mm_loadu_si128((__m128i *)(LinePS + X + 0)), _mm_cvtepu16_epi32(SrcV))); _mm_storeu_si128((__m128i *)(LinePS + X + 4), _mm_add_epi32(_mm_loadu_si128((__m128i *)(LinePS + X + 4)), _mm_cvtepu16_epi32(_mm_srli_si128(SrcV, 8)))); _mm_storeu_si128((__m128i *)(LinePP + X + 0), _mm_add_epi32(_mm_loadu_si128((__m128i *)(LinePP + X + 0)), _mm_cvtepu16_epi32(PowerV))); _mm_storeu_si128((__m128i *)(LinePP + X + 4), _mm_add_epi32(_mm_loadu_si128((__m128i *)(LinePP + X + 4)), _mm_cvtepu16_epi32(_mm_srli_si128(PowerV, 8)))); } for (int X = Block * BlockSize; X < ValidW; X++) { LinePS[X] += LinePT[X]; LinePP[X] += LinePT[X] * LinePT[X]; } LinePT += StrideS; } } Index++; } } }
那麼目前演算法研究到這一步，其實我後續一直想攻克的就是形狀模型的創建速度和模型文件的大小問題，在Halcon中，我們會發現形狀模型創建的速度特別快，而且模型文件也非常小。內部的機理我想無非就是他是在創建時只保存了為旋轉和縮放的模板的不同金字塔層的特徵，然後在匹配的時候進行特徵的旋轉。 而我們現在都是創建的時候旋轉影像，然後再計算出個角度的特徵。這個計算量就特別大了，如果同時還考慮縮放，那基本上模板圖稍微大一點，就會造成速度奇慢和記憶體暴漲。所以我現在一直沒有做即帶旋轉有帶縮放的匹配，不是技術上實現不了，而是實現的實用性夠嗆。
目前，關於這個，我也一直在構思，是不是可以通過亞像素的canny來實現類似的功能呢，期待吧，也許將來補救就會有突破，相信自己。
最新版的一個測試DEMO: 帶蒙版的模板匹配。
如果想時刻關注本人的最新文章，也可關注公眾號或者添加本人微信 laviewpbt
