CN101872425B - 基于经验模态分解获取图像特征并测量相应物理参数方法 - Google Patents

Abstract

基于经验模态分解获取图像特征并测量相应物理参数方法，属于图像处理领域，本发明为解决现有通过图像分割提取图像特征的方法由于成像对象不均匀性以及图像噪声，无法分离特征层，导致提取图像特征及测量相关参数的精度低的问题，本发明：一、对图像进行自适应灰度拉伸，形成高对比度的图像，二、进行经验模态分解，获取IMF1，三、对IMF1进行梯度变换和分水岭分割，以得到封闭连续的特征轮廓曲线，四、两次扫描，获取上边界和下边界的采样点，五、用最小二乘法拟合，完成被测图像的特征层的提取，六、根据对拟合后的曲线进行横向扫描，均匀取多个采样点，获取上、下边界差值这一物理参数。

Description

基于经验模态分解获取图像特征并测量相应物理参数方法

技术领域

本发明涉及基于经验模态分解获取图像特征并测量相应物理参数的方法，属于图像处理领域。

背景技术

数据、公式、图表和图像等是描述事物或现象特性与特征的重要手段和方法，从图像中提取特征是模式识别的重要内容之一，从获取的特征中实现物理参数的自动测量可以达到从定性到定量分析的目的。从宏观的天体图像，到实用的医学图像，再到微观的显微图像，对图像进行分析、理解和特征提取，可以从中获取很多有用信息。以实用的医学图像为例，医学成像技术已经成为当前一种流行的诊断技术，比如CT图像、MRI图像以及超声图像等，它们通过不同的成像手段获得人体器官成像和软组织结构的详细资料，进而诊断多种疾病。如CT成像，利用X射线进行断层成像，得到人体切面的解剖结构；一般CT图像的分辨率高，诊断结果可信度高，但是X射线对人体有害，所以只有在必要的情况下才进行CT成像；与CT成像不同，超声成像具有无创、便携式、多功能等优点，而且不产生任何有害辐射。因此，不同的成像方式对于不同的诊断有着一定的适应性。

对于不同的成像模式来说，虽然成像手段不同，但是最终都是以图像为介质来描述事物或现象特性与特征的，医学图像即为病灶，从医学图像中实现物理参数的测量是医学诊断的重要依据，医生可以通过自动测量技术获得必要的物理参数信息，进而实现疾病诊断。比如超声图像中对于血管内膜厚度（IMT）的测量就是一种重要的医疗诊断技术，医生通过测量IMT的异常变化来判断血管疾病的发生；再如肝硬化CT图像中通过对肝实质进行形态分级，并结合CT肝脾容积测量，可以对肝硬化进行定量诊断。医学图像的定量测量，需要依据图像特征，通过对医学图像特征的分割及提取进行物理参数的自动测量。

如今国内外通过图像分割提取图像特征的方法主要包括基于模式识别技术、基于模型的方法、基于跟踪的方法、基于人工智能的方法、基于神经网络的方法、复杂管状结构检测方法等六大类。但很多图像由于其成像机制的限制，图像质量不高，特别是由于成像对象的不均匀性以及图像噪声带来的一般图像特征，如医学成像中器官或组织结构的不均匀性，一些微小的变化不能为图像所分辨，使得无法分离特征层，对图像的分割和处理更加困难，提取图像特征及测量相关参数的精度低。

发明内容                               

本发明目的是为了解决现有通过图像分割提取图像特征的方法由于成像对象不均匀性以及图像噪声，无法分离特征层，导致提取图像特征及测量相关参数的精度低的问题，提供了一种基于经验模态分解获取图像特征并测量相应物理参数方法。

本发明包括如下步骤：

步骤一、对图像进行自适应灰度拉伸，形成高对比度的图像，

步骤二、对步骤一形成的高对比度图像进行经验模态分解，获取一阶本征模态函数分量，

步骤三、对所述一阶本征模态函数分量进行梯度变换和分水岭分割，以得到封闭连续的特征轮廓曲线，所述特征轮廓曲线围成封闭特征区域，

步骤四、对所述封闭特征区域进行两次扫描，获取所述封闭特征区域上边界的采样点和下边界的采样点，

步骤五、用最小二乘法对所述封闭特征区域上边界的采样点和下边界的采样点分别进行拟合，以去除特征区域边界的错误的多余采样点，进而得到精确的特征区域的边缘曲线，完成被测图像的特征层的提取，

步骤六、对拟合后的边缘曲线进行横向扫描，均匀取多个采样点，计算每个相同横向坐标处的上、下边界采样点的纵向坐标之差，并计算多个所述纵向坐标之差的平均值，进而获取上、下边界差值这一物理参数。

本发明的优点：

1）本发明的目的在于提出一种基于图像经验模态分解（Image Empirical Mode Decomposition，IEMD）的目标分割特征提取及重要物理参数测量方法，它解决了当前有些图像特征提取方法无法分离特征层以改善提取精度的问题。

2）它对斑点噪声和对比度的不规则具有很强的抑制作用和可靠性，同时全自动运行，不需要人为参与。适用于不同尺寸、形状和斑块图像的特征提取与参数测量。

附图说明

图1为基于经验模态分解获取图像特征方法流程图；

图2为经验模态分解的流程图；

图3和图4为二次扫描确定特征区域；

图5至图7为三幅颈动脉图像；

图8为图5所述颈动脉图像的直方图统计结果；

图9为图6所述颈动脉图像的直方图统计结果；

图10为图7所述颈动脉图像的直方图统计结果； 

图11为图5所述颈动脉图像在[20 130]阈值时的图像拉伸结果；

图12为图6所述颈动脉图像在[30 210]阈值时的图像拉伸结果；

图13为图7所述颈动脉图像在[30 200]阈值时的图像拉伸结果；

图14为待处理的颈动脉内膜原图；

图15为经过灰度拉伸后的效果图；

图16为一阶IMF分量；

图17为二阶IMF分量；

图18为三创IMF分量；

图19至图22为对一阶IMF分量进行梯度变换和分水岭分割，得到的封闭连续的特征轮廓曲线的全过程；

图23为第一次粗定位扫描获得的图；

图24为第二次精定位扫描获得的图；

图25为颈动脉内膜上下边界的50对采样点； 

图26为拟合的曲线；

图27至图29为将一组具体的图像数据实验：作为原图像的图27按本发明方法提取，中间过程为图28，特征层结果为图29；

图30至图32为将另一组具体的图像数据实验：作为原图像的图30按本发明方法提取，中间过程为图31，特征层结果为图32。

具体实施方式

具体实施方式一：下面结合图1和图2说明本实施方式，

经验模态分解(Empirical Mode Decomposition，简称EMD)方法是由美国NASA的黄锷博士提出的一种信号分析方法。它依据数据自身的时间尺度特征来进行信号分解，无须预先设定任何基函数。这一点与建立在先验性的谐波基函数和小波基函数上的傅里叶分解与小波分解方法具有本质性的差别。正是由于这样的特点， EMD 方法在理论上可以应用于任何类型的信号的分解，因而在处理非平稳及非线性数据上，具有非常明显的优势。

利用信号内部时间尺度的变化做能量与频率的解析，将信号展开成数个本征模态函数（Intrinsic Mode Function，IMF），再利用希尔伯特变换（Hilbert Transform，HT）获得IMF的瞬时频率和振幅，上述过程总称为希尔伯特黄转换（Hilbert-Huang Transform，HHT）。

EMD是HHT算法的重要步骤，不同于使用固定形态窗口为分界基底函数的传统方法，EMD的基底函数是从信号中提取得到的，即使用IMF作基底。而IMF必须满足下列条件：

1）在整个函数中，极值点的数目与穿越零点的数目相等或者相差1；

2）在任何时刻，由局部极值包络线所定义的包络线局部均值为零。

其中，第一个条件与传统高斯平稳过程中窄频宽的要求类似。第二个条件是一个新的想法：将整体性要求改变为局部性要求，使得瞬时频率不会因为不对称波形的存在而导致不必要的晃动。依托这两个条件构建起来的EMD及HHT被认为是强有力地求解非线性、非平稳信号的自适应方法，是近年来对以傅立叶变换为基础的线性及稳态谱分析的重大突破，并得到了广泛的应用。

HHT从瞬时频率的定义和解法出发，定义了EMD方法和IMF的概念，通过EMD方法可以将任意信号分解为从高频到低频的IMF分量的信号叠加，对于图像信号相当于将信号分解成由不同频率信号组成的图像层。通过筛选得到所要的图像特征层为图像的特征提取以及提高提取精度提供了新思路。

本实施方式方法通过以下技术方案实现的：首先对图像进行自适应灰度拉伸后对经过图像拉伸的图像进行经验模态分解，得到扩展的图像特征函数，用本征模态函数表示，对通过筛选出的能表征图像特征的像素层进行梯度变换和分水岭分割，以得到特征轮廓曲线，通过对特征区域的二次扫描来进行区域特征定位，用最小二乘法对特征曲线进行拟合得到精确的目标轮廓曲线，具体方法包括以下步骤：

步骤一、对图像进行自适应灰度拉伸，形成高对比度的图像，

步骤二、对步骤一形成的高对比度图像进行经验模态分解，获取一阶本征模态函数分量，

步骤三、对所述一阶本征模态函数分量进行梯度变换和分水岭分割，以得到封闭连续的特征轮廓曲线，所述特征轮廓曲线围成封闭特征区域，

步骤四、对所述封闭特征区域进行两次扫描，获取所述封闭特征区域上边界的采样点和下边界的采样点，

步骤五、用最小二乘法对所述封闭特征区域上边界的采样点和下边界的采样点分别进行拟合，以去除特征区域边界的错误的多余采样点，进而得到精确的特征区域的边缘曲线，完成被测图像的特征层的提取，

步骤六、对拟合后的边缘曲线进行横向扫描，均匀取多个采样点，计算每个相同横向坐标处的上、下边界采样点的纵向坐标之差，并计算多个所述纵向坐标之差的平均值，进而获取上、下边界差值这一物理参数。

步骤一中对图像进行自适应灰度拉伸，分别搜索从低灰度级到高灰度级以及反向从高灰度级到低灰度级的灰度水平。将第一个最大像素灰度级记为背景灰度级，第二个最大灰度级作为前景灰度级。对经过图像拉伸的图像进行经验模态（IEMD）分解，得到扩展的图像特征函数，用本征模态函数表示；

本步骤首先采用灰度自适应拉伸法来提高关注区域的灰度动态范围，灰色拉伸又称为对比度扩展，它是图像增强过程中一种基本方法。给出的具体实施像中提出了直方图统计的方法。通常情况下，处于某一灰度级别的图像像素越多，该灰度级越重要，影响越大。因此，可以分别搜索从低灰度级到高灰度级以及反向从高灰度级到低灰度级的灰度水平。将第一个最大像素灰度级记为背景灰度级，第二个最大灰度级作为前景灰度级。

步骤二对步骤一形成的高对比度图像进行经验模态分解获取一阶本征模态函数分量的具体过程为：

设定高对比度图像输入信号为                                                

，

，

步骤21、IMF分解过程初始化：，且满足关系式

成立，其中为第

次分解后剩余的残差函数；

步骤22、筛选过程初始化，

，且满足关系式

成立，其中

为第次本征模态函数分解中经过第

次筛选后的剩余函数；

步骤23、根据筛选程序获取残差函数中经过第

次筛选后的剩余函数

，

所述残差函数是输入的待处理曲线

经过第

次本征模态函数分解的剩余的残差函数；

步骤24、采用标准偏差准则判断步骤23获得的剩余函数

是否满足本征模态函数的条件，即

是否小于阈值

，

；

判断结果为是，执行步骤25，判断结果为否，则

，然后执行步骤23，

步骤25、提取一阶本征模态函数分量IMF：；和获取高对比度图像输入信号

经过第1次本征模态函数分解的剩余的残差函数

。

其中，步骤23根据筛选程序获取输入信号

经过第

次本征模态函数分解的剩余趋势函数中经过第次筛选后的剩余函数

的过程为：

步骤31、利用三次样条函数获取输入信号

经过第次本征模态函数分解的剩余趋势函数中经过第

次筛选后的剩余函数

的上、下包络，

步骤32、计算所述剩余函数

上、下包络曲线在各个

的均值

，

步骤33、获取输入信号

经过第

次本征模态函数分解的剩余趋势函数中经过第

次筛选后的剩余函数。

步骤二提取的一阶本征模态函数分量IMF即是我们要用的能表征图像特征的像素层。

步骤三中的梯度变换采用Sobel算子。

在众多图像处理算法中，梯度变换算法可以减少斑点噪声的影响，产生更好的分割结果，常用于图像增强中。其中作为梯度变换算法中的一种，Sobel算子在边缘检测的过程中是一种常用的方法。它具有两种形式。我们一方面可以检测横向水平的边缘，另一方面可以用于检测垂直边缘。相对于其他一些算子，Sobel算子由于沿水平和垂直方向和噪声平滑效果更好，所以是一种边缘检测常用的方法。同时，和诸如拉普拉斯和对数算子相比Sobel算子具有更好的方向保持性。

同时，采用分水岭分割用于准确提取图像特征边缘。分水岭算法对微弱边缘反应灵敏，同时能够保证连续封闭边缘。此外，分水岭算法能获得了一块封闭的区域，这为图像的区域特征分析提供了便利。

步骤四中对所述封闭特征区域进行两次扫描，获取所述封闭特征区域上边界的采样点和下边界的采样点的过程为：

第一次扫描为粗定位，用于确定所述封闭特征区域的基本区域，所述基本区域的边界包括整个目标边界以及一些多余的组织边缘，

第二次扫描为精细定位，对所述基本区域进行重新扫描，用于定位目标边界的上限和下限，进而获得基本区域的上边界的采样点和下边界的采样点。

通过之前步骤的图像检测后，图像中目标组织边缘和其他多余的边缘被保留了下来。我们通过接下来的阶段来消除多余的边缘同时得到分割出的轮廓的确切位置。

虽然一般的扫描方式可以通过得到一系列的采样点来计算被分割目标的边缘，但如果图像中还有一些其他组织边缘，这样采样的最高点将远远偏离目标的上限或下限范围，导致不准确的测量结果。因此，一种改进扫描方式被提了出来，共包括两次扫描。第一次是用来确定基本区域，它包括了整个目标边界以及一些多余的组织边缘如图3所示，我们把它叫做粗定位。第二次对基本区域进行重新扫描，用于定位目标边界的上限和下限，如4所示。由于第二次扫描是在第一次扫描得到的基本区域上进行处理，得到了上边界和下边界的采样点，因此，称之为精细定位。经过两次扫描，目标边界得到完全的确定。

通过扫描，可以提取得到整个被分割目标边缘的多对样本采样点。但由于整个扫描过程可能会有一些错误采样点，有些不属于目标上边界和下边界错误样本点会被提取出来。因此，我们在这里通过对若干对采样点进行拟合从而建立目标的一个连续的边界。在我们的方法中采用多项式拟合法，通过最小二乘法对数据进行拟合。具体为步骤五所述。

步骤五中的用最小二乘法对所述封闭特征区域上边界的采样点和下边界的采样点分别进行拟合的过程相同，上边界的采样点拟合出特征区域上边界，下边界的采样点拟合出特征区域下边界，以下将上边界的采样点和下边界的采样点统称为采样点x，按如下公式获取特征区域的边缘曲线：

,

其中，

、

和

为系数。

假设给定数据m和矢量

，其中，

，是一个阶次不超过n的多项式函数。系数

和拟合多项式由下面的最小二乘法确定。

，                             

 ,                     

在我们的方法中，我们采样二次多项式进行拟合，多项式

可以写为

。

具体实施方式二，下面结合图1至图32进行说明，本实施方式为了评估利用经验模态分解获取图像特征并测量相应物理参数的方法,以颈动脉医学超声图像的特征提取和内膜厚度测量为例进行说明：

 执行步骤一：对图5至图7所示的三幅颈动脉血管内膜彩色超声图像进行图像拉伸处理，分别搜索从低灰度级到高灰度级以及反向从高灰度级到低灰度级的灰度水平。将第一个最大像素灰度级记为背景灰度级，第二个最大灰度级作为前景灰度级，图8至图10显示了三副动脉超声图像和相应的直方图统计结果。图11为图5所述颈动脉图像在[20 130]阈值时的图像拉伸结果；图12为图6所述颈动脉图像在[30 210]阈值时的图像拉伸结果；图13为图7所述颈动脉图像在[30 200]阈值时的图像拉伸结果。从图11至图13可以看出，颈动脉内膜中层超声图像区域已经得到了加强，背景灰度与前景灰度得到明显区分。

执行步骤二、对步骤一形成的高对比度图像进行经验模态分解，获取一阶本征模态函数分量。图14为待处理的颈动脉内膜原图；图15为经过灰度拉伸后的效果图；图16为经过图像拉伸的图像进行经验模态（IEMD）分解，获取一阶本征模态函数分量IMF1；如进行多次分解可获取多个IMF分量，图17为二阶IMF分量；图18为三创IMF分量，可以看出血管内膜区域得到了很好的加强和锐化，轮廓清晰可见。

执行步骤三：对通过筛选出的能表征图像特征的像素层（IMF1）进行梯度变换和分水岭分割，以得到颈动脉血管内膜的特征轮廓曲线。图19至图22所示为颈动脉内膜的检测分割全过程。

执行步骤四：通过对特征区域的两次扫描来进行区域特征定位，用最小二乘法对特征曲线进行拟合得到精确的颈动脉血管内膜曲线。

通过之前步骤二和三的图像检测后，内膜组织和其他多余的边缘被保留了下来。我们通过接下来的阶段来消除多余的边缘同时得到内膜的确切位置。图23为第一次粗定位扫描获得的图；显示了区域定位方式。扫描方向设置为从上到下。每次扫描过程中，沿扫描方向挑选前两个采样点。同时，以横坐标作为扫描位置灰度强度作为纵坐标绘制灰度强度曲线，如图24所示。

执行步骤五、用最小二乘法对所述封闭特征区域上边界的采样点和下边界的采样点分别进行拟合，以去除特征区域边界的错误的多余采样点，进而得到精确的特征区域的边缘曲线，完成被测图像的特征层的提取。

步骤六、对拟合后的曲线进行横向扫描，均匀取多个采样点，计算每个相同横向坐标处的上、下边界采样点的纵向坐标之差，并计算多个所述纵向坐标之差的平均值，进而获取上、下边界差值这一物理参数。

图25显示了颈动脉内膜上下边界的50对采样点；图26为拟合的曲线,从图25-26可以看出，整个拟合过程可为我们的应用提供准确动脉内膜边缘。

为了更加清晰的显示经过处理后的效果，我们将处理图像与原图像进行对比。同时进行了多幅图像的测试。图27至图29为将一组具体的图像数据实验：作为原图像的图27按本发明方法提取，中间过程为图28，特征层结果为图29；图30至图32为将另一组具体的图像数据实验：作为原图像的图30按本发明方法提取，中间过程为图31，特征层结果为图32。可以看出一副是有明确的边界的，另一幅边界不是很明显。无论是图像的质量如何，我们提出的方法可以产生很好的效果。

图像质量的优劣对超声图像的分割影响很大。这使得完成图像分割的任务变的很复杂。表1显示了整个实验数据的统计结果。在实验中，样品的采样点对数设置为50，即每重新扫描一张图片需扫描50对采样点。通过整个图像数据结果可以得出，无论图像的质量如何，我们可以得到准确的膜厚度测量结果。在最后一栏我们列出所用的时间不到0.6s，它能满足实时性的要求。

表 1 图像数据分析结果

 

Claims (7)

1.基于经验模态分解获取图像特征并测量相应物理参数方法，其特征在于，它包括如下步骤：

步骤一、对图像进行自适应灰度拉伸，形成高对比度的图像，

步骤二、对步骤一形成的高对比度图像进行经验模态分解，获取一阶本征模态函数分量，

步骤三、对所述一阶本征模态函数分量进行梯度变换和分水岭分割，以得到封闭连续的特征轮廓曲线，所述特征轮廓曲线围成封闭特征区域，

步骤四、对所述封闭特征区域进行两次扫描，获取所述封闭特征区域上边界的采样点和下边界的采样点，

步骤五、用最小二乘法对所述封闭特征区域上边界的采样点和下边界的采样点分别进行拟合，以去除特征区域边界的错误的多余采样点，进而得到精确的特征区域的边缘曲线，完成被测图像的特征层的提取，

步骤六、对拟合后的边缘曲线进行横向扫描，均匀取多个采样点，计算每个相同横向坐标处的上、下边界采样点的纵向坐标之差，并计算多个所述纵向坐标之差的平均值，进而获取上、下边界差值这一物理参数。

2.根据权利要求1所述的基于经验模态分解获取图像特征并测量相应物理参数方法，其特征在于，步骤二对步骤一形成的高对比度图像进行经验模态分解获取一阶本征模态函数分量的具体过程为：

设定高对比度图像输入信号为x(t)，t＝1，2，---，N，

步骤21、本征模态函数分解过程初始化：n＝1，且满足关系式r_n-1(t)＝x(t)成立，其中r_n-1(t)为第(n-1)次分解后剩余的残差函数；

步骤22、筛选过程初始化，k＝1，且满足关系式h_n(k-1)(t)＝r_n-1(t)成立，其中h_n(k-1)(t)为第n次本征模态函数分解中经过第(k-1)次筛选后的剩余函数；

步骤23、根据筛选程序获取残差函数中经过第k次筛选后的剩余函数h_nk(t)，

所述残差函数是输入的待处理曲线x(t)经过第n次本征模态函数分解的剩余的残差函数；

步骤24、采用标准偏差准则判断步骤23获得的剩余函数h_nk(t)是否满足本征模态函数的条件，即

是否小于阈值H_SD，0.2≤H_SD≤0.3；

判断结果为是，执行步骤25，判断结果为否，则k＝k+1，然后执行步骤23，

步骤25、提取一阶本征模态函数分量IMF1：c₁(t)＝h_1k(t)；和获取高对比度图像输入信号x(t)经过第1次本征模态函数分解的剩余的残差函数r₁(t)＝x(t)-c₁(t)。

3.根据权利要求2所述的基于经验模态分解获取图像特征并测量相应物理参数方法，其特征在于，步骤23根据筛选程序获取输入信号x(t)经过第n次本征模态函数分解的剩余的残差函数中经过第k次筛选后的剩余函数h_nk(t)的过程为：

步骤31、利用三次样条函数获取输入信号x(t)经过第n次本征模态函数分解的剩余趋势函数中经过第k-1次筛选后的剩余函数h_n(k-1)(t)的上、下包络，

步骤32、计算所述剩余函数h_n(k-1)(t)上、下包络曲线在各个t的均值

步骤33、获取输入信号x(t)经过第n次本征模态函数分解的剩余趋势函数中经过第k次筛选后的剩余函数 $h_{\mathrm{nk}}\left(t\right)=h_{n(k-1)}\left(t\right)-{\stackrel{\‾}{m}}_{n(k-1)}\left(t\right).$

4.根据权利要求2所述的基于经验模态分解获取图像特征并测量相应物理参数方法，其特征在于，步骤24的中H_SD＝0.25。

5.根据权利要求1所述的基于经验模态分解获取图像特征并测量相应物理参数方法，其特征在于，步骤三中的梯度变换采用Sobel算子。

6.根据权利要求1所述的基于经验模态分解获取图像特征并测量相应物理参数方法，其特征在于，步骤四中对所述封闭特征区域进行两次扫描，获取所述封闭特征区域上边界的采样点和下边界的采样点的过程为：

第一次扫描为粗定位，用于确定所述封闭特征区域的基本区域，所述基本区域的边界包括整个目标边界以及一些多余的组织边缘，

第二次扫描为精细定位，对所述基本区域进行重新扫描，用于定位目标边界的上限和下限，进而获得基本区域的上边界的采样点和下边界的采样点。

7.根据权利要求1或6所述的基于经验模态分解获取图像特征并测量相应物理参数方法，其特征在于，步骤五中的用最小二乘法对所述封闭特征区域上边界的采样点和下边界的采样点分别进行拟合的过程相同，上边界的采样点拟合出特征区域上边界，下边界的采样点拟合出特征区域下边界，以下将上边界的采样点和下边界的采样点统称为采样点x，按如下公式获取特征区域的边缘曲线：

p_n(x)＝a₀+a₁x+a₂x²

其中，a₀、a₁和a₂为系数。

Images