CN102073994A - 基于多尺度各项异性扩散的超声医学图像散斑噪声抑制方法 - Google Patents

Abstract

基于多尺度各项异性扩散的超声医学图像散斑噪声抑制方法，它涉及一种医学超声图像噪声抑制方法。它为解决将现有各项异性扩散技术应用于超声图像时，用梯度作为边缘检测算子难以准确区分图像边缘和噪声的问题而提出。一：对图像进行小波分解，计算标准模值；二：若为首次迭代，估计标准模值的分布参数，否则跳到步骤四；三：分类得到各尺度下的噪声平稳区域；四：用噪声区域内标准模值的均值确定扩散阈值；并代入扩散函数得到扩散参数；五：用扩散参数对小波系数进行加权后，进行小波逆变换；六：比较相邻两次迭代中的扩散阈值，判断是否终止。本发明在多尺度下对图像进行处理，结果更加精确，可广泛适用于各种需要对超声图像进行处理的场合。

Description

基于多尺度各项异性扩散的超声医学图像散斑噪声抑制方法

技术领域

本发明涉及一种超声图像散斑噪声抑制方法。

背景技术

超声成像技术现广泛用于对诸如心脏的体内结构进行成像。由于成像的物理因素，相干波的叠加形成散斑噪声，严重影响了图像质量，尤其是使边缘和细节变得模糊，降低了图像的可观性和分割精度，限制了超声医学图像在临床诊断中的进一步应用。对超声图像进行降噪处理，可以为图像分割、斑点追踪等后续处理提供更有利的条件。因此，研究在保持图像边缘和细节信息的基础上有效抑制散斑噪声，有助于提高临床诊断的准确性，具有重要意义。

从目前在该领域内对超声医学图像进行噪声抑制的方法来看，主要有空域滤波、小波变换和各向异性扩散算法等。空域滤波主要指邻域平均法和中值滤波等技术，其本质是根据像素点邻域窗口内的信息确定一个新值来取代原值，在抑制噪声的同时往往会引起图像的模糊和重要细节的丢失，限制了处理的效果。基于小波变换的方法将超声图像变换到小波域，通过阈值方法去除某些小波系数，再逆变换以去除噪声，该方法主要用于加性噪声的去除，对于如散斑噪声等乘性噪声阈值难以确定。基于各向异性扩散方程的降噪技术使用选择性扩散方式，在图像的噪声处有较大的扩散系数，有利于图像的平滑，而在边缘处有较小的扩散系数，保持了图像的细节，但超声图像中，由于乘性噪声的影响，难以从灰度变化上准确的区分图像的边缘和噪声，加大了算法的难度。

各项异性扩散的过程可通过小波变换的方式来实现，过程如图7所示，将分解得到的小波系数逐一加权后再进行逆变换输出结果，相当于在多个尺度下进行各向异性扩散。该方法能够在多尺度下对图像进行处理，结果更加精确，同时继承了各向异性扩散能够保持图像细节的特点，不失为超声图像处理的好办法。该方法中，如何针对超声医学图像的特点选用合适的边缘检测器，以及如何调整扩散参数是关系到降噪质量好坏的最关键问题。

发明内容

本发明为了解决现有基于各向异性扩散方程的降噪技术使用梯度检测边缘的方法易受乘性噪声影响，难以从灰度变化上准确区分图像边缘和噪声的问题，而提出的基于多尺度各项异性扩散的超声医学图像散斑噪声抑制方法。

基于多尺度各项异性扩散的超声医学图像散斑噪声抑制方法，所述方法由如下步骤完成：

步骤一：使用Mallat定义的二次样条小波滤波器对图像进行非抽样小波分解，得到三个尺度下的近似分量，水平方向的小波系数和竖直方向的小波系数，计算小波系数标准模值M′_jf；

步骤二：判断是否为首次迭代；是，则用混合瑞利模型来模拟图像扇形扫描区内小波系数标准模值M′_jf的分布，并设定出每个尺度下模型参数的值，执行步骤三；否，则跳到步骤四；

步骤三：使用极大似然分类的方法得出各尺度下的噪声平稳区域；

步骤四：通过计算每个尺度下噪声区域内M′_jf的均值确定扩散阈值；并将该扩散阈值代入扩散函数求得扩散参数p_j；

步骤五：用步骤四求得的扩散参数p_j对步骤一分解得到的小波系数进行加权计算，然后对图像进行小波逆变换，完成一次多尺度各向异性扩散过程；

步骤六：比较相邻两次迭代过程中的扩散阈值，判断它们之间相差的百分比是否小于阈值th；是，则抑制完成；否，则返回步骤一进行下次迭代。

本发明排除了乘性噪声影响，可以从图像灰度变化上更加准确地区分图像边缘和噪声。本发明所述方法在多尺度下对图像进行处理，结果更加精确，同时继承了各向异性扩散能够保持图像细节的特点。本发明可广泛适用于各种需要对超声图像进行处理的场合。

附图说明

图1为本发明所述方法的流程图；图2为实际图像中相同尺度下标准模值的分布直方图与混合瑞利分布模型曲线的比较示意图，图中实心点为标准模值的分布直方图，曲线为混合瑞利分布概率密度曲线；图3为扩散阈值随迭代次数的变化情况；图4为原始未处理的超声图像；图5为采用本发明所述方法当k＝0.5时输出的图像；图6为采用本发明所述方法当k＝2时输出的图像；图7为用小波变换实现一次多尺度各向异性扩散的示意框图。

具体实施方式

具体实施方式一：结合图1至图7至说明本实施方式，所述方法由如下步骤完成：

步骤一：使用Mallat定义的二次样条小波滤波器对图像进行非抽样小波分解，得到三个尺度下的近似分量，水平方向的小波系数和竖直方向的小波系数，计算小波系数标准模值M′_jf；

步骤二：判断是否为首次迭代；是，则用混合瑞利模型来模拟图像扇形扫描区内小波系数标准模值M′_jf的分布，并设定出每个尺度下模型参数的值，执行步骤三；否，则跳到步骤四；

步骤三：使用极大似然分类的方法得出各尺度下的噪声平稳区域；

步骤四：通过计算每个尺度下噪声区域内M′_jf的均值确定扩散阈值；并将该扩散阈值代入扩散函数求得扩散参数p_j；

步骤五：用步骤四求得的扩散参数p_j对步骤一分解得到的小波系数进行加权计算，然后对图像进行小波逆变换，完成一次多尺度各向异性扩散过程；

步骤六：比较相邻两次迭代过程中的扩散阈值，判断它们之间相差的百分比是否小于阈值th；是，则抑制完成；否，则返回步骤一进行下次迭代。

步骤一中小波系数标准模值M′_jf采用如下算法得到：

首先计算小波系数模：

$M_{j}f=\sqrt{{\left(W_j^{1}f\right)}^2+{\left(W_j^{2}f\right)}^2}$ 公式1

式中参数和

分别代表j尺度下水平方向的小波系数和竖直方向的小波系数；然后计算标准模值：

$M_j^{\′}f=\frac{M_{j}f}{\sqrt{{\mathrm{\μ}}_s}}$ 公式2

式中μ_s是图像中一活动窗口内的灰度均值，各尺度下窗的大小定义为D_j＝2^j-1(D₀-1)+1，参数D₀＝3为j＝1时窗的大小。

步骤一采用的小波滤波器的系数参加下表：

步骤二中步骤二中用混合瑞利模型来模拟图像扇形扫描区内小波系数标准模值M′_jf分布的方法为：边缘和噪声处的标准模值M′_jf分别服从方差为和

的瑞利分布，概率密度函数分别由p_e(x)和p_n(x)表示：

$p_e\left(x\right)=\frac{x}{{\mathrm{\σ}}_e^2}\exp (-\frac{x^2}{2{\mathrm{\σ}}_e^2}),$ $p_n\left(x\right)=\frac{x}{{\mathrm{\σ}}_n^2}\exp (-\frac{x^2}{2{\mathrm{\σ}}_n^2})$

则图像扇形扫描区内的M′_jf服从混合瑞利分布：

p(x)＝ω_np_n(x)+(1-ω_n)p_e(x)             公式3

式中参数

和ω_n的值通过期望最大化方法估计得出。实际分布曲线与标准混合瑞利模型曲线如图2所示。

步骤三使用极大似然分类的方法得出各尺度下的噪声平稳区域的方法如下：图像扫描区可分为两部分：细节区域和噪声区域；利用极大似然分类，错误概率最小时得到分类阈值：

$T=\sqrt{\frac{2(\log \frac{{\mathrm{\σ}}_n^2}{{\mathrm{\σ}}_e^2}+\frac{{\mathrm{\ω}}_n}{1-{\mathrm{\ω}}_n})}{|\frac{1}{{\mathrm{\σ}}_e^2}-\frac{1}{{\mathrm{\σ}}_n^2}|}}$ 公式4

按照从大尺度到小尺度的方法对图像进行分类：

公式5

当参数U＝1的区域表示噪声平稳区域，且在最大尺度j＝J下认为U_J＝1；K为可调节参数，K较小时(如K＝0.5)可用于增强图像可读性，K较大时(如K＝2)可用于后续分割。

步骤四中扩散参数p_j的计算方法如下：选择非线性扩散函数：

$c\left(\mathrm{\η}\right)=\left\{\begin{array}{cc}1,& \mathrm{\η}\≤0\\ 1-\exp [-\frac{-3.315}{{(\mathrm{\η}/\mathrm{\λ})}^4}],& \mathrm{\η}>0\end{array}\right.$ 公式6

令η＝M′_jf；扩散阈值λ_j取为对应尺度下噪声平稳区域内标准模值M′_jf的均值与系数的加权：

${\mathrm{\λ}}_j=\frac{\mathrm{Mean}\left(U_j{M}_j^{\′}f\right)}{{\sqrt{2}}^{j^{\′}}},$ $\left\{\begin{array}{cc}{j}^{\′}=0,& j=1\\ j^{\′}=j,& j\≤2\end{array}\right.$ 公式7

扩散参数p_j由下式确定：

p_j＝1-c_j(η)              公式8

步骤五中通过扩散参数p_j对小波系数进行加权计算的方法如下：

将小波系数W_jf乘以扩散参数p_j，得到新的小波系数W′_jf；

$W_j^{\′1}f=p_j\·W_j^{1}f,$ $W_j^{\′2}=f=p_j\·W_j^{2}f$ 公式9。

随着迭代次数的增加，扩散阈值λ逐渐趋于稳定，比较前后两次迭代过程中的扩散阈值，相差百分比小于指定阈值th时停止迭代；设λⁿ和λⁿ⁺¹分别为由公式7算得的第n次和第n+1次迭代中的扩散阈值，当所有尺度下的λ_j(j＝1，…，J)均满足以下条件时迭代停止：

$\frac{|{\mathrm{\&lambda;}}^{n+1}-{\mathrm{\&lambda;}}^n|}{{\mathrm{\&lambda;}}^n}<\mathrm{th}$ 公式10

阈值th一般取0.5％，λ在各尺度下随迭代次数的变化情况如图3所示。

本发明采用小波系数标准模值M′_jf作为超声图像边缘检测器的原理：

原始的散斑噪声是一种乘性噪声，幅值具有瑞利分布的特性。在医学超声成像系统中，为改善质量图像在输出前往往要经过同态滤波等预处理，输出图像数学模型可表示如下：

$f\left(x\right)=g\left(x\right)+\sqrt{g\left(x\right)}n\left(x\right)$ 公式11

其中，参数g(x)为原始无噪图像。设图像中一移动窗口内原始灰度均值为μ，则该窗口内原始图像可定义为g(x)＝μR(x)，则公式11可写为：

$f\left(x\right)=\mathrm{\&mu;R}\left(x\right)+\sqrt{\mathrm{\&mu;R}\left(x\right)}n\left(x\right)$ 公式12

对上式进行小波变换：

公式13

窗口很小时可认为R(x)＝1，由小波变换性质易知

则上式中等号右侧第一项为零，即

故将小波系数模值除以局部窗内灰度均值可去除原始图像自身灰度的影响。同理，小波系数模值除以局部窗内灰度均值得到的标准模值，去除了原始图像自身灰度的影响，使它相对于现有的使用梯度检测边缘的方法更准确更有效，成为对超声医学图像更准确更有效的边缘检测器。

以上内容是结合具体的优选实施方式对本发明所作的进一步详细说明，不能认定本发明的具体实施只局限于这些说明。对于本所属技术领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干简单推演或替换，都应当视为属于本发明所提交的权利要求书确定的专利保护范围。

Claims (6)

1.基于多尺度各项异性扩散的超声医学图像散斑噪声抑制方法，其特征在于所述方法由如下步骤完成：

步骤一：使用Mallat定义的二次样条小波滤波器对图像进行非抽样小波分解，得到三个尺度下的近似分量，水平方向的小波系数和竖直方向的小波系数，计算小波系数标准模值M′_jf；

步骤二：判断是否为首次迭代；是，则用混合瑞利模型来模拟图像扇形扫描区内小波系数标准模值M′_jf的分布，并设定出每个尺度下模型参数的值，执行步骤三；否，则跳到步骤四；

步骤三：使用极大似然分类的方法得出各尺度下的噪声平稳区域；

步骤四：通过计算每个尺度下噪声区域内M′_jf的均值确定扩散阈值；并将该扩散阈值代入扩散函数求得扩散参数p_j；

步骤五：用步骤四求得的扩散参数p_j对步骤一分解得到的小波系数进行加权计算，然后对图像进行小波逆变换，完成一次多尺度各向异性扩散过程；

步骤六：比较相邻两次迭代过程中的扩散阈值，判断它们之间相差的百分比是否小于阈值th；是，则抑制完成；否，则返回步骤一进行下次迭代。

2.根据权利要求1所述的基于多尺度各项异性扩散的超声医学图像散斑噪声抑制方法，其特征在于步骤一中归一化小波系数标准模值M′_jf采用如下算法得到：

首先计算小波系数模：

$M_{j}f=\sqrt{{\left(W_j^{1}f\right)}^2+{\left(W_j^{2}f\right)}^2}$ 公式1

式中参数

和

分别代表j尺度下水平方向的小波系数和竖直方向的小波系数；然后计算标准模值：

$M_j^{\&prime;}f=\frac{M_{j}f}{\sqrt{{\mathrm{\&mu;}}_s}}$ 公式2

式中μ_s是图像中一活动窗口内的灰度均值，各尺度下窗的大小定义为D_j＝2^j-1(D₀-1)+1，参数D₀＝3为j＝1时窗的大小。

3.根据权利要求2所述的基于多尺度各项异性扩散的超声医学图像散斑噪声抑制方法，其特征在于步骤二中用混合瑞利模型来模拟图像扇形扫描区内小波系数标准模值M′_jf分布的方法为：边缘和噪声处的标准模值M′_jf分别服从方差为

和

的瑞利分布，概率密度函数分别由p_e(x)和p_n(x)表示：

$p_e\left(x\right)=\frac{x}{{\mathrm{\&sigma;}}_e^2}\exp (-\frac{x^2}{2{\mathrm{\&sigma;}}_e^2}),$ $p_n\left(x\right)=\frac{x}{{\mathrm{\&sigma;}}_n^2}\exp (-\frac{x^2}{2{\mathrm{\&sigma;}}_n^2})$

则图像扇形扫描区内的M′_jf服从混合瑞利分布：

p(x)＝ω_np_n(x)+(1-ω_n)p_e(x)          公式3

式中参数

和ω_n的值通过期望最大化方法估计得出。

4.根据权利要求3所述的基于多尺度各项异性扩散的超声医学图像散斑噪声抑制方法，其特征在于步骤三使用极大似然分类的方法得出各尺度下的噪声平稳区域的方法如下：利用极大似然分类，错误概率最小时得到分类阈值：

$T=\sqrt{\frac{2(\log \frac{{\mathrm{\&sigma;}}_n^2}{{\mathrm{\&sigma;}}_e^2}+\frac{{\mathrm{\&omega;}}_n}{1-{\mathrm{\&omega;}}_n})}{|\frac{1}{{\mathrm{\&sigma;}}_e^2}-\frac{1}{{\mathrm{\&sigma;}}_n^2}|}}$ 公式4

然后，按照从大尺度到小尺度的方法对图像进行分类：

公式5

其中，参数U＝1的区域表示噪声平稳区域。

5.根据权利要求4所述的基于多尺度各项异性扩散的超声医学图像散斑噪声抑制方法，其特征在于步骤四中扩散参数p_j的求得方法如下：选择非线性扩散函数：

$c\left(\mathrm{\&eta;}\right)=\left\{\begin{array}{cc}1,& \mathrm{\&eta;}\&le;0\\ 1-\exp [-\frac{-3.315}{{(\mathrm{\&eta;}/\mathrm{\&lambda;})}^4}],& \mathrm{\&eta;}>0\end{array}\right.$ 公式6

令η＝M′_jf；扩散阈值λ_j取为对应尺度下噪声平稳区域内标准模值M′_jf的均值与系数的加权：

${\mathrm{\&lambda;}}_j=\frac{\mathrm{Mean}\left(U_j{M}_j^{\&prime;}f\right)}{{\sqrt{2}}^{j^{\&prime;}}},$ $\left\{\begin{array}{cc}{j}^{\&prime;}=0,& j=1\\ j^{\&prime;}=j,& j\&le;2\end{array}\right.$ 公式7

扩散参数p_j由下式确定：

p_j＝1-c_j(η)               公式8。

6.根据权利要求5所述的基于多尺度各项异性扩散的超声医学图像散斑噪声抑制方法，其特征在于步骤五中通过扩散参数p_j对小波系数进行加权计算的方法如下：

将小波系数W_jf乘以扩散参数p_j，得到新的小波系数W′_jf；

$W_j^{\&prime;1}f=p_j\&CenterDot;W_j^{1}f,$ $W_j^{\&prime;2}=f=p_j\&CenterDot;W_j^{2}f$ 公式9。

Images