# Image Approximation with Orthogonal Bases¶


This numerical tour uses several orthogonal bases to perform non-linear image approximation.

Important: Scilab user have to call the function |extend_stack_size(4)| before starting the tour to be able to handle large images.

In [2]:
addpath('toolbox_signal')


## Best $M$-terms Non-linear Approximation¶

This tours makes use of an orthogonal base $\Bb = \{ \psi_m \}_{m=0}^{N-1}$ of the space $\RR^N$ of the images with $N$ pixels.

The best $M$-term approximation of $f$ is obtained by a non-linear thresholding $$f_M = \sum_{ \abs{\dotp{f}{\psi_m}}>T } \dotp{f}{\psi_m} \psi_m,$$ where the value of $T>0$ should be carefully selected so that only $M$ coefficients are not thresholded, i.e. $$\abs{ \enscond{m}{ \abs{\dotp{f}{\psi_m}}>T } } = M.$$

The goal is to use an ortho-basis $\Bb$ so that the error $\norm{f-f_M}$ decays as fast as possible when $M$ increases, for a large class of images.

This tour studies several different orthogonal bases: Fourier, wavelets (which is at the heart of JPEG-2000), cosine, local cosine (which is at the heart of JPEG).

First we load an image of $N = n \times n$ pixels.

In [3]:
n = 512;
f = rescale( load_image('lena', n) );


Display it.

In [4]:
clf;
imageplot(f);


## Fourier Approximation¶

The discrete 2-D Fourier atoms are defined as: $$\psi_m(x) = \frac{1}{\sqrt{N}} e^{ \frac{2i\pi}{n} ( x_1 m_1 + x_2 m_2 ) },$$ where $0 \leq m_1,m_2 < n$ indexes the frequency.

The set of inner products $\{ \dotp{f}{\psi_m} \}_m$ is computed in $O(N \log(N))$ operations with the 2-D Fast Fourier Transform (FFT) algorithm (the Matlab function is |fft2|).

Compute the Fourier transform using the FFT algorithm. Note the normalization by $1/\sqrt{N}$ to ensure orthogonality (energy conservation) of the transform.

In [5]:
fF = fft2(f)/n;


Display its magnitude (in log scale). We use the function |fftshift| to put the low frequency in the center.

In [6]:
clf;
imageplot(log(1e-5+abs(fftshift(fF))));


An image is recovered from a set of coefficients $c_m$ using the inverse Fourier Transform (Matlab function |ifft2|) than implements the formula $$f_M = \sum_m c_m \psi_m.$$

Performs a thresholding.

In [7]:
T = .3;
c = fF .* (abs(fF)>T);


Inverse the Fourier transform.

In [8]:
fM = real(ifft2(c)*n);


Display the approximation.

In [9]:
imageplot(clamp(fM));


Exercise 1

Compute a best $M$-term approximation in the Fourier basis of $f$, for $M \in \{N/100, N/20\}$. Compute the approximation using a well chosen hard threshold value $T$.

In [10]:
exo1()

In [11]:
%% Insert your code here.


The best $M$-term approximation error is computed using the conservation of energy as $$\epsilon[M]^2 = \norm{f-f_M}^2 = \sum_{ \abs{\dotp{f}{\psi_m}} \leq T } \abs{\dotp{f}{\psi_m}}^2.$$

If one denotes by $\{ c_R[k] \}_{k=0}^{N-1}$ the set of coefficients magnitudes $\abs{\dotp{f}{\psi_m}}$ ordered by decaying magnitudes, then this error is easily computed as $$\epsilon[M]^2 = \sum_{k=M}^{N-1} c_R[k]^2 = \norm{f}^2 - \sum_{k=0}^{M-1} c_R[k]^2.$$ This means that $\epsilon^2$ is equal to $\norm{f}^2$ minus the discrete primitive of $c_R^2$.

Exercise 2

Compute and display in log scales the ordered coefficients $c_R$. Hint: a discrete primitive can be computed using the function |cumsum|.

In [12]:
exo2()

In [13]:
%% Insert your code here.


Exercise 3

Compute and display in log-scale the non-linear approximation error $\epsilon[M]^2$. Store the values of $\epsilon[M]^2$ in a vector |err_fft|.

In [14]:
exo3()

In [15]:
%% Insert your code here.


## Wavelet Approximation¶

The Wavelet basis of continuous 2-D functions is defined by by scaling and translating three mother atoms $\{\psi^H,\psi^V,\psi^D\}$: $$\psi_{j,n}^k(x) = \frac{1}{2^j}\psi^k\pa{\frac{x-2^j n}{2^j}}$$

Non-linear wavelet approximation is a the heart of the JPEG-2000 compression standard.

The set of inner products $\{ \dotp{f}{\psi_m} \}_m$ is computed in $O(N)$ operations with the 2-D Fast Wavelet Transform algorithm.

Perform a wavelet transform. Here we use a daubechies wavelet transform.

In [16]:
Jmin = 1;
options.h = compute_wavelet_filter('Daubechies',10);
fW = perform_wavortho_transf(f,Jmin,+1, options);


Display the coefficients.

In [17]:
clf;
plot_wavelet(fW,Jmin);
title('Wavelet coefficients');