多相图像分割的目的是根据图像特征将其自动划分为互相毗邻且不重叠的多个区域，是对两相图像分割问题研究的扩展，在医学、交通、遥感等领域的图像处理与分析中具有大量实际应用(Li等，2016；Spencer等，2019)。近几年，多相图像分割的变分模型还为高维数据的分析与分类(郑世秀，2017)奠定了重要基础。

变分图像分割方法使用连续或离散水平集函数构造不同区域的特征函数，后者又称为离散标签函数。Potts模型(Potts，1952)是多相图像分割变分模型的基础，Samson等人(2000)提出了使用n个连续水平集函数构造n个区域特征函数的多相图像分割方法，但必须添加单纯形约束以避免间隙和重叠问题。为了自动满足单纯形约束条件，潘振宽等人(2009)提出使用n-1个水平集函数设计n个特征函数的区域划分方法。为减少水平集函数的个数，Vese和Chan(2002)基于四色定理，使用m个水平集函数针对不同区域设计2^m个特征函数。为了进一步降低模型的规模，Chung和Vese(2005)提出用1个水平集函数划分多相区域的分割模型，但特征函数的次数较高，不利于优化求解。Lie等人(2006)将连续水平集函数推广到离散水平集函数，提出用m个二值标签函数定义2^m个特征函数的多相分割模型。Li等人(2010)还提出使用n个离散水平集函数划分n个区域的分割方法。近年来，快速分裂Bregman方法(split Bregman method, SBM)、对偶方法(dual method, DM)、增广拉格朗日方法(augmented Lagrangian method, ALM)应用于全变差(total variation, TV)模型中(Wu和Tai，2010)，避免计算由TV规则项产生的曲率项，并有效提高计算效率。但使用离散水平集表示特征函数，变分模型变为非凸的，通常需要通过标签函数凸松弛(Baxter等，2017)将离散问题转化为连续问题，进一步通过凸优化及阈值化求解能量泛函的最优解。

为了提高分割计算效率，Greig等人(1989)引入一种高效的图割(graph cut，GC)算法，将能量泛函的极值问题转化为最大流/最小割问题。Boykov等人(2001)，Bae和Tai(2009)将最小割方法应用到离散标签函数的图像分割模型。根据最大流和最小割定理，最小割问题可以转换为从源点到汇点的网络最大流问题。离散最大流方法存在度量误差，并且不能并行执行。为了避免该局限性，Yuan等人(2010a)通过引入对偶变量，提出与图割方法等价的连续最大流(continuous max-flow method，CMF)方法求解分割模型(Bae等，2014)。为了提高分割精度，杨晓艺等人(2013)提出空间流约束变量为非负函数的连续最大流模型及算法。为了提高多相图像分割的计算效率，Yuan等人(2010b)将CMF方法拓展到以Potts模型为基础的能量泛函极值问题。

由上述分析可见，多相图像分割Vese-Chan模型所用水平集函数的数量较少，定义的特征函数次数较低，从而模型的计算效率适中、优化计算性态较好，并得到了大量推广。本文拟基于标签函数表达的多相图像分割Vese-Chan模型，设计其快速的连续最大流方法。为此，首先引入多个对偶变量，将能量泛函的极小值问题转化为基于源流的能量泛函极大值问题。然后，引入Lagrange乘子和惩罚参数，为该模型设计交替方向乘子方法(alternating direction method of multipliers，ADMM)(Goldstein等，2014；Xiu等，2019)。最后，本文对所设计的算法进行数值实验，并验证其分割效果与计算效率。

Chan和Vese(2001)提出了基于连续水平集函数的Chan-Vese模型，其定义为

$ \mathop {\min }\limits_{{u_1}, {u_2}, \varphi } \left\{ {\begin{array}{*{20}{l}} {E\left({{u_1}, {u_2}, \varphi } \right) = {\alpha _1}\int_\mathit{\boldsymbol{ \boldsymbol{\varOmega} }} {{Q_1}} H(\varphi){\rm{d}}x + }\\ {{\alpha _2}\int_\mathit{\boldsymbol{ \boldsymbol{\varOmega} }} {{Q_2}} (1 - H(\varphi)){\rm{d}}x + }\\ {\gamma \int_\mathit{\boldsymbol{ \boldsymbol{\varOmega} }} {\left| {\nabla H(\varphi)} \right|} {\rm{d}}x} \end{array}} \right\} $

(1)

式中，$E\left(\cdot \right)$表示模型的能量泛函，等号右边的前两项是数据项，${Q_1} = {\left({f - {u_1}} \right)^2}, {Q_2} = {\left({f - {u_2}} \right)^2}$，$f$为要分割的图像强度，${u_1}, {u_2}$为分割后的分段常值近似值，$\varphi \left(x \right)$是用符号距离函数表示的连续水平集函数，满足Eikonal等式

$ \left| {{\nabla _\varphi }(x)} \right| = 1 $

(2)

$H\left({\varphi \left(x \right)} \right)$为Heaviside函数，式(1)第3项是TV规则项，$\nabla $为梯度算子，${\alpha _1}, {\alpha _2}, \gamma $为惩罚参数，$\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}$表示图像区域。

Chan等人(2006)使用凸方法分割连续图像区域，将离散标签函数$\phi \in \left\{ {0, 1} \right\}$凸松弛为$\phi \in \left[ {0, 1} \right]$，式(1)重写为

$ \mathop {\min }\limits_{{u_1}, {u_2}, \phi \in [0, 1]} \left\{ {\begin{array}{*{20}{l}} {E\left({{u_1}, {u_2}, \phi } \right) = {\alpha _1}\int_\mathit{\boldsymbol{ \boldsymbol{\varOmega} }} {{Q_1}\phi {\rm{d}}x + } }\\ {{\alpha _2}\int_\mathit{\boldsymbol{ \boldsymbol{\varOmega} }} {{Q_2}(1 - \phi){\rm{d}}x + } }\\ {\gamma \int_\mathit{\boldsymbol{ \boldsymbol{\varOmega} }} {|\nabla \phi |{\rm{d}}x} } \end{array}} \right\} $

(3)

对于图像任意一点$x \in \mathit{\boldsymbol{ \boldsymbol{\varOmega} }}$，Yuan等人(2010a)引入3个对偶变量：从源点$s$流向$x$的源流${q_s}$，从$x$流向汇点$t$的汇流${q_t}$，以及空间流$\mathit{\boldsymbol{p}}$，式(3)可表示为

$ \mathop {\min }\limits_{\phi \in [0,1]} {\mkern 1mu} \mathop {\max }\limits_{\begin{array}{*{20}{c}} {{q_t} \le {Q_1},}\\ {\begin{array}{*{20}{c}} {{q_s} \le {Q_2},}\\ {|\mathit{\boldsymbol{p}}| \le \gamma {\rm{ }}} \end{array}} \end{array}} {\mkern 1mu} \left\{ {\begin{array}{*{20}{l}} {E\left( {\phi ,{q_t},{q_s},\mathit{\boldsymbol{p}}} \right) = \int_{\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}{\rm{ }}} {{q_t}} \phi {\rm{d}}x + }\\ {\int_{\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}{\rm{ }}} {{q_s}} (1 - \phi ){\rm{d}}x + \int_{\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}{\rm{ }}} \phi \nabla \cdot \mathit{\boldsymbol{p}}{\rm{d}}x = }\\ {\int_{\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}{\rm{ }}} {{q_s}} {\rm{d}}x + \int_{\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}{\rm{ }}} \phi \left( {{q_t} - {q_s} + \nabla \cdot \mathit{\boldsymbol{p}}} \right){\rm{d}}x} \end{array}} \right\} $

(4)

式中，$\nabla \cdot$为散度算子。由上式的极小值问题得到流守恒条件

$ {q_t} - {q_s} + \nabla \cdot \mathit{\boldsymbol{p}} = 0 $

(5)

由极大值问题可将Chan-Vese模型的能量泛函极小值转化为源流的极大值

$ \mathop {\max }\limits_{{q_t}{Q_1}} \mathop {\max }\limits_{{q_s}{Q_2}} \mathop {\max }\limits_{|\mathit{\boldsymbol{p}}|\gamma } \int_\mathit{\boldsymbol{ \boldsymbol{\varOmega} }} {{q_s}{\rm{d}}x} $

(6)

Potts模型是变分多相图像分割的基础，其图像区域$\mathit{\boldsymbol{ \boldsymbol{\varOmega} = }}\left\{ {{\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}_i}, i = 1, \cdots, n} \right\}$满足

$ \mathit{\boldsymbol{ \boldsymbol{\varOmega} }} = \bigcup\limits_{i = 1}^n {{\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}_i}, {\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}_i}} \bigcap\limits_{i \ne j} {{\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}_j} = } \emptyset $

(7)

其能量泛函的极值问题为

$ \mathop {\min }\limits_{\left\{ {{u_i}} \right\}_{i = 1}^n} \mathop {\min }\limits_{\left\{ {{\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}_i}} \right\}_{i = 1}^n} \left\{ {\sum\limits_{i = 1}^n \alpha \int_{{\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}_i}} {{Q_i}} {\rm{d}}x + \sum\limits_{i = 1}^n \gamma \left| {\partial {\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}_i}} \right|} \right\} $

(8)

式中，${Q_i} = {\left({f - {u_i}} \right)^2}, {u_i}\left({i = 1, \cdots, n} \right)$是各区域${\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}_i}$的分段常值近似值，$\alpha, \gamma $为惩罚参数，$\left| {\partial {\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}_i}} \right|$表示各区域的边界长度或曲面面积。

对于多相分割而言，为每个区域分配一个特征函数${\chi _i}\left(x \right) \in \left\{ {0, 1} \right\}, \left({i = 1, \cdots, n} \right)$，例如，当$x \in {\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}_i}$时，${\chi _i}\left(x \right) = 1$；当$x \notin {\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}_i}$时，${\chi _i}\left(x \right) = 0$, 式(8)重写为

$ \mathop {\min }\limits_{\left\{ {{u_i}} \right\}_{i = 1}^n} \left\{ {\mathop \sum \limits_{i = 1}^n \alpha \int_\mathit{\boldsymbol{ \boldsymbol{\varOmega} }} {{Q_i}{\chi _i}{\rm{d}}x + \mathop \sum \limits_{i = 1}^n \gamma \int_\mathit{\boldsymbol{ \boldsymbol{\varOmega} }} {\left| {{\nabla _{{\chi _i}}}} \right|{\rm{d}}x} } } \right\} $

(9)

为了满足式(7)，需要满足单纯性约束

$ \mathop \sum \limits_{i = 1}^n {\chi _i}(x) = 1 $

(10)

Yuan等人(2010b)提出基于Potts模型的连续最大流方法，给定连续区域$\mathit{\boldsymbol{ \boldsymbol{\varOmega} = }}\left\{ {{\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}_i}, i = 1, \cdots, n} \right\}$，源流${{q_s}}$在每一个${{\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}_i}}$都是相同的，汇流${{q_i}}$为从${{\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}_i}}$的$x$流向汇点$t$的流，空间流${{p_i}}$互不相同，其连续最大流模型为

$ \mathop {\max }\limits_{{q_s}, {q_i}, {\mathit{\boldsymbol{p}}_i}} \int_\mathit{\boldsymbol{ \boldsymbol{\varOmega} }} {{q_s}{\rm{d}}x} $

(11)

且满足以下两个流约束

$ \begin{array}{*{20}{c}} {{q_i}{Q_i}, \;\;\left| {{\mathit{\boldsymbol{p}}_i}} \right|\gamma }\\ {{q_i} - {q_s} + \nabla \cdot {\mathit{\boldsymbol{p}}_i} = 0} \end{array} $

(12)

Vese和Chan提出基于变分水平集框架的多相分割Vese-Chan模型，使用$m$个连续水平集函数${\varphi _j}\left({j = 1, \cdots, m} \right)$定义$n = {2^m}$个特征函数进行区域划分，其能量泛函的极值问题为

$ \mathop {\min }\limits_\varphi \left\{ {E(\varphi) = \mathop \sum \limits_{i = 1}^{{2^m}} \alpha \int_\mathit{\boldsymbol{ \boldsymbol{\varOmega} }} {{Q_i}{\chi _i}{\rm{d}}x + \mathop \sum \limits_{j = 1}^m \gamma } \int_\mathit{\boldsymbol{ \boldsymbol{\varOmega} }} {\left| {\nabla H\left({{\varphi _j}} \right)} \right|{\rm{d}}x} } \right\} $

(13)

式中，特征函数${{\chi _i}}$满足

$ \begin{array}{*{20}{c}} {{\chi _i} = \prod\limits_{j = 1}^m \mathit{\Psi } \left({{\varphi _j}} \right)}\\ {\mathit{\Psi }\left({{\varphi _j}} \right) = \left\{ {\begin{array}{*{20}{l}} {H\left({{\varphi _j}} \right)}&{x \in {\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}_j}}\\ {1 - H\left({{\varphi _j}} \right)}&{x \notin {\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}_j}} \end{array}} \right.} \end{array} $

(14)

连续水平集函数${{\varphi _j}}$满足约束条件式(2)。

Lie等人(2006)使用$m$个二值标签函数${\phi _j} \in \left\{ {0, 1} \right\}\left({j = 1, \cdots, m} \right)$定义$n = {2^m}$个特征函数，将${\phi _j}$凸松弛为${\phi _j} \in \left[ {0, 1} \right]$，基于标签函数的Vese-Chan模型为

$ \mathop {\min }\limits_{\phi \in [0, 1]} \left\{ {E(\phi) = \mathop \sum \limits_{i = 1}^{{2^m}} \alpha \int_\mathit{\boldsymbol{ \boldsymbol{\varOmega} }} {{Q_i}{\chi _i}{\rm{d}}x + \mathop \sum \limits_{j = 1}^m \gamma } \int_\mathit{\boldsymbol{ \boldsymbol{\varOmega} }} {\left| {\nabla {\phi _j}} \right|{\rm{d}}x} } \right\} $

(15)

式中，特征函数${\chi _i}$为

$ \begin{array}{l} {\chi _i} = \prod\limits_{j = 1}^m \mathit{\Psi } \left({{\phi _j}} \right)\\ \mathit{\Psi }\left({{\phi _j}} \right) = \left\{ {\begin{array}{*{20}{l}} {{\phi _j}}&{x \in {\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}_j}}\\ {1 - {\phi _j}}&{x \notin {\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}_j}} \end{array}} \right. \end{array} $

(16)

王琦等人(2010)根据要划分的区域编号与其二进制表达的关系建立特征函数${\chi _i}$的通用表达形式

$ {\chi _i} = \mathop \prod \limits_{j = 1}^m \left({b_j^{i - 1} + {{(- 1)}^{b_j^{i - 1}}}{\phi _j}} \right) $

(17)

式中，区域的编号为$i = 1, \cdots, {2^m}$，第$i$个区域的二进制表达式为$b_1^{i - 1}, b_2^{i - 1}, \cdots, b_m^{i - 1}$，其中$b_j^{i - 1} = \left\{ {0, 1} \right\}$。以$m = 2$为例，特征函数可表示为

$ \left\{ {\begin{array}{*{20}{l}} {{\chi _1} = {\phi _1}{\phi _2}}\\ {{\chi _2} = {\phi _1}\left({1 - {\phi _2}} \right)}\\ {{\chi _3} = \left({1 - {\phi _1}} \right){\phi _2}}\\ {{\chi _4} = \left({1 - {\phi _1}} \right)\left({1 - {\phi _2}} \right)} \end{array}} \right. $

(18)

为了提高计算效率，Tan等人(2015)为Vese-Chan模型设计ADMM方法，引入辅助变量${\mathit{\boldsymbol{w}}_j} = {\nabla _{{\phi _j}}}$、拉格朗日乘子${\lambda _j}$和惩罚参数$\theta $，式(15)可表示为

$ \mathop {\min }\limits_{\phi \in [0, 1], \mathit{\boldsymbol{w}}} \left\{ {\begin{array}{*{20}{l}} {E(\phi, \mathit{\boldsymbol{w}}) = }\\ {\sum\limits_{i = 1}^{{2^m}} \alpha \int_\mathit{\boldsymbol{ \boldsymbol{\varOmega} }} {{Q_i}} {\chi _i}{\rm{d}}x + \sum\limits_{j = 1}^m \gamma \int_\mathit{\boldsymbol{ \boldsymbol{\varOmega} }} {\left| {{\mathit{\boldsymbol{w}}_j}} \right|} {\rm{d}}x + }\\ {\sum\limits_{j = 1}^m {\int_\mathit{\boldsymbol{ \boldsymbol{\varOmega} }} {{\lambda _j}} } \cdot \left({{\mathit{\boldsymbol{w}}_j} - \nabla {\phi _j}} \right){\rm{d}}x + }\\ {\frac{\theta }{2}\sum\limits_{j = 1}^m {\int_\mathit{\boldsymbol{ \boldsymbol{\varOmega} }} {{{\left| {{\mathit{\boldsymbol{w}}_j} - \nabla {\phi _j}} \right|}^2}} } {\rm{d}}x} \end{array}} \right\} $

(19)

变量${\phi _j}$和${\mathit{\boldsymbol{w}}_j}$的优化计算分别为

$ \left\{ {\begin{array}{*{20}{c}} {\sum\limits_{i = 1}^{{2^m}} \alpha {Q_i}\frac{{\partial {\chi _i}}}{{\partial {\phi _j}}} + \nabla \cdot \lambda _j^k + }&{x \in \mathit{\boldsymbol{ \boldsymbol{\varOmega} }}}\\ {\theta \;\nabla \cdot \left({\mathit{\boldsymbol{w}}_j^k - \nabla {\phi _j}} \right) = 0}&{}\\ { - \left({\lambda _j^k + \mathit{\boldsymbol{w}}_j^k - \nabla {\phi _j}} \right) \cdot \mathit{\boldsymbol{n}} = 0}&{x \in \partial \mathit{\boldsymbol{ \boldsymbol{\varOmega} }}} \end{array}} \right. $

(20)

$ \mathit{\boldsymbol{w}}_j^{k + 1} = \max \left({\left| {\nabla \phi _j^{k + 1} - \frac{{\lambda _j^k}}{\theta }} \right| - \frac{\gamma }{\theta }, 0} \right)\frac{{\nabla \phi _j^{k + 1} - \frac{{\lambda _j^k}}{\theta }}}{{{{\left| {\nabla \phi _j^{k + 1} - \frac{{\lambda _j^k}}{\theta }} \right|}_\zeta }}} $

(21)

式中，$k$表示迭代次数，$\zeta $为极小正数值以防止分母为零。

连续最大流方法应用到两相图像分割模型中，当扩展到多相分割时，需要引入大量的对偶变量。从1.1节基于Chan-Vese模型的连续最大流方法中提取关键思想，根据式(17)重新设计特征函数

$ {\chi _i} = \left\{ {\begin{array}{*{20}{l}} {{\eta _{ij}}{\phi _j}}&{b_j^{i - 1} = 0}\\ {{\eta _{ij}}\left({1 - {\phi _j}} \right)}&{b_j^{i - 1} = 1} \end{array}} \right. $

(22)

式中，${\eta _{ij}}\left({i = 1, \cdots, {2^m}, j = 1, \cdots, m} \right)$表示$m + 1$个二值标签函数的乘积(不包含${\phi _j}$)，定义如下

$ {\eta _{ij}} = \mathop \prod \limits_{\begin{array}{*{20}{c}} {l = 1}\\ {l \ne j} \end{array}}^m \left({b_l^{i - 1} + {{(- 1)}^{b_l^{i - 1}}}{\phi _l}} \right) $

(23)

以2个二值标签函数为例，二进制区域编号与${\eta _{ij}}$的关系如表 1所示。

当$j$取$1, \cdots, m$中某一值时，根据第$i$个区域二进制编号$b_j^{i - 1}$的取值为0或1，将函数${\eta _{ij}}$划分为$\eta _{ij}^0$和$\eta _{ij}^1$两类。当$b_j^{i - 1} = 0$时，将函数${\eta _{ij}}$表示为$\eta _{ij}^0$；当$b_j^{i - 1} = 1$时，将函数${\eta _{ij}}$表示为$\eta _{ij}^1$，因此多相分割模型中的数据项被划分为两类。将Vese-Chan模型转化为交替优化的两相分割模型，其能量泛函的极值问题如下所示

$ \mathop {\min }\limits_{\phi \in [0, 1]} \left\{ {\begin{array}{*{20}{l}} {E(\phi) = \mathop \sum \limits_{i = 1}^{{2^m}} \alpha \int_\mathit{\boldsymbol{ \boldsymbol{\varOmega} }} {{Q_i}\eta _{ij}^0{\phi _j}{\rm{d}}x + } }\\ {\mathop \sum \limits_{i = 1}^{{2^m}} \alpha \int_\mathit{\boldsymbol{ \boldsymbol{\varOmega} }} {} {Q_i}\eta _{ij}^1\left({1 - {\phi _j}} \right){\rm{d}}x + \gamma \int_\mathit{\boldsymbol{ \boldsymbol{\varOmega} }} {|\overline \nabla \phi |{\rm{d}}x} } \end{array}} \right\} $

(24)

式中，${\bar \nabla }$表示对标量参数的矢量化算子，${{\bar \nabla }_\phi }$的表达式为

$ {{\bar \nabla }_\phi } = \left[ {\begin{array}{*{20}{c}} {\nabla {\phi _1}}\\ {\nabla {\phi _2}}\\ \vdots \\ {\nabla {\phi _m}} \end{array}} \right] $

(25)

为了形象地表示两相分割模型，令${Q^0} = \mathop \sum \limits_{i = 1}^{{2^m}} \alpha {Q_i}\eta _{ij}^0$，${Q^1} = \mathop \sum \limits_{i = 1}^{2m} \alpha {Q_i}\eta _{ij}^1$，整理式(24)得到

$ \mathop {\min }\limits_{\phi \in \left[ {0, 1} \right]} \left\{ {E\left(\phi \right) = \int_\mathit{\boldsymbol{ \boldsymbol{\varOmega} }} {{Q^0}} {\phi _j}{\rm{d}}x + \int_\mathit{\boldsymbol{ \boldsymbol{\varOmega} }} {{Q^1}\left({1 - \phi } \right){\rm{d}}x + \gamma \int_\mathit{\boldsymbol{ \boldsymbol{\varOmega} }} {\left| {{{\overline \nabla }_\phi }} \right|{\rm{d}}x} } } \right\} $

(26)

以$m = 2$为例，当$j = 1$时

$ \begin{array}{*{20}{c}} {E(\phi) = \alpha {Q_1}\eta _{11}^0{\phi _1} + \alpha {Q_2}\eta _{21}^0{\phi _1} + }\\ {\alpha {Q_3}\eta _{31}^1\left({1 - {\phi _1}} \right) + \alpha {Q_4}\eta _{41}^1\left({1 - {\phi _1}} \right) = }\\ {\left({\alpha {Q_1}{\phi _2} + \alpha {Q_2}\left({1 - {\phi _2}} \right)} \right){\phi _1} + }\\ {\left({\alpha {Q_3}{\phi _2} + \alpha {Q_4}\left({1 - {\phi _2}} \right)} \right)\left({1 - {\phi _1}} \right)} \end{array} $

(27)

当$j = 2$时

$ \begin{array}{*{20}{c}} {E(\phi) = \alpha {Q_1}\eta _{12}^0{\phi _2} + \alpha {Q_2}\eta _{22}^1\left({1 - {\phi _2}} \right) + }\\ {\alpha {Q_3}\eta _{32}^0{\phi _2} + \alpha {Q_4}\eta _{42}^1\left({1 - {\phi _2}} \right) = }\\ {\left({\alpha {Q_1}{\phi _1} + \alpha {Q_3}\left({1 - {\phi _1}} \right)} \right){\phi _2} + }\\ {\left({\alpha {Q_2}{\phi _1} + \alpha {Q_4}\left({1 - {\phi _1}} \right)} \right)\left({1 - {\phi _2}} \right)} \end{array} $

(28)

由此可见，$j$取不同的值，$E\left(\phi \right)$的表达式是一致的。

式(26)已将多相分割模型转化为交替优化的两相分割模型，因此，引入矢量化对偶变量${\mathit{\boldsymbol{q}}_t}, {\mathit{\boldsymbol{q}}_s}, \mathit{\boldsymbol{p}}$，采用连续最大流方法求解Vese-Chan模型能量泛函的极值问题，转化为关于源流的极大值问题

$ \mathop {\max }\limits_{\left| {{\mathit{\boldsymbol{q}}_t}} \right|{Q^0}|} \mathop {\max }\limits_{\left| {{\mathit{\boldsymbol{q}}_s}} \right|{Q^1}} \mathop {\max }\limits_{|\mathit{\boldsymbol{p}}|\gamma } \left\{ {E\left({{\mathit{\boldsymbol{q}}_t}, {\mathit{\boldsymbol{q}}_s}, \mathit{\boldsymbol{p}}} \right) = \int_\mathit{\boldsymbol{ \boldsymbol{\varOmega} }} {{\mathit{\boldsymbol{q}}_s}{\rm{d}}x} } \right\} $

(29)

并且满足流约束条件

$ {\mathit{\boldsymbol{q}}_t} - {\mathit{\boldsymbol{q}}_s} + \overline \nabla \cdot \mathit{\boldsymbol{p}} = 0 $

(30)

式中，$\overline \nabla \cdot $表示对矢量参数的矢量化算子，$\overline \nabla \cdot \mathit{\boldsymbol{p}}$的表达为

$ \overline \nabla \cdot \mathit{\boldsymbol{p}} = \left[ {\begin{array}{*{20}{c}} {\nabla \cdot {\mathit{\boldsymbol{p}}_1}}\\ {\nabla \cdot {\mathit{\boldsymbol{p}}_2}}\\ \vdots \\ {\nabla \cdot {\mathit{\boldsymbol{p}}_m}} \end{array}} \right] $

(31)

为了提高计算效率，为本文模型设计ADMM方法，引入拉格朗日乘子$\lambda $和惩罚参数$\theta $，得到能量泛函的交替优化形式

$ \begin{array}{l} \left({\mathit{\boldsymbol{q}}_t^{k + 1}, \mathit{\boldsymbol{q}}_s^{k + 1}, {\mathit{\boldsymbol{p}}^{k + 1}}} \right) = \\ \arg \mathop {\max }\limits_{\left| {{\mathit{\boldsymbol{q}}_t}} \right| \le {Q^0}|} \mathop {\max }\limits_{\left| {{\mathit{\boldsymbol{q}}_s}} \right| \le {Q^1}} \mathop {\max }\limits_{|\mathit{\boldsymbol{p}}| \le \gamma } \left\{ {\begin{array}{*{20}{l}} {E\left({{\mathit{\boldsymbol{q}}_t}, {\mathit{\boldsymbol{q}}_s}, \mathit{\boldsymbol{p}}} \right) = \int_\mathit{\boldsymbol{ \boldsymbol{\varOmega} }} {{\mathit{\boldsymbol{q}}_s}} {\rm{d}}x + }\\ {\left. {\int_\mathit{\boldsymbol{ \boldsymbol{\varOmega} }} {{\lambda ^k}} \cdot \left({{\mathit{\boldsymbol{q}}_t} - {\mathit{\boldsymbol{q}}_s} + \bar \nabla \cdot \mathit{\boldsymbol{p}}} \right){\rm{d}}x - } \right\}}\\ {\frac{\theta }{2}\int_\mathit{\boldsymbol{ \boldsymbol{\varOmega} }} {{{\left| {{\mathit{\boldsymbol{q}}_t} - {\mathit{\boldsymbol{q}}_s} + \bar \nabla \cdot \mathit{\boldsymbol{p}}} \right|}^2}} {\rm{d}}x} \end{array}} \right\} \end{array} $

(32)

交替优化的主要思想是固定其他变量求解一个变量的极值问题，因此，式(32)可转化为以下子优化问题

$ \left\{ {\begin{array}{*{20}{l}} {\mathit{\boldsymbol{q}}_t^{k + 1} = \arg \mathop {\max }\limits_{\left| {{\mathit{\boldsymbol{q}}_t}} \right|{Q^0}} E\left({{\mathit{\boldsymbol{q}}_t}, \mathit{\boldsymbol{q}}_s^k, {\mathit{\boldsymbol{p}}^k}} \right)}\\ {\mathit{\boldsymbol{q}}_s^{k + 1} = \arg \mathop {\max }\limits_{\left| {{\mathit{\boldsymbol{q}}_s}} \right|{Q^1}} E\left({\mathit{\boldsymbol{q}}_t^{k + 1}, {\mathit{\boldsymbol{q}}_s}, {\mathit{\boldsymbol{p}}^k}} \right)}\\ {{\mathit{\boldsymbol{p}}^{k + 1}} = \arg \mathop {\max }\limits_{|\mathit{\boldsymbol{p}}|\gamma } E\left({\mathit{\boldsymbol{q}}_t^{k + 1}, \mathit{\boldsymbol{q}}_s^{k + 1}, \mathit{\boldsymbol{p}}} \right)} \end{array}} \right. $

(33)

对于式(33)，固定${\mathit{\boldsymbol{q}}_s^k}$和${{\mathit{\boldsymbol{p}}^k}}$，得到关于${\mathit{\boldsymbol{q}}_t^{k + 1}}$的Euler-Lagrange方程

$ {\mathit{\boldsymbol{\lambda }}^k} - \theta \left({{\mathit{\boldsymbol{q}}_t} - \mathit{\boldsymbol{q}}_s^k + \overline \nabla \cdot {\mathit{\boldsymbol{p}}^k}} \right) = 0 $

(34)

${\mathit{\boldsymbol{q}}_t^{k + 1}}$的解析解为

$ \left\{ {\begin{array}{*{20}{l}} {\mathit{\boldsymbol{q}}_t^{k + 1} = \frac{{{\mathit{\boldsymbol{\lambda }}^k}}}{\theta } + \mathit{\boldsymbol{q}}_s^k - \bar \nabla \cdot {\mathit{\boldsymbol{p}}^k}}\\ {\mathit{\boldsymbol{q}}_t^{k + 1} = \max \left({\min \left({\left| {\mathit{\boldsymbol{q}}_t^{k + 1}} \right|, {Q^0}} \right), 0} \right)} \end{array}} \right. $

(35)

固定$\mathit{\boldsymbol{q}}_t^{k + 1}$和${\mathit{\boldsymbol{p}}^k}$，得到关于$\mathit{\boldsymbol{q}}_s^{k + 1}$的Euler-Lagrange方程

$ \mathit{\boldsymbol{1}} - {\mathit{\boldsymbol{\lambda }}^k} + \theta \left({\mathit{\boldsymbol{q}}_t^{k + 1} - {\mathit{\boldsymbol{q}}_s} + \overline \nabla \cdot {\mathit{\boldsymbol{p}}^k}} \right) = 0 $

(36)

$\mathit{\boldsymbol{q}}_s^{k + 1}$的解析解为

$ \left\{ {\begin{array}{*{20}{l}} {\mathit{\boldsymbol{q}}_s^{k + 1} = \frac{{1 - {\mathit{\boldsymbol{\lambda }}^k}}}{\theta } + \mathit{\boldsymbol{q}}_t^{k + 1} + \bar \nabla \cdot {\mathit{\boldsymbol{p}}^k}}\\ {\mathit{\boldsymbol{q}}_s^{k + 1} = \max \left({\min \left({\left| {\mathit{\boldsymbol{q}}_s^{k + 1}} \right|, {Q^1}} \right), 0} \right)} \end{array}} \right. $

(37)

固定$\mathit{\boldsymbol{q}}_t^{k + 1}$和$\mathit{\boldsymbol{q}}_s^{k + 1}$得到${\mathit{\boldsymbol{p}}^{k + 1}}$的Euler-Lagrange方程

$ \left\{ {\begin{array}{*{20}{c}} { - \overline \nabla {\mathit{\boldsymbol{\lambda }}^k} + \theta \overline \nabla \left({\mathit{\boldsymbol{q}}_t^{k + 1} - \mathit{\boldsymbol{q}}_s^{k + 1} + \overline \nabla \cdot \mathit{\boldsymbol{p}}} \right) = 0\;\;\;x \in \mathit{\boldsymbol{ \boldsymbol{\varOmega} }}}\\ {{\mathit{\boldsymbol{\lambda }}^k} - \theta \left({\mathit{\boldsymbol{q}}_t^{k + 1} - \mathit{\boldsymbol{q}}_s^{k + 1} + \overline \nabla \cdot \mathit{\boldsymbol{p}}} \right) = 0\;\;\;\;\;\;\;\;x \in \partial \mathit{\boldsymbol{ \boldsymbol{\varOmega} }}} \end{array}} \right. $

(38)

并通过投影公式计算${\mathit{\boldsymbol{p}}^{k + 1}}$的精确值

$ {\mathit{\boldsymbol{p}}^{k + 1}} = \frac{{{\mathit{\boldsymbol{p}}^{k + 1}}}}{{\max \left({1, \frac{{\left| {{\mathit{\boldsymbol{p}}^{k + 1}}} \right|}}{\gamma }} \right)}} $

(39)

拉格朗日乘子${{\mathit{\boldsymbol{\lambda }}^{k + 1}}}$被视为二值标签函数$\phi $的近似值，其更新计算为

$ \left\{ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{\lambda }}^{k + 1}} = {\mathit{\boldsymbol{\lambda }}^k} - \theta \left({\mathit{\boldsymbol{q}}_t^{k + 1} - \mathit{\boldsymbol{q}}_s^{k + 1} + \overline \nabla \cdot {\mathit{\boldsymbol{p}}^{k + 1}}} \right)}\\ {{\mathit{\boldsymbol{\lambda }}^{k + 1}} = \max \left({\min \left({\left| {{\mathit{\boldsymbol{\lambda }}^{k + 1}}} \right|, 1} \right), 0} \right)} \end{array}} \right. $

(40)

为了表达分割图像的边界，需要将凸松弛后的二值标签函数$\phi \in \left[ {0, 1} \right]$转化为$\phi \in \left\{ {0, 1} \right\}$，对$\phi $的阈值化如下所示

$ \phi = \left\{ {\begin{array}{*{20}{c}} {1\;\;\phi \ge \tau }\\ {0\;\;\phi < \tau } \end{array}} \right. $

(41)

式中，$\tau \in \left({0, 1} \right)$，参数$\tau $通常取$\frac{1}{2}$。当能量误差公式满足

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

(42)

时, 交替优化过程结束，其中${E^k}$为第$k$步能量泛函的值, $\varepsilon $为极小正数值。

上述多相图像分割Vese-Chan模型的连续最大流方法描述如下:

定义:${u^0} = {f^0}, \left({\mathit{\boldsymbol{q}}_t^0, \mathit{\boldsymbol{q}}_t^0, {\mathit{\boldsymbol{p}}^0}, {\mathit{\boldsymbol{\lambda }}^0}} \right) = 0$, step=500;for$\left({k = 0;k < = step;k + + } \right)$

根据式(35)计算$\mathit{\boldsymbol{q}}_t^{k + 1}$;

根据式(37)计算$\mathit{\boldsymbol{q}}_s^{k + 1}$;

根据式(38)(39)计算${\mathit{\boldsymbol{p}}^{k + 1}}$;

根据式(40)计算${\mathit{\boldsymbol{\lambda }}^{k + 1}}$;

若满足式(42), 优化结束

end for

Duan等人(2014)将SBM方法、DM方法和ALM方法应用到Vese-Chan分割模型中，并验证以上方法具有近似而且快速的计算效率，而ADMM方法等价于ALM方法，因此，对灰度图像和彩色图像分别采用ADMM方法(Tan等，2015)和本文提出的连续最大流方法进行分割实验，对比两者的分割效果及计算效率。为了客观比较实验结果，对每个实验中的两种方法使用相同的初始化位置，相同的参数$\alpha $和$\lambda $，以保证模型的一致性，并固定参数$\varepsilon = {10^{ - 4}}$以获得相同精度的分割结果。所有实验在PC (Intel(R) Core(TM) i5-4590 CPU @ 3.30 GHz，4.00 GB内存)上完成，编程环境为MATLAB R2016a。

对图 1所示的4幅灰度图像进行实验，式(24)中数据项的表达式为

$ {Q_i} = {\left({f - {u_i}} \right)^2}, i = 1, \cdots, {2^m} $

(43)

图 1(a)-(c))包含3相区域, 图 1 (d)包含4相区域，因此，对以上图像均使用2个二值标签函数进行演化，初始化二值标签函数为2个圆，在圆内$\phi = 1$，否则$\phi = 0$。

图 2是对图 1中4幅图像的分割结果对比, 参数$\alpha = 1$。图 2第1行参数$\gamma = 1$，从分割结果对比看出，本文方法分割更细致，可以获得更精确的大脑CT图像分割结果。图 2第2行参数$\gamma = 1$，比较图中黑色杯子的手柄处，本文方法实现杯子与背景更好的分离。图 2第3行参数$\gamma = 50$，图 2第4行参数$\gamma = 5$，从视觉效果看，两种方法得到一致的分割效果。可以得出，本文方法和ADMM方法对组合图像和几何图像的分割具有通用性，而在医学图像及自然图像的处理上，前者具有更精确的分割结果。

表 2为对图 1实验的迭代次数及运算时间。当模型能量泛函满足收敛条件时，迭代停止，两种方法的收敛条件是一致的。相比较ADMM方法，本文方法的迭代次数均有减少，表明该方法提高收敛速度。从运算时间看，本文方法所需运算时间少于ADMM方法，提高计算效率。图 1(c)的图像大小在前3幅3相图像中最大，运算时间最长，运算时间加速比最大；图 1(d)是简单的几何图像，迭代次数和运算时间较少。可以看出，随着图像尺寸及图像复杂度的增加，本文方法的计算效率提高越明显。

对图 3所示的3幅灰度图像进行实验。这3幅图像的区域相数均大于4，因此，使用3个二值标签函数进行演化。

图 4是对图 3中3幅图像的分割结果对比，参数$\alpha = 10$。图 4第1行参数$\gamma = 1$，图 4第2行参数$\gamma = 10$，从两个简单几何图像的分割结果看，ADMM方法和本文方法具有相同的效果。图 4第3行是从分割结果得到的灰度图，参数为$\gamma = 5$，比较灰度图的右上角区域，可以看出，本文方法的分割结果更加精确，与原图像更加相近。表 3是对图 3进行分割实验的迭代次数及运算时间。图 3(a)是5相图像，图 3(b)是6相图像，且图 3(b)的图像尺寸大于其余两幅图像，所以运算时间是最长的，运算时间加速比最大，因此，本文方法计算效率的提高也与图像所包含的区域相数有关。

对图 5所示的2幅彩色图像进行实验，由于彩色图像包括3层，式(24)中数据项的表达式为

$ {Q_i} = \sum\limits_{c = 1}^3 {{{\left({f - {u_i}} \right)}^2}, i = 1, \cdots, {2^m}} $

(44)

式中，$c$表示彩色图像的分层数，$u_i^c$表示第$c$层第$i$个区域的分段常值近似值。图 5每种颜色代表不同的区域相数，这2幅图像的区域相数均大于8，因此，使用4个二值标签函数进行演化。

图 6是对图 5中2幅图像的分割结果对比，参数$\alpha = 5$。图 6第1行是根据两种方法得到的近似分割结果，参数为$\gamma = 10$，从图中左下角黑色影子部分可以看出，本文方法更加精确。图 6第2行参数为$\gamma = 50$，ADMM方法和本文方法具有相同的分割效果。表 4为对图 5进行分割实验的迭代次数及运算时间对比。

ADMM方法通过引入辅助变量和拉格朗日乘子，避免复杂曲率项的计算，并通过对变量的交替优化提高计算效率。本文方法引入对偶变量，转化为连续最大流模型，保持高效性及高精度性。从表 2—表 4的数据可以得出，分割的迭代次数与运算时间随着标签函数的增加而增加，并且相比较ADMM方法，本文方法的运算时间加速比有明显提高。除此之外，计算效率的提高程度也与图像的复杂程度及尺寸大小有关。

为了进一步提高多相图像分割模型的计算效率，本文提出基于凸松弛的多相图像分割Vese-Chan模型的连续最大流方法。首先，为Vese-Chan模型设计新的特征函数，根据特征函数将模型的数据项整理为两项，从而将多相分割模型转化为交替优化的两相分割模型，便于理解与实现。其次，引入对偶变量，将Vese-Chan模型能量泛函的极小值问题转化为源流的极大值问题，为了提高计算效率，并为该模型设计ADMM方法，将能量泛函的极值问题转化为一系列简单的子优化问题。最后，对不同的多相灰度图像和彩色图像进行数值实验，对比ADMM方法和本文方法的分割效果及计算效率。实验结果表明，相比较ADMM方法，本文方法迭代次数少，具有快速的收敛速度和计算效率，并且对医学图像、遥感图像等复杂图像具有更精确的分割结果。

但是本文方法对于复杂图像的计算效率及分割精度需要进一步提高，今后可以将本文方法应用到医学图像的3维重建等方面的研究中。