CN101452574B - 一种提升超声图像质量的方法 - Google Patents

Abstract

本发明提供一种提升超声图像质量的方法，用于超声成像系统中对超声扫描图像显示数据的优化处理，包括步骤：采用高斯-拉普拉斯金字塔分解方法将超声图像数据分解若干层，对分解结果中的高斯层数据进行各向异性扩散滤波处理，然后对处理后的超声图像数据逆向重构，经过若干次的迭代处理，提升了超声图像质量。其中，各向异性扩散滤波算法根据超声图像的局部结构信息计算出不同方向上的扩散系数，对图像数据进行滤波操作。结果是图像的组织结构区域信息得到增强，而非组织结构区域噪声得到有效抑制。

Description

一种提升超声图像质量的方法

技术领域

本发明涉及超声成像技术领域，特别是超声成像技术中的数字图像处理技术，提升超声图像质量的方法。

背景技术

一个典型的超声成像系统如图1所示，它包括一个超声波成像装置，主控制器控制整个系统的正常运转。在主控制器的控制下，探头发射出超声波，延时一定时间后接收从组织反射回来的超声波。波束形成器将各通道反射回波信号进行聚焦延迟、加权、求和，形成一条或多条扫描线，再经过滤波、检波等处理获得扫描回波的包络信号。回波包络送到数字扫描变换模块完成坐标变换，形成图像数据。所形成的图像数据再送到显示模块，叠加上字符后输出到显示器显示。

由于超声束的空间分辨力有限，组织的反射面不光滑，不同的反射源产生的回波信号相位也不一样，在波束形成的时候可能相叠加也可能相消，及电子电路的噪声等因素导致超声图像中存在许多斑点噪声。斑点噪声的存在将会掩盖超声图像中的部分有用信息，对医生的诊断造成干扰，影响医生的视觉感受甚至影响到医生的诊断信心。

有很多的方法抑制斑点噪声，例如多帧图像平均，扫描线之间的平滑，图像像素的邻域平均等。另外为了突出有用的组织结构特征，有的超声成像设备采用锐化图像的方式，使得超声图像中的组织结构、边缘信息得到加强。

上述常用的办法虽然抑制了噪声或增强了组织结构特征，但在抑制了斑点噪声的同时，却模糊了组织结构或损失了动态特性，或在增强了组织结构信息的同时，却将噪声也一起加强了，无法真正提高超声图像质量。

与解决上述问题相关的是，GE公司在美国专利US 6,208,763和US 6,592,523中公开了一种基于梯度分割的图像增强方法。该方法依据梯度信息，将图像分割为结构区域与非结构区域。对结构区域进行基于灰度加权方向性二阶倒数的各向异性锐化处理，以增强图像边缘的对比度；对非结构区域进行各向同性的平滑处理以消除斑点噪声。

深圳迈瑞生物医疗电子股份有限公司在中国专利200510035913.8和200510036841.9中公开了一种基于梯度分割的图像增强方法。该方法依据梯度信息，将图像分割为结构区域与非结构区域。对结构区域进行基于中间像素的加权平均各向异性增强处理；对非结构区域进行各向同性的平滑处理以消除斑点噪声。

发明内容

上述方法对特征区域处理时，不能建立有效的数学模型增强图像的边缘信息和结构信息。本发明要解决的就是上述方法的不足，即利用热扩散方程模型，达到在抑制斑点噪声的同时，不损失反而增强超声图像中组织结构信息的目的，从而提升超声诊断图像质量。

本发明的目的在于克服现有技术的缺点，提供一种提升超声图像质量的方法。

为了在抑制斑点噪声的同时增强组织结构特征，本发明的基本思想：在现有超声图像基础上，通过后处理的方式，利用由结构张量表示的图像像素的局部特征，将图像划分成组织结构区域，不含组织结构的平坦区域以及组织结构的边缘区域。对三类区域分别采用不同的处理方式，从而达到抑制斑点噪声的同时增强组织结构信息。为了抑制各种大小范围的斑点噪声，增强不同尺度的组织结构特征，所述的图像处理分别在多个尺度范围内进行。

为实现上述目的，本发明提出了一种提升超声图像质量方法，该方法增强超声图像结构信息与抑制斑点噪声，用于超声成像系统对扫描图像显示数据的优化处理。其的步骤如下：

A.利用高斯-拉普拉斯金字塔分解算法分解超声显示图像；

B.对分解产生的高斯层图像数据进行各向异性扩散滤波；

C.对扩散结果进行锐化处理；

D.利用高斯-拉普拉斯金字塔重构算法将各向异性扩散滤波结果与拉普拉斯层数据重构回去；

E.A、B、C、D三个步骤迭代若干次。

金字塔分解可以获得不同分辨尺度下原始图像的表达，金字塔重构可以把不同分辨尺度下的图像几乎没有任何信息丢失重构回去。本发明正是基于这种思想，将图像分解到不同的尺度，然后在不同的尺度下对图像做扩散滤波处理，将不同尺度下处理过的图像重构回去，从而获得最终的结果。

各向异性扩散是基于热传导模型，实现对图像的平滑的同时可以保持甚至增强图像的边界或结构信息。采用高斯核与输入图像做卷积的图像滤波方法，等同于扩散系数为常值的热传导方程的解，即各向同性扩散。由于这种基于高斯核的图像滤波方法没有考虑图像特征在空间上的分布，因此在平滑图像的同时也滤除了图像的一些重要的特征信号。本发明考虑到各向同性扩散滤波的缺点，采用各向异性扩散滤波，根据图像的局部特征计算扩散系数。图像的局部特征用结构张量的特征向量和特征值来表征。特征值间的差值越大，表明这个区域的结构特征越强，大特征值方向采用较小的扩散系数，而小特征值方向采用较大的扩散系数；反之如果特征值间的差值越小，表明这个区域有结构信息的可能性就越小，即越平坦，两个特征方向的扩散系数基本相同。如此一来，扩散结果就是，有结构信息的地方沿着结构边缘扩散，而垂直于结构的边缘方向则几乎没有扩散，从而实现了结构信息的保留甚至增强。对于扩散方程的求解，本发明采用托马斯求解方法，以最高的效率完成扩散滤波算法，以便于实现实时处理超声扫描显示图像。

上述方案中，所述步骤B，对由输入图像分解出来的高斯层进行4个方向的各向异性扩散滤波，由以下几个步骤组成：

a.计算4个方向扩散系数，4个方向分别是直角坐标系中的-45°、0°、45°和90°方向；

b.在上述4个方向进行扩散滤波。

所述过程a，由以下几个步骤组成：

1.高斯平滑；

2.计算梯度；

3.计算结构张量；

4.计算扩散张量；

5.裂解为4个方向的扩散系数。

有益效果

本发明达到在抑制斑点噪声的同时，不损失反而增强超声图像中组织结构信息的目的，从而提升超声诊断图像质量。

附图说明

图1是典型的超声成像系统图。

图2是超声成像系统中本发明的图像处理模块图。

图3是提升超声图像质量的方法总体流程图。

图4是高斯-拉普拉斯金字塔分解流程图。

图5是高斯-拉普拉斯金字塔重构流程图。

图6是各向异性扩散滤波算法流程图。

图7是进行扩散的4个方向示意图。

具体实施方式

图1是典型的超声成像系统图。图中由主控制器、探头、波束形成器、数字扫描变换(DSC，Digital Scan Conversion)、显示器组成，在主控制器的控制下，探头发射出超声波，延时一定时间后接收从组织反射回来的超声波。波束形成器将各通道反射回波信号进行聚焦延迟、加权、求和，形成一条或多条扫描线，再经过滤波、检波等处理获得扫描回波的包络信号。回波包络送到数字扫描变换模块完成坐标变换，形成图像数据。所形成的图像数据再送到显示模块，叠加上字符后输出到显示器显示。

图2是超声成像系统中本发明的图像处理模块图。

本发明提升超声图像质量的方法，抑制超声图像中的斑点噪声，增强超声图像中组织结构信息，用于在超声成像系统中对超声扫描图像显示数据的优化处理，其在系统中的位置如图2所示。是在DSC与显示器之间加入一个图像处理模块，图中虚线所示。该图像处理模块是提升超声图像质量方法的核心，主要包括：图像数据输入接口，处理参数输入接口，图像显示数据输出接口，高斯-拉普拉斯金字塔分解重构算法，结构张量求解算法，扩散张量求解算法，各向异性扩散滤波算法，加性裂解算法，数字图像锐化算法。

图3是提升超声图像质量的方法总体流程图。

本发明提升超声图像质量的方法，其总体实现如图3所示。首先把将要送出显示的图像数据利用高斯-拉普拉斯金字塔分解算法，将其分解成若干层。接着对数据量最小的高斯层进行各向异性扩散滤波处理，再对处理结果做锐化处理。将处理过的高斯层数据与未处理的拉普拉斯层数据利用高斯-拉普拉斯重构算法进行重构。以上四个步骤经过若干次循环即可实现提升超声图像质量的目标，图像斑点噪声得到有效抑制，同时组织结构信息得到增强。

图4、图5是高斯-拉普拉斯金字塔分解、重构流程图。

本发明为了减少整个算法的计算量，以及实现在多尺度下对图像进行处理，对输入图像进行了高斯-拉普拉斯金字塔分解，将图像中的能量集中于数据量较少的高斯层数据。所述高斯-拉普拉斯金字塔分解与重构算法，为PETER J.BURT等在1983年提出来的。分解算法的流程如图4所示，重构算法的流程如图5所示。

高斯金字塔或拉普拉斯金字塔的分解和重构均由两个重要的算子组成，分别为Reduce和Expand，Reduce是数据经过数字滤波器后再下采样的过程，Expand与Reduce相反，是对数据上采样再滤波的过程。

分解过程中，高斯金字塔的求法：

Gauss₀＝x；

Gauss₁＝Reduce(Gauss_1-1)。

而拉普拉斯金字塔的求法：

Laplace₁＝Gauss₁-Expand(Gauss₁₊₁)

重构算法：

Gauss₁＝Expand(Gauss₁₊₁)+Laplace₁

其中x为输入阵列，Gauss₁为第1层高斯层，Laplace₁为第L层拉普拉斯层。

本发明在金字塔分解算法中的滤波器，采用若干阶FIR滤波器。下采样前的滤波器和上采样后的滤波器具有特定的关系。

图6、图7是各向异性扩散滤波算法、进行扩散流程图。

本发明通过高斯-拉普拉斯金字塔分解算法对输入图像分解后，对最底层的高斯层进行各向异性扩散滤波，目的之一是减少计算量，目的之二是不同层数的金字塔分解，其最底层高斯层尺度也不同。因而通过多次循环可以实现在不同尺度下处理图像。本发明采用四个方向扩散滤波，分别是直角坐标系中的-45°、0°、45°和90°方向。而要进行各向异性扩散滤波，必须求出相应每个像素的4个方向的扩散系数。

本发明提升超声图像质量的方法，为了抑制超声图像中的噪声，增强组织结构信息，在计算扩散系数的时候将图像的局部特征考虑进去。本发明采用结构张量来表示局部的特征，某个像素的结构张量为 $\left(\begin{array}{cc}{I}_x^2& I_x{I}_y\\ I_x{I}_y& I_y^2\end{array}\right),$ 其中，I_x为x方向的梯度，I_y为y方向的梯度。结构张量可以分解为两个正交的向量。利用这两个向量的方向和特征值(u₁，u₂)来表征局部特性。两个特征值的差值越大，说明这个区域的结构特征越强，大特征值方向的扩散小，小特征值方向的扩散大，反之如果两个特征值的差值越小，说明这个地方越没有结构信息，即越平坦，两个特征方向的扩散可基本相同。如此一来，扩散结果就是，有结构信息的地方沿着结构边缘有扩散，而垂直于结构的边缘方向则几乎没有扩散，从而实现了结构信息的保留甚至增强。

结构张量的计算要考虑到噪声的影响，因此在计算梯度的时候，首先采用若干阶FIR滤波器平滑图像，平滑之后的图像噪声得到抑制。基于平滑之后的图像求出每个像素的x，y方向的梯度，获得结构张量的组成部分，即I_x ²，I_y ²，I_xI_y，并分别对这三个数据组成的阵列再次作平滑处理，进一步减小噪声对整个算法的影响，得到最终需要的结构张量 $\left(\begin{array}{cc}{j}_{11}& j_{12}\\ j_{12}& j_{22}\end{array}\right).$

结构张量 $\left(\begin{array}{cc}{j}_{11}& j_{12}\\ j_{12}& j_{22}\end{array}\right)$ 的特征值和特征向量为：

${\mathrm{\μ}}_1=\frac{1}{2}(j_{11}+j_{22}+\sqrt{{(j_{11}-j_{22})}^2+{4j}_{12}^2}),{\mathrm{\μ}}_2=\frac{1}{2}(j_{11}+j_{22}-\sqrt{{(j_{11}-j_{22})}^2+{4j}_{12}^2}),$

$\stackrel{\→}{{\mathrm{\ω}}_1}\left|\right|\left(\begin{array}{c}{2j}_{12}\\ j_{22}-j_{11}+\sqrt{{(j_{11}-j_{22})}^2+{4j}_{12}^2}\end{array}\right),$

$\stackrel{\→}{{\mathrm{\ω}}_2}\left|\right|\left(\begin{array}{c}{2j}_{12}\\ j_{22}-j_{11}-\sqrt{{(j_{11}-j_{22})}^2+{4j}_{12}^2}\end{array}\right).$

扩散张量 $\left(\begin{array}{cc}{d}_{11}& d_{12}\\ d_{12}& d_{22}\end{array}\right)$ 的计算：由结构张量的特征值，计算出扩散张量的特征值λ₁和λ₂，再由新的特征值结合结构张量的特征向量计算出扩散张量。

${({\mathrm{\μ}}_1-{\mathrm{\μ}}_2)}^2={(j_{11}-j_{22})}^2+4*j_{12}^2$

dd＝(λ₂-λ₁)(j₁₁-j₂₂)/(μ₁-μ₂)

d₁₁＝0.5*[(λ₂+λ₁)+(λ₂-λ₁)(j₁₁-j₂₂)/(μ₁-μ₂)]

d₁₂＝(λ₁-λ₂)j₁₂/(μ₁-μ₂)

d₂₂＝0.5*[(λ₂+λ₁)-(λ₂-λ₁)(j₁₁-j₂₂)/(μ₁-μ₂)]

4个扩散系数的获得。再由扩散张量转化为4个方向( ${\mathrm{\β}}_{-1}=-\frac{\mathrm{\π}}{4},$ β0＝0， ${\mathrm{\β}}_2=\frac{\mathrm{\π}}{4},$ ${\mathrm{\β}}_2=\frac{\mathrm{\π}}{2}$ )的扩散系数：

${\mathrm{\α}}_{-1}=\frac{\left|d_{12}\right|-d_{12}}{2},$

α₀＝d₁₁-|d₁₂|，

${\mathrm{\α}}_1=\frac{\left|d_{12}\right|+d_{12}}{2},$

α₂＝d₂₂-|d₁₂|。

扩散滤波采用托马斯算法，以达到处理的高效率。

本发明的一次处理循环包括两次各项异性扩散滤波，第一次的扩散系数计算方法为：

${\mathrm{\λ}}_2=\left\{\begin{array}{cc}\mathrm{\α}\·(1-\frac{{({\mathrm{\μ}}_1-{\mathrm{\μ}}_2)}^2}{s^2})& \mathrm{if}{({\mathrm{\μ}}_1-{\mathrm{\μ}}_2)}^2\≤s^2\\ 0& \mathrm{else}\end{array}\right.$

λ₂＝α

初步实现在结构信息较少的区域做平滑，而在结构信息较多的地方保留其结构信息。

第二次扩散系数的计算方法为：

 if(μ₁-μ₂)²＜ss_low

   λ₁＝α

else

 λ＝0

if(μ₁-μ₂)²＜ss_low or(μ₁-μ₂)²＞ss_high

λ₂＝α

else

λ₂＝0

ss_low，ss_high为预设的两个阈值，(μ₁-μ₂)²值低于ss_low认为是平坦区域，高于ss_high认为是强结构区域。第二次各向异性扩散滤波进一步平滑平坦区域，而在组织结构的边界上沿着边缘方向扩散，使组织结构信息增强。

两次扩散滤波之后，对结果进行锐化处理，以避免扩散过程可能带来的图像模糊，同时增强组织结构信息。

本发明在完成一次处理循环之后将继续在下一个尺度下对图像做类似的处理。不同尺度下判定区域是否平坦或是包含丰富结构信息，所用的阈值与所在的尺度相关，从而保证处理后的图像质量大幅提升。

Claims (1)

1.一种提升超声图像质量方法，用于超声成像系统对扫描图像显示数据的优化处理，其特征在于，步骤如下：

A.利用高斯-拉普拉斯金字塔分解算法分解显示的超声图像，将超声图像中的能量集中于数据量较少的高斯层；

B.对分解产生的高斯层图像数据进行各向异性扩散滤波；

C.对扩散结果进行锐化处理；

D.利用高斯-拉普拉斯金字塔重构算法将锐化处理过的高斯层图像数据与未处理的拉普拉斯层图像数据重构回去；

E.A、B、C、D四个步骤循环若干次；

在所述步骤B中，对分解产生的高斯层图像数据进行各向异性扩散滤波采用4个方向的各向异性扩散滤波，由以下几个步骤组成：

a.计算4个方向扩散系数，4个方向分别是直角坐标系中的-45°、0°、45°和90°方向；

b.在上述4个方向进行各项异性扩散滤波；

所述步骤a，由以下几个步骤组成：

(1)高斯平滑；(2)计算梯度；(3)计算结构张量；(4)计算扩散张量；(5)裂解为4个方向的扩散系数，具体为：

采用结构张量来表示局部的特征，某个像素的结构张量为

其中，I_x为x方向的梯度，I_y为y方向的梯度，结构张量分解为两个正交的向量，利用这两个向量的方向和特征值(u₁，u₂)来表征局部特性；

结构张量的计算要考虑到噪声的影响，计算梯度的时候，首先采用若干阶FIR滤波器平滑图像，平滑之后的图像噪声得到抑制，基于平滑之后的图像求出每个像素的x，y方向的梯度，获得结构张量的组成部分，即

并分别对这三个数据组成的阵列再次作平滑处理，进一步减小噪声对整个算法的影响，得到最终需要的结构张量

结构张量

的特征值和特征向量为：

${\mathrm{\μ}}_1=\frac{1}{2}(j_{11}+j_{22}+\sqrt{{(j_{11}-j_{22})}^2+{4j}_{12}^2}),$ ${\mathrm{\μ}}_2=\frac{1}{2}(j_{11}+j_{22}-\sqrt{{(j_{11}-j_{22})}^2+{4j}_{12}^2}),$

$\stackrel{\→}{{\mathrm{\ω}}_1}\left|\right|\left(\begin{array}{c}{2j}_{12}\\ j_{22}-j_{11}+\sqrt{{(j_{11}-j_{22})}^2+{4j}_{12}^2}\end{array}\right),$

$\stackrel{\→}{{\mathrm{\ω}}_2}\left|\right|\left(\begin{array}{c}{2j}_{12}\\ j_{22}-j_{11}-\sqrt{{(j_{11}-j_{22})}^2+{4j}_{12}^2}\end{array}\right),$

扩散张量的计算：由结构张量

的特征值，计算出扩散张量

的特征值λ₁和λ₂，再根据下式计算出扩散张量

${({\mathrm{\μ}}_1-{\mathrm{\μ}}_2)}^2={(j_{11}-j_{22})}^2+4*j_{12}^2,$

dd＝(λ₂-λ₁)(j₁₁-j₂₂)/(μ₁-μ₂)，

d₁₁＝0.5*[(λ₂+λ₁)+(λ₂-λ₁)(j₁₁-j₂₂)/(μ₁-μ₂)]，

d₁₂＝(λ₁-λ₂)j₁₂/(μ₁-μ₂)，

d₂₂＝0.5*[(λ₂+λ₁)-(λ₂-λ₁)(j₁₁-j₂₂)/(μ₁-μ₂)]，

再由扩散张量

转化为上述4个方向的扩散系数：

${\mathrm{\α}}_{-1}=\frac{\left|d_{12}\right|-d_{12}}{2},$

α₀＝d₁₁-|d₁₂|，

${\mathrm{\α}}_1=\frac{\left|d_{12}\right|+d_{12}}{2},$

α₂＝d₂₂-|d₁₂|；

其中，一次处理循环包含两次各向异性扩散滤波，第一次各向异性扩散滤波是初步实现在结构信息较少的区域作平滑，而在结构信息较多的地方保留其结构信息，第二次各向异性扩散滤波是进一步平滑平坦区域，而在组织结构的边界上沿着边缘方向扩散，使得组织结构信息增强；在第一次各向异性扩散滤波中，扩散张量

的特征值λ₁和λ₂的计算方法为：

${\mathrm{\λ}}_1=\left\{\begin{array}{cc}\mathrm{\α}\·(1-\frac{{({\mathrm{\μ}}_1-{\mathrm{\μ}}_2)}^2}{s^2})& \mathrm{if}{({\mathrm{\μ}}_1-{\mathrm{\μ}}_2)}^2\≤s^2\\ 0& \mathrm{else}\end{array}\right.,$

λ₂＝α，

在第二次各向异性扩散滤波中，扩散张量

的特征值λ₁和λ₂的计算方法为：

如果(μ₁-μ₂)²＜ss_low，

λ₁＝α，

否则，

λ₁＝0；

如果(μ₁-μ₂)²＜ss_low or(μ₁-μ₂)²＞ss_high，

λ₂＝α，

否则，

λ₂＝0；

其中，ss_low，ss_high为预设的两个阈值，(μ₁-μ₂)²值低于ss_low是平坦区域，高于ss_high是强结构区域，介于两个阈值之间的区域是结构组织，在所述平坦区域、强结构区域和结构组织三种区域采用不同的扩散系数计算算法。

Images