Category Archives: Image Processing

Dilation

In the previous blog, we discussed erosion operation. In this blog, we will cover another morphological operation – Dilation which is the dual of erosion. Dual in the sense that dilating the object region is equivalent to eroding the background region and vice versa. So, let’s get started.

Dilation

As clear from the name, this operation dilates or expands the object region. This is just opposite of erosion. In this, we ask the simple question of whether the structuring element hits the object or not? (“hits” here means that at least one of the image pixels underlying the structuring element (SE) should have the same value as that of the corresponding SE). If it hits, that pixel is set to 1 else set to 0. This is all the concept behind the dilation operation. Now, let’s formulate this in terms of the set operation.

In general, the dilation of the binary image A by some SE B is defined as

That is the set of all values of z such that the intersection of B (translated by z and reflected about its origin) and A is non-empty. In other words, we place the SE over the image so that the origin of the SE coincides with the input pixel position and compare the underlying image pixels with the pixels of the corresponding SE. If the SE shares at least one common element with its underlying image pixels, then the central image pixel is set to 1 else 0.

Thus this increases the size of the object. If some holes or pepper noise is present in the object, this results in bridging the gaps or removing the noise similar to what we discussed in the low pass filtering. The extent of thickening is controlled by the shape and size of the SE.

For binary images, this can be simply done by taking the maximum of the neighborhood defined by the SE. Now, let’s see how to do this using OpenCV-Python. OpenCV provides a builtin function for this as shown below.

cv2.dilate(src, kernel[, anchor[, iterations[, borderType[, borderValue]]]])

1	cv2.dilate(src, kernel[, anchor[, iterations[, borderType[, borderValue]]]])

Here, src is the input image with any number of channels( all will be processed independently) and the kernel is the structuring element whose origin is defined by the anchor (default (-1,-1) i.e at the center of the SE). You can create the SE using cv2.getStructuringElement() or simply using numpy. Iterations specify how many times to repeat the dilation process. It is sometimes useful to pad the image to account for the boundary pixels or if the image is of non-regular shape and this can be done using the “borderType” and “borderValue” arguments. Below is an example where we dilate the image with the rectangular SE.

import cv2
import numpy as np

# load the image
img = cv2.imread('kang.png', 0)

# Following are the two methods by which we can construct
# the structuring element

# Method-1: using numpy
kernel_1 = np.ones((3,1), dtype='uint8')

# Method-2: using opencv getStructuringElement() function
kernel_2 = cv2.getStructuringElement(cv2.MORPH_RECT,(1,3))

# Apply the dilation function
dilation = cv2.dilate(b, kernel_2, iterations = 1)

import cv2

import numpy as np

# load the image

img = cv2.imread('kang.png', 0)

# Following are the two methods by which we can construct

# the structuring element

# Method-1: using numpy

kernel_1 = np.ones((3,1), dtype='uint8')

# Method-2: using opencv getStructuringElement() function

kernel_2 = cv2.getStructuringElement(cv2.MORPH_RECT,(1,3))

# Apply the dilation function

dilation = cv2.dilate(b, kernel_2, iterations = 1)

Similar to erosion, this can also be used to remove noise, detect the object boundary, etc. Although neither erosion nor dilation alone is effective in reducing noise. A more efficient approach is erosion followed by dilation or opening operation in general. Most of the morphological algorithms which we will discuss in the next blogs are also based on dilation and erosion. 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.

Erosion

Leave a reply

In the previous blog, we touched on the introduction of morphological operations. In this blog and the next blog, we will discuss two of the most fundamental morphological operations – Erosion and Dilation. All other morphological operations can be defined in terms of these two basic operations. So, let’s get started.

Erosion

As clear from the name, this operation erodes or remove the pixels from the object boundary. In this, we ask the simple question of whether the structuring element fits the object or not? (“Fits” here means that all the image pixels underlying the structuring element (SE) should have the same value as that of the corresponding SE) If the image pixel fits, it is assigned 1, otherwise eroded (assigned 0). Thus if we use a square SE (say of size 3×3), then all the object boundary pixels will be eroded away. Now, let’s understand this in terms of the set operation.

In general, the erosion of the binary image A by some SE B is defined as

That is the set of all values of z such that B translated by z, is the subset of A or is contained in A. In other words, you shift the SE over the image and set the positions, where the SE doesn’t share any common element with the background, to 1 and erode all the remaining positions.

Thus this results in a decrease in the object area. If some holes are present in the object, this operation tends to increase the hole area. For binary images, this can be simply applied by taking the minimum of the neighborhood defined by the SE. Now, let’s see how to do this using OpenCV-Python. OpenCV provides a builtin function for this as shown below.

cv2.erode(src, kernel[, anchor[, iterations[, borderType[, borderValue]]]])

1	cv2.erode(src, kernel[, anchor[, iterations[, borderType[, borderValue]]]])

Here, src is the input image with any number of channels( all will be processed independently) and the kernel is the structuring element whose origin is defined by the anchor (default (-1,-1)). You can create the SE using cv2.getStructuringElement() or simply using numpy. Iterations specify how many times to repeat the erosion process. It is sometimes useful to pad the image to account for the boundary pixels or if the image is of non-regular shape and this can be done using the “borderType” and “borderValue” arguments. Below is an example where we erode the image with the rectangular SE.

import cv2
import numpy as np

# load the image
img = cv2.imread('kang.png', 0)

# Following are the two methods by which we can construct
# the structuring element

# Method-1: using numpy
kernel_1 = np.ones((3,3), dtype='uint8')

# Method-2: using opencv getStructuringElement() function
kernel_2 = cv2.getStructuringElement(cv2.MORPH_RECT,(3,3))

# Apply the erosion function
erosion = cv2.erode(b, kernel_2, iterations = 1)

import cv2

import numpy as np

# load the image

img = cv2.imread('kang.png', 0)

# Following are the two methods by which we can construct

# the structuring element

# Method-1: using numpy

kernel_1 = np.ones((3,3), dtype='uint8')

# Method-2: using opencv getStructuringElement() function

kernel_2 = cv2.getStructuringElement(cv2.MORPH_RECT,(3,3))

# Apply the erosion function

erosion = cv2.erode(b, kernel_2, iterations = 1)

Different structuring elements, whether in the terms of shape or size, will produce different results. Mostly people prefer disc-shaped structuring element. If the size of the SE exceeds the size of the object, then the entire object will be eroded away. Erosion can be useful in removing noise (subjected to some conditions), detecting the object boundaries (subtracting the eroded image from the original one), separating the connected components or structures of a certain shape or size, etc.

In the next blog, we will discuss another morphological operation known as Dilation in greater detail. 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.

Morphological Image Processing

Leave a reply

In the previous blogs, we discussed various thresholding algorithms like otsu, adaptive, BHT, etc. All these resulted in a binary image which in general are distorted by noise, holes, etc. Thus there is a need to process these images so as to remove the imperfections. And sometimes, we also need to extract image features such as boundaries, etc that are useful in the representation of the object of interest. This all is done using Morphological image processing that applies non-linear transformations based on the image shape. Now, let’s discuss why the name morphology?

In general, morphology stands for the study of the form and the structure of the things. Because as a result of conversion to the binary image, we have lost the intensity information. Thus the only information that remains is the spatial location or the structure of the image. That’s why known as morphological image processing.

Morphological image processing was originally developed for binary images but later this was also extended to the grayscale images. The watershed algorithm is an outcome of this generalization.

Let’s understand the concept behind the morphological image processing (MIP) wrt. convolution operation that we studied earlier.

Remember in convolution, we have a filter/window and we move this filter over the image. The output values are then computed by some linear operation between the filter weights and the pixel values. Similarly, in MIP we have a structuring element and we move this over the entire image. The output values are computed by applying non-linear operations like set operations (intersection, union, etc.) between the structuring element and the underlying pixel values. These non-linear operations are known as morphological operations.

So, the above paragraph contains two terms – structuring element and the morphological operations. So, let’s understand what’s a structuring element?

The structuring element is a binary image (consisting of 0’s and 1’s) that is used to probe an image for finding the region of interest. For instance, if we want to detect lines in an image, we create a linear structuring element. The pattern of 1’s and 0’s specifies the shape of the structuring element. Below is an example of elliptical SE.

Mostly the dimensions of the SE are odd with the origin at the center. OpenCV provides a builtin function for creating SE as shown below.

cv2.getStructuringElement(shape, ksize[, anchor]) → retval

1	cv2.getStructuringElement(shape, ksize[, anchor]) → retval

Here, shape refers to the SE shape. This can take one of the following values

cv2.MORPH_RECT – creates a rectangular SE.
cv2.MORPH_ELLIPSE – creates an elliptical SE.
cv2.MORPH_CROSS – cross-shaped SE.

The “ksize” specifies the size of the SE and anchor specifies the origin position (default is (-1,-1) i.e at the center of the SE). Below is an example that creates an elliptical SE.

cv2.getStructuringElement(cv2.MORPH_ELLIPSE,(5,5))

1	cv2.getStructuringElement(cv2.MORPH_ELLIPSE,(5,5))

To obtain the output, morphological operations are performed between the SE and the underlying pixel values. Morphological operations are nothing but basic set operations like union, intersection, etc. For instance, an example of morphological operations can be the set of all values such that at least one of the SE pixel values is equal to the underlying image pixel values. This operation usually leads to an increase in the size of the object and fills the holes if present in the object. Below figure shows this morphological operation.

If you use other SE, the result would be different. So, select the shape of the SE according to your application. In the next blog, we will discuss various morphological operations in detail. 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.

Adaptive Thresholding

Leave a reply

In the previous blog, we discussed how global thresholding can be a tedious task when dealing with images having non-uniform illumination. This is because you need to ensure that while subdividing an image, each sub-image histogram is bimodal. Otherwise, the segmentation task will fail.

In this blog, we will discuss adaptive thresholding that works well for varying conditions like non-uniform illumination, etc. In this, the threshold value is calculated separately for each pixel using some statistics obtained from its neighborhood. This way we will get different thresholds for different image regions and thus tackles the problem of varying illumination.

The whole procedure can be summed up as:

For each pixel in the image
- Calculate the statistics (such as mean, median, etc.) from its neighborhood. This will be the threshold value for that pixel.
- Compare the pixel value with this threshold

Now, let’s discuss the OpenCV function for adaptive thresholding.

cv2.adaptiveThreshold(src, maxValue, adaptiveMethod, thresholdType, blockSize, C)

1	cv2.adaptiveThreshold(src, maxValue, adaptiveMethod, thresholdType, blockSize, C)

src: 8-bit greyscale image
thresholdType: This tells us what value to assign to pixels greater/less than the threshold. Must be either THRESH_BINARY or THRESH_BINARY_INV. (You can read more about it here).
maxValue: This is the value assigned to the pixels after thresholding. This depends on the thresholding type. If the type is cv2.THRESH_BINARY, all the pixels greater than the threshold are assigned this maxValue.
adaptiveMethod: This tells us how the threshold is calculated from the pixel neighborhood. This currently supports two methods:
- cv2.ADAPTIVE_THRESH_MEAN_C: In this, the threshold value is the mean of the neighborhood area.
- cv2.ADAPTIVE_THRESH_GAUSSIAN_C: In this, the threshold value is the weighted sum of the neighborhood area. This uses Gaussian weights computed using getGaussiankernel() method. You can read more about it here.
blockSize: This is the neighborhood size.
C: a constant which is subtracted from the threshold.

As discussed OpenCV only provides mean and weighted mean to serve as the threshold. But don’t limit yourself to these two statistics. Try other statistics like standard deviation, median, etc. by writing your own helper function. Let’s see how to use this.

import cv2
# Load the image
img1 = cv2.imread("D:/downloads/adap1.jpg",0)
# Apply Otsu method
ret, thres = cv2.threshold(img2,127,255,cv2.THRESH_BINARY+cv2.THRESH_OTSU)
# Apply adaptive threshold
th3 = cv2.adaptiveThreshold(img2,255,cv2.ADAPTIVE_THRESH_GAUSSIAN_C,cv2.THRESH_BINARY,5,2)
# Display the result
cv2.imshow('original',img1)
cv2.imshow('otsu', thres)
cv2.imshow('adaptive', th3)
cv2.waitKey(0)

import cv2

# Load the image

img1 = cv2.imread("D:/downloads/adap1.jpg",0)

# Apply Otsu method

ret, thres = cv2.threshold(img2,127,255,cv2.THRESH_BINARY+cv2.THRESH_OTSU)

# Apply adaptive threshold

th3 = cv2.adaptiveThreshold(img2,255,cv2.ADAPTIVE_THRESH_GAUSSIAN_C,cv2.THRESH_BINARY,5,2)

# Display the result

cv2.imshow('original',img1)

cv2.imshow('otsu', thres)

cv2.imshow('adaptive', th3)

cv2.waitKey(0)

See how effective adaptive thresholding is in the case of non-uniform illumination. 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.

Balanced histogram thresholding

Leave a reply

In the previous blogs, we discussed different methods for automatically finding the global threshold for an image. For instance, the iterative method, Otsu’s method, etc. In this blog, we will discuss another very simple approach for automatic thresholding – Balanced histogram thresholding. As clear from the name, this method tries to automatically find the threshold by balancing the image histogram. Let’s understand 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.

Concept

Suppose you have a perfectly balanced histogram i.e. a histogram where the distribution of the background and the roi is the same. If you place such a histogram over the lever, it will be balanced. And the optimum threshold will be at the center of the lever as shown in the figure below

This is the main idea behind the Balanced Histogram Thresholding. This method tries to balance the image histogram and then infer the threshold value from that.

But in real-life situations, we don’t encounter images with such perfectly balanced histograms. So, let’s see how this method balances the unbalanced histograms.

First, it places the histogram over the lever and calculates the center point.
Then this calculates the left side and right side weights from the center point.
Removes weight from the heavier side and adjust the center.
Repeat the above two steps until the starting and the endpoints are equal to the center.

The whole procedure can be summed up in the below gif (taken from Wikipedia)

**Credits: By Power3d – Own work, CC BY-SA 3.0**

Below is the python code for this. Here, i_s, i_e are the starting and the endpoints of the histogram and i_m is the center

def balanced_hist_thresholding(b):
    # Starting point of histogram
    i_s = np.min(np.where(b[0]>0))
    # End point of histogram
    i_e = np.max(np.where(b[0]>0))
    # Center of histogram
    i_m = (i_s + i_e)//2
    # Left side weight
    w_l = np.sum(b[0][0:i_m+1])
    # Right side weight
    w_r = np.sum(b[0][i_m+1:i_e+1])
    # Until starting point not equal to endpoint
    while (i_s != i_e):
        # If right side is heavier
        if (w_r > w_l):
            # Remove the end weight
            w_r -= b[0][i_e]
            i_e -= 1
            # Adjust the center position and recompute the weights
            if ((i_s+i_e)//2) < i_m:
                w_l -= b[0][i_m]
                w_r += b[0][i_m]
                i_m -= 1
        else:
            # If left side is heavier, remove the starting weight
            w_l -= b[0][i_s]
            i_s += 1
            # Adjust the center position and recompute the weights
            if ((i_s+i_e)//2) >= i_m:
                w_l += b[0][i_m+1]
                w_r -= b[0][i_m+1]
                i_m += 1
    return i_m

def balanced_hist_thresholding(b):

# Starting point of histogram

i_s = np.min(np.where(b[0]>0))

# End point of histogram

i_e = np.max(np.where(b[0]>0))

# Center of histogram

i_m = (i_s + i_e)//2

# Left side weight

w_l = np.sum(b[0][0:i_m+1])

# Right side weight

w_r = np.sum(b[0][i_m+1:i_e+1])

# Until starting point not equal to endpoint

while (i_s != i_e):

# If right side is heavier

if (w_r > w_l):

# Remove the end weight

w_r -= b[0][i_e]

i_e -= 1

# Adjust the center position and recompute the weights

if ((i_s+i_e)//2) < i_m:

w_l -= b[0][i_m]

w_r += b[0][i_m]

i_m -= 1

else:

# If left side is heavier, remove the starting weight

w_l -= b[0][i_s]

i_s += 1

# Adjust the center position and recompute the weights

if ((i_s+i_e)//2) >= i_m:

w_l += b[0][i_m+1]

w_r -= b[0][i_m+1]

i_m += 1

return i_m

The above function takes the image histogram as the input and returns the optimum threshold. Let’s take an example to check how this works.

import cv2
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
np.random.seed(7)

# Create a sample image
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
b1 = plt.hist(img.ravel(),256,[0,256])
plt.show()

import cv2

import numpy as np

import matplotlib.pyplot as plt

%matplotlib inline

np.random.seed(7)

# Create a sample image

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

b1 = plt.hist(img.ravel(),256,[0,256])

plt.show()

Below is the histogram of the image constructed.

Now, let’s apply the Balanced Histogram thresholding method to check what threshold value this outputs.

thresh_value = balanced_hist_thresholding(b1)
>>> 87

1 2	thresh_value = balanced_hist_thresholding(b1) >>> 87

87 looks like a reasonable threshold, check the image histogram above. So, that’s all for this time. 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.

Optimum Global Thresholding using Otsu’s Method

1 Reply

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 σ²₀ and σ²₁. Then the weighted within-class variance is given by

where w₀(t) and w₁(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.

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()

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.

# 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

# 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.

gray = cv2.imread('kang.jpg',0)
retval, thresh = cv2.threshold(gray,0,255,cv2.THRESH_BINARY+cv2.THRESH_OTSU)

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
1. Update the weights and the mean
2. 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.

# 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)

# 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.

Improving Global Thresholding

Leave a reply

In the previous blog, we discussed otsu’s method for automatic image thresholding. Then we also discussed the limitations of the otsu’s method. In this blog, we will discuss how to handle these limitations so as to produce satisfactory thresholding results. So, let’s get started.

Case-1: When the noise is present in the image

If the noise is present in the image, then this tends to change the modality of the histogram. The sharp valleys between the peaks of the bimodal histogram start degrading. In that case, the otsu’s method or any other global thresholding method will fail. So, in order to find the global threshold, one should first remove the noise using any smoothing filters like Gaussian, etc. and then apply any automatic thresholding method like otsu, etc.

Case-2: When the object area is small compared to the background area

In this case, the image histogram will be dominated by a large background area. This will increase the probability of any pixel belonging to the background. So, the histogram will no longer exhibit bimodality and thus otsu will result in segmentation error. To prevent this, one should only consider pixels that lie on or near the edges between the objects and the background. Doing so will result in an image histogram with peaks of approximately the same size. Then we can apply any automatic thresholding method like otsu, etc. Below are the steps to implement the above procedure.

Calculate the edge image using any high pass filter like Sobel, Laplacian, etc.
Select any threshold value (T).
Threshold the above edge image to produce a binary mask.
Apply the mask image on the input image using any bitwise operations or any other method.
This results in only those pixels where the mask image was white.
Compute the histogram of only those pixels
Finally, apply any automatic global thresholding method like otsu, etc.

Case-3: When the image is taken under non-uniform illumination conditions

In this case, the histogram no longer remains bimodal and thus we will not be able to segment the image satisfactorily. One of the simplest approaches is to subdivide the image into non-overlapping images/rectangles. The size of these rectangles is chosen such that the illumination is nearly constant in each of these rectangles. Then we will apply any global thresholding technique like otsu for each of these rectangles.

The above procedure only works when the size of the object and the background are comparable in the rectangle. This is quite intuitive as only then we will have a bimodal histogram. Taking care of the background and the object sizes in each rectangle is a tedious task.

So, in the next blog, we will discuss adaptive thresholding that works pretty well for the above conditions. That’s all for this blog. 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.

Global Thresholding

Leave a reply

In the previous blog, we discussed image thresholding and when to use this for image segmentation. We also learned that thresholding can be global or adaptive depending upon how the threshold value is selected.

In this blog, we will discuss

global thresholding
OpenCV function for global thresholding
How to choose threshold value using the iterative algorithm

In global thresholding, each pixel value in the image is compared with a single (global) threshold value. Below is the code for this.

def global_threshold(image, thres_value, val_high, val_low):
    img = image.copy()
    for i in range(image.shape[0]):
        for j in range(image.shape[1]):
            if image[i,j] > thres_value:
                img[i,j] = val_high
            else:
                img[i,j] = val_low
    return img

def global_threshold(image, thres_value, val_high, val_low):

img = image.copy()

for i in range(image.shape[0]):

for j in range(image.shape[1]):

if image[i,j] > thres_value:

img[i,j] = val_high

else:

img[i,j] = val_low

return img

Here, we assign a value of “val_high” to all the pixels greater than the threshold otherwise “val_low“. OpenCV also provides a builtin function for thresholding the image. So, let’s take a look at that function.

OpenCV

cv2.threshold(src, thresh, maxval, type) → retval, dst

1	cv2.threshold(src, thresh, maxval, type) → retval, dst

This function returns the thresholded image(dst) and the threshold value(retval). Its arguments are

src: input greyscale image (8-bit or 32-bit floating point)
thresh: global threshold value
type: Different types that decide “val_high” and “val_low“. In other words, these types decide what value to assign for pixels greater than and less than the threshold. Below figure shows different thresholding types available.
maxval: maximum value to be used with THRESH_BINARY and THRESH_BINARY_INV. Check the below image.

To specify the thresholding type, write “cv2.” as the prefix. For instance, write cv2.THRESH_BINARY if you want to use this type. Let’s take an example

import cv2
# Load an image in the greyscale
img = cv2.imread('D:/downloads/opencv_logo1.PNG',cv2.IMREAD_GRAYSCALE)
# threshold the image
ret, thresh = cv2.threshold(img,75,255,cv2.THRESH_BINARY)

import cv2

# Load an image in the greyscale

img = cv2.imread('D:/downloads/opencv_logo1.PNG',cv2.IMREAD_GRAYSCALE)

# threshold the image

ret, thresh = cv2.threshold(img,75,255,cv2.THRESH_BINARY)

Similarly, you can apply other thresholding types to check how they work. Till now we discussed how to threshold an image using a global threshold value. But we didn’t discuss how to get this threshold value. So, in the next section, let’s discuss this.

How to choose the threshold value?

As already discussed, that global thresholding is a suitable approach only when intensity distributions of the background and the ROI are sufficiently distinct. In other words, there is a clear valley between the peaks of the histogram. We can easily select the threshold value in that situation. But what if we have a number of images. In that case, we don’t manually want to first check the image histogram and then deciding the threshold value. We want something that can automatically estimate the threshold value for each image. Below is the algorithm that can be used for this purpose.

Below is the code for the above algorithm.

def thres_finder(img, thres=20,delta_T=1.0):
    
    # Step-2: Divide the images in two parts
    x_low, y_low = np.where(img<=thres)
    x_high, y_high = np.where(img>thres)
    
    # Step-3: Find the mean of two parts
    mean_low = np.mean(img[x_low,y_low])
    mean_high = np.mean(img[x_high,y_high])
    
    # Step-4: Calculate the new threshold
    new_thres = (mean_low + mean_high)/2
    
    # Step-5: Stopping criteria, otherwise iterate
    if abs(new_thres-thres)< delta_T:
        return new_thres
    else:
        return thres_finder(img, thres=new_thres,delta_T=1.0)

def thres_finder(img, thres=20,delta_T=1.0):

# Step-2: Divide the images in two parts

x_low, y_low = np.where(img<=thres)

x_high, y_high = np.where(img>thres)

# Step-3: Find the mean of two parts

mean_low = np.mean(img[x_low,y_low])

mean_high = np.mean(img[x_high,y_high])

# Step-4: Calculate the new threshold

new_thres = (mean_low + mean_high)/2

# Step-5: Stopping criteria, otherwise iterate

if abs(new_thres-thres)< delta_T:

return new_thres

else:

return thres_finder(img, thres=new_thres,delta_T=1.0)

Now, let’s take an example to check how’s this working.

import cv2
# Load an image in the greyscale
img = cv2.imread('D:/downloads/opencv_logo1.PNG',cv2.IMREAD_GRAYSCALE)
# apply threshold finder
vv1 = thres_finder1(img, thres=30,delta_T=1.0)
# threshold the image
ret, thresh = cv2.threshold(img,vv1,255,cv2.THRESH_BINARY)
# Display the image side by side
out = cv2.hconcat([img,thresh])
cv2.imshow('threshold',out)
cv2.waitKey(0)

import cv2

# Load an image in the greyscale

img = cv2.imread('D:/downloads/opencv_logo1.PNG',cv2.IMREAD_GRAYSCALE)

# apply threshold finder

vv1 = thres_finder1(img, thres=30,delta_T=1.0)

# threshold the image

ret, thresh = cv2.threshold(img,vv1,255,cv2.THRESH_BINARY)

# Display the image side by side

out = cv2.hconcat([img,thresh])

cv2.imshow('threshold',out)

cv2.waitKey(0)

That’s all for this blog. In the next blog, we will discuss how to perform optimum global thresholding using 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.

Image Thresholding

Leave a reply

Image Segmentation is the process of subdividing an image into its constituent regions or objects. In many computer vision applications, image segmentation is very useful to detect the region of interest. For instance, in medical imaging where we have to locate tumors, or in object detection like self-driving cars have to detect pedestrians, traffic signals, etc or for video surveillance, etc. There are a number of methods available to perform image segmentation. For instance, thresholding, clustering methods, graph partitioning methods, and convolutional methods to mention a few.

In this blog, we will discuss Image Thresholding which is one of the simplest methods for image segmentation. In this, we partition the images directly into regions based on the intensity values. So, let’s discuss image thresholding in greater detail.

Concept

If the pixel value is greater than a threshold value, it is assigned one value (maybe white), else it is assigned another value (maybe black).

In other words, if f(x,y) is the input image then the segmented image g(x,y) is given by

If the threshold value T remains constant over the entire image, then this is known as global thresholding. When the value of T changes over the entire image or depends upon the pixel neighborhood, then this is known as adaptive thresholding. We will cover both these types in greater detail in the following blogs.

Applicability Condition

Thresholding is only guaranteed to work when a good contrast ratio between the region of interest and the background exists. Otherwise, the thresholding will not be able to fully detect the region of interest. Let’s understand this by an example.

Suppose we have two images from which we want to segment the square region (our region of interest) from the background.

Let’s plot the histogram of these two images.

Clearly as expected for “A“, the histogram is showing two peaks corresponding to the square and the background. The separation between the peaks shows that the background and ROI have a good contrast ratio. By choosing a threshold value between the peaks, we will be able to segment out the ROI. While for “B”, the intensity distribution of the ROI and the background is not that distinct. Thus we may not be able to fully segment the ROI.

Thresholded images are shown below (How to choose a threshold value will be discussed in the next blog).

So, always plot the image histogram to check the contrast ratio between the background and the ROI. Only if the contrast ratio is good, choose the thresholding method for image segmentation. Otherwise, look for other methods.

In the next blog, we will discuss global thresholding and how to choose the threshold value using the iterative 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.

Laplacian of Gaussian (LoG)

2 Replies

In the previous blog, we discuss various first-order derivative filters. In this blog, we will discuss the Laplacian of Gaussian (LoG), a second-order derivative filter. So, let’s get started

Mathematically, the Laplacian is defined as

Unlike first-order filters that detect the edges based on local maxima or minima, Laplacian detects the edges at zero crossings i.e. where the value changes from negative to positive and vice-versa.

Let’s obtain kernels for Laplacian similar to how we obtained kernels using finite difference approximations for the first-order derivative.

Adding these two kernels together we obtain the Laplacian kernel as shown below

This is called a negative Laplacian because the central peak is negative. Other variants of Laplacian can be obtained by weighing the pixels in the diagonal directions also. Make sure that the sum of all kernel elements is zero so that the filter gives zero response in the homogeneous regions.

Let’s now discuss some properties of the Laplacian

Unlike first-order that requires two masks for finding edges, Laplacian uses 1 mask but the edge orientation information is lost in Laplacian.
Laplacian gives better edge localization as compared to first-order.
Unlike first-order, Laplacian is an isotropic filter i.e. it produces a uniform edge magnitude for all directions.
Similar to first-order, Laplacian is also very sensitive to noise

To reduce the noise effect, image is first smoothed with a Gaussian filter and then we find the zero crossings using Laplacian. This two-step process is called the Laplacian of Gaussian (LoG) operation.

But this can also be performed in one step. Instead of first smoothing an image with a Gaussian kernel and then taking its Laplace, we can obtain the Laplacian of the Gaussian kernel and then convolve it with the image. This is shown below where f is the image and g is the Gaussian kernel.

Now, let’s see how to obtain LoG kernel. Mathematically, LoG can be written as

The LoG kernel weights can be sampled from the above equation for a given standard deviation, just as we did in Gaussian Blurring. Just convolve the kernel with the image to obtain the desired result, as easy as that.

Select the size of the Gaussian kernel carefully. If LoG is used with small Gaussian kernel, the result can be noisy. If you use a large Gaussian kernel, you may get poor edge localization.

Now, let’s see how to do this using OpenCV-Python

OpenCV-Python

OpenCV provides a builtin function that calculates the Laplacian of an image. You can find it here. Below is the basic syntax of what this function looks like

cv2.Laplacian(src, ddepth[, ksize[, scale[, delta[, borderType]]]])
# src - input image
# ddepth - Desired depth of the destination image.
# ksize - kernel size

cv2.Laplacian(src, ddepth[, ksize[, scale[, delta[, borderType]]]])

# src - input image

# ddepth - Desired depth of the destination image.

# ksize - kernel size

Steps for LoG:

Apply LoG on the image. This can be done in two ways:
- First, apply Gaussian and then Laplacian or
- Convolve the image with LoG kernel directly
Find the zero crossings in the image
Threshold the zero crossings to extract only the strong edges.

Let’s understand each step through code

import cv2
import numpy as np

# Load the image in greyscale
img = cv2.imread('D:/downloads/clouds.jpg',0)

# Apply Gaussian Blur
blur = cv2.GaussianBlur(img,(3,3),0)

# Apply Laplacian operator in some higher datatype
laplacian = cv2.Laplacian(blur,cv2.CV_64F)

import cv2

import numpy as np

# Load the image in greyscale

img = cv2.imread('D:/downloads/clouds.jpg',0)

# Apply Gaussian Blur

blur = cv2.GaussianBlur(img,(3,3),0)

# Apply Laplacian operator in some higher datatype

laplacian = cv2.Laplacian(blur,cv2.CV_64F)

Since zero crossings is a change from negative to positive and vice-versa, so an approximate way is to clip the negative values to find the zero crossings.

# But this tends to localize the edge towards the brighter side.
laplacian1 = laplacian/laplacian.max()

cv2.imshow('a7',laplacian1)
cv2.waitKey(0)

# But this tends to localize the edge towards the brighter side.

laplacian1 = laplacian/laplacian.max()

cv2.imshow('a7',laplacian1)

cv2.waitKey(0)

Another way is to check each pixel for zero crossing as shown below

def Zero_crossing(image):
    z_c_image = np.zeros(image.shape)
    
    # For each pixel, count the number of positive
    # and negative pixels in the neighborhood
    
    for i in range(1, image.shape[0] - 1):
        for j in range(1, image.shape[1] - 1):
            negative_count = 0
            positive_count = 0
            neighbour = [image[i+1, j-1],image[i+1, j],image[i+1, j+1],image[i, j-1],image[i, j+1],image[i-1, j-1],image[i-1, j],image[i-1, j+1]]
            d = max(neighbour)
            e = min(neighbour)
            for h in neighbour:
                if h>0:
                    positive_count += 1
                elif h<0:
                    negative_count += 1


            # If both negative and positive values exist in 
            # the pixel neighborhood, then that pixel is a 
            # potential zero crossing
            
            z_c = ((negative_count > 0) and (positive_count > 0))
            
            # Change the pixel value with the maximum neighborhood
            # difference with the pixel

            if z_c:
                if image[i,j]>0:
                    z_c_image[i, j] = image[i,j] + np.abs(e)
                elif image[i,j]<0:
                    z_c_image[i, j] = np.abs(image[i,j]) + d
                
    # Normalize and change datatype to 'uint8' (optional)
    z_c_norm = z_c_image/z_c_image.max()*255
    z_c_image = np.uint8(z_c_norm)

    return z_c_image

def Zero_crossing(image):

z_c_image = np.zeros(image.shape)

# For each pixel, count the number of positive

# and negative pixels in the neighborhood

for i in range(1, image.shape[0] - 1):

for j in range(1, image.shape[1] - 1):

negative_count = 0

positive_count = 0

neighbour = [image[i+1, j-1],image[i+1, j],image[i+1, j+1],image[i, j-1],image[i, j+1],image[i-1, j-1],image[i-1, j],image[i-1, j+1]]

d = max(neighbour)

e = min(neighbour)

for h in neighbour:

if h>0:

positive_count += 1

elif h<0:

negative_count += 1

# If both negative and positive values exist in

# the pixel neighborhood, then that pixel is a

# potential zero crossing

z_c = ((negative_count > 0) and (positive_count > 0))

# Change the pixel value with the maximum neighborhood

# difference with the pixel

if z_c:

if image[i,j]>0:

z_c_image[i, j] = image[i,j] + np.abs(e)

elif image[i,j]<0:

z_c_image[i, j] = np.abs(image[i,j]) + d

# Normalize and change datatype to 'uint8' (optional)

z_c_norm = z_c_image/z_c_image.max()*255

z_c_image = np.uint8(z_c_norm)

return z_c_image

Depending upon the image you may need to apply thresholding and median blurring to suppress the noise.

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.

TheAILearner

Mastering Artificial Intelligence

Category Archives: Image Processing

Dilation

Dilation

Erosion

Erosion

Morphological Image Processing

Adaptive Thresholding

Balanced histogram thresholding

Concept

Optimum Global Thresholding using Otsu’s Method

OpenCV

A Faster Approach

Limitations

Improving Global Thresholding

Case-1: When the noise is present in the image

Case-2: When the object area is small compared to the background area

Case-3: When the image is taken under non-uniform illumination conditions

Global Thresholding

OpenCV

How to choose the threshold value?

Image Thresholding

Concept

Applicability Condition

Laplacian of Gaussian (LoG)

OpenCV-Python

Steps for LoG: