卷積其實(shí)是圖像處理中最基本的操作，我們常見(jiàn)的一些算法比如：均值模糊、高斯模糊、銳化、Sobel、拉普拉斯、prewitt邊緣檢測(cè)等等一些和領(lǐng)域相關(guān)的算法，都可以通過(guò)卷積算法實(shí)現(xiàn)。只不過(guò)由于這些算法的卷積矩陣的特殊性，一般不會(huì)直接實(shí)現(xiàn)它，而是通過(guò)一些優(yōu)化的手段讓計(jì)算量變小。但是有些情況下卷積矩陣的元素值無(wú)甚規(guī)律或者有特殊要求，無(wú)法通過(guò)常規(guī)手段優(yōu)化，這個(gè)時(shí)候只能通過(guò)原始的方式實(shí)現(xiàn)。因此，如何快速的實(shí)現(xiàn)圖像的任意卷積矩陣操作也有必要做適當(dāng)?shù)难芯俊?/p>

? ? ? 目前，通過(guò)友人共享或自己搜索找到的一片關(guān)于任意核算法優(yōu)化的文章有：?Reshuf?ing: A Fast Algorithm for Filtering?with Arbitrary Kernels，改文章稱能夠提高原始程序速度的40%左右，但是原始的程序是如何寫的也還不明白。

? ? ? 在matlab中有幾個(gè)函數(shù)都與圖像卷積有關(guān)，比如imfilter就可以實(shí)現(xiàn)卷積，或者 conv2也行，他們的速度都是相當(dāng)快的，比如3000*3000的灰度圖，卷積矩陣大小為15*15，在I5的CPU上運(yùn)行時(shí)間只要170ms左右，相當(dāng)?shù)慕o力。

? ? ? 在Celery的博客中，也提到了他的優(yōu)化后的conv2和matlab相當(dāng)甚至快于matlab，詳見(jiàn)http://blog.csdn.net/celerychen2009/article/details/38852105。

? ? ? 由于matlab的代碼中使用到了IPL庫(kù)進(jìn)行加速，目前我寫的Conv2函數(shù)還無(wú)法做到和其相當(dāng)，對(duì)于任何核速度約為matlab的一半。

? ? ? 簡(jiǎn)單的記錄下我在做卷積過(guò)程中用到的優(yōu)化吧。

? ? ? 原始的卷積的實(shí)現(xiàn)需要四重循環(huán)，簡(jiǎn)單的表達(dá)如下：

for (Y = 0; Y < Height; Y++) {for (X = 0; X < Width; X++){Index = .....;Sum = 0;for (XX = 0; XX < ConvW; XX++){for (YY = 0; YY < ConvH; YY++){Index1 = ..... ;Index2 = ..... ;Sum += Conv[Index1] * Pixel[Index2];}}Dest[Index] = Sum / Weight;} }

　　當(dāng)卷積矩陣較大時(shí)，計(jì)算量將會(huì)很大，而且由于程序中的內(nèi)存訪問(wèn)很頻繁，cache miss現(xiàn)象比較嚴(yán)重，因此效率極為低下。

? ? ?我的優(yōu)化方法主要包括以下幾個(gè)方面：?

? ? ?一：使用SSE進(jìn)行乘法計(jì)算，由于SSE可以一次性進(jìn)行4個(gè)單精度浮點(diǎn)數(shù)的計(jì)算，因此可以有明顯的速度提升。

? ? ?二：通過(guò)適當(dāng)?shù)奶幚矸绞?#xff0c;對(duì)每個(gè)取樣點(diǎn)周邊的卷積矩陣內(nèi)的元素進(jìn)行集中，使得每移動(dòng)一個(gè)像素點(diǎn)不會(huì)需要從內(nèi)存中進(jìn)行大量的搜索工作。

? ? ?具體來(lái)說(shuō)實(shí)現(xiàn)過(guò)程如下：

? ? ? ? ? ?1、為了使用SSE的優(yōu)勢(shì)，首先將卷積矩陣進(jìn)行調(diào)整，調(diào)整卷積矩陣一行的元素個(gè)數(shù)，使其為不小于原始值的4的整數(shù)倍，并且讓新的卷積矩陣的內(nèi)存布局符合SSE相關(guān)函數(shù)的16字節(jié)對(duì)齊的要求。

　　實(shí)現(xiàn)代碼如下：

float *Conv16 = (float *)_mm_malloc(PadConvLine * ConvH * sizeof(float), 16); // 保存16字節(jié)對(duì)齊的卷積矩陣，以方便使用SSEfor(Y = 0; Y < ConvH; Y++) {memcpy (Conv16 + Y * PadConvLine, Conv->Data.F + Y * ConvW , ConvW * sizeof(float)); // 復(fù)制卷積矩陣的數(shù)據(jù)memset(Conv16 + Y * PadConvLine + ConvW, 0, (PadConvLine - ConvW) * sizeof(float)); // 把冗余部分的卷積數(shù)據(jù)設(shè)置為0 }

? ? 其中PadConvLine = Pad4(ConvW) 以及Pad4的原型為：?#define Pad4(bits) (((bits) + 3) / 4 * 4)；

? ??注意_mm_malloc函數(shù)分配的內(nèi)存中的值是隨機(jī)值，對(duì)于擴(kuò)展的部分一定要填充0，否則就會(huì)破壞卷積的結(jié)果。

? ? 那么如果我們也同時(shí)獲得了需要被卷積的部分?jǐn)?shù)據(jù)的話（卷積核肯定和卷積矩陣一樣大小，且也應(yīng)該是16字節(jié)對(duì)齊的），可以用如下的SSE的代碼進(jìn)行乘法計(jì)算：

float MultiplySSE(float *Kernel, float *Conv, int Length) {int Block; const float *Data; // 將SSE變量上的多個(gè)數(shù)值合并時(shí)所用指針.float Sum = 0;if (Length > 16) // 可以進(jìn)行四次SSE計(jì)算，測(cè)試表明，這個(gè)還是快些的 {const int BlockWidth = 4 * 4; // 塊寬. SSE寄存器能一次處理4個(gè)float，然后循環(huán)展開(kāi)4次.Block = Length / BlockWidth; // 塊數(shù). float *KernelP = Kernel, *ConvP = Conv; // SSE批量處理時(shí)所用的指針.__m128 Sum0 = _mm_setzero_ps(); // 求和變量。SSE賦初值0__m128 Sum1 = _mm_setzero_ps();__m128 Sum2 = _mm_setzero_ps();__m128 Sum3 = _mm_setzero_ps();for(int I = 0; I < Block; I++){Sum0 = _mm_add_ps(Sum0, _mm_mul_ps(_mm_load_ps(KernelP), _mm_load_ps(ConvP))); // SSE單精浮點(diǎn)緊縮加法Sum1 = _mm_add_ps(Sum1, _mm_mul_ps(_mm_load_ps(KernelP + 4), _mm_load_ps(ConvP + 4)));Sum2 = _mm_add_ps(Sum2, _mm_mul_ps(_mm_load_ps(KernelP + 8), _mm_load_ps(ConvP + 8)));Sum3 = _mm_add_ps(Sum3, _mm_mul_ps(_mm_load_ps(KernelP + 12), _mm_load_ps(ConvP + 12)));KernelP += BlockWidth;ConvP += BlockWidth;}Sum0 = _mm_add_ps(Sum0, Sum1); // 兩兩合并(0~1).Sum2 = _mm_add_ps(Sum2, Sum3); // 兩兩合并(2~3).Sum0 = _mm_add_ps(Sum0, Sum2); // 兩兩合并(0~2).Data = (const float *)&Sum0;Sum = Data[0] + Data[1] + Data[2] + Data[3];Length = Length - Block * BlockWidth; // 剩余數(shù)量.}if (Length != 0){const int BlockWidth = 4; // 程序已經(jīng)保證了數(shù)量必然是4的倍數(shù) Block = Length / BlockWidth; float *KernelP = Kernel, *ConvP = Conv; __m128 Sum0 = _mm_setzero_ps(); for(int I = 0; I < Block; I++){Sum0 = _mm_add_ps(Sum0, _mm_mul_ps(_mm_load_ps(KernelP), _mm_load_ps(ConvP))); KernelP += BlockWidth;ConvP += BlockWidth;}Data = (const float *)&Sum0;Sum += Data[0] + Data[1] + Data[2] + Data[3];}return Sum; }

　　當(dāng)卷積矩陣（擴(kuò)充后）的元素?cái)?shù)量大于16時(shí)，我們采用了4路并行的SSE乘法實(shí)現(xiàn)，我在I3的CPU上測(cè)試時(shí)，2路SSE和4路SSE已經(jīng)沒(méi)有啥大的區(qū)別了，而在I5的CPU上則4路還是有較為明顯的提高，因此采用4路SSE同時(shí)運(yùn)行。當(dāng)然1路SSE肯定還是比2路慢。另外，如果元素的數(shù)量少于16或者大于16但不能被16整除，那么余下的部分由于先前的擴(kuò)充，剩余元素?cái)?shù)量也肯定是4的倍數(shù)，因此可以用單路的SSE實(shí)現(xiàn)。 這也是編碼上的技巧。

? ? ?2、前面提到了需要被卷積的部分?jǐn)?shù)據(jù)，這部分如何快速的獲取呢。觀察最原始的4重循環(huán)，其內(nèi)部的2重即為獲取需要被卷積的部分，但是這里其實(shí)有很多問(wèn)題。第一：由于卷積取樣時(shí)必然有部分取樣點(diǎn)的坐標(biāo)在原始圖像的有效范圍外，因此必須進(jìn)行判斷，耗時(shí)。第二：同樣為了使用SSE,也必須把取樣的數(shù)據(jù)放在和擴(kuò)充的卷積矩陣一樣大小的內(nèi)存中。這里我先貼出我的代碼在進(jìn)行解釋具體的實(shí)現(xiàn)：

IS_RET __stdcall Conv2(TImage *Src, TMatrix *Conv, TImage *Dest, EdgeMode Edge) {if (Src == NULL || Dest == NULL || Conv == NULL) return IS_RET_ERR_PARA;if (Src->Width != Dest->Width || Src->Height != Dest->Height || Src->BitCount != Dest->BitCount || Src->Stride != Dest->Stride) return IS_RET_ERR_PARA;if (Src->Scan0 == NULL || Dest->Scan0 == NULL || Conv->Data.F == NULL) return IS_RET_ERR_MEM;if (Conv->Width < 1 || Conv->Height < 1) return IS_RET_ERR_PARA;int Width = Src->Width, Height = Src->Height, Stride = Src->Stride;int ConvW = Conv->Width, ConvH = Conv->Height;unsigned char *PtSrc = Src->Scan0, *PtDest = Dest->Scan0;if (Src->BitCount == 24){}else{int Left = ConvW / 2, Top = ConvH / 2, Right = ConvW - Left - 1, Bottom = ConvH - Top - 1, ExpHeight = Height + ConvH - 1; // 注意核中心那個(gè)元素不用擴(kuò)展，比如核的寬度為3，則只要左右各擴(kuò)展一個(gè)像素就可以了int PadConvLine = Pad4(ConvW), Length = PadConvLine * ConvH;int X, Y, IndexD, IndexE, IndexK, ExpStride;float *CurKer, Inv, Sum = 0;unsigned char *PtExp, *PtDest;TImage *Expand;IS_RET Ret = GetPadImage(Src, &Expand, Left, Right, Top, Bottom, Edge); // 得到擴(kuò)展后的數(shù)據(jù)，可以提速和方便編程，但是多占用一份內(nèi)存if (Ret != IS_RET_OK) return Ret;PtExp = Expand->Scan0; PtDest = Dest->Scan0; ExpStride = Expand->Stride;for (X = 0; X < ConvH * ConvW; X ++) Sum += Conv->Data.F[X];Inv = (Sum == 0 ? 1: 1 / Sum); // 如果卷積舉證的和為0，則設(shè)置為1float *Conv16 = (float *)_mm_malloc(PadConvLine * ConvH * sizeof(float), 16); // 保存16字節(jié)對(duì)齊的卷積矩陣，以方便使用SSEfloat *Kernel = (float *)_mm_malloc(PadConvLine * ExpHeight * sizeof(float), 16); // 保存16字節(jié)對(duì)齊的卷積核矩陣，以方便使用SSEfor(Y = 0; Y < ConvH; Y++) {memcpy (Conv16 + Y * PadConvLine, Conv->Data.F + Y * ConvW , ConvW * sizeof(float)); // 復(fù)制卷積矩陣的數(shù)據(jù)memset(Conv16 + Y * PadConvLine + ConvW, 0, (PadConvLine - ConvW) * sizeof(float)); // 把冗余部分的卷積數(shù)據(jù)設(shè)置為0}for (Y = 0; Y < ExpHeight; Y++){IndexE = Y * ExpStride;CurKer = Kernel + Y * PadConvLine; // 計(jì)算第一列所有像素將要取樣的卷積核數(shù)據(jù)for (X = 0; X < ConvW; X++){CurKer[X] = PtExp[IndexE++];}}for (X = 0 ; X < Width ; X ++){if (X != 0) // 如果不是第一列，需要更新卷積核的數(shù)據(jù){memcpy(Kernel, Kernel + 1, (PadConvLine * ExpHeight - 1) * sizeof(float)); // 往前移動(dòng)一個(gè)數(shù)據(jù)IndexK = ConvW - 1 ;IndexE = IndexK + X;for (Y = 0; Y < ExpHeight; Y++){Kernel[IndexK] = PtExp[IndexE]; // 只要刷新下一個(gè)元素IndexK += PadConvLine;IndexE += ExpStride;}}CurKer = Kernel; IndexD = X;for (Y = 0; Y < Height; Y ++) // 沿列的方向進(jìn)行更新{PtDest[IndexD] = Clamp((int)( MultiplySSE(Conv16, CurKer, Length) * Inv + 0.5)); // 直接把函數(shù)放在這里也沒(méi)有啥提速的，注意改函數(shù)不會(huì)被內(nèi)聯(lián)的CurKer += PadConvLine;IndexD += Stride;}}_mm_free(Conv16);_mm_free(Kernel);FreeImage(Expand);return IS_RET_OK;} }

? ? ?對(duì)于第一個(gè)問(wèn)題，解決的方式很簡(jiǎn)答，即用空間換時(shí)間，新建一副（Width + ConvW - 1, Height + ConvH -1)大小的圖像，然后四周的ConvW及ConvH的像素用邊緣的值或者邊緣鏡像的值填充，正中間的則用原來(lái)的圖復(fù)制過(guò)來(lái)，這樣操作后進(jìn)行取樣時(shí)不再原圖取樣，而在這福擴(kuò)展的圖中取樣，就避免了坐標(biāo)判斷等if語(yǔ)句的跳轉(zhuǎn)耗時(shí)了，上GetPadImage即實(shí)現(xiàn)了改功能。

? ? ?第二個(gè)問(wèn)題則需要有一定的實(shí)現(xiàn)技巧，我們分配一塊PadConvLine * (Height + ConvH - 1) 大小的內(nèi)存，然后計(jì)算原圖第一列像素串聯(lián)起來(lái)的需要卷積的部分的數(shù)據(jù)，這一部分代碼如上述44-52行所示。有了這樣的數(shù)據(jù)，如果需要計(jì)算第一列的卷積結(jié)果，則很簡(jiǎn)單了，每跳過(guò)一列則把被卷積的數(shù)據(jù)起點(diǎn)增加PadConvLine個(gè)元素，在調(diào)用上述MultiplySSE函數(shù)獲得卷積結(jié)果。接著則計(jì)算第二列像素的卷積值，此時(shí)需要整體更新這一列像素串聯(lián)起來(lái)的需要被卷積的數(shù)據(jù)，更新也很簡(jiǎn)單，就是把原來(lái)的數(shù)據(jù)整體向左移動(dòng)一個(gè)像素，這個(gè)可以用memcpy快速實(shí)現(xiàn)，然后在填充入新進(jìn)來(lái)的那個(gè)元素，就ok了，接著就是再次調(diào)用MultiplySSE函數(shù)，如此重復(fù)下去。

? ? ?經(jīng)過(guò)編碼測(cè)試，對(duì)于3000*3000的灰度圖，15*15的核在I5的CPU上的測(cè)試平均結(jié)果為360ms,比matlab的慢了一半。

? ? ?最后說(shuō)明一點(diǎn)，很多人都說(shuō)用FFT可以快速的實(shí)現(xiàn)卷積，并且是O(1)的，我比較同意后半句，但是前面半句是絕對(duì)的有問(wèn)題的，至少在核小于50*50時(shí)，FFT實(shí)現(xiàn)的卷積不會(huì)比直接實(shí)現(xiàn)塊。要知道FFT的計(jì)算量其實(shí)是很大的。

? ? ?http://www.cnblogs.com/Imageshop/p/4126753.html

http://blog.csdn.net/celerychen2009/article/details/38852105

總結(jié)

以上是生活随笔為你收集整理的图像处理之卷积---任意卷积核的快速实现的全部?jī)?nèi)容，希望文章能夠幫你解決所遇到的問(wèn)題。

如果覺(jué)得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。