0 引言

图像分割一直是研究人员关注的焦点，在许多领域中均有着较为广泛的应用。图像存在噪声点、密度不均等问题，使得图像分割在准确率、稳定性、实时性等方面仍然存在着很大的提升空间。随着人工智能的推广，脉冲耦合神经网络(Pulse Coupled Neural Network, PCNN)不需要学习和训练以及具有同步脉冲发放和全局耦合的特性，使其在图像分割领域得到了广泛的运用。Eckhorn等^[1]提出了基于猫的视觉皮层神经元信号传导特性的神经网络模型；针对PCNN模型参数的设置决定了其性能好坏的问题，Li等^[2]提出了一种免疫克隆算法来优化PCNN的参数调试问题，实现了PCNN参数的自适应调整；Zhou等^[3-5]提出了迭代分割由粗到细的策略，简化了PCNN的输入和动态阈值变化策略；于江波等^[6-9]确保了参数的正确调整和迭代结果的自动控制。这些研究使得PCNN能够自适应地去进行图像分割，为PCNN的应用提供了基础。

图像分割的实时性和稳定性均会影响相关人员的判断，并且实时性和稳定性也是图像分割的关键因素^[10]。景晓军等^[11]提出快速二维最大类间方差法(OTSU)，减小了时间复杂度。为了进一步提高算法对对比度和信噪比都较低的图像的分割效果，景晓军等^[12]又提出快速递推的三维最大类间方差分割算法，提高了算法的鲁棒性；范九伦等^[13]给出了三维最大类间方差分割算法的递推公式，在原来的基础上减少了计算时间。

上述算法完成了图像分割的基本要求，但是在很多领域中，图像分割更注重实时性、稳定性以及图像分割的准确率，所以在一定抗噪性能的基础上，如何有效地从原图中获取更多的信息量尤为重要。本文首先针对传统一维最大类间方差法(One-dimensional OTSU, 1OTSU)、2OTSU和3OTSU的图像分割结果携带原图信息量较少的问题，提出一维结合熵的最大类间方差法(One-dimensional OTSU-entropy, 1OTSU-H)、2OTSU-H和3OTSU-H准则，该准则可以获取原图较大的类间方差，拥有原图更多的信息量，具有较好的抗噪性能；其次为了进一步提高算法在图像分割中的实时性，提出了用粒子群优化(Particle Swarm Optimization, PSO)算法、布谷鸟搜索(Cuckoo Search, CS)算法、萤火虫算法(Firefly Algorithm, FA)、遗传算法(Genetic Algorithm, GA)四种群智能算法^[14-17]进行优化，将OTSU-H准则的目标函数作为四种群智能算法的适应度函数；最后将所提出的算法引入PCNN中，作为PCNN自动获取循环迭代次数的新准则，将搜索的全局最优阈值作为PCNN的动态阈值，进而实现基于区域的最佳自动图像分割。本文算法既保证了获得较大的类间方差的基础之上拥有原图更多的信息量，又减小了算法时间消耗，提高了图像分割的准确性。通过与原始的OTSU、最大熵准则以及近三年研究者们提出的五种最新图像分割算法对比，验证了本文算法在保证实时性的同时有着较好的图像分割效果。

1 OTSU-H算法

传统的OTSU算法的目标是获取原图最大的类间方差，保证图像中目标和背景有最大的区分度。但是该算法没有考虑图像熵，而在一些特定的领域(医学、航空和船舶等)图像分割中熵值和抗噪性能均是比较重要的因素，因此在保证拥有较大类间方差的基础上获取原图更多的信息量是本文算法研究的重点。

1.1 一维OTSU-H算法

针对1OTSU准则下图像分割携带原图信息量不足的问题，将一维图像熵引入到一维的目标函数中，这种结合了一维OTSU和一维图像熵的目标函数称为1OTSU-H算法。

一维熵H_s定义为：

$ {{H}_{s}}=(-\sum\limits_{s=0}^{L-1}{{{p}_{s}}\operatorname{l}\rm{b}\ }{{p}_{s}}) $

(1)

w₀(s)和w₁(s)分别为一维图像分割中目标区域和背景区域的灰度概率：

$ {{w}_{0}}\left( \mathit{s} \right)=\sum\limits_{i=0}^{s}{{{p}_{i}}} $

(2)

$ {{w}_{1}}\left( \mathit{s} \right)=\sum\limits_{i=s+1}^{L-1}{{{p}_{i}}} $

(3)

$ {{w}_{0}}\left( \mathit{s} \right)+{{w}_{1}}\left( \mathit{s} \right)=1 $

(4)

u₀(s)、u₁(s)和u_T(s)分别为一维图像分割中目标区域的均值、背景区域相应的均值和总体均值：

$ {{u}_{0}}\left( \mathit{s} \right)=\sum\limits_{i=0}^{s}{i{{p}_{i}}/{{w}_{0}}}\left( \mathit{s} \right) $

(5)

$ {{u}_{1}}\left( \mathit{s} \right)=\sum\limits_{i=s+1}^{L-1}{{}}i{{p}_{i}}/{{w}_{1}}\left( \mathit{s} \right) $

(6)

$ {{u}_{T}}\left( \mathit{s} \right)={{w}_{0}}\left( \mathit{s} \right){{u}_{0}}\left( \mathit{s} \right)+{{w}_{1}}\left( \mathit{s} \right){{u}_{1}}\left( \mathit{s} \right) $

(7)

一维图像分割中目标区域的类间方差、背景区域的类间方差和总的类间方差分别为式(8)、(9)和(10)：

$ {{\sigma }_{0}}\left( \mathit{s} \right)={{w}_{0}}\left( \mathit{s} \right){{({{u}_{0}}\left( \mathit{s} \right)-{{u}_{T}}\left( \mathit{s} \right))}^{2}} $

(8)

$ {{\sigma }_{1}}\left( \mathit{s} \right)={{w}_{1}}\left( \mathit{s} \right){{({{u}_{1}}\left( \mathit{s} \right)-{{u}_{T}}\left( \mathit{s} \right))}^{2}} $

(9)

$ {{t}_{r}}{{\sigma }_{b}}\left( s \right)={{\sigma }_{0}}\left( \mathit{s} \right)+{{\sigma }_{1}}\left( \mathit{s} \right) $

(10)

本文提出的1OTSU-H算法的目标函数及其对数表达式分别为式(11)和(12)：

$ h({{s}_{0}})=\max \left\{ {{({{t}_{r}}{{\sigma }_{b}}\left( s \right))}^{N_{\rm{var}}}}\times {{({{H}_{s}})}^{{{N}_{h}}}} \right\} $

(11)

$ \ln \left[h({{s}_{0}}) \right]=\max \left\{ {{N}_{\operatorname{var}}}\ln \left[{{t}_{r}}{{\sigma }_{b}}(s) \right]+{{N}_{h}}\ln \left[{{H}_{s}} \right] \right\} $

(12)

1.2 2OTSU-H算法

一维OTSU-H算法仅利用了图像灰度的分布信息，而没有考虑图像灰度的相关性，而许多图像的图像像元点之间是存在一定相关性的，因此需要构造二维观测空间，此时算法的实时性变差，基于以上考虑，本文又提出了2OTSU-H快速迭代算法。结合二维OTSU和二维熵的特性形成2OTSU-H算法。

设图像I大小为M*N，灰度级区间为(0, L-1)，像素点(s, t)的灰度值为f(s, t)，邻域均值为g(s, t)，g(s, t)的定义为：

$ g(s, t)=\frac{1}{{{k}^{2}}}\sum\limits_{i=-(k-1)/2}^{(k-1)/2}{\sum\limits_{j=-(k-1)/2}^{(k-1)/2}{f(s+i, t+j)}} $

(13)

g(s, t)满足：

$ 0\le g(s, t)\le \frac{1}{^{\mathop{k}^{2}}}\times \mathop{k}^{2}\times L=L $

(14)

即g(s, t)和f(s, t)有相同的灰度值区间。

由f(s, t)和g(s, t)构成二元组(i, j)，设p_ij为I的灰度联合概率分布，即：

$ {{p}_{ij}}={{n}_{ij}}/\left( M\times N \right) $

(15)

二维观测空间背景和目标的概率分别为式(16)、(17)：

$ {{w}_{0}}(s, t)=\sum\limits_{i=0}^{s}{\sum\limits_{j=0}^{t}{{{p}_{ij}}}} $

(16)

$ {{w}_{1}}(s, t)=\sum\limits_{i=s+1}^{L-1}{\sum\limits_{j=t+1}^{L-1}{{{p}_{ij}}}} $

(17)

阈值(s, t)把Cameraman图划分为如图 1(b)所示的四个部分; 图 1(c)是图像Ⅰ的二维灰度直方图，从图 1(c)可以看出，I灰度信息集中于区域1和区域3, 区域2和4是图像中存在的边缘和噪声信息，而且像素点很少可忽略不计。

$ {{\mathit{\boldsymbol{w}}}_{0}}+{{\mathit{\boldsymbol{w}}}_{1}}=1 $

(18)

二维直方图目标类均值、背景类均值及总体均值分别为式(19)、(20)和(21)：

$ {{\mathit{\boldsymbol{u}}}_{0}}={{({{u}_{0i}}, {{u}_{0j}})}^{\rm{T}}}={{\left[\sum\limits_{i=0}^{s}{\sum\limits_{j=0}^{t}{\frac{i{{p}_{ij}}}{{{w}_{0}}(s, t)}}}, \sum\limits_{i=0}^{s}{\sum\limits_{j=0}^{t}{\frac{j{{p}_{ij}}}{{{w}_{0}}(s, t)}}} \right]}^{\rm{T}}} $

(19)

$ {{\mathit{\boldsymbol{u}}}_{1}}={{({{u}_{1i}}, {{u}_{1j}})}^{\rm{T}}}={{\left[\sum\limits_{i=s+1}^{L-1}{\sum\limits_{j=s+1}^{L-1}{\frac{i{{p}_{ij}}}{{{w}_{1}}(s, t)}}}, \sum\limits_{i=s+1}^{L-1}{\sum\limits_{j=s+1}^{L-1}{\frac{j{{p}_{ij}}}{{{w}_{1}}(s, t)}}} \right]}^{\rm{T}}} $

(20)

$ \begin{align} & {{\mathit{\boldsymbol{u}}}_{T}}={{({{u}_{Ti}}, {{u}_{Tj}})}^{\rm{T}}}={{\left[\sum\limits_{i=0}^{s}{\sum\limits_{j=0}^{t}{i{{p}_{ij}}}}, \sum\limits_{i=0}^{s}{\sum\limits_{j=0}^{t}{j{{p}_{ij}}}} \right]}^{\rm{T}}}\approx \\ & \ \ \ \ \ \ \ \ \ \ {{w}_{0}}\left( \mathit{s}, \mathit{t} \right){{\mathit{\boldsymbol{u}}}_{0}}+{{w}_{1}}\left( \mathit{s}, \mathit{t} \right){{\mathit{\boldsymbol{u}}}_{1}} \\ & \\ \end{align} $

(21)

考虑到该算法实时性较差的问题，可用二维OTSU递推公式，表述为：

$ {{w}_{0}}(0, 0)={{p}_{00}} $

(22)

$ {{w}_{0}}(s, 0)={{w}_{0}}(s-1, 0)+{{p}_{s0}} $

(23)

$ {{w}_{0}}(0, t)={{w}_{0}}(0, t-1)+{{p}_{0t}} $

(24)

$ {{w}_{0}}(s, t)={{w}_{0}}(s, t-1)+{{w}_{0}}(s-1, t)-{{w}_{0}}(s-1, t-1)+{{p}_{st}} $

(25)

$ {{u}_{i}}(0, 0)=0 $

(26)

$ {{u}_{i}}(s, 0)={{u}_{i}}(s-1, 0)+s\cdot {{p}_{s0}} $

(27)

$ {{u}_{i}}(0, t)={{u}_{i}}(0, t-1)+s\cdot {{p}_{0t}} $

(28)

$ \begin{align} &{{u}_{i}}(s, t)={{u}_{i}}(s, t-1)+{{u}_{i}}(s-1, t)- \\ &\ \ \ \ \ \ \ \ \ \ \ \ \ \ {{u}_{i}}(s-1, t-1)+s\cdot {{p}_{st}} \\ \end{align} $

(29)

由式(22)~(25)可得w₀, 由式(26)~(29)可得u_i。同理可得：w₁和u_j。

二维观测空间中类间方差的迹为：

$ \begin{align} & {{t}_{r}}{{\sigma }_{b}}(s, t)= \\ & \frac{\left[{{({{w}_{0}}(s, t){{u}_{\mathit{T}i}}-{{u}_{i}}(s, t))}^{2}}+{{({{w}_{0}}(s, t){{u}_{\mathit{T}j}}-{{u}_{j}}(s, t))}^{2}} \right]}{{{w}_{0}}(s, t){{w}_{1}}(s, t)} \\ \end{align} $

(30)

二维熵为：

$ {{H}_{st}}=(-\sum\limits_{s=0}^{L-1}{\sum\limits_{t=0}^{L-1}{{{p}_{st}}\operatorname{l}{\rm{b}}\ \ {{p}_{st}}}}) $

(31)

本文提出的2OTSU-H算法的目标函数及其对数表达式分别为式(32)和(33):

$ h({{s}_{0}}, {{t}_{0}})=\max \left\{ {{\left( {{t}_{r}}{{\sigma }_{b}}(s, t) \right)}^{{{N}_{\operatorname{var}}}}}\times {{\left( {{H}_{st}} \right)}^{{{N}_{h}}}} \right\} $

(32)

$ \ln \left[h({{s}_{0}}, {{t}_{0}}) \right]=\max \left\{ {{N}_{\operatorname{var}}}\ln \left[{{t}_{r}}{{\sigma }_{b}}(s, t) \right]+{{N}_{h}}\ln \left[{{H}_{st}} \right] \right\} $

(33)

2OTSU-H快速迭代算法的时间复杂度由以前的O(L⁴)降到O(L²)。

1.3 3OTSU-H算法

在相同的参数下，2OTSU-H抗噪性能要优于1OTSU-H，而二维观测空间区域1及区域3概率和为1的假设普遍性不够，需要建立三维观测空间，而三维观测空间进行图像分割会提高计算复杂度，因此提出3OTSU-H算法以及3OTSU-H快速迭代算法。

任选一幅图像，像素点(s, t)对应的f(s, t)、g(s, t)和h(s, t)构成三元组(i, j, k)。其灰度联合概率分布为：

$ {{p}_{ijk}}={{n}_{ijk}}/\left( M\times N \right) $

(34)

背景和目标的概率分别为式(35)和式(36)：

$ {{w}_{0}}(s, t, q)=\sum\limits_{i=0}^{s}{\sum\limits_{j=0}^{t}{\sum\limits_{k=0}^{q}{{{p}_{ijk}}}}} $

(35)

$ {{w}_{1}}\left( s, t, q \right)=\sum\limits_{i=s+1}^{L-1}{\sum\limits_{j=t+1}^{L-1}{\sum\limits_{k=q+1}^{L-1}{{{p}_{ijk}}}}} $

(36)

阈值(s, t, q)把三维观测空间划分为如图 2(a)~(d)所示的8个区域，图像的灰度信息主要分布于区域1和8；边缘信息和噪声信息主要分布于其他6个区域。除了区域1和8之外，其余6个区域中像素点很少，可假设区域2~7上中像素点概率和为0。

三维直方图目标类均值、背景类均值及总体均值分别为式(37)、(38)和(39)：

$ \begin{align} & {{\mathit{\boldsymbol{u}}}_{0}}={{({{u}_{0i}}, {{u}_{0j}}, {{u}_{0k}})}^{\rm{T}}}=\left[\sum\limits_{i=0}^{s}{\sum\limits_{j=0}^{t}{\sum\limits_{k=0}^{q}{\frac{i{{p}_{ijk}}}{{{w}_{0}}\left( s, t, q \right)}}}}, \right. \\ & \;\;\;\sum\limits_{i=0}^{s}{\sum\limits_{j=0}^{t}{\sum\limits_{k=0}^{q}{\frac{j{{p}_{ijk}}}{{{w}_{0}}\left( s, t, q \right)}}}}, {{\left. \sum\limits_{i=0}^{s}{\sum\limits_{j=0}^{t}{\sum\limits_{k=0}^{q}{\frac{k{{p}_{ijk}}}{{{w}_{0}}\left( s, t, q \right)}}}} \right]}^{\rm{T}}} \\ \end{align} $

(37)

$ \begin{align} & {{\mathit{\boldsymbol{u}}}_{1}}={{({{u}_{1i}}, {{u}_{1j}}, {{u}_{1k}})}^{\rm{T}}}=\left[\sum\limits_{i=s+1}^{L-1}{\sum\limits_{j=t+1}^{L-1}{\sum\limits_{k=q+1}^{L-1}{\frac{i{{p}_{ijk}}}{{{w}_{1}}\left( s, t, q \right)}}}}, \right. \\ & \;\;\;\sum\limits_{i=s+1}^{L-1}{\sum\limits_{j=t+1}^{L-1}{\sum\limits_{k=q+1}^{L-1}{\frac{j{{p}_{ijk}}}{{{w}_{1}}\left( s, t, q \right)}}}}, {{\left. \sum\limits_{i=s+1}^{L-1}{\sum\limits_{j=t+1}^{L-1}{\sum\limits_{k=q+1}^{L-1}{\frac{k{{p}_{ijk}}}{{{w}_{1}}\left( s, t, q \right)}}}} \right]}^{\rm{T}}} \\ \end{align} $

(38)

$ \begin{align} & {{\mathit{\boldsymbol{u}}}_{T}}={{({{u}_{Ti}}, {{u}_{Tj}}, {{u}_{Tk}})}^{\rm{T}}}=\left[\sum\limits_{i=0}^{L-1}{\sum\limits_{j=0}^{L-1}{\sum\limits_{k=0}^{L-1}{i{{p}_{ijk}}}}} \right., \sum\limits_{i=0}^{L-1}{\sum\limits_{j=0}^{L-1}{\sum\limits_{k=0}^{L-1}{j{{p}_{ijk}}}}},\\ & \;\;\;\; {{\left. \sum\limits_{i=0}^{L-1}{\sum\limits_{j=0}^{L-1}{\sum\limits_{k=0}^{L-1}{k{{p}_{ijk}}}}} \right]}^{\rm{T}}}\approx {{w}_{0}}\left( \mathit{s}, \mathit{t}, \mathit{q} \right){{\mathit{\boldsymbol{u}}}_{0}}+{{w}_{1}}\left( \mathit{s}, \mathit{t}, \mathit{q} \right){{\mathit{\boldsymbol{u}}}_{1}} \\ \end{align} $

(39)

传统3OTSU时间复杂度高，可由三维递推公式表述为：

$ {{w}_{0}}\left( 0, 0, 0 \right)={{p}_{000}} $

(40)

$ {{w}_{0}}\left( s, 0, 0 \right)={{w}_{0}}\left( s-1, 0, 0 \right)+{{p}_{s00}} $

(41)

$ {{w}_{0}}\left( 0, t, 0 \right)={{w}_{0}}\left( 0, t-1, 0 \right)+{{p}_{0t0}} $

(42)

$ {{w}_{0}}\left( 0, 0, q \right)={{w}_{0}}\left( 0, 0, q-1 \right)+{{p}_{00q}} $

(43)

$ \begin{align} &\ \ \ \ \ \ \ \ \ \ \ \ \ \ {{w}_{0}}\left( s, t, 0 \right)= \\ &\ \ \ {{w}_{0}}\left( s-1, t, 0 \right)+{{w}_{0}}\left( s, t-1, 0 \right)- \\ &\ \ \ {{w}_{0}}\left( s-1, t-1, 0 \right)+{{p}_{st0}} \\ \end{align} $

(44)

$ \begin{align} &{{w}_{0}}\left( s, 0, q \right)={{w}_{0}}\left( s-1, 0, q \right)+{{w}_{0}}\left( s, 0, q-1 \right)- \\ &\ \ \ \ \ \ {{w}_{0}}\left( s-1, 0, q-1 \right)+{{p}_{s0q}} \\ \end{align} $

(45)

$ \begin{align} &{{w}_{0}}\left( 0, t, q \right)={{w}_{0}}\left( 0, t-1, q \right)+{{w}_{0}}\left( 0, t, q-1 \right)- \\ &\ \ \ \ \ \ \ \ \ {{w}_{0}}\left( 0, t-1, q-1 \right)+{{p}_{0tq}} \\ \end{align} $

(46)

$ \begin{align} &{{w}_{0}}\left( s, t, q \right)={{w}_{0}}\left( s-1, t, q \right)+{{w}_{0}}\left( s, t-1, q \right)+ \\ &\ \ {{w}_{0}}\left( s, t, q-1 \right)-{{w}_{0}}\left( s, t-1, q-1 \right)-{{w}_{0}}\left( s-1, t, q-1 \right)- \\ &\ \ {{w}_{0}}\left( s-1, t-1, q \right)-{{w}_{0}}\left( s-1, t-1, q-1 \right)+{{p}_{stq}} \\ \end{align} $

(47)

$ {{u}_{i}}\left( 0, 0, 0 \right)=0 $

(48)

$ {{u}_{i}}\left( s, 0, 0 \right)={{u}_{i}}\left( s-1, 0, 0 \right)+s\cdot {{p}_{s00}} $

(49)

$ {{u}_{i}}\left( 0, t, 0 \right)={{u}_{i}}\left( 0, t-1, 0 \right)+0\cdot {{p}_{0t0}} $

(50)

$ {{u}_{i}}\left( 0, 0, q \right)={{u}_{i}}\left( 0, 0, q-1 \right)+0\cdot {{p}_{00q}} $

(51)

$ \begin{align} &{{u}_{i}}\left( s, t, 0 \right)={{u}_{i}}\left( s-1, t, 0 \right)+{{u}_{i}}\left( s, t-1, 0 \right)- \\ &\ \ \ \ {{u}_{i}}\left( s-1, t-1, 0 \right)+s\cdot {{p}_{st0}} \\ \end{align} $

(52)

$ \begin{align} &{{u}_{i}}\left( s, 0, q \right)={{u}_{i}}\left( s-1, 0, q \right)+{{u}_{i}}\left( s, 0, q-1 \right)- \\ &\ \ \ \ \ \ \ \ \ {{u}_{i}}\left( s-1, 0, q-1 \right)+s\cdot {{p}_{s0q}} \\ \end{align} $

(53)

$ \begin{align} &{{u}_{i}}\left( 0, t, q \right)={{u}_{i}}\left( 0, t-1, q \right)+{{u}_{i}}\left( 0, t, q-1 \right)- \\ &\ \ \ \ \ \ \ \ {{u}_{i}}\left( 0, t-1, q-1 \right)+s\cdot {{p}_{0tq}} \\ \end{align} $

(54)

$ \begin{align} &{{u}_{i}}\left( s, t, q \right)={{u}_{i}}\left( s-1, t, q \right)+{{u}_{i}}\left( s, t-1, q \right)+ \\ &\ \ {{u}_{i}}\left( s, t, q-1 \right)-{{u}_{i}}\left( s, t-1, q-1 \right)-{{u}_{i}}\left( s-1, t, q-1 \right)- \\ &\ \ {{u}_{i}}\left( s-1, t-1, q \right)+{{u}_{i}}\left( s-1, t-1, q-1 \right)+s\cdot {{p}_{stq}} \\ \end{align} $

(55)

由式(40)~(47)可得w₀, 由式(48)~(55)可得u_i。同理可得：w₁、u_j、u_k。

三维观测空间类间方差的迹如式(56)：

$ \begin{align} & {{t}_{r}}{{\sigma }_{b}}\left( s, t, q \right)=\left[{{\left( {{w}_{0}}\left( s, t, q \right){{u}_{Ti}}-{{u}_{i}}\left( s, t, q \right) \right)}^{2}}+ \right. \\ &\;\; {{\left( {{w}_{0}}\left( s, t, q \right){{u}_{Tj}}-{{u}_{j}}\left( s, t, q \right) \right)}^{2}}+ \\ &\;\; \left. {{\left( {{w}_{0}}\left( s, t, q \right){{u}_{Tk}}-{{u}_{k}}\left( s, t, q \right) \right)}^{2}} \right]/\left[{{w}_{0}}\left( s, t, q \right){{w}_{1}}\left( s, t, q \right) \right] \\ \end{align} $

(56)

三维熵为：

$ {{H}_{stq}}=(-\sum\limits_{s=0}^{L-1}{\sum\limits_{t=0}^{L-1}{\sum\limits_{q=0}^{L-1}{{{p}_{stq}}}}}\rm{lb}\ {{\mathit{p}}_{\mathit{stq}}}) $

(57)

本文提出的3OTSU-H算法的目标函数及其对数表达式分别为式(58)、(59)：

$ h\left( {{s}_{0}}, {{t}_{0}}, {{q}_{0}} \right)=\max \left\{ {{\left( {{t}_{r}}{{\sigma }_{b}}\left( s, t, q \right) \right)}^{{{N}_{\operatorname{var}}}}}\times {{\left( {{H}_{stq}} \right)}^{{{N}_{h}}}} \right\} $

(58)

$ \ln \left[h\left( {{s}_{0}}, {{t}_{0}}, {{q}_{0}} \right) \right]=\max \left\{ {{N}_{\operatorname{var}}}\ln \left[{{t}_{r}}{{\sigma }_{b}}\left( s, t, q \right) \right]+\\ \;\;{{N}_{h}}\ln \left[{{H}_{stq}} \right] \right\} $

(59)

3OTSU-H快速迭代算法的时间复杂度由以前的O(L⁶)降到O(L³)。

2 四种群智能算法优化的OTSU-H

虽然2OTSU-H和3OTSU-H快速递推算法与传统的2OTSU和3OTSU算法相比实时性有了很大程度的提升，但是在实际应用中实时性和准确率的问题仍待解决。群智能算法的出现，不仅可以提高算法的实时性，而且能搜索到OTSU-H算法的全局最优值，避免陷入OTSU-H算法的局部最优值，提高算法的准确率。因此本文提出用四种群智能算法分别一维、二维及三维OTSU-H算法进行优化, 其中对3OTSU-H算法的优化过程如下：

2.1 布谷鸟搜索算法优化OTSU-H(CS-OTSU-H)

CS算法根据布谷鸟以巢寄生方式提高自己的繁殖率，并结合Lévy飞行行为更新鸟巢位置，获得全局最优值。用该算法优化时，不需要重新匹配大量参数，需要初始化算法的基本参数，将OTSU-H算法的目标函数作为布谷鸟算法的适应度函数。

算法1   CS-3OTSU-H算法。

输入  鸟巢数量n=50，鸟巢替换的概率Pa=0.25，搜索范围[Lb, Ub]为[1,255], 参数N_var=1和N_h=1；

输出  最优鸟巢(s₀, t₀, q₀)。

3OTSU-H目标函数F_s0t0q0=ln[h(s₀, t₀, q₀)];

初始化n个鸟巢位置s_it_iq_i(i=1, 2, …, n)；

while (t＜最大迭代次数)

  通过Lévy飞行得到鸟巢i；

  从n个鸟巢中随机选择一个鸟窝j；

  评估适应度F_sitiqi、F_sjtjqj；

  if (F_sitiqi＞F_sjtjqj)

   用新鸟巢i替换j；

  end

  抛弃并更新较差鸟窝的Pa；

  找出当前最优鸟巢；

end

输出CS-3OTSU-H的最优阈值(s₀, t₀, q₀)。

2.2 粒子群优化算法优化OTSU-H(PSO-OTSU-H)

PSO算法依据鸟类在捕食过程中，个体食物源位置以及离食物的距离到底有多远，为了找到食物，搜索食物源附近的个体周围区域，实现了信息的共享，每个个体建立自身的信念，当发现其他个体的适应度函数值高于自己时，个体会调整自己的速度和位置，从而以提高自身的搜索能力。该算法在优化1OTSU-H、2OTSU-H和3OTSU-H时，将OTSU-H算法的目标函数作为PSO的适应度函数。

算法2   PSO-3OTSU-H算法流程。

输入  粒子个数n=50, bird_setp= 20, 相关参数c₁=2、c₂=2、w=0.9，产生的位置和速度, 参数N_var=1和N_h=1；

输出  最优种群位置(s₀, t₀, q₀)。

3OTSU-H目标函数F_s0t0q0=ln[h(s₀, t₀, q₀)];

初始化n个种群位置和速度；

while(t＜最大迭代次数)

  计算F_stq；

  记录每个个体取最好适应度值的位置；

  与历史最优个体适应度值进行比较；

  if(F_stq＞=F_max)

   记录F_stq的位置，作为全局最优解；

  end

  找出当前最优种群位置(s₀, t₀, q₀)；

end

输出PSO-3OTSU-H的最优阈值(s₀, t₀, q₀)。

2.3 萤火虫算法优化OTSU-H(FA-OTSU-H)

FA算法反映每个萤火虫被绝对亮度比它大的萤火虫所吸引，从而更新更新自身的位置。在优化OTSU-H时，初始化算法参数，将OTSU-H的目标函数作为FA的适应度函数。

算法3   FA-3OTSU-H算法。

输入  粒子个数n=50, MaxGeneration=20, 相关参数α=0.5, β=0.2, γ=1, 参数N_var=1和N_h=1；

输出  最优种群位置(s₀, t₀, q₀)。

3OTSU-H目标函数F_s0t0q0=ln[h(s₀, t₀, q₀)];

随机产生n只萤火虫，确定每个萤火虫的绝对亮度I(x_i)；

while(t＜最大迭代次数)

 for i=1:n

  for j=1:n

   if I(x_i)＞I(x_j)

    计算i和萤火虫j的距离r_ij;

    计算i和萤火虫j的引力β;

    把萤火虫j向i移动;

    计算更新新的萤火虫的亮度；

   end

  end

 end

end

输出FA-3OTSU-H的最优阈值(s₀, t₀, q₀)。

2.4 遗传算法优化OTSU-H(GA-OTSU-H)

GA是模仿“物竞天择，适者生存”的演化法则。将解空间数据表示成遗传空间中染色体及染色体上的基因，基因在染色体上数据的不同组合构成不同的点, 通过对位置的编码、选择、交叉、变异和解码找到问题的最优解。在优化OTSU-H的过程中，需要初始化遗传代数、种群数目、选择、交叉和变异的概率，将OTSU-H的目标函数作为GA的适应度函数。

算法4   GA-3OTSU-H算法。

输入  选择、交叉、变异的概率分别为0.9，0.7和0.05，染色体数目40，最大遗传代数为500, 参数N_var=1和N_h=1；

输出  最优染色体(s₀, t₀, q₀)。

3OTSU-H目标函数F_s0t0q0=ln[h(s₀, t₀, q₀)];

设置最大迭代次数，初始化n个染色体s_i(i=1, 2, …, n)；

while(t＜最大迭代次数)

  从n个染色体中随机选择一个染色体；

  评估适应度F_sitiqi、F_sjtjqj；

  if(F_sitiqi>F_sjtjqj)

   用新染色体i替换j；

  end

  找出当前最优染色体(s₀, t₀, q₀)；

end

输出GA-3OTSU-H的最优阈值(s₀, t₀, q₀)。

3 PCNN模型中应用

在其他参数固定的前提下，PCNN循环迭代次数的选取影响着图像的分割效果。在现有的研究中，马义德等^[5]提出了基于图像熵准则的自动图像分割算法。于江波等^[6]研究了PCNN模型中除了循环迭代次数以外的其他参数在图像处理中的确定问题。最大熵准则能获取原图最大的信息量，但是抗噪性能确不理想；最大类间方差准的抗噪能力较强，但是从原图中获取的信息量不足：因此最大熵和最大类间方差准则下图像分割效果均不理想。本文的2OTSU-H、二维最大熵和2OTSU算法的分割效果如图 3所示。

基于以上考虑，将群智能算法优化后的OTSU-H准则作为PCNN循环迭代次数选取的准则。

PCNN简化模型其公式如下：

$ {{F}_{ij}}\left( n \right)={{I}_{ij}} $

(60)

$ {{L}_{ij}}\left( n \right)=\sum\limits_{k,l}{{{W}_{kl}}{{Y}_{kl}}\left( n-1 \right)} $

(61)

$ {{U}_{ij}}\left( n \right)={{F}_{ij}}\left( n \right)\left( 1+b{{L}_{ij}}\left( n \right) \right) $

(62)

$ {{T}_{ij}}\left( n \right)={{\rm{e}}^{-a}}{{T}_{ij}}\left( n-1 \right)+V{{Y}_{ij}}\left( n \right) $

(63)

$ {Y_{ij}}\left( n \right) = \left\{ \begin{array}{l} 1, \;\;\;\;\;\;\;\;\;\;\;{U_{ij}}\left( n \right){\rm{ > }}{T_{ij}}\left( {n - 1} \right)\\ 0, \;\;\;\;\;\;\;\;\;\;\;{\rm{其他}} \end{array} \right. $

(64)

其中参数$ \mathit{\boldsymbol{W}}=\left( \begin{matrix} \sqrt{1/2} & 1 & \sqrt{1/2} \\ 1 & 1 & 1 \\ \sqrt{1/2} & 1 & \sqrt{1/2} \\ \end{matrix} \right)$; a=0.1; b=0.2;V=0.2。

任选一幅图像，像素(i, j)的原始灰度值是I_ij, F_ij为反馈域输入通道, L_ij为链接域通道, b为链接的强度系数, U_ij位内部状态信号, T_ij为门限阈值, Y_ij为脉冲输出。简化的PCNN模型中，b是可变的，b值越大，图像分解得越精细，图像的轮廓越明显，反之图像的细节越丰富。

4 实验结果及分析

为了验证本文算法的性能, 选取了大量图像在Matlab R2016a实验平台上进行测试。计算机的配置为3.5 GHz Intel Core i3-4150 CPU和8.0 GB内存。

4.1 OTSU-H的抗噪性能

对cameraman图分别加入高斯噪声、椒盐噪声和混合噪声，在2OTSU-H和3OTSU-H准则下进行抗噪性能测试, 测试结果如图 4所示。从图 4中可以看出，在噪声类型不同的情况下，3OTSU-H的抗噪性能均优于2OTSU-H，对大部分噪声都有着较好的抑制能力，因此在相同的参数下3OTSU-H的抗噪性能优于2OTSU-H，同时2OTSU-H和3OTSU-H都拥有较好的抗噪性能。

4.2 OTSU-H下的主客观评价

在二维观测空间和三维观测空间中，分别用最大类间方差准则、OTSU-H和最大熵准则对cameraman和lena进行图像分割，分割结果如图 5~6所示。

在表 1中，参数为(1, 0)的情况下表示最大类间方差准则；参数为(0, 1)情况下表示最大熵准则; 参数为(1, 1)和(1, 2)是结合了最大类间方差和图像熵的OTSU-H准则。在OTSU-H准则下，分割之后的图像所携带的信息量均有一定的提升。从图 5~6中可以看出OTSU-H的图像分割效果要优于OTSU和最大熵准则下的图像分割效果。

4.3 群智能算法优化后OTSU-H的主客观评价

为了进一步降低OTSU-H的时间复杂度，提出了用四种群智能算法分别对2OTSU-H和3OTSU-H进行优化，统计各算法的时间复杂度、阈值和灰度熵, 不同参数下的结果如表 2所示。

表 2中数据是在cameraman图中测试得到，在二维算法优化中，PSO使2OTSU-H算法的实时性提高了41.8%，CS使2OTSU-H算法的实时性提高了28.2%，FA使2OTSU-H算法的实时性提高了26.3%，GA使2OTSU-H算法的实时性提高了19.5%。在三维算法优化中，PSO使3OTSU-H算法的实时性提高了44.4%，CS使3OTSU-H算法的实时性提高了43.8%，FA使3OTSU-H算法的实时性提高了24.3%，GA使3OTSU-H算法的实时性提高了22.1%。因此，二维图像分割算法中PSO-2OTSU-H算法的实时性最强，三维图像分割算法中PSO-3OTSU-H和CS-3OTSU-H算法的实时性最强。综上所述，基于粒子群算法和布谷鸟算法对2OTSU-H、3OTSU-H的优化性能最好，实时性最强。

4.4 OTSU-H与PCNN模型融合

对cameraman和lena进行分割，比较在三种准则下PCNN的分割效果。从图 7~8可以看出，OTSU-H与PCNN融合的的一维、二维和三维图像分割效果优于最大类间方差及最大熵准则。

4.5 OTSU-H准则的实用性验证

采用最大熵准则、OTSU和OTSU-H准则分别在二维观测空间中和三维观测空间中，对任意两幅图像(一幅医学图像，一幅航空图)进行图像分割，分割结果如图 9所示。

从图 9可以看出，最大熵准则和最大类间方差准则都不同程度上出现了错分，而本文算法在医学和航空图像中分割效果比最大类间方差和最大熵分割效果好。

从表 2可以看出，基于粒子群算法和布谷鸟算法优化后的OTSU-H拥有最好的实时性。因此PSO-OTSU-H及CS-OTSU-H实时性和准确率的提高，可以给医生和飞行员提供实时性的判断。

4.6 本文算法与其他算法对比实验

图像分割方法主要有图论分割、像素的聚类分割、基于候选区域的语义分割^[18]。Li等^[19]提出基于图论的线性谱聚类(Linear Spectral Clustering, LSC)算法，迭代地使用k-means算法在高维特征空间的聚类作用。Achanta等^[20]改进了更为复杂的距离计算公式，提出了基于像素聚类的自适应简单线性迭代聚类(Adaptively Simple Linear Iterative Clustering, ASLIC)算法。Long等^[21]提出了全卷积网络(Fully Convolutional Network, FCN)。Hariharan等^[22]借用图像金字塔的思想，提出了Hypercolumn算法。Zheng等^[23]提出了RNN(Recurrent Neural Network)算法。测试的数据集为模式分析统计建模和计算学习视觉对象类(Pattern Analysis Statistical modeling and Computational Learning Visual Object Classes, PASCAL VOC)^[24]。不同算法的分割效果对比如图 10所示，相应的性能评价结果如表 3所示。

图 10为本文算法以及LSC、ASLIC、FCN、Hypercolumn和RNN等算法的分割效果。从图 10可以看出，图 10(b)、(d)、(e)和(f)分割效果最好。

从表 3可以看出，本文算法、FCN和Hypercolumn进行图像分割时拥有原图较多的信息量，虽然本文算法拥有的信息量略低于FCN，但是本文算法模型简单，不需要训练，而最新算法如FCN、Hypercolumn和RNN都是基于卷积神经网络，模型复杂，训练时间长，且测试时间也相对较长。因此可以得出，本文算法在保证实时性的同时有着较好的图像分割效果。

5 结语

本文提出了基于群智能算法优化的OTSU-H与PCNN融合的图像分割算法。该算法根据OTSU-H准则来确定PCNN的最佳迭代次数，并且运用群智能算法来优化融合算法的目标函数，减少算法的时间损耗。在算法的实现过程中，首先构造了观测空间；其次给出了融合算法的快速递推公式；最后用群智能算法优化目标函数。与OTSU算法和最大熵准则相比，本文算法既有着较大的类间方差也携带原图较多的信息量。与图论分割、像素的聚类分割和基于候选区域的语义分割的分割效果相比，本文算法时间损耗少，携带信息量较大，同时稳定性好。在实验过程中，对本文算法的实时性进行了评测对实用性进行了验证。实验结果表明：本文算法在低对比度和低信噪比的图像中，仍具有较好的图像分割效果, 同时本文算法减少了所需的计算量，节约了计算机的存储空间，具有较强的抗干扰能力，在保证实时性的同时有着较好的图像分割效果，可以运用于医学、航空等领域。下一步将研究，在保证一定实时性的同时，与基于语义的图像分割算法相结合，进行多信息融合的图像分割，进一步提高算法的精度。