%
% HistEQ.m
%
%   Perform histogram equalization on a gray scale byte-per-pixel image.
%
%   Usage: 
%   >> K = HistEQ(x);
%
%   Input parameters:
%     x          double array, holds the original image
%
%   Output parameters:
%     K          double array, holds the equalized image
%
% 4/3/03 jph
% 2/3/15: improvements for faster run time. jph
% 

function [K] = HistEQ(x)

%
% Compute the histogram of the input image
%
InputHist = sum(hist(x,0:255)');

%
% Compute the normalized histogram of the input image
%
[Nrows,Ncols] = size(x);
Npix = Nrows * Ncols;
One_on_Npix = 1 / Npix;
p_I = InputHist * One_on_Npix;

%
% Compute the CDF (distribution)
%
P_I = zeros(1,256);
P_I(1) = p_I(1);
for i=2:256
  P_I(i) = P_I(i-1) + p_I(i);
end

%
% Compute the intermediate image J.  It is the CDF of x.
%
J = zeros(Nrows,Ncols);
for row=1:Nrows
  for col=1:Ncols
    J(row,col) = P_I(x(row,col)+1);
  end
end

%
% The output image K is obtained by applying a
%     full-scale contrast stretch to J.
%
K = stretch(J);
