王 玉,王 宏
(1.東北大學(xué) 中荷生物醫(yī)學(xué)與信息工程學(xué)院,沈陽(yáng) 110819;2.東北大學(xué) 機(jī)械工程與自動(dòng)化學(xué)院,沈陽(yáng)110004)
近些年來(lái),隨著CT等影像設(shè)備的不斷發(fā)展和升級(jí)換代,使得患者影像數(shù)據(jù)量越來(lái)越大,傳統(tǒng)的二維顯示方式已無(wú)法滿足實(shí)際臨床需要。醫(yī)學(xué)三維可視化技術(shù)近些年來(lái)成為被關(guān)注的領(lǐng)域。
臨床應(yīng)用中除了對(duì)三維可視化重建圖像質(zhì)量提出較高的要求外,對(duì)交互繪制速度也提出了較高的要求。基于光線投射方法可以生成高質(zhì)量的三維圖像,能夠滿足臨床診斷的需要。但它的一個(gè)突出缺點(diǎn)就是繪制速度慢,如果沒(méi)有一種有效的加速方法對(duì)其繪制速度進(jìn)行改善,其很難適應(yīng)實(shí)際應(yīng)用和醫(yī)學(xué)圖像技術(shù)發(fā)展的要求。
光線投射方法提出后,相應(yīng)的加速技術(shù)的研究就從來(lái)沒(méi)有停止過(guò)。目前,針對(duì)光線投射的軟件加速技術(shù)主要有以下幾種[1-2]:包圍盒技術(shù)(bounding box)、空間剖分(space subdivision)技術(shù)和光線相關(guān)性(ray coherence)技術(shù)。
包圍盒技術(shù)[3-5]的基本原理是在物體外圍生成一個(gè)與物體外接的最小多面體,當(dāng)采用光線投射方法進(jìn)行繪制時(shí),直接跳過(guò)多面體以外的數(shù)據(jù)而只在該多面體內(nèi)進(jìn)行采樣和圖像合成,則可以大大提高圖像繪制的速度。
空間剖分技術(shù)主要是利用數(shù)據(jù)空間的相關(guān)性,通過(guò)對(duì)數(shù)據(jù)空間進(jìn)行剖分,將連續(xù)的空體元?jiǎng)澐值揭欢ǖ淖訁^(qū)間,當(dāng)光線在數(shù)據(jù)空間穿行時(shí),碰到這些子區(qū)間只需進(jìn)行簡(jiǎn)單的求交運(yùn)算,就可以略過(guò)整個(gè)子區(qū)間,從而減少光線投射的步數(shù),達(dá)到降低繪制時(shí)間的目的。空間網(wǎng)格剖分有兩種方式:一種是將空間均勻地分割成大小相同的小立方體網(wǎng)格的方法,具有代表性的是Fujimoto等[6]提出的ARTS(Accelrated Ray Tracing System)方法和Yagel等[7]提出的離散光線投射(Discrete Ray Tracing)方法。另外一種網(wǎng)格剖分方式是采用分層的空間剖分方法。例如二叉剖分(BSP)方法、八叉樹(shù)(Octree)方法[8]和Kd樹(shù)方法[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)。
以上提到的一些加速方法,主要都是基于忽略無(wú)效計(jì)算的思想減少時(shí)間消耗,應(yīng)該說(shuō)都取得了很好的加速效果,但對(duì)于數(shù)據(jù)量越來(lái)越大的醫(yī)學(xué)影像三維可視化,還是很難達(dá)到實(shí)時(shí)的顯示效果。
本文基于NVIDIA顯卡通用計(jì)算架構(gòu)CUDA技術(shù),將GPU加速技術(shù)應(yīng)用于醫(yī)學(xué)三維可視化體顯示中,三維可視化重建速度得到了大幅的提高。同時(shí),利用MFC擴(kuò)展動(dòng)態(tài)庫(kù)及導(dǎo)出類(lèi)技術(shù),在架構(gòu)設(shè)計(jì)上避免了大量代碼移植工作。針對(duì)特定的GeForce 9600 GT顯卡,對(duì)GPU同時(shí)執(zhí)行的線程數(shù)與計(jì)算時(shí)間做了統(tǒng)計(jì),找到了計(jì)算時(shí)間最少的最優(yōu)線程數(shù)。
早在1984年,Tuy等[11]就提出了光線追蹤的基本思想,并實(shí)現(xiàn)了對(duì)三維物體表面的直接繪制。Levoy[12]和Sabella[13]分別從插值方法、可視化映射、明暗處理、圖像合成等方面對(duì)該方法進(jìn)行了完善和擴(kuò)充,形成目前體繪制中的標(biāo)準(zhǔn)光線投射繪制方法。
光線追蹤的基本原理如圖1所示。從屏幕的每個(gè)像素點(diǎn)出發(fā),沿著視線的方向投射一條光線,以一定的步長(zhǎng)在三維數(shù)據(jù)場(chǎng)中穿行。在它行進(jìn)的過(guò)程中,不斷進(jìn)行重采樣及顏色合成,直到阻光度足夠大或光線已經(jīng)穿過(guò)整個(gè)體數(shù)據(jù)空間為止。當(dāng)屏幕上所有像素的光線投射過(guò)程都完成后,得到最終的顯示圖像。由于光線投射的重采樣和圖像合成是按照屏幕上每個(gè)像素上發(fā)出的光線逐個(gè)進(jìn)行的,因而屬于圖像空間掃描的體繪制方法。
圖1 光線追蹤步長(zhǎng)Fig.1 Ray tracing steps
基于光線追蹤的體顯示圖像合成公式由下式給出[14-15]:
式中:cout,λ(u,v)為光線離開(kāi)體元時(shí)的顏色值; cin,λ(uv)為光線進(jìn)入體元前的顏色值;α(xi,yj,zk)為當(dāng)前體元的阻光度;cλ(xi,yj,zk)為當(dāng)前體元的顏色值。
從屏幕像素cλ(um,vn)出發(fā),沿光線方向?qū)ψ韫舛燃邦伾至坎蓸雍铣?,最終得到屏幕像素的顏色值,由下式給出:
式中:cλ(xi,yj,z0)=cbkg,λ;α(xi,yj,z0)=1。
由式(2)可知,合成屏幕像素的光線追蹤過(guò)程對(duì)于每個(gè)象素都是獨(dú)立的,每個(gè)像素的光線追蹤過(guò)程都可以由GPU中的一個(gè)thread完成,具有非常好的并行計(jì)算特性,非常適用GPU加速。
由于現(xiàn)代的顯示芯片具有高度的可程序化能力及相當(dāng)高的內(nèi)存帶寬和大量的執(zhí)行單元,可以利用顯示芯片幫助進(jìn)行一些計(jì)算工作,即GPGPU。CUDA是NVIDIA的GPGPU模型,它以C語(yǔ)言為基礎(chǔ),可以直接使用C語(yǔ)言編寫(xiě)可在顯示芯片上執(zhí)行的程序。
在CUDA的架構(gòu)下,一個(gè)程序分為兩個(gè)部分:host端和device端。host端是指在CPU上執(zhí)行的部分,而device端則是在顯示芯片上執(zhí)行的部分。CUDA架構(gòu)下,顯示芯片執(zhí)行時(shí)的最小單位是 thread。數(shù)個(gè) thread可以組成一個(gè) block。執(zhí)行相同程序的block,可以組成grid。
VS6沒(méi)有集成CUDA程序編譯,為此我們創(chuàng)建了一個(gè)接口導(dǎo)出類(lèi)CGPUSpeedUI,在VS2008下,將其編譯成相應(yīng)的lib文件。原產(chǎn)品VS6工程通過(guò)該lib文件進(jìn)行編譯。這樣 device端的kernel程序可以在VS2008下的動(dòng)態(tài)庫(kù)執(zhí)行,不需要對(duì)原有工程做大量的代碼移植工作,GPU加速架構(gòu)設(shè)計(jì)如圖2所示。
圖2 GPU加速架構(gòu)設(shè)計(jì)Fig.2 GPU speedup architecture
CUDA中,GPU不能直接存取主內(nèi)存,只能存取顯卡上的顯示內(nèi)存。因此,需要將數(shù)據(jù)從主內(nèi)存先復(fù)制到顯卡內(nèi)存中,進(jìn)行運(yùn)算后,再將結(jié)果從顯卡內(nèi)存中復(fù)制到主內(nèi)存中。為此,需要在 host端CPU程序中將劑量計(jì)算需要的數(shù)據(jù)準(zhǔn)備好,然后通過(guò)調(diào)用接口類(lèi)CGPUSpeedUI成員函數(shù)將數(shù)據(jù)復(fù)制到顯卡內(nèi)存中。device端kernel程序執(zhí)行結(jié)束后再通過(guò)調(diào)用接口類(lèi) CGPUSpeedUI成員函數(shù)將計(jì)算好的數(shù)據(jù)復(fù)制到同一主內(nèi)存中。
這些被復(fù)制到顯卡內(nèi)存中的數(shù)據(jù)包括在device端 kernel程序需要定義的一些數(shù)組及變量。在CUDA內(nèi)存管理中,紋理內(nèi)存具有只讀屬性,全局內(nèi)存可以進(jìn)行讀寫(xiě)操作。但紋理內(nèi)存空間具有高速緩存,紋理拾取僅在高速緩存未命中時(shí),才會(huì)從設(shè)備內(nèi)存中讀取數(shù)據(jù),所以紋理內(nèi)存與全局內(nèi)存相比,具有更高的讀取效率。針對(duì)以上特性,我們對(duì)只讀數(shù)組采用紋理內(nèi)存方式管理,對(duì)于可讀寫(xiě)數(shù)組聲明為全局變量
圖1中屏幕像素經(jīng)過(guò)離散處理后,由式(2)可知,對(duì)于屏幕任意離散點(diǎn)對(duì)應(yīng)像素顏色值Cλ(μm,vn)的計(jì)算可以在顯示芯片執(zhí)行最小單位thread中獨(dú)立進(jìn)行。體顯示流程圖如圖3所示。
圖3 體顯示GPU計(jì)算流程圖Fig.3 Volume rendering GPU computing process flowchart
在試驗(yàn)過(guò)程中,我們選擇了五例臨床患者CT序列影像數(shù)據(jù)進(jìn)行試驗(yàn)。程序運(yùn)行軟件環(huán)境為Windows XP操作系統(tǒng)。硬件配置為臺(tái)式機(jī),CPU配置為Intel(R)Core(TM)2 Duo E8400@ 3.00 Hz,內(nèi)存2 GB。GPU配置為GeForce 9600 GT。
在試驗(yàn)過(guò)程中,對(duì)GPU同時(shí)執(zhí)行的線程數(shù)與計(jì)算時(shí)間做了統(tǒng)計(jì),找到了計(jì)算時(shí)間最少的最優(yōu)線程數(shù)。具體如圖4所示數(shù)據(jù)曲線,從圖4中可以看出,針對(duì)GeForce 9600 GT.顯卡,同時(shí)并行執(zhí)行200個(gè)左右的threads時(shí),計(jì)算時(shí)間最短。且當(dāng)threads個(gè)數(shù)大于220后,計(jì)算時(shí)間會(huì)有明顯增加,主要原因是同時(shí)進(jìn)行的 thread數(shù)目過(guò)多,導(dǎo)致一部分的數(shù)據(jù)儲(chǔ)存在顯卡內(nèi)存中,降低了執(zhí)行效率。
圖4 并行執(zhí)行threads個(gè)數(shù)與計(jì)算時(shí)間關(guān)系Fig.4 Threads number and computing time relationship
在試驗(yàn)過(guò)程中,對(duì)CPU計(jì)算時(shí)間和GPU計(jì)算時(shí)間也做了對(duì)比。對(duì)比結(jié)果如表1所示。
表1 CPU與GPU計(jì)算時(shí)間比較Table 1 CPU and GPU Com putation Time Comparison
從表1數(shù)據(jù)可以得出,GeForce 9600 GT顯卡計(jì)算速度可提高6倍左右,使三維可視化重建時(shí)間在1 s以內(nèi)。
本文提出的基于GPU加速技術(shù)也可以用于其他科學(xué)計(jì)算領(lǐng)域。而且考慮到VS6沒(méi)有集成CUDA的編譯環(huán)境,本文采用2.2節(jié)的程序架構(gòu)設(shè)計(jì),可以最大限度地減少代碼重構(gòu)的工作量,只需要將GPU device端運(yùn)行的代碼移植到更高版本的VS版本即可,這對(duì)在VS6環(huán)境開(kāi)發(fā)的軟件改造非常重要。同時(shí),對(duì)于相對(duì)配置較低GeForce 9600 GT顯卡即可達(dá)到6倍左右的加速效果,如果在更高配置的顯卡上運(yùn)行,將會(huì)得到更好的加速效果,具有非常高的應(yīng)用價(jià)值。
[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]彭群生,鮑虎軍,金小剛.計(jì)算機(jī)真實(shí)感圖形的算法基礎(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.