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

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

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

　　我在某篇博客里看到說(shuō)intrinsics語(yǔ)法雖然簡(jiǎn)化了SSE編程的難度，但是他無(wú)法直接控制XMM0-XMM7寄存器，很多指令中間都會(huì )用內存做中轉，所以我就想我如果直接用匯編寫(xiě)效率肯定還能有進(jìn)一步的提高，于是我首先嘗試把GaussBlurFromLeftToRight_SSE優(yōu)化，仔細觀(guān)察這個(gè)函數，如果弄得好，確實(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è)寄存器都充分利用了，在我的預想中應該能有速度的提升的，可是一執行，真的好悲劇，和原先相比速度毫無(wú)變化，這是怎么回事呢。

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

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

　　這個(gè)嘗試純粹是隨意而為，誰(shuí)知道居然非常有效果，具體而言就是在GaussBlurFromLeftToRight_SSE和GaussBlurFromRightToLeft_SSE函數的Y循環(huán)內部，一次性處理二行代碼，我們以L(fǎng)eftToRight為例，示意代碼如下：

　　__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)重復數據

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

　　}

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

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

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

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

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

　　可見(jiàn)，上述過(guò)程存在的一定的精度損失，因為在行方向的處理完成后的浮點(diǎn)到字節數據的轉換丟失了部分數據。但是考慮到是模糊，這種丟失對于結果在視覺(jué)上是基本察覺(jué)不到的。因此，是可以接受的，測試表明，純C版本的這種做法和純C版本的標準做法在速度上基本相當。

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

　　但是垂直方向呢，如果按照上述方式處理，就無(wú)法利用到SSE的優(yōu)勢了，因為上述方式要求每次都是隔行取樣，Cache miss的可能性太高，那么還能不能利用我們在高斯模糊算法的全面優(yōu)化過(guò)程分享(一)提高的那種方式呢。

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

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

　　還有一個(gè)問(wèn)題就是，這樣做會(huì )增加很多次自己數據到浮點(diǎn)數據的轉換，這種轉換的耗時(shí)是否對最后的結果又重要的影響呢，只有實(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行緩沖區

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

　　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); // 分開(kāi)來(lái)執行速度比寫(xiě)在一起有要快些

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

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

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

　　}

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

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

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

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

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

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

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

　　在上一篇文章中，我們提到了由于float類(lèi)型的精度問(wèn)題，當模糊的半徑較大時(shí)，算法的結果會(huì )出現很大的瑕疵，一種方式就是用double類(lèi)型來(lái)解決，還有一種方式就是可以用Deriche濾波器來(lái)解決，為了完美解決這個(gè)問(wèn)題，我還是恨著(zhù)頭皮用SSE實(shí)現了Deriche濾波器，這里簡(jiǎn)要說(shuō)明如下：

　　Deriche濾波器和高斯濾波器有很多類(lèi)似的地方：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.

　　按照維基的解釋?zhuān)珼eriche濾波器可按照如下的步驟執行：詳見(jiàn)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)他們也是行列可分離的算法。

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

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

　　第六：多線(xiàn)程

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

　　第七：比較

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

　　總結：有心就有成績(jì)

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

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

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

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

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

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

　　拒絕無(wú)腦索取代碼。