- Published on

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