%%%
% read the image tire.tif
% The image is a part of standard Image Processing Toolbox distribution
%%%
IM    = imread('tire.tif');
[r,c] = size(IM);

scrsz = get(0,'ScreenSize');	% parameters of the screen


%%%
% first equalized by the built-in matlab function histeq
%%%
IMeq1 = histeq(IM);

%%%
% use cumulative distribution
%%%
[H,intensity] = imhist(IM);
cumsumIM	  = cumsum(H);						% cumulative distribution
cumsumIM	  = round(255*cumsumIM/(r*c));		% normalization
IMeq2 = uint8(cumsumIM(double(IM)+1));			% transformation

%%%
% plot the original image and its histogram
%%%
figure(1)
set(1,'Position',[1 (scrsz(4)/2)-75 scrsz(3)/3 scrsz(4)/2]);
clf
subplot(211)
imshow(IM)
title('original image')
subplot(212)
imhist(IM)

%%%
% plot equalized images
% and their histograms
%%%
figure(2)
set(2,'Position',[scrsz(3)/3 (scrsz(4)/2)-75 scrsz(3)/3 scrsz(4)/2]);
clf
subplot(211)
imshow(IMeq1)
title('equalized image by histeq')
subplot(212)
imhist(IMeq1)

figure(3)
set(3,'Position',[2*scrsz(3)/3 (scrsz(4)/2)-75 scrsz(3)/3 scrsz(4)/2]);
clf
subplot(211)
imshow(IMeq2)
title('equalized image by cumsum')
subplot(212)
imhist(IMeq2)

%%%
% compare histograms and their cumulative distributions
%%%
[H1,intensity] = imhist(IMeq1);
cumsumIMeq1	   = cumsum(H1);
cumsumIMeq1	   = round(255*cumsumIMeq1/(r*c));
[H2,intensity] = imhist(IMeq2);
cumsumIMeq2	   = cumsum(H2);
cumsumIMeq2	   = round(255*cumsumIMeq2/(r*c));

figure(4)
set(4,'Position',[1 0 scrsz(3)/3 (scrsz(4)/2)-(2*75)]);
clf
plot(cumsumIM,'-k')
hold on
plot(cumsumIMeq1,'-r')
plot(cumsumIMeq2,'-b')
grid on
title('Cumulative distribution functions of image histograms')
text(120,60,'\bf Original image','Color','black')
text(120,40,' Equalized image by \bf histeq','Color','red')
text(120,20,' Equalized image by \bf cumsum','Color','blue')

%%%
% How do the equalized histograms differ?
% Both have very similar cumulative distribution funcitions
% as it is shown in figure(4)
%
% BUT the number of intensities in images vary
%%%
numOfIntensitiesOrig = sum(H>0)
numOfIntensitiesEq1	 = sum(H1>0)
numOfIntensitiesEq2	 = sum(H2>0)



