使用CUDA實作image convolution

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 coordinate
4:   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 range
14:     if (BiasY   j < 0 || BiasY   j >= dataH) 
15:       continue;
16: 
17:     for(int i = 0; i < KERNEL_LENGTH;   i)
18:     {
19:       //out of image range
20:       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 thread
6:   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 coordinate
12:   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:   else
18:     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 wise
29:     for (int i = -KERNEL_RADIUS; i <= KERNEL_RADIUS;   i) 
30:     {
31:       // col wise
32:       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 thread
6:   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 coordinate
12:   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:   else
17:     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 wise
26:     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 thread
6:   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 coordinate
13:   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:   else
18:     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 wise
27:     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，這些程式都還需要再加以改進。