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

摘要

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

主权项

一种基于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)[\frac{\▿f}{{|\▿f|}^2+{\mathrm{\α}}^2{(f-m)}^2}+\frac{\▿m}{{|\▿m|}^2+{\mathrm{\α}}^2{(f-m)}^2}]---\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{\δ}}\⊗(u_n+\frac{(m-f)\▿f}{{|\▿f|}^2+{\mathrm{\α}}^2{(m-f)}^2}+\frac{(m-f)\▿m}{{|\▿m|}^2+{\mathrm{\α}}^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{\ρ}}_{X,Y}=\frac{E\left(\mathrm{XY}\right)-E\left(X\right)E\left(Y\right)}{\sqrt{E\left(X^2\right)-{\left(E\left(X\right)\right)}^2}\sqrt{E\left(Y^2\right)-{\left(E\left(Y\right)\right)}^2}}\∈[-\mathrm{1,1}]---\left(3\right)$${\mathrm{\ρ}}_{X,Y}=\frac{\underset{i}{\mathrm{\Σ}}(X_i-\stackrel{\‾}{X})(Y_i-\stackrel{\‾}{Y})}{\sqrt{\underset{i}{\mathrm{\Σ}}{(X_i-\stackrel{\‾}{X})}^2\underset{i}{\mathrm{\Σ}}{(Y_i-\stackrel{\‾}{Y})}^2}}\∈[-\mathrm{1,1}]---\left(4\right)$其中，ρ_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。