CN108416819B - 一种基于curvelet-fista的压缩采样磁共振图像重建方法 - Google Patents

Abstract

一种基于curvelet‑fista的压缩采样磁共振图像重建方法，其特征在于包括如下步骤：步骤(1)初始化：步骤(2)Curvelet稀疏变换：步骤(3)对图像进行快速迭代收缩阈值处理，利用快速迭代收缩阈值算法计算系数矩阵Θ_k。本发明提出一种基于Curvelet变换的MR图像压缩感知重建方法，利用Curvelet稀疏表示的优越特性，对图像平滑部分和边缘部分进行表达。对应边缘的大系数能起到很好的能量集中效果，并且采用快速迭代收缩阈值进行截断，使得在重建过程中得到了边缘及细节部分重建效果的提升，图像整体重建质量也得到了提升。

Description

一种基于curvelet-fista的压缩采样磁共振图像重建方法

技术领域

本发明涉及一种基于curvelet-fista的压缩采样磁共振图像重建方法。

背景技术

随着软硬件技术的进步，磁共振成像^[1](Magnetic resonance imaging,MRI)在临床上的应用日益广泛，相比于计算机断层扫描^[2](Computed Tomography,CT)技术的优势是对人体没有电离辐射等伤害，而且可以对人体的各种软组织进行成像。较慢的数据扫描速度一直是制约其进一步发展的关键问题，由此而引发的其他问题也有很多。因此，合理的采用快速采集技术不仅可以缩短MRI的检查时间，还可以提高检查质量。

压缩感知^[3](Compressed Sensing，CS)论作为一个新兴的信号采集、处理领域理论，突破了奈奎斯特(Nyquist)采样定理中对采样频率的限制。它指出，只要信号是可压缩的或在某个变换域是稀疏的，那么就可以用一个与变换基不相关的观测矩阵^[4]将变换所得高维信号投影到一个低维空间上，然后通过求解一个优化问题就可以从这些少量的投影中以高概率重建出原信号，可以证明这样的投影包含了重建信号的足够信息。

目前，基于CS理论的MRI技术尚处于理论研究阶段，很多关键问题需要解决。现阶段，压缩感知主要采用正交小波进行图像的稀疏表示^[5]。小波是表示具有点奇异性目标函数的最优基，一般的小波函数^[6]是变量可分离的二元函数，这样通过一维多尺度分析导出二维多尺度分析，进而导出二维小波空间和小波函数，常用的二维小波是一维正交小波的张量积，仅有水平、垂直和对角面3个方向，且各向同性，轮廓、边缘和纹理更高维的几何特征表达能力有限。为了进一步表示多维信号中更为普遍的曲面型奇异性，发展出了Curvelet变换^[7]。与小波变换不同，Curvelet包含尺度、位移和方向3个参量，具有更好的方向辨识能力，对图像的边缘，如曲面和直面的几何特征的表达更优于小波变换。

Curvelet变换是基于傅里叶变换和小波变换的一种改进，其特点是有高度的各向异性，具有良好表达图形沿边缘的信息的能力，对于恢复形状的沿边缘的主要结构和抑制周边噪声有其特有优势。

参考文献

[1]Eutsler E P,Khanna G.Whole-body magnetic resonance imaging inchildren:technique and clinical applications[J].Pediatric Radiology,2016,46(6):858.

[2]Jarrett C.Computed Tomography:Fundamentals,System Technology,ImageQuality,Applications[Book Review][J].IEEE Engineering in Medicine&BiologyMagazine,2007,26(2):12-0.

[3]Blaszczyk L.Compressed sensing in MRI–mathematical preliminariesand basic examples[J].Nukleonika,2016,61(1):41-43.

[4]Wang B,Ma S X.Improvement of Gaussian Random Measurement Matricesin Compressed Sensing[J].Advanced Materials Research,2011,301-303:245-250.

[5]Caballero J,Rueckert D,Hajnal J V.Dictionary Learning and TimeSparsity in Dynamic MRI[M].Medical Image Computing and Computer-AssistedIntervention-MICCAI2012.Springer Berlin Heidelberg,2012:256-263.

[6]Jiang W,Yang J J.A New Wavelet-Based Compressive Sensing for ImageCompression[J].Advanced Materials Research,2013,756-759:1684-1690.

[7]Emmanuel Candès,Laurent Demanet,David Donoho,et al.Fast DiscreteCurvelet Transforms[J].Multiscale Modeling&Simulation,2006,5(3):861--899.

[8]Wang S.An improved reconstruction method for CS-MRI based onexponential wavelet transform and iterative shrinkage/thresholding algorithm[J].Journal of Electromagnetic Waves&Applications,2014,28(18):2327-2338.

[9]王开,刘郁林,和继威.准Toeplitz测量矩阵的有限等距性质分析[J].计算机应用研究,2011,28(4):1512-1514.

[10]方晟,吴文川,应葵,等.基于非均匀螺旋线数据和布雷格曼迭代的快速磁共振成像方法[J].物理学报,2013,62(4):502-508.

[11]Lustig M,Donoho D J.Sparse MRI:The application of compressedsensing for rapid MR imaging[J].Journal of the Chemical Fertilizer Industry,2007,58(6):1182.

[12]Jiao L C,Tan S,Liu F.Ridgelet Theory:from Ridgelet Transform toCurvelet[J].Chinese Journal of Engineering Mathematics,2005,22(5):761-773.

[13]Hammouche A M,El-Bakry H M,Mostafa R R.Image Contrast EnhancementUsing Fast Discrete Curvelet Transform via Unequally Spaced Fast FourierTransform(FDCT-USFFT)[J].International Journal of Electronics Communicationand Computer Engineering.2016,7(2):88-93.

[14]Chambolle A,Dossal C.On the Convergence of the Iterates of the“Fast Iterative Shrinkage/Thresholding Algorithm”[J].Journal of OptimizationTheory&Applications,2015,166(3):968-982.

[15]Wang Y,Lian Q S,Li K.MRI reconstruction based on compoundregularizers and compressed sensing[J].Optical Technique,2010,36(3):350-355.

[16]裴营.基于曲波变换的图像压缩算法研究[D].天津师范大学,2009:1-50.

[17]Tan H,Meyer C H.Estimation of k-space trajectories in spiral MRI[J].Magnetic Resonance in Medicine,2009,61(6):1396–1404.

[18]Ravishankar S,Bresler Y.MR Image Reconstruction From HighlyUndersampled k-Space Data by Dictionary Learning[J].IEEE Transactions onMedical Imaging,2011,30(5):1028-1041.

[19]王斌.MATLAB实现数字图像增强处理[J].佳木斯大学学报(自然科学版),2005,23(1):31-34.

[20]李国燕,侯向丹,周博君,等.基于离散剪切波的压缩感知MRI图像重建[J].计算机应用研究,2013,30(6):1895-1898.

发明内容

为解决压缩感知MRI图像重建问题，克服目前压缩感知算法中存在的小波变换方向选择性差、细节重建效果不好的缺点^[16]，本发明提供一种基于curvelet-fista的压缩采样磁共振图像重建方法，本发明将Curvelet变换融入快速迭代收缩阈值算法中，与快速迭代收缩阈值算法^[8](Fast Iterative Shrinkage-Thresholding Algorithm,FISTA)相比，具有更好的恢复效果。通过Curvelet-FISTA算法对MR图像的重建，可以提高图像恢复效果和收敛速度。

为解决上述问题，本发明提供以下技术方案：

一种基于curvelet-fista的压缩采样磁共振图像重建方法，采用了Curvelet变换结合快速迭代收缩阈值算法实现图像重建；基于Curvelet的快速迭代收缩阈值算法，包括如下步骤：

步骤(1)初始化：给定迭代停止阈值ε、图像分块观测矩阵Φ、图像观测矩阵Y，设定迭代计数器、Curvelet变换自适应阈值收缩算子及分别为k＝1、t＝0、λ＝λ₀、μ＝0.9；

步骤(2)Curvelet稀疏变换：

Θ_k+1＝C(X_k+1)

(5-1)

其中，C为Curvelet稀疏变换；

步骤(3)对图像进行快速迭代收缩阈值处理，利用快速迭代收缩阈值算法计算系数矩阵Θ_k；

步骤(4)更新参数t,利用快速迭代收缩阈值算法计算t_k+1，利用快速迭代收缩阈值算法计算z_k+1；

步骤(5)利用快速迭代收缩阈值算法计算迭代终止函数G_k(Θ_k,Θ_k-1)，如果G_k(Θ_k,Θ_k-1)≤ε，停止迭代，直接执行步骤(8)，如果G_k(Θ_k,Θ_k-1)＞ε，则继续往下执行步骤(6)；

步骤(6)收缩正则化参数λ，λ＝0.9×λ；

步骤(7)更新迭代次数，令k＝k+1，返回到步骤(2)；

步骤(8)对迭代得到的系数矩阵进行Curvelet反变换，得到重建图像X：

X_k+1＝C^T(Θ_k+1) (5-2)

其中，C^T为Curvelet反变换；

步骤(9)显示重建图像X，算法结束。

所述快速迭代收缩阈值算法包括以下步骤：

由拉格朗日约束条件：

所示，

令

g(Θ)＝||Θ||₁；

arg min f(Θ)+λg(Θ) (4-1)

其中，正则化参数λ在式(4-1)中的作用是平衡数据f(Θ)与g(Θ)的比重。求解式(4-1)中f(Θ)的一个简单的方法就是引入梯度法，通过梯度法在每一次迭代中，不断修正Θ，具体如下(a>0)：

这等价于：

用同样的方法求解式(2-7)，迭代计算式如下：

省略掉常数项后，上式可化为：

由于l₁范数是可分离的，||Θ||₁为它的所有元素的绝对值之和，式(4-1)中g(Θ)的求解可以简化为求Θ的每一个最小化问题，这可以通过阈值收缩求得，则式(4-5)可转化为：

其中，shrink是软阈值算子，即：

shrink(x,β)＝sign(x)·max{|x|-β,0}

(4-7)

为加快收敛速度，引入参数t和参数z，结合前2次的迭代值对其进行更新：

其中，Θ_k-1、Θ_k为Θ前2次的迭代值，此外，t的更新公式为：

由式(4-6)求得Θ_k+1，计算过程如下：

Θ_k+1＝shrink((z_k+1-a▽f(z_k+1)),λa)

(4-10)

迭代的终止条件是通过Θ相邻2次迭代值的相对误差来设定的，终止函数定义为：

终止阈值设为ε，如果G_k(Θ_k,Θ_k-1)≤ε时，则停止迭代；

所述步骤(3)对图像进行快速迭代收缩阈值处理，利用式(4-10)计算系数矩阵Θ_k；

所述步骤(4)更新参数t,利用式(4-9)计算t_k+1，根据式(4-8)计算z_k+1；

所述步骤(5)按式(4-11)计算迭代终止函数G_k(Θ_k,Θ_k-1)，如果G_k(Θ_k,Θ_k-1)≤ε，停止迭代，直接执行步骤(8)，如果G_k(Θ_k,Θ_k-1)＞ε。

所述的一种基于curvelet-fista的压缩采样磁共振图像重建方法还包括离散Curvelet变换的步骤：

Curvelet变换是在频域内实现的，采用频域窗函数U来表示Curvelet基函数

的傅里叶变换。定义了窗函数W(r)和角度窗函数V(t)，均满足可容许条件：

设U_j(ω)＝ψ_j(ω₁)V_j(ω)，其中ψ_j(ω₁)＝ψ(2^-jω₁)，ψ(ω₁)是一个带通滤波器，对于θ_j∈[-π/4,π/4)，有：

U_j,l(ω)＝ψ_j(ω₁)V_j(S_θlω)＝U_j(S_θlω)

(3-2)

在构造上

则对于ω₁＞0，根据式(3-1)有：

当-1≤tanθ_l＜1，即

时，有：

又因为：

则根据式(3-6)满足：

通过对称和旋转才做可以将整个笛卡尔平面进行局部化分割，由

可以得到Curvelet基函数在时域的形式，则Curvelet变换表示：

本发明所提出的重建算法采用的是USFFT(Unequally Spaced fast FourierTransform)算法，实现步骤如下：

第一步，对于给定的笛卡尔坐标中的二维图像f[n₁,n₂]进行2DFFT变换，得到二维的频域表示为：

第二步，在频域，对于每一对(j,l)，对

重新采样或插值得到符合Curvelet频域局部化的形状

第三步，将

与拟合窗U_j,l相乘，得到

在拟合窗上的傅里叶变换

这就是在频域的Curvelet变换；

第四步，对每个乘积

做二维逆傅里叶变换，即得到离散Curvelet系数C^D(j,l,k)。

有益效果：

为了提高MR图像的重建效果和缩短MR图像重建时间，本发明提出一种基于curvelet-fista的压缩采样磁共振图像重建方法，本发明是一种基于Curvelet变换的MRI图像快速迭代收缩阈值重建算法。Curvelet变换具有的多尺度、各向奇异性、更高稀疏表示性能等特性。利用Curvelet稀疏表示和快速迭代收缩阈值算法，以此更好地保持重建图像的细节信息并解决信号重建噪声问题。阈值函数的使用促进了图像的稀疏性是压缩感知(CS)图像重建的一个关键因素。快速迭代收缩阈值算法FISTA与迭代收缩阈值算法ISTA相比，不仅保留了原先的简洁的优势，而且提高了最优化梯度的速率，已经证明该Curvelet-FISTA算法可有效恢复完全采样图像从核磁共振成像中的欠采样数据。重建图像的质量以峰值信噪比(PSNR)、均方误差(MSE)、结构相似性度(SSIM)来衡量，与其他方法相比，该方法显示了较好的重建效果和较快的收敛速度。

本发明提出一种基于Curvelet变换的MR图像压缩感知重建方法，利用Curvelet稀疏表示的优越特性，对图像平滑部分和边缘部分进行表达。对应边缘的大系数能起到很好的能量集中效果，并且采用快速迭代收缩阈值进行截断，使得在重建过程中得到了边缘及细节部分重建效果的提升，图像整体重建质量也得到了提升。实验结果表明，本发明的FCISTA算法重建所得图像可满足人眼视觉的需要，相比于传统的MR图像重建算法，其重建图像质量及细节信息有了明显的提升，对于疾病的检查和判断有着积极的作用。本发明能利用更少的有效信息、更集中的数据，以较大概率重建图像，推进图像压缩感知重建研究的发展。

附图说明

图1为压缩感知观测向量的矩阵表示示意图。

图2为实施例的原始lena图。

图3为实施例的原始脑部MR图。

图4为实施例的spiral采样轨迹图。

图5为实施例的不同图像在FCISTA、FISTA算法下重建的比较示意图。

图6为实施例的FCISTA、FISTA算法下的差值图像示意图。

图7为实施例的不同采样率下FCISTA算法与FISTA算法性能比较示意图。

具体实施方式

下面结合说明书附图和实施例，对本发明的具体实施例做进一步详细描述：

如图1至图7，一种基于curvelet-fista的压缩采样磁共振图像重建方法，采用了Curvelet变换结合快速迭代收缩阈值算法实现图像重建；基于Curvelet的快速迭代收缩阈值算法，包括如下步骤：

步骤(1)初始化：给定迭代停止阈值ε、图像分块观测矩阵Φ、图像观测矩阵Y，设定迭代计数器、Curvelet变换自适应阈值收缩算子及分别为k＝1、t＝0、λ＝λ₀、μ＝0.9；

步骤(2)Curvelet稀疏变换：

Θ_k+1＝C(X_k+1)

(5-1)

其中，C为Curvelet稀疏变换；

步骤(3)对图像进行快速迭代收缩阈值处理，利用快速迭代收缩阈值算法计算系数矩阵Θ_k；

步骤(4)更新参数t,利用快速迭代收缩阈值算法计算t_k+1，利用快速迭代收缩阈值算法计算z_k+1；

步骤(5)利用快速迭代收缩阈值算法计算迭代终止函数G_k(Θ_k,Θ_k-1)，如果G_k(Θ_k,Θ_k-1)≤ε，停止迭代，直接执行步骤(8)，如果G_k(Θ_k,Θ_k-1)＞ε，则继续往下执行步骤(6)；

步骤(6)收缩正则化参数λ，λ＝0.9×λ；

步骤(7)更新迭代次数，令k＝k+1，返回到步骤(2)；

步骤(8)对迭代得到的系数矩阵进行Curvelet反变换，得到重建图像X：

X_k+1＝C^T(Θ_k+1) (5-2)

其中，C^T为Curvelet反变换；

步骤(9)显示重建图像X，算法结束。

所述快速迭代收缩阈值算法包括以下步骤：

由拉格朗日约束条件：

所示，

令

g(Θ)＝||Θ||₁；

arg min f(Θ)+λg(Θ) (4-1)

其中，正则化参数λ在式(4-1)中的作用是平衡数据f(Θ)与g(Θ)的比重。求解式(4-1)中f(Θ)的一个简单的方法就是引入梯度法，通过梯度法在每一次迭代中，不断修正Θ，具体如下(a>0)：

这等价于：

用同样的方法求解式(2-7)，迭代计算式如下：

省略掉常数项后，上式可化为：

由于l₁范数是可分离的，||Θ||₁为它的所有元素的绝对值之和，式(4-1)中g(Θ)的求解可以简化为求Θ的每一个最小化问题，这可以通过阈值收缩求得，则式(4-5)可转化为：

其中，shrink是软阈值算子，即：

shrink(x,β)＝sign(x)·max{|x|-β,0} (4-7)

为加快收敛速度，引入参数t和参数z，结合前2次的迭代值对其进行更新：

其中，Θ_k-1、Θ_k为Θ前2次的迭代值，此外，t的更新公式为：

由式(4-6)求得Θ_k+1，计算过程如下：

迭代的终止条件是通过Θ相邻2次迭代值的相对误差来设定的，终止函数定义为：

终止阈值设为ε，如果G_k(Θ_k,Θ_k-1)≤ε时，则停止迭代；

所述步骤(3)对图像进行快速迭代收缩阈值处理，利用式(4-10)计算系数矩阵Θ_k；

所述步骤(4)更新参数t,利用式(4-9)计算t_k+1，根据式(4-8)计算z_k+1；

所述步骤(5)按式(4-11)计算迭代终止函数G_k(Θ_k,Θ_k-1)，如果G_k(Θ_k,Θ_k-1)≤ε，停止迭代，直接执行步骤(8)，如果G_k(Θ_k,Θ_k-1)＞ε。

所述的一种基于curvelet-fista的压缩采样磁共振图像重建方法还包括离散Curvelet变换的步骤：

Curvelet变换是在频域内实现的，采用频域窗函数U来表示Curvelet基函数

的傅里叶变换。定义了窗函数W(r)和角度窗函数V(t)，均满足可容许条件：

设U_j(ω)＝ψ_j(ω₁)V_j(ω)，其中ψ_j(ω₁)＝ψ(2^-jω₁)，ψ(ω₁)是一个带通滤波器，对于θ_j∈[-π/4,π/4)，有：

在构造上

则对于ω₁＞0，根据式(3-1)有：

当-1≤tanθ_l＜1，即

时，有：

又因为：

则根据式(3-6)满足：

通过对称和旋转才做可以将整个笛卡尔平面进行局部化分割，由

可以得到Curvelet基函数在时域的形式，则Curvelet变换表示：

本发明所提出的重建算法采用的是USFFT(Unequally Spaced fast FourierTransform)算法，实现步骤如下：

第一步，对于给定的笛卡尔坐标中的二维图像f[n₁,n₂]进行2DFFT变换，得到二维的频域表示为：

第二步，在频域，对于每一对(j,l)，对

重新采样或插值得到符合Curvelet频域局部化的形状

第三步，将

与拟合窗U_j,l相乘，得到

在拟合窗上的傅里叶变换

这就是在频域的Curvelet变换；

第四步，对每个乘积

做二维逆傅里叶变换，即得到离散Curvelet系数C^D(j,l,k)。

压缩感知理论：

假设原始信号X＝[x₁,x₂,...,x_N]^T∈R^N可由某一组原子

线性表示为：

将式(2-1)改写成矩阵形式，可以得到：

X＝ψΘ (2-2)

其中，ψ为N×N维的稀疏矩阵，Θ是在一个稀疏变换域中的投影系数，若Θ中仅有k个非零元素，则称X在ψ域上为k稀疏。设向量空间R^N中的离散实值信号X在变换Ψ下是稀疏的或是可压缩的，则可将X投影到随机测量矩阵Φ上，得到观测值向量Y，即：

Y＝ΦX＝ΦΨΘ＝A^CSΘ (2-3)

其中，CS信息算子A^CS＝ΦΨ，压缩感知观测向量的矩阵表示如图1所示的压缩感知观测向量的矩阵表示示意图。

测量值向量Y包含了稀疏信号X在变换Φ下的M个线性测量值，要精确重建原始信号，首先要保证观测矩阵满足RIP条件^[9]，在此前提下，重建X就是求解范数下的最优化问题，即求解l₀范数最小化，具体描述如下：

由于观测值Y的维数M远小于稀疏系数Θ的维数N，所以这是一个NP-Hard问题，为了解决这个问题，则必须将稀疏系数Θ中非零值的

种可能一一列举，使用穷举法求解，然而这种方式会耗费大量的时间。由于最小l₁范数和最小l₀范数能够在一定条件下互相转化，具有等价性。那么式(2-4)可转化为l₁最小范数下的最优化问题^[10]：

对于不确定的系统，Y是一个复矢量，考虑到由成像过程引起的噪音，所以式(2-5)可约束成以下形式：

MRI重建问题还包括了稀疏性和k空间数据一致性，由下面的拉格朗日约束条件^[11]给出：

其中，λ为正则化参数，用来平衡式前后两项所占的比重，利用优化算法求解式(2-7)，得到最优解Θ：

将

代入式(2-2)便可重建出原始信号。

离散Curvelet变换：

小波变换在某些应用中长期受到沿边缘信息表达能力不足的困扰，为克服这一不足，科学家们摆脱了对Ridgelet变换的依赖^[12]提出第二代离散Curvelet变换并且构造了Curvelet的紧框架。利用各向异性的基函数对曲线逼近。Curvelet变换是在频域内实现的，采用频域窗函数U来表示Curvelet基函数

的傅里叶变换。定义了窗函数W(r)和角度窗函数V(t)，均满足可容许条件：

设U_j(ω)＝ψ_j(ω₁)V_j(ω)，其中ψ_j(ω₁)＝ψ(2^-jω₁)，ψ(ω₁)是一个带通滤波器，对于θ_j∈[-π/4,π/4)，有：

在构造上

则对于ω₁＞0，根据式(3-1)有：

当-1≤tanθ_l＜1，即

时，有：

又因为：

则根据式(3-6)满足：

通过对称和旋转才做可以将整个笛卡尔平面进行局部化分割，由

可以得到Curvelet基函数在时域的形式，则Curvelet变换表示：

关于Curvelet离散变换的实现方法，E.Candes和Demanet^[13]给出了两种快速算法，因为本发明所提出的重建算法采用的是USFFT(Unequally Spaced fast FourierTransform)算法，所以本发明仅简要介绍USFFT的实现步骤，如下：

第一步，对于给定的笛卡尔坐标中的二维图像f[n₁,n₂]进行2DFFT变换，得到二维的频域表示为：

第二步，在频域，对于每一对(j,l)，对

重新采样或插值得到符合Curvelet频域局部化的形状

第三步，将

与拟合窗U_j,l相乘，得到

在拟合窗上的傅里叶变换

这就是在频域的Curvelet变换

第四步，对每个乘积

做二维逆傅里叶变换，即得到离散Curvelet系数C^D(j,l,k)

快速迭代收缩阈值算法：

快速迭代收缩阈值算法^[14]主要是利用前两次迭代的值和不断更新的参数t及不断收缩的正则化参数λ来获得新的迭代值，为重建MR二维图像，需求解l₁范数最小化：

由式(2-7)所示，令

g(Θ)＝||Θ||₁

arg min f(Θ)+λg(Θ) (4-1)

其中，正则化参数λ在式(4-1)中的作用是平衡数据f(Θ)与g(Θ)的比重。求解式(4-1)中f(Θ)的一个简单的方法就是引入梯度法^[15]，通过梯度法在每一次迭代中，不断修正Θ，具体如下(a>0)：

这等价于：

用同样的方法求解式(2-7)，迭代计算式如下：

省略掉常数项后，上式可化为：

由于l₁范数是可分离的，||Θ||₁为它的所有元素的绝对值之和，式(4-1)中g(Θ)的求解可以简化为求Θ的每一个最小化问题，这可以通过阈值收缩求得，则式(4-5)可转化为：

其中，shrink是软阈值算子，即：

shrink(x,β)＝sign(x)·max{|x|-β,0} (4-7)

为加快收敛速度，引入参数t和参数z，结合前2次的迭代值对其进行更新：

其中，Θ_k-1、Θ_k为Θ前2次的迭代值，此外，t的更新公式为：

由式(4-6)求得Θ_k+1，计算过程如下：

迭代的终止条件是通过Θ相邻2次迭代值的相对误差来设定的，终止函数定义为：

终止阈值设为ε，如果G_k(Θ_k,Θ_k-1)≤ε时，则停止迭代。

基于Curvelet变换的MR图像压缩感知重建算法：

采用了Curvelet变换结合快速迭代收缩阈值算法实现图像重建。基于Curvelet的快速迭代收缩阈值算法(Fast CurveletIterative Shrinkage-Thresholding Algorithm，FCISTA)简要步骤如下：

步骤(1)初始化：给定迭代停止阈值ε、图像分块观测矩阵Φ、图像观测矩阵Y，设定迭代计数器、Curvelet变换自适应阈值收缩算子及分别为k＝1、t＝0、λ＝λ₀、μ＝0.9；

步骤(2)Curvelet稀疏变换：

Θ_k+1＝C(X_k+1) (5-1)

其中，C为Curvelet稀疏变换；

步骤(3)对图像进行快速迭代收缩阈值处理，利用式(4-10)计算系数矩阵Θ_k；

步骤(4)更新参数t,利用式(4-9)计算t_k+1，根据式(4-8)计算z_k+1；

步骤(5)按式(4-11)计算迭代终止函数G_k(Θ_k,Θ_k-1)，如果G_k(Θ_k,Θ_k-1)≤ε，停止迭代，直接执行步骤(8)，如果G_k(Θ_k,Θ_k-1)＞ε，则继续往下执行步骤(6；

步骤(6)收缩正则化参数λ，λ＝0.9×λ；

步骤(7)更新迭代次数，令k＝k+1，返回到步骤(2)；

步骤(8)对迭代得到的系数矩阵进行Curvelet反变换，得到重建图像X：

X_k+1＝C^T(Θ_k+1) (5-2)

其中，C^T为Curvelet反变换；

步骤(9)显示重建图像X，算法结束。

实施例实验结果与分析：

本发明所提出的算法选择在自然图像(Lena图像，512×512像素)和医学图像(脑部MR图像，512×512像素)进行实验，如图2图3所示，其中，图3脑部MR图像来自江苏省人民医院磁共振成像检测部门的3T MRI扫描仪上完全采样所得的数据。

在固定值20％采样率的条件下，对图2、图3两幅图像进行重建仿真。N＝512×512，M为观测值Y的长度，采样率为M/N，实验所采用的是spiral欠采样填充轨迹^[17]对k空间^[18]中心区域采用多次采样，如图4所示。在相同的观测矩阵和稀疏基的基础上，对原始图像利用本发明提出的基于Curvelet变换的快速迭代收缩阈值算法进行仿真实验，并与基于Harr小波变换的快速迭代收缩阈值算法的仿真结果进行比较。FCISTA的参数设定为：δ＝0.005，ε＝10^-5，μ＝0.9，FISTA算法的参数设定与FCISTA算法一致。实验环境：一台PC机(CPU为3.6GHz，内存为4GB)，Matlab版本是R2015a。

1、重建图像细节

在选用spiral采样轨迹及0.2023采样率的基础上，将所提出的方法的结果与基于Harr小波变换的快速迭代收缩阈值算法进行比较，如图5所示：

图5中(a1)-(b1)表示在基于Curvelet变换的快速迭代收缩阈值算法下恢复的lena图像和脑部MR图像，(a2)-(b2)表示在基于Harr小波变换的快速迭代收缩阈值算法下恢复的lena图像和脑部MR图像，(a3)-(b3)与(a4)-(b4)分别表示在FCISTA算法与FISTA算法重建图像中提取的细节部分。对比相同部位细节图可以看出：FCISTA及FISTA算法均可恢复出绝大部分原始图像信息，但基于Harr小波变换的FISTA算法重建图像边缘都略显模糊且粗糙，本发明提出的基于Curvelet变换的快速迭代收缩阈值算法能较好地恢复图像纹理及细节信息，重建图像很接近原图像，所得重建结果视觉效果更佳。

2、差值图像

在相同的迭代终止条件ε＜10^-5的基础上，利用差值图像可以更好地比较FCISTA、FISTA算法重建的图像与原图的差别，因为FCISTA算法与FISTA算法的原始差值图像较不清晰，所以在原始差值图像基础上进行了亮度增强处理^[19]，正如图5所示：

图6中(a5)、(b5)为基于FCISTA算法重建出的图像与相应的原始图像相减处理后的差值图像，(a6)、(b6)为基于FISTA算法重建出的图像与相应的原始图像相减处理后的差值图像。由图6中差值图像所示，FCISTA算法的差值图像相比与FISTA算法的差值图像更为模糊，可以更好的说明本发明所提出的FCISTA算法可以更好地保存结构和恢复重建图像的更多信息，相比FISTA算法重建图像质量明显有所提升，体现了其在处理边缘信息以及细节上的优势。

3、质量评估参数

为了更进一步比较不同算法对图像的重建效果，图7与表1给出了FCISTA算法与FISTA算法重建的lena图像和脑部MR图像在均方误差(MSE)、峰值信噪比(PSNR)和结构化相似度(SSIM)之间^[20]的关系。

表1相同采样率下FCISTA算法与FISTA算法性能比较

参数	图像	FCISTA	FISTA
				MSE	lena图像	16.55121791	33.35619108
	脑部MR图像	4.790353206	29.09848574
				PSNR	lena图像	41.99709966	38.95363469
	脑部MR图像	47.38172386	39.54669533
				SSIM	lena图像	0.99980888	0.99899364
	脑部MR图像	0.999982976	0.99980483

由图7的6幅曲线图中可以看出，在不同采样率的情况下，FCISTA算法和FISTA算法均会随着采样率的提高而提升其恢复效果，然而在相同MSE、PSNR、SSIM值时，FCISTA算法可以在更低的采样率，恢复出FISTA算法相同的效果，以达到提高效率的目的。并且FCISTA算法重建图像的质量明显高于FISTA算法，对于图7中脑部MRI折线图来说，在采样率比较低时，FCISTA算法与FISTA算法重建图像的PSNR值相差不明显，当采样率M/N＞0.1时，FCISTA算法重建图像的峰值信噪比明显加大，而且随着采样率的增大，差距在加大。

由表1的数据可以看出，在相同采样率0.2023的条件下，本发明提出的FCISTA算法在处理lena图像和脑部MR图像过程中MSE、PSNR、SSIM方面均优于FISTA算法。由上述仿真结果可以看出，Curvelet变换具有更好的稀疏特性，基于Curvelet变换的快速迭代收缩阈值算法性能优于传统Harr小波变换快速迭代收缩阈值算法。

结论：本发明提出一种基于Curvelet变换的MR图像压缩感知重建方法，利用Curvelet稀疏表示的优越特性，对图像平滑部分和边缘部分进行表达。对应边缘的大系数能起到很好的能量集中效果，并且采用快速迭代收缩阈值进行截断，使得在重建过程中得到了边缘及细节部分重建效果的提升，图像整体重建质量也得到了提升。实验结果表明，本发明的FCISTA算法重建所得图像可满足人眼视觉的需要，相比于传统的MR图像重建算法，其重建图像质量及细节信息有了明显的提升，对于疾病的检查和判断有着积极的作用。本发明能利用更少的有效信息、更集中的数据，以较大概率重建图像，推进图像压缩感知重建研究的发展。

以上所述，仅是本发明的较佳实施例而已，并非对本发明的技术范围作出任何限制，故凡是依据本发明的技术实质对以上实施例所作的任何细微修改、等同变化与修饰，均仍属于本发明的技术方案的范围内。

Claims (1)

1.一种基于curvelet-fista的压缩采样磁共振图像重建方法，该方法采用了Curvelet变换结合快速迭代收缩阈值算法实现图像重建；基于Curvelet的快速迭代收缩阈值算法，其特征在于包括如下步骤：

步骤(1)初始化：给定迭代停止阈值ε、图像分块观测矩阵Φ、图像观测矩阵Y，设定迭代计数器、Curvelet变换自适应阈值收缩算子及分别为k＝1、t＝0、λ＝λ₀、μ＝0.9；

步骤(2)Curvelet稀疏变换：

Θ_k+1＝C(X_k+1) (5-1)

其中，C为Curvelet稀疏变换；

步骤(3)对图像进行快速迭代收缩阈值处理，利用快速迭代收缩阈值算法计算系数矩阵Θ_k；

步骤(4)更新参数t,利用快速迭代收缩阈值算法计算t_k+1，利用快速迭代收缩阈值算法计算z_k+1；

步骤(5)利用快速迭代收缩阈值算法计算迭代终止函数G_k(Θ_k,Θ_k-1)，如果G_k(Θ_k,Θ_k-1)≤ε，停止迭代，直接执行步骤(8)，如果G_k(Θ_k,Θ_k-1)＞ε，则继续往下执行步骤(6)；

步骤(6)收缩正则化参数λ，λ＝0.9×λ；

步骤(7)更新迭代次数，令k＝k+1，返回到步骤(2)；

步骤(8)对迭代得到的系数矩阵进行Curvelet反变换，得到重建图像X：

X_k+1＝C^T(Θ_k+1) (5-2)

其中，C^T为Curvelet反变换；

步骤(9)显示重建图像X，算法结束；所述快速迭代收缩阈值算法包括以下步骤：由拉格朗日约束条件：

所示，

令

arg min f(Θ)+λg(Θ) (4-1)

其中，正则化参数λ在式(4-1)中的作用是平衡数据f(Θ)与g(Θ)的比重；求解式(4-1)中f(Θ)的一个简单的方法就是引入梯度法，通过梯度法在每一次迭代中，不断修正Θ，具体如下(a>0)：

这等价于：

迭代计算式如下：

省略掉常数项后，上式可化为：

由于l₁范数是可分离的，||Θ||₁为它的所有元素的绝对值之和，式(4-1)中g(Θ)的求解可以简化为求Θ的每一个最小化问题，这可以通过阈值收缩求得，则式(4-5)可转化为：

其中，shrink是软阈值算子，即：

shrink(x,β)＝sign(x)·max{|x|-β,0} (4-7)

为加快收敛速度，引入参数t和参数z，结合前2次的迭代值对其进行更新：

其中，Θ_k-1、Θ_k为Θ前2次的迭代值，此外，t的更新公式为：

由式(4-6)求得Θ_k+1，计算过程如下：

迭代的终止条件是通过Θ相邻2次迭代值的相对误差来设定的，终止函数定义为：

终止阈值设为ε，如果G_k(Θ_k,Θ_k-1)≤ε时，则停止迭代；

所述步骤(3)对图像进行快速迭代收缩阈值处理，利用式(4-10)计算系数矩阵Θ_k；

所述步骤(4)更新参数t,利用式(4-9)计算t_k+1，根据式(4-8)计算z_k+1；

所述步骤(5)按式(4-11)计算迭代终止函数G_k(Θ_k,Θ_k-1)，如果G_k(Θ_k,Θ_k-1)≤ε，停止迭代，直接执行步骤(8)，如果G_k(Θ_k,Θ_k-1)＞ε；还包括离散Curvelet变换的步骤：

Curvelet变换是在频域内实现的，采用频域窗函数U来表示Curvelet基函数

的傅里叶变换；定义了窗函数W(r)和角度窗函数V(t)，均满足可容许条件：

设U_j(ω)＝ψ_j(ω₁)V_j(ω)，其中ψ_j(ω₁)＝ψ(2^-jω₁)，ψ(ω₁)是一个带通滤波器，对于θ_j∈[-π/4,π/4)，有：

在构造上

则对于ω₁＞0，根据式(3-1)有：

当-1≤tanθ_l＜1，即

时，有：

又因为：

则根据式(3-6)满足：

通过对称和旋转才做可以将整个笛卡尔平面进行局部化分割，由

可以得到Curvelet基函数在时域的形式，则Curvelet变换表示：

所提出的重建算法采用的是USFFT(Unequally Spaced fast Fourier Transform)算法，实现步骤如下：

第一步，对于给定的笛卡尔坐标中的二维图像f[n₁,n₂]进行2DFFT变换，得到二维的频域表示为：

第二步，在频域，对于每一对(j,l)，对

重新采样或插值得到符合Curvelet频域局部化的形状

第三步，将

与拟合窗U_j,l相乘，得到

在拟合窗上的傅里叶变换

这就是在频域的Curvelet变换；

第四步，对每个乘积

做二维逆傅里叶变换，即得到离散Curvelet系数C^D(j,l,k)。

Images