CN107301630B - 一种基于排序结构组非凸约束的cs-mri图像重构方法 - Google Patents

Abstract

本发明公开了一种基于排序结构组非凸约束的CS‑MRI图像重构方法。属于数字图像处理技术领域。它是一种利用结构组排序提高固定字典表示能力并利用log‑sum范数来对结构组进行非凸约束的图像重构方法。首先找到目标图像块的相似图像块集合即结构组，然后对目标图像块建立排序模型并将获得的排序矩阵用于结构组的排序，以提高固定字典对结构组的表达能力，最后利用log‑sum范数对稀疏系数进行非凸约束；本发明对结构组内部进行排序，提高了结构组稀疏表示的性能，并采用一种快速阈值算子对系数的非凸优化进行求解，使估计出的系数更接近真实值，通过本发明得到的图像整体更加清晰，并保留了图像的大量细节信息，恢复的准确度更高，因此可用于医学图像的重构。

Description

一种基于排序结构组非凸约束的CS-MRI图像重构方法

技术领域

本发明属于数字图像处理技术领域，它特别涉及在变换域对图像进行增强稀疏表示以及非凸约束来对图像重构的方法，用于医学图像高质量恢复。

背景技术

磁共振成像(MRI)由于其高分辨率和无创性在临床医学诊断中有着非常广泛的应用。传统的磁共振成像需要对原始数据进行奈奎斯特采样，所需的时间较长，成本较高，因此一定程度上限制了该技术在医学上的应用。

随着近年来压缩感知(CS)理论的提出，磁共振成像在减少成像时间上有了重大突破的可能。压缩感知理论提出可以利用信号的稀疏特性，在远小于奈奎斯特采样率的情况下，将采样信号利用非线性算法进行重建。磁共振图像本身存在数据的冗余，而在变换域上具有稀疏特性，这就给压缩感知能从降采样的原始数据K空间样本精确重建出原始信号提供了必要条件。而如何利用该降采样样本重构出更清晰的磁共振图像也成为这一过程中一个非常关键的环节，这也是图像重构中的一个研究热点。

由于图像在变换域的稀疏度越高，重构的图像效果越好。一些图像重构方法选取不同的稀疏字典(如全变分，离散余弦，小波等)来对图像进行稀疏表示，但这些固定字典缺乏对图像的自适应性，对图像的细节保留能力有限。随后基于学习的字典(如KSVD)被用到MRI图像重构中，且取得了良好的结果，但这种基于全局冗余学习对图像块进行编码的字典复杂度较高，且容易出现过拟合现象。

最近的研究表明，图像存在着非局部相似性，即不同区域之间有很多相似的结构。为利用这一特性进一步提高稀疏度，需要以相似图像块集合为处理对象进行稀疏表示，使重构的图像细节更加丰富。

发明内容

本发明的目的在于针对现有CS-MRI图像重构方法存在的不足，提出一种基于排序结构组非凸约束的CS-MRI图像重构方法。该方法充分考虑了图像在变换域的稀疏特性，在考虑图像块间非局部相似性的基础上，将图像块内像素进行排序，增强了字典对结构组的表示能力，同时在对结构组稀疏系数进行估计的过程中，使用一种快速阈值算子，使估计出的稀疏系数更接近真实值，从而使整个图像更清晰，同时能更多保留图像的细节信息。具体包括以下步骤：

(1)输入一幅MRI原始K空间观测数据，对输入数据y进行传统压缩感知初始重构，得到初始重构图像x⁽⁰⁾；

(2)在重构图像x⁽⁰⁾内对每个目标图像块x_i进行相似图像块搜索，并使目标图像块与其相似度最高的S-1个图像块组成对应的相似图像块集合即结构组X_i＝[x_i,0,x_i,1,…x_i,S-1]，其中x_i,0＝x_i；

(3)对每个结构组的目标图像块x_i建立排序模型，通过该模型得到其对应的排序矩阵P_i；

(4)对每个结构组X_i内的图像块的像素利用排序矩阵P_i进行排序，再利用排序结构组具有的稀疏性和非局部相似性进行离散小波和离散余弦变换，并将变换后的系数矩阵向量化以获得稀疏系数z_i；

(5)以结构组为处理对象建立非凸约束下的MRI图像重构模型：

其中F_U为降采样傅里叶编码矩阵，R_i为抽取图像块矩阵，表示从整个图像中抽取对应结构组x_i，B_i为对结构组向量

进行两次变换的矩阵，即

Φ表示离散小波变换，Ψ表示离散余弦变换，

表示克罗内克积，

为B_i的共轭转置，N为结构组的数量，λ和β为正则化参数，C(z_i,ε)为log-sum的非凸约束项，ε为避免数值不稳定问题而引入的一个较小的正数，然后逐渐增大β值，并利用交替方向算法(ADM)对整个重构模型进行求解：

(5a)对于模型中的变量

给定x和β，则重构模型变为求解关于每个结构组稀疏系数z_i的子问题：

该模型可利用一种快速阈值算子来求解；

(5b)在得到结构组稀疏系数估计值z_i后，关于x的重构模型为：

该模型为最小二乘模型，可用共轭梯度法来求解得到

(6)重复步骤(2)～(5)，直到估计图像满足条件或迭代次数达到预设上限。

本发明的创新点是利用图像局部稀疏性和非局部相似性对结构组在变换域进行稀疏表示；为了增强字典对结构组的表示性能，在稀疏表示之前对结构组进行排序；使用log-sum非凸约束项来对稀疏系数进行约束，来更好地估计真实图像的稀疏系数；在估计稀疏系数的过程中采用一种快速阈值算子，使估计出的稀疏系数更接近真实值，并将该方法应用于核磁共振图像(MRI)重构。

本发明的有益效果：对结构组进行排序，增强了字典对结构组的稀疏表示能力，充分利用结构组具有的图像块内块间相似性，对结构组进行左右两次变换以提高稀疏度；采用log-sum对稀疏系数进行非凸约束，并采用一个快速阈值算子实现稀疏系数的估计，因此最终估计的图像不仅整体视觉效果良好，还保留了图像内部大量细节，使整个估计结果更接近真实值。

本发明主要采用仿真实验的方法进行验证，所有步骤、结论都在MATLAB8.0上验证正确。

附图说明

图1是本发明的工作流程框图；

图2是本发明仿真中使用的MRI人体心脏图像原图；

图3是用各方法(依次为PBDW方法、PANO方法、NLR方法和本发明方法)对采样率为10％的人体心脏图像的重构结果；

图4是对应的各方法(依次为PBDW方法、PANO方法、NLR方法和本发明方法)对采样率为10％的人体心脏图像的重构结果的误差。

具体实施方式

参照图1，本发明是基于排序结构组非凸约束的CS-MRI图像重构方法，具体步骤包括如下：

步骤1，对图像进行初始恢复，并建立每个图像块对应的结构组。

(1a)输入一幅MRI原始K空间观测数据y，利用总变分法对其进行初始重构，得到初始重构图像x⁽⁰⁾；

(1b)将图像按照大小为

的图像块进行抽取，并对每个目标图像块x_i与在搜索范围内的其他图像块进行欧式距离比较；

(1c)取出与目标图像块x_i欧氏距离最小的S-1个图像块，并与目标图像块组成结构组X_i＝[x_i,0,x_i,1,…x_i,S-1]，其中x_i,0＝x_i。

步骤2，由于结构组内所有相似图像块均与目标图像块高度相似，为寻找结构组的排序矩阵，可以结构组内的目标图像块为排序对象，并将其获得的排序矩阵作为整个结构组的排序矩阵，再对排序后的结构组进行稀疏表示。

(2a)以结构组内目标图像块x_i为对象建立排序模型：

其中P_i为排序矩阵，

表示每次迭代后目标图像块x_i的重构结果，μ为平衡这两项的参数，

为P_i的共轭转置，Φ是一个离散小波字典；

(2b)对于(2a)中的α_i，可直接用硬阈值求解得到：

其中l为迭代次数，硬阈值算子Η_μ为：

(2c)通过(2b)求得α_i后，关于P_i的子问题为：

将该目标函数展开可得：

其中

为常量，

第三项的上界为：

其中

p和q为分别将b和

按元素数值大小排序的序列，故排序矩阵P_i对应的排序序列为θ＝q(p^-1)，由此可求得P_i。

(2d)求得排序矩阵P_i后，将结构组X_i向量化得到

然后对其进行变换得到稀疏系数z_i：

其中Φ表示离散小波变换，Ψ表示离散余弦变换，

表示克罗内克积。

步骤3，建立非凸约束模型并估计系数。

(3a)以结构组为处理对象建立非凸约束下的MRI图像重构模型：

其中F_U为降采样傅里叶编码矩阵，R_i为抽取图像块矩阵，表示从整个图像中抽取对应结构组x_i，N为结构组的数量，λ和β为正则化参数，C(z_i,ε)为log-sum的非凸约束项，ε为避免数值不稳定问题而引入的一个较小的正数；

(3b)对于(3a)中的变量

给定x和β，则模型变为求解关于每个结构组稀疏系数z_i的子问题：

该模型可利用一种快速阈值算子来求解：

3b1)令

τ＝1/β，可知

的下界为：

3b2)将3b1)中的

取到下界时的

代入关于每个结构组稀疏系数z_i的子问题可得：

3b3)对于3b2)中任意一项

可用以下快速阈值算子求解：

其中阈值函数为：

以及

其中

由此可获得整个模型的求解结果；

(3c)通过(3b)估计出稀疏系数

后，(3a)中的目标函数变为：

该模型为最小二乘模型，可用共轭梯度法来避免求解逆矩阵得到

步骤4，重复步骤(2)～(3)，直到得到的估计图像满足条件或迭代次数达到预设上限。

本发明的效果可以通过以下仿真实验进一步说明：

一、实验条件和内容

实验条件：实验使用伪径向采样矩阵；实验图像采用真实人体心脏图像，如图2所示；实验结果评价指标采用峰值信噪比PSNR和高频误差范数HFEN来客观评价重构结果，其中HFEN定义为：

其中LoG(·)表示拉普拉斯高斯滤波器，来抽取代表图像细节信息的高频部分，

为重构结果，x为原始图像。

实验内容：在上述实验条件下，重构结果使用目前在CS-MRI图像重构领域具有代表性的PBDW方法、PANO方法和NLR方法与本发明方法进行对比。

实验1：用本发明方法和PBDW方法、PANO方法和NLR方法分别对图2采样后的图像进行重构。其中PBDW方法首先利用方向小波对图像块进行排序，并采用l₁范数进行系数约束以实现MRI图像重构，其重构结果为图3(a)，重构误差为图4(a)；PANO方法一种典型的对结构组进行三维小波变换并用l₁范数约束稀疏系数的重构方法，其重构结果为图3(b)，重构误差为图4(b)；而NLR方法利用结构组的低秩特性，并采用logdet(·)作为结构组的非凸约束项，其重构结果为图3(c)，重构误差为图4(c)。实验中本发明方法设置图像块大小

结构组内图像块个数S＝32，最大迭代次数Τ＝100，迭代终止系数η＝5×10^-8；最终重构结果为图3(d)，重构误差为图4(d)。

从图3各方法的重构结果及局部区域放大图可以看出，对比PBDW方法、PANO方法、NLR方法与本发明方法可以看出，本发明方法在重构结果的细节部分高于其他对比方法，图4的重构误差可以更直观地表现各方法的重构差异，其中如颜色参照图所示，误差图中颜色越暖表示误差越大，颜色越冷表示误差越小，从各方法的误差图中，可以直观地看出本发明方法的重构误差最小，因此重构结果更好。

表1不同重构方法的PSNR指标

图像	PBDW方法	PANO方法	NLR方法	本发明方法
					心脏图	32.29	32.75	33.34	34.29

表1给出了各方法重构结果的PSNR指标情况，其中PSNR值越高表示重构效果越好；可以看出本发明方法对比其他方法均有较大提高，此结果与重构效果图相吻合。

表2不同重构方法的HFEN指标

图像	PBDW方法	PANO方法	NLR方法	本发明方法
					心脏图	1.1504	1.1599	1.1278	0.9989

表2给出了各方法重构结果的HFEN指标情况，其中HFEN值越低表示重构结果细节保留更好；可以看出本发明方法相较其他方法，本发明方法能保留更多的细节信息，此结果与重构效果图相吻合。

上述实验表明，本发明重构方法不仅还原效果明显，而且重构后图像内容丰富，同时视觉效果及客观评价指标都较好，由此可见本发明对医学图像重构是有效的。

Claims (4)

1.一种基于排序结构组非凸约束的CS-MRI图像重构方法，包括以下步骤：

(1)输入一幅MRI原始K空间观测数据，对输入数据y进行传统压缩感知初始重构，得到初始重构图像x⁽⁰⁾；

(2)在重构图像x⁽⁰⁾内对每个目标图像块x_i进行相似图像块搜索，并使目标图像块与其相似度最高的S-1个图像块组成对应的相似图像块集合即结构组X_i＝[x_i,0,x_i,1,…x_i,S-1]，其中x_i,0＝x_i；

(3)对每个结构组的目标图像块x_i建立排序模型

其中P_i为排序矩阵，

表示每次迭代后目标图像块x_i的重构结果，μ为平衡这两项的参数，

为P_i的共轭转置，Φ是一个离散小波字典，α_i是目标图像块x_i的稀疏系数，通过该模型得到排序矩阵P_i；

(4)对每个结构组X_i内的图像块的像素利用排序矩阵P_i进行排序，再利用排序结构组具有的稀疏性和非局部相似性进行离散小波和离散余弦变换，并将变换后的系数矩阵向量化以获得稀疏系数z_i；

(5)以结构组为处理对象建立非凸约束下的MRI图像重构模型：

其中F_U为降采样傅里叶编码矩阵，R_i为抽取图像块矩阵，表示从整个图像中抽取对应结构组x_i，B_i为对结构组向量

进行两次变换的矩阵，即

Φ表示离散小波变换，Ψ表示离散余弦变换，

表示克罗内克积，

为B_i的共轭转置，N为结构组的数量，λ和β为正则化参数，C(z_i,ε)为log-sum的非凸约束项，ε为避免数值不稳定问题而引入的一个较小的正数，然后逐渐增大β值，并利用交替方向算法(ADM)对整个重构模型进行求解：

(5a)对于模型中的变量

给定x和β，则重构模型变为求解关于每个结构组稀疏系数z_i的子问题：

该模型可利用一种快速阈值算子来求解；

(5b)在得到结构组稀疏系数估计值z_i后，关于x的重构模型为：

该模型为最小二乘模型，可用共轭梯度法来求解得到

(6)重复步骤(2)～(5)，直到估计图像满足条件或迭代次数达到预设上限。

2.根据权利要求1所述的一种基于排序结构组非凸约束的CS-MRI图像重构方法，其特征在于，步骤(3)中对目标图像块x_i建立排序模型得到对应排序矩阵P_i，通过求解以下模型得到：

其中

表示每次迭代后目标图像块x_i的重构结果，μ是用来平衡这两项的参数，

为P_i的共轭转置，Φ是一个离散小波字典，该模型可转化为分别求解关于P_i和α_i的两个子问题：

(3a)对于α_i，可直接用硬阈值求解得到：

其中l为迭代次数，硬阈值算子Η_μ为：

(3b)求得α_i后，关于P_i的子问题为：

将该目标函数展开可得：

其中

为常量，

的上界为：

其中

p和q为分别将b和

按元素绝对值大小的排序序列，故排序矩阵P_i对应的排序序列为θ＝q(p^-1)，由此可求得P_i。

3.根据权利要求1所述的一种基于排序结构组非凸约束的CS-MRI图像重构方法，其特征在于，步骤(5)中的log-sum非凸约束项，其具体形式为：

其中m为图像块内元素个数，这种log-sum非凸约束项能更好逼近l₀范数，使约束的结果更精确。

4.根据权利要求1所述的一种基于排序结构组非凸约束的CS-MRI图像重构方法，其特征在于，步骤(5a)中的模型求解问题，可按照以下步骤得到：

5a1)令

τ＝1/β，可知

的下界为：

5a2)将5a1)中

取到下界时的

代入关于每个结构组稀疏系数z_i的子问题可得：

5a3)对于5a2)中任意一项

可用以下快速阈值算子求解：

其中阈值函数为：

以及

其中

由此可实现整个模型的求解。

Images