CN103065299B - 超声图像边缘提取方法和装置 - Google Patents

Abstract

本发明提供了一种超声图像边缘提取方法和装置。所述方法包括：对超声图像进行二维经验模态分解得到二维本征模态函数；解析所述二维本征模态函数得到所述超声图像的谱特征；对所述谱特征进行图像分割得到对应的图像边缘信息，并结合所述谱特征对应的图像边缘信息得到所述超声图像的边缘图像。所述系统包括：分解模块，用于对超声图像进行二维经验模态分解得到二维本征模态函数；特征解析模块，用于解析所述二维本征模态函数得到所述超声图像的谱特征；处理模块，用于对所述谱特征进行图像分割得到对应的图像边缘信息，并结合所述谱特征对应的图像边缘信息得到所述超声图像的边缘图像。采用本发明能提高准确性。

Description

超声图像边缘提取方法和装置

技术领域

本发明涉及图像处理技术，特别是涉及一种超声图像边缘提取方法和装置。

景技术

超声图像具有分辨率低、边缘模糊的特点，传统的边缘提取算法，例如，差分算子Roberts、Laplace、Sobel、Canny等边缘提取算法均存在着边缘像素宽、噪声干扰严重的缺点，有待进一步改善。

另一方面，近年来在新的理论基础上，如小波变换、分开、遗传、模糊数学和数学形态学等基础上所发展起来的新算法在超声图像中所提取得到的边缘模糊，受操场干扰大，也不适用于超声图像中边缘的提取。

发明内容

基于此，有必要针对超声图像的边缘提取受噪声干扰大，边缘提取不完整的缺陷的问题，提供一种能提高准确性的超声图像边缘提取方法。

另外，还有必要提供一种能提高准确性的超声图像边缘提取装置。

一种超声图像边缘提取方法，包括如下步骤：

对超声图像进行二维经验模态分解得到二维本征模态函数；

解析所述二维本征模态函数得到所述超声图像的谱特征；

对所述谱特征进行图像分割得到对应的图像边缘信息，并结合所述谱特征对应的图像边缘信息得到所述超声图像的边缘图像。

一种超声图像边缘提取装置，包括：

分解模块，用于对超声图像进行二维经验模态分解得到二维本征模态函数；

特征解析模块，用于解析所述二维本征模态函数得到所述超声图像的谱特征；

处理模块，用于对所述谱特征进行图像分割得到对应的图像边缘信息，并结合所述谱特征对应的图像边缘信息得到所述超声图像的边缘图像。

上述超声图像边缘提取方法和装置，进行二维经验模态分解得到二维本征模态函数，并解析二维本征模态函数得到超声图像的谱特征，进而通过二维经验模态分解和解析得到谱特征的作用下结合超声图像本身的性质得到连续清晰的边缘，由于谱特征能够准确地反映超声图像内部的客观信息，因此，实现了超声图像边缘的准确定位，进而提高了超声图像边缘提取的准确性，抗噪能力强。

附图说明

图1为一个实施例中超声图像边缘提取方法的流程图；

图2为图1中对超声图像进行二维经验模态分解得到二维本征模态函数的方法流程图；

图3为图1中解析二维本征模态函数得到超声图像的谱特征的方法流程图；

图4为图1中对谱特征进行图像分割得到对应的图像边缘信息，并结合谱特征对应的图像边缘信息得到超声图像的边缘图像的方法流程图；

图5为一个实施例中超声图像边缘提取装置的结构示意图；

图6为图5中分解模块的结构示意图；

图7为图5中特征解析单元的结构示意图；

图8为一个实施例中的超声骨骼图像；

图9为图8中的第一层二维本征模态函数示意图；

图10为通过图9得到的水平瞬时频率图像；

图11为通过图9得到的垂直瞬时频率图像；

图12为提取得到的边缘图像。

具体实施方式

如图1所示，在一个实施例中，一种超声图像边缘提取方法，包括：

步骤S110，对超声图像进行二维经验模态分解得到二维本征模态函数。

本实施例中，二维经验模态分解是将图像基于局部空间尺度自适应地分解为一系列空间尺度由小到大、频率由高到低的二维本征模态函数和一个残余项。

获取超声图像，对获取的超声图像进行二维经验模态分解以得到一系列的二维本征模态函数，换而言之，由一系列的二维本征模态函数分解的逆过程即可重构得到超声图像。

如图2所示，在一个实施例中，上述步骤S110的具体过程为：

步骤S111，在超声图像的趋势项中通过八邻域法搜索极值点得到超声图像中的极大值点和极小值点，该趋势项的初始值即为超声图像。

本实施例中，在进行二维经验模态分解的过程中，首先进行初始化，将二维本征模态函数指数的初始值设置为1，图像参数的初始值为超声图像，超声图像的趋势项为图像参数，即超声图像。

在超声图像的趋势项，即超声图像中搜索所有的极值点，其中，该极值点为极大值点和极小值点。

进一步的，通过八邻域法搜索极值点，在超声图像中，极大值点是灰度值比周围八个相邻像素点灰度值都高的点，极小值点是灰度值比周围八个相邻像素点灰度值都低的点。然而，对于超声图像中的边界点，将只会在一半的邻域内寻找极值点，对于位于四角处的边界点，将只考虑四分之一的邻域内。

如表1所示，超声图像中的像素点I(i，j)位于超声图像内部，需要与周围的八个像素点进行灰度值的比较；如表2所示，超声图像的像素点I(i，j)位于超声图像的边界，需要与周围的5个像素点进行灰度值的比较；如表3所示，超声图像的像素点I(i，j)是图像触点，位于超声图像边界上的四角处，需要与周围的3个像素点进行灰度值的比较，进而得到超声图像中的所有极值点。

I(i-1，j-1)	I(i-1，j)	I(i-1，j+1)
			I(i，j-1)	I(i，j)	I(i，j+1)2 -->
I(i+1，j-1)	I(i+1，j)	I(i+1，j+1)

表1

I(i-1，j)

I(i-1，j+1)

I(i，j)	I(i，j+1)
		I(i+1，j)	I(i+1，j+1)

表2

I(i，j)	I(i，j+1)
		I(i+1，j)	I(i+1，j+1)

表3

为实现上述三种情况下八邻域法的极值搜索，将对整个超声图像进行逐行扫描，并分上述三种情况进行灰度值的比较，进而确定超声图像中的极大值点和极小值点。

步骤S113，根据极大值点和极小值点得到超声图像的均值曲面。

本实施例中，在搜索得到超声图像中的极大值点和极小值点之后，将对散乱分布的极大值点和极小值点进行处理得到相应的曲面，即极大值曲面和极小值曲面，进而将极大值曲面和极小值曲面进行平均得到均值曲面。

进一步的，可通过基于三角剖分的立方插值法得到极大值曲面和极小值曲面。对于一组散乱分布的数据点，例如，极大值点或极小值点，如何将各点之间以三角形朴素连接，形成一张三角网格，并使网格质量较优，实质是以三角风格反映各个数据点与其邻近各点之间的拓扑连接关系。

在通过基于三角剖分的立方插值过程中，首先通过三角剖分得到三角形区域，并结合立方插值得到整个插值曲面，由于三角剖分能够较好地逼近不规则区域，因此使得通过三角剖分和立方插值所得到的极大值曲面和极小值曲面更为精确。

在分别进行立方插值得到极大值曲面和极小值曲面之后，将极大值曲面和极小值曲面平均得到均值曲面E_m，即E_m＝(E_u+E_l)/2，其中，E_u为极大值曲面，E_l为极小值曲面，实际上，均值曲面E_m是超声图像中的低频部分，表示超声图像中信号的整体趋势。

步骤S115，通过均值曲面对超声图像进行抽取得到高频部分。

本实施例中，由于均值曲面为超声图像中的低频部分，因此，通过均值曲面对超声图像进行抽取可直接得到超声图像中的高频部分Res，即Res＝Res-E_m，其中，高频部分代表了超声图像中信号的细节变化信息和噪声。

步骤S117，判断高频部分是否为二维本征模态函数，若是，则进入步骤S118，若否，则返回步骤S111。

本实施例中，由于上述步骤S111至步骤S115是一个不断循环的过程，若在当前的循环中高频部分被视为二维本征模态函数，则将二维本征模态函数置为高频部分，将二维本征模态函数指数j加1，将图像参数和二维本征模态函数之间的差值置为新的图像参数值，即BIMF(j)＝Res,j＝j+1,I＝I-Res，并通过步骤S118进一步判断当前所进行的二维经验模态分解是否达到了预设的分解层数，若是，则停止循环.

若判断到高频部分不可视为二维本征模态函数，则应当返回步骤S111得到新的高频部分，以再次进行二维本征模态函数的判断。

进一步的，可通过标准偏差（SD）的方式来检验高频部分是否被视为二维本征模态函数，即

$\mathrm{SD}=\frac{\underset{x=1}{\overset{M}{\mathrm{\Σ}}}\underset{y=1}{\overset{N}{\mathrm{\Σ}}}{|{\mathrm{Res}}_{j+1}(x,y)-{\mathrm{Res}}_j(x,y)|}^2}{\underset{x=1}{\overset{M}{\mathrm{\Σ}}}\underset{y=1}{\overset{N}{\mathrm{\Σ}}}{\left|{\mathrm{Res}}_j(x,y)\right|}^2}$

其中，(x,y)为超声图像中像素点在x和y方向上对应的数值，M为x方向上的最大值，N为y方向的最大值。

在计算得到的标准偏差SD处于预设的数值范围时，将认为高频部分为二维本征模态函数。在优选的实施例中，预设的数值范围是0.2~0.3。

步骤S118，将高频部分置为趋势项，并进一步判断二维经验模态分解是否达到预设的分解层数，若是，则结束，若否，则返回步骤S111。

本实施例中，二维经验模态分解达到预设的分解层数将作为上述步骤S111至117停止循环的条件，若当前所进行的二维经验模态分解还未达到预设的分解层数，则需要返回步骤S111中继续进行高频部分的运算。在优选的实施例中，考虑到二维经验模态分解的实时性和超声图像本身的特性，将认为超声图像的大部分边缘细节信息都包含在第一层二维本征模态函数中，因此，二维经验模态分解所预设的分解层数为两层，以得到第一层二维本征模态函数和趋势项。

步骤S130，解析二维本征模态函数得到超声图像的谱特征。

本实施例中，由于单纯通过经验模态分解实现超声图像的边缘提取也存在着噪声干扰大，边缘提取不完整的缺陷，因此需要结合超声图像的谱特征实现边缘的提取。

对二维本征模态函数进行信号解析以得到超声图像的谱特征，该谱特征反映了超声图像中的特征，包括超声图像的瞬时振幅、瞬时相位、二维瞬时频率等。在优选的实施例中，解析得到的谱特征为二维瞬时频率，以保证能够快速、精确地获知超声图像中的特征。

如图3所示，在一个实施例中，上述步骤S130的具体过程为：

步骤S131，对二维本征模态函数分别进行部分希尔伯特变换和整体希尔伯特变换得到部分变换信号和整体变换信号。

本实施例中，在分解得到了一系列二维本征模态函数之后，将对每一个二维本征模态函数分别进行部分希尔伯特变换和整体希尔伯特变换，以得到每一个二维本征模态函数所对应的部分变换信号和整体变换信号。

进一步的，对于超声图像中的二维本征模态函数进行的部分希尔伯特变换和整体希尔伯特变换可定义为如下三种形式：

定义1，设f为二维信号，即超声图像中的二维本征模态函数，它在x方向和y方向上的部分希尔伯特变换是：

$f_{\mathrm{Hx}}\left(x\right)=f\left(x\right)*\left\{\frac{\mathrm{\δ}\left(x\right)}{\mathrm{\πx}}\right\}$

$f_{\mathrm{HY}}\left(x\right)=f\left(x\right)*\left\{\frac{\mathrm{\δ}\left(y\right)}{\mathrm{\πy}}\right\}$

其中，f_Hx(x)为x方向上进行部分希尔伯特变换所得到的部分变换信号，f_Hy(x)为y方向上进行部分希尔伯特变换所得到的部分变换信号，δ(x)为水平方向冲激函数，δ(y)为垂直方向的冲激函数，π为圆周率，*表示二维卷积，频率域上定义为F_Hx(u)＝iF(u)sign(u)和F_Hy(u)＝iF(u)sign(v)， $\mathrm{sign}\left(u\right)=\left\{\begin{array}{c}1u>0\\ 0u=0\\ -1u<0\end{array}\right.,$ sign(v)与sign(u)类似，二维解析信号定义为e_⊥是单位向量，和e下次，即e_⊥·e＝0。

定义2，二维信号f在x方向和y方向上的整体希尔伯特变换为其中，*表示二维卷积，频率域上定义为F_H(u)＝-F(u)*sign(u)sign(v)，在此基础上，二维解析信号定义为 $f_A\left(x\right)=f\left(x\right)*[{\mathrm{\&delta;}}^2\left(x\right)+\frac{i}{{\mathrm{\&pi;}}^2\mathrm{xxy}}].$

定义3，单象限希尔伯特变换中解析信号的定义为 $f_A=f\left(x\right)*[\mathrm{\&delta;}\left(x\right)+\frac{i}{\mathrm{\&pi;x}}][\mathrm{\&delta;}\left(y\right)+\frac{i}{\mathrm{\&pi;y}}]=f\left(x\right)-{\mathrm{if}}_H\left(x\right)+i(f_{\mathrm{Hx}}\left(x\right)+f_{\mathrm{Hy}}\left(x\right)),$ 其中，f_H(x)是整体希尔伯特变换所得到的整体变换信号，f_Hx(x)和f_Hy(x)是部分希尔伯特变换所得到的部分变换信号。

然而，对于定义1而言，定义的解析信号只是一种简单的信号，即沿一个方向的信号，该方向可根据图像变化的方向进行选取，但是，本质上由于定义1所定义的解析信号是一维的，因而不同方向所得到的结果也不唯一。而定义2和定义3也并不完全符合解析信号在二维上应用的性质。

步骤S133，通过部分变换信号和整体变换信号进行四元数解析得到超声图像的谱特征。

本实施例中，由于定义1、定义2和定义3中的解析信号均不适用于二维上的应用，因此，需要通过四元数解析来弥补解析信号的不足。

一个四元数可写成超复数的形式，即q＝a+bi+cj+dk，其中，a、b、c和d是实数，i、j和k是复数单位，并满足运算规则：i²＝j²＝-1，ij＝ji＝k。

四元数由实数部分和虚数部分组成，其中，实数部分定义为a，虚数部分由三个成分组成，因此，虚数部分可以写成V(q)＝bi+cj+dk的形式，四元数可以写成q＝S(q)+V(q)的形式，若四元数的实数部分为零，则该四元数为纯四元数。

四元数的模和为： $\left|q\right|=\sqrt{a^2+b^2+c^2+d^2},$ 共轭定义为： $\stackrel{\&OverBar;}{q}=a-\mathrm{bi}-\mathrm{cj}-\mathrm{dk},$ 模是1为四元数称为单元四元数，通过和分别表示四元数的特征轴和特征角，进而定义四元数解析信号为：其中，x＝(x,y),n^T＝(i,j,k)^T和f_HT是向量形式，包括整体变换信号f_H和部分变换信号f_Hx(x)和f_Hy(x)，即

$f_{\mathrm{HT}}\left(\mathrm{\&chi;}\right)=\left\{\begin{array}{c}{f}_{\mathrm{Hx}}\left(\mathrm{\&chi;}\right)\\ f_{\mathrm{Hy}}\left(\mathrm{\&chi;}\right)\\ f_H\left(\mathrm{\&chi;}\right)\end{array}\right.$

在一个实施例中，上述步骤S133的具体过程为：通过极坐标的方式将二维本征模态函数、部分变换信号和整体变换信号定义为四元数解析信号，并通过四元数解析信号进行运算得到瞬时相位，根据瞬时相位计算得到二维瞬时频率。

本实施例中，在包含整体变换信号和部分变换信号的四元数解析信号中，通过极坐标的方式对四元数解析信号进行定义，进而对四元数解析信号进行运算得到瞬时相位，并根据瞬时相位计算得到二维瞬时频率。

对定义的四元数解析信号采用极坐标的方式表示，令其中， $A=\sqrt{f{\left(\mathrm{\&chi;}\right)}^2+f_{\mathrm{Hx}}{\left(\mathrm{\&chi;}\right)}^2+f_H{\left(\mathrm{\&chi;}\right)}^2+f_{\mathrm{Hy}}{\left(\mathrm{\&chi;}\right)}^2}$ 表示瞬时振幅，表示瞬时相位， $u=\frac{f_{\mathrm{Hx}}\left(\mathrm{\&chi;}\right)i+f_{\mathrm{Hy}}\left(\mathrm{\&chi;}\right)j+f_H\left(\mathrm{\&chi;}\right)k}{|f_{\mathrm{Hx}}\left(\mathrm{\&chi;}\right)i+f_{\mathrm{Hy}}\left(\mathrm{\&chi;}\right)j+f_H\left(\mathrm{\&chi;}\right)k|}=u_{1}i+u_{2}j+u_{3}k$ 表示向量部分的方向，u₁、u₂和u₃是其三个成分。

在此基础上，进一步计算得到垂直方向和水平方向上的瞬时频率：进而得到瞬时频率向量[diff1,diff2]。

通过上述极坐标的方式得到了瞬时相位和两个方向上的瞬时频率，这两个方向的瞬时频率构成了超声图像中的二维瞬时频率。

上述四元数解析将提取得到超声图像中二维瞬时频率所对应的图像，即水平瞬时频率图像和垂直瞬时频率图像。

步骤S150，对谱特征进行图像分割得到对应的图像边缘信息，并结合谱特征对应的图像边缘信息得到超声图像的边缘图像。

本实施例中，针对谱特征进行图像分割以得相应的图像边缘信息，在优选的实施例中，谱特征为二维瞬时频率，包括了水平瞬时频率图像和垂直瞬时频率图像。

具体的，针对谱特征进行二值化处理，以得到边缘信息，对于所得到的边缘信息，还可根据实际情况应用形态学细化边缘。

如图4所示，在一个实施例中，上述步骤S150的具体过程为：

步骤S151，根据分割阈值对谱特征进行二值化处理得到谱特征对应的图像边缘信息。

本实施例中，针对水平瞬时频率图像和垂直瞬时频率图像分别应用OTSU法（大津算法），自动搜索分割阈值，对水平瞬时频率图像和垂直瞬时频率图像中大于阈值的设置为逻辑“1”，小于阈值的设置为逻辑“0”，实现图像的二值化得以对应的图像边缘信息。

步骤S153，将谱特征对应的图像边缘信息相加得到超声图像的边缘图像。

本实施例中，将包含了“1”和“0”的图像边缘信息相加即可得到超声图像的边缘图像。

如图5所示，在一个实施例中，一种超声图像边缘提取装置，包括分解模块110、特征解析模块130和处理模块150。

分解模块110，用于对超声图像进行二维经验模态分解得到二维本征模态函数。

本实施例中，二维经验模态分解是将图像基于局部空间尺度自适应地分解为一系列空间尺度由小到大、频率由高到低的二维本征模态函数和一个残余项。

分解模块110获取超声图像，对获取的超声图像进行二维经验模态分解以得到一系列的二维本征模态函数，换而言之，由一系列的二维本征模态函数分解的逆过程即可重构得到超声图像。

如图6所示，在一个实施例中，上述分解模块110包括搜索单元111、极值处理单元113、抽取单元115和判断单元117。

搜索单元111，用于在超声图像的趋势项中通过八领域法搜索极值点得到超声图像中的极大值点和极小值点，该趋势项的初始值即为超声图像。

本实施例中，在进行二维经验模态分解的过程中，搜索单元111首先进行初始化，将二维本征模态函数指数的初始值设置为1，图像参数的初始值为超声图像，超声图像的趋势项为图像参数，即超声图像。

在超声图像的趋势项，即超声图像中搜索所有的极值点，其中，该极值点为极大值点和极小值点。

进一步的，搜索单元111通过八邻域法搜索极值点，在超声图像中，极大值点是灰度值比周围八个相邻像素点灰度值都高的点，极小值点是灰度值比周围八个相邻像素点灰度值都低的点。然而，对于超声图像中的边界点，搜索单元111将只会在一半的邻域内寻找极值点，对于位于四角处的边界点，搜索单元111将只考虑四分之一的邻域内。

超声图像中的像素点I(i，j)位于超声图像内部，需要搜索单元111与周围的八个像素点进行灰度值的比较；超声图像的像素点I(i，j)位于超声图像的边界，需要搜索单元111与周围的5个像素点进行灰度值的比较；超声图像的像素点I(i，j)是图像触点，位于超声图像边界上的四角处，需要搜索单元111与周围的3个像素点进行灰度值的比较，进而得到超声图像中的所有极值点。

为实现上述三种情况下八邻域法的极值搜索，搜索单元111将对整个超声图像进行逐行扫描，并分上述三种情况进行灰度值的比较，进而确定超声图像中的极大值点和极小值点。

极值处理单元113，用于根据极大值和极小值得到超声图像的均值曲面。

本实施例中，在搜索得到超声图像中的极大值点和极小值点之后，极值处理单元113将对散乱分布的极大值点和极小值点进行处理得到相应的曲面，即极大值曲面和极小值曲面，进而将极大值曲面和极小值曲面进行平均得到均值曲面。

进一步的，极值处理单元113可通过基于三角剖分的立方插值法得到极大值曲面和极小值曲面。对于一组散乱分布的数据点，例如，极大值点或极小值点，如何将各点之间以三角形朴素连接，形成一张三角网格，并使网格质量较优，实质是以三角风格反映各个数据点与其邻近各点之间的拓扑连接关系。

在通过基于三角剖分的立方插值过程中，极值处理单元113首先通过三角剖分得到三角形区域，并结合立方插值得到整个插值曲面，由于三角剖分能够较好地逼近不规则区域，因此使得通过三角剖分和立方插值所得到的极大值曲面和极小值曲面更为精确。

在分别进行立方插值得到极大值曲面和极小值曲面之后，极值处理单元113将极大值曲面和极小值曲面平均得到均值曲面E_m，即E_m＝(E_u+E_l)/2，其中，E_u为极大值曲面，E_l为极小值曲面，实际上，均值曲面E_m是超声图像中的低频部分，表示超声图像中信号的整体趋势。

抽取单元115，用于通过均值曲面对超声图像进行抽取得到高频部分。

本实施例中，由于均值曲面为超声图像中的低频部分，因此，抽取单元115通过均值曲面对超声图像进行抽取可直接得到超声图像中的高频部分Res，即Res＝Res-E_m，其中，高频部分代表了超声图像中信号的细节变化信息和噪声。

判断单元117，用于判断高频部分是否为二维本征模态函数，若是，则停止执行，若否，则将高频部分置为趋势项，并通知搜索单元111。

本实施例中，由于上述搜索单元111至抽取单元115之间的处理过程是一个不断循环的过程，若在当前的循环中高频部分被视为二维本征模态函数，则判断单元117将二维本征模态函数置为高频部分，将二维本征模态函数指数j加1，将图像参数和二维本征模态函数之间的差值置为新的图像参数值，即BIMF(j)＝Res,j＝j+1,I＝I-Res，并进一步判断当前所进行的二维经验模态分解是否达到了预设的分解层数，若是，则停止循环.

若判断到高频部分不可视为二维本征模态函数，则应当通知搜索单元111得到新的高频部分，以再次进行二维本征模态函数的判断。

进一步的，可通过标准偏差（SD）的方式来检验高频部分是否被视为二维本征模态函数，即

$\mathrm{SD}=\frac{\underset{x=1}{\overset{M}{\mathrm{\&Sigma;}}}\underset{y=1}{\overset{N}{\mathrm{\&Sigma;}}}{|{\mathrm{Res}}_{j+1}(x,y)-{\mathrm{Res}}_j(x,y)|}^2}{\underset{x=1}{\overset{M}{\mathrm{\&Sigma;}}}\underset{y=1}{\overset{N}{\mathrm{\&Sigma;}}}{\left|{\mathrm{Res}}_j(x,y)\right|}^2}$

其中，(x，y)为超声图像中像素点在x和y方向上对应的数值，M为x方向上的最大值，N为y方向的最大值。

在计算得到的标准偏差SD处于预设的数值范围时，将认为高频部分为二维本征模态函数。在优选的实施例中，预设的数值范围是0.2~0.3。

上述判断单元117还用于进一步判断二维经验模态分解是否达到预设的分解层数，若是，则停止执行，若否，则通知搜索单元111。

本实施例中，二维经验模态分解达到预设的分解层数将作为上述搜索单元111至判断单元117之间的循环过程停止的条件，若判断单元117判断到当前所进行的二维经验模态分解还未达到预设的分解层数，则需要通知搜索单元111中继续进行高频部分的运算。在优选的实施例中，考虑到二维经验模态分解的实时性和超声图像本身的特性，将认为超声图像的大部分边缘细节信息都包含在第一层二维本征模态函数中，因此，二维经验模态分解所预设的分解层数为两层，以得到第一层二维本征模态函数和趋势项。

特征解析模块130，用于解析二维本征模态函数得到超声图像的谱特征。

本实施例中，由于单纯通过经验模态分解实现超声图像的边缘提取也存在着噪声干扰大，边缘提取不完整的缺陷，因此需要结合超声图像的谱特征实现边缘的提取。

特征解析模块130对二维本征模态函数进行信号解析以得到超声图像的谱特征，该谱特征反映了超声图像中的特征，包括超声图像的瞬时振幅、瞬时相位、二维瞬时频率等。在优选的实施例中，特征解析模块130解析得到的谱特征为二维瞬时频率，以保证能够快速、精确地获知超声图像中的特征。

如图7所示，在一个实施例中，上述特征解析单元130包括变换单元131和四元数解析单元133。

变换单元131，用于对二维本征模态函数分别进行部分希尔伯特变换和整体希尔伯特变换得到部分变换信息和整体变换信息。

本实施例中，在分解得到了一系列二维本征模态函数之后，变换单元131将对每一个二维本征模态函数分别进行部分希尔伯特变换和整体希尔伯特变换，以得到每一个二维本征模态函数所对应的部分变换信号和整体变换信号。

进一步的，变换单元131对于超声图像中的二维本征模态函数进行的部分希尔伯特变换和整体希尔伯特变换可定义为如下三种形式：

定义1，设f为二维信号，即超声图像中的二维本征模态函数，它在x方向和y方向上的部分希尔伯特变换是：

$f_{\mathrm{Hx}}\left(x\right)=f\left(x\right)*\left\{\frac{\mathrm{\&delta;}\left(x\right)}{\mathrm{\&pi;x}}\right\}$

$f_{\mathrm{HY}}\left(x\right)=f\left(x\right)*\left\{\frac{\mathrm{\&delta;}\left(y\right)}{\mathrm{\&pi;y}}\right\}$

其中，f_Hx(x)为x方向上进行部分希尔伯特变换所得到的部分变换信号，f_Hy(x)为y方向上进行部分希尔伯特变换所得到的部分变换信号，δ(x)为水平方向的冲激函数，δ(y)为垂直方向的冲激函数，π为圆周率，*表示二维卷积，频率域上定义为F_Hx(u)＝iF(u)sign(u)和F_Hy(u)＝iF(u)sign(v)， $\mathrm{sign}\left(u\right)=\left\{\begin{array}{c}1u>0\\ 0u=0\\ -1u<0\end{array}\right.,\mathrm{sign}\left(v\right)$ 与sign(u)类似，二维解析信号定义为e_⊥是单位向量，和e下次，即e_⊥·e＝0。

定义2，二维信号f在x方向和y方向上的整体希尔伯特变换为其中，*表示二维卷积，频率域上定义为F_H(u)＝-F(u)*sign(u)sign(v)，在此基础上，二维解析信号定义为 $f_A\left(x\right)=f\left(x\right)*[{\mathrm{\&delta;}}^2\left(x\right)+\frac{i}{{\mathrm{\&pi;}}^2\mathrm{xxy}}].$

定义3，单象限希尔伯特变换中解析信号的定义为 $f_A=f\left(x\right)*[\mathrm{\&delta;}\left(x\right)+\frac{i}{\mathrm{\&pi;x}}][\mathrm{\&delta;}\left(y\right)+\frac{i}{\mathrm{\&pi;y}}]=f\left(x\right)-{\mathrm{if}}_H\left(x\right)+i(f_{\mathrm{Hx}}\left(x\right)+f_{\mathrm{Hy}}\left(x\right)),$ 其中，fH(x)是整体希尔伯特变换所得到的整体变换信号，f_Hx(x)和f_Hy(x)是部分希尔伯特变换所得到的部分变换信号。

然而，对于定义1而言，定义的解析信号只是一种简单的信号，即沿一个方向的信号，该方向可根据图像变化的方向进行选取，但是，本质上由于定义1所定义的解析信号是一维的，因而不同方向所得到的结果也不唯一。而定义2和定义3也并不完全符合解析信号在二维上应用的性质。

四元数解析单元133，用于通过部分变换信号和整体变换信号进行四元数解析得到超声图像的谱特征。

本实施例中，由于定义1、定义2和定义3中的解析信号均不适用于二维上的应用，因此，需要四元数解析单元133通过四元数解析来弥补解析信号的不足。

一个四元数可写成超复数的形式，即q＝a+bi+cj+dk，其中，a、b、c和d是实数，i、j和k是复数单位，并满足运算规则：i²＝j²＝-1，ij＝-ji＝k。

四元数由实数部分和虚数部分组成，其中，实数部分定义为a，虚数部分由三个成分组成，因此，虚数部分可以写成V(q)＝bi+cj+dk的形式，四元数可以写成q＝S(q)+V(q)的形式，若四元数的实数部分为零，则该四元数为纯四元数。

四元数的模和为： $\left|q\right|=\sqrt{a^2+b^2+c^2+d^2},$ 共轭定义为： $\stackrel{\&OverBar;}{q}=a-\mathrm{bi}-\mathrm{cj}-\mathrm{dk},$ 模是1为四元数称为单元四元数，通过和分别表示四元数的特征轴和特征角，进而定义四元数解析信号为：其中，χ＝(x,y),n^T＝(i,j,k)^T和f_HT是向量形式，包括整体变换信号f_H和部分变换信号f_Hx(χ)和f_Hy(χ)，即

$f_{\mathrm{HT}}\left(\mathrm{\&chi;}\right)=\left\{\begin{array}{c}{f}_{\mathrm{Hx}}\left(\mathrm{\&chi;}\right)\\ f_{\mathrm{Hy}}\left(\mathrm{\&chi;}\right)\\ f_H\left(\mathrm{\&chi;}\right)\end{array}\right.$

在一个实施例中，上述四元数解析单元133还用于通过极坐标方式将二维本征模态函数、部分变换信号和整体变换信号定义为四元数解析信号，并对四元数解析信号进行运算得到瞬时数位，根据瞬时相位计算得到二维瞬时频率。

本实施例中，在包含整体变换信号和部分变换信号的四元数解析信号中，四元数解析单元133通过极坐标的方式对四元数解析信号进行定义，进而对四元数解析信号进行运算得到瞬时相位，并根据瞬时相位计算得到二维瞬时频率。

对定义的四元数解析信号采用极坐标的方式表示，令其中， $A=\sqrt{f{\left(\mathrm{\&chi;}\right)}^2+f_{\mathrm{Hx}}{\left(\mathrm{\&chi;}\right)}^2+f_H{\left(\mathrm{\&chi;}\right)}^2+f_{\mathrm{Hy}}{\left(\mathrm{\&chi;}\right)}^2}$ 表示瞬时振幅，表示瞬时相位， $u=\frac{f_{\mathrm{Hx}}\left(\mathrm{\&chi;}\right)i+f_{\mathrm{Hy}}\left(\mathrm{\&chi;}\right)j+f_H\left(\mathrm{\&chi;}\right)k}{|f_{\mathrm{Hx}}\left(\mathrm{\&chi;}\right)i+f_{\mathrm{Hy}}\left(\mathrm{\&chi;}\right)j+f_H\left(\mathrm{\&chi;}\right)k|}=u_{1}i+u_{2}j+u_{3}k$ 表示向量部分的方向，u₁、u₂和u₃是其三个成分。

在此基础上，四元数解析单元133进一步计算得到垂直方向和水平方向上的瞬时频率：进而得到瞬时频率向量[diff1,diff2]。

四元数解析单元133通过上述极坐标的方式得到了瞬时相位和两个方向上的瞬时频率，这两个方向的瞬时频率构成了超声图像中的二维瞬时频率。

上述四元数解析将提取得到超声图像中二维瞬时频率所对应的图像，即水平瞬时频率图像和垂直瞬时频率图像。

处理模块150，用于对谱特征进行图像分割得到对应的图像边缘信息，并结合谱特征对应的图像边缘信息得到超声图像的边缘图像。

本实施例中，处理模块150针对谱特征进行图像分割以得相应的图像边缘信息，在优选的实施例中，谱特征为二维瞬时频率，包括了水平瞬时频率图像和垂直瞬时频率图像。

具体的，处理模块150针对谱特征进行二值化处理，以得到边缘信息，对于所得到的边缘信息，还可根据实际情况应用形态学细化边缘。

在另一个实施例中，上述处理模块150还用于根据分割阈值对谱特征进行二值化处理得到谱特征对应的图像边缘信息，将谱特征对应的图像边缘信息相加得到超声图像的边缘图像。

本实施例中，处理模块150针对水平瞬时频率图像和垂直瞬时频率图像分别应用OTSU法（大津算法），自动搜索分割阈值，对水平瞬时频率图像和垂直瞬时频率图像中大于阈值的设置为逻辑“1”，小于阈值的设置为逻辑“0”，实现图像的二值化得以对应的图像边缘信息。

处理模块150将包含了“1”和“0”的图像边缘信息相加即可得到超声图像的边缘图像。

以超声骨骼图像为例应用本发明进行边缘提取，其中，超声骨骼图像如图8所示，对超声骨骼图像进行二维经验模态分解所得到的第一层二维本征模态函数如图9所示，解析第一层二维本征模态函数所得到的超声图像的谱特征，即水平瞬时频率图像，如图10所示，和垂直瞬时频率图像，如图11所示，在对谱特征进行图像分割得到对应的图像边缘信息，并结合谱特征对应的图像边缘信息所得到的超声图像的边缘图像如图12所示。

上述超声图像边缘提取方法和装置，进行二维经验模态分解得到二维本征模态函数，并解析二维本征模态函数得到超声图像的谱特征，进而通过二维经验模态分解和解析得到谱特征的作用下结合超声图像本身的性质得到连续清晰的边缘，由于谱特征能够准确地反映超声图像内部的客观信息，因此，实现了超声图像边缘的准确定位，进而提高了超声图像边缘提取的准确性，抗噪能力强。

本领域普通技术人员可以理解实现上述实施例方法中的全部或部分流程，是可以通过计算机程序来指令相关的硬件来完成，所述的程序可存储于一计算机可读取存储介质中，该程序在执行时，可包括如上述各方法的实施例的流程。其中，所述的存储介质可为磁碟、光盘、只读存储记忆体（Read-OnlyMemory，ROM）或随机存储记忆体（RandomAccessMemory，RAM）等。

以上所述实施例仅表达了本发明的几种实施方式，其描述较为具体和详细，但并不能因此而理解为对本发明专利范围的限制。应当指出的是，对于本领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干变形和改进，这些都属于本发明的保护范围。因此，本发明专利的保护范围应以所附权利要求为准。

Claims (6)

1.一种超声图像边缘提取方法，包括如下步骤：

在所述超声图像的趋势项中通过八邻域法搜索极值点得到超声图像中的极大值点和极小值点，所述趋势项的初始值即为超声图像；

根据所述极大值点和极小值点通过三角剖分的立方插值法得到所述超声图像的极大值曲面和极小值曲面，将所述极大值曲面和极小值曲面进行平均得到均值曲面；

通过所述均值曲面对超声图像进行抽取得到高频部分；

通过标准差SD是否处于预设数值范围的方式判断所述高频部分是否为二维本征模态函数，若是，则结束，若否，则将所述高频部分置为趋势项，并返回所述在所述超声图像的趋势项中通过八邻域法搜索极值点得到超声图像中的极大值和极小值的步骤，所述标准差SD计算公式为：其中(x,y)为超声图像中像素点在x和y方向上对应的数值，M为x方向上的最大值，N为y方向的最大值，j为高频部分被视为二维本征模态函数BIMF(j)时对应的二维本征模态函数指数，其中j满足BIMF(j)＝Res，Res为所述高频部分，Res_j(x,y)为二维本征模态函数指数为j时，像素点(x,y)对应的高频部分，Res_j+1(x,y)为二维本征模态函数指数为j+1时，像素点(x,y)对应的二维本征模态函数值；

解析所述二维本征模态函数得到所述超声图像的谱特征；

根据分割阈值对所述谱特征进行二值化处理得到所述谱特征对应的图像边缘信息，针对水平瞬时频率图像和垂直瞬时频率图像分别应用大津算法自动搜索分割阈值，将所述谱特征对应的图像边缘信息相加得到所述超声图像的边缘图像。

2.根据权利要求1所述的超声图像边缘提取方法，其特征在于，所述解析所述二维本征模态函数得到所述超声图像的谱特征的步骤包括：

对所述二维本征模态函数分别进行部分希尔伯特变换和整体希尔伯特变换得到部分变换信号和整体变换信号；

通过所述部分变换信号和整体变换信号进行四元数解析得到所述超声图像的谱特征。

3.根据权利要求2所述的超声图像边缘提取方法，其特征在于，所述通过所述部分变换信号和整体变换信号进行四元数解析得到所述超声图像的谱特征的步骤为：

通过极坐标的方式将所述二维本征模态函数、部分变换信号和整体变换信号定义为四元数解析信号，并对所述四元数解析信号进行运算得到瞬时相位，根据所述瞬时相位计算得到二维瞬时频率。

4.一种超声图像边缘提取装置，其特征在于，包括：

分解模块，包括：

搜索单元，用于在所述超声图像的趋势项中通过八邻域法搜索极值点得到超声图像中的极大值点和极小值点，所述趋势项的初始值即为超声图像；

极值处理单元，用于根据所述极大值点和极小值点通过三角剖分的立方插值法得到所述超声图像的极大值曲面和极小值曲面，将所述极大值曲面和极小值曲面进行平均得到均值曲面；

抽取单元，用于通过所述均值曲面对超声图像进行抽取得到高频部分；

判断单元，用于通过标准差SD是否处于预设数值范围的方式判断所述高频部分是否为二维本征模态函数，若是，则结束，若否，则将所述高频部分置为趋势项，并通知所述搜索单元，所述标准差SD计算公式为：其中(x,y)为超声图像中像素点在x和y方向上对应的数值，M为x方向上的最大值，N为y方向的最大值，j为高频部分被视为二维本征模态函数BIMF(j)时对应的二维本征模态函数指数，其中j满足BIMF(j)＝Res，Res为所述高频部分，Res_j(x,y)为二维本征模态函数指数为j时，像素点(x,y)对应的高频部分，Res_j+1(x,y)为二维本征模态函数指数为j+1时，像素点(x,y)对应的二维本征模态函数值；

特征解析模块，用于解析所述二维本征模态函数得到所述超声图像的谱特征；

处理模块，用于根据分割阈值对谱特征进行二值化处理得到所述谱特征对应的图像边缘信息，针对水平瞬时频率图像和垂直瞬时频率图像分别应用大津算法自动搜索分割阈值，将谱特征对应的图像边缘信息相加得到超声图像的边缘图像。

5.根据权利要求4所述的超声图像边缘提取装置，其特征在于，所述特征解析模块包括：

变换单元，用于对所述二维本征模态函数分别进行部分希尔伯特变换和整体希尔伯特变换得到部分变换信号和整体变换信号；

四元数解析单元，用于通过所述部分变换信号和整体变换信号进行四元数解析得到所述超声图像的谱特征。

6.根据权利要求5所述的超声图像边缘提取装置，其特征在于，所述四元数解析单元还用于通过极坐标方式将所述二维本征模态函数、部分变换信号和整体变换信号定义为四元数解析信号，并对所述四元数解析信号进行运算得到瞬时相位，根据所述瞬时相位计算得到二维瞬时频率。