CN104156994B - 一种压缩感知磁共振成像的重建方法 - Google Patents

Abstract

一种压缩感知磁共振成像的重建方法，涉及图像处理。提供可提升重建图像的主观视觉效果的一种压缩感知磁共振成像的重建方法。在基于非局部相似块构成的低秩矩阵下构造MRI图像重建的目标函数；根据构造出的目标函数，利用变量替换方法将目标函数的优化求解问题转化为低秩矩阵的去噪与目标图像重建的求解问题；对于低秩矩阵的去噪问题，采用奇异值分解方法将聚类的低秩矩阵进行分解，对分解得到的特征值进行软阈值处理，获得去噪后的低秩矩阵；将得到的低秩矩阵带入到目标函数中，经过近似优化，利用最小二乘算法得到最终重建后的MR图像。

Description

一种压缩感知磁共振成像的重建方法

技术领域

本发明涉及图像处理，具体是涉及一种压缩感知磁共振成像的重建方法。

背景技术

磁共振成像技术(MRI)是利用磁场共振原理成像，能够使人们无损伤地获取活体器官和组织的详细诊断图像，避免了不必要的手术痛苦以及探查性手术所带来的副损失及并发症。由于其能够给医生提供清晰、精细、分辨率高、对比度好、信息量大的人体结构医学图像，目前已普遍应用于临床，已成为一些疾病诊断必不可少的检查手段。

然而MRI的不足之处就在于成像速度慢，人在成像过程中必须保持静止状态，即使是轻微的运动也会使成像产生伪影，影响临床的诊断。压缩感知(CS)理论的出现使得人们在获取低数据量的情况下也能够获得高分辨率的图像。CS理论认为在K空间对原始数据进行下采样(采集的数据可以远小于全采样的数据)，只要满足获得的MR图像在某一变换域上是稀疏的，则可以从这一小部分K空间数据中精确恢复出原图像，从而解决了在成像速度快的条件下也能够获得分辨率高的MR图像。

传统的压缩感知MRI图像的重建方法，是利用MR图像的稀疏性来重建的。随着稀疏表示理论在MRI领域的应用，越来越多的方法开始被提出来，大部分方法可以归结为两类：一类是利用图像在某一变换域下的稀疏表示重建图像，例如Lusting(Lusting et al.,Sparse mri:The application of compressed sensing for rapid mr imaging,Magnetic Resonance in Medicine,vol.58,no.6,pp.1182–1195,2007)以及屈小波(X.Quet al.,Undersampled mri reconstruction with patchbased directional wavelets,Magnetic Resonance Imaging,vol.30,no.7,pp.964–977,2012)提出的小波域稀疏方法；另一类是根据样本学习字典来重建数据，例如S.Ravishankar(S.Ravishankar et al.,Mrimage reconstruction from highly undersampled k-space data bydictionarylearning,IEEE Trans.on Medical Imaging,vol.30,no.5,pp.1028–1041,2011)提出的KSVD方法来训练字典，然后重构MR图像。最近，已有一些学者对图像的非局部先验(特殊的稀疏性)进行了研究，Kostadin(Kostadin et al.,Image denoising bysparse 3d transform-domain collaborative filtering,IEEE Trans.on ImageProcessing,vol.16,no.8,pp.2080–2095,2007)提出基于协同滤波的BM3D方法，董伟生(W.Dong,Nonlocal image restoration with bilateral variance estimation:a low-rank approach,IEEE Trans.Image Processing,vol.22,no.2,pp.700–711,2013)提出的一种基于非局部先验的低秩方法。但是这些方法没有充分利用非局部以及低秩先验知识来提升MR图像重建算法性能。

发明内容

本发明的目的是针对现有的MRI图像重建方法的不足，提供可提升重建图像的主观视觉效果的一种压缩感知磁共振成像的重建方法。

本发明包括如下步骤：

A、在基于非局部相似块构成的低秩矩阵下构造MRI图像重建的目标函数；

B、根据构造出的目标函数，利用变量替换方法将目标函数的优化求解问题转化为低秩矩阵的去噪与目标图像重建的求解问题；

C、对于低秩矩阵的去噪问题，采用奇异值分解方法将聚类的低秩矩阵进行分解，对分解得到的特征值进行软阈值处理，获得去噪后的低秩矩阵；

D、将得到的低秩矩阵带入到目标函数中，经过近似优化，利用最小二乘算法得到最终重建后的MR图像。

在步骤A中，所述在基于非局部相似块构成的低秩矩阵下构造MRI图像重建的目标函数的具体步骤如下：

利用低秩矩阵的特性，构造基本MRI重构模型：

其中，y表示经过磁共振扫描仪获得的观测数据，F_u为部分傅立叶变换操作算子，x为重建后的图像；已知，a_i表示在图像x中以像素点i为中心的图像块，图像块集合为N表示图像总的像素点个数，A_i表示以图像块a_i为参考块，在其一定领域内寻找与a_i最相似的M块图像块将这些相似块拉成列向量，构成的相似块矩阵；rank(A_i)表示对矩阵A_i求秩操作；

式(1)是一个非凸优化问题，目前还没有一个很好的解决手段，可以放松约束条件，用核范数最小替代秩最小，所以式(1)可以改写成:

其中，矩阵A_i的核范数||A_i||_*定义为λ_i,j为矩阵A_i的特征值，r表示矩阵A_i的秩大小，λ＞0是权重参数。

在步骤B中，所述根据构造出的目标函数，利用变量替换方法将目标函数的优化求解问题转化为低秩矩阵的去噪与目标图像重建的求解问题的具体方法如下：

利用变量替换的方法将目标模型(2)变换为如下形式：

其中B_i是在求解MR图像过程中产生的中间变量，称之为辅助变量，然后再将这个目标函数写成非约束形式：

其中，定义为：x_i,j表示矩阵X中坐标(i,j)的像素值，ω为常数；对于式(4)，采用交替最小化过程将目标函数转化为低秩矩阵的复原以及目标图像重建问题。

在步骤C中，所述对于低秩矩阵的去噪问题，采用奇异值分解方法将聚类的低秩矩阵进行分解，对分解得到的特征值进行软阈值处理，获得去噪后的低秩矩阵的具体步骤如下：

对于目标函数(4)，当x已知的情况下，求解未知变量B_i，得到下式：

对低秩矩阵A_i进行一个阈值去噪，得到新的低秩矩阵B_i，具体可分为两步：

a、对低秩矩阵A_i进行奇异值分解(SVD)：

(U_i,Σ_i,V_i)＝svd(A_i)

其中，Σ_i表示对矩阵A_i分解获得的特征值矩阵，而U_i与V_i为分解得到的酉矩阵；

b、对特征值矩阵进行软阈值处理阈值其中，r表示矩阵A_i秩大小，估计是根据矩阵A_i所包含的噪声大小来确定的，噪声方差估计如下：

其中，S表示矩阵A_i的特征值由第1个至第k个特征值之和小于A_i总能量的85％所对应的最大k值，1＜k＜r，n、p分别表示A_i的行数与列数，min(n-1,p)表示取n-1与p之间较小的一方的值，特征值阈值为

通过阈值去噪后得到

在步骤D中，所述将得到的低秩矩阵带入到目标函数中，经过近似优化，利用最小二乘算法得到最终重建后的MR图像的具体步骤如下：

对于目标函数(4)，当B_i已知的情况下，求解未知变量x，得到下式：

其中，λ₁＝λ·ω，为了能够把重构MR图像x写成闭型解的形式，需要对式(6)进行改写以便能够得到一个最小二乘的解，由于可以写成向量的和的形式，所以可以得到：

其中，N是图像x像素个数，a_i表示图像x中第i个像素点对应的图像块向量，b_i表示与a_i对应的中间变量，O_i表示中出现的次数，因此，必然满足：

其中，O_min和O_max分别表示的最小值和最大值，因此可以对取一个近似的估计(ρ>0)；因此式(6)可以写为

λ₂是一个权重常数，由于是复原的所有图像块的集合，通过图像块加权平均，得同样，通过块加权平均得到图像x，因此，可以将式(8)转变成式(9)：

λ₃是大于0的常数；

在式(9)的基础上，利用ADMM算法在图像域以及空间域的约束项中同时引入噪声回加过程，可得到式(10)：

式(10)是一个典型的最小二乘问题，从而能够得到重构后的x'，其中，E_x与E_k是噪声回加变量，更新的E_x'与E_k'分别为：E_k′＝E_k+F_ux′-y。

本发明利用非局部块的相似性以及相似块构成的低秩矩阵这两种先验知识，提出了一种新的算法，利用自适应的参数估计以及最小二乘算法来很好的重建MR图像。通过实验证明了该算法模型相比于传统的方法在对MR图像的细节边缘处理上具有很好的重建效果。

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

第一，本发明利用非局部块相似、低秩两种先验，提出了一种新的算法模型；

第二，本发明在对特征值阈值处理过程中采取不同特征值对应不同的阈值，并且其参数是自适应的；

第三，本发明在解法上采用一种近似最小二乘的思路来重建MR图像，相比于现有技术在边缘及边界区域复原更加清晰。

附图说明

图1为本发明仿真实验所用的磁共振MRI血管(circle)图像及对比实验结果。在图1中，(a)为原始图像，(b)为PBDW，(c)为BPFA+TV，(d)为本发明所提方法。

图2本发明针对血管图像在笛卡尔采样条件下重建结果的PSNR。

图3本发明针对血管图像在径向采样条件下重建结果的PSNR。

具体实施方式

本发明是基于非局部、低秩先验模型的压缩感知MRI图像重建方法，其实现步骤如下：

步骤1，对输入磁共振MRI原始K空间观测数据y进行初始化处理，将其进行反傅立叶变换，获得填零重建后的初始图像x0。

步骤2，创建初始化图像x0的相似块矩阵，生成相似块矩阵索引集合G_index

(2a)将初始化的图像x0按步长为1进行取块操作，得到全部图像块的集合标注每个图像块在初始化图像中对应的位置

(2b)以像素点i为中心，取图像块a_i作为参考块，在i点为中心的邻近区域内寻找与其相似的M个相似块创建相似块矩阵A_i，A_i＝[a_i,1,a_i,2,…,a_i,M]，其中,表示M个相似块的位置索引(g_i,j表示相似块a_i,j的索引位置)。

(2c)根据相似块矩阵索引生成所有块相似矩阵索引集合G_index：

G_index＝{g_i,j},i＝1,2,...,N；j＝1,2,...,M。

步骤3，使用基于低秩矩阵优化算法重构磁共振MRI图像x

(3a)设置循环迭代次数iternum＝1,2,...,T,T表示最大迭代次数，在实验中T＝300

(3b)对相似块矩阵A_i进行去噪处理，即

首先，将相似块矩阵A_i进行奇异值分解(SVD)：(U_i,Σ_i,V_i)＝svd(A_i)，其中，svd表示奇异值分解符号，酉矩阵U_i、V_i中的每一列均是与特征值矩阵Σ_i中对应的特征值有关的特征向量。Σ_i为对角线矩阵，λ_i,j为对角线上得特征值，r表示矩阵A_i秩大小。

其次，对特征值矩阵进行软阈值处理阈值估计是根据矩阵A_i所包含的噪声大小来确定的，噪声估计如下：

其中，S表示矩阵A_i的特征值由第1个至第k个(1＜k＜r)特征值之和小于A_i总能量的85％所对应的最大k值。n、p分别表示A_i的行数与列数，min(n-1,p)表示取n-1与p之间较小的一方的值，特征值阈值为

通过阈值去噪后得到

重复步骤(3b)，待将块相似索引集G_index中所指引的相似块矩阵全部处理完为止。

(3c)将图像的全局保真项与低秩矩阵恢复的非局部保真项结合，得到目标函数：

其中，λ₁＞0是一个正则化参数；为了能够把重构MR图像x写成闭型解的形式，需要对式(12)进行改写以便能够得到一个最小二乘的解。由于可以写成向量的和的形式，所以可以得到其中N是图像像素个数，O_i表示在中出现的次数。因此，必然满足：

O_min和O_max分别表示的最小值和最大值。因此可以对取一个近似的估计(ρ>0)。因此式(12)可以写为

λ₂是一个常数。由于是复原的所有图像块的集合，通过图像块加权平均，同样，通过块加权平均得到图像x。因此，可以将式(13)转变成式(14)：

(3d)在式(14)的基础上，引入交替方向乘子算法(ADMM)，利用ADMM算法在图像域以及空间域的约束项中同时引入噪声回加过程，可得到式(15)：

其中E_x与E_k是噪声回加变量，更新的E_x'与E_k'分别为：

E_k′＝E_k+F_ux′-y

x'表示重建后的MR图像。ADMM算法放松了对参数λ₃'的约束，对于给定一个参数，在交替迭代的过程中必然收敛，在实验中，设λ₃'＝0.2。

(3e)求解步骤(3d)目标函数，得到优化的图像x：

其中，等式中H表示共轭转置操作，I表示单位阵。通过直接求解式(16)得到优化的图像x将涉及到非常大的计算量(因为在求解x过程中需要求解大小为N×N这样一个非常庞大的逆矩阵)，因此需要进一步的优化。将等式两边同时取傅立叶变换：

其中，F表示傅立叶变换操作算子，矩阵FF_u ^HF_uF^H是一个对角阵，对角阵上的数值是由0和1组成，1表示在K空间上数据被采到的点，0表示未采到的点的位置(实际上就是K空间的下采样模板)。

令可得：

Ω表示采到的K空间子集，(k_x,k_y)表示K空间所对应的坐标，当(k_x,k_y)∈Ω，则对x进行更新，否则保持不变。最后去反傅立叶变换得到重建后的x'。

(3f)重复步骤2与步骤3共T次，得到最终的优化图像x',x'即为低秩矩阵的压缩感知MRI图像重建方法的最终结果。

本发明同目前国际上几种主流的MRI图像重建算法进行比较，它们分别是Lustig(Lusting et al.,Sparse mri:The application of compressed sensing for rapid mrimaging,Magnetic Resonance in Medicine,vol.58,no.6,pp.1182–1195,2007)提出的SparseMRI、RaVishankar(S.Ravishankar et al.,Mr image reconstruction fromhighly undersampled k-space data by dictionary learning,IEEE Trans.on MedicalImaging,vol.30,no.5,pp.1028–1041,2011)提出的DLMRI、屈小波(X.Qu et al.,Undersampled mri reconstruction with patchbased directional wavelets,MagneticResonance Imaging,vol.30,no.7,pp.964–977,2012)提出的PBDW及丁兴号(X.Ding etal.,Compressed sensing mri with Bayesian dictionary learning,IEEEInternational Conference on Image Processing(ICIP),2013)提出的BPFA+TV，选取一幅大小为512×512的血管图像，在配置为Pentium(R)Dual-Core E53002.60GHz内存4GB，Matlab7.5.0的PC机上进行实验对比。从图1的重建结果可以看出BPFA+TV的重建结果远好于PBDW的重建结果，但是同本发明的重建结果进行对比可以发现，本发明提出的方法在边缘及边界区域重建结果更加清晰。从图2和3的结果可以看出无论是在笛卡尔采样还是在径向采样下，本发明所提出的方法的峰值信噪比整体都要好于几种经典的方法。

Claims (4)

1.一种压缩感知磁共振成像的重建方法，其特征在于包括如下步骤：

A、在基于非局部相似块构成的低秩矩阵下构造MRI图像重建的目标函数，具体步骤如下：

利用低秩矩阵的特性，构造基本MRI重构模型：

$\underset{x}{\mathrm{argmin}}\left|\right|F_{u}x-y|{|}_2^2+rank\left(A_i\right);i=1,...,N---\left(1\right)$

其中，y表示经过磁共振扫描仪获得的观测数据，F_u为部分傅立叶变换操作算子，x为重建后的图像；已知，a_i表示在图像x中以像素点i为中心的图像块，图像块集合为N表示图像总的像素点个数，A_i表示以图像块a_i为参考块，在其邻域内寻找与a_i最相似的M块图像块将这些相似块拉成列向量，构成的相似块矩阵；rank(A_i)表示对矩阵A_i求秩操作；

式(1)是一个非凸优化问题，用核范数最小替代秩最小，式(1)改写成:

$\underset{x}{\mathrm{argmin}}\left|\right|F_{u}x-y|{|}_2^2+\mathrm{\λ}{\Σ}_i|\left|A_i\right|{|}_{*}---\left(2\right)$

其中，矩阵A_i的核范数||A_i||_*定义为λ_i,j为矩阵A_i的特征值，r表示矩阵A_i的秩大小，λ＞0是权重参数；

B、根据构造出的目标函数，利用变量替换方法将目标函数的优化求解问题转化为低秩矩阵的去噪与目标图像重建的求解问题；

C、对于低秩矩阵的去噪问题，采用奇异值分解方法将聚类的低秩矩阵进行分解，对分解得到的特征值进行软阈值处理，获得去噪后的低秩矩阵；

D、将得到的低秩矩阵带入到目标函数中，经过近似优化，利用最小二乘算法得到最终重建后的MRI图像。

2.如权利要求1所述一种压缩感知磁共振成像的重建方法，其特征在于在步骤B中，所述根据构造出的目标函数，利用变量替换方法将目标函数的优化求解问题转化为低秩矩阵的去噪与目标图像重建的求解问题的具体方法如下：

利用变量替换的方法将目标模型(2)变换为如下形式：

$\begin{array}{cc}\underset{x}{\mathrm{argmin}}\left|\right|F_{u}x-y|{|}_2^2+\mathrm{\λ}{\Σ}_i|\left|B_i\right|{|}_{*}& s.t.A_i=B_i\end{array}---\left(3\right)$

其中B_i是在求解MRI图像过程中产生的中间变量，称之为辅助变量，然后再将这个目标函数写成非约束形式：

$\underset{x}{\mathrm{argmin}}\left|\right|F_{u}x-y|{|}_2^2+{\mathrm{\λ\Σ}}_i\{\mathrm{\ω}\left|\right|B_i-A_i|{|}_F^2+|\left|B_i\right|{|}_{*}\}---\left(4\right)$

其中，定义为：x_i,j表示矩阵X中坐标(i,j)的像素值，ω为常数；对于式(4)，采用交替最小化过程将目标函数转化为低秩矩阵的复原以及目标图像重建问题。

3.如权利要求2所述一种压缩感知磁共振成像的重建方法，其特征在于所述对于低秩矩阵的去噪问题，采用奇异值分解方法将聚类的低秩矩阵进行分解，对分解得到的特征值进行软阈值处理，获得去噪后的低秩矩阵的具体步骤如下：

对于目标函数(4)，当x已知的情况下，求解未知变量B_i，得到下式：

$\underset{B_i}{\mathrm{argmin}}{\mathrm{\Σ}}_i\left\{\mathrm{\ω}\right||B_i-A_i|{|}_F^2+\left|\right|B_i\left|{|}_{*}\right\}---\left(5\right)$

对低秩矩阵A_i进行一个阈值去噪，得到新的低秩矩阵B_i，具体分为两步：

a、对低秩矩阵A_i进行奇异值分解(SVD)：

(U_i,Σ_i,V_i)＝svd(A_i)

其中，Σ_i表示对矩阵A_i分解获得的特征值矩阵，而U_i与V_i为分解得到的酉矩阵；

b、对特征值矩阵进行软阈值处理阈值其中，r表示矩阵A_i秩大小，估计是根据矩阵A_i所包含的噪声大小来确定的，噪声方差估计如下：

${{\widehat{\mathrm{\σ}}}_i}^2=\frac{{\Σ}_{j=S+1}^{\min (n-1,p)}{\mathrm{\λ}}_{i,j}}{np-p-nS-pS+S+S^2}$

其中，S表示矩阵A_i的特征值由第1个至第k个特征值之和小于A_i总能量的85％所对应的最大k值，1＜k＜r，n、p分别表示A_i的行数与列数，min(n-1,p)表示取n-1与p之间较小的一方的值，特征值阈值为

${\mathrm{\τ}}_{i,j}=\frac{np{{\widehat{\mathrm{\σ}}}_i}^2}{min(n-1,p){\mathrm{\λ}}_{i,j}}$

通过阈值去噪后得到

4.如权利要求3所述一种压缩感知磁共振成像的重建方法，其特征在于所述将得到的低秩矩阵带入到目标函数中，经过近似优化，利用最小二乘算法得到最终重建后的MRI图像的具体步骤如下：

对于目标函数(4)，当B_i已知的情况下，求解未知变量x，得到下式：

$\underset{x}{\mathrm{argmin}}\left|\right|F_{u}x-y|{|}_2^2+{\mathrm{\λ}}_1{\mathrm{\Σ}}_i||A_i-B_i|{|}_F^2---\left(6\right)$

其中，λ₁＝λ·ω，为了能够把重构MRI图像x写成闭型解的形式，需要对式(6)进行改写以便能够得到一个最小二乘的解，由于可以写成向量的和的形式，所以可得到：

${\mathrm{\Σ}}_i\left|\right|A_i-B_i|{|}_F^2={\mathrm{\Σ}}_i{O}_i||a_i-b_i|{|}_2^2,i=1...N$

其中，N是图像x像素个数，a_i表示图像x中第i个像素点对应的图像块向量，b_i表示与a_i对应的中间变量，O_i表示在中出现的次数，因此，必然满足：

$O_{\min }{\mathrm{\Σ}}_i\left|\right|a_i-b_i|{|}_2^2\≤{\mathrm{\Σ}}_i||A_i-B_i|{|}_F^2\≤O_{max}{\mathrm{\Σ}}_i\left|\right|a_i-b_i|{|}_2^2---\left(7\right)$

其中，O_min和O_max分别表示的最小值和最大值，因此可以对取一个近似的估计(ρ>0)；因此式(6)可以写为

$\underset{x}{\mathrm{argmin}}\left|\right|F_{u}x-y|{|}_2^2+{\mathrm{\λ}}_2{\mathrm{\Σ}}_i||a_i-b_i|{|}_2^2---\left(8\right)$

λ₂是一个权重常数，由于是复原的所有图像块的集合，通过图像块加权平均，得同样，通过块加权平均得到图像x，因此，可以将式(8)转变成式(9)：

$\underset{x}{\arg \min }\left|\right|F_{u}x-y|{|}_2^2+{\mathrm{\λ}}_3||x-{\hat{x}}_{PCA}|{|}_2^2---\left(9\right)$

λ₃是大于0的常数；

在式(9)的基础上，利用ADMM算法在图像域以及空间域的约束项中同时引入噪声回加过程，得到式(10)：

$\underset{x}{\mathrm{argmin}}\left|\right|F_{u}x-y+E_k|{|}_2^2+{\mathrm{\λ}}_3||x-{\hat{x}}_{PCA}+E_x|{|}_2^2---\left(10\right)$

式(10)是一个典型的最小二乘问题，从而能够得到重构后的x'，其中，E_x与E_k是噪声回加变量，更新的E_x'与E_k'分别为：E_k′＝E_k+F_ux′-y。