基于图割优化的局部主动轮廓图像分割方法
【技术领域】
[0001] 本发明属于图像处理技术领域，涉及图像分割方法，可用于紫外极光图像的分割 以及灰度不均匀图像的分割。
【背景技术】
[0002] 极光是来自磁层的高能粒子沉降到高层大气并与中性成分碰撞激发的一种大气 发光现象，它主要出现在以地磁极为中心的环带状区域内，该区域又叫极光卵。极光卵的赤 道向边界和极向边界是重要的地球物理参数，其与太阳风、地磁活动有着密切的关系。而且 极光卵的边界随着地磁活动而变化，对其研究有助于进一步了解日地耦合过程，认知空间 气候变化规律。因此，有效确定极光卵边界是非常重要的。但由于极光卵就是紫外极光图 像的前景，而紫外极光图像的对比度较低，因此，用传统的紫外极光图像分割方法来确定极 光卵的边界都具有一定程度的局限性。
[0003] 在早期的紫外极光图像分割方法中主要有两种形式：一种是无形状先验的，如 2003年Hung等人提出的基于直方图的K均值算法HKM和2004年Li等人提出的自适应最 小误差阈值转换法AMET;另一种是加入形状先验的，如2009年Cao等人提出的基于随机霍 夫变换的线性最小二乘法LLSRHT。而前面两种无形状先验的方法都是基于像素亮度值的分 割方法，仅对于对比度高的紫外极光图像可以得到相对完整的极光卵。但是大多紫外极光 图像的对比度都比较低，所以，无形状先验的方法无法得到完整的极光卵。LLSRHT方法加入 了形状先验，虽然在紫外极光图像对比度较低的情况下也可以获得完整的极光卵，但该方 法获得的极光卵的边界过于平滑，与实际情况不符。
[0004] 近几年，主动轮廓方法非常流行，并且已被广泛应用到图像分割上。主动轮廓主要 分为基于边缘的主动轮廓和基于区域的主动轮廓。基于边缘的主动轮廓方法由于其对噪声 和初始轮廓的位置较敏感而逐渐被基于区域的主动轮廓方法代替。最典型的基于区域的 主动轮廓方法就是2001年ChanandVese提出的CV模型，该模型用梯度下降流来最小化 能量函数容易陷入局部最小。所以有人将CV模型和图割结合，用图割方法来优化CV模型 能量函数，提出了基于图割优化的主动轮廓方法ACBGC，该方法可以找到能量函数的全局最 小。但是该方法不能分割灰度不均匀的图像和对比度极低的紫外极光图像。
【发明内容】
[0005] 本发明的目的在于克服上述方法的不足，提出了一种基于图割优化的局部主动轮 廓图像分割方法，以实现对对比度极低的紫外极光图像和灰度不均匀图像的分割。
[0006] 实现本发明目的技术方案是：给定初始轮廓，并把图像轮廓内的部分当作目标，标 号为1，把图像轮廓外的部分当作背景，标号为〇 ;取图像的子图，在子图内用图割优化主动 轮廓能量函数，更新子图内所有像素的标号；通过迭代使整幅图像像素的标号得到更新，直 至曲线到达目标的边界时完成图像分割。其实现步骤如下：
[0007] (1)在输入图像I上定义初始轮廓曲线C，并对图像上的所有像素进行标号，即把 轮廓曲线内像素当作目标，标号为1，把轮廓曲线外像素当作背景，标号为O;
d(x，y) <a得到一个窄带区域N，a为一个常数，取值在1到2之间；
[0009] (3)以窄带区域内的第i个像素点为中心，r为半径做一个圆域，用这个圆域的外 接正方形为边界的区域作为图像I的一个子图I1，其中r的的值在8到40之间；
[0010] (4)在子图I1上分别计算轮廓内像素的灰度均值ci、轮廓外像素的灰度均值(：2和 能量函数E(C,C1,C2):
[0014] 其中，y为控制轮廓曲线光滑程度的参数，取y> 0,VXq为子图中相邻两像素P，q的标号，取值只能为〇或1 ;
[0015] (5)将步骤（3)得到的子图以央射为一个包含两种顶点和两种边的几何图，这两 种顶点分别是普通顶点和终端顶点，其中普通顶点对应于子图I1中的每个像素，终端顶点 对应于子图I1中前景部分和背景部分，即把前景作为第一终端顶点S，把背景作为第二终端 顶点T，两种边分别是n-links和t-links，其中n-links边对应于每两个相邻的普通顶点 的连线，且以步骤（4)得到的能量函数E((^c1,c2)为权值，t-links边对应于每个普通顶点 和2个终端顶点的连线，其权值按以下规则设置：
[0016] 如果第一终端顶点S与第二终端顶点T连接的普通顶点p满足（I(P)-C1)2 <(I(p)-c2)2,则将p与S连成的t-links边的权值设为9*109,这里"9*109"表示9和IO9 相乘；
[0017] 如果第一终端顶点S与第二终端顶点T连接的普通顶点p满足（I(p)-Cl)2> (I(p)-c2)2，则将p与T连成的t-links边的权值设为9*109;
[0018] (6)用最大流\最小割算法找到步骤（5)中几何图的最小割，以将几何图分成两 部分，一部分由第一终端顶点S以及S连接的部分普通顶点组成，这部分普通顶点被标号为 1，另一部分由第二终端顶点T以及T连接的部分普通顶点组成，这部分普通顶点被标号为 〇,以对子图I1的能量函数E(C，Cl，C2)进行更新，从而驱动子图I1内的轮廓曲线向前景的边 界演化；
[0019] (7)取窄带区域内的其它像素，重复步骤⑶到（6)中的操作，直至取完窄带区域 内所有像素，使图像I内的轮廓曲线C向前景的边界靠拢；
[0020] (8)重复步骤⑵到（7)中的操作，直至轮廓曲线C到达目标边界，这时图像在轮 廓曲线内的部分就是前景，在轮廓曲线外的部分就是背景，完成图像分割。
[0021] 本发明具有如下优点：
[0022] 1)本发明由于在分割紫外极光图像时，用基于随机霍夫变换的线性最小二乘方法 得到的初步分割结果的边界作为初始轮廓曲线，即加入了极光卵的形状先验信息，因此，能 够得到既完整又相对准确的极光卵边界，且由于初始轮廓曲线更接近极光卵的真实边缘， 使轮廓曲线演化到极光卵边界所用时间大大减少，从而加快了分割的速度；
[0023] 2)本发明基于图割优化的局部主动轮廓图像分割方法不仅克服了CV模型用梯度 下降流优化能量函数容易陷入局部最小的问题，同时，克服了CV模型和基于图割优化的主 动轮廓方法ACBGC不能分割灰度不均匀图像以及低对比度图像的问题。
【附图说明】
[0024] 图1是本发明的实现流程图；
[0025] 图2是用本发明对紫外极光图像进行分割的仿真结果图；
[0026] 图3是用本发明与现有方法对灰度不均匀图像进行分割的对比结果图；
[0027] 图4是本发明对同一图像选取不同r值的进行分割的对比结果图。
【具体实施方式】
[0028] 以下结合附图对本发明的内容和效果进行进一步描述。
[0029] 参照图1，本实例的具体实现步骤如下：
[0030] 步骤1 :在输入图像I上定义初始轮廓曲线C，并对图像上的所有像素进行标号。
[0031] la)在输入图像上定义初始轮廓曲线：
[0032] 若输入图像紫外极光图像，则先用已有的基于随机霍夫变换的线性最小二乘方法 得到初步的分割结果，然后提取初步分割结果的边界作为本发明的初始轮廓曲线；
[0033] 若输入图像为其他类型的图像，则用户自定义初始轮廓，但初始轮廓曲线必须要 与前景接触。
[0034] Ib)对图像上的像素进行标号时，把输入图像看作一个矩阵，其中这个矩阵的行数 就是图像的长，列数就是图像的宽，矩阵中的元素就是图像上对应位置的像素灰度值，再把 矩阵中对应于图像轮廓曲线内像素的元素换成1，即把图像轮廓曲线内的像素标号为1，把 矩阵中对应于图像轮廓曲线外像素的元素换成0,即把轮廓曲线外的像素标号为0。
[0035] 步骤2 :根据初始轮廓C构造一个窄带区域N。
[0036] 2a)构造符号距离函数p(.v,j')：
[0037] 以图像的左下角的像素为坐标原点，取图像上横坐标为x，纵坐标为y的像素 (x，y)，记（x，y)到轮廓曲线C的距离为符号距离函数夢其中若（x，y)在轮廓曲线 内，则记，若（X，y)在轮廓曲线外，则记，得到符号距离函数
[0038] 2b)令d(x，y) <a，即取图像上到轮廓曲线C的距离小于等于a的所有像素，用这 些像素构成一个窄带区域N，其中，a为一个常数，取值在1到2之间。
[0039] 步骤3 :以窄带区域内的像素i为中心，在输入图像上取子图L。
[0040] 3a)以窄带区域内的第i个像素点为中心，r为半径做一个圆域，
[0041] 其中r的大小根据图像灰度的不均匀性以及图像的背景复杂性选取，当图像灰度 较不均匀，背景也较复杂时，r的值应当选的小一些，当图像灰度比较均匀，背景也较单一 时，r的值应当选的大一些，但r的值在8-40之间，参照图4,r的大小对本发明的实验结果 有很大影响，所以要根据上述原则恰当的选取r的大小；
[0042] 3b)用这个圆域的外接正方形为边界的区域作为图像I的一个子图L。
[0043] 步骤4 :在子图上I1分别计算轮廓内像素的灰度均值C1、轮廓外像素的灰度均值C2 和能量函数E(C,C1,C2):
[0047] 其中，y为控制轮廓曲线光滑程度的参数，取y> 0,VXq为子图中相邻两像素 P、q的标号，取值为〇或1。
[0048] 步骤5 :将步骤⑶得到的子图I1映射为一个包含两种顶点和两种边的几何图。
[0049] 5a)确定几何图的顶点；
[0050] 这个几何图包括两种顶点，第一种顶点为普通顶点，即子图I1I的每个像素；第二 种顶点为终端顶点，即子图I1的前景和背景，其中，把前景作为第一终端顶点S，把背景作为 第二终端顶点T;
[0051] 5b)确定几何图的边；
[0052] 这个几何图包括两种边，第一种边为每两个普通顶点的连线，又叫n-links，第二 种边为每个普通顶点与2个终端顶点S和T的连线，又叫t-links;
[0053] 5c)确定几何图的边的权值；
[0054] 将n-1inks的权值设置为能量函数E(C，Cl，C2);
[0055] 根据以下规则设置t-1inks的权值则：
[0056] 如果第一终端顶点S与第二终端顶点T连接的普通顶点p满足（I(P)-C1)2 <(I(p) -C2)2，则将P与S连成的t-1inks边的权值设为9*109,这是因为如果p满足 (I(p)-C1) 2<(I(p)-C2)2,则P的灰度值与轮廓内像素的灰度均值C1更接近，就可以认为P 属于前景，即P和代表前景的第一终端顶点S有很高的相似度，则p与S连成的边的权值是 一个很大的数，这里设成9*109，"9*109"表示9和IO9相乘；
[0057] 如果第一终端顶点S与第二终端顶点T连接的普通顶点p满足（I(p)-Cl)2> (I(P)-C2)2，则将P与T连成的t-links边的权值设为9*109。这是因为如果p满足（I(P)-C1)2 >(I(p)-C2)2，则P的灰度值与轮廓外像素的灰度均值C2更接近，就