1 Introduction

Multiphase image segmentation under variational framework has found a lot of applications including multi-target detection and recognition, 3D segmentation and reconstruction in medical images, remote sensing images, etc. [1, 2], due to its property of multiple cue integration. The aim of multiphase image segmentation is to partition images into different regions without any overlaps and without any unlabeled region (called in the sequel vacuum) automatically. It is a natural extension of the two-phase image segmentation based on the variational image analysis paradigm.

The Mumford-Shah model [3] is fundamental to variational image segmentation; it is a region-based model which approximates an image to a piecewise smooth one and edges. To circumvent the difficulty of its implementation, Chan and Vese [4] proposed the classical Vese-Chan model under variational level set framework [5] based on reduced Mumford-Shah model with piecewise constant image assumption. The Chan-Vese model introduced one-level set function to construct two characteristic functions for two regions. Using the same concept, [6] introduced n-level set functions to partition n regions for Potts model [7], but a simplex constraint must be added to avoid vacuum and overlapping problem. In order to improve computational efficiency from the viewpoint of model, Pan et al. [8] proposed to use (n−1)-level set functions to design n characteristic functions satisfying the simplex constraint naturally. Vese and Chan [9] proposed a strategy using m-level set functions to design n=2m characteristic functions for different regions with less evolution equations to be solved. Furthermore, Chung and Vese [10] proposed a more efficient scheme using only one smooth function to partition different regions.

Motivated by the relationship between the Heaviside function of a level set function and a binary label function, the piecewise constant level set function method was adopted to two-phase image segmentation [11] combined with convex relaxation and thresholding techniques with high efficiency. The models for multiphase image segmentation using variational level set method as mentioned in the previous paragraph have been extended to the counterparts using piecewise constant level set function method successively, such as, [12] used n binary label functions to partition n regions, [13] used m binary label functions for n=2m regions, and [14] used one piecewise constant level set function for all regions. Comparatively, Vese-Chan type scheme uses less label functions than [12], and its characteristic functions have lower degree than [14]. Obviously, lower degree characteristic functions are convenient to achieve minimum of the energy functional. To achieve higher computation efficiency for the optimization problems, Goldstein and Osher [15] proposed the split Bregman method, and Duan et al. [16] proposed some fast projection methods without re-initialization.

Our research in this paper starts from Vese-Chan model via binary label functions [11], and the goal is to improve the computational efficiency through modifying its original model with fewer parameters to be estimated. Using Vese-Chan model, the number n of regions is determined by binary label function number m (i.e., n=2m). For instance, if we partition 9 regions, we need 4 binary functions to construct 16 characteristic functions, and thus 7 regions are empty; however, the parameters in these empty regions must be estimated when using the original model. To avoid the parameter calculations in these redundant regions, [17] proposed a modified scheme for characteristic functions after 2m−1 regions that include two parts to discard the last empty regions, which can overcome the problem of redundant parameter estimations.

In this paper, we use a unified characteristic function expression for all regions including redundant regions, but we add some simple area constraints of redundant regions to avoid estimations of redundant parameters, thus reducing costs. Although we transform the original model into a constrained optimization problem, we use ADMM [18–20] to solve it easily and systematically without additional constructions of characteristic functions for empty regions.

The outline of the paper is as follows. In Section 2, we present some classical models of variational image segmentation, covering the Potts model and Vese-Chan model for multiphase image segmentation using binary label functions. In Section 3, we propose the modified Vese-Chan model and design its ADMM. In Section 4, experimental results of gray images and color images are given to illustrate the efficiency of the proposed model in this paper. Concluding remarks are drawn finally.

2 Some previous works of variational image segmentation

2.1 Image segmentation under variational framework

The Mumford-Shah model and the Potts model are fundamental to variational image segmentation. The former one is to partition an image f(x):Ω→R into piecewise smooth parts u(x):Ωs→R and edges Γ such that Ω=Ωs∪Γ and Ωs∩Γ=Ø, its energy functional minimization problem is stated as

$$ \begin{aligned} {\underset{u,\Gamma}{\text{Min}}} &\left\{ E\left({u,\Gamma} \right) = \alpha \int_{{\Omega^{s}}} {{{\left({f - u} \right)}^{2}}dx}\right.\\ &\left. + \beta \int_{{\Omega^{s}}} {{{\left| {\nabla u} \right|}^{2}}dx} + \gamma \oint_{\Gamma} {ds} \right\}, \end{aligned} $$
(1)

where α,β,γ are penalty parameters. The second one is to partition the image domain Ω into

$$ \Omega = \overset{n}{\underset{i = 1}{\cup}} {\Omega_{i}} \quad \text{and} \quad {\Omega_{i}} \cap {\Omega_{j}} = \text{\O}. $$
(2)

Its energy functional minimization problem is

$$ \begin{aligned} \underset{\left\{ {u_{i}} \right\}_{i = 1}^{n}}{\text{Min}} \underset{\left\{ {\Omega_{i}} \right\}_{i = 1}^{n}}{\text{Min}} &\left\{ {\alpha} \sum\limits_{i = 1}^{n} { \int_{{\Omega_{i}}} {{{\left({f - {u_{i}}} \right)}^{2}}dx}}\right.\\ &\left. + \gamma \sum\limits_{i = 1}^{n} {\left| {\partial {\Omega_{i}}} \right|} \right\}, \end{aligned} $$
(3)

where |∂Ωi| stands for the boundary length of each sub-region.

In order to solve these problems in image domain, one approach is assigning a characteristic function for each region and making use of its total variation to replace length terms. If zero level set {x:ϕ(x)=0} of a level set function ϕ(x) is used to describe the curve for region partition implicitly, its Heaviside function H(ϕ(x)) can be used to define characteristic functions of regions. ϕ(x) and H(ϕ(x)) are defined respectively as

$$ \phi \left(x \right) = \left\{ {\begin{array}{cc} { > 0} & {x \ \text{inside} \ \Gamma} \\ 0 & {x \ \text{on} \ \Gamma} \\ { < 0} & {x \ \text{outside} \ \Gamma} \\ \end{array}} \right., $$
(4)
$$ H\left({\phi \left(x \right)} \right) = \left\{ {\begin{array}{ll} 1 & {\phi \left(x \right) \ge 0} \\ 0 & {\text{otherwise}} \\ \end{array}} \right.. $$
(5)

Then, the length of a closed curve Γ is

$$ \oint_{\Gamma} {ds} = \int_{\Omega} {\left| {\nabla H\left({\phi \left(x \right)} \right)} \right|} dx. $$
(6)

Using level set method, Chan and Vese [4] proposed the Vese-Chan model for two-phase segmentation under variational level set framework based on reduced Mumford-Shah model with piecewise constant image assumption, and it is stated as

$$ \begin{aligned} \underset{{u_{1}},{u_{2}},\phi}{\text{Min}} &\left\{ E\left({{u_{1}},{u_{2}},\phi} \right) = {\alpha_{1}}\int_{\Omega} {{{\left({f - {u_{1}}} \right)}^{2}}H\left(\phi \right)} dx \right.\\ & + {\alpha_{2}}\int_{\Omega} {{{\left({f - {u_{2}}} \right)}^{2}}\left({1 - H\left(\phi \right)} \right)dx} \\ & + \left. {\gamma \int_{\Omega} {\left| {\nabla H\left(\phi \right)} \right|} dx} \right\}. \end{aligned} $$
(7)

If ϕ(x) is defined as a signed distance function, it must fulfill the following Eikonal equation

$$ \left| {\nabla \phi} \right| = 1. $$
(8)

For the case of multiphase segmentation, if each region is assigned a characteristic function χi(x)∈{0,1},(i=1,...,n), (i.e., χi(x)=1 for x∈Ωi,χi(x)=0 for x∉Ωi), the Potts model (3) can be rewritten as

$$ \begin{aligned} \underset{\left\{ {{u_{i}}} \right\}_{i = 1}^{n}}{\text{Min}} \underset{\left\{ {{\Omega_{i}}} \right\}_{i = 1}^{n}}{\text{Min}} &\left\{ {\alpha}\sum\limits_{i = 1}^{n} {\int_{{\Omega_{i}}} {{{\left({f - {u_{i}}} \right)}^{2}}{\chi_{i}}dx}} \right.\\ &\left. + \gamma \sum\limits_{i = 1}^{n} {\int_{{\Omega_{i}}} {\left| {\nabla {\chi_{i}}} \right|dx}} \right\}. \end{aligned} $$
(9)

In order to fulfill the condition (2), χi(x) must satisfy the following simplex constraint

$$ \sum\limits_{i = 1}^{n} {{\chi_{i}}\left(x \right)} = 1. $$
(10)

Using variational level set method, each region is described by one-level set function ϕi(x), then, its characteristic function can be denoted by the Heaviside function of ϕi(x) (i.e., χi(x)=H(ϕi(x))). Replacing χi(x) with H(ϕi(x)) in Eqs. (9) and (10) leads to the Potts model under variational level set framework.

Motivated by the Four-Color Theorem, Vese and Chan [9] proposed a strategy to partition n regions using m-level set functions with n=2m, which can avoid the problem of vacuum and overlaps, and where the model is expressed as

$$ \begin{aligned} \underset{u,\phi}{\text{Min}} \Bigg\{& E\left({u,\phi} \right) = {\alpha}\sum\limits_{i = 1}^{{2^{m}}} {\int_{\Omega} {{{\left({f - {u_{i}}} \right)}^{2}}{\chi_{i}}dx}} \\ & + {\gamma}\sum\limits_{j = 1}^{m} {\int_{\Omega} {\left| {\nabla H\left({{\phi_{j}}} \right)} \right|} dx} \Bigg\}, \end{aligned} $$
(11)

where χi is characteristic function with the form

$$ {\chi_{i}} = \prod\limits_{j = 1}^{m} {\Psi ({\phi_{j}})}, \Psi ({\phi_{j}}) = H({\phi_{j}})\ or \ 1{{ - }}H({\phi_{j}}). $$
(12)

For instance, if m=2,χ1=H(ϕ1)H(ϕ2),χ2=H(ϕ1)(1−H(ϕ2)),χ3=(1−H(ϕ1))H(ϕ2),χ4=(1−H(ϕ1))(1−H(ϕ2)).

2.2 Variational convex model for multiphase image segmentation

Taking into account the relationship of level set function and characteristic function, [11] rewrote the Chan-Vese model via a binary label function ϕ∈{0,1}

$$ \begin{aligned} \underset{{u_{1}},{u_{2}},\phi \in \left\{ {0,1} \right\}}{\text{Min}} \Bigg\{& E\left({{u_{1}},{u_{2}},\phi} \right) = {\alpha_{1}}\int_{\Omega} {{{\left({f - {u_{1}}} \right)}^{2}}\phi} dx \\ & + {\alpha_{2}}\int_{\Omega} {{{\left({f - {u_{2}}} \right)}^{2}}\left({1 - \phi} \right)dx} \\ & + \left. {\gamma \int_{\Omega} {\left| {\nabla \phi} \right|} dx} \right\} \Bigg.. \end{aligned} $$
(13)

During automatic computation, ϕ is relaxed to a convex version ϕ∈[0,1], which is recovered to its binary label function finally. For the case of multiphase image segmentation, [13] introduced m binary label functions ϕj∈{0,1},(j=1,...,m) to define n=2m characteristic functions for region partitioning, leading to the new Vese-Chan model in convex form

$$ \begin{aligned} \underset{u,\phi}{\text{Min}} \Bigg\{& E\left({u,\phi} \right) = {\alpha}\sum\limits_{i = 1}^{{2^{m}}} {\int_{\Omega} {{{\left({f - {u_{i}}} \right)}^{2}}{\chi_{i}}dx}} \\ & + {\gamma}\sum\limits_{j = 1}^{m} {\int_{\Omega} {\left| {\nabla {\phi_{j}}} \right|} dx} \Bigg\}, \end{aligned} $$
(14)

where the form of characteristic function is

$$ {\chi_{i}} = \prod\limits_{j = 1}^{m} {\Psi ({\phi_{j}})}, \quad \Psi ({\phi_{j}}) = {\phi_{j}}\ or \ 1{{ - }}{\phi_{j}}. $$
(15)

Whenever the number of image regions is taken as N∈(2m−1,2m), the problem of the automatic estimation of redundant parameters ui arises, and this in turn increases the computation cost due to inherent redundancy. Assuming that the number of redundant regions start from N+1, Li et al. [17] proposed a scheme to avoid redundant parameter estimation using various characteristic function formulas. Their energy functional minimization problem is as follows

$$ \begin{aligned} \underset{u,\phi}{\text{Min}} \Bigg\{& E\left({u,\phi} \right) = \alpha \sum\limits_{i = 1}^{N} {\int_{\Omega} {{{\left({f - {u_{i}}} \right)}^{2}}\chi_{i}^{N}dx}} \\ & + \sum\limits_{j = 1}^{m} {\int_{\Omega} {\left| {\nabla {\phi_{j}}} \right|} dx} \Bigg\}, \end{aligned} $$
(16)

where \({u = \left \{ {{u_{i}}} \right \}_{i = 1}^{N}}, {\phi = \left \{ {{\phi _{j}}} \right \}_{j = 1}^{m}}, {\chi _{i}^{N}}\) is characteristic function of the ith region. If \(N = {2^{m}}, {\chi _{i}^{N}}\) is defined evenly

$$ \chi_{i}^{{2^{m}}} = {\left({ - 1} \right)^{s_{i}^{m}}} \overset{m}{\underset{j = 1}{\Pi}} \left({{\phi_{j}} - b_{j}^{m,i - 1}} \right), $$
(17)

if 2m−1<N<2m,

$$ \left\{ \begin{array}{l} \chi_{i}^{N} = {\left({ - 1} \right)^{s_{i}^{m}}} \overset{m}{\underset{j = 1}{\Pi}} \left({{\phi_{j}} - b_{j}^{m,i - 1}} \right)\\ \quad for \ i = 1, \cdots,2{{i_{0}}} \\ \chi_{i}^{N} = {\left({ - 1} \right)^{s_{{i_{1}}}^{{m_{1}}}}} \overset{m_{1}}{\underset{j = 1}{\Pi}} \left({{\phi_{j}} - b_{j}^{{m_{1}},{i_{1}} - 1}} \right)\\ \quad for \ i = 2{{i_{0}}} + 1, \cdots,N \\ \end{array} \right., $$
(18)

where \({m_{1}}=m - 1, {{i_{0}}} = N - {2^{{m_{1}}}}, {{i_{1}} = i - {i_{0}}}, s_{i}^{m} = \sum \limits _{j = 1}^{m} {b_{j}^{m,i - 1}}, b_{j}^{m,i - 1} = \{ 0,1\} \). This scheme can handle different situations for different numbers of image regions. However, it destroys the original symmetric forms with complicated different characteristic functions. For example, there are different forms of characteristic function for N=5 and N=7, which are as follows

$$\begin{array}{*{20}l} & \left\{\begin{array}{l} {\chi_{1}^{5} =\phi_{1} \phi_{2} \phi_{3},} \\ { \chi_{2}^{5}=\phi_{1} \phi_{2}\left(1-\phi_{3}\right)} \\ {\chi_{3}^{5}=\phi_{1}\left(1-\phi_{2}\right)}\\ {\chi_{4}^{5}=\left(1-\phi_{1}\right) \phi_{2}} \\ {\chi_{5}^{5}=\left(1-\phi_{1}\right)\left(1-\phi_{2}\right)} \end{array}\right., \end{array} $$
(19a)
$$\begin{array}{*{20}l} & \left\{ \begin{array}{l} {\chi_{1}^{7} =\phi_{1} \phi_{2} \phi_{3}} \\ { \chi_{2}^{7}=\phi_{1} \phi_{2}\left(1-\phi_{3}\right)} \\ {\chi_{3}^{7}=\phi_{1}\left(1-\phi_{2}\right) \phi_{3}}\\ {\chi_{4}^{7}=\phi_{1}\left(1-\phi_{2}\right)\left(1-\phi_{3}\right)} \\ {\chi_{5}^{7}=\left(1-\phi_{1}\right) \phi_{2} \phi_{3}}\\ { \chi_{6}^{7}=\left(1-\phi_{1}\right) \phi_{2}\left(1-\phi_{3}\right)} \\ {\chi_{7}^{7}=\left(1-\phi_{1}\right)\left(1-\phi_{2}\right)} \end{array}\right., \end{array} $$
(19b)

which will reduce the generality of characteristic function for different region numbers and increase the computation time. In the next section, we propose another modified Vese-Chan model to overcome these problems.

3 Methods

3.1 Characteristic functions based on binary decomposition

According to the original scheme by [9], the intersection of m binary functions can generate n=2m characteristic functions for region partition. Let the binary expression of characteristic functions χi be \(b_{1}^{i - 1}b_{2}^{i - 1}\cdots b_{m}^{i - 1}\), where the value of \(b_{j}^{i - 1}\) is related to the binary expression about the subtraction of index by one (i.e., i−1) for the ith region Ωi, then the characteristic function of Ωi is

$$ {\chi_{i}} = \prod\limits_{j = 1}^{m} \left({b_{j}^{i - 1} + {{\left({ - 1} \right)}^{b_{j}^{i - 1}}}{\phi_{j}}} \right), $$
(20)

where ϕj∈{0,1}, \(b_{j}^{i - 1} = \left \{ {0,1} \right \}, i=1,\cdots,2^{m}, {j = 1,\cdots,m}\). Taking m=3 for example, 8 regions are partitioned as shown in Fig. 1, and the relevant characteristic functions are listed in Table 1.

Fig. 1
figure 1

Eight regions partitioned by three binary label functions

Full size image
Table 1 Relationship between binary region number and characteristic function
Full size table

3.2 The classic optimization scheme of the Vese-Chan model

Equation (14) stated a multiphase segmentation Vese-Chan model using m binary label functions to define 2m characteristic functions, which is usually solved by alternative optimization method.

The idea of alternative optimization is to solve minimum problem of a variable by fixing others. Firstly fix ϕ to optimize u, we obtain

$$ u_{i}=\frac{\int_{\Omega} f \chi_{i} d x}{\int_{\Omega} \chi_{i} d x}, (i=1,\cdots,2^{m}). $$
(21)

When u is fixed, we compute ϕ by gradient descent method, which is shown below

$$ \frac{\partial \phi_{j}}{\partial t}=\gamma \nabla \cdot\left(\frac{\nabla \phi_{j}}{\left|\nabla \phi_{j}\right|}\right)-\alpha \sum_{i=1}^{2^{m}} \left(f-u_{i}\right)^{2} \frac{\partial \chi_{i}}{\partial \phi}, $$
(22)

where \(\nabla \cdot \left (\frac {\nabla \phi }{\left |\nabla \phi \right |}\right)\) is the curvature term of ϕ. We use the upwind differential scheme [16] to compute this term

$$ \begin{aligned} & \nabla \cdot\left(\frac{\nabla \phi_{x, y}}{\left|\nabla \phi_{x, y}\right|}\right) \\ & =\nabla^{+} \cdot\left(\frac{\nabla^{-} \phi_{x, y}}{ | \nabla^{-} \phi_{x, y}}\right) \\ & =d_{1 x, y} \phi_{x+1, y}+d_{2 x, y} \phi_{x-1, y} \\ & +d_{3 x, y} \phi_{x, y+1}+d_{4 x, y} \phi_{x, y-1}-D_{x, y} \phi_{x, y} \end{aligned} $$
(23)

where

$$ \begin{aligned} & {d_{1 x, y} =\frac{1}{\sqrt{\left(\frac{\phi_{x+1, y}-\phi_{x, y}}{h}\right)^{2}+\left(\frac{\phi_{x,y+1}-\phi_{x,y-1}}{2 h}\right)^{2}+\epsilon}}} \\ & {d_{2 x, y} =\frac{1}{\sqrt{\left(\frac{\phi_{x, y}-\phi_{x-1, y}}{h}\right)^{2}+\left(\frac{\phi_{x-1, y+1}-\phi_{x-1, y-1}}{2 h}\right)^{2}+\epsilon}}} \\ & {d_{3 x, y} =\frac{1}{\sqrt{\left(\frac{\phi_{x, y+1}-\phi_{x, y}}{h}\right)^{2}+\left(\frac{\phi_{x+1, y}-\phi_{x-1, y}}{2 h}\right)^{2}+\epsilon}}} \\ & {d_{4 x, y} =\frac{1}{\sqrt{\left(\frac{\phi_{x, y}-\phi_{x, y-1}}{h}\right)^{2}+\left(\frac{\phi_{x+1, y-1}-\phi_{x-1, y-1}}{2 h}\right)^{2}+\epsilon}}} \\ & {D_{x y} =d_{1 x, y}+d_{2 x, y}+d_{3 x, y}+d_{4 x, y}}. \end{aligned} $$
(24)

In Eq. (24), x and y represent coordinates in the image, h is spatial step, and ε is a small positive value.

For the optimization of variable ϕ, the curvature term of ϕ will be obtained if the alternative optimization method is used directly, causing complex difference scheme and low computational efficiency. However, the ADMM method can avoid computing the curvature term and simplify the solving process by introducing auxiliary variables and Lagrange multipliers.

3.3 The modified model for gray images and its ADMM

For a gray value image f(x):Ω→R with N<2m regions, there exist 2m−N redundant regions in the range [N+1,2m]. To avoid the redundant parameter estimations, we add area constraints on these regions

$$ \int_{\Omega} {{\chi_{l}}\left({\phi} \right)dx} = 0, \left({l = N + 1,...,{2^{m}}} \right). $$
(25)

We cast the Vese-Chan model into

$$ \begin{aligned} \underset{u,\phi}{\text{Min}} \Bigg\{& {E_{1}}\left({u,\phi} \right) = {\alpha}\sum\limits_{i = 1}^{N} {\int_{\Omega} {{{\left({f - {u_{i}}} \right)}^{2}}{\chi_{i}}\left(\phi \right)dx}} \\ &+ {\gamma}\sum\limits_{j = 1}^{m} {\int_{\Omega} {\left| {\nabla {\phi_{j}}} \right|} dx} \Bigg\} \end{aligned} $$
(26)

with constraints (25), where \({u = \left \{ {{u_{i}}} \right \}_{i = 1}^{N}}\) consists of parameters need to be estimated, \({\phi = \left \{ {{\phi _{j}}} \right \}_{j = 1}^{m}}\) consists of binary label functions.

In order to improve the efficiency, ADMM is used in our method. We introduce auxiliary variable w=∇ϕ, Lagrange multipliers λj,λl and positive penalty parameters θ,μ. We add constraints (25) into Eq. (26), which can be rewritten as the following form of alternative optimization

$$ \begin{aligned} \left({u^{k + 1},\phi^{k + 1},{\boldsymbol{w}}^{k + 1}} \right) = \arg \underset{u,\phi,{\boldsymbol{w}}}{\text{Min}} \\ \left\{ \begin{array}{ll} {E_{1}}\left({u,{\phi},{\boldsymbol{w}}} \right) &= {\alpha}\sum\limits_{i = 1}^{N} {\int_{\Omega} {{{\left({f - {u_{i}}} \right)}^{2}}{\chi_{i}}\left(\phi \right)} dx} \\ &+ {\gamma}\sum\limits_{j = 1}^{m} {\int_{\Omega} {\left| {{{\boldsymbol{w}}_{j}}} \right|} dx} \\ &+ \sum\limits_{j = 1}^{m} {\int_{\Omega} {{\boldsymbol{\lambda }}_{j}^{k} \cdot ({{\boldsymbol{w}}_{j}} - \nabla {\phi_{j}})dx}} \\ &+ {\frac{\theta}{2}\sum\limits_{j = 1}^{m} \int_{\Omega} {{{\left| {{{\boldsymbol{w}}_{j}} - \nabla {\phi_{j}}} \right|}^{2}}dx}} \\ &+ \sum\limits_{l = N + 1}^{{2^{m}}} {\int_{\Omega} {\lambda_{l}^{k}{\chi_{l}}\left(\phi \right)dx}} \\ &+ {\frac{\mu}{2}\sum\limits_{l = N + 1}^{{2^{m}}} \int_{\Omega} {{{\left({{\chi_{l}}\left(\phi \right)} \right)}^{2}}dx}} \\ \end{array} \right\}, \end{aligned} $$
(27)

where k is the iteration step.

Using the idea of alternative optimization, Eq. (27) is transformed into the following sub-problems of optimization in each loop

$$\begin{array}{*{20}l} \left({u^{k + 1}} \right) &= \arg \underset{u}{\text{Min}} {E_{1}}\left({u,\phi^{k},{\boldsymbol{w}}^{k}} \right), \end{array} $$
(28a)
$$\begin{array}{*{20}l} \left({\phi^{k + 1}} \right) &= \arg \underset{\phi}{\text{Min}} {E_{1}}\left({u^{k + 1},{\phi},{\boldsymbol{w}}^{k}} \right), \end{array} $$
(28b)
$$\begin{array}{*{20}l} \left({{\boldsymbol{w}}^{k + 1}} \right) &= \arg \underset{\boldsymbol{w}}{\text{Min}} {E_{1}}\left({u^{k + 1},\phi^{k + 1},{{\boldsymbol{w}}}} \right). \end{array} $$
(28c)

Using variational method to Eq. (28a), fix ϕk and wk to optimize uk+1, we obtain

$$ u_{i}^{k + 1} = \frac{{\int_{\Omega} {f{\chi_{i}}\left(\phi \right)} dx}}{{\int_{\Omega} {{\chi_{i}}\left(\phi \right)} dx}}, (i = 1, \cdots,N). $$
(29)

For the problem (28b), since ϕ∈{0,1} is a non-convex function, we first transform ϕ∈{0,1} to ϕ∈[0,1] with convex relaxation and then fix uk+1 and wk, and the sub-problem of optimization with respect to ϕk+1 is as follows

$$ {{}\begin{aligned} &\left({\phi^{k + 1}} \right) = \arg {\underset{{\phi} \in \left[ {0,1} \right]}{\text{Min}}} \\ &\left\{ \begin{array}{ll} {E_{1}}\left({u^{k + 1},{\phi},{\boldsymbol{w}}^{k}} \right) &= {\alpha}\sum\limits_{i = 1}^{N} {\int_{\Omega} {{{\left({f - u_{i}^{k + 1}} \right)}^{2}}{\chi_{i}}\left(\phi \right)} dx} \\ &+ \sum\limits_{j = 1}^{m} {\int_{\Omega} {{\boldsymbol{\lambda }}_{j}^{k} \cdot ({{\boldsymbol{w}}_{j}^{k}} - \nabla {\phi_{j}})dx}} \\ &+ {\frac{\theta}{2}\sum\limits_{j = 1}^{m} \int_{\Omega} {{{\left| {{{\boldsymbol{w}}_{j}^{k}} - \nabla {\phi_{j}}} \right|}^{2}}dx}} \\ &+ \sum\limits_{l = N + 1}^{{2^{m}}} {\int_{\Omega} {\lambda_{l}^{k}{\chi_{l}}\left(\phi \right)dx}} \\ &+ {\frac{\mu}{2}\sum\limits_{l = N + 1}^{{2^{m}}} \int_{\Omega} {{{\left({{\chi_{l}}\left(\phi \right)} \right)}^{2}}dx}} \\ \end{array} \right\}, \\ \end{aligned}} $$
(30)

we obtain the Euler-Lagrange equation of \({\phi ^{k + 1}_{j}}\)

$$ \left\{ {\begin{array}{*{20}{l}} {A + B + C = 0} & {\quad x \in \Omega} \\ { - \left({{\boldsymbol{\lambda }}_{j}^{k} + {\boldsymbol{w}}_{j}^{k} - \nabla {\phi_{j}}} \right) \cdot {\boldsymbol{n}} = 0} & {\quad x \in \partial \Omega} \\ \end{array}} \right. \\ $$
(31)

where \(A=\sum \limits _{i = 1}^{N} {{\alpha }{{\left ({f - u_{i}^{k + 1}} \right)}^{2}}\frac {{\partial {\chi _{i}}\left (\phi \right)}}{{\partial \phi }}}, B=\sum \limits _{j = 1}^{m} \left (\nabla \cdot {\boldsymbol {\lambda }}_{j}^{k} + {\theta }\nabla \cdot ({\boldsymbol {w}}_{j}^{k} - \nabla {\phi _{j}}) \right), C = \sum \limits _{l = N + 1}^{{2^{m}}} {\left ({\lambda _{l}^{k} + {\mu }{\chi _{l}}\left (\phi \right)} \right)\frac {{\partial {\chi _{l}}\left (\phi \right)}}{{\partial \phi }}}\). \({\phi _{j}^{k + 1}}\) can be solved by Gauss-Seidel iterative method, then

$$ \phi_{j}^{k + 1} = \text{Max}\left(\text{Min}\left({\phi_{\mathrm{j}}^{\mathrm{k} + 1},1} \right),0 \right). $$
(32)

For Eq. (28c), to fix uk+1 and ϕk+1 to optimize wk+1, we obtain the Euler-Lagrange equation

$$ {\gamma}\frac{{{{\boldsymbol{w}}_{j}}}}{{\left| {{{\boldsymbol{w}}_{j}}} \right|}} + {\boldsymbol{\lambda }}_{j}^{k} + {\theta}({{\boldsymbol{w}}_{j}} - \nabla \phi_{j}^{k + 1}) = 0, $$
(33)

when |wj|≠0. Equation (33) can be expressed as a following generalized soft thresholding formula in analytical form

$$ {\boldsymbol{w}}_{j}^{k + 1} = \max \left({\left| {\nabla \phi_{j}^{k + 1} - \frac{{{\boldsymbol{\lambda }}_{j}^{k}}}{{{\theta}}}} \right| - \frac{{{\gamma }}}{{{\theta }}},0} \right)\frac{{\nabla \phi_{j}^{k + 1} - \frac{{{\boldsymbol{\lambda }}_{j}^{k}}}{{{\theta}}}}}{{{{\left| {\nabla \phi_{j}^{k + 1} - \frac{{{\boldsymbol{\lambda }}_{j}^{k}}}{{{\theta}}}} \right|}_{\epsilon} }}}, $$
(34)

The Lagrange multipliers \({\boldsymbol {\lambda }}_{j}^{k + 1}\) and \(\lambda _{l}^{k + 1}\) can be updated as the following

$$ \left\{ \begin{array}{l} {\boldsymbol{\lambda}}_{j}^{k + 1} = {\boldsymbol{\lambda}}_{j}^{k} + {\theta}\left({\boldsymbol{w}}_{j}^{k + 1} - \nabla \phi_{j}^{k + 1}\right) \\\\ \lambda_{l}^{k + 1} = \lambda_{l}^{k} + {\mu}{\chi_{l}}\left(\phi \right) \end{array} \right.. $$
(35)

In order to represent the boundary of segmented images, it is necessary to transform ϕj after convex relaxation to binary label function as the following

$$ {\phi_{j}} = \left\{ {\begin{array}{ll} 1 & \quad {{\phi_{j}} \ge \tau} \\ 0 & \quad {{\phi_{j}} < \tau} \\ \end{array}} \right., $$
(36)

where τ∈(0,1). The parameter τ is usually chosen as \(\tau = \frac {1}{2}\). The stopping criterion is based on the relative energy error formula |E1k+1−E1k|E1k≤ε, where ε is a small positive value.

The ADMM of modified model for gray images can be described as Algorithm 1.

3.4 The modified model for color images and its ADMM

Different from a scalar image, a color image consists of three layers, so parameter estimation is needed for each layer in different regions. Let f=(f1,f2,f3)∈R3 be a color image, ui=(ui1,ui2,ui3) be the piecewise constant parameters to be estimated in the ith region, ϕ consists of m binary functions as the previous sub-section, and the variational model for multiphase color image segmentation is

$$ \begin{aligned} \underset{{\boldsymbol{u}},\phi}{\text{Min}} \Bigg\{& {E_{2}}\left({{\boldsymbol{u}},\phi} \right) = {\alpha}\sum\limits_{i = 1}^{N} {\sum\limits_{p = 1}^{3} {\int_{\Omega} {{{\left({{f_{p}} - {u_{ip}}} \right)}^{2}}{\chi_{i}}\left(\phi \right)dx}} } \\ & + {\gamma}\sum\limits_{j = 1}^{m} {\int_{\Omega} {\left| {\nabla {\phi_{j}}} \right|} dx} \Bigg\}, \end{aligned} $$
(37)

it is also subjected to the constraints (25).

To design its ADMM method, we introduce auxiliary variable w=∇ϕ, Lagrange multipliers λj,λl and positive penalty parameters θ,μ. Then we add constraints (25) into Eq. (37) and transform it into the following form of alternative optimization

$$ \begin{aligned} &\left({\boldsymbol{u}^{k + 1},\phi^{k + 1},{\boldsymbol{w}}^{k + 1}} \right) = \arg \underset{{\boldsymbol{u}},\phi,{\boldsymbol{w}}}{\text{Min}} \\ &\left\{ \begin{array}{ll} {E_{2}}\left({\boldsymbol{u},{\phi},{{\boldsymbol{w}}}} \right) & = {\alpha}\sum\limits_{i = 1}^{N} {\sum\limits_{p = 1}^{3} {\int_{\Omega} {{{\left({{f_{p}} - {u_{ip}}} \right)}^{2}}{\chi_{i}}\left(\phi \right)} dx}} \\ &+ {\gamma}\sum\limits_{j = 1}^{m} {\int_{\Omega} {\left| {{{\boldsymbol{w}}_{j}}} \right|} dx} \\ &+ \sum\limits_{j = 1}^{m} {\int_{\Omega} {{\boldsymbol{\lambda }}_{j}^{k} \cdot ({{\boldsymbol{w}}_{j}} - \nabla {\phi_{j}})dx}} \\ &+ {\frac{\theta}{2}\sum\limits_{j = 1}^{m} \int_{\Omega} {{{\left| {{{\boldsymbol{w}}_{j}} - \nabla {\phi_{j}}} \right|}^{2}}dx}} \\ &+ \sum\limits_{l = N + 1}^{{2^{m}}} {\int_{\Omega} {\lambda_{l}^{k}{\chi_{l}}\left(\phi \right)dx}} \\ &+ {\frac{\mu}{2}\sum\limits_{l = N + 1}^{{2^{m}}} \int_{\Omega} {{{\left({{\chi_{l}}\left(\phi \right)} \right)}^{2}}dx}} \\ \end{array} \right\}. \\ \end{aligned} $$
(38)

Using the same procedure as the last sub-section in each loop of optimization, we get successively

$$ u_{ip}^{k + 1} = \frac{{\int_{\Omega} {{f_{p}}{\chi_{i}}\left(\phi \right)} dx}}{{\int_{\Omega} {{\chi_{i}}\left(\phi \right)} dx}}. $$
(39)

Fixing \({u_{ip}^{k + 1}}\) and wk, we obtain the Euler-Lagrange equation of ϕk+1 as the following

$$ \left\{ {\begin{array}{ll} {A_{i} + B + C = 0} & {x \in \Omega} \\ { - \left({{\boldsymbol{\lambda }}_{j}^{k} + {\boldsymbol{w}}_{j}^{k} - \nabla {\phi_{j}}} \right) \cdot {\boldsymbol{n}} = 0 \quad } & {x \in \partial \Omega} \\ \end{array}} \right., $$
(40)

where \(A_{i}={\alpha }\sum \limits _{i = 1}^{N} {\sum \limits _{p = 1}^{3} {{\left ({{f_{p}} - u_{ip}^{k + 1}} \right)^{2}}\frac {{\partial {\chi _{i}}\left (\phi \right)}}{{\partial \phi }}} }\). We design Gauss-Seidel iterative method from Eq. (40) to compute \({\phi _{j}^{k + 1}}\).

\({\boldsymbol {w}}_{j}^{k + 1}\) can be solved from Eq. (34). The Lagrange multipliers \({\boldsymbol {\lambda }}_{j}^{k + 1}\) and \(\lambda _{l}^{k + 1}\) can be updated as Eq. (35).

The ADMM of the modified model for color images can be described as Algorithm 2.

4 Results and discussion

In this section, we present two groups of numerical experiments for segmented images including gray images and color images to study the effectiveness and efficiency of our method. For Eqs. (27) and (38), if the number of regions N=2m, this method represents traditional 2m-phase segmentation with ADMM. If N<2m, it represents the proposed method in this paper. The initial binary label function ϕ is initialized as m circles, there are two cases, in the circle ϕ0=1 otherwise ϕ0=0. In our experiments, We set α=1,μ=1, h=1,ε=10−6, ε=10−5. All experiments are performed using MATLAB R2016a on a Window 7 platform with an Intel(R) Core(TM) i5-4590 CPU at 3.30GHz and 4.00GB memory.

4.1 Numerical experiments for gray image segmentation

Firstly, we compare the Vese-Chan model without redundant regions by using classic optimization method shown in Section 3.2 and ADMM method in Section 3.3 to study advantages of the latter method. Besides, we compare Algorithm 1 with the original Vese-Chan model in [9] and TV regularization method in [17] for images with redundant regions.

4.1.1 Comparisons between the classic optimization method and ADMM

Figure 2 a shows a synthetic image with four phases, so we use two binary label functions to design characteristic functions. Figure 2b and c are segmentation results for classic optimization method and ADMM method, respectively, which contains four meaningful phases (Fig. 2d–g, h–k). We conclude that the above two methods obtain the same segmentation effectiveness.

Fig. 2
figure 2

Comparison of the classic optimization method and ADMM. a The input synthetic image. b The segmentation result of classic optimization method. c The segmentation result of ADMM. d–g Different regions obtained via the classic optimization method. h–k Different regions obtained via ADMM

Full size image

In Fig. 3, we compare the classic optimization method and ADMM method for a brain MR image. The four phases of two methods are displayed in the second row and the third row, respectively. The ADMM method can give the better segmentation result from the comparison of Fig. 3f and g.

Fig. 3
figure 3

Four-phase segmentation of the classic optimization method and ADMM. a The input brain MR image. b The segmentation result of classic optimization method. c The segmentation result of ADMM. d–g Different regions obtained via the classic optimization method. h–k Different regions obtained via ADMM

Full size image

Table 2 demonstrates the iterations and computational time for Figs. 2a and 3a. ADMM method consumes less computational time and iterations than the classic optimization method. We can conclude that using ADMM method can improve computational efficiency for image segmentation.

Table 2 Iterations and computational time for Figs. 2a and 3a
Full size table

4.1.2 Comparisons between Algorithm 1 and other methods

Figure 4 a shows a picture consists of a little bottle, a cup, and background, which occupy three regions. Using Vese-Chan scheme for division, we need two binary label functions to design characteristic functions. The original Vese-Chan method uses N=4 with one empty region, TV regularization method and our method use N=3 without empty regions. Three methods get the same segmentation results as shown in Fig. 4b–d, which contain three meaningful regions (Fig. 4e–g, i–k, and m–o) and one empty region (Fig. 4h, l, and p). The differences of their computation efficiency are listed lately together with other experiments.

Fig. 4
figure 4

Comparison of original Vese-Chan method, TV regularization method, and our method using two label functions. a The input nature image. b The segmentation of original Vese-Chan method. c The segmentation of TV regularization method. d The segmentation of our method. e–h Different regions obtained via the original Vese-Chan method. i–l Different regions obtained via TV regularization method. m–p Different regions obtained via our method

Full size image

Figure 5 presents a three-phase segmentation problem also, but the image is a synthetic one with noises as shown in Fig. 5a. Figure 5 b–d show the segmentation results via the original Vese-Chan method, TV regularization method, and our method, respectively. Comparison between four characteristic functions of every method (Fig. 5e–p) shows that three methods obtain the same segmentation result.

Fig. 5
figure 5

Three-phase segmentation of original Vese-Chan method, TV regularization method, and our method. a The input synthetic image. b The segmentation of original Vese-Chan method. c The segmentation of TV regularization method. d The segmentation of our method. e–h Different regions obtained via the original Vese-Chan method. i–l Different regions obtained via TV regularization method. m–p Different regions obtained via our method

Full size image

In Fig. 6, we compare the original Vese-Chan method, TV regularization method, and our method for a three-phase CT image slice of brain, and the segmentation results using these three methods are given in Fig. 6b–d, respectively. Different regions of three methods are shown in the last three rows. From experiment results, we conclude that three methods can obtain the similar segmentation effectiveness.

Fig. 6
figure 6

Three-phase segmentation for a CT image slice of brain. a The input CT image slice of brain. b The segmentation of original Vese-Chan method. c The segmentation of TV regularization method. d The segmentation of our method. e–h Different regions obtained via the original Vese-Chan method. i–l Different regions obtained via TV regularization method. m–p Different regions obtained via our method

Full size image

Figure 7 a shows a synthetic image consists of six regions, and we need three binary label functions to design characteristic functions. The original Vese-Chan method uses N=8 with two empty regions, TV regularization method, and our model use N=6 without empty regions. Three methods get the same segmentation results as shown in Fig. 7b–d, which can obtain the same meaningful phases. Therefore, we only show six different regions obtained via our method in Fig. 7e–j.

Fig. 7
figure 7

Six-phase segmentation of original Vese-Chan method, TV regularization method and our method. a The input synthetic image. b The segmentation of original Vese-Chan method. c The segmentation of TV regularization method. d The segmentation of our method. e–l Different regions obtained via our method

Full size image

We record iterations and computational time of three methods for Figs. 4a, 5a, 6a, and 7a in Table 3. From Table 3, we reach that TV regularization method and our method both consume less computational time than original Vese-Chan method, and the reason of which is the first two methods reduce estimation of redundant parameters. Our method can obtain the highest computational efficiency. Next, we will study the reasons for this result.

Table 3 Iterations and computational time for Figs. 4a, 5a, 6a, and 7a
Full size table

TV regularization method uses an acceleration method, and our method also uses ADMM acceleration method. In order to further explore the efficiency of our method, we compare four methods, including TV regularization without acceleration method (TV-RWAM), our method without ADMM, TV regularization method, and our method. The first two methods use upwind differential scheme as shown in Section 3.2 to compute label function ϕ. Figure 8 presents a six-phase segmentation problem, the image is a remote sensing one of coastline as shown in Fig. 8a. Figure 8 b–e show the segmentation results via the above four methods, respectively. Figure 8f–i are gray images according to four final segmentations. The iterations and computational time of four methods are showed in Table 4.

Fig. 8
figure 8

Comparison of TV regularization method and our method both solved without/with acceleration method. a The remote sensing image of coastline. b–e The segmentation of TV-RWAM, our method without ADMM, TV regularization method and our method. f–i Gray images obtained via the final segmentation results

Full size image
Table 4 Iterations and computational time for Fig. 8a
Full size table

From Table 4, our method without ADMM consumes less iterations and computational time than TV regularization without acceleration method, because we add area constraints of redundant regions and simplify the expression of characteristic functions. Besides, our method consumes less computational time than TV regularization method, and the reason of which is we use ADMM method to compute energy optimization, improving the computational efficiency.

4.2 Numerical experiments for color image segmentation

Color image segmentation is an extension of gray image segmentation. We compare the original Vese-Chan method and TV regularization method with our method presented in Section 3.4 for color images.

Figure 9 presents a three-phase segmentation problem, the color image is a part of flower as shown in Fig. 9a. Figure 9b–d show the segmentation results via the above three methods respectively. Figure 9e–h show different regions via the original Vese-Chan method. Because the TV regularization method and our method obtain the same result, we only show different regions about our method (Fig. 9i–l). Comparison between empty regions (Fig. 9h and l) shows that our method and TV regularization method give better segmentation result.

Fig. 9
figure 9

Three-phase segmentation of a color image. a The input color image. b The segmentation of original Vese-Chan method. c The segmentation of TV regularization method. d The segmentation of our model. e–h Different regions obtained via the original Vese-Chan method. i–l Different regions obtained via our method

Full size image

Figure 10a presents a six-phase synthetic color image, the segmentation results using the original Vese-Chan method, TV regularization method and our method are given in Fig. 10b–d, respectively, which obtain the same meaningful regions. Therefore we only show six meaningful regions via our method in Fig. 10e–j.

Fig. 10
figure 10

Six-phase segmentation of a color image. a The input color synthetic image. b The segmentation of original Vese-Chan method. c The segmentation of TV regularization method. d The segmentation of our method. e–l Different regions obtained via our method

Full size image

Figures 11a, 12a, and 13a are synthetic color images with nine, ten and eleven different regions respectively, so we need four binary label functions to design characteristic functions. We use three methods: the original Vese-Chan method, TV regularization method and the proposed method. The first one uses N=16 for three images, causing seven, six and five empty regions separately. Our method uses N=9,N=10 and N=11 without empty regions respectively. For Fig. 11a, results are shown in Fig. 11b–d. Three methods obtain the same segmentation results. As a representative, Fig. 11e–m show nine different regions obtained via our method. Figure 12b–d show segmentation results for a ten-phase image. The first eight different regions for Fig. 12 are the same as Fig. 11, and two other regions as shown in Fig. 12e–f. The segmentation results for Fig. 13a are shown in Fig. 13b–d, and its 11 different regions include Figs. 11e–l, 12e, and 13e–f.

Fig. 11
figure 11

Nine-phase segmentation of a color image. a The input color synthetic image. b The segmentation of original Vese-Chan method. c The segmentation of TV regularization method. d The segmentation of our model. e–m Nine meaningful regions obtained via our method

Full size image
Fig. 12
figure 12

Ten-phase segmentation of a color image. a The input color synthetic image. b The segmentation of original Vese-Chan method. c The segmentation of TV regularization method. d The segmentation of our model. e–f Two meaningful regions obtained via our method

Full size image
Fig. 13
figure 13

Eleven-phase segmentation of a color image. a The input color synthetic image. b The segmentation of original Vese-Chan method. c The segmentation of TV regularization method. d The segmentation of our model. e–f Two meaningful regions obtained via our method

Full size image

We record total iterations and computational time of the above three methods for Figs. 9a, 10a, 11a, 12a, and 13a in Table 5, from which we can conclude that compared to the original Vese-Chan method and TV regularization method, and our method can improve computational efficiency for color images. With increase of the empty regions of images, it is more obvious that the computational efficiency is improved. Besides, the iterations and computational time are positively correlated with the image sizes.

Table 5 Iterations and computational time for Figs. 9a, 10a, 11a, 12a, and 13a
Full size table

From the above experiments, we can conclude that the computational cost of our method includes two parts: parameter computations for ϕ,w,λ and parameter estimations for u, both of which are relative to the image sizes. Besides, the former one is also relative to the number of binary label function m. When there exist redundant regions in the images, (i.e., N<2m), parameter estimations are relevant to N.

The number of regions N can be detected according to the number of target objects in the image. In the paper, the segmented images are simple and easy to determine the value of N. However, if the image is complex, the value of N is difficult to determine, which is a limit that needs to be solved in the future.

5 Conclusions

In order to improve the computation efficiency of the Vese-Chan model for multiphase image segmentation in a systematic form, we have designed a modified Vese-Chan model by introducing some simple area constraints. The forms of characteristic functions are unified, and the redundant parameter estimations are not needed, therefore, the cost of computation is reduced. Obviously, the computational efficiency is higher with the increase of the empty regions. Additionally, we formulate the ADMM for the proposed model to further improve efficiency. Some numerical examples for gray image and color image multiphase segmentation are presented to demonstrate that the proposed model has the same or better segmentation effects and higher efficiency. Our method can be applied into motion segmentation, surface segmentation, and 3D reconstruction in the future work.