基于方向加权中值滤波的超声图像处理方法

摘要

本发明公开了 基于方向加权中值滤波的超声图像处理方法，该方法在中值滤波的基础上，考虑了图像的边缘方向和边缘幅度信息，针对某个具体的当前像素点，只选取邻域窗口中与当前点边缘方向相近的点做中值滤波，而且中值的排序选取是在用边缘幅度加权后的像素值上进行的，本发明的方法能够很好的区分超声图像的噪声和边缘信息，使得在去噪的同时有效地保护了超声图像的边缘细节信息，更加有利于后期的诊断与分析。而且本发明采用的基础算法是中值滤波，算法简单，时间复杂度低，处理耗时少，更加符合超声成像系统实时处理显示图像的要求。

说明书

技术领域

本发明涉及超声成像中的数据后处理技术，特别涉及提高超声图像质量的图 像处理技术，更具体地说，涉及一种用于超声图像去噪和增强的方向加权中值 滤波方法。

背景技术

在超声图像形成过程中，当超声波长与照射物体表面粗糙度相当时，就会产 生斑点噪声，这一现象可以用随机散射模型来解释。这些噪声的存在，使得超 声图像的清晰度不高，这也是超声成像的主要缺点之一。超声图像特有的散斑 噪声不仅使得超声图像的质量较差，尤其是掩盖和降低了图像某些细节信息， 也使得对图像细节的识别和分析变得更加困难。为以后的病情诊断及定量分析， 图像特征提取和识别造成不利的影响。因此，抑制这些散斑噪声，改善图像质 量是超声图像分析和识别的重要预处理环节。

超声图像的斑点噪声去除问题一直是国内外超声成像技术的重要课题之 一。超声图像去噪一般要求有效抑制散斑噪声，同时要保留对后期分析和诊断 有用的图像细节信息。超声图像去噪的主要难点在于：1)散斑噪声可以大致看 作为一种乘性噪声；2)噪声的随机性质比较复杂；3)噪声易与图像细节相混， 而图像细节又复杂且多样。

目前，针对抑制超声图像散斑噪声，改善图像质量，人们已经提出了许多 方法。1)图像平均方法利用在不同时间、不同频率或不同扫描方向得到同一目 标的一系列图像，将它们平均形成一幅复合图像以提高图像的信噪比(SNR)。尽 管这种方法简单、快速，然而它受到一些限制：需要严格控制系列图像的形成， 而且需要图像的配准；由于图像模糊的影响，一些小的细节(例如小的血管、 纹理等)会丢失，因此降低了空间分辨能力。2)自适应加权中值滤波方法对图 像中的每一个像素点的值，用其局部邻域窗口的加权中值代替。设定窗口中心 点为(i_n，i_n)，窗口大小为(2w+1)×(2w+1)，则窗口内各点加权系数计算为： weight(i，j)＝[w(i_o，i_o)-a*d*σ²/m]

其中，[·]为取整运算，a为系数，d为点(i，j)到中心点(i_n，i_n)的距离，σ2 和m分别为局部邻域窗口内的方差和均值。取邻域窗口内的加权中值代替中心 点的像素值。其中，邻域窗口内各点的权值能够根据图像的局部统计特征来自 动选取，邻域窗口的大小也能够根据其局部信噪比自动调节。相对于单纯的中 值滤波，这种自适应加权处理在保留细节方面取得了一定的效果。但是对窗口 的选择很敏感，限制了处理效果，在去除噪声的同时，也会导致一些细微细节 的损失。3)小波阈值收缩是一类重要的去除散斑噪声的方法，这种方法主要基 于Donoho提出的小波软阈值(soft-thresholding)收缩去噪：首先将图像小波 分解，设定一个阈值，对于小于阈值的小波系数用零代替，而对于大于阈值的 小波系数用它减去阈值代替，被处理后的小波系数作小波反变换，就可以得到 重建图像。同时，还有许多对于小波阈值的改进算法。然而，上述方法主要是 针对高斯分布的噪声，解决其它分布的噪声时，效果不甚理想；同时，小波阈 值的选取也是一个需要着重解决的问题，既不能选的太大而丢失一些边缘和局 部细节，也不能选的太小而对噪声的抑制不充分。如何选择小波变换的尺度和 阈值尚无确定方法。4)基于各向异性扩散模型的去噪方法。各向异性扩散实际 上是一个非线性的偏微分方程，由图像的梯度来决定扩散速度，能够同时兼顾 噪声消除和特征保持两个方面。目前，以Perona-Malik(PM)模型为代表的各向 异性扩散方程(Anisotropic Diffusion，AD)得到了广泛的应用。Yu等将各向 异性扩散方程应用到散斑抑制中，提出了去除散斑噪声的各向异性扩散模型 (Speckle Reducing Anisotropic Diffusion，SRAD)。Yu等修正了扩散系数， 使扩散方程能够根据图像噪声的情况而调整扩散系数，并且能够对图像的细节 信息更加敏感。但SRAD模型也存在明显的缺陷：模型中尺度函数是由初始图像 中尽可能大的均匀区域计算得到的，模型的关键是如何选取图像中一个尽可能 大的，合适的均匀区域，这一区域的选取往往很大程度上影响着扩散结果，给 实验结果带来较大的偶然性。此外，SRAD模型所使用的局部统计信息实际上是 各向同性的，这也背离了各向异性扩散算法的本质。而且各向异性扩散需要多 次迭代处理才能取得比较好的效果，迭代又是非常耗时的，所以该算法很难满 足超声实时处理的需求。

总结一下，现有的去除超声图像散斑噪声的方法有两个比较关键的不足需 要提高：一是在对图像平滑去除噪声的同时，对细节和边缘的模糊比较严重， 导致了一些边缘细节的损失，影响了后期的诊断与分析；二是处理的时间复杂 度比较高，难以满足超声成像系统实时处理显示图像的需求。

发明内容

本发明的目的是提供基于方向加权中值滤波的超声图像处理方法，该处理 方法在去噪的同时有效地保护了超声图像的边缘细节信息，而且本处理方法采 用的基础算法是中值滤波，算法简单，时间复杂度低。

为达到上述目的，本发明采用了以下技术方案：基于方向加权中值滤波的超 声图像处理方法，包括以下步骤：

(1)计算所述超声图像中每个像素点的边缘幅度和边缘方向；

(2)取所述图像中的一像素点作为当前点；

(3)根据所述当前点的边缘幅度或边缘方向判断该当前点是否为边缘点， 如果不是边缘点，则对该当前点进行自适应加权中值滤波处理；如果是边缘点， 则对该当前点进行下面步骤(4)(5)(6)的处理；

(4)取该当前点的M*N邻域，选择该M*N邻域中与该当前点的边缘方向相 关的点作为处理点；

(5)取所述处理点的像素值和边缘幅度值，对所述处理点的像素值用相应 的边缘幅度值加权，对加权后的值进行排序，找出中间值；

(6)用所述中间值对应的像素值替换该当前点的像素值；

(7)取下一个像素点作为当前点进行步骤(3)的处理，直到处理完所述图 像中的所有像素点，最后输出处理后的图像。

步骤(1)的一种具体实施方式，所述像素点的边缘幅度和边缘方向的计算 方法包括：

(1a)计算该像素点的水平方向和垂直方向的梯度，设定该像素点的位置为 (i，j)，该像素点的像素值为I(i，j)，则该像素点的水平方向的梯度Gx(i， j)、垂直方向的梯度Gy(i，j)的计算方式如下：

Gx(i，j)＝I(i，j+1)-I(i，j-1)    Gy(i，j)＝I(i+1，j)-I(i-1，j)；

(1b)计算该像素点的梯度幅度GradientAm和梯度方向GradientAn，计算 方式如下：

GradientAn＝atan(Gy/Gx)

(1c)判断该像素点的梯度幅度GradientAm是否大于预先设定的阈值，如 果该梯度幅度GradientAm大于或等于所述阈值，则取该梯度幅度GradientAm、 该梯度方向GradientAn为该像素点的边缘幅度和边缘方向，表示该像素点为边 缘点；如果该梯度幅度GradientAm小于所述阈值，则设定该像素点的边缘幅度 和边缘方向均为负值，表示该像素点不是边缘点。

步骤(1)的另一种具体实施方式，所述像素点的边缘幅度和边缘方向的计 算方法包括：

(1a)选取不同方向的边缘检测算子；

(1b)将每个方向的边缘检测算子中的数值与该像素点的邻域中的对应像素 值相乘后求和，求和后取绝对值作为输出值；

(1c)在不同方向的输出值中，将最大值减去最小值得差值，如果该差值大 于或等于预先设定的阈值，则取所述不同方向的输出值中的最大值为该像素点 的边缘幅度值，取该最大值对应的边缘检测算子的方向为该像素点的边缘方向， 表示该像素点为边缘点；如果该差值小于所述阈值，则设定该像素点的边缘幅 度和边缘方向均为负值，表示该像素点不是边缘点。

步骤(4)的第一种具体实施方式，从所述M*N邻域中选择与所述当前点的 边缘方向相关的点作为处理点的方法包括：

(4a)选取不同方向的取点模板；

(4b)从所述不同方向的取点模板中，取与所述当前点的边缘方向最接近的 方向对应的取点模板；

(4c)在所述M*N邻域中，取被所述最接近的方向对应的取点模板击中的点 作为处理点。

步骤(4)的第二种具体实施方式，从所述M*N邻域中，选择边缘方向与所 述当前点的边缘方向的差值小于预先设定的阈值的点作为处理点。

步骤(5)具体实施时，所述处理点的像素值用相应的边缘幅度值进行加权 处理指得是该像素值与相应的边缘幅度值进行相乘处理。

由于上述技术方案的运用，本发明与现有技术相比具有下列优点：本发明 的基于方向加权中值滤波的超声图像处理方法，在中值滤波的基础上，考虑了 图像的边缘方向和边缘幅度信息，针对某个具体的当前像素点，只选取邻域窗 口中与当前点边缘方向相近的点做中值滤波，而且中值的排序选取是在用边缘 幅度加权后的像素值上进行的，本发明的方法能够很好的区分超声图像的噪声 和边缘信息，使得在去噪的同时有效地保护了超声图像的边缘细节信息，更加 有利于后期的诊断与分析。而且本发明采用的基础算法是中值滤波，算法简单， 时间复杂度低，处理耗时少，更加符合超声成像系统实时处理显示图像的要求。

附图说明

附图1为本发明的流程示意图；

附图2为本发明中通过梯度计算图像像素点的边缘幅度和边缘方向的 流程示意图；

附图3为本发明中通过边缘检测算子计算图像像素点的边缘幅度和边 缘方向的流程示意图；

附图4显示了边缘检测算子作用的2*2邻域的取法示意图；

附图5显示了四个不同方向的边缘检测算子；

附图6为本发明中通过当前点的边缘方向选择处理点的流程示意图；

附图7显示了四个不同方向的取点模板；

附图8为本发明中通过当前点与邻域中的点的边缘方向差来选择处理 点的流程示意图。

具体实施方式

下面结合附图来进一步阐述本发明。

本发明提出了一种方向加权中值滤波算法，该算法的基本思路如图1所示， 以超声图像为背景，在中值滤波的基础上，考虑了图像的边缘方向和边缘幅度 信息，针对某个具体的当前点，只选取滤波窗口中与当前点边缘方向相近的点 做中值滤波，而且中值的排序选取是在用边缘幅度加权后的像素值上进行的。 方向加权中值滤波算法考虑到边缘不同于背景，其灰度分布存在明显的差异， 而且这种差异具有方向性，因此采用边缘幅度和边缘方向结合的方法。

基于方向加权中值滤波算法，用于超声图像的去噪和增强。该算法的流程 图参见图1，首先计算图像I中每个像素点的边缘幅度和边缘方向；然后在此基 础上，逐点处理，针对每一像素点，根据边缘方向选择其M*N邻域中方向相关 的点作为处理点；最后，针对邻域中选出的处理点，用边缘幅度进行加权中值 滤波处理。在导入超声图像数据后，主要包括以下步骤：

(1)计算超声图像中每个像素点的边缘幅度EdgeAm和边缘方向EdgeAn；

(2)取图像中的一像素点作为当前点；

(3)根据当前点的边缘幅度或边缘方向判断该当前点是否为边缘点，如果 不是边缘点，则对该当前点进行自适应加权中值滤波处理，该自适应加权中值 滤波处理的方法采用的是现有技术，具体参见“背景技术”；如果是边缘点，则 对该当前点进行下面步骤(4)(5)(6)的处理；

(4)取该当前点的M*N邻域，选择该M*N邻域中与该当前点的边缘方向相 关的点作为处理点；

(5)取上述处理点的像素值和边缘幅度值，对处理点的像素值用相应的边 缘幅度值加权，对加权后的值进行排序，找出中间值；

(6)用上述中间值对应的像素值替换当前点的像素值；

(7)取下一个像素点作为当前点进行步骤(3)的处理，直到处理完图像中 的所有像素点，最后输出处理后的图像。

该算法中，通过边缘幅度控制减少了噪声对中值滤波的影响，边缘方向信 息的加入又使得该算法对边缘方向特别敏感，保留细节能力强，减少了边缘模 糊的程度。因此，这种改进的中值滤波算法，在去除散斑噪声能力、保持边缘 能力以及处理速度等指标上均取得了比较好的效果。

上述算法中，步骤(1)计算图像I中每点的边缘幅度和方向，有很多种方 法能够实现。这里给出两种具体的实施方法：1.通过图像的梯度幅度和方向得 到图像的边缘幅度和方向；2.用不同方向的边缘检测算子得到图像的边缘幅度 和方向。但不仅限于这两种方法。

方法1：通过图像的梯度幅度和方向得到图像的边缘幅度和方向的方法，其 流程图参见图2。首先，根据水平方向的梯度Gx(i，j)和垂直方向的梯度Gy (i，j)计算图像的梯度幅度GradientAm和梯度方向GraditenAn，计算方式如 下：

设定图像像素点的位置为(i，j)，该像素点的像素值为I(i，j)，则该像 素点的Gx(i，j)、Gy(i，j)：

Gx(i，j)＝I(i，j+1)-I(i，j-1)    Gy(i，j)＝I(i+1，j)-I(i-1，j)

则该像素点的梯度幅度GradientAm和梯度方向GraditenAn：

GradientAn＝atan(Gy/Gx)

其次，判断该像素点的梯度幅度值是否大于预先设定的阈值GradientAmThresh。 如果大于或等于阈值，则表示该像素点为边缘点，把梯度幅度看作是边缘幅度， 梯度方向看作是边缘方向；如果小于阈值，则认为无边缘(表示该像素点不是 边缘点)，把边缘幅度和方向都设成负值(比如-1)，以便在后面处理中区分像 素点是否为边缘点。即：

$\left(\begin{array}{cc}\left(\begin{array}{c}\mathrm{EdgeAm}(i,j)=\mathrm{GradientAm}(i,j)\\ \mathrm{EdgeAn}(i,j)=\mathrm{GranientAn}(i,j)\end{array}\right)& \mathrm{GradientAm}(i,j)\ge \mathrm{GradientAmThresh}\\ \left(\begin{array}{c}\mathrm{EdgeAm}(i,j)=-1\\ \mathrm{EdgeAn}(i,j)=-1\end{array}\right)& \mathrm{GradientAm}(i,j)<\mathrm{GradientAmThresh}\end{array}\right)$

对图像I的所有像素点逐点进行上述处理，则得到每个像素点的边缘幅度 EdgeAm和边缘方向EdgeAn。

方法2：通过不同方向的边缘检测算子得到图像边缘幅度和方向的方法， 参见图3，首先选取一组有不同方向的边缘检测算子；其次把这组算子分别作用 于图像上的一像素点的邻域上，得到每个方向上算子的输出；比较各个方向上 的输出，如果差别比较大，选取最大的输出作为当前像素点的边缘幅度，最大 输出对应的算子的方向为边缘方向；如果差别不大，则认为当前像素点不是边 缘点。

举例说明，选取如图5所示的四个方向(90度，0度，45度，135度)上 的边缘检测算子，使用该组边缘检测算子计算图像上每点的边缘幅度和方向， 参见图3，首先，针对图像中的一像素点(i，j)，取如图4所示的2*2邻域上的 点；然后，把图5所示的四个方向的边缘检测算子分别作用于图4的2*2邻域， 取每个算子模板与邻域模板中的对应数值相乘后求和的绝对值，作为输出值 0peratorAm＝[am90，am0，am45，am135]，相应的方向0peratorAn＝[90，0，45，135]； 接着，比较四个方向上的输出值，如果差别比较大，也就是四个方向的输出值 中最大值减去最小值的差大于或等于预先设定的阈值0peratorAmThresh，则表 示该像素点为边缘点，当前像素点(i，j)的边缘幅度为输出的最大值，边缘方向 为最大值对应的边缘检测算子的方向，假设输出的最大值为am45，则边缘幅度 为am45，边缘方向为45。如果输出差别很小，也就是最大值减去最小值的差小 于阈值，则认为无边缘(表示当前像素点不是边缘点)，当前像素点的边缘幅度 和方向都设成负值(比如-1)。用公式表示如下：

$\left(\begin{array}{c}\left(\begin{array}{cc}\mathrm{EdgeAm}(i,j)=\max \left(\mathrm{OperatorAm}\right)& \max \left(\mathrm{OperatorAm}\right)-\min \left(\mathrm{OperatorAm}\right)\\ \mathrm{EngeAn}(i,j)=\mathrm{OperatorAn}\left(\mathrm{index}\right)& \ge \mathrm{OperatorAmThresh}\end{array}\right)\\ \left(\begin{array}{cc}\mathrm{EdgeAm}(i,j)=-1& \max \left(\mathrm{OperatorAm}\right)-\min \left(\mathrm{OperatorAm}\right)\\ \mathrm{EdgeAn}(i,j)=-1& <\mathrm{OperatorAmThresh}\end{array}\right)\end{array}\right)$

其中，index为max(OperatorAm)对应的序号。

对图像的所有像素点都按照上述步骤处理，就能得到每个像素点的边缘幅度 EdgeAm和方向EdgeAn。

上述基于方向加权中值滤波算法的步骤(4)中，如何选择M*N邻域中与当 前点(i，j)的边缘方向相关的点作为处理点，这里也给出了两种具体实现方法： 1.以当前点的边缘方向为准，选择邻域中这个方向上的点为处理点；2.选择与 当前点的边缘方向差小于预设阈值的点为处理点。但不仅限于这两种方法。

方法1：以当前点的边缘方向为准，选择邻域中这个方向上的点为处理点， 方法1的流程图如图6所示。首先选取不同方向的取点模板，如图7，图7显示 了四个方向(90度，0度，45度，135度)的取点模板，判断当前点(i，j)的边 缘方向EdgeAn(i，j)与哪个方向(90度，0度，45度，135度)最接近，也就是 说EdgeAn(i，j)落在如下哪个角度区间里面：[67.5，112.5]， [0，22.5]U[157.5，180]，[22.5，67.5]，[112.5，157.5]；然后，按照如下规 则选取邻域中沿着这个方向上的点：

如果与90度最接近，也就是落在角度区间[67.5，112.5]中，选择图7中 模板2a)或者2b)击中的点作为处理点；

如果与0度最接近，也就是落在区间[0，22.5]U[157.5，180]中，选择图7 中模板3a)或者3b)击中的点作为处理点；

如果与45度最接近，也就是落在区间[22.5，67.5]中，选择图7中模板4a) 或者4b)击中的点作为处理点；

如果与135度最接近，也就是落在区间[112.5，157.5]中，选择图7中模 板5a)或者5b)击中的点作为处理点。

假设选出的处理点的个数为k，边缘幅值为

LocalEdgeAm＝[ea1，ea2，ea3，…，eak]，像素值为

LocalPixel＝[p1，p2，p3，…，pk]。

需要说明的是，为了直观起见，图7中只给出了5*5大小的取点模板2a)-5a) 或者2b)-5b)，实际上这里的模板大小应该是M*N大小的。其中，2a)-5a)是针 对边缘宽度是单点的情况下的取点模板，2b)-5b)是针对边缘宽度比较宽的情况 下的取点模板。如果选择模板2a)-5a)，容易出现误把噪声点增强，如果选择模 板2b)-5b)，容易把一些细小的边缘抹掉，所以要根据具体脏器的边缘宽窄情况， 选择合适的取点模板。

方法2：选择与当前点(i，j)的边缘方向的差小于预设阈值的点为处理点， 方法2的流程图如图8所示。在当前点(i，j)的邻域中的N*M个点中，选择边缘 方向EdgeAn与当前点(i，j)边缘方向EdgeAn(i，j)的差小于设定的阈值 EdgeAnThresh的点为处理点。假设选中的处理点有k个，它们的边缘幅度值为 LocalEdgeAm＝[ea1，ea2，ea3，…，eak]，像素值为 LocalPixel＝[p1，p2，p3，…，pk]。

上述基于方向加权中值滤波算法中，步骤(5)(6)中涉及到的，在选出 了处理点并拿到了它们的边缘幅度值和像素值后，进行加权中值滤波处理，首 先，用边缘幅度对像素值加权，即：LEAP＝[ea1*p1，ea2*p2，ea3*p3，…，eak*pk]； 其次，把加权后的值LEAP排序，取中间值median(LEPA)所对应的像素值px代 替当前点(i，j)处的像素值，如果中间值为ea2*p2，则用p2代替点(i，j)的像素 值。

本发明的基于方向加权中值滤波的超声图像处理方法，能够很好的区分超 声图像的噪声和边缘信息，使得在去噪的同时有效地保护了超声图像的边缘细 节信息，更加有利于后期的诊断与分析。而且本发明采用的基础算法是中值滤 波，算法简单，时间复杂度低，处理耗时少，更加符合超声成像系统实时处理 显示图像的要求。