R Quiroz

T Alves

B Poussot

J.-M Laheurte

Combined RFID tag antenna for recipients containing liquids

In UHF RFID (radio frequency identification), passive tag antennas can be fixed on boxes or recipients where the nature and contents might vary. These variations strongly affect the antenna performance. A combined antenna is proposed to ensure an effective read-range for a plastic recipient containing water or not. Two separate antennas are first designed for the filled and unfilled cases, respectively, then combined for a correct working in both configurations.

Introduction:

The close environment of UHF RFID tags, i.e. the medium on which the inlay is fixed (plastic, metal, cardboard) and its contents (liquids, highly dielectric, metallised,…), degrades the radiation pattern, the antenna matching to the integrated circuit and the overall tag efficiency. This, added to the obstructions and multipath effects in indoor environments, can dramatically reduce the read-range between the reader and the tag.

In this Letter, the objective is to design a tag antenna attached to a plastic recipient that may be empty or filled with a liquid. The originality of the structure relies on a small-loop based module which is used to excite two dipole antennas through inductive coupling. Each antenna is designed to work either for the filled or unfilled case at 868 MHz (UHF band).

Module description:

The MuTRAK640 module manufactured by Tagsys [1] is essentially the series connection of a small loop antenna with a UHF RFID chip, the Impinj Monza 4 [2]. The measured chip read sensitivity and input impedance are -14.7 dBm and 1100Ω// 2.11pF, respectively [START_REF] Saba | A Resistant textile tag antenna for RFID UHF frequency band[END_REF]. Let Z c = R c + jX c be the series equivalent impedance of the chip load (Z c = 6.8Ω -j86 Ω at 868 MHz). The circuit is encapsulated in a rectangular housing made of FR4 epoxy. A key point is that the small loop dimensions yield a radiation resistance smaller than 1Ω. Therefore, the loop only allows short-range reading distances (few centimetres) restricted to the reactive near-field region.

A useful feature of the module is that it can be used as a primary excitation for larger tag antennas. The idea is that a small device basically developed for short reading distances because of its low radiation resistance also performs well at distances up to 10 m when coupled to a dipole-like antenna. Typically, the magnetic field normal to the loop surface turns around the dipole located at a close proximity of the module and generates a current in the dipole. The dipole boosts the module radiation by increasing its radiation efficiency.

Inductive coupling and antenna matching: The excitation of a dipole by an inductively-coupled loop yields the following impedance seen from the chip terminals [START_REF] Son | Design of RFID tag antennas using an inductively coupled feed[END_REF]:

Z antenna = Z loop + (Mv) 2 Z dipole (1)
The above expression includes the mutual coupling factor M caused by the proximity between the dipole and the small loop, and the dipole impedance Z dipole besides the loop impedance Z loop . From (1), it can be observed that the series resonance of a half-wavelength dipole is transformed into a parallel resonance at the loop terminals. Equation (2) gives the power wave reflection coefficient of the antenna normalised to the antenna resistance [5]:

G = Z c -Z * antenna Z c + Z antenna (2)
Antenna design: The maximum power transfer between the antenna and the chip will occur for Z antenna = Z c *. Apart from the compensation of the chip capacitance by the antenna inductance, it is also necessary to keep a low antenna resistance. A dipole with spiral inductors at its ends is coupled to the encapsulated circuit MuTRAK as shown in Fig. 1. Spiral dipoles are compact and more efficient than zig-zag or meandered dipoles as large gaps between the parallel segments are introduced. The dipole was made of a thin copper wire which is fixed on a piece of paper to keep the antenna shape and rigidity. The dimensions of the dipole are given in Fig. 1. The distance between the module edge and segment L1 is 2mm. The loop-coupled spiral dipole antenna is designed to resonate at 868 MHz. All simulations are performed with the 4NEC2 code based on the method of moments. The simulated impedance is Z antenna = (5.1 + j62) Ω for the required frequency (Fig. 2). The estimated gain is 0dBi. Coupling the dipole to the module introduces a parallel resonant at the chip terminals at 868 MHz as observed in Fig. 2. As a result, the antenna impedance shows a higher radiation resistance and a flatter reactance around 868 GHz compared to the loop impedance. However, the reactance mean value is essentially the loop reactance which is too low at 868 MHz to exactly compensate for the -86Ω reactance of the chip. Therefore, even though the comparable antenna and chip resistances result in a Γ minimum, the total reactance (approximately -24Ω) is large compared to the total resistance (approximately 12Ω) which leads to Γ = 0.89. It means that the read-range which would be obtained for a perfect match is divided by 1/ 1 -|G| 2 ≃ 2.2. This mismatch is inherent to the module resonance centred at 915 MHz, not to the design procedure. The following tag resonance is due to the cancellation of the total reactance and occurs above 1000 MHz, i.e. much higher than the module resonance as the dipole coupling reduces the reactance of the isolated loop. Measurements of initial tag: Measurements of the dipole resonant frequency are performed with the King's shielded loop [START_REF] King | The loop antenna for transmission and reception[END_REF] used as a near-field sensor in the vicinity of the dipole. This particular probe avoids external sheet currents on the measurement cable. The dipole resonance is observed on a vector network analyser at the King's loop terminals. The dipole length is adjusted to resonate at 868 MHz, yielding the dimensions given in Fig. 1. Then, the MuTRAK module is implemented and read-range measurements are performed with the help of the Voyantic Tagformance Measurement System [START_REF] O'connor | Voyantic helps companies put RFID tags to the test[END_REF]. The tag read-range measured in the 800-900 MHz band (Fig. 3) shows a maximum of 3.8 m at 868 MHz. This is in agreement with the estimated read-range using the modified Friis formulation given in [5]. Measurements with empty and filled plastic recipient: The spiral antenna described previously is placed on an empty plasticpolypropylene (ε r = 4) recipient which has a size of 23 × 15.5 × 14.5 cm. Due to the plastic wall, a 60 MHz negative shift of the resonant frequency is measured. To adjust the resonance frequency the wire length is reduced, by gradually and symmetrically shortening the wire until the correct resonant frequency is obtained. A total length of 11.5 mm is finally removed from the dimensions given in Fig. 1 (L7 = 0, L6′ = 5.5mm). The maximum read-range drops to 3.7 m.

ELECTRONICS LETTERS 14th

Once the recipient is filled with fresh water (ε r = 80, σ = 0.1 S/m), tag detection is not possible even at a few centimetres from the reader. Using the King's shielded loop, a 550 MHz negative frequency shift of the dipole resonance is first determined. An additional 19 mm length reduction is necessary to resonate at 868 MHz (L7, L6 and L5 equal to zero and L4′ = 13mm). In the presence of water, the maximum read-range is 34 cm.

Combined antenna: The final tag combines both proposed antennas and the MuTRAK module as shown in Fig. 4. As the dipole resonant lengths are very different, no direct coupling between the dipoles is observed and each dipole behaves as if it were alone. Dipole 1 and dipole 2 resonate with the empty and filled recipient, respectively. Without water, the maximum read-range remains around 3.7 m. In presence of water, the read-range is reduced to 31 cm, i.e. 3cm shorter than for the single dipole. In any case, this is a strong improvement compared to the single antenna designed for the empty recipient where no detection is observed for any distance.

Fig. 1

 1 Fig. 1 Spiral dipole coupled to MuTRAK module (loop + chip) Initial dimensions spiral dipole working in free-space: L1 = 40mm, L2 = 20mm, L3 = 13.5mm, L4 = 16.5mm, L5 = 10mm, L6 = 13.5mm, L7 = 3.5mm, wire radius = 0.15mm. Total length: 181 mm

Fig. 2

 2 Fig. 2 Impedance of loop-coupled spiral dipole

Fig. 3

 3 Fig. 3 Read-range measurements with Voyantic setup

Fig. 4

 4 Fig. 4 Combined antenna structure


T Alves

B Poussot

J.-M Laheurte

Non-intrusive estimation of loss resistance of planar antenna over human body

A new method is presented to estimate the loss resistance of antennas located slightly above a lossy medium. This technique is more specifically developed for the context of body area networks. This method is based on the estimation of the antenna resistance above the medium and above a perfect electric conductor.

Introduction:

The 'Wheeler cap method' [START_REF] Wheeler | The radiansphere around a small antenna[END_REF] is the most commonly used method to estimate the efficiency of small antennas. This consists in enclosing the antenna with a metallic box that is large enough to not disturb the reactive near field and cancel the antenna radiation (typically cap dimensions should be greater than λ/2π). The antenna efficiency is then calculated from the knowledge of the antenna impedance with and without the Wheeler cap.

However, in the presence of any lossy medium, the induced currents into it can be seen as an extension of the antenna body [START_REF] Monteath | The effect of the ground constants, and of an earth system, on the performance of a vertical medium-wave aerial[END_REF]. As a result, the total reactive field is contained in a nearly λ/2 diameter sphere [START_REF] Trainotti | Short low-and medium-frequency antenna performance[END_REF] which, in turn, involves a bigger Wheeler cap.

In the context of body area networks (BANs), antennas are in the vicinity of the body, a lossy medium whose composition [START_REF] Alves | Conception d'Antennes pour le Réseau BAN et Modélisation du Canal de Propagation[END_REF][START_REF] Christ | The dependence of electromagnetic far-field absorption on body tissue composition in the frequency range from 300 MHz to 6 GHz[END_REF] may largely vary. In [START_REF] Salim | Efficiency measurement of antennas for on-body communications[END_REF], the Wheeler cap method was applied to BAN antennas with small ground planes partially including the phantom inside the cap. Other authors have introduced the 'body-worn efficiency' [START_REF] Kellomäki | Effects of the human body on single-layer wearable antennas[END_REF] which is the ratio of the in-body 'radiated power' to the free space radiated power.

The methods described above either require a field measurement inside a body phantom or an unpractical Wheeler cap overlapping body tissues. This Letter proposes a simple method based on the compensation theorem to estimate the loss resistance of BAN antennas. From the simulated or measured values of the input resistance of the antenna, first above the body tissues, then above a metallic ground plane, it is possible to extract the loss resistance above the body tissues. Knowing the loss resistance R L and the total input resistance R in , it is further possible to deduce the antenna efficiency η = 1 -R L / R in in the presence of the body.

Compensation theorem: The compensation theorem applied in electromagnetics was first demonstrated by Monteath [START_REF] Monteath | Application of the compensation theorem to certain radiation and propagation problems[END_REF]. It gives the impedance variation ΔZ between the impedance Z of the antenna located above a lossy medium and the impedance Z PEC of the same antenna located above a perfect electric conductor (PEC) [START_REF] Wait | On the input impedance of a hertzian dipole over a flat surface[END_REF][START_REF] Wait | Impedance characteristics of electric dipoles over a conducting half-space[END_REF]

DZ = Z -Z PEC = 1 I 2 0 S H 2 // Z s dS (1)
† Z s : surface impedance modelling of the lossy medium † Z: impedance of the antenna over the lossy medium † Z PEC : impedance of the antenna above a PEC † ΔZ: impedance increment † H // : tangential magnetic field over a PEC † I 0 : root-mean-square current value at the antenna feeding point † S: boundary surface exhibiting a surface impedance Z s .

Power considerations: Comparing the power delivered to the antenna over a PEC and over a lossy medium, we exploit (1) to show that an antenna over a lossy medium consumes an additional power that is

DP = <{DZ}I 2 0 = < S H 2 // Z s dS (2a)
This leads to an increment of power given by

DP = S H // 2 Z s | |cos w s + 2w h dS (2b)
where w s is the phase of the surface impedance and w h is the phase of the tangential magnetic field. On the other hand, the power P L absorbed by any lossy medium is given by the well-known formula

P L = S H // 2 Z s | |cos w s dS (3)
From (2b) and (3), it is clear that P L and ΔP are equal when w h is small. This is verified when the antenna is close to the medium [START_REF] Wait | On the input impedance of a hertzian dipole over a flat surface[END_REF][START_REF] Wait | Impedance characteristics of electric dipoles over a conducting half-space[END_REF]. Therefore, we conclude that when an antenna is close to a lossy medium, the additional consumed power is mainly due to the medium losses. Under this condition, the loss resistance R L of the antenna is given by

R L ≃ <{DZ} (4)
The method presented in this Letter is based on this assumption.

Estimation of loss resistance:

To control the validity of (4), we compare the R L values determined by two methods: a. R L determined from (4) and 4NEC2 simulations [START_REF]4NEC2 is a free interface developed by Arie Voors to use the NEC2 engine[END_REF]: † Computation of the antenna impedance Z over a lossy medium † Computation of the antenna impedance Z PEC over a PEC † Use of (4):

<{DZ} = < Z -Z PEC { }=R L .
It is important to note that this method can also be applied experimentally.

b. R L determined from the following relation using the HFSS calculations of the antenna impedance Z and the efficiency η in the presence of a lossy medium:

R L = < Z { }(1 -h) (5)
The comparison between these methods is realised at 2. agreement is observed between the two methods. However, it must be stressed that the relative differences are 2 and 8% for h = λ 0 /40, and 17 and 25% for h = λ 0 /10 with the IIFA and the dipole, respectively. This increasing divergence for larger heights flows from the approximations in (4) that are not valid for large h/λ 0 , i.e. the phase of the tangential magnetic field is no longer negligible on the body surface in the vicinity of the antenna. We conclude that the R L estimation based on the compensation theorem is valid in the body context as long as the height above the interface is lower that λ 0 /10, approximately (1.25 cm at 2.4 GHz), which is sufficient in most BAN scenarios. Measurements are realised with the IIFA, first on the torso phantom, then on a large copper plate. Special care is taken to avoid cable radiations from the flow of external current with baluns made of a pair of λ/4 stubs. Tunable plastic spacers are used to fix the antenna heights above the phantom surface. R L is then deduced from the procedure described above and compared with the value obtained with the same procedure using the HFSS impedance values. In Fig. 3, the R L curves follow the same trend. Discrepancies can be attributed to the difficulty in measuring the correct impedance in small antennas where the antenna body extends partially to the cable and the balun structure. In addition, the presence of the antenna connector and cable in the vicinity of the interface perturbs the measurements. E-mail: laheurte@univ-mlv.fr

4 Fig. 1 Fig. 2

 412 Fig. 1 Details of IIFA Dimensions: w = 1 mm, s = 2 mm, l = 10.5 mm, d = 3.5 mm FR4 substrate (ɛ r = 4.4, tan δ = 0.02); thickness = 1.6 mm

Fig. 3

 3 Fig. 3 Estimation of IIFA loss resistance based on compensation theoremComparison of experimental and simulated resultsConclusion:A new technique has been described to estimate the loss resistance of an antenna close to a lossy interface. The method is


Dobrina Boltcheva

Bruno Lévy

Simple and Scalable Surface Reconstruction

Keywords: Descriptors (according to ACM CCS): I.3.5 [Computer Graphics]: Computational Geometry and Object Modeling-Boundary representations & Curve, surface, solid, and object representations

Figure 1: Our method applied to a pointset with 33 million points. The total computation took 360 seconds (49 sec. for 2 iterations of smoothing and normal estimation, 78 sec. reconstruction, 61 sec. surface extraction, 114 sec. post-processing and 58 sec. loading and saving)

Introduction

Over the past two decades there has been a significant amount of effort dedicated to the problem of surface reconstruction. The problem of surface reconstruction may be formulated as follows: given a sampling of points from a surface, recover the original surface from which those points were sampled. This general problem is motivated by numerous applications in reverse-engineering, prototyping, visualisation, or computer vision since a growing variety of scanning devices nowadays provides measurements of objects in the form of point sets.

In the surface reconstruction problem we are given only a finite sampling P ⊂ R 3 of an unknown surface S. The goal is to build a model of the surface S from the sampling P which is referred to as a reconstruction of S from P. It is generally represented as a triangulated surface mesh that can be directly used by computer programs from further processing. The reconstruction should match the original surface in terms of geometric and topological properties. In general, surface reconstruction is highly under-constrained and there are multiple triangulated surfaces that might fulfil these criteria. The difficulty of meeting geometric and topological criteria depends on properties of the sampling and on properties of the sampled surface. In particular, sparsity, redundancy, occlusion, noise, non-smoothness and boundaries make surface reconstruction a challenging problem.

Surface mesh reconstruction methods can be divided into methods that approximate the points by an implicit function and methods that connect the points to form a surface mesh. The first approach is widely spread because it generates smooth and closed meshes, usually by extracting the zero level of some potential field. The second type of approach is useful whenever it is desired to keep the input points, for instance when they represent measurements. In geological applications, the points are usually obtained by complex drilling simulations and they are expected to be contained in the final mesh for further processing. In this paper, we are interested in the second kind of approach, which aims at building a mesh by connecting the input points.

We present a new reconstruction algorithm that generates a surface triangulation from an input point set. Our algorithm is designed to process a clean and accurately registered point set, and does not attempt to average out noise, outliers or residual registration errors. The algorithm has several interesting properties: it is very simple to implement and it is time and memory efficient.It is "embarrassing parallel" and outperforms the concurrent algorithms (by two order of magnitude in some cases). Our algorithm uses only one data structure (a kd-tree), and accommodates evenly sampled irregular points sets, as demonstrated in the results. It handles large pointsets -up to 36 million points in minutes on an off-the-shelf PC.

Related work

There are two classes of approaches which differ as to whether they approximate or interpolate the input points.

The approximation methods aim to fit the best surface approximating the sampling and they usually ignore the sampling noise in the input data. A recent state-of-the-art and a benchmark of this family of algorithms can be found in [BLN * 13, BTS * 14]. These methods usually fit the points using the zero set of an implicit function, such as a signed distance function [HDD * 92, CL96, BC02], a sum of radial bases [CBC * 01], piecewise polynomial functions [ABCO * 03, OBA05, NOS09] or an indicator function [KBH06, KH13, MPS08, ACSTD07]. In particular, the Poisson surface reconstruction method solves for an approximate indicator function of the inferred solid, whose gradient best matches the input normals. The output scalar function, represented in an adaptive octree, is then iso-contoured using an adaptive marching cubes [START_REF] Lorensen | Marching cubes: A high resolution 3d surface construction algorithm[END_REF]. Recently, a screened Poisson algorithm [KH13] has been introduced which resolves the over-smoothing problem and offers the ability to reconstruct meshes with boundaries.

Most of these methods estimate a signed-distance function. If the input points are unoriented, an important step is to correctly infer the sign of the resulting distance field [ACSTD07, MdGD * More precisely, our algorithm belongs to the category of the Restricted Delaunay based surface reconstruction methods. The main idea behind these methods is to filter out a subset of the Delaunay triangulation Del(P) of the sampling P, by restricting it to some subset of R 3 (thickening of the point set) which is a good approximation of the unknown surface S and can be efficiently computed from P.

This idea is motivated by the fact that, if the input meets prescribed sampling assumptions on its quality (ε-sampling), it has been proven that the output mesh is homeomorphic to S (see the theorems in [ES94, AB98]). The general ε-sampling framework proposed by Amenta Bern [AB98] defines a real-valued function f , called the local feature size on a smooth surface S (as the distance to the medial axis of S) and, essentially, requires that the sample set P leaves no ball centered at a surface point p and of radius ε × f (p) unsampled.

Alpha shapes are one of the first attempt to use this idea for surface reconstruction, [START_REF] Edelsbrunner | Three-dimensional alpha shapes[END_REF]. These methods define the dilation of the pointset as the union of balls of radius α: P = ∪Bα(p i ). Then they build the α-complex which is the restriction of the Delaunay triangulation of P to the union of α-balls : Del(P)| ∪ Bα(p i ).

The Ball Pivoting algorithms [BMR * 99, DADSG11] are optimized variants of Alpha Shapes where the computation of the 3D Delaunay triangulation is avoided. They generate the surface incrementally, starting from a single triangle and adding new triangles to the surface by pivoting an α-ball. Since all the computations are local, these algorithms exhibit linear complexity and are efficiently parallerized to handle large data sets [START_REF] Digne | An Analysis and Implementation of a Parallel Ball Pivoting Algorithm[END_REF]. The method introduced by Digne et al. [DMSL11] improves the noise-resilience of the BPA by applying it to a scale-space version of the point set.

The Cocone algorithms, fist designed by Amenta and Bern [ACDL00], build on the same idea. Here, the thickening of the pointset is defined as the union of co-cones placed at each point (double cones with opening angle of π/8 with respect to the estimated normal). For the Cocone algorithm the theoretical guarantees come with a sampling density at least 0.18 × l f s(x) and for θ = π/8. Intuitively, the Cocone angle θ absorbs the deviation between the estimated normal vector and the true normal of the surface.

One of the main limitations of the Cocone algorithm is its computation cost. For large datasets, the data structures employed exceed the system memory and processing speed becomes very slow. Even if efficient and robust implementations of the 3D Delaunay triangulation are used [CGAL, GEOGRAM], these algorithms do not scale up well and cannot accommodate large data sets comprising dozens of million of points within a reasonable space and time. Moreover, even on well sampled smooth surfaces, the worstcase complexity of the 3D Delauany triangulation can be quadric O(n 2 log n), where n = P is the size of the sampling [Eri01].

Local Cocone variants were developed to mitigate this problem. These algorithms are local in sense that the restricted Voronoi cells of each point are computed using a small neighbourhood around the point. Our method belongs to this category and in the following we detail the state of the art of these algorithms.

Funke and Ramos [FR02] demonstrated that the Cocone algorithm can be modified so that no global computation of the 3D Delaunay triangulation of the entire pointset is required. They introduced an additional local uniformity assumption on the sampling which requires a minimum separation between the sample points (there is some δ ∈ {0, 1} such that any p, q ∈ S must be separated y a distance of at least δεl f s(p)). This involves decimating the input pointset so that the working point set is a locally uniform subset of the original input. Within this theoretical framework they showed that the cocone triangles can be computed locally. Indeed, the resulting restricted cocone of each point p is a local object in a sense that it is completely determined by nearby samples, when the sampling conditions are met. When the algorithm expands the calculation of the 3D Voronoi cell for each point p, it stops when the distance between p and the new samples exceeds 3 × dmax, where the maximum distance of a point on a Voronoi edge of Vp within Cp. Thus, the algorithm makes sure that the cocone region of each 3D Voronoi cell (with respect to the estimated normal) is computed exactly since further nearest neighbours are added until it is sure that the cocone region is not affected any more. The complexity of the resulting theoretical algorithm is almost linear O(n log n), but it has not been implemented. Amenta and Kil [START_REF] Kil | Gpu-assisted surface reconstruction on locally-uniform samples[END_REF] employed similar techniques to design a very simple and efficient reconstruction algorithm which processes the sample points in parallel on the GPU. However, all these algorithms require a noise free locally uniform ε-sampling which is not reasonable to expect from raw pointsets and the problem of building such sampling from an arbitrary input pointset is still open.

Alternatively, Dey et al. [DDW11] developed an octree-based version of the Cocone algorithm such that the 3D Delaunay triangulation is only computed on small clusters of the total point set. They relaxed the local uniformity criterion and adopted the assumption that the local sampling density is bounded by a constant, [ABL03]. Under this assumption, the Localized Cocone maintains the theoretical guarantees of the original Cocone algorithm.

The Localized Cocone has been implemented using the CGAL library, to facilitate the computation of the 3D Delaunay triangulation on each octree leaf. The reported results show that it reconstructs inputs of 5 million points in more than 1 hour and uses more than 2000MB of memory. Our algorithm described below takes less than a minute on comparable point sets. In addition, in the localized cocone algorithm, when the input pointset does not meet the stringent sampling requirements, the manifold extraction step is not particularly robust. If the sampling density is insufficient to ensure a full umbrella of cocone triangles around a vertex, then in the subsequent trimming stage of the algorithm, all cocone triangles are removed. In practice, it is necessary to introduce heuristics to prevent catastrophic failure to accommodate real sample sets. But in this case, there is no longer a guarantee that the output will be topologically correct.

The SSSR algorithm

Our algorithm is strongly related to the theoretical algorithm of Funke and Ramos [FR02]. In particular, it uses similar ideas as in the Localized Cocone algorithm [DDW11] and the Vorpaline algorithm [LB12] but it differs from previous Cocone algorithms in two main respects:

• Our algorithm computes the Delaunay triangulation of the pointset restricted to a set of disks. Instead of computing the Delaunay triangulation then restricting it, we directly compute the restricted Voronoi cells (intersections between the Voronoi cells and the disks), and deduce the restricted Delaunay triangles from the combinatorial information of the restricted Voronoi vertices (see Section 3.4, Figure 3, Figure 5 below), this results in an algorithm that is simple and "embarassingly parallel". • We also present a new manifold surface extraction method which avoids the trimming step (the most critical because it can end up with an empty mesh). Starting from a subset of candidate triangles which forms a clean and orientable manifold with boundaries, we add the remaining triangles, one by one, only if they do not break the topological properties of the initial mesh. The algorithm is guaranteed to terminate with a non-empty manifold mesh.

The input is an unorganized 3D point set, possibly with normal attributes (unoriented or oriented). The output of the algorithm is a triangular mesh connecting the input points which approximates the surface. Specifically, the output triangulation is a compact 2pseudo-manifold (see Section 3.7), it is orientable (no Moebius The algorithm uses two pre-processing steps which allows us to accommodate raw data sets. First, it uses a smoothing algorithm to reduce the noise in the input data. Second, if input data does not come with normals, it uses a normal estimation algorithm. Note, however, that this algorithm does not need to consistently orient the normals.

The general workflow of our algorithm in summarised in Figure 2. We give details of each of these steps in the following sections.

Kd-tree construction

During this initial step the data is loaded and prepared for the remaining operations. The vertices P are organized in a geometric search data structure, such that for any p i we can efficiently compute the list of nearest neighbours p j sorted by increasing distance to p i , [Sam05].

To do this step, one can use the Approximate Nearest Neighbours (ANN) library [MA97] which with ε set to 0 computes the exact nearest neighbours for each point. We have developed our own version of the same algorithm that only uses contiguous arrays in memory. To save one integer per vertex, we compute a balanced binary tree. As a consequence, the links are completely implicit in the tree.

Smoothing

Since our algorithm interpolates the input points, it requires a clean and accurately registered point set. Thus, if the input set is noisy such as raw scanning data, it has to be smoothed before the reconstruction.

Our current implementation uses projection onto the best approximating plane of the k-nearest neighbours, [HDD * 92]. Specifically, for a point p and its k nearest neighbours {p i } k i=1 , we find the fitting plane n T x = c for p by minimizing the error term e(n, c) = ∑ k i=1 (n T p ic) 2 under the constraint n T n = 1. Then, the point p is projected onto this plane. We compute the best fit plane for each point in parallel. Once all the points are projected, we update the kd-tree.

Parameters: This step uses two parameters : the number of nearest neighbours nb_neighbors (by default 30) and the number of smoothing iterations smooth_iter (by default 0). We have used the default values for the small data sets shown in the results, while, for the large ones, we have done two iterations

Normal estimation

If the input points come without normals, our algorithm estimates them locally with the best fitting plane following exactly the same formulation as in the smoothing step, [HDD * 92].

Notice that the algorithm dose not need to orient the normals since we only need an approximation of the (unoriented) normal direction at each point which is used to create the tangent disk, in the next step of the algorithm.

Parameters: This step uses one parameter: the number of nearest neighbours nb_neighbors set by default to 30. Given a pointset {p i } n i=1 with estimated normals {n i } n i=1 and a radius r, our algorithm computes the intersection between the Voronoi cells of the points Vor(p i ) and the disks D r i of radius r centered on the points and orthogonal to the normals. Figure 4: The Voronoi cell Vor(p i ) is the result of clipping the entire space by all the bisectors (continuous lines) of the segments [p i , p j ] (dashed lines). Among the bisectors, some are contributing (black) and some are non-contributing (gray).

Candidate triangles

During this step, we fist build a disk D r i at each point p i ∈ P orthogonal to n i and with user-given radius r . In our implementation, we use a polygon (with 10 vertices) that approximates the disk. In contrast with previous algorithms, we explicitly build the restriction of the Voronoi cell Vor(p i ) to the tangent disk D r i at p, Vor(p i ) ∩ D r i , as shown on Fig. 3 and Fig. 5. Thus we collect a set of restricted Voronoi vertices.

The set T of candidate triangles is then defined by the triangles which are dual to the restricted Voronoi vertices (more on this below).

The basic operation in this step of the algorithm consists of clipping the disk D r i with the Voronoi cell Vor(p i ) of each point p i . Here we use ideas that are similar to those presented in the Localized Cocone algorithm [DDW11] which allows us to build locally the candidates triangles by using only the k-nearest neighbours. More precisely, we use the clipping algorithm and the local characterization of the restricted Voronoi cells introduced in [LB12] which we recall in the following in order to make the explanation self-contained.

Clipping a disk

Let us recall that a Voronoi cell is a convex polytope, that can be defined as the intersection of halfspaces : Vor(p i ) = j =i Π + (i, j), where Π + (i, j) denotes the half-space bounded by the bisector of (p i , p j ) that contains p i . Note that the bisectors between p i and all the other points are involved in the definition above, whereas only a small subset of them corresponds to actual Voronoi edges (see Figure 4). As a consequence, we can classify the bisectors Π + (i, j) into two sets : non-contributing if Vor(p i ) ⊂ Π + (i, j) and contributing otherwise. In other words, clipping a Voronoi cell by a non-contributing bisector does not change the result. Therefore, the Voronoi cell can be computed by intersecting a superset of the contributing bisectors.

Let p j1 , p j2 , . . . p jn-1 denote the vertices sorted by increasing distance from p i . Let V k denote the intersection of the k first halfspaces i (in blue) by the bisectors of the neighboring points p ji . The so-defined restricted Voronoi vertex q i j1 j2 and the associated restricted Delaunay triangle (p i , p j1 , p j2 ) are highlighted.

and let R k denote its p i -centered radius :

V k (p i ) = k l=1 Π + (i, j l ) R k = max{d(p i , x)|x ∈ V k (p i )}.
We can then make the following simple observation: For all j such that d(p i , p j ) > 2R k , the bisector

Π + (i, j) is non-contributing, i.e. V k (p i ) ⊂ Π + (i, j). Consider p ∈ V k (p i ) and p j such that d(p i , p j ) > 2R k . By defi- nition of R k , d(p, p i ) < R k . We have d(p i , p) + d(p, p j ) ≥ (p i , p j ) (triangular inequality) and d(p i , p j ) > 2R k , therefore d(p, p j ) > R k > d(p, p i ) and p ∈ Π + (i, j)
As a direct consequence of the previous observation:

d(p i , p j k+1 ) > 2R k ⇒ V k = Vor(p i ).
We call radius of security the first value of R k that satisfies this condition. Note that this observation does not have a practical value in the case of the (unrestricted) Voronoi diagram, since some cells are unbounded and have infinite R k therefore, for an infinite cell, all the bisectors are considered, leading to prohibitive computation time. However, the observation can be clearly used to compute the restricted Voronoi cells Vor(p i ) ∩ D r i , since they are finite (they are contained by D r i ). We use Algorithm 1 to build the restricted Voronoi cell for each point, which is the intersection between the tangent disk and the Voronoi cell of the point.

Note that each restricted Voronoi vertex q is the intersection between the disk D r i and the bisectors Π i, j1 and Π i, j2 . When we build the clipped disk, we memorize the indices j 1 and j 2 which generate every restricted Voronoi vertex. This allows us to easily output the corresponding dual triangle i, j 1 , j 2 , see Fig. 5.

Set of candidate triangles

Once we have computed the restricted Voronoi cell D r i ∩Vor(p i ) for each point (p i ), we are able to find all the triangles i, j 1 , j 2 using the stored indices. We call T the set of candidate triangles, that are dual Algorithm 1 Computes a Restricted Voronoi Cell, i.e. the intersection between a disk and the Voronoi cell of a point. Data: the index i of the point x i , the disk Dr centred on it, the number of neighbours n and the set of points P Result: Restricted Voronoi Facet at p i :

Vor(p i ) ∩ D r i V ← D r i R k ← max{d(p, p i )|p ∈ V } k ← 1 while d(x i , x j k ) < 2R k f and k < n do V ← V ∩ Π + (i, j k ) R k ← max{d(p, p i )|p ∈ V } k ← k + 1 end
to every restricted vertex. From this set T of candidate triangles, that may contain a non-manifold configuration, we extract a clean surface, in the next steps of the algorithm. Let us further analyse the structure of T :

Figure 6: A: the set of triangles seen from three restricted Voronoi cells (T 3 ), with some holes (yellow). B: the set of triangles seen from one or two Voronoi cells (T 1,2 , in red) fills some holes, but have non-manifold configurations (e.g., the "sliver" formed by the four red triangles in the right). C: result after carefully inserting the T 1,2 triangles one by one in a way that preserve manifold topology. D: the remaining holes are suppressed by a post-processing step.

The set T of candidate triangles is composed of the restricted Delaunay triangles T 3 plus some other triangles, as explained below.

The set T 3 of restricted Delaunay triangles is defined as:

T 3 = {(i, j, k) | (Vor(p i )∩D r i )∩(Vor(p j )∩D r j )∩(Vor(p k )∩D r k ) = ∅}. (1)
In other words, they correspond to triples of restricted Voronoi cells (i, j, k) that are mutually in contact (see Figures 5 and3). These triangles can be easily found, by generating all the index triples that correspond to the restricted Voronoi vertices, sorting them (using for instance the lexicographic order) and keeping the triples that appear three times in the sorted sequence. The socomputed T 3 triangles are displayed in Figure 6-A. Note that in the sorted sequence, there are also triples that appear once or twice. The corresponding set of triangles is referred to as T 1,2 and is defined by:

T 12 = {(i, j, k) | D r i ∩Vor( j) ∩Vor(k) = ∅} -T 3 . (2)
In other words, a triangle (i, j, k) is in T 12 whenever the restricted Voronoi cell of p i "sees" p j and p k but not both conversely. A 2d example is shown in Figure 7(a): in the Voronoi cell of p j , D r j (symbolized as a thick line) touches a side that is adjacent to a cell Vor(p k ) that is different from Vor(p i ). Other examples are shown in Figure 7(b), that correspond to nearly co-spherical configurations. The so-defined T 12 triangles are displayed in red in Figure 6-B. As expected, they are not as "reliable" as the T 3 triangles, in the sense that they are likely to generate non-manifold configurations (see for instance the four interconnected red triangles that form a "sliver" in the top right part of the neck in Figure 6-B). However, one can observe that they contain interesting information, and could be used to fill the holes in the T 3 triangles (Figure 6-A). Thus, the idea is to start from the T 3 triangles, then carefully insert the triangles from the T 12 set, one triangle at a time, making sure that the reconstructed mesh remains manifold with the additional triangle inserted. If it is not the case, the triangle is rejected. The algorithm that implements this idea is detailed in the next section.

Parameters: This step uses one parameter : the disk radius (in % of the bounding box diagonal). Note that this is not a "real" parameter as it is the case for any Ball Pivoting algorithm, for example. Indeed, when the radius r is sufficiently large (i.e. larger than the width of the Voronoi cells), the disks D r touch all the sides of their Voronoi cells (like on the bottom-right of Figure 3), and increasing the radius no longer changes the result. Thus, the radius can always be set to some large value. In our experimental results, we used r=5% of the bounding box diagonal and we decrease the size in order to accelerate the clipping for the very large data sets (r=0.5% bbox diagonal when number of vertices > 10M).

Manifold mesh extraction

The original "prunning and walking" approach [ACDL00] consists of first deleting all triangles incident to sharps edges and then extracting the outer boundary of the set of remaining triangles by a depth-first walk along each of its connected components. As mentioned earlier, this algorithm is not particularly robust and uses many heuristics for the practical implementation and may end up with an empty mesh.

Here we take a different strategy which allows us to avoid this problem. We first build a clean orientable 2-manifold mesh with a subset of the candidate triangles. Then we fill the holes iteratively, with the remaining triangles, if they do not jeopardize the topological properties of the output.

More precisely, the mesh structure that we build incrementally is a manifold mesh where each edge can only be adjacent to two facets, edges with only one incident triangle are allowed, as well as triangles adjacent only by one vertex (this is a non-manifold vertex whose neighbourhood is not a topological 2-sphere). The output mesh (also known as a 2-pseudo-manifold [Mun84]) is orientable and can have several connected components.

We initialize the output mesh with all T 3 triangles. Then we remove the triangles which exhibit non-manifold edges and nonmanifold vertices by excess. Recall that a non-manifold edge is incident to more than 2 triangles, while a 'by-excess' non-manifold vertex has a closed loop of triangles (clean umbrella) in its neighbours plus additional triangles (see Fig. 8(a)) We also check that there are no Moebius strips within the initial mesh (see Fig. 8(b)). Note, that the initial mesh can however exhibit triangles incident to each other only by a vertex. We then fill the holes of the initial mesh with the remaining triangles from T 12 . The T 12 triangles are tested against the following criteria, from the simplest to the most time consuming. Whenever a criterion is not satisfied, the considered triangle is rejected. Only the triangles that pass all the test are inserted in the final surface.

• Geometric test: This criterion ensures that the normals of the neighbouring triangles agree, which means that they do not make a sharp edge (less than π/4) In the following, we give the overview of the manifold extraction algorithm :

Algorithm 2 Extracting the manifold mesh. Data: A set of candidate triangles T Result: An orientable 2-manifold mesh built from T

(1) Compute the subsets T 3 and T 12 of triangles (2) Initialize the mesh M with the subset T 3 (3) Remove from M : -triangles incident to non-manifold edges (Test II) -triangles which are incident to non-manifold vertices by excess (Test III) -triangles inducing orientability problems (Moebius strip) (Test IV) (4) Tentatively add into M each triangle t from T 12 if it successfully passes all following tests: -Geom. Test : the normals to t and its neighbour should not point to opposite directions -Test I : t is incident to at least two edges of existing triangles, -Test II : t does not create non-manifold edges -Test III : t does not create non-manifold vertex by excess -Test IV : t maintains the global orientability of the mesh The output of this algorithm is an orientable triangular mesh which has no non-manifold edges but can have non-manifold vertices, such as a pinched torus. It possibly has several connected components and boundaries.

Parameters: This step uses one parameter which is the maximum deviation angle between the triangles normals max_N_angle. We always use the default value set to 60 degrees.

Post-processing

At this point, after the manifold extraction test, the obtained mesh may still contain a number of isolated holes, caused mainly by insufficiently sampled sharp features (see e.g. the top of the ears of the horse in Figure 12 C). To eliminate these, we add an extra postprocessing step, that fills some of the holes.

The algorithm takes an user-defined maximum hole size parameter, in terms of number of edges. Then it looks for simply connected holes and fills them using a classical loop-split algorithm [Lie03]. This algorithm splits the hole recursively until it gets holes composed of exactly 3 edges and fills them with the corresponding triangle.

In order to ensure that the holes are simple loops, we first search and remove the configurations where the holes contain "bridges"(see Fig. 9). We call bridge triangles those which are visited more than one time when we walk around the border edges of a hole. Note that this definition dose not capture bridges wider than two triangles but we did not encounter any, in practice. Parameters: This step uses two of the following four parameters which give the size of the holes to be filled and the connected components to be removed. In our experiments, we have always used the default values:

• max_hole_area (=5%) : Fill holes smaller than (in % total area) • max_hole_edges (=500) : Fill holes with a smaller number of edges • min_comp_area (=0.01%) : Remove small components (in % total area) • min_comp_facets (=10) : Remove small components (in facet number)

Experimental results and Discussion

Our implementation is available in C++, in the supplemental material, together with datasets and comparison data for ScaleSpace reconstruction [Dig15], screened Poisson [KH13] and SuperCocone [DDW11]. The implementation is straightforward. A naive implementation fits within 200 lines of code (100 for the Kd-Tree, 50 for the polygon clipping routine and 50 for putting everything together). Our parallel implementation uses 2000 lines (the additional lines are for manifold extraction, book-keeping, handling several options, and for the parallel implementation).

We tested our algorithm on a database of examples available from the world-wide-web [AIM, Sta, Far, EE] and compared the results with the three state-of-the-art methods cited above. Our database of examples is decomposed into two subsets SMALL (up to 1 Million of points) and LARGE (from 1M up to 36M).

Table 1 lists the models and reports the timings for the four algorithms. We used implementations written by the authors of the cited methods, and the same computer for all the tests [Processor: Intel Core i7-4900MQ CPU @ 2.80GHz, RAM:16Gb, Disk: SSD LITEONIT LCS-256M6S 2.5 7mm 256GB, OS:Linux 4.4.0-1-amd64 -Debian distribution]. Our algorithm takes 5 minutes to reconstruct a surface from 36.2 million points. Table 1: Comparison of processing times for the 4 methods. Timings are in seconds and model sizes in million vertices. We used implementations written by the authors of the cited methods, and the same computer for all the tests [Intel Core i7-4900MQ CPU @ 2.80GHz, 16GB]. For the three biggest models, ScaleSpace crashed on all of them, SuperCocone ran out of memory on some models, and Poisson failed on some models (probably due to incoherent normals). KO means that the result has too many holes or is mostly unrelated with the data or the computation time explodes as compared to data of similar size.

The implementation is very simple, with no dependencies, and thus can be easily ported to any architecture. We ported it to an Android phone [HTC M7, Processeur: ARMv7 NeonFPU 4 cores, RAM: 2Gb, Android 4.4 (kitkat)]. The algorithm exploits the 4 cores of the phone. It is approximately 10 times slower than on the PC. It took 1.5 seconds to reconstruct 25K points, 7.5 sec. for 70K points, 14 sec. for 200K points and 55 sec. for 900K points.

The timing breakdown of the different phases of our algorithm, computed on the entire database are summarized in the charts in Figure 10 and size-time curves for all methods and datasets are shown in Figure 11.

We also show some images with example results from our dataset. To be able to analyse and compare the behaviour of the algorithm, we deactivated the post-processing step of our algorithm, so that what we compare with the state of the art in the images is the output of the manifold extraction of our algorithm. Figure 12 shows a horse sampled with 100000 points with a varying sampling density. ScaleSpace has difficulties finding the correct scale in this dataset, and takes more than 3000 seconds. The result has holes when the triangles are larger than the estimated scale. Su-perCocone divides the input points into several subsets using an octree. It sometimes fails to recover triangles that cross octree cells boundaries. This behaviour of ScaleSpace and SuperCocone is also observed on the hand dataset in Figure 13. In this figure, one can observe that Poisson depends on a correct orientation of the normals.

In Figure 14, we show the behaviour of the algorithms for a noisy pointset(original registered pointset data for the Stanford bunny). For this data, Poisson both cancels the noise and recovers fine details. ScaleSpace manages to extract a manifold mesh that exactly follows the details of the data, even when it is very bumpy. Our algorithm outputs a result that is very similar to that of SuperCocone. The same dataset with one iteration of smoothing is shown in Figure 15. All the algorithms give a similar result. Note the aligned missing triangles in SuperCocone at boundaries between octree cells. Figure 16 shows the results obtained from another scanner dataset. Similar behaviour is observed. On this dataset too, ScaleSpace had difficulties estimating the right scale, resulting in some missing triangles (but processing time, 300s, remains reasonable as compared to the irregularly horse dataset). Our algorithm computes the result in 27 seconds.

Figure 17 compares the result on a clean pointset that has bumpy features (the feathers of the eagle). Poisson faithfully reconstructs Figure 16: Experimental results with a scanned pointset of intermediary size (2.6 M vertices). Poisson performs well on this data (it takes 45 seconds) ScaleSpace encounters some difficulties, probably due to the structured patterns and varying sampling, and misses some triangles (647 seconds). SuperCocone (379 s.) and our method (27 s.) give a similar result. Note the aligned missing triangles in SuperCocone result, at the junction between octree cells.

Figure 17: This dataset is clean (no much noise) but has smallscale bumpy features. The behaviour of the methods is similar to Figure 14: Poisson creates a detailed isosurface, and ScaleSpace manages to extract a manifold mesh. SuperCocone and our method create a smaller number of triangles. However, in our result, missing triangles form small holes, easily filled by our post-processing. the features. ScaleSpace reconstructs a continuous manifold mesh that also follows the bumpy features. SuperCocone and our algorithm give a similar result, with more missing triangles than ScaleSpace. However, the missing triangles are mostly small holes that can be removed by our post-processing step.

Figure 18: A dataset with 11 million vertices, extracted from a reconstruction (courtesy of the Digital Michelangelo project). Since it is not raw scanned data, the only difficulty is the number of points. Poisson creates a nice isosurface, that both respects the details and fills-in the holes (time: 230 seconds). SuperCocone gives a good result. It misses some triangles, probably in junctions between octree cells (time: 1806 seconds). Our algorithm gives a result that is similar to the output of SuperCocone in 110 seconds. ScaleSpace crashed on this dataset.

The 'night' dataset in Figure 18 is a set of vertices from an already reconstructed surface, therefore it has no noise, and the only challenge is its size. Poisson has no problem scaling-up, since its representation (implicit function on an octree) is decoupled from the data. Moreover, it very nicely fills all the holes in the input data. SuperCocone scales up, but still misses some triangles at junctions between octree cells. Scale Space crashed on this dataset.

Finally, we compare and analyse the result of our method for the three largest datasets, with up to 36 million points. For such datasets, existing Delaunay-based method (ScaleSpace and Su-perCocone) start encountering memory limitations. ScaleSpace crashed on all examples. SuperCocone gave an incorrect result on one of them ('facade-val-de-grace') since many triangles were missing and crashed on the others ('chateau-rives' and 'rosetta'). The results are displayed in Figures 19 andFigure 1. A rendering of 'rosetta' is available in the supplemental material. Results of the 'facade-val-de-grace' are shown in Figure 19. For this dataset, Poisson had some problems with some incoherent normals, causing bubbles to appear. It may not be too difficult to fix, but the need for having coherent normals is a difficulty for using Poisson. SuperCocone missed many triangles, we suspect this comes form their implementation that starts to encounter memory problems.

We also show in Figure 1 on the first page, the 'chateau-rives' example with 36 million points. For this dataset, only our algorithm outputs a valid result.

Figure 19: A scanned pointset with 30 million vertices. On this dataset, some normals were not coherent, causing some bubbles in Poisson reconstruction (time: 856 seconds). SuperCocone gave a result in 5011 seconds, but has many missing triangles (maybe it starts reaching memory limits). Our algorithm took 340 seconds.

Conclusion

Based on a simple observation about the restricted Voronoi cell, our algorithm has a trivial implementation, is embarassingly parallel and economic in terms of memory resources, therefore it scales up quite well to datasets with tens of millions vertices while keeping computation times smaller than 5 minutes. Its simplicity makes it possible to embed it into devices such as smartphones and cameras with a CPU (see Android ARMV7 timings in Section 4).

Compared to the state of the art, it is not as good as Poisson for filling the holes, and when ScaleSpace can estimate the correct scale, ScaleSpace has a better set of candidate triangles (see e.g. the eagle dataset). Regarding timings, our algorithm is the fastest, by two order of magnitudes in some cases. However, one needs to take into account that the implementation of SuperCocone that we obtained from its authors is sequential, and could probably be easily parallelized. We think that our algorithm will still be faster since it does not need to compute any Delaunay triangulation (it directly computes the Restricted Delaunay Triangulation).

Observing the relative timings of the different steps of our algorithm, sadly, manifold extraction and post-processing clearly dominate the "core" of the algorithm. We think there is room for improvement. Our implementation of hole filling is not well optimized, not parallel, and we think it could probably be made two or three times faster.

We designed this algorithm mainly with practical considerations in mind. A theoretical study would be interesting. For now, the only property that we can guarantee is that what we compute is exactly the Delaunay triangulation of the input points restricted to the union of the D r i disks (see the observation in Section 3.5). In the results, we observed that the behaviour of our algorithm is very similar to the Cocone algorithm. Therefore, it would be interesting to characterize the configurations where the outputs of both algorithms match. We conjecture that it will be some variant of ε-sampling conditions, that we plan to study in future work.

The algorithm is available in our open source GEOGRAM library [geo].

  10]. These approaches are well-equipped to handle various imperfections in the data. However, normal estimation in the presence of imperfect data remains an open problem [MN03, DLS05]. The interpolation methods elaborate upon Voronoi-Delaunay concepts and come with theoretical guarantees under sampling conditions. Surveys of these algorithms can be found in [CG06, Dey06]. Among these methods there are Tangent plane methods: [Boi84,GK02,CSD04,KY05,BG14], Inside/Outside labelling methods: Power Crust [ACK01], Natural neighbors [BC02], Wrap [Ede03,RS07], Convection [ACA07], Delaunay refinement [BO05] and Restricted Delaunay based methods.

  Dumitriu et al. [DFKM08] used this locally uniform ε-sampling framework to design a surface reconstruction algorithm with correctness guarantees. Although the algorithm is quite complicated, Dumitriu et al. have produced an implementation [DFKM10].

Figure 2 :

 2 Figure 2: Overview of the algorithm.

Figure 3 :

 3 Figure 3: We compute the Voronoi diagram of the input pointset restricted to a union of disks centered on the points and orthogonal to the estimated normals. From left to right: the radius of the disks is increased from 0.2% to 4% of the bounding box diagonal. All the restricted Voronoi cells (colored polygons) are computed in parallel. Whenever three restricted Voronoi cells meet, the three corresponding input points will be connected with a triangle.

Figure 5 :

 5 Figure5: Clipping a disk D r i (in blue) by the bisectors of the neighboring points p ji . The so-defined restricted Voronoi vertex q i j1 j2 and the associated restricted Delaunay triangle (p i , p j1 , p j2 ) are highlighted.

Figure 7 :

 7 Figure 7: (a) Configuration resulting into a T 12 triangle . (b) Configurations of 4 co-circular where T 12 triangles appear.

Figure 8 :

 8 Figure 8: (a) Combinatorial test III. A non-manifold vertex 'by excess'. (b) Combinatorial test IV. Moebius band.

Figure 9 :

 9 Figure 9: (a) Simple bridge. (b) Double bridge.

Figure 10 :

 10 Figure 10: Timing breakdown of the different steps of our algorithm. (a) Small data. (b) Large data with two iterations of smoothing.

Figure 11 :

 11 Figure 11: Performance of the tested algorithms. Up, timings obtained on the Small data. Bottom, timing obtained on the Large data.

Figure 12 :

 12 Figure 12: Comparison of the algorithms with a small dataset (100K vertices) that has high variations of sampling density.

Figure 13 :

 13 Figure 13: Poisson needs consistent normal orientations. ScaleSpace had difficulties estimating the scale, and lost features (finger tips).

Figure 14 :

 14 Figure 14: A noisy dataset (original scanner data of the Stanford Bunny). Poisson both filters-out the noise and recovers the fine features. ScaleSpace successfully extracts a continuous manifold mesh. SuperCocone and our algorithm recover a smaller number of triangles.

Figure 15 :

 15 Figure 15: The same dataset as in Figure 14 with one iteration of smoothing applied. The methods give similar results.

•

  Combinatorial test I: This connectivity test ensures that the triangle is properly connected to the current mesh. Every new triangle should be either incident to two edges of existing triangles, or to one edge of an existing triangles and one isolated point. • Combinatorial test II: This test checks for non-manifold edges and tests whether the three candidate edges are manifold. • Combinatorial test III: This criterion checks that inserting the new triangle do not generate 'by-excess' non-manifold vertices. • Combinatorial test IV: This global orientability test checks if in the neighbourhood of the triangle the same connected component appears with two opposite orientations, then connecting the triangle would create a Moebius strip. The triangle is rejected if it is incident to the same connected component with two different orientations.