卷積其實是圖像處理中最基本的操作，我們常見的一些算法比如：均值模糊、高斯模糊、銳化、Sobel、拉普拉斯、prewitt邊緣檢測等等一些和領域相關的算法，都可以通過卷積算法實現。只不過由於這些算法的卷積矩陣的特殊性，一般不會直接實現它，而是通過一些優化的手段讓計算量變小。但是有些情況下卷積矩陣的元素值無甚規律或者有特殊要求，無法通過常規手段優化，這個時候只能通過原始的方式實現。因此，如何快速的實現圖像的任意卷積矩陣操作也有必要做適當的研究。
目前，通過友人共享或自己搜索找到的一片關於任意核算法優化的文章有： Reshuf?ing: A Fast Algorithm for Filtering with Arbitrary Kernels，改文章稱能夠提高原始程序速度的40%左右，但是原始的程序是如何寫的也還不明白。
在matlab中有幾個函數都與圖像卷積有關，比如imfilter就可以實現卷積，或者 conv2也行，他們的速度都是相當快的，比如3000*3000的灰度圖，卷積矩陣大小為15*15，在I5的CPU上運行時間只要170ms左右，相當的給力。
在Celery的博客中，也提到了他的優化後的conv2和matlab相當甚至快於matlab，詳見http://blog.csdn.net/celerychen2009/article/details/38852105。
由於matlab的代碼中使用到了IPL庫進行加速，目前我寫的Conv2函數還無法做到和其相當，對於任何核速度約為matlab的一半。
簡單的記錄下我在做卷積過程中用到的優化吧。
原始的卷積的實現需要四重循環，簡單的表達如下：
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;
}
}
當卷積矩陣較大時，計算量將會很大，而且由於程序中的內存訪問很頻繁，cache miss現象比較嚴重，因此效率極為低下。
我的優化方法主要包括以下幾個方面：
一：使用SSE進行乘法計算，由於SSE可以一次性進行4個單精度浮點數的計算，因此可以有明顯的速度提升。
二：通過適當的處理方式，對每個取樣點周邊的卷積矩陣內的元素進行集中，使得每移動一個像素點不會需要從內存中進行大量的搜索工作。
具體來說實現過程如下：
1、為了使用SSE的優勢，首先將卷積矩陣進行調整，調整卷積矩陣一行的元素個數，使其為不小於原始值的4的整數倍，並且讓新的卷積矩陣的內存布局符合SSE相關函數的16字節對齊的要求。
實現代碼如下：
float *Conv16 = (float *)_mm_malloc(PadConvLine * ConvH * sizeof(float), 16); // 保存16字節對齊的卷積矩陣，以方便使用SSE
for(Y = 0; Y < ConvH; Y++)
{
memcpy (Conv16 + Y * PadConvLine, Conv->Data.F + Y * ConvW , ConvW * sizeof(float)); // 復制卷積矩陣的數據
memset(Conv16 + Y * PadConvLine + ConvW, 0, (PadConvLine - ConvW) * sizeof(float)); // 把冗余部分的卷積數據設置為0
}
其中PadConvLine = Pad4(ConvW) 以及Pad4的原型為： #define Pad4(bits) (((bits) + 3) / 4 * 4)；
注意_mm_malloc函數分配的內存中的值是隨機值，對於擴展的部分一定要填充0，否則就會破壞卷積的結果。
那麼如果我們也同時獲得了需要被卷積的部分數據的話（卷積核肯定和卷積矩陣一樣大小，且也應該是16字節對齊的），可以用如下的SSE的代碼進行乘法計算：
float MultiplySSE(float *Kernel, float *Conv, int Length)
{
int Block;
const float *Data; // 將SSE變量上的多個數值合並時所用指針.
float Sum = 0;
if (Length > 16) // 可以進行四次SSE計算，測試表明，這個還是快些的
{
const int BlockWidth = 4 * 4; // 塊寬. SSE寄存器能一次處理4個float，然後循環展開4次.
Block = Length / BlockWidth; // 塊數.
float *KernelP = Kernel, *ConvP = Conv; // SSE批量處理時所用的指針.
__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單精浮點緊縮加法
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; // 剩余數量.
}
if (Length != 0)
{
const int BlockWidth = 4; // 程序已經保證了數量必然是4的倍數
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;
}
當卷積矩陣（擴充後）的元素數量大於16時，我們采用了4路並行的SSE乘法實現，我在I3的CPU上測試時，2路SSE和4路SSE已經沒有啥大的區別了，而在I5的CPU上則4路還是有較為明顯的提高，因此采用4路SSE同時運行。當然1路SSE肯定還是比2路慢。另外，如果元素的數量少於16或者大於16但不能被16整除，那麼余下的部分由於先前的擴充，剩余元素數量也肯定是4的倍數，因此可以用單路的SSE實現。 這也是編碼上的技巧。
2、前面提到了需要被卷積的部分數據，這部分如何快速的獲取呢。觀察最原始的4重循環，其內部的2重即為獲取需要被卷積的部分，但是這裡其實有很多問題。第一：由於卷積取樣時必然有部分取樣點的坐標在原始圖像的有效范圍外，因此必須進行判斷，耗時。第二：同樣為了使用SSE,也必須把取樣的數據放在和擴充的卷積矩陣一樣大小的內存中。這裡我先貼出我的代碼在進行解釋具體的實現：
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; // 注意核中心那個元素不用擴展，比如核的寬度為3，則只要左右各擴展一個像素就可以了
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); // 得到擴展後的數據，可以提速和方便編程，但是多占用一份內存
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，則設置為1
float *Conv16 = (float *)_mm_malloc(PadConvLine * ConvH * sizeof(float), 16); // 保存16字節對齊的卷積矩陣，以方便使用SSE
float *Kernel = (float *)_mm_malloc(PadConvLine * ExpHeight * sizeof(float), 16); // 保存16字節對齊的卷積核矩陣，以方便使用SSE
for(Y = 0; Y < ConvH; Y++)
{
memcpy (Conv16 + Y * PadConvLine, Conv->Data.F + Y * ConvW , ConvW * sizeof(float)); // 復制卷積矩陣的數據
memset(Conv16 + Y * PadConvLine + ConvW, 0, (PadConvLine - ConvW) * sizeof(float)); // 把冗余部分的卷積數據設置為0
}
for (Y = 0; Y < ExpHeight; Y++)
{
IndexE = Y * ExpStride;
CurKer = Kernel + Y * PadConvLine; // 計算第一列所有像素將要取樣的卷積核數據
for (X = 0; X < ConvW; X++)
{
CurKer[X] = PtExp[IndexE++];
}
}
for (X = 0 ; X < Width ; X ++)
{
if (X != 0) // 如果不是第一列，需要更新卷積核的數據
{
memcpy(Kernel, Kernel + 1, (PadConvLine * ExpHeight - 1) * sizeof(float)); // 往前移動一個數據
IndexK = ConvW - 1 ;
IndexE = IndexK + X;
for (Y = 0; Y < ExpHeight; Y++)
{
Kernel[IndexK] = PtExp[IndexE]; // 只要刷新下一個元素
IndexK += PadConvLine;
IndexE += ExpStride;
}
}
CurKer = Kernel; IndexD = X;
for (Y = 0; Y < Height; Y ++) // 沿列的方向進行更新
{
PtDest[IndexD] = Clamp((int)( MultiplySSE(Conv16, CurKer, Length) * Inv + 0.5)); // 直接把函數放在這裡也沒有啥提速的，注意改函數不會被內聯的
CurKer += PadConvLine;
IndexD += Stride;
}
}
_mm_free(Conv16);
_mm_free(Kernel);
FreeImage(Expand);
return IS_RET_OK;
}
}
對於第一個問題，解決的方式很簡答，即用空間換時間，新建一副（Width + ConvW - 1, Height + ConvH -1)大小的圖像，然後四周的ConvW及ConvH的像素用邊緣的值或者邊緣鏡像的值填充，正中間的則用原來的圖復制過來，這樣操作後進行取樣時不再原圖取樣，而在這福擴展的圖中取樣，就避免了坐標判斷等if語句的跳轉耗時了，上GetPadImage即實現了改功能。
第二個問題則需要有一定的實現技巧，我們分配一塊PadConvLine * (Height + ConvH - 1) 大小的內存，然後計算原圖第一列像素串聯起來的需要卷積的部分的數據，這一部分代碼如上述44-52行所示。有了這樣的數據，如果需要計算第一列的卷積結果，則很簡單了，每跳過一列則把被卷積的數據起點增加PadConvLine個元素，在調用上述MultiplySSE函數獲得卷積結果。接著則計算第二列像素的卷積值，此時需要整體更新這一列像素串聯起來的需要被卷積的數據，更新也很簡單，就是把原來的數據整體向左移動一個像素，這個可以用memcpy快速實現，然後在填充入新進來的那個元素，就ok了，接著就是再次調用MultiplySSE函數，如此重復下去。
經過編碼測試，對於3000*3000的灰度圖，15*15的核在I5的CPU上的測試平均結果為360ms,比matlab的慢了一半。
最後說明一點，很多人都說用FFT可以快速的實現卷積，並且是O(1)的，我比較同意後半句，但是前面半句是絕對的有問題的，至少在核小於50*50時，FFT實現的卷積不會比直接實現塊。要知道FFT的計算量其實是很大的。
****************************基本上我不提供源代碼，但是我會盡量用文字把對應的算法描述清楚或提供參考文檔************************
*************************************因為靠自己的努力和實踐寫出來的效果才真正是自己的東西，人一定要靠自己*******************
****************************作者： laviewpbt 時間： 2014.11.27 聯系QQ: 33184777 轉載請保留本行信息**********************
