CN103400402A

CN103400402A - 基于低秩结构稀疏的压缩感知mri图像重建方法

Info

Publication number: CN103400402A
Application number: CN2013102921374A
Authority: CN
Inventors: 董伟生; 石光明; 王春美
Original assignee: Xidian University
Current assignee: Xidian University
Priority date: 2013-07-12
Filing date: 2013-07-12
Publication date: 2013-11-20
Anticipated expiration: 2033-07-12
Also published as: CN103400402B

Abstract

本发明公开了一种基于低秩结构稀疏的压缩感知MRI图像重建方法，主要解决现有技术难以精确恢复磁共振MRI图像的问题。其实现步骤是：采用传统压缩感知初始恢复磁共振MRI图像，并在该图像中寻找相似块矩阵，形成索引集合；对相似块矩阵进行奇异值分解并计算阈值，利用该阈值对奇异值进行阈值计算得到阈值后奇异值；利用阈值后奇异值优化磁共振MRI图像，即对该磁共振MRI图像循环进行更新相似块矩阵及其索引、相似块矩阵奇异值分解以及算阈值和奇异值阈值计算的过程，得到最后的恢复图像。本发明所恢复的磁共振MRI图像更清晰，图像边缘更锐利，可用于对医学影像的处理。

Description

基于低秩结构稀疏的压缩感知MRI图像重建方法

技术领域

本发明涉及医学影像处理技术领域，具体涉及一种磁共振成像MRI的图像重建方法，主要用于对医学图像的清晰快速恢复。

背景技术

磁共振成像MRI，因其低损害性及高诊断意义而得到医学界的广泛应用，磁共振MRI建立在磁共振原理的基础上。磁共振涉及的基本物理概念主要包括：原子的自旋和磁矩，自旋磁矩在外磁场中的能量状态，产生磁共振的条件，拉莫进动，磁化强度矢量，以及射频场对磁化强度矢量和弛豫过程。

传统的磁共振成像MRI需要对原始数据K空间按照奈奎斯特采样定理进行密集采样，然后对所采集密集数据通过逆傅立叶变换重建磁共振MRI图像，这将导致重建磁共振成像所需数据量大采样时间长的难题。压缩感知理论的出现允许对原始数据K空间进行降采样，采集的原始数据K空间样本可以远小于传统磁共振MRI成像中所需采集的原始数据K空间样本数，从而大大减少采样所需数据，节约采样时间。而如何从降采样的原始数据K空间样本重建高清晰的磁共振MRI图像是压缩感知磁共振MRI成像方法成功的一个关键因素，也是近年来研究的热点。

现有磁共振压缩感知CS-MRI图像重建方法，是利用磁共振MRI图像的稀疏性来重建磁共振MRI图像，稀疏基的选取多利用全变差、DCT、小波基、以及学习得到的冗余字典对磁共振MRI图像进行稀疏表示。相比全变差、DCT以及小波基，基于学习的冗余字典有明显的优越性，能更好地刻画磁共振MRI图像中的边缘结构，但是现有的这种基于字典学习方法学习得到全局的字典难以有效表征各种磁共振MRI图像的局部结构。

此外，从欠采样的原始数据K空间恢复磁共振MRI图像的稀疏表示系数是一个病态逆问题。通过对稀疏系数存在的结构相关性进行约束可以有效提升图像稀疏分解的精度，这种结构稀疏模型已被证明是一个更加鲁棒、更精确的稀疏模型。但是，现有CS-MRI图像重建方法由于主要是利用磁共振MRI图像的稀疏性实现对图像重建，并没有利用稀疏表示系数之间存在的结构相关性，因而难以精确重建出原始真实的磁共振MRI图像，导致医学诊断困难。

发明内容

本发明的目的在于针对现有CS-MRI重建方法的不足，提出一种基于低秩结构稀疏的压缩感知磁共振MRI图像重建方法，以提升重建图像的质量。

实现本发明目的技术思路是：利用结构聚类技术对局部图像块进行结构聚类，将图像内的相似块聚为一类，对于每一个子类，利用结构稀疏模型对相似图像块进行联合稀疏编码；通过局部PCA字典学习方法构造局部自适应稀疏字典，并利用矩阵低秩逼近的方法来优化求解PCA字典学习联合结构稀疏编码问题。具体步骤包括如下：

（1）输入磁共振MRI原始K空间观测数据y，对该输入数据y进行传统压缩感知初始恢复，得到初始恢复图像x⁽¹⁾，x⁽¹⁾∈C^N，C^N表示N维复数空间；

（2）在初始恢复图像x⁽¹⁾内以第i个像素为中心，取大小为n的块

采用k近邻方式在块

周围找到与其相似的m个相似块，创建相似块矩阵

表示在初始恢复图像x⁽¹⁾内的第j_is块与

相似，生成相似块矩阵索引集合G_i=[j_i1,j_i2,...,j_im]，s=1,2,...,m，i=1,2,...,N；

（3）初始化循环次数t=1,2,...,M，M为迭代次数最大值，初始化块更新次数L，L<M，使用低秩结构稀疏算法优化重构磁共振MRI图像Z：

（3a）对相似块矩阵

进行奇异值分解，即

其中，

表示第t次奇异值分解所得相似块矩阵的左酉矩阵，

表示第t次奇异值分解所得相似块矩阵

的右酉矩阵，

Λ_{i}^{(t + 1 / 2)} = diag [λ_{j_{i 1}}^{(t + 1 / 2)}, λ_{j_{i 2}}^{(t + 1 / 2)}, . . ., λ_{j_{ip}}^{(t + 1 / 2)}, . . ., λ_{j_{ir}}^{(t + 1 / 2)}]

表示第t次奇异值分解所得相似块矩阵

的奇异值矩阵，

表示奇异值矩阵第p个奇异值元素，p=1,2,...,r，r=min(n,m)，svd表示奇异值分解符号，diag表示对角线矩阵；

（3b）计算第t次奇异值

的阈值

τ_{j_{ip}} = \frac{2 \sqrt{2} σ_{n}^{2}}{\sqrt{m} {\hat{θ}}_{j_{ip}}},

其中，

{\hat{θ}}_{j_{ip}} = \sqrt{\max (λ_{j_{ip}}^{(t + 1 / 2)} / m {- σ}_{n}^{2}, 0)},

p=1,2,...,r，σ_n为噪声方差；

（3c）利用阈值

对奇异值矩阵

内的奇异值元素进行软阈值计算，得到软阈值奇异值元素为：

λ_{j_{ip}}^{(t + 1)} = \{\begin{matrix} λ_{j_{ip}}^{(t + 1 / 2)}, & λ_{j_{ip}}^{(t + 1 / 2)} {> τ}_{j_{ip}} \\ 0, & λ_{j_{ip}}^{(t + 1 / 2)} \leq τ_{j_{ip}} \end{matrix};

（3d）由软阈值奇异值元素

作为对角元素，得到软阈值奇异值矩阵

为：

Λ_{i}^{(t + 1)} = diag [λ_{j_{i 1}}^{(t + 1)}, λ_{j_{i 2}}^{(t + 1)}, . . ., λ_{j_{ip}}^{(t + 1)}, . . ., λ_{j_{ir}}^{(t + 1)}];

（3e）根据左酉矩阵

软阈值奇异值矩阵

右酉矩阵

为构建优化图像的如下目标函数：

x^{(t + 1)} = \underset{x}{\arg \min} {| | y - Hx | |}_{2}^{2} + β \underset{i}{Σ} \underset{j &Element; G_{i}}{Σ} {| | R_{j} x_{j} - U_{i}^{(t + 1)} Λ_{j}^{(t + 1)} {(v_{j}^{(t + 1)})}^{T} | |}_{2}^{2},

其中，表示目标函数取得最小值时x所取的值，

表示2范数，H表示随机观测矩阵，x为未知真实图像，β为可调参数，R_j表示取图像x^t第j个块的矩阵，即取块矩阵，x_j为图像x^t的第j个块，

为右酉矩阵

的第j列转置；

（3f）求解步骤（3e）目标函数，得到优化的图像x^(t+1)为：

x^{(t + 1)} = \frac{[H^{T} y + β \underset{i}{Σ} \underset{j &Element; G_{i}}{Σ} R_{j}^{T} U_{i}^{(t + 1)} Λ_{j}^{(t + 1)} {(v_{j}^{(t + 1)})}^{T}]}{(H^{T} H + β \underset{i}{Σ} \underset{j &Element; G_{i}}{Σ} R_{j}^{T} R_{j})},

其中，H^T表示随机观测矩阵H的转置，

表示取块矩阵R_j的转置；

（3g）由左酉矩阵

右酉矩阵和软阈值奇异值矩阵得到优化图像块：

X_{i}^{t + 1} = U_{i}^{t + 1} Λ_{i}^{t + 1} V_{i}^{t + 1};

（3h）重复步骤（3a）-（3g）共L次，按步骤（2）方法在优化图像x^(t+1)中重新寻找相似块矩阵

更新索引G_i；

（3i）重复步骤（3a）-（3g）共M次，得到最终的优化图像x^(M+1)，x^(M+1)即为低秩结构稀疏算法优化重构的磁共振MRI图像Z。

本发明与现有技术相比具有以下优点：

第一，本发明由于在图像重建中，利用到了稀疏系数的结构相关性，能更加精确恢复磁共振MRI图像；

第二，本发明由于在参数调解中，对正则参数自适应的选取，使得重构方法更加鲁棒性；

第三，本发明提出了局部PCA字典学习方法，并提出利用矩阵低秩逼近的方法来求解PCA字典学习和联合结构稀疏编码问题，使得重构图像更加精确。

附图说明

图1为本发明的总流程图；

图2为本发明仿真实验所用的磁共振MRI人脑图像Head原图；

图3为用SparseMRI方法对人脑图像Head在观测码率为0.2时的恢复结果；

图4为用TVCMRI方法对人脑图像Head在观测码率为0.2时的恢复结果；

图5为用RecPF方法对人脑图像Head在观测码率为0.2时的恢复结果；

图6为用TV_CSA方法对人脑图像Head在观测码率为0.2时的恢复结果；

图7为用本发明方法对人脑图像Head在观测码率为0.2时的恢复结果。

具体实施方式

参照图1，本发明是基于低秩结构稀疏的压缩感知MRI图像重建方法，其实现步骤如下：

步骤1，对输入磁共振MRI原始K空间观测数据y进行传统压缩感知初始恢复，求得初始恢复图像x⁽¹⁾。

（1a）设置传统压缩感知恢复阈值Q>0，对输入磁共振MRI原始K空间观测数据y进行初始阈值计算，得初始阈值计算结果

f_{i}^{1} = \{\begin{matrix} y_{i}, & y_{i} > Q \\ 0, & y_{i} \leq Q \end{matrix} i = 1,2,3, . . ., N,

其中，y_i为观测数据y内的元素；

（1b）由初始阈值计算结果

得到初始观测数据：

（1c）利用传统压缩感知算法对磁共振MRI图像进行恢复：

（1c1）初始化最大循环次数E=300,设置传统压缩感知恢复循环次数：

e=1,2,...,E；

（1c2）生成随机观测矩阵H，设置传统压缩感知迭代步长ξ>0,本实例取ξ=1.2；

（1c3）对观测数据f^e进行反投影，得到第e次传统压缩感知恢复数据g^e：

g^e=f^e+ξ(H^T(y-A))，

其中，

为恢复数据

内的元素，H^T为随机观测矩阵H的转置，A为对观测数据f^e再观测所得数据，A=Hf^e；

（1c4）对反投影所得恢复数据g^e内的元素进行阈值计算，得第e次阈值计算结果

f_{i}^{e + 1} = \{\begin{matrix} g_{i}^{e}, & g_{i}^{e} > Q \\ 0, & g_{i}^{e} \leq Q \end{matrix};

（1c5）由阈值计算结果

得到优化的恢复数据

f^{e + 1} = [f_{1}^{e + 1}, f_{2}^{e + 1}, . . ., f_{N}^{e + 1}];

（1c6）循环步骤（1c3）-（1c5）共E次，得到最终优化的恢复数据

（1d）生成N维傅里叶变换基D，将上述恢复数据f^E+1与傅里叶基D相乘，得到磁共振MRI初始恢复图像x⁽¹⁾为：

x⁽¹⁾=Df^E+1。

步骤2，创建初始恢复图像x⁽¹⁾的相似块矩阵

生成相似块矩阵索引集合G_i。

（2a）将初始恢复图像按步长1分成大小为n的块C=[x₁,x₂,x₃,...]，标注每块在初始恢复图像x⁽¹⁾中的位置索引A=[1,2,3,...]；

（2b）在初始恢复图像x⁽¹⁾内以第i个像素为中心，取大小为n的块

（2c）采用k近邻方式在块

周围找到与其相似的m个相似块，创建相似块矩阵：

X_{i}^{1} = [x_{j_{i 1}}, x_{j_{i 2}}, . . ., x_{j_{is}}, . . ., x_{j_{im}}],

s=1,2,...,m；

其中，

表示在初始恢复图像x⁽¹⁾内的第j_is块与

相似，

j_is为属于相似块矩阵

的位置索引，j_is∈A；

（2d）根据相似块矩阵索引j_is，生成相似块矩阵

的索引：

G_i=[j_i1,j_i2,...,j_is,...,j_im]。

步骤3，使用低秩结构稀疏算法优化重构磁共振MRI图像Z。

（3a）设置循环次数t=1,2,....,M，设置相似块矩阵更新次数1≤L≤M，设置最大循环次数M≥200，本实例取L=18，M=300；

（3b）对相似块矩阵

进行奇异值分解，即

其中，

表示第t次奇异值分解所得相似块矩阵

的左酉矩阵，表示第t次奇异值分解所得相似块矩阵

的右酉矩阵，

Λ_{i}^{(t + 1 / 2)} = diag [λ_{j_{i 1}}^{(t + 1 / 2)}, λ_{j_{i 2}}^{(t + 1 / 2)}, . . ., λ_{j_{ip}}^{(t + 1 / 2)}, . . ., λ_{j_{ir}}^{(t + 1 / 2)}]

表示第t次奇异值分解所得相似块矩阵

的奇异值矩阵，

表示奇异值矩阵

第p个奇异值元素，p=1,2,...,r，r=min(n,m)，svd表示奇异值分解符号，diag表示对角线矩阵；

（3c）计算第t次奇异值

的阈值

τ_{j_{ip}} = \frac{2 \sqrt{2} σ_{n}^{2}}{\sqrt{m} {\hat{θ}}_{j_{ip}}},

其中，

{\hat{θ}}_{j_{ip}} = \sqrt{\max (λ_{j_{ip}}^{(t + 1 / 2)} / m {- σ}_{n}^{2}, 0)},

p=1,2,...,r，σ_n为噪声方差；

（3d）利用阈值

对奇异值矩阵

内的奇异值元素

进行软阈值计算，得到软阈值奇异值元素

为：

λ_{j_{ip}}^{(t + 1)} = \{\begin{matrix} λ_{j_{ip}}^{(t + 1 / 2)}, & λ_{j_{ip}}^{(t + 1 / 2)} {> τ}_{j_{ip}} \\ 0, & λ_{j_{ip}}^{(t + 1 / 2)} \leq τ_{j_{ip}} \end{matrix};

（3e）由软阈值奇异值元素作为对角元素，得到软阈值奇异值矩阵

为：

Λ_{i}^{(t + 1)} = diag [λ_{j_{i 1}}^{(t + 1)}, λ_{j_{i 2}}^{(t + 1)}, . . ., λ_{j_{ip}}^{(t + 1)}, . . ., λ_{j_{ir}}^{(t + 1)}];

（3f）根据左酉矩阵软阈值奇异值矩阵

右酉矩阵

构建优化图像的目标函数x^(t+1)：

（3f1）由磁共振MRI图像的非局部自相似性，构建压缩感知图像恢复的先验项为：

\underset{i}{Σ} \underset{j &Element; G_{i}}{Σ} {| | R_{j} x_{j} - U_{i}^{(t + 1)} Λ_{j}^{(t + 1)} {(v_{j}^{(t + 1)})}^{T} | |}_{2}^{2};

其中，

表示2范数，R_j表示取图像x^t第j个块的矩阵，即取块矩阵，x_j为图像x^t的第j个块，

为右酉矩阵

的第j列转置；

（3f2）将先验项与压缩感知图像恢复的似然项

结合，得到目标函数：

x^{(t + 1)} = \underset{x}{\arg \min} {| | y - Hx | |}_{2}^{2} + β \underset{i}{Σ} \underset{j &Element; G_{i}}{Σ} {| | R_{j} x_{j} - U_{i}^{(t + 1)} Λ_{j}^{(t + 1)} {(v_{j}^{(t + 1)})}^{T} | |}_{2}^{2},

其中，

表示目标函数取得最小值时x所取的值，x为未知真实图像，β为可调参数，本实例设置β=1.3；

（3g）求解步骤（3f2）的目标函数，得到优化的图像x^(t+1)为：

x^{(t + 1)} = \frac{[H^{T} y + β \underset{i}{Σ} \underset{j &Element; G_{i}}{Σ} R_{j}^{T} U_{i}^{(t + 1)} Λ_{j}^{(t + 1)} {(v_{j}^{(t + 1)})}^{T}]}{(H^{T} H + β \underset{i}{Σ} \underset{j &Element; G_{i}}{Σ} R_{j}^{T} R_{j})},

其中，H^T表示随机观测矩阵H的转置，

表示取块矩阵R_j的转置；

（3h）由左酉矩阵

右酉矩阵和软阈值奇异值矩阵

得到优化图像块:

X_{i}^{t + 1} = U_{i}^{t + 1} Λ_{i}^{t + 1} V_{i}^{t + 1};

（3i）重复步骤（3b）-（3h）共L次，按步骤（2）方法在优化图像x^(t+1)中重新寻找相似块矩阵更新索引G_i；

（3j）重复步骤（3b）-（3h）共M次，得到最终的优化图像x^(M+1)，x^(M+1)即为低秩结构稀疏算法优化重构的磁共振MRI图像Z。

本发明的效果可以通过如下仿真实验具体说明：

1.仿真条件：

1）仿真实验观测矩阵采用随机傅里叶观测矩阵；

2）仿真实验所用编程平台为MatlabR2010b；

3)仿真实验所用人脑图像Head来源于飞利浦公司，如图2所示；

4）仿真试验中，采用峰值信噪比PSNR指标来评价压缩感知实验结果，峰值信噪比PSNR定义为：

PSNR = {10 \log}_{10} (\frac{255^{2}}{Σ {| | x - Z | |}^{2}})

其中，x为原图像，Z为恢复出来的图像。

2.仿真内容：

仿真1，采用SparseMRI方法，对人脑图像Head在观测码率为0.2时进行恢复，其恢复结果如图3所示；

仿真2，采用TVCMRI方法，对人脑图像Head在观测码率为0.2时进行恢复，其恢复结果如图4所示；

仿真3是采用RecPF方法，对人脑图像Head在观测码率为0.2时进行恢复，其恢复结果如图5所示；

仿真4，采用TV_CSA方法，对人脑图像Head在观测码率为0.2时进行恢复，其恢复结果如图6所示；

仿真5，采用本发明方法对人脑图像Head在观测码率为0.2时进行恢复。其恢复结果如图7所示。

从图3-图7所显示的人脑图像Head的恢复结果可以看出，本发明的低秩结构稀疏的压缩感知磁共振MRI重建方法恢复出来的图像比其他方法恢复出来的图像更干净，清晰，图像边缘更锐利，视觉效果更好。

将SparseMRI方法、TVCMRI方法、RecPF方法、TV_FCSA方法和本发明方法分别对人脑图像Head进行压缩感知图像重建仿真，得到的峰值信噪比PSNR结果见表1。

表1恢复图像的峰值信噪比PSNR值(单位dB)

从表1可以看出，本发明的峰值信噪比PSNR比SparseMRI和TVCMRI在观测码率为0.2的时候要平均高出15.86B和12.19dB，比TV_FCSA高出5.45dB。

Claims

1.一种基于低秩结构稀疏的压缩感知MRI图像重建方法，包括如下步骤：

采用k近邻方式在块

周围找到与其相似的m个相似块，创建相似块矩阵

表示在初始恢复图像x⁽¹⁾内的第j_is块与

（3a）对相似块矩阵进行奇异值分解，即

svd (X_{i}^{t}) = (U_{i}^{(t + 1)}, Λ_{i}^{(t + 1 / 2)}, V_{i}^{(t + 1)}),

其中，

表示第t次奇异值分解所得相似块矩阵

的左酉矩阵，

表示第t次奇异值分解所得相似块矩阵

的右酉矩阵，

Λ_{i}^{(t + 1 / 2)} = diag [λ_{j_{i 1}}^{(t + 1 / 2)}, λ_{j_{i 2}}^{(t + 1 / 2)}, . . ., λ_{j_{ip}}^{(t + 1 / 2)}, . . ., λ_{j_{ir}}^{(t + 1 / 2)}]

表示第t次奇异值分解所得相似块矩阵

的奇异值矩阵，表示奇异值矩阵

（3b）计算第t次奇异值的阈值

τ_{j_{ip}} = \frac{2 \sqrt{2} σ_{n}^{2}}{\sqrt{m} {\hat{θ}}_{j_{ip}}},

其中，

{\hat{θ}}_{j_{ip}} = \sqrt{\max (λ_{j_{ip}}^{(t + 1 / 2)} / m {- σ}_{n}^{2}, 0)},

p=1,2,...,r，σ_n为噪声方差；

（3c）利用阈值

对奇异值矩阵

内的奇异值元素

进行软阈值计算，得到软阈值奇异值元素

为：

λ_{j_{ip}}^{(t + 1)} = \{\begin{matrix} λ_{j_{ip}}^{(t + 1 / 2)}, & λ_{j_{ip}}^{(t + 1 / 2)} {> τ}_{j_{ip}} \\ 0, & λ_{j_{ip}}^{(t + 1 / 2)} \leq τ_{j_{ip}} \end{matrix};

（3d）由软阈值奇异值元素

作为对角元素，得到软阈值奇异值矩阵为：

Λ_{i}^{(t + 1)} = diag [λ_{j_{i 1}}^{(t + 1)}, λ_{j_{i 2}}^{(t + 1)}, . . ., λ_{j_{ip}}^{(t + 1)}, . . ., λ_{j_{ir}}^{(t + 1)}];

（3e）根据左酉矩阵

软阈值奇异值矩阵右酉矩阵为构建优化图像的如下目标函数：

x^{(t + 1)} = \underset{x}{\arg \min} {| | y - Hx | |}_{2}^{2} + β \underset{i}{Σ} \underset{j &Element; G_{i}}{Σ} {| | R_{j} x_{j} - U_{i}^{(t + 1)} Λ_{j}^{(t + 1)} {(v_{j}^{(t + 1)})}^{T} | |}_{2}^{2},

其中，

表示目标函数取得最小值时x所取的值，

为右酉矩阵

的第j列转置；

（3f）求解步骤（3e）目标函数，得到优化的图像x^(t+1)为：

x^{(t + 1)} = \frac{[H^{T} y + β \underset{i}{Σ} \underset{j &Element; G_{i}}{Σ} R_{j}^{T} U_{i}^{(t + 1)} Λ_{j}^{(t + 1)} {(v_{j}^{(t + 1)})}^{T}]}{(H^{T} H + β \underset{i}{Σ} \underset{j &Element; G_{i}}{Σ} R_{j}^{T} R_{j})},

其中，H^T表示随机观测矩阵H的转置，表示取块矩阵R_j的转置；

（3g）由左酉矩阵右酉矩阵

和软阈值奇异值矩阵

得到优化图像块：

X_{i}^{t + 1} = U_{i}^{t + 1} Λ_{i}^{t + 1} V_{i}^{t + 1};

更新索引G_i;

2.根据权利要求1所述的重建方法，其特征在于所述步骤（1）中对输入数据y进行传统压缩感知初始恢复，得到初始恢复图像x⁽¹⁾，通过如下步骤进行：

（1a）求解输入数据y在傅里叶基D下的稀疏系数

\hat{a} = \arg \min_{a} {| | y - HDa | |}_{2}^{2} + λ {| | a | |}_{1},

其中，

表示目标函数取得最小值时a所取的值，a为初始恢复磁共振MRI图像真实傅里叶稀疏系数，

表示2范数，H表示随机观测矩阵，λ为可调正则参数。

（1b）将上述稀疏系数

与傅里叶基D相乘，得到磁共振MRI初始恢复图像x⁽¹⁾为：

x^{(1)} = D \hat{a} .