CN111754598B - 基于变换学习的局部空间邻域并行磁共振成像重构方法 - Google Patents

基于变换学习的局部空间邻域并行磁共振成像重构方法 Download PDF

Info

Publication number
CN111754598B
CN111754598B CN202010593729.XA CN202010593729A CN111754598B CN 111754598 B CN111754598 B CN 111754598B CN 202010593729 A CN202010593729 A CN 202010593729A CN 111754598 B CN111754598 B CN 111754598B
Authority
CN
China
Prior art keywords
representing
iteration
matrix
image
auxiliary variable
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.)
Active
Application number
CN202010593729.XA
Other languages
English (en)
Other versions
CN111754598A (zh
Inventor
段继忠
刘畅
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Kunming University of Science and Technology
Original Assignee
Kunming University of Science and Technology
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by Kunming University of Science and Technology filed Critical Kunming University of Science and Technology
Priority to CN202010593729.XA priority Critical patent/CN111754598B/zh
Publication of CN111754598A publication Critical patent/CN111754598A/zh
Application granted granted Critical
Publication of CN111754598B publication Critical patent/CN111754598B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T11/00—2D [Two Dimensional] image generation
    • G06T11/003—Reconstruction from projections, e.g. tomography
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2210/00—Indexing scheme for image generation or computer graphics
    • G06T2210/41—Medical

Landscapes

  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • Theoretical Computer Science (AREA)
  • Magnetic Resonance Imaging Apparatus (AREA)

Abstract

本发明公开一种基于变换学习的局部空间邻域并行磁共振成像重构方法，属于磁共振图像处理技术领域。并行磁共振成像重构一直是近几年的研究热点，它是基于接收阵列的空间灵敏度的一种加速磁共振成像方法局部空间邻域的低秩模型(low‑rank modeling of local k‑space neighborhoods，LORAKS)，是近几年提出的一种利用图像的线性位移不变性，将k‑space数据映射到高维矩阵中，使用低秩正则化来重构图像的并行磁共振成像重构的方法。本发明基于LORAKS重构框架，提出一种结合变换学习(Transform Learning,TL)的高质量磁共振成像重构算法，命名为TL‑LORAKS。利用ADMM(Alternating Direction Method Of Multipliers)技术和共轭梯度法进行求解，实验结果表明，提出的算法有着更好的重构质量，同时拥有更好的边缘和细节。

Description

基于变换学习的局部空间邻域并行磁共振成像重构方法
技术领域
本发明涉及一种基于变换学习的局部空间邻域并行磁共振成像重构方法，属于磁共振图像处理技术领域。
背景技术
磁共振成像(Magnetic Resonance Imaging，MRI)在医学诊断中应用广泛，然而，高磁场对人体的影响仍然是未知的，通常MRI的扫描时间需要15到90分钟，且噪声巨大，可能使患者的听力受到损伤，且对体内有磁金属或起搏器的特殊病人不适用。MRI需要大量的时间来对数据进行采样，直接导致了检查时间较长和成像的成本较高。由于受到一些物理上和非物理上的限制，大多数的MRI成像技术都需要获取大量的数据，这些数据中的每一个序列都是k空间中的一部分，由获取的k空间数据来恢复图像。传统的k空间数据的采样方式受到奈奎斯特采样准则(Nyquist criterion)的限制，违反奈奎斯特采样准则会导致重构的图像中出现伪影，然而，围绕这个准则进行采样，会花费巨大的时间成本。
提高磁共振(Magnetic Resonance，MR)成像的扫描速度是近些年研究的热点，如何能在低于奈奎斯特采样频率的情况下采样，且能恢复出高质量的MR图像受到了极大的关注。并行MR成像技术使用多个线圈来接收数据，每个线圈有着不同的空间灵敏度，以此达到缩短MR数据的获取时间多通道的数据获取往往比单通道的数据获取多出一些信息，这意味着，可以在不满足奈奎斯特采样频率的情况下重构MR图像。早期的MR成像利用线圈的不同的灵敏度来获取k-space数据，以此减少数据的获取量，如同时获取空间谐波(Simultaneous acquisition of spatial harmonics,SMASH),灵敏度编码(SENSitivityEncoding，SENSE)等。SENSE解决了SMASH中一些关于校准的问题：如线圈的配置，切片的几何形状等；而这些方法的重构策略受到线性编码的缘故，重构方法也必须是线性的。一般自校准部分并行采样(GeneRalized Autocalibarating Partially Parallel Acquisition，GRAPPA)和迭代自一致性并行成像重构(Iterative Self-consistentparallel imagingreconstruction，SPIRiT)方法是一类基于自动校准的重构方法，只需要利用少量中心全采样信号(ACS)进行校准，无需显式使用线圈灵敏度，也就避免了SMASH和SENSE中灵敏度估计困难的问题。这一方法最主要的限制是要获取足够的自校准信号通常是不可行的，且对于高加速因子这一方法是无效的。
由于上述方法的不足，Uecker等提出可估计多组线圈灵敏度的模型ESPIRiT，当视场(Field of View,FOV)受限时可以采用多组线圈灵敏度进行重构，从而提高重构质量。研究者还提出了基于自适应稀疏表示的重构方法：基于字典学习的MR成像重构算法和基于稀疏变换学习的重构算法。字典学习MR成像(Dictionary Learning MR imaging，DLMRI)在重构的同时进行字典学习，分别使用K-SVD(Singular value decomposition，SVD)和正交匹配追踪(Orthogonal Matching Pursuit OMP)进行字典学习和稀疏编码。然而，字典学习和稀疏编码均需要迭代求解，故复杂度非常高。随后出现了一种方法变换学习(TransformLearning MR imaging,TLMRI),该方法的重要步骤变换学习和稀疏编码均有解析解，避免了DLMRI的迭代求解，因此提高了效率，又保证了重构质量。
最近，一种新的模型局部空间邻域的低秩模型(Low-Rank modeling of local K-Space neighborhoods,LORAKS)进入了视野，基于空间支持有限的图像或相位缓慢变化的图像在k-space中的线性依赖关系，利用一种新颖的方法将低维矩阵映射到高维矩阵中来解决问题。该方法可以很好的移植到并行MR成像，命名为P-LORAKS，有着不错的重构质量，且足够灵活，能与多种正则项相结合，进一步提高重构质量。
虽然，TLMRI方法相比于其他方法能有效的提高图像的质量和处理速度，但该方法在多线圈MR成像中的图像质量和处理速度上还有一定的改善空间；同样，P-LORAKS算法的成像质量还有提高的空间，可以尝试与其它正则项结合来提高成像质量。
发明内容
本发明的目的在于克服现有技术的不足，提供一种基于变换学习的局部空间邻域并行磁共振成像重构方法。
本发明的目的是通过以下技术方案来实现的：基于变换学习的局部空间邻域并行磁共振成像重构方法，它包括以下步骤：
S0：初始化，令B1,W1,x1,f1,ux 1＝0,uB 1＝0,i＝1；
其中，
Figure BDA0002556749910000021
表示
Figure BDA0002556749910000022
图像块的稀疏变换矩阵，n表示二维图像块的总像素点数；
Figure BDA0002556749910000031
为辅助变量，
Figure BDA0002556749910000032
Bj表示辅助变量B的第j个分量，L表示线圈图像的总个数，
Figure BDA0002556749910000033
为第j个矩阵提取线性算子，表示从每个线圈提取
Figure BDA0002556749910000034
图像块再组合成
Figure BDA0002556749910000035
的矩阵，Q＝nh×ny表示单幅图像的像素点个数，nh和ny表示单幅图像的长和宽；
Figure BDA0002556749910000036
为辅助变量，F表示傅里叶变换，F-1表示傅里叶逆变换，
Figure BDA0002556749910000037
为第c个线圈图像的k空间数据，
Figure BDA0002556749910000038
和
Figure BDA0002556749910000039
表示与算法相关的拉格朗日乘子，
Figure BDA00025567499100000310
i是循环变量，表示数据的第i次迭代，上标1表示第1次迭代。
S1：计算Wi+1，计算公式如下：
Wi+1＝VUH
式中，Wi+1表示第i+1次迭代的图像块的方形稀疏变换，定义U∑VH由
Figure BDA00025567499100000311
奇异值分解得到，U，∑，V是奇异值分解因子，上标H表示矩阵的共轭转置，
Figure BDA00025567499100000312
表示第i次迭代的拉格朗日乘子；
Figure BDA00025567499100000313
表示第i次迭代的辅助变量xi中抽取LQ个图像块并组成n×LQ的矩阵。
S2：初始化j＝1；
S3：计算第i+1次迭代的辅助变量Bi+1，计算公式如下：
Figure BDA00025567499100000314
式中，Bj i+1表示第i+1次迭代的辅助变量Bi+1的第j个分量，与Wi+1Pj(xi)对应，
Figure BDA00025567499100000315
表示第i次迭代的拉格朗日乘子
Figure BDA00025567499100000316
的第j个分块，HJ(δ,θ)是硬阈值函数，
Figure BDA00025567499100000317
δ表示输入矩阵，θ表示阈值，α＞0和μ2＞0为正则化参数。
S4:判断是否完成B的所有图像块分块的更新，若是，则进入步骤S5；否则，对j进行加一操作后返回S3；
S5：初始化，c＝1；
S6：计算第c个线圈的辅助变量
Figure BDA0002556749910000041
计算公式如下：
Figure BDA0002556749910000042
其中，
Figure BDA0002556749910000043
表示第i+1次的第c个线圈的辅助变量，
Figure BDA0002556749910000044
表示第i+1次迭代的辅助变量Bi+1的第j个分块的
Figure BDA0002556749910000045
的第c列，
Figure BDA0002556749910000046
表示第i次迭代的拉格朗日乘子
Figure BDA0002556749910000047
的第c列，
Figure BDA0002556749910000048
表示第i次迭代的拉格朗日乘子
Figure BDA0002556749910000049
的第j个分块的第c列，
Figure BDA00025567499100000410
表示第i+1次迭代的待重构的第c个线圈图像的k空间数据，μ1为正则化参数，WHW＝I是一个单位矩阵，
Figure BDA00025567499100000411
为对角矩阵。
S7：计算待重构的第c个线圈图像的k空间数据
Figure BDA00025567499100000412
使用共轭梯度法求解下述问题得到
Figure BDA00025567499100000413
Figure BDA00025567499100000414
式中，A表示欠采样矩阵，
Figure BDA00025567499100000415
表示第i次迭代的拉格朗日乘子
Figure BDA00025567499100000416
的第c列，dc表示第c个线圈的欠采样数据，λ为正则化参数，Nr(·)是一个构造矩阵的算子，构造出的矩阵的列是近似零空间的标准正交基，
Figure BDA00025567499100000417
是一个构造高维矩阵的线性算子，能够将k空间数据f(nh,ny)构造成如下矩阵：
Figure BDA00025567499100000418
矩阵RS(f)是近似低秩的。其中，矩阵
Figure BDA00025567499100000419
分别表示为：
Figure BDA0002556749910000051
Figure BDA0002556749910000052
Figure BDA0002556749910000053
Figure BDA0002556749910000054
k＝1,...,K，m＝1,...,NR，
Figure BDA0002556749910000055
和
Figure BDA0002556749910000056
分别代表
Figure BDA0002556749910000057
的实部和虚部，上标～表示该数据是频域数据。
Figure BDA0002556749910000058
表示K个不同的k空间位置，(nh,ny)满足-Nh+R≤nh≤Nh-R，-Ny+R≤ny≤Ny-R，Nh和Ny表示k空间度量区域的正整数，R为傅里叶截断半径，当
Figure BDA0002556749910000059
时，
Figure BDA00025567499100000510
使用
Figure BDA00025567499100000511
表示集合
Figure BDA00025567499100000512
中不同元素的有序组合，p和q表示截断傅里叶半径内的图像横纵坐标，下表m表示第m个像素点，NR表示ΛR的基数。
S8：判断是否完成所有线圈，若是，进入步骤S9；否则，对c进行加一操作后返回S6；
S9：计算拉格朗日乘子ux i+1＝ux i+(F-1fi+1-xi+1)。
其中，ux i+1是第i+1次迭代的辅助变量xi+1对应的拉格朗日乘子，fi+1表示第i+1次迭代的待重构的所有线圈图像的k空间数据；
S10：计算拉格朗日乘子
Figure BDA00025567499100000513
其中，
Figure BDA00025567499100000514
表示第i+1次迭代的辅助变量Bi+1对应的拉格朗日乘子，
Figure BDA00025567499100000515
表示从第i+1次迭代的辅助变量xi+1中抽取LQ个图像块并组成n×LQ的矩阵。
S11：判断是否达到最大迭代次数iter；若不大于则对i做加一操作之后返回步骤S1，否则进入步骤S12；
S12：输出重构的k空间数据f，进行傅里叶反变换之后再使用平方和的平方根方法来形成幅度图像，计算公式如下：
Figure BDA0002556749910000061
本发明的有益效果是：基于LORAKS模型，本发明提出了一种基于变换学习的局部空间邻域并行磁共振成像重构方法。首先，使用ADMM算法将稀疏约束正则项进行简化，再采用OS技术将简化后的重构问题转化成一个变换更新问题和一个稀疏编码问题和一个去噪问题，其次，去噪问题应用ADMM算法进行求解。理论分析和成像实验表明，与传统的重构算法相比，提出的新的并行成像算法可以更有效地提高重构图像的质量。
附图说明
图1为本发明方法流程图；
图2为使用了八通道的头部线圈进行完全采样的来获取受试者的体内人脑数据(即data1)图；
图3为8通道的膝盖数据(即data2)图；
图4为加速比约为6的泊松亚采样示意图；
图5-7为数据集data1下，分别使用P-LORAKS，JTV-LORAKS，TL-LORAKS重建并行MR图像(具有6×加速度和24×24中心校准区)；
图8-10为数据集data2下，分别使用P-LORAKS，JTV-LORAKS，TL-LORAKS重建并行MR图像(具有6×加速度和24×24中心校准区)；
图11-13为数据集data1下，分别使用P-LORAKS，JTV-LORAKS，TL-LORAKS重建的误差图(具有6×加速度和24×24中心校准区)；
图14-16为数据集data2下，分别使用P-LORAKS，JTV-LORAKS，TL-LORAKS重建的误差图(具有6×加速度和24×24中心校准区)。
具体实施方式
下面结合附图和具体实施例进一步详细描述本发明的技术方案。
实施例1：本发明是基于P-LORAKS框架提出的一种基于变换学习的局部空间邻域并行磁共振成像重构方法。
在P-LORAKS模型中，当线圈图像在k空间中具有线性依赖关系或者缓慢变化的相位时，全采样的k空间数据可以映射到一个低秩矩阵，这就使得利用低秩正则化来重构欠采样数据或者噪声数据成为可能，用下面这个式子表示MRI的采样过程：
Af＝d (1)
其中，
Figure BDA0002556749910000071
为第c个线圈图像的k空间数据，
Figure BDA0002556749910000072
是欠采样矩阵，
Figure BDA0002556749910000073
表示欠采样数据。
S矩阵可通过对k空间数据f(nh,ny)选取K个不同的像素点进行线性运算的线性算子
Figure BDA0002556749910000074
来构造，具体形式如下：
Figure BDA0002556749910000075
Figure BDA0002556749910000076
矩阵RS(f)是近似低秩的。其中，矩阵
Figure BDA0002556749910000077
分别表示为：
Figure BDA0002556749910000078
Figure BDA0002556749910000079
Figure BDA00025567499100000710
Figure BDA00025567499100000711
k＝1,...,K，m＝1,...,NR，
Figure BDA00025567499100000712
和
Figure BDA00025567499100000713
分别代表
Figure BDA00025567499100000714
的实部和虚部。
Figure BDA00025567499100000715
表示K个不同的k空间位置，(nh,ny)满足-Nh+R≤nh≤Nh-R，-Ny+R≤ny≤Ny-R，Nh和Ny表示k空间度量区域的正整数，R为傅里叶截断半径，当
Figure BDA00025567499100000716
时，
Figure BDA00025567499100000717
使用
Figure BDA00025567499100000718
表示集合
Figure BDA00025567499100000719
中不同元素的有序组合，p和q表示截断傅里叶半径内的图像横纵坐标，下表m表示第m个像素点，NR表示ΛR的基数。
P-LORAKS重建问题可以用下面的公式表示：
Figure BDA0002556749910000081
Figure BDA0002556749910000082
是每个线圈全采样k空间数据fc的集合；
Figure BDA0002556749910000083
是每个线圈欠采样数据dc的集合。
其中，λ为正则化参数，Jr(·)是非凸正则化函数，它迫使矩阵参数的秩小于等于r，r为一个整数，假设
Figure BDA0002556749910000084
则：
Figure BDA0002556749910000085
其中σk是矩阵X的第k个奇异值，IQ-r是一个(Q-r)×(Q-r)的单位矩阵。T是X的最优低秩近似，G是一个零填充矩阵。
P-LORAKS模型为并行成像提供了一种新的重构方法。尽管P-LORAKS指出了结合了JTV正则项的LORAKS模型重构算法可有效改善重构质量，但仍有较大的改进余地。为了进一步提高重构质量，本文将联合变换学习(TransformLearning,TL)与LORAKS模型相结合，提出TL-LORAKS算法，得到如下最优化问题：
Figure BDA0002556749910000086
其中，
Figure BDA0002556749910000087
表示
Figure BDA0002556749910000088
图像块的稀疏变换矩阵；
Figure BDA0002556749910000089
为第j个矩阵提取线性算子，表示从每个线圈提取
Figure BDA00025567499100000810
图像块再组合成
Figure BDA00025567499100000811
的矩阵，α为正则化参数。
算法推导：
引入辅助变量
Figure BDA00025567499100000812
以及拉格朗日乘子
Figure BDA00025567499100000813
则使用交替方向乘子法(Alternating Direction Method Of Multipliers，ADMM)技术可将问题(10)转化成如下形式的求解：
Figure BDA0002556749910000091
引入辅助变量
Figure BDA0002556749910000092
将Bj组合成
Figure BDA0002556749910000093
即B＝[B1,...,Bj,...,BQ]，并引入拉格朗日乘子
Figure BDA0002556749910000094
则式(11)可转化成：
Figure BDA0002556749910000095
使用ADMM算法，问题(12)可分解为以下迭代子问题：
Figure BDA0002556749910000096
Figure BDA0002556749910000097
Figure BDA0002556749910000098
Figure BDA0002556749910000099
Figure BDA00025567499100000910
Figure BDA00025567499100000911
在(13)-(18)中，W表示图像的方形稀疏变换，Wi+1表示第i+1次迭代的图像块的方形稀疏变换，
Figure BDA00025567499100000912
为辅助变量，
Figure BDA00025567499100000913
Bj表示辅助变量B的第j个分量，L表示线圈图像的总个数，
Figure BDA00025567499100000914
为第j个矩阵提取线性算子，表示从每个线圈提取
Figure BDA00025567499100000915
图像块再组合成
Figure BDA00025567499100000916
的矩阵，Q＝nh×ny表示单幅图像的像素点个数，nh和ny表示单幅图像的长和宽，
Figure BDA00025567499100000917
为辅助变量，表示第c个线圈的辅助变量，F表示傅里叶变换，F-1表示傅里叶逆变换，
Figure BDA00025567499100000918
为第c个线圈图像的k空间数据，xcBi+1第i+1次迭代的辅助变量，Bj i+1表示第i+1次迭代的辅助变量Bi+1的第j个分量，与Wi+1Pj(xi)对应，
Figure BDA0002556749910000101
和
Figure BDA0002556749910000102
表示与算法相关的拉格朗日乘子，
Figure BDA0002556749910000103
表示第i次迭代的拉格朗日乘子，
Figure BDA0002556749910000104
表示第i次迭代的拉格朗日乘子
Figure BDA0002556749910000105
的第j个分块，μ2＞0为正则化参数。
Figure BDA0002556749910000106
表示第i+1次的第c个线圈的辅助变量，
Figure BDA0002556749910000107
表示第i+1次迭代的辅助变量Bi+1的第j个分块的
Figure BDA0002556749910000108
的第c列，
Figure BDA0002556749910000109
表示第i次迭代的拉格朗日乘子
Figure BDA00025567499100001010
的第c列，
Figure BDA00025567499100001011
表示第i+1次迭代的待重构的第c个线圈图像的k空间数据，μ1为大于0的正则化参数，WHW＝I是一个单位矩阵，上标H表示矩阵的共轭转置，A表示欠采样矩阵，
Figure BDA00025567499100001012
表示第i次迭代的拉格朗日乘子
Figure BDA00025567499100001013
的第c列，dc表示第c个线圈的欠采样数据，λ为正则化参数，
Figure BDA00025567499100001014
是一个构造高维矩阵的线性算子。
对上述几个子问题进行计算：
变换更新：令
Figure BDA00025567499100001015
表示第i次迭代的辅助变量xi中抽取LQ个图像块并组成n×LQ的矩阵。则问题(13)可以重写成：
Figure BDA00025567499100001016
对
Figure BDA00025567499100001017
进行奇异值分解得到U∑VH,U，∑，V是奇异值分解因子，则W＝VUH。
稀疏编码：问题(14)可以写成：
Figure BDA00025567499100001018
该问题可以使用硬阈值法求解得到，
Figure BDA00025567499100001019
μ2为正则化参数，HJ(δ,θ)的定义为：
Figure BDA00025567499100001020
δ表示输入矩阵，θ表示阈值。
令问题(15)的目标函数对x的导数为0，则可得：
Figure BDA0002556749910000111
WHW＝I是一个单位矩阵，
Figure BDA0002556749910000112
为对角矩阵。
问题(16)的求解：问题(16)的目标函数非凸，计算相对比较复杂，故引入一种近似方法，引入函数Jr(Rx(f))在点fi处的一个majorizer：
Figure BDA0002556749910000113
假设
Figure BDA0002556749910000114
那么可以将Nr(X)定义为生成一个q×(q-r)矩阵的算子，这个矩阵的列等于与X的广义奇异值分解中最小的q-r或零奇异值下相关的右奇异向量；换句话说，Nr(X)的列构成了X的q-r维近似零空间的标准正交基。
利用式(23)，问题(16)可转化成：
Figure BDA0002556749910000115
问题(24)使用共轭梯度法求解。综上所述，子问题(13)(14)(15)(16)均可求解。
具体流程如图1所示，其中步骤如下：
S0：初始化，令B1,W1,x1,f1,ux 1＝0,uB 1＝0,i＝1；
B1表示WP(x)对应辅助变量的初始值，W1表示图像块的方形稀疏变换第一次迭代的初始值，x1表示F-1f对应辅助变量第一次迭代的初始值，f1表示待重构图像k空间数据第一次迭代的初始值，ux 1为辅助变量x对应的拉格朗日乘子第一次迭代的初始值，uB 1为辅助变量B对应的拉格朗日乘子第一次迭代的初始值，i为迭代次数的初始值。
S1：计算Wi+1，计算公式如(19)
S2：初始化j＝1；
S3：计算Bi+1，计算公式如(20)
S4：判断是否完成B的所有图像块分块的更新，若是，则进入步骤S5；否则，对j进行加一操作后返回S3；
S5：初始化，c＝1；
S6：计算
Figure BDA0002556749910000121
计算公式如(22)；
S7：计算
Figure BDA0002556749910000122
公式如(24)；
S8：判断是否完成所有线圈，若是，进入步骤S9；否则，对c进行加一操作后返回S6；
S9：计算拉格朗日乘子ux i+1，公式如(17)；
S10：计算拉格朗日乘子
Figure BDA0002556749910000123
公式如(18)；
S11：判断是否达到最大迭代次数iter；若不大于则对i做加一操作之后返回步骤S1，否则进入步骤S12；
S12：输出重构图像f，输出重构的k空间数据f，进行傅里叶反变换之后再使用平方和的平方根方法来形成幅度图像，计算公式如下：
Figure BDA0002556749910000124
下面结合具体的实验对本发明做进一步的说明。
在下面的实验中，比较了本发明提出的TL-LORAKS算法(问题(10)),添加了联合全变分(Joint Total Variation,JTV)正则项的JTV-LORAKS和传统的P-LORAKS的算法性能。实验中的所有算法都是使用Matlab(MathWorks，Natick，MA)实现的。
为了比较各算法的重构性能，使用了1张人类脑部切片图和1张人类膝盖图来进行仿真实验，命名为data1,data2(如图2，图3)。data1为八通道的头部切片，data2为八通道的膝盖图。为了生成测试数据集，使用加速因子为R×的笛卡尔泊松亚采样，对全采样的数据集进行人工采样并用于仿真重建，如图4所示。在以下实验中，对于所有比较的算法，使用尺寸为24×24校准区域。
所有的实验都是在配置Intel(R)Core(TM)i5-10210U@2.6GHz处理器,12GB内存，Windows 10操作系统(64位)的笔记本电脑上进行的。比较算法的正则化参数是针对SNR最优进行调整的。在以下实验中，信噪比(Signal Noise Ratio，SNR)和均方根误差(RootMean Square Error，NRMSE)用于定量评估重建图像的质量，分别定义为：
Figure BDA0002556749910000131
MSE表示通过使用来自全采样数据集的平方和的平方根(SRSoS)方法计算的重建图像和参考图像之间的均方误差，Var表示参考图像的方差。
Figure BDA0002556749910000132
其中x是参考图像，是
Figure BDA0002556749910000133
重构图像。
首先，对加速因子为6时，四个数据集下的各个算法的重构图像进行了视觉上的比较，如图5-图10所示。图5、图6和图7展示了在数据集data1下，算法P-LORAKS,算法JTV-LORAKS和算法TL-LORAKS的重构图像。从图5可以看出，使用P-LORAKS算法的重构图像有着明显的模糊伪影。图6使用JTV-LORAKS算法的重构图像质量稍好于图5，但是也存在一些伪影，图7使用TL稀疏正则项的TL-LORAKS算法的重构图像质量明显优于使用P-LORAKS算法和JTV-LORAKS算法的重构质量。在数据集data2下,如图8-10所示，TL-LORAKS算法比JTV-LORAKS算法和P-LORAKS算法更为有效的去除重构图像中的alias伪影，并且保持了物体锐利的边缘，拥有更好的细节。
为了进一步比较提出的新算法的重构性能，在图11-图16中展示了两种算法在数据集data1和data2下的误差图像(白点越多，误差越大)。从图11，图14可以明显看出，P-LORAKS重构图像有较大的误差。图12，图15能有效地改善重构误差，不过还是存在一些误差。而图13，图16能更有效的改善重构误差，拥有更好重构质量。综上所述，在两个测试集中，使用TL稀疏正则项的TL-LORAKS算法总是能够获得质量最好的重构图像，而JTV-LORAKS重构图像质量较差，P-LORAKS的重构质量是三种算法中最差的。
本发明将自适应的变换学习与P-LORAKS模型结合，提出了一种新的并行MR成像方法，利用变换学习的方式增强图像的稀疏性，使用P-LORAKS模型来重构图像。进行了多组试验后，实验结果显示，两种方式结合的成像结果比分别使用两种方法的成像结果要好得多。
以上结合附图对本发明的具体实施方式作了详细说明，但是本发明并不限于上述实施方式，在本领域普通技术人员所具备的知识范围内，还可以在不脱离本发明宗旨的前提下作出各种变化。

Claims (1)

1.一种基于变换学习的局部空间邻域并行磁共振成像重构方法，其特征在于：包括以下步骤：
S0：初始化，令B1,W1,x1,f1,ux 1＝0,uB 1＝0,i＝1；
其中，
Figure FDA0002556749900000011
表示
Figure FDA0002556749900000012
图像块的稀疏变换矩阵，n表示二维图像块的总像素点数；
Figure FDA0002556749900000013
为辅助变量，
Figure FDA0002556749900000014
Bj表示辅助变量B的第j个分量，L表示线圈图像的总个数，
Figure FDA0002556749900000015
为第j个矩阵提取线性算子，表示从每个线圈提取
Figure FDA0002556749900000016
图像块再组合成
Figure FDA0002556749900000017
的矩阵，Q＝nh×ny表示单幅图像的像素点个数，nh和ny表示单幅图像的长和宽；
Figure FDA0002556749900000018
为辅助变量，F表示傅里叶变换，F-1表示傅里叶逆变换，f＝[f1,...,fc,..., fL]，
Figure FDA0002556749900000019
为第c个线圈图像的k空间数据，
Figure FDA00025567499000000110
和
Figure FDA00025567499000000111
表示与算法相关的拉格朗日乘子，
Figure FDA00025567499000000112
i是循环变量，表示数据的第i次迭代，上标1表示第1次迭代；
S1：计算Wi+1，计算公式如下：
Wi+1＝VUH
式中，Wi+1表示第i+1次迭代的图像块的方形稀疏变换，定义U∑VH由
Figure FDA00025567499000000113
奇异值分解得到，U，∑，V是奇异值分解因子，上标H表示矩阵的共轭转置，
Figure FDA00025567499000000114
表示第i次迭代的拉格朗日乘子；
Figure FDA00025567499000000115
表示第i次迭代的辅助变量xi中抽取LQ个图像块并组成n×LQ的矩阵；
S2：初始化j＝1；
S3：计算第i+1次迭代的辅助变量Bi+1，计算公式如下：
Figure FDA00025567499000000116
式中，Bj i+1表示第i+1次迭代的辅助变量Bi+1的第j个分量，与Wi+1Pj(xi)对应，
Figure FDA0002556749900000021
表示第i次迭代的拉格朗日乘子
Figure FDA0002556749900000022
的第j个分块，HJ(δ,θ)是硬阈值函数，
Figure FDA0002556749900000023
δ表示输入矩阵，θ表示阈值，α＞0和μ2＞0为正则化参数；
S4:判断是否完成B的所有图像块分块的更新，若是，则进入步骤S5；否则，对j进行加一操作后返回S3；
S5：初始化，c＝1；
S6：计算第c个线圈的辅助变量
Figure FDA0002556749900000024
计算公式如下：
Figure FDA0002556749900000025
其中，
Figure FDA0002556749900000026
表示第i+1次的第c个线圈的辅助变量，
Figure FDA0002556749900000027
表示第i+1次迭代的辅助变量Bi+1的第j个分块的
Figure FDA0002556749900000028
的第c列，
Figure FDA0002556749900000029
表示第i次迭代的拉格朗日乘子
Figure FDA00025567499000000210
的第c列，
Figure FDA00025567499000000211
表示第i次迭代的拉格朗日乘子
Figure FDA00025567499000000212
的第j个分块的第c列，
Figure FDA00025567499000000213
表示第i+1次迭代的待重构的第c个线圈图像的k空间数据，μ1为正则化参数，WHW＝I是一个单位矩阵，
Figure FDA00025567499000000214
为对角矩阵；
S7：计算待重构的第c个线圈图像的k空间数据
Figure FDA00025567499000000215
使用共轭梯度法求解下述问题得到
Figure FDA00025567499000000216
Figure FDA00025567499000000217
式中，A表示欠采样矩阵，
Figure FDA00025567499000000218
表示第i次迭代的拉格朗日乘子
Figure FDA00025567499000000219
的第c列，dc表示第c个线圈的欠采样数据，λ为正则化参数，Nr(·)是一个构造矩阵的算子，构造出的矩阵的列是近似零空间的标准正交基，
Figure FDA00025567499000000220
是一个构造高维矩阵的线性算子，能够将k空间数据f(nh,ny)构造成如下矩阵：
Figure FDA0002556749900000031
矩阵RS(f)是近似低秩的，其中，矩阵
Figure FDA0002556749900000032
分别表示为：
Figure FDA0002556749900000033
Figure FDA0002556749900000034
Figure FDA0002556749900000035
Figure FDA0002556749900000036
k＝1,...,K，m＝1,...,NR，
Figure FDA0002556749900000037
和
Figure FDA0002556749900000038
分别代表
Figure FDA0002556749900000039
的实部和虚部，上标～表示该数据是频域数据，
Figure FDA00025567499000000310
表示K个不同的k空间位置，(nh,ny)满足-Nh+R≤nh≤Nh-R，-Ny+R≤ny≤Ny-R，Nh和Ny表示k空间度量区域的正整数，R为傅里叶截断半径，当
Figure FDA00025567499000000311
时，
Figure FDA00025567499000000312
使用
Figure FDA00025567499000000313
表示集合
Figure FDA00025567499000000314
中不同元素的有序组合，p和q表示截断傅里叶半径内的图像横纵坐标，下表m表示第m个像素点，NR表示ΛR的基数；
S8：判断是否完成所有线圈，若是，进入步骤S9；否则，对c进行加一操作后返回S6；
S9：计算拉格朗日乘子ux i+1＝ux i+(F-1fi+1-xi+1)；
其中，ux i+1是第i+1次迭代的辅助变量xi+1对应的拉格朗日乘子，fi+1表示第i+1次迭代的待重构的所有线圈图像的k空间数据；
S10：计算拉格朗日乘子
Figure FDA00025567499000000315
其中，
Figure FDA00025567499000000316
表示第i+1次迭代的辅助变量Bi+1对应的拉格朗日乘子，
Figure FDA00025567499000000317
表示从第i+1次迭代的辅助变量xi+1中抽取LQ个图像块并组成n×LQ的矩阵；
S11：判断是否达到最大迭代次数iter；若不大于则对i做加一操作之后返回步骤S1，否则进入步骤S12；
S12：输出重构的k空间数据f，进行傅里叶反变换之后再使用平方和的平方根方法来形成幅度图像，计算公式如下：
Figure FDA0002556749900000041
CN202010593729.XA 2020-06-27 2020-06-27 基于变换学习的局部空间邻域并行磁共振成像重构方法 Active CN111754598B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010593729.XA CN111754598B (zh) 2020-06-27 2020-06-27 基于变换学习的局部空间邻域并行磁共振成像重构方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010593729.XA CN111754598B (zh) 2020-06-27 2020-06-27 基于变换学习的局部空间邻域并行磁共振成像重构方法

Publications (2)

Publication Number Publication Date
CN111754598A CN111754598A (zh) 2020-10-09
CN111754598B true CN111754598B (zh) 2022-02-25

Family

ID=72677368

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010593729.XA Active CN111754598B (zh) 2020-06-27 2020-06-27 基于变换学习的局部空间邻域并行磁共振成像重构方法

Country Status (1)

Country Link
CN (1) CN111754598B (zh)

Families Citing this family (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112991483B (zh) * 2021-04-26 2022-03-01 昆明理工大学 一种非局部低秩约束的自校准并行磁共振成像重构方法
CN113628298B (zh) * 2021-08-10 2023-09-08 昆明理工大学 基于特征向量的自一致性和非局部低秩的并行mri重构方法
CN114004764B (zh) * 2021-11-03 2024-03-15 昆明理工大学 一种基于稀疏变换学习的改进灵敏度编码重建方法

Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102930567A (zh) * 2012-09-25 2013-02-13 电子科技大学 多核加权最小二乘支撑向量机的磁共振并行成像重建方法
CN103793711A (zh) * 2014-01-17 2014-05-14 首都医科大学 一种基于脑部核磁共振图像的多维度纹理提取方法
CN104282036A (zh) * 2013-07-09 2015-01-14 韦伯斯特生物官能(以色列)有限公司 由稀疏样本对心脏进行的基于模型的重构
CN107730516A (zh) * 2017-10-16 2018-02-23 江南大学 一种基于模糊聚类的脑mr影像分割方法
CN107991636A (zh) * 2017-11-24 2018-05-04 哈尔滨工业大学 一种基于适应性结构低秩矩阵的快速磁共振图像重建方法
CN108198193A (zh) * 2018-01-16 2018-06-22 北京航空航天大学 一种利用改进直觉模糊聚类算法分割红外舰船图像的方法
CN109920017A (zh) * 2019-01-16 2019-06-21 昆明理工大学 基于特征向量的自一致性的联合全变分Lp伪范数的并行磁共振成像重构方法
CN109934884A (zh) * 2019-01-22 2019-06-25 昆明理工大学 一种基于变换学习和联合稀疏性的迭代自一致性并行成像重构方法

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2016065159A1 (en) * 2014-10-23 2016-04-28 Washington University Systems and methods for measuring cardiac strain

Patent Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102930567A (zh) * 2012-09-25 2013-02-13 电子科技大学 多核加权最小二乘支撑向量机的磁共振并行成像重建方法
CN104282036A (zh) * 2013-07-09 2015-01-14 韦伯斯特生物官能(以色列)有限公司 由稀疏样本对心脏进行的基于模型的重构
CN103793711A (zh) * 2014-01-17 2014-05-14 首都医科大学 一种基于脑部核磁共振图像的多维度纹理提取方法
CN107730516A (zh) * 2017-10-16 2018-02-23 江南大学 一种基于模糊聚类的脑mr影像分割方法
CN107991636A (zh) * 2017-11-24 2018-05-04 哈尔滨工业大学 一种基于适应性结构低秩矩阵的快速磁共振图像重建方法
CN108198193A (zh) * 2018-01-16 2018-06-22 北京航空航天大学 一种利用改进直觉模糊聚类算法分割红外舰船图像的方法
CN109920017A (zh) * 2019-01-16 2019-06-21 昆明理工大学 基于特征向量的自一致性的联合全变分Lp伪范数的并行磁共振成像重构方法
CN109934884A (zh) * 2019-01-22 2019-06-25 昆明理工大学 一种基于变换学习和联合稀疏性的迭代自一致性并行成像重构方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
基于深度学习的快速磁共振成像方法研究;肖韬辉;《中国优秀博硕士学位论文全文数据库(硕士)信息科技辑》;20190215(第02期);I140-57 *

Also Published As

Publication number Publication date
CN111754598A (zh) 2020-10-09

Similar Documents

Publication Publication Date Title
CN108734659B (zh) 一种基于多尺度标签的亚像素卷积图像超分辨率重建方法
CN108460726B (zh) 一种基于增强递归残差网络的磁共振图像超分辨重建方法
CN111754598B (zh) 基于变换学习的局部空间邻域并行磁共振成像重构方法
WO2018099321A1 (zh) 一种基于广义树稀疏的权重核范数磁共振成像重建方法
CN104933683B (zh) 一种用于磁共振快速成像的非凸低秩重建方法
CN110148215B (zh) 一种基于平滑约束和局部低秩约束模型的四维磁共振图像重建方法
CN107274462B (zh) 基于熵和几何方向的分类多字典学习磁共振图像重建方法
CN107991636B (zh) 一种基于适应性结构低秩矩阵的快速磁共振图像重建方法
CN107993194B (zh) 一种基于平稳小波变换的超分辨率重建方法
CN112819949B (zh) 一种基于结构化低秩矩阵的磁共振指纹图像重建方法
CN112991483B (zh) 一种非局部低秩约束的自校准并行磁共振成像重构方法
CN111487573B (zh) 一种用于磁共振欠采样成像的强化型残差级联网络模型
CN109920017B (zh) 基于特征向量的自一致性的联合全变分Lp伪范数的并行磁共振成像重构方法
CN109934884B (zh) 一种基于变换学习和联合稀疏性的迭代自一致性并行成像重构方法
CN112617798B (zh) 一种基于Lp范数联合全变分的并行磁共振成像重构方法
CN112734869A (zh) 基于稀疏复数u型网络的快速磁共振成像方法
CN109188327B (zh) 基于张量积复小波紧框架的磁共振图像快速重构方法
He et al. Dynamic MRI reconstruction exploiting blind compressed sensing combined transform learning regularization
CN114140404A (zh) 基于人工智能的肺部多核mri双域超分辨率重建方法
CN114004764B (zh) 一种基于稀疏变换学习的改进灵敏度编码重建方法
CN113674379A (zh) 一种基于参考支撑集的共稀疏分析模型的mri重构方法、系统及计算机可读存储介质
CN108346167B (zh) 一种基于正交字典下同时稀疏编码的mri图像重构方法
CN113628298B (zh) 基于特征向量的自一致性和非局部低秩的并行mri重构方法
CN113034639B (zh) 一种基于可分离汉克尔矩阵的磁共振成像的图像重建方法
CN113892938B (zh) 一种基于非局部低秩约束的改进灵敏度编码重构方法

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant