Implementation of Canny Edge Detection using Python

Task Description: Understand and implement Canny Edge Detection Algorithm from scratch using Python

Note: This is a note to self while learnign Computer Vision basics.

The edge of the image contains rich information. Canny Edge Detector is one of the most effective edge detectors whose input is the output of Sobel. In a nutshell, it finds the angle and applies the threshold to the output of Sobel.

Apply Sobel
Find gradient direction
Apply Non-max suppression
Apply Double Thresholding and Hysteresis

Sobel Filter

$G_x$ = $\begin{bmatrix} -1 & 0 & +1\\ -2 & 0 & +2\\ -1 & 0 & +1 \end{bmatrix}$ $* I$

$G_y$ = $\begin{bmatrix} -1 & -2 & -1\\ 0 & 0 & 0\\ +1 & +2 & +1 \end{bmatrix}$ $* I$

Magnitude:

$G = \sqrt{G_x^2 + G_y^2}$

Non-max suppression

The Non-max suppression algorithm selects high scoring and deletes close-by less confident neighbours. So, ultimately it creates thin edges from the thick edges created by Sobel Filter.

In order to implement this we have approaches like rounding the angle and interpolation.

Rounding

Let q be the pixel that we are working on and $\theta$ be the angle. Let $p_1$ , $p_2, p_3, p_4, p_5, p_6, p_7$ and $p_8$ be the neighbours of $q$ in right, right bottom, bottom, left bottom, left, left top, top and top right respectively. Then we can determine q by comparing the value with the neighbour.

For example let $\theta = 20\degree$ , q is compared with $p_1$ and $p_5$

Then q is suppressed if the value at $q$ is less than values at $p_1$ and $p_5$ .

Interpolation

Let us assume the following matrix. Here also we look at the gradient and determine in between which angles does the gradient lie.

Let us suppose we are looking at a pixel value at 1,1. Let $\theta$ be between 0 and 45. So we get the following image. So here we have to interpolate between two points $y_0$ and $y_1$ . Since Y axis is parallel to this, we find the sin component of the difference in the pixel value. After this adding $y_0$ offset to the obtained value gives the pixel value ( $Y_1$ ) that we are looking for. Similarly we find $Y_2$ on the other side (in between 180 and -135). Finally compare these values at interpolant with the value at 1,1. If value at 1,1 is smaller than both $Y_1$ and $Y_2$ , then it is suppressed (set to zero).

This process is carried out in every pixel of the image.

Double thresholding

This is simple using double threshold values and finding strong, weak and non-relevant pixels.

Strong pixels are pixels that have an intensity so high that we are sure they contribute to the final edge.
Weak pixels are pixels that have an intensity value that is not enough to be considered as strong ones, but yet not small enough to be considered as non-relevant for the edge detection.
Other pixels are considered as non-relevant for the edge.

So we can see that

High threshold is used to identify the strong pixels (intensity higher than the high threshold)
Low threshold is used to identify the non-relevant pixels (intensity lower than the low threshold)
All pixels having intensity between both thresholds are flagged as weak and the Hysteresis mechanism (next step) will help us identify the ones that could be considered as strong and the ones that are considered as non-relevant.

Edge tracking by Hysteresis

Finally we use hysteresis for edge tracking. Based on the threshold results, the hysteresis transforms weak pixels to strong ones, if and only if at least one of the pixels around the one being processed is strong.

Code

# Canny Edge Detection From Scratch
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.image as mpimg
from scipy import ndimage

def rgb2gray(rgb):
    return np.dot(rgb[...,:3], [0.299, 0.587, 0.144])
img = mpimg.imread('res/hs.jpg')
gray = rgb2gray(img)
plt.imshow(gray, cmap = plt.get_cmap('gray'))
plt.show()

camera_blurred = ndimage.gaussian_filter(gray, sigma=1.4)
print(camera_blurred.shape)
plt.imshow(camera_blurred, cmap = plt.get_cmap('gray'))
plt.show()

def sobelFilter(img, filterDirection):
    '''
    SobelFilter(image, gradient)

    colvolves the image with Sobel filters and returns ndimage.
    '''
    if(filterDirection == 'x'):
        Gx = np.array([[-1, 0, +1], [-2, 0, +2], [-1, 0, +1]])
        convolved = ndimage.convolve(img, Gx)
    elif(filterDirection == 'y'):
        Gy = np.array([[-1, -2, -1], [0, 0, 0], [+1, +2, +1]])
        convolved = ndimage.convolve(img, Gy)
    return convolved

# Normalize the pixel to <= 1
def Normalize(img):
    img = img/np.max(img)
    return img

# Apply Sobel filter in X - direction (Gx)
vertical_edge = sobelFilter(camera_blurred, 'x')
vertical_edge = Normalize(vertical_edge)
plt.imshow(vertical_edge, cmap = plt.get_cmap('gray'))
plt.show()

# Apply Sobel filter in Y - direction (Gy)
horizontal_edge = sobelFilter(camera_blurred, 'y')
horizontal_edge = Normalize(horizontal_edge)
plt.imshow(horizontal_edge, cmap = plt.get_cmap('gray'))
plt.show()

# Calculate the magnitude G = sqrt(gx^2 + gy^2)
magnitude = np.hypot(horizontal_edge, vertical_edge)
magnitude = Normalize(magnitude)
plt.imshow(magnitude, cmap= plt.get_cmap('gray'))
plt.show()

# Calculate angle
# theta = arctan(Gy / Gx)
# arctan2 return arctan(Gy / Gx) randian in between -pi to pi and is converted to degree
theta = np.degrees(np.arctan2( horizontal_edge, vertical_edge ))

# Perform Non-Maximum supression without interpolation
def NonMaxSupWithoutInterpol(g_mag, grad):
    '''
    NonMaxSupWithoutInterpol(g_mag, grad)

    supreses the pixel that is not grater than it's neighbour.
    neighbours are selected based on the gradient angle.

    returns numpy array
    '''
    NMS = np.zeros(g_mag.shape)
    for i in range(1, g_mag.shape[0] - 1):
        for j in range(1, g_mag.shape[1] - 1):
            # Left, Right
            if( (grad[i, j] >= -22.5 and grad[i, j] <= 22.5) or (grad[i, j] >= 157.5 and grad[i, j] <= -157.5) ):
                if( g_mag[i, j] > g_mag[i, j+1] and g_mag[i, j] > g_mag[i, j-1] ):
                    NMS[i, j] = g_mag[i, j]
                else:
                    # Set black
                    NMS[i, j] = 0
            # Top right, Bottom left
            elif( (grad[i, j] >= 22.5 and grad[i, j] <= 67.5) or (grad[i, j] >= -157.5 and grad[i, j] <= -112.5) ):
                if( g_mag[i, j] > g_mag[i-1, j+1] and g_mag[i, j] > g_mag[i+1, j-1] ):
                    NMS[i, j] = g_mag[i, j]
                else:
                    NMS[i, j] = 0
            # Top, bottom
            elif( (grad[i, j] >= 67.5 and grad[i, j] <= 112.5) or (grad[i, j] >= -112.5 and grad[i, j] <= -67.5) ):
                if( g_mag[i, j] > g_mag[i-1, j] and g_mag[i, j] > g_mag[i+1, j] ):
                    NMS[i, j] = g_mag[i, j]
                else:
                    NMS[i, j] = 0
            # Top left, bottom right
            elif( (grad[i, j] >= 112.5 and grad[i, j] <= 157.5) or (grad[i, j] >= -67.5 and grad[i, j] <= -22.5) ):
                if( g_mag[i, j] > g_mag[i-1, j-1] and g_mag[i, j] > g_mag[i+1, j+1] ):
                    NMS[i, j] = g_mag[i, j]
                else:
                    NMS[i, j] = 0
    return NMS

NMS = NonMaxSupWithoutInterpol(magnitude, theta)
NMS = Normalize(NMS)
plt.imshow(NMS, cmap = plt.get_cmap('gray'))
plt.show

# NMS with interpolation
def NonMaxSuppressionWithInterpol(g_mag, grad, g_x, g_y):
    '''
    NonMaxSuppressionWithInterpol(g_mag, grad, g_x, g_y)

    supreses the pixel that is not grater than it's neighbour.
    neighbour values are selected with interpolation.

    returns numpy array
    '''
    NMS = np.zeros(g_mag.shape)

    for i in range(1, g_mag.shape[0] - 1):
        for j in range(1, g_mag.shape[1] - 1):
            # 0 <= gradient < 45 and  180 <= gradient < -135
            if( (grad[i, j] >= 0 and grad[i, j] <= 45) or (grad[i, j] >= 180 and grad[i, j] < -135) ):
                # y = [y0, y1]
                # y0 = y[a,b], if it is at top or bottom or left or right
                # y1 = y[a,b], otherwise
                y_bot = np.array( [ g_mag[i, j+1], g_mag[i+1, j+1] ] )
                y_top = np.array( [ g_mag[i, j-1], g_mag[i-1, j-1] ] )
                # As we are interpolating in Y direction
                # sin (theta) = y / magnitude
                x_est = np.absolute( g_y[i, j] / g_mag[i, j] )

                #                   (magnitude along Y)sin(theta) + offset
                if(g_mag[i,j] >= ( (y_bot[1]-y_bot[0])*x_est + y_bot[0] ) and g_mag[i, j] >= ( (y_top[1]-y_top[0])*x_est + y_top[0] ) ):
                    NMS[i,j] = g_mag[i,j]
                else:
                    NMS[i,j] = 0
            # 45 <= gradient < 90 and  -135 <= gradient < -90
            if( (grad[i, j] >= 45 and grad[i, j] <= 90) or (grad[i, j] >= -135 and grad[i, j] < -90) ):
                y_bot = np.array( [ g_mag[i+1, j], g_mag[i+1, j+1] ] )
                y_top = np.array( [ g_mag[i-1, j], g_mag[i-1, j-1] ] )
                # As we are interpolating in X direction
                # cos(theta) = x / magnitude
                x_est = np.absolute( g_x[i, j] / g_mag[i, j] )

                #                   (magnitude along X)cos(theta) + offset
                if(g_mag[i,j] >= ( (y_bot[1]-y_bot[0])*x_est + y_bot[0] ) and g_mag[i, j] >= ( (y_top[1]-y_top[0])*x_est + y_top[0] ) ):
                    NMS[i,j] = g_mag[i,j]
                else:
                    NMS[i,j] = 0
             # 90 <= gradient < 135 and  -90 <= gradient < -45
            if( (grad[i, j] >= 90 and grad[i, j] <= 135) or (grad[i, j] >= -90 and grad[i, j] < -45) ):
                y_bot = np.array( [ g_mag[i+1, j], g_mag[i+1, j-1] ] )
                y_top = np.array( [ g_mag[i-1, j], g_mag[i-1, j+1] ] )
                # cos(theta) = x / magnitude
                x_est = np.absolute( g_x[i, j] / g_mag[i, j] )
                if(g_mag[i,j] >= ( (y_bot[1]-y_bot[0])*x_est + y_bot[0] ) and g_mag[i, j] >= ( (y_top[1]-y_top[0])*x_est + y_top[0] ) ):
                    NMS[i,j] = g_mag[i,j]
                else:
                    NMS[i,j] = 0
             # 135 <= gradient < 180 and  -45 <= gradient < 0
            if( (grad[i, j] >= 135 and grad[i, j] <= 180) or (grad[i, j] >= -45 and grad[i, j] < 360) ):
                y_bot = np.array( [ g_mag[i, j-1], g_mag[i+1, j-1] ] )
                y_top = np.array( [ g_mag[i, j+1], g_mag[i-1, j+1] ] )
                x_est = np.absolute( g_y[i, j] / g_mag[i, j] )
                if(g_mag[i,j] >= ( (y_bot[1]-y_bot[0])*x_est + y_bot[0] ) and g_mag[i, j] >= ( (y_top[1]-y_top[0])*x_est + y_top[0] ) ):
                    NMS[i,j] = g_mag[i,j]
                else:
                    NMS[i,j] = 0
    return NMS

NMS = NonMaxSuppressionWithInterpol(magnitude, theta, vertical_edge ,horizontal_edge)
NMS = Normalize(NMS)
plt.imshow(NMS, cmap = plt.get_cmap('gray'))
plt.show

def double_threshold(img, low_threshold=0.05, high_threshold=0.09):
    high_threshold = img.max() * high_threshold_ratio
    low_threshold = high_threshold * low_threshold_ratio

    res = np.zeros(img.shape)

    weak = np.int32(25)
    strong = np.int32(255)

    strong_row, strong_col = np.where(img >= high_threshold)
    zeros_row, zeros_col = np.where(img < low_threshold)

    weak_row, weak_col = np.where((img <= high_threshold) & (img >= low_threshold))

    res[strong_row, strong_col] = strong
    res[weak_row, weak_col] = weak

    return res

def hysteresis(img, weak=25, strong=255):
    M, N = img.shape
    for i in range(1, M-1):
        for j in range(1, N-1):
            if (img[i,j] == weak):
                try:
                    if ((img[i+1, j-1] == strong) or (img[i+1, j] == strong) or (img[i+1, j+1] == strong)
                        or (img[i, j-1] == strong) or (img[i, j+1] == strong)
                        or (img[i-1, j-1] == strong) or (img[i-1, j] == strong) or (img[i-1, j+1] == strong)):
                        img[i, j] = strong
                    else:
                        img[i, j] = 0
                except IndexError as e:
                    pass
    return img

final_image = double_thresholding(NMS)
final_image = hysteresis(final_image)
final_image = Normalize(final_image)
plt.imshow(final_image, cmap = plt.get_cmap('gray'))
plt.show()

Acknowledgement

The author would like to thank Ekbana Solution and the AI team.