Histogram equalisation process and errors

Iam working on implenting a histogram equalisation of an image without using inbuilt functions specifcally imhist and histeq, I have, for the most part understood what do during a histogram equalisation i.e. find the histogram, normalize the histogram values, find the cummulative probabillity and then map the old image values to the new ones.

Can someone look at this code and help me fix what is wrong? Iam allowed to use cumsum and Hist or Histcounts. I have a test function that checks my code for different values that I've added for reference. Thank you.

This is my code:

 %input image data is assumed to be in range 0..1

I = imread(img);
image =I;
[m,n] = size(image);
L = 256;
H = histcounts(image(:),(0:256));
H = H.';
[counts] = H;
x = 0:255;
myCDF = cumsum(counts)/(m*n);
equalizedI = (L-1)*myCDF(double(image)+1);
equalizedI = uint8(equalizedI);
histMyOut = histcounts(equalizedI,256);

eq_img = histMyOut
return

My error:

Array indices must be positive integers or logical values.
Error in histeq_contrast (line 22)
equalizedI = (L-1)*myCDF(double(image)+1);
Error in histeq_test (line 16)
I1eq = histeq_contrast(I1);

For reference, my test function is:

I1 = imread('vintage_postcard.tif');
%I1 = imread('pout.tif');
%I1 = imread('office_5.jpg'); I1=I1(:,:,2);
%I1 = imread('office_1.jpg'); I1=I1(:,:,2);
%I1 = imread('low_light.tif');
I1=im2double(I1);

% damage contrast
I1=0.8*I1;

I1eq = histeq_contrast(I1);
figure
subplot(1,2,1);
imagesc(I1); caxis([0 1]); title('Test Image 1'); axis equal tight
subplot(1,2,2);
imagesc(I1eq); caxis([0 1]); title('Histeq Result');  axis equal tight
colormap(gray);

fprintf(1,'Min/max of input image  1: %3.3f %3.3f\n', min(I1(:)),max(I1(:)) );
fprintf(1,'Min/max of output image 1: %3.3f %3.3f\n', min(I1eq(:)),max(I1eq(:)) );

% damage contrast
I2 = I1*0.25 + 0.25;

I2eq = histeq_contrast(I2);
figure
subplot(1,2,1);
imagesc(I2); caxis([0 1]); title('Test Image 2');  axis equal tight
subplot(1,2,2);
imagesc(I2eq); caxis([0 1]); title('Histeq Result');  axis equal tight
colormap(gray);

fprintf(1,'Min/max of input image  2: %3.3f %3.3f\n', min(I2(:)),max(I2(:)) );
fprintf(1,'Min/max of output image 2: %3.3f %3.3f\n', min(I2eq(:)),max(I2eq(:)) );

Answer 1

[DGM]

Given the way you're trying to do 1D lookup for the final remapping transformation, your input must be integer-valued and nonzero. If your inputs are unit-scale float, direct lookup by array indexing won't work. You'd need to actually do the interpolation.

That said, so long as your input is strictly uint8, your example works. It can be cleaned up significantly:

% setup
inpict = imread('cameraman.tif');

% transform the image
counts = histcounts(inpict,'binmethod','integers','binlimits',[0 255]);
myCDF = cumsum(counts)/numel(inpict);
outpict = 255*myCDF(double(inpict)+1);
outpict = uint8(outpict);

% test the results
opref = histeq(inpict,256); % the reference
immse(outpict,opref)

ans = 2.8904

imshow([inpict opref outpict])

In this example, we're assuming:

• input is single-channel
• input is uint8 class
• number of input and output bins is 256
• target output distribution is uniform

Note that the example and reference outputs are visually similar, but not identical. They differ by only a few LSB, why? The fact is that because histeq() is designed to be more flexible than this integer-only code, it does the transformation calculations and output interpolation differently.