劉俞辰

(中北大學(xué)信息與通信工程學(xué)院，山西 太原 030051)

近年來(lái)，工業(yè)CT 技術(shù)高速發(fā)展，應(yīng)用也越來(lái)越廣泛。三維錐束CT 具有掃描速度快、射線(xiàn)利用效率高和重建圖像軸向分辨率高［1］等特點(diǎn)，受到了越來(lái)越多的關(guān)注。錐束CT圖像重建的數(shù)據(jù)量大，公式比較復(fù)雜，計(jì)算量較大。傳統(tǒng)的CPU 重建方法，耗時(shí)巨大，無(wú)法滿(mǎn)足工業(yè)上實(shí)時(shí)、快速、準(zhǔn)確的重建要求［2］。2007 年NVIDIA 公司推出了基于CUDA 的GPU，CUDA 不需要借助于圖形學(xué)API，采用了比較容易掌握的類(lèi)C 語(yǔ)言進(jìn)行開(kāi)發(fā)，大大降低了開(kāi)發(fā)難度。Scherl 等［3］使用CUDA 技術(shù)實(shí)現(xiàn)了FDK 算法的加速，取得了比基于圖形學(xué)API 加速算法更好的加速效果。Mueller 等［4］通過(guò)利用緩存優(yōu)化，取得了加速比的進(jìn)一步提高;Papenhausen 等［5］提出的精細(xì)優(yōu)化策略，更好地發(fā)揮了GPU 在通用計(jì)算方面的優(yōu)勢(shì)。

本文使用CUDA 流和優(yōu)化的濾波方法以及合理的線(xiàn)程分配等方法，更大程度發(fā)揮GPU 在并行計(jì)算方面的優(yōu)勢(shì)，提高了FDK 算法的效率。實(shí)驗(yàn)結(jié)果表明，本文提出的優(yōu)化方法在實(shí)際工業(yè)中有廣泛的應(yīng)用前景。

1 基礎(chǔ)知識(shí)

1.1 CUDA 基礎(chǔ)

CUDA 是NVIDIA 公司2007 年推出的一種GPU 架構(gòu)模式，該模式集合了數(shù)據(jù)并行計(jì)算的軟/硬件體系和開(kāi)發(fā)工具?；贑UDA 的編程方式由四層組成:應(yīng)用程序編程接口、硬件驅(qū)動(dòng)程序以及它的運(yùn)行庫(kù)、兩個(gè)高級(jí)數(shù)學(xué)函數(shù)庫(kù)(CUFFT和CUBLAS)。CUDA 還提供了CUDA profiler 用于分析CUDA 程序的性能，可以提供程序中所有核函數(shù)和內(nèi)存拷貝占用的時(shí)間、流多處理器的占用率、活動(dòng)block 數(shù)量、靜態(tài)共享內(nèi)存大小以及全局存儲(chǔ)器的訪(fǎng)問(wèn)情況等信息，對(duì)開(kāi)發(fā)人員分析和優(yōu)化程序有很大的幫助。

1.2 FDK 算法原理

在眾多的三維重建算法中，考慮到運(yùn)算量和工程實(shí)現(xiàn)難度，F(xiàn)DK 算法最為實(shí)用。圖1 為該算法的掃描方式，o－xyz為掃描對(duì)象坐標(biāo)系，o－ts 為射線(xiàn)源旋轉(zhuǎn)坐標(biāo)系，o－tz 為探測(cè)器坐標(biāo)系。S 射線(xiàn)源通過(guò)射線(xiàn)源中心并與探測(cè)器平面垂直，D 為射線(xiàn)到旋轉(zhuǎn)中心O 的距離，β 為射線(xiàn)源旋轉(zhuǎn)的角度。

圖1 FDK 掃描原理圖

FDK 重建算法［6］如下:

其中:f(t，s，z)為被重建函數(shù)，h(·)為濾波函數(shù)，Rβ(p，β)為物體轉(zhuǎn)動(dòng)到β 角度下的二維投影，p1、ξ 為反投影地址。FDK的重建步驟可歸納為:

2)對(duì)二維投影數(shù)據(jù)Rβ(p，β)沿行方向進(jìn)行濾波，濾波函數(shù)為h(·);

3)計(jì)算出投影地址p1、ξ，進(jìn)行反投影重構(gòu)。

從式(1)可知，F(xiàn)DK 算法在濾波和反投影過(guò)程是對(duì)每個(gè)角度下的數(shù)據(jù)單獨(dú)處理，且不同角度下的數(shù)據(jù)是獨(dú)立的，對(duì)于同一角度下的數(shù)據(jù)，不同的重建點(diǎn)也是相互獨(dú)立的。在計(jì)算過(guò)程中，不依賴(lài)相鄰像素值，也不依賴(lài)前一計(jì)算結(jié)果，因此具有很強(qiáng)的并行性，易于用CUDA 實(shí)現(xiàn)。

2 加速方法

2.1 使用GPU 并行化FDK 重建方法

2.1.1 CUDA 流加速

由于錐束CT 重建的數(shù)據(jù)量很大，目前的解決方法都是把重建數(shù)據(jù)分成N 塊，分步傳入GPU 中，數(shù)據(jù)經(jīng)處理后，再傳回CPU。因?yàn)镚PU 的運(yùn)算速度遠(yuǎn)大于數(shù)據(jù)傳輸速度，導(dǎo)致GPU 編譯器大部分時(shí)間處于等待狀態(tài)，浪費(fèi)了大量時(shí)間。如今大部分GPU 中包含一個(gè)或多個(gè)內(nèi)存復(fù)制操作引擎，以及一個(gè)核函數(shù)執(zhí)行引擎，這些引擎彼此獨(dú)立地對(duì)操作進(jìn)行排隊(duì)。因此本文利用CUDA 流在GPU 上同時(shí)執(zhí)行多個(gè)任務(wù)，減少等待狀態(tài)的時(shí)間浪費(fèi)來(lái)加速重建過(guò)程。

CUDA 流表示一個(gè)GPU 操作隊(duì)列，該隊(duì)列中的操作將以指定的順序執(zhí)行。本文使用兩個(gè)流進(jìn)行加速。每個(gè)流都包含三個(gè)操作:將重建數(shù)據(jù)復(fù)制到GPU，執(zhí)行核函數(shù)，將計(jì)算結(jié)果復(fù)制回CPU。如果同時(shí)調(diào)度某個(gè)流的所有操作，那么很容易阻塞另一個(gè)流的復(fù)制操作或者核函數(shù)執(zhí)行，影響GPU運(yùn)行效率。為了解決這個(gè)問(wèn)題，本文不采用深度優(yōu)先方式，而是寬度優(yōu)先方式，即不是首先添加第0 個(gè)流的所有三個(gè)操作，然后繼續(xù)添加第1 個(gè)流的所有三個(gè)操作，而是將這兩個(gè)流之間的操作交叉添加。

圖2 執(zhí)行時(shí)間線(xiàn)

執(zhí)行時(shí)間線(xiàn)如圖2 所示，其中箭頭表示引擎間的依賴(lài)性，空的方框表示一個(gè)流正在等待執(zhí)行某個(gè)操作的時(shí)刻，這個(gè)操作不能與其它流的操作互相重疊。

從圖中可以看出，第0 個(gè)流對(duì)處理后數(shù)據(jù)的復(fù)制操作將不會(huì)阻塞第1 個(gè)流對(duì)處理前數(shù)據(jù)的內(nèi)存復(fù)制操作。這使得GPU 能夠并行地執(zhí)行復(fù)制操作和核函數(shù)，從而使應(yīng)用程序的運(yùn)行速度顯著加快。若執(zhí)行復(fù)制操作的引擎和執(zhí)行核函數(shù)的引擎始終保持運(yùn)行，那么可以在性能上獲得近2 倍的提升。

2.1.2 濾波加速方法

在FDK 算法步驟，對(duì)投影數(shù)據(jù)濾波相當(dāng)于求卷積。CUDA 的CUFFT 函數(shù)庫(kù)中自帶了FFT 函數(shù)，因此可以使用GPU并行運(yùn)算來(lái)加速濾波過(guò)程。

由于CUFFT 庫(kù)中的FFT 函數(shù)只能對(duì)復(fù)數(shù)進(jìn)行傅里葉變換，因此先通過(guò)CUDA 流把一部分投影數(shù)據(jù)傳輸?shù)紾PU 顯存后，再把相鄰兩排的投影數(shù)據(jù)一個(gè)作為復(fù)數(shù)的實(shí)部，另一個(gè)作為復(fù)數(shù)的虛部組成新的一行，每行的投影數(shù)據(jù)從M 個(gè)變成2 ×M 個(gè)，總行數(shù)從N 行變成了N/2 行，這樣就能同時(shí)對(duì)兩行數(shù)據(jù)進(jìn)行濾波。然后調(diào)用CUFFT 庫(kù)中的FFT 函數(shù)對(duì)重排后的數(shù)據(jù)按行進(jìn)行一維快速傅里葉變換，與濾波核函數(shù)卷積計(jì)算完后，再進(jìn)行一次快速傅里葉反變換運(yùn)算。最后將傅里葉反變換得到的投影數(shù)據(jù)重排到對(duì)應(yīng)的兩行中。

2.1.3 反投影加速方法

反投影過(guò)程是FDK 算法中計(jì)算量最大的部分，主要包括各個(gè)角度下每個(gè)像素投影坐標(biāo)的計(jì)算、距離加權(quán)與求和。

GPU 在進(jìn)行并行運(yùn)算時(shí)，會(huì)以block 為單位，把一個(gè)個(gè)的block 分配給SM 進(jìn)行運(yùn)算;而block 中的thread，又會(huì)以warp 為單位，把thread 分組來(lái)計(jì)算。GPU 的計(jì)算效率取決于grid 和block 中線(xiàn)程的分配方式，而線(xiàn)程分配是否恰當(dāng)取決于SM 的占用率的高低。這是因?yàn)镚PU 是通過(guò)線(xiàn)程間的切換來(lái)隱藏長(zhǎng)延時(shí)操作(訪(fǎng)問(wèn)全局存儲(chǔ)器，線(xiàn)程間同步等)的，當(dāng)一個(gè)warp 中的線(xiàn)程進(jìn)行長(zhǎng)延時(shí)操作時(shí)，另一個(gè)活動(dòng)warp 中的線(xiàn)程便會(huì)自動(dòng)開(kāi)始運(yùn)算。但SM 占用率也不是越高越好，因?yàn)镾M 占用率越高則每個(gè)線(xiàn)程占用的warp 資源越少，即每個(gè)線(xiàn)程完成的計(jì)算量越少，而每個(gè)SM 上的最大活動(dòng)warp 數(shù)量是一定的，這就有可能導(dǎo)致所有活動(dòng)warp 中線(xiàn)程同時(shí)進(jìn)入長(zhǎng)延時(shí)操作，便不能充分地隱藏延遲。因此，在SM 占用率和每個(gè)線(xiàn)程分配的計(jì)算量之間有一個(gè)最佳點(diǎn)，使得GPU 的性能發(fā)揮最大。通過(guò)借助CUDA profiler 進(jìn)行反復(fù)實(shí)驗(yàn)得出，當(dāng)block 的大小恒定為1616 時(shí)，加速效果最好。

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

為了驗(yàn)證該方法的可行性，本文采用實(shí)際投影數(shù)據(jù)進(jìn)行測(cè)試。掃描工件為圓柱型，掃描模式為圓軌跡。射線(xiàn)源到旋轉(zhuǎn)中心的距離為800 mm，物體中心到探測(cè)器的距離為150 mm，射線(xiàn)為ISOVOLT 450，探測(cè)器為paxscan2520，探元大小0.127 mm，電壓110 kV，電流2.5 mA，投影數(shù)據(jù)大小為570460，中心坐標(biāo)(285，228)，旋轉(zhuǎn)角度為360 度，掃描間隔為1 度。其中，實(shí)驗(yàn)平臺(tái)采用的計(jì)算機(jī)配置如下，CPU:I7-4770，GPU:GTX760。軟件環(huán)境為Visual studio 2012 +CUDA 5.5。

表1 重建過(guò)程時(shí)間對(duì)比表(單位:s)

表1 給出了具體重建過(guò)程中每個(gè)步驟的時(shí)間對(duì)比表。從中可以看出，本文方法與原GPU 方法相比，提高了約2倍，這是由于GPU 顯存具有較高的訪(fǎng)存延遲，經(jīng)CUDA 流優(yōu)化后重建速度明顯提升。與原CPU 方法相比，加速比達(dá)到了70 多倍，重建速度得到明顯改善，充分說(shuō)明了GPU 在大規(guī)模運(yùn)算中的優(yōu)勢(shì)。

圖3 原CPU 方法和本文方法重建結(jié)果對(duì)比

針對(duì)Z=100 層圖像，圖3 給出了CPU 與本文方法重建結(jié)果對(duì)比圖?？梢钥闯觯瑑煞N方法的重建結(jié)果基本一致，在精度允許范圍之內(nèi)，這是因?yàn)楸疚乃惴ǖ膶?shí)現(xiàn)方式不包含任何圖形編程內(nèi)容，只是將運(yùn)行于CPU 上的代碼改寫(xiě)為GPU上執(zhí)行的代碼。這進(jìn)一步驗(yàn)證了本文加速方法的實(shí)用性和有效性。

4 結(jié)論

針對(duì)錐束CT 三維重建耗時(shí)長(zhǎng)的問(wèn)題，本文利用GPU 中的CUDA 流和高速緩存等技術(shù)共同加速錐束CT 重建整個(gè)過(guò)程。實(shí)驗(yàn)結(jié)果表明，本文提出的方法在比原GPU 方法在重建速度上提高了近3 倍，比CPU 方法提高了約70 倍，基本滿(mǎn)足工業(yè)實(shí)時(shí)性需求。

［1］馬俊峰.三維錐束CT 圖像重建加速技術(shù)研究［D］.濟(jì)南:山東大學(xué)，2011.

［2］Zhao Y，Zou X，Yu W.Iterative Reconstruction Algorithm for Half-Cover Dual Helical Cone-Beam Computed Tomography［J］.Insight-Non-Destructive Testing and Condition Monitoring，2013，55(5):232－236.

［3］Scherl H，Keck B，Kowarschik M，et al.Fast GPU-based CT Reconstruction Using the Common Unified Device Architecture (CUDA)［C］//Nuclear Science Symposium Conference Record，2007.NSS＇07.IEEE.IEEE，2007，6:4464－4466.

［4］Zheng Z，Mueller K.Cache-aware GPU Memory Scheduling Scheme for CT Back-projection［C］.Nuclear Science Symposium Conference Record (NSS/MIC)，2010 IEEE.IEEE，2010:2248－2251.

［5］Papenhausen E，Zheng Z，Mueller K.GPU-accelerated Back-Projection Revisited:Squeezing Performance by Careful Tuning［C］//Workshop on High Performance Image Reconstruction，2011，1:1－4.

［6］張全紅，路宏年，楊民.錐束工業(yè)CT 中Feldkamp 重建算法的快速實(shí)現(xiàn)［J］.計(jì)算機(jī)工程與設(shè)計(jì)，2006(6):931－933.