Louis Cuel

Jacques-Olivier Lachaud

Quentin Mérigot

Boris Thibert

Robust Geometry Estimation using the Generalized Voronoi Covariance Measure * †

Keywords: Geometric Inference, Normal estimation, Curvature estimation, Distance to a measure, Voronoi Covariance Measure, Power Diagram

The Voronoi Covariance Measure of a compact set K of R d is a tensor-valued measure that encodes geometrical information on K and which is known to be resilient to Hausdorff noise but sensitive to outliers. In this article, we generalize this notion to any distance-like function δ and define the δ-VCM. Combining the VCM with the distance to a measure and also with the witnessed-k-distance, we get a provably good tool for normal estimation that is resilient to Hausdorff noise and to outliers. We present experiments showing the robustness of our approach for normal and curvature estimation and sharp feature detection.

Introduction

The estimation of normal vectors and other differential quantities, such as principal curvature directions and direction of sharp features, has many application in computer graphics, geometry processing or computational geometry. The output of most acquisition devices is a raw point cloud, often corrupted with noise. For applications such as surface reconstruction, it is therefore very useful to be able to perform these estimations directly on such data. Many methods have been proposed for estimating normal and curvature directions directly from a point cloud, such as principal component analysis [START_REF] Hoppe | Surface reconstruction from unorganized points[END_REF], local fitting of polynomials quantities [START_REF] Cazals | Estimating differential quantities using polynomial fitting of osculating jets[END_REF], integral estimates [14], moving least squares [START_REF] Fleishman | Robust moving least-squares fitting with sharp features[END_REF], statistical normal clustering [START_REF] Boulch | Fast and robust normal estimation for point clouds with sharp features[END_REF], to name but a few. Our work is in the continuity of Voronoi-based methods for normal estimations, which have been introduced in [START_REF] Amenta | Surface reconstruction by voronoi filtering[END_REF] and refined by many authors [START_REF] Tamal | Normal and feature approximations from noisy point clouds[END_REF][START_REF] Alliez | Symposium on Geometry processing, Figure 12: Mean absolute curvature on "hand" data estimated with d w P,k -VCM (yellow is low curvature while red is high curvature)[END_REF][START_REF] Mérigot | Voronoi-based curvature and feature estimation from point clouds[END_REF]. These methods are robust to Hausdorff noise but not outliers. Our contribution is to generalize the Voronoi-based approach, and combine it to the notion of distance to a measure [START_REF]Geometric inference for probability measures[END_REF][START_REF] Guibas | Witnessed k-distance[END_REF] so as to gain resilience to outliers.

Distance to a measure All the aforementioned Voronoi-based methods for normal estimation rely directly or indirectly on the notion of distance function. Recall that the distance function to a compact subset K of R d is the function on R d defined by the formula d K (x) = min p∈K px . The stability properties of geometric inference methods based on the distance function is derived from the stability property of the distance function, namely that if K and L are Hausdorff-close, then d K and d L are uniformly close. Unfortunately, geometric data is usually corrupted with outliers, and the hypothesis of Hausdorff noise is not realistic in practice. To make things worse, even the addition of a single outlier to a point cloud can perturb the distance function to that point cloud drastically. To counter this difficulty a robust variant of the notion of distance function to a compact set, called the distance to the measure, was proposed in [START_REF]Geometric inference for probability measures[END_REF]. This new definition of distance is able to cope with the presence of some outliers. Moreover, the distance to a measure is distance-like: this means that it possesses the regularity properties of distance functions to compact sets which makes them amenable to geometric inference.

Contributions

• We extend the definition of Voronoi covariance measure of a compact set. More precisely, we associate to any distance-like function δ, a family of tensor-valued measures called the δ-Voronoi covariance measure (δ-VCM).

• We show the stability of the δ-VCM. Our main general theorem (Theorem 2.1) asserts that if a distance-like function δ well approximates the distance function to a compact set K, then the δ-VCM is close to the VCM of K. When applied to a point cloud P approximating a surface S of R 3 , this implies that one can recover the normal vectors of S very accurately (Proposition 3.1). This estimation is Hausdorff stable and robust to outliers.

• The distance to a measure of a point cloud being not computable in practice, we replace it by the witnessed k-distance [START_REF] Guibas | Witnessed k-distance[END_REF]. We show that the associated VCM still well approximates the VCM of the underlying surface (Proposition 4.1), which opens the door to practical computations.

• We show on various examples that the δ-VCM provides a robust normal estimator resilient to Hausdorff noise and to outliers. For the experiments, we introduce another distancelike function, the median-k-distance. Although we do not have any guarantee for the VCM based on the median-k-distance, it gives very robust estimations in practice. We also use the δ-VCM to estimate curvatures and sharp features. Our estimators improve the results based on the VCM [START_REF] Mérigot | Voronoi-based curvature and feature estimation from point clouds[END_REF] or on the Jet Fitting [START_REF] Cazals | Estimating differential quantities using polynomial fitting of osculating jets[END_REF], even when there are no outliers. They are also compared favorably to the normal classifier of Boulch et al. [START_REF] Boulch | Fast and robust normal estimation for point clouds with sharp features[END_REF].

Notation

In the following we denote by . the Euclidean norm of R d . We will call tensor a square matrix. The tensor product v ⊗ w of two vectors v, w is the d × d matrix whose entries are defined by (v ⊗ w) ij = v i w j . The d-dimensional ball of center x and radius r is denoted by B(x, r).

δ-Voronoi Covariance Measure

As remarked in [START_REF]Geometric inference for probability measures[END_REF], many inference results rely only on two properties of the distance function to a compact set. Any functions satisfying these properties is called distance-like; in particular, the usual distance function to a compact set is distance-like. The goal of this section is to extend the definition of Voronoi Covariance Measure of a compact set, introduced in [START_REF] Mérigot | Voronoi-based curvature and feature estimation from point clouds[END_REF] for the purpose of geometric inference. We associate to any distance-like function δ a tensor-valued measure called the δ-Voronoi covariance measure or δ-VCM. Our main theoretical result is Theorem 2.1, which asserts in a quantitative way that if a distance-like function δ is uniformly close to the distance function to a compact set, then the δ-VCM is close to the VCM of this compact set.

Informally, this theorem shows that one can recover geometric information about a compact set using an approximation of its distance function by a distance-like function.

δ-Voronoi Covariance Measure

In this paragraph, we introduce the definitions necessary for the precise statement of the main theorem. We start by the definition of distance-like function, following [START_REF]Geometric inference for probability measures[END_REF]. Note that we used the remark following Proposition 3.1 in [START_REF]Geometric inference for probability measures[END_REF] to drop the 1-Lipschitz assumption, which follows from the two other assumptions.

Definition 2.1 (Distance-like function). A function

δ : R d → R + is called distance-like if • δ is proper, i.e. lim x →∞ δ(x) = ∞. • δ 2 is 1-semiconcave, that is δ 2 (.) -. 2 is concave.
The typical examples of distance-like functions that we will consider are power distances. Given a point cloud P and a family of non-negative weights (ω p ) p∈P , we call power distance to P the distance-like function δ P defined by

δ P (x) := min p∈P x -p 2 + ω p 1/2 . ( 1
)
Note that when the weights are all zero, the power distance is nothing but the distance function to the point cloud P .

Definition 2.2 (δ-VCM). The δ-Voronoi Covariance Measure is a tensor-valued measure.

Given any non-negative probe function χ, i.e. an integrable function on R d , we associate a positive semi-definite matrix defined by

V δ,R (χ) := δ R n δ (x) ⊗ n δ (x).χ (x -n δ (x)) dx, (2)
where δ R := δ -1 ((-∞, R]) and where n δ (x) := 1 2 ∇δ 2 (x). Note that this vector field n δ is defined at almost every point in R d by the 1-concavity property of distance-like functions.

The tensor V δ,R (χ) is a convolution of a tensor with the function χ that localizes the calculation on the support of χ. Intuitively, if χ is the indicatrix of a ball B ⊂ R d , then the VCM V δ,R (χ) is the integral of n δ ⊗ n δ over the set of points of δ R that are "projected" into the ball B.

Example 2.1 (δ-VCM of a distance function). When considering the distance function to a compact set K, the set d R K coincides with the offset of K with radius R, thus explaining our choice of notation. Moreover, the vector field n d K has an explicit expression in term of the projection function p K on K, that is the application that maps a point to its closest neighbor in K:

n d K (x) = x -p K (x). ( 3
)
Comparing Eq. (2) to Eq. (4.2) in [START_REF] Mérigot | Voronoi-based curvature and feature estimation from point clouds[END_REF] and using this remark, one sees that the d K -VCM defined here matches the original definition of the VCM of K:

V d K ,R (χ) = K R (x -p K (x)) ⊗ (x -p K (x)).χ(p K (x))dx.
Consider a smooth compact surface S of R 3 , with exterior unit normal n. If R is chosen small enough the following expansion holds as r → 0, where . op is the operator norm [START_REF] Mérigot | Voronoi-based curvature and feature estimation from point clouds[END_REF]:

V d S ,R (1 B(p,r) ) - 2π 3 R 3 r 2 [n(p) ⊗ n(p)] op = O(r 3 ). (4)
This equation shows that one can recover local information about differential quantities from the VCM of a surface. Note that curvature information is also present in the matrix.

Example 2.2 (δ-VCM of a power distance). Let us give a closed form expression for the VCM of a power distance, which we will use in the computations. Each power distance defines a decomposition of the space into a family of convex polyhedra, called power cells, on which the function δ 2 P is quadratic. The power cell of the point p in P is defined by

Pow P (p) = {x ∈ R d ; ∀q ∈ P, x -p 2 + ω p ≤ x -q 2 + ω q }. (5)
When the weight vector ω vanishes, we recover the notion of Voronoi cell. The following Lemma generalizes Eq. (2.1) in [START_REF] Mérigot | Voronoi-based curvature and feature estimation from point clouds[END_REF], and shows that computing the VCM of a power distance amounts to computing the covariance matrix of the intersection of each power cell with a ball (see also Algorithm 1).

Lemma 2.1. Let (P, ω) be a weighted point cloud. Given a probe function χ,

V δ P ,R (χ) = p∈P χ(p) M p , (6)
where M p is the covariance matrix of

C p := Pow P (p) ∩ B(p, (R 2 -ω p ) 1/2 ), M p := Cp (x -p) ⊗ (x -p)dx. (7)
Figure 1: Integration domain of δ P -VCM. The weighted point cloud P is represented by a union of circles whose radius are the weights. We suppose χ to be the indicatrix of the ball B(y, r) of center y and radius r. The number of power cells is less than the cardinal of the point cloud since one cell is empty.

Proof. Since for every point x in the interior of the power cell Pow P (p),

n δ P (x) = x -p, hence χ(x -n δ P (x)) = χ(p
) is constant and one can decompose the integral over power cells

V δ P ,R (χ) = p∈P χ(p) Pow P (p)∩δ R P (x -p) ⊗ (x -p)dx.

Furthermore a computation gives Pow

P (p) ∩ δ R P = Pow P (p) ∩ B p, (R 2 -ω p ) 1/2 , if we consider that R 2 -ω p < 0 defines the empty ball.

Stability of the δ-VCM

We are now able to state our main theorem, which is a generalization of the stability theorem for VCM proven in [START_REF] Mérigot | Voronoi-based curvature and feature estimation from point clouds[END_REF]. Our theorem asserts that if the distance function to a compact set K is well approximated by a distance-like function δ, then the VCM of K is also well approximated by the δ-VCM. The hypothesis of this theorem is satisfied for instance under the sampling model considered in [START_REF] Guibas | Witnessed k-distance[END_REF].

The uniform norm of a function

f on R d is denoted f ∞ = sup R d |f |. Given a Lipschitz function f with Lipschitz constant Lip(f ), we introduce its bounded-Lipschitz norm f BL = f ∞ + Lip(f ).
Theorem 2.1. Let K be a compact set and δ a distance-like function. For any bounded Lipschitz function χ : R d -→ R + , one has

V δ,R (χ) -V d K ,R (χ) op ≤ C 1 χ BL δ -d K 1 2 ∞ ,
where C 1 is a constant that only depends on R, d and diam(K).

In practice choosing a probe function χ supported in a small ball allows one to recover local information from the δ-VCM. Remark 1. A notable feature of this theorem is that the constant in the upper bound depends only on the diameter of K and not on its local geometry or on its regularity.

We now recall a simplified version of Davis-Kahan Theorem [START_REF] Davis | The rotation of eigenvectors by a perturbation[END_REF].

Theorem 2.2 (Davis-Kahan). Let M and M be two symmetric matrices, λ an eigenvalue of M and δ λ (M ) be the minimum distance between λ and any other eigenvalue of M . Then for every eigenvector v of M with eigenvalue λ, there exists an eigenpair (λ , v ) of M such that

|λ -λ | ≤ √ 2 M -M op and v -v ≤ √ 2 M -M op δ λ (M ) , provided that M -M op ≤ √ 2 δ λ (M ).
Remark 2. When the compact set K is a smooth surface S and χ r p is the indicatrix of a ball centered at a point p of S, Equation ( 4) implies that the eigenvector associated to the highest eigenvalue of V d K ,R (χ r p ) gives an estimation of the normal vector to S at p. Indeed, λ = 1 is the only non-vanishing eigenvalue of the matrix M := n(p) ⊗ n(p) and δ λ(M ) = 1. Hence, by Davis-Kahan theorem, and under reasonable assumptions, an eigenvector associated to the highest eigenvalue of V δ,R (χ r p ) thus gives an estimation of the normal direction of S at the point p.

Stability of gradients of distance-like functions

We mention in this subsection an intermediary result, Corollary 2.1, that guarantees the L 1stability of gradients of distance-like functions. It is a consequence of Theorem 3.5 of [START_REF] Chazal | Boundary measures for geometric inference[END_REF] which gives the stability of gradients of convex functions. For any function f : R d → R k and any set E of R d , we put

f 1,E = E f (x) dx and f ∞,E = sup x∈E f (x) .
If E is rectifiable, we denote by H n (E) its n-dimensional Hausdorff measure. If f is differentiable almost everywhere, one puts ∇f (E) := {∇f (x), x ∈ E ∩Ω f } where Ω f is the set of points where f is differentiable. The diameter of a set X ⊂ R d is given by diam(X) := sup x,y∈X xy . We first recall the following theorem.

Theorem 2.3 (Theorem 3.5 of [START_REF] Chazal | Boundary measures for geometric inference[END_REF]). Let E be an open subset of R d with rectifiable boundary, and f , g be two convex functions from

E to R such that diam(∇f (E) ∪ ∇g(E)) ≤ D. Then ∇f -∇g 1,E ≤ C 2 H d (E) + (D + f -g 1 2 ∞,E )H d-1 (∂E) f -g 1 2 ∞,E ,
where C 2 is a constant that only depends on the dimension d.

Corollary 2.1. Let E be a set of R d with rectifiable boundary and ε > 0. For any distance-like functions δ and δ such that δδ ∞,E ≤ ε, one has

n δ -n δ 1,E ≤ C 3 [H d (E) + (2diam(E) + 4R + √ 2Rε) H d-1 (∂E)] √ 2Rε, where R = max( δ ∞,E , δ ∞,E
) and C 3 is a constant depending only on d.

We introduce ψ δ (x) = x 2 -δ 2 (x) and ψ δ (x) = x 2 -δ 2 (x). For almost every x, one has

n δ (x) -n δ (x) = 1 2 (∇δ 2 (x) -∇δ 2 (x)) = 1 2 (∇ψ δ (x) -∇ψ δ (x)).
Using the convexity of ψ δ and ψ δ and Theorem 2.3, Corollary 2.1 follows from the following lemma.

Lemma 2.2. Under the assumptions of Corollary 2.1, one has

(i) ∀x ∈ E, |ψ δ (x) -ψ δ (x)| ≤ 2Rε (ii) diam(∇ψ δ (E) ∪ ∇ψ δ (E)) ≤ 2diam(E) + 4R.
Proof. For every point x in E, one has

|ψ δ (x) -ψ δ (x)| = |δ 2 (x) -δ 2 (x)| = |δ(x) -δ (x)| × |δ(x) + δ (x)| ≤ 2Rε.
Let now X and X be points in ∇ψ δ (E) and ∇ψ δ (E) respectively. There exist x, x in E such that X = 2x -2δ(x)∇δ(x) and X = 2x -2δ (x )∇δ (x ). Then

X -X || ≤ 2x -2x + 2δ (x )∇δ (x ) -2δ(x)∇δ(x) ≤ 2 diam(E) + 2R ∇δ ∞,E + 2R ∇δ ∞,E ≤ 2 diam(E) + 4R.
In the last inequality, we used the fact that a distance-like function is 1-Lipschitz (see Proposition 3.1 in [START_REF]Geometric inference for probability measures[END_REF]). The result still holds if X, X both belong to ∇ψ δ (E) or to ∇ψ δ (E).

Proof of Theorem 2.1

This proof follows the proof of the stability theorem in [START_REF] Mérigot | Voronoi-based curvature and feature estimation from point clouds[END_REF]. The idea is to compare the two integrals on the common set

E = K R-ε with ε = δ -d K ∞,E
and to show that remaining parts are negligible. By Proposition 4.2 of [START_REF] Chazal | Boundary measures for geometric inference[END_REF], the set ∂E is rectifiable. For every R ∈ R + , we denote by N (∂K, R) the covering number of ∂K with a radius parameter R, namely the minimal number of balls of radius R needed to cover ∂K. Let us first suppose that ε < R 2 . We have

V d K ,R (χ) = E n d K (x) ⊗ n d K (x)χ(x -n d K (x))dx + K R \E n d K (x) ⊗ n d K (x)χ(x -n d K (x))dx. For every x ∈ K R , one has n d K (x) = x -p K (x) ≤ R. The fact that χ ∞ ≤ χ BL implies V d K ,R (χ) - E n d K (x) ⊗ n d K (x)χ(x -n d K (x))dx op ≤ R 2 • χ BL • H d (K R \K R-ε ). (8)
We proceed similarly for the δ-VCM. By definition, we have

n δ (x) = δ(x)∇δ(x) ≤ |δ(x)| ≤ R, since δ is 1-Lipschitz. Moreover, one has K R-ε ⊂ δ R ⊂ K R+ε . Hence, the volume of δ R \ K R-ε is less than the volume of K R+ε \K R-ε and V δ,R (χ) - E n δ (x) ⊗ n δ (x)χ(x -n δ (x))dx op ≤ R 2 • χ BL • H d (K R+ε \K R-ε ). ( 9
)
We now bound the volume of K R+ε \K R-ε by using Proposition 4.2 of [START_REF] Chazal | Boundary measures for geometric inference[END_REF]. We set ω n (t) to be the volume of the n-dimensional ball of radius t.

H d (K R+ε \K R-ε ) = R+ε R-ε H d-1 (∂K t ) dt ≤ R+ε R-ε N (∂K, t)ω d-1 (2t) dt ≤ N (∂K, R -ε)ω d-1 (2(R + ε))2ε ≤ 2 N (∂K, R 2 )ω d-1 (3R) ε. (10)
We now need to bound the operator norm of the integral

E [A d K (x) -A δ (x)
]dx, where

A d K (x) := n d K (x) ⊗ n d K (x)χ(x -n d K (x)),
and A δ is defined similarly. The difference of these two terms is decomposed as:

A d K (x) -A δ (x) = χ(x -n d K (x))(n d K (x) ⊗ n d K (x) -n δ (x) ⊗ n δ (x)) + [χ(x -n d K (x)) -χ(x -n δ (x))]n δ (x) ⊗ n δ (x)
From the facts that

n d K (x) ⊗ n d K (x) -n δ (x) ⊗ n δ (x) ≤ n d K (x) ⊗ (n d K (x) -n δ (x)) + (n δ (x) -n d K (x)) ⊗ n δ (x) ≤ 2R n d K (x) -n δ (x) ,
and also

[χ(x -n d K (x)) -χ(x -n δ (x))]n δ (x) ⊗ n δ (x) ≤ Lip(χ) n d K (x) -n δ (x) R 2 one has A d K (x) -A δ (x) op ≤ χ BL (R 2 + 2R) n d K (x) -n δ (x) ,
which by integration leads to

E [A d K (x) -A δ (x)]dx op ≤ χ BL (R 2 + 2R) n d K -n δ 1,E .
Corollary 2.1 implies that the norm n d K -n δ 1,E is bounded by

C 2 × [H d (E) + (2diam(E) + 4R + √ 2Rε) × H d-1 (∂E)] × √ 2Rε.
Since diam(E) ≤ diam(K) + 2R and ε ≤ R 2 , the two last equations imply

E [A d K (x) -A δ (x)]dx op ( 11
)
≤ C 2 χ BL (R 2 + 2R)[H d (K R ) + (2diam(K) + 9R)H d-1 (∂K R-ε )] √ 2R √ ε.
Again by Proposition 4.2 of [START_REF] Chazal | Boundary measures for geometric inference[END_REF], we bound

H d-1 (∂K R-ε ) by N (∂K, R 2 )•ω d-1 (3R). Furthermore, since the diameter of K R is less than diam(K) + 2R, its volume H d (K R ) is bounded by a constant involving diam(K), R and d. Remark that N (∂K, R
2 ) is also bounded by a constant depending on the same quantities. Hence, the constants involved in Equations ( 11) and ( 10) also depend solely on diam(K), R and d. We conclude the proof of Theorem 2.1 by combining [START_REF] Guibas | Witnessed k-distance[END_REF] with [START_REF] Davis | The rotation of eigenvectors by a perturbation[END_REF] and [START_REF] Tamal | Normal and feature approximations from noisy point clouds[END_REF].

The case where ε ≥ R 2 is trivial (and not interesting in practice) since both V δ,R (χ) and V d K ,R (χ) can be upper bounded by a quantity depending on diam(K), R and d, which can be put into the constant C 1 of the Theorem.

Figure 2: Distance to a measure. We suppose here that the measure µ is supported on a compact set K. The quantity δ µ,m (x) is the minimal radius r such that B(x, r) ∩ K has a mass m.

VCM using the distance to a measure

The distance to a measure is a distance-like function that is known to be resilient to outliers. It is therefore natural to consider the δ-VCM in the particular case where δ is a distance to a measure. First we recall the definition and some stability results of the distance to a measure. Then we study the δ-VCM when δ is a distance to a measure.

Distance to a measure

The distance to a measure has been introduced in [START_REF]Geometric inference for probability measures[END_REF] and is defined for any probability measure µ on R d . We denote in the following supp(µ) the support of µ. Definition 3.1. Let µ be a probability measure on R d and m 0 a regularization parameter in (0, 1). The distance function to the measure µ is defined for every point x in R d by

d µ,m0 (x) := 1 m 0 m0 0 δ 2 µ,m (x)dm 1/2 , ( 12
)
where δ µ,m (x) = inf{r ≥ 0, µ(B(x, r)) ≥ m}.

The formula defining the function δ µ,m mimics a similar formula for the distance function to a set K: d K (x) = inf{r ≥ 0, K ∩ B(x, r) = ∅}, as shown on Figure 2. However it turns out that the function δ µ,m is not distance-like [START_REF]Geometric inference for probability measures[END_REF], while the regularization defined by Eq. ( 12) is distance-like. Furthermore, the distance to the measure has been shown to be resilient to outliers, and more precisely, Theorem 3.5 of [START_REF]Geometric inference for probability measures[END_REF] states that

d µ,m0 -d µ ,m0 ∞ ≤ 1 √ m 0 W 2 (µ, µ ),
where W 2 is the 2-Wasserstein distance between measures. For more details on Wasserstein distances, which are also known under the name Earthmover distances in image retrieval and computer science, one can refer to [15]. To give an intuition, when µ K and µ K are uniform probability measures on two point clouds K and K with the same number of points, the distance W 2 (µ K , µ K ) is the square root of the cost of a least square assignment between K and K . Point cloud case Let P ⊂ R d be a finite point set with n points, k ∈ (0, n) a real number and m 0 = k/n. Following [START_REF] Guibas | Witnessed k-distance[END_REF], we call k-distance to P and denote d P,k the distance to the measure for the uniform measure on P for the parameter m 0 . In the particular case where k is an integer, a simple calculation [START_REF]Geometric inference for probability measures[END_REF] shows that for every point x in R d ,

d 2 P,k (x) = 1 k pi∈N P,k (x) x -p i 2 ,
where N P,k (x) are the k nearest neighbors of x in P . Furthermore, the k-distance d P,k is a power distance [11, Proposition 3.1]. More precisely, if we denote Bary P,k the set of isobarycenters of k distinct points in P , one has

∀x ∈ R d d 2 P,k (x) = min b∈Bary P,k x -b 2 + ω b ,
where the weight ω b = 1 k pi∈N P,k (b) bp i 2 . Figure 3 illustrates the stability of the kdistance of a point cloud with respect to outliers. Remark that the level sets of the k-distance (second picture), are much smoother and faithful than the level sets of the distance function (first picture).

Stability

We state here a stability result obtained by combining our result with a stability result for the distance to a measure established in [START_REF] Guibas | Witnessed k-distance[END_REF]. Although this result can hold for measures supported in R d , we choose to state it for uniform measures on surfaces of R 3 , so as to give a better intuition.

Let S be a two dimensional surface and µ S denote the uniform probability measure on S. An application of Günther-Bishop theorem which can be found explicitely on page 749 of [START_REF]Geometric inference for probability measures[END_REF] is that there exists a constant α S > 0 such that

∀p ∈ S, ∀r ≤ diam(S), µ S (B(p, r)) ≥ α S r 2 . ( 13
)
Proposition 3.1. Let P be a point cloud of cardinal n, S a surface of R 3 and let µ P and µ S be the uniform probability measures on P and S respectively. If we define k = W 2 (µ S , µ P ) √ α S n, then we have

V d P,k ,R (χ) -V d S ,R (χ) op ≤ C 4 α -1/8 S χ BL W 2 (µ S , µ P ) 1 4 ,
where the constant C 4 depends on diam(S) and R.

Proof. Let k ∈ (0, n) and m 0 = k/n. By Theorem 1 of [START_REF] Guibas | Witnessed k-distance[END_REF] (Section 4.2), one has

d P,k -d K ∞ ≤ m -1 2 0 W 2 (µ S , µ P ) + α -1 2 S m 1 2 0 .
The right hand-side is minimal when m 0 satisfies m

-1 2 0 W 2 (µ S , µ P ) = α -1 2 S m 1 2 0 , namely m 0 = W 2 (µ S , µ P ) √ α S . If k = W 2 (µ S , µ P )
√ α S n, the bound in the previous equation becomes

W 2 (µ S , µ P ) 1 2 α -1/4 S
. We conclude by Theorem 2.1.

Remark this result can be easily extended to any dimension d and also to the case where S is a compact set with a dimension at most l (i.e. there exists a constant α S such that µ S (B(p, r)) ≥ α S r l ).

Tractable variants for the d P,k -VCM

We have seen in the previous section that the d P,k -VCM is resilient to outliers. Unfortunately, it is not practically computable, because the computation of the k-distance is not. We therefore study in this section variants of the d P,k -VCM, by using variants of the k-distance. In Section 4.1, we propose an efficient relaxation of the d P,k -VCM for which we have a stability result. In Section 4.2, we investigate another δ-VCM interesting for its resilience to outliers.

Witnessed-k-distance and d w P,k -VCM

Let P ⊂ R d be a point cloud, k an integer. We have seen in the previous section that the k-distance is a power distance on the set Bary P,k of k points of P . Lemma 2.1 then allows us in theory to calculate the d P,k -VCM. Unfortunately, it is not computable in practice. The set Bary P,k is indeed huge since its cardinal is of the order of n k , if n is the cardinal of P . To overcome this problem, we use the witnessed k-distance that was introduced in [START_REF] Guibas | Witnessed k-distance[END_REF].

Relaxation with the witnessed k-distance A point is said to be a k witness of P if it is the barycenter of k points p 1 , . . . , p k of P such that p 2 , . . . , p k are the (k -1)-nearest neighboors of p 1 in P \ {p 1 }. We denote by Bary w P,k this set of points. The witnessed k-distance is then defined as the following power distance

d w P,k (x) := min b∈Bary w P,k x -b 2 + ω b 1/2
,

where ω b = 1 k pi∈N P,k (b) b -p i 2 .
We then use Lemma 2.1 for the computation of its VCM. Note that the number of cells in the power diagram is bounded by the number of points in P , and can therefore be computed efficiently using e.g. CGAL [START_REF] Cgal | [END_REF]. We also have a stability result.

→ R + V d w P,k ,R (χ) -V S,R (χ) op ≤ C 5 α -1/8 S χ BL W 2 (µ P , µ S ) 1 4 ,
where the constant C 5 depends on diam(S) and R.

The proof of this proposition is similar to the one of Proposition 3.1, and is using Theorem 4 of [START_REF] Guibas | Witnessed k-distance[END_REF].

Median-k-distance and d m P,k -VCM

We introduce here the median-k-distance which is a distance-like function derived from the witnessed-k-distance. We do not have yet stability results for this function, but numerical experiments in the next section show that the corresponding δ-VCM is very stable in practice (see also the visual intuition given by Figure 2). Let P be a point cloud and k > 1 be an integer. A point is said to be a median k-witness of P if it is the geometric median of k points p 1 , . . . , p k of P such that p 2 , . . . , p k are the (k -1)-nearest neightbors of p 1 in P \ {p 1 }. We denote by Bary m P,k this set of points.

Definition 4.1. The median-k-distance is the power distance defined by

d m P,k (x) := min b∈Bary m P,k x -b 2 + ω b 1/2
,

where

ω b = 1 k pi∈N P,k (b) b -p i 2 .
The key idea is to replace the barycenter of k points involved in the witnessed-k-distance by the geometric median. This can be seen as replacing the L 2 -norm by the L 1 -norm. Indeed, it is well-known that the barycenter of p 1 , . . . , p k is the point b that minimizes

d i=1 b -p i 2 .
Similarly, a geometric median is a point b that minimizes d i=1 b-p i . Note that the geometric median is unique when the k points are not colinear.

As for the witnessed-k-distance, one can compute the d m P,k -VCM by using Lemma 2.1.

Computation and Experiments

Computation of the VCM for a power distance

We describe here our algorithm to compute an approximation of the VCM of a power distance. Let thus (P, ω) be a weighted point cloud that defines a distance-like function δ P . For the purpose of normal and curvature estimation, we have to compute for each point q in P the covariance matrix V δ P ,R (χ r q ), where χ r q is the indicatrix of the ball B(q, r). Using Lemma 2.1, we have

V δ P ,R (χ r q ) = p∈P ∩B(q,r) M p ,
where the matrix M p is defined by

C p = Pow P (p) ∩ B p, (R 2 -ω p ) 1/2 M p = Cp (x -p) ⊗ (x -p)dx
The main difficulty is to compute these covariance matrices, and in practice we approach it by replacing the ball in the definition of C p by a convex polyhedron. The input of our algorithm (summarized in Algorithm 1) is a weighted point cloud, a radius R and an approximation of the unit ball by a convex polyhedron B. In all examples, the polyhedron B was a dodecahedron circumscribed to the unit ball. To compute the approximate covariance matrices (M B p ) p∈P , we first build the power diagram of the weighted point cloud. Then, we proceed in two steps for each point:

Intersection

We compute an approximation of the cell C p by

C B p := Pow P (p) ∩ (p + (R 2 -ω p ) 1/2 B) (14)
To perform this computation, we gather the half-spaces defining both polyhedra and we compute their intersection. By duality, this is equivalent to the computation of a convex hull.

Integration

The boundary of the polyhedron C B p is decomposed as a union of triangles, and we consider the tetrahedron ∆ 1 p , . . . , ∆ kp p joining these triangles to the centroid of C p . We compute

M B p = C B p (x -p) ⊗ (x -p)dx (15)
by summing the signed contribution of each tetrahedron, which can be evaluated exactly using the same formulas as in [START_REF] Alliez | Symposium on Geometry processing, Figure 12: Mean absolute curvature on "hand" data estimated with d w P,k -VCM (yellow is low curvature while red is high curvature)[END_REF]. We implemented this algorithm using CGAL [START_REF] Cgal | [END_REF] for the power diagram calculation and the intersection step. We report in Table 1 some running times for the witnessed-k-distance and the median-k-distance.

Algorithm 1 Computation of V δ P ,R .

Require: P ⊆ R d point cloud, R > 0, B = approximation of B(0, 1) Computation of the power diagram (Pow P (p)) p∈P for all p ∈ P do

C B p ← Pow P (p) ∩ (p + (R 2 -ω p ) 1/2 B) { §5.1.1} ∆ 1 p , . . . , ∆ kp p ← decomposition of C B p into tetrahedra { §5.1.2} M B p ← kp i=1 ∆ i p (x -p) ⊗ (x -p)dx { §5.1.2} end for return (M B p ) p∈P .

d w P,k -VCM evaluation

We show in this section that the d w P,k -VCM provides in practice a robust estimator for normal vectors and curvatures, as well as for feature detection. It is also resilient to both Hausdorff noise and outliers. Note that in our experiments, the point cloud P is the set of vertices of a mesh. The mesh itself is only used for visualisation, our algorithm being only based on the point cloud P . The diameter of the point cloud in all our examples is D = 2. We say that the Hausdorff noise is ε if every point p of P can be randomly moved at a distance less than ε. We will speak about outliers if in addition a certain amount of points can be moved much further.

Normal estimation

As suggested by Remark 2, we define a normal at each point p of a point cloud P as the eigenvector associated to the largest eigenvalue of V δ P ,R (χ r p ). Comparison with other methods We compare the accuracy of our method to the VCM method [START_REF] Mérigot | Voronoi-based curvature and feature estimation from point clouds[END_REF], to the Jet Fitting method [START_REF] Cazals | Estimating differential quantities using polynomial fitting of osculating jets[END_REF] and to the Randomized Hough Transform method [START_REF] Boulch | Fast and robust normal estimation for point clouds with sharp features[END_REF]. Figure 4 displays the average angular deviation between the estimated normal and the ground truth in the case of a noisy ellipsoid with outliers. Our method gives better result than the others.

Sensitivity to parameters To measure the dependence of our method on the different parameters, we computed the average deviation of the normal on a noisy ellipsoid for different choices of parameters k and R. As can be seen in Figures 5 and6, the results are essentially stable. Figure 5 shows that we should choose k larger when the noise is large. We observe that the value k = 30 gives good results in both presence and absence of noise. Therefore all our experiments (except of course the ones of Figure 5) are done with the value k = 30. Note that the classical VCM method coincides exactly with the d w P,k -VCM method when k = 1.

Visualisation of normal estimation

We test the quality of our normal estimator on standard shapes, by rendering meshes according to their estimated normals. In Figure 7, we notice

Estimation of curvatures and features detection

The covariance matrix also carries curvature information along other eigendirections [START_REF] Mérigot | Voronoi-based curvature and feature estimation from point clouds[END_REF]. We denote by λ 0 ≥ λ 1 ≥ λ 2 the three eigenvalues of V δ P ,R (χ r p ) at a point p. Up to a multiplicative constant, λ 1 and λ 2 correspond to the absolute value of the minimal and maximal curvatures respectively [START_REF] Mérigot | Voronoi-based curvature and feature estimation from point clouds[END_REF]. We call the corresponding eigenvectors respectively minimal and maximal principal directions.

We compare our method with the Jet Fitting method [START_REF] Cazals | Estimating differential quantities using polynomial fitting of osculating jets[END_REF] in Figure 9. Experiments have been done for a large choice of different parameters, and we find better principal directions with our method. In the presence of many outliers, we plot in Figure 10 the minimal principal direction estimation of d w P,k -VCM projected on the initial mesh. To illustrate the performance of curvature estimation by d w P,k -VCM, the estimated mean absolute curvature λ 1 +λ 2 is displayed for the standard "caesar" and "hand" models (Figures 11 and12). On Figure 11, we also render the triangulation before any processing to better illustrate the amplitude of noise.

We also detect sharp features, using the same criteria as in [START_REF] Mérigot | Voronoi-based curvature and feature estimation from point clouds[END_REF]: we say that a point p belongs to a sharp feature if λ 1 /(λ 0 + λ 1 + λ 2 ) ≥ T , for a given threshold T . Figure 13 shows results on fandisk data with increasing noise (Hausdorff+outliers). It is clear that both witnessed-k-distance VCM and median-k-distance VCM are much less sensitive to outliers than VCM. More precisely our methods classify an outlier as a feature only if it is close to a sharp feature, while the VCM fails. This behavior is reinforced when the noise level increases. [15] Cédric Villani, Optimal transport: old and new, vol. 338, Springer, 2008.

Figure 3 :

 3 Figure 3: Level sets of distance functions of a noisy point cloud. From left to right: usual distance function to the point cloud; k-distance; witnessed-k-distance; median-witnessed-kdistance

Proposition 4 . 1 .

 41 Let µ S be the uniform measure on a surface S ⊂ R 3 and µ P denotes the uniform probability measure on a point cloud P . Let k = W 2 (µ P , µ S ) √ α S n/4, where n is the number of points of P . Then for any bounded Lipschitz function χ : R d

  VCM R = 0.1, r = 0.2 VCM, R = 0.1, r = 0.1 d w P,k -VCM R = 0.1, r = 0.2, k = 30 d w P,k -VCM R = 0.1, r = 0.1, k = 30 Jet fitting k = 50 Jet fitting k = 80 Hough Transform k = 400 Hough Transform k = 700

Figure 4 :

 4 Figure 4: Ellipsoid normal estimation. Left: the rendering illustrates the quality of the normal estimation (parameters: k = 30, R = 0.2D, r = 0.2D, Hausdorff noise = 0.4D), where D is the diameter of the original shape. Right: Comparison of the average normal deviation for Jet Fitting, VCM, d w P,k -VCM and the Randomized Hough Transform method [4] for different values of Hausdorff noise.

Figure 7 :

 7 Figure 7: Rendering of "caesar" data using triangle normal (top row) and estimated d w P,k -VCM normal (bottom row) with Phong shading (parameters R = 0.04D, r = 0.04D, k = 30, where D is the diameter of the original shape). From left to right, the Hausdorff noise is 0.02D, 0.04D and 0.06D.

Figure 8 :

 8 Figure 8: Rendering comparison of "hand" data: (left) using our d w P,k -VCM normal estimation with parameters R = 0.04D, r = 0.02D, k = 30 (D is the diameter of the original shape), and (right) the Randomized Hough Transform method [4] with standard parameters: neighbors number = 500 and the others equal to default values. Triangles are displayed with Phong shading and "hand" is perturbated by a Hausdorff noise equal to 0.04D.

Figure 9 :

 9 Figure 9: Principal direction estimation on Bimba datas (no noise). We put a line-segment aligned at every point to the minimal principal direction. Left: Jet-fitting method (with a polynomial of degree 4 × 4 and k = 100 neighbors for each point). Right: d w P,k -VCM method (parameters R = 0.04D, r = 0.08D, k = 30, where D is the diameter of the original shape.)

Figure 10 :

 10 Figure 10: Principal direction estimation on Bimba datas of diameter D. Left: input data with outliers (80% of the points are moved at a distance at most 0.02D, 10% are moved at a distance between 0.02D and 0.1D and 10% are outliers taken randomly in the bounding box). Right: for every vertex, we project the minimal principal direction estimation of the d w P,k -VCM on the initial mesh.

Figure 11 :Figure 13 :

 1113 Figure 11: Mean absolute curvature on "ceasar" data (of diameter D) estimated with d w P,k -VCM (yellow is low curvature while red is high curvature). Parameters are set as R = 0.04D, r = 0.04D, k = 30. From left to right: input data with Hausdorff noise = 0.02D; corresponding curvature estimation; input data with Hausdorff noise = 0.06D; corresponding curvature estimation.

Figure 15 :

 15 Figure 15: Comparison between d w P,k -VCM and d m P,k -VCM (for the same parameters R = 0.1D, r = 0.1D, k = 30, where D is the diameter of the original shape).Hausdorff noise is 0.2D. The percentage of "outliers" corresponds to the percentage of points that have been moved at a distance in-between 0.05D and 0.1D. Left: input mesh with 50% of outliers. Right: evolution of the angle error when the proportion of outliers increases.

This research has been supported in part by the ANR grants DigitalSnow ANR-11-BS02-009, KIDICO ANR-2010-BLAN-0205 and TopData ANR-13-BS01-0008.

† An extended abstract of this work has been presented at EuroCG 2014.

that the rendering done with normals computed with our method is much better than the rendering done with normals induced by the geometry of the underlying mesh, at the same time robust to noise while keeping intact significant features. In Figure 8, we compare this rendering with the rendering obtained when the normal is calculated by the Randomized Hough Transform method [START_REF] Boulch | Fast and robust normal estimation for point clouds with sharp features[END_REF], which is a statistic method known to be resilient to noise and outliers. Our method achieves a much smoother rendering.