Image Approximation with Fourier and Wavelets

\[ \newcommand{\NN}{\mathbb{N}} \newcommand{\CC}{\mathbb{C}} \newcommand{\GG}{\mathbb{G}} \newcommand{\LL}{\mathbb{L}} \newcommand{\PP}{\mathbb{P}} \newcommand{\QQ}{\mathbb{Q}} \newcommand{\RR}{\mathbb{R}} \newcommand{\VV}{\mathbb{V}} \newcommand{\ZZ}{\mathbb{Z}} \newcommand{\FF}{\mathbb{F}} \newcommand{\KK}{\mathbb{K}} \newcommand{\UU}{\mathbb{U}} \newcommand{\EE}{\mathbb{E}} \newcommand{\Aa}{\mathcal{A}} \newcommand{\Bb}{\mathcal{B}} \newcommand{\Cc}{\mathcal{C}} \newcommand{\Dd}{\mathcal{D}} \newcommand{\Ee}{\mathcal{E}} \newcommand{\Ff}{\mathcal{F}} \newcommand{\Gg}{\mathcal{G}} \newcommand{\Hh}{\mathcal{H}} \newcommand{\Ii}{\mathcal{I}} \newcommand{\Jj}{\mathcal{J}} \newcommand{\Kk}{\mathcal{K}} \newcommand{\Ll}{\mathcal{L}} \newcommand{\Mm}{\mathcal{M}} \newcommand{\Nn}{\mathcal{N}} \newcommand{\Oo}{\mathcal{O}} \newcommand{\Pp}{\mathcal{P}} \newcommand{\Qq}{\mathcal{Q}} \newcommand{\Rr}{\mathcal{R}} \newcommand{\Ss}{\mathcal{S}} \newcommand{\Tt}{\mathcal{T}} \newcommand{\Uu}{\mathcal{U}} \newcommand{\Vv}{\mathcal{V}} \newcommand{\Ww}{\mathcal{W}} \newcommand{\Xx}{\mathcal{X}} \newcommand{\Yy}{\mathcal{Y}} \newcommand{\Zz}{\mathcal{Z}} \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} \renewcommand{\phi}{\varphi} \renewcommand{\th}{\theta} \newcommand{\om}{\omega} \newcommand{\Om}{\Omega} \renewcommand{\epsilon}{\varepsilon} \newcommand{\Calpha}{\mathrm{C}^\al} \newcommand{\Cbeta}{\mathrm{C}^\be} \newcommand{\Cal}{\text{C}^\al} \newcommand{\Cdeux}{\text{C}^{2}} \newcommand{\Cun}{\text{C}^{1}} \newcommand{\Calt}[1]{\text{C}^{#1}} \newcommand{\lun}{\ell^1} \newcommand{\ldeux}{\ell^2} \newcommand{\linf}{\ell^\infty} \newcommand{\ldeuxj}{{\ldeux_j}} \newcommand{\Lun}{\text{\upshape L}^1} \newcommand{\Ldeux}{\text{\upshape L}^2} \newcommand{\Lp}{\text{\upshape L}^p} \newcommand{\Lq}{\text{\upshape L}^q} \newcommand{\Linf}{\text{\upshape L}^\infty} \newcommand{\lzero}{\ell^0} \newcommand{\lp}{\ell^p} \renewcommand{\d}{\ins{d}} \newcommand{\Grad}{\text{Grad}} \newcommand{\grad}{\text{grad}} \renewcommand{\div}{\text{div}} \newcommand{\diag}{\text{diag}} \newcommand{\pd}[2]{ \frac{ \partial #1}{\partial #2} } \newcommand{\pdd}[2]{ \frac{ \partial^2 #1}{\partial #2^2} } \newcommand{\dotp}[2]{\langle #1,\,#2\rangle} \newcommand{\norm}[1]{\| #1 \|} \newcommand{\normi}[1]{\norm{#1}_{\infty}} \newcommand{\normu}[1]{\norm{#1}_{1}} \newcommand{\normz}[1]{\norm{#1}_{0}} \newcommand{\abs}[1]{\vert #1 \vert} \newcommand{\argmin}{\text{argmin}} \newcommand{\argmax}{\text{argmax}} \newcommand{\uargmin}[1]{\underset{#1}{\argmin}\;} \newcommand{\uargmax}[1]{\underset{#1}{\argmax}\;} \newcommand{\umin}[1]{\underset{#1}{\min}\;} \newcommand{\umax}[1]{\underset{#1}{\max}\;} \newcommand{\pa}[1]{\left( #1 \right)} \newcommand{\choice}[1]{ \left\{ \begin{array}{l} #1 \end{array} \right. } \newcommand{\enscond}[2]{ \left\{ #1 \;:\; #2 \right\} } \newcommand{\qandq}{ \quad \text{and} \quad } \newcommand{\qqandqq}{ \qquad \text{and} \qquad } \newcommand{\qifq}{ \quad \text{if} \quad } \newcommand{\qqifqq}{ \qquad \text{if} \qquad } \newcommand{\qwhereq}{ \quad \text{where} \quad } \newcommand{\qqwhereqq}{ \qquad \text{where} \qquad } \newcommand{\qwithq}{ \quad \text{with} \quad } \newcommand{\qqwithqq}{ \qquad \text{with} \qquad } \newcommand{\qforq}{ \quad \text{for} \quad } \newcommand{\qqforqq}{ \qquad \text{for} \qquad } \newcommand{\qqsinceqq}{ \qquad \text{since} \qquad } \newcommand{\qsinceq}{ \quad \text{since} \quad } \newcommand{\qarrq}{\quad\Longrightarrow\quad} \newcommand{\qqarrqq}{\quad\Longrightarrow\quad} \newcommand{\qiffq}{\quad\Longleftrightarrow\quad} \newcommand{\qqiffqq}{\qquad\Longleftrightarrow\qquad} \newcommand{\qsubjq}{ \quad \text{subject to} \quad } \newcommand{\qqsubjqq}{ \qquad \text{subject to} \qquad } \]

This numerical tour overviews the use of Fourier and wavelets for image approximation.

Contents

Installing toolboxes and setting up the path.

You need to download the following files: signal toolbox and general toolbox.

You need to unzip these toolboxes in your working directory, so that you have toolbox_signal and toolbox_general in your directory.

For Scilab user: you must replace the Matlab comment '%' by its Scilab counterpart '//'.

Recommandation: You should create a text file named for instance numericaltour.sce (in Scilab) or numericaltour.m (in Matlab) to write all the Scilab/Matlab command you want to execute. Then, simply run exec('numericaltour.sce'); (in Scilab) or numericaltour; (in Matlab) to run the commands.

Execute this line only if you are using Matlab.

getd = @(p)path(p,path); % scilab users must *not* execute this

Then you can add the toolboxes to the path.

getd('toolbox_signal/');
getd('toolbox_general/');

Note: to measure the error of an image \(f\) with its approximation \(f_M\), we use the SNR measure, defined as \[ \text{SNR}(f,f_M) = -20\log_{10} \pa{ \frac{ \norm{f-f_M} }{ \norm{f} } }, \] which is a quantity expressed in decibels (dB). The higer the SNR, the better the quality.

Image Loading and Displaying

First we load an image \( f \in \RR^N \) of \( N = N_0 \times N_0 \) pixels.

name = 'lena';
n0 = 512;
f = rescale( load_image(name,n0) );

Display the original image.

clf;
imageplot( f, 'Image f');

Display a zoom in the middle.

clf;
imageplot( crop(f,64), 'Zoom' );

An image is a 2D array, that can be modified as a matrix.

clf;
imageplot(-f, '-f', 1,2,1);
imageplot(f(n0:-1:1,:), 'Flipped', 1,2,2);

Blurring is achieved by computing a convolution \(f \star h\) with a kernel \(h\).

Compute the low pass kernel.

k = 9; % size of the kernel
h = ones(k,k);
h = h/sum(h(:)); % normalize

Compute the convolution \(f \star h\).

fh = perform_convolution(f,h);

Display.

clf;
imageplot(fh, 'Blurred image');

Fourier Transform

The Fourier orthonormal basis is defined as \[ \psi_m(k) = \frac{1}{\sqrt{N}}e^{\frac{2i\pi}{N_0} \dotp{m}{k} } \] where \(0 \leq k_1,k_2 < N_0\) are position indexes, and \(0 \leq m_1,m_2 < N_0\) are frequency indexes.

The Fourier transform \(\hat f\) is the projection of the image on this Fourier basis \[ \hat f(m) = \dotp{f}{\psi_m}. \]

The Fourier transform is computed in \( O(N \log(N)) \) operation using the FFT algorithm (Fast Fourier Transform). Note the normalization by \(\sqrt{N}=N_0\) to make the transform orthonormal.

F = fft2(f) / n0;

We check this conservation of the energy.

disp(strcat(['Energy of Image:   ' num2str(norm(f(:)))]));
disp(strcat(['Energy of Fourier: ' num2str(norm(F(:)))]));

Energy of Image:   255.9831
Energy of Fourier: 255.9831

Compute the logarithm of the Fourier magnitude \( \log(\abs{\hat f(m)} + \epsilon) \), for some small \(\epsilon\).

L = fftshift(log( abs(F)+1e-1 ));

Display. Note that we use the function fftshift is useful to put the 0 low frequency in the middle.

clf;
imageplot(L, 'Log(Fourier transform)');

Linear Fourier Approximation

An approximation is obtained by retaining a certain set of index \(I_M\) \[ f_M = \sum_{ m \in I_M } \dotp{f}{\psi_m} \psi_m. \]

Linear approximation is obtained by retaining a fixed set \(I_M\) of \(M = \abs{I_M}\) coefficients. The important point is that \(I_M\) does not depend on the image \(f\) to be approximated.

For the Fourier transform, a low pass linear approximation is obtained by keeping only the frequencies within a square. \[ I_M = \enscond{m=(m_1,m_2)}{ -q/2 \leq m_1,m_2 < q/2 } \] where \( q = \sqrt{M} \).

This can be achieved by computing the Fourier transform, setting to zero the \(N-M\) coefficients outside the square \(I_M\) and then inverting the Fourier transform.

Number \(M\) of kept coefficients.

M = n0^2/64;

Exercice 1: (the solution is exo1.m) Perform the linear Fourier approximation with \(M\) coefficients. Store the result in the variable fM.

exo1;

Compare two 1D profile (lines of the image). This shows the strong ringing artifact of the linea approximation.

clf;
subplot(2,1,1);
plot(f(:,n0/2));
axis('tight'); title('f');
subplot(2,1,2);
plot(fM(:,n0/2));
axis('tight'); title('f_M');

Non-linear Fourier Approximation

Non-linear approximation is obtained by keeping the \(M\) largest coefficients. This is equivalently computed using a thresholding of the coefficients \[ I_M = \enscond{m}{ \abs{\dotp{f}{\psi_m}}>T }. \]

Set a threshold \(T>0\).

T = .2;

Compute the Fourier transform.

F = fft2(f) / n0;

Do the hard thresholding.

FT = F .* (abs(F)>T);

Display. Note that we use the function fftshift is useful to put the 0 low frequency in the middle.

clf;
L = fftshift(log( abs(FT)+1e-1 ));
imageplot(L, 'thresholded Log(Fourier transform)');

Inverse Fourier transform to obtained \(f_M\)

fM = real( ifft2(FT)*n0 );

Display.

clf;
imageplot(clamp(fM), ['Non-linear, Fourier, SNR=' num2str(snr(f,fM), 4) 'dB']);

Given a \(T\), the number of coefficients is obtained by counting the non thresholded coefficients \( \abs{I_M} \).

m = sum(FT(:)~=0);
disp(['M/N = 1/'  num2str(round(n0^2/m)) '.']);

M/N = 1/32.

Exercice 2: (the solution is exo2.m) Compute the value of the threshold \(T\) so that the number of coefficients is \(M\). Display the corresponding approximation \(f_M\).

exo2;

Wavelet Transform

A wavelet basis \( \Bb = \{ \psi_m \}_m \) is obtained over the continuous domain by translating an dilating three mother wavelet functions \( \{\psi^V,\psi^H,\psi^D\} \).

Each wavelet atom is defined as \[ \psi_m(x) = \psi_{j,n}^k(x) = \frac{1}{2^j}\psi^k\pa{ \frac{x-2^j n}{2^j} }. \] The scale (size of the support) is \(2^j\) and the position is \(2^j(n_1,n_2)\). The index is \( m=(k,j,n) \) for \{ j \leq 0 \}.

The wavelet transform computes all the inner products \( \{ \dotp{f}{\psi_{j,n}^k} \}_{k,j,n} \).

Set the minimum scale for the transform to be 0.

Jmin = 0;

Perform the wavelet transform, fw stores all the wavelet coefficients.

fw = perform_wavelet_transf(f,Jmin,+1);

Display the transformed coefficients.

clf;
plot_wavelet(fw);

Wavelet Approximation

Linear wavelet approximation with \(M=2^{-j_0}\) coefficients is obtained by keeping only the coarse scale (large support) wavelets: \[ I_M = \enscond{(k,j,n)}{ j \geq j_0 }. \]

It corresponds to setting to zero all the coefficients excepted those that are on the upper left corner of fw.

Exercice 3: (the solution is exo3.m) Perform linear approximation with \(M\) wavelet coefficients.

exo3;

A non-linear approximation is obtained by keeping the \(M\) largest wavelet coefficients.

As already said, this is equivalently computed by a non-linear hard thresholding.

Select a threshold.

T = .2;

Perform hard thresholding.

fwT = fw .* (abs(fw)>T);

Display the thresholded coefficients.

clf;
subplot(1,2,1);
plot_wavelet(fw);
title('Original coefficients');
subplot(1,2,2);
plot_wavelet(fwT);

Perform reconstruction.

fM = perform_wavelet_transf(fwT,Jmin,-1);

Display approximation.

clf;
imageplot(clamp(fM), strcat(['Approximation, SNR=' num2str(snr(f,fM),3) 'dB']));

Exercice 4: (the solution is exo4.m) Perform non-linear approximation with \(M\) wavelet coefficients by chosing the correct value for \(T\). Store the result in the variable fM.

exo4;

Compare two 1D profile (lines of the image). Note how the ringing artifacts are reduced with respec to the Fourier approximation.

clf;
subplot(2,1,1);
plot(f(:,n0/2));
axis('tight'); title('f');
subplot(2,1,2);
plot(fM(:,n0/2));
axis('tight'); title('f_M');

Copyright (c) 2010 Gabriel Peyre