Hervé Delingette
email: herve.delingette@inria.fr

General Object Reconstruction based on Simplex Meshes

Keywords: Image Segmentation, Deformable Models, 3D Reconstruction, Medical Imaging, Range Images

In this paper, we propose a general tridimensional reconstruction algorithm of range and volumetric images, based on deformable simplex meshes. Simplex meshes are topologically dual of triangulations and have the advantage of permitting smooth deformations in a simple and e cient manner. Our reconstruction algorithm can handle surfaces without any restriction on their shape or topology. The di erent tasks performed during the reconstruction include the segmentation of given objects in the scene, the extrapolation of missing data, and the control of smoothness, density, and geometric quality of the reconstructed meshes. The reconstruction takes place in two stages. First, the initialization stage creates a simplex mesh in the vicinity of the data model either manually or using an automatic procedure. Then, after a few iterations, the mesh topology can be modi ed by creating holes or by increasing its genus. Finally, a n i t e r a t i v e re nement algorithm decreases the distance of the mesh from the data while preserving high geometric and topological quality. Several reconstruction examples are provided with quantitative and qualitative results.

Introduction

General Object Reconstruction

From laser range nders to volumetric medical imagery, the development of three-dimensional acquisition devices have stressed the need for general shape reconstruction techniques. Object reconstruction is a task that consists in building a geometric model from a tridimensional dataset obtained from a scanning device. Such tridimensional datasets may vary substantially in terms of accuracy, resolution, or data structure. In all cases, the reconstruction task should create a ge-ometric model corresponding to the real object. The nature of this reconstructed model is dependent on the type of high level task performed a posteriori (see gure 1), such as visualization, object recognition, or scienti c computing (mechanical or ow analysis).

In this paper, we address the problem of general object reconstruction as opposed to sensordependent reconstruction techniques. Ideally, a general reconstruction system should handle volumetric images as well as tridimensional range data, with varying noise level and resolution. Furthermore, it should reconstruct smooth objects as well as polyhedral shapes with a control over the mesh density and the closeness of t. The gure recovery of objects from tridimensional datasets requires several tasks including: Segmentation: the scene described in the tridimensional dataset, is usually made-up of several objects lying close to each other. The segmentation task consists in isolating the object points from data points corresponding to neighboring objects or outliers.

Filtering: real tridimensional datasets always include some amount o f n o i s e . Therefore, it is often necessary to smooth the reconstructed model in order to reduce the e ect of noise without removing salient features such as sharp edges.

Extrapolation of missing data in many cases, 3D datasets do not entirely describe a given object. Therefore, it is necessary to handle those missing data points by extrapolating the surface in the most "intuitive" manner.

Mesh Density Control: when dense datasets are provided, it is often required, especially for visualization purposes, to greatly reduce the amount of information stored in the geometric model. On the contrary, when sparse data is provided, it is necessary to increase the mesh density b y re ning the geometric model.

Mesh Quality Control: many scienti c computation algorithms (continuum mechanics, ow analysis,: : : ) require the de nition of nite elements on meshes of high geometric quality. The quality of a mesh may be de ned in several terms 2]. For instance, on triangulated models, the geometric quality m a y be measured using the minimum, median or average angle of triangles.

Few existing algorithms address the "general reconstruction problem". For instance, isosurfac-ing 20] is a commonly used technique for reconstructing surfaces from volumetric or range images 10]. However, because it does not handle missing and noisy data, it cannot be considered as a general reconstruction technique. Similarly, Delaunay triangulation 1, 3] is widely used for serial slice reconstruction, or reconstruction from unorganized points but it cannot extrapolate missing data points.

Deformable modeling is well-suited for general object reconstruction because it makes few assumptions about the shape to recove r a n d i t c a n deal with missing and noisy data. There are several existing frameworks for deformable models, but a common approach consists of formalizing the deformation as a variational problem involving an internal energy ensuring the geometric continuity of the model, and an external energy controlling the closeness of t. Many researchers have proposed reconstruction systems based on deformable models [START_REF] Cohen | Using deformable surfaces to segment 3-d images and infer di erential structures[END_REF][START_REF] Delingette | Shape representation and image segmentation using deformable surfaces[END_REF][START_REF] Vemuri | Constructing intrinsic parameters with active models for invariant surface reconstruction[END_REF][START_REF] Mcinerney | A nite element model for 3d shape reconstruction and nonrigid motion tracking[END_REF][START_REF] Chen | Surface description of complex objects from multiple range images[END_REF] (for a survey of deformable models in medical image analysis see 25]). However, few systems address the problem of general surface reconstruction, including the ve tasks previously listed.

Deformable Models

A key issue for a general reconstruction system is related to the choice of a good surface representation. More precisely, s u c h representations should be well-suited for all surfaces independent of their geometry and topology. In this section, we review the existing geometric representations of deformable models and introduce our original surface representation.

The Parameterization Problem

Most reconstruction systems are based on parametric representations, such as splines or nite elements. A parametric representation provides a continuous transformation between a parameter space embedded in the Euclidean plane IR 2 and a tridimensional surface. The existence of a continuous representation enables the de nition of geometric quantities, such as a normal vector or curvature information, everywhere on the surface model. This is often required for high level tasks such as CAD design.

However, parametric representations su er from two problems. The rst problem is related to the representation of complex shapes. Because parameterizing a shape is equivalent to mapping a subset of the Euclidean plane onto that shape, problems occur when the object is not of planar, cylindrical, or toroidal topology. For spherical shapes for instance, at least one degenerate point or pole is created by the mapping of a plane onto a sphere. This point leads to many problems in the surface deformation because normal vector and curvature cannot be computed in a stable manner at that vertex. A lot of work has been done in the eld of computer aided geometric design to overcome these topological problems. A popular approach consists of sewing several parametric patches and ensuring that proper geometric continuity is realized between the patches. For instance Hoppe and Eck 15] use the construction scheme of Peters 29] to build G 1 continuous surfaces of arbitrary topology while other reconstruction systems such as Loop and De Rose 19] rely on di erent p a t c h corners such as Sabin nets. An important constraint in order to represent e c i e n tly general deformable models with spline patches is that the G 1 continuity equation across patches must be linear in the control points. Those continuity constraints across patches usually entail the addition of a large number of control points and furthermore tend to break the homogeneity o f t h e mesh 18].

The second problem lies in the nature of the parameterization which greatly in uences the deformation scheme. The in uence of the parameterization originates from the nature of the fairness functionals that are usually de ned in terms of parametric-dependent quantities such as partial derivatives. This is the case, for instance, of the widely used thin plate functional, based on rst order derivatives. When the deformation is constrained by parametric-dependent functionals, the reconstructed geometric model is dependent on the nature of the initial parameterization. In gure 2, we demonstrate the in uence of parameterization in a simple case where we are minimizing the following criterion:

E(r) = Z 1 0 k dr du k 2 du + kP ; r(u 0 )k 2 (1)
The choice of the parameter u 0 of the attachment point r(u 0 ) has a drastic e ect on the resulting approximation, better visual results being obtained when the initial parameterization is close to arc length. On the other hand, when the smoothness functional is intrinsic, i.e. independent of parameterization, the solution of the approximation problem is no longer dependent o f t h e c hoice of u 0 . In gure 2(right), we use the length of the curve a s the fairness functional of the approximation problem:

E intrinsic (r) = Z 1 0 k dr du kdu + kP ; r(u 0 )k 2 (2)
Details of these two a p p r o ximation problems are provided in appendix A.1. This example reveals the trade-o in choosing between parameterdependent and parameter-independent smoothness functionals. Parameter dependent functionals are easier to handle numerically, but the quality of the approximation depends on the metric distortion of its parameterization. On the other hand, intrinsic functionals guarantee a high level of smoothness but lead to more complex nonlinear problems.

In most research papers, little attention is paid to the errors introduced by metric distortions between the parameter space and the surface in ap-proximation problems. This is mainly because they limit themselves to objects with simple topology (planar or cylindrical). However, Eck and Hoppe 15] have addressed this reparameterization problem with the use of harmonic maps. Brechbuhler 5] reparameterizes closed surfaces in order to perform e cient object recognition.

As a conclusion, building parametric deformable models of arbitrary topology is a di cult task due to the constraints of having G 1 continuity across patches as well as keeping a parameterization with little distortion. It is not well suited for the task of general surface reconstruction due to the computational expense of performing topological changes such as creating holes or local re nement.

Non Parametric Deformable Models

To overcome the topological restriction of parametric deformable models, researchers have d e v eloped alternate representations. For instance, implicit deformable models, where the surface model is de ned through an implicit function f(x y z) = 0, h a ve the advantage of not restricting the topology, connectivity, and the geometry of the surface. In most cases, the function is de ned over a regularly tessellated domain. For instance, Mc Inerney and Terzopoulos 24] de ned topologically adaptable snakes through reparameterization on regular simplicial domain. Maladi and Vemuri 21] propose to deform a two dimensional and tridimensional implicit model by propagating fronts based on the formulation of Osher and Sethian. Muraki 27] and Gascuel et al. 38] propose implicit models based on radial functions to t range or medical data. Taubin et al. 37] t algebraic curves and surfaces of high degree on range data. Curless and Levoy 10] use distance maps to build a single model from multiple range images. If implicit model recovery has many advantages in terms of topological exibility, it leads to a computationally expensive scheme since it deforms the embedding space rather than the manifold itself. Furthermore, implicit models cannot represent surfaces with holes such as surfaces of planar topology.

Unstructured meshes, where no underlying parameterization is de ned, have been used as de-formable surface representation as well. The problem with using unstructured meshes lies in the definition of a stable and meaningful internal energy, or internal force. Vasilescu and Terzopoulos 39] use non-linear springs and masses on triangulated meshes for the reconstruction of range images. Even if this representation allows for local re nement and the detection of discontinuities, the numerical stability of spring models with non-zero rest length, is very sensitive to the mesh topology. In particular, if not enough springs are attached a vertex, the system becomes under-constrained and several rest shapes are possible. On the contrary, if too many springs are attached a vertex, the system is over-constrained and the bending of the spring model is very limited.

Other methods for regularizing a triangulation have been proposed. A traditional method for fairing triangulations is to apply Laplacian smoothing where each vertex is moved towards the center of its neighbors. Taubin 36] reduces the shrinkage of Gaussian smoothing by applying a l o w-pass lter. Recently, B o yer 4] proposed to minimize locally the area of the triangles adjacent to a vertex. Mallet 22] uses a smoothness energy on a triangulation de ned as a quadratic form of the neighboring vertex positions. These three internal energies provide intuitive deformation, but are only stabilizers of degree 1 and therefore tend to atten the curved parts of the triangulation. Welch and Witkin 41] use more sophisticated intrinsic functionals on triangulations by tting a local coordinate frame at each vertex, and then minimizing the functional locally. In fact, this method proposes new nite di erences expressions on an unstructured mesh. Despite the generality of this approach, it does not seem to bewell-suited for our purpose because it involves non linear optimization and numerical instability may occur when performing intensive re nements or large deformations.

Subdivision surface techniques 16, 33] allow t o describe smooth surfaces by iteratively subdividing an unstructured mesh. However, since the limit surface cannot be recovered explicitly but only approximately as a function of the control points, it often leads to computationally expensive s c hemes for approximating real data.

Oriented particles proposed by Szelinski 35] do not make a n y assumptions about the connectivity, or topology of the reconstructed model. However, it is not clear how robust is this method against outliers and narrow objects.

Simplex Meshes as Deformable Models

We propose an original surface representation, called simplex meshes as the basis of our reconstruction system. Simplex meshes can represent surfaces of all topologies, just as triangulation meshes. Furthermore, the geometry of simplex meshes enables to de ne at each v ertex, discrete geometric quantities such as mean curvature or normal vectors. We have de ned regularizing forces on simplex meshes allowing a high order of geometric continuity (C 2 continuity) in a simple and e cient manner. In addition, local and global topological operators are de ned to re ne a simplex mesh or change its genus, in a simple fashion. Because simplex meshes are not parametric models, it is not necessary to update a rigidity m a t r i x or the parameterization mapping when performing topological transformations.

In this paper, we p r e s e n t a shape recovery system based on this surface representation. In particular, we discuss the problem of initializing a simplex mesh, often considered as the weak point of deformable models. Also, we propose a semiautomatic algorithm for changing the mesh topology by creating holes or increasing the surface genus. Finally, w e describe a re nement algorithm based on the minimization of a geometric criterion based on the distance to the data or the local curvature.

The outline of the paper is as follows. In section 2, we brie y introduce the simplex mesh both in terms of topology as well as geometry. In section 3, we introduce the internal and external forces applied on these meshes in presence of range data or volumetric images. In section 4, we describe the various components of our reconstruction system including the initialization and topology control. In section 6, we present several examples of shape recovery and we conclude in section 7.

Simplex Meshes

The properties and de nitions of simplex meshes are only summarized in this section. Some com-plementary de nitions have been added as appendices and additional information can be found in 11]. Some preliminary results have been published in 12].

Definition of Simplex Meshes

The de nitions of simplex meshes and triangulations are closely related. More precisely, their underlying graphs are dual of each other. Another important property of simplex meshes is their constant vertex connectivity. In this section, we o n l y i n troduce the topological properties of simplex meshes. First, we g i v e a general de nition of k-simplex meshes embedded in a Euclidean space IR d of dimension d. Then, we consider 1 and 2-simplex meshes of IR 3 as representations of contours and surfaces. We de ne a k-simplex as a union of p-cells, 0 p k. Since these cells are p-simplex meshes, the de nition of a cell is recurrent: De nition 1. We de ne a 0-cell of IR d as a point P of IR d and a 1-cell of IR d as an edge of IR d , i . e . an unordered pair of distinct vertices (P M). We recursively de ne a p-cell (p 2) C of IR d as a union of (p ; 1)-cells such t h a t :

Definition of Cells and Simplex Meshes

1. Every vertex belonging to C belongs to p distinct (p ; 1)-cells. 2. The intersection of two (p ; 1)-cells is either empty o r i s a (p ; 2)-cell. Proposition 1. A k-simplex mesh is a (k + 1 )connected mesh: each vertex has (k + 1 ) neighboring vertices.

The constant connectivity b e t ween vertices implies a simple relation between the numberofvertices and the number of edges. Table 1 summarizes the connectivity b e t ween vertices, edges, faces, and cells of a k-simplex mesh. If a k-simplex mesh is (k + 1 ) -connected, all (k + 1 ) -connected meshes are not necessary simplex meshes. For instance, a 3-connected mesh where two faces intersect along two edges cannot be a 2-simplex mesh.

We write the set of n vertices of M as V (M) and its connectivity function as N(M). If P i is a vertex of a k-simplex mesh M then (P N0(i) P N1(i) : : : P N k (i) ) are its (k + 1 ) neighbors. A k-triangulation is composed of p-simplices (0 p k) which are the p-faces of the triangulation. 0-faces are the vertices, 1-faces the edges and 2-faces the triangles.

We can de ne a topological transformation that associates a k-simplex mesh to a k-triangulation. This transformation is pictured in gure 5 and considers di erently the vertices and edges located at the boundary of the triangulation from those located inside . Basically, this duality transformation associates a p-face of a k-triangulation with a (k ; p)-cell of a k-simplex mesh. In a triangulation, there is a real notion of boundary since a triangle may not be surrounded by three triangles. In a simplex mesh, a hole is simply an "empty" cell, since each v ertex is surrounded by (k + 1 ) k-cells.

There exists other duality transformations that have been de ned on k-triangulations. The most commonly studied has been the duality between triangulations and cellular complexes through the duality between Delaunay triangulation and Vorono diagrams 31]. Vorono diagrams are cellular complexes and the duality relation with Delaunay triangulation is geometric because it depends on the position of its vertices. On the contrary, the duality b e t ween triangulations and simplex meshes is purely topological since there are no geometric bijections between simplex meshes and triangulations.

The inexistence of geometric duality can be proven by considering the relative n umber of degrees of freedom of simplex meshes and triangulations. The geometry of a non-degenerate ksimplex mesh or a non-degenerate k-triangulation is determined by the set of coordinates of their vertices. However, for k > 1 the number of vertices V sm of a k-simplex mesh is di erent from the number of vertices V tr of a k-triangulation. For k = 2 for instance, for a triangulation without Table 3. Duality b e t ween a k-triangulation and a k-simplex mesh for boundary faces.

1-Tr

=) 1;SM 2-Tr =) 2-SM 1-SM =) 1;Tr 2-SM =) 2-Tr p = 0 Vertex =) Edge Vertex =) Face Vertex =) (nil) Vertex =) (nil) Vertex =) Vertex Vertex =) Edge p = 1 Edge =) Edge Edge =) (nil) Edge =) Vertex
holes having genus number g, the Euler relation gives:

V tr ; V sm 2 = 2 ( 1 ; g) Therefore, we cannot build an geometric homeomorphism between triangulations and simplex meshes since the vectorial space representing their space of possible con gurations have di erent d imensions. Only for 1-simplex meshes, is it possible to build a geometric dual 1-triangulation since they have the same number of vertices. Geometric equivalence with triangulations can only exist if we consider degenerate k-simplex meshes, for instance by restraining to convex polytopes (convex simplex meshes with planar faces).

The inexistence of a geometric transformation between simplex meshes and triangulations implies that the geometric deformation of a simplex mesh is not equivalent to the geometric deformation of a triangulation. This is the reason why simplex meshes are a surface representation distinct from triangulations.

2.1.3. Contours Contours de ned on 2-simplex meshes, are 1-simplex meshes, i.e. closed polygonal curves. They are simply de ned as a set of neighboring vertices such that a contour vertex has two and only two neighbors belonging to the contour (see gure 6). Contours can be de ned around any simplex-mesh faces especially faces corresponding to a hole. Contours are considered as deformable models moving independently from the mesh. The surface m e s h i s t h us attached to contour vertices which s e t the boundary conditions for the mesh deformation. The de nition of contours allows a greater control of a simplex mesh shape.

Mesh Transformation

Simplex meshes as triangulations are unstructured meshes and therefore can be locally re ned or decimated. In addition, simplex meshes can be cut along contours and surface handles can be created.

We de ne four basic topological operators acting on a simplex mesh, T 2 1 T 2 2 T 2 3 T 2 4 , described in gure 7. The rst two operators, are Eulerian since they do not modify the mesh genus. On the other hand, T 2 3 and T 2

Geometry of simplex meshes

In this section, we present the main geometric relations existing in a 2-simplex mesh of IR 3 . Similar results exist for 1-simplex meshes of IR 3 (i.e. tridimensional contours) and are described in appendix A.3. The main result consists of a simple equation (equation A1 and equation 7) giving the position of a vertex relatively to its neighborsand some geometric quantities: the simplex angle and the metric parameters.

On a tridimensional 2-simplex mesh M 2 IR 3 , a vertex P i is surrounded by three vertices (P N1(i) P N2(i) P N3(i) ). These three points de ne the tangent plane at P i with its normal vector n i :

n i = P N 1 (i) ^PN 2 (i) +P N 2 (i) ^PN 3 (i) +P N 3 (i) ^PN 1 (i) kP N 1 (i) ^PN 2 (i) +P N 2 (i) ^PN 3 (i) +P N 3 (i) ^PN 1 (i) k
We introduce the circle S 1 of center C i and radius r i circumscribed to the triangle (P N1(i) P N2(i) P N3(i) ). We introduce also the sphere S 2 of center O i and radius R i , circumscribed to the four vertices (P i P N1(i) P N2(i) P N3(i) ).

The simplex angle ' i = 6 (P i P N1(i) P N2(i) P N3(i) ) at P i is de ned by the two equations: ' i 2 ; ] : The simplex angle has many original properties that are developed in 11]. In particular, there is a simple relationship between the simplex angle ' i and the curvature H i = 1 Ri , also called the mean curvature at vertex P i :

sin(' i ) = ri Ri sign((P N1(i) ; P i ) n i ) cos(' i ) = kCi;Oik Ri sign((C i ; O i ) n i ) (3) O i P N2(i) F i C i r i P N1(i) P N2(i) n i P i R i (a) n i F i O i C i i d P i φ i r i R i (b)
H i = sin ' i r i (4)
From a vertex P i , w e i n troduce the orthogonal projection F i of P i onto the neighboring triangle (P N1(i) P N2(i) P N3(i) ). The metric parameters at a v ertex P i are the barycentric coordinates of F i with respect to the triangle (P N1(i) P N2(i) P N3(i) ):

F i = 1i P N1(i) + 2i P N2(i) + 3i P N3(i) (5) 1i + 2i + 3i = 1 (6)
The simplex angle with two metric parameters (since the three metric parameters are linked by equation 6 ) represent the position of P i with respect to its three neighbors: P i = 1i P N1(i) + 2i P N2(i) + 3i P N3(i) +L(r i d i ' i )n i [START_REF] Cohen | Using deformable surfaces to segment 3-d images and infer di erential structures[END_REF] where r i is the radius of the circumscribed circle at the triangle (P N1(i) P N2(i) P N3(i) ). d i is the distance between F i and the center C i of the circumscribed circle. L(r i d i ' i ) is de ned as:

L(r i d i ' i ) = (r 2 i ; d 2 i ) tan(' i ) p r 2 i + ( r 2 i ; d 2 i ) tan 2 (' i ) + r i (8) with: = 1 if j' i j < = 2 = ;1 if j' i j > = 2
The simplex angle and the metric parameters de ne the local shape around a given vertex. The simplex angle controls the local mean curvature, i.e, the height o f a v ertex with respect to the tangent plane. The metric parameters control the vertex position in the tangent plane with respect to its three neighbors. In another words, the metric parameters change the local parameterization whereas the simplex angle changes the extrinsic curvature of the surface.

Deformable Simplex Meshes

Law of Motion

In this section, we describe the law of motion for simplex meshes. Unlike parametric deformable models, the deformation of simplex meshes is not based on the notion of partial derivatives, but on the relative position of a vertex with respect of its neighbors, i.e. in terms of the simplex angle and the metric parameters.

As in most deformable model schemes, all vertices of a simplex mesh are considered as a physical mass submitted to a Newtonian law of motion including internal and external forces: m d 2 P i dt 2 = ; dP i dt + F int + F ext [START_REF] Cohen | A nite element m e t h o d a pplied to new active contour models and 3d reconstruction from cross sections[END_REF] where m is the vertex mass and is the damping factor. F int is the internal force assuring some level of geometric continuity. F ext is the external force constraining the distance between the mesh and some tridimensional dataset.

We discretized the time t in order to compute the evolution of the simplex mesh under the law of motion [START_REF] Cohen | A nite element m e t h o d a pplied to new active contour models and 3d reconstruction from cross sections[END_REF]. Using central nite di erences with an explicit scheme, the law of motion is discretized as:

P t+1 i = P t i + ( 1 ; )(P t i ;P t;1 i ) + i F int + i F ext (10)
Both forces F int and F ext are computed at time t. This explicit scheme is conditionally stable and therefore the two parameters i and i must belong to a given interval to guarantee a stable scheme. In equation 10, the internal and external forces have the dimension of a displacement.

Internal Force Computation

We propose an original formulation of internal forces based on the geometry of simplex meshes. The large extent of shape control is the main motivation for using simplex meshes rather than triangulations as a deformable surface representation (see section 5). The proposed internal force expressions are both geometrically signi cant and computationally e cient. However, they do not derive from the minimization of a global functional. Indeed, the geometric parameters (mean curvature, simplex angle, and metric parameters) at a vertex are related to the vertex position with a complex non-linear relationship. Therefore, the minimization of a global functional expressed in term of the geometric parameters, would lead to unnecessarily complex force expressions. We h a ve c hosen to use simpli ed regularizing force formulae at the expense of not having a global functional for guiding the minimization.

The internal force applied on the simplex mesh vertices can be decomposed into a tangential force F Tangent and a normal force F Normal .

F int = F Tangent + F Normal

Tangential Internal Force

The goal of the tangential force is to control the vertex position with respect to its three neighbors in the tangent plane. In general, it is preferable to have v ertices uniformly spread over the surface of a simplex mesh as in gure 9 (a). In some other cases, it is advantageous to have a concentration of vertices in places with high curvature as in gure 9 (b). At e a c h v ertex P i , w e h a ve de ned in section 2.2 the three metric parameters ( 1i 2i 3i ). In addition, we de ne three reference metric parameters (~ 1i ~ 2i ~ 3i ) corresponding to the prescribed value of the metric parameters. These two sets of metric parameters correspond to the barycentric coordinates of the two points F i and Fi . The tangential force F Tangent acting on a vertex P i is simply a spring force between F i and Fi :

F Tangent = Fi ; F i = ( ~ 1i ; 1i )P N1(i) + (~ 2i ; 2i )P N2(i) + ( ~ 3i ; 3i )P N3(i)
To obtain uniformly spread vertices, the reference metric parameters should all be equal ~ 1i = ~ 2i = ~ 3i = 1 =3. Otherwise, any strictly positive value of the reference metric parameters is possible for a stable deformation.

Normal Internal Force

The force acting in the normal direction constrains the mean curvature of the surface through the simplex angle (see equation 4). The general expression of F Normal is governed by the reference simplex angle 'i : regularizing force: the normal component is related to the variation of mean curvature i n a n e i g h borhood of vertex P i .

F Normal = ( L(r i d i 'i ) ; L(r i d i ' i )) n i F i N (i) 1 P i P P N (i) 3 F Tangent Fi F Normal 2 P N (i) (a)
The expression of 'i depends on the level of geometric continuity that should be enforced:

Position Continuity or C 0 constraint: We s e t 'i = ' i . In this case, F Normal = 0 and the surface can freely bend around vertex P i . This is useful to represent surface orientation discontinuities.

Surface Orientation Continuity or C 1 constraint: We choose 'i = 0 , w h i c h corresponds to moving each vertex towards its center of curvature. The total internal force is then equivalent to the mean curvature motion corresponding to the minimization of the total surface area. By using equation [START_REF] Cohen | Using deformable surfaces to segment 3-d images and infer di erential structures[END_REF], the internal force can be written as: F int = ( 1i P N1(i) + 2i P N2(i) + 3i P N3(i) ; P i ).

When metric parameters are equal to 1=3, the vertex P i is attracted towards the center of its neighbors and the internal force is equivalent to Laplacian Smoothing.

Shape Constraint: To constrain the local curvature of a vertex, we s e t 'i = ' 0 i , w h e r e = ' 0 i is a constant. When this constraint is applied on every vertex of a simplex mesh, this entails a global shape constraint up to a rotation, translation, and scale transformation.

Delingette

Simplex Angle Continuity or C 2 constraint : 'i is dened as an average of the simplex angles at neighboring vertices:

'i = X j2Q s i (Pi)
' j [START_REF] Delingette | Simplex meshes: a general representation for 3d shape reconstruction[END_REF] The neighborhood Q si (P i ) is recursively de ned as the union of Q si;1 (P i ) with the neighbors of Q si;1 (P i ) (we s e t Q 0 (P i ) = fP i g as shown in gure 11). The neighborhood size s i corresponds intuitively to the notion of rigidity, or deformation scale. In a highly rigid model, the model is smoothed over a large scale, and therefore an external constraint would entail a small de ection but over a large extent of the surface. The rigidity also has an e ect on the dynamics of the deformable model. This is because the scale of deformation in uences the speed of constraint propagation on the surface. A exible mesh with a small rigidity needs more iterations to reach its stable position than a mesh with high rigidity. In practice, we are using the C 2 smoothness constraint for shape recovery. Indeed, with the C 0 smoothness constraint, the simplex mesh is equivalent to spring-mass models and therefore only guarantees that a vertex stays close to its three neighbors. The C 1 smoothness constraint entails a continuity of surface normals but shrinks the surface towards its center of curvature. This shrinking e ect is a well known limitation of linear ltering (F int is then a linear combination of the vertex positions) and results in a bias estimate of the surface reconstruction. Indeed, the shrinking e ect implies that the mesh has a privileged direction of deformation. To l i m i t this bias, sev-eral algorithms has been proposed. For instance, a balloon force 8] has been added that arti cially pushes each vertex outwards. Another approach consists in creating an initial mesh that encloses the object to reconstruct. Instead, we propose to use a C 2 smoothness constraint for extrapolating missing data since it performs a good visual continuity of the mesh without any shrinkage. The maximum value of i is related to the stability criterion of the di erential equation [START_REF] Curless | A v olumetric method for building complex models from range images[END_REF].

When applying the C 0 or C 1 constraint, equation 10 can be seen as a straightforward Jacobi relaxation scheme and the maximum value of i is 1:0. The expression of the C 2 internal force is non-linear and we have found experimentally that i value below 1=2, ensures a stable iterative scheme.

Contour Deformation

Contours de ned on simplex meshes are deformable as well, and vertices belonging to those contours follow the same law of motion as in equation ( 9), but with a di erent expression of the internal force. The internal forces of tridimensional contours are described in appendix A.4 and are closely related to the formulae of section 3.2.1 and 3.2.2 ranging from position continuity t o c u r v ature continuity.

Because contour vertices are deformed independently of mesh vertices, they act as boundary conditions for the surface mesh. There are two t ypes of surface-contour constraints:

Simplex Angle Condition: The simplex angles at contour vertices are set to a given value.

Tangent Condition: The angle between the mesh normal and the contour normal measured around the contour tangent vector, intuitively corresponds to the angle between the surface and the contour.

Examples

The de nitions of C 0 , C 1 , C 2 and shape constraints on simplex mesh vertices and contour vertices, provide a powerful framework for reconstructing tridimensional objects. The variational modeling ability of simplex meshes proved to be a key feature for reconstructing surfaces of complex geometry and topology.

In gure 12, we s h o w a smooth junction between 4 four cylinders with C 1 and C 2 smoothness constraints. The junction has been created from a Fig. 12. Example of smoothness constraints applied to simplex meshes. We constrain a surface with four holes (left) to smoothly connect 4 cylinders with surface orientation continuity (center) and simplex angle continuity ( r i g h t).

sphere with four holes. The holes have been manually connected with the extremities of each cylinder.

In gure 13, we s h o w three examples of variational modeling based on simplex meshes and contours. The rst gure is an example of a minimal surface (zero mean curvature everywhere) based on simplex meshes. In gures 13 (b) and (c), we demonstrate the di erent boundary conditions between simplex meshes and contours.

External Force Computation

The expression of the external force F ext is dependent on the nature of the dataset and is given in sections 3.3.1 and 3.3.2 for range data and volumetric images. In all cases, the external force is directed along the normal direction n i at the vertex P i where the force is applied. The collinearity of F ext with n i is important f o r t wo reasons.

First, it entails smooth deformations with attraction forces varying smoothly along the mesh. Since we use an iterative process, it is important to guarantee a stable and smooth behavior over time. Indeed, tangential displacements could create vertices with negative metric parameters or even a self-intersecting mesh.

Second, the normal direction of the external force ensures the resulting mesh shape of the mesh will be smooth even in the presence of sparse data as shown in gure 14.

The complexity of the computation of F ext is linear in the number of mesh vertices, unlike Metaxas 17] and Terzopoulos 23] who compute the external force of all vertices with a complexity linear with the number of data points since they compute the closest mesh vertex from each data point. This method is not appropriate for general surface reconstruction since it makes the double assumption that all data points belong to the same object and that there are no outliers in the data. Furthermore, in presence of dense data, this method becomes too computationally expensive.

We propose an expression of the external force that deals with sparse as well as dense datasets and that can take into account outliers. We describe the computation of F ext for range data and volumetric images.

External Force Computation on Range Data

The computation of the external force is dependent on the notion of closest points. For every vertex P i of the mesh, we look for the closest data point M Cl (i) and the force F ext is then computed as: [START_REF] Delingette | Simplex meshes: a general representation for 3d shape reconstruction[END_REF] where i is a weight parameter ( 1:0), n i is the normal vector at P i , G(x) is the derivative o f The computation of M Cl (i) depends on the nature of the dataset and can be achieved with four algorithms:

F ext = i G (M Cl (i) ; P i ) n i D ref n i
1. Projection Method : On dense range images extracted from triangulation principles 34] or stereo-vision 14], the search for the closest point is theoretically of complexity O(m 2 ) for an m m image. However, when the calibration matrix is known, we propose a method in O(1), that approximates the search for the closest point. Indeed, we restrict the search along a two-dimensional segment in the image, projection of the tridimensional line passing through P i and directed along n i (see gure 15 (right)). This segment is centered around the projection of P i and is only a few pixels long, depending on the value of D ref .

Once the closest point along the segment h a s been found, we s e a r c h for the closest point i n a 5 5 window around that point.

Projection Method and Normal Orientation

: We can further restrict the choice of the closest point by considering the normal orientation of the range data in addition to the Euclidean distance. We then check if the dot product between the normal n i at the mesh and the normal at the range images is strictly positive.

3. KD-tree : When no calibration matrix is available, we compute the closest point with a kdtree 31]. This data structure gives the closest point inside a sphere centered around P i and of radius D ref , in nearly constant time.

Because of memory limitation, it is not always possible to store all data points inside a kd-tree. In the rst stage of the deformation where D ref is relatively large, we usually subsample the whole dataset in order to guarantee a fast computation of F ext . During the second stage, we discard all the data points that are far away from the mesh, and store all the remaining points in the kd-tree. The size of the kd-tree may still be large, but since the value D ref is substantially lower, the compu- tation time is still moderate.

External Force Computation on Volumetric

Images On volumetric images, the reconstruction task usually consists in isolating regions of consistent intensity values. Therefore, gradient intensity is the main information on which the external force is based. As in 7] and 23], we combine both gradient i n tensity and edge information for the computation of F ext : F ext = F grad + F edge [START_REF] Delingette | Shape representation and image segmentation using deformable surfaces[END_REF] The gradient i n tensity is used for local de ection of the mesh towards the voxels of maximum variation of intensity. The edge information, on the other hand, corresponds to gradient maxima and entails larger deformations of the mesh. The edge image is a binary image extracted from the gradient i n tensity image through thresholding. Previous work 23] has derived the computation of F grad as the gradient of the potential eld krIk 2 , t h us generalizing the approach o f a ctive contours. However, this formulation has the drawback of creating oscillations of the deformable model across the voxels with high gradient, due to the choice of the time step and the space discretization. Cohen 9] has proposed to normalize the gradient v ector in order to limit those oscillations within a pixel. However, this has the drawback of limiting the range of in uence of the gradient force to an arbitrary value.

We propose an expression of the gradient force that avoid all oscillations and that is not disturbed by local maxima of gradient i n tensity. The force F grad at vertex P i relies on the search in a neigh-borhoodofP i , for the voxel of maximum gradient intensity. If V is the voxel containing P i , then we inspect all voxels in an m m m window around V for the voxel center G i , of highest gradient intensity (see gure 16 (a)). The force expression is then:

F grad = grad i ((G i ; P i ) n i )n i (14)
where grad i is a weighting parameter with 0 grad i 1. Since the force is proportional to the true de ection vector G i ;P i (and not the gradient vector), this force does not entail any oscillations of the mesh. The computation of this force is actually fast since we can pre-compute the voxel center G i given the voxel V and a xed window size.

The gradient force can be made more speci c by incorporating additional constraints:

Gradient Direction Constraint: The voxel G i must have the highest gradient i n tensity with a gradient direction consistent with the normal direction at P i : rI(G i ) n i > 0. This constraint is natural since the gradient is always directed inwards or outwards from a region of constant i n tensity v alue.

Intensity Value:

The voxel G i must have a n i n tensity v alue in a given range value I min I max ]. This constraint helps to discriminate against neighboring regions of high contrast with di erent i n tensity values.

The edge image is created by thresholding the maxima of the gradient norm in the gradient direction. The removal of small connected components may help to obtain reliable edge information. The traditional approach for the computation of the edge force F edge has been to use distance maps built from the binary edge image 7, 23, 28]. Distance maps provide a simple algorithm for computing the direction of the edge force even at a long range. However, this approach has several limitations. First, it entails very large deformations away from the edge voxels, which causes an unstable behavior of the mesh model. Furthermore, distances are ambiguous at voxels equidistant from two e d g e v oxels. Our method uses the normal orientation at vertex P i to solve this ambiguity. Our approach consists in nding the closest edge voxel in the normal direction of the mesh.

N i P i F ext G i Voxel ν (a) P i i E Normal L Edge Voxel (b)
At vertex P i , we nd the closest voxel V, and a tridimensional line of voxels is scanned in the direction of n i (see gure 16 (b)). The maximum numb e r o f v oxels scanned is given by the reference distance D ref , determined as a percentage of the overall radius of the edge image. In general, less than 30 voxels are scanned for each v ertex. The tridimensional Bresenham line drawing algorithm is used to e ciently scan the voxel image along the line. If E i is the closest edge voxel along the normal line, then the edge force is given by: F edge = edge i (E i ; P i ) [START_REF] Eck | Automatic reconstruction of b-spline surfaces of arbitrary topological type[END_REF] The vector E i ; P i is collinear with n i and edge corresponds to the sti ness of the edge force, 0 edge i

1.

As with the gradient force, we can add constraints for nding the closest edge voxel : gradient direction constraint: The voxel E i must have the highest gradient i n tensity with a gradient direction consistent with the normal direction at P i : rI(G i ) n i > 0.

intensity value: The voxel E i must have a n i n tensity v alue in a given range value I min I max ].

Reconstruction Algorithm

Modeling Scheme

We propose a reconstruction scheme that is independent of the dataset type and that requires little a priori knowledge about the scene. The scheme is presented in gure 17 (a) and consists of two main stages.

The rst stage corresponds to the mesh deformation from an initial position to a close approximation of the dataset shape. The mesh is initialized with a spherical, cylindrical, or planar topology. The user can choose between a manual initialization or an automatic initialization algorithm described in section 4.2, if the dataset has few outliers.

In all cases, the mesh topology can be changed after the rst stage of deformation. It consists in creating holes, eventually followed by an increase of the surface genus, as seen in the diagram 17. All topological changes are performed semi-automatically with an automatic detection of topological problems, but a m a n ual validation of all operations.

Once the mesh has the desired topology, the second stage of deformation consists in providing a more accurate geometric representation of the object by re ning or adapting the mesh. Renement a n d adaptation are two complementary tasks that respectively add vertices and move v ertices towards parts of high curvature. Both tasks ensure that the nal mesh will be at a given distance from the data while having good geometric and topological properties.

Only the rigidity parameter r i and the external force parameter i are modi ed between the rst and second stage. During the rst stage, we set r i to high values (r i 10) and i to low values ( i 0:1) in order to obtain smooth and large scale deformations on the mesh. However, during the second stage, we u s e l o w v alues of r i (r i 1) and high values of i ( i 0:5) since a control of the t is required. rst stage of deformation: the mesh is brought from its initial shape to a rough approximation of the data (c) The second stage of deformation: the mesh is re ned or adapted in order to create an optimal representation of the object.

Initialization

One of the main problems with reconstruction systems based on deformable models is that the recovered shapes depend on the initial position of the mesh. In general, these models need to be initialized in a close vicinity of the dataset. When modeling complex shapes, manual initialization may not give proper results, despite the ability of changing the mesh topology.

We propose an initialization algorithm that can only handle a limited number of topologies (three in the current implementation) but that performs well even in the presence of missing data. Initial models can be created from both range and volumetric images. The algorithm is not sensitive t o noise but may be in uenced by outliers.

The creation of a spherical mesh takes place by rst computing the centroid G of the dataset, and then by nding the spherical surface centered on G enclosing the whole dataset. To do so, we consider a tessellated unit sphere centered on G, a n d then project each data point or voxel onto that sphere. We store on each surface element of the sphere, the maximum radial distance of projected points (see gure 18). When there are no points projecting on a sphere element, the mean value of the maximum projected radial distances is used instead. The variance of the maximum projected radial distances on the unit sphere determines the numb e r o f v ertices in the initial simplex mesh. Finally, the initialization of planar meshes is similar to the cylindrical case. Instead of the minimum inertia axis, we de ne the maximum inertia plane, passing through the centroid G and directed by the direction of maximum inertia (see gure 20). We project data points onto the plane, and compute for each planar element, the average height from the plane. Similarly, the variance of the height determines the number of vertices in the initial mesh. ). For the skull, holes has been created around the the orbits, the nose and the foramen. Other zones have not been cut due to their small size.

The face example shows why this procedure cannot be made fully automatic. Two zones have indeed been created: one around the upper part of the head and one at the neck level. The zone located at the hair level corresponds to missing data due to the limitation of the acquisition technology. On the other hand, the zone located around the neck does not correspond to any existing physical points and should be removed. This is why the hole creation algorithm requires some user interaction based on his a priori knowledge of the tridimensional scene.

Genus Modification

The genus of a surface is a topological characteristic de ned as the number of surface handles. When dealing with incomplete data, the initialization stage proposed in section 4.2 provides meshes with zero genus (sphere, cylinder, or plane). It is therefore necessary to provide an algorithm that increases the mesh genus.

We propose an algorithm that consist in connecting, with the topological operator T 2 4 , two contours surrounding two holes. Among all pairs of contours surrounding holes, we select those verifying the following criteria:

1. The centers of the two contours must be close to each other. 2. The diameters of the two contours must be similar. 3. The two contours must be located inside nearly parallel planes.

In gure 23, we initialize a mesh of spherical topology around a vertebra. The algorithm described in the previous section, automatically creates two holes on each side of the handle. If those contours verify the three criteria above listed, the genus is increased by applying the topological operator T 2 4 . The mesh is subsequently re ned to closely t the model. The Adaptation Algorithm Since simplex meshes are discrete representation of surfaces, the mesh goodness of t not only depends on the distance of vertices to the data, but also on the relative location of vertices on the surface object. A good strategy is to concentrate vertices towards parts of high curvature into order to optimize the shape description.

The reference metric parameters ~ i = f~ 1i ~ 2i ~ 3i g control the relative position a vertex with respect to its three neighboring vertices. Our algorithm updates the values of the reference metric parameters depending on the mean curvature values of neighboring vertices.

Basically, the algorithm operates as follows: vertices of low mean curvature migrate towards neighboring vertices of relatively higher mean curvature whereas vertices of constant mean curvature have metric parameters close to one third.

The reference metric parameters are updated every p iterations (p 1) b y the following relaxation algorithm:

~ t+p i = ~ t i + 1 2 (~ opt i ; ~ t i ))
where ~ t i is the reference metric parameters at iteration t and ~ opt i is the optimal metric parameter value based on the local variation of the mean curvature at P i .

To get a stable estimate of the mean curvature variation at a given vertex, we rst consider the mean value of the absolute mean curvature value jH i j:

jH i j = ( jH N1(i) j + jH N2(i) j + jH N3(i) j)=3
We then compute the relative mean curvature deviation vector jHj i as:

jHj i = 0 B B @ jH N 1 (i) j; jHij jHij jH N 2 (i) j; jHij jHij jH N 3 (i) j; jHij jHij 1 C C A
The vector jHj i is related to the variation of absolute mean curvature between one vertex and its three neighbors. The vector norm is null if the absolute mean curvature is locally constant. The optimal reference metric parameter value is related to jHj i with the following relationship:

~ opt i = 1 3 + i jHj i ( 16
)
where i is a constant that controls the extent o f the vertex adaptation and is usually chosen between 0.03 and 0.25. However, since reference metric parameters ~ ki must belong to the interval 0 1], w e m a y compute a value of i substantially lower than 0.03. In practice, we c hoose to update the metric parameters every p = 1 0 iterations in order to stabilize the mesh before re-evaluating the mean curvature over the mesh. Example of mesh adaptation is provided in gure 9 (a) and (b) on a mesh representing a cube. Figure 9 (a) corresponds to a mesh deformation without adaptation (reference metric parameters equal to one third) whereas gure 9 (b) corresponds to a mesh adaptation with i = 0 :20.

The Refinement Algorithm

Since the shape description of a simplex mesh is determined by the number of its vertices, we introduce a face re nement algorithm that controls the mesh resolution. Our algorithm is based on a re nement measure describing whether a face is a good candidate for re nement. Several geometric criteria may b e c o m bined to estimate the re nement m e asure.

More precisely, w e de ne the following four geometric criteria as re nement measures:

The relative face area : It is computed as the ratio of a face area to the total mesh area. The face elongation : It is computed as 1:0 plus the di erence of length between the longest edge and the smallest edge of a face, divided by the median value of face edges length. The face Gaussian curvature : It is evaluated as the area of the spherical polygon described by the face normal vectors on the Gaussian sphere. The face distance to a 3D dataset : It is based on the ratio between the smallest distance to a 3D dataset with the cut-o distance D ref .

The re nement algorithm works in an iterative manner (see gure 24 (a)) at every q iterations (q = 5 in general). We rst sort all faces according to their re nement measure such that faces having high re nement measures are rst processed. Each face having a re nement measure above a given threshold, is re ned by creating a new edge with the topological operator T 2

2 . The mesh is then locally smoothed around each newly created edge and neighboring edges are eventually swapped in order to keep the numberofvertices per face close to six. Finally, w e recompute the re nement m e asure for all re ned faces and iterate until all faces have a re nement measure below the threshold. This algorithm is fairly e cient because all operations are performed locally. Figure 24 shows an example of mesh re nement based on the distance to a synthetic cross dataset.

We often choose the curvature and distance to 3D dataset as the adequate criteria of re nement. However, when the 3D dataset is very noisy, or has many outliers, the distance criterion should n o t b e u s e d s i n c e i t w ould result in actually tting the noise. Finally, it is often useful to restrict the re nement algorithm to given parts of the mesh in order to optimize the total number of vertices of the simplex mesh.

Simplex meshes versus Triangulated meshes

Simplex meshes and triangulated meshes are topologically dual of each other. This duality implies that all topological operations de ned on a simplex mesh can be de ned on a triangulation and inversely. However, in terms of data structure, simplex meshes are slightly easier to work with because of the constant vertex connectivity and because borders are considered as empty faces of a simplex mesh. Note that borders of triangulations do not correspond to triangles and therefore must be handled in a separate data structure.

The main di erence between simplex meshes and triangulation is geometric in nature. O n a v ertex of a simplex mesh, surface normal, mean curvature, simplex angle, and metric parameters are unambiguously de ned. We have devised a simple and e cient i n ternal force on simplex meshes solely based on these geometric quantities. Efciency is an important issue for shape recovery because it allows some level of user interaction.

Furthermore, we h a ve found that the C 2 smoothness constraint is the appropriate constraint for extrapolating missing data since the C 1 smoothness constraint tends to shrink surfaces towards their center of curvature.

Several simple and e cient s c hemes have been proposed for deforming triangulations under the C 1 continuity constraint (a good example can be found in 30]). However, as seen in section 1.4, enforcing the C 2 continuity constraint on nonparametric triangulated representations requires either sophisticated spline formulations 26] or subdivision surfaces [START_REF] Hoppe | Piecewise smooth surface reconstruction[END_REF][START_REF] Stam | Exact evaluation of catmull-clark subdivision surfaces at arbitrary parameter values[END_REF]. Those representations lead to more computationally intensive algorithms and are not well-suited for a fast update of their data structure after a local (mesh re nement) or global (genus increase) topological change. Indeed, a key advantage of simplex meshes lies in their internal regularizing force that is only based on the vertices position without any additional degrees of freedom. Therefore, all topological operations can be done e ciently with a simple data structure.

However, a disadvantage of using simplex meshes over triangulations is that they must be triangulated in order to be displayed or transfered to other high level tasks. In general, for triangulating a simplex mesh, the dual triangulation is not a good choice because it strongly underestimates the vertex curvature. Instead, it is preferable to triangulate each face independently for instance by adding a central vertex. In this section, we compare the e ciency of the Newtonian law of motion with the Lagrangian law of motion. We recall that the discrete Newtonian law of motion is: P t+1 i = P t i + ( 1 ; )(P t i ;P t;1 i ) + F int +F ext [START_REF] Koh | Hierarchical shape representation using locally adaptive nite elements[END_REF] We note that if = 1, then equation 17 corresponds to a Lagrangian law of motion. As the damping factor decreases towards zero, the effect of the acceleration increases. For = 0 , t h e deformable model behaves as a perfect oscillator. We have tested the e ect of the damping factor on the deformation of a simplex mesh on a range image of a foot. The mesh is initialized as an ellipsoid and we perform 30 iterations of deformation.

The initial mesh position and the deformed mesh with a value of = 0 :80 are shown in gure 25. In gure 26, we h a ve studied the displacement of vertices during the deformation for several values of the damping factor. The displacement is de ned as the sum of vertex motion between two iterations. It appears that as decreases, the convergence speed increases. However, when is smaller than 0:20, large oscillations prevent the mesh from converging towards the shape of the range data. For = 0:20, we observe that the resulting mesh self-intersects. The results show that the Newtonian law o f m otion is more e cient than the Lagrangian law o f motion. However, the choice of the damping factor must be governed by the trade o between e ciency and stability. If is too small, the mesh may not converge towards the right shape, especially when the mesh is far away from the data. If is too close from 1.0, then the e ciency of the deformation process is poor. In practice, we choose a value of = 0 :65.

Performance Evaluation of Closest-Point Algorithms

We have introduced two methods for nding the closest data point to a given vertex. The projection method uses the calibration matrices of camera-based range-nders to nd the closest point along the normal line. This method does not guarantee that we will nd the true close s t d a t a p o i n t, but it is fairly computationally ecient. On the other hand, the kd-tree data structure, gives the true closest point from a vertex, but at a higher memory and computation cost.

In this section, we compare the accuracy and the computational e ciency of these two methods at three di erent stages of deformation: initialization, the rst deformation and after re nement. The three meshes are shown respectively in gure 25 (a) and gures 27 (a) and (b).

Since the kd-tree method gives the true closest point to the data, we e v aluate the accuracy of the projection method in comparison with the distance given by the kd-tree method. All distances are measured along the vertex normal direction, in order to remove the e ect of sparse data. In gure 4, we can see that the accuracy of the projection method, increases drastically as the mesh is closer to the range data. The maximum errors occur at parts where the range data is incomplete (such as the tip of the foot) and therefore does not have an important in uence on the overall shape of the deformable model.

The computational cost of both methods are compared in gure 4. It appears that the eciency of the kd-tree method is linked with the radius D ref of search around the vertex which i s not the case for the projection method. When the mesh is far from the data, the kd-tree method is 35 times slower than the projection method whereas it is only 2:5 times slower for lower value of D ref .

In all cases, the time needed to build a kd-tree from a 512 512 range image was 4:03s.

From this experiment, we conclude that when the calibration parameter of the range sensor are known, the projection method is better-suited than the kd-tree method because of its low computational and memory requirement. Furthermore, the poor accuracy of the projection method when far from the data is balanced by the low v alue of the sti ness parameter i ( i = 0 :1).

Effect of the Refinement on the Distance Error

In this section, we study the e ect of re nement o n the overall distance error. Starting from a rough estimate shown in gure 25 (a), the mesh position after applying the rst stage of deformation is shown in gure 27 (a). This model has 1280 vertices and over-smoothes the shape of the original data.

We re ne this mesh based on its distance to the dataset. In this example, we choose two di erent re nement l e v els set by two re nement measures. The former re nement measure has been automatically computed in order to re ne 25% of the original number of faces whereas the latter corresponds to a re nement level of 40%. Since the re nement measure is linked to the maximum distance to the data, choosing a re nement m e asure is equivalent t o c hoosing a maximum value for that distance. The rst re nement measure corresponds to a maximum distance of 22.10 mm and the second to a distance of 14.25 mm. Note that this distance is computed not only for all mesh vertices, but also for the all face centers. Meshes before and after re nement are shown in gure 27. Quantitative results on the distance to the original data, are given in gure 5. The re nement procedure is called every 5 iterations and stops when there are less than 5 faces to re ne. We limit the number of faces re ned per procedure to 100, in order to obtain a more homogeneous spacing of vertices. The number of faces re ned per call is described in gure 28 (a). The total number of iterations increases from 45 to 130 as the required re nement l e v el is modi ed. In gure 28 (b), we can see the change of the distance error histogram before and after re nement. As shown in gure 5, the maximum error is within the bounds set by t h e re nement criterion. Note that the median error does not vary substantially between the two l e v els of re nement since it approximately corresponds to the range data resolution.

The computation times were evaluated on a DEC Alphastation 200/233. We indicate in gure 5 t h e time needed for the rst stage of deformation (where the mesh is deformed from an ellipsoid) as well as for the re nement stage. Within the re nement stage, most of the time is spent o n evaluating the re nement criterion for all faces.

Qualitative Results

Changing the Connectivity

In this section, we demonstrate that our reconstruction scheme can take i n to account t h e c hange of connectivity. More precisely, starting from a simplex mesh, the creation of holes, as described in section 4.3.1, can cut a simplex mesh into two independent meshes. The change of connectivity occurs when the part of the zone interpolated has more than one border.

In this example, we use a pre-segmented volumetric image consisting of manually segmented contours. The di culty lies in building a smooth geometric model that closely ts those contours. Direct reconstruction techniques, such as Delaunay triangulation 3] or the Marching Cubes al-gorithm, extrapolate the contours without taking into account a n y smoothness constraints.

In our framework, the deformable mesh is attracted by the volumetric image and submitted to regularizing forces. To compute the attraction force, we simply consider the contours as edge voxels. We start the reconstruction with an automatic mesh initialization of spherical topology around the data. We obtain a mesh surrounding both lungs (see gure 29 (a)). The rst deformation stage creates a smooth model where we can correctly detect the part which has been interpolated. Since this interpolated part has two borders, its removal entails the creation of a second simplex mesh model ( gure 29 (c)). Those meshes are then re ned and the nal reconstruction is pictured in gure 29(f).

Reconstruction from Volumetric Images

In this example, we reconstruct the ventricles, atria and the pericardium from a T1-weighted MRI volumetric image. The image corresponds to a routine acquisition of the abdomen and consists of only 10 slices without any contrast agent. Because of the low-contrast, deformable models are the only alternative to manual segmentation. Furthermore, we cannot automatically initialize our models because it is impossible to isolate the ventricles and atria in the image. Instead, the ventricles and atria models are initialized as ellipsoids ( gure 30 (a)) and positioned manually within the volumetric image. The pericardium is initialized as a cylinder. All models are attracted towards edges with the additional constraint that the surface normal must be coherently oriented with the image gradient d irection. The gradient information was computed as a 2D gradient since the inter-slice distance is six times larger than the pixel size.

First, the di erent models are deformed with a high rigidity parameter and, in a second stage, the rigidity parameter is decreased to a minimum while slightly re ning the meshes based on area. Because of the rough mesh initialization, we allow the user to locally guide the mesh in order to remove it from a local minimum. However, we found that using the coherence between the gradient direction and the orientation of the deformable model greatly helps the segmentation. The right and left ventricles have been easily segmented ( gure 30 (b) and (c)) whereas the recovery of the atria is more problematic. In particular, the right and left atrium intersect each other because the septum separating the two cavities is missing in the image ( gure 30(h)). In gure 30(i), we s h o w the slices of the di erent reconstructed models. A better t could be achieved by increasing the resolution of the models at the expense of computation time. The number of vertices for the right/left ventricles, the right/left atrium and the pericardium is respectively 1920, 1280, 1320, 1592 and 366.

Reconstruction by Parts

When the object to recover is strongly non-convex or non starshaped, it may be necessary to recover sub-parts of that object and subsequently join the di erent parts.

For instance, we have reconstructed a hand model by reconstructing independently the palm and the 5 ngers. Those models have been initialized with a sphere and 5 cylinders, as seen in gure 31 (a). After recovering each piece, we m a n ually connect them together with several T 2 3 operators. The junctions between meshes are then automatically smoothed to ensure a continuity of normal orientation. The nal model has about 8000 vertices and uses the texture of both range data. The exibility of the simplex mesh representation allows us to perform the reconstruction by parts in a straightforward manner.

Conclusion

In this paper, we h a ve p r e s e n ted a general reconstruction framework based on deformable simplex meshes. It di ers from previous approaches by using a non-parametric representation of surfaces and by the addition of deformable contours lying on the surface model. The deformable models can be re ned, adapted, and the level of smoothness can be controlled in a computationally e cient manner. Because simplex meshes can represent surfaces of all topologies, we h a ve proposed an algorithm for adapting the mesh topology to the topology of the object. Finally, w e h a ve proposed new formulations of the external forces from range data or volumetric images.

In the current implementation, we rely on an initialization algorithm that is limited to only three di erent topologies. Furthermore, it requires that objects of interest have been previously isolated from other objects. We i n tend to develop a more general initialization technique that does not require any a priori knowledge of the object topology. Such an approach m ust be coupled with a low-level segmentation algorithm that can extract a region of interest. Finally, when the object to recover is known a priori, the robustness of segmentation and reconstruction could be improved by including in the deformable model a notion of shape statistics. This would speci cally be bene cial for segmenting medical images where databases of volumetric images and anatomic structures could be used.

Fig. 1 .

 1 Fig. 1. The general reconstruction scheme

Fig. 2 .

 2 Fig. 2. The parameterization problem : (left) The curve minimizing the approximation functional of equation 1 with u 0 = 0 :2 and = 2 :0 (center) Same as (left) with u 0 = 0 :8 (right) The curve minimizing the intrinsic approximation functional of equation 2 with = 2 :0.In gure 2, we demonstrate the in uence of parameterization in a simple case where we are minimizing the following criterion:

Fig. 3 .

 3 Fig. 3. Examples of p-cells 0 p 3.

A 2 -Fig. 4 .

 24 Fig. 4. Four generic 2-simplex meshes with di erent v alues of genus and number of holes.

2. 1 . 2 .

 12 Duality with Triangulations It is important to stress the dual nature between ksimplex meshes and k-triangulations. The ktriangulations, also called k-manifolds 32], are sets of p-simplices (0 < p k) that follow strict topological rules such as the Euler-Poincar relation. The k-triangulations are actually a subset of more general sets of p-simplices, called k-simplicial complexes.

Fig. 5 .

 5 Fig. 5. Duality b e t ween k-triangulations (k = 1 2) d r a wn with solid lines circles, and k-simplex meshes drawn with dashed lines and squares. The rst two gures correspond to triangulations without boundaries whereas the last two have boundaries.

Fig. 6 .

 6 Fig. 6. Two c o n tours de ned on a 2-simplex mesh. The rightmost contour surrounds a face of the mesh.

Fig. 7 . 2 1 T 2 2 2 3

 7222 Fig. 7. (a) The two Eulerian operators T 2 1 T 2 2 de ned on 2-simplex meshes (b) The two meta-operators T 2 3 and T 24 ( c ) The edge swap operator T2 7 .

Fig. 8 .

 8 Fig. 8. (a) The circumscribed sphere S 2 of radius R i and the circumscribed circle S 1 of radius r i . (b) Projection of gure (a) onto the plane (O i , C i ,P i ). The simplex angle c a n b e i n terpreted as an angle of planar geometry.

Fig. 9 .

 9 Fig. 9. (a) Deformable simplex mesh with reference metric parameters equal to 1=3 (b) Deformable simplex after adaptation of its reference metric parameter in order to increase the resolution in places with high curvature (c) De nition of the tangential force as ( Fi ; F i ).

Fig. 10 .

 10 Fig. 10. (a) Geometric interpretation of the C 1 regularizing force: in this case the internal force can be written as F int = Fi ; P i (b) Geometric interpretation of the C 2

Fig. 11 .

 11 Fig. 11. Description of a neighborhoodQ s i (P i ) around a vertex P i with s i = 0 (see (a)) , s i = 1 (see (b)) and s i = 2 (see (c))

Fig. 13 .

 13 Fig. 13. (a) a simplex mesh under C 1 continuity constraint (with null mean curvature everywhere) having the symmetry of a dodecahedron (b) a simplex mesh under C 2 continuity constraint and four contours the surface-contour constraints of the two i n termediate contours correspond to curvature constraint (zero mean curvature) and tangent constraint (c) a vase created from a cylinder and ve contours

Fig. 15 .

 15 Fig. 15. (left) The function G(x) corresponding to the derivative of an M-estimator (right) Search for the closestpoint on a range image using the projection method Pi and ñi are the projection of P i and n i on the image plane.

Fig. 16 .

 16 Fig. 16. (a) Computation of the gradient force F gradient : G i is the center of the voxel with the highest gradient i ntensity in the neighborhood ofvoxel (b) Search for the closest edge voxel along the normal line for the computation of F edge .

Fig. 17 .

 17 Fig. 17. (a) The general reconstruction scheme (b) Therst stage of deformation: the mesh is brought from its initial shape to a rough approximation of the data (c) The second stage of deformation: the mesh is re ned or adapted in order to create an optimal representation of the object.

Fig. 18 .Fig. 19 .

 1819 Fig. 18. (a) For initializing a spherical mesh, each data point is projected on a sphere (b) The initial simplex mesh build from a parametric surface on spherical coordinates

Fig. 20 . 4 . 3 . 1 .

 20431 Fig. 20. (a) The initialization of a planar mesh is based on the projection of data points on the plane of maximum inertia (b) The initial simplex mesh build from the elevation map

Fig. 21 .

 21 Fig. 21. Hole creation on a skull model: (a) Mesh after the rst stage of deformation (b) Color coding of the distance of mesh vertices to the dataset (c) Zones to be cut (d) Mesh after the hole creation algorithm We present t wo examples of hole creation with the reconstruction of a skull model (see gure 21) and face model (see gure 22). For the skull, holes has been created around the the orbits, the nose and the foramen. Other zones have not been cut due to their small size.The face example shows why this procedure cannot be made fully automatic. Two zones have indeed been created: one around the upper part of the head and one at the neck level. The zone located at the hair level corresponds to missing data due to the limitation of the acquisition technology. On the other hand, the zone located around the neck does not correspond to any existing physical points and should be removed. This is why the hole creation algorithm requires some user interaction based on his a priori knowledge of the tridimensional scene.

Fig. 22 .

 22 Fig. 22. Holes creation on a face model : (a) Simplex mesh after the rst stage of deformation (b) The zone corresponding to the missing data around the hair (c) Zone corresponding to the neck (d) Mesh after cutting the neck zone

Fig. 23 .

 23 Fig. 23. (a) Deformation of a mesh initialized as a sphere on a vertebra dataset (b) Color coding of the distance to the data (c) The hole creation algorithm generates two holes (d) The two c o n tours are merged and the mesh has a toroidal topology

Fig. 24 .

 24 Fig. 24. (a) Flow c hart of the re nement algorithm (b) A simplex mesh t on a synthetic cross without re nement (1620 vertices) (b) Mesh with a low level of re nement (1970 vertices) (c) Mesh with a high level of re nement (4020 vertices)

1 .

 1 Influence of the Damping Factor Our deformable model framework is based on a Newtonian law of motion (see equation9), that includes a damping factor in order to prevent oscillations of the system. Most deformable model systems, on the contrary, are based on a Lagrangian law o f motion, where the speed of each v ertex is equal to the sum of the internal and external forces.

Fig. 25 .

 25 Fig. 25. (a) Initialization of the mesh around the foot range data (b) Mesh after 30 iterations with = 0 :80

Fig. 26 .

 26 Fig. 26. (a) Chart of the mean vertex displacement a s a function of the iteration number for several values of the damping factor (b) The mean and median distance to the range data for several values of the damping factor after 30 iterations.

Fig. 27 .

 27 Fig. 27. (a) The mesh of gure 25 (a) after the rst stage of deformation (b) Mesh after ne re nement.

Fig. 28 .

 28 Fig. 28. (a) The number of faces re ned for each call of the re nement procedure (b) The distance error histogram before and after re nement.

Fig. 29 .

 29 Fig. 29. (a) Initial mesh (b) A zone detected as interpolated (c) Change of connectivity (d) and (e) Reconstructed model of the rst and second lung (f) and (g) Final reconstruction (h) and (i) Slices of the lung models.

Fig. 30 .Fig. 31 .

 3031 Fig. 30. (a) Initial meshes (b) and (c) The right and left reconstructed ventricles (d) The right and left atrium (e) (f) and (g) The 5 recovered models from the MRI image (h) and (i) Slices of the di erent models.

Table 1 .

 1 The connectivity relations of a k-simplex mesh.

		Edges / Vertex Faces / Vertex Faces / Edge
	k=1 k=2	2 3	3	2

Table 2

 2

	summarizes the transformation between a p-face of a k-triangulation and a (k ; p)-cell of a k-simplex mesh. For cells or faces belonging to a triangulation boundary or a simplex mesh boundary, the dual transformation applies di er-ently. Vertices, edges, and triangles belonging to a k-triangulation boundary are associated with two cells of a k-simplex mesh: one (k ;p)-cell and one (k ; p ; 1)-cell (0 p < k ) (see table 3).

Table 2 .

 2 Duality between a k-triangulation and a ksimplex mesh for internal faces.

	1;Tr () 1;SM 2;Tr () 2;SM
	p = 0 Vertex () Edge Vertex () Face p = 1 Edge () Vertex Edge () Edge p = 2 Triangle () Vertex

Table 4 .

 4 Accuracy and computational cost of the projection and kd-tree method.

	Accuracy Rate	Mean Error	Maximum Error	D ref	Computation Cost of the Projection Method	Computation Cost of the KD-tree
	Initial Mesh Deformed Mesh 95.2 % 33 % Re ned Mesh 98.7 %	36.93 mm 371 mm 380 mm 0.38 mm 45.89 mm 95 mm 0.04 mm 14.32 mm 95 mm	0.79 ms 0.81 ms 0.81 ms	27.20 ms 2.21 ms 2.02 ms

Table 5 .

 5 Distance error and computation time for three meshes at di erent re nement levels.

	Mean Distance Distance Distance Median Maximum Numberof Numberof Computation Vertices Iterations Time
	Deformed Mesh Coarsely Re ned Mesh 3.56 mm 2.44 mm 21.98 mm 16.01 mm 14.53 mm 75.70 mm Finely Re ned Mesh 2.94 mm 2.27 mm 14.07 mm	1280 1518 2020	30 20 50	23.77 s 35.73 s 124.43 s

are meta-operators because they can break a mesh into two pieces, or change its genus (number of surface handles). All topological transformations can be decomposed into a set of those four operators.

Acknowledgements

The author would like to thank Y. Watanabe from NTT Yokosuka Research Laboratory for providing the Cyberware range images as well as Pr Darcourt from the Centre Antoine Lacassagne in Nice (France) for providing the MRI heart images. Also, the author would like to thank all the members of the Epidaure project for interesting discussions regarding the topic of this paper.

From the Euler-Lagrange equation associated with equation 1, it can be proved that the curve minimizing equation ( 1) is piecewise linear. If we note A and B the extremities of the curve, P the data point and r 0 the vertex of parameter u 0 , t h e n the problem lies in nding r 0 that minimizes: E(r 0 ) = kA ; r 0 k 2 u 0 + kB ; r 0 k 2 1 ; u 0 + kr 0 ; Pk 2 which leads to the following equation: r 0 = u 0 (1 ; u 0 )P + u 0 B + ( 1 ; u 0 )A 1 + u 0 (1 ; u 0 )

Replacing the sum of square derivatives with the total curve length, removes u 0 in the expression of the optimal curve: E intrinsic (r 0 ) = kA ; r 0 k + kB ; r 0 k + kr 0 ; Pk 2 which leads to the following non-linear equation that does not have a n y closed form solution: r 0 = 2 kB ; r 0 kkA ; r 0 kP + kA ; r 0 kB + kB ; r 0 kA kB ; r 0 k + kA ; r 0 k + 2 kB ; r 0 kkA ; r 0 k

Geometry of Planar 1-Simplex Meshes

The geometry of a planar 1-simplex mesh M 2 IR 2 is described by its metric parameters and its simplex angle.

We rst de ne the local tangent t i and normal vector n i around a vertex P i as:

t i = P i+1 ; P i;1 kP i+1 ; P i;1 k n i = t ?

The local curvature i is de ned as the inverse of the radius of the circumscribed circle at (P i;1 P i P i+1 ) (see gure A.2). We d e n e a s w ell r i as the half distance between P i;1 and P i+1 : r i = kP i+1 ; P i;1 k=2. We t h e n de ne the metric parameters and the simplex angle at P i .

De nition 3.

Let P i be a vertex of a planar 1-simplex mesh. Let F i be the orthogonal projection of P i onto the line P i+1 ; P i;1 ]. Then, the metric parameters at P i , ( 1i 2i ) are de ned as Delingette the barycentric coordinates of F i relative t o P i;1 and P i+1 :

De nition 4. The simplex angle ' i at P i is the oriented angle existing between the two adjacent segments P i;1 P i ] and P i P i+1 ].

The fundamental relation gives the position of vertex P i as a function of its two neighbors P i;1 and P i+1 and the three shape parameters ( 1i 2i ' i ):

The values of the metric parameters and the simplex angle describe the shape of the mesh up to a translation, rotation, and scale transformation.

Geometry of Tridimensional 1-Simplex Mesh

Contours are tridimensional 1-simplex mesh M 2 IR 3 , where we de ne a local frame made of the tangent t i , normal n i and binormal b i : t i = Pi+1;Pi;1 kPi+1;Pi;1k The metric parameters and the simplex angle are de ned as for planar meshes. However, we need to introduce another angle i as additional shape parameter. This parameter is de ned through vector r i : r i = t i ^((P i;1 ; P i;2 ) ^(P i+2 ; P i+1 )) kt i ^((P i;1 ; P i;2 ) ^(P i+2 ; P i+1 ))k (r i and t i ^ri ) de nes the normal plane, orthogonal to t i . Therefore, we de ne i as the angle between n i and r i : n i = c o s ( i )r i + sin( i )t i ^ri vectors. F i is the foot of P i on P i;1 P i+1 .

The shape of a tridimensional simplex mesh is then fully described by its metric parameters, simplex angle and angle i :

with r i = 1 2 kP i+1 ; P i;1 k and d i = j2 1i ; 1jr i .

Internal Force Applied on Contours.

The internal force applied on a tridimensional contour C can also be decomposed into a tangential and normal force.

F int = F Tangent + F Normal For the tangential force, we de ne a set of reference metric parameters (~ 1i ~ 2i ) which corresponds to the position in the tangent plane of a vertex with respect to its 2 neighbors. If P i is a vertex of C and if P i;1 , P i+1 are its two neighbors, then the tangential force is expressed as: F Tangent = ( ~ 1i ; 1i )P i;1 + ( ~ 2i ; 2i )P i+1 With the notation of appendix A.3, the normal force is written as:

The choice of ~ i and ' determines the level of continuity constraint.

Position Continuity or C 0 Constraint: . We set 'i = ' i and ~ i = i . The normal force is then null and the contour can freely bend around each v ertex.

General Object Reconstruction Based on Simplex Meshes 141 Normal Continuity or C 1 Constraint: We set 'i = 0 and ~ i = 0. The force then simply writes as F int = ~ 1i P i;1 + ~ 2i P i+1 ; P i which corresponds to Laplacian smoothing .

Curvature Continuity or C 2 Constraint: We set ~ i = 0 and 'i = 'i;1+'i+'i+1 3 . Shape constraint: We set 'i = ' 0 i and ~ i = 0 i where ' 0 i and 0 i are two constants corresponding to the reference shape.