Multi-phase image segmentation by the Allen–Cahn Chan–Vese model

In [19], we give the definition of the inhomogeneous graph Laplacian operator for a grayscale image, which can be easily extended to a color image by summation of information on different channels.

Let $\Omega$ be a 2D discrete image domain, and $I$ be an image defined on it with $M_{1}\times M_{2}$ pixels. For pixel $x_{0}=(i,j)\in\Omega$, the generalized inhomogeneous graph Laplacian operator $L$ for a color image is given by

where $I_{i,j,r}^{\ell}$ is intensity value of the $\ell$-th neighbour point of $I_{i,j}$ in the $r$-th channel. The relevant 8 neighbors are

and the weighted coefficient

which is between $0$ and $1$. $\kappa$ in (3.2) is a given non-negative parameter.

Subsequently, we choose a threshold value $\sigma$ to locate the zero-cross points of the generalized inhomogeneous Laplacian operator by condition $\left|L(x_{0})\right|>\sigma$. Finally, all zero-cross points are collected in the set $\mathcal{S}$ to represent rough initial edges.

This part is devoted to dividing the rough initial edges obtained in the last stage into different phases. We utilize the K-means clustering algorithm proposed in [22], minimizing the sum of distances from each object to its cluster centroid for all clusters. The edges of each phase have different intensities as well as the interior region. Consequently, the K-means clustering algorithm can help segment the rough initial edges into $m$ phases.

Due to the influence of noise, the denoising process is necessary. Similar to [19], the diagonal connectivity is used as a guiding principle to remove most noise points from the rough initial color contours. The edge of an objective should have connectivity, which means that points of edges connect with each other, while the noise doesn’t possess this property.

To eliminate noise in the rough initial color contours, we keep all the diagonally connected points and remove the other points. The denoising process needs to be repeated $M$ times where $M$ is a pre-setting small integer.

Now combining all above treatments, we restate the Multi-IGLIM as follows.

From the Multi-IGLIM in Algorithm 1, we get the denoising edge detection results $v_{i}^{0}$ for $m$ phases. However, the number of the initial values for the minimization process for the ACCV model (2.1) is $n=\lceil\log_{2}{m}\rceil$. In this part, $m=4,n=2$. We first sort the results of the Multi-IGLIM, $v^{0}_{i}$, in ascending order by the perimeter for each phase. Then the initial contours for two phase variables are

In the case of overlapped edges, $u^{0}_{i}$ will be beyond $1$ somewhere, for which we will force them to equal to $1$ directly. The integration method (3.3)-(3.4) is not unique, and other reasonable treatment can also be considered.

This part introduces the alternating direction minimization method (ADMM) to minimize the energy functional (2.2). The minimization procedure can be written as

with $\mathbf{u}=(u_{1},u_{2}).$

Since the Heaviside function $H(u)$ in (2.3) is discontinuous, we give $H_{\epsilon_{1}}$ as a regularization of the Heaviside function $H(u)$ in practice

Here, the value of $\epsilon_{1}$ is requested to be less than $\frac{1}{2}$. The derivative of $H_{\epsilon_{1}}(u)$ is an approximation of Dirac function $\displaystyle\delta(u)=\frac{\mathrm{d}}{\mathrm{du}}H(u)$ given by [15]

For a fixed $u_{i},i=1,2$ and $r=1,\cdots,\omega$, the variational calculation gives

For (3.5), we only need to solve the following two Allen–Cahn equations to the steady state,

with the homogeneous Neumann boundary condition. Here,

The function $w(u)=2u(2u-1)(u-1)$ is the derivative of $W(u)$ with respect to $u$.

Considering the properties of maximum boundedness and energy decay for the classical Allen–Cahn equation, we would next explore the relevant properties for the proposed Allen–Cahn equations (3.8) and (3.9). According to the criterion in [8], we can verify that the proposed Allen–Cahn equations (3.8) and (3.9) satisfy the following maximum bound principle

for homogeneous Neumann boundary condition. Besides, the energy decay property for (2.2) holds obviously. In the following part, we will design efficient numerical schemes preserving the corresponding properties, which are stable and reliable to get the steady state solution. In [7, 10], it is shown that the ETD method can maintain the maximum bound principle and energy stability for the Allen–Cahn equation. As a consequence, we will utilize the ETD method as temporal discretization to construct the numerical schemes for the proposed Allen–Cahn equations (3.8) and (3.9).

Using the central finite difference method to discretize (3.8) and (3.9) in space, we obtain the following ODE system for numerical solution $U_{i}=[U_{i,\ell}]^{M_{1}M_{2}}_{\ell=1}\in\mathbb{R}^{M_{1}M_{2}}$,

where

Here, $D_{h}$ is a 2D discrete Laplacian matrix under the homogeneous Neumann boundary condition,

where $I_{M_{i}}$ is an $M_{i}\times M_{i}$ identity matrix and

$I_{d}\in\mathbb{R}^{M_{1}M_{2}\times M_{1}M_{2}}$ is an identity matrix. The positive constant $S$ is called the stabilizer.

Multiplying (3.10) with an integrating factor $e^{L_{h}t}$, and integrating from $t^{n}$ to $t^{n+1}$ give

Then two approximations for the nonlinear term $N(U(t^{n}+s))$ are proposed.

• •

  ETD1 : Approximating $N(U_{i}(t^{n}+s))$ by $N(U_{i}(t^{n}))$, we have a first-order scheme

  $\displaystyle U_{i}^{n+1}=\phi_{0}(L_{h}\Delta t)U_{i}^{n}+\Delta t\phi_{1}(L_{h}\Delta t)N(U_{i}^{n}).$

  (3.11)

• •

  ETDRK2: Approximating $N(U_{i}(t^{n}+s))$ with $\displaystyle(1-\frac{s}{\Delta t})N(U_{i}(t^{n}))+\frac{s}{\Delta t}N(\hat{U}_{i}^{n+1})$, wherein $\hat{U}_{i}^{n+1}$ is a first-order approximation obtained from ETD1. As a consequence, the ETDRK2 scheme is given by

  $$\left\{\begin{split}&\hat{U}_{i}^{n+1}=\phi_{0}(L_{h}\Delta t)U_{i}^{n}+\Delta t\phi_{1}(L_{h}\Delta t)N(U^{n}_{i}),\\ &{U}_{i}^{n+1}=\phi_{0}(L_{h}\Delta t)\hat{U}_{i}^{n+1}+\Delta t\phi_{2}(L_{h}\Delta t)(N(\hat{U}^{n+1}_{i})-N(U^{n}_{i})).\\ \end{split}\right.$$

  (3.12)

In (3.11) and (3.12), $\phi_{j},j=0,1,2$ are defined as

The calculation process of $\phi_{i}(L_{h}\Delta t)$ can be simplified by 2D discrete cosine transform (DCT) efficiently. The time complexity is $\mathcal{O}((M_{1}M_{2})^{2}log(M_{1}M_{2}))$ per time step. For more detailed illustrations, please refer to [11, 14].

We summarize the algorithm described above as follows.

In this part, we characterize the stabilizer $S$ such that the discrete maximum bound principle holds for the proposed ETD scheme (3.11) (or (3.12)). We assume the phase variable $u_{i}$ in the interval $[0,1]$, and project $u_{i}\in[0,1]$ to $\tilde{u}_{i}\in[-1,1]$ by the operator $\mathcal{P}$

Therefore the original image $I$ and the average intensity $\mathbf{C}$ are also rescaled in $[-1,1]$ and denoted by $\tilde{I},\tilde{\mathbf{C}}$, correspondingly. Substituted $\tilde{u}_{i},\tilde{I},\tilde{\mathbf{C}}$ into the Allen–Cahn equation (3.8) and (3.9), the equations become

where $\tilde{w}(u)=\frac{1}{2}u(u-1)(u+1)$. Then the ODE system (3.10) turns to

with

Following the idea in [23], we turn to the proof of discrete maximum bound principle for the ETD numerical schemes of (3.16).

With the estimate in Lemma 3.1, we can deduce the discrete maximum bound principle for ETD1 and ETDRK2 schemes of (3.16) under $\tilde{S}\geq\gamma$, referring to Theorems 3.4 and 3.5 in [7]. That is, the numerical solution $\tilde{U_{i}}^{k}$ satisfies $\left\|\tilde{U}_{i}^{k}\right\|_{\infty}\leq 1$ if the initial value have $\left\|\tilde{U}_{i}^{0}\right\|_{\infty}\leq 1$. $\tilde{U}_{i}$ and $U_{i}$ also satisfy the projection relations in (3.13). Thus, the numerical solution of the ETD scheme (3.11) (or (3.12)) for equations (3.10) preserves the discrete maximum bound principle, as shown in the following theorem.

Theorem 3.2 demonstrates that for fixed $\mathbf{C}$, each $U_{i}^{k}$ can be kept in $[0,1]$, which causes the interval which $\mathbf{C}$ belong to is invariant, after updating $\mathbf{C}$ in Step 2 of Algorithm 2. It guarantees that the discrete maximum bound principle always holds for the whole iteration.

For given $\mathbf{C}$, the Allen–Cahn equations (3.8) and (3.9) are energy stable based on (2.2), viz.,

Subsequently, we will prove that the ETD scheme (3.11) or (3.12) for (3.10) has discrete energy stability in the sense that the following discrete energy is non-increase in time:

The proof for the ETD1 scheme (3.11) can be obtained with a similar argument in [7]. For the ETDRK2 scheme (3.12), the proof in [7] is not optimal, and we can adopt a similar proof in [10] to get the unconditional energy stability in Lemma 3.3. In both proofs, the vital step lies in the boundedness of $\left\|F^{\prime\prime}_{\epsilon_{1}}\right\|$ in (3.24), which can be proved by the boundedness of numerical solutions guaranteed by Theorem 3.2 and estimates on the regularized Heaviside function $H_{\epsilon_{1}}(u)$ and the regularized Dirac function $\delta_{\epsilon_{1}}(u)$. The proofs are quite similar to those in [7] and [10]. So we omit the proof of Lemma 3.3.