使用CUDA實作image convolution

Image convolution是影像處理上常用的方法之一,可應用在影像模糊、影像強化、擷取邊緣…等等。在數學上,convolution代表兩個函數(ex: f and g)的重疊運算,結果會產生第3個函數(h) :

f * g = h

應用在影像處理時,函數f為原影像,另一個函數g則為filter,而h可視為原影像被convolution後的結果。convolution在2維discrete domain的運算過程,可用下圖來表示:image

要計算source image上每一個像素的convolution結果,需要以該像素為中心取出與filter相同大小的區域(window),然後該區域的每一個像素與filter相對應的位置做點對點的相乘,最後計算相乘結果的總和,即為該像素convolution的結果。

在影像上計算每一像素的convolution時,會遇到邊界像素的window超出影像範圍的問題,超出的部分就稱為apron,如下圖所示:

image

而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:

image

Parallelization with OpenMP:

image

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。

image

效能上比平行化的CPU程式快1.29倍。

image

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區域,如下圖所示:

image

在此仍然讓一個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倍。

image

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,如下圖:

image

如此一來,不用考慮Y方向的apron,即可增加每個block中用來計算的thread。另外,也可減少影像上同一區域被重複讀進shared memory的次數,原本2D convolution中,同一區域最多會被讀取9次,separable convolution方法中最多只會被讀取6次。

效能上比上一個方法快4.63倍。

image

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在適當的位置。

image

從範例程式中得知,這個位置的所在是從白色區域往前算的第16個資料的倍數,當appron寬度小於16時,就往前多載入16個資料,若寬度大於16,則載入16的倍數的資料,如此一來讀取的範圍既可涵蓋appron又滿足位置的aligned。接下來談第二項條件,連續讀取記憶體的部分,先了解GPU的運作方式,當一個thread去存取記憶體並等待資料時,GPU會先切換到下個thread並執行,因此存取記憶體的順序應該是:

image

所以下圖讓thread 0存取紅色部分,thread 1與2各別存取黃色與綠色部分,並沒有滿足連續存取記憶體。

image

而是要如下圖所示,才能滿足連續存取記憶體。

image

nvidia的程式透過memory coalescence的最佳化,並將某些迴圈以loop unrolling的方式展開,來提高執行效率,結果如下:

效能上比上一個方法快1.24倍。

image

測試結果

下圖為本篇所有測試的效能比較圖。

image

最後結果,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,這些程式都還需要再加以改進。

1 thought on “使用CUDA實作image convolution”

  1. 博主,你好!最近在学习基于CUDA的卷积优化方法,你的这篇《使用CUDA實作image convolution》对我很有帮助,但有一个小小的问题,这篇博客中的图解都down掉了,很影响整体的阅读和理解,可否将这篇博文电子档发我一份,以便参考学习之用,多谢楼主?
    Email:hpliu51@163.com

發佈留言

發佈留言必須填寫的電子郵件地址不會公開。