Canny edge detector implementation in Python¶
Faisal Qureshi
Professor
Faculty of Science
Ontario Tech University
Oshawa ON Canada
http://vclab.science.ontariotechu.ca
Copyright information¶
© Faisal Qureshi
License¶
This work is licensed under a Creative Commons Attribution-NonCommercial 4.0 International License.
Lesson Plan¶
- Sobel filters
- Image gradient magnitude
- Using color to depict edge orientation
- Non-maxima suppression
- Hysteresis and edge-linking
- Difference of Gaussian
In [1]:
import cv2 as cv
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
Developing our own Canny edge detector¶
In [2]:
filename = 'data/lena.png'
In [3]:
weak = 64
strong = 255
hysteresis_value = 128
weak_color=[0.0,0.3,0.6]
strong_color=[1,0,0]
hysteresis_color=[0,1,0]
lowThresholdRatio=0.05
highThresholdRatio=0.09
In [4]:
img_bgr = cv.imread(filename)
img_rgb = cv.cvtColor(img_bgr, cv.COLOR_BGR2RGB)
plt.figure(figsize=(7,7))
plt.imshow(img_rgb)
plt.title(filename)
plt.xticks([])
plt.yticks([]);
In [5]:
img_gray = cv.cvtColor(img_rgb, cv.COLOR_RGB2GRAY)
img = img_gray.astype('float32') / 255.0
plt.figure(figsize=(10,10))
plt.imshow(img, cmap='gray')
plt.yticks([])
plt.xticks([]);
Sobel filters to compute image derivatives¶
In [6]:
sobelx = cv.Sobel(img,cv.CV_64F,1,0,ksize=5)
sobely = cv.Sobel(img,cv.CV_64F,0,1,ksize=5)
plt.figure(figsize=(20,10))
plt.subplot(121)
plt.imshow(sobelx, cmap='gray')
plt.title('$I_x$')
plt.xticks([])
plt.yticks([])
plt.subplot(122)
plt.title('$I_y$')
plt.imshow(sobely, cmap='gray')
plt.xticks([])
plt.yticks([]);
Computing image gradient magnitude¶
In [7]:
g = np.sqrt(np.square(sobelx)+np.square(sobely))
plt.figure(figsize=(10,10))
plt.xticks([])
plt.yticks([])
plt.imshow(g, cmap='gray');
In [8]:
plt.figure(figsize=(30,10))
plt.subplot(131)
plt.xticks([])
plt.yticks([])
plt.imshow(img_rgb)
plt.subplot(132)
plt.imshow(img, cmap='gray')
plt.xticks([])
plt.yticks([])
plt.subplot(133)
plt.xticks([])
plt.yticks([])
plt.imshow(g, cmap='gray');
Using color to depict edge orientation¶
In [9]:
def colorize_gradient_orientation(g, sobelx, sobely):
hsv = np.zeros([g.shape[0],g.shape[0],3], dtype=np.uint8)
hsv[..., 1] = 255
mag, ang = cv.cartToPolar(sobelx[...], sobely[...])
#print(hsv.shape)
#print(mag.shape)
#print(ang.shape)
hsv[..., 0] = ang * 180 / np.pi / 2
hsv[..., 2] = cv.normalize(mag, None, 0, 255, cv.NORM_MINMAX)
gc = cv.cvtColor(hsv, cv.COLOR_HSV2BGR)
return gc
In [10]:
gc = colorize_gradient_orientation(g, sobelx, sobely)
plt.figure(figsize=(10,10))
plt.imshow(gc)
plt.xticks([])
plt.yticks([]);
Non-maxima suppression¶
Use non-maxima suppression for better edge localization. Notice thick edges in the zoomed in picture on the right.
In [11]:
def crop_and_highlight(img, crop):
rs = crop[0]
re = crop[1]
cs = crop[2]
ce = crop[3]
assert(rs <= re)
assert(cs <= ce)
assert(rs >=0 and re <= img.shape[0])
assert(cs >=0 and ce <= img.shape[1])
img_cropped = img[rs:re, cs:ce, ...]
if re == img.shape[0]:
re = re - 1
if ce == img.shape[1]:
ce = ce - 1
x = [cs, cs, ce, ce, cs]
y = [rs, re, re, rs, rs]
return img_cropped, x, y
In [12]:
crop1 = (0, 150, 512-150, 512)
In [13]:
g_crop1, x, y = crop_and_highlight(g, crop1)
gc_crop1, _, _ = crop_and_highlight(gc, crop1)
sobelx_crop1, _, _ = crop_and_highlight(sobelx, crop1)
sobely_crop1, _, _ = crop_and_highlight(sobely, crop1)
plt.figure(figsize=(20,10))
plt.subplot(121)
plt.imshow(gc, cmap='gray')
plt.plot(x, y, 'r-', linewidth=2)
plt.xticks([])
plt.yticks([])
plt.subplot(122)
plt.xticks([])
plt.yticks([])
plt.imshow(gc_crop1, cmap='gray');
In [14]:
crop2 = (50, 75, 75, 100)
In [15]:
g_crop2, x, y = crop_and_highlight(g_crop1, crop2)
gc_crop2, _, _ = crop_and_highlight(gc_crop1, crop2)
sobelx_crop2, _, _ = crop_and_highlight(sobelx_crop1, crop2)
sobely_crop2, _, _ = crop_and_highlight(sobely_crop1, crop2)
plt.figure(figsize=(20,10))
plt.subplot(121)
plt.imshow(gc_crop1, cmap='gray')
plt.plot(x, y, 'r-')
plt.xticks([])
plt.yticks([])
plt.subplot(122)
plt.imshow(gc_crop2, cmap='gray', interpolation='None')
plt.xticks(np.linspace(0,24,25))
plt.yticks(np.linspace(0,24,25));
The key idea for non-maxima suppression is: for each edge pixel, look across the edge and only keep pixel that has the highest intensity, zeroing out all other pixels.
In [16]:
def line_endpoints(x, y, dx, dy, scale=0.1):
xx = [x - scale*dx, x + scale*dx]
yy = [y - scale*dy, y + scale*dy]
return xx, yy
In [17]:
def sample_nn_along_line(img, x, y, dx, dy, width=3):
try:
values = np.empty([1, 2*width+1, img.shape[2]], dtype=img.dtype)
except:
values = np.empty([1, 2*width+1], dtype=img.dtype)
l = np.linalg.norm(np.array([dx, dy]))
for i in range(-width, width+1):
x1 = int(x + i*dx/l)
y1 = int(y + i*dy/l)
values[0, i + width, ...] = img[x1, y1, ...]
return values
In [18]:
x, y, scale, = 15, 14, 0.2
dx, dy = sobelx_crop2[x, y], sobely_crop2[x, y]
xx, yy = line_endpoints(x, y, dx, dy, scale)
plt.figure(figsize=(10,10))
plt.subplot(111)
plt.imshow(g_crop2, interpolation='None', cmap='gray')
plt.xticks(np.linspace(0,24,25))
plt.xlabel('x')
plt.ylabel('y')
plt.yticks(np.linspace(0,24,25))
plt.plot(xx, yy, 'r-')
plt.arrow(x, y, scale*dx, scale*dy, head_width=1, head_length=1, fc='r', ec='r')
plt.scatter(x, y, color='blue');
#plt.subplot(121)
In [19]:
width = 4
bgr_vals = sample_nn_along_line(gc_crop2, x, y, dx, dy, width)
g_vals = sample_nn_along_line(g_crop2, x, y, dx, dy, width)
m = np.argmax(g_vals) # Just for generating picture for slides +1
#print(m)
#print(g_vals1.shape)
g_vals1 = np.copy(g_vals)
#print(g_vals1)
if not m == width:
g_vals1[0,m] = 0
plt.figure(figsize=(10,5))
plt.subplot(311)
plt.imshow(g_vals, cmap='gray')
plt.scatter([width],[0])
plt.yticks([])
plt.xticks([])
plt.subplot(312)
plt.xticks([])
plt.plot(g_vals[0], 'r-')
plt.scatter([width], g_vals[0,width])
plt.subplot(313)
plt.plot(g_vals1[0], 'r-')
#plt.scatter([m], g_vals1[0,m], color='black')
plt.text(0,(np.max(g_vals) - np.min(g_vals))/2,'After non-maxima suppression');
In [20]:
x, y, scale, = 12, 2, 0.2
dx, dy = sobelx_crop2[x, y], sobely_crop2[x, y]
xx, yy = line_endpoints(x, y, dx, dy, scale)
plt.figure(figsize=(10,10))
plt.subplot(111)
plt.imshow(g_crop2, interpolation='None', cmap='gray')
plt.xticks(np.linspace(0,24,25))
plt.xlabel('x')
plt.ylabel('y')
plt.yticks(np.linspace(0,24,25))
plt.plot(xx, yy, 'r-')
plt.arrow(x, y, scale*dx, scale*dy, head_width=1, head_length=1, fc='r', ec='r')
plt.scatter(x, y, color='blue');
#plt.subplot(121)
In [21]:
width = 4
bgr_vals = sample_nn_along_line(gc_crop2, x, y, dx, dy, width)
g_vals = sample_nn_along_line(g_crop2, x, y, dx, dy, width)
m = np.argmax(g_vals)
#print(m)
#print(g_vals1.shape)
g_vals1 = np.copy(g_vals)
#print(g_vals1)
if not m == width:
g_vals1[0,width] = 0
plt.figure(figsize=(10,5))
plt.subplot(311)
plt.imshow(g_vals, cmap='gray')
plt.scatter([width],[0])
plt.yticks([])
plt.xticks([])
plt.subplot(312)
plt.xticks([])
plt.plot(g_vals[0], 'r-')
plt.scatter([width], g_vals[0,width])
plt.subplot(313)
plt.plot(g_vals1[0], 'r-')
#plt.scatter([m], g_vals1[0,width], color='black')
plt.text(0,(np.max(g_vals) - np.min(g_vals))/2,'After non-maxima suppression');
Figures for course slides
In [22]:
gc_crop1, x1, y1 = crop_and_highlight(gc, crop1)
sobelx_crop1, _, _ = crop_and_highlight(sobelx, crop1)
sobely_crop1, _, _ = crop_and_highlight(sobely, crop1)
gc_crop2, x2, y2 = crop_and_highlight(gc_crop1, crop2)
sobelx_crop2, _, _ = crop_and_highlight(sobelx_crop1, crop2)
sobely_crop2, _, _ = crop_and_highlight(sobely_crop1, crop2)
plt.figure(figsize=(30,90))
plt.subplot(131)
plt.imshow(gc)
plt.plot(x1, y1, 'r-', linewidth=2)
plt.xticks([])
plt.yticks([])
plt.subplot(132)
plt.imshow(gc_crop1)
plt.plot(x2, y2, 'r-', linewidth=2)
plt.xticks([])
plt.yticks([])
plt.subplot(133)
x, y, scale = 15, 14, 0.2
dx, dy = sobelx_crop2[x, y], sobely_crop2[x, y]
xx, yy = line_endpoints(x, y, dx, dy, scale)
plt.imshow(gc_crop2)
plt.xticks(np.linspace(0,24,25))
plt.xlabel('x')
plt.ylabel('y')
plt.yticks(np.linspace(0,24,25))
plt.plot(xx, yy, 'r-')
plt.arrow(x, y, scale*dx, scale*dy, head_width=1, head_length=1, fc='r', ec='r')
plt.scatter(x, y, color='blue');
Implementation¶
In [23]:
import math
def sample_neighbours(grad, x, y, dx, dy, width):
if math.isclose(dx, 0.0, rel_tol=1e-5) and math.isclose(dy, 0.0, rel_tol=1e-5):
return []
M, N = grad.shape
#print(M, N)
neighbours = []
l = np.linalg.norm(np.array([dx, dy], dtype=np.float32))
for i in range(-width, width+1):
x1 = int(x + i*dx/l)
y1 = int(y + i*dy/l)
if x1 >= 0 and x1 < M and y1 >= 0 and y1 < N:
neighbours.append((x1, y1, grad[x1, y1]))
return neighbours
Testing sampling
In [24]:
width = 5
x, y = 15, 14
dx, dy = sobelx_crop2[x, y], sobely_crop2[x, y]
neighbours = sample_neighbours(g_crop2, x, y, dx, dy, width)
for i in range(len(neighbours)):
print(i, neighbours[i])
0 (18, 17, np.float64(0.8009413437795406)) 1 (18, 16, np.float64(0.24615457891281892)) 2 (17, 15, np.float64(8.292644617706165)) 3 (16, 15, np.float64(16.33219594008025)) 4 (15, 14, np.float64(28.119370722233324)) 5 (15, 14, np.float64(28.119370722233324)) 6 (14, 13, np.float64(15.594596325917706)) 7 (13, 12, np.float64(4.536588559993903)) 8 (12, 12, np.float64(2.2308175381431625)) 9 (11, 11, np.float64(0.08318902812231753)) 10 (11, 10, np.float64(0.35537286538155094))
Testing find maximum
In [25]:
def find_maximum(neighbours):
if len(neighbours) == 0:
return -1
ind = 0
for i in range(len(neighbours)):
if neighbours[ind][2] < neighbours[i][2]:
ind = i
return ind
In [26]:
#neighbours = [(1,2,3),(3,4,5),(6,7,1)]
#print(neighbours[0][0])
ind = find_maximum(neighbours)
print(ind)
print(neighbours[ind])
4 (15, 14, np.float64(28.119370722233324))
In [27]:
def non_max_suppression1(grad, sobelx, sobely, width):
M, N = grad.shape
Z = np.zeros(grad.shape, dtype=grad.dtype)
# rows
for i in range(M):
# columns
for j in range(N):
dx, dy = sobelx[i,j], sobely[i,j]
neighbours = sample_neighbours(grad, i, j, dx, dy, width)
ind = find_maximum(neighbours)
if ind > -1:
x, y, _ = neighbours[ind]
# Keep this pixel if this the maxima, otherwise set it to 0
if i == x and j == y:
Z[i, j] = grad[i, j]
return Z
In [28]:
width = 2
g_thinned = non_max_suppression1(g, sobelx, sobely, width)
In [29]:
g_thinned_crop1, x1, y1 = crop_and_highlight(g_thinned, crop1)
g_thinned_crop2, x2, y2 = crop_and_highlight(g_thinned_crop1, crop2)
plt.figure(figsize=(30,90))
plt.subplot(131)
plt.imshow(g_thinned, cmap='gray')
plt.plot(x1, y1, 'r-', linewidth=2)
plt.xticks([])
plt.yticks([])
plt.subplot(132)
plt.imshow(g_thinned_crop1, cmap='gray')
plt.plot(x2, y2, 'r-', linewidth=2)
plt.xticks([])
plt.yticks([])
plt.subplot(133)
plt.imshow(g_thinned_crop2, cmap='gray')
plt.xticks(np.linspace(0,24,25))
plt.xlabel('x')
plt.ylabel('y')
plt.yticks(np.linspace(0,24,25));
In [30]:
gc_thinned = np.copy(gc)
gc_thinned[g_thinned == 0] = 0
gc_thinned_crop1, x1, y1 = crop_and_highlight(gc_thinned, crop1)
gc_thinned_crop2, x2, y2 = crop_and_highlight(gc_thinned_crop1, crop2)
plt.figure(figsize=(30,90))
plt.subplot(131)
plt.imshow(gc_thinned, cmap='gray')
plt.plot(x1, y1, 'r-', linewidth=2)
plt.xticks([])
plt.yticks([])
plt.subplot(132)
plt.imshow(gc_thinned_crop1, cmap='gray')
plt.plot(x2, y2, 'r-', linewidth=2)
plt.xticks([])
plt.yticks([])
plt.subplot(133)
plt.imshow(gc_thinned_crop2, cmap='gray')
plt.xticks(np.linspace(0,24,25))
plt.xlabel('x')
plt.ylabel('y')
plt.yticks(np.linspace(0,24,25));
Dividing pixels into weak and strong edge pixels¶
In [31]:
# See also https://towardsdatascience.com/canny-edge-detection-step-by-step-in-python-computer-vision-b49c3a2d8123
def threshold(grad, lowThresholdRatio, highThresholdRatio, weak, strong):
"""
Set pixels as follows:
pixels whose grad magnitude < lowThresholdRatio*grad.max() --> 0
pixels whose grad magnitude > highThresholdRatio*grad.max() --> strong
otherwise --> weak
"""
highThreshold = grad.max() * highThresholdRatio;
lowThreshold = highThreshold * lowThresholdRatio;
M, N = grad.shape
res = np.zeros((M,N), dtype=np.int32)
weak = np.int32(weak)
strong = np.int32(strong)
strong_i, strong_j = np.where(grad >= highThreshold)
zeros_i, zeros_j = np.where(grad < lowThreshold)
weak_i, weak_j = np.where((grad <= highThreshold) & (grad >= lowThreshold))
res[strong_i, strong_j] = strong
res[weak_i, weak_j] = weak
return (res, weak, strong)
In [32]:
def colorize_thresholded(g_thinned_thresholded, weak, strong, hysteresis_value, weak_color=weak_color, strong_color=strong_color, hysteresis_color=hysteresis_color):
g = g_thinned_thresholded
c = np.zeros((g.shape[0], g.shape[1], 3), dtype=np.float32)
c[g == weak] = np.array(weak_color)
c[g == strong] = np.array(strong_color)
c[g == hysteresis_value] = np.array(hysteresis_color)
return c
In [33]:
g_thinned_thresholded, _, _ = threshold(g_thinned, lowThresholdRatio, highThresholdRatio, weak, strong)
gc_thinned_thresholded = colorize_thresholded(g_thinned_thresholded, weak, strong, hysteresis_value)
plt.figure(figsize=(20,40))
plt.subplot(121)
plt.xticks([])
plt.yticks([])
plt.imshow(g_thinned, cmap='gray')
plt.subplot(122)
plt.xticks([])
plt.yticks([])
plt.imshow(gc_thinned_thresholded)
plt.title('Weak and strong pixels');
In [34]:
gc_thinned_thresholded_crop1, x1, y1 = crop_and_highlight(gc_thinned_thresholded, crop1)
gc_thinned_thresholded_crop2, x2, y2 = crop_and_highlight(gc_thinned_thresholded_crop1, crop2)
# gc_thinned_thresholded_crop1, x1, y1 = crop_and_highlight(gc_thinned_thresholded, crop1)
# gc_thinned_thresholded_crop2, x2, y2 = crop_and_highlight(gc_thinned_thresholded_crop1, crop2)
plt.figure(figsize=(30,90))
plt.subplot(131)
plt.imshow(gc_thinned_thresholded, cmap='gray')
plt.plot(x1, y1, 'r-', linewidth=2)
plt.xticks([])
plt.yticks([])
plt.subplot(132)
plt.imshow(gc_thinned_thresholded_crop1, cmap='gray')
plt.plot(x2, y2, 'r-', linewidth=2)
plt.xticks([])
plt.yticks([])
plt.subplot(133)
plt.imshow(gc_thinned_thresholded_crop2, cmap='gray')
plt.xticks(np.linspace(0,24,25))
plt.xlabel('x')
plt.ylabel('y')
plt.yticks(np.linspace(0,24,25));
Hysteresis and edge linking¶
The goal is to start edges at strong edge pixels. Continue edges to weak edge pixels. Do not start edges at weak edge pixels. The process of edge linking is actually quite computation intensive.
An approximation to hyesteresis¶
Turn every weak pixel to a strong pixel, if a neighbouring pixel is strong pixel. Strong pixels are our edge pixels. Notice that no real edge linking is taking place here.
In [35]:
def get_sub_region(img, r, c, w):
R, C = img.shape[0], img.shape[1]
rs, re = np.max([r-w, 0]), np.min([r+w, R])+1
cs, ce = np.max([c-w, 0]), np.min([c+w, C])+1
return img[rs:re, cs:ce], (rs, re, cs, ce)
In [36]:
def perform_hysteresis(grad, weak, strong, hysteresis_value, width):
"""
Set weak pixels equal to 0, if no pixel around it in area with radius
"width" is strong.
We use hysteresis value to identify pixels that will be set to strong
due to hysteresis operation where a weak pixel is set to strong if it is
next to a strong pixel. We use this value to display the effects of this
and edge linking operation.
"""
img = np.copy(grad)
M, N = img.shape
for i in range(1, M-1):
for j in range(1, N-1):
if (img[i,j] == weak):
sub_img, _ = get_sub_region(img, i, j, width)
if np.any(sub_img == strong):
img[i, j] = hysteresis_value
else:
img[i, j] = 0
return img
In [37]:
e = perform_hysteresis(g_thinned_thresholded, weak, strong, hysteresis_value, width)
ec = colorize_thresholded(e, weak, strong, hysteresis_value)
plt.figure(figsize=(20,40))
plt.subplot(121)
plt.xticks([])
plt.yticks([])
plt.imshow(gc_thinned_thresholded)
plt.subplot(122)
plt.xticks([])
plt.yticks([])
plt.imshow(ec);
In [38]:
crop3 = (350, 500, 150, 300)
crop4 = (45, 70, 95, 120)
g_thinned_thresholded_crop3, x3, y3 = crop_and_highlight(g_thinned_thresholded, crop3)
g_thinned_thresholded_crop4, x4, y4 = crop_and_highlight(g_thinned_thresholded_crop3, crop4)
gc_thinned_thresholded_crop3, _, _ = crop_and_highlight(gc_thinned_thresholded, crop3)
gc_thinned_thresholded_crop4, _, _ = crop_and_highlight(gc_thinned_thresholded_crop3, crop4)
ec_crop3, _, _ = crop_and_highlight(ec, crop3)
ec_crop4, _, _ = crop_and_highlight(ec_crop3, crop4)
plt.figure(figsize=(30,20))
#plt.suptitle('Performing hysteresis using strong and weak edge pixels')
plt.subplot(231)
plt.imshow(gc_thinned_thresholded)
plt.plot(x3, y3, 'c-', linewidth=2)
plt.xticks([])
plt.yticks([])
plt.subplot(232)
plt.imshow(gc_thinned_thresholded_crop3)
plt.plot(x4, y4, 'c-', linewidth=2)
plt.xticks([])
plt.yticks([])
plt.subplot(233)
plt.imshow(gc_thinned_thresholded_crop4)
plt.xticks(np.linspace(0,24,25))
plt.xlabel('x')
plt.ylabel('y')
plt.yticks(np.linspace(0,24,25))
#plt.title('Blue represents weak edge pixels, and red represents strong edge pixels.')
plt.subplot(234)
plt.imshow(ec)
plt.plot(x3, y3, 'c-', linewidth=2)
plt.xticks([])
plt.yticks([])
plt.subplot(235)
plt.imshow(ec_crop3)
plt.plot(x4, y4, 'c-', linewidth=2)
plt.xticks([])
plt.yticks([])
plt.subplot(236)
plt.imshow(ec_crop4)
plt.xticks(np.linspace(0,24,25))
plt.xlabel('x')
plt.ylabel('y')
plt.yticks(np.linspace(0,24,25));
#plt.title('After performing hysteresis. Red represents strong edge pixels, and green represents weak edge pixels that have been relabelled as strong edge pixels.')
Edgelinking, starting from strong pixels.¶
We will start with a strong edge pixel and we check for other pixels along this edge. If there is a weak pixel, we will set it to to a strong pixel. Otherwise, we stop exploring. Use edge tangents pixel (x,y) to link neighbouring pixels. We can rotate gradient vector by 90 degrees to get the tangent. Recall that gradient is always normal to an edge.
In [39]:
# Rotate [sobelx, sobely] by 90 degrees
tangentx = np.copy(sobely)
tangenty = -np.copy(sobelx)
# A quicker way to normalize tangents
# tangent = np.stack([tangentx, tangenty], axis=2)
# tangent_norm = np.linalg.norm(tangent, axis=2)
# ind = tangent_norm != 0
# tangentx[ind] = tangentx[ind] / tangent_norm[ind]
# tangenty[ind] = tangenty[ind] / tangent_norm[ind]
tangentx_crop3, _, _ = crop_and_highlight(tangentx, crop3)
tangenty_crop3, _, _ = crop_and_highlight(tangenty, crop3)
tangentx_crop4, _, _ = crop_and_highlight(tangentx_crop3, crop4)
tangenty_crop4, _, _ = crop_and_highlight(tangenty_crop3, crop4)
#gc_thinned_crop3, _, _ = crop_and_highlight(gc_thinned, crop3)
#gc_thinned_crop4, _, _ = crop_and_highlight(gc_thinned_crop3, crop4)
In [40]:
i, j = 16, 11
plt.figure(figsize=(30,90))
plt.subplot(131)
plt.imshow(gc_thinned_thresholded, cmap='gray')
plt.plot(x3, y3, 'c-', linewidth=2)
plt.xticks([])
plt.yticks([])
plt.subplot(132)
plt.imshow(gc_thinned_thresholded_crop3, cmap='gray')
plt.plot(x4, y4, 'c-', linewidth=2)
plt.xticks([])
plt.yticks([])
plt.subplot(133)
plt.imshow(gc_thinned_thresholded_crop4, cmap='gray')
plt.xticks(np.linspace(0,24,25))
plt.xlabel('x')
plt.ylabel('y')
plt.yticks(np.linspace(0,24,25))
x, y, scale = j, i, 3
dx, dy = tangentx_crop4[x, y], tangenty_crop4[x, y]
xx, yy = line_endpoints(x, y, dx, dy, scale)
plt.plot(xx, yy, 'c-', linewidth=2)
#plt.arrow(x, y, scale*dx, scale*dy, head_width=0, head_length=0, fc='r', ec='r', linewidth=5)
plt.scatter(x, y, color='white');
In [41]:
i, j = 15, 2
plt.figure(figsize=(30,90))
plt.subplot(131)
plt.imshow(gc_thinned_thresholded, cmap='gray')
plt.plot(x3, y3, 'c-', linewidth=2)
plt.xticks([])
plt.yticks([])
plt.subplot(132)
plt.imshow(gc_thinned_thresholded_crop3, cmap='gray')
plt.plot(x4, y4, 'c-', linewidth=2)
plt.xticks([])
plt.yticks([])
plt.subplot(133)
plt.imshow(gc_thinned_thresholded_crop4, cmap='gray')
plt.xticks(np.linspace(0,24,25))
plt.xlabel('x')
plt.ylabel('y')
plt.yticks(np.linspace(0,24,25))
x, y, scale = j, i, 3
dx, dy = tangentx_crop4[x, y], tangenty_crop4[x, y]
xx, yy = line_endpoints(x, y, dx, dy, scale)
plt.plot(xx, yy, 'c-', linewidth=2)
#plt.arrow(x, y, scale*dx, scale*dy, head_width=0, head_length=0, fc='r', ec='r', linewidth=5)
plt.scatter(x, y, color='white');
Edge linking requires that we explore neighbouring pixels along the edge
In [42]:
def find_neighbours_along_edge(g_thinned_thresholded, i, j, tx, ty, flags, weak, strong, hysteresis_value):
"""
Returns a list of neighbouring pixels coordinates that should be explored.
These pixels are neighbours of (i, j) along the edge as determined by tangent
(tx, ty) at pixel (i, j). The pixel values are either weak or strong.
Please remember that it returns (x,y) coordinates where x -> columns and y -> rows.
"""
I = g_thinned_thresholded
M, N = I.shape[0], I.shape[1]
t = [ tx[i,j], ty[i,j] ]
t = t / np.linalg.norm(t)
x, y, dx, dy = j, i, t[0], t[1]
c = []
fns = [np.floor, np.ceil]
for s in [1.414, -1.414, 1, -1]:
for f in fns:
x1, y1 = f(x + s*dx), f(y+s*dy)
x1, y1 = int(x1), int(y1)
if x1 >= 0 and x1 < N and y1 >= 0 and y1 < M:
if (x1 == x+1 or x1 == x or x1 == x-1) and (y1 == y+1 or y1 == y or y1 == y-1):
if not I[y1, x1] == 0:
if not (x1 == x and y1 == y) and flags[y1,x1] == False:
c.append( (x1, y1) )
flags[y1, x1] = True
return c
In [43]:
i, j = 16, 11
scale = 2
M, N = g_thinned_thresholded_crop4.shape[0], g_thinned_thresholded_crop4.shape[1]
flags_crop4 = np.empty((M,N), dtype=bool)
flags_crop4[:] = False
n = find_neighbours_along_edge(g_thinned_thresholded_crop4, i, j, tangentx_crop4, tangenty_crop4, flags_crop4, weak=64, strong=255, hysteresis_value=200)
x, y = j, i
dx, dy = tangentx_crop4[i, j], tangenty_crop4[i, j]
plt.figure(figsize=(10,10))
plt.imshow(gc_thinned_thresholded_crop4, cmap='gray')
xx, yy = line_endpoints(x, y, dx, dy, scale)
plt.plot(xx, yy, 'w-', linewidth=2)
plt.plot(x, y, 'c.')
plt.scatter(x, y, color='white')
for c in n:
plt.plot(c[0], c[1], 'c.')
plt.xticks(np.linspace(0,g_thinned_thresholded_crop4.shape[0]-1,g_thinned_thresholded_crop4.shape[0]))
plt.yticks(np.linspace(0,g_thinned_thresholded_crop4.shape[1]-1,g_thinned_thresholded_crop4.shape[1]));
In [44]:
i, j = 15, 12
scale = 2
M, N = g_thinned_thresholded_crop4.shape[0], g_thinned_thresholded_crop4.shape[1]
flags_crop4 = np.empty((M,N), dtype=bool)
flags_crop4[:] = False
n = find_neighbours_along_edge(g_thinned_thresholded_crop4, i, j, tangentx_crop4, tangenty_crop4, flags_crop4, weak=64, strong=255, hysteresis_value=200)
x, y = j, i
dx, dy = tangentx_crop4[i, j], tangenty_crop4[i, j]
plt.figure(figsize=(10,10))
plt.imshow(gc_thinned_thresholded_crop4, cmap='gray')
xx, yy = line_endpoints(x, y, dx, dy, scale)
plt.plot(xx, yy, 'w-', linewidth=2)
plt.plot(x, y, 'c.')
plt.scatter(x, y, color='white')
for c in n:
plt.plot(c[0], c[1], 'c.')
plt.xticks(np.linspace(0,g_thinned_thresholded_crop4.shape[0]-1,g_thinned_thresholded_crop4.shape[0]))
plt.yticks(np.linspace(0,g_thinned_thresholded_crop4.shape[1]-1,g_thinned_thresholded_crop4.shape[1]));
In [45]:
i, j = 24, 1
scale = .2
M, N = g_thinned_thresholded_crop4.shape[0], g_thinned_thresholded_crop4.shape[1]
flags_crop4 = np.empty((M,N), dtype=bool)
flags_crop4[:] = False
n = find_neighbours_along_edge(g_thinned_thresholded_crop4, i, j, tangentx_crop4, tangenty_crop4, flags_crop4, weak=64, strong=255, hysteresis_value=200)
x, y = j, i
dx, dy = tangentx_crop4[i, j], tangenty_crop4[i, j]
plt.figure(figsize=(10,10))
plt.imshow(gc_thinned_thresholded_crop4, cmap='gray')
xx, yy = line_endpoints(x, y, dx, dy, scale)
plt.plot(xx, yy, 'w-', linewidth=2)
plt.plot(x, y, 'c.')
plt.scatter(x, y, color='white')
for c in n:
plt.plot(c[0], c[1], 'c.')
plt.xticks(np.linspace(0,g_thinned_thresholded_crop4.shape[0]-1,g_thinned_thresholded_crop4.shape[0]))
plt.yticks(np.linspace(0,g_thinned_thresholded_crop4.shape[1]-1,g_thinned_thresholded_crop4.shape[1]));
In [46]:
def link_edges_using_hysteresis(g_thinned_thresholded, tx, ty, weak, strong, hysteresis_value):
I = np.copy(g_thinned_thresholded)
M, N = I.shape[0], I.shape[1]
flags = np.empty((M,N), dtype=bool)
flags[:] = False
for i in range(M):
for j in range(N):
#print('Processing ({},{})'.format(i,j))
if I[i,j] == strong and flags[i,j] == False:
flags[i, j] = True
n = find_neighbours_along_edge(I, i, j, tx, ty, flags, weak=weak, strong=strong, hysteresis_value=hysteresis_value)
while len(n)>0:
p = n.pop()
i1, j1 = p[1], p[0]
#print('Popping ({},{})'.format(i1, j1))
if not I[i1, j1] == strong:
I[i1, j1] = hysteresis_value
nx = find_neighbours_along_edge(I, i1, j1, tx, ty, flags, weak, strong, hysteresis_value)
n.extend(nx)
# if flags[24,1] == True:
# print(i,j)
# return I
I[flags == False] = 0
return I
Testing edge linking
In [47]:
el_crop4 = link_edges_using_hysteresis(g_thinned_thresholded_crop4, tangentx_crop4, tangenty_crop4, weak, strong, hysteresis_value)
elc_crop4 = colorize_thresholded(el_crop4, weak, strong, hysteresis_value)
plt.imshow(elc_crop4);
In [48]:
el = link_edges_using_hysteresis(g_thinned_thresholded, tangentx, tangenty, weak=weak, strong=strong, hysteresis_value=hysteresis_value)
elc = colorize_thresholded(el, weak, strong, hysteresis_value)
elc_crop3, _, _ = crop_and_highlight(elc, crop3)
elc_crop4, _, _ = crop_and_highlight(elc_crop3, crop4)
plt.figure(figsize=(30,30))
#plt.suptitle('Performing hysteresis using strong and weak edge pixels')
plt.subplot(331)
plt.imshow(gc_thinned_thresholded)
plt.plot(x3, y3, 'c-', linewidth=2)
plt.xticks([])
plt.yticks([])
plt.subplot(332)
plt.imshow(gc_thinned_thresholded_crop3)
plt.plot(x4, y4, 'c-', linewidth=2)
plt.xticks([])
plt.yticks([])
plt.subplot(333)
plt.imshow(gc_thinned_thresholded_crop4)
plt.xticks(np.linspace(0,24,25))
plt.xlabel('x')
plt.ylabel('y')
plt.yticks(np.linspace(0,24,25))
#plt.title('Blue represents weak edge pixels, and red represents strong edge pixels.')
plt.subplot(334)
plt.imshow(ec)
plt.plot(x3, y3, 'c-', linewidth=2)
plt.xticks([])
plt.yticks([])
plt.subplot(335)
plt.imshow(ec_crop3)
plt.plot(x4, y4, 'c-', linewidth=2)
plt.xticks([])
plt.yticks([])
plt.subplot(336)
plt.imshow(ec_crop4)
plt.xticks(np.linspace(0,24,25))
plt.xlabel('x')
plt.ylabel('y')
plt.yticks(np.linspace(0,24,25))
#plt.title('After performing hysteresis. Red represents strong edge pixels, and green represents weak edge pixels that have been relabelled as strong edge pixels.')
plt.subplot(337)
plt.imshow(elc)
plt.plot(x3, y3, 'c-', linewidth=2)
plt.xticks([])
plt.yticks([])
plt.subplot(338)
plt.imshow(elc_crop3)
plt.plot(x4, y4, 'c-', linewidth=2)
plt.xticks([])
plt.yticks([])
plt.subplot(339)
plt.imshow(elc_crop4)
plt.xticks(np.linspace(0,24,25))
plt.xlabel('x')
plt.ylabel('y')
plt.yticks(np.linspace(0,24,25));
#plt.title('After performing hysteresis. Red represents strong edge pixels, and green represents weak edge pixels that have been relabelled as strong edge pixels.')
In [49]:
plt.figure(figsize=(30,20))
#plt.suptitle('Performing hysteresis using strong and weak edge pixels')
plt.subplot(231)
plt.imshow(gc_thinned_thresholded)
plt.plot(x3, y3, 'c-', linewidth=2)
plt.xticks([])
plt.yticks([])
plt.subplot(232)
plt.imshow(gc_thinned_thresholded_crop3)
plt.plot(x4, y4, 'c-', linewidth=2)
plt.xticks([])
plt.yticks([])
plt.subplot(233)
plt.imshow(gc_thinned_thresholded_crop4)
plt.xticks(np.linspace(0,24,25))
plt.xlabel('x')
plt.ylabel('y')
plt.yticks(np.linspace(0,24,25))
#plt.title('Blue represents weak edge pixels, and red represents strong edge pixels.')
plt.subplot(234)
plt.imshow(elc)
plt.plot(x3, y3, 'c-', linewidth=2)
plt.xticks([])
plt.yticks([])
plt.subplot(235)
plt.imshow(elc_crop3)
plt.plot(x4, y4, 'c-', linewidth=2)
plt.xticks([])
plt.yticks([])
plt.subplot(236)
plt.imshow(elc_crop4)
plt.xticks(np.linspace(0,24,25))
plt.xlabel('x')
plt.ylabel('y')
plt.yticks(np.linspace(0,24,25));
#plt.title('After performing hysteresis. Red represents strong edge pixels, and green represents weak edge pixels that have been relabelled as strong edge pixels.')
Putting it all together.¶
In [50]:
e_final = np.copy(e)
e_final[e == hysteresis_value] = strong
el_final = np.copy(el)
el_final[el == hysteresis_value] = strong
e_canny = cv.Canny(img_gray, threshold1=60, threshold2=120, apertureSize=3, L2gradient=True)
plt.figure(figsize=(20,5))
plt.subplot(141)
plt.xticks([])
plt.yticks([])
plt.imshow(img_gray, cmap='gray')
plt.subplot(142)
plt.xticks([])
plt.yticks([])
plt.title('Ours w/o edge linking')
plt.imshow(e_final, cmap='gray')
plt.subplot(143)
plt.imshow(el_final, cmap='gray')
plt.xticks([])
plt.yticks([])
plt.title('Ours w edge linking')
plt.subplot(144)
plt.imshow(e_canny, cmap='gray')
plt.xticks([])
plt.yticks([])
plt.title('OpenCV');
In [51]:
e_canny = cv.Canny(img_gray, threshold1=60, threshold2=120, L2gradient=True)
plt.figure(figsize=(20,10))
plt.subplot(121)
plt.imshow(img_rgb)
plt.xticks([])
plt.yticks([])
plt.subplot(122)
plt.imshow(e_canny, cmap='gray')
plt.xticks([])
plt.yticks([])
plt.title('OpenCV');
Difference of Gaussians¶
You can also use Difference of Gaussians to identify edge pixels for subsequent processing --- thresholding, edge linking, etc. You will still need to compute image gradients for thinning and edge linking.
In [52]:
blur5 = cv.GaussianBlur(img,(9,9),0)
blur3 = cv.GaussianBlur(img,(3,3),0)
dog = blur5 - blur3
plt.figure(figsize=(21,7))
plt.subplot(131)
plt.title(r'$9 \times 9$ filter')
plt.imshow(blur5, cmap='gray')
plt.subplot(132)
plt.title(r'$3 \times 3$ filter')
plt.imshow(blur3, cmap='gray')
plt.subplot(133)
plt.title('Difference of Gaussians')
plt.title('Diff. of Gaussians')
plt.imshow(dog, cmap='gray');
In [ ]: