多源γ 射束立體定向放射治療計(jì)劃系統(tǒng)的劑量計(jì)算加速算法研究

【作 者】曹國(guó)剛

上海世鵬實(shí)驗(yàn)室科技發(fā)展有限公司，上海，200072

1951年，瑞典著名神經(jīng)外科專家Leksell教授提出立體定向放射外科(Stereotactic Radiosurgery)的概念。1968年，Leksell教授與Larsson、Rexed等核物理專家以及放射治療學(xué)專家、解剖學(xué)家通力合作，在瑞典的Uppsala大學(xué)首先研制成功世界上第一臺(tái)伽瑪?shù)叮_了立體定向放射神經(jīng)外科的新紀(jì)元[1-2]。立體定向放射治療是指采用立體定向手段進(jìn)行的高精度放射治療的技術(shù)。與傳統(tǒng)的放射治療相比，它在增加腫瘤照射劑量的同時(shí)，分散了對(duì)健康組織的照射劑量，大大提高了腫瘤放射治療的臨床療效。經(jīng)過(guò)各國(guó)專家?guī)资甑呐?，立體定向放射治療現(xiàn)已形成了靜態(tài)式、旋轉(zhuǎn)式和弧形等中心三大類方式，出現(xiàn)了頭部伽瑪?shù)?、旋轉(zhuǎn)式頭部和體部伽瑪?shù)逗皖^部和體部X刀等多種裝置。

放射治療的質(zhì)量取決于治療時(shí)的吸收劑量分布，而劑量計(jì)算是目前能全面給出吸收劑量分布的唯一方法，因此高精度地計(jì)算人體的三維吸收劑量分布就成為評(píng)價(jià)治療計(jì)劃優(yōu)劣的先決條件。放射治療計(jì)劃系統(tǒng)中，劑量計(jì)算總是在計(jì)算速度和計(jì)算精度間尋找平衡點(diǎn)?；贛onte Carlo粒子輸運(yùn)理論的劑量計(jì)算是目前最精確的算法，但計(jì)算所需的時(shí)間是驚人的，極少能用于臨床[3-4]；基于Collapsed Cone卷積疊加模型的劑量計(jì)算方法，常用于X射線計(jì)劃系統(tǒng)[5-6]；而查表插值法等近似模型則適合 γ - 射線的性質(zhì)及照射方法[7]。

本文以陀螺刀（一種陀螺旋轉(zhuǎn)式多源 γ 射束立體定向放射治療系統(tǒng)）治療計(jì)劃系統(tǒng)為例，以查表插值法近似模型為基礎(chǔ)，介紹一種基于錐形坐標(biāo)系和多核并行化的劑量計(jì)算加速算法。該治療計(jì)劃系統(tǒng)由上海世鵬實(shí)驗(yàn)室科技發(fā)展有限公司自主研發(fā)，已通過(guò)國(guó)家相關(guān)機(jī)構(gòu)檢驗(yàn)。

1 多源 γ 射束立體定向放射治療計(jì)劃系統(tǒng)

多源 γ 射束立體定向放射治療系統(tǒng)(Multi-Source Stereotactic Radiotherapy System with Gamma Beam)，俗稱伽瑪?shù)?Gamma Knife)，采用幾何聚焦的原理，利用立體定向裝置、CT模擬定位、計(jì)算機(jī)立體成像與三維治療計(jì)劃系統(tǒng)等技術(shù)手段，將由多個(gè)鈷-60放射源經(jīng)準(zhǔn)直器形成的很細(xì)的γ-射線在三維空間聚焦于一點(diǎn)，要求該點(diǎn)與病灶完全重合，形成一個(gè)高劑量分布區(qū)，稱之為靶點(diǎn)(shot)。由于靶點(diǎn)上的能量極高，可以在短時(shí)間內(nèi)將靶點(diǎn)擊毀；而靶點(diǎn)之外的健康組織所受的照射劑量卻很小，從而達(dá)到了比手術(shù)切除更好的治療效果[2,8]。

治療計(jì)劃系統(tǒng)（treatment planning system,TPS）是伽瑪?shù)兜闹匾M成部分，是一種融合了諸如圖像處理和分析、數(shù)值計(jì)算和可視化、計(jì)算機(jī)圖形學(xué)、軟件工程、放射物理學(xué)以及腫瘤學(xué)等多種學(xué)科的復(fù)雜軟件系統(tǒng)，醫(yī)生藉此可以方便迅速地制定出合理的治療計(jì)劃方案。治療計(jì)劃系統(tǒng)一般由以下幾個(gè)功能模塊組成：劑量學(xué)數(shù)據(jù)管理、病例數(shù)據(jù)管理、病例數(shù)據(jù)導(dǎo)入和導(dǎo)出、框架定位、圖像配準(zhǔn)與融合、輪廓勾畫、虛擬人體建模、計(jì)劃規(guī)劃與優(yōu)化、劑量計(jì)算與計(jì)劃評(píng)估、治療計(jì)劃輸出與存檔等[9]。

目前的伽瑪?shù)吨委熡?jì)劃系統(tǒng)均采用正向設(shè)計(jì)、手工規(guī)劃的方法，計(jì)劃者負(fù)責(zé)選擇靶點(diǎn)相關(guān)的參數(shù)，由計(jì)算機(jī)進(jìn)行劑量計(jì)算。這是一個(gè)反復(fù)試誤(try and error)的過(guò)程，在找到合適的方案之前要一直進(jìn)行下去。這種方法中劑量計(jì)算速度決定計(jì)劃設(shè)計(jì)時(shí)間[10]，因此提高劑量計(jì)算速度成為計(jì)劃系統(tǒng)的關(guān)鍵。

2 治療計(jì)劃系統(tǒng)中的劑量計(jì)算

所謂劑量計(jì)算，就是計(jì)算粒子或光子射線照射到人體（或體模）中能量吸收的空間分布。由于粒子或光子線與人體組織的相互作用比較復(fù)雜，射線的種類、能量和照射方式不同，其形成的劑量分布會(huì)有很大的差別。在伽瑪?shù)吨委熡?jì)劃系統(tǒng)中，根據(jù) γ - 射線的性質(zhì)及照射方法，采用劑量計(jì)算近似模型，用實(shí)際測(cè)量的劑量學(xué)數(shù)據(jù)進(jìn)行查表插值的方法計(jì)算劑量分布的近似估計(jì)。該方法主要包括單源劑量計(jì)算、多源劑量分布疊加和多靶點(diǎn)劑量分布疊加。

2.1 單源劑量計(jì)算方法

單個(gè)靜態(tài) γ - 射線源對(duì)病人體內(nèi)某一點(diǎn)貢獻(xiàn)的劑量可以用下式表示：

式中：D(x,y,z,p,n)為當(dāng)源在位置p，使用第n個(gè)準(zhǔn)直器照射時(shí)，體內(nèi)某點(diǎn)(x,y,z)的劑量值；x,y,z為病人體內(nèi)某點(diǎn)的坐標(biāo)；p為源位置；n為準(zhǔn)直器編號(hào)；DF(n)為第n個(gè)準(zhǔn)直器的標(biāo)定劑量率；TMR(deff,n)為第n個(gè)準(zhǔn)直器，深度為deff時(shí)的組織最大劑量比(tissue maximum ratio,TMR)；deff為計(jì)算點(diǎn)的有效深度(包含非均勻組織的電子密度校正）；OAR(d,r,n)為第n個(gè)準(zhǔn)直器，離軸距離為r，深度為d的離軸比(off axis ratio,為深度校正幾何因子；SAD為源軸距，即源至等中心的距離；△d為計(jì)算點(diǎn)在中心軸上投影距離等中心的距離，遠(yuǎn)離源方向時(shí)取正，反之取負(fù)。

2.2 多源劑量計(jì)算方法

對(duì)于多個(gè)放射源或放射源為運(yùn)動(dòng)的情況，采用離散累加的方式計(jì)算劑量，所以計(jì)算點(diǎn)的實(shí)際劑量是多個(gè)放射源位置的劑量貢獻(xiàn)之和：

其中，D(x,y,z,n)為使用第n個(gè)準(zhǔn)直器照射時(shí)，體內(nèi)某點(diǎn)(x,y,z)的劑量值；Wp是源在位置p的權(quán)重。

2.3 多靶點(diǎn)劑量計(jì)算方法

如果有m個(gè)靶點(diǎn)，病人體內(nèi)的劑量是多個(gè)靶點(diǎn)的加權(quán)和，即：

其中，D(x,y,z)為體內(nèi)某點(diǎn)(x,y,z)的劑量值，wi是第i個(gè)靶點(diǎn)的權(quán)重，ni為第i個(gè)靶點(diǎn)的準(zhǔn)直器編號(hào)。

3 劑量計(jì)算加速算法

傳統(tǒng)劑量計(jì)算方法通常在直角坐標(biāo)系中進(jìn)行，不同計(jì)算點(diǎn)需重復(fù)計(jì)算深度、有效深度、離軸距離等中間結(jié)果多次，且這些中間結(jié)果無(wú)法重復(fù)使用，導(dǎo)致計(jì)算速度慢。因此，我們考慮建立一個(gè)“錐形”坐標(biāo)系，重復(fù)利用中間結(jié)果，提高計(jì)算速度。此外，劑量計(jì)算過(guò)程中，需多次計(jì)算源在不同位置的劑量分布，是一個(gè)典型的計(jì)算密集型過(guò)程，十分耗時(shí)。隨著多核處理器的可用性和性價(jià)比的不斷提高，越來(lái)越多的應(yīng)用將采用多核架構(gòu)，多核并行化的劑量計(jì)算的研究更具有實(shí)用價(jià)值。

3.1 基于錐形坐標(biāo)系的劑量計(jì)算

首先建立錐形坐標(biāo)系。通過(guò)空間幾何變換將機(jī)器坐標(biāo)系變換到以源為原點(diǎn)的直角坐標(biāo)系(x,y,z)，錐形坐標(biāo)系(d,φ,θ)與上述直角坐標(biāo)系關(guān)系如圖1。兩個(gè)坐標(biāo)系的轉(zhuǎn)換公式如下：

其中，當(dāng)z=0時(shí)，θ=0，φ=0；當(dāng)x=0，y>0時(shí),θ=90；當(dāng)x=0，y<0時(shí)，θ=270。

在錐形坐標(biāo)系中進(jìn)行劑量計(jì)算時(shí)，可以減少對(duì)deff、d、r和△d等中間結(jié)果的大量重復(fù)計(jì)算，加快計(jì)算速度。在該坐標(biāo)系下進(jìn)行劑量計(jì)算后，只需將其結(jié)果轉(zhuǎn)換回直角坐標(biāo)系下即可。

圖1 錐形坐標(biāo)系與直角坐標(biāo)系Fig.1 Cone coordingate system and cartesian coordinate system

3.2 劑量計(jì)算的多核并行化

多核技術(shù)是近年來(lái)全球計(jì)算機(jī)技術(shù)發(fā)展的重要內(nèi)容，自從英特爾在2006年底發(fā)布了全球第一款主流四核處理器至強(qiáng)5300后，計(jì)算機(jī)行業(yè)宣告正式進(jìn)入多核時(shí)代。多核計(jì)算將成為一種普及的計(jì)算模式，并將深刻影響企業(yè)和消費(fèi)者的使用模式[11-12]。為了提高劑量計(jì)算速度，本文提出一種多核并行化的劑量計(jì)算算法，以加快劑量計(jì)算速度。

多核并行化將每個(gè)單源劑量計(jì)算分配給一個(gè)虛擬的處理器(核)，所有處理器共享輸入數(shù)據(jù)，形成數(shù)據(jù)并行的算法。每個(gè)虛擬的處理器處理對(duì)應(yīng)的計(jì)算過(guò)程，無(wú)需關(guān)心其他的計(jì)算過(guò)程，但需在對(duì)共享數(shù)據(jù)進(jìn)行讀寫時(shí)，采用同步處理保證對(duì)共享數(shù)據(jù)的互斥訪問(wèn)。

執(zhí)行共享數(shù)據(jù)的多核并行計(jì)算時(shí)，線程不是孤立的，各線程間存在依賴關(guān)系，線程間對(duì)數(shù)據(jù)的共享可能會(huì)導(dǎo)致錯(cuò)誤的行為和結(jié)果。因此，應(yīng)用多線程實(shí)現(xiàn)并行計(jì)算必需引入同步機(jī)制協(xié)調(diào)各線程，以便安全共享數(shù)據(jù)。線程的同步機(jī)制包括臨界區(qū)(Critical sections)、互斥量(Mutex)、信號(hào)量(Semaphores)等。其中，臨界區(qū)和互斥量用來(lái)串行化多個(gè)線程共享數(shù)據(jù)的訪問(wèn)，其思想是每個(gè)獨(dú)占性地訪問(wèn)一個(gè)資源的線程，可在訪問(wèn)那個(gè)資源之前鎖定臨界區(qū)，訪問(wèn)之后解除鎖定。如果另一個(gè)線程也試圖鎖定該臨界區(qū)，則該線程將被阻塞直到該臨界區(qū)空閑。信號(hào)量保存有代表可用資源數(shù)量的資源數(shù)，鎖定信號(hào)量會(huì)減少資源數(shù)，釋放信號(hào)量則增加資源數(shù)。只有在線程試圖鎖定資源數(shù)為0的信號(hào)量時(shí)，線程才會(huì)被阻塞。信號(hào)量可以用來(lái)同步化同一進(jìn)程中的線程，也可同步化不同進(jìn)程中的線程。當(dāng)系統(tǒng)中有多個(gè)同種類型的資源共享時(shí)，就可利用信號(hào)量在多個(gè)線程間分配這些資源。

基于多核處理器的并行化劑量計(jì)算是共享輸入數(shù)據(jù)的并行化算法，其并行程序與串行程序十分相似，僅在同時(shí)計(jì)算累加劑量分布數(shù)據(jù)時(shí)用臨界區(qū)保證互斥訪問(wèn)。

4 實(shí)驗(yàn)結(jié)果與分析

4.1 實(shí)驗(yàn)條件

實(shí)驗(yàn)采用直徑160 mm的有機(jī)玻璃球模，并將根據(jù)該球模產(chǎn)生的仿真CT數(shù)據(jù)輸入TPS；采用標(biāo)稱尺寸為12 mm，標(biāo)稱吸收劑量率為200cGy/min的聚焦野劑量學(xué)數(shù)據(jù)；將靶點(diǎn)置于球模中心，滾筒旋轉(zhuǎn)120o，持續(xù)照射5分鐘。采用不同方法進(jìn)行劑量計(jì)算，比較不同方法的點(diǎn)劑量誤差、面積重合率和計(jì)算速度等參數(shù)。

實(shí)驗(yàn)環(huán)境：HP xw6600圖形工作站；2顆英特爾至強(qiáng)E5405四核處理器，4G內(nèi)存；Windows XP Professional SP3。

4.2 評(píng)價(jià)方式

分別在直角坐標(biāo)系和錐形坐標(biāo)系下計(jì)算劑量分布，對(duì)比劑量分布中的點(diǎn)劑量誤差和面積重合率。用點(diǎn)劑量誤差評(píng)價(jià)某個(gè)固定點(diǎn)在兩種計(jì)算方法中絕對(duì)劑量的差異，而面積重合率評(píng)價(jià)兩種計(jì)算方法中相對(duì)劑量分布的差異。

點(diǎn)劑量絕對(duì)誤差和相對(duì)誤差定義分別如下：

其中，△D為某點(diǎn)的劑量絕對(duì)誤差，△DR為該點(diǎn)的劑量相對(duì)誤差，D1和D2分別為在直角坐標(biāo)系和錐形坐標(biāo)系下計(jì)算所得的該點(diǎn)劑量值。

面積重合率定義如下：

其中，△S為某指定等劑量線的面積重合率，S1和S2分別為在直角坐標(biāo)系和錐形坐標(biāo)系下計(jì)算所得的該等劑量線所包圍區(qū)域的面積，SC為S1和S2兩部分面積的重合部分。

4.3 不同坐標(biāo)系下計(jì)算結(jié)果與分析

根據(jù)過(guò)球模中心CT層面的劑量分布結(jié)果，統(tǒng)計(jì)點(diǎn)劑量誤差和面積重合率。兩種坐標(biāo)系下劑量計(jì)算結(jié)果的點(diǎn)劑量誤差和面積重合率如表1和表2，在兩種坐標(biāo)系下靶點(diǎn)計(jì)算量與計(jì)算時(shí)間結(jié)果如表3，表中各符號(hào)含義同4.2節(jié)。

從上述數(shù)據(jù)可以看出，兩種計(jì)算方法在點(diǎn)劑量誤差和面積重合率兩個(gè)方面的差異均非常小，可以忽略不計(jì)；且新方法可以極大提高計(jì)算速度，與原有方法相比，該方法可以節(jié)約65%的時(shí)間。

表1 點(diǎn)劑量誤差Tab.1 Error of point dose

表2 面積重合率Tab.2 Overlap rate of area

表3 靶點(diǎn)劑量和計(jì)算時(shí)間Tab.3 Point dose of shot center and calculation time

表4 多核并行計(jì)算時(shí)間與速度Tab.4 Time and speed in multi-core based parallelized dose calculation

4.4 多核并行化劑量計(jì)算結(jié)果與分析

表4列出了兩種計(jì)算方法使用不同數(shù)量核計(jì)算所需的時(shí)間和速度。其中“傳統(tǒng)”表示不采用多核并行的傳統(tǒng)方法，表中分別將兩種方法的傳統(tǒng)計(jì)算方法速度定義為1.0。

圖2為計(jì)算速度與核數(shù)量之間的關(guān)系，圖中橫坐標(biāo)為劑量計(jì)算使用核的數(shù)量，縱坐標(biāo)為計(jì)算速度，虛線、點(diǎn)虛線、實(shí)線分別代表直角坐標(biāo)系、錐形坐標(biāo)系和理想情況。

圖2 多核并行化劑量計(jì)算中速度與核數(shù)量的關(guān)系Fig.2 Relationship of speed and number of cores in multi-core based parallelized dose calculation

從上述圖表可以看出，多核并行化的劑量計(jì)算只使用一個(gè)核進(jìn)行計(jì)算時(shí)，與不采用多核并行的傳統(tǒng)計(jì)算方法耗時(shí)相當(dāng)，無(wú)明顯差異；多核并行化的劑量計(jì)算速度與所使用核的數(shù)量接近正比，但比理想情況稍差，原因是多線程運(yùn)算時(shí)，線程間切換和數(shù)據(jù)同步需要占用少量的時(shí)間。因此，多核并行化的劑量計(jì)算方法更進(jìn)一步提高了劑量計(jì)算速度，其速度提升與使用的核接近正比，當(dāng)采用8核進(jìn)行并行計(jì)算時(shí)，其計(jì)算速度可提升7-8倍。

5 小結(jié)

本文提出一種基于錐形坐標(biāo)系和多核并行化的劑量計(jì)算加速算法?；阱F形坐標(biāo)系，該方法在保證計(jì)算精度的前提下可以極大提高劑量計(jì)算速度；同時(shí)利用多核并行化計(jì)算方法，還可進(jìn)一步提高劑量計(jì)算的速度，其速度提升與使用的核接近正比。從實(shí)驗(yàn)結(jié)果和分析看出，該方法可以將計(jì)算時(shí)間從傳統(tǒng)計(jì)算方法的584.3秒減少到30.3秒，節(jié)約95%的時(shí)間，即速度提升到傳統(tǒng)方法的20倍左右。在實(shí)際應(yīng)用中，還可以通過(guò)減小計(jì)算區(qū)域和降低計(jì)算網(wǎng)格精度等方法，得到劑量計(jì)算的快算結(jié)果，在可接受范圍內(nèi)降低計(jì)算精度，進(jìn)一步提高計(jì)算速度。

[1]陳炳桓.立體定向放射[M].北京: 北京出版社,1994

[2]嚴(yán)玉龍.立體定向放射療法的治療規(guī)劃——適形放療與立體定向放射外科治療規(guī)劃的理論與技術(shù)研究[D].東南大學(xué)博士論文,1997

[3]Kawrakow I,Fippel M.Investigation of variance reduction techniques for Monte Carlo photon dose calculation using XVMC[J].Phys.Med.Biol,2000,45: 2163-2183.

[4]Kawrakow I,Fippel M.VMC++,a MC algorithm optimized for electron and photon beam dose calculations for RTP[C].Proceedings of the 22th annual EMBS international conference,2000,1490-1493.

[5]Miften M,Wiesmeyer M,Monthofer S,et al.Implement at ion of FFT convolution and multigrid superposition models in the FOCUS RTP system[J].Phys.Med.Biol.,2000,45(4): 817- 833.

[6]Zhengdong Zhou,Wei Song.A Study of Dose Calculation Based on Collapsed Cone Convolution Algorithm[C].BMEI2010,2010,1331-1335.

[7]羅浩,郭吉豐,劉涵.CT機(jī)器人伽瑪?shù)吨谢谌S區(qū)域生長(zhǎng)算法的劑量快速搜索計(jì)算研究[J].中國(guó)醫(yī)療器械雜志,2005,29(3):161-163.

[8]于金明,殷蔚伯,李寶生主編.腫瘤精確放射治療學(xué)[M].山東科學(xué)技術(shù)出版社,2004

[9]於文雪,李松毅,羅立民.鈷-60三維精確適形放射治療計(jì)劃系統(tǒng)的研制.中國(guó)醫(yī)療器械雜志[J],2007,31(6): 407-410.

[10]周正東,羅立民,舒華忠,等.放射治療計(jì)劃的優(yōu)化方法[J].中國(guó)醫(yī)療器械雜志,2007,31(6): 391-394.

[11]CAO Guogang,LUO Limin,RONG Chengcheng.Multicore-based parallelized differential evolution for medical image registration[C].MIPPR2009,2009,7497: 74972L1-74972L7.

[12]曹國(guó)剛.計(jì)算機(jī)輔助手術(shù)導(dǎo)航關(guān)鍵技術(shù)研究[D].東南大學(xué)博士論文,2010.