Histogram Equalization

Concept

Generate a processed image which has a uniform histogram. If there is a function which can transform image to have a uniform histogram instead of leaned to certain range of intensity histogram, this function can increase the Dynamic range of image to enhance the contrast.

## Principle

• Find a function f(x) that maps x to y.
• f(x) is a mapping function so it have to meet 2 conditions.
  • f(x) must be one-to-one function. So it must be strictly monotonically increasing.
  • variable x and y has same domain range.

Derivation

p\left(y\right)\left|dy\right|=p\left(x\right)|dx|\left(\because f\left(x\right)\ is\ one-to-function\right)

\int_{0}^{y}{p\left(y\right)dy=\int_{0}^{x}p\left(x\right)dx\left(\because d x,dy>0\right)} \\ \frac{1}{L-1}\times\int_{0}^{y}dy=\int_{0}^{x}p\left(x\right)dx \\ \frac{1}{L-1}y=P\left(x\right)-P\left(0\right)\\ \therefore y=f(x)=\left(L-1\right)p\left(x\right)

Actually p(x) is cumulative distribution frequency so transform function can be written in psudo code like this.

Intensity_{out}=round(\left(L-1\right)cdf(Intensity_{in}))

(Round function is used because we handle discrete image.)

Python Implementation

• Python implementation without opencv

import numpy as np

def histogram(data, bin_cnt, bin_range, normalized=False):
    r
    Compute the histogram of a set of data.
    Simply use uniformed bin size

    Parameters
    ----------
    data (1D ndarray) : Ravelled data array
    bin_cnt(int) : total count of bins
    bin_range (tuple) (min, max) range of data value.

    Returns
    -------
    hist (1D ndarray) : 1D histogram of data.
    bin_edges (1D ndarray) : bin edges 

    Error
    --------
    data, bin_cnt, bin_range must not be None

    See Also 
    --------

    Examples 
    --------

    assert type(data) != None, 'data cannot not be None'
    assert type(bin_cnt) != None, 'bin_cnt cannot be None'
    assert type(bin_range) != None, 'bin_range cannot be None'

    bin_edges = np.linspace(bin_range[0], bin_range[1], bin_cnt)
    hist = np.zeros(bin_cnt, dtype=float)

    for datum in data:
        hist[np.where(bin_edges == datum)] += 1

    if normalized:
        hist /= len(data)

    return hist, bin_edges


def cumulative_histogram(hist):
    r
    Compute the cumulative histogram from given histogram.

    Parameters
    ----------
    hist (1D ndarray) : histogram

    Returns
    -------
    cumulative_hist (1D ndarray) : return None if 'hist' is None

    Error
    --------

    See Also 
    --------

    Examples 
    --------

    if hist is None:
        return None

    bin_cnt = len(hist)
    cumulative_hist = np.zeros(bin_cnt)
    prev_sum = 0

    for i in range(0, bin_cnt):
        prev_sum += hist[i]
        cumulative_hist[i] += prev_sum

    return cumulative_hist


def histogram_equalization(img):
    r
    Return histogram equalized img.

    Parameters
    ----------
    img (NxMx1 ndarray) : 1 channel img

    Returns
    -------
    result_img (NxMx1 ndarray) : Equalized 1 channel img
    Error
    --------

    See Also 
    --------
    histogram(), cumulative_histogram()

    Examples 
    --------

    result_img = np.empty(img.shape)
    ravelled_img = img.ravel()

    hist, bin_edges = histogram(ravelled_img, 256, (0, 255), True)
    c_hist = cumulative_histogram(hist)

    for i in range(0, img.shape[0]):
        for j in range(0, img.shape[1]):
            result_img[i, j] = np.around(255 * c_hist[np.where(bin_edges == img[i, j])])
    return result_img


if __name__ == __main__:
    import PIL.Image
    import cv2
    from matplotlib import pyplot as plt

    img = np.asarray(PIL.Image.open('../test_imgs/lena.jpg').convert('L'), 'uint8')
    ravelled_img = img.ravel()

    hist, bin_edges = histogram(ravelled_img, 256, (0, 255), True)
    c_hist = cumulative_histogram(hist)

    equalized_img = histogram_equalization(img)
    cv2.imwrite(../test_imgs/eqj_lena.jpg, equalized_img)
    e_hist, _ = histogram(equalized_img.ravel(), 256, (0, 255), True)
    e_c_hist = cumulative_histogram(e_hist)

    plt.subplot(2, 2, 1)
    plt.xlabel(Width)
    plt.ylabel(Height)
    plt.title(Image)
    plt.imshow(img, cmap='gray')

    plt.subplot(2, 2, 3)
    plt.xlabel(Width)
    plt.ylabel(Height)
    plt.title(Equalized Image)
    plt.imshow(equalized_img, cmap='gray')

    plt.subplot(4, 2, 2)
    plt.xlabel(Intensity)
    plt.ylabel(Frequency)
    plt.title(Histogram)
    plt.xticks(np.arange(0, 263, 8))
    plt.bar(bin_edges, hist)

    plt.subplot(4, 2, 4)
    plt.xlabel(Intensity)
    plt.ylabel(Frequency)
    plt.title(Cumulative Histogram)
    plt.xticks(np.arange(0, 263, 8))
    plt.bar(bin_edges, c_hist)

    plt.subplot(4, 2, 6)
    plt.xlabel(Intensity)
    plt.ylabel(Frequency)
    plt.title(Equalized Histogram)
    plt.xticks(np.arange(0, 263, 8))
    plt.bar(bin_edges, e_hist)

    plt.subplot(4, 2, 8)
    plt.xlabel(Intensity)
    plt.ylabel(Frequency)
    plt.title(Equalized Cumulative Histogram)
    plt.xticks(np.arange(0, 263, 8))
    plt.bar(bin_edges, e_c_hist)

    plt.tight_layout()
    plt.show()