伽玛刀治疗计划.docx

1、伽玛刀治疗计划大连理工大学第十届数学建模竞赛参赛作品选题：A.伽玛刀治疗计划组别：第61组队员：计算机系 0102班 刘力协计算机系 0101班 杨 威计算机系 0101班 冯 维2003年4月20日伽玛刀治疗计划摘要我们在对剂量分布问题的充分研究的基础上，提出利用中轴变换抽取出病灶的骨架（3D Skeleton），迅速掌握了病灶的特征，然后利用双重标准进行填充的计算机自动生成治疗方案的模型。我们利用二次积分和微分的概念求出了靶点附近剂量的分布函数，在此基础上讨论了有效剂量的选择对周围正常脑细胞的影响及不同准直器对周围脑细胞的影响，进一步讨论了两个靶点组合进行治疗时，互相影响所引起的等剂量曲线

2、的变化。通过对上面三个问题的分析，我们提出了用低剂量相互叠加，有效等剂量曲线和用来进行病灶覆盖的等剂量曲线（即我们在进行填充时的参考形状）的双重标准，我们认为70%等剂量曲线是有效剂量，而50%的剂量曲线是用来进行病灶覆盖的等剂量曲线。这样就能够在放置靶点的同时，将靶点之间的相互影响的问题也同时考虑进去，从而避免了计算量很大的迭代，能够一步到位，而且使的我们的模型很有弹性。这样就避免了纯粹用几何来解决时遇到的几何和剂量分布差别较大的缺点。在具体建立模型过程中，我们在一个50X50X50（625，000个像素）的正方形网格中讨论问题，这是能够满足一般情况下的技术精度的。因为没有现成的病灶，我们利

3、用随机和分形相结合方法生成了一个简陋的病灶（图五）。在图形处理方面，我们利用现在在图形保存和压缩方面比较常用的骨架抽取的方法，具体通过利用“带符号的距离地图”（8SSED, Sequential Euclidean Distance mapping algorithms6812）的思想来实现，通过记录病灶内部的像素到边界像素的最近距离，找到图形的中轴。抽取出病灶的骨架后，我们利用骨架的特征（最大内球球心的集合）来进行靶点对病灶进行从“角”到中心类似贪心算法的覆盖，在形成治疗方案后，再对其进行评价，利用我们模型的弹性化设计再进行一些局部优化，将得到一个满意的结果。我们的模型有很强的可操作性和弹性

4、，并且算法复杂度低，稳定性好，完全实现计算机自动处理。一 问题重述： 伽玛刀是世界上最先进的治疗仪器之一，主要通过伽玛刀组发射的高剂量的射线来达到摧毁肿瘤的目的。伽玛刀组有固定在头罩上的201个钴-60单元来组成，单一一条射线是不足以杀死肿瘤细胞的，而201束射线同时在同一中心交汇，形成一个大致球体的剂量分布（称为靶点），并达到有效的剂量水平，杀死肿瘤细胞。可以通过四个可互换的外部瞄准仪头罩（有4，8，14，18mm四种射线束通道）来照射出不同的剂量分布尺寸。如果目标的体积较大，可以通过多个靶点的组合来覆盖整个目标，在实际中，大多需要一到十五个靶点。现在治疗方案大多是凭医师的经验和一个迭代算法

5、，结果的满意程度在很大程度上取决于医师的经验。我们的模型就是要找到一种有效的方法，使计算机针对肿瘤的形状和周围健康组织的情况来自动生成一套治疗方案。二 假设条件：1 由核磁共振(MRI)或CT提供的图像是足够精度的，我们暂不考虑其精度问题对我们模型的影响。2 大脑组织是均匀的，且对射线的传播没有影响，即射线在其中传播能能量不会衰减。3 每条射线是一个理想的圆柱形，且各处能量密度相同，同时由于伽玛射线的波长很短且各放射源相互独立，因此认为相交时不发生干涉。4 从不同准直器出来的光线的能量密度是相同的，在我们的模型中认为是不可改变的。4 射线组的焦点是严格准确的，不考虑实际中存在的焦点偏差。5 各

6、个靶点的照射时间是相同的，不考虑因为此而带来的剂量问题。6 靶点是一个正球体。7 某一点累积的剂量和这些剂量同时作用于这一点所起的效果是相同的。三 名词解释1 有效剂量：能杀死肿瘤细胞的最低剂量。2 能量密度：单位体积介质中的射线的能量。3 IDL：isodose lines,等剂量曲线，在我们的模型中，说50%IDLs或选取50%等剂量线，是指用50%等剂量曲线的形状来作为靶点的参考形状。四 模型建立和算法设计我们的模型有三个重要的组成部分：1靶点的剂量问题。靶点剂量分布是有201条射线叠加而形成的，在我们的假设条件下通过计算讨论剂量的分布函数，多个靶点之间剂量的影响等问题，为下面的优化和评

7、价提供基础。2治疗计划的生成。我们用三维抽取骨架算法（3D Skeleton）来求出肿瘤的骨架，继而通过填充算法来对靶点进行操作并进行自评价，最终生成一个较优化的治疗方案。3对治疗计划的评价。主要考虑剂量的覆盖面，剂量的分布是否均匀，及周围组织受到得剂量大小等。（一）靶点的剂量分布靶点的附近剂量曲线根据我们的假设，射线呈圆柱形，不衰减，且不发生干涉。叠加效果如图一所示,根据这一假设，我们可以求出靶点周围的剂量分布。先考虑一条射线，设这条射线的能量密度为Po, 截面半径为a（可以为2mm/4mm/7mm/9mm）。一个以靶心的中心为球心，半径为R的球，截得的圆柱体积为Vcut(R)，根据二次积分

8、公式可以求得Vcut(R)=。 （1） 图一：光线叠加的结果其中所包含的能量P(R)=Vcut(R)*Po。下面同时考虑201条射线，因为201条射线射线同时照射一个中心，可以认为以该中心为球心的球面上各处的能量密度是相同的。半径为R(Ra)的球面上的能量密度为 （2）(其中Vsp(R)表示半径为R的球的体积，表示距靶点中心R处的能量密度)进一步可以求得 （3）从此公式我们可以得到剂量分布如图二所示，从中可以明显看出从1-0.5曲线下降得很厉害，以后逐渐变得平缓，因此选取不同的等剂量线作为有效剂量造成的对周围脑组织的影响是不同的，越低越大。为了对周围脑组织造成较小的影响，相对剂量水平一定要达到

9、50%以上，国际上较通用做法也是选取50%等剂量曲线作为有效剂量覆盖病灶周边，目的是利用立体定向放射外科靶点外剂量成梯度锐减的特点，减少正常脑组织的损伤。在考虑到一些其他因素后，我们将做出一点改进，这在下面的有效剂量的选择中将重点讨论。不同直径靶点对周围组织的影响的比较我们将其中两种准直器的靶点剂量分布做了比较，如图三所示，从中可以看出，18mm准直器的射线束直径大，剂量从聚焦中心向外周的递减比小直径射线束较弥散，对周围脑组织影响较大，术后并发症的危险性也相应增大，为了达到使用最小的靶点目的，势必要多选取较大直径的靶点，然而这样也会增大对周边脑组织的影响，因此为了达到总体的优化，应该从总体上达

10、到一个平衡，在我们的模型中也考虑了这一因素。图二：8mm准直器下靶点的剂量分布。右图，横坐标表示距靶点中心的距离，纵坐标为相对剂量。左图为靶点的等剂量线，最大圆为90%等剂量线。图三：图中心的空白处为靶点的大小，线条表示射线束。大直径的准直器靶点外射线束较弥散，而小直径准直器射线较集中，对靶点外正常组织影响小。两个靶点之间剂量的相互影响当所用到的靶点超过一个时，必然会涉及到两个靶点周边的重叠，而且这种影响是很显著的，因此有必要对这个问题进行深入的讨论。通过我们上面所得到的剂量分布公式，我们计算了两个2mm准直器靶点相互作用时他们中心线上的剂量分布，如图四所示，从中我们可以清楚看到，在70%-8

11、0%叠加时，等剂量线要比50%60%变化快，中央形成一个很大的陡坡，这不是我们想要的结果，我们要使最大剂量附近剂量变化尽可能平缓，即要使剂量的变化率尽可能小。还可以看出，70%-80%等剂量线叠加时，等剂量线将发生显著的变化（易位），这种变化是在我们进行计算机优化是应该考虑的。一般的模型先将靶点看成是刚性小球，根据病灶形状进行填充，然后再考虑剂量的影响，这往往使先前的工作变成徒劳，因为先前纯几何的讨论往往不能达到剂量分布的要求，一旦有多个靶点同时作用时，等剂量线将发生很大的易位，而调整起来又比较困难。我们的模型将这种因素融进去，做了一些改进，在作几何填充的同时将剂量因素考虑进去，这也是我们之所

12、以花很多精力来研究靶点的剂量问题。我们在剂量方面的模型从上面的讨论我么可以看到，因为靶点不是理想的刚性球，当多个靶点共同配合来治疗时，后加入的靶点将影响到在此以前的剂量分布，使以前的等剂量线发生图四：在不同IDLs（等剂量线）处叠加时中心线上的剂量分布，注意纵坐标的不同很大的变形，和病灶的外形就不会很好的匹配，因此常常需要一个迭代的算法。这种用几何一次性填充想要达到很好匹配的想法是不够合理的，付出的代价就是不断的迭代。我们想出一个解决的办法就是有新的靶点加入时尽量保持原来的等剂量线保持不变。要达到这个目的，我们想通过低剂量的叠加实现，即让每个靶点只在他们的低剂量区域开始进行叠加，这样对原来的剂

13、量分布造成的影响就小的多，而且能过满足有效剂量的要求。在我们的模型中将有效剂量定为70%以上，而不是通常的50%；我们将50%等剂量曲线作为我们覆盖病灶使达到最好匹配的参考形状（即我们在进行几何覆盖计算时使用的尺寸）。这样就可以有效的避免不断的迭代，同时在这一范围内里，剂量变化是较迅速的，因此不会对周边脑组织造成大的影响。这样做的同时还有一个好处，就是在这一段距离内，病灶的放置方式可以有一定程度的自由，我们认为两个靶点在50%-60%等剂量曲线组成的球壳区域内是可以有重叠的，因此我们的模型是有一定的弹性。同时若要考虑关键组织的情况，我们可以先通过对关键组织周围病灶进行匹配，通过这种办法来达到对

14、其的保护。（二）治疗方案生成要达到和病灶的最佳匹配，必须同时考虑靶点的大小和位置。最容易想到的就是将靶点放到病灶的中心，然后尽可能的靶点放进去，以达到最大的覆盖率，这种想法不是我们模型的思想，但求最大球的思路是可以肯定的。我们在查阅资料中发现一个这方面的一个有效的算法，就是三维骨架抽取（或称中轴变换）算法（3D Skeleton）6-10，这使得我们的想法有了实现的方法。这部分我们主要进行了三个主要的工作：（一）通过利用类似分形的算法来生成一个病灶。（二）利用距离地图（Distance Map/SED）的方法来抽象出病灶的骨架。（三）利用我们的算法来实现治疗方案的最后设计。病灶的生成我们的模型

15、是在一个50mmX50mmX50mm正方体内（共有625，000个像素），通过向里面放进一些随机产生的“种子”（像素），然后让它们根据自己的“邻居”的情况遵守一定规则来向各个方向自动“生长”，直到有一个方向达到边界。一些定义：我们认为像素是一系列的正方体，它们紧密排列。如果两个正方体共用一个面，我们称它们为面邻居（F-neighbors）。如果两个正方体共用一条棱，我们称它们为边邻居(E-neighbors)。如果两个正方体共用一个顶点，我们称它们为点邻居(V-neighbors)。一步：是指向某个方向“生长”一个像素。规则1如果一个像素没有邻居，它可以向四周的面邻居和边邻居“生长”一步。2如

16、果一个像周边已经被堵死，即在两步之内找不到能出去的路，则将周围填满。3如果一个像素有邻居但还没有堵死，则只向面邻居生长（如果原来已经有了，就不用再生成）。 图五：用我们算法生成的三维肿瘤三维图像的骨架抽取（3D Skeleton）一个图形的骨架（有的称中轴）是这样一个集合，这个集合中的元素是有一系列的点和它所对应的最大内圆（二维）或球（三维）组成的。一个最大内圆至少和图形的边界相切于两点，并且不能被别的最大内圆完全包含。骨架抽取在过去的年里已经取得了许多成果，在我们所找到的资料中，主要有两种主要算法（一）细化方法。即逐次去掉物体的边界点，同时必须保证区域的连通性不被破坏，最后保留的点组成物体的

17、中轴。这种方法多采用模板方式来实现。抽取的中轴可以保留物体的拓扑结构。（二）中轴转化法。即先对图像进行距离变换，找出每个像素与最近边界点的距离，将距离相对最大的点提取出来，作为中轴点。这种方法易于实现，效率高。但由此方法得到的中轴点离散分布可能不连通，会出现中轴断裂的情况。图六很清楚的说明了它们的区别图六：上面两个图形的骨架是通过细化方法生成的，中间的点或线就是所谓的中轴或骨架。下面的两个是通过中轴变换生成的。我们可以看出用中轴变换生成的骨架内圆半径有一个很大的近似连续变化的范围，这就就为我们用提供更大的选择空间。这在我们的模型中将用到。我们在算法中用第二种方法，主要是针对我们的问题并不需要过

18、多的保留拓扑结构，因为我们能用的球只有四种，按我们的模型，四种球可以在很小的范围上变化。即使这样，也不可能完全按曲率刚好匹配，所以我们需要在更大范围上有所选择。基于这一点和效率的考虑，我们采用算法复杂度较小的带符号的距离地图算法（8SSED, Sequential Euclidean Distance mapping algorithms）6812来实现。治疗方案设计以上算法生成的骨架相对应的最大内球，和我们的靶点唯一的区别就是，最大内球是可以任意叠加的，根据我们的模型，靶点只能在一定程度上重叠，因此我们要想出一种规则，来达到这个目的。基于上面算法生成的骨架，我们提出以下规则来对病灶进行覆盖（

19、填充），以达到对病灶的最大覆盖。定义：“角”和“分支点”： 如上页图七所示。分支的尽头： 如果将骨架看成是一棵树的话，那么根节点或叶子节点我们成为分支的尽头。规则：（一）靶点的中心在生成的中轴（骨架）上。（二）靶点只允许在低剂量区域（70%）重叠。（三）靶点从骨架的病灶的一个“角”开始填充，也就是从骨架的分支的尽头开始。（四）已经填充的区域不在下一次填充的考虑范围。（五）一次将病灶所有的“角”都填充好，然后再计算剩余部分的中轴，。具体算法（参照示意图:图七）：第一步：通过“带符号的距离地图”求出图形的骨架或称中轴。第二步：按我们定的规则进行填充，对于骨架的每一个分支的尽头，沿这条分支寻找，在不

20、通过分支点（包括分支点）的情况下，找到一个适合四种靶点之一且是相对最大靶点的位置并填充。如果没有找到适合靶点的位置，则可以超过分支点，类似第二步进行填充，直到找到合适位置。第三步：当所有的分支的尽头都已经填充好后，重新计算剩余部分的中轴，即将剩余部分看成是一个新问题再调用我们的算法。注意： 这里所说的剩余部分是指在原来图形中只除去已经填充进去的靶点部分。 图七_A （第一步） 图七_B（第二步1） 图七_C（第二步2） 图六：我们算法的示意图。括号写出了对应算法的步数，这只是我们算法开始的两次迭代。（图胚来源5） 图七_D（又第一步） 图七_E（又第二步）（三）治疗计划的评价我们最终得到的治疗

21、方案要达到一定的临床和生物学上一些要求。1剂量的变化率我们用等剂量线的疏密来表示。2有效剂量的覆盖面可以用70%等剂量曲面包围的体积占病灶总体积的百分比来表 示3对周边脑组织的影响我们的模型是选取一个最低相对剂量（比如20%），低于这个水平的剂量认为对脑组织没有影响。然后可以通过数值积分的方法求出。4优先考虑关键组织。5如果上面生成的治疗方案有小的偏差，可以通过适当的调整来实现，因为我们放置的时候就在不影响治疗效果的前提下留有一定的调整程度。五 模型及算法评价我们的模型的优点主要表现在以下几个方面：1我们对剂量分布进行了深入的分析，考虑了各个靶点之间的相互影响，及由于这种影响而造成的纯粹用几何

22、匹配所带来的一系列问题，并得到了一些结论。因此我们的模型在放置靶点的同时，也将靶点之间的影响问题同时考虑，制定了有效剂量曲面和用来覆盖病灶的等剂量曲面的分别考虑。从而避免了计算量很大的迭代，能够一步到位。2我们采用了在图形存储和压缩方面常用的骨架抽取思想，能够迅速抓住病灶的特征，为后面治疗计划的生成提供了有利的前提条件，虽然骨架抽取本身也是一个难题，但是在实际应用中已经有了一些很好的近似算法（线性的）和实际应用，这将为我们的模型提供一个有力的后盾，也将有力控制我们模型的复杂度，同时使我们的模型有更大的可操作性。3根据我们在二维图形上的尝试，我们的模型能够得到理想的结果。我们的模型也存在一些不足

23、和需要改进之处：1我们的模型将射线想象成完美的圆柱体，并且没有考虑到射线的衰减及身体组织对其的影响。这和实际是有一定出入的，从而也会造成治疗方案的偏差，但我们认为关于射线的形状并不能影响到我们模型的实用价值。只要适当调整有效剂量曲线和用来覆盖病灶的等剂量曲线就可以解决。2在实际中，靶点绝对剂量是可以通过射线的照射时间和堵住一些源来改变，这就对模型提出了更高的要求，并且可以更好的解决关键组织免遭破坏的更多方案,由于时间关系，我们没能够就这一点进一步研究。3根据我们的填充算法，抽取骨架的算法还可以更简单，因为我们是从病灶的“角”开始填充的，因此求出图形所有的骨架是一种浪费，理想状态时只要能确定靶点

24、的位置和大小就可以了，这是我们将来要简化的，将两者集中在一个算法里也不是没有可能的。4我们的遗憾之处就是我们发现骨架提取算法的时间较晚，并且其本身也相当复杂，对于三维图形更是如此。我们的三个主要算法：(1)病灶的生成算法，（2）骨架的提取算法（8SSED）,(3)病灶的填充算法 由于时间的限制，我们顺利的完成了第一个，其余我们只能完成部分，因此不能得到完整的结果，在写论文的时间，我的同伴正在调试，争取能够让程序正常运行，只要有充足的时间。参考文献：1Clinical Implementation of an Automated Planning System for Gamma Knife R

25、adiosurgery. David M.Shepard etc.2 Treatment Planning for Image-Guided Robotic Radiosurgery，Rhea Tombropoulos, John Adler, and Jean-Claude Latombe，Computer Science Department，Neurosurgery Department，Stanford University3 王恩敏 周良辅 潘力 不同准直器对大病灶中心剂量和周围脑组织受照容积影响的研究功能性和立体定向神经外科杂志 1996年第九卷第四期4A Study and Au

26、tomatic Solution for Multishot Treatment Planning for the gamma-knife Q.jackie Wu, Ph.D. and J.Daniel Bourland,Ph.D. Journal of Radiosurgery,Vol.3，No.2,20005Computerized Treatment Planning for Stereotactic Radiosurgery D.M.Shepard , University of Maryland School of Medicine, M.C.Ferris ，University o

27、f Wisconsin-Madison6 Wu Q, Bourand JD, Robb RA: Fast 3d medial axis transformation to reduce computation and complexity in radiosurgery treatment planning. SPIE2710:562-571,19967 Efficient Computation of A Simplified Medial Axis Mark Foskey, Ming C.Liu, and Dinesh Manocha Departments of Radiology an

28、d Computer Science, UNC-Chapel Hill forkey,lin,dmcs.unc.edu8Extraction of the Euclidean skeleton based on a connectivity criterion， Wai-Pak Choi, Kin-man Lam, Wan-Chi Sui ， Patterm Recognition 36(2003) 721-7299Fast and exact Signed Euclidean distance transformation with linear complexityO.Cuisenaire

29、 and B.Macq10Three-dimensional skeleton and centerline generation based on an approximate minimum distance field Yong Zhou,Arie Kaufman,Arthur W.Toga the Visual Computer (1998)14:303-31411. 韩西安杨海涛等基于切片图像的血管三维重建方法 装备指挥技术学院学报Vol.13 NO.3 June 200212. Q.Z.Ye, The signed Euclidean distance transform and its applications.in:Proceedings of the ninth International Conference on Pattern Recognition, Rome, 1988, pp.495-499.13Gamma Knife Dosimetry & Treatment Planning Jurgen Arndt,Karolinska Hospital Skockholm Sweden, AAPM 1999