% histeqdemo.m: histogram equalization step by step demonstration
%
% x0: M x N input image, each pixel has 8 bits
% y0: histogram equalized image, each pixel has 8 bits
% L: [0:L-1] # of gray levels. default L = 256
% 
% Algorithm:
% 1. compute the PDF, and then CDF (T(r))of x0
% 2. point-wise gray-level transform according to T(r)
% (C) copyright 2005 by Yu Hen Hu

% function y0=histeq(x0,L)
% Usage: y0=histeq(x0,L)
% histogram equalization algorithm
% (C) 2002 by Yu Hen Hu
% created: 10/1/2002
%
% x0: M x N input image, each pixel has 8 bits
% y0: histogram equalized image, each pixel has 8 bits
% L: [0:L-1] # of gray levels. default L = 256
% 
% Algorithm:
% 1. compute the PDF, and then CDF (T(r))of x0
% 2. point-wise gray-level transform according to T(r)

% if nargin==1, L=256; end

clear all, close all
load p64int.txt; % 64 x 64 image, L = 256
x0=p64int; 
L=256;

[m,n]=size(x0); 
len=m*n;
x=reshape(x0,len,1);

xpdf=hist(x,[0:L-1]); % pdf, 1 x L
tr=round(xpdf*triu(ones(L))*(L-1)/len);  % cdf, range from 0 to L-1

y0=zeros(m,n);
for i=1:L,
   if xpdf(i)>0,
      y0=y0+[x0==i-1]*tr(i);
   end
end
ypdf=hist(reshape(y0,len,1),[0:L-1]); % pdf of y, 1 x L

figure(1),subplot(211),stem([0:L-1],xpdf),title('histogram, original')
axis([0 256 0 500])

subplot(212),stem([0:L-1],ypdf),title('histogram, equalized'),axis([0 256 0 500])

% figure(2),subplot(121),image(x0),colormap('gray')
% subplot(122),image(y0),colormap('gray')
figure(2),subplot(121),imshow(uint8(x0)),title('before')
subplot(122),imshow(uint8(y0)),title('after')

figure(3)
stairs([0:L-1],tr),title('transformation'),axis([0 256 0 256])