GPU加速技術在醫學三維可視化中的應用

王 玉，王 宏

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

近些年來，隨著CT等影像設備的不斷發展和升級換代，使得患者影像數據量越來越大，傳統的二維顯示方式已無法滿足實際臨床需要。醫學三維可視化技術近些年來成為被關注的領域。

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

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

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

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

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

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

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

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

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

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

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

本文基于NVIDIA顯卡通用計算架構CUDA技術，將GPU加速技術應用于醫學三維可視化體顯示中，三維可視化重建速度得到了大幅的提高。同時，利用MFC擴展動態庫及導出類技術，在架構設計上避免了大量代碼移植工作。針對特定的GeForce 9600 GT顯卡，對GPU同時執行的線程數與計算時間做了統計，找到了計算時間最少的最優線程數。

1 體顯示

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

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

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

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

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

從屏幕像素cλ(um，vn)出發，沿光線方向對阻光度及顏色分量采樣合成，最終得到屏幕像素的顏色值，由下式給出:

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

2 程序實現

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

2.1 GPU加速技術介紹

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

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

2.2 程序架構設計

VS6沒有集成CUDA程序編譯，為此我們創建了一個接口導出類CGPUSpeedUI，在VS2008下，將其編譯成相應的lib文件。原產品VS6工程通過該lib文件進行編譯。這樣 device端的kernel程序可以在VS2008下的動態庫執行，不需要對原有工程做大量的代碼移植工作，GPU加速架構設計如圖2所示。

圖2 GPU加速架構設計Fig.2 GPU speedup architecture

2.3 數據交換

CUDA中，GPU不能直接存取主內存，只能存取顯卡上的顯示內存。因此，需要將數據從主內存先復制到顯卡內存中，進行運算后，再將結果從顯卡內存中復制到主內存中。為此，需要在 host端CPU程序中將劑量計算需要的數據準備好，然后通過調用接口類CGPUSpeedUI成員函數將數據復制到顯卡內存中。device端kernel程序執行結束后再通過調用接口類 CGPUSpeedUI成員函數將計算好的數據復制到同一主內存中。

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

2.4 程序設計

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

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

3 結果與比較

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

3.1 GPU線程優化

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

圖4 并行執行threads個數與計算時間關系Fig.4 Threads number and computing time relationship

3.2 CPU與GPU計算時間比較

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

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

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

4 結束語

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

［1］Hu Ying，Hou Yue，Xu Xin－he.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］彭群生，鮑虎軍，金小剛.計算機真實感圖形的算法基礎［M］.北京:科學出版社，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.