周正東 宋 威

(南京航空航天大學(xué)核科學(xué)與工程系，南京 210016)

引言

放射治療是目前治療惡性腫瘤的主要手段之一，隨著計(jì)算機(jī)、醫(yī)學(xué)影像、放射物理以及放射生物等技術(shù)的發(fā)展，以X射線為放射源的適形和調(diào)強(qiáng)放射治療技術(shù)迅速發(fā)展。放射治療計(jì)劃系統(tǒng)作為實(shí)現(xiàn)放射治療精確定位、精確計(jì)劃、精確治療目標(biāo)的必要手段具有重要的意義。放射治療計(jì)劃的目的是盡可能選擇最好的照射方案，以達(dá)到對靶區(qū)施行高劑量照射，在產(chǎn)生不可恢復(fù)性摧毀的同時(shí)，使周圍正常組織及關(guān)鍵或敏感組織所受影響最小。劑量計(jì)算是放射治療計(jì)劃系統(tǒng)的核心，collapsed cone卷積/疊加劑量計(jì)算方法[1]是一種基于點(diǎn)核模型的劑量計(jì)算方法，該方法較常用的筆形束劑量計(jì)算方法計(jì)算精度高，但計(jì)算速度相對較慢，自從該方法提出后，有關(guān)此算法的研究一直受到眾多學(xué)者的關(guān)注[2-8]。我們對collapsed cone卷積/疊加劑量計(jì)算方法進(jìn)行了深入研究，在此基礎(chǔ)上開發(fā)了一套光子束放射治療計(jì)劃系統(tǒng)，并運(yùn)用蒙特卡羅方法對劑量計(jì)算的準(zhǔn)確性進(jìn)行了驗(yàn)證。

1 Collapsed Cone劑量計(jì)算方法

卷積/疊加算法是目前放射治療計(jì)劃系統(tǒng)常用的光子束劑量計(jì)算方法，它是一種基于核模型的算法。常用的卷積核包括點(diǎn)核(point kernel)，筆形束核(pencil beam kernel)、面核(planar kernel)等。與筆形束核模型相比較，點(diǎn)核模型可以處理電子失衡區(qū)域的劑量計(jì)算問題[9]，計(jì)算精度更高，不足之處是計(jì)算速度較慢，但是通過算法的改進(jìn)，在保證一定計(jì)算精度的前提下，能夠滿足臨床對計(jì)算速度的要求。以點(diǎn)核為模型的卷積算法分為兩步[10]，首先由質(zhì)量衰減系數(shù)計(jì)算原射線穿過人體組織時(shí)單位質(zhì)量釋放的總能量 (total energy released per unit mass，TERMA)，然后根據(jù)事先計(jì)算得到的點(diǎn)核，經(jīng)過組織密度不均勻修正，與TERMA卷積，最終得到人體組織中的劑量分布。

對于具有n3個(gè)體素的三維數(shù)據(jù)，以點(diǎn)核為模型的標(biāo)準(zhǔn)卷積劑量計(jì)算算法的時(shí)間復(fù)雜度為O(n6)，如果再考慮組織密度的修正，計(jì)算復(fù)雜度將達(dá)到O(n7)，難以滿足快速計(jì)算的要求。Ahnesj?提出了collapsed cone卷積算法[1]，此算法的時(shí)間復(fù)雜度為O(mn3)，其中 m為能量沉積核的沉積方向數(shù)。collapsed cone卷積算法將散射劑量的積分劃分為具有一定立體角的同軸錐形筒串，位于筒串軸的散射體積單元釋放的能量沿筒串軸線性傳遞、衰減、沉積。

式(1)為點(diǎn)核 hm(r)的表達(dá)式，其中 am、bm、Am，Bm為某一能量沉積方向的沉積系數(shù)，可以通過對蒙特卡羅程序計(jì)算得到的能量沉積核進(jìn)行擬合得到。

對各個(gè)能量沉積方向，按照式(2)分別計(jì)算原射線劑量和散射線劑量并求和

式中，m、n分別為相對于射束入射方向的天頂角和方位角的離散方向索引，天頂角和方位角的角增量分別記為Δσθ和Δφ，r為劑量計(jì)算區(qū)域中的任意一點(diǎn)的位置向量，分別表示原射線和散射線劑量。

對于每條能量沉積射線上體素的原射線劑量和散射線劑量，可分別通過式(3)和式(4)迭代求出[1]，此兩式等號右側(cè)第一項(xiàng)為某一方向能量沉積線上前方體素ri-1對當(dāng)前計(jì)算體素ri的劑量貢獻(xiàn)，第二項(xiàng)為當(dāng)前體素對自身的劑量貢獻(xiàn)。

式中，Ti為當(dāng)前體素單位質(zhì)量釋放的總能量，Ωmn為某一能量沉積方向筒串的立體角，ηi為當(dāng)前體素相對水的電子密度，用于對非均勻介質(zhì)的密度修正。

在collapsed cone算法中用平行于劑量沉積方向的一系列射線，來替代每一個(gè)體素的三維卷積方向進(jìn)行加速計(jì)算。為了快速計(jì)算平行的能量沉積射束與病人體數(shù)據(jù)網(wǎng)格的交點(diǎn)以及相交線段的長度，Ahnesj?提出了種子平面的概念[3]，穿過種子平面的平行射束與病人網(wǎng)格的交點(diǎn)，可以簡化為求其中一條射線與網(wǎng)格的交點(diǎn)坐標(biāo)，然后與射線在種子平面上的入射點(diǎn)坐標(biāo)相加獲得。整個(gè)劑量計(jì)算有以下步驟。

步驟1：計(jì)算單位質(zhì)量釋放的總能量(TERMA)

步驟1-1：根據(jù)機(jī)床角、機(jī)架角、機(jī)頭角、射野的形狀計(jì)算待計(jì)算區(qū)域的矩形包圍盒;

步驟1-2：計(jì)算包圍盒內(nèi)各體素的 CT值、質(zhì)量衰減系數(shù)、質(zhì)量密度、電子密度;

步驟1-3：根據(jù)射野的形狀和計(jì)算分辨率，計(jì)算射束內(nèi)射線的離散數(shù)目n;

步驟1-4：計(jì)算每條射線在包圍盒內(nèi)沉積的TERMA。

步驟2：計(jì)算劑量

步驟2-1：根據(jù)設(shè)定的天頂角增量Δθ和方位角增量 Δφ，將方位角和天頂角離散為(θi，φj);

步驟2-2：對于每一個(gè)能量沉積方向(θi，φj)計(jì)算種子平面的大小、能量沉積射線的數(shù)量、過原點(diǎn)的能量沉積射線與包圍盒內(nèi)體數(shù)據(jù)網(wǎng)格的交點(diǎn)和相交線段長度，然后根據(jù)坐標(biāo)平移關(guān)系，計(jì)算每一條能量沉積射線與體數(shù)據(jù)網(wǎng)格的交點(diǎn)坐標(biāo)，進(jìn)而計(jì)算每一條能量沉積射線在包圍盒內(nèi)沉積的劑量，并累加。

2 軟件系統(tǒng)結(jié)構(gòu)及功能

為了驗(yàn)證上述劑量計(jì)算方法的準(zhǔn)確性，我們在VS 2005平臺上開發(fā)了一套光子束放射治療計(jì)劃系統(tǒng)，其軟件結(jié)構(gòu)如圖1所示。主要包括DICOM數(shù)據(jù)的讀取、圖像分割、三維組織重建、計(jì)劃設(shè)計(jì)、劑量計(jì)算、劑量評估、數(shù)據(jù)庫管理七個(gè)模塊，每個(gè)模塊又可以細(xì)分為若干子模塊。系統(tǒng)首先讀入DICOM序列，以位圖格式顯示在軟件界面上;通過手動分割或者智能剪刀(live wire)分割方法[11-13]提取靶區(qū)和關(guān)鍵組織的輪廓線后，使用輪廓線重建算法、marching cubes表面重建算法或體繪制算法進(jìn)行三維重建;根據(jù)勾畫出的靶區(qū)，設(shè)定加速器及射野參數(shù)等，進(jìn)行劑量計(jì)算;將計(jì)算得到劑量數(shù)據(jù)以等劑量面、等劑量線或者以偽彩色的方式顯示在病人斷層數(shù)據(jù)上，并計(jì)算靶區(qū)及敏感組織的劑量體積直方圖，以便進(jìn)行劑量評估。最后把輪廓線數(shù)據(jù)、計(jì)劃設(shè)計(jì)的參數(shù)、劑量計(jì)算的參數(shù)及結(jié)果等存放到指定的路徑下，并把文件路徑保存到數(shù)據(jù)庫中，軟件系統(tǒng)界面示例圖如圖2所示。

圖1 放射治療計(jì)劃系統(tǒng)軟件結(jié)構(gòu)Fig.1 Software structure of radiation treatment planning system

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

為了驗(yàn)證上述劑量計(jì)算方法的準(zhǔn)確性，采用蒙特卡羅程序DOSXYZnrc進(jìn)行比較，在相同能譜、射野大小、水體模等條件下進(jìn)行了模擬計(jì)算，比較射線能量分別為6 MV和15 MV條件下的百分深度劑量(PDD)，及20 cm×20 cm方野6 MV射束條件下不同深度的離軸比。源皮距SSD=100 cm，水體模體積為30 cm×30 cm×30 cm，射野大小為3 cm×3 cm、5 cm ×5 cm、10 cm ×10 cm、20 cm ×20 cm，DOSXYZnrc程序的歷史次數(shù)為108，collapsed cone算法中天頂角增量和方位角增量(Δθ、Δφ)取(7.5°、22.5°)、(3.75°、22.5°)、(3.75°、15°)等 3 種情況，用來分析離散角對計(jì)算精度的影響。劑量計(jì)算平均誤差見表1，可以看出：射束能量增加誤差減少;射野增大誤差增加;計(jì)算分辨率提高誤差降低，當(dāng) Δθ、Δφ 取 3.75°、15°時(shí)劑量計(jì)算誤差最小，圖 3給出了該情況下，本算法與蒙特卡羅程序DOSXYZnrc計(jì)算結(jié)果比較的示例圖，圖中“DOSXYZnrc”表示蒙特卡羅程序 DOSXYZnrc計(jì)算的結(jié)果，“CC”表示本 collapsed cone算法計(jì)算得到的結(jié)果。圖3(a)是射線能量為6 MV，照射野分別為 3 cm ×3 cm、5 cm ×5 cm、10 cm ×10 cm、20 cm ×20 cm條件下的百分深度劑量。圖3(b)是射線能量為15 MV，照射野分別為 3 cm×3 cm、5 cm×5 cm、10 cm×10 cm、20 cm×20 cm條件下的百分深度劑量。圖3(c)是射線能量為6 MV，照射野大小為 20 cm ×20 cm，深度為 1.5、10、30 cm 條件下的離軸比。從圖3中的(a)和(b)可以看出：1)隨著射野或能量增大，百分深度劑量增加;2)當(dāng)深度較小時(shí)，collapsed cone算法的計(jì)算結(jié)果與蒙特卡羅模擬的結(jié)果一致，從圖3(c)亦可看出并可得出此結(jié)論，但隨著深度的增加，由于射束隨深度增加而硬化，導(dǎo)致能量沉積核發(fā)生變化，因而誤差略有增加;3)相比于高能射線，低能射線硬化效應(yīng)相對明顯，因而深度較大時(shí)的誤差相對較大;4)在射野較大的情況下，由于每一能量沉積方向需要計(jì)算的距離增大，因而相對于小射野，射野大時(shí)誤差略大。

表1 深度劑量計(jì)算平均誤差(%)表Tab.1 Mean calculated dose error

圖3 與蒙特卡羅程序DOSXYZnrc的劑量比較示例圖。(a)6 MV射線下百分深度劑量;(b)15 MV射線下百分深度劑量;(c)6 MV射線20 cm×20 cm射野下離軸比Fig.3 Comparision of dose with DOSXYZnrc.(a)PDDs for 6 MV beam;(b)PDDs for 15 MV beam;(c)OARs for 6 MV beam with field size of 20 cm×20 cm

4 結(jié)論與討論

對collapsed cone劑量計(jì)算的步驟以及加速算法進(jìn)行了深入研究，在VS 2005平臺上開發(fā)了一套功能相對完善的光子束放射治療計(jì)劃系統(tǒng)，運(yùn)用蒙特卡羅方法和水體模對該劑量計(jì)算方法的準(zhǔn)確性進(jìn)行驗(yàn)證。與蒙特卡羅方法比較結(jié)果表明，本文研究開發(fā)的collapsed cone劑量計(jì)算方法準(zhǔn)確可靠，具有重要的臨床應(yīng)用價(jià)值。百分深度劑量與射線的能量有關(guān)，臨床應(yīng)根據(jù)病灶的具體情況選擇合適能量的射線進(jìn)行治療。與常用的筆形束劑量計(jì)算方法相比，collapsed cone劑量計(jì)算方法具有更高的計(jì)算精度，其劑量計(jì)算精度與劑量計(jì)算分辨率(即collapsed cone算法的離散參數(shù))及射束硬化等因素有關(guān)。隨著計(jì)算分辨率的提高，計(jì)算時(shí)間相應(yīng)增加，計(jì)算復(fù)雜度增加、速度降低;射束硬化會導(dǎo)致劑量計(jì)算誤差增加。對于射束硬化引起的計(jì)算誤差，通過對能譜進(jìn)行硬化校正，可以減小計(jì)算誤差，提高劑量計(jì)算的精度;由于劑量需要在三維空間進(jìn)行計(jì)算，且在放射治療計(jì)劃優(yōu)化的過程中需要反復(fù)進(jìn)行劑量計(jì)算，因此在確保劑量計(jì)算精度的前提下，劑量計(jì)算的速度也至關(guān)重要，GPU具有高效的并行計(jì)算能力且成本低，近幾年來在諸多領(lǐng)域得到眾多研究人員的關(guān)注，對于劑量計(jì)算同樣具有較高的實(shí)用價(jià)值。今后我們將研究基于GPU的劑量計(jì)算方法，以提高劑量計(jì)算的速度，同時(shí)研究劑量計(jì)算中的射束硬化校正方法，進(jìn)一步提高劑量計(jì)算的精度，并完善軟件系統(tǒng)功能。

[1]Ahnesj? A.Collapsed cone convolution of radiant energy for photon dose calculation in heterogeneous media[J].Med Phys，1989，16(4)：577-592.

[2]Hobant PW，Murray DC，Round WH.Photon beam convolutionusing polyenergetic energy deposition kernels[J].Phys Med Biol，1994，39(4)：669-685.

[3]Ahnesj? A.Cone discretization for the collapsed cone algorithm[C]//Proceedings of International Conference on the Use of Computers in Radiation Therapy.XII ICCR.Madison：Leavitt DD and Starkschall G Medical Physics Publishing，1997：114-116.

[4]Maspradakis M， Morrison RH， Richmond N， et al.Experimental verification of convolution/superposition photon dose calculations for radio therapy treatment planning[J].Phys Med Biol，2003，48：2873-2893.

[5]Lu Weiguo， Olivera1 GH， Chen Mingli， et al. Accurate convolution/superposition for multi-resolution dose calculation using cumulative tabulated kernels[J].Phys Med Biol，2005，50(4)：655-680.

[6]Zhou B，Hu XS，Chen DZ，et al.Hardware acceleration for 3-D radiation dose calculation[C]//Proceedings of 18 th IEEE International Conference on Application-specific Systems，Architectures and Processors.Montreal：IEEE，2007：290-296.

[7]Carlsson TA，Ahnesj? A.Optimization of the computational efficiency of a 3 D，collapsed cone dose calculation algorithm for brachytherapy[J].Med Phys，2008，35(4)：1611-1618.

[8]Hissoiny S， Ozell B. Fast convolution-superposition dose calculation on graphics hardware[J].Med Phys，2009，36(6)：1998-2005.

[9]陳朝斌，黃群英，吳宜燦，等.蒙特卡羅方法在放療計(jì)劃中的應(yīng)用[J].核技術(shù)，2006，29(1)：22-28.

[10]Ahnesj? A，Aspradakis MM.Dose calculations for external photon beams in radiotherapy[J].Phys Med Biol，1999，44(11)：99-155.

[11]Mortensen EN， Barrett WA. Interactive segmentation with intelligent scissors[J].Graphical Models and Image Processing，1998，60(5)：349-384.

[12]高新波，雷云，姬紅兵.一種改進(jìn)的 Live-Wire交互式圖像分割算法[J].系統(tǒng)工程與電子技術(shù)，2003，25(8)：915-918.

[13]楊金柱，趙姝穎，胡英，等.序列醫(yī)學(xué)圖像三維分割的一種方法[J].系統(tǒng)仿真學(xué)報(bào)，2005，17(12)：2896-2900.