Lesson 16: Wavelets and Gabor Filters¶

The Fourier transform (Lesson 15) tells you which frequencies are present in a signal, but not where — a sine wave basis function extends across the entire image, so a localized feature (an edge, a texture patch) gets smeared across the whole spectrum. Wavelets and Gabor filters are two different fixes for this: both are built from small, spatially localized oscillations instead of infinite sinusoids, giving a joint sense of where and what frequency.

In [1]:

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

The Haar wavelet: the simplest possible wavelet¶

The Haar wavelet transform splits a signal into an approximation (local averages) and a detail (local differences), computed on non-overlapping pairs of samples:

$$a_k = \frac{x_{2k} + x_{2k+1}}{\sqrt{2}}, \qquad d_k = \frac{x_{2k} - x_{2k+1}}{\sqrt{2}}$$

This should look familiar: it's essentially the same "blur + keep the residual" idea as the Laplacian pyramid (Lesson 13), just computed pairwise instead of with a Gaussian kernel, and without overlap between neighborhoods.

In [2]:

signal = np.array([1, 3, 5, 11, 7, 9, 2, 4], dtype=np.float64)

approx = (signal[0::2] + signal[1::2]) / np.sqrt(2)
detail = (signal[0::2] - signal[1::2]) / np.sqrt(2)

cA_ref, cD_ref = pywt.dwt(signal, 'haar')
print('manual matches pywt.dwt?', np.allclose(approx, cA_ref) and np.allclose(detail, cD_ref))
print('approximation (cA):', np.round(approx, 2))
print('detail (cD):       ', np.round(detail, 2))

manual matches pywt.dwt? True
approximation (cA): [ 2.83 11.31 11.31  4.24]
detail (cD):        [-1.41 -4.24 -1.41 -1.41]

From Haar to Daubechies: smoother wavelets¶

The Haar wavelet is discontinuous (a hard step), which gives it poor frequency localization — its own frequency content is spread out, the opposite of what we want. Daubechies wavelets (Daubechies, 1988) use longer, smoother filters with more vanishing moments, trading a wider spatial support for much better frequency behavior. db2 (sometimes called "D4" for its 4 filter taps) is the next step up in smoothness from Haar.

In [3]:

db2 = pywt.Wavelet('db2')
print('db2 low-pass (scaling) filter :', np.round(db2.dec_lo, 4))
print('db2 high-pass (wavelet) filter:', np.round(db2.dec_hi, 4))

fig, axes = plt.subplots(1, 2, figsize=(8, 3))
axes[0].stem(db2.dec_lo)
axes[0].set_title('db2 scaling filter (low-pass)')
axes[1].stem(db2.dec_hi)
axes[1].set_title('db2 wavelet filter (high-pass)')
plt.tight_layout()
plt.show()

db2 low-pass (scaling) filter : [-0.1294  0.2241  0.8365  0.483 ]
db2 high-pass (wavelet) filter: [-0.483   0.8365 -0.2241 -0.1294]

No description has been provided for this image

The Haar transform above is just a convolve-and-downsample-by-2 operation with the 2-tap filters $[\tfrac{1}{\sqrt2}, \tfrac{1}{\sqrt2}]$ and $[\tfrac{1}{\sqrt2}, -\tfrac{1}{\sqrt2}]$. Daubechies wavelets follow the exact same recipe with longer filters. We reproduce db2's decomposition from scratch, as a periodic convolution, and check it against pywt.

In [4]:

def dwt_periodic(x, filt):
    """One level of a periodic discrete wavelet transform with an arbitrary filter."""
    n, L = len(x), len(filt)
    k = np.arange(n // 2)[:, None]
    i = np.arange(L)[None, :]
    idx = (2 * k + 2 - i) % n
    return (filt[None, :] * x[idx]).sum(axis=1)

longer_signal = np.array([1, 3, 5, 11, 7, 9, 2, 4, 6, 8], dtype=np.float64)

cA_mine = dwt_periodic(longer_signal, np.array(db2.dec_lo))
cD_mine = dwt_periodic(longer_signal, np.array(db2.dec_hi))
cA_ref, cD_ref = pywt.dwt(longer_signal, db2, mode='periodization')

print('matches pywt (periodization mode)?', np.allclose(cA_mine, cA_ref) and np.allclose(cD_mine, cD_ref))

reconstructed = pywt.idwt(cA_mine, cD_mine, db2, mode='periodization')
print('perfect reconstruction?', np.allclose(reconstructed, longer_signal))

matches pywt (periodization mode)? True
perfect reconstruction? True

2D wavelet decomposition of an image¶

Like the separable Gaussian filter in Lesson 10, a 2D wavelet transform is applied as two 1D passes: rows, then columns. One level of decomposition splits an image into four subbands:

LL: low-pass both directions — a coarser, half-resolution copy of the image (like one Gaussian pyramid level)
LH: low-pass rows, high-pass columns — responds to horizontal edges
HL: high-pass rows, low-pass columns — responds to vertical edges
HH: high-pass both directions — responds to diagonal detail and corners

In [5]:

img = np.zeros((200, 200), dtype=np.float64)
cv2.rectangle(img, (40, 40), (160, 160), 200, -1)
cv2.circle(img, (100, 100), 40, 100, -1)

LL, (LH, HL, HH) = pywt.dwt2(img, 'db2', mode='periodization')

fig, axes = plt.subplots(1, 5, figsize=(14, 3))
for ax, im, title in zip(axes, [img, LL, LH, HL, HH],
                          ['Original', 'LL (approx.)', 'LH (horiz. edges)', 'HL (vert. edges)', 'HH (diagonal)']):
    ax.imshow(im, cmap='gray')
    ax.set_title(title, fontsize=9)
    ax.axis('off')
plt.tight_layout()
plt.show()

Just like the Laplacian pyramid, this decomposition is exactly invertible: pywt.idwt2 reconstructs the original from the four subbands with no loss.

In [6]:

reconstructed_img = pywt.idwt2((LL, (LH, HL, HH)), 'db2', mode='periodization')
print('max reconstruction error:', np.abs(reconstructed_img - img).max())

max reconstruction error: 1.4210854715202004e-13

Gabor filters: localized oscillations tuned to orientation¶

A Gabor filter is a sinusoidal grating multiplied by a Gaussian envelope — a wave that's localized in space, tuned to a specific frequency and orientation. Unlike Daubechies wavelets (built for compact, orthogonal, invertible multi-resolution decomposition), Gabor filters are used more for feature extraction: detecting oriented texture and edges at a chosen scale.

Gabor filters also have a striking biological connection. Recording the neurons in a cat's visual cortex, Hubel and Wiesel (Hubel and Wiesel, 1959) found "simple cells" that fire selectively for a bar or edge at one specific orientation and position, and barely at all for other orientations — a foundational discovery in visual neuroscience. It was later shown (Marcelja, 1980; Daugman, 1985) that a 2D Gabor function is a remarkably good mathematical model of these simple-cell receptive fields, which is part of why Gabor filters became a standard, biologically-motivated tool for orientation-selective feature extraction in computer vision.

In [7]:

orientations_deg = [0, 45, 90, 135]
wavelengths = [4, 8, 12, 16]

fig, axes = plt.subplots(len(wavelengths), len(orientations_deg), figsize=(11, 11))
for i, lambd in enumerate(wavelengths):
    for j, t in enumerate(orientations_deg):
        kernel = cv2.getGaborKernel((25, 25), sigma=3, theta=np.radians(t), lambd=lambd, gamma=0.5, psi=0)
        ax = axes[i, j]
        ax.imshow(kernel, cmap='gray')
        ax.set_xticks([])
        ax.set_yticks([])
        if i == 0:
            ax.set_title(fr'$\theta$ = {t} deg', fontsize=9)
        if j == 0:
            ax.set_ylabel(fr'$\lambda$ = {lambd}', fontsize=9)
plt.tight_layout()
plt.show()

Orientation selectivity, demonstrated quantitatively¶

We draw four bars at four different orientations and filter the image with a Gabor kernel tuned to each orientation, then measure the average response magnitude near each bar. A filter should respond most strongly to the bar matching its tuning and weakly to the others.

In [8]:

def draw_bar(image, center, angle_deg, length=60, thickness=6, value=255):
    angle = np.radians(angle_deg)
    dx, dy = length / 2 * np.cos(angle), length / 2 * np.sin(angle)
    p1 = (int(center[0] - dx), int(center[1] - dy))
    p2 = (int(center[0] + dx), int(center[1] + dy))
    cv2.line(image, p1, p2, value, thickness)

bar_orientations = [0, 45, 90, 135]
centers = [(50, 50), (150, 50), (50, 150), (150, 150)]

bars_img = np.zeros((200, 200), dtype=np.float64)
for c, ang in zip(centers, bar_orientations):
    draw_bar(bars_img, c, ang)

plt.imshow(bars_img, cmap='gray')
plt.title('Bars at 0, 45, 90, 135 degrees')
plt.axis('off')
plt.show()

In [9]:

responses = np.zeros((4, 4))
for fi, bar_angle in enumerate(bar_orientations):
    # Note: cv2's theta parameter is the orientation of the sinusoidal stripes inside the
    # kernel, which runs *perpendicular* to the bar it responds to -- hence the +90 offset.
    kernel = cv2.getGaborKernel((25, 25), sigma=4, theta=np.radians(bar_angle + 90), lambd=10, gamma=0.5, psi=0)
    response = cv2.filter2D(bars_img, cv2.CV_64F, kernel)
    for ci, c in enumerate(centers):
        region = response[c[1] - 20:c[1] + 20, c[0] - 20:c[0] + 20]
        responses[fi, ci] = np.abs(region).mean()

print(f'{"filter tuned for":>18}', ' '.join(f'{a:>6}' for a in bar_orientations), '  <- bar orientation')
for fi, ang in enumerate(bar_orientations):
    print(f'{ang:>18}', ' '.join(f'{v:>6.0f}' for v in responses[fi]))

best_match = [bar_orientations[i] for i in responses.argmax(axis=1)]
print('\neach filter peaks at its own bar orientation?', best_match == bar_orientations)

  filter tuned for      0     45     90    135   <- bar orientation
                 0   2580    612    343    612
                45    495   3277    495    428
                90    343    612   2580    612
               135    495    428    495   3277

each filter peaks at its own bar orientation? True

The response matrix is strongly diagonal-dominant: each filter's largest response lands squarely on the bar it was tuned for, exactly the orientation selectivity Hubel and Wiesel observed biologically.

Exercise¶

Try db4 (pywt.Wavelet('db4')) instead of db2 for the 2D image decomposition. How does the LH/HL/HL subband appearance change, given db4's longer, smoother filters?
Apply pywt.dwt2 a second time to the LL subband from the image decomposition above, to get a second, coarser level — this is a wavelet pyramid, directly analogous to the Gaussian/Laplacian pyramids of Lessons 10 and 12.
Increase lambd (the sinusoid's wavelength) in the Gabor kernel while keeping sigma fixed. What happens to the number of visible stripes inside the Gaussian envelope, and how would you expect that to change which real-image texture scale the filter responds to?

← Lesson 15: The Fourier Transform and Frequency-Domain Filtering

Lesson 17: Image Compression — Lossless and Lossy →