CN104739410B - 一种磁共振图像的迭代重建方法 - Google Patents

Abstract

一种磁共振图像的迭代重建方法，涉及磁共振成像。提供计算简单、运算速度快、占用计算机内存少的一种磁共振图像的迭代重建方法。首先提出一个解空间受限的稀疏重建模型，再通过迭代运算求解重建模型并获得重建磁共振图像，其中在每次迭代中利用相邻两次迭代的重建图像对图像进行修正，用来加快图像重建的速度。这种重建方法的参数少，重建的迭代运算复杂度低，运算速度快。同时，图像重建过程不需要额外计算变换域系数，占用少量计算机内存。

Description

一种磁共振图像的迭代重建方法

技术领域

本发明涉及磁共振成像，尤其是涉及一种磁共振图像的迭代重建方法。

背景技术

磁共振成像在医学临床诊断中有着广泛的应用，但在传统的磁共振成像中，获取完整的k空间数据这个过程耗时较长。成像速度慢使得磁共振成像在功能磁共振成像、心脏和腹部磁共振成像等应用中的图像质量不理想。一种加快磁共振成像的方法就是通过对k空间数据的欠采样(Sub-Nyquist Sampling)来减少采集的数据量，达到缩短成像时间的目的。欠采样的k空间需要进行图像的稀疏重建得到完整的磁共振图像(M.Lustig,D.Donoho,and J.M.Pauly,"Sparse MRI:The application of compressed sensing for rapid MRimaging,"Magnetic Resonance in Medicine,vol.58,pp.1182-1195,2007)。

不同的图像重建方法在重建的速度和重建图像的质量上存在明显差异。这些重建方法通常基于两种不同的图像重建模型：综合模型和分解模型。对于综合模型，A.Beck等(Afast iterative shrinkage-thresholding algorithm for linear inverse problems,SIAM Journal on Imaging Sciences,vol.2,pp.183-202,2009)提出快速迭代软阈值方法来解决综合模型问题。这种方法运算复杂度低且算法速度快。但信号处理的研究表明，当使用框架来稀疏磁共振图像的时候，综合型方法的图像重建误差大于分解型的的图像重建误差。对于分解模型，Qu等(Undersampled MRI reconstruction with patch-baseddirectional wavelets,Magnetic Resonance Imaging,vol.30,pp.964-977,2012.)和Afonso等(,J.M.Bioucas-Dias,and M.A.T.Figueiredo,"Fast Image Recovery UsingVariable Splitting and Constrained Optimization,"IEEE Transactions on ImageProcessing,vol.19,pp.2345-2356,2010.)提出交替方向法来重建磁共振图像。交替方向法的优点是图像重建速度很快和重建误差小。但是这类方法需要设置的参数比较多，当使用框架变换来进行图像重建时占用大量计算机内存。

发明内容

本发明的目的在于提供计算简单、运算速度快、占用计算机内存少的一种磁共振图像的迭代重建方法。

本发明包括以下步骤：

1)提出解空间受限的稀疏重建模型：测量到的欠采样的k空间数据y表示为y＝UFx，其中x表示待重建的图像，F表示对磁共振图像x进行傅里叶变换，Fx就是完整的k空间数据，U表示对完整的k空间数据进行欠采样；用Ψ表示任意的框架变换，Ψ^*表示框架变换的复共轭，Φ表示框架Ψ的对偶框架变换，对偶框架变换Φ通过矩阵Ψ^*Ψ的逆(Ψ^*Ψ)^-1乘以Ψ^*来表示，也就是Φ＝(Ψ^*Ψ)^-1Ψ^*，其中符号“-1”表示矩阵求逆；将N个像素的图像x表示成复数域中的一个向量图像x在框架变换Ψ的变换域系数α构成的空间记为Range(Ψ)，也就是提出的解空间受限的稀疏重建模型为：

其中正则化参数λ用于权衡||Ψx||₁和两项的重要性，α∈Range(Ψ)表示变换域系数α是系数空间Range(Ψ)中的一个向量。

2)重建图像：解空间受限的稀疏重建模型(1)中的变换域系数α通过公式(2)的迭代运算来求解：

α_k+1＝T_γλ(ΨΦα_k+γΦ^*F^*U^T(y-UFΦα_k)) (2)

其中α_k和α_k+1是在第k和k+1次迭代时得到的变换域系数，U^T表示U的转置，F^*表示反傅里叶变换，γ是步长，其取值范围为γ＞0；步长γ与正则化参数λ的乘积为γλ，T_γλ表示对给定向量z＝(ΨΦα_k+γΦ^*F^*U^T(y-UFΦα_k))的所有元素进行软阈值操作，软阈值操作的定义为：给定向量z，若向量z的第j个元素z_j的绝对值|z_j|≤γλ，则z_j＝0；若|z_j|＞γλ，则z_j＝sgn(z_j)(|z_j|-γλ)，其中sgn(z_j)是符号函数。

由于待重建图像x用对偶框架Φ与变换域系数α的乘积来表示，即x＝Φα，因此图像x可按照公式(3)的迭代运算获得

x_k+1＝ΦT_γλ(Ψx_k+γΦ^*F^*U^T(y-UFx_k)) (3)

其中x_k和x_k+1是在第k和k+1次迭代时得到的图像。

3)加快图像重建的速度：引入一组因子{t_k}，其中k是大于等于0的整数，初始化的因子t₀＝1，已知这组因子中的第k项因子t_k，通过表达式可以计算第k+1项因子t_k+1；在第k+1次迭代重建图像后，利用第k和k+1次迭代时得到的图像x_k和x_k+1对迭代的图像进行修正，从而加快图像重建的收敛速度，修正的方法是：

公式(4)中的符号表示将公式(4)中右边的结果赋值给公式(4)左边的x_k+1中。

迭代过程中，当相邻两次的重建图像x_k和x_k+1的L2范数误差||x_k+1-x_k||₂/||x_k||₂缩小到指定的范围时，停止迭代运算。将最后一步迭代的重建结果x_k+1输出为最终的重建结果图像。

本发明的有益效果是：图像重建的参数少，仅引入步长γ和正则化参数λ；重建的迭代运算复杂度低，运算速度快。同时，图像重建过程不需要额外计算变换域系数，占用少量计算机内存。

附图说明

图1是实施例中进行k空间欠采样的采样模板。

图2是根据图1的采样模板得到欠采样的k空间数据。

图3是迭代过程中的人脑图像和全采样人脑图像的L2范数误差。

图4是使用本发明方法从图2中所示的k空间数据中重建出的人脑图。

具体实施方式

本发明实施例使用磁场强度为3特斯拉的磁共振成像仪对自愿者的大脑进行成像。本实施例使用的序列为T2加权快速自旋回波，序列的回波时间TE＝99ms，重复时间TR＝610ms，选层厚度Δd＝3mm，视野FOV＝256mm×256mm。

具体实施过程如下：

第一步：提出解空间受限的稀疏重建模型

对于实施例中重建256×256个像素的人脑磁共振图像x，根据图1中的采样模板控制频率和相位编码，使得采样模板中白色对应位置的k空间数据采集到。欠采到的k空间数据y(如图2所示)总共有26368点。欠采样的k空间数据y与重建图像x的关系式为：y＝UFx。其中F表示对磁共振图像x进行傅里叶变换，预期得到的65536个傅里叶系数就是待重建的完整的k空间数据，U表示对完整的k空间数据根据图1中的采样模板进行欠采样。采用尺度数为4的冗余Daubechies小波基构成框架变换Ψ。用Ψ^*表示框架变换Ψ的复共轭，Ψ对应的对偶框架变换Φ通过矩阵Ψ^*Ψ的逆(Ψ^*Ψ)^-1乘以Ψ^*来求解，也就是Φ＝(Ψ^*Ψ)^-1Ψ^*(符号“-1”表示矩阵求逆)。提出的解空间受限的稀疏重建模型为：

其中α是变换域系数，Range(Ψ)是所有65536个像素图像x在框架Ψ投影下的变换域系数构成的系数空间，α∈Range(Ψ)表示变换域系数α是系数空间Range(Ψ)中的一个向量，正则化参数λ用于权衡||Ψx||₁和的重要性。实施例人脑图像的重建中，设定正则化参数λ＝0.01。

第二步：重建图像

人脑图像x通过公式(6)的迭代运算来获得

x_k+1＝ΦT_γλ(Ψx_k+γΦ^*F^*U^T(y-UFx_k)) (6)

其中x_k和x_k+1是在第k和k+1次迭代时得到的图像，U^T表示U的转置，F^*表示反傅里叶变换，γ是步长，实施例中设置迭代步长γ＝1。下标步长γ与正则化参数λ的乘积为γλ＝0.01。T_γλ表示对给定向量z＝(Ψx_k+γΦ^*F^*U^T(y-UFx_k))的所有元素进行软阈值操作，软阈值操作的定义为：给定向量z，若向量z的第j个元素z_j的绝对值|z_j|≤γλ，则z_j＝0，若|z_j|＞γλ，则z_j＝sgn(z_j)(|z_j|-γλ)，其中sgn(z_j)是符号函数。

第三步：加快图像重建的速度

在第k+1次迭代时，利用第k和k+1次迭代时得到的图像x_k和x_k+1对迭代的图像进行修正，加快迭代重建的收敛速度。修正的方法是

第k项因子t_k和第k+1项因子t_k+1满足关系式其中初始的第k项的因子t₀＝1。公式(7)中的符号表示将公式(7)中右边的结果赋值给公式(7)左边的x_k+1中。

迭代过程中，当相邻两次的重建图像x_k和x_k+1的L2范数误差||x_k+1-x_k||₂/||x_k||₂小于10^-4时重建图像稳定，停止迭代运算。将最后一步迭代的重建结果x_k+1输出为最终的重建结果图像。

图3显示了迭代过程中的人脑图像和全采样人脑图像的L2范数误差，在计算时间15s后，L2范数误差的值趋于稳定，重建结果很快地收敛。本实施例重建出的人脑图像如图4所示。

Claims (1)

1.一种磁共振图像的迭代重建方法，其特征在于包括以下步骤：

1)提出解空间受限的稀疏重建模型：测量到的欠采样的k空间数据y表示为y＝UFx，其中x表示待重建的磁共振图像，F表示对待重建的磁共振图像x进行傅里叶变换，Fx就是待重建磁共振图像的完整的k空间数据，U表示对完整的k空间数据进行欠采样；用ψ表示任意的框架变换，ψ^*表示框架变换ψ的复共轭，Φ表示与框架变换ψ对应的对偶框架变换，对偶框架变换Φ通过矩阵Ψ*Ψ的逆(Ψ*Ψ)^-1乘以Ψ^*来表示，也就是Φ＝(Ψ*Ψ)^-1Ψ^*，其中符号“-1”表示矩阵求逆；将N个像素的待重建的磁共振图像x表示成复数域中的一个向量待重建的磁共振图像x在框架变换Ψ的变换域系数α构成的空间记为Range(Ψ)，也就是提出的解空间受限的稀疏重建模型为：

$\underset{\mathrm{\α}\∈Range\left(\mathrm{\Ψ}\right)}{min}\mathrm{\λ}\left|\right|\mathrm{\α}|{|}_1+\frac{1}{2}||y-UF\mathrm{\Φ}\mathrm{\α}|{|}_2^2$

其中正则化参数λ用于权衡||Ψx||₁和两项的重要性，α∈Range(Ψ)表示变换域系数α是系数空间Range(Ψ)中的一个向量；

2)重建磁共振图像：解空间受限的稀疏重建模型(1)中的变换域系数α通过公式(2)的迭代运算来求解：

α_k+1＝T_γλ(ΨΦα_k+γΦ^*F^*U^T(y-UFΦα_k)) (2)

其中α_k和α_k+1分别是在第k和k+1次迭代时得到的变换域系数，U^T表示U的转置，F^*表示反傅里叶变换，γ是步长，其取值范围为γ＞0；步长γ与正则化参数λ的乘积为γλ，T_γλ表示对给定向量z＝(ΨΦα_k+γΦ^*F^*U^T(y-UFΦα_k))的所有元素进行软阈值操作，软阈值操作的定义为：给定向量z，若向量z的第j个元素z_j的绝对值|z_j|≤γλ，则z_j＝0；若|z_j|＞γλ，则z_j＝sgn(z_j)(|z_j|-γλ)，其中sgn(z_j)是符号函数；

由于待重建的磁共振图像x用对偶框架Φ与变换域系数α的乘积来表示，即x＝Φα，因此待重建的磁共振图像x按照公式(3)的迭代运算获得

x_k+1＝ΦT_γλ(Ψx_k+γΦ^*F^*U^T(y-UFx_k)) (3)

其中x_k和x_k+1分别是在第k和k+1迭代时得到的磁共振图像；

3)加快图像重建的速度：引入一组因子{t_k}，其中k是大于等于0的整数，初始化的因子t₀＝1，已知这组因子中的第k项因子t_k，通过表达式可以计算第k+1项因子t_k+1；在第k+1次迭代重建图像后，利用第k和k+1迭代时分别得到的磁共振图像x_k和x_k+1对迭代的磁共振图像进行修正，从而加快磁共振图像重建的收敛速度，修正的方法是：

$x_{k+1}\⇐x_{k+1}+\left(\frac{t_k-1}{t_{k+1}}\right)(x_{k+1}-x_k)---\left(4\right)$

公式(4)中的符号表示将公式(4)中右边的结果赋值给公式(4)左边的x_k+1中；

迭代过程中，当相邻两次的重建图像x_k和x_k+1的L2范数误差||x_k+1-x_k||₂/||x_k||₂缩小到指定的范围时，停止迭代运算；将最后一步迭代的重建结果x_k+1输出为最终重建的磁共振图像。