基于GPU并行計算實現(xiàn)快速顯微CT重構(gòu)

沈 飛 陳榮昌 肖體喬

1 (中國科學(xué)院上海應(yīng)用物理研究所 上海 201800)

2 (中國科學(xué)院研究生院 北京 100049)

近些年來，三維CT在醫(yī)學(xué)和工業(yè)領(lǐng)域得到廣泛應(yīng)用。隨著探測器分辨率的改進，CT投影圖像采樣分辨率也不斷提高，對高質(zhì)量輸出容積體期望值的不斷提升，使存儲容量和數(shù)據(jù)的處理量越來越大，嚴重影響CT重建速度。目前，上海光源成像線站采用串行的 CT重建程序并通過 CPU完成計算，計算數(shù)據(jù)容積為2 048×2 048×1 500，就需40多小時，迫切需要進行CT重建的加速。

圖像處理器(GPU, graphics processing unit)是高度并行化的流式處理器，并行計算能力強，運算速度快，在通用計算(GPGPU)領(lǐng)域有極大優(yōu)勢。文獻[1]利用 GPU 實現(xiàn) FDK 算法(The algorithm of Feldkamp, Davis and Kress)的加速，最高可達15倍左右的加速比。Scherl[2]利用 CUDA在 NVDIA Geforce8800實現(xiàn)了FDK算法，取得5倍左右的加速比。

使用FDK算法和迭代算法進行CT重建的GPU加速研究，已取得較佳效果[3?5]。而基于濾波反投影算法的GPU加速的研究并不多見，該法的主要優(yōu)點是重建速度快。本文采用濾波反投影(FBP)算法實現(xiàn)CT重建，由于反投影時每個位置的像素值重建相互獨立，可同時進行重建計算，故可將其并行化設(shè)計。

1 GPU并行化實現(xiàn)CT重建

1.1 CUDA編程模型

CUDA編程模型[6]將CPU作為主機(Host)，GPU作為協(xié)處理器(co-processor)或設(shè)備(Device)。一個系統(tǒng)中可存在一個主機和若干個設(shè)備，采用單指令、多線程執(zhí)行模式(SIMT)。

模型中CPU和GPU協(xié)同工作，各司其職。CPU負責(zé)進行邏輯性強的事務(wù)處理和串行計算，GPU則專注于執(zhí)行高度線程化的并行處理任務(wù)。CPU和GPU各自擁有相互獨立的存儲器地址空間：主機端的內(nèi)存和設(shè)備端的顯存。

一旦確定了程序中的并行部分，即可考慮把這部分計算工作交給GPU。將運行在GPU上的并行計算函數(shù)稱為kernel(內(nèi)核函數(shù))。一個kernel函數(shù)并不是一個完整的程序，而是整個CUDA程序中一個可被并行執(zhí)行的步驟(圖1)。一個完整的CUDA程序，是由一系列設(shè)備端kernel函數(shù)并行步驟和主機端串行處理步驟共同組成。這些處理步驟會按照程序中的相應(yīng)語句順序依次執(zhí)行，滿足順序一致性。

圖1 CUDA編程模型Fig.1 Programming model of CUDA.

進行GPU代碼設(shè)計時，對全局現(xiàn)存的訪問是否符合合并訪問條件是對 CUDA程序性能影響最明顯的因素之一，可能產(chǎn)生一個量級的差距。因為，若符合一定的訪問條件，只需一次傳輸即可完成一個Half-warp所有線程的數(shù)據(jù)訪問[6]。

1.2 CT重建kernel的設(shè)計

濾波反投影重建是一種解析方法，先獲取某一層所有角度的投影數(shù)據(jù)得到正弦圖，再進行濾波降噪處理，最后進行反投影重建得到斷層每個位置的像素值，其優(yōu)點是重建速度快[7]。

1.2.1 正弦圖

濾波反投影前，需獲取重建切片的正弦圖，由于直接獲取的正弦圖含有背景噪聲，需去除之，我們采用的方法如式(1)：

式中，Pi(x,y)為第 i層切片(x,y)處的像素值，D(x,i)為無光照時第 i層切片(x,i)處的像素值，F(xiàn)(x,i)為有光照無樣品時第i層切片(x,i)處的像素值。

由式(1)，每個坐標處的計算互不相干，可并行執(zhí)行，則可在CPU代碼中獲取正弦圖數(shù)據(jù)，并編寫kernel函數(shù)去除背景噪聲，可見kernel執(zhí)行過程的數(shù)據(jù)訪問符合合并訪問的條件。

1.2.2 濾波處理

由于濾波反投影算法會帶來星狀偽跡，故需進行濾波消除偽跡。本文采用Sheep Logan濾波函數(shù)進行濾波處理，用濾波函數(shù)重建圖像，其振蕩響應(yīng)減小，對于含有噪聲的投影數(shù)據(jù)，重建效果較好[8]，其函數(shù)式為：

式中，n為取值范圍(?width/2,width/2)，width為圖像的寬度；d取1。其圖形如圖2所示。

圖2 S_L濾波函數(shù)的離散表示Fig.2 Discrete representation of S_L filter function.

將圖像函數(shù)和濾波函數(shù)進行卷積操作，即：

由于CUDA有FFT和IFFT的API，本文采用傅里葉變換方式求卷積，濾波函數(shù)的取值范圍和圖像的寬度相一致。傅里葉變換后求乘積即位置相對應(yīng)的點相乘，可將其并行化處理，而且兩個函數(shù)的寬度一樣，所以kernel在執(zhí)行時其數(shù)據(jù)訪問符合合并訪問的條件。

1.2.3 反投影

CT重建的計算量主要在反投影部分，反投影部分的加速效果直接影響重建過程的加速比。反投影過程即為累加求和的過程，在進行累加求和前需進行射束計算和線性內(nèi)插。圖像重建時指定切片的像素為N×N，如圖(3)所示。

圖3 圖像重建時所用的坐標系Fig.3 Coordinate system of image reconstruction.

濾波反投影公式為：

由式 4，圖像重建時每個點的重建互不干擾，且每個點重建的方法相同，只是重建時所需的數(shù)據(jù)不一樣。所以，單個切片的每個坐標處的重建是可以并行執(zhí)行的。

由于重建時數(shù)據(jù)的訪問不符合合并訪問的條件，帶有隨機性，所以采用將數(shù)據(jù)映射到紋理存儲器。紋理存儲器通過緩存利用數(shù)據(jù)的局部性提高效率，使用紋理時未嚴格遵守合并訪問條件，也能獲得較高的帶寬[6,8]。

2 實驗結(jié)果與分析

實驗平臺的配置為：中央處理器Intel i5-650，圖形處理器NVIDIA GTX295，內(nèi)存Kingston DDR3 4GB，操作系統(tǒng)Windows XP Professional SP2 X86，顯卡驅(qū)動NV driver 190.20 + CUDA 2.3 SDK。

開發(fā)工具為 Visual Studio2005，主機端使用C++進行編程，設(shè)備端采用CUDA標準的類C語言進行編程。

我們在CPU中進行CT重建時，采用單精度浮點型數(shù)據(jù)進行保存，而GPU支持雙精度浮點運算，其重建結(jié)果符合實驗要求(圖4)。

圖4 CPU和GPU重建結(jié)果的比較Fig.4 Comparison of reconstruction between CPU and GPU.

由于GTX295擁有30個流多處理器(SM)，而每個SM中擁有8個SIMD的流處理器(SP)，總計算核心達 240個。實驗中對尺寸為 704×301×1 200的數(shù)據(jù)進行了重建的加速(704和301是圖像的高度與寬度，1 200是投影數(shù))，其中對單張切片的正弦圖處理、濾波、反投影部分采用GPU加速后的時間進行了比較，結(jié)果見表1。

由表1，除去數(shù)據(jù)傳輸?shù)难訒r，使用GPU加速，其正弦圖處理和濾波部分可得到200倍以上的加速比，而在反投影部分，將數(shù)據(jù)映射到紋理存儲器中，其加速比增加一個量級。本實驗還對尺寸為1 016× 1 016×1 620和2 048×2 048×1 500的物體進行重建加速的比較，計算部分的加速比都在150倍以上。

表1 某一層數(shù)據(jù)重建的時間(物體尺寸：704×301, 投影數(shù)：1 200)Table 1 Reconstruction time of a slice(object size, 704×301; projection number, 1 200).

隨著數(shù)據(jù)量的不斷增加，CPU讀取硬盤數(shù)據(jù)及CPU內(nèi)存與GPU顯存之間數(shù)據(jù)傳輸?shù)臅r間不斷增加。對不同容積的數(shù)據(jù)，讀取硬盤數(shù)據(jù)和 CPU與 GPU間數(shù)據(jù)傳輸并將數(shù)據(jù)寫存到硬盤的總時間見圖5(a)；圖像重建的加速比見圖5(b)。

由圖5(b)，當(dāng)重建的數(shù)據(jù)大小達到一定程度時，其整體的加速比會降低，這是因為當(dāng)數(shù)據(jù)量增至一定程度，其計算部分的加速比達到飽和，而此時CPU的內(nèi)存和GPU的顯存之間以及CPU內(nèi)存和硬盤之間數(shù)據(jù)傳輸?shù)臅r間隨著數(shù)據(jù)尺寸而增大，從而使整體的加速比降低。

圖5 不同容積的數(shù)據(jù)傳輸時間(a)和加速比(b)的比較Fig.5 Comparison of different data sizes in transmission time (a) and speedup ratio (b).

3 結(jié)語

本文研究了通過 GPU加速使用濾波反投影原理的CT重建方法，得到了比較理想的結(jié)果。實驗表明，支持單精度的GPU重建圖像質(zhì)量完全符合要求，其單張切片的重建可得到150多倍的加速比。GPU強大的計算能力可用于快速CT重構(gòu)，是一種性價比較高的重構(gòu)算法。若采用多GPU進行CT重建的加速，則加速比更高，甚至直接用于實時 CT重構(gòu)。

1 FANG Xu. Real-time 3D computed tomographic reconstruction using commodity graphics hardware[J]. Phys Med Biol, 2007, 52(12): 3405–3419

2 Scherl H. Fast GPU-based CT reconstruction using the Common Unified Device Architecture (CUDA)[Z]. IEEE Nucl Sci Symp, Med Imaging Conf, 2007, (6): 4464–4466

3 Mueller K, Yagel R. Rapid 3-D cone-beam reconstruction with the Simultaneous Algebrmc Reconstruction Technique(SART) using 2-D texture mapping hardware[J]. IEEE Trans Med Imag, 2000, 19(12): 1227–1237

4 Sharp G C. GPU-based streaming architectures for fast cone-beam CT image reconstruction and demons deformable registration[J]. Phys Med Biol, 2007, 52(19): 5771–5783

5 胡君杰. 利用 GPU實現(xiàn)三維圖像重建方法研究[J]. 世界科技研究與發(fā)展，2008, 30(1): 24–26

HU Junjie. The research of 3-D image reconstruction in GPU[J]. Res Dev world, 2008, 30(1): 24–26

6 張 舒. GPU高性能運算之CUDA[M]. 北京: 中國水利水電出版社，2009: 16–168

ZHANG Shu. High performance computing of GPU–CUDA[M]. Beijing: China Water Power Press, 2009: 16–168

7 莊天戈. CT原理與算法[M]. 上海: 上海交通大學(xué)出版社，1992: 30–60

ZHANG Tiange. The principle and algorithm of CT[M]. Shanghai: Shanghai Jiaotong University Press, 1992: 30–60

8 馬車平. GPU多重紋理加速三維ct圖像重建[J]. 計算機工程與應(yīng)用，2008, 44(7): 227–230

MA Cheping. Accelerating the 3-D CT image reconstruction with multiple texture of GPU[J]. Comput Eng Appl, 2008, 44(7): 227–230