Lab 9.4 — Image 2-D Convolution on the Embedded Target

← Course 3 syllabus · Bonus Module 9 · Prev: « Lab 9.3 · Next: Lab 9.5 »

Goal

Take the streaming discipline into two dimensions: the host sends a real grayscale image to the device, which runs 2-D convolution — a separable Gaussian blur and a Sobel edge detector — on it, block by block, and streams the result back for the host to save as a PNG and check against scipy.ndimage / OpenCV. On the memory-tight STM32 the key trick is streaming the image row-by-row through a small sliding window of rows, never holding the whole picture, which is exactly how line-buffer image pipelines work in real camera/vision firmware; larger images move to the Pi 5 / Jetson. The verification is quantitative: PSNR and max-abs-error vs. a reference implementation, so “the blur looks right” becomes a number. This is the image-signal-processing lab of the module and the bridge to the video pipeline in Lab 9.5.

Equipment & parts

  • STM32 Nucleo-L476RG for small images (see the SRAM budget below), or Raspberry Pi 5 / Jetson Orin Nano for larger images.
  • Host Mac with numpy scipy pillow (and optionally opencv-python).
  • A grayscale test image — a PGM (trivial to parse) or a PNG loaded via Pillow and converted to 8-bit gray. Start small (e.g. 128×128) for the STM32.

Wiring & bench setup

Same single-cable hookup as Lab 9.1 (Mac USB → ST-LINK Micro-B, USART2 VCP); nothing else on the bench. The blocks are now image rows:

flowchart LR
  PNG["media/in/test.png"]
  HOST["Host: Pillow → 8-bit gray<br/>Lab 9.1 framing, one row/block"]
  MCU["STM32 3-row window<br/>separable blur + Sobel"]
  OUT["media/out/blur.png<br/>media/out/edges.png"]
  REF["scipy.ndimage reference<br/>PSNR + max-abs-error"]
  PNG --> HOST
  HOST -- "header + rows" --> MCU
  MCU -- "output rows" --> HOST
  HOST --> OUT -- "compare" --> REF

flowchart LR
  PNG["media/in/test.png"]
  HOST["Host: Pillow → 8-bit gray<br/>Lab 9.1 framing, one row/block"]
  MCU["STM32 3-row window<br/>separable blur + Sobel"]
  OUT["media/out/blur.png<br/>media/out/edges.png"]
  REF["scipy.ndimage reference<br/>PSNR + max-abs-error"]
  PNG --> HOST
  HOST -- "header + rows" --> MCU
  MCU -- "output rows" --> HOST
  HOST --> OUT -- "compare" --> REF

Safety & don’t-break-it

Data lab; the constraints are memory and numeric scaling:

  • The L476RG has 128 KB SRAM — you cannot hold a big image. A 512×512 8-bit image is 256 KB, already twice the SRAM. This is why you stream row-by-row and keep only a 3-row sliding window (a few KB) plus one output row. Sizing the whole frame into RAM is the mistake that won’t even link.
  • Border handling must be decided, not accidental. Reading img[-1] or past the last row silently reads garbage. Choose a border policy (zero-pad, replicate/clamp, or reflect) and implement it explicitly at the top/bottom/left/right edges — and use the same policy in the host reference or the PSNR will be wrong at the borders.
  • Fixed-point scaling overflows fast in 2-D. A 3×3 kernel sums nine products; an unnormalized Sobel or an integer Gaussian can exceed 8 bits. Accumulate in int16/int32, normalize (divide by the kernel sum for the blur), and saturate back to uint8 — never let the byte wrap.
  • Keep the Lab 9.1 flow control: a row is a block; don’t let convolution compute stall the RX DMA.

Project & environment setup

Firmware — reuse firmware/m9-media/ (Lab 9.1). No new peripherals and no CMSIS-DSP needed (the 3×3 passes are integer C); confirm USART2 921600 + DMA + 80 MHz clock (setup essentials) as before.

Host — course venv plus the image libraries this lab adds:

source venv/bin/activate
pip install pillow opencv-python   # Pillow: PNG → 8-bit gray (PIL.Image); OpenCV optional cross-check
mkdir -p labs/lab-9-4/host labs/lab-9-4/captures

Library roles: PIL.Image (load/convert/save), scipy.ndimage.convolve1d (matched separable reference, as in Part C; scipy.signal.convolve2d is the non-separable cross-check), numpy (PSNR/max-abs-error math), cv2.GaussianBlur/cv2.Sobel (optional second opinion — remember they use their own kernels/normalization, so score against the matched convolve1d reference, not these). Transport is labs/lab-9-1/host/harness.py; the row-streamer + scoring script go in labs/lab-9-4/host/ (you write them). Input media/in/test.png — the snippet’s in.png.

Keep this lab’s reconciliation in labs/lab-9-4/host/analysis.ipynb — the notebook convention — and export final figures next to it.

Where results go:

Artifact Path
Bench note (PSNR / max-abs-error, SRAM budget) labs/lab-9-4/notes.md
Device blur / edge outputs media/out/blur.png · media/out/edges.png
Reference images labs/lab-9-4/host/ref-blur.png · labs/lab-9-4/host/ref-edges.png
Score logs (PSNR, max-abs-error, per border mode) labs/lab-9-4/captures/scores.txt

Background

2-D convolution. For an image \(I\) and kernel \(K\) of size \((2a{+}1)\times(2b{+}1)\):

\[(I * K)[i,j] \;=\; \sum_{u=-a}^{a}\sum_{v=-b}^{b} K[u,v]\,I[i-u,\,j-v].\]

A direct \(M\times N\)-pixel convolution with a \(k\times k\) kernel costs \(O(MNk^2)\) multiply-adds.

Separability. A Gaussian kernel factors as an outer product of two 1-D Gaussians,

\[K = g\,g^{\mathsf T}, \qquad g[t] = \tfrac{1}{Z}\,e^{-t^2/(2\sigma^2)},\]

so the 2-D convolution becomes two 1-D passes — convolve each row with \(g\), then each column with \(g\):

\[I * K \;=\; \big(I *_{\text{rows}} g\big) *_{\text{cols}} g.\]

This drops the cost from \(k^2\) multiply-adds per pixel to \(2k\) — for the \(k=3\) kernel here, \(9 \to 6\) (a \(1.5\times\), i.e. 33%, saving; the asymptotic factor is \(k/2\)) — and, in Course 1 Part I terms, is just the rank-1 factorization of the kernel matrix. The 3-row sliding window on the STM32 is exactly what the column pass needs: a \(3\times3\) kernel touches only rows \(i{-}1,i,i{+}1\), so you buffer three input rows, emit one output row, slide.

Sobel edges. The Sobel operator estimates the image gradient with two separable \(3\times3\) kernels:

\[G_x = \begin{bmatrix} 1 & 0 & -1\\ 2 & 0 & -2\\ 1 & 0 & -1\end{bmatrix} = \begin{bmatrix}1\\2\\1\end{bmatrix}\begin{bmatrix}1&0&-1\end{bmatrix}, \qquad G_y = G_x^{\mathsf T},\]

and the edge (gradient-magnitude) image is

\[G[i,j] = \sqrt{G_x[i,j]^2 + G_y[i,j]^2}.\]

Both \(G_x\) and \(G_y\) are separable (a smoothing \([1,2,1]\) and a difference \([1,0,-1]\)), so the same line-buffer machinery applies.

PSNR. To compare the device output \(\hat{I}\) against the reference \(I_{\text{ref}}\) (8-bit, peak \(=255\)):

\[\text{MSE} = \frac{1}{MN}\sum_{i,j}\big(\hat{I}[i,j]-I_{\text{ref}}[i,j]\big)^2, \qquad \text{PSNR} = 10\log_{10}\frac{255^2}{\text{MSE}}\ \text{dB}.\]

Identical integer math gives \(\text{MSE}=0\) (PSNR \(=\infty\)); a difference of a few LSB from rounding/border choices gives a high but finite PSNR you can predict.

Theory — Images: Convolution, Filtering, and Restoration

This section is the image-processing theory formerly carried as a lesson of the Course 1 booklet, moved here because this lab (with Lab 9.5 and Lab 9.6) is where it is used; the edge, corner, and motion half lives in Lab 9.6’s theory section.

Images as arrays

A continuous scene \(f(s,t)\) becomes a digital image by sampling and quantization — Course 1 Lesson 29’s story, run once per axis. The result is an \(M \times N\) array \(f(x,y)\) with the origin at the top left and intensities on \(L=2^k\) integer levels. Aliasing is also Lesson 29 twice over: a band-limited image is recoverable iff the sampling rate exceeds twice the bandwidth per axis, jaggies and moiré are the 2-D aliases, and the rules are the same — anti-alias before sampling (optical blur; there is no software cure afterward), and resampling is sampling again, so smooth before shrinking.

Spatial filtering: correlation and convolution

A linear spatial filter replaces each pixel by a sum of products of an \(m \times n\) kernel \(w(s,t)\) (\(m = 2a+1\), \(n = 2b+1\), odd) with the neighborhood:

\[ g(x,y) = \sum_{s=-a}^{a} \sum_{t=-b}^{b} w(s,t)\, f(x+s,\, y+t). \]

This is correlation; convolution pre-rotates the kernel by \(180^\circ\). The litmus test: convolving with a unit impulse yields an exact copy of the kernel (Course 1 Lesson 25’s sifting property, twice over); correlating yields it rotated. The two coincide only for symmetric kernels — which is why the distinction is invisible until it bites. Borders need padding (\(a\) rows, \(b\) columns minimum; zero, replicate, and mirror padding trade artifacts), and convolution’s algebra — commutative, associative, distributive — is what the transform theory below requires.

Smoothing and sharpening kernels

Box kernels (all ones, normalized) are cheapest and worst: not isotropic (their transform is a 2-D sinc with directional lobes). Gaussian kernels \(G(s,t) = K e^{-(s^2+t^2)/2\sigma^2}\) are the only circularly symmetric separable kernels; size \(\lceil 6\sigma \rceil\) (odd) captures essentially all the mass — bigger is waste, smaller truncates. Cascaded Gaussians obey \(\sigma^2 = \sigma_1^2 + \sigma_2^2\), and repeated box filtering tends Gaussian — the central limit theorem, now with kernels. Sharpening is differentiation: the digital Laplacian

\[ \nabla^2 f = f(x{+}1,y) + f(x{-}1,y) + f(x,y{+}1) + f(x,y{-}1) - 4f(x,y), \]

kernel \([0\ 1\ 0;\ 1\ {-4}\ 1;\ 0\ 1\ 0]\), gives \(g = f - \nabla^2 f\); unsharp masking (\(g = f + k(f - \bar f)\), subtract a blur, add the difference back) is the same idea in disguise.

Separability

A rank-1 kernel \(w = \mathbf{v}\mathbf{w}^\mathsf{T}\) filters in two 1-D passes: \((m+n)\) multiplies per pixel instead of \(mn\) — for a \(15 \times 15\) kernel, a \(7.5\times\) saving. Box, Gaussian, and Sobel kernels all qualify, and the two-pass result equals the one-pass result to machine precision. On an embedded target this is the difference between a frame rate and a slide show (this lab).

The 2-D DFT and frequency-domain filtering

\[ F(u,v) = \sum_{x=0}^{M-1} \sum_{y=0}^{N-1} f(x,y)\, e^{-j2\pi(ux/M + vy/N)} \]

— Course 1 Lesson 11’s DFT run on rows then columns (the kernel is separable), so the FFT gives \(O(MN \log MN)\). Translation changes phase only; real images give conjugate-symmetric transforms; \(F(0,0) = MN \bar f\); magnitude says how much of each frequency, phase says where. Products of DFTs implement circular convolution, so linear filtering demands padding (typically to \(2M \times 2N\)): unpadded, a blur wraps around and borrows from the opposite edge. Filtering is \(g = \mathrm{IDFT}[H F]\) with a real, centered, symmetric \(H\); with \(D(u,v)\) the distance from the center of the (padded) frequency rectangle, the ideal lowpass (\(H = 1\) for \(D \le D_0\), else \(0\)) rings — its spatial kernel is a 2-D sinc, so this is Gibbs — while the Gaussian \(H = e^{-D^2/2D_0^2}\) is ring-free, and highpass filters come free as \(1 - H_{LP}\).

Restoration: why the inverse filter explodes, and Wiener

Restoration models the degradation: \(g = h*f + \eta\), i.e. \(G = HF + N\). Direct inversion \(\hat F = G/H = F + N/H\) blows up wherever \(H\) is small — and real blurs guarantee such places: uniform linear motion gives a sinc-shaped \(H\) with periodic zeros, and the streaky “curtain of noise” in inverse-filtered images is exactly \(N/H\) erupting there. The principled fix is the Wiener filter of Lab 6.8’s theory section with \(D(u,v)\) in place of \(\omega\): minimizing \(\mathbb{E}[(f-\hat f)^2]\) gives

\[ \hat F = \left[ \frac{1}{H}\, \frac{|H|^2}{|H|^2 + S_\eta/S_f} \right] G \;\approx\; \left[ \frac{1}{H}\, \frac{|H|^2}{|H|^2 + K} \right] G, \]

inverse filtering tempered by the noise-to-signal ratio: \(\to 1/H\) where signal dominates, \(\to 0\) where noise does, with the constant-\(K\) form tuned interactively when the spectra are unknown. (Tikhonov/constrained least squares replaces \(S_\eta/S_f\) by \(\gamma|P(u,v)|^2\) with \(P\) the Laplacian’s transform.) For impulse (salt-and-pepper) noise no linear filter wins: the median of the window discards impulses outright — probability’s tool, not Fourier’s.

NoteConnection

Everything here is Course 1 Part V promoted to two indices: the filter sum is Lesson 25’s convolution with two subscripts, the 2-D DFT is Lesson 11’s unitary transform on rows then columns, and Wiener deconvolution is the orthogonality principle of Lesson 19 in the frequency plane. This lab implements exactly this section’s 2-D convolution on the embedded target (separability and fixed-point kernels decide whether the frame loop closes in time); Lab 9.5 runs per-frame Sobel / thresholding; Lab 9.6’s motion estimation is the block-matching paragraph of its theory section; and Lab 8.3’s learned denoiser is benchmarked against the Wiener baseline above. For computer vision this is the classical prerequisite layer its geometry consumes.

Worked by hand

P1. Apply the \(-4\)-center Laplacian kernel to the center of the patch \([10\ 10\ 10;\ 10\ 20\ 10;\ 10\ 10\ 10]\) and interpret the sign of the response.

P3. At a frequency where \(|H| = 0.1\), compare the inverse filter’s gain with the Wiener gain for \(K = 0.01\); then show the constant-\(K\) Wiener formula reduces to \(1/H\) as \(K \to 0\).

(P2, the Sobel gradient computation, is in Lab 9.6.)

Theory exercises

Worked by hand, like every exercise on this site; statements only.

Theory exercise 1 [Hand] — (G&W 3.27) An image is filtered four times with a \(3 \times 3\) Gaussian kernel of standard deviation 1.0. By associativity, the same result comes from a single Gaussian kernel. (a) What is its size? (b) What is its standard deviation?

Theory exercise 2 [Hand] — (G&W 4.47) A \(3 \times 3\) kernel averages the four closest neighbors of a point \((x,y)\) but excludes the point itself from the average. (a) Find the equivalent frequency-domain transfer function \(H(u,v)\). (b) Show that it is a lowpass filter transfer function.

Theory exercise 3 [Hand] — (G&W 5.30) A motion-blurred image corrupted by additive Gaussian noise is inverse-filtered. The blurring itself is corrected, but the restored image shows a strong streak pattern that is not apparent in the blurred image. Explain how this pattern originated.

Procedure

Part A — Host: send the image row-major.

  1. Load and convert to 8-bit gray, send dimensions in a header frame, then stream one row per block through the Lab 9.1 harness:
from PIL import Image
import numpy as np
img = np.asarray(Image.open("in.png").convert("L"), dtype=np.uint8)   # H×W
send_header(img.shape)                       # H, W, kernel id, border policy
for row in img:                              # one block per row
    send_block(row.tobytes())

Part B — Firmware: separable convolution over a 3-row window.

  1. Keep a ring of three input rows. On each new row, run the row (horizontal) 1-D pass into the window, then the column (vertical) 1-D pass across the three buffered rows to produce one output row. Illustrative kernel (structure only):
/* rows[3][W] hold i-1, i, i+1 after the horizontal pass.  Gaussian [1 2 1]/4 sep. */
for (int j = 0; j < W; ++j) {
    int32_t acc = rows[0][j] + 2*rows[1][j] + rows[2][j];   // vertical pass
    out_row[j] = (uint8_t)(acc >> 2);                        // normalize /4, saturate
}
/* Sobel: Gx = sep([1 2 1]^T, [1 0 -1]);  Gy = transpose. mag = sqrt(gx*gx+gy*gy) */
  1. Handle top/bottom rows with the chosen border policy (replicate is simplest). Stream each finished output row back.

Part C — Host: reassemble and score.

  1. Reassemble output rows into an array, save out.png, and score against the reference with the same border mode:
from scipy.ndimage import convolve1d
# Reference MUST be the SAME integer kernel the device runs — a [1,2,1]/4 binomial
# blur is NOT a true Gaussian of any sigma, so scoring against gaussian_filter(sigma)
# would guarantee a nonzero "error" that is really just a different kernel.
k = np.array([1.0, 2.0, 1.0]) / 4.0
tmp = convolve1d(img.astype(np.float64), k, axis=1, mode="nearest")   # horizontal pass
ref = convolve1d(tmp,                    k, axis=0, mode="nearest")   # vertical pass
mse = np.mean((out.astype(np.float64) - ref)**2)
psnr = 10*np.log10(255.0**2 / mse) if mse > 0 else np.inf
maxerr = np.max(np.abs(out.astype(int) - np.rint(ref).astype(int)))

Deliverable & expected results

  • blur.png and edges.png from the device, the reference images, and the PSNR / max-abs-error for each.

For a 3×3 separable Gaussian (\([1,2,1]/4\) each pass) and Sobel on an 8-bit image, integer device math vs. a matched integer reference:

Quantity Predicted Measured
Separable vs. direct multiply count (3×3) \(2k\) vs \(k^2\) per pixel = 6 vs 9 → 33% fewer …
Blur PSNR vs. reference (matched rounding) \(\gtrsim 45\) dB (few-LSB rounding) or \(\infty\) if bit-exact …
Blur max abs error \(\le 1\)–2 LSB …
Sobel edge-map max abs error \(0\) (integer, matched border) …
STM32 SRAM for 3-row window (W=512) \(\approx 3\times512\) B + out row ≈ 2 KB …
Full 512×512 image in SRAM 256 KB > 128 KB → must stream …

Analysis & reconciliation

If the interior of the blur matches the reference but the borders are off, your device border policy and the reference mode= disagree — align them and the PSNR jumps. A uniform small error across the whole image is rounding: integer >>2 truncates where gaussian_filter rounds — switch to round-to-nearest ((acc + 2) >> 2) and watch the max-abs-error drop to 0. A structured error (e.g. every column shifted) is an indexing/transpose bug in the separable column pass, not quantization. Reconcile the measured multiply count against the \(2k\)-vs-\(k^2\) (6-vs-9) separable prediction — that speedup is the whole reason to factor the kernel (Course 1 Part I’s rank-1 view). For Sobel, integer math with a matched border should be exactly the reference, so any nonzero max-abs-error there is a bug to find, not a bound to accept.

Going further

  • Add reflect and zero border modes and quantify their effect on edge PSNR near the frame boundary.
  • Do the Gaussian at Q15 fixed point and compare its PSNR to the integer version — the 2-D analog of Lab 9.2’s fixed-point study.
  • Move a 1024×1024 image to the Pi 5 or Jetson (over the Lab 9.1 Jetson transport) and hold the whole frame; compare wall-clock to the STM32’s streamed version and note where the memory-vs-throughput tradeoff flips. On the Jetson, add the GPU rung — cupyx.scipy.ndimage.gaussian_filter/sobel on the full frame — and carry the three-way number (STM32 streamed / SBC CPU / Jetson GPU) into Lab 9.5, which runs this per-frame at video rate.
  • Chain blur → Sobel (edge detection on a denoised image) entirely in the streaming pipeline — the exact per-frame operation you’ll run on video in Lab 9.5.