[快速閱讀八] Matlab中bwlookup的實現及其在計算二值影像的尤拉數、面積及其他morph變形中的應用。

　　以前看過matlab的bwlookup函式，但是總感覺有點神秘，一直沒有去仔細分析，最近在分析計算二值影像的尤拉數時，發現自己寫的程式碼和matlab的總是對不少，於是又去翻了下matlab的原始碼，看到了matlab裡實現尤拉數的程式碼非常簡單，如下所示：

if n==4
    lut = 4*[0 0.25 0.25 0 0.25 0  .5 -0.25 0.25  0.5  0 -0.25 0 ...
             -0.25 -0.25 0] + 2;
else
    lut = 4*[0 0.25 0.25 0 0.25 0 -.5 -0.25 0.25 -0.5  0 -0.25 0 ...
             -0.25 -0.25 0] + 2;
end
% Need to zero-pad the input
b = padarray(a,[1 1],'both');

weights = bwlookup(b,lut);
if coder.isColumnMajor
    e = (sum(weights(:),'double') - 2*numel(b)) / 4;
else
    e = (sum(sum(weights,2,'double'),1,'double') - 2*numel(b)) / 4;
end

　　在內部就是呼叫了bwlookup函式，那沒辦法了，就仔細看了M的幫助文件，發現原來這個函式真的非常簡單，我們貼下M的幫助文件：

　　A = bwlookup(BW,lut)

　　The bwlookup function performs these steps to determine the value of each pixel in the processed image A:

• Locate the pixel neighborhood in input image BW based on the coordinates of the target pixel in A. The function zero-pads border pixels of image BW when the neighborhood extends past the edge of BW.

• Calculate an index, idx, based on the binary pixel pattern of the neighborhood.

• Set the target pixel in A as the value of the lookup table at the index idx, in other words, the value of lut(idx).

2-by-2 Neighborhood Lookup

For 2-by-2 neighborhoods, there are four pixels in each neighborhood. Each binary pixel has two possible states, therefore the total number of permutations is 24 and the length of the lookup table lut is 16.

To find the value of the target output pixel at (row, column) coordinate (r,c), bwlookup uses the 2-by-2 pixel neighborhood in the input binary image BW whose top left pixel is at coordinate (r,c). The index idx into the lookup table is the weighted sum of the four pixels in the neighborhood, plus 1.

3-by-3 Neighborhood Lookup

For 3-by-3 neighborhoods, there are nine pixels in each neighborhood. Each binary pixel has two possible states, therefore the total number of permutations is 29 and the length of the lookup table lut is 256.

To find the value of the target output pixel at (row, column) coordinate (r,c), bwlookup uses the 3-by-3 pixel neighborhood in the input binary image BW centered on coordinate (r,c). The index idx into the lookup table is the weighted sum of the nine pixels in the neighborhood, plus 1.

　　簡單的說，這個查詢表只有2種可能，一個是表的元素是16個，一種是由512個元素，其中16個元素時，查表的依據是以當前點為起點，向右向下各擴充套件一個畫素範圍，得到一個2*2領域4個畫素，4個畫素，每個畫素也就只有2種可能的取值，這樣就有16種可能的組合，對應16個元素。

　　當表的元素是512個時，查表的依據是以當前點為中心點，向左向右，向上向下各擴充套件一個畫素，得到一個3*3的領域9個畫素，同樣每個畫素也只有2種可能取值，所以共有1<<9即512種取值。

　　因為涉及到領域，所有在原圖邊緣的地方會有越界的影像，這個時候的原則是填充0，而不是複製最近的畫素。

　　這樣，當表中有不同的內容時，就可以實現不同的結果。

　　我們以16個元素的表為例，假定 lut[16] ={ 6 3 16 11 7 14 8 5 15 1 2 4 13 9 10 12},

　　對於下面的一個4*4的二值資料：

   0   0   1   1
   0   0   1   1
   1   1   0   0
   1   1   0   0

　　如果當前的點位於第二行第一列，則其2*2的領域範圍的畫素資訊為：

   0   0
   1   1
　　按照位置資訊即對應的值計算索引為:  1 + 0 * 1 + 1 * 2 + 0 * 4 + 1 * 8 = 11， 對用的查表結果為 lut[11] = 2。
　　使用類似的計算方法可以得到其他位置經過查表後的結果為：

  　　6    13    12    11
     2     8    14     3
    12    11     6     6
    14     3     6     6

　注意，以上所有的下標起點都是1(matlab內部的約定），而且行列順序和我們常用的C語言的二維陣列也是掉了邊的。 
  除了剛剛說的尤拉數的計算（其實沒有搞懂尤拉數里的那個查詢表是如何設計的，且8領域的尤拉數也是用的2*2的表，而不是3*3的），我們以matlab的bwarea函式為例，說明這個函式的價值。

  matlab中關於bwarea函式的定義為：

　　bwarea 透過對影像中每個畫素的面積求和來估計影像中所有 on 畫素的面積。單個畫素的面積是透過觀察其 2×2 鄰域來確定的。有六種不同的情形，每種情形表示一個不同面積：

• 具有零個 on 畫素的情形（面積 = 0）

• 具有一個 on 畫素的情形（面積 = 1/4）

• 具有兩個相鄰 on 畫素的情形（面積 = 1/2）

• 具有兩個對角 on 畫素的情形（面積 = 3/4）

• 具有三個 on 畫素的情形（面積 = 7/8）

• 具有所有四個 on 畫素的情形（面積 = 1）

　　用圖形化的表達方式上述六種情況對應如下（為了顯示，使用綠色表示黑色值，紅色表示白色值）：

　 　　　　　　 　　　　　　　　　　　　　　　　　　　　

索引值 　 0　　 　 1　　 　 2　　　　3　　 4　　　 5　　　 6　　　 7　　　 8　　 9　　 10　　 11　　 12　　 13　　 14　　 15　　

面積值　　　　 0　　 　1/4　　 1/4　　　1/2　　 1/4　　 1/2　　 3/4　　　 7/8　　　1/4　　 3/4　　 1/2　　 7/8　　 1/2　　 7/8　　 7/8　　 1

把面積值放大八倍得到面積值依次為： 0 2 2 4 2 4 6 7 2 6 4 7 4 7 7 8，我們翻看matlab的bwarea函式，可以得到其程式碼如下：

lut = [0     2     2     4     2     4     6     7     2     6     4     7 ...
       4     7     7     8];

% Need to zero-pad the input.
bb = false(size(b)+2);
bb(2:end-1,2:end-1) = b;

weights = bwlookup(bb,lut);
a = sum(weights(:),[],'double')/8;

　　查詢表和這裡的一摸一樣。

　　使用查詢表的優點是把領域所有的組合提前計算出結果來，而不用到程式裡進行一些列複雜的判斷，這樣在很多情況下是可以提高速度的，比如剛剛這個面積計算的東西，如果在迴圈內部做，則要做N種判斷，那如果是涉及到3*3的領域，那就更為複雜了。

　　樸素的來說，有了bwlookup的原理，要把這個演算法移植到C++中，就是一個比較簡單的過程了，假如是處理16元素的表，也就是涉及到了4個畫素，假定他們的值分別為P0\P1\P2\P3，則最基本的C程式碼如下所示：

int Index = 0;
if (P0 == 255)    Index += 1;
if (P1 == 255)    Index += 2;
if (P2 == 255)    Index += 4;
if (P3 == 255)    Index += 8;
LinePD[X] = Lut[Index];

　　處理512個元素的過程也是類似的，只不過索引值多加了幾個（注意，這裡用255表示白色，而不是1）。

　　如果考慮指令集最佳化，一個自然的翻譯過程如下所示：

    __m128i Index = _mm_setzero_si128();
    Index = _mm_blendv_epi8(Index, _mm_add_epi8(Index, _mm_set1_epi8(1)), _mm_cmpeq_epi8(P0, _mm_set1_epi8(255)));
    Index = _mm_blendv_epi8(Index, _mm_add_epi8(Index, _mm_set1_epi8(2)), _mm_cmpeq_epi8(P1, _mm_set1_epi8(255)));
    Index = _mm_blendv_epi8(Index, _mm_add_epi8(Index, _mm_set1_epi8(4)), _mm_cmpeq_epi8(P2, _mm_set1_epi8(255)));
    Index = _mm_blendv_epi8(Index, _mm_add_epi8(Index, _mm_set1_epi8(8)), _mm_cmpeq_epi8(P3, _mm_set1_epi8(255)));
    _mm_storeu_si128((__m128i*)(LinePD + X), _mm_shuffle_epi8(Lut128, Index));

　　一次性處理16個位元組，最為關鍵的一點是，這個查表可以方便的直接使用_mm_shuffle_epi8非常高效的實現，這個在我前期的文章裡多次提到（16個位元組表的查詢，我們假定lut的型別為位元組型別，這個在實際的應用中也足夠了），

　　這個程式碼的效能提升相對於C++來說是非常可觀的，大概又4倍多的速度提升。

　　仔細上述程式碼其實還有提高的地方，我們使用_mm_blendv_epi8來進行結果的混合， 而混合的標誌是值是否等於255，這裡用_mm_cmpeq_epi8做判斷，而_mm_cmpeq_epi8的結果就是如果P等於255，則返回255，否則返回0，那這個不就是P本省嗎，所以根本就不用做這個判斷，這樣程式碼就變為：

    __m128i Index = _mm_setzero_si128();
    Index = _mm_blendv_epi8(Index, _mm_add_epi8(Index, _mm_set1_epi8(1)), P0, _mm_set1_epi8(255));
    Index = _mm_blendv_epi8(Index, _mm_add_epi8(Index, _mm_set1_epi8(2)), P1, _mm_set1_epi8(255));
    Index = _mm_blendv_epi8(Index, _mm_add_epi8(Index, _mm_set1_epi8(4)), P2, _mm_set1_epi8(255));
    Index = _mm_blendv_epi8(Index, _mm_add_epi8(Index, _mm_set1_epi8(8)), P3, _mm_set1_epi8(255));
    _mm_storeu_si128((__m128i*)(LinePD + X), _mm_shuffle_epi8(Lut128, Index));

　　另外，再觀察這個C語言的特殊性，我們可以把上述C程式碼修改為：

    int Index = 0;
    Index += (1 & P0);
    Index += (2 & P1);
    Index += (4 & P2);
    Index += (8 & P3);
    LinePD[X] = Lut[Index];

　　即直接使用位運算替換掉那個判斷，這樣對應的SSE指令也可以進行修改為如下程式碼：

    __m128i Index = _mm_setzero_si128();
    Index = _mm_add_epi8(Index, _mm_and_si128(_mm_set1_epi8(1), P0));
    Index = _mm_add_epi8(Index, _mm_and_si128(_mm_set1_epi8(2), P1));
    Index = _mm_add_epi8(Index, _mm_and_si128(_mm_set1_epi8(4), P2));
    Index = _mm_add_epi8(Index, _mm_and_si128(_mm_set1_epi8(8), P3));
    //    可以直接接用shuffle實現查詢表
    _mm_storeu_si128((__m128i*)(LinePD + X), _mm_shuffle_epi8(Lut128, Index));

　　相比於原始的SSE程式碼，這個改動也約有30%的效能提升的。

　　對於512個元素的表，情況會有所不同，因為索引值大於了255，所以已經無法用位元組型別來儲存累加後得到的索引了，至少的用ushort型別，但是呢，前8個位置的索引相加還是不會超出位元組型別的，所以前8個位置還是可以一次性處理16個畫素，到最後一個位置時，單獨轉換為16位，然後再接用2次16資料的處理指令，就可以一次性得到16個位元組對應的查詢表索引。

　　注意，512個元素的查表如果是純C++程式碼， 最後一個的 & 操作，一定要注意資料範圍，不能寫成 Index += (256 & P9);, 而是 Index += (256 & (65280 + P9));

　　但是這個時候，由於表的大小時512，所以已經無法使用SSE去最佳化這個查表功能了，只能把索引值單個的提取出來，然後用普通的C語言去執行程式碼。如果CPU能支援AVX2，那麼AVX2倒是又有關指令能直接查表，不過也要做很多的改造工程。

　　我們測試，對於512個元素的表，最佳化後的SSE指令處理一副 3000*2000的二值圖，大概耗時在2ms不到，速度還是相當的快的。

　　搜尋matlab的程式碼，除了bweuler及bwarea內使用了bwlookup，另外就是bwmorph裡也大量的使用了bwlookup，而且都是使用的512個元素的表，也就是說使用3*3的領域，我們看bwmorph的幫助文件，有很多相關內容：

　　這些對應的查詢表在 MATLAB\R2023b\toolbox\images\images\+images\+internal 目錄下以lut開頭的檔案中可以找到，比如說，這個clean的表內容如下：

　　

　　這些表的內容都是提前設計好的，然後可以多次重複的呼叫，從而得到某種效果。

　　一般來說，這種3*3的領域運算元，在迭代了一定的次數後效果就不會有變化了。

　　 我們在morph中也能看到一些常用的運算元，比如這個remove就可以得到二值影像的邊界，這個majority可以平滑噪音（和我部落格裡專門講的那個majority效果還是有差異的）。

　　　　　　　　原圖　　　　　　　　　　　　　　　　　　remove　　　　　　　　　　　　　　majority

　　　　 fatten 　　　　　　　　　　　　　　　　　 skeleton　　　　　　　　　　　　　　 　　thin

　　 個人感覺這個方法在處理一些比較複雜的領域資訊時還是很有幫助的，特別是不方便用判斷條件一個一個的寫的，如果能提前把這個表弄出來，那就效率能得到很大的提升的。

　　 在個人的SSE Demo裡，也整合了這個演算法，其內容在Binary(二值處理）--》Processing(後處理)-->> Morph(形態學)中，有興趣的朋友可以看看。

　　 SSE Demo下載地址： https://files.cnblogs.com/files/Imageshop/SSE_Optimization_Demo.rar

翻譯

搜尋

複製