# Fast and Efficient Numerical Finite Difference Method for Multiphase Image Segmentation

Authors: Yibao Li; Sungha Yoon; Jian Wang; Jintae Park; Sangkwon Kim; Chaeyoung Lee; Hyundong Kim; Junseok Kim
Journal: Mathematical Problems in Engineering
 (2021)

---

## Abstract

We present a simple numerical solution algorithm for a gradient flow for the Modica–Mortola functional and numerically investigate its dynamics. The 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.


---

## Body

## 1. Introduction

We consider a practical and efficient numerical method to solve the following equation:(1)∂ϕx,t∂t=2εΔϕx,t−πεsin2πϕx,t,x∈Ω,t>0,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 ℝdd=1,2. We apply the homogeneous Neumann boundary condition n⋅∇ϕ=0on∂Ω, where n denotes the unit normal vector to ∂Ω. Equation (1) depicts a gradient flow equation in the L2Ω-norm for the Modica–Mortola functional [1, 2]:(2)Eϕ≔∫Ωε∇ϕ2+1εsin2πϕdx.Figure1 shows the multiple periodic well potential, 1/εsin2πϕ, given by the second term in equation (2).Figure 1
Multiple periodic well potential1/εsin2πϕ. Here, ϵ=1 is used.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–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. The 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–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 phase-field model. They 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. The detailed process of the proposed explicit segmentation scheme is described in Section 3.3. The main purpose of this paper is to present a simple explicit numerical solution algorithm to the gradient flow for the Modica–Mortola functional.This paper is organized as follows. In Section2, 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.


## 2. Numerical Solution Algorithm

We discretize equation (1) in two-dimensional space, i.e., Ω=a,b×c,d. The one-dimensional case is defined in the same manner. Let Nx and Ny be positive integers, h=b−a/Nx be the uniform mesh size, and Ωh=xi,yj:xi=a+i−0.5h,yj=c+j−0.5h,1≤i≤Nx,1≤j≤Ny be the discrete domain. Let ϕijn be approximations of ϕxi,yj,nΔt, where Δt is the time step. Let ϕ2=∑i=1Nx∑j=1Nyϕij2/NxNy be the discrete l2-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(3)∂ϕx,t∂t=f1ϕx,t+f2ϕx,t,where f1ϕx,t=2ϵΔϕx,t and f2ϕx,t=−π/εsin2πϕx,t. Then, OSM is as follows: we solve equation (4) with ψx,0=ϕx,nΔt for a given value ϕx,nΔt to get ψx,Δt.(4)∂ψx,t∂t=f1ψx,t.Next, using the solution obtained in the first stepψx,Δt, we solve equation (5) with φx,0=ψx,Δt to get φx,Δt.(5)∂φx,t∂t=f2φx,t.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.(6)ϕij∗−ϕijnΔt=2εΔhϕijn,1≤i≤Nx,1≤j≤Ny,where Δhϕijn=ϕi−1,jn+ϕi+1,jn−4ϕijn+ϕi,j−1n+ϕi,j+1n/h2. Here, the homogeneous Neumann boundary condition [31] is applied. Second, we solve equation (5) using the explicit Euler method.(7)ϕijn+1−ϕij∗Δt=−πεsin2πϕij∗,1≤i≤Nx,1≤j≤Ny.For the stability of the scheme, we require a constraint for the time step. From equation (6), we have Δt<h2/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<ϕijn+1<ϕij∗. The inequality ϕijn+1<ϕij∗ is consistently satisfied for any Δt, as sin2πϕij∗>0 for 0<ϕij∗<0.5. The inequality 0<ϕijn+1 implies(8)Δt<εϕij∗πsin2πϕij∗,which is a monotonically increasing function with respect to ϕij∗ in 0,0.5. If we take limit as ϕij∗⟶0 to the right hand side term in equation (8), then we obtain(9)limϕij∗⟶0εϕij∗πsin2πϕij∗=ε2π2.Therefore, we can assume the time step to beΔt=sminh2/8ε,ε/2π2, where 0<s<1 is a safety factor. For the one-dimensional case, we have Δt=sminh2/4ε,ε/2π2.


## 3. 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 denoteN 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 l2-norm of the relative difference between two consecutive solutions is less than a given tolerance level, i.e.,(10)1NxNy∑i,jϕijn−ϕijn+1ϕijn<tol,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.


### 3.1. Phase Evolution in One-Dimensional Computational Domain

Unless specified otherwise, we fixh=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 ϕm∞≤0.05<ϕm+1∞, such that xα=xm+h0.05−ϕm/ϕm+1−ϕm is shown in 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(b). Here, ϕ∞xα≈0.05, and ϕ∞xβ≈0.95. Figure 2(c) shows the numerical steady states with different ε=2h,4h,…,20h values. The initial condition in the computational domain Ω=0,1 is ϕx,0=0 if x<0.5, and ϕx,0=1 otherwise. The tolerance tol=1 e−16 is used.Figure 2
(a) Schematic illustration of linear interpolation evaluatingxα. (b) Schematic illustration of transition thickness, L=xβ−xα. (c) Equilibrium solutions with different ε=2h,4h,…,20h values.
(a)(b)(c)Table1 lists the thickness L=xβ−xα of the interface transition layer for various ϵ values.Table 1
ThicknessL=xβ−xα of interface transition layer for various ϵ values.
Caseϵ=2hϵ=4hϵ=6hϵ=8hϵ=10hL2.9903h6.5364h9.7668h12.9192h16.2230hCaseϵ=12hϵ=14hϵ=16hϵ=18hϵ=20hL19.4323h22.6640h25.9192h29.1198h32.3832hFigure3 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(11)εL=0.6158L+0.0004.Figure 3
Plot of interface transition thicknessL=xβ−xα with respect to ϵ.Note that we employ the linear regression based on the least-squares method in order to find such a linear function. Therefore, if we want to haveL as the interface thickness, then we can use equation (11) to define the value of ϵ.Figure4 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 4
Numerical steady states with initial conditionsϕx,0=0 when x<0.5, and ϕx,0=m otherwise. (a), (b), (c), and (d) are the results with m=1,2,3, and 4, respectively. ε=ε4h is used.
(a)(b)(c)(d)Figure5 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 5
Numerical steady states with initial conditionsϕx,0=0 when x<0.5, and ϕx,0=m otherwise. (a), (b), (c), and (d) are the results with m=1,2,3, and 4, respectively. ε=ε8h is used.
(a)(b)(c)(d)Figure6 shows the temporal evolution of ϕ with an initial condition on a computational domain Ω=0,1:(12)ϕx,0=2randx+5x−0.5,where rand x is a random number between −1 and 1. The parameters h=0.01, ε=ε4h, and Δt with s=0.1 are used. In Figure 6, the dotted, solid, and thick solid lines are the initial profile, intermediate, and final solution after 100 iterations, respectively. While initial profile and intermediate solutions have random values, we can observe some flat regions with integer values at the final solution.Figure 6
Temporal evolution of numerical solutions with initial condition,ϕx,0=2randx+5x−0.5.


### 3.2. Phase Evolution in Two-Dimensional Computational Domain

In this section, we begin our experiment in a two-dimensional 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. The initial condition in Ω is given by(13)ϕx,y,0=0.5randx,y+5x+y−1,where rand x,y is a random number between −1 and 1. ε=ε2h and Δt with s=0.99 are used.Figure 7
Temporal evolution of numerical solution: (a), (b), and (c) represent snapshots att=0,Δt, and 7Δt, respectively.
(a)(b)(c)Figure8 shows the temporal evolution of ϕ: Figures 8(a), 8(b), and 8(c) represent snapshots at t=0,Δt, and 4Δt, respectively. The initial condition is defined by perturbed integers from 0, 1, 3, and 6 in 2 by 2 subdomains of Ω. The maximum amplitude of the random number is 0.5. ε=ε2h and Δt with s=0.99 are used.Figure 8
Temporal evolution ofϕ with s=0.99: (a), (b), and (c) are snapshots at t=0,Δt, and 4Δt, respectively.
(a)(b)(c)In addition, we investigate the same problem in Figure8 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. Therefore, one might use the value of s as large as possible that is guaranteed to converge in order to gain fast convergence.Figure 9
Temporal evolution ofϕ with s=0.11: (a), (b), and (c) are the snapshots at t=0,9Δt, and 36Δt, respectively.
(a)(b)(c)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. The initial condition is defined by perturbed random integers between 0 and 26 in 5 by 5 subdomains of Ω. The maximum amplitude of the random number is 0.5. The space step h=0.01, ε=ε2h, 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.Figure 10
Temporal evolution ofϕ: (a), (b), and (c) are snapshots at t=0,Δt, and 2Δt, respectively.
(a)(b)(c)


### 3.3. Image Segmentation Application

Image segmentation plays a significant role in image unperceiving, computer vision, medical image processing, and pattern identification [32–36]. We want to partition an image into several regions and make each region have characteristics such as edges, color, and object. The 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 data-fitting 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:(14)Eϕ=∫Ωε∇ϕ2+1εsin2πϕ+GKϕ,f0dx.The first termε∇ϕ2 makes a smooth transition between phases and penalizes a large amount of oscillation. The 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. The 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. The third term is a fitting term, which is defined as follows:(15)GKϕ,f0=λ2∑k=0KCk−f02sinc2ϕ−k,where f0 is the given image, λ is a positive parameter which indicates degree of fidelity, and Ck is the average of f0 in the k-level, where k=0,1,…,K, i.e.,(16)Ck=∫Ωf0xsinc2ϕx−kdx∫Ωsinc2ϕx−kdx.Note that as the value ofλ increases, the influence of the edge-stopping function GK 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 sinc2ϕ. Note that the square of the sincϕ plays the same role as the delta function in k-level.Figure 11
Functionssincϕ and sinc2ϕ.Subsequently, we obtain the following gradient descent flow equation:(17)ϕt=2εΔϕ−πεsin2πϕ−λ∑k=0KCk−f02sin2πϕ−k2πϕ−k2−sin2πϕ−kπ2ϕ−k3,where we assume that Ck 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:(18)∂ϕx,t∂t=−λ∑k=0KCk−f02×sin2πϕx,t−k2πϕx,t−k2−sin2πϕx,t−kπ2ϕx,t−k3,∂ϕx,t∂t=−πεsin2πϕx,t,∂ϕx,t∂t=2εΔϕx,t.First, we solve the above three equations by applying the fully explicit method with a known valueϕn, that is,(19)ϕijn+1/3−ϕijnΔt=−λ∑k=0KCkn−f0,ij2×πϕijn−ksin2πϕijn−k−2sin2πϕijn−k2π2ϕijn−k3+δ,(20)ϕijn+2/3−ϕijn+1/3Δt=−πεsin2πϕijn+1/3,(21)ϕijn+1−ϕijn+2/3Δt=2ϵΔhϕijn+2/3,where we added a small value δ in the denominator to avoid singularities.Note that if we replaceϕijn+1 by ϕijn+2/3 on the right side of equation (21), it becomes an implicit method, which generally needs an iterative method to solve it.NowCkn is defined as follows:(22)Ckn=∑i=1Nx∑j=1Nyf0,ijsinc2ϕijn−k∑i=1Nx∑j=1Nysinc2ϕijn−k.In the numerical experiments of image segmentation, we normalize the selected imagef as f0=f−fmin/fmax−fmin, where fmin and fmax are the minimum and maximum values of the selected image, respectively. For the initial phase-field value, we use the following formula:(23)ϕij0=n−1+f0,ij−Tn−1Tn−Tn−1,iff0,ij∈Tn−1,Tn,where Tn for n=0,1,…,N are target levels of the image. For example, if we have four target values, T0=0,T1=0.3,T2=0.8, and T3=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–39] is a traditional method and has been used so far.Figure 12
Initialization of multiphase image segmentation.


### 3.4. Synthetic Image 1

We test image segmentation using the image with 5 phases of240×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=sminh2/8ε,ε/2π2 are used. Additionally, we use 5 target values: T0=0, T1=0.25, T2=0.5, T3=0.75, and T4=1. We evaluate the discrete value of the Mumford–Shah energy functional as follows:(24)Εϕn=∑i,jεϕi+1,jn−ϕijn2+ϕi,j+1n−ϕijn2h2+1εsin2πϕijn+λ2∑k=0KCkn−f0,ij2sinc2ϕijn−kh2,where Ckn is described earlier in equation (22).Figure 13
Image segmentation process with 5 phases. (a) Initial gray-scale image with 5 uniformly distributed target values. (b) Mesh plot of the initial image. (c) Converged phase field after 5 iterations. (d) Energy decrease plot and corresponding snapshots.
(a)(b)(c)(d)As shown in Figure13(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<tole. This threshold tole is set equal to 1e−4 for the tests with the synthetic image unless otherwise specified. Therefore, 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. Therefore, uniformly distributed target values used in Figure13 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 the5% 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. The 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.Figure 14
Image segmentation process of the noise-corrupted image with 5 phases. (a) Initial gray-scale image corrupted by5% random noise from the uniform distribution with 5 nonuniformly distributed target values. (b) Mesh plot of the initial image. (c) Converged phase field after 15 iterations. (d) Energy decrease plot and corresponding snapshots.
(a)(b)(c)(d)Subsequently, we generate a kind of impulsive noise, which is called salt and pepper noise, in the synthetic image. Figure15 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.Figure 15
Image segmentation process of the noise-corrupted image with 5 phases. (a) Initial gray-scale image corrupted by salt and pepper noise with 5 nonuniformly distributed target values. (b) Segmented image after 150 iterations. (c) Converged phase field after 150 iterations. (d) Energy decrease plot and corresponding snapshots.
(a)(b)(c)(d)Since relatively smalls is used, the number of iterations has to increase to achieve a similar level of segmentation as shown in Figure 14. The 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. The 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. This 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.Figure 16
Image segmentation process of noise-corrupted image with 5 phases. (a) Initial gray-scale image anisotropically corrupted by20% random noise from the standard normal distribution. We use 5 uniformly distributed target values with buffer. Note that the noise occurs in red boxes. (b) Segmented image after 150 iterations. (c) Initial gray-scale image corrupted by periodic noise with 5 nonuniformly distributed target values. (d) Segmented image after 150 iterations.
(a)(b)(c)(d)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.Figure 17
Energy plot of the Mumford–Shah energy functional obtained with proposed method equations (19)–(21) with the initial condition in Figure 14.According to Figure17, our method has long-term energy stability, which asserts nonoscillation asymptotic behaviors.Though 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. Figure18 shows original and predicted bounding boxes with respect to different ϵ values, and Table 2 lists the value of IoU for each case.Figure 18
Bounding boxes of both original and segmented images: (a)ε=εh, (b) ε=ε2h, (c) ε=ε4h, and (d) ε=ε8h. Note that the yellow boxes are original, and the green boxes are predicted.
(a)(b)(c)(d)Table 2
IoU for differentϵ values (note that the values are rounded up from the fifth decimal place).
εHH=hH=2hH=4hH=8hIoU0.98320.94300.89800.8980In this test, we use the parameter valuess=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. Therefore, IoU might be similar for ϵ values exceeding a certain scale in diffuse-interface models due to the lack of the number of pixels. This is the reason why we have similar values for ε4h and ε8h.Moreover, we verify the dullness of edge of our proposed method. Figure19 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. The 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.Figure 19
Dullness of edge of our proposed method with respect toεm, where m=h,2h,4h,8h.


### 3.5. Complex Image

For the following step, we use a more complex image. We fixtol=1 e−6 and tole=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. Therefore, 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 ik,jk for k=0,1,…,K, then we define Ck as the local mean of the near 8 cells and itself.(25)Ck=19∑a,b=−1,0,1f0,ik+a,jk+b.Note that we fix these defined values until completion of the segmentation process. Here, we uses=0.3, h=1, ε=ε2h, λ=1, Δt=sminh2/8ε,ε/2π2, and K=3. The size of image is 290×320. Figures 20(a) and 20(b) show the original image and the image with Ck positions for k=0,1,2,3, respectively. Figures 20(c) and 20(d) show mesh plots of the original brain MRI and the segmented phase field, respectively. Figures 20(e)–20(h) show the sectional segmentation images for each k-level.Figure 20
(a) Initial image of brain MRI. (b) Locations ofCk for k=0,1,2,3. (c) Mesh plot of initial image. (d) Segmentation results. Sectional image segmentation process of brain MRI: (e) ϕ=0, (f) ϕ=1, (g) ϕ=2, and (h) ϕ=3.
(a)(b)(c)(d)(e)(f)(g)(h)Note that the expression in black and white colors is just to clarify the visualization. Therefore, 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 meanCk is presented in Figure 21. Here, we use all the same parameters listed above for Figure 20, except for the number of iterations (500).Figure 21
Energy plot of the Mumford–Shah energy functional using our proposed method with local meanCk.


### 3.6. Parameter Test

We perform parameter tests on image with low contrast and high-level noise in this section. We employ the132×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.Figure 22
(a) Original image of vertebral MRI with low contrast near boundaries of object. (b) Perturbed image of vertebral MRI with±30% noise from the uniform distribution. (c) Segmented image using (b) with default parameters.
(a)(b)(c)We have targeted the spines within the image to segment, andK=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 tole=1 e−3, respectively.Figure 23
(a)–(c) verify the effect ofs on the image segmentation. (d)–(f) verify the effect of h on the image segmentation. (g)–(i) verify the effect of λ on the image segmentation. (j)–(l) verify the effect of ϵ on the image segmentation.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)As depicted in Figures23(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. The 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. Therefore, these have sharp-interface limits asymptotically and are smoothly smeared. Figures 23(j)–23(l) describe the effect of ϵ. The small value of ϵ makes the transition width to be narrow, resulting in jagged edges due to inconsistent separation depending on the shape of boundary while larger ϵ values make the boundaries too smooth, making them more blunt than the actual shapes.


### 3.7. 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. Figure24 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.Figure 24
(a)–(d) Segmented images using our proposed method. (e)–(h) Segmented images using implicit method followed by Li and Kim [27]. (i)–(l) Segmented images using quadratic regularization method followed by Kang and March [45]. (m)–(p) Segmented images using the superlinear model followed by Kang and March [45].
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)(m)(n)(o)(p)We adopt the parametersϵ=0.01, Δt=0.25h2, 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 Δtref=1e−5. We fix ε=ε2h for both the schemes. We use the following initial value ϕ0, ϕij0=Kf0,ij, where K=3,4 and target values as T0=0.01, T1=0.35, T2=0.66, T3=0.92 in the case of K=3 and T0=0.004, T1=0.23, T2=0.52, T3=0.77, T4=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 Δtref, 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. This result implies that we can achieve fast computation speed via the proposed method with the moderate time step proposed in Section 2.Table 3
Elapsed times taken by the proposed method and the reference method [27].
CasesProposed methodΔtProposed methodΔtrefReferenceΔtrefK=30.315 s2.297 s2.843 sK=40.403 s2.896 s3.062 s


## 3.1. Phase Evolution in One-Dimensional Computational Domain

Unless specified otherwise, we fixh=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 ϕm∞≤0.05<ϕm+1∞, such that xα=xm+h0.05−ϕm/ϕm+1−ϕm is shown in 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(b). Here, ϕ∞xα≈0.05, and ϕ∞xβ≈0.95. Figure 2(c) shows the numerical steady states with different ε=2h,4h,…,20h values. The initial condition in the computational domain Ω=0,1 is ϕx,0=0 if x<0.5, and ϕx,0=1 otherwise. The tolerance tol=1 e−16 is used.Figure 2
(a) Schematic illustration of linear interpolation evaluatingxα. (b) Schematic illustration of transition thickness, L=xβ−xα. (c) Equilibrium solutions with different ε=2h,4h,…,20h values.
(a)(b)(c)Table1 lists the thickness L=xβ−xα of the interface transition layer for various ϵ values.Table 1
ThicknessL=xβ−xα of interface transition layer for various ϵ values.
Caseϵ=2hϵ=4hϵ=6hϵ=8hϵ=10hL2.9903h6.5364h9.7668h12.9192h16.2230hCaseϵ=12hϵ=14hϵ=16hϵ=18hϵ=20hL19.4323h22.6640h25.9192h29.1198h32.3832hFigure3 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(11)εL=0.6158L+0.0004.Figure 3
Plot of interface transition thicknessL=xβ−xα with respect to ϵ.Note that we employ the linear regression based on the least-squares method in order to find such a linear function. Therefore, if we want to haveL as the interface thickness, then we can use equation (11) to define the value of ϵ.Figure4 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 4
Numerical steady states with initial conditionsϕx,0=0 when x<0.5, and ϕx,0=m otherwise. (a), (b), (c), and (d) are the results with m=1,2,3, and 4, respectively. ε=ε4h is used.
(a)(b)(c)(d)Figure5 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 5
Numerical steady states with initial conditionsϕx,0=0 when x<0.5, and ϕx,0=m otherwise. (a), (b), (c), and (d) are the results with m=1,2,3, and 4, respectively. ε=ε8h is used.
(a)(b)(c)(d)Figure6 shows the temporal evolution of ϕ with an initial condition on a computational domain Ω=0,1:(12)ϕx,0=2randx+5x−0.5,where rand x is a random number between −1 and 1. The parameters h=0.01, ε=ε4h, and Δt with s=0.1 are used. In Figure 6, the dotted, solid, and thick solid lines are the initial profile, intermediate, and final solution after 100 iterations, respectively. While initial profile and intermediate solutions have random values, we can observe some flat regions with integer values at the final solution.Figure 6
Temporal evolution of numerical solutions with initial condition,ϕx,0=2randx+5x−0.5.


## 3.2. Phase Evolution in Two-Dimensional Computational Domain

In this section, we begin our experiment in a two-dimensional 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. The initial condition in Ω is given by(13)ϕx,y,0=0.5randx,y+5x+y−1,where rand x,y is a random number between −1 and 1. ε=ε2h and Δt with s=0.99 are used.Figure 7
Temporal evolution of numerical solution: (a), (b), and (c) represent snapshots att=0,Δt, and 7Δt, respectively.
(a)(b)(c)Figure8 shows the temporal evolution of ϕ: Figures 8(a), 8(b), and 8(c) represent snapshots at t=0,Δt, and 4Δt, respectively. The initial condition is defined by perturbed integers from 0, 1, 3, and 6 in 2 by 2 subdomains of Ω. The maximum amplitude of the random number is 0.5. ε=ε2h and Δt with s=0.99 are used.Figure 8
Temporal evolution ofϕ with s=0.99: (a), (b), and (c) are snapshots at t=0,Δt, and 4Δt, respectively.
(a)(b)(c)In addition, we investigate the same problem in Figure8 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. Therefore, one might use the value of s as large as possible that is guaranteed to converge in order to gain fast convergence.Figure 9
Temporal evolution ofϕ with s=0.11: (a), (b), and (c) are the snapshots at t=0,9Δt, and 36Δt, respectively.
(a)(b)(c)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. The initial condition is defined by perturbed random integers between 0 and 26 in 5 by 5 subdomains of Ω. The maximum amplitude of the random number is 0.5. The space step h=0.01, ε=ε2h, 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.Figure 10
Temporal evolution ofϕ: (a), (b), and (c) are snapshots at t=0,Δt, and 2Δt, respectively.
(a)(b)(c)


## 3.3. Image Segmentation Application

Image segmentation plays a significant role in image unperceiving, computer vision, medical image processing, and pattern identification [32–36]. We want to partition an image into several regions and make each region have characteristics such as edges, color, and object. The 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 data-fitting 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:(14)Eϕ=∫Ωε∇ϕ2+1εsin2πϕ+GKϕ,f0dx.The first termε∇ϕ2 makes a smooth transition between phases and penalizes a large amount of oscillation. The 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. The 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. The third term is a fitting term, which is defined as follows:(15)GKϕ,f0=λ2∑k=0KCk−f02sinc2ϕ−k,where f0 is the given image, λ is a positive parameter which indicates degree of fidelity, and Ck is the average of f0 in the k-level, where k=0,1,…,K, i.e.,(16)Ck=∫Ωf0xsinc2ϕx−kdx∫Ωsinc2ϕx−kdx.Note that as the value ofλ increases, the influence of the edge-stopping function GK 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 sinc2ϕ. Note that the square of the sincϕ plays the same role as the delta function in k-level.Figure 11
Functionssincϕ and sinc2ϕ.Subsequently, we obtain the following gradient descent flow equation:(17)ϕt=2εΔϕ−πεsin2πϕ−λ∑k=0KCk−f02sin2πϕ−k2πϕ−k2−sin2πϕ−kπ2ϕ−k3,where we assume that Ck 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:(18)∂ϕx,t∂t=−λ∑k=0KCk−f02×sin2πϕx,t−k2πϕx,t−k2−sin2πϕx,t−kπ2ϕx,t−k3,∂ϕx,t∂t=−πεsin2πϕx,t,∂ϕx,t∂t=2εΔϕx,t.First, we solve the above three equations by applying the fully explicit method with a known valueϕn, that is,(19)ϕijn+1/3−ϕijnΔt=−λ∑k=0KCkn−f0,ij2×πϕijn−ksin2πϕijn−k−2sin2πϕijn−k2π2ϕijn−k3+δ,(20)ϕijn+2/3−ϕijn+1/3Δt=−πεsin2πϕijn+1/3,(21)ϕijn+1−ϕijn+2/3Δt=2ϵΔhϕijn+2/3,where we added a small value δ in the denominator to avoid singularities.Note that if we replaceϕijn+1 by ϕijn+2/3 on the right side of equation (21), it becomes an implicit method, which generally needs an iterative method to solve it.NowCkn is defined as follows:(22)Ckn=∑i=1Nx∑j=1Nyf0,ijsinc2ϕijn−k∑i=1Nx∑j=1Nysinc2ϕijn−k.In the numerical experiments of image segmentation, we normalize the selected imagef as f0=f−fmin/fmax−fmin, where fmin and fmax are the minimum and maximum values of the selected image, respectively. For the initial phase-field value, we use the following formula:(23)ϕij0=n−1+f0,ij−Tn−1Tn−Tn−1,iff0,ij∈Tn−1,Tn,where Tn for n=0,1,…,N are target levels of the image. For example, if we have four target values, T0=0,T1=0.3,T2=0.8, and T3=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–39] is a traditional method and has been used so far.Figure 12
Initialization of multiphase image segmentation.


## 3.4. Synthetic Image 1

We test image segmentation using the image with 5 phases of240×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=sminh2/8ε,ε/2π2 are used. Additionally, we use 5 target values: T0=0, T1=0.25, T2=0.5, T3=0.75, and T4=1. We evaluate the discrete value of the Mumford–Shah energy functional as follows:(24)Εϕn=∑i,jεϕi+1,jn−ϕijn2+ϕi,j+1n−ϕijn2h2+1εsin2πϕijn+λ2∑k=0KCkn−f0,ij2sinc2ϕijn−kh2,where Ckn is described earlier in equation (22).Figure 13
Image segmentation process with 5 phases. (a) Initial gray-scale image with 5 uniformly distributed target values. (b) Mesh plot of the initial image. (c) Converged phase field after 5 iterations. (d) Energy decrease plot and corresponding snapshots.
(a)(b)(c)(d)As shown in Figure13(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<tole. This threshold tole is set equal to 1e−4 for the tests with the synthetic image unless otherwise specified. Therefore, 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. Therefore, uniformly distributed target values used in Figure13 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 the5% 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. The 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.Figure 14
Image segmentation process of the noise-corrupted image with 5 phases. (a) Initial gray-scale image corrupted by5% random noise from the uniform distribution with 5 nonuniformly distributed target values. (b) Mesh plot of the initial image. (c) Converged phase field after 15 iterations. (d) Energy decrease plot and corresponding snapshots.
(a)(b)(c)(d)Subsequently, we generate a kind of impulsive noise, which is called salt and pepper noise, in the synthetic image. Figure15 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.Figure 15
Image segmentation process of the noise-corrupted image with 5 phases. (a) Initial gray-scale image corrupted by salt and pepper noise with 5 nonuniformly distributed target values. (b) Segmented image after 150 iterations. (c) Converged phase field after 150 iterations. (d) Energy decrease plot and corresponding snapshots.
(a)(b)(c)(d)Since relatively smalls is used, the number of iterations has to increase to achieve a similar level of segmentation as shown in Figure 14. The 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. The 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. This 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.Figure 16
Image segmentation process of noise-corrupted image with 5 phases. (a) Initial gray-scale image anisotropically corrupted by20% random noise from the standard normal distribution. We use 5 uniformly distributed target values with buffer. Note that the noise occurs in red boxes. (b) Segmented image after 150 iterations. (c) Initial gray-scale image corrupted by periodic noise with 5 nonuniformly distributed target values. (d) Segmented image after 150 iterations.
(a)(b)(c)(d)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.Figure 17
Energy plot of the Mumford–Shah energy functional obtained with proposed method equations (19)–(21) with the initial condition in Figure 14.According to Figure17, our method has long-term energy stability, which asserts nonoscillation asymptotic behaviors.Though 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. Figure18 shows original and predicted bounding boxes with respect to different ϵ values, and Table 2 lists the value of IoU for each case.Figure 18
Bounding boxes of both original and segmented images: (a)ε=εh, (b) ε=ε2h, (c) ε=ε4h, and (d) ε=ε8h. Note that the yellow boxes are original, and the green boxes are predicted.
(a)(b)(c)(d)Table 2
IoU for differentϵ values (note that the values are rounded up from the fifth decimal place).
εHH=hH=2hH=4hH=8hIoU0.98320.94300.89800.8980In this test, we use the parameter valuess=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. Therefore, IoU might be similar for ϵ values exceeding a certain scale in diffuse-interface models due to the lack of the number of pixels. This is the reason why we have similar values for ε4h and ε8h.Moreover, we verify the dullness of edge of our proposed method. Figure19 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. The 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.Figure 19
Dullness of edge of our proposed method with respect toεm, where m=h,2h,4h,8h.


## 3.5. Complex Image

For the following step, we use a more complex image. We fixtol=1 e−6 and tole=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. Therefore, 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 ik,jk for k=0,1,…,K, then we define Ck as the local mean of the near 8 cells and itself.(25)Ck=19∑a,b=−1,0,1f0,ik+a,jk+b.Note that we fix these defined values until completion of the segmentation process. Here, we uses=0.3, h=1, ε=ε2h, λ=1, Δt=sminh2/8ε,ε/2π2, and K=3. The size of image is 290×320. Figures 20(a) and 20(b) show the original image and the image with Ck positions for k=0,1,2,3, respectively. Figures 20(c) and 20(d) show mesh plots of the original brain MRI and the segmented phase field, respectively. Figures 20(e)–20(h) show the sectional segmentation images for each k-level.Figure 20
(a) Initial image of brain MRI. (b) Locations ofCk for k=0,1,2,3. (c) Mesh plot of initial image. (d) Segmentation results. Sectional image segmentation process of brain MRI: (e) ϕ=0, (f) ϕ=1, (g) ϕ=2, and (h) ϕ=3.
(a)(b)(c)(d)(e)(f)(g)(h)Note that the expression in black and white colors is just to clarify the visualization. Therefore, 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 meanCk is presented in Figure 21. Here, we use all the same parameters listed above for Figure 20, except for the number of iterations (500).Figure 21
Energy plot of the Mumford–Shah energy functional using our proposed method with local meanCk.


## 3.6. Parameter Test

We perform parameter tests on image with low contrast and high-level noise in this section. We employ the132×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.Figure 22
(a) Original image of vertebral MRI with low contrast near boundaries of object. (b) Perturbed image of vertebral MRI with±30% noise from the uniform distribution. (c) Segmented image using (b) with default parameters.
(a)(b)(c)We have targeted the spines within the image to segment, andK=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 tole=1 e−3, respectively.Figure 23
(a)–(c) verify the effect ofs on the image segmentation. (d)–(f) verify the effect of h on the image segmentation. (g)–(i) verify the effect of λ on the image segmentation. (j)–(l) verify the effect of ϵ on the image segmentation.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)As depicted in Figures23(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. The 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. Therefore, these have sharp-interface limits asymptotically and are smoothly smeared. Figures 23(j)–23(l) describe the effect of ϵ. The small value of ϵ makes the transition width to be narrow, resulting in jagged edges due to inconsistent separation depending on the shape of boundary while larger ϵ values make the boundaries too smooth, making them more blunt than the actual shapes.


## 3.7. 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. Figure24 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.Figure 24
(a)–(d) Segmented images using our proposed method. (e)–(h) Segmented images using implicit method followed by Li and Kim [27]. (i)–(l) Segmented images using quadratic regularization method followed by Kang and March [45]. (m)–(p) Segmented images using the superlinear model followed by Kang and March [45].
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)(m)(n)(o)(p)We adopt the parametersϵ=0.01, Δt=0.25h2, 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 Δtref=1e−5. We fix ε=ε2h for both the schemes. We use the following initial value ϕ0, ϕij0=Kf0,ij, where K=3,4 and target values as T0=0.01, T1=0.35, T2=0.66, T3=0.92 in the case of K=3 and T0=0.004, T1=0.23, T2=0.52, T3=0.77, T4=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 Δtref, 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. This result implies that we can achieve fast computation speed via the proposed method with the moderate time step proposed in Section 2.Table 3
Elapsed times taken by the proposed method and the reference method [27].
CasesProposed methodΔtProposed methodΔtrefReferenceΔtrefK=30.315 s2.297 s2.843 sK=40.403 s2.896 s3.062 s


## 4. Geometric Active Segmentation

For approximation of the geometric active segmentation for multiphase image, the governing equation is given as follows:(26)∂ϕx,t∂t=gf0x2εΔϕx,t−πεsin2πϕx,t+λg2f0xϕx,t1εsin2πϕx,t,x∈Ω,0<t≤T,where f0x for x=x,y is a given image, g is an edge-stopping function, ϕx,t is an order parameter, ϵ is a positive constant, and λ is a constant parameter. Ω, in particular, is a domain bounded in ℝdd=1,2. Note that we restrict ϕx,t≥0 for convenience. In equation (26), gf0x is the edge-stopping function and is defined as(27)gf0x=11+∇Gσ∗f0xp,where Gσ∗f0x is the convolution of the given image f0 and the Gaussian function(28)Gσ=12πσ2e−x2+y2/2σ2.Here, we will usep=2. The function gf0x 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:(29)∂ϕx,t∂t=λg2f0xϕx,t1εsin2πϕx,t,∂ϕx,t∂t=−πεgf0xsin2πϕx,t,∂ϕx,t∂t=2εgf0xΔϕx,t.First, we solve the above three equations by applying the fully explicit method with a known valueϕn, that is,(30)ϕijn+1/3−ϕijnΔt=λgij2ϕijn1εsin2πϕijn,ϕijn+2/3−ϕijn+1/3Δt=−πεgijsin2πϕijn+1/3,ϕijn+1−ϕijn+2/3Δt=2εgijΔhϕijn+2/3.Also, the edge functiongf0 with k=2 can be numerically computed by(31)gf0ij=11+Gσ∗f0x,ij2+Gσ∗f0y,ij2,where Gσ∗f0x,ij=Gσ∗f0i+1,j−Gσ∗f0i−1,j/2h and Gσ∗f0y,ij=Gσ∗f0i,j+1−Gσ∗f0i,j−1/2h. Here, Gσ∗f0ij can be computed using a 3×3 smoothing kernel as(32)Gσ∗f0ij=∑m=i−1i+1∑n=j−1j+1f0mn2πσ2e−i−m2+j−n2h2/2σ2.As mentioned above, the homogeneous Neumann boundary condition is applied to the domain.


### 4.1. Synthetic Image 2

We apply our proposed method to a synthetic image. The synthetic image consists of three rectangles and one circle. That is, we try to separate into four phases. The initial contour is set to the median of two consecutive numbers of0,1,2,3, and 4. The size of the image is 256×256. The parameters are chosen as ε=ε3h, Δt=5 e−6, T=10000Δt, and λ=1 e2. The results of the synthetic image segmentation can be seen in Figure 25. Figure 25(a) shows the initial contour-shaped rectangular. While Figure 25(b) shows the result at t=2000Δt, Figure 25(c) shows the result at final time, T=10000Δt. Figure 25(d) shows the result of multiphase image segmentation in a mesh plot.Figure 25
Snapshots of contours are shown in (a), (b), and (c) att=0, t=2000Δt, and t=10000Δt, respectively. Mesh plot of (c) is shown in (d).
(a)(b)(c)(d)


### 4.2. 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.Figure 26
Segmentation with the real-world images chosen from the CMU-Cornell iCoseg Dataset: (a) color images, (b) results of the proposed method, and (c) ground truth of the color images.
(a)(b)(c)


### 4.3. Measurement of Quality of the Segmentation

As quantitative measures of the ability to restore damaged fingerprints to the proposed algorithm, we use peak signal-to-noise ratio (PSNR) and structural similarity index map (SSIM), which are defined as follows:(33)PSNR=10log10Imax2∑i=1Nx∑j=1NyIij−Kij2/NxNy,SSIM=2μIμK+c12σIK+c2μI2+μK2+c1σI2+σK2+c2,where Iij=fij/fmax−fmin, Kij=0.5ϕij+1; μI,μK denote the average of I and K, respectively; and σI2,σK2 are the variances of I and K, respectively. Here, σIK is the covariance of I and K, c1=0.01L2,andc2=0.03L2 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 Figure27. From the result, we can see that the curve is increasing, which indicates that the proposed method works well for the image segmentation.Figure 27
Time evolution of a quantitative measure of the segmentation of the proposed method.


### 4.4. X-Ray Example

We perform the multiphase image segmentation using a knee X-ray as the original image. The size of given image is128×128. The original image and scalar fields of function g and g2 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.Figure 28
(a) Original image. (b) Edge-stopping functiong. (c) Edge-stopping function g2.
(a)(b)(c)Figure 29
Snapshots of contours at (a)t=0, (b) t=500000Δt, and (c) t=1400000Δt. Mesh plot of (c) is shown in (d).
(a)(b)(c)(d)The results yield that the edge-stopping functionsg and g2 work well and one can easily extract each activated contour only using multiple-well potential.


## 4.1. Synthetic Image 2

We apply our proposed method to a synthetic image. The synthetic image consists of three rectangles and one circle. That is, we try to separate into four phases. The initial contour is set to the median of two consecutive numbers of0,1,2,3, and 4. The size of the image is 256×256. The parameters are chosen as ε=ε3h, Δt=5 e−6, T=10000Δt, and λ=1 e2. The results of the synthetic image segmentation can be seen in Figure 25. Figure 25(a) shows the initial contour-shaped rectangular. While Figure 25(b) shows the result at t=2000Δt, Figure 25(c) shows the result at final time, T=10000Δt. Figure 25(d) shows the result of multiphase image segmentation in a mesh plot.Figure 25
Snapshots of contours are shown in (a), (b), and (c) att=0, t=2000Δt, and t=10000Δt, respectively. Mesh plot of (c) is shown in (d).
(a)(b)(c)(d)


## 4.2. 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.Figure 26
Segmentation with the real-world images chosen from the CMU-Cornell iCoseg Dataset: (a) color images, (b) results of the proposed method, and (c) ground truth of the color images.
(a)(b)(c)


## 4.3. Measurement of Quality of the Segmentation

As quantitative measures of the ability to restore damaged fingerprints to the proposed algorithm, we use peak signal-to-noise ratio (PSNR) and structural similarity index map (SSIM), which are defined as follows:(33)PSNR=10log10Imax2∑i=1Nx∑j=1NyIij−Kij2/NxNy,SSIM=2μIμK+c12σIK+c2μI2+μK2+c1σI2+σK2+c2,where Iij=fij/fmax−fmin, Kij=0.5ϕij+1; μI,μK denote the average of I and K, respectively; and σI2,σK2 are the variances of I and K, respectively. Here, σIK is the covariance of I and K, c1=0.01L2,andc2=0.03L2 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 Figure27. From the result, we can see that the curve is increasing, which indicates that the proposed method works well for the image segmentation.Figure 27
Time evolution of a quantitative measure of the segmentation of the proposed method.


## 4.4. X-Ray Example

We perform the multiphase image segmentation using a knee X-ray as the original image. The size of given image is128×128. The original image and scalar fields of function g and g2 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.Figure 28
(a) Original image. (b) Edge-stopping functiong. (c) Edge-stopping function g2.
(a)(b)(c)Figure 29
Snapshots of contours at (a)t=0, (b) t=500000Δt, and (c) t=1400000Δt. Mesh plot of (c) is shown in (d).
(a)(b)(c)(d)The results yield that the edge-stopping functionsg and g2 work well and one can easily extract each activated contour only using multiple-well potential.


## 5. 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. The results indicate that proper phase separations are achieved via our simple explicit method. The 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. Therefore, 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].



---

*Source: 2414209-2021-11-20.xml*