Image Filtering: 20-Bio-B4 Robotics, December 2014, Q2

Question 2: Image Filtering (20 marks)

Question text not reproduced: the examination questions are © Engineers and Geoscientists BC. Open the official past paper (linked at the top of this page) to read the question, then follow the worked solution below.

Given. An image $f$ of size $M_1 \times N_1$ convolved with a kernel $h$ of size $M_2 \times N_2$.

Find. The three standard output sizes, a direct-convolution algorithm, an FFT-based algorithm, their computational complexities, and the invertibility condition on $h$.

(a) The Three Convolution Output Sizes

Two-dimensional discrete convolution slides the (flipped) kernel over the image; how far the kernel is allowed to overhang the image border fixes the output size:

The three standard 2-D convolution output sizes for an $M_1 \times N_1$ image and an $M_2 \times N_2$ kernel.

"Full" convolution keeps every non-zero overlap position, including where the kernel only partially overlaps $f$: size $(M_1+M_2-1) \times (N_1+N_2-1)$.
"Same" convolution crops the full result (or equivalently zero-pads $f$ before convolving) so the output matches the input size: $M_1 \times N_1$. This is the size used by most filtering functions by default.
"Valid" convolution keeps only positions where the kernel lies entirely inside $f$ (no zero-padding assumptions at the border): size $(M_1-M_2+1) \times (N_1-N_2+1)$.

(b) Direct (Spatial-Domain) Pseudo-Code

Direct convolution flips the kernel and, at every output position, sums the elementwise products of the kernel with the underlying image patch (here shown producing the "same"-size output, with zero-padding at the border):

function g = convolve2D(f, h):
    // f is M1 x N1, h is M2 x N2 (M2, N2 assumed odd)
    a = floor(M2 / 2); b = floor(N2 / 2)
    g = zeros(M1, N1)
    for m = 0 .. M1-1:
        for n = 0 .. N1-1:
            acc = 0
            for i = -a .. a:
                for j = -b .. b:
                    if (m-i, n-j) is inside f:
                        // h is flipped: index (a-i, b-j) is h(-i,-j)
                        acc += f(m-i, n-j) * h(a-i, b-j)
                    // else: f is treated as 0 (zero-padding)
            g(m, n) = acc
    return g

(c) Fast Convolution via the Fourier Transform

The convolution theorem states that convolution in the spatial domain is elementwise multiplication in the frequency domain. To avoid the wraparound error that a naive equal-size DFT would introduce (circular convolution ≠ linear convolution), both $f$ and $h$ are first zero-padded to at least $(M_1+M_2-1) \times (N_1+N_2-1)$:

Zero-pad $f$ and $h$ to a common size $P \times Q \ge (M_1+M_2-1) \times (N_1+N_2-1)$ (rounded up to a size efficient for the FFT, e.g. a power of two).
$F = \text{FFT}(f_{\text{padded}}), \quad H = \text{FFT}(h_{\text{padded}})$
Multiply elementwise in the frequency domain: $G(u,v) = F(u,v)\,H(u,v)$.
Inverse transform and crop: $g = \text{IFFT}(G)$, then crop to the desired "full"/"same"/"valid" size from (a).

This, the same sequence convolved via zero-padded forward DFT → pointwise multiply → inverse DFT — confirming the convolution theorem invoked above.

(d) Computational Complexity and the Crossover

Approach. Count multiplications for each method as a function of the image and kernel sizes, then compare on a concrete example.

Direct method. Each of the $M_1 N_1$ output pixels performs $M_2 N_2$ multiplications: $$C_{\text{direct}} = M_1 N_1 M_2 N_2$$
FFT-based method. Padding to $P \times Q$ (with $P,Q = O(M_1+M_2)$), each 2-D FFT/IFFT costs $O(PQ\log_2(PQ))$; there are three transforms (forward $f$, forward $h$, one inverse) plus $PQ$ pointwise multiplications: $$C_{\text{FFT}} \approx 3\,PQ\log_2(PQ) + PQ$$
Crossover. For a $256\times256$ image with a small $5\times5$ kernel, $C_{\text{direct}} = 256^2\cdot25 = 1{,}638{,}400$, against $C_{\text{FFT}} \approx 1.44\times10^7$ (both $f$ and $h$ must be padded to at least $260\times260$ to avoid wraparound, which rounds up to a $512\times512$ FFT) — the direct method wins by roughly a factor of 9. For a large kernel approaching the image size (e.g. $129\times129$ on the same image), $C_{\text{direct}}$ grows to $\approx1.09\times10^9$, while the required padding is still $256+129-1=384$, rounding up to the same $512\times512$ transform, so $C_{\text{FFT}}$ is unchanged at $\approx1.44\times10^7$ — now the FFT method wins by roughly a factor of 75. The comparison shows the FFT-based cost is set almost entirely by which power-of-two padded transform size the kernel forces, not by the kernel area directly, whereas the direct cost scales linearly with $M_2N_2$ regardless of padding.

Case	Direct multiplications	FFT-based multiplications (approx.)	Winner
$256\times256$ image, $5\times5$ kernel	1,638,400	≈1.44×10⁷	Direct
$256\times256$ image, $129\times129$ kernel	≈1.09×10⁹	≈1.44×10⁷	FFT

Answer to (d): when $M_2, N_2 \ll M_1, N_1$, the direct method is more efficient, because its cost scales with the small kernel area $M_2N_2$, whereas the FFT approach pays a near-fixed overhead set by the (much larger) padded image size regardless of how small the kernel is.

(e) Conditions on $h$ for Deconvolution

Recovering $f$ from $g = f * h$ means, in the frequency domain, dividing $G(u,v)$ by $H(u,v) = \text{FFT}(h)$. This is only well-posed if:

$H(u,v) \neq 0$ for every $(u,v)$ in the band of interest — any frequency where $H$ has a zero (or is very small) makes the inverse filter $1/H(u,v)$ blow up, and any noise at that frequency is amplified without bound.
$h$ (equivalently $H$) must be known (non-blind deconvolution) or reliably estimable (blind deconvolution, which additionally requires a statistical/parametric model of $h$ and/or $f$ to make the problem determined).
In the presence of noise, exact inversion is numerically unstable even where $H \neq 0$ but is merely small; a regularized inverse (e.g. the Wiener filter, $\hat{F} = \dfrac{H^{*}G}{|H|^2 + K}$ with $K$ set by the noise-to-signal power ratio) is required in practice rather than the naive $G/H$.

Concept & Theory: Convolution, the Convolution Theorem, and Deconvolution

Topic: 2-D discrete convolution, fast convolution via the DFT/FFT, and the ill-posedness of inverse filtering.
Key relations: $g = f * h \Leftrightarrow G = F \cdot H$ (convolution theorem); direct cost $M_1N_1M_2N_2$ vs. FFT cost $O(PQ\log PQ)$; Wiener deconvolution $\hat{F} = H^*G/(|H|^2+K)$.
Why this works: the DFT diagonalizes the (circulant) convolution operator, turning an $O(n^2)$ spatial operation into an $O(n\log n)$ elementwise one, at the cost of needing zero-padding to convert the DFT's implicit circular convolution back into linear convolution.
Common pitfall: forgetting to zero-pad before the FFT (producing wraparound/circular-convolution artifacts at the image border), and naive division by a near-zero $H(u,v)$ during deconvolution, which catastrophically amplifies noise.
Source: Gonzalez & Woods, Digital Image Processing, ch. 4 (frequency-domain filtering, the convolution theorem) and ch. 5 (image restoration and inverse/Wiener filtering).

Question 2 of 6: Image Filtering