Image Processing Fundamentals

A digital image is a discrete, finite approximation of a continuous 2D signal. Physically, a camera sensor divides the scene into a regular grid of photosensitive cells called pixels (picture elements). Each pixel records light intensity over a small spatial area, and the sensor's analog-to-digital converter maps that intensity to an integer in some fixed range.

Grayscale images are 2D arrays of integers: I[y][x] gives the brightness at row y, column x. For an 8-bit image, values range from 0 (black) to 255 (white). A 1920×1080 grayscale image is simply a matrix with about 2 million entries.

Color images add a third dimension. The most common representation is RGB: three separate 2D channels (red, green, blue), so the image is a 3D array of shape (H, W, 3). Other color spaces are often more useful for processing:

• HSV (Hue-Saturation-Value): separates chromatic content from brightness, making color-based segmentation easier.
• YCbCr: separates luminance (Y) from color difference (Cb, Cr); JPEG compression stores Y at full resolution and subsamples the color channels because the human eye is more sensitive to luminance variations.
• *LAB (CIE Lab)**: a perceptually uniform space where Euclidean distance approximates perceived color difference.

Sampling and quantization introduce two fundamental constraints. The spatial resolution (pixels per inch) determines the finest detail you can represent — this is the sampling rate. The bit depth determines how many distinct intensity levels are available — this is quantization. Both create information loss: aliasing from undersampling, and banding/posterization from coarse quantization.

In code, NumPy makes image arithmetic natural. OpenCV loads images as uint8 NumPy arrays by default (BGR channel order, not RGB — a persistent source of confusion).

import numpy as np
from PIL import Image

# Load as numpy array: shape (H, W, 3), dtype uint8
img = np.array(Image.open('photo.jpg'))   # RGB order

# Grayscale conversion (ITU-R BT.601 luma coefficients)
# Human vision weights green most, red second, blue least
gray = 0.299*img[:,:,0] + 0.587*img[:,:,1] + 0.114*img[:,:,2]
gray = gray.astype(np.uint8)  # shape (H, W)

# Normalise to [0,1] float for arithmetic
gray_f = gray.astype(np.float32) / 255.0

print(f'Resolution: {img.shape[1]}×{img.shape[0]}  Channels: {img.shape[2]}')

Convolution is the fundamental operation of image filtering. Intuitively, it is a sliding weighted average: you place a small weight matrix (called a kernel or filter) over each pixel, multiply the kernel weights by the corresponding pixel values, sum the products, and write the result at the center position. You then slide the kernel one pixel, repeat, and build up an output image of the same size.

Mathematically, the 2D discrete convolution of image I with kernel K is:

(I * K)[y,x] = Σ_m Σ_n  I[y-m, x-n] · K[m, n]

Note the flipping: the kernel is indexed by (−m, −n), meaning it is flipped both horizontally and vertically before the dot product. This is the technically correct definition that ensures commutativity and associativity. In practice, most image-processing frameworks implement cross-correlation instead (no flip), because filters are usually symmetric and the distinction disappears. NumPy's np.correlate and PyTorch's F.conv2d use cross-correlation.

The kernel encodes what feature you extract:

• A kernel of all equal weights → box blur (local average).
• A kernel shaped like a Gaussian → Gaussian blur (weighted average, smooth rolloff).
• A kernel with a positive center and negative surround → edge detection.
• A kernel with a large positive center and small negative ring → sharpening (unsharp mask).

Border handling is a practical choice. When the kernel overlaps the image boundary, you need values for out-of-bounds pixels. Options: zero-pad (cheapest), reflect the image, or replicate the edge pixel. Zero-padding can create dark halos; reflection is usually the best default.

Separability is a key optimization. If a 2D kernel can be written as the outer product of two 1D vectors — K = v·h^T — you can apply h horizontally and v vertically in two passes. This reduces complexity from O(K²) to O(2K) multiplications per pixel, where K is the kernel side length. Gaussian kernels are separable; Sobel kernels are separable.

import numpy as np
from scipy.ndimage import convolve

def convolve2d_naive(image, kernel):
    """Slow reference implementation for teaching."""
    H, W = image.shape
    kH, kW = kernel.shape
    pH, pW = kH // 2, kW // 2
    # Reflect-pad the image
    padded = np.pad(image, ((pH, pH), (pW, pW)), mode='reflect')
    out = np.zeros_like(image, dtype=np.float32)
    for y in range(H):
        for x in range(W):
            patch = padded[y:y+kH, x:x+kW]
            out[y, x] = np.sum(patch * kernel)
    return out

# Box blur kernel (3x3 average)
box = np.ones((3, 3), dtype=np.float32) / 9.0

# Apply to grayscale image (float32, values in [0,1])
blurred = convolve(gray_f, box, mode='reflect')

The Gaussian filter is the gold standard for smoothing and noise reduction. Its kernel weights are samples of a 2D Gaussian function:

G(x,y) = (1 / 2πσ²) · exp(−(x²+y²) / 2σ²)

The parameter σ (sigma) controls the spatial extent of the blur: larger σ means more smoothing. You truncate the kernel at ±3σ (beyond that the weight is negligible), giving a kernel of side length roughly 6σ+1.

Why is the Gaussian special? Several deep reasons:

1. Optimal time-frequency tradeoff: Among all smoothing functions, the Gaussian simultaneously minimizes spread in both spatial and frequency domains (it is the eigenfunction of the Fourier transform). No other kernel can be as compact in both domains at once.
2. The diffusion equation: Gaussian blurring satisfies the heat equation ∂I/∂t = ∇²I. Increasing σ is equivalent to letting heat diffuse in the image — a physical interpretation that makes the mathematics feel natural.
3. Scale space: A sequence of Gaussian-blurred images at increasing σ forms a Gaussian scale space, which is the mathematically principled way to analyze an image at multiple resolutions. The SIFT feature detector, which revolutionized computer vision in 2004, is built on Gaussian scale space.
4. Separability: G(x,y) = G(x)·G(y), so the 2D Gaussian is exactly separable into two 1D Gaussians. A 31×31 Gaussian blur (961 multiplies per pixel) becomes two 31×1 passes (62 multiplies). This makes large-σ blurs practical.

Noise model motivation: Salt-and-pepper noise is best removed with a median filter (nonlinear: replace each pixel with the median of its neighborhood), which preserves edges far better than Gaussian blur. Gaussian noise (from camera sensors in low light) is well-matched to Gaussian blur because the filter is optimal for signals corrupted by additive Gaussian noise in a least-squares sense.

import numpy as np

def gaussian_kernel_1d(sigma, truncate=3.0):
    """Build a 1D Gaussian kernel (for separable application)."""
    radius = int(truncate * sigma + 0.5)
    x = np.arange(-radius, radius + 1)
    kernel = np.exp(-x**2 / (2 * sigma**2))
    return kernel / kernel.sum()   # normalise so DC gain = 1

def gaussian_blur(image, sigma):
    from scipy.ndimage import convolve1d
    k = gaussian_kernel_1d(sigma)
    # Separable: apply horizontally then vertically
    tmp = convolve1d(image, k, axis=1, mode='reflect')
    return convolve1d(tmp,   k, axis=0, mode='reflect')

# sigma=1 → mild softening; sigma=5 → heavy blur
blur1 = gaussian_blur(gray_f, sigma=1.0)
blur5 = gaussian_blur(gray_f, sigma=5.0)

Sharpening is the complement of blurring: instead of attenuating details, we amplify them. The most principled technique is unsharp masking, a name that comes from the darkroom era of photography: photographers would make a blurry (unsharp) negative and subtract it from the original to isolate the fine details, then add those details back with enhanced strength.

Mathematically:

sharp = original + λ · (original − blurred)
       = (1+λ) · original − λ · blurred

where λ > 0 controls sharpening strength. The term (original − blurred) is called the detail signal or high-frequency residual — it contains edges, textures, and noise.

Rearranging: sharp = original + λ · detail. If λ = 1, this is equivalent to convolving with the kernel:

  0  −1   0
 −1   5  −1
  0  −1   0

This is the Laplacian sharpening kernel: the center adds four copies of the pixel and the cross subtracts neighbors, amplifying center-surround differences. A more aggressive version amplifies all diagonals too:

 −1  −1  −1
 −1   9  −1
 −1  −1  −1

A critical practical warning: sharpening amplifies noise as much as signal. Always apply some noise reduction (Gaussian blur) before sharpening if the input is noisy. The sharpening filter cannot tell the difference between a real edge and a noisy pixel; it indiscriminately amplifies both.

High-boost filtering generalizes unsharp masking by using a constant A instead of 1+λ:

high-boost = A · original − blurred

When A=1 this is an ordinary high-pass filter (edge detector). When A>1 it is a sharpener that retains the original image content (so the output still looks like an image, not just edges).

import numpy as np
from scipy.ndimage import gaussian_filter

def unsharp_mask(image, sigma=1.0, strength=1.5):
    """
    Sharpen by adding back a scaled high-frequency residual.
    strength=1.5 means we boost detail by 150%.
    """
    blurred = gaussian_filter(image, sigma=sigma)
    detail  = image - blurred          # high-frequency residual
    sharp   = image + strength * detail
    return np.clip(sharp, 0.0, 1.0)    # keep values valid

# Laplacian sharpening kernel (equivalent to strength=1)
laplacian_sharpen = np.array([
    [0, -1,  0],
    [-1, 5, -1],
    [0, -1,  0]
], dtype=np.float32)

from scipy.ndimage import convolve
sharp_lap = np.clip(convolve(gray_f, laplacian_sharpen, mode='reflect'), 0, 1)

Edges are the places in an image where intensity changes rapidly — the boundaries between objects, the silhouettes, the creases. Mathematically, edges correspond to large values of the image gradient.

For a continuous image I(x,y), the gradient is the vector ∇I = (∂I/∂x, ∂I/∂y). Its magnitude |∇I| is large at edges. In the discrete setting, we approximate these partial derivatives with finite differences, implemented as convolution with derivative kernels.

First-order (gradient) detectors:

The Sobel operator uses 3×3 kernels that combine a derivative estimate in one direction with a 1D Gaussian smoothing in the other direction:

Sx = [-1  0  1]     Sy = [-1 -2 -1]
     [-2  0  2]          [ 0  0  0]
     [-1  0  1]          [ 1  2  1]

Note Sx = [1 2 1]^T · [-1 0 1] (separable!). Apply both kernels, then:

magnitude = sqrt(Sx_response² + Sy_response²)
direction = atan2(Sy_response, Sx_response)   # angle in radians

Second-order (Laplacian) detectors:

The Laplacian ∇²I = ∂²I/∂x² + ∂²I/∂y² measures curvature. Edges appear as zero crossings of the Laplacian (where it transitions from positive to negative). The 3×3 Laplacian kernel is:

 0   1   0
 1  -4   1
 0   1   0

Because the Laplacian amplifies noise aggressively, in practice you first smooth with a Gaussian, then apply the Laplacian. This combination is called the Laplacian of Gaussian (LoG) or Marr-Hildreth operator. The combined filter can be precomputed as a single kernel. The LoG can be approximated by the Difference of Gaussians (DoG) — the difference of two Gaussian-blurred images at different scales — which is far cheaper to compute and underlies SIFT.

Canny edge detector (1986, widely considered the best classical edge detector):

1. Gaussian smooth with chosen σ.
2. Compute Sobel gradient magnitude and direction.
3. Non-maximum suppression (NMS): for each pixel, check whether its gradient magnitude is a local maximum along the gradient direction. If not, suppress it (set to 0). This thins edges to 1-pixel-wide curves.
4. Double thresholding + hysteresis: use a high threshold to find definite edges; use a low threshold to extend edges by linking adjacent pixels that exceed the low threshold. This gives connected edge curves and suppresses isolated noise responses.

Canny is an optimal edge detector in a formal sense: it maximizes signal-to-noise ratio, localizes edges accurately, and returns minimal spurious responses.

import numpy as np
from scipy.ndimage import convolve, gaussian_filter

def sobel_edges(image):
    Sx = np.array([[-1, 0, 1], [-2, 0, 2], [-1, 0, 1]], dtype=np.float32)
    Sy = np.array([[-1,-2,-1], [ 0, 0, 0], [ 1, 2, 1]], dtype=np.float32)
    gx = convolve(image, Sx, mode='reflect')
    gy = convolve(image, Sy, mode='reflect')
    magnitude  = np.sqrt(gx**2 + gy**2)
    direction  = np.arctan2(gy, gx)   # in radians, range [-π, π]
    return magnitude, direction

# For Canny, use OpenCV's well-optimised implementation:
import cv2
img_u8 = (gray_f * 255).astype(np.uint8)
edges_canny = cv2.Canny(img_u8,
                         threshold1=50,   # low hysteresis threshold
                         threshold2=150)  # high hysteresis threshold

mag, dirn = sobel_edges(gray_f)
# Threshold at 0.1 for a simple binary edge map
binary_edges = (mag > 0.1).astype(np.uint8) * 255

The Fourier transform answers a profound question: What mixture of sinusoids makes up this signal? For images, the signal is 2D, and the sinusoids are 2D spatial waves — patterns of light and dark stripes at various orientations and frequencies.

Why sinusoids? They are the eigenfunctions of the convolution operator. If you convolve a sinusoid with any filter, the output is the same sinusoid multiplied by a (complex) scalar. This means convolution, which is complex and global in the spatial domain, becomes pointwise multiplication in the frequency domain — dramatically simpler.

2D Discrete Fourier Transform (DFT): For an M×N image I[y,x], the DFT is:

F[v,u] = Σ_{y=0}^{M-1} Σ_{x=0}^{N-1}  I[y,x] · exp(−2πi(uy/M + vx/N))

F[v,u] is a complex number whose magnitude |F[v,u]| says how strongly that spatial frequency is present, and whose phase arg(F[v,u]) says where the pattern is aligned.

Interpreting the spectrum: The center of the Fourier spectrum (after fftshift) corresponds to zero frequency (DC) — the average brightness. Moving outward increases spatial frequency. High frequencies correspond to fine details, textures, and edges; low frequencies correspond to the coarse overall structure and gradual brightness variations.

A natural image has a characteristic 1/f spectrum: energy falls off as 1/f² approximately, meaning low frequencies dominate. Edges appear as high-frequency ridges perpendicular to their orientation in the spectrum.

The Convolution Theorem is the payoff:

fft2(I * K) = fft2(I) · fft2(K)

Convolution in space = multiplication in frequency. This lets you:

1. Apply huge kernels cheaply (O(MN log MN) via FFT vs. O(MNK²) in spatial domain).
2. Design filters directly by specifying what frequencies to keep or attenuate.
3. Understand what any kernel does by looking at its Fourier transform.

The Fourier transform of a Gaussian is a Gaussian — this is why Gaussian blur is the canonical low-pass filter. Its frequency-domain interpretation is clean: it is a smooth bell curve that attenuates high frequencies and leaves low frequencies intact, with no ringing artifacts.

import numpy as np
import matplotlib.pyplot as plt

def show_spectrum(image):
    """Compute and display log-magnitude Fourier spectrum."""
    F = np.fft.fft2(image)
    F_shifted = np.fft.fftshift(F)       # move DC to center
    magnitude = np.abs(F_shifted)
    # Log scale for visualisation (dynamic range is huge)
    log_mag = np.log1p(magnitude)
    return log_mag

# Applying a filter in the frequency domain
def ideal_low_pass(image, cutoff_fraction=0.1):
    """Zero out all frequencies above cutoff (fraction of Nyquist)."""
    F = np.fft.fft2(image)
    F_s = np.fft.fftshift(F)
    H, W = image.shape
    cy, cx = H // 2, W // 2
    radius = cutoff_fraction * min(H, W) / 2
    Y, X = np.ogrid[:H, :W]
    mask = (X - cx)**2 + (Y - cy)**2 <= radius**2  # circular mask
    F_filtered = F_s * mask
    return np.real(np.fft.ifft2(np.fft.ifftshift(F_filtered)))

# This gives heavy blur + ringing (Gibbs phenomenon) because
# the ideal low-pass has a sinc impulse response in the spatial domain.

The frequency-domain perspective unifies all the filters we have seen and opens up new designs.

Filter design vocabulary:

• Low-pass filter (LPF): passes low frequencies, attenuates high frequencies. Effect in space: blurring, smoothing.
• High-pass filter (HPF): passes high frequencies, attenuates low frequencies. Effect in space: edges, sharpening. HPF = 1 − LPF (their impulse responses sum to a delta function).
• Band-pass filter: passes a range of frequencies; used to isolate texture at a specific scale.
• Notch filter: removes specific periodic noise (e.g., a scanner interference pattern, visible as a spike in the Fourier spectrum).

The Gibbs phenomenon and filter design: An ideal brick-wall low-pass filter (hard cutoff in frequency domain) has a sinc function as its spatial-domain impulse response. Sincs have infinite tails, meaning the filter would require an infinitely large kernel. When truncated to a finite kernel, the abrupt cutoff produces ringing artifacts (oscillations near edges), called the Gibbs phenomenon. The Gaussian filter avoids this: it has no sharp cutoff, and its infinite spatial tails decay so rapidly (exponentially) that a small kernel captures almost all the energy.

Window functions (Hann, Hamming, Blackman) are used in 1D signal processing to taper the kernel and reduce ringing. In image processing, Gaussian kernels serve the same purpose naturally.

Practical Fourier filtering workflow:

1. Compute fft2(image) and fftshift.
2. Inspect the log-magnitude spectrum. Identify noise (sharp spikes at specific frequencies).
3. Construct a frequency-domain mask H[v,u].
4. Multiply: F_filtered = fftshift(fft2(image)) * H.
5. Reconstruct: result = real(ifft2(ifftshift(F_filtered))).

Phase matters, magnitude matters. A famous demonstration: swap the phase of one image with that of another while keeping each other's magnitudes. The reconstructed image looks like the phase donor, not the magnitude donor. Phase encodes structural information (where things are); magnitude encodes spectral power (how much of each frequency). For human perception, phase is dominant.

import numpy as np
from scipy.ndimage import gaussian_filter

def gaussian_low_pass(image, sigma):
    """Gaussian LPF — smooth in both space and frequency, no ringing."""
    return gaussian_filter(image, sigma=sigma)

def gaussian_high_pass(image, sigma):
    """HPF = original - LPF. Extracts fine details / edges."""
    return image - gaussian_low_pass(image, sigma)

def frequency_domain_notch(image, notch_freqs):
    """
    Remove periodic noise spikes at given (fy, fx) frequency pairs.
    notch_freqs: list of (fy, fx) in units of cycles/image,
                 where (0,0) is DC after fftshift.
    """
    F = np.fft.fftshift(np.fft.fft2(image))
    H, W = image.shape
    cy, cx = H // 2, W // 2
    mask = np.ones((H, W), dtype=complex)
    r = 5   # notch radius in pixels
    for (fy, fx) in notch_freqs:
        Y, X = np.ogrid[:H, :W]
        # Zero out a small disc around each noise frequency
        mask[(X-(cx+fx))**2 + (Y-(cy+fy))**2 <= r**2] = 0
        mask[(X-(cx-fx))**2 + (Y-(cy-fy))**2] = 0  # conjugate
    return np.real(np.fft.ifft2(np.fft.ifftshift(F * mask)))

# Example: separate a face image into structure and texture
lpf = gaussian_low_pass(gray_f, sigma=4)   # coarse structure
hpf = gaussian_high_pass(gray_f, sigma=4)  # fine details + edges

When we photograph a scene, we sample the continuous light field at the pixel grid. When we resize an image, we resample it. Sampling introduces a fundamental constraint: you can only faithfully represent frequencies up to half your sampling rate.

Nyquist–Shannon sampling theorem: To reconstruct a signal that contains frequencies up to f_max Hz, you must sample at a rate of at least 2·f_max samples per second. The Nyquist frequency is f_N = f_sample / 2. Any frequency above f_N, if present in the signal, will be aliased — it will masquerade as a lower frequency in the sampled signal, and you cannot tell the alias from genuine low-frequency content.

Visual aliasing examples:

• A fence photographed at an angle can produce moiré patterns where the fence spacing approaches the pixel spacing.
• When you zoom out (downsample) an image naively (dropping every other pixel), checkerboard textures create false large-scale patterns.
• A car wheel in a video can appear to spin backward when the frame rate and wheel rotation rate satisfy a particular relationship.

Anti-aliasing: The solution is to low-pass filter before sampling. When downsampling an image by a factor of k:

1. Apply a Gaussian blur with σ ≈ k/2 (or roughly, cutoff ≈ new_Nyquist).
2. Then subsample (take every k-th pixel).

Skipping step 1 causes aliasing. OpenCV's cv2.resize with INTER_AREA implements this correctly. Many naive implementations use bilinear interpolation without pre-blurring, which introduces subtle aliasing at high downsampling ratios.

In deep learning: convolutional neural network downsampling layers (stride-2 convolutions, max pooling) should apply anti-aliasing before subsampling. A 2019 paper by Richard Zhang showed that ignoring this makes CNNs unexpectedly sensitive to small translations — a failure mode that has direct roots in sampling theory.

import numpy as np
from scipy.ndimage import gaussian_filter
import cv2

def downsample_naive(image, factor):
    """Wrong: skips anti-aliasing, produces aliasing artifacts."""
    return image[::factor, ::factor]

def downsample_correct(image, factor):
    """
    Correct: low-pass filter first, then subsample.
    sigma = factor/2 is a rule of thumb for the anti-alias blur.
    """
    sigma = factor / 2.0
    blurred = gaussian_filter(image, sigma=sigma)
    return blurred[::factor, ::factor]

# Observe the difference on a checkerboard
H, W = 256, 256
y, x = np.mgrid[:H, :W]
checkerboard = ((x // 8 + y // 8) % 2).astype(np.float32)

aliased  = downsample_naive(checkerboard, factor=4)
correct  = downsample_correct(checkerboard, factor=4)
# aliased shows strong moiré; correct is approximately uniform gray

The Laplacian of Gaussian (LoG), also called the Mexican-hat filter, is one of the most beautiful operators in image processing: it is simultaneously a second-derivative edge detector and a blob detector, and it forms the mathematical foundation of scale-space theory.

The LoG combines Gaussian smoothing with the Laplacian:

LoG(x,y) = ∇²G(x,y) = (x² + y² − 2σ²) / σ⁴ · exp(−(x²+y²)/2σ²)

In practice, LoG(image) = Laplacian(Gaussian_blur(image, σ)) (applying them in sequence is equivalent by linearity of convolution).

As an edge detector: The LoG responds strongly at edges, and edge locations correspond to its zero crossings (sign changes). Zero-crossing detection gives subpixel edge localization — you can interpolate the exact position where LoG transitions from positive to negative.

As a blob detector: The LoG has a maximum (negative) response at a circular blob of radius r ≈ √2·σ. By computing the LoG at multiple scales σ and finding local maxima in the 3D (x, y, σ) scale-space, you detect blobs at their characteristic scale. This is the core of the SIFT (Scale-Invariant Feature Transform) detector, though SIFT approximates LoG with the Difference of Gaussians (DoG) for efficiency:

DoG(x,y,σ) = G(x,y,kσ) − G(x,y,σ)  ≈  (k−1)σ² · LoG(x,y,σ)

Subtracting two Gaussian-blurred images at nearby scales is much cheaper than computing the Laplacian, and it is a good approximation.

Gaussian scale space is the principled answer to the question: at what scale should I analyze this image? The idea (Lindeberg, 1994) is that a camera can photograph a scene from any distance, and a principled multiscale representation should not invent or destroy features as you change scale. The Gaussian is the unique kernel that guarantees this causality property — a feature at scale σ must have existed at finer scales.

This deep connection between differential geometry, partial differential equations, and image analysis is one of the most intellectually satisfying corners of computer vision.

import numpy as np
from scipy.ndimage import gaussian_filter, laplace

def log_filter(image, sigma):
    """Laplacian of Gaussian at scale sigma."""
    blurred = gaussian_filter(image, sigma=sigma)
    return laplace(blurred)   # scipy's discrete Laplacian

def dog_filter(image, sigma, k=1.6):
    """Difference of Gaussians: approximates (k-1)*sigma^2 * LoG."""
    return gaussian_filter(image, sigma=k*sigma) - gaussian_filter(image, sigma=sigma)

# Multi-scale blob detection: find local maxima in scale-space
sigmas = [1, 2, 4, 8, 16]
scale_space = np.stack([log_filter(gray_f, s) * s**2   # normalise for scale
                         for s in sigmas], axis=0)  # shape: (len(sigmas), H, W)
# Local maxima in this 3D volume are scale-space blobs.
# scipy.ndimage.maximum_filter on scale_space would find them.

The mathematical elegance of convolution meets real engineering constraints when you process images at video rate or on embedded hardware. Understanding the performance landscape makes you a better practitioner.

Complexity comparison: For an M×N image and a K×K kernel:

• Naive 2D convolution: O(M·N·K²) multiply-adds per image.
• Separable kernel (two 1D passes): O(M·N·2K).
• FFT-based convolution: O(M·N·log(M·N)), independent of kernel size.

The FFT method wins for large kernels (K ≳ 30), but has overhead from the FFT itself. For small kernels (K ≤ 15), spatial-domain convolution with SIMD vectorization is faster. For separable kernels at any size, the two-pass approach is almost always the practical choice.

Integral images (summed-area tables): For box filters specifically, you can compute the sum over any rectangular region in O(1) time using a precomputed integral image:

I_int[y,x] = Σ_{y'≤y, x'≤x} I[y',x']
box_sum(r1,c1,r2,c2) = I_int[r2,c2] - I_int[r1-1,c2] - I_int[r2,c1-1] + I_int[r1-1,c1-1]

This makes a box blur O(M·N) regardless of kernel size. OpenCV's integral() function builds this in linear time. The technique underpins the Viola-Jones face detector (2001), which evaluates Haar-like features in constant time.

GPU and SIMD: Modern image processing exploits data-level parallelism heavily. GPUs process thousands of pixels simultaneously (each thread handles one or a few output pixels). CPU SIMD instructions (SSE, AVX, NEON on ARM) process 8–32 values per instruction. OpenCV calls highly optimized SIMD paths automatically for standard operations.

Numerical precision matters: Operations on uint8 (0–255) values can overflow easily during convolution (e.g., summing 9 pixels of value 255 gives 2295, exceeding uint8). Always convert to float32 before applying non-trivial filters, then clamp back to [0,1] or [0,255] before saving.

import numpy as np

def integral_image(img):
    """Build 2D prefix-sum (summed area table) in O(H*W)."""
    return img.cumsum(axis=0).cumsum(axis=1)

def box_filter_integral(img, half_k):
    """Box filter using integral image — O(H*W) regardless of kernel size."""
    H, W = img.shape
    pad = half_k
    # Pad so boundary checks are trivial
    padded = np.pad(img.astype(np.float64), pad, mode='reflect')
    II = integral_image(padded)
    # Insert zero row and column for easy boundary arithmetic
    II = np.pad(II, ((1, 0), (1, 0)), mode='constant')
    # Output pixel (y,x) sums over [y:y+2k+1, x:x+2k+1] in padded
    k = half_k
    r1, c1 = np.arange(H), np.arange(W)
    r2, c2 = r1 + 2*k + 1, c1 + 2*k + 1   # broadcast over grid
    # Using np.ix_ for 2D broadcasting
    r1g, r2g = r1[:, None] + 1, (r1 + 2*k + 2)[:, None]
    c1g, c2g = c1[None, :] + 1, (c1 + 2*k + 2)[None, :]
    out = (II[r2g, c2g] - II[r1g-1, c2g]
           - II[r2g, c1g-1] + II[r1g-1, c1g-1])
    out /= (2*k+1)**2
    return out