Histogram Equalization without using built-in histogram methods in python

Question

I wrote the code below and what I get is the output below. What I suppose to do is write an histogram equalization function(without built in methods) I get no error, however output is not what it should to be. I could not notice any logic mistakes im my code. Although, while writing the loop for calculating cdf and/or mapping I couldn't follow what happens behind it exactly, maybe the problem is there but I am not sure.

def my_float2int(img):
    
    img = np.round(img * 255, 0)
    img = np.minimum(img, 255)
    img = np.maximum(img, 0)
    img = img.astype('uint8')
    
    return img

def equalizeHistogram(img):
    
    img_height = img.shape[0]
    img_width = img.shape[1]
    histogram = np.zeros([256], np.int32) 
    
    
    # calculate histogram 
    for i in range(0, img_height):
        for j in range(0, img_width):
            histogram[img[i, j]] +=1
            
    # calculate pdf of the image
    pdf_img = histogram / histogram.sum()
    
    ### calculate cdf 
    # cdf initialize . 
    cdf = np.zeros([256], np.int32) 

    # For loop for cdf
    for i in range(0, 256):
        for j in range(0, i+1):
            cdf[i] += pdf_img[j]
     
    cdf_eq = np.round(cdf * 255, 0) # mapping, transformation function T(x)
    
    imgEqualized = np.zeros((img_height, img_width))
    
    # for mapping input image to s.
    for i in range(0, img_height):
        for j in range(0, img_width):
            r = img[i, j] # feeding intensity levels of pixels into r. 
            s = cdf_eq[r] # finding value of s by finding r'th position in the cdf_eq list.
            imgEqualized[i, j] = s # mapping s thus creating new output image.
            
    # calculate histogram equalized image here
    # imgEqualized = s # change this
    
    return imgEqualized

# end of function


# 2.2 obtain the histogram equalized images using the above function
img_eq_low = equalizeHistogram(img_low)
img_eq_high = equalizeHistogram(img_high)

img_eq_low = my_float2int(img_eq_low)
img_eq_high = my_float2int(img_eq_high)

# 2.3 calculate the pdf's of the histogram equalized images
hist_img_eq_low = calcHistogram(img_eq_low)
hist_img_eq_high = calcHistogram(img_eq_high)

pdf_eq_low = hist_img_eq_low / hist_img_eq_low.sum()
pdf_eq_high = hist_img_eq_high / hist_img_eq_high.sum()

# 2.4 display the histogram equalized images and their pdf's
plt.figure(figsize=(14,8))
plt.subplot(121), plt.imshow(img_eq_low, cmap = 'gray', vmin=0, vmax=255)
plt.title('Hist. Equalized Low Exposure Image'), plt.xticks([]), plt.yticks([])
plt.subplot(122), plt.imshow(img_eq_high, cmap = 'gray', vmin=0, vmax=255)
plt.title('Hist. Equalized High Exposure Image'), plt.xticks([]), plt.yticks([])
plt.show()
plt.close()

My output:

Expected output: with the built-in methods.

Answer 1

I found two minor bugs, and one efficiency issue:

• Replace cdf = np.zeros([256], np.int32) with cdf = np.zeros([256], float)
  In the loop, you are putting float elements in cdf, so the type should be float instead of int32.
• Replace img = np.round(img * 255, 0) with img = np.round(img, 0) (in my_float2int).
  You are scaling img by 255 twice (the first time is in cdf_eq = np.round(cdf * 255, 0)).

You may compute cdf more efficiently.
Your implementation:

for i in range(0, 256):
    for j in range(0, i+1):
        cdf[i] += pdf_img[j]

Suggested implementation (more efficient way for computing "accumulated sum"):

cdf[0] = pdf_img[0]
for i in range(1, 256):
    cdf[i] = cdf[i-1] + pdf_img[i]

It's not a bug, but a kind of academic issue (regarding complexity).

---

Here is an example for corrected code (uses only img_low):

import numpy as np
import cv2

def my_float2int(img):
    
    # Don't use *255 twice
    # img = np.round(img * 255, 0)
    img = np.round(img, 0)
    img = np.minimum(img, 255)
    img = np.maximum(img, 0)
    img = img.astype('uint8')
    
    return img

def equalizeHistogram(img):
    
    img_height = img.shape[0]
    img_width = img.shape[1]
    histogram = np.zeros([256], np.int32) 
    
    
    # calculate histogram 
    for i in range(0, img_height):
        for j in range(0, img_width):
            histogram[img[i, j]] +=1
            
    # calculate pdf of the image
    pdf_img = histogram / histogram.sum()
    
    ### calculate cdf 
    # cdf initialize . 
    # Why does the type np.int32?
    #cdf = np.zeros([256], np.int32)
    cdf = np.zeros([256], float)

    # For loop for cdf
    for i in range(0, 256):
        for j in range(0, i+1):
            cdf[i] += pdf_img[j]

    # You may implement the "accumulated sum" in a more efficient way:
    cdf = np.zeros(256, float)
    cdf[0] = pdf_img[0]
    for i in range(1, 256):
        cdf[i] = cdf[i-1] + pdf_img[i]
     
    cdf_eq = np.round(cdf * 255, 0) # mapping, transformation function T(x)
    
    imgEqualized = np.zeros((img_height, img_width))
    
    # for mapping input image to s.
    for i in range(0, img_height):
        for j in range(0, img_width):
            r = img[i, j] # feeding intensity levels of pixels into r. 
            s = cdf_eq[r] # finding value of s by finding r'th position in the cdf_eq list.
            imgEqualized[i, j] = s # mapping s thus creating new output image.
            
    # calculate histogram equalized image here
    # imgEqualized = s # change this
    
    return imgEqualized

# end of function


# Read input image as Grayscale
img_low = cv2.imread('img_low.png', cv2.IMREAD_GRAYSCALE)

# 2.2 obtain the histogram equalized images using the above function
img_eq_low = equalizeHistogram(img_low)
img_eq_low = my_float2int(img_eq_low)

# Use cv2.imshow (instead of plt.imshow) just for testing.
cv2.imshow('img_eq_low', img_eq_low)
cv2.waitKey()

---

Result: