基于小波与剪切波的压缩感知纵向MR成像
康瑞瑞1, 曹斌2, 王斌3, 闫龙4, 渠刚荣1,5
1.北京交通大学 理学院,北京100044
2.湘潭大学 信息工程学院,湖南 湘潭 411105
3.北京航空航天大学 精密光机电一体化技术教育部重点实验室,北京 100191
4. 延安市人民医院MRI、CT诊断科,陕西 延安 716000
5. 北京数学与信息交叉学科中心,北京 100048
通讯作者:渠刚荣(1961—),男,内蒙古呼和浩特人,教授.email: grqu@bjtu.edu.cn.

第一作者:康瑞瑞(1984—),女,陕西清涧人,博士生.研究方向为图像重建.email:11118381@bjtu.edu.cn.

摘要

针对传统二维小波变换不能最优表示MR图像和MR图像数据采集缓慢的问题,基于相似性与联合稀疏变换对纵向压缩感知磁共振成像(LCS-MRI)提出新的正则化模型与相应的重建算法(frsLCS-MRI).为了验证联合稀疏变换与算法frsLCS-MRI的有效性,利用两个医学数据集与基于小波变换和相似性的算法LACS-MRI作了比较.数值实验表明:frsLCS-MRI算法可重建更高精度和更高信噪比的MR图像.使用联合稀疏变换比使用小波变换或剪切波变换可更好地重建MR图像的细节信息,如边界、拐角、轮廓、二维奇异曲线等.此外,利用相似性先验信息可节省数据采集时间.

关键词: 图像处理; 磁共振成像; 压缩感知; 相似性; 小波; 剪切波
中图分类号:TN957.52;TP391.41;R445.2 文献标志码:A 文章编号:1673-0291(2017)03-0127-06
Compressed sensing longitudinal MRI based on wavelet and shearlet
KANG Ruirui1, CAO Bin2, WANG Bin3, YAN Long4, QU Gangrong1,5
1. School of Science, Beijing Jiaotong University, Beijing 100044,China
2. College of Information Engineering, Xiangtan University, Xiangtan Hunan 411105, China
3. Key Laboratory of Precision Opto-mechatronics Technology of Education Ministry,Beihang University, Beijing 100191, China
4. MRI, CT Diagnosis Branch, Yan'an People's Hospital, Yan'an Shaanxi 716000, China
5. Beijing Center for Mathematics and Information Interdisciplinary Sciences, Beijing 100048, China
Abstract

In order to overcome the deficiencies that the traditional two-dimensional wavelet transform can not provide the optimal representation for MR images and to reduce the scanning time of MR images, a new regularization model and the corresponding reconstruction algorithm (frsLCS-MRI) are proposed based on the joint sparse transform and similarity prior. In order to verify the effectiveness of the joint sparse transform and the proposed algorithm frsLCS-MRI, we compare it with the LACS-MRI algorithm based on single sparse transform and similarity by using two medical datasets. Numerical results indicate that the proposed method frsLCS-MRI can reconstruct MR images with higher accuracy and higher signal to noise ratio. Compared with the method LACS-MRI using only wavelet transform, or shearlet transform, more detail information such as boundaries, corners, contours, can be reconstructed by the proposed method frsLCS-MRI using the joint sparse transform. Moreover, exploiting similarity based on reference image saves acquisition time.

Keyword: image processing; magnetic resonance imaging (MRI); compressed sensing; similarity; wavelet; shearlet

磁共振成像(Magnetic Resonance Imaging, MRI)是一种非侵入性和非电离的成像技术, 它提供了多种对比机制, 很好地实现了解剖结构和生理功能的可视化.但是, 由于MRI中的数据是在时间线上连续地采集 K⁃空间(图像的Fourier变换空间)获得的, 因此MRI速度缓慢, 影响了临床吞吐量与成像质量.目前, 如何精确重建MR图像中的细节信息, 包括边缘、纹理、拐角、曲线等成为MRI的一大挑战.

最近发展起来的压缩感知(Compressed Sensing, CS)理论表明, 从高度欠采样的K-空间中可以以压倒性的概率精确重建MR图像[1], 这可以极大地减少MR图像的扫描时间和减轻病人的不适感.此外, 在许多临床MRI中, 利用已经存在的图像信息可以进一步缩短数据采集时间并提高重建图像的信噪比.这里把已经得到的采集图像称为参考图像, 参考图像与被重建的目标图像之间具有相似性[2, 3].Lior Weizman等[2]利用这种相似性与MR图像在小波变换下的稀疏性, 对纵向压缩感知磁共振成像(Longitudinal Compressed Sensing Magnetic Resonance Imaging, LCS-MRI), 提出一种自适应采样和加权重建算法LACS-MRI, 进一步地减少了MR图像扫描时间, 并提高了重建MR图像的信噪比.

虽然文献[2]利用时间域的相似性减少了采样, 但由于其利用了小波变换作为空间域的稀疏变换, 使得重建的图像依然存在类似于曲线状的伪影.这主要是因为二维小波变换是由一维小波变换张量积构成的, 小波变换可以有效地逼近具有点奇异性的一维图像特征但不能高效逼近具有细节信息等高维奇异性的图像特征.与小波变换相比, 剪切波变换[4]能够更加有效地逼近含有如边界、拐角、曲线、尖峰等几何信息的逐段光滑图像特征, 具有多尺度方法提取图像几何信息的强大能力.但是剪切波变换不能有效表示类似于点状的图像特征.很自然的想法是将两者相结合.小波变换与剪切波变换的线性组合构成新的稀疏变换(称为联合稀疏变换), 可以避免重建图像的点状伪影与曲线状的伪影, 表示更丰富的几何图像信息, 既弥补了单一方法的不足又保留了各自的优点.

本文作者利用相似性先验信息, 结合小波变换和剪切波变换, 提出一种新的正则化模型.并为了处理模型中不光滑的正则项和实现快速收敛的目的, 给出一种新的求解算法frsLCS-MRI.数值实验表明, 使用联合变换比只使用小波变换或剪切波变换可以重建更多几何信息和更高信噪比的MR图像, 而且使用相似性减少了MR图像的扫描时间.

1 提出的模型与算法

CS是一种有效的数学框架, 只要信号是稀疏的或者测量是不相干的, 就可以使用非常少的, 远低于奈奎斯特采样定率所要求的采样重建高维信号.压缩感知磁共振成像(Compressed Sensing Magnetic Resonance Imaging, CS-MRI)重建模型为

MinxBx-y22+αΨx1(1)

式中: xCN是所要重建图像的向量形式, 其中 CNN维复数域; yCM是观测向量; ΨCM×N是一个稀疏变换矩阵, BCM×N( MN)为部分Fourier变换矩阵 FuCM×N, 即 B=Fu; 正则化参数 α是用来调节保真项与稀疏性之间的平衡.

在许多MRI的应用中, 参考图像与被重建的目标图像之间具有相似性.记参考图像为 x0CN.在MRI中参考图像一般是同一个扫描切片在某个对比度下的图像或者是同一被诊断者不同时期的扫描图像.

由于 x0x在大部分图像域是相似的[2, 3, 5], 因此差分 x0-x可以模式化为稀疏性, 把时间相似性 x0-x嵌入到CS-MRI重建模型(1)中就得到

MinxFux-y22+αΨx1+βx-x01(2)

其中 αβ是正则化参数, 用来平衡变换域与时间域之间的稀疏性.

尽管利用相似性先验信息可以加速采样、减少扫描时间、提高信噪比.但是, 由于有些扫描切片之间的时间间隔较长, 给LCS-MRI带来很大的挑战.这一方面是因为时间间隔较长的扫描切片之间的相似性可能是不存在的, 如有的扫描切片之间发生了大的病理改变或外科技术的改变.另一方面是因为在每个时间点都要求重建高质量的MR图像, 时间域欠采样会导致重建的MR图像不理想.Lior Weizman等[2]根据相似性的相似程度提出一种自适应采样和加权重建算法LACS-MRI.用公式表示如下:

MinxFux-y22+αW1Ψx1+βW2(x-x0)1(3)

其中 W1RN×NW2RN×N是加权矩阵, RN×NN×N维实数域, 并有

Wk=diag([wk1, wk2, , wkN]), k=1, 2,

且对 i=1, 2, , N,

w1i=11+[|Ψ1x˙|]iw2i=11+[|x˙-x0|]i(4)

式中: x˙是被重建图像x的一个估计; ·i表示方括号中向量的第 i个元素, 且有 0wki1.

1.1 提出的模型

为了重建MR图像的更多方向特征与高阶奇异性, 本文将小波与剪切波的线性组合作为新的稀疏变换, 并结合相似性先验提出的新模型如下:

MinxFux-y22+αaW1aΨax1+αbW1bΨbx1+βW2(x-x0)1(5)

式中: Ψa, ΨbCM×N分别是小波变换与剪切波变换矩阵; αa, αbβ是正则化参数; W1a, W1bRN×NW2RN×N是加权矩阵, 并有

Wk=diag([wk1, wk2, , wkN]), k=1a, 1b, 2,

且对 i=1, 2, , N,

w1ai=1/(1+[|Ψax˙|]i)w1bi=1/(1+[|Ψbx˙|]i)w2i=1/(1+[|x˙-x0|]i)(6)

1.2 frsLCS-MRI算法

将基于相似性与联合稀疏变换提出快速LCS-MRI成像算法记为frsLCS-MRI, 简称算法1, 详细过程见表1.

表1 frsLCS-MRI算法 Tab.1 frsLCS-MRI algorithm

注意到算法1中并不假设相似性存在, 因此在第2步加权矩阵初始化为 W1a=W1b=IW2=0.

1.3 rsSFISTA算法

快速阈值迭代算法(FISTA)[6]原本是求解下面的问题:

MinxBx-y22+αx1(7)

最近 Tan等[7]基于FISTA提出一个光滑快速阈值迭代算法(SFISTA), 求解(7)的加权形式问题

MinxBx-y22+αDx1(8)

其中 D是分析字典, Dx是稀疏的.

我们对算法SFISTA进行扩展, 提出算法rsSFISTA, 简称算法2, 用来求解最小化问题(5).为方便记述, 我们给出下面形式记法

f(x)=Fux-y(9)

g1(W1aΨax)=minu{αau1+12μu-W1aΨax22}(10)

g1(W1bΨbx)=minu{αbu1+12μu-W1bΨbx22}(11)

g2μ(W2(x-x0))=minu{βu1+12μu-W2(x-x0)22}(12)

以及

Hμ(x)=Fu-y22+λag1(W1aΨax)+λbg1(W1bΨbx)+g2μ(W2(x-x0))(13)

式中: g1, g1, g2μ分别是 αau1, αbu1, βu1的莫罗包络(Moreau envelope); 算子 f, g1, g1g2μ的梯度分别记为 f, g1, g1g2μ; ·2表示矩阵的最大奇异值; 算子 Γαμ(z)是作用在向量 z元素上的软阈值算子, 定义为

Γαμ(z)=|zi|-αμzizi,   |zi|> αμ0,         否则(14)

其中 αμ是阈值且 α> 0, μ> 0.矩阵 H的共轭转置记为 H* , 若 H是实值矩阵, 则 H* 就为 H的转置.

rsSFISTA算法的具体过程见表2.

表2 求解最小化问题(5)的rsSFISTA算法 Tab.2 rsSFISTA algorithm for solving minimization problem (5)

算法2中两个稀疏性的平衡性由 λaαa+λbαbβ之间的比率决定, 可通过函数 Γ(·)实现, 算法2的收敛性是由 μ来控制的.

2 数值实验

本节从不完全测量数据中重建高质量的MR图像来验证本文所提出的模型与算法.我们在两个临床纵向MRI真实数据集上进行了测验, 并与基于相似性与小波变换的算法LACS-MRI作了比较.其中真实数据与部分源码来自于文献[8].两个数据集的部分K-空间数据都是通过对一个全采样的K-空间进行降采样获得的.为了简单起见, 本文使用变密度随机采样机制和快速扫描机制.变密度随机采样机制是根据数据的概率密度函数进行采样, 概率密度函数是根据距离K-空间中心的远近来定义的, 也就是距离K-空间中心越近采样越多, 距离中心越远采样越少[2, 3].

本文使用Daubechies 4小波变换和ShearLab 3D[9, 10]的离散不可分离剪切波变换作为稀疏变换.其中离散不可分离剪切波变换使用默认设置, 也就是对0, 1, 2, 3尺度分别设有4, 4, 8, 8方向滤波的剪切滤波.总共有25个子带.在下面两个实例中, 参数 λa, λb, μ的值都分别设为 λa=0.31, λb=0.69, μ=0.0364.

2.1 同一切片在T2-加权与FLAIR下的相似性

同一被诊断者的同一个切片在多个对比下都具有相似性, 特别是T2-加权对比与流体衰减反相恢复(Fluid-Attenuated Inversion Recovery, FLAIR)对比的图像在非流体区域具有高度的相似性, 利用这种相似性可以减少采样时间并提高重建MR图像的信噪比.本实例(例1)是基于一个扫描切片在不同对比度下的相似性的MRI, 具体就是基于非流体区域在T2-加权与FLAIR对比下具有高度相似性, 使用FLAIR的15%采样重建FLAIR图像.其中全采样的T2-加权对比扫面切片作为参考图像, 加权矩阵 W1a, W1bW2按照算法1来迭代计算.对只采用小波变换作为稀疏变换的正则化参数设为 αa=0.025, αb=0, β=0.03; 对只采用剪切波变换作为稀疏变换的参数设为 αa=0, αb=0.025, β=0.03.对采用联合稀疏变换作为稀疏变换的正则化参数设为 αa=0.0077, αb=0.0172, β=0.03.

重建的FLAIR图像如图1所示.从图1的(e)和(f)看出小波能够有效地表示点奇异性图像特征而不能够有效表示曲线奇异性几何特征.从图1的(g)和(h)看出剪切波能够更好地表示边界与曲线, 而对点的奇异性表示则差一些.从图1的(i)和(j)看出, 使用小波与剪切波的线性组合作为稀疏变换比只使用单个稀疏变换(小波变换或剪切波变换)作为稀疏变换, 可以重建更多FLAIR图像细节特征, 包括点奇异性、拐角、边界轮廓曲线等, 并且还抑制了只用小波变换引起的方块状伪影和只用剪切波变换引起的曲线状伪影.

图1 采样方式和使用不同稀疏变换从K-空间的15%数据中重建的图像Fig.1 Sampling pattern and the images reconstructed from 15% of K⁃space data by using different sparse transforms

2.2 同一被诊断者不同时间扫描切片间的相似性

患有脑部肿瘤的病人经常每隔几周或几个月就会进行一次脑部磁共振扫描.同一个病人在时间线上的不同扫描切片之间具有相似性, 利用这种相似性能够加速后续扫描切片的重建.本实例(例2)基于同一个被诊断者在时间(如数周或数月)线上不同扫描切片之间的相似性, 重建后续MR图像.这里把先前扫描的切片作为参考图像.对使用小波变换或剪切波变换作为稀疏变换时, 正则化参数的取值范围设为 αa, αb, β[0.001, 0.9]; 对使用小波变换与剪切波变换的线性组合作为稀疏变换时, 正则化参数取值范围设为 αa, αb[3.1×10-4, 0.621], β[0.001, 0.9].

重建的后续MR图像如图2所示.从图2看出, 使用小波变换与剪切波变换的线性组合作为稀疏变换, 比使用小波或剪切波单个变换作为稀疏变换可重建更高质量和更高精度的MR图像.只用单个稀疏变换不能重建使用联合变换重建的某些几何图像特征, 这是因为小波与剪切波都有一定局限性, 在抑制某些伪影的同时会引起另一些伪影.两者相结合可以弥补各自的不足, 使得重建的脑部图像有效地抑制了正方形状伪影和曲线状伪影, 并提高了重建图像的峰值信噪比(Peak Signal-to-Noise Ratio, PSNR)与质量.

图2 采样方式与使用不同稀疏变换从K-空间的6%数据中重建的图像Fig.2 Sampling pattern and the images reconstructed from 6% of K⁃space data by using different sparse transforms

为了定量说明本文所提出方法的性能, 我们考察了两个数据集的PSNR, 定义如下:

P=20lg(M/V)(15)

式中:P为峰值信噪比; M是图像中最大像素值; V是原始图像 x与所要重建图像 x˙的均方误差.

表3 基于相似性使用不同稀疏变换重建图像的PSNR对比 Tab.3 PSNR contrast of images reconstructed using different sparse transforms based on similarity

表3的结果与上述两个例子的结论一致, 即使用联合稀疏变换比使用单个稀疏变换重建的MR图像具有更丰富的图像特征和更高的PSNR.此外, 与第二个例子相比, 第一个例子的PSNR值的变化范围相对比较小, 这解释了T2-加权对比与FLAIR对比之间的相似性只是存在于部分区域而不是整个图像.

3 结论

本文利用MRI已有的图像与所要重建图像之间的相似性, 结合小波与剪切波提出了新的正则化模型.为了处理模型中的不可微项和达到快速收敛的目的, 对算法SFISTA进行了扩展, 提出frsLCS-MRI算法.在两个临床纵向MRI真实数据集上进行了测验, 与算法LACS-MRI作了比较, 结果表明, 使用联合变换比使用单个变换作为稀疏变换可以更好地重建图像的细节、纹理与边界等特征, 避免了由于只使用单个变换引起的伪影.另外利用相似性减少了MR图像的扫描时间.

The authors have declared that no competing interests exist.

参考文献
[1] 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-1195. [本文引用:1]
[2] WEIZMAN L, EIADR Y C, EILAM A, et al. Fast reference based MRI[C]//37th Annual International Conference of the IEEE, 2015: 7486-7489. [本文引用:6]
[3] WEIZMAN L, ELDAR Y C, BASHAT D B. Compressed sensing for longitudinal MRI: An adaptive-weighted approach[J]. Medical Physics, 2014, 42(9): 5195-5208. [本文引用:3]
[4] EASLEY G, LABATE D, WANG Q L. Sparse directional image representations using the discrete shearlet transform[J]. Applied and Computational Harmonic Analysis, 2008, 25(1): 25-46. [本文引用:1]
[5] GAMPER U, BOESIGER P, KOZERKE S. Compressed sensing in dynamic MRI[J]. Magnetic Resonance in Medicine, 2008, 59(2): 365-373. [本文引用:1]
[6] BECK A, TEBOULLE M. A fast iterative shrinkage-thresholding algorithm for linear inverse problems[J]. SIAM Journal on Imaging Sciences, 2009, 2(1): 183-202. [本文引用:1]
[7] ZHAO T, ELDAR Y C, BECK A, et al. Smoothing and decomposition for analysis sparse recovery[J]. IEEE Transactions on Signal Processing, 2014, 62(7): 1762-1774. [本文引用:1]
[8] WEIZMAN L, ELDAR Y C. Reference based MRI[EB/OL]. [2016-03-22]. http://webee.technion.ac.il/people/YoninaEldar/software_det13.php. [本文引用:1]
[9] KUTYNIOK G, WANG Q L, REISENHOFER R. ShearLab 3D: Faithful digital shearlet transforms based on compactly supported shearlets[J]. ACM Transactions on Mathematical Software, 2016, 42(1): 5. [本文引用:1]
[10] KUTYNIOK G, SHAHRAM M, ZHUANG X S. ShearLab: A rational design of a digital parabolic scaling algorithm[J]. SIAM Journal on Imaging Sciences, 2012, 5(4): 1291-1332. [本文引用:1]