CN109934884B - 一种基于变换学习和联合稀疏性的迭代自一致性并行成像重构方法 - Google Patents

一种基于变换学习和联合稀疏性的迭代自一致性并行成像重构方法 Download PDF

Info

Publication number
CN109934884B
CN109934884B CN201910056437.XA CN201910056437A CN109934884B CN 109934884 B CN109934884 B CN 109934884B CN 201910056437 A CN201910056437 A CN 201910056437A CN 109934884 B CN109934884 B CN 109934884B
Authority
CN
China
Prior art keywords
representing
matrix
iteration
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
CN201910056437.XA
Other languages
English (en)
Other versions
CN109934884A (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 CN201910056437.XA priority Critical patent/CN109934884B/zh
Publication of CN109934884A publication Critical patent/CN109934884A/zh
Application granted granted Critical
Publication of CN109934884B publication Critical patent/CN109934884B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Landscapes

  • Magnetic Resonance Imaging Apparatus (AREA)

Abstract

本发明涉及一种基于变换学习和联合稀疏性的迭代自一致性并行成像重构方法，属于医学磁共振成像技术领域。本发明基于迭代自一致性并行成像重构问题，提出了一种结合变换学习和联合稀疏正则项的笛卡尔迭代自一致性并行磁共振成像重构方法，利用变量分离(Variable Splitting，VS)和交替方向乘子法(Alternating Direction Method of Multipliers，ADMM)技术进行求解。对两个实际数据集的仿真实验表明，与其他比较方法相比，本发明提出的新算法能获得更好的重构质量。

Description

一种基于变换学习和联合稀疏性的迭代自一致性并行成像重 构方法
技术领域
本发明涉及一种基于变换学习和联合稀疏性的迭代自一致性并行成像重构方法，属于医学磁共振成像技术领域。
背景技术
并行磁共振(Magnetic Resonance，MR)成像是一种众所周知的加速成像方法。其优点是通过多线圈阵列接收器接收空间灵敏度信息来减少MR成像的采样时间。在过去的二十年里，许多并行成像重构方法被提出，这些方法因灵敏度信息的使用不同方式而有所不同；如：灵敏度编码(SENSitivity Encoding，SENSE)方法是利用显式灵敏度信息来进行重构的。这一方法在实际的应用中最主要的限制是很难准确的测量灵敏度信息。另一类是一般自校准部分并行采样(GeneRalized Autocalibarating Partially ParallelAcquisition，GRAPPA)和迭代自一致性并行成像重构(Iterative Self-consistentparallel imaging reconstruction，SPIRiT)方法，它们利用隐式灵敏度信息进行重构，避免了难以准确测量灵敏度信息的限制。除此之外，还有一些不需要校准灵敏度信息的MR成像重构方法，如：无校准多线圈MR图像重构。但是相比于基于SENSE和SPIRiT的重构方法，无校准的MR图像重构方法在使用相同正则项的情况下，重构效果通常不理想。
最近，一些研究人员利用稀疏正则项模型来改善GRAPPA校准内核的重构质量。基于GRAPPA自校准方法，SPIRiT是一种在k空间网格中为每一个点(采样区和未采样区)执行校准一致性的重构方法。为了提高重构质量，一些研究者提出了联合L1(L1,2)正则项来促进多线圈图像在一些变换域中的联合稀疏性，并且在频率域或图像域中使用基于插值算子的凸集投影(Projection Over Convex Sets，POCS)方法求解SPIRiT重构问题。Weller等引入联合总变分(Joint TotalVariation，JTV)和L1,2正则项，提出了非笛卡尔SPIRiT重构方法，并通过基于先验条件的变量分离(Variable Splitting，VS)和交替方向乘子法(Alternating Direction Method of Multipliers，ADMM)进行求解。在前期工作中，提出了含JTV和L1,2正则项的笛卡尔SPIRiT重构问题的快速求解方法，即：基于SPIRiT的含联合稀疏的有效算子分裂方法(Efficient operator splitting algorithm for jointsparsity-regularized，EOSJS)。
为了进一步改善重构质量，一些研究工作者提出了基于稀疏变换(或字典)学习的重构方法。通过对参考图像进行全采样，提出了基于块的K-SVD(Singularvaluedecomposition，SVD)字典的变换学习方法。然而，在一些新图像的新特征中，从原图像获取的字典学习块通常不是有效稀疏的，从而不能有效的对图像进行重构。为此，Otaza等研究了一种在MR成像中利用K-SVD的一维(one dimensional，1D)字典的重构方法，能有效的将学习字典块进行稀疏化，但是在图像的一些二维局部结构中，1D字典法也不能使学习的字典块有效的稀疏。Ravishankar等提出了一种有效的字典学习MR成像(DictionaryLearning MR imaging，DLMRI)重构方法，它是通过从亚采样k空间方案和同时采样的MR图像中自适应学习的稀疏变换字典，该字典是利用K-SVD方法得到的。为了更有效地重构MR成像，Ravishankar等进一步提出了一种稀疏变换学习MR成像(Transform Learning MRimaging，TLMRI)重构方法。这种方法也是一种可以从训练数据中通过学习条件良好的平方稀疏变换的新重构方法。该方法可以更新学习的稀疏变换，从而得到问题的近似解。
相比于传统的固定稀疏变换的压缩感知MR成像，TLMRI方法可以有效的提高图像的质量。另外，相比于字典是利用迭代K-SVD方法求解的DLMRI方法的合成字典法，TLMRI方法是通过引入一个变换学习字典来对求出的近似解进行分析，从而能在单线圈MR成像重构中有效的节约采样时间。虽然，TLMRI方法相比于其他方法能有效的提高图像的质量和处理速度，但该方法在多线圈MR成像中的图像质量和处理速度上还有一定的改善空间。
发明内容
本发明要解决的技术问题是提供了一种基于变换学习和联合稀疏性的迭代自一致性并行成像重构方法，在于克服现有技术的不足，提高MR成像的质量。
本发明采用的技术方案是：一种基于变换学习和联合稀疏性的迭代自一致性并行成像重构方法，所述方法的具体步骤如下：
S0：初始化，令Γ-1＝((G-I)H(G-I)+μ1I)-1，x0＝AHy，
Figure BDA0001952623020000021
k＝1；
其中Γ是一个具有对角结构的矩阵，能够直接进行求逆，G表示从自校准线中获得的图像域中的k空间自一致性卷积，I表示一个(N·Nc)×(N·Nc)的单位矩阵，μ1＞0是正则项参数，(·)H表示共轭转置；
Figure BDA0001952623020000031
是待重构的所有线圈图像，x＝(x1,...,xc,...,xNc)，
Figure BDA0001952623020000032
表示x的第c列，也就是第c个线圈图像的向量化表示，N＝nx×ny表示单幅图像的像素点个数，nx和ny分别表示单幅图像的行数和列数，Nc表示线圈图像的总个数，x0表示待重构的所有线圈图像的初始值；A是系统编码矩阵，即部分傅里叶变换矩阵，且
Figure BDA0001952623020000033
是由N×N的单位矩阵IN的M行构成的亚采样矩阵，M表示单位矩阵IN中的M行，即采样数据的点数，
Figure BDA0001952623020000034
是一个N×N的矩阵，
Figure BDA0001952623020000035
分别是nx点和ny点的离散傅里叶变换矩阵，
Figure BDA00019526230200000318
表示Kronecker积；
Figure BDA0001952623020000036
是所有线圈图像的部分傅里叶测量数据，y＝[y1,...,yc,...,yNc]，
Figure BDA0001952623020000037
表示y的第c列，也就是第c个线圈图像的部分傅里叶测量数据；按照重叠步幅为1且边界环绕方式，以单线圈二维图像的每个像素为左上顶点生成N个
Figure BDA0001952623020000038
的重叠二维图像块，n表示二维图像块的总像素点数，1≤j≤N表示与重叠二维图像块相关的索引变量，
Figure BDA0001952623020000039
表示第j个
Figure BDA00019526230200000310
的二维图像块抽取矩阵，Pj由N×N的单位矩阵IN的n行构成，Pj与
Figure BDA00019526230200000311
的向量化单线圈图像相乘，则可抽取得到
Figure BDA00019526230200000312
的二维图像块的n点列向量表示；Pj与
Figure BDA00019526230200000313
的向量化多线圈图像相乘，则可抽取得到Nc个
Figure BDA00019526230200000314
的二维图像块的列向量表示组成的n×Nc矩阵，(·)T表示矩阵转置；W表示图像块的方形稀疏变换，Wk表示第k次迭代的图像块的方形稀疏变换，k为迭代次数，W0表示图像块的方形稀疏变换的初始值；由于按照重叠步幅为1且边界环绕方式将二维图像分割为重叠二维图像块，故Dc是对角矩阵，c表示线圈图像的索引变量；
Figure BDA00019526230200000315
表示
Figure BDA00019526230200000316
点离散余弦变换矩阵；
Figure BDA00019526230200000317
为辅助变量z对应的对偶变量，
Figure BDA00019526230200000319
表示uz的第c列，
Figure BDA0001952623020000041
为待重构的所有线圈图像x的辅助变量，z＝(z1,...,zc,...,zNc)，
Figure BDA0001952623020000042
表示z的第c列；
Figure BDA0001952623020000043
为辅助变量，B＝[B1,:,...,Bj,:,...,BN,:]，
Figure BDA0001952623020000044
表示B的第j个分块，
Figure BDA0001952623020000045
表示Bj,:的第c列；
Figure BDA0001952623020000046
为辅助变量B对应的对偶变量，uB＝[(uB)1,:,...,(uB)j,:,...,(uB)N,:]，
Figure BDA0001952623020000047
表示uB的第j个分块，
Figure BDA0001952623020000048
Figure BDA0001952623020000049
表示(uB)j,:的第c列；
Figure BDA00019526230200000410
为辅助变量z对应的对偶变量的初始值；
Figure BDA00019526230200000411
为辅助变量B对应的对偶变量的初始值；
S1：计算辅助变量zk，计算公式如下：
Figure BDA00019526230200000412
其中zk表示第k次迭代的待重构的所有线圈图像xk的辅助变量，xk-1表示第k-1次迭代的待重构的所有线圈图像，
Figure BDA00019526230200000413
是第k-1次迭代的辅助变量zk-1的对偶变量；vec(·)将矩阵操作数的所有列堆积成一个列向量，unvec(·)是将vec(·)转化成的列向量转化成原来形式的矩阵；
S2：初始化，j＝1；
S3：计算辅助变量的第j个分块
Figure BDA00019526230200000414
公式如下：
Figure BDA00019526230200000415
式中，Wk-1表示第k-1次迭代的图像块的方形稀疏变换，
Figure BDA00019526230200000416
表示从第k-1次迭代的待重构的所有线圈图像xk-1中抽取Nc个图像块并组成n×Nc的矩阵；
Figure BDA00019526230200000417
表示第k次迭代的辅助变量Bk的第j个分块，与
Figure BDA00019526230200000418
相对应，
Figure BDA00019526230200000419
表示第k-1次迭代的辅助变量Bk-1的对偶变量，
Figure BDA00019526230200000420
表示第k-1次迭代的对偶变量
Figure BDA00019526230200000421
的第j个分块，μ2＞0是正则项参数；HJ(b,θ)是硬阈值处理函数，其第一个输入参数
Figure BDA00019526230200000422
b＝[b1,...,bc,...,bNc]，
Figure BDA00019526230200000423
为b的第c列，第二个参数θ＞0为一个常数，HJ(b,θ)定义为：
Figure BDA00019526230200000518
其中，
Figure BDA00019526230200000519
表示矩阵的逐点乘积，
Figure BDA0001952623020000051
为一个中间变量，计算公式为：
Figure BDA0001952623020000052
repmat(s,[1,Nc])是MATLAB中的函数，功能为将第一个参数s按照第二个参数[1,Nc]的维度进行复制拓展，拓展之后得到
Figure BDA0001952623020000053
的矩阵；
S4：判断是否完成B的所有图像块分块的更新，若是，则进入步骤S5；否则，对j进行加一操作后返回S3；
S5：更新第k次迭代的图像块的方形稀疏变换Wk，计算公式如下：
Figure BDA0001952623020000054
其中：In表示n×n的单位矩阵，
Figure BDA0001952623020000055
是一个分解因子，L＝(Xk-1(Xk-1)H+0.5λIn)1/2，
Figure BDA0001952623020000056
是一个中间变量矩阵，其值为Xk-1＝[P1xk-1,...,Pjxk-1,...,PNxk-1]，xk-1表示第k-1次迭代的待重构所有线圈图像，且
Figure BDA0001952623020000057
的完整奇异值分解为UΣVH，
Figure BDA0001952623020000058
是奇异值分解因子，Bk表示第k次迭代的辅助变量，λ＞0表示正则项参数；
S6：初始化，c＝1；
S7：计算待重构的第c个线圈图像
Figure BDA0001952623020000059
公式如下：
Figure BDA00019526230200000510
其中：
Figure BDA00019526230200000511
表示第k次迭代的待重构的第c个线圈图像，
Figure BDA00019526230200000512
表示第k次迭代的辅助变量zk的第c列，
Figure BDA00019526230200000513
表示第k-1次迭代的对偶变量
Figure BDA00019526230200000514
的第c列，
Figure BDA00019526230200000515
表示第k次迭代的辅助变量Bk的第j个分块
Figure BDA00019526230200000516
的第c列，
Figure BDA00019526230200000517
表示第k-1次迭代的对偶变量
Figure BDA0001952623020000061
的第j个分块
Figure BDA0001952623020000062
的第c列，α＞0表示正则项参数；
S8：判断是否完成所有线圈，若是，进入步骤S9；否则，对c进行加一操作后返回S7；
S9：计算对偶变量
Figure BDA0001952623020000063
公式如下：
Figure BDA0001952623020000064
其中，
Figure BDA00019526230200000611
是第k次迭代的辅助变量zk对应的对偶变量，xk表示第k次迭代的待重构的所有线圈图像；
S10：计算对偶变量
Figure BDA0001952623020000065
公式如下：
Figure BDA0001952623020000066
其中，
Figure BDA0001952623020000067
表示第k次迭代的辅助变量Bk对应的对偶变量，
Figure BDA0001952623020000068
是一个中间变量矩阵，其值为Xk＝[P1xk,...,Pjxk,...,PNxk]，
Figure BDA0001952623020000069
表示第k次迭代的待重构的所有线圈图像xk中抽取Nc个图像块并组成n×Nc的矩阵；
S11：判断是否达到最大迭代次数K，K是最大迭代次数，若是，输出重构的所有线圈图像x；否则，对k进行加一操作后返回步骤S1；
S12：得到重构的所有线圈图像x后，使用平方和的平方根方法来形成幅度图像，计算公式如下：
Figure BDA00019526230200000610
本发明的有益效果是：在本发明中，提出了一种将自适应稀疏变换学习和联合稀疏正则项与笛卡尔SPIRiT并行MR成像相结合的方法，以提高重建质量。通过使用VS和ADMM技术，可以将重建问题转换为具有近似形式解的一些子问题。对两个活体数据集的实验仿真表明，所提出的方法优于其他比较的方法，能有效的提高重构图像的视觉质量。
附图说明
图1为本发明方法流程图；
图2在data1下，使用ADMMJTVLI方法重构的图像。
图3在data1下，使用DLSPIRiT方法重构的图像。
图4在data1下，使用JTLSPIRiT方法重构的图像。
图5在data1下，使用ADMMJTVLI方法重构的误差图像。
图6在data1下，使用DLSPIRiT方法重构的误差图像。
图7在data1下，使用JTLSPIRiT方法重构的误差图像。
图8在data2下，使用ADMMJTVLI方法重构的图像。
图9在data2下，使用DLSPIRiT方法重构的图像。
图10在data2下，使用JTLSPIRiT方法重构的图像。
图11在data2下，使用ADMMJTVLI方法重构的误差图像。
图12在data2下，使用DLSPIRiT方法重构的误差图像。
图13在data2下，使用JTLSPIRiT方法重构的误差图像。
具体实施方式
下面结合附图进一步详细描述本发明的技术方案：一种基于变换学习和联合稀疏性的迭代自一致性并行成像重构方法，所述方法的具体步骤如下：
假设
Figure BDA0001952623020000071
表示待重构的所有线圈图像，x＝(x1,...,xc,...,xNc)，
Figure BDA0001952623020000072
表示x的第c列，也就是第c个线圈图像的向量化表示，N＝nx×ny表示单幅图像的像素点个数，nx和ny分别表示单幅图像的行数和列数，Nc表示线圈图像的总个数。所有线圈图像的部分傅里叶测量数据或k空间亚采样数据
Figure BDA0001952623020000073
由下式给定：
y＝Ax (1)
其中，
Figure BDA0001952623020000074
Figure BDA0001952623020000075
表示y的第c列，也就是第c个线圈图像的部分傅里叶测量数据；A是一个部分傅里叶变换矩阵，且
Figure BDA0001952623020000076
Figure BDA0001952623020000077
是由N×N的单位矩阵IN的M行构成的亚采样矩阵，M表示单位矩阵IN中的M行，即采样数据的点数。
Figure BDA0001952623020000078
是一个N×N的矩阵，
Figure BDA0001952623020000079
分别是nx点和ny点的离散傅里叶变换矩阵，
Figure BDA00019526230200000710
表示Kronecker积。
SPIRiT内核一致性项是由通过在图像域中的矩阵乘积和相应k空间中的卷积块的基形成的。在图像域中的SPIRiT校准一致性方程如下所示：
vec(x)＝G·vec(x) (2)
其中G表示从自校准线中获得的图像域中的k空间自一致性卷积，vec(·)将矩阵操作数的所有列堆积成一个列向量，unvec(·)是将vec(·)转化成的列向量转化成原来形式的矩阵。
考虑到线圈图像之间小波变换域的相似性，Murphy等引入了联合稀疏来提高重建质量。因此，并行成像重建问题可以表述为以下最小化问题：
Figure BDA0001952623020000081
其中Ψ是逐线圈小波变换，且||ω||1,2定义为：
Figure BDA0001952623020000082
其中ω表示一个中间变量，ωrc表示第c个线圈的第r个k空间位置的中间变量。c是线圈图像索引变量，r是k空间位置索引变量。
Murphy等提出了一种有效的POCS方法求解问题(3)。Lustig等和Vasanawala等利用了k空间自一致性卷积，也通过使用POCS方法来解决类似问题。为了提高并行MR图像的质量，Weller等提出了利用JTV和L1,2正则项来代替问题(3)中的L1,2正则项来进行图像重构，相应的无约束问题如下：
Figure BDA0001952623020000083
其中
Figure BDA0001952623020000084
表示Frobenius范数，||·||JTV表示联合总变分正则项。I表示一个(N·Nc)×(N·Nc)的单位矩阵，α1，α2和γ均是正则化参数。Weller等提出了利用VS和ADMM方法求解问题(5)。
为了进一步改善图像重构的质量，本发明提出了一种引入多线圈图像的稀疏变换学习(Transform Learning，TL)和联合稀疏正则项模型，将重构问题转换为：
Figure BDA0001952623020000085
其中
Figure BDA0001952623020000086
为中间变量，
Figure BDA0001952623020000087
Figure BDA0001952623020000091
表示w的第c列。
Figure BDA0001952623020000092
表示线圈图像的方形稀疏变换，按照重叠步幅为1且边界环绕方式，以单线圈二维图像的每个像素为左上顶点生成N个
Figure BDA00019526230200000920
的重叠二维图像块，n表示二维图像块的总像素点数，1≤j≤N表示与重叠二维图像块相关的索引变量，
Figure BDA0001952623020000093
表示第j个
Figure BDA0001952623020000094
的二维图像块抽取矩阵，Pj由N×N的单位矩阵IN的n行构成，Pj与
Figure BDA0001952623020000095
的向量化单线圈图像相乘，则可抽取得到
Figure BDA0001952623020000096
的二维图像块的n点列向量表示；Pj与
Figure BDA0001952623020000097
的向量化多线圈图像相乘，则可抽取得到Nc个
Figure BDA0001952623020000098
的二维图像块的列向量表示组成的n×Nc矩阵，
Figure BDA0001952623020000099
表示从待重构的所有线圈图像x中抽取Nc个图像块并组成n×Nc的矩阵。函数
Figure BDA00019526230200000910
是W变换的一个正则项。与之前的变换学习的方法类似，定义
Figure BDA00019526230200000911
在目标函数中作为正则项函数，det(·)表示矩阵的行列式。α，λ均表示正则项参数。
引入辅助变量z＝x和B＝[B1,:,...,Bj,:,...,BN,:]＝[WP1x,...,WPjx,...,WPNx]，以及相应的对偶变量uz和uB。
Figure BDA00019526230200000912
表示从多线圈图像x抽取Nc个图像块组成的n×Nc的矩阵的方形稀疏变换，
Figure BDA00019526230200000913
表示B的第j个分块，
Figure BDA00019526230200000914
表示Bj,:的第c列。利用ADMM求解问题(6)，子问题的迭代解如下：
Figure BDA00019526230200000915
Figure BDA00019526230200000916
Figure BDA00019526230200000917
Figure BDA00019526230200000918
Figure BDA00019526230200000919
Figure BDA0001952623020000101
在(7)-(12)中，
Figure BDA0001952623020000102
为待重构的所有线圈图像x的辅助变量，z＝(z1,...,zc,...,zNc)，
Figure BDA0001952623020000103
表示z的第c列，zk表示第k次迭代的所有线圈图像x的辅助变量，k为迭代次数；
Figure BDA0001952623020000104
为辅助变量z对应的对偶变量，
Figure BDA0001952623020000105
表示uz的第c列，
Figure BDA0001952623020000106
是第k-1次迭代的辅助变量zk-1对应的对偶变量；xk-1表示第k-1次迭代的待重构的所有线圈图像。
Figure BDA0001952623020000107
表示第k次迭代的辅助变量Bk的第j个分块，与
Figure BDA0001952623020000108
相对应；Wk-1表示第k-1次迭代的线圈图像的方形稀疏变换；
Figure BDA0001952623020000109
表示第k-1次迭代的待重构的所有线圈图像xk-1中抽取Nc个图像块并组成n×Nc的矩阵；
Figure BDA00019526230200001010
为辅助变量B对应的对偶变量，
Figure BDA00019526230200001011
表示uB的第j个分块，
Figure BDA00019526230200001012
表示
Figure BDA00019526230200001013
的第c列，
Figure BDA00019526230200001014
表示第k-1次迭代的对偶变量
Figure BDA00019526230200001015
的第j个分块。Wk表示第k次迭代的线圈图像的方形稀疏变换。
Figure BDA00019526230200001016
是一个中间变量矩阵，其值为
Figure BDA00019526230200001017
Bk表示第k次迭代的辅助变量。
Figure BDA00019526230200001018
表示第k次迭代的待重构的第c个线圈图像，
Figure BDA00019526230200001019
表示第k次迭代的辅助变量zk的第c列，
Figure BDA00019526230200001020
表示第k-1次迭代的对偶变量
Figure BDA00019526230200001021
的第c列，
Figure BDA00019526230200001022
表示从向量化的第c个线圈图像xc中抽取一个向量化为n点的图像块，
Figure BDA00019526230200001023
表示第k次迭代的辅助变量Bk的第j个分块
Figure BDA00019526230200001024
的第c列，
Figure BDA00019526230200001025
表示第k-1次迭代的对偶变量
Figure BDA00019526230200001026
的第j个分块
Figure BDA00019526230200001027
的第c列。
Figure BDA00019526230200001028
是第k次迭代的辅助变量zk对应的对偶变量，xk表示第k次迭代的待重构所有线圈图像。
Figure BDA00019526230200001029
表示第k次迭代的辅助变量Bk对应的对偶变量，
Figure BDA00019526230200001030
是一个中间变量矩阵，其值为
Figure BDA00019526230200001031
表示第k次迭代的待重构的所有线圈图像xk中抽取Nc个图像块并组成n×Nc的矩阵。μ1，μ2均表示大于零的正则项参数。
设Γ＝(G-I)H(G-I)+μ1I，子问题(7)可以通过下式进行求解：
Figure BDA00019526230200001112
式中，Γ是一个具有对角结构的矩阵，(·)H表示共轭转置，由于SPIRiT一致性算子(G-I)的对角块结构仅在线圈之间耦合，因此可以直接对Γ中的每一个Nc×Nc块进行求逆。
问题(8)的解如下：
Figure BDA0001952623020000111
其中HJ(b,θ)是硬阈值处理函数，其第一个输入参数
Figure BDA0001952623020000112
Figure BDA0001952623020000113
为b的第c列，第二个参数θ＞0为一个常数，HJ(b,θ)定义为：
Figure BDA00019526230200001113
其中，
Figure BDA00019526230200001114
表示矩阵的逐点乘积，
Figure BDA0001952623020000114
为一个中间变量，计算公式为：
Figure BDA0001952623020000115
repmat(s,[1,Nc])是MATLAB中的函数，功能为将第一个参数s按照第二个参数[1,Nc]的维度进行复制拓展，拓展之后得到
Figure BDA0001952623020000116
的矩阵。
子问题(9)可以由下式求解：
Figure BDA0001952623020000117
其中In表示n×n的单位矩阵，
Figure BDA0001952623020000118
是一个分解因子，L＝(Xk-1(Xk-1)H+0.5λIn)1/2，且
Figure BDA0001952623020000119
具有UΣVH的全部奇异值分解(Singular value decomposition，SVD)，
Figure BDA00019526230200001110
是奇异值分解因子。
由于按照重叠步幅为1且边界环绕方式将二维图像分割为重叠二维图像块，故
Figure BDA00019526230200001111
是一个对角矩阵，其中(·)T表示矩阵转置。因此，子问题(10)能有效的求解：
Figure BDA0001952623020000121
现在，得到了子问题(7)-(10)的求解方法(如公式(13)-(17)所示)，一个基于联合稀疏和变换学习的迭代自一致性并行MR成像重建方法(Joint sparsity and TransformLearning based SPIRiT parallel MR imaging reconstruction algorithm，JTLSPIRiT)，具体流程如图1所示，其中步骤如下：
S0：初始化，令Γ-1＝((G-I)H(G-I)+μ1I)-1，x0＝AHy，
Figure BDA0001952623020000122
x0表示待重构的所有线圈图像的初始值；W0表示图像块的方形稀疏变换的初始值；
Figure BDA0001952623020000123
表示
Figure BDA0001952623020000124
离散余弦变换矩阵；
Figure BDA0001952623020000125
为辅助变量z对应的对偶变量的初始值；
Figure BDA0001952623020000126
为辅助变量B对应的对偶变量的初始值。
S1：计算辅助变量zk，计算公式如(13)；
S2：初始化，j＝1；
S3：计算辅助变量的第j个分块
Figure BDA0001952623020000127
公式如(14)；
S4：判断是否完成B的所有图像块分块更新，若是，则进入步骤S5；否则，对j进行加一操作后返回S3；
S5：更新第k次迭代的图像块的方形稀疏变换Wk，计算公式如(16)；
S6：初始化，c＝1；
S7：计算待重构的第c个线圈图像
Figure BDA0001952623020000128
公式如(17)；
S8：判断是否完成所有线圈，若是，进入步骤S9；否则，对c进行加一操作后返回S7；
S9：计算对偶变量
Figure BDA0001952623020000129
公式如(11)；
S10：计算对偶变量
Figure BDA0001952623020000131
公式如(12)；
S11：判断是否达到最大迭代次数K，K是最大迭代次数，若是，输出重构的所有线圈图像x；否则，对k进行加一操作后返回步骤S1；
S12：得到重构的所有线圈图像x后，使用平方和的平方根方法来形成幅度图像，计算公式如下：
Figure BDA0001952623020000132
在以下实验中，比较了所提出的方法JTLSPlRiT与使用JTV和L1正则项的SPIRiT并行成像重构方法(ADMMJTVL1)和基于K-SVD字典学习的SPIRiT并行成像重构方法(DLSPIRiT)。所有比较的方法都在Matlab中实现。
本发明使用了两个由八通道线圈获得的完整视场(Field of View，FOV)的体内大脑数据集，即data1和data2。为了生成测试数据集，采用加速因子为R×(不包括自校准采样(Self-calibrated sampling，ACS)线)和24×24ACS线的笛卡尔泊松盘亚采样模式对完全采样的数据集进行人工亚采样。此外，在所有比较方法的实验中使用5×5SPIRiT校准内核，并且从完全采样数据集中使用公式(18)计算参考图像。
对于DLSPIRiT，使用大小为6×6(n＝36)的图像块，并且利用从14400个随机选择块中用K-SVD迭代方法(10次迭代)来学习1.5倍的超完备合成词典。对于JTLSPIRiT，使用6×6(n＝36)的图像块。手动调整方法的正则化参数以获得更好的重建质量。当xk和xk-1之间的相对误差低于公差tol＝2e-4时，所有比较的方法都被终止。本发明所有的实验都是在配备Intel Core i7-7500U@2.7GHz处理器，8GB内存和Windows 10操作系统(64位)的笔记本电脑上进行的。
图2-13展示了采用加速因子为6×和24×24ACS线的泊松盘亚采样模式，在数据集data1和data2下所有比较方法的并行MR成像的重构图像及其相应的误差图。在图2中，可以明显看出ADMMJTVL1方法的重构的图像出现了较多的伪影。在图3中，DLSPIRiT方法可以稍微缓解图像重构的伪影，但在过平滑的区域中，其重构的图像存在一些细节上的缺失。相比之下，在图4中，提出的JTLSPIRiT方法有效的抑制了图像重构中的伪像，并且保留了其他方法中缺失的细节。图5-7分别是ADMMJTVLI，DLSPIRiT和JTLSPIRiT的重建误差图像，从图7中可以看出，提出的JTLSPIRiT重构图像的误差算子较小(黑色的部分越多，误差算子越小)，能够重构出图像的更多细节。因此JTLSPIRiT在所有比较方法中能得到最佳的视觉质量。对于图8-13，也可以得出同样的结论。
在实验仿真中，还对所有比较的方法和本发明提出的方法的运行时间进行了比较。对于data1，ADMMJTVL1,DLSPIRiT和JTLSPIRiT的平均运行时间分别为20秒，1600秒和260秒。对于data2，ADMMJTVL1,DLSPIRiT和JTLSPIRiT的平均运行时间分别为68秒，1950秒和410秒。综上可以看出，与ADMMJTVL1相比，提出的JTLSPIRiT算法尽管需要更多的重构时间，但重构质量有较大的提高。而与DLSPIRiT相比，JTLSPIRiT重构图像的SNR有所增加，视觉质量有较大改善，而且重构速度是DLSPIRiT的5倍以上。
以上结合附图对本发明的具体实施方式作了详细说明，但是本发明并不限于上述实施方式，在本领域普通技术人员所具备的知识范围内，还可以在不脱离本发明宗旨的前提下作出各种变化。

Claims (1)

1.一种基于变换学习和联合稀疏性的迭代自一致性并行成像重构方法，其特征在于：具体步骤如下：
S0：初始化，令Γ-1＝((G-I)H(G-I)+μ1I)-1，x0＝AHy，
Figure FDA0001952623010000011
k＝1；其中Γ是一个具有对角结构的矩阵，能够直接进行求逆，G表示从自校准线中获得的图像域中的k空间自一致性卷积，I表示一个(N·Nc)×(N·Nc)的单位矩阵，μ1＞0是正则项参数，(·)H表示共轭转置；
Figure FDA0001952623010000012
是待重构的所有线圈图像，
Figure FDA00019526230100000118
Figure FDA0001952623010000013
表示x的第c列，也就是第c个线圈图像的向量化表示，N＝nx×ny表示单幅图像的像素点个数，nx和ny分别表示单幅图像的行数和列数，Nc表示线圈图像的总个数，x0表示待重构的所有线圈图像的初始值；A是系统编码矩阵，即部分傅里叶变换矩阵，且
Figure FDA0001952623010000014
是由N×N的单位矩阵IN的M行构成的亚采样矩阵，M表示单位矩阵IN中的M行，即采样数据的点数，
Figure FDA0001952623010000015
是一个N×N的矩阵，
Figure FDA0001952623010000016
分别是nx点和ny点的离散傅里叶变换矩阵，
Figure FDA0001952623010000017
表示Kronecker积；
Figure FDA0001952623010000018
是所有线圈图像的部分傅里叶测量数据，
Figure FDA0001952623010000019
Figure FDA00019526230100000110
表示y的第c列，也就是第c个线圈图像的部分傅里叶测量数据；按照重叠步幅为1且边界环绕方式，以单线圈二维图像的每个像素为左上顶点生成N个
Figure FDA00019526230100000111
的重叠二维图像块，n表示二维图像块的总像素点数，1≤j≤N表示与重叠二维图像块相关的索引变量，
Figure FDA00019526230100000112
表示第j个
Figure FDA00019526230100000113
的二维图像块抽取矩阵，Pj由N×N的单位矩阵IN的n行构成，Pj与
Figure FDA00019526230100000114
的向量化单线圈图像相乘，则可抽取得到
Figure FDA00019526230100000115
的二维图像块的n点列向量表示；Pj与
Figure FDA00019526230100000116
的向量化多线圈图像相乘，则可抽取得到Nc个
Figure FDA00019526230100000117
的二维图像块的列向量表示组成的n×Nc矩阵，(·)T表示矩阵转置；W表示图像块的方形稀疏变换，Wk表示第k次迭代的图像块的方形稀疏变换，k为迭代次数，W0表示图像块的方形稀疏变换的初始值；由于按照重叠步幅为1且边界环绕方式将二维图像分割为重叠二维图像块，故Dc是对角矩阵，c表示线圈图像的索引变量；
Figure FDA0001952623010000021
表示
Figure FDA0001952623010000022
点离散余弦变换矩阵；
Figure FDA0001952623010000023
为辅助变量z对应的对偶变量，
Figure FDA0001952623010000024
Figure FDA0001952623010000025
表示uz的第c列，
Figure FDA0001952623010000026
为待重构的所有线圈图像x的辅助变量，
Figure FDA0001952623010000027
Figure FDA0001952623010000028
表示z的第c列；
Figure FDA0001952623010000029
为辅助变量，
Figure FDA00019526230100000210
Figure FDA00019526230100000211
表示B的第j个分块，
Figure FDA00019526230100000212
Figure FDA00019526230100000213
表示Bj,:的第c列；
Figure FDA00019526230100000214
为辅助变量B对应的对偶变量，uB＝[(uB)1,:,...,(uB)j,:,...,(uB)N,:]，
Figure FDA00019526230100000215
表示uB的第j个分块，
Figure FDA00019526230100000216
Figure FDA00019526230100000217
表示(uB)j,:的第c列；
Figure FDA00019526230100000218
为辅助变量z对应的对偶变量的初始值；
Figure FDA00019526230100000219
为辅助变量B对应的对偶变量的初始值；
S1：计算辅助变量zk，计算公式如下：
Figure FDA00019526230100000220
其中zk表示第k次迭代的待重构的所有线圈图像xk的辅助变量，xk-1表示第k-1次迭代的待重构的所有线圈图像，
Figure FDA00019526230100000221
是第k-1次迭代的辅助变量zk-1的对偶变量；vec(·)将矩阵操作数的所有列堆积成一个列向量，unvec(·)是将vec(·)转化成的列向量转化成原来形式的矩阵；
S2：初始化，j＝1；
S3：计算辅助变量的第j个分块
Figure FDA00019526230100000222
公式如下：
Figure FDA00019526230100000223
式中，Wk-1表示第k-1次迭代的图像块的方形稀疏变换，
Figure FDA00019526230100000224
表示从第k-1次迭代的待重构的所有线圈图像xk-1中抽取Nc个图像块并组成n×Nc的矩阵；
Figure FDA0001952623010000031
表示第k次迭代的辅助变量Bk的第j个分块，与Wk-1Pjxk-1相对应，
Figure FDA0001952623010000032
表示第k-1次迭代的辅助变量Bk-1的对偶变量，
Figure FDA0001952623010000033
表示第k-1次迭代的对偶变量
Figure FDA0001952623010000034
的第j个分块，μ2＞0是正则项参数；HJ(b,θ)是硬阈值处理函数，其第一个输入参数
Figure FDA0001952623010000035
Figure FDA0001952623010000036
Figure FDA0001952623010000037
为b的第c列，第二个参数θ＞0为一个常数，HJ(b,θ)定义为：
Figure FDA0001952623010000038
其中，
Figure FDA0001952623010000039
表示矩阵的逐点乘积，
Figure FDA00019526230100000310
为一个中间变量，计算公式为：
Figure FDA00019526230100000311
repmat(s,[1,Nc])是MATLAB中的函数，功能为将第一个参数s按照第二个参数[1,Nc]的维度进行复制拓展，拓展之后得到
Figure FDA00019526230100000312
的矩阵；
S4：判断是否完成B的所有图像块分块的更新，若是，则进入步骤S5；否则，对j进行加一操作后返回S3；
S5：更新第k次迭代的图像块的方形稀疏变换Wk，计算公式如下：
Figure FDA00019526230100000313
其中：In表示n×n的单位矩阵，
Figure FDA00019526230100000314
是一个分解因子，L＝(Xk-1(Xk-1)H+0.5λIn)1/2，
Figure FDA00019526230100000315
是一个中间变量矩阵，其值为Xk-1＝[P1xk-1,...,Pjxk-1,...,PNxk-1]，xk-1表示第k-1次迭代的待重构所有线圈图像，且
Figure FDA00019526230100000316
的完整奇异值分解为UΣVH，U,Σ,
Figure FDA00019526230100000317
是奇异值分解因子，Bk表示第k次迭代的辅助变量，λ＞0表示正则项参数；
S6：初始化，c＝1；
S7：计算待重构的第c个线圈图像
Figure FDA00019526230100000318
公式如下：
Figure FDA0001952623010000041
其中：
Figure FDA0001952623010000042
表示第k次迭代的待重构的第c个线圈图像，
Figure FDA0001952623010000043
表示第k次迭代的辅助变量zk的第c列，
Figure FDA0001952623010000044
表示第k-1次迭代的对偶变量
Figure FDA0001952623010000045
的第c列，
Figure FDA0001952623010000046
表示第k次迭代的辅助变量Bk的第j个分块
Figure FDA0001952623010000047
的第c列，
Figure FDA0001952623010000048
表示第k-1次迭代的对偶变量
Figure FDA0001952623010000049
的第j个分块
Figure FDA00019526230100000410
的第c列，α＞0表示正则项参数；
S8：判断是否完成所有线圈，若是，进入步骤S9；否则，对c进行加一操作后返回S7；
S9：计算对偶变量
Figure FDA00019526230100000411
公式如下：
Figure FDA00019526230100000412
其中，
Figure FDA00019526230100000413
是第k次迭代的辅助变量zk对应的对偶变量，xk表示第k次迭代的待重构的所有线圈图像；
S10：计算对偶变量
Figure FDA00019526230100000414
公式如下：
Figure FDA00019526230100000415
其中，
Figure FDA00019526230100000416
表示第k次迭代的辅助变量Bk对应的对偶变量，
Figure FDA00019526230100000417
是一个中间变量矩阵，其值为
Figure FDA00019526230100000418
Figure FDA00019526230100000419
表示第k次迭代的待重构的所有线圈图像xk中抽取Nc个图像块并组成n×Nc的矩阵；
S11：判断是否达到最大迭代次数K，K是最大迭代次数，若是，输出重构的所有线圈图像x；否则，对k进行加一操作后返回步骤S1；
S12：得到重构的所有线圈图像x后，使用平方和的平方根方法来形成幅度图像，计算公式如下：
Figure FDA00019526230100000420
CN201910056437.XA 2019-01-22 2019-01-22 一种基于变换学习和联合稀疏性的迭代自一致性并行成像重构方法 Active CN109934884B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201910056437.XA CN109934884B (zh) 2019-01-22 2019-01-22 一种基于变换学习和联合稀疏性的迭代自一致性并行成像重构方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201910056437.XA CN109934884B (zh) 2019-01-22 2019-01-22 一种基于变换学习和联合稀疏性的迭代自一致性并行成像重构方法

Publications (2)

Publication Number Publication Date
CN109934884A CN109934884A (zh) 2019-06-25
CN109934884B true CN109934884B (zh) 2022-05-24

Family

ID=66985005

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201910056437.XA Active CN109934884B (zh) 2019-01-22 2019-01-22 一种基于变换学习和联合稀疏性的迭代自一致性并行成像重构方法

Country Status (1)

Country Link
CN (1) CN109934884B (zh)

Families Citing this family (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111754598B (zh) * 2020-06-27 2022-02-25 昆明理工大学 基于变换学习的局部空间邻域并行磁共振成像重构方法
CN112834971A (zh) * 2020-12-31 2021-05-25 苏州朗润医疗系统有限公司 一种基于奇异值分解的mri迭代自校准并行成像算法
CN112991483B (zh) * 2021-04-26 2022-03-01 昆明理工大学 一种非局部低秩约束的自校准并行磁共振成像重构方法

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105184755A (zh) * 2015-10-16 2015-12-23 西南石油大学 基于自一致性的含联合全变分的并行磁共振成像高质量重构方法
CN106019189A (zh) * 2016-05-18 2016-10-12 西南石油大学 一种基于自一致性的并行磁共振成像快速重构方法

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US9734601B2 (en) * 2014-04-04 2017-08-15 The Board Of Trustees Of The University Of Illinois Highly accelerated imaging and image reconstruction using adaptive sparsifying transforms

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105184755A (zh) * 2015-10-16 2015-12-23 西南石油大学 基于自一致性的含联合全变分的并行磁共振成像高质量重构方法
CN106019189A (zh) * 2016-05-18 2016-10-12 西南石油大学 一种基于自一致性的并行磁共振成像快速重构方法

Non-Patent Citations (7)

* Cited by examiner, † Cited by third party
Title
Transductive Zero-Shot Learning with a Self-training dictionary approach;Yunlong Yu等;《IEEE ACCESS》;20170527;第1-10页 *
基于字典学习和加权TV的MRI重构算法;宋长明等;《电子技术应用》;20171231(第01期);第148-151页 *
基于自一致性的磁共振并行成像高效重构算法;段继忠等;《天津大学学报(自然科学与工程技术版)》;20140531;第47卷(第5期);第414-419页 *
新一代视频编码技术HEVC算法分析及比较;刘昱等;《电视技术》;20121017(第20期);第51-55页 *
约束恒模医学CT图像盲均衡算法;孙云山等;《计算机应用》;20110601(第06期);第129-131页 *
采用稀疏变换和拉普拉斯金字塔的数字水印算法;顾宇鑫等;《计算机辅助设计与图形学学报》;20180515(第05期);第157-166页 *
面向目标检测的稀疏表示方法研究进展;高仕博等;《电子学报》;20150215(第02期);第114-126页 *

Also Published As

Publication number Publication date
CN109934884A (zh) 2019-06-25

Similar Documents

Publication Publication Date Title
Wang et al. DIMENSION: dynamic MR imaging with both k‐space and spatial prior knowledge obtained via multi‐supervised network training
Schlemper et al. A deep cascade of convolutional neural networks for dynamic MR image reconstruction
CN104933683B (zh) 一种用于磁共振快速成像的非凸低秩重建方法
Ravishankar et al. Data-driven learning of a union of sparsifying transforms model for blind compressed sensing
Loktyushin et al. Blind multirigid retrospective motion correction of MR images
CN109934884B (zh) 一种基于变换学习和联合稀疏性的迭代自一致性并行成像重构方法
CN111754598B (zh) 基于变换学习的局部空间邻域并行磁共振成像重构方法
CN107274462A (zh) 基于熵和几何方向的分类多字典学习磁共振图像重建方法
CN112991483B (zh) 一种非局部低秩约束的自校准并行磁共振成像重构方法
KR102584166B1 (ko) 리스케일링과 인공신경망을 적용한 자기 공명 영상 처리 장치 및 그 방법
Kelkar et al. Prior image-constrained reconstruction using style-based generative models
CN109920017B (zh) 基于特征向量的自一致性的联合全变分Lp伪范数的并行磁共振成像重构方法
Fan et al. An interpretable MRI reconstruction network with two-grid-cycle correction and geometric prior distillation
Ke et al. CRDN: cascaded residual dense networks for dynamic MR imaging with edge-enhanced loss constraint
Zhang et al. T2LR-Net: An unrolling network learning transformed tensor low-rank prior for dynamic MR image reconstruction
Guo et al. A Joint Group Sparsity-based deep learning for multi-contrast MRI reconstruction
CN112617798A (zh) 一种基于Lp范数联合全变分的并行磁共振成像重构方法
CN114004764B (zh) 一种基于稀疏变换学习的改进灵敏度编码重建方法
He et al. Dynamic MRI reconstruction exploiting blind compressed sensing combined transform learning regularization
Liu et al. Highly undersampling dynamic cardiac MRI based on low-rank tensor coding
Pan et al. Iterative self-consistent parallel magnetic resonance imaging reconstruction based on nonlocal low-rank regularization
CN112634385B (zh) 一种基于深度拉普拉斯网络的快速磁共振成像方法
CN113674379A (zh) 一种基于参考支撑集的共稀疏分析模型的mri重构方法、系统及计算机可读存储介质
Ke et al. Deep low-rank prior in dynamic MR imaging
Huang et al. Orthogonal tensor dictionary learning for accelerated dynamic MRI

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