Submanifold Reconstruction Jean-Daniel Boissonnat DataShape, INRIA - - PowerPoint PPT Presentation

Submanifold Reconstruction

Jean-Daniel Boissonnat DataShape, INRIA http://www-sop.inria.fr/datashape

Algorithmic Geometry Submanifold reconstruction J-D. Boissonnat 1 / 36

Geometric data analysis

Images, text, speech, neural signals, GPS traces,...

Geometrisation : Data = points + distances between points Hypothesis : Data lie close to a structure of “small” intrinsic dimension Problem : Infer the structure from the data

Algorithmic Geometry Submanifold reconstruction J-D. Boissonnat 2 / 36

Submanifolds of Rd

A compact subset M ⊂ Rd is a submanifold without boundary of (intrinsic) dimension k < d, if any p ∈ M has an open (topological) k-ball as a neighborhood in M

W U Rm φ RN M

Intuitively, a submanifold of dimension k is a subset of Rd that looks locally like an open set of an affine space of dimension k A curve a 1-dimensional submanifold A surface is a 2-dimensional submanifold More generally, manifolds are defined in an intrinsic way,

d

Algorithmic Geometry Submanifold reconstruction J-D. Boissonnat 3 / 36

Triangulation of a submanifold

We call triangulation of a submanifold M ⊂ Rd a simplicial complex ˆ M such that ˆ M is embedded in Rd its vertices are on M it is homeomorphic to M Submanifold reconstruction The problem is to construct a triangulation ˆ M of some unknown submanifold M given a finite set of points P ⊂ M

Algorithmic Geometry Submanifold reconstruction J-D. Boissonnat 4 / 36

Triangulation of a submanifold

We call triangulation of a submanifold M ⊂ Rd a simplicial complex ˆ M such that ˆ M is embedded in Rd its vertices are on M it is homeomorphic to M Submanifold reconstruction The problem is to construct a triangulation ˆ M of some unknown submanifold M given a finite set of points P ⊂ M

Algorithmic Geometry Submanifold reconstruction J-D. Boissonnat 4 / 36

Issues in high-dimensional geometry

Dimensionality severely restricts our intuition and ability to visualize data ⇒ need for automated and provably correct methods methods Complexity of data structures and algorithms rapidly grow as the dimensionality increases ⇒ no subdivision of the ambient space is affordable ⇒ data structures and algorithms should be sensitive to the intrinsic dimension (usually unknown) of the data Inherent defects : sparsity, noise, outliers

Algorithmic Geometry Submanifold reconstruction J-D. Boissonnat 5 / 36

Issues in high-dimensional geometry

Dimensionality severely restricts our intuition and ability to visualize data ⇒ need for automated and provably correct methods methods Complexity of data structures and algorithms rapidly grow as the dimensionality increases ⇒ no subdivision of the ambient space is affordable ⇒ data structures and algorithms should be sensitive to the intrinsic dimension (usually unknown) of the data Inherent defects : sparsity, noise, outliers

Algorithmic Geometry Submanifold reconstruction J-D. Boissonnat 5 / 36

Issues in high-dimensional geometry

Dimensionality severely restricts our intuition and ability to visualize data ⇒ need for automated and provably correct methods methods Complexity of data structures and algorithms rapidly grow as the dimensionality increases ⇒ no subdivision of the ambient space is affordable ⇒ data structures and algorithms should be sensitive to the intrinsic dimension (usually unknown) of the data Inherent defects : sparsity, noise, outliers

Algorithmic Geometry Submanifold reconstruction J-D. Boissonnat 5 / 36

Looking for small and faithful simplicial complexes

Need to compromise Size of the complex

◮ can we have dim ˆ

M = dim M ?

Efficiency of the construction algorithms and of the representations

◮ can we avoid the exponential dependence on d ? ◮ can we minimize the number of simplices ?

Quality of the approximation

◮ Homotopy type & homology

(Cech and α complexes, persistence)

◮ Homeomorphism

(Delaunay-type complexes)

Algorithmic Geometry Submanifold reconstruction J-D. Boissonnat 6 / 36

Sampling and distance functions

[Niyogi et al.], [Chazal et al.]

Distance to a compact K : dK : x → infp∈K x − p

Stability

If the data points C are close (Hausdorff) to the geometric structure K, the topology and the geometry of the offsets Kr = d−1([0, r]) and Cr = d−1([0, r]) are close

Algorithmic Geometry Submanifold reconstruction J-D. Boissonnat 7 / 36

Distance functions and triangulations

ˇ

Nerve theorem (Leray)

The nerve of the balls (Cech complex) and the union of balls have the same homotopy type (same result for the α-complex)

Algorithmic Geometry Submanifold reconstruction J-D. Boissonnat 8 / 36

Questions

+ The homotopy type of a compact set X can be computed from the ˘ Cech complex of a sample of X + The same is true for the α-complex – The ˘ Cech and the α-complexes are huge (O(nd) and O(n⌈d/2⌉)) and difficult to compute in high dimensions – Both complexes are not in general homeomorphic to X (i.e. not a triangulation of X) – The ˘ Cech complex cannot be realized in general in the same space as X

Algorithmic Geometry Submanifold reconstruction J-D. Boissonnat 9 / 36

˘ Cech and Rips complexes

The Rips complex is easier to compute but still very big, and less precise in approximating the topology

α

Algorithmic Geometry Submanifold reconstruction J-D. Boissonnat 10 / 36

An example where no offset has the right topology !

• 1. Manifold + small noise assumption
• 2. Call persistent homology at rescue !

Algorithmic Geometry Submanifold reconstruction J-D. Boissonnat 11 / 36

The curses of Delaunay triangulations in higher dimensions

Complexity depends exponentially on the ambient dimension. Robustness issues become very tricky Higher dimensional Delaunay triangulations are not thick even if the vertices are well-spaced The restricted Delaunay triangulation is no longer a good approximation of the manifold even under strong sampling conditions (for d > 2)

Algorithmic Geometry Submanifold reconstruction J-D. Boissonnat 12 / 36

3D Delaunay Triangulations are not thick even if the vertices are well-spaced

Each square face can be circumscribed by an empty sphere This remains true if the grid points are slightly perturbed therefore creating thin simplices

Algorithmic Geometry Submanifold reconstruction J-D. Boissonnat 13 / 36

Badly-shaped simplices

Badly-shaped simplices lead to bad geometric approximations

which in turn may lead to topological defects in Del|M(P)

[Oudot]

see also [Cairns], [Whitehead], [Munkres], [Whitney]

Algorithmic Geometry Submanifold reconstruction J-D. Boissonnat 14 / 36

Tangent space approximation

M is a smooth k-dimensional manifold (k > 2) embedded in Rd

Bad news

[Oudot 2005]

The Delaunay triangulation restricted to M may be a bad approximation of the manifold even if the sample is dense

u v w p0 c0 t = ∆ x y z t p t = ∆ + δ/2 c x y z t

Algorithmic Geometry Submanifold reconstruction J-D. Boissonnat 15 / 36

Thickness and tangent space approximation

Lemma

[Whitney 1957]

If σ is a j-simplex whose vertices all lie within a distance η from a hyperplane H ⊂ Rd, then sin ∠(aff (σ), H) ≤ 2j η D(σ)

Corollary

If σ is a j-simplex, j ≤ k, vert (σ) ⊂ M, ∆(σ) ≤ δ rch(M) ∀p ∈ σ, sin ∠(aff(σ), Tp) ≤ δ Θ(σ)

(η ≤

∆(σ)2 2 rch(M) by the Chord Lemma) Algorithmic Geometry Submanifold reconstruction J-D. Boissonnat 16 / 36

The assumptions

M is a differentiable submanifold of positive reach of Rd The dimension k of M is small P is an ε-net of M, i.e.

◮

∀x ∈ M, ∃ p ∈ P, x − p ≤ ε rch(M)

◮ ∀p, q ∈ P, p − q ≥ ¯

η ε

ε is small enough

Algorithmic Geometry Submanifold reconstruction J-D. Boissonnat 17 / 36

The tangential Delaunay complex

[B. & Ghosh 2010]

p Tp M

1

Construct the star of p ∈ P in the Delaunay triangulation DelTp(P)

• f P restricted to Tp

2

DelTM(P) =

p∈P star(p)

Algorithmic Geometry Submanifold reconstruction J-D. Boissonnat 18 / 36

+ DelTM(P) ⊂ Del(P) + star(p), DelTp(P) and therefore DelTM(P) can be computed without computing Del(P) – DelTM(P) is not necessarily a triangulated manifold

Algorithmic Geometry Submanifold reconstruction J-D. Boissonnat 19 / 36

Construction of DelTp(P)

Given a d-flat H ⊂ R, Vor(P) ∩ H is a weighted Voronoi diagram in H

pi pj x p′

i

p′

j

H

x − pi2 ≤ x − pj2 ⇔ x − p′

i2 + pi − p′ i2 ≤ x − p′ j2 + pj − p′ j2

Corollary: construction of DelTp

Ψp(pi) = (p′

i, −pi − p′ i2)

(weighted point)

1

project P onto Tp which requires O(Dn) time

2

construct star(Ψp(pi)) in Del(Ψp(pi)) ⊂ Tpi

3

star(pi) ≈ star(Ψp(pi)) (isomorphic )

Algorithmic Geometry Submanifold reconstruction J-D. Boissonnat 20 / 36

Construction of DelTp(P)

Given a d-flat H ⊂ R, Vor(P) ∩ H is a weighted Voronoi diagram in H

pi pj x p′

i

p′

j

H

x − pi2 ≤ x − pj2 ⇔ x − p′

i2 + pi − p′ i2 ≤ x − p′ j2 + pj − p′ j2

Corollary: construction of DelTp

Ψp(pi) = (p′

i, −pi − p′ i2)

(weighted point)

1

project P onto Tp which requires O(Dn) time

2

construct star(Ψp(pi)) in Del(Ψp(pi)) ⊂ Tpi

3

star(pi) ≈ star(Ψp(pi)) (isomorphic )

Algorithmic Geometry Submanifold reconstruction J-D. Boissonnat 20 / 36

Inconsistencies in the tangential complex

A simplex is not in the star of all its vertices τ ∈ star(pi) ⇔ Tpi ∩ Vor(τ) = ∅ ⇔ B(cpi(τ) ∩ P = ∅ τ ∈ star(pj) ⇔ Tpj ∩ Vor(τ) = ∅ ⇔ B(cpj(τ) ∩ P ∋ p

Algorithmic Geometry Submanifold reconstruction J-D. Boissonnat 21 / 36

Inconsistency (k + 1)-trigger

pi pj τ Bpj(τ) Bpi(τ) p Tpi ∈ Vor(τ) ∈ aff(Vor(τ)) cpi(τ) Tpj cpj(τ) M iφ

Bpi(τ) : ball circumscribing τ centered on Tpi, cpi its center Inconsistency : Bpi(τ) ∩ P = ∅ and Bpj(τ) ∩ P = ∅ pl ∈ Bij, first point hit by (1 − λ)Bpi + λBpj, λ : 0 → 1 Trigger τ † : (k + 1)-simplex τ ⋆ pl ∈ Del((P))

Algorithmic Geometry Submanifold reconstruction J-D. Boissonnat 22 / 36

Inconsistency (k + 1)-triggers are flakes

pi pj τ Bpj(τ) Bpi(τ) p Tpi ∈ Vor(τ) ∈ aff(Vor(τ)) cpi(τ) Tpj cpj(τ) M iφ

If τ is small and thick, then Tpi ≈ Tpj ≈ aff(τ)

(sample density)

cpi − cpj small ⇒ Bij := Bpi(τ) \ Bpj(τ) = ∅ is small

(τ thick)

the trigger τ † = τ ⋆ pl is not thick

Algorithmic Geometry Submanifold reconstruction J-D. Boissonnat 23 / 36

Bound on the diameter of the simplices of DelTM(P)

(i) Vor(p) ∩ Tp ⊆ B(p, α0 rch(M)) where α0 ≈ ε (ii) ∀σ ∈ star(p), Rp(σ) ≤ α rch(M) (iii) ∀σ ∈ DelTM(P), ∆(σ) ≤ 2α rch(M).

Algorithmic Geometry Submanifold reconstruction J-D. Boissonnat 24 / 36

Proof of (i)

x′ x x′′ p

x ∈ Vor(p) ∩ Tp, p − x = α rch(M) x′ be the point of M closest to x and x′′ = Πp(x′) p − x′ ≤ p − x + x − x′ ≤ 2p − x ⇒ x′ − x′′ ≤ p−x′2

2rch(M) ≤ 2α2 rch(M) (Chord Lemma)

x′ − x′′ = x − x′ cos φ, where φ = ∠(Tx′, Tp) and cos φ ≥ 1 − 8α2 ⇒ x − x′ ≤ 2α2 rch(M) 1 − 8α2 assuming α ≤ √ 2 4 P is an ε-dense sample : ∃q ∈ P, x′ − q ≤ ε rch(M) x−p = α rch(M) ≤ x−q ≤ x−x′+x′−q ≤

• 2α2

1 − 8α2 + ε

• rch(M)

Algorithmic Geometry Submanifold reconstruction J-D. Boissonnat 25 / 36

Bound on the diameter of inconsistency triggers

τ an inconsistent k-simplex, τ † a trigger, θ = maxp∈τ ∠(aff(τ), Tp) Lemma sin θ ≤

∆(τ) Θ(τ)) rch(M)

and R(τ †) ≤ R(τ)

cos θ

Tpi cpi c(τ) pi τ R(τ) ω pl

Proof d(pl, Tp) ≤

∆2(τ) 2rch(M) (Chord Lemma)

sin ∠(aff(τ), Tp) ≤

2 ∆2(τ)

2rch(M)

Θ(τ) ∆(τ) = ∆(τ) Θ(τ) rch(M) (Whitney’s angle bound)

Rpi(τ) = pi − cpi ≤ R(τ)

cos θ

and R(τ †) ≤ i(τ †) − pi ≤ R(τ)

cos θ

Algorithmic Geometry Submanifold reconstruction J-D. Boissonnat 26 / 36

Bound on the thickness of inconsistency triggers

Lemma Θ(τ †) ≤

∆(τ †) 2(k+1)rch(M)

• 1 +

2 Θ(τ)

• Tpi

cpi c(τ) pi τ R(τ) ω pl

Proof Let q ∈ τ. D(pl, τ †) = pl − q sin ∠(pl − q, aff(τ)) ≤ ∆(τ †) (sin ∠(pl − q, Tq) + sin ∠(Tq, aff(τ))) ≤ ∆(τ †) ∆(τ †) 2rch(M) + ∆(τ) Θ(τ) rch(M)

• (Chord + previous Lemmas)

≤ ∆2(τ †) 2rch(M)

• 1 +

2 Θ(τ)

• Hence, if τ is thick, τ † cannot be so : we say that τ † is a flake

Algorithmic Geometry Submanifold reconstruction J-D. Boissonnat 27 / 36

Reconstruction of smooth submanifolds

1

For each vertex v, compute the star star(p) of p in Delp(P)

2

Remove inconsistencies among the stars by weighting the points

3

Stitch the stars to obtain a triangulation of P

v

Algorithmic Geometry Submanifold reconstruction J-D. Boissonnat 28 / 36

Algorithm hypotheses

Known quantities in red M = a differentiable submanifold of positive reach of dim. k ⊂ Rd P = an (ε, δ)-sample of M ε ≤ ε0 ε/δ ≤ η0 we can estimate the tangent space Tp at any p ∈ P

Algorithmic Geometry Submanifold reconstruction J-D. Boissonnat 29 / 36

Removing inconsistencies by removing flakes

Another application of Moser-Tardos algorithmic LLL

Input: P, {Tp, p ∈ P}, ˜ w0, Θ0 Initialize all weights to 0 and compute DelTM( ˆ P) while there are Θ0-flakes or inconsistencies in DelTM( ˆ P) do while there is a Θ0-flake σ in DelTM( ˆ P) do resample σ, i.e. reweight the vertices of σ update DelTM( ˆ P) if there is an inconsistent simplex σ in DelTM( ˆ P) then compute a trigger simplex σ† associated to σ resample the flake σ ⊂ σ† update DelTM( ˆ P) Output: A weighting scheme on P and DelTM( ˆ P) DelTM( ˆ P) is Θ0-thick and has no inconsistency

Algorithmic Geometry Submanifold reconstruction J-D. Boissonnat 30 / 36

Summary

Termination

◮ If ¯

η 2 ≥ ¯

ρ ≥ f(Θ0), the algorithm terminates and returns a complex ˆ M that has no inconsistent configurations

Complexity

◮ No d-dimensional data structure ⇒ linear in d ◮ exponential in k

Approximation

◮ ˆ

M is a PL simplicial k-manifold

◮ ˆ

M ⊂ tub(M, ε)

◮ ˆ

M is homeomorphic to M

Algorithmic Geometry Submanifold reconstruction J-D. Boissonnat 31 / 36

ˆ M is a PL simplicial k-manifold

Lemma Let P be an ε-sample of a manifold M and let p ∈ P. The link of any vertex p in ˆ M is a topological (k − 1)-sphere Proof :

• 1. Since ˆ

M contains no inconsistencies, ∀p ∈ vert( ˆ M), star(p, ˆ M) = star(p, Delp(P))

• 2. Delp(P) ⊂ Rd ≈ Del(Ψp(P)) ⊂ Tp

⇒ star(p) ≈ starp(p)

• 3. starp(p) is a k-dimensional triangulated topological ball (general position)
• 4. p cannot belong to the boundary of starp(p)

(the Voronoi cell of p = Ψp(p) in Vor(Ψp(P)) is bounded)

Algorithmic Geometry Submanifold reconstruction J-D. Boissonnat 32 / 36

ˆ M is a triangulation of M

Theorem M : a connected compact k-submanifold of Rd without boundary ˆ M : combinatorial k-manifold without boundary, embedded in Rd s.t.

1

P = vert( ˆ M) ⊂ M

2

∀σ ∈ ˆ M, ∆(σ) ≤ δ0 rch(M) where δ0 ≤ 1

5 (longest edge)

L(σ) ≥ λ0 rch(M)

(shortest edge)

Θ(σ) ≥ 8.13 δ2

0/λ0 (thickness)

3

∃pi ∈ ˆ M s.t. Π−1(pi) = {pi} Then

1

ˆ M is a triangulation of M

2

The Hausdorff distance between ˆ M and M is at most 2δ2

0 rch(M)

3

If σ is a k-simplex of ˆ M and p one of its vertices, we have sin ∠(aff(σ), Tp) ≤ δ0 Θ0

Algorithmic Geometry Submanifold reconstruction J-D. Boissonnat 33 / 36

Triangulating a Riemann surface

a complex curve in the projective plane (parameterized as a real surface of R8)

Algorithmic Geometry Submanifold reconstruction J-D. Boissonnat 34 / 36

Triangulating the conformational space of C8H16

Algorithmic Geometry Submanifold reconstruction J-D. Boissonnat 35 / 36

Applications and extensions

Discrete metric sets (see the previous lecture on the witness complex) Anisotropic mesh generation Non euclidean embedding space (e.g. statistical manifolds)

Algorithmic Geometry Submanifold reconstruction J-D. Boissonnat 36 / 36