[发明专利]一种基于curvelet-fista的压缩采样磁共振图像重建方法

权利要求书

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(Θ)的一个简单的方法就是引入梯度法，通过梯度法在每一次迭代中，不断修正Θ，具体如下(a0)：

这等价于：

迭代计算式如下：

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

由于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)。