高斯模糊算法的全面優(yōu)化過(guò)程分享（一）
加入技術(shù)交流群
這里的高斯模糊采用的是論文《Recursive implementation of the Gaussian filter》里描述的遞歸算法。
本文引用地址：http://dyxdggzs.com/article/201703/345016.htm仔細觀(guān)察和理解上述公式，在forward過(guò)程中,n是遞增的，因此，如果在進(jìn)行forward之前，把in數據先完整的賦值給w，然后式子(9a)就可以變?yōu)椋?/p>
w[n] = B w[n] + (b1 w[n-1] + b2 w[n-2] + b3 w[n-3]) / b0; ---------> (1a)
在backward過(guò)程中，n是遞減的，因此在backward前，把w的數據完整的拷貝到out中，則式子(9b)變?yōu)椋?/p>
out[n] = B out[n] + (b1 out[n+1] + b2 out[n+2] + b3 out[n+3]) / b0 ; <--------- (1b)
從編程角度看來(lái)，backward中的拷貝是完全沒(méi)有必要的，因此 (1b)可以直接寫(xiě)為：
w[n] = B w[n] + (b1 w[n+1] + b2 w[n+2] + b3 w[n+3]) / b0 ; <--------- (1c)
從速度上考慮，浮點(diǎn)除法很慢，因此，我們重新定義b1 = b1 / b0, b2 = b2 / b0, b3 = b3 / b0, 最終得到我們使用的遞歸公式：
w[n] = B w[n] + b1 w[n-1] + b2 w[n-2] + b3 w[n-3]; ---------> (a)
w[n] = B w[n] + b1 w[n+1] + b2 w[n+2] + b3 w[n+3] ; <--------- (b)
上述公式是一維的高斯模糊計算方法，針對二維的圖像，正確的做法就是先對每個(gè)圖像行進(jìn)行模糊處理得到中間結果，再對中間結果的每列進(jìn)行模糊操作，最終得到二維的模糊結果，當然也可以使用先列后行這樣的做法。
另外注意點(diǎn)就是，邊緣像素的處理，我們看到在公式(a)或者(b)中，w[n]的結果分別依賴(lài)于前三個(gè)或者后三個(gè)元素的值，而對于邊緣位置的值，這些都是不在合理范圍內的，通常的做法是鏡像數據或者重復邊緣數據，實(shí)踐證明，由于是遞歸的算法，起始值的不同會(huì )將結果一直延續下去，因此，不同的方法對邊緣部分的結果還是有一定的影響的，這里我們采用的簡(jiǎn)單的重復邊緣像素值的方式。
由于上面公式中的系數均為浮點(diǎn)類(lèi)型，因此，計算一般也是對浮點(diǎn)進(jìn)行的，也就是說(shuō)一般需要先把圖像數據轉換為浮點(diǎn)，然后進(jìn)行高斯模糊，在將結果轉換為字節類(lèi)型的圖像，因此，上述過(guò)程可以分別用一下幾個(gè)函數完成：
CalcGaussCof // 計算高斯模糊中使用到的系數
ConvertBGR8U2BGRAF // 將字節數據轉換為浮點(diǎn)數據
GaussBlurFromLeftToRight // 水平方向的前向傳播
GaussBlurFromRightToLeft // 水平方向的反向傳播
GaussBlurFromTopToBottom // 垂直方向的前向傳播
GaussBlurFromBottomToTop // 垂直方向的反向傳播
ConvertBGRAF2BGR8U // 將結果轉換為字節數據
我們依次分析之。
CalcGaussCof，這個(gè)很簡(jiǎn)單，就按照論文中提出的計算公式根據用戶(hù)輸入的參數來(lái)計算，當然結合下上面我們提到的除法變乘法的優(yōu)化，注意，為了后續的一些標號的問(wèn)題，我講上述公式(a)和(b)中的系數B寫(xiě)為b0了。
void CalcGaussCof(float Radius, float &B0, float &B1, float &B2, float &B3)
{
float Q, B;
if (Radius >= 2.5)
Q = (double)(0.98711 * Radius - 0.96330); // 對應論文公式11b
else if ((Radius >= 0.5) && (Radius < 2.5))
Q = (double)(3.97156 - 4.14554 * sqrt(1 - 0.26891 * Radius));
else
Q = (double)0.1147705018520355224609375;
B = 1.57825 + 2.44413 * Q + 1.4281 * Q * Q + 0.422205 * Q * Q * Q; // 對應論文公式8c
B1 = 2.44413 * Q + 2.85619 * Q * Q + 1.26661 * Q * Q * Q;
B2 = -1.4281 * Q * Q - 1.26661 * Q * Q * Q;
B3 = 0.422205 * Q * Q * Q;
B0 = 1.0 - (B1 + B2 + B3) / B;
B1 = B1 / B;
B2 = B2 / B;
B3 = B3 / B;
}
由上述代碼可見(jiàn)，B0+B1+B2+B3=1，即是歸一化的系數，這也是算法可以遞歸進(jìn)行的關(guān)鍵之一。
接著(zhù)為了方便中間過(guò)程，我們先將字節數據轉換為浮點(diǎn)數據，這部分代碼也很簡(jiǎn)單：
void ConvertBGR8U2BGRAF(unsigned char *Src, float *Dest, int Width, int Height, int Stride)
{
//#pragma omp parallel for
for (int Y = 0; Y < Height; Y++)
{
unsigned char *LinePS = Src + Y * Stride;
float *LinePD = Dest + Y * Width * 3;
for (int X = 0; X < Width; X++, LinePS += 3, LinePD += 3)
{
LinePD[0] = LinePS[0]; LinePD[1] = LinePS[1]; LinePD[2] = LinePS[2];
}
}
}
大家可以嘗試自己把內部的X循環(huán)再展開(kāi)看看效果。
水平方向的前向傳播按照公式去寫(xiě)也是很簡(jiǎn)單的，但是直接使用公式里面有很多訪(fǎng)問(wèn)數組的過(guò)程(不一定就慢)，我稍微改造下寫(xiě)成如下形式：
void GaussBlurFromLeftToRight(float *Data, int Width, int Height, float B0, float B1, float B2, float B3)
{
//#pragma omp parallel for
for (int Y = 0; Y < Height; Y++)
{
float *LinePD = Data + Y * Width * 3;
float BS1 = LinePD[0], BS2 = LinePD[0], BS3 = LinePD[0]; // 邊緣處使用重復像素的方案
float GS1 = LinePD[1], GS2 = LinePD[1], GS3 = LinePD[1];
float RS1 = LinePD[2], RS2 = LinePD[2], RS3 = LinePD[2];
for (int X = 0; X < Width; X++, LinePD += 3)
{
LinePD[0] = LinePD[0] * B0 + BS1 * B1 + BS2 * B2 + BS3 * B3;
LinePD[1] = LinePD[1] * B0 + GS1 * B1 + GS2 * B2 + GS3 * B3; // 進(jìn)行順向迭代
LinePD[2] = LinePD[2] * B0 + RS1 * B1 + RS2 * B2 + RS3 * B3;
BS3 = BS2, BS2 = BS1, BS1 = LinePD[0];
GS3 = GS2, GS2 = GS1, GS1 = LinePD[1];
RS3 = RS2, RS2 = RS1, RS1 = LinePD[2];
}
}
}
不多描述，請大家把上述代碼和公式(a)對應一下就知道了。
有了GaussBlurFromLeftToRight的參考代碼，GaussBlurFromRightToLeft的代碼就不會(huì )有什么大的困難了，為了避免一些懶人不看文章不思考，這里我不貼這段的代碼。
那么垂直方向上簡(jiǎn)單的做只需要改變下循環(huán)的方向，以及每次指針增加量更改為Width * 3即可。
那么我們來(lái)考慮下垂直方向能不能不這樣處理呢，指針每次增加Width * 3會(huì )造成嚴重的Cache miss，特別是對于寬度稍微大點(diǎn)的圖像，我們仔細觀(guān)察垂直方向，會(huì )發(fā)現依舊可以按照Y -- X這樣的循環(huán)方式也是可以的，代碼如下：
void GaussBlurFromTopToBottom(float *Data, int Width, int Height, float B0, float B1, float B2, float B3)
{
for (int Y = 0; Y < Height; Y++)
{
float *LinePD3 = Data + (Y + 0) * Width * 3;
float *LinePD2 = Data + (Y + 1) * Width * 3;
float *LinePD1 = Data + (Y + 2) * Width * 3;
float *LinePD0 = Data + (Y + 3) * Width * 3;
for (int X = 0; X < Width; X++, LinePD0 += 3, LinePD1 += 3, LinePD2 += 3, LinePD3 += 3)
{
LinePD0[0] = LinePD0[0] * B0 + LinePD1[0] * B1 + LinePD2[0] * B2 + LinePD3[0] * B3;
LinePD0[1] = LinePD0[1] * B0 + LinePD1[1] * B1 + LinePD2[1] * B2 + LinePD3[1] * B3;
LinePD0[2] = LinePD0[2] * B0 + LinePD1[2] * B1 + LinePD2[2] * B2 + LinePD3[2] * B3;
}
}
}
就是說(shuō)我們不是一下子就把列方向的前向傳播進(jìn)行完，而是每次只進(jìn)行一個(gè)像素的傳播，當一行所有像素都進(jìn)行完了列方向的前向傳播后，在切換到下一行，這樣處理就避免了Cache miss出現的情況。
注意到列方向的邊緣位置，為了方便代碼的處理，我們在高度方向上上下分別擴展了3個(gè)行的像素，當進(jìn)行完中間的有效行的行方向前向和反向傳播后，按照前述的重復邊緣像素的方式填充上下那空出的三行數據。
同樣的道理，GaussBlurFromBottomToTop的代碼可由讀者自己補充進(jìn)去。
最后的ConvertBGRAF2BGR8U也很簡(jiǎn)單，就是一個(gè)for循環(huán)：
void ConvertBGRAF2BGR8U(float *Src, unsigned char *Dest, int Width, int Height, int Stride)
{
//#pragma omp parallel for
for (int Y = 0; Y < Height; Y++)
{
float *LinePS = Src + Y * Width * 3;
unsigned char *LinePD = Dest + Y * Stride;
for (int X = 0; X < Width; X++, LinePS += 3, LinePD += 3)
{
LinePD[0] = LinePS[0]; LinePD[1] = LinePS[1]; LinePD[2] = LinePS[2];
}
}
}
在有效的范圍內，上述浮點(diǎn)計算的結果不會(huì )超出byte所能表達的范圍，因此也不需要特別的進(jìn)行Clamp操作。
最后就是一些內存分配和釋放的代碼了，如下所示：
void GaussBlur(unsigned char *Src, unsigned char *Dest, int Width, int Height, int Stride, float Radius)
{
float B0, B1, B2, B3;
float *Buffer = (float *)malloc(Width * (Height + 6) * sizeof(float) * 3);
CalcGaussCof(Radius, B0, B1, B2, B3);
ConvertBGR8U2BGRAF(Src, Buffer + 3 * Width * 3, Width, Height, Stride);
GaussBlurFromLeftToRight(Buffer + 3 * Width * 3, Width, Height, B0, B1, B2, B3);
GaussBlurFromRightToLeft(Buffer + 3 * Width * 3, Width, Height, B0, B1, B2, B3); // 如果啟用多線(xiàn)程，建議把這個(gè)函數寫(xiě)到GaussBlurFromLeftToRight的for X循環(huán)里，因為這樣就可以減少線(xiàn)程并發(fā)時(shí)的阻力
memcpy(Buffer + 0 * Width * 3, Buffer + 3 * Width * 3, Width * 3 * sizeof(float));
memcpy(Buffer + 1 * Width * 3, Buffer + 3 * Width * 3, Width * 3 * sizeof(float));
memcpy(Buffer + 2 * Width * 3, Buffer + 3 * Width * 3, Width * 3 * sizeof(float));
GaussBlurFromTopToBottom(Buffer, Width, Height, B0, B1, B2, B3);
memcpy(Buffer + (Height + 3) * Width * 3, Buffer + (Height + 2) * Width * 3, Width * 3 * sizeof(float));
memcpy(Buffer + (Height + 4) * Width * 3, Buffer + (Height + 2) * Width * 3, Width * 3 * sizeof(float));
memcpy(Buffer + (Height + 5) * Width * 3, Buffer + (Height + 2) * Width * 3, Width * 3 * sizeof(float));
GaussBlurFromBottomToTop(Buffer, Width, Height, B0, B1, B2, B3);
ConvertBGRAF2BGR8U(Buffer + 3 * Width * 3, Dest, Width, Height, Stride);
free(Buffer);
}
正如上所述，分配了Height + 6行的內存區域，主要是為了方便垂直方向的處理，在執行GaussBlurFromTopToBottom之前按照重復邊緣的原則復制3行，然后在GaussBlurFromBottomToTop之前在復制底部邊緣的3行像素。
至此，一個(gè)簡(jiǎn)單而又高效的高斯模糊就基本完成了，代碼數量也不多，也沒(méi)有啥難度，不曉得為什么很多人搞不定。
按照我的測試，上述方式代碼在一臺I5-6300HQ 2.30GHZ的筆記本上針對一副3000*2000的24位圖像的處理時(shí)間大約需要370ms，如果在C++的編譯選項的代碼生成選項里的啟用增強指令集選擇-->流式處理SIMD擴展2(/arch:sse2)，則編譯后的程序大概需要220ms的時(shí)間。
我們嘗試利用系統的一些資源進(jìn)一步提高速度，首先我們想到了SSE優(yōu)化，關(guān)于這方面Intel有一篇官方的文章和代碼：
IIR Gaussian Blur Filter Implementation using Intel? Advanced Vector Extensions [PDF 513KB]
source code: gaussian_blur.cpp [36KB]
我只是簡(jiǎn)單的看了下這篇文章，感覺(jué)他里面用到的計算公式和Deriche濾波器的很像，和本文參考的Recursive implementation 不太一樣，并且其SSE代碼對能處理的圖還有很多限制，對我這里的參考意義不大。
我們先看下核心的計算的SSE優(yōu)化，注意到 GaussBlurFromLeftToRight 的代碼中， 其核心的計算部分是幾個(gè)乘法，但是他只有3個(gè)乘法計算，如果能夠湊成四行，那么就可以充分利用SSE的批量計算功能了，也就是如果能增加一個(gè)通道，使得GaussBlurFromLeftToRight變?yōu)槿缦滦问剑?/p>
void GaussBlurFromLeftToRight(float *Data, int Width, int Height, float B0, float B1, float B2, float B3)
{
//#pragma omp parallel for
for (int Y = 0; Y < Height; Y++)
{
float *LinePD = Data + Y * Width * 4;
float BS1 = LinePD[0], BS2 = LinePD[0], BS3 = LinePD[0]; // 邊緣處使用重復像素的方案
float GS1 = LinePD[1], GS2 = LinePD[1], GS3 = LinePD[1];
float RS1 = LinePD[2], RS2 = LinePD[2], RS3 = LinePD[2];
float AS1 = LinePD[3], AS2 = LinePD[3], AS3 = LinePD[3];
for (int X = 0; X < Width; X++, LinePD += 4)
{
LinePD[0] = LinePD[0] * B0 + BS1 * B1 + BS2 * B2 + BS3 * B3;
LinePD[1] = LinePD[1] * B0 + GS1 * B1 + GS2 * B2 + GS3 * B3; // 進(jìn)行順向迭代
LinePD[2] = LinePD[2] * B0 + RS1 * B1 + RS2 * B2 + RS3 * B3;
LinePD[3] = LinePD[3] * B0 + AS1 * B1 + AS2 * B2 + AS3 * B3;
BS3 = BS2, BS2 = BS1, BS1 = LinePD[0];
GS3 = GS2, GS2 = GS1, GS1 = LinePD[1];
RS3 = RS2, RS2 = RS1, RS1 = LinePD[2];
AS3 = AS2, AS2 = AS1, AS1 = LinePD[3];
}
}
}
則很容易就把上述代碼轉換成SSE可以規范處理的代碼了。
而對于Y方向的代碼，你仔細觀(guān)察會(huì )發(fā)現，無(wú)論是及通道的圖，天然的就可以使用SSE進(jìn)行處理，詳見(jiàn)下面的代碼。
好，我們還是一個(gè)一個(gè)的來(lái)分析：
第一個(gè)函數 CalcGaussCof 無(wú)需進(jìn)行任何的優(yōu)化。
第二個(gè)函數 ConvertBGR8U2BGRAF按照上述分析需要重新寫(xiě)，因為需要增加一個(gè)通道，新的通道的值填0或者其他值都可以，但建議填0，這對有些SSE函數很有用，我把這個(gè)函數的SSE實(shí)現共享一下：
void ConvertBGR8U2BGRAF_SSE(unsigned char *Src, float *Dest, int Width, int Height, int Stride)
{
const int BlockSize = 4;
int Block = (Width - 2) / BlockSize;
__m128i Mask = _mm_setr_epi8(0, 1, 2, -1, 3, 4, 5, -1, 6, 7, 8, -1, 9, 10, 11, -1); // Mask為-1的地方會(huì )自動(dòng)設置數據為0
__m128i Zero = _mm_setzero_si128();
//#pragma omp parallel for
for (int Y = 0; Y < Height; Y++)
{
unsigned char *LinePS = Src + Y * Stride;
float *LinePD = Dest + Y * Width * 4;
int X = 0;
for (; X < Block * BlockSize; X += BlockSize, LinePS += BlockSize * 3, LinePD += BlockSize * 4)
{
__m128i SrcV = _mm_shuffle_epi8(_mm_loadu_si128((const __m128i*)LinePS), Mask); // 提取了16個(gè)字節，但是因為24位的，我們要將其變?yōu)?2位的，所以只能提取出其中的前12個(gè)字節
__m128i Src16L = _mm_unpacklo_epi8(SrcV, Zero);
__m128i Src16H = _mm_unpackhi_epi8(SrcV, Zero);
_mm_store_ps(LinePD + 0, _mm_cvtepi32_ps(_mm_unpacklo_epi16(Src16L, Zero))); // 分配內存時(shí)已經(jīng)是16字節對齊了，然后每行又是4的倍數的浮點(diǎn)數，因此，每行都是16字節對齊的
_mm_store_ps(LinePD + 4, _mm_cvtepi32_ps(_mm_unpackhi_epi16(Src16L, Zero))); // _mm_stream_ps是否快點(diǎn)?
_mm_store_ps(LinePD + 8, _mm_cvtepi32_ps(_mm_unpacklo_epi16(Src16H, Zero)));
_mm_store_ps(LinePD + 12, _mm_cvtepi32_ps(_mm_unpackhi_epi16(Src16H, Zero)));
}
for (; X < Width; X++, LinePS += 3, LinePD += 4)
{
LinePD[0] = LinePS[0]; LinePD[1] = LinePS[1]; LinePD[2] = LinePS[2]; LinePD[3] = 0; // Alpha最好設置為0，雖然在下面的CofB0等SSE常量中通過(guò)設置ALPHA對應的系數為0，可以使得計算結果也為0，但是不是最合理的
}
}
}
稍作解釋?zhuān)琠mm_loadu_si128一次性加載16個(gè)字節的數據到SSE寄存器中，對于24位圖像，16個(gè)字節里包含了5 + 1 / 3個(gè)像素的信息，而我們的目標是把這些數據轉換為4通道的信息，因此，我們只能一次性的提取處16/4=4個(gè)像素的值，這借助于_mm_shuffle_epi8函數和合適的Mask來(lái)實(shí)現，_mm_unpacklo_epi8/_mm_unpackhi_epi8分別提取了SrcV的高8位和低8位的8個(gè)字節數據并將它們轉換為8個(gè)16進(jìn)制數保存到Src16L和Src16H中， 而_mm_unpacklo_epi16/_mm_unpackhi_epi16則進(jìn)一步把16位數據擴展到32位整形數據，最后通過(guò)_mm_cvtepi32_ps函數把整形數據轉換為浮點(diǎn)型。
可能有人注意到了 int Block = (Width - 2) / BlockSize; 這一行代碼，為什么要-2操作呢，這也是我在多次測試發(fā)現程序總是出現內存錯誤時(shí)才意識到的，因為_(kāi)mm_loadu_si128一次性加載了5 + 1 / 3個(gè)像素的信息，當在處理最后一行像素時(shí)(其他行不會(huì ))，如果Block 取Width/BlockSize, 則很有可能訪(fǎng)問(wèn)了超出像素范圍內的內存，而-2不是-1就是因為那個(gè)額外的1/3像素起的作用。
接著(zhù)看水平方向的前向傳播的SSE方案：
void GaussBlurFromLeftToRight_SSE(float *Data, int Width, int Height, float B0, float B1, float B2, float B3)
{
const __m128 CofB0 = _mm_set_ps(0, B0, B0, B0);
const __m128 CofB1 = _mm_set_ps(0, B1, B1, B1);
const __m128 CofB2 = _mm_set_ps(0, B2, B2, B2);
const __m128 CofB3 = _mm_set_ps(0, B3, B3, B3);
//#pragma omp parallel for
for (int Y = 0; Y < Height; Y++)
{
float *LinePD = Data + Y * Width * 4;
__m128 V1 = _mm_set_ps(LinePD[3], LinePD[2], LinePD[1], LinePD[0]);
__m128 V2 = V1, V3 = V1;
for (int X = 0; X < Width; X++, LinePD += 4) // 還有一種寫(xiě)法不需要這種V3/V2/V1遞歸賦值，一次性計算3個(gè)值，詳見(jiàn) D:程序設計正在研究CoreShockFilterConvert里的高斯模糊，但速度上沒(méi)啥區別
{
__m128 V0 = _mm_load_ps(LinePD);
__m128 V01 = _mm_add_ps(_mm_mul_ps(CofB0, V0), _mm_mul_ps(CofB1, V1));
__m128 V23 = _mm_add_ps(_mm_mul_ps(CofB2, V2), _mm_mul_ps(CofB3, V3));
__m128 V = _mm_add_ps(V01, V23);
V3 = V2; V2 = V1; V1 = V;
_mm_store_ps(LinePD, V);
}
}
}
和上面的4通道的GaussBlurFromLeftToRight_SSE比較，你會(huì )發(fā)現基本上語(yǔ)法上沒(méi)有任何的區別，實(shí)在是太簡(jiǎn)單了，注意我沒(méi)有用_mm_storeu_ps，而是直接使用_mm_store_ps，這是因為，第一，分配Data內存時(shí)，我采用了_mm_malloc分配了16字節對齊的內存，而Data每行元素的個(gè)數又都是4的倍數，因此，每行起始位置處的內存也是16字節對齊的，所以直接_mm_store_ps完全是可以的。
同理，在垂直方向的前向傳播的SSE優(yōu)化代碼就更直接了：
void GaussBlurFromTopToBottom_SSE(float *Data, int Width, int Height, float B0, float B1, float B2, float B3)
{
const __m128 CofB0 = _mm_set_ps(0, B0, B0, B0);
const __m128 CofB1 = _mm_set_ps(0, B1, B1, B1);
const __m128 CofB2 = _mm_set_ps(0, B2, B2, B2);
const __m128 CofB3 = _mm_set_ps(0, B3, B3, B3);
for (int Y = 0; Y < Height; Y++)
{
float *LinePS3 = Data + (Y + 0) * Width * 4;
float *LinePS2 = Data + (Y + 1) * Width * 4;
float *LinePS1 = Data + (Y + 2) * Width * 4;
float *LinePS0 = Data + (Y + 3) * Width * 4;
for (int X = 0; X < Width * 4; X += 4)
{
__m128 V3 = _mm_load_ps(LinePS3 + X);
__m128 V2 = _mm_load_ps(LinePS2 + X);
__m128 V1 = _mm_load_ps(LinePS1 + X);
__m128 V0 = _mm_load_ps(LinePS0 + X);
__m128 V01 = _mm_add_ps(_mm_mul_ps(CofB0, V0), _mm_mul_ps(CofB1, V1));
__m128 V23 = _mm_add_ps(_mm_mul_ps(CofB2, V2), _mm_mul_ps(CofB3, V3));
_mm_store_ps(LinePS0 + X, _mm_add_ps(V01, V23));
}
}
}
對上面的代碼不想做任何解釋了。
最有難度的應該算是ConvertBGRAF2BGR8U的SSE版本了，由于某些原因，我不太愿意分享這個(gè)函數的代碼，有興趣的朋友可以參考opencv的有關(guān)實(shí)現。
最后的SSE版本高斯模糊的主要代碼如下：
void GaussBlur_SSE(unsigned char *Src, unsigned char *Dest, int Width, int Height, int Stride, float Radius)
{
float B0, B1, B2, B3;
float *Buffer = (float *)_mm_malloc(Width * (Height + 6) * sizeof(float) * 4, 16);
CalcGaussCof(Radius, B0, B1, B2, B3);
ConvertBGR8U2BGRAF_SSE(Src, Buffer + 3 * Width * 4, Width, Height, Stride);
GaussBlurFromLeftToRight_SSE(Buffer + 3 * Width * 4, Width, Height, B0, B1, B2, B3); // 在SSE版本中，這兩個(gè)函數占用的時(shí)間比下面兩個(gè)要多,不過(guò)C語(yǔ)言版本也是一樣的
GaussBlurFromRightToLeft_SSE(Buffer + 3 * Width * 4, Width, Height, B0, B1, B2, B3); // 如果啟用多線(xiàn)程，建議把這個(gè)函數寫(xiě)到GaussBlurFromLeftToRight的for X循環(huán)里，因為這樣就可以減少線(xiàn)程并發(fā)時(shí)的阻力
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));
GaussBlurFromTopToBottom_SSE(Buffer, Width, Height, B0, B1, B2, B3);
memcpy(Buffer + (Height + 3) * Width * 4, Buffer + (Height + 2) * Width * 4, Width * 4 * sizeof(float));
memcpy(Buffer + (Height + 4) * Width * 4, Buffer + (Height + 2) * Width * 4, Width * 4 * sizeof(float));
memcpy(Buffer + (Height + 5) * Width * 4, Buffer + (Height + 2) * Width * 4, Width * 4 * sizeof(float));
GaussBlurFromBottomToTop_SSE(Buffer, Width, Height, B0, B1, B2, B3);
ConvertBGRAF2BGR8U_SSE(Buffer + 3 * Width * 4, Dest, Width, Height, Stride);
_mm_free(Buffer);
}
對于同樣的3000*2000的彩色圖像，SSE版本的代碼耗時(shí)只有145ms的耗時(shí)，相對于普通的C代碼有約2.5倍左右的提速，這也情有可原，因為我們在執行SSE的代碼時(shí)時(shí)多處理了一個(gè)通道的計算量的，但是同編譯器自身的SSE優(yōu)化220ms，只有1.5倍的提速了，這說(shuō)明編譯器的SSE優(yōu)化能力還是相當強的。
進(jìn)一步的優(yōu)化就是我上面的注釋掉的opemp相關(guān)代碼，在ConvertBGR8U2BGRAF / GaussBlurFromLeftToRight / GaussBlurFromRightToLeft / ConvertBGRAF2BGR8U 4個(gè)函數中，直接加上簡(jiǎn)單的#pragma omp parallel for就可以了，至于GaussBlurFromTopToBottom_SSE、 GaussBlurFromBottomToTop_SSE則由于上下行之間數據的相關(guān)性，是無(wú)法實(shí)現并行計算的，但是測試表示他們的耗時(shí)要比水平方向的少很多。
比如我們指定openmp使用2個(gè)線(xiàn)程，在上述機器上(四核的)，純C版本能優(yōu)化到252ms，而純SSE的只能提高到100ms左右，編譯器自身的SSE優(yōu)化速度大約是150ms，基本上還是保持同一個(gè)級別的比例。
對于灰度圖像，很可惜，上述的水平方向上的優(yōu)化方式就無(wú)論為力了，一種方式就是把4行灰度像素混洗成一行，但是這個(gè)操作不太方便用SSE實(shí)現，另外一種就是把水平方向的數據先轉置，然后利用垂直方向的SSE優(yōu)化代碼處理，完成在轉置回去，最后對轉置的數據再次進(jìn)行垂直方向SSE優(yōu)化，當然轉置的過(guò)程是可以借助于SSE的代碼實(shí)現的，但是需要額外的一份內存，速度上可能和普通的C相比就不會(huì )有那么多的提升了，這個(gè)待有時(shí)間了再去測試。
前前后后寫(xiě)這個(gè)大概也花了半個(gè)月的時(shí)間，寫(xiě)完了整個(gè)流程，在倒過(guò)來(lái)看，真的是非常的簡(jiǎn)單，有的時(shí)候就是這樣。
本文并沒(méi)有提供完整的可以直接執行的全部代碼，需者自勤，提供一段C#的調用工程供有興趣的朋友測試和比對(未使用OPENMP版本)。
http://share.eepw.com.cn/share/download/id/383730
后記：后來(lái)測試發(fā)現，當半徑參數較大時(shí)，無(wú)論是C版本還是SSE版本都會(huì )出現一些不規則的瑕疵，感覺(jué)像是溢出了，后來(lái)發(fā)現主要是當半徑變大后，B0參數變得很小，以至于用float類(lèi)型的數據來(lái)處理遞歸已經(jīng)無(wú)法保證足夠的精度了，解決的辦法是使用double類(lèi)型，這對于C語(yǔ)言版本來(lái)說(shuō)是秒秒的事情，而對于SSE版本則是較大的災難，double時(shí)換成AVX可能改動(dòng)量不大，但是AVX的普及度以及.....，不過(guò)，一般情況下，半徑不大于75時(shí)結果都還是正確的，這對于大部分的應用來(lái)說(shuō)是足夠的，半徑75時(shí)，整個(gè)圖像已經(jīng)差不多沒(méi)有任何的細節了，再大，區別也不大了。
解決上述問(wèn)題一個(gè)可行的方案就是使用Deriche濾波器，我用該濾波器的float版本對大半徑是不會(huì )出現問(wèn)題的，并且也有相關(guān)SSE參考代碼。
評論