SSE圖像算法優化系列二：高斯模糊算法的全面優化過程分享（二）。

時間 2019-11-08

標籤 sse 圖像算法優化系列高斯模糊全面過程分享简体版

原文原文鏈接

在高斯模糊算法的全面優化過程分享（一）一文中咱們已經給出了一種至關高性能的高斯模糊過程，可是優化沒有終點，通過上一個星期的發憤圖強和測試，對該算法的效率提高又有了一個新的高度，這裏把優化過程當中的一些心得和收穫用文字的形式記錄下來。算法

第一個嘗試直接使用內聯彙編替代intrinsics代碼（無效）

我在某篇博客裏看到說intrinsics語法雖然簡化了SSE編程的難度，可是他沒法直接控制XMM0-XMM7寄存器，不少指令中間都會用內存作中轉，因此我就想我若是直接用匯編寫效率確定還能有進一步的提升，因而我首先嚐試把GaussBlurFromLeftToRight_SSE優化，仔細觀察這個函數，若是弄得好，確實能有效的利用這些寄存器，有關代碼以下：編程

void GaussBlurFromLeftToRight_SSE(float *Data, int Width, int Height, float B0, float B1, float B2, float B3)
{
    float *MemB3 = (float *)_mm_malloc(4 * sizeof(float), 16);
    MemB3[0] = MemB3[1] = MemB3[2] = MemB3[3] = B3;
    int Stride = Width * 4 * sizeof(float);
    _asm
    {
        mov     ecx, Height
        movss   xmm0, B0
        shufps  xmm0, xmm0, 0            //    xmm0 = B0
        movss   xmm1, B1
        shufps  xmm1, xmm1, 0            //    xmm1 = B1
        movss   xmm2, B2
        shufps  xmm2, xmm2, 0            //    xmm2 = B2
        mov     edi, MemB3
    LoopH24 :
        mov     esi, ecx
        imul    esi, Stride
        add     esi, Data                //    LinePD = Data + Y * Width * 4
        mov     eax, Width
        movaps  xmm3, [esi]              //    xmm3 = V1
        movaps  xmm4, xmm3                //  xmm4 = V2 = V1
        movaps  xmm5, xmm3                //     xmm5 = V3 = V1
    LoopW24 :
    movaps  xmm6, [esi]                //    xmm6 = V0
        movaps  xmm7, xmm3                //    xmm7 = V1
        mulps   xmm5, [edi]                //    xmm5 = V3 * B3
        mulps   xmm7, xmm1                //    xmm7 = V1 * B1
        mulps   xmm6, xmm0                //    xmm6 = V0 * B0
        addps   xmm6, xmm7                //    xmm6 = V0 * B0 + V1 * B1
        movaps  xmm7, xmm4                //    xmm7 = V2
        mulps   xmm7, xmm2                //    xmm7 = V2 * B2
        addps   xmm5, xmm7                //    xmm5 = V3 * B3 + V2 * B2
        addps   xmm6, xmm5                //    xmm6 = V0 * B0 + V1 * B1 + V3 * B3 + V2 * B2
        movaps  xmm5, xmm4                //    V3 = V2            
        movaps  xmm4, xmm3                //    V2 = V1
        movaps [esi], xmm6
        movaps  xmm3, xmm6                //    V1 = V0
        add     esi, 16
        dec     eax
        jnz     LoopW24
        dec     ecx
        jnz     LoopH24
    }
    _mm_free(MemB3);
}

　　看上面的代碼，基本上把XMM0-XMM7這8個寄存器都充分利用了，在個人預想中應該能有速度的提高的，但是一執行，真的好悲劇，和原先相比速度毫無變化，這是怎麼回事呢。緩存

後來我反編譯intrinsics的相關代碼，發現編譯器真的很厲害，他的彙編代碼和我上面的基本一致，只是寄存器的利用順序有所不一樣而已，後面又看了其餘的幾個函數，發現編譯器的彙編碼都寫的很是高效，基本上咱們是超不過他了，並且編譯器還能充分調整指令執行的順序，使得有關指令還能實現指令層次的並行，而若是咱們本身寫ASM，這個對功力的要求就更高了，因此說網絡上的說法也不能夠徹底相信，而若是不是有特別強的彙編能力，也不要去挑戰編譯器。網絡

第二個嘗試水平方向的模糊一次執行二行（15%提速）

這個嘗試純粹是隨意而爲，誰知道竟然很是有效果，具體而言就是在GaussBlurFromLeftToRight_SSE和GaussBlurFromRightToLeft_SSE函數的Y循環內部，一次性處理二行代碼，咱們以LeftToRight爲例，示意代碼以下：多線程

    __m128 CofB0 = _mm_set_ps(0, B0, B0, B0);
    __m128 CofB1 = _mm_set_ps(0, B1, B1, B1);
    __m128 CofB2 = _mm_set_ps(0, B2, B2, B2);
    __m128 CofB3 = _mm_set_ps(0, B3, B3, B3);
    __m128 V1 = _mm_load_ps(LineP1);                //    起點重複數據
    __m128 W1 = _mm_load_ps(LineP2);
    __m128 V2 = V1, V3 = V1;
    __m128 W2 = W1, W3 = W1;
    for (int X = 0; X < Length; X++, LineP1 += 4, LineP2 += 4)            
    {
        __m128 V0 = _mm_load_ps(LineP1);
        __m128 W0 = _mm_load_ps(LineP2);
        __m128 V01 = _mm_add_ps(_mm_mul_ps(CofB0, V0), _mm_mul_ps(CofB1, V1));
        __m128 W01 = _mm_add_ps(_mm_mul_ps(CofB0, W0), _mm_mul_ps(CofB1, W1));
        __m128 V23 = _mm_add_ps(_mm_mul_ps(CofB2, V2), _mm_mul_ps(CofB3, V3));
        __m128 W23 = _mm_add_ps(_mm_mul_ps(CofB2, W2), _mm_mul_ps(CofB3, W3));
        __m128 V = _mm_add_ps(V01, V23);
        __m128 W = _mm_add_ps(W01, W23);
        V3 = V2;    V2 = V1;    V1 = V;
        W3 = W2;    W2 = W1;    W1 = W;
        _mm_store_ps(LineP1, V);
        _mm_store_ps(LineP2, W);
    }

　　就是把原來的代碼複製一份，在稍微調整一下，固然注意這個時候Y份量一次要遞增2行了，還有若是Height是奇數，還要對最後一行作處理。這些活都是細活，稍微注意就不會出錯了。ide

就這麼樣的簡單的一個調整，通過測試性能竟然能有15%的提高，真是意想不到，分析具體的緣由，我以爲Y循環變量的計數耗時的減小在這裏是可有可無的，核心可能仍是這些intrinsics內部寄存器的一些調整，是的更多的指令能並行執行。函數

可是，在垂直方向的SSE代碼用相似的方式調整彷佛沒有性能的提高，還會到底代碼的可讀性較差。oop

　　第三種嘗試：不使用中間內存實現的近似效果（80%提速）post

之前我在寫高斯模糊時考慮到內存佔用問題，採用了一種近似的方式，在水平方向計算時，只須要分配一行大小的浮點數據，而後每一行都利用這一行數據作緩存，當一行數據的水平模糊計算完後，就把這些數據轉換爲字節數據保存到結果圖像中，當水平方向都計算完成後，在進行列方向的處理。列方向也是隻分配高度大小的一列中間浮點緩存數據，而後進行高度方向處理，每列處理完後，把浮點的結果轉換成字節數據。

可見，上述過程存在的必定的精度損失，由於在行方向的處理完成後的浮點到字節數據的轉換丟失了部分數據。可是考慮到是模糊，這種丟失對於結果在視覺上是基本察覺不到的。所以，是能夠接受的，測試代表，純C版本的這種作法和純C版本的標準作法在速度上基本至關。

咱們考慮這種作法的SSE優化，第一，是水平方向的處理，想一想很簡單，核心的代碼和以前的是沒有區別的，固然咱們也應該帶上咱們的兩行一次性處理這種訣竅的。

可是垂直方向呢，若是按照上述方式處理，就沒法利用到SSE的優點了，由於上述方式要求每次都是隔行取樣，Cache miss的可能性過高，那麼還能不能利用咱們在高斯模糊算法的全面優化過程分享（一）提升的那種方式呢。

仔細看看（一）中的過程，很明顯他一次性只會利用到4行的數據，同時，相鄰兩行的處理數據有3行是重疊的，那麼這就爲咱們的低內存佔用同時又能高效的利用SSE提供了可能性，咱們只須要分配4行的浮點緩存區，而後每次交換行行之間的指針，對垂直方向的處理就能利用相同的SIMD優化算法了。

可是這樣作又會帶來另一個小問題，就是在Top到Bottom的處理過程當中，每一行處理完後又會有一個浮點到字節數據的精度丟失，這種丟失通過測試也是能夠接受的。

還有一個問題就是，這樣作會增長不少次本身數據到浮點數據的轉換，這種轉換的耗時是否對最後的結果又重要的影響呢，只有實測才知道。咱們待會再分析，這裏貼出這種近似的優化的有關代碼：

 void GaussBlur_FastAndLowMemory_SSE(unsigned char *Src, unsigned char *Dest, int Width, int Height, int Stride, float Radius)
 {
     float B0, B1, B2, B3;                            
     float *Line0, *Line1, *Line2, *Line3, *Temp;
     int Y = 0;
     CalcGaussCof(Radius, B0, B1, B2, B3);

     float *Buffer = (float *)_mm_malloc(Width * 4 * 4 * sizeof(float), 16);                //    最多須要4行緩衝區

     //    行方向的優化，這個是沒有啥精度損失的
     for (; Y < Height - 1; Y += 2)                                                            //    兩行執行的代碼比單行快
     {
         ConvertBGR8U2BGRAF_Line_SSE(Src + (Y + 0) * Stride, Buffer, Width);
         ConvertBGR8U2BGRAF_Line_SSE(Src + (Y + 1) * Stride, Buffer + Width * 4, Width);            //    讀取兩行數據
         GaussBlurLeftToRight_TwoLine_SSE(Buffer, Width, B0, B1, B2, B3);                    //    分開來執行速度比寫在一塊兒有要快些
         GaussBlurRightToLeft_TwoLine_SSE(Buffer, Width, B0, B1, B2, B3);
         ConvertBGRAF2BGR8U_Line_SSE(Buffer, Dest + (Y + 0) * Stride, Width);                    //    浮點轉換爲字節數據
         ConvertBGRAF2BGR8U_Line_SSE(Buffer, Dest + (Y + 1) * Stride, Width);
     }
     for (; Y < Height; Y++)                                                                //    執行剩下的單行
     {
         ConvertBGR8U2BGRAF_Line_SSE(Src + Y * Stride, Buffer, Width);
         GaussBlurLeftToRight_OneLine_SSE(Buffer, Width, B0, B1, B2, B3);
         GaussBlurRightToLeft_OneLine_SSE(Buffer, Width, B0, B1, B2, B3);
         ConvertBGRAF2BGR8U_Line_SSE(Buffer, Dest + Y * Stride, Width);
     }

     //    列方向考慮優化，多了一次浮點到字節類型的轉換，有精度損失
     ConvertBGR8U2BGRAF_Line_SSE(Dest, Buffer + 3 * Width * 4, Width);
     memcpy(Buffer + 0 * Width * 4, Buffer + 3 * Width * 4, Width * 4 * sizeof(float));            //    起始值取邊界的值
     memcpy(Buffer + 1 * Width * 4, Buffer + 3 * Width * 4, Width * 4 * sizeof(float));
     memcpy(Buffer + 2 * Width * 4, Buffer + 3 * Width * 4, Width * 4 * sizeof(float));

     Line0 = Buffer + 0 * Width * 4;    Line1 = Buffer + 1 * Width * 4;
     Line2 = Buffer + 2 * Width * 4;    Line3 = Buffer + 3 * Width * 4;
     for (Y = 0; Y < Height; Y++)
     {
         ConvertBGR8U2BGRAF_Line_SSE(Dest + Y * Stride, Line3, Width);                                //    轉換當前行到浮點緩存
         GaussBlurTopToBottom_LowMemory_SSE(Line0, Line1, Line2, Line3, Width, B0, B1, B2, B3);    //    垂直方向處理
         ConvertBGRAF2BGR8U_Line_SSE(Line3, Dest + Y * Stride, Width);                                //    又再次轉換爲字節數據
         Temp = Line0;    Line0 = Line1;    Line1 = Line2;    Line2 = Line3;    Line3 = Temp;            //    交換行緩存
     }    

     ConvertBGR8U2BGRAF_Line_SSE(Dest + (Height - 1) * Stride, Buffer + 3 * Width * 4, Width);        //    重複邊緣像素
     memcpy(Buffer + 0 * Width * 4, Buffer + 3 * Width * 4, Width * 4 * sizeof(float));
     memcpy(Buffer + 1 * Width * 4, Buffer + 3 * Width * 4, Width * 4 * sizeof(float));
     memcpy(Buffer + 2 * Width * 4, Buffer + 3 * Width * 4, Width * 4 * sizeof(float));

     Line0 = Buffer + 0 * Width * 4;    Line1 = Buffer + 1 * Width * 4;
     Line2 = Buffer + 2 * Width * 4;    Line3 = Buffer + 3 * Width * 4;
     for (Y = Height - 1; Y > 0; Y--)                                                            //    垂直向上處理
     {
         ConvertBGR8U2BGRAF_Line_SSE(Dest + Y * Stride, Line0, Width);
         GaussBlurBottomToTop_LowMemory_SSE(Line0, Line1, Line2, Line3, Width, B0, B1, B2, B3);
         ConvertBGRAF2BGR8U_Line_SSE(Line0, Dest + Y * Stride, Width);
         Temp = Line3;    Line3 = Line2;    Line2 = Line1;    Line1 = Line0;    Line0 = Temp;
     }
     _mm_free(Buffer);
 }

　　上述代碼中的ConvertBGR8U2BGRAF_Line_SSE和ConvertBGRAF2BGR8U_Line_SSE是以前的相關函數的單行版。

通過測試，上述改進後的算法在一樣配置的電腦上，針對3000*2000的彩色圖像耗時約爲86ms，和以前的145ms相比，提速了近一倍，而基本不佔用額外的內存，但是爲何呢，彷佛代碼中還增長了不少浮點到字節和字節到浮點數據的轉換代碼，總的計算量應該是增長的啊。按照個人分析，我認爲這是這裏分配的輔助內存很小，被分配到一級緩存或者二級緩存或其餘更靠近CPU的位置的內尺寸區域的可能性更大，而初版本的內存因爲過大，只可能分配堆棧中，同時咱們算法裏有着大量訪問內存的地方，這樣雖然總的轉換量增長了，可是內存訪問節省的時間已經超越了轉換增長的時間了。

　　第四種嘗試：列方向直接使用BGR而不是BGRA的SSE優化（100%提速）

在高斯模糊算法的全面優化過程分享（一）中，爲了解決水平方向上的SSE優化問題，咱們將BGR數據轉換爲了BGRA格式的浮點數後再進行處理，這樣在列方向處理時一樣須要處理A的數據，可是在通過第三種嘗試後，在垂直方向的處理咱們還有必要處理這個多餘的A嗎，固然沒有必要，這樣垂直方向總體上又能夠減小約25%的時間，耗時只有75ms左右了，實現了約100%的提速。

第五種嘗試：算法穩定性的考慮和最終的妥協

　在上一篇文章中，咱們提到了因爲float類型的精度問題，當模糊的半徑較大時，算法的結果會出現很大的瑕疵，一種方式就是用double類型來解決，還有一種方式就是能夠用Deriche濾波器來解決，爲了完美解決這個問題，我仍是恨着頭皮用SSE實現了Deriche濾波器，這裏簡要說明以下：

　　Deriche濾波器和高斯濾波器有不少相似的地方：The Deriche filter is a smoothing filter (low-pass) which was designed to optimally detect, along with a derivation operator, the contours in an image (Canny criteria optimization). Besides, as this filter is very similar to a gaussian filter, but much simpler to implement (based on simple first order IIR filters), it is also much used for general image filtering.

按照維基的解釋，Deriche濾波器可按照以下的步驟執行：詳見https://en.wikipedia.org/wiki/Deriche_edge_detector。

It's possible to separate the process of obtaining the value of a two-dimensional Deriche filter into two parts. In first part, image array is passed in the horizontal direction from left to right according to the following formula:

y_{{ij}}^{1}=a_{1}x_{{ij}}+a_{2}x_{{ij-1}}+b_{1}y_{{ij-1}}^{1}+b_{2}y_{{ij-2}}^{1}

and from right to left according to the formula:

y_{{ij}}^{2}=a_{3}x_{{ij+1}}+a_{4}x_{{ij+2}}+b_{1}y_{{ij+1}}^{2}+b_{2}y_{{ij+2}}^{2}

The result of the computation is then stored into temporary two-dimensional array:

\theta _{{ij}}=c_{1}(y_{{ij}}^{1}+y_{{ij}}^{2})

The second step of the algorithm is very similar to the first one. The two-dimensional array from the previous step is used as the input. It is then passed in the vertical direction from top to bottom and bottom-up according to the following formulas:

y_{{ij}}^{1}=a_{5}\theta _{{ij}}+a_{6}\theta _{{i-1j}}+b_{1}y_{{i-1j}}^{1}+b_{2}y_{{i-2j}}^{1}

y_{{ij}}^{2}=a_{7}\theta _{{i+1j}}+a_{8}\theta _{{i+2j}}+b_{1}y_{{i+1j}}^{2}+b_{2}y_{{i+2j}}^{2}

\Theta _{{ij}}=c_{2}(y_{{ij}}^{1}+y_{{ij}}^{2})

　　可見他們也是行列可分離的算法。

一樣爲了節省內存，咱們也使用了相似上述第三種和第四重嘗試的方式，可是考慮到Deriche的特殊性（主要是 $\Theta _{{ij}}=c_{2}(y_{{ij}}^{1}+y_{{ij}}^{2})$ 這裏），他仍是須要一份中間內存的，爲了效率和內存，咱們再次以犧牲精度爲準備，中間使用了一份和圖像同樣的字節數據內存。

因爲計算量較原先的高斯有所增長，這裏最後的優化完成的耗時約爲100ms。

第六：多線程

在水平方向計算時，各行之間的計算時獨立的，所以是能夠並行處理的，可是垂直方向因爲是先後依賴的，是沒法並行的。好比用openmp作2個線程的並行，大概速度能提升到（高斯）到60ms，可是這個東西在不是本文這裏的重點。

　第七：比較

同標準的高斯濾波相比，Deriche濾波器因爲其特性，沒法支持In-Place操做，也就是說Src和Dest不能相同，若是必定要相同，就只有經過一箇中間內存來過渡了，而標準高斯是能夠的。第二就是高斯是能夠不佔用太多額外的內存就能夠實現的，而Deriche須要一份一樣大小的內存。第三就是標準高斯速度仍是要快一點。第四Deriche濾波器的精度在float類型時精度要比標準高斯高。綜合選擇，我以爲仍是之後選擇Deriche代替標準的高斯模糊。

總結：有心就有成績

同opencv的cvsmooth相比，一樣的機器上一樣的3000*2000大小的彩圖，Ksize我取100時，須要1200多ms，比我這裏慢了10倍，可是cvsmooth彷佛對ksize參數敏感，他並非與核大小無關的，ksize較小時還會很快的，不過除了一些特效外，在不少場合咱們其實須要大半徑的高斯的（好比圖像加強、銳化上）。

作完了在回頭看看優化的過程，以爲和看書是一個道理，先是越看越厚，通了就像一張薄紙同樣。

最後總結下，就是一件事情，只要你有時間和信心，就能有進步，堅持是勝利的必要條件。

提供一個測試的Demo: http://files.cnblogs.com/files/Imageshop/FastGaussBlur.rar

由測試Demo能夠測試出，當選擇低內存近似版本或者準確版本時，當半徑較大時，若是連續的拖動滾動條，圖像會出現閃爍，而若是選擇Deriche時，則圖像變換很平緩，而當半徑特別大時，若是選擇低內存近似版本或者準確版本，則圖像有可能會出現線條或者色塊，只有Deriche濾波的效果是完美的。

　　高斯模糊的優化到此結束，若是有誰有用GPU實現的，還請告訴我下大概的耗時。

拒絕無腦索取代碼。