Reconstruction Method Using GPU (CUDA)

疊代式的影像重建演算法 SART 與 OS-SART 計算過程需要非常龐大的 計算量，所以利用 CPU 在實作上非常耗時；而且當投影資料(Sinogram)角度 和射線數目愈多，則需要花費的時間就愈龐大。

因此我們藉由圖形處理器(GPGPU)平行處理與適合大量資料的運算能 力，來加快疊代重建過程需花費的時間。分析疊代重建演算法的步驟中，可 以利用平行處理特性運算包含兩個部分[11]：估計投影資料(Re-projection) 與反投影(Back projection)。

3.4.1 Parallel Re-projection Computing Method

在 Re-projection 步驟中，利用圖形處理器的架構下，我們平行計算 角 度中整排射線的取樣過程，可以同時得到該 角度下全部射線的估計投影值，

概念如圖 3-9。

圖 3-9：平行計算應用於 Re-projection 概念圖

利用 CUDA 開發環境實作 Re-projection 平行化，必須將其參考矩陣、

投影角度 、單位向量(unit vector)與射線進入點座標由 CPU 端複製到圖形 處理器端。我們配置一個執行緒(Thread)代表該角度的每一條射線，並且讓 每個 Thread 執行 GPGPU 中「kernel」程式。於程式中計算每條射線的進入 點(Entry point)，並且每間隔單位向量長度 1 同時對參考矩陣取樣；利用雙 線性內插的計算方式，將每一個取樣點的估計投影值依序累加至 ，即為每 條射線所得到的估計投影值。存於一維陣列中，最後將此一維陣列的資料複 製回 CPU 端進行運算。

由於每條射線在影像區域中必須取樣非常多次，因此在圖形處理器記憶 體的管理中，我們使用 Texture Memory 來存放物體影像矩陣。CUDA 中的 Texture Memory 的特性為，限制只能讀取不能寫入、支援 Spatial Locality 存 取方式下快速存取且頻寬較大、以及可以自動計算貼圖中座標位置的內插取 值，相較於利用 Global Memory 來存放物體影像矩陣，改善讀取記憶體過程 的時間延遲。

3.4.2 Parallel Back Projection Computing Method

在原先 Back projection 的實作中，演算方法為將該射線的修正值反向投 影回物體的影像矩陣中。而我們在圖形處理器的架構下，將此步驟的演算方 法設計為適合平行處理的計算方式。

3.4.1 Re-projection 可以同時得到某 角度中所有射線的修正值。因此我 們可以設計平行處理演算方法為，以影像矩陣中的每個元素為出發點，反向 對修正值進行內插取值，並填回矩陣元素中成為影像像素的 Intensity 值；平 行計算概念如圖 3-10。

圖 3-10：平行計算應用於 Back projection 概念圖

利用 CUDA 開發環境實作 Back projection 的平行化，首先將修正值、

影像矩陣與投影角度 由 CPU 端複製到圖形處理器端。修正值矩陣的記憶

體配置上，我們選擇存放於 Share Memory 中，讓同一個 block 中的所有 thread 快速的進行取值。

由於 CUDA 中 1 個 block 只能包含 512 條 thread，通常二維影像矩陣遠 大於這個大小限制；因此我們必須依照影像的尺寸來配置 block 與 thread 的 大小。例如：若影像矩陣尺寸為 256x256，則我們可以令二維 block 維度 16，

即索引值為(0,0),(0,1),...,(16,16)，總個數為 256 個；而每個 block 裡面包含二 維 thead 數 目 為 (256/16) (256/16) 個 ， 即 索 引 為 (0,0),(0,1),...,((256/16),(256/16))。利用這個方式，我們可以計算 block 與 thread 的索引編號，來代表影像矩陣中的每個元素；概念如圖 3-11 所示。

圖 3-11：CUDA 配置二維矩陣索引概念圖

GPGPU 中 Back projection 的「kernel」程式，程式中利用上述方法得到 (a,b)，用來代表影像矩陣中每個元素的索引值。將原先的座標系統依照投影 角度 旋轉，得到原座標(a,b)對於新座標系統的 x 座標， cos

sin 。利用此x 座標，即可對該角度的修正值進行內插取值。平行計 算影像矩陣中所有元素(a,b)，同時得到該像素的 Intensity 值；再將此計算完 畢的影像矩陣複製到 CPU 端，即完成該 角度的 Back projection 步驟。

利 用 圖 形 處 理 器 平 行 計 算 Re-projection 與 Back projection ， 不 用 ray-by-ray 對參考矩陣與影像矩陣作內插取值計算，大幅降低演算方法的時 間複雜度。因此我們可以簡化 OS-SART 影像重建演算法如下：

OS-SART algorithm implementation:

假設角度 0 ~π之間，CT Scanner 對物體影像進行 p 次投影，且每個投影 角度皆包含q 條射線；subset_num 為 OS-SART 中 subset 的個數。

Input : sinogram[p][q];

Array : estimate_projection[q], correction_array[q], reference_image[q][q];

Output : CT_image[q][q];

for i = 1 to p

do estimate_projection[q] Å Re-projection using GPU(*reference_image);

for j = 1 to q

do correction_array[q] Å sinogram[i][j] - estimate_projection[j];

do CT_imageÅ Back-projection using GPU(*correction_array);

if i = subset_num

then reference_image Å CT_image