Lesson 13: Second Derivatives and the Laplacian¶
Lesson 12 found edges as peaks in the first derivative (gradient magnitude). This lesson uses the second derivative instead, where edges show up as zero crossings — the basis of the classic Marr-Hildreth edge detector. We then reuse the same blur-and-difference idea to build a Laplacian pyramid, which fixes the information-loss problem from Lesson 11's Gaussian pyramid.
import numpy as np
import cv2
import matplotlib.pyplot as plt
Why zero crossings? A 1D intuition¶
Consider a single step edge along one row. The first derivative is a spike at the edge; the second derivative swings from positive to negative (or vice versa) and crosses exactly zero right at the edge location.
x = np.linspace(0, 10, 400)
step = 1 / (1 + np.exp(-6 * (x - 5))) # a smoothed step, to keep derivatives well-defined
first_deriv = np.gradient(step, x)
second_deriv = np.gradient(first_deriv, x)
fig, axes = plt.subplots(3, 1, figsize=(6, 6), sharex=True)
for ax, y, title in zip(axes, [step, first_deriv, second_deriv],
['Intensity (step edge)', 'First derivative (peak at edge)', 'Second derivative (zero crossing at edge)']):
ax.plot(x, y)
ax.axvline(5, color='gray', linestyle='--', linewidth=1)
ax.set_title(title, fontsize=10)
if True:
axes[2].axhline(0, color='red', linewidth=0.8)
plt.tight_layout()
plt.show()
The Laplacian: a 2D second derivative¶
The Laplacian sums the unmixed second partial derivatives of a 2D image:
$$\nabla^2 I = \frac{\partial^2 I}{\partial x^2} + \frac{\partial^2 I}{\partial y^2}$$
Unlike the gradient, it's a single scalar at each pixel (no direction) — it just measures how much a pixel differs from the average of its neighbors. Two common discrete kernels:
$$K_4 = \begin{bmatrix}0&1&0\\1&-4&1\\0&1&0\end{bmatrix} \qquad K_8 = \begin{bmatrix}1&1&1\\1&-8&1\\1&1&1\end{bmatrix}$$
$K_4$ uses only the 4-connected neighbors; $K_8$ also includes the diagonals for a stronger, more isotropic response.
img = np.zeros((150, 150), dtype=np.float64)
cv2.rectangle(img, (30, 30), (100, 100), 200, -1)
cv2.circle(img, (110, 110), 25, 120, -1)
k4 = np.array([[0, 1, 0], [1, -4, 1], [0, 1, 0]], dtype=np.float64)
lap = cv2.filter2D(img, cv2.CV_64F, k4)
print('matches cv2.Laplacian?', np.allclose(lap, cv2.Laplacian(img, cv2.CV_64F, ksize=1)))
fig, axes = plt.subplots(1, 2, figsize=(7, 3.5))
axes[0].imshow(img, cmap='gray')
axes[0].set_title('Original')
axes[1].imshow(lap, cmap='coolwarm', vmin=-200, vmax=200)
axes[1].set_title('Laplacian response\n(red=positive, blue=negative)')
for ax in axes:
ax.axis('off')
plt.tight_layout()
plt.show()
Every edge produces a positive lobe on one side and a negative lobe on the other — the edge itself sits at the zero crossing between them.
The Marr-Hildreth detector: Laplacian of Gaussian (LoG)¶
The Laplacian, like any derivative, amplifies high-frequency noise — and a second derivative amplifies it even more than a first derivative does. Marr and Hildreth's fix (1980): smooth first, then take the Laplacian. Since convolution is associative, this is equivalent to convolving directly with a single combined kernel, the Laplacian of Gaussian (LoG):
$$\text{LoG}_\sigma(I) = \nabla^2(G_\sigma * I) = (\nabla^2 G_\sigma) * I$$
Edges are then found as zero crossings of the LoG response, rather than by thresholding its magnitude directly.
rng = np.random.default_rng(1)
disk = np.zeros((150, 150), dtype=np.float64)
cv2.circle(disk, (75, 75), 50, 255, -1)
noisy_disk = np.clip(disk + rng.normal(0, 10, disk.shape), 0, 255)
sigma = 2.0
log_response = cv2.Laplacian(cv2.GaussianBlur(noisy_disk, (0, 0), sigmaX=sigma), cv2.CV_64F, ksize=3)
def zero_crossings(response, threshold=4.0):
"""Mark a pixel as an edge if it's adjacent to a sign change of at least `threshold` magnitude."""
edges = np.zeros_like(response, dtype=np.uint8)
shifts = [(0, 1), (1, 0), (1, 1), (1, -1)]
for dy, dx in shifts:
shifted = np.roll(np.roll(response, dy, axis=0), dx, axis=1)
sign_change = (response * shifted < 0) & (np.abs(response - shifted) > threshold)
edges |= sign_change.astype(np.uint8)
return edges * 255
log_edges = zero_crossings(log_response)
fig, axes = plt.subplots(1, 3, figsize=(10, 3.5))
for ax, im, title in zip(axes, [noisy_disk, log_response, log_edges],
['Noisy input', 'LoG response', 'Zero crossings\n(Marr-Hildreth edges)']):
ax.imshow(im, cmap='gray' if title != 'LoG response' else 'coolwarm')
ax.set_title(title, fontsize=9)
ax.axis('off')
plt.tight_layout()
plt.show()
The cheaper approximation: Difference of Gaussians (DoG)¶
Computing a true LoG kernel is more expensive than blurring. A widely used shortcut: the difference of two Gaussian blurs at slightly different scales approximates a scaled LoG:
$$\text{DoG} = G_{\sigma} - G_{k\sigma} \;\approx\; -(k-1)\sigma^2 \cdot \text{LoG}_\sigma$$
This is the same building block later reused by SIFT for keypoint detection at multiple scales.
k = 1.6
g_sigma = cv2.GaussianBlur(noisy_disk, (0, 0), sigmaX=sigma)
g_ksigma = cv2.GaussianBlur(noisy_disk, (0, 0), sigmaX=sigma * k)
dog = (g_ksigma - g_sigma) / ((k - 1) * sigma**2)
correlation = np.corrcoef(log_response.ravel(), dog.ravel())[0, 1]
print(f'correlation between LoG and scaled DoG: {correlation:.3f}')
fig, axes = plt.subplots(1, 2, figsize=(7, 3.5))
axes[0].imshow(log_response, cmap='coolwarm', vmin=-30, vmax=30)
axes[0].set_title('LoG')
axes[1].imshow(dog, cmap='coolwarm', vmin=-30, vmax=30)
axes[1].set_title('DoG (scaled)')
for ax in axes:
ax.axis('off')
plt.tight_layout()
plt.show()
The two are highly correlated but not identical — DoG is an approximation, not an exact substitute, but a much cheaper one (two blurs and a subtraction vs. an explicit second-derivative kernel).
Laplacian pyramids: recovering what Gaussian pyramids throw away¶
Lesson 11 showed that a Gaussian pyramid loses information: blurring and downsampling repeatedly, then trying to upsample back, doesn't reconstruct the original. A Laplacian pyramid fixes this by storing, at each level, exactly the detail that would otherwise be lost — the difference between a Gaussian level and the coarser level upsampled back to match it. That difference is a discrete approximation of the Laplacian, which is why it shares the name.
def gaussian_pyramid(image, num_levels):
pyramid = [image]
current = image
for _ in range(num_levels - 1):
current = cv2.pyrDown(current)
pyramid.append(current)
return pyramid
def laplacian_pyramid(gauss_pyramid):
lap_pyramid = []
for i in range(len(gauss_pyramid) - 1):
size = (gauss_pyramid[i].shape[1], gauss_pyramid[i].shape[0])
upsampled = cv2.pyrUp(gauss_pyramid[i + 1], dstsize=size)
detail = gauss_pyramid[i].astype(np.int16) - upsampled.astype(np.int16)
lap_pyramid.append(detail)
lap_pyramid.append(gauss_pyramid[-1].astype(np.int16)) # the smallest level: no finer detail to subtract
return lap_pyramid
def reconstruct_from_laplacian(lap_pyramid):
current = lap_pyramid[-1]
for i in range(len(lap_pyramid) - 2, -1, -1):
size = (lap_pyramid[i].shape[1], lap_pyramid[i].shape[0])
upsampled = cv2.pyrUp(current, dstsize=size)
current = upsampled + lap_pyramid[i]
return np.clip(current, 0, 255).astype(np.uint8)
photo = np.zeros((256, 256, 3), dtype=np.uint8)
photo[:] = (40, 40, 40)
cv2.circle(photo, (128, 128), 90, (255, 120, 30), -1)
cv2.rectangle(photo, (30, 30), (110, 110), (30, 200, 255), -1)
g_pyr = gaussian_pyramid(photo, num_levels=5)
l_pyr = laplacian_pyramid(g_pyr)
fig, axes = plt.subplots(1, len(l_pyr), figsize=(13, 3))
for ax, level in zip(axes, l_pyr):
display = np.clip(level.astype(np.int32) + 128, 0, 255).astype(np.uint8) # shift for visibility
ax.imshow(display)
ax.set_title(f'{level.shape[1]}x{level.shape[0]}', fontsize=9)
ax.axis('off')
plt.tight_layout()
plt.show()
Each level (shown shifted by +128 gray levels so negative values are visible) is mostly flat gray except right at edges — exactly where detail is lost by blurring and downsampling. The final, smallest level stores actual image content rather than a difference, since there's nothing coarser left to compare it to.
Exact reconstruction¶
Because each level stores exactly what its Gaussian-pyramid counterpart discarded, summing back up the pyramid reconstructs the original image exactly (up to integer rounding).
reconstructed = reconstruct_from_laplacian(l_pyr)
diff = np.abs(reconstructed.astype(int) - photo.astype(int))
print(f'max reconstruction error = {diff.max()} gray levels')
print(f'mean reconstruction error = {diff.mean():.4f} gray levels')
fig, axes = plt.subplots(1, 2, figsize=(6, 3.5))
axes[0].imshow(photo)
axes[0].set_title('Original')
axes[1].imshow(reconstructed)
axes[1].set_title('Reconstructed from Laplacian pyramid')
for ax in axes:
ax.axis('off')
plt.tight_layout()
plt.show()
Compare this to Lesson 11's pyrDown + pyrUp result, which had a visibly nonzero reconstruction error: the Laplacian pyramid stores just enough extra information at each level to make the process perfectly reversible, at the cost of needing to keep all the levels around (not just the smallest one).
Exercise¶
- Try
threshold=1.0andthreshold=15.0inzero_crossings. How does the number of detected edge pixels change, and why does a higher threshold on a second-derivative sign-change criterion behave differently than a magnitude threshold on a first derivative? - Increase the noise level in
noisy_diskand see how largesigmaneeds to be before the zero-crossing edge map stops being dominated by spurious noise loops. - Laplacian pyramids are the classic tool behind seamless image blending (e.g. blending two photos along a mask, level by level). Sketch out — in words or code — how you'd blend two Laplacian pyramids (e.g. average them with a spatially-varying weight per level) before reconstructing, and why blending in this representation avoids sharp seams that blending the original images directly would produce.