在工業應用場合，有着非常廣泛使用場景的是多目标多角度的模闆比對（基本無縮放或輕微縮放)，本人嘗試基于NCC模仿halcon的一些方法和步驟，加以複現，基本實作了具有一定價值的快速比對。本文對實作過程中一些細節予以描述。
願意寫代碼的人一般都不太願意去寫文章，因為代碼方面的藝術和文字中的美學往往很難兼得，兩者都兼得的人通常都已經被西方極樂世界所收羅，我也是隻喜歡寫代碼，讓那些字母組成美妙的歌曲，然後自我沉浸在其中自得其樂。而今天，在清明之際，在踏青時節，我還是忍不住停下來歇歇腳，稍微共享一下最近一直研究的一個非常基礎的算法和應用 - 多目标多角度的模闆比對。
模闆比對，這是一個幾十年來一直為業界所重點研究和處理的算法，存在于各種不同的機器視覺庫中，如果哪一個沒有提供這個功能，那麼他将無法擷取大家的認可，也就失去了最基本的活力。可以說模闆比對基于機器視覺就相當于數組在程式設計語言中一樣，基礎但是不可或缺。
在2004年時，我的畢業設計中一個很重要的部分也是模闆比對，當時用模闆比對找到每個量杯中黃色的油的位置，現在看來用那個算法也是醉了，不過能順利畢業還考得就是他。
在我的早期部落格中，有一篇文章已經談到了這個算法，詳見：标準的基于歐式距離的模闆比對算法優源碼化和實作(附源代碼)， 但是這個是個非常慢的過程，而且是單目标無旋轉的實作，在實際應用中，這個基本沒有啥實際的價值。
在工業應用場合，有着非常廣泛使用場景的是多目标多角度的模闆比對（基本無縮放或輕微縮放），這方面實作的比較好的有halcon、海康、康耐視等，國内也有一些小機關有做研究，而且效果不錯。在網絡上其實也有比較多的文章談到了多目标模闆比對，基本上都是基于Opencv實作，良心的說也談到了一些核心技術，但是還是皮毛，基本都是一帶而過，而且實作的效率也基本是沒有什麼實用價值的，可能是怕說多了别人學會了吧。
雖然在我的實作中，也參考了不少網絡上的文章，但是大部分的細節還是靠的自己的思考和朋友的一些指導，為了尊重他人，我也不打算特别深入講解我的實作，但是還是把一些具有一定深度的問題提出來，也算是回報網絡吧。
1、概述
這裡先提工業界最為常用，也是最為基本的模闆比對方式，基于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這個指令）。
其原理如下：
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、進一步提升
有了上一步的分析，也許你就準備開始動工了，千裡之行始于足下嘛，但是作為過來人，我奉勸你除非是為了看效果，否則，你還是等一等吧，下面的更精彩。
如果我們按照公式（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取樣會越界）呢？
3、多目标
現實中通常需要在一副原圖中尋找多個比對的目标，此時，我們的難度就有所上升了，金字塔我們依舊還可以使用，但是如何區分多目标呢。
在單目标時，我們對最高層金字塔進行了全圖比對後，隻需要取最大的比對值作為候選點，這也是理所當然的，當問題來到多目标時，一個自然的想法就是對比對值進行排序，然後取前N個最大值作為候選點。但是，這種直接的想法确實非常的不科學的，因為通常在某個最大值附近，由于鄰域的相似性，還存在很多和最大值非常接近的得分點，而他們實際上都是對應同一個目标。
在百度搜尋多目标模闆比對，大部分情況你會找到如下的一段代碼：
這是一段比較好的參考代碼，但是僅僅是比較好，還不能解決很多問題，但是我們可以從中學習到一些東西。
代碼中輸入參數中有一個參數是前一次的最小值的坐标，然後在這個坐标附近一定的矩形範圍内（上述代碼是模闆圖像的1/3尺寸），将得分值修改為某個很大的值，接着再進行全範圍的最大和最小值定位，此時，肯定就定位到了離輸入最小值有一定距離的另外一個最小值，而輸入最小值附近的那些次小值就不會幹擾結果。
注意上面代碼是最小值，因為他用的檢測名額是CV_TM_SQDIFF_NORMED，而非NCC，對于NCC，則需要歸為最大值。
這裡的templatW / 3和templatW / 3有點重疊範圍的意思了，但是還完全不夠。
在最頂層金字塔中找到了多個目标的粗糙位置後，就可以和前所述的一樣的方式一步一步的向下一層金字塔進行細化，直到處理到頂層金字塔為止。
4、多目标+多角度
當問題來到這裡時，整個的難度就是階躍式的提高了一個檔次。
如果目标存在旋轉，為了能找到發生旋轉的物體，我們可以建立多個方向的旋轉對象，也就是說，将搜尋空間離散化，此時，有兩個可選的方式：一個是旋轉搜尋圖像，然後用模闆在旋轉後的圖像中搜尋，二是旋轉模闆，用旋轉後的模闆在搜尋圖像中定位。我們說，第一種方式基本不可取，原因有三。
(1)、搜尋圖像一般來說都是較大的圖，對其進行旋轉耗時比較可觀。、
(2)、實際情況需要多個角度的旋轉，對原圖旋轉記憶體方面也會有過多的消耗
(3)、工業應用時，一般模闆比較固定，而搜尋圖像總是時刻變化的。
當選擇第二種方法時，對于較小的模闆圖像，是可以在執行搜尋前把相關旋轉資訊提前準備好，在搜尋時刻直接使用，而無需做無謂的耗時。
此時，在金字塔的最頂層，需要做的計算工作也有所增加，我們需要對每個角度的模闆都做一個全圖的比對，得到比對的結果，然後對每個可能點，選擇比對度最大的那個角度作為頂層的候選點。
類似的，在向下一層金字塔映射時，不僅僅需要映射比對點的X和Y坐标，還需要映射角度資訊，同理，為了保證角度方面的精度，也需要适當的擴大角度的搜尋區間。
5、更多的問題
實際上，為了實作多目标和多角度的比對，還存在很多很多的細節問題，需要取研究，這些方面的細節有些已經有了部分結論，但是大部分在網絡上鮮有探讨和方向，這裡列出一些問題和大家共同探讨一番。
問題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指令進行處理，尤其是全圖比對的計算過程。
6、公關成果
大概前前後後再一起有折騰了大概一年的時間吧，初步還是有了一定的成果，無論是在速度上，還是在準确度中，還是能解決一定的工程問題了。
共享一些測試圖，和大家一起比較。
模闆圖 待搜尋圖 搜尋結果
模闆圖 待搜尋圖 搜尋結果
基于NCC的速度方面和很多參數有關，比如角度的搜尋範圍，金字塔層數，模闆的大小（一般模闆大，速度反而快些，特别小的模闆則非常耗時，知道原因嗎？），重疊的大小等等都有關。一些簡單的測試資料如下：
雖然曆經千辛萬苦，在磨砺了很久之後，也對這個初有小成，基本實作了這樣的一些過程，但是和halcon相比，無論是從穩定性還是效率方面都還是有一定的差距的，是以标題中的無限接近 就是一句诳語而已。
如果你希望有一個簡單的可視化測試界面，可以從如下連結中擷取，但是請注意這個Demo本身是有一些BUG的（不影響測試使用），請不要将其直接應用到工業環境中，以免造成不必要的損失。
可視化測試Demo: https://files.cnblogs.com/files/Imageshop/TemplateMatching.rar
7、更多疑惑
弄得越多，發現不了解的也越多，特别是在研究比對Halcon的過程中，還發現了一些暫時沒有弄清楚的事情，比如：
①、在不勾選亞像素時，Halcon傳回的坐标值很多情況下也是帶小數點（特别是非0度時的結果），這個作何解釋。
②、當模闆的理論角度是0度時，即使按照AngleStart和AngleStep的值依次計算，取樣的角度值不會準确的取為0，而且也未勾選亞像素，他也能正确的傳回0值，難道他對0度做了特别的處理嗎？
這些還需要慢慢的取探索。
本文的姊妹篇： 工程應用二】 多目标多角度的快速模闆比對算法（基于邊緣梯度)
如果想時刻關注本人的最新文章，也可關注公衆号：