CN105631820A - 基于小波变换和三边滤波器的医学超声图像去噪方法 - Google Patents

Abstract

基于小波变换和三边滤波器的医学超声图像去噪方法，包括如下步骤：步骤1)建立医学超声图像模型；步骤2)对第一步得到的对数变换后的图像进行小波分解，得到四个频域(LL¹、LH¹、HL¹和HH¹)。对低频域LL¹继续进行小波分解，再得到四个频域(LL²、LH²、HL²和HH²)；然后重复这个步骤，直到分解最大层数J；步骤3)对每一层的高频部分(LH^j、HL^j和HH^j，j＝1,2,...,J)的小波系数进行阈值法收缩处理；步骤4)利用三边滤波器对最后一层的低频部分(LL^J)中的小波系数做滤波处理；步骤5)作小波逆变换处理，得到去噪后的医学超声图像。

Description

基于小波变换和三边滤波器的医学超声图像去噪方法

技术领域

本发明应用于医学超声图像去噪领域，设计一种适用于医学超声图像的基于小波变换的去噪方法。

背景技术

随着科技的发展，在医学成像领域，超声成像、CT、MRI等成像技术已应用于医学临床诊断中。由于超声成像技术凭借高分辨率，操作简单，即时性等优点快速发展，具有无创、无放射性损害、快捷方便等特性，已经成为一种广泛使用且高度安全的医疗诊断技术。尤其对人体的身体器官检查及肿瘤组织，超声成像技术的使用更为重要。

由于超声成像机理的限制，斑点噪声的存在严重影响了超声图像的质量，导致了超声图像质量较差。斑点噪声的产生是由于超声成像中的基本分辨单元内存在大量的随机散射现象，在图像上表现为空间域内相关的形状各异的小斑点，它将掩盖那些灰度差别很小的图像特征。对于临床医生而言，斑点噪声对他们的准确诊断造成了很大的干扰，特别是对于经验不是很丰富的医生造成的影响更大。因此，从临床应用的角度出发，需要研究去除斑点噪声的算法，为医生做出更准确的诊断提供技术支持，降低人工诊断的风险。

由于医院资源的局限性，特别是医生每天进行人工诊断病人的数量无法满足社会整个阶层的需求，即面临着病人多医生少的情况。因此，各种自动诊断仪器的需求越来越大，自动诊断仪器的出现，一方面可以节约医生资源，另一方面可以方便更多的病人进行诊断。随着当今社会经济的飞跃发展，人们自身健康情况却不容乐观，所以人们对家用型医疗自动诊断仪器的需求也非常大，例如家用超声图像自动诊断仪等。但是超声图像自动诊断仪同样面临着图像质量不高的问题，并且自动诊断仪需要对超声图像做后期的智能分析，如特征提取、边缘检测和图像分类识别等。因此，从自动化诊断技术的角度出发，需要研究去除斑点噪声的方法，为图像的后期智能处理提供技术保障，促进自动诊断技术的发展。

综上所述，研究医学超声图像去噪方法具有非常重要的意义：

(1)提高医学超声图像的质量，改善视觉效果；

(2)方便医生更加准确地针对病灶区域做出判断，降低辅助诊断的风险；

(3)促进超声图像自动化诊断技术的发展，具有不可估量的价值。

在数字图像处理领域，滤波常用来修改或增强图像，对图像的某些特征，如轮廓、边缘、细节和对比度等进行锐化，提高图像的视觉质量。由于抑制斑点噪声具有非常重要的意义，众多科研工作者在此问题上投入了大量的精力。近几十年来出现的医学超声图像去噪方法，可以简单分为5种类型：自适应去噪方法、各向异性扩散去噪方法、非局部均值去噪方法、小波变换去噪方法和混合型去噪方法。通过实验虽然自适应滤波方法的复杂度低，但是往往会模糊图像的细节部分，对于斑点噪声的抑制效果不是很理想。各向异性扩散去噪方法，具有很强的去噪能力，但是结果可能会出现过度平滑的现象。非局部去噪算法对于斑点噪声的抑制效果比较理想，但是这类去噪方法的复杂度较高，不易满足医学超声成像系统的实时性要求，往往用于医学超声图像的后期去噪处理。虽然双边滤波器能够很好地保持边缘信息得到大家的认可，但是由于双边滤波器运算复杂度较高并且存在“梯度失真”现象，并不能满足医学超声成像系统的实时性要求。

发明内容

本发明要克服现有技术的上述缺点，结合斑点噪声的模型的特点和医学超声图像的处理需求提出了一种新的去噪方法，即一种基于小波变换和三边滤波器的医学超声图像去噪方法。

小波变换具有时频分析和多尺度分析等优越性，并已在图像处理领域得到了广泛的应用。在处理加性噪声问题时，小波的去噪效果较好，能够满足一般的产品需求。然而，仅仅利用小波变换的去噪方法对医学超声图像中斑点噪声的抑制效果不好。虽然双边滤波器在去噪的过程中也能很好保持边缘信息，但是其效率较低，运行时间将会很长并且存在“梯度失真”现象，难以用于实时系统。随着图像的分辨率越来越大，这在很大的程度上限制了双边滤波的应用空间。因此本发明利用三边滤波器替换掉双边滤波器，极大提高去噪性能和运算效率。对于三边滤波器，一方面可以去除低频部分的脉冲噪声和斑点噪声，另一方面保持图像边缘细节，并且可以解决梯度失真的问题。具体思路如下：在传统的小波去噪方法的基础之上，根据小波域内超声图像及斑点噪声的统计特性，提出了改进的自适应小波收缩算法，能够更有效地去除高频部分的斑点噪声；然后通过提出的三边滤波器去除低频噪声。

基于小波变换和三边滤波器的医学超声图像去噪方法，包括如下步骤：

步骤1)建立医学超声图像模型

如果认为超声成像系统能够对那些影响声波功率的因素做出恰当的动态补偿，则超声成像系统采集的包络信号由两部分组成，一是有意义的体内组织的反射信号，另一部分是噪声信号。其中噪声信号可分为相乘噪声与相加噪声。相乘噪声与超声信号成像的原理有关，主要来源于随机的散射信号。相加噪声认为是系统噪声，如传感器的噪声等。对于包络线的通用模型模型如下

s(x,y)＝r(x,y)n(x,y)(1)

这里，(x,y)代表2D空间坐标，r(x,y)表示无噪声信号，n(x,y)表示相乘噪声。

为了适应超声成像系统显示屏幕的动态显示范围，对超声成像系统采集到的包络信号进行对数压缩处理。此时相乘的式(1)模型将变为相加的模型，如下

log(s(x,y))＝log(r(x,y))+log(n(x,y))(2)

此时，得到的信号log(s(x,y))即是通常看到的医学超声图像。

步骤2)对第一步得到的对数变换后的图像进行小波分解，得到四个频域(LL¹、LH¹、HL¹和HH¹)。

LL1分量是对原始信号LL0的列和行进行小波分解后得到的低频分量，即一级小波分解后近似部分，它包含了原始图像最多的低频信息；

LH1是一次小波分解后的垂直方向上的高频分量，即它包含了图像水平方向上的近似信息和垂直方向上的边缘等高频信息；

HL1是一次小波分解后的水平方向上的高频分量，即它包含了图像垂直方向上的近似信息和水平方向上的边缘等高频信息；

HH1是一次小波分解后对角方向上的高频分量，即它包含了图像水平和垂直方向上的边缘等高频信息；

对低频域LL¹继续进行小波分解，再得到第二层四个频域(LL²、LH²、HL²和HH²)。然后重复这个步骤，直到分解最大层数J。

由于小波变换是线性变换，因此式(3)模型经过二维离散小波变换后得到下面模型：

其中和分别表示含有噪声图像的小波系数、无噪声图像的小波系数和斑点噪声的小波系数。其中上标j为小波变换的分解层数，下标(l,k)为小波域内的坐标。

经过小波分解后的无噪信号的小波系数服从广义拉普拉斯分布，其概率分布如下

式中，是伽马函数，b为形状参数，a为尺度参数，

同时斑点噪声的小波系数服从零均值高斯分布

式中σ_N为小波域内噪声的标准差。

步骤3)对每一层的高频部分(LH^j、HL^j和HH^j，j＝1,2,...,J)的小波系数进行阈值法收缩处理。

在小波去噪方法中，阈值函数的选择会直接影响到最终的图像去噪结果。当阈值选择较小时，一部分大于该阈值的噪声系数会被当作有用信号保留下来，这就导致去噪后的图像依然存在大量噪声；当阈值选择较大时，会将很多系数很小的有用信息当作噪声而置零，这将使得去噪后的图像变得很平滑，损失很多细节信息。因此选择恰当的小波阈值函数非常重要。

Donoho等人提出了一种典型的阈值选取方法，并且从理论上证明了该阈值与噪声的标准差成正比，改阈值函数又称为统一阈值函数，其公式如下

其中，M即是对应小波域内小波系数的总体个数，σ_n是噪声的标准差。在这种阈值函数中，阈值T受小波系数的个数影响较大，即当M过大时，较大的阈值可能会平滑掉那些系数较小的有用信息。

基于式(7)的基础之上，提出了一种更加适合超声图像的阈值函数，其公式如下

其中，T_j是新的阈值函数，σ_n是噪声的标准差，t_j代表j层的自适应参数。这是种常见的阈值改进的方法，t_j的选取是根据实验决定的，在小波分解后，在不同层分解的小波系数具有不同的分布，由此t_j的选择基于j层的选择，但这种选择不是最佳的，如果适当的选择，所提出的方法将反射更多的优越性。

在小波去噪方法中，首先选定一个给定阈值，然后按照一定的规则对小波系数进行收缩，便完成了对小波系数的去噪。即给定一个阈值，所有绝对值小于这个阈值的系数被当作噪声，然后对其作置零处理；对绝对值大于阈值的小波系数用一定的方法进行缩减，然后得到缩减后的新值。

无噪信号的小波系数服从广义拉普拉斯分布，小波域内的斑点噪声部分服从高斯分布；选择b＝1，则式(4)变为拉普拉斯分布

经典的小波收缩方法有软阈值法和硬阈值法，但是在软阈值法中，较大的小波系数总是被阈值缩减，因此收缩后的信号的数学期望与收缩之前不同，所以处理后的图像相对平滑一些。硬阈值法的缺点是在零值域附近的小波系数被突然置零，导致了小波数据的不连续性，并且这使得信号的方差更大了，这些变换对于图像中的细节影响较大。但是在实际应用中，特别是噪声水平很高时，硬阈值法处理后的图像在不连续点周围会产生震荡，影响图像的去噪效果。

由于经典的阈值收缩方法不能满足对医学超声图像去噪的要求，所以本发明对收缩方法做了改进。

为了得到小波域内的信号估计值，使用贝叶斯最大后验估计的方法。在后验概率的计算过程中，使用贝叶斯公式如下

其中，P_(R|S)(r|s)表示小波域内的信号估计值，R表示无噪声图像小波量，S表示有噪声图像小波量。

将式(5)、式(8)带入上式(9)，得到

为了得到最大后验概率，将ln(p_R|S(r|s))对r求一次导数的方程置零，最后得到

为r的估计。这样就得到新的收缩方法

小波收缩函数在曲线图像上表现的更加平滑，尤其当小波系数大于小波阈值的区间范围内。

步骤4)利用三边滤波器对最后一层的低频部分(LL^J)中的小波系数做滤波处理

一般基于小波的去噪方法，即保留低频域(LL)的小波系数不变，仅对高频域(LH、HL、HH)的小波系数做阈值处理。然而，此方法应用于医学超声图像去噪时表现不佳。经过多次实验，发现低频域内的小波系数依然具有很多斑点噪声，传统的双边滤波器的优点是边缘保留，但通常会导致梯度失真，更重要是无法解决脉冲噪声。为了更有效地滤除低频域内的斑点噪声，本发明选择三边滤波器对低频域内的小波系数作滤波处理。我们提出了一个新的想法，在原始的双边滤波器的高斯距离权重和灰度权重基础上，增加了“脉冲”权重，也称为质量权重。此外，由于脉冲量在一定程度上表示了图像的梯度信息，三边滤波器解决了梯度失真的问题。

三边滤波器由传统的双边滤波器结构发展而来，结构如下

其中，f(ζ)表示原来被噪声污染的图像矩阵，x＝(x₁,x₂)是当前点的位置,ζ表示Ω域内的一个点，Ω_x(N):＝{x+(i,j):-N≤i,j≤N}表示一种连接关系；在实际效果中，选择Ω＝Ω_x(1)；w(x,ζ)加权函数为

w(x,ξ)＝w_S(x,ξ)w_R(x,ξ)(14)

w_S(x,ξ)表示区域滤波器，w_R(x,ξ)表示值域滤波器。

为实现三边滤波器，我们需要用加权函数计算出出图像中的噪声点。

首先引用函数f_m(x)来估计像素x是边缘点还是噪声点，d(x,ξ)表示x和ξ之间的像素差的绝对值，

d(x,ξ)＝|f(x)-f(ξ)|(15)

f_m(x)设为如下式：

g_i(x)为除d(x,ξ)外第i_th个最小值。

其想法为：如果一个像素点是为图像中边缘点，则其邻域中至少有一半左右的点和其灰度值差不多，从而有比较小的f_m(x)；否则，若其为被脉冲噪声污染的点，则其他点和这点灰度值差别较大，故有比较大的f_m(x)函数值。所以在增加脉冲量后，双边滤波器的加权函数改为

w'(x,ξ)＝w_S(x,ξ)w_R(x,ξ)^1-H(x,ξ)w_I(ξ)^H(x,ξ)(17)

表示脉冲权重函数，根据可知，当x为边缘点，H(x,ξ)≈0，当x为噪声点，H(x,ξ)≈1。参数σ_I表示补偿f_m(x)高值的近似阈值，参数σ_H表示控制函数H(x,ξ)的形状。

综上所述，经三边滤波器去噪后的图像h(x)被表示为

步骤5)作小波逆变换处理，得到去噪后的医学超声图像。

经过阈值收缩处理和三边滤波器处理就可以得到去噪后的小波系数，为了得到去噪后的超声图像，需要对小波系数进行小波逆变换，从而可以得到利于医师分析的去噪后的图像，通过实验也验证了本发明确实可以满足对于医学超声图像去噪的要求。

与现有技术相比，本发明的优点在于：

本发明一方面对通用的阈值函数和收缩方法进行了改进，对高频域的细小噪声具有很强的去噪能力；另一方面，由于本发明利用三边滤波器对低频部分进行了过滤处理，因此对于颗粒较大的斑点噪声(存在于低频部分)同样具有很强的抑制能力。由于本发明的时间的消耗很大部分是由于引入的用于处理低频部分噪声的滤波器，所以在本发明中三边滤波器的引入不仅能很好的提高去噪的性能并且大大提高了处理的效率。同时针对医学超声图像的特点，这种结合的方法不仅能很好的抑制斑点噪声，同时还能够保留图像中病灶边缘等的细节部分，能更好的帮助医师进行病情分析。

附图说明

图1是本发明的方法流程图

图2是本发明的实验过程示意图

图3是仿真图像，其中图3a是无噪图像，大小400×400，图3b是噪声图像，斑点噪声方差σ²＝0.1。

图4是仿真图像的去噪效果对比图，其中图4a是Lee算法的效果图，图4b是SRBF算法的效果图，图4c是Frost算法的效果图，图4d是三边滤波算法的效果图，图4e是DPAD算法的效果图，图4f是OBNLM算法的效果图，图4g是小波收缩算法的效果图，图4h是小波与双边的结合算法的效果图，图4i是本发明的效果图。

图5是肝脏临床超声图像。

图6是肝图像的去噪效果对比图，其中图6a是Lee算法的效果图，图6b是SRBF算法的效果图，图6c是Frost算法的效果图，图6d是三边滤波算法的效果图，图6e是DPAD算法的效果图，图6f是OBNLM算法的效果图，图6g是小波收缩算法的效果图，图6h是小波与双边的结合算法的效果图，图6i是本发明的效果图。

具体实施方式

为使本发明的目的、技术方案和优点更加清晰，下面结合附图就对本发明的技术方案作进一步描述。

基于小波变换和三边滤波器的医学超声图像去噪方法，包括如下步骤：

步骤1)建立医学超声图像模型

如果认为超声成像系统能够对那些影响声波功率的因素做出恰当的动态补偿，则超声成像系统采集的包络信号由两部分组成，一是有意义的体内组织的反射信号，另一部分是噪声信号。其中噪声信号可分为相乘噪声与相加噪声。相乘噪声与超声信号成像的原理有关，主要来源于随机的散射信号。相加噪声认为是系统噪声，如传感器的噪声等。对于包络线的通用模型模型如下

s(x,y)＝r(x,y)n(x,y)(1)

这里，(x,y)代表2D空间坐标，r(x,y)表示无噪声信号，n(x,y)表示相乘噪声。

为了适应超声成像系统显示屏幕的动态显示范围，对超声成像系统采集到的包络信号进行对数压缩处理。此时相乘的式(1)模型将变为相加的模型，如下

log(s(x,y))＝log(r(x,y))+log(n(x,y))(2)

此时，得到的信号log(s(x,y))即是通常看到的医学超声图像。

步骤2)对第一步得到的对数变换后的图像进行小波分解，得到四个频域(LL¹、LH¹、HL¹和HH¹)。

LL1分量是对原始信号LL0的列和行进行小波分解后得到的低频分量，即一级小波分解后近似部分，它包含了原始图像最多的低频信息；

LH1是一次小波分解后的垂直方向上的高频分量，即它包含了图像水平方向上的近似信息和垂直方向上的边缘等高频信息；

HL1是一次小波分解后的水平方向上的高频分量，即它包含了图像垂直方向上的近似信息和水平方向上的边缘等高频信息；

HH1是一次小波分解后对角方向上的高频分量，即它包含了图像水平和垂直方向上的边缘等高频信息；

对低频域LL¹继续进行小波分解，再得到第二层四个频域(LL²、LH²、HL²和HH²)。然后重复这个步骤，直到分解最大层数J。

由于小波变换是线性变换，因此式(3)模型经过二维离散小波变换后得到下面模型：

其中和分别表示含有噪声图像的小波系数、无噪声图像的小波系数和斑点噪声的小波系数。其中上标j为小波变换的分解层数，下标(l,k)为小波域内的坐标。

本发明认为经过小波分解后的无噪信号的小波系数服从广义拉普拉斯分布，其概率分布如下

式中，是伽马函数，b为形状参数，a为尺度参数，

同时斑点噪声的小波系数服从零均值高斯分布

式中σ_N为小波域内噪声的标准差。

步骤3)对每一层的高频部分(LH^j、HL^j和HH^j，j＝1,2,...,J)的小波系数进行阈值法收缩处理。

在小波去噪方法中，阈值函数的选择会直接影响到最终的图像去噪结果。当阈值选择较小时，一部分大于该阈值的噪声系数会被当作有用信号保留下来，这就导致去噪后的图像依然存在大量噪声；当阈值选择较大时，会将很多系数很小的有用信息当作噪声而置零，这将使得去噪后的图像变得很平滑，损失很多细节信息。因此选择恰当的小波阈值函数非常重要。

Donoho等人提出了一种典型的阈值选取方法，并且从理论上证明了该阈值与噪声的标准差成正比，改阈值函数又称为统一阈值函数，其公式如下

其中，M即是对应小波域内小波系数的总体个数，σ_n是噪声的标准差。在这种阈值函数中，阈值T受小波系数的个数影响较大，即当M过大时，较大的阈值可能会平滑掉那些系数较小的有用信息。

基于式(7)的基础之上，提出了一种更加适合超声图像的阈值函数，其公式如下

其中，T_j是新的阈值函数，σ_n是噪声的标准差，t_j代表j层的自适应参数。这是种常见的阈值改进的方法，t_j的选取是根据实验决定的，在小波分解后，在不同层分解的小波系数具有不同的分布，由此t_j的选择基于j层的选择，但这种选择不是最佳的，如果适当的选择，所提出的方法将反射更多的优越性。

在小波去噪方法中，首先选定一个给定阈值，然后按照一定的规则对小波系数进行收缩，便完成了对小波系数的去噪。即给定一个阈值，所有绝对值小于这个阈值的系数被当作噪声，然后对其作置零处理；对绝对值大于阈值的小波系数用一定的方法进行缩减，然后得到缩减后的新值。

无噪信号的小波系数服从广义拉普拉斯分布，小波域内的斑点噪声部分服从高斯分布；选择_b＝1，则式(4)变为拉普拉斯分布

经典的小波收缩方法有软阈值法和硬阈值法，但是在软阈值法中，较大的小波系数总是被阈值缩减，因此收缩后的信号的数学期望与收缩之前不同，所以处理后的图像相对平滑一些。硬阈值法的缺点是在零值域附近的小波系数被突然置零，导致了小波数据的不连续性，并且这使得信号的方差更大了，这些变换对于图像中的细节影响较大。但是在实际应用中，特别是噪声水平很高时，硬阈值法处理后的图像在不连续点周围会产生震荡，影响图像的去噪效果。

由于经典的阈值收缩方法不能满足对医学超声图像去噪的要求，所以本发明对收缩方法做了改进。

为了得到小波域内的信号估计值，使用贝叶斯最大后验估计的方法。在后验概率的计算过程中，使用贝叶斯公式如下

其中，P_(R|S)(r|s)表示小波域内的信号估计值，R表示无噪声图像小波量，S表示有噪声图像小波量。

将式(5)、式(8)带入上式(9)，得到

为了得到最大后验概率，将ln(p_R|S(r|s))对r求一次导数的方程置零，最后得到

为r的估计。这样就得到新的收缩方法

本发明改进的小波收缩函数在曲线图像上表现的更加平滑，尤其当小波系数大于小波阈值的区间范围内。

步骤4)利用三边滤波器对最后一层的低频部分(LL^J)中的小波系数做滤波处理

一般基于小波的去噪方法，即保留低频域(LL)的小波系数不变，仅对高频域(LH、HL、HH)的小波系数做阈值处理。然而，此方法应用于医学超声图像去噪时表现不佳。经过多次实验，发现低频域内的小波系数依然具有很多斑点噪声，传统的双边滤波器的优点是边缘保留，但通常会导致梯度失真，更重要是无法解决脉冲噪声。为了更有效地滤除低频域内的斑点噪声，本发明选择三边滤波器对低频域内的小波系数作滤波处理。我们提出了一个新的想法，在原始的双边滤波器的高斯距离权重和灰度权重基础上，增加了“脉冲”权重，也称为质量权重。此外，由于脉冲量在一定程度上表示了图像的梯度信息，三边滤波器解决了梯度失真的问题。

三边滤波器由传统的双边滤波器结构发展而来，结构如下

其中，f(ζ)表示原来被噪声污染的图像矩阵，x＝(x₁,x₂)是当前点的位置,ζ表示Ω域内的一个点，Ω_x(N):＝{x+(i,j):-N≤i,j≤N}表示一种连接关系；在实际效果中，选择Ω＝Ω_x(1)；w(x,ζ)加权函数为

w(x,ξ)＝w_S(x,ξ)w_R(x,ξ)(14)

w_S(x,ξ)表示区域滤波器，w_R(x,ξ)表示值域滤波器。

为实现三边滤波器，我们需要用加权函数计算出出图像中的噪声点。

首先引用函数f_m(x)来估计像素x是边缘点还是噪声点，d(x,ξ)表示x和ξ之间的像素差的绝对值，

d(x,ξ)＝|f(x)-f(ξ)|(15)

f_m(x)设为如下式：

g_i(x)为除d(x,ξ)外第i_th个最小值。

其想法为：如果一个像素点是为图像中边缘点，则其邻域中至少有一半左右的点和其灰度值差不多，从而有比较小的f_m(x)；否则，若其为被脉冲噪声污染的点，则其他点和这点灰度值差别较大，故有比较大的f_m(x)函数值。所以在增加脉冲量后，双边滤波器的加权函数改为

w'(x,ξ)＝w_S(x,ξ)w_R(x,ξ)^1-H(x,ξ)w_I(ξ)^H(x,ξ)(17)

表示脉冲权重函数，根据可知，当x为边缘点，H(x,ξ)≈0，当x为噪声点，H(x,ξ)≈1。参数σ_I表示补偿f_m(x)高值的近似阈值，参数σ_H表示控制函数H(x,ξ)的形状。

综上所述，经三边滤波器去噪后的图像h(x)被表示为

步骤5)作小波逆变换处理，得到去噪后的医学超声图像。

经过阈值收缩处理和三边滤波器处理就可以得到去噪后的小波系数，为了得到去噪后的超声图像，需要对小波系数进行小波逆变换，从而可以得到利于医师分析的去噪后的图像，通过实验也验证了本发明确实可以满足对于医学超声图像去噪的要求。

本方法的整体步骤的示意图如下图1所示。

实验验证

为了客观地评价本发明提出的去噪方法，以峰值信噪比(PSNR)、结构相似度(SSIM)、FoM(Pratt’sFigureofMerit)和运行时间Time作为图像质量评价标准。峰值信噪比的计算公式如下

式中，为信号X的估计值，MSE由下面公式计算得到

这里的M，N分别表示二维信号X的长度与宽度。

结构相似度能够量化两幅图像在结构上的差异，公式定义如下

式中，μ_X、和分别是参考图像和估计图像的均值和方差。是X和的协方差，c₁和c₂为常量。当c₁和c₂都选择为正数时，SSIM的取值范围为[01]，其中1为最好结果，表示两幅图的结构相同。

FoM能够客观地比较去噪图像的边缘检测质量，公式定义如下

式中，N_X和分别表示理想的和实际检测到的边缘像素个数。α为常数(通常取α＝1/9)，d_i表示为第i边缘像素点到最近理想边缘像素点的距离。FoM的取值范围为[01]，其中1为最好结果，表示为检测到的图像边缘和理想的图像边缘一致。这里检测边缘像素时使用的是Canny检测算法(高斯滤波器的标准差取值σ＝3)。

在本文中，为了让本发明更有说服力并且更好展现其优势，不仅将实验分成两个部分，一个是斑点噪声仿真实验，另一个是真实的临床医学超声图像(肝图像)；而且还做了与其他8种经典方法的对比实验，并结合上述四种量化指标来清晰地评价出本发明的优势，实验示意图如下图2所示：

为了能更定量且直观地验证本发明的优越性，先用仿真图(如图3)实验，并与8种经典算法对比，分别为Lee,SRBF,Frost,三边滤波,DPAD,OBNLM,小波收缩算法以及小波与双边的结合算法。实验结果见如下图4，对比指标量值见如下表1。

表1九大算法的去噪性能比较

算法	FOM	PSNR/dB	SSIM	Time/s
					噪声图像	0.3170	18.7605	0.4708	—
Lee	0.79695	22.20069	0.85591	3.51231
					SRBF	0.64657	24.94255	0.75256	11.08316
Frost	0.81464	23.54707	0.84193	11.37657
					三边滤波	0.82163	25.71782	0.82474	9.27613
DPAD	0.70546	25.37754	0.72347	32.76338
					OBNLM	0.87622	22.02718	0.70989	49.87029
小波收缩	0.33209	24.27656	0.71628	1.33192
					小波_双边	0.44828	25.56923	0.81025	9.44184
本发明方法	0.88251	26.31741	0.84792	9.51625

备注:红色下标指9种算法位于所在列的指标排序。

仿真图像是通过对无噪图像加入斑点噪声所得。通过以上实验得到的定量数据可看出，经典的Lee滤波速度快且对于原图像边缘的破坏性很小，但去噪效果最差。而SRBF虽然比Lee略胜一筹，但图像较模糊，而且为突出的优势指标。Frost除了边缘保持得不错外，对噪声的抑制能力较弱。DPAD在边缘保持和去噪上效果均不突出。OBNLM可以看到清晰地边缘，然而正相反，图像整体模糊不清。本发明方法的原分支——三边滤波和小波收缩，单独去噪的效果不佳。仅用前者，整体在抑制斑点噪声中的效果并不好；仅用后者，低频部分仍存在斑点噪声，所以需要引入滤波器对低频部分做滤波处理。如果引入的是双边滤波器，效果虽然优于原生小波，但梯度失真较明显，且性能不高。但通过比较发现，本发明在低频引用三边滤波，与小波变换相结合，不仅边缘保持得很好，而且解决了梯度失真问题，极大地提高了去噪后的图像质量。

再利用临床超声图像进行试验验证，选择的是如图5的肝脏临床超声图像。实验结果如图6。

由于现实中并不存在无噪图像，因此以上所列的PSNR等四大质量指标在此无法有效地使用，通过调研和文献检索，本文引入另一个指标对比，即无参考图像质量指标(NIQE)。评价指标对比结果见表2。

表2九大算法的NIQE指标对比

由以上定量数据可直观看出，本发明方法在实际应用于医学超声图像的过程和结果与应用与仿真图像时如出一辙，不仅去噪效果得到了显著提高，而且更好地保留图像边缘信息，并无梯度失真，从而达到医学超声图像对于去噪的要求。

Claims (1)

1.基于小波变换和三边滤波器的医学超声图像去噪方法，包括如下步骤：

步骤1)建立医学超声图像模型；

如果认为超声成像系统能够对那些影响声波功率的因素做出恰当的动态补偿，则超声成像系统采集的包络信号由两部分组成，一是有意义的体内组织的反射信号，另一部分是噪声信号；其中噪声信号可分为相乘噪声与相加噪声；相乘噪声与超声信号成像的原理有关，主要来源于随机的散射信号；相加噪声认为是系统噪声，如传感器的噪声等；对于包络线的通用模型模型如下：

s(x,y)＝r(x,y)n(x,y)(1)

这里，(x,y)代表2D空间坐标，r(x,y)表示无噪声信号，n(x,y)表示相乘噪声；

为了适应超声成像系统显示屏幕的动态显示范围，对超声成像系统采集到的包络信号进行对数压缩处理；此时相乘的式(1)模型将变为相加的模型，如下

log(s(x,y))＝log(r(x,y))+log(n(x,y))(2)

此时，得到的信号log(s(x,y))即是通常看到的医学超声图像；

步骤2)对第一步得到的对数变换后的图像进行小波分解，得到四个频域(LL¹、LH¹、HL¹和HH¹)；

LL1分量是对原始信号LL0的列和行进行小波分解后得到的低频分量，即一级小波分解后近似部分，它包含了原始图像最多的低频信息；

LH1是一次小波分解后的垂直方向上的高频分量，即它包含了图像水平方向上的近似信息和垂直方向上的边缘等高频信息；

HL1是一次小波分解后的水平方向上的高频分量，即它包含了图像垂直方向上的近似信息和水平方向上的边缘等高频信息；

HH1是一次小波分解后对角方向上的高频分量，即它包含了图像水平和垂直方向上的边缘等高频信息；

对低频域LL¹继续进行小波分解，再得到第二层四个频域(LL²、LH²、HL²和HH²)；然后重复这个步骤，直到分解最大层数J；

由于小波变换是线性变换，因此式(3)模型经过二维离散小波变换后得到下面模型：

$S_{l,k}^j=R_{l,k}^j+N_{l,k}^j,j=1,2,...,J.(l,k)\∈Z^2---\left(3\right)$

其中和分别表示含有噪声图像的小波系数、无噪声图像的小波系数和斑点噪声的小波系数；其中上标j为小波变换的分解层数，下标(l,k)为小波域内的坐标；

经过小波分解后的无噪信号的小波系数服从广义拉普拉斯分布，其概率分布如下

$p_R\left(r\right)=\frac{b}{2a\mathrm{\Γ}\left(\frac{1}{b}\right)}\exp (-|\frac{r}{a}{|}^b),a,b>0---\left(4\right)$

式中，是伽马函数，b为形状参数，a为尺度参数，

同时斑点噪声的小波系数服从零均值高斯分布

$p_N\left(n\right)=\frac{1}{\sqrt{2\mathrm{\π}}{\mathrm{\σ}}_N}\exp (-\frac{n^2}{2{\mathrm{\σ}}_N^2})---\left(5\right)$

式中σ_N为小波域内噪声的标准差；

步骤3)对每一层的高频部分(LH^j、HL^j和HH^j，j＝1,2,...,J)的小波系数进行阈值法收缩处理；

基于式(7)的改阈值函数，提出了一种更加适合超声图像的阈值函数，其公式如式(8)：

$T={\mathrm{\σ}}_n\sqrt{2\log M}---\left(6\right)$

其中，M即是对应小波域内小波系数的总体个数，σ_n是噪声的标准差；

$T_j=t_j{\mathrm{\σ}}_N\sqrt{2\log M}---\left(7\right)$

其中，T_j是新的阈值函数，σ_n是噪声的标准差，t_j代表j层的自适应参数；这是种常见的阈值改进的方法，t_j的选取是根据实验决定的，在小波分解后，在不同层分解的小波系数具有不同的分布，由此t_j的选择基于j层的选择；

在小波去噪方法中，首先选定一个给定阈值，然后按照一定的规则对小波系数进行收缩，便完成了对小波系数的去噪；即给定一个阈值，所有绝对值小于这个阈值的系数被当作噪声，然后对其作置零处理；对绝对值大于阈值的小波系数用一定的方法进行缩减，然后得到缩减后的新值；

无噪信号的小波系数服从广义拉普拉斯分布，小波域内的斑点噪声部分服从高斯分布；选择b＝1，则式(4)变为拉普拉斯分布

$p_R\left(r\right)=\frac{1}{2a}\exp (-|\frac{r}{a}\left|\right),a>0---\left(8\right)$

为了得到小波域内的信号估计值，使用贝叶斯最大后验估计的方法；在后验概率的计算过程中，使用贝叶斯公式如下

$p_{R|S}\left(r\right|s)=\frac{1}{p_S\left(s\right)}{p}_N(s-r)p_R\left(r\right)---\left(9\right)$

其中，P_(R|S)(r|s)表示小波域内的信号估计值，R表示无噪声图像小波量，S表示有噪声图像小波量。

将式(5)、式(8)带入上式(9)，得到

$p_{R|S}\left(r|s\right)=\frac{1}{2\sqrt{2\mathrm{\π}}{\mathrm{a\σ}}_N{p}_S\left(s\right)}\exp \left(-\frac{2{\mathrm{\σ}}_N^2\left|r\right|+a{\left(s-r\right)}^2}{2{\mathrm{a\σ}}_N^2}\right)---\left(10\right)$

为了得到最大后验概率，将ln(p_R|S(r|s))对r求一次导数的方程置零，最后得到

$\hat{r}=\mathrm{sgn}\left(s\right)max\left(\right|s|-\frac{{\mathrm{\σ}}_N^2}{a},0)---\left(11\right)$

为r的估计；这样就得到新的收缩方法

$\hat{r}=\left\{\begin{array}{cc}0& s\≤T_j\\ \mathrm{sgn}\left(s\right)max\left(\right|s|-\frac{{\mathrm{\σ}}_N^2}{a},0)& s>T_j\end{array}\right.---\left(12\right)$

小波收缩函数在曲线图像上表现的更加平滑，尤其当小波系数大于小波阈值的区间范围内；

步骤4)利用三边滤波器对最后一层的低频部分(LL^J)中的小波系数做滤波处理；

三边滤波器的结构如下

$h\left(x\right)=\frac{\underset{\mathrm{\ξ}\∈\mathrm{\Ω}}{\Σ}w(x,\mathrm{\ξ})f\left(\mathrm{\ξ}\right)}{\underset{\mathrm{\ξ}\∈\mathrm{\Ω}}{\Σ}w(x,\mathrm{\ξ})}---\left(13\right)$

其中，f(ζ)表示原来被噪声污染的图像矩阵，x＝(x₁,x₂)是当前点的位置,ζ表示Ω域内的一个点，Ω_x(N):＝{x+(i,j):-N≤i,j≤N}表示一种连接关系；在实际效果中，选择Ω＝Ω_x(1)；w(x,ζ)加权函数为

w(x,ξ)＝w_S(x,ξ)w_R(x,ξ)(14)

w_S(x,ξ)表示区域滤波器，w_R(x,ξ)表示值域滤波器；

为实现三边滤波器，需要用加权函数计算出出图像中的噪声点；

首先引用函数f_m(x)来估计像素x是边缘点还是噪声点，d(x,ξ)表示x和ξ之间的像素差的绝对值，

d(x,ξ)＝|f(x)-f(ξ)|(15)

f_m(x)设为如下式：

$f_m\left(x\right)=\underset{i=1}{\overset{m}{\Σ}}{g}_i\left(x\right)---\left(16\right)$

g_i(x)为除d(x,ξ)外第i_th个最小值；

如果一个像素点是为图像中边缘点，则其邻域中至少有一半左右的点和其灰度值差不多，从而有比较小的f_m(x)；否则，若其为被脉冲噪声污染的点，则其他点和这点灰度值差别较大，故有比较大的f_m(x)函数值；所以在增加脉冲量后，双边滤波器的加权函数改为

w'(x,ξ)＝w_S(x,ξ)w_R(x,ξ)^1-H(x,ξ)w_I(ξ)^H(x,ξ)(17)

表示脉冲权重函数，根据可知，当x为边缘点，H(x,ξ)≈0，当x为噪声点，H(x,ξ)≈1；参数σ_I表示补偿f_m(x)高值的近似阈值，参数σ_H表示控制函数H(x,ξ)的形状；

综上所述，经三边滤波器去噪后的图像h(x)被表示为

$h\left(x\right)=\frac{\underset{\mathrm{\ξ}\∈\mathrm{\Ω}}{\Σ}{w}^{\′}(x,\mathrm{\ξ})f\left(\mathrm{\ξ}\right)}{\underset{\mathrm{\ξ}\∈\mathrm{\Ω}}{\Σ}{w}^{\′}(x,\mathrm{\ξ})}---\left(18\right)$

步骤5)作小波逆变换处理，得到去噪后的医学超声图像；

经过阈值收缩处理和三边滤波器处理就可以得到去噪后的小波系数，为了得到去噪后的超声图像，需要对小波系数进行小波逆变换，从而可以得到利于医师分析的去噪后的图像。