CN100538740C - 医学超声图像各向异性扩散降噪方法 - Google Patents

Abstract

本发明提供了一种医学超声图像各向异性扩散降噪方法。首先计算每个像素的边缘信赖度，并由此确定每个像素的扩散参数取值；然后进行各向异性扩散算法降噪，即根据像素的梯度和扩散参数确定扩散程度，同时更新像素的灰度值；最后计算此次图像扩散降噪结果与扩散前的图像的非线性相关系数，如果该系数低于指定值则图像降噪处理结束，否则回到上一步，即根据像素新的梯度和扩散参数确定扩散程度，并更新像素的灰度值，直至相邻两次图像的非线性相关系数低于指定值。本发明显著提高了扩散算法进行图像降噪的适应性，有着广泛的应用前景。

Description

医学超声图像各向异性扩散降噪方法

(一)技术领域

本发明属于图像处理技术领域，具体涉及一种医学超声图像处理技术。

(二)背景技术

在医学超声成像中，由于相干波叠加造成的散斑噪声对超声图像影响严重，降低了超声成像系统对目标的分辨能力，影响了人们对图像的理解和自动特征的检测与提取，限制了医学超声图像在临床诊断中的进一步应用。对超声图像进行降噪处理，是医学超声图像处理的重要一环。超声图像降噪算法要在有效抑制散斑噪声的同时，保留对后期的分析和诊断有用的图像细节信息。

图像降噪技术一般分为空间域和频域处理两类。由于医学超声图像处理往往要求实时性比较好，而频率域的傅立叶变换数据量大且运算时间长，因此通常采用空间域降噪方法。常用的空间域线性降噪技术是邻域平均法，即将邻域各点的某种加权平均作为该点的滤波结果，其实质就是根据所取邻域内的信息确定一个新值来取代原值，以使它符合某种预定义的规则。为了避免引起亮度的偏移，该方法在抑制噪声的同时会引起图像的模糊和重要细节的丢失，不适于医学超声图像处理。还有一种多图像平均方法，若噪声是互不相关的零均值加性噪声，可用该方法来达到降噪的目的，该方法在医学超声图像应用中的最大困难是难以取得多幅精确定位的人体剖面图像。空间域线性降噪算法理论基础完善、数学处理简单且易于实现，但是会模糊图像的细节，因此它不是图像滤波的主流，尤其是对于医学超声图像更是如此。

由于线性算法在医学超声图像处理上的缺点，非线性医学超声图像降噪算法得到了充分的研究和发展，目前医学超声图像的降噪方法主要是中值滤波、小波变换和各向异性扩散算法等。基于中值滤波的方法根据图像的局部统计特征来自动选取滤波窗口内各点权值及窗口的大小和形状，尽管在保留图像细节方面取得了一定效果，但这类基于窗口的平滑方法对窗口的选择很敏感，限制了处理效果。基于小波变换的方法将超声图像变换到小波域，利用小波阈值处理将某些系数丢弃，再逆变换以去除噪声，但在该方法中很难选择合适的尺度。

基于各向异性扩散方程的医学超声图像降噪技术用选择性扩散方式代替了高斯平滑这种各向同性的扩散方式，即在图像的噪声处有较大的扩散系数，对图像平滑较多，而在图像的边界处扩散系数较小，使图像细节得到保持。各向异性扩散算法为

$\left\{\begin{array}{c}{I}_{x,y}^{t+1}=I_{x,y}^t+\mathrm{\λ}{[C_N{\▿}_{N}I+C_S{\▿}_{S}I+C_E{\▿}_{E}I+C_W{\▿}_{W}I]}_{x,y}^t\\ C_{\mathrm{Dir}}\left({\▿}_{\mathrm{Dir}}I\right)=\exp (-{\left(\frac{{\▿}_{\mathrm{Dir}}I}{K}\right)}^2),\mathrm{Dir}=N,S,E,W\end{array}\right.$

其中为第t次扩散处理的图像，原图像可表示为，x、y为图像像素的坐标，

为像素在Dir方向上的梯度，即东E、西W、南S和北N。

各向异性扩散算法降噪效果很大程度上决定于算法中对图像的噪声和图像的边界的灰度变化的不同扩散程度的处理。另外，怎样选取一个合适的扩散截止时间也是各向异性扩散算法的一个难题，如果对原图像进行无限制扩散，最终图像将是初始图像灰度平均值的单色图像。

针对超声图像的特点，如何调整扩散参数，何时终止扩散是需要解决的问题。

(三)发明内容

本发明的目的在于提出一种显著提高了图像降噪的快速性的使用边缘信赖度的扩散参数设置方法，并提出了采用非线性相关系数来衡量扩散算法降噪程度的医学超声图像各向异性扩散降噪方法。

本发明是通过以下技术方案实现的：首先计算每个像素的边缘信赖度，并由此确定每个像素的扩散参数取值；然后进行各向异性扩散算法降噪，即根据像素的梯度和扩散参数确定扩散程度，同时更新像素的灰度值；最后计算此次图像扩散降噪结果与扩散前的图像的非线性相关系数，如果该系数低于指定值则图像降噪处理结束，否则回到上一步，即根据像素新的梯度和扩散参数确定扩散程度，并更新像素的灰度值，直至相邻两次图像的非线性相关系数低于指定值。

本发明包括如下技术特征：

第一步，计算每个像素的边缘信赖度，并由此确定每个像素的扩散参数取值：

图像I_x，y的边缘信赖度EC_x，y的计算方法为，将一个局部窗口内的图像像素分布看成一个矢量，将水平和垂直方向的梯度模板看成图像的两个基，图像梯度就是窗口矢量在这两个基上的投影，从而估计像素所在位置的边缘方向。随后通过计算像素矢量与同边缘方向的理想边缘模板的相似程度来估计窗口中心像素的边缘信赖度。边缘信赖度在0和1间取值，越大表示是边缘的可能性越大。若理想边缘模板为A_ref，以像素P为中心的图像块在零均值和归一化处理后用A表示，函数vec()表示将图像块矢量化，则P的边缘信赖度为：

EC_x，y＝|<vec(A)，vec(A_ref)>|

边缘信赖度0.5是边缘和噪声的分界点，即拥有大于0.5的边缘信赖度的像素被认为是边缘，希望扩散参数小些以保持图像的细节，反之边缘信赖度小于0.5的像素被视为噪声，要增大其扩散程度。

定义每个像素的扩散参数：

$\left\{\begin{array}{c}{K}_{x,y}={\mathrm{Ke}}^{4(0.5-{\mathrm{EC}}_{x,y})}\\ {\mathrm{\&lambda;}}_{x,y}=\min (\mathrm{\&lambda;}{e}^{(0.5-{\mathrm{EC}}_{x,y})/2},0.25)\end{array}\right.$

此时K_x，y在0.1353到7.3891倍K之间，λ_x，y在0.7788到1.2840倍λ之间变化，但不大于0.25(这是保证各向异性扩散算法稳定的条件)。这里的K和λ是参数的全局控制参数，决定了参数的变化范围，也是边缘信赖度为0.5的像素的扩散参数值。

第二步，各向异性扩散算法降噪。

利用每个像素的梯度和对应的扩散参数，各向异性扩散算法可改写为：

$\left\{\begin{array}{c}{I}_{x,y}^{t+1}=I_{x,y}^t+{\mathrm{\&lambda;}}_{x,y}{[C_N{\&dtri;}_{N}I+C_S{\&dtri;}_{S}I+C_E{\&dtri;}_{E}I+C_W{\&dtri;}_{W}I]}_{x,y}^t\\ C_{\mathrm{Dir}}=\exp (-{\left(\frac{{\&dtri;}_{\mathrm{Dir}}I}{K_{x,y}}\right)}^2),\mathrm{Dir}=N,S,E,W\end{array}\right.$

第三步，计算

和

的非线性相关系数，若大于指定值则停止，否则转到第二步。

非线性相关系数的处理如下：

对两个离散变量X和Y，变量取值个数为N，变量可取的状态数为b，状态的分布由以下的方式来确定：

1)分别将变量X和Y的数据按从小到大的顺序排列；

2)最小的N/b个值为第一个状态，接下来的N/b个值为第二个状态，如此类推，称每个状态的最小值和最大值为状态阈值；

3)根据状态阈值将变量对(x，y)放到b×b个二维状态格中。

此时对于变量X或Y的状态i，其概率为 $p_i=\frac{N/b}{N}=\frac{1}{b}$ 变量X和Y的联合概率 $p_{\mathrm{ij}}=\frac{n_{\mathrm{ij}}}{N},$ 其中n_ij是在第(i，j)个二维状态格中变量X和Y对的个数。N个数据对在状态格中的分配蕴含着两个变量间统计意义上的普遍相关性。

非线性相关系数定义为：

I^r(X；Y)＝H^r(X)+H^r(Y)-H^r(X，Y)

H^r(X)为改进过的变量X的熵，H^r(X，Y)为改进的变量X和Y联合熵：

$H^r\left(X\right)=-\underset{i=1}{\overset{b}{\mathrm{\&Sigma;}}}{p}_i{\log }_b{p}_i$

$H^r(X,Y)=-\underset{i=1}{\overset{b}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{b}{\mathrm{\&Sigma;}}}{p}_{\mathrm{ij}}{\log }_b{p}_{\mathrm{ij}}$

则非线性相关系数可以重写为：

$I^r(X,Y)=2*(-\underset{i=1}{\overset{b}{\mathrm{\&Sigma;}}}\frac{1}{b}{\log }_b\frac{1}{b})+\underset{i=1}{\overset{b^2}{\mathrm{\&Sigma;}}}{p}_i{\log }_b{p}_i=2+\underset{i=1}{\overset{b}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{b}{\mathrm{\&Sigma;}}}{p}_{\mathrm{ij}}{\log }_b{p}_{\mathrm{ij}}$

指定值一般根据经验取0.6左右。

本发明中，引入边缘信赖度，解决扩散参数调整问题。具体方法为边缘信赖度结合了模型匹配和图像基函数两种思想的概念，将一个局部窗口内的图像像素分布看成一个矢量，将水平和垂直方向的梯度模板看成图像的两个基，图像梯度就是窗口矢量在这两个基上的投影，从而估计像素所在位置的边缘方向。随后通过计算像素矢量与同边缘方向的理想边缘模板的相似程度来估计窗口中心像素的边缘信赖度。边缘信赖度在0和1间取值，越大表示是边缘的可能性越大。与梯度相比，边缘信赖度受噪声的影响较小，且基本能指示出正确的边缘位置和强度。

引入非线性相关系数，解决扩散何时终止问题。具体方法为线性相关系数基于互信息的定义，以[0，1]之间的值来表征两个变量间非线性相关程度的强弱，0表示最弱，1表示最强。本发明将针对医学超声图像的特点，结合边缘信赖度来调整各向异性扩散算法中参数的设置，使每一个像素在初始时就处于不同的扩散条件中，增强了各向异性扩散算法的局部适应性，并利用非线性相关系数来衡量相邻两幅图像的相似性，并以此作为扩散算法的停止标准，避免了图像被过度平滑。

结合图1-图4，本发明在各向异性扩散算法中增加了反映图像局部信息的基于边缘信赖度的扩散参数设置使得降噪结果获得明显提高，并提出了基于非线性相关系数的扩散算法停止判定，使算法可根据运算结果自动停止。本发明给出了扩散参数调整具体方法和扩散终止判别方法，提高了超声图像降噪的各向异性扩散算法的实用性，克服了以往的各向异性扩散算法中参数设置适应性不强的问题。

(四)附图说明

图1-图4为本发明的实验结果，其中图1为带有散斑噪声的医学超声原图像，图2为对图像进行本发明的方法处理后的边缘提取结果，图3为对图像进行各向异性扩散算法处理后的边缘提取结果，图4为对图像进行中值滤波处理后的边缘提取结果。

图5为本发明的过程框图。

(五)具体实施方式

下面结合附图对本发明作进一步描述：

结合图5，本实施例的步骤为：第一步，计算每个像素的边缘信赖度，并由此确定每个像素的扩散参数取值。第二步，各向异性扩散算法降噪。第三步，计算

和

的非线性相关系数，若大于指定值则停止，否则转到第二步。

结合图1，基于本发明的扩散方法和各向异性扩散算法的K和λ相同，扩散算法的参数设置为K＝30，λ＝0.15，非线性相关系数的阈值为0.61。由于边缘信赖度对噪声不敏感，只在扩散算法开始时计算边缘信赖度，在循环过程中不再计算，这样方法的计算量只增加了少许。本发明的扩散方法处理过程中，非线性相关系数的指定值为0.61，迭代1次后自动停止。常规各向异性扩散算法计算过程中非线性相关系数＝0.59、0.62，迭代2次后停止。从而可以看出，本发明的处理方法降噪速度比常规各向异性扩散算法快。比较边缘检测出来的结果，本发明的扩散方法在保留了边缘的前提下，降噪能力好于常规各向异性扩散算法和中值滤波算法。

综上所述，对于医学超声图像的降噪，本发明的处理方法能够得到更好的结果图像。

Claims (1)

1、一种医学超声图像各向异性扩散降噪方法，其特征在于，首先计算每个像素的边缘信赖度，并由此确定每个像素的扩散参数取值；然后进行各向异性扩散算法降噪，即根据图像梯度和扩散参数确定扩散程度，同时更新像素的灰度值；最后计算此次图像扩散降噪结果与扩散前的图像的非线性相关系数，如果该系数低于指定值则图像降噪处理结束，否则回到上一步，即根据像素新的梯度和扩散参数确定扩散程度，并更新像素的灰度值，直至相邻两次图像的非线性相关系数低于指定值；

所述的计算每个像素的边缘信赖度的方法为：假设图像I_x，y的边缘信赖度为EC_x，y，理想边缘模板为A_ref，以像素P为中心的图像块在零均值和归一化处理后用A表示，函数vec()表示将图像块矢量化，则像素P的边缘信赖度为：EC_x，y＝|<vec(A)，vec(A_ref)>|；

所述的确定每个像素的扩散参数取值的方法为：定义每个像素的扩散参数K和λ：

$\left\{\begin{array}{c}{K}_{x,y}={\mathrm{Ke}}^{4(0.5-{\mathrm{EC}}_{x,y})}\\ {\mathrm{\&lambda;}}_{x,y}=\min ({\mathrm{\&lambda;e}}^{(0.5-{\mathrm{EC}}_{x,y})/2},0.25)\end{array}\right.$

K_x，y在0.1353到7.3891倍K之间，λ_x，y在0.7788到1.2840倍λ之间变化，但不大于0.25；

所述的非线性相关系数的计算步骤为：

对两个离散变量X和Y，变量取值个数为N，变量可取的状态数为b，状态的分布由以下的方式来确定：

1)分别将变量X和Y的数据按从小到大的顺序排列；

2)最小的N/b个值为第一个状态，接下来的N/b个值为第二个状态，如此类推，称每个状态的最小值和最大值为状态阈值；

3)根据状态阈值将变量对(x，y)放到b×b个二维状态格中，此时对于变量X或Y的状态i，其概率为 $p_i=\frac{N/b}{N}=\frac{1}{b},$ 变量X和Y的联合概率 $p_{\mathrm{ij}}=\frac{n_{\mathrm{ij}}}{N},$ 其中n_ij是在第(i，j)个二维状态格中变量X和Y对的个数，非线性相关系数定义为：

$I^r(X,Y)=2*(-\underset{i=1}{\overset{b}{\mathrm{\&Sigma;}}}\frac{1}{b}{\log }_b\frac{1}{b})+\underset{i=1}{\overset{b^2}{\mathrm{\&Sigma;}}}{p}_i{\log }_b{p}_i=2+\underset{i=1}{\overset{b}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{b}{\mathrm{\&Sigma;}}}{p}_{\mathrm{ij}}{\log }_b{p}_{\mathrm{ij}}.$

Images