Geodesic Medial Axsis

Important: Please read the installation page for details about how to install the toolboxes. $\newcommand{\dotp}[2]{\langle #1, #2 \rangle}$ $\newcommand{\enscond}[2]{\lbrace #1, #2 \rbrace}$ $\newcommand{\pd}[2]{ \frac{ \partial #1}{\partial #2} }$ $\newcommand{\umin}[1]{\underset{#1}{\min}\;}$ $\newcommand{\umax}[1]{\underset{#1}{\max}\;}$ $\newcommand{\umin}[1]{\underset{#1}{\min}\;}$ $\newcommand{\uargmin}[1]{\underset{#1}{argmin}\;}$ $\newcommand{\norm}[1]{\|#1\|}$ $\newcommand{\abs}[1]{\left|#1\right|}$ $\newcommand{\choice}[1]{ \left\{ \begin{array}{l} #1 \end{array} \right. }$ $\newcommand{\pa}[1]{\left(#1\right)}$ $\newcommand{\diag}[1]{{diag}\left( #1 \right)}$ $\newcommand{\qandq}{\quad\text{and}\quad}$ $\newcommand{\qwhereq}{\quad\text{where}\quad}$ $\newcommand{\qifq}{ \quad \text{if} \quad }$ $\newcommand{\qarrq}{ \quad \Longrightarrow \quad }$ $\newcommand{\ZZ}{\mathbb{Z}}$ $\newcommand{\CC}{\mathbb{C}}$ $\newcommand{\RR}{\mathbb{R}}$ $\newcommand{\EE}{\mathbb{E}}$ $\newcommand{\Zz}{\mathcal{Z}}$ $\newcommand{\Ww}{\mathcal{W}}$ $\newcommand{\Vv}{\mathcal{V}}$ $\newcommand{\Nn}{\mathcal{N}}$ $\newcommand{\NN}{\mathcal{N}}$ $\newcommand{\Hh}{\mathcal{H}}$ $\newcommand{\Bb}{\mathcal{B}}$ $\newcommand{\Ee}{\mathcal{E}}$ $\newcommand{\Cc}{\mathcal{C}}$ $\newcommand{\Gg}{\mathcal{G}}$ $\newcommand{\Ss}{\mathcal{S}}$ $\newcommand{\Pp}{\mathcal{P}}$ $\newcommand{\Ff}{\mathcal{F}}$ $\newcommand{\Xx}{\mathcal{X}}$ $\newcommand{\Mm}{\mathcal{M}}$ $\newcommand{\Ii}{\mathcal{I}}$ $\newcommand{\Dd}{\mathcal{D}}$ $\newcommand{\Ll}{\mathcal{L}}$ $\newcommand{\Tt}{\mathcal{T}}$ $\newcommand{\si}{\sigma}$ $\newcommand{\al}{\alpha}$ $\newcommand{\la}{\lambda}$ $\newcommand{\ga}{\gamma}$ $\newcommand{\Ga}{\Gamma}$ $\newcommand{\La}{\Lambda}$ $\newcommand{\si}{\sigma}$ $\newcommand{\Si}{\Sigma}$ $\newcommand{\be}{\beta}$ $\newcommand{\de}{\delta}$ $\newcommand{\De}{\Delta}$ $\newcommand{\phi}{\varphi}$ $\newcommand{\th}{\theta}$ $\newcommand{\om}{\omega}$ $\newcommand{\Om}{\Omega}$

This tour studies the computation of the medial axis using the Fast Marching.

In [2]:
addpath('toolbox_signal')
addpath('toolbox_general')
addpath('toolbox_graph')
addpath('solutions/shapes_6_medialaxis')

Voronoi Diagram

The Voronoi diagram is the segmentation of the image given by the region of influence of the set of starting points.

Load a distance map.

In [3]:
n = 200;
W = load_image('mountain', n);
W = rescale(W,.25,1);

Select seed points.

In [4]:
pstart = [[20;20] [120;100] [180;30] [60;160]];
nbound = size(pstart,2);

Display the map and the points.

In [5]:
ms = 20;
clf; hold on;
imageplot(W);
h = plot(pstart(2,:), pstart(1,:), '.r'); set(h, 'MarkerSize', ms);

Compute the geodesic distant to the whole set of points.

In [6]:
[D,S,Q] = perform_fast_marching(W, pstart);

Display the geodesic distance.

In [7]:
clf; hold on;
imageplot(convert_distance_color(D, W));
h = plot(pstart(2,:), pstart(1,:), '.r'); set(h, 'MarkerSize', ms);

Display the Voronoi Segmentation.

In [8]:
clf; hold on;
imageplot(Q);
h = plot(pstart(2,:), pstart(1,:), '.r'); set(h, 'MarkerSize', ms);
colormap jet(256);

Medial Axis from the Voronoi Map

The medial axis is difficult to extract from the singularity of the distance map. It is much more robust to extract it from the discontinuities in the Voronoi index map |Q|.

Compute the derivative, the gradient.

In [9]:
G = grad(Q);

Take it modulo |nbound|.

In [10]:
G(G<-nbound/2) = G(G<-nbound/2) + nbound;
G(G>nbound/2) = G(G>nbound/2) - nbound;

Compute the norm of the gadient.

In [11]:
G = sqrt(sum(G.^2,3));

Compute the medial axis by thresholding the gradient magnitude.

In [12]:
B = 1 - (G>.1);

Display.

In [13]:
clf; hold on;
imageplot(B);
h = plot(pstart(2,:), pstart(1,:), '.r'); set(h, 'MarkerSize', ms);

Skeleton of a Shape

The sekeleton, also called Medial Axis, is the set of points where the geodesic distance is singular.

A binary shape is represented as a binary image.

In [14]:
n = 200;
name = 'chicken';
M = load_image(name,n);
M = perform_blurring(M,5);
M = double( rescale( M )>.5 );
if M(1)==1 
    M = 1-M;
end

Compute its boundary, that is going to be the set of starting points.

In [15]:
pstart = compute_shape_boundary(M);
nbound = size(pstart,2);

Display the metric.

In [16]:
lw = 2;
clf; hold on;
imageplot(-M);
h = plot(pstart(2,:), pstart(1,:), 'r'); set(h, 'LineWidth', lw); axis ij;

Parameters for the Fast Marching: constant speed |W|, but retricted using |L| to the inside of the shape.

In [17]:
W = ones(n);
L = zeros(n)-Inf; L(M==1) = +Inf;

Compute the fast marching, from the boundary points.

In [18]:
options.constraint_map = L;
[D,S,Q] = perform_fast_marching(W, pstart, options);
D(M==0) = Inf;

Display the distance function to the boundary.

In [19]:
clf;
hold on;
display_shape_function(D);
h = plot(pstart(2,:), pstart(1,:), 'r'); set(h, 'LineWidth', lw); axis ij;

Display the index of the closest boundary point.

In [20]:
clf;
hold on;
display_shape_function(Q);
h = plot(pstart(2,:), pstart(1,:), 'r'); set(h, 'LineWidth', lw); axis ij;

Exercise 1

Compute the norm of the gradient |G| modulo |nbound|. Be careful to remove the boundary of the shape from this indicator. Display the thresholded gradient map. radient ompute the norm of the gadient. emove the boundary to the skeletton.

In [21]:
exo1()
In [22]:
%% Insert your code here.

Exercise 2

Display the Skeleton obtained for different threshold values.

In [23]:
exo2()
In [24]:
%% Insert your code here.

Skeletton