Laplacian of Gaussian

Jacky Baltes
National Taiwan Normal University
Taipei, Taiwan
jacky.baltes@ntnu.edu.tw

25 October 2022

import matplotlib.pyplot as plt
import numpy as np
import math
import random

import cv2
from skimage import io

url = "https://upload.wikimedia.org/wikipedia/commons/thumb/8/86/Red-tailed_Tropicbird3.jpg/200px-Red-tailed_Tropicbird3.jpg" #@param {type:"string"}
img = io.imread( url )
img = cv2.resize(img, (320, 240 ) )

gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)

grad_x = cv2.Sobel( gray, cv2.CV_16S, 1, 0, 3 )

grad_y = cv2.Sobel( gray, cv2.CV_16S, 0, 1, 3 )

abs_grad_x = cv2.convertScaleAbs( grad_x )
#print( grad_x[0:5,0:5])
#print( abs_grad_x[0:5,0:5])

abs_grad_y = cv2.convertScaleAbs( grad_y )

fig = plt.figure( figsize=(10,10) )

ax1 = fig.add_subplot(1,2,1)
ax1.set_title('Sobel X')
ax1.imshow(abs_grad_x, interpolation = 'bicubic', cmap = "gray")
ax1.set_xticks([])
ax1.set_yticks([])  

ax2 = fig.add_subplot(1,2,2)
ax2.set_title('Sobel Y')
ax2.imshow(abs_grad_y, interpolation = 'bicubic', cmap = "gray")
ax2.set_xticks([])
ax2.set_yticks([])  

bird21 = addJBFigure("bird21", 0, 0, fig )
plt.close()

Sobel Edge Detection

Sobel edge detector uses first derivative to find edges (i.e., large changes in brightness to the neighboring pixels)

grad = cv2.addWeighted( abs_grad_x, 0.5, abs_grad_y, 0.5, 0 )

ret, threshold = cv2.threshold( grad, 100, 255, cv2.THRESH_BINARY )
#print( threshold[0:5, 0:5] )

ax1 = fig.add_subplot(1,2,1)
ax1.set_title('Threshold')
ax1.imshow(threshold, interpolation = 'bicubic', cmap = "gray")
ax1.set_xticks([])
ax1.set_yticks([])  

ax2 = fig.add_subplot(1,2,2)
ax2.set_title('Threshold (Zoom)')
ax2.imshow(threshold[100:140,150:190], interpolation = 'bicubic', cmap = "gray")
ax2.set_xticks([])
ax2.set_yticks([])  

b23 = addJBFigure("b23", 0, 0, fig )
plt.close()

grad = cv2.addWeighted( abs_grad_x, 0.5, abs_grad_y, 0.5, 0 )

ret, threshold = cv2.threshold( grad, 100, 255, cv2.THRESH_BINARY )
#print( threshold[0:5, 0:5] )

ax1 = fig.add_subplot(1,2,1)
ax1.set_title('Threshold')
ax1.imshow(threshold, interpolation = 'bicubic', cmap = "gray")
ax1.set_xticks([])
ax1.set_yticks([])  

ax2 = fig.add_subplot(1,2,2)
ax2.set_title('Threshold (Zoom)')
ax2.imshow(threshold[100:140,150:190], interpolation = 'bicubic', cmap = "gray")
ax2.set_xticks([])
ax2.set_yticks([])  

b23 = addJBFigure("b23", 0, 0, fig )
plt.close()

Thresholding

Threshold the absolute gradient image to classify pixels as edge/non-edge pixels

Problems with Sobel Edge Detection

In the zoomed image, we can see that there are multiple pixels per edge. Should be just a single pixel

Some spurious edge pixels in the background and the bird

First Order Operators and Edges

Sobel edge detection uses the first derivative to detect an edge (first order operator)

Two adjacent edges may have a high first derivative, hence we get multiple edge pixels for one edge

The second derivative allows us to distinguish between a maximum, minimum, or saddle point

These are called second order operators

Laplacian of Gaussian (LoG)

Check exact point where 2nd derivative is 0. May not lie on an exact pixel

Look for zero crossing of the 2nd derivative

url = "https://upload.wikimedia.org/wikipedia/commons/thumb/8/86/Red-tailed_Tropicbird3.jpg/200px-Red-tailed_Tropicbird3.jpg" #@param {type:"string"}
img = io.imread( url )
img = cv2.resize(img, (320, 240 ) )

fig = plt.figure( figsize=(10,10) )
ax1 = fig.add_subplot(1,1,1)
ax1.set_title('Original Image')

ax1.imshow(img, interpolation = 'bicubic')
ax1.set_xticks([])
ax1.set_yticks([])  
plt.show()

Derivation of the Laplacian Operator

\[ I(x,y) = \left[ \begin{array}{ccccc} p_{-2,-2}, p_{-1,-2}, p_{0,-2}, p_{1,-2}, p_{2,-2} \\ p_{-2,-1}, p_{-1,-1}, p_{0,-1}, p_{1,-1}, p_{2,-1} \\ p_{-2,0}, p_{-1,0}, \underline{p_{0,0}}, p_{1,0}, p_{2,0} \\ p_{-2,1}, p_{-1,1}, p_{0,1}, p_{1,1}, p_{2,1} \\ p_{-2,2}, p_{-1,2}, p_{0,2}, p_{1,2}, p_{2,2} \\ \end{array} \right] \]

Calculating the First Derivative

First derivative is the change in value of the image function \(I(x,y)\)

Since we are dealing with samples measured along the x and y direction

Abbreviated as \( I_{x}(x,y), I_{y}(x,y) \)

\[ I_x(x,y) = \frac{p_{1,0} - p_{-1,0}}{2}\\ I_y(x,y) = \frac{p_{0,1} - p_{0,1}}{2} \]

Calculating the Second Derivative

First derivative is the change in the first derivative, I need to calculate two first derivatives

\[ I^{+}_{x}(x,y) = p_{1,0} - p_{0,0} \\ I^{-}_{x}(x,y) = p_{0,0} - p_{-1,0} \\ \]

\[ \begin{array}{c} I^{+}_{x}(x,y) - I^{-}_{x}(x,y)\\ ( p_{1,0} - p_{0,0} ) - ( p_{0,0} - p_{-1,0} ) \\ p_{1,0} - 2 p_{0,0} + p_{-1,0} \\ \end{array} \]

Convolution mask for x and y direction simulataneously \[ \begin{array}{ccc} 0 & 1 & 0 \\ 1 & -4 & 1 \\ 0 & 1 & 0 \\ \end{array} \]

Laplacian Operator

Usually, also check for diagonal edges

\[ \begin{array}{ccc} 1 & 1 & 1 \\ 1 & -8 & 1 \\ 1 & 1 & 1 \\ \end{array} \]

kernel = 7 #@param ["3", "5", "7", "9", "11"] {type:"raw"}
fig = plt.figure()
ax1 = fig.add_subplot(1,1,1)
ax1.set_title('Grayscale Image')

gray = cv2.cvtColor( img,  cv2.COLOR_RGB2GRAY );

ax1.imshow(gray, interpolation = 'bicubic', cmap='gray')
ax1.set_xticks([])
ax1.set_yticks([])  
bird24 = addJBFigure( "bird24", 0, 0, fig )
plt.close()

Grayscale Conversion using OpenCV

%matplotlib notebook
laplacian = cv2.Laplacian(gray,cv2.CV_64F)
sobelX = cv2.Sobel(gray,cv2.CV_64F,1,0,ksize=3)  # x
sobelY = cv2.Sobel(gray,cv2.CV_64F,0,1,ksize=3)  # y

fig = plt.figure(figsize=(15,15))
ax1 = fig.add_subplot(1,2,1)
ax1.set_title('Laplacian')
#ax1.set_xticks([])
#ax1.set_yticks([])

ax1.imshow( laplacian, cmap = 'gray')

ax2 = fig.add_subplot(1, 2, 2)
ax2.set_title('SobelXY')
ax2.set_xticks([])
ax2.set_yticks([])
ax2.imshow( sobelX + sobelY, cmap = 'gray' )

f1 = addJBFigure( "f1", 0, 0, fig )

%matplotlib notebook
roi = gray[100:200,125:225]

laplacian = cv2.Laplacian(roi,cv2.CV_64F)
sobelX = cv2.Sobel(roi,cv2.CV_64F,1,0,ksize=5)  # x
sobelY = cv2.Sobel(roi,cv2.CV_64F,0,1,ksize=5)  # y

fig = plt.figure( figsize=(15,15) )
ax1 = fig.add_subplot(1,2,1)
ax1.set_title('Laplacian')
ax1.set_xticks([])
ax1.set_yticks([])

ax1.imshow( laplacian, cmap = 'gray')

ax2 = fig.add_subplot(1, 2, 2)
ax2.set_title('SobelXY')
ax2.set_xticks([])
ax2.set_yticks([])
ax2.imshow( sobelX + sobelY, cmap = 'gray' )

f2 = addJBFigure( "f2", 0, 0, fig )

gray2 = gray + 3 * np.random.normal( 0.0, 5, (gray.shape) )

fig = plt.figure( figsize=(15,15) )
ax1 = fig.add_subplot(1,2,1)
ax1.set_title('Original Image')

ax1.imshow(gray, interpolation = 'bicubic', cmap='gray')
ax1.set_xticks([])
ax1.set_yticks([])  

ax2 = fig.add_subplot(1,2,2)
ax2.set_title('Image distorted with noise')

ax2.imshow(gray2, interpolation = 'bicubic', cmap='gray')
ax2.set_xticks([])
ax2.set_yticks([])  


f10 = addJBFigure( "f10", 0, 0, fig )

roi = gray[100:200,125:225]
roi2 = gray2[100:200, 125:225]

fig = plt.figure( figsize=(15,15) )
ax1 = fig.add_subplot(1,2,1)
ax1.set_title('Original Image')

ax1.imshow(roi, interpolation = 'bicubic', cmap='gray')
ax1.set_xticks([])
ax1.set_yticks([])  

ax2 = fig.add_subplot(1,2,2)
ax2.set_title('Image distorted with noise')

ax2.imshow(roi2, interpolation = 'bicubic', cmap='gray')
ax2.set_xticks([])
ax2.set_yticks([])  

f11 = addJBFigure("f11", 0, 0, fig )

laplacian = cv2.Laplacian(roi2,cv2.CV_64F)
sobelX = cv2.Sobel(roi2,cv2.CV_64F,1,0,ksize=3)  # x
sobelY = cv2.Sobel(roi2,cv2.CV_64F,0,1,ksize=3)  # y

ret,thresh1 = cv2.threshold(laplacian,200,255,cv2.THRESH_BINARY)
ret,thresh2 = cv2.threshold(laplacian,60,255,cv2.THRESH_BINARY_INV)

fig = plt.figure( figsize=(15,15) )
ax1 = fig.add_subplot(1,2,1)
ax1.set_title('Laplacian')
ax1.set_xticks([])
ax1.set_yticks([])

ax1.imshow( thresh1 + thresh2, cmap = 'gray')

ax2 = fig.add_subplot(1, 2, 2)
ax2.set_title('SobelXY')
ax2.set_xticks([])
ax2.set_yticks([])
ax2.imshow( sobelX + sobelY, cmap = 'gray' )

f12 = addJBFigure("f12", 0, 0, fig )

Laplacian with Noise

Laplacian operator is very sensitive to noise

kernel = 9 #@param ["3", "5", "7", "9"] {type:"raw"}
roi2 = gray2[100:200, 125:225]

k = cv2.getGaussianKernel( kernel,0)

blur = cv2.GaussianBlur(roi2, (kernel,kernel), 0)

laplacian = cv2.Laplacian(blur,cv2.CV_64F)
sobelX = cv2.Sobel(blur,cv2.CV_64F,1,0,ksize=5)  # x
sobelY = cv2.Sobel(blur,cv2.CV_64F,0,1,ksize=5)  # y

ret,thresh1 = cv2.threshold(laplacian,8,255,cv2.THRESH_BINARY)
ret,thresh2 = cv2.threshold(laplacian,0,255,cv2.THRESH_BINARY_INV)

fig = plt.figure( figsize=(15,15) )
ax1 = fig.add_subplot(1,2,1)
ax1.set_title('Laplacian')
ax1.set_xticks([])
ax1.set_yticks([])

ax1.imshow( thresh1, cmap = 'gray')

ax2 = fig.add_subplot(1, 2, 2)
ax2.set_title('SobelXY')
ax2.set_xticks([])
ax2.set_yticks([])
ax2.imshow( sobelX + sobelY, cmap = 'gray' )

f13 = addJBFigure( "f13", 0, 0, fig )

blur = cv2.GaussianBlur(gray2, (kernel,kernel), 0)

laplacian = cv2.Laplacian(blur,cv2.CV_64F)
sobelX = cv2.Sobel(blur,cv2.CV_64F,1,0,ksize=5)  # x
sobelY = cv2.Sobel(blur,cv2.CV_64F,0,1,ksize=5)  # y

ret,thresh1 = cv2.threshold(laplacian,8,255,cv2.THRESH_BINARY)
ret,thresh2 = cv2.threshold(laplacian,0,255,cv2.THRESH_BINARY_INV)

fig = plt.figure( figsize=(15,15) )
ax1 = fig.add_subplot(1,2,1)
ax1.set_title('LoG')
ax1.set_xticks([])
ax1.set_yticks([])

ax1.imshow( thresh1, cmap = 'gray')

ax2 = fig.add_subplot(1, 2, 2)
ax2.set_title('SobelXY')
ax2.set_xticks([])
ax2.set_yticks([])
ax2.imshow( sobelX + sobelY, cmap = 'gray' )
f14 = addJBFigure("f14", 0, 0, fig )

blur = cv2.GaussianBlur(gray2, (kernel,kernel), 0)

laplacian = cv2.Laplacian(blur,cv2.CV_64F)
sobelX = cv2.Sobel(blur,cv2.CV_64F,1,0,ksize=5)  # x
sobelY = cv2.Sobel(blur,cv2.CV_64F,0,1,ksize=5)  # y

ret,thresh1 = cv2.threshold(laplacian,8,255,cv2.THRESH_BINARY)
ret,thresh2 = cv2.threshold(laplacian,0,255,cv2.THRESH_BINARY_INV)

fig = plt.figure( figsize=(15,15) )
ax1 = fig.add_subplot(1,2,1)
ax1.set_title('LoG')
ax1.set_xticks([])
ax1.set_yticks([])

ax1.imshow( thresh1, cmap = 'gray')

ax2 = fig.add_subplot(1, 2, 2)
ax2.set_title('SobelXY')
ax2.set_xticks([])
ax2.set_yticks([])
ax2.imshow( sobelX + sobelY, cmap = 'gray' )
f14 = addJBFigure("f14", 0, 0, fig )

Differences between Sobel and LoG

Best choice of edge detector depends on your application

Sobel (First order operators): + robust to noise, + complete outlines, - multiple pixels per edge, - extra edge pixels

Laplacian (Second order operators): + single pixel edges, - sensitive to noise (Gaussian blur), - holes in the outline

Differences between Sobel and LoG

Best choice of edge detector depends on your application

Sobel (First order operators): + robust to noise, + complete outlines, - multiple pixels per edge, - extra edge pixels

Laplacian (Second order operators): + single pixel edges, - sensitive to noise (Gaussian blur), - holes in the outline

Differences between Classification Algorithms

Two error classes: False positives, False negatives

In practical applications, the penalty for an error are not the same

Machine vision system to detect errors in a painting job, fit circles to rivet holes. False positives were a bigger problem, since they would need to check the part.

System to detect bad potato chips on an assembly line, then false negatives were a bigger problem. If you put a moist potato chip in the bag, then the customer will never buy your potato chips again

Error and Their Consequences

Assume that you want to buil a vision system to detect if someone is carrying a gun in carry-on luggage

Two types of errors, but false negatives may cause people to die in a terrorist attack

Even more, the priors are very different - millions of normal travellers, luckily very small number of terrorists. This causes a problem, since any algorithm that always says no gun will be 99.99..% accurate

Search everybody, but then that would be a huge invasion in privacy and costs a lot of time and money