0 引言

图像分割是将图像中的感兴趣区域提取出来的图像处理过程。目前有多种分割方法, 其中活动轮廓模型是公认的一种高效的分割方法^[1-2], 其基本思想是最小化一个关于活动轮廓曲线的能量泛函, 最终得到目标物体的边界。活动轮廓模型可以归结为两类:一类是基于边缘的活动轮廓分割模型^[3-8]; 一类是基于区域的活动轮廓模型^[9-15]。前者的基本思想是用边缘检测函数驱使活动轮廓朝图像边界移动, 最初由文献[6]在snake模型中提出来的, 而基于区域的分割方法的思想则主要集中在用全局信息识别目标边缘。

全局图像分割方法是将图像中所有物体按照模型设定全部分割出来的方法。然而, 在一些图像分割问题中, 只想分割其中的某一个物体而非全部物体, 这就是所谓的选择性分割。比如在一幅肾脏CT图像中, 只需要分割其中的一颗肾脏, 则全局分割方法将不再适用, 因为这些模型会将其它有相同特征(纹理、颜色)的区域一并分割出来。但是给予一些几何约束条件后, 本研究提出的方法可以实现对其中一颗肾脏的分割。

文献[16]提出了灰度图像的单水平集选择性分割模型, 该模型用变分思想求能量泛函极小值; 与双水平集方法相比, 单水平集方法分割速度更快。考虑到对彩色图像分割的研究具有更广阔的应用前景, 本研究基于Lavdie-Chen的模型提出带几何约束的彩色图像分割方法。该方法沿用单水平集的思想, 对于彩色图像, 传统的处理方法是将三个通道分别处理, 但是对于基于边缘的分割方法, 若分别对彩色图像三通道进行分割可能得不到统一的分割结果。因此, 该模型将彩色图像作为整体求图像梯度和边缘检测函数; 边缘检测函数中引入一个可以控制分割速度的参数; 然后利用边缘检测函数和关于靠近目标物体边界的约束点的距离函数以及约束点形成的多边形内外面积共同控制曲线演化, 最终实现彩色图像的选择性分割。试验仿真表明, 本研究提出的方法不仅对有接触的物体可以进行选择性分割, 对有重叠区域的图像同样能进行选择性分割; 此外, 对医学彩色CT图像中边界不清晰的器官也能进行选择性分割。

1 灰度图像单水平集选择性分割模型

对于一幅定义在矩形域Ω上的灰度图像u(x, y), 在感兴趣区域边界附近给定n个约束点, 用A={(x_i, y_i)∈Ω, 1≤i≤n}⊂Ω表示。Lavdie-Chen提出灰度图像的单水平集选择性分割模型

$ \begin{array}{*{20}{c}} {\mathop {\min }\limits_{\mathit{\Gamma },{c_2}} F\left( {\mathit{\Gamma }{c_2}} \right) = \mu \int\limits_\mathit{\Gamma } {g\left( {\left| {\nabla u\left( {x,y} \right)} \right|} \right){\rm{d}}s} + {\lambda _1}\int\limits_{{\rm{inside}}\left( \Gamma \right)} {{{\left| {u\left( {x,y} \right) - {c_1}} \right|}^2}{\rm{d}}x{\rm{d}}y} + }\\ {{\lambda _2}\int\limits_{{\rm{outside}}\left( \mathit{\Gamma } \right)} {{{\left| {u\left( {x,y} \right) - {c_2}} \right|}^2}{\rm{d}}x{\rm{d}}y} + \nu \left[ {{{\left( {\int\limits_{{\rm{inside}}\left( \mathit{\Gamma } \right)} {{\rm{d}}x{\rm{d}}y} - {A_1}} \right)}^2} + {{\left( {\int\limits_{{\rm{outside}}\left( \mathit{\Gamma } \right)} {{\rm{d}}x{\rm{d}}y} - {A_2}} \right)}^2}} \right],} \end{array} $

(1)

式中: μ, λ₁, λ₂, ν为非负加权参数; g(|∇u(x, y)|)为边缘检测函数; c₁, c₂分别为曲线Γ内外区域的灰度值; Α₁, Α₂分别为约束点确定的多边形内外区域面积。边缘检测函数的数学表达式如下所示:

$ g\left( {\left| {\nabla u\left( {x,y} \right)} \right|} \right) = \frac{1}{{1 + {{\left| {\nabla u\left( {x,y} \right)} \right|}^2}}} = \frac{1}{{1 + \left( {u_x^2 + u_y^2} \right)}}。$

(2)

2 带约束的彩色图像选择性分割模型

对于定义在矩形域Ω上的彩色图像u(x, y, 3), 给定n个约束点, 即A={(x_i, y_i)∈Ω, 1≤i≤n}⊂Ω, 约束点的数量与位置的选取取决于所要分割的目标物体, 一般情况下数量选3~6个即可, 若目标物体较为模糊, 初始约束点选取在接近边缘位置的分割效果要优于其它位置; 若目标物体较为清晰, 初始点可任意选择在待分割物体上。据此提出彩色图像选择性分割方法

$ \begin{array}{*{20}{c}} {\mathop {\min }\limits_{\mathit{\Gamma },c_1^{\left( i \right)},c_2^{\left( i \right)}} F\left( {\mathit{\Gamma },c_1^{\left( i \right)},c_2^{\left( i \right)}} \right) = \mu \int\limits_\mathit{\Gamma } {d\left( {x,y} \right){g_{{\text{color}}}}{\text{d}}s} + \frac{1}{3}\iint\limits_{{\text{inside}}\left( \Gamma \right)} {\sum\limits_{i = 1}^3 {\lambda _1^{\left( i \right)}{{\left| {{u^{\left( i \right)}} - c_1^{\left( i \right)}} \right|}^2}{\text{d}}x{\text{d}}y} } + } \\ {\frac{1}{3}\iint\limits_{{\text{outside}}\left( \mathit{\Gamma } \right)} {\sum\limits_{i = 1}^3 {\lambda _2^{\left( i \right)}{{\left| {{u^{\left( i \right)}} - c_2^{\left( i \right)}} \right|}^2}{\text{d}}x{\text{d}}y} } + \nu \left[ {{{\left( {\iint\limits_{{\text{inside}}\left( \Gamma \right)} {{\text{d}}x{\text{d}}y} - {A_1}} \right)}^2} + {{\left( {\iint\limits_{{\text{outside}}\left( \Gamma \right)} {{\text{d}}x{\text{d}}y} - {A_2}} \right)}^2}} \right],} \end{array} $

(3)

式中: μ, λ₁, λ₂, ν为非负加权参数; c₁⁽ⁱ⁾, c₂⁽ⁱ⁾, i=1, 2, 3分别为曲线Γ内外区域的三通道图像灰度值; Α₁, Α₂分别为约束点确定的多边形内外区域的面积。g_color为彩色图像边缘检测函数, d(x, y)^[17]为距离函数, 定义如下:

$ d\left( {x,y} \right) = \prod\limits_{i = 1}^n {\left( {1 - \exp \left( { - \frac{{{{\left( {x - {x_i}} \right)}^2}}}{{2{\sigma ^2}}}} \right)\exp \left( { - \frac{{{{\left( {y - {y_i}} \right)}^2}}}{{2{\sigma ^2}}}} \right)} \right)} ,\;\;\forall \left( {x,y} \right) \in \mathit{\Omega }。$

(4)

2.1 彩色图像边缘检测函数

根据黎曼几何的思想, 彩色图像u(x, y, 3)可以看成是Euclidean空间中以(x, y)为参数的曲面, 则曲面上任意方向的弧长微元^[16] ${\rm{d}}u\left( {\frac{{\partial {u^{\left( 1 \right)}}}}{{\partial x}}{\rm{d}}x + \frac{{\partial {u^{\left( 1 \right)}}}}{{\partial y}}{\rm{d}}y, \frac{{\partial {u^{\left( 2 \right)}}}}{{\partial x}}{\rm{d}}x + \frac{{\partial {u^{\left( 2 \right)}}}}{{\partial y}}{\rm{d}}y, \frac{{\partial {u^{\left( 3 \right)}}}}{{\partial x}}{\rm{d}}x + \frac{{\partial {u^{\left( 3 \right)}}}}{{\partial y}}{\rm{d}}y} \right) $, 则

$ {\left| {{\rm{d}}u} \right|^2} = \left\langle {{\rm{d}}u,{\rm{d}}u} \right\rangle = \sum\limits_{i = 1}^3 {\left( {\frac{{\partial {u^{\left( i \right)}}}}{{\partial x}}{\rm{d}}x + \frac{{\partial {u^{\left( i \right)}}}}{{\partial y}}{\rm{d}}y} \right)} = {C_1}{\left( {{\rm{d}}x} \right)^2} + 2{C_2}{\rm{d}}x{\rm{d}}y + {C_3}{\left( {{\rm{d}}y} \right)^2}, $

式中:

$ {C_1} = \sum\limits_{i = 1}^3 {{{\left( {\frac{{\partial {u^{\left( i \right)}}}}{{\partial x}}} \right)}^2}} ,{C_2} = \sum\limits_{i = 1}^3 {\left( {\frac{{\partial {u^{\left( i \right)}}}}{{\partial x}}} \right)\left( {\frac{{\partial {u^{\left( i \right)}}}}{{\partial y}}} \right)} ,{C_3} = \sum\limits_{i = 1}^3 {{{\left( {\frac{{\partial {u^{\left( i \right)}}}}{{\partial y}}} \right)}^2}} 。$

(5)

式(5)可改写成二次型形式${\left| {{\rm{d}}u} \right|^2} = {\left[\begin{array}{l} {\rm{d}}x\\ {\rm{d}}y \end{array} \right]^{\rm{T}}}\left[{\begin{array}{*{20}{c}} {{C_1}}&{{C_2}}\\ {{C_2}}&{{C_3}} \end{array}} \right]\left[\begin{array}{l} {\rm{d}}x\\ {\rm{d}}y \end{array} \right] = {\left[\begin{array}{l} {\rm{d}}x\\ {\rm{d}}y \end{array} \right]^{\rm{T}}}\mathit{\boldsymbol{A}}\left[\begin{array}{l} {\rm{d}}x\\ {\rm{d}}y \end{array} \right] $, 可求得矩阵A的特征值为

$ {\sigma _{1,2}} = \left( {{C_1} + {C_3} \pm \sqrt {{{\left( {{C_1} - {C_3}} \right)}^2} + 4C_2^2} } \right)/2,{\sigma _1} > {\sigma _2}。$

(6)

则彩色图像边缘检测函数

$ {g_{{\rm{color}}}} = g\left( {{\sigma _1} - {\sigma _2}} \right) = 1/\left[ {1 + {{\left( {\left( {{\sigma _1} - {\sigma _2}} \right)/k} \right)}^2}} \right],k > 0。$

(7)

2.2 模型的水平集方程

采用水平集方法^[18-21], 可以将式(3)中的积分扩展到整个图像区域Ω上。设ϕ:Ω→R为Lipschitz函数, 其零水平集为Γ, 即Γ={(x, y)∈Ω:ϕ(x, y)=0}, 规定曲线内部ϕ < 0, 曲线外部ϕ>0。则引入水平集函数后, 式(3)变为

$ \begin{array}{*{20}{c}} {\mathop {\min }\limits_{\phi \left( {x,y} \right),c_1^{\left( i \right)},c_2^{\left( i \right)}} F\left( {\phi \left( {x,y} \right),c_1^{\left( i \right)},c_2^{\left( i \right)}} \right) = \mu \iint\limits_\mathit{\Omega } {d\left( {x,y} \right){g_{{\text{color}}}}\left| {\nabla H\left( \phi \right)} \right|{\text{d}}x{\text{d}}y} + } \\ {\frac{1}{3}\iint\limits_\mathit{\Omega } {\sum\limits_{i = 1}^3 {\lambda _1^{\left( i \right)}{{\left| {{u^{\left( i \right)}} - c_1^{\left( i \right)}} \right|}^2}H\left( \phi \right){\text{d}}x{\text{d}}y} } + \frac{1}{3}\iint\limits_\mathit{\Omega } {\sum\limits_{i = 1}^3 {\lambda _2^{\left( i \right)}{{\left| {{u^{\left( i \right)}} - c_2^{\left( i \right)}} \right|}^2}\left( {1 - H\left( \phi \right)} \right){\text{d}}x{\text{d}}y} } + } \\ {\nu \left[ {{{\left( {\iint\limits_\mathit{\Omega } {H\left( \phi \right){\text{d}}x{\text{d}}y} - {\mathit{\boldsymbol{A}}_1}} \right)}^2} + {{\left( {\iint\limits_\mathit{\Omega } {\left( {1 - H\left( \phi \right)} \right){\text{d}}x{\text{d}}y} - {\mathit{\boldsymbol{A}}_2}} \right)}^2}} \right],} \end{array} $

(8)

式中:H为一维Heaviside函数; ∫_Ω|∇H(ϕ(x, y))|dxdy为曲线Γ的长度。

由于标准Heaviside函数在原点处不可微, 因此用规则化H_ε代替式(8)中的标准H函数^[22]

$ {H_{1\varepsilon }}\left( x \right) = \left\{ \begin{array}{l} 0,\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;x < - \varepsilon ,\\ \frac{1}{2}\left[ {1 + \frac{x}{\varepsilon } + \frac{1}{{\rm{ \mathsf{ π} }}}\sin \left( {\frac{{{\rm{ \mathsf{ π} }}x}}{\varepsilon }} \right)} \right],\left| x \right| \le \varepsilon ,{H_{2\varepsilon }}\left( x \right) = \frac{1}{2}\left( {1 + \frac{2}{{\rm{ \mathsf{ π} }}}\arctan \left( {\frac{x}{\varepsilon }} \right)} \right)。\\ 1,\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;x > \varepsilon , \end{array} \right. $

(9)

式(8)变为

$ \begin{array}{*{20}{c}} {\mathop {\min }\limits_{\phi \left( {x,y} \right),c_1^{\left( i \right)},c_2^{\left( i \right)}} {F_\varepsilon }\left( {\phi ,c_1^{\left( i \right)},c_2^{\left( i \right)}} \right) = \mu \int\limits_\mathit{\Omega } {d\left( {x,y} \right){g_{{\rm{color}}}}{\delta _\varepsilon }\left( \phi \right)\left| {\nabla \phi } \right|{\rm{d}}x{\rm{d}}y} + }\\ {\frac{1}{3}\int\limits_\mathit{\Omega } {\sum\limits_{i = 1}^3 {\lambda _1^{\left( i \right)}{{\left| {{u^{\left( i \right)}} - c_1^{\left( i \right)}} \right|}^2}H_\varepsilon\left( \phi \right){\rm{d}}x{\rm{d}}y} } + \frac{1}{3}\int\limits_\mathit{\Omega } {\sum\limits_{i = 1}^3 {\lambda _2^{\left( i \right)}{{\left| {{u^{\left( i \right)}} - c_2^{\left( i \right)}} \right|}^2}\left( {1 - H_\varepsilon\left( \phi \right)} \right){\rm{d}}x{\rm{d}}y} } + }\\ {\nu \left[ {{{\left( {\int\limits_\mathit{\Omega } {{H_\varepsilon }\left( \phi \right){\rm{d}}x{\rm{d}}y} - {\mathit{\boldsymbol{A}}_1}} \right)}^2} + {{\left( {\int\limits_\mathit{\Omega } {\left( {1 - {H_\varepsilon }\left( \phi \right)} \right){\rm{d}}x{\rm{d}}y} - {\mathit{\boldsymbol{A}}_2}} \right)}^2}} \right]。} \end{array} $

(10)

保持ϕ(x, y)不变, 式(10)最小化得c₁⁽ⁱ⁾, c₂⁽ⁱ⁾的计算公式如下:

$ c_1^{\left( i \right)} = \frac{{\int\limits_\mathit{\Omega } {{u^{\left( i \right)}}{H_\varepsilon }\left( {\phi \left( {x,y} \right)} \right){\rm{d}}x{\rm{d}}y} }}{{\int\limits_\mathit{\Omega } {{H_\varepsilon }\left( {\phi \left( {x,y} \right)} \right){\rm{d}}x{\rm{d}}y} }},i = 1,2,3, $

(11)

如果$\int\limits_\mathit{\Omega } {{H_\varepsilon }\left( {\phi \left( {x, y} \right)} \right)} {\rm{d}}x{\rm{d}}y > 0 $, 则曲线Γ内部非空;

$ c_2^{\left( i \right)} = \frac{{\int\limits_\mathit{\Omega } {{u^{\left( i \right)}}\left[ {1 - {H_\varepsilon }\left( {\phi \left( {x,y} \right)} \right)} \right]{\rm{d}}x{\rm{d}}y} }}{{\int\limits_\mathit{\Omega } {\left[ {1 - {H_\varepsilon }\left( {\phi \left( {x,y} \right)} \right)} \right]{\rm{d}}x{\rm{d}}y} }},i = 1,2,3; $

(12)

如果$\int\limits_\mathit{\Omega } {\left[{1-{H_\varepsilon }\left( {\phi \left( {x, y} \right)} \right)} \right]} {\rm{d}}x{\rm{d}}y > 0 $, 则曲线Γ外部非空。保持c₁⁽ⁱ⁾, c₂⁽ⁱ⁾不变, 关于ϕ(x, y)最小化式(10)得Euler-Lagrange方程如下:

$ \begin{array}{*{20}{c}} {{\delta _\varepsilon }\left( \phi \right)\mu \nabla \cdot \left( {G\left( {x,y} \right)\frac{{\nabla \phi }}{{\left| {\nabla \phi } \right|}}} \right) - \frac{1}{3}{\delta _\varepsilon }\left( \phi \right)\left( {\sum\limits_{i = 1}^3 {\lambda _1^{\left( i \right)}{{\left( {{u^{\left( i \right)}} - c_1^{\left( i \right)}} \right)}^2}} - \sum\limits_{i = 1}^3 {\lambda _2^{\left( i \right)}{{\left( {{u^{\left( i \right)}} - c_2^{\left( i \right)}} \right)}^2}} } \right) - }\\ {\nu \left[ {\left( {\int\limits_\mathit{\Omega } {{H_\varepsilon }\left( \phi \right){\rm{d}}x{\rm{d}}y} - {A_1}} \right) - \left( {\int\limits_\mathit{\Omega } {\left( {1 - {H_\varepsilon }\left( \phi \right)} \right){\rm{d}}x{\rm{d}}y} - {A_2}} \right)} \right] = 0,} \end{array} $

(13)

边界条件为$ \frac{{\partial \phi }}{{\partial \mathit{\boldsymbol{\vec n}}}} = 0, \mathit{\boldsymbol{\vec n}}$为边界∂Ω的单位外法向量, G(x, y)=d(x, y)g_color。在式(13)中加入气球力^{[16, 20]}αG(x, y)|∇ϕ|, α>0可提高模型收敛速度。式(13)引入时间变量t, 得到下式的梯度下降流公式

$ \begin{array}{*{20}{c}} {\frac{{\partial \phi }}{{\partial t}} = \mu {\delta _\varepsilon }\left( \phi \right)\nabla \cdot \left( {G\left( {x,y} \right)\frac{{\nabla \phi }}{{\left| {\nabla \phi } \right|}}} \right) - \frac{1}{3}{\delta _\varepsilon }\left( \phi \right)\left( {\sum\limits_{i = 1}^3 {\lambda _1^{\left( i \right)}{{\left( {{u^{\left( i \right)}} - c_1^{\left( i \right)}} \right)}^2}} } - \sum\limits_{i = 1}^3 {\lambda _2^{\left( i \right)}{{\left( {{u^{\left( i \right)}} - c_2^{\left( i \right)}} \right)}^2}}\right) - }\\ {\nu \left[ {\left( {\int\limits_\mathit{\Omega } {{H_\varepsilon }\left( \phi \right){\rm{d}}x{\rm{d}}y} - {A_1}} \right) - \left( {\int\limits_\mathit{\Omega } {\left( {1 - {H_\varepsilon }\left( \phi \right)} \right){\rm{d}}x{\rm{d}}y} - {A_2}} \right)} \right] + \alpha G\left( {x,y} \right)\left| {\nabla \phi } \right|,} \end{array} $

(14)

满足纽曼边界条件。

2.3 加法算子分裂算法

为了降低计算的复杂程度, 同时避免过度分割的问题^[23], 式(14)可简化为:

$ \frac{{\partial \phi }}{{\partial t}} = \mu {\delta _\varepsilon }\left( \phi \right)\nabla \cdot \left( {F\nabla \phi } \right) + \varphi = \mu {\delta _\varepsilon }\left( \phi \right)\left( {{\partial _x}\left( {F{\partial _x}\phi } \right) + {\partial _y}\left( {F{\partial _y}\phi } \right)} \right) + \varphi , $

(15)

式中:

$ \begin{array}{*{20}{c}} {F = \frac{G}{{{{\left| {\nabla \phi } \right|}_\beta }}},\varphi = - \frac{1}{3}{\delta _\varepsilon }\left( \phi \right)\left( {\sum\limits_{i = 1}^3 {\lambda _1^{\left( i \right)}{{\left( {{u^{\left( i \right)}} - c_1^{\left( i \right)}} \right)}^2}} - \sum\limits_{i = 1}^3 {\lambda _2^{\left( i \right)}{{\left( {{u^{\left( i \right)}} - c_2^{\left( i \right)}} \right)}^2}} } \right) - }\\ {\nu \left[ {\left( {\int\limits_\mathit{\Omega } {{H_\varepsilon }\left( \phi \right){\rm{d}}x{\rm{d}}y} - {\mathit{\boldsymbol{A}}_1}} \right) - \left( {\int\limits_\mathit{\Omega } {\left( {1 - {H_\varepsilon }\left( \phi \right)} \right){\rm{d}}x{\rm{d}}y} - {\mathit{\boldsymbol{A}}_2}} \right)} \right] + \alpha G\left( {x,y} \right)\left| {\nabla \phi } \right|。} \end{array} $

(16)

加法算子分裂算法^[24-25]是将一个m维的空间算子分裂成m个一维算子的和。本研究是将二维空间算子分裂成两个一维问题。下面以x方向为例

$ \begin{array}{*{20}{c}} {\frac{{\phi _{i,j}^{k + 1} - \phi _{i,j}^k}}{{\Delta t}} = \mu {\delta _\varepsilon }\left( \phi \right)\frac{{F_{i + 1/2,j}^k{{\left( {{\partial _x}{\phi ^{k + 1}}} \right)}_{i + 1/2,j}} - F_{i - 1/2,j}^k{{\left( {{\partial _x}{\phi ^{k + 1}}} \right)}_{i - 1/2,j}}}}{{{h_x}}} + {\varphi _{i,j}} = }\\ {\mu {\delta _\varepsilon }\left( \phi \right)\left( \begin{array}{l} \frac{{F_{i + 1,j}^k + F_{i,j}^k}}{{2h_x^2}}\left( {\phi _{i + 1,j}^{k + 1} - \varphi _{i,j}^{k + 1}} \right)\\ - \frac{{F_{i,j}^k + F_{i - 1,j}^k}}{{2h_x^2}}\left( {\phi _{i,j}^{k + 1} - \phi _{i - 1,j}^{k + 1}} \right) \end{array} \right) + {\varphi _{i,j}},} \end{array} $

(17)

即

$ \phi _{i,j}^{k + 1} = \phi _{i,j}^k + \Delta t\left( {{l_1}\phi _{i + 1,j}^{k + 1} - {l_2}\phi _{i,j}^{k + 1} + {l_3}\phi _{i - 1,j}^{k + 1}} \right) + \Delta t{\varphi _{i,j}}, $

(18)

式中:

$ {l_1} = \mu {\delta _\varepsilon }\left( {\phi _{i,j}^k} \right)\frac{{F_{i + 1,j}^k + F_{i,j}^k}}{{2h_x^2}};{l_2} = \mu {\delta _\varepsilon }\left( {\phi _{i,j}^k} \right)\frac{{F_{i + 1,j}^k + 2F_{i,j}^k + F_{i - 1,j}^k}}{{2h_x^2}};{l_3} = \mu {\delta _\varepsilon }\left( {\phi _{i,j}^k} \right)\frac{{F_{i,j}^k + F_{i - 1,j}^k}}{{2h_x^2}}。$

同理, 可得y方向递推方程

$ \phi _{i,j}^{k + 1} = \phi _{i,j}^k + \Delta t\left( {{{l'}_1}\phi _{i,j + 1}^{k + 1} - {{l'}_2}\phi _{i,j}^{k + 1} + {{l'}_3}\phi _{i,j - 1}^{k + 1}} \right) + \Delta t{\varphi _{i,j}}, $

(19)

式中: $ l{'_1} = \mu {\delta _\varepsilon }\left( {{\phi _{i, j}}} \right)\frac{{F_{i, j + 1}^k + F_{i, j}^k}}{{2h_y^2}};{l_1} = \mu {\delta _\varepsilon }\left( {{\phi _{i, j}}} \right)$$\frac{{F_{i, j + 1}^k + 2F_{i, j}^k + F_{i, j-1}^k}}{{2h_y^2}};l{'_1} = \mu {\delta _\varepsilon }\left( {{\phi _{i, j}}} \right)\frac{{F_{i, j}^k + F_{i, j-1}^k}}{{2h_y^2}} $。

本研究取h_x=h_y=1。因此, 式(18)(19)可合并改写为下列所示的矩阵形式

$ \left( {\mathit{\boldsymbol{I}} - 2\Delta t{\mathit{\boldsymbol{A}}_m}\left( {{\phi ^k}} \right)} \right)\phi _m^{k + 1} = {{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over \varphi } }^k},\;\;\left( {m = 1,2} \right),\\{\phi ^{k + 1}} = \frac{1}{2}\sum\limits_{m = 1}^2 {\phi _m^{k + 1}} ,\;\;k = 0,1,2, \cdots ,{{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over \varphi } }^k} = {\phi ^k} + \Delta t{\varphi ^k}。$

(20)

式中:I为单位矩阵; ϕ₁^k+1, ϕ₂^k+1分别为x方向和y方向更新的ϕ值; A_m(m=1, 2)为两个三角矩阵。

3 试验仿真结果

本研究分别对合成图像和真实图像进行了分割试验, 以验证本研究提出方法的有效性。由于合成图像的边界非常清晰, 因此在合成图像分割时, 均采用函数H_1ε; 在真实图像的分割中, 根据图像差异采用H_1ε或H_2ε。在进行图像分割时, 初始轮廓均根据初始约束点的情况初始化为圆:

$ \phi \left( {x,y} \right) = r - c - \sqrt {{{\left( {x - {x_0}} \right)}^2} + {{\left( {y - {y_0}} \right)}^2}} ,c \ge 0, $

式中: $ r = \min \left( {\left\| {p-{p_0}} \right\|} \right);{p_0}\left( {{x_0}, {y_0}} \right);$$ {x_0} = \sum\limits_{i = 1}^n {{x_i}/n} ;y_0=\sum\limits_{i = 1}^n {{y_i}/n} ;p \in \mathit{\Omega ;n}$为约束点个数。

(1) 合成图像选择性分割

首先将本研究应用到合成图像的分割中, 其中一幅图像为物体彼此接触, 另一幅图像则为物体间有重叠区域, 两幅图像中的物体颜色差异不大。

同一目标不同位置约束点分割效果如图 1所示。原始图像为2个圆和1个正方形, 3个物体的颜色非常接近, 并且紧挨着。其中(a)和(c)的区别在于约束点的选取, 导致初始曲线的半径和位置均不相同, 选取ε=2, μ=100, k=4, λ₁⁽ⁱ⁾=λ₂⁽ⁱ⁾=1, c=7分别对(a)和(c)进行分割, (b)和(d)为其对应分割结果, 由图 1可看出, 边界分割效果良好。由此可见, 在模型中加入区域信息后, 可以克服基于边缘分割模型对初始曲线位置敏感的缺点。另外由(a)和(c)可以看出, 无论约束点是否在目标物体边界上, 本研究方法仍然能得到理想的分割效果, 但前提是必须保证不能将约束点取在临近的其他物体内部。

不同目标分割效果如图 2所示。图像为颜色相差不大的一个正方形和一个三角形, 并且两个物体有重叠部分。(a)和(c)分别表示目标物体为三角形和矩形, 选取α=0.009, Δt=25, σ=4, ε=6, μ=12, k=9, λ₁⁽ⁱ⁾=λ₂⁽ⁱ⁾=1, 对(a)进行分割得到(b); 选取α=0.003, Δt=23, σ=4, ε=2, μ=100, k=4, c=13, λ₁⁽ⁱ⁾=λ₂⁽ⁱ⁾=1, 对(c)进行分割得到(d), 综合结果可以看出, 本研究提出的方法能够对有重叠区域的图像分别进行正确的选择性分割, 且分割边界平均精度较高。

(2) 真实图像的选择性分割

图 3为枫叶采集图像分割效果, 目标物体为3片真实枫叶中的一片叶子, 其形状较为规则。图(a)给出任意3个约束点, 选取α=0.002, Δt=25, σ=4, ε=7, μ=100, k=15, c=4, λ₁⁽ⁱ⁾=λ₂⁽ⁱ⁾=1, 对(a)进行选择性分割, 通过仿真实验发现其分割效果(b)能较好地对感兴趣区域进行选择性分割。

图 4为来自医学研究的肾脏的CT图像, 图片结构较为复杂, 需分别对左右两颗肾脏进行分割。其中(a)和(c)分别是左右两颗肾脏靠近边缘处的3个约束点及初始水平集曲线。选取α=0.004, Δt=30, σ=4, ε=6, μ=100, k=3, c=4, λ₁⁽ⁱ⁾=λ₂⁽ⁱ⁾=1, 得到(b)为左边肾脏分割效果图, 结果较理想; 选取α=0.007, Δt=30, σ=2, ε=3, μ=100, k=5, c=4, λ₁⁽ⁱ⁾=λ₂⁽ⁱ⁾=1, 得到(d)为右边肾脏分割效果图, 结果同样较为理想。由此可见, 本研究提出的方法能对此类背景较为复杂的目标物体进行正确分割, 并且取得较高的边界分割精度。

图 5为来自医学图像脑积水CT图像, 分别对脑积水CT图像的左边和右边进行选择性分割。这类图像边缘有一定的模糊地带, 其中(a)和(c)分别为左右脑积水及位于边缘处的4个约束点确定的初始水平集曲线, 选取α=0.009, Δt=30, σ=5, ε=6, μ=200, k=10, c=3, λ₁⁽ⁱ⁾=λ₂⁽ⁱ⁾=1时, 得到(b)为左边部分脑积水的分割效果; 选取α=0.009, Δt=25, σ=3, ε=7, μ=200, k=12, λ₁⁽ⁱ⁾=λ₂⁽ⁱ⁾=1, c=3时, 得到(d)为右边部分脑积水分割效果。依据分割效果可知:虽然该图的目标物边缘存在模糊地带, 但是本研究提出的方法仍能得到较高精度的分割结果。

图 6为医学检验的红细胞图像, (a)标记初始选择点处为待分割的目标物体, 即细胞核周围给定5个约束点, 当α=0.003, Δt=30, σ=4, ε=2, μ=1, k=2, c=4, λ₁⁽ⁱ⁾=λ₂⁽ⁱ⁾=1时, 得到如(b)所示的分割结果, 边界精度较理想。

通过上述几类图片分别进行分割仿真试验可知:在待分割目标物体上选取初始点, 若待分割部分边界较为简单且清晰, 3个初始点基本上能满足要求, 位置可以在目标物体上任意选取; 若待分割物体边界模糊或轮廓复杂, 初始点数量一般扩充为4~6个, 位置尽可能的选取靠近目标物体边缘处, 这样获得的分割边界精度将更为理想。

4 结论

本研究基于Lavdie-Chen的灰度图像单水平集选择性分割方法, 提出带几何约束的彩色图像选择性分割方法。该方法将彩色图像作为整体得到图像梯度和边缘检测函数; 在边缘检测函数中引入可以控制分割速度的参数; 图像边缘梯度、距离函数以及约束点形成的多边形内外面积共同控制着曲线演化的走向; 水平集的引入将曲线上的积分拓展到图像平面上; 加法算子分裂算法的运用加快了分割过程。针对合成图像与真实图像分别进行分割仿真试验, 依次选取不同位置、不同目标、不同清晰程度的待分割物体进行试验, 结果发现:本研究提出的方法对边缘不清晰、灰度不均匀的图像进行选择性分割, 其获取的边界精度均较为理想, 体现了较好的鲁棒性。

本研究的创新点:(1)距离函数与彩色图像边缘函数的充分结合, 共同控制曲线的演化过程, 增强曲线演化的稳定性; (2)新的分割方法有效地采用了正则化优化方案, 克服了原有模型在凹陷处分割效果不理想的缺点, 实现了演化曲线的自动收缩与扩张, 控制着曲线演化的方向性; (3)快速算法的引进降低了计算的复杂度, 加快了迭代运算的速度, 加强了曲线朝着目标边缘演化的快速性。