Installing toolboxes and setting up the path.

You need to unzip these toolboxes in your working directory, so that you have toolbox_signal, toolbox_general, toolbox_graph and toolbox_wavelet_meshes in your directory.

For Scilab user: you must replace the Matlab comment '%' by its Scilab counterpart '//'.

Recommandation: You should create a text file named for instance numericaltour.sce (in Scilab) or numericaltour.m (in Matlab) to write all the Scilab/Matlab command you want to execute. Then, simply run exec('numericaltour.sce'); (in Scilab) or numericaltour; (in Matlab) to run the commands.

Triangulated Mesh

We consider a 3-D surface discretized using a triangular mesh with \(N\) vertices. We denote \( \{x_i\}_{i=1}^{N} \) the set
of vertices, where \(x_i \in \RR^3\). We denote \(\Ff \subset \{1, \ldots, N\}^3\) the set of faces.

Control Formulation of the Eikonal Equation

Our goal is to compute the geodesic distance to a set of starting points.

We consider an isotropic metric on the mesh \(W_i > 0\).

Given a set \(I\) of starting point, we want to compute the discrete geodesic distance map \( U \in \RR^N \) to the starting
points \( \{x_i\}_{i \in I} \).

Let's consider a geodesic distance \(d(x,y)\) on a manifold. A general mathematical result is that the geodesic distance \[
d(x,A) = \umin{a \in A} d(x,a) \] to a set \(A\) is the unique solution to the following equation \[ \forall \, x \notin A,
\quad d(x,A) = \umin{ y \in B(x) } d(y,A) + d(x,y) \] where \( B(x) \) is a small disk around \(x\) that does not contain
\(A\).

In a discrete setting, one can use this property with \(B(x_i)\) being the one ring, made of the edges connecting vertices
that are around \(x_i\). Using a piecewise linear interpolation of \(U\) on the mesh, one obtain that \(U\) is a solution
to a fixed point equation \[ U = \Ga(U) \] where the update operator \( \Ga : \RR^N \rightarrow \RR^N \) is defined as \[
\Ga(U)_i = \umin{ (i,j,k) \in \Ff } d_{i,j,k} \] where \[ d_{i,j,k} = \umin{0 \leq t \leq 1} t U_{j} + (1-t) U_k + W_i \norm{
t x_j + (1-t) x_k - x }. \] The quantities \( d_{i,j,k} \) can in fact be computed by solving a second order polynomial,
as we describe in the following section.

Geodesic distance Computation Using Jacobi Iterations

The mapping \(\Ga\) is monotone \[ U \leq V \qarrq \Ga(U) \leq \Ga(V). \] If one uses an initialization \(U^{(0)}\) such that
\(\Ga(U^{(0)}) \geq U^{(0)} \) (for instance setting \(U^{(0)}=0\)), then the iterates \[ U^{(\ell+1)} = \Ga(U^{(\ell)}) \]
are increasing. One can prove that they are bounded, and hence they converge to the (unique) fixed point solution of \(\Ga(U)=U\).

Define the metric \(W_i > 0\) on the mesh. We start with the uniform (constant) metric.

W = ones(n,1);

Indexes \(I \subset \{1,\ldots,N\}\) of starting points. Warning: this list works for the mesh nefertitit refined twice. For other meshes, use the function select3dtool to retrive the indexes of a few vertices.

To ensure correctness of the scheme, the update should come from within the triangle (which correspond to the condition \[t
\in [0,1]\] in the original definition of \(d_{i,j,k}\) using the control formulation). This corresponds to having \( \al_1<=0
\) and \(\al_2 \leq 0\) where \( \al = S ( u - d \mathbb{I} ) \).

Compute \(\al \in \RR^2\).

alpha = Mult( S, u - repmat(d, 2, 1) );

For the index \(i \in J\) where the update comes from outside the edge \( [x_j,x_k] \), we use an update along the edge \[
d = \min\pa{ U_j + \norm{x-x_j}W_i, U_k + \norm{x-x_k}W_i }. \]

Now that d stores the value of d_{i,j,k} for a bunch of faces assign in U1(i) the mininimum between the previous U1(i) and all the entries in d that corresponds to a face \((i,j,k)\).

U1 = accumarray(i', d, [n 1], @min); U1(U1==0) = Inf;

Enforce the boundary conditions \(\forall i \in I, \, U_i = 0\).

U1(I) = 0;

Update the solution.

U = U1;

Exercice 1: (check the solution) Perform the full iterative algorithm until convergence. Store in err(l) the fixed point error \( \norm{ U^{(\ell+1)} - U^{(\ell)} } \). Note: you might want to put outside of the loop all the quantities that do not depend on \(u\), e.g. \(S, a, \) etc.

Geodesic Path Extraction

A geodesic minimal path from any point on the mesh to the starting point can be computed using a gradient descent. Since the
computed geodesic distance is affine on each triangle, it makes sense to define a discretized path that is a line segment
on each face.

Extract minimal paths from all over the mesh, starting from the boundary.

Exercice 3: (check the solution) Display the convergence of the computed geodesic distance to the the true geodesic distance (which is the Euclidean distance
\( \norm{x_i} \)) as \(n\) increases. Note: the triangulation with increasing number of points should be refining (i.e. a finer triangulation should contains all the
other ones).