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.

Note that \(L(\gamma)\) is invariant under re-parameterization of the curve \(\gamma\).

A geodesic curve \(\gamma\) between two points \(x_0\) and \(x_1\) has minimum length among curves joining \(x_0\) and \(x_1\),
\[ \umin{\ga(0)=x_0, \ga(1)=x_1} L(\ga). \] A shortest curve thus tends to pass in areas where \(W\) is small.

The geodesic distance between the two points is then \(d(x_0,x_1)=L(\gamma)\) is the geodesic distance according to the metric
\(W\).

Define start and end points \(x_0\) and \(x_1\) (note that you can use your own points).

x0 = [14;161];
x1 = [293;148];

The metric is defined according to \(f\) in order to be low at pixel whose value is close to \(f(x)\). A typical example is
\[ W(x) = \epsilon + \abs{f(x_0)-f(x)} \] where the value of \( \epsilon>0 \) should be increased in order to obtain smoother
paths.

epsilon = 1e-2;
W = epsilon + abs(f-f(x0(1),x0(2)));

Display the metric \(W\).

clf;
imageplot(W);

Set options for the propagation: infinite number of iterations, and stop when the front hits the end point.

options.nb_iter_max = Inf;
options.end_points = x1;

Perform the propagation, so that \(D(a,b)\) is the geodesic distance between the pixel \(x_1=(a,b)\) and the starting point
\(x_0\). Note that the function perform_fast_marching takes as input the inverse of the metric \(1/W(x)\).

[D,S] = perform_fast_marching(1./W, x0, options);

Display the propagated distance map \(D\). We display in color the distance map in areas where the front has propagated, and
leave in black and white the area where the front did not propagate.

Exercice 1: (check the solution) Using options.nb_iter_max, display the progressive propagation. This corresponds to displaying the front \( \enscond{x}{D(x) \leq t} \) for various
arrival times \(t\).

exo1;

Geodesic Curve Extraction

Once the geodesic distance map \(D(x)\) to a starting point \(x_0\) is computed, the geodesic curve between any point \(x_1\)
and \(x_0\) extracted through gradient descent \[ \ga'(t) = - \eta_t \nabla D(\ga(t)), \] 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}\).

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]);

Display \(G\).

clf;
imageplot(G);
colormap jet(256);

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.

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 2: (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\).

Exercice 4: (check the solution) Perform the shortest path extraction for various images such as 'cavern' or 'mountain'.

exo4;

Edge-based Geodesic Metric

It is possible to extract the boundary of an object using shortest paths that follows region of high gradient.

First we load an image \(f\).

n = 256;
name = 'cortex';
f = rescale( sum(load_image(name,n),3) );

Display it.

clf;
imageplot(f);

An edge-attracting potential \(W(x)\) should be small in regions of high gradient. A popular choice is \[ W(x) = \frac{1}{\epsilon
+ G_\si \star G(x)} \qwhereq G(x) = \norm{\nabla f(x)}, \] and where \(G_\si\) is a Gaussian kernel of variance \(\si^2\).

Compute the gradient norm \(G(x)\).

G = grad(f,options);
G = sqrt( sum(G.^2,3) );

Smooth it by \(G_\si\).

sigma = 3;
Gh = perform_blurring(G,sigma);

Display the smoothed gradient \( G \star G_\si \).

clf;
imageplot(Gh);

Compute the metric.

epsilon = 0.01;
W = 1./( epsilon + Gh );

Display it.

clf;
imageplot(W);

Set two starting point \( \Ss = \{x_0^1,x_0^2\} \) (you can use other points).

Exercice 5: (check the solution) Extract the set of points that are along the boundary of the Voronoi region. This corresponds for instance to the points
of the region \( \enscond{x}{Q(x)=1} \) that have one neighbor inside the region \( \enscond{x}{Q(x)=2} \). Compute the geodesic
distance \(D(x)\) at these points, and choose two points \(a\) and \(b\) on this boundary that have small values of \(D\).

exo5;

Exercice 6: (check the solution) Extract the geodesics joining \(a\) and \(b\) to the two starting points (this makes 4 geodesic curves). Use them to perform
segmentation.

exo6;

Vessel Segmentation and Centerline Extraction

One can extract a network of geodesic curve starting from a central point to detect vessels in medical images.

Load an image. This image is extracted from the DRIVE database of retinal vessels.

n = 256;
name = 'vessels';
f = rescale(load_image(name, n));

Display it.

clf;
imageplot(f);

We clean the image by substracting the smoothly varying background \[ f_1 = f - G_\si \star f, \] where \(G_\si\) is a Gaussian
kernel of variance \(\si^2\). Computing \(f_1\) corresponds to a high pass filtering.