Image Derivatives

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

© Faisal Qureshi

License

Creative Commons Licence
This work is licensed under a Creative Commons Attribution-NonCommercial 4.0 International License.

Lesson Plan

  • Why do we care about image gradients?
  • Computing image derivatives
  • Sobel filters
  • Gradient magnitude and directions
  • Visualizing image gradients

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

  • 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

Speeding up

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

Kernel for computing image gradients

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>)

Gradient magnitude

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.imshow(square_grad_mag, cmap='gray', interpolation='none')

plt.xticks([])
plt.yticks([])
plt.title('Gradient magnitude')
Out[14]:
Text(0.5, 1.0, 'Gradient magnitude')

Gradient direction

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.title('Gradient')
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.imshow(square_grad_mag, cmap='gray', interpolation='none')
plt.xticks([])
plt.yticks([])
plt.title('Gradient magnitude')
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([])
plt.title('Gradient')
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')

Gradient magnitude

In [21]:
circle_grad_mag = np.sqrt(np.square(Hx_circle) + np.square(Hy_circle))
plt.figure(figsize=(10,10))
plt.imshow(circle_grad_mag, cmap='gray', interpolation='none')
plt.xticks([])
plt.yticks([])
plt.title('Gradient magnitude')
Out[21]:
Text(0.5, 1.0, 'Gradient magnitude')

Gradient direction

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.title('Gradient')
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.imshow(circle_grad_mag, cmap='gray', interpolation='none')
plt.xticks([])
plt.yticks([])
plt.title('Gradient magnitude')

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.title('Gradient')
# 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.title('Gradient')
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>)