2018, Vol. 35 
2. 上海卡勒幅磁共振技术有限公司, 上海 200062;
3. 西门子医疗东北亚科研合作部, 上海 201318
2. Shanghai Colorful Magnetic Resonance Technology Co., Ltd., Shanghai 200062, China;
3. MR Collaboration NE Asia, Siemens Healthcare, Shanghai 201318, China
磁共振成像（MRI）在场强较低或采用一些特殊的成像序列（如高b值扩散加权成像）扫描时，图像的信噪比（SNR）往往比较差.在临床应用中，一般通过增加采集次数，将多幅图像进行累加平均以提高图像的信噪比.然而，由于心脏和腹部主动脉的波动、呼吸运动和肠的蠕动等多种因素的影响，多次采集的图像之间存在偏移或形变，直接进行累加平均会导致图像模糊.采用心电门控和呼吸门控的方法，可以在一定程度上缓解这一问题，但并不能完全解决问题.
提高磁共振图像的信噪比是经典的图像去噪问题[1].用于磁共振图像去噪的方法包括基于小波变换的滤波方法[2-4]、基于稀疏性的去噪方法[5, 6]及字典学习的方法[7, 8]等.2005年，Buades等人[9]基于具有相似邻域模式的像素具有相似灰度值的假设，提出了非局域均值去噪算法（NLM）.由于这种算法操作简单，在降低图像噪声的同时对图像细节影响比较小，因此得到了广泛的使用[10-16].
针对低信噪比图像累加平均变模糊的问题，本研究组曾基于各帧图像在发生偏移过程中，图像的局部结构仍保持不变的假设，利用NLM算法计算图像的局部相似性方法，对图像间的平移进行修正后再进行加权平均，以提高图像的信噪比，即匹配加权平均（MWA）方法[17].但该方法受局部模式相似性计算方法的限制，不能处理图像之间存在局部旋转的图像序列.在此基础上，本文提出一种旋转不变的非局域均值（RINLM）算法，利用圆形邻域区域取代方形邻域区域，并将其以中心像素划分为一系列等面积的圆环，模式的相似性不再受旋转因素影响，充分利用图像中发生旋转的相似模式，提高算法的去噪效果.将RINLM算法应用于图像序列的平均和去噪中，可以克服局部运动的旋转成分对图像局部偏移量校正的影响，进一步提高图像的信噪比.我们利用模拟和真实的磁共振图像进行实验，将本文提出的RINLM算法与NLM算法，以及结合RINLM的MWA方法（MWA+RINLM）与MWA方法进行了比较.
1 理论与方法 1.1 NLM算法NLM算法[9]利用了图像中大量存在的冗余信息，假设具有相似邻域区域的像素具有相似的灰度值，利用整幅图像中所有像素加权平均得到目标像素的估计值.每个像素对应的权重，由该像素对应邻域区域与目标邻域区域之间的相似性大小决定.
假设实际采集到的带噪图像为v，带噪图像中的噪声为加性噪声.则NLM算法对带噪图像中像素i的估计值NL(i)为
| $ NL(i) = \sum\nolimits_{j \in v} {w(i, j)} v(j) $ | (1) |
| $ w(i, j) = \frac{1}{{Z(i)}}{{\rm{e}}^{ - \frac{{||v({N_i}) - v({N_j})||_{2, \alpha }^2}}{{{h^2}}}}} $ | (2) |
| $ Z(i) = \sum\nolimits_{j \in v} {{{\rm{e}}^{ - \frac{{||v({N_i}) - v({N_j})||_{2, \alpha }^2}}{{{h^2}}}}}} $ | (3) |
其中，w(i, j)为像素j对像素i的权重，取决于两个像素邻域区域之间的相似程度；Ni和Nj分别为以像素i和j为中心的正方形邻域区域；α表示在邻域区域上所施加高斯核的标准差，起到突出中心像素的作用；h为指数函数的衰减因子，控制指数函数的衰减快慢；Z(i)为归一化系数.
1.2 RINLM算法在NLM算法中，邻域区域之间的相似程度取决于两者之间加权欧几里德距离的大小，当邻域区域发生相对旋转时，邻域区域之间的加权欧几里德距离就不能真实反映两者的相似程度.所以，NLM算法不能有效利用图像中发生旋转的相似邻域区域，不适用于局部旋转现象较多的图像处理.本文选取圆形邻域区域代替原来的正方形邻域区域，并将其划分为一系列以中心像素为圆心的等面积圆环，每个圆环代表一个环形体素，再计算邻域模式之间的相似性，使算法不受旋转的影响，更真实地反映邻域区域之间的相似性.
如图 1所示，将以像素i为中心的邻域区域Ci划分为s个等面积的圆环，这里以s=9为例.邻域区域的第n个圆环对应的灰度值为：
| $ {u_i}(n) = \sum\nolimits_{m \in {N_i}} {{R_{n, m}} \cdot {v_i}(m)} $ | (4) |
|
图 1 以像素i为中心的圆形邻域区域Ci Figure 1 Circular neighborhood Ci |
其中，vi(m)表示像素i的邻域中像素m的灰度值；Rn, m表示像素m对第n个圆环的权重，具体计算方法为
| $ {R_{n, m}} = \frac{{{S_{n, m}}}}{{\sum\nolimits_{m \in {N_i}} {{S_{n, m}}} }} $ | (5) |
其中，Sn, m表示第n个圆环与像素m重叠的面积，如图 1中网格标记的面积.重叠面积越大，对应的权重值就越大.
RINLM算法中，像素j对像素i的权重计算方法变为
| $ w(i, j) = \frac{1}{{{Z_i}}}{{\rm{e}}^{ - \frac{{||v({C_i}) - v({C_j})||_{2, \alpha }^2}}{{{h^2}}}}} $ | (6) |
| $ Z(i) = \sum\nolimits_{j \in v} {{{\rm{e}}^{ - \frac{{||v({C_i}) - v({C_j})||_{2, \alpha }^2}}{{{h^2}}}}}} $ | (7) |
其中，α为所施加高斯核的标准差，为了与方形邻域区域高斯核的施加方法保持一致，这里，每个圆环上所施加的权重为其内环半径对应的高斯函数值.
2 结果与讨论 2.1 邻域区域之间的相似性为了RINLM算法选择相似邻域模式的正确性，选用256×256的Shepp-Logan图像，选取在图像中处于两种组织边界、位置为（100, 75）的像素作为参考像素，并以其为中心选取边长为5的正方形邻域区域，作为参考邻域区域[图 2(a)].图 2(b)和2(c)分别为采用(2)式和(6)式计算得到的图像中每个像素邻域区域与参考像素邻域区域之间的相似程度，计算值越大，表示相似程度越高，在图中显示为更亮.可以看出，与NLM算法相比，RINLM算法可以搜索到图像中发生相对旋转的相似邻域区域.
|
图 2 (a)无噪声Shepp-Logan图像，图中白色方框表示选取参考邻域区域；(b)利用NLM算法计算的无噪声图的相似性；(c)利用RINLM算法计算的无噪声图的相似性；(d)带噪Shepp-Logan图像，图中白色方框表示选取参考邻域区域；(e)利用NLM算法计算的带噪图的相似性；(f)利用RINLM算法计算出的带噪图的相似性 Figure 2 (a) Shepp-Logan image without noise. The white box is the reference neighborhood; (b) Similarity distribution of Shepp-Logan image without noise treated with NLM; (c) Similarity distribution of Shepp-Logan image without noise treated with RINLM; (d) Shepp-Logan image with noise. The white box is reference neighborhood; (e) Similarity distribution of Shepp-Logan image with noise treated with NLM; (f) Similarity distribution of Shepp-Logan image with noise treated with RINLM |
为了进一步验证RINLM算法在带噪图像中的表现，在图 2(a)上添加均值为0、标准差为0.03的高斯噪声，模拟带噪图像.选择相同位置的像素作为参考像素，以其为中心选取相同大小参考邻域区域[图 2(d)].图 2(e)和2(f)图分别为利用(2)式和(6)式计算得到的图像中每个像素邻域区域与参考像素邻域区域之间的相似性，计算的值越大，表示相似程度越高，在图中显示为更亮.结果显示，带噪图像中，与NLM算法相比，RINLM算法也可以寻找到带噪图像中发生相对旋转的相似邻域区域.
2.2 单幅图像去噪NLM算法中，邻域区域和搜索窗的大小均与图像中的噪声水平有关[18-20]，同样，在RINLM算法中，图像中噪声水平越高，需要更大的邻域区域使其足够强健，同时需要更大的搜索窗，以保证可以找到足够多的相似邻域模式.邻域区域所施加高斯核的标准差和圆环数目需要与邻域区域的尺寸相匹配.衰减因子的取值与图像噪声水平成正比，不同的图像取值不同.
为了验证RINLM算法对单幅带噪图像的去噪效果，首先利用BrainWeb[21]中提供的无噪声T1加权磁共振脑图（大小为217×181）作为原始图像[图 3(a)].在其上面添加均值为0、标准差为0.03的高斯噪声，模拟带噪图像[图 3(d)].图 3(b)为对图 3(d)使用NLM算法去噪的结果（邻域区域边长为5，高斯加权核标准差为1.5，衰减因子为0.033，搜索窗边长为21）；图 3(c)图为对图 3(d)使用RINLM算法去噪的结果（最佳参数：邻域区域直径为5，高斯加权核标准差为1.5，衰减因子为0.021，搜索窗边长为15，圆环数目为11）；图 3(e)和3(f)图分别为图 3(b)和3(c)与原始图像之间的5倍差异图.可以看出，图 3(e)中包含大量明显的边缘细节信息；与之相比，图 3(f)中虽然也包含部分结构信息，但边缘部分不明显.说明RINLM算法由于可以利用图像中发生旋转的相似邻域区域，在去噪过程中可以更好的保持图像的边缘信息.以峰信噪比（PSNR）定量[22]的评价两种算法的去噪效果，RINLM算法去噪结果（PSNR= 37.025 4）优于NLM算法（PSNR= 36.874 3）.
|
图 3 (a)无噪声T1加权磁共振脑图；(b)使用NLM算法去噪处理图像（PSNR= 36.874 3）；(c)使用RINLM算法去噪处理图像（PSNR= 37.025 4）；(d)带噪的T1加权磁共振脑图；(e)和(f)分别表示图(b)和图(c)与图(a)的5倍差异图 Figure 3 (a) Raw T1 weighted MR brain image without noise; (b) Image treated with NLM (PSNR=36.874 3); (c) Image treated with RINLM (PSNR=37.025 4); (d) T1 weighted MR brain image with noise; (e) 5 times residual of image treated with NLM; (f) 5 times residual of image treated with RINLM |
为了验证RINLM算法在不同噪声水平下，对不同类型带噪图像的去噪效果，利用BrainWeb[21]中提供的无噪声T1、T2和质子密度（PD）加权磁共振脑图（大小为217×181）分别作为原始图像.在其上面分别添加标准差为0.01、0.03、0.05、0.07和0.09的0均值高斯噪声，模拟不同噪声水平的带噪图像.图 4表示使用NLM和RINLM算法对不同噪声水平的T1、T2和PD加权带噪图像去噪结果的PSNR.图 4中第1幅和第3幅图显示，随着图像中噪声的增大，NLM和RINLM两种算法处理结果之间的差异逐渐减小，但RINLM算法去噪结果整体优于NLM算法.第2幅图显示，当图像中噪声水平较低时，RINLM算法去噪结果远好于NLM算法；当图像中噪声大于0.05时，NLM算法去噪结果优于RINLM算法，即RINLM算法对低噪声带噪图像去噪效果更好，在去噪过程中能更好的保持图像边缘信息.
|
图 4 使用NLM和RINLM算法对不同噪声水平的T1、T2和PD加权磁共振脑图带噪图像的去噪效果 Figure 4 PSNR of T1, T2 and PD weighted MR images with different noise levels treated with NLM and RINLM |
我们使用的临床真实数据来自1.5 T西门子磁共振仪器Espree，采集序列为扩散加权成像（DWI）序列，观察视野（FOV）= 350 mm×380 mm；层厚（slice thickness）= 8 mm；层数（slice count）= 1；回波时间（TE）= 60 ms；脉冲间隔重复时间（TR）= 350 ms，采样矩阵（matrix）= 98×122.采集的是健康被试的弥散加权肝脏磁共振图像，包括两组数据：一组是在采集过程中添加了自由呼吸心电门控；另一组在采集过程中未加入自由呼吸心电门控.
对于在采集过程中添加了自由呼吸心电门控的弥散加权肝脏图像序列，图 5(a)表示参考图像，这里选择的是第3幅图像；图 5(b)表示累加平均得到的结果；图 5(c)表示使用MWA方法处理结果[17]；图 5(d)表示将RINLM与MWA方法结合应用于图像序列对应的处理结果（最优参数：邻域区域直径为7，圆环数目为19，高斯加权核标准差为2，匹配窗的边长为21，衰减因子为0.02，搜索窗的大小为11，容差为2）；图 5(e)~(h)为图 5(a)~(d)的局部放大图。结果显示，累加平均得到的图像信噪比有所提高，但其中包含的细节结构与参考图像不一致（箭头2）；使用MWA方法处理得到的图像是与参考图像相对应的高信噪比图像，其中包含的细节结构信息与参考图像一致；而使用RINLM与MWA方法结合则可以更好的保留参考图像中的细节结构信息，边缘细节更加锐利，视觉效果也优于前面两种方法.
|
图 5 有自由呼吸门控弥散加权肝脏图像.(a)参考图像；(b)简单相干累加结果；(c)使用MWA方处理结果；(d)使用MWA+RINLM方法处理结果；(e)~(h)为图(a)~(d)的局部放大图 Figure 5 SWI of liver with ECG. (a) Reference image; (b) Image treated with simple average; (c) Image treated with MWA; (d) Image treated with MWA+RINLM; (e)~(h) Partial enlarged details from (a)~(d) |
对于无自由呼吸门控的肝脏弥散加权图像序列，图 6(a)为参考图像，这里选择的是图像序列中的第8幅图像；图 6(b)表示累加平均得到的结果；图 6(c)表示使用MWA算法处理结果[17]；图 6(d)表示RINLM方法结合MWA方法后的处理结果（最优参数：邻域区域直径为7，圆环数目为19，高斯加权核标准差为2，匹配窗的边长为19，衰减因子的为0.03，搜索窗的边长为13，容差为2）；图 6(e)~(h)为图 6(a)~(d)的局部放大图.可以看出，采集过程中不加自由呼吸门控时，图像中噪声比较大.同样，累加平均处理结果中的细节结构与参考图像中不完全一致（箭头2）；MWA方法处理结果中包含的细节结构信息与参考图像一致；相比较而言，本文提出的RINLM算法与MWA方法的结合更好的保留了参考图像中的细节结构信息，边缘细节更加锐利.
|
图 6 无自由呼吸门控肝脏弥散加权图像.(a)参考图像；(b)简单相干累加结果；(c)使用MWA方法处理结果；(d)使用MWA+RINLM方法处理结果；(e)~(h)为图(a)~(d)的局部放大图 Figure 6 SWI of liver (no ECG). (a) Reference image; (b) Image treated with simple average; (c) Image treated with MWA; (d) Image treated with MWA+RINLM; (e)~(h) Partial enlarged details from (a)~(d) |
进行MRI扫描时，如果图像信噪比比较低常常需要进行多次扫描，简单的累加平均不能很好地处理图像间发生运动的情况.我们曾提出通过计算图像间的局部相似性获得图像间的局部偏移量，得到了比简单累加平均更好的结果.但这种方法不能充分利用图像中发生相对旋转的局部相似结构.在此基础上，本文进一步提出RINLM算法，选取圆形邻域区域替代正方形邻域区域，并以中心像素为圆心划分为一系列等面积的圆环，再计算邻域模式之间的相似性.该算法可以更好地利用图像中的冗余信息，提高去噪效果.将RINLM算法应用于图像降噪中，并结合MWA方法应用于图像序列的平均和去噪中，实验结果表明：与NLM算法方法相比，RINLM算法可以进一步提高图像的质量；与MWA方法相比，RINLM与MWA方法的结合可以进一步提高图像的质量.
| [1] | MOHAN J, KRISHNAVENI V, GUO Y H, et al. A survey on the magnetic resonance image denoising methods[J]. Biomed Signal Proces, 2014, 9: 56-69. DOI: 10.1016/j.bspc.2013.10.007. |
| [2] | DONOHO D L, JOHNSTONE I M. Ideal spatial adaptation by wavelet shrinkage[J]. Biometrika, 1994, 81(3): 425-455. DOI: 10.1093/biomet/81.3.425. |
| [3] | DONOHO D L. De-noising by soft-thresholding[J]. IEEE T Inform Theory, 1995, 3(41): 613-627. |
| [4] | PIZURICA A, WINK A M, VANSTEENKISTE E, et al. A review of wavelet denoising in MRI and ultrasound brain imaging[J]. Current Med Imaging Rev, 2006, 2(2): 247-260. DOI: 10.2174/157340506776930665. |
| [5] | MANJON J V, COUPEB P, BUADES A, et al. New methods for MRI denoising based on sparseness and self-similarity[J]. Med Image Anal, 2012, 16(1): 18-27. |
| [6] | BAO L J, ROBINI M, LIU W Y, et al. Structure-adaptive sparse denoising for diffusion-tensor MRI[J]. Med Image Anal, 2013, 17(4): 442-457. DOI: 10.1016/j.media.2013.01.006. |
| [7] | AHARON M, ELAD M, BRUCKSTEIN A. rmK-SVD:an algorithm for designing over complete dictionaries for sparse representation[J]. IEEE T Signal Proces, 2006, 54(11): 4311-4322. DOI: 10.1109/TSP.2006.881199. |
| [8] | ELAD M, AHARON M. Image denoising via sparse and redundant representations over learned dictionaries[J]. IEEE Trans Image Process, 2006, 15(12): 3736-3745. DOI: 10.1109/TIP.2006.881969. |
| [9] | BUADES A, COLL B, MOREL J M, et al. A review of image denoising algorithms, with a new one, Multiscale Model[J]. Siam Journal on Multiscale Modeling & Simulation, 2005, 4(2): 490-530. |
| [10] | MANJÓN J V, COUPÉ P, BUADES A, et al. New methods for MRI denoising based on sparseness and self-similarity[J]. Med Image Anal, 2012, 16(1): 18-27. |
| [11] | CAI B, LIU W, ZHENG Z, et al. A new similarity measure for non-local means denoising[M]. ZHA H, CHEN X, WANG L, et al. Computer vision. communications in computer and information science. Heidelberg: Springer, 2015, 546: 306-316. |
| [12] | LI H J, SUEN C Y. A novel Non-local means image denoising method based on grey theory[J]. Pattern Recogn, 2016, 49: 237-248. DOI: 10.1016/j.patcog.2015.05.028. |
| [13] | PRASATH V B S, KALAVATHI P. Adaptive nonlocal filtering for brain MRI restoration[M]. THAMPI S, BANDYOPADHYAY S, KRISHNAN S, et al. Advances in signal processing and intelligent recognition systems. Switzerland: Springer International Publishing, 2016, 425: 571-580. |
| [14] | MANJON J V, COUPE P, MARTI-BONMATI L, et al. Adaptive non-local means denoising of MR images with spatially varying noise levels[J]. J Magn Reson Imaging, 2010, 31(1): 192-203. DOI: 10.1002/jmri.22003. |
| [15] | COUPE P, YGER P, PRIMA S, et al. An optimized blockwise nonlocal means denoising filter for 3-D magnetic resonance images[J]. IEEE Trans Med Imaging, 2008, 27(4): 425-441. DOI: 10.1109/TMI.2007.906087. |
| [16] | MAHMOUDI M, SAPIRO G. Fast image and video denoising via nonlocal means of similar neighborhoods[J]. IEEE Signal Proc Let, 2005, 12(12): 839-842. DOI: 10.1109/LSP.2005.859509. |
| [17] |
LI W J, XIE H B, YAN X, et al. MR image average based on local offset correction[J].
Chinese J Magn Reson, 2017, 34(3): 294-301.
李文静, 谢海滨, 严序, 等. 基于局部位移校正的磁共振图像相干平均[J]. 波谱学杂志, 2017, 34(3): 294-301. |
| [18] | HE L L, GREENSHIELDS I R. A nonlocal maximum likelihood estimation method for rician noise reduction in MR images[J]. IEEE Trans Med Imaging, 2009, 28(2): 165-172. |
| [19] | SALMON J. On two parameters for denoising with non-local means[J]. IEEE Signal Proc Let, 2010, 17(3): 269-272. DOI: 10.1109/LSP.2009.2038954. |
| [20] | VAN DE VILLE D, KOCHER M. Nonlocal means with dimensionality reduction and sure-based parameter selection[J]. IEEE T Image Process, 2011, 20(9): 2683-2690. DOI: 10.1109/TIP.2011.2121083. |
| [21] | http://brainweb.bic.mni.mcgill.ca/brainweb/selection_normal.html[OL] |
| [22] | WANG Z, SIMONCELLI E P, BOVIK A C, et al. Multiscale structural similarity for image quality assessment[C]. Pacific Grove: Conference Record of the Thirty-Seventh Asilomar Conference on Signals, Systems and Computers, 2004. DOI: 10.1109/ACSSC.2003. 1292216 |