使用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:

1. 每一個half-wrap的第一個thread有aligned在適當的記憶體位置 (該位置必需是每個thread存取的資料大小的16倍)
   
2. "連續"地讀取記憶體

為了滿足第一項條件，以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，這些程式都還需要再加以改進。