高斯模糊算法的全面優化過程分享（二）

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

　　第一個嘗試 直接使用內聯匯編替代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是奇數，還要對最后一行做處理。這些活都是細活，稍微注意就不會出錯了。

　　就這么樣的簡單的一個調整，經過測試性能居然能有15%的提升，真是意想不到，分析具體的原因，我覺得Y循環變量的計數耗時的減少在這里是無關緊要的，核心可能還是這些intrinsics內部寄存器的一些調整，是的更多的指令能并行執行。

　　但是，在垂直方向的SSE代碼用類似的方式調整似乎沒有性能的提升，還會到底代碼的可讀性較差。

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

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

　　可見，上述過程存在的一定的精度損失，因為在行方向的處理完成后的浮點到字節數據的轉換丟失了部分數據。但是考慮到是模糊，這種丟失對于結果在視覺上是基本察覺不到的。因此，是可以接受的，測試表明，純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:

　　and from right to left according to the formula:

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

　　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:

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

　　同樣為了節省內存，我們也使用了類似上述第三種和第四重嘗試的方式，但是考慮到Deriche的特殊性(主要是這里)，他還是需要一份中間內存的，為了效率和內存，我們再次以犧牲精度為準備，中間使用了一份和圖像一樣的字節數據內存。

　　由于計算量較原先的高斯有所增加，這里最后的優化完成的耗時約為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://share.eepw.com.cn/share/download/id/383731

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

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

　　拒絕無腦索取代碼。