How does Histogram Equalization work?

I need to be able to perform histogram equalization without using the built-in MATLAB function. This would not be an issue if I could understand what it is and how it works.

I thought it was multiplying each value in the grayscale image array by the corresponding probability mass function in that array based off the histogram bin and probability density functions.

Is that correct? If so, how do I code that? If not (more likely in my opinion), what am I not understanding? What would best help without using the built-in function?

I attached my entire project code just so you could have access to my data, if it helps. The section under "The Part I Need Help With" is the relevant portion.

[DGM]

I can't open .mlx files in any practical way, so this isn't really a direct answer. I'm just going to point out that histeq() can always be opened and inspected. Alternatively, MIMT contains a replacement that may (or may not) be easier to read, though it relies on other MIMT tools. It should still be fairly easy to read between the lines though.

function [outpict TF] = histeqFB(inpict,varargin)
%   OUTPICT=HISTEQFB(INPICT,{NBINS})
%   OUTPICT=HISTEQFB(INPICT,{HGRAM})
%   [OUTPICT TF]=HISTEQFB( ... )
%   Adjust the contrast of an image such that its intensity distribution conforms to that described
%   by a given histogram.  If no histogram is given explicitly, a flat (uniform) histogram of either
%   a default or user-specified number of bins is used.
% 
%   Optionally, the 1D transfer function TF mapping the input image to the output image can be returned.
%
%   This is a passthrough to the IPT function histeq(), with an internal fallback implementation 
%   to help remove the dependency of MIMT tools on the Image Processing Toolbox. As with other fallback 
%   tools, performance without IPT may be degraded or otherwise different due to the methods used.  
%
%   INPICT is a single-channel image of any standard image class. 
%      While the IPT tools support indexed images, the fallback tools do not. 
%   NBINS optionally specifies how many bins should be used for flat histograms (default 64)
%   HGRAM optionally specifies the target output distribution explicitly
%
%  Output image class is inherited from INPICT
%  TF class is double
%
% Webdocs: http://mimtdocs.rf.gd/manual/html/histeqFB.html
% See also: histeq, imhist, imhistFB, imlnc

% IF IPT IS INSTALLED, USE IT
if hasipt()
    [outpict TF] = histeq(inpict,varargin{:});
    return;
end

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

nin = 256; % internal default

if numel(varargin) == 0
    % HISTEQFB(I)
    nout = 64;
    hgram = ones([1 nout])*numel(inpict)/nout;
else
    if size(varargin{1},2) == 3
        % HISTEQFB(IX,MAP)
        % HISTEQFB(IX,MAP,HGRAM)
        error('HISTEQFB: fallback methods to not support indexed images')
    elseif numel(varargin{1}) == 1
        % HISTEQFB(I,NBINS)
        nout = round(varargin{1});
        hgram = ones([1 nout])*numel(inpict)/nout;
    else
        % HISTEQFB(I,HGRAM)
        hgram = round(varargin{1});
        hgram = hgram*numel(inpict)/sum(hgram); % normalize wrt sum
        nout = numel(hgram);
        if ~isvector(hgram)
            error('HISTEQFB: HGRAM is supposed to be a vector')
        end
    end
end

% get input histogram
hgramin = imhistFB(inpict,nin)'; % same as IPT imhist()

% calculate cumulative hgram sums
incsum = cumsum(hgramin); 
outcsum = cumsum(hgram*numel(inpict)/sum(hgram));

% generate TF mapping inpict to outpict
numin = numel(inpict);
tol = ones([nout 1])*min([hgramin(1:nin-1) 0; 0 hgramin(2:nin)])/2;
err = (outcsum(:)*ones([1 nin])-ones([nout 1])*incsum(:)')+tol;
orv = find(err < -numin*sqrt(eps));
if ~isempty(orv)
    err(orv) = numin*ones(size(orv));
end
[~,idx] = min(err); % location of minimum error
TF = (idx-1)/(nout-1); % normalize

% this replicates the functionality of grayxformmex
[inpict inclass] = imcast(inpict,'double'); % imcast() replaces im2double(),im2uint8(), etc
inpict = min(max(inpict,0),1);
outpict = interp1(linspace(0,1,numel(TF)),TF,inpict,'nearest');
outpict = imcast(outpict,inclass);

I doubt you need to make your implementation mimic all the different functionality of histeq(), so some simplification is probably due.

Answer 1

[Udit06]

Hi Aiden,

I understand that you want to perform histogram equalization without using the inbuilt MATLAB function and want to know the steps to perform the same. Let's assume that the image has L total intensity levels from 0 to L-1 (for an 8-bit image, L=256). Let r be the input pixel intensities in the original image, and s be the output pixel intensities in the equalized image.

The transformation T that maps the input pixel value r to the output pixel value s is given by:

where CDF(r) is the cumulative distribution function of the input image's histogram at the pixel intensity level r.

You can refer to the page 28 of the following pdf resource to understand more about histogram equalization.

https://www.imageprocessingplace.com/downloads_V3/dipum3e_downloads/dipum3e_sample_book_material/dipum3e_Sample%20Book%20Material.pdf

I hope this helps.