GPU加速技術(shù)在醫(yī)學(xué)三維可視化中的應(yīng)用

王 玉，王 宏

(1.東北大學(xué) 中荷生物醫(yī)學(xué)與信息工程學(xué)院，沈陽 110819;2.東北大學(xué) 機械工程與自動化學(xué)院，沈陽110004)

近些年來，隨著CT等影像設(shè)備的不斷發(fā)展和升級換代，使得患者影像數(shù)據(jù)量越來越大，傳統(tǒng)的二維顯示方式已無法滿足實際臨床需要。醫(yī)學(xué)三維可視化技術(shù)近些年來成為被關(guān)注的領(lǐng)域。

臨床應(yīng)用中除了對三維可視化重建圖像質(zhì)量提出較高的要求外，對交互繪制速度也提出了較高的要求?；诠饩€投射方法可以生成高質(zhì)量的三維圖像，能夠滿足臨床診斷的需要。但它的一個突出缺點就是繪制速度慢，如果沒有一種有效的加速方法對其繪制速度進(jìn)行改善，其很難適應(yīng)實際應(yīng)用和醫(yī)學(xué)圖像技術(shù)發(fā)展的要求。

光線投射方法提出后，相應(yīng)的加速技術(shù)的研究就從來沒有停止過。目前，針對光線投射的軟件加速技術(shù)主要有以下幾種［1－2］:包圍盒技術(shù)(bounding box)、空間剖分(space subdivision)技術(shù)和光線相關(guān)性(ray coherence)技術(shù)。

包圍盒技術(shù)［3－5］的基本原理是在物體外圍生成一個與物體外接的最小多面體，當(dāng)采用光線投射方法進(jìn)行繪制時，直接跳過多面體以外的數(shù)據(jù)而只在該多面體內(nèi)進(jìn)行采樣和圖像合成，則可以大大提高圖像繪制的速度。

空間剖分技術(shù)主要是利用數(shù)據(jù)空間的相關(guān)性，通過對數(shù)據(jù)空間進(jìn)行剖分，將連續(xù)的空體元劃分到一定的子區(qū)間，當(dāng)光線在數(shù)據(jù)空間穿行時，碰到這些子區(qū)間只需進(jìn)行簡單的求交運算，就可以略過整個子區(qū)間，從而減少光線投射的步數(shù)，達(dá)到降低繪制時間的目的?？臻g網(wǎng)格剖分有兩種方式:一種是將空間均勻地分割成大小相同的小立方體網(wǎng)格的方法，具有代表性的是Fujimoto等［6］提出的ARTS(Accelrated Ray Tracing System)方法和Yagel等［7］提出的離散光線投射(Discrete Ray Tracing)方法。另外一種網(wǎng)格剖分方式是采用分層的空間剖分方法。例如二叉剖分(BSP)方法、八叉樹(Octree)方法［8］和Kd樹方法［9］等。

基于光線相關(guān)性的加速方法主要基于以下一種或幾種原理［10］:

(1)像空間相關(guān)性(pixel－space coherency);

(2)物空間相關(guān)性(object－space coherence);

(3)光線間相關(guān)性(inter－ray coherency);

(4)空間跳躍(space－leaping);

(5)序列間相關(guān)性(sequence coherence)。

以上提到的一些加速方法，主要都是基于忽略無效計算的思想減少時間消耗，應(yīng)該說都取得了很好的加速效果，但對于數(shù)據(jù)量越來越大的醫(yī)學(xué)影像三維可視化，還是很難達(dá)到實時的顯示效果。

本文基于NVIDIA顯卡通用計算架構(gòu)CUDA技術(shù)，將GPU加速技術(shù)應(yīng)用于醫(yī)學(xué)三維可視化體顯示中，三維可視化重建速度得到了大幅的提高。同時，利用MFC擴展動態(tài)庫及導(dǎo)出類技術(shù)，在架構(gòu)設(shè)計上避免了大量代碼移植工作。針對特定的GeForce 9600 GT顯卡，對GPU同時執(zhí)行的線程數(shù)與計算時間做了統(tǒng)計，找到了計算時間最少的最優(yōu)線程數(shù)。

1 體顯示

早在1984年，Tuy等［11］就提出了光線追蹤的基本思想，并實現(xiàn)了對三維物體表面的直接繪制。Levoy［12］和Sabella［13］分別從插值方法、可視化映射、明暗處理、圖像合成等方面對該方法進(jìn)行了完善和擴充，形成目前體繪制中的標(biāo)準(zhǔn)光線投射繪制方法。

光線追蹤的基本原理如圖1所示。從屏幕的每個像素點出發(fā)，沿著視線的方向投射一條光線，以一定的步長在三維數(shù)據(jù)場中穿行。在它行進(jìn)的過程中，不斷進(jìn)行重采樣及顏色合成，直到阻光度足夠大或光線已經(jīng)穿過整個體數(shù)據(jù)空間為止。當(dāng)屏幕上所有像素的光線投射過程都完成后，得到最終的顯示圖像。由于光線投射的重采樣和圖像合成是按照屏幕上每個像素上發(fā)出的光線逐個進(jìn)行的，因而屬于圖像空間掃描的體繪制方法。

圖1 光線追蹤步長Fig.1 Ray tracing steps

基于光線追蹤的體顯示圖像合成公式由下式給出［14－15］:

式中:cout，λ(u，v)為光線離開體元時的顏色值; cin，λ(uv)為光線進(jìn)入體元前的顏色值;α(xi，yj，zk)為當(dāng)前體元的阻光度;cλ(xi，yj，zk)為當(dāng)前體元的顏色值。

從屏幕像素cλ(um，vn)出發(fā)，沿光線方向?qū)ψ韫舛燃邦伾至坎蓸雍铣桑罱K得到屏幕像素的顏色值，由下式給出:

式中:cλ(xi，yj，z0)=cbkg，λ;α(xi，yj，z0)=1。

2 程序?qū)崿F(xiàn)

由式(2)可知，合成屏幕像素的光線追蹤過程對于每個象素都是獨立的，每個像素的光線追蹤過程都可以由GPU中的一個thread完成，具有非常好的并行計算特性，非常適用GPU加速。

2.1 GPU加速技術(shù)介紹

由于現(xiàn)代的顯示芯片具有高度的可程序化能力及相當(dāng)高的內(nèi)存帶寬和大量的執(zhí)行單元，可以利用顯示芯片幫助進(jìn)行一些計算工作，即GPGPU。CUDA是NVIDIA的GPGPU模型，它以C語言為基礎(chǔ)，可以直接使用C語言編寫可在顯示芯片上執(zhí)行的程序。

在CUDA的架構(gòu)下，一個程序分為兩個部分:host端和device端。host端是指在CPU上執(zhí)行的部分，而device端則是在顯示芯片上執(zhí)行的部分。CUDA架構(gòu)下，顯示芯片執(zhí)行時的最小單位是 thread。數(shù)個 thread可以組成一個 block。執(zhí)行相同程序的block，可以組成grid。

2.2 程序架構(gòu)設(shè)計

VS6沒有集成CUDA程序編譯，為此我們創(chuàng)建了一個接口導(dǎo)出類CGPUSpeedUI，在VS2008下，將其編譯成相應(yīng)的lib文件。原產(chǎn)品VS6工程通過該lib文件進(jìn)行編譯。這樣 device端的kernel程序可以在VS2008下的動態(tài)庫執(zhí)行，不需要對原有工程做大量的代碼移植工作，GPU加速架構(gòu)設(shè)計如圖2所示。

圖2 GPU加速架構(gòu)設(shè)計Fig.2 GPU speedup architecture

2.3 數(shù)據(jù)交換

CUDA中，GPU不能直接存取主內(nèi)存，只能存取顯卡上的顯示內(nèi)存。因此，需要將數(shù)據(jù)從主內(nèi)存先復(fù)制到顯卡內(nèi)存中，進(jìn)行運算后，再將結(jié)果從顯卡內(nèi)存中復(fù)制到主內(nèi)存中。為此，需要在 host端CPU程序中將劑量計算需要的數(shù)據(jù)準(zhǔn)備好，然后通過調(diào)用接口類CGPUSpeedUI成員函數(shù)將數(shù)據(jù)復(fù)制到顯卡內(nèi)存中。device端kernel程序執(zhí)行結(jié)束后再通過調(diào)用接口類 CGPUSpeedUI成員函數(shù)將計算好的數(shù)據(jù)復(fù)制到同一主內(nèi)存中。

這些被復(fù)制到顯卡內(nèi)存中的數(shù)據(jù)包括在device端 kernel程序需要定義的一些數(shù)組及變量。在CUDA內(nèi)存管理中，紋理內(nèi)存具有只讀屬性，全局內(nèi)存可以進(jìn)行讀寫操作。但紋理內(nèi)存空間具有高速緩存，紋理拾取僅在高速緩存未命中時，才會從設(shè)備內(nèi)存中讀取數(shù)據(jù)，所以紋理內(nèi)存與全局內(nèi)存相比，具有更高的讀取效率。針對以上特性，我們對只讀數(shù)組采用紋理內(nèi)存方式管理，對于可讀寫數(shù)組聲明為全局變量

2.4 程序設(shè)計

圖1中屏幕像素經(jīng)過離散處理后，由式(2)可知，對于屏幕任意離散點對應(yīng)像素顏色值Cλ(μm，vn)的計算可以在顯示芯片執(zhí)行最小單位thread中獨立進(jìn)行。體顯示流程圖如圖3所示。

圖3 體顯示GPU計算流程圖Fig.3 Volume rendering GPU computing process flowchart

3 結(jié)果與比較

在試驗過程中，我們選擇了五例臨床患者CT序列影像數(shù)據(jù)進(jìn)行試驗。程序運行軟件環(huán)境為Windows XP操作系統(tǒng)。硬件配置為臺式機，CPU配置為Intel(R)Core(TM)2 Duo E8400@ 3.00 Hz，內(nèi)存2 GB。GPU配置為GeForce 9600 GT。

3.1 GPU線程優(yōu)化

在試驗過程中，對GPU同時執(zhí)行的線程數(shù)與計算時間做了統(tǒng)計，找到了計算時間最少的最優(yōu)線程數(shù)。具體如圖4所示數(shù)據(jù)曲線，從圖4中可以看出，針對GeForce 9600 GT.顯卡，同時并行執(zhí)行200個左右的threads時，計算時間最短。且當(dāng)threads個數(shù)大于220后，計算時間會有明顯增加，主要原因是同時進(jìn)行的 thread數(shù)目過多，導(dǎo)致一部分的數(shù)據(jù)儲存在顯卡內(nèi)存中，降低了執(zhí)行效率。

圖4 并行執(zhí)行threads個數(shù)與計算時間關(guān)系Fig.4 Threads number and computing time relationship

3.2 CPU與GPU計算時間比較

在試驗過程中，對CPU計算時間和GPU計算時間也做了對比。對比結(jié)果如表1所示。

表1 CPU與GPU計算時間比較Table 1 CPU and GPU Com putation Time Comparison

從表1數(shù)據(jù)可以得出，GeForce 9600 GT顯卡計算速度可提高6倍左右，使三維可視化重建時間在1 s以內(nèi)。

4 結(jié)束語

本文提出的基于GPU加速技術(shù)也可以用于其他科學(xué)計算領(lǐng)域。而且考慮到VS6沒有集成CUDA的編譯環(huán)境，本文采用2.2節(jié)的程序架構(gòu)設(shè)計，可以最大限度地減少代碼重構(gòu)的工作量，只需要將GPU device端運行的代碼移植到更高版本的VS版本即可，這對在VS6環(huán)境開發(fā)的軟件改造非常重要。同時，對于相對配置較低GeForce 9600 GT顯卡即可達(dá)到6倍左右的加速效果，如果在更高配置的顯卡上運行，將會得到更好的加速效果，具有非常高的應(yīng)用價值。

［1］Hu Ying，Hou Yue，Xu Xin－h(huán)e.Fast volume rendering formedical image［C］//Proceedings of XI International Congress for Stereology，Beijin，2003.

［2］Szirmay－Kalos L，Havran V，Balazs B，et al.On the efficiency of ray－shooting acceleration schemes［C］//Proceedings of SCCG'02 conference，2002:89－98.

［3］Chang A Y.A survey ofgeometric data structures for ray tracing［D］.Brooklyn USA:Polytechnic University，2001:12－75.

［4］Clark JH.Hierarchical geometricmodels for visible surface algorithm［J］.Communication of the ACM，1976，19(10):547－554.

［5］彭群生，鮑虎軍，金小剛.計算機真實感圖形的算法基礎(chǔ)［M］.北京:科學(xué)出版社，1999.

［6］Fujimoto A，Tanata T，Iwata K.ARTS:Accelerated raytracing system［J］.IEEEComputer Graphics and Applications，1986，6(4):16－26.

［7］Yagel R，Cohen D，Kaufman A.Discrete ray tracing［J］.IEEE Computer Graphics and Applications，1992，12(5):19－28.

［8］Fuchs H，Kedem Z M，Naylor B F.On visible surface generation by a priori tree structures［J］.Computer Graphics，1980，14(3):124－133.

［9］Glassner A S.Space subdivision for fast ray tracing［J］. IEEE Computer Graphics and applications，1984，4 (10):15－22.

［10］Bentley J L.Multidimensional binary search trees used for associative searching［C］//Communications of the ACM，1975，18(9):509－517.

［11］Tuy H K，Tuy L T.Direct 2－D display of 3－D objects［J］.IEEE Computer Graphics and Applications，1984，4(10):29－33.

［12］Levoy M.Display of surfaces from volume data［J］.IEEE Computer Graphics and Applications，1988，8(5):29－37.

［13］Sabella P.A Rendering algorithm for visualizing 3D scalar field［J］.Computer Graphics，1988，22(4):51－58.

［14］Westover L.Footprint evaluation for volume rendering［J］.Computer Graphics，1990，24(4):367－376.

［15］Frieder G，Gordon D，Reynolds R A.Back－to－front display of voxel－based objects［J］.IEEE Computer Graphics and Applications，1985，5(1):52－60.