2. 河南师范大学 计算机与信息工程学院, 河南 新乡 453007
2. College of Computer and Information Engineering, Henan Normal University, Xinxiang Henan 453007, China
图像分割是图像处理的重要组成部分,它可以把图像分割成若干个感兴趣区域,为目标识别、编码等更高层图像处理服务。在图像分割算法中,主动轮廓方法能够完美地解决拓扑不变性的问题且不需人工干预,具有强大的竞争优势,它可分为基于边缘[1-2]和基于区域[3-4]两种模型。前者利用边缘信息执行分割,但对于弱边缘图像,分割曲线易穿越目标轮廓造成目标泄露;后者利用区域的统计信息,相对于前者鲁棒性更强,对于初始轮廓的依赖性更弱。主动轮廓方法中的无边缘主动轮廓模型(C-V模型)[1]假设图像分为前景和背景两个区域,其各自具有不同的强度均值且缺乏明显边缘信息。C-V模型是经典的全局区域型模型,但由于某些图像在获取过程中存在一些问题,导致了图像灰度缓慢变化(非同质或者不均匀),这种图像若利用分段常数的C-V模型效果较差。为此很多学者提出纠正灰度偏差的方法,最早可以追溯到1986年,但是利用主动轮廓纠正图像偏差直到近几年才被Li等[5-6]提出,此方法假设前景和背景都具有正态分布,并结合局部信息和全局信息自适应地对非均匀图像进行分割,效果较好。
若图像灰度统计分布未知或者前景、背景混叠较为严重,特别是目标非单一灰度均值的情况[7],利用训练样本的先验知识应是最有效的解决方法。从机器学习的角度而言,主动轮廓同样也可以分为有监督和无监督两种类型。利用自组织神经网络对训练样本进行学习从而逼近前景或者背景分布,再利用该分布结合主动轮廓引导曲线进化,这种有监督的主动轮廓模型可以弥补目标统计分布未知的缺陷,达到有效分割的目的。本文提出利用核自组织映射(Kernel Self-Organizing Map,KSOM)对前景和背景的先验数据进行学习从而匹配概率统计分布,然后利用该分布结合主动轮廓达到自动分割图像的目的。
1 无监督C-V模型 1.1 全局C-V模型[8]C-V模型的能量函数E(C)可以表示为:
$\begin{align} & E(C)=\mu \cdot \text{Length}(C)+\nu \cdot Area(in(C)) \\ & \text{ +}{{\lambda }_{1}}\int_{\text{For}(C)}{{{(I(x)-{{C}_{1}})}^{2}}dx} \\ & \text{ +}{{\lambda }_{2}}\int_{\text{Bac}(C)}{{{(I(x)-{{C}_{2}})}^{2}}dx} \\ \end{align}$ | (1) |
其中:C为前景目标轮廓;I(x)∈R是像素x点处的灰度值; μ≥0是正则参数控制曲线的光滑度;ν≥0是另外一个正则参数惩罚大的前景区域;For(C)和Bac(C)分别代表前景和背景区域,其均值分别是C1和C2,其能量项的控制参数分别是λ1和λ2。若用变分水平集的形式[9],式(1)变为:
$\frac{\partial \phi }{\partial t}={{\delta }_{\varepsilon }}(\phi )[\mu \kappa -\nu -{{\lambda }_{1}}{{(I-{{C}_{1}})}^{2}}+{{\lambda }_{2}}{{(I-{{C}_{2}})}^{2}}]$ | (2) |
其中:φ是水平集函数;δε(φ)是狄拉克函数;κ为曲线曲率。
${{\delta }_{\varepsilon }}(\phi )=\frac{1}{\pi }\frac{\varepsilon }{{{\varepsilon }^{2}}+{{\phi }^{2}}}$ | (3) |
参数ε可以控制δε(φ)的有效宽度。
水平集的初始化值φ0为:
${{\phi }_{0}}=\left\{ \begin{align} & \alpha ,\text{前景} \\ & -\alpha ,\text{背景} \\ \end{align} \right.\text{ }$ | (4) |
参数α必须满足α≈2ε。
若Ω是整个图像区域,曲线C可表示为:
$C=\{x\in \Omega :\phi (x)=0\}$ | (5) |
前景和背景可分别表示为:
$\left\{ \begin{align} & \text{For}(C)=\{x\in \Omega :\phi (x)>0\} \\ & Bac(C)=\{x\in \Omega :\phi (x)<0\} \\ \end{align} \right.$ | (6) |
全局C-V模型仅仅考虑全局信息,假设前景和背景的所有像素都分别服从各自的正态分布,因此在解决非同质或者由多个不同强度组成的目标分割方面效果欠佳。Liu等[10-11]提出基于局部区域的C-V模型,它用以下函数代替全局模型下的C1和C2:
${{C}_{1}}(x,C)=\frac{\int_{\text{For(C)}}{{{g}_{\sigma }}(x-y)I(y)dy}}{\int_{\text{For(C)}}{{{g}_{\sigma }}(x-y)dy}}$ | (7) |
${{C}_{2}}(x,C)=\frac{\int_{\text{Bac(C)}}{{{g}_{\sigma }}(x-y)I(y)dy}}{\int_{\text{Bac(C)}}{{{g}_{\sigma }}(x-y)dy}}$ | (8) |
其中:gσ是核宽为σ高斯核函数;C1(x,C)和C2(x,C)是像素x的σ邻域局部像素均值,把其代入式(1)和(2)即可以得到局部模型。
但是这种局部主动轮廓仍然局限于前景和背景的统计假设,并未探索利用已知先验样本对C-V模型进行监督。
2 核自组织映射自组织映射可以表征输入模式的拓扑映射结构,即用网络神经元的空间位置坐标表示输入模式包含的内在统计特征。由于自组织映射具有估计精度不足和最优目标函数缺乏两个局限,VAN Hulle[12]提出了基于核的自组织映射方式,它可以改善拓扑映射结构。
考虑由l个神经元组成的网格P,这些神经元由相应核集刻画:k(x,ωi,σi),目标函数是定义在第i个神经元输出之上:
${{y}_{i}}=k(x,{{\omega }_{i}},{{\sigma }_{i}}),\text{ }i=1,2,...,l$ | (9) |
获胜神经元被定义为:
$i(x)=\arg \underset{i}{\mathop{\max }}\,{{y}_{j}}(x),\text{ } j\in P$ | (10) |
和自组织映射(Self-Organizing Map,SOM)相同,邻域函数hj,i(x)为距获胜神经元i(x)的网格距离的单调减函数,ra为邻域函数hj,i(x)的范围。
${{h}_{j,i}}(x)=\exp (-\frac{{{\left\| {{x}_{j}}-{{\omega }_{i}} \right\|}^{2}}}{2{{r}_{a}}^{2}}),\text{ } j\in P$ | (11) |
前景和背景两种样本分别通过各自的KSOM进行训练,分别得到各自的分布特征,对于标量图像而言,特征映射图如图 1所示。
图 1中,标量图像的输入为一维像素灰度空间,输出为神经元的权值向量。把前景先验图像输入后经过训练可以得到前景图像神经元的值WF。把背景先验图像输入后可以得到背景图像神经元的值WB。利用此先验值可以在水平集部分引导曲线自动分割前景目标。测试图像在每一次迭代时都会得到一个前景轮廓C,此时前景中的每一个像素x的映射值为WF中获胜的神经元[13]:
${{\omega }_{F}}(x,C)=\{{{\omega }_{F,j}}|\arg \underset{j}{\mathop{\min }}\,\left\| {{\omega }_{F,j}}-\frac{\int_{\text{For(C)}}{{{g}_{\sigma }}(x-y)I(y)dy}}{\int_{\text{For(C)}}{{{g}_{\sigma }}(x-y)dy}} \right\|\}$ | (12) |
背景中的每一个像素映射值为WB中获胜的神经元:
${{\omega }_{B}}(x,C)=\{{{\omega }_{B,j}}|\arg \underset{j}{\mathop{\min }}\,\left\| {{\omega }_{B,j}}-\frac{\int_{\text{Bac(C)}}{{{g}_{\sigma }}(x-y)I(y)dy}}{\int_{\text{Bac(C)}}{{{g}_{\sigma }}(x-y)dy}} \right\|\}$ | (13) |
即使用最小距离法确定像素的映射值都为各自网络最匹配或者获胜的神经元值。
则此时能量函数变为:
$\begin{align} & E(C)=\mu \cdot \text{Length}(C)+\nu \cdot Area(in(C)) \\ & \text{ +}{{\lambda }_{1}}\int_{\text{For}(C)}{{{(I(x)-{{\omega }_{F}}(x,C))}^{2}}dx} \\ & \text{ +}{{\lambda }_{2}}\int_{\text{Bac}(C)}{{{(I(x)-{{\omega }_{B}}(x,C))}^{2}}dx} \\ \end{align}$ | (14) |
引入水平集函数φ后,式(14)变为:
$\begin{align} & E(C)=\mu \cdot \text{Length}(C)+\nu \cdot Area(in(C)) \\ & \text{ +}{{\lambda }_{1}}\int_{\Omega }{{{(I(x)-{{\omega }_{F}}(x,C))}^{2}}H(\phi (x))dx} \\ & \text{ +}{{\lambda }_{2}}\int_{\Omega }{{{(I(x)-{{\omega }_{B}}(x,C))}^{2}}(1-H(\phi (x)))dx} \\ \end{align}$ | (15) |
其中正则化Heaviside函数为H(φ):
${{H}_{\varepsilon }}(\phi )=\frac{1}{2}[1+\frac{2}{\pi }\arctan (\frac{\phi }{\varepsilon })]$ | (16) |
设
$\begin{align} & {{e}_{F}}(x,C)={{(I(x)-{{\omega }_{F}}(x,C))}^{2}} \\ & {{e}_{B}}(x,C)={{(I(x)-{{\omega }_{B}}(x,C))}^{2}} \\ \end{align}$ |
迭代方程式(2)变为:
$\frac{\partial \phi }{\partial t}={{\delta }_{\varepsilon }}(\phi )[\mu \kappa -\nu -{{\lambda }_{1}}{{e}_{F}}+{{\lambda }_{2}}{{e}_{B}}]$ | (17) |
KSOM训练好后,背景和前景的训练精度会有差别,如果不进行修正,那么测试图像的分割效果就会变差,假设背景和前景训练图像的平均误差分别为eBT和eFT。
${{e}_{BT}}=\sum\limits_{x\in {{U}_{B}}}{{{(I(x)-{{\omega }_{B}}(x))}^{2}}}/{{N}_{B}}$ | (18) |
${{e}_{FT}}=\sum\limits_{x\in {{U}_{F}}}{{{(I(x)-{{\omega }_{F}}(x))}^{2}}}/{{N}_{F}}$ | (19) |
其中:UB和UF分别是训练图像中的背景和前景像素集合;NB和NF分别为背景和前景的像素个数;ωB和ωF分别是上述集合像素x映射后的获胜神经元。
另外,eF和eB的差是驱动曲线进化的关键因素,它们前面的参数λ1和λ2,以前文章中均采用固定参数。本文在做了多次实验后,用启发式的方法,给出了λ1和λ2随着迭代次数n变化的公式。无论是前景或背景,它们的区域面积越小,误差就越小,仅仅用误差相减驱动曲线进化是不合理的,应该用误差除以各自面积,用单位面积误差相减进行驱动,所以为了达到前景和背景的平衡,上述两个参数应该与区域面积成反比,令λ2(n)始终为1,则λ1(n)为:
${{\lambda }_{1}}(n)=\frac{{{A}_{B}}(n)}{{{A}_{F}}(n)}$ | (20) |
其中:AF(n)和AB(n)分别为第n次迭代时前景和背景的面积。
经过以上分析,式(17)变为:
$\begin{align} & \phi (n+1)=\phi (n)+\Delta t{{\delta }_{\varepsilon }}(\phi )\cdot \\ & \text{ }[\mu \kappa -\nu -{{\lambda }_{1}}(n){{e}_{F}}/{{e}_{FT}}+{{e}_{B}}/{{e}_{BT}}] \\ \end{align}$ | (21) |
首先对前景区域进行训练,步骤如下:
1) 初始化。首先对初始权值ωi(0)和核宽σi(0)(i=1,2,…,lF)选择随机值,其中lF为前景区域网络结构中神经元个数,初始化学习率参数和最大迭代步数N。
2) 取样。从前景区域中按一定概率取出一个像素值输入网络。
3) 相似性匹配。在算法的第n次迭代,按照式(10)确定获胜神经元i(x)。
4) 自适应。用式(22)~(23)调整权值向量和每个核的宽度:
${{\omega }_{j}}(n+1)={{\omega }_{j}}(n)+\frac{{{\eta }_{w}}{{h}_{j,i}}(x)}{\sigma _{j}^{2}}(x(n)-{{\omega }_{j}}(n))$ | (22) |
${{\sigma }_{j}}(n+1)={{\sigma }_{j}}(n)+\frac{{{\eta }_{\sigma }}{{h}_{j,i}}(x)}{{{\sigma }_{j}}(n)}[\frac{{{\left\| x(n)-{{\omega }_{j}}(n) \right\|}^{2}}}{m\sigma _{j}^{2}(n)}-1]$ | (23) |
其中:ηw和ησ为学习率参数;m=1为输入样本像素的维数。
5) 若达到最大迭代次数结束迭代过程。
背景区域的训练过程和前景区域相同。通过以上过程即可得到前景图像神经元的值ωF和背景图像神经元的值ωB,同时对以上两个网络评估训练精度,得到误差eBT和eFT。
3.2.2 测试过程步骤如下:
1) 初始化水平集函数。在测试图像域中选取一块区域,区域内部φ>0,边缘φ=0,其他位置φ<0,初始化最大迭代次数M。
2) 根据式(12)和(13)计算测试图像的映射输出,并通过式(20)计算λ(n),利用式(21)更新水平集函数φ。
3) 利用φ(n+1)=φ(n)+Δt·▽2φ(n)对水平集函数正则化。
4) 若曲线不再变化或者达到最大迭代次数结束进化。
4 实验结果和分析 4.1 仿真结果和定性分析为验证本文算法核自组织映射主动轮廓(Kernel Self-Organizing map Active Contour,KSOAC)的适用性,在处理器为Intel Core i7-3770 CPU PC中,利用Matlab7.11.0 对多幅强度不均匀图像进行仿真,并和C-V[1]、LGD(Local Gaussian Distribution)[9]、SOM主动轮廓(Self-Organizing map Active Contour,SOAC)算法进行实验对比和理论分析。文中采用小汽车图像、飞机图像和心脏图像进行说明,其原始图像如图 2所示,三幅图均为强度不均匀图像。飞机图像的背景像素强度不均匀,四个尖端区域边缘模糊,小汽车图像整体强度不均匀,并且前景区域有多个目标强度,心脏图像前景区域和背景区域灰度差别较小。三幅图像均采用图 1网络结构,其中飞机图和小汽车图均为3个神经元,心脏图为5个神经元,KSOM的训练参数如表 1所示。
KSOAC和SOAC的训练误差以及分割目标所需时间对比如表 2,训练后神经元权值如表 3所示。从表 2可以看出对于训练样本无论背景还是前景KSOM相对于SOM的误差更小。从表 3可看出飞机和小汽车背景和前景灰度值差别较为分明,而心脏图的背景神经元最大权值和前景神经元最小权值较为接近,因此最可能造成错分;同时心脏背景神经元的权值112和114,5和7差别较小,可以考虑优化神经元个数从而降低计算复杂度。表 4为本文算法水平集分割时所需的参数。
图 3为三幅图像前景和背景的训练区域,其中细线矩形区域均为背景,粗线区域均为前景;图 4为水平集初始轮廓;图 5为三幅图像不同算法的对比。从飞机图像实验结果可看出:LGD的分割效果较差,C-V在边缘处分割效果较差,本文算法分割的轮廓相对于SOAC较为完整,特别是在4个尖端部分,这是因为本文算法相对于SOAC算法训练误差更小,并且在能量函数中考虑了训练误差的影响。从小汽车图像实验结果可看出:C-V分割出了较多的虚假目标,效果较差;LGD在汽车右轮子处和车顶部分不如本文算法,但是在车底部的分割效果较好;而SOAC最终分割效果最差,陷入局部极值停止迭代。从心脏图像实验结果可看出:C-V、LGD和SOAC均未达到分割前景和背景的目的,而本文算法右端目标分割不够完整,这是因为此部分前景灰度值和背景神经元的最大权值比较接近,造成了错分。图 6为本文算法心脏图分割的中间过程,展现了随着迭代次数的增加,分割效果的变化以及能量项的控制参数数值变化,其中iter为迭代次数。
图 7为采用本文算法,在小汽车分割过程中λ1(n)随着迭代次数变化的曲线。从图 7可以看出初始时刻背景面积与前景面积比大约为19,经过约20次迭代后,λ1(n)趋于稳定,即此时背景和前景区域的比达到定值不再变化。
利用Matlab 7.11.0 软件,在程序代码中设置tic和toc得到程序消耗的时间,本文对比了SOAC和KSOAC算法,结果如表 2最后一列。可以看出,本文算法相对于SOAC要慢一点。这主要是因为KSOM需要调整核宽,且水平集迭代时,每次迭代都需要计算能量项的控制参数λ1(n)。
4.2 定量分析为评价KSOAC算法的分割精度,本文采用三个变量准确率P、查全率R和F参数定量分析,定义如下:
$\left\{ \begin{align} & P=\frac{TP}{TP+FP} \\ & R=\frac{TP}{TP+FN} \\ & F=\frac{2PR}{P+R} \\ \end{align} \right.$ | (24) |
其中:TP为被判定为前景样本,事实上是前景样本;FP为被判定前景样本,事实上是背景样本;FN为被判定为背景样本,事实上是前景样本;P为判定为前景样本中真正前景样本的比重,与过分割程度关联较大;R为被判定为前景样本占总前景样本的比重,与漏分割程度关联较大;F是总体精确度评价指标,值越大,分割效果越好。真实的目标轮廓如图 8所示,基于不同算法的分割效果度量参数如表 5所示。对于LGD算法,小汽车图的F值最大,分割效果最好,但飞机图的F值最小,分割效果最差,这可能是因为参数调节的问题,特别是能量项控制参数没有选好,因此LGD算法对于能量项控制参数的选取比较敏感。SOAC算法在小汽车图上的分割效果最差,同样是因为参数调节的问题,同时它欠缺考虑训练误差的影响。C-V算法的F值较稳定但同样存在过分割的情况,在小汽车图和心脏图上面有较多的虚假目标。总体来说KSOAC的准确率、查全率和F值最稳定,都超过了0.9,且平均值最大,所以KSOAC的稳健性更好。但是小汽车图、飞机图和心脏图三幅图比较,KSOAC在心脏图上的分割效果较差,有漏检现象,这主要是因为前景和背景灰度值过于接近。
本文研究利用核自组织神经映射和主动轮廓对目标进行分割,由于核自组织映射相对于自组织映射改善了拓扑映射结构,因此训练精度更高和对输入的逼近误差更小,更能精确地表征输入信号的分布规律;并且在此基础上还提出了将网络的训练误差加入能量函数,同时自适应地获得能量控制参数λ1来引导曲线收敛于目标轮廓。此种有监督的主动轮廓模型对于目标概率分布未知、图像灰度不均匀和多像素目标等分割问题具有强大的理论意义和工程价值。以后的工作将围绕如何自动优化神经元个数,以及如何解决背景和前景灰度相差较小的图像分割方面展开研究,同时快速算法以及多目标分割问题也是尚待研究的领域。
[1] | CHAN T F, VESE L A. Active contours without edges[J]. IEEE Transactions on Image Processing, 2001, 10 (2) : 266-277. doi: 10.1109/83.902291 (0) |
[2] | LI C M, KAO C Y, GORE J C, et al. Minimization of region-scalable fitting energy for image segmentation[J]. IEEE Transactions on Image Processing, 2008, 17 (10) : 1940-1949. doi: 10.1109/TIP.2008.2002304 (0) |
[3] | CASELLES V, KIMMEL R, SAPIRO G. Geodesic active contours[J]. International Journal of Computer Vision, 1997, 22 (1) : 61-79. doi: 10.1023/A:1007979827043 (0) |
[4] | KIMMEL R, AMIR A, BRUCKSTEIN A M. Finding shortest paths on surfaces using level sets propagation[J]. IEEE Transactions on Pattern Analysis and Machine Intelligence, 1995, 17 (6) : 635-640. doi: 10.1109/34.387512 (0) |
[5] | LI C, GORE J C, DAVATZIKOS C. Multiplicative Intrinsic Component Optimization (MICO) for MRI bias field estimation and tissue segmentation[J]. Magnetic Resonance Imaging, 2014, 32 (7) : 913-923. doi: 10.1016/j.mri.2014.03.010 (0) |
[6] | ZHANG K, ZHANG L, LAM K M, et al. A local active contour model for image segmentation with intensity inhomogeneity [EB/OL]. [2015-10-10]. https://arxiv.org/abs/1305.7053. (0) |
[7] | WANG P, SUN K, CHEN Z. Local and global intensity information integrated geodesic model for image segmentation[C]//Proceedings of the 2012 International Conference on Computer Science and Electronics Engineering. Piscataway, NJ: IEEE, 2012, 2:129-132. (0) |
[8] | 敖谦, 朱燕平, 江少锋. 基于混合水平集的脑组织自动提取方法[J]. 计算机应用, 2013, 33 (7) : 2014-2017. ( AO Q, ZHU Y P, JIANG S F. Automatic brain extraction method based on hybrid level set model[J]. Journal of Computer Applications, 2013, 33 (7) : 2014-2017. doi: 10.3724/SP.J.1087.2013.02014 ) (0) |
[9] | WANG L, HE L, MISHRA A, et al. Active contours driven by local Gaussian distribution fitting energy[J]. Signal PRocessing, 2009, 89 (12) : 2435-2447. doi: 10.1016/j.sigpro.2009.03.014 (0) |
[10] | LIU S, PENG Y. A local region-based Chan-Vese model for image segmentation[J]. Pattern Recognition, 2012, 45 (7) : 2769-2779. doi: 10.1016/j.patcog.2011.11.019 (0) |
[11] | 林亚忠, 李新, 张会奇, 等. 快速的局部二值拟合优化分割算法[J]. 计算机应用, 2013, 33 (2) : 491-494. ( LIN Y Z, LI X, ZHANG H Q, et al. Fast local binary fitting optimization approach for image segmentation[J]. Journal of Computer Applications, 2013, 33 (2) : 491-494. doi: 10.3724/SP.J.1087.2013.00491 ) (0) |
[12] | VAN HULLE M M. Kernel-based topographic map formation achieved with normalized Gaussian competition[C]//Proceedings of the 2002 12th IEEE Workshop on Neural Networks for Signal Processing. Piscataway, NJ: IEEE, 2002:169-178. (0) |
[13] | ABDELSAMEA M M, GNECCO G, GABER M M. An efficient self-organizing active contour model for image segmentation[J]. Neurocomputing, 2015, 149 (PB) : 820-835. (0) |