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

沈 飛 陳榮昌 肖體喬

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

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

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

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

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

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

1.1 CUDA編程模型

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

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

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

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

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

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

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

1.2.1 正弦圖

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

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

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

1.2.2 濾波處理

由于濾波反投影算法會(huì)帶來(lái)星狀偽跡，故需進(jìn)行濾波消除偽跡。本文采用Sheep Logan濾波函數(shù)進(jìn)行濾波處理，用濾波函數(shù)重建圖像，其振蕩響應(yīng)減小，對(duì)于含有噪聲的投影數(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ù)進(jìn)行卷積操作，即：

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

1.2.3 反投影

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

圖3 圖像重建時(shí)所用的坐標(biāo)系Fig.3 Coordinate system of image reconstruction.

濾波反投影公式為：

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

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

2 實(shí)驗(yàn)結(jié)果與分析

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

開(kāi)發(fā)工具為 Visual Studio2005，主機(jī)端使用C++進(jìn)行編程，設(shè)備端采用CUDA標(biāo)準(zhǔn)的類(lèi)C語(yǔ)言進(jìn)行編程。

我們?cè)贑PU中進(jìn)行CT重建時(shí)，采用單精度浮點(diǎn)型數(shù)據(jù)進(jìn)行保存，而GPU支持雙精度浮點(diǎn)運(yùn)算，其重建結(jié)果符合實(shí)驗(yàn)要求(圖4)。

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

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

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

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

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

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

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

3 結(jié)語(yǔ)

本文研究了通過(guò) GPU加速使用濾波反投影原理的CT重建方法，得到了比較理想的結(jié)果。實(shí)驗(yàn)表明，支持單精度的GPU重建圖像質(zhì)量完全符合要求，其單張切片的重建可得到150多倍的加速比。GPU強(qiáng)大的計(jì)算能力可用于快速CT重構(gòu)，是一種性?xún)r(jià)比較高的重構(gòu)算法。若采用多GPU進(jìn)行CT重建的加速，則加速比更高，甚至直接用于實(shí)時(shí) 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實(shí)現(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高性能運(yùn)算之CUDA[M]. 北京: 中國(guó)水利水電出版社，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 馬車(chē)平. GPU多重紋理加速三維ct圖像重建[J]. 計(jì)算機(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