CN106485764A - Mri图像的快速精确重建方法 - Google Patents

Abstract

本发明公开了一种MRI图像的快速精确重建方法，包括：在目标视野中进行扫描得到K空间全采样数据；利用径向采样对所述K空间全采样数据进行采样得到欠采样数据；对欠采样数据进行多线圈数据恢复得MRI图像；利用双密度双树复小波变换作为稀疏基，根据MRI图像在该双密度双树复小波变换下的先验信息，与FISTA算法相结合进行图像重构。该重建方法，能够更好地刻画图像细节以及信息，提高恢复图像的质量。

Description

MRI图像的快速精确重建方法

技术领域

本发明涉及医学图像处理技术领域，尤其涉及一种MRI图像的快速精确重建方法。

背景技术

磁共振成像(Magnetic Resonance Imaging，MRI)无电离辐射，具有极高的软组织对比度以及空间分辨率，利用氢核在外界磁场作用下的共振原理进行成像，已得到医学检测中得到广泛应用。然而，MRI由于硬件扫描系统和传统奈奎斯特采样定理的限制，磁共振成像在K空间的采样时间过长，成像速度慢，且其较长的成像时间会引入运动伪影和原始图像的对比度失真，制约其在临床医学中发展。

2004年，Candes等提出的压缩感知(Compressed sensing)理论突破了传统奈奎斯特采样定理的限制，它充分运用了信号本身或者在变换域上可以稀疏表示这一先验信息，利用随机投影实现了在远低于奈奎斯特频率的采样频率下直接采样少量的数据点，在满足Tao等提出的受限等距性条件下，进而利用非线性重建算法可以恢复出原始信号。但是，MRI由于硬件采样系统限制，不能做到完全的随机采样。

2007年，Lustig等人利用变密度随机采样，选取传统的离散小波变换进行图像的稀疏变换，利用L1范数最小化的非线性重建方法，首次将压缩感知应用到MRI领域。自此，为了提高MRI图像恢复质量，研究者大都在K空间采样轨迹，稀疏基，非线性重建算法三个方面进行研究拓展。

K空间采样轨迹是实现采样的具体方式，分为笛卡尔采样和非笛卡尔采样。现有成像设备的MRI脉冲序列大都是基于笛卡尔采样的。非笛卡尔采样可分为径向采样，螺旋采样，高斯随机采样。相比于笛卡尔采样，非笛卡尔采样更满足不相干要求，支持高度的欠采样。

针对MRI图像选择合适的稀疏基是现有磁共振压缩感知CS-MRI重建算法的基础和前提，合适的稀疏基有利于提高信号的采集速度，减少存储信号所占用的资源。传统的离散小波变换具有对输入数据敏感缺少平移不变性，方向选择性差，在±45°方向有混叠现象以及没有相位空间信息等缺点，恢复MRI图像时，不能更好地刻画图像的边缘成分，缺少图像的轮廓信息。轮廓波作为稀疏基时，它只能很好地表示曲线信息而不能描述图像中奇异点的信息。双树复小波提供±15°，±45°和±75°等6个主方向上的信息描述，但对于一些复杂的图像，由于其有限的方向选择性的限制，在边缘有衍生物产生，不能精确描述原始图像的细节信息。

此外，现有的恢复CS-MRI重建算法共轭梯度法，迭代收缩阈值法，两步迭代收缩阈值法，快速迭代收缩阈值法(fast iterative shrinkage/soft algorithm，FISTA)等，虽然能够恢复图像，但是采样时间长，且都是以小波作为稀疏基，不能描述原始图像细节信息。

鉴于此，有必要进行深入研究，以实现MRI图像的快速精确重建。

发明内容

本发明的目的是提供一种MRI图像的快速精确重建方法，能够更好地刻画图像细节以及信息，提高恢复图像的质量。

本发明的目的是通过以下技术方案实现的：

一种MRI图像的快速精确重建方法，包括：

在目标视野中进行扫描得到K空间全采样数据；

利用径向采样对所述K空间全采样数据进行采样得到欠采样数据；

对欠采样数据进行多线圈数据恢复得MRI图像；

利用双密度双树复小波变换作为稀疏基，根据MRI图像在该双密度双树复小波变换下的先验信息，与FISTA算法相结合进行图像重构。

所述在目标视野中进行扫描得到K空间全采样数据的步骤包括：

选择脉冲扫描序列，并设置扫描参数后对目标视野中的扫描对象进行扫描；

收集由所述扫描对象受脉冲序列激发所产生的磁共振信号，获取K空间全采样数据。

所述利用径向采样对所述K空间全采样数据进行采样得到欠采样数据包括：

利用径向采样方式根据K空间全采样数据能量分布规律，生成径向采样轨迹，从而得到欠采样数据；

所述径向采样是沿辐射式方式采集K空间全采样数据，所采样得到的数据中条线含有等量的低频到高频的信息。

所述对欠采样数据y进行压缩感知初始恢复得MRI图像的步骤包括：

构建多线圈数据重建模型：

其中，m为线圈编号，M为总的线圈数，I_m(i,j)为第m个线圈的图像，为第m个线圈的噪声方差，I(i,j)为重建的二维MRI图像；

对每个线圈的采样数据利用逆傅里叶变换，单独重建得到上述1)中的I_m(i,j)；

对上述2)中获得的所有线圈I_m(i,j)用其平方和的平方根获得最终的二维MRI图像I(i,j)。

所述利用双密度双树复小波变换作为稀疏基，根据MRI图像x在该双密度双树复小波变换下的先验信息，与FISTA算法相结合进行图像重构的步骤包括：

步骤a、构建基于双密度双数复小波的压缩感知MRI图像重建的优化数学模型：

F(x)＝f(x)+g(x)

其中，f(x)是光滑连续凸函数，g(x)是连续凸函数；

上式中，Φ_i表示双密度双树复小波变换的N个高频子带系数构成的小波基，i＝1,2,…,N，R为欠采样的傅里叶变换，τ为可调正则化参数，g(x)为惩罚项，Φ_ix要尽量稀疏，为待求解的F(x)的最小值，也即最终重建结果；

步骤b、设置如下参数：最大迭代次数maxIter，迭代次数k＝1,2,…,maxIter，迭代系数t＝1，，初始阈值threshold，给定k＝1时近似函数起始图像z₁＝x₀＝I(i,j)，I(i,j)为恢复得的MRI图像；

步骤c、获取连续可微函数f(x)在z_k处的梯度为

其中，z_k为近似函数起始图像；y为欠采样数据；R^T为欠采样的傅里叶逆变换；

步骤d、利用z_k沿着负梯度方向下降，获取局部优化图像x_g：

x_g＝z_k-ρR^T(Rz_k-y)；

其中，ρ为固定下降步长；

步骤e、对大小为M*N的MRI局部优化图像x_g，以双密度双树复小波变换为稀疏基，进行L层双密度双树复小波变换，得到各个尺度上的高频子带系数w_h；

w_h＝τ||Φ_ix_g||₁,i＝1,2,…,N

步骤f、由于L1范数的可分离性，使用收缩阈值算法进行计算，对各个尺度上的高频子带系数w_h进行修正：

w_h'＝shrink(w_h,threshold)；

其中shrink为收缩阈值算子，阈值选取为threshold，w_h'为修正后的高频子带系数；

步骤g、利用步骤f中修正后的高频子带系数w_h'作双密度双树复小波反变换，获取第k次迭代的图像x_k：

x_k＝Φ_i ^-1(w_h')，i＝1,2,…,N

步骤h、更新阈值，进行收缩阈值松弛：

threshold'＝threshold*β

其中，β为松弛系数，当松弛后的阈值threshold'＜ε时，截断迭代，，输出图像x_k作为最终重建结果否则，判断k是否达到最大迭代次数maxIter，若是，则输出图像x_k作为最终重建结果若否，则令threshold＝threshold'，并转入步骤i；

步骤i、更新迭代系数：

步骤j、第k+1次迭代使用前两次迭代过程的结果x_k，x_k-1，以及经过迭代系数t_k与步骤i所获得的t_k+1的线性组合相应的近似函数起始图像z_k+1，再转入步骤c：

双树双密度复小波是基于两个不同的尺度函数 和四个不同的小波函数 构成；其中一对小波函数构成希尔伯特变换对是一种基2冗余且不随尺度变化的紧框架变换；

且双密度双树复小波对输入图像并行使用4个过采样的二维双密度滤波器组实现，先对行滤波再对列滤波，并且行和列采用不同的滤波器组，每层获得4个低频子带图和32个高频子带图，共描述16个不同的主方向信息。

由上述本发明提供的技术方案可以看出，采用径向采样，对于伪影具有更低的敏感度，并且能够高度欠采样，减少采样时间；同时，采用双密度双树复小波变换作为稀疏基，结合FISTA算法，收敛速度快且可以准确捕捉MRI图像的细节及边缘轮廓信息，能更好的描述图像的方向属性。

附图说明

为了更清楚地说明本发明实施例的技术方案，下面将对实施例描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域的普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他附图。

图1为本发明实施例提供的一种MRI图像的快速精确重建方法的流程图。

图2为本发明实施例提供的在目标视野中进行扫描得到K空间全采样数据的流程图；

图3为本发明实施例提供的采样率为20％的径向采样生成的径向采样轨迹示意图；

图4为本发明实施例提供的二维二层双密度双树复小波变换的滤波器组结构示意图；

图5为本发明实施例提供的本发明方法和传统方法对MRI脑部图像image1的重建效果对比图；

图6为本发明实施例提供的本发明方法和传统方法对MRI脑部图像image2的重建效果对比图；

图7为本发明实施例提供的本发明方法和传统方法为MRI颈椎图像image3的重建效果对比图。

具体实施方式

下面结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明的保护范围。

图1为本发明实施例提供的一种MRI图像的快速精确重建方法的流程图。如图1所示，其主要包括如下步骤：

步骤11、在目标视野中进行扫描得到K空间全采样数据。

步骤12、利用径向采样对所述K空间全采样数据进行采样得到欠采样数据。

步骤13、对欠采样数据进行多线圈数据重建恢复得MRI图像。

步骤14、利用双密度双树复小波变换作为稀疏基，根据MRI图像在该双密度双树复小波变换下的先验信息，与FISTA算法相结合进行图像重构。

本发明实施例的上述方案，是为了提高恢复磁共振图像的细节信息，加快磁共振成像的速度。该方案继承了FISTA算法收敛速度快的优点，采用双密度双树复小波作为稀疏基，能够更好地刻画图像细节以及信息，提高恢复图像的质量。

为了便于理解，下面针对上述方案的各个步骤做详细介绍。

1、上述步骤11中描述了，在目标视野中进行扫描得到K空间全采样数据，该步骤可以通过如图2所示方式实现：

1)选择脉冲扫描序列，并设置扫描参数后对目标视野中的扫描对象进行扫描。

2)收集由所述扫描对象受脉冲序列激发所产生的磁共振信号，获取K空间全采样数据。

本领域技术人员可以理解，上述过程所选择的脉冲扫描序列，以及设置的扫描参数可以根据需求或者经验来选定。

2、上述步骤12中描述了，利用径向采样对所述K空间全采样数据进行采样得到欠采样数据该步骤可以通过如下方式实现：

利用径向采样方式根据K空间全采样数据能量分布规律(K空间全采样数据中心区域包含图像大量信息采样密度高，四周区域信息较少采样密度低)，生成径向采样轨迹(如图3所示)，从而得到欠采样数据；

所述径向采样是沿辐射式方式采集K空间全采样数据，所采样得到的数据中条线含有等量的低频到高频的信息。

3、上述步骤13中描述了，对欠采样数据y进行多线圈数据初始重建得MRI图像，该步骤可以通过如下所示方式实现：

1)构建多线圈数据重建模型：

其中，m为线圈编号，M为总的线圈数，I_m(i,j)为第m个线圈的图像，为第m个线圈的噪声方差，I(i,j)为重建的二维MRI图像；

2)对每个线圈的采样数据利用逆傅里叶变换，单独重建得到上述1)中的I_m(i,j)。

3)对上述2)中获得的所有线圈I_m(i,j)用其平方和的平方根获得最终的二维MRI图像I(i,j)。

4、上述步骤14中描述了，利用双密度双树复小波变换作为稀疏基，根据MRI图像x在该双密度双树复小波变换下的先验信息，与FISTA算法相结合进行图像重构，该步骤可以通过如下所示方式实现：

步骤a、构建基于双密度双数复小波的压缩感知MRI图像重建的优化数学模型：

F(x)＝f(x)+g(x)

其中，f(x)是光滑连续凸函数，g(x)是连续凸函数，不一定光滑。

上式中，Φ_i表示双密度双树复小波变换的N个高频子带系数构成的小波基，i＝1,2,…,N，R为欠采样的傅里叶变换，τ为可调正则化参数，g(x)为惩罚项，Φ_ix要尽量稀疏，即非零元素个数要尽可能少，为待求解的F(x)的最小值，也即最终重建结果。

步骤b、设置如下参数：最大迭代次数maxIter,迭代次数k＝1,2,…,maxIter,迭代系数t＝1，初始阈值threshold，给定k＝1时近似函数起始图像z₁＝x₀，x₀＝I(i,j)。

步骤c、获取连续可微函数f(x)在z_k处的梯度为

其中，z_k为近似函数起始图像；R^T为欠采样的傅里叶逆变换；y为欠采样数据。

步骤d、利用z_k沿着负梯度方向下降(此时下降最快)，获取局部优化图像x_g：

x_g＝z_k-ρR^T(Rz_k-y)；

其中，ρ为固定下降步长；

步骤e、对大小为M*N的MRI局部优化图像x_g，以双密度双树复小波变换为稀疏基，进行L层双密度双树复小波变换，得到各个尺度上的高频子带系数w_h；

w_h＝τ||Φ_ix_g||₁,i＝1,2,…,N

步骤f、由于L1范数的可分离性，使用收缩阈值算法进行计算，对各个尺度上的高频子带系数w_h进行修正：

w_h'＝shrink(w_h,threshold)；

其中，shrink为收缩阈值算子，阈值选取为threshold，w_h'为修正后的高频子带系数。

步骤g、利用步骤f中修正后的高频子带系数w_h'作双密度双树复小波反变换，获取第k次迭代的图像x_k：

x_k＝Φ_i ^-1(w_h')，i＝1,2,…,N

步骤h、更新阈值，进行收缩阈值松弛：

threshold'＝threshold*β；

其中，β为松弛系数，当松弛后的阈值threshold'＜ε时，截断迭代，输出图像x_k作为最终重建结果否则，判断k是否达到最大迭代次数maxIter，若是，则输出图像x_k作为最终重建结果若否，则令threshold＝threshold'，并转入步骤i。

步骤i、更新迭代系数：

步骤j、第k+1次迭代使用前两次迭代过程的结果x_k，x_k-1，以及经过迭代系数t_k与步骤i所获得的t_k+1的线性组合生成相应的近似函数起始图像z_k+1，再转入步骤c：

示例性的，上述步骤所设定的参数可以采用如下数值：最大迭代次数maxIter＝200；迭代系数t＝1；初始阈值threshold＝0.015；N＝16；ρ＝1；L＝3，4，5；β＝0.9；ε＝0.001，均加入标准差为0.01的高斯白噪声。需要说明的是，上述具体数值仅为举例并非构成限制，在实际应用中，用户可根据实际情况或者经验来改变上述参数的具体数值。

本发明实施例中，所述双树双密度复小波是基于两个不同的尺度函数 和四个不同的小波函数 构成；其中一对小波函数构成希尔伯特变换对是一种基2冗余且不随尺度变化的紧框架变换；且双密度双树复小波对输入图像并行使用4个过采样的二维双密度滤波器组实现，先对行滤波再对列滤波，并且行和列采用不同的滤波器组，每层获得4个低频子带图和32个高频子带图，共描述16个不同的主方向信息。二维二层双密度双树复小波的滤波器组结构如图4所示。Lp和Hp构成了第一层的滤波器组结构，分别代表1个尺度滤波器和8个小波滤波器，Lo,Hi表示第二层及其以后各层的分解滤波器。

所述双密度双树复小波具有近似平移不变性、良好的方向选择性，有限的冗余度，更加近似连续的小波变换。它的16个主方向更能准确地描述复杂MRI图像的细节，能反应出图像在不同分辨率上沿16个方向的变化情况，能更好的描述图像的方向属性。

为了验证本发明的可行性和有效性，利用Matlab平台进行仿真实验，使用GEMR750仪器的真实采集K空间数据，通过仿真结果可以证明本发明方法相对于传统的以小波基作为稀疏基的MRI图像恢复方法，更能具体描绘图像的细节：

1、仿真条件

1)仿真实验采用模拟的径向采样作为采样方式

2)仿真实验所用编程平台为MatlabR2015b

3)仿真实验的实验对象分别采用离散小波变换结合FISTA算法(FISTA-DWT)，双树双密度复小波变换结合FISTA算法(FISTA-CDDDT)。

2、仿真实验所用MRI脑部图像image1，MRI脑部图像image2，MRI颈椎图像image3来自志愿者在GE MR750仪器扫描的扫描数据。

本发明采用的恢复图像质量的评价指标为：

1)相对误差(Relative error)，表示为：

式中，x表示原始图像，表示重建图像，相对误差越小，反映了重建图像与原始图像之间的相似度越高。

2)峰值信噪比(Peak signal-to-noise ratio,PSNR)

式中，M和N分别表示图像行和列的大小，MAX表示输入图像最大值，峰值信噪比越高越好。

3、仿真实验结果：

表1恢复图像image1，image2，image3的PSNR(单位dB)和Relative error值

从表1中可以看出，本发明重建方法FISTA-CDDDT比已有重建方法FISTA-DWT，在采样率为20％时对image1和image2以及在采样率为20％时对image3，进行图像重建时，平均要高3dB。

利用本发明重建方法和已有重建方法FISTA-DWT对MRI脑部图像image1对比图如图5所示。其中图5(a)为MRI脑部图像image1；图5(b)为采用FISTA-DWT得到的重建图像；图5(c)为采用本发明方法得到的重建图像；图5(d)，(e)，(f)分别是图5(a)，(b)，(c)中白色方框的局部细节放大图。

利用本发明重建方法和已有重建方法FISTA-DWT对MRI脑部图像image2对比图如图6所示。其中图6(a)为MRI真实图像image2；图6(b)为采用FISTA-DWT得到的重建图像；图6(c)右图为采用本发明方法得到的重建图像；图6(d)，(e)，(f)分别是图5(a)，(b)，(c)中白色方框的局部细节放大图。

利用本发明重建方法和已有重建方法FISTA-DWT对MRI颈椎图像image3对比图如图7所示。其中图7(a)为MRI真实图像image3；图7(b)为采用FISTA-DWT得到的重建图像；图7(c)为采用本发明方法得到的重建图像。

参照图5，图6和图7可见，采用FISTA-DWT方法得到的恢复图像，图5(e)，图6(e)以及图7(e)产生了伪影，对中心细节恢复效果差。采用本发明提出的方法FISTA-CDDDT：图5(f)，图6(f)和图7(f)恢复图像清晰，准确刻画了细节和边缘信息，减少伪影，恢复图像的均匀区域较为平滑，恢复图像质量更好。

综上所述，本发明方法相较于现有FISTA-DWT重建方法，对MRI仿体(加有琼脂和明胶，水脂含量约为7:3)，MRI脑部图像和MRI颈椎图像有良好的重建效果，在能够快速重建实现MRI图像同时保证恢复图像的细节，减少伪影，这是一种可行的MRI快速精确重建算法。

以上所述，仅为本发明较佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明披露的技术范围内，可轻易想到的变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应该以权利要求书的保护范围为准。

Claims (6)

1.一种MRI图像的快速精确重建方法，其特征在于，包括：

在目标视野中进行扫描得到K空间全采样数据；

利用径向采样对所述K空间全采样数据进行采样得到欠采样数据；

对欠采样数据进行多线圈数据恢复得MRI图像；

利用双密度双树复小波变换作为稀疏基，根据MRI图像在该双密度双树复小波变换下的先验信息，与FISTA算法相结合进行图像重构。

2.根据权利要求1所述的一种压缩感知MRI图像重建方法，其特征在于，所述在目标视野中进行扫描得到K空间全采样数据的步骤包括：

选择脉冲扫描序列，并设置扫描参数后对目标视野中的扫描对象进行扫描；

收集由所述扫描对象受脉冲序列激发所产生的磁共振信号，获取K空间全采样数据。

3.根据权利要求1所述的一种压缩感知MRI图像重建方法，其特征在于，所述利用径向采样对所述K空间全采样数据进行采样得到欠采样数据包括：

利用径向采样方式根据K空间全采样数据能量分布规律，生成径向采样轨迹，从而得到欠采样数据；

所述径向采样是沿辐射式方式采集K空间全采样数据，所采样得到的数据中条线含有等量的低频到高频的信息。

4.根据权利要求1所述的一种压缩感知MRI图像重建方法，其特征在于，所述对欠采样数据y进行压缩感知初始恢复得MRI图像的步骤包括：

构建多线圈数据重建模型：

$I(i,j)=\sqrt{\underset{m=1}{\overset{M}{\Σ}}\frac{\left|I_m(i,j)\right|}{{\mathrm{\δ}}_m^2}};$

其中，m为线圈编号，M为总的线圈数，I_m(i,j)为第m个线圈的图像，为第m个线圈的噪声方差，I(i,j)为重建的二维MRI图像；

对每个线圈的采样数据利用逆傅里叶变换，单独重建得到上述1)中的I_m(i,j)；

对上述2)中获得的所有线圈I_m(i,j)用其平方和的平方根获得最终的二维MRI图像I(i,j)。

5.根据权利要求1或4所述的一种压缩感知MRI图像重建方法，其特征在于，所述利用双密度双树复小波变换作为稀疏基，根据MRI图像x在该双密度双树复小波变换下的先验信息，与FISTA算法相结合进行图像重构的步骤包括：

步骤a、构建基于双密度双数复小波的压缩感知MRI图像重建的优化数学模型：

F(x)＝f(x)+g(x)

其中，f(x)是光滑连续凸函数，g(x)是连续凸函数；

$f\left(x\right)=\frac{1}{2}\left|\right|Rx-y|{|}_2^2,g\left(x\right)=\mathrm{\τ}|\left|{\mathrm{\Φ}}_{i}x\right|{|}_1,i=1,2,\mathrm{...},N$

$\hat{x}=\underset{x}{\arg \min }\left\{\frac{1}{2}\right||Rx-y|{|}_2^2+\mathrm{\τ}\left|\right|{\mathrm{\Φ}}_{i}x|{|}_1,i=1,2,\mathrm{...},N\}$

上式中，Φ_i表示双密度双树复小波变换的N个高频子带系数构成的小波基，i＝1,2,…,N，R为欠采样的傅里叶变换，τ为可调正则化参数，g(x)为惩罚项，Φ_ix要尽量稀疏，为待求解的F(x)的最小值，也即最终重建结果；

步骤b、设置如下参数：最大迭代次数maxIter，迭代次数k＝1,2,…,maxIter，迭代系数t＝1，，初始阈值threshold，给定k＝1时近似函数起始图像z₁＝x₀＝I(i,j)，I(i,j)为恢复得的MRI图像；

步骤c、获取连续可微函数f(x)在z_k处的梯度为▽f(z_k)：

▽f(z_k)＝R^T(Rz_k-y)

其中，z_k为近似函数起始图像；y为欠采样数据；R^T为欠采样的傅里叶逆变换；

步骤d、利用z_k沿着负梯度方向下降，获取局部优化图像x_g：

x_g＝z_k-ρR^T(Rz_k-y)；

其中，ρ为固定下降步长；

步骤e、对大小为M*N的MRI局部优化图像x_g，以双密度双树复小波变换为稀疏基，进行L层双密度双树复小波变换，得到各个尺度上的高频子带系数w_h；

w_h＝τ||Φ_ix_g||₁,i＝1,2,…,N

步骤f、由于L1范数的可分离性，使用收缩阈值算法进行计算，对各个尺度上的高频子带系数w_h进行修正：

$shrink(w_h,threshold)=\frac{w_h}{\left|w_h\right|}*max\left(\right|w_h|-threshold,0);$

w_h'＝shrink(w_h,threshold)；

其中shrink为收缩阈值算子，阈值选取为threshold，w_h'为修正后的高频子带系数；

步骤g、利用步骤f中修正后的高频子带系数w_h'作双密度双树复小波反变换，获取第k次迭代的图像x_k：

x_k＝Φ_i ^-1(w_h')，i＝1,2,…,N

步骤h、更新阈值，进行收缩阈值松弛：

threshold'＝threshold*β

其中，β为松弛系数，当松弛后的阈值threshold'＜ε时，截断迭代，，输出图像x_k作为最终重建结果否则，判断k是否达到最大迭代次数maxIter，若是，则输出图像x_k作为最终重建结果若否，则令threshold＝threshold'，并转入步骤i；

步骤i、更新迭代系数：

$t_{k+1}=\frac{1+\sqrt{1+4t_k^2}}{2};$

步骤j、第k+1次迭代使用前两次迭代过程的结果x_k，x_k-1，以及经过迭代系数t_k与步骤i所获得的t_k+1的线性组合相应的近似函数起始图像z_k+1，再转入步骤c：

$z_{k+1}=x_k+\frac{t_k-1}{t_{k+1}}(x_k-x_{k-1}).$

6.根据权利要求1所述的一种压缩感知MRI图像重建方法，其特征在于，

双树双密度复小波是基于两个不同的尺度函数和四个不同的小波函数p＝1,2构成；其中一对小波函数构成希尔伯特变换对是一种基2冗余且不随尺度变化的紧框架变换；

且双密度双树复小波对输入图像并行使用4个过采样的二维双密度滤波器组实现，先对行滤波再对列滤波，并且行和列采用不同的滤波器组，每层获得4个低频子带图和32个高频子带图，共描述16个不同的主方向信息。