田林琳，李 瑩

(沈陽(yáng)工學(xué)院 信息與控制學(xué)院，遼寧 撫順 113122)

1 概 述

放射治療計(jì)劃系統(tǒng)(TPS)是利用計(jì)算機(jī)軟件技術(shù)、劑量計(jì)算方法、控制手段等，解決治療過(guò)程中的精確定位、劑量計(jì)算的快速精確、保證目標(biāo)劑量的準(zhǔn)確實(shí)施以及治療計(jì)劃的優(yōu)化等問(wèn)題[1]。數(shù)字重建放射影像(DRR)是TPS系統(tǒng)中的關(guān)鍵功能，是通過(guò)模擬X光線的衰減過(guò)程，從CT斷層數(shù)據(jù)中生成類似X光片的成像效果，這對(duì)于放療過(guò)程中的質(zhì)量控制具有非常重要的意義[2]。DRR可在任意三維空間角度下實(shí)現(xiàn)對(duì)病患的“透視”模擬，能夠隨意觀察放射靶區(qū)、特定組織或器官，或靶區(qū)與周圍器官間的空間關(guān)系，因?yàn)闆](méi)有治療床的限制，可利用DRR得到模擬定位機(jī)難以拍攝的圖像。同時(shí)省略了使用模擬放射攝影圖像的成像步驟，既節(jié)省了臨床費(fèi)用，也減少了對(duì)病患的輻射[3]。

在TPS的臨床應(yīng)用中，DRR具有重要作用，為了達(dá)到臨床應(yīng)用的目標(biāo)，要求DRR能夠準(zhǔn)確、快速地反映放射治療中的真實(shí)情況，并能夠滿足動(dòng)態(tài)性、交互性的需要[4]。這也是三維圖像可視化領(lǐng)域中，人們一直試圖解決的兩個(gè)問(wèn)題：一是實(shí)時(shí)地進(jìn)行圖像繪制，在DRR的實(shí)際使用中，要求可以快速、多角度地生成圖像；二是真實(shí)地顯示結(jié)果，這要求計(jì)算結(jié)果能夠真實(shí)、客觀地反映其解剖結(jié)構(gòu)在三維空間上的細(xì)節(jié)與特征[5]。

典型的DRR算法是以光線跟蹤算法為基礎(chǔ)，利用虛擬攝影機(jī)視圖追蹤光線至投影圖像中的每一個(gè)像素，每一條光線所經(jīng)過(guò)的體元數(shù)據(jù)進(jìn)行累加，累加值與像素亮度值成正比，使用這種方法實(shí)現(xiàn)重建出的DRR圖像顯示精度和圖像質(zhì)量都非常好，滿足了顯示結(jié)果真實(shí)性的要求[6-7]。但由于此算法需要跟蹤每條從虛擬光源發(fā)出的光線，涉及到大量的光線與CT三維體數(shù)據(jù)進(jìn)行求交測(cè)試的運(yùn)算，運(yùn)算量非常大。傳統(tǒng)的基于CPU技術(shù)的實(shí)現(xiàn)方式通常不能滿足交互式DRR算法的需求，隨著GPU(graphic processing unit)在并行計(jì)算能力、存儲(chǔ)容量和可編程能力等方面的發(fā)展，使DRR算法的性能提升成為可能[8]。

英偉達(dá)(NVIDIA)公司于2007年提出的CUDA架構(gòu)(compute unified device architecture)使得開發(fā)者能夠以GPU的強(qiáng)大并行計(jì)算能力為基礎(chǔ)，建立一種更快更高效的數(shù)據(jù)計(jì)算解決方案。但CUDA架構(gòu)的一個(gè)最主要問(wèn)題是通用性，人們只能使用NVIDIA系列產(chǎn)品對(duì)GPU進(jìn)行開發(fā)。而作為一個(gè)軟件系統(tǒng)，TPS是運(yùn)行在醫(yī)院提供的計(jì)算機(jī)上，各醫(yī)院的硬件環(huán)境千差萬(wàn)別，所以使用CUDA加速的局限性比較高。而基于異構(gòu)平臺(tái)的OpenCL通用計(jì)算規(guī)范，允許開發(fā)者以相同的代碼在任何支持該規(guī)范的平臺(tái)上運(yùn)行，包括NVIDIA和AMD顯卡，甚至包括支持OpenCL的CPU，這大大降低了程序?qū)τ布脚_(tái)的依賴性，具有較強(qiáng)的理論意義和實(shí)踐價(jià)值。

使用OpenCL對(duì)TPS系統(tǒng)中的DRR算法進(jìn)行并行加速，主要出于三個(gè)目的：

(1)通過(guò)對(duì)DRR算法的加速，滿足算法交互性和實(shí)時(shí)性的需要，這不僅能減少患者的等待時(shí)間，也提高了TPS的使用效率。

(2)OpenCL作為通用計(jì)算標(biāo)準(zhǔn)，與CUDA架構(gòu)相比，其性能缺乏評(píng)判，通過(guò)使用OpenCL技術(shù)在NVIDIA顯卡上實(shí)現(xiàn)對(duì)DRR算法的加速，軟件人員可以更好地審視其性能，以便于在今后的開發(fā)中進(jìn)行平臺(tái)選擇。

(3)通過(guò)將數(shù)據(jù)存儲(chǔ)在不同的硬件存儲(chǔ)器上的性能對(duì)比，開發(fā)人員可以深入了解不同存儲(chǔ)器在并行優(yōu)化中的影響，便于在使用OpenCL技術(shù)時(shí)進(jìn)行數(shù)據(jù)的合理布局。

2 OpenCL概述

OpenCL是一個(gè)可以在異構(gòu)平臺(tái)上進(jìn)行并行編程的開放框架標(biāo)準(zhǔn)，包括一組底層程序接口和一個(gè)高效、快速、可移植的抽象層，為開發(fā)者提供了一個(gè)有效的并行開發(fā)環(huán)境、一組與平臺(tái)無(wú)關(guān)的工具和豐富的中間層。

OpenCL框架運(yùn)行時(shí)使用主機(jī)端程序?qū)λ兄С諳penCL的設(shè)備進(jìn)行統(tǒng)一管理。對(duì)于每個(gè)支持OpenCL規(guī)范的計(jì)算設(shè)備，其在內(nèi)部又進(jìn)一步劃分為多個(gè)計(jì)算單元(compute unit)，每個(gè)計(jì)算單元又劃分為多個(gè)處理單元(processing element)，每個(gè)處理單元以SIMD或SPMD的方式運(yùn)行[9]。

為了支持代碼的可移植性，OpenCL定義了一個(gè)抽象的內(nèi)存模型，開發(fā)者可以根據(jù)該模型編寫代碼，硬件廠商可以將該模型映射到他們自己實(shí)際的硬件上[10-11]。OpenCL定義的內(nèi)存空間如圖1所示，這些內(nèi)存空間與OpenCL程序密切相關(guān)。

圖1 OpenCL內(nèi)存模型

OpenCL內(nèi)存模型按照訪問(wèn)權(quán)限劃分，依次分為全局內(nèi)存和常量?jī)?nèi)存(Kernel范圍)、本地內(nèi)存(Workgroup范圍)、私有內(nèi)存(Work-item范圍)。其中紋理內(nèi)存和常量?jī)?nèi)存的數(shù)據(jù)位于全局內(nèi)存，但其擁有高速緩存[12]。

3 DRR算法的優(yōu)化過(guò)程

3.1 DRR算法原理

DRR圖像是通過(guò)模擬X光線衰減過(guò)程生成的，X光線穿過(guò)的組織不同，被吸收的劑量也不同，所以采用X光線衰減模型進(jìn)行模擬，有如下公式：

(1)

其中，I0為X射線的初始強(qiáng)度；μi為組織i的線性衰減系數(shù)；L為X射線的長(zhǎng)度。

式1的離散形式如式2所示：

I=I0×e-∑μili

(2)

其中，I0和μi同式1；li表示X射線穿過(guò)組織i的距離。式2為生成DRR圖像的關(guān)鍵公式[13]。

3.2 編寫串行程序

編寫一個(gè)能夠正確運(yùn)行、沒(méi)有存儲(chǔ)器資源泄露的CPU串行程序，作為并行優(yōu)化的正確性和性能測(cè)試的基準(zhǔn)。這樣做的原因是：(1)由于在目前的OpenCL版本上，使用C++只能實(shí)現(xiàn)一些OpenCL編程中主機(jī)部分的函數(shù)，包括OpenCL設(shè)備初始化、CPU和GPU之間的數(shù)據(jù)通信、注銷OpenCL計(jì)算環(huán)境等，而真正在GPU上執(zhí)行的kernel函數(shù)和device函數(shù)都是以C語(yǔ)言函數(shù)的形式編寫的，所以首先把TPS系統(tǒng)中DRR算法部分由C++程序轉(zhuǎn)換為C程序。(2)OpenCL device上運(yùn)行的代碼不易調(diào)試，所以在串行程序中去除不易并行的因素，使得程序中的循環(huán)或者迭代之間沒(méi)有或者較少的數(shù)據(jù)依賴，這會(huì)大大減少并行程序出錯(cuò)的概率，因此接下來(lái)在第1步串行程序的基礎(chǔ)上改寫以符合OpenCL的并行操作，例如函數(shù)中的靜態(tài)變量和全局變量在程序運(yùn)行過(guò)程中，會(huì)被多個(gè)線程同時(shí)訪問(wèn)，但其值一直保持不變，會(huì)使并行程序中線程之間存在依賴性，不利于程序的并行化，所以根據(jù)算法的實(shí)現(xiàn)過(guò)程改造為局部變量。還有將程序中代表CT體數(shù)據(jù)的數(shù)據(jù)指針由二級(jí)指針修改為一級(jí)指針，這樣CT體數(shù)據(jù)成為一個(gè)連續(xù)的內(nèi)存數(shù)據(jù)塊，方便了數(shù)據(jù)從CPU端向設(shè)備端的復(fù)制和在綁定紋理時(shí)對(duì)體數(shù)據(jù)內(nèi)存排列的要求。

3.3 確定并行化方案

支持OpenCL計(jì)算的GPU通常具有數(shù)個(gè)多處理器，每個(gè)多處理器都具有SIMD架構(gòu)，支持在同一時(shí)鐘周期內(nèi)處理不同的數(shù)據(jù)，這使得OpenCL主要針對(duì)基于數(shù)據(jù)的并行計(jì)算方式[14]。光線追蹤算法中每一條光線的計(jì)算過(guò)程都是相互獨(dú)立的，不同光線之間沒(méi)有依賴關(guān)系，因此可以方便地進(jìn)行并行化，只需將每一條光線設(shè)計(jì)成一個(gè)獨(dú)立的計(jì)算線程就能實(shí)現(xiàn)基于光線追蹤的DRR算法并行化。這是利用OpenCL實(shí)現(xiàn)基于光線追蹤的DRR算法的可能性。圖2是DRR算法的并行化方案。

圖2 DRR算法的并行化方案

使用OpenCL進(jìn)行并行程序設(shè)計(jì)，其中一個(gè)主要特點(diǎn)是針對(duì)不同的存儲(chǔ)器類型進(jìn)行數(shù)據(jù)存儲(chǔ)。GPU中有全局內(nèi)存、常量?jī)?nèi)存、紋理內(nèi)存和共享內(nèi)存。DRR算法的輸入數(shù)據(jù)CT體數(shù)據(jù)少則80～150張，多則200～400張甚至更多。以150張CT圖像為例，其需要的內(nèi)存為150*0.5 MB=75 MB，所以只有全局內(nèi)存和紋理內(nèi)存可供選擇。紋理內(nèi)存緩存在芯片上，因此它能夠減少對(duì)內(nèi)存的請(qǐng)求并提供更高效的內(nèi)存帶寬[15]。紋理緩存是專門為那些在內(nèi)存訪問(wèn)模式中存在大量空間局部性(spatial locality)的圖形應(yīng)用程序而設(shè)計(jì)的[16]。文中對(duì)CT體數(shù)據(jù)存儲(chǔ)在全局內(nèi)存和紋理內(nèi)存的性能進(jìn)行了對(duì)比。

常量?jī)?nèi)存中的數(shù)據(jù)位于全局內(nèi)存，但擁有局部緩存加速，常量?jī)?nèi)存的空間較小(只有64 K)，因此將虛擬光源的點(diǎn)坐標(biāo)、CT圖像的層厚、圖像大小、像素大小等在計(jì)算過(guò)程中恒定不變的少量參數(shù)放在常量?jī)?nèi)存中。

對(duì)于不同CT值數(shù)據(jù)的衰減表，大小是4 096，類型是單精度浮點(diǎn)，可以選擇將數(shù)據(jù)存儲(chǔ)到全局內(nèi)存、紋理內(nèi)存或者常量?jī)?nèi)存，文中也對(duì)三種情況下的程序性能進(jìn)行了對(duì)比。

3.4 并行化實(shí)現(xiàn)步驟

OpenCL實(shí)現(xiàn)DRR算法的步驟主要有：

(1)調(diào)用clGetPlatformIDs獲得可用的平臺(tái)，設(shè)置使用OpenCL設(shè)備；

(2)調(diào)用clCreateContext在指定平臺(tái)上創(chuàng)建上下文；

(3)調(diào)用clCreateCommandQueue函數(shù)，在已經(jīng)創(chuàng)建好的上下文上創(chuàng)建命令隊(duì)列；

(4)在GPU上分配三維image對(duì)象，并將CT體數(shù)據(jù)和CT衰減表從CPU讀取至設(shè)備的全局內(nèi)存或者紋理內(nèi)存中；

(5)在GPU上為算法的其他參數(shù)分配設(shè)備內(nèi)存或者常量設(shè)備內(nèi)存，如虛擬光源的點(diǎn)坐標(biāo)、CT圖像的層厚、圖像大小、像素大小等，并將CPU端的數(shù)據(jù)復(fù)制到對(duì)應(yīng)的設(shè)備內(nèi)存中；

(6)創(chuàng)建一個(gè)網(wǎng)格，并且設(shè)置好workgroup和workitem的大小；

(7)創(chuàng)建內(nèi)核對(duì)象，使用clCreateKernel創(chuàng)建一個(gè)kernel對(duì)象，調(diào)用clEnqueueNDRangeKernel啟動(dòng)kernel函數(shù)；

(8)kernel函數(shù)執(zhí)行，先判斷當(dāng)前的workitem要處理的像素是否在DRR圖像范圍內(nèi)，如果不是，則不進(jìn)行計(jì)算直接返回；否則執(zhí)行DRR算法內(nèi)核函數(shù)，在內(nèi)核函數(shù)的執(zhí)行過(guò)程中更新DRR圖像關(guān)聯(lián)的全局設(shè)備內(nèi)存；

(9)kernel函數(shù)執(zhí)行完畢，把更新后全局內(nèi)存數(shù)據(jù)復(fù)制回CPU內(nèi)存；

(10)在屏幕上描畫CPU內(nèi)存上的結(jié)果；

(11)釋放OpenCL設(shè)備上的空間以及CPU上的內(nèi)存。

3.5 并行參數(shù)設(shè)置

基于OpenCL的DRR算法中，用一個(gè)線程來(lái)計(jì)算DRR圖像中一個(gè)像素的值，如DRR圖像的分辨率為m*n，則在一個(gè)kernel中就包括使用m*n個(gè)線程，每個(gè)workgroup中的workitem設(shè)置為32*16，則一個(gè)kernel中的workgroup數(shù)目為(m/32)*(n/16)。

4 實(shí) 驗(yàn)

4.1 實(shí)驗(yàn)環(huán)境

針對(duì)相同的算法復(fù)雜度，將CT體數(shù)據(jù)分別存儲(chǔ)在全局內(nèi)存和紋理內(nèi)存并分別與串行程序進(jìn)行了對(duì)比，另外對(duì)CT衰減表分別存儲(chǔ)于全局內(nèi)存、紋理內(nèi)存和常數(shù)內(nèi)存并與串行程序進(jìn)行了對(duì)比，以評(píng)估OpenCL將數(shù)據(jù)存儲(chǔ)在不同存儲(chǔ)器上的性能。

使用的平臺(tái)如下：

GPU：NVIDIA GeForce GTX760 2 G顯存；

CPU：Intel i5 4590 4核4線程 3.3 GHz；

內(nèi)存：DDR3 1 600 MHz 8 G。

4.2 結(jié)果分析

文中分別對(duì)分辨率為2 048*2 048和1 024*1 024的DDR圖像進(jìn)行了CPU和GPU的計(jì)算結(jié)果和性能對(duì)比，使用GPU運(yùn)算的結(jié)果和CPU相比，運(yùn)行結(jié)果誤差是1e-2，結(jié)果可以接受，如圖3所示。

圖3 CPU和GPU的運(yùn)算結(jié)果對(duì)比

性能對(duì)比如表1～3所示。

表1 體數(shù)據(jù)在全局內(nèi)存時(shí)的性能對(duì)比

表2 體數(shù)據(jù)在紋理內(nèi)存時(shí)的性能對(duì)比

表3 CT衰減表在不同存儲(chǔ)器上的性能對(duì)比

由表1～3中的對(duì)比數(shù)據(jù)可見：經(jīng)過(guò)OpenCL優(yōu)化的DRR算法速度比CPU串行版本明顯快很多；在OpenCL編程中使用紋理內(nèi)存可以大大提高程序的加速比；在相同條件下，紋理內(nèi)存的速度略大于常數(shù)內(nèi)存，并且遠(yuǎn)大于全局內(nèi)存；OpenCL并行程序和串行程序的加速比與基于CUDA的DRR圖像生成算法相近[13]。

5 結(jié)束語(yǔ)

通過(guò)利用DRR算法的可并行性，使用OpenCL技術(shù)對(duì)基于光線追蹤的DRR算法進(jìn)行優(yōu)化，并在GPU上進(jìn)行了并行算法調(diào)優(yōu)；通過(guò)對(duì)OpenCL中的各種存儲(chǔ)器的優(yōu)化使用，顯著提高了程序的運(yùn)行效率。目前，OpenCL的使用仍然有很多局限，在諸如復(fù)雜的分支和邏輯判斷等方面仍然不如CPU靈活，但隨著GPU的快速發(fā)展，相信OpenCL的應(yīng)用范圍將更加廣闊。

參考文獻(xiàn):

[1] 閆 鋒.數(shù)字重建放射影像生成方法及其應(yīng)用研究[D].合肥：合肥工業(yè)大學(xué)，2010.

[2] 戴建榮，胡逸民.圖像引導(dǎo)放療的實(shí)現(xiàn)方式[J].中華放射腫瘤學(xué)雜志，2006,15(2):132-135.

[3] 吳宜燦，李國(guó)麗，陶聲祥，等.精確放射治療系統(tǒng)ARTS的研究與發(fā)展[J].中國(guó)醫(yī)學(xué)物理學(xué)雜志，2005,22(6):683-690.

[4] FERNANDO R,KILGARD M J.The Cg tutorial:the definitive guide to programmable real-time graphics[M].[s.l.]:Addison-Wesley,2003.

[5] 汪 俊，吳章文，張從華，等.基于射線追蹤的數(shù)字重建影像增強(qiáng)技術(shù)[J].生物醫(yī)學(xué)工程學(xué)雜志,2005,22(1):125-128.

[6] HUANG Jianbing,CARTER M B.Interactive transparency rendering for large CAD models[J].IEEE Transactions on Visualization and Computer Graphics,2005,11(5):584-595.

[7] LEVORY M.Volume rendering by adaptive refinement[J].Visual Computer,1990,6(1):2-7.

[8] RUSSAKOFF D B,ROHLFING T,MORI K,et al.Fast generation of digitally reconstructed radiographs using attenuation fields with application to 2D-3D image registration[J].IEEE Transactions on Medical Imaging,2005,24(11):1441-1454.

[9] 李 森，李新亮，王 龍，等.基于OpenCL的并行方腔流加速性能分析[J].計(jì)算機(jī)應(yīng)用研究,2011,28(4):1401-1403.

[10] 賈海鵬，張?jiān)迫垏?guó)平，等.基于OpenCL的拉普拉斯圖像增強(qiáng)算法優(yōu)化研究[J].計(jì)算機(jī)科學(xué),2012,39(5):271-277.

[11] NVIDIA Corporation. NVIDIA OpenCL JumpStart guide[M]．America:NVIDIA,2009．

[12] MUNSHI A．The OpenCL specification[S]．America:Khronos OpenCL Working Group，2009.

[13] 杜曉剛，黨建武，王陽(yáng)萍.基于CUDA的數(shù)字重建影像生成算法[J].計(jì)算機(jī)科學(xué),2015,42(2):301-305.

[14] 劉 磊.基于GPU的醫(yī)學(xué)圖像三維重建及可視化技術(shù)研究[D].廣州:南方醫(yī)科大學(xué),2008.

[15] 劉 鵬，高 軍，雷勛祖,等.一種基于光場(chǎng)的數(shù)字重建影像快速生成算法[J]. 南方醫(yī)科大學(xué)學(xué)報(bào),2007,27(10):1537-1539.

[16] NAGY Z,KLEIN R.Depth-peeling for texture-based volume rendering[C]//Proceedings of the 11th Pacific conference on computer graphics and applications.Washington,DC,USA:IEEE Computer Society,2003:429-433.