Geodesic Triangulation for Image Compression¶

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 triangulations to perform image compression.

In [2]:

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

Image Approximation with Triangulation¶

It is possible to approximate an image over a triangulation using piecewise linear splines.

Load an image.

In [3]:

name = 'cameraman';
n = 256;
M = rescale( load_image(name, n) );

Number of points used to compute the approximation. The more points, the smallest the error.

In [4]:

m = 400;

Seeds random points, include the corners into these points.

In [5]:

vertex = floor( rand(2,m-4)*(n-1) ) +1;
vertex(:,end+1:end+4) = [[1;1] [1;n] [n;n] [n;1]];

Compute a Delaunay triangulation.

In [6]:

faces = compute_delaunay(vertex);

A first way to perform an approximation with |m| triangles is to interpolate the image at the sampling points |vertex|.

In [7]:

vinterp = interp2(M, vertex(2,:), vertex(1,:));

Each |vinterp(i)| is the value of the approximating image at the |vertex(:,i)|. We compute a spline interpolation.

In [8]:

Minterp = compute_triangulation_interpolation(faces,vertex,vinterp, n);

Display the approximation.

In [9]:

clf;
subplot(1,2,1);
plot_triangulation(vertex,faces, M);
title('Triangulation');
subplot(1,2,2);
imageplot(clamp(Minterp), ['Interpolation, SNR=' num2str(snr(Minterp,M)) 'dB']);

Another, better way to compute the approximation is to compute coefficients |vapprox| that performs the best L2 approximation with linear spline.

In [10]:

vapprox = compute_orthoproj_triangulation(vertex, faces, M);
Mapprox = compute_triangulation_interpolation(faces,vertex,vapprox, n);

Compare interpolation and approximation.

In [11]:

clf;
imageplot(clamp(Minterp), ['Interpolation, SNR=' num2str(snr(Minterp,M),3) 'dB'], 1,2,1);
imageplot(clamp(Mapprox), ['Approximation, SNR=' num2str(snr(Mapprox,M),3) 'dB'], 1,2,2);

Isotropic Metrics for Image Approximation¶

It is possible to compute optimized sampling location |vertex| by using the farthest point sampling algorithm with a well chosen metric |W| so that more points are put in areas of strong gradient.

The metric will be of the form

|W(x) = (norm(grad_x(M))+epsilon|)^alpha|

where |epsilon| and |alpha| control the density variation strength.

Parameters for the metric.

In [12]:

alpha = .7;
epsilon = 1e-2;

Exercise 1

In [13]:

exo1()

In [14]:

%% Insert your code here.

Exercise 2

Perform farthest points sampling to compute sampling location |vertex| and the corresponding geodesic Delaunay triangulation |faces|. ompute the Delaunay triangulation

In [15]:

exo2()

In [16]:

%% Insert your code here.

Perform approximation over the triangulation.

In [17]:

vgeod = compute_orthoproj_triangulation(vertex, faces, M);
Mgeod = compute_triangulation_interpolation(faces,vertex,vgeod, n);

Compare interpolation and approximation.

In [18]:

clf;
subplot(1,2,1);
plot_triangulation(vertex,faces, M);
subplot(1,2,2);
imageplot(clamp(Mgeod), ['SNR=' num2str(snr(Mgeod,M),3) 'dB']);