Gabor filters: runnable experiments

Companion notebook to the Gabor filters topic. We look at the kernels, use a Gabor bank for orientation and texture analysis, and detect blobs with a Difference of Gaussians. Everything imports from the clean module gabor.py.

import numpy as np
import matplotlib.pyplot as plt
from skimage import data
from skimage.filters import gabor_kernel

from gabor import gabor_response, gabor_feature_stack, dominant_orientation, dog

1. The kernels: even and odd, at several orientations

The real (cosine) part is a bar detector; the imaginary (sine) part is an edge detector. Rotating theta steers both.

fig, axes = plt.subplots(2, 4, figsize=(11, 5.5))
for col, deg in enumerate([0, 45, 90, 135]):
    k = gabor_kernel(frequency=0.2, theta=np.deg2rad(deg), bandwidth=1.0)
    axes[0, col].imshow(k.real, cmap='gray'); axes[0, col].set_title(f'real, {deg}deg')
    axes[1, col].imshow(k.imag, cmap='gray'); axes[1, col].set_title(f'imag, {deg}deg')
    for r in (0, 1):
        axes[r, col].set_xticks([]); axes[r, col].set_yticks([])
fig.suptitle('Gabor kernels (even = real, odd = imaginary)')
fig.tight_layout(); plt.show()

2. Oriented energy on a natural image

Each Gabor channel lights up where the image has structure at that orientation and scale. This oriented-energy stack is the feature a texture classifier starts from.

img = data.camera().astype(float) / 255.0
freqs = [0.1, 0.2, 0.3]
thetas = np.deg2rad([0, 45, 90, 135])
stack, params = gabor_feature_stack(img, freqs, thetas)
print('feature stack:', stack.shape, '=', len(freqs), 'scales x', len(thetas), 'orientations')

fig, axes = plt.subplots(1, 4, figsize=(11, 3.2))
for ax, deg in zip(axes, [0, 45, 90, 135]):
    ax.imshow(gabor_response(img, 0.2, theta=np.deg2rad(deg)), cmap='magma')
    ax.set_title(f'{deg}deg'); ax.set_xticks([]); ax.set_yticks([])
fig.suptitle('Oriented energy at frequency 0.2 cyc/pixel')
fig.tight_layout(); plt.show()

3. Recovering the orientation of a texture

Summing Gabor energy over the image, the peak channel reports the dominant orientation. We test it on a grating with a known angle.

def grating(n=128, frequency=0.1, theta=0.0):
    y, x = np.mgrid[0:n, 0:n]
    xr = x * np.cos(theta) + y * np.sin(theta)
    return np.cos(2 * np.pi * frequency * xr)

true_deg = 60
g = grating(theta=np.deg2rad(true_deg), frequency=0.12)
orients = np.deg2rad(np.arange(0, 180, 10))
energy = [gabor_response(g, 0.12, theta=t).sum() for t in orients]
best = np.rad2deg(dominant_orientation(g, 0.12, orients))
print(f'true {true_deg} deg, recovered {best:.0f} deg')

plt.figure(figsize=(9, 3.2))
plt.plot(np.rad2deg(orients), energy / np.max(energy), 'o-')
plt.axvline(true_deg, color='C3', ls='--', label='true')
plt.xlabel('channel orientation [deg]'); plt.ylabel('normalised energy')
plt.title('Dominant-orientation readout'); plt.legend(); plt.grid(True, alpha=0.3)
plt.tight_layout(); plt.show()

4. Blob detection with Difference of Gaussians

The center-surround DoG is isotropic: it finds spots and blobs of a chosen size, regardless of orientation, and discards smooth gradients.

coins = data.coins().astype(float) / 255.0
d = dog(coins, sigma1=2.0, sigma2=5.0)

fig, (a1, a2) = plt.subplots(1, 2, figsize=(10, 4.5))
a1.imshow(coins, cmap='gray'); a1.set_title('input'); a1.axis('off')
m = np.abs(d).max()
a2.imshow(d, cmap='gray', vmin=-m, vmax=m); a2.set_title('DoG (sigma 2 vs 5)'); a2.axis('off')
fig.tight_layout(); plt.show()

The bright rings trace coin edges and the centres respond as blobs. See the embedded page for running this same 2-D convolution on hardware from an 8-bit AVR up to a Cortex-M33 with an NPU.