Eric Ferley
email: eric.ferley@imag.fr

Marie-Paule Cani-Gascuel
email: marie-paule.gascuel@imag.fr

Dominique Attali
email: dominique.attali@imag.fr

Skeletal Reconstruction of Branching Shapes

We present a new method to reconstruct an implicit representation of a branching object from a set of data points scattered on its surface. The method is based on the computation of a geometric skeleton inside the data set. This skeleton is simplified in order to filter noise and converted into skeletal elements -a graph of interconnected curves -that generate an implicit surface. We use Bézier triangles as extra skeletal elements to perform bulge free blends between branches while controlling the blend extent. The result is a smooth reconstruction of the object, that can be computed whatever its topology. The skeleton offers compact storage, and provides an underlying structure for the reconstructed object, making it easier to edit in a modeling or animation environment.

Introduction

With the recent developement of advanced rangeimaging sensors, automatic reconstruction of realworld objects from scattered data points has become an important issue in Computer Graphics. Applications include medical imaging, inverse engineering, simulation, and semi-automatic modeling for design or animation. For the first set of applications, the aim is to provide an improved and enhanced visualization of the data. In other cases, the reconstructed object is going to be edited, deformed and animated. Here, providing a structured reconstruction is essential.

In the area of shape recognition, researchers structure a data set by computing its skeleton, defined as the locus of the centers of maximal spheres inside the data [START_REF] Blum | A transformation for extracting new descriptors of shape[END_REF][START_REF] Arcelli | Euclidean skeleton via centre-of-maximal-disc extraction[END_REF] . This skeleton, a thin centered structure that throughout this paper will be called geometric skeleton, provides a compact representation of both the topology and the geometry of the shape. Such a notion of a geometric skeleton appears attractive in the † iMAGIS is a joint project of CNRS, INRIA, Institut National Polytechnique de Grenoble and Université Joseph Fourier.

area of implicit surface modeling. The questions to be asked are: does this skeleton represent the same notion as the skeletal elements used to generate an implicit surface [START_REF] Bloomenthal | Interactive techniques for implicit modeling[END_REF][START_REF] Bloomenthal | Convolution surfaces[END_REF][START_REF]Bulge elimination in convolution surfaces[END_REF] ? Can we use it to construct a smooth implicit representation of an object? Compared with previous approaches for reconstruction using implicit surfaces, the introduction of a geometric skeleton to generate the surface gives important benefits:

• The implicit surface can be computed in a purely geometric way, without the need of optimization processes such as in [START_REF] Muraki | Volumetric shape description of range data using blobby model[END_REF][START_REF] Bittar | Automatic reconstruction of unstructured 3d data: Combining medial axis and implicit surfaces[END_REF] . These processes are often computationally intensive and difficult to control. • The skeleton, a graph of interconnected curve segments or surface elements, defines a structure for the reconstructed object. It can be used to edit the shape in an natural way. Such a structure is not provided when objects are reconstructed by directly computing an implicit function [START_REF] Whitaker | Algorithms for implicit deformable models[END_REF][START_REF] Velho | Approximate conversion of parametric to implicit surfaces[END_REF] , or by merely generating a list of skeletal points [START_REF] Muraki | Volumetric shape description of range data using blobby model[END_REF][START_REF] Bittar | Automatic reconstruction of unstructured 3d data: Combining medial axis and implicit surfaces[END_REF] .

In consequence, an approach based on skeletal reconstruction seems more adapted to modeling and animation applications.

This paper focuses on the reconstruction of free form branching shapes from the geometric skeleton of data points scattered on the object's surface. Conversion from geometric skeleton to parametric representation has already been widely studied [START_REF] Gelston | Boundary surface recovery from skeleton curves and surfaces[END_REF] . The aim of this paper is different since we want to build a smooth implicit surface representing the shape. We call "branching shapes" these objects for which the geometric skeleton can be represented as a graph of interconnected curves, without any surface element. The topology of these objects is arbitrary Section 2 reviews the definition of the geometric skeleton of an object and presents the semi-continuous approach that was used to compute it from a data set. Section 3 points out that this skeleton cannot directly be used to generate an implicit surface, and lists the problems to solve. Section 4 presents a method to convert each branch of the geometric skeleton into a skeletal curve that generates a smooth implicit surface of varying radius in the cross-section. Section 5 develops a new method, based on Bézier triangular skeletal elements, to smoothly blend branches at joints without generating bulges or creases. Results are presented in Section 6. We conclude and focus on future work in Section 7.

Geometric Skeleton

Definitions

The geometric skeleton of an object is the locus of the centers of maximal spheres inside this object. A sphere included in an object is said to be maximal if it is not included in any other sphere included in the object.

The geometric skeleton is a thin structure, represented as a graph of interconnected curve segments or surface elements, centered in the object (see Figure 1). It stores in a compact way both the topological properties of the object -given by its graph structureand the geometrical information, i.e. the distance to the surface. Note that the geometric skeleton is based on the same principles as the "medial axis" used in discrete geometry [START_REF] Lam | Thinning methodologies -a comprehensive survey[END_REF] . However, the medial axis of the data set is built as a volume made of an unstructured set of voxels. In contrast, the graph structure of the geometric skeleton gives useful neighboring information.

A large amount of work has been done in image analysis on the computation of the skeleton [START_REF] Blum | A transformation for extracting new descriptors of shape[END_REF][START_REF] Arcelli | Euclidean skeleton via centre-of-maximal-disc extraction[END_REF][START_REF] Brandt | Continuous skeleton computation by Voronoi diagram[END_REF] . Skeletonization methods are quite different depending on the way objects are represented. If the object is defined as a collection of voxels (or pixels in 2D), then the skeleton can be computed either by simulating the propagation of a wire initialized at the boundary of the object, by topological thinning, or by first extracting the medial axis from a distance map. In this paper, as the object is only known through a set of unorganized points {pi} n i=1 located on its boundary, a different approach must be considered.

Computation of the geometric skeleton

Our approach exploits the relationship between the skeleton of an object and the Voronoi graph of its boundary points. The Voronoi graph of a set of points (called seeds) divides the space into regions [START_REF] Aurenhammer | Voronoi diagrams -a survey of a fundamental geometric data structure[END_REF] . Each region is the set of points closer to a particular seed than to any other seed. The skeleton may be approximated as the subset of the Voronoi elements (the edges of the Voronoi diagramm) that are completely included in the object [START_REF] Brandt | Continuous skeleton computation by Voronoi diagram[END_REF] (see Figure 2). The better the sample of the object, the better the approximation of the skeleton. The computation of the skeleton requires to determine which part of the Voronoi graph is located inside the object. As the initial data is a mere set of boundary points with no structure, there is no immediate way to decide whether a point is inside or outside the object. Therefore, an heuristic must be used. We first compute the Delaunay tetrahedrization of the boundary points pi, which is the dual of the Voronoi graph. Then, Delaunay tetrahedra are progressively merged until space is partitioned into two sets that represent c The Eurographics Association 1997 the interior and the exterior of the object 14 (see Figure 4(a)). After this step, the skeleton can easily be retrieved from the tesselation of the object.

The result is a thin shape made up of straight-line segments in 2D and polygons in 3D (Figures 3(a) and 4(b)). Each Voronoi vertex located on the geometric skeleton is the center of a sphere passing through at least three seeds and having no seed in its iterior. The radius of this sphere approximates the distance to the surface at each vertex.

Pruning branches due to noise

A drawback of the skeleton transformation is its sensibility to noise. Noise on the boundary of an object may significantly change its skeleton (Figure 3 The skeleton is simplified by iteratively removing those Voronoi vertices that are located on the border of the skeleton and verify a removing criterion. By construction, the simplified skeleton is a subset of the initial skeleton. Moreover, the simplification method does not change the homotopy type of the skeleton, and thus the topology of the object remains unchanged.

In this paper, we apply a removing criterion based on the dual relationship between the object and its skeleton 15 16 . Since the center and the radius of each maximal sphere is known, the skeleton stores the information on the boundary of the object. Hence, the amount of change in this boundary gives a criterion to guide the simplification process: we remove a Voronoi vertex from the skeleton if the volume of its dual Delaunay simplex is less than a given threshold. This prunes branches due to noise, and in the 3D case, can lead to wireframe skeletons, which are "natural skeletons" for branching shapes -made up of generalized cylinders (see Figure 4(c)). By construction, the order of all branching points is three.

Problems Raised by Skeletal Implicit Modelling

People working in the area of implicit modeling use a set of geometric primitives -a skeleton -associated with distance varying field functions to generate a surface. Ideally, a smooth surface representing an object should be directly computed from any skeletal description. Unfortunately, using the geometric skeleton described in Section 2 for iso-surface generation is not so trivial, since this skeleton is a graph of very small line segments.

This section first reviews methods used in skeletal implicit modeling, and then points out the difficulties raised by the use of the geometric skeleton computed in Section 2.

Implicit surfaces generated by skeletons

Implicit surfaces have been defined as iso-surfaces of field functions (Sometimes called "implicit function" although they are defined by an explicit equation) generated by a set of geometric primitives that define the "skeleton" of the object [START_REF] Bloomenthal | Interactive techniques for implicit modeling[END_REF][START_REF] Bloomenthal | Convolution surfaces[END_REF] .

The field contribution fi of each skeletal primitive Si is a decreasing function of the distance to this element. The implicit surface S is defined as an isosurface of isovalue c for the sum of skeletal contributions:

S = P | i fi(P ) -c = 0 (1)
This formulation is a generalization of the original Blinn's objects [START_REF] Blinn | A generalization of algebraic surface drawing[END_REF] , metaballs [START_REF] Nishimura | Objects modeling by distribution function and a method of image generation (in japanese)[END_REF] and soft objects [START_REF] Wyvill | Data structure for soft objects[END_REF] , where only skeletal points were used.

If the summation is replaced by a maximum in equation ( 1), the implicit volume is the union of the volumes created by each skeletal primitive (see Figure 5(c)). The advantage of summation is that it defines a smooth blend between these contributions. However, the way each fi varies affects the blend so that generating a smooth shape is not an easy task. In particular, bulges (see Figure 5(b)), characterized by a "cross section that exhibits negative, then positive, then negative curvature with respect to the underlying skeleton" [START_REF]Bulge elimination in convolution surfaces[END_REF] , should be avoided.

The next paragraphs review the different models that have been proposed as field functions.

Distance surfaces

Distance surfaces 3 define fi as a decreasing C 1 continuous function F of the distance to Si (defined as the distance to the closest point on Si):

fi(P ) = F (d(P, Si)) (2)
The choice of a function F with local support limits the influence of each skeletal component, thus providing local control of the shape and reducing computations. Anisotropic field functions depending on the location of the closest point on Si may be used to generate more complex shapes [START_REF] Kacic-Alesic | Controlled blending of procedural implicit surfaces[END_REF] .

Distance surfaces can be safely used only with convex skeletons, otherwise creases may appear. Dividing non-convex skeletons into convex parts does not give any solution, since bulges at joints would be produced in this case. Consequently, modeling a ramification without any bulge or crease appears difficult (see Figure 5).

Convolution surfaces

Convolution surfaces 4 define fi as the integral over Si of contributions from individual points. Here, the resulting surface does not depend on the way a skeleton c The Eurographics Association 1997 is divided into skeletal elements, which is more satisfactory. Moreover, the resulting surface is smooth, even for a non-convex skeleton. However, J. Bloomenthal showed that the convolution method generates bulges at a joint of one dimensionnal skeletal pieces 5 : one must thicken the skeletal pieces and handle 2D (or 3D for circular crosssections) skeletons to prevent that phenomenon.

Procedural methods

Procedural field functions were introduced as an example of complex implicit definition 3 , in the specific case of a skeleton made up of two branching curves S1 and S2. In order to compute the field contribution at a given point P , the points P1 and P2 closest to P on each of the two skeletal curves Si are computed, and then fi(P ) is defined by the distance to the closest point on the segment [P1, P2]:

fi(P ) = F (d(P, [P1, P2]))
Unfortunately, this approach that seems to offer an elegant solution to the bulge problem cannot be applied in the general case:

• It does not work for sharp angles between branches, as shown in Figure 6, since the surface tangent changes suddenly when a projected point P1 or P2 is located on the extremity of the line segment. • There is no control of the blend extent: it is restricted by the projective domains of points P on Si. In consequence, two consecutive ramifications can not be modeled with this method.

Problems to be solved

Applying one of the previous skeletal modeling techniques to reconstruct a smooth shape from the geometric skeleton of a data set is not straightforward:

Each branch of the geometric skeleton computed in Section 2 is made up of many segments. Due to noise in the data, the resulting polyline is not smooth. The distance surface generated by this polyline would exhibit creases, while the summation of each segment's contribution would generate bulges. As we want to exploit the simplicity of the 1D skeletal parts computed, the use of convolution surfaces does not fit our problem. None of the previous methods of field computation can directly be applied to the skeleton we obtained. We have to find a way to smoothly blend consecutives line-segments and smooth junctions between several branching curves.

The rest of this paper presents our solutions to these problems:

P1 Pt P2 P -4 -3 -2 -1 0 1 2 3 -4 -2 0 2 4 (a) -4 -3 -2 -1 0 1 2 3 -3 -2 -1 0 1 2 3 (b)
Figure 6: Effects of the procedural method to blend the contributions of two line segments: (a) The method works when the angle between the segments is greater than π 2 . (b) Strange features appear for sharp angles.

In Section 4, we approximate each polyline by a spline curve, and define a field function of varying strength along it in order to generate an implicit generalized cylinder of varying radius.

Our solution to the branching problem is described in Section 5. It was inspired by the procedural method reviewed above: in this method, the locus of points Pt (see Figure 6(a)) can be seen as an extra skeletal element -a curve in 2D, a surface in 3D -that is used to compute a distance surface. The procedural approach did not offer any chance of controlling the shape of this element nor the smoothness of its junction with the skeletal curves, which resulted in the problems listed above. Our basic idea is thus to define such an element explicitly. We use a triangular patch whose corners fit with the branching curves. The non-manifold skeleton made of the union of this patch and of the curves generates a junction without any bulge or crease.

Implicit Reconstruction along a Branch

The method for computing the geometric skeleton described in Section 2 outputs a graph of connected line segments, the distance to the surface being stored at each vertex (this distance is the radius of the maximal sphere associated with the vertex).

The first step for converting the geometric skeleton into a set of skeletal elements that can be used in isosurface generation is to split the graph into a set of polylines (or "branches"), connected by joints. Then, we simplify each polyline and define the field function that generates an implicit generalized cylinder around it.

Smoothing a polyline

The algorithm of Section 2 computes a skeleton whose number of vertices is of the same order as the number of initial data points. Consequently, each branch is an irregular polyline made up of a large number of very small segments (see Figure 4).

As the easiest way to generate an implicit generalized cylinder is to use the distance surfaces along a smooth curve, we approximate the polyline by a spline curve. This is done in two steps:

• We first reduce the number of vertices along the polyline using Pavlidis algorithm 21 : the polyline is first approximated by the line segment between its extremities, and then this coarse representation is progressively refined by adding some of the inner initial vertices. The algorithm to refine a segment is:

1. Compute the distance between the segment and each of the initial vertices along this branch; 2. Stop the refinement if the maximal distance is smaller than the level of detail threshold; 3. Otherwise, split the segment into two parts by adding the inner vertex which is the farther from the line segment.

As we need to take into account the variations of the radius of the cylinder as well as the geometry of its axis, we use the above algorithm in a 4D space, where the fourth coordinate of a vertex is the distance to the surface (i.e. the radius of the associated maximal sphere). • Then, we use the few remaining vertices as control points for a spline curve that approximates the initial polyline (we currently use cubic B-Spline curves).

An important remark is that the distance surface generated by a curve is smooth if and only if the radius of the cross section at each point of the curve is smaller than the local radius of curvature (see Figure 7). Otherwise, a crease will be generated. In our current implementation, we apply more smoothing to the polyline when the crease problem appears (our hypothesis is that the shape to reconstruct was smooth).

Figure 7: Constraint on the curvature of a skeletal curve: it should remain smaller than the distance to the surface. Otherwise, a crease will appear.

In future work, we plan to automatically compute the level of detail threshold to be used in step 2 of the above algorithm from the local constraint of curvature.

Field contribution along a curve

The skeleton gives precise information on the geometry of the object, since the distance to the surface is stored for each vertex. The reconstructed implicit surface should take this information into account. In P Figure 8: Implicit contribution along a curve. particular, using distance to the closest point on the curve to define the implicit contribution at point P , as in equation ( 2), no longer works (see Figure 8). Euclidean distance must be replaced by a generalized distance that takes the radius along the curve into account.

Computing the field contribution generated at a given point P by a cubic curve S(u) can be done as follows:

1. Compute the set of curve parameters (ui)i=1..p that define local minima of the distance to the point P , by finding the roots of a 5th degree univariate polynomial [START_REF] Schneider | Solving the nearest-point-on-curve problem[END_REF][START_REF] Glassner | Skeletal design of natural forms[END_REF] . 2. For i = 1..p, Compute the object thickness r(ui) from the distances to the surface at the curve's control points (we use again a cubic B-spline approximation to interpolate the radii associated with the vertex). This definition of the field function could be used for directly ray-tracing the implicit surface [START_REF] Blinn | A generalization of algebraic surface drawing[END_REF] . In practice, the images illustrating this paper have rather been rendered from a polygonalisation of the implicit surface, computed from the parametrization of the skeletal curves [START_REF] Gelston | Boundary surface recovery from skeleton curves and surfaces[END_REF] .

Figure 9 shows the conversion of the geometric skeleton of the chromosome into a graph of B-spline curves, and the surface -with creases -obtained as the union of the generalized cylinders defined along each branch.

Bulge Free Blend Between Implicit Branches

The last step of the implicit reconstruction algorithm is to generate a smooth blend between the field contributions of connected skeletal curves.

As noticed in Section 3.1, using a mere summation of field contributions gives satisfactory results when the thicknesses of connected branches are very different (see Figure 5). In that case, the blend extent can be controlled by carefully choosing the support of the function F used for each branch. However, if the branches have approximately the same size, a bulge appears near the junction. Since none of the previous approach solves this problem (see Section 3), we propose a new technique, specifically designed to model a joint between three branches (as explained in Section 2, only joints of order three are generated during the computation of the geometric skeleton). Possible extensions to higher order joints are discussed below.

Blending three branches

Our solution consist in smoothing the joint by locally replacing the three curves by a small surface patch. Then, the union of the curves and this patch, which is a smooth non-manifold entity, will be considered as a single skeletal element generating a distance surface.

In the case of three skeletal curves, we use a cubic Bézier triangle [START_REF] Farin | Curves and Surfaces for Computer Aided Geometric Design[END_REF] to smooth the joint. This triangle, defined by 10 control points, is represented by the

S(u, v, w) = u 3 P0 + v 3 P1 + w 3 P2 + 6uvwP9 + 3u 2 vP3 + 3uv 2 P4 + 3v 2 wP5 + 3vw 2 P6 + 3uw 2 P7 + 3u 2 wP8
Points P0, P1 and P2 are interpolated, while points P3 to P8 define the tangent directions (see Figure 10).

To generate a bulge/crease-free blend between three curves, we use the Bézier triangle as follows: we put P0, P1 and P2 on the curves to blend. Their position along the curve control the blend extent. We choose them according to the thickness of the solid at the joint. The same property as in Section 4.1 must be satisfied: at each point on the Bézier triangle the local radius of curvature must be larger than the surface radius, otherwise no smooth surface can be defined. be sharp to fit the curves, our first attempt consisted in merging pairs of secondary control points. These points were positioned along the tangent direction to the curve at the nearest primary control point, in order to satisfy the skeleton G 1 continuity (same tangent direction). However, when the three tangents to the curves are not coplanar, the triangle radius of curvature reaches zero near the primary control points (see Figure 11(a)); so the curvature condition (similar to Section 4.1) can never be satisfied, which causes the creases in Figure 11(b).

To avoid this problem, we use a composition of two Bézier triangles, one in the 2D parametric domain of the other one (which is in 3D):

• As before, points P0, P1 and P2 of the 3D triangle are positioned on the curves to blend. Then, pairs of adjacent secondary control points are set in the plane that passes through the closest primary control point and includes both the tangent to the curve at that point and the vector formed by the two other primary control points (see Figure 12(a) and (b)).

Since the corners of this triangle tend to be located in a plane, there is no more problem with high curvature.

• Then, we restrict the parametric domain of this patch to ensure the G 1 -continuity with the curves: we use a 2D Bézier triangle whose pairs of secondary control points are merged to map the parametric domain of the patch into a smaller domain having sharp angles at corners (see Figure 12(c)).

The resulting offset surface presents no more creases, as depicted on Figures 11(c) and 11(d).

Implicit contribution along the Bézier triangle

When the surface has a varying radius along the three skeletal curves, the distance to the implicit surface c The Eurographics Association 1997 must vary smoothly along the Bézier triangle T . We use the same approach as in Section 4.2: the equation of the Bézier triangle is used again to interpolate the distances to the surface at control points, in order to compute the desired radius r(u, v, w) at the point of barycentric coordinates u, v, w in T . Then, the field contribution along the surface can be computed in the same way than along a curve (Section 4.2), using the generalized distance: D(P, S(u, v, w)) = (d(P, S(u, v, w))r(u, v, w)) in the field computation.

Again, our implementation is based on a parametric method to display the corresponding offset surface (inspired from 10 ).

The G 1 continuity of the resulting surface is ensured by: 1. the C 1 variation of both the skeletal curves and triangles, and the offset radius along these elements. 2. the G 1 connection between these elements. (this holds also for the radius, considered as the fourth dimension of the skeleton) 3. the observance of the curvature condition.

Figure 13 shows the results for 3D curves along which the thickness is varying.

Generalizations

Although not yet implemented, extensions of this blending approach to more general branchings seems promising.

In the case of a flat branching between n curves, the method would first project the curves to a plane that is "tangent" to the joint, in order to define an order between the curves. Then, an n-sided Bézier patch 25 could be used instead of a Bézier triangle to thicken the junction.

In the general branching case, the n branches may go in any direction in space, so there is no way of flattening the junction. Then, the generalization of our method would mean embedding the joint in a Bézier volume interpolating each of the n branches (see Figure 14).

The faces of this Bézier volume would be Bézier triangles joining three of the n curves, so computing the closest point of this volume would involve closest point computations for Bézier triangles only.

Results

Figure 15 shows the final reconstruction for the chromosome example of Figure 4.

Conclusion and Future Work

This paper has developed a new approach to the reconstruction of 3D free-form surfaces, based on the combination of geometric skeletons and implicit surfaces. The main features of the method are:

• a purely geometric way to perform the reconstruction, without any optimization process to define the skeletal elements that generate the surface. • a structured reconstruction: the skeleton that generates the implicit surface is a graph of interconnected curves, so it gives an intuitive and compact representation of the object. Editing this skeleton to model deformations or articulations also appears feasible. Therefore, the method seems promising for further applications in modeling and animation.

Perfecting the reconstruction method has led us to solve a problem inherent in skeletal implicit design: the bulge free blending of surfaces generated by skeletal curves. Instead of trying to define more complex field contributions, as was done with convolution surfaces or with procedural field functions, our solution is based on the local thickening of the skeleton near a joint to prevent the creases on the corresponding distance surface. This approach offers the advantage of an explicit control on the blend extent near the joint. In a modeling system, the surface or volume element used to thicken the skeleton could be computed in an automatic way, thus being transparent to the user.

Future work includes the implementation of a smooth blend between any number of skeletal curves, using Bézier volumes. We are also looking for bounds on the curvature of Bézier triangles in order to fix their size automatically according to the local thickness of the object at a joint. This paper has studied the reconstruction of "branching shapes" for which the skeleton was a graph of interconnected curve segments. In the general case the geometric skeleton of an object includes both curves and surface elements. We are planning to use either interconnected surface patches or curve segments with anisotropic implicit contributions to model these surface elements.

Figure 1 :

 1 Figure 1: The geometric skeleton of an object.

Figure 2 :

 2 Figure 2: Relation between the Voronoi graph (a) and the geometric skeleton (b) of a set of point.

Figure 3 :

 3 Figure 3: Sensitivity of the geometric skeleton to noise in the data.

Figure 4 :

 4 Figure 4: The data set of a chromosome, its geometric skeleton (including surface elements), and the wireframe skeleton obtained at the end of the simplification process.

Figure 5 :

 5 Figure 5: Distance surface blend for a skeleton made of two segments: (a) The blend is smooth when the thicknesses are very different. (b) A bulge appears when the branches have about the same thickness. (c) Union surface produced by replacing sum by maximum in equation (1).

3 .

 3 Select the parameter u k , for k ∈ {1..p}, that minimizes the generalized distance: D(P, S(u)) = (d(P, S(u))r(u)) 4. Then, equation 2 is replaced by: fi(P ) = F (D(P, S(u k )))(3)where F is a C 1 continuous function of local support satisfying F (0) = c.

Figure 9 :

 9 Figure 9: Implicit reconstruction along branches for the chromosome example: (a) Conversion of the polylines into cubic B-spline curves of 51 control points. (b) Union surface composed of implicit surfaces of varying radius generated along each curve.

Figure 10 :Figure 11 :

 1011 Figure 10: Definition of a Bézier triangular patch: P0, P1 and P2 are the primary control points; P3, P4, P5, P6, P7 and P8 are the secondary control points.

Figure 12 :

 12 Figure 12: Definition of an extra-skeletal element.

  Figure15shows the final reconstruction for the chromosome example of Figure4. Figure16illustrates the

Figure 13 :

 13 Figure 13: Bulge free blend between three skeletal curves

Figure 14 :Figure 15 :

 1415 Figure 14: Embedding of a joint in a Bézier volume

Figure 16 :

 16 Figure 16: Reconstruction of a free form branching shape.

c The Eurographics Association 1997

Acknowledgements

Our understanding of the bulge free blend problem was really deepened by discussions with Mathieu Desbrun, Agata Opalach and Jules Bloomenthal. Many thanks to Mathieu for his help during the implementation and to Agata for carefully rereading this paper.


Eric Ferley
email: eric.ferley@imag.fr

Marie-Paule Cani Gascuel
email: marie-paule.gascuel@imag.fr

Dominique Attali
email: dominique.attali@imag.fr

Skeletal Reconstruction of Branching Shapes

We present a new method for the implicit reconstruction of branching shapes from a set of scattered data points. The method is based on the computation of a geometric skeleton inside the data set. This skeleton is simplified in order to filter noise and converted into skeletal elements -a graph of interconnected curves -that generate an implicit surface.

We use Bézier triangles as extra skeletal elements to perform bulge free blends between branches while controlling the blend extent. This leads to a smooth implicit representation of the shape, directly computed in a purely geometric way.

Introduction

With the recent developement of advanced range-imaging sensors, automatic reconstruction of real-world objects from scattered data points has become an important issue in Computer Graphics. Applications include medical imaging, inverse engineering, simulation, and semi-automatic modeling for design or animation. For the first set of applications, the aim is to provide an improved and enhanced visualization of the data. In other cases, the reconstructed object is going to be edited, deformed and animated. Here, providing a structured reconstruction is essential.

In the area of shape recognition, researchers structure a data set by computing its skeleton, defined as the locus of the centers of maximal spheres inside the data [START_REF] Blum | A transformation for extracting new descriptors of shape[END_REF][START_REF] Arcelli | Euclidean skeleton via centre-of-maximaldisc extraction[END_REF]. This skeleton, a thin centered structure that throughout this paper will be called geometric skeleton, provides a compact representation of both the topology and the geometry of the shape. Such a notion of a geometric skeleton appears attractive in the area of implicit surface modeling. The questions to be asked are: does this skeleton represent the same notion as the 1 skeletal elements used to generate an implicit surface [START_REF] Bloomenthal | Interactive techniques for implicit modeling[END_REF][START_REF] Bloomenthal | Convolution surfaces[END_REF][START_REF] Bloomenthal | Bulge elimination in implicit surface blend[END_REF]? Can we use it to construct a smooth implicit representation of an object?

Compared with previous approaches for reconstruction using implicit surfaces, the introduction of a geometric skeleton to generate the surface gives important benefits:

• The implicit surface can be computed in a purely geometric way, without the need of optimization processes such as in [START_REF] Muraki | Volumetric shape description of range data using blobby model[END_REF][START_REF] Bittar | Automatic reconstruction of unstructured 3d data: Combining medial axis and implicit surfaces[END_REF]. These processes are often computationally intensive and difficult to control.

• The skeleton, a graph of interconnected curve segments or surface elements, defines a structure for the reconstructed object. It can be used to edit the shape in an intuitive way. Such a structure is not provided when objects are reconstructed by directly computing an implicit function [START_REF] Whitaker | Algorithms for implicit deformable models[END_REF][START_REF] Velho | Approximate conversion of parametric to implicit surfaces[END_REF], or by merely generating a list of skeletal points [START_REF] Muraki | Volumetric shape description of range data using blobby model[END_REF][START_REF] Bittar | Automatic reconstruction of unstructured 3d data: Combining medial axis and implicit surfaces[END_REF].

In consequence, an approach based on skeletal reconstruction seems more adapted to modeling and animation applications.

This paper focusses on the reconstruction of free form branching shapes from the geometric skeleton of data points scattered on the object's surface. Conversion from geometric skeleton to boundary representation has already been widely studied for parametric surfaces [START_REF] Gelston | Boundary surface recovery from skeleton curves and surfaces[END_REF]. The aim of this paper is different since we want to build a smooth implicit surface representing the shape. We call "branching shapes" these objects for which the geometric skeleton can be represented as a graph of interconnected curves. Extensions to skeletons that include surface elements are discussed in the conclusion.

Section 2 reviews the definition of the geometric skeleton of an object and presents the semi-continuous approach that was used to compute it from a data set. Section 3 points out that this skeleton cannot directly be used as an input to generate an implicit surface, and lists the problems to solve. Section 4 presents a method to convert each branch of the geometric skeleton into a skeletal curve that generates a smooth implicit surface of varying radius in the cross-section. Section 5 develops a new method, based on Bézier triangular skeletal elements, to smoothly blend branches at joints without generating bulges or creases. Results are presented in Section 6. We conclude and focus on future work in Section 7.

Geometric Skeleton

Definitions

The geometric skeleton of an object is the locus of the centers of maximal spheres inside this object. A sphere included in an object is said to be maximal if it is not included in any other sphere included in the object.

The geometric skeleton is a thin structure, represented as a graph of interconnected curve segments or surface elements, centered in the object (see Figure 1). It stores in a compact way both the topological properties of the object -given by its graph structure -and the geometrical information, i.e. the distance to the surface. Note that the geometric skeleton is based on the same principles as the "medial axis" used in discrete geometry [START_REF] Lam | Thinning methodologies -a comprehensive survey[END_REF]. However, the medial axis of a data set is built as a volume made of an unstructured set of voxels.

In contrast, the graph structure of the geometric skeleton gives useful neighboring information. The next section explains how the geometric skeleton, which may initially include surface elements, is approximated by a set of polylines.

Computation of the geometric skeleton

The geometric skeleton of an object is computed using the continuous approach described in [START_REF] Brandt | Continuous skeleton computation by Voronoi diagram[END_REF]. The input is a set of points {p i } n i=1 located on the surface of the object. The output is a subgraph of the Voronoi graph of the boundary points p i .

The Voronoi graph is a well known data structure in computational geometry [START_REF] Aurenhammer | Voronoi diagrams -a survey of a fundamental geometric data structure[END_REF]. In summary, the Voronoi graph of a set of points (called seeds) divides the space into regions. Each region is the set of points closer to a particular seed than to any other seed.

In [START_REF] Brandt | Continuous skeleton computation by Voronoi diagram[END_REF], the geometric skeleton is approximated with the Voronoi elements completely included in the object (see figure 2). The result is a thin shape made up of straight-line segments in 2D and polygons in 3D (Figures 3(a) and 4(b)). Each Voronoi vertex located on the geometric skeleton is the center of a sphere passing through at least three seeds and having no seed in its interior. The radius of this sphere approximates the distance to the surface at each vertex.

Pruning branches due to noise

A drawback of the skeleton transformation is its sensibility to noise. Noise on the boundary of an object may significantly change the aspect of its skeleton (Figures 3(b)). A simplification algorithm is therefore necessary to remove peripheral branches having no perceptual relevance. Attali and Montanvert [START_REF] Attali | Semicontinuous skeletons of 2D and 3D shapes[END_REF] propose a method of simplification of continuous skeletons, that are made of both line segments and surface elements in the general case (see Figure 4). Peripheral branches and polygons are removed one after the other while they satisfy a removing criterion. By construction, the simplified skeleton is a subset of the initial skeleton that has the same class of homotopy. In 3D, the result is a wireframe figure. The order of all branching points is three.

Problems Raised by Skeletal Implicit Modelling

For people working in the area of implicit modeling, the use of a skeletonie. a set of geometric primitives admitting well-defined distance functions -to generate an implicit surface is an intuitive concept. Ideally, a smooth surface representing an object should be directly computed from any skeletal description.

Unfortunately, using the geometric skeleton described in Section 2 for isosurface generation is not trivial. This section first reviews methods for skeletal implicit modeling, and then points out the difficulties raised by the use of the geometric skeleton computed in Section 2.

Implicit surfaces generated by skeletons

Implicit surfaces have been defined as iso-surfaces of field functions 1 generated by a set of geometric primitives that define the "skeleton" of the object [START_REF] Bloomenthal | Interactive techniques for implicit modeling[END_REF][START_REF] Bloomenthal | Convolution surfaces[END_REF].

The implicit contribution f i of each skeletal primitive S i is a decreasing function of the distance to this element. The implicit surface S is defined as an iso-surface of isovalue c for the sum of skeletal contributions:

S = P/ i f i (P ) -c = 0 (1)
This formulation is a generalization of the original Blinn's objects [START_REF] Blinn | A generalization of algebraic surface drawing[END_REF], metaballs [NHK + 85] and soft objects [START_REF] Wyvill | Data structure for soft objects[END_REF], where only skeletal points were used.

If the summation is replaced by a maximum in equation (1), the implicit volume is the union of the volumes created by each skeletal primitive (see Figure 5(c)). The advantage of summation is that it defines a smooth blend between these contributions.

However, the way each f i varies affects the blend and generating a smooth shape is not easy. In particular, bulges (see figure 5(b)), characterized by a "cross section that exhibits negative, then positive, then negative curvature with respect to the underlying skeleton" 2 , should be avoided.

The next paragraphs review the different models that have been proposed for field functions.

Distance surfaces

Distance surfaces [START_REF] Bloomenthal | Interactive techniques for implicit modeling[END_REF] define f i as a decreasing C 1 continuous function F of the distance 3 to S i :

f i (P ) = F (d(P, S i )) (2)
The choice of a function F with local support limits the influence of each skeletal component, thus providing local control of the shape and reducing computations. Anisotropic field functions depending on the location of the closest point on S i may be used to generate more complex shapes [START_REF] Kacic | Controlled blending of procedural implicit surfaces[END_REF].

Distance surfaces can only be used for convex skeletons, otherwise creases will appear. Dividing non-convex skeletons into convex parts is not a solution, 1 Sometimes called "implicit function" although they are defined by an explicit equation. 2 This definition was proposed by Jules Bloomenthal in [START_REF] Bloomenthal | Bulge elimination in implicit surface blend[END_REF].

3 Defined as the distance to the closest point on Si since bulges at joints would be produced in this case. Consequently, modeling a ramification without any bulge or crease is difficult (see Figure 5).

Procedural methods

Procedural field functions were introduced in [BW90] as an example of complex implicit definition. For a specific example of two branching curves S 1 and S 2 , the two points P 1 and P 2 respectively closest to P on the two skeletal curves S i are computed, and then f i (P ) is defined by the distance to the closest point on the segment [P 1 , P 2 ]:

f i (P ) = F (d(P, [P 1 , P 2 ]))
Unfortunately, this approach that seems to offer an elegant solution to the bulge problem cannot be applied in the general case:

• It does not work for sharp angles between branches, as shown in Figure 6, since the surface tangent changes suddenly when the projected point P 1 (respectively P 2 ) is located on the extremity of the line segment.

• There is no control of the blend extent: it is restricted by the projective domains of points P on S i .

Two consecutive ramifications cannot be modeled with this method.

Convolution surfaces

Convolution surfaces [START_REF] Bloomenthal | Convolution surfaces[END_REF] define f i as the integral over S i of contributions from individual points. Here, the resulting surface does not depend on the way a skeleton is divided into skeletal elements, which is more satisfactory. Moreover, the resulting surface is smooth, even for a non-convex skeleton. Convolution surfaces solve the bulge problem if certain constraints are applied to the skeleton. Previously, these constraints have been applied to skeletons consisting of surfaces or volumes. No solution has been proposed in the

P1 Pt P2 P -4 -3 -2 -1 0 1 2 3 -4 -2 0 2 4 (a) -4 -3 -2 -1 0 1 2 3 -3 -2 -1 0 1 2 3 (b)
Figure 6: Effects of the procedural method to blend the contributions of two line segments: (a) The method works when the angle between the segments is greater than π 2 . (b) Strange features appear for sharp angles. case of one dimensional primitives: a bulge will appear near the joint in the 3D case. The solution that uses a thick volume skeleton instead of curves is not applicable in the general case, since convolution of generalized cylinders would be difficult to compute.

Problems to be solved

Applying one of the previous skeletal modeling techniques to reconstruct a smooth shape from the geometric skeleton of a data set is not straightforward:

1. Each branch of the geometric skeleton computed in Section 2 is made up of many segments. Due to noise in the data, the resulting polyline may not be very smooth. The distance surface generated by this polyline would exhibit creases, while the summation of each segment's contribution would generate bulges. The use of convolution surfaces would be computationally intensive regarding the number and size of the segments.

2. In order to create a smooth reconstructed surface, we must be able to define a smooth blend between connected branches of the skeleton. None of the previously mentioned implicit techniques gives a direct answer to this problem.

Section 4 describes a solution to the first problem, based on the conversion of skeleton branches into smooth skeletal curves.

Section 5 presents a new method to perform a bulge-free blend of the implicit surfaces generated by branching curves.

Implicit Reconstruction along a Branch

Polyline Simplification

The continuous approach of Section 2 produces a skeleton whose number of vertices is of the same order of magnitude as the number of the data points. be replaced by a generalized distance that takes the radius along the curve into account.

The implicit contribution f i (P ) of varying strength along a cubic curve S(u) can be for instance computed as in [START_REF] Bloomenthal | Skeletal design of natural forms[END_REF], by computing the zeros of a five degree polynomial, and taking the one that minimizes the distance normalized by the thickness.

In practice, our figures were generated by another method inspired by the parametric approach of [START_REF] Gelston | Boundary surface recovery from skeleton curves and surfaces[END_REF].

Figure 8 shows the conversion of the geometric skeleton of the chromosome into a graph of B-spline curves, and the union surface obtained by using each of these curves in isolation to generate an implicit surface of varying radius.

Bulge free blend between implicit branches

The last step of the implicit reconstruction algorithm is to generate a smooth blend between implicit contributions of connected skeletal curves.

Branches of very different thicknesses

As noticed in Section 3.1, the direct use of blended distance surfaces gives good results when the thicknesses of connected branches are very different (see Figure 5). In this case, the locality of the blend can be controlled by carefully choosing the support of the function F used for each branch. However, if branches have approximately the same size, a bulge will appear near the junction. As explained in Section 3, none of the previous implicit modeling techniques gives a solution to this specific problem in the case of one dimensional skeletal primitives.

The next sections present a new solution to bulge free blending, based on the local thickening of skeletons at joints. The method has already been implemented for joints between three skeletal curves. We are currently extending it to more general junctions.

Blending three branches

Our solution proposes to smooth the joint by locally replacing the three curves by a small surface patch. Then, the union of the curves and this patch, which is a smooth non-manifold entity, will be considered as a single skeletal element from which distance surfaces will be generated.

In the case of three skeletal curves, we use a cubic Bézier triangle [START_REF] Farin | Curves and Surfaces for Computer Aided Geometric Design[END_REF] to smooth the joint. This triangle, defined by 10 control points, is represented by the parametric equation: S(r, s, t) = r 3 P 0 + s 3 P 1 + t 3 P 2 + 6rstP 9 + 3r 2 sP 3 + 3rs 2 P 4 + 3s 2 tP 5 + 3st 2 P 6 + 3rt 2 P 7 + 3r 2 tP 8 Points P 0 , P 1 and P 2 are interpolated, while points P 3 to P 8 define the tangent directions (see Figure 9).

The points where the triangle is attached to the curves control the blend extent. We choose them according to the thickness of the solid at the joint. The same property as in Section 4.2 must be verified: the radius of curvature of the Bézier triangle must be larger than the surface radius, otherwise no smooth distance surface can be defined.

This restrained us from merging the tangents to the triangle at its extremities (by merging P 3 and P 8 , P 4 and P 5 and P 6 and P 7 ): that would have ensured the G 1 -continuity between the curves and the triangular patch, but then the radius of curvature would have reached zero near their connexion. So, we proceed in two steps:

• First, we position a triangular patch satisfying the curvature condition (see figure 10(a) and (b)).

• Then, we restrict the patch to ensure the G 1 -continuity between the curves and the patch. This is done by composition of Bézier triangles: a 2D triangle defines the parameter space of the previously positioned triangle (see figure 10(c)).

Implicit contribution along the Bézier triangle

When the surface has a varying radius along the three skeletal curves, the distance to the implicit surface must vary smoothly along the Bézier triangle T . We use the same approach as in Section 4.3: the equation of the Bézier triangle is used again to interpolate the distances to the surface at control points, in order to compute the desired radius r(r, s, t) at the point of barycentric coordinates r, s, t in T . Then, the implicit contribution along its surface can be computed in the same way than along a curve (section 4.3). Again, our implementation is based on an explicit method to display the corresponding offset surface (inspired from [START_REF] Gelston | Boundary surface recovery from skeleton curves and surfaces[END_REF]).

Figure 11 shows the results for 3D curves along which the thickness is varying.

Generalizations

Although not yet implemented, extensions of this blending approach to more general branchings seem promising.

In the case of a flat branching between n curves, the method would first project the curves to a plane that is "tangent" to the joint, in order to define In the general branching case, the n branches may go in any direction in space, so there is no way of flattening the junction. Then, the generalization of our method would mean embedding the joint in a Bézier volume interpolating each of the n branches (see figure 12).

Conclusion and Future Work

This paper has developed a new approach to the reconstruction of 3D free-form solids, based on the combination of geometric skeletons and implicit surfaces. The main features of the method are:

• a purely geometric way to perform the reconstruction, without any optimization process to define the skeletal elements that generate the surface.

• a structured reconstruction: the skeleton that generates the implicit surface is a graph of interconnected curves, so it gives an intuitive and compact representation of the object. Editing this skeleton to model deformations or articulations appears also feasible. Therefore, the method is potentially adapted to further applications in modeling and animation.

Perfecting the reconstruction method has led us to solve a problem inherent in skeletal implicit design: the bulge free blending of surfaces generated by skeletal curves. Instead of trying to define more complex implicit contributions, as was done with convolution surfaces or with procedural field functions, our solution is based on the local thickening of the skeleton near a joint to prevent the creases on the corresponding distance surface. This approach offers the advantage of an explicit control on the blend extent near the joint. In a modeling system, the surface or volume element used to thicken the skeleton could be computed in an automatic way, thus being transparent to the user.

Future work includes the implementation of a smooth blend between any number of skeletal curves, using Bézier volumes. We are also looking for bounds on the curvature of Bézier triangles in order to fix their size automatically according to the local thickness of the object at a joint.

This paper has studied the reconstruction of "branching shapes" for which the skeleton was a graph of interconnected curve segments. In the general case the geometric skeleton of a solid includes both curves and surface elements. We are planning to use either interconnected surface patches or curve segments with anisotropic implicit contributions to model these surface elements.

Figure 1 :

 1 Figure 1: The geometric skeleton of an object.

Figure 2 :

 2 Figure 2: Relation between the Voronoi graph (a) and the geometric skeleton (b) of a set of point.

Figure 3 :

 3 Figure 3: Sensitivity of the geometric skeleton to noise in the data.

Figure 4 :

 4 Figure 4: The data set of a chromosome, its geometric skeleton (including surface elements), and the wireframe skeleton obtained at the end of the simplification process.

Figure 5 :

 5 Figure 5: Distance surface blend for a skeleton made of two segments: (a) The blend is smooth when the thicknesses are very different. (b) A bulge appears when the branches have about the same thickness. (c) Union surface produced by replacing sum by maximum in equation (1).

PFigure 7 :

 7 Figure 7: Implicit contribution along a curve.

Figure 8 :

 8 Figure 8: Implicit reconstruction along branches for the chromosome example: (a) Conversion of the polylines into cubic B-spline curves of 51 control points. (b) Union surface composed of implicit surfaces of varying radius generated along each curve.

Figure 9 :

 9 Figure 9: Definition of a Bézier triangular patch.

Figure 10 :

 10 Figure 10: Definition of an extra-skeletal element.

Figure 11 :

 11 Figure 11: Bulge free blend between three skeletal curves

Figure 12 :

 12 Figure 12: Embedding of a joint in a Bézier volume

Figure 14 :

 14 Figure 14: Reconstruction of a free form branching shape.

Acknowledgements

Our understanding of the bulge free blend problem was really deepened by heated debates with Mathieu Desbrun, Agata Opalach and Jules Bloomenthal. Many thanks to Mathieu for his help during the implementation and to Agata for carefully rereading this paper.

Consequently, each branch of the geometric skeleton is an irregular polyline made up of a large number of very small segments (see Figure 4).

A first step before using these polylines in skeletal implicit modeling is to reduce their number of vertices, according to the local smoothness of the polyline and the desired level of detail.

We use the recursive algorithm of Pavlidis [START_REF] Pavlidis | Segmentation of plane curves[END_REF]: each polyline is first approximated by the line segment between its extremities, and then this coarse representation is recursively refined by adding some of the inner initial vertices. The algorithm to refine a segment is:

1. Compute the distance between the segment and each of the initial vertices along this branch;

2. Stop the refinement if the maximal distance is smaller than the level of detail threshold;

3. Otherwise, split the segment into two parts by adding the most distant of the initial vertices, and recursively apply the same process to the two resulting segments.

For our application, we need to take into account the variations of the object's thickness along a polyline to decide how much it should be simplified. Thus, we use the above algorithm in a 4D space, where the last coordinate of each skeleton vertex is the distance to the surface (given by the radius of the associated maximal sphere).

Conversion to a smooth curve

Polylines computed in the previous section cannot directly be used as skeletal elements generating a distance surface, since creases would appear in concavities. In order to get a smooth distance surface, we smooth the branches of the skeleton by converting the polyline into curves.

At this point, we may remark that the distance surface generated by a curve is smooth provided that the thickness remains locally smaller than the local radius of curvature (see Figure 7(a)). Under the hypothesis that we are trying to reconstruct a smooth shape, this property should always be satisfied.

Considering this remark, we approximate the polyline by a cubic B-spline curve, whose control points are the vertices of the simplified polyline. At the moment, if the curvature property is not verified, we increase the polyline simplification.

Implicit contribution along a curve

The skeleton gives precise information on the geometry of the object, since the distance to the surface is stored for each vertex. The reconstructed implicit surface should take this information into account. In particular, using distance to the closest point on the curve to define the implicit contribution at point P , as in equation (2), no longer works (see Figure 7(b)). Euclidean distance must


Nicolas Tsingos
email: nicolas.tsingos@imag.fr

Eric Bittar
email: eric.bittar@imag.fr

Marie-Paule Gascuel
email: marie-paule.gascuel@imag.fr

Implicit Surfaces for Semi-Automatic Medical Organs Reconstruction

Keywords: Implicit surfaces, medical applications, shape reconstruction

This paper proposes a new method based on implicit surfaces for medical organs reconstruction. The original data are noisy scattered points of an organ, that may be arranged in any non-uniform repartition. The knowledge of the normal vectors is not needed. Organs of any geometry and topology can be reconstructed, as for instance a vertebra that is characterized by a hole and by several branchings.

The proposed method uses implicit iso-surfaces generated by skeletons, that are a particularly compact way of giving a smooth representation of a surface. Validity of the reconstructed shape is insured, since implicit surfaces bound a well defined volume. As in a previous work [Mur91] the principle is the minimization of an energy that represents the distance between the set of points and the iso-surface. Skeletons generating field functions are automatically subdivided when needed. But the algorithm we use is novel. Firstly, we introduce local control on the reconstructed shape thanks to local and less expensive field functions. Secondly, we propose a much more efficient skeleton subdivision process, based on a notion of "skeleton energy" that gives a robust criterion for choosing the next skeleton to subdivide. Another optimization consists in splitting the data space into several overlapping windows, where only local data points are used for reconstruction. Implemented as a semi-automatic process, the reconstruction module enables the user to take benefits of his initial knowledge of the surface to guide computations. We use an interactive system for visualizing the data, initially positioning some skeletons, defining local reconstruction areas, and visualizing intermediate results. We have successfully applied the method to both non-medical and medical data.

Introduction

The problem of shape reconstruction of an object is a center of interest of many researchers in various fields, such as medical applications, computer vision, computer animation, or CAD. The medical community needs to reconstruct the human organs for diagnostics, as well as for surgery planning or simulation. The information can be contained into images, volumes, or can merely be a set of scattered points. Either coming from a range sensor or being the result of images segmentation, these points are the basis of most of the reconstruction algorithm. The large interest for the problem explains why so many different approaches have been developed: Boissonnat [START_REF] Boissonnat | Geometric structures for three-dimensional shape representation[END_REF][START_REF] Boissonnat | Three-dimensional reconstruction of complex shapes based on the Delaunay triangulation[END_REF] proposes to use the Delaunay triangulation of the set of points. Edelsbrunner and Mücke [START_REF] Edelsbrunner | Three-dimensional alpha shapes[END_REF] introduce a parameter α that enables to tune the degree of detail while deriving a shape from the Delaunay Triangulation. Attali [START_REF] Attali | Using polyballs to approximate shape and skeletons[END_REF] computes the Voronoï graph of the points to build the "skeleton" of the object and reconstruct its shape. An other approach, introduced in [HDDW92] and completed in [START_REF] Hoppe | Piecewise Smooth Surface Reconstruction[END_REF] by a mesh optimization, consists in using graph traversal techniques to calculate a signed distance function from the points, and build an iso-surface from this function.

A different philosophy consists in deforming a surface or a volume to fit the data points. The introduction of energy-minimizing curves, known as "snakes", is due to Kass, Witkin and Terzopoulos [START_REF] Kass | Snakes: Active contour models[END_REF]. The idea was then widely used and developed [START_REF] Terzopoulos | Constraints on deformable models: Recovering 3D shape and nonrigid motion[END_REF], and the active contour turned to an active surface, deformed by internal forces, such as pressure, and by external forces based on the data to fit [START_REF] Cohen | Deformable Models for 3D Medical Images using Finite Element and Balloons[END_REF]. The implementation relies most of the time on splines, but Miller also proposes a Geometrically Deformed Model based on triangles [MBE + 91]. Another option is to use deformable implicit volumes. Sclaroff and Pentland [START_REF] Pentland | Generalized implicit functions for computer graphics[END_REF] fit a super-quadric by use of modal deformations and displacement maps. Terzopoulos and Metaxas choose a physically based approach for local and global deformations of super-quadrics [START_REF] Terzopoulos | Dynamic 3D models with local and global Deformable superquadrics[END_REF].

Most of the developed methods only reconstruct surfaces of fixed topological type. However, a few papers describe models that are able to change their topology: Delingette [START_REF] Delingette | Modélisation et reconnaissance d'objets tridimensionnels 'a l'aide de maillages simplexes[END_REF] introduces the "Simplex Meshes" for which the user can interactively adapt topology. Leitner uses an adaptative spline surface [START_REF] Leitner | Segmentation dynamique d'images tridimensionnelles[END_REF], and Lachaud and Bainville a model derived from triangulation [START_REF] Lachaud | A discrete adaptative model following topological modifications of volumes[END_REF].

Muraki [START_REF] Muraki | Volumetric shape description of range data using blobby model[END_REF] is the first one who uses implicit iso-surfaces generated by skeletons for shape reconstruction. An initial skeleton is positioned at the center of mass of the data points, and divided until the process reaches a correct approximation level. But Muraki does not take a full advantage of the potentialities of implicit surfaces, and the reconstruction process he proposes is very slow. This paper presents a new method for reconstruction with implicit surfaces generated by skeletons. We introduce local control on the reconstructed shape thanks to a local field function, that enables the definition of local energy terms associated with each skeleton. This leads to a much more efficient skeleton subdivision process, since we get a robust criterion telling which skeleton should be divided next. Unlike in Muraki's work, the knowledge of the normal vectors at the data points is not needed. The method works as a semi-automatic process: the user can visualize the data, initially position some skeletons thanks to an interactive implicit surfaces editor, and further optimize the process by specifying several "reconstruction windows", that slightly overlap, and where surface reconstruction follows a local criterion. If needed, different precisions of reconstruction can be defined in each window. The shapes to reconstruct can be of any topology and geometry, and for instance include holes and branchings. We show reconstruction experiments from noisy medical data, for which scattered points are arranged in non-uniform repartition.

The remainder of this paper develops as follows: Section 2 introduces implicit surfaces generated by skeletons and discusses Muraki's reconstruction method. Section 3 details the amelioration we propose to the original algorithm, that offer both local control on the reconstructed shape and a robust heuristic for choosing the skeleton to split. The user interface we offers for semi-automatic reconstruction is described in Section 4, and results are presented in Section 5. Section 6 concludes and discusses work in progress.

2 Shape reconstruction with implicit surfaces 2.1 Implicit surfaces generated by skeletons Developed up to now as a tool for free form modeling [WMW86, WW89, BW90, BS91], implicit surfaces generated by skeletons constitute a good alternative to traditional implicit surfaces directly given by an analytical equation. Indeed, manipulating simple geometric primitives such as skeletons facilitates the design of complex shapes.

An implicit surface S generated by a set of skeletons S i (i = 1..n) with associated "field functions" f i is defined by:

S = {P ∈ IR 3 / f (P ) = iso} where f (P ) = n i=1 f i (P )
• iso is a given scalar iso-value.

• The skeletons S i can be any geometric primitive admitting a well defined distance function.

• The field functions f i are monotonically decreasing functions of the distance to the associated skeleton. They can be defined, for instance, by exponential functions [START_REF] Blinn | A generalization of algebraic surface drawing[END_REF], by pieces of polynomials [NHK + 85, WMW86], or by parametrized procedural functions [START_REF] Kacic | Controlled blending of procedural implicit surfaces[END_REF].

The surface surrounds a volume defined by f (P ) ≥ iso, and normal vectors are directed along the gradient of f .

Muraki's reconstruction method

Muraki [START_REF] Muraki | Volumetric shape description of range data using blobby model[END_REF] presents an automatic method for generating an implicit shape description from range data that include both points and normal vectors.

The surface used for reconstruction is defined by Blinn's "blobby model" [START_REF] Blinn | A generalization of algebraic surface drawing[END_REF]. Only point-skeletons are used, and fields are expressed as exponential functions of the distance d(P, S i ):

f i (P ) = b i e -aid(P,Si)
where b i can be either positive of negative1 .

The principle of the reconstruction algorithm consists in minimizing an energy that represents the "distance" between the current implicit surface and the data. Muraki proposes the use of three terms:

• E value fits the surface to the data points:

E value = M j=1 (f (P j ) -iso) 2
• E normal uses normal vectors to insure that the external side of the isosurface corresponds to the external side of the object:

E normal = M j=1 |nj -N (P i)| ||N (P i)|| 2
• E shrink expresses constraints on the relative values of the parameters a i and b i of the n skeletons that avoid surface degeneration in low-sampled areas:

E shrink = ( n i=1 a i -3/2 |bi|) 2
The total energy E is given by the expression:

E = 1 M (E value + α E normal ) + β E shrink
where M is the number of data points.

Directly using n skeletons to reconstruct the surface would be a very computer intensive process, since it would require solving a non-linear minimization problem with 5n unknowns (the 5 parameters a i , b i , x i , y i , z i for each skeleton). Muraki performs instead a progressive refinement process. The main steps of the algorithm are:

1. Position a single skeleton at the barycenter of the data points.

Make trials, independently for each skeleton S i :

• Split S i into two skeletons S ′ i and S" i : • Minimize energy by modifying only S ′ i and S" i 's parameters ; • Store the final energy.

3. Keep the modification of step 2 that gave the lowest energy (all the other trials are not used anymore).

4. If the energy is sufficiently low, stop the process ; else, go back to step 2.

Discussion

Using implicit surfaces generated by skeletons for reconstructing smooth shapes seems to be a very good idea, since these surfaces are particularly convenient for defining an energy that gives a current "distance" to the data, and since they can be stored in a very compact way (position and field parameters for each skeleton).

However, several problems can be identified, that seem to forbid the direct use of Muraki's method for our goal of medical organs reconstruction:

• The choice of an exponential field function doesn't give any locality on the reconstructed shape. Every data point is influenced by all the skeletons, even those that are distant. In consequence, each refinement changes the shape everywhere, even in areas where the reconstruction was already sufficient.

• The initialization is performed with a single skeleton, positioned at the barycenter of the data-points. This gives quite good results for reconstructing heads (the only examples given in the paper), for which the progressive refinement of an implicit sphere is well adapted. But one can wonder how the process would work for arbitrary shapes, with holes and branchings, such as a vertebra. In particular, the initial skeleton could be located inside a hole, probably perturbating the reconstruction.

• The very important computational time is really limitative: the given examples, that are said to take "a few days", correspond to 2893 data points. Direct use of the method on the echography of an organ (about 25000 scattered points) would be intractable.

The most expensive process seems to be the way the next skeleton to split is selected. The author admits that trying every skeleton for just keeping the solution that decreases the energy the most is very slow, and should be changed. In practice, some of the examples in [START_REF] Muraki | Volumetric shape description of range data using blobby model[END_REF] are computed by merely splitting each skeleton successively, in the arbitrary order given by the current list of skeletons. As this method doesn't seem very satisfactory as well, better heuristics need to be found.

A novel method for reconstruction

This paper proposes a more efficient method for shape reconstruction with implicit surfaces, that we are experimenting on medical data. This section describes our main contributions for improving Muraki's technique. The main idea we are developing is to take benefits of locally defined field functions to optimize the reconstruction process. This is done by defining a new energy criterion telling which skeleton must be splitted next, and by offering the possibility of processing reconstruction in a more local way, through the use of several "reconstruction windows".

Choice of a local field function

Local field function, that become zero with their derivatives at a certain radius of influence R have been introduced in the literature for optimizing field computations during the design of complex objects [START_REF] Wyvill | Data structure for soft objects[END_REF]. Another benefits is that these fields offer local control on the implicit surface, which is particularly important for our reconstruction process: splitting and optimizing skeletons in an area should not modify a good reconstruction of the data already obtained in another area. Contrary to Muraki that uses both positive and negative field functions, we want to reconstruct objects with positive fields only. Indeed, for medical applications that may include a future physically-based simulation of organs deformations, the exclusive use of positive skeletons seems more appropriate.

The chosen field function will be used in a computer-intensive minimization process. So we need a model that is both controlled by very few parameters (in order to limit the dimension of the search space), and that provides efficient computations. In consequence, we are currently experimenting field contributions f i that are composed of a linear and a quadratic curve segments, and that are parametrized by two parameters (see Figure 1):

• The radius e i of the sphere created by a skeleton S i alone (e i is characterized by f i (e i ) = iso) ;

• The stiffness k i in e i , that defines the blending properties of the surface.

The function is expressed by:

f i (P ) = -k i r + k i e i + 1 if r ∈ [0, e i ] f i (P ) = k i e i (e i -R i ) + 3e i -R i (e i -R i ) 3 (r -R i ) 2 if r ∈ [e i , R i ]
f i (P ) = 0 elsewhere Where r = d(P, S i ) and R i = (e i -2 ki ) is the radius of influence. As in Muraki's work, there are 5 parameters to optimize for each skeleton: e i , k i , x i , y i , z i .

Basics for semi-automatic reconstruction

The idea is to keep the main principles described in Section 2.2, but with several modifications.

First of all, most medical data not include any description of normal vectors, that were used by Muraki to insure that the internal volume was reconstructed. However, we are defining a semi-automatic reconstruction process, where the user takes benefits of his knowledge by interactively positioning the initial skeletons with respect to the data. If initial skeletons are positioned inside the volume to reconstruct, the method will work without the need of normal vectors.

In consequence, we redefine the energy to minimize by:

E = 1 M   M j=1 (f (P j ) -iso) 2   + 1 n α 1 n i=1 e -β1ei + α 2 n i=1 e -β2ki
The two last terms are basically equivalent to E shrink in Muraki's method: they prevent the skeletons from degenerating to negative radius or stiffness (n is the current number of skeletons).

As stressed in Section 2.2, the main problem with Muraki's method is the very intensive computation process, due to the method used for selecting the skeleton to divide next: every solution is tested (which requires running the full minimization) and only the best one is used. The next section describes the heuristic we have developed for accelerating this process.

An efficient subdivision algorithm

Selecting the skeleton to divide can be done without any extra minimization, by using the new local properties of the field functions. Indeed, the best candidate for splitting is the skeleton whose area of influence corresponds to the zone of the surface were the reconstruction is the worst. See Figure 2.

C i = 1 m i   mi j=1 (f (P j ) -iso) 2  
where the sum is made only on the m i point that are inside S i 's sphere of influence.

At each iteration of the algorithm, the skeleton with greatest C i is replaced by two skeletons, whose 10 parameters are first optimized while keeping the skeletons and fields constant. Then, we process a few steps of global optimization before going to the next skeleton subdivision (this global optimization step reduces the final number of skeletons: in some cases, just adjusting all the current skeletons at the same time all is sufficient for decreasing the energy to the desired threshold).

Reconstruction windows

For large scattered points databases, using all the points for every step of energy minimization can be really intractable, since the energy includes a sum of field terms over the points. On an other hand, many objects to reconstruct are easy to decompose into several specific zones, called "reconstruction windows", where reconstruction can be processed separately. For instance, 3 different windows can be defined for the object in Figure 3. The energy to minimize in a 3D window will be calculated by only using the data-points located in this window, and the global optimization step after a refinement will only be performed on the skeletons affected to this window. Processing local reconstructions is justified by our use of local field functions: when we move a skeleton, only local points are affected, so recomputing the entire energy is not really useful. Local windows optimize computations since less points are considered in energy terms, and minimization is performed in a space of smaller dimension. The algorithm we use for reconstructing data when several reconstruction windows are defined develops as follows. At each reconstruction step:

1. Select the window were:

W k = 1/w k ( w k j=1 (f (P j ) -iso) 2
is the highest (w k is the number of points in this window).

2. Use the criterion of Section 3.3 for selecting a skeleton to divide in this window.

3. Split the skeleton and optimize its two "children", everything else being kept constant. The energy to minimize is computed from the w k points and the n k skeletons of window:

E k = 1 m k   m k j=1 (f (P j ) -iso) 2   + 1 n k α 1 n k i=1 e -β1ei + α 2 n k i=1 e -β2ki
4. Process a few steps of optimization in the window in order to further decrease E k (and W k ).

5. If the global value defined by: W global = W i is smaller than the desired threshold, stop the process. Else, go back to step 1.

In practice, the list of skeletons and points included into each window are maintained to accelerate the computations. The user is provided with an interactive interface for visualizing the data, defining reconstruction windows and positioning some initial skeletons. We describe it in the next section.

Interface for Semi-Automatic reconstruction

The semi-automatic reconstruction process we have presented is implemented within the system F abule, a platform for interactive modeling and animation that we are developing at iMAGIS [START_REF] Gascuel | Fabule : un environnement de recherche pour l'animation et la simulation[END_REF]. It includes an object-oriented library, implemented in C++ and built on Open-Inventor, and an interface with the interpreted language Tcl and OSF-Motif for the easy specification of demonstrators and interfaces.

The user visualizes the data and creates reconstruction windows thanks to 3D widgets provided by Inventor. He also interactively specifies skeletons and associated field functions that will initialize reconstruction. Intermediate reconstruction results are visualized within the same interactive system. We offer an interactive graphic interface for modeling skeletons, visualizing implicit surfaces, and editing the field function parameters (see Figure 4). Implicit surfaces are visualized in real time thanks to a new sampling technique on seeds-migration, that adapts at interactive rates to any surface deformation. Sample points are not directly used for display: the user can choose between displaying "scales" laying in each sample point tangent plane, or a piece-wise polygonization of the surface, where a polygon set is associated with each skeleton (this display mode has been used in Figure 8). The full system for interactive modeling and visualization of implicit surfaces is described in [START_REF] Tsingos | Un modeleur interactif d'objets définis par des surfaces implicites[END_REF].

Experimental results

Validation of the algorithm on implicitly-defined objects

We have first validated the reconstruction algorithm on objets that had initially been created with implicit surfaces.

The object depicted in Figure 5 has been reconstructed with a single initial skeleton (thus validating the automatic subdivision process), and a single reconstruction window, including all the scene. Results show that the original skeletons and fields are correctly reconstructed by the process, without the creation of any extra skeleton.

The use of reconstruction windows

The example depicted in Figure 6 illustrates the use of reconstruction windows on an object with hole and branches, also created with implicit surfaces. We have performed two different reconstructions based on the same initial skeleton: one with a single large window surrounding the whole scene (a), and the other with several reconstruction windows (c). Figure 6 (a) shows the first result: an approximation with 16 skeletons obtained after more than two hours of computations on an Indigo 2, R4400, 150 MHz workstation.

Reconstruction of medical organs

A kidney

The difficulty for reconstructing this kidney was due to a very noisy data-base, obtained from an echography. Figure 7 shows the data-points, arranged in a non-uniform repartition, and the reconstructed shape we get.

A vertebra (work in progress)

A vertebra is a challenging example for validating the reconstruction process, due to it its very complicated geometry (a central hole and several small branchings).

Conclusion

We have presented a semi-automatic reconstruction method, that can be used on noisy scattered points of a medical organ. The method is based on implicit iso-surfaces generated by skeletons, that provides a smooth and compact representation of the surface. The basics is the minimization of an energy, that represents the distance between the data and the surface, during a series of skeleton refinements. We propose an efficient automatic way of selecting the skeleton to subdivide next, and a new method for optimizing reconstruction through the use of "reconstruction windows". The method is validated on both medical and non-medical data. The user can guide the reconstruction by initializing some skeletons and their reconstruction windows, thus taking benefits of his initial knowledge of the data.

An easy extension would be to enable the use of other kinds of skeletons generating field functions (such as segments) during the reconstruction process.

For some applications, an automatic generation of the initial skeletons may be more convenient that the interactive initialization process we have described. We are currently studying an approach based on "medial axis" for the automatic positioning of initial skeletons and reconstruction windows.

Figure 1 :

 1 Figure 1: Local field function with two parameters.

Figure 2 :

 2 Figure 2: Selection of the next skeleton to subdivide So we define the splitting criteria as:

Figure 3 :

 3 Figure 3: Defining several windows for local shape reconstruction

  Figure 4: Interface for editing skeletons and field functions

Figure 5 :

 5 Figure 5: Reconstruction of an iso-surface initially created with 3 punctual skeletons

  Figure 6 (b), computed much more efficiently thanks to the local windows, is an approximation with 19 skeletons obtained in half an hour.

  (a) Reconstruction with a single large (about 2 hours) (b) Interactive positioning of 6 local windows in the data space.(c) Reconstruction with local reconstruction windows (about 30 mn).

Figure 6 :

 6 Figure 6: The use of several windows for accelerating reconstruction

Figure 7 :

 7 Figure 7: Reconstruction of a kidney from a very noisy database

Using a "negative skeleton" ie. a skeleton with a negative field reduces the implicit volume, while positive skeletons extend it.