基于贝叶斯非局部平均滤波的超声图像斑点噪声抑制算法

方宏道，周颖玥，林茂松

(1.西南科技大学 信息工程学院，四川 绵阳 621010； 2.特殊环境机器人技术四川省重点实验室(西南科技大学)，四川 绵阳 621010)

0 引言

超声成像是观察人体内部组织器官(例如腹部、乳房、肝脏、肾脏和肌肉骨骼)的强大技术。与其他医学诊断方法相比，超声成像具有对人体无害、非放射性、低成本、高效等优点[1]；因此，超声成像技术在临床应用上显得更加重要，特别表现在观察孕妇胎儿生长状况和诊断腹部器官损伤这些方面。然而，像所有相干成像方法一样，超声图像容易被斑点噪声[2]污染。斑点噪声是一种粒状图案，由随机分布的结构散射的能量干扰所引起的，其分辨率超过图像系统的能力。由于这种斑点噪声的过多存在，大大降低了超声图像的质量，严重影响了医生对病人的疾病作出正确的判断，因此，有必要研究如何去除超声图像中的斑点噪声。

近些年来，许多去除斑点噪声的方法已经被广泛研究。这些方法可以简单地分为三类，即局部自适应滤波、基于偏微分方程滤波和小波滤波。常用的自适应滤波，例如Lee滤波[3]、Kuan滤波[4]、Fronst滤波[5]、中值滤波器[6-7]等，这些方法都是将斑点噪声当作乘性噪声对待的。在20世纪90年代初，Lopes等[8]建议通过对像素进行分类来改进Lee滤波和Frost滤波，以便于对不同的类采用特定的处理。基于这一思想，自适应斑点减少(Adaptive Speckle Reduction, ASR)滤波器[9]就是利用局部图像统计来确定特定区域被进一步处理。虽然局部自适应滤波器具有较低的算法复杂度，但它们往往会使图像模糊，并且不能得到满意的去噪结果。基于偏微分方程的方法有各向异性扩散(Anisotropic Diffusion, AD)滤波[10-11]和全变分最小模型(Total Variation, TV)[12]，这些方法都可以用来进行超声图像去噪处理。改进的各向异性扩散(Speckle Reducing Anisotropic Diffusion, SRAD)[13-14]是一个非线性迭代滤波器，推导出瞬时方差系数作为边缘检测器。与前面方法相比，基于偏微分方程的方法在保护图像边缘的同时，对图像进行了平滑处理。然而，有意义的结构细节在迭代的过程中不幸被去除。一些基于小波的方法[15]开始用于斑点减少的研究，基于这一假设，斑点噪声可以用对数运算转换为加性高斯噪声。然而，作为非线性运算，对数变换不仅完全改变了斑点的统计，而且偏离了噪声的平均值。为了放宽这种限制性假设，首先利用小波系数构建基本边缘图，然后使用前向和后向扩散技术进行图像去噪[16]。Achim等[17]利用贝叶斯框架来确定适应于非高斯信号的统计的小波阈值。Khare等[18]提出使用复尺度小波系数的虚分量来检测强边，然后在非边缘点处对小波系数的幅度进行收缩。另外，基于稀疏表示的图像去噪技术在图像高斯噪声的去除中表现很好[19]，根据噪声图像所训练出的自适应过完备字典能更加稀疏地表示图像块信号，使得去噪效果进一步提升。随着深度学习技术在信号处理领域的广泛应用，一些学者尝试利用神经网络、深度学习的方法来去除图像中的高斯噪声[20]，也取得了非常不错的效果。

非局部平均(Non-Local Means, NLM)滤波算法[21-23]是一种简单有效的算法。该算法充分考虑了图像中充满着丰富的冗余信息，采用非局部加权平均的方式对图像去噪，在加性高斯噪声的去除上取得了超越传统算法的结果。为了将该算法应用于非高斯的噪声模型，有研究提出了一种贝叶斯非局部平均滤波器[24]，极大地扩展了原算法在其他噪声模型上的应用。本文以被斑点噪声污染的超声图像为研究对象，结合斑点噪声的模型，利用贝叶斯公式推导出两种新的统计距离来比较图像块的相似性，并用于非局部平均滤波器中，提高了算法对于噪声模型的鲁棒性。为了降低算法的复杂度，本文采用了块预选择的方法来选择最相关的图像块，从而避免无用权重的计算，在一定程度上提升了去噪效果，并且加速了算法。另外，在参数的设置上，本文根据多次实验总结出一种定量估算滤波参数h最优值的方法。在算法实现方面，本文基于VS(Visual Studio)以及OpenCV(Open source Computer Vision library)开发了超声图像去斑应用程序，并有良好的用户操作界面，每幅图像的去噪控制在很短时间以内。在实验中，将所提出的方法与其他几个经典方法进行比较，结果表明，该方法不仅具有很强的去斑能力，而且保持了图像细节，如边缘等，从而提高了医学超声图像的质量。

1 斑点噪声模型和非局部平均滤波算法

1.1 斑点噪声模型

在过去的几年中，人们已经提出了一些统计模型来描述超声图像中的斑点噪声。通常认为斑点噪声是一种乘性噪声，并且服从瑞利分布[2]。然而，为了提高图像的质量，需要在输出之前对射频(Radio Frequency, RF)信号进行一系列的标准处理，例如非线性放大、对数压缩、低通滤波、插值等运算，这些操作可能改变了原始信号的统计特性。最近的研究[7]表明，斑点噪声近似地满足Gamma分布或Fisher-Tippett分布。目前，最通用的斑点噪声模型如下：

u(x)=v(x)+vγ(x)η(x)

(1)

其中：u(x)和v(x)分别代表着观测图像和原始图像，η(x)～N(0,σ2)是均值为0的高斯噪声。因数γ的值依赖于超声装置和图像形成的附加处理，通常情况下，取为0.5比较适合[7]。

1.2 非局部平均滤波技术

假设一个灰度级噪声图像u=(u(xi))xi∈Ω定义在一个有界域Ω⊂R2上，u(xi)∈R+是像素点xi∈Ω的噪声观察灰度值。在非局部平均滤波原始的方法[25]中，每个像素xi的修复值NL(u)(xi)是通过对图像中所有像素的灰度值加权平均得到的：

(2)

(3)

NLM滤波算法的基本思想是自然图像中都包含有重复的结构信息。像素与像素之间存在着相关性，而这种相关性并不局限于邻域内，而是涉及到整幅图像；并且，与其他去噪方法不同，NLM算法是将图像块之间进行比较，而不是基于单个像素的对比。因此，它在去除高斯噪声和保护图像边缘上表现良好，但这种方法却不能直接应用于去除斑点噪声，原因在于斑点噪声和高斯噪声有很大不同。下面就来探讨如何将非局部平均滤波器应用于斑点噪声的去除中。

2 基于贝叶斯NLM的超声图像去斑算法

2.1 贝叶斯NLM

(4)

(5)

其中:p(u(N)|v(N))表示在v(N)条件下u(N)的分布。因为p(u(N)|v(N))和p(u(N))不能够仅从图像的一个场景中估计，换个思路，这些概率密度函数可以从图像块Ni的搜索窗口Δi中获取另一个图像块Nj进行估计。根据文献[24]和[28]，接下来对上式分子分母的近似值作转变，分别如下所示：

(6)

(7)

(8)

其中:p(u(Ni)|v(Nj))表示在干净未知的图像块v(Nj)的条件下u(Ni)的概率密度函数。由于真实值v(Nj)是未知的，所以用观测值u(Nj)来代替。因此，得到贝叶斯非局部平均滤波模型：

(9)

2.2 Pearson距离

基于贝叶斯NLM式(9)和斑点噪声模型(1)，对每一个像素，有：

u(x)|v(x)～N(v(x),v(x)2γσ2)

(10)

则：

(11)

给定一个图像块Ni，P为该图像块包含的像素数目，概率计算如下(条件独立性假设):

(12)

因此，可引出斑点噪声干扰下图像块之间的距离为：

(13)

将dPearson代替式(3)中的高斯加权欧氏距离，即为适应于斑点噪声模型的NLM，取名PearsonNLM(Pearson Non-Local Means)。

2.3 根距离

(14)

(15)

则可引出一个新的距离dsqrt用来描述受到斑点噪声污染的图像块之间的差异，定义如下：

dsqrt(u(Ni),u(Nj))=

(16)

因为v是未知的，一个简单但经典的估计是将v替换成u，这样最终距离公式变为：

(17)

将dsqrt代替式(3)中的高斯加权欧氏距离，即为另一种适应于斑点噪声的NLM，取名SqrtNLM(Sqrt Non-Local Means)。

2.4 快速算法

NLM滤波器已经被证明是从图像中去除噪声的最有效技术之一，然而，由于计算负担的原因阻碍了它的广泛应用。如果待处理的图像大小为I×I，相应的匹配区域半径为d，搜索区域的半径为M，则NLM算法的整个计算复杂度将为I2×(2d+1)2×(2M+1)2。对于正在处理的每个像素，搜索区域中的所有像素都要被计算，并且相应邻域之间的差异也要被计算。为了加速算法，本文采用了图像块预选择[29]的方法。在搜索区域Δi中对像素进行预选择，以此避免无用权重的计算。其主要思想是在Δi中选择权重高的像素xj，而不是使用搜索区域中的所有像素去计算距离。在本文中，使用u(Ni)和u(Nj)的平均值比值和方差比值进行预选择，那些比值在固定范围内的图像块才用于权重计算。预选择方式如下所示：

w(xi,xj)=

(18)

其中:μ1和σ12是两个固定阈值，控制着最相关图像块的选择。

通过这种预选择，不仅对算法进行了加速，还能使NLM滤波器更好地保护细节区域。因为在平坦区域中，由于存在大量相似的像素，这时增加像素的数量倾向于改善去噪；与之相反，在杂乱的区域，因为有很少相似的像素，所以像素数量的增加往往会在平滑过程中消除细节。

另外，为了加速算法，本文基于VS和OpenCV对算法进行了实现，将处理一幅超声斑点噪声图片的时间控制在很短时间之内。在实验中，对于所有基于NLM的方法，如果将匹配窗口和搜索窗口分别设为3×3和17×17，则PearsonNLM方法用时1.77 s，而SqrtNLM方法由于根号操作的存在时间大约31.03 s，并有良好的用户操作界面,如图1所示。

图1 超声图像恢复界面 Fig. 1 Ultrasonic image recovery interface

3 实验结果与分析

图2 滤波参数h=Cσ对PSNR的影响 Fig. 2 Influence of smoothing parameter h=Cσ on PSNR

3.1 幻影图像去噪实验

在这个实验中，幻影图像如图3(a)所示，它是由Matlab软件生成的，是8位的灰度图像，大小为256×256。斑点噪声的模拟大致如下：首先产生和该图像大小一样的实部和虚部都是均值为0方差为σ2的复数矩阵；然后对该复数矩阵使用3×3的滑动窗口进行均值滤波；最后使用式(1)的噪声模型对生成的幻影图像进行加噪处理来模拟超声图像的斑点噪声，这里γ取值为0.5。在污染图像时，本文设置了四种噪声级别σ={5,10,15,20}进行测试，对比算法的参数设置如表1所示。

为了评估去噪算法，采用两种评价指标。一是峰值信噪比(Peak Signal-to-Noise Ratio, PSNR)[22]值，单位dB，用式(9)来计算：

(19)

表2列出了所有评估方法的PSNR和MSSIM值，这些实验都是在不同程度的斑点噪声污染的幻影图像上进行的，对应σ=10的噪声图像和去噪图像如图3所示。从该实验中，可以看出基于两种新距离的NLM去噪算法在去除斑点噪声上都有很好的表现，不仅在去噪图像的PSNR值和MSSIM值上有很大的提升，同时在视觉上保留了有意义的边缘细节。

表1 用于模拟斑点图像的参数最佳设置Tab. 1 Optimal parameter setting for speckle image simulation

表2 γ=0.5时不同噪声水平下的去斑方法的性能比较Tab. 2 Performance comparison for despeckling methods under different levels of noise when γ=0.5

图3 模拟斑点图像不同滤波器的去斑结果(σ=10) Fig. 3 Despeckling results obtained with different filters for simulated speckle image when σ=10

3.2 真实超声图像去噪实验

为了评估算法对于真实超声图像去噪的效果，进行了如图4和图5所示的实验。由于没有任何本真图像，因此很难用数字来评估真实超声图像的去斑效果。在视觉上，基于两种新距离的NLM滤波器有效地去除了斑点成分，同时增强边缘并保留了图像的结构细节。而其他的去噪方法要么斑点噪声去除得不干净，要么使图片过度平滑，一些重要细节丢失,这些缺陷都不利于医生对病情的诊断。

图4 不同滤波器对体内显示超声图像的去斑结果 Fig. 4 Despeckling results obtained with different filters for in vivo displayed ultrasonic image

图5 不同滤波器对真实超声图像的去斑结果 Fig. 5 Despeckling results obtained with different filters for real ultrasonic image

4 结语

在本文中，通过贝叶斯非局部平均滤波模型引入Pearson距离和根距离进行图像块相似度的计算，从而得到适用于滤除超声图像斑点噪声的NLM滤波器。同时，本文采用了图像块预选择的方案来选择最相关的图像块进行权重计算，这在一定程度上提升了去噪效果。由于本文算法是在VS2013平台下用OpenCV视觉库实现的，所以算法运行的速度比在Matlab平台下实现的算法快得多。为测试本文算法的去噪性能，分别在幻影图像和真实超声图像上进行了实验。实验结果表明，相比一些现有的方法，基于两种距离的NLM算法在去除斑点噪声上都有很好的表现，同时保留了重要的结构细节，这无疑有利于图像后处理工作(如图像配准或图像分割)进一步得到展开。

参考文献(References)

[1] LOIZOU C P, PATTICHIS C S, CHRISTODOULOU C I, et al. Comparative evaluation of despeckle filtering in ultrasound imaging of the carotid artery [J]. IEEE Transactions on Ultrasonics Ferroelectrics and Frequency Control, 2005, 52(10):1653-1669.

[2] 刘春明,张相芬,陈武凡.基于小波的医学超声图像斑点噪声抑制方法[J].中国医学物理学杂志,2006,23(5):364-367.(LIU C M, ZHANG X F, CHEN W F. Wavelet-based method for speckle reduction in medical ultrasound image [J]. Chinese Journal of Medical Physics, 2006, 23(5): 364-367.)

[3] LEE J S. Digital image enhancement and noise filtering by use of local statistics [J]. IEEE Transactions on Pattern Analysis and Machine Intelligence, 1980, 2(2):165-168.

[4] FROST V S, STILES J A, SHANMUGAN K S, et al. A model for radar images and its application to adaptive digital filtering of multiplicative noise [J]. IEEE Transactions on Pattern Analysis and Machine Intelligence, 1982, 4(2):157-166.

[5] KUAN D T, SAWCHUK A A, STRAND T C, et al. Adaptive noise smoothing filter for images with signal-dependent noise [J]. IEEE Transactions on Pattern Analysis and Machine Intelligence, 1985, 7(2):165-177.

[6] 赵高长,张磊,武风波.改进的中值滤波算法在图像去噪中的应用[J].应用光学,2011,32(4):678-682.(ZHAO G C, ZHANG L, WU F B. Application of improved median filtering algorithm to image de-noising [J]. Journal of Applied Optics, 2011, 32(4):678-682.)

[7] 胡静波.改进的中值滤波去噪算法分析[J].信息技术,2011(8):32-33.(HU J B. Analysis of improved median filtering de-noising algorithm [J]. Information Technology,2011(8):32-33.)

[8] LOPES A, TOUZI R, NEZRY E. Adaptive speckle filters and scene heterogeneity [J]. IEEE Transactions on Geoscience and Remote Sensing,1990, 28(6): 992-1000.

[9] KARAMAN M, KUTAY M A, BOZDAGI G. An adaptive speckle suppression filter for medical ultrasonic imaging [J]. IEEE Transactions on Medical Imaging, 1995, 14(2): 283-292.

[10] PERONA P, MALIK J. Scale-space and edge detection using anisotropic diffusion [J]. IEEE Transactions on Pattern Analysis and Machine Intelligence, 2002, 12(7): 629-639.

[11] RAMOS-LLORDEN G, VEGAS-SANCHEZ-FERRERO G, MARTIN-FERNANDEZ M, et al. Anisotropic diffusion filter with memory based on speckle statistics for ultrasound images [J]. IEEE Transactions on Image Processing, 2015, 24(1): 345-358.

[12] RUDIN L I, OSHER S, FATEMI E. Nonlinear total variation based noise removal algorithms [J]. Physica D: Nonlinear Phenomena, 1992, 60(1/2/3/4): 259-268.

[13] 洪志强,张立亭,陈竹安,等.基于改进的各向异性扩散图像去噪算法研究[J].江西科学,2017,35(3):433-437.(HONG Z Q, ZHANG L T, CHEN Z A, et al. Research of images denoising based on improved anisotropic diffusion algorithm [J]. Jiangxi Science, 2017, 35(3):433-437.)

[14] YU Y, MOLLOY J A, ACTON S T. Three-dimensional speckle reducing anisotropic diffusion [C]// Proceedings of the 2003 Conference Record of the 37th Asilomar Conference on Signals, Systems and Computers. Piscataway, NJ: IEEE, 2003: 1987-1991.

[15] RABBANI H, VAFADUST M, ABOLMAESUMI P, et al. Speckle noise reduction of medical ultrasound images in complex wavelet domain using mixture priors [J]. IEEE Transactions on Biomedical Engineering, 2008, 55(9): 2152-2160.

[16] NI C, LI Q, XIA L Z. A novel method of infrared image denoising and edge enhancement [J]. Signal Processing, 2008, 88(6): 1606-1614.

[17] ACHIM A, BEZERIANOS A, TSAKALIDES P. Novel Bayesian multiscale method for speckle removal in medical ultrasound images [J]. IEEE Transactions on Medical Imaging, 2001, 20(8): 772-783.

[18] KHARE A, KHARE M, JEONG Y, et al. Despeckling of medical ultrasound images using Daubechies complex wavelet transform [J]. Signal Processing, 2010, 90(2): 428-439.

[19] 赵井坤,周颖玥,林茂松.基于稀疏表示与非局部相似的图像去噪算法[J].计算机应用,2016,36(2):551-555.(ZHAO J K, ZHOU Y Y, LIN M S. Image denoising algorithm based on sparse representation and nonlocal similarity [J]. Journal of Computer Applications, 2016, 36(2):551-555.)

[20] HARMELING S. Image denoising: can plain neural networks compete with BM3D? [C]// Proceedings of the 2012 IEEE Conference on Computer Vision and Pattern Recognition. Washington, DC: IEEE Computer Society, 2012: 2392-2399.

[21] 张兆伦.基于非局部均值图像去噪算法研究[D].南京:南京邮电大学,2015:19-24.(ZHANG Z L. The research on image denoising algorithm based on non-local mean [D]. Nangjing: Nangjing University of Posts and Telecommunications, 2015: 19-24.)

[22] 周颖玥,臧红彬,赵井坤,等.基于非局部平均滤波的冲击噪声图像恢复算法[J].计算机应用研究,2016,33(11):3489-3494.(ZHOU Y Y, ZANG H B, ZHAO J K, et al. Image recovering algorithm for impulse noise based on nonlocal means filter [J]. Application Research of Computers, 2016, 33(11): 3489-3494.)

[23] SUDEEP P V, PALANISAMY P, RAJAN J, et al. Speckle reduction in medical ultrasound images using an unbiased non-local means method [J]. Biomedical Signal Processing and Control, 2016, 28(6): 1-8.

[24] KERVRANN C, BOULANGER J, COUPÉ P. Bayesian non-local means filter, image redundancy and adaptive dictionaries for noise removal [C]// International Conference on Scale Space and Variational Methods in Computer Vision, LNCS 4485. Berlin: Springer, 2007: 520-532.

[25] BUADES A, COLL B, MOREL J M. A review of image denoising algorithms, with a new one [J]. SIAM Journal on Multiscale Modeling and Simulation, 2010, 4(2): 490-530.

[26] BUADES A, COLL B, MOREL J. On image denoising methods [EB/OL]. [2017- 03- 12]. http://www.ece.lsu.edu/gunturk/EE7700/OnImageDenoisingMethods.pdf.

[27] COUPE P, HELLIER P, KERVRANN C, et al. Nonlocal means-based speckle filtering for ultrasound images [J]. IEEE Transactions on Image Processing, 2009, 18(10): 2221-2229.

[28] SPJØTVOLL E. A nonlinear gaussian filter applied to images with discontinuities [J]. Journal of Nonparametric Statistics, 1997, 8(1): 21-43.

[29] MAHMOUDI M, SAPIRO G. Fast image and video denoising via nonlocal means of similar neighborhoods [J]. IEEE Signal Processing Letters, 2005, 12(12): 839-842.

This work is partially supported by the National Natural Science Foundation of China (61401379), the General Project of Educational Commission of Sichuan Province (14ZB0107).

FANGHongdao, born in 1991,M.S. candidate. His research interests include digital image processing.

ZHOUYingyue, born in 1983, Ph.D., associate research fellow. Her research interests include image processing and analysis.

LINMaosong, born in 1964, Ph.D., professor. His research interests include computer graphics and image processing, visualization in scientific computation.