高斯模糊算法的全面優(yōu)化過程分享（二）
加入技術交流群
在高斯模糊算法的全面優(yōu)化過程分享(一)一文中我們已經(jīng)給出了一種相當高性能的高斯模糊過程，但是優(yōu)化沒有終點，經(jīng)過上一個星期的發(fā)憤圖強和測試，對該算法的效率提升又有了一個新的高度，這里把優(yōu)化過程中的一些心得和收獲用文字的形式記錄下來。
本文引用地址：http://butianyuan.cn/article/201703/345017.htm第一個嘗試 直接使用內聯(lián)匯編替代intrinsics代碼(無效)
我在某篇博客里看到說intrinsics語法雖然簡化了SSE編程的難度，但是他無法直接控制XMM0-XMM7寄存器，很多指令中間都會用內存做中轉，所以我就想我如果直接用匯編寫效率肯定還能有進一步的提高，于是我首先嘗試把GaussBlurFromLeftToRight_SSE優(yōu)化，仔細觀察這個函數(shù)，如果弄得好，確實能有效的利用這些寄存器，有關代碼如下：
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個寄存器都充分利用了，在我的預想中應該能有速度的提升的，可是一執(zhí)行，真的好悲劇，和原先相比速度毫無變化，這是怎么回事呢。
后來我反編譯intrinsics的相關代碼，發(fā)現(xiàn)編譯器真的很厲害，他的匯編代碼和我上面的基本一致，只是寄存器的利用順序有所不同而已，后面又看了其他的幾個函數(shù)，發(fā)現(xiàn)編譯器的匯編碼都寫的非常高效，基本上我們是超不過他了，而且編譯器還能充分調整指令執(zhí)行的順序，使得有關指令還能實現(xiàn)指令層次的并行，而如果我們自己寫ASM，這個對功力的要求就更高了，所以說網(wǎng)絡上的說法也不可以完全相信，而如果不是有特別強的匯編能力，也不要去挑戰(zhàn)編譯器。
第二個嘗試 水平方向的模糊一次執(zhí)行二行(15%提速)
這個嘗試純粹是隨意而為，誰知道居然非常有效果，具體而言就是在GaussBlurFromLeftToRight_SSE和GaussBlurFromRightToLeft_SSE函數(shù)的Y循環(huán)內部，一次性處理二行代碼，我們以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); // 起點重復數(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);
}
就是把原來的代碼復制一份，在稍微調整一下，當然注意這個時候Y分量一次要遞增2行了，還有如果Height是奇數(shù)，還要對最后一行做處理。這些活都是細活，稍微注意就不會出錯了。
就這么樣的簡單的一個調整，經(jīng)過測試性能居然能有15%的提升，真是意想不到，分析具體的原因，我覺得Y循環(huán)變量的計數(shù)耗時的減少在這里是無關緊要的，核心可能還是這些intrinsics內部寄存器的一些調整，是的更多的指令能并行執(zhí)行。
但是，在垂直方向的SSE代碼用類似的方式調整似乎沒有性能的提升，還會到底代碼的可讀性較差。
第三種嘗試：不使用中間內存實現(xiàn)的近似效果(80%提速)
以前我在寫高斯模糊時考慮到內存占用問題，采用了一種近似的方式，在水平方向計算時，只需要分配一行大小的浮點數(shù)據(jù)，然后每一行都利用這一行數(shù)據(jù)做緩存，當一行數(shù)據(jù)的水平模糊計算完后，就把這些數(shù)據(jù)轉換為字節(jié)數(shù)據(jù)保存到結果圖像中，當水平方向都計算完成后，在進行列方向的處理。列方向也是只分配高度大小的一列中間浮點緩存數(shù)據(jù)，然后進行高度方向處理，每列處理完后，把浮點的結果轉換成字節(jié)數(shù)據(jù)。
可見，上述過程存在的一定的精度損失，因為在行方向的處理完成后的浮點到字節(jié)數(shù)據(jù)的轉換丟失了部分數(shù)據(jù)。但是考慮到是模糊，這種丟失對于結果在視覺上是基本察覺不到的。因此，是可以接受的，測試表明，純C版本的這種做法和純C版本的標準做法在速度上基本相當。
我們考慮這種做法的SSE優(yōu)化，第一，是水平方向的處理，想想很簡單，核心的代碼和之前的是沒有區(qū)別的，當然我們也應該帶上我們的兩行一次性處理這種訣竅的。
但是垂直方向呢，如果按照上述方式處理，就無法利用到SSE的優(yōu)勢了，因為上述方式要求每次都是隔行取樣，Cache miss的可能性太高，那么還能不能利用我們在高斯模糊算法的全面優(yōu)化過程分享(一)提高的那種方式呢。
仔細看看(一)中的過程，很明顯他一次性只會利用到4行的數(shù)據(jù)，同時，相鄰兩行的處理數(shù)據(jù)有3行是重疊的，那么這就為我們的低內存占用同時又能高效的利用SSE提供了可能性，我們只需要分配4行的浮點緩存區(qū)，然后每次交換行行之間的指針，對垂直方向的處理就能利用相同的SIMD優(yōu)化算法了。
但是這樣做又會帶來另外一個小問題，就是在Top到Bottom的處理過程中，每一行處理完后又會有一個浮點到字節(jié)數(shù)據(jù)的精度丟失，這種丟失經(jīng)過測試也是可以接受的。
還有一個問題就是，這樣做會增加很多次自己數(shù)據(jù)到浮點數(shù)據(jù)的轉換，這種轉換的耗時是否對最后的結果又重要的影響呢，只有實測才知道。我們待會再分析，這里貼出這種近似的優(yōu)化的有關代碼：
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)化，這個是沒有啥精度損失的
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); // 浮點轉換為字節(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)化，多了一次浮點到字節(jié)類型的轉換，有精度損失
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); // 又再次轉換為字節(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); // 重復邊緣像素
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是之前的相關函數(shù)的單行版。
經(jīng)過測試，上述改進后的算法在同樣配置的電腦上，針對3000*2000的彩色圖像耗時約為86ms，和之前的145ms相比，提速了近一倍，而基本不占用額外的內存，可是為什么呢，似乎代碼中還增加了很多浮點到字節(jié)和字節(jié)到浮點數(shù)據(jù)的轉換代碼，總的計算量應該是增加的啊。按照我的分析，我認為這是這里分配的輔助內存很小，被分配到一級緩存或者二級緩存或其他更靠近CPU的位置的內尺寸區(qū)域的可能性更大，而第一版本的內存由于過大，只可能分配堆棧中，同時我們算法里有著大量訪問內存的地方，這樣雖然總的轉換量增加了，但是內存訪問節(jié)省的時間已經(jīng)超越了轉換增加的時間了。
第四種嘗試：列方向直接使用BGR而不是BGRA的SSE優(yōu)化(100%提速)
在高斯模糊算法的全面優(yōu)化過程分享(一)中，為了解決水平方向上的SSE優(yōu)化問題，我們將BGR數(shù)據(jù)轉換為了BGRA格式的浮點數(shù)后再進行處理，這樣在列方向處理時同樣需要處理A的數(shù)據(jù)，但是在經(jīng)過第三種嘗試后，在垂直方向的處理我們還有必要處理這個多余的A嗎，當然沒有必要，這樣垂直方向整體上又可以減少約25%的時間，耗時只有75ms左右了，實現(xiàn)了約100%的提速。
第五種嘗試：算法穩(wěn)定性的考慮和最終的妥協(xié)
在上一篇文章中，我們提到了由于float類型的精度問題，當模糊的半徑較大時，算法的結果會出現(xiàn)很大的瑕疵，一種方式就是用double類型來解決，還有一種方式就是可以用Deriche濾波器來解決，為了完美解決這個問題，我還是恨著頭皮用SSE實現(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é)省內存，我們也使用了類似上述第三種和第四重嘗試的方式，但是考慮到Deriche的特殊性(主要是這里)，他還是需要一份中間內存的，為了效率和內存，我們再次以犧牲精度為準備，中間使用了一份和圖像一樣的字節(jié)數(shù)據(jù)內存。
由于計算量較原先的高斯有所增加，這里最后的優(yōu)化完成的耗時約為100ms。
第六：多線程
在水平方向計算時，各行之間的計算時獨立的，因此是可以并行處理的，但是垂直方向由于是前后依賴的，是無法并行的。比如用openmp做2個線程的并行，大概速度能提高到(高斯)到60ms，但是這個東西在不是本文這里的重點。
第七：比較
同標準的高斯濾波相比，Deriche濾波器由于其特性，無法支持In-Place操作，也就是說Src和Dest不能相同，如果一定要相同，就只有通過一個中間內存來過渡了，而標準高斯是可以的。第二就是高斯是可以不占用太多額外的內存就可以實現(xiàn)的，而Deriche需要一份同樣大小的內存。第三就是標準高斯速度還是要快一點。第四Deriche濾波器的精度在float類型時精度要比標準高斯高。綜合選擇，我覺得還是以后選擇Deriche代替標準的高斯模糊。
總結：有心就有成績
同opencv的cvsmooth相比，同樣的機器上同樣的3000*2000大小的彩圖，Ksize我取100時，需要1200多ms，比我這里慢了10倍，但是cvsmooth似乎對ksize參數(shù)敏感，他并不是與核大小無關的，ksize較小時還會很快的，不過除了一些特效外，在很多場合我們其實需要大半徑的高斯的(比如圖像增強、銳化上)。
做完了在回頭看看優(yōu)化的過程，覺得和看書是一個道理，先是越看越厚，通了就像一張薄紙一樣。
最后總結下，就是一件事情，只要你有時間和信心，就能有進步，堅持是勝利的必要條件。
提供一個測試的Demo: http://share.eepw.com.cn/share/download/id/383731
由測試Demo可以測試出，當選擇低內存近似版本或者準確版本時，當半徑較大時，如果連續(xù)的拖動滾動條，圖像會出現(xiàn)閃爍，而如果選擇Deriche時，則圖像變換很平緩，而當半徑特別大時，如果選擇低內存近似版本或者準確版本，則圖像有可能會出現(xiàn)線條或者色塊，只有Deriche濾波的效果是完美的。
高斯模糊的優(yōu)化到此結束，如果有誰有用GPU實現(xiàn)的，還請告訴我下大概的耗時。
拒絕無腦索取代碼。
評論