Geodesic Farthest Point Sampling¶

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 geodesic computation to perform image sampling.

In [2]:

addpath('toolbox_signal')
addpath('toolbox_general')
addpath('toolbox_graph')
addpath('solutions/fastmarching_5_sampling_2d')

Voronoi Segmentation¶

A geodesic Voronoi segementation is obtained by computing a geodesic distance from multiple starting points.

Compute an image with bumps.

In [3]:

n = 256; % size of the image
sigma = n/8; % width of the bumps
[Y,X] = meshgrid(1:n,1:n);
x = n/4; y = n/4;
M = exp( -( (X-x).^2 + (Y-y).^2 )/(2*sigma^2) );
x = 3*n/4; y = 3*n/4;
M = M + exp( -( (X-x).^2 + (Y-y).^2 )/(2*sigma^2) );

Compute a metric by rescaling |M|.

In [4]:

W = rescale(M,1e-2,1);

Create random starting points.

In [5]:

m = 20; 
pstart = floor( rand(2,m)*(n-1) ) +1;

Perform the propagation using the Fast Marching.

In [6]:

[D,Z,Q] = perform_fast_marching(1./W, pstart);

Display the distance map.

In [7]:

clf;
subplot(1,2,1);
hold on;
imageplot(perform_hist_eq(D,'linear'));  title('Geodesic distance'); 
plot(pstart(2,:), pstart(1,:), 'r.');
subplot(1,2,2);
hold on;
imageplot(Q); title('Voronoi'); 
plot(pstart(2,:), pstart(1,:), 'r.');
colormap jet(256);

Geodesic Delaunay Triangulation¶

A geodesic Delaunay triangulation is obtained by linking starting points whose Voronoi cells touch. This is the dual of the original Voronoi segmentation.

Exercise 1

Using |Q|, compute the faces |faces| of the Delaunay triangulation. To that end, extract each quad of values |Q(i,j),Q(i+1,j),Q(i+1,j+1),Q(i,j+1)|, and add a new face when three of these four values are different (this corresponds to a Voronoi point). Display the obtained triangulation. ompute quad of values

In [8]:

exo1()

In [9]:

%% Insert your code here.

Display the obtained triangulation.

In [10]:

clf;
subplot(1,2,1);
hold on;
imageplot(Q, 'Voronoi'); axis ij;
plot(pstart(2,:), pstart(1,:), 'r.');
subplot(1,2,2);
plot_triangulation(pstart,faces, M);
colormap jet(256);

Compare this Geodesic Delaunay triangulation with the Euclidean triangulation.

In [11]:

faces_euc = compute_delaunay(pstart);

Display.

In [12]:

clf;
subplot(1,2,1);
plot_triangulation(pstart,faces, M); title('Geodesic');
subplot(1,2,2);
plot_triangulation(pstart,faces_euc, M); title('Euclidean');
colormap jet(256);

Farthest Point Sampling¶

To sample point uniformly according to the geodesic distance, one can use an iterative farthest point sampling scheme.

Construct a metric that is large in area where we want to sample more points. The front should move slowly in high density sampling region.

In [13]:

W = rescale(M,3*1e-1,1);

Choose the initial point.

In [14]:

vertex = [1;1];

Compute the geodesic distance.

In [15]:

[D,Z,Q] = perform_fast_marching(1./W, vertex);

Choose the second point as the farthest point.

In [16]:

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

Display distance and points.

In [17]:

clf;
subplot(1,2,1);
hold on;
imageplot(W, 'Metric'); 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( perform_hist_eq(D, 'linear'), 'Distance'); axis ij;
plot(vertex(2,1), vertex(1,1), 'r.');
plot(vertex(2,2), vertex(1,2), 'b.');
colormap jet(256);

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

In [18]:

options.constraint_map = D;
[D1,Z,Q] = perform_fast_marching(1./W, vertex(:,end), options);

display old/new

In [19]:

clf;
imageplot( convert_distance_color(D,M,1), 'Update distance', 1,3,1 );
imageplot( convert_distance_color(D1,M,1), 'Update distance', 1,3,2 );
imageplot( convert_distance_color(min(D,D1),M,1), 'New distance', 1,3,3 );

update

In [20]:

D = min(D,D1);

Exercise 2

Iterate the sampling process to add more and more points.

In [21]:

exo2()

In [22]:

%% Insert your code here.

Exercise 3

Display the geodesic Delaunay triangulation corresponding to the sampling

In [23]:

exo3()

In [24]:

%% Insert your code here.

Lloyd Relaxation¶

The farthest point sampling strategy is greedy and does not move the positions of the points once they are seeded.

To enhance the sampling, it is possible to relocate iteratively the points at the center of the Voronoi cells. This corresponds to the Lloyd algorithm, first developped for vector quantization.

First with a constant metric.

In [25]:

n = 512;
W = ones(n);

Seed random points.

In [26]:

p = 40;
vertex = floor(rand(2,p)*n-1)+1;

Compute Voronoi partition.

In [27]:

[D,Z,Q] = perform_fast_marching(1./W, vertex);

Display Vornoi cells.

In [28]:

clf; hold on;
imageplot(Q');
h = plot(vertex(1,:), vertex(2,:), 'k.');
set(h, 'MarkerSize', 15);
colormap(jet(256));

Re-center each point at the barycenter of its cell.

In [29]:

for i=1:p
    [x,y] = ind2sub(size(W), find(Q==i));
    vertex(:,i) = [mean(x);mean(y)];
end

Display updated partitions.

In [30]:

[D,Z,Q] = perform_fast_marching(1./W, vertex);
clf; hold on;
imageplot(Q');
h = plot(vertex(1,:), vertex(2,:), 'k.');
set(h, 'MarkerSize', 15);
colormap(jet(256));

Exercise 4

Perform the Lloyd iterative algorithm. etric eed random points.

In [31]:

exo4()

In [32]:

%% Insert your code here.

Now we define a non-constent metric.

In [33]:

n = 256;
x = linspace(-1,1,n);
[Y,X] = meshgrid(x,x);
sigma = .3;
W = exp( -(X.^2+Y.^2)/(2*sigma^2) );
W = rescale(W,.02,1);

Display.

In [34]:

clf;
imageplot(W);

Exercise 5

Perform the Lloyd iterative algorithm. eed random points.

In [35]:

exo5()

In [36]:

%% Insert your code here.