Thresholding Image with Unimodal Histogram

Calculating Threshold of Unimodal histogram¶

There are many approaches of finding perfect point threshold the image, automatically, like, Otsu's method, T-point etc.. Here, I have tried to simplify the process by just finding the difference along the array of histogram, and find the maximum difference point. This would give us the point of maximum difference, which we would consider to be a Thresholding point $$T = max(H_i-H_{i-1})$$

In [1]:

import numpy as np
import cv2
import matplotlib.pyplot as plt

Just for trying this, let us take a Guassian Distribution curve

In [2]:

x = np.arange(250)
y = np.exp(-(x-100)**2/(2*16**2))/np.sqrt(2*np.pi*16**2)
y=y/np.amax(y)
plt.plot(x,y)
plt.title("Example - Distribution")
plt.show()

we use np.diff() to calculate differences. It takes the array of $N$ elements as argument and returns array of $N-1$ elements, that undergo $$D_i = A_{i+1}-A_i$$

In [3]:

dy = np.diff(y)
max_change_point = np.amin(np.nonzero(dy==np.amax(dy))[0])
print("Maximum changing point = ",max_change_point)
plt.plot(x,y)
plt.axvline(x=max_change_point,color='red',linestyle = '--')
plt.title("Indicating the point")
plt.show()

Maximum changing point =  83

In [4]:

def get_histogram(img):
    max_val = np.amax(img)
    k=1
    while k<max_val:
        k=k*2
    histogram = np.zeros(k)
    for i in range(img.shape[0]):
        for j in range(img.shape[1]):
            histogram[img[i][j]] = histogram[img[i][j]]+1
    return histogram

def thresholder(img,t):
    for i in range(img.shape[0]):
        for j in range(img.shape[1]):
            if img[i][j] > t:
                img[i][j] = 255
            else:
                img[i][j] = 0
    return img

In [5]:

def unimodal_threshold(histogram):
    H = np.diff(histogram)
    return np.amin(np.nonzero(H==np.amax(H))[0])

I have created a sample image, that has a Unimodal Distribution, and demonstrated the same logic

In [6]:

img = cv2.imread('unimodal_image.png')
img = cv2.cvtColor(img,cv2.COLOR_BGR2GRAY)

histogram = get_histogram(img.copy())
new_thres = unimodal_threshold(histogram)
img_t = thresholder(img.copy(),new_thres)

plt.subplot(2,2,1)
plt.imshow(img,cmap='gray')
plt.title("Original Image")
plt.yticks([],[])
plt.subplot(2,2,2)
plt.imshow(img_t,cmap='gray')
plt.title("Threshold = "+str(new_thres))
plt.yticks([],[])

plt.subplot(2,1,2)
plt.plot(np.arange(len(histogram)),histogram)
plt.axvline(x=new_thres,color='red',linestyle='--')
plt.show()

If at all we have a very noisy histogram, we can smoothen it by using Moving Average Filter which makes our estimation of threshold more accurate and not altering the original image too.

Why is this required?¶

• Few commonly used methods are Otsu's method, which works best with Bimodal histograms
• Computation is very simple in this case One cannot surely say, which method is required where, everything can be altered and used according to the requirement</br>

Applications:¶

Such methods can be used in images used for Astronomical studies (Galaxies, etc) where, one cannot use Otsu's method as it can crush many details in the image.