J.-P Pons

J.-D Boissonnat

A Lagrangian Approach to Dynamic Interfaces through Kinetic Triangulation of the Ambient Space

Keywords: curve, interface, deformation, Lagrangian, topology, collision, triangulation, kinetic data structure, exact computation, multiple junction ACM CCS: I.3.5 Computer Graphics: Computational Geometry and Object Modeling -Curve, surface, solid, and object representations, 6.8 Simulation and Modeling Types of simulation -Animation

In this paper, we propose a robust and efficient Lagrangian approach for modeling dynamic interfaces between different materials undergoing large deformations and topology changes, in two dimensions. Our work brings an interesting alternative to popular techniques such as the level set method and the particle level set method, for two-dimensional and axisymmetric simulations. The principle of our approach is to maintain a two-dimensional triangulation which embeds the one-dimensional polygonal description of the interfaces. Topology changes can then be detected as inversions of the faces of this triangulation. Each triangular face is labeled with the type of material it contains. The connectivity of the triangulation and the labels of the faces are updated consistently during deformation, within a neat framework developed in computational geometry: kinetic data structures. Thanks to the exact computation paradigm, the reliability of our algorithm, even in difficult situations such as shocks and topology changes, can be certified. We demonstrate the applicability and the efficiency of our approach with a series of numerical experiments in two dimensions. Finally, we discuss the feasibility of an extension to three dimensions.

Introduction

Modeling dynamic interfaces between several materials undergoing large deformations is a ubiquitous task in science and engineering. The existing techniques roughly fall into two categories: Eulerian and Lagrangian formulations. It is commonly admitted that both viewpoints have strengths and weaknesses and that only a hybrid approach can overcome the limitations of both. In this paper, we propose a method which has the particularity of being purely Lagrangian and of retaining all the associated advantages, while achieving topological adaptivity with a comparable robustness and efficiency to Eulerian methods.

Eulerian methods

The Eulerian formulation casts deformation as a time variation of quantities defined over a fixed grid. The inter-faces have to be represented implicitly, since the grid does not conform to them. In computational physics, this is also known as the front capturing method. Two notable front capturing techniques are the level set method, introduced by Osher and Sethian [START_REF] Osher | Fronts propagating with curvature-dependent speed: Algorithms based on Hamilton-Jacobi formulations[END_REF], and the volumeof-fluid (VOF) method, pioneered by Hirt and Nichols [START_REF] Hirt | Volume of fluid ({VOF}) method for the dynamics of free boundaries[END_REF].

Hereafter, we will mainly focus on the level set method because it is an established technique in fluid animation. Basically, this method consists in representing the interface as the zero level set of a higher-dimensional scalar function. The movement of the interface can be cast as an evolution of the embedding level set function by an Eulerian PDE (partial differential equation). We refer the reader to some good reviews [START_REF] Sethian | Level Set Methods and Fast Marching Methods: Evolving Interfaces in Computational Geometry[END_REF][START_REF] Osher | Level set methods: an overview and some recent results[END_REF] for all the details about the theory, the recent developments, the implementation and the applications of the level set method.

On the one hand, this approach has several advantages over an explicit Lagrangian representation of the interface: no parameterization is needed, topology changes are handled automatically, intrinsic geometric properties such as normal or curvature can be computed easily from the level set function. Last but not least, the theory of viscosity solutions provides robust numerical schemes and strong mathematical results to deal with the evolution PDE. These advantages explain the popularity of the level set method for multi-phase fluid flow simulation [START_REF] Sussman | A level set approach for computing solutions to incompressible two-phase flow[END_REF] in CFD (computational fluid dynamics), as well as for computer animation of fluids with free surfaces [FF01; EFFM02; EMF02; LGF04; ELF05].

On the other hand, several serious shortcomings limit the applicability of the level set method. First, the higher dimensional embedding makes the level set method much more expensive computationally than explicit representations. Much effort has been done to alleviate this drawback, leading to the narrow band methodology [START_REF] Adalsteinsson | A fast level set method for propagating interfaces[END_REF] and to the PDE-based fast local level set method [PMO * 99]. More recently, octree decompositions have been proposed [LGF04; ELF05; LFO06; COQ06] to circumvent the typically fixed uniform sampling of the level set method, in order to reach high resolution (typically an effective resolution of 512 3 ) while keeping the computational and memory cost sustainable. However, these methods somewhat lose the simplicity of the original level set method, as an efficient implementation of such tree-based methods turns out to be a tricky task.

Second, as discussed and numerically demonstrated by Enright et al. in [START_REF] Enright | A hybrid particle level set method for improved interface capturing[END_REF], the level set method is strongly affected by mass loss, smearing of high curvature regions and inability to resolve very thin parts. These limitations have motivated the development of some hybrid Eulerian-Lagrangian methods, such as the particle level set method outlined by Foster and Fedkiw [START_REF] Foster | Practical animation of liquids[END_REF] and later improved by Enright and coworkers [EFFM02; EMF02; ELF05]. While the latter method yields state-of-the-art results, an objection could be the large number of parameters controlling the particle reseeding strategy included in this approach. To some extent, also related is the contouring method [Str01; SJ02; BGOS06], which converts back and forth from a Lagrangian mesh to an implicit function evaluated on a regular Eulerian grid using distance computations and isocontour extraction algorithms.

Third, purely Eulerian formulation is not very appropriate for tracking interface properties such as color or texture coordinates, as may be needed in computer graphics applications. An approach based on a coupled system of Eulerian PDEs was recently proposed by Pons et al. to overcome this limitation within the level set framework [START_REF] Pons | Maintaining the point correspondence in the level set framework[END_REF], but this capability comes at a significant additional computational cost.

Previous Lagrangian methods

The Lagrangian formulation adopts a more "natural" point of view. It explicitly tracks the interfaces between the different materials with some points advected by the motion. There are mainly two classes of Lagrangian techniques: meshbased methods and particle-based methods. We first tackle the mesh-based approach, also known as the front tracking method in computational physics [TBE * 01, UT92]. When dealing with large deformations, this approach is hampered by distortion and entanglement of the mesh, source of numerical instabilities or even of breakdowns of the simulation. For instance, if corners and cusps develop in the evolving front, front tracking methods usually form "swallowtail" solutions. These defective parts of the interface must be detected, then removed through intricate delooping procedures.

Another major shortcoming of the mesh-based Lagrangian approach is that a fully automatic, robust and efficient handling of topology changes remains an open issue, despite many heuristic solutions proposed in various fields of application (e.g. [MT99, MT00; LM99] in medical imaging, [TBE * 01, SJ02] in computational fluid dynamics, . . .).

McInerney and Terzopoulos [MT99, MT00] propose topology adaptative deformable curves and meshes, called

T-snakes and T-surfaces. During the evolution, the model is periodically resampled by computing its intersections with a regular simplicial decomposition of space. A labeling of the vertices of the simplicial grid as inside or outside of the model is maintained. Besides being computationally expensive, this procedure loses the desirable adaptivity of the Lagrangian formulation, by imposing a fixed uniform spatial resolution.

Lachaud and Montanvert [START_REF] Lachaud | Deformable model with a complexity independent from image resolution[END_REF][START_REF] Lachaud | Deformable meshes with automated topology changes for coarse-tofine 3D surface extraction[END_REF] use the concept of δ-triangulation. A length parameter δ is used to control the sampling of the mesh and to detect self-intersections, by monitoring the distance between pairs of neighbor and non-neighbor vertices. Besides being only approximate, this approach for topological flexibility is costly, even when optimizing the pairwise distance computations with an octree structure.

This major shortcoming of mesh-based methods have gained popularity to the particle-based approach [TPF89; DG95; FM96; DC96; CD97; MCG03; PTB * 03; MKN * 04; PKA * 05] for representing dynamic interfaces undergoing complex topology changes. However, it is generally admitted that with the meshless approach, the localization of the interface and the computation of interface properties such as normal and curvature gets cumbersome.

For sake of completeness, let us also mention a recent approach for fluid simulation based on dynamic meshes [START_REF] Klingner | Fluid animation with dynamic meshes[END_REF]. This method bears some resemblance with ours, in the sense that it maintains a triangulation of the computational domain during motion, which conforms to a moving object. However, their purpose, fluid/rigid body coupling, is different from ours, fluids with free interfaces: in their case, the geometry and the topology of the front is fixed.

Novelty of our method

In this paper, we present a purely Lagrangian approach, which combines the advantages of front tracking and front capturing methods, while discarding their respective drawbacks. Our work brings a robust and efficient solution to remeshing and topological adaptivity for Lagrangian dynamic interfaces in two dimensions. The principle of our approach is to maintain a two-dimensional triangulation of space which embeds the one-dimensional polygonal description of the interfaces. Interestingly, this increase of dimension bears similarity to the spirit of the level set method. Topology changes can then be detected as inversions of the triangular faces (called faces for short) of the triangulation. Each face is labeled with the type of material it contains. This embedding triangulation enforces directly watertight interfaces free of loops, swallowtails or self intersections at all times.

The connectivity of the triangulation and the labels of the faces are updated consistently during deformation, within the kinetic data structures framework, introduced by Basch et al. [START_REF] Basch | Data structures for mobile data[END_REF]. Rather than trying to repair the possibly entangled triangulation after each iteration of the simulation, we smoothly interpolate the coordinates of interface points between two consecutive time steps, and we modify the connectivity of the triangulation and the labels of the faces exactly as and when it is required. As a result, our lazy kinetic triangulation has the desirable property not to introduce any unnecessary perturbation of the interfaces in the absence of topology changes, and to handle them exactly and efficiently when they occur.

Being purely Lagrangian, our method does not suffer from mass loss which plagues the level set method, and can track material properties such as color or texture coordinates during motion at no additional cost, while being free of the localization problem of particle-based approaches.

Moreover, our method is not limited to two-phase motion. It seamlessly accommodates any number of materials, whereas this requires special care in most existing Lagrangian and Eulerian methods (cf for example [START_REF] Smith | A projection method for motion of triple junctions by level sets[END_REF] on the treatment of triple junctions with the level set method).

The remainder of this paper is organized as follows. Section 2 gives some background on kinetic data structures, and particularly on kinetic triangulations and their application to efficient collision detection. The different elements of our lazy kinetic triangulation are described in Section 3. In Section 4, we give the details of our implementation with CGAL (Computational Geometry Algorithms Library) [START_REF] Boissonnat | Triangulations in CGAL[END_REF] and we demonstrate the applicability and the efficiency of our approach with a series of numerical experiments. Finally, we discuss an extension to three dimensions in Section 5.

Background

Kinetic data structures

Kinetic data structures, introduced by Basch et al. [START_REF] Basch | Data structures for mobile data[END_REF], are a general framework to efficiently maintain a discrete attribute of a set of moving objects in time. Typically, this attribute consists of a combinatorial description such as the convex hull [START_REF] Basch | Data structures for mobile data[END_REF] or the Delaunay triangulation [AGMR98, GR04] of a set of moving points. This framework was motivated by the important limitations of the method of widest practical use, which is a naive incremental update of the attribute at some predefined (usually uniformly distributed) time samples. Because the distribution in time of the actual changes of the attribute is far from uniform, the choice of the time step is problematic. Any time step is bound to oversample the system sometimes, then performing useless computations, and to undersample it sometimes, then missing some possibly important intermediate states of the attribute.

In contrast, kinetic data structures take advantage of the knowledge of the motion to exactly determine the times when events occur, i.e. when the attribute changes. No computation is made between two consecutive events. This is achieved by maintaining a set of elementary geometric predicates, which is logically equivalent to the current state of the attribute. For instance, the validity of a 2D kinetic Delaunay triangulation [AGMR98, GR04] is certified by one InCircle predicate for each triangle, plus some Orientation predicates for the edges of the convex hull. The next change of the attribute is then scheduled at the earliest time for which one of the predicates fails. This is possible in practice whenever the trajectories of the objects are described by low-degree polynomials, so that the failure times of the predicates are roots of polynomials. Moreover, the efficient processing of the events in chronological order is achieved by storing the failure times in a global priority queue. When an event occurs, the attribute is properly modified, then some predicates affected by the change are created, rescheduled or deleted, while the priority queue is updated accordingly.

We refer the reader to some recent reviews [AGE * 02, Gui04] for a thorough description of the kinetic data structures framework and its applications to the maintenance of classical combinatorial structures such as convex hulls, Voronoi diagrams and Delaunay triangulations, closest pairs, minimum spanning trees, etc.

Kinetic triangulations

Maintaining a triangulation of a set of moving points is of particular interest for our problematic of modeling materials undergoing large deformations and topology changes. There has been a significant amount of work on time-dependent triangulations within the kinetic data structures framework, but most of it was focused on algorithmic complexity.

An important result in this area was obtained by Agarwal et al. [ABdB * 99]: maintaining a triangulation with a kinetic data structure requires at least a quadratic number of events in the worst case, even if the addition of a linear number of Steiner (i.e. additional) points is allowed. So far, no algorithm was shown to attain this lower-bound result. The algorithm described in [ABG * 02] processes O(N 7/3 ) events in the worst case. More recently, in [START_REF] Agarwal | A 2D kinetic triangulation with near-quadratic topological changes[END_REF], an algorithm was proposed that only requires N 2 2 O( √ log N log log N ) events.

In the particular case of the Delaunay triangulation, the gap between theoretical bounds and existing algorithms is wider, with a quadratic best lower bound (i.e. the same as for an arbitrary triangulation) and a nearly cubic best known upper bound [START_REF] Albers | Voronoi diagrams of moving points[END_REF].

The spirit of our work significantly differs from that of the above references. We have several reasons to believe that the quadratic lower bound reported in [ABdB * 99] is not relevant to our problem. First, our context is far from the worstcase scenario, since the boundary of the materials where the points are distributed, as well as the motion, typically exhibit some regularity properties. In [START_REF] Attali | Complexity of the Delaunay triangulation of points on surfaces: the smooth case[END_REF], it was shown that such a smoothness assumption reduces the complexity of the static 3D Delaunay triangulation of N points from O(N 2 ) to O(N log N ). We expect a similar complexity drop for a kinetic triangulation too, although we are not aware of a theoretical result in this direction.

Second, we do not resort to a canonically defined triangulation, such as a Delaunay triangulation [AGMR98, GR04] or a constrained fan triangulation [START_REF] Agarwal | A 2D kinetic triangulation with near-quadratic topological changes[END_REF]. In other words, the state of our kinetic triangulation does not depend only on the positions of the points, but also on the history of the motion. As a consequence, it is very difficult to obtain formal statements on the complexity of our approach.

Last, our kinetic triangulation has the particularity of being coupled with a numerical simulation which computes the positions of the points at evenly distributed time instants. At each time step, all the trajectories are updated and the whole set of predicates is recomputed. With reasonably small time steps, few events occur during an interval of the simulation. Under these conditions, the computational expense of our kinetic data structure is mostly related to its compactness, as defined in [AGE * 02, Gui04], i.e. the number of predicates used to certify the triangulation. It is exactly N -2, the number of triangles, with our simple algorithm. In the following, we will discuss possible improvements on this number.

Collision detection

An efficient and exact detection of collisions between several moving polygonal objects can be achieved by maintaining a tiling of the free space between them with a kinetic data structure. A triangulation is the more immediate candidate for such a tiling. However, several authors advocate the use of a more complex structure, namely a pseudo-triangulation [BEG * 99; KS02; ABG * 02], that is to say a tiling of space with pseudo-triangles. A pseudo-triangle is a simple polygon composed of three concave chains joining at their endpoints. A kinetic pseudo-triangulation is much harder to implement than a kinetic triangulation, but has significant advantages in some situations.

Compared to a kinetic triangulation, a kinetic pseudotriangulation requires a smaller number of elements and adapts to the free space between moving objects with fewer events. Also, it is very compact, as defined previously, under some assumptions on the type of motion, on the separation between the objects or on their convexity. However, these methods do not apply to non-convex objects under non-rigid motion, with the exception of [ABG * 02]. Unfortunately, the latter approach does not guarantee a better compactness than our simple kinetic triangulation in the case of nonrigid motion, which makes the choice of a kinetic pseudotriangulation questionable for our application. Also, the reader should note that the problem that we address in this paper is less restrictive than collision detection: our method allows the simulation to continue after a collision, by properly handling the fusion or the splitting of the different materials.

Methods

The principle of our approach is to couple the numerical solver which outputs the position of the vertices at some (typically evenly-spaced) time instants, with a kinetic triangulation. The numerical solver is specific to each application and is not discussed here. Some examples will be briefly presented in Section 4. In this section, we describe the particularity of our lazy kinetic triangulation and explain its pertinence for the modeling of dynamic interfaces.

Embedding triangulation

A determinant feature of the triangulation we maintain is that a physical interpretation is attached not only to its vertices, but also to its faces, namely the type of material that it contains. For instance, in a multi-phase fluid simulation, each triangle is mapped to one of the phases. This information, that we call for short the label of a triangle, can be conveniently represented by an integer value. The interfaces of interest are embedded in this triangulation. They are composed of the edges adjacent to two triangles having different labels. This representation has many virtues. It automatically enforces watertight interfaces. It seamlessly accommodates any number of materials. It makes the request of the type of material at some point very efficient. An example of this embedding triangulation, including three materials with some triple junctions, is displayed in Figure 1.

Lazy approach

We now describe how the connectivity of the triangulation and the labels of the faces are consistently updated during motion, in order to prevent the breakdown of the simulation, and to faithfully reflect the deformations and the topology changes of the different materials.

Our kinetic triangulation is termed "lazy" because it undergoes no connectivity change as long as it remains a valid geometric triangulation, that is to say as long as the orientations of the triangles do not change and the triangulation covers the convex hull of the vertices. In other words, the triangulation is modified only when a triangle becomes flat or two consecutive border edges become collinear. The corresponding set of predicates is an Orientation predicate for each triangle, plus an Orientation predicate for each triplet of consecutive border vertices. The events are processed with elementary local connectivity modifications, namely edge flips and edge collapses, that we will describe in detail later.

Contrarily to a Delaunay triangulation, the state of our triangulation is not canonically defined: it does not depend only on the positions of the points, but also on the history of the motion. As a result, it does not suffer from an instability and a multiplication of meaningless events when the points are very close to a configuration change, like it happens with Delaunay when four or more points are nearly co-circular.

Another advantage of an arbitrary triangulation is that the predicates are cheaper to compute than for Delaunay. For polynomial trajectories of degree d, the failure times of the Orientation predicates are roots of polynomials of degree 2d, against degree 4d for the InCircle predicates.

Note that the lazy strategy may not be efficient on the long term, as the total number of events may be higher than with a more foreseeing strategy like Delaunay. However, the actual motivation for this choice is not efficiency, but fidelity to deformations and topology changes of the interfaces. The lazy behavior guarantees that the triangulation remains a faithful embedding of the polygonal interfaces. A connectivity change for the convenience of the kinetic data structure would constitute an unnecessary perturbation of the location of the interfaces.

An objection to the lazy kinetic triangulation could be the absence of control on the quality of the triangles. In some other contexts where the numerical robustness directly depends on the quality of the worst-shaped element of the triangulation, this can impair the simulation. This limitation does not apply to our approach, because we are only interested in the interfaces between the different materials. The tiling of each material is not used in the simulation, it is only available internally to the kinetic data structure.

Predicates and motion model

The convex hull can be kinetically maintained by considering an orientation predicate for each triplet of consecutive border vertices. However, for the sake of simplicity, we assume here that the convex hull remains constant during motion. This can be easily achieved by adding to the triangulation fixed points that delimit a bounding box of the simulation. As a result, the certification of our lazy kinetic triangulation only requires one orientation predicate for each triangle.

As for the motion model of the vertices between two consecutive steps of the simulation, it can be chosen arbitrarily with a canonically-defined structure such as the Delaunay triangulation, as long as it respects the continuity of the trajectories. In this case, indeed, the intermediates states of the triangulation can be safely disregarded. In [START_REF] Guibas | An empirical comparison of techniques for updating Delaunay triangulations[END_REF], Guibas and Russel empirically compare the efficiency of three different motion models for the 3D kinetic Delaunay triangulation: a linear motion of all the points simultaneously, a linear motion of one point at a time, and a linear motion of all the points along one coordinate at a time.

In contrast, the state of our lazy kinetic triangulation depends on the history of the motion, and the choice of the motion model conditions the fidelity to the simulated physical phenomenon. Thus, the "one point at a time" and "one coordinate at a time" motion models would be questionable, despite their higher efficiency. In most applications, linear or quadratic trajectories are relevant. Quadratic trajectories faithfully model accelerated motion, but are significantly more expensive computationally, because they require to manipulate roots of polynomials of degree 4.

Processing of the events

First, we consider the generic event of a flat triangle whose vertices have three distinct positions. We use the classical Lawson edge flip [START_REF] Lawson | Software for C 1 surface interpolation[END_REF] to restore a valid geometric triangulation. The longest edge of the flat triangle is flipped, as shown in Figure 2. An edge flip requires the two adjacent triangles to form a convex quad. Note that this condition is automatically fulfilled in our case, since three vertices are collinear.

Also, consistent labels must be assigned to the two new triangles generated by the edge flip: these two triangles are given the label of the original non-flat triangle involved in the flip (cf Figure 2), so that the location of the interfaces is not altered. Finally, the predicates of these two new triangles are built and their failure times are placed in the priority queue.

As a result of these changes in the connectivity and in the labels, some vertices of the triangulation may no longer participate in the interface, i.e., all the adjacent triangles may now have the same label. It typically happens in the case of a topology change of the materials or of a shock, as we will show in our experiments in Section 4. For the sake of efficiency, these vertices are removed from the triangulation and their holes are retriangulated, again without affecting the position of the interfaces. There is no ambiguity in assigning labels to the newly created triangles.

We now describe how to cope with degenerate cases, in order to make our algorithm perfectly reliable. Flat triangles with two or three vertices having the same position are handled with the edge collapse operation. It consists in collapsing the two vertices of an edge into a single vertex. One of the vertices of the edge, the edge itself and the two triangles adjacent to the edge are deleted. In the edges and triangles involving the deleted vertex, the latter is replaced by the other vertex of the collapsed edge. In general, an edge collapse may generate triangles of reverse orientation, and make a geometric triangulation invalid. Note that it is not possible in our case, since the vertices of the collapsed edge have the same position. As regards the labels of the triangles, no change is needed after an edge collapse. The handling of simultaneous events also requires special care. It happens when several triangles become flat exactly at the same time. The kinetic data structures framework does not prescribe the order in which simultaneous events must be processed. This order depends on the particular attribute maintained, and turns out to be important for our lazy kinetic triangulation. Indeed, processing the events in a wrong order can compromise the termination of the algorithm. This is due to the fact that, in some degenerate cases, the processing of an event can generate one of several new flat triangles. For example, it happens when an edge between two flat triangles is flipped. We give an example of such a situation in Figure 3. There are two possible edge flips for these two adjacent flat triangles. If the shortest of the two possible flipping edges is flipped first, two new flat triangles are created, and the final configuration is qualitatively the same as before the flip. If the algorithm constantly chooses the wrong flipping edge, it enters an infinite loop.

In contrast, in Figure 4 we show a configuration in which an edge flip between two flat triangles makes some progress. The difference between these two degenerate configurations inspires the following strategy to handle simultaneous events. We perform in priority any needed edge collapse. Then, we process any needed edge flip, except that we skip the flipping edges adjacent to two flat triangles if the flipping edge is not the largest edge of both triangles. After an edge flip, new edge collapses may become possible: we perform them before considering the next edge flip. Now, we give an elementary proof, in three arguments, that the algorithm terminates with this strategy, and hence that it restores a valid geometric triangulation. First, it is trivial that the sequence of edge collapses terminates, since each edge collapse deletes a vertex. Second, among the possible edge flips, there is always an acceptable one, for instance the one with the largest flipping edge. Third, for each acceptable edge flip: r either it strictly decreases the number of flat triangles and of flipping edges, as shown in Figure 2, r or it strictly decreases the length of the possible flipping edges, while keeping the number of flat triangles constant, as illustrated in Figure 4.

Throughout this section, we have assumed that the convex hull remains constant. However, handling specifically the events involving border vertices and border edges suffices to extend the above algorithm to the case of a varying convex hull.

Resampling of the interfaces

In most applications, the resolution of the interfaces has to be corrected during deformation, to maintain the desired localization accuracy. Although this task is highly applicationdependent, we propose a basic technique to control the distribution of the lengths of interface edges. This resampling step is performed at the end of each step of the simulation, after a valid geometric triangulation has been restored by the kinetic data structure.

Basically, the edges of the interfaces, hence not all the edges of the embedding triangulation, are tested against a maximum and a minimum threshold. Long edges are bisected by adding their center point to the triangulation, then replacing each of the triangle adjacent to the edge by two new triangles with identical labels. Note that this modification neither changes the location of the interface, nor jeopardizes the validity of the triangulation.

In contrast, collapsing short edges requires special care, because it alters the position of the interface and it can sometimes generate inverted triangles. We emphasize that, although very short edges are an unnecessary computational burden, they do not affect the stability of the simulation, as our approach is not subject to the formation of swallowtails.

Numerical Experiments

Implementation issues

By using CGAL (Computational Geometry Algorithms Library, homepage: www.cgal.org) [BDTY00], we have been able to implement our approach with only 1500 lines of C++ code. We have taken advantage of the flexibility of this library to derive our labeled triangulation data structure from the existing triangulation data structure. In particular, we have overloaded the edge flip and edge collapse operations to consistently handle the labels of the triangles.

We have developed our own implementation of the kinetic data structure. Our priority queue is implemented with a pairing heap [START_REF] Fredman | The pairing heap: A new form of self-adjusting heap[END_REF], an efficient self-adjusting heap data structure. In the future, we will probably consider conforming to the CGAL package implementing the computational framework of Guibas, Karavelas and Russel [START_REF] Guibas | A computational framework for handling motion[END_REF].

In order to make our implementation robust, the exact computation paradigm [Yap97] can be followed. The coordinates of the vertices may be represented with arbitrary precision rational numbers provided by the GMP (GNU Multiple Precision) arithmetic library [GMP]. In this case, the failure times of the predicates are the roots of polynomials with rational coefficients. The method of Emiris and Tsigaridas [START_REF] Emiris | Comparing real algebraic numbers of small degree[END_REF], based on Sturm sequences, may be used to sort them in the priority queue. The latter method also allows to exactly determine the sign of a polynomial at the root of another polynomial. This can be used to robustly compare the lengths of the different edges of a flat triangle, to choose what edge collapse and what edge flip must be performed.

However, in general, the exact computation paradigm does not apply to the numerical solver. Consequently, in our implementation, for efficiency purposes, we have used a fixedprecision floating-point arithmetic both for the solver and for the lazy kinetic triangulation. This avoids converting the coordinates of the vertices back and forth between floats and exact number representations. We have not encountered any robustness issues in our numerical experiments.

Mass loss

In our first numerical experiment, we compare the behavior of the level set method and of our method regarding mass loss, using an experimental methodology initially proposed by LeVeque [START_REF] Leveque | High-resolution conservative algorithms for advection in incompressible flow[END_REF], and later popularized in computer graphics by Enright and coworkers [EFFM02; EMF02; ELF05].

The spatial domain is the unit square [0, 1] × [0, 1], the initial interface is a circle with radius 0.15 centered at (0.5, 0.75) and the evolution is driven by an analytical incompressible velocity field with non-constant vorticity that varies sinusoidally in time with a period T. Specifically, the velocity field v is defined by

v(x, y, t) = cos πt T sin 2 (π x) sin(2πy) -sin(2π x) sin 2 (πy) . ( 1
)
This experiment is challenging because the flow considerably stretches the interface, as can be seen in Figure 5. Note that the exact solution of this evolution is not available at all time, so at first sight we cannot measure the error. But the velocity reverses at time T /2, so that the initial interface should be recovered at time T. This provides a convenient way to evaluate the mass loss caused by each method. In this experiment, we take T = 4.

Due to the gap between Lagrangian and Eulerian viewpoints, a direct comparison between the two methods for the same parameter values is not possible. We have done our best in adjusting these parameters so as to obtain the fairest possible comparison. In both experiments, we use a simple Euler time stepping scheme with a time step dt = 0.05/128. We use a highly-optimized home-made implementation of the level set method used in our previous work [START_REF] Pons | Maintaining the point correspondence in the level set framework[END_REF], with a 128 × 128 resolution. To fix ideas, extracting the isocontour from the initial level set function generates 156 vertices. We set to 100 the initial number of vertices in our kinetic triangulation.

t = 0 t = 0.5 t = 1 t = 1.5 t = 2 t = 2.5 t = 3 t = 3.5 t = 4
Figure 5 displays the evolution obtained with the level set method (top) and with our method (bottom) at different times. As expected, the level set method fails to preserve the thin parts of the interface, despite a quite high resolution and the use of a high-order spatial differencing scheme (namely a fifth-order WENO scheme [START_REF] Jiang | Efficient implementation of weighted ENO schemes[END_REF][START_REF] Jiang | Weighted ENO schemes for Hamilton-Jacobi equations[END_REF]). As a result, a significantly higher mass loss is observed with the level set method.

Example of coalescence

In our second experiment, we demonstrate how our approach handles topology changes, in the simple case of the coalescence of two expanding objects. Here, the number of vertices is very low in order to ease the visualization of the embedding triangulation. In this simulation, the numerical solver only advects the vertices with a unit speed in the direction of the estimated normal. To be more precise, we obtain this velocity field by a variational paradigm, as the gradient ascent of the area enclosed by the interfaces.

Figure 6 shows different stages of the evolution. We invite the reader to have a close look at the connectivity and the labels of the embedded triangulation just before and after the two objects merge.

Three-phase combustion simulation

In our third experiment, we simulate the combustion of two different materials plunged in a third ambient material, to demonstrate the ability of our approach to accommodate any number of materials, and in particular multiple junctions. This type of motion generates complex phenomena called shocks and rarefactions in the study of conservation laws [START_REF] Leveque | Numerical Methods for Conservation Laws[END_REF]. Whereas shocks are extremely difficult to simulate with a traditional front-tracking method, due to the repeated formation of swallowtails, they are handled efficiently and robustly by our approach.

Figure 7 shows this simulation at various time samples. The combustion speed of the blue material is four times that of the red material, which creates a complex and interesting series of topology changes. In this figure, the tiling inside the materials is not represented, due to the large number of vertices. Anyway, we again emphasize that this tiling is only relevant to the internal functioning of the kinetic lazy triangulation.

Brain segmentation from medical images

In our last experiment, we give a glimpse of the possible applications of our approach in medical imaging. Deformable curves and surfaces are a widely used technique for shape reconstruction in image processing. We refer the interested reader to two reviews on this topic [START_REF] Mcinerney | Deformable models in medical image analysis: a survey[END_REF][START_REF] Xu | Medical image segmentation using deformable models[END_REF].

Here, we address the automatic delineation of the brain from magnetic resonance imaging (MRI) of the head. Figure 8 shows the 2D magnetic resonance image whose intensities drive the motion, the different stages of the evolution and the embedding triangulation of the final interface. Note how the initial seeds grow and progressively merge to fit the complex shape of the brain.

Quantitative evaluation

Table 1 gathers some quantitative measurements on our different experiments: r the number of vertices, at the beginning and at the end of the evolution, r the total number of events, r the minimum, maximum and average computation time required by one iteration of our algorithm. The time spent in the numerical solver and the time needed to update Although the numerical solvers used in our experiments are particularly simple and computationally inexpensive, the computation time required to maintain our lazy kinetic triangulation keeps lower. Our more complex simulation, brain segmentation, runs in about 300 milliseconds on a standard workstation.

Discussion and Future Work

We have proposed a robust and efficient Lagrangian approach for modeling dynamic interfaces between different materials undergoing large deformations. Thanks to the kinetic data structures framework and the exact computation paradigm, the reliability of our algorithm, even in difficult situations such as shocks and topology changes, can be certified. Being purely Lagrangian, our method does not suffer from mass loss which plagues the level set method, and can track material properties such as color or texture coordinates during motion at no additional cost, while being free of the localization problem of particle-based approaches. Whereas our work in its current progress does not claim state-of-the-art results, we have demonstrated in a series of numerical experiments that it brings a promising rigorous alternative to existing techniques for two-dimensional and axisymmetric simulations.

We now discuss the feasibility of an extension of our approach to three dimensions. As for the embedding triangulation, it immediately extends to a tetrahedralization, whose cells are labeled with the type of material they contain. All the advantages of this representation still hold, such as the automatic enforcement of watertight surfaces and the seamless handling of multiple junctions. However, as expected, maintaining this tetrahedralization is significantly more complex than in two dimensions. First, the predicates needed to certify the kinetic lazy triangulation are one orientation predicate for each tetrahedron. It requires to manipulate polynomials of degree 3, assuming a linear motion model, while a quadratic motion model now becomes prohibitively expensive.

But actually, the main difficulty is that more complex events must be processed to faithfully track the location of the interfaces. In particular, the addition of Steiner points cannot be avoided in some situations, which is consistent with the well known fact that many polyhedra have no tetrahedralization [START_REF] Ruppert | On the difficulty of triangulating three-dimensional nonconvex polyhedra[END_REF]. Thus, we cannot handle some events with the classical edge collapse, 2 → 3 and 3 → 2 tetrahedra flips. For instance, when two opposite edges of a flat tetrahedron intersect, it is sometimes mandatory to add the intersection point to the triangulation. Our future work includes investigating whether the number of these "spontaneous" Steiner points becomes prohibitive for the simulation.

Figure 1 :

 1 Figure 1: Embedding labeled triangulation.

Figure 2 :

 2 Figure 2: Example of the removal of a flat triangle by flipping its longest edge. The fill pattern represents the labels of the triangles.

Figure 3 :

 3 Figure 3: Example of an ineffective edge flip: it replaces two flat triangles with two new flat triangles in a similar configuration.

Figure 4 :

 4 Figure 4: Example of an acceptable edge flip: it replaces two flat triangles with two new flat triangles, but strictly decreases the length of the possible flipping edges.

Figure 5 :

 5 Figure5: Comparison between the level set method (top) and our method (bottom) for a circle stretched and distorted by a "vortex-in-a-box" velocity field.

Figure 6 :

 6 Figure 6: From left to right, different stages of the coalescence of two expanding objects.

  c 2007 The Authors Journal compilation c 2007 The Eurographics Association and Blackwell Publishing Ltd.

Figure 7 :

 7 Figure 7: (From left to right and top to bottom, different stages of a three-phase combustion simulation.

Figure 8 :

 8 Figure 8: Application to brain segmentation: (a) MR image, (b-g) different stages of the evolution, (h) labeled triangulation of the final shape.

Table 1 :

 1 Number of vertices, total number of events and measurements of the computation time in our different experiments.

	Comput. time / it. (ms)

c 2007 The Authors Journal compilation c 2007 The Eurographics Association and Blackwell Publishing Ltd.


Jean-Philippe Pons
email: pons@certis.enpc.fr

GEOMETRICA Jean-Daniel Boissonnat
email: boissonnat@sophia.inria.fr

Delaunay Deformable Models: Topology-Adaptive Meshes Based on the Restricted Delaunay Triangulation

In this paper, we propose a robust and efficient Lagrangian approach, which we call Delaunay Deformable Models, for modeling moving surfaces undergoing large deformations and topology changes. Our work uses the concept of restricted Delaunay triangulation, borrowed from computational geometry. In our approach, the interface is represented by a triangular mesh embedded in the Delaunay tetrahedralization of interface points. The mesh is iteratively updated by computing the restricted Delaunay triangulation of the deformed objects. Our method has many advantages over popular Eulerian techniques such as the level set method and over hybrid Eulerian-Lagrangian techniques such as the particle level set method: localization accuracy, adaptive resolution, ability to track properties associated to the interface, seamless handling of triple junctions. Our work brings a rigorous and efficient alternative to existing topology-adaptive mesh techniques such as Tsnakes.

Introduction

Deformable models, also known in the literature as snakes, active contours/surfaces, deformable contours/surfaces, constitute a widely used computerized technique to address various shape reconstruction problems in image processing. They have been initially proposed for the purpose of image segmentation by Kass, Witkin and Terzopoulos in [START_REF] Kass | Snakes: Active contour models[END_REF], but they have proven successful in many other contexts, in computer vision and in medical imaging, including region tracking and shape from X. More generally, modeling dynamic interfaces between several materials undergoing large deformations is a ubiquitous task in science and engineering. Let us mention computer aided design, physics simulation and computer graphics.

The existing techniques roughly fall into two categories: Eulerian and Lagrangian formulations. It is commonly ad-mitted that both viewpoints have strengths and weaknesses and that only a hybrid approach can overcome the limitations of both. In this paper, we propose a method which has the particularity of being purely Lagrangian with all the associated advantages, while achieving topological adaptivity with a comparable robustness and efficiency to Eulerian methods.

Eulerian methods

The Eulerian formulation casts deformation as a time variation of quantities defined over a fixed grid. The interfaces have to be represented implicitly, since the grid does not conform to them. In computational physics, this is also known as the front capturing method. Two notable front capturing techniques are the level set method, introduced by Osher and Sethian [START_REF] Osher | Fronts propagating with curvaturedependent speed: Algorithms based on Hamilton-Jacobi formulations[END_REF], and the volume-of-fluid (VOF) method, pioneered by Hirt and Nichols [START_REF] Hirt | Volume of fluid (VOF) method for the dynamics of free boundaries[END_REF].

Hereafter, we will mainly focus on the level set method because it is an established technique in computer vision. Basically, this method consists in representing the interface as the zero level set of a higher-dimensional scalar function. The movement of the interface can be cast as an evolution of the embedding level set function by an Eulerian PDE (partial differential equation). We refer the reader to some good reviews [START_REF] Osher | Level set methods: an overview and some recent results[END_REF][START_REF] Sethian | Level Set Methods and Fast Marching Methods: Evolving Interfaces in Computational Geometry, Fluid Mechanics[END_REF] for all the details about the theory, the recent developments, the implementation and the applications of the level set method.

On the one hand, this approach has several advantages over an explicit Lagrangian representation of the interface: no parameterization is needed, topology changes are handled automatically, intrinsic geometric properties such as normal or curvature can be computed easily from the level set function. Last but not least, the theory of viscosity solutions provides robust numerical schemes and strong mathematical results to deal with the evolution PDE. These advantages explain the popularity of the level set method not only in computer vision but also for multi-phase fluid flow simulation [START_REF] Sussman | A level set approach for computing solutions to incompressible two-phase flow[END_REF] in CFD (computational fluid dynamics), as well as for computer animation of fluids with free surfaces [START_REF] Enright | A hybrid particle level set method for improved interface capturing[END_REF][START_REF] Enright | A fast and accurate semi-Lagrangian particle level set method[END_REF][START_REF] Enright | Animation and rendering of complex water surfaces[END_REF][START_REF] Foster | Practical animation of liquids[END_REF][START_REF] Losasso | Simulating water and smoke with an octree data structure[END_REF].

On the other hand, several serious shortcomings limit the applicability of the level set method:

First, the higher dimensional embedding makes the level set method much more expensive computationally than explicit representations. Much effort has been done to alleviate this drawback, leading to the narrow band methodology [START_REF] Adalsteinsson | A fast level set method for propagating interfaces[END_REF] and to the PDE-based fast local level set method [START_REF] Peng | A PDE-based fast local level set method[END_REF]. More recently, octree decompositions have been proposed [START_REF] Cecil | Simplex free adaptive tree fast sweeping and evolution methods for solving level set equations in arbitrary dimension[END_REF][START_REF] Enright | A fast and accurate semi-Lagrangian particle level set method[END_REF][START_REF] Losasso | Spatially adaptive techniques for level set methods and incompressible flow[END_REF][START_REF] Losasso | Simulating water and smoke with an octree data structure[END_REF] to circumvent the typically fixed uniform sampling of the level set method, in order to reach high resolution (typically an effective resolution of 512 3 ) while keeping the computational and memory cost sustainable. However, these methods somewhat lose the simplicity of the original level set method, as an efficient implementation of such tree-based methods turns out to be a tricky task.

Second, as discussed and numerically demonstrated by Enright et al. in [START_REF] Enright | A hybrid particle level set method for improved interface capturing[END_REF], the level set method is strongly affected by mass loss, smearing of high curvature regions and inability to resolve very thin parts. These limitations have motivated the development of some hybrid Eulerian-Lagrangian methods, such as the particle level set method outlined by Foster and Fedkiw [START_REF] Foster | Practical animation of liquids[END_REF] and later improved by Enright and coworkers [START_REF] Enright | A hybrid particle level set method for improved interface capturing[END_REF][START_REF] Enright | A fast and accurate semi-Lagrangian particle level set method[END_REF][START_REF] Enright | Animation and rendering of complex water surfaces[END_REF]. While the latter method yields state-of-the-art results, an objection could be the large number of parameters controlling the particle reseeding strategy included in this approach.

Third, purely Eulerian formulation is not very appropriate for tracking interface properties such as color or texture coordinates, as may be needed in computer graphics applications. Some approaches based on a coupled system of Eulerian PDEs were recently proposed to overcome this limitation within the level set framework [START_REF] Adalsteinsson | Transport and diffusion of material quantities on propagating interfaces via level set methods[END_REF][START_REF] Pons | Maintaining the point correspondence in the level set framework[END_REF][START_REF] Xu | An Eulerian formulation for solving partial differential equations along a moving interface[END_REF], but this capability comes at a significant additional computational cost.

Previous Lagrangian methods

The Lagrangian formulation adopts a more "natural" point of view. It explicitly tracks the interfaces between the different materials with some points advected by the motion. There are mainly two classes of Lagrangian techniques: mesh-based methods and particle-based methods. We first tackle the mesh-based approach, which corresponds to the so-called "snake" methodology in computer vision. and which is also known as the front tracking method in computational physics. When dealing with large deformations, this approach is hampered by distortion and entanglement of the mesh, source of numerical instabilities or even of breakdowns of the simulation. For instance, if corners or cusps develop in the evolving front, front tracking methods usually form "swallowtail" solutions. These defective parts of the interface must be detected, then removed through intricate delooping procedures.

Another major shortcoming of the mesh-based Lagrangian approach is that a fully automatic, robust and efficient handling of topology changes remains an open issue, despite several heuristic solutions proposed in computer vision [START_REF] Bredno | A general discrete contour model in two, three, and four dimensions for topology-adaptive multichannel segmentation[END_REF][START_REF] Delingette | General object reconstruction based on simplex meshes[END_REF][START_REF] Delingette | Shape and topology constraints on parametric active contours[END_REF][START_REF] Duan | A subdivision-based deformable model for surface reconstruction of unknown topology[END_REF][START_REF] Lachaud | Deformable meshes with automated topology changes for coarse-to-fine 3D surface extraction[END_REF][START_REF] Lachaud | Deformable model with a complexity independent from image resolution[END_REF][START_REF] Mcinerney | Topology adaptive deformable surfaces for medical image volume segmentation[END_REF][START_REF] Mcinerney | T-snakes: Topology adaptive snakes[END_REF].

McInerney and Terzopoulos [START_REF] Mcinerney | Topology adaptive deformable surfaces for medical image volume segmentation[END_REF][START_REF] Mcinerney | T-snakes: Topology adaptive snakes[END_REF] propose topology adaptive deformable curves and meshes, called T-snakes and T-surfaces. During the evolution, the model is periodically resampled by computing its intersections with a regular simplicial decomposition of space. A labeling of the vertices of the simplicial grid as inside or outside of the model is maintained. This procedure loses the desirable adaptivity of the Lagrangian formulation, by imposing a fixed uniform spatial resolution. Also, not all motions are admissible: this approach only works when the model inflates or deflates everywhere, which considerably restricts the range of applications.

Several authors have proposed alternatives to the Tsnakes approach: Lachaud and coworkers [START_REF] Lachaud | Deformable meshes with automated topology changes for coarse-to-fine 3D surface extraction[END_REF][START_REF] Lachaud | Deformable model with a complexity independent from image resolution[END_REF], Bredno et al. [START_REF] Bredno | A general discrete contour model in two, three, and four dimensions for topology-adaptive multichannel segmentation[END_REF], Duan and Qin [START_REF] Duan | A subdivision-based deformable model for surface reconstruction of unknown topology[END_REF], Delingette and Montagnat [START_REF] Delingette | General object reconstruction based on simplex meshes[END_REF][START_REF] Delingette | Shape and topology constraints on parametric active contours[END_REF]. Basically, these approaches consist in detecting self-intersections in the evolving mesh and in merging the colliding regions using a set of heuristic remeshing rules. Unfortunately, the detection of intersection is computationally expensive, even when optimizing pairwise distance computations with an octree structure. In practice, it requires the most part of total computation time. Also, these approaches lack a systematic and provably correct remeshing strategy. Consequently, they are very likely to break in some complex or degenerate practical cases.

This major shortcoming of mesh-based methods have gained popularity to the particle-based approach [START_REF] Cani | Animation of deformable models using implicit surfaces[END_REF][START_REF] Desbrun | Smoothed particles: A new paradigm for animating highly deformable bodies[END_REF][START_REF] Desbrun | Animating soft substances with implicit surfaces[END_REF][START_REF] Foster | Realistic animation of liquids[END_REF][START_REF] Müller | Point based animation of elastic, plastic and melting objects[END_REF][START_REF] Pauly | Meshless animation of fracturing solids[END_REF][START_REF] Premoze | Particle-based simulation of fluids[END_REF][START_REF] Szeliski | Surface modeling with oriented particle systems[END_REF][START_REF] Terzopoulos | Heating and melting deformable models (from goop to glop)[END_REF] for representing dynamic interfaces undergoing complex topology changes. However, it is generally admitted that with the meshless approach, the localization of the interface and the computation of interface properties such as normal and curvature gets cumbersome.

Novelty of our method

In this paper, we present a purely Lagrangian approach, which combines the advantages of front tracking and front capturing methods, while discarding their respective drawbacks. Our work brings a robust and efficient solution to remeshing and topological adaptivity for Lagrangian dynamic interfaces, thanks to the concept of restricted Delaunay triangulation, borrowed from computational geometry. In our approach, the interface is represented by a triangular mesh embedded in the Delaunay tetrahedralization of interface points. Interestingly, this increase of dimension bears similarity to the spirit of the level set method. This embedding tetrahedralization enforces directly watertight in-terfaces free of loops, swallowtails or self intersections at all times. The mesh is iteratively updated by computing the restricted Delaunay triangulation of the deformed objects. Thus, our approach does not need heuristic rules to update the mesh connectivity. The proximity structure encoded in the Delaunay triangulation allows to detect topology changes efficiently and to handle them naturally.

Moreover, improving the geometric quality and adapting the resolution of the deformable model is far easier than in previous approaches, since the adequate connectivity is automatically obtained through the restricted Delaunay triangulation.

Being purely Lagrangian, our method does not suffer from mass loss which plagues the level set method, and can track material properties such as color or texture coordinates during motion at no additional cost, while being free of the localization problem of particle-based approaches.

Finally, our method is not limited to two-phase motion. It seamlessly accommodates any number of materials, whereas this requires special care in most existing Lagrangian and Eulerian methods (cf for example [START_REF] Smith | A projection method for motion of triple junctions by level sets[END_REF] on the treatment of triple junctions with the level set method). To achieve this, it suffices to label each Delaunay tetrahedron with the type of material it contains.

The remainder of this paper is organized as follows. Section 2 gives some background on the basic computational geometry concepts needed in our approach: Voronoi diagrams, Delaunay triangulations and restricted Delaunay triangulations. Our algorithm is described in Section 3. In Section 4, we demonstrate the applicability and the efficiency of our approach with some numerical experiments, including the multi-label segmentation of real medical images.

Background

Voronoi diagram and Delaunay triangulation

Voronoi diagrams are versatile structures which encode proximity relationships between objects. They are particularly relevant to perform nearest neighbors search and motion planning (e.g. in robotics), and to model growth processes (e.g. crystal growth in materials science). Delaunay triangulations, which are geometrically dual to Voronoi diagrams, are a classical tool in the field of mesh generation and mesh processing due to its optimality properties.

Most of the following definitions are taken from [START_REF] Boissonnat | Provably good sampling and meshing of surfaces[END_REF]. We also refer the interested reader to some recent computational geometry textbooks [START_REF] Boissonnat | Algorithmic Geometry[END_REF][START_REF] Berg | Computational Geometry, Algorithms and Applications[END_REF].

In the sequel, we call k-simplex the convex hull of k + 1 affinely independent points. For example, a 0-simplex is a point, a 1-simplex is a line segment, a 2-simplex is a triangle and a 3-simplex is a tetrahedron.

Let E = {p 1 , . . . , p n } be set of points in R d . Note that in this work, we are mainly interested in d = 2 and d = 3. Voronoi region, or Voronoi cell, denoted by V (p i ), associated to a point p i is the region of space that is closer from p i than from all other points in E:

V (p i ) = {p ∈ R d : ∀j, p -p i ≤ p -p j } . (1)
V (p i ) is the intersection of n-1 half-spaces bounded by the bisector planes of segments [p i p j ], j = i. V (p i ) is therefore a convex polytope, possibly unbounded. The Voronoi diagram of E, denoted by Vor(E), is the partition of space induced by the Voronoi cells V (p i ).

See Figure 1 for a two-dimensional example of a Voronoi diagram In two dimensions, the edges shared by two Voronoi cells are called Voronoi edges and the points shared by three Voronoi cells are called Voronoi vertices. Similarly, in three dimensions, we term Voronoi facets, edges and vertices the geometric objects shared by one, two and three Voronoi cells, respectively. The Voronoi diagram is the collection of all these k-dimensional objects, with 0 ≤ k ≤ d, which we call Voronoi objects. In particular, note that Voronoi cells V (p i ) correspond to d-dimensional Voronoi objects.

The Delaunay triangulation Del(E) of E is defined as the geometric dual of the Voronoi diagram: there is an edge between two points p i and p j in the Delaunay triangulation if and only if their Voronoi cells V (p i ) and V (p j ) have a non-empty intersection. It yields a triangulation of E, that is to say a partition of the convex hull of E into ddimensional simplices (i.e. into triangles in 2D, into tetrahedra in 3D and so on).

The fundamental property of the Delaunay triangulation is called the empty circle (resp. empty sphere in 3D) property: in 2D (resp. in 3D), a triangle (resp. tetrahedron) belongs to the Delaunay triangulation if and only if its circumcircle (resp. circumsphere) does not contain any other points of E.

The algorithmic complexity of the Delaunay triangulation of n points is O(n log n) in 2D, and O(n 2 ) in 3D. Fortunately, as was recently proven in [START_REF] Attali | Complexity of the Delaunay triangulation of points on surfaces: the smooth case[END_REF], the complexity in 3D drops to O(n log n) when the points are distributed on a smooth surface, which is the case of interest here.

Restricted Delaunay triangulation

Each k-simplex in the Delaunay triangulation is dual to a (d -k)-dimensional Voronoi object. For instance, in 3D, the dual of a Delaunay tetrahedron is the Voronoi vertex which coincides with the circumcenter of the tetrahedron, the dual of a Delaunay facet is a Voronoi edge, the dual of a Delaunay edge is a Voronoi facet, and the dual of a Delaunay vertex p i is the Voronoi cell V (p i ).

Given a subset Ω ∈ R d , typically a manifold of dimension k ≤ d, we call the Delaunay triangulation of E restricted to Ω, and we note Del| Ω (E) the subcomplex of Del(E) composed of the Delaunay simplices whose dual Voronoi objects intersect Ω. For example, in 2D, as illustrated in Figure 2 (left), the Delaunay triangulation restricted to a curve C is composed of the Delaunay edges whose dual Voronoi edges intersect C. Similarly, as shown in Figure 2 (right), the Delaunay triangulation restricted to a region R is composed of the Delaunay triangles whose circumcenters are contained in R. The attentive reader may have noticed that in both cases the restricted Delaunay triangulation forms a good approximation of the object.

Actually, this is a general property of the restricted Delaunay triangulation. It can be shown that, under some assumptions, and especially if E is a "sufficiently dense" sample of Ω, in some sense defined in [START_REF] Amenta | Surface reconstruction by voronoi filtering[END_REF], Del| Ω (E) is a good approximation of Ω, both in a topological and in a geometric sense: as regards topology, Del| Ω (E) is homeomorphic to Ω; as regards geometry, the Hausdorff distance between Del| Ω (E) and Ω can be made arbitrarily small; normals and curvatures of Ω can be consistently approximated from Del| Ω (E).

Based on these approximation properties, a family of provably correct algorithms for mesh generation and mesh reconstruction from point clouds have been designed in the last decade. We refer the reader to [START_REF] Boissonnat | Provably good sampling and meshing of surfaces[END_REF] and references therein for more details.

Methods

Embedding labeled Delaunay triangulation

A determinant feature of our approach is that the deforming objects and the interfaces between them are modeled as a subset of the Delaunay simplices of a point sample E.

More specifically, each tetrahedron of Del(E) is labeled with the type of material it contains. For instance, in a multi-phase fluid simulation, each tetrahedron would be mapped to one of the phases. This information, that we call for short the label of a tetrahedron, can be conveniently represented by an integer value. The interfaces of interest are embedded in the Delaunay triangulation: they are composed of the triangular facets adjacent to two tetrahedra having different labels.

This representation has many virtues. It automatically enforces watertight interfaces. It seamlessly accommodates any number of materials. It makes the request of the type of material at some point inexpensive, since there exist very efficient location algorithms dedicated to triangulations [START_REF] Berg | Computational Geometry, Algorithms and Applications[END_REF]. An example of this embedding Delaunay triangulation in 2D, including three materials with some triple junctions, is displayed in Figure 3.

Algorithm

We now describe how our embedding labeled Delaunay triangulation is consistently updated during motion, in order to faithfully reflect the deformations and the topology changes of the different materials. For sake of simplicity, the following description only deals with the case of two materials (inside and outside), the most common in practice. However, the described algorithm extends to any number of materials with a few minor changes.

Let us denote by E n the point sample at iteration n, and by D n the subset of the embedding Delaunay triangulation which is labeled as inside, i.e. foreground objects, at the corresponding time. The core idea of our approach is to obtain D n by taking the Delaunay triangulation of E n restricted to some region Ω n :

D n = Del| Ωn (E n ) . (2)
Given the approximation properties of the restricted Delaunay triangulation evoked in the previous section, this approach is valid as soon as E n is a good point sample of the interfaces and Ω n faithfully represents the interior of objects.

With this construction on hand, we only need to compute the point sample E n+1 at the subsequent iteration and the domain Ω n+1 after deformation, given the current E n and Ω n . As for the connectivity of interface meshes, it is automatically determined through Equation [START_REF] Adalsteinsson | Transport and diffusion of material quantities on propagating interfaces via level set methods[END_REF]. No heuristic and intricate remeshing rule is needed, as opposed to previous work [START_REF] Bredno | A general discrete contour model in two, three, and four dimensions for topology-adaptive multichannel segmentation[END_REF][START_REF] Delingette | General object reconstruction based on simplex meshes[END_REF][START_REF] Delingette | Shape and topology constraints on parametric active contours[END_REF][START_REF] Duan | A subdivision-based deformable model for surface reconstruction of unknown topology[END_REF][START_REF] Lachaud | Deformable meshes with automated topology changes for coarse-to-fine 3D surface extraction[END_REF][START_REF] Lachaud | Deformable model with a complexity independent from image resolution[END_REF].

E n+1 is obtained simply by advecting the points in E n with the desired velocity field, then by adapting the sampling resolution if needed. Specifically, for each pair of points of E n corresponding to an interface edge in D n :

• if the distance between the two displaced points is too small, the points are discarded and replaced by their midpoint in E n+1 ;

• otherwise, the displaced points are included in E n+1 ;

• if the distance between the two displaced points is too large, their midpoint is also added to E n+1 .

Of course, the aforementioned length thresholds may be dependent on the input data or on the geometry of the interface itself, e.g. on the structure tensor of the image as in [START_REF] Lachaud | Deformable model with a complexity independent from image resolution[END_REF], or on the curvature of the object. Also, we emphasize that, in contrast with existing approaches, the above resampling procedure can be applied without taking care of the connectivity of interface meshes. Again, the adequate connectivity after resampling is recovered through Equation ( 2). The ability of our approach to deal with topology changes is intimately related to the computation of the domain of objects after deformation. Ω n+1 is obtained as the union of tetrahedra in D n , after advecting their vertices with the velocity field, and after discarding tetrahedra which have become inverted. This can be written as

Ω n+1 = { T | T ∈ D n , orientation( T ) > 0} , (3)
where T denotes a tetrahedron and T is the same tetrahedron after displacement. For sake of clarity, the different substeps of an iteration of our algorithm are illustrated in Figure 4 (a-f), in two dimensions:

• (a) shows the embedding triangulation D n before deformation.

• (b) shows the triangulation after deformation, with its connectivity kept constant. Note that it is no longer a valid geometric triangulation, since there are some inverted triangles, dotted with yellow edges.

• (c) shows the domain Ω n+1 , i.e. the union of displaced triangles still having a correct orientation. Note that this domain accounts for the change in topology.

• (d) shows the Delaunay triangulation of the new point sample E n+1 along with the associated circumcenters, dotted with blue crosses.

• (e) shows the Delaunay triangulation of E n+1 restricted to the domain Ω n+1 . It is obtained by labeling Delaunay triangles as inside if the circumcenters displayed in (d) are inside the domain shown in (c).

• In some cases, especially after a topology change, some vertices of the new triangulation may be surrounded by simplices all sharing the same label. In the example shown in (e), four vertices are in this case. These vertices are of no use for representing the interfaces, so they are removed from the triangulation, yielding the final embedding Delaunay triangulation to be used in the next iteration, shown in (f).

As regards removing useless points in the last substep of each iteration, using a Delaunay triangulation rather than an arbitrary triangulation plays an important role in 3D. Indeed, although it may be impossible to tetrahedralize the hole created by removing a point from an arbitrary 3D triangulation [START_REF] Ruppert | On the difficulty of triangulating three-dimensional nonconvex polyhedra[END_REF], it is always possible in a 3D Delaunay triangulation, and there exist efficient algorithms to do it [START_REF] Devillers | Perturbations and vertex removal in a 3d delaunay triangulation[END_REF].

Numerical experiments

By using CGAL (Computational Geometry Algorithms Library, homepage: www.cgal.org) [START_REF] Boissonnat | Triangulations in CGAL[END_REF], we have been able to implement our approach with only 1000 lines of C++ code. CGAL defines all the needed geometric primitives and provides an excellent algorithm to compute the Delaunay triangulation in 3D: it is robust to degenerate configurations and floating-point error, thanks to the use of exact geometric predicates, while being able to process millions of points per minute on a standard workstation. As regards the efficient computation of the restricted Delaunay triangulation, which requires to test if a point intersects a "soup" of tetrahedra, we use a segment tree data structure [START_REF] Berg | Computational Geometry, Algorithms and Applications[END_REF], also implemented in CGAL.

Imbricated toruses

In our first experiment, we demonstrate that our approach naturally handles topology changes and resolution adaptivity, by segmenting a synthetic image of two imbricated toruses, starting from a sphere. The obtained evolution is shown in Figure 5. Note that we have purposely enforced different resolutions in the lower and upper parts of the object. Whereas this is straightforward to achieve with our approach, it is intricate when using Eulerian techniques such as the level set method.

In this experiment, the number of iterations to reach convergence is 120, the number of vertices ranges from 8,000 to about 22,000 along the evolution, and the total computation time is 171 s. Thus, the efficiency of our algorithm is comparable to the recent work of Lachaud and Taton [START_REF] Lachaud | Deformable model with a complexity independent from image resolution[END_REF]: roughly one second per iteration for 10,000 vertices. Comparing these timings to other related methods [START_REF] Bredno | A general discrete contour model in two, three, and four dimensions for topology-adaptive multichannel segmentation[END_REF][START_REF] Delingette | General object reconstruction based on simplex meshes[END_REF][START_REF] Delingette | Shape and topology constraints on parametric active contours[END_REF][START_REF] Duan | A subdivision-based deformable model for surface reconstruction of unknown topology[END_REF][START_REF] Mcinerney | Topology adaptive deformable surfaces for medical image volume segmentation[END_REF][START_REF] Mcinerney | T-snakes: Topology adaptive snakes[END_REF] turned out to be a tricky task, due to missing information such as number of iterations, number of vertices, or hardware specifications. Moreover, as the complexity of these algorithms is not linear, a fair comparison would suppose to use the same number of vertices, if not the same dataset.

Medical image segmentation

In our second experiment, we apply our approach to the segmentation of real 3D medical data. In order to demonstrate the ability of our approach to seamlessly handle any number of materials, we use three different labels: air, soft tissues (skin, muscles, . . . ) and bone. We set the velocity of the model to a combination of mean curvature motion and of a region-based data attachment term (see [START_REF] Mcinerney | Deformable models in medical image analysis: a survey[END_REF][START_REF] Xu | Medical image segmentation using deformable models[END_REF] for more details). Figure 6 shows the input CT image that drives the motion and the different stages of the evolution. Note how the initial random seeds corresponding to the different tissues grow and progressively merge, while creating triple junctions, until they successfully partition the head. Some movies of these experiments are available as supplemental material. Also, we plan to make our code available at the time of the conference.

Conclusion and future work

We have proposed a robust and efficient Lagrangian approach for modeling dynamic interfaces between different materials undergoing large deformations and topology changes, based on the rigorous concept of restricted Delaunay triangulation, borrowed from computational geometry. We have demonstrated the applicability and the efficiency of our approach with some numerical experiments, including the multi-label segmentation of real medical images. As our algorithm easily extends to any number of dimensions, our future work includes investigating its effectiveness in 4D, e.g. for dynamic scene reconstruction from multi-view image sequences, or for spatio-temporal MRI segmentation.

Figure 1 .

 1 Figure 1. Left: Voronoi diagram of a set of points in the plane. Right: Its dual Delaunauy triangulation.

Figure 2 .

 2 Figure 2. Left: The Delaunay triangulation restricted to the blue curve is plotted with a thick red line. Right: The Delaunay triangulation restricted to the region bounded by the blue curve is composed of the filled red triangles, whose circumcenters (blue crosses) are inside the region.

Figure 3 .

 3 Figure 3. Embedding labeled Delaunay triangulation, including three materials and triple junctions. The different materials are color-coded. Thick red edges indicate the interfaces between different materials.

Figure 4 .

 4 Substeps of one iteration of our algorithm (see text).

Figure 5 .Figure 6 .

 56 Figure 5. From left to right and top to bottom, different stages of the segmentation of two imbricated toruses, starting from a sphere, using a non-uniform sampling resolution.