Image Approximation with Orthogonal Bases
This numerical tour uses several orthogonal bases to perform image approximation.
Contents
Installing toolboxes and setting up the path.
You need to download the general purpose toolbox and the signal toolbox.
You need to unzip these toolboxes in your working directory, so that you have toolbox_general/ and toolbox_signal/ 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 commands 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(path,p); % scilab users must *not* execute this
Then you can add these toolboxes to the path.
% Add some directories to the path getd('toolbox_signal/'); getd('toolbox_general/');
Fourier Approximation
Approximation with Fourier basis is not efficient because it is too global, and also it treats the boundary of the image using periodic condition.
First we load an image.
n = 256;
M = rescale( load_image('lena', n) );
Compute the Fourier transform using the FFT algorithm.
MF = fft2(M);
Display its magnitude. We use the function fftshift to put the low frequency in the center.
clf; imageplot(M, 'Original image', 1,2,1); imageplot(log(1e-5+abs(fftshift(MF))), 'log(FFT)', 1,2,2);
Exercice 1: (the solution is exo1.m) Compute a best m-term approximation in the Fourier basis of M (use a well chosen hard thresholding). You can use several m values, for instance m=.02*n^2 and m=.1*n^2. Use the function ifft2 to recover an image from the coefficients.
exo1;
Best m-term approximation error err_fft(m) by summing the v(i+1:end) where v is the vector of decreasing ordered squared coefficients.
v = sort(abs(MF(:)).^2); if v(n^2)<v(1) v = reverse(v); % be sure it is in reverse order end err_fft = reverse( cumsum(v) ); err_fft = err_fft/err_fft(1);
Cosine Approximation
To avoid periodicity approximation artifacts, one can symmetrize the image along the horizontal and vertical directions, and perform a Fourier approximation. Since the image is symmetric, one can only consider 1/4 of the (2*n)^2 sets of symmetric coefficients, which corresponds to the decomposition in a Discrete Cosinde orthogonal basis.
Compute the Cosine transform using the fast DCT algorithm.
MC = dct2(M);
Display its magnitude. Note that the low frequencies are now in the upper-left corner.
clf; imageplot(M, 'Original image', 1,2,1); imageplot(log(1e-5+abs(MC)), 'log(DCT)', 1,2,2);
Exercice 2: (the solution is exo2.m) Compute a few best m-term approximations in the DCT basis of M. Use the function ifft2 to recover an image from the coefficients.
exo2;
Exercice 3: (the solution is exo3.m) Record in err_dct(m) the error of the m-term DCT approximation.
exo3;
Local Cosine Approximation
To improve the global DCT approximation, one can approximate independantly small patches in the image. This corresponds to a decomposition in a local cosine basis, which is at the heart of the JPEG image compression standard.
The only parameter of the transform is the size of the square.
w = 16;
Initialize at zero the transformed image in the local DCT basis.
ML = zeros(n,n);
Example of patch index.
i = 5; j = 7;
For a given path index (i,j), we extract an (w,w) patch.
seli = (i-1)*w+1:i*w; selj = (j-1)*w+1:j*w; P = M(seli,selj);
Compute the Cosine transform of the patch using the fast DCT algorithm.
ML(seli,selj) = dct2(P);
Display the patch and its coefficients. We removed the low frequency of P for display purpose only.
clf; imageplot(P,'Patch',1,2,1); imageplot(dct2(P-mean(P(:))),'DCT',1,2,2);
Exercice 4: (the solution is exo4.m) Compute the local DCT transform ML by transforming each patch.
exo4;
Display the coefficients.
clf; imageplot(M,'Image',1,2,1); imageplot(min(abs(ML),.005*w*w),'Local DCT',1,2,2);
Exercice 5: (the solution is exo5.m) Compute the inverse local DCT transform M1 of ML by inverse transforming each patch using the function idct2.
exo5;
Error |M-M1|/|M| = 3.7263e-16
Exercice 6: (the solution is exo6.m) Compute a few best m-term approximations in the Local DCT basis of M.
exo6;
Exercice 7: (the solution is exo7.m) Record in err_ldct(m) the error of the m-term DCT approximation.
exo7;
Wavelet Transform of an Image
Local cosine approximation creates blocking artifact for aggressive approximation. This can be avoided using a wavelet basis which provide a multiscale decomposition of the image.
Perform a wavelet transform. Here we use a 7/9 biorthogonal wavelet transform, which is nearly orthogonal and is used in the JPEG-2000 image compression standard.
Jmin = 2; MW = perform_wavelet_transf(M,Jmin,+1);
Display the coefficients.
clf; subplot(1,2,1); imageplot(M, 'Image'); subplot(1,2,2); plot_wavelet(MW,Jmin); title('Wavelet coefficients');
Exercice 8: (the solution is exo8.m) Compute a few best m-term approximations in the wavelet basis of M. Use the function perform_wavelet_transf(MWT,Jmin,-1) to recover an image from the thresholded coefficients MWT.
exo8;
Exercice 9: (the solution is exo9.m) Record in err_wav(m) the error of the m-term wavelet approximation.
exo9;
Comparison of the Orthogonal Bases
An orthognoal bases is better than an other one for approximating a given image if the m-term approximation error decay faster.
Exercice 10: (the solution is exo10.m) Display in log-log plot a zoom on a sub-set of the error decay curve for the different orthogonal bases.
exo10;
Comparison of Wavelet Approximations of Several Images
An image is more complicated than an other one for a given orthogonal basis if its approximation error decays more slowly.
First load several high resolution images.
n = 512;
name = {'regular3','phantom','lena','mandrill'};
Mlist = rescale( load_image(name,n) );
Display them.
clf; imageplot(Mlist,name,2,2);
Exercice 11: (the solution is exo11.m) Compare the approximation error decay for those images.
exo11;
