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 and toolbox_graph 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.

Fast Marching

The Dijstra algorithm suffers from a strong metrization problem, and it actually computes the \(\ell^1\) distance on the grid.

The Fast Marching algorithm replace the graph update by a local resolution of the Eikonal equation. This reduces significantly
the grid bias, and can be shown to converge to the underlying geodesic distance when the grid step size tends to zero.

Exercice 3: (check the solution) Compute the distance map to these starting point using the FM algorithm.

exo3;

Once the geodesic distance map to \(\Ss\) has been computed, the geodesic curve between any point \(x_1\) and \(\Ss\) extracted
through gradient descent \[ \ga'(t) = - \eta_t \nabla D(\ga(t)) \qandq \ga(0)=x_1 \] where \(\eta_t>0\) controls the parameterization
speed of the resulting curve.

To obtain unit speed parameterization, one can use \(\eta_t = \norm{\nabla D(\ga(t))}^{-1}\) (one need to be careful when
\(\ga\) approaches \(\Ss\) since \(D\) is not smooth on \(\Ss\)).

Normalize the gradient to obtained \(G(x) = G_0(x)/\norm{G_0(x)}\), in order to have unit speed geodesic curve (parameterized
by arc length).

G = G0 ./ repmat( sqrt( sum(G0.^2, 3) ), [1 1 2]);

The geodesic is then numerically computed using a discretized gradient descent, which defines a discret curve \( (\ga_k)_k
\) using \[ \ga_{k+1} = \ga_k - \tau G(\ga_k) \] where \(\ga_k \in \RR^2\) is an approximation of \(\ga(t)\) at time \(t=k\tau\),
and the step size \(\tau>0\) should be small enough.

Step size \(\tau\) for the gradient descent.

tau = .8;

Initialize the path with the ending point.

x1 = round([.9;.88]*n);
gamma = x1;

Define a shortcut to interpolate \(G\) at a 2-D points. Warning: the interp2 switches the role of the axis ...

Compute the gradient at the last point in the path, using interpolation.

g = Geval(G, gamma(:,end));

Perform the descent and add the new point to the path.

gamma(:,end+1) = gamma(:,end) - tau*g;

Exercice 4: (check the solution) Perform the full geodesic path extraction by iterating the gradient descent. You must be very careful when the path become
close to \(x_0\), because the distance function is not differentiable at this point. You must stop the iteration when the
path is close to \(x_0\).