賴佳路 宋 瑩 周 莉 白 雪 侯 氫

1（四川大學華西醫(yī)院放療科 成都610041）2（四川大學原子核科學技術研究所輻射物理及技術教育部重點實驗室 成都610064）

如今調強技術（Intensity Modulated Radiation Therapy，IMRT）已經廣泛應用于臨床，調強計劃制定的質量和速度受限于劑量計算的精度與速度。目前，放療中光子線劑量計算算法主要有蒙特卡羅算法［1］（Monte Carlo，MC）、卷 積/積 分 算 法［2-5］（Convolution/Superposition，CS）以及有限筆形束算法［6］（Finite Size Pencil Beam，F(xiàn)SPB）。MC算法被稱為劑量計算的“金標準”，它能計算各種復雜條件下的劑量分布［7］。由于MC計算時間太長，該方法極少應用于臨床。FSPB是目前應用于臨床劑量計算的方法，其計算速度快，但不能對三維非均勻介質進行修正，當介質非均勻梯度很大時候帶來誤差較大［8］。CS算 法 模 型 來 源 于Boyer［9］、Mohan［10］和Mackie［11］等在各自研究中提出計算模型。CS算法能進行三維非均勻性修正，其計算精度僅次于MC。盡管CS算法計算速度遠遠快于MC，但該算法計算速度要完全應用于臨床較困難。1989年Ahnesj?［4］提出了筒串卷積疊加算法（Collapsed Cone Convolution，CCC），該算法對CS進行了簡化，縮短了計算時間。目前CCC算法已經應用于Pinnacle、RayStation治療計劃系統(tǒng)。由于簡化后的CCC算法計算精度稍遜于CS算法，本文針對的是原始的CS算法。目前，圖形處理器（Graphics Processing Unit，GPU）已經廣泛應用于并行計算及醫(yī)學影像處理［12］。王先良等［13］借助單顆GPU對CS算法進行了加速，與僅用CPU計算相比速度提升了60倍，單野計算時間達到1 min左右。雖然加速效果不錯，但該速度應用于調強計劃顯然不足。基于此，本文詳細研究了采用1～7塊GPU對CS算法的加速情況，探討通過CPU+多GPU異構計算模型讓CS算法應用于臨床調強計劃的可行性。

1 CS算法原理和CPU+多GPU環(huán)境搭建

1.1 CS算法原理

CS算法以點核為能量沉積核（Energy Deposit Kernel，EDK），主要涉及兩個步驟［14］。首先，計算原射線在單位質量內釋放的總能量（Total Energy Released Per Unit Mass，Terma）。再 將Terma和EDK進行卷積/積分得到最終沉積的能量。其原理如圖1所示，均勻介質中任意劑量沉積點r處沉積的能量可通過式（1）計算。

圖1 卷積/積分算法原理圖Fig.1 Schematic diagram of convolution/superposition algorithm

式中：r′為光子源到光子散射點的矢量；r0為光子源到模體表面的距離；r為光子源到能量沉積點的位置矢量；D(r)為r處的劑量；T(r′)為原射線在位移處r′處產生的Terma；keff(r-r′)為有效能量沉積核，它代表原始射線在r′處與介質發(fā)生作用，產生的次級電子在位移r處單位體積產生的能量沉積分數。Terma可通過式（2）和式（3）計算［4］得到：

式中：Ψ(Ei)為能量Ei對應的能量注量；uˉ(r)為點r處的平均線性衰減系數；u(Ei，r)為點r對應能量為Ei的線性衰減系數。

1.2 CUDA架構

計算機統(tǒng)一設備架構（Compute Unified Device Architecture，CUDA）為GPU提供了一整套通用并行計算解決方案［15］。該架構由“中央”處理器（Central Processing Unit，CPU）和“協(xié)”處理器GPU組成。CPU稱為主機端，GPU稱為設備端。在CUDA架構中，GPU上運行的函數叫內核函數（Kernel）。Kernel函數以線程網格（Grid）的形式組織，Grid為多個可以被并行執(zhí)行的線程塊（Block）集合。Block包含大小相同的線程（Thread），Thread是函數執(zhí)行的最基本單元。GPU的計算核心是流多處理器（SM），每個SM以線程束（warp）為單位執(zhí)行，同一個warp的線程執(zhí)行不同數據的相同指令，這種執(zhí)行方式能有效隱藏線程的延遲和等待問題［16］。在CUDA編程模型中，編程者可以通過三個預定義變量［17］blockDim、blockIdx和threadIdx實現(xiàn)數組元素與線程的對應，從而實現(xiàn)線程間并行計算。在多GPU架構中，主機只有一個，設備可以有多個。借助CUDA 4.0工具包，可以通過cudaSetDevice（）函數在各個GPU之間切換。

2 CPU+多GPU實現(xiàn)CS算法加速

CS算法的程序實現(xiàn)過程主要涉及6個步驟：

1）相關計算參數讀取；

2）創(chuàng)建虛擬模體；

3）計算Terma；

4）計算卷積/積分；

5）劑量插值；

6）計算所有射束產生的總劑量；。其中：步驟2）～5）這4個模塊（下文簡稱為M2、M3、M4、M5）計算數據量大、數據相關性低、數據有相同的執(zhí)行程序。該4個模塊適合用GPU進行加速。下面將對4個模塊進行詳細介紹。

2.1 創(chuàng)建虛擬模體

CS算法中的EDK來自MC方法在水中模擬得到的數據。模擬過程中，射束垂直入射到模體。因此計算過程中我們需要創(chuàng)建一個與射束中心軸方向垂直的虛擬體模，如圖2所示。

圖2 創(chuàng)建虛擬模體原理圖Fig.2 Schematic diagram of creating virtual phantom

式中：θG、θT、θC分別代表機架、治療床以及準直器角度，Matrixv→b(θG，θT，θC)為射束坐標系到人體坐標系的轉換矩陣。為了得到更精確的CT密度，需要細化虛擬模體的體元。細化后的體元多達數百萬以上。與CPU循環(huán)計算方式不同，借助GPU可以讓一個線程負責一個體元的計算。

2.2 計算Terma

計算Terma采用射線追蹤的方法［18］。射線從點源出發(fā)沿著固定的方向依次穿過多個體元，在射野平面加入射束的強度信息后計算每個體元的Terma。為了保留射束發(fā)散的物理性質，把射束看成連續(xù)的在入射表面等面積的錐形束，這樣可以得到更精確的計算結果。為了保證追蹤過程每個體元都至少有一條射線穿過，需要將射束細化，細化后的射束數目較多，且相互獨立，該過程非常適合GPU并行計算。

2.3 計算卷積/積分

卷積計算中將要計算的體元作為中心，計算所有位于積分核半徑Radcut之內對其有影響的Terma產生的劑量總和。卷積/積分中計算時間和精度受Radcut大小影響較大。參考陳炳周的研究結果［19］，本文積分核半徑選擇3 cm。計算過程中將需要計算的體元r作為追蹤的起點，沿著r′-r方向進行射線追蹤，然后得到對應等效長度的積分核，對于密度不均勻的介質，需要對其做等效厚度修正，計算等效厚度使用射線追蹤的方法，三維空間中點(x，y，z)的等效深度d(x，y，z)可以通過下面公式進行計算：

式中：ρ(x′，y′，z′)是體元(x′，y′，z′)相對于水的CT值；l(x′，y′，z′)表示射線穿過體元(x′，y′，z′)的幾何長度。CS算法需要將點核離散化［20］，為了實現(xiàn)精度與速度的平衡，本文中，頂點角φ和方位角θ都取48個。在計算積分之前，用光線追蹤方法［21］建立一個體元密度和坐標無關的查找表，在等效厚度計算的時候，只需要按角度和步數Nsteps就可以簡單地得到在x、y、z、r方向的增量，將r方向的增量Δr與體元的電子密度ρ相乘便是射線穿過體元的等效長度。該過程計算量大，非常適合GPU進行加速。

2.4 劑量插值

由于機架角度變化會改變射束坐標系，為了處理不同射束坐標系計算得到的劑量，計算得到的虛擬模體中劑量還應插值到初始模體中去，便于劑量的分析、顯示等。該過程類似于模塊2，同樣適合并行計算。

3 測試結果與分析

本文實現(xiàn)了CS算法基于CUDA架構的多GPU并行計算模型。測試平臺參數如表1所示。測試模體是一個大小為30 cm×30 cm×30 cm密度非均勻模體，體元密度值是相對于水為0～1.2之間的隨機數，體元尺寸大小設定為0.2 cm×0.2 cm×0.2 cm，總的體元數目約為337萬，該體元數目遠遠超過我們臨床中常見部位的體元數目，如頭部約160萬、胸部約239萬、腹部約250萬。虛擬模體體元的大小和初始模體體元大小設置相同。射野的大小為10 cm×10 cm，野的分辨率為0.2 cm×0.2 cm。測試機架角度為0°和45°。光子能量采用6 MV，GPU采用1～7個。本文定義了加速比TiG/C，表示第i個模塊用CPU計算花費時間與使用GPU花費的時間之比，式如（6）所示，為總的加速因子。

本文GPU計算時間包括在各個GPU上分配內存以及將數據拷貝到各個GPU上進行并行計算所用時間，由于本實驗中內核啟動極其快，所以本文忽略掉了內核啟動時間。

表2和表3分別給出了機架角度為0°和45°時，CS算法4個模塊單用CPU和采用1～7個GPU的計算時間及加速比。

多GPU并行計算的環(huán)境比傳統(tǒng)的單GPU復雜，本文測試環(huán)境中使用的是共享系統(tǒng)多GPU，即一個單獨的系統(tǒng)中包含多個GPU，這些GPU通過共享一個CPU的RAM進行連接，各個GPU數據之間不能直接傳輸，必須通過CPU內存進行中轉。隨著GPU使用個數的增加，需要在更多的GPU上為核函數所需數據分配內存，同時需要將CPU數據分別傳向更多的GPU。當多個GPU掛載到總線上后，總線傳輸壓力增大［22］，由于本文研究內容需要數據在CPU和GPU之間頻繁進行交換，數據傳輸時間延遲便成了本文研究過程中的一個關鍵問題。特別是對于在CPU上計算時間本身較短的情況。表2和表3中的M2、M3、M5三個模塊采用多個GPU計算時間說明了這種情況。該三個模塊在CPU計算時間均在1 s之內，隨著GPU使用個數的增加一方面數據傳輸總時間增加，另一方面內核啟動總時間也會相應增加，Dana Schaa與David Kaeli對共享總線結構的單GPU和多GPU的傳輸速率做過測試，結果表明單GPU的吞吐率甚至會高于多GPU［23］。針對本文研究情況，建議M2、M3、M5三塊模塊使用1～2個GPU進行計算。對于M4模塊，隨著GPU使用個數的增加加速倍數越來越高，分析原因是數據傳輸延遲和內核啟動時間與模塊整個計算時間相比完全可以忽略不計，當采用7個GPU時候加速倍數已經高達180倍以上。

表2 機架角度為0°時，2～5四個模塊計算時間和加速比Table 2 Execution time and speedup factors for 2～5 module when gantry angle was 0°

表3 機架角度為45°時，2～5四個模塊計算時間和加速比Table 3 Execution time and speedup factors for 2～5 module when gantry angle was 45°

對于模塊M4經過代碼優(yōu)化后，加速比達到了198倍以上，主要原因是研究過程中取消了M4中大量的邏輯判斷語句。該模塊采用CPU計算時，為了減少計算時間會加入一些邏輯判斷語句來減少循環(huán)次數。GPU與CPU相比并不擅長邏輯運算，考慮到GPU中擁有大量線程，因此沒必要通過邏輯判斷來減少循環(huán)次數。通過測試不難發(fā)現(xiàn)，提高CS算法計算加速比可以通過兩種方式：一是對程序代碼進行優(yōu)化，盡量減少邏輯判斷語句；二是對程序中各個模塊的GPU使用個數進行優(yōu)化。表4是優(yōu)化代碼和各個模塊GPU使用個數后計算結果，對模塊M2和M5采用1個GPU，對模塊M3采用兩個GPU，對模塊M4采用7個GPU且取消M4中邏輯判斷語句。通過優(yōu)化組合，最終加速倍數達到了207倍，計算時間縮短至9 s。上文提到臨床中頭部、胸部和腹部擁有更少的體元數目，這些部位若采用本文的方法將花費更少的計算時間，因此加速后的CS算法能夠滿足臨床要求。

表4 通過優(yōu)化代碼和GPU使用個數優(yōu)化后四個模塊總的加速比Table 4 The speedup factors for 4 modules after optimizing code and number of GPU used

由于不是所有的CUDA函數都滿足IEEE 754標準，因此有必要對結果準確性進行驗證。驗證過程中將CPU計算的劑量結果作為參考劑量，將采用不同GPU個數計算后的劑量結果作為評估劑量，對二者進行γ因子分析［24］。分析過程中距離標準選擇3 mm，百分劑量差標準選擇3%。經過評估，無論采用多少個GPU進行計算，GPU計算結果的γ通過率都是100%。說明使用不同個數GPU計算后的結果滿足精度要求。值得一提的是勾成俊等［25］借助單個GPU讓電子束劑量計算加速了98倍，在后續(xù)工作中，我們將研究多GPU在電子束劑量計算中加速情況。

4 結語

本文研究采用不同個數GPU對CS算法的加速結果表明：并行計算的加速倍數并非完全與GPU使用個數呈線性關系，還應充分考慮數據傳輸問題。對于在CPU上計算時間較短的程序代碼，僅用一兩個GPU進行加速即可。此外將CPU執(zhí)行代碼改為GPU執(zhí)行代碼過程中還應該考慮到邏輯運算問題，在編寫GPU計算代碼過程中應該盡量減少邏輯判斷語句，而不是簡單地對CPU代碼進行移植。總之，本文研究結果表明：通過程序代碼和GPU使用個數的優(yōu)化，CS算法能完全應用于臨床調強計劃。