Fast and Efficient Numerical Finite Difference Method for Multiphase Image Segmentation

We present a simple numerical solution algorithm for a gradient flow for the Modica–Mortola functional and numerically investigate its dynamics.-e proposed numerical algorithm involves both the operator splitting and the explicit Euler methods. A time step formula is derived from the stability analysis, and the goodness of fit of transition width is tested. We perform various numerical experiments to investigate the property of the gradient flow equation, to verify the characteristics of our method in the image segmentation application, and to analyze the effect of parameters. In particular, we propose an initialization process based on target objects. Furthermore, we conduct comparison tests in order to check the performance of our proposed method.


Introduction
We consider a practical and efficient numerical method to solve the following equation: where ϕ(x, t) is a phase-field function in space x and time t and ϵ is a positive constant. Here, Ω is a bounded domain in R d (d � 1, 2). We apply the homogeneous Neumann boundary condition n · ∇ϕ � 0 on zΩ, where n denotes the unit normal vector to zΩ. Equation (1) depicts a gradient flow equation in the L 2 (Ω)-norm for the Modica-Mortola functional [1,2]: (2) Figure 1 shows the multiple periodic well potential, (1/ε)sin 2 (πϕ), given by the second term in equation (2).
Phase-field models have been actively studied until recently in various fields [3,4]. Bogosel et al. [5] proposed an efficient phase-field method based on a multiphase Γ-convergence. Oudet and Santambrogio [6] presented the Modica-Mortola approximation for branched transport. Conti et al. [7] investigated elastic shape optimization. Pegon et al. [8] studied fractal shape optimization in branched transport. Kang et al. [9] proposed a new clustering model using regularized k-means. Furthermore, there are numerous studies addressing image processing based on partial differential equations (PDEs) using the variational method [10][11][12][13][14][15]. Li [16] proposed a new multiphase image segmentation considering the images' morphological diversity. In [17], the authors introduced a new image segmentation method using multiphase segmentation based on binary image segmentation. Hu et al. [18] presented the minimum barrier superpixel segmentation which is capable of achieving cutting-edge performance. Liu et al. [19] proposed a graph-cut-based minimization method for image segmentation. e authors in [20] proposed a parametric snake model for image segmentation. Image segmentation methods based on fractional calculus are popular emerging methods [21,22]. Recently, deep learning-based segmentation methods have been actively studied in the image segmentation [23][24][25][26].
Several research studies have been conducted for multiphase image segmentation. Inspired by the work of Jung et al. [1], the authors in [27] presented a hybrid numerical method for multiphase image segmentation using a phasefield model. ey used a periodic quartic polynomial instead of a sinusoidal function. A simplified Modica-Mortola functional was introduced for estimating piecewise-constant solutions for Fredholm integral equations of the first kind and a damped Newton method was employed to minimize the functional [28]. Huska et al. [29] presented an extension of the Mumford-Shah model for multiphase image segmentation using alternating directions methods of multipliers.
Most of the aforementioned methods use implicit schemes which update a solution by solving a discretized equation involving both the current state (n time level) of a given system and later one (n + 1 time level). However, in this study, we propose an explicit segmentation scheme which directly updates the state of the given system at a later time level from the state of the system at the current time level without applying any iterative methods. e detailed process of the proposed explicit segmentation scheme is described in Section 3.3. e main purpose of this paper is to present a simple explicit numerical solution algorithm to the gradient flow for the Modica-Mortola functional.
is paper is organized as follows. In Section 2, we propose the numerical solution algorithm. In Section 3, we perform several characteristic numerical experiments to demonstrate the efficiency of the proposed numerical method. Finally, conclusions are provided in Section 5.

Numerical Solution Algorithm
We discretize equation (1) in two-dimensional space, i.e., Ω � (a, b) × (c, d). e one-dimensional case is defined in the same manner. Let N x and N y be positive integers, h � (b − a)/N x be the uniform mesh size, and Ω h � (x i , y j ): where Δt is the time step. Let be the discrete l 2 -norm. In [1], the authors developed a convex-splitting algorithm for minimization of the nonconvex energy functional. We propose the following operator splitting method (OSM) [30] for equation (1). Let us rewrite equation (1) as where . en, OSM is as follows: we solve equation (4) with Next, using the solution obtained in the first step ψ(x, Δt), we solve equation (5) Finally, we set ϕ(x, (n + 1)Δt) � φ(x, Δt). Now, we describe the above procedure in discrete approximation. First, we solve equation (4) using the explicit Euler method. where Here, the homogeneous Neumann boundary condition [31] is applied. Second, we solve equation (5) using the explicit Euler method.
For the stability of the scheme, we require a constraint for the time step. From equation (6), we have Δt < h 2 /(8ε). Because the right hand side term in equation (7) is periodic, we only need to check the stability with 0 < ϕ * ij < 0.5. More precisely, the sine function is 2π-periodic while the sign is not important in the stability, and hence we only need to consider a half cycle. For this ϕ * ij , we require 0 < ϕ n+1 ij < ϕ * ij . e inequality ϕ n+1 ij < ϕ * ij is consistently satisfied for any Δt, Figure 1: Multiple periodic well potential (1/ε)sin 2 (πϕ). Here, ϵ � 1 is used.

Numerical Experiments
In this section, we perform various numerical experiments to demonstrate the efficiency of the proposed method and to investigate the property of the gradient flow equation. We denote N as the number of iterations, h as the size of spatial step, and Δt as the time step scale. We terminate the main loop whence the discrete l 2 -norm of the relative difference between two consecutive solutions is less than a given tolerance level, i.e., for some n. Here, all the numerical experiments are conducted on 2.70 GHz CPU and 4 GB DDR4 memory. We employ MATLAB R2019a as our test tool.

Phase Evolution in One-Dimensional Computational
Domain. Unless specified otherwise, we fix h � 0.01 in this section. Let ϕ ∞ � ϕ n+1 be taken as an equilibrium solution.
We compute x α and x β using a linear interpolation. For example, m depicts the index satisfying the condition Figure 2(a). Similarly, x β can be obtained in a similar way. We define the transition thickness of ϕ ∞ as L � x β − x α , which is schematically illustrated in Figure 2 Figure 2(c) shows the numerical steady states with different ε � 2h, 4h, . . . , 20h values. e initial condition in the computational domain Ω � (0, 1) is ϕ(x, 0) � 0 if x < 0.5, and ϕ(x, 0) � 1 otherwise. e tolerance tol � 1 e− 16 is used. Table 1 lists the thickness (L � x β − x α ) of the interface transition layer for various ϵ values. Figure 3 shows the transition thickness L � x β − x α against ϵ. We observe the linear relationship between L and ϵ. Fitting the data with a linear function, we have Note that we employ the linear regression based on the least-squares method in order to find such a linear function. erefore, if we want to have L as the interface thickness, then we can use equation (11) to define the value of ϵ. Figure 4 shows the numerical steady states with initial conditions ϕ(x, 0) � 0 if x < 0.5, and ϕ(x, 0) � m otherwise. Figures 4(a), 4(b), 4(c), and 4(d) are the results with m � 1, 2, 3, and 4, respectively. ε � ε(4h), tol � 1 e-16, and Δt with s � 0.1 are used. Figure 5 shows the numerical steady states with initial conditions ϕ(x, 0) � 0 if x < 0.5, and ϕ(x, 0) � m otherwise. Figures 5(a), 5(b), 5(c), and 5(d) are the results with m � 1, 2, 3, and 4, respectively. ε � ε(8h), tol � 1 e− 16, and Δt with s � 0.1 are used. Compared to the case with ε � ε(4h), the width of the transition layer is larger in the results with ε � ε(8h). Figure 6 shows the temporal evolution of ϕ with an initial condition on a computational domain Ω � (0, 1): where rand (x) is a random number between − 1 and 1. e parameters h � 0.01, ε � ε(4h), and Δt with s � 0.1 are used. In Figure

Phase Evolution in Two-Dimensional Computational
Domain. In this section, we begin our experiment in a twodimensional computational domain defined by Ω � (0, 1) × (0, 1). To illustrate the robustness of the proposed model, Figure 7 shows the temporal evolution of ϕ: Figures 7(a), 7(b), and 7(c) are snapshots at t � 0, Δt, and 7Δt, respectively. e initial condition in Ω is given by where rand (x, y) is a random number between − 1 and 1. ε � ε(2h) and Δt with s � 0.99 are used. Figure 8 shows the temporal evolution of ϕ: Figures 8(a), 8(b), and 8(c) represent snapshots at t � 0, Δt, and 4Δt, respectively. e initial condition is defined by perturbed integers from 0, 1, 3, and 6 in 2 by 2 subdomains of Ω. e maximum amplitude of the random number is 0.5. ε � ε(2h) and Δt with s � 0.99 are used.
In addition, we investigate the same problem in Figure 8 with different s values to examine the effect of the safety factor s. Figure 9 shows the temporal evolutions of ϕ with s � 0.11. Compared to the previous result, a small value of s brings out the small time step; hence, it takes a long computational time. erefore, one might use the value of s as large as possible that is guaranteed to converge in order to gain fast convergence.
For the next step, we perform a test with a large variation of ϕ and a larger value of s. Figure 10 shows the temporal evolution of ϕ: Figures 10(a), 10(b), and 10(c) are snapshots at t � 0, Δt, and 2Δt, respectively. e initial condition is defined by perturbed random integers between 0 and 26 in 5 by 5 subdomains of Ω.
e maximum amplitude of the random number is 0. 5 and Δt with s � 0.99 are used. Even though we use a large value of s � 0.99, which is close to the critical value to the numerical stability, we have a good result. Furthermore, we get Figure 10(c) only after 2 iterations.

Image Segmentation Application.
Image segmentation plays a significant role in image unperceiving, computer vision, medical image processing, and pattern identification [32][33][34][35][36]. We want to partition an image into several regions and make each region have characteristics such as edges, color, and object. e segmentation results are mainly for better analysis and meaningful interpretation of images. In this section, we apply our method to multiphase image segmentation with a datafitting term. We assume that there are K + 1 phase fields for minimization of the following Mumford-Shah energy functional based on the Modica-Mortola functional: e first term ε|∇ϕ| 2 makes a smooth transition between phases and penalizes a large amount of oscillation. e reason is that minimizing the square of the gradient indicates a diffusion phenomenon and restricts sudden gradient changes, that is, the oscillation is suppressed because a large amount of sharp interface disappears. e second term represents a multiple periodic well potential, which derives the phase-field variable ϕ to the closest integer value. Note that the reason for using the periodic function as a nonlinear reaction term unlike other models such as Mumford-Shah, Chan-Vese, and other snake contour-based methods is that the multiphase separation based on the multiple periodic well potential can be performed with single-order parameter. e third term is a fitting term, which is defined as follows: where f 0 is the given image, λ is a positive parameter which indicates degree of fidelity, and C k is the average of f 0 in the k-level, where k � 0, 1, . . . , K, i.e., Note that as the value of λ increases, the influence of the edge-stopping function G K becomes stronger, and the tendency to regress to the boundary of the original image becomes stronger.
Here, sinc(ϕ) � sin(πϕ)/(πϕ), which is plotted by solid line in Figure 11, and the dashed line represents function sinc 2 (ϕ). Note that the square of the sinc(ϕ) plays the same role as the delta function in k-level.
Subsequently, we obtain the following gradient descent flow equation:

Mathematical Problems in Engineering
where we assume that C k is constant.
We present the operator splitting numerical algorithm to solve the proposed model. We split the original equation (17) into the following three equations:

Mathematical Problems in Engineering
First, we solve the above three equations by applying the fully explicit method with a known value ϕ n , that is, where we added a small value δ in the denominator to avoid singularities.
Note that if we replace ϕ n+1 ij by ϕ n+2/3 ij on the right side of equation (21), it becomes an implicit method, which generally needs an iterative method to solve it. Now C n k is defined as follows: In the numerical experiments of image segmentation, we normalize the where f min and f max are the minimum and maximum values of the selected image, respectively. For the initial phase-field value, we use the following formula: where T n for n � 0, 1, . . . , N are target levels of the image. For example, if we have four target values, T 0 � 0, T 1 � 0.3, T 2 � 0.8, and T 3 � 1, then we have the initial phase-field value as illustrated in Figure 12. Note that the target values can be randomly selected; however, selecting them according to the distribution of the gray scale of the image ensures fast convergence. As our proposed method involves first observing the image and determining the target values, it is intuitive to follow the distribution of the gray scale of the image. Manual target selection [37][38][39] is a traditional method and has been used so far.

Synthetic Image 1.
We test image segmentation using the image with 5 phases of 240 × 240 mesh size [40]. We fix tol � 1 e− 7 unless otherwise specified. Figure 13 shows the image segmentation using uniformly distributed 5 target values.
Here, s � 0.99, h � 1/240, ε � ε(2h), λ � 10, and Δt � s min[h 2 /(8ε), ε/(2π 2 )] are used. Additionally, we use 5 target values: T 0 � 0, T 1 � 0.25, T 2 � 0.5, T 3 � 0.75, and T 4 � 1. We evaluate the discrete value of the Mumford-Shah energy functional as follows: where C n k is described earlier in equation (22). As shown in Figure 13(c), the synthetic image with 5 uniformly distributed target values is segmented to 5 phases. Moreover, we further consider a normalized ratio of E(ϕ n )/E(ϕ 0 ) in the segmentation process. Figure 13(d) shows the energy plot E(ϕ n ), and the criterion for convergence is n � 5, which is determined by |(E(ϕ n ) − E(ϕ n− 1 ))/E(ϕ 0 )| < tol e . is threshold tol e is set equal to 1e− 4 for the tests with the synthetic image unless otherwise specified. erefore, numerical steady state is achieved at the intersection of two types of stopping criteria, relative consecutive errors, and normalized decaying energy ratios.
Note that the important point is that the target values are set as hyperparameters in the case of real image or automatically set through the distribution of gray-scale values of image. erefore, uniformly distributed target values used in Figure 13 are nothing special, i.e., these can be nonuniformly distributed.
To illustrate the robustness of the proposed model with respect to noise, we conduct an image segmentation experiment for the 5% random noise from the uniform distribution added to the image at first. In this test, we choose the initial image with levels 0, 0.2, 0.4, 0.7, and 1. e image corrupted by 5% noise and the mesh plot of the image are shown in Figures 14(a) and 14(b), respectively. After 15 iterations, we confirm the converged phase field as shown in Figure 14(c) by using Figure 14(d). Compared to Figure 13(d), there is no significant difference in the degree of energy decaying; however, 10 more iterations are needed because the random noise is added, and therefore the convergence of relative consecutive errors takes more iterations relatively.
Subsequently, we generate a kind of impulsive noise, which is called salt and pepper noise, in the synthetic image. Figure 15 depicts the segmentation results with salt and pepper noise. Note that we change the corresponding values of parameters, s � 0.1 and λ � 500, in this case.
Since relatively small s is used, the number of iterations has to increase to achieve a similar level of segmentation as  shown in Figure 14. e value of s is concerned with the accuracy of solutions in this case, which is covered in more detail in the parameter test section below. e boundaries tend to be relatively smeared rather than the results of randomly perturbed noise image segmentation as shown in Figure 14. Another type of impulsive noise is anisotropic noise which is concentrated in specific area. In this case, we generate anisotropically the 20% random noise value from the standard normal distribution. Moreover, we conduct a similar examination for periodic noise. Figure 16 shows the segmentation results initially corrupted by anisotropic and periodic noises. For the anisotropic noise case, we adjust the target values since the specified areas occupy a small portion of the entire area, which implies that the outliers are scaled together when scaling process is done; hence, one may create a buffer as 0.2 to both sides, 0 and 1. We take 0.2, 0.35, 0.5, 0.65, and 0.8 as the target values. is is the reason why Figure 16(a) has a slightly different colormap from others in noise-corrupted images. We further change the value of fidelity, λ � 100, and the others are the same as those of salt and pepper noise case. On the other hand, we adopt the same target values described in Figure 14 for periodic noise case and just vary the parameters s � 0.1, λ � 10. Note that we omit mesh plots and energy plots since they are duplicated as the results listed in Figure 15.
As can be seen from the above results, the segmentation results using our proposed method are robust. Furthermore, we present the long time numerical energy stability for equations (19)- (21) in Figure 17. Here, we use all the same parameters listed above for Figure 14, but for the number of iterations to 500, the final time is T � 500Δt ≈ 30.895.
According to Figure 17, our method has long-term energy stability, which asserts nonoscillation asymptotic behaviors.
ough a phase-field modeling itself does not include an optimizing algorithm, we examine the value of intersection over union (IoU) to measure the accuracy of segmentation. Figure 18 shows original and predicted bounding boxes with respect to different ϵ values, and Table 2 lists the value of IoU for each case.
In this test, we use the parameter values s � 0.1, λ � 100, and the others are the same as used in Figure 14. Note that ϵ is concerned with transition width and the number of pixels in this image is 240 for each axis. erefore, IoU might be similar for ϵ values exceeding a certain scale in diffuse-interface models due to the lack of the number of pixels. is is the reason why we have similar values for ε(4h) and ε(8h).
Moreover, we verify the dullness of edge of our proposed method. Figure 19 shows the dullness of edge with respect to various ϵ values. According to this result, the boundaries become smooth enough as ε increases; it may result in a drawback such as a smeared boundary problem. e theoretical details of both diffuse-interface and sharp-interface limits can be checked in the following references [41,42] and can be also found in the references therein.

Complex Image.
For the following step, we use a more complex image. We fix tol � 1 e− 6 and tol e � 1 e− 4 unless otherwise specified. We take a brain MRI image used in [43] as the first example. In medical science, there are cases in which diagnostic results are obtained by tracking and observing areas that have the same color in the image. erefore, the brain MRI image is suitable for checking the performance of our method. We convert the brain MRI image to gray scale, which ranges from 0 to 1, and choose the pixels as target values T as our region of interest. Fast convergence of the explicit scheme is therefore generally ensured, and one may not continuously select the pixels that have almost the same gray-scale distribution values in reality. If we want to segment regions of intensity at (i k , j k ) for k � 0, 1, . . . , K, then we define C k as the local mean of the near 8 cells and itself.
Note that we fix these defined values until completion of the segmentation process. Here, we use Note that the expression in black and white colors is just to clarify the visualization. erefore, the gray texture in the image can be obtained in exactly the same way. Similarly, the numerical energy stability for the proposed method with   local mean C k is presented in Figure 21. Here, we use all the same parameters listed above for Figure 20, except for the number of iterations (500).

Parameter Test.
We perform parameter tests on image with low contrast and high-level noise in this section. We employ the 132 × 641 sized vertebral image listed in [44], which has low contrast near the boundaries of the object. We initially transform the image into gray scale and perturb with random noise about ±30% within [− 1, 1]. Figure 22 shows both original and noise perturbed vertebral MRI images. We have targeted the spines within the image to segment, and K � 3 is adopted for convenience. Recall that K is a hyperparameter, and hence there is an adequate value of K to achieve fast computation as well as the target values. We fix the default set of parameters as s � 0.6, h � 0.1, λ � 200, and ε � ε(12h). Figure 23 depicts the effect of parameters. We set the tolerance levels for relative consecutive errors to tol � 1 e− 5 and for normalized energy ratios to tol e � 1 e− 3, respectively.
As depicted in Figures 23(a)-23(c), the details near the boundaries of object are captured as the value of s is small. Since s affects the scale of time step Δt, this numerical result of the dynamics of phase field with small value of s fits well with theoretical predictions. However, there is a drawback of segmenting the wrong areas around the boundaries of object if the noise is severe. e effect of space step h is shown in Figures 23(d)-23(f ). It affects the time step size like s; however, it simultaneously affects the value of ϵ. Since we analyzed theoretically the derivation of ϵ for h � 0.01, relatively large value of h leads to a torn segmented image. Figures 23(g)-23(i) illustrate the effect of λ. If λ is small, then the segmented image has a smeared boundary, while the image has a sharp interface with large value of λ. It implies that the fidelity force is dominant if the value of λ is large; therefore, unintended boundaries can be generated inside of object. Note that boundaries in phase-field models are represented by an average level of transition region. erefore, these have sharp-interface limits asymptotically

Comparison Test.
We compare the segmented images followed by using different schemes with respect to the results followed by using our method in this section. Figure 24 shows the segmented images using our method and the methods listed in [27,45]. For more detailed information on the methods, refer to [27,45] and the references therein.
We adopt the parameters ϵ � 0.01, Δt � 0.25h 2 , and μ � 100 for both quadratic regularization and superlinear models. As depicted in Figure 24, the segmented results using our proposed method do not differ significantly from the results using existing methods.
Furthermore, we compare the performance of our simple method to the implicit method used in [27] in terms of the CPU time. We measured the elapsed time until both the stopping criteria are satisfied. Two different time step sizes are used: Δt � 0.5sε/π 2 ≈ 0.0187 with a safety factor s � 0.3 and Δt ref � 1e− 5. We fix ε � ε(2h) for both the schemes. We use the following initial value ϕ 0 , ϕ 0 ij � Kf 0,ij , where K � 3, 4 and target values as T 0 � 0.01, T 1 � 0.35, T 2 � 0.66, T 3 � 0.92 in the case of K � 3 and T 0 � 0.004, T 1 � 0.23, T 2 � 0.52, T 3 � 0.77, T 4 � 0.94 in the case of K � 4. Here, K is the number as described in equation (15). Furthermore, we resize the size of image used in Figure 20 to 256 × 256 since the reference method employs the multigrid method. Table 3 shows the result of the CPU time performance test. If we use Δt ref , then we have similar CPU times for both the methods. However, if we use a larger time step, then the CPU time is reduced by one order of magnitude. is result implies that we can achieve fast computation speed via the proposed method with the moderate time step proposed in Section 2.

Geometric Active Segmentation
For approximation of the geometric active segmentation for multiphase image, the governing equation is given as follows: , t)), where f 0 (x) for x � (x, y) is a given image, g is an edgestopping function, ϕ(x, t) is an order parameter, ϵ is a positive constant, and λ is a constant parameter. Ω, in particular, is a domain bounded in R d (d � 1, 2). Note that we restrict ϕ(x, t) ≥ 0 for convenience. In equation (26), g(f 0 (x)) is the edge-stopping function and is defined as where (G σ * f 0 )(x) is the convolution of the given image f 0 and the Gaussian function Here, we will use p � 2. e function g(f 0 (x)) is close to 1 in homogeneous regions and is close to 0, when the gradient of the image is large. We present the operator splitting numerical algorithm to solve the proposed model. We split the original equation (26) into the following three equations: where Here, (G σ * f 0 ) ij can be computed using a 3 × 3 smoothing kernel as As mentioned above, the homogeneous Neumann boundary condition is applied to the domain.

Synthetic Image 2.
We apply our proposed method to a synthetic image. e synthetic image consists of three rectangles and one circle. at is, we try to separate into four phases. e initial contour is set to the median of two consecutive numbers of 0, 1, 2, 3, and 4. e size of the image is 256 × 256. e parameters are chosen as ε � ε(3h), Δt � 5 e− 6, T � 10000Δt, and λ � 1 e2. e results of the synthetic image segmentation can be seen in Figure 25. Figure 25(a) shows the initial contour-shaped rectangular. While Figure 25

Segmentation with Real-World Images.
In this section, we applied the proposed method on the real-world images, which are chosen from the CMU-Cornell iCoseg Dataset [46], to demonstrate the robustness of the segmentation. Figure 26 shows the comparison results. Figure 26(a) shows the real-world images of CMU-Cornell iCoseg Dataset. Figure 26(b) shows the results of the proposed method. Figure 26(c) shows the ground truth of the color images. As can  be seen from the results, the proposed method can preserve the original shape of the object and shows good performance in image segmentation. Comparing ground truth of the images, we can clearly see that details have been preserved by the proposed method.

Measurement of Quality of the Segmentation.
As quantitative measures of the ability to restore damaged fingerprints to the proposed algorithm, we use peak signal-tonoise ratio (PSNR) and structural similarity index map (SSIM), which are defined as follows:  PSNR � 10log 10 where I ij � f ij /(f max − f min ), K ij � 0.5(ϕ ij + 1); μ I , μ K denote the average of I and K, respectively; and σ 2 I , σ 2 K are the variances of I and K, respectively. Here, σ IK is the covariance of I and K, c 1 � (0.01L) 2 , and c 2 � (0.03L) 2 with L � 1. High PSNR and SSIM values indicate a good segmentation of the images. To verify the efficiency of the proposed method, we compute the average PSNR of the proposed method with the CMU-Cornell iCoseg Dataset, which is shown in Figure 27. From the result, we can see that the curve is increasing, which indicates that the proposed method works well for the image segmentation.

X-Ray Example.
We perform the multiphase image segmentation using a knee X-ray as the original image. e size of given image is 128 × 128. e original image and scalar fields of function g and g 2 are given in Figure 28. Figure 29 shows the result of multiphase image segmentation in two ways: contour and mesh plot. Here, we use ε � ε(6h) for the width of interfacial region, Δt � 1 e− 7, T � 1400000Δt, and λ � 1 e2. Initial contours are set to medians between each assigned discrete value of 0, 1, 2, and 3 so that we can get each initialized value (hence only one) between each two successive integers, respectively. e results yield that the edge-stopping functions g and g 2 work well and one can easily extract each activated contour only using multiple-well potential.

Conclusion
In this paper, we present a simple numerical solution algorithm for the gradient flow for the Modica-Mortola functional and numerically investigate its dynamics. We figure out the property of the corresponding gradient flow equation and analyze the effect of parameters. Especially, we analyze the stability of scheme for time step size to our method and examine the goodness of fit of linear relationship for various ϵ values.
e results indicate that proper phase separations are achieved via our simple explicit method. e proposed numerical algorithm can be applied to multiphase image segmentation problems. We performed various numerical experiments with synthetic and actual MRI images to show the characteristics of our proposed method. We compared the proposed method to several other existing methods and confirmed that the segmented results are not quite different from those of existing methods. In particular, the implicit method using the multigrid method has a mesh size restriction in common, which is of a multiple of 2; however, our method has the advantage of being free to the mesh size. In addition, our proposed method produces good results in less time even though it has a limitation of time step size. Our method can segment objects of simple color images as shown in Section 4.2. However, it is difficult to achieve good results when segmenting complex color images. erefore, for future work, we plan to extend the current model to applications, such as complex nature image segmentation [47,48] and video co-segmentation [49].

Data Availability
e data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
e authors declare that there are no conflicts of interest regarding the publication of this paper.