Feature Detection, Part 2: Laplace & Gaussian Operators

Contents

Feature detection is a domain of computer vision that focuses on using tools to detect regions of interest in images. A significant aspect of most feature detection algorithms is that they do not employ machine learning under the hood, making the results more interpretable and even faster in some cases. In the previous article of this series, we learned about edges — zones where the image intensity drastically changes. Edges usually represent local zones where an object of interest can be found. For example, consider an image of a blue sky with a small plane flying over it. The intensity difference across the entire image will be nearly zero, except for the small zone containing the intersection of the plane with the blue sky background. As a result, we can easily detect the zone where the plane is located and its shape using the image derivative or the Sobel filter. Edges are located in zones of substantial pixel intensity change. In this example, edges can be found near the border of the plane and clouds. In other words, we can detect edges by taking the minimum/maximum peaks of the first derivative of the intensity function, which is represented by the image in the middle. Second image derivative Definition Discrete approximation Diagonal edge detection Isotropy Noise Gaussian filter Derivative of Gaussian Laplacian of Gaussian OpenCV Laplacian of Gaussian Zero-crossing detection Visualization Conclusion Resources

Second image derivative

Let us return to the mentioned example and ask ourselves what will happen if we take the second derivative? The answer is shown on the right graph below. In this case, the minimum/maximum peaks of the first derivative will correspond to zeros of the second derivative.

Intensity values across the X-axis (image on the left).
First derivative (image in the middle): local extrema indicate an edge.
Second derivative (image on the right): zero-crossing area indicates an edge.
Image adapted by the author. Source: Laplace Operator (OpenCV documentation).

This means that we can now use zero values of the second derivative as an additional criterion for edge detection. However, it is essential to remain cautious, as other points on the image can also result in second derivatives being equal to 0, yet not be edges themselves. For that, it is recommended to use other filters to eliminate such cases.

Generally, image edges result in very sharp zero-crossing areas.

Definition

The Laplacian is formally defined by the formula below:

As we can see, the Laplacian is simply the sum of second derivatives with respect to x and y directions.

The Laplacian does not provide information about edge direction.

Discrete approximation

Given the Laplacian formula above for the continuous case, let us try to obtain an approximation for the discrete case of images.

For that, let us suppose that for the image part below, we want to calculate the Laplacian for the pixel I₂₂. Let us compute a derivative in the X-axis direction.

In this example, we would like to approximate Laplace for the central pixel I₂₂

For that, we reuse twice the derivative formula from the previous article by choosing two values for Δx: -1 and 1. We get:

For Δx = -1: dI / dx = I₂₃ – I₂₂
For Δx = 1: dI / dx = I₂₂ – I₂₁

To calculate the first derivative, we basically calculate the difference between adjacent pixels. We can proceed with the same logic for the second derivative, which can be informally thought of as the difference of the difference.

Therefore, we can take the two differences we have just calculated and find their difference. More formally, it would be the same if we applied the standard derivative formula by setting Δx to -1. So we get:

d²I / dx = (I₂₃ – I₂₂) – (I₂₂ – I₂₁) = I₂₃ – 2I₂₂ + I₂₁

Great! We have computed the second derivative in the X-direction. We can do the analogous procedure for the second derivative in the Y-direction. We will get the following expression:

d²I / dy = (I₃₂ – I₂₂) – (I₂₂ – I₁₂) = I₃₂ – 2I₂₂ + I₁₂

Finally, the only thing left to do is to sum up both derivatives. We get:

d²I / dx + d²I / dy = I₁₂ + I₂₁ + I₂₃ + I₃₂ – 4I₂₂

The whole Laplacian computation process can be visualised by the diagram below:

Obtaining a formula for the Laplace calculation in the discrete case

Based on the achieved result, we can also construct a 3×3 convolutional kernel for the Laplacian:

As an interesting fact, the sum of the elements of this kernel is 0, which means that if an input image is constant, then applying the filter will produce a matrix with zero elements. It is a logical result since there is no change in intensity.

In comparison to Sobel and Scharr kernels, the Laplace kernel detects intensity changes in both directions. It is sufficient to apply the 3×3 kernel to any image, and the Laplace operator will output the final scalar values representing intensity changes.

As a reminder, with Sobel and Scharr operators, the kernels were applied separately to the X and Y axes, and then their magnitudes were calculated.

Diagonal edge detection

We have deducted a kernel that represents a sum of second derivatives across the X and Y axes. Therefore, this method would be well-suited for detecting both horizontal and vertical edges. But what if an edge in the image has a diagonal direction? We have not taken it into account yet.

For that reason, the kernel above is usually slightly adjusted to account for the diagonal direction. One of the most popular options is to use the following kernel:

An alternative to Laplace kernel which is symmetric horizontally, vertically and diagonally.

Isotropy

We can observe that the Laplace matrix is symmetric, which makes it isotropic. Isotropy is a property according to which the kernel is rotation-invariant, meaning that the output resulting from applying an isotropic filter to an image and its rotated version is the same.

Noise

The Laplace kernel we saw above works well for edge detection. However, we didn’t take into account another aspect that would prevent us from efficiently applying that filter in real life: noise.

We saw several graphs of image intensity change along the X-axis at the beginning of this article, where we plotted the first and second derivatives of the image intensity. In fact, these plots were constructed for a perfect image where there is no noise.

If the noise is present, we might end up in a situation like the one depicted in the diagram below.

Given the fact that noise is present on the input image (on the left), intensity values along the X-axis result in oscillations making the first derivative (on the right) (and as a result, the second one) hard to analyse.

The graph on the left illustrates a more realistic scenario of intensity values along an axis in the image. Although the noise is not very strong, it still fluctuates rapidly in local areas. At the same time, we use derivatives to detect the same rapid changes, which creates a problem.

Ultimately, taking the first derivative would result in the change in intensity, as shown on the graph on the right. Obviously, normal edge detection using derivatives is not possible in the presence of such oscillations.

Gaussian filter

A Gaussian filter is a method for suppressing noise in an image, allowing the Laplacian or another edge detection operator to be applied without constraints.

Formally, the Gaussian distribution is given by the formula below (where μ = 0):

Gaussian formula. μ is the mean (in our case μ = 0), and σ represents the standard deviation, which can be tuned for a given problem.

The Gaussian function g(x) is plotted below on the top-right graph:

Top-left: Intensity function.
Top-right: Gaussian function.
Bottom-left: the product between intensity and the Gaussian functions.
Bottom-right: derivative of a product.

Multiplying the intensity function I(x) by the Gaussian g(x) tends to smooth the intensity in general, which results in easier analysis. The example is shown in the bottom-left diagram.

Given the smooth function form, we can then take the first derivative, whose extrema points will indicate edges (bottom-right diagram).

Derivative of Gaussian

We can notice that that taking the first image derivative is a linear operation, thus we can decompose our computations in the following manner:

Instead of calculating the derivative of a product, we can precompute the derivative of the Gaussian and then multiply it by the intensity function; the result will be the same but will require less computations.

This allows us to pre-compute the first derivative of the Gaussian in advance and then multiply it by the intensity function, which optimizes calculations.

Laplacian of Gaussian

We can apply the same trick to the second derivative, about which we discussed at the beginning of this article, for edge detection. As a result, we can compute the second Gaussian derivative in advance and then multiply it by the intensity function. The result is shown on the bottom-right graph:

Top-left: Intensity function.
Top-right: Gaussian function.
Bottom-left: the second derivative of the Gaussian function (also known as the inverted Mexican sombrero).
Bottom-right: the product between the second Gaussian derivative and the Intensity function.

The resulting function is called the Laplacian of Gaussian and is used extensively in edge detection, while also being resilient to image noise.

The second derivative of a Gaussian (bottom-left image) is often referred to as an inverted “Mexican sombrero” because of its high shape similarity with a sombrero hat.

OpenCV

Having studied edge detection theory using second image derivatives, it is time to look at how Laplace filters can be implemented in OpenCV.

Laplacian of Gaussian

First of all, let us import the necessary libraries and upload an input image:

import cv2
import numpy as np
import matplotlib.pyplot as plt

We are going to use an image example containing several coins:

Let us read the image:

image = cv2.imread('data/input/coins.png')
image = cv2.cvtColor(image, cv2.COLOR_BGR2GRAY)

We have converted the input image to grasycale format to be able to later take derivatives relative to intensity changes.

As discussed above, before using the Laplace operator, we are going to apply the Gaussian filter:

image = cv2.GaussianBlur(image, (7, 7), 0)

The second parameter (7, 7) refers to the kernel size.
The third parameter is used to define the value of standard deviation σ from the Gaussian formula given above. Here, we set it to 0, which means that the standard deviation σ will be chosen automatically by OpenCV based on the provided kernel size.

Now we are ready to apply the Laplace filter:

laplacian_signed = cv2.Laplacian(image, cv2.CV_16S, ksize=3)

Here we specify the output type as cv2.CV_16S, which is equivalent to the short type, with values in the range [-32768, 32767]. We do that since the Laplacian filter can produce values outside the standard pixel interval [0, 255]. If we had not done this, then the output values would have been clipped to [0, 255] and we would have lost information.

OpenCV also provides a useful function to map the raw Laplacian output (or any other resulting values) to the standard range [0, 255]:

laplacian_absolute = cv2.convertScaleAbs(laplacian_signed)

The resulting variable laplacian_absolute has the same shape as laplacian_signed but with values clipped to [0, 255]. The mapping is done using the function cv2.convertScaleAbs() in a way that preserves information about zero-crossings:

The values whose absolute value is greater than 255 are clipped to 255.
For the values between -255 and 255, the absolute value is taken.

cv2.convertScaleAbs() function maps pixel values to the [0, 255] interval.

Zero-crossing detection

Unfortunately, the OpenCV documentation only shows an example of Laplacian application but does not demonstrate how to use its results for zero-cross edge detection. Below is a simple snippet example that implements that:

laplacian_sign = np.sign(laplacian_signed)

zero_crossings = np.zeros_like(laplacian_sign, dtype=bool)
for shift in [-5, 5]:
    zero_crossings |= (np.roll(laplacian_sign, shift, axis=0) * laplacian_sign < 0)
    zero_crossings |= (np.roll(laplacian_sign, shift, axis=1) * laplacian_sign < 0)

threshold = 20
edges = np.uint8(zero_crossings & (np.abs(laplacian_signed) > threshold)) * 255

In simple terms, we create a zero_crossings variable that contains information on whether a pixel at position (x, y) is a zero-crossing candidate. In our code, we consider a pixel a zero-crossing if its sign (+ or -) differs from the sign of any shifted pixel by five positions in either horizontal or vertical direction.

We could have chosen another shift constant (not necessarily 5) or even combine several of them, and/or also take into account diagonal shifts as well.

To exclude non-relevant zero-crossings, we then only keep those whose absolute laplacian value is greater than a defined threshold (20).

Visualization

Finally, let us plot the images we have retrieved throughout our analysis:

plt.figure(figsize=(12, 4))

plt.subplot(1, 3, 1)
im1 = plt.imshow(laplacian_signed, cmap='gray', vmin=laplacian_signed.min(), vmax=laplacian_signed.max())
plt.title("Laplacian (signed)")
plt.axis('off')
plt.colorbar(im1, fraction=0.05, pad=0.05, label='Pixel value')

plt.subplot(1, 3, 2)
im2 = plt.imshow(laplacian_absolute, cmap='gray', vmin=laplacian_absolute.min(), vmax=laplacian_absolute.max())
plt.title("Laplacian (absolute)")
plt.axis('off')
plt.colorbar(im2, fraction=0.05, pad=0.05, label='Pixel value')

plt.subplot(1, 3, 3)
im2 = plt.imshow(edges, cmap='gray', vmin=edges.min(), vmax=edges.max())
plt.title("Edges from zero-crossings")
plt.axis('off')
plt.colorbar(im2, fraction=0.05, pad=0.05, label='Pixel value')

plt.tight_layout()
plt.show()

Here is the output result we get:

Output after applying Laplace filters. The right image is binary (containing only pixels with values of 0 or 255).

As you can see, we successfully detected edges on the right binary image. The next step could, for example, involve applying the OpenCV method cv2.findContours() to detect contours and further analyze them.

Given the simple form of the signed Laplacian output, we could have also used it to detect the borders of coins.

Conclusion

Using knowledge about image derivatives from the previous part, this article analyzes the role of the second derivative and zero-cross detection. It was also a great way to learn about the Laplace and Gaussian filters, which can be very easily implemented with OpenCV for edge detection.

Resources

All images unless otherwise noted are by the author.