Histogram Equalization method without use of histeq

Question

I am new to Matlab and am trying to implement code to perform the same function as histeq without actual use of the function. In my code the image colour I get changes drastically when it should not change that much. The average intensity in the image (ranging between 0 and 255) is 105.3196. The image is of an open source pollen particle.

Any help would be much appreciated. The sooner the better! Please could any help be simplified as my Matlab understanding is limited. Thanks.

clc;
clear all;
close all;

pollenJpg = imread ('pollen.jpg', 'jpg');
greyscalePollen = rgb2gray (pollenJpg);

histEqPollen = histeq(greyscalePollen);

averagePollen = mean2 (greyscalePollen)

sizeGreyScalePollen = size(greyscalePollen);
rowsGreyScalePollen = sizeGreyScalePollen(1,1);
columnsGreyScalePollen = sizeGreyScalePollen(1,2);

for  i = (1:rowsGreyScalePollen)
    for  j = (1:columnsGreyScalePollen)
        if (greyscalePollen(i,j) > averagePollen)
            greyscalePollen(i,j) = greyscalePollen(i,j) + (0.1 * averagePollen);

            if (greyscalePollen(i,j) > 255)
                greyscalePollen(i,j) = 255;
            end

        elseif (greyscalePollen(i,j) < averagePollen)
            greyscalePollen(i,j) = greyscalePollen(i,j) - (0.1 * averagePollen);

            if (greyscalePollen(i,j) > 0)
                greyscalePollen(i,j) = 0;
            end

        end
    end
end

figure;
imshow (pollenJpg);
title ('Original Image');
figure;
imshow (greyscalePollen);
title ('Attempted Histogram Equalization of Image');
figure;
imshow (histEqPollen);
title ('True Histogram Equalization of Image');

Answer 1

Relatively few (5) lines of code can do this. I used a low contrast file called 'pollen.jpg' that I found at http://commons.wikimedia.org/wiki/File%3ALepismium_lorentzianum_pollen.jpg

I read it in using your code, run all the above, then do the following:

% find out the index of pixels sorted by intensity:
[gv gi] = sort(greyscalePollen(:));  

% create a table of "approximately equal" intensity values:
N = numel(gv);
newVals = repmat(0:255, [ceil(N/256) 1]);

% perform lookup: 
% the pixels in sorted order need new values from "equal bins" table:
newImg = zeros(size(greyscalePollen));
newImg(gi) = newVals(1:N);

% if the size of the image doesn't divide into 256, the last bin will have 
% slightly fewer pixels in it than the others

When I run this algorithm, and then create a composite of the four images (original, your attempt, my attempt, and histeq), you get the following:

I think it's convincing. The images are not exactly identical - I believe that is because the matlab histeq routine ignores all pixels with value 0. Since it is fully vectorized it is also pretty fast (although not nearly as fast as histeq by about a factor 15 on my image.

EDIT: a bit of explanation might be in order. The repmat command I use to create the newVals matrix creates a matrix that looks like this:

0  1  2  3  4  ...  255
0  1  2  3  4  ...  255
0  1  2  3  4  ...  255
...
0  1  2  3  4  ...  255

Since matlab stores matrices in "first index first" order, if you read this matrix with a single index (as I do in the line newVals(1:N)), you access first all the zeros, then all the ones, etc:

0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 3 3 3 3 3 3 3 3 ...

So - when I know the indices of the pixels in the order of their intensity (as returned by the second argument of the sort command, which I called gi), then I can easily assign the value 0 to the first N/256 pixels, the value 1 to the next N/256 etc, with the command I used:

newImg(gi) = newVals(1:N);

I hope this makes the code a little easier to understand.

Answer 2

To implement the equalisation algorithm described on the Wikipedia page, follow these these steps:

Decide on a binSize to group greyscale values. (This is a tweakable, the larger the bin, the less accurate the result from the ideal case, but I think it can cause problems if chosen too small on real images).

Then, calculate the probability of a pixel being a shade of grey:

pixelCount = imageWidth * imageHeight
histogram = all zero
for each pixel in image at coordinates i, j
    histogram[floor(pixel / 255 / 10) + 1] += 1 / pixelCount // 1-based arrays, not 0-based
    // Note a technicality here: you may need to 
    // write special code to handle pixels of 255,
    // because they will fall in their own bin. Or instead use rounding with an offset. 

The histogram in this calculation is scaled (divided by the pixel count) so that the values make sense as probabilities. You can of course factor the division out of the for loop.

Now you need to calculate the accumulative sum of this:

histogramSum = all zero // length of histogramSum must be one bigger than histogram

for i == 1 .. length(histogram)
    histogramSum[i + 1] = histogramSum[i] + histogram[i]

Now you have to invert this function and this is the tricky part. The best is to not calculate an explicit inverse, but calculate it on the spot, and apply it on the image. The basic idea is to search for the pixel value in the histogramSum (find the closest index below), and then do a linear interpolation between the index and the next index.

foreach pixel in image at coordinates i, j
   hIndex = findIndex(pixel, histogramSum) // You have to write findIndex, it should be simple
   equilisationFactor = (pixel - histogramSum[hIndex])/(histogramSum[hIndex + 1] - histogramSum[hIndex]) * binSize
   // This above is the linear interpolation step. 
   // Notice the technicality that you need to handle:
   // histogramSum[hIndex + 1] may be out of bounds

   equalisedImage[i, j] = pixel * equilisationFactor

---

Edit: without drilling into the maths, I can't be 100% sure, but I think that division by 0 errors are possible. These can occur if one bin is empty, so consecutive sums are equal. So you need special code to handle this case too. The best you can do is take the value for the factor as halfway between hIndex, hIndex + n, where n is the highest value for which histogramSum[hIndex + n] == histogramSum[hIndex].

---

And that should be it, once you have dealt with all the technicalities.

The above algorithm is slow (especially in the findIndex step). You may be able to optimize this with a special lookup datastructure. But only do that when it's working, and only if necessary.

---

One more thing about your Matlab code: the rows and columns are inverted. Because of the symmetry in the algorithm, the result is the same, but it can cause puzzling bugs in other algorithms, and be very confusing if you examine pixel values during debugging. In the pseudocode above I used them the same as you, though.