Gabriel Peyré Ceremade
email: gabriel.peyre@ceremade.dauphine.fr

Numerical Mesh Processing

These course notes are an introduction to the geometry of surfaces both from a local (differential), global (geodesic) and multiscale point of view. The emphasis is put on the discrete computation of differential and geodesic quantities on 3D meshes. These computational tools can be used to perform various mesh processing tasks such as surface denoising, compression or recognition. All the numerical experiments shown in these notes can be reproduced with a set of Matlab scripts available in the courses section of my webpage.

This chapter exposes the basics of surface approximation with 3D meshes and the way to process such meshes with linear operators. In particular, it studies filtering on 3D meshes and explains how a Fourier theory can be built to analyze these filters.

Surface Discretization with Triangulated Mesh 1.Continuous Geometry of Surfaces

In this course, in order to simplify the mathematical description of surfaces, we consider only globally parameterized surfaces. We begin by considering surfaces embedded in euclidean space M ⊂ R k .

Definition 1 (Parameterized surface). A parameterized surface is a mapping

u ∈ D ⊂ R 2 → ϕ(u) ∈ M.
Of course, most surfaces do not benefit from such a simple parameterization. For instance, a sphere should be split into two parts in order to be mapped on two disks D 1 , D 2 . These topological difficulties require the machinery of manifolds in order to incorporate a set of charts D = {D i } i that overlap in a smooth manner. All the explanations of this course extend seamlessly to this multi-charts setting.

A curve is defined in parameter domain as a 1D mapping t ∈ [0, 1] → γ(t) ∈ D. This curve can be traced over the surface and its geometric realization is γ(t) def.

= ϕ(γ(t)) ∈ M. The computation of the length of γ in ambient k-dimensional space R k follows the usual definition, but to do the computation over the parametric domain, one needs to use a local metric defined as follow.

Definition 2 (First fundamental form). For an embedded manifold M ⊂ R k , the first fundamental form is

I ϕ = ∂ϕ ∂u i , ∂ϕ ∂u j i,j=1,2 .
This local metric I ϕ defines at each point the infinitesimal length of a curve as

L(γ) def. = 1 0 ||γ ′ (t)||dt = 1 0 γ ′ (t) T I ϕ (γ(t))γ ′ (t)dt.
This fundamental form is an intrinsic invariant that does not depends on how the surfaces is isometrically embedded in space (since the length depends only on this tensor field I ϕ ). In contrast, higher order differential quantities such as curvature might depend on the bending of the surface and are thus usually not intrinsic (with the notable exception of invariants such as the gaussian curvature). In this course, we restrict ourselves to first order quantities since we are mostly interested in lengths and the intrinsic study of surfaces.

Example 1 (Isometry and conformality). A surface M is locally isometric to the plane if I ϕ = Id 2 . This is for instance the case for a cylinder. The mapping ϕ is said to be conformal if I ϕ (u) = λ(u)Id 2 . It means that the length of a curve over the plane is only locally scaled when mapped to the surface. In particular, the angle of two interescting curves is the same over the parametric domain and over the surface. This is for instance the case for the stereographic mapping between the plane and a sphere.

Discretization of Surfaces with Triangulations

Mesh Data Structure A triangulated mesh is a discrete structure that can be used to approximate a surface embedded in Euclidean space R k . It is composed of a topological part M = (V, E, F ) and a geometrical realization M = (V, E, F). It is important to make the distinction between these two parts since many algorithms rely only on geometry (point clouds processings such as dimension reduction) or on topology (such as compression).

The topology M of the mesh is composed of Vertices (0D): this is an abstract set of indices V ≃ {1, . . . , n}.

Edges (1D): this is a set of pair of vertices E ⊂ V × V . This set is assumed to be symmetric

(i, j) ∈ E ⇐⇒ i ∼ j ⇔ (j, i) ∈ E.
Faces (2D): this is a collection of 3-tuples of vertices F ⊂ V × V × V , with the additional compatibility condition (i, j, k) ∈ F =⇒ (i, j), (j, k), (k, i) ∈ E.

We further assumer that there is no isolated edges

∀ (i, j) ∈ E, ∃ k, (i, j, k) ∈ F.
The set of edges can be stored in a symmetric matrix A ∈ R n×n such that A ij = 1 if (i, j) ∈ E and A ij = 0 otherwise. This matrix is often stored as a sparse matrix since the number of edges is usually much smaller than n 2 . The set of vertices and edges form a non-oriented graph G = (V, E). Faces are often stored as a matrix A F ∈ {1, . . . , n} 3×m where m is the number of faces and a column ((A F ) i,1 , (A F ) i,2 , (A F ) i,3 ) stores the indices of a face. In a triangulation, the face matrix A F allows to recover the edge incidence matrix A. The face data structure allows to really capture the 2D geometry of surfaces, which is not possible with graphs alone.

The geometric realization M is defined through a spacial localization of the vertices (for instance in 3D space) V def.

= {x i \ i ∈ V } ⊂ R 3 .
This allows to define a piecewise linear mesh

F def. = (i,j,k)∈F Conv(x i , x j , x k ) ⊂ R 3 ,
where the convex envelop Conv(x, y, z) of three points is the Euclidean triangle generated by (x, y, z).

This piecewise linear realization M can be displayed as a 3D surface on a computer screen. This is performed through a perspective projection of the points and a linear interpolation of color and light inside the triangle. Figure 1.1 shows an example of 3D display, with a zoom on the faces of the mesh.

Adjacency Relationships From the basis topological information given by M = (V, E, F ), one can deduce several adjacency data-structures that are important to navigate over the triangulation.

Definition 3 (Vertex 1-ring). The vertex 1-ring of a vertex i ∈ V is

V i def.
= {j ∈ V \ (i, j) ∈ E} ⊂ V. ' (1.1)

The s-ring is defined by induction as

∀ s > 1, V (s) i
= j ∈ V \ (k, j) ∈ E and k ∈ V (s-1) i

.

(1.2) Definition 4 (Face 1-ring). The face 1-ring of a vertex i ∈ V is

F i def. = {(i, j, k) ∈ F \ i, j ∈ V } ⊂ F.
The geometrical realization of a vertex 1-ring is

V i = (i,j,k)∈Vi
Conv(x i , x j , x k ).

A triangulated mesh is a manifold mesh if all the rings V i for i ∈ V are homeomorphic to either a disk (for interior vertices) or to a half disk (for boundary vertices). This ensures that the geometrical mesh really has the topology of a 2D surface embedded in R 3 (possibly with boundaries). In particular, it implies that there is at most two faces connected to each edge

∀ (i, j) ∈ E, # {k \ (i, j, k) ∈ F } 2.
As an application of these local rings, one can compute a normal at each point using a simple rule

∀ f = (i, j, k) ∈ F, -→ n f def. = (x j -x i ) ∧ (x k -x i ) ||(x j -x i ) ∧ (x k -x i )|| .
and where

∀ i ∈ V, -→ n i def. = f ∈Fi -→ n f || f ∈Fi -→ n f || .
These normals are used to define for instance a light intensity I(i) = max( n i , ℓ(i), , 0), where ℓ(i) is the incident light. In practice one uses a infinite light source ℓ(i) = ℓ =constant or a local spot located at position s ∈ R 3 through ℓ(i) = (v i -s)/||v i -s||. This light intensity is interpolated on the whole mesh during display.

Linear Mesh Processing

The light intensity I is a particular example of a function defined at each vertex of the mesh. Mesh processing is intended to process such functions and we thus define carefully vector spaces and operators on meshes.

Functions on a Mesh

In this course, a function is a discrete set of values defined at each vertex location.

Definition 5 (Linear space on a mesh). A function on a mesh is a mapping f ∈ ℓ 2 (V) ≃ ℓ 2 (V ) ≃ R n and can be viewed equivalently as

f : V -→ R x i -→ f (x i ) ⇐⇒ f : V -→ R i -→ f i ⇐⇒ f = (f i ) i∈V ∈ R n .
The linear space of the functions on a mesh is equipped with an Hilbert space structure that allows to quantify approximation error and compute projections of functions.

Definition 6 (Inner product and norm). One defines the following inner product and norm for vector f, g ∈ R n f, g def.

=

i∈V f i g i and ||f || 2 = f, f .

In order to modify (process) functions on a mesh (such as a light intensity I), this course considers only linear operations that are defined through a large matrix.

Definition 7 (Linear operator A). A linear operator A is defined as

A : ℓ 2 (V ) → ℓ 2 (V ) ⇐⇒ A = (a ij ) i,j∈V ∈ R n×n (matrix).
and operate on a function f as follow

(Af )(x i ) = j∈V a ij f (x j ) ⇐⇒ (Af ) i = j∈V a ij f j .
Example 2. If the coordinates of the point of a mesh are written

x i = (x 1 i , x 2 i , x 3 i ) ∈ R 3 , then the X-coordinate defines a function f : i ∈ V → x 1 i ∈ R.
A geometric mesh M is thus 3 functions defined on M . Mesh processing is the task of modifying functions f ∈ ℓ 2 (V ). For instance, one can denoise a mesh M as 3 functions on M . The usual strategy applies a linear operator f → Af . Sometimes, A can computed from M only (for instance for compression) but most of the times it requires both M and M.

Local Operators

In most applications, one can not store and manipulate a full matrix A ∈ R n×n . Furthermore, one is usually interested in exploiting the local redundancies that exist in most usual functions f ∈ R n defined on a mesh. This is why we restrict our attention to local operators that can be conveniently stored as sparse matrices (the zeros are not kept in memory).

Definition 8 (Local operator). A local operator

W ∈ R n×n satisfies w ij = 0 if (i, j) / ∈ E. (W f ) i = (i,j)∈E w ij f j .
A particularly important class of local operators are local smoothings (also called filterings) that perform a local weighted sum around each vertex of the mesh. For this averaging to be consistent, we define a normalized operator W whose set of weights sum to one.

Definition 9 (Local averaging operator).

A local normalized averaging is W = ( wij ) i,j∈V 0 where ∀ (i, j) ∈ E, wij = w ij (i,j)∈E w ij .

It can be equivalently expressed in matrix form as

W = D -1 W with D = diag i (d i ) where d i = (i,j)∈E w ij .
The smoothing property corresponds to W 1 = 1 which means that the unit vector is an eigenvector of W with eigenvalue 1.

Example 3. In practice, we use three popular kinds of averaging operators.

Combinatorial weights: they depends only on the topology (V, E) of the vertex graph

∀ (i, j) ∈ E, w ij = 1.
Distance weights: they depends both on the geometry and the topology of the mesh, but do not require faces information,

∀ (i, j) ∈ E, w ij = 1 ||x j -x i || 2 .
Conformal weights: they depends on the full geometrical realization of the 3D mesh since they require the face information

∀ (i, j) ∈ E, w ij = cot(α ij ) + cot(β ij ).
(1.3) Figure 1.2 shows the geometrical meaning of the angles α ij and β ij α ij = ∠(x i , x j , x k1 ) and

β ij = ∠(x i , x j , x k2 ),
where (i, j, k 1 ) ∈ F and (i, j, k 2 ) ∈ F are the two faces adjacent to edge (i, j) ∈ E. We will see in the next section the explanation of these celebrated cotangent weights.

x i One can use iteratively a smoothing in order to further filter a function on a mesh. The resulting vectors W f, W 2 , . . . , W k f are increasingly smoothed version of f . Figure 1.3 shows an example of such iterations applied to the three coordinates of mesh. The sharp features of the mesh tend to disappear during iterations. We will make this statement more precise in the following, by studying the convergence of these iterations.

x j x k 1 x k 2 α ij β ij

Approximating Integrals on a Mesh

Before investigating algebraically the properties of smoothing operators, one should be careful about what are these discrete operators really approximating. In order for the derivation to be simple, we make computation for a planar triangulation M of a mesh M ⊂ R 2 .

In the continuous domain, filtering is defined through integration of functions over the mesh. In order to descretize integrals, one needs to define a partition of the plane into small cells centered around a vertex or an edge. These Voronoi cells indeed form a partition of the mesh

Definition 10 (Vertices Voronoi). The Voronoi diagram associated to the vertices is

∀ i ∈ V, E i = {x ∈ M \ ∀ j = i, ||x -x i || ||x -x j ||}
M = i∈V E i = e∈E E e .
The following theorem gives the formula for the area of these cells.

Theorem 1 (Voronoi area formulas). For all e = (i, j) ∈ E, ∀ i ∈ V , one has

A e = Area(E e ) = 1 2 ||x i -x j || 2 (cot(α ij ) + cot(β ij )) A i = Area(E i ) = 1 2 j∈Ni A (ij) .
With these areas, one can approximate integrals on vertices and edges using

M f (x)dx ≈ i∈V A i f (x i ) ≈ e=(i,j)∈E A e f ([x i , x j ]).
Of particular interest is the approximation of the so-called Dirichelet energy M ||∇ x f || 2 dx. In order to compute it on a triangular mesh, one can use a finite difference approximation of the gradient of a function at the point x ij = (x i + x j )/2 along an edge (i, j)

∇ xij f, x i -x j ||x i -x j || ≈ f (x i ) -f (x j ) ||x i -x j || .
This leads to the following approximation of the Dirichlet energy

M ||∇ x f || 2 dx ≈ (i,j)∈E A (i,j) ∇ xij f, x i -x j ||x i -x j || 2 ≈ (i,j)∈E A (i,j) |f (x j ) -f (x i )| 2 ||x j -x i || 2 (1.4) = (i,j)∈E w ij |f (x j ) -f (x i )| 2 where w ij = cot(α ij ) + cot(β ij ). (1.5)
This discrete formulation shows that the correct weights to approximate the Dirichlet energy are the cotangent one, already introduced in equation (1.3).

Example on a Regular Grid

A regular grid is an uniform discretization with n points of [0, 1) (in 1D) or [0, 1) 2 (in 2D). One usually assumes periodic boundary conditions, which means that each side of the square is associated with its opposite.

Since the geometry of a regular grid is invariant under translation, local averaging operators can be computed as convolution on D = (Z/pZ) d where n = p d for d the dimension of the domain

(d = 1 or d = 2) ∀ i ∈ D, W f (i) = k∈D f (k) w(i -k),
where the operation + and -should be computed modulo p and w(k) = W (0, k) is the convolution kernel.

Example 4 (Averaging). The uniform averaging filter is defined as

W f (i) = 1 |N | k∈N f (i + k),
where N is the set of neighbors of the point 0 and |N | = 2 d . In this case, in dimension 1, w = (1, 0, 1)/2, where this notation assumes that w is centered at the point 0.

In order to study translation invariant operators like local filtering, one needs to use the discrete Fourier transform that diagonalizes these operators.

Definition 12 (Discrete Fourier transform). The 1D discrete Fourier transform Φ(f ) ∈ C n of the vector f ∈ C n Φ(f )(ω) = f (ω) def. = 1 n k f k e 2iπ n kω .
A similar definition can be given for the 2D discrete Fourier transform. The main property of the Fourier transform is the following diagonalization result.

Theorem 2 (Convolution and Fourier). For any vector f , one has

Φ( W k f ) = Φ( w * . . . * w * f ) =⇒ Φ( W k f )(ω) = w(ω) k f (ω).
The main interest of this tools is that Φ(f ) can be computed in O(n log(n)) operations with the FFT algorithm. Using the following theorem, it gives an alternative expression of a local filtering. This expression in the Fourier domain can be used to speed up the computation of w * f if w has a lot of non zero entries (which is not the case in our setting of local operators). It is also useful to analyze theoretically the behavior of iterated filterings.

Theorem 3 (Convergence). For any function f defined on a regular grid in 1D or 2D, one has

W k f k→+∞ -→ 1 |V | i∈V f i
This Fourier theory can only be developed for domains that have a group structure that enables translation invariant filtering. In particular, it does not carry over easily to an arbitrary surface. In the remaining, we define a corresponding theory for graphs and triangulated surfaces using the eigenvector of Laplacian operators. This Fourier transform on meshes enables the analysis of the convergence of many filtering schemes.

Gradients and Laplacians on Meshes

Gradient operator A gradient operator defines directional derivatives on a triangulation. It maps functions defined on vertices to functions defined on the set of oriented edges

Ē def. = {(i, j) ∈ E \ i > j} .
Definition 13 (Gradient). Given a local averaging W , the gradient operator G is defined as

∀ (i, j) ∈ E, i < j, (Gf ) (i,j) def. = √ w ij (f j -f i ) ∈ R.
This mapping can be viewed equivalently as

G : ℓ 2 (V ) -→ ℓ 2 (E), or G : R n -→ R p where p = |E|, or G ∈ R n×p (a matrix).
The value of (Gf ) e for an edge e = (i, j) can be thought as a derivative along direction --→ x i x j .

Example 5. For the local averaging based on square distances, one has

w ij = ||x i -x j || -2 , (Gf ) (i,j) = f (x j ) -f (x i ) ||x i -x j || .
which is exactly the finite difference discretization of a directional derivative.

One a regular grid, one can note that

Gf discretizes ∇f = ∂f ∂x , ∂f ∂y T . G T v discretizes div(v) = ∂v1 ∂x + ∂v2 ∂y .
Laplacian Operator A Laplacian operator is a discrete version of a second order derivative operator.

Definition 14 (Laplacian). Given a local averaging W , the Laplacian operator D is defined as

L def. = D -W, where D = diag i (d i ), with d i = j w ij .
In the remaining, we also make use of normalized operators, which have an unit diagonal.

Definition 15 (Normalized Laplacian). The normalized Laplacian is defined as

L def. = D -1/2 LD -1/2 = Id n -D -1/2 W D 1/2 = Id n -D 1/2 W D -1/2 .
This normalized Laplacian correspond to the weighted graph Laplacian used in graph theory, see for instance [START_REF] Chung | Spectral graph theory[END_REF].

Linear Mesh Processing

Remark 1. One can note that Laplacians are symmetric operators L, L ∈ R n×n . L acts like a (second order) derivative since L1 = 0. in contrast, the normalized Laplacian is not a real derivative since L1 = 0 in general.

The main interest of the gradient operator is that it factorizes the Laplacian as follow.

Theorem 4 (Laplacian factorization). One has

L = G T G and L = (GD -1/2 ) T (GD -1/2 ).
This theorem proves in particular that L and L are symmetric positive definite operators. The inner product defined by the Laplacian can be expressed as an energy summed over all the edges of the mesh

Lf, f = ||Gf || 2 = (i,j)∈E w ij ||f i -f j || 2 .
In the particular case of the cotangent weights introduced in equation ( 1.3), one can see that the Laplacian norm Lf, f is exactly the finite differences approximation of the continuous Dirichlet energy M |∇ x f |dx derived in equation (1.5). This is why these cotangent weights are the best choice to compute a Laplacian that truly approximates the continuous Laplace Beltrami operator (see definition [START_REF] Khodakovsky | Progressive geometry compression[END_REF].

A similar expression is derived for the normalized laplacian

Lf, f = ||GD -1/2 f || 2 = (i,j)∈E w ij f i √ d i - f j d j 2 .
Of particular interest for the study of filtering on meshes is the behavior of the spectrum of the Laplacian. We can first study its kernel.

Theorem 5 (Kernel of the Laplacian). If M is connected, then ker(L) = span(1) and ker(L) = span(D 1/2 ).

Examples in 1D and 2D

In 1D, all local weights are equivalent since the points are equi-spaced. The corresponding Laplacian is a convolution that can be written as

(Lf ) i = 1 h 2 (2f i -f i+1 -f i-1 ) = 1 h 2 f * (-1, 2, 1) ,
where it is important to remember that the notation (-1, 2, 1) means that the vector is centered around 0.

This discrete 1D Laplacian is the finite difference approximation of the continuous Laplacian on the torus T of the segment [0, 1) modulo 1. Up to a minus sign, this Laplacian is just the second order derivative

L h→0 -→ - d 2 f dx 2 (x i )
One should be careful with our notation that consider positive semi-definite Laplacian, that have the opposite sign with respect to second order derivative operators (which are definite negative).

The gradient operator corresponds to a discretization of the first order derivative f → f ′ (which is anti symmetric). The continuous counterpart of the factorization L = G T G is the integration by part formula on the torus

T f ′′ (x)g(x)dx = - T -f (x)g ′ (x)dx =⇒ T f ′′ (x)f (x)dx = - T |f ′ (x)| 2 0.
The discrete Laplacian on a 2D grid can also be written as a 2D convolution

(Lf ) i = 1 h 2 (4f i -f j1 -f j2 -f j3 -f j4 ) = 1 h 2 f *   0 -1 0 -1 4 -1 0 -1 0  
where {j k } k are the four neighbors of the point i. This operator is the finite difference approximation to the continuous 2D Laplacian

L h→0 -→ - ∂ 2 f ∂x 2 (x i ) - ∂ 2 f ∂y 2 (x i ) = -∆f (x i ).
The factorization Lf = G T Gf corresponds to the decomposition ∆f = div(∇f ).

Example of a Parametric Surface

We recall that a parameterized surface is a mapping u ∈ D ⊂ R 2 → ϕ(u) ∈ M. Whereas the continuous Laplacian is simple to define on the plane using partial derivatives, its definition on a surface requires the intervention of an arbitrary parameterization ϕ which makes its expression cumbersome.

Definition 16 (Laplace-Beltrami). The Laplace-Beltrami operator on a parametric surface M is defined as

√ g∆ M def. = ∂ ∂u 1 g 22 √ g ∂ ∂u 1 - g 12 √ g ∂ ∂u 2 + ∂ ∂u 2 g 11 √ g ∂ ∂u 2 - g 12 √ g ∂ ∂u 1
where g = det(I ϕ ) and

I ϕ = (g ij ) i,j=1,2 .
The Laplacian is however an intrinsic operator that does not depends on the chosen parameterization, as shown by the following approximation theorem.

Remark 2 (Laplacian using averaging).

∆ M f (x) = lim h→0 1 |B h (x)| y∈M f (y)dy where B h (x) = {y \ d M (x, y) h}
where d M is the geodesic distance on M and h = max (i,j)∈E ||x i -x j || is the discretization precision.

Diffusion and Regularization on Surfaces 1.3.1 Heat Diffusion

The main linear PDE for regularization of functions is the heat equation that governs the isotropic diffusion of the values of a function in time.

Definition 17 (Heat diffusion). ∀ t > 0, one defines

F t : M → R solving ∂F t ∂t = -D -1 LF t = -(Id n -W )F t and ∀ i ∈ V, F 0 (i) = f (i)
In order to compute numerically the solution of this PDE, one can fix a time step δ > 0 and use an explicit discretization in time Fk as F 0 = f and

1 δ Fk+1 -Fk = -D -1 L Fk =⇒ Fk+1 = Fk -δD -1 L Fk = (Id -δ) Fk + δ W Fk . (1.6)
If δ is small enough, one hopes that the discrete solution Fk is close to the continuous time solution F t for t = δk. This is indeed the case as proven later in these notes.

Diffusion and Regularization on Surfaces

Remark 3. In order for this scheme to be stable, one needs δ < 1. This is be proven later using the extension of Fourier theory to meshes.

Remark 4. If δ = 1, then the discretization of the Heat equation corresponds to iterative smoothing since Fk = W k f . In this case stability is not guaranteed but only pathological meshes give unstable filtering (see theorem 13).

Instead of using the explicit discretization in time (1.6), one can use an implicit scheme which compute an approximate solution Fk at step k by solving

1 δ Fk+1 -Fk = -D -1 L Fk+1 =⇒ ((δ + 1)Id n -δ W ) Fk+1 = Fk . (1.7)
Computing Fk requires the solution of a sparse linear system at each step k. The implicit scheme (1.7) is thus computationally more involved than the explicit scheme (1.6). We will however see later that the implicit scheme is always stable for any value of δ 1.

Example 6 (Mesh smoothing). In order to smooth a mesh whose points are In practice, mesh smoothing is used to denoise a function f = f 0 + σg where g ∈ R n is a realization of a gaussian white noise (each entry g(i) are independent and follow a gaussian law with unit variance). The difficult task it to find an optimal stopping time t to minimize ||F t -f 0 ||, which is not available since one does not know f 0 . For uniformly smooth surfaces, the theory predicts that a linear filtering such as the heat equation requires a stopping time proportional to the noise level σ. This is however false for more complex surfaces such as the one used in computer graphics. In these case, alternate non linear diffusions such as non-linear PDE or wavelet thresholding usually perform better, see [START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF] for an overview of these methods in image processing. 1 Linear Mesh Processing Other differential equations. One can solve other partial differential equations involving the Laplacian over a 3D mesh M = (V, E, F ). For instance, one can consider the wave equation, which defines, for all t > 0, a vector F t ∈ ℓ 2 (V ) as the solution of

x i = (x 1 i , x 2 i , x 3 i ), one can perform a heat diffusion for each component f i = (x k i ), k = 1, 2, 3.
∂ 2 F t ∂t 2 = -D -1 LF t and F 0 = f ∈ R n , d dt F 0 = g ∈ R n , (1.8)
In order to compute numerically the solution of this PDE, one can fix a time step δ > 0 and use an explicit discretization in time Fk as F 0 = f , F 1 = F 0 + δg and for k > 1

1 δ 2 Fk+1 + Fk-1 -2 Fk = -D -1 L Fk =⇒ Fk+1 = 2 Fk -Fk-1 -δ 2 D -1 L Fk .

Spectral Decomposition

In order to better understand the behavior of linear smoothing on meshes, one needs to study the spectral content of Laplacian operators. This leads to the definition of a Fourier theory for meshes. The decomposition L = (GD -1/2 ) T (GD -1/2 ) of the Laplacian implies that it is a positive semi-definite operator. One can thus introduce the following orthogonal factorization.

Theorem 6 (Eigen-decomposition of the Laplacian). It exists a matrix U, U T U = Id n such that

L = U ΛU T where Λ = diag ω (λ ω ), λ 1 . . . λ n .
The eigenvalues λ ω correspond to a frequency index that ranks the eigenvectors u ω of U = (u ω ) ω . One can first state some bounds on these eigenvalues.

Theorem 7 (Spectral bounds). ∀ i, λ i ∈ [0, 2] and If M is connected then 0 = λ 1 < λ 2 . λ n = 2 if and only if M is 2-colorable.
We recall the definition of a colorable graph next.

Definition 18 (Colorable graph

). A graph (V, E) is k-colorable if it exist a mapping f : V → {1, . . . , k} such that ∀ (i, j) ∈ E, f (i) = f (j).
A 2-colorable graph is also called bi-partite. A 2-colorable mesh is pathological for filtering since one can split the set of vertices into two parts without inner connexions. The filtering process can oscillate by exchanging values between these sets, thus never converging.

The orthogonal eigen-basis U = (u ω ) ω is an orthogonal basis of the space R n ≃ ℓ 2 (V ), which can be written as

u ω : V -→ R i -→ u ω (x i )
The orthogonality means that u ω , u ω ′ = δ ω ′ ω . This basis allows to compute an orthogonal decomposition of any functions f

∀ f ∈ ℓ 2 (V ), f = ω f, u ω u ω .
Having such a tool allows to split a function f in elementary contributions f, u ω with a control in the energy because of orthogonality shows some examples of eigenfunctions depicted using color ranging from blue (negative values of the eigenfunction) to red (positive values). One can see that these functions are oscillating, in a way similar to the traditional Fourier basis. In some sense (made more precise latter), this basis is the extension of the Fourier basis to meshes. A function u ω corresponding to a large spectral value λ ω is highly oscillating and corresponds thus intuitively to a high frequency atom.

||f || 2 = ω | f, u ω | 2 . ω = 2 ω = 4 ω = 8 ω = 12 ω = 16
Extracting numerically eigenvectors from a large matrix is a difficult problem. If the matrix is sparse, a method of choice consists in using iterative powers of a shifted version of the laplacian. One starts from a random initial vector v 0 and iterates

v k+1 = w k+1 ||w k+1 || where w k+1 = ( L -λId n ) -1 v k .
(1.9)

These iterates converges to the eigenvectors corresponding to the eigenvalue the closest to λ, as staten in the following theorem.

Theorem 8 (Inverse iterations). For a given shift λ, lets denote

ω ⋆ = argmin ω |λ -λ ω | and ω + = argmin ω =ω ⋆ |λ -λ ω | If |λ -λ ω ⋆ | < |λ -λ ω + |, then v k k→+∞ -→ u ω ⋆ and Lv k , v k k→+∞ -→ λ ω ⋆ .
The speed of convergence of these inverse iterations is governed by the conditioning of ( L -

λId n ) -1 since ||v k -u ω ⋆ || Cρ(λ) k where ρ(λ) def. = |λ -λ ω ⋆ | |λ -λ ω + | < 1.
The smallest ρ(λ) is, the faster the method converges.

In order to compute an iteration (1.9) of the method, one needs to solve a sparse linear system Aw k+1 = v k whith A = L -λId n . In order to do so, one can use a direct method such as LU factorization. The advantage of such an approach is that the factorization is computed once for all and can be re-used to solve very quickly at each step k. These factorization are however quite slow to compute especially for large matrices. For large problems, one can solve this linear system using an iterative algorithm such as conjugate gradient. These iterative method are attractive for sparse matrices, but a fast convergence requires 1/ρ(λ), the conditioning of L -λId n to be not large, with is contradictory with the constraint for iterations 1.9 to converge fast.

Spectral Theory on a Regular Grid

In the particular case of a 1D or 2D lattice, the eigenfunctions defined earlier correspond exactly to the Fourier basis used in the discrete Fourier transform.

Theorem 9 (Spectrum in 1D). For a 1D regular lattice,

u ω (k) = 1 √ n exp 2iπ n kω and λ ω = 4 sin 2 2π n ω .
Theorem 10 (Spectrum in 2D). For a 2D regular lattice, n = n 1 n 2 , ω = (ω 1 , ω 2 )

u ω (k) = 1 √ n exp 2iπ n k, ω and λ ω = 4 sin 2 2π n 1 ω 1 + sin 2 2π n 2 ω 2 .
As already mentioned, on a mesh, the eigenvectors of L correspond to a extension of the Fourier basis to meshes. The definition of the Fourier transform on meshes requires a little care since a diagonal normalization by D is used as defined next.

Definition 19 (Manifold-Fourier transform). For f ∈ ℓ 2 (V ), Φ(f )(ω) = f (ω) def. = D 1/2 f, u ω ⇐⇒ Φ(f ) = f = U T D 1/2 .
where (u ω ) ω are the eigenvectors of L.

One can note that there is still a degree of freedom in designing this Fourier transform since one can use any local weighting (for instance combinatorial, distance or conformal). Depending on the application, one might need to use weights depending only on the topology of the mesh (combinatorial for mesh compression).

A major theoretical interest of this Fourier transform is that it diagonalizes local averaging operators.

Diffusion and Regularization on Surfaces

Theorem 11 (Spectral smoothing). One has Φ W Φ -1 = Id n -Λ and thus for any function f

W f (ω) = (1 -λ ω ) f (ω)
This diagonalization allows to prove the convergence of iterative smoothing.

Theorem 12 (Convergence of iterated smoothing). If λ n < 2 (i.e. M is not 2-colorable), then for any function f

W k f k→+∞ -→ 1 n i∈V f i .

Spectral Resolution of the Heat Diffusion

Recall that the heat diffusion is defined as

∀ t > 0, ∂F t ∂t = -D -1 LF t = -(Id n -W )F t
Using the manifold Fourier expansion Ft def.

= U T D 1/2 F t , this differential equation can be re-written as

∂ Ft (ω) ∂t = -λ ω Ft (ω) =⇒ Ft (ω) = exp(-λ ω t) f (ω). (1.10)
This allows to study the convergence of the continuous heat equation.

Theorem 13 (Convergence of heat equation). If M is connected,

F t t→+∞ -→ 1 n i∈V f i .
Recall that the heat equation is discretized using the following explicit and implicit schemes, equations (1.6) and (1.7) Fk = (1 -δ) Fk + δ W Fk ,

((1 + δ)Id n -δ W ) Fk+1 = Fk .
These filtering iterations can be re-written over the Fourier domain as

   Fk+1 (ω) = (1 -δλ ω ) Fk (ω), Fk+1 (ω) = 1 (1+δλω) Fk (ω).
This allows to state the stability and convergence of the finite difference discretization.

Theorem 14 (Convergence of discretization). The explicit scheme is stable if δ < 1. The implicit scheme is always stable. One has

Ft/δ δ→0 -→ F t , Ft/δ δ→0 -→ F t .
with the restriction that for the explicit scheme, the mesh must not be 2-colorable.

Other Differential Equations. The manifold Fourier transform can also be used to solve the wave equation (1.8) since

∂ 2
Ft (ω)

∂t 2 = -λ ω Ft (ω) =⇒ Ft (ω) = cos( λ ω t) f (ω) + 1 √ λ ω sin( λ ω t)ĝ(ω).

Quadratic Regularization

Instead of using a PDE for regularization, one can try to find a new function that is both close to the original one f and that is smooth in a certain sense. This leads to the notion of quadratic regularization, where one uses a Laplacian as a smoothness prior on the recovered function.

Definition 20 (Quadratic regularizer). For t > 0, one defines

F q t = argmin g∈R n ||f -g|| 2 + t|| Gg|| 2 where G = GD -1/2 .
This optimization replaces f ∈ ℓ 2 (V ) by F q t ∈ ℓ 2 (V ) with small gradients. This optimization can be found in closed form by inverting a sparse linear system.

Theorem 15 (Solution of quadratic regularization). F q t is unique and

F q t = (Id n + t L) -1 f.
Over the Fourier domain, this inversion reads

F q t (ω) = 1 1 + tλ ω f (ω).
This corresponds to an attenuation of the high frequency content of f , in a way very similar to equation (1.10).

Once again, similarly to the heat equation, the spectral expression of the quadratic regularizer allows to study its convergence for large t.

Theorem 16 (Convergence of quadratic regularization). If M is connected,

F q t t→+∞ -→ 1 n i∈V f i .

Application to Mesh Compression

We have shown how the Fourier basis on meshes can be used to compute in a diagonal fashion filtering, heat diffusion and quadratic regularization. This Fourier transform is however of little interest in practice, since the original filterings (or finite difference approximation of the heat equation) are usually faster to compute directly than over the Fourier domain. The Fourier transform is thus mainly of theoretical interest in these cases since it allows to prove convergence results.

Another class of applications makes use of an orthogonal expansion such as the Fourier one to perform mesh compression. This section shows how to compute a linear M -term approximation in this Fourier basis and to do mesh compression. We refer to the survey [START_REF] Alliez | Recent advances in compression of 3d meshes[END_REF] for more advanced non-linear mesh compression methods.

The orthogonal basis U = (u ω ) ω of ℓ 2 (V ) ≃ R n , where L = U ΛU T allows to define a linear approximation as followed.

Definition 21 (Linear M -term approximation). For any M > 0, the linear M -term approximation of f is

f = n ω=1 f, u ω u ω M -term approx. =⇒ f M def. = M ω=1 f, u ω u ω .
The quality of the approximation is measured using the error decay, which can in turn be estimated using the removed coefficients

E(M ) def. = ||f -f M || 2 = ω>M | f, u ω | 2 .
A good orthogonal basis U is a basis for which E(M ) decays fast on the signals of interest. Equivalently, a fast decay of E with M corresponds to a fast decay of | f, u ω | for large ω. We recall that the Fourier atoms

∀ ω ∈ Z, u ω (x) = 1 √ 2π e iωx
are the eigenvectors of the compact, symmetric, semi-definite negative operator f → f ′′ (that should be defined on the Hilbert space of twice Sobolev derivable functions). This set of function is also an Hilbert basis of the space L 2 (R/(2πZ)) of 2π-periodic square integrable functions and a Fourier coefficient is f (ω)

def.

= f, u ω .

Approximation theory studies this linear error decay for classical functional spaces. One can for instance study the Fourier expansion over euclidean spaces.

Theorem 17 (Fourier in 1D). If f is C α regular on R/(2πZ), | f (ω)| ||f (α) || ∞ |ω| -α .
This result can be proven with a simple integration by parts. A slightly more difficult result shows that the linear approximation error decays like M -α .

Theorem 18 (Fourier approximation). If f is C α on R/(2πZ), then it exist C > 0 such that ω |ω| 2α | f, u ω | 2 < +∞ =⇒ E(M ) CM -α .
This kind of results can be extended to continuous surfaces thanks to the continuous Laplacian. We suppose that M is a surface parameterized by ϕ, and a function f = ϕ • f is defined on it. By definition, this function f is C α if f is C α in euclidean space. For a compact surface M, the Laplace-Beltrami operator ∆ M is symmetric (for the inner product on the surface), is negative semi-definite and has a discrete spectrum ∆ M u ω = -λ ω u ω for ω ∈ N. The functions {u ω } ω are an orthogonal basis for function of finite energy on the surface L 2 (M). The inner product of an arbitrary smooth function f ∈ C α (M) can be bounded using integration by parts

f, u ω = 1 λ k ω ∆ k M f, u ω =⇒ | f, u ω | ||f || C α λ α/2 ω .
1 Linear Mesh Processing

This proves the efficiency of the Fourier basis on surfaces to approximate smooth functions.

When computing the M -term approximation f M of f one removes the small amplitude Fourier coefficients of the orthogonal expansion of f . Figure 1.9 shows some examples of mesh approximation where one retains an increasing number of Fourier coefficients. Mesh compression is only a step further, since one also need to code the remaining coefficients. This requires first quantifying the coefficients up to some finite precision and then binary code these coefficients into a file. 1% 2% 3% 4% 5% 1% 2% 3% 4% Figure 1.9: Examples of spectral mesh compression.

Application to Mesh Parameterization

This section is restricted to the study of meshes that can be globally parameterized on a plane. It means that they are topologically equivalent to a 2D disk. More complex meshes should be first segmented in cells that are equivalent to a disk.

A parameterization of a continuous surface M is a bijection

ψ : M -→ D ⊂ R 2 .
A similar definition applies to a discrete mesh where one computes a 2D position ψ(i) for all the vertices i ∈ V and then interpolates linearly the mapping to the whole piecewise linear geometric mesh. This section explains the basics of linear methods for mesh parameterization. We refer to various surveys [START_REF] Floater | Surface parameterization: a tutorial and survey[END_REF][START_REF] Sheffer | Mesh parameterization methods and their applications[END_REF] for more details on mesh parameterization.

Usually, a 2D mesh is computed from range scanning or artistic modeling, so it does not come with such a parameterization. In order to perform texture mapping or more general mesh deformations, it is however important to use such a parameterization. Since many bijections are possible to layout the mesh in 2D, the mapping ψ has to satisfy additional smoothness assumptions. Classically, one requires that each coordinate of ψ has a vanishing Laplacian (it is thus harmonic) outside a set of constrained vertices that enforce boundary conditions.

Diffusion and Regularization on Surfaces

More precisely, ψ = (ψ 1 , ψ 2 ) is the solution of

∀ i / ∈ ∂M, (Lψ 1 )(i) = (Lψ 2 )(i) = 0 ∀ i ∈ ∂M, ψ(i) = ψ 0 (i) ∈ ∂D,
where ∂M is the boundary of the mesh, which consists in vertices whose face ring is not homeomorphic to a disk but rather to a half disk. This formulation requires the solution of two sparse linear systems (one for each coordinate of ψ).

The boundary condition ψ 0 (i) for i ∈ ∂M describes a 1D piecewise linear curve in the plane, that is fixed by the user. In the following, we will see that this curve should be convex for the parameterization to be bijective.

Remark 5. For such an harmonic parameterization, each point is the average of its neighbors since

∀ i, ψ(i) = 1 j w ij (i,j)∈E w i,j ψ(j).
The powerful feature of this linear parameterization method is that it can be proven to produce a valid (bijective) parameterization as long as the constrained position (boundary values of ψ) are along a convex curve.

Theorem 19 (Tutte theorem). If ∀ (i, j) ∈ E, w ij > 0, and if ∂D is a convex curve, then ψ is a bijection. One is free to use any laplacian (combinatorial, distance or conformal) as long as it produces positive weights. There is a issue with the conformal weights, which can be negative if the mesh contains obtuse triangles. In practice however it leads to the best results. The efficiency of a parameterization can be measured by some amount of distortion induced by the planar mapping. Linear methods cannot hope to cope with large isoperimetric distortions (for instance large extrusions in the mesh) since harmonicity leads to clustering of vertices.

Application to Mesh Flattening

One of the difficulty with linear parameterization methods is that they require to set up the positions of the vertices along the boundary of the mesh. In order to let the boundary free to evolve and find some optimal shape, one can replace the fixed point constraint by a global constraints of unit variance as follow

min ψ1,ψ2∈R n || Gψ 1 || 2 + || Gψ 2 || 2 with    ||ψ i || = 1, ψ 1 , ψ 2 = 0, ψ i , 1 = 0.
This optimization problem also has a simple global solution using eigenvectors of the Laplacian.

Theorem 20 (Mesh flattening solution). The mesh flattening solution is given by

Span(ψ 1 , ψ2) = Span(u 1 , u 2 ) where L = U ΛU T .
In order to compute this flattening, one thus needs to extract 2 eigenvectors from a sparse matrix. Note however that, in contrast to linear parameterization schemes, this flattening is not ensured to be bijective.

Geodesic Mesh Processing

This chapter studies the numerical computation of geodesic distance on Riemannian manifolds. It shows both the algorithms and applications to surface processing, in particular sampling, meshing and comparison of surfaces.

The original book on Fast Marching methods is [START_REF] Sethian | Level Sets Methods and Fast Marching Methods[END_REF]. For other applications to computer graphics and image processing one can see [START_REF] Osher | Level Set Methods and Dynamic Implicit Surfaces[END_REF] and [START_REF] Kimmel | Numerical Geometry of Images: Theory, Algorithms, and Applications[END_REF]. The recent book [3] treats all the details of the geometry of non-rigid surfaces, including geodesic distance computation and shape comparison. One can also see the two books [START_REF] Osher | Geometric Level Set Methods in Imaging, Vision, and Graphics[END_REF][START_REF] Paragios | Handbook of Mathematical Models in Computer Vision[END_REF] that contain review articles with some applications of Fast Marching and geodesic methods.

Manifold Geometry of Surfaces

Riemannian Manifold

We have seen in definition 1 that a parameterized surface is embedded into some Euclidean domain R k , which allows to define a local metric thanks to the first fundamental form I ϕ . It is however possible to consider directly a field of positive definite tensors on a parametric domain D = R s (in practice here s = 2 for surfaces or s = 3 for volumes). With a slight abuse in notations, we will even assimilate the resulting abstract surface M with D. Once again, we consider only surfaces globally mapped some Euclidean domain, but handling generic surfaces requires to split the manifold into overlapping charts.

Definition 22 (Riemannian manifold). A Riemannian manifold is an abstract parametric space M ⊂ R s equipped with a metric x ∈ M → H(x) ∈ R s×s positive definite.

Using the Riemannian metric, one can compute the length of a curve (γ(t)

) T t=0 ∈ M L(γ) def. = T 0 γ ′ (t) T H(γ(t))γ ′ (t)dt.
This length is invariant under a reparameterization of the curve, which means that for any strictly increasing ψ : [0, T ′ ] → [0, T ], one has

L(γ) = L(γ • ψ).
This notion of curve length is well defined for smooth curves γ, but it can be extended to piecewise smooth curves by splitting the integral into each part where the curve is smooth.

At each location x, the Riemannian tensor can be diagonalized as follow

H(x) = λ 1 (x)e 1 (x)e 2 (x) T + λ 2 (x)e 2 (x)e 2 (x) T with λ 1 λ 2 ,
and e 1 , e 2 are two orthogonal eigenvector fields. In fact, e i should be understood as direction (unoriented) field since both e i and -e i are eigenvectors of the tensor. A curve γ passing at location γ(t) = x with speed γ ′ (t) has a shorter local length if γ ′ (t) is collinear to e 1 (x) rather to an other direction. Hence shortest paths (to be defined in the next section) tends to be tangent to the direction field e 1 .

In practice, the Riemannian metric H is given by the problem one wishes to solve. In image processing, the manifold is the image plane M = [0, 1] 2 equipped with a metric derived from the image (for instance its gradient).

Geodesic Distances

The local Riemannian metric H(x) allows to define a global metric on the space M using shortest paths. This corresponds to the notion of geodesic curves.

Definition 23 (Geodesic distance). Given some Riemannian space (M, H) with M ⊂ R s , the geodesic metric is defined as

∀ (x, y) ∈ M 2 , d M (x, y) def.

= min

T >0,γ∈P T (x,y)

L(γ)

where P T (x, y) denotes the set of piecewise smooth curves joining x and y

P T (x, y) def.
= {γ \ γ(0) = x and γ(T ) = y} .

The shortest path between two points according to the Riemannian metric is called a geodesic. If the metric H is well chosen, then geodesic curves can be used to follow salient features on images and surfaces.

Definition 24 (Geodesic curve). A geodesic curve γ(t) ∈ P T (x, y) is such such that L(γ) = d M (x, y).

A geodesic curve between two points might not be unique, think for instance about two anti-podal points on a sphere. A complete surface is a surface for which for any pair of points (x, y) ∈ M 2 , it exists a geodesic in P T (x, y). In this course we consider only compact connected surfaces that are complete. For instance the space M = [0, 1] 2 -{(1/2, 1/2)}, which is not closed, is not a complete manifold for the euclidean metric since there is no shortest path in M between x = (0, 0) and y = (1, 1).

In order to perform the numerical computation of geodesic distances, we fix a starting point x 0 and consider only distance and geodesic curves from that point.

Definition 25 (Distance map). The distance map to a starting point x 0 ∈ M is defined as

U x0 (x) def. = d M (x 0 , x).
The main theorem that characterizes the geodesic distance is the following, that replaces the optimization problem of finding the minimum distance by a non-linear partial differential equation.

Theorem 21 (Eikonal equation).

If the metric H is continuous, then for any x 0 ∈ M, the map U x0 is the unique viscosity solution of the Hamilton-Jacobi equation

||∇ x U x0 || H(x) -1 = 1 with U x0 (x 0 ) = 0, (2.1)
where

||v|| A = √ v T Av.
It is important to notice that, even if the metric x → H(x) is a smooth function, the distance function U x0 might not be continuous. This is why the machinery of viscosity solution is needed to give a sense to the solution of the Hamilton-Jacobi equation. See for instance [START_REF] Crandall | User's guide to viscosity solutions of second order partial differential equations[END_REF] for an introduction to viscosity solutions and a proof of theorem 21.

It is possible to define the distance map U S (x) to a compact set of points S ⊂ M as the shortest path between x and any point in

S d M (S, x) def. = U S (x) def. = min y∈S d M (x, y).
The function U S also satisfy the Eikonal equation (2.1), but with the boundary condition U S (x) = 0 for x ∈ S. In this course, we restrict ourself to initial conditions that are either a single point S = {x 0 } or a finite number of such points S = {x i } m i=1 . Once the distance map U x0 has been computed by solving the Eikonal equation (2.1), one can extract any geodesic joining x 0 using a gradient descent on the function U x0 .

Theorem 22 (Gradient descent). The geodesic curve γ between x 1 and x 0 solves

γ ′ (t) = - H(γ(t)) -1 ∇ γ(t) U x0 ||H(γ(t)) -1 ∇ γ(t) U x0 || with γ(0) = x 1 .
The geodesic curve γ extracted using this gradient descent is parameterized with unit speed since ||γ ′ || = 1. Example 8. For an isotropic metric H(x) = W (x)Id x , the Eikonal equation and the geodesic extraction are

||∇ x U x0 || = W (x) and γ ′ (t) = - ∇ x U x0 ||∇ x U x0 || .

Mapping Between Surfaces

The notion of a parametric surface (M, ϕ) allows to compare the surface M with its parametric domain thanks to the first fundamental form. For instance, a surface is isometric to euclidean space if I ϕ = Id s for some parameterization ϕ.

However, in most situations, one wants to compare two different surfaces, maybe embedded in different euclidean spaces and of course having different parameterizations with the same parametric domain. Given two parameterized surfaces (M 1 , ϕ 1 ) and (M 2 , ϕ 2 ), a matching is a bijective function f : M 1 → M 2 . Such a matching allows to define an induced parameterization of M 2 : φ2 = f • ϕ 1 . This leads to the notion of isometry between surfaces.

Definition 26 (Isometry

). M 1 is isometric to M 2 if it exists a matching f such that ∀ u, I ϕ2 (u) = I φ2 (u).
This definition is local and requires only the comparison of the first fundamental forms. The difficulty is that it requires the use of arbitrary parameterizations ϕ 1 and ϕ 2 , although the notion of isometry is intrinsic to the surface. The following theorem states that isometry can in fact be checked in a global and intrinsic fashion by looking at geodesic lengths.

Theorem 23. f is isometric if and only if ∀ (x, y) ∈ M 2 , d M1 (x, y) = d M2 (f (x), f (y)).
This theorem thus gives a new definition of isometry that can be used in numerical applications to check if a mapping f is isometric, or even to optimize an unknown isometric mapping f between two surfaces.

One can also relax the notion of isometry and ask only for a conservation of angles. This leads to the notion of conformality between surfaces, which generalizes conformal parameterizations.

Definition 27 (Conformality). M 1 conformally equivalent to M 2 if it exists f such that ∀ u, I ϕ2 (u) = λ(u)I φ2 (u) with λ(u) ∈ R.

Numerical Computations of Geodesic Distances

In order to make all the previous definitions effective in practical situation, one needs a fast algorithm to compute the geodesic distance map U x0 . This section explains algorithms based on front propagation that enables to compute the distance map by propagating the distance information from the starting point x 0 .

Front Propagation Algorithms

Depending on the properties of the metric, one needs to consider several algorithms, that all rely on the idea of front propagation. This family of algorithms allows to order to computations in such a way that each point of the discretization grid is visited only once. This ordering is feasible for distance computation because the distance value of a grid point only depends (and can be computed) from a small number of points having only smaller distances. If one can order the grid points with increasing distance, then one gets a coherent ordering of the computations. Of course, this is not that easy since this distance ordering would require the knowledge of the solution of the problem (the distance itself). But depending on the application, it is possible to devise a selection rule that actually select at each step the correct grid point.

A front propagation labels the points of the grid according to a state S(x) ∈ {Computed, F ront, F ar}.

During the iterations of the algorithm, a point can change of label according to

F ar → F ront → Computed.
Computed points S(x) = Computed are those that the algorithm will not consider any more (the computation of U x0 (x) is done for these points). Front points S(x) = F ront are the points that are being processed (the value of U (x) ≈ U x0 (x) is well defined but might change in future iterations). Far points S(x) = F ar are points that have not been yet processed.

In practice, a front propagation algorithm requires three key ingredients: Given a point x in the grid, a local set of neighbors N eigh(x) connected to x. A priority P(x) among points x in the front, that allows to select the point to process at a given iteration. In most application, this priority is computed as the current value of the distance P(x) def.

= U (x). Section 2.3.2 shows how to change this priority in order to speed up computations. A procedure x → U pdate(x) ∈ R that computes the distance value approximating U x0 (x) knowing the value U x0 (x) for computed point and an approximate value for points in the front. This procedure usually solves some kind of equation that discretizes the Eikonal equation (2.1) one wishes to solve.

Listing 1 gives the details of the front propagation algorithm that computes a distance map U approximating U x0 (x) on a discrete grid. The following section details for actual implementations of the U pdate procedure for different metrics.

The numerical complexity of this scheme is O(n log(n)) for a discrete set of n points. This is because all the points are visited (tagged Computed ) once, and the selection of min P from the front points takes at most log(n) operations with a special heap data structure (although in practice it takes much less and the algorithm is nearly linear in time).

1. Initialization: U (x 0 ) ← 0, S(x 0 ) ← F ront, ∀ y = x 0 , S(y) ← F ar. -If S(y) = F ront, then U (y) ← min(U (y), U pdate(y)).

-Recompute the priority P(y).

Stop:

If x = x 1 , go back to 2.
Table 1: Front propagation algorithm.

Geodesic Computation on a Graph

A graph is a discrete data structure that links together points with edges. This is a crude discretization of many problems (such as those involving surfaces). In particular it is poorer than the mesh data structure introduced in section 1.1.2 since it lacks two dimensional connectivity information given by the faces of the mesh. However, in many applications including high dimensional data processing, it is the only way to capture the local connectivity information.

Definition 28 (Graph). A graph is a data structure G = (V, E) where V = {1, . . . , n} and E ⊂ V × V is symmetric.
Similarly to the notion of local averaging weights introduced in section 1.2.2, a metric on a graph is a weight defined on each edge.

Definition 29 (Graph metric).

A metric H on a graph is given by ∀ (i, j) ∈ E, H(i, j) > 0.

It is important to note that the metric H(i, j) is usually very different from the local averaging weights W (i, j) introduced in section 1.2.2. In general, H(i, j) is small to represent the fact that i and j share similarities, and in this case W (i, j) should be large (to produce strong averaging between values at i and j).

Example 9 (Geometric realization). In practice, a graph often comes with some geometric realization, which is a set of positions V = {x 1 , . . . , x n } ⊂ R k . In this case, the metric can be defined as H(i, j)

def.

= ||x i -x j ||.
If the set of points is sampled from some continuous surface, this edge metric is a crude discretization of the metric defined by the embedding of the surface in euclidean space. One can note that in this case, if one uses the local averaging weights W based on distance, one has W (i, j) = H(i, j) -2 .

Definition 30

(Path). A path π is π = (i 1 , . . . , i m ) ⊂ V m with (i k , i k+1 ) ∈ E. The starting point is S(π) = i 1 and the ending point is E(π) = i m .
It is important to note that a path π is a topological structure. In order to compute its length, one just needs to use the metric.

Definition 31 (Length). The length of a path π ⊂ V m is L(π) def. = n i=1 W (v i , v i+1 ).
This path length allows to define a geodesic metric on vertices of the graph.

Definition 32 (Geodesic distance on a graph). The geodesic distance between (x, y) ∈ V 2 is defined as

d G (x, y) = min π {L(π) \ S(π) = x, E(π) = y}
In order to compute the geodesic map U x0 (x) = d G (x 0 , x), one can use the front propagation algorithm, listing 1 with the following update procedure

U pdate(x) = min (x,y)∈E U (y) + H(x, y) .
This corresponds to the Dijkstra propagation that is classical in graph theory, see for instance [START_REF] Cormen | Introduction to Algorithms[END_REF].

Geodesic Computation on a Square Grid

The classical Fast Marching algorithm, introduced by Sethian [START_REF] Sethian | Level Sets Methods and Fast Marching Methods[END_REF], is a fast procedure to solve the Eikonal equation (2.1) for an isotropic metric H(x) = W (x)Id k for a uniform regular grid that discretizes [0, 1] k . We recall this procedure for a planar domain k = 2 although it can be extended to any dimension.

In order to capture the viscosity solution of an Hamilton Jabobi equation, one cannot use standard finite differences because of the apparition of chocks and singularities in the solution of the equation. One needs to choose, at each grid point, the optimal finite difference scheme (differentiation on the left or on the right to approximate d/dx for instance). This optimal differentiation should be chosen in the direction where the solution of the equation decreases. This is called an upwind finite difference scheme, and on a 2D grid it leads to find u = U pdate(x) at a grid point x = x i,j that is the smallest solution of

max(u -U (x i-1,j ), u -U (x i+1,j ), 0) 2 + max(u -U (x i,j-1 ), u -U (x i,j+1 ), 0) 2 = h 2 W (x i,j ) 2 .
Sethian shown that the smallest solution of this equation leads to a stable and convergent scheme that can be safely used in the front propagation algorithm listing 1.

Figure 2.2 shows some examples of front propagation with the Fast Marching. The colored area shows, at some given step of the algorithm, the set of computed points (its boundary being the set of front points). During the iterations, the front propagates outwards until all the grid points are visited. The numerical complexity of this scheme is O(n log(n)) for a grid of n points. S(x 1 ) = Computed to avoid performing useless computations. The idea of using geodesics in order to extract salient curves in images as been introduced in [START_REF] Kimmel | Finding shortest paths on surfaces using level sets propagation[END_REF].

In practice, the difficult task is to design a metric W in order have meaningful geodesics. Here are some examples of possible choices, for image processing with an input image f :

Pixel value based potential: in many application, one simply wishes to extract curves with a constant value c. In this case, one can use a potential like Gradient avoiding potential: for edge detection, one wants geodesics that avoid traversing image edges. In this case, the potential should be small in area of large gradient in the image

W (x) = 1 ε + |f (x) -c| .
W (x) = G σ2 * 1 ε + ||∇ x f * G σ1 || ,
where G σ are gaussian filter. One usually uses a very small σ 1 but a larger σ 2 to cope with noise.

Gradient attracting potential: for road tracking in satellite images, the structures to extract are usually very thin, and one thus would only like to follow regions with elongated gradient. To avoid cancellation in the gradient, only a smoothing on the magnitude of the gradient should be performed The Fast Marching works the same way in any spacial dimension k and in particular can be used to extract shortest paths in 3D volumetric medical data. Such a volume is a discretization of a mapping f : [0, 1] 3 → R. Figure 2.4 shows two alternative ways to explore a 3D cube of a scanned medical data. The first way is simply to visualize 2D slices that are traditional images, for instance (f (x, y, z 0 )) x,y for some fixed z 0 ∈ [0, 1]. An alternative way is to perform a 3D display with a semi-transparent mapping that removes more or less parts of the data. The transparency at point (x, y, z) is defined as ρ(f (x, y, z)) where ρ : [f min , f max ] → [0, 1] is the α-mapping. A third

W (x) = ε + G σ2 * ||∇ x f ||.
(x) = (|f (x) -f (x 0 )| + ε) -1
where x 0 is a point given by the user and supposed to lie inside some vessel. A geodesic follows nicely a vessel since its density is constant and thus the value of f is approximately equal to f (x 0 ) inside the vessel.

Geodesic computation on a triangulation

The classical Fast Marching algorithm is restricted to isotropic metrics on a regular grid. It can thus only be used to compute geodesics on surface parameterized conformaly on a 2D square. This setting is useful for image and volumetric data processing, but in order to deal with arbitrary surfaces embedded in R k , one needs to use algorithms designed for triangulated meshes. Kimmel and Sethian [START_REF] Kimmel | Computing Geodesic Paths on Manifolds[END_REF]have developed a version of the Fast Marching algorithm for surface M ⊂ R n with metric W (x) for x in embedded space R n . In the continuous setting, a parameteric surface (M, ϕ) embedded with a metric W (x) in ambient space correspond to a Riemannian manifold with a metric I ϕ (x)W (ϕ(x)) in parameter space x ∈ R 2 .

The algorithm of Kimmel and Sethian works on a triangulated mesh and treats the triangles of this mesh as locally flat and equipped with an isotropic metric W . In order to compute the update value at a given vertex x 0 , it computes an update value U pdate f (x 0 ) for each triangle f ∈ F x0 in the face 1-ring around x 0 , F x0 = {f 1 , . . . , f k }. The resulting Fast Marching update step is defined as U pdate(x 0 ) = min

f ∈Fx 0 U pdate f (x 0 ).
In order to derive the expression for U pdate f (x 0 ), one considers a planar triangle f = (x 0 = 0, x 1 , x 2 ) and denote as X = (x 1 , x 2 ) ∈ R 2×2 . In the following, we assume, without loss of generality that x = 0.

If only U (x 1 ) or U (x 2 ) is available for computation (meaning that one of the points has not been reached by the front), then the update procedure is simply

U pdate f (x 0 ) = min(U (x 0 ), U (x 1 ) + W (x 0 )||x 1 ||, U (x 2 ) + W (x 0 )||x 2 ||),
which corresponds to a simple Dijkstra-like update (propagation along the 1D edges of the triangle). Otherwise, the known distances are u = (U (x 1 ), U (x 2 ))

T and one wishes to solve for U pdate f (x) = p = U (x).

The discretization of the Eikonal equation leads to consider a linear approximation of the geodesic distance map U inside the current triangle f .

U (x) ≈ g, x + p where ∇ x U ≈ g.
One can evaluate this linear equation at x 1 and x 2 , which leads to the following quadratic equation in the unknown p u = X T g + p1

||g|| 2 = W (x) 2 =⇒ 1 T Q1p 2 + (21 T Qu)p + (u T Qu -W (x) 2 ) = 0,
where we have denoted Q = (X T X) -1 . This matrix Q account for the local deformation of the surface (which corresponds to the anisotropy of the metric). The only admissible solution to this problem is

p = 1Qu + (1 T Qd) 2 + 1 T Q1(u T Qu -W (x) 2 ) 1 T Q1
Some care should be taken to compute the value of the update on triangle f so that the computed value is larger than the available distances at x 1 and x 2

U pdate f (x) = p if p max(U (x 1 ), U (x 2 )), otherwise: min(U (x 0 ), U (x 1 ) + W (x 0 )||x 1 ||, U (x 2 ) + W (x 0 )||x 2 ||).
There is some technical difficulties with this scheme on triangulation that contains obtuse angles that we shall ignore here.

Applications and Extensions of Fast Marching

Applications and Extensions of Fast Marching 2.3.1 Shape Analysis

In order to analyze the shape of planar objects, one can consider the metric space obtained by restricting the plane to the inside of a planar domain.

Definition 33 (2D shape). A 2D shape S is a connected, closed compact set S ⊂ R 2 , with a piecewise-smooth boundary ∂S.

The geodesic distance inside such a shape is obtained by constraining the curve to lie inside S. where P(x, y) ⊂ S are the paths with starting point x and ending point y. The geodesic distance can be used to define several functions on the 2D shape. This section studies the eccentricity of a shape, as introduced by [START_REF] Ion | Shape matching using the geodesic eccentricity transform[END_REF] to perform shape recognition. . These eccentric points are in fact located along the boundary.

Definition 35 (Eccentricity). The eccentricity E

S : M → R is E S (x) def. = max y∈S d S (x, y) = max y∈∂S d S (x, y).

Theorem 24 (Location of eccentric points). One has E(S) ⊂ ∂S.

A more general definition of eccentricity allows to replace the maximum by a weighted average of geodesic distances. Definition 37 (α-eccentricity). The α eccentricity of some shape S is defined as

E α S (x) def. = S d S (x, y) α dy 1/α .
This eccentricity allows to generalize the notion of gravity center to the geodesic setting.

Definition 38 (Euclidean gravity center). The Euclidean gravity center is

argmin x S ||x -y|| 2 dy.
The α-eccentric center is argmin

x E α S (x).

Remark 6. For α = 2, the eccentric center is called geodesic gravity center (and equivalent to the Euclidean center for an uniform metric).

Having defined a function such as E S inside a shape S, one can collect information about the shape using the histogram of that function.

Definition 39 (Descriptors). The eccentricity histogram descriptor h S ∈ R m of a shape is

∀ i = 1, . . . , m, h S (i) = 1 |S| # x ∈ S \ i -1 m E S (x) -min(E S ) max(E S ) -min(E S ) < i m .
In particular, one can compare shapes by measuring the distance between the histograms

δ(h, h) 2 def. = m i=1 (h(i) -h(i)) 2 .
These histograms are invariant if one modifies a shape isometrically. In the plane, geodesic isometry of shapes are not interesting since they are rotations and translations. One can however consider approximate isometries such as articulations, that are useful to model deformations of planar shapes, as defined in [START_REF] Ling | Using the inner-distance for classification of articulated shapes[END_REF].

Definition 40 (ε-articulated object). An articulated object S can be split as

S = m i=1 S i i =j J ij , with diam(J ij ) ε.
Definition 41 (Articulation). An articulation is a mapping between two articulated shapes S, S ′ such that

f : S → S ′ = m i=1 S ′ i i =j J ′ ij is rigid on S i → S ′ i .
The eccentricity is approximately invariant for shapes that are modified by articulation.

Theorem 25 (Articulation and isometry). If f is an articulation, then

|d S (x, y) -d S ′ (x, y)| mε | E S (x) -E S ′ (x)| mε .
Starting from a shape library {S 1 , . . . , S p }, one can use the shape signature h S to do shape retrieval using for instance a nearest neighbor classifier, as shown in pseudo-code 2. Figure 2.10 shows examples of typical shape retrievals.More complex signatures can be constructed out of geodesic distances and un-supervised recognition can also be considered.

Heuristically Driven Propagation

The various implementations of the front propagation algorithm, pseudo-code 1, use a simple priority P(x) = U x0 (x), where U (x) ≈ d M (x 0 , x) is the current value of the distance to the starting point. This strategy leads to an isotropic grow of the front which enforces the exploration of a large area of the computational grid. The advantage of using this priority is that it does not favor any points and thus produces provable valid approximations of geodesic distance (both on a graph with Dijkstra and on a square/triangular grid with Fast Marching).

In order to reduce the computational burden, one could think about using more aggressive ordering of the front that favors some specific direction in the front. The hope is that the front would advance faster in the direction of the goal x 1 one wishes to reach. Ultimately, one would like the front to explore only points along the geodesic γ ∈ P T (x 0 , x 1 ) joining the starting point to the ending point.

If one has an oracle: V (x) ≈ d(x 1 , x) that estimates the remaining geodesic distance from the current point x to the end x 1 , one can use as priority map

P(x) = U (x) + V (x).
The map V is called an heuristic since the exact distance d(x 1 , x) is not available in practice. The value of a good heuristic close to the real distance is revealed by the following theorem.

Theorem 26 (Geodesic segment). The function ψ(x) = d(x 0 , x)+d(x 1 , x) is minimal and constant ψ(x) = d(x 0 , x 1 ) along the geodesic path joining x 0 and x 1 .

In the setting of graph theory, the Dijkstra algorithm can be replaced by the A * (A-star), [START_REF] Nilsson | Problem-solving Methods in Artificial Intelligence[END_REF], which uses an heuristic to speed up computations. The following theorem proves the validity of this approach.

Theorem 27 (A * validity). If the heuristic satisfies V (x) d(x 1 , x), then the curve γ ∈ P(x 0 , x 1 ) extracted from the front propagation, algorithm 1, is a geodesic between x 0 and x 1 .

Over a continuous domain, one can invoke a similar (but weaker) theorem.

Theorem 28 (Explored area). If the heuristic satisfies V (x) d(x 1 , x), then the geodesic γ ∈ P 1 (x 0 , x 1 ) between x 0 and x 1 satisfies

{γ(t) \ t ∈ [0, 1]} ⊂ {x \ P(x) = U (x) + V (x) P(x 1 )} .
This theorem shows why it is important to approximated the geodesic distance by below, since otherwise the region explored by the algorithm might not contain the true geodesic. Many strategies can be used to estimate an heuristic. For instance, on a Riemannian metric (M, H(x)), one could use

V (x) = ρ||x -x 1 || where ρ = min x =0,v =0 ||v|| H(x) .
In this case, ρ is the minimum eigenvalue of all the tensors H(x). This heuristic estimates the geodesic distance with an Euclidean distance and satisfies V (x) d(x 1 , x).

For a propagation on a graph (A * algorithm) that is embedded in Euclidean space according to i ∈ V → x i ∈ R k , one could also define

∀ i ∈ V, V (i) = ||x i1 -x i ||,
where i 1 is the index of the ending point. This heuristic also satisfy V (i) d(i 1 , i).

These Euclidean heuristics performs poorly on spaces that are not relatively flat. In order to compute more accurate heuristic, we use an expression of the geodesic distance as a minimization. This kind of approximation has been used first in graph theory [START_REF] Goldberg | Computing the shortest path: A* search meets graph theory[END_REF] and it is defined in a continuous setting in [START_REF] Peyré | Heuristically driven front propagation for fast geodesic extraction[END_REF]. This leads to an heuristic V (x) = d(x, x 1 ) that has the following properties. In a numerical application that requires the extraction of many geodesics in real time over a large domain, one can pre-compute (off-line) the set of distance maps to the landmarks {d(x, z i )} m i=1 . At run time, this set of distances is used to compute the heuristic and speed up the propagation. Figure 2.12 shows how the quality of the heuristic increases with the number of landmarks. Figure 2.13 shows an application to geodesic extraction on 3D meshes.

Surface Sampling

In order to acquire discrete samples from a continuous surface, or to reduce the number of samples of an already acquired mesh, it is important to be able to seed evenly a set of points on a surface. This is relevant in numerical analysis to have a good accuracy in computational simulations, or in computer graphics to display 3D models with a low number of polygons. In practice, one typically wants to enforce that the samples are approximately at the same distance ones from each others. The numerical computation of geodesic distances is thus a central tool, that we are going to use both to produce the sampling and to estimate the connectivity of a triangular mesh.

Farthest Point Sampling

A sampling of a Riemannian surface M is a set of points {x 1 , . . . , x n } ⊂ M. If the surface is parameterized by ϕ : [0, 1] 2 → M, the easiest way to compute a sampling is to seed points regularly over the parametric domain

∀ (i, j) ∈ {1, . . . , √ n} 2 , x i,j = ϕ(i/ √ n, j/ √ n).
This strategy performs poorly if the mapping ϕ introduces heavy geodesic distortion and the sampling might not be regular any more for the geodesic metric on the surface. In order to ensure the quality of a sampling, one can use the notion of a well separated covering.

Definition 42 (ε-covering). A sampling {x 1 , . . . , x n } ⊂ M is an ε-covering if i B ε (x i ) = M where B ε (x) def. = {y \ d M (x, y) ε} . Definition 43 (ε-separated). A sampling {x 1 , . . . , x n } ⊂ M is ε-separated if max(d M (x i , x j )) ε.
The farthest point sampling algorithm is a simple greedy strategy able to produce quickly a good sampling. This algorithm has been introduced in image processing to perform image approximation [START_REF] Eldar | The farthest point strategy for progressive image sampling[END_REF]. It is used in [START_REF] Peyré | Geodesic remeshing using front propagation[END_REF] together with geodesic Delaunay triangulation (to be defined in the next section) to do surface remeshing.

Pseudo code 3 gives the details of this iterative algorithm. In particular, note that the update of the distance d(x) to the set of already seeded points goes faster at each iteration since the domain of update is smaller when the number of points increases.

Initialization

: x 1 ←random, d(x) ← d M (x 1 , x), set i = 1. 2. Select point: x i+1 = argmax x d(x), ε = d(x i+1 ). 3. Local update of the distance: d(x) ← min(d(x), d M (x i+1 , x)).
This update is restricted to the set of points {x \ d M (x i+1 , x) < d(x)}.

4.

Stop: If i < n or ε > ε 0 , set i ← i + 1 and go back to 2.

Table 3: Farthest point sampling algorithm.

The output sampling of the algorithm enjoys the property of being a well separated covering of the manifold.

Theorem 31 (Farthest seeding properties). The farthest point sampling {x 1 , . . . , x n } is an εcovering that is ε-separated for

ε = max i=1,...,n min j=1,...,n d M (x i , x j ).
Note however that there is no simple control on the actual number of samples n required to achive a given accuracy ε. We refer to [START_REF] Clarkson | Building triangulations using epsilon-nets[END_REF] for an in-depth study of the approximation power of this greedy sampling scheme.

Triangulations

Having computed, for instance with farthest points, a sampling {x i } i∈V ⊂ M, the next step is to compute some connectivity between the samples in order to build a graph, or even better, a triangulation. The problem of surface remeshing has been studied extensively in computer graphics, see the survey [START_REF] Alliez | Recent advances in remeshing of surfaces[END_REF]. This section explains a solution based on the geodesic Delaunay triangulation.

The following definition generalizes the notion of an Euclidean Voronoi diagram, definition 10, to an arbitrary surface.

Definition 44 (Voronoi segmentation). The Voronoi segmentation of a sampling {x

i } i∈V ⊂ M is M = i V i with (V i ) m i=1 def. = Voronoi M ({x i } i )
where

V i def. = {x \ ∀ j = i, d M (x, x i ) d M (x, x j ).}
Each Voronoi cell V i is thus composed of points that are closer to x i than to any other sampling point. The boundary between two adjacent cells V i and V j is thus a piece of curve at equal distance between x i and x j . One can then compute the graph dual to a given partition, which joins together pair of adjacent cells. This leads to the notion of Delaunay graph.

Definition 45 (Geodesic Delaunay graph). The Delaunay graph (V, E) of a samplind {x i } i∈V ⊂ M is defined for V = {1, . . . , n} as

E = {(i, j) ∈ V \ ∂V i ∩ ∂V j = 0} .
The main interest of this Delaunay graph is that, if the number of points is large enough to capture the topology of the surface (for instance at least 4 points are needed on a sphere), then one gets a valid triangulation.

Theorem 32. For a large enough number of points, the Delaunay graph is a valid triangulation.

This theorem means that one can find a set of faces F such that (V, E, F ) is a triangulated mesh. One can see [START_REF] Leibon | Delaunay triangulations and voronoi diagrams for riemannian manifolds[END_REF] for a theoretical study of geodesic Delaunay triangulations. This Delaunay triangulation can thus be used to perform a geodesic meshing or re-meshing of any Riemannian surface, as explained in [START_REF] Peyré | Geodesic remeshing using front propagation[END_REF]. Figure 2. [START_REF] Ion | Shape matching using the geodesic eccentricity transform[END_REF] shows examples of Voronoi segmentations on the plane for various isotropic Riemannian metrics W (x). The Delaunay graph allows to define a planar mesh of points evenly sampled according to the metric. In order to mesh the interior of a planar shape S ⊂ R 2 , one can use the Euclidean metric inside the shape and compute a geodesic Delaunay triangulation. Figure 2.17 shows some examples of shape meshing with this uniform metric. This triangulation is however very close to the usual definition of a planar Euclidean Delaunay triangulation. In contrast, one can use a non-uniform metric W (x) and compute a sampling inside the shape that conforms itself to this density. Figure 2.17 shows a sampling and meshing that uses a metric W (x) = (ε + d(x, ∂S)) -1 that tends to seed more points on the boundary of the shape S. An option to compute this metric is to use a texture mapped on the surface. Starting from some parametric surface:

Metric

ϕ : D ⊂ [0, 1] 2 → M, a texture T is a mapping T : [0, 1] 2 → R.
It allows to define an isotropic metric using for instance an edge adaptive function

∀ x ∈ D, H(x) = ψ T (x)Id 2 .
where the edge-based stopping function is ψ T (x) = (||∇ x T || + ε) -1 . Figure 2.20 shows examples of remeshing with a texture-adapted metric with a decreasing value of ε (increasing adaptivity).

Centroidal Tesselation

For some applications, the sampling quality provided by the greedy farthest point sampling, pseudo code 3, might not be good enough. One can thus try to enhance its quality by a relaxation scheme that updates the position of each sample in order for them to be globally more evenly distributed. To give a precise definition of this sampling quality, one can use the geodesic gravity center, already introduced in definition 38 in the special case of a planar shape.

Definition 46 (Geodesic gravity center). The geodesic center of gravity of

A ⊂ M is g M (A) = argmin x A d M (x, y) 2 dy.
For an Euclidean metric where H(x) = Id, one gets the traditional definition of gravity center

g eucl (A) = argmin x A ||x -y|| 2 dy. =⇒ g eucl (A) = 1 |A| A ydy,
where |A| is the Euclidean area of A.

A good sampling is then defined as a sampling whose Voronoi cells are compact and rounded, which is equivalent to each cell having a gravity center equal to its seed.

Definition 47 (Centroidal tesselation). A sampling {x

i } m i=1 is centroidal if its Voronoi tessellation Voronoi M ({x i } i ) = {V i } m i=1 satisfies x i = g M (V i ).
One can see [START_REF] Du | Centroidal Voronoi tessellations: Applications and algorithms[END_REF] for an overview of the theory and applications of centroidal tesselations with weighted Euclidean metrics. The geodesic centroidal tesselation has been introduced for surface remeshing and segmentation in [START_REF] Peyré | Geodesic remeshing using front propagation[END_REF]. See also [START_REF] Alliez | Recent advances in remeshing of surfaces[END_REF] for other ways to use centroidal tesselations for remeshing.

This definition can be turned into an energy minimization as follow.

Theorem 33 (Variational characterization). A centroidal tesselation is a local minimizer of

E({x i }, {V i } i ) = m i=1 Vi d M (x i , y) 2 dy.
Such a local minimizer can be computed by a gradient descent of E, using

E i (x) = Vi d M (x i , y) 2 dy =⇒ ∇ x E i = Vi d M (x i , y) ---→ n x (y)dy,
where n x (y) = γ ′ (0)/||γ ′ (0)|| is the unit tangent at x to the geodesic γ ∈ P T (x, y) joining x to y.

Semi-

{x i } i ← random. 2. Update regions: ∀ i, {V i } i ← Voronoi M ({x i } i ). 3. Update centers: ∀ i, x i ← g M (V i ).
4. Stop: while not converged, go back to 2. , where ρ : [0, 1] 2 → R + is a user-defined weight function. Not that this definition does not corresponds to a geodesic gravity center with an isotropic geodesic metric H(x) = ρ(x)Id 2 , but this usually gives similar results. In practical situation, one wants to choose a weigth function W (x) or ρ(x) in order for the size of the triangles near a position x to conform to some sizing field s(x). As explained in [START_REF] Du | Centroidal Voronoi tessellations: Applications and algorithms[END_REF], the weights should obey approximately the scaling relation d+2) where d is the dimension of the domain (d = 2 for a 2D surface mesh).

ρ(x) ≈ W (x) ∼ s(x) -(

High Dimensional Data Analysis

This chapter presents various methods to analyze point clouds in high dimension. The emphasis is put on manifold models that allows to parameterize these point clouds with a small number of parameters. Examples in image processing shows that some image ensemble indeed satisfy these manifold models.

Manifold Models of High Dimensional Data

Flattening and Dimension Reduction

In order to perform data analysis, one often has to reduce the dimension of the input set of points. This can be useful to extract salient features from a data set and to speed up processing. This section shows how several related problems all lead to dimensionality reduction using spectral methods. In particular, surface flattening can be re-casted as a matrix approximation problem with dimensionality reduction.

Point Clouds

A data set is a large point clouds X = (x j ) p j=1 ∈ R n×p with x j ∈ R n , n ≪ p. From this set of input data, one can define two matrices that represent in transposed forms the correlation present in the input data (second order statistics).

Definition 48 (Correlation matrix). C = XX T ∈ R n×n , C ij = k x k (i)x k (j). Definition 49 (Gram matrix). K = X T X ∈ R p×p , K ij = ( x i , x j ) i,j .
It is important to remember that these matrices do not carry the same amount of information and depending on the application (flattening, dimension reduction, best approximation, etc) one can use either the correlation or the gram matrix. An important tools from numerical analysis is the singular value decomposition (SVD) that generalizes the orthogonal diagonalization of symmetric matrices to non-symmetric, non-square matrices. Such an orthogonal decomposition is crucial to perform matrix approximation by deleting small spectral values.

Theorem 34 (Singular value decomposition). Any matrix X can be decomposed as

X = U ΣV T = U Σ n V n T with U ∈ O(n), V ∈ O(p),
where O(n) ⊂ R n×n is the group of orthogonal matrices and

Σ = diag(σ i ) i , Σ n = diag(σ 1 , . . . , σ n ), V n = (v 1 , . . . , v n ) ∈ R p×n .
This decomposition satisfies

XX T = U ΛU T X T X = V ΛV T .
The singular values Λ = Σ 2 are ordered λ 1 . . . λ n λ n+1 = 0.

In order to compare point clouds, one can use various matrix norms, among which the most popular are the Frobenius norm and the norm derived from the vectorial norm.

Definition 50 (Matricial norms).

||X|| 2 def. = argmin ||a||=1 ||Xa|| and ||X|| 2 F = ij X 2 ij = tr(X T X) = tr(XX T ).
The properties of these norms are recalled in this theorem.

Theorem 35. These norms are invariant under isometries, since for any orthogonal matrix U ,

||U X|| F = ||X|| F and ||U X|| 2 = ||X|| 2 .
These matricial norms satisfy

||X|| 2 = σ 1 and ||X|| 2 F = i σ 2 i .

Linear Dimension Reduction

In order to perform linear approximation of some data set, one needs to use orthogonal projectors.

Definition 51 (Orthogonal projectors). The set of orthogonal projectors of rank k is defined as

P k def. = P = AA T \ A ∈ R n×k and A T A = Id k
The optimal projection of a given data set can be computed easily using a spectral decomposition.

Theorem 36 (Dimension reduction solution). Up to a rotation, U k = (u 1 , . . . , u k ) satisfies

U k U k T = argmin P ∈P k ||X -P X|| F = argmin P ∈P k ||X -P X|| 2 .
The Principal component analysis (PCA) of a data set X corresponds to reducing the dimensionality using the optimal linear projector

X ∈ R n×p -→ X k = U k T X ∈ R k×p (dim. reduction) X k ∈ R k×p -→ X def. = U X k (reconstruction)
The error that occurs from this linear reduction / reconstruction cycle can be monitored using the singular values that have been removed during reduction Theorem 37 (PCA error). One has ||X -X|| 2 = i>k σ 2 i .

It is important to remember that the PCA approximation is optimal only in a linear setting. for input exemplar x j that corresponds to smooth function (or similarly for data set that are close to gaussian), this linear approximation is very good and is close to the low-frequency Fourier expansion. However for complex data x j such as points on non-smooth meshes or natural images, this linear setting performs poorly. Non-linear methods such as wavelets can overcome this drawback.

Distance and Gram Matrices

An input data set is treated as a points cloud X = (x j ) p j=1 ∈ R n×p with x j ∈ R n . In order to manipulate this points cloud, this section considers two kinds of euclidean invariants: inner products and distances. Inner products matrices are easy to deal with since they form a convex cone that can be mapped back to the point cloud X using the diagonalization of symmetric operators. In contrast, distance matrices form a complex set difficult to handle. To ease the computation with distance matrices, we consider a simple mapping between inner products and distances.

Definition 52 (Gram matrices). The set of Gram matrices is defined as

K n def. = X T X \ X ∈ R n×p ⊂ R p×p and K = 1 n p K n .
The nice property of inner product matrices is that they can be characterized by their spectral content.

Theorem 38 (Kernel factorization). K ∈ K n ⇐⇒ K is positive semi-definite and rank(K) = n.

Remark 7 (Description of the factorization). In particular, an inner product matrix can be mapped back to a point cloud X using diagonalization.

K = V ΛV T rank(K) = n =⇒ K = X T X with X = √ Λ n V n V n = (v 1 , . . . , v n ) .
An other important set of matrices is those that contain pairwise distances between points in Euclidean space.

Definition 53 (Distance matrices). The set of euclidean distance matrices is defined as

D n def. = D = (||x i -x j || 2 ) i,j \ X ∈ R n×p ⊂ R p×p and D = 1 n p D n .
It is not obvious to check wether a given matrix is a distance matrix or not. At least it has to satisfy the following properties of a pre-distance matrix.

Definition 54 (Pre-distance matrices). The set of pre-distance matrices is defined as

D + def. = D \ D 0, D T = D, diag(D) = 0 .
Remark 8. One has the inclusions:

D 1 ⊂ D 2 ⊂ . . . ⊂ D p ⊂ D + .
Using the expansion

||x i -x j || 2 = ||x i || 2 + ||x j || 2 -2 x i , x j = K ii + K jj -2K ij ,
one can define a mapping between kernels and distance matrices as follow

ϕ : K n -→ D n K -→ KI + IK -2K
where I = 11 T .

Remark 9. One can readily note that

ϕ((X + λ) T (X + λ)) = ϕ(X) =⇒ ϕ is not injective. ∀ D ∈ D, ∃K, D = ϕ(K) =⇒ ϕ is surjective.
In order to make this mapping bijective, one has to consider only centered data sets.

Definition 55 (Centered points). The set of centered points clouds is

Pn def. = X ∈ R n×p \ X1 = 0 .
Definition 56 (Centered inner products). The set of centered inner product matrices is

Kn def. = X T X \ X ∈ Pn .
One can center a points cloud and a gram matrix using the following projectors

R n×p -→ Pn X -→ XJ and K n -→ Kn K -→ JKJ .
where J def.

= Id p -I/p.

The following theorem shows that centered inner product and distance matrices are in bijection.

Theorem 39 (Mapping between gram and distance). One has the following bijections

ψ • ϕ = Id Kn , ϕ • ψ = Id Dn , where ψ : D n -→ Kn D -→ -1 2 JKJ.

Manifold Flattening

In order to flatten a manifold M, we forget about the location of the points on the manifold but only consider pairwise geodesic distances using the metric d M . This gives an input geodesic distance matrix D = (d M (x i , xj ) 2 ) i,j ∈ D + , where xi ∈ M.

The flattening corresponds to finding

X = (x i ) p i=1 ∈ R n×p such that ϕ(X T X) def. = D ≈ D.
The underlying idea is to find an euclidean layout of points that respects the geometry of the manifold (the pairwise distances). The input matrix D ∈ D + is not euclidean. To find the layout, one would ideally like to solve min This optimization has a simple spectral solution, that has been exploited in the Isomap algorithm, introduced in [?] (listing ??). This corresponds to the classical multidimensional scaling theory.

Theorem 40 (Isomap solution). One has

K = U Λ n U T where ψ( D) = U T ΛU where Λ = (λ i ) i , |λ 1 | . . . |λ p | and Λ n = diag(λ 1 , . . . , λ n , 0, . . . , 0). 1. Compute geodesic distances D = (d M (x i , xj ) 2 ) i,j .
2. Compute centered kernel K = -J DJ/2 with J = Id -I/p.

Compute spectral decomposition

K = U ΛU T , U = (u j ) j ∈ R p×p .
4. Embedding locations are

X = Λ n U n T = diag( λ 1 , . . . , λ n )(u 1 , . . . , u n ) T ∈ R n×p .
Table 5: Isomap algorithm.

Remark 10 (Equivalence PCA/Isomap). If D is euclidean, i.e. D = ϕ(X T X) ∈ D, then isomap is equivalent to PCA.

A flattening is a mapping f = (f 1 , f 2 ) → R 2 . Figure 1.11, right, shows an example of mesh flattening using Isomap. Such a flattening has been used to perform texture mapping [?]. Figure 1.11 also compares this flattening method with the Laplacian-based methods introduced in section 1.3.8.

Laplacian eigenmaps enforces local smoothness through

f i = argmin ||f ||=1 ||Gf ||. In this case, (f 1 , f 2 ) are eigenvectors (#2,#3) of L = G T G. Isomap enforces global constraints through ||f (x) -f (y)|| ≈ d M (x, y). In this case, (f 1 , f 2 ) are eigenvectors (#1,#2) of -J(d M(xi,xj ) ) ij J.
In laplacian methods, the matrix to process is sparse, which eases computations but may leads to poor global control over the flattening.

Bending Invariants

Instead of performing a reduction in dimension (flattening) one can use the Isomap embedding in order to obtain a representation of the surface that is invariant under geodesic isometries, as introduced by [?].

Definition 57 (Bending invariants). Given a surface M, the bending invariant I M (x) ∈ R 3 is the output of the Isomap algorithm.

A geodesic isometry is formally a mapping ψ : M → M ′ such that

∀ (x, y) ∈ M 2 , d M (x, y) = d M ′ (ψ(x), ψ(y)).
The main properties of bending invariants is the following theorem.

Theorem 41 (Bending invariance). Up to rigid motion, I M is invariant to geodesic isometries:

I M (x) = v + U I M ′ (ψ(x)) where U ∈ O(3) and v ∈ R 3 . Figure ??
shows some examples of bending invariants. These bending invariants can be used to perform shape recognition with isometric invariance, which is a useful property in face recognition for instance, see [3].

Multiresolution Mesh Processing

This chapter shows how computations on a mesh can be performed in a multiscale manner, by considering meshes of increasing resolutions. This leads to the notion of subdivision surfaces and wavelet transform, which are two different tools to interpolates and decompose functions on meshes. Both methods rely on a special kind of meshes whose triangulations can be obtained by applying a regular refinement rule. In order to perform multiscale mesh processing, one needs to pack the vertices V of a topological mesh M = (V, E, F ) in sets of increasing resolution. As explained in section 1.1.2, it is important to remember that this construction is purely combinatorial, in that no geometrical information (such as actual positions of the vertices in R 3 ) is required to build the set of multi-resolution meshes. In fact these multiscale grids can be used to actually process the geometrical realization M of the mesh M as three real valued functions (the three coordinates of the points).

Semi-regular Meshes

We thus consider a set of nested indexes

V 0 ⊂ V -1 ⊂ . . . ⊂ V L = V which are split according to V j = V j+1 ∪ H j+1 .
Next section describes how to actually compute this set of nested grids using a triangular split, but most of the mathematical tools are in fact valid for arbitrary set of indices, as long as they are embedded in one each other through scales.

For mesh processing, an index ℓ ∈ V j corresponds to a vertex x ℓ ∈ V ⊂ R 3 . The signals to be processed are vectors f ∈ R n of size n = |V L | defined on the grid V L . We sometimes write f ∈ ℓ 2 (V L ) instead of f ∈ R n to emphasis the domain on which f is indexed. This chapter describes transforms for signals f ∈ ℓ 2 (V L ) sampled on the finest grid V L . γ(e) The combinatorial structure of a triangular mesh is defined in section 1.1.2. This chapter considers only a certain class of meshes M = (V, E, F ) that can be obtained by a regular split of faces, starting from an initial coarse triangulation. This splitting leads to a set of multiresolution meshes M j = (V j , E j , F j ) for J j 0, where the full mesh is M J = M . Starting from this coarse triangulation, one defines by subdivision a multiscale triangulation (V j , E j , F j ) L j 0 where For each edge e ∈ E j , a central index γ(e) ∈ V j-1 is added to the vertices

Regular 1:4 Subdivision

µ 1 (f ) µ 2 (f ) µ 3 (f ) µ 4 (f ) f subdivision subdivision
j = 0 j = -1 j = -2 j = -3
V j-1 = V j ∪ {γ(e) \ e ∈ E j } .
Each edge is subdivided into two finer edges ∀ e = (a, b) ∈ E j , σ 1 (e) = (a, γ(e)) and σ 2 (e) = (b, γ(e)).

The subdivided set of edges is then

E j-1 = {σ i (e) \ i = 1, 2 and e ∈ E j } . Each face f = (a, b, c) ∈ F j is subdivided into four faces µ 1 (f ) = (a, γ(a, b), γ(a, c)), µ 2 (f ) = (b, γ(b, a), γ(b, c)), µ 3 (f ) = (c, γ(c, a), γ(c, b)), µ 4 (f ) = (γ(a, b), γ(b, c), γ(c, a)).
The subdivided set of faces is then The set of vertices can be classified as Regular vertices are those who belong neither to the coarse mesh V 0 nor to a boundary of a mesh M j . These vertices have always 6 neighbors. Extraordinary vertices are the initial vertices of V 0 . They exhibit arbitrary connectivity. Boundary vertices are those belonging to a mesh boundary. Boundary vertices not in V 0 always have 4 immediate neighbors. Obviously not every meshes can be obtained from such a subdivision process. In practice, an arbitrary mesh, obtained from CAD design or range scanning usually does not have any multiscale structure. It is thus necessary to remesh it in order to modify the connectivity of the mesh. During this process, the position of the vertices in R 3 is modified in order for the geometrical realization to stay close from the original piecewise linear surface. One can see [START_REF] Alliez | Recent advances in remeshing of surfaces[END_REF] fur a survey of various semi-regular remeshing methods.

F j-1 = {µ i (f ) \ i = 1, 2, 3, 4 and f ∈ F j } .

Spherical Geometry Images

Starting from some input surfaces S ⊂ R 3 , one typically wants to compute a semi-regular meshes (M j ) j L that approximate S. In most case, the surface S is actually given as an arbitrary triangulated mesh and this process corresponds to a semi-regular remeshing. Many algorithm have been devised for surface remeshing and we describe here a method [START_REF] Praun | Spherical parametrization and remeshing[END_REF] that works for surfaces that have the topology of a sphere. It means that the surface has genus 0, without boundary and without handles. This methods works by computing several intermediate surface-wise parameterization.

Spherical parameterization: each points of the original triangulation of S is mapped onto the unit sphere. This create a bijective parameterization

ϕ S : S 2 → S.
This is a non-linear process that differ from the planar parameterization introduced in section 1.3.7. We do not give the details of such a process, but it requires minimizing the smoothness of the mapping ϕ -1 S under the constraint that it maps points of S to unit length vectors (point on the sphere S 2 ). The algorithm is explained in details in [START_REF] Praun | Spherical parametrization and remeshing[END_REF]. Spherical-tetraedron flattening: one flatten each quadrant (1/8) of the sphere in order to have a mapping ϕ T : Octaedron → S 2 .

One can use for instance a mapping between spherical barycentric coordinate on each quadrant and Euclidean barycentric coordinates on each face of the octahedron.

Tetraedron unfolding: One maps each equilateral face of the octaedron on a rectangular triangle that corresponds to 1/8th of the square [0, 1] 2 ϕ U : [0, 1] 2 → Octaedron.

Regular sampling: the geometry image is obtained by regularly sampling the square on a uniform grid

x ℓ = ϕ S • ϕ T • ϕ U (ℓ/n) for ℓ i = 0, . . . , n -1.
The mapping ℓ → x ℓ ∈ R 3 is the geometry image, which can be stored as a 3-channel (color) image.

From such a geometry image x ℓ , one can easily compute a semi-regular mesh by simply performing a regular 1:4 subdivision of the octaedron. Figure 3.4 shows the steps of the construction of a geometry image, and the resulting semi-regular mesh.

Subdivision Curves

Before getting into the detail of subdivision surfaces, we describe the subdivision process in the simpler setting of 1D signals. This leads to the construction of subdivision of 1D functions and subdivision curves.

In this 1D setting, the grid point indexes are dyadic sub-grids of Z ∀ j L, V j = ℓ2 j-L \ 0 ℓ < s 0 2 -j , where s 0 = |V 0 | is the size of the initial vector f 0 to be subdivided. Each subdivision steps computes, from a set f j (ℓ) ∈ ℓ 2 (V j ) of coarse values, a refined vector f j-1 ∈ ℓ 2 (V j-1 ) defined by

h 0 h 1 h-1 h1 h0 V j H j V j-
∀ k ∈ H j , f j-1 (k) = t f j ((k -1)/2 + t)h(t), ∀ ℓ ∈ V j , f j-1 (ℓ) = t f j (ℓ + t) h(t).
where the set of weights h and h acts as local averaging operators. This averaging should be corrected at the boundary, and we use here cyclic boundary conditions which identifies 0 and s 0 2 -j in V j . Figure 3.5 shows a graphical display of these averaging operators.

One can write this subdivision steps as convolution by introducing the global set of weights g = [. . . , h(-1), h(0), h(0), h(1), h(1), . . .] since one has f j-1 = (f j ↑ 2) * g where a ↑ 2 = [. . . , 0, a(-1), 0, a(0), 0, a(1), 0, . . .].

This corresponds to the traditional description of the wavelet low-pass filtering [START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF]. One can apply this subdivision of functions to a pair of signals

(X 0 , Y 0 ) : V 0 → R 2
which is a control polygon composed of points located in the plane. The subdivision curve converges to the limiting curve (X j , Y j )

j→-∞ -→ (X(t), Y (t)) 1 t=0 ⊂ R 2 .
An interesting property is that this curve is included in the convex hull of the control polygon (X(t), Y (t)) t ⊂ Conv(X 0 , Y 0 ).

Subdivision Surfaces

Subdivision schemes allows to compute a set of progressively refined vectors on a semi-regular mesh. More precisely, from an initial vector f 0 ∈ R |V0| defined on the coarse mesh M 0 , local interpolation kernels computes iteratively vectors f j ∈ R |Vj | of finer resolution. When applied to 3 function (f i 0 ) i=1,2,3 defining the geometrical position of points in R 3 , this hierarchical construction defines a subdivision surface. These subdivisions surfaces are used extensively in computer aided geometry and computer graphics. One can see [START_REF] Schroeder | Subdivision surfaces in character animation[END_REF] for a survey of subdivision surfaces and their applications.

Interpolation Operators

In order to refine a vector f j ∈ R |Vj | defined on the vertex V j of the mesh M j , one uses two interpolators

P j : ℓ 2 (V j ) -→ ℓ 2 (H j ) and Pj : ℓ 2 (V j ) -→ ℓ 2 (V j ). ( 4.1)
A new refined function f j-1 ∈ R |Vj-1| defined on the vertices V j-1 = V j ∪ H j of M j-1 is defined by applying these two refinement operators:

∀ ℓ ∈ V j-1 , f j-1 (ℓ) = (P j f j )(ℓ) if ℓ ∈ V j , ( Pj f j )(ℓ) if ℓ ∈ H j .
Since V j ⊂ V j-1 , the operator Pj only modify slightly the value at vertex in V j . On the other hand, the operator P j creates new value at the vertices of H j that are inserted between V j and V j-1 .

In practical applications, these interpolating operators are local, meaning that the value of (P j f j )(ℓ) and ( Pj f j )(ℓ) depends only on values f j (ℓ ′ ) for ℓ ′ ∈ V j being close to ℓ ∈ V j-1 , typically in the 1-ring or 2-ring vertex neighborhood.

A particularly important setting for subdivision scheme is when one apply the subdivision steps in parallel to three vectors (X j , Y j , Z j ) starting from three initial vectors describing the position in 3D space of a coarse mesh M 0 . This allows to defines finer and finer spacial localization for the vertex of the refined meshes M j . Figure 3.9 shows an example of such a subdivision surface. In order for the resulting infinitely refined surface to have good properties such as being continuous and even smooth, one needs to design carefully the interpolation operators. Next section gives examples of such operators.

Some Classical Subdivision Stencils

In order to define the interpolation operators P j and Pj of equation (3.1), one needs to use a naming convention for the neighborhoods of vertices.

For a vertex ℓ ∈ V j , the one ring neighborhood V ℓ has already been defined in equation (1.1). It is the set of vertices adjacent to ℓ. In a regular point (that does not belongs to V 0 and not on a boundary of the mesh), its size is |V ℓ | = 6 since a point has 6 neighbors. This 1-ring is used to define Pj .

Subdivision Surfaces

For a vertex k ∈ H j ⊂ V j-1 , the butterfly neighborhood is a set of vertices in V j close to k. This neighborhood is used to define P j . The two immediate neighbors are

(v 1 k , v 2 k ) def. = {v ∈ V j \ (v, k) ∈ E j-1 } .
Two other vertices (w 1 k , w 2 k ) are defined using the two faces adjacent to edge (

v 2 k , v 2 k ) ∈ E j f 1 k = (v 1 k , v 2 k , w 1 k ) ∈ F j and f 2 k = (v 1 k , v 2 k , w 2 k ) ∈ F j .
For edges E j on the boundary of M j , one one face is available, in which case we implicitly assume that f 1 = f 2 (reflecting boundary conditions). The four last vertices are defined using faces adjacent to f 1 and f 2 :

∀ i, j = 1, 2, f i,j k def. = (z i,j k , v j k , w j k ) ∈ F j with f i,j k = f j .
Once again, reflecting boundary condition are applied for faces on the boundary of the mesh. The butterfly neighborhood is depicted on figure 3.8.

k v 1 k v 2 k f 2 k f 1 k w 1 k w 2 k z 1,1 k f 1,1 k f 1,2 k f 2,2 k f 2,1 k z 1,2 k z 2,2 k z 2,1 k Figure 4.8: The butterfly neighborhood of a vertex k ∈ H j .
Linear Interpolating Scheme The simplest subdivision rule compute values along edge mide point using a simple linear interpolation as follow

∀ k ∈ H j , (P j f j )(k) = 1 2 (f (v 1 k ) + f (v 2 k )), ∀ ℓ ∈ V j , ( Pj f j )(ℓ) = f j (ℓ). (4.2)
Since Pj is the identity operator, this scheme is called interpolating. It means that value of f 0 on points of the coarse triangulation are kept during iteration of the subdivision.

Butterfly Interpolating Scheme The linear scheme creates function that are piecewise linear on each face of the coarse triangulation F 0 . In order to create smooth surface, one needs to use more points in the butterfly neighborhood as follow Loop Approximating Scheme In order to gain flexibility in the subdivision design, one can also modify points in V j during the iterations. This means that Pj is not any more the identity, and that all the values will evolves during the iterations. The question of wether these iterated modification actually converge to a limit value is studied in the next section.

∀ k ∈ H j , (P j f j )(k) = 1 2 2 i=1 f (v i k ) + 1 8 2 i=1 f (w i k ) -1 16 2 i,j=1 f (z i,j k ), ∀ ℓ ∈ V j , ( Pj f j )(ℓ) = f j (ℓ). (4.3)
The Loop subdivision rule is defined as

∀ k ∈ H j , (P j f j )(k) = 3 8 2 i=1 f (v i k ) + 1 8 2 i=1 f (w i k ), ∀ ℓ ∈ V j , ( Pj f j )(ℓ) = (1 -|V ℓ |β |V ℓ | )f j (ℓ) + β |V ℓ | ℓ ′ ∈V ℓ f j (ℓ ′ ).
where the weights depends on the number of neighbors and are defined as

β m def. = 1 m 5 8 - 3 8 + 1 4 cos(2π/m) 2 .
Other schemes. It is possible to define subdivision schemes using rules that do not involve a regular 1:4 splitting of each coarse face. For instance, in dual schemes such as the one depicted in figure 3.11, the faces of F j are not included in F j-1 but only in F j-2 .

Invariant Neighborhoods

In order to study the convergence of subdivision schemes, one needs to consider independently each vertex x ∈ V j0(x) , where j 0 (x) is the coarser scale at which x appears

j 0 (x) = max {j \ x ∈ V j } .
Original vertices satisfy j 0 (x) = 0 and are the only one (except boundary vertices) that have a non-regular connectivity.

The vertex x belongs to the mesh M j0(x) which is going to be refined through scales j < j 0 (x). In order to analyze this refinement, one needs to define an invariant neighborhood V x j ⊂ V j of x for each scale j j 0 (x). These neighborhood are the set of points that are required to compute the operators P j and Pj . More precisely, given a vector f ∈ ℓ 2 (V j-1 ), the neighborhoods are required to satisfy

∀ ℓ ∈ V x j-1 ∩ V j , ( Pj f )(ℓ) depends only on V x j ∀ k ∈ V x j-1 ∩ H j , (P j f )(k) depends only on V x j .
We further impose that all the invariant neighborhoods have the same size

∀ j j 0 (x), #V x j = m x .
ℓ , as defined in (1.2).

Thanks to the invariance of these neighborhood systems, one can restrict the predictors around x and define

P x j : V x j -→ V x j-1 ∩ V j and P x j : V x j -→ V x j-1 ∩ H j .
The subdivision matrix S x j ∈ R mx×mx is then defined as matrix of the following mapping ( P x j , P x j ) : V j x -→ V j-1

x .

All the subdivision schemes studied in this chapter are invariant, meaning that the subdivision rule does not change through the scales j. This impose that the subdivision matrices are constant S x j = S x . In fact, in all the examples given in the previous section, they only depends on the number |V x | of neighbors in the one ring of x. .12: Invariant neighborhood V x j and V x j-1 (indexing with red circles) of the Loop subdivision scheme for a vertex of valence |V ℓ | = 0. The number in {0, . . . , 9} refers to the numbering of the vertices in V x j and V x j-1

Convergence of Subdivisions

The value at x ∈ V j0(x) of a function f j ∈ ℓ 2 (V j ) obtained by subdividing at scale j j 0 (x) an initial vector f 0 ∈ ℓ 2 (V 0 ) can be computed as

f j (x) = S x f x j+1 (x) = (S x ) j0(x)-j f x j0(x) (x),
where the vector f x i ∈ R mx is the restriction of f i to the set V x j . In order to analyze the limiting function resulting from an infinite number of subdivision, one can use the eigen vector decomposition of the matrix S x S x = ΦV ΛΦ T where Φ T = Φ-1 , Λ = diag(λ i ), λ 1 λ 2 . . . λ mx .

Since the subdivision matrix S x is not symmetric, some of the eigenvalues might be complex, and we shall ignore this difficulty here. The fact that P j and Pj are predictor implies that the subdivision matrix has to satisfy S x 1 = 1, meaning that φ1 = 1 is an eigenvector associated to the eigenvalue 1. In the following we further makes the following assumption

1 = λ 1 < λ def. = λ 2 = λ 3 < λ 4 . (4.4)
This hypothesis is satisfied by all the subdivision rules introduced in the previous section.

If one write Φ = (ϕ i ) m i=1 and Φ = (ϕ i ) m i=1 , one has the following decomposition of a vector

f ∈ R mx f = mx i=1 f, ϕ i φi and (S x ) k (x) = mx i=1 λ k i f, ϕ i φi .
One thus has the following asymptotic expansion

1 λ k (f -f, ϕ 1 1) = f, ϕ 2 φ2 + f, ϕ 3 φ3 + o(1). (4.5)
This expression describes the asymptotic behavior of the subdivision scheme at zero order (position) and first order (tangents).

Theorem 42 (Convergence of the subdivision scheme). If the subdivision matrix S x of a point x satisfies (3.4) then the subdivision process converges at x to the value

f j (x) j→-∞ -→ f x j0(x) , ϕ 1 .
The smoothness of the resulting function is more difficult to analyze. A particularly important setting is when one computes the subdivision of 3 function p 0 = (X 0 , Y 0 , Z 0 ) ∈ ℓ 2 (V 0 ) 3 corresponding to the position in R 3 (geometrical realization) of a coarse mesh M 0 . In this case, the subdivided functions p j = (X j , Y j , Z j ) gives refined 3D meshes that converge uniforlmy to a continuous surfaces

p(x) = (X(x), Y (x), Z(x)) = ( X x j0 , ϕ 1 , Y x j0 , ϕ 1 , Z x j0 , ϕ 1 ).
Condition (3.4) nearly implies that the resulting surface is smooth. Indeed, the asymptotic expansion (3.5) shows that for a point x ′ near x in the subdivision domain, the differential vector can be well approximated as a projection on a 2D plane

p(x) -p(x ′ ) + o(1) ∈ Span(τ x 2 , τ x 3 ) where τ i (x) def. = ( X x j0 , ϕ i , Y x j0 , ϕ i , Z x j0 , ϕ i ).
If the vectors τ x 2 and τ x 3 are linearly independent, they form a basis of the tangent plane at p(x).

Example of the Loop subdivision. For the Loop interpolation operators defined in equation (3.3.2), the invariant neighborhood V x j correspond to the 2-ring of x in the triangulation G j , as shown in figure 3.12. For a vertex with k neighbors, |V x | = k, the size of these invariant neighborhood is m x = 3k + 1. A particular neighboring for k = 3 is depicted in figure 3.12, together with an indexing in {0, . . . , 3k = 9} of the points in V x j and V x j-1 . For this indexing, the subdivision matrix reads

                7 3 3 3 1 1 1 1 10 1 1 1 1 1 1 1 10 1 1 1 1 1 1 1 10 1 1 3 3 1 1 3 3 1 1 3 3 1 3 1 1 1 1 3 1 1 2 1 3                
where the 0's have been omitted and where the rows should be rescaled to sum to 1. The eigenvalues of this matrix satisfy λ 1 = 1 and λ 2 = λ 3 = 1/3 > λ 4 .

Wavelets on Meshes

Multiscale Biorthogonal Bases on Meshes

The transforms considered in this section are multiscale and indexed by the set of nested grids (V j ) L<j J . This corresponds to computing a set of coefficients (d j ) L<j J ∪ f J from an initial input signal f . These coefficients corresponds to inner products with basis vectors

d j ∈ ℓ 2 (H j ) where ∀ k ∈ H j , d j (k) = f, ψ j,k , f J ∈ ℓ 2 (V J ) where ∀ ℓ ∈ V J , f J (ℓ) = f, ϕ J,ℓ .
By analogy with the wavelet setting, the vectors ψ j,k ∈ R n corresponds to primal wavelets and are intended to capture the details present in the signal f at a scale j, whereas the scaling vectors ϕ J,k ∈ R n capture the missing coarse approximation of f at scale J. This decomposition is stopped at any coarse scale L < J 0.

In order to reconstruct the function f from this set of transformed coefficients, one needs to use a set of bi-orthogonal basis vectors

f = L<j J,k∈Hj d j (k) ψj,k + ℓ∈V J f J (ℓ) φJ,ℓ .
If this reconstruction formula holds for any scale L < J 0, the set of vectors (ψ j,k , ϕ j,ℓ ) L<j 0 k∈Hj ,ℓ∈Vj and ( ψj,k , φj,ℓ ) L<j 0 k∈Hj ,ℓ∈Vj ,

is said to be a pair of primal and dual multiscale bases (together with their scaling functions).

The following paragraph shows how one can modify such a pair of multiscale bases while still maintaining the biorthogonality property. This lifting process is useful to design multiscale bases with various properties on complicated domains.

The Lifting Scheme

The lifting scheme is a construction of multiscale biorthogonal bases introduced by Sweldens [START_REF] Sweldens | The lifting scheme: A custom-design construction of biorthogonal wavelets[END_REF][START_REF] Sweldens | The lifting scheme: A construction of second generation wavelets[END_REF]. It extends the traditional construction of wavelets in two main directions:

As explained in [START_REF] Daubechies | Factoring wavelet transforms into lifting steps[END_REF], it allows to implement already existing filter banks more efficiency by splitting the computation into elementary blocks. This computational gain is described at the end of the section together with the factorization of wavelets into lifting steps. It allows to define multiscale transforms over domains that are not translation invariant. This section gives two examples of such transforms: a non-separable 2D wavelet transform and wavelets on triangulated meshes.

In order to build wavelets on triangulation, one can specialize the lifting scheme to a particular setting where only two lifting steps are applied.

Forward lifting scheme. The forward algorithm performs the transform (f j-1 (ℓ)) ℓ∈Vj-1 -→ (d j (k)) k∈Hj ∪ (f j (ℓ)) ℓ∈Vj by applying the following steps Splitting: this corresponds selecting the coefficient of f j-1 (ℓ) that are in V j or in H j (f j-1 (ℓ)) ℓ∈Vj-1 = (f j (ℓ)) ℓ∈Vj ∪ (f j (ℓ)) ℓ∈Hj .

These two sets of coefficients are treated differently in the two remaining steps of the transform. Predict step: creates wavelets coefficients d j by computing local differences between each coefficient in V j and its neighbors in H j ∀ k ∈ H j , d j (k) = f j-1 (k) -ℓ∈Vj p j (k, ℓ)f j-1 (ℓ).

The coefficients p j (k, ℓ) are weights that determine the predict operator P j : ℓ 2 (V j ) -→ ℓ 2 (H j ) g -→ h = P j g where h(k) = k∈Hj p j (k, ℓ)g(ℓ).

Update step: enhance the properties of each remaining low pass coefficients f j-1 (ℓ) for ℓ ∈ V j by pooling locally the wavelets coefficients d j (k) for k around ℓ ∀ ℓ ∈ V j , f j (ℓ) = f j-1 (ℓ) + k∈Hj u j (ℓ, k)d j (k).

The coefficients u j (ℓ, k) are weights that determine the update operator

U j : ℓ 2 (H j ) -→ ℓ 2 (V j ) h -→ g = U j h
where g(ℓ) = k∈Hj u j (ℓ, k)h(k). One of the main feature of the lifting scheme is that this is achieved by simply reversing the order of the lifting steps and interchanging +/-signs.

Inverse update step:

∀ ℓ ∈ V j , f j-1 (ℓ) = f j (ℓ) -k∈Hj u j (ℓ, k)d j (k).

Inverse predict step:

∀ k ∈ H j , f j-1 (k) = d j (k) + ℓ∈Vj p j (k, ℓ)f j-1 (ℓ).

Merging: makes the union of the coefficients computed in the two previous steps (f j-1 (ℓ)) ℓ∈Vj-1 = (f j (ℓ)) ℓ∈Vj ∪ (f j (ℓ)) ℓ∈Hj . The lifting scheme is more general than the algorithm described in this section since several passes of predict/update steps can be applied to further enhance the properties of the resulting transform. However, the steps beyond the two initial ones are difficult to analyze, except in the notable exception of points sampled evenly on a 1D axes, where a factorization algorithm [START_REF] Daubechies | Factoring wavelet transforms into lifting steps[END_REF] allows to recover traditional wavelet filters.

Imposing vanishing moments.

The operator P j is called a predictor since the values of P j f Vj j-1 should typically be close to f Hj j-1 for the wavelet coefficients d j to be small. Such predictors have already been constructed in equations (3.2), (3.3) and (3.3.2).

The operator U j is called an update operator since the additional term U j d j should enhance the properties of f Vj j-1 . This update steps does not appears in the theory of subdivision surface and this section considers a local update operator which guaranty the conservation of the mean value when switching from f j-1 to f j .

Polynomial vectors. In order to select predict and update operator that have good properties, one follow the insight gained from the analysis of the wavelet approximation of signal on the real line. In order to do so, one need analyze the effect of a lifting wavelet transform on polynomials. The most basic constraint enforces one vanishing moment by imposing orthogonality with the constant vector Φ 0 = 1. This constraint does not require to known the spacial location x ℓ of each index ℓ ∈ V L . In order to impose higher order vanishing moments, one needs to assume some sampling pattern, for instance ∀ ℓ ∈ V L , f (ℓ) = f (x ℓ ) where x ℓ ∈ R q and where f is a function defined on R q . For instance, the points x ℓ might corresponds to a regular sampling of the line (this is the traditional wavelet setting) or to an irregular sampling of a 2D surface embedded in R 3 . The next paragraphs describe several situations with different sampling grids. Once the precise locations of the samples are known, one can for instance select Φ s as some monomials of degree (s 1 , . . . , s q ) over R q .

Vanishing moment and polynomials reproduction. Having defined these polynomial vectors, one requires that the following constraints are fulfilled.

Vanishing moments: the wavelet coefficients of a low order polynomial should be 0, which implies that ∀ k ∈ H j , Φ s , ψ j,ℓ = 0. (4.7)

Polynomial reproduction: coarse coefficients f j computed from a polynomial f f -1 should also be polynomials, which implies that ∀ ℓ ∈ V j , Φ s , ϕ j,ℓ = Φ s (ℓ) (4.8)

In order for the wavelets and scaling function to satisfy conditions (3.7) and (3.8), the predict operator P j and update operator U j should be designed carefully. One can impose these constraint from the fine scale j = L until the coarse scale j = 0. Indeed, if (ϕ j-1,ℓ , ψ j-1,ℓ ) k,ℓ satisfy conditions (3.7) and (3.8), then, for the scale j where Φ A s ∈ ℓ 2 (A) is the restriction of Φ s to A. In contrast, the constraint (3.8) on the update operator P j is more involved and the next section shows how to handle it on a triangulation situations for only one vanishing moment |S| = 1.

Lifted Wavelets on Meshes

The lifted wavelet bases can be used to process signals f ∈ ℓ 2 (V L ) where ℓ ∈ V L index a sampling x ℓ of an arbitrary surface. The construction of biorthogonal wavelets on triangulated mesh has been first proposed by Lounsbery et al. [START_REF] Lounsbery | Multiresolution analysis for surfaces of arbitrary topological type[END_REF] and re-casted into the lifting scheme framework by Schroeder and Sweldens [START_REF] Schröder | Spherical Wavelets: Efficiently Representing Functions on the Sphere[END_REF][START_REF] Schröder | Spherical wavelets: Texture processing[END_REF].

Designing predict operators. The constraints (3.7) on the predictor P j is easily solved. For instance, for each k, one selects only |S| non vanishing weights (p j (k, ℓ)) ℓ and solves a small |S| × |S| linear system. Furthermore, in the case of a regular triangulation with edges of constant length, predictors with several vanishing moments have been already defined in (3.2), (3.3) and (3.3.2). Figure 3.14 shows the weights for these predictors.

Examples of Update Operators

]) erators: λ.

24.

P j 1 = 1. One can choose any of these operators, and creates respectively linear, butterfly and Loop wavelets bases. All these predictors have one vanishing moment since they satisfy P j 1 Hj = 1 Vj . In fact they have more vanishing moments if one consider polynomials Φ s sampled at points x ℓ ∈ R 2 of an hexagonal tiling with constant edge length. In practice, if the triangulation under consideration have edges with smoothly varying length, the resulting predictor are efficient to predict the value of smooth functions on the triangulation.

Designing update operators. In order to ensure the reproduction of constant polynomials, we design the update operator so that it depends only on the direct neighbors in H j of each point in V j ∀ ℓ ∈ V j , V ℓ = {γ(ℓ, ℓ ′ ) \ (ℓ, ℓ ′ ) ∈ E j } .

One then looks for a valid update operator in the following form ∀ h ∈ ℓ 2 (H j ), ∀ ℓ ∈ V j , (U j h)(ℓ) = λ ℓ k∈V ℓ h(k), (4.9)

where each λ ℓ should be fixed in order for condition (3.8) to be satisfied.

In an semi-regular triangulation, |V ℓ | = 6 excepted maybe for some points in the coarse grid ℓ ∈ V 0 . In this setting, the values of λ ℓ can be computed by a recursion through the scales. In an ideal triangulation where |V ℓ | = 6 for all ℓ, one can use a constant weight λ ℓ = λ.

Non-linear Mesh Compression

These wavelets can be used to perform an approximation of a function f ∈ ℓ 2 (V L ) defined on the fine triangulation. For instance a wavelet approximation can be applied to each coordinate f i , i = 1, 2, 3 of the actual position x ℓ = (f 1 (ℓ), f 2 (ℓ), f 3 (ℓ)) ∈ R 3 of the surface points, as done in [START_REF] Guskov | Multiresolution signal processing for meshes[END_REF][START_REF] Khodakovsky | Progressive geometry compression[END_REF]. This leads to a scheme to approximate and compress a 3D surface using the lifted biorthogonal wavelets associated to the semi-regular triangulation. This is possible because these wavelets depend only on the combinatorial grids V j and not on the precise position of the samples x ℓ in 3D. Note that for each coefficient the threshold T is scaled according to the size of the support of the wavelet in order to approximately normalize the wavelets in ℓ 2 (V L ) norm.

Figure 1 . 1 :

 11 Figure 1.1: Example of display of a 3D mesh.

Figure 1 . 2 :

 12 Figure 1.2: One ring around a vertex i, together with the geometrical angles α ij and β ij used to compute the conformal weights.

Figure 1 . 3 :Figure 1 . 4 :

 1314 Figure 1.3: Examples of iterative smoothing of a 3D mesh. Definition 11 (Edges Voronoi). The Voronoi diagram associated to the edges is ∀ e = (i, j) ∈ E, E e = {x ∈ M \ ∀ e ′ = e, d(x, e) d(x, e ′ )}

Figure 1 .

 1 [START_REF] Clarkson | Building triangulations using epsilon-nets[END_REF] shows an example of such a smoothing.

Figure 1 . 5 :

 15 Examples of mesh denoising with the heat equation.

Figure 1 .

 1 Figure 1.6 shows examples of the resolution of the wave equation on 3D meshes.

Figure 1 . 6 :

 16 Figure 1.6: Example of evolution of the wave equation on 3D mesh. The initial condition f is a superposition of small positive and negative gaussians.

Figure 1 . 7 :

 17 Figure 1.7: Examples of eigenvectors u ω of the Laplacian L. The blue colors indicated negative values, red colors positive ones. The black curve is the 0 level set of the eigenvector.

Figure 1 . 7

 17 Figure 1.7 shows some examples of eigenfunctions depicted using color ranging from blue (negative values of the eigenfunction) to red (positive values). One can see that these functions are oscillating, in a way similar to the traditional Fourier basis. In some sense (made more precise latter), this basis is the extension of the Fourier basis to meshes. A function u ω corresponding to a large spectral value λ ω is highly oscillating and corresponds thus intuitively to a high frequency atom.

Figure 1 . 8 :

 18 Figure 1.8: Examples of Fourier spectrum for a smooth and a non-smooth function.

Figure 1 .

 1 Figure 1.10 shows several examples of parameterizations. One is free to use any laplacian (combinatorial, distance or conformal) as long as it produces positive weights. There is a issue with the conformal weights, which can be negative if the mesh contains obtuse triangles. In practice however it leads to the best results. The efficiency of a parameterization can be measured by some amount of distortion induced by the planar mapping. Linear methods cannot hope to cope with large isoperimetric distortions (for instance large extrusions in the mesh) since harmonicity leads to clustering of vertices.

Figure 1 .Figure 1 . 10 :

 1110 [START_REF] Eldar | The farthest point strategy for progressive image sampling[END_REF] shows that for meshes with large distortion, this flattening indeed leads to wrong parameterizations. Examples of mesh parameterizations.MeshCombinatorial Conformal Isomap Figure 1.11: Examples of mesh flattening.

Example 7 .

 7 Figure 2.1 shows some some frequently used geodesic metric spaces: Euclidean space: M = R s and H(x) = Id s . 2D shape: M ⊂ R 2 and H(x) = Id 2 . Parametric surface: H(x) = I ϕ (x) is the first fundamental form. Isotropic metric: H(x) = W (x)Id s , W (x) > 0 being some weight function. Image processing: given an image I : [0, 1] 2 → R, one can use an edge-stopping weight W (x) = (ε + ||∇ x I||) -1. This way, geodesic curves can be used to perform segmentation since they will not cross boundaries of the objects. fMRI imaging: M = [0, 1] 3 , and H(x) is a field of diffusion tensors acquired during a scanning experiment. For fMRI imaging, the direction field e 1 indicates the direction of elongated fibers of the white matter.

Figure 2 . 1 :

 21 Figure 2.1: Examples of metric (top row) and geodesics (bottom row). From left to right: euclidean, isotropic, shape, surface and 3D Riemannian manifold metrics (diffusion tensor). The blue/red colormap indicates the geodesic distance to the starting point.

2 .

 2 Select point: x ←-argmin S(z)=T rial P(z).

3 . 4 .

 34 Tag: S(x) ← Computed. Update neighbors: for all y ∈ N eigh(x), -If S(y) = F ar, then S(y) ← F ront and U (y) ← U pdate(y).

Figure 2 .

 2 Figure 2.3 shows examples of distance functions to a starting point x 0 with the corresponding geodesics γ(t) extracted from some ending point x 1 . The front propagation is stopped when

Figure 2 . 2 :

 22 Figure 2.2: Examples of front propagation. The colormap indicates the values of the distance functions at a given iteration of the algorithm.

Figure 2 . 3 ,

 23 Figure 2.3, left and middle, shows examples of such curves extractions.Gradient avoiding potential: for edge detection, one wants geodesics that avoid traversing image edges. In this case, the potential should be small in area of large gradient in the image

Figure 2 . 3 ,

 23 Figure 2.3, right, shows an example of such a road extraction.

Figure 2 . 3 :

 23 Figure 2.3: Example of distance functions (top row) and geodesics (bottom row).way (applied in figure2.5 where one can see in red the front of the Fast Marching propagation) is to display various isosurfaces which are level sets {(x, y, z) \ f (x, y, z) = c} for some c.

Figure 2 . 4 :

 24 Figure 2.4: Example of volume data visualization with slices (top row) and semi-transparent mapping (bottom row).

Figure 2 .

 2 Figure 2.5 shows some examples of geodesic extraction on a medical image that represents tubular structures (blood vessels) around the heart. The potential W (x) is chosen as W(x) = (|f (x) -f (x 0 )| + ε) -1where x 0 is a point given by the user and supposed to lie inside some vessel. A geodesic follows nicely a vessel since its density is constant and thus the value of f is approximately equal to f (x 0 ) inside the vessel.

Figure 2 . 5 :

 25 Figure 2.5: Example of volumetric Fast Marching evolution (top row) and geodesic extractions (bottom row).

Figure 2 . 6 :

 26 Figure 2.6: Example of Fast Marching propagation on a triangulated mesh.

Figure 2 .

 2 Figure 2.6 shows an example of propagation on a triangulated surface. The colored region corresponds to the point that are computed (its boundary being the front).

Figure 2 . 7 :

 27 Figure 2.7: Examples of geodesic extraction on a mesh with an increasing number of starting points.

  Figure 2.7 shows examples of such distances to several points together with geodesics. A geodesic γ links a point x to its closest point in S.

Definition 34 (

 34 Geodesic distance in S). The geodesic distance in S for the uniform metric is d S (x, y)

Figure 2 . 8 :

 28 Figure 2.8: Geodesics inside a 2D shape.

Figure 2 .

 2 Figure 2.8 shows examples of shapes together with the geodesic distance to a starting point. The geodesic curve is the union of segments inside S and pieces of the boundary ∂S.

Figure 2 . 9 (

 29 Figure 2.9 (top row) shows several examples of eccentricity. The colormap indicates in blue points with small eccentricity. The points for which the minimum in the definition of E S is obtained are called eccentric. Definition 36 (Eccentric points). An eccentric point x ∈ E(S) satisfies ∃y, E S (y) = d(x, y). These eccentric points define regions of influence which perform a segmentation of the shape as follow S = x∈E(S) {y ∈ S \ E S (y) = d(x, y)}

Figure 2 . 9 :

 29 Figure 2.9: Example of eccentricity E S and corresponding histograms h S .

1 .

 1 Dataset:shapes {S 1 , . . . , S p } (binary images). 2. Preprocessing: compute eccentricity descriptors h Si . 3. Input: shape S. 4. Retrival: return i ⋆ = argmin i δ(h S , h Si ).

Table 2 :

 2 Shape retrieval process.

Figure 2 . 10 :

 210 Figure 2.10: Examples of shape recognitions. The shape on the left is the input S and the second shape in each raw is S i ⋆ .

Figure 2 . 11 :

 211 Figure 2.11: Example of propagations with a priority P(x) = U (x) + λV (x) for λ = 0, 0.5, 0.9.

Figure 2 .

 2 Figure 2.11shows examples of heuristics that approximate the true remaining distance by bellow. One can see how the explored area of the propagation progressively shrinks while containing the true geodesic. Such an heuristic is however impossible to use in practice since one does not have direct access to the remaining distance during the propagation.

Theorem 29 (

 29 Reversed triangular inequality). For all (x, y) ∈ M, one has d(x, y) = sup z |d(x, z) -d(z, y)| . If one restricts the minimum to a small subset of landmark points {z 1 , . . . , z n } ⊂ M, one can define the following approximate distance dz1...zn (x, y) = sup k=1...n |d k (x) -d k (y)| ,

Figure 2 . 12 :

 212 Figure 2.12: Heuristically driven propagation in 2D with an increasing number of landmark points.

Figure 2 . 13 :

 213 Figure 2.13: Heuristically driven propagation on a 3D mesh with landmark points.

Figure 2 . 14 :

 214 Figure 2.14: Examples of farthest point sampling (the colormap indicates the distance function to the seeds).

Figure 2 .

 2 Figure 2.14 shows examples of farthest point sampling with an uniform (top row) and a spatially varying isotropic metric W (x) (bottom row). One can see that this scheme seeds more points in areas where the metric W is large. One can thus control the sampling density by modifying the metric W .

Sampling Voronoi Delaunay Figure 2 . 15 :

 215 Examples of sampling and triangulations.

Figure 2 .

 2 Figure 2.16 shows examples of Voronoi cells on a surface embedded in R 3 .

Figure 2 .

 2 Figure 2.18 shows an example of uniform remeshing of a 3D surface with an increasing number of points. Figure2.19 shows how one can adapt the density by defining a non-constant isotropic metric on the surface.

Figure 2 .

 2 Figure 2.18 shows an example of uniform remeshing of a 3D surface with an increasing number of points. Figure2.19 shows how one can adapt the density by defining a non-constant isotropic metric on the surface.

Figure 2 . 16 :

 216 Figure 2.16: Example of Voronoi segmentations for an increasing number of seeding points.

Figure 2 . 17 :

 217 Figure 2.17: Shape meshing with an increasing number of points. Left and center: uniform meshing, right: adaptive meshing.

Figure 2 . 18 :

 218 Figure 2.18: Geodesic remeshing with an increasing number of points. Top: sampling with corresponding distance function, bottom: corresponding geodesic triangulation.

Figure 2 . 19 :

 219 Figure 2.19: Adaptive remeshing with a binary density function (top row) and a linearly increasing density (bottom row).

Figure 2 . 20 :

 220 Figure 2.20: Adaptive remeshing with a density given by a texture. A gradient descent of E i allows to compute the geodesic gravity center of a given region of the surface. It is the central ingredient of the Lloyd algorithm detailed in pseudo code 4, that converges in practice to a local optimum of the energy E and thus a centroidal tesselation. Of course, centroidal tesselation of a given domain are not unique and computing the best one (with the smallest value of E) is a difficult task. In a discrete setting (point clouds in the Euclidean space) the Lloyd algorithm is equivalent to the k-means algorithm used to perform vector quantization and classification. 1. Initialize: {x i } i ← random.

Figure 2 .

 2 Figure 2.21, top row, shows examples of iterations of the Lloyd algorithm on a square with euclidean metric. Middle and bottom row shows the Lloyd algorithm with a weighted euclidean metric, for which the gravity center are defined as

Figure 2 .

 2 Figure 2.22 shows iterations of the algorithm on a 3D surface, equipped with an uniform metric.

Figure 2 . 21 :

 221 Figure 2.21: Iterations of the Lloyd algorithm on a square with euclidean metric (top row) and weighted euclidean metric (middle and bottom rows). The Delaunay triangulation is depicted in red and the corresponding Voronoi diagram is in blue.

Figure 2 . 22 :

 222 Figure 2.22: Iterations of the Lloyd algorithm on two different surfaces.

3. 1 . 1

 11 Point Clouds and Manifold 3.1.2 Graph Approximation and Geodesic Distances 3.1.3 Image Ensembles

D∈Dn

  ||D -D||.This is a difficult optimization task, and in order to simplify it, this optimization can be approximated by replacing the distance matrices by gram matrices as follow min K∈Kn ||K -ψ( D)||.

Figure 3 . 1 :

 31 Figure 3.1: Top: original surfaces, bottom: examples of bending invariants.

4. 1 . 1

 11 Nested Multiscale Grids.

Figure 4 . 1 :

 41 Figure 4.1: Edge-splitting subdivision.

Figure 4 . 2 :

 42 Regular subdivision 1:4 of a single triangle. Regular subdivision of a planar triangulation M 0 .

Figure 3 .

 3 Figure 3.1 shows the notations related to the subdivision process.Figure 3.2 shows an example of recursive splitting of a triangle and a coarse triangulation.Figure 3.3 shows examples of semiregular triangulation using a geometric realization (position of the vertices) to create a 3D surface.

Figure 3 .

 3 Figure 3.1 shows the notations related to the subdivision process.Figure 3.2 shows an example of recursive splitting of a triangle and a coarse triangulation.Figure 3.3 shows examples of semiregular triangulation using a geometric realization (position of the vertices) to create a 3D surface.

Figure 3 .

 3 Figure 3.1 shows the notations related to the subdivision process.Figure 3.2 shows an example of recursive splitting of a triangle and a coarse triangulation.Figure 3.3 shows examples of semiregular triangulation using a geometric realization (position of the vertices) to create a 3D surface.

4. 1 3 Figure 4 . 3 :

 1343 Figure 4.3: Examples of semi-regular meshes (V j ) j for increasing scale j (from left to right).

"

  #$%$&'(!!!!!! )*+,#$-'(!*'#'.,/#$0'/$"&! "-/'+,1#'(!*'#'.,/#$0'/$"&! %,".,/#2!$.'%,!3($/4! #,.,)+,1!%,".,/#2!!! 5$%6#,!78!9,."&)/#'/$"&!":!)*+,#$-'(!*'#'.,/#$0'/$"&!'&1!)6;),<6,&/!#,)'.*($&%!$&/"!'!%,".,/#2!$.'%,=!

Figure 4 . 4 :

 44 Figure 4.4: Spherical geometry image construction, taken from [31].

Figure 4 . 5 :

 45 Figure 4.5: 1D subdivision scheme with filters h and h. The red curve represent the original signal f 0 .

Figure 4 . 6 :

 46 Figure 4.6: 1D subdivision of a signal. Bottom row shows the subdivision from an impulse signal, converging to the scaling function ϕ.

Figure 3 .

 3 Figure 3.6 shows several steps of subdivision, starting from an initial vector of size |V 0 | = 10.

Figure 3 .

 3 Figure 3.7 shows examples of subdivision curves.

Figure 4 . 7 :

 47 Figure 4.7: Two examples of subdivision curves. The red curve represent the original curve (X 0 , Y 0 ).

Figure 4 . 9 :Figure 4 . 10 :

 49410 Figure 4.9: Examples of iterative subdivision using Loop scheme. The points (X 0 , Y 0 , Z 0 ) of the initial coarse mesh M 0 are shown in red.

Figure 4 . 11 :

 411 Figure 4.11: Surface after 0, 1 and 3 step of √ 3 subdivision [20].

Figure 3 .

 3 Figure 3.12 shows an example of invariant neighborhood which corresponds to the 2-ring V

Figure 4

 4 Figure 4.12: Invariant neighborhood V x j and V x j-1 (indexing with red circles) of the Loop subdivision scheme for a vertex of valence |V ℓ | = 0. The number in {0, . . . , 9} refers to the numbering of the vertices in V x j and V x

Figure 3 .

 3 Figure 3.13, top row, shows the block diagram associated to this forward lifting wavelet transform.The iterations of the forward lifting transform can also be written in vector and operator formatd j = f Hj j-1 -P j f Vj j-1 , f j = f Vj j-1 + U j d j = (Id Vj -U j P j )f Vj j-1 + U j f Hj j-1 ,

)

Figure 4 . 13 :

 413 Figure 4.13: Block diagrams for the forward and backward lifting scheme.

Figure 3 .

 3 Figure 3.13, top row, shows the block diagram associated to this backward lifting wavelet transform.

∀ s ∈ S, ( 3 . 7 )

 37 ⇐⇒ P j Φ

Figure 4 . 14 :

 414 Figure 4.14: Predict operators on a triangulation.

For

  the predictors defined in (3.2), (3.3) and (3.3.2), one has P j T 1 Hj = 3 × 1 Vj and U j T 1 Vj = 6λ1 Hj so setting λ ℓ = 1/24 solves equation (3.8). Figure 3.15 shows examples of butterfly wavelets on a planar semi-regular triangulation.

j = - 2 j = - 3 j = - 4 Figure 4 . 15 :

 234415 Example of wavelets ψ j,k on a semi-regular triangulation. The height over the triangle (together with the color) indicates the value of the wavelet vector.In order to perform a wavelet approximation in this biorthogonal basis, one uses a non-linear thresholding atT > 0 f = (j,k)∈I T f, ψ j,k ψj,kwhereI T = (j, k) \ k ∈ H j and | f, ψ j,k | > T | supp(ψ j,k )| -1/2 .

Figure 3 .

 3 Figure 3.16 shows an example of compression of the position of a vertex in 3D spaces as 3 functions defined on a semi-regular mesh. Figure3.17 shows an example of compression of a spherical texture map which is a single function defined at each vertex of a semi-regular mesh obtained by subdividing an icosaedron.

Figure 3 .Figure 4 . 16 :

 3416 Figure 3.16 shows an example of compression of the position of a vertex in 3D spaces as 3 functions defined on a semi-regular mesh. Figure3.17 shows an example of compression of a spherical texture map which is a single function defined at each vertex of a semi-regular mesh obtained by subdividing an icosaedron.

Table 4 :

 4 Max-Lloyd relaxation algorithm.

100%

10% 5% 2% Figure 4.17: Non-linear spherical wavelet compression with a decreasing number of coefficients.