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 inbetween the regions [25]. 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 MumfordShah (MS) model [23],
(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 [8] and the following model can be achieved:
(1.2) 
which is called the twophase piecewise constant MS model, also known as ChanVese 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, phasefield 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 phasefield models have been derived for image segmentation through the convergence [1] for the length term of partition constraints, see e.g., [2, 16, 21, 32] and references therein. Merriman, Bence, and Osher(MBO) [24] developed a threshold dynamic method for the motion of an interface driven by the mean curvature. As pointed in [14]
, 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 phasefield approach of the ChanVese model (1.2) for the twophase segmentation of grayscale images. An alternating minimization method will be used to update the phase variable , the twophase average intensities and iteratively. Two exponential time differencing (ETD) schemes will be proposed to solve the resulted AllenCahn 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 AllenCahn 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 AllenCahn equation. The nonlocal edge detection method derived by using the nonlocal Laplacian operator, which is comparable and even superior to the existing gradientbased 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 widelyused edge detection operators based on the firstorder and secondorder 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 AllenCahn type phasefield 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 gradientbased 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 gradientbased 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 zerocross point of Laplacian. Thus the local extremum of gradients and the zerocross point of Laplacian would be two significant measurements for detecting edge.
Some detection methods are based on the firstorder derivative, e.g., Roberts, Prewitt and Sobel operators [5]. 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 firstorder 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 wellknown Canny edge detector [4] 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 secondorder derivative, the zerocross point of the Laplacian operator will indicate the edge. Applying the Laplacian operator to the original image directly will cause the doubleedge 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 [22]:
where . Bigger gives us more blur of the image. After the above convolution, the zerocross point of indicates the desired edge.
2.2 Nonlocal Laplacian operator
The pixel structure of the image provides a natural twodimensional 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
on the boundary layer is defined as
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]
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
where dimension .
Given the uniform square mesh as , . At any node , we rewrite the nonlocal Laplacian operator as
With a quadraturebased finite difference discretization [10], the nonlocal Laplacian operator becomes
The operator
here denotes the piecewise multilinear interpolation with respect to
, which can be represented asHere, the basis function satisfies
Hence, the discrete nonlocal Laplacian operator is formulated as follows:
where the coefficient is
The operator is selfadjoint and negative semidefinite.
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
where and
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
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 zerocross 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
correspondingly.
The nonlocal Laplacian operator is evaluated with . In Fig. 2.1, we can see that the nonlocal Laplacian also has zerocross 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 zerocross 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 zerocross points pointwisely. However, it is still sensitive to strong noise points. To get an accurate segmentation of a grayscale image, we will introduce a phasefield model in next section, which takes the detected edge as initial data.
3 Image segmentation by the AllenCahn 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 AllenCahn model will be used to segment the given image to two parts with evident intensity change.
3.1 A phasefield approach to the ChanVese model
A phasefield approach to the ChanVese model (1.2) gives the following energy functional:
(3.1) 
where
Here, is a periodic potential function [16], the Heaviside function is given by
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
(3.2)  
(3.3) 
We give a regularization of the Heaviside function by
Here, we take . The derivative of is an approximation of Dirac function given by [18]
For a fixed , the solution for the equation (3.3) can be given by
(3.4) 
Therefore our attention should be paid on solving (3.2) by finding the steady state solution of the following AllenCahn equation:
(3.5) 
where . We consider a homogenous Neumann boundary condition for this AllenCahn 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 AllenCahn 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 AllenCahn 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 23 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 AllenCahn equation (3.5). Discretizing (3.5) in spatial variables by the finite difference method leads to a system of ODEs:
(3.6) 
where
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(3.7) 
Then we propose two approximations of (3.7) as follows.

ETD1 : Approximating with , we obtain
(3.8) 
ETDRK2: Evaluating by ETD1 and approximating with , we obtain
(3.9)
are defined as
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 ,
where
Here
are the eigenvalues of
. Then we have [15]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 AllenCahn equation (3.5) satisfy the maximum bound principle [12]. 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
(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 AllenCahn equation (3.5), we have
(3.11) 
where . The ODE system obtained by the finite difference spatial discretization for equation (3.11) becomes
(3.12) 
with
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
Proof.
The derivative of is given by
(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
i.e.
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
and then we have
when the stabilizer satisfies . On the other hand, the nonlinear term also has property
The monotonicity leads us to the conclusion,
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 [11] 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.
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 AllenCahn equation (3.5) holds the energy stability
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
Here,
Lemma 3.3.
For the fixed constants , the energy stability of the ETD schemes is given by
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 [11], 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.
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 phasefield approach of the ChanVese model. The direct edge detection is sensitive to noise and may causes some broken edges, while the segmentation by the AllenCahn equation can fix these problems to handle more complex images.
The fractional power kernel for the nonlocal Laplacian operator we choose is as follows,
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 gradientbased 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:
where is the sum of absolute value of each elements for the given matrix. is a synthetic exact image segmentation given by a twophase 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 gradientbased edge detection method. The nonlocal radius here is taken as 3.
Roberts  Sobel  Log  Canny  Nonlocal detection  

FPR  0.6703  0.7430  0.7989  0.6449  0.6952 
FNR  2.0328  2.8917  2.5709  1.8159  1.6759 
RSE  0.5041  0.5911  0.6095  0.4759  0.4914 
CPU time  5.67E02  4.66E02  1.53E01  1.04E01  1.81E01 
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 AllenCahn equation
When we solve the AllenCahn 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 AllenCahn 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 phasefield approach of the ChanVese 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 firstorder (ETD1) and secondorder (ETDRK2) schemes, respectively. Compared with the firstorder scheme, the secondorder 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 nonlocaldetection initialization method is a little longer than the threshold method, we will explain the superiority of the nonlocaldetection 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 AllenCahn equation. The energy stability is shown numerically in Figure 4.5. It can be observed that both schemes are energy nonincreasing, which implies that our numerical schemes are energy stable.
Stage 1  Stage 2  Exteriorloop  CPUtime  Error  Maximum bound  

Scheme  min  1max  
1stthreshold  23  12  1  6  23.1  5.33E04  1.39E12  1.39E12 
2ndthreshold  15  8  1  3  23.8  5.29E04  2.85E11  2.84E11 
1stnonlocal  36  12  1  6  33.3  5.67E04  8.68E15  1.52E14 
2ndnonlocal  23  8  1  3  31.1  5.38E04  1.96E11  3.10E11 
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 secondorder scheme can avoid the inaccuracy on the neck of the tiger by the firstorder scheme. In addition, the nonlocaldetection initialization method can ignore the irrelevant area in the topright corner without affecting the primary part. And the shape of grass in bottomright 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 [32]. 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 .
Stage 1  Stage 2  Exteriorloop  CPUtime  Maximum bound  

Scheme  min  1max  
1stthreshold  83  16  1  4  8.5  6.07E14  5.13E14 
2ndthreshold  65  11  1  5  10.3  3.34E15  4.77E15 
1stnonlocal  83  17  1  4  8.4  1.53E14  2.90E14 
2ndnonlocal  64  11  1  3  9.7  1.75E13  3.48E13 
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 AllenCahn 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 oversmoothing 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 lesssmoothing or oversmoothing 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 [32]. The remaining parameters are , for stage 1 , interaction radius and threshold value .
Stage 1  Stage 2  Exteriorloop  CPUtime  Maximum bound  

Scheme  min  1max  
1stthreshold  127  22  1  13  24.0  7.81E14  2.95E14 
2ndthreshold  93  13  1  9  26.0  1.60E14  6.67E15 
1stnonlocal  145  20  1  11  26.7  1.55E12  6.09E13 
2ndnonlocal  91  13  1  4  23.3  1.28E11  4.70E12 
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 gradientbased edge detection methods. Then to get a better twophase segmentation of grayscale images, we adopt a phasefield approach of the ChanVese model. An alternating minimization algorithm is developed to solve this model. For the derived AllenCahn 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 AllenCahn 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. Multiphase 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 RFS20215S03 and GRF grants 15300417 and 15302919. Q. Zhang’s research is supported by the 2019 Hong Kong Scholar Program GYZ2Y.
References
 [1] L. Ambrosio and V. M. Tortorelli. Approximation of functional depending on jumps by elliptic functional via tconvergence. Communications on pure and applied mathematics, 1990, 43(8): 9991036.

[2]
A.L. Bertozzi and A. Flenner. Diffuse interface models on graphs for classification of high dimensional data. SIAM Review, 2016, 58(2): 293328.
 [3] X. Cai, R. Chan and T. Zeng. A twostage image segmentation method using a convex variant of the Mumford–Shah model and thresholding. SIAM J. Imaging Sci., 2013, 6(1): 368C390.
 [4] J. Canny. A computational approach to edge detection. IEEE Transactions on pattern analysis and machine intelligence, 1986 (6): 679698.
 [5] A.K. Cherri and M.A. Karim. Optical symbolic substitution: edge detection using Prewitt, Sobel, and Roberts operators. Applied optics, 1989, 28(21):46448.

[6]
V. Caselles, R. Kimmel and G. Sapiro. Geodesic active contours, in: Processing of IEEE International conference on computer vision 95, Boston, MA, 1995, pp.694C699.
 [7] K.R. Castleman. Digital image processing. Prentice Hall, 1998.
 [8] T.F. Chan and L.A. Vese. Active contours without edges. IEEE Transactions on image processing, 2001, 10(2):266277.
 [9] R. Crandall. Image segmentation using the ChanVese algorithm. Project report from ECE, 2009, 532: 123.
 [10] 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: 607625.
 [11] Q. Du, L. Ju, X. Li and Z. Qiao. Maximum principle preserving exponential time differencing schemes for the nonlocal AllenCahn equation. SIAM Journal on numerical analysis, 2019, 57 (2): 875898.
 [12] 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.
 [13] S. Esedog̃lu and Y.H.R. Tsai. Threshold dynamics for the piecewise constant MumfordShah functional. Journal of computational physics, 2006, 211(1):367384.
 [14] S. Esedog̃lu and F. Otto. Threshold dynamics for networks with arbitrary surface tensions. Communications on pure and applied mathematics, 2015, 68.
 [15] N.J. Higham. Functions of Matrices: Theory and Computation. SIAM, Philadelphia, 2008, https://doi.org/10.1137/1.9780898717778.

[16]
Y. Jung, S. Kang and J. Shen. Multiphase image segmentation via ModicaMortola phase transition. SIAM applied mathematics, 2007, 67, 12131232.
 [17] M. Kass, A. Witkin and D. Terzopoulos. Snakes: Active contour models. International journal of computer vision, 1988, 1(4): 321331.

[18]
C. Li, C. Xu and C. Gui, et al. Level set evolution without reinitialization: A new variational formulation. IEEE Computer society conference on computer vision and pattern recognition, 2005.
 [19] 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.

[20]
B. Liu, H.D. Cheng and J. Huang, et al. Probability density differencebased active contour for ultrasound image segmentation. Pattern recognition, 2010, 43(6): 20282042.
 [21] Y. Li and J. Kim. Multiphase image segmentation using a phasefield model. Computers and mathematics with applications, 2011, 62(2):737745.
 [22] D. Marr and E. Hildreth. Theory of edge detection. Proceedings of the royal society of London, 1980, 207(1167):187217.
 [23] D. Mumford and J. Shah. Optimal approximations by piecewise smooth functions and associated variational problems. Communications on pure and applied mathematics, 1989, 42(5):577685.
 [24] B. Merriman, J.K. Bence and S. Osher. Diffusion generated motion by mean curvature. Department of Mathematics, University of California, Los Angeles, 1992.
 [25] A. Mitiche. Variational and level set methods in image segmentation. Springer Berlin Heidelberg, 2011.
 [26] S. Osher and J.A. Sethian. Fronts propagating with curvaturedependent speed: algorithms based on HamiltonJacobi formulations. 1988, 79(1):1249.
 [27] S. Osher, R. Fedkiw and K. Piechor. Level set methods and dynamic implicit surfaces. Applied mechanics reviews, 2004, 57(3):xiv+273.
 [28] V. Torre and T.A. Poggio. On edge detection. IEEE Transactions on pattern analysis and machine intelligence, 1986 (2): 147163.
 [29] 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):271293.
 [30] D. Wang and X.P. Wang. The iterative convolutionthresholding method (ICTM) for image segmentation. arXiv preprint arXiv:1904.10917, 2019.
 [31] D. Wang, H. Li and X. Wei, et al. An efficient iterative thresholding method for image segmentation. Journal of computational physics, 2017, 350: 657667.
 [32] W. Yang, Z. Huang and W. Zhu. Image segmentation using the CahnHilliard equation. Journal of scientific computing, 2019, 79(2):10571077.
 [33] K. Zhang, H. Song and L. Zhang. Active contours driven by local image fitting energy. Pattern recognition, 2010, 43(4): 11991206.