CN108828483A - 一种磁共振图像的收缩阈值迭代重建算法 - Google Patents

Abstract

本发明提供一种磁共振图像的收缩阈值迭代重建算法，包括提出带约束解空间的稀疏重建模型；提出阈值函数收缩规则；重建图像；加速图像重建速度直至修正收敛输出重建结果。本发明为CS_MRI重建提出了一个快速高效的重建算法，收缩阈值投影算法STPROX，利用l_p‑范数的非凸优化模型代替l₁‑范数凸优化模型，来减少采样点数目，同时采用阈值算法进行求解。迭代阈值算法不仅算法本身简单、高效，而且还可以有效的消除由于K空间欠采样所带来的伪影，为了使阈值函数适用l_p‑范数的非凸优化模型，避免收敛到局部最优解，对软阈值的阈值函数和收缩规则进行了修改以让其收敛到全局最优解。

Description

一种磁共振图像的收缩阈值迭代重建算法

技术领域

本发明涉及一种磁共振图像的收缩阈值迭代重建算法。

背景技术

MRI在医疗临床诊断中始终扮演着举足轻重的角色。尽管MRI是目前少有的对人体无伤害的准确的临床诊断方法，具有多方位、多参数、多模态等优点，但是它的扫描时间、成像速度却不尽人意，引入压缩感知理论，针对稀疏信号或可压缩信号，该方法是在获取信号的同时，就对数据进行适当的压缩，其采样频率可低于传统的奈奎斯特频率，其突出的优点是可减少采样数据，节省存储空间，但包含有足够的信息量。

在CS_MRI中有两种典型的稀疏模型：综合型和分解型。综合型模型是假设图像可以被某一字典中的元素稀疏表示，而分解型模型是指图像在某一变换域下是稀疏的。在正交系统下两种模型是等价的，如果在冗余系统下两种模型就有明显的差别。

大多数磁共振图像重建算法针对的是l₁-范数凸优化问题的求解，然而，实际应用中基于l₁-范数最小化的凸优化算法不能很好地去除数据间较大的冗余，所需的观测点数仍然较理论所需观测点数的最小值还大。

因此有必要为CS_MRI重建提出一个快速高效的重建算法，用以解决上述问题。

发明内容

根据上述提出的技术问题，而提供一种磁共振图像的收缩阈值迭代重建算法。本发明主要利用利用l_p-范数的非凸优化模型代替l₁-范数凸优化模型，再对稀疏重建模型进行阈值收缩和图像重建，从而让其收敛到全局最优解。

本发明采用的技术手段如下：

一种磁共振图像的收缩阈值迭代重建算法，其特征在于，包括如下步骤：

S1、提出带约束解空间的稀疏重建模型：临床磁共振成像从K空间域采集数据组成待重建的图像，通过傅里叶变换得到欠采样K空间数据y，y＝UFx，其中，x为待重建的图像，Fx为对图像x做傅里叶变换得到的完整K空间数据，U表示欠采算子是对完整K空间域进行欠采样；

用Ψ表示任意的紧框架变换，Ψ^H表示紧框架的共轭的转置；将N个像素的待重建的图像x表示成复数域的一个向量，待重建的图像x在紧框架变换Ψ的变换域系数α构成的空间为Range(Ψ)；利用l_p-范数的非凸优化模型代替l_1-范数凸优化模型；提出稀疏重建模型为：

其中，λ表示正则化参数，用来权衡和两项的重要性；

S2、阈值函数收缩规则：

其中，

S3、重建图像：提出的稀疏重建模型(1)中的变换域系数可以由公式(3)收缩阈值迭代运算求解：

α^k+1＝Γ_p(ΨΨ^Hα^k+γΨF^HU^T(y-UFΨ^Hα^k)) (3)

其中，α^k+1、α^k是第k+1和k次迭代的变换域系数；U^T表示U的转置；F^H表示傅里叶反变换；γ表示步长，取值范围γ＞0；Γ_p表示对给定向量z＝ΨΨ^Hα^k+γΨF^HU^T(y-UFΨ^Hα^k)的所有元素按公式(2)进行阈值收缩，阈值收缩的定义为：给定向量z，若向量z的第k个元素的绝对值|z_k|≤τ_p(λ)，则z_k＝0；若|z_k|＞τ_p(λ)，则z_k＝sgn(z)ST_p(|z_k|；λ)，其中sgn(z_k)是符号函数；

由于待重建的图像x是用紧框架Ψ的共轭转置与变换域系数α的乘积来表示，即：x＝Ψ^Hα，因此待重建的图像x按照公式(4)迭代运算获得：

x^k+1＝Ψ^HΓ_p(Ψ(x^k+γF^HU^T(y-UFx^k))) (4)

S4、加速图像重建的速度：引入一组因子{t_k}，初始化t₀＝1，通过表达式对第k次和k+1次迭代得到的图像进行修正，从而加快图像重建的速度，修正的方法：

在迭代的过程中若相邻两次的重建图像x^k和x^k+！的l₂范数缩小到指定的范围，则停止迭代，输出最后重建结果x^k+！。

较现有技术相比，本发明是为CS_MRI重建提出了一个快速高效的重建算法，收缩阈值投影算法STPROX，利用l_p-范数的非凸优化模型代替l₁-范数凸优化模型，来减少采样点数目，同时采用阈值算法进行求解。迭代阈值算法不仅算法本身简单、高效，而且还可以有效的消除由于K空间欠采样所带来的伪影，为了使阈值函数适用l_p-范数的非凸优化模型，避免收敛到局部最优解，对软阈值的阈值函数和收缩规则进行了修改以让其收敛到全局最优解。

本发明的有益效果是：参数少，仅引入了两个可调参数λ,γ，在实验调试中给出最优取值，求解模型是l_p-范数模型，可以减少采样点数，p的取值范围为0＜p＜1，可调节，既能保证算法的收敛速度又能提高图像的重建质量。

附图说明

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

图1是本发明实施例中的人脑图像。

图2是对k空间进行欠采样的采样模板。

图3是不同y值情况下，f(u)的变换图形。

图4是FISTA、ADMM、pFISTA、STPROX四种算法重建图像质量对比图。

图5是当p＝0.5时FISTA、ADMM、pFISTA、STPROX四种算法收敛曲线图。

图6是当p＝0.1时FISTA、ADMM、pFISTA、STPROX四种算法的收敛曲线图。

图7是当p＝0.8时FISTA、ADMM、pFISTA、STPROX四种算法的收敛曲线图。

具体实施方式

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

一种磁共振图像的收缩阈值迭代重建算法，包括如下步骤：

S1、提出带约束解空间的稀疏重建模型：

对于实施例中的图像x是256×256的人脑图像，如图1所示，根据图2的欠采样模板，采样率为30％的径向采样，使得采样模板中黑色对应位置的K空间数据不采，白色对应位置的K空间数据采集到。欠采到的K空间数据y总共有196601数据点。

临床磁共振成像从K空间域采集数据组成待重建的图像，通过傅里叶变换得到欠采样K空间数据y，y＝UFx，其中，x为待重建的图像，F表示对磁共振图像作傅里叶变换得到65536个傅里叶系数构成完整的K空间数据，U表示对完整的K空间数据用图1的采样模板进行欠采。选择4个分解集的冗余Daubechies小波为紧框架变换Ψ，Ψ^H表示紧框架的共轭的转置。提出带约束解空间的稀疏重建模型为：

其中，α是变换域系数，Range(α)是所有65532个像素图像x在紧框架Ψ下的投影系数所构成的空间，α∈Range(Ψ)就是解空间中的一个向量。λ表示正则化参数，用来权衡和的重要性。在实施人脑重建中取λ＝0.001。

S2、阈值函数收缩规则：

具体地，针对步骤02，τ_p(λ)的获得：

对于l_p-范数重建模型近似应用软迭代阈值求解的模目标函数可以定义为：

其中，z＝α-γ▽f(α)，观察到，f(u)的值在(0,+∞)的范围内是随着参数值的不同而不断变化。首先取定p＝0.5,λ＝1，在图3中显示f(u)随不同的取值z的变化情况。当p和λ取定值时，将存在一个特殊的阈值τ_p(λ)，使得当z＜τ_p(λ)时，x＝0是全局最小值；否则，x在一非零点处取得最小值。因此为了求解目标函数的最优值，新的阈值函数和收缩规则就显得尤为重要。

首先，对f(u)分别求一阶导数和二阶导数：

f(u)'＝λpu^p-1+u-z

f(u)”＝λp(p-1)u^p-2+1

令f(u)”＝0，得到拐点不难发现f(u)在区间上是凸函数，而在上是凹函数。为了保证f(u)在区间上有最小值，还需满足从图中我们不难发现，存在一个特殊的值y，使得与之对应的点处，有为了解出正确的阈值τ_p(λ)，可以利用两个临界条件代入得到如下的数学公式：

由公式可得将τ_p(λ)的值代入公式可以得到区间上唯一的解和阈值τ_p(λ)：

由此可以得到新的阈值规则：

收缩阈值算法：

下面给出改进的阈值算法l_p-ST，Γ_p(z；λ)＝l_p-ST(z,λ,p,m)的具体步骤：

(1)初始化，输入z,λ,p,m，m是迭代次数；

(2)比较输入值y与阈值τ_p(λ)的大小，如果|z|≤τ_p(λ)，其中τ_p(λ)为公式的值，则令Γ_p(z；λ)＝0。否则进行步骤(3)

(3)求解非零点出的最小值，令k＝0，u^k＝|z|。进行迭代，k＝0,1,···,m，对y进行阈值收缩得到更新的值u^k+1＝|z|-λp(u^k)^p-1；

(4)达到迭代次数，出处阈值收缩函数Γ_p(z；λ)的值。

S3、重建图像：提出的稀疏重建模型(1)中的变换域系数可以由公式(3)收缩阈值迭代运算求解：

α^k+1＝Γ_p(ΨΨ^Hα^k+γΨF^HU^T(y-UFΨ^Hα^k)) (3)

其中，α^k+1、α^k是第k+1和k次迭代的变换域系数；U^T表示U的转置；F^H表示傅里叶反变换；γ表示步长，取值范围γ＞0；Γ_p表示对给定向量z＝ΨΨ^Hα^k+γΨF^HU^T(y-UFΨ^Hα^k)的所有元素按公式(2)进行阈值收缩，阈值收缩的定义为：给定向量z，若向量z的第k个元素的绝对值|z_k|≤τ_p(λ)，则z_k＝0；若|z_k|＞τ_p(λγ)，则z_k＝sgn(z)ST_p(|z_k|；λ)，其中sgn(z_k)是符号函数；

由于待重建的图像x是用紧框架Ψ的共轭转置与变换域系数α的乘积来表示，即：x＝Ψ^Hα，因此待重建的人脑图像x按照公式(4)迭代运算获得：

x^k+1＝Ψ^HΓ_p(Ψ(x^k+γF^HU^T(y-UFx^k))) (4)

其中，γ＝1。

S4、加速图像重建的速度：引入一组因子{t_k}，初始化t₀＝1，通过表达式对第k次和k+1次迭代得到的图像进行修正，从而加快图像重建的速度，修正的方法：

在迭代的过程中若相邻两次的重建图像x^k和x^k+！的l₂范数缩小到指定的范围，则停止迭代，输出最后重建结果x^k+！。

图4显示了新提出的算法与FISTA、ADMM、pFISTA的图像重建质量的对比，其中p＝0.5。从中可以看出STPROX的重建图像质量最高，其次是pFISTA，而FISTA的重建图像伪影比较明显。

图5是四种算法收敛曲线的对比图，其中p＝0.5。如图所示，无论是在收敛速度上还是重建误差上，新提出的算法都有明显的优势。

图6与图7是为了验证新提出算法的另一个重要优势，对于p的选择，当p取接近于0的值时，STPROX的收敛速度明显加快，如图6所示，其中p＝0.1；当p的取值接近1的时候，STPROX的重建误差会降低更多，图像的重建质量更优，如图7所示，其中p＝0.8。

最后应说明的是：以上各实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述各实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分或者全部技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明各实施例技术方案的范围。

Claims (1)

1.一种磁共振图像的收缩阈值迭代重建算法，其特征在于，包括如下步骤：

S1、提出带约束解空间的稀疏重建模型：临床磁共振成像从K空间域采集数据组成待重建的图像，通过傅里叶变换得到欠采样K空间数据y，y＝UFx，其中，x为待重建的图像，Fx为对图像x做傅里叶变换得到的完整K空间数据，U表示欠采算子是对完整K空间域进行欠采样；

用Ψ表示任意的紧框架变换，Ψ^H表示紧框架的共轭的转置；将N个像素的待重建的图像x表示成复数域的一个向量，待重建的图像x在紧框架变换Ψ的变换域系数α构成的空间为Range(Ψ)；利用l_p-范数的非凸优化模型代替l₁-范数凸优化模型；提出稀疏重建模型为：

其中，λ表示正则化参数，用来权衡和两项的重要性，p的取值为：0＜p＜1；

S2、阈值函数收缩规则：

其中，0＜p＜1；

S3、重建图像：提出的稀疏重建模型(1)中的变换域系数可以由公式(3)收缩阈值迭代运算求解：

α^k+1＝Γ_p(ΨΨ^Hα^k+γΨF^HU^T(y-UFΨ^Hα^k)) (3)

其中，α^k+1、α^k是第k+1和k次迭代的变换域系数；U^T表示U的转置；F^H表示傅里叶反变换；γ表示步长，取值范围γ＞0；Γ_p表示对给定向量z＝ΨΨ^Hα^k+γΨF^HU^T(y-UFΨ^Hα^k)的所有元素按公式(2)进行阈值收缩，阈值收缩的定义为：给定向量z，若向量z的第k个元素的绝对值|z_k|≤τ_p(λ)，则z_k＝0；若|z_k|＞τ_p(λ)，则z_k＝sgn(z)ST_p(|z_k|；λ)，其中sgn(z_k)是符号函数；

由于待重建的图像x是用紧框架Ψ的共轭转置与变换域系数α的乘积来表示，即：x＝Ψ^Hα，因此待重建的图像x按照公式(4)迭代运算获得：

x^k+1＝Ψ^HΓ_p(Ψ(x^k+γF^HU^T(y-UFx^k))) (4)

S4、加速图像重建的速度：引入一组因子{t_k}，初始化t₀＝1，通过表达式对第k次和k+1次迭代得到的图像进行修正，从而加快图像重建的速度，修正的方法：

在迭代的过程中若相邻两次的重建图像x^k和x^k+！的l₂范数缩小到指定的范围，则停止迭代，输出最后重建结果x^k+！。