基于稀疏低秩与Atlas集的3D MRI胰腺分割方法

摘要

本发明提出一种基于稀疏低秩与Atlas集的3D MRI胰腺分割方法，主要解决现有技术中分割胰腺粘连现象严重、分割结果不精确的问题。其步骤如下：1)输入3D MRI数据并归一化处理；2)对归一化处理后的序列图像进行稀疏低秩分解并进行Hessian矩阵增强；3)在增强后的图像上人机交互分割出胰腺标签；4)构造粗分割集；5)输入Atlas集，依次进行刚性配准和弹性配准，并进行标签融合，输出分割结果。本发明运用矩阵增强、矩阵分解、图集法，并利用3D MRI图像的空间信息，进行胰腺分割，减少了胰腺目标与其他组织器官的粘连现象，图像轮廓清晰，提高了分割的准确度，可用于胰腺癌精确放射治疗中的胰腺目标检测。

说明书

技术领域

本发明属于医学图像处理技术领域，更进一步涉及三维核磁共振图像3D MRI胰腺分割，可用于胰腺癌精确放射治疗中的胰腺目标检测。

背景技术

胰腺癌是出现在胰脏的癌症，通常认为胰腺癌是常见肿瘤中恶性程度最高、死亡率最高的癌症，约90％的患者无法以手术根治，五年存活率低于5％。根据世界卫生组织估计，胰腺癌在全世界发生率排第十三，死亡率排第八。随着科技的发展，其他癌症的治疗都有了提高，但胰腺癌却没有明显的突破，死亡率也没有明显的下降趋势，所以必须加强研究，早日攻克胰腺癌对人类的威胁。

传统检查胰腺癌的方式是医生根据经验和观察直接阅片，定性的判断患者的病情，这样易受情绪影响，且随着工作量增加，阅片准确率下降。对胰腺癌的治疗方法中，放射治疗是一种主要的方法，他试图用射线杀死肿瘤，但也会对正常组织造成一次性或永久性伤害。精确放疗通过精确的肿瘤定位、精确的计划设计、精确的剂量计算及在治疗机上精确执行的一种全新的肿瘤放疗技术，有助于减少放射治疗中对健康组织放射造成的伤害。精确分割胰腺组织有助于帮助医生分析胰腺病情，也是精确放疗不可缺少的一步。

目前，对胰腺癌进行检测和治疗的设备主要是核磁共振成像MRI系统，该设备功能强大、技术含量高、软组织图像清晰，而且对人体检测损害小，已成为对胰腺进行分析、诊断和治疗的重要手段，而从MRI图像中精确的分割出胰腺，就是这些手段的重要前提。

传统医学图像分割方法包括基于区域的方法、基于边缘的方法、基于神经网络的方法、基于模糊聚类的方法、基于图谱的方法等。由于这些方法大多数是基于2D图像序列的分割，利用2D图像本身的特性信息完成感兴趣区域的分割，例如灰度、纹理等信息，而没有利用体素的空间信息，因此不能进行准确的胰腺提取；医学MRI图像有噪声强、对比度低、软组织边界不明显的特点，而这些方法大多数对图像噪声和对比度敏感，所以难以定位目标边界，使得分割结果不精确。

发明内容

本发明的目的在于针对以上现有技术的不足，提出一种基于稀疏低秩与Atlas集的3D MRI胰腺分割方法，以提高3D MRI胰腺分割的精度。

实现本发明目的的技术方案如下：

(1)输入三维核磁共振3D MRI序列图像，对其进行归一化处理，得到归一化后的序列图像C；

(2)用交替迭代法对归一化处理后的序列图像C进行稀疏低秩矩阵分解，得到稀疏图像序列A和低秩图像序列B；

(3)对稀疏图像序列A进行Hessian矩阵增强，得到二值增强图像序列A₁；

(4)通过一次人机交互，在二值增强图像序列A₁中确定胰腺的位置区域P，根据胰腺位置区域P，利用形态学原理在二值增强图像序列A₁中分割出胰腺的标签图像序列L；

(5)根据标签图像序列L和步骤(1)中的序列图像C，得到粗分割集Q，其中包含F帧图像；

(6)输入胰腺图像序列W及人工勾画出的胰腺标签序列E，即Atlas集；

(7)选择粗分割集Q中的第i帧图像作为参考图像，记为Q_i，将Atlas集中的胰腺图像序列W作为浮动图像，用基于互信息的刚性配准算法对Q_i、W进行配准,得到刚性配准后的胰腺图像序列W₁和刚性转换公式T₁，根据T₁对胰腺标签序列E进行迁移，得到刚性配准后的胰腺标签序列E₁，其中i的取值范围为1～F；

(8)将刚性配准后的胰腺图像序列W₁作为浮动图像，使用基于Demons的弹性配准算法对参考图像Q_i和该W₁进行配准，得到弹性配准后的胰腺图像序列W₂和弹性转换公式T₂，根据T₂对刚性配准后的胰腺标签序列E₁进行迁移，得到弹性配准后的胰腺标签序列E₂；

(9)计算参考图像Q_i与弹性配准后的胰腺图像序列W₂的相关性序列，根据该相关性序列及弹性配准后的胰腺图像序列W₂和弹性配准后的胰腺标签序列E₂这三者是一一对应的关系，取相关性序列最高的前五个系数对应的标签，并对这五个标签进行融合，得到参考图像Q_i精分割后的结果R_i，其中i的取值范围为1～F；

(10)重复步骤(7)～(9)，对粗分割集Q中所有图像精分割后，得到最终结果精分割集R。

本发明与现有技术相比具有以下优点：

1.本发明采用Hessian矩阵增强的方法，将三维核磁共振图像3D MRI的空间信息运用到胰腺目标的分割中，使分割结果更加准确。

2.本发明采用稀疏低秩矩阵分解方法，将三维核磁共振图像3D MRI分解成低秩序列部分和稀疏序列部分，通过对稀疏序列部分进行操作，有效减少目标与其他器官的粘连现象。

3.本发明运用Atlas集方法，将粗分割结果与胰腺模板库图像配准，再进行标签迁移和融合得到最终分割结果，使得分割目标区域能保持较好的轮廓。

附图说明

图1是本发明的实现流程图；

图2是待分割3D MRI图像序列中的一帧及该帧的标签；

图3是用本发明与现有level set法和graph cut法的分割结果对比图。

具体实施方式

参照图1，本发明的实现步骤如下：

步骤1、输入数据并归一化处理。

输入三维核磁共振图像3D MRI序列C₀，按下列公式对其进行归一化处理，得到归一化处理后的序列图像C：

$C=\frac{C_0}{max\left(C_0\right)}$

其中，max()为提取矩阵最大元素。

步骤2、对归一化处理后的序列图像进行稀疏低秩分解。

用交替迭代法对归一化处理后的序列图像C进行稀疏低秩矩阵分解，得到稀疏图像序列A和低秩图像序列B:

(2a)设置初始值，稀疏图像序列A为零矩阵，低秩图像序列B＝C；

(2b)通过公式：更新稀疏矩阵A，其中，是欧几里得投影，γ为常数，且γ＞0，β指偏离线性约束的惩罚参数，且β＞0,k为迭代次数，B^k为第k次迭代后的低秩矩阵，Z^k为第k次迭代后的约束矩阵；

(2c)对进行奇异值分解产生三个不同的矩阵：左奇异向量矩阵U^k+1、右奇异向量矩阵V^k+1及对角奇异值矩阵

(2d)通过公式：更新低秩矩阵B,其中，T为矩阵的转置，max{}为提取最大元素，diag()为提取对角线元素；

(2e)通过公式Z^k+1＝Z^k-β(A^k+1+B^k+1-C)更新约束性矩阵Z，其中，A^k+1是第k+1次迭代后的稀疏矩阵，Z^k+1是第k+1次迭代后的线性约束矩阵，B^k+1是第k+1次迭代后的低秩矩阵；

(2f)设定迭代次数k＝1000，重复步骤(2a)～(2e)，迭代结束后，得到稀疏图像序列A和低秩图像序列B。

步骤3、对稀疏图像序列进行Hessian矩阵增强。

对稀疏图像序列A进行Hessian矩阵增强，得到二值增强图像序列A₁:

(3a)根据下面公式，求出稀疏图像序列A中每个像素点I的Hessian矩阵H：

$H=\left(\begin{array}{ccc}{I}_{xx}& I_{xy}& I_{xz}\\ I_{yx}& I_{yy}& I_{yz}\\ I_{zx}& I_{zy}& I_{zz}\end{array}\right)$

其中，I表示稀疏图像序列A中每一个像素点，I_xx代表沿x方向的二阶偏微分，I_xy代表先沿x方向再沿y方向的二阶偏微分，I_xz代表先沿x方向再沿z方向的二阶偏微分，I_yx代表先沿y方向再沿x方向的二阶偏微分，I_yy代表沿y方向的二阶偏微分，I_yz代表先沿y方向再沿z方向的二阶偏微分，I_zx代表先沿z方向再沿x方向的二阶偏微分，I_zy代表先沿z方向再沿y方向的二阶偏微分，I_zz代表沿z方向的二阶偏微分；

(3b)根据Cholesky分解法，计算步骤(3a)中每一个Hessian矩阵H的三个特征值，并按绝对值从小到大顺序排序为λ₁，λ₂，λ₃；

(3c)根据下面公式和每一个像素点I的Hessian矩阵特征值，对稀疏图像序列A的每一个像素点进行增强；

$v_1\left(\lambda \right)=\left(\begin{array}{cc}0& \begin{array}{cc}if& {\lambda }_2<0{\mathrm{or\lambda }}_3<0\end{array}\\ (1-\exp (-\frac{R_A^2}{2{\alpha }^2}\left)\right)\exp (-\frac{R_B^2}{2{\tau }^2})(1-\exp (-\frac{S^2}{2c^2}\left)\right)& \end{array}\right)$

其中，α是起调节面状权重R_A作用的常数，α＞0，τ是起调节管状权重R_B作用的常数，τ＞0，c是起调节球状权重S作用的常数，c＞0，v₁(λ)是指每个像素点I增强后的结果。

步骤4、人机交互分割出胰腺标签。

利用matlab上的roipoly函数，在步骤3中得到的二值增强图像序列A₁上手动选择胰腺的位置区域P，利用形态学原理，将A₁中与P组成连通域的部分提取出来，得到胰腺的标签图像序列L:

(4a)对二值增强图像序列A₁进行形态学腐蚀操作，得到腐蚀后的二值增强图像序列A₁′：

(4b)根据胰腺位置区域P从腐蚀后的二值增强图像序列A₁中提取出连通区域L′；

(4c)对提取出来的连通区域L′进行形态学膨胀操作，得到标签图像序列L。

步骤5、构造粗分割集。

根据标签图像序列L和步骤(1)中的序列图像C，按下列公式，得到包含F帧图像的粗分割集Q：

Q＝C.*L

其中，.*表示矩阵对应元素相乘。

步骤6、输入Atlas集。

输入胰腺图像序列W及人工勾画出的胰腺标签序列E，W和E构成Atlas集。

步骤7、基于互信息的刚性配准。

(7a)选择粗分割集Q中的第i帧图像作为参考图像，记为Q_i，将Atlas集中的胰腺图像序列W作为浮动图像，用基于互信息的刚性配准算法对Q_i、W进行配准,得到刚性配准后的胰腺图像序列W₁和刚性转换公式T₁：

(7a1)设初始状态k＝0，初始参数为(0,0,0)，W₁⁰＝W，其中，p₁代表水平变量、p₂代表垂直变量、p₃代表顺时针旋转变量，k为迭代次数,W₁⁰为迭代第0次后的刚性配准胰腺图像序列；

(7a2)用PV插值法统计Q_i和W₁^k的联合直方图，根据该联合直方图计算Q_i和W₁^k的互信息值，其中，W₁^k为迭代第k次后的刚性配准胰腺图像序列；

(7a3)用Powell算法根据步骤(7a2)中得到的互信息值判断参数是否最优，如果不是最优，则进行步骤(7a4),如果是最优，则停止迭代，输出最优参数进行步骤(7a5)，其中，为迭代第k次后的水平变量，为迭代第k次后的垂直变量，为迭代第k次后的顺时针旋转变量；

(7a4)对W₁^k进行平移、旋转，得到W₁^k+1和调整后的参数令k＝k+1,返回到步骤(7a2)，其中，W₁^k+1为迭代第k+1次后的刚性配准胰腺图像序列，为迭代第k+1次后的水平变量，为迭代第k+1次后的垂直变量，为迭代第k+1次后的顺时针旋转变量；

(7a5)根据最优参数(p₁,p₂,p₃)按下列公式得到刚性转换公式T₁，根据T₁将W转换成刚性配准后的胰腺图像序列W₁：

(7b)根据T₁对胰腺标签序列E进行迁移，得到刚性配准后的胰腺标签序列E₁。

步骤8、基于Demons的非刚性配准。

(8a)将刚性配准后的胰腺图像序列W₁作为浮动图像，使用基于Demons的弹性配准算法对参考图像Q_i和该W₁进行配准，得到弹性配准后的胰腺图像序列W₂和弹性转换公式T₂：

(8a1)设初始状态k＝1，在参考图像Q_i上选择所有像素点作为Demons点，给定允许误差ε，其中ε＞0,k为迭代次数,其中为第一次迭代的弹性配准胰腺图像序列；

(8a2)根据Symmetric Demons公式计算中所有像素点的形变向量：

$u^k=\frac{2(W_2^k-Q_i)(▿W_2^k+▿Q_i)}{|▿W_2^k+▿Q_i{|}^2+{\eta }^2{(W_2^k-Q_i)}^2}$

其中，u表示形变向量，u^k表示迭代第k次后的形变向量，表示迭代第k次后弹性配准胰腺图像序列，表示求图像的梯度，||表示取模操作，η表示归一化因子；

(8a3)根据形变向量u^k按下列公式求出弹性转换公式根据将转换为其中，为迭代第k次后的弹性转换公式，为迭代第k+1次后的弹性配准胰腺图像序列：

$T_2^k=\xi ·u^k$

其中，ξ表示项力衰减系数；

(8a4)计算Q_i与的误差平方和，如果该误差平方和小于ε，输出弹性转换公式和弹性配准后的胰腺图像序列结束迭代，否则k＝k+1,返回步骤(8a2)；

(8b)根据T₂对刚性配准后的胰腺标签序列E₁进行迁移，得到弹性配准后的胰腺标签序列E₂。

步骤9、进行标签融合。

计算参考图像Q_i与弹性配准后的胰腺图像序列W₂的相关性序列，根据该相关性序列及弹性配准后的胰腺图像序列W₂和弹性配准后的胰腺标签序列E₂这三者是一一对应的关系，取相关性序列最高的前五个系数对应的标签，并对这五个标签进行融合，得到参考图像Q_i精分割后的结果R_i。

(9a)根据相关性序列最高的前五个系数g₁～g₅和他们对应的标签图像G₁～G₅，计算标签融合后的结果R_i′：

${R_i}^{\prime }=\underset{j=1}{\overset{5}{\Sigma }}\frac{g_j}{{\Sigma }_{j=1}^5{g}_j}{G}_j$

其中j的取值范围为1～5，∑表示求和操作；

(9b)按下列公式对(9a)中标签融合后的结果R_i′进行修正，得到精分割后的结果R_i：

$R_i=\left(\begin{array}{cc}0& {R_i}^{\prime }<0.4\\ 1& {R_i}^{\prime }\ge 0.4\end{array}\right)$

步骤10、重复步骤(7)～(9)，对粗分割集Q中所有图像精分割，得到最终结果精分割集R。

以下结合仿真图对本发明的效果做进一步的说明。

1.仿真条件：

硬件平台为：CPU为Inter Core i7-4790K，主频为4GHz，内存为32GB；

软件平台为：Windows 10专业版64位操作系统，MatlabR2014a。

2.仿真数据：

本例中待分割的数据采用的是通过三维核磁共振仪获取的一名志愿者腹部三维核磁共振图像序列，图像大小为320*284，一共有5帧图像，通过人工勾画出胰腺作为标签，如图2所示，其中，图2a是胰腺图像，图2b是胰腺图像的标签。Atlas集采用的是其他志愿者通过三维核磁共振仪获取的三维核磁共振图像，它们对应的标签图像通过人工勾画获得。

3.仿真内容：

为了验证算法的有效性，一般在进行图像分割以后，会以标签为参考图像计算指标，豪斯多夫距离、平均绝对表面距离、相似性、灵敏度、特异性五个指标作为分割精度的判断标准,豪斯多夫距离、平均绝对表面距离这两个指标数值越大分割结果越精确，相似性、灵敏度、特异性这三个指标数值越小分割结果越精确。

用本发明和现有的两个具有代表性的图像分割方法level set法和graph cut法进行仿真实验对比。

仿真一：分别通过这三种方法对上述仿真数据进行胰腺分割仿真，对三种方法得到的分割结果进行三维重建，结果如图3所示，其中图3a为本发明的分割结果，图3b为level set法的分割结果，图3c为graph cut法的分割结果，可以看出用本发明分割出的胰腺轮廓更好。

仿真二：分别通过这三种方法对上述仿真数据进行胰腺分割仿真，计算分割结果的五个指标：豪斯多夫距离、平均绝对表面距离、相似性、灵敏度、特异性，结果如表1所示：

表1用本发明与现有level set法和graph cut法的分割指标对比表

豪斯多夫距离平均绝对表面距离相似性灵敏性特异性本发明18.0482mm1.6165mm0.79640.76320.9966Level set32.3503mm3.9333mm0.57210.47260.996Graph cut69.7327mm12.3378mm0.44760.64590.9722

从表1可以看出本发明的分割结果在指标上优于其他两种算法，说明本发明分割胰腺结果精度更高。