CN104091337B - 一种基于PCA及微分同胚Demons的变形医学图像配准方法 - Google Patents

Abstract

本发明涉及一种基于PCA及微分同胚Demons算法的变形医学图像配准方法，属于医学图像处理技术领域。该方法首先对采集的图像进行预处理。之后采用多分辨率机制，在设定迭代次数下进行配准循环：计算移动图像的变形场；用高斯滤波器进行平滑；利用求得的变形场矩阵，在指数域对移动图像进行插值；使用PCA提取图像关键特征，结合SSD，Pearson，Spearman及Kendall计算两图像间的相似性测度。与传统的demons方法相比，PCA微分同胚log demons相关配准方法，降低了计算量，能很好的抑制噪声，比单纯使用SSD相似性测度时鲁棒性高，提高了收敛速度。

Description

一种基于PCA及微分同胚Demons的变形医学图像配准方法

技术领域

本发明涉及一种基于PCA及微分同胚Demons算法的变形医学图像配准方法，属于医学图像处理技术领域。

背景技术

计算机辅助手术是依靠图像引导的介入性手术，在术前利用当今医学领域的先进成像设备如CT、MRI、PET等,得到病灶的多模医学图像数据,并由医生制定合理的手术方案,在术中利用计算机和立体定位系统进行图像相关处理及实时监视,利用一定的引导系统,对手术进行干预。影响手术导航系统性能的最主要指标是手术引导的精确度。图像数据、立体定位技术、配准算法及影像漂移等是影响精确度的主要因素。其中医学图像配准是整个导航系统的核心内容,也是研究的重点和难点。配准的目的就是将术前规划数据(病人图像或生理结构模型)与术中实时数据统一到同一坐标系下，达到指导手术的目的。如今，术前数据通常是3D的CT，MRI或2D的X线投影透视图像。术中数据通常是2D的超声图像，X线投影透视图像，CT透视图像以及光学图像或3D的锥形CT和超声图像或3D数字化的点或面集。因此依据成像的维数，配准分为3D/2D、3D/3D或2D/3D。不同维数图像需先转换为同一维度下方可继续进行配准。在此重点讨论同一维度下的图像配准。

算法的精准度及速度是医学图像配准最重要的两个评价指标。配准旨在寻求一个变换T，使得在不同时间、场景、设备等条件下获取的同一或不同场景下的两张图像达到最大程度的相近，即使两图像间误差E达到极小。由于变换T有多个参数，因此配准过程中多参数的优化，变换T及相似性测度E是研究配准的关键问题。配准依据变换T分为刚性及非刚性配准，刚性配准适用于骨骼等不易变形的部位，非刚性配准适用于软组织等易变形的器官。非刚性配准的变换T分为基于物理模型，插值原理，知识形变模型和特征约束的模型。其中基于物理模型的配准包括弹性体模型，粘性流体模型，扩散模型，微分同胚模型和基于曲率的配准。弹性配准模型仅能配准小形变的图像，但无法精确配准具有形变较大的图像；粘性流体配准模型可以配准大形变图像，但可能误配准相连的不同组织结构，易得到病态结果。基于插值原理的配准算法具有易于求解的优点，但是这类算法不能揭示组织变形的物理意义，需要人工干预配准过程，配准精度易受特征提取的精度会影响。

扩散模型中的Demons算法是一类优秀的算法，它运算简单，效率高且配准结果较好。Demons是热力学理论中的概念，于19世纪由Maxwell提出。Demons具有这样一种属性，他能够区分分子类别并促使其产生选择性扩散而最终将其隔离。假设A、B两个容器中存放有AB两类混合粒子。A、B容器中间有一个半透膜，上面布满了Demons粒子，由于Demons粒子具有某种特性能够识别A、B粒子，它只允许A分子扩散到A区，B分子扩散到B区,最后将两种粒子完全分离。Thirion把以位移场表示几何变换的非参数配准看成扩散过程，于1998年提出基于光流场的最初的Demons算法，且Hellier于2003年对常用的6种配准算法进行对比，结果表明原始Demons算法在各种精确度评价指标下都最高。Wang于2005年提出ActiveDemons方法，将来自两幅图像绝对梯度的两个单向力进行叠加变为双向力，从而提高了配准精确度和收敛速度。Roglj于2006年提出Symmetric Demons方法。Vercauteren于2007年将微分同胚的概念引入log-Demons配准中，很好的保留了解剖体的拓扑结构。Peyrat于2010年将多通道的Demons用来配准4D时间序列的心脏图像。Lombaert于2014年提出Spectral-log-demons方法，此算法在图像存在大的变形情况下仍可以获得精准的配准结果。本发明采用微分同胚Demons框架。

PCA(Principal Component Aanlysis)即主成分分析，是图像处理中常用的降维方法。假设一张图像有300个特征点，每个特征点有相应的128维向量描述，则一幅图像可用300*128维向量表示，而经过PCA处理，可降为300*64维。而医学图像中3D图像大小可达217*181*181，用PCA降维后可大大较少运算量且提高收敛速度，同时由于PCA提取了图像中关键像素，使得配准结果仍能保持原有配准精度。本发明将PCA和微分同胚log demons算法相结合，创新性如下：1)使用PCA提取图像的关键像素，保留图像最主要的特征；2)在3D图像中，从x，y，z三个维度的每个切片分别应用PCA，保留了图像像素的空间信息；3)将PCA运算结果与Pearson，Spearman等统计学相关系数相结合组成新的相似性测度；4)将此新的相似性测度与微分同胚log Demons配准框架相结合，在保证相同配准质量的条件下提升了收敛速度。本发明在不同3D及2D MRI脑部图像添加及未添加噪声条件下进行测试，实验结果表明PCA相关方法优于原始微分同胚log demons配准算法。

发明内容

本发明的目的在于，针对传统变形配准方法基于手动标记的特征不能适用于所有图像的缺陷，提出一种基于PCA及微分同胚log demons的方法，并将提取的特征作为SSD，Pearson等相似性测度的输入，实验结果表明PCA-Demons方法鲁棒性更高，并提高了收敛速度。

为实现上述目的，本发明采用的技术方案为一种基于PCA及微分同胚Demons的变形医学图像配准方法，该方法的实现步骤如下，

步骤(1)采集数据

分别采集60个人的三维脑部MRI图像各一张，选取其中一人的图像作为固定参考图像，另外一人的图像作为移动待配准图像。另外采集T1/T2时间下的二维脑部图像各两张。T1/T2是MRI的两类加权成像。T1/T2表示处于偏离平衡状态的氢核，在作用力停止后，自动向平衡状态恢复经历的不同时间。T1为纵向弛豫时间，T2为横向弛豫时间。其中三维图像采用UCLA-LONI实验室及英国帝国理工学院的生物医学图像分析组在网上提供的LPBA40及IXI图像集，二维图像采用Matlab Central在网上提供的T1/T2图像。

步骤(2)数据预处理

针对上述采集的图像做如下处理，

步骤(2.1),对固定及移动图像实施仿射变换初配准；

步骤(2.2),将两幅图像的像素灰度值归一化到0～255之间；

步骤(2.3),采用多分辨率机制，由粗到精，设定级数，下采样图像；

步骤(3)开始配准循环

步骤(3.1)初始化参数，并设定循环次数。

步骤(3.2)计算变形场

步骤(3.2.1)将初始值为零的变形场变换到指数空间；

步骤(3.2.2)依据上一步得到的初始变形场线性插值移动图像；

步骤(3.2.3)计算变形场

Demons算法应用在图像配准中，设M为移动待配准图像，F为固定参考图像。把参考图像全部像素点看成Demons点，移动图像视为可变形的网格。每个网格上的Demons力沿着参考图像的灰度梯度方向使浮动图像想参考图像变形，同时引入移动图像的梯度，直到两图像匹配，公式如下：

其中，u为上述的Demons力即待求的变形场，u_f为固定图像对变形场的贡献，下标f代表固定图像分量，u_m为移动图像对变形场的贡献，下标m代表移动图像分量，f为某一点p处的固定图像灰度值，m为移动图像对应p点像素值，▽代表求梯度，α为归一化因子。得到x、y、z三个方向的变形场分量。

步骤(3.2.4)平滑变形场

由于Demons算法利用局部图像信息来变换图像，为保证变换在全局范围连续且保持图像的拓扑结构，使用高斯滤波器平滑得到的偏移向量，公式如下：

其中，u_n+1为第n+1次迭代时的变形场，G_δ为高斯滤波器，下标δ代表滤波器核函数的均方差，为卷积操作，u_n为第n次迭代时的变形场，其余参数的物理意义参照公式(1)。

步骤(3.2.5)将求得的变形场变换到指数域，并利用变形场对移动待配准图像进行线性插值，得到Mp，即配准后的移动图像。

步骤(3.3)关键特征提取

在降维的同时，为了保持图像空间位置信息，分别对x、y、z三个维度的每个切片进行PCA降维，在三个维度上将所有切片得到的降维后的矩阵累加求和，保留了最重要的特征，同时抑制了噪声等干扰。

步骤(3.4)计算相似性测度并优化配准结果

评判两个序列相关系数的相似性测度种类有很多，不同相似性测度体现的要对比的两个序列间的关系不同。常用的有互相关，互信息，归一化互相关、归一化互信息、模式灰度；在此使用偏差平方和(Sum of Squared Deviations,SSD)，Pearson(Pearson product-moment correlation coefficient，PPMCC)，Spearman及Kendall作为相似性测度，后三者的取值范围为[-1,1]。Pearson，Spearman及Kendall是统计学常用的相似性测度，前两者的计算公式如下：

其中，ρ_X,Y为Pearson或Spearman的相似性测度值，E(X)、E(Y)、E(X²)、E(X²)分别为序列X、Y、X²、Y²的期望，X_i、Y_i为X、Y序列中第i个值，为序列X、Y的均值。

Spearman对变量的分布没有要求。而Pearson要求变量呈正态分布，在此在微分同胚Demons方法中引入log，避免了这一影响。将PCA对固定图像F及变换后的移动图像Mp分别降维后得到的矩阵作为输入，分别与SSD，Pearson，Spearman，Kendall组合成四种新的相似性测度，分别标记为PCA-SSD，PCA-Pearson，PCA-Spearman，PCA-Kendall；而Mp与F图像作为输入与SSD作为最原始的相似性测度，标记为Ori-SSD。

步骤(3.4.1)计算PCA-SSD，PCA-Pearson，PCA-Spearman，PCA-Kendall

依据上述方法得到降维后的矩阵，对矩阵的x、y、z三个维度分别应用SSD，Pearson，Spearman，Kendall，求和后取平均。由于Pearson，Spearman，Kendall所得两序列间关系的取值范围为[-1,1]，为便于比较，利用下列公式进行转换，便得到四类相似性测度的值：

Dpearson＝(1-Cpearson)*100 (5)

其中Cpearson为最终变换后的相似性测度值，Dpearson为直接利用测度公式求得的相似性测度值。

步骤(3.4.2)计算Ori-SSD

利用Mp及F计算SSD，取平均值，得到Ori-SSD。

步骤(3.4.3)最优配准判定

判断相似性测度值是否小于提前设定的标准值，如果是则将此测度值赋予标准值。判断此时相似性测度值是否满足中断条件，如果是则结束循环，如果否则继续执行下面的步骤。

步骤(3.4.4)实时显示配准结果

实时显示配准的过程和结果，显示的图像包括，原始固定参考图像F、移动待配准图像M，配准后的移动图像Mp，固定图像与配准后移动图像的差F-Mp，每一次迭代下相似性测度值曲线energy，计算的变形域(ux,uy,uz)以及指数域的变换(sx,sy,sz)。

在预定循环次数及条件内，循环执行步骤(2.3)及步骤3。

与现有技术相比，本发明具有如下有益效果。

本发明提出了一种新的基于PCA微分同胚Demons的变形医学图像配准算法。该算法将PCA降维思想引入到变形医学图像配准中，并与SSD，Pearson，Spearman等组成新的相似性测度。在微分同胚logdemons配准框架下，与原始的配准方法相比，在保证原有配准精度的条件下，表现出较高的鲁棒性，提高了收敛速度。

附图说明

图1为整个算法流程图；

图2为PCA降维算法流程图；

图3a为LPBA40数据集图像配准过程中各相似性测度的收敛幅值；

图3b为LPBA40数据集图像加入高斯白噪声后配准过程中各相似性测度的收敛幅值；

图4a为IXI数据集图像配准过程中各相似性测度的收敛幅值；

图4b为IXI数据集图像加入高斯白噪声后配准过程中各相似性测度的收敛幅值；

图5a为T1数据集中固定参考图像；

图5b为T1数据集中移动待配准图像；

图5c为T1数据集中配准后的移动图像；

图5d为T1数据集中固定参考图像与配准后的移动图像间的差值图像；

图5e为T1数据集图像配准过程中各相似性测度的收敛幅值；

图6a为T2数据集中固定参考图像；

图6b为T2数据集中移动待配准图像；

图6c为T2数据集中配准后的移动图像；

图6d为T2数据集中固定参考图像与配准后的移动图像间的差值图像；

图6e为T2数据集图像配准过程中各相似性测度的收敛幅值。

具体实施方式

本发明采取的技术方案为：

一种基于PCA微分同胚Demons的变形医学图像配准方法。该方法首先通过仿射变换预配准图像，并将图像灰度值归一化到0～255之间。循环配准过程中，利用微分同胚logdemons方法计算得到变形场，并用高斯滤波器进行平滑。将固定参考图像F与配准后的移动图像Mp进行PCA降维，保留关键特征的同时保证了图像空间位置信息的一致。将降维得到的矩阵输入到SSD，Pearson，Spearman，Kendall计算两序列间的相似度。分别使用原始图像及加入高斯噪声的3D及2D图像序列进行测试。实验结果表明，在保证原有配准精度的条件下，PCA微分同胚Demons相关配准方法比未使用PCA降维的配准方法表现出更高的鲁棒性，提高了收敛速度。

本发明具体步骤如下：

步骤(1)采集数据

采用UCLA-LONI实验室提供的LPBA40三维脑MRI数据集，英国帝国理工学院的生物医学图像分析组网上提供的IXI三维脑部MRI图像集，以及Matlab Central提供的T1/T2二维脑图像作为测试数据。其中LPBA40数据集包含40个研究对象。将对象1作为固定参考图像，其余作为移动待配准图像，随机选取十个作为实验对象，形成十对测试图像对，图像大小为217×181×181。T1/T2图像各包含一张固定参考图像和一张待配准的移动图像，形成两对测试图像对，图像大小为192×192。

本发明的目的是提出一种鲁棒性高的变形医学图像配准方法。在保证原有配准精度的同时，提高了算法收敛速度。

步骤(2)数据预处理

步骤(2.1),初配准

为便于后续配准的进行，将上述数据集中的每一对图像实施仿射变换，完成初始配准；

步骤(2.2),灰度值归一化

医学图像的灰度级数高，先获得灰度图的最大及最小值，计算图像的灰度跨度，再将两幅图像内的所有像素归一化到0～255之间；

步骤(2.3),下采样图像

采用多分辨率机制，利用图像在不同层次上的相似性，可使配准精度从低到高逐步提升。设定最大级数为3，对移动和固定图像进行下采样，下采样频率为2^-(N-1),N＝1,2,3。

步骤(3)开始配准循环

步骤(3.1)初始化参数，并设定循环次数。

定义一个结构体参数opt，并初始化参数,简单列举一下重要参数，如下所示：

opt.sigma_diffusion＝1.0；％高斯滤波均方差

opt.sigma_i＝1.0；opt.sigma_x＝1.0；％计算变形场公式中的系数

opt.niter＝250；％配准最大循环迭代次数

opt.vx＝zeros(size(M))；opt.vx＝zeros(size(M))；opt.vx＝zeros(size(M))；％变形场

步骤(3.2)计算变形场

步骤(3.2.1)将初始值为零的变形场变换到指数空间，输入为vx，vy，vz，输出为sx，sy，sz；

步骤(3.2.2)依据上一步得到的初始变形场sx，sy，sz线性插值移动图像M，输出M_prime；

步骤(3.2.3)计算变形场

首先依据公式diff＝F-M_prime计算两图像插值及两图像梯度矩阵[gx,gy,gz][gx_f,gy_f,gz_f]，之后依据变形场公式计算变形场ux，uy，uz。

步骤(3.2.4)平滑变形场

采用三维高斯滤波器，核函数均方差为opt.sigma_fluid＝1，在x，y，z三个维度计算高斯滤波器方差，范围为[-3:3,-3:3,-3:3]。对计算的变形场ux，uy，uz进行高斯平滑滤波。

步骤(3.2.5)在初始变形场vx，vy，vz基础上累积变形场ux，uy，uz，并将其变换到指数域得到sx，sy，sz，利用变形场对移动待配准图像进行线性插值，得到移位后的移动图像Mp。

步骤(3.3)关键特征提取

在x轴维度方向上，依次对每一切片图像实施二维PCA降维操作得到矩阵pca，并将pca累加到x维度PCA分量矩阵pcax上，对y轴及z轴维度切片实施相同的操作，得到pcay，pcaz。

步骤(3.4)计算相似性测度并优化配准结果

本发明比较了Ori-SSD，PCA-SSD，PCA-Pearson，PCA-Spearman，PCA-Kendall五种相似性测度下算法的性能和结果。

步骤(3.4.1)计算PCA-SSD，PCA-Pearson，PCA-Spearman，PCA-Kendall

将移动及固定图像的降维矩阵[mx,my,mz][fx,fy,fz]作为四种相似性测度的输入，依据如下公式计算四种相似性测度值CSSD,Cpearson,Cspearman,Ckendall。

CSSD＝(sum(diff2x(:))+sum(diff2y(:))+sum(diff2z(:)))/area；

Cpearson＝(PEARSON(fx,mx)+PEARSON(fy,my)+PEARSON(fz,mz))/3；

Cspearman＝(SPEARMAN(fx,mx)+SPEARMAN(fy,my)+SPEARMAN(fz,mz))/3；

Ckendall＝(KENDALL(fx,mx)+KENDALL(fy,my)+KENDALL(fz,mz))/3；

由于Pearson，Spearman，Kendall所得相似性测度取值范围限定在[-1,1]，为便于算法的比较，将上述后三者的值依据下列公式进行转换，得到最终的测度值DpearsonDspearman Dkendall。

Dpearson＝(1-Cpearson)*100；

Dspearman＝(1-Cspearman)*100；

Dkendall＝(1-Ckendall)*100；

步骤(3.4.2)计算Ori-SSD

计算Mp及F的差方和，再取均值，得到Ori-SSD。

步骤(3.4.3)最优配准判定

若计算的相似性测度值e(iter)小于提前设定的标准值e_min，用此相似性测度值更新标准值，其中iter代表迭代次数。如果此时相似性测度值满足中断条件则结束循环，否则继续执行下面的步骤。中断条件如下表示，其中opt.stop_criterium＝1e-4。

iter>1&&abs(e(iter)-e(max(1,iter-5)))<e(1)*opt.stop_criterium

步骤(3.4.4)实时显示配准结果

从上到下，从左至右依次实时显示原始固定参考图像F、移动待配准图像M，配准后的移动图像Mp，固定图像与配准后移动图像的差F-Mp，每一次迭代下相似性测度值曲线energy，计算的变形域(ux,uy,uz)以及指数域的变换(sx,sy,sz)。

在满足预定循环次数及设定条件内，循环执行步骤(2.3)及步骤3。

本发明使用不同分辨率、维度、添加及不添加噪声的图像，比较不同相似性测度在不同迭代次数下的值及其变化趋势，分析出适合MRI脑部图像处理的相似性测度。

为了检验本发明所提出的方法的性能，分别在3D的LPBA40，IXI MRI脑部数据集及2D的T1/T2图像上进行了实验。其中选取LPBA40图像十组，IXI图像8组，T1/T2图像各一组。将利用PCA降维下的各相似性测度结果PCA-SSD，PCA-Pearson，PCA-Spearman，PCA-Kendall与原始Ori-SSD相比较。结果表明，在保证配准精度不变的条件下，PCA微分同胚Demons相关方法表现出较快的收敛速度。

对于3D LPBA40数据集，选取11个对象的数据，将第一个作为固定图像，其余10个作为参考图像，组成10组测试数据。图像进行预处理后，对10组数据依次实施上述五种配准方法。计算每对数据在每次迭代时的收敛幅值并进行归一化，求十组数据在每次迭代的归一化平均值。附录图3展示了实验结果，横坐标代表迭代次数，纵坐标代表归一化后的收敛幅值。其中图a)，b)分别展示了未添加及添加噪声的情况下的实验结果。未添加噪声的情况下，Ori-SSD方法表现出较好的性能。但在测试图像中加入高斯噪声后，原始的Ori-SSD方法表现出了相对噪声鲁棒性较差的特性，本应逐渐收敛的曲线变得逐渐发散，相比之下，PCA-SSD和PCA-Pearson方法表现出了较好的收敛效果。另外，Kendall、Spearman方法也表现出了鲁棒性较差的特点。

对于3D IXI数据集，选取9个对象的数据，组成8组测试数据，附录图4展示了实验结果，可得出类似LPBA40数据的结论。

T1/T2代表了不同弛豫时间下的MRI成像。附录图5展示了2D T1图像对的实验结果。其中图a)至图e)依次表示固定参考图像，移动待配准图像，配准后的移动图像，配准后移动图像与固定图像间的差图像，收敛幅值图。PCA-SSD方法和PCA-Pearson方法比原始的Ori-SSD方法收敛性好。PCA-Kendall方法的收敛性并不均匀，但相比原始Ori-SSD方法，整体性能上还是比其要好。

2DT2图像配准结果见附录图6，可以得到和T1图像类似的结果，但是Kendall方法明显不如Ori-SSD方法收敛性好。

实验结果表明，将PCA应用到微分同胚Demons配准框架中，并与SSD，Pearson等组成新的相似性测度。在处理MRI图像时，PCA微分同胚Demons算法提取图像最主要的特征，降低了计算量；同时能很好的抑制噪声，比单纯使用SSD相似性测度时鲁棒性高；提高了收敛速度。

Claims (2)

1.一种基于PCA及微分同胚Demons的变形医学图像配准方法，其特征在于：该方法的实现步骤如下，

步骤(1)采集数据

分别采集60个人的三维脑部MRI图像各一张，选取其中一人的图像作为固定参考图像，另外一人的图像作为移动待配准图像；另外采集T1/T2时间下的二维脑部图像各两张；T1/T2是MRI的两类加权成像；T1/T2表示处于偏离平衡状态的氢核，在作用力停止后，自动向平衡状态恢复经历的不同时间；T1为纵向弛豫时间，T2为横向弛豫时间；其中三维图像采用UCLA-LONI实验室及英国帝国理工学院的生物医学图像分析组在网上提供的LPBA40及IXI图像集，二维图像采用Matlab Central在网上提供的T1/T2图像；

步骤(2)数据预处理

针对上述采集的图像做如下处理，

步骤(2.1),对固定及移动图像实施仿射变换初配准；

步骤(2.2),将两幅图像的像素灰度值归一化到0～255之间；

步骤(2.3),采用多分辨率机制，由粗到精，设定级数，下采样图像；

步骤(3)开始配准循环

步骤(3.1)初始化参数，并设定循环次数；

步骤(3.2)计算变形场

步骤(3.2.1)将初始值为零的变形场变换到指数空间；

步骤(3.2.2)依据上一步得到的初始变形场线性插值移动图像；

步骤(3.2.3)计算变形场

Demons算法应用在图像配准中，设M为移动待配准图像，F为固定参考图像；把参考图像全部像素点看成Demons点，移动图像视为可变形的网格；每个网格上的Demons力沿着参考图像的灰度梯度方向使浮动图像向参考图像变形，同时引入移动图像的梯度，直到两图像匹配，公式如下：

$u=u_f+u_m=(m-f)\&lsqb;\frac{\&dtri;f}{|\&dtri;f{|}^2+{\mathrm{\&alpha;}}^2{(f-m)}^2}+\frac{\&dtri;m}{|\&dtri;m{|}^2+{\mathrm{\&alpha;}}^2{(f-m)}^2}\&rsqb;---\left(1\right)$

其中，u为上述的Demons力即待求的变形场，u_f为固定图像对变形场的贡献，下标f代表固定图像分量，u_m为移动图像对变形场的贡献，下标m代表移动图像分量，f为某一点p处的固定图像灰度值，m为移动图像对应p点像素值，代表求梯度，α为归一化因子；得到x、y、z三个方向的变形场分量；

步骤(3.2.4)平滑变形场

由于Demons算法利用局部图像信息来变换图像，为保证变换在全局范围连续且保持图像的拓扑结构，使用高斯滤波器平滑得到的偏移向量，公式如下：

$u_{n+1}=G_{\mathrm{\&delta;}}\&CircleTimes;(u_n+\frac{(m-f)\&dtri;f}{|\&dtri;f{|}^2+{\mathrm{\&alpha;}}^2{(m-f)}^2}+\frac{(m-f)\&dtri;m}{|\&dtri;m{|}^2+{\mathrm{\&alpha;}}^2{(m-f)}^2})---\left(2\right)$

其中，u_n+1为第n+1次迭代时的变形场，G_δ为高斯滤波器，下标δ代表滤波器核函数的均方差，为卷积操作，u_n为第n次迭代时的变形场，其余参数的物理意义参照公式(1)；

步骤(3.2.5)将求得的变形场变换到指数域，并利用变形场对移动待配准图像进行线性插值，得到Mp，即配准后的移动图像；

步骤(3.3)关键特征提取

在降维的同时，为了保持图像空间位置信息，分别对x、y、z三个维度的每个切片进行PCA降维，在三个维度上将所有切片得到的降维后的矩阵累加求和，保留了最重要的特征，同时抑制了噪声干扰；

步骤(3.4)计算相似性测度以及图像显示

评判两个序列相关系数的相似性测度种类有很多，不同相似性测度体现的要对比的两个序列间的关系不同；常用的有互相关，互信息，归一化互相关、归一化互信息、模式灰度；在此使用偏差平方和，Pearson，Spearman及Kendall作为相似性测度，后三者的取值范围为[-1,1]；Pearson，Spearman及Kendall是统计学常用的相似性测度，前两者的计算公式如下：

${\mathrm{\&rho;}}_{X,Y}=\frac{E\left(XY\right)-E\left(X\right)E\left(Y\right)}{\sqrt{E\left(X^2\right)-{\left(E\right(X\left)\right)}^2}\sqrt{E\left(Y^2\right)-{\left(E\right(Y\left)\right)}^2}}\&Element;\&lsqb;-1,1\&rsqb;---\left(3\right)$

${\mathrm{\&rho;}}_{X,Y}=\frac{\underset{i}{\&Sigma;}(X_i-\stackrel{\&OverBar;}{X})(Y_i-\stackrel{\&OverBar;}{Y})}{\sqrt{\underset{i}{\&Sigma;}{(X_i-\stackrel{\&OverBar;}{X})}^2\underset{i}{\&Sigma;}{(Y_i-\stackrel{\&OverBar;}{Y})}^2}}\&Element;\&lsqb;-1,1\&rsqb;---\left(4\right)$

其中，ρ_X,Y为Pearson或Spearman的相似性测度值，E(X)、E(Y)、E(X²)、E(Y²)分别为序列X、Y、X²、Y²的期望，X_i、Y_i为X、Y序列中第i个值，为序列X、Y的均值；

Spearman对变量的分布没有要求；而Pearson要求变量呈正态分布，在此在微分同胚Demons方法中引入log，避免了这一影响；将PCA对固定图像F及变换后的移动图像Mp分别降维后得到的矩阵作为输入，分别与SSD，Pearson，Spearman，Kendall组合成四种新的相似性测度，分别标记为PCA-SSD，PCA-Pearson，PCA-Spearman，PCA-Kendall；而Mp与F图像作为输入与SSD作为最原始的相似性测度，标记为Ori-SSD；

步骤(3.4.1)计算PCA-SSD，PCA-Pearson，PCA-Spearman，PCA-Kendall

依据上述方法得到降维后的矩阵，对矩阵的x、y、z三个维度分别应用SSD，Pearson，Spearman，Kendall，求和后取平均；由于Pearson，Spearman，Kendall所得平均后的两序列间关系的取值范围为[-1,1]，为便于比较，利用下列公式进行转换，便得到三类相似性测度的值：

Dpearson＝(1-Cpearson)*100 (5)

Dspearman＝(1-Cspearman)*100；

Dkendall＝(1-Ckendall)*100；

其中Dpearson为最终变换后的相似性测度值，Cpearson为直接利用测度公式求得的相似性测度值；

步骤(3.4.2)计算Ori-SSD

利用Mp及F计算SSD，取平均值，得到Ori-SSD；

步骤(3.4.3)最优配准判定

判断相似性测度值是否小于标准值，如果是则将此测度值赋予标准值；判断此时相似性测度值是否满足中断条件，如果是则结束循环，如果否则继续执行下面的步骤；

步骤(3.4.4)实时显示配准结果

实时显示配准的过程和结果，显示的图像包括，原始固定参考图像F、移动待配准图像M，配准后的移动图像Mp，固定图像与配准后移动图像的差F-Mp，每一次迭代下相似性测度值曲线energy，计算的变形域(ux,uy,uz)以及指数域的变换(sx,sy,sz)；

在预定循环次数及条件内，循环执行步骤(2.3)及步骤3。

2.根据权利要求1所述的一种基于PCA及微分同胚Demons的变形医学图像配准方法，其特征在于：该方法首先通过仿射变换预配准图像，并将图像灰度值归一化到0～255之间；循环配准过程中，利用微分同胚log demons方法计算得到变形场，并用高斯滤波器进行平滑；将固定参考图像F与配准后的移动图像Mp进行PCA降维，保留关键特征的同时保证了图像空间位置信息的一致；将降维得到的矩阵输入到SSD，Pearson，Spearman，Kendall计算两序列间的相似度；分别使用原始图像及加入高斯噪声的3D及2D图像序列进行测试；

本方法具体步骤如下：

步骤(1)采集数据

采用UCLA-LONI实验室提供的LPBA40三维脑MRI数据集，英国帝国理工学院的生物医学图像分析组网上提供的IXI三维脑部MRI图像集，以及Matlab Central提供的T1/T2二维脑图像作为测试数据；其中LPBA40数据集包含40个研究对象；将对象1作为固定参考图像，其余作为移动待配准图像，随机选取十个作为实验对象，形成十对测试图像对，图像大小为217×181×181；T1/T2图像各包含一张固定参考图像和一张待配准的移动图像，形成两对测试图像对，图像大小为192×192；

步骤(2)数据预处理

步骤(2.1),初配准

为便于后续配准的进行，将上述数据集中的每一对图像实施仿射变换，完成初始配准；

步骤(2.2),灰度值归一化

医学图像的灰度级数高，先获得灰度图的最大及最小值，计算图像的灰度跨度，再将两幅图像内的所有像素归一化到0～255之间；

步骤(2.3),下采样图像

采用多分辨率机制，利用图像在不同层次上的相似性，可使配准精度从低到高逐步提升；设定最大级数为3，对移动和固定图像进行下采样，下采样频率为2^-(N-1),N＝1,2,3；

步骤(3)开始配准循环

步骤(3.1)初始化参数，并设定循环次数；

定义一个结构体参数opt，并初始化参数,简单列举一下重要参数，如下所示：

opt.sigma_diffusion＝1.0；％高斯滤波均方差

opt.sigma_i＝1.0；opt.sigma_x＝1.0；％计算变形场公式中的系数

opt.niter＝250；％配准最大循环迭代次数

opt.vx＝zeros(size(M))；opt.vx＝zeros(size(M))；opt.vx＝zeros(size(M))；％变形场

步骤(3.2)计算变形场

步骤(3.2.1)将初始值为零的变形场变换到指数空间，输入为vx，vy，vz，输出为sx，sy，sz；

步骤(3.2.2)依据上一步得到的初始变形场sx，sy，sz线性插值移动图像M，输出M_prime；

步骤(3.2.3)计算变形场

首先依据公式diff＝F-M_prime计算两图像差值及两图像梯度矩阵[gx,gy,gz][gx_f,gy_f,gz_f]，之后依据变形场公式计算变形场ux，uy，uz；

步骤(3.2.4)平滑变形场

采用三维高斯滤波器，核函数均方差为opt.sigma_fluid＝1，在x，y，z三个维度计算高斯滤波器方差，范围为[-3:3,-3:3,-3:3]；对计算的变形场ux，uy，uz进行高斯平滑滤波；

步骤(3.2.5)在初始变形场vx，vy，vz基础上累积变形场ux，uy，uz，并将其变换到指数域得到sx，sy，sz，利用变形场对移动待配准图像进行线性插值，得到移位后的移动图像Mp；

步骤(3.3)关键特征提取

在x轴维度方向上，依次对每一切片图像实施二维PCA降维操作得到矩阵pca，并将pca累加到x维度PCA分量矩阵pcax上，对y轴及z轴维度切片实施相同的操作，得到pcay，pcaz；

步骤(3.4)计算相似性测度以及图像显示

比较Ori-SSD，PCA-SSD，PCA-Pearson，PCA-Spearman，PCA-Kendall五种相似性测度下算法的性能和结果；

步骤(3.4.1)计算PCA-SSD，PCA-Pearson，PCA-Spearman，PCA-Kendall将移动及固定图像的降维矩阵[mx,my,mz][fx,fy,fz]作为四种相似性测度的输入，依据如下公式计算四种相似性测度值CSSD,Cpearson,Cspearman,Ckendall；

CSSD＝(sum(diff2x(:))+sum(diff2y(:))+sum(diff2z(:)))/area；

Cpearson＝(PEARSON(fx,mx)+PEARSON(fy,my)+PEARSON(fz,mz))/3；

Cspearman＝(SPEARMAN(fx,mx)+SPEARMAN(fy,my)+SPEARMAN(fz,mz))/3；

Ckendall＝(KENDALL(fx,mx)+KENDALL(fy,my)+KENDALL(fz,mz))/3；

由于Pearson，Spearman，Kendall所得相似性测度取值范围限定在[-1,1]，为便于算法的比较，将上述后三者的值依据下列公式进行转换，得到最终的测度值DpearsonDspearman Dkendall；

Dpearson＝(1-Cpearson)*100；

Dspearman＝(1-Cspearman)*100；

Dkendall＝(1-Ckendall)*100；

步骤(3.4.2)计算Ori-SSD

计算Mp及F的差方和，再取均值，得到Ori-SSD；

步骤(3.4.3)最优配准判定

若计算的相似性测度值e(iter)小于提前设定的标准值e_min，用此相似性测度值更新标准值，其中iter代表迭代次数；如果此时相似性测度值满足中断条件则结束循环，否则继续执行下面的步骤；中断条件如下表示，其中opt.stop_criterium＝1e-4；

iter>1&&abs(e(iter)-e(max(1,iter-5)))<e(1)*opt.stop_criterium

步骤(3.4.4)实时显示配准结果

从上到下，从左至右依次实时显示原始固定参考图像F、移动待配准图像M，配准后的移动图像Mp，固定图像与配准后移动图像的差F-Mp，每一次迭代下相似性测度值曲线energy，计算的变形域(ux,uy,uz)以及指数域的变换(sx,sy,sz)；

在满足预定循环次数及设定条件内，循环执行步骤(2.3)及步骤3；

本方法使用不同分辨率、维度、添加及不添加噪声的图像，比较不同相似性测度在不同迭代次数下的值及其变化趋势，分析出适合MRI脑部图像处理的相似性测度；

为了检验本方法所提出的方法的性能，分别在3D的LPBA40，IXI MRI脑部数据集及2D的T1/T2图像上进行了实验；其中选取LPBA40图像十组，IXI图像8组，T1/T2图像各一组；将利用PCA降维下的各相似性测度结果PCA-SSD，PCA-Pearson，PCA-Spearman，PCA-Kendall与原始Ori-SSD相比较；结果表明，在保证配准精度不变的条件下，PCA微分同胚Demons相关方法表现出较快的收敛速度；

对于3D LPBA40数据集，选取11个对象的数据，将第一个作为固定图像，其余10个作为参考图像，组成10组测试数据；图像进行预处理后，对10组数据依次实施上述五种配准方法；计算每对数据在每次迭代时的收敛幅值并进行归一化，求十组数据在每次迭代的归一化平均值。