Convection-Driven Dynamic Surface Reconstruction

Remi Allegre Raphaelle Chaine Samir Akkouche

LIRIS CNRS, Universite Claude Bernard Lyon 1, France

{Remi.Allegre, Raphaelle.Chaine, Samir.Akkouche }@liris.cnrs.fr

Reconstruction of a simplified mesh Local dynamic reconstruction refinement

Figure 1. Our reconstruction framework, illustrated on the TRIPLE HECATE model. Starting from a dense input point set, we reconstruct a simplifiedmesh (center). Benefiting from the connectivity of this initial reconstruction, we can make it to evolve dynamically so as to refine the approximationlocally. This refinement can be achieved either in an automatic fashion, for example in order to improve the quality of the elements of the mesh, orinteractively, in order to add or remove sample points. Here, the draped dress has been locally enhanced (right).

Abstract

In this paper, we introduce a flexible framework for thereconstruction of a surface from an unorganized point set,extending the geometric convection approach introduced byChaine [9]. Given a dense input point cloud, we first ex-tract a triangulated surface that interpolates a subset of theinitial data. We compute this surface in an output sensitivemanner by decimating the input point set on-the-fly duringthe reconstruction process. Our simplification procedurerelies on a simple criterion that locally detects and reducesoversampling. If needed, we then operate in a dynamic fash-ion for local refinement or further simplification of the re-constructed surface. Our method allows to locally updatethe reconstructed surface by inserting or removing sam-ple points without restarting the convection process fromscratch. This iterative correction process can be controlledinteractively by the user or automatized given some specificlocal sampling constraints.

Keywords: surface reconstruction, point set simplification,geometric convection, dynamic correction.

1. Introduction

Shape modeling from point-sampled geometry has re-ceived considerable attention in the past few years, due tothe recent advances of 3D digital acquisition and the in-creasing number of application domains. Today’s rangescanning devices are able to produce highly detailed digitalsurface models that can contain millions of sample points.To produce efficient shape representations for interactive vi-sualization or further geometry processing, the complexityof these models has to be reduced.

From a dense input point set, we consider the problem ofcomputing a simplified piecewise linear surface. This goalcan be achieved by first reconstructing an initial mesh, andthen simplifying this mesh. However, the time and mem-ory costs can be prohibitive if connectivity relations haveto be established for all points of the input point set. An-other solution consists in first simplifying the input pointset, and then reconstructing. The goal of point set simpli-fication algorithms is to extract the relevant data of a denseinput point set so as to accelerate a subsequent surface re-construction process. Most of these algorithms are not de-signed to perform both point set simplification and surface

reconstruction in a single stage, and do not compute theirresult in an output sensitive manner. The coarse-to-fine sur-face reconstruction algorithm by Boissonnat and Cazals [5]is a notable exception. One limitation of this work is that therefinement is not achieved in a complete dynamic fashionwhen a consistent orientation for normals has to be deter-mined. For the purpose of multiresolution shape modeling,when transmitting point samples over networks, or duringthe 3D acquisition process, updating the reconstructed sur-face on-the-fly by incorporating or removing points dynam-ically can be useful.

In this paper, we introduce a dynamic framework for thereconstruction of a simplified triangulated surface from aset of unorganized points sampled from a smooth surfacein R3. We build upon Chaine’s geometric convection al-gorithm [9]. As many surface reconstruction algorithms inthe Computational Geometry community [8], this algorithmoutputs a triangulated surface embedded in the 3D Delau-nay triangulation (or DT for short) of the input point set.Unlike in the original method, we do not require a global3D DT of the entire input point set. Given a geometric errortolerance, we produce a piecewise linear approximation ofthe original surface that interpolates only a relevant subsetof the input data. We compute this surface in an output sen-sitive manner by decimating the input point set on-the-flyduring the reconstruction process. Our simplification pro-cedure involves a simple criterion that locally detects andreduces oversampling. If desired, the reconstructed surfacecan be customized in a second stage either to refine the ap-proximation, or to eliminate undesirable features. One keyfeature of our method is the ability to dynamically insert andremove sample points without restarting the reconstructionprocess from scratch, taking benefit from the current recon-structed surface. This iterative process can be automatizedgiven some specific local sampling constraints or controlledinteractively.

1.1. Related Work

Our work is closely related to surface reconstruction andsurface resampling problems. We focus here on point setsimplification techniques and sampling algorithms that out-put a mesh approximation.

There are mainly two kinds of algorithms for simplify-ing a dense point set:subsamplingand resamplingalgo-rithms. Subsampling algorithms output a decimated pointset that is a subset of the original point set. This can beachieved fine-to-coarse by iterative point removal opera-tions [12, 15, 1, 19], or in a coarse-to-fine fashion, usinga Farthest Point Sampling Strategy [16]. Boissonnat andCazals [5] have proposed a coarse-to-fine reconstruction al-gorithm controlled by a signed distance function to the re-constructed surface. This signed distance function is com-

puted from natural neighbor interpolation of a random sub-set of the input point set with oriented normals. This subsetis iteratively enriched till the result fits a geometric error.

Resampling algorithms rely on a global or local esti-mated representation of the true surface to compute new,well chosen point locations. This representation is com-puted from the input point set [13, 18, 17, 19]. Boissonnatand Oudot [6, 7] have recently revisited Chew’s FarthestPoint sampling technique [10] to generate optimal featuresize-dependent triangulations for a fixed implicit or polyhe-dral surface. The authors reconstruct a mesh surface coarse-to-fine from a point set using a Moving Least Squares im-plicit surface representation.

1.2. Overview

Let P = {pi} be a set of sample points that lie on ornear a smooth surfaceS embedded inR3. We suppose thatP is sufficiently dense in the sense this point set forms aε-sample ofS for some constantε > 0 [2]. This point setcan be locally oversampled w.r.t. the local feature size, i.e.the shortest distance to the medial axis. In our framework,we need an (unoriented) normal direction at every samplepoint p ∈ P . If normals are not supplied as part of theinput data, we estimate the normal direction at a pointp byfitting a least squares plane top and itsk nearest neighborsin a preprocessing step. For a reliable estimation, a locallyuniform sampling distribution is required [3]. We recall thegeometric convection algorithm in Section 2. We proceedin two stages, as illustrated in Figure1.• In the first stage (Section 3), we compute a linear ap-

proximation ofS that interpolates a subsetP ′ of P w.r.t.a geometric error toleranceρ > 0 prescribed by the user.Each time a new point is inserted in the reconstructed sur-face, we decimate the input point set in a small neighbor-hood around. We achieve this simplification in a feature-sensitive manner thanks to a normal-based error metric. Theresult is a consistent triangulated surface embedded in the3D DT of P ′. For this reconstruction, we do not requireto compute the 3D DT ofP explicitly. We only need tocompute the 3D DT ofP ′ in prevision of the second stage.• In the second stage (Section 4), corrections can be ap-

plied dynamically to the reconstructed surface. We proposea refinement algorithm to improve the quality of the trian-gles and give the possibility to the user to customize theresult by adding or removing details. For this purpose, weintroduce an algorithm to locally update the reconstructedsurface by inserting or removing sample points. We storethe history of the reconstruction process by maintaining the3D DT of the set of points that actually belongs to the re-constructed surface.

In Section 5, we present some experimental results. Weconclude and discuss future work in Section 6.

(a) (b) (c) (d)

Figure 2. Convection towards a 2D point set. In (a), the pseudo-curve is initialized to the convex hull of the point set. The current edge, enclosedby its (dashed) nonempty half-circle, forms a Delaunay triangle (dark gray) with the dark gray point. This triangle becomes external, and the pseudo-curve is updated consistently in (b), where two other Delaunay triangles have been opened. The pseudo-curve evolves as long as the oriented Gabrielproperty is not met for any of its oriented edges. In (c), an edge is found to block a cavity. The result in (d) contains coupled oriented edges (the’tail’ at the bottom right of the shape).

2. The geometric convection algorithm

The geometric convection algorithm [9] is based on theconvection model introduced by Zhao, Osher and Fed-kiw [20]. From a point setP sampled from a surfaceS,the latter solve the reconstruction problem by computing aclosed surface that minimizes a global distance function tothe input point set. A convection scheme is used to com-pute an initial approximation ofS. This approximation isobtained by shrinking a surfaceS′ that enclosesP . At eachstep, every pointx of S′ evolves along the normal directionn(x) of S′ at pointx, with displacement speed proportionalto−∇d(x) · n(x) whered(x) is the distance betweenx andits closest point inP .

Zhao et al. [20] compute the convection result on a regu-lar grid with a so-called fast tagging algorithm. Chaine [9]translates the convection scheme into Computational Ge-ometry terms, yielding an efficient, purely data-dependentsurface reconstruction algorithm. This work is built on thefollowing theorem.

Theorem (proved in Chaine [9]) Given a closed surfaceS′ enclosing a point setP , the convection ofS′ throughthe velocity field−∇d(x) converges to a closed, piecewiselinear pseudo-surface. All the facets of this pseudo-surfaceare Delaunay facets oriented consistently towards theinterior of the shape that meet anoriented Gabriel property.

An oriented Delaunay facetf is said to meet theori-ented Gabriel propertyif the half of its minimum enclosingsphere located on the positive side off does not containany point ofP in its interior. The termpseudo-surfacemeans that different parts of this surface can be pinchedtogether, i.e. can locally share common geometric infor-mation, while remaining topologically independent. Moreformally, a piecewise-linear pseudo-surface can be definedas the geometric embedding of an orientable manifoldpolyhedral complex such that the geometric images of twovertices, edges of facets are either identical or disjoint. Thisstructure can represent the evolving surface all along the

reconstruction process. It supports topological changes andallows to maintain a consistent mesh at every step of theconvection scheme.

Algorithm From the above theorem is derived an algorithmthat extracts the reconstructed surfaceS′ from the 3D DT ofP . The idea is to shrink an initial piecewise-linear pseudo-surface (or just pseudo-surface for short) through this trian-gulation, which is equivalent to make it convect through thevelocity field−∇d(x). Figure2 illustrates the reconstruc-tion process on a 2D point set.

The pseudo-surface is initialized with the boundary ofthe convex hull ofP , all facets oriented inwards. If an ori-ented facet does not meet the oriented Gabriel property, itis removed from the pseudo-surface to be replaced by threeother ’hidden’, consistently oriented facets of the Delaunaytetrahedron it belongs to. An oriented facet can be openedtowards a point location that is already attached to a vertexof the current pseudo-surface. Two oriented facets of thepseudo-surface with identical geometry are said to becou-pled – they necessarily have opposite orientations. Whenone of two coupled oriented facets of the pseudo-surfacedoes not meet the oriented Gabriel criterion, both are re-moved and are said tocollapse, which can involve topolog-ical changes of the pseudo-surface. All facets are processedseparately, following a breadth-first traversal.

The convection stops as soon as every oriented facetmeets the oriented Gabriel property. In presence of con-cavities larger than the Gabriel half-spheres raised from thepseudo-surface, the basic convection process stops prema-turely. If the point sample correctly reflects the local fea-ture size, the pseudo-surface is further shrunk through thesepocketsby detecting inconsistencies between the size ofblocking facets and local density. A more global solutionfor this problem can be derived from a topological persis-tence criterion [14, 8].

The resulting pseudo-surface is a combinatorial mani-fold, that can contain coupled oriented facets calledthinparts. These thin parts are not always significant, so it isimportant to identify desirable ones. This is achieved by

pursuing the convection process in 2D on these thin parts,starting from their boundary.

Complexity The complexity of this algorithm is dominatedby the time to compute the 3D DT of the input point set. Itcan be computed incrementally in timeO(N log N), whereN is the size of the input point set [11]. For a locally uni-form ε-sample, the number of tetrahedra is almost linearw.r.t. N [4]. Remaining operations correspond to a partialbreadth-first traversal of the 3D DT.

The basic convection process does not require anyglobalDelaunay-related information such as the poles [2], that canbe needed to obtain an estimation of the local feature sizefor example. Only the pocket detection could benefit fromsuch a global information, but the proposed local solutiongives satisfactory results in practice. The 3D DT allows todirectly identify the points towards which oriented facetsopen during the convection process. For redundant inputpoint sets, time and memory are wasted unnecessarily, sincenot all points are relevant. For our purpose of reconstructinga simplified triangulated surface, we investigate a way toavoid the computation of the 3D DT of the entire point set.

3. Reconstruction and simplification

The first stage of our framework consists in constructinga mesh from the input point cloud while simplifying it. Tofit our simplification purpose, we couple the geometric re-construction algorithm with a subsampling procedure. Wetake benefit of the locality property of the convection pro-cess to avoid the computation of the 3D DT ofP . Next,we describe our decimation scheme. Then, we present thecomplete reconstruction algorithm and take a look at someproperties of the reconstructed surface.

3.1. Decimation scheme

Let Pev denote the set of points interpolated by theevolving pseudo-surface at a given time. If a new sam-ple pointp is inserted into the evolving pseudo-surface, weremove some redundant points ofP − Pev located in itssurface neighborhood following an isotropic region grow-ing strategy. We grow the neighboring region aroundp byadding its nearest neighbors in the order of their distanceto p. The idea is to makep a good representative of thegeometric information hold by the sample points to be re-moved in its neighborhood, given an error toleranceρ ≤ 1.A good representativemeans that the set of facets that willbe incident to this point will have to correctly approximatethe surface both geometrically and topologically, at globalscale. We thus define two kinds of sub-sampling criteriato control the decimation process: onegeometric approxi-mationcriterionCgeom, and onetopological diskcriterion

Ctopo. A point that fulfills both criteriaCgeom andCtopo iseliminated, and the growing stops as soon as one of theseconditions is not met. The distance betweenp and the far-thest point that fulfills bothCgeom andCtopo will be calleddecimation radiusof p.

Our geometric criterion (Fig.3) relies on a normal-basederror metric. A pointpi fulfills the geometric approxima-tion criterion if and only if the following condition is met:

(Cgeom) |n(pi) · n(p)| > ρ

where0 < ρ ≤ 1.The topological criterion must guarantee that the ball

B(p) bounding the decimation region aroundp con-tains only sample points that belong to the local surfaceneighborhood ofp on S. This criterion is designed toprevent oriented facets from having an half-sphere that’encroaches’ another part of the surface, which would resultin a topologically incorrect reconstruction. As illustratedin Figure 4, we derive an intrinsic criterion. A pointpi

fulfills the topological criterion if and only if the followingcondition is met:

(Ctopo) |n(pi) ·p− pi

‖p− pi‖| < ρ′

where0 < ρ′ ≤ 1 is a value that depends on the samplingconditions. IfP = S, the idea is that ifB(p) touches another part ofS at a sample pointpi, B(p) becomes tangentto S at pi so that n(pi) andppi are collinear. The trans-position to the discrete setting is straightforward, providedP reflects the local feature size. Depending on the inputtype, theoretical bounds forρ′ could be derived. In prac-tice, we found thatρ′ = 0.95 generally gives satisfactoryresults. Figure5 illustrates the effect of this criterion on thescrewdriver model.

pi

θmax

pi

pn(p)

p

Figure 3. Decimation of the point set in the neighborhood ofp.A point pi with a normal that makes an angle less thanθmax withn(p) satifies the conditionCgeom.

3.2. Reconstruction

We extend the original convection algorithm by in-troducing the previously described decimation scheme toproduce a simplified reconstruction. Our algorithm doesnot compute the 3D DT of the input point set explicitly.We shrink a pseudo-surface mesh aroundP by computingrequired Delaunay tetrahedra on-the-fly, while eliminating

(a) (b) (c) (d)

Figure 6. Convection with simplification towards a 2D point set. In (a), the pseudo-curve is first initialized to the convex hull of the point set. Theblack points have been eliminated while computing the convex hull. The current edge forms a Delaunay triangle with the dark gray point. In (b), thetwo black points in the neighborhood of the selected point are found to lie below the prescribed error tolerance. These points are removed and thepseudo-curve is updated in (c). The pseudo-curve evolves as long as the oriented Gabriel property is not met for any of its oriented edges. The finalresult is shown in (d).

pi

n(p )i

n(p)

p

Figure 4. Illustration of the topological disk criterion. Here,piis the first point that does not belong to the topological disk con-structed fromp. The angle between the normal n(pi) and the edgeppi should be small ifP is sufficiently dense.

Figure 5. Reconstruction of the SCREWDRIVERmodel with sim-plification (k = 6, ρ = 0.98). The reconstruction on the left doesnot uses our topological criterion. On the right, the topological cri-terion is enabled.

irrelevant sample points. We will denote asR the set ofpoints eliminated during the reconstruction process, i.e.R = P − P ′. In the two following paragraphs, we presentour data-structures and detail our algorithm.

Data-structures We represent the evolving surface by apseudo-surface mesh data-structure [9] that is not supportedby an explicit 3D DT of the input point set. We delegatespatial searching to akd-tree data-structure that stores theinput point setP . We require this data-structure to performpoint locations efficiently.

Algorithm We first initialize the pseudo-surfaceSev bycomputing an approximate convex hull ofP using the

Figure 7. Two reconstructions of the APHRODITE model (left,46K points) with different values ofρ. With ρ = 0.9, the surface isapproximated by 3.5K points (center). Withρ = 0.98, the surfaceis approximated by 10.8K points (right).

Quick Hull algorithm. Each time a new sample point isadded to the convex hull, we call our decimation procedureon this sample point. We then start from this mesh for theconvection process.

Let pqr be an oriented facet of the pseudo-surface. Todecide whether this facet should be opened or not towardsa point ofP , it suffices to check whether its meets the ori-ented Gabriel criterion. We first report all points ofP lo-cated in the Gabriel half-sphereH of pqr .• If the criterion is not verified, we check whether the twocoupled oriented facetspqr and qpr both belong toSev.If it is the case, the two coupled oriented facets collapseand the connectivity is restored between their neighboringfacets [9]. Otherwise, we have to find the points∈ H∩(P−R) such thatpqrs forms a Delaunay tetrahedron. This pointis the one that maximizes the radius of the circumsphereS of pqrs. The sphereS is a medial sphere that satisfiesthe Delaunay empty ball criterion. We call the decimationprocedure ons if s does not belong to the pseudo-surface.Then, we just replacepqr by the oriented trianglespqs, qrsandrps, which can be achieved by a simple vertex insertion.• If pqr satisfies the oriented Gabriel property, we rely onnormals to detect and pursue the convection process throughpockets. Let nf denote the unit normal to the oriented facet

f . If the greatest value between|nf · n(p)|, |nf · n(q)|,and|nf · n(r)| is less than0.5, we consider thatf blocks apocket. In this case, we search for the points such that thetetrahedronpqrs forms a Delaunay tetrahedron. We call thedecimation procedure ons if sdoes not already belong to thepseudo-surface, and then inserts into pqr as previously.

The convection stops as soon as every oriented facet ofthe pseudo-surface meets the oriented Gabriel property anddoes not pass the pocket detection test. Our reconstructionalgorithm is schematized in 2D in Figure6. Figure7 il-lustrates our technique on the APHRODITE point set withdifferent values of the tolerance parameterρ.

3.3. Sampling density and regularity of the mesh

The sampling density and regularity of the resultingmesh are controlled both by the radius of the decimationregion around sample points and by the order in which theconvection processes them. We provide here some elementsof explanation, but not a complete study.

Our simplification procedure is closely related to the ge-ometry of the whole surface, but not directly guided by thelocal feature size [2]. The radius of decimation is anothergeometric measure that quantifies the local thickness of thesurface. The growth of a decimation ball around a samplepoint in regions with important curvature changes is stoppedby theCgeom criterion, that only depends on the local sur-face variations. The radius of a decimation ball built arounda sample point on a less curved part is determined byCtopo

for a sufficiently permissive criterionCgeom, i.e. by the dis-tance to the opposite side of the surface.

The final sampling distribution can be explained by theorder in which sample points are incorporated into theevolving pseudo-surface. An oriented facet that does notmatch the oriented Gabriel property tends to open towardsa point that forms a neighborhood edge with one of its threevertices, except in the case where the entered Delaunay cellcorresponds to a branching of the medial axis. The neigh-borhood graph on the pseudo-surface thus results from abreadth-first propagation around existing vertices that is in-duced by the geometric convection process.

This way of processing yields regular connectivity inregions that exhibit no important change in curvature(Fig. 8(a)). In smooth regions, it can also be observedthat retained points are distributed uniformly around a givenone. However, strong curvature variations induced, forexample, by the presence of sharp features such as edgesor corners, produce long and skinny triangles coupled withhigh valence vertices, which can be inadequate for furthermesh processing (Fig.8(b)). When growing the decima-tion region around a sample point, we have no cheap meansto anticipate on these variations. Sharp features could bedetected earlier by precomputing the radius of the decima-

tion region for every input sample point without simplify-ing, but we prefer to compute this radius only for retainedsample points. We improve the quality of the mesh that re-sults from the reconstruction process in a second stage, thatwe develop in the following Section.

p

(a) (b)

p

Figure 8. Optimal (a) vs. non-optimal (b) neighborhood configu-ration. Every sample point is represented with its circular decima-tion region.

4. Dynamic correction

The purpose of the correction stage is to refine the ini-tial reconstructed surface by inserting or removing samplepoints. This functionality first serves the purpose of improv-ing the quality of bad shaped triangles. Second, we want toprovide the user with means of locally changing the levelof detail of the reconstruction. Two kinds of questions nat-urally arise: Which points are good candidates for furtherinsertions or deletions? How to update the pseudo-surfacein a dynamic fashion after these points have been chosen?

For triangle-quality improvement, we devise a simplegreedy refinement algorithm inspired from Chew’s algo-rithm [10], recently revisited by Boissonnat and Oudot [6].Unlike in the latter method, we do not resample a smoothapproximation of the surface, but simply reinsert some sam-ple points that have been eliminated in the first stage till asmooth density gradient is achieved on the whole mesh andan aspect ratio criterion is met for every triangle of the re-constructed surface. For user-controlled detail insertion orsuppression, we propose to restart the convection process ina prescribed region with a value of the parameterρ modu-lated according to a potential field function.

To achieve sample point insertions or deletions in thereconstructed surface, we have designed an efficient algo-rithm that ’reinflates’ the pseudo-surface in an altered re-gion and restarts the convection process only locally. Thisalgorithm requires the computation of the 3D DT ofP ′ dur-ing the first stage, and supports either one or several inser-tion or deletion operations at a time.

In the following paragraphs, we first describe our al-gorithm for local update of the pseudo-surface. Then, wepresent our methods for improving the quality of the trian-gles and for interactive refinement.

4.1. Local update of the pseudo-surface

Let us recall that a pseudo-surfaceS′ that interpolatesa subsetP ′ of the original input point setP is embeddedin the 3D DT ofP ′. At the end of the convection process,every oriented facet ofS′ meets the oriented Gabriel prop-erty w.r.t.P ′. When a pointp is inserted (resp. removed) inthe 3D DT, the latter is modified onlylocally in a connectedregion spanned by the set of cells in conflict withp (resp.incident top). We exploit this locality property to updatethe convection result without restarting the convection pro-cess from the convex hull of the point set. Figure9 gives asimple example on a 2D point set, where one new samplepoint is inserted.

Figure 9. Local update of a pseudo-surface in 2D. The top rowillustrates the initial reconstruction process. In the bottom trow,a new (black) sample point is inserted. The conflict region isbounded by the dark contour. The DT is updated and the finalpseudo-surface is shown bottom-right.

At the beginning of the initial convection process, thepseudo-surface lies on the convex hull of the points ofP ′.All the Delaunay cells ofP ′ areinternal except the infinitecells1 that areexternal. During the convection process, thepseudo-surface evolves and the cells it goes through becomeexternal. An external cellC2 is said to have beendiscoveredfrom a cellC1 if C2 becomes external when the facet inci-dent toC1 andC2, oriented towardsC2, is opened by theconvection process. A cellC2 can be discovered onlyonce.In the case a facet oriented from a cellC3 towards a cellC2

is pushed towards a cellC2 that has already been discov-ered, it means that the pseudo-surface locally has coupledfacets betweenC2 andC3 that collapse. The cellC2 is saidto have beenrediscoveredfrom the cellC3.

The graph that represents the discovering relation be-tween cells is a forest of rooted treesD = ∪Di where eachtreeDi is rooted on an infinite cell. A given forestD is

1These cells are artifically created so that the facets of the convex hullare incident to two cells.

not unique in the sense there may exists several equivalentconfiguration of the discovering relation depending on theorder in which Delaunay tetrahedra are opened. Given acell C2 discovered from a cellC1 and rediscovered froma cell C3, the discovering and the rediscovering relationscan be switched if and only ifC2 is not an ascendant ofC3

(Fig. 10).

C1

C2

C3

C1

C2C3

Discovering relation.

Rediscovering relation.

Figure 10. Equivalence between two cell configurations in 2D.The discovering and rediscovering relations can be switched pro-vided there is no cycle creation.

When points are inserted into or removed from the 3DDT, some cells disappear, that we will call cellsin conflictin both cases. These are replaced bynew cells resultingfrom the local retriangulation. This alters the integrity ofthe discovering relation. The external cells that do not dis-appear and that were previously discovered from a conflictcell become roots of the discovering relation, though theyare not infinite cells. These cells are calledorphancells.To restart the convection process, we reinflate the pseudo-surface so that it encloses the conflict region (Fig.11,(a)).The newly created cells then become internal. The pseudo-surface is composed of facets oriented towards the interiorof the surface, that are of three types:• The restart facets, that are oriented towards new internalcells and that previously encoded discovering or rediscov-ering relations.• Thereversable facets, that are oriented from new internalcells towards orphan cells.• The other facets.

The update process consists in coherently restoring thediscovering relation regarding the new cell configuration in-duced by the updated new point set. This process can bedescribed through the following three steps:1. We first launch the convection process from the restartfacets and restrict it to the conflict region. The pseudo-surface stops temporarily at the boundary of this conflict re-gion, on standby of an opportunity to re-establish a discov-ering or rediscovering relation towards the outside. Thesefacets are calledtemporary facets(Fig. 11(b)).2. This step deals with orphan cells. We have to restorethe connectivity between cells so that these orphan cells getdiscovered, while maintaining a forest. Care must be taken

orphan cellconflict regioninitial pseudo−surface

(a) (b)reversable facets temporary facet

Figure 11. Convection through a conflict region in 2D. The con-flict region is bounded by the dashed curve. In (a), the convectionstarts from two entry points, corresponding to discovering and re-discovering relations. Gray cells are orphan cells. In (b), the con-vection through the conflict region has stopped. One temporaryfacet appears.

to avoid the creation of cycles while re-establishing the dis-covering relation. If that is not possible for a given orphancell, this cell cannot be external with the new point set con-figuration and the pseudo-surface will be reversed locally.

Let C denote an orphan cell. Suppose there exists a cellC1 adjacent toC, outside the evolving pseudo-surface andsuch thatC could be opened fromC1. If C is not an as-cendant ofC1 for the discovering relation, we create a dis-covering relation betweenC1 andC. Now thatC has beendiscovered, the facets that separateC from the interior ofthe pseudo-surface could be discovered by the convectionprocess. These facets are pushed into the temporary set.

If one of the candidate cell to beC1 is a new cell (i.e. in-side the retriangulated region), then it is chosen rather thanthe others. If there is no cell candidate to discoverC, thecell C becomes internal and the pseudo-surface backtracksconsistently. If a cellC2 was discovered fromC, then it be-comes an orphan cell (if not infinite) that will be processedin turn. The connectivity of the discovering relation is thenrestored.3. The last step consists in launching the convection processfrom the remaining facets in the temporary set.

4.2. Triangle quality improvement

Let f be an oriented facet ofS′. We call rmin (resp.rmax) the minimum (resp. maximum) decimation radius ofthe three vertices off . We define the following two refine-ment criteria:

• (Car) f meets the criterion if the ratio between the ra-dius of the circumcircle off and the length of its short-est edge is less than a positive constantβ.

• (Cden) f meets the criterion ifrmin

rmax≥ γ, whereγ is a

positive constant.

CriterionCar aims at guaranteeing a minimal aspect ra-tio for every triangle inS′. Criterion Cden prevents big

facets from being incident to small facets and so, ensuresa smooth gradient of density over the whole mesh.

The idea of the refinement algorithm is to ’break’ everytriangle ofS′ that does not meet one of the above condi-tions. In Chew’s algorithm, an interesting idea for produc-ing well-shaped triangles is to create a new sample pointequidistant from the three vertices of a facet with bad as-pect ratio. In our framework we extend this idea consider-ing the reconstructed pseudo-surfaceS′ and the remainingset of pointsR. Let f denote a facet that does not passCar

or Cden, and letBf be its minimal enclosing sphere. Forevery such facet, our algorithm proceeds by inserting inS′

the pointsof Bf ∩R that is the nearest from the intersectionpoint betweenS′ and the line that supports the dual Voronoıedge off . The algorithm stops either when both conditionsCar andCden are met or whenBf ∩R is empty. Our currentimplementation does not take boundaries of thin parts intoaccount, but it could be easily extended. Figure12 illus-trates one level of refinement for the APHRODITEmodel.

Figure 12. Triangle quality improvement for the APHRODITE

model. Starting from an initial reconstruction withρ = 0.95 (left),the refinement was performed withβ = 0.8 andγ = 0.5 (right).

4.3. Interactive refinement

Our interactive refinement technique consists in rescal-ing the value of the parameterρ locally and to update thereconstructed surface according to the new local samplingconditions. By being more or less restrictive locally on thenormal variation, we allow to add details or remove somefeatures with an intuitive control. We have implemented asimple brush tool based on a potential field function param-eterized by a centerc, a radiusr and a maximum intensityρc at c. For every point samplep is computed a local er-ror toleranceρloc(p) that depends on the distance betweenc andp. We define this local error tolerance as follows:

ρloc(p) = ρ +ρc − ρ

r2n(rn − ‖p− c‖n)2

where n ≥ 2 is an integer constant. This function

smoothly varies in a monotonically fashion betweenρc

(reached ifp = c), andρ (reached if‖p−c‖ ≥ r). If ρc > ρ,more points will be inserted. Otherwise, some points willbe removed. TheC2 nature of this function ensures a con-tinuous density gradient between the altered region and theremaining of the reconstructed surface.

To refine the reconstruction, we reinflate the pseudo-surface in the region of influence of the tool, and restartthe convection process in this region regarding the new lo-cal simplification parameters. Figures1, 13and14illustrateour interactive refinement tool on various point sets.

Figure 13. Interactive refinement of the Isis model. The originalmodel exhibits hieroglyphs on the back that are not captured for atoo small value ofρ (left, center-top). Our interactive refinementtool allows to make them to appear (center-bottom, right).

Figure 14. Interactive customization of the IGEA model. The orig-inal point set has 134K points. The reconstruction on the left hasonly 17K points. The ridge on the left cheek was interactively re-moved to obtain the result on the right.

5. Results and discussion

We have implemented our dynamic reconstructionframework on a Linux platform using the ComputationalGeometry Algorithm Library, CGAL2. We require CGAL’s

2http://www.cgal.org

Figure 15. Reconstruction of a noisy point set. In (a) is shown areconstruction of the original RAM model with 622K points, with-out simplification. In (b) is shown a simplified reconstruction ob-tained with our method, withk = 18 andρ = 0.94 (43K points).Our normal-based error metric based on a local normal estimationacts as a noise filter.

filtered predicates for robust point location and computationof Delaunay tetrahedra.

We demonstrate the effectiveness of our framework onseveral point set models that were obtained from laser rangescanning (Figs.1, 5, 14, 7), including a particularly noisypoint set (Fig.15). If normal directions are not supplied, theuser has to give a value for the size of thek-neighborhood.A value for the geometric error toleranceρ is required.The reconstruction stage then works automatically. Thecorrection stage can be either skipped or run given user-defined parameters. The triangle-quality improvement pro-cedure requires a maximum tolerated aspect ratioβ, and/ora minimum density factorγ. Interactive correction neces-sitates defining the brush tool properties and areas of in-terest picked on the reconstructed surface. Table 1 reportsthe overall timings and final number of points for the ini-tial reconstruction stage and for the correction stage. Allthe results presented here were obtained on a Pentium IV3.2GHz, 1GB RAM workstation. Reconstruction timingstake into account the incremental generation of the 3D DTof the simplified point set. This generation takes less than 5seconds for all the point sets we tested.

Our approach for locating points that form a Delaunaytetrahedron in Gabriel half-spheres is not currently optimal,which results in relatively high computation times. Somefacets, in particular on the convex hull, may have a Gabrielhalf-sphere that contains a great part of the input point set.Since the intersection between these half-spheres are not al-ways empty, some point samples can be tested many timesbefore they become part of the surface or they are elimi-nated. This search could be improved in several ways. Wecould benefit from the 3D DT of the points inserted in thesurface to improve the locality of the point locations or fur-ther exploit the normals to guess the position of the nextcandidate to insertion.

Model Reconstruction Correctionname #points ρ #points time #points time

TRIPLE HECATE 90,180 0.98 28,718 281 34,310 120SCREWDRIVER 27,152 0.98 7,944 49 – –

APHRODITE 46,096 0.90 4,507 67 7,644 42ISIS 187,644 0.95 8,368 96 10,994 38IGEA 134,344 0.98 17,232 102 17,104 24RAM 622,716 0.94 43,498 1,472 – –

Table 1. Performance of our reconstruction framework for variouspoint sets. Computational timings are given in seconds for bothinitial reconstruction and correction steps.

6. Conclusion and future work

In this paper, we have presented a new framework forreconstructing a surface from an unorganized point set thattakes only relevant sample points into account. In a firststage, we construct a triangulated surface that interpolatesonly a relevant subset of the input data. We decimate theinput point set on-the-fly during the reconstruction process.The sampling density is controlled by local geometric andtopological constraints. If needed, we then make correc-tions to the reconstructed surface, which requires the 3D DTof the simplified point set. We improve the quality of thetriangles by a refinement algorithm and enable interactiveinsertion of details or further simplification. These correc-tions are achieved in a dynamic fashion, without restartingthe reconstruction process from scratch, which makes ourreconstruction framework very flexible.

Future work will first include the search for a more effi-cient, dedicated data-structure for point locations. We couldthen extend our algorithm to automatically produce a mul-tiresolution decomposition of an input point set. This de-composition could be used for progressive reconstructionwith our dynamic update procedure. For non-uniformly dis-tributed point sets, it could be interesting to incorporate a re-laxation procedure into our algorithm. The question of giv-ing guarantees on the sampling density of the output pointset also certainly deserves further investigation.

Acknowledgments The authors would like to thankDavid Bourguignon for expressing the need for a dy-namic approach of geometric convection. The TRIPLE

HECATE model is provided courtesy of Christian Lahanierfrom C2RMF (Musee du Louvre, departement des Antiq-uites Etrusques et Romaines, marbre, MA 2594, C2RMF:FZ31452). The other point set models used in this paper arecourtesy of: Cyberware, Polygon Technology, and FaustoBernardini from IBM.

References

[1] M. Alexa, J. Behr, D. Cohen-Or, S. Fleishman, D. Levin, andC. T. Silva. Point set surfaces. InProc. IEEE VisualizationConference, pages 21–28, 2001.

[2] N. Amenta and M. Bern. Surface Reconstruction by VoronoiFiltering. InProc. Discrete Computational Geometry, pages481–504, 1999.

[3] M. Andersson, J. Giesen, M. Pauly, and B. Speckman.Bounds on the k-Neighborhood for Locally Uniform Sam-pled Surfaces. InProc. Symposium on Point-Based Graph-ics, pages 167–171, 2004.

[4] D. Attali, J.-D. Boissonnat, and A. Lieutier. Complexityof the Delaunay triangulation of Points on Surfaces: theSmooth Case. InProc. ACM Symposium on ComputationalGeometry, pages 201–210, 2003.

[5] J.-D. Boissonnat and F. Cazals. Coarse-to-fine surface sim-plification with geometric guarantees. InProc. Eurograph-ics, pages 490–499, 2001.

[6] J.-D. Boissonnat and S. Oudot. Provably Good Surface Sam-pling and Approximation. InProc. of Symposium on Geom-etry Processing, 2003.

[7] J.-D. Boissonnat and S. Oudot. An effective Condition forSampling Surfaces with Guarantees. InProc. Symposium onSolid Modeling and Applications, pages 101–112, 2004.

[8] F. Cazals and J. Giesen. Delaunay Triangulation based Sur-face Reconstruction: Ideas and Algorithms. Technical Re-port 5393, INRIA, November 2004.

[9] R. Chaine. A geometric convection approach of 3-D re-construction. InProc. Symposium on Geometry Processing,pages 218–229, 2003.

[10] L. P. Chew. Guaranteed-Quality Mesh Generation forCurved Surfaces. InProc. of Symposium on ComputationalGeometry, pages 274–280, 1993.

[11] O. Devillers. The Delaunay hierarchy.Internat. J. Found.Comput. Sci., (13):163–180, 2002.

[12] T. K. Dey, J. Giesen, and J. Hudson. Decimating samplesfor mesh simplification. InProc. Canadian Conference onComputational Geometry, pages 85–88, 2001.

[13] T. K. Dey, J. Giesen, and J. Hudson. Sample shuffling forquality hierarchic surface meshing. InProc. 10th Inter-national Meshing Roundatble Conference, pages 143–154,2001.

[14] H. Edelsbrunner. Surface reconstruction by wrapping finitepoint sets in space. In B. Aronov, S. Basu, J. Pach, andS.-V. M. Sharir, editors,Ricky Pollack and Eli GoodmanFestscrift, pages 379–404. Springer-Verlag, 2002.

[15] L. Linsen. Point cloud representation. Technical Report2001-3, Universitat Karlsruhe, Germany, 2001.

[16] C. Moenning and N. A. Dogson. A new point cloud sim-plification algorithm. InProc. Int. Conf. on Visualization,Imaging and Image Processing, pages 1027–1033, 2003.

[17] C. Moenning and N. A. Dogson. Intrinsic point cloud sim-plification. InProc. GraphiCon, 2004.

[18] M. Pauly, M. Gross, and L. P. Kobbelt. Efficient Simplifica-tion of Point-Sampled Surfaces. InProc. IEEE VisualizationConference, pages 163–170, 2002.

[19] J. Wu and L. P. Kobbelt. Optimized Sub-Sampling of PointSets for Surface Splatting. InProc. Eurographics, pages643–652, 2001.

[20] H.-K. Zhao, S. Osher, and R. Fedkiw. Fast Surface Recon-struction using the Level Set Method. InProc. IEEE Work-shop on Variational and Level Set Methods in Computer Vi-sion (VLSM), pages 194–202, 2001.