Gabor Filter

A Gabor filter is a linear filter whose kernel is a sinusoidal plane wave under a Gaussian envelope. It responds strongly where the image contains a grating of one spatial frequency and orientation near one position, and weakly elsewhere, so it measures “what” (frequency and orientation) and “where” at once. Among all functions it is the one that localizes best in both respects together: it attains the lower bound of the uncertainty relation between spread in space and spread in spatial frequency. Banks of Gabor filters at several scales and orientations are a classical front end for texture description, iris codes, the GIST scene descriptor and face analysis, and Goodfellow, Bengio and Courville (2016) note that most deep learning algorithms learn Gabor-like features in the first layer of a convolutional neural network. The 2D functions are also called Gabor wavelets when they are used as a self-similar family, and Gabor functions or Gabor elementary functions in signal analysis.

Definition

A Gabor filter is a convolution kernel equal to a Gaussian function multiplied by a complex exponential (the complex Gabor filter) or by a cosine or sine (its even and odd real parts). Its parameters are the carrier’s spatial frequency and orientation, the widths of the Gaussian envelope along and across the carrier, the carrier’s phase, and the position of the envelope’s centre. The modulus of the complex filter’s response is the local amplitude or energy of image structure in the filter’s frequency band, and its argument is the local phase.

Intuition

A Gaussian filter is a low-pass filter: its frequency response is a Gaussian centred on zero frequency. Multiplying the kernel by a carrier $e^{i 2\pi f_0 x}$ shifts that Gaussian, by the modulation theorem of the Fourier transform, so that it is centred on frequency $f_0$ instead. The Gabor filter is therefore a band-pass filter whose passband is a Gaussian blob around one point of the 2D frequency plane: the distance of that point from the origin is the preferred spatial frequency, and its direction is the preferred orientation. In space, the kernel looks like a few parallel stripes that fade out with distance from the centre. Convolving it with an image correlates every neighbourhood with a short patch of grating.

The width of the envelope sets the trade-off. A wide envelope contains many stripes: it is selective in frequency and orientation but blurs position. A narrow one localizes well but responds to a broad range of frequencies and orientations. No choice escapes the trade-off, and the Gaussian envelope makes it as favourable as possible.

The real, even filter (a cosine under the envelope) responds to bright or dark bars centred on it; the odd filter (a sine) responds to edges. Both responses oscillate as a grating slides under the filter. Their squared sum does not: the even and odd filters form a quadrature pair, and the modulus of the complex filter’s response measures how much grating there is, whatever its phase. Goodfellow, Bengio and Courville (2016, §9.10) describe this as the cartoon model of a “complex cell” in the primary visual cortex, computed from two “simple cells” a quarter cycle out of phase.

Formal Definition

The kernel

Let $(x, y)$ be pixel coordinates and rotate them by the filter’s orientation $\theta$:

$$x' = x\cos\theta + y\sin\theta, \qquad y' = -x\sin\theta + y\cos\theta .$$

The complex Gabor kernel, in the convention scikit-image uses, is

$$g(x, y) = \frac{1}{2\pi\sigma_x\sigma_y} \exp\!\left(-\frac{x'^2}{2\sigma_x^2} - \frac{y'^2}{2\sigma_y^2}\right) \exp\!\big(i\,(2\pi f_0 x' + \psi)\big),$$

with carrier frequency $f_0$ in cycles per pixel (wavelength $\lambda = 1/f_0$), envelope standard deviations $\sigma_x$ along the carrier and $\sigma_y$ across it, and phase offset $\psi$. The orientation $\theta$ is the direction of the carrier’s wave vector, which is normal to the stripes. With $\psi = 0$ the real part is even (cosine) and the imaginary part odd (sine).

The real form used in the computational neuroscience literature and by OpenCV writes the envelope with one width $\sigma$ and a spatial aspect ratio $\gamma$,

$$g(x, y) = \exp\!\left(-\frac{x'^2 + \gamma^2 y'^2}{2\sigma^2}\right)\cos\!\left(2\pi\frac{x'}{\lambda} + \varphi\right),$$

so $\sigma_y = \sigma/\gamma$. The Groningen group’s page that scikit-image’s source cites gives this form, notes that $\varphi = 0$ and $180°$ give symmetric functions and $\pm 90°$ antisymmetric ones, and reports that $\gamma$ has been found to vary between about 0.23 and 0.92 for simple cells (Grigorescu and Petkov). Goodfellow, Bengio and Courville (2016, eq. 9.16) write the same function with a general amplitude, centre and rotation as the model of a V1 simple cell’s weights.

Daugman (2003, eq. 2.4) defines the complex 2D Gabor wavelet with widths $(\alpha, \beta)$ in an $e^{-\pi(\cdot)}$ normalization and a wave vector $(u_0, v_0)$, whose polar coordinates are the preferred spatial frequency and orientation. His lecture notes add that a further degree of freedom, the orientation of the elliptical envelope relative to the carrier, creates cross terms in $xy$ (Daugman, 2021, slide 80); the forms above align the envelope with the carrier.

Fourier transform

With the transform $G(u, v) = \iint g(x,y)\,e^{-i2\pi(ux + vy)}\,dx\,dy$ and $\psi = 0$, the normalized kernel above has

$$G(u, v) = \exp\!\Big(-2\pi^2\big[\sigma_x^2 (u' - f_0)^2 + \sigma_y^2 v'^2\big]\Big), \qquad u' = u\cos\theta + v\sin\theta,\; v' = -u\sin\theta + v\cos\theta .$$

This follows from the transform of a Gaussian and the modulation theorem: a Gaussian of the same functional form, centred on the wave vector $f_0(\cos\theta, \sin\theta)$, with peak value 1 there. Daugman (2003, eq. 2.5) states the same property for his parameterization: the 2D Gabor wavelet is self-Fourier, its transform has the same functional form with the parameters interchanged or inverted, and its power spectrum is a bivariate Gaussian centred on $(u_0, v_0)$. A real (even or odd) Gabor filter has two such lobes, at $\pm f_0(\cos\theta, \sin\theta)$, because a cosine is the sum of two complex exponentials (Daugman, 2021, slide 81).

The response of the complex filter to a real grating $\cos(2\pi \mathbf{k}\cdot\mathbf{x} + \phi)$ with wave vector $\mathbf{k}$ is $\tfrac{1}{2}\big(G(\mathbf{k})e^{i(\cdot)} + G(-\mathbf{k})e^{-i(\cdot)}\big)$. When $G(-\mathbf{k})$ is negligible, which holds whenever the passband stays away from the origin, the modulus is $\tfrac12 G(\mathbf{k})$ at every position: the energy $|r|^2$ is independent of the grating’s phase, and $G$ is directly the filter’s tuning curve.

Bandwidth

Along the carrier, $G$ falls to half its peak at $f_0 \pm \sqrt{\ln 2 / 2}\,/(\pi\sigma_x)$. Expressed as a ratio of the upper to the lower half-magnitude frequency, the bandwidth in octaves is

$$b = \log_2 \frac{\pi\sigma_x f_0 + \sqrt{\ln 2/2}}{\pi\sigma_x f_0 - \sqrt{\ln 2/2}} \quad\Longleftrightarrow\quad \sigma_x f_0 = \frac{1}{\pi}\sqrt{\frac{\ln 2}{2}}\;\frac{2^b + 1}{2^b - 1},$$

the relation the Groningen page gives with $f_0 = 1/\lambda$ and the one scikit-image’s gabor_kernel uses to set $\sigma$ from its bandwidth argument. A one-octave filter has $\sigma \approx 0.56\lambda$, so its envelope spans a little more than two carrier cycles at $\pm 2\sigma$. Narrower envelopes widen the bandwidth without limit as $\sigma_x f_0$ approaches $\sqrt{\ln 2/2}/\pi \approx 0.19$, where the lower half-magnitude point reaches zero frequency.

At the preferred frequency, with an isotropic envelope ($\sigma_x = \sigma_y = \sigma$), the response falls to half at the angle

$$\Delta\theta_{1/2} = \arccos\!\left(1 - \frac{\ln 2}{4\pi^2\sigma^2 f_0^2}\right),$$

which follows from the distance between two points on the circle of radius $f_0$. For one octave this is about $19.2°$, so the full orientation bandwidth is about $38°$. With an elongated envelope ($\sigma_y > \sigma_x$, as with $\gamma < 1$), the orientation tuning narrows while the frequency tuning stays set by $\sigma_x$ (our reading of the formula for $G$).

Properties

Joint localization: the uncertainty bound

For any filter, let $\Delta x$ be the standard deviation of $|g|^2$ along an axis and $\Delta u$ the standard deviation of $|G|^2$ along the corresponding frequency axis. Then $\Delta x\,\Delta u \ge 1/(4\pi)$, with equality only for a Gaussian times a complex exponential; this is the uncertainty principle that the Fourier Transform article states for the Gaussian. Daugman’s lecture notes credit Dennis Gabor with discovering these uncertainty-minimizing elementary functions in one dimension (Daugman, 2021, slide 55; Gabor, 1946), and Daugman (2003) describes the 2D Gabor wavelet as a generalization of them, citing his 1980 paper (Daugman, 1980), and states that it achieves the lower bound of the uncertainty relation for joint resolution in the 2D space and 2D Fourier domains, citing his 1985 paper. In two dimensions the bound applies to each pair of conjugate axes, so with $\Delta x \Delta y$ for the spatial extent and $\Delta u \Delta v$ for the spectral extent the product is at least $1/(16\pi^2)$.

The example below uses a rotation-free form of the 2D bound. With $\Sigma_x$ the $2\times2$ covariance of $|g|^2$ and $\Sigma_u$ that of $|G|^2$,

$$\rho = 16\pi^2\sqrt{\det\Sigma_x\,\det\Sigma_u} \;\ge\; 1 .$$

This form is our derivation: the joint covariance of position and frequency has determinant at least $(4\pi)^{-4}$ (the Robertson–Schrödinger form of the uncertainty relation), and by Fischer’s inequality $\det\Sigma_x \det\Sigma_u$ is at least that determinant. Equality holds for Gabor functions, whose position and frequency are uncorrelated.

The bound concerns the complex filter. A real Gabor filter’s power spectrum has two lobes, so its spectral spread about the mean frequency (zero) includes the distance $f_0$ between the lobes and the origin, and the product is far above the bound.

Response to constant intensity

The complex Gabor kernel, and its even part, do not integrate to zero: the even part’s response to a uniform image is $G$ evaluated at the origin, $\exp(-2\pi^2\sigma_x^2 f_0^2)$ in the normalization above. This is small for narrow bandwidths and grows quickly as the envelope narrows. Kovesi points out that once the standard deviation of the two Gaussians that make up the even filter’s transfer function exceeds about a third of the centre frequency, they overlap excessively at the origin, which limits Gabor filters with a reasonably small DC component to about one octave of bandwidth, so many filters are needed to cover the spectrum. The odd part is antisymmetric and always has zero DC response. Implementations that need a zero-mean even filter subtract the DC term, as the iris-code example in the biometrics article does.

Amplitude and phase

Convolving with the complex filter gives at every pixel a complex number whose modulus is the local contrast amplitude and whose angle is the local phase (Daugman, 2003, eqs. 2.6–2.7). Amplitude is the usual feature for texture and detection: Daugman’s lecture notes show the modulus after convolution with complex 2D Gabor wavelets picking out the eyes and mouth in a face image (slides 86–88) and texture-energy variation across an image as a basis for segmentation (slide 104). Phase carries position: shifting a grating by a fraction of a wavelength changes the phase by the same fraction of a cycle and leaves the amplitude unchanged. Daugman’s iris codes keep only phase, quantized to its quadrant (Daugman, 2003).

Completeness without orthogonality

With dilation, rotation and translation, a 2D Gabor mother wavelet generates a self-similar family that Daugman (2003) describes as closed under these transformations and under convolution, and as a complete image basis, or frame. The members are not orthogonal, so the coefficients that reconstruct an image are not simply the inner products of the image with the wavelets; Daugman (2003) cites his 1988 relaxation network as one solution (Daugman, 1988). His lecture notes show a face reconstructed from a few hundred wavelets with six orientations, two quadrature phases and five frequencies an octave apart (Daugman, 2021, slide 84).

Sampling and truncation

A digital Gabor filter samples the continuous function on the pixel grid and truncates it. Sampling aliases the spectrum, which matters only when the passband approaches the Nyquist frequency of 0.5 cycles per pixel; Kovesi prefers wavelengths of at least 3 pixels for that reason. Truncation multiplies the kernel by a box, which convolves the spectrum with a sinc and changes the tuning. scikit-image’s gabor_kernel truncates by default to a rectangle whose half-width along each image axis is the larger of the projections onto that axis of $3\sigma_x$ along the carrier and $3\sigma_y$ across it (n_stds=3), which is $3\sigma$ at $0°$ but less in oblique directions, and scales the kernel by $1/(2\pi\sigma_x\sigma_y)$; OpenCV’s getGaborKernel uses the same $3\sigma$ rule, rounded rather than rounded up, unless a size is given, returns the real form with no normalization and a default phase of $\pi/2$ (an odd kernel), and stores the kernel rotated by 180°, consistent with filter2D computing correlation rather than convolution. The example measures what the $3\sigma$ default costs.

Filter Banks

A single Gabor filter measures one band. Applications use a bank: several orientations at each of several scales, with one or both quadrature phases. The choices:

• Orientations. Energy is periodic with period $180°$, so $K$ orientations are spaced $180°/K$ apart: six in Daugman’s face example and in Kovesi’s own practice, eight in the GIST descriptor’s default bank (see scene recognition), and, according to Daugman’s lecture notes, usually eight for facial action coding (slide 211).
• Scales. Centre frequencies are spaced geometrically, for example an octave apart in Daugman’s face example (slide 84). For one-octave filters that spacing makes the half-magnitude points of neighbouring filters meet: a filter at $f_0$ falls to half at $\tfrac{4}{3}f_0$, where the filter at $2f_0$ does (from the bandwidth formula). Kovesi frames the choice as a trade-off: enough overlap that the sum of the transfer functions covers the spectrum evenly, but as little as possible so that the outputs of neighbouring filters are not redundant.
• Phase. Even and odd filters together give amplitude and phase; one phase alone confounds them.
• Implementation. Small banks are applied as spatial convolutions; large ones can be built directly as transfer functions and applied by multiplying the image’s spectrum, as the GIST code and Kovesi’s log-Gabor code do. Kovesi notes that a frequency-domain filter with a single lobe on one side of the origin gives, after the inverse transform, the even-filter response in the real part and the odd-filter response in the imaginary part.

Bank outputs are turned into features by pooling the energy over regions (mean and variance per filter are a classical texture feature), by quantizing phase (iris codes), or by sampling the energies on a grid (GIST). Estimating an orientation more finely than the bank’s spacing needs interpolation between filters: since orientation is periodic with period $180°$, the energies are averaged as vectors at twice their angles, $\hat\theta = \tfrac12 \arg \sum_k E_k e^{2i\theta_k}$, rather than taking the filter with the largest energy.

Variants

• Even, odd and complex filters. The even and odd real filters are used alone when only bars or only edges matter; the complex filter, or the pair, gives phase-invariant energy.
• Gabor wavelets. A self-similar family generated from one mother wavelet by dilations in octave steps, rotations and translations, with a $2^{-2m}$ amplitude factor at dilation level $m$ (Daugman, 2003, eqs. 2.1–2.3).
• Log-Gabor filters. Kovesi describes the log-Gabor function, which he attributes to Field (1987), as a filter whose transfer function is Gaussian on a logarithmic frequency axis, $\exp\!\big(-\ln^2(f/f_0) / (2\ln^2\kappa)\big)$ for a constant shape ratio $\kappa$. It has no DC component by construction, has an extended tail toward high frequencies, and can be built with any bandwidth; it has no closed form in space and is designed in the frequency domain. Kovesi gives $\kappa \approx 0.75$ for about one octave and $0.55$ for about two, and reports that the spatial width of log-Gabor filters is near its minimum for bandwidths of about one to three octaves. He also reports that log-Gabor and Gabor functions have almost identical shapes below one octave.
• Learned Gabor layers. Pérez et al. (2020) replace the first convolutional layers of LeNet, AlexNet, VGG16 and WideResNet by Gabor layers: kernels sampled from a real Gabor function whose scale, envelope width, aspect ratio, carrier frequency and phase are learned, at a fixed set of orientations spread uniformly over the circle, applied depth-wise and followed by a $1\times1$ convolution. They report that networks with Gabor layers gain robustness against adversarial examples while keeping test accuracy, and derive a Lipschitz constant for the filters that they use as a regularizer.

Examples

The script measures three things on synthetic gratings, whose frequency, orientation and phase are known exactly. Part A compares a Gabor filter’s measured frequency and orientation tuning with the closed form above, for scikit-image’s default truncation (n_stds=3) and for n_stds=5, and measures the DC response of the even kernel. Part B computes the uncertainty ratio $\rho$ for Gabor filters and for other band-pass filters. Part C estimates the orientation of noisy gratings from the energies of a bank of eight filters.

Conventions: $x$ runs along columns and $y$ down the rows; $\theta$ is the direction of the wave vector measured from the $x$ axis toward $+y$, as in scikit-image; frequencies are in cycles per pixel; filtering is convolution, as in skimage.filters.gabor and scipy.signal.fftconvolve. In Part A the energy at the centre of each grating is averaged over two grating phases, $0$ and $\pi/2$, which cancels the small term that beats with position (the $G(-\mathbf{k})$ term above) exactly, so the expected value is $\big(G(\mathbf{k})^2 + G(-\mathbf{k})^2\big)/4$. As checks, the script compares the hand-written kernel with scikit-image’s gabor_kernel, the centre sum with scipy.ndimage.convolve, and each measured energy with the exact energy computed from the frequency response of the sampled, truncated kernel. The three blocks form one script; it runs in about 3 seconds with NumPy 2.3.5, SciPy 1.16.3 and scikit-image 0.25.2.

import os
os.environ["OMP_NUM_THREADS"] = "1"          # one BLAS thread: identical sums on every run
os.environ["OPENBLAS_NUM_THREADS"] = "1"
os.environ["MKL_NUM_THREADS"] = "1"
import math
import numpy as np
from scipy import ndimage
from scipy.optimize import brentq
from scipy.signal import fftconvolve
from skimage.filters import gabor_kernel


def prefactor(b):
    """sigma * f for a half-magnitude frequency bandwidth of b octaves."""
    return math.sqrt(math.log(2) / 2) / math.pi * (2**b + 1) / (2**b - 1)


def gabor(f, theta, b=1.0, sx=None, sy=None, n_stds=3):
    """Complex Gabor kernel in scikit-image's convention (x = columns, y = rows)."""
    sx = prefactor(b) / f if sx is None else sx
    sy = prefactor(b) / f if sy is None else sy
    ct, st = math.cos(theta), math.sin(theta)
    x0 = math.ceil(max(abs(n_stds * sx * ct), abs(n_stds * sy * st), 1))
    y0 = math.ceil(max(abs(n_stds * sy * ct), abs(n_stds * sx * st), 1))
    y, x = np.mgrid[-y0:y0 + 1, -x0:x0 + 1]
    xr, yr = x * ct + y * st, -x * st + y * ct
    return np.exp(-0.5 * (xr**2 / sx**2 + yr**2 / sy**2) + 2j * np.pi * f * xr) / (2 * np.pi * sx * sy)


def grid(g):
    h, w = g.shape[0] // 2, g.shape[1] // 2
    return np.mgrid[-h:h + 1, -w:w + 1]                 # (y, x) offsets of the kernel taps


# ---------- A. tuning curves against the closed form ----------
F0, TH0 = 1 / 8, math.radians(30)
SIG = prefactor(1.0) / F0                               # 4.497 px for 1 octave at f0 = 1/8

def G_closed(kx, ky):
    """Continuous Fourier transform of the kernel at wave vector (kx, ky), cycles/px."""
    u = kx * math.cos(TH0) + ky * math.sin(TH0) - F0
    v = -kx * math.sin(TH0) + ky * math.cos(TH0)
    return math.exp(-2 * math.pi**2 * SIG**2 * (u**2 + v**2))

def E_closed(f, th):
    kx, ky = f * math.cos(th), f * math.sin(th)
    return (G_closed(kx, ky)**2 + G_closed(-kx, -ky)**2) / 4

def E_measured(g, f, th):
    """Energy |response|^2 at the centre of a grating cos(2 pi f (x cos th + y sin th) + phi),
    averaged over phi = 0 and pi/2 (which cancels the term that beats with position)."""
    y, x = grid(g)
    e = 0.0
    for phi in (0.0, math.pi / 2):
        patch = np.cos(2 * math.pi * f * (x * math.cos(th) + y * math.sin(th)) + phi)
        r = np.sum(g * patch[::-1, ::-1])               # convolution evaluated at the centre
        e += abs(r)**2 / 2
    return e

def E_dtft(g, f, th):
    """Exact energy from the frequency response of the sampled, truncated kernel."""
    y, x = grid(g)
    H = lambda kx, ky: np.sum(g * np.exp(-2j * math.pi * (kx * x + ky * y)))
    kx, ky = f * math.cos(th), f * math.sin(th)
    return (abs(H(kx, ky))**2 + abs(H(-kx, -ky))**2) / 4

g3 = gabor(F0, TH0)
print(f"hand kernel vs skimage gabor_kernel: max abs difference {np.max(np.abs(g3 - gabor_kernel(F0, theta=TH0))):.1e}")
img = np.cos(2 * math.pi * 0.1 * (np.mgrid[0:65, 0:65][1] * math.cos(0.3) + np.mgrid[0:65, 0:65][0] * math.sin(0.3)) + 0.7)
r_nd = ndimage.convolve(img, g3.real)[32, 32] + 1j * ndimage.convolve(img, g3.imag)[32, 32]
y, x = grid(g3)
r_c = np.sum(g3 * img[32 - y.max():33 + y.max(), 32 - x.max():33 + x.max()][::-1, ::-1])
print(f"centre sum vs scipy.ndimage.convolve: difference {abs(r_nd - r_c):.1e}")

freqs = np.arange(0, 0.3001, 0.005)
dths = np.radians(np.arange(-90, 90.01, 2.5))
peak = E_closed(F0, TH0)
print("\nn_stds  kernel  max|meas-closed|/peak (freq, orient)  max|meas-DTFT|  octaves  orient. half-width")
for n in (3, 5):
    g = gabor(F0, TH0, n_stds=n)
    mf = np.array([E_measured(g, f, TH0) for f in freqs])
    mo = np.array([E_measured(g, F0, TH0 + d) for d in dths])
    dev_f = np.max(np.abs(mf - [E_closed(f, TH0) for f in freqs])) / peak
    dev_o = np.max(np.abs(mo - [E_closed(F0, TH0 + d) for d in dths])) / peak
    dev_x = max(np.max(np.abs(mf - [E_dtft(g, f, TH0) for f in freqs])),
                np.max(np.abs(mo - [E_dtft(g, F0, TH0 + d) for d in dths])))
    a = lambda f: math.sqrt(E_measured(g, f, TH0)); a0 = a(F0)
    flo, fhi = brentq(lambda f: a(f) - a0 / 2, 0.02, F0), brentq(lambda f: a(f) - a0 / 2, F0, 0.3)
    dh = brentq(lambda d: math.sqrt(E_measured(g, F0, TH0 + d)) - a0 / 2, 0.01, 1.2)
    print(f"{n:6d}  {g.shape[0]}x{g.shape[1]}   {dev_f:.1e}, {dev_o:.1e}{'':17s}{dev_x:.1e}     {math.log2(fhi / flo):.4f}   {math.degrees(dh):.3f} deg")
c = math.sqrt(math.log(2) / 2) / (math.pi * SIG)
print(f"closed form{'':48s}{math.log2((F0 + c) / (F0 - c)):.4f}   "
      f"{math.degrees(math.acos(1 - math.log(2) / (4 * math.pi**2 * SIG**2 * F0**2))):.3f} deg")
for b in (1, 2, 3):
    s = prefactor(b) / F0
    print(f"DC response of the even (real) kernel, bandwidth {b} octave(s): "
          f"{np.sum(gabor(F0, TH0, b, n_stds=5).real):.6f} (closed form {math.exp(-2 * math.pi**2 * s**2 * F0**2):.6f})")

hand kernel vs skimage gabor_kernel: max abs difference 0.0e+00
centre sum vs scipy.ndimage.convolve: difference 2.8e-17

n_stds  kernel  max|meas-closed|/peak (freq, orient)  max|meas-DTFT|  octaves  orient. half-width
     3  25x25   2.1e-02, 2.1e-02                 1.2e-16     1.0164   19.510 deg
     5  41x41   2.0e-05, 2.0e-05                 1.4e-16     1.0000   19.188 deg
closed form                                                1.0000   19.188 deg
DC response of the even (real) kernel, bandwidth 1 octave(s): 0.001954 (closed form 0.001953)
DC response of the even (real) kernel, bandwidth 2 octave(s): 0.145817 (closed form 0.145816)
DC response of the even (real) kernel, bandwidth 3 octave(s): 0.317965 (closed form 0.317965)

The measured energies equal the sampled kernel’s exact frequency response to about $10^{-16}$, so every difference from the closed form comes from sampling and truncation, not from the measurement. With $5\sigma$ the curves match the continuous formula to $2\times10^{-5}$ of the peak, and the half-magnitude frequency bandwidth and orientation half-width equal their closed forms (1 octave, $19.19°$). scikit-image’s default $3\sigma$ box changes the tuning by 2% of the peak and makes the filter slightly broader than requested: 1.016 octaves instead of 1 and a half-width of $19.51°$ instead of $19.19°$. The 25 × 25 box of the $30°$ filter ends 12 pixels from the centre along the image axes, about $2.7\sigma$, where the envelope is still about 3% of its peak, and cutting it off with a box spreads the passband; this matters when tuning curves themselves are the point, not for most feature extraction. The DC responses follow $\exp(-2\pi^2\sigma^2 f_0^2)$: 0.2% of the peak gain at one octave but 15% at two octaves and 32% at three, the overlap at the origin that limits Gabor filters with small DC response to about one octave.

# ---------- B. joint space-frequency uncertainty ----------
def uncertainty_ratio(g, n_fft=512):
    """16 pi^2 sqrt(det Cov_x det Cov_u) of |g|^2 and |G|^2; >= 1, = 1 only for Gabor functions."""
    y, x = grid(g)
    p = (np.abs(g)**2).ravel(); p /= p.sum()
    X = np.stack([x.ravel(), y.ravel()]).astype(float); X -= X @ p[:, None]
    q = (np.abs(np.fft.fft2(g, (n_fft, n_fft)))**2).ravel(); q /= q.sum()
    fy, fx = np.meshgrid(np.fft.fftfreq(n_fft), np.fft.fftfreq(n_fft), indexing="ij")
    U = np.stack([fx.ravel(), fy.ravel()]); U -= U @ q[:, None]
    return 16 * math.pi**2 * math.sqrt(np.linalg.det((X * p) @ X.T) * np.linalg.det((U * q) @ U.T))

def hann_carrier(f, theta, R=12):
    y, x = np.mgrid[-R:R + 1, -R:R + 1]
    r = np.hypot(x, y)
    w = np.where(r < R, np.cos(math.pi * r / (2 * R))**2, 0.0)
    return w * np.exp(2j * math.pi * f * (x * math.cos(theta) + y * math.sin(theta)))

def log_gabor(f, theta, sigma_on_f=0.75, dtheta_sigma=math.pi / 8 / 1.5, n=256):
    """One-sided log-Gabor built in the frequency domain (Kovesi's parameterization), returned in space."""
    fy, fx = np.meshgrid(np.fft.fftfreq(n), np.fft.fftfreq(n), indexing="ij")
    rad = np.hypot(fx, fy); rad[0, 0] = 1.0
    radial = np.exp(-np.log(rad / f)**2 / (2 * math.log(sigma_on_f)**2)); radial[0, 0] = 0.0
    dth = np.abs(np.angle(np.exp(1j * (np.arctan2(fy, fx) - theta))))
    H = radial * np.exp(-dth**2 / (2 * dtheta_sigma**2))
    g = np.fft.fftshift(np.fft.ifft2(H))
    return g[1:, 1:]                                    # odd size, centred on the origin

rows = [("Gabor, sigma 4.50, f 1/8, 30 deg, untruncated (8 sigma)", gabor(F0, TH0, n_stds=8)),
        ("Gabor, sigma 3 x 6, f 0.15, 60 deg, untruncated (8 sigma)", gabor(0.15, math.radians(60), sx=3.0, sy=6.0, n_stds=8)),
        ("Gabor, sigma 4.50, truncated at 3 sigma (scikit-image default)", gabor(F0, TH0, n_stds=3)),
        ("Hann-windowed complex carrier, radius 12, f 1/8", hann_carrier(F0, TH0)),
        ("log-Gabor, f 1/8, sigma_on_f 0.75 (about 1 octave)", log_gabor(F0, TH0)),
        ("even (real) part of the first Gabor", gabor(F0, TH0, n_stds=8).real.astype(complex))]
print("\nfilter                                                           16 pi^2 sqrt(det Sx det Su)")
for name, g in rows:
    print(f"{name:65s}{uncertainty_ratio(g):.6f}")

filter                                                           16 pi^2 sqrt(det Sx det Su)
Gabor, sigma 4.50, f 1/8, 30 deg, untruncated (8 sigma)          1.000000
Gabor, sigma 3 x 6, f 0.15, 60 deg, untruncated (8 sigma)        1.000000
Gabor, sigma 4.50, truncated at 3 sigma (scikit-image default)   1.006244
Hann-windowed complex carrier, radius 12, f 1/8                  1.068259
log-Gabor, f 1/8, sigma_on_f 0.75 (about 1 octave)               1.141693
even (real) part of the first Gabor                              5.094186

Both untruncated Gabor filters, one isotropic and one elongated and rotated, attain the bound $\rho = 1$ to the printed precision. Everything else exceeds it: the $3\sigma$ truncation slightly (1.006), a Hann window under the same complex carrier by 7%, and the log-Gabor filter by 14%. The even (real) Gabor filter, which is the real part of the first row, gives 5.09, because its power spectrum has two lobes on opposite sides of the origin. The log-Gabor’s larger product is not a defect for its purpose: it was designed for wide bandwidths and zero DC response, not for the minimum joint spread.

# ---------- C. local orientation of noisy gratings from a bank's energy ----------
N_IMG, SZ, LAM, K = 500, 64, 8.0, 8
rng = np.random.default_rng(0)
thetas = rng.uniform(0, np.pi, N_IMG)                  # wave-vector direction, x = columns, y = rows
phis = rng.uniform(0, 2 * np.pi, N_IMG)
yy, xx = np.mgrid[0:SZ, 0:SZ]
clean = np.cos(2 * np.pi / LAM * (xx[None] * np.cos(thetas)[:, None, None]
                                  + yy[None] * np.sin(thetas)[:, None, None]) + phis[:, None, None])
PTS = np.array([20, 28, 36, 44])                         # 4 x 4 sample points, kernels stay inside
TK = np.arange(K) * np.pi / K
bank = [gabor_kernel(1 / LAM, theta=t, bandwidth=1.0) for t in TK]
print(f"\ndata: {N_IMG} gratings {SZ}x{SZ}, wavelength {LAM:.0f} px, theta {np.degrees(thetas.min()):.1f}-"
      f"{np.degrees(thetas.max()):.1f} deg, clean range {clean.min():.3f}..{clean.max():.3f}; "
      f"bank of {K}, largest kernel {max(k.shape[0] for k in bank)} px")

def angular_error(est, true):
    d = np.abs(est - true[:, None]) % np.pi
    return np.degrees(np.minimum(d, np.pi - d))

print("noise sd   argmax mean / p95     vector mean / p95     even-only mean / p95   (degrees)")
for sd in (0.0, 1.0, 2.0, 4.0, 8.0):
    imgs = clean + sd * rng.standard_normal(clean.shape)
    R = np.stack([[fftconvolve(im, k, mode="same")[np.ix_(PTS, PTS)].ravel() for k in bank] for im in imgs])
    E, Ee = np.abs(R)**2, R.real**2                      # (N_IMG, K, 16): quadrature and even-only energy
    est = {"argmax": TK[np.argmax(E, axis=1)],
           "vector": 0.5 * np.angle(np.einsum("nkp,k->np", E, np.exp(2j * TK))) % np.pi,
           "even": 0.5 * np.angle(np.einsum("nkp,k->np", Ee, np.exp(2j * TK))) % np.pi}
    errs = [angular_error(est[m], thetas) for m in ("argmax", "vector", "even")]
    print(f"{sd:7.1f}   " + "   ".join(f"{e.mean():8.2f} / {np.percentile(e, 95):6.2f}" for e in errs))

data: 500 gratings 64x64, wavelength 8 px, theta 0.1-179.5 deg, clean range -1.000..1.000; bank of 8, largest kernel 29 px
noise sd   argmax mean / p95     vector mean / p95     even-only mean / p95   (degrees)
    0.0       5.86 /  10.66       0.31 /   0.62       0.31 /   0.62
    1.0       5.94 /  11.23       1.55 /   3.91       7.38 /  46.13
    2.0       6.26 /  12.69       4.40 /  12.19      14.45 /  67.88
    4.0      17.45 /  74.97      18.90 /  69.46      27.43 /  80.38
    8.0      35.32 /  83.34      35.42 /  82.57      38.41 /  83.83

Without noise, the argmax estimator can only return one of the eight filter orientations, so its errors spread between zero and half the $22.5°$ spacing (95th percentile 10.66°, mean 5.86°, close to the 5.6° of a uniform spread). The doubled-angle vector average interpolates between filters and is accurate to a fraction of a degree; the remaining error comes from the angular tuning curve being sampled at only eight orientations (in our reading). The even-only energy gives the same accuracy without noise, because on a noise-free grating all eight even responses share the same phase factor (apart from the small $G(-\mathbf{k})$ term), which cancels in the angle.

Noise separates the estimators. The unit-amplitude gratings give a response of modulus about 0.5 at the matched filter, while the noise in each filter’s output has standard deviation equal to the pixel noise times the kernel’s norm, about $0.063$ for $\sigma = 4.5$ (this follows from the kernel’s normalization: its squared norm is close to $1/(4\pi\sigma^2)$). At noise standard deviation 1 and 2 the vector average is more accurate than the argmax, as expected. At 4 it is not: its mean error (18.90°) is larger than the argmax’s (17.45°), although its 95th percentile is smaller; at 8 the two are about equal, with mean errors near 35° against 45° for a random guess. In our reading, noise energy enters all eight channels, which shortens the averaged vector and makes its angle noisy, while the argmax needs only the matched channel to stay largest. The even-only estimator is worse than the quadrature estimators at every noise level: already at noise 1 its 95th percentile error is 46°. Its energy is the quadrature energy times the squared cosine of the local phase, so at sample points near a zero crossing of the even response the signal vanishes and noise decides the angle. That is the practical case for quadrature pairs.

History

Daugman’s lecture notes credit Dennis Gabor with discovering the uncertainty-minimizing elementary functions now called wavelets (Daugman, 2021, slide 55), and Daugman (2003) cites Gabor’s 1946 paper “Theory of Communication” for the one-dimensional elementary functions (Gabor, 1946). Daugman (2003) describes the 2D Gabor wavelet as a generalization of those functions, citing his own 1980 paper (Daugman, 1980), and cites his 1980 and 1985 papers for two of its properties: it models the receptive field profiles of simple cells in the mammalian visual cortex, and it achieves the lower bound of the uncertainty relation in 2D (Daugman, 1985). Goodfellow, Bengio and Courville (2016, §9.10) state that reverse-correlation measurements show most V1 cells to have weights described by Gabor functions, and that, since Olshausen and Field (1996) showed that sparse coding learns features with receptive fields similar to those of simple cells, a wide variety of learning algorithms, including most deep learning algorithms in their first layer, have been found to learn Gabor-like features from natural images; they warn that, because so many different algorithms learn such edge detectors, the features alone do not identify the right model of the brain. The image filtering article describes the frequency- and orientation-selective first-layer kernels of the ImageNet classifier of Krizhevsky, Sutskever and Hinton (2012).

In computer vision, Pérez et al. (2020) describe Gabor filters as having become customary, used for texture characterization, character recognition, edge detection and face recognition. Daugman’s iris recognition is built on them, and the GoogLeNet paper article notes that Serre et al. (2007), modelling the primate visual cortex, had used fixed Gabor filters of several sizes.

Common Misconceptions

• “A Gabor filter is an edge detector.” Only the odd filter is edge-like; the even filter responds to bars, and the complex filter’s modulus responds to gratings of the preferred frequency regardless of phase. All three are band-pass, not derivative, filters.
• “Any windowed sinusoid attains the uncertainty bound.” Only a Gaussian window with a complex carrier does; the example shows a Hann window 7% above the bound and truncation of the Gaussian moving it above too.
• “The real Gabor filter is optimally localized.” The bound concerns the complex function. A real filter’s spectrum has two lobes, and the product of spreads measured about its mean is far larger (5.09 times the bound in the example).
• “Gabor filters ignore constant intensity.” The even filter has a nonzero DC response that becomes large at wide bandwidths; only the odd filter, a mean-corrected even filter or a log-Gabor filter is exactly zero-mean.
• “Choosing the strongest filter gives the orientation.” It gives the orientation to within half the spacing of the bank; interpolating the energies with a doubled-angle average is far more accurate on clean data, though, as the example shows, not necessarily under heavy noise.
• “Library Gabor kernels are interchangeable.” scikit-image returns a normalized complex kernel truncated by default at n_stds=3 and parameterized by frequency and bandwidth; OpenCV returns an unnormalized real kernel parameterized by wavelength, $\sigma$ and aspect ratio, with an odd default phase. The same nominal parameters give different filters.

Where It Is Used

• Texture. Energies of a Gabor bank, pooled over regions, describe textures by their frequency and orientation content, and their variation across an image can drive segmentation (Daugman, 2021, slide 104). Texture analysis compares them with co-occurrence statistics and local binary patterns.
• Iris recognition. Daugman’s iris codes quantize the phase of complex 2D Gabor wavelets in polar coordinates to two bits per coefficient; the Biometrics article covers the codes and their matching.
• Scene recognition. The GIST descriptor pools the magnitudes of a bank of Gabor-like filters over a coarse grid; Scene Recognition describes the filters and their use.
• Faces and expressions. Gabor magnitudes at several scales and orientations are classical facial features; facial expression recognition lists a bank of 40 Gabor filters with an SVM among the baselines of RAF-DB.
• Pooled energy features. The fine-grained image classification example pools Gabor energy maps over the image as hand-designed features and shows what global pooling loses.
• Neural networks. The first-layer filters learned by convolutional networks resemble Gabor functions, and Gabor layers build the resemblance in; the receptive field and feature visualization articles discuss what early units respond to.