CN107330953B - 一种基于非凸低秩的动态mri重建方法 - Google Patents

Abstract

本发明公开一种基于非凸低秩的动态MRI重建方法，其将非凸正则项引入动态MRI的重建模型中，减少与真实值之间的误差，即在分离MR图像时，可以得到更精确的前景和背景。其次在视觉效果上，能够保留较多的纹理信息和轮廓信息，更便于病情的诊断。

Description

一种基于非凸低秩的动态MRI重建方法

技术领域

本发明涉及数字图像处理技术领域，具体涉及一种基于非凸低秩的动态MRI(Magnetic Resonance Imaging，核磁共振成像)重建方法。

背景技术

图像重建在医疗领域的应用越来越广泛，且大部分临床诊断依赖于计算机硬件设备，因此图像重建性能的提高对诊断的精确度具有重大意义。当采用MRI对心脏，胃等器官检测时，由于这些器官是动态的，使得最后呈现的图像受此干扰而产生不同程度的伪影，导致质量下降，影响对病情的诊断。因此，构建一个重建质量较高且快速的方法成为解决该问题的有效途径。

针对这一现状，已经有众多文献提出了一些不同解决方案。例如Xu提出的通过非凸低秩矩阵逼近进行动态MRI重建，(参考文献：Xu F，Han J，Wang Y，et al.DynamicMagnetic Resonance Imaging via Nonconvex Low-Rank Matrix Approximation[J].IEEE Access，2017，5:1958-1966.)。该方法利用低秩矩阵的伽马范数和拉普拉斯范数分别逼近低秩矩阵的秩，能够较好的逼近矩阵的秩，但是没有考虑稀疏矩阵的l₁范数和l₀范数之间的误差，使得重建的图像仍有少量伪影。然而，通过对上述医学图像重建方法的分析发现，现有方法在不同程度上存在伪影、细节丢失、轮廓不明显等视觉效果下降，从而影响对病情的诊断问题。

发明内容

本发明所要解决的技术问题是现有动态MRI图像重建方法在不同程度上存在伪影、细节丢失、轮廓不明显等视觉效果下降，从而影响对病情的诊断的问题，提供一种基于非凸低秩的动态MRI重建方法。

为解决上述问题，本发明是通过以下技术方案实现的：

一种基于非凸低秩的动态MRI重建方法，包括如下步骤：

步骤1，输入观测的动态MRI数据d，编码算子E；

步骤2，初始化低秩矩阵L⁽⁰⁾，稀疏矩阵S⁽⁰⁾，惩罚参数μ⁽⁰⁾，权值v⁽⁰⁾，拉格朗日乘子Y⁽⁰⁾，以及迭代次数k＝1；

步骤3，采用奇异值阈值方法求解最小化问题，得到第k次迭代低秩矩阵L^(k)：

步骤4，采用软阈值方法求解最小化问题，得到第k次迭代稀疏矩阵S^(k)：

步骤5，判断第k次迭代低秩矩阵L^(k)与稀疏矩阵S^(k)之和是否均满足收敛条件；如果满足收敛条件，转至步骤8；否则，转至步骤6；

步骤6，更新惩罚参数μ^(k)和拉格朗日乘子Y^(k)，其中

惩罚参数μ^(k)更新为：

μ^(k)＝ρμ^(k-1)，

拉格朗日乘子Y^(k)更新为：

Y^(k)＝Y^(k-1)+μ^(k-1)E^H(E(L^(k)+S^(k))-d)；

步骤7，令迭代次数k加1，并转至步骤3；

步骤8，将最终迭代得到的低秩矩阵

和稀疏矩阵

相加得到重建后的图像

上述各式中，L^(k)表示第k次迭代的低秩矩阵，S^(k)表示第k次迭代的稀疏矩阵，SVT_τ表示奇异值阈值算子，τ表示大于0的常量，

表示软阈值算子，

表示阈值，λ_S表示设定的正则化参数，μ^(k-1)表示第k-1次迭代的惩罚参数，v^(k-1)表示第k-1次迭代的权值，E表示编码算子，H表示共轭转置，d表示动态MRI数据，S^(k-1)表示第k-1次迭代的稀疏矩阵，Y^(k-1)表示第k-1次迭代的拉格朗日乘子，μ^(k)表示第k次迭代的惩罚参数，ρ表示大于1的常数，Y^(k)表示第k次迭代的拉格朗日乘子，k表示迭代次数，k＝1,2,……。

上述步骤1中，编码算子E为傅里叶变换。

上述步骤2中，初始低秩矩阵L⁽⁰⁾＝E^Hd，初始稀疏矩阵S⁽⁰⁾＝0，初始惩罚参数μ⁽⁰⁾>0，初始权值v⁽⁰⁾＝1，初始拉格朗日乘子

其中E表示编码算子，H表示共轭转置，d表示动态MRI数据，σ(·)表示取矩阵的奇异值。

上述步骤5中，收敛条件为：当第k次迭代低秩矩阵L^(k)和稀疏矩阵S^(k)之和L^(k)+S^(k)与第k-1次迭代低秩矩阵L^(k-1)和稀疏矩阵S^(k-1)之和L^(k-1)+S^(k-1)的相对误差小于预设误差值时，则说明解已满足收敛条件。

与现有技术相比，本发明将非凸正则项引入动态MRI的重建模型中，减少与真实值之间的误差，即在分离MR图像时，可以得到更精确的前景和背景。其次在视觉效果上，能够保留较多的纹理信息和轮廓信息，更便于病情的诊断。

附图说明

图1为一种基于非凸低秩模型的动态MRI重建方法的流程图。

图2是不同方法所重建出的动态心脏灌注的重建效果图；

(a-1)表示k-t SENSE方法重建出的第2帧，(b-1)表示LplusS方法重建出的第2帧，(c-1)表示本发明方法重建出的第2帧；

(a-2)表示k-t SENSE方法重建出的第8帧，(b-2)表示LplusS方法重建出的第8帧，(c-2)表示本发明方法重建出的第8帧；

(a-3)表示k-t SENSE方法重建出的第14帧，(b-3)表示LplusS方法重建出的第14帧，(c-3)表示本发明方法重建出的第14帧；

(a-4)表示k-t SENSE方法重建出的第24帧，(b-4)表示LplusS方法重建出的第24帧，(c-4)表示本发明方法重建出的第24帧。

具体实施方式

下面结合本发明具体实施例中的附图，对本发明的技术方案进行详细地描述。

本实施例以动态心脏灌注图像的重建为例所设计的一种基于非凸低秩的动态MRI重建方法，图1所示，其具体包括步骤如下：

步骤1，输入观测的动态MRI数据d，编码算子E。

输入的观测数据d属于k-t空间，编码算子E是指傅里叶变换，在Matlab中直接调用小波工具箱中的函数，其作用是将图像域数据与k-t空间的观测数据对应，即E(L+S)＝d。

步骤2，初始化低秩矩阵L⁽⁰⁾，稀疏矩阵S⁽⁰⁾，惩罚参数μ⁽⁰⁾，权值

拉格朗日乘子Y⁽⁰⁾。

迭代开始时，初始化低秩矩阵L⁽⁰⁾＝E^Hd，L⁽⁰⁾即为图像域的数据与k-t空间的初始观测数据d对应，E^H为E的共轭转置。稀疏矩阵S⁽⁰⁾＝0，惩罚参数μ⁽⁰⁾>0，权重

拉格朗日乘子

σ(·)表示取矩阵的奇异值。

步骤3，采用奇异值阈值(SVT)方法求解最小化问题

其中，||·||_*为核范数，即为矩阵L奇异值的和，E^H为E的共轭转置，λ_L为正则化参数，本实施例中的λ_L＝0.01，右上标(k)表示第k次迭代。

得到低秩矩阵L^(k)：

其中，SVT(X)可以被表示成SVT(X)＝U[diag{∑-τI}]V^T，其中(x)₊＝max(x,0)，τ是一个大于0的常量，本实施例中τ＝0.0025。

步骤4，采用软阈值方法求解最小化问题

其中，n为图像域上图像的宽和高，s_i,j为S中的第i行第j列元素，f(·)＝|·|^p，并且假设非凸函数f:R→R⁺是连续的且在[0,∞)单调递增，

是f(x)在x₀处的导数，

E^H为E的共轭转置，λ_S为正则化参数，右上标(k)表示第k次迭代。

得到稀疏矩阵S^(k)：

其中，软阈值算子被定义为

阈值为

λ_S为正则化参数，本发明取λ_S＝0.03，同时更新权值

i＝1,2,...,n，j＝1,2,...,n，

ε是一个很小的正实数，本实施例中ε＝0.0001。p的范围是(0,1)，通过实验可知，p值的选取严重影响着收敛速度。当p值越大时，函数的收敛速度越快。但当p值越小时，解的精确度越高。因此p的取值要适当，本实施例中的p均取0.2。

步骤5，判断步骤3和步骤4所得到的第k次迭代低秩矩阵L^(k)与稀疏矩阵S^(k)之和是否均满足收敛条件；如果满足收敛条件，转至步骤8；否则，转至步骤6；

步骤6，分别更新惩罚参数μ^(k)和拉格朗日乘子Y^(k)，其中

惩罚参数μ^(k)更新为：

μ^(k)＝ρμ^(k-1)，

拉格朗日乘子Y^(k)更新为：

Y^(k)＝Y^(k-1)+μ^(k-1)E^H(E(L^(k)+S^(k))-d)；

其中，ρ为常数，一般情况下ρ>1，本实施例中ρ＝1.2。

步骤7，令迭代次数k加1，并转至步骤3；

步骤8，将最终迭代得到的低秩矩阵

和稀疏矩阵

相加得到重建后的图像

在本实施例中，收敛条件是L^(k)+S^(k)与前一次的L^(k-1)+S^(k-1)的相对误差小于2.5×10^-3时，说明解已达到收敛。

本发明的效果由以下仿真进一步说明。

仿真条件：仿真的运行环境为Matlab 2014a，Inter Core i5-4690处理器，8GB的内存，win7 64位操作系统。

仿真内容：仿真数据是40帧大小为128×128的心脏灌注数据。另外，客观评价指标分别为相对误差和均方根误差，相对误差作为仿真迭代停止的指标，形式如下：

均方根误差为重建后图像质量的评价指标，形式如下：

其中，d为观测数据，L,S分别为重建图像相应的低秩矩阵和稀疏矩阵。Err越小说明收敛速度越快，当Err小于2.5×10^-3时，迭代停止，RMSE越小说明重建误差越小。用本发明与k-t SENSE，LplusS分别对以上数据进行重建，从以上两方面进行比较。

仿真结果：

本发明比k-t SENSE方法和LplusS方法在视觉上取得较好的重建效果。仿真结果如图2所示。从图2中可以看出，对于这四帧图像，三种方法都呈现出良好的重建效果，但是k-t SENSE和LplusS重建后的较多区域出现模糊，背景不清晰，部分区域被噪点污染的情况。正如图中红框中的部分所示，所提方法重建的主体部分与边缘组织的轮廓相比于这两种方法都较清晰，亮度也较亮。

另外，从客观上比较本发明和k-t SENSE和LplusS的重建效果，如表1所示。

表1

从表1中可以看出，在对动态心脏灌注图像的重建中，本发明方法和LplusS的总耗时都比k-t SENSE的少，LplusS的总耗时最少。而本发明方法的Err和RMSE比k-t SENSE和LplusS的低很多，说明本发明方法的重建效果较好。

本发明在非凸优化问题研究的基础上，引入非凸正则项，增大图像的稀疏性，然后用交替方向乘子法求解模型，得到迭代更新后的低秩矩阵和稀疏矩阵，将两者相加得到重建的MRI图像。本发明的优点是将非凸正则项引入动态MRI的重建模型中，减少与真实值之间的误差，使求得的解更精确。并且利用交替方向乘子法求解模型，大大减少了计算量，使计算更简单快捷，且用该方法求得的解具有良好的收敛性。实验结果在主观视觉和客观评价指标两方面均取得较好的效果，能够保留较多的纹理信息和轮廓信息，更便于病情的诊断。

需要说明的是，尽管以上本发明所述的实施例是说明性的，但这并非是对本发明的限制，因此本发明并不局限于上述具体实施方式中。在不脱离本发明原理的情况下，凡是本领域技术人员在本发明的启示下获得的其它实施方式，均视为在本发明的保护之内。

Claims (3)

1.一种基于非凸低秩的动态MRI重建方法，其特征是，包括如下步骤：

步骤1，输入观测的动态MRI数据d，编码算子E；

步骤2，初始化低秩矩阵L⁽⁰⁾，稀疏矩阵S⁽⁰⁾，惩罚参数μ⁽⁰⁾，权值v⁽⁰⁾，拉格朗日乘子Y⁽⁰⁾，以及迭代次数k＝1；

步骤3，采用奇异值阈值方法求解最小化问题，得到第k次迭代低秩矩阵L^(k)，

步骤4，采用软阈值方法求解最小化问题，得到第k次迭代稀疏矩阵S^(k)，

步骤5，判断第k次迭代低秩矩阵L^(k)与稀疏矩阵S^(k)之和是否均满足收敛条件；即当第k次迭代低秩矩阵L^(k)和第k次迭代稀疏矩阵S^(k)之和L^(k)+S^(k)与第k-1次迭代低秩矩阵L^(k-1)和第k-1次迭代稀疏矩阵S^(k-1)之和L^(k-1)+S^(k-1)的相对误差小于预设误差值时，则说明解已满足收敛条件：如果满足收敛条件，转至步骤8；否则，转至步骤6；

步骤6，更新惩罚参数μ^(k)和拉格朗日乘子Y^(k)，其中

惩罚参数μ^(k)更新为，

μ^(k)＝ρμ^(k-1)，

拉格朗日乘子Y^(k)更新为，

Y^(k)＝Y^(k-1)+μ^(k-1)E^H(E(L^(k)+S^(k))-d)；

步骤7，令迭代次数k加1，并转至步骤3；

步骤8，将最终迭代得到的低秩矩阵

和稀疏矩阵

相加得到重建后的图像

上述各式中，L^(k)表示第k次迭代的低秩矩阵，S^(k)表示第k次迭代的稀疏矩阵，SVT_τ表示奇异值阈值算子，τ表示大于0的常量，

表示软阈值算子，

表示阈值，λ_S表示设定的正则化参数，μ^(k-1)表示第k-1次迭代的惩罚参数，v^(k-1)表示第k-1次迭代的权值，E表示编码算子，H表示共轭转置，d表示动态MRI数据，S^(k-1)表示第k-1次迭代的稀疏矩阵，Y^{(k -1)}表示第k-1次迭代的拉格朗日乘子，μ^(k)表示第k次迭代的惩罚参数，ρ表示大于1的常数，Y^(k)表示第k次迭代的拉格朗日乘子，k表示迭代次数，k＝1,2,……。

2.根据权利要求1所述的一种基于非凸低秩的动态MRI重建方法，其特征是，步骤1中，编码算子E为傅里叶变换。

3.根据权利要求1所述的一种基于非凸低秩的动态MRI重建方法，其特征是，步骤2中，初始低秩矩阵L⁽⁰⁾＝E^Hd，初始稀疏矩阵S⁽⁰⁾＝0，初始惩罚参数μ⁽⁰⁾＞⁰，初始权值v⁽⁰⁾＝1，初始拉格朗日乘子

其中E表示编码算子，H表示共轭转置，d表示动态MRI数据，σ(·)表示取矩阵的奇异值。

Images