毕业论文- 基于Elastix的CT与MRI图像配准研究

湖南大学毕业设计(论文)HUNANUNIVERSITY毕业论文论文题目基于Elastix的CT与MRI医学图像配准研究学生姓名学生学号专业班级自动化2班学院名称电气与信息工程学院指导老师学院院长2015年 月日 湖南大学毕业设计(论文)第II页 湖南大学毕业设计(论文)第II页第I页 湖南大学毕业设计(论文)第I页摘要医学图像是现代医疗诊断的重要方法，由于不同的物理理论不断地应用到医学影像中，医学成像设备不断升级，出现了各种各样的成像类型，常见的成像手段有核磁共振成像（MRI）、正电子发射型计算机断层显像（PET）、计算机断层成像（CT）等等，各类成像方式因为有成像原理上的差异，它们提供的图像具有不同的医学特征信息，且由于病人在不同时间，不同条件下，或是不同的设备（或是同一设备的不同参数）产生的图像也会有所差异，因此各图像的医学特征信息可能呈现片面化，但同时各图像的信息也是互补的，通过图像的配准可以全面地了解病人病情，方便医生进行诊断。本文主要通过在Elastix中进行实验的方式，探究了头颈部肿瘤CT与MRI图像的刚性配准以及弹性配准的性能特点，在此之外还进行了联合测度配准的实验，经过实验证明如何选用算法和参数才能得到较好的配准结果。关键词：图像配准；Elastix；ITK；刚性配准；弹性配准；联合参数配准CT/PET/MRIimageregistrationbasedonElastixAbstractMedicalimageisanimportantmethodofmodernmedicaldiagnosis,duetodifferentphysicaltheoriescontinuetobeappliedtomedicalimaging,andmedicalimagingdevicesescalating,appearedinavarietyofimagetypes,commonimagingmeanshavemagneticresonanceimaging(MRI),positronemissiontomography(PET),computedtomography(CT)andsoon,becausetherearedifferencesbetweenvarioustypesofimagingmodalityimagingprinciple,imagestheyofferdifferentfeaturesmedicalinformation,andbecausedifferentpatientstime,theimageunderdifferentconditions,ordifferentdevices(ordifferentparametersofthesamedevice)generatedwillvary,andthereforethemedicalinformationofeachimagefeaturemayappearone-sided,buteachimageinformationiscomplementarytobyimageregistrationcanfullyunderstandthepatient'scondition,convenientfordoctorstodiagnose.Inthispaper,byperformingexperimentsinElastixthewaytoexploretheperformancecharacteristicsofrigidregistrationandelasticregistration,andinadditiontothisalsocarriedoutjointmeasureswithquasiexperimentafterexperimenttodemonstratehowthechoiceofalgorithmsandparametersinordertogetbetterTheregistrationresults.KeyWords：imageregistration;Elastix;ITK;rigidregistration;deformableregistration;multi-metricregistration第I页 湖南大学毕业设计(论文) 第I页目录TOC\o"1-3"\h\u83261绪论 页1绪论1.1医学图像配准研究背景及其意义自从上世纪初，X射线被用来作为医学诊断手段以来，医学诊断技术一直在更新换代，新的造影技术层出不穷，不断成熟的医学成像技术使得临床诊断更加方便和直观，确保了医生对病人病情的正确判断，为后续的有效治疗提供了坚实的保证。医学图象技术所探究的内容包括图象的获取、重建、配准、分割和分析等方面，而此中的医学图像配准环节是一个非常实用且重要发研究方向领域[10]。伴随着造影技术的升级，出现了各种不同的造影技术，从最早的X射线成像，到如今的计算机断层成像（CT）、核磁共振（MRI）、电子发射型计算机断层显像（PET）和单光子发射计算机断层成像术（SPECT)等技术，因为原理的不同，出现了多种模态的医学图像，各模态图像为医生提供了不同侧面的人体医疗信息。各成象技术根据其所呈现的医学特征信息可分为两个类别：一类是描述人体的解剖特征的解剖结构图像，例如CT、MRI等；还有一类是可以显示人体的新陈代谢信息的图像，例如PET、SPECT等；两种成像手段各有其特点，各有其优劣，解剖结构图像可以获得各器官的解剖结构信息，图像分辨率高，但是对疾病敏感度不高。当检测到疾病时，组织结构已经发生病变，甚至已经到了疾病的中晚期。而功能结构图像可以获取各个器官的新陈代谢信息，能检测到细胞以及分子层面的病变信息，可以进行疾病的早期诊断，特别是对于癌症治疗，疾病的早期诊断是极其重要的。遗憾的是功能结构图像的分辨率低，难以分辨解剖结构。并且同一种类别的造影技术也各有侧重，例如CT图像相对于MRI图像来说侧重于骨组织；而MRI相对于CT图像来说侧重于软组织。由此可见，要想通过医学图像来了解病人的病情全貌，仅仅使用单模态的图像是远远不够的，于是催生了多种模态图像配准的思想。医学图像配准技术通过计算机算法找出不同模态图像的空间对应关系，通过一系列形变，使两幅图像(或者多幅图像）在空间上对齐，进而得到信息整合的综合图像。1.2国内外研究现状随着医疗设备的不断更新换代，计算机的存储容量不断扩张，处理器的速度不断提高，以及图像配准算法的不断改进，图像配准在日新月异地发展中。国外的图像配准技术起步早，对其研究已形成比较完整的体系，各大研究机构以及高校与企业开发了许多图像配准软件和算法平台。当今距离单模态匹配问题的基本完成已经过去了几十年，但对于多模态匹配，因其问题的复杂性，到现在仍是图像处理中备受关注的焦点问题[1]。二十世纪八十年代以来，控制理论的不断发展，以及智能控制理论的出现，为图像处理提供了更多的途径与方法。例如细胞神经网络（CNN）如今已经被广泛应用于图像处理，在图像处理中的线检验、图像分割、图像编码、边缘方向检测、孔填充、运动检验、特征提取、字符识别、去噪等等领域都取得了一定的成果[2]。国内的图像配准技术起步较晚，各研究机构与高校相继开展了相关研究，提出了一些新的可行算法，对图像配准的发展起到了重要的作用。1.3配准研究方法现存的配准算法有很多，不同的算法有各自的优劣，不同的算法适用的场合也不尽相同，于是需要通过对比优化实验，得到最佳的配准算法以及配准参数。通过使用Elastix配准软件，可以快速使用ITK上几乎所有配准算法，且无需进行复杂繁琐的程序编写，只需编写模块化的参数文档，再在CMD.EXE中调用ELASTX参数文档即可完成配准实验，大大提高了配准实验的效率，以便让研究人员把更多精力集中在算法的优化和对比中。模块化的参数文档配置方便研究人员快速尝试各种配准算法，调试各种配准参数，以得到满意的，结合自己研究实际的算法和配准参数文档。本文致力于研究医学图像的刚性和弹性配准，及其联合测度配准。通过设计配准实验，对比刚性配准和弹性配准各自的特点，总结刚性与弹性配准各自所适用的配准对象，通过对结果的配准速度、配准精度等指标进行评估（可视化评估，以及量化评估），最终得到我们需要的，能用于病人实例的配准参数文档。2图像配准基本理论想要完成配准的实验，应当先掌握其基础理论知识，以及对相关软件的熟练运用。本章作为后续章节的铺垫，叙述图像配准中的基本原理以及Elastix的基本使用知识。2.1图像配准的定义在配准的过程中，其中一幅图像是需要进行变换的图像，被称为浮动图象（movingimage），另一张作为基准的无需变换的图像被称为固定图像（fixedimage）。设浮动图像和固定图像的维度为d，在他们各自的空间域中分别有：以及。且设一种变换，使得在空间上与对齐，这就是配准的最终目的[5]。其示意图如下：图2.1图像配准的任务是寻找两幅图像的空间映射关系，左为固定图像，右为浮动图像从数学实质上来说，医学图像配准的过程是一个寻优过程，但是在寻优过程中有很多参数，是一个很复杂的优化过程，这个优化过程的代价函数可如下表示：，（2-1）配准过程可以表示为存在关于参数μ的空间变换使得两图像配准时代价函数最小。2.2医学图像配准的一般流程一般来说，图像配准有四个基本的步骤：2.2.1特征提取这一步进行图像特征的提取，通常有三种特征是经常被使用的[12]：（1）特征点：即为在几何上有有意义的特殊点。特征点又可分为内部特征点以及外部特征点，人体内部固有的解剖结构点被称为内部特征点，而在外部人为标记而成的特征点被称为外部特征点，由于外部特征点需要人为干涉，并且对患者造成痛苦，已较少使用，现在医学图像处理主要基于内部特征点；（2）特征曲线（曲面）：用图像分割的方法，把需要进行处理的部分提取出来，以便后续处理；（3）灰度统计特征，灰度统计特征由于其无需对图像进行预处理，是图像处理自动化的一个重点研究工具，本文所用的互信息测度准则就是基于统计特征来进行计算的。2.2.2空间变换经过特征提取后，图像配准问题就变成了浮动图像的特征空间到固定图像的特征空间的空间变换问题，也即找出浮动图像特征空间与固定图像特征空间的映射关系。凭借各类空间变换性质的不同可以把空间变换分为两个大类：刚性变换（rigidtransformation)与弹性变换(deformabletransformation)，刚性变换允许对图像进行缩放、旋转、位移等操作，但是不能扭曲图像曲线，刚性变换一般用来对图像进行粗配准，把图像的位置对齐，但是无法处理形变区域进行；非线性变换可以处理形变区域，一般用于刚性变换之后进行局部处理，尤其是对浮动图像畸变区域（如果存在的话）进行矫正。2.2.3相似性测度用来评估两图像的匹配程度的度量函数被称为相似性测度，它用来测量图像配准是否达标。相似性测度会明显影响到配准的精确度，是配准过程中的重中之重。2.2.4优化策略选择合适的优化算法在配准中非常重要，因为医学图像配准的过程本质上是一个多参数的寻优过程，选择合适的优化算法，可以迅速找到最优解，而不必遍历整个特征空间。以上的几个步骤可以表示为下图：相似性测度相似性测度优化策略配准结果浮动图像图像插值空间变换固定图像像素 优化结束像素变换参数图2.2配准流程简单示意图2.3Elastix中的配准Elastix是基于著名的图像分割与配准工具包ITK的开源软件，Elastix集成了很多用来解决图像配准问题的常用算法。Elastix使用非常简单，用户只需编写简单的模块化参数文件加以调用就能快速进行配准实验，免去了繁琐的程序编写与调试，让用户集中更多精力在配准算法的对比与优化上。下面将介绍Elastix陪准过程中的各大模块，本文所用的完整的参数文档请参见附录，参数文档中各参数的推荐设置可在http://Elastix.isi.uu.nl/doxygen/parameter.html中查询。2.3.1相似性测度（Metrics）在上节所述，相似性度量是一个通用的配准过程中十分重要的一步，它决定了图像配准的精度和准确度，与配准结果的质量直接相关。Elastix中提供了几种常用的测度：平方和的差（SSD），在Elastix中选用它的语句为(Metric"AdvancedMeanSquares")，SSD的定义如下：，(2-2）其中为固定图像的值域，是固定图像的像素数量。根据已给的变换，SSD能在固定图像中轻松迭代实现。互信息（MI），在Elastix中通过语句（Metric”AdvancedMattesMutualInformation")来选择互信息用于相似性测度计算,其定义如下：，（2-3）其中和表示两图像上某相同位置的灰度值，是联合概率分布，以及是边缘概率分布。互信息也可以定义为：，（2-4）其中和分别为固定图像和浮动图像的熵，为两图像的联合熵，当两图像完全相同时，联合熵为0，此时互信息最大。互信息表示两幅图像的统计相关性，也即一幅图像包含另一幅图像的信息多少。在Elastix中使用的是Mattes等人提出的负的互信息，其值越小，表示两图像相关程度越高[5]。归一化互信息（NMI），在Elastix中用（Metric”NormalizedMutualInformation”)来选择使用NMI作为相似性测度。 归一化互信息可定义为：，（2-5）其中H为熵。对比互信息的熵定义（2-4）：，可以发现归一化互信息是为了避免当两张图像重叠区域太少，导致互信息太小而配准失效的情况，基于归一化互信息的图像配准具有更好的鲁棒性和配准精度。2.3.2图像采样器（Imagesamplers)通常来说，把图像所有像素进行采样可以获得最好的配准精度，但在实际应用中，这是没有必要的，很多时候只需要采样一定的子集即可满足需求。选择什么样的采样方法需要用户根据自己的实际需求进行选择，在参数文档中用语句(ImageSampler"XXX")选择采样器，其中XXX为采样器名称，Elastix中提供了如下几种常用的采样方式以供选择：全采样（Full），即如上述所说的对所有像素进行采样的方法，Elastix中虽然提供了这种方法，但是因为效率过低，且没有必要，通常不使用全采样。网格采样（Grid），网格采样可理解为在固定图像上加上一定密度的网格，随之在网格点上进行采样。相对于全采样，网格采样大大提高了采样速度。用户可以根据自己的实际需要在参数文档中设置网格的大小。随机采样（Random），随机采样会从固定图像中随机地采样出用户设置的像素值。被采样图像中的每一个像素被采样的概率相同，并且可能被多次选中，也可通过参数文档设置避免同一像素被反复采样，采样数量由用户设置，是一种常用的采样方式。随机坐标采样（RandomCoordinate），随机坐标采样和随机采样类似，但随机坐标采样不止限于像素点的采样，两个相邻的像素点之间的位置也可能被选中，和随机采样一样，随机坐标采样也是一种常用的采样方式，本文实验部分采用的就是随机坐标采样。2.3.3插值器(Interpolators)数字图像是由离散的像素组成，在对图像进行处理的过程中，变换后的图像像素坐标可能不是整数，为了使变换后的图像重新落在整数坐标轴上，形成新的数字图像，必须对变换后的图像进行合适的插值操作。Elastix中同样提供了几种常用的插值器供用户选择使用，参数文档中选择插值器的语句为(Interpolator"XXX")，其中XXX为插值器在Elastix中的名称，下面列举几种Elastix中的插值器：最近邻插值（NN），在Elastix中的名称为NearestNeighborInterpolator，最近邻插值是一种最简单也最容易理解的插值方式，即将变换后图像的像素点赋予与它最近的坐标点的值。这种算法的优点是实现简单，计算量小，计算速度快，但是往往精度不高。线性插值（Linear），在Elastix中的名称为LinearInterpolator，也被称为一阶插值，该法用非整数像素点附近的四个点的加权平均值（各点的权重为他们各自与非整数像素点的距离）来得出此位置的灰度值。高阶插值（N-thorderB-spline），也称为N阶B样条插值，在Elastix中的名称为BSplineInterpolator。最近邻插值和线性插值可以分别看成阶数N=0和N=1的高阶插值，通常来说，高阶插值的阶数越高，所得到的图像质量越好，但同时也会消耗更多的计算资源，并且到一定阶数以上之后，得到的图像效果改善不会太大，下图会说明这一现象。图2.3为几种不同的插值器所得结果：图2.3不同插值器结果，（a）为最近邻插值，（b）为线性插值，（c）为2阶B样条插值，（d）为3阶B样条插值，（e）为5阶B样条插值。由上图可见，二阶以上的插值器效果改善不很明显，通常一阶插值就能满足实际需要，并且一阶插值有很好的性价比，所以我们经常在配准过程中使用一阶插值。而在产生最终图像时，特别是最终图像需要进行弹性形变时，通常使用三阶B样条插值完成，以得到高精度的输出图像。产生最终图像的插值器与配准过程中的插值器在语句上有一点区别，最终插值器的语句为(ResampleInterpolator"FinalXXX")，其中XXX为插值器名称，例如(ResampleInterpolator"FinalBSplineInterpolator")2.3.4几何变换（Transfomation）几何变换的方向关于几何变换的方向，经常会有误解。直觉上来看，感觉是浮动图像进行几何变换后，映射到固定图像的特征空间，再与固定图像进行相似性测度计算。而实际上，在Elastix中，几何变换的定义为从固定图像空间映射到浮动图像空间[7]，即，（2-6）下图说明了这一过程：插值采样插值采样图2.4几何映射的方向由图可见，实际中的几何变换，是固定图像接受几何变换后，经过相似性测度的评价，得到浮动图像中的映射坐标，最后浮动图像根据被映射的点，经过插值后得到配准结果。Elastix这样定义变换方向的原因是，配准过程中，固定图像需要进行采样后得到物理坐标系下的图像，再从固定图像的物理坐标系映射到浮动图像的物理物理坐标系（此时需要进行几何变换），被映射得到的浮动再经过插值运算反映射到逻辑坐标系中。如果把几何变换定义为从浮动图像映射到固定图像，则重采样过程中不能保证固定逻辑坐标系图像采样后得到唯一的值，这会导致采样后的图像出现空洞，所以才把几何变换的方向定义为从固定图像映射到浮动图像，以保证与重采样方向的一致性。Elastix中的常用几何变换Elastix中同样集成了几种常用的几何变换以供用户根据实际需求选择使用，在参数文档中选择几何变换的语句为(Transform"XXX")，其中XXX为几何变换在Elastix中的名称，例如(Transform"BSplineTransform")，表示选择使用B样条变换作为配准过程中的几何变换。下面介绍几种在Elastix中的几何变换：（1）刚体变换（Rigid），在Elastix中的名称为EulerTransform，刚体变换可定义为，（2-7）其中R为旋转矩阵，t为位移量，c为旋转中心。在刚体变换中，整个图像被看作一个刚体，可以对其进行平移和旋转，但是不能对其进行缩放。在对同一病人的脑部进行配准时常使用刚体变换，因为同一病人的脑部图像形变不会太大，刚体变换能够满足配准要求。刚体变换也常常用来对图像进行粗配准，在复杂的形变较多的医学图像配准中，可以先用刚体变换，将图像全局对齐，之后再用弹性变换对局部形变进行处理，后一章将要进行的弹性配准以及联合参数就是用刚性配准的结果作为初始化变换来处理病人头颈部图像的配准问题。（2）相似变换（Similarity），在Elastix中的语句为SimilarityTransform，它的定义式如下：，（2-8）相似变换是在刚体变换的基础上加上一个伸缩因子s，它不仅可以进行刚体变换的所有操作，也能进行图像的放大缩小，变换前后的图像是互为相似的。（3）仿射变换（Affine），在Elastix中的语句为AffineTransform，其定义式为，（2-9）其中A为任意矩阵，仿射变换的含义是指变换前空间内两条平行的直线，在变换后的空间内任然保持平行,但是两条直线的夹角变换后角度可能会变化，变化后图像不与变换前图像相似。（4）B样条变换（B-splines），其在Elastix中的名称为BSplineTransform，B样条变换是一种弹性（非刚性）变换，它的定义式如下：，（2-10）其中表示操作点，是系数向量，是B样条操作点空间，是B样条在上的压缩支撑集，操作点被定义在固定图像上的网格上，是三次B样条多项式，其定义式如下：，（2-11）弹性变换会使直线经过变换后变成曲线，这种特性用来对图像中的形变进行矫正是非常有效的，但是需要注意设置多项参数，以免出现不切实际的变形。2.3.5优化算法（Optimiser）由于图像配准问题在本质上是一个多参数寻优问题，所以优化算法的选择至关重要，它直接关系到配准过程的速度以及算法的收敛性，Elastix中同样集成了几种常用的优化算法供用户选择使用，参数文档中选择优化算法的语句为(Optimizer"XXX")，其中XXX为优化算法在Elastix中的名称。下面介绍Elastix中的常用优化算法：梯度下降法（Gradientdescent，GD），在Elastix中的名称为StandardGradientDescent或者RegularStepGradientDescent，梯度下降法每一步的搜索方向是代价函数的负梯度方向，用式子可以表示为：，（2-12）其中，为k位置的梯度。梯度下降法可以处理大部分问题，在实验中使用的是梯度下降法的一个变种，自适应随机梯度下降法，在Elastix中的语句为AdaptiveStochasticGradientDescent，经过实验证明，自适应随机梯度下降法性能良好。进化算法(Evolutionaryalgorithms，EAs），在Elastix中的语句为CMAEvolutionStrategy，进化算法是一种全局优化算法，它的灵感来自自然界中生物的演化，在自然界中，生物总能通过遗传、变异和选择淘汰来与环境适应，进化算法通过选择、重组和变异这三个环节来模仿进化过程，实现优化过程的求解。进化算法是一种具有高自适应性，自组织性，以及高鲁棒性的全局算法。2.3.6多级分辨率（Multi-Resolution）很多医学图像，特别是三维图像的体积是十分巨大的，再加上配准过程中需要对图像进行计算量非常大的相似性测度以及几何变换的迭代计算，使得整个配准过程计算量非常巨大。为了减少计算代价，在Elastix中对图像进行多级分辨率处理，把整个配准过程分为几个子过程，形成图像分辨率金字塔，把分辨率分为从粗到细的几个层次。依次将上一层的配准结果作为下一层的初始条件进行处理，直到达到最高分辨率。下图为多分辨率配准的示意图：图2.5多级分辨率分层处理多分辨率图像配准与直接配准相比，有很多优势：配准的过程是从粗到细处理，每一层都有上一层的配准结果作为初始条件（最顶层是零初始条件），使得整个配准大部分过程都在低分辨率下进行，使得计算量大大减小；多级分辨率在处理过程中逐步包含全局信息，可以避免落入局部最优解;可以利用多级分辨率实现不同分辨率图像的配准，比如CT图像和PET图像的配准，前者分辨率高而后者分辨率低，可以把CT图像为低分辨率图像后再与PET图像配准。在Elastix的参数文档中使用多级分辨率的语句为(FixedImagePyramid"FixedSmoothingImagePyramid")和(MovingImagePyramid"MovingSmoothingImagePyramid")，前者为固定图像设置多级分辨率，后者为浮动图像设定多级分辨率，除此以外还需要为图像金字塔设置层级数量以及每一级的分辨率。例如后面实验会用到的，三维图像的多层分辨率可以如下配置：(NumberOfResolutions3)//分辨率有三个层级；(ImagePyramidSchedule882441111)//分辨率的分配依次为在x和y方向上为原始分辨率的1/8、1/4，最后一层为原始分辨率；在z方向上为原始分辨率的1/2、后面两层均为原始分辨率。2.4配准评估配准完成后，配准结果是否合格，是否能满足临床需求，还需要对其进行评判，一般来说，这是一件很难的事情，我们无法知道哪个像素应该映射到哪一点，但是任有一些配准评估的方法可以供我们使用。首先，可以用可视化工具，把配准后图像与基准图像（固定图像）进行叠加，从直观上观察配准是否合格，常用的医学可视化软件有MeVisLab、ITK-SNAP、ParaView、3DSlicer、VV等，其中除了MeVisLab以外都是开源软件，配准用户可以方便地使用他们快速观察自己的配准结果是否在直观上可行。但是仅仅以目测做为评判标准是远远不够的，还有一个好的方法是观察计算过程中的相似性测度，例如将在后面实验中所用的Mattes互信息，观察测度在配准过程中的曲线，若曲线是收敛的，说明配准是比较成功的，若曲线是发散的，说明配准参数还需调整。3基于Elastix的CT/MRI配准本章论述几组CT/MRI医学图像配准实验研究内容。通过实验主要进行刚性与弹性配准的性能对比，以及单相似性测度配准与联合相似性测度配准的性能对比。本文所述实验均在64位WINDOWS7系统下运行，Elastix版本为最新的4.7，编译器为VC++2010，Elastix的编译和安装使用跨平台的编译软件Cmake来完成。本实验的病人图像从医院获得，有两位病人的CT/MRI图像，对两位病人均进行同样的参数文档实验，以排除参数文档对于单一病人的特殊适用性，保证了配准参数文档一定程度的实用性。本实验的参数设置参考自Elastix官网的参数数据库（http://Elastix.bigr.nl/wiki/index.php/Parameter_file_database），并根据自己的实际需求进行了优化和改进。在Elastix中进行配准实验非常简单，只要编写好参数文档，并在DOS中输入调用语句即可进行实验，在DOS中调用Elastix的语句为Elastix-ffixedimage.ext-mmovingimage.ext-outoutputdirect-pparameterfile.txt，其中-f为固定图像的设置语句，-m为浮动图像的设置语句，-out为输出文件夹的设置语句，-p为本次配准所使用的参数文档的设置语句。还有一些可选的额外语句，在需要的时候加以调用：如-fMaskfixedmaskimage.ext-mMaskmovingmaskimage.ext-t0initialtransformparameterfile.txt，其中-fMask为设置固定图像掩膜的语句，-mMask为设置浮动图像掩膜的语句，-t0为设置初始变换的语句。Elastix可以支持大部分医学图像格式，例如.mhd、.hdr、.nii、.gipl、.dcm，但是DICOM文件夹不能直接在Elastix中使用，本文实验所用的病人图像是DICOM文件夹格式，为了Elastix方便处理，用MRIConvert软件将其转换成了.mhd格式。配准结束后用VV进行直观评估，从视觉上判断图像配准是否成功，同时用EXCEL画出相似性测度在配准中的变化曲线，对其配准过程的相似性测度进行定量分析，以此评判配准是否成功。3.1刚性及弹性配准实验下面将进行刚性与弹性配准的对比实验，从实验的角度来验证刚性配准和弹性配准各自的特点以及适用范围。因为当今医院放射科很多都使用了PET\CT一体机，CT和PET的配准已经在成像设备内部完成，所以本文主要完成的是CT与MRI图像的配准，CT图像对人体骨骼组织的分辨率较高而对软组织的成像效果差，而MRI图像对人体软组织的成像效果好而对软组织的成像效果差，因此对CT\MRI的配准是非常有实用意义的。3.1.1刚性配准实验（1）对第一位病人的刚性配准首先对第一个病人的CT\MRI图像进行刚性配准，500mm×500mm×392.4mm，MRI图像大小为240mm×240mm×216mm。以下是未配准前的CT\MRI图像叠加的图像：图3.1第一位病人的CT\MRI配准前叠加图像图中绿色的部分为MRI图像，紫色的部分为CT图像，可见配准前两图像在空间上有很大的错位，且如前面所介绍，CT图像包含了更多的骨组织信息，而MRI图像包含了更多的软组织信息，由于是头颈部部图像，骨质成分很多，形变较少，因此可以考虑使用刚性变换进行配准。配准过程以多级分辨率形式进行，分为三层分辨率层级，采样器选用随机采样，插值器选用B样条插值（过程中插值为1阶，最终插值为3阶），几何变换为欧拉变换，优化器选用自适应随机梯度下降法，相似性测度采用Mattes互信息，主要的参数设置如下：(Registration"MultiResolutionRegistration")(NumberOfResolutions3)//多分辨率配准，分为三个级层。(Interpolator"BSplineInterpolator")(ResampleInterpolator"FinalBSplineInterpolator")(BSplineInterpolationOrder1)(FinalBSplineInterpolationOrder3)//使用1阶B样条插值器，以及3阶B样条作为最终插值。(Optimizer"AdaptiveStochasticGradientDescent")(MaximumNumberOfIterations500)//优化策略为自适应随机梯度下降法，每个分辨率层级最大迭代次数500次(Transform"EulerTransform")//几何变换采用欧拉变换。(Metric"AdvancedMattesMutualInformation")//相似性测度选用Mattes互信息。完整的参数文档请参见附录A。整个刚性配准程序运行耗时2分40秒，以下为配准后的图像与固定图像（CT图像）叠加的图像：图3.2第一位病人的CT\MRI图像进行配准后的图像图中绿色部分为MRI图像，紫色部分为CT图像，可见配准后的图像在直观上已经对齐了，软组织信息和骨组织信息有效地形成了互补。同时对配准过程中的互信息进行了定量分析，可见于以下图表：图3.3对第一位病人进行刚性配准过程中的互信息图中横坐标为迭代次数，纵坐标为相似性测度，初始互信息为-0.64768，最终互信息为-0.601635。虽然从直观上来看两图像貌似是对齐了，但是从定量分析来看，互信息并没有很好地收敛，究其原因，是因为刚性配准对处理复杂图像的局限性。因为刚性配准只能进行位移和旋转运算，对图像的形变无法处理，实际上在第一位病人的图像中，虽然肉眼比较难分辨，但是放大后仍可见有很多明显形变的地方，这说明仅仅使用刚性配准无法完成符合临床需求的配准工作。对第二位病人的图像的刚性配准为了作对比，以及保证参数文档一定的泛用性，还用了另外一位病人的CT\MRI图像进行配准。第二位病人的CT图像大小为其配准前图像如下：图3.4第二位病人的配准前CT图像图3.5第二位病人配准前的MRI图像观察图像可知，第二位病人的CT图像出现严重的失真变形情况，无疑加大了配准的难度，若仅使用刚性配准，可能得不到理想的结果。使用与第一位病人相同的参数文档进行配准后的结果如下图：图3.6第二位病人刚性配准后的图像从图可见本次配准并不成功，究其原因是因为图像原本的形变比较大，而刚性配准无法进行形变操作，因为刚性配准无法完成本次配准，若想完成第二位病人的配准工作，必须在刚性配准的基础上，再使用弹性配准进行形变操作，下面将用实验证明是否可行。3.1.2弹性配准实验第一位病人的弹性配准实验首先对第一位病人进行弹性配准实验，在第一位病人的刚性配准中，已经几乎完成了配准工作，但是因为有微小形变仍然存在，还需要进行弹性配准进行矫正。弹性配准的参数配置如下：采样器选用随机坐标采样器，样本空间数设置为5000；几何变换使用1阶B样条变换；相似性测度使用Mattes互信息；使用三层图像金字塔；优化器选择自适应随机梯度下降法，每个分辨率层级迭代1000次；插值器与最终插值器分别采用1阶与3阶的B样条插值器；并将上小节的刚性配准结果作为弹性配准的初始化变换，主要的参数配置如下：(Registration"MultiResolutionRegistration")(NumberOfResolutions3)//多分辨率配准，原图像分为三个分辨率层级。(Interpolator"BSplineInterpolator")(ResampleInterpolator"FinalBSplineInterpolator")(BSplineInterpolationOrder1)(FinalBSplineInterpolationOrder3)//使用B样条插值，1阶过程插值，3阶最终插值。(Optimizer"AdaptiveStochasticGradientDescent")(MaximumNumberOfIterations1000)//优化策略选用自适应随机梯度下降法，最大迭代次数1000次。(Transform"BSplineTransform")//几何变换选用B样条变换。(Metric"AdvancedMattesMutualInformation")//测度选用Mattes互信息。实验结果如下：图3.7对第一位病人进行弹性配准的结果本次实验用时12分22秒，初始互信息为-0.599965，最终互信息为-0.665683。从图像上来看，与刚性配准的变化不大，但从互信息的角度可以发现弹性配准的性能较刚性配准确有改进，但同时也有不足的地方：在图像的一些地方出现了不切实际的形变，以俯视图为例，形变情况如下图：图3.8左图为MRI原图像，右图为弹性变换后的图像可见出现了比较严重的非实际形变现象，配准方法和参数任需改进。互信息的变化见下图：图3.8第一位病人弹性配准过程中的互信息从上图来看，互信息呈现一定的收敛趋势，说明弹性配准确实对配准的效果起到了提升作用。（2）第二位病人的弹性配准对第二位病人的刚性配准是失败的，因为第二位病人的图像存在大量形变，下面考察用弹性配准是否能提升配准质量。第二位病人的弹性配准使用与第一位病人相同的配准参数文档。其配准结果如下：图3.9第二位病人的弹性配准结果本次配准耗时14分25秒，初始互信息为-0.704621，最终互信息为-0.764674。从图可见，与刚性配准比较，弹性配准较好地完成了第二位病人的配准任务，但是仍有不少问题存在，比如出现了不切实际的形变，说明配准参数还需改进。对第二位进行弹性配准过程中的互信息变化过程如下图：图3.10第二位病人的弹性配准过程中的互信息变化曲线从图表中可见，互信息整体显示出下降趋势，但不如第一位病人的互信息下降明显，从互信息也可以发现配准不是很完美，还有改善的余地。3.2联合测度配准实验在Elastix中，不仅可以实现两幅图像使用一种相似性测度进行配准，更可以实现多测度联合使用，还能实现多张图像的联合配准，甚至可以进行十分复杂的多图像多测度联合配准实验，下面将进行联合测度配准的实验，以说明联合配准在配准精度上的优越性。在Elastix中，可以组合多种测度进行联合配准（一般是组合两种测度进行配准），它的好处是能提高配准的精度，但同时因为相似性测度是配准过程中计算量最大的环节，使用联合测度会不可避免地降低配准的速度。在使用联合测度时，代价函数变为如下形式：，（3-1）其中为第i个测度的权重，为第i个测度的类型。本文使用Mattes互信息+BEP（bendingenergypenalty）作为联合测度，其中BEP的定义式如下[15]：,(3-2)其中是像素位置，P是像素数量。是变换的第i要素，、、分别为、、方向的坐标。下面将对两位病人进行联合测度的配准实验。第一位病人的联合测度配准实验上面的实验中，对第一位病人的刚性和弹性实验中，虽然都得到了可以接受的结果，但是都存在问题，刚性配准无法处理形变，而弹性配准却产生了不切实际的形变。下面进行的实验着力于增加配准的精度，但是同时牺牲了配准的速度，在时间充裕的情况下，下述方法是十分可靠的。联合测度配准首先需要写入参数文档的语句是(Registration"MultiMetricMultiResolutionRegistration")，说明这是一次联合测度配准， 同时在选择测度时可以填写多个测度，其中用空格分开，例如(Metric"AdvancedMattesMutualInformation""TransformBendingEnergyPenalty")。而其他的参数设置和单测度配准相同，无需赘述。本文使用Mattes互信息+BEP作为联合测度，使用联合测度时需要对每个测度进行权重设置，本实验对两个互信息的权重分别设置为0.17和0.83，同时改进了一些重要参数，例如NumberOfHistogramBins从64改为了50，还对分辨率层级的分配做了调整，此外因为原图像有很多无用因素，比如CT机床架的图像，所以同时还使用了图像掩模来进一步增加配准精度，详细参数文档请参见附录C。本次实验结果如下：图3.11第一位病人的联合参数配准实验结果本次实验耗时较长，程序运行了1小时4分27秒，初始互信息为-0.735267，最终互信息为-0.782261。由图可见，本次配准效果很好，且没有出现肉眼明显可见的非实际形变，可以说本次实验比刚性以及弹性配准都要成功。下图是本次联合参数配置过程中的互信息变化曲线：图3.12第一位病人的联合测度配准过程中互信息变化曲线从互信息曲线可以明显看到互信息强烈收敛，从此可以说本次对第一位病人的配准是比前两次都成功的，这说明了联合测度配准确实对提高配准精度有很大的作用，但也如预期一样，联合测度牺牲了配准的速度，需要很长时间才能完成配准。第二位病人的联合测度弹性配准实验将同样的参数文档用到第二位病人的联合测度配准中，得到了以下的结果：图3.13第二位病人的联合测度配准实验本次实验耗时32分50秒，初始互信息为-0.723749，最终互信息为-0.821640。从图像上来看，本次实验较前两次成功，但是仍有非实际形变产生。下面观察其互信息变化：图3.14第二位病人的联合测度配准过程中互信息变化曲线由互信息变化曲线可以明显看到互信息强烈收敛，说明本次配准确实是比较成功的，但是从观察上来看仍有非实际形变，参数文档仍有改进的空间。4结论本文主要对现有配准中已存算法进行了实验探究，分析刚性、弹性以及联合测度配准的性能特点。通过实验得出以下结论：刚性配准速度快精度可靠，但是无法处理形变，因此使用的范围有限，大多只能作为后续配准的初始化变换。弹性配准可以处理形变，但是需要仔细设置参数，否则会出现非实际形变。联合测度配准可以使配准的精度更上一层楼，但是需要消耗很多计算时间，因此在时间条件充裕的情况下可以通过此法来提高配准精度。由于时间和自身水平有限的关系，本文还存在很多缺陷，难以探究配准的深层次问题，本文的不足之处主要有以下几点：没有联合配准金标准来进行配准结果的评估，而仅仅只用可视化以及互信息测度来分析配准是否成功。因为对配准知识的理解不够透彻，不能从根源理论来指导实验，这使得实验效率不高。软件水平有限，无法进行对开源软件的个性化编写以满足自己的实际需求，譬如设计一个可视化的快速配准窗口。通过本次毕业设计，我收获了很多，学到了很多，首先学会了文献的查阅，以及跨平台开源软件的安装和使用；还学到在实际中遇到问题以及解决问题的方法和手段；最后学到了如何写出一篇规范的论文。这些宝贵经验无疑对我以后的学习和工作都有巨大的帮助。致谢光阴似箭，转眼间已经到了毕业季，大学最后一段时间是在毕业设计中度过的，这段时间是我人生中最宝贵的一段岁月，这次毕业设计我付出了很多，也收获了很多，积累了宝贵的学习和人生经验。但是自己一个人是不可能完成这些的，这都归功于大家的帮助和支持。首先要感谢我的毕业设计导师刘国才教授，从选题开始，到最后的论文的写作，刘国才教授都给予了我极大的帮助，细心指导毕业设计中的每一个环节存在的各种问题，没有刘国才教授全面指导，我相信是无法完成本次毕业设计的。此外还要感谢聂茂学长，在毕业设计的过程中，无论是理论上的问题，还是实验操作中的困难，聂茂学长都给予了细致的讲解和启发，让我在毕业设计过程中少走了不少弯路。同时感谢胡浩，黄睿，洪国凯，康志鹏等同学在毕业设计期间的大力支持和协助，通过积极讨论毕业设计中的各种问题，我们共同进步，共同克服毕业设计中存在的问题。感谢父母和家人，是他们无私的爱支持着我，让我我不断地学习和进步。最后感谢学校的各位领导、老师和专家，感谢你们在百忙之中抽出宝贵的时间对我们的论文进行审批。参考文献[1］俞亚青，田学隆，闫春红．医学图像配准方法分类及现状[J]．重庆大学学报,2003，26(8)：114-118．[2]王苏苏．基于CNN和互信息的PET/CT图像配准方法的研究[D]．青岛：青岛大学，2013．[3]SchlemmerHP,PichlerBJ,SchmandM,etal.SimultaneousMR/PETimagingofthehumanbrain:feasibilitystudy[J].Radiology,2008;248:1028-1035.[4]ChuanxinNiu，MedicalimageregistrationbasedonmutualinformationusingKrigingprobabilitydensityestimation[J].IEEETransactionsOnEngineeringinMedicineandBiologySociety,2005:3097-3099.[5]StefanKlein,MariusStaring,KeelinMurphy,MaxA.Viergever,andJosienP.W.Pluim.Elastix:atoolboxforintensity-basedmedicalimageregistration.IEEETRANSACTIONSONMEDICALIMAGING,2010:196-205.[6]程莹.基于细胞神经网络的应用研究[D].南京航空航天大学,2008.[7]张存利.基于互信息的多模态医学图像CT_MRI配准[D].长春：吉林大学，2011.[8]胡永祥.基于互信息的多模态医学图像非刚性配准研究[D].长沙：中南大学,2012.[9]陈伟卿.基于互信息的医学图像刚性配准研究[D].大连：大连理工大学，2010.[10]别术林.基于互信息的医学图像配准算法研究[D].北京：北京交通大学，2014.[11]张倩.基于互信息的医学图像配准算法研究[D].济南：山东大学，2008.[12]徐显辉.小鼠CT图像非刚性配准方法研究[D].西安：西安电子科技大学，2014.[13]彭文.基于特征的医学图像配准中的若干关键技术的研究：[D].杭州：浙江大学，2007.[14]ValerioFortunati,RenéF.Verhaart,FrancescoAngeloni,AadvanderLugt,WiroJ.Niessen,JifkeF.Veenland,MargarethusM.Paulides,TheovanWalsum.FeasibilityofMultimodalDeformableRegistrationforHeadandNeckTumorTreatmentPlanning[J].InternationalJournalofRadiationOncology•Biology•Physics,

July8,2014,Vol.90,

Issue1,

p85–93.[15]SaraLeibfarth,

DavidMönnich,

StefanWelz,

ChristineSiegel,

NinaSchwenzer,

HolgerSchmidt,

DanielZips,and

DanielaThorwarth.AstrategyformultimodaldeformableimageregistrationtointegratePET-MRintoradiotherapytreatmentplanning[J].ActaOncologica,2013;52:1353–1359.附录附录A刚性配准参数文档============startofParameterFile:Par0023input/Rigid.txt=============(FixedInternalImagePixelType"float")(MovingInternalImagePixelType"float")(UseDirectionCosines"true")//****************MainComponents**************************(Registration"MultiResolutionRegistration")(Interpolator"BSplineInterpolator")(ResampleInterpolator"FinalBSplineInterpolator")(Resampler"DefaultResampler")(FixedImagePyramid"FixedSmoothingImagePyramid")(MovingImagePyramid"MovingSmoothingImagePyramid")(Optimizer"AdaptiveStochasticGradientDescent")(Transform"EulerTransform")(Metric"AdvancedMattesMutualInformation")//*****************Transformation**************************(AutomaticScalesEstimation"true")(AutomaticTransformInitialization"true")(HowToCombineTransforms"Compose")//*******************Similaritymeasure*********************(NumberOfHistogramBins64)(ErodeMask"false")//********************Multiresolution**********************(NumberOfResolutions3)(ImagePyramidSchedule882441111)//*******************Optimizer****************************(MaximumNumberOfIterations500)(MaximumStepLength1.0)(RequiredRatioOfValidSamples0.05)//****************Imagesampling**********************(NumberOfSpatialSamples2000)(NewSamplesEveryIteration"true")(ImageSampler"Random")//*************InterpolationandResampling****************(BSplineInterpolationOrder1)(FinalBSplineInterpolationOrder3)(DefaultPixelValue0)(WriteResultImage"true")(ResultImagePixelType"short")(ResultImageFormat"mhd")============endofParameterFile:Par0023input/Rigid.txt============附录B弹性配准参数文档=========startofParameterFile:Par0023input/Deformable_alter.txt=========(FixedInternalImagePixelType"float")(MovingInternalImagePixelType"float")(UseDirectionCosines"true")//****************MainComponents**************************(Registration"MultiResolutionRegistration")(Interpolator"BSplineInterpolator")(ResampleInterpolator"FinalBSplineInterpolator")(Resampler"DefaultResampler")(FixedImagePyramid"FixedSmoothingImagePyramid")(MovingImagePyramid"MovingSmoothingImagePyramid")(Optimizer"AdaptiveStochasticGradientDescent")(Transform"BSplineTransform")(Metric"AdvancedMattesMutualInformation")//*****************Transformation**************************(FinalGridSpacingInPhysicalUnits16)(HowToCombineTransforms"Compose")(BSplineTransformSplineOrder1)(UseCyclicTransform"true")//*******************Similaritymeasure*********************(NumberOfHistogramBins64)(ErodeMask"false")//*