CN112991483B - 一种非局部低秩约束的自校准并行磁共振成像重构方法 - Google Patents

一种非局部低秩约束的自校准并行磁共振成像重构方法 Download PDF

Info

Publication number
CN112991483B
CN112991483B CN202110452944.2A CN202110452944A CN112991483B CN 112991483 B CN112991483 B CN 112991483B CN 202110452944 A CN202110452944 A CN 202110452944A CN 112991483 B CN112991483 B CN 112991483B
Authority
CN
China
Prior art keywords
representing
iteration
coil
image
reconstructed
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
CN202110452944.2A
Other languages
English (en)
Other versions
CN112991483A (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 CN202110452944.2A priority Critical patent/CN112991483B/zh
Publication of CN112991483A publication Critical patent/CN112991483A/zh
Application granted granted Critical
Publication of CN112991483B publication Critical patent/CN112991483B/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

Landscapes

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

Abstract

本发明涉及一种非局部低秩约束的自校准并行磁共振成像重构方法，属于磁共振成像技术领域。并行磁共振成像重构一直是近几年的研究热点。本发明基于SPIRiT重构框架，提出一种图像非局部低秩约束(nonlocal low‑rank，NLR)的高质量磁共振成像重构模型，命名为NLR‑SPIRiT，利用交替方向乘子法(alternating direction method of multipliers，ADMM)进行求解。实验比较了结合了联合全变分(Joint Total Variation,JTV)正则项的迭代自一致并行成像JTV‑SPIRiT模型。实验结果表明，NLR‑SPIRiT算法重构出的图像有更好的重构质量，能更好的去除欠采样伪影。

Description

一种非局部低秩约束的自校准并行磁共振成像重构方法
技术领域
本发明涉及一种非局部低秩约束的自校准并行磁共振成像重构方法，属于磁共振成像技术领域。
背景技术
磁共振成像(Magnetic Resonance Imaging,MRI)是一种主要的医学成像技术，对于学术研究和临床应用都有很大的价值。压缩感知(Compressed Sensing，CS)与并行成像(Parallel Imaging，PI)技术都能通过减少k空间数据的采集量来加快磁共振成像速度。如何从低于奈奎斯特采样率的采样数据下重构出高质量的图像，一直是磁共振成像领域的研究热点。
在压缩感知重构中，使用正则化约束探索图像先验信息，通过非线性重建算法重构图像。探索图像先验信息的正则化有：利用图像局部特性的全变分(Total Variation，TV)，利用固定变换下存在稀疏性的l1范数，以及一些自适应的变换方法，例如字典学习(Dictionary Learning，DL)和变换学习(Transform Learning，TL)。
除了稀疏性之外，图像还包含非局部结构，例如自相似性：补丁通常与同一图像中的其他非局部结构相似。在压缩感知中使用块匹配(Block Match，BM)方法获取非局部组稀疏性，之后通过低秩方法解决最小化问题，得到了基于非局部低秩正则化的压缩感知模型(Compressive Sensing via Nonlocal Low-Rank Regularization，NLR-CS)，与之前的方法相比具有显著的提高。
并行成像也是一种流行的加速磁共振算法。早期的并行成像利用线圈的不同的灵敏度来获取k空间数据。如基于线圈灵敏度信息的灵敏度编码(SENSitivity Encoding，SENSE)，它是一种直接的利用灵敏度信息进行并行采样图像恢复的方法，在具有准确的灵敏度信息时允许最佳的重构结果。Uecker等人提出特征向量算子的迭代自一致性的并行成像重构(iTerative Self-consistent Parallel Imaging Reconstruction usingEigenvector maps,ESPIRiT)模型，通过间接的方式，应用k空间算子的主要特征向量来表示线圈灵敏度算子，这些灵敏度算子可以通过在图像空间中的使用奇异值分解来快速计算，从而能够从k空间中心的自校准区域中有效的估计线圈灵敏度，是一种有效的图像并行重构方案。
然而，在实际应用中高精度的灵敏度信息测量很困难，迭代自一致性并行成像重构(Iterative self-consistent parallel imaging reconstruction，SPIRiT)方法是一类基于自校准的重构方法，只需要利用少量中心全采样信号(ACS)进行校准，无需直接使用线圈灵敏度，也就避免了灵敏度估计困难的问题。但是现有的方法仅使用了简单的将联合全变分(Joint Total Variation，JTV)和联合L1范数(Joint L1-norm，JL1)正则化项，算法重构图像的图像质量还有待提高。
发明内容
本发明的目的在于克服现有技术的不足，提供一种非局部低秩约束的自校准并行磁共振成像重构方法，能进一步提高磁共振成像重构图像的质量。
本发明采用的技术方案是：一种非局部低秩约束的自校准并行磁共振成像重构方法，包括以下步骤：
S0：初始化，令
Figure BDA0003038957100000021
Z0＝0，B0＝0，D0＝0，z0＝0，d0＝0，b0＝0；
其中，上标“0”表示迭代前的初始值。
Figure BDA0003038957100000022
表示待重构的多线圈图像，X0表示X的初始值，
Figure BDA0003038957100000023
表示欠采样的k空间数据，X与Y的每一列由单线圈数据Xc与Yc按列堆叠得到，
Figure BDA0003038957100000024
表示待重构的第c个线圈图像，
Figure BDA0003038957100000025
表示第c个线圈欠采样k空间数据，c＝1,...,C表示线圈索引变量。N表示待重构单线圈图像的像素点数，M表示单线圈图像的k空间欠采样投影数据个数，C表示线圈数。
Figure BDA0003038957100000026
表示傅里叶变换，
Figure BDA0003038957100000027
表示
Figure BDA0003038957100000028
的共轭转置，“H”表示共轭转置运算。
Figure BDA0003038957100000029
表示辅助变量,
Figure BDA00030389571000000210
和
Figure BDA00030389571000000211
表示辅助变量Z和B的第c列，Z0和B0表示Z和B的初始值。
Figure BDA00030389571000000212
和
Figure BDA00030389571000000213
表示对应于B和Z的拉格朗日乘子，
Figure BDA0003038957100000031
和
Figure BDA0003038957100000032
表示z和b的第c列，z0和b0表示z和b的初始值。D＝[D1,:,D2,:,...,DC-1,:,DC,:]表示辅助变量，Dc,:表示D的第c个分块，
Figure BDA0003038957100000033
表示Dc,:的第i列，D0表示D的初始值。d＝[d1,:,d2,:,...,dC-1,:,dC,:]表示对应于D的拉格朗日乘子，dc,:表示d的第c个分块，
Figure BDA0003038957100000034
表示dc,:的第i列，d0表示d的初始值；
定义Dci＝Vci(Xc)，其中
Figure BDA00030389571000000314
是块匹配运算符，具体地说：对Xc取补丁得到重叠的Np个大小为
Figure BDA0003038957100000035
的参考块，n表示参考块中的像素点个数，对于每一个参考块
Figure BDA0003038957100000036
c＝1,...,C，i＝1,...,Np，i表示参考块索引变量，在大小为s×s的搜索窗口中找到与参考块uci最相似(欧几里得距离最小)的m个块，将提取的相似块按照欧几里得距离升序排列构成
Figure BDA0003038957100000037
的列。
S1：预处理：计算Γ-1，其中Γ定义为：
Γ＝μ1(G-I)H(G-I)+β1I
其中，G是用SPIRiT核卷积整个k空间的插值算子，μ1＞0和β1＞0为正则化参数，I为单位矩阵。对矩阵重新排序会产生由Nc×Nc个块组成的块对角线结构，通过直接求每个Nc×Nc块的逆来计算Γ-1，“-1”表示求逆运算。
S2:初始化，迭代次数k＝1；
S3:初始化，c＝1；
S4:计算第k+1次迭代的辅助变量
Figure BDA0003038957100000038
计算公式如下：
Figure BDA0003038957100000039
其中，
Figure BDA00030389571000000310
表示第k次迭代的待重构第c个线圈图像，
Figure BDA00030389571000000311
表示第k次迭代中对应于
Figure BDA00030389571000000312
的拉格朗日乘子。
S5:初始化，i＝1；
S6:对Xc进行相似块匹配得到
Figure BDA00030389571000000313
S7:对
Figure BDA0003038957100000041
进行奇异值分解，得到
Figure BDA0003038957100000042
其中，∑为奇异值，U和V是特征向量，
Figure BDA0003038957100000043
表示第k次迭代中对应于
Figure BDA0003038957100000044
的拉格朗日乘子，上标“T”表示矩阵的转置运算，∑＝diag(γ1,γ2,...,γJ-1,γJ)，令γj表示Dci的第j个奇异值，j＝1,...,J表示奇异值索引变量，J＝min(m,n)，。
S8:计算权重w＝[w1,w2,...,wJ-1,wJ]；
其中，w包含J个值，令wj表示第j个权重，计算式为
Figure BDA0003038957100000045
式中ε是一个小的常数，用于保证分母不为0。
S9:计算第k+1次迭代的辅助变量
Figure BDA0003038957100000046
计算公式如下：
Figure BDA0003038957100000047
其中，
Figure BDA0003038957100000048
是第k+1次迭代的辅助变量，
Figure BDA0003038957100000049
为奇异值阈值，β2＞0和μ2＞0表示参数，符号：(x)+＝max{x,0}。
S10:判断是否完成每一个线圈中Np个
Figure BDA00030389571000000410
的更新，当i＝Np时进入步骤S11；否则，令i＝i+1后返回S5；
S11：计算第k+1次迭代的辅助变量
Figure BDA00030389571000000411
计算公式如下：
Figure BDA00030389571000000412
其中，
Figure BDA00030389571000000413
是第k+1次迭代的辅助变量，
Figure BDA00030389571000000414
是第k次迭代中对应于
Figure BDA00030389571000000415
的拉格朗日乘子，β2＞0和ν＞0表示参数。
Figure BDA00030389571000000416
表示将去噪得到的
Figure BDA00030389571000000417
使用添加的方法放回块匹配之前的位置，
Figure BDA00030389571000000418
表示Xc的第i个像素点在Vci(Xc)中出现的次数，也称为Vci(Xc)对Xc的第i个像素点引用次数。
S12：计算第k+1次迭代的第c个线圈的待重构图像
Figure BDA00030389571000000419
计算公式如下：
Figure BDA0003038957100000051
其中，
Figure BDA0003038957100000052
表示k空间欠采样算子，
Figure BDA0003038957100000053
为矩阵
Figure BDA0003038957100000054
的共轭，
Figure BDA0003038957100000055
为对角矩阵，
Figure BDA0003038957100000056
表示第c个线圈的k空间采集数据。
S13:初始化，i＝1；
S14:计算拉格朗日乘子
Figure BDA0003038957100000057
其中，
Figure BDA0003038957100000058
是第k+1次迭代中对应于
Figure BDA0003038957100000059
的拉格朗日乘子，η2＞0表示参数。
S15:判断是否完成每一个线圈中Np个
Figure BDA00030389571000000510
的更新，当i＝Np时进入步骤S16；否则，令i＝i+1后返回S13；
S16:计算拉格朗日乘子
Figure BDA00030389571000000511
其中，
Figure BDA00030389571000000512
是第k+1次迭代中对应于
Figure BDA00030389571000000513
的拉格朗日乘子，η1＞0表示参数。
S17：计算拉格朗日乘子
Figure BDA00030389571000000514
其中，
Figure BDA00030389571000000515
是第k+1次迭代中对应于
Figure BDA00030389571000000516
的拉格朗日乘子，η3＞0表示参数。
S18：判断是否完成所有线圈的计算，若c＝C，进入步骤S19，否则，令c＝c+1返回S3；
S19：判断是否达到最大迭代次数K，若k＝K，进入步骤S20；否则，令k＝k+1返回S2；
S20：输出重构的线圈图像Xc，c＝1,...,C，使用平方和的平方根方法来形成最终重构的图像
Figure BDA00030389571000000517
计算公式如下：
Figure BDA00030389571000000518
本发明的有益效果是：迭代自一致性并行成像重建(iterative self-consistentparallel imaging reconstruction,SPIRiT)是一种利用k空间核校准的自校准并行重构方法，将多线圈数据重构问题转化为含有数据一致性与校准一致性的正则化的最小二乘优化问题。基于SPIRiT模型，本发明提出了一种采用图像非局部低秩作为约束项的高质量并行成像重建算法。首先，本发明使用变量分离和ADMM算法将重构问题转化成一个去噪子问题和两个最小二乘子问题。之后，对每个子问题进行求解：关于非局部低秩的去噪子问题，使用加权核范数作为秩代理转化为低秩最小化问题求解；关于SPIRiT校准最小二乘子问题，使用预处理共轭梯度方法进行求解；重构图像的最小二乘子问题再次应用变量分离和ADMM算法进行求解。理论分析和成像实验表明，与之前的重构算法相比，本申请提出的新的并行成像算法可以有效地提高重构图像的质量。
附图说明
图1为本发明方法流程图；
图2为使用了八通道的膝盖线圈图像进行完全采样的来获取受试者的体内人脑数据(即dataset1)图；
图3为具有6×采样率和24×24ACS线的二维泊松圆盘欠采样掩模图；
图4-图5展示了在数据集dataset1下，分别使用算法JTV-SPIRiT和算法NLR-SPIRIT从含有6×加速度和24×24中心校准二维泊松圆盘欠采样数据重构的图像；
图6-图7中展示了对应于图4-图5的误差图像。
具体实施方式
下面结合附图进一步详细描述本发明的技术方案：
实施例1：本发明是基于SPIRiT框架提出的一种高效的重构方法。
假设
Figure BDA0003038957100000061
表示待重构的多线圈图像，
Figure BDA0003038957100000062
表示欠采样的k空间数据。N表示待重构单线圈图像的像素点数，M表示单线圈图像的k空间欠采样投影数据个数，C为线圈数。X与Y的每一列由单线圈图像按列堆叠得到。多线圈图像的欠采样k空间数据Y由下式给出：
Y＝AX (1)
其中，矩阵
Figure BDA0003038957100000063
表示对多线圈数据X进行欠采样投影的算子，
Figure BDA0003038957100000064
为欠采样算子，
Figure BDA0003038957100000065
为傅里叶变换。
SPIRiT对图像网格的每个点及其整个邻域点(无论它们是否被获取)之间执行一致性，所有k空间位置的校准一致性公式为：
X＝GX (2)
其中G是用SPIRiT核gic卷积整个k空间的插值算子。gic是插值核，也称为SPIRiT核，是根据以全采样的ACS数据进行标定的。
尽管结合了联合L1范数(JL1)正则项和联合总变化(JTV)正则项的SPIRiT算法可有效改善重构质量，但仍有较大的改进余地。为了进一步提高重构质量，本发明将非局部低秩正则化(NLR)与SPIRiT模型相结合，提出NLR-SPIRIT算法，得到如下最优化问题：
Figure BDA0003038957100000071
SPIRiT是一种逐线圈求解的方法，为了便于理解，之后介绍对每个线圈数据Xc重构问题的求解：
Figure BDA0003038957100000072
在式(3)和式(4)中，
Figure BDA0003038957100000073
c＝1,...,C表示待重构的第c个线圈的图像，X＝[X1,X2,...,XC-1,XC]，c表示线圈索引。
Figure BDA0003038957100000074
表示第c个线圈图像的欠采样k空间数据，Y＝[Y1,Y2,...,YC-1,YC]。μ1与μ2为正则化参数。
rank(Vci(Xc))表示图像非局部低秩正则项。其中
Figure BDA0003038957100000075
是块匹配运算符，具体地说：对Xc取补丁块得到重叠的Np个大小为
Figure BDA0003038957100000076
的块，n表示参考块中的像素点个数。对于每一个参考块
Figure BDA0003038957100000077
c＝1,...,C，i＝1,...,Np，i表示参考块索引变量。在大小为s×s的搜索窗口中找到与参考块uci最相似(欧几里得距离最小)的m个块，对提取的相似块按照欧几里得距离升序排列构成
Figure BDA0003038957100000078
的列。
算法推导：
引入辅助变量
Figure BDA0003038957100000079
与辅助变量
Figure BDA00030389571000000710
以及对应的拉格朗日乘子
Figure BDA00030389571000000711
和
Figure BDA00030389571000000712
使用ADMM方法，可以将公式(4)转化为以下子公式进行求解：
Figure BDA00030389571000000713
Figure BDA0003038957100000081
Figure BDA0003038957100000082
Figure BDA0003038957100000083
Figure BDA0003038957100000084
在式(5)-式(9)中，变量的上标“k+1”和“k”表示变量为第k+1次迭代与第k次迭代的变量。在式(5)中，辅助变量
Figure BDA0003038957100000085
是第k+1次迭代的辅助变量，
Figure BDA0003038957100000086
是第k次迭代中对应于
Figure BDA0003038957100000087
的拉格朗日乘子，
Figure BDA0003038957100000088
表示第k次迭代的待重构第c个线圈图像，β1＞0为正则项
Figure BDA0003038957100000089
对应的参数。在式(6)中，
Figure BDA00030389571000000810
每个线圈的图像重构中含有Np个Dci，
Figure BDA00030389571000000811
是第k+1次迭代的辅助变量，
Figure BDA00030389571000000812
是第k次迭代中对应于
Figure BDA00030389571000000813
的拉格朗日乘子，β2＞0为正则项
Figure BDA00030389571000000814
对应的参数。在式(7)中，
Figure BDA00030389571000000815
表示第k+1次迭代的待重构第c个线圈的图像。式(8)和(9)中，
Figure BDA00030389571000000816
是第k+1次迭代中对应于
Figure BDA00030389571000000817
的拉格朗日乘子，
Figure BDA00030389571000000818
是第k+1次迭代中对应于
Figure BDA00030389571000000819
的拉格朗日乘子，η1＞0与η2＞0表示更新拉格朗日乘子的参数。
公式(8)和(9)可以直接进行求解，之后对上述其他几个公式进行计算：
关于Zc的子问题：根据式子(5)，得到Zc的解为：
Figure BDA00030389571000000820
使用预处理方法可以简单计算得Z的精确解。具体地说，令：
Γ＝μ1(G-I)H(G-I)+β1I (11)
事实上，矩阵Γ具有对角线块结构，对矩阵进行简单的排序会产生由Nc×Nc个块组成的块对角结构。只要Nc足够小，就可以通过对每个Nc×Nc的块进行翻转得到Γ-1。Γ-1算子可以在循环之前计算得到，那么就有Zc的精确解为：
Figure BDA0003038957100000091
关于Dci的子问题：一般来说，低秩极小化公式(6)是一个非确定多项式困难(nondeterministic polynomial time(NP)-hard)问题，因此不能直接进行求解。为了得到低秩近似，可以使用加权核范数进行代理求解，加权核范数定义如下：
Figure BDA0003038957100000092
其中wj为权重。使用加权核范数作为秩代理，可以将式(6)写为：
Figure BDA0003038957100000093
公式(14)可以使用加权奇异值阈值算法进行求解，求解过程如下：对X进行相似块匹配，得到对应的低秩近似Vci(Xc)，对
Figure BDA0003038957100000094
进行奇异值分解得到U∑VT，其中∑为奇异值，U和V是特征向量，∑＝diag(γ1,γ2,...,γJ-1,γJ)，J＝min(m,n)，令γj表示Vci(Xc)+dci的第j个奇异值。那么，可以得到Dci的最优解可以由加权奇异值算法给出：
Figure BDA0003038957100000095
其中
Figure BDA0003038957100000096
为奇异值阈值，w＝[w1,w2,...,wJ-1,wJ]表示权重，令wj表示第j个权重，则
Figure BDA0003038957100000097
其中ε是一个小的常数，保证分母不为0，符号(x)+＝max{x,0}。
关于Xc的子问题：根据式(7)，Xc的封闭形式解由下式给出：
Figure BDA0003038957100000098
上式无法直接求解，需要再次使用变量分离以及ADMM，引入辅助变量
Figure BDA0003038957100000099
以及对应的对偶变量
Figure BDA00030389571000000910
Xc的子问题可以通过以下式子的迭代求解：
Figure BDA0003038957100000101
Figure BDA0003038957100000102
Figure BDA0003038957100000103
在式(17)中，
Figure BDA0003038957100000104
是第k+1次迭代的辅助变量，
Figure BDA0003038957100000105
是第k次迭代中对应于
Figure BDA0003038957100000106
的拉格朗日乘子，ν＞0表示新引入的正则化
Figure BDA0003038957100000107
的参数。在式(19)中，
Figure BDA0003038957100000108
是第k+1次迭代中对应于
Figure BDA0003038957100000109
的拉格朗日乘子，η3＞0表示参数。
公式(19)可以直接进行求解，之后对上述其他两个公式进行计算：
对于公式(17)，Bc的封闭解为：
Figure BDA00030389571000001010
其中，
Figure BDA00030389571000001011
表示将去噪得到的
Figure BDA00030389571000001012
使用添加的方法放回块匹配之前的位置，
Figure BDA00030389571000001013
表示Xc的第i个像素点在Vci(Xc)中出现的次数，也称为Vci(Xc)对Xc的第i个像素点引用次数。
关于公式(18)：代入
Figure BDA00030389571000001014
得到Xc的精确解为：
Figure BDA00030389571000001015
其中，
Figure BDA00030389571000001016
表示k空间欠采样算子，
Figure BDA00030389571000001017
为矩阵
Figure BDA00030389571000001018
的共轭，
Figure BDA00030389571000001019
为对角矩阵，
Figure BDA00030389571000001020
表示第c个线圈的k空间采集数据。
综上所述，所有子问题均可求解。得到重构的每一个线圈的(22)图像Xc，c＝1,...,C后，使用平方和的平方根方法来形成最终重构的图像
Figure BDA00030389571000001021
计算公式如下：
Figure BDA00030389571000001022
具体流程如图1所示，其中步骤如下：
S0：令
Figure BDA0003038957100000111
Z0＝0，B0＝0，D0＝0，z0＝0，d0＝0，b0＝0；
上标“0”表示迭代前的初始值。
Figure BDA0003038957100000112
Figure BDA0003038957100000113
表示辅助变量,
Figure BDA0003038957100000114
和
Figure BDA0003038957100000115
表示辅助变量Z和B的第c列。
Figure BDA0003038957100000116
和
Figure BDA0003038957100000117
表示对应于B和Z的拉格朗日乘子，
Figure BDA0003038957100000118
和
Figure BDA0003038957100000119
表示z和b的第c列。D＝[D1,:,D2,:,...,DC-1,:,DC,:]表示辅助变量，Dc,:表示D的第c个分块，
Figure BDA00030389571000001110
表示Dc,:的第i列，D0表示D的初始值。d＝[d1,:,d2,:,...,dC-1,:,dC,:]表示对应于D的拉格朗日乘子，dc,:表示d的第c个分块，
Figure BDA00030389571000001111
表示dc,:的第i列，d0表示d的初始值
S1：预处理：计算Γ-1；
S2:初始化，迭代次数k＝1；
S3:初始化，c＝1；
S4:计算
Figure BDA00030389571000001112
公式如(12)；
S5:初始化，i＝1；
S6:对Xc进行相似块匹配得到
Figure BDA00030389571000001113
S7:对
Figure BDA00030389571000001114
进行奇异值分解，得到
Figure BDA00030389571000001115
S8:计算权重w＝[w1,w2,...,wJ-1,wJ]，第j个权重
Figure BDA00030389571000001116
S9:计算
Figure BDA00030389571000001117
公式如(15)；
S10:判断是否完成每一个线圈中Np个
Figure BDA00030389571000001118
的更新，当i＝Np时进入步骤S11；否则，令i＝i+1后返回S5；
S11：计算
Figure BDA00030389571000001119
公式如(20)；
S12：计算
Figure BDA00030389571000001120
公式如(21)；
S13:初始化，i＝1；
S14:计算拉格朗日乘子
Figure BDA0003038957100000121
公式如(9)；
S15:判断是否完成每一个线圈中Np个
Figure BDA0003038957100000122
的更新，当i＝Np时进入步骤S16；否则，令i＝i+1后返回S13；
S16:计算拉格朗日乘子
Figure BDA0003038957100000123
公式如(8)；
S17：计算拉格朗日乘子
Figure BDA0003038957100000124
公式如(19)；
S18：判断是否完成所有线圈的计算，若c＝C，进入步骤S19，否则，令c＝c+1返回S3；
S19：判断是否达到最大迭代次数K，若k＝K，进入步骤S20；否则，令k＝k+1返回S2；
S20：输出重构的线圈图像Xc，c＝1,...,C，使用平方和的平方根方法来形成最终输出的重构图像x，计算公式如(23)。
实验结果：
在下面的实验中，为了验证本申请所提出的方法的性能，本申请将它与结合联合全变分的迭代自一致并行成像JTV-SPIRiT方法进行了比较。所有算法均用MATLAB实现。所有的实验均在笔记本电脑上进行，笔记本电脑的配置为：因特尔内核[email protected]双核,12GB内存，Windows 7***(64bit)。
为了比较各算法的重构性能，本申请使用了1张人类膝盖图来进行仿真实验，命名为dataset1(如图2)。为了生成测试数据集，使用具有R×(R＝3:8)采样率和24×24ACS线的二维泊松圆盘欠采样掩模对完全采样的数据集进行人工欠采样得到CS并行测量值，如图3所示(R＝6)。
比较算法的正则化参数是针对SNR最优进行调整的。在以下实验中，信噪比(Signal Noise Ratio,SNR)和高频误差范数(High-Frequency Error Norm,HFEN)用于定量评估重建图像的质量。SNR数值越高，HFEN数值越低，则图像重构质量越好。对于参考图像x与重构图像
Figure BDA0003038957100000125
SNR(dB)的定义为：
Figure BDA0003038957100000126
其中Var表示参考图像x与重构图像
Figure BDA0003038957100000127
之间的均方误差，MSE为参考图像x的方差。
HEFN定义为：
Figure BDA0003038957100000131
其中filter(·)是一个拉普拉斯高斯滤波器，用于捕捉图像边缘。
本申请比较了加速因子为6时，使用各个算法对所选择的数据集下进行重构的视觉呈现。图4和图5展示了在数据集dataset1下，使用算法JTV-SPIRiT和算法NLR-SPIRIT的重构图像。从图4和图5可以看出，使用JTV-SPIRiT算法的重构图像在膝盖关节处有明显的模糊伪影，使用NLR-SPIRiT的重构图像更加清楚，更接近于原始图像。
为了更直观地比较几种算法重构效果，本申请在图6和图7中展示了JTV-SPIRiT算法和NLR-SPIRiT算法在数据集dataset1下的误差图像(白点表示误差，白点越多，误差越大)。从图6和图7可以明显看出，JTV-SPIRiT的重构图像有较大的误差，NLR-SPIRiT具有较小的重构误差，拥有更好重构质量。
另外，本申请还分析了几种方法的重构图像的SNR与HFEN的数值分析。JTV-SPIRiT方法的重构图像(如图4)的SNR＝13.45dB，HFEN＝0.3796；NLR-SPIRiT方法的重构图像(如图5)的SNR＝14.97dB，HFEN＝0.2813。可以明显看出，在图像指标下，本申请提出的NLR-SPIRIT方法具有最佳的SNR与HFEN。
本发明将非局部低秩约束与SPIRiT模型结合，提出了一种新的并行磁共振成像方法，利用非局部低秩约束探索图像的自相似特性，使用SPIRiT模型的数据校准方法来重构图像，本申请的方法称为非局部低秩约束的迭代自一致性并行成像重构(NLR-SPIRIT)模型。进行了多组试验后，实验结果显示，本申请提出的NLR-SPIRiT方法比JTV-SPIRiT方法的成像结果更好。
综上所述，使用本申请选用的测试集对不同算法进行实验，使用基于非局部低秩正则化的NLR-SPIRiT方法在两种不同的评估指标以及视觉上均优于JTV-SPIRiT方法。
以上结合附图对本发明的具体实施方式作了详细说明，但是本发明并不限于上述实施方式，在本领域普通技术人员所具备的知识范围内，还可以在不脱离本发明宗旨的前提下作出各种变化。

Claims (1)

1.一种非局部低秩约束的自校准并行磁共振成像重构方法，包括以下步骤：
S0：初始化，令
Figure FDA0003038957090000011
Z0＝0，B0＝0，D0＝0，z0＝0，d0＝0，b0＝0；
其中，上标“0”表示迭代前的初始值，
Figure FDA0003038957090000012
表示待重构的多线圈图像，X0表示X的初始值，
Figure FDA0003038957090000013
表示欠采样的k空间数据，X与Y的每一列由单线圈数据Xc与Yc按列堆叠得到，
Figure FDA0003038957090000014
表示待重构的第c个线圈图像，
Figure FDA0003038957090000015
表示第c个线圈欠采样k空间数据，c＝1,...,C表示线圈索引变量，N表示待重构单线圈图像的像素点数，M表示单线圈图像的k空间欠采样投影数据个数，C表示线圈数，
Figure FDA0003038957090000016
表示傅里叶变换，
Figure FDA0003038957090000017
表示
Figure FDA0003038957090000018
的共轭转置，“H”表示共轭转置运算；
Figure FDA0003038957090000019
表示辅助变量,
Figure FDA00030389570900000110
和
Figure FDA00030389570900000111
表示辅助变量Z和B的第c列，Z0和B0表示Z和B的初始值，
Figure FDA00030389570900000112
和
Figure FDA00030389570900000113
表示对应于B和Z的拉格朗日乘子，
Figure FDA00030389570900000114
和
Figure FDA00030389570900000115
表示z和b的第c列，z0和b0表示z和b的初始值，D＝[D1,:,D2,:,...,DC-1,:,DC,:]表示辅助变量，Dc,:表示D的第c个分块，
Figure FDA00030389570900000116
Figure FDA00030389570900000117
表示Dc,:的第i列，D0表示D的初始值，d＝[d1,:,d2,:,...,dC-1,:,dC,:]表示对应于D的拉格朗日乘子，dc,:表示d的第c个分块，
Figure FDA00030389570900000118
Figure FDA00030389570900000119
表示dc,:的第i列，d0表示d的初始值；
定义Dci＝Vci(Xc)，其中
Figure FDA00030389570900000120
是块匹配运算符，具体地说：对Xc取补丁得到重叠的Np个大小为
Figure FDA00030389570900000121
的参考块，n表示参考块中的像素点个数，对于每一个参考块
Figure FDA00030389570900000122
c＝1,...,C，i＝1,...,Np，i表示参考块索引变量，在大小为s×s的搜索窗口中找到与参考块uci欧几里得距离最小的m个块，将提取的相似块按照欧几里得距离升序排列构成
Figure FDA00030389570900000123
的列；
S1：预处理：计算Γ-1，其中Γ定义为：
Γ＝μ1(G-I)H(G-I)+β1I
其中，G是用SPIRiT核卷积整个k空间的插值算子，μ1＞0和β1＞0为正则化参数，I为单位矩阵，对矩阵重新排序会产生由Nc×Nc个块组成的块对角线结构，通过直接求每个Nc×Nc块的逆来计算Γ-1，“-1”表示求逆运算；
S2:初始化，迭代次数k＝1；
S3:初始化，c＝1；
S4:计算第k+1次迭代的辅助变量
Figure FDA0003038957090000021
计算公式如下：
Figure FDA0003038957090000022
其中，
Figure FDA0003038957090000023
表示第k次迭代的待重构第c个线圈图像，
Figure FDA0003038957090000024
表示第k次迭代中对应于
Figure FDA0003038957090000025
的拉格朗日乘子；
S5:初始化，i＝1；
S6:对Xc进行相似块匹配得到
Figure FDA0003038957090000026
S7:对
Figure FDA0003038957090000027
进行奇异值分解，得到
Figure FDA0003038957090000028
其中，∑为奇异值，U和V是特征向量，
Figure FDA0003038957090000029
表示第k次迭代中对应于
Figure FDA00030389570900000210
的拉格朗日乘子，上标“T”表示矩阵的转置运算，∑＝diag(γ1,γ2,...,γJ-1,γJ)，令γj表示Dci的第j个奇异值，j＝1,...,J表示奇异值索引变量，J＝min(m,n)；
S8:计算权重w＝[w1,w2,...,wJ-1,wJ]；
其中，w包含J个值，令wj表示第j个权重，计算式为
Figure FDA00030389570900000211
式中ε是一个小的常数，用于保证分母不为0；
S9:计算第k+1次迭代的辅助变量
Figure FDA00030389570900000212
计算公式如下：
Figure FDA00030389570900000213
其中，
Figure FDA0003038957090000031
是第k+1次迭代的辅助变量，
Figure FDA0003038957090000032
为奇异值阈值，β2＞0和μ2＞0表示参数，符号：(x)+＝max{x,0}；
S10:判断是否完成每一个线圈中Np个
Figure FDA0003038957090000033
的更新，当i＝Np时进入步骤S11；否则，令i＝i+1后返回S5；
S11：计算第k+1次迭代的辅助变量
Figure FDA0003038957090000034
计算公式如下：
Figure FDA0003038957090000035
其中，
Figure FDA0003038957090000036
是第k+1次迭代的辅助变量，
Figure FDA0003038957090000037
是第k次迭代中对应于
Figure FDA0003038957090000038
的拉格朗日乘子，β2＞0和ν＞0表示参数，
Figure FDA0003038957090000039
表示将去噪得到的
Figure FDA00030389570900000310
使用添加的方法放回块匹配之前的位置，
Figure FDA00030389570900000311
表示Xc的第i个像素点在Vci(Xc)中出现的次数，也称为Vci(Xc)对Xc的第i个像素点引用次数；
S12：计算第k+1次迭代的第c个线圈的待重构图像
Figure FDA00030389570900000312
计算公式如下：
Figure FDA00030389570900000313
其中，
Figure FDA00030389570900000314
表示k空间欠采样算子，
Figure FDA00030389570900000315
为矩阵
Figure FDA00030389570900000316
的共轭，
Figure FDA00030389570900000317
为对角矩阵，
Figure FDA00030389570900000318
表示第c个线圈的k空间采集数据；
S13:初始化，i＝1；
S14:计算拉格朗日乘子
Figure FDA00030389570900000319
其中，
Figure FDA00030389570900000320
是第k+1次迭代中对应于
Figure FDA00030389570900000321
的拉格朗日乘子，η2＞0表示参数；
S15:判断是否完成每一个线圈中Np个
Figure FDA00030389570900000322
的更新，当i＝Np时进入步骤S16；否则，令i＝i+1后返回S13；
S16:计算拉格朗日乘子
Figure FDA00030389570900000323
其中，
Figure FDA00030389570900000324
是第k+1次迭代中对应于
Figure FDA00030389570900000325
的拉格朗日乘子，η1＞0表示参数，
S17：计算拉格朗日乘子
Figure FDA0003038957090000041
其中，
Figure FDA0003038957090000042
是第k+1次迭代中对应于
Figure FDA0003038957090000043
的拉格朗日乘子，η3＞0表示参数；
S18：判断是否完成所有线圈的计算，若c＝C，进入步骤S19，否则，令c＝c+1返回S3；
S19：判断是否达到最大迭代次数K，若k＝K，进入步骤S20；否则，令k＝k+1返回S2；
S20：输出重构的线圈图像Xc，c＝1,...,C，使用平方和的平方根方法来形成最终重构的图像
Figure FDA0003038957090000044
计算公式如下：
Figure FDA0003038957090000045
CN202110452944.2A 2021-04-26 2021-04-26 一种非局部低秩约束的自校准并行磁共振成像重构方法 Active CN112991483B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110452944.2A CN112991483B (zh) 2021-04-26 2021-04-26 一种非局部低秩约束的自校准并行磁共振成像重构方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110452944.2A CN112991483B (zh) 2021-04-26 2021-04-26 一种非局部低秩约束的自校准并行磁共振成像重构方法

Publications (2)

Publication Number Publication Date
CN112991483A CN112991483A (zh) 2021-06-18
CN112991483B true CN112991483B (zh) 2022-03-01

Family

ID=76340238

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110452944.2A Active CN112991483B (zh) 2021-04-26 2021-04-26 一种非局部低秩约束的自校准并行磁共振成像重构方法

Country Status (1)

Country Link
CN (1) CN112991483B (zh)

Families Citing this family (3)

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

Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101975936A (zh) * 2010-09-03 2011-02-16 杭州电子科技大学 一种基于cs压缩感知技术的快速磁共振成像方法
CN102930567A (zh) * 2012-09-25 2013-02-13 电子科技大学 多核加权最小二乘支撑向量机的磁共振并行成像重建方法
CN106019189A (zh) * 2016-05-18 2016-10-12 西南石油大学 一种基于自一致性的并行磁共振成像快速重构方法
CN109920017A (zh) * 2019-01-16 2019-06-21 昆明理工大学 基于特征向量的自一致性的联合全变分Lp伪范数的并行磁共振成像重构方法
CN109934884A (zh) * 2019-01-22 2019-06-25 昆明理工大学 一种基于变换学习和联合稀疏性的迭代自一致性并行成像重构方法
CN111754598A (zh) * 2020-06-27 2020-10-09 昆明理工大学 基于变换学习的局部空间邻域并行磁共振成像重构方法

Patent Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101975936A (zh) * 2010-09-03 2011-02-16 杭州电子科技大学 一种基于cs压缩感知技术的快速磁共振成像方法
CN102930567A (zh) * 2012-09-25 2013-02-13 电子科技大学 多核加权最小二乘支撑向量机的磁共振并行成像重建方法
CN106019189A (zh) * 2016-05-18 2016-10-12 西南石油大学 一种基于自一致性的并行磁共振成像快速重构方法
CN109920017A (zh) * 2019-01-16 2019-06-21 昆明理工大学 基于特征向量的自一致性的联合全变分Lp伪范数的并行磁共振成像重构方法
CN109934884A (zh) * 2019-01-22 2019-06-25 昆明理工大学 一种基于变换学习和联合稀疏性的迭代自一致性并行成像重构方法
CN111754598A (zh) * 2020-06-27 2020-10-09 昆明理工大学 基于变换学习的局部空间邻域并行磁共振成像重构方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
基于小波树稀疏结构的磁共振成像快速重构算法;鲍中文 等;《陕西师范大学学报》;20201012;第48卷(第6期);1-9 *

Also Published As

Publication number Publication date
CN112991483A (zh) 2021-06-18

Similar Documents

Publication Publication Date Title
CN112991483B (zh) 一种非局部低秩约束的自校准并行磁共振成像重构方法
WO2018099321A1 (zh) 一种基于广义树稀疏的权重核范数磁共振成像重建方法
CN107274462B (zh) 基于熵和几何方向的分类多字典学习磁共振图像重建方法
CN112150568A (zh) 基于Transformer模型的磁共振指纹成像重建方法
Poddar et al. Manifold recovery using kernel low-rank regularization: Application to dynamic imaging
WO2020114329A1 (zh) 磁共振快速参数成像方法及装置
CN111754598B (zh) 基于变换学习的局部空间邻域并行磁共振成像重构方法
CN112819949B (zh) 一种基于结构化低秩矩阵的磁共振指纹图像重建方法
CN105931242B (zh) 基于字典学习和时间梯度的动态核磁共振图像重建方法
CN109920017B (zh) 基于特征向量的自一致性的联合全变分Lp伪范数的并行磁共振成像重构方法
CN114820849A (zh) 基于深度学习的磁共振cest图像重建方法、装置及设备
CN109934884B (zh) 一种基于变换学习和联合稀疏性的迭代自一致性并行成像重构方法
CN112617798B (zh) 一种基于Lp范数联合全变分的并行磁共振成像重构方法
CN110148193A (zh) 基于自适应正交字典学习的动态磁共振并行重建方法
CN112581385B (zh) 基于多先验约束扩散峰度成像张量估计方法、介质和设备
Liu et al. High-fidelity MRI reconstruction using adaptive spatial attention selection and deep data consistency prior
Luo et al. Generative image priors for MRI reconstruction trained from magnitude-only images
He et al. Dynamic MRI reconstruction exploiting blind compressed sensing combined transform learning regularization
CN109188327A (zh) 基于张量积复小波紧框架的磁共振图像快速重构方法
CN115471580A (zh) 一种物理智能高清磁共振扩散成像方法
Yu et al. Universal generative modeling in dual domains for dynamic MRI
CN113628298B (zh) 基于特征向量的自一致性和非局部低秩的并行mri重构方法
CN113892938B (zh) 一种基于非局部低秩约束的改进灵敏度编码重构方法
CN114004764B (zh) 一种基于稀疏变换学习的改进灵敏度编码重建方法
Zhang et al. T2LR-Net: An unrolling network learning transformed tensor low-rank prior for dynamic MR image reconstruction

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