高斯模糊算法的全面優(yōu)化過程分享（二）

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

　　第一個(gè)嘗試 直接使用內(nèi)聯(lián)匯編替代intrinsics代碼(無效)

　　我在某篇博客里看到說intrinsics語法雖然簡化了SSE編程的難度，但是他無法直接控制XMM0-XMM7寄存器，很多指令中間都會(huì)用內(nèi)存做中轉(zhuǎn)，所以我就想我如果直接用匯編寫效率肯定還能有進(jìn)一步的提高，于是我首先嘗試把GaussBlurFromLeftToRight_SSE優(yōu)化，仔細(xì)觀察這個(gè)函數(shù)，如果弄得好，確實(shí)能有效的利用這些寄存器，有關(guān)代碼如下：

　　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個(gè)寄存器都充分利用了，在我的預(yù)想中應(yīng)該能有速度的提升的，可是一執(zhí)行，真的好悲劇，和原先相比速度毫無變化，這是怎么回事呢。

　　后來我反編譯intrinsics的相關(guān)代碼，發(fā)現(xiàn)編譯器真的很厲害，他的匯編代碼和我上面的基本一致，只是寄存器的利用順序有所不同而已，后面又看了其他的幾個(gè)函數(shù)，發(fā)現(xiàn)編譯器的匯編碼都寫的非常高效，基本上我們是超不過他了，而且編譯器還能充分調(diào)整指令執(zhí)行的順序，使得有關(guān)指令還能實(shí)現(xiàn)指令層次的并行，而如果我們自己寫ASM，這個(gè)對功力的要求就更高了，所以說網(wǎng)絡(luò)上的說法也不可以完全相信，而如果不是有特別強(qiáng)的匯編能力，也不要去挑戰(zhàn)編譯器。

　　第二個(gè)嘗試 水平方向的模糊一次執(zhí)行二行(15%提速)

　　這個(gè)嘗試純粹是隨意而為，誰知道居然非常有效果，具體而言就是在GaussBlurFromLeftToRight_SSE和GaussBlurFromRightToLeft_SSE函數(shù)的Y循環(huán)內(nèi)部，一次性處理二行代碼，我們以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); // 起點(diǎn)重復(fù)數(shù)據(jù)

　　__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);

　　}

　　就是把原來的代碼復(fù)制一份，在稍微調(diào)整一下，當(dāng)然注意這個(gè)時(shí)候Y分量一次要遞增2行了，還有如果Height是奇數(shù)，還要對最后一行做處理。這些活都是細(xì)活，稍微注意就不會(huì)出錯(cuò)了。

　　就這么樣的簡單的一個(gè)調(diào)整，經(jīng)過測試性能居然能有15%的提升，真是意想不到，分析具體的原因，我覺得Y循環(huán)變量的計(jì)數(shù)耗時(shí)的減少在這里是無關(guān)緊要的，核心可能還是這些intrinsics內(nèi)部寄存器的一些調(diào)整，是的更多的指令能并行執(zhí)行。

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

　　第三種嘗試：不使用中間內(nèi)存實(shí)現(xiàn)的近似效果(80%提速)

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

　　可見，上述過程存在的一定的精度損失，因?yàn)樵谛蟹较虻奶幚硗瓿珊蟮母↑c(diǎn)到字節(jié)數(shù)據(jù)的轉(zhuǎn)換丟失了部分?jǐn)?shù)據(jù)。但是考慮到是模糊，這種丟失對于結(jié)果在視覺上是基本察覺不到的。因此，是可以接受的，測試表明，純C版本的這種做法和純C版本的標(biāo)準(zhǔn)做法在速度上基本相當(dāng)。

　　我們考慮這種做法的SSE優(yōu)化，第一，是水平方向的處理，想想很簡單，核心的代碼和之前的是沒有區(qū)別的，當(dāng)然我們也應(yīng)該帶上我們的兩行一次性處理這種訣竅的。

　　但是垂直方向呢，如果按照上述方式處理，就無法利用到SSE的優(yōu)勢了，因?yàn)樯鲜龇绞揭竺看味际歉粜腥?，Cache miss的可能性太高，那么還能不能利用我們在高斯模糊算法的全面優(yōu)化過程分享(一)提高的那種方式呢。

　　仔細(xì)看看(一)中的過程，很明顯他一次性只會(huì)利用到4行的數(shù)據(jù)，同時(shí)，相鄰兩行的處理數(shù)據(jù)有3行是重疊的，那么這就為我們的低內(nèi)存占用同時(shí)又能高效的利用SSE提供了可能性，我們只需要分配4行的浮點(diǎn)緩存區(qū)，然后每次交換行行之間的指針，對垂直方向的處理就能利用相同的SIMD優(yōu)化算法了。

　　但是這樣做又會(huì)帶來另外一個(gè)小問題，就是在Top到Bottom的處理過程中，每一行處理完后又會(huì)有一個(gè)浮點(diǎn)到字節(jié)數(shù)據(jù)的精度丟失，這種丟失經(jīng)過測試也是可以接受的。

　　還有一個(gè)問題就是，這樣做會(huì)增加很多次自己數(shù)據(jù)到浮點(diǎn)數(shù)據(jù)的轉(zhuǎn)換，這種轉(zhuǎn)換的耗時(shí)是否對最后的結(jié)果又重要的影響呢，只有實(shí)測才知道。我們待會(huì)再分析，這里貼出這種近似的優(yōu)化的有關(guān)代碼：

　　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行緩沖區(qū)

　　// 行方向的優(yōu)化，這個(gè)是沒有啥精度損失的

　　for (; Y < Height - 1; Y += 2) // 兩行執(zhí)行的代碼比單行快

　　{

　　ConvertBGR8U2BGRAF_Line_SSE(Src + (Y + 0) * Stride, Buffer, Width);

　　ConvertBGR8U2BGRAF_Line_SSE(Src + (Y + 1) * Stride, Buffer + Width * 4, Width); // 讀取兩行數(shù)據(jù)

　　GaussBlurLeftToRight_TwoLine_SSE(Buffer, Width, B0, B1, B2, B3); // 分開來執(zhí)行速度比寫在一起有要快些

　　GaussBlurRightToLeft_TwoLine_SSE(Buffer, Width, B0, B1, B2, B3);

　　ConvertBGRAF2BGR8U_Line_SSE(Buffer, Dest + (Y + 0) * Stride, Width); // 浮點(diǎn)轉(zhuǎn)換為字節(jié)數(shù)據(jù)

　　ConvertBGRAF2BGR8U_Line_SSE(Buffer, Dest + (Y + 1) * Stride, Width);

　　}

　　for (; Y < Height; Y++) // 執(zhí)行剩下的單行

　　{

　　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);

　　}

　　// 列方向考慮優(yōu)化，多了一次浮點(diǎn)到字節(jié)類型的轉(zhuǎn)換，有精度損失

　　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); // 轉(zhuǎn)換當(dāng)前行到浮點(diǎn)緩存

　　GaussBlurTopToBottom_LowMemory_SSE(Line0, Line1, Line2, Line3, Width, B0, B1, B2, B3); // 垂直方向處理

　　ConvertBGRAF2BGR8U_Line_SSE(Line3, Dest + Y * Stride, Width); // 又再次轉(zhuǎn)換為字節(jié)數(shù)據(jù)

　　Temp = Line0; Line0 = Line1; Line1 = Line2; Line2 = Line3; Line3 = Temp; // 交換行緩存

　　}

　　ConvertBGR8U2BGRAF_Line_SSE(Dest + (Height - 1) * Stride, Buffer + 3 * Width * 4, Width); // 重復(fù)邊緣像素

　　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是之前的相關(guān)函數(shù)的單行版。

　　經(jīng)過測試，上述改進(jìn)后的算法在同樣配置的電腦上，針對3000*2000的彩色圖像耗時(shí)約為86ms，和之前的145ms相比，提速了近一倍，而基本不占用額外的內(nèi)存，可是為什么呢，似乎代碼中還增加了很多浮點(diǎn)到字節(jié)和字節(jié)到浮點(diǎn)數(shù)據(jù)的轉(zhuǎn)換代碼，總的計(jì)算量應(yīng)該是增加的啊。按照我的分析，我認(rèn)為這是這里分配的輔助內(nèi)存很小，被分配到一級緩存或者二級緩存或其他更靠近CPU的位置的內(nèi)尺寸區(qū)域的可能性更大，而第一版本的內(nèi)存由于過大，只可能分配堆棧中，同時(shí)我們算法里有著大量訪問內(nèi)存的地方，這樣雖然總的轉(zhuǎn)換量增加了，但是內(nèi)存訪問節(jié)省的時(shí)間已經(jīng)超越了轉(zhuǎn)換增加的時(shí)間了。

　　第四種嘗試：列方向直接使用BGR而不是BGRA的SSE優(yōu)化(100%提速)

　　在高斯模糊算法的全面優(yōu)化過程分享(一)中，為了解決水平方向上的SSE優(yōu)化問題，我們將BGR數(shù)據(jù)轉(zhuǎn)換為了BGRA格式的浮點(diǎn)數(shù)后再進(jìn)行處理，這樣在列方向處理時(shí)同樣需要處理A的數(shù)據(jù)，但是在經(jīng)過第三種嘗試后，在垂直方向的處理我們還有必要處理這個(gè)多余的A嗎，當(dāng)然沒有必要，這樣垂直方向整體上又可以減少約25%的時(shí)間，耗時(shí)只有75ms左右了，實(shí)現(xiàn)了約100%的提速。

　　第五種嘗試：算法穩(wěn)定性的考慮和最終的妥協(xié)

　　在上一篇文章中，我們提到了由于float類型的精度問題，當(dāng)模糊的半徑較大時(shí)，算法的結(jié)果會(huì)出現(xiàn)很大的瑕疵，一種方式就是用double類型來解決，還有一種方式就是可以用Deriche濾波器來解決，為了完美解決這個(gè)問題，我還是恨著頭皮用SSE實(shí)現(xiàn)了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濾波器可按照如下的步驟執(zhí)行：詳見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:

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

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

　　由于計(jì)算量較原先的高斯有所增加，這里最后的優(yōu)化完成的耗時(shí)約為100ms。

　　第六：多線程

　　在水平方向計(jì)算時(shí)，各行之間的計(jì)算時(shí)獨(dú)立的，因此是可以并行處理的，但是垂直方向由于是前后依賴的，是無法并行的。比如用openmp做2個(gè)線程的并行，大概速度能提高到(高斯)到60ms，但是這個(gè)東西在不是本文這里的重點(diǎn)。

　　第七：比較

　　同標(biāo)準(zhǔn)的高斯濾波相比，Deriche濾波器由于其特性，無法支持In-Place操作，也就是說Src和Dest不能相同，如果一定要相同，就只有通過一個(gè)中間內(nèi)存來過渡了，而標(biāo)準(zhǔn)高斯是可以的。第二就是高斯是可以不占用太多額外的內(nèi)存就可以實(shí)現(xiàn)的，而Deriche需要一份同樣大小的內(nèi)存。第三就是標(biāo)準(zhǔn)高斯速度還是要快一點(diǎn)。第四Deriche濾波器的精度在float類型時(shí)精度要比標(biāo)準(zhǔn)高斯高。綜合選擇，我覺得還是以后選擇Deriche代替標(biāo)準(zhǔn)的高斯模糊。

　　總結(jié)：有心就有成績

　　同opencv的cvsmooth相比，同樣的機(jī)器上同樣的3000*2000大小的彩圖，Ksize我取100時(shí)，需要1200多ms，比我這里慢了10倍，但是cvsmooth似乎對ksize參數(shù)敏感，他并不是與核大小無關(guān)的，ksize較小時(shí)還會(huì)很快的，不過除了一些特效外，在很多場合我們其實(shí)需要大半徑的高斯的(比如圖像增強(qiáng)、銳化上)。

　　做完了在回頭看看優(yōu)化的過程，覺得和看書是一個(gè)道理，先是越看越厚，通了就像一張薄紙一樣。

　　最后總結(jié)下，就是一件事情，只要你有時(shí)間和信心，就能有進(jìn)步，堅(jiān)持是勝利的必要條件。

　　提供一個(gè)測試的Demo: http://share.eepw.com.cn/share/download/id/383731

　　由測試Demo可以測試出，當(dāng)選擇低內(nèi)存近似版本或者準(zhǔn)確版本時(shí)，當(dāng)半徑較大時(shí)，如果連續(xù)的拖動(dòng)滾動(dòng)條，圖像會(huì)出現(xiàn)閃爍，而如果選擇Deriche時(shí)，則圖像變換很平緩，而當(dāng)半徑特別大時(shí)，如果選擇低內(nèi)存近似版本或者準(zhǔn)確版本，則圖像有可能會(huì)出現(xiàn)線條或者色塊，只有Deriche濾波的效果是完美的。

　　高斯模糊的優(yōu)化到此結(jié)束，如果有誰有用GPU實(shí)現(xiàn)的，還請告訴我下大概的耗時(shí)。

　　拒絕無腦索取代碼。