CN113298907B - 一种基于伽马核范数和总变分的核磁图像重建方法 - Google Patents
一种基于伽马核范数和总变分的核磁图像重建方法 Download PDFInfo
- Publication number
- CN113298907B CN113298907B CN202110691450.XA CN202110691450A CN113298907B CN 113298907 B CN113298907 B CN 113298907B CN 202110691450 A CN202110691450 A CN 202110691450A CN 113298907 B CN113298907 B CN 113298907B
- Authority
- CN
- China
- Prior art keywords
- nuclear magnetic
- magnetic image
- matrix
- value
- 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
Links
- 238000000034 method Methods 0.000 title claims abstract description 52
- 238000005070 sampling Methods 0.000 claims abstract description 73
- 238000012545 processing Methods 0.000 claims abstract description 6
- 239000011159 matrix material Substances 0.000 claims description 139
- 230000003190 augmentative effect Effects 0.000 claims description 22
- 238000004364 calculation method Methods 0.000 claims description 13
- 238000005457 optimization Methods 0.000 claims description 13
- 238000000354 decomposition reaction Methods 0.000 claims description 7
- 230000017105 transposition Effects 0.000 claims description 6
- OAICVXFJPJFONN-UHFFFAOYSA-N Phosphorus Chemical compound [P] OAICVXFJPJFONN-UHFFFAOYSA-N 0.000 claims description 5
- 230000008569 process Effects 0.000 claims description 5
- 238000006243 chemical reaction Methods 0.000 claims description 3
- 230000005251 gamma ray Effects 0.000 claims description 3
- 230000005389 magnetism Effects 0.000 claims description 3
- 230000001131 transforming effect Effects 0.000 claims description 3
- 235000009508 confectionery Nutrition 0.000 claims description 2
- 230000006870 function Effects 0.000 description 39
- 230000009286 beneficial effect Effects 0.000 description 5
- 238000003384 imaging method Methods 0.000 description 4
- 238000002595 magnetic resonance imaging Methods 0.000 description 4
- 230000000747 cardiac effect Effects 0.000 description 2
- 230000007547 defect Effects 0.000 description 2
- 238000013461 design Methods 0.000 description 2
- 238000002059 diagnostic imaging Methods 0.000 description 2
- 230000000694 effects Effects 0.000 description 2
- 238000012360 testing method Methods 0.000 description 2
- 238000004458 analytical method Methods 0.000 description 1
- 210000003484 anatomy Anatomy 0.000 description 1
- 238000003759 clinical diagnosis Methods 0.000 description 1
- 238000004195 computer-aided diagnosis Methods 0.000 description 1
- 238000003745 diagnosis Methods 0.000 description 1
- 201000010099 disease Diseases 0.000 description 1
- 208000037265 diseases, disorders, signs and symptoms Diseases 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 230000005865 ionizing radiation Effects 0.000 description 1
- 210000003127 knee Anatomy 0.000 description 1
- 238000011031 large-scale manufacturing process Methods 0.000 description 1
- 238000005259 measurement Methods 0.000 description 1
- 230000007246 mechanism Effects 0.000 description 1
- 238000013421 nuclear magnetic resonance imaging Methods 0.000 description 1
- 210000000056 organ Anatomy 0.000 description 1
- 230000001575 pathological effect Effects 0.000 description 1
- 230000010412 perfusion Effects 0.000 description 1
- 230000010287 polarization Effects 0.000 description 1
- 238000011084 recovery Methods 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 238000004088 simulation Methods 0.000 description 1
- 238000001356 surgical procedure Methods 0.000 description 1
- 230000004083 survival effect Effects 0.000 description 1
- 238000009423 ventilation Methods 0.000 description 1
- 230000000007 visual effect Effects 0.000 description 1
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
本发明公开了一种基于伽马核范数和总变分的核磁图像重建方法，包括以下步骤：S1、对原始的核磁图像帧进行下采样处理，得到下采样核磁图像帧；S2、基于下采样核磁图像帧，建立动态核磁图像重建模型；S3、对动态核磁图像重建模型进行求解，得到重构后的动态核磁图像；本发明解决了L1正则项模型并不能取得最稀疏的结果，只能获得次优解的问题。
Description
技术领域
本发明涉及图像处理技术领域，具体涉及一种基于伽马核范数和总变分的核磁图像重建方法。
背景技术
核磁共振成像采用一种非侵入式的成像模式，对人体没有相应的电离辐射和损害，其成像能够有效地反映人体的生理信息和分子生物信息特征。更重要的是核磁共振成像能够获取人体器官在时域维度上的解剖结构以及生理特征上的连续影像，在医院的临床诊断具有重大的参考作用，尤其在视觉引导手术、计算机辅助诊断以及临床病理机理研究等领域。核磁共振成像虽然具备以上优点，但也存在一些不足，如扫描的速度较慢、成像质量不佳等。扫描数据的缓慢速度一方面影响了成像的视觉质量，导致核磁图像中出现混叠伪影、噪声等，低质量的核磁图像很难辅助医生的疾病诊断；另一方面，扫描的时间过长，对病人的要求相应的提高，即要求病人长时间的保持不动，给病人增加了痛苦，同时也使其规模化的生产应用受到一定的局限。
近年来，基于压缩感知的核磁共振成像(Compressed Sensing based MagneticResonance Imaging,CS-MRI)技术被引入到核磁图像重建中，显著地提升了成像速度。CS只需要采样少数的k空间(也就是傅里叶域)数据，然后从中恢复出完整的磁共振(MR)图像。当前的CS-MRI主要包括离线CS-MRI和在线CS-MRI两类方法。离线CS-MRI方法一般在对信号进行重建前需要先要获得整个图像序列的下采样数据。这一类重建方法主要利用的是帧与帧之间的纹理相关性和整个序列在时空域上的结构稀疏性，从而重建得到整个序列影像。不同于离线重建方法，在线CS-MRI方法一般基于单帧重建，即对整个动态核磁影像序列的恢复是一帧一帧得到的。这类重建方法只利用本帧图像的观测数据及已经重建得到的图像来恢复当前帧。即在线CS-MRI方法可以概括为两步，即第1步是利用已重建好的图像来得到预测帧，第2步是对残差图像进行恢复。由于核磁图像帧间存在较强的纹理相关性，可以利用该特性来估计得到预测帧，而残差帧和参考帧的恢复则要利用整体的稀疏性及所提出的重建模型。
虽然CS-MRI重建方法具有较好的优点，但使用CS-MRI重建方法获得的图像质量将受到两个因素的影响。其中一个因素是感兴趣信号的稀疏性或者在某种表示下的稀疏性对重建性能的影响，另一个因素是采样矩阵的约束等距性常数。在实际应用中，重建方法需要求解一个非线性凸优化问题。当核磁图像数据量较大时，重建的效率将较低。另外，基于以上分析，在重建方法的设计中，大部分正则化方法基于L1范数进行求解。实际上，L1正则项模并不能取得最稀疏的结果，只能获得次优解。比较基于L1范数的正则项模型和基于L0范数的正则项模型，采用L1正则项模型去逼近L0正则项模型，存在一定的偏差，主要是因为L0范数反映的是非零元的个数，而L1范数约束解的绝对值个数不是解的非零元素数目。
现有技术中基于压缩感知的核磁图像重建方法存在的问题：实际上，L1正则项模型并不能取得最稀疏的结果，只能获得次优解。
发明内容
针对现有技术中的上述不足，本发明提供的一种基于伽马核范数和总变分的核磁图像重建方法解决了L1正则项模型并不能取得最稀疏的结果，只能获得次优解的问题。
为了达到上述发明目的，本发明采用的技术方案为：一种基于伽马核范数和总变分的核磁图像重建方法，包括以下步骤：
S1、对原始的核磁图像帧进行下采样处理，得到下采样核磁图像帧；
S2、基于下采样核磁图像帧，建立动态核磁图像重建模型；
S3、对动态核磁图像重建模型进行求解，得到重构后的动态核磁图像。
进一步地，步骤S1中进行下采样处理的公式为：
bt＝UFxt+nt
其中，bt为第t采样时刻的下采样核磁图像帧；U为下采样算子，F为2维离散傅里叶变换矩阵，此时，n为矩阵F的列，m为矩阵F的行，xt为第t采样时刻的原始的核磁图像帧，nt为第t采样时刻的噪声，除了矩阵F外，bt，U，xt，nt中的m均为矩阵行大小，n为矩阵列大小，
为复数域空间。
进一步地，步骤S2中动态核磁图像重建模型的公式为：
B＝[b1,b2,...,bT]H，t＝1,2,...,T，
其中，argmin(·)为对等式右边的目标函数求最小值，将目标函数取得最小值时的极小点赋给等式左边的变量，
为重构后的动态核磁图像，
为矩阵的F范数的平方，tr(·)为求一个矩阵的迹，U为下采样算子，F为2维离散傅里叶变换，X为需要重构的动态核磁图像，||X||γ为需要重构的动态核磁图像X的伽马核范数，Vec(·)为向量化算子，B为由T个下采样核磁图像帧构成的动态下采样核磁图像，T为总的采样时刻数量，bt为第t采样时刻的下采样核磁图像帧，H为转置运算，λ1和λ2为正则化参数，γ为伽马参数，min{P,T}为取P和T的最小值，P＝mn，P为需要重构的动态核磁图像X的行大小，σi(X)表示需要重构的动态核磁图像X的第i个奇异值，||X||TV为需要重构的动态核磁图像X的总变分，
为复数域空间。
上述进一步方案的有益效果为：本发明利用加码范数较好地刻画了核磁图像的全局信息，从而可较大程度地去除核磁图像中的伪影，进一步使得重建后的核磁图像更为清晰。
进一步地，需要重构的动态核磁图像X的总变分||X||TV的计算过程为：
A1、计算每个采样时刻的核磁图像帧的总变分：
其中，||xt||TV为第t采样时刻的核磁图像帧的总变分，xi,j,t表示第t采样时刻下原始核磁图像帧xt在(i,j)位置处的像素值，xi+1,j,t表示第t采样时刻下原始核磁图像帧xt在(i+1,j)位置处的像素值；i＝1,2,...,m—1，j＝1,2,...,n—1表示像素所在的位置下标；
A2、计算T个采样时刻的核磁图像帧的总变分的总和，得到需要重构的动态核磁图像X的总变分||X||TV：
其中，Xt为需要重构的动态核磁图像X的第t个向量，O为一个算子，该算子将需要重构的动态核磁图像X的第t个向量变换为相应的二维核磁图像帧。
上述进一步方案的有益效果为：用总变分方法能充分描述核磁图像帧的局部一致性，这样能较好地保持重构后核磁图像的细节信息。
进一步地，步骤S3包括以下分步骤：
S31、将动态核磁图像重建模型表示为优化问题：
满足M＝X,Y＝M
其中，min(·)为最小化算子，即求括号里目标函数的最优值，U为下采样算子，F为2维离散傅里叶变换，X为需要重构的动态核磁图像，B为由T个下采样核磁图像帧构成的动态下采样核磁图像，H为转置运算，λ1和λ2为正则化参数，γ为伽马参数，σi(X)表示需要重构的动态核磁图像X的第i个奇异值，||X||TV为需要重构的动态核磁图像X的总变分，min{P,T}为取P和T的最小值，P＝mn，m为矩阵行大小，n为矩阵列大小，
为复数域空间，P为需要重构的动态核磁图像X的行大小，tr(·)为求一个矩阵的迹，M、Y均为与需要重构的动态核磁图像X等价的辅助矩阵，||X||γ为重构后的动态核磁图像X的伽马核范数，
为矩阵的F范数的平方；
S32、将优化问题构建为增广拉格朗日函数：
<Λ1,X—M>＝tr(Λ1 H(X-M))
<Λ2,M—Y>＝tr(Λ2 H(X-Y))
其中，
为最小化算子，具体为求关于优化变量X,Y,M,Λ1,Λ2下整个括号里的目标函数的最优值，
为矩阵的F范数的平方，||M||γ为M的伽马核范数，||Y||TV为Y的总变分，Λ1、Λ2为矩阵形式的拉格朗日乘子，<Λ1,X—M>表示求Λ1与(X—M)的内积，<Λ2,M—Y>表示求Λ2与(M-Y)的内积，β1、β2为正的惩罚参数；
S33、对增广拉格朗日函数进行求解，得到重构后的动态核磁图像
上述进一步方案的有益效果为：将模型对应的优化问题采用增广拉格朗日函数进行求解，主要是为了对变量进行分离，即将目标函数分别对每个优化变量进行求解。这样可使得问题简单化，且具有较高的求解效率。
进一步地：步骤S33包括以下分步骤：
S331、基于增广拉格朗日函数，固定辅助矩阵Y和M、矩阵形式的拉格朗日乘子Λ1和Λ2四个变量，对需要重构的动态核磁图像X进行求解，得到第k+1次迭代的动态核磁图像值：
A＝UF
其中，Xk+1为第k+1次迭代的动态核磁图像值，A为中间变量，Mk为第k次迭代的辅助矩阵M的值，
为第k次迭代的拉格朗日乘子Λ1的值，H为转置运算，I为单位矩阵，大小与矩阵AHA的大小相同，β1表示正的惩罚参数；
S332、基于增广拉格朗日函数，固定需要重构的动态核磁图像X、辅助矩阵M、矩阵形式的拉格朗日乘子Λ1和Λ2四个变量，对辅助矩阵Y进行求解，得到第k+1次迭代的辅助矩阵Y的值为：
Yt＝OY
其中，Yk+1为第k+1次迭代的辅助矩阵Y的值，Yt k为第t采样时刻下辅助矩阵Y对应的第k个核磁图像帧，T为总的采样时刻数量，argmin(·)为对等式右边的目标函数求最小值，将目标函数取得最小值时的极小点赋给等式左边的变量，Q为一个辅助的中间变量，OY为将辅助矩阵Y的第t个向量变换为相应的二维核磁图像帧Yt，Qt为第t采样时刻下的Q值，意为每个采样时刻t都有对应的Q值，
为第k次迭代的拉格朗日乘子Λ2的值，||Yt||TV为二维核磁图像帧Yt的总变分；
S333、基于增广拉格朗日函数，固定需要重构的动态核磁图像X、辅助矩阵Y、矩阵形式的拉格朗日乘子Λ1和Λ2四个变量，对辅助矩阵M进行求解，可进一步将增广拉格朗日函数表示为：
为计算方便，引入辅助变量S，并设
进一步转化为：
进而计算得到第k+1次迭代的辅助矩阵M的值为：
Mk+1＝U*diag(σ*)(V*)H
其中，Mk+1为第k+1次迭代的辅助矩阵M的值，U*、V*分别为矩阵S经奇异值分解后的两个正交矩阵，奇异值分解后得到S≈U*diag(σS)(V*)H，diag(σS)为矩阵S经奇异值分解后的对角矩阵，σS为矩阵S的奇异值向量，diag(σ*)为使函数
取得最小值时的极小点σ*组成的对角矩阵，即
σi(M)表示辅助矩阵M的第i个奇异值，μ＝β1+β2为惩罚因子，
为2范数的平方，σ为辅助矩阵M的奇异值向量，
为第t采样时刻辅助矩阵M的奇异值向量的极小点；
S334、采用第k+1次迭代的动态核磁图像值Xk+1、第k+1次迭代的辅助矩阵M的值Mk+1和第k+1次迭代的辅助矩阵Y的值Yk+1，对矩阵形式的拉格朗日乘子Λ1和Λ2进行更新，得到第k+1次迭代的矩阵形式的拉格朗日乘子Λ1的值和Λ2的值：
其中，
为第k+1次迭代的矩阵形式的拉格朗日乘子Λ1的值，
为第k+1次迭代的矩阵形式的拉格朗日乘子Λ2的值；
S335、基于第k+1次迭代的矩阵形式的拉格朗日乘子Λ1的值和Λ2的值，重新执行步骤S331～步骤S334，直到达到设定的迭代次数，将最终得到的动态核磁图像值作为重构后的动态核磁图像
上述进一步方案的有益效果为：在采用增广拉格朗日函数方法对重建模型进行迭代求解，引入相关辅助变量主要是为了简化优化问题。这样将提升重建模型的速度，另一方面进行多次迭代求解可进一步提高重建模型的精度。
综上，本发明的有益效果为：为了克服L1稀疏优化模型的不足，本发明针对基于压缩感知的核磁图像重建方法，设计了一种基于伽马核范数和总变分的核磁图像重建方法。本发明采用伽马核范数来刻画核磁图像的全局信息，同时利用总变分方法来描述每个核磁图像帧的局部一致性。通过结合伽马核范数和总变分方法，本发明充分挖掘了核磁图像重建重的空域和时域冗余性，可进一步增强核磁图像重建的鲁棒性能。
附图说明
图1为一种基于伽马核范数和总变分的核磁图像重建方法的流程图；
图2为下采样的核磁图像；
图3为重建后的核磁图像；
图4是重建后的核磁图像与原始核磁图像的差值图像；
图5为本发明方法与文献1、文献2方法的性能比较图。
具体实施方式
下面对本发明的具体实施方式进行描述，以便于本技术领域的技术人员理解本发明，但应该清楚，本发明不限于具体实施方式的范围，对本技术领域的普通技术人员来讲，只要各种变化在所附的权利要求限定和确定的本发明的精神和范围内，这些变化是显而易见的，一切利用本发明构思的发明创造均在保护之列。
如图1所示，一种基于伽马核范数和总变分的核磁图像重建方法，包括以下步骤：
S1、对原始的核磁图像帧进行下采样处理，得到下采样核磁图像帧；
步骤S1中进行下采样处理的公式为：
bt＝UFxt+nt
其中，bt为第t采样时刻的下采样核磁图像帧；U为下采样算子，F为2维离散傅里叶变换矩阵，此时，n为矩阵F的列，m为矩阵F的行，xt为第t采样时刻的原始的核磁图像帧，nt为第t采样时刻的噪声，除了矩阵F外，bt，U，xt，nt中的m均为矩阵行大小，n为矩阵列大小，
为复数域空间。
S2、基于下采样核磁图像帧，建立动态核磁图像重建模型；
步骤S2中动态核磁图像重建模型的公式为：
B＝[b1,b2,...,bT]H，t＝1,2,...,T，
其中，argmin(·)为对等式右边的目标函数求最小值，将目标函数取得最小值时的极小点赋给等式左边的变量，
为重构后的动态核磁图像，
为矩阵的F范数的平方，tr(·)为求一个矩阵的迹，U为下采样算子，F为2维离散傅里叶变换，X为需要重构的动态核磁图像，||X||γ为需要重构的动态核磁图像X的伽马核范数，Vec(·)为向量化算子，即Vec(xt),t＝1,2,...,T表示将
拉成一列向量，所以P＝mn，B为由T个下采样核磁图像帧构成的动态下采样核磁图像，T为总的采样时刻数量，bt为第t采样时刻的下采样核磁图像帧，H为转置运算，λ1和λ2为正则化参数，γ为伽马参数，min{P,T}为取P和T的最小值，P＝mn，P为需要重构的动态核磁图像X的行大小，σi(X)表示需要重构的动态核磁图像X的第i个奇异值，||X||TV为需要重构的动态核磁图像X的总变分，
为复数域空间。
需要重构的动态核磁图像X的总变分||X||TV的计算过程为：
A1、计算每个采样时刻的核磁图像帧的各向异性总变分：
其中，||xt||TV为第t采样时刻的核磁图像帧的总变分，xi,j,t表示第t采样时刻下原始核磁图像帧xt在(i,j)位置处的像素值，xi+1,j,t表示第t采样时刻下原始核磁图像帧xt在(i+1,j)位置处的像素值；i＝1,2,...,m—1，j＝1,2,...,n—1表示像素所在的位置下标；
A2、计算T个采样时刻的核磁图像帧的各向异性总变分的总和，得到需要重构的动态核磁图像X的总变分||X||TV：
其中，Xt为需要重构的动态核磁图像X的第t个向量，O为一个算子，该算子将需要重构的动态核磁图像X的第t个向量变换为相应的二维核磁图像帧。
在本实施例中：以两组心脏图像Cardiac cine和Cardiac perfusion作为测试数据，每个时刻下采样得到的核磁图像帧大小为256x256。下采样算子U设置为二进制采样矩阵，其中该二进制采样矩阵中非零项的数据决定了下采样率(元素1的个数除以整个矩阵的元素个数)，本实施例中选择采样率为0.20，T设置为30(即有30个核磁图像帧，每个核磁图像帧的大小均为256x256)。
S3、对动态核磁图像重建模型进行求解，得到重构后的动态核磁图像。
在对动态核磁图像重建模型求解时，需初始化相关参数；包括下采样率(本实施例设置为0.25)，核磁图像帧的个数(本实施例设置为30)。正则化参数λ1和λ2，其大小分别设置为0.9,0.01。伽马参数γ设置为0.01。
步骤S3包括以下分步骤：
S31、在进行模型求解时，引入了辅助变量M、Y，M、Y均与X等价，进一步将动态核磁图像重建模型表示为优化问题：
满足M＝X,Y＝M
其中，min(·)为最小化算子，即求括号里目标函数的最优值，U为下采样算子，F为2维离散傅里叶变换，X为需要重构的动态核磁图像，B为由T个下采样核磁图像帧构成的动态下采样核磁图像，H为转置运算，λ1和λ2为正则化参数，γ为伽马参数，σi(X)表示需要重构的动态核磁图像X的第i个奇异值，||X||TV为需要重构的动态核磁图像X的总变分，min{P,T}为取P和T的最小值，P＝mn，m为矩阵行大小，n为矩阵列大小，
为复数域空间，P为需要重构的动态核磁图像X的行大小，tr(·)为求一个矩阵的迹，M、Y均为与需要重构的动态核磁图像X等价的辅助矩阵，||X||γ为重构后的动态核磁图像X的伽马核范数，
为矩阵的F范数的平方；
S32、将优化问题构建为增广拉格朗日函数：
<Λ1,X-M>＝tr(Λ1 H(X-M))
<Λ2,M-Y>＝tr(Λ2 H(X-Y))
其中，
为最小化算子，具体为求关于优化变量X,Y,M,Λ1,Λ2下整个括号里的目标函数的最优值，
为矩阵的F范数的平方，||M||γ为M的伽马核范数，||Y||TV为Y的总变分，Λ1、Λ2为矩阵形式的拉格朗日乘子，<Λ1,X—M>表示求Λ1与(X—M)的内积，<Λ2,M—Y>表示求Λ2与(M—Y)的内积，β1、β2为正的惩罚参数；
S33、对增广拉格朗日函数进行求解，得到重构后的动态核磁图像
步骤S33包括以下分步骤：
在对增广拉格朗日函数求解前，先初始化拉格朗日乘子Λ1、Λ2和惩罚参数β1、β2，在本实施例中初始化Λ1＝Λ2＝0，β1＝0.1，β1＝0.05。k＝1表示第1次循环迭代，初始化总的循环迭代次数K为100。
S331、基于增广拉格朗日函数，固定辅助矩阵Y和M、矩阵形式的拉格朗日乘子Λ1和Λ2四个变量，对需要重构的动态核磁图像X进行求解，得到第k+1次迭代的动态核磁图像值：
A＝UF
其中，Xk+1为第k+1次迭代的动态核磁图像值，A为中间变量，Mk为第k次迭代的辅助矩阵M的值，
为第k次迭代的拉格朗日乘子Λ1的值，H为转置运算，I为单位矩阵，大小与矩阵AHA的大小相同，β1表示正的惩罚参数；
S332、基于增广拉格朗日函数，固定需要重构的动态核磁图像X、辅助矩阵M、矩阵形式的拉格朗日乘子Λ1和Λ2四个变量，对辅助矩阵Y进行求解，得到第k+1次迭代的辅助矩阵Y的值为：
Yt＝OY
其中，Yk+1为第k+1次迭代的辅助矩阵Y的值，Yt k为第t采样时刻下辅助矩阵Y对应的第k个核磁图像帧，T为总的采样时刻数量，argmin(·)为对等式右边的目标函数求最小值，将目标函数取得最小值时的极小点赋给等式左边的变量，Q为一个辅助的中间变量，OY为将辅助矩阵Y的第t个向量变换为相应的二维核磁图像帧Yt，Qt为第t采样时刻下的Q值，意为每个采样时刻t都有对应的Q值，
为第k次迭代的拉格朗日乘子Λ2的值，||Yt||TV为二维核磁图像帧Yt的总变分；
S333、基于增广拉格朗日函数，固定需要重构的动态核磁图像X、辅助矩阵Y、矩阵形式的拉格朗日乘子Λ1和Λ2四个变量，对辅助矩阵M进行求解，可进一步将增广拉格朗日函数表示为：
为计算方便，引入辅助变量S，并设
进一步转化为：
进而计算得到第k+1次迭代的辅助矩阵M的值为：
Mk+1＝U*diag(σ*)(V*)H
其中，Mk+1为第k+1次迭代的辅助矩阵M的值，U*、V*分别为矩阵S经奇异值分解后的两个正交矩阵，奇异值分解后得到S≈U*diag(σS)(V*)H，diag(σS)为矩阵S经奇异值分解后的对角矩阵，σS为矩阵S的奇异值向量，diag(σ*)为使函数
取得最小值时的极小点σ*组成的对角矩阵，即
σi(M)表示辅助矩阵M的第i个奇异值，μ＝β1+β2为惩罚因子，
为2范数的平方，σ为辅助矩阵M的奇异值向量，
为第t采样时刻辅助矩阵M的奇异值向量的极小点；
S334、采用第k+1次迭代的动态核磁图像值Xk+1、第k+1次迭代的辅助矩阵M的值Mk+1和第k+1次迭代的辅助矩阵Y的值Yk+1，对矩阵形式的拉格朗日乘子Λ1和Λ2进行更新，得到第k+1次迭代的矩阵形式的拉格朗日乘子Λ1的值和Λ2的值：
其中，
为第k+1次迭代的矩阵形式的拉格朗日乘子Λ1的值，
为第k+1次迭代的矩阵形式的拉格朗日乘子Λ2的值；
S335、基于第k+1次迭代的矩阵形式的拉格朗日乘子Λ1的值和Λ2的值，重新执行步骤S331～步骤S334，直到达到设定的迭代次数，将最终得到的动态核磁图像值作为重构后的动态核磁图像
为了验证本发明提出方法的核磁图像重建性能，本发明在实施例中选用了一组膝盖图像作为测试数据，每个时刻下采样得到的核磁图像帧大小均为256x256。下采样算子U设置为二进制采样矩阵。本发明中选择采样率为0.20，选择30个核磁图像帧。本发明方法仿真环境为64位Windows10操作系统、CPU2.6G，内存16G以及Matlab2018。
图2表示下采样的核磁图像，图3表示重建后的核磁图像，图4是重建后的核磁图像与原始核磁图像的差值图像。从图2～4的效果图可以看出本发明方法的有效性。
进一步，本发明与文献1：B.Tremoulheac,N.Dikaios,D.Atinson,et al.,DynamicMRI image reconstruction separation from undersampled(k,t)-space via low-rankplus sparse prior.IEEE Transactions on Medical Imaging,2014,33(8):1689-1701和文献2：S,G.Lingala,Y.Hu,E.DiBella,et al.,Accelerated dynamic MRIexpolitingsparsity and low-rank structure:k-t SLR.IEEE Transactions onMedical Imaging,2011,30(5):1042-1054中提出的核磁图像重建方法进行了性能比较。本发明与上述文献中的采样率均选择0.20，核磁图像帧个数均为30)，其他参数的配置根据本发明和文献1和文献2所提方法中相应的建议值进行设置，并将信号误差比作为性能比较的度量准则。其中信号误差比(Signal to Error Ratio,SER)的计算公式为：
其中Io表示无失真的原始核磁图像，
表示重建后得到的核磁图像。
表示矩阵的F范数的平方。如
表示Io的F范数的平方，如
图5给出了本发明方法与文献1、文献2方法的性能比较图，从图5中可以当采样率从0.1变化到0.45时，本发明方法的信号误差比值有明显的提高。进一步验证了本发明方法的信号误差比均比文献1、文献2方法的信号误差比性能要好。
综上所述，本发明的有益效果是通过结合伽马核范数和总变分方法，充分挖掘了核磁图像重建重的空域和时域冗余性，进一步提升了核磁图像重构的质量和鲁棒性能。
Claims (3)
1.一种基于伽马核范数和总变分的核磁图像重建方法，其特征在于，包括以下步骤：
S1、对原始的核磁图像帧进行下采样处理，得到下采样核磁图像帧；
S2、基于下采样核磁图像帧，建立动态核磁图像重建模型；
S3、对动态核磁图像重建模型进行求解，得到重构后的动态核磁图像；
所述步骤S1中进行下采样处理的公式为：
bt＝UFxt+nt
其中，bt为第t采样时刻的下采样核磁图像帧；U为下采样算子，F为2维离散傅里叶变换矩阵，此时，n为矩阵F的列，m为矩阵F的行，xt为第t采样时刻的原始的核磁图像帧，nt为第t采样时刻的噪声，除了矩阵F外，bt，U，xt，nt中的m均为矩阵行大小，n为矩阵列大小，
为复数域空间；
所述步骤S2中动态核磁图像重建模型的公式为：
B＝[b1,b2,...,bT]H，t＝1,2,...,T，
其中，argmin(·)为对等式右边的目标函数求最小值，将目标函数取得最小值时的极小点赋给等式左边的变量，
为重构后的动态核磁图像，
为矩阵的F范数的平方，tr(·)为求一个矩阵的迹，U为下采样算子，F为2维离散傅里叶变换，X为需要重构的动态核磁图像，||X||γ为需要重构的动态核磁图像X的伽马核范数，T为总的采样时刻数量，bt为第t采样时刻的下采样核磁图像帧，H为转置运算，λ1和λ2为正则化参数，γ为伽马参数，min{P,T}为取P和T的最小值，P＝mn，P为需要重构的动态核磁图像X的行大小，σi(X)表示需要重构的动态核磁图像X的第i个奇异值，||X||TV为需要重构的动态核磁图像X的总变分，
为复数域空间；
所述需要重构的动态核磁图像X的总变分||X||TV的计算过程为：
A1、计算每个采样时刻的核磁图像帧的总变分：
其中，||xt||TV为第t采样时刻的核磁图像帧的总变分，xi,j,t表示第t采样时刻下原始核磁图像帧xt在(i,j)位置处的像素值，xi+1,j,t表示第t采样时刻下原始核磁图像帧xt在(i+1,j)位置处的像素值；i＝1,2,...,m-1，j＝1,2,...,n-1表示像素所在的位置下标；
A2、计算T个采样时刻的核磁图像帧的总变分的总和，得到需要重构的动态核磁图像X的总变分||X||TV：
其中，Xt为需要重构的动态核磁图像X的第t个向量，O为一个算子，该算子将需要重构的动态核磁图像X的第t个向量变换为相应的二维核磁图像帧。
2.根据权利要求1所述的基于伽马核范数和总变分的核磁图像重建方法，其特征在于，所述步骤S3包括以下分步骤：
S31、将动态核磁图像重建模型表示为优化问题：
满足M＝X,Y＝M
其中，min(·)为最小化算子，即求括号里目标函数的最优值，U为下采样算子，F为2维离散傅里叶变换，X为需要重构的动态核磁图像，B为由T个下采样核磁图像帧构成的动态下采样核磁图像，H为转置运算，λ1和λ2为正则化参数，γ为伽马参数，σi(X)表示需要重构的动态核磁图像X的第i个奇异值，||X||TV为需要重构的动态核磁图像X的总变分，min{P,T}为取P和T的最小值，P＝mn，m为矩阵行大小，n为矩阵列大小，
为复数域空间，P为需要重构的动态核磁图像X的行大小，tr(·)为求一个矩阵的迹，M、Y均为与需要重构的动态核磁图像X等价的辅助矩阵，||X||γ为重构后的动态核磁图像X的伽马核范数，
为矩阵的F范数的平方；
S32、将优化问题构建为增广拉格朗日函数：
其中，
为最小化算子，具体为求关于优化变量X,Y,M,Λ1,Λ2下整个括号里的目标函数的最优值，
为矩阵的F范数的平方，||M||γ为M的伽马核范数，||Y||TV为Y的总变分，Λ1、Λ2为矩阵形式的拉格朗日乘子，<Λ1,X-M>表示求Λ1与(X-M)的内积，<Λ2,M-Y>表示求Λ2与(M-Y)的内积，β1、β2为正的惩罚参数；
S33、对增广拉格朗日函数进行求解，得到重构后的动态核磁图像
3.根据权利要求2所述的基于伽马核范数和总变分的核磁图像重建方法，其特征在于，所述步骤S33包括以下分步骤：
S331、基于增广拉格朗日函数，固定辅助矩阵Y和M、矩阵形式的拉格朗日乘子Λ1和Λ2四个变量，对需要重构的动态核磁图像X进行求解，得到第k+1次迭代的动态核磁图像值：
A＝UF
其中，Xk+1为第k+1次迭代的动态核磁图像值，A为中间变量，Mk为第k次迭代的辅助矩阵M的值，
为第k次迭代的拉格朗日乘子Λ1的值，H为转置运算，I为单位矩阵，大小与矩阵AHA的大小相同，β1表示正的惩罚参数；
S332、基于增广拉格朗日函数，固定需要重构的动态核磁图像X、辅助矩阵M、矩阵形式的拉格朗日乘子Λ1和Λ2四个变量，对辅助矩阵Y进行求解，得到第k+1次迭代的辅助矩阵Y的值为：
Yt＝OY
其中，Yk+1为第k+1次迭代的辅助矩阵Y的值，Yt k为第t采样时刻下辅助矩阵Y对应的第k个核磁图像帧，T为总的采样时刻数量，arg min(·)为对等式右边的目标函数求最小值，将目标函数取得最小值时的极小点赋给等式左边的变量，Q为一个辅助的中间变量，OY为将辅助矩阵Y的第t个向量变换为相应的二维核磁图像帧Yt，Qt为第t采样时刻下的Q值，意为每个采样时刻t都有对应的Q值，
为第k次迭代的拉格朗日乘子Λ2的值，||Yt||TV为二维核磁图像帧Yt的总变分；
S333、基于增广拉格朗日函数，固定需要重构的动态核磁图像X、辅助矩阵Y、矩阵形式的拉格朗日乘子Λ1和Λ2四个变量，对辅助矩阵M进行求解，可进一步将增广拉格朗日函数表示为：
为计算方便，引入辅助变量S，并设
进一步转化为：
进而计算得到第k+1次迭代的辅助矩阵M的值为：
Mk+1＝U*diag(σ*)(V*)H
其中，Mk+1为第k+1次迭代的辅助矩阵M的值，U*、V*分别为矩阵S经奇异值分解后的两个正交矩阵，奇异值分解后得到S≈U*diag(σS)(V*)H，diag(σS)为矩阵S经奇异值分解后的对角矩阵，σS为矩阵S的奇异值向量，diag(σ*)为使函数
取得最小值时的极小点σ*组成的对角矩阵，即
σi(M)表示辅助矩阵M的第i个奇异值，μ＝β1+β2为惩罚因子，
为2范数的平方，σ为辅助矩阵M的奇异值向量，
为第t采样时刻辅助矩阵M的奇异值向量的极小点；
S334、采用第k+1次迭代的动态核磁图像值Xk+1、第k+1次迭代的辅助矩阵M的值Mk+1和第k+1次迭代的辅助矩阵Y的值Yk+1，对矩阵形式的拉格朗日乘子Λ1和Λ2进行更新，得到第k+1次迭代的矩阵形式的拉格朗日乘子Λ1的值和Λ2的值：
其中，
为第k+1次迭代的矩阵形式的拉格朗日乘子Λ1的值，
为第k+1次迭代的矩阵形式的拉格朗日乘子Λ2的值；
S335、基于第k+1次迭代的矩阵形式的拉格朗日乘子Λ1的值和Λ2的值，重新执行步骤S331～步骤S334，直到达到设定的迭代次数，将最终得到的动态核磁图像值作为重构后的动态核磁图像
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110691450.XA CN113298907B (zh) | 2021-06-22 | 2021-06-22 | 一种基于伽马核范数和总变分的核磁图像重建方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110691450.XA CN113298907B (zh) | 2021-06-22 | 2021-06-22 | 一种基于伽马核范数和总变分的核磁图像重建方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN113298907A CN113298907A (zh) | 2021-08-24 |
| CN113298907B true CN113298907B (zh) | 2022-09-13 |
Family
ID=77329159
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202110691450.XA Active CN113298907B (zh) | 2021-06-22 | 2021-06-22 | 一种基于伽马核范数和总变分的核磁图像重建方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN113298907B (zh) |
Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN104077791A (zh) * | 2014-05-22 | 2014-10-01 | 南京信息工程大学 | 一种多幅动态对比度增强核磁共振图像联合重建方法 |
| CN106780645A (zh) * | 2016-11-28 | 2017-05-31 | 四川大学 | 动态mri图像重建方法及装置 |
| CN107330953A (zh) * | 2017-07-06 | 2017-11-07 | 桂林电子科技大学 | 一种基于非凸低秩的动态mri重建方法 |
| WO2020182033A1 (zh) * | 2019-03-08 | 2020-09-17 | 腾讯科技（深圳）有限公司 | 图像区域定位方法、装置和医学图像处理设备 |
Family Cites Families (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US9504851B2 (en) * | 2011-06-27 | 2016-11-29 | Koninklijke Philips N.V. | Magnetic resonance imaging of bone tissue |
-
2021
- 2021-06-22 CN CN202110691450.XA patent/CN113298907B/zh active Active
Patent Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN104077791A (zh) * | 2014-05-22 | 2014-10-01 | 南京信息工程大学 | 一种多幅动态对比度增强核磁共振图像联合重建方法 |
| CN106780645A (zh) * | 2016-11-28 | 2017-05-31 | 四川大学 | 动态mri图像重建方法及装置 |
| CN107330953A (zh) * | 2017-07-06 | 2017-11-07 | 桂林电子科技大学 | 一种基于非凸低秩的动态mri重建方法 |
| WO2020182033A1 (zh) * | 2019-03-08 | 2020-09-17 | 腾讯科技（深圳）有限公司 | 图像区域定位方法、装置和医学图像处理设备 |
Non-Patent Citations (5)
| Title |
|---|
| Compressed sensing reconstruction of undersampled 3D NOESY spectra: application to large membrane proteins;Bostock, MJ 等;《JOURNAL OF BIOMOLECULAR NMR》;20120930;第54卷(第1期);第15-32页 * |
| 一种优化的核磁共振图像重建方法;郭姿鹬;《电脑知识与技术》;20200205(第04期);第211-212+213页 * |
| 基于向量稀疏和矩阵低秩的压缩感知核磁共振图像重建算法;张红雨;《天津理工大学学报》;20170215;第33卷(第01期);第25-29页 * |
| 基于稀疏与低秩的核磁共振图像重构算法;敬朝阳等;《计算机应用研究》;20150331;第32卷(第03期);第942-945页 * |
| 基于非凸全变差正则的核磁共振图像重构算法;沈马锐等;《计算机应用》;20200810;第40卷(第08期);第2358-2364页 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN113298907A (zh) | 2021-08-24 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN104156994B (zh) | 一种压缩感知磁共振成像的重建方法 | |
| CN111028306A (zh) | 基于AR2 U-Net神经网络的快速磁共振成像方法 | |
| Zhao et al. | SwinGAN: A dual-domain Swin Transformer-based generative adversarial network for MRI reconstruction | |
| Bao et al. | Undersampled MR image reconstruction using an enhanced recursive residual network | |
| Wang et al. | Pyramid convolutional RNN for MRI reconstruction | |
| CN106618571A (zh) | 一种磁共振成像方法和系统 | |
| Shao et al. | SPECTnet: a deep learning neural network for SPECT image reconstruction | |
| CN112991483B (zh) | 一种非局部低秩约束的自校准并行磁共振成像重构方法 | |
| CN105654425A (zh) | 一种应用于医学x光图像的单幅图像超分辨率重建方法 | |
| CN111784792A (zh) | 基于双域卷积神经网络的快速磁共振重建系统及其训练方法与应用 | |
| Lingala et al. | Accelerated first pass cardiac perfusion MRI using improved k− t SLR | |
| Manimala et al. | Convolutional neural network for sparse reconstruction of MR images interposed with gaussian noise | |
| Xiao et al. | SR-Net: a sequence offset fusion net and refine net for undersampled multislice MR image reconstruction | |
| CN113476064B (zh) | 基于bcd-ed的单扫描双示踪剂pet信号分离方法 | |
| Barbano et al. | Steerable conditional diffusion for out-of-distribution adaptation in imaging inverse problems | |
| Fan et al. | An interpretable MRI reconstruction network with two-grid-cycle correction and geometric prior distillation | |
| CN110148193A (zh) | 基于自适应正交字典学习的动态磁共振并行重建方法 | |
| Usman et al. | Motion corrected multishot MRI reconstruction using generative networks with sensitivity encoding | |
| CN113298907B (zh) | 一种基于伽马核范数和总变分的核磁图像重建方法 | |
| Liu et al. | Highly undersampling dynamic cardiac MRI based on low-rank tensor coding | |
| Ahmed et al. | Dynamic MRI using deep manifold self-learning | |
| Yakkundi et al. | Convolutional LSTM: A deep learning approach for dynamic MRI reconstruction | |
| CN106228583A (zh) | 一种腹部图像重建的装置 | |
| Gao | Prior rank, intensity and sparsity model (PRISM): A divide-and-conquer matrix decomposition model with low-rank coherence and sparse variation | |
| Wang et al. | Multi-level temporal information sharing transformer-based feature reuse network for cardiac MRI 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 |