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()
In [31]:
%% Insert your code here.

Anisotropic Image Approximation

One can combine the farthest point sampling strategy with a geodesic Delaunay triangulation to perform image approximation.

Load an image.

In [32]:
n = 256;
name = 'peppers-bw';
M = rescale(load_image(name, n));

Exercise 3

Compute a metric |H| adapted to the approximation of this image.

In [33]:
exo3()
In [34]:
%% Insert your code here.

Exercise 4

Perform farthest point sampling.

In [35]:
exo4()
In [36]:
%% Insert your code here.

Exercise 5

Compute the geodesic Delaunay triangulation of this set of point.

In [37]:
exo5()
In [38]:
%% Insert your code here.

Exercise 6

Perform image approximation using linear splines.

In [39]:
exo6()
In [40]:
%% Insert your code here.