CN103049923A - 磁共振快速成像的方法 - Google Patents

Abstract

本发明属于图像处理技术领域，尤其涉及一种磁共振快速成像的方法。本发明的共振快速成像的方法包括：步骤A：在图像梯度域的水平梯度图像和垂直梯度图像上进行字典学习，建立图像模型；步骤B：利用重建算法交替更新图像块的稀疏表示，恢复水平梯度和垂直梯度，然后在这两个方向梯度重建图像。本发明实施例磁共振快速成像的方法通过自适应学习字典的引入可以克服固定的有限差分变换导致目标图像的块状效应，可以处理结构更复杂的图像，从而产生更精确的重建；另外本发明实施例对梯度图像进行处理，其比原图像更稀疏，因而字典学习可以变得更精确和鲁棒，从而更稀疏地表示图像，且具有更好的保真度，且可以恢复出更多的细节。

Description

磁共振快速成像的方法

技术领域

本发明属于图像处理技术领域，尤其涉及一种磁共振快速成像的方法。

背景技术

为了缩短磁共振图像采集时间，压缩感知理论被成功应用到磁共振成像中。压缩感知理论利用信号在某个基的稀疏性，实现了在非相干采样矩阵下，只需少量采样(远少于奈奎斯特采样理论所需的采样)即可高质量重建原始信号。而且，若信号在某个基越稀疏，那么所需要的采样量则越少。因此，在压缩感知理论中，一个重要的问题就是稀疏基的选取。在过去几年中，从欠采样K空间信号来重建磁共振图像一般使用如下全变分模型(TV模型)：

$\underset{u}{\min }\{\mathrm{\μ}{\left|\right|u\left|\right|}_{\mathrm{TV}}+\frac{1}{2}{\left|\right|F_{p}u-f\left|\right|}_2^2\}---\left(1\right)$

其中， ${\left|\right|u\left|\right|}_{\mathrm{TV}}={\left|\right|\▿u\left|\right|}_1=\underset{i,j}{\mathrm{\Σ}}\left|{\left({\▿}_{x}u\right)}_{i,j}\right|+\left|{\left({\▿}_{y}u\right)}_{i,j}\right|,$ 称为各向异性离散总变分正则化方程，

定义为水平方向和垂直方向的差分算子；μ＞0为罚参数；第二项

作为保真项，用l₂-norm控制误差；我们定义

为需要重建的图像，

表示傅里叶欠采样信号，两者的关系为F_pu＝f，其中表示傅里叶欠采样矩阵。

TV模型有着很好的保存图像边缘的能力。但是在较大的欠采样下可能会产生块状效应，因此把该模型应用于磁共振快速成像上时，图像重建质量会受到一定的影响。为了提高图像重建质量，可以在TV模型基础上添加其它稀疏约束项(例如小波变换)，例如近年来由Yang等人提出的RecPF快速重建方法模型如下：

$\underset{u}{\min }\{{\mathrm{\μ}}_1{\left|\right|u\left|\right|}_{\mathrm{TV}}+{\mathrm{\μ}}_2{\left|\right|\mathrm{\ψu}\left|\right|}_1+\frac{1}{2}{\left|\right|F_{p}u-f\left|\right|}_2^2\}---\left(2\right)$

其中，ψ表示小波变换，μ₁，μ₂＞0，用于权衡前两个正则化项和保真项。

除了引入基于已知固定的稀疏变换的正则化来提高稀疏性，近年来，基于图像块字典学习的图像稀疏表示方法得到了越来越多的研究。其典型的方法描述如下：

对一幅

的图像u，用一个图像块提取算子R对图像进行分块，得到包含L个图像块的集合R(u)＝[R₁u，R₂u，…，R_Lu]，

定义为大小为

的图像块的向量形式。Elad等人提出，对于所有的图像块R_lu，都可以在字典D^*上稀疏的表示，其模型如下：

${\mathrm{\α}}_l^{*}=\underset{{\mathrm{\α}}_l}{\arg \min }{\left|\right|D^{*}{\mathrm{\α}}_l-R_{l}u\left|\right|}_2^{2}s.t.{\left|\right|{\mathrm{\α}}_l\left|\right|}_0\≤T_0,l=\mathrm{1,2},\·\·\·,L---\left(3\right)$

其中，α_l表示第l个图像块在字典D^*上的表示系数，T₀控制表示系数的稀疏度，||·||₀定义为向量中非零元素的个数。

Ravishanker等人把字典学习模型应用到K空间欠采样的磁共振图像重建上，提出了DLMRI模型：

$\underset{u,D,\mathrm{\Γ}}{\min }\{\underset{l}{\mathrm{\Σ}}{\left|\right|D{\mathrm{\α}}_l-R_{l}u\left|\right|}_2^2+v{\left|\right|F_{p}u-f\left|\right|}_2^2\}$

$s.t.{\left|\right|{\mathrm{\α}}_l\left|\right|}_0\≤T_0,\∀l---\left(4\right)$

其中，Γ＝[α₁，α₂，…，α_L]定义为所有图像块对应的稀疏系数矩阵。前一项保证图像块在自适应的学习字典上稀疏表示，后一项为在K空间上的信号保真项，而正则化参数v平衡这两项间的权重。求解该模型一般用两步迭代交替更新的方法：第一步，应用K-SVD算法训练稀疏表示字典来消除混淆和噪声；第二步，通过学习字典和稀疏系数来重建图像。

现有技术的缺陷在于：在TV模型和RecPF模型中，TV正则化及小波稀疏约束等作为非自适应的固定变换，并不能理想地稀疏表示所有的图像；而Elad等人和Ravishanker等人提出的字典学习方法都是在图像域上训练字典，不能使图像更稀疏的表示，导致需要采样的信号多，另外上述的方法在图像域上学习反应细节较差。

发明内容

本发明提供了一种磁共振快速成像的方法，旨在解决现有的磁共振成像方法不能理想地稀疏表示所有的图像，导致需要采样信号多且成像反应细节差的技术问题。

本发明提供的技术方案为：一种磁共振快速成像的方法，包括：

步骤A：在图像梯度域的水平梯度图像和垂直梯度图像上进行字典学习，建立图像模型；

步骤B：利用重建算法交替更新图像块的稀疏表示，恢复水平梯度和垂直梯度，然后在这两个方向梯度重建图像。

本发明的技术方案还包括：所述步骤A还包括：在图像域上进行字典学习，即在图像域和梯度域的水平梯度图像和垂直梯度图像上同时进行字典学习，建立图像模型。

本发明的技术方案还包括：在所述步骤B中，恢复图像及其水平、垂直梯度图像，然后根据图像以及水平、垂直方向梯度重建图像。

本发明的技术方案还包括：在所述步骤A中，建立的图像模型为：

$\underset{u,D^{\left(i\right)},{\mathrm{\Γ}}^{\left(i\right)}}{\min }\{\underset{i=1}{\overset{2}{\mathrm{\Σ}}}\underset{l}{\mathrm{\Σ}}{\left|\right|D^{\left(i\right)}{\mathrm{\α}}_l^{\left(i\right)}-R_l\left({\▿}^{\left(i\right)}u\right)\left|\right|}_2^2+\frac{v_1}{2}{\left|\right|F_{p}u-f\left|\right|}_2^2\}$

$s.t.{\left|\right|{\mathrm{\α}}_l^{\left(i\right)}\left|\right|}_0\≤T_0,\∀l,i,$ 其中，第一项为梯度图像在字典上稀疏表示，第二项保证重建结果与K空间采样信号相匹配；权重v₁＝(λ/σ)，σ为测量噪声的标准差，λ为正常数。

本发明的技术方案还包括：所述步骤B的重建算法包括：引入辅助变量w⁽ⁱ⁾，i＝1，2，应用布雷格曼技术，定义 $\▿=\left[\begin{array}{c}{\▿}^{\left(1\right)}\\ {\▿}^{\left(2\right)}\end{array}\right],$ $b=\left[\begin{array}{c}{b}^{\left(1\right)}\\ b^{\left(2\right)}\end{array}\right],$ $w=\left[\begin{array}{c}{w}^{\left(1\right)}\\ w^{\left(2\right)}\end{array}\right].$

本发明的技术方案还包括：所述步骤B的重建算法还包括：更新梯度图像变量w⁽ⁱ⁾，i＝1，2；更新梯度图像块的稀疏表示；交替更新字典D⁽ⁱ⁾和系数矩阵

在稀疏编码阶段，固定字典D⁽ⁱ⁾，通过贪婪算法正交匹配追踪来更新

在字典更新阶段，固定系数

通过奇异值分解来逐列更新字典的每一列，最小化近似误差。

本发明的技术方案还包括：所述建立的图像模型为：

$\begin{array}{c}\underset{u,D,\mathrm{\ΓB},X}{\min }\{\underset{l}{\mathrm{\Σ}}{\left|\right|D{\mathrm{\α}}_l-R_{l}u\left|\right|}_2^2+\mathrm{\τ}\underset{l}{\mathrm{\Σ}}{\left|\right|B{x}_l-R_l(\▿u)\left|\right|}_2^2+{\mathrm{\μ}}_1{\left|\right|F_{p}u-f\left|\right|}^2\}\\ s.t.{\left|\right|{\mathrm{\α}}_l\left|\right|}_0\≤T_0,{\left|\right|x_l\left|\right|}_0\≤T_0,\∀l\end{array},$ 其中，D，B＝[B₁，B₂]分别为原图像块和梯度图像块对应的字典， ${\mathrm{\α}}_l,x_l=\left[\begin{array}{c}{x}_l^1\\ x_l^2\end{array}\right]$ 为相应的系数；τ，μ₁为权重参数。

本发明的技术方案还包括：所述重建算法包括：引入辅助变量v，应用布雷格曼技术，将图像模型变更为：

$\{u^{k+1},v^{k+1},D^{k+1},{\mathrm{\α}}_l^{k+1},B^{k+1},x_l^{k+1}\}=\arg \underset{u,v,D,\mathrm{\Γ},B,X}{\min }\left\{\begin{array}{c}\underset{l}{\mathrm{\Σ}}{\left|\right|D{\mathrm{\α}}_l-R_{l}u\left|\right|}_2^2+\mathrm{\τ}\underset{l}{\mathrm{\Σ}}{\left|\right|B{x}_l-R_l\left(v\right)\left|\right|}_2^2+{\mathrm{\μ}}_1{\left|\right|F_{p}u-f\left|\right|}^2+{\mathrm{\μ}}_2{\left|\right|b_1^k+\▿u-v\left|\right|}_2^2\\ s.t.{\left|\right|{\mathrm{\α}}_l\left|\right|}_0\≤T_0,{\left|\right|x_l\left|\right|}_0\≤T_0,\∀l\end{array}.\right.$

本发明的技术方案还包括：所述重建算法还包括：更新u，固定其余的变量；更新变量v；更新字典及系数：B和x_l，D和α_l，l＝1，2，…，L，交替更新字典D⁽ⁱ⁾和系数矩阵

在稀疏编码阶段，固定字典D⁽ⁱ⁾，通过贪婪算法正交匹配追踪来更新

在字典更新阶段，则固定系数

通过奇异值分解来逐列更新字典的每一列，最小化近似误差。

本发明的技术方案具有如下优点或有益效果：本发明实施例磁共振快速成像的方法通过自适应学习字典的引入可以克服固定的有限差分变换导致目标图像的块状效应，可以处理结构更复杂的图像，从而产生更精确的重建；另外本发明实施例对梯度图像进行处理，其比原图像更稀疏，因而字典学习可以变得更精确和鲁棒，从而更稀疏地表示图像，且具有更好的保真度，且可以恢复出更多的细节。

附图说明

附图1是本发明第一实施例的磁共振快速成像的方法的流程图；

附图2是本发明第二实施例的磁共振快速成像的方法的流程图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。

请参阅图1，为本发明第一实施例的磁共振快速成像的方法的流程图。本发明第一实施例的磁共振快速成像的方法包括：

步骤100：在图像梯度域的水平梯度图像和垂直梯度图像上进行字典学习，建立图像模型；

其中，自适应学习字典的引入可以克服有限差分变换导致重建图像的块状效应，同时也可以更好地表示结构复杂的图像，从而产生更精确的重构，在稀疏的变换域上学习字典，其稀疏表示效果越好，重建的效果也因此更好。对于图像的一阶梯度，其水平梯度图像和垂直梯度图像比原图像更稀疏。因此，为了同时利用字典学习的优点和磁共振图像在梯度域上的稀疏性质，本发明第一实施例的磁共振快速成像在图像梯度域(包含两个方向，水平方向和垂直方向)上进行字典学习的方法，建立如下新模型(称为TVDL模型)：

$\underset{u,D^{\left(i\right)},{\mathrm{\Γ}}^{\left(i\right)}}{\min }\{\underset{i=1}{\overset{2}{\mathrm{\Σ}}}\underset{l}{\mathrm{\Σ}}{\left|\right|D^{\left(i\right)}{\mathrm{\α}}_l^{\left(i\right)}-R_l\left({\▿}^{\left(i\right)}u\right)\left|\right|}_2^2+\frac{v_1}{2}{\left|\right|F_{p}u-f\left|\right|}_2^2\}---\left(5\right)$

$s.t.{\left|\right|{\mathrm{\α}}_l^{\left(i\right)}\left|\right|}_0\≤T_0,\∀l,i$

其中，令

第一项使梯度图像在字典上稀疏表示，同时，第二项保证重建结果与K空间采样信号相匹配；权重v₁＝(λ/σ)，σ为测量噪声的标准差，λ为正常数。

步骤110：利用重建算法交替更新图像块的稀疏表示，恢复水平梯度和垂直梯度，然后在这两个方向梯度重建图像。

对于本发明第一实施例提出的模型(5)，可以应用Bregman技术来求解，重建算法称为GradDLRec，该算法交替更新图像块的稀疏表示，恢复水平梯度和垂直梯度，然后由这两个方向梯度重建图像。算法的具体描述如下：

引入辅助变量w⁽ⁱ⁾，i＝1，2，方程(5)可以改写如下：

$\underset{u,w,D^{\left(i\right)},{\mathrm{\Γ}}^{\left(i\right)}}{\min }\{\underset{i=1}{\overset{2}{\mathrm{\Σ}}}\underset{l}{\mathrm{\Σ}}{\left|\right|D^{\left(i\right)}{\mathrm{\α}}_l^{\left(i\right)}-R_l\left(w^{\left(i\right)}\right)\left|\right|}_2^2+{v_1\left|\right|F_{p}u-f\left|\right|}_2^2\}---\left(6\right)$

$s.t.{\left|\right|{\mathrm{\α}}_l^{\left(i\right)}\left|\right|}_0\≤T_0,\∀l,i;w^{\left(i\right)}={\▿}^{\left(i\right)}u,\∀i;$

通过应用Bregman方法，定义 $\▿=\left[\begin{array}{c}{\▿}^{\left(1\right)}\\ {\▿}^{\left(2\right)}\end{array}\right],$ $b=\left[\begin{array}{c}{b}^{\left(1\right)}\\ b^{\left(2\right)}\end{array}\right],$ $w=\left[\begin{array}{c}{w}^{\left(1\right)}\\ w^{\left(2\right)}\end{array}\right],$ 可以得到如下子问题：

$\{u^{k+1},w^{k+1},{\left(D^{\left(i\right)}\right)}^{k+1},{\left({\mathrm{\α}}_l^{\left(i\right)}\right)}^{k+1}\}=\underset{u,w,D,\mathrm{\Γ}}{\arg \min }\left\{\begin{array}{c}\underset{i=1}{\overset{2}{\mathrm{\Σ}}}\underset{l}{\mathrm{\Σ}}{\left|\right|D^{\left(i\right)}{\mathrm{\α}}_l^{\left(i\right)}-R_l\left(w^{\left(i\right)}\right)\left|\right|}_2^2+v_1{\left|\right|F_{p}u-f\left|\right|}_2^2+v_2{\left|\right|b^k+\▿u-w\left|\right|}_2^2\\ s.t.{\left|\right|{\mathrm{\α}}_l^{\left(i\right)}\left|\right|}_0\≤T_0,\∀l,i\end{array}\right\}---\left(7\right)$

$b^{k+1}=b^k+\▿u^{k+1}-w^{k+1},---\left(8\right)$

其中v₂定义为正罚参数。应用交替方向法(ADM，Alternating DirectionMethod)来求解问题(7)。问题(7)求解过程如下：

步骤71：更新u，在第k次迭代，固定w，D⁽ⁱ⁾，假设他们的值分别为

去除常量，更新u的目标函数为：

$u^{k+1}=\underset{u}{\arg \min }\{v_1{\left|\right|F_{p}u-f\left|\right|}_2^2+v_2{\left|\right|b^k+\▿u-w^k\left|\right|}_2^2\}---\left(9\right)$

由于(9)是一个简单的最小二乘问题，可以直接求出变量u的解析解：

$u^{k+1}=F^{-1}\left(\frac{F[v_1{F}_p^{T}f+v_2{\▿}^T(w^k-b^k)]}{v_1{\mathrm{FF}}_p^T{F}_p{F}^T+v_{2}F{\▿}^T{F}^{T}F\▿F^T}\right)---\left(10\right)$

其中，

定义为正规化的全采样傅里叶矩阵，即F^TF＝1_N。矩阵

为只包含0和1的对角矩阵，1在对角线上，其位置对应于K空间采样到的位置。

步骤72：更新梯度图像变量w⁽ⁱ⁾，i＝1，2(分别对应水平、垂直方向)；

由于w⁽¹⁾和w⁽²⁾是分离的，因此可以分别求解：

${\left(w^{\left(i\right)}\right)}^{k+1}=\underset{w^{\left(i\right)}}{\arg \min }\{\underset{l}{\mathrm{\Σ}}{\left|\right|{\left(D^{\left(i\right)}\right)}^k{\left({\mathrm{\α}}_l^{\left(i\right)}\right)}^k-R_l\left(w^{\left(i\right)}\right)\left|\right|}_2^2+v_2{\left|\right|{\left(b^{\left(i\right)}\right)}^k+{\left({\▿}^{\left(i\right)}u\right)}^{k+1}-w^{\left(i\right)}\left|\right|}_2^2\}---\left(11\right)$

该最小二乘问题的解析解如下：

${\left(w^{\left(i\right)}\right)}^{k+1}=\frac{v_2[{\left(b^{\left(i\right)}\right)}^k+{\left({\▿}^{\left(i\right)}u\right)}^{k+1}]+\underset{l}{\mathrm{\Σ}}{R}_l^T{\left(D^{\left(i\right)}\right)}^k{\left({\mathrm{\α}}_l^{\left(i\right)}\right)}^k/\mathrm{\β}}{v_2+1}---\left(12\right)$

步骤73：更新梯度图像块的稀疏表示(D⁽ⁱ⁾和

由于问题(7)中关于梯度图像在水平和垂直方向上的字典和系数是分离的，因此可以分别求解两个方向上对应的字典及系数：

$\{{\left(D^{\left(i\right)}\right)}^{k+1},{\left({\mathrm{\α}}_l^{\left(i\right)}\right)}^{k+1}\}=\underset{D^{\left(i\right)},{\mathrm{\Γ}}^{\left(i\right)}}{\arg \min }\left\{\begin{array}{c}\underset{l}{\mathrm{\Σ}}{\left|\right|D^{\left(i\right)}{\mathrm{\α}}_l^{\left(i\right)}-R_l{\left(w^{\left(i\right)}\right)}^{k+1}\left|\right|}_2^2\\ s.t.{\left|\right|{\mathrm{\α}}_l^{\left(i\right)}\left|\right|}_0\≤T_0,\∀l\end{array}\right\},i=\mathrm{1,2}---\left(13\right)$

求解(13)的方法与模型(3)使用的K-SVD方法一样，即交替更新字典D⁽ⁱ⁾和系数矩阵在稀疏编码阶段，固定字典D⁽ⁱ⁾，通过贪婪算法正交匹配追踪(OMP)来更新

在字典更新阶段，则固定系数

通过奇异值分解(SVD)来逐列更新字典的每一列，从而最小化近似误差。

以上即为求解问题(7)的过程。综上所述，本发明实施例提出的完整的GradDLRec算法可归纳如下：

算法1：GradDLRec算法，包括：

1：初始化： ${\left({\mathrm{\Γ}}^{\left(i\right)}\right)}^0=0,{\left(D^{\left(i\right)}\right)}^0,{\left(b^{\left(i\right)}\right)}^0=0,i=\mathrm{1,2};u^0=F_p^{T}f;$

2：For k＝1，2，…直到满足终止准则：

3： ${\left(w^{\left(i\right)}\right)}^{k+1}=\frac{v_2[{\left(b^{\left(i\right)}\right)}^k+{\left({\▿}^{\left(i\right)}u\right)}^{k+1}]+\underset{l}{\mathrm{\Σ}}{R}_l^T{\left(D^{\left(i\right)}\right)}^k{\left({\mathrm{\α}}_l^{\left(i\right)}\right)}^k/\mathrm{\β}}{v_2+1},i=\mathrm{1,2}$

4：更新 $\{{\left(D^{\left(i\right)}\right)}^{k+1},{\left({\mathrm{\α}}_l^{\left(i\right)}\right)}^{k+1}\},i=\mathrm{1,2}$

5： $u^{k+1}=F^{-1}\left(\frac{F[v_1{F}_p^{T}f+v_2{\▿}^T(w^k-b^k)}{v_{1}F{F}_p^T{F}_p{F}^T+v_{2}F{\▿}^T{F}^{T}F\▿F^T}\right)$

6： ${\left(b^{\left(i\right)}\right)}^{k+1}={\left(b^{\left(i\right)}\right)}^k+{\left({\▿}^{\left(i\right)}u\right)}^{k+1}-{\left(w^{\left(i\right)}\right)}^{k+1},i=\mathrm{1,2}$

7：End

8：输出u^k+1

请参阅图2，为本发明第二实施例的磁共振快速成像的方法的流程图。本发明第二实施例的磁共振快速成像的方法包括：

步骤200：在图像域和梯度域的水平梯度图像和垂直梯度图像上同时进行字典学习，建立图像模型；

在步骤200中，在图像及其水平和垂直梯度图像上同时进行字典学习，建立如下图像模型：

$\underset{u,D,\mathrm{\ΓB},X}{\min }\left\{\underset{l}{\mathrm{\Σ}}\right||{\left[\begin{array}{ccc}D& \sqrt{\mathrm{\τ}}{B}_1& \sqrt{\mathrm{\τ}}{B}_2\end{array}\right]\left[\begin{array}{c}{\mathrm{\α}}_l\\ x_l^1\\ x_l^2\end{array}\right]-R_l\left[\begin{array}{c}u\\ {\▿}_{1}u\\ {\▿}_{2}u\end{array}\right]\left|\right|}_2^2+{\mathrm{\μ}}_1{\left|\right|F_{p}u-f\left|\right|}^2\}---\left(14\right)$

$s.t.{\left|\right|{\mathrm{\α}}_l\left|\right|}_0\≤T_0,{\left|\right|x_l\left|\right|}_0\≤T_0,\∀l$

由于 ${\left|\right|\left[\begin{array}{ccc}D& \sqrt{\mathrm{\τ}}{B}_1& \sqrt{\mathrm{\τ}}{B}_2\end{array}\right]\left[\begin{array}{c}{\mathrm{\α}}_l\\ x_l^1\\ x_l^2\end{array}\right]-R_l\left[\begin{array}{c}u\\ {\▿}_{1}u\\ {\▿}_{2}u\end{array}\right]\left|\right|}_2^2\≤\underset{l}{\mathrm{\Σ}}{\left|\right|D{\mathrm{\α}}_l-R_{l}u\left|\right|}_2^2+\mathrm{\τ}\underset{l}{\mathrm{\Σ}}{\left|\right|B{x}_l-R_l(\▿u)\left|\right|}_2^2,$ 因此提出如下新模型(称为HTVDL模型)：

$\underset{u,D,\mathrm{\ΓB},X}{\min }\{\underset{l}{\mathrm{\Σ}}{\left|\right|D{\mathrm{\α}}_l-R_{l}u\left|\right|}_2^2+\mathrm{\τ}\underset{l}{\mathrm{\Σ}}{\left|\right|{\mathrm{Bx}}_l-R_l(\▿u)\left|\right|}_2^2+{\mathrm{\μ}}_1{\left|\right|F_{p}u-f\left|\right|}^2\}---\left(15\right)$

$s.t.{\left|\right|{\mathrm{\α}}_l\left|\right|}_0\≤T_0,{\left|\right|x_l\left|\right|}_0\≤T_0,\∀l$

其中，D，B＝[B₁，B₂]分别为原图像块和梯度图像块对应的字典， ${\mathrm{\α}}_l,x_l=\left[\begin{array}{c}{x}_l^1\\ x_l^2\end{array}\right]$ 为相应的系数；τ，μ₁为权重参数。

步骤210：利用重建算法交替更新图像块的稀疏表示，恢复图像及其水平、垂直梯度图像，然后根据图像以及水平、垂直方向梯度重建图像。

引入辅助变量v，问题(15)可以改写成如下形式：

$\underset{u,v,D,\mathrm{\Γ},B,X}{\min }\{\underset{l}{\mathrm{\Σ}}{\left|\right|D{\mathrm{\α}}_l-R_{l}u\left|\right|}_2^2+\mathrm{\τ}\underset{l}{\mathrm{\Σ}}{\left|\right|B{x}_l-R_l\left(v\right)\left|\right|}_2^2+{\mathrm{\μ}}_1{\left|\right|F_{p}u-f\left|\right|}^2\}---\left(16\right)$

$s.t.{\left|\right|{\mathrm{\α}}_l\left|\right|}_0\≤T_0,{\left|\right|x_l\left|\right|}_0\≤T_0,\∀l;v=\▿u;$

应用Bregman技术，问题(16)可以通过如下方法求解：

$\{u^{k+1},v^{k+1},D^{k+1},{\mathrm{\α}}_l^{k+1},B^{k+1},x_l^{k+1}\}=\arg \underset{u,v,D,\mathrm{\Γ},B,X}{\min }\left\{\begin{array}{c}\underset{l}{\mathrm{\Σ}}{\left|\right|D{\mathrm{\α}}_l-R_{l}u\left|\right|}_2^2+\mathrm{\τ}\underset{l}{\mathrm{\Σ}}{\left|\right|B{x}_l-R_l\left(v\right)\left|\right|}_2^2+{\mathrm{\μ}}_1{\left|\right|F_{p}u-f\left|\right|}^2+{\mathrm{\μ}}_2{\left|\right|b_1^k+\▿u-v\left|\right|}_2^2\\ s.t.{\left|\right|{\mathrm{\α}}_l\left|\right|}_0\≤T_0,{\left|\right|x_l\left|\right|}_0\≤T_0,\∀l\end{array}\right\}---\left(17\right)$

$b_1^{k+1}=b_1^k+\▿u^{k+1}-v^{k+1}---\left(18\right)$

其中，μ₂为正罚参数。

问题(17)的求解过程如下：

步骤171：更新u，固定其余的变量，问题(17)关于变量u的方程如下：

$u^{k+1}=\arg \underset{u}{\min }\{\underset{l}{\mathrm{\Σ}}{\left|\right|D{\mathrm{\α}}_l-R_{l}u\left|\right|}_2^2+{\mathrm{\μ}}_1{\left|\right|F_{p}u-f\left|\right|}_2^2+{\mathrm{\μ}}_2{\left|\right|b_1^k+\▿u-v\left|\right|}_2^2\}---\left(19\right)$

该最小二乘问题的解满足如下等式：

$(\underset{l}{\mathrm{\Σ}}{R}_l^l{R}_l+{\mathrm{\μ}}_1{F}_p^l{F}_p+{\mathrm{\μ}}_2{\▿}^T\▿)u^{k+1}=(\underset{l}{\mathrm{\Σ}}{R}_l^{l}D{\mathrm{\α}}_l^k+{\mathrm{\μ}}_1{F}_p^{l}f+{\mathrm{\μ}}_2{\▿}^l(v^k-b_1^k))---\left(20\right)$

令定义为正规化的全采样傅里叶矩阵，即F^TF＝1_N。对(20)左右两边做傅里叶变换：

$(F\underset{l}{\mathrm{\Σ}}{R}_l^T{R}_l{F}^T+{\mathrm{\μ}}_{1}F{F}_p^T{F}_p{F}^T+{\mathrm{\μ}}_{2}F{\▿}^T{F}^{T}F\▿F^T)F{u}^{k+1}=F(\underset{l}{\mathrm{\Σ}}{R}_l^{T}D{\mathrm{\α}}_l^k+{\mathrm{\μ}}_1{F}_p^{T}f+{\mathrm{\μ}}_2{\▿}^T(v^k-b_1^k))---\left(21\right)$

应用傅里叶变换的卷积定理，可得：

$u^{k+1}=F^{-1}\left(\frac{F[\underset{l}{\mathrm{\Σ}}{R}_l^{T}D{\mathrm{\α}}_l^k+{\mathrm{\μ}}_1{F}_p^{T}f+{\mathrm{\μ}}_2{\▿}^T(v^k-b_1^k)]}{F\underset{l}{\mathrm{\Σ}}{R}_l^T{R}_l{F}^T+{\mathrm{\μ}}_{1}F{F}_p^T{F}_p{F}^T+{\mathrm{\μ}}_{2}F{\▿}^T{F}^{T}F\▿F^T}\right)---\left(22\right)$

步骤172：更新变量v；

问题(17)关于变量v的方程如下：

$v^{k+1}=\arg \underset{v}{\min }\{\underset{l}{\mathrm{\Σ}}{\left|\right|B{x}_l-R_l\left(v\right)\left|\right|}_2^2+{\mathrm{\μ}}_2{\left|\right|b_1^k+\▿u^{k+1}-v\left|\right|}_2^2---\left(23\right)$

解为： $v^{k+1}=\frac{{\mathrm{\μ}}_2(b_1^k+\▿u^{k+1})+\mathrm{\τ}\underset{l}{\mathrm{\Σ}}{R}_l^{T}B{x}_l^k}{{\mathrm{\μ}}_2+\mathrm{\τ}\underset{l}{\mathrm{\Σ}}{R}_l^T{R}_l}---\left(24\right)$

步骤173：更新字典及系数：B和x_l，D和α_l，l＝1，2，…，L

$\{D^{k+1},{\mathrm{\α}}_l^{k+1}\}=\arg \underset{D,\mathrm{\Γ}}{\min }\left\{\begin{array}{c}\underset{l}{\mathrm{\Σ}}{\left|\right|D{\mathrm{\α}}_l-R_l{u}^{k+1}\left|\right|}_2^2\\ s.t.{\left|\right|{\mathrm{\α}}_l\left|\right|}_0\≤T_0,\∀l\end{array}\right\}---\left(25\right)$

$\{B^{k+1},x_l^{k+1}\}=\arg \underset{B,X}{\min }\left\{\begin{array}{c}\underset{l}{\mathrm{\Σ}}{\left|\right|B{x}_l-R_l\left(v^{k+1}\right)\left|\right|}_2^2\\ s.t.{\left|\right|x_l\left|\right|}_0\≤T_0,\∀l\end{array}\right\}---\left(26\right)$

求解(25)(26)的方法与K-SVD、DLMRI一致。

以上即为求解问题(17)的过程。

综上所述，本发明实施例提出的完整的HDLRec算法可归纳如下：

算法2：HDLRec算法，包括：

1：初始化： ${\left({\mathrm{\Γ}}^{\left(i\right)}\right)}^0=0,{\left(D^{\left(i\right)}\right)}^0,{\left(b^{\left(i\right)}\right)}^0=0,i=\mathrm{1,2};u^0=F_p^{T}f;$

2：For k＝1，2，…直到满足终止准则：

3：更新 $\{{\left(D^{\left(i\right)}\right)}^{k+1},{\left({\mathrm{\α}}_l^{\left(i\right)}\right)}^{k+1}\}$

$4:u^{k+1}=F^{-1}\left(\frac{F[\underset{l}{\mathrm{\Σ}}{R}_l^{T}D{\mathrm{\α}}_l^k+{\mathrm{\μ}}_1{F}_p^{T}f+{\mathrm{\μ}}_2{\▿}^T(v^k-b_1^k)]}{F\underset{l}{\mathrm{\Σ}}{R}_l^T{R}_l{F}^T+{\mathrm{\μ}}_{1}F{F}_p^T{F}_p{F}^T+{\mathrm{\μ}}_{2}F{\▿}^T{F}^{T}F\▿F^T}\right)$

5：更新 $\{B^k,x_l^k\}$

6： $v^{k+1}=\frac{{\mathrm{\μ}}_2(b_1^k+\▿u^{k+1})+\mathrm{\τ}\underset{l}{\mathrm{\Σ}}{R}_l^{T}B{x}_l^k}{{\mathrm{\μ}}_2+\mathrm{\τ}\underset{l}{\mathrm{\Σ}}{R}_l^T{R}_l}$

7： $b_1^{k+1}=b_1^k+\▿u^{k+1}-v^{k+1}$

8：End

9：输出u^k+1

本发明实施例磁共振快速成像的方法提出两种新的自适应字典学习模型，第一种为稀疏梯度域自适应字典学习模型(简称TVDL模型)，该模型在图像的水平梯度和垂直梯度上分别进行字典学习，可以看做TV模型的自适应策略的推广。第二种模型为图像域及其一阶梯度域自适应字典学习模型(简称HTVDL模型)，该模型在图像域及其一阶梯度域的两个方向上同时进行字典学习，可以看做是图像域字典学习的高阶模型，该模型内含了TVDL模型，因而保存了TVDL模型的优点。

本发明实施例磁共振快速成像的方法通过自适应学习字典的引入可以克服固定的有限差分变换导致目标图像的块状效应，可以处理结构更复杂的图像，从而产生更精确的重建；另外本发明实施例对梯度图像进行处理，其比原图像更稀疏，因而字典学习可以变得更精确和鲁棒，从而更稀疏地表示图像，且具有更好的保真度，且可以恢复出更多的细节。

以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。

Claims (9)

1.一种磁共振快速成像的方法，包括：

步骤A：在图像梯度域的水平梯度图像和垂直梯度图像上进行字典学习，建立图像模型；

步骤B：利用重建算法交替更新图像块的稀疏表示，恢复水平梯度和垂直梯度，然后在这两个方向梯度重建图像。

2.根据权利要求1所述的磁共振快速成像的方法，其特征在于，所述步骤A还包括：在图像域上进行字典学习，即在图像域和梯度域的水平梯度图像和垂直梯度图像上同时进行字典学习，建立图像模型。

3.根据权利要求2所述的磁共振快速成像的方法，其特征在于，在所述步骤B中，恢复图像及其水平、垂直梯度图像，然后根据图像以及水平、垂直方向梯度重建图像。

4.根据权利要求1所述的磁共振快速成像的方法，其特征在于，在所述步骤A中，建立的图像模型为：

$\underset{u,D^{\left(i\right)},{\mathrm{\Γ}}^{\left(i\right)}}{\min }\{\underset{i=1}{\overset{2}{\mathrm{\Σ}}}\underset{l}{\mathrm{\Σ}}{\left|\right|D^{\left(i\right)}{\mathrm{\α}}_l^{\left(i\right)}-R_l\left({\▿}^{\left(i\right)}u\right)\left|\right|}_2^2+\frac{v_1}{2}{\left|\right|F_{p}u-f\left|\right|}_2^2\}$

$s.t.{\left|\right|{\mathrm{\α}}_l^{\left(i\right)}\left|\right|}_0\≤T_0,\∀l,i,$

其中，第一项为梯度图像在字典上稀疏表示，第二项保证重建结果与K空间采样信号相匹配；权重v₁＝(λ/σ)，σ为测量噪声的标准差，λ为正常数。

5.根据权利要求1所述的磁共振快速成像的方法，其特征在于，所述步骤B的重建算法包括：引入辅助变量w⁽ⁱ⁾，i＝1，2，应用布雷格曼技术，定义 $\▿=\left[\begin{array}{c}{\▿}^{\left(1\right)}\\ {\▿}^{\left(2\right)}\end{array}\right],$ $b=\left[\begin{array}{c}{b}^{\left(1\right)}\\ b^{\left(2\right)}\end{array}\right],$ $w=\left[\begin{array}{c}{w}^{\left(1\right)}\\ w^{\left(2\right)}\end{array}\right].$

6.根据权利要求5所述的磁共振快速成像的方法，其特征在于，所述步骤B的重建算法还包括：更新梯度图像变量w⁽ⁱ⁾，i＝1，2；更新梯度图像块的稀疏表示；交替更新字典D⁽ⁱ⁾和系数矩阵

在稀疏编码阶段，固定字典D⁽ⁱ⁾，通过贪婪算法正交匹配追踪来更新在字典更新阶段，固定系数

通过奇异值分解来逐列更新字典的每一列，最小化近似误差。

7.根据权利要求2所述的磁共振快速成像的方法，其特征在于，所述建立的图像模型为： $\begin{array}{c}\underset{u,D,\mathrm{\ΓB},X}{\min }\{\underset{l}{\mathrm{\Σ}}{\left|\right|D{\mathrm{\α}}_l-R_{l}u\left|\right|}_2^2+\mathrm{\τ}\underset{l}{\mathrm{\Σ}}{\left|\right|B{x}_l-R_l(\▿u)\left|\right|}_2^2+{\mathrm{\μ}}_1{\left|\right|F_{p}u-f\left|\right|}^2\}\\ s.t.{\left|\right|{\mathrm{\α}}_l\left|\right|}_0\≤T_0,{\left|\right|x_l\left|\right|}_0\≤T_0,\∀l\end{array},$ 其中，

D，B＝[B₁，B₂]分别为原图像块和梯度图像块对应的字典， ${\mathrm{\α}}_l,x_l=\left[\begin{array}{c}{x}_l^1\\ x_l^2\end{array}\right]$ 为相应的系数；τ，μ₁为权重参数。

8.根据权利要求7所述的磁共振快速成像的方法，其特征在于，所述重建算法包括：引入辅助变量v，应用布雷格曼技术，将图像模型变更为：

$\{u^{k+1},v^{k+1},D^{k+1},{\mathrm{\α}}_l^{k+1},B^{k+1},x_l^{k+1}\}=\arg \underset{u,v,D,\mathrm{\Γ},B,X}{\min }\left\{\begin{array}{c}\underset{l}{\mathrm{\Σ}}{\left|\right|D{\mathrm{\α}}_l-R_{l}u\left|\right|}_2^2+\mathrm{\τ}\underset{l}{\mathrm{\Σ}}{\left|\right|B{x}_l-R_l\left(v\right)\left|\right|}_2^2+{\mathrm{\μ}}_1{\left|\right|F_{p}u-f\left|\right|}^2+{\mathrm{\μ}}_2{\left|\right|b_1^k+\▿u-v\left|\right|}_2^2\\ s.t.{\left|\right|{\mathrm{\α}}_l\left|\right|}_0\≤T_0,{\left|\right|x_l\left|\right|}_0\≤T_0,\∀l\end{array}.\right.$

9.根据权利要求8所述的磁共振快速成像的方法，其特征在于，所述重建算法还包括：更新u，固定其余的变量；更新变量v；更新字典及系数：B和x_l，D和α_l，l＝1，2，…，L，交替更新字典D⁽ⁱ⁾和系数矩阵

在稀疏编码阶段，固定字典D⁽ⁱ⁾，通过贪婪算法正交匹配追踪来更新在字典更新阶段，则固定系数

通过奇异值分解来逐列更新字典的每一列，最小化近似误差。

Images