3.2.2. Histogram Equalization¶

Histogram equalization is a point operator such that the histogram of the resultant image is constant. Histogram equalization is often used to correct for varying illumination conditions.

In the first column from top to botton: the original image, its histogram and its cumulative histogram. In the second column the histogram equalized version of the image with its histogram and cumulative histogram.

3.2.2.1. Definition¶

Let \(f:\set D\rightarrow[0,1]\in\setR\) be a scalar image with histogram \(h_f\). Histogram equalization constructs an image point operator \(\Psi\) such that the resultant image \(\Psi f\) has a constant histogram, meaning that all scalar values in the resultant image are equally probable.

Constructing the transform \(\Psi\) is easy, assuming that it is a point transform, i.e. \((\Psi f)(\v x) = \psi(f(\v x))\) and that the ordering of scalar values is to be preserved, i.e.:

\[v_1 \leq v_2 \Rightarrow \psi(v_1) \leq \psi(v_2)\]

Because the ordering is preserved the p-th percentile of image \(f\) should be equal to the p-th percentile of \(g = \Psi f\). In terms of the cumulative histograms:

\[H_f(v) = H_{g}(\psi(v))\]

Because the histogram of \(g\) is flat (constant) and the image values run from 0 to 1, we have that \(H_{g}(v)=v\), so:

\[\psi(v) = H_f(v)\]

3.2.2.2. Discretization¶

Now consider the image \(f:\set D\rightarrow[0,255]\in\setZ\). The histogram of \(f\) has at most \(256\) (non empty) bins. The value \(h_f(v)\) for \(v=0,\ldots255\) is the height of the bar in the bar plot of the histogram. Histogram equalization then amounts to move the bars over the entire range while preserving the order. Bars can be put on top of eachother but the left right relation (order) can not be changed. And of course one bar cannot be split up in two bars. The finite number of different values in the image range thus accounts for the fact that a flat histogram of the result image is not feasable.

3.2.2.3. Practical Use¶

Histogram equalization is an important image processing operation in practice for the following reason. Consider two images \(f_1\) and \(f_2\) of the same object but taken under two different illumination conditions (say one image taken on a bright and sunny day and the other image taken on a cloudy day). The difference between these images can be approximated with an order preserving mapping \(\gamma\), i.e. \(f_2 = \gamma( f_1)\).

It is not hard to prove that the histogram equalized versions of both images are equal: \(\Psi_1 f_1 = \Psi_2 f_2\). Here \(\Psi_1\) and \(\Psi_2\) are the histogram equalizing mappings defined above (these mappings are image dependent, that is why we write the subscript to denote the image).

Histogram equalization thus serves as a preprocessing step in vision systems to compensate for the effect of an unknown change in illumination (grey value scaling). For instance in a face recognition system where images of human faces have to be compared with other images of the same face (and most probably taken under different illumination conditions).

3.2.2.4. Code¶

Using the Numpy function histogram to make histograms and the interp function for linear table interpolation, the code for histogram equalization is simple:

def histogramEqualization(f, bins=100):
    his, be = histogram(f, range=(0,1), bins=bins)
    his = his.astype(float)/sum(his)
    return interp(f, be, hstack((zeros((1)), cumsum(his))))

Here we assume that the original image \(f\) has \([0,1]\in\setR\) as range. We leave it as an exercise to the reader to adapt the code to the case the image has range from \(0\) to \(M\).

3.2.2.5. Exercises¶

Equalization. Let \(f:\set D\rightarrow[0,M]\in\setR\). Give the expression for the function \(\psi\) in this case.
Code. In practice scalar images are often 8, 10 or 12 bit per pixel. In such cases it is unnescessary to first make the image into a floating point image (with 32 or even 64 bits per pixel), you would like a histogram equalization algorithm to work on the integer representation instead of a real valued representation. Write a program to do just that.

Hint: an often used way to do this is to work with a look-up table. For each possible scalar value v you calculate the histogram equalized value (and you may use floating point values here) and store the result in an array, say HE. In Numpy the actual table lookup then is as simple as HE[image] where image is the original image array.

Please note that the imread function from matplotlib.pylab always returns a floating point representation of an image even it was stored in the file in an 8 bit per pixel format. You can use the imread function from scipy.ndimage instead. The image ‘cameraman.png’ from the standard images is an 8 bit per pixel image.
Practical Use. Make several pictures of the same scene/object from different points of view. Most likely the overal intensity scaling in the images is not equal due to the automatic lighting correction in your camera. Use histogram equalization to correct for this.
Interpolation. The function interp is actually a simple interpolator for a one dimensional function. Let xs be the array of sample x-values and F the array of sample function values of a function in these sample points. Then interp(x,xs,F) calculates the interpolated function value for given x using a simple linear interpolation scheme. Let the number of samples in xs (and thus in F) be equal to \(N\) what is the order for an efficient algorithm to implement the interp function? Note: you may only assume that the sample x values in xs are ordered, i.e. that all(diff(xs))>0 is True.