# Curve and Surface Smoothing Using a Modified Cahn-Hilliard Equation

Authors: Yongho Choi; Darea Jeong; Junseok Kim
Journal: Mathematical Problems in Engineering
 (2017)

---

## Abstract

We present a new method using the modified Cahn-Hilliard (CH) equation for smoothing piecewise linear shapes of two- and three-dimensional objects. The CH equation has good smoothing dynamics and it is coupled with a fidelity term which keeps the original given data; that is, it does not produce significant shrinkage. The modified CH equation is discretized using a linearly stable splitting scheme in time and the resulting scheme is solved by using a Fourier spectral method. We present computational results for both curve and surface smoothing problems. The computational results demonstrate that the proposed algorithm is fast and efficient.


---

## Body

## 1. Introduction

For many computational purposes such as visualization of scientific data and registration of multimodal medical data, smooth curves and surfaces must be approximated by polygonal curves and polyhedral surfaces, respectively. However, these approximation algorithms generate faceted curves and surfaces [1]. In medical applications, Vollmer et al. [2] presented an improved Laplacian technique for smoothing polygonal surface meshes that avoids the well-known problem of deformation and shrinkage caused by many smoothing methods as shown in Figure 1. Using Bezier-smoothed boundary representation, Anflor and Marczak [3] studied topological optimization of anisotropic heat conducting devices.Figure 1
(a)–(e) Origin, noisy mesh, three iterations of the algorithm in [2]. Reprinted from Vollmer et al. [2] with permission from the Wiley Online Library.
(a)
(b)
(c)
(d)
(e)To prevent the mesh from shrinking, Liu et al. [4] proposed a volume-constrained smoothing algorithm for triangular meshes, which preserves exactly the mesh volume during the smoothing process. Wei et al. [5] proposed an effective mesh smoothing algorithm which preserves morphology on polygonized isosurfaces of inhomogeneous binary volumes. Recently, various smoothing methods have been developed [6, 7].In this paper, we propose a new method using a modified Cahn-Hilliard (CH) equation for smoothing piecewise linear shapes of two- and three-dimensional objects. Letψ(x) be a given mosaic binary data, that is, either plus one or minus one. Here, x=(x,y) or x=(x,y,z). We want to construct smooth contour line or isosurface from the mosaic binary data by using the modified CH equation: (1)ϕtx,t=ΔF′ϕx,t-ε2Δϕx,t+λψx-ϕx,t,where λ is the fidelity parameter and Fϕ=0.25ϕ2-12. When λ=0, then (1) becomes the CH equation which was originally proposed by Cahn and Hilliard [8, 9] as a phenomenological model of phase separation in a binary alloy. See [10] and the references therein for physical, mathematical, and numerical derivations of the CH equation. The various applications of the CH equation can be seen in [11]. Equation (1) was used in image inpainting with a space-dependent λx [12].This paper is organized as follows. In Section2, the discrete equations for the governing equations are presented. In Section 3, we present various computational examples. Finally, in Section 4, conclusions are drawn.


## 2. Algorithm

We use an unconditionally stable Fourier spectral method for the modified CH equation (1) in two-dimensional space Ω=a,b×c,d. Fourier spectral method has been applied in many studies [13, 14]. The three-dimensional discretization is defined analogously. Let h=b-a/Nx=d-c/Ny be the spatial step size, where Nx and Ny are positive even integers. Let ϕmnk and μmnk be approximations of ϕxm,yn,tk and μxm,yn,tk, respectively, where xm,yn=a+m-0.5h,b+n-0.5h and tk=kΔt with Δt the temporal step. Let ϕ^pqk=∑m=1Nx∑n=1Nyϕmnke-iξpxm+ηqyn be the discrete Fourier transform with a given data ϕmnk∣m=1,…,Nxandn=1,…,Ny, where ξp=2πp/b-a and ηq=2πq/d-c. The inverse discrete Fourier transform is (2)ϕmnk=1NxNy∑p=-Nx/2Nx/2-1∑q=-Ny/2Ny/2-1ϕ^pqkeiξpxm+ηqyn.Letϕx,y,kΔt  =  (1/NxNy)∑p=-Nx/2Nx/2-1∑q=-Ny/2Ny/2-1ϕ^pqkei(ξpx+ηqy). Then, we have(3)Δϕx,y,kΔt=-1NxNy∑p=-Nx/2Nx/2-1∑q=-Ny/2Ny/2-1ξp2+ηq2ϕ^pqkeiξpx+ηqy.By using the linearly stabilized splitting scheme [15] to (1), we have (4)ϕmnk+1-ϕmnkΔt=Δ2ϕmnk+1-ε2Δϕmnk+1+fϕmnk+λψmn-ϕmnk+1,where fϕ=ϕ3-3ϕ. Thus, (4) can be transformed into the discrete Fourier space as follows: (5)ϕ^pqk+1-ϕ^pqkΔt=-ξp2+ηq22ϕ^pqk+1+ε2ξp2+ηq2ϕ^pqk+1+f^pqk+λψ^pq-ϕ^pqk+1.Therefore, we obtain the following discrete Fourier transform:(6)ϕ^pqk+1=ϕ^pqk-ξp2+ηq2Δtf^pqk+λψ^pq1+Δtλ+2ξp2+ηq2+ε2ξp2+ηq22.Then, using (2), we get the updated numerical solution ϕmnk+1 as (7)ϕmnk+1=1NxNy∑p=-Nx/2Nx/2-1∑q=-Ny/2Ny/2-1ϕ^pqk+1eiξpxm+ηqyn.Standard fast Fourier transform (FFT) has been used to speed up the computations [16].


## 3. Computational Results and Discussion

### 3.1. Circular Mosaic Shape

In this test problem, we will show the smoothing effect of our proposed method. Therefore, we choose a coarse mesh grid20×20 on the computational domain Ω=-1,1×-1,1. The initial configuration is(8)ϕx,y,0=1ifx2+y2≤0.6,-1otherwise,which represents a circular mosaic shape. The numerical parameters h=0.1,ε=h,Δt=0.1,andλ=1000 are used.Figure2 shows the temporal evolution at t=0, Δt, and 10Δt. Here, the shaded area is the filled contour with ϕ≥0. Figures 2(d)–2(f) show the values of ϕx,y at the part of the red box in Figures 2(a)–2(c). As shown in Figure 2(c), we can see the smooth interface after 10 iterations.Figure 2
(a) Initial circular mosaic shape. (b)-(c) Numerical results after1 and 10 iterations by the proposed method, respectively. (d)–(f) Values of ϕx,y at the red box.
(a)
(b)
(c)
(d)
(e)
(f)


### 3.2. Effect ofλ

Next, we investigate the effect of parameterλ on the smoothing dynamics. On the computational domain Ω=-1,1×-1,1, the initial configuration is given as (9)ϕx,y,0=1ifx2+y2+0.3cos⁡6θ≤0.5,-1otherwise,where θ=tan-1⁡y/x if x>0;θ=π+tan-1⁡y/x otherwise. We use the following parameters: h=2/30, ε=h, and Δt=0.1. With the initial configuration as shown in Figure 3(a), we obtain the various numerical results (see Figures 3(b)–3(d)) at time T=20Δt with respect to different λ. In this test, we use three different λ=106,104, and 102. As λ is large, ϕ(x,y) does not evolve and is close to the initial mosaic profile. However, if λ is too small, then the CH dynamics dominates and the profile becomes a circular shape which has the minimum interfacial length.Figure 3
Effect ofλ on the smoothing dynamics. (a) Initial condition and numerical results at T=20Δt with (b) λ=106, (c) λ=104, and (d) λ=102.
(a)
(b)
(c)
(d)


### 3.3. Effect ofε

In our proposed model, parameterε determines the thickness of the interface of two phases, that is, interfacial region of transition profile. Therefore, it is one of the important factors on the smoothing process. Now, we study the effect of ε on the smoothing dynamics.On the computational domainΩ=-1,1×-1,1, we use the initial condition (9); see Figure 4(a). The other parameters are used as h=2/30,Δt=0.1,andλ=10000.Figure 4
Effect ofε on the smoothing dynamics. (a) Initial condition and numerical results at T=20Δt with (b) ε=4h, (c) ε=h, and (d) ε=h/4.
(a)
(b)
(c)
(d)Figures4(b)–4(d) represent the numerical configurations at T=20Δt with ε=4h,ε=h,andε=h/4, respectively. As shown in Figure 4(b), large value of ε causes wide transition layer. Therefore, it overly smooths the curve. On the other hand, if ε is too small, then the phase-field has abrupt transition and the profile becomes mosaic.


### 3.4. 2D Various Shapes

To demonstrate the robustness of the proposed algorithm, we consider various complex shapes on two-dimensional spaceΩ=-1,1×-1,1. For the numerical test, we use N=30, h=2/N, ε=h, λ=10000, and Δt=0.1. The initial configurations are represented in Figure 5(a) and its corresponding smooth results in Figure 5(b). Through these results, we can see that our proposed algorithm is a good strategy to remove staircases while preserving the initial feature.Figure 5
(a) Three initial configurations. (b) Corresponding smooth results after 10 iterations.
(a)
(b)


### 3.5. 3D Various Shapes

In this section, we consider various morphologies on the three-dimensional space.


### 3.6. Ellipsoidal Shape

We consider an ellipsoidal shape on the computational domainΩ=-3,3×-1.5,1.5×-1.5,1.5. The initial condition is set to (10)ϕx,y,z,0=1ifx24+y2+z2≤1,-1otherwise.In the test problem, we useh=0.075,ε=5h,Δt=0.1,andλ=1000 on an 80×40×40 mesh grid. Figures 6(a)–6(c) show the temporal evolutions of the surface at T=0,Δt,and10Δt. As shown in Figure 6, the proposed method leads to a smooth ellipsoidal surface.Figure 6
(a) Initial condition. Numerical results after (b) 1 and (c) 10 iterations.
(a)
(b)
(c)


### 3.7. Volume Preservation

Laplacian smoothing is a well-known algorithm to smooth a polygonal mesh. However, this method can cause volume shrinkage or shape distortion [5]. However, our proposed method is a feature-preserving approach without severe volume changes. In this section, we demonstrate that our smoothing approach preserves the initial volume. We first define the volume of an object using triangles on the surface. For s=1,2,…,M, let Tris=Xl,Xm,Xn=Xl,Yl,Zl,Xm,Ym,Zm,Xn,Yn,Zn be a triangle consisting of three points on a triangular mesh. Here, M is the number of triangles. We have the following formula for the volume of the polyhedron as (11)VX=16∑s=1NXnYlZm-YmZl-YnXlZm-XmZl+ZnXlYm-XmYl.See [17] for more detailed explanation about the polyhedron volume. Also, we define the percent of relative error of volume VX with the initial volume VX0 as (12)VerrX=VX-VX0VX0×100%.Now, we consider three-dimensional morphologies such as to holed tetrahedron, many spheres, and holed sphere onΩ=-1,1×-1,1×-1,1; see Figure 7(a). On a 40×40×40 mesh grid, we perform the numerical test with h=0.1, ε=h, Δt=0.1, and λ=5000. We obtain smooth results after 10 iteration steps by our proposed method as shown in Figure 7(b). Also, to show the effectiveness of our method, we compute the temporal evolutions of the relative percent error; see Figure 7(c).Figure 7
(a) Three different initial conditions. (b) Smoothed results after 10 iterations by the proposed method. (c) Percentage of relative error of volume against iteration steps.
(a)
(b)
(c)Next, to show the robustness of the proposed method, we perform surface smoothing problem with noise in Figure8. The initial conditions are the same as in Figure 7 with 30 percent noise. The other parameters are the same as the previous test. Figure 8 also shows the smooth numerical results. Even with noise, our smoothing approach has robust and effective results.Figure 8
(a) Three different initial condition with30% noise. (b) Smoothed results after 10 iterations by the proposed method.
(a)
(b)


#### 3.7.1. Complex 3D Morphology

As the final test, we take Armadillo, Stanford bunny, and happy Buddha models as the initial conditions. For all tests, the parametersλ=104, Δt=0.01, and ε=h are used. The others are stated in Table 1.Table 1
Domain size, mesh grid points, and spatial step size used in Armadillo, stanford bunny, and happy buddha models.
Case
(a) Armadillo
(b) Stanford bunny
(c) Happy buddha
Domain size
0,2
×
0,2.4
×
0,2
0,1
×
0,1
×
0,1
0,1
×
0,2
.
2
×
0,1
Mesh grid points
50
×
60
×
50
60×60×50
100
×
220
×
100
Spatial step size
h=0.04
h=0.0167
h=0.01The columns in Figure9 represent the initial condition and numerical results after 1 and 3 iterations. As we expected, the numerical results are smooth and preserve the initial feature.Figure 9
Smoothing morphologies atT=0 (first column), T=Δt (second column), and T=3Δt (third column). Here, (a) Armadillo, (b) Stanford bunny, and (c) happy Buddha models are used.
(a)
(b)
(c)To show the effectiveness of the proposed smoothing approach, we compare results in [5] with those from our proposed method. Figure 10 illustrates CPU time in seconds against the total number of vertices of a model. As shown in Figure 10, under the similar number of vertex points, our proposed method gets smooth results and a reduction of over two orders of magnitude in CPU time.Figure 10
Total number of vertices versus CPU time.Also, in Figure11, we represent the percent of the relative error of volume VerrX which is stated in (12). Through these results, we can see that the proposed method implements the smoothing procedure almost without total volume loss.Figure 11
Temporal percent of relative error of volume with initial volume of Figure9.


## 3.1. Circular Mosaic Shape

In this test problem, we will show the smoothing effect of our proposed method. Therefore, we choose a coarse mesh grid20×20 on the computational domain Ω=-1,1×-1,1. The initial configuration is(8)ϕx,y,0=1ifx2+y2≤0.6,-1otherwise,which represents a circular mosaic shape. The numerical parameters h=0.1,ε=h,Δt=0.1,andλ=1000 are used.Figure2 shows the temporal evolution at t=0, Δt, and 10Δt. Here, the shaded area is the filled contour with ϕ≥0. Figures 2(d)–2(f) show the values of ϕx,y at the part of the red box in Figures 2(a)–2(c). As shown in Figure 2(c), we can see the smooth interface after 10 iterations.Figure 2
(a) Initial circular mosaic shape. (b)-(c) Numerical results after1 and 10 iterations by the proposed method, respectively. (d)–(f) Values of ϕx,y at the red box.
(a)
(b)
(c)
(d)
(e)
(f)


## 3.2. Effect ofλ

Next, we investigate the effect of parameterλ on the smoothing dynamics. On the computational domain Ω=-1,1×-1,1, the initial configuration is given as (9)ϕx,y,0=1ifx2+y2+0.3cos⁡6θ≤0.5,-1otherwise,where θ=tan-1⁡y/x if x>0;θ=π+tan-1⁡y/x otherwise. We use the following parameters: h=2/30, ε=h, and Δt=0.1. With the initial configuration as shown in Figure 3(a), we obtain the various numerical results (see Figures 3(b)–3(d)) at time T=20Δt with respect to different λ. In this test, we use three different λ=106,104, and 102. As λ is large, ϕ(x,y) does not evolve and is close to the initial mosaic profile. However, if λ is too small, then the CH dynamics dominates and the profile becomes a circular shape which has the minimum interfacial length.Figure 3
Effect ofλ on the smoothing dynamics. (a) Initial condition and numerical results at T=20Δt with (b) λ=106, (c) λ=104, and (d) λ=102.
(a)
(b)
(c)
(d)


## 3.3. Effect ofε

In our proposed model, parameterε determines the thickness of the interface of two phases, that is, interfacial region of transition profile. Therefore, it is one of the important factors on the smoothing process. Now, we study the effect of ε on the smoothing dynamics.On the computational domainΩ=-1,1×-1,1, we use the initial condition (9); see Figure 4(a). The other parameters are used as h=2/30,Δt=0.1,andλ=10000.Figure 4
Effect ofε on the smoothing dynamics. (a) Initial condition and numerical results at T=20Δt with (b) ε=4h, (c) ε=h, and (d) ε=h/4.
(a)
(b)
(c)
(d)Figures4(b)–4(d) represent the numerical configurations at T=20Δt with ε=4h,ε=h,andε=h/4, respectively. As shown in Figure 4(b), large value of ε causes wide transition layer. Therefore, it overly smooths the curve. On the other hand, if ε is too small, then the phase-field has abrupt transition and the profile becomes mosaic.


## 3.4. 2D Various Shapes

To demonstrate the robustness of the proposed algorithm, we consider various complex shapes on two-dimensional spaceΩ=-1,1×-1,1. For the numerical test, we use N=30, h=2/N, ε=h, λ=10000, and Δt=0.1. The initial configurations are represented in Figure 5(a) and its corresponding smooth results in Figure 5(b). Through these results, we can see that our proposed algorithm is a good strategy to remove staircases while preserving the initial feature.Figure 5
(a) Three initial configurations. (b) Corresponding smooth results after 10 iterations.
(a)
(b)


## 3.5. 3D Various Shapes

In this section, we consider various morphologies on the three-dimensional space.


## 3.6. Ellipsoidal Shape

We consider an ellipsoidal shape on the computational domainΩ=-3,3×-1.5,1.5×-1.5,1.5. The initial condition is set to (10)ϕx,y,z,0=1ifx24+y2+z2≤1,-1otherwise.In the test problem, we useh=0.075,ε=5h,Δt=0.1,andλ=1000 on an 80×40×40 mesh grid. Figures 6(a)–6(c) show the temporal evolutions of the surface at T=0,Δt,and10Δt. As shown in Figure 6, the proposed method leads to a smooth ellipsoidal surface.Figure 6
(a) Initial condition. Numerical results after (b) 1 and (c) 10 iterations.
(a)
(b)
(c)


## 3.7. Volume Preservation

Laplacian smoothing is a well-known algorithm to smooth a polygonal mesh. However, this method can cause volume shrinkage or shape distortion [5]. However, our proposed method is a feature-preserving approach without severe volume changes. In this section, we demonstrate that our smoothing approach preserves the initial volume. We first define the volume of an object using triangles on the surface. For s=1,2,…,M, let Tris=Xl,Xm,Xn=Xl,Yl,Zl,Xm,Ym,Zm,Xn,Yn,Zn be a triangle consisting of three points on a triangular mesh. Here, M is the number of triangles. We have the following formula for the volume of the polyhedron as (11)VX=16∑s=1NXnYlZm-YmZl-YnXlZm-XmZl+ZnXlYm-XmYl.See [17] for more detailed explanation about the polyhedron volume. Also, we define the percent of relative error of volume VX with the initial volume VX0 as (12)VerrX=VX-VX0VX0×100%.Now, we consider three-dimensional morphologies such as to holed tetrahedron, many spheres, and holed sphere onΩ=-1,1×-1,1×-1,1; see Figure 7(a). On a 40×40×40 mesh grid, we perform the numerical test with h=0.1, ε=h, Δt=0.1, and λ=5000. We obtain smooth results after 10 iteration steps by our proposed method as shown in Figure 7(b). Also, to show the effectiveness of our method, we compute the temporal evolutions of the relative percent error; see Figure 7(c).Figure 7
(a) Three different initial conditions. (b) Smoothed results after 10 iterations by the proposed method. (c) Percentage of relative error of volume against iteration steps.
(a)
(b)
(c)Next, to show the robustness of the proposed method, we perform surface smoothing problem with noise in Figure8. The initial conditions are the same as in Figure 7 with 30 percent noise. The other parameters are the same as the previous test. Figure 8 also shows the smooth numerical results. Even with noise, our smoothing approach has robust and effective results.Figure 8
(a) Three different initial condition with30% noise. (b) Smoothed results after 10 iterations by the proposed method.
(a)
(b)


### 3.7.1. Complex 3D Morphology

As the final test, we take Armadillo, Stanford bunny, and happy Buddha models as the initial conditions. For all tests, the parametersλ=104, Δt=0.01, and ε=h are used. The others are stated in Table 1.Table 1
Domain size, mesh grid points, and spatial step size used in Armadillo, stanford bunny, and happy buddha models.
Case
(a) Armadillo
(b) Stanford bunny
(c) Happy buddha
Domain size
0,2
×
0,2.4
×
0,2
0,1
×
0,1
×
0,1
0,1
×
0,2
.
2
×
0,1
Mesh grid points
50
×
60
×
50
60×60×50
100
×
220
×
100
Spatial step size
h=0.04
h=0.0167
h=0.01The columns in Figure9 represent the initial condition and numerical results after 1 and 3 iterations. As we expected, the numerical results are smooth and preserve the initial feature.Figure 9
Smoothing morphologies atT=0 (first column), T=Δt (second column), and T=3Δt (third column). Here, (a) Armadillo, (b) Stanford bunny, and (c) happy Buddha models are used.
(a)
(b)
(c)To show the effectiveness of the proposed smoothing approach, we compare results in [5] with those from our proposed method. Figure 10 illustrates CPU time in seconds against the total number of vertices of a model. As shown in Figure 10, under the similar number of vertex points, our proposed method gets smooth results and a reduction of over two orders of magnitude in CPU time.Figure 10
Total number of vertices versus CPU time.Also, in Figure11, we represent the percent of the relative error of volume VerrX which is stated in (12). Through these results, we can see that the proposed method implements the smoothing procedure almost without total volume loss.Figure 11
Temporal percent of relative error of volume with initial volume of Figure9.


## 3.7.1. Complex 3D Morphology

As the final test, we take Armadillo, Stanford bunny, and happy Buddha models as the initial conditions. For all tests, the parametersλ=104, Δt=0.01, and ε=h are used. The others are stated in Table 1.Table 1
Domain size, mesh grid points, and spatial step size used in Armadillo, stanford bunny, and happy buddha models.
Case
(a) Armadillo
(b) Stanford bunny
(c) Happy buddha
Domain size
0,2
×
0,2.4
×
0,2
0,1
×
0,1
×
0,1
0,1
×
0,2
.
2
×
0,1
Mesh grid points
50
×
60
×
50
60×60×50
100
×
220
×
100
Spatial step size
h=0.04
h=0.0167
h=0.01The columns in Figure9 represent the initial condition and numerical results after 1 and 3 iterations. As we expected, the numerical results are smooth and preserve the initial feature.Figure 9
Smoothing morphologies atT=0 (first column), T=Δt (second column), and T=3Δt (third column). Here, (a) Armadillo, (b) Stanford bunny, and (c) happy Buddha models are used.
(a)
(b)
(c)To show the effectiveness of the proposed smoothing approach, we compare results in [5] with those from our proposed method. Figure 10 illustrates CPU time in seconds against the total number of vertices of a model. As shown in Figure 10, under the similar number of vertex points, our proposed method gets smooth results and a reduction of over two orders of magnitude in CPU time.Figure 10
Total number of vertices versus CPU time.Also, in Figure11, we represent the percent of the relative error of volume VerrX which is stated in (12). Through these results, we can see that the proposed method implements the smoothing procedure almost without total volume loss.Figure 11
Temporal percent of relative error of volume with initial volume of Figure9.


## 4. Conclusion

In this paper, we have introduced a novel method using the modified CH equation for smoothing piecewise linear shapes of two- and three-dimensional objects without producing significant shrinkage. The main advantages of the proposed method are very fast and easy to implement because we use a fast Fourier transform. Various computational experiments for both curve and surface smoothing problems were presented. The results demonstrate that the proposed algorithm is effective and efficient for smoothing curves and surfaces.



---

*Source: 5971295-2017-10-25.xml*