彩色纹理图像分割在医学、公安、交通等诸多领域存在大量应用。由于其纹理成分复杂，纹理强度不同，采用传统的变分方法基于边缘^[1]、区域^[2-3]、先验形状^[4]等特征的图像分割模型已很难得出准确分割结果。当彩色图像中包含纹理成分时，其纹理特征表达和区域边界描述仍存在诸多困难，如不充分考虑各层图像之间的耦合效应，则将引起边缘以及纹理模糊。

近些年，有许多针对纹理图像分割的研究，文献[5-6]分别通过Gabor、Wavelet滤波器生成多尺度矢量图像，采用矢量图像分割的Chan-Vese模型^[7]实现纹理图像分割，Li^[8]在Gabor变换基础上，用四元数处理多通道信息。滤波器会减少图像中的成分，导致纹理图像区域边缘模糊，分割不准确；Wang等^[9]通过图像分解变分方法^[10]提取纹理图像中的分段常值结构化图像成分，并将图像分解模型与Chan-Vese模型组合实现纹理图像分割。这种模型效率低下且分割效果不佳。Gillboa等^[11]提出了非局部算子在图像处理中的应用，Bresson等^[12]将传统的Mumford-Shah模型^[13]与非局部算子结合为非局部Mumford-Shah模型，但忽略了纹理图像成分隐藏的图像结构化成分，导致目标边缘分割较差。文献[14]通过扩展上述非局部模型为多通道模型并增加Tikhonov规则项，提出了Mumford-Shah-MT(mumford- shah-multichannel tikhonov)模型，但由于彩色图像各层图像边缘处的扩散强度不一致，导致轮廓线无法准确贴近目标边缘。在图像恢复中表现良好的多通道全变差模型^[15]实现了良好的彩色图像层与层之间的耦合效果且边缘不扩散，在非局部Mumford-Shah模型分割中增加多通道全变差规则项可以实现图像边缘分割更加精确。

由上述分析可知，综合非局部Mumford-Shah模型与多通道全变差模型的优点是提出非局部Mumford-Shah多通道全变差模型，该模型可以对纹理成分准确描述，并保持彩色图像层与层之间耦合性，最终实现彩色纹理图像更精确的分割。为了提高计算效率，笔者还设计了所提出模型的ADMM(alternating direction method of multipliers)算法^[16-18]。

1 彩色纹理图像分割模型 1.1 非局部Mumford-Shah模型及其拓展

Mumford-Shah模型^[13]是变分图像分割的基础，其基本思想是用分段光滑图像及其边界近似来分割图像。为了克服在不同维度上优化的困难，Bresson等^[18]直接采用二值标记函数ϕ(x)∈(0, 1)划分区域，用总变差表达分割线长度。提出与Chan-Vese模型^[3]对应的分段光滑图像分割的变分模型为

$ \begin{array}{l} \mathop {\min }\limits_{{u_1},{u_2},\phi } \left\{ {E\left( {{u_1},{u_2},\phi } \right) = {\alpha _1}\int_{{\mathit{\Omega }_1}} {\left( {{{\left( {{u_1} - f} \right)}^2} + {{\left| {\nabla {u_1}} \right|}^2}} \right)\phi {\rm{d}}x} + } \right.\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{\alpha _2}\int_{{\mathit{\Omega }_2}} {\left( {{{\left( {{u_2} - f} \right)}^2} + {{\left| {\nabla {u_2}} \right|}^2}} \right)\left( {1 - \phi } \right){\rm{d}}x} + \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\left. {\gamma \int_\mathit{\Omega } {\left| {\nabla \phi } \right|{\rm{d}}x} } \right\}, \end{array} $

(1)

其中:u₁, u₂分别表示分段常值图像在区域Ω₁和Ω₂的均值; f为原图像。

$ \phi \left( x \right) = \left\{ \begin{array}{l} 1,x \in {\mathit{\Omega }_1}\\ 0,x \in {\mathit{\Omega }_2} \end{array} \right.,\mathit{\Omega } = {\mathit{\Omega }_1} \cup {\mathit{\Omega }_2},{\mathit{\Omega }_1} \cap {\mathit{\Omega }_2} = \emptyset , $

为了使演化方程在求解时得到稳态解，首先将离散的ϕ(x)凸松弛为ϕ(x)∈[0, 1]，最终通过阈值化方法还原为原问题的离散解。

借助纹理图像表达的非局部算子，Bresson，Chan^[12]将式(1)扩展为非局部Mumford-Shah模型

$ \begin{array}{l} \mathop {\min }\limits_{{u_1},{u_2},\phi } \left\{ {E\left( {{u_1},{u_2},\phi } \right) = {\alpha _1}\int_{{\mathit{\Omega }_1}} {\left( {{{\left( {{u_1} - f} \right)}^2} + {\beta _1}{{\left| {{\nabla _{{\rm{NL}}}}{u_1}} \right|}^2}} \right)\phi {\rm{d}}x} + } \right.\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{\alpha _2}\int_{{\mathit{\Omega }_2}} {\left( {{{\left( {{u_2} - f} \right)}^2} + {\beta _2}{{\left| {{\nabla _{{\rm{NL}}}}{u_2}} \right|}^2}} \right)\left( {1 - \phi } \right){\rm{d}}x} + \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\left. {\gamma \int_\mathit{\Omega } {\left| {\nabla \phi } \right|{\rm{d}}x} } \right\}。\end{array} $

(2)

为了使上述模型适用彩色纹理图像分割并发挥纹理图像成分隐藏的结构化图像成分在图像分割中的作用，Wang等^[14]提出了非局部Mumford-Shah-MT模型，即

$ \begin{array}{l} \mathop {\min }\limits_{u,\phi \in \left\{ {0,1} \right\}} \left\{ {E\left( {{u_{ij}},\phi } \right) = \sum\limits_{i = 1}^2 {\sum\limits_{j = 1}^3 {{\alpha _i}\int_\mathit{\Omega } {\left( {{{\left( {{u_{ij}} - {f_j}} \right)}^2} + {\beta _i}{{\left| {\nabla {u_{ij}}} \right|}^2} + {\lambda _i}{{\left| {{\nabla _{{\rm{NL}}}}{u_{ij}}} \right|}^2}} \right){\chi _i}\left( \phi \right){\rm{d}}x} } } + } \right.\\ \left. {\gamma \int_\mathit{\Omega } {\left| {\nabla \phi } \right|{\rm{d}}x} } \right\}, \end{array} $

(3)

其中：χ₁(ϕ)=ϕ, χ₂(ϕ)=1－ϕ分别为分割目标部分与背景部分的特征函数; f(x)=(f₁(x), f₂(x), f₃(x)):Ω∈R³，表示彩色纹理图像，f_j(x)表示彩色图像的3个通道。

$ \sum\limits_{j = 1}^3 {|\nabla {u_j}{|^2}} $虽然考虑到图像中隐藏的结构化成分但并没有考虑到彩色图像层与层之间的耦合，导致边缘以及纹理模糊，影响分割精度。笔者采用考虑层与层之间耦合关系的多通道全变差规则项$ \sqrt {\sum\limits_{j = 1}^3 {|\nabla {u_j}{|^2}} } $替代$ \sum\limits_{j = 1}^3 {|\nabla {u_j}{|^2}} $，以改进分割精度。

1.2 非局部算子

图像科学中的非局部导数、非局部散度等定义由Gillboa等^[11]基于Buades等^[19]的非局部均值概念提出，并在图像恢复^[20-22]与分割^[23-25]等领域表现较好的效果。该方法定义了表示图像片之间相似度的权函数为

$ \mathit{\boldsymbol{w}}\left( {x,y} \right) = {{\rm{e}}^{\frac{{ - G\sigma \times \left\| {f\left( {x + \cdot } \right) - f\left( {y + \cdot } \right)} \right\|_2^2}}{h}}}, $

(4)

其中：x∈Ω，y∈Ω；G_σ是高斯核函数；σ是高斯核函数的标准差；h是控制相似度的阈值；f(x+·)表示以x为中心的小邻域(片)中的像素值；f(y+·)表示以y为中心的小邻域(片)的像素值。为提高计算效率，上述定义可限于搜索窗口S(x)中，如图 1所示。x, y∈S(x)表示搜索窗中的任一像素点。权w(x, y)是基于以点x、y间为中心(片)的广义距离，而非仅2点的局部信息。

基于权w(x, y)可定义x，y点间的非局部导数为

$ {\partial _x}u\left( {x,y} \right) = \left( {u\left( y \right) - u\left( x \right)} \right)\sqrt {\mathit{\boldsymbol{w}}\left( {x,y} \right)} ,x,y \in \mathit{\Omega }, $

(5)

以x为基点的非局部梯度为

$ {\nabla _{{\rm{NL}}}}u\left( {x,y} \right) = \left( {u\left( y \right) - u\left( x \right)} \right)\sqrt {\mathit{\boldsymbol{w}}\left( {x,y} \right)} ,x,y \in \mathit{\Omega }。$

(6)

非局部梯度表示图像中片与片的映射: Ω×Ω→R，非传统意义下的矢量。对于非局部矢量v(x, y): Ω×Ω→R，定义相关非局部点积

$ \left\langle {{\mathit{\boldsymbol{v}}_1} \cdot {\mathit{\boldsymbol{v}}_2}} \right\rangle \left( x \right) = \int_\mathit{\Omega } {{\mathit{\boldsymbol{v}}_1}\left( {x,y} \right){\mathit{\boldsymbol{v}}_2}\left( {x,y} \right){\rm{d}}y} :\mathit{\Omega } \to \mathit{R}。$

(7)

相应非局部矢量的模为

$ \left| \mathit{\boldsymbol{v}} \right|\left( x \right) = \sqrt {\int_\mathit{\Omega } {{{\left( {\mathit{\boldsymbol{v}}\left( {x,y} \right)} \right)}^2}{\rm{d}}y} } :\mathit{\Omega } \to \mathit{R}。$

(8)

根据广义散度原理

$ \left\langle {{\nabla _{{\rm{NL}}}}u,\mathit{\boldsymbol{v}}} \right\rangle = - \left\langle {u,{\rm{di}}{{\rm{v}}_{{\rm{NL}}}}\mathit{\boldsymbol{v}}} \right\rangle 。$

(9)

可求出非局部散度

$ \left\langle {{\nabla _{{\rm{NL}}}} \cdot \mathit{\boldsymbol{v}}} \right\rangle \left( x \right) = \int_\mathit{\Omega } {\left( {\mathit{\boldsymbol{v}}\left( {x,y} \right) - \mathit{\boldsymbol{v}}\left( {y,x} \right)} \right)\sqrt {\mathit{\boldsymbol{w}}\left( {x,y} \right)} {\rm{d}}y} :\mathit{\Omega } \to \mathit{R}。$

(10)

和非局部拉普拉斯算子

$ {\nabla _{{\rm{NL}}}}u\left( x \right) = \frac{1}{2}\left( {{\nabla _{{\rm{NL}}}} \cdot \left( {{\nabla _{{\rm{NL}}}}u} \right)} \right)\left( x \right) = \int_\mathit{\Omega } {\left( {u\left( y \right) - u\left( x \right)} \right)\mathit{\boldsymbol{w}}\left( {x,y} \right){\rm{d}}y} 。$

(11)

基于上述定义，非局部梯度的模可表达为

$ \left| {{\nabla _{{\rm{NL}}}}u} \right|\left( x \right) = \sqrt {\int_\mathit{\Omega } {{{\left( {u\left( y \right) - u\left( x \right)} \right)}^2}\mathit{\boldsymbol{w}}\left( {x,y} \right){\rm{d}}y} } :\mathit{\Omega } \to \mathit{R}。$

(12)

2 非局部Mumford-Shah多通道全变差模型及其ADMM算法

利用纹理图像中隐藏的结构化图像成分，考虑到层与层之间的耦合，用多通道全变差规则项$ \sqrt {\sum\limits_{j = 1}^3 {|\nabla {u_j}{|^2}} } $替代非局部Mumford-Shah-MT模型中Multi-channel Tikhonov规则项$ \sum\limits_{j = 1}^3 {|\nabla {u_j}{|^2}} $，得到如下非局部Mumford-Shah多通道全变差模型

$ \begin{array}{*{20}{c}} {\mathop {\min }\limits_{{u_i},\phi \in \left[ {0,1} \right]} \left\{ {E\left( {{u_i},\phi } \right) = \sum\limits_{i = 1}^2 {\sum\limits_{j = 1}^3 {{\alpha _i}\int_\mathit{\Omega } {\left( {{{\left( {{u_{ij}} - {f_j}} \right)}^2} + {\mu _{ij}}{{\left| {{\nabla _{{\rm{NL}}}}{u_{ij}}} \right|}^2}} \right){\chi _i}\left( \phi \right){\rm{d}}x} } } + } \right.}\\ {\left. {\sum\limits_{i = 1}^2 {{\beta _i}\int_\mathit{\Omega } {\sqrt {\sum\limits_{j = 1}^3 {{{\left| {\nabla {u_{ij}}} \right|}^2}} } {\chi _i}\left( \phi \right){\rm{d}}x} } + \gamma \int_\mathit{\Omega } {\left| {\nabla \phi } \right|{\rm{d}}x} } \right\}。} \end{array} $

(13)

等式右端第1项表示数据项，表达分段光滑图像与原图像的逼近程度，(i=2，3；j=1，2，3)，表示彩色图像3个通道在目标和背景2个区域的分段常值，是惩罚参数，值越大，越接近；第2项使得到的图像为光滑的，惩罚参数越大，图像越光滑；第3项表示彩色纹理图像隐藏的结构化图像成分的梯度，惩罚参数控制着结构化成分的光滑程度；第4项表示分段光滑图像边界的长度，惩罚参数对轮廓线的紧凑程度惩罚。当β_i=0时，该模型是多通道非局部Mumford-Shah模型^[11]；μ_i=0时，该模型是分段光滑彩色图像分割的Chan-Vese模型^[3]；当β_i=0、μ_i=0时，该模型为分段常值图像分割的经典Chan-Vese模型^[2]。

所提出模型与模型(3)对比，从图像边缘保持和图像之间层与层的耦合性具有优越性：全变差规则项在图像边缘保持方面优于Thiknov规则项，分割轮廓线会紧贴目标边缘，过分割和欠分割现象较少。模型(3)中的规则项并没有考虑彩色图像3个通道的耦合，先分别计算3个通道，最后叠加。由文献[15]可知，$ \sqrt {\sum\limits_{j = 1}^3 {|\nabla {u_j}{|^2}} } $考虑层与层之间的耦合，纹理保持较好，不会出现过分割。

当直接采用交替优化方法求解时会得到关于u₁，u₂和ϕ3个变量的曲率，导致差分计算格式复杂。为了提高数值计算效率，通过ADMM ^[16-17]算法，引入辅助变量、Lagrange乘子以及惩罚参数，在交替优化的每一步将原优化问题分解为一系列简单、高效的子优化问题。

研究通过引入辅助变量v_ij=▽u_ij, (i=1, 2;j=1, 2, 3)、w=▽ϕ、Lagrange乘子λ_ij, (i=1, 2;j=1, 2, 3)、σ及惩罚参数θ_i, (i=1, 2)、ω，将式(11)转化为交替优化能量模型

$ \begin{array}{*{20}{c}} {\left( {u_{ij}^{k + 1},{\mathit{\boldsymbol{v}}^{k + 1}}_{ij},{\phi ^{k + 1}},{\mathit{\boldsymbol{w}}^{k + 1}}} \right) = }\\ {\mathop {arg\min }\limits_{{u_{ij}},{\mathit{\boldsymbol{v}}_{ij}},\mathit{\boldsymbol{w}},\phi \in \left[ {0,1} \right]} \left\{ {E\left( {{u_{ij}},{\mathit{\boldsymbol{v}}_{ij}},\mathit{\boldsymbol{w}},\phi } \right) = \sum\limits_{i = 1}^2 {\sum\limits_{j = 1}^3 {{\alpha _i}\int_\mathit{\Omega } {{Q_{ij}}{\chi _i}\left( \phi \right){\rm{d}}x} } } + \sum\limits_{i = 1}^2 {{\beta _i}\int_\mathit{\Omega } {\sqrt {\sum\limits_{j = 1}^3 {{{\left| {{\mathit{\boldsymbol{v}}_{ij}}} \right|}^2}} } {\chi _i}\left( \phi \right){\rm{d}}x} } + } \right.}\\ {\sum\limits_{i = 1}^2 {\sum\limits_{j = 1}^3 {\frac{{{\theta _i}}}{2}\int_\mathit{\Omega } {{{\left| {{\mathit{\boldsymbol{v}}_{ij}} - \nabla {u_{ij}}} \right|}^2}{\rm{d}}x} } } + \sum\limits_{i = 1}^2 {\sum\limits_{j = 1}^3 {\int_\mathit{\Omega } {\lambda _{ij}^k \cdot \left( {{\mathit{\boldsymbol{v}}_{ij}} - \nabla {u_{ij}}} \right){\rm{d}}x} } } + \frac{\mathit{\boldsymbol{\omega }}}{2}\int_\mathit{\Omega } {{{\left| {\mathit{\boldsymbol{w}} - \nabla \phi } \right|}^2}{\rm{d}}x} + \gamma \int_\mathit{\Omega } {\left| {\nabla \phi } \right|{\rm{d}}x} 。} \end{array} $

(14)

$ \left\{ \begin{array}{l} \lambda _{ij}^{k + 1} = \lambda _{ij}^k + {\theta _i}\left( {{\mathit{\boldsymbol{v}}^{k + 1}}_{ij} - \nabla u_{ij}^{k + 1}} \right)\\ {\mathit{\boldsymbol{\sigma }}^{k + 1}} = {\mathit{\boldsymbol{\sigma }}^k} + \mathit{\boldsymbol{\omega }}\left( {{\mathit{\boldsymbol{w}}^{k + 1}} - \nabla {\phi ^{k + 1}}} \right) \end{array} \right., $

(15)

上述模型中:Q_ij=(u_ij－f_j)²+μ_i|▽_NLu_ij|²。

ADMM算法的基本思想是在每一步交替优化时，暂固定其他变量不变，仅对某一变量求解极值，上述优化可分解为如下子优化问题

$ \left( {u_{ij}^{k + 1}} \right) = \arg \mathop {\min }\limits_{{u_{ij}}} E\left( {{u_{ij}},\mathit{\boldsymbol{v}}_{ij}^k,{\mathit{\boldsymbol{w}}^k},{\phi ^k}} \right), $

(16)

$ \left( {\mathit{\boldsymbol{v}}_{ij}^{k + 1}} \right) = \arg \mathop {\min }\limits_{{\mathit{\boldsymbol{v}}_{ij}}} E\left( {u_{ij}^{k + 1},{\mathit{\boldsymbol{v}}_{ij}},{\mathit{\boldsymbol{w}}^k},{\phi ^k}} \right), $

(17)

$ \left( {{\phi ^{k + 1}}} \right) = \arg \mathop {\min }\limits_{\phi \in \left[ {0,1} \right]} E\left( {u_{ij}^{k + 1},{\mathit{\boldsymbol{v}}^{k + 1}}_{ij},{\mathit{\boldsymbol{w}}^k},\phi } \right), $

(18)

$ \left( {{\mathit{\boldsymbol{w}}^{k + 1}}} \right) = \arg \mathop {\min }\limits_\mathit{\boldsymbol{w}} E\left( {u_{ij}^{k + 1},{\mathit{\boldsymbol{v}}^{k + 1}}_{ij},\mathit{\boldsymbol{w}},{\phi ^{k + 1}}} \right), $

(19)

对式(16)采用变分方法得到关于u_ij(i=1, 2;j=1, 2, 3)的Euler-Lagrange方程为

$ \left\{ \begin{array}{l} \left( {{u_{ij}} - {f_j}} \right){\chi _i}\left( \phi \right) + \nabla \cdot \left( {\lambda _{ij}^k - {\theta _i}\nabla {u_{ij}} + {\theta _i}\mathit{\boldsymbol{v}}_{ij}^k} \right) + \\ {\mu _i}{\nabla _{{\rm{NL}}}} \cdot \left( {{\chi _i}\left( \phi \right){\nabla _{{\rm{NL}}}}{u_{ij}}} \right) = 0,\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;x \in \mathit{\Omega },\\ - \left( {\lambda _{ij}^k - {\theta _i}\nabla {u_{ij}} + {\theta _i}\mathit{\boldsymbol{v}}_{ij}^k} \right) \cdot \mathit{\boldsymbol{n}} = 0,\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;x \in \mathit{\Omega }, \end{array} \right. $

(20)

该方程通过采用Gauss-Seidel迭代法求解u_ij^k+1(i=1, 2;j=1, 2, 3)。

同样，对于v_ij(i=1, 2;j=1, 2, 3)的求解极值问题式(17)采用变分法，其解可表达为如下解析形式的广义软阈值公式

$ \begin{array}{*{20}{c}} {{v^{k + 1}}_{ij} = \max \left( {\left| {\mathit{\boldsymbol{B}}_{ij}^{k + 1}} \right| - \frac{{{\beta _i}}}{{{\theta _i}}} \cdot \frac{{\int_\mathit{\Omega } {\left| {{v_{ij}}^k} \right|{\rm{d}}x} }}{{\int_\mathit{\Omega } {\sqrt {\sum\limits_{j = 1}^3 {{{\left| {{v_{ij}}} \right|}^2}} } } }},0} \right)\frac{{\mathit{\boldsymbol{B}}_{ij}^{k + 1}}}{{\left| {\mathit{\boldsymbol{B}}_{ij}^{k + 1}} \right|}},}\\ {0\frac{\mathit{\boldsymbol{0}}}{{\left| \mathit{\boldsymbol{0}} \right|}} = \mathit{\boldsymbol{0}},} \end{array} $

(21)

式中，$ {\mathit{\boldsymbol{B}}}_{ij}^{k + 1} = \nabla u_{ij}^{k + 1} - \frac{{\lambda _{ij}^k}}{{{\theta _i}}} $。

类似式(18)和式(19)的解ϕ、w分别对应如下Euler-Lagrange方程

$ \left\{ \begin{array}{l} \left( {\sum\limits_{i = 1}^2 {\sum\limits_{j = 1}^3 {{Q_{ij}}\left( u \right)} } + \sum\limits_{i = 1}^2 {{\beta _i}\int_\mathit{\Omega } {\sqrt {\sum\limits_{j = 1}^3 {{{\left| {\nabla {u_{ij}}} \right|}^2}} } } } } \right)\frac{{\partial {\chi _i}\left( \phi \right)}}{{\partial \phi }},\\ \nabla \cdot \left( {{\mathit{\boldsymbol{w}}^k} - \nabla \phi + {\sigma ^k}} \right) = 0,\;\;\;\;\;\;\;\;\;x \in \mathit{\Omega },\\ - \left( {{\mathit{\boldsymbol{w}}^k} - \nabla \phi + {\sigma ^k}} \right) \cdot \mathit{\boldsymbol{n}} = 0,\;\;\;\;\;\;\;x \in \mathit{\Omega }, \end{array} \right. $

(22)

$ \phi = \max \left( {\min \left( {\phi ,1} \right),0} \right)。$

(23)

及广义软阈值公式

$ \begin{array}{*{20}{c}} {\mathit{\boldsymbol{w}}k + 1 = \max \left( {\left| {\nabla {\phi ^{k + 1}} - \frac{{{\mathit{\boldsymbol{\sigma }}^k}}}{\mathit{\boldsymbol{\omega }}}} \right| - \frac{\gamma }{\mathit{\boldsymbol{\omega }}},0} \right) \cdot \frac{{\nabla {\phi ^{k + 1}} - \frac{{{\mathit{\boldsymbol{\sigma }}^k}}}{\mathit{\boldsymbol{\omega }}}}}{{\left| {\nabla {\phi ^{k + 1}} - \frac{{{\mathit{\boldsymbol{\sigma }}^k}}}{\mathit{\boldsymbol{\omega }}}} \right|}},}\\ {0\frac{\mathit{\boldsymbol{0}}}{{\left| \mathit{\boldsymbol{0}} \right|}} = \mathit{\boldsymbol{0}},} \end{array} $

(24)

离散方程亦通过Gauss-Seidel迭代求解。当满足

$ \frac{{\left| {{E^{k + 1}} - {E^k}} \right|}}{{{E^k}}} \le \varepsilon , $

(25)

上述交替优化过程结束。其中，E^k为第k步能量泛函的值，ε为迭代误差容限。为准确表达分割区域边界，最终需采用硬阈值公式将凸松弛后的标记函数恢复为二值标记函数，即

$ \phi = \left\{ {\begin{array}{*{20}{c}} {1,}&{\phi \ge \tau ,}\\ {0,}&{{\rm{otherwise}},} \end{array}} \right. $

(26)

其中，τ∈(0, 1)，通常取$ \tau = \frac{1}{2} $。

上述计算过程可综述为

Compute w from (7);u_ij⁰=f⁰；v_ij⁰=▽u_ij; 1⁰=▽ϕ; λ0_ij=0;σ⁰=0;ε=0.001; N=200;E⁰=0;N=200; for(k=0;k＜=N; k++)

1. Compute λ_ij^k+1(i=1, 2;j=1, 2, 3) from (15)

2. Compute u_ij^k+1(i=1, 2;j=1, 2, 3) from (17);

3.Compute v_ij^k+1(i=1, 2;j=1, 2, 3) from (18);

4. Compute ϕ^k+1from (22) and (23);

5. Compute w^k+1 from (21);

6. Compute E^k+1 from (14);

If (25) is fulfilled

  Break

    Else Continue

  Endif

End for

3 实验结果与分析

实验研究分2部分，第一部分通过典型彩色纹理图像的数值实验考察所提出模型对纹理和边缘的保持效果。第二部分的目的在于考察模型及算法与现有模型分割效果的区别。所有实验在(Intel(R) Core(TM).i5-4590 CPU@3.30 GHz，4.00 GB内存)上完成，编程环境为Matlab R2014b。实验数据选自伯克利分割数据集^[26]。

3.1 提出模型对复杂纹理图像处理结果

在图 2中，选取了4幅合成的复杂纹理图像。这4幅图像纹理复杂，并且目标与背景相似，分割难度较大。4幅图像采用像素片大小的轮廓线根据图像的特点在所分割目标附近进行初始化。4幅图像采用的像素片大小均为5×5，前2幅图搜索窗口为9×9，第3、4幅图搜索窗口为7×7。通过实验结果可以看出，模型对于纹理和边缘信息保持较好，可以得到比较理想的分割结果。

在图 3中，轮廓线沿着目标物体演化最终到边缘停止演化，(e)是最终分割结果。4幅图像素片大小均为5，前3幅图像搜索窗口为11，第4幅图搜索窗口为9。第1幅老虎的图像取值轮廓线紧贴老虎边缘，尾巴后面没有精确分割。第2幅飞机图像的分割结果比较精确，保持层与层之间的耦合，并且没有出现过分割。第3幅斑点鱼和第4幅花豹图像目标和背景比较接近，纹理复杂，非局部项▽_NLu_ij可以保持纹理信息，但还有细节分割不精确。总体来说，提出的模型对于不同环境下的自然纹理图像都可以得到比较理想的分割结果。

3.2 与现有方法实验对比

为了进一步证明提出的彩色纹理分割模型的优越性，用客观评价方法F-measure^[27]对模型的分割结果进行评估。通过手工标记的分割标准图像与实验分割结果进行比对计算出准确率(P)和召回率(R)。F-measure计算公式如下

$ F = \frac{{P * R}}{{\xi P + \left( {1 - \xi } \right)R}}, $

(27)

式中：ξ=0.5。F∈[0, 1]，F越接近1代表分割效果越好，如表 1所示。

在图 4中用提出的方法与另外6中彩色纹理图像的分割算法进行对比。伯克利分割数据集提供了分割标准图像。(f)是文献[9]分割结果。从实验结果可以看出，由于纹理信息复杂，JSEG方法^[28]出现了过分割，背景区域也出现过分割，N-cuts^[29]和CTM^[30]也出现同样的过分割现象，并且较严重。MAP-ML方法^[30]目标区域也会出现较少的过分割。文献[9]分割结果出现少量过分割和欠分割。非局部Mumford-Shah-MT模型对于纹理信息的保持很好，没有出现过分割现象，这是因为非局部方法是基于片的相似性而不是点与点之间的，是对局部信息和全局信息的整合，所以纹理信息可以处理好。研究所提出的模型同样使用非局部算子，很好地保持了纹理信息，没有出现过分割现象，由于使用多通道全变差规则项，彩色图像层与层之间的耦合很好，边缘没有出现模糊现象，轮廓线紧贴目标边缘。

表 2是2种非局部方法的消耗时间，根据图 2可知，轮廓线随着迭代次数的增加会逐渐逼近目标边缘，当轮廓线到达目标边缘时，能量函数会收敛，趋于稳态，得到稳态解，分割结束。为了客观评价2种方法的耗时，2个种非局部方法尽可能使用一致的惩罚参数。对公式(13)中的$ \sqrt {\sum\limits_{j = 1}^3 {|\nabla {u_j}{|^2}} } $使用变分时，会产生曲率项，由于曲率项计算复杂，由第3部分可知引入了辅助变量及Lagrange乘子，交替优化迭代，相对于计算$ {\sum\limits_{j = 1}^3 {|\nabla {u_j}{|^2}} } $会消耗较长的时间，但分割精确度提升较大，可以弥补时间上的不足。

4 结语

在综合了经典的Mumford-Shah模型、多通道全变差模型、非局部变分模型基础上提出了彩色纹理图像分割的非局部Mumford-Shah多通道全变差模型，并设计了相应的ADMM算法，将原变分问题分解为一系列可通过Gauss-Seidel迭代和广义软阈值公式快速计算的简单子优化问题。通过大量数值实验比较验证了研究方法较现有方法在纹理图像边缘表达以及图像层与层耦合方面的优越性。未来还可将相关研究推广到曲面纹理图像的分割。