MR图像去噪算法研究进展

王 昊，康晓东，刘玲玲，耿佳佳，郭 军

MR成像技术是利用核磁共振原理对被成像物的整体或局部进行断面或立体成像。现在，无论在临床诊断还是基础研究中，MR成像技术都已成为必不可少的重要工具之一。MRI之所以能被如此广泛的应用，主要是因为它具有2个特点：首先它是一种无创性(noninvasive) 的成像手段，检测过程中不会对被检测对象造成损伤。其次，MRI技术具有灵活多样性和可变通性。在同一次扫描过程中，不同的MR成像和波谱技术可用来获得同一成像体内不同层次、不同类别的信息。与CT、超声成像等传统医学影像技术相比，MRI可以利用多种参数(T1、T2、扩散等)加权成像，能够为临床诊断提供更多、更详细和更准确的信息[1-5]。

尽管MRI有上述优点，然而，由于成像机制的限制，MRI的时间、空间分辨率和信噪比之间存在矛盾[6]。在MR图像中，噪声主要来自硬件电路和被成像对象两个方面，主要是热噪声,另外还有生理学噪声[7]。快速MRI会导致低的信噪比和对比度，因此，在提高图像时间分辨率的时候，图像噪声会大大降低MR图像的质量，使得一些组织的边界变得模糊，细微结构难以辨别，增加了对图像细节识别与分析的难度，影响医学诊断。尽可能地降低图像噪声是所有MR成像设备制造商所追求的目标之一。虽然图像的累加平均能够提高图像的信噪比，却增加了图像采集的时间[8-9]，因此，图像的后处理去噪对于提高图像质量就非常重要。

在MR图像去噪方面，根据实际图像的特点、噪声频谱规律和统计特性，提出了多种图像去噪方法，通常可以分为图像空间域去噪方法和图像变换域去噪方法两大类。图像空间域去噪方法是指直接在空间域上对原始图像的像素点进行处理的技术。它的缺点是对图像中的每一个像素都用相同的方式去进行处理，没有考虑每个像素点自身所固有的特性；优点是算法比较简单，计算速度快。图像变换域去噪方法是先对图像进行某种数学变换，使其从空间域转换到变换域，然后再对变换域中的变换系数进行修正处理，最后进行反变换而得到处理后的图像[10]。

基于上述讨论，笔者将图像去噪算法按照信号域的不同进行划分，从图像空间域去噪算法和图像变换域去噪算法[11]等方面进行综述。

1 空间域去噪算法的研究

空间域去噪方法是指在原始图像上直接进行数据运算，对像素的灰度值进行处理。经典的空间域去噪方法有均值滤波、中值滤波等。这些方法的共同特点是用同一种方式去处理图像中的每一个像素，而不顾每个像素点自身的特性。它们在去除附加的随机噪声方面非常有效，但是在降低噪声的同时，也使图像产生了较为严重的模糊，特别是在图像的边缘和细节处。

空域的方法是直接对这些像素操作的过程，定义公式如下[1]：

式中f(x,y)是输入图像，g(x,y)是处理后的图像，T是对f的一种操作，其定义在点(x,y)的邻域。

1.1 图像平滑

平滑技术用于消除图像中的噪声。经典的平滑噪声方法主要有均值滤波法，中值滤波法等。

均值滤波是空间域线性滤波中最典型的一种滤波技术，也称为线性平滑滤波，其对处理高斯噪声非常有效。均值滤波的思想是对含噪图像中的每一像素点取一个领域，然后计算领域内各像素点灰度值的平均值，作为处理后图像各像素点的值。

均值滤波对噪声的滤除效果与所选取的领域尺寸有关，尺寸越大则去噪的效果越好，但同时也会引起图像边界、轮廓的模糊[3]。

中值滤波是一种能有效抑制噪声的非线性信号处理技术，它最初主要应用于时间序列的分析，其后这种新的滤波思想被引入到了图像处理中，并取得了较好的效果。由于中值滤波是一种非线性的滤波方式，所以具有较好的适应性，使其可以有效地克服线性滤波的一些不足。

中值滤波的基本原理是把数字图像中一点的值用该点领域内各点的中值来代替，设对一个一维序列f1,f2,… ,fi,… 进行中值滤波，取窗口的长度为L(L为奇数)，则其实现的过程为：首先从输入的序列中抽出L个数，然后将这L个数按数值的大小进行排列，最后取中心点的那个数进行输出。用数学表达式表示为：

公式2中yi为像素序列按数值大小排列所确定的中值，Med为取中值操作，式中中括号内为像素序列按数值大小的排序，i∈Z为像素序列属于非零整数集。

若对一个二维图像信号进行中值滤波，其实现过程与一维序列类似。

窗口的选取在中值滤波中具有重要的意义，直接关系到滤波效果的好坏。当图像信息中含有丰富的细节信息时，为了使这些细节信息能够很好的得以保留，会选择较小的滤波窗口。然而在图像去噪中，大的滤波窗口往往会取得好的滤波效果。所以在实际使用窗口时，会选取不同尺寸的窗口，直到得到满意的滤波效果。

空间域维纳(Wiener)滤波是一种重要的滤波方法[12]，它是1949年由Wiener提出来的，是一种最优的线性最小二乘滤波器，广泛应用于预测、估计、内插、信号和噪声的过滤等领域。维纳滤波器是根据有用信号和干扰信号的统计特性(自相关函数或功率谱)，以线性最小均方误差估计准则设计的最优滤波器，它假定噪声是高斯的且是加性的，噪声和信号相互独立，这样能得到最小均方误差意义下的最优滤波。

维纳滤波器的优点是适应面较广，无论平稳随机过程是连续的还是离散的是标量的还是向量的都可应用。对某些问题，还可求出滤波器传递函数的显式解，并进而采用由简单的物理元件组成的网络构成维纳滤波器。维纳滤波器的缺点是，要求得到半无限时间区间内的全部观察数据的条件很难满足，同时它也不能用于噪声为非平稳的随机过程的情况，对于向量情况应用也不方便。因此，维纳滤波在实际问题中应用不多。实现维纳滤波的要求是首先输入过程是广义平稳的；其次输入过程的统计特性是已知的。根据其他最佳准则的滤波器亦有同样要求。然而，由于输入过程取决于外界的信号、干扰环境，这种环境的统计特性常常是未知的、变化的，因而难以满足上述两个要求。这就促使人们研究自适应滤波器。

1.2 图像锐化

锐化处理[13-15]的主要目的是突出图像中的细节或者增强被模糊了的细节。其原理就是和采用空间的微分来完成。微分算子的响应强度域图像与该点(应用了算子) 的突变程度有关。

1.2.1 一阶微分图像去噪(梯度法)

在图像处理中，一阶微分是通过梯度来实现的。函数f(x,y)在其坐标 (x,y)上的梯度模值(梯度) 表示为：

1.2.2 二阶微分算子(拉普拉斯算子)

通过比较一阶和二阶的微分处理的响应，得知二阶微分处理对细节有较强的响应，如细节和孤立点。所以在锐化处理时，为了更加凸显细节，二阶的效果要比一阶的好。一个二元图像f(x,y)的拉普拉斯变换定义为：

2 变换域去噪算法的研究

变换域的去噪方法是另外一类非常有效的去噪方法，其基本思想是：首先对含噪图像进行某种变换，将其从空间域转换到变换域，随后再对变换域中的变换系数进行处理，之后进行反变换，将含噪图像从变换域再转换到原始空间域，以此来达到去除图像噪声的目的[2-4]。

多数MRI系统都是通过正交检波在频域内采集成像对象的数据[8]，再进行反变换之后得到图像数据。关于MR幅度图像中的噪声目前仍存在争议，多数学者认为MR在变换域内的幅度图像中的噪声是Rician噪声，但是，也有些学者认为是高斯白噪[6-8]声。Rician 噪声是一种与信号相关的噪声，在信号强度高的区域表现为高斯噪声，而在信号强度低的区域是Rician噪声。由于MR图像背景区域的信噪比非常低，因此，图像背景区域的噪声是Rician噪声。Rician 噪声不仅可以引起像素值的随机波动，而且可以引入与信号相关的偏差。

所以，很多在空间域中无法进行有效分析的信号，放到变换域中则可以进行有效的分析，其变换后系数的分布特征会较明显，更有利于进行MR图像的去噪处理。常见的将图像从空间域转换到变换域的变换方法很多[16]，如经典的傅立叶变换、小波变换以及近年来发展起来的多尺度几何分析方法等。

2.1 低通滤波[3]

2.1.1 理想低通滤波器

最简单的低通滤波器是“截断”傅立叶变换中所有高频成分，这些成分处在距变换原点的距离比指定的距离D0要远得多的位置。

2.1.2 巴特沃思低通滤波器(BLPF)n级巴特沃思低通滤波器的传递函数的定义为：

式中D(u,v)是 (u,v)点距频率矩形原点的距离，D0为截止频率距原点的距离。

2.1.3 高斯低通滤波器

截止频率距原点为D0的高斯型低通滤波器的传递函数为：

与上面提到的两种滤波器相比，图像会更平滑，即使对微小物体和细小条纹，高斯滤波器也是比较清晰的。

2.2 高通滤波

在频域中进行图像的锐化处理的实质上是加强需要的高频分量。并必须考虑到要在锐化图像的同时抑制噪声。由上面的讨论可知，H(u,v)是滤波器的转移函数，只要选择合适的H(u,v)，让高频分量通过，锐化图像的目的便可以得到。所以称此过程为高通滤波。

2.3 基于小波变换图像去噪方法

长期以来，主要的图像去噪方法往往采用传统的线性滤波以及基于傅里叶变换的去噪方法，但去噪效果不是很理想。随着小波理论的发展和不断完善，其良好的时频特性在信号去噪领域中受到了极大地关注。小波变换进行图像去噪的思想是在小波域中利用相应的规则对含噪图像的小波系数进行处理。含噪图像经过小波变换可以得到图像信号的小波系数和噪声的小波系数，由于这两种系数的特性不同，因此可以在小波域内对这两种系数进行相应的处理，从而使图像信号的小波系数尽可能地被保留，噪声的小波系数尽可能地被剔除。小波图像去噪的步骤是首先对含噪图像进行小波变换，然后在小波域内对系数进行相应的处理，最后对处理后的小波系数进行逆变换[17]。

目前，基于小波变换的图像去噪方法大致可以分为三大类：阈值法、相关法和模极大值法。

2.3.1 阈值法

小波阈值去噪法最早是由Donoho和Johnstone[18]提出，它主要基于以下事实：小波变换是一种去相关变换，经小波变换后，信号能量主要集中于少数幅值较大的小波系数上，而噪声能量则均匀分布于大部分小波系数上，即比较大的小波系数一般以实际信号为主，而比较小的小波系数则很大程度上是噪声。因此，在小波系数中，低频分量中含有大量的有用信息，应该予以保留；同时在高频分量中，一些绝对值较大的重要的小波系数并不是噪声，而是边缘信息，也应予以保留。因此可以通过设定合适的阈值，将绝对值小于阈值的小波系数置零，而保留绝对值大于阈值的小波系数。在小波阈值去噪过程中，阈值和阈值函数的选取是两个非常关键的因素，直接影响到图像的去噪效果。

2.3.2 相关法

信号或图像边缘在各尺度相应位置上的小波系数之间往往具有很强的相关性，而噪声对应的小波系数则具有弱相关或不相关的特点，相关法就是根据此原理来去除噪声的[19]。例如，Xu等[20]提出了一种SSNF(spatiaily selective noise filtration)方法，该方法就是利用了信号和噪声在小波域内的不同相关特性，通过计算和比较相邻尺度间小波系数的相关性大小，来判别小波系数的类型，从而进行取舍的。

2.3.3 模极大值法

模极大值去噪法最早是由Mallat和Hwang[21]提出。由于信号和噪声在小波变换中有着不同的传播特性，即随着尺度的增大，信号所对应的模极大值将增大，而噪声所对应的模极大值将减小。因此，连续做若干次小波变换之后，噪声对应的模极大值已经基本去除或者幅值很小，其余极值点主要由信号控制。模极大值去噪法即是利用上述信号和噪声在不同尺度上模极大值的不同传播特性，从所有小波变换模极大值中选择信号的模极大值而去除噪声的模极大值，然后利用剩余的小波变换模极大值重构出原始信号。

2.4 基于多尺度几何分析图像去噪方法

由于小波变换在低频部分具有较高的频率分辨率和较低的时间分辨率，而在高频部分具有较高的时间分辨率和较低的频率分辨率，因此，小波变换是一种理想的时频分析工具，但由于小波变换在方向性上和各向异性上存在不足，使其不能“最优”的表示二维图像信息，从而需要寻找新的方法来进行表示。也正是这个原因，引起了多尺度几何分析这场变革[22]。

2002年，Do和Vetterli[23]提出了一种“真正”二维图像表示方法，即Contourlet变换。Contourlet变换是一种多方向的、局部的、多分辨率的图像表示方法，其继承了Curvelet变换的一些优点，因此在某种程度上可以认为Contourlet变换是Curvelet变换的另一种实现形式。Contourlet变换的基的支撑区间呈“长条形”结构，并且这种“长条形”结构的长宽比随尺度而变化。Contourlet变换首先用拉普拉斯金字塔滤波器对图像进行多尺度分解，这样可以用来捕捉点奇异，然后用方向滤波器组将同一方向上的奇异点合并成一个系数。Contourlet变换的最后结果实际上是用类似于线段的基结构来描述原图像。2003年Lu和Do[24]提出了一种多尺度、多方向、非冗余的图像表示方法(CRISP-Contourlet变换)。CRISPContourlet变换是从Contourlet变换发展而来的，其思想是利用非可分叠代滤波器组来实现Contourlet变换中方向分析和多尺度分析的分离，且有类似于Contourlet变换的频域划分形式。da Cunha等[25]在2006年提出了一种非下采样Contourlet变换，其主要通过综合非下采样的金字塔分解和非下采样的方向滤波器组来实现，它不但继承了变换多尺度、多方向的特性，同时还具有平移不变性。

基于Contourlet变换的图像去噪算法的主要步骤是首先对含噪图像进行Contourlet变换，然后对变换后的Contourlet系数按照某种规则进行修正，最后对修正后的Contourlet系数进行Contourlet逆变换而得到去噪后的图像。

虽然多尺度几何分析理论还处于完善阶段，但它已经在图像处理的各个领域体现出了极大的优越性。现在人们所提出的多尺度几何分析方法有很多，其主要方法主要有Meyer和Coifiman[26]在1997年提出的Brushlet变换、Candes[27]在 1998年提出的Ridgelet变换，在Ridgelet变换的基础上于1999年提出的Curvelet变换[28]、Pennec和Mallat[29]在2000年提出的Bandelet变换、在2001年提出的Beamlet变换[30]以及Velisavljević等[31]在2004年提出的Directionlet变换等。

3 总结与展望

随着近年来MR成像技术的迅猛发展，带动和促进了医学影像学相关领域的研究和发展，而且MR成像与其他医学成像技术相比，具有对人体无电离辐射损害、软组织高对比度分辨率、提供丰富的诊断信息等优点。虽然扫描技术在空间分辨率、获取速度和信噪比方面有了很大地改进，MR图像依然受到了信号强度不均一、噪声等的影响，这些因素影响了MR图像在医学诊断方面的应用效果。可以看到，MR图像在去噪领域已经取得了很大的进步，但至今还没有一种最好的方法适用于所有情况。所以如何在不损害图像细节信息的前提下，能够有效去除MR图像中所包含的噪声一直是困扰人们的难题，因此MR图像去噪是个不断发展的领域。

首先，传统的MR图像去噪方法，如中值滤波、低通滤波等对去除成像过程中某一特定类型的噪声特别有效，但由于其本身的一些缺陷性和局限性，对于混合噪声的去除可以通过进一步改进，建立诸如有效像素点的统计模型等方法，具有较为深远的意义。

其次，基于小波变换的MR图像去噪及其改进方法一直是研究的热点问题，但如何选取合适的小波基函数，对小波变换后的图像系数建立更加有效的模型、能够更好地刻画不同层之间小波系数的关系都是需要进一步研究的方向。

同时，多尺度几何分析理论目前还处于起始阶段，每类多尺度变换都有自己的特点，其相关理论还很不成熟，并没有一个完整的框架，如何统一这一理论思想用于MR图像的去噪，也是今后一段时期研究的热点方向。

综上所述，MR图像去噪在医学领域具有广泛的应用，它是基于图像在临床诊断和治疗领域重要的基础之一。随着图像处理研究不断深入，在新理论、新工具的支撑下，图像去噪必将更好地服务于图像处理和医学影像的发展。

[References]

[1] Dai JP, Chen HY. MR pulse sequences in the central nervous system: part I. Chin J Magn Reson Imaging, 2010,1(3): 220-226.戴建平, 陈红艳. 磁共振脉冲序列在中枢神经系统中的应用(一). 磁共振成像, 2010, 1(3): 220-226.

[2] Dai JP, Shen HC, Li SW. MR pulse sequences in the central nervous system: part II. Chin J Magn Reson Imaging, 2010, 1(3): 305-310.戴建平, 沈慧聪, 李少武. 磁共振脉冲序列在中枢神经系统中的应用(二). 磁共振成像, 2010, 1(3): 305-310.

[3] Kang XD. Medical imaging processing. Beijing: People’s Medical Publishing House, 2009: 43-60.康晓东.医学影像图像处理.北京:人民卫生出版社,2009: 43-60.

[4] Zhang YJ. Image processing and analysis. Beijing: Higher Education Press, 2002: 94-117.章毓晋.图像处理和分析基础.北京：高等教育出版社,2002: 94-117.

[5] Huang SL, Ye CH, Qiu JQ. Noise reduction of magnetic resonance images based on wavelet transform. Chin J of Magn Reson, 2006, 23(4): 529-541.黄世亮, 叶朝辉, 裘鉴卿. 基于小波变换的磁共振图像去噪声. 波谱学杂志, 2006, 23(4): 529-541.

[6] Wood WL, Bronskill MJ, Mulkern RV, et al. Physical MR desktop data. J Magn Reson Imaging, 1993, 3(Suppl):S19-S24.

[7] Edelstein WA, Glover G, Hardy C, et al. The intrinistic signal-to-noise ratio in NMR imaging. Magn Reson Med,1986, 3(4): 604-618.

[8] Alexander ME, Baumgartner R, Summer AR, et al. A wavelet-based method for improving signal-to-noise ratio and contrast in MR images. Magn Reson Imaging, 2000,18 (2): 169-180.

[9] Sijbers J, den Dekker AJ, der Linden AV, et al. Adaptive anisotropic noise filtering for magnitude MR data. Magn Reson Imaging, 1999, 17(10): 1533-1539.

[10] Yang F. Research on Image denoising based on contourlet transform. Beijing: Beijing Jiaotong University, 2008: 4-5.杨帆.基于Contourlet变换的图像去噪算法研究. 北京：北京交通大学硕士学位论文, 2008: 4-5.

[11] Zhu XC, Liu F, Hu D. Digital image processing and communication. Beijing: Beijing University of Posts and Telecommunications Press, 2002: 58-76.朱秀昌, 刘峰, 胡栋. 数字图像处理与图像通信. 北京：北京邮电大学出版社, 2002: 58-76.

[12] Hamerly G, Elkan C. Altematives to the k2 means algorithm that find better clusterings // Proceedings of the 11th International Conference on Knowledge Management and Knowledge Technologies. New York: ACM Press,2002: 600-607.

[13] Ramponi G. A cubic unsharp masking technique for contrast enhancement. Signal Process, 1998, 67 (2): 211-222．

[14] Polesel A, Ramponi G, Mathews VJ. Image enhancement via adaptive unsharp masking. IEEE Trans Image Process,2000, 9(3):505-510.

[15] Fu SJ, Ruan QQ, Wang WQ, et al. A feature-dependent fuzzy bidirectional flow for adaptive image sharpening.Neurocomputing, 2007, 70 (4-6): 883-895．

[16] Ruan QQ, Peng YH. Digital image processing. Beijing:Publishing House of electronics industry, 2001: 63-211.阮秋琦, 彭玉华. 数字图像处理学. 北京: 电子工业出版社, 2001: 63-211.

[17] Xie JC, Zhang DL, Xu WL. Overview on wavelet image delousing. J Image and Graphics, 2002, 7(3): 209-217.谢杰成, 张大力, 徐文立. 小波图像去噪综述. 中国图像图形学报, 2002, 7(3): 209-217.

[18] Candes EJ. Monoseale ridgelets for the representation of images with edges. California: Department of Statisties,Stanford University, 1999: 1-26.

[19] Johnstone IM, Silverman BW. Wavelet threshold estimators for data with correlated noise. J Royal Statics Society Series, 1997, 59(2): 319-351.

[20] Xu Y, Weaver JB, Healy DM, et al. Wavelet transforms domain filters: a spatially selective noise filtration technique. IEEE Trans Image Process, 1994, 3(6):747-758.

[21] Mallat S, Hwang WL. Singularity detection and processing with wavelets. IEEE Trans, 1992, 38(2):617-643.

[22] Jiao LC, Tan S. Development and prospect of image multiscale geometric analysis. Acta Electronica Sinica,2003, 31(12A): 1975-1981.

[23] Do MN, Vetterli M. Contourlets: a directional multiresolution image representation. Rochester: IEEE International Conference on Image Processingm, 2002:357-360.

[24] Lu Y, Do MN. CRISP-Contourlet:a critically sampled directional multiresolution image representation.SanDiego: ProSPIE Conf on Wavelets X, 2003: 1-10.

[25] da Cunha AL, Zhou J, Do MN. The nonsubsampled contourlet transform: theory, design, and applications.IEEE Trans Image Process, 2006, 15(10): 3089-3101.

[26] Meyer FG, Coifiman RR. Brushlets: a tool for directional image analysis and image compression. Applied and Comput Harmonic Analysis, 1997, 4(2): 147-187.

[27] Candes EJ. Ridgelets: theory and applications. California:Stanford University, 1998: 1-121.

[28] Candes EJ, Donoho DL. Curvelets: a surprisingly effective nonadaptive representation for objects with edges:technical report. California: Stanford University, 2000:1-16.

[29] Pennec EL, Mallat S. Image compression with geometrical wavelets. Vancouver, Canada, September: Proc. of ICIP,2000: 661-664.

[30] Donoho DL, Huo XM. Beamlets and multiscale image analysis. Stanford: Stanford University, 2001: 1-49.

[31] Velisavljević V, Beferull-Lozano B, Vetterli M, et al.Directionlets: anisotropic multidirectional representation with separable filtering. IEEE Trans Image Process, 2006,15(7): 1916-1933.