Derivative Filters: Sobel, Laplacian & LoG

"I don't care what the pixel value is. I care where it's going. Show me change or show me zeros."

An Edge Detector Who Sees Things in Black and White

So far this chapter has produced images for humans: smoothed in Section 3.2, sharpened in Section 3.3. This section produces images for algorithms. A gradient map looks like a ghostly line drawing, not a photograph, but it answers the question downstream code actually asks: where does the scene change, how abruptly, and in what direction?

1. Derivatives on a Pixel Grid Intermediate

An image is a sampled function, so derivatives must become finite differences. For the horizontal partial derivative, three candidates present themselves: forward difference $I(x+1) - I(x)$, backward difference $I(x) - I(x-1)$, and the central difference:

$$ \frac{\partial I}{\partial x} \;\approx\; \frac{I(x+1, y) - I(x-1, y)}{2} $$

The central difference wins on accuracy (its error is second-order in the sample spacing, versus first-order for the others) and on symmetry (it assigns the derivative to the pixel itself, not to a half-pixel offset). As a kernel it is $\frac{1}{2}[\,-1\;\;0\;\;+1\,]$, an antisymmetric row whose weights sum to zero, the signature of every derivative filter: flat regions must map to exactly zero response.

Stacking both partials gives the gradient vector $\nabla I = (\partial I/\partial x,\; \partial I/\partial y)$, summarized by two scalar fields that downstream code consumes constantly:

$$ \text{magnitude:}\;\; \|\nabla I\| = \sqrt{I_x^2 + I_y^2} \qquad \text{orientation:}\;\; \theta = \operatorname{atan2}(I_y,\, I_x) $$

The magnitude is large wherever brightness changes quickly (edges), and the orientation points across the edge, perpendicular to its direction. These two numbers per pixel are the most reused intermediate quantity in classical vision.

2. Sobel, Prewitt, Scharr: Smoothed Derivatives Intermediate

The bare central difference has a fatal flaw: noise. Differentiation is a high-frequency amplifier (in Chapter 4's terms, its frequency response grows with frequency), and pixel noise is precisely high-frequency. Apply $[\,-1\;0\;+1\,]$ to a real photograph and the output is a blizzard of noise responses with edges barely visible inside it. The cure is the chapter's recurring move: smooth first. Smoothing perpendicular to the differentiation direction suppresses noise without diluting the derivative itself, and packaging both into one kernel via the separability logic of Section 3.6 gives the Sobel operator:

$$ S_x = \begin{bmatrix} -1 & 0 & +1 \\ -2 & 0 & +2 \\ -1 & 0 & +1 \end{bmatrix} = \begin{bmatrix} 1 \\ 2 \\ 1 \end{bmatrix} \begin{bmatrix} -1 & 0 & +1 \end{bmatrix}, \qquad S_y = S_x^\top $$

Read the factorization: a vertical $[1\;2\;1]$ smoothing column (a tiny Gaussian approximation) times the horizontal derivative row. Prewitt uses uniform $[1\;1\;1]$ smoothing instead; Scharr re-optimizes the weights to $[3\;10\;3]$ for markedly better rotational accuracy, which matters whenever the orientation $\theta$ feeds geometry. All three are the same idea with different smoothing budgets. For larger noise, OpenCV's ksize parameter grows the Sobel kernel ($5 \times 5$, $7 \times 7$), folding in more smoothing.

Figure 3.4.1 shows the geometry that makes first and second derivatives complementary instruments: across an edge, the first derivative peaks at the transition, while the second derivative crosses zero there, with lobes of opposite sign on either side.

In code, the cardinal rule is to compute derivatives in a signed float type. A Sobel response is negative on dark-to-light transitions read right-to-left, and storing it in uint8 silently destroys half the signal, the second classic dtype bug of this chapter after Section 3.3's wraparound.

import cv2
import numpy as np

gray = cv2.imread("staircase.jpg", cv2.IMREAD_GRAYSCALE)

# Signed float output (CV_32F) is essential: gradients go negative.
gx = cv2.Sobel(gray, cv2.CV_32F, 1, 0, ksize=3)   # d/dx
gy = cv2.Sobel(gray, cv2.CV_32F, 0, 1, ksize=3)   # d/dy

mag, ang = cv2.cartToPolar(gx, gy, angleInDegrees=True)

print(f"magnitude: max {mag.max():.0f}, mean {mag.mean():.1f}")
print(f"strong-edge pixels (mag > 100): {(mag > 100).mean():.1%}")
# Representative output:
# magnitude: max 1040, mean 21.3
# strong-edge pixels (mag > 100): 4.7%
# Typical scenes are mostly flat: a few percent of pixels carry the edges.

# For display only: rescale magnitude into [0, 255].
view = cv2.normalize(mag, None, 0, 255, cv2.NORM_MINMAX).astype(np.uint8)

What are gradient maps for? Thresholding the magnitude with the techniques of Chapter 2 yields a crude but serviceable edge map. Histograms of the orientation channel summarize the dominant directions in a region, an idea that matures into the HOG descriptors of classical recognition and the keypoint orientations of Chapter 10. The example below is among the simplest profitable uses: estimating the skew of a scanned document from its gradient orientations.

3. The Laplacian: Change of Change Intermediate

The Laplacian $\nabla^2 I = \partial^2 I/\partial x^2 + \partial^2 I/\partial y^2$ sums the unmixed second derivatives. Applying the finite-difference recipe twice per axis and adding gives the standard 4-neighbor kernel, with an 8-neighbor variant that includes diagonals:

$$ \nabla^2_{4} = \begin{bmatrix} 0 & 1 & 0 \\ 1 & -4 & 1 \\ 0 & 1 & 0 \end{bmatrix} \qquad \nabla^2_{8} = \begin{bmatrix} 1 & 1 & 1 \\ 1 & -8 & 1 \\ 1 & 1 & 1 \end{bmatrix} $$

Both are rotationally symmetric (unlike Sobel, the Laplacian has no preferred direction; it measures total curvature in one number) and both sum to zero. You met the 4-neighbor version inside Section 3.3's sharpening kernel, and Figure 3.4.1's bottom panel shows its behavior at an edge: a positive lobe, a negative lobe, and a zero crossing at the transition itself. Zero crossings localize edges more precisely than thresholded gradient peaks, and they form closed contours, properties that made them the foundation of an entire edge-detection school. Their weakness is the theme of this section squared: two rounds of differentiation amplify noise twice as enthusiastically as one, so the raw Laplacian of a photograph is unusable without prior smoothing. Which leads directly to the section's final construction.

4. LoG and DoG: Derivatives With a Scale Dial Advanced

Smooth with a Gaussian, then take the Laplacian. By the associativity from Section 3.1, the two convolutions collapse into one kernel, the Laplacian of Gaussian:

$$ \mathrm{LoG}_\sigma(x, y) = \nabla^2 (G_\sigma * I) = (\nabla^2 G_\sigma) * I, \qquad \nabla^2 G_\sigma(r) = \frac{r^2 - 2\sigma^2}{\sigma^4}\, G_\sigma(r) $$

The kernel looks like a sombrero: a negative center pit ringed by a positive annulus (or the inverse, by sign convention). It responds maximally to blobs, roughly circular regions of size matched to $\sigma$, dark-on-light or light-on-dark, and it crosses zero on their boundaries. The parameter $\sigma$ is now a genuine instrument dial: sweep it and the filter detects structure at chosen scales, the idea formalized by Marr and Hildreth in 1980 as a theory of biological edge detection and later industrialized by scale-space theory in Chapter 4.

Because building exact LoG kernels for many sigmas is wasteful, practice substitutes the Difference of Gaussians: $\mathrm{DoG} = G_{k\sigma} * I - G_{\sigma} * I$, which approximates the LoG closely when $k \approx 1.6$ and costs only Gaussian blurs, which Section 3.6 makes nearly free. A stack of DoG responses across scales is precisely the detection engine inside SIFT, the keypoint detector of Chapter 10; the code below builds one level of it.

import cv2
import numpy as np

gray = cv2.imread("cells.png", cv2.IMREAD_GRAYSCALE).astype(np.float32)

# --- LoG: smooth, then Laplacian (associativity makes this exact) ---
sigma = 2.0
smooth = cv2.GaussianBlur(gray, (0, 0), sigma)
log = cv2.Laplacian(smooth, cv2.CV_32F, ksize=3)

# --- DoG: difference of two Gaussians approximates LoG at ~zero cost ---
k = 1.6                                            # canonical scale ratio
dog = cv2.GaussianBlur(gray, (0, 0), k * sigma) \
    - cv2.GaussianBlur(gray, (0, 0), sigma)

# The two responses are proportional; compare after normalization.
corr = np.corrcoef(log.ravel(), dog.ravel())[0, 1]
print(f"LoG vs DoG correlation: {corr:.3f}")
# Representative output:
# LoG vs DoG correlation: 0.987   (DoG is a positive multiple of LoG here)