%
%   Inew = histeq( I ) 
%
%      - performs a histogram equalization on I
%
%      - number of gray levels is estimated as the nearest
%        power of 2 above the highest intensity in I
%
%   Inew = histeq( I, G )
%
%      - G specifies the maximum gray level, eg. 255
%
function Inew = histeq(I, G)

	% If no max gray level is specified, make a guess
	if nargin == 1
		G = 2^ceil(log(max(max(I)))/log(2)) - 1
	end
	
	% Compute the transfer function T
	H = hist(I(:), 0:G);
	CP = cumsum(H)./sum(H);
	T = round((CP-min(CP))/(1-min(CP))*G);
	
	% Use the transfer function to compute the equalized image
	Inew = zeros(size(I));
	Inew(:) = T(round(I(:))+1);
end