超越halcon速度的二值圖像的腐蝕和膨脹，實現目前最快的半徑相關類算法（附核心源碼）。

　　我在兩年前的博客裏曾經寫過 SSE圖像算法優化系列七：基於SSE實現的極速的矩形核腐蝕和膨脹（最大值和最小值）算法  一文，經過SSE的優化把矩形核心的腐蝕和膨脹作到了不只和半徑無關，並且速度也至關的快，當時在被博文的評論裏有博友提出了以下的問題：

#1樓
2018-02-21 20:26 | 胡一譚  
博主的思路很巧妙，只是這個算法自己仍是不夠快，優化效果與商業軟件仍是有比較大差距，4096X8192大小的的灰度圖商業軟件（halcon）只須要33ms， 本文須要250ms，考慮到商業軟件採用多核優化，我測試機器是4核，
一般優化加速比在3倍左右，所以，本文並行化後的理論耗時爲250/3=83.33ms。但我採用OpenMP對本文算法進行優化後達不到3倍的加速比。仍是須要尋找更好的思路。

　　當時看到這個評論後，真的以爲這博友是否是搞錯了，這麼大的圖像，怎麼可能只要33ms就處理完了呢，就是最簡單的一個圖像處理算法，反色（Invert）通過極度優化後也須要大概7/8毫秒的，因此我當時心裏是不承認這個速度的。

　　後續我也在考慮二值圖像的這個特殊性，曾經有考慮過好比膨脹時，遇到有個是白色的像素則中止循環，也考慮過使用直方圖的方式進行優化，畢竟直方圖也只有兩個像素了，可是也仍是達不到上述速度，有些甚至還更慢。因此後續一直也沒有什麼進步。

　　前幾日，網友LQC-Jack忽然又再次提到了這個問題，他認爲針對這個問題確實有更快的方法，畢竟二值得特殊性擺在那裏： 

     

　　其中的「你box濾波的，sum>0當前點就是255」  這個是關鍵，是啊，針對二值圖求局部矩形內的最大值，和求二值圖像的局部均值若是咱們可以創建起聯繫，那麼就能夠藉助於快速的局部均值算法間接的實現腐蝕或膨脹，我在博客裏有多篇文章提到了局部均值的終極優化，特別是SSE圖像算法優化系列十三：超高速BoxBlur算法的實現和優化（Opencv的速度的五倍）一文中提到的方式，效率及其高，針對4096X8192的二值圖也就是30ms左右能搞定。但願燃起。

　　那如何將二者搭橋呢，仔細想一想確實很簡單，若是是求最大值（膨脹），那麼只要局部有一個像素爲255，結果就爲255，此時的局部均值必然大於0 （考慮實際因素，應該是局部累加值，由於考慮最後的整除，不排除某個局部區域，只有一個白點，當局部過大時，整除後的結果可能也爲0），而只有全部局部內的像素都爲0是，最大值才爲0，這個時候 的局部累加值也必然爲0。若是是求最,小值（腐蝕），只要局部有一個像素爲0值，結果就爲0，只有局部全部像素都爲255，結果才爲255，那麼這裏的信息反饋到局部均值就等同於說平均值爲255，則結果爲255，不然結果就爲0（一樣的道理，這裏實際編程時要用局部累加值，而不是平均值）。

　　如此一來，咱們會發現，這種實現過程相比標準的方框模糊來講還少了一些步驟，咱們先貼下我SSE優化方框模糊的核心部分：

 1 int BlockSize = 4, Block = (Width - 1) / BlockSize;
 2 __m128i OldSum = _mm_set1_epi32(LastSum);
 3 __m128 Inv128 = _mm_set1_ps(Inv);
 4 for (int X = 1; X < Block * BlockSize + 1; X += BlockSize)
 5 {
 6     __m128i ColValueOut = _mm_loadu_si128((__m128i *)(ColValue + X - 1));
 7     __m128i ColValueIn = _mm_loadu_si128((__m128i *)(ColValue + X + Radius + Radius));
 8     __m128i ColValueDiff = _mm_sub_epi32(ColValueIn, ColValueOut);                            //    P3 P2 P1 P0                                                
 9     __m128i Value_Temp = _mm_add_epi32(ColValueDiff, _mm_slli_si128(ColValueDiff, 4));        //    P3+P2 P2+P1 P1+P0 P0
10     __m128i Value = _mm_add_epi32(Value_Temp, _mm_slli_si128(Value_Temp, 8));                 //    P3+P2+P1+P0 P2+P1+P0 P1+P0 P0
11     __m128i NewSum = _mm_add_epi32(OldSum, Value);
12     OldSum = _mm_shuffle_epi32(NewSum, _MM_SHUFFLE(3, 3, 3, 3));                              //    從新賦值爲最新值
13     __m128 Mean = _mm_mul_ps(_mm_cvtepi32_ps(NewSum), Inv128);
14     _mm_storesi128_4char(_mm_cvtps_epi32(Mean), LinePD + X);
15 }

　　注意第14及第15行爲求均值並最終保存數據到內存中的過程，其中的NewSum中保存即爲累加的值，注意這裏由於有除法，因此借用了浮點版本的相關指令，同時增長了相關的類型轉換過程。

　　這裏，爲了適應咱們的腐蝕和膨脹的需求，這兩句是不須要的，按照上述分析，好比膨脹效果，只需做以下改動：

LinePD[X + 0] = NewSum.m128i_i32[0] > 0 ? 255 : 0;
LinePD[X + 1] = NewSum.m128i_i32[1] > 0 ? 255 : 0;
LinePD[X + 2] = NewSum.m128i_i32[2] > 0 ? 255 : 0;
LinePD[X + 3] = NewSum.m128i_i32[3] > 0 ? 255 : 0;

　　以上代碼只是示意，若是真的這樣寫，會破壞SSE算法的總體的和諧性，並且這種SSE中穿插普通C代碼會帶來性能上的極大損失，一種處理方法以下所示：

__m128i Flag = _mm_cmpgt_epi32(NewSum, _mm_setzero_si128());
Flag = _mm_packs_epi32(Flag, Flag);
*((int *)(LinePD + X)) = _mm_cvtsi128_si32(_mm_packs_epi16(Flag, Flag));

　　咱們利用SSE中比較運算符的特殊性，產生諸如0XFFFFFFFF這樣的結果，而後在經過有關飽和運算將他們減低到8位，注意上面使用的都是有符號的飽和計算。

　　對於腐蝕的過程，你知道怎麼寫嗎？

　　咱們通過簡單測試，處理一副4096X8192大小的二值圖，任意的半徑大小，耗時基本穩定在24ms左右，比boxblur也快了不少。

　　我也構思過不實用累加和的方式判斷，好比使用或運算或者與運算，可是都是解決不了進出像素的處理問題，所以，總體看來是仍是用累加最爲科學。

　　其實對於半徑比較小時，仍是有更爲快速的方法的，這裏稍微簡單描述下，可是可能不少人看不懂。

　　在咱們上述的實現中，咱們用的是int類型的數據來保存累加值，這是由於半徑稍微大一點累加值就可能超過short類型所能表達的範圍，可是int類型SSE一次只能處理4個，而short類型數據SSE一次能處理8個，所以，若是作適當的變更是否有可能使用short類型呢，是用可能的。

　　由於是二值圖，因此就只有0和255兩個值，0值無所謂，那若是咱們把255這個值修改爲1，那麼在半徑不大於某個數值（64仍是其餘數，能夠本身畫一畫）時，累加值將可控在short類型所能表達的範圍。

　　這是還有個問題就是，255這個值如何變爲1，若是使用_mm_blendv_epi8集合有關判斷語句是能夠實現的，可是這個Blend是比較耗時的，反而得不償失。一個最好的辦法就是充分利用無符號和有符號數之間的特色，當咱們把一個等於255的unsigned char數據類型強制轉換爲signed char時，他的值就等於-1，和咱們要的值1相反， 這個時候咱們本來代碼裏的_mm_add_epi8接可使用_mm_sub_epi8代替，反之亦然。而在SSE裏，這種類型轉換還不須要強制進行，由於他直接操做內存。

　　咱們貼下下面的代碼可能有人就能明白是什麼意思了。

memset(ColValue + Radius, 0, Width * sizeof(unsigned char));
for (int Z = -Radius; Z <= Radius; Z++)
{
    unsigned char *LinePS = Src + ColOffset[Z + Radius] * Stride;
    int BlockSize = 16, Block = Width / BlockSize;
    for (int X = 0; X < Block * BlockSize; X += BlockSize)
    {
        unsigned char *DestP = ColValue + X + Radius;
        __m128i Sample = _mm_loadu_si128((__m128i *)(LinePS + X));                //    255成爲-1
        _mm_storeu_si128((__m128i *)DestP, _mm_sub_epi8(_mm_loadu_si128((__m128i *)DestP), Sample));
    }
    for (int X = Block * BlockSize; X < Width; X++)
    {
        ColValue[X + Radius] += (LinePS[X] == 255 ? 1 : 0);                                            //    更新列數據
    }
}

　　普通的C代碼部分及時直接實現，而SSE部分，並無看到明顯的255到1之間的轉換，一塊兒都在那幾句簡簡單單的代碼中。

　　經過這種相關的優化，大概4096X8192的圖能作到12到13毫秒之間，已經徹底超過了Halcon的速度。

　　halcon中的腐蝕和膨脹也有圓形半徑的，一樣的半徑下圓形半徑在halcon中的耗時大概是矩形半徑的8倍左右，我相信halcon的圓形半徑的算法也是經過EDM算法來實現的，詳見SSE圖像算法優化系列二十五:二值圖像的Euclidean distance map（EDM)特徵圖計算及其優化一文， 而我這裏也差不都是這樣的時間比例。

　　源代碼下載：https://files.cnblogs.com/files/Imageshop/FastBlur.rar

　　極度優化版本工程：https://files.cnblogs.com/files/Imageshop/SSE_Optimization_Demo.rar，見Binary->Processing->Erode/Dilate菜單。