繼去年上半年一鼓作氣研究了幾種不同的模版匹配演算法後，這個方面的工作基本停滯了有七八個月沒有去碰了，因為感覺已經遇到了瓶頸，無論是速度還是效率方面，以當時的理解感覺都到了頂了。年初，公司業務慘淡，也無心向佛，總要找點事情做一做，充實下自己，這裡選擇了前期一直想繼續研究的基於離散夾角餘弦相似度指標的形狀匹配最佳化。
在前序的一些列文章裡，我們也描述了我從linemod模型裡抽取的一種相似度指標用於形狀匹配，個人取名為離散夾角餘弦，其核心是將傳統的基於梯度點積相似度的的指標進行了離散化：
傳統的梯度點積計算公式如下：
對於任意的兩個點，透過各自的梯度方向，按照上述公式可計算出他們的相似度。
那麼離散夾角餘弦的區別是，不直接計算兩個梯度方向的餘弦，而是提前進行一些定點化，比如，早期的linemod就是把360度角分為了8份，每份45度，這樣[0,45]的梯度方向標記為0，[45,90]的梯度方向標記為1，依次類推，直到[315, 360]之間的梯度方向標記為7，這樣共有8個標記，然後提前構建好8個標記之間的一個得分表，比如，下面這樣的一個表：
這個表的意思也很簡單，就是描述不同標記之間的得分，比如兩個點的梯度方向，都為3或者4或者5，那麼他們的得分就比較高，可以取8，如果一個為1，一個為6，則得分就只有2， 一個為0，一個為5，則得分為0（即完全相反的兩個方向）。
很明顯，角度量化的越細，則得到的結果越和傳統的梯度點積結果越接近。但是計算量可能也就越大，
關於這個過程，我去年的版本也有弄過8角度、16角度及32角度，個人覺得，在目前的計算機框架下，16角度應該是既能滿足進度要求，又能在速度方面更為完美的一個選擇。
這裡記錄下最近對基於16角度離散餘弦夾角指標的形狀匹配的進一步最佳化過程。
一、核心的最佳化策略
透過前面的描述，我們知道，這種方法的得分是透過查表獲取的，而且，在大部分的計算中，是沒有涉及到浮點計算的，我們透過適當的構造表的內容，可以透過簡單的整數型別的加減乘除來得到最後的得分。這個的好處有很多，其中一個就是精度問題，在基於梯度點積的計算中，如果採用float型別來累計計算結果，通常或多或少的存在某些情況下的精度丟失，而且還不好定位哪裡有問題。還有個好處就是真的可以加速，當然這個的加速主要是透過一個很特殊，但是又很有效的SSE指令集語句實現的，這個指令就是_mm_shuffle_epi8，其原型及相關介紹如下：
__m128i _mm_shuffle_epi8 (__m128i a, __m128i b)
這是個很牛逼的東西，如果我們把a看成一個16個元素的位元組陣列，b也是一個16個元素的位元組資料，則簡單的理解他就是實現下述功能：
dst[i] = a[b[i]]; 即一個16個元素的查表功能。
很明顯，b[i]要在0到15之間才有效，否則，就查不到元素了，但是該指令還有比較隱藏的功能是，當b[i]大於15時，dst就返回0了。
為什麼說這個指令牛逼，我們看我們前面說的這個獲取離散夾角餘弦的過程，對於兩個等面積的區域，假定一個區域的量化後的離散角度標記儲存在QuantizedAngleT記憶體中，另外一個儲存在QuantizedAngle記憶體中，他們的寬和高分別為ValidW和ValidH，則這兩個區域按照前面定義的標準，其得分可用下述程式碼表示（這裡SimilarityLut是16角度離散的）:
int Score = 0; for (int Y = 0; Y < ValidH; Y++) { int Index = Y * ValidW; for (int X = 0; X < ValidW; X++) { Score += SimilarityLut[QuantizedAngleT[Index + X] * 16 + QuantizedAngle[Index + X]]; // Score += SimilarityLut[QuantizedAngleT[Index + X], QuantizedAngle[Index + X]]; } }
如果把SimilarityLut看成是二維的陣列，上面註釋掉的得分可能看起來更為清晰。
實際上，我們在進行模版匹配的時候大部分都是在進行這樣的得分計算，因此，如果上面的過程能夠提速，那麼整體將提速很多。
通常，查表的演算法是無法進行指令集最佳化的（AVX2的gather雖然有一定效果，但是弄的不好會適得其反），但是，正是因為我們本例的特殊性，使得這個查表反而更有利於演算法的效能提高。
注意到前面有說過我們在進行16角度量化時，量化的標記範圍是[0,15]，意味著什麼，這個正好是_mm_shuffle_epi8 指令裡引數b的有效範圍啊。
但是仔細看上面的SimilarityLut表，他由兩個變數確定索引，這就有點麻煩了，解決的辦法是換位思考，我們能不能固定其中的一個呢，這個就要結合我們的實際應用了。
在形狀匹配中，我們提取了很多特徵點，然後需要使用這些特徵點對影像中有效區域範圍的目標進行得分統計，也就是說任何一個位置，都要計算所有特徵點的得分，並計總和，一個簡單的表示為：
for (int Y = 0; Y < ValidH; Y++) { int Index = Y * ValidW; for (int X = 0; X < ValidW; X++) { int b = 影像位置對應的量化值 for (int Z = 0; Z < Template.PointAmount; Z++) { int a = 模版位置對應的量化值 Score += SimilarityLut16[a, b]; } } }
這種情況a，b在每次獨立的小迴圈中還是變化的，一樣無法使用指令集。
不過，如果我們調換下迴圈的順序，改為以下方式：
for (int Z = 0; Z < Template.PointAmount; Z++) { int a = 模版位置對應的量化值 for (int Y = 0; Y < ValidH; Y++) { int Index = Y * ValidW; for (int X = 0; X < ValidW; X++) { int b = 影像位置對應的量化值 Score += SimilarityLut16[a, b]; } } }
此時，a在內部迴圈裡是不變的，我們透過a之可以定位到SimilarityLut16的a*16地址處，並載入16位元組內容，作為查詢表的內容，而b值也可以一次性載入16個位元組，作為查詢的索引，這樣一次性就能獲得16個位置的得分了。
還有一點，我們在演算法裡有個最小對比度的東西， 這個東西是用來加快速度的，即梯度值小於這個資料，我們不要這個點參與到匹配中，即此時這個點的得分是0，為了標記這樣的點，我們需要再原圖的量化值裡增加不在[0,15]範圍內的東西，一旦有這個東西存在，我們的普通C程式碼裡就需要新增類似這樣的程式碼了：
if (QuantizedAngle[Index + X] != 255) Score += SimilarityLut[QuantizedAngleT[Index + X] * 16 + QuantizedAngle[Index + X]];
這裡我使用了255這個不在[0,15]範圍內的數字來表示這個點不需要參與匹配計算。本來說，如果剛剛那條_mm_shuffle_epi8指令，只是純粹的實現[0,15]索引之間的查表，那有了這個東西，這個指令又就沒法用了，但是恰恰這個指令能實現在索引不在0到15之間，可以返回其他值，而這個其他值恰好又是0，你們說是多麼的巧合和神奇啊。所以，這一切都是為這個指令準備的。我們看看其他shuffle指令，包括_mm_shuffle_epi32，_mm_shuffle_ps，都沒有這個附帶的功能。
其實這裡還是要交代一點，這個演算法，如果遇到那種不能用指令集的機器，或者說用純C語言去實現，效率就比較低了，因為C語言裡只能直接查表，而且還要有那個判斷。
二、特徵點數量的展開即貪婪度引數的捨棄
linemod裡使用8角度的特徵，其兩個特徵之間的得分最大值為8，其內部使用了16位的加法，所以其最大的特徵點數量為8096個，當模型較大時，往往會超出這個數量的特徵點，特別是在最後面基層金字塔的時候，一種方案就是我們限制特徵點的數量，並對特徵點進行合適的提取，這也不失為一種方案。 那想要完美呢，就必須還得是上32位的加法。這裡就存在一個問題，精度和速度如何同時保證，畢竟在SSE指令集的世界裡，16位的加法是要比32位的加法快的。
一種解決方案就是，對特徵點進行分割槽，我們按照可能超出16位能表達的最大範圍時特徵點的數量為區間大小，對特徵點按順序進行多區間劃分，在每個小區間裡還是用16位的指令集加法，完成一次後，把臨時結果在加到32位的資料裡。這樣就精度和速度都能兼顧了，只不過程式碼又複雜了一些。
當我們嘗試了這麼多努力後，我們發現無論是頂層的得分計算，還是後續的每層的區域性更新，其速度都變得飛快，這個時候我們又想到了一個貪婪度引數，這個東西，在論文有提到，可以提前結束迴圈，加快速度。可是，我也想嘗試把這個東西加入到我這個演算法的過程裡，我發現他會破壞我整體的節奏，最終我們選擇捨棄了這個引數，核心理由如下：
1、提前結束迴圈，是需要進行判斷的，而且是每次都要判斷，特別是對於後期的區域性更新判斷，因為大部分候選點都是能滿足最小得分要求的，這部分的判斷一般來說，基本上就無法滿足了，也就是純粹的多了這些無效判斷。
2、因為我們使用了_mm_shuffle_epi8指令，一次性可以處理16個位置的得分，也就是這個粒度是16個畫素，而如果使用SSE進行判斷，也只有當16個位置都不滿足最小得分要求後，才可以跳出迴圈，這個在很多情況下也是得不償失的。
三、目標只有部分在影像內的識別
有些情況下，被識別的目標只有區域性在影像範圍內，而我們也期望能識別他，這個功能，我知道早期版本的halcon是沒有的，他只能識別那些特徵點完全在影像內的目標（不是模版影像邊緣）。我早期版本也麼有這個功能，後期有做過一些擴充套件，擴充套件的方法是透過擴大原始影像合適的範圍，同時為了避免不增加新的邊緣資訊，擴充套件的部分都是用了邊界的畫素值。這樣做在大部分情況下是能夠解決問題的，不過其實也隱藏的一些不合理的地方，這些擴充套件的部分在細節上還是會產生額外的邊緣的，只是不怎麼明顯。因此，這個版本，我也考慮了幾個最佳化，在內部直接實現了邊緣的擴充套件。
這裡有幾個技巧。
1、原圖的金字塔影像還是不動。
2、計算原圖每層金字塔影像的角度量化值時，對這個量化值進行擴充套件，擴充套件的部分的量化值填充前面說的那個不在[0,15]之間的無效值，比如這裡是255，這樣，這些區域的得分就是0。
3、為了能保證在極端情況下這些部分在影像的目標能被識別，擴充套件的大小要以特徵有效值的外接矩形的對角線長度為基準，再進行適當的擴充套件（這個擴充套件也有特別要求）。
4、計算完成後，座標值要進行相應的調整。
透過這種方式，可在內部實現缺失目標的識別，而且在記憶體佔用、速度等方面也有一定的優勢。
四、最小外接矩形識別重疊
halcon有說過其maxoverlap引數是透過計算特徵點的最小外接矩形之間的重疊來實現的，在我以前的版本中，這個功能是透過其他的簡易方法來搞定的。這個主要是以前搞不定最小外接矩形的計算，年初，恰好從opencv裡翻譯可扣取了一些程式碼，起重工就有最小外接矩形的獲取，這個需要透過計算凸包以前其他一些複雜的東西搞定，我沒有看懂原理，只是直接扣取了程式碼，不過CV的程式碼繞來繞去，扣的我也是頭暈腦脹，總算搞出來了。
那麼這裡其實也有蠻多的細節和可選方案，我列舉如下：
1、在建立特徵時，計算好每個旋轉後的特徵的最小外接矩形（勾選了預生成模型資料）。
2、在最後確定的底層金字塔裡所有的候選點出計算每個特徵點對應的外接矩形。
3、只計算底層金字塔0角度是特徵單的最小外接矩形，然後其他底層金字塔的最小外接矩形用他旋轉得到。
我們實際考慮啊，方案一對建立模型不友好，方案二實際測試對執行的效率產生了不良影響，方案3最好，基本不耗時，而且對精度的影響也非常有限，所以可以選擇方案3。
五、其他的一些我未公開的討論的課題
1、16角度SimilarityLut的值如何設計，其實在halcon裡有個metric引數，他有三種選擇，使用極性、忽略全域性極性、忽略區域性極性。如果想前面給出的那種8角度的SimilarityLut，是隻能實現使用極性和忽略全域性極性的。這裡適當擴充套件，就可以實現三個都有了，而且對是速度提升還有好處。
2、5*5區域性得分過程的特別最佳化，尤其是如何高效的載入每行5個位元組，並拼接成合適的形式，使得能快速的使用指令集。
3、也可以使用8*8的區域性區域（非對稱的區域性更新），這樣方便使用指令集，但是因為數量變大，還是沒有最佳化後的5*5快。
4、在最頂層，計算候選點時，可以直接計算，也可以考慮使用ResponseMaps，其中，測試表明還是使用ResponseMaps快一點。
5、還是候選點的選擇問題，在最頂層，目前我還是用的某個角度下的二維得分結果中選擇得分大於最小得分要求，同時是5*5領域的最大值作為候選點，這種方式留下的候選點還是有很多的，對於只有旋轉的匹配，是否可以考慮在3D（X方向、Y方向以及角度方向）的空間裡，選得分大於最小得分要求且是5*5*5領域的最大值呢，這樣候選點肯定會少很多，但是程式碼的編寫似乎變得困難了很多，還有佔用記憶體問題。那如果擴充套件到多尺度的匹配，或者是各項異性的匹配，那就要在4D和5D空間搜尋最大值了，這個感覺就更為複雜了。
6、另外，再有頂層金字塔向下層金字塔迭代更新的過程中的候選點的捨棄問題，也值得探討，到底是隻根據得分是否滿足需求，還是根據一些物理空間或者角度方面的特性做些特別的最佳化和捨棄呢，而且這種捨棄行為本身有的時候可能會帶來效能的下降，因為在搜尋那些目標可以捨棄時，需要一個迴圈，當候選目標有幾千個的時候，兩重這樣的迴圈會對對後續候選點變少帶來的效能加速起到反向作用，這個甚至會超過候選點變少帶來的加速。所以 目前我這個版本為了穩定性，只是對得到的重複的點做了捨棄。
7、原本再想一個最佳化，即我們的特徵點的儲存順序問題，現在只有0角度的特徵點的儲存順序是X及Y方向都是由小向大方向座標排列，這樣訪問的時候和影像的記憶體佈局方向性一致，按理說cache miss要小一些。但是，如果按照順序把0角度的特徵點旋轉後，得到新的位置資訊，一般來說肯定不是按照這樣的順序排列的了，所以訪問時隨機性就差了一些，那是否我再計算前把這些點再按0角度時那樣做個排序後，有利於演算法的效能呢，我做了實際的測試，應該說基本沒啥區別吧。
六、結論
綜合以上各種最佳化手段後，目前經過測試，速度較以前有很大的提高，而且和基於傳統的梯度點積的方法比較，速度更具有優勢，而且精度也不逞多讓。我們測試在2500*2000的灰度中查詢 300 * 150的多目標，大約耗時11ms, 傳統的方式結合貪婪度耗時也需要18ms。這個時間我測試和halcon的比較已經非常接近了。當然這種比較還是要看具體的測試圖。
從演算法精度上看，怎麼上定位也是很準確的，在執行過程中，佔用的記憶體也不大，因此，個人覺得這個方法不失為一個優質的運算元。
本文測試DEMO連結：https://files.cnblogs.com/files/Imageshop/QuantizedOrientations_Matching.rar?t=1710810282&download=true
如果想時刻關注本人的最新文章，也可關注公眾號或者新增本人微信： laviewpbt
