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