2. 广东省医学图像处理重点实验室,广东 广州 510515
2. Guangdong Provincial Key Laboratory of Medical Imaging Processing, Southern Medical University, Guangzhou 510515, China
Gibbs伪影是在磁共振成像中常见的主要存在于磁共振图像边缘处的一种伪影,它是采用部分k空间数据进行图像重建[1-3]时产生的。在临床应用中,考虑到患者的实际情况,通常会忽略一定数量k空间外围相位编码线以缩短采集时间。对有限的数据使用傅里叶变换进行重建,在重建图像中就会出现Gibbs伪影[4-6],最终影响临床磁共振图像的诊断价值。最近有文献报道Gibbs伪影降低了心内膜下缺血的诊断[7-9]以及引起弥散加权图像中定量分析结果的偏差[10-11]。
目前临床上常用消除Gibbs伪影的用技术[12-13]是图像滤波器,通过削弱k空间的高频信息从而达到消除Gibbs伪影目的。这种滤波技术的优点是速度快,消除效果明显。但滤波器参数的优化[14-17]要根据每种应用类型分别设置其参数。在降低Gibbs伪影的其它方法中,Ferreira等[10,18-21]提出的亚体素移位的方法消除Gibbs伪影效果明显,该方法是通过图像域的全局亚体素移位来近似过零点处重采样振荡模式。但是,优化全局移位可能会影响图像的细节信息。考虑到局部亚像素移位的优化问题,Elias等[22]提出了一种新的方法局部亚像素移位(LSS)方法,它基于亚体素移位的思想,并通过最小化局部总变差消除了Gibbs伪影。尽管LSS方法显著地降低了图像的Gibbs伪影,但该方法不能够消除k空间零填充重建图像中的Gibbs伪影。
在本文中,我们提出了两种将LSS应用于欠采样数据的方法。第1种方法为LSS+插值,该方法首先将非零填充的k空间数据傅里叶变换到图像域,其次在图像域使用LSS算法搜索最佳的亚像素移位,最后将得到图像进行内插获得最终图像。第2种方法首先对k空间数据进行零填充,其次将零填充后的k空间数据傅里叶变换[23]到图像域,然后在图像域通过隔行局部变差(iLV)法来搜索最佳的局部亚像素移位,从而消除图像的Gibbs伪影。这项工作在体模、T2加权(T2WI)数据中均进行了验证,我们发现,第2种方法,也就是隔行局部变化(iLV)的方法在保留图像清晰边缘方面比第1种方法表现更突出。
1 方法 1.1 伪影校正在磁共振成像中,Gibbs伪影通常出现在低分辨率图像上,也可以认为是高分辨率图像的k空间乘以一个狭窄的矩形窗,即,对于|N| > ωk为零。在图像域中,该现象可以被看作是对应图像与‘sinc’函数的卷积。然而在高对比度组织边界,如图像边缘的急剧变化,图像域卷积可以看作是一系列‘sinc’函数的和,并会造成“振铃”形振荡的周期性与幅度衰减。非零填充k空间重建数据的振荡周期是两个像素。一般地,非零填充k空间数据重建的振荡信号可以由公式(1)给出:
$ f[n] = \sum\limits_{\Omega = -\frac{N}{2}}^{\frac{N}{2}-1} F(\Omega) e^{j \frac{\Omega}{N} \pi n} $ | (1) |
f [n]表示在n处的重建数据,F(Ω)表示在Ω处非零填充的k空间数据的截断数据,并且当
$ \begin{array}{l} f{[n]_{m, padded}} = \sum\limits_{\Omega = -\frac{N}{2}m}^{m\frac{N}{2} -1} F (\Omega ){e^{jn\frac{\Omega }{{mN}}2\pi }} = \sum\limits_{\Omega = -\frac{N}{2}}^{\frac{N}{2} - 1} F (\Omega ){e^{jn\frac{\Omega }{{mN}}2\pi }}\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; = \sum\limits_{\Omega = - \frac{N}{2}}^{\frac{N}{2} - 1} F (\Omega ){e^{j\frac{\Omega }{N}2\pi \frac{n}{m}}} \end{array} $ | (2) |
在公式(2)中,f [n]m, padded可以看作是非零填充重建振荡模型f [n]的m倍密集采样。
基于上面的描述,我们提出一种新的iLV将LSS扩展到零填充的情况。对于整数倍的零填充,iLV最初将密集样本分为m个独立的振荡模型,然后在每一个独立的振荡模型中用LSS去搜索最小的局部变差。原始LSS和iLV的比较如图 1A和图 1B所示,我们可以看到振荡周期为4个像素。图中的黑色线表示真实值,红色线表示由50% k空间低频部分重建的Gibbs伪影。在A中绿色的点表示目标点之一,蓝色的点表示在LSS内核中计算的对应点。在B中,橙色点表示目标点之一,蓝色点是在iLV内核中计算的对应点。在C中,对于相同的目标点,在LSS和iLV内核计算中的值分别用绿线和橙线表示,显然iLV可以接近最小值,这表明iLV可以很大程度的降低振荡并减弱环状伪影的影响。
对于非整数倍的零填充,我们首先完成额外的零填充以确保零填充的部分可以被非零部分整除,然后再进行iLV操作,最后我们重新调整重建后图像的矩阵大小。现在以30%的零填充为例,k空间数据包括70%的采集数据与30%的零填充数据,为了确保整数倍的零填充,我们需要进行40%的零填充以确保零填充部分与非零部分是整数比的关系。然后我们将重建的振荡模型分为两个独立的振荡模型,这样LSS就能够使用。最后我们重新调整重建后图像的矩阵大小。iLV的流程图如图 2所示。为了更好的比较,我们提出LSS+插值的策略。首先,我们将k空间70%的采集数据进行重建,其次使用LSS操作,最后在图像域上执行内插。本文中的LSS直接作用在k空间零填充的数据中。本文还与汉明窗方法做对比,汉明窗滤波器[24-26]直接在零填充k空间上执行,然后完成傅里叶逆变换到图像域(图 2)。
常用的汉明窗滤波器由公式(3)给出:
$ (\varpi ) = \left\{ {\begin{array}{*{20}{c}} {1, \varpi \in {\varpi _p}}\\ {0.54-0.46\cos \left( {\frac{{|\varpi |-{\varpi _c}}}{{{\varpi _t}}}} \right), \varpi \in {\varpi _t}}\\ {0, \varpi \in {\varpi _c}} \end{array}} \right. $ | (3) |
如公式(3)所示:ϖp是通带,ϖt是过渡带,ϖc是阻带。在本文中,过渡带是从k空间非零部分的边界开始。当明显的振荡被完全消除并且保留有很好的分辨率时,我们认为通带达到最优。
1.3 LSS与iLV内核设置我们采用Elias[22]论文中的优化参数实现原始LSS算法,即LV内核(kernel k)设置为[1,3],并且亚体素移位在iLV中覆盖一个像素,内核计算每个单独振荡模式中的相邻体素,并且移位范围随振荡周期而变化。
1.4 操作我们分别在一维空间、二维空间进行数字仿真,去验证本文提出的消除Gibbs伪影方法的有效性与可行性。此外我们还采用了活体数据去验证本文提出方法的有效性与实用性。本实验使用了“双三次”插值算法,所有处理过程都在计算机MATLAB软件自编程序中实现。
1.5 二维数字模型对矩阵尺寸大小为100×100的体模进行实验,在50% k空间零填充数据中测试本文提出方法的性能。在体模中手动添加高斯噪声(SNR = 50),我们将没有伪影的图像作为参考图像。我们通过在相位编码方向手动欠采样50% k空间线重建得到伪影图像,参考方法在伪影图像上实现了不同内核k = [1,2]和k = [1,3]的方法,实验结果如图 4所示。为了评估方法的鲁棒性,在50%零填充的情况下使用具有不同噪声水平的数据进行验证,数字模型选取内核为k = [1,2]进行实验(图 5)。
T2WI数据是使用8通道接收头部线圈在3T扫描仪(Achieva, Philips,Netherlands)采用二维的自旋回拨序列采集,其主要参数为:TR/TE = 2000/72 ms,NAS = 2,层厚 = 6 mm,分辨率 = 1.04 mm×1.04 mm。本文提出的方法和参考方法在40%零填充的T2WI数据进行测试(图 6)。为了验证本文所提方法的鲁棒性,分别对60%,45%,30%零填充的T2WI数据进行内核参数k = [1,2]的测试(图 7)。
实验结果如图 3所示,黑色线表示真实值;红色线表示在50%零填充下的Gibbs伪影;绿色线表示通过汉明滤波后的结果,此方法使边缘信息得到很大的平滑;蓝色线表示用原始LSS方法处理Gibbs伪影的结果,这种方法依然残留了很大的伪影;粉色线表示用LSS+插值的方法处理Gibbs伪影的结果,橙色线表示用iLV处理的结果,根据图像可知,LSS+插值地方法在边缘截断处比iLV方法更加平滑。iLV在去除Gibbs伪影的同时,还保存了细节信息。
图 4展示了本文提出的方法与参考方法的对比结果。其结果与一维仿真的结果类似。我们可以清楚地看到汉明窗滤波器在降低Gibbs伪影的同时还给图像带来了严重的模糊,造成图像细节信息丢失,原始LSS方法处理后图像中依然残留很多Gibss伪影,LSS+插值的方法可以很好地去除Gibbs伪影,但在保持图像分辨率方面比iLV方法差(绿色箭头所示)。本文提出的方法在选取不同k参数的情况下,k = [1,3]的结果比k = [1, 2]的结果更平滑。
在图 5中,我们可以看出本文提出的方法在不同噪声水平下都能很好地消除Gibbs伪影,且iLV方法在保留图像锐转变(白色箭头所示)方面比LSS+插值方法更好。
2.3 活体实验结果40%零填充T2WI图像(冠状位)的测试结果如图 6所示,A表示参考图象,B表示Gibbs伪影图像,C表示汉明窗滤波后的图像,D表示直接用LSS方法得到的图像,E表示用LSS+插值的方法得到的图像,F表示用iLV方法得到的图像。从图像可知使用汉明窗滤波使图像变得非常模糊;直接使用LSS方法得到的图像依然还残留很多Gibbs伪影;在图像E和F中,Gibbs伪影被很好的消除了。接下来的两行显示了3个局部放大图,在局部放大图中我们观察到LSS+插值的方法在去除Gibbs伪影的同时也给图像带来了极大的平滑,因此iLV方法在去除Gibbs伪影的同时也可以减轻对图像的平滑,更好的保留图像的细节信息。
T2WI在60%、45%、30%零填充(冠状位)情况下执行LSS+插值与iLV两种算法,实验结果如图 7所示:两种方法都能够很好的降低Gibbs伪影,但iLV在保留图像细节方面优于LSS+插值的算法。内核参数设置为k = [1,2]。
3 讨论因此,我们提出一种新的方法iLV,以消除k空间零填充情况下的Gibbs伪影。该伪影通常发生在实际MRI扫描中,在考虑不影响采集时间与图像对比度的情况下减少相位编码的步数。iLV作为LSS的扩展,通过利用Gibbs伪影的周期性并通过在零填充的情况下用逐步亚体素移位法在‘Sinc’函数过零点处重采样来减少Gibbs伪影。LSS+插值方法是另一种对照比较方法,直接对k空间非零部分进行LSS操作,然后执行图像域的插值操作。一维仿真、体模以及活体T2WI的测试结果表明,本文所提出的方法都优于文献中的LSS,可以有效的消除零填充情况下的Gibbs伪影,并且iLV在保持图像细节方面比LSS+插值方法更占优势。
与iLV方法相比,我们可以看出LSS+插值方法得到的图像更加平滑。从理论上讲,图像平滑可以归因于LSS操作,并且上采样会将图像平滑给放大。因此LSS+插值组合将进一步降低图像的分辨率。相反,iLV处理零填充的数据,图像的平滑程度与非零填充下情况下的LSS相似,然后下采样到所需的矩阵大小,这样既不会造成图像平滑,也不会造成图像锐变。
与汉明窗滤波器相比,本文提出的方法可以在零填充情况下以最小的平滑稳健性有效的消除Gibbs伪影,尽管处理速度比滤波器慢,在计算机上需要耗费几秒钟,但是本文提出的方法比滤波器得到的图像有更少的模糊。由于iLV不是迭代算法,因此在将来可以做到从体素到图像的加速。
根据图 4的结果可知,不同的内核宽度可能会导致不同的结果,较大的内核尺度会导致图像模糊,而较小的内核尺度会导致iLV易受噪声的干扰。实际上,参数k的选择还应考虑Gibbs伪影的干扰与图像细节的大小。当两个急剧过渡离得太近并且伪影是两个周期性振荡模式的叠加时,经常会发生Gibbs伪影干扰。在我们的活体内实验中(图 7),为了方便计算,我们选用Elias[22]论文中优化的内核参数。未来的工作将集中在自动化的内核算法上。
本文中所提出的方法仍然局限于笛卡尔采集方案和对称采集方案。在非对称采集中,Gibbs伪影是周期性振荡模式和非周期性模式的叠加,因此本文所提出的方法不适用,需要进一步研究算法。与SENSE等并行成像方法[27-30]的结合将有利于减少成像时间并减轻对重建图像中Gibbs伪影的影响,这些内容我们将在接下来的工作中进行研究。
[1] |
柴青焕, 苏冠群, 聂生东. 磁共振部分k空间重建算法[J]. 中国医学物理学杂, 2018, 35(5): 537-42. DOI:10.3969/j.issn.1005-202X.2018.05.008 |
[2] |
张煜东, 王水花. 磁共振加速方法[J]. 中国医学物理学杂志, 2014, 31(4): 5015-21. DOI:10.3969/j.issn.1005-202X.2014.04.010 |
[3] |
江贵平, 黄鑫, 冯衍秋, 等. 一种适用于真实人体数据能有效消除Gibbs伪影的MR重建新算法[J]. 计算机学报, 2007, 30(11): 2040-7. DOI:10.3321/j.issn:0254-4164.2007.11.018 |
[4] |
彭莹莹, 张煜, 江贵平. 基于逆扩散方法的MR图像Gibbs伪影校正[J]. 南方医科大学学报, 2010, 30(9): 2074-6. |
[5] |
冯前进, 黄鑫, 冯衍秋. 基于Chebyshev多项式的消除Gibbs伪影的快速算法[J]. 中国图象图形学报, 2006, 8(8): 1132-8. DOI:10.3969/j.issn.1006-8961.2006.08.014 |
[6] |
冯衍秋, 陈武凡, 梁斌, 等. 基于Gibbs随机场与模糊C均值聚类的图像分割新算法[J]. 电子学报, 2004, 32(4): 645-7. DOI:10.3321/j.issn:0372-2112.2004.04.027 |
[7] |
Bella E, Parker DL, Sinusas AJ. On the dark rim artifact in dynamic contrast-enhanced MRI myocardial perfusion studies[J]. Magn Reson Med, 2005, 54(5): 1295-9. DOI:10.1002/(ISSN)1522-2594 |
[8] |
Ferreira P, Gatehouse P, Kellman PA, et al. Variability of myocardial perfusion dark rim Gibbs artifacts due to sub-pixel shifts[J]. J Cardiovasc Magn Reson, 2009, 11(5): 17. |
[9] |
Sharif B, Dharmakumar R, Labounty T, et al. Towards elimination of the Dark-Rim artifact in First-Pass myocardial perfusion MRI: removing gibbs ringing effects using optimized radial imaging[J]. Magn Reson Med, 2014, 72(1): 124-36. DOI:10.1002/mrm.24913 |
[10] |
Perrone D, Aelterman J, Pizurica A, et al. The effect of Gibbs ringing artifacts on measures derived from diffusion MRI[J]. Neuroimage, 2015, 120(5): 441-55. |
[11] |
Veraart J, Fieremans E, Jelescu IO, et al. Gibbs ringing in diffusion MRI[J]. Magn Reson Med, 2016, 76(1): 301-14. DOI:10.1002/mrm.v76.1 |
[12] |
Muzi L, Siderius M, Gebbie J. An analysis of beamforming algorithms for passive Bottom reflectopn-loss estimation[J]. J Acoust SocAm, 2018, 144(5): 3046. DOI:10.1121/1.5080258 |
[13] |
Bakir T, Reeves SJ. A filter design method for minimizing ringing in a region of interest in Mr spectroscopic images[J]. IEEE Trans Med Imaging, 2000, 19(6): 585-600. DOI:10.1109/42.870664 |
[14] |
Archibald R, Gelb A. A method to reduce the Gibbs ringing artifact in MRI scans while keeping tissue boundary integrity[J]. IEEE Trans Med Imaging, 2002, 21(4): 305-19. DOI:10.1109/TMI.2002.1000255 |
[15] |
Zha YE, Huang YL, Sun ZC, et al. Bayesian deconvolution for angular super-resolution in forward-looking scanning radar[J]. Sensors (Basel), 2015, 15(3): 6924-46. DOI:10.3390/s150306924 |
[16] |
Stern AS, Hoch JC. A new approach to compressed sensing for NMR[J]. Magn Reson Chem, 2015, 53(11, SI): 908-12. DOI:10.1002/mrc.4287 |
[17] |
刘平, 冯衍秋, 陈武凡. 基于凸集投影的稀疏磁共振图像重建新算法[J]. 中国医学物理学杂志, 2009, 26(3): 1160-2, 1179. |
[18] |
Gottlieb D, Shu CW. On the gibbs phenomenon and its resolution[J]. SIAM Review, 1997, 39(4): 644-68. DOI:10.1137/S0036144596301390 |
[19] |
Sarra SA. Digital total variation filtering as postprocessing for Chebyshev pseudospectral methods for conservation laws[J]. NumericalAlgorithms, 2006, 41(1): 17-33. |
[20] |
Block KT, Uecker M, Frahm J. Suppression of MRI truncation artifacts using total variation constrained data extrapolation[J]. Int J Biomed Imaging, 2008, 25(5): 184123. |
[21] |
Knoll F, Bredies K, Pock T, et al. Second order total generalized variation (TGV) for MRI[J]. Magn Reson Med, 2011, 65(2): 480-91. DOI:10.1002/mrm.22595 |
[22] |
Kellner E, Dhital B, Kiselev VG, et al. Gibbs-ringing artifact removal based on local subvoxel- shifts[J]. Magn Reson Med, 2016, 76(5): 1574-81. DOI:10.1002/mrm.v76.5 |
[23] |
陈文静, 苏显渝. 提高Fourier变换轮廓术测量精度的新方法[J]. 光电工程, 2002, 29(1): 19-22, 68. DOI:10.3969/j.issn.1003-501X.2002.01.006 |
[24] |
李宏恩, 周晋阳, 崔艳斌. 抗混叠滤波器在医学成像系统中的应用研究和设计[J]. 中国医学物理学杂志, 2016, 33(8): 844-7. DOI:10.3969/j.issn.1005-202X.2016.08.018 |
[25] |
黄河, 任超. 基于各向异性高斯滤波器的视网膜新生血管增强算法[J]. 中国医学物理学杂志, 2015, 32(6): 887-91. DOI:10.3969/j.issn.1005-202X.2015.06.027 |
[26] |
王红敏, 聂生东, 王远军. 低场核磁共振信号降噪方法研究进展[J]. 中国医学物理学杂志, 2013, 30(4): 4261-5. DOI:10.3969/j.issn.1005-202X.2013.04.012 |
[27] |
肖智魁, 胡广书. 二维SPACE RIP并行成像技术在三维磁共振成像中的应用[J]. 清华大学学报:自然科学版网络.预览, 2009, 49(3): 439-42. |
[28] |
徐雅洁, 祝祯伟, 田浩然, 等. 非线性梯度场下并行成像技术进展[J]. 中国医学影像技术, 2012, 28(11): 2098-101. |
[29] |
Koolstra K, van Gemert J, Boernert P, et al. Accelerating compressed sensing in parallel imaging reconstructions using an efficient circulant preconditioner for cartesian trajectories[J]. Magn Reson Med, 2019, 81(1): 670-85. DOI:10.1002/mrm.27371 |
[30] |
Zhang J, Feng L, Otazo R, et al. Rapid dynamic contrast-enhanced MRI for small animals at 7T using 3D ultra-short echo time and golden-angle radial sparse parallel MRI[J]. Magn Reson Med, 2019, 81(1): 140-52. DOI:10.1002/mrm.27357 |