CN102930273B - 基于亮度自适应水平集的极光卵分割方法 - Google Patents

Abstract

本发明公开了一种基于亮度自适应水平集的极光卵分割方法，主要解决现有极光卵分割方法分割精度较低、鲁棒性较差和适用面较窄的缺陷。其实现步骤为：①采用形态学成分分析法对紫外极光图像进行预处理；②构建一个形态学显著性图作为极光卵的形状特征；③应用形态学显著性图的边缘曲线对水平集函数初始化；④计算亮度自适应的水平集演化速度和停止函数；⑤根据亮度自适应的水平集演化方程更新水平集函数；⑥提取更新完成后的零水平集曲线，并将其作为极光卵边界输出。本发明避免了现有分割方法出现的结果偏移和边缘泄漏现象，具有分割精度高和鲁棒性强的优点，可适用于各类紫外极光图像的分割。

Description

基于亮度自适应水平集的极光卵分割方法

技术领域

本发明属于图像处理技术领域，涉及紫外波段捕获图像的目标分割，可用于从复杂的紫外极光图像中将极光卵区域分割出来。

背景技术

极光是极地区域空间光学辐射中一种重要的自然发光现象。它是由太阳风和磁尾带电粒子体与大气中的分子、原子在近地空间相互作用激发而成的。它的发生和发展与太阳活动周期、地磁活动强度密切相关。紫外波段成像探测由于可以捕捉到任何其它波段都难以观测到的极其微弱的能量散射，进而对太阳进入地球极区的能量波动进行全天候、全区域的探测，从而成为极光探测的主要手段。紫外极光图像中的极光卵区域是极光粒子以磁极为中心沉降在地球南北极地区形成的椭圆带状区域，它的赤道向边界和极向边界是重要的地球物理参数，因此，有效提取紫外极光图像中的极光卵边界对于研究极光粒子沉降规律，进而理解太阳风-磁层能量耦合过程起着十分重要的作用。

由于紫外极光图像的强噪声、低对比度以及众多的混淆现象，如日晖、彗星、宇宙射线等，所以将极光卵区域从图像中分割出来具有很大的挑战性。目前，对紫外极光图像极光卵区域进行分割的方法主要可以分为两类：无形状知识的方法和有形状知识的方法。前者不需要任何形状先验知识，仅考虑图像自身的灰度和像素间的空间关系信息。例如，Hung等人在文献“C.Hung,G.Germany,K-means anditerative selection algorithms in image segmentation,in Proc.IEEE Southeastcon，2003,(Session 1:Software Development).”中提出的基于直方图的k-均值法(Histogram-basedK-means,HKM)，该方法对极光图像的灰度直方图进行k-均值聚类后，将高于灰度阈值的几个聚类作为极光卵区域；Li等人在文献“X.Li,R.Ramachandran，M.He,S.Movva,J.A.Rushing,S.J.Graves,W.Lyatsky,A.Tan，G.A.Germany,Comparingdifferent thresholding algorithms for segmenting auroras,in Proc.IEEE Conf.InformationTechnology:Coding and Computing,2004,pp.594-601.”中提出的自适应最小误差阈值法(Adaptive Minimum Error Thresholding,AMET)，该方法根据磁地方时(MagneticLocal Time,MLT)将紫外极光图像划分成若干个子区域，然后对每个子区域应用最小误差阈值法进行分割处理；Germany等人在文献“G.A.Germany,G.K.Parks,H.Ranganath,R.Elsen,P.G.Richards,W.Swift,J.P.Spann，M.Brittnacher,Analysis ofauroral morphology:substorm precursor and onset on January 10,1997,GeophysicalResearch Letters 25(15)(1998)3043-3046.”中提出的脉冲耦合神经网络法(Pulse-coupled Neural Network,PCNN)，该方法将图像的每个像素点看作一个分割神经元，利用相邻像素之间的灰度关系，采用神经网络法对极光卵区域进行分割。此类无形状知识的极光卵分割方法对于目标和背景灰度相差较大的极光图像分割效果较好，但是，在众多低对比度图像中会出现噪点大、误分率和漏分率较高等问题。而有形状知识的极光卵分割方法则结合了极光卵的形状特征先验知识，例如Cao等人在文献“C.Cao,T.S.Newman,G.Germany,A new shape-based auroral oval segmentationdrivenby LLS-RHT,Pattern Recognition 42(2009)607-618.”中提出的基于椭圆形状先验知识的线性最小二乘-随机霍夫变换法(Linear least-squares(LLS)-randomized Houghtransform(RHT),LLS-RHT)，该方法结合极光图像灰度信息和形状特征，用一个类似椭圆的图形来拟合极光卵区域，虽然能够对形态为完整型椭圆的极光卵进行较好的分割，但是，对于形态是缺口状椭圆的极光卵，分割结果会出现很大的偏差。

水平集方法是近年来图像分割领域研究的热点之一。其中，Li等人在文献“C.Li,C.Xu,C.Gui,M.D.Fox,Level set evolution without re-initialization:A newvariational formulation，in Proc.IEEE Conf.Comput.Vis.Pattern Recog.,2005,pp.430-436.”中提出的不需要初始化的水平集方法，由于在保证了高分割准确率的同时，较高地提升了分割速度而广泛应用在医学图像分割等领域。然而，Li等人的方法对曲线的初始化是随机的，而极光图像的背景和目标灰度相差较小，若初始曲线离目标边缘较远，分割结果就会出现很大的偏离。同时，该方法中的曲线演化速度和停止函数与曲线周围的灰度值没有相关性，因而在处理极光卵分割这类弱边缘问题时，即使选用了离目标边缘很近的初始化曲线，也会出现较大的边缘泄漏现象。

综上所述，现有的极光卵分割方法和目前常用的水平集图像分割方法对紫外极光图像都不能达到满意的分割效果。

发明内容

本发明的目的在于针对上述极光卵分割方法的不足，提出一种基于亮度自适应水平集的极光卵分割方法，以提高极光卵区域分割的精度和准确率，为空间物理研究提供尽可能精确的极光卵参数信息。

实现本发明的技术思路是：通过采用形态学成分分析法对紫外极光图像进行预处理，达到降噪的目的；通过构建一个形态学显著性图，作为极光卵的形状特征，并以形态学显著性图的边缘曲线作为水平集演化的初始曲线，采用亮度自适应的水平集演化方程对曲线进行演化，使曲线最终逼近目标的边缘，从而实现极光卵的分割，具体实现步骤包括如下：

(1)输入一幅紫外极光图像I，运用形态学成分分析的方法将紫外极光图像I分解为极光形态学子层I_A和噪声形态学子层I_N，去掉噪声形态学子层I_N，将极光形态学子层I_A输出，得到预处理后的图像I_P；

(2)根据磁地方时将预处理后的图像I_A分为24个子区域，对每个子区域进行基于最小平方误差准则的阈值选取，根据每个子区域不同的阈值对预处理后的图像I_A进行二值化操作，得到二值化后的图像I_B；

(3)将二值化后的图像I_B通过一组结构元素逐渐变大的复合形态学滤波器F，进行噪声点去除，得到滤波后的图像I_F；

(4)寻找滤波后的图像I_F中的最大连通域，找出最大连通域外边界曲线围成的图像区域，计算该图像区域内值为0的像素点个数N，设定经验阈值T_A＝1200，判断紫外极光图像的类型，若N大于等于T_A，则该图像属于完整型椭圆极光卵图像，若N小于T_A，则属于缺口状椭圆极光卵图像；

(5)根据紫外极光图像类型计算其形态学显著性图SMM，若图像属于完整型椭圆极光卵图像，则对滤波后的图像I_F进行椭圆拟合，得到椭圆拟合后的图像I_E，若图像属于缺口状椭圆极光卵图像，则对滤波后的图像I_F不做处理，即：

(6)将图像的形态学显著性图SMM的边界作为水平集演化的初始曲线，初始化水平集函数：

${\mathrm{\φ}}^0=\left\{\begin{array}{cc}d& (x,y)\∈\mathrm{SMM}-\∂\mathrm{SMM}\\ 0& (x,y)\∈\∂\mathrm{SMM}\\ -d& (x,y)\∈\mathrm{\Ω}-\mathrm{SMM}\end{array}\right.,$

其中，φ⁰为水平集函数的初始值，d为正的常数，Ω表示整个紫外极光图像区域，(x,y)表示紫外极光图像中第x行第y列的像素点，若该像素点属于形态学显著性图SMM以内的区域，则给φ⁰赋值为d，若该像素点在SMM的边界上，则给φ⁰赋值为0，若该像素点属于SMM以外的区域，则给φ⁰赋值为-d；

(7)计算像素点(x,y)及其八邻域的平均灰度值，将其归一化到[0,1]区间作为归一化的局部灰度平均值I_l(x,y)，计算亮度自适应的水平集演化速度v(x,y)和停止函数g(x,y)：

$v(x,y)=q\frac{1}{1+\exp (-p(I_l(x,y)-0.5))},$ $g(x,y)=\frac{1}{1+I_l(x,y){|\▿G_{\mathrm{\λ}}*I_A|}^m},$

其中，q为控制演化速度大小的常数，p为控制演化速度随I_l(x,y)非线性变化程度的常数，G_λ为标准差为λ的高斯核，*表示卷积运算，表示梯度运算，||表示取模运算，m为常数；

(8)设置水平集函数更新停止条件的阈值ε＝0.01，根据亮度自适应的水平集演化方程更新水平集函数：

(8a)输入t时刻的水平集函数φ^t，根据亮度自适应的水平集演化方程，计算水平集函数的变化值

$\frac{\∂\mathrm{\φ}}{\∂t}=u[\mathrm{\Δ\φ}-\mathrm{div}\left(\frac{\▿\mathrm{\φ}}{|\▿\mathrm{\φ}|}\right)]+\mathrm{\β\δ}\left(\mathrm{\φ}\right)\mathrm{div}\left(g\frac{\▿\mathrm{\φ}}{|\▿\mathrm{\φ}|}\right)+\mathrm{\δ}\left(\mathrm{\φ}\right)\mathrm{vg},$

其中，φ为水平集函数，表示水平集函数随时间的变化值，u和β为控制水平集演化曲线的形状如何变化的常数，Δ为拉普拉斯算子，div表示散度运算，δ(·)为单变量脉冲函数(Dirac函数)，v和g分别为亮度自适应的水平集演化速度和停止函数；

(8b)计算更新后的水平集函数φ^t+1：

${\mathrm{\φ}}^{t+1}={\mathrm{\φ}}^t+\mathrm{\Δt}\frac{\∂\mathrm{\φ}}{\∂t},$

其中，Δt为时间步长；

(8c)判断水平集函数更新前后的数值差异是否满足水平集函数更新的停止条件：

|φ^t+1-φ^t|≤ε，

若满足，则停止水平集函数的更新，输出水平集函数，若不满足，则继续进行水平集函数的更新；

(9)提取更新完成后的水平集函数中值为0的像素点，将其作为极光卵区域的边界，输出极光卵区域。

本发明与现有的极光卵分割方法和目前常用的水平集图像分割方法相比，具有以下优点：

1)本发明通过构建极光图像的形态学显著性图，获得了良好的极光卵形状特征，将该形态学显著性图的边缘作为水平集算法的初始演化曲线，能够实现极光卵区域的精确分割，避免了目前常用水平集分割方法出现的结果偏移现象；

2)本发明结合极光图像的局部灰度信息，设计了亮度自适应的水平集演化速度和停止函数，能够实现极光卵弱边缘的良好分割，避免了目前常用水平集分割方法出现的边缘泄漏现象；

3)本发明与现有的极光卵分割方法相比，不但能够实现目标与背景灰度相差较大的极光图像的精确分割，而且对目标和背景灰度相差不大的极光图像也能获得良好的分割效果，有效降低了误分率和漏分率，并且对完整型椭圆和缺口状椭圆的极光卵图像同时适用，扩大了极光卵分割方法的适用面。

本发明的技术过程和效果可结合以下附图详细说明。

附图说明

图1是本发明的流程图；

图2是本发明中基于磁地方时区域划分的示意图；

图3是现有完整型椭圆极光卵原图像；

图4是本发明与现有五种分割方法对图3的分割效果图；

图5是现有缺口状椭圆极光卵原图像；

图6是本发明与现有五种分割方法对图5的分割效果图。

具体实施方式

参照图1，本发明基于亮度自适应水平集的极光卵分割方法主要包括三个阶段：图像预处理、极光卵形状特征提取、亮度自适应水平集演化，具体步骤如下：

步骤1：对输入的紫外极光图像进行预处理：

1.1)构建紫外极光图像数据库，即将原始紫外极光频道定义CDF数据转化为标签图像文件TIFF数据，每个图像文件为大小是200×228的灰度图；

1.2)输入一幅紫外极光图像I，并对其进行预处理：

1.2.1)构建形态学成分分析方法的稀疏表示字典Φ_A和Φ_N；

1.2.2)设定形态学成分分析方法的条件为：

$I=I_A^{*}+I_N^{*},$

其中，为目标形态学子层，为背景形态学子层；

1.2.3)将目标形态学子层和背景形态学子层作为变量，代入形态学成分分析方法的目标函数：

$\underset{I_A^{*},I_N^{*}}{\min }{\left|\right|{\mathrm{\Φ}}_A^T{I}_A^{*}\left|\right|}_1+{\left|\right|{\mathrm{\Φ}}_N^T{I}_N^{*}\left|\right|}_1,$

即求解在稀疏表示字典Φ_A和Φ_N下，目标形态学子层和背景形态学子层的最优表达，其中，|| ||₁表示1范数运算；

1.2.4)输出目标形态学子层和背景形态学子层的最优表达，将目标形态学子层的最优表达作为极光形态学子层I_A，将背景形态学子层的最优表达作为噪声形态学子层I_N；

1.2.5)去掉噪声形态学子层I_N，将极光形态学子层I_A输出，得到预处理后的图像I_P。

步骤2：提取紫外极光图像中极光卵的形状特征：

2.1)对预处理后的图像I_P进行基于磁地方时自适应门限的二值化操作：

2.1.1)参照图2，根据磁地方时将预处理后的图像I_P分为24个子区域，每个子区域跨越2小时的磁地方时，相邻的两个子区域有1小时磁地方时的重叠；

2.1.2)选取24个子区域中的一块，假设它的目标和背景服从均值为μ_k，标准差为σ_k的正态分布，其中k＝1,2，当k＝1时表示目标，当k＝2时表示背景；

2.1.3)根据均值μ_k和标准差σ_k估计所选子区域的直方图分布H_e(i,t)：

$H_e(i,T)=\underset{k=1}{\overset{2}{\mathrm{\Σ}}}\frac{P_k}{\sqrt{2\mathrm{\π}}{\mathrm{\σ}}_k}\exp [-\frac{{(i-{\mathrm{\μ}}_k)}^2}{2{\mathrm{\σ}}_k^2}],$

其中，i表示该子区域中某一点像素的灰度值，T表示分割目标和背景的灰度阈值门限，P_k为目标和背景所占该子区域面积的比例；

2.1.4)计算所选子区域实际的直方图分布H(i)；

2.1.5)根据最小平方误差准则选取最优阈值门限T_m：

$T_m=\arg \underset{T}{\min }\underset{i=0}{\overset{255}{\mathrm{\Σ}}}{(H_e(i,T)-H\left(i\right))}^2,$

其中，(H_e(i,T)-H(i))为估计值对实际值的误差，T_m为能够使该误差的平方最小的阈值门限T的取值；

2.1.6)对所选子区域进行二值化操作，即将灰度值大于等于T_m的像素点当作目标并设定取值为1，将灰度值小于T_m的像素点当作背景并设定取值为0；

2.1.7)对24个子区域都进行二值化操作后，得到二值化后的图像I_B；

2.2)将二值化后的图像I_B进行多结构元素复合形态学滤波：

2.2.1)选择一组形状为圆形，大小逐步增大的结构元素{B}＝{B₁,B₂,...,B_n}，B_r为这组结构元素中的一个，r＝1,2,…n，n为这组结构元素的个数；

2.2.2)构建多结构元素复合形态学滤波器F：

F＝(((AοB_j)·B_j)οB_j+1)·B_j+1…  j＝1,...,n，

其中，A为将要通过滤波器F的图像，j为结构元素的序号，取值从1到n，ο表示形态学开操作，·表示形态学闭操作；

2.2.3)通过滤波器F对二值化后的图像I_B进行处理，得到输出的图像I_F：

I_F＝I_B*F+I_R*F，

其中，I_R＝(I_B-I_B*F)，其为二值化后的图像I_B与滤波器F卷积后的剩余部分；

2.3)寻找滤波后的图像I_F中的最大连通域，找出最大连通域外边界曲线围成的图像区域，计算该图像区域内值为0的像素点个数N，设定经验阈值T_A＝1200，判断紫外极光图像的类型，若N大于等于T_A，则该图像属于完整型椭圆极光卵图像，若N小于T_A，则属于缺口状椭圆极光卵图像；

2.4)若紫外极光图像属于完整型椭圆极光卵图像，则对滤波后的图像I_F进行椭圆拟合：

2.4.1)对滤波后的图像I_F随机采样l个点{p(x,y)l}，l＝1,2,…L,L≥5，其中(x,y)为采样点的坐标值；

2.4.2)用随机采样点{p(x,y)_l}拟合二次方程式：

A^*x²+B^*y²+C^*xy+D^*x+E^*y+F^*＝0，

其中，A^*为x的二次项系数，B^*为y的二次项系数，C^*为x和y的乘积项系数，D^*为x的一次项系数，E^*为y的一次项系数，F^*为常数项；

2.4.3)判断拟合的二次方程式是否满足如下条件：

$\left|\begin{array}{ccc}{A}^{*}& C^{*}/2& D^{*}/2\\ C^{*}/2& B^{*}& E^{*}/2\\ D^{*}/2& E^{*}/2& F^{*}\end{array}\right|\≠0,$ $\left|\begin{array}{cc}{A}^{*}& C^{*}/2\\ C^{*}/2& B^{*}\end{array}\right|>0,$

$\left|\begin{array}{ccc}{A}^{*}& C^{*}/2& D^{*}/2\\ C^{*}/2& B^{*}& E^{*}/2\\ D^{*}/2& E^{*}/2& F^{*}\end{array}\right|/(A^{*}+B^{*})<0,$

若满足，则进行下一步，若不满足，则重新进行步骤2.4.1)；

2.4.4)将二次方程式变形为椭圆的标准方程式：

$\frac{{[(x-x_c)\cos \mathrm{\&theta;}+(y-y_c)\sin \mathrm{\&theta;}]}^2}{a^2}+\frac{{[(y-y_c)\cos \mathrm{\&theta;}-(x-x_c)\sin \mathrm{\&theta;}]}^2}{b^2}=1,$

其中，(x_c,y_c)是椭圆中心，a和b为椭圆长短轴，θ为椭圆方向角；

2.4.5)根据椭圆中心(x_c,y_c)、椭圆长短轴a和b以及椭圆方向角θ，确定椭圆拟合后的图像I_E；

2.5)若紫外极光图像属于缺口状椭圆极光卵图像，则对滤波后的图像I_F不做处理；

2.6)将紫外极光图像的形态学显著性图SMM作为极光卵的形状特征：

步骤3：对紫外极光图像进行基于亮度自适应的水平集演化：

3.1)初始化水平集函数，即将紫外极光图像的形态学显著性图SMM的边界作为水平集演化的初始曲线，初始化水平集函数：

${\mathrm{\&phi;}}^0=\left\{\begin{array}{cc}d& (x,y)\&Element;\mathrm{SMM}-\&PartialD;\mathrm{SMM}\\ 0& (x,y)\&Element;\&PartialD;\mathrm{SMM}\\ -d& (x,y)\&Element;\mathrm{\&Omega;}-\mathrm{SMM}\end{array}\right.,$

其中，φ⁰为水平集函数的初始值，d为正的常数，Ω表示整个紫外极光图像区域，(x,y)表示紫外极光图像中第x行第y列的像素点，若该像素点属于形态学显著性图SMM以内的区域，则给φ⁰赋值为d，若该像素点在SMM的边界上，则给φ⁰赋值为0，若该像素点属于SMM以外的区域，则给φ⁰赋值为-d；

3.2)计算亮度自适应的水平集演化速度和停止函数：

3.2.1)计算像素点(x,y)及其八邻域的平均灰度值，将其归一化到[0,1]区间作为归一化的局部灰度平均值I_l(x,y)，

3.2.2)计算亮度自适应的水平集演化速度v(x,y)：

$v(x,y)=q\frac{1}{1+\exp (-p(I_l(x,y)-0.5))},$

其中，q为控制演化速度大小的常数，p为控制演化速度随I_l(x,y)非线性变化程度的常数；

3.2.3)计算亮度自适应的水平集停止函数g(x,y)：

$g(x,y)=\frac{1}{1+I_l(x,y){|\&dtri;G_{\mathrm{\&lambda;}}*I_A|}^m},$

其中，G_λ为标准差为λ的高斯核，*表示卷积运算，表示梯度运算，||表示取模运算，m为常数；

3.3)对水平集函数进行更新：

3.3.1)设定水平集函数更新停止条件的阈值ε＝0.01；

3.3.2)输入t时刻的水平集函数φ^t；

3.3.3)根据亮度自适应的水平集演化方程，计算水平集函数的变化值

$\frac{\&PartialD;\mathrm{\&phi;}}{\&PartialD;t}=u[\mathrm{\&Delta;\&phi;}-\mathrm{div}\left(\frac{\&dtri;\mathrm{\&phi;}}{|\&dtri;\mathrm{\&phi;}|}\right)]+\mathrm{\&beta;\&delta;}\left(\mathrm{\&phi;}\right)\mathrm{div}\left(g\frac{\&dtri;\mathrm{\&phi;}}{|\&dtri;\mathrm{\&phi;}|}\right)+\mathrm{\&delta;}\left(\mathrm{\&phi;}\right)\mathrm{vg},$

其中，φ为水平集函数，表示水平集函数随时间的变化值，u和β为控制水平集演化曲线形状变化的常数，Δ为拉普拉斯算子，div表示散度运算，δ(·)为单变量脉冲函数(Dirac函数)，v和g分别为亮度自适应的水平集演化速度和停止函数；

3.3.4)计算更新后的水平集函数φ^t+1：

${\mathrm{\&phi;}}^{t+1}={\mathrm{\&phi;}}^t+\mathrm{\&Delta;t}\frac{\&PartialD;\mathrm{\&phi;}}{\&PartialD;t},$

其中，Δt为时间步长；

3.3.5)判断水平集函数更新前后的数值差异是否满足水平集函数更新的停止条件：

|φ^t+1-φ^t|≤ε，

若满足，则停止水平集函数的更新，输出水平集函数，若不满足，则继续进行水平集函数的更新；

3.4)提取更新完成后的水平集函数中值为0的像素点，将这些像素点作为极光卵区域的边界并输出。

本发明的效果可通过以下仿真实验进一步说明：

实验1：极光卵分割效果比较

实验条件：实验用到的完整型椭圆极光卵原图像如图3所示，其中，图3(a)为原始极光卵图像，图3(b)为极光专家对图3(a)进行手工标定的原图像；实验用到的缺口状椭圆极光卵原图像如图5所示，其中，图5(a)为原始极光卵图像，图5(b)为极光专家对图5(a)进行手工标定的原图像。

实验内容：本实验将比较本发明和现有五种分割方法对紫外极光图像极光卵区域分割的效果，现有五种分割方法为背景技术中提及的HKM、AMET、PCNN、LLS-RHT以及Li等人的水平集方法。实验分成两部分：第一部分为完整型椭圆极光卵图像的分割效果比较，第二部分为缺口状椭圆极光卵图像的分割效果比较。

1.完整型椭圆极光卵图像的分割效果比较：图4是本发明与现有五种分割方法对图3的分割效果图，其中，图4(a)为HKM对图3的分割效果图，图4(b)为AMET对图3的分割效果图，图4(c)为PCNN对图3的分割效果图，图4(d)为LLS-RHT对图3的分割效果图，图4(e)为Li等人的水平集方法对图3的分割效果图，图4(f)为本发明对图3的分割效果图。

从图4中可以看出，HKM、AMET、PCNN对完整型椭圆极光卵图像的分割结果含有较多的噪点，并且误分率和漏分率较高；LLS-RHT的分割结果大体符合完整型椭圆极光卵的形态，但是却没有表现出极光卵内边界的凹凸起伏；Li等人的水平集方法出现了边缘泄漏，分割不完整的现象；而本发明的分割结果则比较精确。

2.缺口状椭圆极光卵图像的分割效果比较：图6是本发明与现有五种分割方法对图5的分割效果图，其中，图6(a)为HKM对图5的分割效果图，图6(b)为AMET对图5的分割效果图，图6(c)为PCNN对图5的分割效果图，图6(d)为LLS-RHT对图5的分割效果图，图6(e)为Li等人的水平集方法对图5的分割效果图，图6(f)为本发明对图5的分割效果图。

从图6中可以看出，HKM、AMET、PCNN对缺口状椭圆极光卵图像的分割结果存在与完整型椭圆极光卵图像相同的问题；LLS-RHT的分割结果出现了很大的偏移现象；Li等人的水平集方法分割结果大体良好，但仍旧没有本发明的结果精确。

实验2：极光卵分割误差比较

实验条件：本实验用到了30幅完整型椭圆极光卵原图像和30幅缺口状椭圆极光卵原图像。

实验内容：比较本发明与现有五种分割方法对极光卵原图像中极光卵区域进行分割的总分割误差，标准边界由极光专家手工标定。

对极光卵原图像中极光卵区域进行分割的总分割误差SE计算如下：

SE＝α^*·FS+β^*·MS+γ^*·BDE，

其中，FS为误分率，表示将背景像素点当作极光卵的误差，MS为漏分率，表示将极光卵像素点当作背景的误差，BDE为边界偏离误差，表示分割结果的极光卵边界与标准边界之间的偏移距离总和，α^*为误分率FS的权值比重，β^*为漏分率MS的权值比重，γ^*为边界偏离误差BDE的权值比重。

误分率FS、漏分率MS和边界偏离误差BDE根据下式进行计算：

$\mathrm{FS}=\frac{1}{A^{\&prime;}}\underset{(x,y)\&NotElement;R}{\mathrm{\&Sigma;}}r(x,y),$ $\mathrm{MS}=\frac{1}{A^{\&prime;}}\underset{(x,y)\&Element;R}{\mathrm{\&Sigma;}}[1-r(x,y)],$ $\mathrm{BDE}=\sqrt{E_{i^{*}=1}^{N^{\&prime;}}{d}^2\left(i^{*}\right)}/A^{\&prime;},$

其中，R为面积是A′的极光卵区域，该区域像素点取值为1，其余背景像素点取值为0，N′为分割结果图中极光卵边界像素点的个数，d²(i^*)表示分割结果的第i^*(i^*＝1,2,…N′)个边界像素点与标准边界像素点的距离，r(x,y)表示极光卵分割结果图中第x行第y列的像素点：

$r(x,y)=\left\{\begin{array}{cc}0,& (x,y)\&NotElement;R\\ 1,& (x,y)\&Element;R\end{array}\right..$

表6为本发明与现有五种分割方法对30幅完整型椭圆极光卵原图像的总分割误差比较。

表6.完整型椭圆极光卵图像分割结果的总分割误差比较

从表6可以看出，本发明的总分割误差的平均值最小，说明本发明对完整型椭圆极光卵图像的分割最精确，同时，本发明的总分割误差的均方差也最小，说明本发明对完整型椭圆极光卵图像的分割最鲁棒。

表7为本发明与现有五种分割方法对30幅缺口状椭圆极光卵原图像的总分割误差比较。

表7.缺口状椭圆极光卵图像分割结果的总分割误差比较

从表7可以看出，本发明的总分割误差的平均值最小，说明本发明对缺口状椭圆极光卵图像的分割最精确，同时，本发明的总分割误差的均方差也最小，说明本发明对缺口状椭圆极光卵图像的分割最鲁棒。

综上，本发明方法对紫外极光图像进行基于亮度自适应水平集的极光卵分割，提高了极光卵分割的准确性，增强了分割效果的鲁棒性，扩大了对不同形态极光卵图像的适用性。

Claims (4)

1.一种基于亮度自适应水平集的极光卵分割方法，包括如下步骤：

(1)输入一幅紫外极光图像I，运用形态学成分分析的方法将紫外极光图像I分解为极光形态学子层I_A和噪声形态学子层I_N，去掉噪声形态学子层I_N，将极光形态学子层I_A输出，得到预处理后的图像I_P；

(2)根据磁地方时将预处理后的图像I_P分为24个子区域，对每个子区域进行基于最小平方误差准则的阈值选取，根据每个子区域不同的阈值对预处理后的图像I_P进行二值化操作，得到二值化后的图像I_B；

(3)将二值化后的图像I_B通过一组结构元素逐渐变大的复合形态学滤波器F，进行噪声点去除，得到滤波后的图像I_F；

(4)寻找滤波后的图像I_F中的最大连通域，找出最大连通域外边界曲线围成的图像区域，计算该图像区域内值为0的像素点个数N，设定经验阈值T_A＝1200，判断紫外极光图像的类型，若N大于等于T_A，则该图像属于完整型椭圆极光卵图像，若N小于T_A，则属于缺口状椭圆极光卵图像；

(5)根据紫外极光图像类型计算其形态学显著性图SMM，若图像属于完整型椭圆极光卵图像，则对滤波后的图像I_F进行椭圆拟合，得到椭圆拟合后的图像I_E，若图像属于缺口状椭圆极光卵图像，则对滤波后的图像I_F不做处理，即：

(6)将紫外极光图像的形态学显著性图SMM的边界作为水平集演化的初始曲线，初始化水平集函数：

${\mathrm{\&phi;}}^0=\left\{\begin{array}{cc}d& (x,y)\&Element;\mathrm{SMM}-\&PartialD;\mathrm{SMM}\\ 0& (x,y)\&Element;\&PartialD;\mathrm{SMM}\\ -d& (x,y)\&Element;\mathrm{\&Omega;}-\mathrm{SMM}\end{array}\right.,$

其中，φ⁰为水平集函数的初始值，d为正的常数，Ω表示整个紫外极光图像区域，(x,y)表示紫外极光图像中第x行第y列的像素点，若该像素点属于形态学显著性图SMM以内的区域，则给φ⁰赋值为d，若该像素点在SMM的边界上，则给φ⁰赋值为0，若该像素点属于SMM以外的区域，则给φ⁰赋值为-d；

(7)计算像素点(x,y)及其八邻域的平均灰度值，将其归一化到[0,1]区间作为归一化的局部灰度平均值I_l(x,y)，计算亮度自适应的水平集演化速度ν(x,y)和停止函数g(x,y)：

$v(x,y)=q\frac{1}{1+\exp (-p(I_l(x,y)-0.5))},g(x,y)=\frac{1}{1+I_l(x,y){|{\&dtri;G}_{\mathrm{\&lambda;}}*I_A|}^m},$

其中，q为控制演化速度大小的常数，p为控制演化速度随I_l(x,y)非线性变化程度的常数，G_λ为标准差为λ的高斯核，*表示卷积运算，表示梯度运算，||表示取模运算，m为常数；

(8)设置水平集函数更新停止条件的阈值ε＝0.01，根据亮度自适应的水平集演化方程更新水平集函数：

(8a)输入t时刻的水平集函数φ^t，根据亮度自适应的水平集演化方程，计算水平集函数的变化值

$\frac{\&PartialD;\mathrm{\&phi;}}{\&PartialD;t}=u[\mathrm{\&Delta;\&phi;}-\mathrm{div}\left(\frac{\&dtri;\mathrm{\&phi;}}{|\&dtri;\mathrm{\&phi;}|}\right)]+\mathrm{\&beta;\&delta;}\left(\mathrm{\&phi;}\right)\mathrm{div}\left(g\frac{\&dtri;\mathrm{\&phi;}}{|\&dtri;\mathrm{\&phi;}|}\right)+\mathrm{\&delta;}\left(\mathrm{\&phi;}\right)\mathrm{vg},$

其中，φ为水平集函数，表示水平集函数随时间的变化值，u和β为控制水平集演化曲线的形状如何变化的常数，Δ为拉普拉斯算子，div表示散度运算，δ(·)为单变量脉冲Dirac函数，ν和g分别为亮度自适应的水平集演化速度和停止函数；

(8b)计算更新后的水平集函数φ^t+1：

${\mathrm{\&phi;}}^{t+1}={\mathrm{\&phi;}}^t+\mathrm{\&Delta;t}\frac{\&PartialD;\mathrm{\&phi;}}{\&PartialD;t},$

其中，Δt为时间步长；

(8c)判断水平集函数更新前后的数值差异是否满足水平集函数更新的停止条件：

|φ^t+1-φ^t|≤ε，

若满足，则停止水平集函数的更新，输出水平集函数，若不满足，则继续进行水平集函数的更新；

(9)提取更新完成后的水平集函数中值为0的像素点，将其作为极光卵区域的边界，输出极光卵区域。

2.根据权利要求1所述的基于亮度自适应水平集的极光卵分割方法，其中步骤(1)所述的运用形态学成分分析的方法将紫外极光图像I分解为极光形态学子层I_A和噪声形态学子层I_N，按如下公式进行：

$\underset{I_A^{*},I_N^{*}}{\min }{\left|\right|{\mathrm{\&Phi;}}_A^T{I}_A^{*}\left|\right|}_1+{\left|\right|{\mathrm{\&Phi;}}_N^T{I}_N^{*}\left|\right|}_1,$

$s.t.I=I_A^{*}+I_N^{*}$

其中，Φ_A为极光形态学子层的稀疏表示字典，Φ_N为噪声形态学子层的稀疏表示字典，为目标形态学子层，为背景形态学子层，||||₁表示1范数运算，形态学成分分析法就是在目标形态学子层和背景形态学子层之和为输入图像I的条件下，寻找稀疏表示字典Φ_A和Φ_N下，目标形态学子层和背景形态学子层的最优表达，并将目标形态学子层的最优表达作为极光形态学子层I_A，将背景形态学子层的最优表达作为噪声形态学子层I_N。

3.根据权利要求1所述的基于亮度自适应水平集的极光卵分割方法，其中步骤(2)所述的对每个子区域进行基于最小平方误差准则的阈值选取，按如下步骤进行：

(3a)选取24个子区域中的一块，假设它的目标和背景服从均值为μ_k，标准差为σ_k的正态分布，其中k＝1,2，当k＝1时表示目标，当k＝2时表示背景；

(3b)估计所选子区域的直方图分布H_e(i,T)：

$H_e(i,T)=\underset{k=1}{\overset{2}{\mathrm{\&Sigma;}}}\frac{P_k}{\sqrt{2\mathrm{\&pi;}}{\mathrm{\&sigma;}}_k}\exp [-\frac{{(i-{\mathrm{\&mu;}}_k)}^2}{{2\mathrm{\&sigma;}}_k^2}],$

其中，i表示该子区域中某一点像素的灰度值，T表示分割目标和背景的灰度阈值门限，P_k为目标和背景所占该子区域面积的比例；

(3c)计算所选子区域实际的直方图分布H(i)；

(3d)根据最小平方误差准则选取最优阈值门限T_m：

$T_m=\arg \underset{T}{\min }\underset{i=0}{\overset{255}{\mathrm{\&Sigma;}}}{(H_e(i,T)-H\left(i\right))}^2,$

其中，(H_e(i,T)-H(i))为估计值对实际值的误差，则T_m为能够使该误差的平方最小的阈值门限T的取值。

4.根据权利要求1所述的基于亮度自适应水平集的极光卵分割方法，其中步骤(3)所述的将二值化后的图像I_B通过一组结构元素逐渐变大的复合形态学滤波器F进行噪声点去除，按如下步骤进行：

(4a)选择一组形状为圆形，大小逐步增大的结构元素{B}＝{B₁,B₂,...,B_n}，B_r为这组结构元素中的一个，r＝1,2,…n，n为这组结构元素的个数；

(4b)构建复合形态学滤波器F：

F＝(((AоB_j)·B_j)оB_j+1)·B_j+1…j＝1,...,n，

其中，A为将要通过滤波器F的图像，j为结构元素的序号，取值从1到n，ο表示形态学开操作，·表示形态学闭操作；

(4c)通过滤波器F对二值化后的图像I_B进行处理，得到输出的图像I_F：

I_F＝I_B*F+I_R*F，

其中，I_R＝(I_B-I_B*F)，为二值化后的图像I_B与滤波器F卷积后的剩余部分。