Image Derivatives

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

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

Image edges and image deriatives

import numpy as np
import matplotlib.pyplot as plt

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

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

[<matplotlib.lines.Line2D at 0x1147aaad0>]

# 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$)')

Text(0, 0.5, 'Height ($y$)')

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

Text(0, 0.5, '$I(x)$')

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

[<matplotlib.lines.Line2D at 0x115b20550>]

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

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

# 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([])

([], [])

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

([], [])

Square

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([])

([], [])

Hx_square = sp.signal.convolve2d(square, Hx, 'same')

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

([], [])

Hy_square = sp.signal.convolve2d(square, Hy, 'same')

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

([], [])

Gradient magnitude

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

Text(0.5, 1.0, 'Y derivative')

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

Text(0.5, 1.0, 'Gradient magnitude')

Gradient direction

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([])

([], [])

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

Circle

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([]);

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([]);

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([]);

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

Gradient magnitude

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

Gradient direction

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([]);

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([]);

Using color to visualize image gradients

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([]);

Effects of noise - code

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

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([]);

plt.figure(figsize=(20,5))
dy = np.convolve(y, Hx, 'same')
plt.plot(x, dy, 'b')
plt.xticks([])
plt.yticks([])

([], [])

Creating a noisy signal

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([]);

Edge finding in noisy signal

dyn = np.convolve(yn, Hx, 'same')
plt.figure(figsize=(20,5))
plt.plot(x, dyn, 'b')
plt.xticks([])
plt.yticks([]);

1D Gaussian - recap

s = 4
ticks = np.linspace(-s,s,2*s+1)

gxx = np.linspace(-40,40,801)
mu = 0
sigma = 5
g = np.exp(- (gxx-mu)**2 / (2*(sigma**2)) ) / ( 2 * np.pi * sigma)

plt.figure(figsize=(20,5))
plt.plot(gxx, g, 'k', linewidth=5)
plt.xlim([-s*sigma, s*sigma])
plt.yticks([])
plt.title(r'Gaussian.  ($\mu=0, \sigma=5$)');

1D Gaussian - First Derivative

dg = np.convolve(g, Hx, 'same')

plt.figure(figsize=(20,5))
plt.plot(gxx, dg, 'k', linewidth=5)
plt.xlim([-s*sigma, s*sigma])
plt.yticks([])
plt.title(r'Gaussian 1st derivative.  ($\mu=0, \sigma=5$)');

1D Guassian - Second Derviative

ddg = np.convolve(dg, Hx, 'same')

plt.figure(figsize=(20,5))
plt.plot(gxx, ddg, 'k', linewidth=5)
plt.xlim([-s*sigma, s*sigma])
plt.yticks([])
plt.title(r'Gaussian 2nd derivative.  ($\mu=0, \sigma=5$)');

Setting up Gaussian and Gaussian derivative

gxx = np.linspace(-20,20,41)
mu = 0
sigma = 7
g = np.exp(- (gxx-mu)**2 / (2*(sigma**2)) ) / ( 2 * np.pi * sigma)
plt.figure(figsize=(20,5))
plt.plot(gxx, g)
plt.xticks([])
plt.yticks([])
plt.title('Gaussian');

dg = np.convolve(g, Hx, 'same')
plt.figure(figsize=(20,5))
plt.plot(gxx, dg)
plt.xticks([])
plt.yticks([])
plt.title('First derivative of Gaussian');

Gaussian smoothing

gyn = np.convolve(yn, g, 'same')
plt.figure(figsize=(20,5))
plt.plot(x[5:-15], gyn[5:-15])
plt.plot(x[5:-15], gyn[5:-15],'r.')
plt.xticks([])
plt.yticks([])
plt.title('Signal smoothed with Gaussian');

Compute derivative of smoothed signal

dgyn = np.convolve(gyn, Hx, 'same')
plt.figure(figsize=(20,5))
plt.plot(x[5:-15], dgyn[5:-15], 'b')
plt.xticks([])
plt.yticks([])
plt.title('Derivative of smoothed signal');

Smooth signal with derivative of Gaussian

Convolve noisy signal with derivative of Gaussian (saving one convolution operation)

dgyn2 = np.convolve(dg, yn, 'same')
plt.figure(figsize=(20,5))
plt.plot(x[5:-15], dgyn2[5:-15], 'b')
plt.xticks([])
plt.yticks([])
plt.title('Signal convolution with first derivative of Gaussian');