Image Derivatives¶

Faisal Qureshi
Professor
Faculty of Science
Ontario Tech University
http://vclab.science.ontariotechu.ca

Lesson Plan¶

• Computing image derivatives
• Sobel filters

Image edges and image deriatives¶

In [1]:
import numpy as np
import matplotlib.pyplot as plt

In [2]:
def sigmoid(x, s=1.0, t=0.0):
return 1./(1.+np.exp(-s*(x-t)))

In [3]:
x = np.linspace(0,10,100)
plt.plot(x, sigmoid(x, s=3, t=5))

Out[3]:
[<matplotlib.lines.Line2D at 0x123580ed0>]
In [4]:
# im = np.zeros([100,200])
# for i in range(100):
#     im[i,:] = np.hstack([sigmoid(x, s=5, t=5), 1-sigmoid(x, s=2.5, t=5)])

# The following does exactly what the for-loop shown above does
import numpy.matlib as matlab
im = matlab.repmat(np.hstack([sigmoid(x, s=5, t=5), 1-sigmoid(x, s=2.5, t=5)]), 100,1)

plt.figure(figsize=(10,7))
plt.imshow(im, cmap='gray')
plt.xlabel('Width ($x$)')
plt.ylabel('Height ($y$)')

Out[4]:
Text(0, 0.5, 'Height ($y$)')
In [5]:
plt.figure(figsize=(10,12))
plt.subplot(211)
plt.imshow(im, cmap='gray')
plt.title('Image')
plt.plot([1,199],[40,40],'r-')
plt.xlabel(r'Width ($x$)')
plt.ylabel(r'Height ($y$)')
plt.subplot(212)
plt.title('Row 40')
plt.plot(im[40,1:199],'r')
plt.xlabel(r'Width ($x$)')
plt.ylabel(r'$I(x)$')

Out[5]:
Text(0, 0.5, '$I(x)$')
In [6]:
Hx = [1,-1]
dx = np.convolve(im[0,:], Hx)
d2x = np.convolve(dx, Hx)
plt.figure(figsize=(12,15))
plt.subplot(311)
plt.title('Intensity')
plt.plot(im[0,:],'r')
plt.xticks([])
#plt.xlabel(r'Width ($x$)')
plt.ylabel(r'$I(x)$')
#plt.plot([100,100],[0,1],'k--')
plt.plot([150,150],[0,1],'k--')
plt.plot([50,50],[0,1],'k--')
plt.subplot(312)
#plt.xlabel(r'Width ($x$)')
plt.ylabel(r"$I'(x)$")
plt.title('Derivative')
plt.plot(dx)
plt.xticks([])
#plt.plot([100,100],[-0.04,0.08],'k--')
plt.plot([50,50],[-0.06,0.13],'k--')
plt.plot([150,150],[-0.06,0.13],'k--')
plt.subplot(313)
plt.xlabel(r'Width ($x$)')
plt.ylabel(r"$I''(x)$")
plt.title('Second Derivative')
plt.plot(d2x)
#plt.plot([100,100],[-0.01,0.01],'k--')
plt.plot([50,50],[-0.03,0.03],'k--')
plt.plot([150,150],[-0.03,0.03],'k--')

Out[6]:
[<matplotlib.lines.Line2D at 0x1247c4710>]

Computing image derivatives¶

Image gradients capture local changes in the image. These are a fundamental image processing operation that is required for a number of subsequent image analysis procedures, such as edge detection, segmentation, feature construction, etc.

Option 1: Reconstruct a continuous function $f$, then compute partial derivatives as follows

$$\frac{\partial f(x,y)}{\partial x} = \lim_{\epsilon \to 0} \frac{f(x+\epsilon,y) - f(x,y)}{\epsilon}$$$$\frac{\partial f(x,y)}{\partial y} = \lim_{\epsilon \to 0} \frac{f(x,y+\epsilon) - f(x,y)}{\epsilon}$$

Option 2: We can use finite difference approximation to estimate image derivatives $$\frac{\partial f(x,y)}{\partial x} \approx \frac{f[x+1,y] - f[x,y]}{1}$$

$$\frac{\partial f(x,y)}{\partial y} \approx \frac{f[x,y+1] - f[x,y]}{1}$$

Observation: We can compute image derivatives using convolutions (linear filtering)

Kernel for computing image derivative w.r.t. $x$ (width)¶

$$H_x = \left[ \begin{array}{rr} 1 & -1 \\ \end{array} \right]$$

Kernel for computing image derivative w.r.t. $y$ (height)¶

$$H_y = \left[ \begin{array}{r} 1 \\ -1 \\ \end{array} \right]$$

Finite difference filters¶

• Sobel filters
$$H_x = \left[\begin{array}{rrr} -1 & 0 & 1 \\ -2 & 0 & 2 \\ -1 & 0 & 1 \end{array}\right] \mathrm{\ and\ } H_y = \left[\begin{array}{rrr} 1 & 2 & 1 \\ 0 & 0 & 0 \\ -1 & -2 & -1 \end{array}\right]$$
• Prewire
$$H_x = \left[\begin{array}{rrr} -1 & 0 & 1 \\ -1 & 0 & 1 \\ -1 & 0 & 1 \end{array}\right] \mathrm{\ and\ } H_y = \left[\begin{array}{rrr} 1 & 1 & 1 \\ 0 & 0 & 0 \\ -1 & -1 & -1 \end{array}\right]$$
• Roberts
$$H_x = \left[\begin{array}{rr} 0 & 1 \\ -1 & 0 \end{array}\right] \mathrm{\ and\ } H_y = \left[\begin{array}{rr} 1 & 0 \\ 0 & -1 \end{array}\right]$$

Image derviatives highlights edge pixels¶

• X derivative is computed by convolving image with $H_x$, and Y derivative is computed by convolving image with $H_y$.
• Pixels sitting on vertical edges are highlighted in X derivative.
• Pixels sitting on horizontal edges are highlighted in Y derivative.

• Image gradient $\nabla I$ points to the direction of most rapid change in intensity.
$$\nabla I = \left[ \begin{array}{rr}\frac{\partial I}{\partial x} & \frac{\partial I}{\partial y} \end{array}\right]^T$$
• Gradient magnitude (can be used to identify edge pixels)
$$\| \nabla I \| = \sqrt{ \left( \frac{\partial I}{\partial x} \right)^2 + \left( \frac{\partial I}{\partial y} \right)^2 }$$
• Gradient direction (is perpendicular to the underlying edge)
$$\theta = \tan^{-1} \left( \frac{\partial I}{\partial y} / \frac{\partial I}{\partial x} \right)$$

Image Derivatives are highly sensitive to noise¶

Use smoothing first to get rid of the high-frequency components as seen below

Code Example¶

In [7]:
import cv2
import numpy as np
import scipy as sp
from scipy import signal
import matplotlib.pyplot as plt
%matplotlib inline


In [8]:
# flipped since convolution is flips

Hx = np.array([[1,-1]], dtype='float32')
Hy = np.array([[1],[-1]], dtype='float32')

plt.imshow(Hx, cmap='gray', interpolation='none')
plt.xticks([])
plt.yticks([])

Out[8]:
([], <a list of 0 Text major ticklabel objects>)
In [9]:
plt.imshow(Hy, cmap='gray', interpolation='none')
plt.xticks([])
plt.yticks([])

Out[9]:
([], <a list of 0 Text major ticklabel objects>)

Square¶

In [10]:
square = np.zeros((32,32), dtype='float32')
square[8:8+16, 8:8+16] = 1.0

plt.imshow(square, cmap='gray', interpolation='none')
plt.xticks([])
plt.yticks([])

Out[10]:
([], <a list of 0 Text major ticklabel objects>)
In [11]:
Hx_square = sp.signal.convolve2d(square, Hx, 'same')

plt.imshow(Hx_square, cmap='gray', interpolation='none')
plt.xticks([])
plt.yticks([])

Out[11]:
([], <a list of 0 Text major ticklabel objects>)
In [12]:
Hy_square = sp.signal.convolve2d(square, Hy, 'same')

plt.imshow(Hy_square, cmap='gray', interpolation='none')
plt.xticks([])
plt.yticks([])

Out[12]:
([], <a list of 0 Text major ticklabel objects>)

In [13]:
plt.figure(figsize=(15,5))
plt.subplot(131)
plt.imshow(square, cmap='gray', interpolation='none')
plt.xticks([])
plt.yticks([])
plt.title('Square')
plt.subplot(132)
plt.imshow(Hx_square, cmap='gray', interpolation='none')
plt.xticks([])
plt.yticks([])
plt.title('X derivative')
plt.subplot(133)
plt.imshow(Hy_square, cmap='gray', interpolation='none')
plt.xticks([])
plt.yticks([])
plt.title('Y derivative')

Out[13]:
Text(0.5, 1.0, 'Y derivative')
In [14]:
square_grad_mag = np.sqrt(np.square(Hx_square) + np.square(Hy_square))
plt.figure(figsize=(10,10))

plt.xticks([])
plt.yticks([])

Out[14]:
Text(0.5, 1.0, 'Gradient magnitude')

In [15]:
square_grad_angle = np.arctan2(Hy_square, Hx_square)

plt.figure(figsize=(10,10))
x_c, y_c = np.linspace(0,31,32), np.linspace(0,31,32)
xx, yy = np.meshgrid(x_c, y_c)
plt.quiver(xx, yy, Hx_square, Hy_square, scale=25, color='red', width=.005)
plt.imshow(square, cmap='gray')
plt.xticks([])
plt.yticks([])

Out[15]:
([], <a list of 0 Text major ticklabel objects>)
In [16]:
plt.figure(figsize=(20,10))
plt.subplot(121)
plt.xticks([])
plt.yticks([])
plt.subplot(122)
plt.quiver(xx, yy, Hx_square, Hy_square, scale=25, color='red', width=.005)
plt.imshow(square, cmap='gray')
plt.xticks([])
plt.yticks([])

Out[16]:
Text(0.5, 1.0, 'Gradient')

Circle¶

In [17]:
s = 128
circle = np.zeros((s,s), dtype='float32')

import math
cx, cy = 128./2, 128./2
for i in range(0, s-1):
for j in range(0, s-1):
if math.sqrt( (cx - i)**2 + (cy - j)**2 ) < 128./4:
circle[i,j] = 1.0

plt.imshow(circle, cmap='gray', interpolation='none')
plt.xticks([])
plt.yticks([])

Out[17]:
([], <a list of 0 Text major ticklabel objects>)
In [18]:
Hx = np.array([[1,0,-1],[2,0,-2],[1,0,-1]])

Hx_circle = sp.signal.convolve2d(circle, Hx, 'same')

plt.imshow(Hx_circle, cmap='gray', interpolation='none')
plt.xticks([])
plt.yticks([])

Out[18]:
([], <a list of 0 Text major ticklabel objects>)
In [19]:
Hy = np.array([[1,2,1],[0,0,0],[-1,-2,-1]])

Hy_circle = sp.signal.convolve2d(circle, Hy, 'same')

plt.imshow(Hy_circle, cmap='gray', interpolation='none')
plt.xticks([])
plt.yticks([])

Out[19]:
([], <a list of 0 Text major ticklabel objects>)
In [20]:
plt.figure(figsize=(15,5))
plt.subplot(131)
plt.imshow(circle, cmap='gray', interpolation='none')
plt.xticks([])
plt.yticks([])
plt.title('Circle')
plt.subplot(132)
plt.imshow(Hx_circle, cmap='gray', interpolation='none')
plt.xticks([])
plt.yticks([])
plt.title('X derivative')

plt.subplot(133)
plt.imshow(Hy_circle, cmap='gray', interpolation='none')
plt.xticks([])
plt.yticks([])
plt.title('Y derivative')

Out[20]:
Text(0.5, 1.0, 'Y derivative')

In [21]:
circle_grad_mag = np.sqrt(np.square(Hx_circle) + np.square(Hy_circle))
plt.figure(figsize=(10,10))
plt.xticks([])
plt.yticks([])

Out[21]:
Text(0.5, 1.0, 'Gradient magnitude')

In [22]:
circle_grad_angle = np.arctan2(Hy_circle, Hx_circle)

plt.figure(figsize=(10,10))
x_c, y_c = np.linspace(0,127,128), np.linspace(0,127,128)
xx, yy = np.meshgrid(x_c, y_c)
plt.quiver(xx, yy, Hx_circle, Hy_circle, scale=500, color='red')
#plt.imshow(circle, cmap='gray')
plt.xticks([])
plt.yticks([])

Out[22]:
([], <a list of 0 Text major ticklabel objects>)
In [23]:
plt.figure(figsize=(30,10))
plt.subplot(131)
plt.xticks([])
plt.yticks([])

plt.plot([64,104],[64,64], 'r')
plt.plot([64,104],[64-40,64-40], 'r')
plt.plot([64,64],[64,64-40], 'r')
plt.plot([104,104],[64,64-40], 'r')

plt.subplot(132)
# x_c, y_c = np.linspace(0,127,128), np.linspace(0,127,128)
# xx, yy = np.meshgrid(x_c, y_c)
# plt.quiver(xx, yy, Hx_circle, Hy_circle, scale=500, color='red')
# plt.imshow(circle, cmap='gray')
# plt.xticks([])
# plt.yticks([])
# plt.subplot(133)
x_c, y_c = np.linspace(64,104,40), np.linspace(64,104,40)
xx, yy = np.meshgrid(x_c, y_c)
plt.quiver(xx, yy, Hx_circle[64:104,64:104], Hy_circle[64:104,64:104], scale=150, color='red')
#plt.imshow(circle[64:104,64:104], cmap='gray')
plt.xticks([])
plt.yticks([])

Out[23]:
([], <a list of 0 Text major ticklabel objects>)

Using color to visualize image gradients¶

In [24]:
hsv = np.zeros([circle.shape[0],circle.shape[0],3], dtype=np.uint8)
hsv[..., 1] = 255

mag, ang = cv2.cartToPolar(Hx_circle[...], Hy_circle[...])
# print(hsv.shape)
# print(mag.shape)
# print(ang.shape)
hsv[..., 0] = ang * 180 / np.pi / 2
hsv[..., 2] = cv2.normalize(mag, None, 0, 255, cv2.NORM_MINMAX)
bgr = cv2.cvtColor(hsv, cv2.COLOR_HSV2BGR)

plt.figure(figsize=(20,10))
plt.subplot(121)
plt.imshow(bgr)
plt.xticks([])
plt.yticks([])

plt.subplot(122)
width = 1000
x = np.linspace(width, -width, width*2+1)
y = np.linspace(width, -width, width*2+1)
xx, yy = np.meshgrid(x, y)

hsv = np.zeros([x.shape[0], y.shape[0], 3], dtype='uint8')
hsv[..., 1] = 255

mag, ang = cv2.cartToPolar(xx[...], yy[...])
# print(hsv.shape)
# print(mag.shape)
# print(ang.shape)
hsv[..., 0] = ang * 180 / np.pi / 2
hsv[..., 2] = cv2.normalize(mag, None, 0, 255, cv2.NORM_MINMAX)
bgr = cv2.cvtColor(hsv, cv2.COLOR_HSV2BGR)
plt.imshow(bgr)
plt.xticks([])
plt.yticks([])

Out[24]:
([], <a list of 0 Text major ticklabel objects>)

Effects of noise - code¶

In [25]:
# 1D kernel for computing derivatives
Hx = np.array([1,-1], dtype='float32')

In [26]:
np.random.seed(0)

x = np.linspace(0,100,101)
y = np.ones(x.shape)*10
y[x <= 50] = 0
plt.figure(figsize=(20,5))
plt.plot(x, y, 'r.')
plt.plot(x, y, 'b')
plt.xticks([])
plt.yticks([])

Out[26]:
([], <a list of 0 Text major ticklabel objects>)
In [27]:
plt.figure(figsize=(20,5))
dy = np.convolve(y, Hx, 'same')
plt.plot(x, dy, 'b')
plt.xticks([])
plt.yticks([])

Out[27]:
([], <a list of 0 Text major ticklabel objects>)

Creating a noisy signal¶

In [28]:
x = np.linspace(0,1000,1001)
yn = np.ones(x.shape)*9
yn[x <= 500] = 0
yn = yn + np.random.rand(yn.shape[0])*8

plt.figure(figsize=(20,5))
plt.plot(x, yn, 'r.')
plt.plot(x, yn, 'b')
plt.xticks([])
plt.yticks([])

Out[28]:
([], <a list of 0 Text major ticklabel objects>)