Anisotropic Fast Marching¶

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 explores the use of geodesic distances for anisotropic metric.

In [8]:

warning off
addpath('toolbox_signal')
addpath('toolbox_general')
addpath('toolbox_graph')
addpath('solutions/fastmarching_3_anisotropy')
warning on

Structure Tensor Field¶

An anisotropic metric is given through a tensfor field |T| which is an |(n,n,2,2)| array, where |T(i,j,:,:)| is a positive definite symmetric matrix that define the metric at pixel |(i,j)|.

Here we extract the tensor field whose main eigenvector field is alligned with the direction of the texture. This can be achieved using the structure tensor field, which remove the sign ambiguity by tensorizing the gradient, and remove noise by filtering.

Load an image that contains an oscillating texture.

In [9]:

name = 'fingerprint';
n = 150;
M = rescale(load_image(name,n));

Compute its gradient.

In [10]:

options.order = 2;
G = grad(M,options);

Compute a rank-1 tensor with main eigenvector aligned in the direction orthogonal to the gradient.

In [11]:

T = zeros(n,n,2,2);
T(:,:,1,1) = G(:,:,2).^2;
T(:,:,2,2) = G(:,:,1).^2;
T(:,:,1,2) = -G(:,:,1).*G(:,:,2);
T(:,:,2,1) = -G(:,:,1).*G(:,:,2);

Smooth the field (blur each entry).

In [12]:

sigma = 12;
T = perform_blurring(T,sigma);

Compute the eigenvector and eigenvalues of the tensor field.

In [13]:

[e1,e2,l1,l2] = perform_tensor_decomp(T);

Display the main orientation field.

In [14]:

clf;
plot_vf(e1(1:10:n,1:10:n,:),M);
colormap(gray(256));

Anisotropic Fast Marching¶

One can compute geodesic distance and geodesics using an anisotropic Fast Marching propagation.

Anisotropy of the tensor field.

In [15]:

anisotropy = .1;

Build the Riemannian metric using the structure tensor direction.

In [16]:

H = perform_tensor_recomp(e1,e2, ones(n),ones(n)*1/anisotropy );

Starting point.

In [17]:

pstart = [n n]/4;

Perform the propagation.

In [18]:

hx = 1/n; hy = 1/n;
[D, dUx, dUy, Vor, L] = fm2dAniso([hx;hy], H, pstart);

Display the result.

In [19]:

clf;
subplot(1,2,1);
imageplot(M, 'Image');
subplot(1,2,2);
hold on;
imageplot(convert_distance_color(D), 'Geodesic distance');
hh = plot(pstart(2),pstart(1), 'r.');
set(hh, 'MarkerSize',15);
axis('ij');
colormap(gray(256));

Exercise 1

Compute the geodesic distance for several anisotropy, and for several starting points.

In [20]:

exo1()

In [21]:

%% Insert your code here.

Farthest Point Sampling¶

It is possible to perform an anisotropic sampling of the image using a Farthest point sampling strategy.

We use a highly anisotropic metric.

In [22]:

anisotropy = .02;
H = perform_tensor_recomp(e1,e2, ones(n),ones(n)*1/anisotropy );

Choose the initial point.

In [23]:

vertex = [1;1];

Compute the geodesic distance.

In [24]:

[D, dUx, dUy, Vor, L] = fm2dAniso([hx;hy], H, vertex);

Choose the second point as the farthest point.

In [25]:

[tmp,i] = max(D(:));
[x,y] = ind2sub([n n],i); 
vertex(:,end+1) = [x;y];

Display distance and points.

In [26]:

clf;
subplot(1,2,1);
hold on;
imageplot(M, 'Image'); axis ij;
plot(vertex(2,1), vertex(1,1), 'r.');
plot(vertex(2,2), vertex(1,2), 'b.');
subplot(1,2,2);
hold on;
imageplot( convert_distance_color(D), 'Distance'); axis ij;
plot(vertex(2,1), vertex(1,1), 'r.');
plot(vertex(2,2), vertex(1,2), 'b.');
colormap gray(256);

Update the value of the distance map with a partial propagation from the last added point.

In [27]:

[D1, dUx, dUy, Vor, L] = fm2dAniso([hx;hy], H, vertex);

Display old/new.

In [28]:

clf;
imageplot( D, 'Old distance', 1,2,1 );
imageplot( D1, 'New distance', 1,2,2 );
colormap jet(256);

Update.

In [29]:

D = D1;

Exercise 2

Perform farthest point sampling.

In [30]:

exo2()