Color Image Processing

This numerical tour explores color image processing.

Installing toolboxes and setting up the path.
Color Images
Linear Color Spaces
Non-linear Color Spaces
Principal Component Analysis of Colors
Color Image Denoising
Color Image Compression

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/');

Color Images

We load a color image, which is an [n n 3] cube of data.

name = 'flowers';
n = 256;
options.nbdims = 3;
M = load_image(name, n, options);
M = rescale(M);

Display.

clf;
imageplot(M, 'Color image');
clf;
imageplot(M(:,:,1), 'Red',   2,2,1);
imageplot(M(:,:,2), 'Green', 2,2,2);
imageplot(M(:,:,3), 'Blue',  2,2,3);

Linear Color Spaces

A linear color space is obtained by applying a (3,3) transformation matrix T to the RGB color channels of the image

An example of transformation is the YUV color space, where Y is the luminance and UV are the chrominance.

T = [.299 .587 .114; ...
    -.14713 -.28886 .436; ...
    .615 -.51499 -.10001];
% display
disp('YUV color matrix:');
disp(T);

YUV color matrix:
    0.2990    0.5870    0.1140
   -0.1471   -0.2889    0.4360
    0.6150   -0.5150   -0.1000

The RGB to YUV conversion is obtained by applying the matrix.

U = reshape(M, [n*n 3]);
U = reshape( U*T, [n n 3] );
% display
clf;
imageplot(U(:,:,1), 'Y channel', 1,3,1);
imageplot(U(:,:,2), 'U channel', 1,3,2);
imageplot(U(:,:,3), 'V channel', 1,3,3);

The image can be modified by modifying the YUV representation.

% we lower the chrominance of the image.
U1 = U;
U1(:,:,2:3) = U1(:,:,2:3)/2;

Exercice 1: (the solution is exo1.m) Recover an image M1 from the transformed YUV representation U1.

exo1;

Non-linear Color Spaces

A non-linear color space is obtained by a polar or conical parameterization of a linear color space. The angular coordinate in the plane orthogonal to the first linear axis (which is usually the luminance) is called the Hue, and the radial coordinates is called the

Using a luminance which is the sum of the 3 coordinates, one obtain a color system that is quite close to the HSV color system (which has a more complicated definition, but leads to similar results).

First we compute the "Lightness" coordinate

L = sum3(M, 3) / sqrt(3);

The we compute the projection on the plane orthogonal to [1 1 1], for instance using the projection on the two orthognoal unit vectors [0 1 -1]/sqrt(2) and [2 -1 -1]/sqrt(6).

A = ( M(:,:,2)-M(:,:,3) )/sqrt(2);
B = ( 2*M(:,:,1) - M(:,:,2) - M(:,:,3) )/sqrt(6);

The (L,A,B) components are obtained from RGB using a transformation with an orthogonal matrix T.

T = [1/sqrt(3) 1/sqrt(3) 1/sqrt(3); 0 1/sqrt(2) -1/sqrt(2); 2/sqrt(6) -1/sqrt(6) -1/sqrt(6)];
disp('Matrix T');
disp(T);

Matrix T
    0.5774    0.5774    0.5774
         0    0.7071   -0.7071
    0.8165   -0.4082   -0.4082

The Hue/Saturation are the polor coordinates within this plane.

S = sqrt( A.^2 + B.^2 );
H = atan2(B,A);

Display.

clf;
imageplot(L, 'Lightness', 1,3,1);
imageplot(H, 'Hue', 1,3,2);
imageplot(S, 'Saturation', 1,3,3);

We can modify the image by modifying the coordinates. For instance we rotate by a small amount (1/20) the hue of the image.

H1 = H + .05*2*pi();

Exercice 2: (the solution is exo2.m) Reconstruct a modified image using the hue H1.

exo2;

Principal Component Analysis of Colors

Instead of using a fixed color space, one can computed a color space adapted to an image to process.

Each pixel M(i,j,:) defines a 3D point. The color image thus corresponds to a 3D point clouds whose structure reflects the colors within the image.

% number of displayed points
p = 5000;
% randomly selected points
H = reshape(M, [n*n 3]);
sel = randperm(n*n); sel = sel(1:p);
H = H(sel,:);
% display the points
clf;
plot3(H(:,1), H(:,2), H(:,3), 'k.');
set_label('R', 'G', 'B');

The shape of this point cloud can be approximately captured with a 3D ellipsoid whose axis are the eigenvectors of the covariance matrix of the points cloud.

H = reshape(M, [n*n 3]);
% average color (center of the point cloud
m = mean(H, 1);
% centers the point cloud
H = H - repmat(m, [n*n 1]);
% covariance matrix
C = (H'*H)/(n*n);

The ellipsoid direction are obtained with a eigen decomposition of the symmetric covariance matrix, so that C = V*D*V' and D is a diagonal matrix.

[V,D] = eig(C);
D = diag(D);
% re-order according to decaying eigenvalues
[D,I] = sort(D);
if D(1)>D(length(D))
    D = reverse(D); I = reverse(I);
end
V = V(:,I);
% display the matrix
disp('PCA color matrix:');
disp(V);

PCA color matrix:
    0.1425   -0.7314   -0.6669
   -0.2838    0.6153   -0.7354
    0.9482    0.2941   -0.1199

The PCA color space is image-dependant and correponds to a whitening of the point cloud, so that the new normalized image has zero mean and diagonal covariance.

H = reshape(M, [n*n 3]);
% remove mean
H = H - repmat(m, [n*n 1]);
% projection onto the components
H = H*V;
% new transfored image
U = reshape(H, [n n 3]);

The image, withing this new color space, can be manipulated. Here we lower the chrominance of the image. The result is slighly less visible than with YUV modification because the color space is more adapted to the image.

U1 = U;
U1(:,:,2:3) = U1(:,:,2:3)/2;

Exercice 3: (the solution is exo3.m) Inverse the PCA transformation in order to retrive an image M1 from the modified PCA representation U1.

exo3;

Color Image Denoising

Color image denoising is more difficult than grayscale denoising because independant denoising of RGB channels introduces color artifact. For image with homogenous colors, it can be avoided by using an appropriate color spacE.

A noisy image color image is corrupted by a color Gaussian noise.

% the clean image
M0 = rescale(M,.05,.95);
% add gaussian noise
sigma = .13;
Mnoisy = M0 + randn(n,n,3)*sigma;
% display
clf;
imageplot(M0, 'Original image', 1,2,1);
imageplot(clamp(Mnoisy), 'Noisy image', 1,2,2);

Exercice 4: (the solution is exo4.m) Compare (translation invariant) wavelet denoising of color image in the RGB and PCA color space (the PCA space should be estimated from the noisy image).

exo4;

Exercice 5: (the solution is exo5.m) For a complicated, non-homogenous color image, compute a local color space for each pixel, by performing a PCA over group of pixels. Use this adaptive color model to perform wavelet denoising.

exo5;

Color Image Compression

Similarely to color image denoising, color image compression is difficult because of color artifacts.

Exercice 6: (the solution is exo6.m) Compare wavelet domain image compression (quantization+coding) over the original RGB space and a more adapted global or local color space.

exo6;