{ "cells": [ { "cell_type": "markdown", "id": "10916a29", "metadata": {}, "source": [ "# Lesson 16: Wavelets and Gabor Filters\n", "\n", "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*." ] }, { "cell_type": "code", "execution_count": null, "id": "f67d5f53", "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import cv2\n", "import pywt\n", "import matplotlib.pyplot as plt" ] }, { "cell_type": "markdown", "id": "b46b8356", "metadata": {}, "source": [ "## The Haar wavelet: the simplest possible wavelet\n", "\n", "The Haar wavelet transform splits a signal into an **approximation** (local averages) and a **detail** (local differences), computed on non-overlapping pairs of samples:\n", "\n", "$$a_k = \\frac{x_{2k} + x_{2k+1}}{\\sqrt{2}}, \\qquad d_k = \\frac{x_{2k} - x_{2k+1}}{\\sqrt{2}}$$\n", "\n", "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." ] }, { "cell_type": "code", "execution_count": null, "id": "74d59e31", "metadata": {}, "outputs": [], "source": [ "signal = np.array([1, 3, 5, 11, 7, 9, 2, 4], dtype=np.float64)\n", "\n", "approx = (signal[0::2] + signal[1::2]) / np.sqrt(2)\n", "detail = (signal[0::2] - signal[1::2]) / np.sqrt(2)\n", "\n", "cA_ref, cD_ref = pywt.dwt(signal, 'haar')\n", "print('manual matches pywt.dwt?', np.allclose(approx, cA_ref) and np.allclose(detail, cD_ref))\n", "print('approximation (cA):', np.round(approx, 2))\n", "print('detail (cD): ', np.round(detail, 2))" ] }, { "cell_type": "markdown", "id": "e9dc0c5b", "metadata": {}, "source": [ "## From Haar to Daubechies: smoother wavelets\n", "\n", "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." ] }, { "cell_type": "code", "execution_count": null, "id": "22efbe98", "metadata": {}, "outputs": [], "source": [ "db2 = pywt.Wavelet('db2')\n", "print('db2 low-pass (scaling) filter :', np.round(db2.dec_lo, 4))\n", "print('db2 high-pass (wavelet) filter:', np.round(db2.dec_hi, 4))\n", "\n", "fig, axes = plt.subplots(1, 2, figsize=(8, 3))\n", "axes[0].stem(db2.dec_lo)\n", "axes[0].set_title('db2 scaling filter (low-pass)')\n", "axes[1].stem(db2.dec_hi)\n", "axes[1].set_title('db2 wavelet filter (high-pass)')\n", "plt.tight_layout()\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "ffda9f02", "metadata": {}, "source": [ "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`." ] }, { "cell_type": "code", "execution_count": null, "id": "86303224", "metadata": {}, "outputs": [], "source": [ "def dwt_periodic(x, filt):\n", " \"\"\"One level of a periodic discrete wavelet transform with an arbitrary filter.\"\"\"\n", " n, L = len(x), len(filt)\n", " k = np.arange(n // 2)[:, None]\n", " i = np.arange(L)[None, :]\n", " idx = (2 * k + 2 - i) % n\n", " return (filt[None, :] * x[idx]).sum(axis=1)\n", "\n", "longer_signal = np.array([1, 3, 5, 11, 7, 9, 2, 4, 6, 8], dtype=np.float64)\n", "\n", "cA_mine = dwt_periodic(longer_signal, np.array(db2.dec_lo))\n", "cD_mine = dwt_periodic(longer_signal, np.array(db2.dec_hi))\n", "cA_ref, cD_ref = pywt.dwt(longer_signal, db2, mode='periodization')\n", "\n", "print('matches pywt (periodization mode)?', np.allclose(cA_mine, cA_ref) and np.allclose(cD_mine, cD_ref))\n", "\n", "reconstructed = pywt.idwt(cA_mine, cD_mine, db2, mode='periodization')\n", "print('perfect reconstruction?', np.allclose(reconstructed, longer_signal))" ] }, { "cell_type": "markdown", "id": "84f23f4e", "metadata": {}, "source": [ "## 2D wavelet decomposition of an image\n", "\n", "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**:\n", "\n", "- **LL**: low-pass both directions — a coarser, half-resolution copy of the image (like one Gaussian pyramid level)\n", "- **LH**: low-pass rows, high-pass columns — responds to horizontal edges\n", "- **HL**: high-pass rows, low-pass columns — responds to vertical edges\n", "- **HH**: high-pass both directions — responds to diagonal detail and corners" ] }, { "cell_type": "code", "execution_count": null, "id": "0e4836d3", "metadata": {}, "outputs": [], "source": [ "img = np.zeros((200, 200), dtype=np.float64)\n", "cv2.rectangle(img, (40, 40), (160, 160), 200, -1)\n", "cv2.circle(img, (100, 100), 40, 100, -1)\n", "\n", "LL, (LH, HL, HH) = pywt.dwt2(img, 'db2', mode='periodization')\n", "\n", "fig, axes = plt.subplots(1, 5, figsize=(14, 3))\n", "for ax, im, title in zip(axes, [img, LL, LH, HL, HH],\n", " ['Original', 'LL (approx.)', 'LH (horiz. edges)', 'HL (vert. edges)', 'HH (diagonal)']):\n", " ax.imshow(im, cmap='gray')\n", " ax.set_title(title, fontsize=9)\n", " ax.axis('off')\n", "plt.tight_layout()\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "b979a113", "metadata": {}, "source": [ "Just like the Laplacian pyramid, this decomposition is exactly invertible: `pywt.idwt2` reconstructs the original from the four subbands with no loss." ] }, { "cell_type": "code", "execution_count": null, "id": "73ef2ada", "metadata": {}, "outputs": [], "source": [ "reconstructed_img = pywt.idwt2((LL, (LH, HL, HH)), 'db2', mode='periodization')\n", "print('max reconstruction error:', np.abs(reconstructed_img - img).max())" ] }, { "cell_type": "markdown", "id": "c5b88fac", "metadata": {}, "source": [ "## Gabor filters: localized oscillations tuned to orientation\n", "\n", "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.\n", "\n", "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." ] }, { "cell_type": "code", "execution_count": null, "id": "eabcd447", "metadata": {}, "outputs": [], "source": [ "orientations_deg = [0, 45, 90, 135]\n", "wavelengths = [4, 8, 12, 16]\n", "\n", "fig, axes = plt.subplots(len(wavelengths), len(orientations_deg), figsize=(11, 11))\n", "for i, lambd in enumerate(wavelengths):\n", " for j, t in enumerate(orientations_deg):\n", " kernel = cv2.getGaborKernel((25, 25), sigma=3, theta=np.radians(t), lambd=lambd, gamma=0.5, psi=0)\n", " ax = axes[i, j]\n", " ax.imshow(kernel, cmap='gray')\n", " ax.set_xticks([])\n", " ax.set_yticks([])\n", " if i == 0:\n", " ax.set_title(fr'$\\theta$ = {t} deg', fontsize=9)\n", " if j == 0:\n", " ax.set_ylabel(fr'$\\lambda$ = {lambd}', fontsize=9)\n", "plt.tight_layout()\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "815d1eb5", "metadata": {}, "source": [ "### Orientation selectivity, demonstrated quantitatively\n", "\n", "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." ] }, { "cell_type": "code", "execution_count": null, "id": "b6dc82be", "metadata": {}, "outputs": [], "source": [ "def draw_bar(image, center, angle_deg, length=60, thickness=6, value=255):\n", " angle = np.radians(angle_deg)\n", " dx, dy = length / 2 * np.cos(angle), length / 2 * np.sin(angle)\n", " p1 = (int(center[0] - dx), int(center[1] - dy))\n", " p2 = (int(center[0] + dx), int(center[1] + dy))\n", " cv2.line(image, p1, p2, value, thickness)\n", "\n", "bar_orientations = [0, 45, 90, 135]\n", "centers = [(50, 50), (150, 50), (50, 150), (150, 150)]\n", "\n", "bars_img = np.zeros((200, 200), dtype=np.float64)\n", "for c, ang in zip(centers, bar_orientations):\n", " draw_bar(bars_img, c, ang)\n", "\n", "plt.imshow(bars_img, cmap='gray')\n", "plt.title('Bars at 0, 45, 90, 135 degrees')\n", "plt.axis('off')\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": null, "id": "8b970790", "metadata": {}, "outputs": [], "source": [ "responses = np.zeros((4, 4))\n", "for fi, bar_angle in enumerate(bar_orientations):\n", " # Note: cv2's theta parameter is the orientation of the sinusoidal stripes inside the\n", " # kernel, which runs *perpendicular* to the bar it responds to -- hence the +90 offset.\n", " kernel = cv2.getGaborKernel((25, 25), sigma=4, theta=np.radians(bar_angle + 90), lambd=10, gamma=0.5, psi=0)\n", " response = cv2.filter2D(bars_img, cv2.CV_64F, kernel)\n", " for ci, c in enumerate(centers):\n", " region = response[c[1] - 20:c[1] + 20, c[0] - 20:c[0] + 20]\n", " responses[fi, ci] = np.abs(region).mean()\n", "\n", "print(f'{\"filter tuned for\":>18}', ' '.join(f'{a:>6}' for a in bar_orientations), ' <- bar orientation')\n", "for fi, ang in enumerate(bar_orientations):\n", " print(f'{ang:>18}', ' '.join(f'{v:>6.0f}' for v in responses[fi]))\n", "\n", "best_match = [bar_orientations[i] for i in responses.argmax(axis=1)]\n", "print('\\neach filter peaks at its own bar orientation?', best_match == bar_orientations)" ] }, { "cell_type": "markdown", "id": "6501acd2", "metadata": {}, "source": [ "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." ] }, { "cell_type": "markdown", "id": "1cfa437b", "metadata": {}, "source": [ "### Exercise\n", "\n", "1. 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?\n", "2. 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.\n", "3. 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?" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "name": "python", "version": "3.x" } }, "nbformat": 4, "nbformat_minor": 5 }