CN110148215B - 一种基于平滑约束和局部低秩约束模型的四维磁共振图像重建方法 - Google Patents

一种基于平滑约束和局部低秩约束模型的四维磁共振图像重建方法 Download PDF

Info

Publication number
CN110148215B
CN110148215B CN201910432660.XA CN201910432660A CN110148215B CN 110148215 B CN110148215 B CN 110148215B CN 201910432660 A CN201910432660 A CN 201910432660A CN 110148215 B CN110148215 B CN 110148215B
Authority
CN
China
Prior art keywords
image
magnetic resonance
iteration
dimensional
constraint
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
CN201910432660.XA
Other languages
English (en)
Other versions
CN110148215A (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.)
Harbin Institute of Technology
Original Assignee
Harbin Institute of 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 Harbin Institute of Technology filed Critical Harbin Institute of Technology
Priority to CN201910432660.XA priority Critical patent/CN110148215B/zh
Publication of CN110148215A publication Critical patent/CN110148215A/zh
Application granted granted Critical
Publication of CN110148215B publication Critical patent/CN110148215B/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
    • G06T17/00—Three dimensional [3D] modelling, e.g. data description of 3D objects
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/10—Image acquisition modality
    • G06T2207/10072—Tomographic images
    • G06T2207/10088—Magnetic resonance imaging [MRI]
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
    • Y02A90/00—Technologies having an indirect contribution to adaptation to climate change
    • Y02A90/30—Assessment of water resources

Landscapes

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

Abstract

一种基于平滑约束和局部低秩约束模型的四维磁共振图像重建方法，属于四维磁共振图像重建领域，为了解决现有重建方法存为的重建图像伪影严重以及重建时间较长的问题。包括：将欠采样K空间数据进行三维逆傅里叶变换得到初始重建图像；建立图像重建模型；引入辅助变量和正则化参数；获得三个子问题；初始化内、外循环最大迭代次数和收敛阈值；在内循环中迭代求解重建图像；判断重建图像是否满足收敛条件，满足则退出内循环，否则返回内循环继续迭代；未达到外循环最大迭代次数则更新正则化系数后再进入内循环迭代，否则终止，获得最终重建图像。实验表明本发明能有效加速磁共振成像且能较好地恢复图像细节，并获得具有临床诊断意义的重建图像。

Description

一种基于平滑约束和局部低秩约束模型的四维磁共振图像重 建方法
技术领域
本发明属于磁共振成像领域的四维磁共振图像重建领域，具体涉及一种基于平滑约束和局部低秩约束模型的四维磁共振图像重建方法。
背景技术
四维磁共振成像技术(Four-dimensional Magnetic Resonance Imaging,4DMRI)可以同时展现受检组织器官的三维解剖结构信息和动态变化信息。例如，四维心脏磁共振(4D CMR)成像可提供可视化三维心脏解剖结构和心脏瓣膜运动信息，为心血管疾病的临床诊断以及房室结构和心肌功能评价提供帮助；四维放射治疗(4D Radiation Therapy,4DRT)借助四维磁共振图像可实现对肿瘤的四维跟踪，避免对健康组织有所损害。
四维磁共振数据采集时需要获取包含几乎连续时刻上的三维数据以致所需扫描时间非常长，在长时间扫描中，患者不自主的生理运动导致图像产生运动伪影，给临床诊断带来困难。同时对一些患心脏疾病或哮喘病的患者及儿童患者，无法进行长时间扫描。所以加速成像过程成为了四维磁共振成像中亟待解决的问题。加速四维磁共振成像的关键一方面在于对硬件条件的改善和快速扫描序列的研究，一方面在于图像重建方法的研究。受磁共振扫描仪器物理特性及患者生理情况的限制，前者加速磁共振成像的发展空间有限，所以对图像重建算法的研究是解决问题的关键。目前大部分重建算法都仅在空间域或时间域进行约束，并未充分利用数据中的先验信息，因此在重建图像质量等方面存在一定缺陷。
发明内容
本发明为了解决现有重建方法图像细节的重建效果差以及重建时间较长的问题，充分挖掘了四维磁共振图像在时空域上的先验信息，提出了一种基于平滑约束和局部低秩约束模型的四维磁共振图像重建方法，达到了高质量重建图像细节和加速四维磁共振成像的目的。
为实现上述目的，本发明所采取的技术方案是：
一种基于平滑约束和局部低秩约束模型的四维磁共振图像重建方法，所述四维包括三维空间域和一维时间域；所述方法包含以下步骤：
(1)利用预先设置的欠采样模板获取参考四维磁共振图像的欠采样K空间数据，即测量数据b；
(2)对测量数据b进行三维逆傅里叶变换将其变换到图像域，获得初始重建图像X(0)＝F-1b，F表示傅里叶变换操作；
(3)结合四维磁共振图像的平滑特性和局部化低秩特性建立如下图像重建模型：
Figure SMS_1
式中，X*用于表示最终的重建图像；X表示重建四维磁共振图像；A＝SF表示欠采样傅里叶变换；第一项
Figure SMS_2
为数据保真项，b表示测量数据；第二项
Figure SMS_3
表示对空间域及时间域的平滑约束，采用l1范数，其中D(·)表示四维有限差分算子；第三项是对P组矩阵的低秩约束，||·||*表示核范数；其中Pi为提取所有时刻上第i个小图像块的操作算子，再利用提取到的所有时刻上的小图像块构建局部化低秩矩阵，用
Figure SMS_4
表示分组构建局部化低秩矩阵的操作算子；λ1、λ2分别表示平衡第二项、第三项的约束项系数；
(4)引入辅助变量T和Wi，将式(1)表示的重建模型改写为如下易于求解的形式：
Figure SMS_5
利用拉格朗日乘子法，将其写为如下无约束形式：
Figure SMS_6
其中β1和β2为正则化参数；
(5)本发明采用内外双重循环的模式；初始化外循环最大迭代次数，内循环最大的迭代次数和收敛阈值；将上述式(3)所示问题分为X、T和Wi三个子问题，在内循环中采用交替的方法迭代求解X、T和Wi三个子问题，在外循环中更新β1和β2；
得到第k次迭代结果为：
Figure SMS_7
Figure SMS_8
Figure SMS_9
(6)判断本次迭代重建图像X(k)是否满足收敛条件或是否达到最大迭代次数，未达到最大迭代次数且不满足收敛条件则返回步骤(5)进行k+1次迭代求解，否则退出内循环；
(7)判断是否达到最大外循环迭代次数，达到最大迭代次数则算法终止，未达到则更新β1和β2，重置内循环迭代次数为k＝1且将本次重建图像作为X(0)，返回步骤(5)。
进一步地，在步骤(1)中，所述参考四维磁共振图像大小为：D1×D2×D3×Nt；D1、D2、D3分别表示三个空间维度尺寸，Nt表示时间维的大小；
所述欠采样模板为径向SoS采样模板(采样模板采用如图2所示的SoS方式)，该采样模板在kx-ky平面进行径向采样，在kz方向进行笛卡尔采样；数据采集过程中先采集所有kx-ky平面内给定旋转角度下的辐条，然后增加旋转角度再进行采集，增加的角度为111.25°(黄金角度)，重复此过程；所需采集的辐条数量根据欠采样倍数确定；kx、ky、kz表示笛卡尔坐系中K空间的三个坐标轴。
进一步地，在步骤(2)中，只在空间维度上进行逆傅里叶变换。
进一步地，在步骤(3)中，在建立的基于平滑约束和局部低秩约束的四维磁共振图像重建模型中，平滑约束项中D(X)定义为如下加权梯度的形式：
Figure SMS_10
其中，由于四维磁共振图像在时间域和空间域上的平滑性不同，且表现为时间域的平滑性强于空间域，所以引入了参数γ增强了在时间维度上的惩罚，γ>1；
x、y、z、t表示在笛卡尔坐标系下的四维图像坐标。
进一步地，在步骤(3)中(按照如图3所示的方法构建局部化低秩矩阵)，
小图像块的划分方式：设置空间域每个分块尺寸为M×N×Q和扫描步长S1,S2,S3，其中1<S1<M,1<S2<N,1<S3<Q，将每一时刻图像按相同的分块尺寸和扫描步长划分为若干个部分重叠的小图像块，同时获得权重系数矩阵Weight，权重系数表示某像素被划分到图像块数量，由于对三维图像的分块方式是相同的，所以每时刻的分块数量P也相同的；将所有时刻对应位置上大小为M×N×Q小图像块视为一组，共P组，
Figure SMS_11
分组构建局部化低秩矩阵的操作算子
Figure SMS_12
具体为：将每组内的三维图像块列化并按时间维顺序排列成大小为M·N·Q×t矩阵，获得P组低秩矩阵；由于局部像素之间的相关性，所以每组矩阵都具有低秩性。
进一步地，在步骤(5)中，为了使式(3)与式(1)所示问题等价，在外循环中，β1与β2的值需要不断增大，达到106数量级。
进一步地，在步骤(5)中，X、T和Wi三个子问题解的获得过程为：
通过保持其中两个变量不变，更新另一个变量的方式来交替求解式(3)，因此包括以下三个步骤：
步骤5-1，固定X和Wi，求解T，得到如下T子问题：
Figure SMS_13
利用软阈值法求解得：
Figure SMS_14
步骤5-2，固定X和T，求解Wi，得到如下Wi子问题:
Figure SMS_15
上式可以看作是对P组损失函数分别进行最小化，所以将式(10)等价于：
Figure SMS_16
考虑到每小块在计算上是独立的，且对
Figure SMS_17
进行奇异值分解有：
Figure SMS_18
则上式用软阈值法求解，得：
Figure SMS_19
步骤5-3，固定Wi和T，求解Wi，得到如下X子问题：
Figure SMS_20
分析上式的第三项，其是对各分块矩阵差值的平方和，所以可以利用权重矩阵Weight将局部化矩阵重构为原四维图像矩阵以便求解，即
Figure SMS_21
则式(13)等价于如下形式：
Figure SMS_22
显然上式是关于X的二次多项式，对上式进行求导并令导数等于零得方程式(15)，满足方程的解即为X子问题的解；
(2ATA+λ1β1DTD+λ2β2)X＝2ATb+λ1β1DT(T)+λ2β2W (15)
同时对等式两边进行傅里叶变换，再简单移项后，可得求解X的一步计算解析式如下：
Figure SMS_23
由式(9)得到迭代形式的公式(4)，由式(12)得到迭代形式的公式(5)，由式(16)得到迭代形式的公式(6)。
进一步地，步骤(6)中包含以下操作步骤：
步骤6-1，当k小于最大迭代次数且不满足迭代收敛条件时，将本次迭代的解作为第k+1次迭代的输入，返回步骤(5)继续迭代；
步骤6-2，当连续两次迭代的相对误差小于允许范围时，算法收敛，退出内循环，收敛条件为：
|ε(k-1)-ε(k)|<0.005ε(k) (17)
Figure SMS_24
ε表示重建误差。
进一步地，步骤(7)中包含以下操作步骤：
步骤7-1，判断是否达到最大外循环迭代次数，是则退出外循环，终止算法，当前X即为最终的重建图像，否则采用连续策略更新正则化参数β1和β2，如式(18)和式(19)所示：
β1＝β1·βinc1 (18)
β2＝β2·βinc2 (19)
且为了获得与原问题式(1)相同的解，在外循环中，β1与β2的值需要不断增大，达到106数量级，βinc1、βinc2分别表示β1、β2的更新倍数；
步骤7-2，重置内循环迭代次数为k＝1且将本次重建图像作为内循环初始迭代图像X(0)，并返回步骤(5)。
本发明的有益效果是：
本发明研究在K空间数据高倍欠采样的情况下，利用少量采集数据重建高分辨率四维磁共振图像的方法。
本发明充分挖掘四维磁共振图像在时空域上的先验信息，建立并求解重建模型，包括以下步骤：利用预先设置的欠采样模板获取参考四维磁共振图像的欠采样K空间数据；将欠采样四维K空间数据在时间帧上分别进行三维逆傅里叶变换得到初始重建图像；构建局部化低秩矩阵，建立图像重建模型；引入辅助变量和正则化参数；将问题分为三个子问题，采用内外双重循坏模式，在内循环中迭代求解获得重建图像；判断重建图像是否满足收敛条件及是否达到最大内循环迭代次数，不满足收敛条件并且未达到最大内循环迭代次数则返回内循环继续迭代，否则退出内循环，进入外循环；未达到外循环最大迭代次数则更新正则化系数，并将本次重建图像作为内循环迭代初始图像再进入内循环迭代，否则算法终止并获得最终重建图像。
实验表明本发明具有能有效加速磁共振成像、较好地重建图像细节及能较好地克服成像伪影的优点，具有临床诊断意义。所提出的算法计算效率较高，CPU运行时间较短，是一种高效的算法；同时利用本发明提出的算法，可以较好地重建磁共振图像，较完整地保留图像中的细节区域，克服成像伪影。实验结果表明，与传统方法相比，本发明能够在较高加速倍数下获得高质量重建图像，同时提高了重建速度，实现了加速磁共振图像重建的目的。实验表明本发明能有效加速磁共振成像且能较好地恢复图像细节，并获得具有临床诊断意义的重建图像。
本发明的技术要点为：
(1)利用预先设置的欠采样模板获取参考四维磁共振图像的欠采样K空间数据；
(2)将欠采样K空间四维数据在每一时间帧上进行三维逆傅里叶变换得到初始重建图像；
(3)设置分块尺寸和扫描步长，按相同的分块尺寸和扫描步长将四维磁共振图像每一时间帧上的三维图像划分为若干个部分重叠的图像块，同时获得权重系数矩阵，并将图像块构建成局部化低秩矩阵；结合图像空间域和时间域平滑先验信息和局部化低秩特性建立图像重建模型；
(4)引入辅助变量和正则化参数β1和β2；
(5)将问题分为易求解的三个子问题；初始化外循环最大迭代次数、内循环的最大迭代次数和收敛阈值；采用内外双重循环模式，在内循环中迭代求解获得重建图像；
(6)判断重建图像是否满足收敛条件及是否达到内循环最大迭代次数，满足收敛条件或达到内循环最大迭代次数则退出内循环，进入外循环，否则返回步骤(5)继续迭代；
(7)判断是否达到外循环最大迭代次数，若达到则获得最终重建图像，否则在外循环更新正则化系数β1和β2后，将本次重建图像作为内循环迭代初始图像返回步骤(5)。
附图说明
图1为本发明方法流程图(基于平滑约束和局部低秩约束模型的四维磁共振图像重建方法)；
图2为径向SoS采样模板示意图；
图3为局部化低秩矩阵构建方法示意图；
图4为6倍欠采样下重建效果图；
图5为8倍欠采样下重建效果图；
图6为12倍欠采样下重建效果图；
图7为16倍欠采样下重建效果图。
具体实施方法
下面结合附图1至7对本发明进行详细说明。
如图1所示，本发明的具体实施方式如下：
(1)利用预先设置的欠采样模板获取参考四维磁共振图像的欠采样K空间数据，即测量数据，用b表示；
(2)对测量数据b进行三维逆傅里叶变换将其变换到图像域，获得初始重建图像X(0)＝F-1b，F表示傅里叶变换操作；
(3)利用四维磁共振图像空间域、时间域平滑特性和局部化低秩特性建立如下图像重建模型；
Figure SMS_25
式中，X*用于表示最终的重建图像；X表示重建四维磁共振图像；A＝SF表示欠采样傅里叶变换；第一项
Figure SMS_26
为数据保真项，b表示测量数据；第二项
Figure SMS_27
表示对空间域及时间域的平滑约束，采用l1范数，其中D(·)表示四维有限差分算子；第三项是对P组矩阵的低秩约束，||·||*表示核范数；其中Pi为提取所有时刻上第i个小图像块的操作算子，再利用提取到的所有时刻上的小图像块构建局部化低秩矩阵，用
Figure SMS_28
表示分组构建局部化低秩矩阵的操作算子。λ1、λ2分别表示平衡第二项、第三项的约束项系数；
(4)引入辅助变量T和Wi，将式(1)变为如下易求解方程：
Figure SMS_29
利用拉格朗日乘子法，将其写为如下无约束形式：
Figure SMS_30
其中β1和β2为正则化参数。
(5)初始化外循环最大迭代次数Initer，内循环最大的迭代次数Outiter和收敛阈值ε；在内循环中对变量X，T和Wi进行交替求解，求得第k次迭代结果为：
Figure SMS_31
Figure SMS_32
Figure SMS_33
式(5)中Ui (k-1)与Vi (k-1)表示对
Figure SMS_34
进行奇异值分解所得的奇异值向量，如：
Figure SMS_35
(6)判断本次迭代重建图像X是否满足收敛条件或是否达到最大迭代次数，未达到最大迭代次数且不满足收敛条件则返回步骤(5)进行k+1次迭代求解，否则退出内循环，进入外循环；
(7)判断是否达到最大外循环迭代次数，达到最大迭代次数则算法终止，未达到则更新β1和β2，重置内循环迭代次数为k＝1且将本次重建图像作为X(0)，返回步骤(5)。
步骤(1)中欠采样模板使用如图2所示的径向SoS采样模板，该采样模板在kx-ky平面进行径向采样，在kz方向进行笛卡尔采样。数据采集过程中先采集所有kx-ky平面内给定旋转角度下的辐条，然后增加旋转角度再进行采集，增加的角度为111.25°(黄金角度)，重复此过程。所需采集的辐条数量根据欠采样倍数确定，最终获得的采集数据记为测量数据b，欠采样倍数ACC为全采样K空间数据量与欠采样K空间数据量的比值。
步骤(1)中，参考四维磁共振图像大小为D1×D2×D3×Nt(3D空间域+1D时间域)步骤(3)具体包括：
步骤3-1，设置分块尺寸M×N×Q和扫描步长S1,S2,S3，其中1<S1<M,1<S2<N,1<S3<Q，将各时刻上的三维图像划分成若干重叠的图像块，同时获得权重系数矩阵Weight，权重系数表示某像素被划分到的图像块数量；
步骤3-2，按照如图3所示的方法构建局部化低秩矩阵。所有时刻对应位置上的大小为M×N×Q图像块视为一组，一共P组，
Figure SMS_36
再将每组内的三维图像块列化并按时间维顺序排列成大小为M·N·Q×T矩阵，得到P组矩阵。由于局部像素之间的相关性，所以每组矩阵都具有低秩性；
步骤3-3，建立基于平滑约束和局部低秩约束的四维磁共振图像重建模型
步骤(3)平滑约束项中D(X)定义为如下加权梯度的形式：
Figure SMS_37
本发明中，步骤(4)引入了辅助变量T和Wi，获得了易求解的目标函数，并引入正则化参数β1和β2将其写成了无约束的形式。
本发明采用内外双重循环模式，在步骤(5)中初始化外循环的最大迭代次数，内循环的最大迭代次数和收敛阈值，并在内循环中迭代交替求解X，T和Wi。具体是通过保持其中两个变量不变，更新另一个变量的方式来交替求解，求解过程具体有以下三个步骤：
步骤5-1，固定X和Wi，求解T，得到如下T子问题:
Figure SMS_38
利用软阈值法求解得：
Figure SMS_39
步骤5-2，固定X和T，求解Wi，得到如下Wi子问题:
Figure SMS_40
上式可以看作是对P组损失函数分别进行最小化，所以将式(10)等价于：
Figure SMS_41
考虑到每小块在计算上是独立的，且对
Figure SMS_42
进行奇异值分解有：
Figure SMS_43
则上式用软阈值法求解，得：
Figure SMS_44
步骤5-3，固定Wi和T，求解Wi，得到如下X子问题:
Figure SMS_45
分析上式的第三项，其是对各分块矩阵差值的平方和，所以利用权重系数矩阵Weight将局部化矩阵重构为原四维图像矩阵以便求解，即
Figure SMS_46
则上式等价于如下形式：
Figure SMS_47
显然上式是关于X的二次多项式，对上式进行求导并令导数等于零得方程式(15)，满足方程的解即为X子问题的解。
(2ATA+λ1β1DTD+λ2β2)X＝2ATb+λ1β1DT(T)+λ2β2W (14)
同时对等式两边进行傅里叶变换，再简单移项后，可得求解X的一步计算解析式如下：
Figure SMS_48
步骤(6)具体包括：
步骤6-1，当k小于最大迭代次数且不满足迭代收敛条件时，将本次迭代的解作为第k+1次迭代的输入，返回步骤(5)继续迭代。
步骤6-2，当连续两次迭代的相对误差小于允许范围时，退出内循环。收敛条件为：
|ε(k-1)-ε(k)|<0.005ε(k) (16)
Figure SMS_49
ε表示重建误差
步骤(7)中具体包括：
步骤7-1，判断是否达到最大外循环迭代次数，是则终止算法，当前X即为最终的重建图像，否则采用连续策略更新正则化参数β1和β2，如式(18)和式(19)所示：
β1＝β1gβinc1 (17)
β2＝β2gβinc2 (18)
为了获得与原问题式(1)相同的解，在外循环中，β1与β2的值需要不断增大，最终达到106数量级。
步骤7-2，重置内循环迭代次数为k＝1且将本次重建图像作为内循环初始迭代图像X(0)，返回步骤(5)。
本发明将通过以下实验进行有效性的说明
1.实验环境：
实验室用台式机参数：CPU为Inter(R)Core(TM)i5-4590，主频为3.30GHz，内存为16G，操作系统为Win10 64位系统，实验软件为MatlabR2016b。
2.实验结果及结果分析
本发明实验所用的图像数据为四维腹部磁共振图像，实验的欠采样倍数分别为6，8，12，16倍。为了验证本方法的有效性，本发明将利用信号误差比(SER)与结构性相似度(SSIM)对重建效果进行定量评估，并与四维全变分方法(4D TV)进行对比。同时，为了说明本方法的可行性，本方法重建图像将与初始重建图像进行比较。
选择一幅四维腹部磁共振图像作为参考图像，大小为：D1×D2×D3×Nt，将其每个像素值减去最小像素值后，再除以最大与最小像素值之差，得到归一化到[0，1]后的四维参考磁共振图像。之后按照如图1所示流程图进行重建实验。
本发明进行了欠采样倍数分别为6，8，12，16倍的四组实验，对实验结果进行定量评估如表1及表2所示。重建图像如图4到图7所示。
选取四维图像一帧的一个二维层面图像进行实验效果展示，图4展示了6倍欠采样下重建效果，图4(a)-(c)依次为初始重建图像，4D TV方法重建图像，本发明重建图像；图4(d)-(f)为重建图像与参考图像的差值图，依次为初始重建图像差值图，4D TV方法重建图像差值图，本发明重建图像差值图。图5展示了8倍欠采样下重建效果，图5(a)-(c)依次为初始重建图像，4D TV方法重建图像，本发明重建图像；图5(d)-(f)为重建图像与参考图像的差值图，依次为初始重建图像差值图，4D TV方法重建图像差值图，本发明重建图像差值图。图6展示了12倍欠采样下重建效果，图6(a)-(c)依次为初始重建图像，4D TV方法重建图像，本发明重建图像；图6(d)-(f)为重建图像与参考图像的差值图，依次为初始重建图像差值图，4D TV方法重建图像差值图，本发明重建图像差值图。图7展示了16倍欠采样下重建效果，图7(a)-(c)依次为初始重建图像，4D TV方法重建图像，本发明重建图像；图7(d)-(f)为重建图像与参考图像的差值图，依次为初始重建图像差值图，4D TV方法重建图像差值图，本发明重建图像差值图。
实例一：本实例欠采样倍数为6倍，初始重建图像SER：15.74dB；SSIM：0.8137；4DTV重建图像SER：24.93dB；SSIM：0.9902；本发明重建图像SER：30.20dB，SSIM：0.9952；
实例二：本实例欠采样倍数为8倍，如表2所示初始重建图像SER：14.00dB；SSIM：0.7645；4D TV重建图像SER：22.94dB；SSIM：0.9854；本发明重建图像SER：27.88dB；SSIM：0.9922；
实例三：本实例欠采样倍数为12倍，初始重建图像SER：11.45dB；SSIM：0.6758；4DTV重建图像SER：20.52dB；SSIM：0.9737；本发明重建图像SER：25.92dB；SSIM：0.9898；
实例四：本实例欠采样倍数为16倍，初始重建图像SER：9.85dB；SSIM：0.6284；4DTV重建图像SER：18.95dB；SSIM：0.9665；本发明重建图像SER：24.12dB；SSIM：0.9859。表1重建图像信号误差比(dB)
Figure SMS_50
表2重建图像结构性相似度
Figure SMS_51
实验结果表明在相同的欠采样倍数下，利用本发明中提出的方法重建出的四维磁共振图像具有较高的信号误差比与结构性相似度，并且从重建图像上看，本发明方法对细节的重建效果好。相比于四维全变分方法，本发明所提方法可以在高倍欠采样的情况下获得更好的重建结果，重建图像细节更为清晰，与参考图像最为接近。

Claims (6)

1.一种基于平滑约束和局部低秩约束模型的四维磁共振图像重建方法，所述四维包括三维空间域和一维时间域；其特征在于，所述方法包含以下步骤：
(1)利用预先设置的欠采样模板获取参考四维磁共振图像的欠采样K空间数据，即测量数据b；
(2)对测量数据b进行三维逆傅里叶变换将其变换到图像域，获得初始重建图像X(0)＝F-1b，F表示傅里叶变换操作；
(3)结合四维磁共振图像的平滑特性和局部化低秩特性建立如下图像重建模型：
Figure FDA0004092862360000011
式中，X*用于表示最终的重建图像；X表示重建四维磁共振图像；A＝SF表示欠采样傅里叶变换；第一项
Figure FDA0004092862360000012
为数据保真项，b表示测量数据；第二项
Figure FDA0004092862360000013
表示对空间域及时间域的平滑约束，采用1范数，其中D(·)表示四维有限差分算子；第三项是对P组矩阵的低秩约束，||·||*表示核范数；其中Pi为提取所有时刻上第i个小图像块的操作算子，再利用提取到的所有时刻上的小图像块构建局部化低秩矩阵，用
Figure FDA0004092862360000014
表示分组构建局部化低秩矩阵的操作算子；λ1、λ2分别表示平衡第二项、第三项的约束项系数；
(4)引入辅助变量T和Wi，将式(1)表示的重建模型改写为如下易于求解的形式：
Figure FDA0004092862360000015
利用拉格朗日乘子法，将其写为如下无约束形式：
Figure FDA0004092862360000016
其中β1和β2为正则化参数；
(5)采用内外双重循环的模式；初始化外循环最大迭代次数，内循环最大的迭代次数和收敛阈值；将上述式(3)所示问题分为X、T和Wi三个子问题，在内循环中采用交替的方法迭代求解X、T和Wi三个子问题，在外循环中更新β1和β2；
得到第k次迭代结果为：
Figure FDA0004092862360000021
Figure FDA0004092862360000022
Figure FDA0004092862360000023
(6)判断本次迭代重建图像X(k)是否满足收敛条件或是否达到最大迭代次数，未达到最大迭代次数且不满足收敛条件则返回步骤(5)进行k+1次迭代求解，否则退出内循环；
(7)判断是否达到最大外循环迭代次数，达到最大迭代次数则算法终止，未达到则更新β1和β2，重置内循环迭代次数为k＝1且将本次重建图像作为X(0)，返回步骤(5)；
在步骤(1)中，所述参考四维磁共振图像大小为：D1×D2×D3×Nt；
D1、D2、D3分别表示三个空间维度尺寸，Nt表示时间维的大小；
所述欠采样模板为径向SoS采样模板，该采样模板在kx-ky平面进行径向采样，在kz方向进行笛卡尔采样；数据采集过程中先采集所有kx-ky平面内给定旋转角度下的辐条，然后增加旋转角度再进行采集，增加的角度为111.25°，重复此过程；所需采集的辐条数量根据欠采样倍数确定；kx、ky、kz表示笛卡尔坐系中K空间的三个坐标轴；
在步骤(3)中，在建立的基于平滑约束和局部低秩约束的四维磁共振图像重建模型中，平滑约束项中D(X)定义为如下加权梯度的形式：
Figure FDA0004092862360000024
其中，由于四维磁共振图像在时间域和空间域上的平滑性不同，且表现为时间域的平滑性强于空间域，所以引入了参数γ增强了在时间维度上的惩罚，γ>1；
x、y、z、t表示在笛卡尔坐标系下的四维图像坐标；
在步骤(3)中，小图像块的划分方式：设置空间域每个分块尺寸为M×N×Q和扫描步长S1,S2,S3，其中1<S1<M,1<S2<N,1<S3<Q，将每一时刻图像按相同的分块尺寸和扫描步长划分为若干个部分重叠的小图像块，同时获得权重系数矩阵Weight，权重系数表示某像素被划分到图像块数量，由于对三维图像的分块方式是相同的，所以每时刻的分块数量P也相同的；将所有时刻对应位置上大小为M×N×Q小图像块视为一组，共P组，
Figure FDA0004092862360000031
分组构建局部化低秩矩阵的操作算子
Figure FDA0004092862360000032
具体为：将每组内的三维图像块列化并按时间维顺序排列成大小为M·N·Q×t矩阵，获得P组低秩矩阵。
2.根据权利要求1所述的一种基于平滑约束和局部低秩约束模型的四维磁共振图像重建方法，其特征在于，在步骤(2)中，只在空间维度上进行逆傅里叶变换。
3.根据权利要求2所述的一种基于平滑约束和局部低秩约束模型的四维磁共振图像重建方法，其特征在于，在步骤(5)中，为了使式(3)与式(1)所示问题等价，在外循环中，β1与β2的值需要不断增大，达到106数量级。
4.根据权利要求3所述的一种基于平滑约束和局部低秩约束模型的四维磁共振图像重建方法，其特征在于，在步骤(5)中，X、T和Wi三个子问题解的获得过程为：
通过保持其中两个变量不变，更新另一个变量的方式来交替求解式(3)，因此包括以下三个步骤：
步骤5-1，固定X和Wi，求解T，得到如下T子问题：
Figure FDA0004092862360000033
利用软阈值法求解得：
Figure FDA0004092862360000034
步骤5-2，固定X和T，求解Wi，得到如下Wi子问题:
Figure FDA0004092862360000035
上式可以看作是对P组损失函数分别进行最小化，所以将式(10)等价于：
Figure FDA0004092862360000036
考虑到每小块在计算上是独立的，且对
Figure FDA0004092862360000039
进行奇异值分解有：
Figure FDA0004092862360000037
则上式用软阈值法求解，得：
Figure FDA0004092862360000038
步骤5-3，固定Wi和T，求解Wi，得到如下X子问题：
Figure FDA0004092862360000041
分析上式的第三项，其是对各分块矩阵差值的平方和，所以可以利用权重矩阵Weight将局部化矩阵重构为原四维图像矩阵以便求解，即
Figure FDA0004092862360000042
则式(13)等价于如下形式：
Figure FDA0004092862360000043
显然上式是关于X的二次多项式，对上式进行求导并令导数等于零得方程式(15)，满足方程的解即为X子问题的解；
(2ATA+λ1β1DTD+λ2β2)X＝2ATb+λ1β1DT(T)+λ2β2W (15)
同时对等式两边进行傅里叶变换，再简单移项后，可得求解X的一步计算解析式如下：
Figure FDA0004092862360000044
由式(9)得到迭代形式的公式(4)，由式(12)得到迭代形式的公式(5)，由式(16)得到迭代形式的公式(6)。
5.根据权利要求4所述的一种基于平滑约束和局部低秩约束模型的四维磁共振图像重建方法，其特征在于，
步骤(6)中包含以下操作步骤：
步骤6-1，当k小于最大迭代次数且不满足迭代收敛条件时，将本次迭代的解作为第k+1次迭代的输入，返回步骤(5)继续迭代；
步骤6-2，当连续两次迭代的相对误差小于允许范围时，算法收敛，退出内循环，收敛条件为：
|ε(k-1)-ε(k)|＜0.005ε(k) (17)
Figure FDA0004092862360000045
ε表示重建误差。
6.根据权利要求5所述的一种基于平滑约束和局部低秩约束模型的四维磁共振图像重建方法，其特征在于，
步骤(7)中包含以下操作步骤：
步骤7-1，判断是否达到最大外循环迭代次数，是则退出外循环，终止算法，当前X即为最终的重建图像，否则采用连续策略更新正则化参数β1和β2，如式(18)和式(19)所示：
β1＝β1·βinc1 (18)
β2＝β2·βinc2 (19)
且为了获得与原问题式(1)相同的解，在外循环中，β1与β2的值需要不断增大，达到106数量级，βinc1、βinc2分别表示β1、β2的更新倍数；
步骤7-2，重置内循环迭代次数为k＝1且将本次重建图像作为内循环初始迭代图像X(0)，并返回步骤(5)。
CN201910432660.XA 2019-05-22 2019-05-22 一种基于平滑约束和局部低秩约束模型的四维磁共振图像重建方法 Active CN110148215B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201910432660.XA CN110148215B (zh) 2019-05-22 2019-05-22 一种基于平滑约束和局部低秩约束模型的四维磁共振图像重建方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201910432660.XA CN110148215B (zh) 2019-05-22 2019-05-22 一种基于平滑约束和局部低秩约束模型的四维磁共振图像重建方法

Publications (2)

Publication Number Publication Date
CN110148215A CN110148215A (zh) 2019-08-20
CN110148215B true CN110148215B (zh) 2023-05-19

Family

ID=67592715

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201910432660.XA Active CN110148215B (zh) 2019-05-22 2019-05-22 一种基于平滑约束和局部低秩约束模型的四维磁共振图像重建方法

Country Status (1)

Country Link
CN (1) CN110148215B (zh)

Families Citing this family (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110490832B (zh) * 2019-08-23 2023-05-05 哈尔滨工业大学 一种基于正则化深度图像先验方法的磁共振图像重建方法
CN112927313B (zh) * 2019-12-05 2022-12-06 上海联影医疗科技股份有限公司 磁共振图像重建方法、装置、计算机设备及可读存储介质
CN111275669B (zh) * 2020-01-13 2022-04-22 西安交通大学 一种先验信息引导的四维锥束ct图像重建算法
CN112037298A (zh) * 2020-08-20 2020-12-04 上海联影医疗科技股份有限公司 图像重建方法、装置、计算机设备和存储介质
CN112213673B (zh) * 2020-09-07 2022-11-22 上海东软医疗科技有限公司 动态磁共振成像方法、装置、重建计算机及磁共振系统
CN112489155B (zh) * 2020-12-07 2024-01-26 中国科学院深圳先进技术研究院 图像重建方法和装置、电子设备以及机器可读存储介质
CN112819949B (zh) * 2021-02-07 2024-03-26 哈尔滨工业大学 一种基于结构化低秩矩阵的磁共振指纹图像重建方法
CN113763501B (zh) * 2021-09-08 2024-02-27 上海壁仞智能科技有限公司 图像重建模型的迭代方法和图像重建方法
CN113866694B (zh) * 2021-09-26 2022-12-09 上海交通大学 一种快速三维磁共振t1定量成像方法、系统及介质

Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103985111A (zh) * 2014-02-21 2014-08-13 西安电子科技大学 一种基于双字典学习的4d-mri超分辨率重构方法
WO2014162300A1 (en) * 2013-04-05 2014-10-09 Isis Innovation Ltd. Acceleration of low-rank mri data acquisition
CN106530258A (zh) * 2016-11-22 2017-03-22 哈尔滨工业大学 基于高阶全变分正则化的快速迭代磁共振图像重建方法
CN106934778A (zh) * 2017-03-10 2017-07-07 北京工业大学 一种基于小波域结构和非局部分组稀疏的mr图像重建方法
CN107991636A (zh) * 2017-11-24 2018-05-04 哈尔滨工业大学 一种基于适应性结构低秩矩阵的快速磁共振图像重建方法
WO2018137199A1 (en) * 2017-01-25 2018-08-02 Tsinghua University Real-time phase-contrast flow mri with low rank modeling and parallel imaging
CN109615675A (zh) * 2018-12-04 2019-04-12 厦门大学 一种多通道磁共振成像的图像重建方法

Family Cites Families (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
EP2059907A2 (en) * 2006-08-29 2009-05-20 Koninklijke Philips Electronics N.V. Reduction of heart motion artifacts in thoracic ct imaging
US9964621B2 (en) * 2013-10-01 2018-05-08 Beth Israel Deaconess Medical Center, Inc. Methods and apparatus for reducing scan time of phase contrast MRI

Patent Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2014162300A1 (en) * 2013-04-05 2014-10-09 Isis Innovation Ltd. Acceleration of low-rank mri data acquisition
CN103985111A (zh) * 2014-02-21 2014-08-13 西安电子科技大学 一种基于双字典学习的4d-mri超分辨率重构方法
CN106530258A (zh) * 2016-11-22 2017-03-22 哈尔滨工业大学 基于高阶全变分正则化的快速迭代磁共振图像重建方法
WO2018137199A1 (en) * 2017-01-25 2018-08-02 Tsinghua University Real-time phase-contrast flow mri with low rank modeling and parallel imaging
CN106934778A (zh) * 2017-03-10 2017-07-07 北京工业大学 一种基于小波域结构和非局部分组稀疏的mr图像重建方法
CN107991636A (zh) * 2017-11-24 2018-05-04 哈尔滨工业大学 一种基于适应性结构低秩矩阵的快速磁共振图像重建方法
CN109615675A (zh) * 2018-12-04 2019-04-12 厦门大学 一种多通道磁共振成像的图像重建方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
"ADAPTING REGULARIZED LOW-RANK MODELS FOR PARALLEL ARCHITECTURES";Driggs Derek等;《SIAM JOURNAL ON SCIENTIFIC COMPUTING》;20190315;第41卷(第1期);第A163-A189页 *
基于向量稀疏和矩阵低秩的压缩感知核磁共振图像重建算法;张红雨;《天津理工大学学报》;20170215;第33卷(第01期);第28-32页 *

Also Published As

Publication number Publication date
CN110148215A (zh) 2019-08-20

Similar Documents

Publication Publication Date Title
CN110148215B (zh) 一种基于平滑约束和局部低秩约束模型的四维磁共振图像重建方法
Qin et al. Convolutional recurrent neural networks for dynamic MR image reconstruction
Armanious et al. Unsupervised medical image translation using cycle-MedGAN
CN111123183B (zh) 基于复数R2U_Net网络的快速磁共振成像方法
Banerjee et al. A completely automated pipeline for 3D reconstruction of human heart from 2D cine magnetic resonance slices
CN109360152A (zh) 基于稠密卷积神经网络的三维医学图像超分辨率重建方法
CN111127521B (zh) 用于生成和跟踪目标的形状的系统和方法
CN107330953B (zh) 一种基于非凸低秩的动态mri重建方法
WO2022183988A1 (en) Systems and methods for magnetic resonance image reconstruction with denoising
CN111754598B (zh) 基于变换学习的局部空间邻域并行磁共振成像重构方法
CN111487573B (zh) 一种用于磁共振欠采样成像的强化型残差级联网络模型
CN111598964B (zh) 一种基于空间自适应网络的定量磁化率图像重建方法
Bai et al. Super-resolution reconstruction of MR brain images
CN109741439B (zh) 一种二维mri胎儿图像的三维重建方法
Barbano et al. Steerable conditional diffusion for out-of-distribution adaptation in imaging inverse problems
CN112617798B (zh) 一种基于Lp范数联合全变分的并行磁共振成像重构方法
Yi et al. Fast and Calibrationless low-rank parallel imaging reconstruction through unrolled deep learning estimation of multi-channel spatial support maps
Gu et al. Cross-modality image translation: CT image synthesis of MR brain images using multi generative network with perceptual supervision
Lv et al. Reconstruction of undersampled radial free‐breathing 3D abdominal MRI using stacked convolutional auto‐encoders
CN117495992A (zh) 一种欠采样核磁共振图像的重建方法
Qiao et al. CorGAN: Context aware recurrent generative adversarial network for medical image generation
CN114004764B (zh) 一种基于稀疏变换学习的改进灵敏度编码重建方法
CN105488757A (zh) 一种脑纤维稀疏重建的方法
CN113567901B (zh) 一种磁共振旋转坐标系下的自旋晶格弛豫成像方法和系统
Xie et al. Non‐invasive reconstruction of dynamic myocardial transmembrane potential with graph‐based total variation constraints

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
CB03 Change of inventor or designer information

Inventor after: Hu Yue

Inventor after: Li Peng

Inventor after: Lin Disi

Inventor after: Zhao Kuangshi

Inventor before: Hu Yue

Inventor before: Lin Disi

Inventor before: Zhao Kuangshi

CB03 Change of inventor or designer information
GR01 Patent grant
GR01 Patent grant