肿瘤靶向治疗中的建模及仿真研究---硕士论文

河北工业大学硕士学位论文学位论文肿瘤靶向治疗中的建模及仿真研究河北工业大学硕士学位论文肿瘤靶向治疗中的建模及仿真研究摘 要目前，随着医疗技术的不断进步，要求医生在诊疗过程中尽可能做到高疗效、低损伤。“微创”治疗是未来肿瘤治疗发展的一大趋势，其目的在于以最小的创伤争取最好的治疗效果，其中以局部细胞灭活为主的靶向治疗，成为二十一世纪肿瘤治疗的主导方向。这种方法可以准确的杀灭靶区局部的肿瘤细胞，等体积靶向切除肿瘤，而最大限度的减少周围正常组织的损伤。因此，肿瘤的靶向治疗受到了国际、国内的普遍关注，并得到了较广泛的应用。但是，肿瘤的靶向治疗尚未规范，疗效很大程度取决于医生临床经验和治疗设备，不但很难取得稳定良好的疗效，而且极易引起医患纠纷。因此，实现科学规范的治疗方法是提高治疗疗效，减少经济损耗的重要手段。（在肿瘤的靶向治疗中，如何科学、规范、合理的选取并确定靶区，以及如何精确定位是治疗的关键。因此，本论文通过对肿瘤靶向治疗中的一种氩氦刀肿瘤靶向治疗系统进行分析和研究，）因此，我们试图借助先进的监测和控制手段规范肿瘤靶向治疗方法，改善疗效及对治疗效果进行科学的评估。本课题中我们针对计算机辅助肿瘤靶向治疗计划和疗效评估的关键技术进行了深入的研究和探讨，提出了靶区模型、氩氦刀模型和定位模型的概念，为实现计算机辅助治疗和治疗的科学化、规范化打下基础。因此，本论文的研究具有重要的理论意义和实用价值。研究针对肿瘤靶向治疗的需要，基于医学图像对肿瘤病灶、氩氦刀精确定位以及氩氦刀消融过程进行分析和建模研究。我们首先研究如何通过肿瘤 CT 图像的预处理，建立肿瘤病灶的三维模型；并考虑氩氦刀不同结构和使用条件，基于神经网络建立了描述氩氦刀冰球形成过程的数学模型；然后，研究了基于 MBS 理论建立氩氦刀靶向精确定位模型的方法，建立了精确定位数学模型。在此基础上，对肿瘤靶向治疗中的建模与仿真系统的开发进行了探索。关键词： 肿瘤，靶向治疗，氩氦刀，图像处理，数学模型，仿真河北工业大学硕士学位论文RESEARCH ON MODEL AND SIMULATION IN TUMOUR TARGETED THERAPYABSTRACTNowadays, with the continuous progress of medical technology, high curative effect and low damage are demanded in the process of diagnosis and treat. So, ” microtrauma” is a future major tendency of tumour treatment development. Its goal lies in striving for the best treatment result by the smallest wound, targeted treatment which is partial cell deactivation primarily becomes the leading direction of the 21th century tumours treatment. This method may accurately destroy target area partial tumour cell, the amount volume tumour is targeted excised, while reducing normal tissues damage to maximum limit. Therefore, the malignant tumour targeted therapy received universal attention internationally and at home, and obtained the widespread application.The malignant tumour treatment is not completely standard at present, which has seriously affected the curative effect. Realizing scientific and standard method of treatment is the important way to enhance curative effect and reduce the economy loses. In tumour targeted therapy, how scientifically, standardly, reasonably to select and define target area, as well as how to localization precisely are the key. Therefore, the paper analysis and research one of tumour targered therapy-cryocare tumor targeted therapy system, propose the concept of target area model, cryocare model and localization model. It builds the foundation for realizing computer assistant treatment and treatment is scientific and standard. Therefore, the paper research has important theory significance and practical value.The paper analysis and research model about tumour focus, cryocares localization precisely and the melt process of cryocare according to the need of cryocare tumour targeted therapy, based on tumour medical image.We firstly study how to establish three-dimension model of tumour focus through tumour CT images pretreatment; Considering cryocare different structure and exploitation conditions, establish mathematical model of cryocare iceball based on nerve network; then, study pinpoint localization model method of cryocare based on MBS theory, establish pinpoint localization mathematical model.On the basis, we carried on exploration to develop the system of model and simulation in tumour targeted therapy.KEY WORDS: tumour, targeted therapy, cryocare, image processing, mathematical model, simulationv肿瘤靶向治疗中的建模及仿真研究河北工业大学硕士学位论文目录第一章绪论 11-1 课题的提出和意义 1 1-2 国内外研究动态 1 1-3 该领域研究目前存在的主要问题 3 1-4 本论文的主要工作 3 1-4-1 论文的主要研究内容 3 1-4-2 论文的总体结构 3第二章 肿瘤靶区模型研究 42-1 引言 4 2-2 靶区建模的研究 4 2-2-1 基于体数据的建模 4 2-2-2 基于点云逆向工程的建模 5 2-3 靶区建模实例分析 10 2-3-1 基于体素建模实例分析 10 2-3-2 基于点云逆向工程建模实例分析 12 2-4 本章小结 17第三章 氩氦刀冰球模型的研究 183-1 引言 18 3-2 影响氩氦刀冰球的因素分析 19 3-3 氩氦刀冰球模型的建立 22 3-3-1 人工神经元网络介绍 23 3-3-2 BP 网络 24 3-3-3 氩氦刀冰球 BP 网络模型的设计 25 3-4 本章小结 28vii肿瘤靶向治疗中的建模及仿真研究第四章 氩氦刀精确定位模型研究 284-1 引言 28 4-2 用于复杂坐标系统定位计算的 MBS 理论 29 4-2-1 MBS 拓扑结构和低序体阵列 29 4-2-2 相邻体间坐标变换矩阵建立 31 4-2-3 典型体上给定点理想位置方程 34 4-3 基于 MBS 理论的氩氦刀精确定位建模方法研究 35 4-3-1 坐标系转换关系 35 4-3-2 坐标系的建立 37 4-3-3 定位模型的建立 38 4-4 本章小结 39第五章：基于 VRML 的氩氦刀肿瘤靶向治疗中的仿真研究 415-1 引言 41 5-2 关于仿真 41 5-2-1 仿真的定义 41 5-2-2 仿真的应用 41 5-2-3 仿真的步骤 41 5-3 关于 VRML42 5-3-1 VRML 简介 42 5-3-2 VRML 特点 43 5-3-3 VRML 基本节点编程 43 5-3-4 事件和路由 44 5-4 基于 VRML 的氩氦刀冰球生长与消融仿真44 5-5 本章小结 47第六章：用于肿瘤治疗的建模与仿真系统 486-1 引言 48 6-2 系统的开发 48 6-2-1 编程平台的选择 48 6-2-2 系统总体流程规划 48 6-2-3 系统组成及主要功能 49 6-3 本章小结 54河北工业大学硕士学位论文第七章：结论与展望 557-1 本课题研究的总结 55 7-2 本课题研究工作的展望 56参考文献 57致谢 59攻读学位期间所取得的相关科研成果 60ix肿瘤靶向治疗中的建模及仿真研究河北工业大学硕士学位论文第一章 绪论1-1 课题的提出和意义肿瘤特别是恶性肿瘤对人类健康的威胁很大，它和心血管疾患已成为医学上的两大难关，在全世界构成死亡原因的头两位。全世界 52 亿人口中，每年有约 700 万人新患癌症，每年约有 500 多万人死于癌症，几乎每 6 秒就有一名癌症患者死亡。我国目前每年平均有 150 万人新患癌症，每年约有 80 万人死于癌症，癌症的死亡率占总死亡原因的第十位。传统的肿瘤治疗方法有手术、放疗、化疗、中医中药治疗、免疫生物疗法等。手术是目前肿瘤治疗最主要的手段，其效果也最好，治愈率达到了 22%，而其他方法的治愈率都在 20%以下。近年来，随着医学技术的不断进步，要求医生在对患者进行诊疗的过程中，尽可能做到高疗效、低损伤。于是，“微创”治疗是未来肿瘤治疗的发展的一大趋势。其目的是在于以最小的创伤争取最好的治疗效果，其中以局部肿瘤细胞灭活为主的靶向治疗，成为二十一世纪治疗肿瘤的主导方向。这种治疗方法可以准确的灭杀靶区局部的肿瘤细胞，等体积靶向切除肿瘤，而最大限度的减少周围正常组织的损伤。因此，恶性肿瘤的靶向治疗受到了国际、国内的普遍关注，并得到了较广泛的应用1。但恶性肿瘤治疗目前尚未完全规范，这严重影响了治疗效果，也成为造成医患纠纷的隐患因素。实现科学规范的治疗方法，是提高治疗疗效，减少经济损耗的重要手段。为了达到靶向治疗的目的，精确定位技术是治疗的关键和瓶颈。如何科学、规范、合理的选取并确定靶区，建立精确的靶区模型和定位模型，是达到预定治疗效果的关键问题。因此，本论文通过对肿瘤治疗现状和前人研究成果的深入分析和研究提出了肿瘤靶向治疗中的建模及仿真研究的课题。肿瘤靶向治疗可以准确的杀灭癌细胞，最大限度的保护周围正常组织，以获得最好的治疗效果，尽可能的延长患者的生命。而其中的关键技术是靶区和定位模型的准确程度，本文通过对肿瘤靶向治疗中的一种氩氦刀肿瘤靶向治疗系统的研究，建立了靶区模型、氩氦刀模型和定位模型的模型，以便提高肿瘤治疗的效果。因此本论文的研究具有重大的理论意义和实用价值。1-2 国内外研究动态目前，国内外肿瘤专家采取了许多有效的措施，尽力提高肿瘤靶向治疗中的定位精度。由于现在的定位技术多依靠专家的经验，对病人的病体信息和治疗环境进行判断和决策。因此，定位不确定因素多，医生心理负担重。而且，除了诸多人为因素，还受病体的复杂性、环境的多变性以及不可预见性因素的影响。因此不能实现肿瘤靶向治疗中的精确定位，很难发挥靶向治疗的优势。现行的肿瘤靶向治疗中，由于缺乏对病灶的尺寸、形状、位置、状态、边界、内容等的精确确定，通常是参考 B 超、CT 或核磁共振等医学影响资料，根据医生以往的经验来判断病灶各种信息属性。因此，主客观因素的影响，都造成了定位方案的差异性。而随着医疗水平的不断发展，要求医生在对患者进行治疗的过程中，要尽可能的做到高疗效、低损伤。因此，今年来国内外很多医学专家及相关学科学者肿瘤靶向治疗中的定位技术研究，其关键技术已经成为国内外专家研究的前沿和热点2。靶区模型的建立主要是根据计算机断层扫描(CT)、核磁共振成像（MRI）、或者超声（US）得到的人体及其器官内部的二维断层图像序列，得到靶区模型。这属于三维重建的范畴，也属于逆向工程的范畴。这类问题是目前的一个研究热点问题，是一个多学科交叉的研究领域，是计算机图形学和图像处理在生物医学工程中的重要应用。它涉及数字图像处理、计算机图形学以及医学领域的相关知识。国外在三维1肿瘤靶向治疗中的建模及仿真研究重建方面的研究最多的是日本，其次是美国和英国。三维重建的研究在发达国家起步较早，研究也比较深入。上述几个国家的许多所大学以及其他科研机构都在进行该技术的相关研究。目前，医学可视化技术发展十分迅速，已经由单纯的辅助诊断治疗向教学研究以及辅助手术治疗方向发展。由医学图像可视化技术生成的三维图像直接引导手术过程，开发了计算机辅助立体定向导航系统 (Computer Assisted Stereotactic Syetem，CASS)，计算机辅助手术系统及模拟系统(Computer Assisted Design & Simulation System，CADS)等，例如，Philippe St-Jean 等人通过图像变换与融合操作，将已经标注了主要功能区的脑组织切片图与病人的脑部核共振图相匹配，并以此为基础生成病人的脑部三维图像，在其上规划手术路径，并可以在实际手术中引导手术过程3。Yoshinobu Sato 等人针对乳腺癌切除手术的特点，手术前采用超声波诊断，并确定病灶位置，并以超声波数据为原始数据生成三维图形，在手术过程中与摄像机所拍摄的图像结合，以此引导手术进程4。在国内目前工作主要集中于可视化算法研究，开发与之相关的医学可视化系统，并正在向辅助手术与模拟手术方向发展。浙江大学、清华大学、东南大学、中科院自动化所等均做了大量的研究，开发了一些试验系统。四川大学制造科学与工程学院的王德麾等人以一组人体头部层去式实体扫描切片为对象，对数字化人体的三维重建进行了研究。用基于颜色空间的方法提取目标图像，完成图像去噪和边缘检测工作；并用 VC+语言编程调用 matlab 引擎并记录目标像素的坐标值，获得点云数据文件；在 CATIA 环境下进行曲面拟合，完成三维数字化实体的重建5。清华大学郑卫国等人针对 CT 图像数据，提出了基于 CT 图像的快速成型数据建模方法。通过 CT 图像滤波、分割、等处理方法以及轮廓提取和采样等建模方法，得到快速成形数据借口 CLI 文件，进而进行原型的快速成形制造6。目前的定位手段大多是通过将定位面膜、体膜、胸膜固定，通过在病人体表划线，然后利用一些辅助定位装置以及 CT、MRI 图像才能算出病灶中心的深度及病灶在体内三维方向的具体轮廓及其方位7。基于 CT 的模拟定位系统在肿瘤靶向治疗中也有了简单的应用。它是一个与 CT 扫描机相连的实时三维立体治疗计划系统，由 CT 扫描系统、多层面 CT 影像显示控制台、3D 治疗计划系统、激光投影定位装置等部分组成，所有配置均由计算机工作站控制，并已成为一种标准的治疗设计工具，将广泛用于多个部位的肿瘤和其他病变的治疗计划与临床治疗。为了提高治疗效果，节省治疗时间，减少患者在治疗时的痛苦，迫切需要适合于肿瘤靶向治疗的精确定位技术和系统。目前的研究已经开始将肿瘤治疗中的定位方法与计算机常规技术方法相结合，研究靶向治疗中的计算机辅助定位。刘丽冰教授在全国第一届肿瘤靶向治疗会议上，提出了肿瘤靶向治疗中的精确定位技术，得到了与会专家和临床医生的高度评价。仿真是一门利用计算机软件模拟实际环境进行科学实验的技术，虚拟现实技术是利用计算机技术建立一种逼真的虚拟环境，它集成了计算机图形学、多媒体、人工智能、传感器、网络、并行处理等技术的最新发展成果。浏览者通过数字手套、立体头盔、立体眼镜和三维鼠标等传感器与计算机发生联系，最终产生一个拟人化的三维逼真的虚拟环境。目前，它与多媒体、网络技术并称为三大前景最好的计算机技术，在医学研究领域中得到越来越广泛的应用。从医学图像开始，发展到虚拟人体、虚拟医疗系统、虚拟实验室和药物研究。虚拟人体是在虚拟环境下实现人体科学计算可视化。美国国家医学图书馆于 1991 年委托克科罗拉多大学建立起一个男人和一个女人人体的全部解剖结构的数据库。一具男性尸体从头到脚作了 CT 和核磁共振扫描，男的间距 1mm，女的间距 0.3mm。然后以同样的间距作组织切片的数码相机摄影。所得数据共 56GB。这些数据称为 Visible Human 数据集，对于人体研究起了巨大的推动作用。虚拟医疗系统是指虚拟现实技术在临床上的应用系统，初期是一些模拟训练系统，寻找最佳的临床治疗方案，逐步发展成为临床应用系统。如外科手术模拟系统、遥控手术（远程医疗）、复杂医学数据库图像和康复医疗、轮椅虚拟等。外科手术手术模拟是典型的手术再现，用户戴上立体头盔，如同深入透视图像中学习人体解剖，同时可以用解剖刀和钳子等工具虚拟外科手术。外科手术工作站配置了三维监视器和三维摄影系统，以及带感觉输入设备的反应控制器和灵敏的压力反馈操纵杆。在遥控手术中心，医生只需对虚拟病人的模型进行手术，通过高速宽带网络将医生的动作传送到网络另一端的手术机2河北工业大学硕士学位论文器人，由机器人对病人实施手术。1-3 该领域研究目前存在的主要问题在肿瘤靶向治疗中存在的问题很多，例如如何建立精确的靶区模型，要求准确的获得靶区的边缘，这就涉及到图像分割和边缘检测，图像分割与边缘检测的研究具有重要的意义。结构分析、运动分析、三维可视化等后续操作，以及图像引导手术、肿瘤放射治疗、治疗评估等应用研究都假设已对图像做了准确分割，或者说都是以图像分割为基础的。医学图像分割是正常组织和病变组织的三维重建、定量分析等后续操作的基础，也是临床医学应用的瓶颈，分割的准确性对医生判断疾病的真实情况并做出正确的诊断计划至关重要。由于医学图像中具有极其复杂的多样性和复杂性，加上目前医学影像设备(CT、MRI、PET 等)成像技术的特点，使得医学图像存在一定的噪声，图像中目标物体部分边缘也有可能局部不清晰，这使得医学图像分割更加困难。因此我们需要在分割之前采取必要的预处理，以减少图像中的噪声。在氩氦刀肿瘤靶向治疗中，冷冻靶区设计不合理、冷冻范围小、残留率高是导致手术疗效差的直接原因。这种情况往往出现在开展冷冻手术的初始阶段，对氩氦刀和冷冻理念的理解不深所致。缺乏三维设计的理念，冷刀排列、介入点过分依赖手术者个人主观设计，也是导致残留率偏高的重要原因。由于不同冷冻技术的设计和临床参数不同，治疗中的临床疗效不同。冷冻靶点定位穿刺不准确致使靶区脱离病灶，使癌灶仅部分冷冻或未完全冷冻直接影响疗效。B 超定位因其一维定位成像较难反映肿瘤的三维立体结构，如果 B 超医师定位技术不熟练，手术医师对 B 超图像理解较差，将直接影响冷冻效果。B 超定位是要多角度探查，三维评价。采用 CT、MRI 定位引导，可使定位、引导穿刺、和术中监测的准确性得到了显著提高，但目前国内开放性 CT、MRI 较少，适时监控存在一定困难，手术中需要影像医生的密切配合。治疗中要注意建立三维适形靶区消融概念。1-4 本论文的主要工作1-4-1 论文的主要研究内容本论文通过对氩氦刀肿瘤靶向治疗系统的研究，根据目前靶向肿瘤治疗中存在的一些问题，如靶区的立体结构认识不准、定位不准等，提出了氩氦刀靶向治疗肿瘤中靶区模型、氩氦刀模型和定位模型的概念，为氩氦刀肿瘤靶向治疗计算机辅助精确定位打下基础。通过对医学 CT 图像的处理建立靶区模型；通过对氩氦刀冰球数据的分析，建立氩氦刀冰球模型；通过对氩氦刀坐标系、靶区坐标系和人体定位坐标系的分析，建立定位模型。根据建立的氩氦刀冰球模型和靶区模型、定位模型进行仿真。最后开发氩氦刀肿瘤靶向治疗的建模与仿真系统。1-4-2 论文的总体结构第一章：绪论第二章：肿瘤靶区模型研究第三章：氩氦刀冰球建模方法研究第四章：氩氦刀精确定位模型研究第五章：基于 VRML 的氩氦刀肿瘤靶向治疗中的仿真研究第六章：氩氦刀肿瘤靶向治疗中建模与仿真系统开发第七章：结论与展望3肿瘤靶向治疗中的建模及仿真研究第二章 肿瘤靶区模型研究2-1 引言所谓靶区，是指氩氦刀靶向治疗中的病体区域，也是通常的感兴趣区域。随着计算机断层扫描(CT： Computed Tomography)，核磁共振成像（MRI：Magnetic Resonance Imaging），超声（US：Ultrasonography）等医学成像技术的发展，人们可以得到人体及其内部器官的二维数字断层图像序列。但是，二维断层图像只是表达某一截面的解剖信息，仅由二维断层图像上某些解剖部位进行简单的坐标叠加，不能准确给出三维影像，造成靶区定位的失真与畸变8。为了提高医疗诊断和治疗规划的准确性与科学性，我们有必要建立靶区的三维模型，由二维断层图像序列转变成为三维的结构和形态，可以计算它的体积甚至重量。医生可以根据这些影像资料中所显示的肿瘤及周围组织结构的三维立体解剖关系，初步讨论制定一合理、完整的冷冻手术方案，并制定出术中、术后可能发生并发症的应对措施。2-2 靶区建模的研究在利用计算机对客观事物进行分析和研究时，需要建立相应得模型来表示实际的或抽象的对象或现象，这个过程称为建模。利用断层序列图像建模有两种思路：一种是基于体数据的建模，一种是基于点云逆向工程的建模。基于体数据的方法又分为两大类：一类是基于表面的方法，一类是直接体视法。2-2-1 基于体数据的建模1.物体表面建模三维体数据中蕴藏着丰富细腻的物体结构。物体表面重建的任务就是要恢复出蕴藏在三维体数据中的物体的几何信息，以便更加方便的对物体进行定量分析、处理和显示。而要从三维体数据中重建物体表面，首先必须根据采样值的分布变化确定物体的边界，然后用适当的数学方法来描述物体表面9。表面表示是表示三维物体形状最基本的方法，它可以提供三维物体形状的全面信息。根据重建过程所操作的对象来分，重建方法可以分为从轮廓重建物体表面（边界轮廓线表示）和从体素重建物体表面（表面曲面表示）两种方法。1.1 从轮廓重建物体表面断层扫描数据广泛的存在于医学、生物、地质无损探伤等应用领域，是一种最简单的三维标量场。从表面轮廓进行物体的三维重建不仅在体视化中有重要的意义，同时也在计算机辅助设计上有着重要的应用价值，因为在 CAD 中也常用一组放样轮廓来描述形体。一般各断层间是相互平行的，每一断层与实体的交线就是实体在该断层上的轮廓线，即是二维平面上一条封闭的无自交的等值线。断层数据的表面重构是从一系列断层面上的轮廓线中推导出相应实体的空间几何结构。整个表面重构的过程分成两步：拓扑重构和几何重构。拓扑重构获得实体的拓扑表示，是对三维断层数据集中 每一断层上的轮廓线进行分类，确定各轮廓线所属实体及同一实体不同断层上轮廓线上的对应关系。几何重构是在拓扑重构的基础上建立对象的几何表示。几何重构可分为面向曲面的重构和面向体的重构。面向曲面的重构中常用的是三角片表面重构；面向体的重构有 Delaunay 四面体重构和六面体重构。在医学及其他领域，切片大多数重建是以断层图像作为处理的对象，首先在每个切片图像上抽取出物体的轮廓，然后用轮廓拼接构造通过这组轮廓的物体表面物体表面。这方面的代表性工作是 Kepple 和 Boissonanat 的“从平面轮廓重建三维形体”。4河北工业大学硕士学位论文最初的表面重建方法采用基于轮廓线的描述方式，即在断层图像中，通过手工或自动方式实现目标轮廓的确定性分割，然后用各层的轮廓线“堆砌”在一起表示感兴趣的边界10。这种轮廓线表示方法简单，数据量小。除了轮廓线表示物体之外，也可以由轮廓重建物体的表面来表示。Boissonna 提出了一种基于表面轮廓的 Delaunay 三角形方法。这种方法首先对给定的离散数据点进行三角剖分，可分为对测量数据投影域的剖分和在空间直接剖分两种类型，目标是使散乱的数据点在空间连接成一个最优的三角网格，尽量接近 Delaunay 三角化；第二步计算三角网格边界条件，构造初始三角曲面；第三步构造整体 G1 连续的差值曲面。但是这种方法得出的只是分片光滑的表面，采用 Lin 的从轮廓出发的 B 样条插值重建算法，可以得到整体光滑的表面。1.2 从体素重建物体表面分割出的待建组织区域以体素形式记录，我们可以对分割提取出的区域构建出其对应的组织或器官的三维几何模型。其方法是在提取得区域中，根据给定阈值，构建出表达重建组织或器官的等值面。等值面的构造是从体数据恢复物体三维几何描述的常用方法之一。如果我们把体数据看成是某个空间区域内关于某种物理属性的采样集合，非采样点上的值以其临近采样点上的采样值的插值来估计，则该空间区域内所有具有某一个相同值的点的集合将定义一个或多个曲面，称之为等值面。因为不同的物质具有不同的物理属性，因此可以选定适当的值来定义等值面，该等值面表示不同物质的交界。也就是说，一个适当定义的等值面可以代表某种物质的表面。最早的体素级重建方法叫做立方块法(Cuberille)，它是用边界体素的拟合等值面，即把边界体素中相互重合的面去掉，只把不重合的面连接起来近似表示等值面11，这种方法的特点是算法简单易行，便于并行处理，因为对每个体素的处理都是独立的；主要问题是出现严重的走样，显示图像给人一种“块状”感觉，尤其是物体边界处锯齿走形特别醒目，而且显示粗糙，不能很好的显示物体的细节。而 Loresen 等人于 1978 年提出的移动立方体法（Marching Cubes），是最有影响的等值面构造方法，一直沿用至今。最初的 Marching Cubes 算法不能保证三角片所构成的等值面的拓扑一致性，会造成等值面上出现空隙。 M.J.Durst 首先提出了 MC 算法中的二义性，后来许多人在 Lorensen 方法的基础上做了许多改进。移动四面体法(Marching Tetrahedra)是在 MC 算法的基础上发展起来的，该方法首先将立方体元剖分成四面体，然后在其中构造等值面。由于四面体是最简单的多面体，其它类型的多面体都可以剖分为四面体。2. 直接体视法体绘制由于直接研究光线通过体数据场与体素的相互关系，无需构造中间面，体素的许多细节信息得以保留，结果的保真性大为提高。2.1 投影法（Projection）首先根据视点确定每一体素的可见性优先级，然后按优先级由低到高或由高到低的次序将所有体素投影到二维像平面上，在投影过程中，利用光学中的透明公式计算当前颜色与阻光度，依投影顺序（即体素可见性优先级）的不同，投影法分为从前至后（Front-to-Back）算法与从后至前（Back-to-Front）算法。一般说来，前一种算法运算速度快，但除需一个颜色缓存区外，还需要一个阻光度缓存区；后一种算法仅需一个颜色缓存区，并在执行过程中产生不同层面的图像，有助医生更好的理解医学图像。2.2 光线跟踪法（Ray-Casting）该方法是在体数据进行分类后， 从象空间的每一体素出发，根据设定的方法反射一条光线，在其穿过各个切片组成的体域的过程中，等间距的进行二次采样，由每个二次采样点 8 个邻域体素用三次线形插值方法得到采样点的颜色和阻光度值，依据光照模型求出各采样点的光亮度值，从而得到三维数据图像。2-2-2 基于点云逆向工程的建模我们是基于医学图像建立靶区的模型，要提取到靶区的点云信息，一般须经过医学图像的预处理。这些图像预处理通常包括滤波、分割、边缘检测以及轮廓提取等方法。5肿瘤靶向治疗中的建模及仿真研究1. 医学图像预处理滤波处理的基本方法图像在生成和传输的过程中常受到各种噪声的干扰和影响，使图像质量下降，图像模糊，特征淹没，给分析带来困难。为了抑制噪声，改善图像质量，必须对图像进行平滑处理，这可在空域或频域中进行。在平滑噪声时应尽量不损害图像中的边沿和各种细节。对于滤除图像中的噪声，人们已经提出了很多的方法。通常，将数字图像的平滑技术划分为两类。一类是全局处理，即对噪声图像的整体或大的块进行校正以得到平滑的图像。例如在变换中使用 Wiener 滤波、最小二乘滤波等。使用这些技术需要知道信号和噪声的统计模型。但对于大多数图像而言，人们不知道或不可能用简单的随机过程精确的描述统计模型，而且，这些技术计算量也相当大。另一类平滑技术是对噪声图像使用局部算子。当对某一像素进行平滑处理时，仅对它的局部小邻域的一些像素加以运算，其优点是计算效率高，而且可以多个像素并行处理12。1.1.1 中值滤波如果既要消除噪声又要保持图像的细节，可以使用中值（median）滤波。中值滤波是一种有效的非线性滤波，常用于消除随机脉冲噪声。它的基本思想是：在图像上滑动一个含有奇数个象素的窗口，对该窗口所覆盖象素的灰度按大小进行排序，处在灰度序列中间的那个灰度值称为中值，用它来代替窗口中心所对应象素的灰度。中值滤波法在实际运算过程中并不需要图像的统计特性，这带来了不少方便，但是其也有一些缺陷，对一些细节特别多，特别是点、线、尖顶细节多的图像不宜采用该方法。1.1.2 邻域平均法对于图像中的每一个象素，取一个以它为中心的区域，用该区域内各象素灰度的加权平均值取代该象素的灰度值，这就是邻域平均法。具体的做法是取一个方形区域，称为平滑窗口（Window）或掩膜（Mask），它是权值的二维阵列。滤波过程是用窗口在图像上滑动，窗口中心对着的象素根据下面式（2.1）来更新其灰度值。当每个象素都被扫描一次之后，对一幅图像的平滑就完成了。这是实平面上对图像进行平滑的一种最简单的方法。设 f (i, j) 是一幅待平滑的图像，平滑窗口为 W 的大小为 (2N +1) (2N +1) ，则平滑后的图像可以表示为NNg(i, j) = wuv f (i + u, j + v)u =- N v=- N（2.1）NN wuvu =- N u =-v其中， wuv 是权值。我们可以对窗口内的权值进行归一化，即使权值之和等于一。NN wuv = 1（2.2）u =-N u =-v于是，公式（2.1）只剩下分子项。1.1.3 保持边缘滤波法前面介绍的两种滤波方法在抑制噪声的同时都或多或少地钝化了边缘，使图像清晰度下降。我们希望在滤除噪声的同时保持边缘不被平滑掉。保持边缘的滤波器就是为了实现这个目的而设计的。如果在滤波之前能够检测出边缘，那么只要对那些没有边缘的图像区域进行平滑。然而，边缘是很难可靠地检测出来的。尤其是在噪声严重的情况下，更是无法区分哪些是真实的边缘哪些是噪声造成的伪边缘。设计保持边缘滤波器的出发点是：对应于同一个目标的象素具有相近的灰度。对于每一个象素，逐个检测它的每一个邻域（包含该象素），其中，邻域的灰度方差最小的那个邻域是最平滑的，其内部的6河北工业大学硕士学位论文象素基本上可以认为是同属一类。因此，用该邻域的象素灰度均值取代该象素原来的灰度值。1.2 图像分割与边缘检测由于医学图像通常对比度较低，组织特征的可变性、不同组织之间或者组织和病灶之间边界的模糊性，使医学图像技术的智能化受到了限制。图像分析专家正面临着如何在计算机辅助下，精确提取影像中包含的信息问题。医学图像处理除了包括普通数字图像处理通常的图像滤波去噪、图像增强、图像复原、图像变换之外，其中最重要的要算图像分割和边缘检测。图像分割就是把图像中具有特殊含义的不同区域分开来，这些区域使互不相交的每一个区域都满足特定区域的一致性。它是图像处理与图像分析中的一个经典问题，也是制约医学图像处理中其它相关技术发展和应用的瓶颈13。目前针对各种具体问题已经提出了许多不同的图像分割算法，对图像分割的效果也有很好的分析结论。但是由于图像分割问题所面向领域的特殊性，至今尚未得到圆满的、具有普适性的解决方法14。边缘检测技术的目的是在有噪声的图像中确出目标的边界，边缘检测的好坏，会直接影响到后续的治疗过程。它在医学图像匹配、肿瘤病灶确定、造影血管检测、冠心病诊断、左心室边缘抽出等方面占有举足轻重的地位。它还广泛的用于颅脑三维重建前的边缘抽取，尘肺的自动侦测，脑灰质脑白质的抽取。边缘至今为止在数学上没有准确的定义，一般的文献都以局部极值点或灰度发生巨变的点即为边缘。它可以粗略的分为阶跃边缘与屋顶边缘。医学图像分割在临床诊断、病理分析以及治疗方面具有重要意义，具体表现在以下几个方面：（1）图像分割的结果常用于生物医学图像的分析，如不同形式图像的配准、融合、解剖结构的测量，获取先验知识用于图像重建以及心脏运动的跟踪等15。（2）用于测量人体器官、组织或病灶的体积 。通过对这些体积在治疗前后的定量测量和分析，可以帮助医生进行诊断、预后和制定或修改对病人的治疗方案16（3）用于医学图像的 3D 重建方面，便于可视化，外科手术方案的制定和仿真，病理研究 ，药物疗效的评估，解剖参考以及放疗计划中的 3D 定位等17。分割问题的困难在于图像数据的模糊和噪声的干扰。实际图像中景物各异，具体问题具体分析，需要根据实际情况选择适合的方法。分割结果的好坏或者正确与否，目前还没有一个统一的评价判断准则，分割的好坏必须从分割的效果和实际的应用场景来判断。不过在人类研究图像的历史中，还是积累了很多经典的图像分割方法。虽然这些分割方法不适合于所有类型的图像分割，但是这些方法却是图像分割进一步发展的基础。事实上，现代一些分割算法恰恰是从经典的分割方法衍生出来的。目前国内外广泛应用的医学图像分割方法有基于边缘的图像分割、基于区域的图像分割、结合特定理论工具的方法。基于边缘的图像分割假设图像分割结果的某个子区域在原来图像中一定会有边缘存在；基于区域的图像分割假设图像分割结果的某个子区域一定会有相同的性质，而不同区域的像素则没有共同的性质。这两种方法都有缺点和优点，有的学者也试图把两者结合起来进行图像分割。基于边缘的分割方法分为并行边界分割和串行边界分割。并行边界分割主要是以灰度的一阶和二阶导数作为边缘检测的依据，为此人们设计出了各种边缘检测算子，如梯度算子、Sobel 算子、Laplacian 算子、Candy算子和 Marr-Hiddreth 算子等，跟踪检测出的边缘点以构成边缘。这些边缘检测器对边缘灰度值过渡比较尖锐且噪声较小的图象，可取得较好的效果。但对于边缘复杂、采光不均匀的图像，效果不够理想，且对噪声较敏感。串行边界分割主要是边界跟踪。基于区域的图像分割方法也分为并行区域分割和串行区域分割。并行区域分割方法是指采用并行的方法对目标区域进行检测实现图像分割的方法，主要有：阈值分割和聚类。串行区域分割是指采用串行处理策略对目标区域直接进行检测来实现分割的方法。串行区域分割的特点是整个处理过程可以分解为顺序的多个步骤依次进行。串行区域分割一般可以分为两种方法：一是区域生长，二是分裂合并。区域生长是指从某个像素出发，按照一定的准则，逐步加入临近像素，当满足一定的条件时，区域生长终止。分裂合并差不多是区域生长的逆过程：从整个图像出发，不断分裂得到各个子区域，然后再把前景区域合并，实现目标提取。结合特定理论工具的方法主要有可变模型法，人工神经网络法，基于小波变化的方法，基于分型的分割方法，基于统计学的分割方法和基于数学形态学方法。比如基于三维可视化系统结合 Fast Marching 算法和 Watershed 变换的医学图像分割方法，能得到快速、准确的分割结果18。人体器官是一个高精度庞大的结构，所建模型还应考虑临7肿瘤靶向治疗中的建模及仿真研究床应用方面的因素以及某些特殊部位的个体差异，针对不同的生理组织应采用不同的分割方法19。根据不同的图像、不同的组织要选用不同的分割方法。在 CT 图像中常见的分割任务和相应得分割方法见表 2.1表 2.1CT 图像的常用分割方法Table 2.1 The General Dividing Methods of CT Image任务常用分割方法备注骨骼的分割阈值法、区域生长法、马尔科夫随机场、分割结果常用于三形变模型、模糊区域生长维重建胸部扫描图像分类统计聚类、区域生长与水线算法结合、例如用于支气管树区域生长与模糊逻辑结合、形变模型状结构的三维重建肺的分割模型拟合、形变模型腹部动脉瘤、胃部、心形变模型脏的分割在 MRI 图像中常见的分割任务和相应的分割方法见表 2.2表 2.2 MRI 图像的常用分割方法Table 2.2 The General Dividing Methods of MRI Image任务常用分割方法备注通常在初始分割步骤中允许丢失一些由于相邻脑部组织和非脑部组织强度值提取脑室脑部组织，然后用形态学滤波器、形变模重叠，使脑室提取比较型或图谱引导方法找到组织困难使用 T1加权或多普数分类器方法、聚类方法、神经网络方法、据分割脑灰质、白质和马尔科夫随机场脑脊髓液胼胝体、皮层下脑部结构、海马状突起：分割特定的脑部组织形变模型、图谱引导法肿瘤和损伤：神经网络，图谱引导法，线性滤波、模糊区域生长和形变模型心脏图像分割区域生长、阈值法、形变模型如描绘出心脏左心室膝盖图像分割马尔科夫随机场模型、形变模型核 磁 血 管 造 影 图马尔科夫随机场模型、形变模型、阈值（MRA）法8河北工业大学硕士学位论文2.逆向工程建模随着测量设备、计算机辅助几何设计(CAGD)技术的发展，逆向工程技术正在成为研究与应用的热点19,20。基于某种可获得的现有产品或物理模型的测量数据来重构设计概念，实现对原样件的逼近或修改，从而达到复制该产品或获得新产品的过程称为“反求工程”或“逆向工程”（Reverse engineering）。例如，在精致的轿车车身设计或人脸类雕塑曲面模型的 CAD 建模中，常先用油泥制作物理模型，再作三维型值点采样，由所获得的数据来建立 CAD 数字化模型；在医学图像可视化中也常用 CT 切片来得到人体骨骼或脏器表面的三维数据点，进而构建其三维数字化模型21,22,23。首先是测量数据结构，获取三维几何数据的方法有多种24，其中能获取外表面的（从外部可见的）位置数据的方法有坐标测量机法、激光三角形法、模式光法、图像分析法等，能获取物体内部以内表面数据的方法有工业 ICT、医学 CT、MRI 扫描以及层析扫描等方法。这些方法在精度和速度上都存在很大差异，但从建模角度出发，我们更关心的是数据反映的特征及其结构组织。下面列举了建模的几种主要的数据结构：（1）散乱点 无论采用何种测量设备，如果使用了多次测量的策略，导致数据在分布上不遵循一定的规则，无明显的几何关系则称为散乱点。（2）栅格点 测量点在扫描平面上呈等距栅格状分布，在激光扫描仪、光学数字化仪中较多采用。数据点具有确定的位置关系。（3）等值线 在 CT、层析设备中出现，数据点由某一方向的一系列面族与目标表面相交形成，如地形等高线、头骨截面线。对于我们基于 CT 图像的建模而言，每层的轮廓线的生成还需要对原始图像进行复杂的预处理，包括边缘提取、分析内外轮廓等工作。（4）特征点和线 规则面越多，就专门测量相关的特征点，用以形成平面、圆柱、圆锥、球等基本几何元素；对于某些由运动学方程形成的复杂曲面，可测量其母线及其运动轨迹；对于那些完全没有生成规则的曲面还可以先按照交线进行分割，在每块区域沿某个方向测量一些脊线，作为建模的依据。数据的不同结构和特点决定建模过程的难易程度以及将要采取的造型策略。在 4 种不同的数据结构中，包含的模型几何信息呈逐步增多的趋势，对建模而言采用靠后的测量数据，建模相对也要方便一些。建模主要有四种方法：线框造型、曲面造型、实体造型和参数化特征造型，因为测量数据获取的主要是物体表面测量点的坐标，所以曲面造型方法是逆向工程的主要思想。目前实体造型方法也逐渐得到重视25。常用的曲面造型手段主要有：（1）放样法(lofting) 对数据云取不同的截面点，构成截面线，再由截面线构成曲面。构造的截面线必须满足同向、节点分布相近的原则才能保证最终形成的曲面不扭曲；（2）旋转法 对于旋转形成类的曲面，首先切取并拟合一根母线，然后旋转成曲面；（3）包封曲面法 在对由四边形区域分割而成的数据点集合，首先进行测量点的参数化，然后拟合出满足方程 S=BP 的参数曲面。其中 S 为曲面型值点坐标矩阵，B 为曲面基函数系数矩阵，P 为控制点矩阵。在构造单个曲面中需要考虑的主要问题是实现型值点与拟合曲面之间的距离最小的优化目标，采用非线性最小平方和算法对型值点的参数值进一步优化，并得到新的构造曲面是有效的方法。如果优化后的误差距离仍然不满足预定误差要求，就只有采用增加曲面节点等方法改变曲面形式，重新拟合计算。针对大规模散乱数据的曲面重构问题，一些研究人员也提出一种两步法的拟合方案，首先用某种简单的插值法，如 Shepard 插值等快速构造插值曲面，然后在此曲面的基础上根据曲面局部的曲率值自适应得抽取采样点，生成矩阵拓扑网络，然后再运用曲面拟合方法得到 NURBS 曲面26。曲面模型描述物体的几何形状，没有体的信息。表面类型有规则曲面、自由曲面、复杂曲面 3 种，在构造曲面前需要做好不同曲面的区域分割工作，通常即使是全自由曲面构成的零件也需要进行分割，以保证局部细节的拟合精度。分割的方式主要有自动分割和手工分割两种，自动分割算法的方法主要有如下 3 种27,28：（1）四叉树法 将整个曲面首先看成一个整体，拟合成单一曲面，然后检验误差是否满足精度要求，不满足就一分为四，然后对每一部分重复上述过程，直至每个子曲面都能满足精度的要求。9肿瘤靶向治疗中的建模及仿真研究（2）基于特征线的方法 需要寻找曲面的尖锐过渡、曲率极值点以及曲面对称中心等特征作为曲面片之间的边界，这样得到的特征网格比前一种方法得到的网格更能反映物体的结构特征，也适合于下一步的曲面模型建立。（3）基于面的方法 这是一种结合曲面拟合的分割方法，然后按由近及远的顺序向外扩展，同时不断检验曲面的拟合程度，直至误差超出预定要求时停止。在进行层析数据的重构时常常采用这类方法。2-3 靶区建模实例分析2-3-1 基于体素建模实例分析我们以 Matlab 为基础采用一组 CT 切片图像进行体素的体绘制说明建模的过程。体绘制实质上是光线在媒介中传播的近似，媒介就是我们所用的体数据，其类型与结构随所研究的具体问题的不同而不同，在流体力学、气象学、地质学等领域，体数据常常是非均匀几何结构。如结构化网格或曲线化网格数据，非结构网格数据及混合型数据。典型的网格化结构是规则网格 G ( xi , yj , zk) 。xi = i L , i = 0,1,., Nxyj = j M , j = 0,1,.Ny（.）zk = k N , k = 0,1,.NzL，M，N 分别为 x，y，z 方向上节点之间的间距。Nx，Ny，Nz 分别为 x，y，z 方向上的网格点个数，相邻的 8 个网格节点构成一个长方体体元。有多种体绘制方法，其中应用最广的是光线投射整体显示方法。当光线穿过体数据时，它被反射、吸收或者被散射。体绘制仅仅绘制出视觉上可见的部分，对那些大量的透明的体素部分不予处理。由于光线在固定的环境中传播，能很快的达到稳定，因此我们认为光线在给定的区域内在一段时间内是稳定的。首先根据体数据 f ( x, y , z) 的值对每一数据点赋予不透明度及颜色值，再根据各数据点所在处的梯度及光照模型计算出各数据点的光照强度，然后将投射到图像平面中同一像素点的各数据的半透明彩色值综合在一起，最后形成图像。Matlab 提供了丰富的矩阵和图像计算工具，我们可以利用 Matlab 工具箱进行 CT 切片图像的建模。具体步骤如下：（1）数据输入对采集到的 n 幅 CT 切片 xy 尺寸的图像进行三维体数据集 D 的构造，得到 x y n 的矩阵。I(:,:,1)=imread(E:课题课题资料数据3.0mmjjj.bmp29.bmp); I(:,:,2)=imread(E:课题课题资料数据3.0mmjjj.bmp30.bmp); I(:,:,3)=imread(E:课题课题资料数据3.0mmjjj.bmp31.bmp); I(:,:,4)=imread(E:课题课题资料数据3.0mmjjj.bmp32.bmp); I(:,:,5)=imread(E:课题课题资料数据3.0mmjjj.bmp33.bmp); I(:,:,6)=imread(E:课题课题资料数据3.0mmjjj.bmp34.bmp); I(:,:,7)=imread(E:课题课题资料数据3.0mmjjj.bmp35.bmp); I(:,:,8)=imread(E:课题课题资料数据3.0mmjjj.bmp36.bmp); I(:,:,9)=imread(E:课题课题资料数据3.0mmjjj.bmp37.bmp); I(:,:,10)=imread(E:课题课题资料数据3.0mmjjj.bmp38.bmp); I(:,:,11)=imread(E:课题课题资料数据3.0mmjjj.bmp39.bmp); I(:,:,12)=imread(E:课题课题资料数据3.0mmjjj.bmp40.bmp);10河北工业大学硕士学位论文I(:,:,13)=imread(E:课题课题资料数据3.0mmjjj.bmp41.bmp);I(:,:,14)=imread(E:课题课题资料数据3.0mmjjj.bmp42.bmp);I(:,:,15)=imread(E:课题课题资料数据3.0mmjjj.bmp43.bmp);I(:,:,16)=imread(E:课题课题资料数据3.0mmjjj.bmp44.bmp);I(:,:,17)=imread(E:课题课题资料数据3.0mmjjj.bmp45.bmp);I(:,:,18)=imread(E:课题课题资料数据3.0mmjjj.bmp46.bmp);I(:,:,19)=imread(E:课题课题资料数据3.0mmjjj.bmp47.bmp); I(:,:,20)=imread(E:课题课题资料数据3.0mmjjj.bmp48.bmp);（2） 数据预处理和平滑构造的体数据集数量大，可根据实际情况对体数据集进行预处理，以降低所处理的数据量，进行平滑处理。Ds = smooth3(I);（3）体数据的绘制利用 isosurface 函数计算体数据集在显示平面累计投影。isosurface(Ds,5)；（4）构造结果图像碎片利用 patch 函数来对图像子区域进行分类，可以定义结果图像的颜色，光线等信息。hiso = patch(isosurface(Ds,5),.FaceColor,1,.75,.65,.EdgeColor,none);hcap = patch(isocaps(I,5),.FaceColor,interp,.EdgeColor,none);（5）计算碎片的法线方向利用 isonormals 函数计算图形子区域顶部法线方向。isonormals(Ds,hiso)；（6）计算三维图像的几何边界利用 isocaps 函数计算三维图像的几何边界。isocaps(I,5)；（7）图像显示利用 view,daspect,colormap,camlight,lighting 等函数显示图像。其中 view 函数定义观察者视角， daspect 定义 x 轴、y 轴和 z 轴的比例，colormap 定义图像的颜色集，lighting 定义显现图像的光线阴影。view(45,30)axis tightdaspect(1,1,.4)lightangle(45,30);set(gcf,Renderer,zbuffer);lighting phong利用 Matlab 进行建模得到如图.结果：11肿瘤靶向治疗中的建模及仿真研究图.Matlab 三维重建结果Fig. 2.1 3D Reconstruction of Matlab由此可见，利用 Matlab 进行 CT 切片图像的三维建模的效果还是比较理想的，但是由于数据量达，运算时间长，不能满足实时性的要求。2-3-2 基于点云逆向工程建模实例分析医学 CT 图像对于不同的组织和病变呈现为不同的像素灰度分布，可以帮助诊断者了解器官的立体形状，观察、判断、分析和定位病变组织。要根据一组医学图像进行靶区建模，首先要把图像中的靶区提取出来，我们需要将每张 CT 图像中靶区的轮廓提取出来，由于 CT 图像中存在干扰，因此必须在进行图像分割以前进行图像的预处理，包括滤波去噪以滤除图像中的干扰和噪声，增强图像中的有用信息。然后进行图像分割和边缘检测，得到靶区的轮廓。由于我们可以知道图像像素间距与实际空间距离比以及断层图像之间的层厚，可以经过简单处理得道与图像像素点一一对应的像素空间坐标信息29，并以一定的文件格式保存。读入图像数据文件图像像素的空间信息，重建一个逻辑的三维空间场，原来 CT 序列图像中的各个像素以结构网格的形式在其中排列，网格节点的值节点就对应像素的灰度值。因此，我们以第一张 CT 切片图像的左上角像素点为空间坐标点的原点，以从左向右为 X 轴正向，从上向下为 Y 轴正向，垂直屏幕向外为 Z 轴正向，建立坐标系。因此我们就可以知道这组 CT 切片图像中每个像素的灰度信息和坐标信息，当然我们也就得到了靶区轮廓的坐标信息。我们以一组 CT 切片图像为例，介绍靶区轮廓提取的大致过程：首先我们对图像进行预处理以增强图像的有用信息；然后对图像进行分割，得到目标组织的像素区域；由于进行图像分割后得到的像素不连续，存在干扰，为了滤除干扰信息，得到更好的效果，我们对图像进行了开操作处理，接着检测目标边缘，最后将得到的边缘轮廓以一定的文件格式保存。流程如图 2.2 所示：12河北工业大学硕士学位论文开始提示用户设置分割阈值对图像进行预处理对图像进行分割对图像进行开操作对图像进行边缘检测合并所有轮廓文件结束图.图象处理流程图Fig. 2.2 Flow Chart of Image Processing大致步骤如下：图 2.3 原始 CT 图像Fig. 2.3 Original CT Image13肿瘤靶向治疗中的建模及仿真研究（1）滤波。由于 CT 图像的形成不可避免的会引入不同的噪声，所以我们需要在对图像进行平滑和去噪预处理。通常的滤波方法主要有均值滤波、中值滤波和高斯滤波。中值滤波可以有效地去除尖峰信号，削平振动信号，而且对阶梯和陡的下降沿则照样保留，因此能有效地抑制噪声，同时保留有效信号。因此我们选用中值滤波。图 2.4 中值滤波结果Fig. 2.4 Result of Median Filter（2）图像分割。因为不同的组织和器官在 CT 图像呈现不同的灰度值，因此可以用阈值分割法将图像分为不同灰度的目标区域和背景区域。图 2.5 选取阈值为 50 时的分割结果Fig. 2.5 Segmentation Result of Threshold 5014河北工业大学硕士学位论文图 2.6 选取阈值为 100 时的分割结果Fig. 2.6 Segmentation result of Threshold 100从上面的图中可以看出，选取阈值为 100 时周围存在干扰，因此我们可以通过开操作来滤除干扰。（3）开操作。在图像二值化之后，获得了目标和背景两块灰度区域。但是由于图像本身分辨率导致图像二值化的结果不一定能完全区分开目标和背景，从而在背景和目标交界的地方会产生一些过渡像素，使得目标区域的边界处出现参差不齐的现象。在进行轮廓检测时容易对实际轮廓的检测造成影响，所以引入了图像开操作的图像形态运算子。图像开操作包括图像腐蚀和图像膨胀两步操作，属于图像处理中形态运算子之一。通过控制开操作进行腐蚀或者膨胀时的像素范围，可控制一些小的或者狭长物体的消除，达到平滑感兴趣的较大图像区域边界的目的，同时不会显著改变感兴趣的区域形态。图 2.7 开运算结果Fig. 2.7 Result of Open Operation（4）边缘检测。15肿瘤靶向治疗中的建模及仿真研究图 2.8 边缘检测结果Fig.2.8 Result of Edge Detect图 2.9 边缘检测结果Fig. 2.9 Result of Edge Detect（5）保存轮廓文件目标图像的轮廓是以图像形式存在的，必须将轮廓上的点转换成矢量形式才能使用通用的 CAD 软件进行三维重建工作。矢量化的主要工作有：1. 设定坐标系：采用笛卡儿直角坐标系来确定点的空间位置，以图片矩阵的行号为 Y 坐标，以列号为 X 坐标，Z 坐标为切片的序号，间距为 0.5mm。理论上，以这组图片中的任何一个固定的点都可以设定为坐标原点，文中以图片左上角为坐标原点。2. 确定长度单位转换系数：直接从图片上得到的点 X,Y 的坐标是像素的个数，要转换成标准的长度单位（毫米）则需乘上一个系数。借助照相时设置的标尺，可以用标尺上的刻度值与像素数的比例来确定单位转换的系数。若一厘米的标尺长度对应 n 个像素，则长度转换系数为 10/n（毫米/像素）。只要找出每个像素所在的行列数，再乘以单位转换的系数，就可以确定出点的空间坐标，进而用标准长度单位表达模型实际大小。如果还要对模型进行放大或者缩小，在此系数上再乘以缩放系数就可以。3. 程序算法设计：算法的关键在于追踪起点和追踪法则的确定。在目标轮廓图片中，所有白色的点都16河北工业大学硕士学位论文是边缘，因此只要搜索到第一个值为 1 的点就可以当作追踪的起始点。以此点为中心，依次检索它的 8 连通区域的点，判断是否为 1，则该点为有效点。（6）三维实体模型的重建用点云数据反映的三维模型是由离散的点组成的点集，不符合人们的视觉习惯，达不到工程使用的要求；但是点云是包含了大量反映了物体表面的数字信息；是目前逆向工程中唯一的生成实体的方法。由点云数据转换成三维实体一般要通过复杂的数学运算，包括：1.误差点处理：去除离散的噪音点；2.数据插补：补齐因意外丢失的特征点；3.数据平滑：消除数据毛刺，提高点云的质量；4.数据网格化：建立离散点，尤其是层之间数据的联系；5.曲面拟合：由网格化的点构建实体曲面；2-4 本章小结本章首先研究了用于靶区建模的各种方法，根据一组医学 CT 图象进行靶区建模的方法有：基于体数据三维重建和基于点云逆向工程三维重建，分析了其建模的优缺点和过程。然后，通过两组 CT 图像进行了模型的建立，一组用 Matlab 实现了基于体数据三维重建中的体绘制，一组通过图象处理得到了靶区轮廓的点云信息，进行了基于点云逆向工程的三维重建。17肿瘤靶向治疗中的建模及仿真研究第三章 氩氦刀冰球模型的研究3-1 引言Cryocare 氩氦刀靶向手术系统配备有 4 或 8 个氩氦刀探头（细针）。控制器能单独控制或个冷热绝缘氩氦刀探头。氩氦刀探头的直径为 1.58 毫米，目前临床使用的探头直径分别为 2、3、5、8 毫米，当输出高压常温氩气时，氩气快速超低温制冷技术可使氩气在刀尖内急速释放，沿刀尖前端向后，由内向外，形成可控制的倒梨状冷冻区，使靶区病变组织内温度在几十秒内降至-130-175,使靶区组织消融坏死。B 超、CT、MRI、X 线监测可以看到清楚的低密度冰球靶区的形状，氩氦靶向治疗协作组专家的临床资料证实，CT 显示的消融靶区与氩氦刀在肺癌中形成的冰球相符合（如图 3.1 图 3.2 所示）。冰球的形成受到多种因素的影响，我们分析了影响冰球生长的种种因素，并分析了它们之间的关系，以此为依据建立了氩氦刀冰球模型，以便准确的控制氩氦刀冰球的形成，为精确的切除肿瘤，保护正常组织打下基础。图.氩氦刀形成的冰球Fig. 3.1 The Iceball of Cryocare18河北工业大学硕士学位论文图.氩氦刀形成的冰球Fig. 3.2 The Iceball of Cryocare3-2 影响氩氦刀冰球的因素分析根据氩氦刀靶向治疗协作组专家的资料，在临床治疗中冷冻消融靶区的分布受一下几种因素的影响：（1）氩氦刀探头直径的大小2mm、3mm、5mm、8mm 的氩氦刀探头可分别产生最大长径范围为 23cm、56cm、78cm、910cm 的消融靶区（冰球）。临床及基础研究表明冷冻靶区的组织消融温度范围氩氦刀消融温度范围为-38到-140分布区域，氩氦刀冰球靠近边缘 8mm 以内为不可逆消融区。因而临床临床治疗时要充分考虑靶区边缘的设定，根据病人情况常规冷冻区应该超过消融靶区的 12cm。根据不同器官肿瘤情况，部分学者提出在肝癌、肺癌治疗时，如果局部病情允许冰球边缘可以超过肿瘤边缘 13cm。临床治疗时根据病人肿瘤靶区的大小选择合适的氩氦刀探头。Endocare 临床医学专家建议：当一个氩氦刀探头的消融靶区不能完全覆盖肿瘤组织时，可以根据肿瘤大小及形状，采用二或三个氩氦刀探头适形组合应用（如图 3.3 所示）图.多探头组合示意图Fig. 3.3 The Sketch Map of Multi-Probe Compounding19肿瘤靶向治疗中的建模及仿真研究（2）氩气氦气的输出压力美国 Cryocare 氩氦刀靶向手术系统配备氩气容器的最大压力为 6000psi（41.4 兆帕），产生理想冰球冷冻效应所需要的工作压力是 2000psi，当压力低于 2000psi 时将影响冰球的形状、大小及冰球的形成速度。如果氩气压力不能达到 2000psi 冷冻灭活的工作压力，或者氩气达不到理想的纯度要求，可能造成临床靶区治疗计划达不到理想的治疗目的，也是低温治疗临床出现残留的常见原因。（3）冷冻时间冷冻时间无疑是影响冰球生长的因素之一，随着冷冻时间的增加，冰球也逐渐增大，但从理论上推测，冰球不会无限制的增大，当到一定程度，即冷热交换达到一个动态平衡时，冰球就不会再增大，即应有一冷冻时间的极限。根据氩氦刀专家的观察结果分析表明，在冷冻 5 分钟和 10 分钟时，以及冷冻 10分钟和 15 分钟时冰球的增加明显，具有统计学差异，20 分钟时的冰球较 15 分钟时仍有增大，但增幅已很微弱，无统计学差异。（4）冷冻温度冷冻温度也是一个影响冰球生长的因素，一般来说，冰球冷冻温度越低，冰球增长就越缓慢，反之，冰球增长就会快一些。为了具体准确的研究冰球的生长和氩氦刀直径、输出压力、冷冻时间和冷冻温度的关系，我们在相关医疗单位的帮助下做了很多的实验。表 3.1 是 2004 年 9 月 16 号在珠江医院肿瘤中心所作的离体猪肝冷冻试验，设备采用 Endocare CryocareTM Surgical System 冷冻消融治疗机、直径 2、3、5mm 冷冻刀各一把、武汉江堤气体厂生产“江气牌”氩气 2 瓶（纯度99.999%、压力5000Psi）表 3.1 2mm 刀具在相同时间、不同压力条件下所形成的冰球Table 3.1 The Iceball of 2mm Subject to the Same Time But Different Pressure刀体直径 mm气体压力 kPa冷冻时间分秒冷冻温度冰球大小 mm24500501373.84.5235005111382.74.522500501352.54.0215005101352.34.0表 3.22、3、5mm 刀具在相同时间、相同压力条件下所形成的冰球Table 3.2 The Iceball of 2mm,3mm,5mm Subject to the Same Time Period and Pressure刀体直径 mm气体压力 kPa冷冻时间分秒冷冻温度冰球大小 mm254004500501373.84.5354004500551303.55.7553004700521553.76.52450036001001363.45.520河北工业大学硕士学位论文3450036001041283.76.35470038001021514.58.02360027001501354.05.53360027001521284.56.85380030001531505.38.022003300024002001295.07.15200为了研究冷冻时间与冰球大小的关系，从而为临床确定最佳冷冻时间及冷冻靶区的设计提供理论依据宁波解放军 113 医院肿瘤科以 18 例接受氩氦刀治疗的患者为例，CT 实时监测，分别测量 5 分钟、10分钟、15 分钟、20 分钟时的冰球长径与横径，进行统计分析。表 3.3 2mm 氩氦刀冷冻 5、10、15、20 分钟所形成的冰球Table 3.3 The Iceball of 2mm Subject to 5,10,15,20 Minutes病冷冰球大小（cm）例肿 瘤冻5 分钟10 分钟15 分钟20 分钟编类 型编号号纵径横径纵径横径纵径横径纵径横径1恶性神经13.82.34.93.15.63.95.63.9鞘瘤23.92.55.03.25.94.25.94.22肝转移癌33.62.14.42.44.22.54.92.443.72.44.12.94.23.04.73.33肺癌53.62.34.52.75.03.65.13.765.33.15.43.25.63.55.63.64肝癌73.81.73.92.14.12.44.22.583.31.94.12.34.32.54.32.65肝癌93.52.14.72.94.93.450.3.5103.32.04.83.04.93.65.13.8113.62.35.13.35.23.65.43.8123.62.45.23.55.54.65.64.7133.22.34.63.55.23.65.33.9143.42.44.73.65.34.05.34.36肺癌154.82.35.12.85.43.25.53.3164.82.55.03.15.53.45.63.57肺癌175.12.65.23.05.23.25.23.3185.12.85.33.25.43.55.43.58肺癌194.52.55.23.35.33.55.33.6204.32.65.13.15.33.55.33.521肿瘤靶向治疗中的建模及仿真研究9肝癌214.02.55.63.46.24.56.54.8224.13.26.14.27.15.77.15.7233.82.55.33.45.54.25.94.4243.82.65.84.16.95.66.95.710肺癌253.72.44.53.04.83.05.03.1264.02.64.83.55.84.56.24.911肺癌274.22.95.03.85.94.66.24.812肉瘤284.92.35.53.66.54.26.64.3294.53.05.33.56.24.16.24.213肺癌303.82.34.32.94.83.24.83.214肺癌313.42.14.73.14.93.54.93.6323.32.24.52.84.93.55.03.515肺癌334.62.55.63.95.84.25.84.3344.22.35.43.65.94.15.94.116肺癌354.12.05.42.65.52.95.62.917肺癌363.82.14.72.74.83.44.93.5373.51.94.53.04.73.44.83.618肺癌383.42.24.12.74.83.24.93.3表 3.4统计结果：平均值Table 3.4 Statistical Result:Average5 分钟10 分钟15 分钟20 分钟纵径4.44.85.15.1横径2.43.03.63.7我们因此可以得出结论：氩氦刀冷冻时间随时间延长，冰球逐渐增大，其中在 5 分钟到 15 分钟时增大较为明显，此时冰球大小为 5.1*3.6*3.6，15 分钟到 20 分钟增幅趋缓。根据上面的试验和数据，我们知道氩氦刀冰球的生长受到多种因素的影响，就试验数据来看，影响的主要因素是刀径、压力、时间、温度等，它们之间有着非常复杂的关系，因此我们的任务就是要找到它们之间的这种关系。3-3 氩氦刀冰球模型的建立（1）简化模型22河北工业大学硕士学位论文要建立氩氦刀的冰球模型，往往需要作一定的简化处理。我们知道一般单个氩氦刀探针所形成的冰球呈一个规则的倒梨形状，我们可以近似将其看成一个圆球，圆心即为探针的顶点。因此，在建立坐标系时，我们以该点为圆心，建立坐标系，以进刀方向为 Z 轴的正方向，过圆心垂直于 Z 轴的平面为 XY 平面（如图 3.4 所示）。探针冰球XY 平面Z 轴图 3.4 简单氩氦刀坐标系和模型Fig. 3.4 A Simple Cryocare Model & Coordinate System（2）精确模型上面建立的模型作了很大的简化处理，这样建立起来的数学模型往往不能完全真实地反映出实际情况。氩氦刀的生长环境复杂多变，干扰多，不确定因素多，它的冰球模型是严重非线性的。因此，氩氦刀冰球模型与其环境一起构成了一个十分复杂的系统，其中包括许多作用因素和环节，它们之间相互影响、相互耦合，使得系统的输入输出之间呈现十分复杂的相互关系。它的各个环节、各个部分都存在许多非线性因素、不确定因素和无法描述的因素，要用数学表达式进行精确的描述，是十分困难的，甚至是根本不可能的。氩氦刀的冰球模型非常复杂，无法得到精确的数学模型。氩氦刀冰球模型的辨识方法主要有两种：一是通过理论分析的方法确定模型的结构形式，然后采用卡尔曼（Kalman）滤波等参数辨识方法确定模型的结构参数，从而确定系统的数学模型；二是利用神经网络能无限逼近任意有理函数的特性，通过网络的自学习逐步逼近系统模型，从而得到系统模型。通常第一种方法用于线性系统建模，而第二种方法则用于非线性系统建模。针对氩氦刀冰球模型的特点，本文采用非线性系统辨识的方法（人工智能识别方法）30,31进行建模。利用神经网络的自学习、自适应功能，通过系统的输入输出数据，直接辨识其输入与输出之间的关系，从而达到对这一复杂系统进行整体特性的辨识。3-3-1 人工神经元网络介绍人工神经网路的结构和工作机理基本上是以人脑的组织结构（即大脑神经元网络）和活动规律为背景的，它采用物理可实现的器件或计算机来模仿生物体中神经网络的某些结构和功能，并反过来用于工程界和其他领域，其着眼点是采纳生物体神经细胞网络中可利用的部分来克服目前计算机或其他系统不能解决的问题，如学习、识别、控制、专家系统等。它是由大量的、简单的处理单元（称为神经元）广泛的互相连接而形成的复杂网络系统，反映了人脑功能的许多基本特征，是一个高度复杂的非线性动力学系统。神经网络具有大规模并行、分布式存储和处理、自组织、自适应和自学习能力，特别适合处理需要同时考虑许多因素和条件的、不精确的和模糊的信息处理问题。它具有高度的非线性映射能力，事先不需要假设输出变量与输入变量之间的关系，而是通过样本的学习，实现输入与输出之间的非线性映射。在最抽象的层次上，神经网络可以被当作是一个黑箱，数据从一边输入，经过神经网络处理后产生23肿瘤靶向治疗中的建模及仿真研究相应的输出。输出与目标值比较后，产生的误差用于调整神经网络内部的连接权重。如图 3.5 显示了神经网络的功能图。图.神经网络功能图Fig.3.5 Function Chart of Nerve Network目前已有数十种人工神经网络模型。人工神经网络包括前向网络、反馈网络和自组织网络。其中前向网络包括 BP 网、感知器、径向基函数（RBF）网络和 GMDH 网络。反馈网络包括 Hopfield 网络、 Cohen-Grossberg 网络、盒中脑（BSB）模型、双向联想记忆（BAM）网络、回归 BP 网络、Boltzmann机网络。自组织网络包括自适应共振理论(ART)、自组织特征映射、CPN 映射模型。3-3-2 BP 网络BP 网络，即反向传播网络（Back-Propagation Network），是人工神经网络实际应用中使用最广泛的一种，8090 的人工神经网络模型是采用 BP 网络或它的变化形式32。它是前向网络的核心，体现了人工神经网络最精华的部分，它的应用也最为成熟。辽宁工程技术大学的白瑞才、李建刚、郭旭颖等人利用神经网络工具箱开发基于图形用户界面的 BP 网络方法和设计步骤，以露天矿卡车路段行程时间的预测，进行实例仿真试验，说明了该方法的可用性和便利性33。南京航空航天大学的张臣等人将 BP 神经网络的理论和算法应用于球头刀具铣削力建模的研究中。采用 LM 算法建立了铣削力预测的神经网络模型，模型中考虑了影响铣削力的加工参数，选取铣削力试验数据对神经网络模型进行训练，用训练好的神经网络模型对铣削力进行仿真。仿真结果表明，用 BP 神经网络方法建立的铣削力模型能够对铣削力进行准确的预测34。海军工程大学的戴余良等人分析了潜艇在风浪中运动的复杂性、不确定性和运动模型的严重非线性等特点，难以建立准确的数学模型，因此提出了利用神经网络的非线性映射能力，进行潜艇运动模型辨识模型的新方法，通过仿真试验验证了利用 BP 网络建立潜艇在风浪中运动模型的可行性35。上海交通大学的王强等人利用人工神经元网络实现函数逼近36。因此本文选择 BP 网络作为氩氦刀冰球模型。标准的 BP 网络模型，由三层神经元组成37：接受输入信号的单元层称为输入层，输出信号的单元层称为输出层，不直接与输入输出发生联系的单元称为中间层（又称隐含层）。各层次之间的神经元通过权重形成全互联系，各层次内的神经元之间则没有连接。形式如图 3.6 所示：24河北工业大学硕士学位论文图.BP 网络模型Fig. 3.6 BP network model输入层的神经元将输入信息 X i 分配到隐含层的各神经元，隐含层中的神经元 j 的输入加权和为W ji Xi ，其中Wji 为输入信息 X i 与神经元 j 之间连接权值，神经元 j 的输出是这一加权和的函数，即Y j = f ( W ji Xi )(3.1)式中： f ( ) 为作用函数，它一般为 S 型函数，即f ( s ) = 1/(1 + e-s )(3.2)在学习算法中，规定在神经元 i 和 j 之间连接权值的变化为DW ji =hdi Xi(3.3)式中：h 为学习率，di 为一个因子，它随神经元所处的位置不同而不同。为了加快 BP 算法的学习过程，通常在上述算式中增加一个动量项，即DW ji ( k + 1) = hd i X i + mDW ji ( k)(3.4)式中： DW ji ( k +1) 、 DW ji ( k) 分别为 K+1 和 K 时刻的权值的变化量， m 为动量系数。3-3-3 氩氦刀冰球 BP 网络模型的设计（1）网络结构确定理论上已经证明：具有偏差和至少一个 S 型隐含层加上一个线性输出层的网络，能够逼近任何有理函数。根据这一结论及我们的实际使用经验，通常只需具有一个隐含层的三层网络便可实现大多数应用。而且实践证明，对许多实际问题，一个中间层已经足够了。随着中间层的增加，训练时间会急剧增加。这是因为中间层越多，误差传播过程计算就越复杂，使训练时间急剧增加。且中间层增加后，局部极小误差也会增加，使网络在训练过程中，易陷入局部极小误差而无法摆脱。但是，在 BP 神经网络的训练过程中发现，当采用一个中间层时，需要用较多的处理单元，这时如果采用两个中间层，每层的处理单元就会大大减少，反而会取得较好的训练效果。因此，在建立多层 BP 神经网络时，首先考虑选用单隐含层；在这种情况下，如果增加了处理单元数还不能得到满意结果，可以增加隐含层数目，但应减少总的处理单元数。因此本文选取三层 BP 网络结构，即网络由输入层、一个隐含层和输出层三层组成。输入量分别为刀径、输出压力、冷冻时间和冷冻温度，输出量为冰球的横径和纵径。所以本文选取的 BP 网络输入层有 4 个节点，输出层有 2 个节点。采用适当的中间层处理单元数是非常重要的，单元数太少，网络难以处理较复杂的问题；而中间层处理单元数目过多，将使网络训练时间急剧增加，而且过多的处理单元易使网络训练过度，使网络难以分辨数据所代表的真正模式。隐含层节点的选取至今没有统一的规定，只能根据经验和反复的训练选取一个最合适的节点数，根据样本训练精度和学习效率的要求，我25肿瘤靶向治疗中的建模及仿真研究们经过反复试验，最终隐含层选取了 4 个节点。因此氩氦刀冰球模型为 4-4-2 网络结构。隐含层传递函数为非线性 logsig，输出层为线性 purelin，训练函数为 trainlm。（2）训练样本集的选择神经网络部分训练样本集列表如表 3.1、表 3.2 所示。通过修正的 BP 算法利用样本对神经网络进行训练，训练好的 BP 网络作为氩氦刀冰球模型的参数模型。（3）输入输出数据预处理为了加快网络的训练速度，往往对神经网络的输入输出进行一定的预处理。而归一化处理是使用最多的预处理方法之一。它将每组数据都量化到范围0，1内，对于输入向量 p 和目标输出向量 t 进行归一化处理的语句为：pn,minp,maxp,tn,mint,maxt=premnmx(p,t)。网络训练是采用归一化之后的数据，即 net=train(net,pn,tn)，训练结束后可以使用函数 postmnmx 将输出还原成最初目标所用的单位：an=sim(net,pn)；a=postmnmx(an,mint,maxt)。当用训练好的网络对新数据 pnew 进行预测时，也应作相应的预处理：pnew=tramnmx(pnew,minp,maxp)；anew=sim(net,pnew)；anew=postmnmx(anew,mint,maxt)。（4）训练数据的导入方法根据数据的多少，数据文件的格式等获取样本的方法有：1） 直接输入数据：通过采用元素列表方式输入。适用于样本数目较少时。本文采用这种数据输入方式。2） 创建数据文件：通过 MATLAB 提供的装载数据函数，从数据文件中读取。3） 输入向导读取：通过输入向导(Import Wizard)从文件或剪贴板中读取数据。该方法不适合从M-File 文件中读取数据。（5）网络的建立训练通过建立 M-File 文件，将以下代码输入，构建一个神经网络。p=222223523 5235;4500 350025001500495049505000 405040504250 3150 31503400;55 55 55 510 101015 1515;-137-138 -135-135 -137-130 -155 -136 -128-151 -135 -128 -150;t=3.82.72.52.33.83.53.7 3.4 3.7 4.54.04.55.3;4.54.54.04.04.55.76.5 5.5 6.3 8.05.56.88.0;net=newff(minmax(p),4,2,logsig,purelin,trainlm);net.trainParam.show = 500;net.trainParam.epochs =10000;net.trainParam.goal = 0.002;pn,minp,maxp,tn,mint,maxt = premnmx(p,t);net=train(net,pn,tn);%net,tr=train(net,p,t);%a = sim(net,p)%b=sim(net,2;4000;8;-140)an = sim(net,pn)a = postmnmx(an,mint,maxt)训练该网络，得到如图.所示结果：26河北工业大学硕士学位论文图 3.7 训练结果Fig. 3.7 Training result因此可以看出网络经过 1065 次学习达到收敛，达到了指定的精度和速度要求，完成了网络训练的过程。（6）网络的输出表 3.5 BP 网络输出Table 3.5 BP network output刀体直径 mm气体压力 kPa冷冻时间分秒冷冻温度冰球大小 mm24500501373.75744.4353235005111382.67704.364522500501352.47374.1523215005101352.36203.9414254004500501373.84984.6084354004500551303.47115.530027肿瘤靶向治疗中的建模及仿真研究553004700521553.76126.58292450036001001363.40595.53613450036001041283.71776.50085470038001021514.43097.84722360027001501353.96445.48183360027001521284.46356.78065380030001531505.35818.0247与目标输出相比较可以看出，所得输出已经非常接近目标输出，因此，训练所得的神经网络模型是有效的。3-4 本章小结氩氦刀治疗肿瘤的过程中，冰球的生长规律是影响手术计划的重要因素，根据分析，冰球生长因素主要受到刀径、工作压力、冷冻时间和冷冻温度的影响，而它们之间绝对不是简单的线性关系，而是非常复杂的非线性关系。因此要找到它们之间的规律就必须借助于人工神经网络，而 BP 网络是人工神经网络的精华，它具有非常强的非线性映射能力，可以完成任何复杂函数的拟合和逼近，因此我们用 Matlab 神经网络工具箱建立了 BP 神经网络，完成了氩氦刀冰球模型的建立。第四章 氩氦刀精确定位模型研究4-1 引言氩氦刀肿瘤靶向治疗中，准确的定位是取得好的治疗效果的关键。目前的定位方法是：依据 B 超、 CT 及 X 线影像学所显示的肿瘤大小、形态及与周围脏器、组织结构的关系，经全面比较，仔细测算确定氩氦刀介入瘤体的层面以及在同一层面和不同层面内介入氩氦刀的型号、数量和方位。并最终确定氩氦刀头端进入瘤体内的靶点。如肿瘤较大可用以 3mm 冷刀为主的多刀组合冷冻方案，周围情况允许时，选用 5mm 冷到效果更好。设计原则要求多刀组合冷冻所形成的冰球应尽可能将瘤组织包容其内。冰球冷冻范围应大于肿瘤边缘 1cm。肿瘤临近重要组织结构部位时，应该确定测温探针的位置以监控靶区的界28河北工业大学硕士学位论文限。确定瘤体消融靶点后，确定体表皮肤穿刺进针点。皮肤穿刺点选择以手术操作方便，易于呼吸监测管理，穿刺针经过正常组织的路径最短为基本原则。靶点与体表进针点所连直线即为穿刺针、氩氦刀进入靶点所经过的路径。路径内应避开神经、主支气管、心包及大血管。所有设计的体表进针点和瘤内的靶点确定后，在 CT 和 B 超定位片上分别测量显示出各点的进针方向、进针深度和进针角度。并根据定位的坐标数据，在进针点皮肤上粘贴一金属或塑料标记，再次经扫描显示标记物所在的方位与设计要求一致时，则标记物所在点极为穿刺进针点。这是目前人工定位的方法。目前的研究已经开始将肿瘤治疗中的定位方法与计算机常规技术方法相结合，研究靶向治疗中的计算机辅助定位。所谓计算机靶向精确定位技术就是充分利用医学影像的数字信息，建立肿瘤靶向治疗中的全定位数学模型，通过基于智能算法的决策系统，实现肿瘤靶向治疗中的精确定位。本章通过研究氩氦刀肿瘤靶向治疗系统，建立了三个坐标系，分别为氩氦刀坐标系、靶区坐标系和人体定位坐标系。其中氩氦刀坐标系和靶区坐标系是通过计算机建立的虚拟坐标系，在治疗过程中处于人体内部，并没有也不可能实际标记出来，是不可见的坐标系。人体定位坐标系是可以直接得到人体定位的坐标系，是可见的，可以直接用于治疗过程中。我们要研究的就是这几个坐标系之间的关系，已达到用通过可见的坐标系控制不可见的坐标系的目的。多体理论是解决多坐标系复杂系统的有力工具，因此我们应用多体理论建立了氩氦刀定位模型。4-2 用于复杂坐标系统定位计算的 MBS 理论多体系统（MBS）是指多个刚体或柔体通过某种方式联结而成的复杂机械系统38，是对工程实际中出现的复杂工程对象的高度概括和抽象，通过概括和抽象使人更广泛、更深刻地把握事物本质，从而更好的指导人们对工程对象进行研究、设计和分析。以 J.Wittenberrg、A.A.Shabana 以及 T.R.Kane， R.L.Huston 等为代表的学派，在多体系统的理论和方法的研究方面做出了卓越的贡献，特别是凯恩（Kane）发展了一种新的动力学的普遍方法，在 Kane 方程的基础上，休斯敦（Huston）建立了一套完整的算法，包括引入低续体阵列描述系统的拓扑结构，用欧拉参数描述体的方位，将动力学方程系数表达为偏速度和偏角速度及其导数的形式，是一种面向计算机编程的设计方法。该理论方法具有很强的概39,40括性、通用性和系统性，已用于机器人学、机构学以及航天技术等多个领域。4-2-1 MBS 拓扑结构和低序体阵列多个物体通过某种特定的形式联接起来构成 MBS 系统。基于 MBS 理论，可将复杂的系统用 MBS 拓扑结构描述。MBS 拓扑结构的描述方法，是 MBS 理论的基本问题。对于任意 MBS，可用低序体阵列对系统拓扑结构进行数字化描述。29肿瘤靶向治疗中的建模及仿真研究图 4.1 MBS 拓扑结构Fig. 4.1 MBS topology framework首先任选一物体为 B1,然后沿远离 B1 的方向，沿自然增长的数列依次标定每个物体的序号，从系统的一个分支到另一个分支，直到全部物体都标定完毕。如记惯性参考系 R 为 B0 体，则该系统的低序体阵列如表 4.1 所示，系统中每个物体的较低序号物体，用 L(K)表示，L 表示低序列体算子，K 表示物体的序号，它满足：Ln(K)=L(L(n-1)(K)(n,K 为正整数)，且补充定义 L0(K)=K,L(0)=0，则多体系统的任何一个物体都可以通过低序体阵列追溯到它与惯性系的关系。表 4.1 MBS 低序体阵列Table 4.1 MBS low-order body array低序体阵列 L(K)描述了开环多体系统的拓扑构造特点，在 L1(K)这一行中末出现的序号对应末端体，如 B3,B5,B7：在 Ln(K)(n1)中重复出现的序号为分支，如 B1 体，对于典型体 Bk，低序体阵列给出它及其所在分支的所有低序体序号。30河北工业大学硕士学位论文4-2-2 相邻体间坐标变换矩阵建立np3Opnp2Pknp1Sknj3nk3OjQkOknj1nj2Rjnk1Rknk2BjRno3RoBkOno2no1图 4.2 MBS 中的典型体与相邻低序体Fig 4.2 The Classic System and The Low-order System in MBS图 4.2 表示理想情况下，多体系统中的典型体及其相邻低序体。图中，Bk 为典型体，B j 为其相邻低序体。 n0 为惯性坐标系， n j 和 nk 分别为固联在 B j 和 Bk 体上的体坐标系(动坐标系)；O 为惯性系坐标原点， O j 、 Ok 为 B j 和 Bk 体坐标系原点； n p 为典型体 Bk 运动参考坐标系， O p 为典型体 Bk 运动参考坐标系坐标原点。 Pk 为典型体运动参考坐标系坐标原点在 B j 体坐标系中的位置矢量，Pk = Pkx Pky Pkz T ； sk 为典型体件坐标系相对其运动参考坐标系的位移矢量，sk = skx sky skz T ；Qk 为 Bk 体坐标系原点在 B j 体坐标系上的位置矢量，Qk = Qkx Qky Qkz T ； rk 、rj 、ro 分别是典型体 Bk 上给定点 R 在体坐标系 nk 、n j 及惯性坐标系 n0 中的矢径，rk 、rj 、ro 分别为rkx rky rkz T 、rjx rjy rjz T 、roxroy roz T 41。设 Qk = 0 ，且相邻体坐标系 nk 相对 n j 转过的相对方位角为a K 、 bK 、g K 。则典型体 Bk 及其相31肿瘤靶向治疗中的建模及仿真研究邻低序体 B j 体坐标系间的关系可用下式表示。n j = (n j nk1 )nk1 + (n j nk 2 )nk 2 + (n j nk 3 )(4.1)n j = SJK nk (4.2)其中SJK 称为相邻体坐标系之间方位变换矩阵。当a K 、bK 、g K 非常小时，SJK 可简化为下式：1 00 10bK 1- g K0 1- g KbKg KSJK = 01-aK 010g K10= 1-aK 0 aK1- bK01 001- bKaK1(4.3)当 Qk 不为零时，典型体上任意点 R 在其相邻体坐标系中的位置可用下式表示:R j = Qk + SJK Rk (4.4)上式用 4 4 阶 Denavit-Hartenberg 齐次变换矩阵为R R (4.5)j= AJK k114 4 阶矩阵AJK称为理想情况下，MBS 中典型体与其相邻低序体坐标间相邻体间变换矩阵。 SJK Q AJK = k(4.6)01式中， 3 3 矩阵SJK描述了 Bk 体坐标系 nk 相对于 B j 体坐标系 n j 的转动， 3 1向量Qk 描述了 Bk 体坐标系在 B j 体坐标系中的位置。设 AJK = AJK p AJK s SJK P SJK S = 0pk sk(4.7)1 01式中 AJK p SJK P = pk(4.8)01 AJK s SJK S = sk(4.9)01定义：相邻体间 4 4 变换矩阵 AJK 的各组成矩阵 AJK p 、 AJK s 为坐标系特征矩阵，简称特征矩阵，则 AJK p 、 AJK s 为理想状态下，相邻体间特征矩阵。32河北工业大学硕士学位论文其中， AJK p 为理想状态下， Bk 体运动参考坐标系相对其低序体坐标系的特征变换矩阵，由典型体运动参考坐标系在相邻低序体坐标系中的方位矩阵 SJK p 和位置矢量Pk 确定； AJK s 为理想状态下，Bk 体坐标系相对其运动参考坐标系的特征变换矩阵，由 Bk 体坐标系相对其运动参考坐标系的运动方位矩阵 SJK s 和移动矢量Sk 确定。因此，只要确定了相邻体间各坐标系特征矩阵，即可建立相邻体间变换矩阵。实际上，我们已经把建立 MBS 系统中相邻体间变换矩阵问题转换成了建立坐标系特征矩阵问题。下面讨论理想情况下，基于坐标系的特征矩阵形成规律。当典型体 Bk 体运动参考坐标系相对于 B j 体坐标系的位置和方位改变时，其主特征为特征矩阵1-g kpbkppkx AJK g1-ap = kpa1kpky (4.10)p - bp kpkpkz 0001 式中，akp 、 bkp 、g kp 为 Bk 体运动参考坐标系相对 B j 体坐标系转过的方位角； pkx 、 pky 、 pkz 为典型体运动参考坐标系原点在 B j 体坐标系上位置矢量的量化。当典型体 Bk 运动参考坐标系 np 与 B j 体坐标系 n j 平行时，akp = bkp = g kp =0，转化为 100 AJK = 010p 00100 0pkx p ky (4.11)p kz 1 当典型体 Bk 运动参考坐标系 np 相对 B j 体坐标系 n j 只有方位变化时， pkx = pky = pkz =0，转化为1-g kpbkp0 AJK g1-a0= kpa1kp(4.12)p - bkpkp00001若典型体 Bk 运动参考坐标系与 B j 体坐标系重合， akp = bkp = g kp =0 ， pkx = pky = pkz =0，则 AJK p =I，此时典型体体坐标系与其低序体体坐标系重合。理想状态下，相邻体间运动坐标系特征变换矩阵 AJK s 可按下列公式确定。当典型体相对其相邻低序体运动坐标系方向特征为 X（或 Y、Z）时，33肿瘤靶向治疗中的建模及仿真研究 100 x AJK ( x) = 0100(4.13) 0010001 0 1000 AJK ( y) = 010y(4.14) 0010001 0 1000 AJK ( z) = 0100(4.15) 001z 00 01当典型体相对其低序体的运动特征为 R，方向为绕 X 轴转动，转动角度为 （或绕 Y 轴转动 、Z轴转动 ）时，相邻体间理想运动特征变换矩阵 1000 0cos(a )-sin(a)0 AJK (a) = sin(a )cos(a)(4.16) 0000 01 cos( b )0-sin( b )0 AJK ( b ) =0100(4.17) sin( b )0cos( b )00001 cos(g )sin(g ) 0 0 - sin(g )cos(g )00 AJK (g ) = 001(4.18)000014-2-3 典型体上给定点理想位置方程理想状态下，典型体上给定点 R 在惯性坐标系中的位置矢量 R0 可用下式表示。R R (4.19)0= AOK k11式中， AOK 为理想状态下，典型体 Bk 相对于惯性坐标系 R 的总变换矩阵。 AOK 可用下式表示:34河北工业大学硕士学位论文0 AOK = AJK (4.20)u典型体上给定点 R 在惯性坐标系中的位置矢量 R0 可用相邻体间变换矩阵表示:R R0 R(4.21)0= AOK k= AJK k11u1从上式可以看出，典型体给定点可通过相邻体间变换矩阵按低序体阵列顺序连乘，追溯到惯性系 R。4-3 基于 MBS 理论的氩氦刀精确定位建模方法研究4-3-1 坐标系转换关系氩氦刀靶向肿瘤治疗系统是典型的 MBS。运用 MBS 基本思想解决氩氦刀肿瘤靶向治疗中的精确定位问题，可以进一步提高氩氦刀肿瘤靶向治疗的定位精度。本节根据前面对 MBS 基本理论建模方法的研究，探讨氩氦刀肿瘤靶向治疗中的精确定位建模方法。我们可以建立多体系统（MBS）相邻体间坐标变换矩阵。如图 4.3 所示建立三个坐标系，A 系相对于基准坐标系的旋转矩阵为 M A ，平移矢量为 S A ；B 系相对于基准坐标系的旋转矩阵及平移矢量分别为 M B 和 S B 。若在 B 系中有一点 PB 和一方向矢量 N B ，求该点或该矢量在 A 系中的表示。由于 A 系和 B 系之间没有直接的转换通路，因此需要借助他们的公用中介基准坐标系来建立二者之间的关系。首先求得 PB 和 N B 在基准坐标系 Or X r Yr 中的数学表达式，如下：P = M T P + SrBBB N r= M BT N B（4.22）然后用 A 系表达基准坐标系中的测点 P 和方向矢量，得出表达式，P = MA(P - S )ArAN A = M A N r（4.23）YrYBNBPBPOBYAXAPAXBSBSAOAOXr35肿瘤靶向治疗中的建模及仿真研究图 4.3 两个坐标系之间的数据转换Fig. 4.3 Data Transformation of Two Coordinate Systems将式（4.23）代入式（4.22）可得如下公式，实现两坐标系之间的数据转换。P = MA(M TP + SB- S )ABBA N A = M A M BT N B（4.24）这样，B 系中一点 PB 和一方向矢量 N B 在 A 系中就可以表示出来了。对于本课题来说，其结构如图 4.4 所示，其中 O 为人体定位坐标系原点， O j 为氩氦刀坐标系原点，由氩氦刀结构参数确定， Ok 为靶区坐标系原点，是我们定义的靶区原点，我们将图像序列的第一张图像的左上角定义为靶区坐标系的原点。 rk 为靶区中某一点在靶区坐标系中的矢径，需要求得该点在氩氦刀坐标系当中的矢径 rj 。氩氦刀坐标系相对于人体定位坐标系的旋转矩阵为 M A ,平移矢量为S A ，靶区坐标系相对于人体定位坐标系的旋转矩阵为 M B ，平移矢量为 S B 。因此在靶区坐标系中某点的矢径 rk 在人体定位坐标系中的矢径为:r = M T r + S(4.25)oB kB而在氩氦刀坐标系当中rj = M A ( ro - SA )(4.26)因此可以求得:r = M( M T r + S- S)jAB kBA(4.27)36河北工业大学硕士学位论文nk3nk1nj3Oknk2氩氦刀坐标系nj2rkROjrjnj1靶区坐标系n03roOn02n01人体定位坐标系（床身坐标系）图 4.4 氩氦刀系统坐标系Fig. 4.4 Cryocare Systems Coordinate System4-3-2 坐标系的建立（1）靶区坐标系的建立靶区坐标系主要是指在病灶部分建立的坐标系，还包括重要器官及周围环境。通常是根据医生根据所获得的 CT 数据来确定坐标系和原点的。该坐标系的作用就是医生可以根据建立的靶区坐标系，来确定靶区的位置，以达到较好的治疗效果。但是在常规治疗中，该坐标系一般只是在二维图像中建立 ，三维坐标只能抽象的建立在医生的头脑当中，这样就会造成治疗当中的定位不精确，影响治疗的效果。而在计算机辅助精确定位系统中，我们可以根据医学图像序列建立三维空间坐标，使其当中的每个像素点都有精确的空间位置。由于图像信息中包含了图像像素间距和实际空间距离比以及断层图像之间的层厚，可以经过简单处理得到与像素像素点一一对应的像素空间坐标信息，我们可以以序列图像的第一张左上角第一个像素为坐标原点，以向右为 X 轴正向，向下为 Y 轴正向，根据右手定则确定 Z 轴方向。我们可以利用图像处理技术提取出靶区的边缘，这样靶区边缘点的空间位置坐标就都确定了。（2）氩氦刀坐标系的建立我们可以以氩氦刀探针的顶点作为原点，建立氩氦刀坐标系。当计算机控制氩氦刀进行治疗时，探针的位置是可以不断改变的，其坐标应该可以通过计算机数据查找出来。（3）人体定位坐标系的建立以上建立的两组坐标系是通过计算机建立的虚拟坐标系，在治疗过程中均处于人体内部，并没有也不可能实际标记出来。因此，我们要想实现准确的定位就必须还有一组基准坐标系也就是人体的坐标系以供参考。目前在实际治疗中，人体定位的方法有很多种，一般较常用的有体内或体外标记点定位、通过人体自身的结构定位以及设置传感器等定位方法。这些定位方法的特点就是可以直接得到人体定位的坐标系，并且该坐标系是可见的，可以直接应用于治疗过程中。37肿瘤靶向治疗中的建模及仿真研究4-3-3 定位模型的建立在 CT 片中，我们很容易知道两个进刀点在病体坐标系中的坐标，利用病体坐标系与基准坐标体之间的结构参数进行坐标转换的计算，就可以得到这两点在基准坐标系中的坐标，又已知基准坐标系与氩氦刀坐标系之间的结构参数，在进行坐标转换，这样，就可以得出这两点在氩氦刀坐标系下的坐标，从而可控制氩氦刀针尖运动到这两点，进行肿瘤消融治疗。坐标系的变换包括坐标系的平移和旋转，如图所示建立两个坐标系。XYZ和 XYZ 分别为新旧坐标系，并且已知 OX、OY和 OZ三个矢量在旧坐标系中的方向余弦如下：ZZYOOYXX图 4.5 坐标系的平移与旋转Fig. 4.5 The Translation and Rotation of the Coordinate SystemO X = l1 , m1 , n1OY = l2 , m2 , n2 O Z = l3 , m3 , n3 (4.28)新坐标系任意一点 P(x, y, z) 与旧坐标系中的任意一点 P(x, y, z) 的关系用矩阵和矢量形式表示可简化为P = M (P - S)（4.29）其中，P = x yz xxl01P = yS = y0 M = l2，z zl ，0 ，3mn11m2n2 m3n3令 CT 片中标记点 Op 为坐标原点，沿以上坐标轴方向建立病体坐标系，并已知在该坐标系下 Pa 点坐标为（Xa,Ya,Za）, Pb 点坐标为（Xb,Yb,Zb）,又已知病体坐标系的三个坐标轴方向的矢量 OpXp,OpYp 和 OpZp 在基准坐标系下的方向余弦为：38河北工业大学硕士学位论文O p X p = l1, m1, n1O p Y p = l 2, m 2, n2（4.30）O p Z p = l 3, m3, n3且病体坐标系原点 Op 在基准坐标系下的坐标为（x0,y0,z0）,该坐标可根据测量得出，那么 Pa 点在基准坐标系下的坐标（x1,y1,z1）为： x1l1 m1 n1Xa x0 y1= l2 m2 n2（-1）* Ya + y0(4.31) z1l3 m3 n3 Za z0又已知氩氦刀坐标系原点 Od 在基准坐标系下坐标为（Xd,Yd,Zd），该坐标也可根据测量得出，氩氦刀坐标系三个坐标轴方向矢量在基准坐标系下的方向余弦为：Od X d = l1 , m1 , n1O Y= l 2, m 2 , n2(4.32)d dO Zd= l 3 , m3 , n3d那么基准坐标系下一点 Q(Xq,Yq,Zq)在氩氦刀坐标系中的坐标 Q(Xq,Yq,Zq)为：Xq l1m1n1 Xq - Xd m2n2（4.33） Yq = l2* Yq -Ydm3n3 Zq l3 Zq - Zd综上所述，设病体坐标系 Pa 点在氩氦刀坐标系中坐标为（x,y,z），且病体坐标系原点 Op 在基准坐标系下坐标为（x0,y0,z0）, 氩氦刀坐标系原点 Od 在基准坐标系下坐标为（Xd,Yd,Zd），则： x l1 m1 n1 l1m1n1m2n2( -1) y = l2 m2 n2 l2m3n3 z l3 m3 n3 l3 Xa * Ya + Za x 0 y 0 z 0Xd （4.34）- Yd Zd 同理，设病体坐标系下 Pb 点在氩氦刀坐标系中坐标为（x,y,z）,其余参数同上，则 x l1 m1 n1 l1m1n1n2( -1) y= l2 m2 n2 l2 m2m3n3 z l3 m3 n3 l3 Xa * Ya + Za x 0 y 0 z 0Xd （4.35）- Yd Zd 这样，病体坐标系中由 Pa,Pb 两点确定的矢量 (Xb-Xa,Yb-Ya,Zb-Za)T，在氩氦刀坐标系下对应的矢量为(x-x,y-y,z-z) T，则 x - x l1 m1 n1 l1m1n1 Xb - Xa m2n2( -1)*（4.36） y - y = l2 m2 n2* l2 Yb -Ya m3n3 z - z l3 m3 n3 l3 Zb - Za 4-4 本章小结本章分析了基于多体理论的建模方法和特点，提出了基于多体理论建模的规律性，并利用 MBS 理论建立了氩氦刀肿瘤靶向治疗的精确定位模型。主要内容如下：39肿瘤靶向治疗中的建模及仿真研究（1）首先介绍了多体系统的概念及如何建立多体系统的拓扑结构和低序体阵列。（2）详细阐述了基于多体系统理论相邻体间的齐次变换矩阵特征变换矩阵。（3）针对氩氦刀肿瘤靶向治疗系统建立坐标系，并建立它们之间的特征矩阵，根据特征矩阵建立定位模型。40河北工业大学硕士学位论文第五章：基于 VRML 的氩氦刀肿瘤靶向治疗中的仿真研究5-1 引言本章通过对氩氦刀肿瘤靶向治疗的研究，在前面几章建立了靶区模型、氩氦刀冰球模型和定位模型的基础上，对氩氦刀冰球模型和手术消融过程进行仿真。氩氦刀肿瘤靶向治疗手术仿真使医生可以事先发现手术中的问题，并帮助医生合理、定量的制定手术方案，选择最佳的手术路径，减少对周围健康组织的破坏，还可以实现对医生的手术培训。这对提高定位精度、开展复杂的外科手术和提高手术成功率具有重要意义42。目前，虚拟手术仿真正