書接上回，經過自動檢測和裁剪之後，已經有很多切割整齊的硬幣照片了，再來看看相似檢測的方法。
一開始的思路極其簡單，既然前面做了邊緣提取，而對於硬幣的種類，人眼最關心的也是和輪廓、花紋相關的資訊。那就利用邊緣資訊來檢測，怎麼檢測呢？
拿待檢測的圖的邊緣，和所有標準圖的邊緣一個一個疊起來比較，和誰重合的面積越大，就和誰最像。雖然這個思路簡單到小學二年級的同學都會笑出聲，但這個方法在使用得當的時候，準確率可以高達40%（不行還是太弱了）
設待檢測的影像的邊緣資訊為
，n個標準影像的邊緣資訊為
，這裡要求
和
影像尺寸均一致。這樣，與待測影像最接近的標準影像為。
說人話：就是看看標準圖和測試圖重合區域的面積，哪個面積越接近測試影像的面積，就是最相像的標準影像（為了防止混亂，這裡強調一下，“面積”指的是圖片裡白色區域的面積，或者說白色畫素塊的數量。這裡涉及的圖都是隻有黑色和白色兩個顏色的）。
老潘用手工方式標註了200個硬幣配對，在這上面做點簡單測試，這種方法準確率只有8%，經過觀察大部分硬幣圖片都被識別成了這兩個
這時候我有一個猜測，如果某個標準影像非常極端的全部為白色，那麼待測影像和標準影像的重合範圍肯定是非常大的（重合範圍就是100%的待測影像），所以上面的方法需要改一下。必須同時考慮測試圖和標準圖。這裡將
的判別方法改寫為
說人話，就是看看標準圖測試圖重合面積、測試圖面積以及標準圖的面積情況，只有同時和兩者都接近的才能被選作相似的圖片配對。
這個方法的準確率迅猛提升到了40%，但是該方法有一個致命的bug，就是一旦硬幣的角度稍有旋轉（實際這個情況卻很容易發生），識別率顯然會受到影響。而且資料量增大也會使準確率和速度逐漸降低。速度降低很顯然，至於準確率下降，老潘在測試的時候確實發現有這個現象，如果想想為什麼這樣，個人猜測一個可能的解釋應該是：資料分佈的空間有限（和圖片的尺寸有關，固定大小的圖片，畫素點取值是有限種可能的情況，而且我們用邊緣資訊，圖中畫素點只取0-1）。資料少的時候，資料分佈更可能散的比較遠，中間的空隙會很大（相當於畫一個圈，大一些小一些歪一些都能把目標點套住），而資料越來越多就會讓各個資料更擁擠，每個資料點之間的距離縮小了，這個時候隨便畫個圈就不是那麼容易套住想要的東西了，很有可能套到別的點上。而這個小學二年級方法顯然不能如此精準地在擁擠的點中選擇目標，所以準確率會下降。（如果有大佬能給出規範的解釋老潘給你投幣三連）
所以有沒有哪種方法可以在旋轉或者光照發生變化時，對同一個硬幣圖案的描述比較穩定，對不同硬幣的描述差別較大呢？顯然這馬上使我們想到了使用LBP方法來嘗試。LBP（Local Binary Patterns）可以看做一類方法的統稱。大致思想就是在圖中選一個區域，把這個區域的變化情況用數字分類表達。這裡老潘找到一個專門針對硬幣分類的LBP思想方法的論文《Image-based coin recognitionusing rotation-invariant region binary patterns based on gradient magnitudes》[1]。論文思路十分簡單，這裡進行簡單的分析與復現。
演算法本身並不困難，總體思路如圖
首先把硬幣圖片切割好（這我們早已做好了），然後求硬幣的梯度圖Gradient Magnitudes（公式2），之後用蜘蛛網狀分割槽，分幾個環和幾個扇（取值多少屬於調參時要考慮的，呼喚調參俠），對每個小塊裡的值，都有特別的計算方法（一個是求區域裡的均值，一個是公式4），再變換到intra RBP和inter RBP（公式6和7）。看每一個環，他的各個扇區一圈就代表了一個二值模式，上面的圓有3圈，所以有3個二值模式，下面的圓同理。
對於每個二值模式編碼編成多少？在同一個作者的《Rotation andflipping robust region binary patterns for video copy detection》[2]裡介紹了。對於所有可以通過翻轉和旋轉重合的模式，就算為一組，如圖
比如說RFR是②的那一排，一圈上只有一個黑點的，都是一樣的模式（可以旋轉和翻轉互相重合）所以作為一組。這裡面的組順序沒什麼特別要求，只要相同模式歸為一組並且編個號就行。
對於評判兩類模式之間的距離，作者用的方法是，計算兩個組之間每組所有模式兩兩的漢明距離（Hamming Distance），取最小的為兩類之間距離。
思路已經捋清楚了，我們甚至可以當堂完成這個作業。這一次老潘用了MATLAB，沒有用python（不知道為什麼涉及到這種純影像的問題我更習慣先開啟MATLAB）。寫程式很快，一晚上就可以寫完甚至不需要熬夜，但是優化這個東西卻耗了老潘一下午。
先來簡單的。第一步求影像的梯度影像，太簡單了不解釋。（我突然覺得這裡也可以優化的）（我寫推送的時候意識到好像可以用sobel運算元）
M = zeros(I_size - 2);
for x = 2:I_size - 1
for y = 2:I_size - 1
M(x-1, y-1) = sqrt((I(x+1, y) - I(x-1, y))^2 + (I(x, y+1) - I(x,y-1))^2);
end
end
然後做影像的扇區，求出每個扇區下的平均值。這裡稍稍有點繞，我選擇先計算環，再計算角度，再取交集，就像下圖所示。
r = (I_size - 2) / 2;
center_point = r;
axis_x = zeros(I_size - 2);
for row = 1:I_size - 2
axis_x(row,:) = [1:I_size - 2];
end
axis_y = zeros(I_size - 2);
for column = 1:I_size - 2
axis_y(:,column) = [1:I_size - 2];
end
condition_r = abs(sqrt((axis_x-center_point).^2+ (axis_y-center_point).^2));
condition_angle =atan2(axis_y-center_point, axis_x-center_point)/pi*180+180;
pic_num = 1;
m = zeros(rings, sub_regions);
for ring = 1:rings
for region = 1:sub_regions
area_r = condition_r < ring*(r/rings) & condition_r >=(ring-1)*(r/rings);
area_angle = condition_angle < region*(360/sub_regions) &condition_angle >= (region-1)*(360/sub_regions);
area_need = area_r .* area_angle;
area_gradient = M .* area_need;
插個嘴做一些解釋：這裡的area_r和area_angle相當於蒙版，我們只要直接和圖片點乘就可以選出要的區域（1乘一個數還是這個數，0乘一個數一直是0）
取兩個蒙版的交集也是同理，直接互相點乘即可，得到的矩陣就是交集area_need。
註釋掉的地方可以生成一個動態的GIF圖片，展示了每圈每個扇區下的梯度影像情況。
之後計算每個區域下的均值，存下來即可。
% % save a gif to observe
% imshow(area_gradient, [], 'border','tight', 'initialmagnification','fit');
% F=getframe(gcf);
% I=frame2im(F);
% [I,map]=rgb2ind(I,256);
% if pic_num == 1
% imwrite(I,map,'RFR_log/area.gif','gif','Loopcount',inf,'DelayTime',0.1);
% else
% imwrite(I,map,'RFR_log/area.gif','gif','WriteMode','append','DelayTime',0.1);
% end
% pic_num = pic_num + 1;
m(ring, region) = sum(area_gradient(:)) / sum(area_need(:));
end
end
之後的公式3至7的計算略，沒什麼好解釋的。貼出來也只會令人頭暈。
好了，就差構造一個查詢表IND和一個距離矩陣了，但是具體做的時候就沒那麼簡單了。
1
先說構造IND，這個表用來把一個010101這樣的模式換成他所在類的號碼。因為組比較多，所以用一次算一次比較慢，不如用空間換時間。於是第一天晚上掛機一宿，發現這個表還沒羅列好。遂優化。一開始的想法是用一個表格當做查詢表，行數就是類別號，每行儲存這類的二值模式0101碼。但是問題一，表裡會有很多冗餘資料，比如全0的這類，由於陣列尺寸必須一樣，需要把後面空缺的都填滿。就存了20個全0的模式。問題二，時間很慢而且會越來越慢。因為每一個模式都要去表裡查詢是否出現過，表越大查的越慢。
於是想了一招，採用Map格式儲存，首先要宣告一下變數
IND = containers.Map;
代表我們要用Map儲存IND表格，之後，為了提高查詢效率，不用類號做索引，用二值模式做索引，這樣一來可以很快查到表裡有沒有出現過這個二值模式，二來在使用的時候，把圖的二值模式轉成類別號也會快。這時候遍歷所有的二值模式，對於這一組二值模式，如果有一個沒在之前出現過，那就把這一組全存下，並且給一個統一的組號（這一組是用某一個二值模式經過迴圈移位和翻轉出來的所有結果）。如果在之前出現過，就代表這一組的二值模式都已經在這個表裡有過了（通過旋轉和移位終究會變成一樣的）。判斷和儲存的方法為
if ~isKey(IND, dec2bin(CP_FP_dec(1),sub_regions))
for index_map = 1:length
IND(dec2bin(CP_FP_dec(index_map), sub_regions)) = ind;
end
ind = ind + 1;
end
（只記錄了判斷鍵值是否出現過，整體程式碼見“閱讀原文”）
於是一宿未完的問題，10分鐘就可以解決了。其中5分鐘用來把資料儲存到磁碟上用來下次直接查詢該換固態硬碟了。
2
下一個問題是計算兩類之間的距離，原論文仍然是做了一個矩陣，可以直接在裡面查兩類的距離，但是我做的時候發現這個矩陣太大，記憶體會爆（計算一下這個矩陣只佔729Mb，MATLAB居然放不下）。我們也不知道作者用了什麼黑科技（也許是轉了C語言），我想了想，乾脆在判斷圖片的時候直接計算，也差不了多長時間。
還是挺佔時間的，實際用的時候發現只判斷兩個圖片就要將近半分鐘。這要是再一個一個比對過去，天荒地老。遂優化。最開始就是粗暴的兩個for迴圈，然後一個一個的取兩個類別號，然後計算。發現不行，有點慢，於是改用矩陣計算。先做兩個類別號（其實是兩組，一個圖片會有一組類別號）的排列組合矩陣
% reference:https://www.ilovematlab.cn/thread-321138-1-1.html
[m, n] = meshgrid(RFR_i_list_dec,RFR_j_list_dec');
[res(:,1), res(:,2)] = deal(reshape(m, [],1), reshape(n, [], 1));
然後計算這個矩陣按行的漢明距離，取這個矩陣的min就是兩個類的距離。
DH = dec2bin(bitxor(res(:,1), res(:,2)));
[length, ~] = size(DH);
DH2 = uint8(zeros(length, 1));
for i = 1:length
DH2(i) = sum(DH(i,:)=='1');
end
min_Hamming = min(DH2);
之後發現還是很佔時間，偶然瞥到，在每一輪迴圈裡都有這麼個操作
key = keys(IND);
value = uint32(cell2mat(values(IND)));
把類號和二值模式從Map格式裡面拿出來，這個拿的動作很耗時。所以就一開始取出，直接作為引數傳給後面的一系列函式。原本兩個圖片對比時間要大概半分鐘，現在只需半秒。速度提升了60倍。
不過速度還是比較慢，不知道那個線上識別硬幣的網站用的到底是什麼方法。如果選擇做出來距離對照表的話，拿空間換時間，應該速度就會很快了。但是我估算目前的小需求下，做查詢表的時間和直接做計算的時間，好像做表不是太划算。
話說回來，用這個方法的檢測準確率可以達到45.5%，遠不及論文中所達到的準確率（89.1%）。我使用了和作者一致的引數，但是資料集不一樣。作者使用的是遠古版本的MUSCLE-CIS資料集，這是在2006年辦的一個硬幣識別比賽用的資料集，原連結已經失效，在網上也很難搜尋到相關的資料，只有在偶爾一些論文裡有零星展示，在論文[3]裡給了兩個資料集內的圖片。
可以看到這個影像和直接用手機拍的畫面差別挺大的，這個花紋的邊緣非常清晰，有可能是用了特殊的光源和相機拍攝的。在這種優良乾淨的資料上，我相信小學生方法也能有不錯的效果。但是在具體生活場景上，如何在光照不佳，資料量小的情況下提高檢測效率，還是一個有待研究的問題。
經過簡單的進一步調研，在文章[4]中提到，MUSCLE-CIS解析度低並且經過了很強的預處理，實際場景很難得到這樣的圖片。但是硬幣資料集還有很多，例如文章中給出的這個例子，可以看到實際生活中主要面對的還是後面這幾類影像。這個文章著眼於提取硬幣表面的文字，更深刻的分析了硬幣識別裡面的各種技術問題，而且作者他們也自己建立了一個資料集。論文線上網址：
https://www.frontiersin.org/articles/10.3389/fict.2017.00009/full
（這個作者也研究過硬幣分類的問題，用的貌似是類似SIFT那種方法，看起來好像很吊的樣子，這篇文章[5]的方法作為作業下節課找同學上臺講解）
這次老潘用的自建資料集也會放在GitHub上，想研究的同學可以點選閱讀原文拿去玩一玩。
參考文獻：
[1] Kim S, Lee S H, Ro Y M, et al.Image-based coin recognition using rotation-invariant region binary patternsbased on gradient magnitudes[J]. Journal of Visual Communication and Image Representation, 2015: 217-223.
微信掃一掃
關注老潘的公眾號，第一時間獲取新文章