基于GPU并行計(jì)算的X射線動態(tài)顯微CT的快速重構(gòu)

張 園 謝紅蘭 杜國浩 許明偉 薛艷玲 肖體喬

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

2（中國科學(xué)院大學(xué) 北京100049）

3（中國科學(xué)院上海高等研究院上海同步輻射光源 上海201204）

X射線動態(tài)顯微CT成像具有較高穿透性、高空間分辨和快時間分辨的特點(diǎn)，為開展不透明樣品三維微觀結(jié)構(gòu)演變動力學(xué)機(jī)制的探究提供了研究平臺［1-3］，在生物醫(yī)學(xué)、材料科學(xué)、考古等領(lǐng)域已有大量應(yīng)用，如科學(xué)家們利用X射線動態(tài)顯微CT原位觀察載荷下材料和結(jié)構(gòu)的動力學(xué)過程［4-6］、原位研究昆蟲馬蛉呼吸過程中的氣囊收縮特性［7］、原位分析基拉韋厄峰火山碎屑在不同壓力下氣泡產(chǎn)生的復(fù)雜變化［8］、原位觀察超聲熱浪導(dǎo)致的組織根部斷裂和枝晶的斷裂和增殖過程［9］等。近年來，隨著第三代同步輻射光源的廣泛使用，光源通量得到了極大的改善，為實(shí)現(xiàn)X射線動態(tài)顯微CT成像提供了條件［2］。上海光源X射線成像組基于長工作距離顯微透鏡系統(tǒng)和互補(bǔ)金屬氧化物半導(dǎo)體（Complementary Metal-Oxide-Semiconductor，CMOS）數(shù)字高速相機(jī)開發(fā)了快速X射線成像探測器，在彎鐵白光下，動態(tài)顯微CT成像時間分辨率可達(dá)25 Hz（每秒采集25套CT數(shù)據(jù)）［1］。目前X射線成像線站的CT重構(gòu)采用基于CPU串行重構(gòu)程序。對于一套800個投影的CT數(shù)據(jù)，每幀投影大小為1 024×1 024，按目前的重構(gòu)程序重構(gòu)需要24 min，一天的用戶實(shí)驗(yàn)產(chǎn)生的CT數(shù)據(jù)量（約500套CT）則需要重構(gòu)200 h，CT重構(gòu)速度遠(yuǎn)比原始數(shù)據(jù)的產(chǎn)生速度要慢。因此，CT重構(gòu)速度慢造成的原始數(shù)據(jù)積壓問題成為上海光源X射線成像線站遇到的瓶頸，如何實(shí)現(xiàn)快速顯微CT重構(gòu)是本文研究的重點(diǎn)。

CT重構(gòu)屬于計(jì)算與數(shù)據(jù)密集型問題，這類問題可以分解成較小的獨(dú)立單元并并行計(jì)算［10］。2010年以前，國內(nèi)外是通過依靠并行硬件平臺來加速CT重構(gòu)，其中使用廣泛的有專用集成電路（Application Specific Integrated Circuit，ASIC）、現(xiàn)場可編程門陣列（Field Programmable Gate Array，F(xiàn)PGA）［11］、數(shù)字信號處理器（Digital Signal Processor，DSP）［12］等，但此類方案硬件花費(fèi)昂貴且軟件開發(fā)難度大。2010年后，圖形處理器（Graphics Processing Unit，GPU）因其多核架構(gòu)、可大規(guī)模并行化的特點(diǎn)，在高性能計(jì)算領(lǐng)域得到廣泛應(yīng)用，隨后，NVidia的CUDA［13］和Khronos的OpenCL［14］等開發(fā)API大大簡化了基于GPU并行開發(fā)的難度。很多國內(nèi)外團(tuán)隊(duì)借助設(shè)備相對低廉的優(yōu)勢，利用開發(fā)平臺和工具API實(shí)現(xiàn)了基于GPU的CT重構(gòu)［15-17］，常用的重構(gòu)開源框架有TIGRE［18］、TomoPy［19］、ASTRA［20］、UFO［21］、PyHST［22］等。

上海光源隨著新建的快速X射線成像線站即將投入使用以及未來X射線成像探測器的升級，動態(tài)顯微CT成像將產(chǎn)生巨大的數(shù)據(jù)量，目前基于CPU串行計(jì)算的CT重構(gòu)系統(tǒng)將不能應(yīng)對海量數(shù)據(jù)的計(jì)算任務(wù)。本文借助開源框架UFO［21］，基于CUDA編程標(biāo)準(zhǔn)，將濾波反投影（Filtered Back Projection，F(xiàn)BP）的CT重構(gòu)算法的反投影運(yùn)算實(shí)現(xiàn)了并行化，并在NVDIA RTX2080的GPU上成功實(shí)現(xiàn)了CT重構(gòu)的并行計(jì)算。經(jīng)過系統(tǒng)性能測試和實(shí)驗(yàn)驗(yàn)證，相比于CPU串行計(jì)算重構(gòu)，GPU并行計(jì)算重構(gòu)獲得了200倍左右的加速比，重構(gòu)一套CT數(shù)據(jù)僅需7 s左右，大大減少了用戶處理實(shí)驗(yàn)數(shù)據(jù)的時間，且重構(gòu)結(jié)果滿足用戶實(shí)驗(yàn)要求。

1 基于GPU并行計(jì)算的CT重構(gòu)實(shí)現(xiàn)方法

CT重構(gòu)算法包括解析算法和代數(shù)迭代算法，當(dāng)數(shù)據(jù)完整時，解析算法的重構(gòu)質(zhì)量令人滿意，且重構(gòu)速度快，先濾波后反投影的FBP算法是解析算法中知名度最高和最普及的算法，因此本文采用對FBP算法并行化編程實(shí)現(xiàn)X射線動態(tài)顯微CT的快速重構(gòu)。

在三維的成像系統(tǒng)中，假設(shè)z軸代表X射線光束方向，用f(x，y，z)表示物體模型，用θ表示物體的旋轉(zhuǎn)角度，用p(x"，y")表示物體旋轉(zhuǎn)θ之后下物體的投影，投影數(shù)據(jù)由成像探測器采集。二維原始圖像f(x，z)的投影p(x")的傅里葉變換P(ω)等于f(x，z)的傅里葉變換F(ωx，ωz)沿與探測器平行的方向過原點(diǎn)的片段。如果當(dāng)物體在θ?[0，2π)內(nèi)旋轉(zhuǎn)，投影P(ω)將覆蓋整個二維傅里葉平面，由式（1）知，原始圖像f(x，z)可由濾波反投影算法計(jì)算得到［23］。

具體地說，F(xiàn)BP算法可按照三個步驟來實(shí)現(xiàn)：1）獲得投影數(shù)據(jù)p(x"，y")，并對其做傅里葉變換得到F（ωx，ωz）；2）對F（ωx，ωz）進(jìn)行濾波處理；3）對所有旋轉(zhuǎn)角度θ的p(x"，y")做反投影，即得到二維原始圖像f（x，z）。

如圖1所示，沿著y軸方向，物體可看作由很多層切片（即二維原始圖像f(x，z)）構(gòu)成，且每一層切片的重構(gòu)互不干擾，獨(dú)立執(zhí)行，這意味著CT重構(gòu)比較適用GPU的并行設(shè)計(jì)方案，可以把每個獨(dú)立的切片當(dāng)成單個重構(gòu)任務(wù)去執(zhí)行。反投影是圖像重構(gòu)中最復(fù)雜的過程，也是最耗時的過程。但是，重構(gòu)物體每一個切片之間的重構(gòu)是獨(dú)立的，可以通過并行設(shè)計(jì)有效地減少重構(gòu)時間。投影數(shù)據(jù)是按照不同投影方向（不同角度θ）得到的，反投影圖像是對不同投影方向并行求和計(jì)算得出。

圖1 基于GPU并行計(jì)算重構(gòu)的FBP算法中反投影計(jì)算過程示意圖Fig.1 Schematic diagram of back-projection calculation process in FBP algorithm based on GPU parallel computing reconstruction

本文基于GPU并行計(jì)算通過并行化FBP算法反投影過程實(shí)現(xiàn)加速。我們借助NVidia的CUDA的API接口開發(fā)運(yùn)行在GPU的反投影內(nèi)核函數(shù)［16］。CUDA是一種統(tǒng)一計(jì)算設(shè)備架構(gòu)，它可以使擅長并行計(jì)算的GPU和邏輯運(yùn)算的CPU協(xié)同工作。運(yùn)行在GPU上的程序稱為Kernel（內(nèi)核函數(shù)），一個內(nèi)核函數(shù)包含兩個層次的并行，即Grid中Block間并行和Block中Thread間并行。每個線程都有自己的BlockID和ThreadID以區(qū)別于其他線程，每個線程按照指令順序執(zhí)行Kernel函數(shù)。本文根據(jù)投影圖像大小進(jìn)行Kernel設(shè)計(jì)，線程任務(wù)劃分如圖2所示。1個Grid網(wǎng)格開啟8個Block，每個Block的尺寸是32×32，即每個Block由1 024個線程組成，每個Block負(fù)責(zé)投影圖像128行數(shù)據(jù)的重構(gòu)，每行的一個像素點(diǎn)的重構(gòu)由一個線程負(fù)責(zé)。

圖2 根據(jù)投影圖像大小對FBP算法中反投影計(jì)算過程做并行任務(wù)劃分示意圖Fig.2 Schematic diagram of divided parallel tasks in the back projection process of FBP algorithm according to the size of the projection image

在CT技術(shù)中，圖像重構(gòu)算法是決定圖像重構(gòu)質(zhì)量的關(guān)鍵因素，重構(gòu)算法的優(yōu)劣直接影響到圖像重構(gòu)結(jié)果的準(zhǔn)確性。濾波器的選取對重構(gòu)圖像質(zhì)量的影響也很明顯。反投影算法具有較低的時間復(fù)雜度，其重構(gòu)圖像的空間分辨率較高，由于反投影圖像是對不同投影方向求和計(jì)算得出，原本像素值為0的點(diǎn)的像素值變?yōu)榉?，因此反投影重構(gòu)的圖像會有偽影［24］。為了消除偽影，需要在重構(gòu)算法中加入濾波器，常用的濾波器有Hamming濾波器、Shepp-Logan濾波器等。Hamming濾波器的窗函數(shù)采用Hamming窗函數(shù)，該窗函數(shù)主瓣頻寬窄，重構(gòu)圖像空間分辨率高，旁瓣幅值小，吉布斯效應(yīng)弱［25］。對于含有噪聲的投影數(shù)據(jù)，Hamming濾波器的重構(gòu)圖像質(zhì)量比Shepp-Logan濾波器要好，因此本文中CPU串行重構(gòu)算法和GPU并行重構(gòu)算法中均選擇Hamming濾波器。

2 基于GPU并行計(jì)算的CT重構(gòu)的性能測試

本節(jié)測試基于GPU的FBP重構(gòu)算法的性能，測試設(shè)備為裝有CenOS 7系統(tǒng)的DELL T5820工作站，并配備英特爾酷睿i9-9900X處理器（3.5 GHz）和RTX 2080圖形處理器（8 GB），內(nèi)存Kingston DDR4（64 GB）。RTX 2080擁有2 944個流處理器，顯存搭配8 GB GDDR6，位寬256 bit，帶寬為448 GB·s-1。基于GPU的FBP算法的CT重構(gòu)流程主要包括以下4個步驟：從硬盤讀取投影數(shù)據(jù)；在CPU上對投影數(shù)據(jù)p(x")的一維傅里葉變換P(ω)做濾波處理；在GPU上實(shí)現(xiàn)反投影運(yùn)算；寫重構(gòu)數(shù)據(jù)到硬盤。

測試需要的實(shí)驗(yàn)數(shù)據(jù)來自10 Hz時間分辨率的CT數(shù)據(jù)采集系統(tǒng)，包括黃粉蟲體內(nèi)氣泡CT數(shù)據(jù)和黃粉蟲第一軀干CT數(shù)據(jù)。把其中一套黃粉蟲體內(nèi)氣泡CT數(shù)據(jù)分成三份進(jìn)行重構(gòu)時間的測試，分別是256張、512張和1 024張，數(shù)據(jù)像素大小為1 024×1 024?；贕PU并行計(jì)算的CT重構(gòu)算法總的重構(gòu)時間以及各個步驟的重構(gòu)時間如表1所示。

表1 在不同大小CT數(shù)據(jù)下基于GPU并行計(jì)算重構(gòu)的FBP算法的總的重構(gòu)時間以及各個步驟的重構(gòu)時間Table 1 The total reconstruction time and the reconstruction time of each step of FBP algorithm based on GPU parallel computing reconstruction with different sizes of CT data

相比于投影讀取、切片寫入和濾波運(yùn)算，反投影算法的內(nèi)存、時間、指令消耗較大，導(dǎo)致時間復(fù)雜度也很大。當(dāng)輸入規(guī)模按照同比例增加時，時間復(fù)雜度越大的算法運(yùn)行時間增長越快。因此，通過表1可以看出，隨著投影數(shù)據(jù)量的增加，反投影耗時占總時間的比重越來越大。投影數(shù)進(jìn)行四倍增長時，反投影所需時間占總的重構(gòu)時間比例由25.1%增加到49.2%。此外，總的重構(gòu)時間中，磁盤的讀寫包括投影讀取和切片寫入也占用了一半左右的時間。

像素?cái)?shù)目為1 024×1 024的投影數(shù)據(jù)分別利用基于CPU串行和基于GPU并行的FBP算法重構(gòu)，通過不同大小的CT數(shù)據(jù)來比較兩種方案的重構(gòu)時間。圖3顯示了上述兩種方案的重構(gòu)時間隨CT數(shù)據(jù)大小的變化曲線。

圖3 在不同大小CT數(shù)據(jù)下基于CPU串行計(jì)算重構(gòu)和基于GPU并行計(jì)算重構(gòu)的所需時間對比Fig.3 Comparison of reconstruction time based on CPU serial computing and GPU parallel computing with CT data of different sizes

從圖3可以看到，基于GPU并行計(jì)算的CT重構(gòu)算法縮短了CT重構(gòu)的運(yùn)行時間。當(dāng)投影數(shù)據(jù)為800幀時，基于GPU并行計(jì)算的重構(gòu)算法的運(yùn)行速度比基于CPU串行計(jì)算的重構(gòu)算法提升了近200倍。反投影過程所需時間從基于CPU串行重構(gòu)的分鐘級減少到了基于GPU并行重構(gòu)的秒級。

對于基于CPU串行程序，數(shù)據(jù)加載路徑是從硬盤到本地內(nèi)存，對于基于GPU并行程序，數(shù)據(jù)加載路徑額外包括從本地內(nèi)存到GPU顯存。隨著數(shù)據(jù)量的不斷增加，CPU讀取硬盤數(shù)據(jù)及本地內(nèi)存與GPU顯存數(shù)據(jù)傳輸時間不斷增加。對于不同大小的投影數(shù)據(jù)，基于GPU并行和基于CPU串行的重構(gòu)時間加速比隨CT數(shù)據(jù)大小的變化曲線如圖4所示。

圖4 相比于基于CPU串行計(jì)算重構(gòu)，基于GPU并行計(jì)算重構(gòu)獲得的加速比隨CT數(shù)據(jù)大小的變化曲線Fig.4 Curve of speedup of GPU based parallel reconstruction compared with CPU based serial reconstruction vs.CT data size

如圖4所示，投影數(shù)據(jù)量當(dāng)增加到一定規(guī)模時，圖像重構(gòu)的加速比將會隨著數(shù)據(jù)量的增加而減少。其主要原因是基于GPU并行計(jì)算的CT重構(gòu)算法在反投影計(jì)算部分起到了加速效果，而在投影讀取和切片寫入部分并沒有得到加速，由表1可知，隨著投影數(shù)據(jù)的增加，投影讀取和切片寫入部分的時間占總重構(gòu)時間的50%左右，其中，硬盤和本地內(nèi)存之間的投影讀取和切片寫入時間占主要部分，這是CT重構(gòu)總體加速比接近飽和狀態(tài)的主要原因。提高硬件能夠進(jìn)一步加快CT重構(gòu)算法中反投影計(jì)算部分的速度，但不能改變投影讀取和切片寫入部分的速度。在重構(gòu)過程中，提高數(shù)據(jù)傳輸速度成為解決問題的關(guān)鍵。

3 基于GPU并行計(jì)算的CT重構(gòu)的實(shí)驗(yàn)驗(yàn)證

X射線動態(tài)顯微CT實(shí)驗(yàn)的樣品為一只活的黃粉蟲，將其置于直徑2 mm的聚酰亞胺管中。實(shí)驗(yàn)在上海同步輻射光源的BL09B測試線站進(jìn)行，采用通量較高的彎鐵光源白光，輸出X射線的光子能量在6 keV以上。在樣品前2 m處放置1 mm厚鋁箔片以吸收低能X射線光子，降低樣品受到的X射線輻射。實(shí)驗(yàn)樣品以及實(shí)驗(yàn)裝置如圖5所示，實(shí)驗(yàn)樣品臺由三部分組成：從上到下依次是微型位移臺（SmarAct SLS-5252，分辨率小于1 nm）、氣浮轉(zhuǎn)臺（LAB RT150S，最高轉(zhuǎn)速800 rpm）和六軸位移臺（PI H-825，最小位移達(dá)0.25μm）。微型位移臺用于調(diào)整樣品的旋轉(zhuǎn)中心，氣浮轉(zhuǎn)臺用于樣品的高精度穩(wěn)定旋轉(zhuǎn)實(shí)現(xiàn)CT掃描，六軸六足位移臺用于調(diào)節(jié)視場中樣品的位置。上海光源X射線成像組開發(fā)了基于長工作距離顯微鏡頭系統(tǒng)和高速CMOS相機(jī)的快速X射線成像探測器，通過優(yōu)化長距離顯微鏡頭系統(tǒng)，增大其數(shù)值孔徑（8倍鏡頭的數(shù)值孔徑為0.5，4倍鏡頭的數(shù)值孔徑為0.22），提高了耦合效率，為了匹配鏡頭數(shù)值孔徑，閃爍體選擇厚度為50μm，可實(shí)現(xiàn)最佳成像分辨率和最佳成像襯度。為了達(dá)到最優(yōu)的實(shí)驗(yàn)效果，仔細(xì)選擇了實(shí)驗(yàn)參數(shù)，光子能量設(shè)置為18 keV，樣品到探測器的距離設(shè)置為46 cm，氣浮轉(zhuǎn)臺每分鐘旋轉(zhuǎn)300圈，相機(jī)每秒拍攝8 000幀圖片，每套CT數(shù)據(jù)為800幀，每套CT的采集時間為0.1 s。實(shí)驗(yàn)數(shù)據(jù)來自10 Hz時間分辨率的CT數(shù)據(jù)采集系統(tǒng)，包括黃粉蟲體內(nèi)氣泡CT數(shù)據(jù)和黃粉蟲第一軀干CT數(shù)據(jù)，黃粉蟲體內(nèi)氣泡CT數(shù)據(jù)用4倍鏡頭的探測器采集，有效像素5μm，探測器視場為5 mm，實(shí)驗(yàn)拍攝時間為5 s。黃粉蟲第一軀干CT數(shù)據(jù)用8倍鏡頭的探測器采集，有效像素2.5μm，探測器視場為2.5 mm，實(shí)驗(yàn)拍攝時間也為5 s。

圖5 X射線動態(tài)顯微CT實(shí)驗(yàn)樣品(a)以及實(shí)驗(yàn)裝置圖(b)Fig.5 Photographs of X-ray dynamic microscopic CT experimental sample(a)and experimental setup(b)

3.1 基于CPU串行計(jì)算和基于GPU并行計(jì)算的CT重構(gòu)圖像質(zhì)量對比

基于CPU串行計(jì)算和基于GPU并行計(jì)算分別對黃粉蟲體內(nèi)氣泡CT數(shù)據(jù)進(jìn)行重構(gòu)，選擇同一層CT重構(gòu)切片做線型分析，結(jié)果如圖6所示。由于CPU串行重構(gòu)算法和GPU并行重構(gòu)算法均采取了Hamming濾波器，重構(gòu)切片的偽影均得到有效的消除。重構(gòu)切片中的圈表示線型分析的路徑，GPU并行計(jì)算重構(gòu)切片在明暗交界處像素點(diǎn)灰度值的變化明顯，其襯度值大于CPU串行計(jì)算重構(gòu)結(jié)果。

圖6 基于CPU串行計(jì)算重構(gòu)(a)和基于GPU并行計(jì)算重構(gòu)(b)得到的黃粉蟲體內(nèi)氣泡的CT重構(gòu)切片及其線型分析比較Fig.6 Comparisons of CT reconstruction slices and line profile analyses of bubbles in the body of the tenebrio molitor sample obtained from CPU serial reconstruction(a)and GPU parallel reconstruction(b)

基于CPU串行計(jì)算和基于GPU并行計(jì)算分別對黃粉蟲第一軀干CT數(shù)據(jù)進(jìn)行重構(gòu)。選擇同一層CT重構(gòu)切片計(jì)算其灰度值分布，結(jié)果如圖7所示。

通過CT重構(gòu)切片的灰度值可以計(jì)算其襯度，襯度是指圖像上不同區(qū)域間存在的明暗程度差異，數(shù)值越大，越易于分辨圖像細(xì)節(jié)。其中襯度值計(jì)算公式如下：

式中：C表示襯度值；G表示灰度值。

結(jié)合圖7和式（2）計(jì)算可知，基于CPU串行計(jì)算重構(gòu)結(jié)果的襯度值為25.4%，基于GPU并行計(jì)算重構(gòu)結(jié)果的襯度值為33.3%。根據(jù)襯度值定義，圖像襯度值越大，圖像明暗變化越大，圖像細(xì)節(jié)信息越明顯，從圖8可以看出，基于GPU并行計(jì)算的三維CT重構(gòu)結(jié)果除了可以看到黃粉蟲體內(nèi)直腸外，黃粉蟲的直腸周圍細(xì)節(jié)信息也能清晰的觀察到，經(jīng)測量，直腸直徑為47.5μm。

圖7 基于CPU串行計(jì)算重構(gòu)(a)和基于GPU并行計(jì)算重構(gòu)(b)得到的黃粉蟲第一軀干CT重構(gòu)切片及灰度值分布圖的比較Fig.7 Comparisons of the CT reconstruction slice and the gray-scale value distribution map of the first trunk of the tenebrio molitor sample based on CPU serial computing reconstruction(a)and GPU parallel computing reconstruction(b)

圖8 基于CPU串行計(jì)算重構(gòu)(a)和基于GPU并行計(jì)算重構(gòu)(b)獲得的黃粉蟲第一軀干的三維CT重構(gòu)結(jié)果比較Fig.8 Comparisons of the 3D CT reconstruction results of the first trunk of the tenebrio molitor sample based on CPU serial computing reconstruction(a)and GPU parallel computing reconstruction(b)

基于CPU串行重構(gòu)算法和基于GPU并行重構(gòu)算法的像素點(diǎn)灰度值計(jì)算方式不同，造成了重構(gòu)切片襯度值的不同，最終使得樣品三維CT重構(gòu)結(jié)果的內(nèi)部細(xì)節(jié)信息在觀察上存在差異。圖像反相算法是指對圖像中的像素點(diǎn)做灰度值計(jì)算，其計(jì)算方法如式（3）所示，重構(gòu)切片位深為8位，i0為初始像素值，i1為圖像反相后的像素值。

對CPU串行計(jì)算重構(gòu)切片做圖像反相計(jì)算，其襯度值由25.4%提高到31.5%，結(jié)果如圖9所示。因此在基于CPU串行計(jì)算重構(gòu)中，可以進(jìn)一步采用圖像反相算法來提高重構(gòu)圖像的襯度，增強(qiáng)圖像中吸收X射線較弱的信號，這對用戶觀察和分析樣品內(nèi)部結(jié)構(gòu)中的弱信號有著一定實(shí)際意義。

圖9 對CPU串行計(jì)算重構(gòu)切片做反相處理得到的黃粉蟲第一軀干CT重構(gòu)切片灰度值分布圖(a)以及三維CT重構(gòu)結(jié)果(b)Fig.9 The gray-scale value distribution of CT reconstruction slice(a)and 3D CT reconstruction result(b)of the first trunk of the tenebrio molitor sample obtained by inverse processing of serial calculation reconstruction slice of CPU

3.2 基于GPU并行計(jì)算在不同濾波器下CT重構(gòu)質(zhì)量的對比

為了比較不同濾波器的性能差異，基于GPU并行計(jì)算重構(gòu)分別選擇不同的濾波器處理投影數(shù)據(jù)。圖10（a、c）是采用Shepp-Logan濾波器處理得到的黃粉蟲第一軀干的三維CT重構(gòu)結(jié)果的正視圖，圖10（b、d）采用Hamming濾波器處理得到的黃粉蟲第一軀干的三維CT重構(gòu)結(jié)果的俯視圖?？梢钥闯觯瑘D10（b、d）偽影更少，圖像重構(gòu)質(zhì)量較好，采用Hamming濾波器的CT重構(gòu)算法實(shí)驗(yàn)效果更好。

3.3 基于GPU并行計(jì)算的動態(tài)顯微CT的快速重構(gòu)

基于GPU并行計(jì)算的CT重構(gòu)算法大大縮短了CT重構(gòu)的運(yùn)行時間，選擇Hamming濾波器可以減少反投影計(jì)算產(chǎn)生的偽影。最后通過設(shè)計(jì)動態(tài)顯微CT實(shí)驗(yàn)對基于GPU并行計(jì)算和Hamming濾波器的CT重構(gòu)算法進(jìn)行了驗(yàn)證。50套黃粉蟲體內(nèi)氣泡CT數(shù)據(jù)完成重構(gòu)總共耗時5.93 min，而基于CPU串行重構(gòu)耗時20.4 h。選取不同時間點(diǎn)的氣泡CT重構(gòu)數(shù)據(jù)，可以在時間維度上觀察到氣泡的動態(tài)變化過程。圖11是基于GPU并行計(jì)算實(shí)現(xiàn)黃粉蟲體內(nèi)氣泡在不同時刻的運(yùn)動過程的三維CT重構(gòu)結(jié)果，CT重構(gòu)圖像沒有偽影產(chǎn)生，可以清晰觀察到黃粉蟲體內(nèi)氣泡在自上而下運(yùn)動，通過計(jì)算，氣泡半徑為0.624 mm，氣泡平均運(yùn)動速度為1.08 mm·s-1。

圖10 基于GPU并行計(jì)算的CT重構(gòu)中，分別采用Shepp-Logan濾波器和Hamming濾波器獲得的黃粉蟲第一軀干的三維CT重構(gòu)結(jié)果比較 (a、c)為采用Shepp-Logan濾波器的三維重構(gòu)結(jié)果的正視圖和俯視圖，(b、d)為采用的Hamming濾波器的三維重構(gòu)結(jié)果的正視圖和俯視圖Fig.10 Comparisons of the 3D CT reconstruction results of the first trunk of the tenebrio molitor sample using Shepp-Logan filter and Hamming filter in GPU parallel computing reconstruction (a,c)The front view and the top view of the results using Shepp-Logan filter,(b,d)The front view and the top view of the results using Hamming filter

圖11 基于GPU并行計(jì)算重構(gòu)實(shí)現(xiàn)黃粉蟲體內(nèi)氣泡在不同時刻的運(yùn)動過程的三維CT重構(gòu)結(jié)果Fig.11 The 3D CT reconstruction results of the movement process of bubbles in the body of the tenebrio molitor sample at different times based on GPU parallel computing reconstruction

4 結(jié)語

本文基于GPU通過并行化反投影實(shí)現(xiàn)了FBP重構(gòu)算法，借助GPU強(qiáng)大的計(jì)算能力可以實(shí)現(xiàn)快速顯微CT重構(gòu)，相比于CPU串行重構(gòu)，GPU并行計(jì)算獲得了200倍左右的加速比。選擇同一層重構(gòu)切片做X射線線型分析并進(jìn)行襯度計(jì)算，基于GPU并行計(jì)算重構(gòu)結(jié)果的襯度值優(yōu)于基于CPU串行計(jì)算重構(gòu)結(jié)果。在基于CPU串行計(jì)算重構(gòu)中，進(jìn)一步采用圖像反相算法可以提高重構(gòu)圖像的襯度。相比Shepp-Logan濾波器，Hamming濾波器能減少重構(gòu)圖像的偽影，改善重構(gòu)圖像的質(zhì)量。基于GPU并行計(jì)算并采用Hamming濾波器，在6 min內(nèi)實(shí)現(xiàn)了50套X射線動態(tài)顯微CT實(shí)驗(yàn)數(shù)據(jù)的快速重構(gòu)，并能夠清晰地看到黃粉蟲體內(nèi)氣泡在不同時刻的運(yùn)動過程。實(shí)驗(yàn)結(jié)果表明，基于GPU并行計(jì)算并采用Hamming濾波器的CT重構(gòu)算法能夠?qū)崿F(xiàn)海量CT實(shí)驗(yàn)數(shù)據(jù)的快速重構(gòu)，將很大程度上解決上海光源X射線成像線站的大量用戶實(shí)驗(yàn)數(shù)據(jù)積壓的瓶頸問題，顯著提升實(shí)驗(yàn)效率。