多GPU異構(gòu)模型實(shí)現(xiàn)放射治療中卷積/積分算法的快速計(jì)算

賴佳路 宋 瑩 周 莉 白 雪 侯 氫

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

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

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

1.1 CS算法原理

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

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

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

式中：Ψ(Ei)為能量Ei對(duì)應(yīng)的能量注量；uˉ(r)為點(diǎn)r處的平均線性衰減系數(shù)；u(Ei，r)為點(diǎn)r對(duì)應(yīng)能量為Ei的線性衰減系數(shù)。

1.2 CUDA架構(gòu)

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

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

CS算法的程序?qū)崿F(xiàn)過程主要涉及6個(gè)步驟：

1）相關(guān)計(jì)算參數(shù)讀??；

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

3）計(jì)算Terma；

4）計(jì)算卷積/積分；

5）劑量插值；

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

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

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

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

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

2.2 計(jì)算Terma

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

2.3 計(jì)算卷積/積分

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

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

2.4 劑量插值

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

3 測(cè)試結(jié)果與分析

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

本文GPU計(jì)算時(shí)間包括在各個(gè)GPU上分配內(nèi)存以及將數(shù)據(jù)拷貝到各個(gè)GPU上進(jìn)行并行計(jì)算所用時(shí)間，由于本實(shí)驗(yàn)中內(nèi)核啟動(dòng)極其快，所以本文忽略掉了內(nèi)核啟動(dòng)時(shí)間。

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

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

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

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

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

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

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

4 結(jié)語

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