专利名称：一种多叶准直器静态调强叶片序列算法的制作方法
技术领域：
本发明属于多叶准直器算法领域，是一种搜索多叶准直器静态调强叶片序列最 优路径的方法，能对强度分布图进行自主搜索最优路径并得出叶片运动的最佳序列。
背景技术：
多叶准直器用于控制加速器或者其他放射源的电子或光子射线的强度与形状， 其一般由20 - 120对叶片组成，每个叶片在等中心平面上的投影宽度约为0. 3cm至 1. 2cm不等。微型多叶准直器的每个叶片投影宽度可达1. 6mm。各组叶片相向排列安 装于托架上，叶片可以在托架上运动，其运动的方向与射线束轴线方向垂直，每个 叶片可独立驱动。准直器叶片的大小、片数和位置决定了照射野大小和形状近似的 精确度。计划系统输出多叶准直器各叶片的位置、运动速度等参数，通过调节各照 射野形状和照射时间，实现调强。输出叶片位置和速度等参数的策略统称多叶准直 器的叶片序列算法，由于其可以体现许多数学理论的概念如最优化方法、路径最优 规划、贪心算法等，因此有许多研究机构对其进行了理论研究。
多叶准直器的适形调强方式主要有多叶准直器静态调强、动态调强。静态调强 的原理是将射野要求的强度分布进行分级，利用多叶准直器形成的多个子野进行分 布照射。其特征是每个子野照射完毕后，照射切断，多叶准直器调整到另一个子野， 再继续照射，直到所有子野照射完毕，所有子野的流强相加，形成要求的强度分布。 多叶准直器静态调强的优点是控制简单，只需控制叶片的位置，不需要控制速度， 可调制成任意阶梯形的剂量分布。静态调强的研究起步较早，由于其便于剂量验证 和机械控制，因而成为调强的主流。
现有基于多叶准直器的控制系统，生成叶片序列的过程比较固定，即按照一定 理论的进行矩阵分解并将子野实施。如果计算出的子野个数过多，使子野数目过多 就舍去最小的子野。这类方法缺少智能性和自主性，对于最优路径的选择缺乏认知 能力，在实际应用中有一定的局限性。
发明内容
本发明的目的在于针对现有的叶片序列输出的基础子野分割算法进行改进，以便使子野的个数(Number of Segment, NS )尽可能的少，所需要的子野跳数(Total Number of Monitor Unit, TNMU)尽可能少，而提出一种多叶准直器静态调强叶片序 列算法。在各种几何目标约束情况下计算最优双向匹配和最小费用流，从而得到效 率更高的多叶准直器叶片动作序列。 本发明的技术方案如下
一种多叶准直器静态调强叶片序列算法，其特征在于包括以下步骤 (1)通过计划系统计算出的射束剖面强度分布优化结果将影响到剂量的分布， 但一个好的剂量分布对应的初始强度分布并不一定是最可取的，因为剂量适形程度 越高，要求的强度分布就越复杂。对于复杂的分布情形，初始射束剖面强度将变得 非常复杂，子野化后必将产生很多子野，甚至很多的小野。如此以来，将会大大延 长照射的时间，同时随着子野^:目的增加，特别是4艮多小子野的存在，射野利用率 将降低，同时会增加射野漏射(Leakage )、机头散射(Scat ter )和MLC的凸凹槽效 应(Tongue and Groove Effects )等因素带来的不确定性。假定从计划系统获得的 m行n列二维强度矩阵A为已知，通过Savitzky-Golay滤波、矩形域的最小二乘曲面拟 合算法、最佳一致逼近的里米兹方法等将A在保留基本特征的情况下，得到A'。如图
运用滤波算法到强度通量图A上会得到不运用滤波算法相近的剂量分布A,可以 使需要使用的TNMU和NS大量减少。RMS ( root-mean-squared )偏离是判断标准。滤 波算法的运用也有2阶，4阶等不同阶数的区别，图4的b、 c。
(2 )并已通过辐照区域离散化使矩阵A成为一个m行n列的非负强度矩阵A'，矩 阵A'的元素A' ( i， j)代表辐照区域的辐照强度值，0代表没有强度，即是没有辐照， 辐照强度值由小至大，则是表示辐照强度由小至大，即是需要照射单位的少与多； 一般说来，元素A' (i， j)必须满足如下的条件才能被大照射野多叶准直器所实现，
存在正整数A _1，"' +1,其具有如下特性
《)
6其中(+1各自代表了左右叶片的位置，所以强度矩阵A'可以分解为二维矩 阵Si的和
其中Si是MLC可以形成的每一个子野。S(i, j )是Si二维数组中第i行第j列的 元素值。
以子野的个数NS、子野跳数T画U最小化为目标，也就是以SUi值与k值各自都达 到最小值；采用Xia and Verhey算法、Galvin—Chen and Smith算法、Bortfeld算法、 Engel算法、Kalinowski算法、Kamath算法、Boland， Hamacher and Lenzen算法、 Baatar and Hamacher算法、Langer, Thai and Papiez算法、Dai and Zhu算法等 计算结果比较T腿U与NS值的大小后，选择能获得最小TN而与NS值的算法，并按此算 法将A分解而成的二维矩阵一系列的子野S,。
(3) S,对应与MLC的叶片序列关系，如图3。根据大照射野的硬件条件不同，将 射野划分成为多个射束元(beamlet),如下面的某种MLC的技术指标
等中心射野约81 x 270mm;
叶片对数27对(3mm，等中心处投影)；
最小辐射野尺寸3mm x 3mm;
最大过等中心距离90mm。
根据技术指标将上述MLC的调强区域转换为27 x 60的二维数组。将数学关系与实 际结合到了 一起。
所述的滤波算法，其特征在于Savitzky-Golay滤波、矩形域的最小二乘曲面拟 合算法、最佳一致逼近的里米兹方法。
所述的子野分割算法比较，其特征在于采用Xia and Verhey算法、Galvin_Chen and Smith算法、Bortfeld算法、Engel算法、Kalinowski算法、Kamath算法、 Boland, Hamacher and Lenzen算法、Baatar and Hamacher算法、Langer， Thai andPapiez算法、Dai and Zhu算法。在上述算法范围内达到子野数目或机器跳数 最小，减少加速器和多叶准直器机械损耗。
所述的滤波算法在不改变强度通量图主要特性的基础上减少子野数目，从而减 少了叶片机械损耗，并缩短了照射时间。保留的主要特性的评价标准，其特征在于 RMS ( root—mean-squared )偏离。根据权利要求l所述的滤波算法，其特征在于可以根据用户需要按照子野重要性逐步 减少相对不重要的子野。
本发明在数学上最优性的证明可以通过文献(Boyer A L and Strait J P. Delivery of intensity modulated treatments with dynamic multileaf collimators. In Proc. XII Int. Conf. on use of Computers in Radaiation Therapy. Madison, WI: Medical Physics Publishing, 1997. 13—15.)获4寻。
通过Xia and Verhey ，Galvin， Siochi等算法比较输出得2D强度之后，可以 利用图论知识将该2D矩阵转换为一个带权有向图，G= (V， E, W)以及左右两个指定的 顶点v和w (代表左右叶片)，问题转换为要寻找从v到w的最短路径问题。所谓 最短路径问题，就是给定一个赋权有向图。即给了 一个有向图D= (V, E),对每一个弧 a=(vi, vj)，相应地有权重w(a)=wij。又给定D中的两个顶点vs， vt。设P是D中从 vs到vt的一条路径，定义路径P的权是P中所有弧的权之和，记为w(P)。最短路径 就是要在所有从vs到vt的路中，求一条权重最小的路径，即求一条从vs到vt的路 径PO,使w(PO)-Min(w(PO))。这个算法将按照顶点离v增加的次序找从v到另外一 个顶点的最短路。当所计算的顶点到达w时算法停止。寻找最小路径的算法较多， 采用Floyd算法、动态规划算法、贪心算法。以Floyd算法为例
Floyd算法的基本思想是从vs到vt的最短路径是以下各种可能路径中的长 度最小的那条。
若〈vs, vt〉存在，则存在路径{vs, vt}。(路径中不含有其它顶点)
若〈vs, vl>, <vl, vP存在，则存在if各径{vs, vl, vt)。(路径中所含顶点序号不大于
1)
若〈vs，…，v2〉，〈v2,…vt〉存在，则存在一条最短;洛径(vs，…v2…vj)。(路径中所 含顶点序号不大于2)
依次类推，则vs到vt的最短路径应是上述这些路径中路径长度最小者。
关键思想编程，用伪代码表示
从F1 oyd算法的基本思想可以得出其算法的步骤如下
St印l:初始化距离矩阵D (0) = (di j (0))、序号矩阵S (0) = (si j (0))和执行次数变 量k的值。其中距离矩阵D (0)和序号矩阵S (0)分别定义为
8若vi邻接到vj,则dij (0)=wij;否则dij (0) = ~。
si j (O)-j (i, j=l， 2…n),即元素si j (0)的值等于它所在的列数。图中顶点的个 数为n。 k=l
St印2:当k《n时，D(k)中第k行与第k列元素保持与D(k-1)的相应元素相同。 其他元素按下式进行计算
dij (k) =Min(dij (k—1), dik (k—1)+dkj (k—1)) 序号矩阵S (0)的元素取值规则为
若 dij(k)=dij(k-l),则 sij(k)=sij(k-l);若 dij(k)<dij(k-l),则 sij(k)=sik(k)。
经过计算，若dii (k)<0,则结束整个过程的计算。说明图中存在一条含有顶点vi 的负回if各，由序号矩阵中的sij (k)可以找出此回if各，否
则k-k+l，再执行St印2。
St印3:当1[=11时，终止整个过程的计算。若^"1{)=+00,则说明D中不存在从vs 到vt的最短路径；否则di j (k)的值就是vi到vj的最短距离。其相应的路径可由序 号矩阵中的sik(k)找出。(i,j=l，2...n)
从以上介绍的具体算法中可以看出Floyd算法是从邻接矩阵出发，而且最短路 径可以从序号矩阵中找出来。最短路径的权值从距离矩阵中得出。
得到的Ol数组，根据要求再采用Savitzky-Golay算法进行子野过滤(可选)， 进一步减少子野数目，最后按多叶准直器的几何参数转换为叶片序列文件控制多叶 准直器。
本发明的有益成果为以便使子野的个数(Number of Segment, NS)尽可能的 少，所需要的子野跳数(Total Number of Monitor Unit, TNMU)尽可能少，大大减 少了硬件损耗和尽可能缩短了照射的时间。
图l建立多叶准直器叶片序列输出过程示意图。 图2为算法流程图。
图3为2D数组与对应MLC形状，中间有洞的多边形。图4滤波算法使用效果图；(a)是原数组，(b) (c)是二阶、四阶滤波效果图。
在实际应用的具体实施实例 实例
有如下强度矩阵A:
爿=
2111113、
75625441
14566657
67627662
27724557
55544563
54456777
,11213311,
步骤l: Xia and Verhey， Galvin， Siochi, Bortfeld算法等分解结果如下:
Xia and VerheyGalvinSiochiBortfeld
子野数目(NS)78712
表l子野数目的计算结果
按照表l的计算结果，采用Xia and Verhey算法取初始强度
u = INT(log2(I』-1 = 4。
00000000 1110 1110 0 1111111
1110 1110 0 110 111 1
1111 1110 11111111 00000000
将此结果按照1的连续性分为MIX可以执行的形状。对应如图3。按照F1 oyd 算法对2D矩阵进行计算，寻找最小路径。
10编程方法(伪代码)如下 For k=l to n For i=l to n For j=l to n If i=k Then
For c=l to n
dkc (k) =dkc (k-1) #D (k)中第k行与与D(k-l)的第k行的元 素对应相同
dck(k)=dck(k—l)#D(k)中第k列与与D(k-l)的第k列的元 素对应相同 Next
Else
If j#k Then dij (k)=Min(dij (k—1)， dik(k-l)+dkj (k-1) If dij (k) = (dij (k-l)Then sij(k)=sij(k-1)
Else
sij (k)=sik(k)
End End
Next
Next Next
该算法的时间复杂性为0(n3+n2)。计算结果为3,比预期的Xia算法子野 分割数目(2D强度下)要小。因为子野数目已经比较小，所以不需要进行滤波。 如果滤波程序就可以将子野数目按权重能力过滤掉小子野。
权利要求
1、一种多叶准直器静态调强叶片序列算法，其特征在于包括以下步骤(1)通过计划系统计算出的射束剖面强度分布优化结果将影响到剂量的分布，但一个好的剂量分布对应的初始强度分布并不一定是最可取的，因为剂量适形程度越高，要求的强度分布就越复杂；对于复杂的分布情形，初始射束剖面强度将变得非常复杂，子野化后必将产生很多子野，甚至很多的小野；如此以来，将会大大延长照射的时间，同时随着子野数目的增加，特别是很多小子野的存在，射野利用率将降低，同时会增加射野漏射(Leakage)、机头散射(Scatter)和MLC的凸凹槽效应(Tongue and Groove Effects)等因素带来的不确定性；假定从计划系统获得的m行n列二维强度矩阵A为已知，通过Savitzky-Golay滤波、矩形域的最小二乘曲面拟合算法、最佳一致逼近的里米兹方法等将A在保留基本特征的情况下，得到A′；运用滤波算法到强度通量图A上会得到不运用滤波算法相近的剂量分布A，可以使需要使用的TNMU和NS大量减少；RMS(root-mean-squared)偏离是判断标准；滤波算法的运用也有2阶，4阶等不同阶数的区别；(2)并已通过辐照区域离散化使矩阵A成为一个m行n列的非负强度矩阵A′，矩阵A′的元素A′(i，j)代表辐照区域的辐照强度值，0代表没有强度，即是没有辐照，辐照强度值由小至大，则是表示辐照强度由小至大，即是需要照射单位的少与多；一般说来，元素A′(i，j)必须满足如下的条件才能被大照射野多叶准直器所实现，存在正整数li-1，ri+1，其具有如下特性1≤li≤ri+1≤n+1 (i∈[m])，其中li-1，ri+1各自代表了左右叶片的位置，所以强度矩阵A′可以分解为二维矩阵Si的和其中Si是MLC可以形成的每一个子野。S(i，j)是Si二维数组中第i行第j列的元素值；以子野的个数NS、子野跳数TNMU最小化为目标，也就是以∑Ui值与k值各自都达到最小值；采用Xia and Verhey算法、Galvin_Chen and Smith算法、Bortfeld算法、Enge1算法、Kalinowski算法、Kamath算法、Boland，Hamacher and Lenzen算法、Baatar and Hamacher算法、Langer，Thai and Papiez算法、Dai and Zhu算法等计算结果比较TNMU与NS值的大小后，选择能获得最小TNMU与NS值的算法，并按此算法将A分解而成的二维矩阵一系列的子野Si；(3)Si对应与MLC的叶片序列关系根据大照射野的硬件条件不同，将射野划分成为多个射束元(beamlet)，如下面的某种MLC的技术指标等中心射野约81×270mm；叶片对数27对(3mm，等中心处投影)；最小辐射野尺寸3mm×3mm；最大过等中心距离90mm；根据技术指标将上述MLC的调强区域转换为27×60的二维数组，将数学关系与实际结合到了一起。
2、 根据权利要求l所述的多叶准直器静态调强叶片序列算法，其特征在于所述的滤 波算法，是Savitzky-Golay滤波、矩形域的最小二乘曲面拟合算法、最佳一致 逼近的里米兹方法。
3、 根据权利要求l所述的多叶准直器静态调强叶片序列算法，其特征在于所述的子 野分割算法比较，采用Xia and Verhey算法、Galvin—Chen and Smith算法、 Bortfeld算法、Engel算法、Kalinowski算法、Kamath算法、Boland， Hamacher and Lenzen算法、Baatar and Hamacher算法、Langer， Thai and Papiez算 法、Dai and Zhu算法，在上述算法范围内达到子野^:目或机器跳数最小，减 少加速器和多叶准直器机械损耗。
4、 根据权利要求l所述的多叶准直器静态调强叶片序列算法，其特征在于所述的滤 波算法在不改变强度通量图主要特性的基础上减少子野数目，从而减少了叶片 机械损耗，并缩短了照射时间；保留的主要特性的评价标准，其特征在于RMS(root-mean-squared )偏离。
5、 根据权利要求l所述的多叶准直器静态调强叶片序列算法，其特征在于所述的滤波算法，根据用户需要按照子野重要性逐步减少相对不重要的子野。
全文摘要
本发明公开了一种多叶准直器静态调强叶片序列算法，综合比较现有算法得出2d强度分布，再利用图论分解2d强度分布并加入子野过滤的叶片序列算法的整体解决方案。在获得靶区强度矩阵后，对其分群并利用多种算法(Galvin，Xia andverhey，Siochi，Bortfeld等)进行子野分割的结果比较总加速器跳数(TNMU)和子野分割数目(NS)得出最优值，后选定子野分割方案，得到2d强度分布，并以图论中最小路径法计算叶片序列，并可根据需要按MU权重过滤子野，以期达到符合需求的优化叶片序列算法。
文档编号G06F17/10GK101499326SQ200910116119
公开日2009年8月5日 申请日期2009年1月22日 优先权日2009年1月22日
发明者兰海洋, 吴宜灿, 耀 孟, 佳 景, 曹瑞芬, 贵 李, 曦 裴, 郑华庆 申请人:中国科学技术大学