[发明专利]一种基于PCA及微分同胚Demons的变形医学图像配准方法

权利要求书

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=uf+um=(m-f)[▿f|▿f|2+α2(f-m)2+▿m|▿m|2+α2(f-m)2]---(1)]]>

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

步骤(3.2.4)平滑变形场

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

un+1=Gδ⊗(un+(m-f)▿f|▿f|2+α2(m-f)2+(m-f)▿m|▿m|2+α2(m-f)2)---(2)]]>

其中，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是统计学常用的相似性测度，前两者的计算公式如下：

ρX,Y=E(XY)-E(X)E(Y)E(X2)-(E(X))2E(Y2)-(E(Y))2∈[-1,1]---(3)]]>

ρX,Y=Σi(Xi-X‾)(Yi-Y‾)Σi(Xi-X‾)2Σi(Yi-Y‾)2∈[-1,1]---(4)]]>

其中，ρ_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。

2.该方法首先通过仿射变换预配准图像，并将图像灰度值归一化到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]，为便于算法的比较，将上述后三者的值依据下列公式进行转换，得到最终的测度值Dpearson Dspearman 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组数据依次实施上述五种配准方法；计算每对数据在每次迭代时的收敛幅值并进行归一化，求十组数据在每次迭代的归一化平均值。