Image Edge Detection
The Canny edge detector is an edge detection operator that uses a multi-stage algorithm to detect a wide range of edges in images. It was developed by John F. Canny in 1986. Canny also produced a computational theory of edge detection explaining why the technique works. (Wikipedia)
The Canny edge detection algorithm is composed of 5 steps:
- Noise reduction (Gaussian blur) - to smooth out edges;
- Gradient calculation (Sobel Operators) - to detect edges;
- Gradient Magnitude Thresholding - to get rid of spurious response to edge detection;
- Double threshold - to determine potential edges;
- Edge Tracking by Hysteresis - suppress weak edges not connected to strong edges.
We've covered steps 1 and 2 in the Image Convolution section. Gradients using Sobel Kernels did a pretty good job of detecting edges in an image all on their own. Canny's algorithm, does a preprocessing step (Gaussian blur), and three post processing steps to create definite edges. Let's go through each step below.
def rgb2gray(rgb):
r, g, b = rgb[:,:,0], rgb[:,:,1], rgb[:,:,2]
gray = 0.2989 * r + 0.5870 * g + 0.1140 * b
return gray[0]
Step 1. Noise reduction (Gaussian blur)
Edge detection is affected by noise, so first we smooth out the image. The selection of the size of the Gaussian kernel with affect the performance of the edge detector. The larger the size, the lower the sensitivity to noise. A \(5\times 5\) kernel is an acceptable size for most cases.
def gaussian_kernel(kernel_size, sigma=1):
size = int(kernel_size) // 2
x, y = np.mgrid[-size:size+1, -size:size+1]
normal = 1 / (2.0 * np.pi * sigma**2)
g = np.exp(-((x**2 + y**2) / (2.0*sigma**2))) * normal
return g
# apply Gaussian Blur
img_gb = ndimage.filters.convolve(img,gaussian_kernel(5))
# display image
plt.imshow(img_gb, cmap="gray")
Step 2. Gradient calculation (Sobel Operators)
Applying Sobel operators we return the intensity gradient-norm matrix, and the gradient-angle matrix. Recall our axes are oriented as follows: the positive x-axis is down, and the positive y-axis is to the right. This is due to indexing. The angle for the gradient vector is calculated as usual by \(\theta = \arctan{\left(\frac{y}{x}\right)}\). See the diagram for a visual interpretation of \(\theta\), and the second diagram for basic angles we'll be considering.
def sobel_filters(img):
''' input: a gray scale image
output: the gray scale edge detection image, and direction matrix theta'''
Kx = np.array([[1, 2, 1], [0, 0, 0], [-1, -2, -1]], np.float32)
Ky = np.array([[1, 0, -1], [2, 0, -2], [1, 0, -1]], np.float32)
Ix = ndimage.filters.convolve(img, Kx)
Iy = ndimage.filters.convolve(img, Ky)
G = np.hypot(Ix, Iy)
G = G / G.max() * 255
theta = np.arctan2(Iy, Ix) # x is vertical, y is horizontal - angle measured off x axis.
return (G, theta)
# apply Sobel Filters
img_sobel, theta = sobel_filters(img_gb)
# display image
plt.show(img_sobel, cmap="gray")
Step 3. Gradient magnitude thresholding or lower bound cut-off suppression
Next we use an edge thinning technique known as minimum cut-off suppression of gradient magnitudes, or lower bound thresholding, or non-max suppression.
The algorithm for each pixel in the gradient image is:
-
Compare the edge strength of the current pixel with the edge strength of the pixel in the positive and negative gradient directions.
-
If the edge strength of the current pixel is the largest compared to the other pixels in the mask with the same direction (e.g., a pixel that is pointing in the y-direction will be compared to the pixel above and below it in the vertical axis), the value will be preserved. Otherwise, the value will be suppressed.
We discretize the directions to simplify the calculations. The direction of the gradient is rounded to the nearest angle: \(0^{\circ}\), \(45^{\circ}\), \(90^{\circ}\), \(135^{\circ}\), \(180^{\circ}\). (Since we only care about the positive or negative gradient direction, we don't need to worry about negative angles, instead we rotate them \(180^{\circ}\) to make them positive.) Therefore, we only need to consider which region in the right hand side of the following diagram does the \(\theta\) value live in.
Once we know the region \(\theta\) is in (i.e. we know the direction of the gradient), then for that particular pixel we compare the pixel intensity with that of it's two neighbours in the direction of the gradient. If either has larger intensity then we zero-out the current pixels intensity (i.e. we suppress that pixel as an edge pizel). Otherwise, we keep its values. The following diagram shows the two neighbours we consider in the direction of the gradient. This is used in selecting indices for \(q\) and \(r\) in the code below.
def edge_thinning(img, D):
'''input: Intensity Matrix img, and direction matrix D'''
M, N = img.shape
Z = np.zeros((M,N), dtype=np.int32)
angle = D * 180. / np.pi # convert direction matrix from radians to degrees
angle[angle < 0] += 180 # map all angles from 0 to 180
for i in range(1,M-1):
for j in range(1,N-1):
try:
q = 255
r = 255
#angle 0
if (0 <= angle[i,j] < 22.5) or (157.5 <= angle[i,j] <= 180):
q = img[i+1, j]
r = img[i-1, j]
#angle 45
elif (22.5 <= angle[i,j] < 67.5):
q = img[i-1, j-1]
r = img[i+1, j+1]
#angle 90
elif (67.5 <= angle[i,j] < 112.5):
q = img[i, j+1]
r = img[i, j-1]
#angle 135
elif (112.5 <= angle[i,j] < 157.5):
q = img[i+1, j-1]
r = img[i-1, j+1]
if (img[i,j] >= q) and (img[i,j] >= r):
Z[i,j] = img[i,j]
else:
Z[i,j] = 0
except IndexError as e:
pass
return Z
# apply edge thinning
img_edge = edge_thinning(img_sobel, theta)
# diplay image
plt.imshow(img_edge, cmap="gray")
Step 4. Double threshold
Filter out edge pixels with a weak gradient value and preserve edge pixels with a high gradient value.
def dbl_threshold(img, lowThresholdRatio=0.05, highThresholdRatio=0.09):
highThreshold = img.max() * highThresholdRatio;
lowThreshold = highThreshold * lowThresholdRatio;
M, N = img.shape
res = np.zeros((M,N), dtype=np.int32)
weak = np.int32(25)
strong = np.int32(255)
strong_i, strong_j = np.where(img >= highThreshold)
zeros_i, zeros_j = np.where(img < lowThreshold)
weak_i, weak_j = np.where((img <= highThreshold) & (img >= lowThreshold))
res[strong_i, strong_j] = strong
res[weak_i, weak_j] = weak
return (res, weak, strong)
Step 5. Edge Tracking by Hysteresis
Only keep weak edge pixels that are neighbours with a strong edge pixel. Otherwise, the weak edge pixel is likely there due to noise or colour variations, and is removed from consideration as an edge.
def hysteresis(img, weak, 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
Let's show all the steps in the full Canny Algorithm:
# load and covert image to gray scale
img_col = plt.imread("cube.jpg")
img = rgb2gray(img_col)
# apply gaussian blur
img_gb = ndimage.filters.convolve(img, gaussian_kernel(5))
# apply Sobel filters
img_sobel, theta = sobel_filters(img_gb)
# apply Edge thinning
img_nms = edge_thinning(img_sobel,theta)
# apply Double Threshold
img_dt, weak, strong = dbl_threshold(img_nms)
# edge tracking by hysteresis
img_h = hysteresis(img_dt, weak)
# display images
fig = plt.figure(figsize=(8,8))
fig.add_subplot(1,2,1)
plt.imshow(img, cmap="gray")
plt.title("Original")
fig.add_subplot(1,2,2)
plt.imshow(img_h, cmap="gray")
plt.title("Edge Detection")
