Grayscale and Color Histogram Manipulations
This tour explore the manipulation (equalization, interpolation) of 1D and 3D histograms.
Contents
This tour considers discrete 1D or 3D data, that are viewed as point clouds in 1D and 3D. Such a point clouds should be understood as a discrete statistical distribution, and should be interpreted as a particular sample from an (unknown) continuous 1D or 3D distribution. In 1D, a useful representation for display is a quantized histogram, but we will always manipulate the original point cloud.
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/');
1D Distributions
A (unordered) 1D set of values is represented using a 1D histogram.
Load an image.
n = 256;
f0 = rescale( load_image('lena', n) );
Display it.
clf; imageplot(f0);
Number of bins.
p = 50;
Compute the histogram.
[h,t] = hist(f0(:), p);
Each entry of h is the number of samples that fails into a bin of size 1/p. To make this curve an approximation of a continuous distribution, it should be normalized by p/n^2.
Display this normalized histogram.
clf;
bar(t,h*p/n^2); axis('tight');
Exercice 1: (the solution is exo1.m) Compute and display the histogram for an increasing number of bins.
exo1;
Load another image.
f1 = rescale( mean(load_image('fingerprint', n),3) );
Display it.
clf; imageplot(f1);
Exercice 2: (the solution is exo2.m) Compare the two histograms.
exo2;
1D Histogram Matching
Histogram matching corresponds to computing the orthogonal projection of a point cloud (a set of values) onto the constraint that its distribution should be equal to a given one.
The name "histogram" equalization is thus a slight abuse since one actually equalize directly the set of values and not quantized histogram. But of course since the value are equals after equalization, the quantized histograms are also equal.
The projected point cloud is equal to the given one up to a permutation. One thus compute the optimal permutation so that the original point cloud is as close as possible (as measured using the L2 norm) from the given one.
It turns out that this projection can be computed very easily for 1D distribution by simply sorting the values.
[tmp,I0] = sort(f0(:)); [tmp,I1] = sort(f1(:)); f0eq(I0) = f1(I1); f0eq = reshape(f0eq, [n n]);
Check the new histogram.
clf; [h,t] = hist(f0eq(:), p); bar(t,h*p/n^2);
Compare before/after equalization.
clf; imageplot(f0, 'Before', 1,2,1); imageplot(f0eq, 'After', 1,2,2);
One can compute a partial matching by linearly interpolating between the matched values.
The interpolation parameter rho in [0,1]. For rho=0.5, this corresponds to a midway equalization.
rho = 0.5;
Perform the linear blending.
f0eq1(I0) = rho*f1(I1) + (1-rho)*f0(I0); f0eq1 = reshape(f0eq1, [n n]);
The progressive equalization.
clf; imageplot(f0, '\rho=0', 1,3,1); imageplot(f0eq1, '\rho=0.5', 1,3,2); imageplot(f0eq, '\rho=1', 1,3,3);
Exercice 3: (the solution is exo3.m) Display the progression of the interpolation of the histograms.
exo3;
2D Distributions
Here we study a simple setting of 2D sets of points, that can be easily visualized.
Number of points.
q = 300;
First distribution of points, inside a square.
d0 = rand(2,q)-.5;
Second distribution of points, inside an anulus.
theta = 2*pi*rand(1,q); r = .8 + .2*rand(1,q); d1 = [cos(theta).*r; sin(theta).*r];
Shortcut for displaying point clouds.
plotp = @(x,c)plot(x(1,:), x(2,:), c);
Display.
clf; hold on; plotp(d0, 'r.'); plotp(d1, 'b.'); axis('off'); axis('equal');
2D Histogram Matching
On the contrary to 1D distribution, matching 2D distributions is difficult because of the lack of ordering.
An iterative algorithm, initially proposed by Marc Bernot, progressively moves the first point cloud toward the second using 1D matching along randomized directions. The final point clouds at convergence realizes the matching.
This algorithm is detailed in
Wasserstein Barycenter and its Application to Texture Mixing Julien Rabin, Gabriel Peyré, Julie Delon and Marc Bernot, Preprint Hal-00476064, April 2010.
Initialize the iterations, the point cloud d1new will be moved toward d2.
d1new = d1;
Compute an orthogonal coordinate system for this iteration.
[U,R] = qr(randn(2));
Transform the point cloud in the coordinate system.
d0U = U*d0; d1U = U*d1;
Perform the equalization along each axis.
for c=1:size(d1,1) d1U(c,:) = perform_hist_eq(d1U(c,:),d0U(c,:)); end
Gradient step size. Should be less than 1. A small step size provides a better matching, but requires more iterations until convergence.
tau = .2;
Retrieve back the position in the original coordonate system.
d1new = (1-tau)*d1new + tau*U'*d1U;
Display the old, new and target position.
clf; hold on; plotp(d0, 'r.'); plotp(d1, 'b.'); plotp(d1new, 'bo'); axis('off'); axis('equal');
Exercice 4: (the solution is exo4.m) Iterate many time the randomized assignement until convergence of d1new toward d2. The random projector U should be re-newed at each iteration.
exo4;
Display final configuration with the matching.
clf; hold on; plotp(d0, 'r.'); plotp(d1, 'b.'); plotp(d1new, 'bo'); plot([d1(1,:);d1new(1,:)], [d1(2,:);d1new(2,:)], 'k'); axis('off'); axis('equal');
An average point cloud is obtained by interpolation.
rho = .5; d1av = (1-rho)*d1 + rho*d1new;
Display.
clf; hold on; plotp(d0, 'r.'); plotp(d1, 'b.'); plotp(d1av, 'g.'); axis('off'); axis('equal');
Exercice 5: (the solution is exo5.m) Show the progressive interpolation.
exo5;
3D Histogram Matching
Load a color image.
n = 256;
f0 = rescale( load_image('hibiscus', n) );
Display.
clf; imageplot(f0);
The distribution of 3D points associated with the image.
d0 = reshape(f0, [n*n 3])';
Shortcut for 3D display of point cloud.
p = 5000;
sel = randperm(n*n); sel = sel(1:p);
plotp3 = @(x,c)plot3(x(1,sel),x(2,sel),x(3,sel), '.');
Display as a point cloud a sub-sample.
clf; plotp3(d0); axis([0 1 0 1 0 1]);
Load a second image.
f1 = rescale( load_image('flowers', n) );
d1 = reshape(f1, [n*n 3])';
Display the image.
clf; imageplot(f1);
Display the point cloud.
clf; plotp3(d1); axis([0 1 0 1 0 1]);
Exercice 6: (the solution is exo6.m) Perform the equalization of each of the coordinate independantly of d1 with d0 to obtain d1new.
exo6;
Retrieve the image from the new point cloud.
f1new = reshape(d1new', [n n 3]);
Display it.
clf; imageplot(f1new);
Equalizing channel by channel does not leads to an image whose distribution is equal to the target distribution (only the projection on the axis are equal).
To obtain a real matching, one can use the iterative projection algorithm.
Exercice 7: (the solution is exo7.m) Perform the iterative randomized projection, display the point cloud after convergence.
exo7;
Retrieve and display the corresponding image.
f1new = reshape(d1new', [n n 3]); clf; imageplot(f1new);
