Last time we stretched the intensities of a specific range of values.  The method for determining the bounds was somewhat arbitrary and we simply spread the intensities equally.  Stretching uses more of the intensity spectrum and, hence, increases contrast, but it may not be the best approach to preserve information.  Could we be a little more intelligent about how we distribute the intensities?

Uniform Distribution

In an image with perfect, uniformly distributed intensities, the probability of any intensity would be the same as the probability of any other intensity.  This generally isn’t possible on a per intensity basis, however, if we consider ranges of intensities, we can get close.  For instance, a random pixel in a high contrast image would have roughly the same probability of being in the light or dark ranges of intensities.  This is achievable through, what I like to call, mathemagic.

Probability

When we roll a die we have six possibilities: 1,2,3,4,5, or 6.  What are the odds that 1 will come up?  1 in 6 right?  In this case our probability is P(1)= .1667.  What about 2? P(2)=.1667 as well, right?  What is the probability of 1 or 2 coming up?  P(1 or 2) = P(1) + P(2) = .3333.  Note that the chances of any side coming up = P(1)+P(2)+P(3)+P(4)+P(5)+P(6)=1.

Back to our situation, what is the probability of a randomly selected pixel having a specific intensity?   Using the histogram, this is trivial to calculate.  Simply take the total number of pixels at the intensity divided by the total number of pixels.

 

Okay, so what?  This is where things get a little tricky.

Probability Density Function

If you were to divide every value in a histogram by the number of pixels you would end up with the probability for every intensity.  This is a discrete probability density function (PDF).

The PDF for our wintery tree from last time:

Cumulative Distribution Function

If we have a given intensity k, what is the probability that a pixel will have an intensity between 0 and k?  I hope it is fairly apparent that P(intensity≤k)=P(0)+P(1)+…+P(k).  In our case, the possible intensities range from 0 to 255. 

So, for 0≤k≤255 the CDF(k) = P(intensity≤k).

Why does this matter?  Well, if we map the current intensities to the CDF, it will evenly distribute our intensities across the spectrum which is exactly what we want.

The Code

private double[] CalculateCDF(WriteableBitmap grayscale)

{

    int[] histogramData = new int[256];

 

    // here we are just calculating the histogram

    // this should be old hat by now

    for (int pixelIndex = 0;

        pixelIndex < grayscale.Pixels.Length;

        pixelIndex++)

    {

        byte intensity = (byte)(grayscale.Pixels[pixelIndex] & 0xFF) ;

        histogramData[intensity]++;

    }

 

 

    double[] cumulativeDistributionFunction = new double[256];

    //the first value is just the probability of the first intensity

    //notice we divide the histogram value by the total number of pixels

    //to get the probability

    cumulativeDistributionFunction[0] =

        (((double)histogramData[0]) / grayscale.Pixels.Length);

    for (int intensity = 1; intensity < 256; intensity++)

    {

        //now we just add the current probability

        //to the previous CDF

        cumulativeDistributionFunction[intensity] =

            (((double)histogramData[intensity]) / grayscale.Pixels.Length)

            + cumulativeDistributionFunction[intensity – 1];

    }

 

    return cumulativeDistributionFunction;

}

 

public WriteableBitmap Equalize(

    WriteableBitmap grayscale,

    double[] cumulativeDistFunction)

{

    WriteableBitmap equalized = new WriteableBitmap(grayscale.PixelWidth, grayscale.PixelHeight);

    int[] histogramData = new int[256];

    for (int pixelIndex = 0;

        pixelIndex < grayscale.Pixels.Length;

        pixelIndex++)

    {

        byte intensity = (byte)(grayscale.Pixels[pixelIndex] & 0xFF);

        histogramData[intensity]++;

    }

 

    histogramData = new int[256];

    int maxCount = 0;

    for (int pixelIndex = 0;

        pixelIndex < grayscale.Pixels.Length;

        pixelIndex++)

    {

        //this is where the magic happens

        //first we get the old intensity

        //then we map the old intensity to

        //the CDF and scale it up to 255

        byte oldintensity = (byte)(grayscale.Pixels[pixelIndex] & 0xFF);

        byte newIntensity =

            (byte)(cumulativeDistFunction[oldintensity] * 255 + .5);

 

        equalized.Pixels[pixelIndex] =

            255 << 24

            | newIntensity << 16

            | newIntensity << 8

            | newIntensity;

 

        histogramData[newIntensity]++;

               

        if (histogramData[newIntensity] > maxCount)

        {

            maxCount = histogramData[newIntensity];

        }

    }

 

    PlotHistogram(histogramData, maxCount);

    return equalized;

}

 

Notice how the CDF is now a straight line?  This means that we have done a pretty good job of distributing the intensities.  If we were to run this procedure again it would not change.  While these results are similar to our contrast stretched tree notice how the spacing is not equal.  In fact higher occurring intensities are spaced further apart than lower occurring intensities.

Summary

The math involved here is a little more intense than the other image enhancement techniques we performed.  This is one procedure that I did not find well explained on the internet.  If you have any difficulties understanding it get in touch with me @azzlsoft on Twitter or rich@azzlsoft.com.

Download Code

http://cid-88e82fb27d609ced.office.live.com/embedicon.aspx/Blog%20Files/PhoneVision/PhoneVision%2010%20-%20Histogram%20Equalization.zip

Up next: Intro to Spatial Filters