Pierre Alliez

Ron Kimmel

Guillermo Sapiro

Jacques Blanc-Talon

Tamy Boubekeur

Frédéric Chazal

Claire Lartigue

Yves Meyer

Freddy Bruckstein

Steve Oudot

José-Luis Lisani

Tony Buades

Adina Aude

Ayman

Benjamin

Eric Bruno

Julien, Marc, Mauricio Miguel Morgan

Zhongwei Merci

Anne Pascal

Keywords: Surfaces

Trois ans de travail sont résumés dans les pages qui vont suivre. Cela semble presque trop peu (les lecteurs me contrediront sans doute sur ce point).

Je voudrais remercier Jean-Michel Morel de m'avoir dirigée pendant ces trois ans. Cette thèse a été riche en apprentissages grâce à sa grande disponibilité et sa patience.

Notations

• M denotes a smooth surface • M S denotes a set of samples on M • p is a point of the smooth surface

• n is the normal to the surface at p • B r and C r are two different types of neighborhoods of p (that will be defined in section ..)

• r is the radius of the neighborhood

• m is a point of the neighborhood B r or C r

• o is the barycenter of the B r or C r

• k 1 and k 2 are the principal curvatures of the surface at p

• t 1 , t 2 are the principal directions of the surface at p

• H is the mean curvature of the surface at p

• V is a mesh and v, v i vertices of V , T is the set of triangles of V .

• g is a polynomial tel-00557589, version 1 -19 Jan 2011 tel-00557589, version 1 -19 Jan 2011

Chapter 1

Introduction Context

This thesis started as a part of a collaboration between two laboratories of the Ecole Normale Supérieure de Cachan, the CMLA (Centre de Mathématiques et Leurs Applications) and the LURPA (Laboratoire Universitaire de Recherche en Productique Automatisée). This project was named Géométrie Inverse Pour l'Industrie (Inverse Geometry for the Industry) and aimed at building a closed loop acquisition system. This thesis deals with the data analysis side of this project.

The LURPA high precision scanner is an experimental device delivering very large high precision raw data point sets with up to 35 million points, usually made of about 300 different scan sweeps.

In the course of the project it was progressively realized that the size of the raw point sets, and their precision (nominally about 20µm) raised questions which did not seem to have been addressed before. A first serious problem was that most state of the art techniques were actually too complex to be even tested on the data in reasonable time. Furthermore, the current meshing and analysis techniques were not designed for the conservation and visualization of high precision clouds, and usually lost accuracy and re-sampled or under-sampled the data. Smoothing and re-sampling were acceptable with noisy and only approximately registered clouds obtained by older generation scanners. Yet it seems that new opportunities for the visual exploration and conservation of valuable objects are opening up with future more precise scanners. As usual with better tools the first question is linked to the data proliferation.

This situation led us to a complete rewriting of the processing chain starting with a raw data set consisting of many scans with very irregular sampling and holes. While the end result remains a meshed surface, the steps leading to it have obeyed new requirements aimed at preserving the high accuracy throughout the processing chain:

• to be able to safely orient the raw cloud before meshing (raw orientation);

• to be able to mesh the raw point could itself, thus obtaining a raw mesh (raw meshing);

tel-00557589, version 1 -19 Jan 2011

• to be therefore able to visualize all scanning artifacts in the raw mesh (raw visualization);

• to correct the merging artifacts with minimal resolution loss (merging of raw scans) ;

• to delineate all missing parts (holes) to guide the laser head toward them, or, at the very least, to detect and fill them up automatically;

• to propose a raw cloud scale space analysis and segmentation method, based on the very same process leading to the surface construction.

A last problem arising in the course of this work and worth mentioning has been the scarcity of a shared data base of raw data point sets on which experiment sharing and performance benchmarks could be performed. Most available point clouds are meshed, which usually means that they are already smoothed and re-sampled. Thus most experiments were made on industrial objects and molded copies of archeological objects with fine texture and detail scanned at LURPA. Fortunately, we were also able to test our algorithms on a few raw data point sets provided by Stanford University (Forma Urbis Romae Project).

The main contributions of the thesis are reviewed in the next sections. For the sake of clarity, the contributions of the thesis are listed dryly, without the necessary bibliographic analysis. Of course, the bibliographic analysis and experimental comparison will be extensive in the other chapters.

Thesis summary

The problem of surface reconstruction from raw points is considered as essentially solved. Indeed, many commercial devices acquire objects and reconstruct their surface. Nevertheless, the reconstructed surface is usually very smooth as a result of the scanner internal de-noising process and the fusion of different scans. The precision loss is obvious in that all pieces look polished and glossy, having visibly lost their grain and texture. By looking at the academic data sets available on the web, a similar conclusion can be drawn. Surface reconstruction is, in fact, far from being adequately addressed in the context of high precision data.

This thesis has adopted the somewhat extreme conservative position, not to loose nor alter any raw sample throughout the whole processing pipeline. There are two main reasons to try to see directly the raw input samples. First, this could and will serve to control the acquisition loop. The input raw points have to be visualized, to point out all surface imperfections (holes, offsets) and correct them immediately. Second, since high precision data can capture the slightest surface variation, any smoothing, any sub-sampling can loose textural detail. Although tel-00557589, version 1 -19 Jan 2011  of course further compressing steps may become necessary, building and seeing first the highest quality reference data set before any compression seems to be a sound precaution. It also raises the hope of getting a real D microscope, enlarging considerably D objects and revealing their details.

The thesis attempts to prove that one can triangulate the raw point cloud with almost no sample loss, and solves the exact visualization problem for large data sets. Although this triangulation process performs remarkably well on accurate data sets, we shall see that it actually can be applied safely on any data set, even on extremely noisy point clouds.

As mentioned, the typical point sets treated here contain around 15 million points (up to 35 million points) and are made of some 300 different scan sweeps. To the best of our knowledge, the bibliographical analysis and numerical comparison will show that no reliable method existed yet for the exact visualization of all the input raw points: most methods either build a mesh by discarding a large part of the input samples, or simply create new vertices and loose the input samples. Large data sets have already been acquired and considered in the literature (Stanford Digital Michelangelo and Fragment Urbis Romae Project). Yet the available models resulting from these acquisitions are smoothed and re-sampled. They are generally not the initial scans.

To achieve the high precision triangulation which is the main goal of this thesis, two major problems will be solved. The first is the orientation of the complete raw point set, previous to meshing. The second one is the correction of tiny scan misalignments leading, nonetheless, to strong high frequency aliasing and hampering completely a direct visualization. Another outcome is the accurate delineation of scanning holes. The final result aimed at is an accurate visualization of the surface containing all raw points, with low frequencies slightly corrected to avoid aliasing effects, and intact high frequencies.

The second development of the thesis is a general low frequency-high frequency decomposition algorithm for any point cloud. This permits D extensions of the level set tree and of the MSER representations, which are classic image analysis tools, and an intrinsic mesh segmentation method.

The underlying mathematical development focuses on an analysis of discrete differential operators acting on raw point clouds which have been proposed in the literature. By considering the asymptotic behavior of these operators on a smooth surface, a classification by their underlying curvature operators is obtained.

Of particular interest is the discrete operator which ends up being the numerical spinal chord of the chain. It is, in one sentence, the iterated projection of each point on the local regression plane of its Euclidean neighborhood. This iterated operator will be proved to be consistent with the mean curvature motion (the intrinsic heat equation). It therefore defines a remarkably simple and robust numerical scale space analysis. By this intrinsic heat equation (using its numerical reversibility), all of the tel-00557589, version 1 -19 Jan 2011 above mentioned different and difficult problems (point set orientation, raw point set triangulation and scan merging), usually solved by separated techniques, are solved in a unified framework. Inasmuch as they can handle large amounts of points, the existing methods have been compared to the new algorithms (see for example fig .).

The next sections review one by one the chapters of the thesis and their contributions.

The thesis by chapters

State of the art (chapter )

This chapter gives an overview of the whole D surface processing field from the early acquisition to the final model rendering and post processing that can be done on the mesh models. It also reviews the multiple ways of acquiring D surfaces, acquisition methods being mainly divided into passive light acquisition methods and active light acquisition methods.

High precision point set orientation and meshing (chapter )

As stated above, the first goal is to visualize the initial raw point cloud including its very small details and its imperfections. All existing methods either smooth the shape or can simply not handle large data sets. Two major problems for D data set processing are solved using a new scale space implementation: the raw point set orientation, and the mesh construction. The result is a mesh whose vertices are almost all (in all examples more than .%) of the input samples. It follows that the input data can be properly visualized. The main tool that will be used to solve these problems is the operator projecting each point into the local regression plane of its Euclidean neighborhood. An asymptotic interpretation of this operator will be proved: This operator is consistent with the mean curvature motion tel-00557589, version 1 -19 Jan 2011

 ∂p ∂t = Hn
where p is a point of the surface, and H and n are the mean curvature and normal to the surface at p. The iterated projection provides a perfect reversible numerical scale space, able to push the points toward a smooth surface or, once reversed, to transport back on the initial cloud a structure which was easily computed on the smoothed one. Orienting a raw point set is not considered an easy problem: although normal directions are easy to compute, the ambiguity lies in finding a coherent orientation for the whole surface. Each point must be oriented consistently with its neighbors normals. This chapter proves the scale space method to work directly and reliably with the simplest implementation. The whole process is completely summarized in the following sentences. The point cloud is first smoothed by the scale space iteration, then the orientation of a seed point is chosen and its neighbors are oriented consistently with this normal. The orientation is then propagated from neighborhood to neighborhood. The final normals are transported back to the original positions of the points. Since sharp features are removed by scale space iterations, the propagation of normals is numerically straightforward.

Once a consistent orientation is found for the raw point cloud, the next step is to build a mesh from the scans. A mesh is the right shape model to generate efficient visualization. Chapter  shows that, again, the scale space method is selfsufficient to the task. The smoothed point cloud is easily meshed using any direct meshing algorithm such as the Ball Pivoting Algorithm ([BMR * ]). The resulting connectivity information (faces and edges of the mesh) can be transported back by reversing the scale space to the original point positions. As shown on figure ., this method permits to recover the tiniest details of an object acquired by a high precision laser scanner, which are lost by other methods. Since this operator recovers the exact initial point positions, the final mesh is not smoothed and scan superposition artifacts can also be visualized.

These results will be published in Scale Space Meshing of Raw Data Point Sets [DMMSL].

High fidelity scan merging (chapter )

Using the scale space meshing method, even the slightest scan misalignment is immediately visible (see fig .(b)). Smoothing the shape, for example by a bilateral filter, would suppress such artifacts, yet it would also remove detail and texture and make the whole meshing method useless, since it was devised precisely to preserve and see exactly the initial data. A new way for merging different scans without smoothing them was therefore necessary. Again, the scale space method proves sufficient to the task. Indeed, in a spirit similar to the low-high frequency decompositions in signal and image processing, the scale space operator can be used to decompose the shape into a low frequency base, the shape filtered by mean curvature motion (see fig. .(c)), and a high frequency component (the scale space motion vectors). This decomposition can also be seen as a smooth base + height decomposition: the smooth base is the shape obtained after smoothing and the height function associates to each point its displacement vector.

Thus, it is shown that all scans can be decomposed separately but consistently into their smooth bases and their intrinsic high frequencies. But they can also receive a common base by applying the scale space to their joint cloud. Adding the scans high frequencies to the common smooth base yields a final merging where all details of all scans are preserved and merged together. This process preserves integrally the high frequencies, but slightly moves the smooth bases (see Fig .(d)). The method proposed in this chapter is not a registration method: it only merges scans which are already well registered. But this process is essential to avoid very strong aliasing effect arising with even the slightest offset between scans.

The results of this chapter are published in High Fidelity Scan Merging [DMAL].

Filling holes in scale space meshes (chapter )

This chapter makes a review of mesh hole filling methods and deduces an adaptation of the Liepa [Lie] method to scale space meshes. Again, the method is based on the same scale space framework. It first detects and fills in the holes in the coarse scale mesh (the mesh built after the scale space iterations). In the regular scale space process, the backward operator should then be applied to get the fine scale mesh. Yet the backward scale space operator is not defined for the patch points. The patches can then be very conspicuous because of their smooth aspect.

tel-00557589, version 1 -19 Jan 2011



To avoid this, a texture noise is added to the patch points in their normal directions. The texture properties are deduced from the motion amplitude of the first scale space iteration. Indeed, as stated above, the first iteration captures the local variations and this is what these patches lack. Implicitly the hypothesis is that the patches contain no feature but are textured as the rest of the shape.

The level set tree on meshes (chapter )

A third contribution deals with the problem of mesh segmentation and feature extraction. Most feature detection methods focus on crest lines (also called ridge and valley lines), yet this kind of feature detection requires the computation of an order 3 surface derivative and leads to a local detection similar to the edge detection paradigm in image analysis. To get a real segmentation method, a more global mesh feature is introduced : the level lines of a proper function defined on the mesh.

Here the analogy with image analysis is useful. Given a real function defined on its triangles, the mesh can be seen as the support of an image, the triangles as its pixels, and the function as the image values. Then level lines and level sets of this function can be considered as mesh features and be used to segment it. To extract the most significant image level lines, the Maximally Stable Extremal Region algorithm was extended to meshes, yielding a time-efficient reliable level set selection. But the obvious question is: which real function can be chosen to intrinsically characterize the mesh? It turns out that the mesh curvature is the simplest lower order local intrinsic function that can be deduced from the mesh itself. By (again!) the scale space method, this curvature can be computed at any scale of smoothing and transported back on the raw pixels. Thus the image defined on the raw mesh simply is its own curvature at some fixed scale. As this chapter shows, segmenting this image defined on the mesh yields a mesh segmentation into smooth parts and high curvature parts. The scale space operator therefore permits to define an intrinsic mesh segmentation.

Experimental results show the segmentation and line extraction working as well on mechanical parts (fig .) as on fine arts objects (fig .).

These results are published in The Level Set Tree on Meshes [DMAMS].

A numerical analysis of raw point cloud smoothing (chapter )

For the sake of continuity in the exposition, this chapter comes late and is placed after the thorough description of the D scale space method and its application to the synthesis and analysis of raw point clouds.

But the chapter contains, nonetheless, the main mathematical results and the mathematical method to analyze the local discrete operators proposed in the literature and in the present thesis. Indeed, the definition and mathematical analysis of a raw numerical scale space contributed in Chapter  leads to a general methodology for analyzing discrete methods. These methods attempt to compute differential operators on irregularly sampled surfaces by discrete schemes based on the local point cloud statistics. Many such methods have been proposed for meshed D data. However, computing directly differential operators on the raw point clouds as they are acquired (e.g.) by triangulation scanners is crucial because this can be done before any mesh resampling and can in particular bring decisive information into a meshing loop.

tel-00557589, version 1 -19 Jan 2011

This chapter proposes a method to analyze and characterize these raw point cloud local operators. It reviews a half dozen basic discrete algorithms which have been proposed to compute discrete curvature-like shape indicators on raw point clouds. It shows that all of them can actually be analyzed mathematically in a unified framework by computing their asymptotic form when the size of the neighborhood tends to zero. Assuming that the underlying manifold is smooth enough, the asymptotic form of these operators is obtained in terms of the principal curvatures or of higher order intrinsic differential operators which characterize the discrete operators. This analysis, completed with numerical experiments, permits to classify the discrete raw point cloud operators, to rule out several of them, and to single out others.

tel-00557589, version 1 -19 Jan 2011  Furthermore, by analyzing asymptotically two very simple moving least squares operators, namely the operator M LS1(p) that projects each point p to the local regression plane (our favorite scale space operator) and the very common operator M LS2(p) that projects each point onto the local degree 2 polynomial regression surface, it is shown that only the first is consistent with the mean curvature motion. More precisely, one has:

< M LS1(p) -p, n >= r 2 8 (k 2 + k 1 ) + O(r 3 ),
where r is the radius of the neighborhood used for the surface regression, k 1 and k 2 are the principal curvatures and n is the normal to the surface at p, oriented towards the concavity.

On the contrary, one has < M LS2(p)p, n >= -r 4 48 (3a 04 + a 22 + 3a 40 ) + O(r 5 )

where

a 40 = 1 4! ∂ 4 f ∂x 4 , a 04 = 1 4! ∂ 4 f ∂y 4 , a 22 = 1 4! ∂ 4 f
∂x 2 ∂ 2 y are the fourth derivatives of the graph function of the manifold in the intrinsic coordinates system at p and x, y are the coordinates along the principal directions.

These results will be published in [DM] Color cloud visualization (chapter ) Publications linked to the thesis

Outline of the thesis

The present thesis is divided as follows:

• Chapter  reviews state of the art work related to data acquisition, representation and extraction of characteristics.

• Chapter  describes the first applications of the new scale space: the robust orientation of surface point clouds and the construction of a high fidelity mesh whose vertices are the raw points.

tel-00557589, version 1 -19 Jan 2011



• Chapter  introduces a method to merge scans to prevent the appearance of D aliasing due to scan misalignment.

• Chapter  describes a postprocessing step to detect and fill the holes in the built meshes.

• Chapter  describes the extension of the image level set tree theory to surface meshes.

• Chapter  reviews and analyzes mathematically the standard ways of computing the curvatures and curvature directions for raw point clouds.

• Chapter  describes another application of the projection filter: the filtering of color clouds extracted from color images.

. D Acquisition

There are multiple ways of acquiring D surfaces. We shortly explain how these data can be acquired. Acquisition methods are mainly divided into passive light acquisition methods and active light acquisition methods. Passive light acquisition methods do not project any light onto the object. This is the case of stereo acquisition for example. On the other side, active light scanner use light projection to acquire the geometry. There exist contact D scanners that acquire a surface by physical touch, yet here we will only consider non-contact D scanners.

.. Stereoscopy

Stereo acquisition is the principal passive light acquisition method: it creates a disparity map by considering two images in epipolar geometry and using the y axis difference between matching points to deduce a depth map. This is very useful in satellite imaging for creating digital elevation models for example (see tel-00557589, version 1 -19 Jan 2011

.. Triangulation laser scanner

Triangulation laser scanner are named this way because of the triangle formed by the laser emitter, the camera optic center and the laser impact point on the surface. Despite the name, triangulation laser scanners do not automatically produce a triangulation. They initially produce data points in a D coordinate system without any other information than the coordinates. Normal information could be deduced easily from the laser emitter position, yet being able to reliably process un-oriented points widens the application field of the algorithms. The most famous acquisition projects, Digital Michelangelo and Forma Urbis Romae (see section .), have been done using triangulation laser scanner. For example, the laser system used for the David range acquisition " employed a  mW -nanometer laser diode, a  x  pixel CCD sensor, and a fixed triangulation angle" ([LPC * ]).

Figure .: Scheme of a triangulation laser scanner acquisition system (copyright Wikipedia)

Efforts have been made to combine information of various kind during the acquisition process to improve the reconstruction quality: [NRDR] combines normals computed by photometry and range images and [MYF] integrates shape from shading and range images (see also [Jas], [DSGV]).

tel-00557589, version 1 -19 Jan 2011

. Some acquisition projects

Over the past years, some acquisition projects have been completed yielding new data sets for the geometry processing community. In particular, the Stanford Digital Michelangelo Project  aimed at building D models from some of Michelangelo's most important sculptures, including the David ([LPC * ], [BRM * ]).

The Stanford Forma Urbis Romae Project  dealt with a marble map of Rome designed in the early rd. century. This marble map was broken into several pieces, some of which were lost or came to us only through drawings. Trying to solve the jigsaw puzzle and putting back all the pieces together is a challenging task. The remaining pieces were therefore acquired by triangulation laser, which yielded another data set that is partly available for research ([KTN * ]). This data set will be used for experimentation in this thesis.

Most data used in this thesis come from the acquisition done by another laboratory: the Laboratoire Universitaire de Recherche en Productique Automatisée (LURPA). The next section shortly describes the acquisition system.

. LURPA acquisition system

The LURPA acquisition system is a triangulation laser scanner with high precision (around 20µm). It is composed of a granite table, a revolving arm and a laser scanner head (fig .). The system projects a laser pencil on the surface whose position is captured by a CCD camera. The device is calibrated so that it can translate position coordinates given in the CCD coordinate system into the machine coordinate system. This process of registering scans together will be detailed in ... Scans are registered to a very high precision but there is always a remnant offset which, as we shall see on the raw mesh, creates strong aliasing effects (see chapter ). The output of this device is a set of unorganized un-oriented points given only by their D coordinates as can be seen on Fig. ..

. Representation of D data

The acquired shape comes usually as a set of points sampled with more or less precision on the object surface. A natural question is then to find a way to represent the surface in a more handy way than a list of D coordinates. The usual choice is to build a mesh: i.e., a set of connected triangles. This makes a piecewise linear approximation of the surface. The more vertices, the finer the approximation. Other surface representation include splines and NURBS surfaces  http://www-graphics.stanford.edu/data/dmich-public/  http://formaurbis.stanford.edu/docs/FURproject.html tel-00557589, version 1 -19 Jan 2011 ( [PT], [MK], [Pet]). An alternative to meshes was presented in [LP], where the shape is represented by a set of triangle fans. A triangle fan is a set of triangles sharing a vertex. This representation allows for faster computations because of the localness of the fans.

In this thesis, we show that meshes can be built on the original data points, thus loosing no information whatsoever on the raw data point set. In that way all scanning artifacts are revealed and can be carefully corrected, in a way that does not loose the initial point accuracy.

There are previous methods attempting to keep the point cloud information all along the process and to perform the final rendering on the point set itself, without meshing [LMR], yet this is not our choice here. The whole theory of Moving Least Squares Surfaces [Lev] also aimed at processing surfaces by local surface regression, without any previous meshing (see chapter ). [MMS * ] proposed a meshless subdivision framework based upon the idea that point-based surface processing relies on intrinsic surface properties instead of intermediary representations. This idea of raw point processing is also at the bottom of our methods. Defining geodesics on point clouds has been investigated in [MS] and [MSa] where theoretical results are given for computing the distances on point clouds. [MSb] and [MS] use the Gromov-Hausdorff distances to compare manifolds given by point clouds. Another formulation was proposed in [BBK]. [MS] use geometric distance distributions on point clouds to recognize shapes.

Other works on point clouds include [BSW] where an algorithm for building a Laplacian operator on point clouds is given and [MTSM] which computes the visibility of a point cloud from a viewpoint.

The next section discusses how to go from an input raw point set to a mesh.

. From a point cloud to a mesh

We assume that the scans are already registered, i.e., that all the scans are given in a global coordinate system, obtained by the device calibration or by previous registration algorithm (in the case of most data processed here the global coordinate system is given by the scanning device itself). Other acquisition systems do not necessarily provide this initial registration. This is why much work has been done on registering shapes to merge multiple view (fig .).

.. Merging multiple views

The registration of multiple views is usually done by finding correspondences between the views and by finding a rigid or non-rigid transform mapping the set of points from the first view to the corresponding set of points in the second view.

tel-00557589, version 1 -19 Jan 2011 Finding these initial correspondences requires a good shape descriptor. This topic has generated a huge literature. Famous shape descriptors include spin images ( [JH], [JH], [Joh]), snapshot descriptors ([Mal]), regional point descriptors [FHK * ], spherical representations ( [HID], [KFR]), harmonic maps ( [ZH]), point signatures ( [CJ], [YF]) and D point's fingerprints ( [SA]). The heat kernel was also used to produce a point signature ( [SOG]).

Registration has also been studied in the shape matching context (e.g. to compare a shape to models in a database) [BBK * ], [BBK], [BBBK]. This is particularly difficult when dealing with deforming surfaces, for example, to recognize faces in spite of expression change. Specific surface descriptors had to be introduced ([BBK], [BBKa], [BBKb], [BBK]).

Once the initial registration has been found, the optimization is done by iteratively finding the registration that minimizes the distance between the overlapping parts of the scans and updating the point positions ([BM], [RL], [Hor]).

Since scans may contain warps, a rigid transform might not be enough to generate a good model, this is why non rigid transforms were considered ([ZZW * ], [CR], [BR], [BR]). More specifically the transformation is not a rotation and translation anymore, but a thin plate spline allowing for distortions inside each scan.

Once scans are registered in the same coordinate system, points must be turned into a single mesh. The next subsection reviews the details of meshing.

.. Meshing

As soon as a set of points is given and we look for a suitable model, the interpolation/approximation dilemma comes in. Should the model include the initial data (i.e. the mesh vertices remain the input data set) or should it approximate the shape (the mesh is close to the initial points but does not necessarily contain the data)?

First, a mesh can be built by creating triangles between input samples. This is tel-00557589, version 1 -19 Jan 2011 the case of all Voronoi/Delaunay based reconstructions. Indeed by building the Delaunay tetrahedralization of a point set and labeling inside and outside tetrahedra, a set of frontier triangles can be extracted which are "on" the surface ( [ACK],

[Boi], [ABK], [AB]). Yet these methods do not handle properly large input point sets. Edelsbrunner and Mecke introduced the concept of α-shapes ( [EM], [Ede], [BB]). The α-shapes are based on the Delaunay diagram of the input samples and are used to build a sub-complex of the Delaunay simplicial complex.

A triangle of the Delaunay triangulation is kept iff there exists a ball with empty interior or the complement of a ball of radius α that contains the whole pointset and such that the three vertices are on its frontier. The Ball-pivoting algorithm [BMR * ]) is based on this principle but avoids building the Delaunay triangulation. It simply rotates a ball of radius r on a set of points and builds triangles whenever the ball has three vertices on its frontier and none inside. For α-shapes as well as for ball-pivoting reconstruction, the problem is setting the ball radius which is a good compromise between detail loss and hole filling.

Over the past few years, another family of shape reconstruction methods has taken over the field. It is based on finding the level set of a function that more or less corresponds to the distance to the underlying shape (see Fig  .). Some methods assume that the points are consistently oriented, and use various function family to approximate the signed distance field: Fourier basis or radial basis functions for example ([HDD * ], [Kaz], [KBH], [ACSTD], [CBC * ], [OBA]), while more recent methods deal with unoriented point sets ([MDGD * ]). An interesting method uses the Voronoi diagram in the level set framework: the distance function is approximated by radial basis functions centered at the voronoi cell centers ( [SAAY]). Once the distance function is built, the level set is extracted using the Marching Cubes algorithm ( [LC]) or its extension ( [KBSS]). Other possibilities include [WH], [CA] and [LGS].

tel-00557589, version 1 -19 Jan 2011

.. Remeshing

Meshes built from raw scans might have a high number of vertices. Therefore, reducing the number of vertices (compressing the mesh data) is of crucial importance for practical uses. Though remeshing will not be addressed in this thesis (it will only be used in the hole filling process), we summarize briefly the remeshing field (a review of remeshing techniques can be found in [AUGA]). Decreasing the number of vertices and triangles while keeping the shape as close as possible to the original was investigated in [HDD * ], [PGK] and [CSAD] among others. The aim of remeshing can also be to improve the mesh regularity ( [AMD], [BK], [AVDI]). Usually the remeshing is done to satisfy a remeshing criterion and many such criteria exist (see [ABE]). Another popular remeshing topic is to transform a triangle mesh into a quadrangle mesh. Indeed quadrangles are attractive because of their tensor-product like nature, which is useful for mapping textures (e.g.) [ACSD * ], [AUGA], [SLI] or NURBS patch fitting in reverse engineering. Other methods have been proposed to redistribute vertices according to a density function ( [PC], [PC], [PC]).

In a related domain, [LCOLTE] proposed a projection algorithm to project an arbitrary point set onto an input point-set which can be used to resample surfaces by point sets.

. Rendering

Once the mesh is built, the object must be rendered to be visualized. This can be done by Ray Tracing methods. The idea is to trace a ray from the optic center of a virtual camera through each pixel in a virtual screen, and computing the color of the object that this ray intersects (see fig .). It requires computing mesh-ray intersections which can be very slow in case of large meshes. Models can be textured and various types of lights can be set to get realistic renderings.

. Post-processing and line detection

When surfaces are built, their geometry can be analyzed. For example, we might want to extract high curvature lines. This is always a very delicate part since high curvature lines, crest and ridge lines would actually need computing zero crossing of principal curvatures derivatives. This is an order 3 derivative of a surface we know only by a few sampled points. In case of a very smooth surface it could actually work, but it proves very unstable. Chapter  will introduce a new, global way for extracting geometry from a shape. In accordance with the edge detection paradigm in image processing, it is common to perform a D shape analysis by extracting the crest lines (the real edges) on meshes or point clouds. Ridge lines are the loci of points where the maximal curvature takes a positive maximum along its curvature line. Valley lines are the loci of points where the minimal principal curvature attains a negative minimum along its curvature line. These points can be linked to form lines (see among others [OBS], [BA], [LFM], [YBS], [SF]). Most methods use a quadric or polynomial regression. In [GWM], the lines are detected by neighborhood covariance analysis. Indeed, from a point neighborhood, the centroid and centered covariance can be computed. Comparing the ratios of the covariance matrix eigenvalues gives the geometry of the neighborhood (see also [MOG]). In [HMG], edges of a mesh are first classified according to their importance (this importance is an increasing function of the adjacent faces angle).

tel-00557589, version 1 -19 Jan 2011

A multiscale approach was proposed in [PKG]. Nearby feature points are first detected. In the neighborhood of these points surfaces are fitted, and depending on the number of fitted surfaces, points are projected to the nearest surface. Intersection points of these surfaces are finally classified as edge or corner points. By increasing the processing radius, one could track feature lines and keep only the ones at a given scale. Though dealing with scales, this method does not introduce a scale space framework. A similar idea for points classification and point projection was used in [DIOHS].

Although these papers introduce a ridge/valley line detection, none of them proposes a ridge and valley segmentation. In [IFP] the idea was suggested, tel-00557589, version 1 -19 Jan 2011 though: points lying near ridges or near valleys were labeled and this labeling was used to obtain a better rendering of the ridge and valley lines. But crest lines as defined by these methods require the computation of degree three surface derivatives.

Chapter  focuses on other interesting and well defined line features: namely the curvature level lines and level sets, analogous to the image gray level lines. Of particular interest are the zero-crossings of the curvature, which are technically similar to the zero-crossings of the Laplacian in image processing. These zero-crossings define inflexion lines, easy to compute from the raw data point set. They reliably segment the surface into ridges and valleys.

. The difference between Image Processing and D surface processing

Both subjects are very related and indeed share common ideas. Nevertheless, there is a major difference: there is no equivalent of the Shannon sampling theory for D surfaces embedded in three dimensions, no Fourier analysis, and no notion of frequency domain. This is why so many algorithms cannot be easily adapted to surfaces. For example, the very powerful non local means algorithm [BCM] was considered for adaptation to meshes by Yoshizawa et al. ([YBS]) yet the sampling problem was not really handled. It was once more adapted in [WZZY], more properly so that neighborhoods were better defined, yet it made the whole algorithm much more complicated. As will be seen in Chapter , however, a (non linear) high frequencylow frequency decomposition of the surface is actually possible thanks to the scalespace method developed in the present thesis. It provides a way to create for each surface a smooth base on which the high frequency component can be defined as a refinement field. This (simple) decomposition will be used for both the processing (chapter ) and the analysis of the surface (chapter ). tel-00557589, version 1 -19 Jan 2011

. Introduction

A growing number of applications involve creating numerical models for existing objects acquired by triangulation laser scanner or other devices. Commercial scanners can directly produce a direct triangulation of points sampled on the surface, but this triangulation is derived from a raw set of points with no connectivity information. Only raw input data will be considered here, namely sets of unorganized and non-oriented points given by their x, y, z coordinates. The proposed method orients and meshes directly the complete raw data set, thus allowing for the visualization of the finest surface details, and an accurate delineation of the scanning holes. The processed point data have a typical acquisition error of 20µ, allowing in principle to recover the finest texture and details.

The main tool introduced here is a raw data set point smoothing operator consistent with the intrinsic heat equation. The intrinsic heat equation, or mean curvature motion (MCM), is the simplest intrinsic method to smooth out a surface. The mean curvature motion writes

dp dt = Hn (.)
where H is the mean curvature at p (whose sign depends on the normal orientation), and n the normal. This motion will be given a robust implementation working directly on raw data points, which can be summarized in few words: it is the iterated projection of each point on the regression plane of its radial neighborhood.

tel-00557589, version 1 -19 Jan 2011

Mathematical and experimental arguments will show that this iterated planar regression consistently implements the MCM and actually permits to compute an accurate denoised curvature. Indeed, Theorem  states that, by these iterations, each raw data set point moves forward at the speed of the surface mean curvature in the direction of the surface normal.

By the iterated projection algorithm each initial raw data point can be tracked forward and backward in the surface smoothing process. As a consequence, the surface structure established at a smooth scale can be transported back on the raw data set point. This back transportation yields a topologically faithful orientation at each raw point, and subsequently a mesh whose vertices are almost all raw data points. It also permits an accurate detection of holes in the raw data, useful for further scanning attempts. Comparative experiments will show that a direct meshing gives poor results, while the back transported mesh allows for the uttermost accurate rendering of the surface, the mesh vertices being more than % of all initial raw points. Obviously, such a complete mesh is not economical, but it permits an accurate rendering of fine art or archaeological pieces at 20µ precision and a detection by visual inspection of the tiniest scanning defects.

The use of the mean curvature motion, forward and backward, is a direct D extension of the scale space paradigm in image processing introduced in the founding Witkin paper [Wit]. It consists of applying the heat equation ∂u ∂t = ∆u to the image u, which entails a rapid image simplification. The main image features (for example the edges) are detected at a coarse scale (large t) and then tracked back to their fine scale position. The next subsection reviews the methods computing curvatures and normals on raw data.

.. Building a mesh

Given an initial oriented point cloud, most meshing methods start by defining a signed distance field to the inferred surface [HDD * ], [KBH]. The signed distance function can be estimated at any point by computing the distance between the point and the regression plane of its k-nearest neighbors [HDD * ]. Since the neighbors are assumed previously oriented, the sign of this distance is straightforward. Other successful methods approximate the distance function using its decomposition on a local radial basis functions [KBH]. Once the distance function is defined, extracting the surface corresponds to extracting the zero level set of the distance function. This can be done using the marching cubes algorithm [LC] or any other other contouring algorithm.

These methods yield meshes that approach well the shape, but the approximation entails an implicit surface smoothing and the loss of fine texture. Acquisition holes are also filled in by those methods, the signed distance function giving a natural close up of the surface. Nonetheless, for scanning applications, the acquisition tel-00557589, version 1 -19 Jan 2011 holes should be detected rather than filled in. The smoothing can be desirable if there are noise and scanning artifacts. However, in the cases we shall examine, texture actually dominates noise. A guarantee that no detail is lost is granted only when almost all raw data points are mesh vertices. [AMD] introduces a remeshing method based on mappings of the mesh onto a plane. Meshing a planar projection will also be used here, but this projection will only be local.

In a way the scale space meshing method is not far from [GKS], where the triangulation is found by locally projecting the points onto a regression plane and performing a D triangulation. Our method will also consider meshing a simpler point set yet it uses a D triangulation method and is done in the scale space framework defined below.

.. Raw data point set processing

Yet, it has been termed impossible to mesh directly the raw data point set. The literature has therefore considered more and more sophisticated smoothing and interpolation methods. The "Moving Least Square Surface" (MLS) introduced in [Lev] is defined as the set of stationary points of an operator projecting each raw point on a local quadratic regression of the surface. The order n MLS algorithm estimates at each point a degree n polynomial surface from a set of weighted neighbors. The obtained surface can be used to project the point on the MLS surface, or to sub-sample the surface by removing step by step the points with least influence [ABCO * ]. Variations of the MLS algorithm for denoising point sampled surfaces and preserving edges were proposed in [FCOS], [GTE * ], [OGG], [LCOL], [GG]. Interpolatory point set surfaces can be achieved using a singular weighting function ( [OGG], [AA], [SOS]), but extracting the isosurface via marching cubes will loose the input point positions At first sight applying MCM to a data point set requires the separate computations of the surface intrinsic Laplacian (mean curvature) and of the normal. For meshes, the standard discretization of the Laplacian operator is done through the cotangent formula [MDSB]. For point clouds, [BSW] proposed the construction of a laplacian operator for functions defined on point clouds (yet no result on real noisy shapes was presented). In [PKG], the curvature is either estimated by a polynomial regression or by projection on a fitted least square surface (in other terms, by MLS). The reverse operator is built by storing the displacements of each point at each step. A similar scale space approach will be used here, but with quite different scopes: in [PKG], the proposed applications were shape morphing and shape editing.

In [UH], another raw data point set MCM discretization was proposed. The surface Laplacian is computed by building an operator A θ at each point position and for every direction θ in the tangent plane. A θ moves a point p proportionally tel-00557589, version 1 -19 Jan 2011 to the curvature H θ of the section curve in direction θ. By integrating over θ, it yields a mean curvature motion.

.. Computing curvatures

Computing the principal curvatures reliably on a given surface is crucial for various applications, in particular the anisotropic filtering preserving sharp edges [HP], [MDSB], or the sampling methods adapting the density to the surface curvatures [PGK]. On meshes, the curvature estimation problem has already been investigated in [MDSB] where the cotangent formula is proven and extended. [Taua] derived an analytic expression for estimating the directional curvatures in the edge directions. In [Rus], [TRZS], the tensor curvature was estimated on each face of a mesh surface. Other mesh curvature computation techniques include the use of the normal cycle theory [CSM]. For a summary and comparison of mesh curvature estimation methods, see [MSR]. It is also possible to estimate curvatures by building curves contained in the surface and passing through the considered point [Tan].

To determine the curvature of a given point, direct methods fit a surface (a polynomial or a quadric) locally to each neighborhood and then compute the fundamental forms in their explicit form. This permits to compute the Weingarten map whose eigenvalues and eigenvectors are the principal curvatures and principal directions ( [LFM] among others). In [BC] the principal curvatures are computed from an oriented raw data set without surface fitting by expressing the fundamental forms of a D surface as covariance matrices. Indeed, the covariance matrix of the point normals projected on the regression plane yields the principal curvatures and their directions. Other approaches avoiding surface regression include the computation of integral invariants [PWHY],[PWY * ]. They are based on the idea that differentiation is not robust in a discrete and potentially noisy data set, whereas integration is much more resistent to noise. The proofs link the computation of the area of the intersection of the surface with a ball to the principal curvatures. Another possibility is to adapt the curvature estimation of [Taua] to the case of point clouds, like in [LP]. Instead of considering the edge direction, since no edge information is given for the point cloud, they consider all directions from the center point to one of its neighbors. MLS surfaces were also used to derive analytic expressions for the curvatures of point set surfaces [YQ]. As far as meshes are concerned, a comparison of various curvature estimations can be found in [SMS * ].

Mathematical results are given in chapter  proving the consistency of the proposed scale space algorithm. This chapter is divided as follows: section . analyzes the discretization problem. Sections ., . describe the two main applications of the scale space: a point cloud orientation method and a faithful mesh construction

. Continuous theory

This section investigates a new way of implementing the mean curvature motion by the iteration of a planar surface regression. The surface M supporting the data point set is assumed to be at least C 2 . The samples on the surface M are denoted by M S .

Let p(x p , y p , z p ) be a point of the surface M. At each non umbilical point p, consider the principal curvatures k 1 and k 2 linked to the principal directions t 1 and t 2 , with k 1 > k 2 where t 1 and t 2 are orthogonal vectors. (At umbilical points, any orthogonal pair ( t 1 , t 2 ) can be taken.) Set n = t 1 × t 2 so that (t 1 , t 2 , n) is an orthonormal basis. The quadruplet (p, t 1 , t 2 , n) is called the local intrinsic coordinate system. In this system we can express the surface as a C 2 graph z = f (x, y). By Taylor expansion,

z = f (x, y) = 1 2 (k 1 x 2 + k 2 y 2 ) + o(x 2 + y 2 ). (.)
Notice that the sign of z depends on the arbitrary surface orientation.

.. Spherical neighborhoods vs cylindrical neighborhoods

Consider two kinds of neighborhoods in M for p defined in the local intrinsic coordinate system (p, t 1 , t 2 , n):

• a neighborhood B r = B r (p) ∩ M is the set of all points m of M with coor- dinates (x, y, z) satisfying (x -x p ) 2 + (y -y p ) 2 + (z -z p ) 2 < r 2 • a cylindrical neighborhood C r = C r (p) ∩ M is the set of all points m(x, y, z) on M such that (x -x p ) 2 + (y -y p ) 2 < r 2 .
For the forthcoming proofs the cylindrical neighborhood will prove much handier than the spherical one. Proof. The surface area element of a point m(x, y, z(x, y)) on the surface M, expressed as a function of x, y, dx and dy is dm(x, y) = 1 + z 2 x + z 2 y dxdy. One has

z x = k 1 x + O(r 2 ) and z y = k 2 y + O(r 2 ). Thus dm(x, y) = (1 + k 2 1 x 2 + k 2 2 y 2 + O(r 3 ))dxdy which yields dm(x, y) = (1 + O(r 2 ))dxdy. (.)
Using (.), the integrals we are interested in become

Cr f (x, y)dm = (1 + O(r 2 )) Br f (x, y)dxdy; (.) Br f (x, y)dm = (1 + O(r 2 )) Cr f (x, y)dxdy (.) = (1 + O(r 2 ))
x 2 +y 2 <r 2 f (x, y)dxdy.

This last form is amenable to analytic computations. Consider polar coordinates (ρ, θ) such that x = ρ cos θ and y = ρ sin θ with -r ≤ ρ ≤ r and 0 ≤ θ ≤ π.

Then for m(x, y, z) belonging to the surface M, we have z = 1 2 ρ 2 (k 1 cos 2 θ + k 2 sin 2 θ) + O(r 3 ). Fixing θ we obtain a curve with equation z = 1 2 ρ 2 k(θ) + O(r 3 ), where k(θ) = k 1 cos 2 θ + k 2 sin 2 θ. The condition that (x, y, z) belongs to the neighborhood B r (p) can therefore be rewritten as ρ 2 + z 2 < r 2 , that is

ρ 2 + 1 4 k(θ) 2 ρ 4 < r 2 + O(r 5 )
Computing the boundaries ±ρ(θ) of this neighborhood yields ρ(θ

) 2 + 1 4 k(θ) 2 ρ(θ) 4 -r 2 + O(r 5 ) = 0. Thus ρ(θ) 2 = -1 + 1 + k(θ) 2 (r 2 + O(r 5 )) 1 2 k(θ)
f (x, y)dxdy = [0,2π] [0,ρ(θ)] f (x, y)ρdρdθ = [0,2π] [0,r] f (x, y)ρdρdθ -E = Cr f (x, y)dxdy -E, with E =: [0,2π] [ρ(θ),r] f (x, y)ρdρdθ. Thus |E| ≤ π 4 sup x 2 +y 2 ≤r 2 |f (x, y)|k(θ) 2 r 4 , which yields |E| ≤ π|k 1 | 2 4 sup x 2 +y 2 ≤r 2 |f (x, y)|r 4 . In particular if f (x, y) = O(r n ), then |E| ≤ O(r 4+n ).

.. Curvature estimation

Theorem  states that projecting a point onto the neighborhood barycenter approximates the mean curvature motion. We shall discuss later on why it cannot be used for implementing the mean curvature motion.

Theorem . In the local intrinsic coordinate system, the barycenter of a neighborhood B r where p is the origin of the neighborhood has coordinates

x o = o(r 2 ), y o = o(r 2 ) and z o = Hr 2 4 + o(r 2 ), where H = k 1 +k 2 2
is the mean curvature at p.

Proof. By Lemma  applied to the numerator and denominator of the following fraction, we have

z o = Br zdm Br dm = x 2 +y 2 <r 2 z(x, y)dxdy + O(r 5 ) x 2 +y 2 <r 2 dxdy + O(r 3 ) = x 2 +y 2 <r 2 1 2 (k 1 x 2 + k 2 y 2 ) + o(x 2 + y 2 ) dxdy x 2 +y 2 <r 2 dxdy + O(r 3 ) = 1 2πr 2 r ρ=0 2π θ=0 ρ 2 (k 1 cos 2 θ + k 2 sin 2 θ)ρdρdθ + o(r 2 ) = r 2 8π (k 1 π + k 2 π) + o(r 2 ) = Hr 2 4 + o(r 2 ).
A similar but simpler computation yields the estimates of x o and y o .

tel-00557589, version 1 -19 Jan 2011

.. Surface motion induced by projections on the regression plane

The main tool of the proposed scale space will be simple projection of each surface point p on the local regression plane. This PCA regression plane is defined as the plane orthogonal to the least eigenvector of the centered local covariance matrix, and passing through the centroid of the neighborhood. The projection of p on this plane will be called p ′ . The next lemma compares the normal to the PCA regression plane with the normal to the surface, n.

Lemma . The normal v to the PCA regression plane at p ∈ M is equal to the surface normal at point p, up to a negligible factor: v = n + O(r).

Proof. The local PCA regression plane of point p is characterized as the plane passing through the barycenter of the neighborhood B r and with normal v minimizing:

I(v) = Br | v, pm | 2 dm s.t. v = 1
Denoting by (v x , v y , v z ) the coordinates of v,

I(v) = Br (v x x + v y y + v z 1 2 (k 1 x 2 + k 2 y 2 ) + o(r 2 )) 2 dxdy.
Considering the particular value v = (0, 0, 1) shows that the minimal value I min of I(v) satisfies I min ≤ O(r 6 ). In consequence the minimum

(v x , v y , v z ) satisfies v x ≤ O(r) and v y ≤ O(r). Thus v z ≥ 1 -O(r) and therefore v = n + O(r).
By Lemma , projecting p onto the regression plane induces a motion which is asymptotically in the normal direction: p ′ p is almost parallel to n. The simple projection of each surface point p onto its local regression plane approximates a D scale space (mean curvature motion) as shown in the next theorem.

Theorem . Let T r be the operator defined on the surface M transforming each point p into its projection p ′ = T r (p) on the local regression plane. Then

T r (p) -p = Hr 2 4 n + o(r 2 ). (.)
Remark More extensive formal calculations (too long to be included here) prove that the error term in Theorems  and  is in fact O(r 4 ).

tel-00557589, version 1 -19 Jan 2011

Proof. By Theorem  the barycenter o of B r has local coordinates -→ po = (o(r 2 ), o(r 2 ), Hr 2 4 + o(r 2 )). On the other hand -→ pp ′ is proportional to the normal to the regression plane, v. Thus by Lemma  -→ pp ′ = λ(O(r), O(r), 1 -O(r)). To compute λ, we use the fact that p ′ is the projection on the regression plane of p, and that o belongs to this plane by definition. This implies that -→ pp ′ ⊥ -→ op ′ and therefore

λ 2 O(r 2 ) + λ(1 -O(r))(H r 2 4 + o(r 2 ) + λ(1 -O(r))) = 0, which yields λ = Hr 2 4 + o(r 2 ) and therefore -→ pp ′ = (O(r 3 ), O(r 3 ), Hr 2 4 + o(r 2 )) = Hr 2 4 n + o(r 2 ).

. The discrete algorithm

The previous theorems assume that the surface is a uniform Lebesgue measure.

A constant sampling density is therefore necessary. This constant density will be approximated on discrete data by weighting each point by a weight inversely proportional to its initial density. More precisely, let p be a point and B r = M s ∩ B r (p). Each point m should ideally have a weight 0 ≤ w(m) ≤ 1 such that m∈Br w(Q) = 1. This amounts to solve a huge linear system. For this reason, we shall be contented with ensuring m∈Br w(Q) ≃ 1 by taking w(p) = 1 ♯(Br) , as proposed in [UH]. Let o be the weighted barycenter of this neighborhood. In R 3 , the coordinates are written with superscripts, e.g. the coordinates of a point u are (u 1 , u 2 , u 3 ). Thus, for i = 1, 2, 3,

o i = 1 m∈Br w(m) m∈Br w(m)m i . The cen- tered covariance matrix Σ = (m ij ) i,j=1,••• ,3 is defined as m ij = m∈Br w(m)(m i - o i ) • (m j -o j ) for i, j = 1, 2, 3. Let λ 0 ≤ λ 1 ≤ λ 2 be the eigenvalues of Σ with corresponding eigenvectors v 0 , v 1 , v 2 . For k = 0, 1, 2, λ k = m∈Br w(m) (m -o), v k 2 . (.)
Each eigenvalue gives the variance of the point set in the direction of the corresponding eigenvector. Since v 1 and v 2 are the vectors that capture most variations, they define the PCA regression plane. The normal v to this plane is the direction v

minimizing m∈Br w(m) (m -o), v 2 .
Effectiveness of Theorems  and . Both Theorems permit a priori to implement the mean curvature motion on the raw data point set without any previous orientation. Nevertheless, the numerical application of these theorems depends on the tel-00557589, version 1 -19 Jan 2011 assumption that a uniform Lebesgue measure on the surface is well represented by a uniform sample density. This is not true for the barycenter method of Theorem . Iterating the barycenter method with a small neighborhood and a slightly varying sample density leads to a local clustering of the samples. Indeed, the barycenter method provokes a normal motion, but also a non negligible tangential motion to the surface. This motion is precisely the one used in the Mean Shift method [Che] for data clustering. This undesirable clustering effect is illustrated in fig. Back propagation A normal motion by mean curvature can be defined for every point p 0 on the initial surface as a solution of (.) ( dp dt = Hn) considered as an ordinary differential equation with initial point p 0 . Thus, the backward scale space is trivial, provided the forward MCM implementation actually implements the evolution of each raw data set point p 0 . Let us consider a point p t and its evolution p t+1 at steps t and t + 1. Now, we can build the sequence d p (t) = p t+1p t and the reverse scale space operator P -1 t (p t+1 ) = p t+1d p (t), this operator allows to go backward in the scale space evolution from step t + 1 to 0. This is exactly the construction proposed in [PKG]. If we only need to go from step t to the initial data 0, without any intermediate step, the operator is even simpler to build, since we only need to store for each point p t its initial position P -1 t (p t ) = p 0 . This reverse

.. Higher order regression surfaces

The authors of [CP] proved that a degree n polynomial fitting estimates all k th order differential quantities to accuracy O(h n-k+1 ). In [PKG] an order  Moving Least Squares (MLS) method projecting the point onto the locally fitted least squares surface was actually proposed as a scale space operator. Yet these iterated projections cannot be consistent with MCM, because by definition they do not evolve degree two surfaces. Furthermore, the first step of MLS is always to compute a regression plane, which gives the estimate of the normal. We have just shown (Theorem ) that this computation by itself is sufficient to get directly not only the normal, but also the surface motion by mean curvature.

Can iterated MLS give a better estimate of the curvature than the projection filter? Comparative experiments were performed on a randomly and uniformly sampled sphere with added gaussian noise. The point samples were filtered four times by T r . By Theorem , each filtering step gives an estimate of the mean curvature. The same sampled sphere was filtered by MLS, the surface mean curvature being computed as the mean curvature of the approximating surface at each point. This estimate is very exact, since the difference on a C 4 surface between a point and its MLS estimate can be proved to be O(r 4 ). Both mean curvature estimates are compared by their mean and standard variations in the table of fig. .. The result shows that when the noise level increases the planar projection yields a much more stable computation (see the fast decay of the standard variation for the curvature estimate). This experiment is also coherent with the MCM consistency theorem. Indeed, the planar projection yields a point set with (slowly) increasing curvatures (once the noise is removed, i.e., once the standard variation is stable).

Fig. .

 is another illustration in D of the same phenomena: a circle with radius 1 and added gaussian noise with variance 0.05 is denoised by iterated T r and by an iterated MLS projection using the same neighborhood radius. In D, T r becomes a simple line regression and the MLS surface a degree  polynomial curve. The simplest MLS method is used: it merely performs a weighted least squares polynomial regression on the local neighborhood. The neighbors weights are equal to G(d) where G is a gaussian and d is the distance between the neighbor and the center point. The standard variation of G is equal to the neighborhood radius.

tel-00557589, version 1 -19 Jan 2011

Noise 0.01 0.05 0.1 Plane 1 1.00/1.95 1.15/5.57 1.27/4.76 Plane 2 0.99/0.07 1.02/2.16 1.17/4.89 Plane 3 1.00/0.02 1.01/0.16 1.05/2.10 Plane 4 1.01/0.01 1.01/0.05 1.02/0.27 Plane 10 1.04/0.01 1.05/0.01 1.09/0.04 MLS 1 0.94/0.22 0.11/2.58 -0.42/2.99 MLS 2 1.01/0.13 1.02/0.49 0.62/1.36 MLS 3 1.01/0.10 1.02/0.36 1.06/0.68 MLS 4 1.00/0.08 1.02/0.30 1.05/2.19 MLS 10 1.00/0.04 1.01/0.14 1.02/0.74 Figure .: Comparison of mean curvature estimates on a noisy sphere with radius 1 (mean/standard variations) given by iterated planar projection and iterated MLS regression. The curvature is evaluated at all points as the displacement along the normal induced by the planar projection (as stated in Thm ) for the planar case, and by the explicit computation of the MLS surface mean curvature in the MLS case. The same radius is used for both iterations and both regressions.

. First application: scale space raw data point orientation

Given an initial non oriented raw point cloud the surface orientation is a much needed information before meshing. The eigenvector of the least eigenvalue of the local covariance matrix is a classic approximation of the normal. We must then pick one of two possible orientations, and this choice must be globally coherent.

The idea is to start by picking a random orientation for one point and to propagate it to the neighboring points. Now, sharp edges or a noisy surface could fool such a propagation. If, however, the surface is smooth enough, the propagation of the normal is safe. Thus the overall technique to orient the raw data set will be to smooth it by the scale space, to orient the smoothed surface, and to transport back this coherent orientation to the initial data points.

The first tool to realize this program is a simple propagation method for a point p whose neighborhood B r contains some previously oriented points. The orientation is transmitted from one point to the next if their normal directions are similar, tel-00557589, version 1 -19 Jan 2011 i.e., their angle is below a certain threshold (algorithm ).

Algorithm : OrientateFromNeighbors(p,r,t)

Data: p an unoriented point, a threshold 0 < t < 1, a radius r, the set B r of p's neighbors within radius r Result: true if the point was oriented, false otherwise The input parameters for the scale space orientation (algorithm ) are the radius r and a threshold 0 ≤t≤ 1. Steps from  to the end are necessary because adding neighbors of points to the stack might not be enough to cover the entire cloud due to sampling irregularities. Once this procedure is over, there might remain non oriented points. These points are usually isolated points, and we choose to ignore them. In all our experiments the number of remaining non oriented points was tel-00557589, version 1 -19 Jan 2011 below 0.1%. At each step the radius is multiplied by an α > 1 factor. In step , the radius r is changed. Thus all normals are not computed with the same radius. This is why we must reverse the scale space to come back to the original point cloud. At scale 0, the normal direction is recomputed by local PCA for all points and the chosen orientation is the one which has positive scalar product with the previous normals. This is an application of the scale space paradigm, where the information is computed at a coarse scale and propagated back to the finest scale.

 Compute p's normal direction n by local PCA;  n ← normalized mean of already oriented neighbors' normals;  if (n • n) 2 > t then  if n • n > 0 then  n(p) = n;  else  n(p) = -n;
A result of this orientation algorithm is shown on a car point set (Fig. .).

. Second application: scale space meshing

We now discuss how to build a mesh on the raw point cloud. Direct meshing is not possible because of the surface oscillation due to fine texture and noise (see

Fig. .(c)).
The idea is again to perform meshing on the smoothed surface and to transport this mesh back on the original point cloud. An efficient triangulation technique, the ball pivoting method [BMR * ] is used in all experiments. The crucial faithfulness requirement is that the final vertices of the mesh must be an overwhelming majority of the raw data set points. This conservative requirement, incompatible with level set methods ( [KBH], [HDD * ], [LC]) is described in Algorithm . puting the octree to sort the points. Indeed the root of the octree is the bounding box of all points. Let us call L max the length of its largest side. Then, each cell represents a D cube with size L max /2 d where d is the depth of the cell. Counting the number of points in that cell gives an approximation of the number of neighbors of a point contained in this cell for a spherical neighborhood of radius r d = L max /2 d+1 . Performing this approximation in all non empty cells at the same depth gives an approximation of the number of neighbors for spherical neighborhoods with radius r d . The projection filter requires at least three neighbors per point to estimate a regression plane, but a robust estimate is experimentally attained with 30 neighbors. Of course, since the same radius is used for all points, it may occur that there are not enough neighbors to perform the plane regression. Those points must be eliminated, but in all the experiments less than 0.1% of points were removed this way. These points are mostly outliers, or isolated points in folds of the acquired object. Although their relative number is low, nonetheless this represents some thousands points that are eliminated.

Once the minimal number of neighborhood points has been fixed (and it has been fixed once and for all on all experiments to ), the radius is also fixed and the tel-00557589, version 1 -19 Jan 2011 meshing scale space only depends on the number of scale space iterations. When setting the radius automatically as described above, it was found that four iterations were always enough to smooth the point cloud and build the mesh. Thus, the scale space is conceived as a very local motion securing a reliable tangent space. In all experiments the points barely moved (less than 40µ for the Tanagra point set). The scale space and the ball pivoting reconstruction use the same ball radius. The parameter of the Poisson reconstruction (namely the octree depth) was set to be the largest allowed by our computing equipment (namely a 8 3Ghz processors computer with 48 Go RAM).

Transporting back the connectivity information (step 3) can in theory lead to a self intersecting mesh. Indeed, if two points lie too close to each other they may "switch position" in the scale space iterations, leading to a complicated surface topology. This problem can be solved by detecting all pairs of intersecting triangles. Then any remeshing algorithm can solve the problem by switching edges in quadrilaterals. However, this additional step was not implemented for two good reasons. First, the existence of a few intersecting triangles would be no serious visual inconvenience. Second, no such crossing was found in any of about twenty experiments on very large data point sets.

. Comparative experiments on high resolution data

The algorithms were devised for highly accurate point clouds acquired by a laser scanner. A typical example of the scanned objects is a mould of a fourth century B.C. Tanagra figurine acquired at the Museum of Cycladic Arts, Athens (Fig.

.(a)).

It is cm high and the point cloud contains 6 • 10 6 points.

Thanks to a very accurate calibration of the laser scanner device, the output is a well registered non oriented point cloud containing a negligible warp. Tests were also made on objects of the Stanford Fragment Urbis Romae database. In that case a registration of the raw sweeps is needed to have a point cloud representing the whole object. Preferring not to address the sweep registration problem in this chapter, we will use single sweeps for our meshing experiments and show that considerable texture information can be recovered from each sweep. Figs. . and . show the application of scale space meshing with a mesh rendering at fine and coarse scale. We can see on Figs. . and . that the surface texture is lost at a coarse scale, but completely and accurately recovered by scale space meshing. Comparing the back projected mesh to the result of a direct meshing of the initial samples (Fig.

.)

shows that the scale space triangulation is much more precise. In fact, a direct meshing is not applicable. It creates, among other artifacts, many spurious triangles. T r has been proposed as the simplest smoothing operation implementing the mean curvature motion. It may be objected that the surface could also be tel-00557589, version 1 -19 Jan 2011 directly approximated by the classic order  moving least square method (MLS). The most objective way to compare T r and M LS2 was to implement them with exactly the same neighborhood radius. . shows a comparison between the reconstruction obtained by the VRIP reconstruction method (see [CL]) and scale space meshing. The scale space method produces a significantly more precise mesh, as can be seen on the close up of GG]) and RIMLS ( [OGG]). APSS builds an implicit function by evaluating the distance between each evaluation point and an algebraic spherical fit of the surface. Though the method is not explicitely deviced for meshing, the isosurface can be extracted using the marching cubes. RIMLS is another modification of the standard MLS procedure. It is based on minimizing an objective function that gives less weight to spatial and normal outliers (i.e., sparse points and features). Here, marching cubes are explicitely mentioned for extracting the surface. For both methods, the resolution depends on the marching cubes grid resolution: it was set so that increasing it would not change much the visual aspect. Though both methods are visually close to scale space meshing, our method is much simpler and does not use an isosurface extraction at all. It is also the only method which preserves input samples and does not add additional vertices (both APSS and RIMLS use more than twice the number of input samples). Another problem is that the isosurface extraction by marching cubes introduces strong artefacts which are avoided by scale space meshing (Fig.

.).

Fig.

. displays the many acquisition holes at the bottom of the Tanagra figurine, in the folds of the tunic or near the right foot. By the scale space meshing these holes are not filled in and can be detected. Since the ball pivoting algorithm is used for triangulation, no triangle larger than a given threshold has been created. Indeed, to form a triangle, three points must lie on a sphere of given radius r. Thus, low density areas are considered as holes. simple shapes. Test point sets were built by sampling perfect geometric shapes (for example a sinusoidal surface). The root mean square distance of the triangle barycenters of the mesh to the real surface were compared for each meshing method. This distance is computed by the Newton-Raphson method. The first surface "Wave " has equation z = 0.2 cos(5x), "Wave " has equation z = 0.2 cos(5x) * cos(5y), the third surface is a regularly sampled sphere and the last one is a sum of two close and narrow Gaussians z =exp -(x-0.1) 2 0.01 exp -(x+0.1) 2 0.01 . The RMSE results are shown in the table of fig. .. It is obvious from these results that the Poisson reconstruction or any level set method cannot be applied to recover a surface with very thin details. On shapes containing no sharp edges, direct BPA and scale space meshing perform comparably. On the thin structure created by adding two very close Gaussians, the loss of precision of BPA is clear. This phenomenon is similar to the one observable in

. Complexity analysis and computation time measures

One scale-space projection requires the following operations: look for neighbors within radius r, build their covariance matrix and their centroid, perform PCA of this 3 × 3 covariance matrix. Therefore, once the neighbors are found, they are sequentially scanned in order to build the covariance matrix and the centroid. This yields 6 multiplications and additions per point for the covariance matrix update and 3 additions per point for the centroid update. The PCA complexity does not depend on the number of neighbors: it requires 9 operations. Knowing the least eigenvector, the projection is only 12 operations. There is one list scan (9 operations per processed point) and 21 operations once the covariance and centroid are built.

Assuming we have 30 neighbors, this yields a total of 200 operations per point. Finally finding the neighbors in an octree structure is O(logN ) (average) and one scale space iteration therefore is O(N (logN + 200)) operations, where N is the total number of points in the point cloud.

The computation time needed for meshing the Tanagra point set with six millions points was as follows: Sorting the points in the octree takes 1.2s. The scale space iterations requires 3 min. This leads to a total computation time of 19min for the scale space point cloud orientation and of 27min for the whole scale space meshing on an 8 3Ghz processors computer with 48 Go RAM. The maximum memory usage was less than 2Go. These figures should be compared with the time required for directly meshing the oriented point set by the ball pivoting method without any scale space iterations, which took 25min. Therefore, only a two additional minutes were used to get a much more faithful mesh.

tel-00557589, version 1 -19 Jan 2011

. Conclusion

The increasing accuracy of D triangulation scanners requires an effort to reconsider the whole rendering chain, and to obtain high quality visualization, actually better than what is obtained by photography. The present chapter has proposed a strategy to mesh the raw original surface, therefore ensuring a faithful rendering. Future work will test a closed scanning loop with our experimental scanner to scan the detected holes. Abstract: For each scanned object D triangulation laser scanners deliver multiple sweeps corresponding to multiple laser motions and orientations. The problem of aligning these scans has been well solved by using rigid and, more recently, non-rigid transformations. Nevertheless, there are always residual local offsets between scans which forbid a direct merging of the scans, and force to some preliminary smoothing. Indeed, the tiling and aliasing effects due to the tiniest normal displacements of the scans can be dramatic. This chapter proposes a general method to tackle this problem. The algorithm decomposes each scan into its high and low frequency components and fuses the low frequencies while keeping intact the high frequency content. It produces a mesh with the highest attainable resolution, having for vertices all raw data points of all scans. This exhaustive fusion of scans maintains the finest texture details. The method is illustrated on several high resolution scans of archeological objects.

. Introduction

Recent high precision triangulation laser scanners can scan surfaces of medium size objects with a precision of less than 10µ. Yet, although each scan has a very high precision, this precision can be lost again when merging multiple scans and meshing them together. This loss of precision entails a loss of visible texture, which explains the smooth and glassy aspect of most rendered scanned objects. On the other hand the merging of the multiple scans (often called super-resolution) is absolutely necessary. A patch of the object may well be acquired tens and even hundreds of times on well exposed parts. Indeed, many sweeps with varying trajectories are necessary to acquire the less exposed parts of the object. The main goal of the merging considered here is not to gain more detail and texture or to denoise the data point cloud by super-resolution: recent triangulation scanners yield scan sweeps with excellent quality. Unfortunately this quality is at risk of being damaged by the merging procedure itself. Thus, more trivially, the goal is to secure that the texture of each scan is not lost again due to slight matching errors which force a smoothing before a joint meshing. provokes strong tiling and aliasing effects, due to very small local offsets between both scans, in spite of the fact that they have been globally well registered. The challenge is therefore to merge both scans in such a way that the rendering quality does not decrease. The numerical problem is made more acute by two facts. First, not just two, but up to hundred scans may overlap in some region. Second, scans tel-00557589, version 1 -19 Jan 2011 boundaries appear everywhere, as illustrated in fig. . and make the fusion near these boundaries still more problematic. Each point of each scan has three-dimensional coordinates given either in a global coordinate system if the acquisition device is calibrated, or in a local coordinate system if the device is not calibrated. In the case of non-calibrated devices, the scans must be registered in a common coordinate system, and the registration problem becomes a rigid transform estimation. This problem has been widely investigated and has found efficient solutions [BM] [RL]. Yet if the scans had some internal local warping (which is usually the case), the rigid transform framework is not sufficient. A whole theory of non-rigid scan registration has therefore been developed [BR], [BR]. If the acquisition device is well calibrated the delivered scans are well registered, up to a given precision. Yet, as we already mentioned, a tiny residual mismatch can provoke strong artifacts similar to aliasing patterns (see fig. .) and forbids a direct meshing of the union of all data point clouds. This problem is generally solved by applying a method which meshes an implicit zero level set of a distance function to the raw points. The distance function is approximated by its Fourier coefficients [Kaz] or by radial basis functions [KBH]. The problem is that these methods result in a serious loss of accuracy when the final result is compared to each scan separately.

This chapter experiments on sets of scans of an object that have been either previously optimally registered by rigid or non-rigid methods, or registered through a high precision calibration of the acquisition tool. To demonstrate that no texture content will be lost, the goal is to mesh the entire point cloud. This means that all raw acquired points of all scans will ideally be vertices of the mesh. This requirement guarantees a complete preservation of all the acquired information, including noise and fine textures. Of course such a mesh is not numerically economic, but it is necessary for two goals: first to get high quality rendering of complex shapes such as archeological objects, and second to precisely explore all remanent artifacts such as the holes, inherent in any scanning process. For scanning control purposes, it is anyway quite rewarding to be able to see exactly what has been scanned.

. Previous Work

.. Rigid Scan registration

When the scans contain no warp, the registration problem sums up to estimating a rigid transform between scan coordinate systems by minimizing the distance between the reference scan and the transformed scan. The estimation of the transform in the quaternion form was proposed in [Hor]. The Iterative Closest Point (ICP) Registration procedure introduced in [BM] and [CM]

tel-00557589, version 1 -19 Jan 2011

from one scan to the closest point of the other scan and computes the transform based on these matchings. ICP was very successful and many variants were introduced ( [RL], [BS]). A study of the optimal sampling for the ICP algorithm was proposed in [GIRL].

Since ICP converges to a local minimum, the initial scan position must be close to optimal. Thus, some robust initial matches are needed to initialize the algorithm. This robust match search has also been investigated. The common idea is to find features easily identifiable on the scans (usually linked with highly curved points) and to match them. Spin Images, introduced in [JH], represent potential feature points by images. The spin image associated with point p is the function that maps a D point to (α, β), where α is the distance to the tangent plane and β the distance to the line parallel to n and passing through p. A linear correlation coefficient is used as a similarity coefficient between images to find the correspondences. In [SA] another feature carrier is introduced. It is based on computing geodesic circles around a point and projecting them onto the tangent plane. This gives a 1D 2π-periodic function parameterized by the angle. Feature matching is done by sampling the contours and computing a similarity measure between the contours. Other popular descriptors were described in [YF], [VSR], [ZH] or [KFR].

Another approach was proposed in [AMCO]: all coplanar -point sets that are approximately congruent are extracted on both shapes and matched using the fact that the distance ratios relatively to the intersection point are invariant to rigid motions.

Rigid scan registration methods assume that the scans must fit perfectly using only rigid transforms. Yet, if the scans have some warping, the method does not apply. Another problem is that the registration error might cumulate when registering multiple scans (see [BR] for examples of bad registrations). These considerations led to modeling the registration transform by a non rigid transform as will be shown in the next subsection.

.. Non-rigid scan registration

The method for non rigid registration using thin plate splines [BR] first applies a hierarchical ICP to find good features: the source mesh is iteratively divided through the middle of its longest axis and each half is realigned separately. This yields good feature correspondences at the cost of substantial discontinuities in the source mesh. These point matches are then used to compute the thin plate spline that best approximates all pairs of points. The spline being continuous, any discontinuity introduced in the scan splitting process is removed by the spline approximation. Once the scans are registered they must be merged. This is done using the VRIP method [CL], which will be described below.

An extension of this method was introduced in [BR]: the same thin plate tel-00557589, version 1 -19 Jan 2011  Chapter . High Fidelity Scan Merging spline non rigid deformation model is applied using an improved earlier point matching: features are found and matched in a process which rejects outliers. Other methods include [CR] where the non rigid registration is done using as input soft assignment between point pairs. This yields a functional minimization comprising a term of soft assignment.

.. Super-resolution from several scans

Recently, the problem of achieving super-resolution from multiple scans has been raised, mostly for improving range image resolutions ( [KMA]) where various low resolution scans with the same depth direction are acquired and registered by ICP. Depth values at each point of a high resolution grid are then interpolated with depth values of points falling in the neighboring cells. [AKSA] used rotating scans around the depth axis to build for each orientation high resolution range images. These high resolutions range images were registered and combined by considering the image gradient and the angle between the baseline and the image gradient to weight each point. Other methods building hybrid scanners to achieve high resolution include [NRDR] where positions and normals are acquired and used to improve the resolution.

Once scans have been computed and registered, a mesh must be built to allow a fast visualization of the surface model. The next section reviews methods to build the mesh.

.. Meshing

Meshing methods can be divided into two categories: methods that approximate the point cloud and methods that mesh directly the point set. Approximating methods usually build a function defined on R 3 whose 0-level set is the shape surface. A mesh is then built on the 0 level-set by the marching cubes algorithm [LC]. These methods include [KBH], [Kaz], among others. A very interesting variant of these methods is the VRIP algorithm ( [CL]). VRIP considers an implicit function taking into account not only point positions but also their reliability. Nonetheless, two drawbacks common to the mentioned methods are the automatic filling in of holes, and the implicit low pass filtering performed by the level set method. These methods usually compute the distance to the surface as an average of the signed distance of the point to its k-nearest neighbors. Thus initial points are discarded and de facto replaced by local averages. This removes noise in the cloud, but also loses fine details and textures.

Direct meshing methods include [ABK] or [AB]. We shall use the incremental ball-pivoting method [BMR * ], which is fast and does not fill holes. The method is based on pivoting a ball of fixed radius r around edges. Three points tel-00557589, version 1 -19 Jan 2011 are triangulated if they lie on a ball with radius r and empty interior. The ball is then pivoted around all three edges of the triangle, until it meets a point and has still empty interior. If no such point is met then the edge is the boundary of a hole. The parameter r is a bound for the creation of triangles: no triangle edge can have length above 2r. Thus low density areas are considered as holes.

In short, the dilemma is as follows: the approximating methods are not sensitive to a slight registration error, but can loose detail and texture. On the contrary, aliasing-like artifacts become visible when a direct meshing method attempts to keep all raw points of several scans. The best choice is to preserve the raw points and to apply a direct meshing method. But this requires eliminating all traces of inaccurate registration.

. Scan Merging

.. Principle

The general idea behind the scan merging method experimented here is to preserve point positions in non overlapping areas, and to make a fusion of the scans on overlapping regions while keeping all raw points. The fusion involves a smoothbase/height-function decomposition for each scan. The decomposition of a surface as the sum of a smooth base and of a height function was proposed for a different purpose in [KST], and [ZTS], where the height function was used to segment the mesh and extract features as contours of the height function. The underlying idea is that a surface S can be decomposed into a smooth base B and a height function h, so that:

S = B + h
B can be seen as the low frequency surface and h can be seen as the high frequency term. Given several surfaces

S 1 = B 1 + h 1 , S 2 = B 2 + h 2 , • • • , S N = B N + h N ,
the idea is to fuse the bases, but to keep exactly the h i terms, thus preserving all fine details. In other terms, a common basis B for all surfaces must be found, the high frequencies of all scans adopting this common basis thereafter. This strategy is comparable to the one used for morphing applications in [PKG]. In a way, the idea is similar to [SCOT], where the high frequency error due to quantization was transformed into a low frequency error much less noticable. The data merging using a high/low frequency decomposition has long been a classic method in image processing [BA]. This article introduced the idea of separating each image into various frequency bands by a Gaussian pyramid. The low frequency bands were merged separately to obtain a smooth blending of different images. The method has been successfully used to create panoramas from multiple images [BL] and texture D models [Bau]. Two major differences are tel-00557589, version 1 -19 Jan 2011 that in [BA] all frequency bands are merged, whereas the method proposed here only merges the low frequencies while keeping the high frequencies intact. Another important difference is the usage of a nonlinear heat equation instead of a linear frequency decomposition.

The next section addresses the robust decomposition of a surface into a base and a height.

.. Low/High frequency surface decomposition

Since the pioneering article [Taub] it is known that mesh high frequencies are removed by the application of the intrinsic heat equation ∂p ∂t = ∆p. Yet, our scanned surfaces are given as point clouds and not as meshes. A numerical scheme of the heat equation for raw point clouds must be used. This question has been addressed in [BSW] and [PKG]. We shall use the simple implementation of the intrinsic heat equation proposed in [DMMSL]: this chapter proves that the intrinsic heat equation can be implemented by iterating a projection of each point onto the local cloud regression plane. Consider the projection operator T r that projects each point p onto the regression plane of the neighbors of p enclosed in a ball of radius r. Then it can be proven that this motion is tangent to the intrinsic heat equation. The iteration of T r yields a scale space (a representation of the shape at various smoothing scales). In all experiments r is set so that the ball B r centered at p contains about 30 neighbors at almost all points, and the number of scale space iterations n is set to 4. The first parameter () is fixed so that a reliable regression plane is always computed. The second parameter, namely the number of iterations , is chosen to guarantee a smooth enough basis in all cases. It can be increased without damage. When iterating the projection operator with an initial surface S 0 , the surface S t is iteratively smoothed. To each point p t of S t corresponds a point p 0 of S 0 , and the height function can be taken to be the vector h(p t ) = p tp 0 .

An alternative definition for the height would be the scalar function h(p t ) = (p tp 0 ) • n, where n is the normal to S t at p t . Yet, the results with both height variants being fairly identical, the simplest definition was kept: it separates each data point into a smooth base point and a high frequency vector.

.. Finding a common smooth basis for all surfaces

Choosing a common basis for all scans is the next question. A natural constraint on the method is to keep fixed the points belonging to regions where only one scan is available. Finding the common basis then becomes straightforward: It is enough to apply the same number n of iterations of T r with the same parameter r to all the sets after they have been put together. This global filtering assumes that the high tel-00557589, version 1 -19 Jan 2011

.. Algorithm

The method is summarized in Algorithm . The algorithm is based on two applications of the intrinsic heat equation scheme (here the iterated projection on the regression plane) with the same parameters and the same number of iterations. All registered scans are given in the same global coordinate system. The first application (Line ) is done on the separate scans yielding the intrinsic high frequencies of each scan. The second application (Line ) is done on all scans together. When filtering all scans together (lines  and ) the registration error is suppressed and we get a common low scale registration or basis, the set of points b(p). Adding back to them the high frequency component ---→ b i (p)p restores all details from all scans.

An important feature of the method is that each region A of the shape that has been acquired by one scan only is not altered. Indeed, inside such a region, applying the separate scale space or the common scale space is strictly equivalent, since there is only one scan in the neighborhood of the points of A. Then the point is first filtered to b i (p) = b(p), and therefore moved back to its original position p at Line . So in areas with only one scan, point positions are not changed. The only effect of the algorithm is the merging of overlapping scans.

.. One-dimensional study

It is easy to illustrate the method in 1-D on simple 1D shapes. Our goal was to check that the proposed method superimposes two simulated scans without any smoothing effect. To do so, two noisy straight lines A and B were synthesized from the same model and then merged by the algorithm. The noise of each set A, B, A ∪ B was estimated as the root mean square error to their regression lines before and after merging. The results in Tab. . show that the merging did not cause any denoising. Indeed, the RMSE does not decrease by the merging procedure.

. Implementation and Results

The inputs of the merging algorithm are the outputs of our laser triangulation scanner. This device being accurately calibrated, the scans are in principle already registered so that no extra software registration is needed. Nevertheless, the ICP algorithm is applied to see if it can remove the aliasing and tiling artifacts. In this case ICP failes: the positions computed by ICP oscillate around the input scan positions, and the resulting meshes are no better. The registration process was implemented on a 1.5 Ghz processor with 48GB RAM. An octree structure was first built to allow for fast access to the neighbors of each given point. ,. and . present the results on these data. For all point sets, two different renderings are displayed: the first one is a ball pivoting [BMR * ] meshing of all raw scan points without any merging. The scans were preregistered by the calibrated acquisition device and no software post-registration was needed. The second rendering is again a ball-pivoting meshing, but applied to the merged point set. The rendering was made using the POV-RAY ray-tracer. The conclusion is common to all experiments: even if the scans are actually very accurately registered, the tiny warps of the grids always create some aliasing visible as grid or tiling effects. After the merging procedure (which only slightly affects the low frequencies), these undesirable effects disappear almost completely. In the procedure more than .% of the raw points were kept. Thus, the final result indeed is highly faithful to the raw scan. Yet a careful attention shows some remains of aliasing (Fig. ., last column). The area of these is actually small, being inferior to the area of the holes. They could easily be removed by a selective local smoothing. Some of the bigger pieces, like Nefertiti, show no defect at all. Comparison To better judge the texture preservation, the rendering of a scan alone (ground truth) was computed and compared to the rendering of all scans in the same region on Fig  .. This shows that the visual information loss after scale space merging is very low compared to the one due to a simple joint meshing.

tel-00557589, version 1 -19 Jan 2011

It is crucial to compare the raw merging method results with results obtained by the level set reconstruction method of the unmerged scans point set. The result of the level set method applied to the Tanagra head (fig. . b), obviously introduces an important smoothing and loses texture in comparison to the merging result (fig.

., a). But even with that smoothing the result still keeps several artifact lines due to the scan offsets: these offsets become visible at the scans boundaries. See the nearly straight long lines on the surface, mostly vertical and horizontal. It can also be asked if an efficient denoising method could actually restore the raw set. Fig.

.-c, shows the result of the application of the bilateral filter [FDCO] to the union of the scans. This iterated filtering was applied up to the point where aliasing artifacts were no more visible. Clearly, this entails a much too strong smoothing of detail and texture.

The scan merging is a very local method which is therefore computationally efficient (see Tab. .). Yet, if the input data are not already well registered the merging could obviously fail. The method corrects the slight misalignments only in the normal direction. A tangential drift in the original registration could therefore cause a loss of sharpness or a loss of small details. Nevertheless, this degradation seems to pass unnoticed. Indeed, for last generation triangulation scanners like the one used for the experiments in this chapter, the registration error is very small. For a point cloud with side-length 99mm the observed average point offset after merging was 0.081mm, with standard deviation 0.012. The tangential offset could not be measured. The explanation of the relative visual success of the method is that even a tiny normal offset causes a dramatic change in triangles orientation, tel-00557589, version 1 -19 Jan 2011 and therefore completely jeopardizes the visual quality of the triangulation. An equally small tangential offset seems to be visually undetectable. Thus, the merging method corrects the normal error, and makes the tangential error unnoticeable.

The proposed merging can be seen as a local non rigid registration. Therefore it can be compared to the result given by state of the art non rigid registration methods [BR]. To perform the comparison, the problem arose that the scans did not systematically contain strongly identified features. Most scans of the mask point set were simply rejected by the non rigid registration method described in [BR]. In order to perform a serious comparison anyway, two sweeps of the fragment u of the Stanford FUR project were used. The computation times were, however, considerably different: it took more than 2h30 to register non rigidly these meshes. On the same computer, using only the raw points and not the meshes, the merging took only 84s. The final meshes were built using Poisson Reconstruction [KBH] in both cases. The registration artifacts (two horizontal lines limiting the overlap area, fig .) are much less visible with the scan merging than with the non rigid registration.

. Conclusion

The main conclusion of the study is that it is possible to fuse multiple raw scans with minimal accuracy loss, provided an accurate previous registration has been performed. Future work will focus on the detection and handling of remaining holes, and on the automatic assessment of surface quality to replace a visual inspection. ) but also to the resolution variation. Indeed the borders of the object were acquired using multiple orientations while the middle was acquired using only one orientation. This is why we have such a precision difference that is not fixed by the algorithm.

tel-00557589, version 1 -19 Jan 2011

tel-00557589, version 1 -19 Jan 2011 Abstract: In this chapter, a method to fill holes in scale space meshes is presented. The method uses two key features of the scale space meshing method: first the fact that it permits to detect accurately the holes, second that it permits to work at a coarse smooth scale and then to transport back the result onto the meshed raw cloud. The goal is to fill holes so that the missing patches blend nicely with the rest of the object. The holes are first detected in coarse scale meshes, the filling patches are refined, and a mesh fairing method is applied so as to blend the patch in the surrounding mesh. The back propagating operator is eventually used to transport the patches to the fine scale mesh.

. Previous Work

The problem of filling holes in meshes is crucial because all D acquisition systems have a limited access to hollow parts of objects. Even a very careful scanning process leaves behind a number of holes corresponding to shadowed parts of the object. Thus many solutions have been proposed to fill in the holes. First, one should note that meshes built by approximation of a signed distance function ( [KBH], [Kaz], [ACSTD] and [HDD * ] among others) automatically fill the holes, so that no postprocessing hole filling technique is required. This is actually a problem:

The holes are not really detected. In consequence no exploration and no improvement of the quality of the restoration in these parts is possible. In agreement with the methodology developed in this thesis, given a mesh with holes the goal will be to fill in the holes without altering the other parts of the mesh.

The existing methods to actually fill in detected holes can be divided in two categories: the volumetric methods and the triangle patch methods.

.. Volumetric methods

These methods use ideas similar to the level set mesh reconstruction methods ( [KBH],[HDD * ]...) to determine in unknown areas (holes) a plausible surface.

In [DMGL], the problem of filling holes is treated by building a signed distance function, whose zero level set is the surface. The goal is to preserve the geometry when it exists and to create smooth transitions to plausible geometry in unobserved areas. The surface is first converted to a volumetric representation: a D grid of values of a signed distance function defined only in a narrow band around the surface. An associated weight function measuring the confidence in the distance value is built. The distance values are diffused from known areas to adjacent unobserved areas. Once the diffusion is complete the zero set is extracted via marching tel-00557589, version 1 -19 Jan 2011 cubes ( [LC]). The diffusion is done by convolution with a low pass filter (see also [GLWYT]). Many approaches based on a similar voxelization have been proposed. For instance [CBC * ] uses Radial Basis Functions (RBF) to approximate the distance and extract a watertight surface. The method in [VCBS] extends the inpainting technique over the voxel grid to "inpaint" the implicit representation of the surface and extract the level set. [NT] gives value 0 and 1 to voxels inside and outside the surface. Determining whether a cube is inside or outside the surface is done by projecting the shape onto k random planes. Each plane votes for the inside or outside position and a majority count yields the classification. In [Ju] the cube labeling is deduced from the definition of a dual surface. In [PR], besides the voxels labeling, cubes intersecting the surface are split in tetrahedra and the tetrahedra are also labeled. Then a graph-cuts algorithm is used to deduce the final labels.

In [ZGL], a method that mixes patch filling and volumetric methods is presented. First holes are identified and filled with triangles using the advancing front mesh technique. Normals of the newly created triangles are estimated and used as vector field for solving a Poisson equation and obtaining the new coordinates of every vertex. Positions are then updated and the result is a smooth based patch.

An advantage of volumetric methods is that complicated hole topologies are not a problem. On the contrary patch filling methods mostly consider holes that do not contain islands. The next section reviews patch filling methods.

.. Patch filling methods

Detecting the holes in a mesh is a trivial task as their boundaries are closed loops of boundary edges (edges adjacent to only one triangle). All methods start by such a detection. Triangulating a D polygon is a problem with many solutions proposed in the past years. Nonetheless there is no simple method to fill a D polygon with triangles, indeed, [BDE] proved that the problem of filling a D polygon is NPcomplete. It also proved that if a perspective projection of a -dimensional polygon P from some point onto some plane Q is simple then P is triangulable. An algorithm for finding a suitable projection plane is proposed. For all pairs of edges, if the projection direction is inside the sphere region that sees the two edges intersecting, then the direction is invalid. A region of admissible projections can then be computed.

Many heuristics have been developed to fill the D mesh holes.

[BS] fills cracks in meshes by matching points on the crack boundary one with another and then stitching the borders with triangles (see also [TL]). In [BNK], vertexvertex and vertex-edge correspondences are found on the hole boundary. The distance between the correspondences is used as an error to sort the pairs into a priority queue. The queue is processed in descending order: pairs are linked, new tel-00557589, version 1 -19 Jan 2011 correspondences created, and the priority queue is maintained up to date.

Since triangulating a D polygon is quite easy, the idea of projecting the polygon on a plane and triangulating has been extensively used. In [PSM * ], the surface is projected onto planes and standard image hole filling techniques are applied to the projection. Then the surface is back projected to D. Yet, a plane selection is possible because this method considers only range images. In [Jun], a method was introduced to divide a complex hole (a hole that self intersects when being projected on a plane) into several simple sub-holes. If two projected edges intersect, the intersection points on both edges are found and projected back on real  edges. These points are used as splitting points and connected, separating the hole in two parts. Finally the inner mesh is refined to obtain regular triangles. In [LYZ], a modification of [Jun] is presented. The PCA regression plane is found and the hole is projected onto it. If the hole is not self-intersecting, it is simply filled. Otherwise, intersecting edges are found, and a parallelepiped-like region around the intersecting edge and orthogonal to the neighboring triangle is built. Inside this region, the algorithm looks for boundary points and selects one of them to create a new triangle adjacent to the orthogonal edge. The selection is done for example by comparing the angles between the initial triangle and the created triangle. The new triangle provides a natural splitting of the hole. In [BWS * ], a different method to unfold the hole boundary on a plane is introduced. The hole unfolding aims at making 4 consecutive vertices planar. Thus, for a loop of n points v i , it minimizes the energy E = n-1 i=0 ∠t i , t (i+1) mod n where

t i = (v i -v i+1 mod n )×(v i+2 mod n -v i+1 mod n ).
The minimization is done subject to the constraint that the minimum distance between two segments of the boundary loop should remain above a threshold ε. Simulated annealing is used to perform the minimization. After the unfolding the obtained polygon is triangulated by constrained Delaunay triangulation. The mesh is embedded in R 3 by moving each of the boundary edge back to their original D position. The resolution of the inner mesh is then refined simply by standard Delaunay refinement algorithms. Inner vertices are added so that they minimize the area of the inner mesh. This is done by applying repeatedly a Laplacian smoothing.

Hole triangulation is not the only problem that arises. Indeed patches must blend nicely into the shape. For example in [PS], holes are filled using splines.

In [Lie] a popular geometric method for filling holes in meshes is proposed. Holes are filled by D polygon triangulation minimizing a given weight. Instead of using an area minimization weight, a weight that takes into account the dihedral angle between adjacent triangles is used. Thus, large dihedral angles are penalized. The patching mesh is refined so that its density matches the one of the surrounding mesh. Mesh fairing is applied to get a smooth enough surface.

In [PMV], a topological grid is inserted into the hole (by mapping the grid of a disk onto the hole), and the inserted mesh is deformed to satisfy blending criteria. A bar network coupled to a geometric model is used to deform the inner mesh while tel-00557589, version 1 -19 Jan 2011 keeping the outside mesh identical. The mesh is deformed so that it minimizes the curvature variation while satisfying additional user-supplied constraints.

Only one of the patch filling methods explicitly addresses the problem of holes with islands: [LMW] first determines the hole position and an "influence" area. If there is an island in the hole, the influence area and possible island is searched for feature points. Those feature points are linked accross the hole, creating feature curves that divide the hole into sub-holes. Each sub-hole is then filled. Yet the computation of the influence area is not explained.

Finally, Moving Least Square surfaces were also used to fill holes. In [WO], planes are fitted to the set of the hole boundary points and a height for each boundary point is found. Finally new samples are interpolated by MLS interpolation on the height field. The use of MLS guarantees that new samples will fit smoothly into the mesh. Yet if the hole is not planar then it may fail dramatically.

Extending the idea of image inpainting has also been considered. The idea is that hole filling patches should look like other parts of the surface. The next section reviews the few surface inpainting methods that have been proposed.

.. Inpainting for surfaces, hole filling based on similarity

A review of mesh hole filling can be found in [BPK * ]. In [SACO], the aim is to complete missing parts of the surface by transferring appropriately sampled and fitting regions of the shape. The surface is approximated by fitting a trivariate low degree polynomial. The space is divided into cells. A surface cell is called valid if it contains at least m points (there are enough points for surface representation in the cell) and invalid otherwise. For each invalid surface cell, the idea is to import and paste the content of a valid surface cell that matches the surface approximation in and around the invalid cell. The best matching cell is found by computing a vector valued signature for each cell and comparing them.

In [PMG * ], shape completion is made by combining geometric information from different context models. This is done by database retrieval, non rigid alignment of the candidate models and blending the aligned segmented models to get a consolidated model. For each missing region the best fitting model in the database has to be found. The shape is aligned to the shape database and information must be extrapolated from the models to region where data is missing. This is done by copying vertices from the model patches. The warping function is then recomputed for each model in the database to conform with the enhanced point set. After updating the alignment context patches are enlarged and the process is iterated until the hole is closed.

In [BSKa] and [BSKb], a multiscale inpainting method for surfaces is presented. Holes and a scale space is built from the input surface. By projecting the neighborhood onto its regression plane and resampling the neighbor's height on a tel-00557589, version 1 -19 Jan 2011 rectangular grid, a neighborhood descriptor is built. Thus the comparison between neighborhoods is easier. The -level inpainting is then defined as finding for every bounding point b an appropriate candidate point from the surface and to copypaste its neighborhood. Multi level inpainting is done by considering the next level filling surface as a guidance surface for the filling of the current level surface.

[VCBS] (see section ..) is also based on inpainting, yet the inpainting is performed on a implicit surface formulation.

.. Point clouds

Holes in point clouds have not been really investigated yet. A noticeable effort was made in [CJ]. First, hole boundaries need to be found. In this work boundary points are characterized as being much farther from the barycenter of their neighborhood than other points. Since this would also detect points on sharp edges, another criterion is added taking into account the point distribution on the local tangent plane. Then a triangulation of the boundary curve is found and a set of points sampled on this triangulated surface is generated. In [BSK], several criteria are combined to give to each point a boundary probability and extract boundary loops using this probability and coherence properties of the boundary.

. A hole filling algorithm

.. Filling the hole: finding an initial patch

In [Lie], an algorithm for finding the initial patch is proposed that aims at minimizing the angles between the patch triangles. Yet it requires storing a n 3 array for a hole with n border points. Some holes in our meshes contain more than 5000 border points, which makes this method intractable. Therefore, we shall use a simpler triangulation by projecting the contours on a D plane and detecting self intersections due to projections. This kind of triangulation can produce sharp angles between triangles. Better initial patch finding algorithms could be found, yet it is a first scale space mesh hole filling investigation. The proposed algorithm is summed up in Alg. . (The definitions of terms used in the algorithm are given below). Note that this algorithm works only for connected meshes. This means that the input mesh should contain a single connected component (the mesh is easily cleaned from remaining "islands" after singular triangles have been removed). This can lead to the loss of some samples.

A few explanations are necessary. Let a "border point" (resp. "border edge") be a point (resp. and edge) of a hole contour. Edges that are not border edges will be called inner edges.

tel-00557589, version 1 -19 Jan 2011

Algorithm : Mesh Hole Filling Data: A mesh with holes, each point is endowed with an oriented normal Result: A hole-free mesh • Singularities: (line ), all triangles having a inner edge with two border endpoints should be removed.

• Connected mesh: a mesh will be called connected if there is a path from one triangle to another composed of edge-adjacent triangles.

• Hole contour: Since we know the normal at each point position, we can easily find the triangle normal orientation. A contour point ordering is defined by setting that two successive points of the contour v i , v i+1 adjacent to triangle

t = (v i , v i+1 , v ′ ) with oriented normal ⋉(t) should verify: ( --→ v ′ v ′′ ⊗ ---→ v i+1 v i ) • n(t) > 0 where v ′′ is the projection of v ′ on (v i , v i+1 ).
• D self-intersection: (line ) In the case of manifold meshes, no edge can intersect with another edge except at the endpoints. In this case removing the singularities means checking that the hole border contains each point only one time. If this is not the case the hole is split into two holes at the intersection point.

• Filling a D hole: (line ) There exist many methods triangulating a D polygon. Here we used Seidel's method ( [Sei]). It simply finds triangles that cover the interior of the polygon. Triangle vertices are vertices of the boundary polygon so that no inner vertex is added in the initial patch finding step.

tel-00557589, version 1 -19 Jan 2011

• D self-intersection: (line ) When projecting a hole on its regression plane the boundary might self-intersect. Thus the hole should be split into parts that project without intersection. Call v i v i+1 and v j v j+1 the edges that intersect when projecting the hole boundary. Triangles (v i , v i+1 , v j+1 ) and (v j , v j+1 , v i+1 ). This two triangles split the hole boundary in two closed loops. At the end all the patches of the split holes, as well as the splitting triangles, are added back to the initial self intersecting D hole patch.

• Uplift: (line ) Uplifting the patch only means creating a D patch from the set of D triangles. Each D point is the projection of a D point, thus we simply create a D triangle whenever the vertices projections are linked by a D triangle.

.. Refining the patch

The created patches have very large triangles (since no vertex is added) and are not smooth. Therefore in a second step the patch is refined. The same refinement method as in [Lie] is applied. We briefly sum up this refinement in algorithm .

The idea is to continue adding inner vertices, while the density of the patch is below the density of the surrounding mesh. To do that a density value σ is computed for each point (called the scale attribute), σ(v) is simply the average length of all edges adjacent to v. Then for each triangle the centroid c is computed along with its estimated scale attribute σ(c) as the average of the scale attributes of the triangle's vertices.

Our experiments led us to consider α = 1.75 as a good value. As a matter of fact, we consider this value as fixed and will not change it anymore.

The algorithm makes use of a well-known algorithm which is the edge-swap or edge-flip method. Consider a patch composed of two triangles (pqr) and (pqs). Edge pq is adjacent to two triangles: (pqr) and (pqs). Flipping (swapping ) pq means destroying those triangles and creating (prs) and (qrs).

Relaxing an edge (Algorithm , lines  and ) means that for two adjacent triangles, we check that the non-mutual vertices are outside the opposite triangle circumsphere. If this is not the case and if an edge swap can fix this, the edge is swapped destroying two triangles and an edge and creating two new triangles and the swapped edge. This edge relaxation is necessary to avoid triangles with sharp angles (though there will be triangles like this near the hole boundary, since we do not add any additional vertices).

tel-00557589, version 1 -19 Jan 2011

 if σ(c) < α v m -c and σ(v m ) < α v m -c for m = i, j, k then  changed ← true;  Split the triangles into (c, v i , v k ), (c, v i , v j ) and (c, v j , v k );  Relax edges (v i , v k ),(v i , v j ),(v j , v

.. Giving shape to the patch

Finally the patch should blend nicely into the surrounding surface. Therefore a standard mesh fairing is applied to the patch inner vertices so that the final mesh is smoothed. [Lie] proposed to minimize a thin plate spline energy to obtain C 2 continuity. This approach was first proposed in [KCVS]. The idea is to minimize the energy:

E T P (f ) = f 2 uu + 2f 2 uv + f 2 vv .
Let the umbrella operator U be defined for each vertex v by

U(v) = 1 #N (v) v i ∈N (v) (v i -v)
where N (v) is the one ring neighborhood of v.

[KCVS] showed that this umbrella operator is a discrete analogon to the surface Laplacian.

Recursively this gives

U 2 (v) = 1 #N (v) v i ∈N (v) (U(v i ) -U(v))
tel-00557589, version 1 -19 Jan 2011



Chapter . Filling Holes in Scale Space Meshes

A way of solving for the thin plate energy minimization is to iterate the following

v ′ = v - 1 ν U 2 (v) where ν = 1 + 1 #N (v) v i ∈N (v) 1 #N (v)
. This leads to an inflation or deflation of the patch to match the normal of the boundary triangles.

.. Special holes

In some cases we do not want to fill the holes, in particular the base holes (the base of a statue for example). This is why in some cases we exclude the base from the set of holes to fill. This base is detected by finding the hole with lowest coordinate in one specified direction. In the case of our acquisition device, this direction was always z.

. Results on synthetic data

.. Sphere example

The following is an implementation of Liepa's method ( [Lie]) with the single difference that a different initial patch finding algorithm is used as discussed in section ...

. Scale space meshes

Scale space meshes are a special case. Indeed they have (at least) two resolutions: the coarse scale mesh and the fine scale mesh. At which scale should the hole filling take place?

Since the hole filling technique assumes a C 2 continuity, it seems wiser to apply the hole filling to the coarse scale mesh and then to use the back propagating operator to go back to the fine scale. Yet no such operator is defined on the hole filling patches.

Each patch point is consistently oriented with the rest of the mesh (this normal computation and orientation is done when adding the vertices). An estimation of the fine scale variations is easily obtained by the scale space decomposition: the amplitude of the first iteration motion gives exactly the local texture and noise of the shape. If we assume that the missing part of the surface do not include any sharp feature then we can just add noise to the patch points in their normal direction. The noise addition is an estimation of the backward scale space operator for patch points. It implicitly assumes that holes do not contain sharp features but only fine scale variations.

The amplitude of the directional noise is gaussian and it has the variance mean of the first iteration motion. This is enough for the patch to "blend in" the rest of the surface. Figures ., ., ., .,., .,., ., . show results of the hole filling method. It appears that it works well for small holes, yet is not good enough to fill large holes satisfyingly (see the bottom of the mask, fig. .). tel-00557589, version 1 -19 Jan 2011

Conclusion

This chapter is clearly exploratory and it is not published. Its goal was to make a first review of hole filling algorithms and to define a first method to detect and fill holes taking advantage of the scale space meshing method. Several improvements are still necessary. In particular, a better initial patch should be found. An important improvement would also be to handle holes with islands.

tel-00557589, version 1 -19 Jan 2011

. Introduction

Segmenting meshes or point clouds into their significant parts is a basic but still challenging problem. Robust shape descriptors are required for most applications such as facet classification, shape registration, mesh simplification or shape retrieval. To this aim, the "crest lines" detection plays a role analogous to edge detection in image analysis. Crest lines are usually defined as the loci of directional extrema of curvature. Ridge (resp. valley) points can be defined as points where the maximum principal curvature takes a positive maximum (resp. negative minimum) along the line tangent to its eigenvector. This notion being closely linked to the surface curvature, robust curvature estimators are needed. Curvatures can be computed by surface regression [CP], by discrete schemes using the mesh triangle geometry [MDSB], or by computing the surface tensor [Rus]. One of the most used curvature estimation method is [Taua], where curvatures are estimated by drawing curves on the mesh surface (see [LP] for an adaptation to point clouds). Other mesh based methods use the angles between adjacent mesh triangles to determine the curvature [DVVR], [HMG].

Most crest extraction methods detect and link the points where the derivative of the curvature crosses zero. These methods are usually defined on meshes, but can be adapted to raw point clouds [LFM], [BA], [YBS], [OBS], [DIOHS], [HPW]. In [DIHOS], potential features are extracted by regressing surfaces near those points, estimating the number of fitted surfaces, and deducing the feature type. The idea of using feature lines for surface segmentation was investigated in [SF] and suggested in [IFP]. Regression free methods include [GWM], where the analysis of the surface local covariance matrix leads to point classification, or [PKG] where a multi-scale approach is introduced, yielding good results for mechanical shapes. Indeed, texture features are detected at fine scale, whereas at coarse scale the features describe shape geometry and are more robust. Recent research has also made huge progress in the rendering of viewpoint tel-00557589, version 1 -19 Jan 2011 dependent apparent edges [DFRS], [JDA].

Crest lines have, however, two limitations. The first is that these lines are often obtained by computing degree three derivatives, which is not an easy task for noisy or textured surfaces. The second is that crest points must be connected by some heuristic linkage. The final crest lines being often open or broken, they do not provide a surface segmentation.

Analogous problems arise in image analysis with edge detection. Image analysis therefore also uses segmentation methods dividing the image into regions separated by closed curves. One of the most reliable ways to define such closed curves is to extract the contrasted level lines. The level lines are the boundaries of connected components of upper (or lower) level sets. They inherit from these connected components an inclusion tree structure (fig. .). This structure was called Tree of Shapes in [BCM], and a fast algorithm computing them, the Fast Level Set Transform is given in [MG].

For an image I defined on a domain D and with values in R the level sets with level λ are F λ = {x ∈ D|I(x) ≥ λ} (upper level set)

F λ = {x ∈ D|I(x) ≤ λ} (lower level set) If λ ′ > λ, then F λ ′ ⊂ F λ and each connected component of F λ ′ is contained in one connected component of F λ .
The set of connected components of upper (resp. lower) level sets partially ordered by inclusion is therefore a tree. The shape tree proposed in [MG] is a fusion of both trees.

The level sets can be represented by their borders ∂F λ and ∂F λ which are unions of closed Jordan curves, the image level lines. Several methods have been proposed to select the relevant level lines. A definition of meaningful level lines is given in [CMS], [DMM]. More recently the MSER method introduced the same objects with different names: the connected components of upper or lower level sets are called extremal regions (ER). The ones with best contrasted level lines are called maximally stable extremal regions (MSER) [MCUP]. The extraction of significant level lines to segment data is so useful that it has been extended to D medical imaging to extract meaningful level surfaces [CSA], [MZFC], and to video analysis [DB], where extremal regions are tracked from frame to frame.

To the best of our knowledge these level line techniques have not yet been extended to meshes. The reason could be the lack of straightforward intrinsic scalar functions linked to a mesh, (such as the grey level for images). But there are actually such functions on meshes, the simplest one being the mean curvature. Several methods have already considered the curvature level lines and the umbilical points, but mostly from a theoretical point of view [KNSS], [CPZ] and [MWP]. Until now, curvature level lines have not been studied as valuable feature lines, or used for surface segmentation. The goal of the present chapter is to describe an algo-tel-00557589, version 1 -19 Jan 2011 rithm computing all conspicuous curvature level lines, and to give experimental evidence that the method detects valuable mesh features.

In a way, the present work extends [KST] and [ZTS]. In these works, the surface is (implicitly) decomposed into a smooth base and a height value in the normal direction. Then the edges or iso-contours of this height are extracted. The Mesh-MSER framework considers this same situation in a more general setting : a surface with any scalar function defined on it. As shown on Fig. ., results comparable to the results of [ZTS] can be obtained by a significantly simpler and more general procedure.

The remainder of this chapter is divided as follows: Section . recalls the image MSER method, discusses its adaptation to meshes, and gives the algorithm building the level set component tree. Section . describes the algorithm extracting maximally stable extremal regions from this tree. Section . shows results on various simple and complex scanned objects, discusses strategies for level line selection, and compares D Mesh-MSER results to D MSER results on pictures of the scanned objects.

. Mesh Extremal Regions

.. Definition of MSER for D images [MCUP]

Let I be a real function defined on an image domain D ⊂ Z 2 . MSER needs an adjacency relation A for elements of I, and usually chooses a  or  connectivity. The boundary ∂N of any set N ⊂ D is defined as {q ∈ D \ N | ∃p ∈ N, pAq}. An Extremal Region N is a region such that for every p ∈ N and every q ∈ ∂N one has I(p) > I(q) (maximum extremal region) or I(p) < I(q) (minimum extremal region). To define Maximally Stable Extremal Regions (MSER), consider a sequence (N i ) i of nested extremal regions (N i+1 ⊂ N i ). A region N i * in the sequence is maximally stable iff its area change rate q(i)

= |N i-δ |-|N i+δ | |N i |
has a local minimum at i * . The small variation δ > 0 is a parameter of the method.

The detection of MSER proceeds by: . sorting pixels by intensity;

. iteratively placing pixels in the image and updating the list of connected components and their areas;

. selecting intensity levels that are local minima of the area change rate as thresholds, producing MSERs.

In a more formal way, the method uses the fact that upper level sets F k = {p|I(p) ≥ k} are ordered by inclusion:

F k+1 ⊂ F k .
F i such that i < k and N k ⊂ N i .
Thus the set of extremal regions is a rooted tree [BCM] called (upper) level set component tree. Its dual tree is the (lower) level set component tree. The fastest component tree method seem to be [NC]. It is its tree structure that allows the fast selection of MSERs [DB].

Recently the question of the interest of MSER features was raised in [KZBB]. An improvement of the stability was proposed that took the perimeter of MSER nodes into account. This improvement avoids favoring regions with rounder contours during MSER selection. A similar computation could be done here. Yet the simple extension proved to work well without this improvement.

We now adapt these definitions to the case of meshes.

.. Level set trees: an extension to D mesh surfaces

Let (V, T ) be a set of vertices and triangles sampled from a 2-manifold M embedded in the D Euclidean space R 3 . Points v ∈ V are linked to other points of V by edges forming triangles t ∈ T . We will assume that each edge is adjacent to either one or two triangles, so that each point belongs to at least one triangle. This means that there is no orphan edge and no orphan point, and that the mesh has no edge adjacent to three triangles. In other terms, the mesh is a manifold. To define a level set tree we need a real function defined on the mesh. The function H: V → R associating with each vertex v its mean curvature H(v) is chosen as our example throughout the chapter. Mesh regions will be defined as unions of mesh triangles. A level set tree requires a topology and therefore an adjacency relation on the mesh. Two triangles will be called adjacent if they share an edge. With this definition, analogous to -connectivity on D images, two regions sharing a vertex but no edge are disconnected. The main differences with the twodimensional case are that the mesh itself can be disconnected, and that it usually contains scanning holes. If the mesh is disconnected the component tree is a forest. The algorithm will process independently each tree. Section .. explains how to handle scanning holes.

As for images, connected components of upper level sets can contain topological holes which are themselves connected components of lower level sets (Fig.

.).

Monasse [MG] therefore proposed to build a shape tree, which is a fusion of the upper level set component tree with the lower level set component tree. Since an upper component tree is faster to build, and since it is the appropriate object to perform MSER extraction on the mesh, we limit ourselves to the upper component tree.

tel-00557589, version 1 -19 Jan 2011

.. Building the component tree for meshes

Set Λ = {h 1 < h 2 < • • • < h n } a set
of quantized values of H on V . The triangles t ∈ T will be ordered by setting Level(t) = max{k| min v∈t (H(v)) ≥ h k }. Referring to the D case, we can say that the triangles play the roles of pixels and that the real function used for building the component tree is Level(t). Algorithm  describes the construction of the tree of level set components on T .

Algorithm : Building the Component Tree Forest

Data: A list of the mesh triangles

F

Sketch of Proof .

For simplicity, in the algorithm and in the comments below, k denotes the level h k . Only the order of the levels matters, and it is reflected by k.

Algorithm  is strongly based on the triangle processing order from higher to lower levels. Indeed, while the algorithm is building nodes at level k, the only previously built nodes have a higher level. This fact is used for expanding the node: the expansion of a level k node is performed by successively adding the not previously added neighboring triangles of level k to the node triangle set. When expanding a new node N new , we may encounter triangles already processed (line ) and belonging to a node N at level l (then necessarily l > k; otherwise the node would not have yet been created). By going up in the parent-child relation, we can retrieve N a (line ), the last created ancestor of N (which can be N itself if it has no parent). For the same reason as before, N a 's level must be superior or equal to k. If N a has the same level as k, then it belongs to the same node, and one can merge both nodes (line ) by merging their triangle sets and adding N a 's children to the set of N new 's children (and then deleting N a ). If N a has a level larger than k, then N a is a child of N new (line ), and we can set N a 's father to be N new and add N a to the set of N new 's children.

During the construction of each node, track is kept of the triangles which are contained in this node, but not in its children. This way, each triangle belongs to a single node. The algorithm also keeps track of the area A (sum of the areas of triangles belonging to the node and to its descendants) while building the tree. This information is used in Mesh-MSER.

Algorithm  being similar to [NC], the complexity of building the forest is also quasi-linear. Indeed, it starts by sorting the triangles, an O(N log N ) step. When expanding a node, the computationally demanding case is when the encountered triangle belongs to an already created node. This requires finding the last created ancestor, which entails some traversing node operations, merge triangle lists (constant time) and add areas (1 operation). Since each triangle is processed only once, the total complexity is roughly N log N .

There are as many trees in the forest as nodes with no parent. By arguments similar to [MG] one can prove: Theorem . Any quantized scalar function H on a manifold mesh (V, T ) can be reconstructed from its component tree. This means that no information is lost on H in the component tree.

Assume the quantized levels of H are h k , k = 1, • • • , n. We want to extract local maximal elements (in the MSER sense) in the level set tree of H. Call δ > 0 tel-00557589, version 1 -19 Jan 2011 the level step used for computing the stability coefficient of a node N h k at level h k . The set of test levels h 1 < • • • < h n must therefore be complemented with all levels (h i + δ) and (h iδ). For each node N h i , going up in the tree hierarchy yields its descendant N h i +δ with level h i + δ, and going down yields its ascendant N h i -δ with level h iδ. Once these nodes are at hand, the stability coefficient is simply

q(N h i ) = |N h i -δ |-|N h i +δ | |N h i |
. For simplicity, we assumed in the previous relation that N h i only has one descendant with level h i + δ, which is not necessarily true. In the case of multiple descendants with level h i + δ, the sum of their areas is used instead of |N h i +δ |. Remark that N h i -δ is not necessarily the father of N h i , since there is no condition on the size of δ relatively to min i (h i+1h i ). For the same reason N h i +δ is not necessarily a son of N h i . This is why we must go up and down father-son relations in the component tree.

. Extracting maximally stable regions from the component tree

.. Mesh-MSER: the algorithm

Algorithm  describes how to extract maximally stable extremal regions (MSERs).

Starting from the component tree, this is an easy task. The tree structure gives a quick access to the area variations between N h i +δ ⊂ N h i ⊂ N h i -δ . When computing the stability coefficient, topological changes are authorized by the algorithm, whereas the original image-MSER technique only compares nodes on branches with no bifurcation. This way more lines are found than in the original D image method. The node merging procedure () is also standard. It generates a subtree whose nodes are exclusively the maximally stable regions. Merging a node N h i+1 into its father N h i requires a) removing N h i+1 from the set of N h i 's children, b) adding all of N h i+1 's children to N h i 's children, c) setting their father to be N h i , d) merging the list of triangles and e) updating the areas accordingly.

.. Triangle classification

Each selected node being given an index l, algorithm  yields a triangle classification L which, with any triangle t of the mesh, associates the label of the highest node containing the triangle (i.e., the label of the node with highest level containing the triangle). Because of the tree structure, it may occur that two triangles with label l are not connected by triangles with label l. Then the node must be split into different parts with different labels.

tel-00557589, version 1 -19 Jan 2011

Algorithm : Mesh-MSER Data: δ a step for computing stability coefficients and a set of levels

h 1 < • • • < h n Result:
A labeling of the triangles and the borders of the detected MSERs  Build the component tree with levels h i ± δ;  for each tree node N h i (where h i is the node level) do  Look for all descendants of N h i with level h i + δ and the sum of their areas A h i +δ ;  Look for the ascendant of N h i with level h iδ and its area

A h i -δ ;  q(N h i ) = A h i -δ -A h i +δ A(N h i )
;

 for each node N h i do  Get q h i+1 the minimal stability of the descendants of N h i with level h i+1 ;  Get q h i-1 the stability of the N h i 's ascendant with level h i-1 ;

 if q(N h i ) < q h i+1 and q(N h i ) < q h i-1 then  Select Node N h i ;
 for all non selected tree nodes N do

 

Merge the node N with its father;

 Associate to each triangle the index of the node with highest level containing it; tel-00557589, version 1 -19 Jan 2011 Figure .: Extracting the border of a region (the region border with no interpolation at all is shown in the blue, and the interpolated region border is the dotted green line). Notice that a part of both borders coincide with the mesh border

.. Border extraction from selected region

Extracting the borders of each selected region is a simple task using the available list of triangles for each node. From this list one can extract those which share an edge with a triangle of lesser level. This yields a set of border edges which can be linked. Since we are dealing with meshes built on raw point sets with a triangulation method which does not fill in the scanning holes, the component border line may encounter a scanning hole. To extract the line, we first extract the set B of edges belonging to at least one triangle of the connected component which is either a hole border or a component border. Starting from an edge of B which is not a hole border, a line is extended by linking edges from B. This way only closed contour lines are built, which are not hole borders, but can partially coincide with hole borders.

. Implementation and results

In the experiments herewith, the function H on the mesh will be the mean curvature. Choosing the levels h i is another question. Since curvatures are real numbers and are estimated only up to a given estimation error depending on the curvature estimation method, using all levels for the h i would not be a good solution. The adopted solution is to use equally spaced bins and δ equal to the quantization step.

The Mesh-MSER algorithm was implemented in C++. On a 1.5 GHz dual core laptop, without any particular effort on code optimization, the whole process lasts tel-00557589, version 1 -19 Jan 2011 for less than one minute for a 500000 vertices mesh. All shapes presented in this section (with the exception of the Stanford Urbis Romae pieces) were acquired using a triangulation laser scanner which yields a high precision dense point cloud (with some acquisition holes though, as can be noticed on fig. .). The nonoriented point cloud was first oriented, its curvature computed at all raw points, and an interpolating mesh built using the method described in [DMMSL].

.. Results on mechanical and geometrical shapes

The first experiment is a sanity check on simple a diamond shaped volumetric pattern (fig. .) with 150k vertices and 400k triangles. Mesh-MSER surrounds all geometrically relevant areas, namely the facets and a single connected region containing all vertices. It could be objected that a single threshold on the curvature would have sufficed to obtain the facets. But, even in that simple case, it was not obvious to predict the right curvature threshold. Furthermore, a simple curvature threshold would have delivered many small extremal regions due to noise inside the facets, which are actually fused to the facets by Mesh-MSER. The second industrial example is a mesh acquired from a water pump (2.5 million vertices, 4 million triangles), whose the mesh was again built directly on the raw data. This object has many acquisition holes (see fig. .). The final classification gives 200 regions. For better visualization random colors were given to the regions. The algorithm automatically separates edges from plane or curved parts. The segmentation of such a huge cloud into only  regions promises to enable a further model analysis, facet by facet.

.. Archeological pieces: comparative results

The next test (fig. .) was performed on a good quality mould of an archaeological piece, which was subjected to a massive scan followed by Mesh-MSER. Mesh tel-00557589, version 1 -19 Jan 2011 and curvatures were obtained using [DMMSL]. This small prehistorical figurine (. B.C.), "La Dame de Brassempouy" is only  centimeters tall. The mesh has approximately 300k vertices and 500k triangles. Notice how each detail of the hair dress is segmented out. This preliminary exploration of the capabilities of Mesh-MSER was continued on the Stanford Forma Urbis Romae database, containing hundreds of archaeological artifacts coming from a broken stone map of Roma (see [KTN * ]). A challenge of this project is to solve the jigsaw puzzle and rebuild the map. It is a crucial test for the Mesh-MSER method to check whether or not it extracts the engraved symbols and drawings figuring the town map, and whether it does it better than what can be done with D-MSER or with a Canny edge detector from simple photographs. Pictures of fragments g and u are given along with the result of MSER extraction on figs . and .. The experiment of fig. .(b) shows Mesh-MSER working on these engravings with a high performance, comparable to the best D image MSER performance on pictures containing high contrasted trademarks and logos [MCUP]. Indeed, all visible symbols and all features of the map plan are faithfully extracted, with very few outliers.

This experiment can be pushed further. Indeed, the pieces being rather flat, a direct comparison of D-and Mesh-MSER on their main facet makes sense. The Mesh-MSER result compares advantageously to D level line or edge extraction methods applied to a picture of the same object (fig. . (a)). The comparison shows that it is far more reliable to detect boundaries on the D mesh. Finally, several strategies for extracting curvature level lines on meshes are compared on fig. . with fragment g of the Stanford Urbis Romae database. Here again the Mesh-MSER results seem to be complete, accurate, and without outliers. The experiment compares the choice of curvature level lines made by MSER with a simple threshold based on level line length. Although this choice indeed removes noisy level lines, it also loses many meaningful ones, and adds anyway a spurious threshold parameter.

To compare the obtained results with those presented in [ZTS], the same data set point (a fragment of antique vase) was used. The mentioned reference is very similar in scopes to MSER: it proposes a two scale analysis on a mesh by defining a "base" and "height function". The lines shown in [ZTS] are level lines of the height function. The base is implicitly defined by its gradient, by a sophisticated variational procedure. Here we used a similar height function to get a relevant comparison. The height function is defined as the difference between the surface and its smoothed out version by a large scale mean curvature motion. The Mesh MSER method can then be directly applied (fig. .) on this scalar function.

. Using the same approach with other functions

Other functions can be used for extracting valuable informations on surfaces. Let us for example consider the case of a digital elevation model. There are two ways of considering the data: either as an image and process it as gray-valued image using D-MSER or as a set of grid points with height values giving a set of D coordinates with a mesh structure. Then the process takes into account a much finer information since it adds precise triangle areas (which depend on the height) instead of adding a constant 1 area to compute node areas.

. Conclusion

This chapter introduced a fast level set tree method, Mesh-MSER, applicable to any scalar function defined on a mesh. This method is a direct extension of classic D image analysis tools building trees of level sets components or of level lines. Using the fact that the curvature is the most straightforward scalar function defined from and on a mesh, the method was used to segment meshes into maximally stable extremal regions (MSERs) of the curvature. Future work will focus on the exploitation of this structure. Indeed the experiments clearly point out the possibility of using the detected curves and regions to perform pattern recognition of complex objects such as the Urbis Romae fragments. On the other hand the method provides automatic segmentations of industrial objects into edge parts and parts with constant or slowly varying curvature, for which spline or conical models should easily be estimated.

Introduction

The output of laser scanners or any surface acquisition system is usually a set of points sampled with more or less precision on the surface. In some cases the result comes as a mesh, i.e., as a set of triangles linking points. In other cases, no such information is given and the machine data is simply an unorganized point cloud.

In this chapter we focus on the mathematical processing of the surfaces defined as point clouds. We also explain in terms of differential operators what happens in some of the most common surface regularization processes. The field of numerical surface analysis has been widely studied over the fifteen past years, due to the development of computer graphics. Yet, in most cases the starting surface representation is a mesh, for the simple reason that meshes are much easier to handle than a point cloud. Indeed, numerically processing a surface always involves finding the neighbors of a surface sample, and a mesh is oriented and has a surface topology. The neighbors of a vertex are all points linked by a specified number of edges to the center vertex. Yet in case of highly irregular meshes the tel-00557589, version 1 -19 Jan 2011 processing can be problematic because then, triangle areas vary a lot. The most common mesh reconstruction methods from a raw point cloud define a signed function over R 3 representing the distance to the object, and then extract the 0 level set which approximates the object surface. See (e.g.) [CL], [HDD * ], [CBC * ], [Kaz], [ACSTD], and the most popular such method [KBH], which solves a Poisson equation to build the indicator field. These methods vary in the approach to compute the distance function, but all extract the 0-level set by using the marching cubes algorithm ( [LC], [KBSS]). In such a meshing process, the initial raw points are irremediably lost, and this incurs into a loss of resolution. This explains why it may be relevant to process directly an unstructured raw point cloud. If no mesh is given as input data, then looking for neighbors means looking for points within a given distance of the query point. This neighborhood definition raises other problems: if the neighborhood is too large then we risk including two different parts of the surface in the neighborhood, and in the normal case of an irregular sampling empty neighborhoods can happen. Therefore the radius choice is crucial in the numerical implementation of any local operator.

The reminder of this chapter is divided as follows: section . reviews surface operators and motions previously defined on meshes and on point clouds, section . gives the necessary definitions and tools. Section . analyzes regression free curvature estimates. Section . and . analyze and compare the motions given by the projection on simple regression surfaces: a plane and a degree 2 polynomial surface. Finally, section . shows comparative numerical experiments.

. Previous work

.. Curvature estimation and surface motions defined on meshes

Reviews for curvature estimation on surfaces can be found in [Rus], [MSR] or [MSR]. Curvature tensor estimation methods were pioneered by Taubin [Taua] who presented a simple approximation for computing the directional curvature in any tangent direction. Then the curvature is computed in all incident edge directions and a covariance matrix of the edge direction weighted by their directional curvatures and the area of the two incident triangles is built. Eigenvectors and eigenvalues of this covariance matrix yield a simple expression of the principal curvatures and curvature directions.

Other curvature tensor estimation methods include [Rus] where the tensor is estimated by building a linear system over the tensor coefficients. This linear system expresses the constraints that multiplying the tensor by an edge direction should give the difference of the edge's endpoints normals. The same method is tel-00557589, version 1 -19 Jan 2011 applied to find the curvature derivatives. Normals are also used in [TRZS] to give a piecewise linear curvature estimation (see also [TT])

To avoid computations of derivatives, a new kind of curvature estimation method has been proposed, based on domain integration. In [YLHP], [PWY * ] and [PWHY], the intersection of the surface with either a sphere or a ball centered at a vertex is analyzed: the covariance matrix of this domain is computed and eigenvalues are expressed in terms of principal curvatures. By increasing the neighborhood radius, the curvature estimate can be made multiscale. A very interesting feature of these methods is that they do not rely on high order derivatives and are therefore more stable. Other methods to compute the curvature include the normal cycle theory [CSM].

Surface motions have also been studied as a part of a mesh fairing process. A key method was introduced by Taubin in [Taub] who considered a discrete Laplacian for a mesh V with vertices

v i ∂V ∂t = λL(V ), L being a discretization of the Laplacian L(v i ) = 1 card N (v i ) j ∈ N (v i )(v j -v i )
where N (v i ) is the set of vertices linked by an edge to v i (1-ring neighborhood). This formulation is widely used. Indeed [DMSB] uses a similar "umbrella" operator. [GH] also computes the discrete Laplacian for all mesh vertices and the eigenvector and eigenvalues of the Laplacian are computed. By removing the smallest eigenvalues, a fair mesh is obtained.

A well known formulation of the Laplace Beltrami operator is the famous cotangent formula ( [MDSB]). This formula writes:

∆v i = 1 2 j∈N v i (cotanα ij + cotanβ ij )
where v i is a vertex of the mesh, N 1 (v i ) its one ring neighborhood, α ij and β ij are the angle opposite to edge v i v j in the two triangles adjacent to v i v j . This has been used to compute the surface intrinsic equation.

.. Curvature estimation and surface motions defined on point clouds

We now discuss the rare approaches dealing directly with point clouds.

In [UH], a scale space decomposition method for point clouds is introduced. At a point p, the intersection of the surface with a plane containing the normal to p is a curve. The curvature of these curves can be easily computed. An operator is defined that moves each point in the normal direction by a factor equal to this curve curvature. Averaging over all possible planes yields a mean curvature motion. The non-uniform sampling problem is then solved by using a density normalized kernel which removes the dependence of the result on variations in sampling density. The scale space is then used to select "scale-space extrema". At each scale the point tel-00557589, version 1 -19 Jan 2011 motion is considered. Introducing a scalar function on the displacement norm, the authors claim that they recover the characteristic scales of the surface (the introduced function is extreme for the characteristic scales).

More theoretical work which is not tested on real raw surfaces can be found in [MOG] where a method to compute curvatures and normals based on voronoi covariance matrix is introduced.

In [Tan], the proposed framework for curvature estimation at a particular point is based on a set of curves representing the local neighborhood of the point under consideration. By considering the neighbors of p, one can have the set of triplets (p i , p j , p) and those triplets can be used to define parametric space curves by quadratic polynomial interpolation with p(0) = p i , p(1) = p j and p(t) = p where t = |p-p i | |p-p i |+|p j -p| . This allows to approximate maximum and minimum curvature values as the minimum and maximum normal curvature values for all possible point triplets. This method can be used either on meshes or point clouds.

In [KSNS], the authors propose a statistical estimation of curvature of point sampled surfaces based on IRLS (iteratively reweighted least squares) and Mestimators. Position difference vector ∆p and normal difference vector ∆n are used to define a linear system yielding a first estimate of the curvature tensor. Then residuals are computed and used to weight the samples and the objective function is minimized by iterative reweighting of point samples. This yields the final curvature tensor estimate. This procedure is rather time-consuming.

Finally in [BSW], an algorithm to compute the Laplacian of a function defined on point clouds in R d is proposed along with convergence proofs. Yet the model is not tested on real surfaces.

Neighborhood covariances being used already for normal estimation, the idea to express fundamental forms as covariances matrices was introduced. Next section reviews the different covariance techniques considered in the literature.

.. Curvature estimation using covariance techniques

There are few covariance approaches and none of them has been analyzed mathematically yet. Nevertheless, covariance methods can be an elegant alternative to surface regression. This chapter performs this analysis linking the various covariances to surface curvatures. Three papers use covariance matrices for curvature estimation: [BC], [LT] and [PGK].

[BC] considers the neighbors (p i ) of a point p. The second fundamental form is then defined as the covariance matrix of vectors pp i projected onto the tangent plane of the surface at p. An equivalent of the Gauss map is also introduced: it is the covariance matrix of the neighbors unit normals projected onto the surface tangent plane at point p. The eigenvectors are said to give the principle directions. In fact these two covariance matrices are inspired from [LT]. Indeed, [LT] first tel-00557589, version 1 -19 Jan 2011 proposed to compute the covariance matrix of the normals of point p's neighbors, to extract the principal eigenvalues which correspond to the principal curvatures of the surface at p. The last covariance method, introduced in [PGK], is not claimed to be explicitly linked to surface curvatures or fundamental forms, yet it is used to account for the surface geometric variations. Consider the covariance matrix of vectors p i where is the barycenter of the neighborhood of p. The surface variation is then defined as the ratio of the least eigenvalue over the sum of all eigenvalues of this covariance matrix. This quantity has the nice property that it is bounded between 0 (flat case) and 1/3 (isotropic case).

All these methods will be described in more details and analyzed in section .. Another interesting technique for surfaces defined by point clouds is the Moving Least Squares Surfaces. [LS] first introduced MLS surfaces as follows. Given a data set of points {p i } i (possibly acquired by a D scanning device), a smooth surface M based on the input points is defined. The goal is to replace the points p defining M with a reduced set R = {r i } defining an MLS M ′ surface which approximates M. The surface given by the point cloud is expected to be a -manifold, C ∞ smooth. The authors also define a bounding error ε such that d(M, M ′ ) < ε, where d is the Hausdorff distance.

.. Moving Least Squares Surfaces

The projection of a point on MLS surface is defined as follows: given a point p, find a local reference domain (plane) for p. The local plane H is computed so as to minimize a local weighted sum of square distances of the points p i to the plane. (Thus is is a weighted regression plane). The weights attached to p i are defined as functions of the distance from p i to the projection of p on plane H, rather than the distance from p i to p.

Assume q is the projection of p onto H, then H is found by locally minimizing

N i=1 (< n, p i > -D) 2 θ( p i -q ),
where θ is a smooth, monotone decreasing function, which is positive on the whole space. We can set q = p + tn for some t ∈ R, which gives

N i=1 (< n, p i -p -tn >) 2 θ( p i -p -tn )
The local reference domain is then given by an orthonormal coordinate system on H so that q is the origin of this system. The reference domain for p is used to tel-00557589, version 1 -19 Jan 2011 compute a local bivariate polynomial approximation to the surface in a neighborhood of p. Let q i be the projection of p i onto H, and f i =< n, p iq i >. In this local coordinate system, let (x i , y i ) be the coordinates of q i on H. The coefficients of the polynomial are computed by minimizing the least square error:

N i=1 (g(x i , y i ) -f i ) 2 θ( p i -q )
The projection of p onto M is defined by the polynomial value at the origin, i.e., q + g(0, 0)n = p + (t + g(0, 0))n. Thus, given a point p and its neighborhood, its projection onto the MLS surface can be computed.

The approximation power of MLS surfaces was shown in [Lev] and the first applications were introduced in [ABCO * ], [AK] and [Lev].

MLS surfaces turned out to be not only theoretically powerful but provided fine implementations for either rendering or up-sampling and down-sampling point sets ([ABCO * ], [PGK]). Variants of MLS were proposed mostly for a better preservation of sharp edges in surfaces defined by point clouds ( [OGG], [LCOL], [FCOS]).

The same framework to build a scale space for point clouds in [PKG]. The surface is evolved through a diffusion process ∂p ∂tλ • ∆p = 0, where p is a point of the surface, λ a diffusion parameter and ∆p = Hn is the Laplace Beltrami Operator (H is the curvature and n the normal at point p, this is the decomposition process). By remembering the set of displacements D i (p) of each point p we have a reconstruction operator. The choice of the Laplacian discretization is very important: a first possibility is to use the standard mesh Laplacian techniques ( [Taub]) adapted for point clouds using the k-nearest neighbors instead of the one ring neighborhood. Another possibility is to use the weighted least squares projection ( [HG], [KG]): the surface is iteratively projected onto the plane defined by the weighted barycenter o and the normal estimated using the weighted neighborhood covariance matrix. Weights are a simple gaussian ponderation on the distance to p and the size of the gaussian kernel is a parameter that controls the amount of smoothing. This projection process is in fact an order 1 projection motion (MLS) that is analyzed in this chapter. To make the projection more efficient, [PKG] proposed to subsample the point cloud. It yields a scale space decomposition where at each level the surface is smoothed and sub-sampled. The scale space decomposition is then applied to the multi-scale freeform deformation and to the morphing problem, with satisfactory results.

The moving least squares (MLS) were used to estimate curvatures. For example, in [YQ], the authors use the MLS framework to build a closed form solution for curvature estimation. Indeed, surfaces implied by point clouds can be seen as the zero level set of an implicit function f . The gradient and Hessian Matrix of f is built. Finally, using formulas for the Gaussian and the mean curvature depending on the Hessian and gradient of f , those curvatures are computed.

In [CP], the problem of estimating differential quantities on point clouds is tel-00557589, version 1 -19 Jan 2011

recast to that of fitting the local representation of the manifold by a jet. A jet is simply a truncated Taylor expansion. A n jet is a Taylor expansion truncated at order n.

A jet of order n contains differential information up to the n-th order. In particular it is stated that a polynomial fitting of degree n estimates any k th order differential quantity to accuracy O(h n-k+1 ). This implies that the coefficient of the first fundamental form and unit vector normal are estimated with O(h n ) precision and the coefficients of the second fundamental form and shape operator are approximated with accuracy O(h n-1 ), and so are the principal directions. In order to characterize curvature properties, the method resorts to the Weingarten map A of the surface, also called the shape operator, that is the tangent map of the Gauss map. Recall that the first and second fundamental forms I, II and A satisfy II(t, t) = I(A(t), t) for any vector t of the tangent space. Second order derivatives are computed by building the Weingarten map of the osculating jet whose eigenvalues are the principal curvatures. Note that the described methods can be used either with a mesh or with a point cloud. Jets are in fact very related to MLS surfaces. Indeed, to estimate differential quantities a polynomial fitting of degree n is done, which is exactly what MLS does. Therefore the analyzis given in section . giving the equation governing MLS and MLS motions are valid for the jets too.

In this chapter, we give the exact partial differential equation that governs the MLS projection motion for order 1 and 2 MLS surfaces. We then link these PDEs to the surface curvatures.

Next section provides some tools to analyze numerically point cloud motions. The following analysis is in spirit very close to the image filter analysis performed in [BCM].

. Tools for numerical analysis of point cloud surface motions 

We always assume the existence of a smooth surface M supporting the point set. These surfaces are the boundaries of solid objects and can therefore be assumed to be locally Lipschitz graphs. However, for a mathematical analysis of smoothing algorithms and curvature estimations on the surface, we shall always assume that the surface is a C ∞ embedded manifold, known from its samples denoted by M S . This is not a limitation, in the sense that any finite sample set can be anyway interpolated by an arbitrarily smooth surface.

z = f (x, y) = 1 2 (k 1 x 2 + k 2 y 2 ) + o(x 2 + y 2 ).  (.)
Notice that the sign of the pair (k 1 , k 2 ) depends on the arbitrary surface orientation.

Points where k 1 and k 2 have the same sign are called parabolic, and points where they have opposite signs are hyperbolic. Consider two kinds of neighborhoods in M for p defined in the local intrinsic coordinate system (p, t 1 , t 2 , n):

• the spherical neighborhood B r = B r (p) ∩ M is the set of all points m of M with coordinates (x, y, z) satisfying (xx p ) 2 + (yy p ) 2 + (zz p ) 2 < r 2

• the cylindrical neighborhood C r = C r (p)∩M is the set of all points m(x, y, z) on M such that (xx p ) 2 + (yy p ) 2 < r 2 .

The spherical neighborhood in the sampled surface M 2 is the only neighborhood to which there is a direct numerical access, and the discrete operators is defined with it. Nevertheless, for the forthcoming asymptotic numerical analysis, the cylindrical neighborhood will prove much handier than the spherical one. The next technical lemma justifies its use in theoretical calculations.

= (1 + O(r 2 ))
x 2 +y 2 <r 2 , (x,y,z)∈M f (x, y)dxdy.

This last form is amenable to analytic computations. Consider polar coordinates (ρ, θ) such that x = ρ cos θ and y = ρ sin θ with -r ≤ ρ ≤ r and 0 ≤ θ ≤ π. Then for m(x, y, z) belonging to the surface M, we have z = 1 2 ρ 2 (k 1 cos 2 θ + k 2 sin 2 θ) + O(r 3 ). Fixing θ we obtain z = 1 2 ρ 2 k(θ) + O(r 3 ), where k(θ) = k 1 cos 2 θ + k 2 sin 2 θ. The condition that (x, y, z) belongs to the neighborhood B r can therefore be rewritten as ρ 2 + z 2 < r 2 , that is

ρ 2 + 1 4 k(θ) 2 ρ 4 < r 2 + O(r 5 )
Computing the boundaries ±ρ(θ) of this neighborhood yields ρ(θ

) 2 + 1 4 k(θ) 2 ρ(θ) 4 -r 2 + O(r 5 ) = 0. Thus ρ(θ) 2 = -1 + 1 + k(θ) 2 (r 2 + O(r 5 )) 1 2 k(θ) 2 .
This yields ρ(θ) = r -1 2 k(θ) 2 r 3 + O(r 4 ). We shall use this estimate for the error term E appearing in This lemma will prove very useful for the rest of the chapter and in particular in the next section where analysis are given for various curvature estimates.

Lemma , as well as all theorems in the remainder of this chapter will assume that the surface is a uniform Lebesgue measure. A constant sampling density is therefore necessary. This constant density is approximated on discrete data by weighting each point by a weight inversely proportional to its initial density. More precisely, let p be a point and N r (p) = M s ∩ B r (p). Each point q should ideally have a weight 0 ≤ w(q) ≤ 1 such that q∈Nr(p) w(q) = 1. This amounts to solve a huge linear system. For this reason, we shall be contented with ensuring q∈Nr(p) w(q) ≃ 1 by taking w(p) = 1 ♯(Bp(r)) , as proposed in [UH].

. Regression-free curvature estimates

This section finds the form of the differential operators underlying four different discrete schemes based on local cloud point statistics, and proposing discrete analogues of the "second fundamental forms" or of the "principal curvatures". These discrete schemes have very simple and robust form, being based on the computation of local moments and eigenvalues of the point cloud. We shall see that they actually compute nonlinear differential operators linked to the principal curvatures.

.. A discrete "second fundamental form" [BC]

Let(p i ) i∈1•••N be the set of neighbors of a point p with normal n. This paper proposes to build the "second fundamental form matrix" as follows. (Although this covariance matrix is not, as we shall see, consistent with the second fundamental form, it is thus called in this paper, and actually has, as we shall see, the principal directions as eigenvectors.) Let s i = (p ip) T • n, let t 1 , t 2 be two orthonormal vectors of p's tangent plane, and

α i = s i • (p i -p) • t 1 (p i -p) • t 2 = ((p i -p) T • n) • (p i -p) • t 1 (p i -p) • t 2 .
tel-00557589, version 1 -19 Jan 2011

The α i are projections of the vectors (p ip) onto the tangent plane to p, weighted by their distance to this plane. The "second fundamental form matrix" is the covariance of these vectors, namely

Σ d = N i=1 (α i -α m ) • (α i -α m ) T (.)
where α m = 1 N N i=1 α i and in Σ d the d stands for "discrete". To compute the underlying differential operators, two assumptions are made throughout this chapter. The first one is that the surface sampling is uniform with respect to the area measure on the surface. The second one is that this sampling is dense enough, so that the averages taken on neighborhoods can be interpreted as integrals. Under this interpretation, we can reinterpret the sum in (.) as an integral on a cylindrical neighborhood of p, assuming the data point set to be a locally smooth manifold. In the local intrinsic surface coordinate system at point p, (p, t 1 , t 2 , n), the surface can be written as a graph z = 1 2 (k 1 x 2 + k 2 y 2 ) + o(r 2 ). Thus the vectors α i are replaced by a continuous vector α(x, y) defined by

α(x, y) = 1 2 (k 1 x 2 + k 2 y 2 ) • x y = 1 2 • k 1 x 3 + k 2 y 2 x k 1 x 2 y + k 2 y 3 + o(r 3 ). (.)
Under the interpretation taken above the "second fundamental matrix" rewrites

Σ = Br (α(x, y) -α m ) • (α(x, y) -α m ) T dm(x, y) (.)
where α m = 1 meas(B r ) Br α(x, y)dm(x, y).

(.)

The proposition made in [BC] is to extract the surface principal curvatures and their corresponding directions at p from this covariance matrix, as its eigenvalues and eigenvectors. The next theorem checks if this works asymptotically in the continuous model.

Theorem . The eigenvectors of the "second fundamental form matrix" Σ give the principal directions with error o(r 8 ). But the eigenvalues of Σ are not the principal curvatures as they satisfy

λ 1 = 1 4 r 8 8 π 8 (5k 2 1 + 2k 1 k 2 + k 2 2 ) + o(r 8 ) and λ 2 = 1 4 r 8 8 π 8 (k 2 1 + 2k 1 k 2 + 5k 2 2 ) + o(r 8 )
where k 1 and k 2 are the principal curvatures at p tel-00557589, version 1 -19 Jan 2011

Proof. In the continuous model α m therefore is close to zero because the integrated function is odd on a symmetric domain. More precisely, using Lemma  in (.) and writing α m = (α mx , α my ),

α mx = 1 2πr 2 x 2 +y 2 <r 2 (k 1 x 3 + k 2 y 2 x i + o(r 5 ))dxdy = o(r 3 )
and similarly

α my = o(r 3 ).
By Lemma  again, the covariance matrix (.) satisfies Σ = x 2 +y 2 <r 2 α(x, y) • α(x, y) T dxdy + o(r 8 ), and, using (.), we can calculate its four components as follows.

Σ 11 = 1 4 2 +y 2 <r 2 (k 1 x 3 + k 2 y 2 x) 2 dxdy + o(r 8 ) = 1 4 x 2 +y 2 <r 2 (k 2 1 x 6 + 2k 1 k 2 x 4 y 2 + k 2 2 y 4 x 2 )dxdy + o(r 8 ) = 1 4 θ∈[0,2π],ρ∈[0,r] ρ 6 (k 2 1 cos(θ) 6 + 2k 1 k 2 cos(θ) 4 sin(θ) 2 + k 2 2 cos(θ) 2 sin(θ) 4 )ρdρdθ + o(r 8 ) = 1 4 r 8 8 θ∈[0,2π] (k 2 1 cos(θ) 6 + 2k 1 k 2 cos(θ) 4 sin(θ) 2 + k 2 2 cos(θ) 2 sin(θ) 4 )dθ + o(r 8 ) = 1 4 r 8 8 (k 2 1 5π 8 + 2k 1 k 2 π 8 + k 2 2 π 8 ) + o(r 8 ) = 1 4 r 8 8 π 8 (5k 2 1 + 2k 1 k 2 + k 2 2 ) + o(r 8 )
By exchanging the roles of k 1 , k 2 , and x, y respectively, we get

Σ 22 = 1 4 r 8 8 π 8 (k 2 1 + 2k 1 k 2 + 5k 2 2 ) + o(r 8 ).
Σ being a symmetric matrix, Σ 12 = Σ 21 and the integrated function being odd,

Σ 12 = 1 4 x 2 +y 2 <r 2 (k 1 x 3 + k 2 y 2 x)(k 1 x 2 y + k 2 y 3 )dxdy + o(r 8 ) = 1 4 x 2 +y 2 <r 2 (k 2 1 x 5 y + 2k 1 k 2 x 3 y 3 + k 2 2 y 5 x)dxdy + o(r 8 ) = o(r 8 ).

tel-00557589, version 1 -19 Jan 2011

Thus, Σ is equivalent to a diagonal matrix whose principal directions are t 1 and t 2 , which validates the theoretical requirements, t 1 and t 2 being the principal directions at point p. However, the corresponding eigenvalues are

λ 1 = 1 4 r 8 8 π 8 (5k 2 1 + 2k 1 k 2 + k 2 2 ) + o(r 8 ) λ 2 = 1 4 r 8 8 π 8 (k 2 1 + 2k 1 k 2 + 5k 2 2 ) + o(r 8
) which are definitely different from λ 1 = k 1 and λ 2 = k 2 . Only the absolute values of k 1 and k 2 can actually be deduced from Σ.

.. Another discrete "second fundamental form"

Another method is also introduced in [BC] which, in a nutshell, computes the covariance matrix of the unit normal vectors projections onto the local tangent plane. By applying again the continuous asymptotic analysis of section .., we shall see in Theorem  that this method actually computes discrete approximations of the squares of the principal curvatures. The discrete algorithm is as follows. Let M be a C 2 surface and p be a point of M. Let (p i ) i be the neighbors of p in a ball neighborhood of radius r. Denote by n i the normal at p i and define v i as the projection of n i onto the tangent plane at p, then the computed "curvatures" are defined as the eigenvalues of the covariance matrix of the vectors v i . The vector v i being the projection of n i onto the tangent plane, we have:

v i = n i • t 1 n i • t 2 Set v m = 1 N N i=1 v i . Then this new discrete covariance matrix writes Σ d = N i=1 (v i -v m ) • (v i -v m ) T .
In the continuous framework, the local points on the surface have coordinates m(x, y) = (x, y, 1 2 (k 1 x 2 + k 2 y 2 ) + o(r 2 )) and the normal vector to this surface is ∂m ∂x (x, y)

∧ ∂m ∂y (x, y) = (-k 1 x, -k 2 y, 1) + o(r). It follows that v(x, y) = 1 1 + k 2 1 x 2 + k 2 2 y 2 -k 1 x -k 2 y + o(r), v m = 1 meas(B r ) v(x, y)dm(x, y),
and the continuous covariance matrix is

Σ := Br (v(x, y) -v m ) • (v(x, y) -v m ) T .
tel-00557589, version 1 -19 Jan 2011

Theorem . The eigenvalues of the covariance matrix Σ of the vectors v(x, y) in the spherical neighorhood B r are

k 2 1 r 4 π 4 + o(r 4 ) and k 2 2 r 4 π 4 + o(r 4 ).
Proof. Let us compute the mean v m of v(x, y) on the spherical neighborhood. By Lemma  the integral on a spherical neighborhood is asymptotically equivalent to the integral on a cylindrical neighborhood and more precisely,

(πr 2 )v m • t 1 = x 2 +y 2 <r 2 -k 1 x 1 + k 2 1 x 2 + k 2 2 y 2 dxdy + o(r 3 ) = x 2 +y 2 <r 2 -k 1 x(1 - 1 2 k 2 1 x 2 - 1 2 k 2 2 y 2 )dxdy + o(r 3 ) = - θ∈[0,2π],ρ∈[0,r] k 1 ρ cos(θ)(1 - ρ 2 2 (k 2 1 cos 2 θ + k 2 2 sin 2 θ))ρdρdθ + o(r 3 ) = - θ∈[0,2π],ρ∈[0,r] k 1 ρ 2 cos(θ) + Aρ 4 dρdθ + o(r 3 ) where A is a bounded constant = - k 1 3 θ∈[0,2π] r 3 cos(θ)dθ + o(r 3 ) = o(r 3 )
and similarly

(πr 2 )v m • t 2 = x 2 +y 2 <r 2 -k 2 y 1 + k 2 1 x 2 + k 2 2 y 2 dxdy + o(r 4 ) = o(r 3 )
Thus the coefficients of Σ satisfy, again by Lemma ,

Σ 11 = x 2 +y 2 <r 2 k 2 1 x 2 1 + k 2 1 x 2 + k 2 2 y 2 dxdy + o(r 4 ) = x 2 +y 2 <r 2 k 2 1 x 2 (1 -k 2 1 x 2 -k 2 2 y 2 )dxdy + o(r 4 ) = k 2 1 θ∈[0,2π],ρ∈[0,r] ρ 3 cos 2 θ(1 -ρ 2 (k 2 1 cos 2 θ + k 2 2 sin 2 θ))dρdθ + o(r 4 ) = k 2 1 r 4 4 θ∈[0,2π] cos 2 θdθ + r 6 A + o(r 4 ) with A bounded = k 2 1 r 4 π 4 + o(r 4 )
tel-00557589, version 1 -19 Jan 2011

 Chapter . A numerical analysis of raw point cloud smoothing

Similarly,

Σ 22 = x 2 +y 2 <r 2 k 2 2 y 2 1 + k 2 1 x 2 + k 2 2 y 2 dxdy + o(r 4 ) = k 2 2 r 4 π 4 + o(r 4 )
and

Σ 12 = Σ 21 = x 2 +y 2 <r 2 k 1 k 2 xy 1 + k 2 1 x 2 + k 2 2 y 2 dxdy + o(r 4 ) = x 2 +y 2 <r 2 k 1 k 2 xy(1 -k 2 1 x 2 -k 2 2 y 2 )dxdy + o(r 4 ) = k 1 k 2 θ∈[0,2π],ρ∈[0,r] ρ 3 sin(θ) cos(θ) -ρ 5 sin(θ) cos(θ)(k 2 1 cos 2 θ + k 2 2 sin 2 θ)dρdθ + o(r 4 ) = k 1 k 2 4 θ∈[0,2π] sin(θ) cos(θ)r 4 dθ + Ar 6 + o(r 4 ) where A is bounded = r 4 k 1 k 2 4 θ∈[0,2π] sin(θ) cos(θ)dθ + o(r 4 ) = o(r 4 )
Thus at order 4, Σ can be considered diagonal and its eigenvalues are Σ 11 and Σ 22 in the principal directions t 1 and t 2 . They asymptotically give an approximation of each of the squared principal curvatures, but not their sign.

.. A third discrete "fundamental form"

The methods analyzed in sections .. and .. are akin to the original method introduced in [LT]. Indeed, in [LT] it was proposed to compute the covariance matrix of the normal vectors of the neighborhood (without projecting them in the local regression plane) and therefore get a 3 × 3 matrix instead of a 2 × 2 matrix. This is actually the simplest imaginable method and we shall see that it gives a result similar to section ... Theorem . Let M be a C 2 surface, let p be a point of M. Then the three eigenvalues of the covariance matrix C of the unit normals in a neighborhood of radius r around p are asymptotically respectively equal to  and to the squares of the principal curvatures at p.

Proof. A normal vector writes

tel-00557589, version 1 -19 Jan 2011

N = 1 1 + k 2 1 x 2 + k 2 2 y 2   -k 1 x -k 2 y 1   + o(r).
As in the previous sections, we easily obtain by Lemma , N mx = o(r), N my = o(r), N mz = 1 + o(r). Thus again by Lemma ,

C= x 2 +y 2 ≤r 2 1 1 + k 2 1 x 2 + k 2 2 y 2   k 2 1 x 2 k 1 k 2 xy k 1 x(1 -N mz ) k 1 k 2 xy k 2 2 y 2 k 2 y(1 -N mz ) k 1 x(1 -N mz ) k 2 y(1 -N mz ) (1 -N mz ) 2   dxdy+o(r 4 )
and, by calculations exactly analogous to Section ..,

C 11 = k 2 1 r 4 π 4 + o(r 4 ), C 22 = k 2 2 r 4 π 4 + o(r 4 ), C 12 = C 21 = k 1 k 2 x,y xydxdy = o(r 4 ), C 13 = C 31 = C 23 = C 32 = C 33 = o(r 4 ).
Thus the eigenvalues are asymptotically equal to k 2 1 r 4 π 4 and k 2 2 r 4 π 4 , which also gives back the squares of the principal curvatures of the surface, but not their sign.

.. A fourth discrete fundamental form: the surface variation

We shall now analyze a last variant introduced in [PGK], the so called surface variation. It is again based on a local covariance analysis. Unlike the previous methods, the surface variation was not claimed to be a curvature estimate, but to be a measure of the neighborhood shape. This subsection establishes again a link between this discrete quantity and the principal curvatures of the surface.

Let p be a point with given neighborhood B r . Let o be the barycenter of the neighborhood. In R 3 , the coordinates are written with superscripts e.g. the coordinates of a point u are (u 1 , u 2 , u 3 ). Thus, for i = 1, 2, 3,

o i = 1 card Br p k ∈Br p i k . The centered covariance matrix Σ = (m ij ) i,j=1,••• ,3 is defined as m ij = p k ∈Br (p i k - o i ) • (p j k -o j ) for i, j = 1, 2, 3. Let λ 0 ≤ λ 1 ≤ λ 2 be the eigenvalues of Σ with corresponding eigenvectors v 0 , v 1 , v 2 . For k = 0, • • • , 2, λ k = p i ∈Br (p i -o), v k 2 . (.)
Each eigenvalue gives the variance of the point set in the direction of the corresponding eigenvector. Since v 1 and v 2 are the vectors that capture most variations, they define the PCA regression plane. The normal v 0 to this plane is the direction v minimizing p i ∈Br (p io), v 2 . [PGK] defines the surface variation by

σ = λ 0 λ 0 + λ 1 + λ 2 . (.)

tel-00557589, version 1 -19 Jan 2011

This quantity measures the ratio of variance along the normal to the total variance. If the neighborhood is highly curved, its surface variation will be high and if the neighborhood is flat the surface variation will be small. This quantity has the property to be bounded between 0 (flat case) and 1/3 (isotropic distribution case).

Lemma . In the local intrinsic coordinate system, the barycenter of a neighborhood B r of point p has coordinates

x o = o(r 2 ), y o = o(r 2 ) and z o = Hr 2 4 + o(r 2 ), where H = k 1 +k 2 2
is the mean curvature at p. Proof. By Lemma  applied to the numerator and denominator of the following fraction, we have

z o = Br zdm Br dm = x 2 +y 2 <r 2 z(x, y)dxdy + O(r 5 ) x 2 +y 2 <r 2 dxdy + O(r 3 ) = x 2 +y 2 <r 2 1 2 (k 1 x 2 + k 2 y 2 ) + o(x 2 + y 2 ) dxdy x 2 +y 2 <r 2 dxdy + O(r 3 ) = 1 2πr 2 r ρ=0 2π θ=0 ρ 2 (k 1 cos 2 θ + k 2 sin 2 θ)ρdρdθ + o(r 2 ) = r 2 8π (k 1 π + k 2 π) + o(r 2 ) = Hr 2 4 + o(r 2 ).
A similar but simpler computation yields the estimates of x o and y o .

Theorem . In the local coordinate system the surface variation σ satisfies

σ = r 2 16 k 2 1 + k 2 2 2 - 1 3 k 1 k 2 + o(r 2 ) (.)
Proof. We need to explain what the covariance eigenvalues stand for. Each eigenvector v i and associated eigenvalue λ i represent a principal direction and the variation along this principal direction:

λ i = m∈Br om, v i 2 dm
Since we have λ 0 ≤ λ 1 ≤ λ 2 , we can see that λ 0 is associated to the direction with the least variation namely the normal direction to the surface oz. Since the eigenvectors form an orthonormal basis, we have

λ 0 + λ 1 + λ 2 = m∈Cr om, v 0 2 + om, v 1 2 + om, v 2 2 dm = m∈Cr om 2 dm

tel-00557589, version 1 -19 Jan 2011

This yields

λ 0 + λ 1 + λ 2 =
x 2 +y 2 <r 2

x 2 + y 2 + (zz o ) 2 dxdy and

λ 0 + λ 1 + λ 2 = πr 4 2 + λ 0 + o(r 6 ) (.)
We first compute λ 0 , applying again Lemma  to get back to the easy cylindrical neighborhood.

λ 0 = x 2 +y 2 ≤r 2 (z -z o ) 2 dxdy = x 2 +y 2 ≤r 2 z 2 dxdy + z 2 o x 2 +y 2 ≤r 2 dxdy -2z o x 2 +y 2 ≤r 2 zdxdy = x 2 +y 2 ≤r 2 z 2 dxdy + H 2 r 4 16 * πr 2 -2 Hr 2 4 r 4 4 πH + o(r 6 ) = 1 4 (k 2 1 x 2 +y 2 ≤r 2 x 4 dxdy + k 2 2 x 2 +y 2 ≤r 2 y 4 dxdy + 2k 1 k 2 x 2 +y 2 ≤r 2 x 2 y 2 dxdy - H 2 r 6 16 π + o(r 6 ) = 1 4 r 6 6 ( 3π 4 (k 2 1 + k 2 2 ) + k 1 k 2 π 2 ) - H 2 r 6 16 π + o(r 6 )
where H = k 1 +k 2 2 is the mean curvature. Thus

λ 0 = πr 6 32 ( k 2 1 + k 2 2 2 - 1 3 k 1 k 2 ) + o(r 6 ) (.)
Using . and . we get:

σ = r 2 16 ( k 2 1 +k 2 2 2 -1 3 k 1 k 2 ) + o(r 2 ) 1 + r 2 16 ( k 2 1 +k 2 2 2 -1 3 k 1 k 2 ) + o(r 2 )
which finally yields:

σ = r 2 16 k 2 1 + k 2 2 2 - 1 3 k 1 k 2 + o(r 2 )
The formula of the surface variation given by Theorem  indeed measures a sort of curvature. To interpret it we can notice that tel-00557589, version 1 -19 Jan 2011

• the surface variation is symmetric in k 1 ,k 2 ;

• in the case of a point lying on a sphere, k 1 = k 2 = k so σ sphere = r 6 24 k 2 ; • in the case of a saddle point k

1 = k = -k 2 , σ saddle = r 6
12 k 2 so σ sphere < σ saddle ;

• in the case of a cylinder k 1 = k, k 2 = 0, σ cylinder = r 6 32 k 2 . It follows from that the surface variation is not a discriminant information on the surface curvature, being unable to discriminate between very different local shapes.

. The MLS projection

We shall now explore a more efficient way of computing locally at least the mean curvature, and to use it as a scale space. The main tool of the scale space proposed in [DMMSL] is a simple projection of each surface point p on the local regression plane. This PCA regression plane is defined as the plane orthogonal to the least eigenvector of the centered local covariance matrix, and passing through the centroid of the neighborhood. The projection of p on this plane will be called p ′ . This projection method is the simplest instance of the moving least square method (MLS) by which each point of a surface is projected to a polynomial regression. The local barycenter can be considered as an MLS of order , MLS, and the present projection on a plane is an order  MLS, which we shall denote by MLS. The next lemma compares the normal to the PCA regression plane with the normal to the surface, n at p. Lemma . The normal v to the PCA regression plane in a spherical neighborhood B r at p ∈ M is equal to the surface normal at point p, up to a negligible factor: v = n + O(r).

Proof. The local PCA regression plane of point p is characterized as the plane passing through the barycenter of the neighborhood B r and with normal v minimizing:

I(v) = Br | v, pp ′ | 2 dp ′ s.t. v = 1 Denoting by (v x , v y , v z ) the coordinates of v, I(v) = Br (v x x + v y y + v z 1 2 (k 1 x 2 + k 2 y 2 ) + o(r 2 )) 2 dxdy.
Considering the particular value v = (0, 0, 1) shows that the minimal value I min of I(v) satisfies I min ≤ O(r 6 ). In consequence the minimum

(v x , v y , v z ) satisfies v x ≤ O(r) and v y ≤ O(r). Thus, since ||v|| = 1, v z ≥ 1 -O(r) and therefore v = n + O(r).

tel-00557589, version 1 -19 Jan 2011

By Lemma , projecting p onto the regression plane induces a motion which is asymptotically in the normal direction: p ′ p is almost parallel to n. Lemma  stated that sending each point onto the barycenter of its neighborhood approximates the mean curvature motion. We shall discuss later on why this result cannot be used for implementing the mean curvature motion. Nevertheless, the simple projection of each surface point p onto its local regression plane approximates a D scale space (mean curvature motion) as shown in the next theorem, and this discrete approximation will show very efficient.

→ pp ′ = λ(O(r), O(r), 1 -O(r)).
To compute λ, we use the fact that p ′ is the projection on the regression plane of p, and that o belongs to this plane by definition. This implies that -→ pp ′ ⊥ -→ op ′ and therefore

λ 2 O(r 2 ) + λ(1 -O(r))(H r 2 4 + o(r 2 ) + λ(1 -O(r))) = 0,
which yields λ = Hr 2 4 + o(r 2 ) and therefore

-→ pp ′ = (O(r 3 ), O(r 3 ), Hr 2 4 + o(r 2 )) = Hr 2 4 n + o(r 2 ).

. Asymptotics of MLS and MLS

There is some particular interest in MLS, because the recent meshing method uses it as a very simple and reversible smoothing tool for point clouds ( [DMMSL] and chapter ). On the other hand many cloud point processing methods involve some variant of the MLS method to smooth, interpolate, or sub-sample a point cloud. MLS and MLS are smoothing operators and therefore could be used as scale spaces, that is, as iterative smoothing operators. But, following [DMMSL] MLS indeed is a scale space and MLS is not, as illustrated in the experiments of Section .. MLS implements very accurately a mean curvature motion, and second that MLS is insensitive to first, second, and third order intrinsic derivatives and has an order  difference to the original surface. Thus, it is to be expected that by MlS a smooth surface will have only a very small motion. This is verified in the experimental section. We start with the asymptotic analysis of MLS.

Theorem . Consider a smooth manifold and its intrinsic coordinates at a point p(0, 0), so that the Taylor expansion in a neighborhood of p satisfies

z = f (x, y) = 1 2 k 1 x 2 + 1 2 k 2 y 2 + f 3 (x, y) + f 4 (x, y) + f 5 (x, y) + O(r 6 )
where f i are homogeneous polynomials in x, y of global degree i. The order  MLS approximation M LS1(p) of p in a radial neighborhood of p with radius r satisfies

< M LS1(p) -p, n >= r 2 8 (k 2 + k 1 ) + O(r 3 )
where x and y are the coordinates, k 1 and k 2 are the principal curvatures and n is the normal to the surface at p, oriented towards the concavity.

Lemma . One can choose the coordinates x and y in the regression plane at p so that, z being the coordinate in the direction of the normal plane, the equation of the manifold around p has the form z = f (x, y) = 5 i,j=1 a ij x i y j + o(|x 2 + y 2 | 3 ), and in addition a ij = ãij (1 + O(r)), where z = f (x, ỹ) = 5 i,j=1 a ij xi ỹj + o(|x 2 + ỹ2 | 3 ) is the equation of the manifold in the coordinates (x, ỹ, z) defined by the normal at p and the directions of the principal curvatures.

Proof. Consider (x, ỹ, z) the coordinates in the intrinsic frame such that x and ỹ are the coordinates associated with the principal curvatures at p, and the plane xpỹ is the tangent plane. Consider now coordinates (x, y, z) associated with the regression plane in a spherical neighborhood. Because the normal at the regression plane tends to the real normal when the spherical neighborhood shrinks, we can choose the coordinate axes (x, y) in the regression plane so that the rotation R which sends one frame to the other is close to the identity, namely Consider now the asymptotic expansion of z as a function of x, ỹ. (We assume the manifold to be at least C 5 ):

z -g(x, ỹ) -O((x 2 + ỹ2 ) 5 2 ) = 0.
Because of the relation (.) we can consider the above equation as an implicit equation in z, x, y, R, namely

Q(x, y, z, R) -O((x 2 + y 2 + z 2 ) 5 2 ) = 0. (.)
However, by the chain rule we have ∂Q ∂z (0, 0, 0, Id) = 1. Thus by the implicit function theorem, there is a function h of class C 5 such that in a neighborhood of (0, 0, 0, Id), (.) is equivalent to

z = h(x, y, R).
Since h is C 5 we can make a Taylor expansion and therefore get

z = g(x, y, R) + O(||R -Id|| 5 + (x 2 + y 2 ) 5 2 ).
In particular for R = Id we obtain by identification of the terms with degree lower than  that g(x, y, Id) = f (x, y). Thus, all monomials a ij (R)x i y j in the expansion of f with respect to x, y satisfy a i,j (R) = ãi,j (Id)+O(I -R), which by (.) yields a i,j (r) = ãi,j + O(r).

Proof of Theorem  In the local coordinate system, the surface can be described as the graph

z = f (x, y) = f 1 (x, y) + f 2 (x, y) + f 3 (x, y) + f 4 (x, y) + f 5 (x, y) + o(|x 2 + y 2 | 5/2 )
where f 1 (x, y) = a 10 x + a 01 y, f 2 (x, y) = a 20 x 2 + a 11 xy + a 02 y 2 f 3 (x, y) = a 30 x 3 + a 21 x 2 y + a 12 xy 2 + a 03 y 3 , f 4 (x, y) = a 40 x 4 + a 31 x 3 y + a 22 x 2 y 2 + a 13 xy 3 + a 04 y 4 ; f 5 (x, y) = a 50 x 5 + a 41 x 4 y + a 32 x 3 y 2 + a 23 x 2 y 3 + a 14 xy 4 + a 05 y 5 .

Let us find the order 1 polynomial g that best fits this surface in the least squares sense. Set g(x, y) = αx + βy + θ Then we must find the parameters Θ = α β θ) which minimize

x 2 +y 2 <r 2 (g(x, y)f (x, y)) 2 dxdy = where X = x y 1 . This is a quadratic minimization and the derivation yields

x 2 +y 2 <r 2 X T (XΘ Tf (x, y))dxdy = 0, and therefore

Θ T = x 2 +y 2 <r 2 (X T X) -1 x 2 +y 2 <r 2 (X T f (x, y)). Let us compute the matrix M = x 2 +y 2 <r 2 (X T X). M = x 2 +y 2 <r 2   x 2 xy x xy y 2 y x y 1   dxdy = πr 4 4   1 0 0 0 1 0 0 0 4 r 2   ; M -1 = 4 πr 4   1 0 0 0 1 0 0 0 r 2 4   .
We also have

X T f 1 (x, y) = πr 4 4   a 10 a 01 0   , X T f 2 (x, y) = πr 4 4   0 0 a 20 + a 02   X T f 3 (x, y) = πr 4 4   r 2 6 (3a 30 + a 12 ) r 2 6 (a 21 + 2a 03 ) 0   , X T f 4 (x, y) = πr 4 4   0 0 r 2 6 (3a 40 + a 22 + 3a 04 )   , X T f 5 (x, y) = πr 4 4  

tel-00557589, version 1 -19 Jan 2011

In this frame ã20 = 1 2 k 1 and ã02 = 1 2 k 2 . So that the final formula for the projection on the regression plane is:

( r 2 8 (k 2 + k 1 ) + r 4 24 (3ã 04 + ã22 + 3ã 40 ) + O(r 6 ))(1 + O(r)) = r 2 8 (k 2 + k 1 ) + O(r 3 )

.. The asymptotic behavior of MLS

We need to specify a simple version of MLS, prone to an asymptotic analysis, but in agreement with the literature. The following version seems to be a sort of common denominator of all versions. The first difficulty of MLS is the fact that a first reference frame must be found, and that then the mean square approximation by order  polynomials is made in this reference frame. The most natural frame is found by applying MLS, and the coordinates (x, y) is therefore the coordinates in the regression plane in a spherical neighborhood B r . The second step is to find the mean square closest order  polynomial in the spherical neighborhood. Because of Lemma  we can specify without loss of generality or precision that the minimization is made in the cylindrical neighborhood C r . In that way, all integrals computed in the approximation process are integrals on the disk x 2 +y 2 ≤ r 2 , which is numerically and formally convenient. Thus the MLS algorithm specified for the analysis works in the two steps:

. compute the regression plane of the manifold in the spherical neighborhood B r = B r (p) ∩ M;

. call (x, y) the reference coordinates in the regression plane. Consider the restriction of the smooth manifold to the disk D r := x 2 + y 2 ≤ r 2 , z = f (x, y). Then find the order  polynomial g(x, y) that best approximates f for the L 2 (D r ) distance;

. set (in the reference frame) M LS2(p) := (0, 0, g(0, 0)).

The next theorem shows that unlike MLS which reveals the mean curvature, the difference between a point smoothed by MLS and its original position is very small (of order ) and actually only reveals a fourth order intrinsic operator. The result is actually interesting because the revealed fourth order operator is a kind of bilaplacian. The evolution by MLS is a fourth order equation that is intuitively well-posed at least for short times. Viewed that way, the next theorem finds a very simple and efficient numerical scheme for at least one fourth order intrinsic evolution.

tel-00557589, version 1 -19 Jan 2011

Theorem . Consider a smooth manifold and its intrinsic coordinates at a point p(0, 0), so that the Taylor expansion in a neighborhood of p satisfies

z = f (x, y) = 1 2 k 1 x 2 + 1 2 k 2 y 2 + f 3 (x, y) + f 4 (x, y) + f 5 (x, y) + O(r 6 )
where f i are homogeneous polynomials in x, y of global degree i.

(x, y) = f 1 (x, y) + f 2 (x, y) + f 3 (x, y) + f 4 (x, y) + f 5 (x, y) + o(|x 2 + y 2 | 5/2 ) where f 1 (x, y) = a 10 x + a 01 y, f 2 (x, y) = a 20 x 2 + a 11 xy + a 02 y 2 f 3 (x, y) = a 30 x 3 + a 21
x 2 y + a 12 xy 2 + a 03 y 3 , f 4 (x, y) = a 40 x 4 + a 31 x 3 y + a 22 x 2 y 2 + a 13 xy 3 + a 04 y 4 , f 5 (x, y) = a 50 x 5 + a 41 x 4 y + a 32 x 3 y 2 + a 23 x 2 y 3 + a 14 xy 4 + a 05 y 5 .

We look for the order 2 polynomial g that best fits this surface in the least squares sense, g(x, y) = αx 2 + βy 2 + γxy + δx + εy + θ.

We therefore must find the parameters Θ = α β γ δ ε θ) which minimize

x 2 +y 2 <r 2 (g(x, y) -f (x, y)) 2 dxdy = x 2 +y 2 <r 2 (XΘ T -f (x, y)) 2 dxdy
where X = x 2 y 2 xy x y 1 . This is a quadratic minimization and differentiating this integral with respect to Θ yields

x 2 +y 2 <r 2 X T (XΘ Tf (x, y))dxdy = 0.

Writing M = x 2 +y 2 <r 2 X T X, finally:

Θ T = x 2 +y 2 <r 2 (X T X) -1
x 2 +y 2 <r 2 (X T f (x, y));

tel-00557589, version 1 -19 Jan 2011

Θ T = M -1
x 2 +y 2 <r 2 X T (f 1 (x, y)+f 2 (x, y)+f 3 (x, y)+f 4 (x, y)+f 5 (x, y)+O((x 2 +y 2 ) 3 );

X T X =        
x 4 x 2 y 2 x 3 y x 3 x 2 y x 2 x 2 y 2 y 4 xy 3 xy 2 y 3 y 2 x 3 y xy

. Numerical experiments

This section shows experiments with the most significant algorithms described in the previous sections. A simulated randomly sampled sphere plays the role of numerical pattern. In particular we evaluate the curvatures given by MLS on the sphere, and the curvatures given by MLS projection, by polynomial regression. We also compute the curvature estimated by the method described in [BC] and by the surface variation of [PGK]. The results are compared by giving the mean tel-00557589, version 1 -19 Jan 2011 estimated curvature and its standard variation. The input data is a randomly sampled sphere with radius 2 corrupted with added centered Gaussian noise of variance 0.1. By comparing the values in tab. ., two conclusions can be drawn: first, MLS is much more stable than the MLS projection, which can be observed by the standard variation on the estimates. Second, the MLS projection does not implement a mean curvature motion. The MLS projection instead yields an increase of the mean curvature (i.e., the sphere radius decreases, which is expected from a mean curvature motion). The standard deviation of the radius of the sphere gets smaller by MLS than by MLS, meaning that MLS is a better curvature estimator. The surface variation of this point set is 0.1419 with standard variation 0.009. Finally, the mean curvature computed by normal covariance analysis is 2.4565 with variance 0.1026.

Iteration

On real surfaces (Figs ., . and .), we present various curvature distribution and surface variation distributions. Intuitively, a mean curvature motion corresponds well to the MLS projection.

To better judge of the smoothing effect of the projection on the regression plane operator, we show the following experiment. First a consistently oriented point set is built (see [DMMSL] or chapter  for an efficient way of doing so). This normal orientation yields the sign of the mean curvature, by computing the scalar product of the oriented normal and the displacement vector. Each point is then plotted in a different color according to its sign, blue for positive and red for negative (see Fig .). This experiment shows that, at the beginning, the sign captures noise and small textures (i.e., small scale variations), and after some iterations, the shape is smoothed and the sign captures the geometry of the shape (large scale variations).

tel-00557589, version 1 -19 Jan 2011

Conclusion

In this chapter, we analyzed several proposition of discrete operators computing curvatures or curvature-like operators computed on surfaces defined by raw point tel-00557589, version 1 -19 Jan 2011 sum of weights is constant over all discrete Euclidean neighborhoods with a fixed size? To the best of our knowledge, there is no answer to this question yet, although this problem appears crucial to the field of numerical surface motion. Finally the simplest analysis performed in this chapter (the analysis of the MLS projection) yields a mean curvature motion, and this operator, proven very robust to irregular sampling, has already multiple practical applications (see chapters  and  which are detailed versions of [DMMSL] and [DMAL]).

Introduction

Color images are functions defined on Ω ⊂ R 2 with values commonly clipped in [0, 255] 3 in the RGB color space. By dissociating pixel values from pixel positions, we get a set of D values in the color cube [0, 255] 3 . These values constitute the color histogram of the image, a "color cloud" that will be filtered by the algorithm.

The idea behind the projection on the regression plane is that the underlying structure of the color cloud is essentially 2-dimensional, otherwise the whole procedure would be pointless. The question of the dimensionality of color clouds has been investigated: [OW] proposed a linear model to account for the color clouds structures. In [BLM], it is shown that the D model accounted much better for the colors in the image. By a statistical analysis and numerical comparison between a D filtering and D filtering, this paper confirms that the color clouds are essentially D, roughly up to a % error. In other terms the clouds can be modeled as a (thick) D manifold. By the projection filter, the color clouds can be filtered to become really D surfaces, whose visualization is no more obscured by a few stray color points. Here, thanks to a careful numerical analysis, we were able to test this D hypothesis on far larger color images (2048 × 1536 or 3648 × 2736 pixels) than those tel-00557589, version 1 -19 Jan 2011

.. Color clouds



studied in [BLM]. In addition, the filtered color clouds could then be meshed by scale space meshing to get a visualization of the color sheets at any smoothing scale.

. Color clouds

Each pixel stores a five dimensional information (2 parameters for space information and 3 parameter for color information: red, green and blue channel). As stated before, if we consider only the color information, then we have a D point cloud contained in the ([0, 255] 3 ) cube. An example of a natural image and its color cloud can be seen on figure ..

. Implementation

The main difference with usual point cloud surfaces coming from triangulation scanners is that such data bring us close to the complexity worst case. Indeed the points are much more spread out in the color cube than they are in a good quality scan. This means that many more octree cells are actually occupied. As a result an acceleration of the code must be introduced. A difference is that the octree size stays always the same: S = 255 (in fact we set it to 265 to allow for a small error). Given a filtering radius r specified by the user (which according to [BLM] goes from 5, 10 to 20), an octree depth d can be inferred so that the leaf size l is:

S 2 l+2 < r < S 2 l+1
tel-00557589, version 1 -19 Jan 2011

Because of the order of magnitude of the radius, this would mean an octree depth of less than 4, and would therefore make the neighbor search intractable. In order to accelerate the search, we use the fact that we do not need the neighbors individually. Indeed, we only need to compute their covariance matrix and barycenter. If each cell (be it leaf or not) stores its barycenter, number of points and un-centered covariance, then each time the cell is entirely included inside the ball centered at the query point p and with radius r, the cell barycenter and coefficients should be considered instead of computing the covariance from the set of neighboring points. More specifically, when computing the local regression plane, we only need to compute the neighborhood covariance Σ and the barycenter o from the set of N r neighboring points. The covariance being a symmetric matrix, we only store its upper coefficients.

When building a tree, values must be stored in each cell. Each time a point p is added to the octree, each cell C it traverses when going down in the tree until reaching the leaf, updates part of the un-centered covariance matrix A and the unnormalized barycenter q: a 00 = a 00 + p. and unnormalized barycenter coefficients q.x = q.x + p.x q.y = q.y + p.y q.z = q.z + p.z, and the number of points N = N + 1. When a query point p is given and we look at cells intersecting its r-ball, then the following alternative decisions can be made

• If the cell is included in the ball neighborhood, the coefficients of A (respectively q) are added to the coefficients of Σ (resp o) and the number of neighbors N r = N r + N ;

• If the cell is not included in the ball neighborhood and does not intersect the ball neighborhood, discard it;

tel-00557589, version 1 -19 Jan 2011

.. Results



• If the cell is not included in the ball neighborhood but intersects the ball then:

-If the cell is a leaf then look individually at the points it contains and update Σ, o, N r accordingly;

-If the cell is not a leaf then process all of its children cells using the same method.

At the end we have the upper coefficients of the un-centered covariance matrix Σ, the sum of all neighbor coordinates O and the number of neighbors and the final values can be deduced as Remembering that Σ is symmetric, we obtain the final centered covariance and barycenter. This simplification explains why when using a high depth value, the computation time can be reduced. Its drawback is that high depth values require a larger memory use (by increasing the octree depth by 1, 8 times more pointers are created!).

. Results

Result

Conclusion

This chapter verified on numerous examples that applying the regression plane projection operator to a color cloud did not change a lot the perception of the image tel-00557589, version 1 -19 Jan 2011

Conclusion

In this thesis, a whole point cloud processing pipeline from the acquisition to the final mesh segmentation and extraction was studied. Yet the field of high precision shape representation is still wide open. Here we summarize some of the possible improvements and future work.

The method introduced proved to be very versatile because it was not only able to process high precision surfaces, but was also able to deal with more noisy data such as color clouds. In the continuation of this work, different types of data and improvements could be made to get a still better virtual representation of the acquired object. It seemed to us that, little by little, this work was clearing the road toward the high precision visualization and exploration of fine arts and archaeological objects acquired by a high precision laser. The field of cultural heritage digitization is already a very active (through the d-coform project  for example), but not to the precision we now get. This work opens, in some sense, the hope of a digital conservation and of an improved and more versatile display to amateurs and professionals of these sometimes fragile small pieces In particular the digital magnification seems to be an efficient visualization tool for high precision scans.

. Another type of data: D from stereo

The existence of passive light acquisition device was already mentioned in this work. From a set of stereo-rectified images, one can extract points on the surface of the object, either by calibrating one camera and rectifying several images, or by using directly a stereo camera. An interesting goal would be to build the stereo reconstruction for the same objects that we used with the laser scanner. Then precisions of active and passive light acquisition devices could and should be compared and, hopefully, lead to technologically very needed cross-calibration tools of D devices.

. High precision registration

Our laser scanner was somewhat limited in terms of the size of the acquired objects. In order to be still able to acquire larger objects with the same precision, Inverse Geometry: From the raw point cloud to the D surface Theory and Algorithms Abstract:

Many laser devices acquire directly D objects and reconstruct their surface. Nevertheless, the final reconstructed surface is usually smoothed out as a result of the scanner internal de-noising process and the offsets between different scans.

This thesis, working on results from high precision scans, adopts the somewhat extreme conservative position, not to loose or alter any raw sample throughout the whole processing pipeline, and to attempt to visualize them. Indeed, it is the only way to discover all surface imperfections (holes, offsets). Furthermore, since high precision data can capture the slightest surface variation, any smoothing and any sub-sampling can incur in the loss of textural detail.

The thesis attempts to prove that one can triangulate the raw point cloud with almost no sample loss. It solves the exact visualization problem on large data sets of up to 35 million points made of 300 different scan sweeps and more. Two major problems are addressed. The first one is the orientation of the complete raw point set, an the building of a high precision mesh. The second one is the correction of the tiny scan misalignments which can cause strong high frequency aliasing and hamper completely a direct visualization.

The second development of the thesis is a general low-high frequency decomposition algorithm for any point cloud. Thus classic image analysis tools, the level set tree and the MSER representations, are extended to meshes, yielding an intrinsic mesh segmentation method.

The underlying mathematical development focuses on an analysis of a half dozen discrete differential operators acting on raw point clouds which have been proposed in the literature. By considering the asymptotic behavior of these operators on a smooth surface, a classification by their underlying curvature operators is obtained.

This analysis leads to the development of a discrete operator consistent with the mean curvature motion (the intrinsic heat equation) defining a remarkably simple and robust numerical scale space. By this scale space all of the above mentioned problems (point set orientation, raw point set triangulation, scan merging, segmentation), usually addressed by separated techniques, are solved in a unified framework.

Keywords: D surface reconstruction, high precision point clouds, point cloud scale space tel-00557589, version 1 -19 Jan 2011

  Figure .: Comparison between several meshing methods on a 1 cm high logo. See chapter  for more details

  Figure .: The scan merging problem and its solution . 38 scans cover the Tanagra head, represented in a different colors in .(a).

Figure

  Figure .: Line extraction (left) and segmentation (right) obtained by Mesh-MSER

Figure

  Figure .: An image and two views of its filtered color cloud

  Figure .: A Peugeot car acquired by a range laser scanner (acquisition made by Délégation Générale de l'Armement).

Figure

  Figure .: Picture of the LURPA scanner laser and its output: a raw point set

Figure

  Figure .: Registering multiple views (Image from David Laser Scanner Website)

  Figure .: Level set method for reconstructing a surface.

Figure

  Figure .: Scheme of ray tracing (copyright Henrik, Wikipedia)

  Figure .: Comparison between cylindrical and spherical neighborhoods

  Finally we have Br f (x, y)dxdy = Cr f (x, y)dxdy + O(r 4+n ). (.) Combining (.), (.) and (.) yields (.).

  Figure .: Comparison of the iterated barycenter and of the iterated projection filter on a randomly sampled sphere. Both motions are consistent with the mean curvature motion, but the iterated barycenter provokes clustering. See text for an explanation.

Figure

  Figure .: Denoising a noisy circle with 1, 50, 100 iterations of T r (top) and MLS (bottom). Even after 100 iterations the oscillations removed by T r persist with MLS. The sphere radius decreases with T r , which is consistent with the mean curvature. tel-00557589, version 1 -19 Jan 2011



  Scale space OrientationData: A point cloud M S , a radius r, an update parameter α > 1 Iterate the projection filter T r and keep track of each raw data point sample (mean curvature motion); Find a point p 0 in a flat area, pick its orientation and mark it as oriented. Add its neighbors to the stack S;  while S is not empty or S does not become constant do  Take p 0 the first point in S;  if orientateFromNeighbors(p,r,t) then  Mark the point as oriented and remove p from S;  Add p's neighbors to S;  Add all remaining unoriented points to S;  while S is not empty and ♯S does not become constant do  r = αr;  for p in S do  Perform orientateFromNeighbors(p,r,t);

Figure

  Figure .: A raw point set (left) and its orientation (right). Points in the right figure are simply given a gray value depending on the scalar product of their normal and the lighting orientation



  Figure .: D example of the steps performed by the scale space meshing algorithm

  Fig. .(f) shows the comparative result on one of the finest details of the Tanagra data set. The results are similar in terms of detail quality, yet the computation time was doubled, and we have seen (Fig. .) that MLS does not deliver a scale space and keeps the smoothed out noise. Fig.

  Fig .. Fig. . shows the scale space reconstruction of one scan of a fine scale object (i.e. the mesh back-projected at all scales).Fig . compares the mesh reconstruction by several meshing algorithm with scale space meshing. The experiment clearly rules out both Ball Pivoting algorithm and Poisson Reconstruction. Two MLS methods are also tested. APSS ([

  Figure .: Multi-resolution mesh reconstruction of the Tanagra point set (22 cm high) illustrating the recovery of fine texture. All back propagated textures are present on the original tel-00557589, version 1 -19 Jan 2011

  Figure .: Comparison between several meshing methods. The width of this logo is approximately 1cm. The direct mesh (.(c)) creates many spurious triangles. The Poisson reconstruction [KBH] clearly smoothes out all details (.(d)). Filtering the logo by order two MLS and meshing the points by the ball pivoting algorithm (.(e)) also creates a smooth mesh. Compare the details in .(b) and .(c). Fig .(f) shows the result of applying the same scale space strategy with the projection on the order 2 MLS surface instead of the regression plane projection operator. The result is alike in detail quality but the computation time is double. See figs .-. for an explanation of this difference. tel-00557589, version 1 -19 Jan 2011

  Figure .: Comparison of the Rosette reconstruction (Picture (a)) using Ball Pivoting Algorithm (b), Poisson Reconstruction (c), RIMLS (d), APSS (e), and scale space meshing (f). APSS and RIMLS yield results that are really close to ours, yet both methods need an isosurface extraction done with the marching cubes, which creates strong artefacts (see a closeup Fig .). Besides, RIMLS and APSS meshes contain around 268500 vertices whereas the scale space mesh contains 132203 vertices. Notice also that APSS and RIMLS introduce some denoising (visible especially in the nearly flat parts). Scale space meshing is the only method that preserves exactly the input data. tel-00557589, version 1 -19 Jan 2011

  Figure .: Detail of the mesh built using (from left to right) RIMLS, APSS and scale space meshing. Notice the horizontal artefacts created by the isosurface extraction for both APSS and RIMLS methods.

Figure

  Figure .: Comparison on a piece of the Fragment Urbis Romae (FUR) database. Texture and details are better recovered on the back-propagated mesh (middle). Compare with the VRIP reconstruction method available on the FUR website (right)

  Fig. . illustrates the problem. With two overlapping shifted scan grids, as seen in (a), the aliasing risk is high. Meshing each scan separately yields two almost identical surfaces and textures (b, c). Nevertheless, a joint meshing (d)

  Figure .: Example of two overlapping scans (a), points of each scan are first meshed ((b)-(c)) separately. The result can be compared to the meshing of points of both scans together (d)

Figure

  Figure .: The unmerged head with aliasing artifacts (left), its smooth base (middle) and the merged result (right)

  Figure .: Noise estimates on each separate scan A and B before and after their merging.

  Figure .: Two noisy sine functions before and after merging.

  Figure .: Two noisy edges before and after merging.

Figure

  Figure .: Merging of the mask scans seen from the back side (top row) and left side (bottom row). Left: picture, middle: without merging, right: with merging tel-00557589, version 1 -19 Jan 2011

Figure

  Figure .: Comparisons of the merging (a) with a level set reconstruction method ([KBH]) of the unmerged scans point set (b) and a filtering of the unmerged scans point set (c).The level set method obviously introduces a serious smoothing, yet does not eliminate the scanning boundary lines. The bilateral filter, applied until all aliasing artifacts have been eliminated, over-smoothes some parts of the shape.

Figure

  Figure .: Comparison of registration of two scans (colored in different colors on the first figure) using Global Non Rigid Alignment [BR] and scale space merging.Meshes were reconstructed using[KBH].

  Figure .: Comparison of the mesh obtained by merging and by Poisson reconstruction on a detail of nefertiti's cheek. In this case, Poisson Reconstruction suppresses registration artefacts but smoothes out the details.

FigureFilling

  Figure .: Details of the Rosetta object without merging (top row) and with fusion (bottom row) tel-00557589, version 1 -19 Jan 2011

Algorithm  : 

 : Patch refinement Data: A set of holes with patches Result: Refined patches  changed ← true;  while changed do  changed ← f alse; Compute the scale attributes of hole boundary points;  for each triangle t = (v i , v j , v k ) in the patch do  Compute the triangle centroid c and its scale attribute σ(c);

Figure

  Figure .: Initial hole (left), intial patch (middle) and refined patch (right).

  Figure .: Thin plate spline minimization iteration , and .

Figure

  Figure .: Thin plate spline minimization iteration ,  and .

Figure

  Figure .: Thin plate spline minimization iteration ,  and .

  Figure .: first test object, a mire (left) and the same object with filled holes (right), here the base is not removed

Figure

  Figure .: Dancer with Crotales: fine scale mesh (left), fine scale mesh with filled holes (right) , here the base hole patch is removed tel-00557589, version 1 -19 Jan 2011

Figure

  Figure .: Mask: coarse scale mesh (left), coarse scale mesh with filled holes (middle), fine scale mesh with filled holes (right) . Here the base hole patch is removed.

Figure

  Figure .: Computation times for various models . All computation times are given for filling the holes in the fine scale mesh

  Figure .: Example of an image (left) and its level set tree (right). Node A is the father of B and B contains a hole tel-00557589, version 1 -19 Jan 2011

Figure

  Figure .: Classification by Mesh-MSER Analysis of a diamond shaped pattern (k triangles).

Figure

  Figure .: Classification by MSER Analysis (From top to bottom: picture of the object; obtained mesh; MSER segmentation; MSER borders)

  Figure .: Maximally stable curvature level lines of "La dame de Brassempouy". Some lines on the above figure appear to be open because of the cropping into front and back part)

  Figure .: Picture (a) of fragment u of Stanford FUR database, D mesh (b),Mesh-MSER result (d). This result can be compared with the D-MSER result (b) on a good quality picture of the same fragment[DMM], and with the result of a Canny edge detector applied to the same picture (c). The comparison gives a sweeping advantage to Mesh-MSER. Indeed, D-MSER misses parts and keeps noisy level lines. Canny's detector has many outliers and yields anyway no segmentation. Similar experiments on artificial renderings of the mesh gave no better results tel-00557589, version 1 -19 Jan 2011

Figure

  Figure .: Rendering of a digital elevation model .

Figure

  Figure .: Result of Mesh-MSER on the pump mesh and comparison with the method by Zatzarinni et al. ([KST]).

Lemma  .

 . Integrating on M any function f (x, y) such that f (x, y) = O(r n ) on a cylindrical neighborhood C r instead of a spherical neighborhood B r introduces an o(r n+4 ) error. More precisely:Br f (x, y)dm = x 2 +y 2 <r 2 f (x, y)dxdy + O(r 4+n ). (.)  We could use z = f (x, y) = -1 2 (k 1 x 2 + k 2 y 2 ) + o(x 2 + y 2) at the cost of changing the orientation and sign of k 1 ,k 2 .tel-00557589, version 1 -19 Jan 2011 Chapter . A numerical analysis of raw point cloud smoothingProof. The surface area element of a point m(x, y, z(x, y)) on the surface M, expressed as a function of x, y, dx and dy is dm(x, y) = 1 + z 2x + z 2 y dxdy. One hasz x = k 1 x + O(r 2 ) and z y = k 2 y + O(r 2 ). Thus dm(x, y) = (1 + k 2 1 x 2 + k 2 2 y 2 + O(r3))dxdy which yields dm(x, y) = (1 + O(r 2 ))dxdy. (.) Using (.), the integrals we are interested in become Cr f (x, y)dm = (1 + O(r 2 )) Br f (x, y)dxdy; (.) Br f (x, y)dm = (1 + O(r 2 )) Cr f (x, y)dxdy (.)

  y)ρdρdθ -E = Cr f (x, y)dxdy -E, tel-00557589, version 1 -19 Jan 2011 .. Regression-free curvature estimates  with E=: [0,2π] [ρ(θ),r] f (x, y)ρdρdθ. Thus |E| ≤ 2π sup x 2 +y 2 ≤r 2 |f (x, y)| max(k 2 1 , k 2 2 )r 4 + O(r 5 ).In particular if f (x, y) = O(r n ), then |E| ≤ O(r 4+n ). Finally we have Br f (x, y)dxdy = Cr f (x, y)dxdy + O(r 4+n ). (.) Combining (.), (.) and (.) yields (.).

Theorem  .

 . Let T r be the operator defined on the surface M transforming each point p into its projection p ′ = T r (p) on the local regression plane. Then T r (p)p = Hr 2 4 n + o(r 2 ). (.) Proof. By Lemma  the barycenter o of B r has local coordinates -→ po = (o(r 2 ), o(r 2 ), Hr 2 4 + o(r 2 )). On the other hand -→ pp ′ is proportional to the normal to the regression plane, v. Thus by Lemma  -

  (x, ỹ, z) = R(x, y, z) (.) is close to the identity: R → Id when r → 0. More precisely, by Lemma , the normal v(r) to the PCA regression plane in a spherical neighborhood B r at p ∈ M is equal to the surface normal at point p, up to a negligible factor: v(r) = n + O(r). Thus we also have R = R(r) = I + O(r). (.) tel-00557589, version 1 -19 Jan 2011

  x 2 +y 2 <r 2(XΘ Tf (x, y)) 2 dxdy,tel-00557589, version 1 -19 Jan 2011  Chapter . A numerical analysis of raw point cloud smoothing

Figure

  Figure .: Curvature evolution by iterative projection on MLS

  Figure .: Other curvature estimates .

Figure

  Figure .: Evolution of the motion direction with projection iterations.

  Figure .: Example of a color cloud .

x 2 a

 2 01 = a 01 + p.x • p.y a 02 = a 02 + p.x • p.z a 11 = a 11 + p.y 2 a 12 = a 12 + p.y • p.z a 22 = a 22 + p.z 2

  o.x = o.x/N r o.y = o.y/N r o.z = o.z/N r Σ 00 = Σ 00 -N r • o.x 2 Σ 01 = Σ 01 -N r • o.x • o.y Σ 02 = Σ 02 -N r • o.x • o.z Σ 11 = Σ 11 -N r • o.y 2 Σ 12 = Σ 12 -N r • o.y • o.z Σ 22 = Σ 22 -N r • o.z 2

  of color cloud filtering are shown on figs ., ., ., ., ., ., ., ., .. One can see on fig . that the perception of the image barely changes, over the iterations. By looking closely, however, one can notice that some colors disappeared (fig. .).

  Figure .: Evolution of image "Turin": original image and image evolution (top row) and two views of its color cloud (two last rows) after 1 and 10 iterations.therefore comforting the theory of[BLM]. The filter is very fast when applied to small size images (below 1000 × 1000), yet when the image size increases, a good compromise must be found between time efficiency and memory use. A useful addition to[BLM] is that now the noisy color cloud can even be meshed at all scale by scale space meshing. Experimental evidence shows that the mesh is already visible after only one filtering iteration. tel-00557589, version 1 -19 Jan 2011

FigureFigure

  Figure .: Evolution of a color image with iterative projection filter (radius ) (top row) and color cloud (bottom row)

FigureFigureFigure

  Figure .: Comparison of the filtering between  iterations of the projection filter with radius 10 and 20.

  The next technical lemma justifies its use.

	Lemma . Integrating on M any function f (x, y) such that f (x, y) = O(r n ) on
	a cylindrical neighborhood C r instead of a spherical neighborhood B r introduces an
	O(r n+4 ) error. More precisely:
	f (x, y)dm =	f (x, y)dxdy + O(r 4+n ).	(.)
	Br	x 2 +y 2 <r 2
	tel-00557589, version 1 -19 Jan 2011

 Chapter . Scale Space Meshing of Raw Data Point Sets

	This yields ρ(θ) = r -1 8 k(θ) 2 r 3 + O(r 3 ). We shall use this estimate for the error
	term E appearing in
	Br
	tel-00557589, version 1 -19 Jan 2011
	2	.

  Table. gives the computation times for various shapes of various sizes with varying numbers of scans. Notice the high number of scans necessary to get a good covering of the object. It entails that several dozens of scans have to be merged on the more exposed parts.

	Point set	points	scans Time(s) height
	Dancer	5, 524, 627	94	321	17cm
	Mime	8, 611, 522	102	140	11cm
	Greek Mask	8, 961, 736	78	106	12cm
	Nefertiti	15, 554, 528 115	819	18cm
	Tanagra	17, 496, 999 160	1258	22cm
	Rosetta	36, 201, 537	32	45min 30cm
	Figure .

: Computation time for the proposed merging. It is significantly faster than the scanning time itself Figs .

Model Number of Number of Number of Computation time points

			triangles	holes
	Crotales 5, 496, 386 10, 988, 071	27	246.13s
	Mask	8, 938, 478 17, 874, 326	108	548.66s
	Mime	8, 571, 061 17, 139, 634	30	270.08s
	Mire 	1, 436, 303	2, 866, 460	117	76.87s
	Tanagra 16, 211, 676 32, 385, 734	630	37min20s

  One calls extremal region any connected tel-00557589, version 1 -19 Jan 2011 component of some level set. For each connected component N k of the level set E k , either k = k min , in which case N k is the whole image domain, or there exists a connected component N i of the upper level set

  tagged with their levels Result: The component tree forest {F(M)} by Algorithm  gives the graph of connected components of level sets of H a mapping defined on M. This graph is a forest (meaning that the constructed graph has no cycles).

 Compute the levels of all triangles;  Sort F in decreasing level order;  Set all triangle markers to active;  for t ∈ F and t is active do

 k = t → level;  Create an empty node N new ;  E = {t};  while E is not empty do  Remove and return the first element t of E;  Get t 1 , t 2 , t 3 the neighbors of t;  for i = 1 • • • 3 do  if t i is active and t i → level == k then  Add t i to the set E;  if t i → level > k then  Get the node N containing t i ;  Get N a the last built ancestor of N ;  if N a → level = k then  Merge node N a into N new ;  else  N a .f ather = N new ;  Add N a to N new 's children;  |N new | ← |N new | + |N a |;

 Mark t as inactive;  Add t to the set of triangles of N new ;  |N new | ← |N new | + Area(t); Theorem . Assume that the mesh (V, T ) is a manifold (meaning that each edge is shared by at most two triangles). Then the list of nodes and father-child relations tel-00557589, version 1 -19 Jan 2011 created

Tools for numerical analysis of point cloud surface motions 

  Let p(x p , y p , z p ) be a point of the surface M. At each non umbilical point p, consider the principal curvatures k 1 and k 2 linked to the principal directions t 1 and t 2 , with k 1 > k 2 where t 1 and t 2 are orthogonal vectors. (At umbilical points, any orthogonal pair (t 1 , t 2 ) can be taken.) Set n = t 1 × t 2 so  Some theorems are in fact already in chapter , we rewrite them for clarity. Comparison between cylindrical and spherical neighborhoods that (t 1 , t 2 , n) is an orthonormal basis. The quadruplet (p, t 1 , t 2 , n) is called the local intrinsic coordinate system. In this system we can express locally the surface as a C 2 graph z = f (x, y). By Taylor expansion,

	Cylindrical Neighborhood
	M
	P	Regression Plane
	Spherical Neighborhood
	Figure .:
	tel-00557589, version 1 -19 Jan 2011

..

  12 + 3a 30 ) + r 4 16 (5a 50 + a 32 + a 14 ) + O(r 5 ) a 01 + r 2 6 (3a 03 + a 21 ) + r 4 16 (a 41 + a 23 + 5a 05 ) + O(r 5 ) (a 20 + a 02 ) r 2 4 + r 4 24 (3a 04 + a 22 + 3a 40 ) + O(r 6 )

					r 4 16 (5a 50 + a 32 + a 14 ) r 4 16 (a 41 + a 23 + 5a 50 )	 
					0
	Thus the parameter Θ satisfies
	Θ = M -1 X T	 a 10 + r 2 6 (a  f (x, y) =   .
		x 2 +y 2 <r 2
	This means that the projection on the plane generates a motion of amplitude
		r 2 4	(a 20 + a 02 ) +	r 4 24	(3a 04 + a 22 + 3a 40 ) + O(r 6 ),

Finally, Lemma  permits to replace a ij by ãij (1 + O(r)) where ãij are the coefficients in the frame defined by the principal curvature eigenvectors and the normal.

  ∂ 2 y are the fourth derivatives of the intrinsic equation at p in the directions of x, y and x, y respectively, and n is the normal to the surface at p, oriented towards the concavity.

			The order  MLS
	approximation M LS2(p) of p in a cylindrical neighborhood of p with radius r satisfies
	< M LS2(p) -p, n >= -	r 4 48	(3a 04 + a 22 + 3a 40 )) + O(r 5 )
	where x and y are the coordinates of principal curvatures, a 40 = 1 4! 1 4! ∂ 4 f ∂ 4 f ∂y 4 , a 22 = 1 4! ∂x 2	∂ 4 f ∂x 4 , a 04 =

Proof. Let us write f

  3 x 2 y 2 x 2 y xy 2 xy x 22 + 6a 40 ) + O(r 4 ) a 02 + r 2 8 (a 22 + 6a 04 ) + O(r 4 ) a 11 + 3r 2 8 (a 13 + a 31 ) + O(r 4 ) a 10 + r 2 6 (a 12 + 3a 30 ) + r 4 16 (5a 50 + a 32 + a 14 ) + O(r 5 ) a 01 + r 2 6 (3a 03 + a 21 ) + r 4 16 (a 41 + a 23 + 5a 05 ) + O(r 5 ) -r 4 48 (3a 04 + a 22 + 3a 40 ) + O(r 6 )

							
							
							
							
				3	xy 2 x 2 y x 2 xy x	 
				x 2 y	y 3	xy 2 xy y 2 y	 
				x 2	y 2	xy	x	y	1
			a 20 + r 2 8 (a 
							
							
							
							
							
							
							
							
	so that the MLS projection satisfies
	g(0, 0) = -(3a Finally Lemma  permits to replace r 4 48
			g(0, 0) = -	r 4 48	(3a 04 + a 22 + 3a 40 ) + O(r 6 ).
	by
	-(	r 4 48	(3ã 04 + ã22 + 3ã 40 ) + O(r 6 ))(1 + O(r)) = -	r 4 48	(3ã 04 + ã22 + 3ã 40 ) + O(r 5 )

04 + a 22 + 3a 40 ) + O(r 6 ).

  Comparison of the curvature estimation by iteration of the MLS projection and iterations of the MLS projection

			MLS		MLS
		mean standard variation mean standard variation
	0	0.5828	2.8609	0.052	1.2879
	1	0.5158	1.2434	0.4920	1.0053
	2	0.5079	0.3196	0.5083	0.1259
	3	0.5102	0.0253	0.5073	0.1001
	4	0.5136	0.0189	0.5068	0.0855
	5	0.5171	0.0165	0.5065	0.0749
	10	0.5356	0.0156	0.5058	0.0489
	Figure .:

  http://www.dcoform.eu/ tel-00557589, version 1 -19 Jan 2011 [BLM] Buades T., Lisani J., Morel J.: On the distribution of colors in natural images. preprint MAP, February . , , ,  [BM] Besl P. J., McKay N. D.: A method for registration of -d shapes. IEEE Trans. Pattern Anal. Mach. Intell. ,  (), -. ,  [BMR * ] Bernardini F., Mittleman J., Rushmeier H., Silva C., Taubin G.: The ball-pivoting algorithm for surface reconstruction. IEEE TVCG  (), -. , , , ,  [BNK] Borodin P., Novotni M., Klein R.: Progressive gap closing for mesh repairing. In Advances in Modelling, Animation and Rendering, Vince J., Earnshaw R., (Eds.). Springer Verlag, July , pp. -.  [Boi] Boissonnat J.-D.: Geometric structures for three-dimensional shape representation. ACM Trans. Graph. ,  (), -.  [BPK Bernardini F., Rushmeier H., Martin I., Mittleman J., Taubin G.: Building a digital model of michelangelo's florentine pieta. -. printed (dossier recalage).  [BS] Barequet G., Sharir M.: Filling gaps in the boundary of a polyhedron. Comput. Aided Geom. Des. ,  (), -.  [BS] Benjemaa R., Schmitt F.: Fast global registration of d sampled surfaces using a multi-z-buffer technique. In Image and Vision Computing (), pp. -.  [BSKa] Bendels G. H., Schnabel R., Klein R.: Detail-preserving surface inpainting. In The th International Symposium on Virtual Reality, Archaeology and Cultural Heritage (VAST) (Nov. ), Eurographics Association, Eurographics, pp. -.  tel-00557589, version 1 -19 Jan 2011 [ZZW * ] Zeng W., Zeng Y., Wang Y., Yin X., Gu X., Samaras D.: d nonrigid surface matching and registration based on holomorphic differentials. In ECCV ': Proceedings of the th European Conference on Computer Vision (Berlin, Heidelberg, ), Springer-Verlag, pp. -. 

	tel-00557589, version 1 -19 Jan 2011 tel-00557589, version 1 -19 Jan 2011

 * ] Botsch M., Pauly M., Kobbelt L., Alliez P., Lévy B.: Geometric modeling based on polygonal meshes. In Eurographics Tutorial ().  [BR] Brown B., Rusinkiewicz S.: Non-rigid range-scan alignment using thin-plate splines. In DPVT' (). printed. , ,  [BR] Brown B., Rusinkiewicz S.: Global non-rigid alignment of -D scans. ACM Transactions on Graphics (Proc. SIGGRAPH) ,  (Aug. ). , , , ,  [BRM * ]

(a) Without Merging (b) With Merging

Remerciements

Acknowledgements The authors thank the Stanford Digital Forma Urbis Romae Project (http://formaurbis.stanford.edu) and professor Marc Levoy for the fragments data used in figures ., . a joint property of the Sovraintendenza of Rome and of Stanford University. tel-00557589, version 1 -19 Jan 2011 Acknowledgements Pictures and scans of the Stanford Forma Urbis Romae are used with permission of Professor Marc Levoy (Stanford Digital Forma Urbis Romae Project). The vase model is property of Laboratory of Computer Graphics & Multimedia at the Technion and the Zinman Institute of Archaeology, University of Haifa. tel-00557589, version 1 -19 Jan 2011 tel-00557589, version 1 -19 Jan 2011

we would need to move the object between different acquisitions. All the acquired parts should then be registered. The problem of shape registration has been studied already, yet in most cases the mesh was then built by the Poisson method, which actually hides the registration artifacts rather than it corrects them. A reasonable question is: can we register point sets well enough so that our reconstruction method can be used? Many good descriptors have been proposed to perform the registration and this could lead to test their precision.

. Color reprojection

If the final goal is a high precision digital reproduction of art pieces, in other terms a digital Tanagra process, the visual information (the object's texture) will have to be fused with its D shape. It is therefore a stimulating problem to realize a very accurate reprojection of the color on the virtual object from a set of pictures. However, the pictures should be first registered to the D object, leading to another type of registration: the D-D registration, which has already been investigated but never for data of this accuracy. The color information should also not be perturbed by shadows nor by geometric textures (such as those are acquired by the laser).

A. Classifying points

All along this work an efficient fixed-radius neighbor finding algorithm was needed. Indeed this is one of the trickiest problems in surface processing: given a point p find the coordinates of all samples lying within a given distance of p. Because of sampling irregularities it may happen that this neighborhood contains either no other point than p or the whole set.

The very naive implementation would need to traverse once the whole set of N points for each of the N points, yielding an algorithm with N 2 complexity. In case of point sets with more than a million points it yields prohibitive computation time. This is why, tree structures are usually used to partition the space into cells. In a nutshell, this partition permits to avoid traversing cells that are too far away from the query point ( [Sam]).

Here an octree was used: it simply divides the bounding box of side L into eight cells of size L/2. The splitting point is of course the center of the parent cell. All points are included in leaf cells, i.e. cells with no child that have depth d, where d is the octree depth specified by the user. This way, only nodes actually containing points are created.

A. Iterating over the octree

Iterating over the octree means from a given point position find a set of neighbors.

In this work radius fixed neighborhoods are used meaning that given a point p and a point set M S , the set of p's neighbors is the subset of M S such that:

The iterating method is basically the one from [FP]. Given a point p, a radius r and a set M S it proceeds as follows:

• In case it is not given, the leaf cell C containing p is found.

• Going upwards in the tree, it finds the cell C 0 containing C and whose length L is such that L 4 < r ≤ L

2

• Then the set of neighbors is included in C 0 and C 0 's neighboring cells at the same depth.

• For each of those cells:

-Look at all non-empty children nodes -If the hypercube maximum distance to p is below r then all points included in its descendants are added to the set of neighbors -If the hypercube maximum distance to p is above r and the minimum distance is below r then, if the node is not a leaf, the node's children are explored recursively. If it is a leaf then all the point distances are computed and compared to r.

-If the hypercube minimum distance to p is above r the cell is not explored any further.

The whole neighboring cell search is made faster using Location Codes (see [FP] for details).

tel-00557589, version 1 -19 Jan 2011

A. A single octree to deal with an evolving point set

Since the pointset evolves with filtering iterations, an evolving structure has to be built. This is done by simply considering that each cell contains more than 1 set of points. Since the initial pointset must be remembered, it is necessary to store  sets per leaf: the set of initial points, the set of points filtered at iteration N and the set of points filtered at step N + 1. Then when performing neighbor-search, the appropriate set should be chosen.

• When filtering set 0, the results are stored in set 1

• When filtering set 1, neighbor-search is done in set 1, the results are stored in set 2 and set 1 is emptied

• When filtering set 2, neighbor-search is done in set 2, the results are stored in set 1 and set 2 is emptied Each point links to the point of set 0 from which it originated.

A. Implementation

All the algorithms presented in this thesis were implemented using C++ and the Standard Template Library (STL). It uses the Template Numerical Toolkit and JAMA-C++ libraries  for Principal Component Analysis and system solving.

Reading/Writing images for color cloud filtering for example was done using the CImg library  .

Finally, reading and writing mesh in Stanford PLY format used rply  . For perenity reasons, files are always written in ascii PLY and not binary PLY, in order to be able to access the information on any computer in the next years, at the cost of bigger files.

Throughout this thesis renderings have been computed using the POV-RAY program  . No texture was used (the triangles were colored in white) and a simulated diffuse light was set.