2015, Vol. 32文章信息
- 肖 洒, 吕植成, 孙献平, 叶朝辉, 周 欣
- XIAO Sa, LV Zhi-cheng,SUN Xian-ping, YE Chao-hui, ZHOU Xin
- CS-MRI中评价随机欠采样矩阵的新方法
- A New Method for Evaluation of Random Undersampling Matrix in Compressed Sensing-MRI
- 波谱学杂志, 2015, 32(4): 584-595
- Chinese Journal of Magnetic Resonance, 2015, 32(4): 584-595
- http://dx.doi.org/10.11938/cjmr20150404
-
文章历史
- 收稿日期：2015-02-15
- 收修改稿日期：2015-10-22
2. 中国科学院大学，北京 100049
2. University of Chinese Academy of Sciences, Beijing 100049, China
压缩感知(Compressed Sensing，CS)作为一种新兴的采样理论，可以在远小于Nyquist采样频率的条件下，使用欠采样方式采集信号数据，然后通过非线性重建算法很好的重建信号．2006年，Donoho、Candes和陶哲轩等人奠定了压缩感知的基础[4]，在此之后压缩感知在信息技术领域引起了广泛关注和研究．磁共振成像(Magnetic Resonance Imaging，MRI)是一种重要的医学影像技术，具有无放射性，非侵入性等特点，但是和其他影像学方法相比其成像速度较慢、成本较高．压缩感知-磁共振成像(Compressed Sensing-Magnetic Resonance Imaging，CS-MRI)是一种将压缩感知技术与磁共振成像相结合，提高成像速度的新技术．通过CS-MRI技术，不仅可以节省采样时间，从而降低磁共振成像的使用成本；而且可以结合现有磁共振成像技术发展许多新的技术，例如快速磁共振功能成像[5]，动态磁共振成像[6]，超极化气体肺部磁共振成像[7, 8]等．
虽然这项技术还处于起步阶段，但已显示出了良好的发展和应用前景．
2007年，Lustig等人[7, 9]将CS技术用于MRI，在显著降低采样比率的情况下，得到了具有较好重建质量的图像，初步建立了CS-MRI方法的框架．在其提出的框架中，CS-MRI方法的主要步骤为：首先为重建图像选择合适的稀疏变换域，研究表明磁共振医学图像在诸如小波变换、离散余弦变换和有限差分变换域内都有较好的稀疏性[7, 8]，此外，还有一些基于新颖的多尺度变换和自适应的稀疏表示[15]，可以较好的降低图像变换后的稀疏度，明显提高图像质量；随后选取合理的随机欠采样矩阵对磁共振图像的K空间进行欠采样实现加快采样速度的目的，随机欠采样矩阵的选取应遵循其和稀疏变换矩阵相关性最小的原则；最后对采样结果采用非线性共轭梯度重建算法[16, 17]重建图像，可以使用的其它重建算法还包括基追踪算法[15, 21]、迭代加权最小二乘法[7, 8]、迭代软阈值法[25]和内点法[26, 27]等．在所有步骤中，合理选取随机欠采样矩阵是决定重建图像质量好坏的关键一步， 这是由于欠采样方式不同，即使是采样比率相同，重建的误差也有明显的差别．Lustig等给出了一种评价随机欠采样矩阵可能产生伪影大小的方法——点扩散函数(Point Spread Function，PSF)法，通过计算点扩散函数在欠采样情况下可能产生伪影的最大值来评价随机欠采样矩阵．点扩散函数是一种天然的评价非相干性的工具，在CS-MRI中其能够很好的反映出由于欠采样产生的伪影状况[22, 24]．许多 后续的CS-MRI研究工作[1]沿用了这种方法和标准．但是，最大值只能反映伪影的最坏情况，并不能反映出伪影的平均状况和波动程度等统计学信息，而这些统计学信息对于重建算法消除伪影的效果十分重要．因此本文认为将平均值(Mean Value，MV)与标准差(Standard Deviation，SD)两个参数也作为评价参数可能会得到较好的重建图像．本文分别在7.0 T和1.5 T磁共振成像仪上获取了小鼠脑部和人体脑部的MRI数据，首先采用不同百分比的最大小波变换系数对小鼠脑部和人体脑部的图像进行了重建，研究图像在小波变换域内的稀疏度；随后分别采用最大值、平均值和标准差作为评价参数在不同的采样比率下生成不同的随机欠采样矩阵；最后对磁共振图像原始数据进行了重采样和重建，并对不同随机欠采样矩阵所获得的重建图像与全采样图像之间的平均绝度误差(Mean Absolute Error，MAE)和均方误差(Mean Squared Error，MSE)进行了比较和评价．
1 实验理论及方法根据CS理论，信号重建效果的好坏以及可以重建信号的最低采样比率与信号在稀疏域内的稀疏度密切相关，依据CS研究中的经验公式：
其中m为欠采样比率，r为图像稀疏变换后的稀疏度[2]．根据(1)式CS采样比率低于4倍的稀疏度时，重建得到的图像质量将会显著下降，甚至无法重建图像．可以从(1)式中看到采样信号经过稀疏变换后越稀疏，CS重建时需要的采样比率就越低．因此本文首先对磁共振全采样图像在小波变换域内的稀疏度进行了分 析，并在实验中设定随机欠采样矩阵的采样比率分别为<、=和>4倍的稀疏度．
点扩散函数是一种天然的评价非相关性的工具，考虑的是欠采样对图像域中单点脉冲的影响．令Fu为由随机欠采样矩阵决定的欠采样傅里叶变换算子，${e_{ij}}$表示一个在(i,j)位置=1，其他位置都为零的脉冲矩阵，那么
PSF表示在(i,j)位置等于1，其位置都为0的矩阵经过傅里叶欠采样再进行傅里叶逆变换后(i,j)位置对(m,n)位置的影响．在满采样的条件下Fu为完整的傅里叶算子，PSF对图像重建不会有影响，即$PSF{(i,{\rm{ }}j;{\rm{ }}m,{\rm{ }}n)_{i,j \ne m,n}} = 0$，但在欠采样的情况下将会产生影响，即在(i,j)位置之外的其他位置也将出现非0点．这些非零点表现为由于欠采样而产生的伪影．在现有技术中，一般通过$PSF{(i,{\rm{ }}j;{\rm{ }}m,{\rm{ }}n)_{i,j \ne m,n}}$最大值的大小来评价伪影的强弱．CS理论是建立在稀疏分解理论之上的，CS理论出现之前稀疏分解最普遍的用途就是图像降噪[3]．稀疏分解正是利用了图像信号和噪声信号在稀疏域的差别实现降噪功能．
CS-MRI中由于欠采样产生的伪影可以认为是一种噪声，最大值虽然体现了欠采样产生伪影的最坏情况，但是不能完全反映噪声的整体情况．平均值和标准差是衡量噪声状况的两个重要指标：平均值反映了伪影的总体大小；标准差反映了伪影的波动状况．
为了评价欠采样重建图像的质量，本文计算了欠采样图像与满采样图像之间的MAE和MSE．MAE和MSE是数字图像处理中常用的两个图像质量评价参数，在CS- MRI中也有着广泛的应用[4]，可以很好地反映欠采样图像与满采样图像之间的差异，
其中，图像矩阵大小为M × N，${f_{ij}}$代表欠采样图像中位置为(i,j)的像素点的像素值，代表满采样图像中位置为(i,j)的像素点的像素值．
本文首先在7 T小动物磁共振成像仪中获取了小鼠脑部的磁共振图像，扫描参数为：重复时间(TR) = 2 ms，回波时间(TE) = 8.667 ms，视野(FOV) = 24 mm×24 mm，扫描层数为10层，层厚=2 mm．
接着在1.5 T人体磁共振成像仪中获取了人体脑部的磁共振图像，扫描参数为：TR = 5 000 ms，TE = 82 ms，FOV = 348 mm×512 mm，扫描层数为5层，层厚=5 mm．
随后的实验中通过计算点扩散函数，分别以最大值、平均值和标准差为评价参数，通过等概率随机的方式，生成不同采样比率下的最佳随机欠采样矩阵，随机次数为100万次．最后对获取的MRI欠采样数据进行重建，图像重建工具为4核3.4 G、内存8 GB、操作系统Windows 7的台式计算机，编程环境为MATLAB 2010a，稀疏变换域为Daubechies小波变换，重建算法采用非线性共轭梯度算法，优化对象中重建图像的部分傅里叶变换与欠采样K空间数据之间的L2范数的系数为1，重建图像小波变换后的L1范数的系数为0.005，重建图像的TV范数的系数为0.002．
得到CS图像后计算了每一幅欠采样图像与全采样图像之间的MAE和MSE，并对结果进行了分析．
2 结果与讨论 2.1 小鼠脑部和人体脑部图像的稀疏度 2.1.1 小鼠脑部图像稀疏度为了研究小鼠脑部图像在小波变换域的稀疏度，采用不同百分比的最大小波变换系数对小鼠脑部图像进行了重建，重建效果图如图 1(a)所示，重建图像与全采样图像的差值图如图 1(b)所示．由此得出的图像稀疏度将在下一步实验中确定采样百分比．
|
| 图 1 (a) 采用不同百分比的最大小波变换系数重建得到的小鼠脑部图像，(b) 重建图像与全采样图像的差值图 Fig. 1 (a) Images of mice brain reconstructed from different percentages of the largest wavelet transform coefficients,(b) Difference between reconstructed images and fully sampled image |
通过观察可以发现，使用5%的最大小波变换系数重建的图像与全采样图像已经几乎没有差异，所以可以根据(1)式确定采样比率为0.2，同时，添加了采样比率0.1和0.3作为参照，研究采样比率与稀疏度关系对重建质量的影响．
2.1.2 人体脑部图像稀疏度采用不同百分比的最大小波变换系数对人体脑部的图像进行了重建，重建效果图如图 2(a)所示，重建图像与全采样图像的差值图如图 2(b)所示．
|
| 图 2 (a) 采用不同百分比最大小波变换系数重建得到的人体脑部图像，(b) 重建图像与全采样图像的差值图 Fig. 2 (a) Images of human brain reconstructed from different percentages of largest wavelet transform coefficients,(b) Difference between reconstructed images and fully sampled image |
通过观察可以发现，使用11%的最大小波变换系数重建的图像与全采样图像几乎已经没有差异，所以可以根据(1)式确定采样比率为0.45，与小鼠情况类似，添加了采样比率0.3和0.6作为参照，研究采样比率与稀疏度关系对重建质量的影响．
2.2 不同随机欠采样矩阵的CS重建效果 2.2.1 小鼠脑部图像结果对于小鼠脑部MRI数据，在采样比率为0.1的情况下分别采用最大值、平均值和标准差作为评价参数选取随机欠采样矩阵对10层磁共振图像的原始数据进行了重采样和重建（由于图像结果类似，限于篇幅，以下展示的小鼠脑部重建效果图都是第一层的图像），结果如图 3所示．
|
| 图 3 采样比率为0.1：全采样图像、最大值图像、平均值图像和标准差图像（本图及以下所有图中Fully sampled对应全采样图像，Max、MV、SD分别对应以最大值、平均值和标准差为评价参数的重建效果图） Fig. 3 Sampling ratio is 0.1: Fully sampled image,Maximum image,MV image and SD image(In this figure and all the following figures,fully sampled stands for the fully sampled image and Max,MV,SD stands for images reconstructed by using maximum value,mean value and the standard deviation as evaluation parameter respectively) |
分别对10层CS图像计算与全采样图像的MAE和MSE，结果如图 4所示．
|
| 图 4 采样比率为0.1: (a) 重建图像与全采样图像之间的MAE值，(b) 重建图像与全采样图像的之间的MSE值 Fig. 4 Sampling ratio is 0.1: (a) MAE between reconstructed images and fully sampled image,(b) MSE between reconstructed images and fully sampled image |
通过观察采样比率为0.1的小鼠脑部图像结果可知：采用3种评价参数得到的重建图像出现了大片的伪影，导致重建图像难以辨认，无法用于诊断和研究．重建图像与全采样图像之间的MAE和MSE都比较接近且数值较大，这是由于采样比率已经低于稀疏度的4倍界限，MRI的原始数据信息丢失严重，重建图像时产生的伪影过多，难以通过CS方法有效消除．这从一个方面验证了CS重建时要求采样比率不小于4倍稀疏度的规律．
在采样比率为0.2的情况下分别采用最大值、平均值和标准差作为评价参数进行图像重建，结果如图 5所示．
|
| 图 5 采样比率为0.2：全采样图像、最大值图像、平均值图像和标准差图像 Fig. 5 Sampling ratio is 0.2: Fully sampled image,Maximum image,MV image and SD image |
分别对10层CS图像计算与全采样图像的MAE和MSE，结果如图 6所示．
|
| 图 6 采样比率为0.2: (a) 重建结果与全采样图像之间的MAE值，(b) 重建结果与全采样图像的之间的MSE值 Fig. 6 Sampling ratio is 0.2: (a) MAE between reconstructed images and fully sampled image,(b) MSE between reconstructed images and fully sampled image |
在采样比率为0.3的情况下分别采用最大值、平均值和标准差作为评价参数进行图像重建，结果如图 7所示．
|
| 图 7 采样比率为0.3：全采样图像、最大值图像、平均值图像和标准差图像 Fig. 7 Sampling ratio is 0.3: Fully sampled image,Maximum image,MV image and SD image |
分别对10层CS图像计算与全采样图像的MAE和MSE，结果如图 8所示．
|
| 图 8 采样比率为0.3: (a) 重建结果与全采样图像之间的MAE值，(b) 重建结果与全采样图像的之间的MSE值 Fig. 8 Sampling ratio is 0.3: (a) MAE between reconstructed images and fully sampled image,(b) MSE between reconstructed images and fully sampled image |
通过观察采样比率为0.2、0.3重建小鼠脑部图像结果可知：重建图像的清晰程度比采样比率为0.1时有了较大提高，伪影也大幅减小，可以用于诊断和研究，这从另一个方面验证了CS重建时要求采样比率不小于4倍稀疏度的规律．同时在采样比率为0.2和0.3时，采用平均值作为评价参数取得了较好的结果，与以最大值为评价参数相比MAE平均下降了47%，MSE平均下降了71%，而采用标准差为参数的效果在采样率为0.2时最差，在0.3时与平均值接近，表现得并不稳定．实验结果产生的原因可能是，图像重建时对伪影总量的控制是非常重要的，平均值均衡的考虑了图像伪影总量的程度，而最大值只是考虑了伪影最坏的可能情况，对伪影总量的大小没有考虑，类似的标准差虽然考虑了伪影波动的情况但是同样缺乏对伪影总量的控制．
2.2.2 人体脑部图像结果对于人体脑部MRI数据，在采样比率为0.3的情况下分别采用最大值、平均值和标准差作为评价参数选取随机欠采样矩阵对5层磁共振图像的原始数据进行了重采样和重建（由于图像结果类似，限于篇幅，以下展示的人体脑部重建的效果图都是第一层的图像），重建结果如图 9所示．
|
| 图 9 采样比率为0.3：全采样图像、最大值图像、平均值图像和标准差图像 Fig. 9 Sampling ratio is 0.3: Fully sampled image,Maximum image,MV image and SD image |
分别对5层CS图像计算与全采样图像的MAE和MSE，结果如图 10所示．
|
| 图 10 采样比率为0.3: (a) 重建结果与全采样图像之间的MAE值，(b) 重建结果与全采样图像的之间的MSE值 Fig. 10 Sampling ratio is 0.3: (a) MAE between reconstructed images and fully sampled image,(b) MSE between reconstructed images and fully sampled image |
通过观察采样比率为0.3的人体脑部图像结果可知：采用最大值作为评价参数重建的图像与采用平均值和标准差相比，具有更小的MAE和MSE．但是MAE和MSE的数值仍然很大，与小鼠情况类似，这同样是由于采样比率低于稀疏度4倍的界限，过低的采样比率导致伪影过多难以通过CS重建来消除，最终重建图像中出现大片的伪影，无法用于诊断和研究．这在对小鼠脑部重建结果的基础上进一步验证了CS重建时采样比率要求不小于4倍稀疏度的规律．
在采样比率为0.45的情况下分别采用最大值，平均值和标准差作为评价参数进行图像重建，结果如图 11所示．
|
| 图 11 采样比率为0.45：全采样图像、最大值图像、平均值图像和标准差图像 Fig. 11 Sampling ratio is 0.45: Fully sampled image,Maximum image,MV image and SD image |
分别对5层CS图像计算与全采样图像的MAE和MSE，结果如图 12所示．
|
| 图 12 采样比率为0.45: (a) 重建结果与全采样图像之间的MAE值，(b) 重建结果与全采样图像的之间的MSE值 Fig. 12 Sampling ratio is 0.45: (a) MAE between reconstructed images and fully sampled image,(b) MSE between reconstructed images and fully sampled image |
在采样比率为0.6的情况下分别采用最大值、平均值和标准差作为评价参数进行图像重建，结果如图 13所示．
|
| 图 13 采样比率为0.6：全采样图像、最大值图像、平均值图像和标准差图像 Fig. 13 Sampling ratio is 0.6: Fully sampled image,Maximum image,MV image and SD image |
分别对5层CS图像计算与全采样图像的MAE和MSE，结果如图 14所示．
|
| 图 14 采样比率为0.6: (a) 重建结果与全采样图像之间的MAE值，(b) 重建结果与全采样图像的之间的MSE值 Fig. 14 Sampling ratio is 0.6: (a) MAE between reconstructed images and fully sampled image,(b) MSE between reconstructed images and fully sampled image |
通过观察采样比率为0.45、0.6重建的人体脑部图像结果可知：图像较采样比率为0.3时清晰程度大幅提高，大片的伪影也得以消除，在此基础上有利于研究和诊断的进一步开展，进一步验证了要求采样比率不小于4倍稀疏度的规律．其中，采用平均值作为评价参数效果最优，而采用标准差的效果与最大值接近．通过比较3个随机欠采样矩阵的参数发现，采用标准差和最大值挑选出的随机欠采样矩阵计算PSF时的平均值也非常接近，进一步验证了平均值在挑选随机欠采样矩阵中的优势以及控制伪影总量对于重建图像的重要性．分析采用平均值作为评价参数的图像结果时发现，与采用最大值相比，MAE平均下降了58%，MSE平均下降了74%，这与小鼠脑部图像的重建结果相比，MAE和MSE下降更多，这可能是因为人脑比小鼠大脑结构更为复杂，在信号中体现为具有更多的高频信息，因此随机欠采样矩阵对图像重建质量的影响也越大．
3 结论在本文的工作中，验证了小鼠脑部和人体脑部的MRI图像在小波变换域内的稀疏性，以及采样比率不小于4倍稀疏度条件在CS-MRI中的通用性．在现有生成随机欠采样矩阵技术的基础上，提出了平均值和标准差作为随机欠采样矩阵的两种新的统计学评价参数，小鼠脑部和人体脑部的MRI数据重建结果表明，在采样比率≥4倍稀疏度时，以平均值为评价参数的CS-MRI比通过其他两种参数获得的图像质量都更优．
| [1] | Candes E J, Romberg J K, Tao T. Robust uncertainty principles: Exact signal reconstruction from highly incompletefrequency information[J]. IEEE T Inform Theory, 2006, 52(2): 489-509. |
| [2] | Donoho D L. Compressed sensing[J]. IEEE T Inform Theory, 2006, 52(4): 1 289-1 306. |
| [3] | Candes E J, Tao T. Near-optimal signal recovery from random projections: Universal encoding strategies[J]. IEEE T InformTheory, 2006, 52(12): 5 406-5 425. |
| [4] | Wech T, Lemke A, Medway D, et al. Accelerating cine-MR imaging in mouse hearts using compressed sensing[J]. J MagnReson Imaging, 2011, 34(5): 1 072-1 079. |
| [5] | Gao Ming-sheng(高明生), Xie Hai-bin(谢海滨), Yan Xu(严序), et al. Selective dual-direction sequential compressedsensing for dynamic MR imaging(选择性双向顺序压缩感知重建动态磁共振成像)[J]. Chinese J Magn Reson(波谱学杂志), 2013, 30(2): 194-203. |
| [6] | Li Hai-dong(李海东), Zhang Zhi-ying(张智颖), Han Ye-qing(韩叶清), et al. Lung MRI using hyperpolarized gases(超极化气体肺部磁共振成像)[J]. Chinese J Magn Reson(波谱学杂志), 2014, 31(3): 307-320. |
| [7] | Lustig M, Donoho D L, Pauly J M. Sparse MRI: The application of compressed sensing for rapid MR imaging[J]. MagnReson Med, 2007, 58(6): 1 182-1 195. |
| [8] | Lustig M, Donoho D L, Santos J M, et al. Compressed sensing MRI[J]. IEEE Signal Proc Mag, 2008, 25(2): 72-82. |
| [9] | Ajraoui S, Lee K J, Deppe M H, et al. Compressed sensing in hyperpolarized He-3 lung MRI[J]. Magn Reson Med, 2010,63(4): 1 059-1 069. |
| [10] | Qu X B, Zhang W R, Guo D, et al. Iterative thresholding compressed sensing MRI based on contourlet transform[J].Inverse Probl Sci En, 2010, 18(6):737-758. |
| [11] | Ravishankar S, Bresler Y. MR image reconstruction from highly undersampled k-space data by dictionary learning[J].IEEE T Med Imaging, 2011, 30(5): 1 028-1 041. |
| [12] | Qu X B, Guo D, Ning B D, et al. Undersampled MRI reconstruction with patch-based directional wavelets[J]. Magn ResonImaging, 2012, 30(7): 964-977. |
| [13] | Song Y, Zhu Z, Lu Y, et al. Reconstruction of magnetic resonance imaging by three-dimensional dual-dictionary learning[J].Magn Reson Med, 2014, 71(3): 1 285-1 298. |
| [14] | Qu X B, Hou Y K, Lam F, et al. Magnetic resonance image reconstruction from undersampled measurements using apatch-based nonlocal operator[J]. Med Image Anal, 2014, 18(6): 843-856. |
| [15] | Chen S S B, Donoho D L, Saunders M A. Atomic decomposition by basis pursuit[J]. Siam J Sci Comput, 1998, 20(1): 33-61. |
| [16] | Ye J C, Tak S, Han Y, et al. Projection reconstruction MR imaging using FOCUSS[J]. Magn Reson Med, 2007, 57(4): 764-775. |
| [17] | Donoho D L, Elad M, Temlyakov V N. Stable recovery of sparse overcomplete representations in the presence of noise[J].IEEE T Inform Theory, 2006, 52(1): 6-18. |
| [18] | Daubechies I, Defrise M, De Mol C. An iterative thresholding algorithm for linear inverse problems with a sparsityconstraint[J]. Commun Pur Appl Math, 2004, 57(11): 1 413-1 457. |
| [19] | Starck J L, Elad M, Donoho D L. Image decomposition via the combination of sparse representations and a variationalapproach[J]. IEEE T Image Process, 2005, 14(10): 1 570-1 582. |
| [20] | Elad M, Matalon B, Zibulevsky M. Coordinate and subspace optimization methods for linear least squares withnon-quadratic regularization[J]. Appl Comput Harmon Anal, 2007, 23(3): 346-367. |
| [21] | Kim S J, Koh K, Lustig M, et al. An interior-point method for large-scale l(1)-regularized least squares[J]. IEEE J STSP,2007, 1(4): 606-617. |
| [22] | Ajraoui S, Parra-Robles J, Wild J M. Incorporation of prior knowledge in compressed sensing for faster acquisition ofhyperpolarized gas images[J]. Magn Reson Med, 2013, 69(2): 360-369. |
| [23] | Imai H, Kimura A, Hori Y, et al. Hyperpolarized 129Xe lung MRI in spontaneously breathing mice with respiratory gatedfast imaging and its application to pulmonary functional imaging[J]. NMR Biomed, 2011, 24(10): 1 343-1 352. |
| [24] | Pang Y, Zhang X. Interpolated compressed sensing for 2D multiple slice fast MR imaging[J]. PLOS ONE, 2013, 8(2):e56098. |
| [25] | Candes E J, Wakin M B. An introduction to compressive sampling[J]. IEEE Signal Proc Mag, 2008, 25(2): 21-30. |
| [26] | Sun Yu-bao(孙玉宝), Wei Zhi-hui(韦志辉), Wu Min(吴敏), et al. Image poisson denoising using sparse representations(稀疏性正则化的图像泊松去噪算法)[J]. Acta Electronica Sinica(电子学报), 2011, 39(2): 285-290. |
| [27] | Fang L, Li S, Nie Q, et al. Sparsity based denoising of spectral domain optical coherence tomography images[J].Biomedical Opt Express, 2012, 3(5): 927-942. |