Keeping it together: a phase field version of pathconnectedness and its implementation
Abstract.
We describe the implementation of a topological constraint in finite element simulations of phase field models which ensures pathconnectedness of preimages of intervals in the phase field variable. Two main applications of our method are presented. First, a discrete steepest decent of a phase field version of a bending energy with spontaneous curvature and additional surface area penalty is shown, which leads to disconnected surfaces without our topological constraint but connected surfaces with the constraint. The second application is the segmentation of an image into a connected component and its exterior. Numerically, our constraint is treated using a suitable geodesic distance function which is computed using Dijkstra’s algorithm.
Key words and phrases:
Willmore Energy, Phase Field Approximation, Topological Constraint, Connectedness, Dijkstra’s Algorithm2010 Mathematics Subject Classification:
49M30, 90C591. Introduction
In this article we describe how to incorporate a topological constraint into a phase field simulation for certain geometric functionals. In three dimensions, the prototypical example of such an energy is Willmore’s energy, i.e., the integral of mean curvature squared, restricted to the class of manifolds embedded into a bounded domain such that the embedding has surface area and is connected. A number of more general functionals are also admissible—the precise theoretical setting for our methods is presented in section 2.1. We furthermore consider the case of functionals controlling only the perimeter of sets contained in a bounded domain . To simplify the presentation, we mostly focus on the threedimensional case with a control of Willmore’s energy.
Our method to enforce this connectedness constraint for diffuse interfaces is based on a functional introduced in [DLW17] and given by
where are continuous functions such that
for some . For a heuristic interpretation of the functional, see section 2.2.
Despite the intimidating appearance of the functional as a double integral coupled to a geodesic distance, an efficient algorithm to treat this term can be implemented. We rely on a decomposition into connected components and a Dijkstra or fast marchingtype method to approximate the geodesic distance as well as its variation.
In our simulations, we compute a discrete time gradient flow (approximately) minimizing in each time step the functional
(1.1) 
where is an admissible phase field approximation of a geometric functional and is a constant. Note that the case of course provides no penalty for disconnectedness.
The article is structured as follows. In section 2.1, we recall precise statements regarding the sharp interface limit of functionals of the form , where is a term controlling a diffuse version of Willmore’s energy as well as the perimeter [BM10, RS06]. In the remainder of section 2, we give a heuristic explanation of our topological functional and briefly review the background material on the distance function on graphs which will be required below. In section 3, we describe an implementation of the connectedness constraint. In section 4, we present simulations with and without the topological constraint. Section 5 provides an outlook to the aforementioned application to image segmentation.
2. Preliminaries
2.1. Geometric Energies and Phasefield Connectedness
We note the following convergence results.
Theorem 2.1.
[DLW17] Let and a sequence of functionals such that for every there exists independent of such that
where
is the usual ModicaMortola approximation of the perimeter functional, and
is an approximation of Willmore’s energy due to de Giorgi. If is a Radon measure such that the diffuse area measures
converge to in the weak* sense for a sequence with , then is connected.
Remark 2.2.
According to [RS06], we have at boundaries for all sequences such that also remains bounded. Furthermore, if is uniformly bounded, then up to a subsequence we know that there exist and a Radon measure supported in such that in and weakly in the sense of Radon measures. Thus our assumption that a limit exists is not detrimental to the generality of the statement.
Remark 2.3.
We note that furthermore, if
at boundaries and the recovery sequence for is given by the usual optimal profile construction, then for any , the functional satisfies
at connected boundaries. The inequality here is obvious since . If is a connected boundary with associated optimal profile sequence , then as demonstrated in [DLW17] and the inequality follows as well. This shows that our topological functional does not change existing limits, apart from enforcing connectedness.
Our main numerical example for this result is the treatment of the functional
(2.1) 
with and in three ambient space dimensions. Like above, the normalizing constant is given by . A modification of the diffuse Willmore functional , this energy includes a spontaneous curvature which the surface would prefer to take. This is a more realistic model for many biological membranes. Theorem 2.1 clearly applies to this , but the hypotheses are also satisfied by an approximation of Helfrich’s Energy as given in [BM10] or if a volume constraint is included, see [DLW17, Chapter 5]. We do note, however, that the limit of this functional is not simply given by Willmore’s energy with spontaneous curvature and area penalty, as that is not a lowersemicontinuous energy [GB93].
2.2. The Topological Term
When , has a continuous representative, so the following notions are welldefined. We can think of the interface as the preimage of any interval , for a precise statement see [DW17, Theorem 2.20]. Thus it is possible to introduce a quantitative notion of pathconnectedness of the interface at two points through a geodesic distance function
with a weight satisfying
In particular, if is (path)connected, then for all . If, however, there are multiple connected components of with a positive spatial separation, then any curve should have uniformly positive length between and if the two points are in different connected components. We now measure the total disconnectedness of by a doubleintegral of the quantitative pathdisconnectedness of at and over the entire interface with respect to both and
where is a bump function
The dependence of on the choices of and , and therefore on and , should be kept in mind, but for notational convenience we will not make it explicit in the remainder of this article. Along the usual optimal profile recovery sequence for a connected, smooth, embedded manifold, vanishes identically. The bound on enforces a strong mode of convergence for to away from which suffices for to detect a disconnected interface in the sense that if has more than one connected component, see [DLW17, DW17]. The normalization factor is used since an interface has width proportional to .
2.3. Discretizing the Geodesic Distance
Let be a finite connected (undirected) graph with vertices and edges that are assigned weights . The distance of two vertices is defined by the length of the shortest path in the graph connecting and where the length of an edge is measured by the weight , i.e.,
In our finite element setting we consider a sequence of quasiuniform triangulations (or tetrahedralizations in 3D) of our domain with a spatial grid scale for . We furthermore assume that is a given continuous function on . Let now be the dual graph associated to a triangulation , i.e., a vertex of corresponds bijectively to a triangle and that if shares an edge (or a face in 3D) with , then the corresponding vertices in are connected by an edge. To such an edge we associate the weight , where . This gives rise to a discrete distance depending on a phase field function .
The triangulations may force a minimal path to zigzag to connect two points, so the distance on the graph may not approximate the distance function
as , but, due to the nondegeneracy assumptions on the triangulation, we note that the discrete distance is equivalent to in an almost biLipschitz sense uniformly in , i.e.,
for all and such that with suitable constants . Such a modification clearly does not change the effect of the topological term on connectedness as only a nonzero lower bound and a vanishing upper bound are required for the inequality and construction, respectively.
As a treatment of the timestep minimization problem will require a variation of the distance as well, so we note that if there exists a unique shortest curve between and then
This identity will be postulated for the procedure below and we approximate the variation of the geodesic distance on the graph as we vary in direction by the discrete term
(2.2)  
where is a shortest connecting path between triangles and containing and , respectively. The above variation is in the spirit of [BCPS10, BLS15] where it was derived for the fast marching method.
3. The Algorithm
In this section, we describe how to include the topological term in an explicit fashion in given finite element code. We note that the limiting effect on the timestep size of this explicit treatment is moderate, as the term does not depend on spatial gradients of the phase field function .
The description is given in the twodimensional case assuming that the finite element space corresponds to a triangulation of with grid length scale . Dimension three and more general basis element shapes can be treated by the same method.
In the set up of the simulation, we create the dual graph corresponding to the finite element triangulation as described in section 2.3 (with the edge weights left unassigned for the time being, as they will change in each time step). Each triangle is furthermore associated with its volume and diameter .
Given a Galerkin space function in time step , do the following.

For all triangles in the triangulation, compute the average integral

For each edge in corresponding to two adjacent triangles and , compute its weight as
The second factor can be replaced by a generic grid length scale if the triangles are sufficiently uniform.

Create a list of all interface elements, i.e. all triangles such that

Separate the elements in into connected components , where two triangles belong to the same component if . If there is only one connected component, the algorithm can be terminated here as our approximations of both and its variation vanish.

For calculate

For , calculate the component distances as well as the shortest connecting paths between components , , . These computations can be performed using Dijkstra’s algorithm on , see [Dij59].

The approximate topological energy can now be computed as
We note that, compared to the original double integral term, this is a major simplification which is due to the specific choice of vanishing identically on the support of .

Our discrete approximation of the variation of with respect to a finite element basis function is then given by
where for any , is given by (2.2).
This algorithm can be added to a given finite element implementation. We may compute the timestep from to with any scheme
that treats the topological term explicitly. Here is the timestep size and is an explicit, implicit or mixed approximation of the variation . In simulations, it has proven useful to use the sum of two functionals of type : one, to keep the portion of the interface close to connected (i.e., with and close to ) and one to keep the portion of the interface close to connected (i.e., with and close to ).
4. Numerical Results
We compare the discrete gradient flow (1.1) of , where is given in (2.1) and either or . The computational domain is a discretization of the three dimensional unit ball by approximately tetrahedral finite elements. Our timestepping algorithm is a simple first order fully implicit Euler scheme coupled to an explicit treatment of as described in section 3. The time step size is given by (with a smaller time step if for the first several hundred time steps during the fast convergence to the optimal profile). The higher spatial gradient in the energy is treated using a CiarletRaviartMonk mixed formulation [CR74, Mon87] with clamped boundary conditions on and on and we note that on our convex domain this variational crime is only a misdemeanor [GSS12].
The parameters in the simulation are , , and . The initial condition is an approximation of the characteristic function of an ellipsoid with principle axes , , and . As seen in Figure 1, in the simulation of the case , without topological constraint, the surface undergoes a pinchoff and the final steady state is given by two spheres of radius approximately equal . A similar result of pinchoff was observed in [DLW06]^{*}^{*}*The simulation in [DLW06] was performed for . We note that our simulation produces virtually the same results for , but we present the case of positive since none of the analytic results in [RS06] are admissible unless remains uniformly bounded as ..
The simulation results for the case , i.e., including the topological constraint, are shown in Figure 2. One can clearly see that the pinchoff into two components has been suppressed and a dumbbelllike shape is the final result. We conjecture that this shape is in fact a stationary point of the energy as no further motion was observed in the simulation even on a longer timescale.
As mentioned at the end of section 3, we implement the connectedness constraint using the sum of two functionals of type , in this specific case one with , (in order to keep the part of the transition layer close to the phase connected) and another one with , (in order to keep the part of the transition layer close to the phase connected). The functions and in are given by
(4.1)  
(4.2) 
respectively, with and chosen such that and such that . We note that in our case, at the pinchoff, the function dips below and thus the start of the interface becoming disconnected is detected by the algorithm. For an illustration see Figure 3.
The equilibrium was reached in approximately 9 hours of walltime using 8 cores of a computer server equipped with two Intel Xeon E52690 v4 processors. We note that the cputime spent computing the topological constraint is negligible (less than 0.1% of the total cputime)—the only computational downside may thus be the aforementioned restriction on the timestep size due to the necessary explicit treatment of . Implicit treatment of is analytically questionable due to the lower regularity of .
5. An application to image segmentation
To conclude, we give an outlook to a different application of the topological functional introduced in this article. We consider the energy
with and given and . Again, the case is that of no disconnectedness penalty, and for disconnectedness of the set is penalized.
For and using a usual doublewell energy (in this case with minima at and and normalization constant ), the functional is a typical imagesegmentation functional with perimeter regularization and a fidelity term.
For we now choose and and ^{†}^{†}†The functions and , as well as the constants and are picked as before, the constant is irrelevant, noting that in the simulation it always holds that . The functional, is given by , since the integrals are now performed over open sets and not over boundary layers.: this means that the topological term will enforce pathconnectedness of the set and so the grey scale image given by should be segmented into a connected segment and its exterior.
In the present example, is the unit square, and . The discretization is given by approximately triangular finite elements. We again compute a timediscrete gradient flow (using now a semiimplicit first order Euler scheme with only the linear highest gradient term being treated implicitly) of the energy until a stationary state has been reached. In all cases, the initial condition is given by the characteristic function of the set in polar coordinates.
In Figures 4 and 5 the results are shown. On the left, the source image to be segmented is displayed (in our simple example, itself only takes values in and is the characteristic function of two separated disks). The two other pictures show the computed phase field minimizer , first without connectedness penalty () then with connectedness penalty (). One can clearly see that without this penalty, the source image is simply reproduced.
large disks  small disks  

Radii  0.16  0.11 
Distance between centers of disks  0.6  0.6 
Twice the distance between disks  0.56  0.76 
Fidelity penalty for removing one disk  0.84  0.40 
In Figure 4 connectedness is established by adding a doublelayer between the two disks in the image . With the disks being smaller in Figure 5 and our choice of parameters, the double layer would be more energetically costly than the error being made in the fidelity term, so one of the disks is simply being ignored. The respective values for the energies have been assembled in Table 1. We do note, however, that a global optimum cannot generally be found by such a gradient flow as local minimal of the energy exist (in our setting, the global minimizers are given by everywhere). A rigorous analysis of this imagesegmentation energy with connectedness constraint will be presented in [DNWW].
Acknowledgements
PWD gratefully acknowledges partial support by the WissenschaftlerRückkehrprogramm GSO/CZS as well as inspiring discussions with Benedikt Wirth (Münster), Sebastian Reuther (Dresden), and Douglas N. Arnold (Minneapolis).
References
 [BCPS10] F. Benmansour, G. Carlier, G. Peyre, and F. Santambrogio. Derivatives with respect to metrics and applications: subgradient marching algorithm. Numerische Mathematik, 116(3):357–381, 2010.
 [BLS15] M. Bonnivard, A. Lemenant, and F. Santambrogio. Approximation of length minimization problems among compact connected sets. SIAM J. Math. Anal., 47(2):1489–1529, 2015.
 [BM10] G. Bellettini and L. Mugnai. Approximation of the Helfrich’s functional via diffuse interfaces. SIAM J. Math. Anal., 42(6):2402–2433, 2010.
 [CR74] P. G. Ciarlet and P.A. Raviart. A mixed finite element method for the biharmonic equation. In Mathematical aspects of finite elements in partial differential equations, pages 125–145. Elsevier, 1974.
 [Dij59] E. W. Dijkstra. A note on two problems in connexion with graphs. Numerische Mathematik, 1(1):269–271, 1959.
 [DLW06] Q. Du, C. Liu, and X. Wang. Simulating the deformation of vesicle membranes under elastic bending energy in three dimensions. Journal of Computational Physics, 212(2):757–777, 2006.
 [DLW17] P. W. Dondl, A. Lemenant, and S. Wojtowytsch. Phase Field Models for Thin Elastic Structures with Topological Constraint. Arch. Ration. Mech. Anal., 223(2):693–736, 2017.
 [DNWW] P. W. Dondl, M. Novaga, B. Wirth, and S. Wojtowytsch. A phasefield approach to connected perimeters in the plane. In preparation.
 [DW17] P. W. Dondl and S. Wojtowytsch. Uniform convergence of phasefields for Willmore’s energy. Calc. Var. PDE, 56(4), 2017.
 [GB93] K. GroßeBrauckmann. New surfaces of constant mean curvature. Mathematische Zeitschrift, 214(1):527–565, 1993.
 [GSS12] T. Gerasimov, A. Stylianou, and G. Sweers. Corners Give Problems When Decoupling Fourth Order Equations Into Second Order Systems. SIAM Journal on Numerical Analysis, 50(3):1604–1623, January 2012.
 [Mon87] P. Monk. A mixed finite element method for the biharmonic equation. SIAM Journal on Numerical Analysis, 24(4):737–749, 1987.
 [RS06] M. Röger and R. Schätzle. On a modified conjecture of De Giorgi. Math. Z., 254(4):675–714, 2006.