沈 樂(lè) 邢宇翔,2

1 (清華大學(xué)工程物理系 北京 100084)

2 (清華大學(xué)粒子技術(shù)與輻射成像教育部重點(diǎn)實(shí)驗(yàn)室 北京 100084)

錐束螺旋X射線(xiàn)CT近年來(lái)受到廣泛關(guān)注，其重建算法是當(dāng)前CT領(lǐng)域的研究熱點(diǎn)。目前，解析重建方法有近似重建方法和精確重建方法，后者是大家追求的目標(biāo)[1,2]。2004年，ZOU 等[3]提出了微分反投影濾波重建方法(DBPF)，能在最少數(shù)據(jù)條件下進(jìn)行精確重建，很多相關(guān)工作(包括軌道的推廣和局部重建等)隨之展開(kāi)。但其具體實(shí)現(xiàn)面臨計(jì)算復(fù)雜度和離散化等問(wèn)題，至今，商用系統(tǒng)仍未使用此類(lèi)方法。因此，尋找一種好的算法實(shí)現(xiàn)策略以及如何對(duì)算法進(jìn)行加速使其應(yīng)用于實(shí)際非常關(guān)鍵。

利用顯卡硬件的并行能力對(duì)CT重建算法加速是目前加速效率最高的一種方式。由于重建算法具有數(shù)據(jù)相互獨(dú)立性及可高度并行處理特點(diǎn)，在GPGPU上設(shè)計(jì)并行重建算法，能獲得很好的性能提升。對(duì)于平行束、扇形束的FBP算法，以及錐束FDK算法，已獲得良好結(jié)果[4–6]。對(duì)于PI線(xiàn)BPF重建算法，文獻(xiàn)[7,8]用 Cg語(yǔ)言處理重建的反投影過(guò)程，重建速度提升了100倍。

我們對(duì)DBPF方法的實(shí)施技巧和硬件加速進(jìn)行了研究，給出了一種新的實(shí)現(xiàn)方法和加速策略，以在最少不均勻采樣基礎(chǔ)上提供可控分辨率重建，并最大程度地利用 GPU的并行加速性能全面提高容積重建速度。

1 DBPF算法和PI線(xiàn)離散方法

螺旋錐束CT的掃描坐標(biāo)系和參數(shù)的定義如圖1所示，射線(xiàn)源和探測(cè)器按螺旋軌道u運(yùn)v動(dòng)，射線(xiàn)源的位置在直角坐標(biāo)系{x,y,z}中表示為Rsinλ,hλ/2π)，其中，R為射線(xiàn)源旋轉(zhuǎn)的半徑，λ為旋轉(zhuǎn)角度，h為螺距，即射線(xiàn)源旋轉(zhuǎn)一周相對(duì)于置物臺(tái)的軸向移動(dòng)距離。定義探測(cè)器旋轉(zhuǎn)坐標(biāo)系{u,v,w}，其基向量為：

圖1 螺旋錐束CT坐標(biāo)系定義Fig.1 Definition of coordinates for the helical conebeam CT.

2004年芝加哥大學(xué)潘曉川教授[3]提出的基于PI線(xiàn)的反投影濾波算法的步驟如下：

(1) 對(duì)投影數(shù)據(jù)求導(dǎo)：

(2) 沿PI線(xiàn)進(jìn)行加權(quán)反投影：由求導(dǎo)后的數(shù)據(jù)G{λ,u,v}在 PI線(xiàn)所在角度范圍λ1–λ2進(jìn)行積分：

(3) 對(duì)加權(quán)反投影結(jié)果進(jìn)行Hilbert變換，得到PI線(xiàn)上的重建結(jié)果：

這里，附加項(xiàng)xA,xB為PI線(xiàn)與重建區(qū)域的交點(diǎn)，P{u,v,λ}為PI線(xiàn)所在位置的射線(xiàn)投影，我們?nèi){u,v,λ} = [P{u,v,λ1}+P{u,v,λ2}]/2。

(4) 將PI線(xiàn)重建結(jié)果按三線(xiàn)性插值方式重新采樣成為均勻網(wǎng)格的三維重建圖像。在BPF重建中，常規(guī)的 PI線(xiàn)選取為從一個(gè)光源角度采樣點(diǎn)到其它多個(gè)光源角度采樣點(diǎn)形成的一個(gè)扇形曲面形式，這不利于建立重建圖像分辨率與PI線(xiàn)密度和PI線(xiàn)上采樣密度的關(guān)系。如圖2，使一組PI線(xiàn)(圖2中實(shí)線(xiàn))在俯視圖中(圖2c)互相平行且距離為DG，采樣點(diǎn)在XY平面上投影的距離為DP。我們把這樣一組PI線(xiàn)構(gòu)成的曲面稱(chēng)為一個(gè)PI面，圖2的虛線(xiàn)所示構(gòu)成另一個(gè)PI面。在選擇軸方向上相鄰的PI面時(shí)，令它們的PI線(xiàn)在XY平面上的投影成一角度Dθ，由重建需要的層厚確定。取兩個(gè)PI面的中心點(diǎn)(0,0,z)作為基準(zhǔn)點(diǎn)(也可選其它點(diǎn)作為基準(zhǔn)點(diǎn)，如圖2b所示，效果相同)，則Dθ由Z軸采樣間隔DZ和螺距的線(xiàn)性關(guān)系確定。如此得到的PI線(xiàn)采樣點(diǎn)在重建區(qū)域內(nèi)具有全局一致性，可方便地控制最后的重建分辨率。如DG=DP=DZ，可得到各向同性的具有均勻分辨率效果的 PI線(xiàn)重建采樣點(diǎn)。用式(4)得到的重建結(jié)果是沿PI線(xiàn)的數(shù)據(jù)，須按照規(guī)則網(wǎng)格重新采樣為三維體數(shù)據(jù)。由于我們的PI線(xiàn)選取方式得到全局基本均勻的重建點(diǎn)，可方便地使用三線(xiàn)性插值方式重采樣得到重建圖像。

731 Comparison of different separation protocols for exosomes derived from urine of patients with aortic dissection

圖2 PI 線(xiàn)選取和采樣點(diǎn) (a)全景圖，(b)側(cè)視圖，(c)俯視圖Fig.2 PI-lines and sampling point selection. (a) Perspective view, (b) Top view, (c) Side view.

式中，(xΠ,yΠ,zΠ)表示 PI線(xiàn)上像素的坐標(biāo)，ax,ay,az為該像素與均勻網(wǎng)格像素之間距離的三個(gè)分量，D為均勻網(wǎng)格的層內(nèi)離散化步距，Ds為均勻網(wǎng)格的層厚，Neighbor(i,j,k)表示(i,j,k)鄰域(圖 3)。

圖3 三線(xiàn)性插值重采樣的示意圖Fig.3 Illustration of tri-linear interpolation resampling.

2 CUDA加速優(yōu)化策略

隨著半導(dǎo)體技術(shù)及芯片制造工藝的不斷突破，圖形處理器的流式核心數(shù)量、浮點(diǎn)運(yùn)算能力不斷提升。近年來(lái)快速發(fā)展的通用圖形處理器(GPGPU)除圖形圖像的顯示處理外，還可完成很多原本由CPU完成的計(jì)算任務(wù)。目前的高端圖形處理器(如美國(guó)NVIDIA公司的 GF100核心和美國(guó) AMD公司的RV870核心)的浮點(diǎn)運(yùn)算能力已超過(guò)2 TFLOPS(1012浮點(diǎn)運(yùn)算/s)，已大大超越多核CPU(如美國(guó)Intel公司芯片)的計(jì)算能力(表 1)。圖形處理器本身系并行處理器，對(duì)于可并行算法，通用圖形處理器的執(zhí)行速度遠(yuǎn)高于CPU。

CUDA(Compute Unified Device Architecture)[9]是美國(guó)NVIDIA公司提出的并行架構(gòu)，可利用圖形處理器的并行處理器能力大幅度提升計(jì)算性能。CUDA包含一系列硬件架構(gòu)、一套指令集和編程模型，針對(duì)C/C++/Fortran語(yǔ)言的應(yīng)用程序編程接口和編譯器。開(kāi)發(fā)人員利用相應(yīng)軟件工具，可在CUDA架構(gòu)的GPU上開(kāi)發(fā)并行程序應(yīng)用及加速算法。

我們利用CUDA架構(gòu)，在GPU上完成：(1) 對(duì)投影數(shù)據(jù)求導(dǎo)；(2) 沿 PI線(xiàn)加權(quán)反投影；(3) 有限Hilbert變換；(4) 將PI線(xiàn)重建結(jié)果按圖像網(wǎng)格重采樣。具體步驟如下：

(1) 求導(dǎo)：將投影數(shù)據(jù)傳入顯卡的紋理顯存中，將探測(cè)器面陣所得的投影數(shù)據(jù)映射到一個(gè)二維線(xiàn)程組中，每一個(gè)線(xiàn)程對(duì)應(yīng)計(jì)算某一點(diǎn)的全導(dǎo)數(shù)。

(2) 反投影：在一個(gè)PI面內(nèi)的一組PI線(xiàn)，不僅起點(diǎn)終點(diǎn)角度不同、長(zhǎng)短不同，PI線(xiàn)上的采樣點(diǎn)數(shù)目也不同，因此無(wú)法將這些PI線(xiàn)上的采樣點(diǎn)直接映射為規(guī)則的二維線(xiàn)程組。若采用最長(zhǎng)PI線(xiàn)的采樣點(diǎn)數(shù)作為一個(gè)映射維度，則其他PI線(xiàn)映射的線(xiàn)程又不能充分利用，對(duì)于顯卡并行模型是不經(jīng)濟(jì)的。CUDA并行編程中應(yīng)盡可能使顯卡上的線(xiàn)程達(dá)到負(fù)載平衡。為此，我們沿Z方向選取長(zhǎng)度相同(圖4)，即起始和終止角度差Dλ=λ2?λ1相等的一組 PI線(xiàn)并發(fā)執(zhí)行反投影操作。不同的PI線(xiàn)僅λ1、λ2的參數(shù)不同，長(zhǎng)度及采樣點(diǎn)數(shù)量都相同。實(shí)驗(yàn)證明此法能有效匹配負(fù)載平衡，提高流處理器的利用率，大大加速反投影運(yùn)算。

(3) 有限Hilbert變換：將一組PI線(xiàn)反投影結(jié)果映射到二維線(xiàn)程組中，每一線(xiàn)程計(jì)算一個(gè)采樣點(diǎn)的卷積結(jié)果。計(jì)算前將反投影數(shù)據(jù)和卷積核綁定到一維紋理中，以利用紋理緩存命中提高數(shù)據(jù)讀取速度。

(4) 重采樣：將一組PI線(xiàn)的上述結(jié)果映射到二維線(xiàn)程組中，按圖3方式將結(jié)果分配到采樣點(diǎn)周?chē)?個(gè)鄰域像素上，并累加8個(gè)像素所占權(quán)重值。當(dāng)重建全部完成后，使用累加的權(quán)重值對(duì)結(jié)果歸一化處理，得到最終圖像。

圖4 反投影的并行方式示意圖Fig.4 Illustration of paralleled back-projection.

為提高算法效率，我們提出了優(yōu)化手段，包括打包存儲(chǔ)重建所需常量參數(shù)，應(yīng)用支持線(xiàn)性插值方式的三維紋理存儲(chǔ)器保存投影數(shù)據(jù)等，取得了較好效果。

首先，在加權(quán)反投影及Hilbert變換中可能用到的重建參數(shù)超過(guò)20個(gè)。若全部作為函數(shù)調(diào)用參數(shù)來(lái)傳遞，則kernel的寄存器占用數(shù)會(huì)增加27%，流處理器占用率明顯減少。對(duì)于計(jì)算能力<1.3的硬件，可能會(huì)因硬件資源不足導(dǎo)致kernel無(wú)法啟動(dòng)。為減少寄存器占用，我們將重建所需的參數(shù)打包存于數(shù)組中，預(yù)先輸入顯卡的常數(shù)存儲(chǔ)器或共享存儲(chǔ)器中以供調(diào)用，有效解決了每個(gè)kernel硬件資源占用過(guò)多的問(wèn)題。

其次，CUDA的紋理存儲(chǔ)器是一種只讀存儲(chǔ)器，具有緩存機(jī)制，同時(shí)支持硬件線(xiàn)性插值。直接訪問(wèn)三維浮點(diǎn)坐標(biāo)的數(shù)據(jù)即可得按三線(xiàn)性插值的結(jié)果，無(wú)需編程實(shí)現(xiàn)。緩存機(jī)制還可提高內(nèi)存的實(shí)際訪問(wèn)帶寬，降低延遲。

最后，在并行重采樣時(shí)，若PI線(xiàn)上的采樣點(diǎn)間隔較小，相鄰的采樣點(diǎn)可能具有公共的鄰域像素。此時(shí)若并行進(jìn)行重采樣，會(huì)引起不同線(xiàn)程對(duì)同一內(nèi)存地址進(jìn)行寫(xiě)操作。為避免線(xiàn)程間的訪問(wèn)沖突，在反投影濾波結(jié)束后，將PI線(xiàn)重建數(shù)據(jù)傳回內(nèi)存，由CPU進(jìn)行串行重采樣。這一步可在子線(xiàn)程上進(jìn)行，而CPU主線(xiàn)程則控制GPU進(jìn)行下一輪反投影濾波運(yùn)算，如此可充分利用GPU計(jì)算時(shí)CPU的空閑時(shí)間。只要在數(shù)據(jù)傳輸前設(shè)置同步，就可保證訪問(wèn)是安全的。這種CPU+GPU同時(shí)計(jì)算的方法，在重建數(shù)據(jù)較小、PI線(xiàn)采樣步距較大時(shí)可有效隱藏重采樣時(shí)間。

對(duì)于大數(shù)據(jù)量處理，同步導(dǎo)致的空閑時(shí)間明顯增加。為此，我們?cè)O(shè)計(jì)了第二種GPU上執(zhí)行的方法：分批并行重采樣，在加速的同時(shí)避免寫(xiě)操作的沖突。在二維線(xiàn)程組映射時(shí)，把 PI線(xiàn)上的采樣點(diǎn)分批處理。如圖5所示，圓點(diǎn)標(biāo)記的PI線(xiàn)采樣點(diǎn)映射為一個(gè)批次的線(xiàn)程組，完成后把方形標(biāo)記PI線(xiàn)采樣點(diǎn)映射為另一批次的線(xiàn)程組進(jìn)行計(jì)算，如此繼續(xù)直至處理完所有的PI線(xiàn)上采樣點(diǎn)。選取每一個(gè)批次內(nèi)處理的PI線(xiàn)采樣點(diǎn)的合適間距，可保證采樣點(diǎn)不會(huì)有公共的被重建鄰域像素，避免了線(xiàn)程訪問(wèn)沖突。

圖5 并行重采樣方式示意圖Fig.5 Illustration of paralleled resampling.

3 結(jié)果與分析

仿真實(shí)驗(yàn)中，我們采用三維Shepp Logan 模型生成螺旋錐束投影數(shù)據(jù)，光源到旋轉(zhuǎn)軸和探測(cè)器的距離分別為650 mm和1000 mm，每個(gè)2π采樣1080個(gè)角度，探測(cè)器為 256×80，每個(gè)單元尺寸為 1.8 mm×1.8 mm，掃描螺距46.7 mm。算法使用Visual Studio 2005，CUDA Toolkit 3.0編程設(shè)計(jì)，程序在一臺(tái)AMD Athlon II X4 620 (2.60 GHz)計(jì)算機(jī)上進(jìn)行測(cè)試，具有DDR2 800 MHz 2 GB RAM和NVIDIA Tesla C1060 GPU。我們分別使用了純CPU執(zhí)行和GPU加速測(cè)試重建256×256×64的結(jié)果，D = 1 mm,采樣間隔DG=DP=DZ=0.8D。有效重建區(qū)域覆蓋52層(3.4×106像素)，在CPU上重建時(shí)間787 s，使用GPU加速后的重建時(shí)間為2.469 s，重建加速比達(dá)到318倍。圖6為重建圖像比較，顯示CPU與GPU重建結(jié)果非常吻合。

圖6 第28層重建結(jié)果比較(a) CPU重建結(jié)果, (b) CUDA加速后的結(jié)果, (c) 圖像(a)(b)第92行的水平剖面線(xiàn)Fig.6 Reconstructions of Slice 28.(a) CPU result, (b) CUDA accelerated result, (c) Horizontal profiles at Row 92.

在以上的重建中，重建速度及質(zhì)量與 Hilbert卷積核選取、附加項(xiàng)計(jì)算、采樣間隔等參數(shù)密切相關(guān)。為避免混疊，卷積核長(zhǎng)度通常為2Np?1，Np表示 PI線(xiàn)上的采樣點(diǎn)數(shù)目?？s短卷積核長(zhǎng)度可加快Hilbert變換的運(yùn)算速度，但會(huì)引入偽影。固定最后的重建分辨率Δ條件下，重建過(guò)程中PI線(xiàn)間的距離DG，PI線(xiàn)上的采樣間距ΔP，Z軸方向采樣間隔DZ三個(gè)參數(shù)的取值會(huì)對(duì)重建結(jié)果及速度產(chǎn)生影響。在需要重建圖像各向同性分辨率情況下，通常取(但不嚴(yán)格要求)DG=DP=DZ。顯然，采樣間隔越小，重建圖像質(zhì)量越好，重建時(shí)間越長(zhǎng)。根據(jù)經(jīng)驗(yàn)，一般可設(shè)采樣間隔(DG, DP, DZ)～0.8D，即使 PI線(xiàn)的采樣略細(xì)于最后的重建分辨率要求即可得到滿(mǎn)意的重建圖像和高的重建速度。圖7顯示不同采樣參數(shù)組的重建效果。圖8是不同參數(shù)情況下的重建時(shí)間變化情況，這里我們?cè)O(shè)置ΔG= ΔP。

圖7 第20層不同采樣參數(shù)的重建圖像Fig.7 Reconstructions of Slice 20 with different parameters.

圖8 不同PI線(xiàn)選取和采樣網(wǎng)格下的重建時(shí)間比較(ΔG=ΔP的情況)Fig.8 Reconstruction time comparison using different sampling grids for PI line setting (ΔG = ΔP).

4 結(jié)論

本文從DBPF重建算法的原理和目前通用圖形處理器與并行重建算法的技術(shù)入手，介紹了一種新的實(shí)現(xiàn)方法，以及使用CUDA平臺(tái)下算法進(jìn)行高效加速的策略和優(yōu)化手段。充分利用了GPU的并行處理與CPU+GPU混合計(jì)算的計(jì)算能力，達(dá)到了很好的加速效果。

該算法仍有進(jìn)一步改進(jìn)的可能，重建過(guò)程中的計(jì)算精度由于使用內(nèi)置快速數(shù)學(xué)函數(shù)會(huì)帶來(lái)一定誤差。2010年4月美國(guó)NVIDIA公司發(fā)布的最新圖形處理核心 GF-100具有更高的單精度和雙精度浮點(diǎn)運(yùn)算能力，并且支持IEEE 754-2008單精度浮點(diǎn)數(shù)標(biāo)準(zhǔn)。在最新架構(gòu)下，可嘗試使用高精度的標(biāo)準(zhǔn)數(shù)

學(xué)函數(shù)，以獲取更優(yōu)的圖像質(zhì)量。此外，反投影過(guò)程是所有步驟中耗時(shí)最長(zhǎng)的操作，有可能通過(guò)繼續(xù)提高反投影過(guò)程中的紋理緩存命中率，減少分支判斷而提高速度。我們將在后續(xù)工作中繼續(xù)研究這些技術(shù)。

1 Tuy H, SIAM J. Appl Math, 1983, 43(3): 546–552

2 Katsevich A. Phys Med Biol, 2002, 47: 2583–2597

3 ZOU Yu, PAN Xiaochuan. Phys Med Biol, 2004, 49:941–959

4 Wang Y J, Hu H F, Xing Y X. Strategy for GPU acceleration of massive data cone beam CT reconstruction.The 10thinternational meeting on fully 3D image reconstruction in radiology and nuclear medicine, Beijing,2009. 57–60

5 Bi W Y, Xing Y X, Chen Z Q,et al. 2008 IEEE Nucl Sci Symp and Med Imag Conf (2008 NSS/MIC). 2009, 1–9:2605–2608

6 Xu F, Mueller K. Phys Med Biol, 2007, 52: 3405–3419

7 Zheng H, Yu Y, Kang Y,et al. Proc SPIE, 2010, 7622:76222G-1–76222G-12

8 Zheng H, Kang Y, Liu J,et al. Implementation of Helical Cone-Beam Back-Projection Filtered Reconstruction Algorithm on GPU. The 10thinternational meeting on fully 3D image reconstruction in radiology and nuclear medicine, Beijing, 2009. 45–48

9 NVIDIA Corporation. NVIDIA CUDA Best Practices Guide [EB/OL]. http://developer. download. nvidia.com/compute/cuda/3_0/toolkit/docs/NVIDIA_CUDA_BestPra cticesGuide.pdf, 2010