基于ITK的人脑图像配准算法实现

南京邮电大学 毕 业 设 计（论 文）题 目基于ITK的人脑图像配准算法实现专 业生物医学工程学生姓名曾宪波班级学号B B指导教师杜若瑜指导单位南京邮电大学地理与生物信息学院 日期：2015年1月15日 至 2015年6月7日毕业设计（论文）原创性声明本人郑重声明：所提交的毕业设计（论文），是本人在导师指导下，独立进行研究工作所取得的成果。除文中已注明引用的内容外，本毕业设计（论文）不包含任何其他个人或集体已经发表或撰写过的作品成果。对本研究做出过重要贡献的个人和集体，均已在文中以明确方式标明并表示了谢意。 论文作者签名： 日期： 年 月 日摘 要伴随着科技的发展与进步，医学图像成像技术也在飞速发展，图像处理在医学研究与临床医学中的应用越来越广泛，涌现出大量的医学图像。但是，大量的医学图像因为成像设备、环境、时序、位置等的差异带来了医学图像信息的不一致性，如何消除里面的差异获取一致性的信息成为一个函待解决的问题。本文介绍了医学图像配准构架中的各个部件以及常见的配准算法。同时设计了2组实验：使用2组人脑MRI图像序列进行配准研究。利用ITK中的居中变换刚性配准进行2D刚性配准；将序列重建为3D体数据，利用Versor 3D刚性配准进行3D刚性配准。通过对2组实验结果进行比较与分析，得出2D刚性配准在配准精度上优于3D刚性配准，但配准耗费的时间比3D刚性配准要多的结论。关键词：ITK；人脑图像；图像配准；刚性配准；MRIABSTRACTWith the development of science, medical imaging technology is expanding rapidly. As image processing is applied widely in medical research and clinical medicine, a large number of medical images are emerging. However, because of the differences in the imaging device, the environment, the timing, location, a majority of medical images has led to the inconsistency of medical image information. Therefore, how to eliminate the differences and obtain consistency of information has become a problem to be addressed.Thisarticleintroduceseachcomponentofmedicalimagingregistrationframeworkandsomecommonregistrationalgorithms.Inthemeantime,itdesignstwocomparativeexperimentstodoregistrationresearchbasedontwosetsofbrainMRIimagesequence.Firstly,itrebuildsthesequenceas3DdataandthenmakesuseofthecentertransformationandrigidregistrationofITKtoconduct2Drigidregistration,aswellasVersor3Dtoconduct3Drigidregistration.Throughanalyzingandresearchingtheresultsoftwoexperiments,itconcludesthat2Drigidregistrationisbetterthan3Drigidregistrationontheaccuracyofregistration,butweakerontheconsumingtimeofregistration.Keywords：ITK; brain image; image registration; rigid registration; MRI目 录第一章 绪论11.1 课题意义11.2 研究现状简述11.2.1 国外研究现状11.2.2 国内研究现状11.3 医学图像配准的问题21.4 本文的主要工作2第二章 ITK简介与配准框架32.1 ITK的安装32.1.1 获取安装资源32.1.2 ITK安装32.2 ITK的背景和理念32.3 ITK核心概念32.3.1 数据表达32.3.2 泛编程式42.3.3 数据处理管道42.3.4 内存管理与智能指针52.3.5 事件对象62.3.6 空间对象62.4 DICOM文件结构62.4.1 DICOM文件头62.4.2 DICOM数据集合7第三章 医学图像配准83.1 图像配准的分类83.2 医学图像配准的基本过程93.3 医学图像配准框架各部分理论介绍93.3.1 空间变换93.3.2 相似性测度（度量）113.3.3 优化器133.3.4 图像插值133.4 医学图像配准的评估14第四章 配准实验154.1 实验准备164.1.1 实验条件164.1.2 DICOM文件的读取164.1.3 序列DICOM文件转换为3D体文件164.1.4 3D体文件沿水平方向切片174.2 2D刚性配准174.2.1 2D刚性配准介绍174.2.2 输入数据194.2.3 实验结果194.3 3D刚性配准204.3.1 Versor3D刚性配准介绍214.3.2 输入数据224.3.3 实验结果224.4 对比研究24结束语25致谢26参考文献27南京邮电大学2015届本科毕业设计（论文）第一章 绪论1.1 课题意义随着我国现代化医学影像设备技术的发展，高科技高技术含量的医学影像设备被大量的应用于临床诊断和治疗中。在实际的临床应用中，单一模态的图像一般无法提供医生所需要的足够的信息，一般需要把不同模态的图像融合在一起，得到更丰富的信息来了解病变组织的综合信息，从而做出准确的判断1。所以，在临床上一般需要对同一病人进行多种模式或同一模式的多次成像，从而能够同时从多幅图像进行综合分析，获得有价值的信息。对多幅不同的医学图像进行分析，首先要解决的问题就是对这几幅图像的严格对齐问题，即是医学图像的配准问题，医学图像配准就是指对一幅医学图像进行一种空间变换，让它与另一幅医学图像上的对应点达到空间位置的一致，即：通过变换将图像与另外一幅图像在相应空间点的位置达到相同，几幅经过配准后的多模态图像可以为后期的医学图像融合做准备2。1.2 研究现状简述1.2.1 国外研究现状国外早在70年代就已经开展了医学图像配准的研究，提出了各种各样的配准方法。根据Institute for Scientific Information的调查数据表明，在1999到2008年十年里至少有不少于3900多篇的学术论文在研究图像配准问题。关于图像配准和相关的专利申请，在美国2006年网上公布的申请专利项目中，就多达50多项3。很多大公司已经意识到这方面的前景，都开始设立专门的工作组从事于图像配准方面的研究。针对医疗图像配准广阔的应用前景，在发达国家的高校和企业均加大投入和研究的力度，开发了很多包含配准功能的高质量医学影像软件和算法平台，使图像配准能够实现自动化和模板化。例如德国汉堡大学开发的Voxel-Man项目实现了高质量人体结构医学图像分析及可视化。德国癌症研究中心开发了一个免费的开源的医学影像交互工具包（MITK），能够实现多模医学图像的配准工作4。不过应用较为广泛的，当属ITK（洞察分割与配准工具），1989年，美国国家医学图书馆进行人体可视化项目的研究，想开发一个分割与配准的程序包，作为visible human项目的一个工具。ITK就是源于这个项目而形成的。ITK是一个开放源码的类库，所以它既可以自己维护，又可以供使用者和开发者不断交流，从而不断扩充该软件。ITK实现了大部分的配准和分割算法，形成一个算法仓库，避免重复劳动，减少开发者的障碍。1.2.2 国内研究现状21世纪随着科技的不断进步，医学影像设备也被广泛的应用，巨大的商机必然引起各个相关行业的重视。虽然我国在医学影像方面起步较晚，基础薄弱，但是随着我国综合科技水平的提高，各大科研机构和院校都在医学影像方面做了大量有益的尝试。如中国科学院自动化研究所人工智能实验室朱付平、田捷、林瑶、葛行飞等人提出了基于Level Set方法的医学图像分割5；四川大学建筑学院生物医学工程研究所袁华、滕奇志、综述、袁支润、陶德元等人提出了边缘检测及其在医学图像中的应用6；曹会志、王晨、罗述虚等人提出结合蚁群算法的Snake模型的医学图像分割的方法7；林瑶、田捷、张晓鹏等人提出了基于模糊连接度的FCM分割方法在医学图像分析中的应用8；中科学自动化研究所的医学影像处理研究所开发了C+类库-MITK，并给予MITK编写了3D医学图像处理分析系统-3DMed,并配套建立了技术支持和社区交流的专用网站，是医学影像算法研究与应用的典范9。1.3 医学图像配准的问题医学图像配准基础理论研究庞杂并且分散。医学图像配准涉及医学成像技术、图像预处理、图像插值、空间几何变换、解空间搜索技术、相似性度量等多方面的知识。虽然这些相关知识都有各自相关专著进行相应的论述，但是具有针对性的统一论述任然非常匮乏；医学图像配准理论相关专著大多是有英文写成，作为图像配准研究领域的基础，这些问题对计算机专业研究者入门会造成非常大的障碍。医学图像配准选用不同的方法配准的结果并不相同，所以应当根据医学图像的特性选择合适的算法来实现配准功能，使图像达到配准的最佳效果。如何选择配准算法以及设置算法中的参数是医学图像中的重要问题。1.4 本文的主要工作本文主要工作在以下几个方面：进行有针对性的理论总结。通过大量阅读国内外该领域的资料，本文分析了国内外医学图像配准的现状，并对医学图像配准框架以及医学图像配准理论进行了总结。介绍2D刚性配准与3D刚性配准原理与步骤，设计实验并进行编程：选取人脑切片为实验对象，完成对人脑图像配准的实验，分别使用2D刚性配准和3D刚性配准进行人脑切片图像的配准处理。并对2D刚性配准与3D刚性配准的实验结果进行了比较和评价，分析了两种配准算法的优劣。第二章 ITK简介与配准框架2.1 ITK的安装2.1.1 获取安装资源（1）Cmake安装资源CMake是一个跨平台的安装（编译）工具，可以用简单的语句来描述所有平台的安装(编译过程)。用于对ITK进行编译，可以在ITK官网/download/进行下载，本文选用的是在Windows平台的cmake-3.2.2-win32-x86.exe版本。（2）ITK安装资源ITK安装资源可以在ITK官网/进行下载，最新版本为4.7.1.，本文采用最新版本资源InsightToolkit-4.7.1.tar.gz。上述的安装基于Windows 7 系统和 Visual Studio 2012 开发平台。2.1.2 ITK安装本文参考/lj/article/details/教程进行安装，ITK安装路径 D:ITK，Cmake安装路径 C:Program Files (x86)CMake 。2.2 ITK的背景和理念a.ITK的设计理念主要基于泛型编程思想；能够提供用来分割和配准的数据表示和算法，主要针对医学应用，但是同样也可以用来处理其它类型的数据；为任意维度图像的数据表示提供统一的格式；不提供可视化或者图形化用户界面，可视化或者图形化用户界面由其他软件完成（如 VTK 等）；开发基于极限编程的原则；支持多线程并行处理；对象工厂用来实例化对象，它允许实时系统扩展；使用Command/Obeserver设计模式处理事件；组织围绕着数据流结构，即它的数据是用数据对象表示的，由一系列处理对象轮流处理。数据对象和处理对象由管道状结构组织在一起。管道可以根据用户定义的存储限制分块处理数据。b.ITK 的应用理念跨平台,支持Unix、Windows和MacOSX系统；由于使用了泛型编程思想，大量的C+模板代码可能会极大增加编译的负担；ITK支持多语言绑定，包括 Tcl,Python,和Java等，绑定是由自动封装机制自动生成的，使用CableSwing进行代码封装；使用了智能指针；使用CMake构建，它不依赖于特定系统和编译器，根据操作系统和所用编译器来有针对性的构建所需文件。2.3 ITK核心概念2.3.1 数据表达在ITK中，有2种基本的数据表达类型：图像和网格。这一功能是通过类Image和类Mesh实现的，Image和Mesh类均为itk:DateObject的子类。在ITK中没数据对象是那些分送系统和参与数据流管道的类10。若要完整的表示一幅图像，则需要提供图像的相应的拓扑结构和像素值的信息。因图像的拓扑结构是规则的，因此为表示一幅图像的结构，至少需要提供图像的原点坐标、沿各个坐标轴方向上的点的间隔和尺寸。图像的信息包括：图像的维数、起点坐标、像素的间隔、图像的尺寸和像素值。 其中坐标、尺寸和间隔的数据类型都是固定的，分别是整形、浮点型和整形。因此剩下的图像维数和像素类型就唯一缺点了一个具体的图像类型。ITK中关于图像的一个非常重要的概念是区域，它是一个连续的、矩形的图像块。区域用来指定图像中处理的部分，如多线程或阿伯村在内存中的部分。ITK支持多种像素类型，除了基本数据类型外，还支持复合类型，如有三个分量的RGB类型和任意分量的Normals、Vectors类型等，其中各分量的类型又可以是任何基本数据类型。网格类表示一个n维无结构的格子。它的拓扑布局是由一套单元表示的，单元是由一个类型和连通性列表来定义的；该连通性列表依次设计点，网格的几何结构是通过n维的点组合起来相关的单元插补函数来进行定义11。点的维数，和点、单元和单元边界相关联的一个像素类型及一个网格的特征模板参数决定了网格的形式。2.3.2 泛编程式ITK在它的实现中使用了泛型编程的技术。范型编程的基本思想是：保持数据的容器、储存数据的迭代器以及范型算法使用容器和迭代器来创建有效的基本算法。该方法的优点是能够通过定义相近的模板类型可以支持几乎无限的数据类型。缺点是有很多编译器仍然不支持这些更新的概念，而且不能在 ITK 中编译，即使它们支持这些概念，也极有可能产生完全不可辨认的错误信息，哪怕仅仅是因为最简单的语法错误造成的。2.3.3 数据处理管道数据对象（如图像和网格）是用来代表数据的，过程对象是用来操作数据对象并且产生新的数据对象的类。像filter object、sourse和mappers都是过程对象12。 Sourse（例如readers）生成数据，filter object载入数据并经过处理产生新的数据，而mappers接收数据并输出到一个文件或其他系统中。数据处理管道的作用是连接数据对象（如图像和网格）和过程对象。管道支持一个自动更新机制。这个机制在当且仅当它的输入或内在状态改变时才会唤起一个滤波器来运行。另外，数据管道还支持streaming，可以自动将数据分成许多小部分，对这些小部分进行逐一处理，并重新集合处理后的数据产生最终结果。ITK的数据处理管道作业方式如图所示。如图示，数据流从最左边读入图像，接着依次从左到右经过管道，传递到最右边作为输出类型。通常使用SetInput( )和GetOutput( )方式连接数据对象和过程对象，方法代码如下：typedef itk:Image FloatImage2DType;itk:RandomImageSource:Pointer random;random = itk:RandomImageSource:New( );random-SetMin(0.0);random-SetMax(1.0);itk:ShrinkImageFilter:Pointer shrink;shrink = itk:ShrinkImageFilter:New( );shrink-SetInput(random-GetOutput( );shrink-SetShrinkFactors(2);itk:ImageFileWriter:Pointer writer;writer = itk:ImageFileWriter:Pointer:New( );writer-SetInput (shrink-GetOutput( );writer-SetFileName( test.raw );writer-Update( );在该例中，源代码对象itk:RandomImageSource同itk:ShrinkImageFilter相连接，shrink滤波器和mapperitk:ImageFileWriter相连。当调用Update( )方式时，数据处理管道便按顺序使用这些滤波器，并且将最终结果储存到文件中。如图2.1所示：图2.1 数据处理管道2.3.4 内存管理与智能指针系统面向对象的、天然的特性通过各种各样的对象类型或类来描述和操作数据。当对一个特定的类实例化生成这个类的实例时，系统将自动分配内存来存储数据属性值和方法指针。其他类或数据结构便可以通过普通程序操作来引用这个对象。一般地，在程序执行过程中，当实例被删除或内存恢复时，应用将会消失。但是，确定实例在什么时间段删除是非常困难的。删除实例过早会造成程序崩溃；而删除时间太晚会出现过多消耗内存或内存泄露。这种内存分配和释放的过程就是内存管理。在ITK中，通过引用计数来实现内存管理。在引用计数中，每个实例的应用次数将被记忆。当对象的引用计数为零时进行释放。在垃圾回收中，一个后台处理系统将自动识别并释放不再被系统引用的对象实例13。与垃圾回收相关的问题是内存删除的确切时间。当对象数据很大时是不可接受的。一旦对象的所有引用消失引用计数将迅速删除内存。通过一个Register( )/Delete( )成员函数界面来实现引用计数。引用计数初始值为0，当创建了一个对象或者某一对象被其它对象应用的时候，调用一个Register( )来对实例的引用计数增加1，当实例引用消失时，调用一个Delete ( )方式来对实例的引用计数减少1。当引用计数为0时，实例将消失。使用一个称为itk:SmartPointer的类来进行引用计数可以大大的简化引用计数操作操作，这样的类称为智能指针。与ITK中大部分实例不同，智能指针可以指向程序块，当智能指针生成的范围关闭时会自动删除。因此，应该尽量少Register( )和Delete( )在ITK中的使用。2.3.5 事件对象ITK中的使用Subject/Observer设计模式来执行事件处理。在该方式中，对象通过注册它们正在观察的实例，来关注这个特定实例所调用的一个特定的事件。如：ITK中的滤波器周期性地调用itk:ProgressEvent。事件发生时对象已记录下在这一事件中通报的信息。在记录过程中，通过调用一个命令来发出通告是特定的。2.3.6 空间对象图像分割和配准的理论是ITK空间对象框架支持对象处理的任务。ITK空间对象提供了一个界面，用于储存对象的物理位置、几何性质以及对象之间的关系，储存是在独立于用来表示这些对象的形式的情况下进行的。空间对象框架提供的功能支持对象配准、分割、表面/立体视图以及其他分析和显示功能。空间对象框架扩充了对计算机视图工具包很普遍的一个“scene graph”概念，因此也支持这些新函数功能。目前已经实现的空间对象类型包括：Ellipse、Blob、Image、Group、Line、21Tube和Surface。对象itk:Scene用于控制一个个依次含有子对象的空间对象列表。每一个空间对象都用一个色彩特征标记，每个空间对象都有自己的功能。ITK中的空间对象和它们的方法有限，但它们的数量正在增长且潜力很大。2.4 DICOM文件结构DICOM文件是由一个头文件与一个图像数据体组成。头文件中包括了标准和自由形成域。标准化域的集合叫作公共DICOM词典，自由域的列表也叫作shadow dictionary。通常一个简单的DICOM文件能够包含多框架，允许存储体和生理活动14。图像数据可以使用许多种标准压缩，包括LZW、JPEG和RLE。2.4.1 DICOM文件头DICOM文件头(DICOM File Meta Information)包含了标识数据集合的相关信息。每个DICOM文件都必须包括该文件头。文件头的最开始是文件前言，它由128个00H字节组成，接下来是DICOM前缀，它是一个长度为4字节的字符串“DICM”，可以根据该值来判断一个文件是不是DICOM文件。文件头中还包括其它一些非常有用的信息，如文件的传输格式、生成该文件的应用程序等等。a. 除了128字节的文件前言和4字节的DICOM前缀外，所有其它的文件头元素都必须采用显示格式编码，各个数据元素排列的顺序按照标签数值从小到大的传输格式(Little Endian)编码。b. 每个文件头元素的长度必须为偶数，否则应该按照规定补充一个字节。c. 有(0002，*)类的标签都为DICOM所保留。为了兼容后续版本，如果发现文件中有目前尚未规定的(0002，*)类标签，则应该忽略它。2.4.2 DICOM数据集合a. DICOM文件主要组成部分就是数据集合。这不仅包括医学图像，还包括许多和医学图像有关的信息，如病人姓名、图像大小等。b. DICOM数据集合是由DICOM数据元素按照指定的顺序依次排列组成的。对于DICOM文件，一般采用显式传输，数据元素按标签从小到大顺序排列，即DICOM PS 3.5规定的Explicit VR Little En-dian Transfer Syntax。c.在DIOCM标准的PS 3.3部分(Information Object Defini-tions)中，定义了各种类型的图像文件必须包括和可选的DICOM数据元素，在制定自己的DICOM文件结构时，必须严格遵照该部分规定。第三章 医学图像配准3.1 图像配准的分类到目前为止，医学图像配准方法还没有形成一个统一的标准，可以根据不同的分类标准将同一配准分类到不同的分类里面。（1） 根据医学图像的模态分类根据医学图像的模态分类可将配准方法分为多模态配准和单模态配准。来自同一成像设备的图像配准一般为单模态图像配准。常用于观察病灶。来自不同成像设备的图像配准一般为多模态配准，主要应用在诊断方面。（2） 按照配准维数分类按照配准图像的维数，配准可以分为2D与2D，2D与3D，3D与3D配准。在断层扫描数据的配准上一般是2D与2D配准；在2D片层扫描数据与3D空间数据的配准时一般使用2D与3D的配准；3D图像的配准一般都使用3D与3D的配准方法。当配准的维数变大时，配准的复杂性将会大大提高。（3） 根据变换性质分类主要可以分为仿射变换、刚体变换、投影变换和曲线变换四种。仿射变换就是将坐标变换和尺度变换结合；刚体变换主要包括坐标轴的平移和旋转；投影变换是将直线映射成直线，平行性和垂直性不能保证；曲线变换是将直线映射为曲线。（4） 根据图像来源进行分类根据图像的来源将配准方法分为同一对象的配准和不同对象的配准以及对象与图谱的配准。待配准图像都来自同一对象的配准称为同一对象的配准；待配准图像来自不同的对象称为不同对象的配准；待配准图像一幅来自对象，一幅来自图谱称为对象与图谱的配准。（5） 按照配准的基准进行分类按照配准的基准进行分类可以分为基于像素的配准和基于特征的图像配准。基于像素的配准无需进行特征的提取与图像分割，且配准过程中不需要人工进行干预，但是计算量大，运行速度低，适应性不佳。基于特征的图像配准又可以细分为基于基于内部特征的配准和基于外部特征的配准。利用体内的解剖特征或集合特征进行空间变换进行配准的称为基于内部特征的配准；在患者身体上设置一些标记点，这些标记点能在不同的影像模式中显示，使用不同的方法用标记对图像进行配准称为基于外部特征的图像配准。（6） 根据用户交互性的多少进行分类根据人工在配准中的参与度，可将配准分为自动配准、半自动配准和交互配准。自动配准不需要人为干预自动进行；半自动配准需要人为设定配准参数；交互配准是软件提供了如读入图像等的交互按钮，人工选择后实现按照设定的方向进行配准。（7） 根据配准过程中变换参数确定的方式进行分类一般分为：能够直接利用公式得到变换参数的配准与参数空间中寻求目标函数的最优解得到空间变换参数的配准。基于特征信息的配准一般使用利用公式得到变换参数的配准，求目标函数的最优解得到空间变换参数的配准一般用于将配准函数看做一个能量极值求解的配准中，并且极值能够用坐标算法求解出来。目前常用的优化算法有梯度下降法、无约束Powell算法、模拟退火法等。3.2 医学图像配准的基本过程医学图像配准的基本流程有以下几个步骤，如图3.1所示：图3.1 医学图像配准的基本流程由图可知，配准框架的基本成员：两个输入图像、一个变换、一个相似性度量、一个插值器、一个优化器。在整个配准过程当中，不断地搜索空间变换参数，使评价配准结果的度量达到最优，这个度量就是相似性度量，通过优化搜索使得相似性度量达到最大，这时得到的配准参数就是最优参数。（1）输入两幅图像参考图像，将一幅固定不变的图像定义为参考图像f(X)，另外一幅图像我们定义为待配准图像m(X)，其中X表示在N维空间中的一个位置，将提取图像的特征信息构成特征空间；（2）通过提取出黎的特征空间来确定一种空间变换（T），并且通过插值得到待配准图像在区域T1(X)的取值m(T1(X)；（3）在相似性度量模块计算插值图像m(T1(X)和参考图像f(X)能够达到定义的相似性度量，即f（X）=T（m(T1(X)）；（4）在不断的空间变化过程中，加入一些优化方法搜索到适合的配准参数使相似性度量更好更快地达到最优解。重复2和步骤3直到相似性度量达到最大，终止迭代循环。其中变换T表示参考图像中的点和配准图像中对应点的空间映射。插值器用来评估待配准图像在非网格位置上的像素的灰度值。相似性度量是评价匹配结果好坏的标准。按照这个标准优化器在本变换参数定义的搜索空间进行优化，寻找符合这个标准的最优函数值。3.3 医学图像配准框架各部分理论介绍3.3.1 空间变换配准问题也可以总结为一个图像空间中的点映射到另一个图像空间中的点的方法。一般由一组参数来确定空间变换。基于所用的配准标准找出最优变换参数是配准的目标。一般来说，自由度越高或交换的次数越多，寻找最优解的的问题就越困难。对于2D与2D或者2D与3D配准问题，空间变换可能是刚性的、非刚性的或者仿射的。在刚性变换中，只能够进行转动和平移的操作。仿射变换只能够在转动和平移基础上的进行倾斜和缩放操作。而非刚性变换是自由映射，通常需要设置一定的规则来限制解空间的范围。在获取结构图像的过程中没有发生形变的情况下一般使用刚性变换，如人脑图像，因为有颅骨的结构避免了非刚性变化的发生。而对于缩放系数未知的刚体结构，一般使用仿射变换。目前刚体变换的技术已经相对成熟；而非刚性变换技术因为它的复杂性，相对刚性变换来说目前并不成熟。a.刚体变换如第一幅图像中的两点间的距离变换到第二幅图像中后依然保持不变，那么这种变换称为刚体变换。刚体变换可分解为平移、旋转和反转。在2D空间中，点(x，y)经过刚体变换到点(x，Y)的变换公式为：（3.1）其中为旋转角，为平移向量。本文以人脑图像为研究对象进行配准，因为人脑一般由颅骨保护不易发生形变，故选用刚性变换。b.仿射变换如经过变换后第一幅图像上的直线映射到第二幅图像依仍为直线，且保持平衡关系，这样的变换称为仿射变换。仿射变换可以分解为线性变换和平移变换。在2D空间中，变换公式为：（3.2）其中为实矩阵。c.投影变换如经过变换后第一幅图像上的直线映射到第二幅图像上依仍为直线，但平行关系基本不保持，这样的变换称为投影变换。投影变换可以用高维空间上的线性变换来表示。变换公式为：(3.3)d.非线性变换非线性变换可把直线变换为曲线。在2D空间中，可以用以下公式表示: (3.4)式中，F表示把第一幅图像映射到第二幅图像上的任意的一种函数形式。典型的非线性变换如多项式变换，在2D空间中，多项式函数可写成如下形式：(3.5)非线性变换一般适用于那些具有全局性形变的图像配准问题，以及整体近似刚体但局部有形变的配准情况。在医学图像配准中，对刚体变换的应用较多，尤其是骨骼及脑部的配准。其他的变换形式多用于软组织的配准中。3.3.2 相似性测度（度量）评价配准结果好坏的标准是相似性测度，一般定义为某种代价函数的形式，在给定的特征空间或变换模型下得到有关变换的参数。在ITK中，类 itk:ImageToImageMetric从数量上通过比较图像亮度的灰度值来待配准图像与参考图像的质量差别配准框架中最为关键的部分为相似性测度，应针对不同的配准问题选择不同类型的相似性测度。它的基本输入是一是：参考图像和待配准图像、一个变换函数、一个插值器。相似性测度有多样的形式，可以分为基于特征的和基于体素的。比较经典的相似性测度包括相关函数、互信息、均方误差等。在ITK中应用较多的是以下基于体素的相似性度量：a. 均方法相关尺度均方相似性测度计算图像A和B中用户定义的区域的方形像素在亮度上不同程度的均方值，公式如下： (3.6)其中是图像A的第i个像素；是图像B的第i个像素；N是涉及的像素的数目。相似性测度的理想值是零。图像A与B搭配较差的情况会使metric的值较大。均方相似性测度依靠在两幅图像中相同的对应点的亮度是相同的假设，因而它的用途是限制相同模态的图像。另外，任何在亮度上的线性改变都会导致较差的搭配。通过不断的平移交换，比较均方相似性测度的大小获取两幅图像的最佳搭配。b.标准化相关度这种度量计算是基于像素的配准图像间的相关度，然后分别利用各个图像相关的平方根来把它标准化。当两幅图像完全一致的时候，这个度量等于1。图像的非正确比使得这个度量值小于1。该度量方法对图相间的常量乘数不敏感。它的数学表达式如下： (3.7)式子乘以-1的目的是为了能够使这个度量可以与最小值搜寻优化器配合使用。当映射位置落在移动图像域外面时，这个方法和均方法一样，把这样的像素对的作用因数抛弃。关于坐标系中度量的曲线，除在零平移位置有微小差异外，这个度量的曲线和均方法得到的结果很相似。c.在差分密度度量方法中，每个像素的影响度用计算得到，其中d是固定图像和移动图像之间图像强度的差分。方程是钟形的，在d=0的时候取得最大值1，在d=的时候取得最小值为0。其中常量控制下降的速度，较小的的将会导致较快的下降速度。这里的对应的是下降50%时候的强度差分。对的选择依赖配准中两幅图像的强度。差分密度的数学表达式如下： (3.8)该度量的最优值为N,较差的比将会对应比较小的度量值。这个度量可以解释为大致一致的差分图像的概率密度函数，最大值取于0值的为止。由于这个度量基于两幅图像的强度差分，因此它对线性差分不敏感。这个度量曲线的轮廓有很尖的峰。尖峰的好处是它可以快速缩小优化范围，然而由于度量倒数在峰值处非常大，这可能会影响优化算法在接近最大值的性能，尤其是对那些依赖于度量导数减量的算法。d.互信息法互信息(Mutual Information)是信息论里一种有用的信息度量，它可以看成一个随机变量中包含的关于另一个随机变量的信息量，或是一个随机变量由于已知另一个随机变量而减少的不肯定性。互信息一般用熵来表示，根据熵的连锁规则，有（3.9）因此（3.10）在图像配准中，固定图像与移动图像的强度都是随机变量。使用互信息法的主要好处：并不需要具体给定两个随机变量之间的相互关系。所以，它可以较好的对两幅图像之间的复杂映射建模，这就使互信息法适用于多模配准中。在医学图像配准中，即使两幅图像可能是来自于不同的时间或者不同的成像设备，但是他们都是基于同一人体解剖信息。因此当两幅图像的空间位置达到完全匹配时，其中待配准图像表达的关于参考图像的信息，即对应像素灰度的互信息为最大。3.3.3 优化器配准过程可以看作是一个参数寻找最优解的问题，优化器被用来找到变换参数去最优化相似性测度。相似性测度就是在搜索空间求最小解的估价函数，得到的解即空间变换的参数。最优化过程从初始参数开始，通过搜索空间不断的循环性的变换位置、评估代价方程来寻求最优解。在ITK中有两种主要类型的优化器。在第一种类型中我们发现优化器适合处理返回单值的cost function。这些确实是最普遍类型的cost function，著名的有Single Valued函数。第二种类型是那些适合于管理每次阶跃中返回多值的cost function。这些cost function普遍用在模型搭配问题，比较著名的有Multi Valued或者Multivariate函数。在ITK中相应的优化器被叫作MultipleValued优化器选择合适的优化器是影响配准速率快慢，结果好坏的关键。常见的优化算法有：梯度下降法、Powell法、模拟退火法、遗传算法、黄金分割法、下山单纯形法等。应当更具不同的情况选择适应的优化算法，选择合适的优化算法有利于提高配准的速度和精确度。3.3.4 图像插值在配准过程中，通过一个变换把一个点从一个空间映射到另外一个空间的时候，一般情况下该点不会被刚好映射到网格区域，而一般被映射到非网格区域，因此我们需要使用插补法评估在映射位置的图像亮度。如图3.2所示：图3.2 参考图像的网格位置映射到待配准图像的非网格位置在配准中，差补操作在一次最优化循环里就要运行数千次。因此，选择差补方法时我们必须在计算简化和最优化的平滑度之间进行衡量。常用的插值方法有以下一些：a. 最近点差值最近点差值简单的用离指定方位最近的网格点上的像素强度作为自己的像素强度。这相当于是假定潜在的图像是块状分布的，在每个网格点之间强度过度是跳跃式进行的。这种插值策略代价比较低，而且不需要进行浮点运算。b.线性差值线性插值中，潜在图像被假定为线性块状分布的。即，网格位置之间的强度按线性规则变化。插值后的图像是空间连续但非可微的。线性插值是个邻居的加权和，公式为：（3.11)这里的n是图像空间的维度。权值是通过计算指定位置和邻居之间的距离得到：（3.12)c.B样条插值在B条样插值中，潜在的图像是使用基于B样条的方程表示的。B样条的系数可用快速递归过滤技术计算得到。在非网格位置上的像素强度可以用B样条系数和所求位置支持区域中的B样条影响函数核相乘得到:（3.13)B 样条影响函数核是可以拆分的:（3.14)对于序大于 1 的样条，插值后的图像和它的导数都是空间连续的。样条的序越高，需要用来计算插值的像素数越多。三阶 B 样条影响函数核比较常用，因为它可以较好的平衡平滑度和计算负担之间的关系。三阶B样条影响函数核的等式如下:（3.15)3.4 医学图像配准的评估医学图像配准所涉及的范围很广泛，使用到搜索算法、数字图像处理等相关知识，而且影响图像配准的因素比较多，不同的应用对配准的具体要求也不一样，这让医学图像配准技术研究和配准结果的评价都是一项很有挑战性的工作。因为在配准中所读入的医学图像数据一般是在不同的医学成像设备或者成像条件下获取到的，使得单模或者多模图像配准选用的评价配准结果好坏的标准并不具有唯一性。根据具体不同的应用，充分考虑精度、稳定性、准确度、可靠性、算法的复杂性和临床可用性等等影响配准的结果，选出特点条件下的最优标准，常用的评估方法有下面几种：a.体模体模的两种常用标准分别是硬件体模和软件体模。硬件体模一般利用CT扫描机等硬件设备进行扫描成像，通过对所得到的影响进行研究和分析，从而获得影像或者成像设备质量的量化评价。软件体模是计算机图像合成的结果。体模法是一种较为粗略的评估方法，所以只能够对医学图像配准的结果做初步的评价。b.准标多模医学图像配准结果精度的评估一般使用人工标记做准标。使用立体定向框架系统的体积图像数据来作为准标的数据，使用该数据可以用来评估所采用的配准方法的适用性。而立体定向框架系统一般是由立体定向图像获取、立体定向参考框架、探针或者手术器械导向几部分组成的，它定位准确，能够有效的避免因为图像畸变而可能带来的问题。c.目测检验直接使用目测观察区别，以人的目测的主观判断来评判两幅医学图像的配准的质量，优点是直观、快捷，缺点是不够客观、精确。d.最小均方误差（MSE）该方法就是在特定的约束条件下求均方误差最小。在医学图像配准中便是利用计算配准后的两幅医学图像的最小均方误差，以此为依据来评价配准结果的好坏。最小均方误差越小表示配准的效果越好。计算公式如下：（3.16)R(i,j)为配准后图像，O（i,j）为参考图像。第四章 配准实验4.1 实验准备前几章已经介绍了医学图像配准的基本知识和ITK系统结构，下面通过具体的实验来分析不同ITK配准方法的不同之处。4.1.1 实验条件本次实验的实验条件均为：CPU：英特尔 第二代酷睿 i5-2410M 2.30GHz 双核内存：8 GB ( 金士顿 DDR3 1600MHz )操作系统：Windows 7 旗舰版 64位 SP1运行环境：命令行ITK版本：InsightToolkit-4.7.1Cmake版本：cmake-3.2.2-win32-x86本次使用的数据为2组MRI人脑图像，每组为160张沿着H方向切片的人脑图像。图像来源：获取使用一台西门子T1加权脉冲序列采集三重奏扫描仪。使用以下设置，扫描平面的矢状面；切片方向；H；视野，220x220毫米；像素矩阵256x256：切片数量；160片；切片厚度为1.2mm，距离因子50%；TR = 2300ms，TE = 2.94ms；GRAPPA因子= 2；扫描时间= 5:17。4.1.2 DICOM文件的读取我们需要对2组DICOM序列图像进行处理，第一步为读入序列DCIOM文件。用户层面：ImageFileReader或ImageSeriesReader(读）ImageFileWriter或ImageSeriesWriter（写）内部实现：由内部ImageIO对象具体负责图像文件读写操作，该对象通过对象工厂根据用户输入文件类型生成相应的ImageIO对象优点：a.使用方便，用户无须关注内部实现细节b.扩展方便，扩展支持新的图像读取而无须修改用户接口，只需添加相应的工厂类和IO类对于序列图像的读取，本文使用ImageSeriesReader来读取序列图像，我们仅仅需要通过压缩了与GDCM类交换信息的ITK类调用这些功能。为方便起见，我们仅仅需要把 DICOM 切 片 存 储 的 目 录 名 传 递 给 这 个 类 。 这 些 可 以 交 给 SetDirectory( ) 。GDCMSeriesFileNames将会搜索这个目录并生成DICOM文件名的序列，达到读取序列图像的功能。4.1.3 序列DICOM文件转换为3D体文件因为给定的序列为沿着垂直平行于X轴方向进行的切片，而3D配准需要输入数据为体数据，因此需要首先将DICOM图像序列重建为3D体数据文件。由于ImageIO工厂机构，我们现在保存测定体积的图像到另一个由用户指定的文件中，仅仅需要通过文件扩展名识别文件格式。typedef itk:ImageFileWriter WriterType;WriterType:Pointer writer = WriterType:New( );writer-SetFileName( argv2 );writer-SetInput( reader-GetOutput( ) );写入图像的程序通过调用writer的Upda