【工程應用一】多目標多角度的快速模板匹配算法（基於NCC，效果無限接近Halcon中........)

時間 2021-04-06

標籤 html 算法編程數組網絡編程語言 post 學習測試優化欄目 HTML 简体版

原文原文鏈接

願意寫代碼的人通常都不太願意去寫文章，由於代碼方面的藝術和文字中的美學每每很難兼得，二者都兼得的人一般都已經被西方極樂世界所收羅，我也是隻喜歡寫代碼，讓那些字母組成美妙的歌曲，而後自我沉浸在其中自得其樂。而今天，在清明之際，在踏青時節，我仍是忍不住停下來歇歇腳，稍微共享一下最近一直研究的一個很是基礎的算法和應用 - 多目標多角度的模板匹配。html

模板匹配，這是一個幾十年來一直爲業界所重點研究和處理的算法，存在於各類不一樣的機器視覺庫中，若是哪個沒有提供這個功能，那麼他將沒法獲取你們的承認，也就失去了最基本的活力。能夠說模板匹配基於機器視覺就至關於數組在編程語言中同樣，基礎可是不可或缺。算法

在2004年時，個人畢業設計中一個很重要的部分也是模板匹配，當時用模板匹配找到每一個量杯中黃色的油的位置，如今看來用那個算法也是醉了，不過能順利畢業還考得就是他。編程

在個人早期博客中，有一篇文章已經談到了這個算法，詳見：標準的基於歐式距離的模板匹配算法優源碼化和實現(附源代碼)，可是這個是個很是慢的過程，並且是單目標無旋轉的實現，在實際應用中，這個基本沒有啥實際的價值。數組

在工業應用場合，有着很是普遍使用場景的是多目標多角度的模板匹配（基本無縮放或輕微縮放），這方面實現的比較好的有halcon、海康、康耐視等，國內也有一些小單位有作研究，並且效果不錯。在網絡上其實也有比較多的文章談到了多目標模板匹配，基本上都是基於Opencv實現，良心的說也談到了一些核心技術，可是仍是皮毛，基本都是一帶而過，並且實現的效率也基本是沒有什麼實用價值的，多是怕說多了別人學會了吧。網絡

雖然在個人實現中，也參考了很多網絡上的文章，可是大部分的細節仍是靠的本身的思考和朋友的一些指導，爲了尊重他人，我也不打算特別深刻講解個人實現，可是仍是把一些具備必定深度的問題提出來，也算是回報網絡吧。編程語言

一、概述post

這裏先提工業界最爲經常使用，也是最爲基本的模板匹配方式，基於NCC的灰度模板匹配。學習

NCC，全稱爲Normalized Cross correlation，即歸一化互相關係數，在模板匹配中使用的很是很是普遍，也是衆多模板匹配方法中很是耀眼的存在，這個匹配的理論核心基礎公式以下：測試

（1）優化

該方法也存在於Opencv的matchTemplate中，較之其餘的CV提供的匹配方法，該算法對於光照、噪音等等的影響，穩定性更佳，也是halcon等商用軟件內嵌的基於像素的模板匹配標準方法。

他的理論匹配度範圍是[-1，1]，爲-1時表示2副圖像的極性徹底相反（原圖和反色後的圖），爲1則表示兩幅圖徹底同樣。通常咱們在計算NCC的時候都是取的絕對值，所以，一般NCC的取值爲[0,1]，值越大，表示兩幅圖像越類似。

實際編程實現時，千萬不要直接用這個公式，若是你使用，那你離砸電腦已經不遠了，請必定要相信我。

實際中，咱們都用下面的式子來實現編碼（不要問我裏面的符號的意思，兩個圖來自不一樣的資料，裏面的字母也不同，可是要研究的這個的人都應該能看懂）：

（2）

　這個式子看上去更爲複雜，可是實質上和公式（2）和公式（1）就是同一個東西。公式（2）咱們能夠把他拆解爲7個部分。咱們一1道來。

①、這個留到最後在說。

②、T表明的是模板，那麼②對於固定的模板來講就是一個定值，在匹配前能夠直接計算好，無需擔心耗時問題。

③、I 表示的搜索圖像中的和模板同樣大的一個子塊，很明顯這個累加有多重方法能夠快速的實現，好比比較原始的積分圖技術，或者個人BoxBlur裏的那種更爲快速的實現，這一項也是和參數無關的。

④、第四項處理方式同 ②項，無需多言。

⑤、第五項徹底同第二項，同時四和五項做爲一個總體也能夠提早計算好，不參與匹配過程的計算。

⑥、第六項處理方式同第三項，也無需多言。

⑦、第七項徹底同第三項，直接使用。

前面的分析代表，第二至第七項要買能夠做爲常量提早計算好，要麼就能夠經過某種技術實現O(1)的快速計算，那麼如今咱們再回過頭來在看第①項，他是模板圖像和搜索圖像同面積區域像素的一個卷積，這個是沒法用某種優化技巧去實現和模板大小無關的快速實現的，註定了他就是NCC計算式中最爲耗時的部分。

有人說卷積能夠有FFT實現優化，沒錯，很是贊成您的觀點，可是，朋友，FFT雖然其第一個F表明了Fast，可是呢他在傅里葉的世界是快的，在咱們模板匹配的空間內他受到了一種無形的壓迫，在工業界仍是沒法接受的。

所以，咱們注意到在本例中，這個卷積其實都是字節類型的計算，對於一個N*M大小的模板圖，這個卷積須要N*M次乘法以及N*M-1次加法，因爲是整數計算，自己運算速度還算比較快的，並且若是在PC上我曾經說起過有一種使用SIMD指令的提速方法，大概能有10倍的運力提升，核心是使用_mm_madd_epi16（對應PMADDWD這個指令）。

　　其原理以下：

_m128i _mm_madd_epi16(__m128i a, __m128i b)

Multiplies the 8 signed 16-bit integers from a by the 8 signed 16-bit integers from b Adds the signed 32-bit integer results pairwise and packs the 4 signed 32-bit integer results.

該指令能夠一次性執行8個16位數據的乘法及4次加法，而咱們只要提早把8爲字節數據轉換爲16爲數據就能夠了，一般這能夠由_mm_unpacklo_epi8或者_mm_ cvtepu8_epi16實現。

二、進一步提高

　有了上一步的分析，也許你就準備開始動工了，千里之行始於足下嘛，可是做爲過來人，我奉勸你除非是爲了看效果，不然，你仍是等一等吧，下面的更精彩。

若是咱們按照公式（2）就開始霸王硬上弓，大部分狀況下你須要耐心的盯着你的屏幕鼠標在那裏轉圈轉圈轉個10來秒吧，而後就能夠見證奇蹟了，也許做爲初學者，你會心動，而那些須要靠這個吃飯的，可能就是心痛了.........

爲了獲得更快的搜索速度，一個最容易想到的策略就是圖像金字塔，圖像金字塔每增長一層，圖像的點數和模板的點數救減小4倍（理論數據，實際因爲非2的寬度和高度，不必定正好是4倍），那麼不考慮Cache等其餘的因素，理論上每增長一層金字塔救能夠提速16倍，所以，若是咱們構建了一個4層的金字塔，那麼在第4層金字塔上的一次完整匹配，其計算的次數和原始的數據相比，就能減小4096倍。

咱們先以無旋轉單目標爲例進行簡單的說明，當咱們在金字塔最高層進行一次完整的匹配後，咱們能夠找一個全局的極值點，這就是在頂層匹配時的最佳匹配位置，此時，咱們能夠將頂層匹配的結果映射到金字塔的下一層中，簡單的說就是將找到的匹配點座標的X/Y值乘以2，那麼爲了保證在下一層中的精度，此時的搜索區域須要適當的增大，好比選擇匹配點周圍5*5或者7*7的一個區域，而後在這感興趣的區域中再進行一個局部的匹配計算，此時只須要計算25或者49次小匹配的計算，當計算完畢後，再次提取出這個小區域的極值，做爲本層的最終種子點，重複這個過程，直到金字塔的最底層（即原始搜索圖像）爲止。

　　稍微分析下，假定上述咱們的搜索的局部範圍爲5*5大小，金字塔取爲4層，那麼總體的計算量是原始直接計算的多少呢，這樣評價：

式中M和N分別表示圖像的寬度和高度。

通常來講M和N都是至少以百爲單位的，所以上述計算式的結果至關小，速度能夠獲得極大的提高。

這裏針僅對一個問題進行展開，即金字塔構建時採用何種下采樣算法，討論以下：

①、最近鄰，這個結果太粗糙，不利於算法穩定性，能夠直接Pass掉。

②、雙線性插值，這個兼顧速度和效果，是個能夠考慮的選項。

③、三次立方插值，這個東西在圖像放大時是個不錯的選項，而金字塔得創建是縮小過程。

④、蘭索斯插值，這個相似於三次立方插值。

⑤、採用嚴格的高斯金字塔採樣矩陣，一個5*5的矩陣，以下所示：

對於縮小，其實③④⑤都不是很好，③內部是一個4*4的取樣，④是一個8*8的取樣，我舉一個簡單的例子。

　　　　原圖（注意紅色框部分的效果）最近鄰　　　　　　　　　　雙線性　　　　　　　　三次立方　　　　　　　　　蘭索斯

你們可用大一點的屏幕去觀察，能夠看到紅色方框內在原圖部分爲很是光滑的，而四個插值中，最近鄰有所模糊，三次立方和蘭索斯在風車葉片邊緣出現了鋸齒，只有雙線性完美的保存了葉片的邊緣光滑性。

那麼通過個人實測，一種更好的方式是直接使用2*2均值下采樣，也就是使用2*2區域內的全部像素的平均值，2*2均值濾波器有一個很是好的特性，他沒有頻率響應問題，而較大的濾波器均存在該問題。同時，使用該濾波器還有一個優異的特性，即他能夠以很是高效的方式實現。

當圖像的寬度和高度都爲2的整數倍時，若是選用雙線性插值創建下一層金字塔，此時的雙線性就退化爲了2*2均值濾波器。

這裏留一個問題你們共同探討吧：

問題1：使用2*2均值濾波器時，對於非偶數的寬度和高度，好比W0 = 101,下一層金字塔的寬度W1究竟是取50，仍是取51呢，若是取51，那麼第51個像素如何獲取結果（2*2取樣會越界）呢？

三、多目標

現實中一般須要在一副原圖中尋找多個匹配的目標，此時，咱們的難度就有所上升了，金字塔咱們依舊還可使用，可是如何區分多目標呢。

在單目標時，咱們對最高層金字塔進行了全圖匹配後，只須要取最大的匹配值做爲候選點，這也是理所固然的，當問題來到多目標時，一個天然的想法就是對匹配值進行排序，而後取前N個最大值做爲候選點。可是，這種直接的想法確實很是的不科學的，由於一般在某個最大值附近，因爲鄰域的類似性，還存在不少和最大值很是接近的得分點，而他們實際上都是對應同一個目標。

在百度搜索多目標模板匹配，大部分狀況你會找到以下的一段代碼：

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; }

　　這是一段比較好的參考代碼，可是僅僅是比較好，還不能解決不少問題，可是咱們能夠從中學習到一些東西。

代碼中輸入參數中有一個參數是前一次的最小值的座標，而後在這個座標附近必定的矩形範圍內（上述代碼是模板圖像的1/3尺寸），將得分值修改成某個很大的值，接着再進行全範圍的最大和最小值定位，此時，確定就定位到了離輸入最小值有必定距離的另一個最小值，而輸入最小值附近的那些次小值就不會干擾結果。

注意上面代碼是最小值，由於他用的檢測指標是CV_TM_SQDIFF_NORMED，而非NCC，對於NCC，則須要歸爲最大值。

這裏的templatW / 3和templatW / 3有點重疊範圍的意思了，可是還徹底不夠。

在最頂層金字塔中找到了多個目標的粗糙位置後，就能夠和前所述的同樣的方式一步一步的向下一層金字塔進行細化，直處處理到頂層金字塔爲止。

四、多目標+多角度

當問題來到這裏時，整個的難度就是階躍式的提升了一個檔次。

若是目標存在旋轉，爲了能找到發生旋轉的物體，咱們能夠建立多個方向的旋轉對象，也就是說，將搜索空間離散化，此時，有兩個可選的方式：一個是旋轉搜索圖像，而後用模板在旋轉後的圖像中搜索，二是旋轉模板，用旋轉後的模板在搜索圖像中定位。咱們說，第一種方式基本不可取，緣由有三。

(1)、搜索圖像通常來講都是較大的圖，對其進行旋轉耗時比較可觀。、

(2)、實際狀況須要多個角度的旋轉，對原圖旋轉內存方面也會有過多的消耗

(3)、工業應用時，通常模板比較固定，而搜索圖像老是時刻變化的。

當選擇第二種方法時，對於較小的模板圖像，是能夠在執行搜索前把相關旋轉信息提早準備好，在搜索時刻直接使用，而無需作無謂的耗時。

此時，在金字塔的最頂層，須要作的計算工做也有所增長，咱們須要對每一個角度的模板都作一個全圖的匹配，獲得匹配的結果，而後對每一個可能點，選擇匹配度最大的那個角度做爲頂層的候選點。

相似的，在向下一層金字塔映射時，不只僅須要映射匹配點的X和Y座標，還須要映射角度信息，同理，爲了保證角度方面的精度，也須要適當的擴大角度的搜索區間。

五、更多的問題

實際上，爲了實現多目標和多角度的匹配，還存在不少不少的細節問題，須要取研究，這些方面的細節有些已經有了部分結論，可是大部分在網絡上鮮有探討和方向，這裏列出一些問題和你們共同探討一番。

問題2：金字塔多少層比較合適？

前面提到了金字塔層數越多，所需的計算量就越少，可是同時帶來一個問題，就是計算的精度會愈來愈差，這是由於，隨着金字塔層數的增長，由於二次採樣圖像會不斷的獲得平滑，在圖像分辨率變得過低時一些基本的特徵已經徹底丟失，各類原本不想關的信息已經徹底融合在一塊兒了。所以，一個合適的金字塔層數就尤其重要。一些成熟的商業軟件通常可提供用戶自行輸入金字塔層數，或者自動肯定。對於大部分用戶而言，他們不懂得更多的技術細節，自動設置顯得尤其重要。

咱們以工業界最爲出色Halcon軟件爲例，通過屢次測試，他的金字塔層數的自動設置是很是智能化的，基本上自動能夠保證速度最優同時效果穩定，一般，咱們認爲金字塔的層數只和模板圖像的尺寸有關，可是，一些仔細的測試代表，兩幅一樣大小的模板圖，halcon也有可能但不必定返回不一樣的值，估計這其中還用到了圖像的一些類似性或者結構方面的一些信息做出的綜合判斷，不過，一個最基本的規律，仍是能夠分享下，那就是若是金字塔某層的短邊像素小於或者等於8了，大家這確定是此模板金字塔層數的極限了。

問題3：角度離散化的間距如何設置最爲合理？

相似於金字塔，角度離散化的間距越大，須要計算的旋轉方面的信息就越少，速度則能夠更快，可是一樣，精度就越差。

通常來講，若是模板越大，離散化的間距則須要越小，這是由於較大的模板可以區別更小角度的變化。一般，對於大小100像素的模板，離散的角度步幅可設置爲1度。一樣，如何自動的設置參數也應該是一個成熟軟件的標誌，一個可行的建議是離散的步長鬚要保證旋轉時兩個相鄰模板之間的長或寬的最小差別值不小於1個像素。

問題4：模板的旋轉如何處理？

一般能夠用雙線性插值或者三次立方插值，來獲取旋轉後的數據，不建議用最近鄰算法，可是不一樣的旋轉算法，最後獲得的匹配結果會有所不一樣，同時這也就說呢，其實帶角度的模板匹配，理論上很難獲取精確解，由於你畢竟不知道原始的旋轉算法是何種，好比我從一個未旋轉圖像中扣一個矩形出來做爲模板，而後把圖像旋轉各10度，用halcon對模板進行匹配，獲得的結果哪怕選擇亞像素也不會是精確的10度。

問題5：旋轉後無效數據的部分怎麼辦？

在對模板進行旋轉時，除非是幾個特殊的角度，好比0、90、270、360度，不會產生額外的信息，其餘的角度，都會有一些未知的區域存在，以下圖所示：

　　　　　　原圖旋轉必定的角度（雙線性插值）邊緣處局部放大

中間爲旋轉後的模板，紅色部分爲新出現的未知區域，若是說咱們這個圖做爲一個總體，去和原圖像進行匹配，也就是說讓紅色部分參與了匹配，這很明顯會獲得一個得分較低的錯誤結果。

一種也許可行的方式是，咱們把這些紅色的區域剔除在匹配有效的數據以外，這時，又會面臨另一個新的問題，在圖像的邊緣部分如何處理。在上面邊緣處局部的方法圖中，咱們能夠看到，因爲插值的特性，邊緣處未能在原始圖像中採集到足夠的採樣點，所以，選擇了紅色背景色做爲融合的基色，此時的結果像素就不徹底是屬於原始圖像了，怎麼辦？我的的一個建議是既然這部分像素也被污染了，那就不該該參與後續的匹配。

　此時，在編程方面就須要繼續克服下一個困難了，前面概述方面講的一些優化方式又不能直接使用了，真是他媽的痛苦，因此，借用某一本書裏的經典勸退名言吧：聰明的機器視覺用戶都依賴標準軟件包來提供這些功能而不會試着本身實現這些算法。

問題6：各層金字塔的角度離散值如何分配？

一般，在金子塔的最底層（和原圖同樣大小那一層），可按照前述的自動角度幅值來一步一步的旋轉圖像，而後隨着金字塔的層數增長，根據模板在每層金字塔中都會縮小2倍的這個事實，在相應金字塔上模板的角度步幅也能夠增長2倍，所以，若是在金字塔最底層上使用的角度步幅喂1度，那麼在第四層上可使用8度做爲角度步幅。

這樣作帶來的好處有不少，由於，一般，咱們須要在最頂層作多角度的全圖的匹配，這個計算的相對來講比較大，若是角度步幅在該層得以以指數級別減小，則計算量也會有量級的下降。

問題7：各層金字塔的最低得分值如何肯定？

前面一直沒提這個得分的問題，在單目標中通常是不存在這個問題的，當有多個目標時，由於可能不是完整的匹配，所以通常須要客戶肯定一個最小的得分值，固然這個得分是指的在最底層時的值。當有了金字塔後，由於下采樣的一些因素，爲了有更高的容錯率，通常是建議每增長一層金字塔，最小的得分值須要適當的下降。好比，用戶設置的最小得分爲0.8，後續各層的得分能夠爲：

　　　　 0.8 　　0.8*0.9　　　　0.8*0.9*0.9　　　　　0.8*0.9*0.9 * 0.9　　

　　固然，其餘的調整方式也何嘗不能夠，可是，整體的區域需保證最小得分愈來愈低。　　

　　問題8：MaxOverlap是什麼鬼，內部是如何操做的？

在Halcon中，還有很是重要的參數-MaxOverlap，一個介於0和1之間的參數，前面也一直沒有談到過。其實，在真正的應用中，存在着一些目標之間有必定程度重疊的狀況，這個時候，若是按照前面的那些處理方式，通常就只能獲取到重疊對象的某一個，而丟掉了其餘的信息。當有了MaxOverlap參數時，咱們就能夠根據這些對象的重疊重複，來決定是否剔除掉某一個不須要的目標。

可是，這也是一個比較難以琢磨的對象，Halcon幫助文檔中的說法是取的對消的最小外接矩形中的重疊比例。說實在的，編程時，這個規則應該還不夠。好比，若是有4個對象，他們都有所重疊，請問這個時候，這個MaxOverlap是指的全部的重疊量合併在一塊兒呢，仍是最大的重疊，抑或是按照得分順序第一個和之重疊的呢。目前，我也是對這個比較搖擺。

同時，另一個比較難以肯定問題是，這個重疊是在金字塔的最頂層就進行判斷仍是如何呢，若是在金字塔的最頂層不進行判斷，那麼金字塔頂層中得分大於MinScore的則有不少個點，若是把這些點都直接向上一層進行映射，那個計算量也是至關客觀的。

問題9：亞像素座標和角度是一塊兒執行的嗎，仍是分開的？

沒有亞像素的模板匹配是沒有靈魂的，特別是帶有角度的匹配。由於，正如前面所述，咱們對角度採用了離散化。那麼這個時候計算的最終角度結果必然是離散化後的序列裏的某個值，這樣的精度有的時候是不夠的。

經過亞像素技術，尋找局部區域裏的最大值，課得到更高的位置精度和角度精度，則能夠有效的得到精準的定位。

可是目前仍是不清楚座標+角度組成的4維的空間的亞像素計算公式是如何的。

問題10：速度優化？

以上談的一些優化基本上是結構上的優化，或者說原理上的優化，在編碼上還能夠考慮時候用SIMD指令進行處理，尤爲是全圖匹配的計算過程。

六、公關成果

大概前先後後再一塊兒有折騰了大概一年的時間吧，初步仍是有了必定的成果，不管是在速度上，仍是在準確度中，仍是能解決必定的工程問題了。

共享一些測試圖，和你們一塊兒比較。

　　模板圖　　　　　　　　　　　　　　　　　待搜索圖　　　　　　　　　　　　　　　　　　　　　　　　　　搜索結果

基於NCC的速度方面和不少參數有關，好比角度的搜索範圍，金字塔層數，模板的大小（通常模板大，速度反而快些，特別小的模板則很是耗時，知道緣由嗎？），重疊的大小等等都有關。一些簡單的測試數據以下：