Image convolution是影像處理上常用的方法之一，可應用在影像模糊、影像強化、擷取邊緣...等等。在數學上，convolution代表兩個函數(ex: f and g)的重疊運算，結果會產生第3個函數(h) :
f * g = h
應用在影像處理時，函數f為原影像，另一個函數g則為filter，而h可視為原影像被convolution後的結果。 convolution在2維discrete domain的運算過程，可用下圖來表示:
要計算source image上每一個像素的convolution結果，需要以該像素為中心取出與filter相同大小的區域(window)，然後該區域的每一個像素與filter相對應的位置做點對點的相乘，最後計算相乘結果的總和，即為該像素convolution的結果。
在影像上計算每一像素的convolution時，會遇到邊界像素的window超出影像範圍的問題，超出的部分就稱為apron，如下圖所示:
而apron的尺寸會等於filter半徑的寬度。當apron所處的位置已超出影像範圍，必需特別給定此apron的值，才能與filter進行運算，而最常見的方法是補0，或者重複邊界像素的值。 進行2維convolution時，一般需要2維影像與2維的filter進行運算，然而某些2維的filter可由2個1維的filter來組成，當此filter可分離時，2維的計算可等效於連續使用2個1維的filter進行運算，如此一來，最大的好處是減少了計算量，例如當filter尺寸為n X m時，原本每個像素需要進行n*m次乘法，分離filter後，只需要n+m次乘法。
測試平台
- OS: Win XP x64
- CPU: Core 2 Extreme X9650 3.0GHz (4 cores)
- VGA: Geforce GTX280
- nVidia CUDA 3.1
使用CPU計算Image convolution
先使用CPU計算的程式來測試，使用3張不同尺寸的圖片，filter使用高斯濾波器，設定filter kernel的尺寸為17*17(半徑為8)，並分離成2個一維的filter來計算。設定超出影像範圍的apron其值為0。效能測試中，皆進行50次運算，再計算平均花費時間。
CPU:
Parallelization with OpenMP:
CUDA程式實作
在nvidia CUDA SDK中，有convolution的範例程式convolutionSeparable及convolutionTexture，這兩支程式做了最佳化而有很好的效能，但它們並沒有實際接受影像輸入，而是用亂數產生float type的資料，模擬一個channel的影像。接下來，我們先不考慮最佳化，而用最直接的方法實作，再一步步改進程式，最後比較各種方法的差異。
The most naive approach
最基本的方法是將影像資料送到device上的global memory，再讓每個thread在計算convolution的kernel中存取它。整張影像將被分成許多block，如下圖所示，而每一個thread負責計算一個像素的結果，因此在這個應用中，將沒有任何idle的thread。設定每一個block的size為16 X 16。
效能上比平行化的CPU程式快1.29倍。
convolution kernel的code如下:
1: __global__ void my2DConvKernel(float *d_Result, float *d_Data, int dataW, int dataH)2: {3: // original image based coordinate4: int y = blockIdx.y * blockDim.y + threadIdx.y;5: int x = blockIdx.x * blockDim.x + threadIdx.x;6:7: int BiasY = y - KERNEL_RADIUS;8: int BiasX = x - KERNEL_RADIUS;9:10: float sum = 0;11: for(int j = 0; j < KERNEL_LENGTH; ++j)12: {13: //out of image range14: if (BiasY + j < 0 || BiasY + j >= dataH)15: continue;16:17: for(int i = 0; i < KERNEL_LENGTH; ++i)18: {19: //out of image range20: if (BiasX + i < 0 || BiasX + i >= dataW)21: continue;22:23: sum += d_Data[(BiasY + j) * dataW + BiasX + i] *24: c_Kernel[KERNEL_LENGTH * j + i];25: }26: }27:28: d_Result[y * dataW + x] = sum;29: }30:
其中24行的c_Kernel為filter，在呼叫kernel前就先把它傳至constant memory中，雖然constant memory是唯讀，但具有快取，適合將唯讀的資料放於此，加速存取。
2D convolution using shared memory
上一個程式，只把影像資料存到global memory，每讀一次原圖的值都要存取一次，然而存取global memory所需要的時間(latency)是很長的，它也沒有快取，很容易讓運算的瓶頸卡在資料傳輸上。改善的方法是利用shared memory，存取的速度上會快很多，甚至與存取暫存器的速度相同，但它有大小限制，每一個multiprocessor只有16KB的shared memory。接下來，我們把一個block送進shared memory中，其中含要處理的像素以及apron區域，如下圖所示:
在此仍然讓一個thread只處理一個像素的結果，在load資料至shared memory時，每個thread負責一個像素的傳輸，而在計算階段，只有負責傳輸白色區域的thread需要計算出結果，負責傳輸apron的thread會idle，因此雖然利用到了shared memory，但卻減少能同時計算的thread數量，有得也有失。然而在此測試中，filter半徑為8時，若白色區域尺寸維持16 X 16，會使block尺寸達到32 X 32，如此一來會超過每個block容許512個thread的上限(顯卡硬體限制，因不同顯卡而異)，所以在此設定白色區域尺寸為4 X 4。
效能上比上一個方法"慢"1.96倍。
convolution kernel的code如下:
1: __global__ void NaiveSharedMemoryGPUKernel(float *d_Result, float *d_Data, int dataW, int dataH)2: {3: __shared__ float data[KERNEL_RADIUS * 2 + ACTIVE_T_W][KERNEL_RADIUS * 2 + ACTIVE_T_H];4:5: // global mem address for this thread6: const int gLoc = blockIdx.y * ACTIVE_T_H * dataW +7: blockIdx.x * ACTIVE_T_W +8: threadIdx.y * dataW + threadIdx.x -9: KERNEL_RADIUS * dataW - KERNEL_RADIUS;10:11: // original image based coordinate12: int x = blockIdx.x * ACTIVE_T_W - KERNEL_RADIUS + threadIdx.x;13: int y = blockIdx.y * ACTIVE_T_H - KERNEL_RADIUS + threadIdx.y;14:15: if (x < 0 || x > dataW - 1 || y < 0 || y > dataH - 1 )16: data[threadIdx.x][threadIdx.y] = 0;17: else18: data[threadIdx.x][threadIdx.y] = d_Data[gLoc];19:20: __syncthreads();21:22: float sum = 0;23: if(threadIdx.x >= KERNEL_RADIUS &&24: threadIdx.x < KERNEL_RADIUS + ACTIVE_T_W &&25: threadIdx.y >= KERNEL_RADIUS &&26: threadIdx.y < KERNEL_RADIUS + ACTIVE_T_H)27: {28: // row wise29: for (int i = -KERNEL_RADIUS; i <= KERNEL_RADIUS; ++i)30: {31: // col wise32: for (int j = -KERNEL_RADIUS; j <= KERNEL_RADIUS; ++j)33: {34: sum += data[threadIdx.x + i][threadIdx.y + j] *35: c_Kernel[KERNEL_RADIUS + i] *36: c_Kernel[KERNEL_RADIUS + j];37: }38: }39:40: d_Result[gLoc] = sum;41: }42: }
以上的ACTIVE_T_W與ACTIVE_T_H為白色區的長與寬，皆設定為4。第3行利用__shared__來宣告一個使用shared memory的變數，data。第20行的__syncthreads()為CUDA函式，讓同一個block內的thread在此同步，才能繼續執行，在此因為要把data的資料都填完才能做後續的計算，因此需要做同步的動作。
Separable convolution using shared memory
此方法先將2D的filter分解成2個一維的filter，再接連對影像做X方向及Y方向的convolution。當只考慮X方向的convolution時，每個block只需包含X方向的apron，如下圖:
如此一來，不用考慮Y方向的apron，即可增加每個block中用來計算的thread。另外，也可減少影像上同一區域被重複讀進shared memory的次數，原本2D convolution中，同一區域最多會被讀取9次，separable convolution方法中最多只會被讀取6次。
效能上比上一個方法快4.63倍。
convolution kernel的code如下:
1: __global__ void mySeparableRowKernel(float *d_Result, float *d_Data, int dataW, int dataH)2: {3: __shared__ float data[KERNEL_RADIUS * 2 + ACTIVE_T_W][ACTIVE_T_H];4:5: // global mem address for this thread6: const int gLoc = blockIdx.y * ACTIVE_T_H * dataW +7: blockIdx.x * ACTIVE_T_W * dataW +8: threadIdx.x -9: KERNEL_RADIUS + threadIdx.y;10:11: // original image based coordinate12: int x = blockIdx.x * ACTIVE_T_W - KERNEL_RADIUS + threadIdx.x;13:14: if (x < 0 || x > dataW - 1)15: data[threadIdx.x][threadIdx.y] = 0;16: else17: data[threadIdx.x][threadIdx.y] = d_Data[gLoc];18:19: __syncthreads();20:21: float sum = 0;22: if(threadIdx.x >= KERNEL_RADIUS &&23: threadIdx.x < KERNEL_RADIUS + ACTIVE_T_W)24: {25: // row wise26: for (int i = -KERNEL_RADIUS; i <= KERNEL_RADIUS; ++i)27: {28: sum += data[threadIdx.x + i][threadIdx.y] *29: c_Kernel[KERNEL_RADIUS + i];30: }31:32: d_Result[gLoc] = sum;33: }34: }
1: __global__ void mySeparableColKernel(float *d_Result, float *d_Data, int dataW, int dataH)2: {3: __shared__ float data[ACTIVE_T_W][KERNEL_RADIUS * 2 + ACTIVE_T_H];4:5: // global mem address for this thread6: const int gLoc = blockIdx.y * ACTIVE_T_H * dataW +7: blockIdx.x * ACTIVE_T_W +8: threadIdx.y * dataW +9: threadIdx.x -10: KERNEL_RADIUS * dataW;11:12: // original image based coordinate13: int y = blockIdx.y * ACTIVE_T_H - KERNEL_RADIUS + threadIdx.y;14:15: if (y < 0 || y > dataH - 1)16: data[threadIdx.x][threadIdx.y] = 0;17: else18: data[threadIdx.x][threadIdx.y] = d_Data[gLoc];19:20: __syncthreads();21:22: float sum = 0;23: if(threadIdx.y >= KERNEL_RADIUS &&24: threadIdx.y < KERNEL_RADIUS + ACTIVE_T_H)25: {26: // col wise27: for (int j = -KERNEL_RADIUS; j <= KERNEL_RADIUS; ++j)28: {29: sum += data[threadIdx.x][threadIdx.y + j] *30: c_Kernel[KERNEL_RADIUS + j];31: }32:33: d_Result[gLoc] = sum;34: }35: }
nvidia convolution application: convolutionSeparable
nvidia SDK中的範例程式convolutionSeparable，使用了一些最佳化技巧來提升程式執行效率，在此介紹一下。首先在記憶體存取上，存取global memory資料時有符合"coalescence"，滿足的條件有2:
- 每一個half-wrap的第一個thread有aligned在適當的記憶體位置 (該位置必需是每個thread存取的資料大小的16倍)
- "連續"地讀取記憶體
為了滿足第一項條件，以row的convolution來說，每個block除了包含需計算結果的區域及apron區外，還多讀取了一些額外的資料，如下圖的綠色區域，其目的就是要讓讀取記憶體的初始位置能align在適當的位置。
從範例程式中得知，這個位置的所在是從白色區域往前算的第16個資料的倍數，當appron寬度小於16時，就往前多載入16個資料，若寬度大於16，則載入16的倍數的資料，如此一來讀取的範圍既可涵蓋appron又滿足位置的aligned。 接下來談第二項條件，連續讀取記憶體的部分，先瞭解GPU的運作方式，當一個thread去存取記憶體並等待資料時，GPU會先切換到下個thread並執行，因此存取記憶體的順序應該是:
所以下圖讓thread 0存取紅色部分，thread 1與2各別存取黃色與綠色部分，並沒有滿足連續存取記憶體。
而是要如下圖所示，才能滿足連續存取記憶體。
nvidia的程式透過memory coalescence的最佳化，並將某些迴圈以loop unrolling的方式展開，來提高執行效率，結果如下:
效能上比上一個方法快1.24倍。
測試結果
下圖為本篇所有測試的效能比較圖。
最後結果，nvidia的方法比沒有平行化的CPU程式還快上12.63倍，也比平行化的CPU程式快上3.78倍。其中2D convolution使用shared memory那一項，因為每個block只有4X4個thread能執行計算，比其他的測試少，因此效能甚至比直接存取global memory還不好。但整體而言，使用shared memory仍對效能有所幫助。
這次的測試程式(包含nvidia code)，其實只能在一些理想條件下運作，例如影像只能有單channel、影像的長與寬必需是block尺寸的倍數、filter的尺寸不能太大以免單一block超過512個thread…等等。所以若要能處理各類型的影像與不同大小filter進行convolution，這些程式都還需要再加以改進。
