 # Two-phase image segmentation by the Allen-Cahn equation and a nonlocal edge detection operator

Based on a nonlocal Laplacian operator, a novel edge detection method of the grayscale image is proposed in this paper. This operator utilizes the information of neighbor pixels for a given pixel to obtain effective and delicate edge detection. The nonlocal edge detection method is used as an initialization for solving the Allen-Cahn equation to achieve two-phase segmentation of the grayscale image. Efficient exponential time differencing (ETD) solvers are employed in the time integration, and finite difference method is adopted in space discretization. The maximum bound principle and energy stability of the proposed numerical schemes are proved. The capability of our segmentation method has been verified in numerical experiments for different types of grayscale images.

## Authors

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

Image segmentation is the process of dividing an image domain into some disjoint regions according to a characterization of the image within or in-between the regions . The characterization can be color, shape, edge, texture, or any feature that can distinguish each region from others. Image segmentation is involved in many fields, such as remote sensing, medical image recognition, robotics, and visual field monitoring.

Among many segmentation methods, variational methods have attracted considerable attention. A typical variational method for image segmentation is based on minimizing an objective energy functional, for instance, the Mumford-Shah (MS) model ,

 EMS(u,I)=∫Ω∖Γ|∇u|2dx+μLength(Γ)+λ∫Ω(u−I)2dx, (1.1)

which is associated with the partition determined by original image and a union of closed edges . The function in the MS model is a piecewise smooth approximation to . and are positive constants. Under many circumstances, the function can degenerate to piecewise constant  and the following model can be achieved:

 ECV(Γ,C1,C2)=Per(Γ;Ω)+λ1∫Γ(C1−I)2dx+λ2∫Ω∖Γ(C2−I)2dx, (1.2)

which is called the two-phase piecewise constant MS model, also known as Chan-Vese model considered previously by Chan and Vese with level set formulation in [8, 29]. Here, denotes the perimeter of the closed curve between two regions and , and are positive constants, and are the average intensities of two phases respectively, which change along with the image intensity. Generally, the objective energy functional contained two parts, including image and partition constraints. Image constraints yield functional terms called fitting terms, which measure how close an approximation fits the original image. Partition constraints give rise to partition priors which usually describe a geometric aspect of the edges, such as their length and smoothness. These two constraints have a wide variety, based on the region [18, 19, 33] or edge [6, 17].

Different approximations combined with different constraints result in various energy functionals. To solve these models for image segmentation, the level set method, phase-field method, and threshold dynamic method have been successfully used. Wherein, level set methods were proposed by Osher and Sethian in [26, 27]. The closed curves describing a segmentation region can split or merge flexibly to implement segmentation with complex structures. Subsequently, it has been found that the energy functional can be replaced by another form and derived the variational level set method [18, 25], which has advantages on numerical stability. Inspired by the variational level set method, various phase-field models have been derived for image segmentation through the -convergence  for the length term of partition constraints, see e.g., [2, 16, 21, 32] and references therein. Merriman, Bence, and Osher(MBO)  developed a threshold dynamic method for the motion of an interface driven by the mean curvature. As pointed in 

, the length term of partition constraints can be estimated by the convolution of a heat kernel. As a consequence, a kind of threshold dynamic methods

[13, 30, 31] had been constructed and obtained effective segmentation results.

In this paper, we plan to adopt a phase-field approach of the Chan-Vese model (1.2) for the two-phase segmentation of grayscale images. An alternating minimization method will be used to update the phase variable , the two-phase average intensities and iteratively. Two exponential time differencing (ETD) schemes will be proposed to solve the resulted Allen-Cahn equation for . The discrete maximum bound principle and energy stability will be proved under certain conditions. Notably, there needs an initial value to solve the Allen-Cahn equation. In general, the initial value is determined by taking an threshold value of the original image [3, 16, 21, 32], or selecting a part of the image [8, 18, 19, 33]. An improper selection of initial values sometimes causes unsatisfactory segmentation results. Therefore, the setup of initial value is also a significant step for image segmentation. In this work, we will develop a nonlocal edge detection method, and use the detected edge as the initial value for solving the Allen-Cahn equation. The nonlocal edge detection method derived by using the nonlocal Laplacian operator, which is comparable and even superior to the existing gradient-based edge detection method [5, 7, 28]. With the nonlocal edge detection method, more detailed information and structure can be retained in the initial value so that a more delicate image segmentation can be achieved.

The rest of this paper is organized as follows. In Section 2, we first review some widely-used edge detection operators based on the first-order and second-order derivatives, respectively. Then, we introduce a novel nonlocal edge detection method. Section 3 is devoted to illustrating an iteration algorithm for the image segmentation by solving an Allen-Cahn type phase-field model, where we use the result of nonlocal edge detection as the initial value to start the iteration process. Plenty of numerical experiments are given in Section 4 to demonstrate the effectiveness and efficiency of our proposed methods. Section 5 concludes the paper.

## 2 Nonlocal edge detection method

In this section, we first review some mature gradient-based edge detection methods, and then introduce the nonlocal Laplacian operator. Based on the nonlocal Laplacian operator, we develop a novel nonlocal edge detection method.

### 2.1 Some gradient-based edge detection operators

When the image intensity of a grayscale image changes intensely, the corresponding gradient will emerge the local extremum, and it also associates with the zero-cross point of Laplacian. Thus the local extremum of gradients and the zero-cross point of Laplacian would be two significant measurements for detecting edge.

Some detection methods are based on the first-order derivative, e.g., Roberts, Prewitt and Sobel operators . They use threshold to estimate the local extremum, after obtaining the corresponding approximate derivatives. The difference of aforementioned three methods lies in different weights of neighbor pixels. Though the consideration of neighbor pixels depresses the irrelevant noise, the edge detected by the first-order derivative edge detection methods is wider than expected, which makes the detected edge less accurate. Besides, the choice of the threshold value is by running trials, which causes the nondeterminacy of the detected edge.

The well-known Canny edge detector  is looking for local extremum of the gradient. After a Gaussian filter on the gradient, this method uses two thresholds to detect edges, including weak edges in the output if they are connected to strong edges. By using two thresholds, the Canny method is less likely than the other methods to be disturbed by noise, and more likely to detect true weak edges.

For edge detection methods based on the second-order derivative, the zero-cross point of the Laplacian operator will indicate the edge. Applying the Laplacian operator to the original image directly will cause the double-edge phenomenon. Another deficiency of the Laplacian operator is the sensitivity to the noise. To fix these problems, Marr proposed the Laplacian of a Gaussian operator as a remedy :

 L=∇2[G(x,y)∗I]=[∇2G(x,y)]∗I,G(x,y)=12πς2exp(−x2+y22ς2),

where . Bigger gives us more blur of the image. After the above convolution, the zero-cross point of indicates the desired edge.

### 2.2 Nonlocal Laplacian operator

The pixel structure of the image provides a natural two-dimensional discretization on a fixed rectangular mesh, therefore we give the elaboration of the 2D nonlocal operator in the subsequent description. For the sake of convenience, the boundary condition is assumed to be the homogeneous Neumann boundary

 ∂I∂x=0

on the boundary layer is defined as

 ΩI={x∉Ω| dist(x,Ω)≤δ}.

Here, is the area where the original image belongs to. Without special illustration, we consider the grayscale image as a 2D matrix whose element denotes the pixel in the corresponding region. Considering that the coordinate of image is always an integer, the interaction radius is also taken as an integer.

#### 2.2.1 General representation

We consider the following nonlocal Laplacian operator [10, 11]

 Lδu(x)=12∫Bδ(0)ρδ(∥s∥)(u(x+s)−2u(x)+u(x−s))ds, x∈Ω,

where is the ball with radius in plane centered at the origin, is a nonnegative kernel function, and stands for the usual Euclidean norm. is further assumed to satisfy

 ∫Bδ(0)∥s∥2ρδ(∥s∥)ds=2d,

where dimension .

Given the uniform square mesh as , . At any node , we rewrite the nonlocal Laplacian operator as

 Lδu(xi)=12∫Bδ(0)u(xi+s)−2u(xi)+u(xi−s)∥s∥2∥s∥1∥s∥2∥s∥1ρδ(∥s∥)ds.

With a quadrature-based finite difference discretization , the nonlocal Laplacian operator becomes

 Lδ,hu(xi)=12∫Bδ(0)Ih(u(xi+s)−2u(xi)+u(xi−s)∥s∥2∥s∥1)∥s∥2∥s∥1ρδ(∥s∥)ds.

The operator

here denotes the piece-wise multi-linear interpolation with respect to

, which can be represented as

 Ihv(s)=∑sjv(sj)ϕj(s).

Here, the basis function satisfies

 ϕj(si)={0,if i≠j,1,if i=j.

Hence, the discrete nonlocal Laplacian operator is formulated as follows:

 Lδ,hu(xi)=∑0≠sj∈Bδ(0)u(xi+sj)−2u(xi)+u(xi−sj)∥∥sj∥∥2∥∥sj∥∥1βδ(sj), xi∈Ω

where the coefficient is

The operator is self-adjoint and negative semi-definite.

#### 2.2.2 Implementation on edge detection

Let be the pixel value of the given image at the mesh point and is a positive integer. The formulation of the nonlocal Laplacian operator is simplified as

 Lδ,hIi,j=r∑p=0r∑q=0cp,q(Ii+p,j+q+Ii−p,j+q+Ii+p,j−q+Ii−p,j−q−4Ii,j)

where and

 cp,q=p+q(p2+q2)h∫∫B+δ(0)ϕp,q(x,y)ρδ(√x2+y2)x2+y2x+ydxdy.

is the bilinear basis function located at the point . The coefficient is symmetric which can be precomputed to reduce the complexity. In the process of image segmentation, the mesh size is taken as 1. We summarize our proposed nonlocal Laplacian edge detection method in Algorithm 1 below.

Algorithm 1 Step 1: Given , we obtain the nonlocal Laplacian

Step 2: Find the desired edge by choosing a suitable threshold value ,

In Algorithm 1, there are two parameters that need to be adjusted. We make the interaction radius , which guarantees that the detection considers the neighbors’ information as much as possible, meanwhile does not cause low contrast by too large . And the threshold value is taken slightly larger than zero to determine the zero-cross points approximately and detect the desired edge. In general, these two parameters vary with specific images, we will discuss it in Section 4.

To illustrate the algorithm, a synthetic image profile is used as an example. We assume the profile consists of a weak edge, a noise point, a jump edge, and a stair edge

 I1=[5 5 4.5 4.5 weak edge0 0 0 7 0 0 0 0 noise point2 5 2 2 jump edge% 0 0 0 0 7 7 7 7 stair edge]/10

correspondingly.

The nonlocal Laplacian operator is evaluated with . In Fig. 2.1, we can see that the nonlocal Laplacian also has zero-cross points around the points where the image intensity changes sharply. is set to be 0.02 for this example. The intersection of and the nonlocal Laplacian is the desired edge . The nonlocal effect avoids the appearance of zero-cross points at the weak edge and inapparent noise. We can decide whether a weak edge shows or not, by adjusting the threshold value according to the change of the image intensity. The threshold method also reduces the complexity of finding the zero-cross points point-wisely. However, it is still sensitive to strong noise points. To get an accurate segmentation of a grayscale image, we will introduce a phase-field model in next section, which takes the detected edge as initial data.

## 3 Image segmentation by the Allen-Cahn equation

The segmentation of a grayscale image from a simple edge detection is inaccurate and the process is very sensitive to noise. In this section, an Allen-Cahn model will be used to segment the given image to two parts with evident intensity change.

### 3.1 A phase-field approach to the Chan-Vese model

A phase-field approach to the Chan-Vese model (1.2) gives the following energy functional:

 Eϵ=∫Ω(ϵ|∇u|2+1ϵW(u)+F(u))dx, (3.1)

where

 F(u)=λ1[(C1−I)2H(u−12)]+λ2[(C2−I)2(1−H(u−12))].

Here, is a periodic potential function , the Heaviside function is given by

 H(u)={0,for u<0,1,for u≥0,

and is the fitting term for the last two terms of (1.2). The diffusion parameter and the weighted coefficients are positive. The phase variable eventually evolves to two phases 0 and 1. We assume the interface, i.e. the edge, lies at the contour of .

### 3.2 The numerical scheme

The energy functional (3.1) depends on the variables and . We adopt an alternating minimization scheme to minimize the energy functionals . Here . For the propose of illustration, we write the iteration process as

 uk+1=argmin Eϵ(Ck,uk), (3.2) Ck+1=argmin Eϵ(uk+1,Ck). (3.3)

We give a regularization of the Heaviside function by

 Hϵ1(u)=⎧⎪ ⎪⎨⎪ ⎪⎩12ϵ1(u+ϵ1πsin(πuϵ1))+12,  for|u|≤ϵ1,1,  for u>ϵ1,0,  for u<−ϵ1.

Here, we take . The derivative of is an approximation of Dirac function given by 

 δϵ1(u)={12ϵ1(1+cos(πuϵ1)),  for|u|≤ϵ1,0,  for |u|>ϵ1.

For a fixed , the solution for the equation (3.3) can be given by

 C1=∫ΩHϵ1(u−12)Idx∫ΩHϵ1(u−12)dx,  C2=∫Ω(1−Hϵ1(u−12))Idx∫Ω(1−Hϵ1(u−12))dx. (3.4)

Therefore our attention should be paid on solving (3.2) by finding the steady state solution of the following Allen-Cahn equation:

 ut=2ϵΔu−1ϵw(u)−[λ1(C1−I)2−λ2(C2−I)2]δϵ1(u−12), (3.5)

where . We consider a homogenous Neumann boundary condition for this Allen-Cahn equation.

We summarize the algorithm described above as

Algorithm 2 Step 0: Give the initial values and , wherein by the nonlocal edge detection method, and . Step 1: Set a suitable . With the initial value , solve the Allen-Cahn equation by the numerical scheme (3.8) (or (3.9)) till the steady state to get an initial segmentation of the original image. Update by (3.4) with and set a smaller . Step 2: With , solve the Allen-Cahn equation by the numerical scheme (3.8) (or (3.9)) till the steady state to get . Step 3: Calculate the new by (3.4) with . Step 4: Set . Repeat Steps 2-3 until the convergence of .

###### Remark 3.1.

As mentioned in Section 2, the nonlocal edge detection is still sensitive to strong noises. To remove the noise effect, we usually set a larger for Step 1 to get an initial segmentation of the original image. This setup can eliminate the irrelevant broken lines or points detected by the nonlocal edge detection method, which underlines the primary part of the given image. We call this procedure as Stage 1 of the Algorithm 2. In Stage 2, we use a smaller in the remaining iterations to get a relatively sharp interface of two phases. As a consequence of Step 4, the boundary of the two parts is specified gradually as the process of iteration goes.

#### 3.2.1 Exponential time differencing methods

In this part, we propose two ETD methods to solve the Allen-Cahn equation (3.5). Discretizing (3.5) in spatial variables by the finite difference method leads to a system of ODEs:

 Ut+LhU=N(U) (3.6)

where

 Lh=−2ϵDh+SId,N(U)=SU−1ϵw(U)−[λ1(C1−I)2−λ2(C2−I)2]δϵ1(U−12)

with discrete Laplacian matrix , which comes from the central finite different difference discretization of

with the homogeneous Neumann boundary condition, and identity matrix

. The positive constant is a stabilizer. Multiply (3.6) an integrating factor , and integrate from to gives

 U(tn+1)=e−LhΔtU(tn)+e−LhΔt∫Δt0eLhsN(U(tn+s))ds. (3.7)

Then we propose two approximations of (3.7) as follows.

• ETD1 : Approximating with , we obtain

 Un+1=ϕ0(LhΔt)Un+Δtϕ1(LhΔt)N(Un). (3.8)
• ETDRK2: Evaluating by ETD1 and approximating with , we obtain

 {^Un+1=ϕ0(LhΔt)Un+Δtϕ1(LhΔt)N(Un),Un+1=ϕ0(LhΔt)^Un+1+Δtϕ2(LhΔt)(N(^Un+1)−N(Un)). (3.9)

are defined as

 ϕ0(a)=e−a,ϕ1(a)=1−e−aa,ϕ2=e−a−1+aa2.

The key process of calculating from the scheme (3.8) or (3.9) is the efficient implementation of the actions of the exponentials . Since comes from the finite difference discretization of with the homogeneous Neumann boundary condition, we could employ the 2D discrete cosine transform (DCT) in the calculation. We use a square domain to illustrate the implementation, where we have pixels. We introduce as the DCT operator, and then for any ,

 Lh=D^LhD−1,

where

 (^LhV)k,l=λk,lVk,l,1≤k,l≤N.

Here

 λk,l=8ϵsin2kπN+8ϵsin2lπN+S,1≤k,l≤N

are the eigenvalues of

. Then we have 

 ϕi(LhΔt)=Dϕi(^LhΔt)D−1, (ϕi(^LhΔt)V)k,l=ϕi(λk,lΔt)Vk,l, i=0,1,2.

The actions of can be implemented by 2D DCT and its inverse transform, respectively. The computational complexity is per time step.

#### 3.2.2 Discrete maximum bound principle

The solutions of the Allen-Cahn equation (3.5) satisfy the maximum bound principle . Next, we will discuss the requirement for the stabilizer such that the solutions of ETD1 (3.8) and ETDRK2 (3.9) preserve the discrete maximum bound principle.

The initial value lies in the interval . Taking the projection , we have

 ~u=P(u)=2(u−12), u=P−1(~u)=12~u+12. (3.10)

Since we rescale the phase variable , the original image and average intensities need to be changed accordingly. Using the same transformation (3.10), we get the new lying in and the new average intensities . Substituting into the Allen-Cahn equation (3.5), we have

 ~ut−2ϵΔ~u+2ϵ~w(~u)+2[λ1(~C1−~I)2−λ2(~C2−~I)2]δϵ1(12~u)=0, (3.11)

where . The ODE system obtained by the finite difference spatial discretization for equation (3.11) becomes

 ~Ut+~Lh~U=~N(~U), (3.12)

with

 ~Lh=−2ϵDh+~SId,  ~N(~U)=~S~U−2ϵ~w(~U)−2[λ1(~C1−~I)2−λ2(~C2−~I)2]δϵ1(12~U).

Subsequently, we prove the solutions of ETD1 (3.8) and ETDRK2 (3.9) for the transformed equation (3.11) preserve the discrete maximum bound principle.

###### Lemma 3.1.

The nonlinear term is bounded by , i.e. , provided that

 ~S≥~L≜2π2ϵ+2λπϵ21,λ≜max(λ1,λ2).
###### Proof.

The derivative of is given by

 ~N′(ξ)=~S−2ϵ~w′(ξ)−[λ1(~C1−~I)2−λ2(~C2−~I)2]δ′ϵ1(12ξ). (3.13)

For the second term of (3.13), is bounded is . As for the third term, the coefficient of fitting term can be estimated by

 −4λ2≤[λ1(~C1−~I)2−λ2(~C2−~I)2]≤4λ1

i.e.

 ∣∣λ1(~C1−~I)2−λ2(~C2−~I)2∣∣≤4max(λ1,λ2)≜4λ,

since represents the rescaling image whose elements are in , and the average intensities are also bounded by . And we know . Hence, the third term of is bounded by . Let

 ~L≜2π2ϵ+2λπϵ21,

and then we have

 ~N′(ξ)≥~S−~L≥0,

when the stabilizer satisfies . On the other hand, the nonlinear term also has property

 ~N(−1)=−~S,~N(1)=~S.

The monotonicity leads us to the conclusion,

 ∣∣~N(ξ)∣∣≤~S,ξ∈[−1,1]

under the condition . ∎

With the bound of the nonlinear term in equation (3.12), the discrete maximum bound principle for the ETD1 (3.8) and ETDRK2 (3.9) can be derived by Theorem 3.4 and Theorem 3.5 in  when the stabilizer . That is, the numerical solution satisfies when the initial value is bounded by and . and also satisfy the relations in (3.10), so numerical solutions of ETD1 (3.8) and ETDRK2 (3.9) for equation (3.6) preserve the discrete maximum bound principle, which is given by the following theorem.

###### Theorem 3.2.

If the initial value satisfies , then numerical solutions of ETD1 (3.8) and ETDRK2 (3.9) for equation (3.6) is also in the interval , for all , provided that the stabilizer satisfies

 S≥L,  L=~L. (3.14)

Theorem 3.2 demonstrates that for fixed , our resolution for preserves the discrete maximum bound principle. It shows that the initial value of the next step in the iteration still lies in . After updating in Step 3 of Algorithm 2, the interval, which belong to, is invariant. It guarantees that the discrete maximum bound principle always holds for the whole iterative process.

#### 3.2.3 Discrete energy stability

The Allen-Cahn equation (3.5) holds the energy stability

 dEϵdt≤0.

In this section, we show that the proposed ETD schemes could preserve this property in the discrete sense. For a rectangular image with pixels, the solution of (3.6), . The discretized energy is defined as defined by

 Eϵh=MN∑i=1(1ϵW(U)+Fϵ1(U))−ϵUTDhU,∀U∈RMN.

Here,

###### Lemma 3.3.

For the fixed constants , the energy stability of the ETD schemes is given by

 ETD1 (???):Eϵ(Un+1)≤Eϵ(Un), ETDRK2 (???):Eϵ(Un+1)≤Eϵ(Un)+C(1+Δt)2

under the condition . The constant is independent of and .

###### Remark 3.2.

Referring to the proofs of Theorem 5.1 and Theorems 5.2 in , we can prove the discrete energy stability for the ETD1 and EDTRK2 schemes in a similar way. The essence of this proof lies in the estimate of . From a similar argument for the bounds of the second term and third term of (3.13), we can have .

###### Theorem 3.4.

The iterative process (3.2)-(3.3) holds unconditional energy stability under the condition ,

 ETD1 :Eϵ(uk+1,Ck+1)≤Eϵ(uk,Ck), (3.15) ETDRK2 :Eϵ(uk+1,Ck+1)≤Eϵ(uk,Ck)+~C, (3.16)

where the constant is independent of and .

###### Proof.

We give the proof for the iterative process (3.2)-(3.3) with ETD1. The process with ETDRK2 is the same.

For the first stage (3.2), Lemma 3.3 derives

 Eϵ(uk+1,Ck)≤Eϵ(uk,Ck).

And then the second stage (3.3) implies

 Eϵ(uk+1,Ck+1)≤Eϵ(uk+1,Ck).

According to the monotonicity of for , we obtain the results of Theorem 3.4. ∎

###### Remark 3.3.

The energy stability in (3.15) for the ETDRK2 scheme cannot be obtained theoretically. However, numerical experiments in Section 4 will show that the energy stability (3.15) holds numerically for both ETD1 and ETDRK2.

## 4 Numerical experiments

We divided this section into two parts. The first one is the edge detection by the nonlocal Laplacian operator. The other one is the image segmentation by the phase-field approach of the Chan-Vese model. The direct edge detection is sensitive to noise and may causes some broken edges, while the segmentation by the Allen-Cahn equation can fix these problems to handle more complex images.

The fractional power kernel for the nonlocal Laplacian operator we choose is as follows,

 ρδ(r)=2(4−α)πδ4−αrαχ(0,δ](r),α∈[0,4).

In this section, the parameter is taken as 1 for the integrality. Actually, there is no big difference when changes. And the parameter in Stage 2 of Algorithm 2 is fixed as . The parameter involved with is given by 0.5 for simplicity. The stabilizer is set according to the bound given in (3.14). All examples were performed on a Windows desktop system using an Intel Core i5 processor, and programmed in Matlab R2019a.

### 4.1 Examples on edge detection

In this subsection, we make the comparison between the nonlocal edge detection method in Section 2.2 and some gradient-based edge detection methods in Section 2.1. To quantified the error, some criteria are given to measure the image segmentation. Referring to [9, 20], the definitions of three kinds of ratios are listed as follows:

 False positive ratio (FPR)=∥S1−S2∥1/∥S1∥1, False negative ratio (FNR)=∥S1−S2∥1/∥S2∥1,

where is the sum of absolute value of each elements for the given matrix. is a synthetic exact image segmentation given by a two-phase image. And is a numerical segmentation result of the corresponding segmentation method. It is clearly seen that the smaller those criteria are, the more accurate the corresponding segmentation we have. The nonlocal interaction radius changes slightly according to different images, and the threshold value is taken as 0.05.

###### Example 4.1.

We begin with a synthetic color image (size [1024,1024]) to test our nonlocal Laplacian operator. Because of the three channels of a color image, we only need one channel to extract the edge. Without losing generality, we just exhibit the result of channel 1 in Table 4.1 and Figure 4.2. It is shown that our nonlocal detection method can achieve smaller segmentation error than Roberts, Sobel and Log operators. Our nonlocal detection method is comparable to the Canny operator which is a pretty delicate gradient-based edge detection method. The nonlocal radius here is taken as 3.

###### Example 4.2.

We choose a realistic image ‘cameraman’(size [256,256]), which is used frequently to test the edge detection operator. In this example, the result segmented by the Canny operator is taken as the reference. The segmentation in Figure 4.3 shows that the right hand of the ‘cameraman’ can be recognized by the Canny operator, where our nonlocal operator has some defects. Whereas, the recognition of the building and the grassland makes our nonlocal edge detection method outperform other detection methods. Here is set to be 4.

### 4.2 Examples on image segmentation by the Allen-Cahn equation

When we solve the Allen-Cahn equation (3.5) to segment the given image, the whole iteration consists of two stages. An initial segmentation is generated in Stage 1 with larger . Stage 2 with smaller specifies the desired edge. We denote as the steps to reach the steady state solution of the Allen-Cahn equation in Stage 1, and for each steady state in Stage 2. The number of total loops for updating is . In the simulation, we record the whole CPU time and verify the discrete maximum bound principle of our methods.

###### Example 4.3.

This example considers the segmentation of a synthetic image (size ) polluted by Gaussian noise with zero mean and standard derivation . The image with such kind of noise can not be segmented by the nonlocal edge detection method directly. The phase-field approach of the Chan-Vese model is applied to segment this image in Figure 4.4. The red line is our detected edge plotted as the contour . The corresponding parameters are listed as , for Stage 1, interaction radius and threshold value for the initial nonlocal detection.

The comparison of two initialization methods, the threshold method and nonlocal detection method, is contained in Table 4.2 for first-order (ETD1) and second-order (ETDRK2) schemes, respectively. Compared with the first-order scheme, the second-order scheme can reduce validly, likewise the number of exterior loop . Because this image is syntectic, the exact segmentation is available. The error between the exact segmentation and the numerical result for all four cases are similar and quite accurate. Although the time taken by the nonlocal-detection initialization method is a little longer than the threshold method, we will explain the superiority of the nonlocal-detection initialization method with some real images later. The maximum bound principle is confirmed in the last two column of Table 4.2. In summary, the nonlocal detection method is comparable to the threshold method to generate initial data for solving the Allen-Cahn equation. The energy stability is shown numerically in Figure 4.5. It can be observed that both schemes are energy non-increasing, which implies that our numerical schemes are energy stable.

###### Example 4.4.

This example aims at segmenting the image ‘tiger in the lake’ (size [294,426]). We use two initialization methods to obtain the initial segmentation . For the threshold method, we made the threshold value . Then in Figure 4.6, four plots are listed for solutions from ETD1 and ETDRK2 with two initialization methods, respectively. To underline the segmentation, we use the original image as the background in all four plots. The second-order scheme can avoid the inaccuracy on the neck of the tiger by the first-order scheme. In addition, the nonlocal-detection initialization method can ignore the irrelevant area in the top-right corner without affecting the primary part. And the shape of grass in bottom-right corner is recognized more clearly. As for the most difficult part ‘the tail’, our segmentation is consecutive without any discontinuity. More comparisons of this image ‘tiger’ can be found in . More specifically, some relevant data are stated in Table 4.3. The ETDRK2 scheme proceeding with the nonlocal edge detection initialization method can result in better segmentation with less CPU time. The corresponding parameters are listed as , for Stage 1, interaction radius and threshold value .

###### Example 4.5.

This example illustrates the sensitivity of threshold value by image ‘vessel’ (size [344,462]). In Figure 4.7, we use three different threshold values to give the initial value of the Allen-Cahn equation at Step 1 of Algorithm 2. We can see that the final segmentation highly depends on the initialization, i.e., the threshold value . When is taken as 0.4, some noise areas are detected, which make the segmentation inaccurate. Conversely, the image is over-smoothing that some small branches of this vessel can not be recognized with . Consequently, we narrow down the range of between 0.4 and 0.5, then the final segmentation with is better than the above two. If we shrink the range further, some detailed information should be detected more accurately. However, a bigger will remove the small structures, while a smaller will show more details but have irrelevant noises. That is, the threshold initial value always exists the less-smoothing or over-smoothing issue. By contrast, the nonlocal edge detection can capture the detailed texture of the given image. By adjusting the parameters and , more information of the image ‘vessel’ can be detected as an initial value, from which we can get a more delicate segmentation. In Figure 4.8, the nonlocal edge detection method is applied to the initialization, which led to a more specific segmentation. In Table 4.4, we give the performance comparison of ETD1 and ETDRK2 schemes with threshold initialization and nonlocal edge detection initialization. It is noticed that the ETDRK2 scheme initialized by the nonlocal detection is more efficient than other methods. More comparisons of this image ‘vessel’ can be also found in . The remaining parameters are , for stage 1 , interaction radius and threshold value .

## 5 Conclusion

In this paper, we first propose a novel nonlocal edge detection method by using the nonlocal Laplacian operator. In practice, the detection results of this method are similar even superior to some widely used gradient-based edge detection methods. Then to get a better two-phase segmentation of grayscale images, we adopt a phase-field approach of the Chan-Vese model. An alternating minimization algorithm is developed to solve this model. For the derived Allen-Cahn equation, it is solved by exponential time integration and finite difference discretization in space. The nonlocal edge detection gives the initial value for solving the Allen-Cahn equation at the first step. For the developed ETD1 and ETDRK2 schemes, the discrete maximum bound principle and energy stability have been proved theoretically and verified numerically. A variety of numerical experiments have been given to demonstrate the effectiveness of our proposed methods. It is notable that the initial value provided by the nonlocal edge detection method can generate better segmentations than those from the traditional threshold initialization method. And we could also observe that the ETDRK2 method with the nonlocal edge detection initialization usually gives better results with less CPU time in the simulation. One future direction of this research is to study the segmentation of images with intensity inhomogeneity. Multi-phase image segmentation and segmentation of color images are also of our interests.

## Acknowledgments

Z. Qiao’s work is partially supported by the Hong Kong Research Council RFS grant RFS2021-5S03 and GRF grants 15300417 and 15302919. Q. Zhang’s research is supported by the 2019 Hong Kong Scholar Program G-YZ2Y.

## References

•  L. Ambrosio and V. M. Tortorelli. Approximation of functional depending on jumps by elliptic functional via t-convergence. Communications on pure and applied mathematics, 1990, 43(8): 999-1036.
• 

A.L. Bertozzi and A. Flenner. Diffuse interface models on graphs for classification of high dimensional data. SIAM Review, 2016, 58(2): 293-328.

•  X. Cai, R. Chan and T. Zeng. A two-stage image segmentation method using a convex variant of the Mumford–Shah model and thresholding. SIAM J. Imaging Sci., 2013, 6(1): 368-C390.
•  J. Canny. A computational approach to edge detection. IEEE Transactions on pattern analysis and machine intelligence, 1986 (6): 679-698.
•  A.K. Cherri and M.A. Karim. Optical symbolic substitution: edge detection using Prewitt, Sobel, and Roberts operators. Applied optics, 1989, 28(21):4644-8.
• 

V. Caselles, R. Kimmel and G. Sapiro. Geodesic active contours, in: Processing of IEEE International conference on computer vision 95, Boston, MA, 1995, pp.694-C699.

•  K.R. Castleman. Digital image processing. Prentice Hall, 1998.
•  T.F. Chan and L.A. Vese. Active contours without edges. IEEE Transactions on image processing, 2001, 10(2):266-277.
•  R. Crandall. Image segmentation using the Chan-Vese algorithm. Project report from ECE, 2009, 532: 1-23.
•  Q. Du, Z. Tao, X.C. Tian and J. Yang. Asymptotically compatible discretization of multidimensional nonlocal diffusion models and approximation of nonlocal Green’s functions. IMA Journal of numerical analysis, 2019, 39: 607-625.
•  Q. Du, L. Ju, X. Li and Z. Qiao. Maximum principle preserving exponential time differencing schemes for the nonlocal Allen-Cahn equation. SIAM Journal on numerical analysis, 2019, 57 (2): 875-898.
•  Q. Du, L. Ju, X. Li and Z. Qiao. Maximum bound principles for a class of semilinear parabolic equations and exponential time differencing schemes. SIAM Review, 2021, In Press.
•  S. Esedog̃lu and Y.H.R. Tsai. Threshold dynamics for the piecewise constant Mumford-Shah functional. Journal of computational physics, 2006, 211(1):367-384.
•  S. Esedog̃lu and F. Otto. Threshold dynamics for networks with arbitrary surface tensions. Communications on pure and applied mathematics, 2015, 68.
•  N.J. Higham. Functions of Matrices: Theory and Computation. SIAM, Philadelphia, 2008, https://doi.org/10.1137/1.9780898717778.
• 

Y. Jung, S. Kang and J. Shen. Multiphase image segmentation via Modica-Mortola phase transition. SIAM applied mathematics, 2007, 67, 1213-1232.

•  M. Kass, A. Witkin and D. Terzopoulos. Snakes: Active contour models. International journal of computer vision, 1988, 1(4): 321-331.
• 

C. Li, C. Xu and C. Gui, et al. Level set evolution without re-initialization: A new variational formulation. IEEE Computer society conference on computer vision and pattern recognition, 2005.

•  C.M. Li, C. Kao, J. Gore and Z. Ding. Implicit active contours driven by local binary fitting energy, in: IEEE Conference on computer vision and pattern recognition, 2007.
• 

B. Liu, H.D. Cheng and J. Huang, et al. Probability density difference-based active contour for ultrasound image segmentation. Pattern recognition, 2010, 43(6): 2028-2042.

•  Y. Li and J. Kim. Multiphase image segmentation using a phase-field model. Computers and mathematics with applications, 2011, 62(2):737-745.
•  D. Marr and E. Hildreth. Theory of edge detection. Proceedings of the royal society of London, 1980, 207(1167):187-217.
•  D. Mumford and J. Shah. Optimal approximations by piecewise smooth functions and associated variational problems. Communications on pure and applied mathematics, 1989, 42(5):577-685.
•  B. Merriman, J.K. Bence and S. Osher. Diffusion generated motion by mean curvature. Department of Mathematics, University of California, Los Angeles, 1992.
•  A. Mitiche. Variational and level set methods in image segmentation. Springer Berlin Heidelberg, 2011.
•  S. Osher and J.A. Sethian. Fronts propagating with curvature-dependent speed: algorithms based on Hamilton-Jacobi formulations. 1988, 79(1):12-49.
•  S. Osher, R. Fedkiw and K. Piechor. Level set methods and dynamic implicit surfaces. Applied mechanics reviews, 2004, 57(3):xiv+273.
•  V. Torre and T.A. Poggio. On edge detection. IEEE Transactions on pattern analysis and machine intelligence, 1986 (2): 147-163.
•  L.A. Vese and T.F. Chan. A multiphase level set framework for image segmentation using the Mumford and Shah model. International Journal of computer vision, 2002, 50(3):271-293.
•  D. Wang and X.P. Wang. The iterative convolution-thresholding method (ICTM) for image segmentation. arXiv preprint arXiv:1904.10917, 2019.
•  D. Wang, H. Li and X. Wei, et al. An efficient iterative thresholding method for image segmentation. Journal of computational physics, 2017, 350: 657-667.
•  W. Yang, Z. Huang and W. Zhu. Image segmentation using the Cahn-Hilliard equation. Journal of scientific computing, 2019, 79(2):1057-1077.
•  K. Zhang, H. Song and L. Zhang. Active contours driven by local image fitting energy. Pattern recognition, 2010, 43(4): 1199-1206.