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.

Managing level set functions

In the level set formalism, the evolution of some curve \( (\ga(t))_{t=0}^1 \) is computed by evolving the zero level of a
function \(\phi : \RR^2 \rightarrow \RR \) \[ \enscond{\ga(s)}{ s \in [0,1] } = \enscond{x \in \RR^2}{\phi(x)=0}. \] This
corresponds to replacing the parameteric representation \(\ga\) of the curve by an implicit representation. This requires
an additional dimension (and hence more storage) but ease the handling of topological change of the curve during the evolution.

Discretazion size \(n \times n\) of the domain \([0,1]^2\).

n = 200;
[Y,X] = meshgrid(1:n,1:n);

One can create a circular shape by using the signed distance function to a circle \[ \phi_1(x) = \sqrt{ (x_1-c_1)^2 + (x_2-c_2)^2
} - r \] where \(r>0\) is the radius and \(c \in \RR^2\) the center.

Radius \(r\).

r = n/3;

Center \(c\).

c = [r r] + 10;

Distance function \(\phi_1\).

phi1 = sqrt( (X-c(1)).^2 + (Y-c(2)).^2 ) - r;

Exercice 1: (check the solution) Load a square shape \(\phi_2\) at a different position for the center.

Exercice 2: (check the solution) Compute the intersection and the union of the two shapes. Store the union in \(\phi_0\) (variable phi0), that we will use in the remaining part of the tour.

exo2;

Mean Curvature Motion.

The mean curvature motion corresponds to the minimizing flow of the length of the curve \[ \int_0^1 \norm{\ga'(s)} d s. \]

It is implemeted in a level set formalism by a familly \(\phi_t\) of level set function parameterized by an artificial time
\(t \geq 0\), that satisfies the following PDE \[ \pd{\phi_t}{t} = -G(\phi_t) \qwhereq G(\phi) = -\norm{\nabla \phi} \text{div}
\pa{ \frac{\nabla \phi}{\norm{\nabla \phi}} } \] and where \(\nabla \phi_t(x) \in \RR^2\) is the spacial gradient.

This flow is computed using a gradient descent \(\phi^{(0)} = \phi_0\) and This is implemented using a gradient descent scheme.
\[ \phi^{(\ell+1)} = \phi^{(\ell)} - \tau G(\phi^{(\ell)}), \] where \(\tau>0\) is small enough time step.

Levelset Re-distancing

During PDE resolution, a level set function \(\phi\) might become ill-conditionned, so that the zero crossing is not sharp
enough. The quality of the level set function is restored by computing the signed distance function to the zero level set.

This corresponds to first extracting the zero level set \[ \Cc = \enscond{x \in \RR^2 }{\phi(x)=0}, \] and then solving the
following eikonal equation PDE on \(\tilde \phi\) (in viscosity sense) \[ \norm{\nabla \tilde \phi(x)} = 1 \qandq \forall
y \in \Cc, \tilde\phi(y)=0. \] The one can replace \(\phi\) by \(\text{sign}(\phi(x))\tilde \phi(x)\) which is the signed
distance function to \(\Cc\).

We set \(\phi=\phi_0^3\) so that they are both valid level set function of the same curve, but \(\phi\) is not the signed
distance function..

Edge-based Segmentation with Geodesic Active Contour

Geodesic active contours compute loval minimum of a weighted geodesic distance that attract the curve toward the features
of the background image.

Note: these active contours should not be confounded with the geodesic shortest paths, that are globally minimizing geodesics between
two points. Here the active contour is a close curve progressively decreasing a weighted geodesic length that is only a local
minimum (the global minimum would be a single point).

Size of the image.

n = 200;

First we load an image \(f_0 \in \RR^{n \times n}\) to segment.

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

Given a background image \(f_0\) to segment, one needs to compute an edge-stopping function \(W\). It should be small in area
of high gradient, and high in area of large gradient.

Exercice 5: (check the solution) Compute and store in G the gradient \(G(\phi)\) (right hand side of the PDE) using the current value of the distance function \(\phi\).

exo5;

Do the descent step.

phi = phi - tau*G;

Once in a while (e.g. every 30 iterations), perform re-distancing of \(\phi\).

phi = perform_redistancing(phi);

Exercice 6: (check the solution) Implement the geodesic active contours gradient descent. Do not forget to do the re-distancing.

exo6;

Region-based Segmentation with Chan-Vese

Chan-Vese active contours corresponds to a region-based energy that looks for a piecewise constant approximation of the image.

The energy to be minimized is \[ \umin{\phi} L(\phi) + \la \int_{\phi(x)>0} \abs{f_0(x)-c_1}^2 d x + \la
\int_{\phi(x)<0} \abs{f_0(x)-c_2}^2 d x \] where \(L\) is the length of the zero level set of \(\phi\). Note that here \((c_1,c_2)
\in \RR^2\) are assumed to be known.

Exercice 7: (check the solution) Compute an initial level set function \(\phi_0\), stored in phi0, for instance many small circles.