基于最优低秩约束的动态心脏磁共振图像重建算法
doi: 10.11918/202309027
段继忠 , 赵蕾 , 黄欢
昆明理工大学信息工程与自动化学院,昆明 650504
基金项目: 国家自然科学基金(61861023) ; 云南省基础研究计划项目(202301AT070452)
Dynamic cardiac magnetic resonance image reconstruction algorithm based on optimal low-rank constraints
DUAN Jizhong , ZHAO Lei , HUANG Huan
Faculty of Information Engineering and Automation, Kunming University of Science and Technology, Kunming 650504 , China
摘要
动态心脏磁共振成像(CMRI)是无创评估心血管疾病的重要工具。在动态 CMRI 中,常运用低秩张量恢复方法来探索动态磁共振图像的稀疏性,然而,张量的不同模态具有不同的低秩属性。研究发现,基于张量的非局部自相似性模态最能提升动态 CMRI 的重建质量。因此,在非局部低秩(NLR)方法的基础上,将高维图像中提取的每一组相似块视作一个矩阵,提出一种具有矩阵稀疏性的最优低秩矩阵恢复(OLRMR)模型。该模型使用加权 Schatten p-范数作为秩代理函数,采用交替方向乘子法(ADMM)和快速软阈值迭代算法求解。基于心脏数据集的实验结果表明,OLRMR 算法比 BCS、k-t SLR 和 k-t LRTC 算法更能有效提升重建图像的质量,更好保留图像细节与边缘轮廓信息。实验还表明,OLRMR 的重建速度比 k-t LRTC 提升了 2. 6 ~ 3 倍。
Abstract
Dynamic cardiac magnetic resonance imaging (CMRI) is an important tool for noninvasive assessment of cardiovascular disease. In dynamic CMRI, a low-rank tensor recovery method is usually employed to explore the sparsity of dynamic magnetic resonance images; however, different modes along the tensor have different low-rank properties. The studies have found that the nonlocal self-similarity mode along the tensor can best improve the reconstruction quality of dynamic CMRI. Therefore, this paper proposes an optimal low-rank matrix recovery (OLRMR) model with matrix sparsity based on the nonlocal low-rank (NLR) method by treating each set of similar blocks extracted from a high-dimensional image as a matrix. The model uses the weighted Schatten p-norm as the rank proxy function and was solved using the alternating direction multiplier method (ADMM) and a fast softthreshold iterative algorithm. Experimental results based on the cardiac dataset show that the OLRMR algorithm is more effective in improving the quality of the reconstructed image than the BCS, k-t SLR, and k-t LRTC algorithms and can better keep the detail of the image and the edge contour information intact. The experimental results also show that OLRMR improves the reconstruction speed by a factor of 2. 6-3 over k-t LRTC.
心血管疾病( cardiovascular disease,CVD)在 60 岁以上人群中较为常见,是全球主要死亡原因之一。心血管磁共振成像( cardiac magnetic resonance imaging,CMRI),是一种评估心血管系统的功能与结构的非侵入式医学成像技术。它运用磁共振成像(magnetic resonance imaging,MRI)基本原理,针对心血管系统成像的特殊困难,优化现有的 MRI 成像技术,使其具有临床价值。然而,要获得具有高时空分辨率、多对比度且全心脏覆盖的图像,动态 CMRI 需长时间采集大量数据。采集过程中还需要考虑呼吸引起的心脏运动,进一步延长了扫描时间[1]。此外,CMRI 扫描期间的呼吸运动可能导致重建图像出现伪影,传统 MRI 为避免呼吸运动产生伪影,通常要求检查者在整个成像期间屏住呼吸,这对于心脏病患者是比较困难的。采集时间过长也会影响图像质量,降低临床诊断的准确性。因此,减少动态磁共振的采样时间并开发能快速重建高质量图像的算法,成为动态 CMRI 重建领域的研究热点。脏病患者是比较困难的。采集时间过长也会影响图像质量,降低临床诊断的准确性。因此,减少动态磁共振的采样时间并开发能快速重建高质量图像的算法,成为动态 CMRI 重建领域的研究热点。
压缩感知( compressed sensing,CS)[2-3] 已广泛应用于磁共振成像,用于加速数据采集过程。特别是在动态 MRI 数据欠采样时,使用 CS-MRI 模型进行图像重建具有重要价值。 CS 还与灵敏度编码(sensitivity encoding,SENSE)[4]等并行 MRI 技术结合[5-6],通过多个接收线圈收集更多数据,从而提升重建图像的时空分辨率。另外,基于字典学习的 CS 技术也被用于 CMRI 重建。与传统的压缩感知技术不同,盲压缩感知(blind compressive sensing,BCS)[7] 从欠采样数据本身自适应地学习稀疏表示和字典,其优势是通过学习针对当前图像内容的专用字典[8],可以更好适应具体数据特征。
低秩性是稀疏性在矩阵上的拓展[9],矩阵秩最小化主要是指利用原始数据矩阵的低秩性进行矩阵重建,这涉及到最小化矩阵的秩函数。联合低秩与稀疏先验知识的图像压缩感知重建已被应用到高光谱图像重建[10]和磁共振成像研究。 Lingala 等[11] 提出 k-t SLR 方法将稀疏性和低秩约束相结合,以加速动态心脏磁共振成像。近年来,基于高阶张量的低秩近似模型研究日益增多,该模型可直接利用多维数据相关性,在动态 MRI 重建中表现出更优性能。然而,在这些基于高阶张量的模型中,动态 MRI 图像序列被简单地建模为具有低秩约束的全局张量,这易导致图像局部细节过度平滑。 Liu 等[12] 提出的 k-t LRTC 方法利用图像的非局部自相似性,引入基于非局部低秩( nonlocal low-rank,NLR)的张量编码算法来约束高维数据的稀疏性,利用高阶奇异值分解( higher-order singular value decomposition,HOSVD)[13-14]处理高维张量,挖掘数据的高维结构信息,提高了动态 CMRI 的重建质量。
在上述方法的基础上,本文在低秩张量建模中发现,张量的不同模态具有不同的低秩属性[15],即动态 CMRI 在时间帧(Nt)、相似块(Patch)和非局部自相似性(Nss)模态下的低秩属性存在差异,且各模态固有的低秩相关性对最终图像重建结果的贡献不同。本文首先通过高阶奇异值分解细致分析了构建的三阶张量沿各模态的秩性质,并通过实验发现沿张量的非局部自相似性模态能提升动态 CMRI 的重建质量; 随后,沿非局部模式进行重建,将高维图像中提取的每一组相似块视作一个矩阵,提出用于动态 CMRI 重建的最优低秩矩阵恢复( optimal lowrank matrix recovery,OLRMR)方法。该算法应用于动态磁共振成像,并基于心脏数据集进行实验,结果表明,最终重建效果优于 BCS、k-t SLR 和 k-t LRTC; 在重建时间上,本文算法较 k-t LRTC 更短。虽然本文重点针对动态心脏磁共振成像,但该算法同样适用于大多数动态 MRI 重建。
1 模型与方法
1.1 动态 MRI 问题
动态心脏磁共振图像(dynamic cardiac magnetic resonance imaging,dCMRI)重建模型可表示为
y=Ax+ε
(1)
式中:xCNx×Ny×Nt为待重建的动态图像序列; NxNy分别为一帧动态图像的水平方向和垂直方向的像素点数; Nt为动态图像序列的帧数; yCMxyNt.为 k-t 空间的测量数据; Mxy为一帧动态图像的欠采样数据点数; A=PFCMxyNt×NxNyNt.为傅里叶欠采样算子,其中P=INtP0RMxyNt×NxNyNt为动态图像序列的欠采样算子,P0RMxy×NxNy为单帧欠采样矩阵,INtNt×Nt的单位矩阵,为 Kronecker 积,F=INtUxUyCNxNyNt×NxNyNt为动态图像序列的傅里叶变换算子,UxCNx×NxUyCNy×Ny分别为 Nx 点和 Ny 点的离散傅里叶变换矩阵; εCMxyNt为测量噪声。研究目标是从测量的 k-t 空间样本 y 中恢复出信号 x。
1.2 最优低秩属性
传统的张量稀疏度度量常通过简单加和各模态的秩,将二阶稀疏度度量扩展到高阶情况( 如 k-t LRTC)。然而,此类方法忽略了一个事实,即高阶张量不同模态的秩性质并不一定相同,表明每个模态的秩与其固有的低秩子空间密切相关。为验证这一点,通过 HOSVD 探索三阶张量χC190×90×10经块匹配后各模态的低秩性质。如图1所示,从左到右依次是原始图像、k-t LRTC 模式以及沿非局部自相似性(Nss)、沿时间帧(Nt)和沿图像块(Patch)的低秩张量模式的图像重建结果。可以看出,沿非局部自相似性模式展开得到的结果比其他方法要好得多,重建结果真实地反映了各模态结构相关性的内在差异。 CMRI 在非局部自相似性模式上的低秩性明显强于时间帧和相似块模式,其次是 k-t LRTC 模式。因此,非局部自相似性是影响 CMRI 恢复性能的关键特性。
1各种低秩张量模式的重建图像
Fig.1Reconstructed images of various low-rank tensor patterns
图2是三阶张量χC190×90×10进行块匹配后的各种低秩模式在不同图像帧上的 PSNR 值。对于基于单模式和组合模式的低秩先验,非局部自相似性模式的性能最好,进一步验证了本文结论:沿非局部自相似性模式的结构相关性远强于时间帧或相似块模式。图2中,Patch、Nt、Nss 分别表示图像块、时间帧和非局部自相似性模式。
综上分析,非局部自相似性模式是促进动态 CMR 图像恢复的关键属性。由此,本文对动态 CMRI 进行重建时,不再通过 HOSVD 分解简单加和各模态的秩,而是直接沿非局部自相似性模式进行重建。
2各种低秩模式下不同图像帧的 PSNR 值对比
Fig.2Comparison of PSNR values of different image frames taken by various low-rank modes
1.3 基于 WSNM 的 OLRMR 模型
基于加权 Schatten p-范数(weighted schatten pnorm,WSNM),给定动态磁共振图像序列xCNx×Ny×Nt及其 k-t 空间的欠采样数据yCMxy×Nt则 CS-dCMRI 重建问题可以归纳为如下非局部低秩[16](NLR)约束最小化问题:
x=argminx 12Ax-y22+βi=1Np rankWi(x)
(2)
其中,第一项为数据保真项,第二项为包含重建图像先验信息的正则化函数。参数 β 是低秩编码正则项的权重; Wi:CNx×Ny×NtCmNt×c表示块匹配算子,即从动态图像中提取出m×m×Nt的图像块,在其局部邻域窗口中搜索出最相似(即欧几里得距离最小)的 c 个图像块,并重塑成mNt×c的矩阵; Np 是图像块的数量; i=1Np表示进行低秩矩阵近似时块索引。
利用近似点梯度法[17],可得公式在固定点 x j处的二次近似为
x=argminx 12Axj-y22+x-xj,AHAxj-y+L2x-xj22+βi=1Np rankWi(x)
(3)
其中L12Ax-y22的梯度的 Lipschitz 常数。忽略与固定点 x j相关的常数项并配方,则公式可转化为
x=argminx L2x-xj-1LAHAxj-y22+βi=1Np rankWi(x)
(4)
可以看出,近似点梯度法实质上就是将问题的光滑部分线性展开,再加上二次项并保留非光滑部分,然后求极小值作为每一步的估计。近似点梯度法的收敛速度为 O(1 / k)。
Beck 等提出的 FISTA [18-19]方法,对近似点梯度法进一步加速,收敛速度能达到 O(1 / k2 )。结合 FISTA 加速算法,公式可以转化为:
zj+1=x-j-1LAHAx-j-y
(5)
xj+1=argminx 12x-zj+122+μi=1Np rankWi(x)
(6)
tj+1=1+1+4tj22
(7)
x-j+1=xj+1+tj-1tj+1xj+1-xj
(8)
式中:zCNx×Ny×Nt; t为 FISTA 迭代的中间变量; x-为FISTA 加速中xjxj+1的线性组合矩阵; j 为迭代次数; μ = β / L。式(5)为梯度下降的迭代公式,式(6)为非局部低秩恢复子问题,式(7)和式(8)用来对xjxj+1进行线性组合得到新点x-j+1以便下一次迭代时在新点处做近似点梯度迭代。
对于非凸最优化问题,可以使用交替方向乘子法(alternating direction multiplier method,ADMM)[20-21]求解局部极值。引入辅助变量Bi=WixCmNi×e公式(6)可转化为以下优化问题:
minx,Bi 12x-zj+122+μi=1Np rankBi s. t. Bi=Wi(x)
(9)
引入 Bi 对应的拉格朗日乘子biCmNt×c则公式的增广拉格朗日函数为
Lx,Bi;bij=12x-zj+122+μi=1Np rankBi+α2i=1Np Bi-Wi(x)-bijF2-α2i=1Np bijF2
(10)
利用 ADMM 技术最小化与 Bix 相关的增广拉格朗日函数LxBi; bij从而得到 Bi x 的解,并更新拉格朗日乘子 bi。迭代子问题如下:
Bij+1=argminBi α2Bi-Wixj+bijF2+μrankBi
(11)
xj+1=argminx x-zj+122+αi=1Np Wi(x)-Bij+1-bijF2
(12)
bij+1=bij+ηWixj+1-Bij+1
(13)
式中:F为 Frobenius 范数; α >0 为惩罚项参数; η >0为乘子更新的参数。
1)首先求解公式。由于秩函数具有高度的非凸和非线性性质,是典型的 NP hard 问题。通常人们利用核范数逼近秩函数来优化求解,然而核范数往往过度收缩秩元素并且平等的对待每个秩元素,因此加权核范数[22] 被提出来赋予奇异值不同的权重,对其进行不同程度的收缩,使得更多有用的信息被保留。近几年,提出了一种加权 Schatten p-范数最小化[23-24],该方法不仅可以更好地逼近原始秩函数,而且还考虑了不同秩元素的重要性,所以加权 Schatten p-范数更适用于描述秩最小化问题。
为了使问题易于求解,采用约束于 Bi 上的非凸函数来增强奇异值的稀疏性,采用加权 Schatten p-范数代替秩函数。因此,问题可以重新表述为
Bij+1=argminBi 12Bi-Wixj+bijF2+τBiw,Spp
(14)
式中:τ=μ/α; wSpp为加权 Schatten p-范数,定义为Biwspp=k=1min{lm} wkσkpBi其中 0 <p <1,wk 为非负权重,表示为
wk=22cσk(1/p)+γ
(15)
式中:σk 为矩阵 Bi 的第 k 个奇异值; γ 是一个较小的值,用于防止分母为 0。公式是加权 Schatten p-范数最小化问题,可以根据加权奇异值阈值算法求解。设UΣV=Wixj+bij是矩阵Wixj+bij的奇异值分解,其中Σ=diagσ1σrr=min{lm}; 正交矩阵 UV 分别为Wixj+bij进行 SVD 分解的左、右奇异值矩阵,满足UUT=IVVT=I那么可以得到公式的解为
Bij+1=UsoftΣ,τwkVT
(16)
式中 soft(·)是软阈值算子,softΣτwk=maxΣ-diagτwk0wk=w1wminlm是权值。
2)求解公式,得到 x 的解为
xj+1=zj+1+αi=1Np Wi*Bij+1-bijI+αi=1Np Wi*Wi
(17)
式中:Wi*:CmNt×cCNx×Ny×NtWi的伴随算子,作用是将去噪后的矩阵块Bij+1放回原位置; i=1Np Wi*WCNxNyNt×NxNyNt是对角矩阵,其对角线元素等于对应像素点出现在所有块Bix中的总次数。
综上所述,通过对公式(5)~(8)的交替求解,最后得到重建图像,算法 1 为 OLRMR 算法的实现步骤。
2 实验结果与分析
2.1 实验设置
本文所有实验均在配置为 Intel Core i7-5500U@ 2.40GHz CPU、8GB 内存和 Windows 10 64 位操作系统的笔记本电脑上完成,运行环境为 MATLAB R2018b。
将所提 OLRMR 方法与经典的 BCS、k-t SLR 以及近年提出的 k-t LRTC 方法进行对比。在高度欠采样的动态 CMRI 数据上的实验表明,所提方法在重建质量上优于现有经典算法,在时间上优于 k-t LRTC 算法。参考 Liu 等的实验参数设置,k-t LRTC 执行 25 次内迭代和 10 次外迭代,张量块尺寸设定为 4 × 4 × Nt,正则化参数 μ 的初始值为 10-5; BCS 设置为 50 次内迭代和 25 次外迭代; k-t SLR 设置为 25 次内迭代和 25 次外迭代。对于本文提出的 OLRMR 算法,设定为 200 次迭代。正则化参数根据经验选择,并通过在各自数据集上进行手动调优以获得最佳结果。
本文方法是在不同的欠采样率( undersampling rate,UR)下重建动态 CMR 图像,并通过多种评价方法估计实验结果。本文通过视觉质量评估结合了图像常用的评价指标,包括峰值信噪比(peak signal-tonoise ratio,PSNR)、结构相似性(structural similarity,SSIM)[25] 和高频误差范数( high-frequency error norm,HFEN)[26],全面评估重建的动态 CMR 图像。其中 PSNR 和 SSIM 的数值越高,或者 HFEN 数值越低,则表示图像重建质量越好。
对于全采样的参考图像x^与重建图像xVPSNR定义为
VPSNR=10log10Vmax2VMSE
(18)
式中VMSEx^ x 的均方误差,Vmaxx^的最大像素值。
VSSIM定义为
VSSIM=2uxux^+c12σx^x^+c2ux2+ux2+c1σx2+σx^2+c2
(19)
式中:uxσx2分别为 x 的均值和方差; ux^σx^2分别为x^的均值和方差; σxx^xx^的协方差; c1c2 是常数。
VHFEN定义为
VHFEN = filter (x^)-filter(x)2filter(x^)2
(20)
式中 filter(·)为拉普拉斯-高斯滤波器,用于提取图像边缘信息。
2.2 实验数据
为了验证提出方法的性能,采用4 个动态 CMRI 数据集进行实验:
1)心肌灌注 MRI 数据集,该数据集来自犹他大学西门子 3T MRI 扫描的受试者,采用饱和恢复序列获取 70 帧 k 空间数据,单个切片的数据在 k 空间矩阵尺寸为 90 × 190(相位编码 × 频率编码)的笛卡尔网格上获得;
2)全采样的心脏电影 MRI 数据集,由延世大学医学中心的1.5T飞利浦扫描仪提供[27],该数据集由 25 个时间帧的 256 × 256 图像切片组成;
3)开放获取的心血管 MRI 数据集[28],该数据集原为 30 个线圈的多线圈数据,使用仿真单线圈(ESC)[29]方法得到单线圈数据,数据集由 18 帧尺寸为 112 × 288 的图像组成;
4)生理改进的非均匀心脏躯干(PINCAT)数值模型[30-31],数据集有 50 个时间帧,空间矩阵尺寸为 128 × 128。
所有数据集都使用伪径向欠采样,采用角间距均匀的径向轨迹掩模对每帧进行 k 空间欠采样,其中辐条之间的角间距根据指定的加速度因子设置。四个数据集的所有真实图像都使用笛卡尔全采样。
2.3 实验结果
首先,将提出的方法与经典的 BCS 和 k-t SLR 以及最近提出的 k-t LRTC 在不同伪径向欠采样率(UR = 4%~16%)的心肌灌注 MRI 重建中进行了比较。实验选取了部分心肌灌注图像序列进行重建比较,重建结果根据评价指标和视觉质量进行评估,如表1图3所示。表1给出了心肌灌注 MRI 数据集在不同欠采样率下不同算法重建的 PSNR、 SSIM 和 HFEN 的均值。从表1可以看出,与其他算法相比,OLRMR 算法的 PSNR 值和 SSIM 值都有了明显的提高,HFEN 值也减小了。 OLRMR 的 PSNR 值在 BCS、k-t SLR 和 k-t LRTC 的基础上平均提升了 6.43 dB、4.21 dB 和 2.44 dB。
图3中,对5 帧心肌灌注 MRI 进行了视觉比较。图中从上到下依次显示了全采样图像、BCS、k-t SLR、 k-t LRTC 和本文算法在 10% 欠采样率下的重建结果。其中,(a)是 5 帧全采样图像及误差图像,(b)~(e)分别是 BCS、 k-t SLR、 k-t LRTC 和 OLRMR 的 5 帧重建图像以及重建图像与参考图像之间的误差图像。从图3重建结果来看,在视觉效果上,本文方法比其他现有方法重建的图像更清晰,与原始图像最接近。从图3可以明显看到,BCS 和 k-t SLR 重建结果都有比较严重的模糊,尤其是在心肌区域。 k-t LRTC重建的图像存在细节和边缘模糊,虽然 OLRMR 重建的图像也存在一些细节模糊,但 k-t LRTC重建的图像有更多的重建伪影。从误差图像中也可以看出 OLRMR 算法的优越性,其他几种算法的误差图中存在大量较亮的点,说明误差较大。
3欠采样率为 10%的 5 帧心肌灌注 MRI 数据的重建图像与误差图像
Fig.3Reconstructed and error images of 5-frame myocardial perfusion MRI data with undersampling rate of 10%
为了客观评价各种重建算法的重建质量,图4中,(a)给出了心肌灌注 MRI 在欠采样率为 4%~24% 的采样模式下的平均 PSNR 值,(b)是心肌灌注 MRI 在 4% 伪径向欠采样率下 10 个时间帧的 PSNR 值。可以明显看出,在不同的欠采样率下,OLRMR 的 PSNR 都高于 BCS、k-t SLR 和 k-t LRTC 方法。虽然在高于24% 欠采样率情况下,OLRMR 的 PSNR 值与 k-t LRTC 算法很接近,但是随着欠采样率降低,OLRMR 算法表现出更大的优势,说明数据采集的越少,本文算法反而比其他算法重建效果更好。从图4(b)可以看出,对所有的图像帧,OLRMR 的重建效果都优于 BCS、k-t SLR 和 k-t LRTC 算法。
4心肌灌注 MRI 数据在不同欠采样率下和不同帧的平均 PSNR 值
Fig.4Mean PSNR of myocardial perfusion MRI data at different undersampling rates and at different frames
1不同欠采样率下心肌灌注 MRI 重建结果的平均 PSNR、SSIM 及 HFEN 对比
Tab.1Comparison of mean PSNR, SSIM, and HFEN for myocardial perfusion MRI reconstruction under different undersampling rates
图5分别显示了 4% 径向欠采样时心肌灌注 MRI 数据集与绿色虚线对应的时间剖面图和相应的误差图像。在图5中加入绿色虚线对应的时间剖面,分析心肌灌注图像随时间的运动情况,并评价重建图像的准确性。从图5可以看出,BCS 的时间轮廓出现严重模糊,并且时间轮廓存在较大的运动误差。 k-t SLR 和 k-t LRTC 的时间轮廓也存在模糊,且 k-t SLR 比 k-t LRTC 存在更多的运动细节误差。 OLRMR 能更好地保留运动图像的细微细节和边缘信息,且模糊程度较低,运动重建结果更准确。因此,OLRMR 方法在时间方向上的保真度优于其他 3 种方法。
为了充分验证所提出方法的鲁棒性和泛化性。采用心脏电影 MRI 数据集进行对比实验。表2给出了心脏电影 MRI 数据集在不同欠采样率下(UR = 4%~11%)各算法重建的平均 PSNR、SSIM 和 HFEN 值。从表2可见,OLRMR 的 PSNR 值比 BCS、 k-t SLR 和 k-t LRTC 平均提升了 5.11 dB、 1.73 dB和 1.62 dB。随着欠采样率减小,OLRMR 方法重建图像的质量提升更明显,表明只需要采集较少的数据即可重建出质量较好的图像。
5欠采样率为 4%的心肌灌注 MRI 数据集的时间剖面图与误差图
Fig.5Time profiles and error plots of myocardial perfusion MRI datasets with undersampling rate of 4%
2不同欠采样率重建心脏电影 MRI 时不同算法重建图像的平均 PSNR、SSIM 和 HFEN 对比
Tab.2Comparison of mean PSNR, SSIM, and HFEN of images reconstructed by different algorithms for cardiac cine MRI at different undersampling rates
图6选取了欠采样率为 4% 的重建结果进行对比。为了便于视觉比较,对全采样图像和重建图像的心脏区域(图6中的绿框)进行了放大,并绘制了重建图像与参考图像之间的误差图像。选取了 2 帧心脏电影 MRI 对应的重建图像和误差图像。通过比较发现,k-t SLR 和 k-t LRTC 的结果优于 BCS。从误差图中可以看出,BCS 重建后的图像中存在明显误差。与 k-t SLR 和 k-t LRTC 相比,OLRMR 重建的图像更平滑,边缘模糊程度较小,心脏部分保留了更多的图像细节信息。
6欠采样率为 4%的 2 帧心脏电影 MRI 数据的重建图像与误差图像
Fig.6Reconstructed images and error images of 2-frame cardiac cine MRI data with an undersampling rate of 4%
为了更加客观地评价所提算法的重建性能,采用欠采样率为 4%~16% 的心血管 MRI 数据集进行对比实验。表3给出了心血管 MRI 数据集在不同重建算法下的平均 PSNR、SSIM 和 HFEN 值。从表3可以得出,在不同欠采样率下,OLRMR 方法的 PSNR 和 SSIM 值均有一定提升,HFEN 值也有一定程度减小,总体上均优于其他重建算法,这与重建结果的视觉比较结果一致。 OLRMR 的 PSNR 值比 BCS、k-t SLR 和 k-t LRTC 平均提升了 6.7 dB、3.4 dB 和 4.2 dB。
3不同欠采样率重建心血管 MRI 时不同算法重建图像的平均 PSNR、SSIM 以及 HFEN 对比
Tab.3Comparison of mean PSNR, SSIM, and HFEN of images reconstructed by different algorithms for cardiovascular MRI at different undersampling rates
图7是心血管 MRI 数据集在欠采样率为 7.3% 下的重建结果,其中,(a)是 5 帧全采样图像及误差图像,(b)~( e)分别是 BCS、k-t SLR、k-t LRTC 和 OLRMR 的 5 帧重建图像以及重建图像与参考图像之间的误差图像。从图7重建结果来看,在视觉效果上,本文方法比其他现有方法重建的图像更清晰,与原始图像最接近。 BCS 重建结果有比较严重的模糊,k-t SLR 和 k-t LRTC 重建的图像虽比 BCS 清晰很多,但相比 OLRMR 重建的图像还是差很多。从误差图像也可以看出 OLRMR 算法的优越性,其他几种算法的误差图中存在大量较亮的点,说明误差较大。
7欠采样率为 7.3%的 5 帧心血管 MRI 数据的重建图像与误差图像
Fig.7Reconstructed and error images of 5-frame cardiovascular MRI data with undersampling rate of 7.3%
最后,用生理改进的非均匀心脏躯干(PINCAT)数据集再次验证所提算法性能,同样采用不同的欠采样率(6%~30%)对不同算法进行对比实验。如图8所示,选取具有代表性的 2 帧 PINCAT 数据集的重建结果,图( a)~( e)依次显示了全采样图像、 BCS、k-t SLR、k-t LRTC 和本文提出的 OLRMR 算法在 11% 欠采样率下的重建图像和误差图像。从图8可以看出,BCS 算法的重建图像存在大量的混叠伪影; k-t SLR 和 k-t LRTC 相较于 BCS 重建效果明显提升,但在一些空间细节和边缘部分仍存在模糊; 而 OLRMR 算法在视觉效果上与原始图像最为接近,从对应的误差图也可以看出,OLRMR 误差图中的亮点明显减少,说明误差最小,最大程度保留了图像细节,图像边缘部分也很清晰。
8欠采样率为 11%的 2 帧 PINCAT 数据的重建图像与误差图像
Fig.8Reconstructed image and error image of 2-frames of PINCAT data with undersampling rate of 11%
此外,实验还验证了本文算法在重建速度上的优势。由于本文算法是基于 k-t LRTC 算法改进而来,所以将本文算法和 k-t LRTC 方法的重建速度进行了比较,结果如图9所示。图9(a)为欠采样率为 11% 的心脏电影 MRI 数据集的重建结果对比,图9(b)是心肌灌注 MRI 数据集在 24% 欠采样率下的重建结果。可以看出,OLRMR 在提升重建质量的同时,其重建速度比 k-t LRTC 提升了 2.6 倍(心肌灌注 MRI 数据集)和 3 倍(心脏电影 MRI 数据集),表明本文方法可以实现更快速度的重建。
9不同算法不同数据集的重建速度比较
Fig.9Comparison of reconstruction speed of different algorithms for different datasets
2.4 实验参数分析
本文提出的 OLRMR 方法包含 μαp 3 个参数,通过遍历调整这 3 个参数可使欠采样数据集重建质量达到最佳。在对 4 个数据集进行图像重建时发现,对于不同的数据集,参数 α 取值较大时,实验结果几乎不再受其影响,可固定 α = 900。以 10% 欠采样率下的心肌灌注 MRI 数据集的重建为例,分析参数 μp 变化时 OLRMR 方法重建图像的 PSNR 变化情况,如图10所示。从图10可以看出,当 μ = 0. 007 和 p = 0.2 时,重建图像的 PSNR 达到最优。当 μ∈[0. 005,0. 01],p∈[0. 05,0.25]时,重建图像的 PSNR 值与最优 PSNR 值差值小于 0.1。
1010%欠采样率的心肌灌注 MRI 进行重建时参数 μp 的变化对 PSNR 值的影响
Fig.10PSNR value for different parameters μ and p when reconstruction of myocardial perfusion MRI data with undersampling rate of 10%
3 结语
低秩张量约束方法在动态磁共振图像重建中被广泛应用,但是哪种模式最有利于动态磁共振图像的恢复尚未有明确结论。本研究通过系统分析低秩张量的结构相关性,并结合实验验证,发现沿非局部自相似性模态施加低秩约束能有效提升动态心脏磁共振图像(CMRI)的重建质量。基于此,本文提出了一种用于动态 CMRI 重建的最优低秩矩阵恢复方法———OLRMR,并基于心脏数据集进行了实验。实验结果表明,与 BCS、k-t SLR 和 k-t LRTC 等方法相比,OLRMR 能重建出更清晰的图像,重建误差也最小,并且在 PSNR、SSIM 以及 HFEN 指标上均取得最佳的整体性能。在重建速度方面,本文算法相比最近提出的 k-t LRTC 也具有明显优势。实验进一步显示,欠采样率越低,OLRMR 对重建质量的提升效果越显著,说明该方法有助于加快动态磁共振成像数据采集的速度,同时也能有效保证重建图像的质量。
本文算法目前针对单线圈 CMRI 设计,未来可以扩展到多线圈采集模式,以进一步提升动态 CMRI 的重建效率。另外,OLRMR 在计算过程中需为每个参考块匹配与其相似的块,并对相似块矩阵进行 SVD 分解,计算开销较大; 由于块匹配和矩阵分解步骤相互独立,该算法非常适合使用 GPU 进行并行计算。后续研究将考虑使用 GPU 优化重建算法,以进一步提升动态 CMRI 的重建速度。
1各种低秩张量模式的重建图像
Fig.1Reconstructed images of various low-rank tensor patterns
2各种低秩模式下不同图像帧的 PSNR 值对比
Fig.2Comparison of PSNR values of different image frames taken by various low-rank modes
3欠采样率为 10%的 5 帧心肌灌注 MRI 数据的重建图像与误差图像
Fig.3Reconstructed and error images of 5-frame myocardial perfusion MRI data with undersampling rate of 10%
4心肌灌注 MRI 数据在不同欠采样率下和不同帧的平均 PSNR 值
Fig.4Mean PSNR of myocardial perfusion MRI data at different undersampling rates and at different frames
5欠采样率为 4%的心肌灌注 MRI 数据集的时间剖面图与误差图
Fig.5Time profiles and error plots of myocardial perfusion MRI datasets with undersampling rate of 4%
6欠采样率为 4%的 2 帧心脏电影 MRI 数据的重建图像与误差图像
Fig.6Reconstructed images and error images of 2-frame cardiac cine MRI data with an undersampling rate of 4%
7欠采样率为 7.3%的 5 帧心血管 MRI 数据的重建图像与误差图像
Fig.7Reconstructed and error images of 5-frame cardiovascular MRI data with undersampling rate of 7.3%
8欠采样率为 11%的 2 帧 PINCAT 数据的重建图像与误差图像
Fig.8Reconstructed image and error image of 2-frames of PINCAT data with undersampling rate of 11%
9不同算法不同数据集的重建速度比较
Fig.9Comparison of reconstruction speed of different algorithms for different datasets
1010%欠采样率的心肌灌注 MRI 进行重建时参数 μp 的变化对 PSNR 值的影响
Fig.10PSNR value for different parameters μ and p when reconstruction of myocardial perfusion MRI data with undersampling rate of 10%
1不同欠采样率下心肌灌注 MRI 重建结果的平均 PSNR、SSIM 及 HFEN 对比
Tab.1Comparison of mean PSNR, SSIM, and HFEN for myocardial perfusion MRI reconstruction under different undersampling rates
2不同欠采样率重建心脏电影 MRI 时不同算法重建图像的平均 PSNR、SSIM 和 HFEN 对比
Tab.2Comparison of mean PSNR, SSIM, and HFEN of images reconstructed by different algorithms for cardiac cine MRI at different undersampling rates
3不同欠采样率重建心血管 MRI 时不同算法重建图像的平均 PSNR、SSIM 以及 HFEN 对比
Tab.3Comparison of mean PSNR, SSIM, and HFEN of images reconstructed by different algorithms for cardiovascular MRI at different undersampling rates
BUSTIN A, FUIN N, BOTNAR R M,et al. From compressedsensing to artificial intelligence-based cardiac MRI reconstruction[J]. Frontiers in Cardiovascular Medicine,2020,7(2):1. DOI:10.3389/fcvm.2020.00017
LUSTIG M, DONOHO D, PAULY J M. Sparse MRI: The application of compressed sensing for rapid MR imaging[J]. Magnetic Resonance in Medicine,2007,58(6):1182. DOI:10.1002/mrm.21391
TOLOUEE A, ALIREZAIE J, BABYN P. Reconstruction of cardiac perfusion MRI with motion compensated compressed sensing[C]//2020 IEEE 5th Middle East and Africa Conference on Biomedical Engineering(MECBME). Marrakech, Morocco,2020:1. DOI:10.1109/MECBME47393.2020.9265153
PRUESSMANN K P, WEIGER M, SCHEIDEGGER M B,et al. SENSE: Sensitivity encoding for fast MRI[J]. Magnetic Resonance in Medicine,1999,42(5):952. DOI:10.1002/(SICI)1522-2594(199911)42:5 < 952:: AID-MRM16 > 3.0. CO;2-S
ARIANI F, BASARI B. Parameter analysis on sensitivity encoding(SENSE)algorithm for parallel imaging of magnetic resonance imaging[C]//2021 IEEE 5th International Conference on Information Technology, Information Systems and Electrical Engineering(ICITISEE). Yogyakarta, Indonesia,2021:35. DOI:10.1109/ICITISEE53823.2021.9655770
HU Yuhan, ZHANG Xinlin, CHEN Dicheng,et al. Spatiotemporal flexible sparse reconstruction for rapid dynamic contrast-enhanced MRI[J]. IEEE Transactions on Biomedical Engineering,2021,69(1):229. DOI:10.1109/TBME.2021.3091881
HE Ning, WANG Ruolin, WANG Yixue. Dynamic MRI reconstruction exploiting blind compressed sensing combined transform learning regularization[J]. Neurocomputing,2020,392(3):160. DOI:10.1016/j.neucom.2018.12.087
WANG Y H, YING L. Compressed sensing dynamic cardiac cine MRI using learned spatiotemporal dictionary[J]. IEEE Transactions on Biomedical Engineering,2014,61(4):1109. DOI:10.1109/TBME.2013.2294939
彭义刚, 索津莉, 戴琼海, 等. 从压缩传感到低秩矩阵恢复: 理论与应用[J]. 自动化学报,2013,39(7):981.PENG Yigang, SUO Jinli, DAI Qionghai,et al. Recovery from compressed sensing to low-rank matrices: Theory and applications[J]. Acta Automatica Sinica,2013,39(7):981. DOI:10.3724/SP. J.1004.2013.00981
孙玉宝, 吴泽彬, 吴敏, 等. 联合低秩与稀疏先验的高光谱图像压缩感知重建[J]. 电子学报,2014,42(11):2219.SUN Yubao, WU Zebin, WU Min,et al. Compressed sensing reconstruction of hyperspectral images with joint low-rank and sparse prior[J]. Acta Electronica Sinica,2014,42(11):2219. DOI:10.3969/j.issn.0372-2112.2014.11.014
LINGALA S G, HU Yue, DIBELLA E,et al. Accelerated dynamic MRI exploiting sparsity and low-rank structure:k-t SLR[J]. IEEE Transactions on Medical Imaging,2011,30(5):1042. DOI:10.1109/TMI.2010.2100850
LIU Die, ZHOU Jinjie, MENG Miaomiao,et al. Highly undersampling dynamic cardiac MRI based on low-rank tensor coding[J]. Magnetic Resonance Imaging,2022,89(4):12. DOI:10.1016/j.mri.2022.01.013
LATHAUWER L D, MOOR B, VANDEWALLE J. Multilinear singular value tensor decompositions[J]. SIAM Journal on Matrix Analysis and Applications,2000,24(4):1253. DOI:10.1137/S0895479896305696
BERGQVIST G, LARSSON E. The higher-order singular value decomposition: Theory and an application[lecture notes][J]. IEEE Signal Processing Magazine,2010,27(3):151. DOI:10.1109/MSP.2010.936030
CHANG Yi, YAN Luxin, CHEN Bingling,et al. Hyperspectral image restoration: Where does the low-rank property exist[J]. IEEE Transactions on Geoscience and Remote Sensing,2020,59(8):6869. DOI:10.1109/TGRS.2020.3024623
潘婷, 段继忠. 基于非局部低秩约束的改进灵敏度编码重建算法[J]. 数据采集与处理,2023,38(1):193.PAN Ting, DUAN Jizhong. An improved sensitivity encoding reconstruction algorithm based on nonlocal low-rank constraints[J]. Journal of Data Acquisition and Processing.2023,38(1):193. DOI:10.16337/j.1004-9037.2023.01.017
刘浩洋, 户将, 李勇锋, 等. 最优化: 建模、算法与理论[M]. 北京: 高等教育出版社,2020.LIU Haoyang, HU Jiang, LI Yongfeng,et al. Optimization: Modeling,algorithm and theory[M]. Beijing: Higher Education Press,2020
BECK A, TEBOULLE M. A fast iterative shrinkage-thresholding algorithm for linear inverse problems[J]. SIAM Journal on Imaging Sciences,2009,2(1):183. DOI:10.1137/080716542
BECK A, TEBOULLE M. Fast gradient-based algorithms for constrained total variation image denoising and deblurring problems[J]. IEEE Transactions on Image Processing,2009,18(11):2419. DOI:10.1109/TIP.2009.2028250
YI Xinlei, ZHANG Shengjun, YANG Tao,et al. Sublinear and linear convergence of modified ADMM for distributed nonconvex optimization[J]. IEEE Transactions on Control of Network Systems,2022,10(1):75. DOI:10.1109/TCNS.2022.3186653
侯榆青, 金明阳, 贺小伟, 等. 基于随机变量交替方向乘子法的荧光分子断层成像[J]. 光学学报,2017,37(7):195.HOU Yuqing, JIN Mingyang, HE Xiaowei. Fluorescence molecular tomography using a stochastic variant of alternating direction method of multipliers[J]. Acta Optica Sinica.2017,37(7):195. DOI:10.3788/AOS201737.0717001
LIU Xuya, HAO Caiyan, SU Zezhao,et al. Image inpainting algorithm based on tensor decomposition and weighted nuclear norm[J]. Multimedia Tools and Applications,2022,82(3):1. DOI:10.1007/s11042-022-12635-3
WANG Li, HOU Wensheng, WU Xiaoying,et al.3D MR image denoising using a modified adaptive high order singular value decomposition method[C]//2020 42nd Annual International Conference of the IEEE Engineering in Medicine and Biology Society(EMBC). Montreal, QC, Canada,2020:1580. DOI:10.1109/EMBC44109.2020.9175418
蒋明峰, 陆亮, 吴龙, 等. 基于加权 Schatten p-范数最小化的磁共振图像重构方法研究[J]. 电子学报,2019,47(4):784.JIANG Mingfeng, LU Liang, WU Long,et al. The Research of MRI reconstruction method by using weighted Schatten P-norm minimization[J]. Acta Electronica Sinica,2019,47(4):784. DOI:10.3969/j.issn.0372-2112.2019.04.003
WANG Z, BOVIK A, SHEIKH H,et al. Image quality assessment: From error visibility to structural similarity[J]. IEEE Transactions on Image Processing,2004,13(4):600. DOI:10.1109/TIP.2003.819861
RAVISHANKAR S, BRESLER Y. MR image reconstruction from highly undersampled k-space data by dictionary learning[J]. IEEE Transactions on Medical Imaging,2011,30(5):1028. DOI:10.1109/TMI.2010.2090538
JUNG H, SUNG K, NAYAK K S,et al.k-t FOCUSS: A general compressed sensing framework for high resolution dynamic MRI[J]. Magnetic Resonance in Medicine,2009,61(1):103. DOI:10.1002/mrm.21757
CHEN Chong, LIU Yingmin, SCHNITER P. OCMR(v1.0)-openaccess multi-coil k-space dataset for cardiovascular magnetic resonance imaging[J]. Image and Video Processing,2020:1. DOI:10.48550/arXiv.2008.03410
TYGERT M, ZBONTAR J. Simulating single-coil MRI from the responses of multiple coils[J]. Image and Video Processing,2018,15(2):115. DOI:10.2140/camcos.2020.15.1
SEGARS W P, TSUI B M. Study of the efficacy of respiratory gating in myocardial SPECT using the new 4-D NCAT phantom[J]. IEEE Transactions on Nuclear Science,2002,49(3):675. DOI:10.1109/TNS.2002.1039548
SHARIF B, BRESLER Y. Adaptive real-time cardiac MRI using PARADISE: Validation by the physiologically improved NCAT phantom[C]//2007 4th IEEE International Symposium on Biomedical Imaging: From Nano to Macro. Arlington, VA, USA,2007:1020. DOI:10.1109/ISBI.2007.357028