一種改進(jìn)的多葉準(zhǔn)直器射野分割算法

張鵬程 朱 健，2 李元金，3 舒華忠

(1東南大學(xué)影像科學(xué)與技術(shù)實(shí)驗(yàn)室，南京210096)

(2山東省腫瘤醫(yī)院，濟(jì)南250117)

(3滁州學(xué)院計(jì)算機(jī)科學(xué)與技術(shù)系，滁州239000)

在調(diào)強(qiáng)放療(intensity-modulated radiotherapy，IMRT)計(jì)劃的實(shí)施過(guò)程中，可以采用多葉光柵靜態(tài)調(diào)強(qiáng)和多葉光柵動(dòng)態(tài)調(diào)強(qiáng)2種模式.前者按照射野要求的強(qiáng)度分布進(jìn)行分割，利用多葉光柵(multileaf collimator，MLC)形成的多個(gè)子野進(jìn)行分布照射，其特征是在每個(gè)子野照射完畢后切斷照射，并將MLC調(diào)到另一個(gè)子野位置，繼續(xù)照射，直到所有子野照射完畢，將所有子野的流強(qiáng)相加，即可形成要求的強(qiáng)度分布［1-2］.后者則利用多葉光柵相對(duì)應(yīng)的一對(duì)葉片的相對(duì)運(yùn)動(dòng)，實(shí)現(xiàn)對(duì)射野強(qiáng)度的調(diào)節(jié)，其特征是葉片運(yùn)動(dòng)過(guò)程中，射線一直處于“Beam on”的位置，劑量率不變［2］.與動(dòng)態(tài)調(diào)強(qiáng)模式相比，靜態(tài)調(diào)強(qiáng)模式具有硬件實(shí)施容易控制、可靠性高等優(yōu)點(diǎn)［3-4］.前期研究在對(duì)比了2種調(diào)強(qiáng)模式后認(rèn)為，靜態(tài)調(diào)強(qiáng)模式在降低正常器官受照射劑量方面優(yōu)于動(dòng)態(tài)調(diào)強(qiáng)模式［5］.因此，本文主要研究靜態(tài)調(diào)強(qiáng)模式中的多葉準(zhǔn)直器射野分割算法.

在靜態(tài)模式中，過(guò)多的總子野個(gè)數(shù) (total number of segments)將延長(zhǎng)治療時(shí)間，過(guò)多的總機(jī)器跳數(shù)(total number of monitor units)將增加葉片間的射線泄露.因此，如何在保證投射劑量的前提下減少總子野個(gè)數(shù)和總機(jī)器跳數(shù)是優(yōu)化子野分割的兩大目標(biāo)［6］.與此同時(shí)，為了減少葉片間的凸凹槽效應(yīng)、提高射線的利用率，在滿足上述目標(biāo)的基礎(chǔ)上，子野的開口面積也應(yīng)設(shè)計(jì)得盡可能大.

本文主要考慮的是相鄰葉片允許交疊情況下子野個(gè)數(shù)最少的射野分割算法.針對(duì)這類射野分割問(wèn)題，目前的解決方法主要是自適應(yīng)算法［7-9］和規(guī)劃類算法［10-11］.但對(duì)于大規(guī)模的矩陣分割，這些方法會(huì)影響計(jì)算速度.在Luan算法中，首先得到不同的多葉準(zhǔn)直器開口形狀和機(jī)器跳數(shù)對(duì)，然后將每行中不同多葉準(zhǔn)直器的開口形狀整合成相應(yīng)的子野，在保證每行MLC開口最大的情況下減少子野個(gè)數(shù)［12］.Engel算法是一種自適應(yīng)算法，在控制總機(jī)器跳數(shù)的同時(shí)，所需的子野個(gè)數(shù)幾乎最少［7］.本文結(jié)合Engel算法所需較少總子野個(gè)數(shù)的優(yōu)點(diǎn)，對(duì)Luan算法進(jìn)行改進(jìn).改進(jìn)后的算法在保證每行中準(zhǔn)直器開口最大的同時(shí)，減少了所需子野個(gè)數(shù)和總機(jī)器跳數(shù)，且所需計(jì)算量較Engel算法更少.

1 算法描述

將m×n的二維矩陣A用于描述來(lái)自于放療逆向計(jì)劃的射野強(qiáng)度圖.A可被分割為包含不同機(jī)器跳數(shù)的子野，即

式中，N為總子野個(gè)數(shù);uk＞0為第k個(gè)子野的強(qiáng)度值;Sk為第k個(gè)子野的形狀，其元素由0，1構(gòu)成，分別表示該位置已被、未被MLC遮擋，且每一行中的所有1必須連續(xù).Sk中元素的值由該子野中每行MLC左、右葉片的位置li，ri決定，即

例如，矩陣A被分割為4個(gè)子野，總機(jī)器跳數(shù)為5+2+1+1=8，則A可以表述為

本文只考慮相鄰葉片允許交疊情況下的射野分割算法.為了描述方便，在矩陣A中加入2列零元素，即

并定義其m×(n+1)維橫向差分矩陣D，其中元素表示為

矩陣A中第i行機(jī)器跳數(shù)復(fù)雜度定義為

矩陣A的機(jī)器跳數(shù)復(fù)雜度為所有行中機(jī)器跳數(shù)復(fù)雜度的最大值，即

第i行的復(fù)雜度間隔為

1.1 Luan 算法

在保證每行MLC開口最大的情況下，可利用Luan算法來(lái)減少子野個(gè)數(shù).算法的具體步驟如下.

①將矩陣A按行分割成m個(gè)行向量Xj，其中j∈［1，m］.

②將每一個(gè)行向量Xj分割成多個(gè)MLC開口形狀和機(jī)器跳數(shù)對(duì)(Sj，i，vj，i)，并將機(jī)器跳數(shù)值 vj，i以∑ai2i的形式分割成多個(gè)MLC形狀和機(jī)器跳數(shù)對(duì)(Sj，i，20)，(Sj，i，21)，(Sj，i，22)，…，(Sj，i，2m).例如，可將(Sj，i，5)分割為(Sj，i，21)，(Sj，i，22).然后，將分割后的MLC開口形狀和機(jī)器跳數(shù)對(duì)(Sj，i，v)按機(jī)器跳數(shù)值重新排序?yàn)?Sj，i，vj，i).

③將每行不同的MLC形狀整合成相應(yīng)的子野.算法描述如下:

Biedl等［13-14］將該算法擴(kuò)展為按 ∑ai3i和∑ai4i的指數(shù)和形式分割，以進(jìn)一步減少子野個(gè)數(shù).

1.2 Engel算法

在保證總機(jī)器跳數(shù)最小的情況下，利用Engle算法可以得到近似最優(yōu)的子野個(gè)數(shù).當(dāng)矩陣A的復(fù)雜度為0時(shí)，射野分割結(jié)束.該算法中，為了減少子野個(gè)數(shù)，每次分割時(shí)盡可能取大的機(jī)器跳數(shù)值uk，以盡可能多地減少矩陣A的復(fù)雜度，即

式(6)中，uk取值為所有行中機(jī)器跳數(shù)值u'k的最小值，u'k為每一行任意準(zhǔn)直器葉片處最大可取的機(jī)器跳數(shù)值.對(duì)于機(jī)器跳數(shù)復(fù)雜度最大行，任意準(zhǔn)直器葉片位置處u'k的取值須滿足以下條件:

式中，i*表示機(jī)器跳數(shù)復(fù)雜度最大行.

對(duì)于其他行，u'k的取值要滿足該次射野分割后該行機(jī)器跳數(shù)復(fù)雜度不大于機(jī)器跳數(shù)復(fù)雜度最大行的復(fù)雜度，即

每次射野分割后，第i行的機(jī)器跳數(shù)復(fù)雜度為

式(9)可以變換為

根據(jù)每行中左、右葉片處的強(qiáng)度差分值，可以分以下2種情況來(lái)確定第i行中的最大機(jī)器跳數(shù)取值 ki(r，l):

1)當(dāng)該行準(zhǔn)直器開口為空時(shí)，則

2)當(dāng)該行準(zhǔn)直器開口非空時(shí)，則

① 若 ki(r，l)≤min{di，li，－ di，ri+1}，則不等式(10)恒成立;

② 若 min{di，li，－ di，ri+1}≤ki(r，l)≤max{di，li，－ di，ri+1}，則

③ 若 ki(r，l)≥max{di，li，－ di，ri+1}，則

為了保證分割后矩陣A中所有元素非負(fù)，每行中的最大機(jī)器跳數(shù)值ki(r，l)應(yīng)小于等于矩陣A中所有元素的最小值.第i行中準(zhǔn)直器開口內(nèi)矩陣A中所有元素的最小值為

當(dāng)前行中最大允許的機(jī)器跳數(shù)取值為

矩陣A中最大允許的機(jī)器跳數(shù)取值為

已知矩陣A中最大允許的機(jī)器跳數(shù)取值uk后，可根據(jù)子野開口面積最大、總機(jī)器跳數(shù)最小、子野開口面和對(duì)應(yīng)機(jī)器跳數(shù)的乘積最大等不同準(zhǔn)則來(lái)確定對(duì)應(yīng)MLC開口形狀Sk.Engle算法在減少子野個(gè)數(shù)的同時(shí)，增加了部分行的復(fù)雜度，而且在這些行中生成的MLC開口都較小.

1.3 改進(jìn)算法

Luan算法和Biedl等［13-14］提出的改進(jìn)算法都存在一個(gè)共同的問(wèn)題，即對(duì)某些數(shù)值的處理反而導(dǎo)致子野個(gè)數(shù)增加.例如，當(dāng)前子野的機(jī)器跳數(shù)為7時(shí)，如果按∑ai2i形式進(jìn)行分次照射，該子野應(yīng)該繼續(xù)分割為機(jī)器跳數(shù)分別為4，2，1的3個(gè)子野，與直接按機(jī)器跳數(shù)7進(jìn)行照射相比，增加了2個(gè)子野.為此，本文結(jié)合Engel算法的優(yōu)點(diǎn)，對(duì)Luan算法提出進(jìn)一步改進(jìn).改進(jìn)算法的具體步驟如下.

①將矩陣A按行分割成m個(gè)行向量Xj，其中j∈［1，m］.

②將每一個(gè)行向量Xj分割成若干個(gè)MLC開口形狀和機(jī)器跳數(shù)對(duì)(Sj，i，vj，i).

③將每行不同的MLC開口形狀整合成相應(yīng)的子野.算法描述如下:

與Engel算法相比，改進(jìn)算法在減少子野個(gè)數(shù)的同時(shí)，盡可能保證了每行的MLC開口面積最大.與Luan算法相比，改進(jìn)算法不再受指數(shù)函數(shù)的限制，減少了部分不必要增加的子野，降低了子野數(shù)，而且沒(méi)有增加機(jī)器總跳數(shù).

2 實(shí)驗(yàn)結(jié)果

利用本文算法對(duì)1 000個(gè)15×15不同最大射野強(qiáng)度的隨機(jī)矩陣進(jìn)行分割，并將結(jié)果與經(jīng)典的Bortfeld 算法［6］、Xia 算法［4］、Engel算法以及 Luan算法得到的結(jié)果進(jìn)行比較(見(jiàn)表1).表1中，T表示1 000個(gè)隨機(jī)矩陣的平均總子野個(gè)數(shù);M表示平均總機(jī)器跳數(shù);L表示隨機(jī)矩陣的最大值.由表可知，與Luan算法相比，本文算法在減少子野個(gè)數(shù)的同時(shí)，降低了總機(jī)器跳數(shù);與Engel算法相比，本文算法沒(méi)有增加總機(jī)器跳數(shù)，僅平均增加了1個(gè)子野.究其原因，Engel算法為了保證在每次分割中盡可能多地減少?gòu)?qiáng)度矩陣的復(fù)雜度，增加了部分行的子野個(gè)數(shù)，而且增加的這些子野的MLC開口都比較小;而本文算法在進(jìn)行分割的時(shí)候，為了確保每行中的MLC開口盡量大，減小了對(duì)應(yīng)的機(jī)器跳數(shù)值，因而增加了總子野個(gè)數(shù).

射野分割算法的運(yùn)行時(shí)間與強(qiáng)度矩陣的大小以及每個(gè)像素的強(qiáng)度值有關(guān).矩陣和強(qiáng)度值越大，所需分割時(shí)間越長(zhǎng).根據(jù)實(shí)際準(zhǔn)直葉片對(duì)數(shù)(如美國(guó)Varian公司用于常規(guī)放療的MLC葉片對(duì)數(shù)有26，40或60對(duì))，對(duì)最大強(qiáng)度值為26、大小為26×26，40×40和60×60的隨機(jī)矩陣的分割時(shí)間進(jìn)行比較.在內(nèi)存為2 GB的雙核英特爾酷睿臺(tái)式機(jī)上運(yùn)行Engel算法、Luan算法及本文算法，結(jié)果見(jiàn)表2.由表可知，與Engel算法相比，本文算法和Luan算法所需計(jì)算時(shí)間較短.而且隨著矩陣的變大，本文算法在縮短計(jì)算時(shí)間上的優(yōu)勢(shì)更加明顯.

表1 本文算法與經(jīng)典算法的計(jì)算結(jié)果比較

表2 本文算法與經(jīng)典算法的計(jì)算時(shí)間比較 s

3 結(jié)語(yǔ)

本文通過(guò)結(jié)合Engel算法的優(yōu)點(diǎn)，對(duì)Luan算法進(jìn)行了改進(jìn).改進(jìn)后的算法不再受指數(shù)函數(shù)的限制，進(jìn)一步減少了射野分割所需的子野個(gè)數(shù).與傳統(tǒng)算法相比，本文算法能在減少子野個(gè)數(shù)和機(jī)器總跳數(shù)的同時(shí)，保證每行中準(zhǔn)直器的開口最大.與Engel算法相比，本文算法所需計(jì)算量更小.

考慮到臨床實(shí)際需求，IMRT計(jì)劃產(chǎn)生的強(qiáng)度矩陣遠(yuǎn)大于15×15，由此帶來(lái)的時(shí)間和劑量學(xué)上的優(yōu)勢(shì)將更為顯著.影像引導(dǎo)放射治療的廣泛應(yīng)用需要在患者不離開治療床的情況下在線快速修改及重新優(yōu)化放療計(jì)劃，這也為本文算法提供了廣闊的發(fā)展空間和良好的應(yīng)用前景.

References)

［1］Boyer A L，Yu C X.Intensity-modulated radiation therapy with dynamic multileaf collimators［J］.Semin Radiat Oncol，1999，9(1):48-59.

［2］胡逸民.腫瘤放射物理學(xué)［M］.北京:原子能出版社，1999:550.

［3］Dai J，Zhu Y.Minimizing the number of segments in a delivery sequence for intensity-modulated radiation therapy with multileaf collimator［J］.Med Phys，2001，28(10):2113-2120.

［4］Xia P，Verhey L J.Multileaf collimator leaf sequencing algorithm for intensity modulated beams with multiple static segments［J］.Med Phys，1998，25(8):1424-1434.

［5］朱健，劉敏，李建彬，等.多葉光柵靜態(tài)調(diào)強(qiáng)技術(shù)對(duì)鼻咽癌放療中危及器官受量的作用［J］.中華放射醫(yī)學(xué)與防護(hù)雜志，2009，26(2):171-173.Zhu Jian，Liu Min，Li Jianbin，et al.The role of the organ at risk in radiotherapy for nasopharyngeal carcinoma using multileaf in collimator static intensity modulated［J］.Chin J Radio Med Prot，2009，26(2):171-173.(in Chinese)

［6］Kalinowski T.Realization of intensity modulated radiation fields using multileaf collimators［C］//Lecture Notes in Computer Science.Berlin:Springer-Verlag，2006，4123:1010-1055.

［7］Engel K.A new algorithm for optimal multileaf collimator field segmentation［J］.Discrete Appl Math，2005，152(1/2/3):35-51.

［8］Kalinowski T.A duality based algorithm for multileaf collimator field segmentation with interleaf collision constraint［J］.Discrete Appl Math，2005，152(1/2/3):52-88.

［9］Siochi R A C.Variable depth recursion algorithm for leaf sequencing［J］.Med Phys，2007，34(2):664-672.

［10］Cambazard H，O'Mahony E，O'Sullivan B.A shortest path-based approach to the multileaf collimator sequencing problem［J］.Discrete Appl Math，2012，160(1/2):81-99.

［11］Taskin Z，Smith J，Romeijn H，et al.Optimal multileaf collimator leaf sequencing in IMRT treatment planning［J］.Operations Research，2010，58(3):674-690.

［12］Luan S，Saia J，Young M.Approximation algorithms for minimizing segments in radiation therapy［J］.Information Processing Letters，2007，101(6):239-244.

［13］Biedl T，Durocher S，Hoos H.et al.A note on improving the performance of approximation algorithms for radiation therapy［J］.Information Processing Letters，2011，111(7):326-333.

［14］Biedl T，Durocher S，Engelbeen C，et al.Faster optimal algorithms for segment minimization with small maximal value［C］//Proceedings of the 12th International Conference on Algorithms and Data Structures.New York，USA，2011:86-97.