Task

Write a program that performs so-called canny edge detection on an image.



A possible algorithm consists of the following steps:

 Noise reduction.   May be performed by Gaussian filter. 
  
 Compute intensity gradient   (matrices 




G

x




{\displaystyle G_{x}}

 and 




G

y




{\displaystyle G_{y}}

)   and its magnitude   



G


{\displaystyle G}

:
           



G
=



G

x


2


+

G

y


2






{\displaystyle G={\sqrt {G_{x}^{2}+G_{y}^{2}}}}


May be performed by convolution of an image with Sobel operators. 
  
 Non-maximum suppression.   
For each pixel compute the orientation of intensity gradient vector:   



θ
=


a
t
a
n
2



(

G

y


,


G

x


)



{\displaystyle \theta ={\rm {atan2}}\left(G_{y},\,G_{x}\right)}

.     
Transform   angle 



θ


{\displaystyle \theta }

   to one of four directions:   0, 45, 90, 135 degrees.     
Compute new array   



N


{\displaystyle N}

:     if         



G

(

p

a


)

<
G

(
p
)

<
G

(

p

b


)



{\displaystyle G\left(p_{a}\right)<G\left(p\right)<G\left(p_{b}\right)}


where   



p


{\displaystyle p}

   is the current pixel,   




p

a




{\displaystyle p_{a}}

   and   




p

b




{\displaystyle p_{b}}

   are the two neighbour pixels in the direction of gradient,   
then     



N
(
p
)
=
G
(
p
)


{\displaystyle N(p)=G(p)}

,       otherwise   



N
(
p
)
=
0


{\displaystyle N(p)=0}

.   
Nonzero pixels in resulting array correspond to local maxima of   



G


{\displaystyle G}

   in direction   



θ
(
p
)


{\displaystyle \theta (p)}

. 
  
 Tracing edges with hysteresis.   
At this stage two thresholds for the values of   



G


{\displaystyle G}

   are introduced:   




T

m
i
n




{\displaystyle T_{min}}

   and   




T

m
a
x




{\displaystyle T_{max}}

.   
Starting from pixels with   



N
(
p
)
⩾

T

m
a
x




{\displaystyle N(p)\geqslant T_{max}}

,   
find all paths of pixels with   



N
(
p
)
⩾

T

m
i
n




{\displaystyle N(p)\geqslant T_{min}}

   and put them to the resulting image.


#!/bin/python
import numpy as np
from scipy.ndimage.filters import convolve, gaussian_filter
from scipy.misc import imread, imshow
 
def CannyEdgeDetector(im, blur = 1, highThreshold = 91, lowThreshold = 31):
	im = np.array(im, dtype=float) #Convert to float to prevent clipping values
 
	#Gaussian blur to reduce noise
	im2 = gaussian_filter(im, blur)
 
	#Use sobel filters to get horizontal and vertical gradients
	im3h = convolve(im2,[[-1,0,1],[-2,0,2],[-1,0,1]]) 
	im3v = convolve(im2,[[1,2,1],[0,0,0],[-1,-2,-1]])
 
	#Get gradient and direction
	grad = np.power(np.power(im3h, 2.0) + np.power(im3v, 2.0), 0.5)
	theta = np.arctan2(im3v, im3h)
	thetaQ = (np.round(theta * (5.0 / np.pi)) + 5) % 5 #Quantize direction
 
	#Non-maximum suppression
	gradSup = grad.copy()
	for r in range(im.shape[0]):
		for c in range(im.shape[1]):
			#Suppress pixels at the image edge
			if r == 0 or r == im.shape[0]-1 or c == 0 or c == im.shape[1] - 1:
				gradSup[r, c] = 0
				continue
			tq = thetaQ[r, c] % 4
 
			if tq == 0: #0 is E-W (horizontal)
				if grad[r, c] <= grad[r, c-1] or grad[r, c] <= grad[r, c+1]:
					gradSup[r, c] = 0
			if tq == 1: #1 is NE-SW
				if grad[r, c] <= grad[r-1, c+1] or grad[r, c] <= grad[r+1, c-1]:
					gradSup[r, c] = 0
			if tq == 2: #2 is N-S (vertical)
				if grad[r, c] <= grad[r-1, c] or grad[r, c] <= grad[r+1, c]:
					gradSup[r, c] = 0
			if tq == 3: #3 is NW-SE
				if grad[r, c] <= grad[r-1, c-1] or grad[r, c] <= grad[r+1, c+1]:
					gradSup[r, c] = 0
 
	#Double threshold
	strongEdges = (gradSup > highThreshold)
 
	#Strong has value 2, weak has value 1
	thresholdedEdges = np.array(strongEdges, dtype=np.uint8) + (gradSup > lowThreshold)
 
	#Tracing edges with hysteresis	
	#Find weak edge pixels near strong edge pixels
	finalEdges = strongEdges.copy()
	currentPixels = []
	for r in range(1, im.shape[0]-1):
		for c in range(1, im.shape[1]-1):	
			if thresholdedEdges[r, c] != 1:
				continue #Not a weak pixel
 
			#Get 3x3 patch	
			localPatch = thresholdedEdges[r-1:r+2,c-1:c+2]
			patchMax = localPatch.max()
			if patchMax == 2:
				currentPixels.append((r, c))
				finalEdges[r, c] = 1
 
	#Extend strong edges based on current pixels
	while len(currentPixels) > 0:
		newPix = []
		for r, c in currentPixels:
			for dr in range(-1, 2):
				for dc in range(-1, 2):
					if dr == 0 and dc == 0: continue
					r2 = r+dr
					c2 = c+dc
					if thresholdedEdges[r2, c2] == 1 and finalEdges[r2, c2] == 0:
						#Copy this weak pixel to final result
						newPix.append((r2, c2))
						finalEdges[r2, c2] = 1
		currentPixels = newPix
 
	return finalEdges
 
if __name__=="__main__":
	im = imread("test.jpg", mode="L") #Open image, convert to greyscale
	finalEdges = CannyEdgeDetector(im)
	imshow(finalEdges)