In the previous blog, we discussed global thresholding and how to find the global threshold using the iterative approach. In this blog, we will discuss Otsu’s method, named after Nobuyuki Otsu, that automatically finds the global threshold. So, let’s discuss this method in detail.
Note: This method assumes that the image histogram is bimodal and a reasonable contrast ratio exists between the background and the region of interest.
In simple terms, Otsu’s method tries to find a threshold value which minimizes the weighted within-class variance. Since Variance is the spread of the distribution about the mean. Thus, minimizing the within-class variance will tend to make the classes compact.
Let’s say we threshold a histogram at a value “t”. This produces two regions – left and right of “t” whose variance is given by σ20 and σ21. Then the weighted within-class variance is given by
where w0(t) and w1(t) are the weights given to each class. Weights are total pixels in a thresholded region (left or right) divided by the total image pixels. Let’s take a simple example to understand how to calculate these.
Suppose we have the following histogram and we want to find the weighted within-class variance corresponding to threshold value 1.
Below are the weights and the variances calculated for left and the right regions obtained after thresholding at value 1.
Similarly, we will iterate over all the possible threshold values, calculate the weighted within-class variance for each of the thresholds. The optimum threshold will be the one with the minimum within-class variance.
Now, let’s see how to do this using python.
|
1 2 3 4 5 6 7 8 9 10 11 12 13 |
import cv2 import matplotlib.pyplot as plt import numpy as np # Create a sample image np.random.seed(7) img = np.random.normal(40,10,size=(500,500)).astype('uint8') img[img>100]=40 img[100:400,100:400] = np.random.normal(150,20,size=(300,300)).astype('uint8') # plot the histogram hist = plt.hist(img.ravel(),256,[0,256]) plt.show() |
The image histogram is shown below
Now, let’s calculate the within-class variance using the steps which we discussed earlier.
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 |
# Set minimum value to infinity final_min = np.inf # total pixels in an image total = np.sum(hist[0]) for i in range(256): # Split regions based on threshold left, right = np.hsplit(hist[0],[i]) # Splt intensity values based on threshold left_bins, right_bins = np.hsplit(hist[1],[i]) # Only perform thresholding if neither side empty if np.sum(left) !=0 and np.sum(right) !=0: # Calculate weights on left and right sides w_0 = np.sum(left)/total w_1 = np.sum(right)/total # Calculate the mean for both sides mean_0 = np.dot(left,left_bins)/np.sum(left) mean_1 = np.dot(right,right_bins[:-1])/np.sum(right) # right_bins[:-1] because matplotlib has uses 1 bin extra # Calculate variance of both sides var_0 = np.dot(((left_bins-mean_0)**2),left)/np.sum(left) var_1 = np.dot(((right_bins[:-1]-mean_1)**2),right)/np.sum(right) # Calculate final within class variance final = w_0*var_0 + w_1*var_1 # if variance minimum, update it if final<final_min: final_min = final thresh = i print(thresh) # 95 |
The gif below shows how the within-class variance (blue dots) varies with the threshold value for the above histogram. The optimum threshold value is the one where the within-class variance is minimum.
OpenCV also provides a builtin function to calculate the threshold using this method.
OpenCV
You just need to pass an extra flag, cv2.THRESH_OTSU in the cv2.threshold() function which we discussed in the previous blog. The optimum threshold value will be returned by this along with the thresholded image. Let’s see how to use this.
|
1 2 |
gray = cv2.imread('kang.jpg',0) retval, thresh = cv2.threshold(gray,0,255,cv2.THRESH_BINARY+cv2.THRESH_OTSU) |
A Faster Approach
We all know that minimizing within-class variance is equivalent to maximizing between-class variance. This maximization operation can be implemented recursively and is faster than the earlier method. The expression for between-class variance is given by
Below are the steps to calculate recursively between-class variance.
- Calculate the histogram of the image.
- Set up weights and means corresponding to the “0” threshold value.
- Loop through all the threshold values
- Update the weights and the mean
- Calculate the between-class variance
- The optimum threshold will be the one with the max variance.
Below is the code in Python that implements the above steps.
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 |
# Calculate the histogram hist = plt.hist(img1.ravel(),256,[0,256]) # Total pixels in the image total = np.sum(hist[0]) # calculate the initial weights and the means left, right = np.hsplit(hist[0],[0]) left_bins, right_bins = np.hsplit(hist[1],[0]) # left weights w_0 = 0.0 # Right weights w_1 = np.sum(right)/total # Left mean mean_0 = 0.0 weighted_sum_0 = 0.0 # Right mean weighted_sum_1 = np.dot(right,right_bins[:-1]) mean_1 = weighted_sum_1/np.sum(right) def recursive_otsu1(hist, w_0=w_0, w_1=w_1, weighted_sum_0=weighted_sum_0, weighted_sum_1=weighted_sum_1, thres=1, fn_max=-np.inf, thresh=0, total=total): if thres<=255: # To pass the division by zero warning if np.sum(hist[0][:thres+1]) !=0 and np.sum(hist[0][thres+1:]) !=0: # Update the weights w_0 += hist[0][thres]/total w_1 -= hist[0][thres]/total # Update the mean weighted_sum_0 += (hist[0][thres]*hist[1][thres]) mean_0 = weighted_sum_0/np.sum(hist[0][:thres+1]) weighted_sum_1 -= (hist[0][thres]*hist[1][thres]) if thres == 255: mean_1 = 0.0 else: mean_1 = weighted_sum_1/np.sum(hist[0][thres+1:]) # Calculate the between-class variance out = w_0*w_1*((mean_0-mean_1)**2) # # if variance maximum, update it if out>fn_max: fn_max = out thresh = thres return recursive_otsu1(hist, w_0=w_0, w_1=w_1, weighted_sum_0=weighted_sum_0, weighted_sum_1=weighted_sum_1, thres=thres+1, fn_max=fn_max, thresh=thresh, total=total) # Stopping condition else: return fn_max,thresh # Check the results var_value, thresh_value = recursive_otsu1(hist, w_0=w_0, w_1=w_1, weighted_sum_0=weighted_sum_0, weighted_sum_1=weighted_sum_1, thres=1, fn_max=-np.inf, thresh=0, total=total) print(var_value, thresh_value) |
This is how you can implement otsu’s method recursively if you consider maximizing between-class variance. Now, let’s discuss what are the limitations of this method.
Limitations
Otsu’s method is only guaranteed to work when
- The histogram should be bimodal.
- Reasonable contrast ratio exists between the background and the roi.
- Uniform lighting conditions are there.
- Image is not affected by noise.
- Size of the background and the roi should be comparable.
There are many modifications done to the original Otsu’s algorithm to address these limitations such as two-dimensional Otsu’s method etc. We will discuss some of these modifications in the following blogs.
In the following blogs, we will also discuss how to counter these limitations so as to get satisfactory results with otsu’s method. Hope you enjoy reading.
If you have any doubt/suggestion please feel free to ask and I will do my best to help or improve myself. Good-bye until next time.