Shape Retrieval with Geodesic Descriptors

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 numerical tour explores the use of geodesic distances within shapes to perform shape retrieval.

This tour is mostly inspired from the following work:

Matching 2D and 3D Articulated Shapes using Eccentricity, A. Ion, N. M. Artner, G. Peyre, W. G. Kropatsch and L. Cohen, Preprint Hal-00365019, January 2009.

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

Geodesic Distances Within a Binary Shape

By restricting shortest path to lie within a shape, one create a geodesic metric that is different from the Euclidean one if the shape is not convex.

A binary shape is represented as a binary image.

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

Display the shape.

In [4]:
clf;
imageplot(-M);

Compute its boundary

In [5]:
bound = compute_shape_boundary(M);
nbound = size(bound,2);

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

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

Initial point for the geodesic computation.

In [7]:
start_points = [95; 20];

Compute the geodesic distance without constraint using Fast Marching. It is simply the Euclidean distance.

In [8]:
options.constraint_map = [];
D0 = perform_fast_marching(W, start_points, options);
D0(M==0) = Inf;

Display Euclidean distance.

In [9]:
clf;
options.display_levelsets = 1;
options.pstart = start_points;
options.nbr_levelsets = 30;
display_shape_function(D0, options);

Compute the geodesic distance with constraints using Fast Marching.

In [10]:
options.constraint_map = L;
D = perform_fast_marching(W, start_points, options);

Display geodesic distance.

In [11]:
clf;
options.nbr_levelsets = 60;
display_shape_function(D, options);

Exercise 1

Using |options.nb_iter_max| display the progression of the Fast Marching.

In [12]:
exo1()
In [13]:
%% Insert your code here.

Compute a geodesic curve.

In [14]:
end_points = [27;112];
p = compute_geodesic(D,end_points);

Display the path.

In [15]:
ms = 30; lw = 3;
clf; hold on;
imageplot(1-M);
h = plot(end_points(2),end_points(1), '.b'); set(h, 'MarkerSize', ms);
h = plot(start_points(2),start_points(1), '.r'); set(h, 'MarkerSize', ms);
h = plot( p(2,:), p(1,:), 'g' ); set(h, 'LineWidth', lw);
axis ij;

Exercise 2

Compute curves joining the start point to several points along the boundary.

In [16]:
exo2()
In [17]:
%% Insert your code here.

Local Geodesic Descriptors

In order to build shape signatures, we compute geodesic distance to all the points on the boundary. We then retrieve some caracteristics from these geodesic distance map.

Select a uniform set of points on the boundary.

In [18]:
nb_samples = 600;
sel = round(linspace(1,nbound+1,nb_samples+1)); sel(end) = [];
samples = bound(:,sel);

Exercise 3

Build a collection |E| of distance maps, so that |E(:,:,i)| is the geodesic distance to |samples(:,i)|.

In [19]:
exo3()
In [20]:
%% Insert your code here.

normalize distances.

In [21]:
E = E/mean(E(:));

Display some locations

In [22]:
points = [[80;20] [95;112] [156;42]];
col = {'r', 'g', 'b', 'k'};
clf; hold on;
imageplot(-M);
for i=1:3
    h = plot(points(2,i), points(1,i), [col{i} '.']);
    set(h, 'MarkerSize', 40);
end
axis('ij');

Display three different features at some locations.

In [23]:
clf;
col = {'r', 'g', 'b', 'k'};
for i=1:3
    subplot(3,1,i);    
    d = E(points(1,i),points(2,i), :); 
    u = hist(d(:), 15); axis tight;
    h = bar(u, col{i}); axis('tight');
    set(gca, 'XTickLabel', []);
end

Global Geodesic Descriptors

One can retain a single statistic from the local descriptors, such as the min, max, mean or median values. The histogram of these values are the global descriptors.

Compute several statistics.

In [24]:
clear A;
A{1} = max(E,[],3);
A{2} = min(E,[],3);
A{3} = mean(E,3);
A{4} = median(E,3);
titles = {'Max', 'Min', 'Mean', 'Median'};

Display as images.

In [25]:
nbr = [20 5 30 30];
options.pstart = [];
clf;
for i=1:4
    subplot(2,2,i);
    options.nbr_levelsets = nbr(i);
    display_shape_function(A{i}, options);
    title(titles{i});
end
colormap jet(256);

Display histograms of the statistics.

In [26]:
clf;
for i=1:4
    u = A{i}(M==1); u = u(u>0);
    subplot(4,1,i);
    hist(u, 40); axis('tight');
    title(titles{i});
end

Shape Retrieval using Geodesic Historams.

One can use the histograms of Eccentricity for shape retrieval.

Exercise 4

Load a library of shapes. Compute the different histograms for these shapes.

In [27]:
exo4()
In [28]:
%% Insert your code here.

Exercise 5

Perform the retrieval by comparing the histogram. Test diffetent metrics for the retrieval.

In [29]:
exo5()
In [30]:
%% Insert your code here.