From Pixels to Matches: A Visual Journey Through Classical Computer Vision

1. The Convolution Machine

An image is a grid of numbers. But a single pixel, alone, tells you almost nothing — is it part of a sky, an edge, a face? The interesting information in an image — edges, textures, objects — emerges only when you consider pixels together with their neighbours.

Convolution is the fundamental operation that turns a small grid of weights into neighbourhood updates on an image: sliding that grid can blur, sharpen, or detect edges — and this one pattern underpins almost everything in classical computer vision.

1.1 Pixels Have Neighbours

Pick any pixel in an image. Now look at the pixels immediately surrounding it — the 3×3 patch centred on your chosen pixel. This patch is called the pixel's neighbourhood.

Figure code: code/ch01_convolution.py.

A 3×3 neighbourhood is highlighted on a real image; the zoomed inset shows the raw intensity values (0–255) that the kernel will operate on.

The figure above shows a real image with one such neighbourhood highlighted. On the right, you can see the actual intensity values — just numbers, typically 0–255 for an 8-bit grayscale image. Individually, these numbers are unremarkable. But their pattern — how they change from one pixel to the next — encodes structure.

Key insight: almost all useful image operations work by examining neighbourhoods, not individual pixels. Blurring, sharpening, edge detection — they all ask the same question: "What's happening in this local patch?"

Why neighbourhoods? A single pixel has no spatial context. You can't tell whether pixel value 128 is part of a smooth gradient, a sharp edge, or random noise. But if you look at a 3×3 patch and see [40, 42, 41; 39, 128, 43; 40, 41, 42], that large jump at the centre screams "edge!" The neighbourhood provides the context that a lone pixel cannot.

1.2 Convolution: The Sliding Window

So we want to combine a pixel with its neighbours. But how? The answer is beautifully simple: use a kernel (also called a filter) — a small matrix of weights — and slide it across the image. At each position:

1. Overlay the kernel on the image patch.
2. Multiply each kernel weight by the corresponding pixel value.
3. Sum all the products.
4. Write the result to the output image.

This operation is called convolution, and the animation below shows it in action:

Figure code: code/ch01_convolution.py.

The kernel (orange grid) slides one pixel at a time. At each position, element-wise products of kernel weights and covered pixel values are summed to produce one output pixel.

Watch carefully: the kernel visits every pixel in the image (except the borders, for now). At each stop, it computes a single number — a weighted combination of the local neighbourhood — and that number becomes one pixel in the output.

Formally, for a kernel $K$ of size $(2a+1) \times (2b+1)$ applied to image $I$, the output at position $(x, y)$ is:

$$ O(x, y) = \sum_{i=-a}^{a} \sum_{j=-b}^{b} K(i, j) \cdot I(x - i, , y - j) $$

For the common 3×3 case ($a = b = 1$), this sums up $3 \times 3 = 9$ products — one for each cell in the neighbourhood.

Writing the sliding window out in plain NumPy makes the mechanism completely transparent:

import numpy as np

def convolve_naive(image, kernel):
    """Manual convolution: for each pixel, extract the neighbourhood,
    multiply element-wise by the kernel, and sum."""
    pad = kernel.shape[0] // 2
    img = image.astype(np.float64)
    out = np.zeros_like(img)
    for y in range(pad, img.shape[0] - pad):
        for x in range(pad, img.shape[1] - pad):
            window = img[y - pad:y + pad + 1, x - pad:x + pad + 1]
            out[y, x] = np.sum(window * kernel)  # the heart of convolution
    return out

In code, you don't need to write the loops yourself. OpenCV's filter2D does it all:

import cv2
import numpy as np

kernel = np.array([[0, -1, 0],
                   [-1, 5, -1],
                   [0, -1, 0]], dtype=np.float64)

output = cv2.filter2D(image, cv2.CV_64F, kernel)

Correlation vs. Convolution — the flip. You may have noticed the minus signs in $I(x - i, y - j)$. True convolution flips the kernel (rotates it 180°) before sliding. Correlation is the same operation without the flip: $\sum K(i,j) \cdot I(x+i, y+j)$. For symmetric kernels (like a Gaussian or box blur), flipping changes nothing — correlation and convolution give identical results. For asymmetric kernels (like Sobel), they differ. Most deep learning frameworks actually implement correlation and call it "convolution" — a harmless abuse of notation, as long as the network learns the right weights either way.

Why does convolution flip? The flip is what makes convolution commutative and associative — two properties that correlation lacks. Section 1.4 explains why those properties let you merge kernels and cut down full-image passes.

Let's work through a concrete example. Consider the 5×5 image patch $I$ and the 3×3 kernel $K$:

$$ I = \begin{bmatrix} 10 & 20 & 30 & 40 & 50 \ 20 & 40 & 60 & 80 & 100 \ 30 & 60 & 90 & 120 & 150 \ 40 & 80 & 120 & 160 & 200 \ 50 & 100 & 150 & 200 & 250 \end{bmatrix}, \quad K = \begin{bmatrix} 0 & -1 & 0 \ -1 & 5 & -1 \ 0 & -1 & 0 \end{bmatrix} $$

To compute the output at position $(1, 1)$ (0-indexed), we overlay $K$ on the 3×3 neighbourhood centred at $(1, 1)$:

$$ \text{patch} = \begin{bmatrix} 10 & 20 & 30 \ 20 & 40 & 60 \ 30 & 60 & 90 \end{bmatrix} $$

Element-wise multiply and sum (using correlation for simplicity — the kernel is symmetric under 180° rotation about its non-zero pattern):

$$ O(1,1) = 0 \cdot 10 + (-1) \cdot 20 + 0 \cdot 30 + (-1) \cdot 20 + 5 \cdot 40 + (-1) \cdot 60 + 0 \cdot 30 + (-1) \cdot 60 + 0 \cdot 90 $$

$$ = 0 - 20 + 0 - 20 + 200 - 60 + 0 - 60 + 0 = 40 $$

The output pixel is 40. The sharpening kernel has amplified the centre relative to its neighbours.

1.3 A Zoo of Kernels

The magic of convolution is that different kernels produce wildly different effects — all using the same sliding-window-and-sum mechanism. Here are six classic kernels and what they do to an image:

Figure code: code/ch01_convolution.py.

Six classic 3×3 kernels and their effects. Each column shows the kernel weights as a colour-coded heatmap, the input image, and the convolution output.

Let's meet each one:

In Python, these are simply NumPy arrays:

KERNELS = {
    "Identity":    np.array([[0, 0, 0], [0, 1, 0], [0, 0, 0]]),
    "Box Blur":    np.ones((3, 3)) / 9,
    "Gaussian 3x3": np.array([[1, 2, 1], [2, 4, 2], [1, 2, 1]]) / 16,
    "Sharpen":     np.array([[0, -1, 0], [-1, 5, -1], [0, -1, 0]]),
    "Sobel X":     np.array([[-1, 0, 1], [-2, 0, 2], [-1, 0, 1]]),
    "Emboss":      np.array([[-2, -1, 0], [-1, 1, 1], [0, 1, 2]]),
}

for name, kernel in KERNELS.items():
    result = cv2.filter2D(image, cv2.CV_64F, kernel)

Why do these kernels work? The secret lies in their weight patterns:

• Box blur assigns equal weight to all 9 neighbours. This averages out local variations — noise gets smoothed, but so do edges. It's the simplest low-pass filter.
• Gaussian blur assigns higher weight to the centre and less to corners. This produces a more natural-looking blur because nearby pixels are more relevant than distant ones. Gaussians have a distinctive set of properties in space, frequency, and separability — a key reason they replace box filters in practice.
• Sharpen has a centre weight of 5 and four neighbours at −1. The weights sum to $5 - 4 = 1$, so in flat regions (where all pixels are equal), the output equals the input. But at an edge, where the centre differs from its neighbours, the difference is amplified. This is the "positive centre, negative surround" pattern — it enhances contrast.
• Sobel X has positive weights on the right column and negative on the left. It computes (roughly) the horizontal derivative: if the right side is brighter than the left, the output is positive; if darker, negative. The middle row has weight ±2 for a slight smoothing effect — a weighted average along the vertical direction before differencing horizontally, which suppresses noise.
• Emboss has negative weights in the top-left and positive in the bottom-right, creating an asymmetric response that simulates light coming from one direction.

1.4 Properties of Convolution

Convolution isn't just a useful operation — it has elegant mathematical properties that make it a powerful tool for building complex image processing pipelines.

Commutativity: The order doesn't matter.

$$f * g = g * f$$

This means convolving kernel $A$ with kernel $B$ gives the same result as convolving $B$ with $A$.

Associativity: You can group operations however you like.

$$(f * g) * h = f * (g * h)$$

This is the big one for practical use. Suppose you want to blur an image and then sharpen it. You could apply blur to the image, then apply sharpen to the result — two full passes over the image. Or you could first convolve the blur kernel with the sharpen kernel to get a single combined kernel, then make one pass over the image. Same result, half the computation.

Distributivity over addition:

$$f * (g + h) = f * g + f * h$$

The practical payoff: In a pipeline with multiple filtering stages, associativity lets you pre-combine all the kernels into one. Instead of $N$ full-image convolutions, you do $N-1$ small kernel-kernel convolutions (cheap) plus one image convolution. This is a real optimisation used in production systems.

The figure below demonstrates associativity visually:

Figure code: code/ch01_convolution.py.

Associativity: blurring then sharpening (left path) gives the same result as applying the single pre-composed blur⊛sharpen kernel (right path).

On the left, we first blur the image, then sharpen the result. On the right, we first convolve the blur kernel with the sharpen kernel to get a combined kernel, then apply it to the image in a single pass. The results are identical — pixel for pixel.

You can verify this directly in code:

import numpy as np
import cv2
from scipy.signal import convolve2d

k_blur  = np.ones((3, 3)) / 9                                            # box blur kernel
k_sharp = np.array([[0, -1, 0], [-1, 5, -1], [0, -1, 0]], dtype=float)  # sharpen kernel

# Path A: combine kernels first, then apply once
combined = convolve2d(k_blur, k_sharp, mode="full")   # 3×3 ∗ 3×3 → 5×5 combined kernel
result_a = cv2.filter2D(image, cv2.CV_64F, combined)

# Path B: apply kernels sequentially
tmp      = cv2.filter2D(image, cv2.CV_64F, k_sharp)
result_b = cv2.filter2D(tmp,   cv2.CV_64F, k_blur)

# Verify associativity: maximum difference should be essentially zero
print(np.abs(result_a - result_b).max())  # ≈ 0 (floating-point rounding only)

Recall the definitions. For 1D signals $f$ and $g$:

Convolution:

$$ (f * g)(x) = \sum_{k} f(k) , g(x - k) $$

Correlation:

$$ (f \star g)(x) = \sum_{k} f(k) , g(x + k) $$

The only difference is the sign: $x - k$ vs $x + k$. Let's check commutativity.

Convolution:

$$ (f * g)(x) = \sum_{k} f(k) , g(x - k) $$

Substitute $m = x - k$ (so $k = x - m$):

$$ = \sum_{m} f(x - m) , g(m) = (g * f)(x) \quad \checkmark $$

Correlation:

$$ (f \star g)(x) = \sum_{k} f(k) , g(x + k) $$

Substitute $m = x + k$ (so $k = m - x$):

$$ = \sum_{m} f(m - x) , g(m) $$

This equals $(g \star f)(-x)$, not $(g \star f)(x)$. Correlation is not commutative — swapping $f$ and $g$ introduces a spatial reversal. This is why the "flip" in convolution isn't just a historical quirk; it's what makes the algebra work.

Associativity proof sketch. Using the Fourier convolution theorem: $\mathcal{F}{f * g} = \mathcal{F}{f} \cdot \mathcal{F}{g}$. Since pointwise multiplication of Fourier transforms is associative (it's just multiplication of complex numbers), convolution inherits associativity:

$$ (f * g) * h ;;\xleftrightarrow{\mathcal{F}};; (\hat{f}\hat{g})\hat{h} = \hat{f}(\hat{g}\hat{h}) ;;\xleftrightarrow{\mathcal{F}^{-1}};; f * (g * h) \quad \checkmark $$

Distributivity follows similarly: $\hat{f}(\hat{g} + \hat{h}) = \hat{f}\hat{g} + \hat{f}\hat{h}$.

1.5 Practical Notes: Output Size, Padding, and Stride

When you slide a 3×3 kernel across a $W \times H$ image, the kernel can't be centred on the border pixels (it would hang off the edge). So the output is smaller than the input. In general:

$$ \text{Output size} = \frac{W - K + 2P}{S} + 1 $$

where:

• $W$ = input width (or height)
• $K$ = kernel size
• $P$ = padding (how many extra pixels we add around the border)
• $S$ = stride (how many pixels the kernel moves per step)

For the default case ($K = 3$, $P = 0$, $S = 1$): the output is $W - 2$ pixels wide. Each border loses 1 pixel.

Padding adds extra pixels around the image border so that the kernel can be centred on edge pixels too. With $P = 1$ (one pixel of padding on each side), the output size equals the input size — this is called "same" padding. The most common padding strategies are:

• Zero padding: fill with 0s (default in most frameworks)
• Replicate: repeat the border pixel values
• Reflect: mirror the image at the border

Stride controls how far the kernel jumps between positions. With stride $S = 1$, it moves one pixel at a time (the default). With $S = 2$, it skips every other pixel, producing an output roughly half the size. Stride > 1 is rarely used in classical filtering but becomes important in convolutional neural networks, where it serves as a form of downsampling.

1.6 Key Takeaways

• An image is a grid of numbers. Useful information emerges from relationships between neighbouring pixels, not from individual values.
• Convolution = sliding window + weighted sum. A kernel slides over the image; at each position, element-wise multiply and sum produces one output pixel.
• Different kernels, different effects. The same convolution mechanism can blur (uniform weights), sharpen (positive centre, negative surround), or detect edges (positive on one side, negative on the other) — all by changing the kernel weights.
• Convolution is commutative and associative. This means you can combine multiple kernels into one before applying to the image — a powerful optimisation.
• Output size depends on kernel size, padding, and stride: $(W - K + 2P)/S + 1$.
• Correlation vs. convolution: convolution flips the kernel; for symmetric kernels, they're identical. The flip gives convolution its nice algebraic properties.

2. The Gaussian and the Idea of Scale

A box blur — equal weights everywhere — can smooth an image. It's simple, it works, and it's easy to understand. So why does virtually every blurring operation in computer vision use something different?

The answer is the Gaussian blur, and it isn't just "another kernel." It has a unique combination of mathematical properties that make it the natural choice for smoothing. Scale is explicit in a single parameter σ: it sets how far the filter looks — and therefore what level of detail survives the blur.

You will understand why the Gaussian kernel is the default choice in both theory and engineering, what σ really means, how to exploit the Gaussian's structure for speed, and how to handle image borders.

2.1 The Gaussian Function

What shape should the weights in a blur kernel have? The box blur says "everyone counts equally." But that's a strange assumption — should a pixel 5 steps away really matter as much as the pixel right next door?

A more natural idea: nearby pixels should contribute more, and far-away pixels should contribute less. The weight should decay smoothly with distance. The function that does this most elegantly is the Gaussian (the bell curve).

In one dimension, the Gaussian is:

$$ G(x) = \frac{1}{\sqrt{2\pi},\sigma} , e^{-x^2 / 2\sigma^2} $$

The parameter σ (sigma) controls the width of the bell. A small σ gives a narrow, peaked curve; a large σ gives a wide, gentle one. The key property: the weights are highest at the centre and fall off symmetrically — exactly the "nearby pixels matter more" intuition we wanted.

For images, we need a 2D Gaussian. Fortunately, the extension is straightforward — just use the distance from the centre in both directions:

$$ G(x, y) = \frac{1}{2\pi\sigma^2} , e^{-(x^2 + y^2) / 2\sigma^2} $$

Translating this formula directly into code gives the weight matrix:

import numpy as np

def gaussian_1d(x, sigma):
    """1D Gaussian: G(x) = exp(−x²/2σ²) / (√(2π) σ)"""
    return np.exp(-x**2 / (2 * sigma**2)) / (np.sqrt(2 * np.pi) * sigma)

def gaussian_2d(size, sigma):
    """Discrete 2D Gaussian kernel, normalised so all weights sum to 1."""
    ax = np.arange(size) - size // 2   # integer axis centred at 0
    xx, yy = np.meshgrid(ax, ax)
    kernel = np.exp(-(xx**2 + yy**2) / (2 * sigma**2))
    return kernel / kernel.sum()        # normalise: preserve average brightness

# Example: 13×13 kernel for σ=2 (covers ±3σ)
kernel = gaussian_2d(size=13, sigma=2.0)
print(kernel.sum())   # 1.0 — unit gain

The figure below shows all three views: the 1D bell curve, the 2D surface (a smooth dome), and the 2D heatmap (bright at centre, dark at edges):

Figure code: code/ch02_gaussian.py.

The Gaussian function in 1D (left: bell curve with σ and 3σ marked), as a 3D surface (centre), and as a 2D kernel heatmap (right). Brighter centre = higher weight; weights taper smoothly to zero.

Core takeaway: the Gaussian kernel assigns the most weight to the centre pixel and smoothly decreasing weight to pixels farther away. This produces the most natural-looking blur — no abrupt transitions, no blocky artifacts.

Why Gaussian and not some other bell shape? There are many functions that peak at the centre and decay outward — a triangle, a cosine lobe, a Cauchy distribution. The Gaussian is special for several reasons:

1. It's separable — a 2D Gaussian can be decomposed into two 1D passes (Section 2.4), cutting computation dramatically.
2. It's the only kernel that doesn't introduce new frequencies — blurring with a Gaussian is a clean low-pass filter in the frequency domain.
3. Convolving two Gaussians gives another Gaussian: $G_{\sigma_1} * G_{\sigma_2} = G_{\sqrt{\sigma_1^2 + \sigma_2^2}}$. No other kernel family has this closure property. This becomes critical for multi-scale image representations such as pyramids.
4. The Central Limit Theorem: if you repeatedly convolve any kernel with itself, the result converges to a Gaussian. The Gaussian is the "attractor" of convolution.

In short, the Gaussian isn't just a reasonable choice — it's mathematically the choice for smoothing.

The continuous 2D Gaussian centred at the origin is:

$$ G(x, y) = \frac{1}{2\pi\sigma^2} , e^{-(x^2 + y^2) / 2\sigma^2} $$

The factor $\frac{1}{2\pi\sigma^2}$ is the normalisation constant — it ensures the function integrates to 1 over the entire plane:

$$ \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} G(x, y) , dx , dy = 1 $$

Proof. Because the exponent factors (see Section 2.4), the double integral separates:

$$ \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \frac{1}{2\pi\sigma^2} e^{-(x^2+y^2)/2\sigma^2} dx,dy = \left(\frac{1}{\sqrt{2\pi}\sigma} \int_{-\infty}^{\infty} e^{-x^2/2\sigma^2} dx\right)^2 $$

Each 1D integral is the well-known Gaussian integral, equal to $\sqrt{2\pi}\sigma$. So each factor evaluates to 1, and the product is $1 \times 1 = 1$. $\checkmark$

Discrete approximation. In practice, we sample $G(x, y)$ at integer coordinates within a finite window, then renormalise so the discrete weights sum to 1. For a $(2r+1) \times (2r+1)$ kernel with $r = \lceil 3\sigma \rceil$:

$$ K[i, j] = \frac{G(i, j)}{\sum_{m=-r}^{r}\sum_{n=-r}^{r} G(m, n)}, \quad -r \leq i, j \leq r $$

The renormalisation compensates for the truncated tails.

2.2 σ and the Idea of Scale

The Gaussian has one free parameter: σ (sigma), the standard deviation. But σ isn't just a number you tweak — it encodes a profound idea.

σ controls the scale at which you observe the image. A small σ preserves fine details — individual hairs, texture grain, thin edges. A large σ washes those away, leaving only large-scale structures — the shape of a face, the outline of a building, the horizon.

Watch what happens as σ increases from 0.5 to 8:

Figure code: code/ch02_gaussian.py.

σ increases from 0.5 to 8. Left: the kernel heatmap grows with σ. Right: the blurred result loses progressively finer detail as σ increases.

As σ grows, small features vanish first (texture, noise), then medium features (thin edges, small objects), until only the coarsest structures remain. This is the idea of scale: different values of σ reveal different levels of structure in the same image.

The static comparison below makes it easier to compare directly:

Figure code: code/ch02_gaussian.py.

Static comparison at four σ values. Kernel size (annotated above each panel) grows with σ; fine textures vanish at higher σ while coarse structure is preserved.

At σ = 0.5, the image looks nearly identical to the original — the kernel is tiny and barely averages anything. By σ = 4, fine textures are completely gone, but the overall composition — where the bright regions are, the rough shape of objects — is perfectly clear.

Core takeaway: σ is not just "amount of blur." It's the scale at which you choose to examine the world. This idea stays central for multi-scale image representations (pyramids) and for detecting and describing structure across scales.

The 3σ Rule and Kernel Size

How big should the kernel be for a given σ? The Gaussian never truly reaches zero — it extends to infinity. But it decays very fast. The 3σ rule from statistics tells us:

• Within ±1σ: 68.3% of the total weight
• Within ±2σ: 95.4% of the total weight
• Within ±3σ: 99.7% of the total weight

Beyond 3σ, the weights are negligible (< 0.3% of the total). So we truncate the kernel at radius $r = \lceil 3\sigma \rceil$, giving a kernel of size:

$$ K = 2 \lceil 3\sigma \rceil + 1 $$

For σ = 1, this gives $K = 7$. For σ = 2, $K = 13$. For σ = 4, $K = 25$.

In code:

import cv2
import numpy as np

sigma = 2.0
ksize = int(np.ceil(sigma * 6)) | 1  # round up to next odd number
blurred = cv2.GaussianBlur(img, (ksize, ksize), sigma)

The | 1 trick ensures an odd kernel size (required so the kernel has a well-defined centre pixel). OpenCV's GaussianBlur handles everything else — sampling the Gaussian, normalising the weights, and applying the filter.

σ as a frequency cutoff. There's another way to think about σ: in the frequency domain, a Gaussian with standard deviation σ in space corresponds to a Gaussian with standard deviation $1/(2\pi\sigma)$ in frequency. Larger σ in space → narrower passband in frequency → more high frequencies removed. This is why Gaussian blur is a low-pass filter: it suppresses high-frequency content (edges, noise, texture) and preserves low-frequency content (smooth gradients, large shapes).

Concretely, a Gaussian with σ = 2 in space attenuates frequencies above roughly $f = 1/(2\pi \times 2) \approx 0.08$ cycles/pixel — meaning structures smaller than about 12 pixels start getting significantly suppressed.

Very large σ → box blur? As σ → ∞, the Gaussian becomes increasingly flat within any finite kernel. If you also let the kernel grow without bound, the weights become nearly uniform — approaching a box blur. In practice, though, nobody uses σ = 100; at such scales, downsampling first is far more efficient.

The Fourier transform of a Gaussian is another Gaussian. For a 1D Gaussian $G(x) = \frac{1}{\sqrt{2\pi}\sigma} e^{-x^2/2\sigma^2}$, its Fourier transform is:

$$ \hat{G}(\omega) = e^{-2\pi^2 \sigma^2 \omega^2} $$

This is a Gaussian in frequency with standard deviation $\sigma_f = \frac{1}{2\pi\sigma}$.

Key relationships:

The inverse relationship $\sigma \cdot \sigma_f = \frac{1}{2\pi}$ is the spatial-frequency uncertainty principle: you cannot be simultaneously localised in both space and frequency. The Gaussian achieves the minimum uncertainty product — it's the optimal trade-off between spatial and frequency localisation.

Cascading Gaussians. If you convolve two Gaussians with parameters $\sigma_1$ and $\sigma_2$:

$$ G_{\sigma_1} * G_{\sigma_2} = G_{\sqrt{\sigma_1^2 + \sigma_2^2}} $$

Proof: In the frequency domain, convolution becomes multiplication: $e^{-2\pi^2\sigma_1^2\omega^2} \cdot e^{-2\pi^2\sigma_2^2\omega^2} = e^{-2\pi^2(\sigma_1^2 + \sigma_2^2)\omega^2}$, which is the Fourier transform of a Gaussian with $\sigma = \sqrt{\sigma_1^2 + \sigma_2^2}$.

This property is essential for building Gaussian pyramids: to go from scale $\sigma_1$ to scale $\sigma_2 > \sigma_1$, you only need to apply an incremental blur with $\sigma_\Delta = \sqrt{\sigma_2^2 - \sigma_1^2}$.

2.3 Separability: A Computational Gift

Here's a practical problem. As σ increases, the kernel size $K$ grows (recall $K \approx 6\sigma + 1$). For σ = 4, the kernel is 25×25 — that's 625 multiplications and additions per pixel. For a 1920×1080 image, that's over 1.2 billion multiply-adds for a single blur. Can we do better?

Yes — dramatically — thanks to a special property of the Gaussian: separability.

The 2D Gaussian factors into the product of two 1D Gaussians:

$$ G(x, y) = G(x) \cdot G(y) $$

This means the 2D kernel can be written as the outer product of a column vector and a row vector:

$$ \mathbf{K}_{2D} = \mathbf{g} \cdot \mathbf{g}^T $$

where $\mathbf{g}$ is the 1D Gaussian sampled as a column vector.

Figure code: code/ch02_gaussian.py.

A 2D Gaussian kernel equals the outer product of a column vector and a row vector. This allows two cheap 1D passes instead of one expensive 2D convolution, reducing cost from K² to 2K per pixel.

The payoff: instead of applying the full 2D kernel (cost: $K^2$ operations per pixel), you apply two 1D passes — first along rows, then along columns (or vice versa). The cost drops to $2K$ operations per pixel.

For σ = 8, separability gives a 24.5× speedup — the difference between real-time and sluggish.

In code, you can see this decomposition directly:

# Build 2D kernel from 1D via outer product
g1d = cv2.getGaussianKernel(ksize=13, sigma=2.0)  # column vector (13×1)
g2d = g1d @ g1d.T  # outer product → 13×13 kernel

# Separable application: row pass then column pass
row_pass = cv2.filter2D(img, cv2.CV_64F, g1d.T)  # 1×K row kernel
col_pass = cv2.filter2D(row_pass, cv2.CV_64F, g1d)  # K×1 column kernel
# col_pass ≈ cv2.filter2D(img, cv2.CV_64F, g2d)  — same result!

Core takeaway: the Gaussian is separable, turning an $O(K^2)$ operation into $O(2K)$. OpenCV's GaussianBlur automatically exploits this — you get the speedup for free.

Why does separability work? The 2D convolution with a separable kernel $K = \mathbf{a}\mathbf{b}^T$ can be rewritten as two nested 1D convolutions. First convolve each row of the image with $\mathbf{b}$ (the row vector), producing an intermediate image. Then convolve each column of that intermediate with $\mathbf{a}$ (the column vector). The result is identical to the full 2D convolution.

Visually, think of it this way: the 2D Gaussian assigns weight to each pixel based on its distance from the centre. But because the Gaussian's distance formula $x^2 + y^2$ separates into $x^2$ and $y^2$, the weight at position $(x, y)$ is just the product of the weight for displacement $x$ and the weight for displacement $y$. Horizontal distance and vertical distance contribute independently.

Not all kernels are separable. The Sobel filter is separable (it factors as a smoothing vector times a derivative vector), but most arbitrary kernels are not. Separability is a special algebraic property — another reason the Gaussian earns its place as the default blur.

The 2D Gaussian is:

$$ G(x, y) = \frac{1}{2\pi\sigma^2} , e^{-(x^2 + y^2) / 2\sigma^2} $$

The exponent factors:

$$ e^{-(x^2 + y^2) / 2\sigma^2} = e^{-x^2 / 2\sigma^2} \cdot e^{-y^2 / 2\sigma^2} $$

And the normalisation constant factors:

$$ \frac{1}{2\pi\sigma^2} = \frac{1}{\sqrt{2\pi}\sigma} \cdot \frac{1}{\sqrt{2\pi}\sigma} $$

Therefore:

$$ G(x, y) = \underbrace{\frac{1}{\sqrt{2\pi}\sigma} e^{-x^2/2\sigma^2}}{G_x(x)} \cdot \underbrace{\frac{1}{\sqrt{2\pi}\sigma} e^{-y^2/2\sigma^2}}{G_y(y)} $$

In matrix form, if $\mathbf{g} = [G(-r), \ldots, G(0), \ldots, G(r)]^T$ is the sampled 1D Gaussian (column vector), then the 2D kernel is:

$$ \mathbf{K} = \mathbf{g} , \mathbf{g}^T $$

This is a rank-1 matrix — the hallmark of separability.

Convolution consequence. For any image $I$:

$$ I * \mathbf{K} = I * (\mathbf{g},\mathbf{g}^T) = (I * \mathbf{g}) * \mathbf{g}^T $$

where the first convolution is along columns (with $\mathbf{g}$) and the second is along rows (with $\mathbf{g}^T$). The order can be swapped (by commutativity).

Operation count. For an $N$-pixel image with kernel size $K$:

• 2D convolution: $N \times K^2$ multiply-adds
• Separable: $N \times K$ (row pass) + $N \times K$ (column pass) = $2NK$
• Ratio: $K^2 / 2K = K/2$

For $K = 25$ (σ = 4), the separable version is 12.5× faster.

2.4 Boundary Handling

There's a detail we've been quietly ignoring: what happens when the kernel reaches the edge of the image?

When the kernel is centred on a border pixel, part of it hangs off the edge — there are no image pixels to multiply with. We need to invent values for these missing pixels. The choice matters more than you might think, especially for large kernels.

There are four common strategies:

Figure code: code/ch02_gaussian.py.

Four border-handling strategies with a large Gaussian blur: zero padding darkens borders; replicate smears edges; reflect (best default) preserves local structure; wrap suits tiling textures only.

Look at the edges carefully in the figure above. With zero padding, the borders darken — the kernel averages real pixels with black zeros, dragging the values down. Replicate avoids the darkening but can create a subtle "smeared" look along the edges. Reflect is usually the best default — the mirrored extension preserves local structure and produces the least visible artifacts. Wrap is appropriate only for images that truly tile (like textures), and produces strange discontinuities on normal photos.

import cv2

p = 20  # padding size

# Zero padding (default)
cv2.copyMakeBorder(img, p, p, p, p, cv2.BORDER_CONSTANT, value=0)

# Replicate border pixels
cv2.copyMakeBorder(img, p, p, p, p, cv2.BORDER_REPLICATE)

# Reflect at border
cv2.copyMakeBorder(img, p, p, p, p, cv2.BORDER_REFLECT)

# Wrap around
cv2.copyMakeBorder(img, p, p, p, p, cv2.BORDER_WRAP)

Core takeaway: image borders need padding, and the choice of padding method affects the result. Use BORDER_REFLECT as a safe default. Be aware of BORDER_CONSTANT (zero padding) when you see dark halos near edges — that's the padding talking, not a bug.

When does padding matter? For small σ (say, σ = 1, kernel size 7), the padding only affects pixels within 3 pixels of the edge — negligible for most applications. But for large σ (σ = 8, kernel size 49), the affected border region is 24 pixels wide. On a small image or a tightly cropped region, this can be a significant fraction of the image.

OpenCV's default. cv2.GaussianBlur uses BORDER_REFLECT_101 by default (reflection that doesn't duplicate the border pixel). This is almost always the right choice. You only need to think about padding when doing something unusual — like applying very large blurs, processing tile-based textures, or implementing frequency-domain filtering where wrap-around is assumed.

2.5 Gaussian vs. Box Blur

With the Gaussian kernel in hand, a direct comparison question is: why not just use a box blur?

The most important practical difference is visual quality. A box blur treats all pixels within the window equally, which can create visible "plateaus" and "steps" in smooth gradients. The Gaussian's smooth falloff avoids these artifacts entirely.

A fun fact about repeated box blurs. If you apply a box blur multiple times (say, 3–4 passes), the result converges to a Gaussian blur — this is the Central Limit Theorem in action. Some real-time graphics systems exploit this: three passes of a small box blur can be faster than one pass of a large Gaussian, while giving nearly indistinguishable results. This technique is called a stacked box blur or box blur cascade.

2.6 Core Takeaways

• The Gaussian kernel assigns the highest weight to the centre pixel and smoothly decreasing weights outward. This produces natural, artifact-free blurring.
• σ controls scale. Small σ preserves fine detail; large σ reveals only coarse structure. The kernel size should be approximately $6\sigma + 1$ (the 3σ rule).
• The Gaussian is separable. A 2D Gaussian is the outer product of two 1D Gaussians, reducing computation from $O(K^2)$ to $O(2K)$ per pixel.
• Boundary handling matters. Zero padding darkens edges; reflect padding is the safest default.
• Gaussian > box blur for almost all applications: smoother results, cleaner frequency response, and the same separability benefit.

Appendix: Reproducing All Figures

Every figure in this tutorial is generated from Python code. Download the self-contained package below — unzip, install dependencies, and run one command to regenerate all 9 figures from scratch:

Download code.zip (contains source code + input image, ~160 KB)

unzip code.zip -d code && cd code
pip install -r requirements.txt
python generate.py

See code/README.md for the full file listing and details.