\[ \newcommand{\NN}{\mathbb{N}} \newcommand{\CC}{\mathbb{C}} \newcommand{\GG}{\mathbb{G}} \newcommand{\LL}{\mathbb{L}} \newcommand{\PP}{\mathbb{P}} \newcommand{\QQ}{\mathbb{Q}} \newcommand{\RR}{\mathbb{R}} \newcommand{\VV}{\mathbb{V}} \newcommand{\ZZ}{\mathbb{Z}} \newcommand{\FF}{\mathbb{F}} \newcommand{\KK}{\mathbb{K}} \newcommand{\UU}{\mathbb{U}} \newcommand{\EE}{\mathbb{E}} \newcommand{\Aa}{\mathcal{A}} \newcommand{\Bb}{\mathcal{B}} \newcommand{\Cc}{\mathcal{C}} \newcommand{\Dd}{\mathcal{D}} \newcommand{\Ee}{\mathcal{E}} \newcommand{\Ff}{\mathcal{F}} \newcommand{\Gg}{\mathcal{G}} \newcommand{\Hh}{\mathcal{H}} \newcommand{\Ii}{\mathcal{I}} \newcommand{\Jj}{\mathcal{J}} \newcommand{\Kk}{\mathcal{K}} \newcommand{\Ll}{\mathcal{L}} \newcommand{\Mm}{\mathcal{M}} \newcommand{\Nn}{\mathcal{N}} \newcommand{\Oo}{\mathcal{O}} \newcommand{\Pp}{\mathcal{P}} \newcommand{\Qq}{\mathcal{Q}} \newcommand{\Rr}{\mathcal{R}} \newcommand{\Ss}{\mathcal{S}} \newcommand{\Tt}{\mathcal{T}} \newcommand{\Uu}{\mathcal{U}} \newcommand{\Vv}{\mathcal{V}} \newcommand{\Ww}{\mathcal{W}} \newcommand{\Xx}{\mathcal{X}} \newcommand{\Yy}{\mathcal{Y}} \newcommand{\Zz}{\mathcal{Z}} \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} \renewcommand{\phi}{\varphi} \renewcommand{\th}{\theta} \newcommand{\om}{\omega} \newcommand{\Om}{\Omega} \renewcommand{\epsilon}{\varepsilon} \newcommand{\Calpha}{\mathrm{C}^\al} \newcommand{\Cbeta}{\mathrm{C}^\be} \newcommand{\Cal}{\text{C}^\al} \newcommand{\Cdeux}{\text{C}^{2}} \newcommand{\Cun}{\text{C}^{1}} \newcommand{\Calt}[1]{\text{C}^{#1}} \newcommand{\lun}{\ell^1} \newcommand{\ldeux}{\ell^2} \newcommand{\linf}{\ell^\infty} \newcommand{\ldeuxj}{{\ldeux_j}} \newcommand{\Lun}{\text{\upshape L}^1} \newcommand{\Ldeux}{\text{\upshape L}^2} \newcommand{\Lp}{\text{\upshape L}^p} \newcommand{\Lq}{\text{\upshape L}^q} \newcommand{\Linf}{\text{\upshape L}^\infty} \newcommand{\lzero}{\ell^0} \newcommand{\lp}{\ell^p} \renewcommand{\d}{\ins{d}} \newcommand{\Grad}{\text{Grad}} \newcommand{\grad}{\text{grad}} \renewcommand{\div}{\text{div}} \newcommand{\diag}{\text{diag}} \newcommand{\pd}[2]{ \frac{ \partial #1}{\partial #2} } \newcommand{\pdd}[2]{ \frac{ \partial^2 #1}{\partial #2^2} } \newcommand{\dotp}[2]{\langle #1,\,#2\rangle} \newcommand{\norm}[1]{|\!| #1 |\!|} \newcommand{\normi}[1]{\norm{#1}_{\infty}} \newcommand{\normu}[1]{\norm{#1}_{1}} \newcommand{\normz}[1]{\norm{#1}_{0}} \newcommand{\abs}[1]{\vert #1 \vert} \newcommand{\argmin}{\text{argmin}} \newcommand{\argmax}{\text{argmax}} \newcommand{\uargmin}[1]{\underset{#1}{\argmin}\;} \newcommand{\uargmax}[1]{\underset{#1}{\argmax}\;} \newcommand{\umin}[1]{\underset{#1}{\min}\;} \newcommand{\umax}[1]{\underset{#1}{\max}\;} \newcommand{\pa}[1]{\left( #1 \right)} \newcommand{\choice}[1]{ \left\{ \begin{array}{l} #1 \end{array} \right. } \newcommand{\enscond}[2]{ \left\{ #1 \;:\; #2 \right\} } \newcommand{\qandq}{ \quad \text{and} \quad } \newcommand{\qqandqq}{ \qquad \text{and} \qquad } \newcommand{\qifq}{ \quad \text{if} \quad } \newcommand{\qqifqq}{ \qquad \text{if} \qquad } \newcommand{\qwhereq}{ \quad \text{where} \quad } \newcommand{\qqwhereqq}{ \qquad \text{where} \qquad } \newcommand{\qwithq}{ \quad \text{with} \quad } \newcommand{\qqwithqq}{ \qquad \text{with} \qquad } \newcommand{\qforq}{ \quad \text{for} \quad } \newcommand{\qqforqq}{ \qquad \text{for} \qquad } \newcommand{\qqsinceqq}{ \qquad \text{since} \qquad } \newcommand{\qsinceq}{ \quad \text{since} \quad } \newcommand{\qarrq}{\quad\Longrightarrow\quad} \newcommand{\qqarrqq}{\quad\Longrightarrow\quad} \newcommand{\qiffq}{\quad\Longleftrightarrow\quad} \newcommand{\qqiffqq}{\qquad\Longleftrightarrow\qquad} \newcommand{\qsubjq}{ \quad \text{subject to} \quad } \newcommand{\qqsubjqq}{ \qquad \text{subject to} \qquad } \]

Shape Retrieval with Geodesic Descriptors

Shape Retrieval with Geodesic Descriptors

This numerical tour explores the use of geodesic distances within shapes to perform shape retrieval.

Contents

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.

Installing toolboxes and setting up the path.

You need to download the following files: signal toolbox, general toolbox and graph toolbox.

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.

Execute this line only if you are using Matlab.

getd = @(p)path(p,path); % scilab users must *not* execute this

Then you can add the toolboxes to the path.

getd('toolbox_signal/');
getd('toolbox_general/');
getd('toolbox_graph/');

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.

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.

clf;
imageplot(-M);

Compute its boundary

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.

W = ones(n);
L = zeros(n)-Inf; L(M==1) = +Inf;

Initial point for the geodesic computation.

start_points = [95; 20];

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

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

Display Euclidean distance.

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.

options.constraint_map = L;
D = perform_fast_marching(W, start_points, options);

Display geodesic distance.

clf;
options.nbr_levelsets = 60;
display_shape_function(D, options);

Exercice 1: (check the solution) Using options.nb_iter_max display the progression of the Fast Marching.

exo1;

Compute a geodesic curve.

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

Display the path.

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;

Exercice 2: (check the solution) Compute curves joining the start point to several points along the boundary.

exo2;

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.

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

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

exo3;

normalize distances.

E = E/mean(E(:));

Display some locations

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.

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.

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.

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.

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.

Exercice 4: (check the solution) Load a library of shapes. Compute the different histograms for these shapes.

exo4;

Exercice 5: (check the solution) Perform the retrieval by comparing the histogram. Test diffetent metrics for the retrieval.

exo5;