超越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優化方框模糊的核心部分：

int BlockSize = 4, Block = (Width - 1) / BlockSize;
__m128i OldSum = _mm_set1_epi32(LastSum);
__m128 Inv128 = _mm_set1_ps(Inv);
for (int X = 1; X < Block * BlockSize + 1; X += BlockSize)
{
    __m128i ColValueOut = _mm_loadu_si128((__m128i *)(ColValue + X - 1));
    __m128i ColValueIn = _mm_loadu_si128((__m128i *)(ColValue + X + Radius + Radius));
    __m128i ColValueDiff = _mm_sub_epi32(ColValueIn, ColValueOut);                            //    P3 P2 P1 P0                                                
    __m128i Value_Temp = _mm_add_epi32(ColValueDiff, _mm_slli_si128(ColValueDiff, 4));        //    P3+P2 P2+P1 P1+P0 P0
    __m128i Value = _mm_add_epi32(Value_Temp, _mm_slli_si128(Value_Temp, 8));                 //    P3+P2+P1+P0 P2+P1+P0 P1+P0 P0
    __m128i NewSum = _mm_add_epi32(OldSum, Value);
    OldSum = _mm_shuffle_epi32(NewSum, _MM_SHUFFLE(3, 3, 3, 3));                              //    重新賦值為最新值
    __m128 Mean = _mm_mul_ps(_mm_cvtepi32_ps(NewSum), Inv128);
    _mm_storesi128_4char(_mm_cvtps_epi32(Mean), LinePD + X);
}

　　注意第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/SSE_Optimization_Demo.rar，見Binary->Processing->Erode/Dilate菜單。