CN104933683B - 一种用于磁共振快速成像的非凸低秩重建方法 - Google Patents
一种用于磁共振快速成像的非凸低秩重建方法 Download PDFInfo
- Publication number
- CN104933683B CN104933683B CN201510310636.0A CN201510310636A CN104933683B CN 104933683 B CN104933683 B CN 104933683B CN 201510310636 A CN201510310636 A CN 201510310636A CN 104933683 B CN104933683 B CN 104933683B
- Authority
- CN
- China
- Prior art keywords
- image
- low
- rank
- convex
- matrix
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Expired - Fee Related
Links
- 238000000034 method Methods 0.000 title claims abstract description 54
- 238000003384 imaging method Methods 0.000 title claims abstract description 19
- 238000004422 calculation algorithm Methods 0.000 claims abstract description 65
- 239000011159 matrix material Substances 0.000 claims abstract description 60
- 238000000354 decomposition reaction Methods 0.000 claims abstract description 16
- 238000000926 separation method Methods 0.000 claims abstract description 3
- 238000012545 processing Methods 0.000 claims description 2
- 230000015572 biosynthetic process Effects 0.000 claims 1
- 238000006467 substitution reaction Methods 0.000 claims 1
- 238000013178 mathematical model Methods 0.000 abstract description 7
- 230000008569 process Effects 0.000 abstract description 6
- 230000009466 transformation Effects 0.000 abstract description 3
- 238000002595 magnetic resonance imaging Methods 0.000 description 114
- 238000005070 sampling Methods 0.000 description 54
- 230000006870 function Effects 0.000 description 17
- 230000000694 effects Effects 0.000 description 6
- 238000004458 analytical method Methods 0.000 description 5
- 238000004364 calculation method Methods 0.000 description 5
- 238000011161 development Methods 0.000 description 5
- 238000005259 measurement Methods 0.000 description 5
- 210000003484 anatomy Anatomy 0.000 description 4
- 238000011160 research Methods 0.000 description 4
- 238000003745 diagnosis Methods 0.000 description 3
- 238000005516 engineering process Methods 0.000 description 3
- 238000002474 experimental method Methods 0.000 description 3
- 238000013459 approach Methods 0.000 description 2
- 238000005457 optimization Methods 0.000 description 2
- 230000003044 adaptive effect Effects 0.000 description 1
- 239000000654 additive Substances 0.000 description 1
- 230000000996 additive effect Effects 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000003247 decreasing effect Effects 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 238000013461 design Methods 0.000 description 1
- 238000002059 diagnostic imaging Methods 0.000 description 1
- 201000010099 disease Diseases 0.000 description 1
- 208000037265 diseases, disorders, signs and symptoms Diseases 0.000 description 1
- 238000001914 filtration Methods 0.000 description 1
- 238000007429 general method Methods 0.000 description 1
- 230000006872 improvement Effects 0.000 description 1
- 230000010365 information processing Effects 0.000 description 1
- 230000014759 maintenance of location Effects 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000013421 nuclear magnetic resonance imaging Methods 0.000 description 1
- 210000000056 organ Anatomy 0.000 description 1
- 230000005855 radiation Effects 0.000 description 1
- 238000011084 recovery Methods 0.000 description 1
- 238000012552 review Methods 0.000 description 1
- 238000012772 sequence design Methods 0.000 description 1
- 238000004088 simulation Methods 0.000 description 1
- 238000001228 spectrum Methods 0.000 description 1
- 238000012549 training Methods 0.000 description 1
Landscapes
- Magnetic Resonance Imaging Apparatus (AREA)
Abstract
一种用于磁共振快速成像的非凸低秩重建方法，本发明建立基于非局部相似图像块低秩先验信息的MR图像数据重构数学模型，釆用交替方向迭代方法对模型进行迭代求解：通过泰勒一阶近似和奇异值分解迭代求解带低秩先验信息的非局部图像模型的低秩矩阵的非凸p范数，得到相似图像块，再通过增加辅助变量和分离变量法迭代求解重建图像。本发明利用图像先验信息，将图像块的非局部相似性与低秩特性结合，并利用傅里叶变换和低秩矩阵的特性简化计算过程，降低算法复杂度，提高了部分K空间数据重建MRI图像的性能，在更少的扫描测量下更精确地重建图像，减少重建图像的伪影，实现快速磁共振成像。
Description
技术领域
本发明属于医学成像领域，尤其涉及磁共振成像。
背景技术
磁共振成像(MRI)技术是一种能取得活体器官和组织详细诊断图像的医疗诊断技术，具有无损伤无辐射等优点，在医学临床与科研领域得到广泛应用。核磁共振成像技术能给医生提供更加清晰和更高对比度的人体结构医学图像，其一诞生就受到临床医生的欢迎，迅速在临床上得到应用，成为一些疾病诊断不可或缺的检查方式。然而MRI具有成像速度慢的不足，为减少MRI成像时间，目前主要有两种途径：一类是对硬件设备进行改进，如多线圈并行成像、快速成像梯度序列设计等；一类是通过减少K空间(Fourier频谱数据组成)采样数据量，再由相应的软件算法重建，即部分K空间重建方法。部分K空间重建由于无需对硬件进行改进，仅通过软件算法即可提高成像速度，近年来备受关注。
压缩感知(Compressed Sensing,CS)最早于2004年提出。压缩感知理论刚提出就受到学术界和工业界的广泛关注，在信号与信息处理领域得到广泛应用，并于2007被美国科技评论评为当年度的十大科技进展之一。压缩感知理论是利用信号的稀疏性，即在稀疏域中，少部分信号为非零，其余信号值为零或接近零，所以只存储少量的非零信号就可以完全重构原信号。为了取得稀疏系数，需要利用传感矩阵对信号进行变换，即可取得观测信号，观测信号的数据量比原信号小很多。因此为从观测信号中重构原始信号，需构造目标函数求解原信号。简言之压缩感知理论是指在某些变换域里具有稀疏性质的图像可以从远小于奈奎斯特采样率的测量数据中得到精确重建。特别地，在磁共振领域中的压缩感知理论被称为CSMRI。
随着CS理论的发展，研究人员开始更多地关注MR图像的不同稀疏表示方式，希望将稀疏表示领域的相关研究应用到MRI。MR图像的稀疏性是重构K空间数据的一个重要前提，而在压缩感知MRI系统的设计中稀疏表示问题可以归结为怎样表达MR图像的稀疏先验。常见的MR图像稀疏先验有MR图像在小波域、离散余弦域、总变差域的稀疏性，其他还有MR图像在自适应训练字典域的稀疏性。大部分方法可以归结为两类：一种是基于已有的数学模型建立一个稀疏表示的字典来表示数据，例如小波基底和离散余弦基底。另一种是根据样本学习一个字典来重建数据，例如2011年，Ravishankar等人利用K-SVD方法来训练字典，然后重构MR图像。DLMRI的重建效果比Lustig等人的SparseMRI方法有较大改善，但计算量非常大。近年来随着对低秩矩阵的稀疏性研究不断深入，越来越期望挖掘MR图像非局部相似性以得到某种形式的稀疏性。
A Buades等人提出非局部均值去噪方法(NLM-based method)后，学者们开始重视图像块的稀疏性及相似性。K.Dabov等人则提出了另一种基于块稀疏的研究方法，他们在对局部块相似性研究的基础上，提出了一种新的块匹配重建算法(BM3D)，该算法利用局部块的相似性对局部块聚类，采用滤波的方法对图像进行稀疏重建。Akcakaya.M等人则提出了基于块匹配的MRI重建算法模型对中心采样的MR图像进行重建，并实现了在4倍下采样的MR图像的精确重建(LOST)。近年来发展起来的非局部均值方法重构MRI下采样数据，同样利用块之间的相似性及冗余性，对下釆样的MRI数据进行重构，实验表明这些采用非局部均值方法的结果在MRI图像细节的保真方面会做的比传统CS重建方法好。
Avishankar等人提出了一种两步迭代方法，将字典学习模型用到K空间欠采样的磁共振图像重建上，提出字典学习重建磁共振图像(DLMRI)模型：
此处，Γ＝[α1，α2，...，αL]为所有图像块对应的稀疏稀疏矩阵。DLMRI直接通过正交匹配追踪方法解决l0最小化问题，虽然这些数据的学习方法比以前预定义为基础字典的方法有很大的改善，但l0最小化问题是非凸的NP难问题。模型中前一项用于图像块在自适应学习字典上的稀疏表示，后一项是在图像数据的拟合保真项。求解该模型的两步迭代方法，第一步是自适应字典学习；第二步是从高度欠采样的K-空间数据重建图像。
在磁共振成像中，现有方法对图像的块稀疏和局部块相似性进行了研究，但是这些方法没有充分考虑到这些块的非局部相似性和相似块组成矩阵的低秩先验信息。非局部相似块组成的低秩矩阵在经过稀疏变换之后的非零系数位置不像一般稀疏的非零系数那样是随机分布的，而是趋向于聚类。现有技术的缺陷在于求解低秩矩阵恢复问题，是一个NP难问题，目前还没有一个很好的解决手段。一般方法是将最低秩问题用核范数最小替代(矩阵核范数的定义为矩阵的特征值之和)，松弛到核范数求解。因此业界需要一种更好的求解带有非局部先验和低秩先验信息的MRI模型的算法，以更好地描述解剖结构，重建更精确的磁共振图像。
发明内容
本发明的目的是提出一种基于非凸p范数低秩约束并迭代求解的磁共振超欠采样K数据成像方法(NLR_MRI)。
本发明基于非凸p范数低秩约束并迭代求解的磁共振超欠采样K数据成像方法中建立基于非局部相似图像块低秩先验信息的MR图像数据重构数学模型，利用傅里叶变换和低秩矩阵的特性简化计算过程，降低算法复杂度。并釆用交替方向法(ADMM)对模型进行迭代求解。相似图像块经过聚类变成低秩矩阵后，可以用更少的稀疏系数和更少的计算量来重构图像。本发明对非局部相似块的低秩特性进行深入研究，充分考虑MR图像的局部和非局部信息，提出基于低秩先MR图像稀疏先验，并建立基于图像低秩先验信息的低采样率MRI重构模型。仿真和实际MRI数据的实验表明本发明具有较强的细节保持能力，同时与其他方法的对比实验也表明所提出的方法能较好地重构，并且可以在更少的测量下更精确的重建图像，减少重建伪迹，恢复更多的图像细节。
本发明通过以下技术方案的步骤实现：
步骤(a)：初始化，设置迭代步数为k，循环迭代总次数为K，k取值为1到K的整数，k＝1。
步骤(b)：在稀疏表示和低秩约束的基础上结合图像块的非局部相似性和低秩特性，得到相似图像块组成的低秩矩阵，并基于非局部相似图像块的低秩先验信息构建磁共振图像稀疏表示模型。
步骤(c)：先固定MR图像，通过泰勒一阶近似和奇异值分解迭代求解基于带低秩先验信息的非局部图像稀疏表示模型的低秩矩阵的非凸p范数，得到相似图像块。
步骤(d)：固定低秩矩阵，通过增加辅助变量和分离变量法迭代求重建图像。当迭代步数k<K时，k＝k+1，转步骤(c)；当迭代步数k≥K时，停止迭代，得到重建图像。
进一步地说，本发明所述步骤(b)为：在稀疏表示和低秩约束的基础上用带非凸p范数相似图像块进行低秩约束，相对于凸1范数约束，将非凸p范数用于相似块的非凸低秩约束可以更接近于0范数，更好地得到非局部相似图像块的低秩先验模型(也就是带低秩先验信息的非局部图像稀疏表示模型)，因此进一步地稀疏表示图像，得到磁共振图像稀疏表示模型。
进一步地说，本发明所述步骤(c)为：先固定MR图像，求解步骤(b)中相似块的非凸p范数，并用局部最小化泰勒一阶近似和奇异值分解算法进行迭代求解，得到加权阈值处理后的相似图像块。
进一步地说，本发明所述步骤(b)为：在稀疏表示和低秩约束的基础上结合图像块的非局部相似性和低秩特性，融入基于非局部相似图像块的低秩先验信息，建立的磁共振图像稀疏表示模型为：
其中，y表示磁共振K空间欠采样数据，Φ表示磁共振K空间欠采样矩阵，表示磁共振重建图像，表示相似块低秩矩阵；模型中第一个子问题的后两项为图像进行组群稀疏和相似块低秩约束的正则项，第一项保证重建结果与K空间欠采样数据保持匹配约束，λ表示图像块的稀疏水平，η表示图像块相似程度拟合的权重。
进一步地说，本发明所述步骤(c)为：先固定MR图像，求解步骤(b)中相似块的非凸函数，并用局部最小化泰勒一阶近似和奇异值分解算法迭代求解带非凸p范数先验信息约束的目标函数。分离出与Li相关项如下：
对式中的第二项非凸p范数函数进行泰勒一阶近似，在(k+1)步迭代中的Li可以通过下式对奇异值分解进行加权阈值处理得到重建图像块：
其中表示Xi的奇异值分解。
进一步地说，本发明所述步骤(d)为：固定低秩矩阵Li，得到求解图像模型：
进一步地说，本发明所述步骤(d)引入辅助变量后，将图像模型变为：
其中是一个辅助参量，β是一个正标量，为拉格朗日乘子，η为图像块相似程度拟合权重。采用分离变量法，将图像模型变为：
μ(l+1)＝μ(l)+β(l)(x(l+1)-z(l+1))，
β(l+1)＝ρβ(l)，
轮换地更新一个变量，同时固定其它变量：固定辅助变量zl，求解xl+1；固定图像xl，通过最小化二次多项式更新zl+1。
进一步地说，本发明所述步骤(d)为对于图像x的子问题，图像x通过求解最小二乘解析问题得：
进一步地说，本发明所述步骤在稀疏表示和低秩约束的基础上融入带非凸p范数的低秩约束进行相似块间的稀疏约束，寻找相似块，进行相似块的组群稀疏和低秩约束，求解带非凸p范数低秩约束的目标函数，重建磁共振超欠采样K数据，得到最终成像结果。
本发明的技术方案具有以下的优点或有益效果：本发明实施例基于非凸低秩算法的磁共振超欠采样K数据成像方法中建立基于非局部相似图像块低秩先验信息的MR图像数据重构数学模型。本技术方案发明方法相对于标准核范数的低秩约束磁共振重建算法(LR_MRI)，使用非凸p范数约束可以进一步的利用MR图像的局部和非局部信息，进一步地稀疏表示图像，在更少的测量数据下更精确的重建图像，恢复更多的图像细节。
附图说明
图1为本发明算法步骤的流程图。
图2为同一采样率下同采样轨迹模板的NLR_MRI的重建结果的峰值信噪比(PSNR)值随P值的变化情况。
图3为P取p＝0.2,0.4,0.6NLR_MRI的重建情况，其中：(a)为原图，(b)(c)(d)分别为p＝0.2,0.4,0.6时NLR_MRI的重建结果图，(e)(f)(g)为p＝0.2,0.4,0.6时NLR_MRI的在的重建误差图。
图4为随机采样轨迹下DLMRI，LR_MRI和NLR_MRI三种算法重建图像的峰值信噪比(PSNR)值随采样因子(Sampling Factor)的变化情况。
图5为随机采样轨迹下DLMRI，LR_MRI和NLR_MRI三种算法重建图像的高频误差(HFEN)值随采样因子(Sampling Factor)的变化情况。
图6为采样率为0.21的随机采样轨迹下DLMRI，LR_MRI，NLR_MRI三种算法的重建性能分析情况。其中：(a)为原图，(b)(c)(d)分别为DLMRI，LR_MRI和NLR_MRI三种算法0.21倍欠采样率下重建结果图。(e)(f)(g)分别为DLMRI，LR_MRI和NLR_MRI三种算法0.21倍欠采样率下重建误差图。
图7为随机采样轨迹和径向采样轨迹下DLMRI，LR_MRI和NLR_MRI三种算法重建结果的峰值信噪比(PSNR)值随迭代次数的变化情况。
图8为随机采样轨迹和径向采样轨迹下DLMRI，LR_MRI和NLR_MRI三种算法重建结果的高频误差(HFEN)值随迭代次数的变化情况。
图9为采样率为0.21的随机采样轨迹下DLMRI，LR_MRI，NLR_MRI三种算法的重建性能分析情况。其中：(a)为原图，(b)(c)(d)分别为DLMRI，LR_MRI和NLR_MRI三种算法0.21倍欠采样率下重建结果图。(e)(f)(g)分别为DLMRI，LR_MRI和NLR_MRI三种算法0.21倍欠采样率下重建误差图。
图10为模拟径向采样轨迹下不同高斯白噪声标准偏差(Standard-Deviation)下DLMRI，LR_MRI和NLR_MRI三种重建算法的峰值信噪比(PSNR)值，
图11为模拟径向采样轨迹下不同高斯白噪声标准偏差(Standard-Deviation)下DLMRI，LR_MRI和NLR_MRI三种重建算法的高频误差(HFEN)值，
图12为模拟径向采样轨迹下不同复杂程度的高斯白噪声下DLMRI，LR_MRI和NLR_MRI三种算法的重建结果。其中：(a)为原图，(b)(c)(d)为DLMRI，LR_MRI和NLR_MRI三种算法在σ＝10时的重建结果图。(e)(f)(g)为DLMRI，LR_MRI和NLR_MRI三种算法在σ＝10时的重建误差图。
图13为磁共振系统评估物理模型的比较。其中：(a)为全采样重建图像，(b)为采样率为0.21的随机欠采样模板；(c)(d)(e)分别为DLMRI，LR_MRI和NLR_MRI三种算法的重建图像，其中NLR_MRI重建时p＝0.5；(f)为(a)(c)(d)(e)的放大图像。
图14为采样率为0.29的径向采样模板下三种重建算法重建性能分析。其中：(a)为全采样重建图像，(b)为采样率为0.29的径向采样模板；(c)(d)(e)分别为DLMRI，LR_MRI和NLR_MRI三种算法的重建图像，其中NLR_MRI重建时p＝0.5；(f)(g)(h)分别为DLMRI，LR_MRI和NLR_MRI三种算法的重建误差图像。
具体实施方式
为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施案例，对本发明进行进一步的详细说明。此处所描述的具体实施例仅用于解释本发明技术方案，并不限于本发明。
参见示出本发明的附图，下文将更详细地描述本发明。现参考附图1描述根据本发明的用于磁共振快速成像的非凸低秩重建算法。根据本发明的方法，本发明技术方案在基于MR图像块的非局部相似性和低秩特性，建立基于图像低秩先验的低采样率MRI重构模型，利用局部最小化泰勒一阶近似和奇异值分解迭代求解非凸P范数低秩约束函数，得到更好的图像重建效果。本发明实施例在基于非局部相似图像块低秩先验的MR图像数据重构数学模型上，釆用交替方向迭代方法对模型进行迭代求解，并利用傅里叶变换和低秩矩阵的特性简化计算过程，降低算法复杂度。相似图像块经过聚类变成低秩矩阵后，可以用更少的稀疏系数和更少的计算量来重构图像，并且可以在更少的测量K数据下更精确的重建图像，恢复更多的图像细节。本发明技术方案算法的重建结果可以避免混叠效应，获得更清晰的图像对比度以及更精确的解剖结构描述。
步骤S101：初始化，设置迭代步数为k，循环迭代总次数为K，k取值为1到K的整数，k＝1。
步骤S102：基于非局部相似图像块的低秩先验信息构建图像稀疏表示模型，结合图像块的非局部相似性和低秩特性，得到相似图像块组成的低秩矩阵。将低秩矩阵的求解用非凸p范数近似求解。建立磁共振图像稀疏表示模型：
y表示磁共振K空间欠采样数据，Φ表示磁共振K空间欠采样矩阵，表示磁共振重建图像，表示相似块低秩矩阵；后两项作为为图像进行组群稀疏和相似块低秩约束的正则项，如下式所示：
(1)式中第一项保证重建结果与K空间欠采样数据保持匹配约束，λ表示图像块的稀疏水平，η表示图像块相似程度拟合的权重。
在磁共振成像中，图像数据经过傅里叶变换到频率域得到一系列变换系数K数据，当使用凸优化l1范数近似l0范数时，会使重建的磁共振图像出现拖尾现象，产生混叠效应，重建图像解剖结构不清晰，对比度不明显。
假设在磁共振图像中每一块样本块在位置i处，大小为有很多的相似块。在此假设下，设置一个门限，对每个样本块进行k邻域搜索
式中T是预定义的临界值，Gi表示相似块位置的集合。再把相似块放到一起组成一个数据矩阵如上假设下，得到的数据矩阵Xi具有低秩的性质。事实上，Xi会被一些噪声腐蚀。为了更好地重建图像，将Xi分解成两个部分，即Xi＝Li+Wi，其中Li和Wi分别是低秩矩阵和高斯噪声。因此低秩问题可以表示为以下式子：
其中表示Fobenius范数，表示加性高斯噪声的方差。
本发明中通过求解以下局部最小化泰勒一阶近似和奇异值分解求解低秩问题
E(X，ε)＝lp det(X+εI) (5)其中lp det(X)＝λXp，(X≥0)，X表示非局部相似图像块组成的矩阵，I表示单位矩阵，λ表示图像块的稀疏水平，ε是一个很小的参量。而对于一般的矩阵把n≤m，把代入式(5)
其中，U表示奇异值分解得到的特征矩阵，∑是的特征值矩阵，即no＝min(n，m)，σj(Li)表示Li的第j个奇异值，并且∑1/2是一个对角矩阵，它对角上的元素是矩阵Li的奇异值。把式(6)代入到式(4)中可以得到下式：
将上式带约束不等式转换成下面无约束等式：
(8)步骤S103：先固定MR图像，求解步骤(102)中相似块的非凸函数，并用局部最小化泰勒一阶近似和奇异值分解算法求解带非凸p范数先验信息约束的目标函数。对于计算Li子问题，可以通过解下面最小化函数得到
(9)式中，η表示图像块相似程度拟合的权重，表示对图像形成非局部相似图像块的算子，x表示要重建的磁共振图像；
将式(6)代入到式(9)中，得到
其中定义是一个非凸函数，用局部最小化一阶泰勒级数方法求解。将f(σ)展开得
其中σ(k)表示在第k步迭代中σ的值。因此式(11)可以转换成下式
其中p＜1是lp函数中的参数，并忽略式(11)中的常数部分。为简化式子，将式(12)重新写成下式
其中τ＝λp/2η，
表示权重函数。
根据权重核范数的邻近算法，在(k+1)步迭代中的Li可以通过下式得到：
其中表示Xi的奇异值分解，(x)+＝max{x，0}。虽然它不是全局最小解，但它在局部总是使目标函数值单调递减。
步骤S104：固定低秩矩阵,通过最小二乘法迭代求重建图像，进行图像更新，得到重建图像。在求解得到Li后：
当矩阵Φ是局部傅里叶变换矩阵时，式(15)用增加辅助变量和分离变量法快速地求解出来。将式(15)分解成闭形解的形式：
其中是一个辅助参量，β是一个正标量，为拉格朗日乘子，η为图像块相似程度拟合权重。式(16)的最优解可以通过下式得到：
其中ρ＞1是一个常数。上标表示迭代序号；先固定x(l)，μ(l)和β(l)的值，再计算z(l+1)的近似解
注意到是一个对角矩阵，对角元素上的值是重叠块的个数，并且是相似块的平均值。接下来固定z(l+1)，μ(l)和β(l)的值来计算x子问题：
把Φ＝DF代入到式(19)中并且对等式两边同时取傅里叶变换，得到下面式子：
其中D是下采样矩阵，F表示傅里叶变换矩阵。对式(20)进行简化得
注意到(DTD+β(l))是一个很容易求解的对角矩阵。然后对上式进行傅里叶反变换，得
最后，当迭代步数k<K时，k＝k+1，转步骤103更新低秩矩阵；当迭代步数k≥K时，停止迭代，得到重建图像。
综上所述在稀疏表示和低秩约束的基础上融入带非凸p范数的低秩约束进行相似块间的稀疏约束及相似块的组群稀疏和低秩约束，近似求解带非凸p范数的目标函数，重建磁共振超欠采样K数据，得到最终成像结果。
具体而言，本发明实施例基于非凸p范数低秩约束并迭代求解的磁共振超欠采样K数据成像方法中建立基于非局部相似图像块低秩先验信息的MR图像数据重构数学模型，釆用交替方向乘子迭代方法(ADMM)对模型进行迭代求解再同时利用傅里叶变换和低秩矩阵的特性简化计算过程，降低算法复杂度。本发明对非局部相似块的低秩特性进行深入研究，充分考虑MR图像的局部和非局部信息，建立基于图像低秩先验信息的低采样率MRI重构模型。本发明具有较强的细节保持能力，同时对比与其他方法，实验也表明所提出的方法可以在更少的测量下更精确的重建图像，减少重建图像的伪影，恢复更多的图像细节。
综上所述，本发明的实施例提出完整的NLR_MRI算法可归纳如下：
(1)：初始化:ωi＝[1，1，...，1]T，x(1)，μ(1)＝0，λ，η，p，τ＝λp/2η,β,K,J,L；
(2)：当k＝1,2,...,K，进行(3)-(11)；
(3)：形成由相似块x(k)组成的矩阵Xi和集合
(4)：当j＝1,2,...,J，进行(5)-(6)；
(5)：若(k＞K0)，更新权重
(6)：通过式(14)计算Li，当j＝J时输出
(7)：当l＝1,2,...,L，进行(8)-(10)；
(8)：通过式(17)和式(22)分别计算z(l+1)和x(l+1)；
(9)：相应的μ(l+1)＝μ(l)+β(l)(x(l+1)-z(l+1))，β(l+1)＝ρβ(l)；
(10)：当l＝L时输出x(k)＝x(l+1)；
(11)：输出
本发明技术方案采用各种不同的欠采样因子和噪声对所提出方法的性能进行评估。本发明提出的NLR_MRI方法与DLMRI和LR_MRI相比，DLMRI直接通过正交匹配追踪方法解决l0最小化问题，LR_MRI则致力于l1解决低秩约束稀疏最小化问题，NLR_MRI方法则通过利用迭代求解带非凸p范数低秩先验信息约束的目标函数。实验过程用实数图像和复数图像来评估本发明提出的方法的性能。各种参数的标准值分别设置如下：图像的大小均为256×256，图像块大小为7×7，对每个样本块寻找50个相似块组成数据矩阵。重建图像的质量通过使用峰值信噪比(PSNR)和高频误差(HFEN)来衡量。
图1为示出本发明算法步骤的流程图。
图2为同一采样率下同采样轨迹模板的NLR_MRI的重建结果的峰值信噪比(PSNR)值随P值的变化情况。
图3为P取p＝0.2,0.4,0.6NLR_MRI的重建情况，其中：(a)为原图，(b)(c)(d)分别为p＝0.2,0.4,0.6时NLR_MRI的重建结果图，(e)(f)(g)为p＝0.2,0.4,0.6时NLR_MRI的在的重建误差图。
图4为随机采样轨迹下DLMRI，LR_MRI和NLR_MRI三种算法重建图像的峰值信噪比(PSNR)值随采样因子(Sampling Factor)的变化情况。从中可以看出LR_MRI与NLR_MRI两种算法重建的PSNR值之间的差距，NLR_MRI算法有较好的重建性能，可以更好的重建出图像细节。
图5为随机采样轨迹下DLMRI，LR_MRI和NLR_MRI三种算法重建图像的高频误差(HFEN)值随采样因子(Sampling Factor)的变化情况。特别地，当采样因子减小到0.1以下，非凸优化的优势降低，三种算法都不能从较少的K数据采样中较好的重建图像。
图6为采样率为0.21的随机采样轨迹下DLMRI，LR_MRI，NLR_MRI三种算法的重建性能分析情况。其中：(a)为原图，(b)(c)(d)分别为DLMRI，LR_MRI和NLR_MRI三种算法0.21倍欠采样率下重建结果图。(e)(f)(g)分别为DLMRI，LR_MRI和NLR_MRI三种算法0.21倍欠采样率下重建误差图。
图7为随机采样轨迹和径向采样轨迹下DLMRI，LR_MRI和NLR_MRI三种算法重建结果的峰值信噪比(PSNR)值随迭代次数的变化情况。
图8为随机采样轨迹和径向采样轨迹下DLMRI，LR_MRI和NLR_MRI三种算法重建结果的高频误差(HFEN)值随迭代次数的变化情况。
图9为采样率为0.21的随机采样轨迹下DLMRI，LR_MRI，NLR_MRI三种算法的重建性能分析情况。其中：(a)为原图，(b)(c)(d)分别为DLMRI，LR_MRI和NLR_MRI三种算法0.21倍欠采样率下重建结果图。(e)(f)(g)分别为DLMRI，LR_MRI和NLR_MRI三种算法0.21倍欠采样率下重建误差图。
图10为模拟径向采样轨迹下不同高斯白噪声标准偏差(Standard-Deviation)下DLMRI，LR_MRI和NLR_MRI三种重建算法的峰值信噪比(PSNR)值。
图11为模拟径向采样轨迹下不同高斯白噪声标准偏差(Standard-Deviation)下DLMRI，LR_MRI和NLR_MRI三种重建算法的高频误差(HFEN)值。
图12为模拟径向采样轨迹下不同复杂程度的高斯白噪声下DLMRI，LR_MRI和NLR_MRI三种算法的重建结果。其中：(a)为原图，(b)(c)(d)为DLMRI，LR_MRI和NLR_MRI三种算法在σ＝10时的重建结果图。(e)(f)(g)为DLMRI，LR_MRI和NLR_MRI三种算法在σ＝10时的重建误差图。从三种重建图像可以看出DLMRI算法比其他两种方法有更严重的重建伪影，本实施例NLR_MRI方法有更好的重建结果。
图13为磁共振系统评估物理模型的比较。其中：(a)为全采样重建图像，(b)为采样率为0.21的随机欠采样模板；(c)(d)(e)分别为DLMRI，LR_MRI和NLR_MRI三种算法的重建图像，其中NLR_MRI重建时p＝0.5；(f)为(a)(c)(d)(e)的放大图像。可以看出NLR_MRI算法可以更好的消除混叠伪影。
图14为在采样率为0.29的径向采样轨迹下DLMRI，LR_MRI和NLR_MRI三种算法的重建情况。其中：(a)为全采样重建图像，(b)为采样率为0.29的径向采样模板；(c)(d)(e)分别为DLMRI，LR_MRI和NLR_MRI三种算法的重建图像，其中NLR_MRI重建时p＝0.5；(f)(g)(h)分别为DLMRI，LR_MRI和NLR_MRI三种算法的重建误差图像。可以看出NLR_MRI在径向采样轨迹下的重建结果比DLMRI呈现更少的混叠伪影，恢复的图像有更好的保真度。
发明实施例基于非凸p范数低秩约束并迭代求解的磁共振超欠采样K数据成像方法在基于MR图像块的非局部相似性和低秩特性，建立基于图像低秩先验的低采样率MR图像重构模型，得到更好的图像重建效果。本发明实施例在基于非局部相似图像块低秩先验的MR图像数据重构数学模型，釆用交替方向迭代方法对模型进行迭代求解，并利用傅里叶变换和低秩矩阵的特性简化计算过程，降低算法复杂度。本发明实施例中相似图像块经过聚类变成低秩矩阵后，可以用更少的稀疏系数和更少的计算量来重构图像，因此进一步地稀疏表示图像，并更好地利用MR图像的局部信息和非局部相似性，在更少的测量数据下得到更精确的重建图像，减少重建图像的伪影，恢复更多的图像细节。本发明实施例算法重建结果可以避免混叠效应，获得更清晰的图像对比度以及更精确的解剖结构描述，恢复的图像有更好的保真度。
因本技术领域的技术人员应理解，本发明可以以许多其他具体形式实现而不脱离本发明的精神或范围。尽管业已描述了本发明的实施例，应理解本发明不应限制为这些实施例，本技术领域的技术人员可如所附权利要求书界定的本发明精神和范围之内作出变化和修改。
Claims (3)
1.一种用于磁共振快速成像的非凸低秩重建方法，其特征是按以下步骤：
步骤(a)：初始化，设置迭代步数为k，循环迭代总次数为K，k取值为1到K的整数，k＝1；
步骤(b)：在稀疏表示和低秩约束的基础上结合图像块的非局部相似性和低秩特性，用带非凸p范数相似图像块进行低秩约束，相对于凸1范数约束，将非凸p范数用于相似块的非凸低秩约束更接近于0范数，得到非局部相似图像块的低秩先验模型，得到相似图像块组成的低秩矩阵，并基于非局部相似图像块的低秩先验信息构建磁共振图像稀疏表示模型；
步骤(c)：先固定MR图像，通过泰勒一阶近似和奇异值分解迭代求解基于非局部相似图像块的低秩先验模型的低秩矩阵的非凸p范数，并用局部最小化泰勒一阶近似和奇异值分解算法进行迭代求解，得到加权阈值处理后的相似图像块；
步骤(d)：固定低秩矩阵，通过增加辅助变量和分离变量法迭代求重建图像；当迭代步数k<K时，k＝k+1，转步骤(c)；当迭代步数k≥K时，停止迭代，得到重建图像；
所述步骤(c)为：先固定MR图像，求解步骤(b)中相似块的非凸函数，并用局部最小化泰勒一阶近似和奇异值分解算法迭代求解带非凸p范数先验信息约束的目标函数：
E(X,ε)＝Ipdet(X+εI) (5)
其中，X表示非局部相似图像块组成的矩阵，I表示单位矩阵，Ipdet(X)＝λXp，(X≥0)，λ表示图像块的稀疏水平，ε是一个参量；而对于矩阵Li∈Cn×m，n≤m，把代入式(5)，得到
(6)式中，U表示奇异值分解得到的特征矩阵，∑是的特征值矩阵，即n0＝min(n,m)，σj表示奇异值，σj(Li)表示Li的第j个奇异值，并且∑1/2是一个对角矩阵，它对角上的元素是矩阵Li的奇异值；
分离出与Li相关项如下：
(9)式中，η表示图像块相似程度拟合的权重，表示对图像形成非局部相似图像块的算子，x表示要重建的磁共振图像；
将式(6)代入到式(9)中，得到：
其中，
定义是一个非凸函数，用一阶泰勒级数方法求解，将f(σ)展开得：
其中σ(k)表示在第k步迭代中σ的值；
对(9)式中的第二项非凸p范数函数进行泰勒一阶近似，在(k+1)步迭代中的Li可以通过下式对奇异值分解进行加权阈值处理得到重建图像块：
其中p<1是lp函数中的参数，并忽略式中的常数部分；
根据权重核范数的邻近算法，在(k+1)步迭代中的Li可以通过下式得到：
其中表示Xi的奇异值分解，τ＝λp/2η，
表示权重函数；
2.根据权利要求1所述的一种用于磁共振快速成像的非凸低秩重建方法，其特征是所述步骤(b)为：在稀疏表示和低秩约束的基础上结合图像块的非局部相似性和低秩特性，融入基于非局部相似图像块的低秩先验信息，建立的磁共振图像稀疏表示模型为：
其中，y表示磁共振K空间欠采样数据，Φ表示磁共振K空间欠采样矩阵，
表示磁共振重建图像，表示相似块低秩矩阵；模型中第一个子问题的后两项为图像进行组群稀疏和相似块低秩约束的正则项，第一项保证重建结果与K空间欠采样数据保持匹配约束，λ表示图像块的稀疏水平，η表示图像块相似程度拟合的权重。
3.根据权利要求1所述的一种用于磁共振快速成像的非凸低秩重建方法，其特征为：固定低秩矩阵Li，得到求解图像模型：
引入辅助变量后，将图像模型变为：
其中z∈CN是一个辅助参量，β是一个正标量，μ∈CN为拉格朗日乘子，η为图像块相似程度拟合权重；
采用分离变量法，将图像模型变为：
μ(l+1)＝μ(l)+β(l)(x(l+1)-z(l+1)),
β(l+1)＝ρβ(l),
ρ表示常数，上标表示迭代序号；轮换地更新一个变量，同时固定其它变量：固定辅助变量zl，求解xl+1；固定图像x1，通过最小化二次多项式更新zl+1。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201510310636.0A CN104933683B (zh) | 2015-06-09 | 2015-06-09 | 一种用于磁共振快速成像的非凸低秩重建方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201510310636.0A CN104933683B (zh) | 2015-06-09 | 2015-06-09 | 一种用于磁共振快速成像的非凸低秩重建方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN104933683A CN104933683A (zh) | 2015-09-23 |
| CN104933683B true CN104933683B (zh) | 2018-12-04 |
Family
ID=54120838
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201510310636.0A Expired - Fee Related CN104933683B (zh) | 2015-06-09 | 2015-06-09 | 一种用于磁共振快速成像的非凸低秩重建方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN104933683B (zh) |
Families Citing this family (31)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN105551000B (zh) * | 2015-12-11 | 2019-09-20 | 中国科学院遥感与数字地球研究所 | 基于参考影像约束和非凸低秩约束的遥感图像重建方法 |
| CN106023274B (zh) * | 2016-01-16 | 2019-03-15 | 中国科学院遥感与数字地球研究所 | 一种结合专家场滤波器稀疏约束的压缩感知图像重建方法 |
| CN105741240A (zh) * | 2016-01-18 | 2016-07-06 | 中国科学院遥感与数字地球研究所 | 基于参考影像纹理约束和非凸低秩约束的遥感图像重建方法 |
| CN105957117B (zh) * | 2016-04-26 | 2018-09-14 | 东南大学 | 并行磁共振的图像重建方法、装置及并行磁共振成像系统 |
| CN106093814A (zh) * | 2016-06-02 | 2016-11-09 | 浙江理工大学 | 一种基于多尺度低秩模型的心脏磁共振成像方法 |
| CN106780372B (zh) * | 2016-11-30 | 2019-06-18 | 华南理工大学 | 一种基于广义树稀疏的权重核范数磁共振成像重建方法 |
| CN106934778B (zh) * | 2017-03-10 | 2019-11-29 | 北京工业大学 | 一种基于小波域结构和非局部分组稀疏的mr图像重建方法 |
| CN107301631B (zh) * | 2017-06-28 | 2020-09-18 | 重庆大学 | 一种基于非凸加权稀疏约束的sar图像降斑方法 |
| CN107330953B (zh) * | 2017-07-06 | 2020-06-26 | 桂林电子科技大学 | 一种基于非凸低秩的动态mri重建方法 |
| CN107993205A (zh) * | 2017-11-28 | 2018-05-04 | 重庆大学 | 一种基于学习字典与非凸范数最小化约束的mri图像重构方法 |
| CN108577840A (zh) * | 2018-02-11 | 2018-09-28 | 南京邮电大学 | 一种动态磁共振稳健pca成像方法 |
| GB201802942D0 (en) * | 2018-02-23 | 2018-04-11 | Univ Leuven Kath | Reconstruction method |
| CN108872292B (zh) * | 2018-05-25 | 2019-11-08 | 中国石油大学(北京) | 基于稀疏采样重建的多维拉普拉斯磁共振方法与装置 |
| CN109145738B (zh) * | 2018-07-18 | 2021-07-27 | 浙江工业大学 | 基于加权非凸正则化和迭代重约束低秩表示的动态视频分割方法 |
| CN109215025B (zh) * | 2018-09-25 | 2021-08-10 | 电子科技大学 | 一种基于非凸秩逼近极小化的红外弱小目标检测方法 |
| CN111090068B (zh) * | 2018-10-23 | 2021-03-23 | 中国科学院深圳先进技术研究院 | 磁共振正对比成像方法及装置 |
| CN109920017B (zh) * | 2019-01-16 | 2022-06-21 | 昆明理工大学 | 基于特征向量的自一致性的联合全变分Lp伪范数的并行磁共振成像重构方法 |
| CN111818345B (zh) * | 2019-04-11 | 2021-10-26 | 四川大学 | 一种联合先验模型和细节增强的图像去压缩方法 |
| CN110505487B (zh) * | 2019-07-16 | 2020-11-10 | 珠海城市职业技术学院 | 基于非凸低秩模型图像量化去块方法、系统及存储介质 |
| CN110675331B (zh) * | 2019-08-13 | 2023-01-31 | 南京人工智能高等研究院有限公司 | 图像去噪方法、装置、计算机可读存储介质及电子设备 |
| CN110728624B (zh) * | 2019-09-29 | 2021-07-23 | 厦门大学 | 一种高分辨率扩散加权图像重建方法 |
| CN111260571B (zh) * | 2020-01-11 | 2022-09-23 | 浙江工业大学 | 一种基于非凸低秩低梯度的深度图像修复方法 |
| CN111314000B (zh) * | 2020-02-18 | 2022-04-15 | 五邑大学 | 一种基于低秩数据驱动权重模型的射频层析成像方法 |
| CN111714124B (zh) * | 2020-06-18 | 2023-11-03 | 中国科学院深圳先进技术研究院 | 磁共振电影成像方法、装置、成像设备及存储介质 |
| CN112710975A (zh) * | 2021-01-25 | 2021-04-27 | 东北林业大学 | 基于稀疏和局部低秩矩阵分解的磁共振扩散图像重建方法 |
| CN112819949B (zh) * | 2021-02-07 | 2024-03-26 | 哈尔滨工业大学 | 一种基于结构化低秩矩阵的磁共振指纹图像重建方法 |
| CN113034635B (zh) * | 2021-03-04 | 2022-08-23 | 重庆不贰科技(集团)有限公司 | 一种抑制工业成像阶梯伪影的图像重建方法 |
| CN113009560B (zh) * | 2021-03-23 | 2022-03-29 | 中国地质大学（武汉） | 一种地震数据重建方法、装置、设备及存储介质 |
| CN113628298B (zh) * | 2021-08-10 | 2023-09-08 | 昆明理工大学 | 基于特征向量的自一致性和非局部低秩的并行mri重构方法 |
| CN113892938B (zh) * | 2021-10-09 | 2023-10-27 | 昆明理工大学 | 一种基于非局部低秩约束的改进灵敏度编码重构方法 |
| CN113971706B (zh) * | 2021-10-15 | 2024-04-30 | 厦门大学 | 一种快速磁共振智能成像方法 |
Citations (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN104156994A (zh) * | 2014-08-14 | 2014-11-19 | 厦门大学 | 一种压缩感知磁共振成像的重建方法 |
Family Cites Families (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US8989465B2 (en) * | 2012-01-17 | 2015-03-24 | Mayo Foundation For Medical Education And Research | System and method for medical image reconstruction and image series denoising using local low rank promotion |
| WO2014162300A1 (en) * | 2013-04-05 | 2014-10-09 | Isis Innovation Ltd. | Acceleration of low-rank mri data acquisition |
-
2015
- 2015-06-09 CN CN201510310636.0A patent/CN104933683B/zh not_active Expired - Fee Related
Patent Citations (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN104156994A (zh) * | 2014-08-14 | 2014-11-19 | 厦门大学 | 一种压缩感知磁共振成像的重建方法 |
Non-Patent Citations (1)
| Title |
|---|
| Compressive Sensing via Nonlocal Low-Rank Regularization;Weisheng Dong 等;《IEEE TRANSACTIONS ON IMAGE PROCESSING》;20140606;第23卷(第8期);摘要，第3619页左栏第-2段、右栏第1段、右栏第3-4段，第3620页左栏第1-3段、右栏第2段，第3620-3624页第IV节；第3621页表1 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN104933683A (zh) | 2015-09-23 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN104933683B (zh) | 一种用于磁共振快速成像的非凸低秩重建方法 | |
| Schlemper et al. | A deep cascade of convolutional neural networks for dynamic MR image reconstruction | |
| Armanious et al. | Unsupervised medical image translation using cycle-MedGAN | |
| Wang et al. | Deep learning for fast MR imaging: A review for learning reconstruction from incomplete k-space data | |
| Ravishankar et al. | Data-driven learning of a union of sparsifying transforms model for blind compressed sensing | |
| WO2018099321A1 (zh) | 一种基于广义树稀疏的权重核范数磁共振成像重建方法 | |
| CN104156994B (zh) | 一种压缩感知磁共振成像的重建方法 | |
| CN106934778B (zh) | 一种基于小波域结构和非局部分组稀疏的mr图像重建方法 | |
| Kelkar et al. | Compressible latent-space invertible networks for generative model-constrained image reconstruction | |
| Kelkar et al. | Prior image-constrained reconstruction using style-based generative models | |
| Nguyen-Duc et al. | Frequency-splitting dynamic MRI reconstruction using multi-scale 3D convolutional sparse coding and automatic parameter selection | |
| CN112991483B (zh) | 一种非局部低秩约束的自校准并行磁共振成像重构方法 | |
| CN111754598A (zh) | 基于变换学习的局部空间邻域并行磁共振成像重构方法 | |
| Manimala et al. | Convolutional neural network for sparse reconstruction of MR images interposed with gaussian noise | |
| Cammarasana et al. | A universal deep learning framework for real-time denoising of ultrasound images | |
| CN109934884B (zh) | 一种基于变换学习和联合稀疏性的迭代自一致性并行成像重构方法 | |
| Manimala et al. | Sparse MR image reconstruction considering Rician noise models: A CNN approach | |
| Liu et al. | High-fidelity MRI reconstruction using adaptive spatial attention selection and deep data consistency prior | |
| CN108460810A (zh) | 一种全变分的并行磁共振图像快速重构方法 | |
| Liu et al. | Sparse and dense hybrid representation via subspace modeling for dynamic MRI | |
| CN114004764B (zh) | 一种基于稀疏变换学习的改进灵敏度编码重建方法 | |
| Golbabaee et al. | Compressive MRI quantification using convex spatiotemporal priors and deep auto-encoders | |
| Wang et al. | LANTERN: Learn analysis transform network for dynamic magnetic resonance imaging | |
| Heckel et al. | Deep learning for accelerated and robust MRI reconstruction: a review | |
| Luo et al. | Towards High-Quality MRI Reconstruction With Anisotropic Diffusion-Assisted Generative Adversarial Networks And Its Multi-Modal Images Extension |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| C06 | Publication | ||
| PB01 | Publication | ||
| C10 | Entry into substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| CB03 | Change of inventor or designer information |
Inventor after: Liu Qiegen Inventor after: Lu Hongyang Inventor after: Wu Xinfeng Inventor after: Long Chengzhi Inventor after: Wang Yuhao Inventor before: Lu Hongyang Inventor before: Liu Qiegen Inventor before: Wu Xinfeng Inventor before: Long Chengzhi Inventor before: Wang Yuhao |
|
| CB03 | Change of inventor or designer information | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant | ||
| TR01 | Transfer of patent right |
Effective date of registration: 20230316 Address after: 330006 Wannianqing Science Park, No. 399, Jingdong Avenue, Nanchang High tech Industrial Development Zone, Nanchang City, Jiangxi Province Patentee after: Nanchang Yunying Medical Technology Co.,Ltd. Address before: 999 No. 330031 Jiangxi province Nanchang Honggutan University Avenue Patentee before: Nanchang University |
|
| TR01 | Transfer of patent right | ||
| CF01 | Termination of patent right due to non-payment of annual fee |
Granted publication date: 20181204 |
|
| CF01 | Termination of patent right due to non-payment of annual fee |