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