Multi-class segmentation of satellite imagery with 72 labeled samples and no deep network: a classical restoration chain, geodesic Voronoi superpixels, per-region statistics and an XGBoost classifier over 6 cartography classes. A segmentation network unbundled into explicit, sample-efficient pieces.
Segmentation has become a reflex: see the problem, grab a U-Net, find labels, train. This toy project asked whether the reflex is necessary. With 72 annotated samples, 6 classes, and no path to the thousands of labels a network wants, I unbundled what a segmentation network does implicitly (denoise, propose regions, classify regions) and implemented each stage explicitly. Prototyping happened in Julia; the final pipeline was rewritten in C++.
The whole pipeline:
(All other figures are illustrative reconstructions on a synthetic toy scene, not screenshots of the actual project data.)
01The data, and why it was worse than it sounds
Satellite imagery is usually multispectral, a dozen+ bands including near-infrared where water and vegetation nearly classify themselves (hence NDVI). I had none of that: plain RGB, plus heavy compression artifacts (block boundaries, ringing, quantization banding).
Compression artifacts are structured noise: fake edges on an 8x8 grid that destroy texture distinguishing, say, bare soil from a construction site. Any region-growing or edge-aware algorithm will latch onto block boundaries if allowed to, which is why the cleaning stack here is unusually heavy: damage control, not perfectionism.
The 72 samples came in varying dimensions, fragments of a smaller number of scenes, so the first step reconstructed the 8 original satellite tiles. From those, adaptive tiling generated ~12k patches of 128x128, grid chosen per tile so patches land at native resolution and avoid bilinear downsampling (a low-pass filter that would smear what compression already ate). The fixed 128x128 size also caps per-patch compute.
02Cleaning: four algorithms, four different problems
No single filter fixes compressed RGB. Four stages, each aimed at one degradation, order not negotiable.
1. Non-Local Bayes + SOS boosting, against stochastic noise. NL-Bayes (Lebrun, Buades and Morel, 2013) groups similar patches and denoises each group with a Bayesian estimate under a Gaussian patch model; satellite scenes are extremely self-similar. SOS boosting (Romano and Elad, 2015) feeds the denoiser its own previous estimate added back to the input, then subtracts it, recovering detail the first pass over-smoothed. Denoising comes first because deconvolution amplifies whatever noise it receives.
2. Short-term wavelet shrinkage, against residual grain. Patch models leave a floor of low-amplitude residue. A light soft-threshold on the finest wavelet scales (Donoho and Johnstone, 1994) removes it without touching structural scales.
3. Richardson-Lucy deconvolution, against optical and atmospheric blur. The imaging chain convolves ground truth with a blur kernel. Richardson (1972) and Lucy (1974) give the classic iterative maximum-likelihood deblur under Poisson statistics, matching sensor photon-counting. RL amplifies noise (stages 1-2 cover that) and rings at borders (one more reason to deconvolve reconstructed tiles, not fragments); iterations are capped early.
4. CLAHE, for the segmenter, not the eye. Contrast Limited Adaptive Histogram Equalization (Zuiderveld, 1994) equalizes local windows with a clip limit. This is conditioning, not restoration: it makes weak-but-real boundaries (dirt road against sand, shallow water edge) produce gradients strong enough to stop a geodesic front. It runs last so its nonlinear remapping doesn't violate the earlier stages' statistical assumptions.
03Superpixels: a geodesic Voronoi diagram with a shortcut
Superpixels replace the learned encoder, turning 16k pixels into a few hundred coherent regions.
Following Zeng et al., Structure-Sensitive Superpixels via Geodesic Distance (ICCV 2011): each pixel belongs to the seed with smallest geodesic distance, path cost accumulating local color variation. Crossing an edge is expensive, so region borders snap onto image structure. A plain Euclidean Voronoi diagram, or spatially-dominated clustering like SLIC, cuts straight through a shoreline when the spatial term wins; the geodesic one doesn't.
The original paper's expensive part is the loop (place seeds, partition, relocate, repeat), mostly recovering from a bad start. So I initialized better instead of iterating more:
- Quantize the cleaned patch to 16 colors, dynamically, by coverage of the color space rather than pixel count: a mass-weighted quantizer would burn ten shades of beige on a 70%-sand scene and starve roads and rooftops. This way every distinct material gets a bin regardless of area.
- Flood-fill the quantized map into connected regions of near-constant color.
- Lay a grid of seeds inside each region, density proportional to area.
Seeds are already color-aware before the first geodesic pass: none straddles a strong edge, water gets sparse coverage, industrial areas get dense coverage. The partition needs dramatically fewer iterations than the paper's uniform-grid init, close to one pass in practice, and boundary adherence is better since iterative relocation can drift centroids next to edges.
04From regions to labels
Each superpixel is summarized by five statistics per RGB channel: min, max, mean, stddev, skewness. Fifteen features, deliberately spartan: with 6 classes and 72 annotations, feature count is a regularization decision, and these moments already separate classes well: water is dark, low-variance, negatively skewed; sand is bright and flat; built-up areas have high variance and heavy tails from roof/shadow mixtures.
The classifier is XGBoost (Chen and Guestrin, 2016), argmax over 6 classes per superpixel. Gradient-boosted trees fit this well: strong on small tabular data, indifferent to feature scaling, axis-aligned splits matching the thresholded nature of the problem.
05The hardest part: turning predictions back into a map
Everything up to here is the part people expect to be hard, and wasn't. The genuinely hard part starts after the classifier speaks.
The classifier gives one prediction per superpixel by design: max() over the segmentation, argmax of class scores, done. Regions are internally consistent by definition, so reconstruction becomes "assemble a few hundred labeled puzzle pieces per patch" instead of reconciling 16k independent pixel opinions. But assembling those pieces splits into three sub-problems:
Fusing superpixels back to tile scale. Predictions are per 128x128 patch; the deliverable is a labeled tile. Adjacent masks must be merged, same-label neighbors fused into single units, geometry realigned by inverting the tiling transform. None of it is conceptually deep, all of it is where bugs concentrate: off-by-one alignments at patch borders look exactly like classifier errors and aren't.
Smoothing the label field. Mispredictions are sparse, isolated superpixels, not large wrong areas, because errors are close to spatially uncorrelated while true land cover is heavily autocorrelated. A robust spatial smoothing pass (majority reassignment of isolated regions from their neighborhood) exploits that mismatch, converting many isolated mistakes into correct labels essentially for free: the cheapest accuracy gain in the pipeline.
The road not taken: keeping the full distribution. Argmax throws information away. The alternative keeps the entire class-probability vector per superpixel and predicts all of them jointly, letting neighboring distributions negotiate (a 55/45 water-vs-sand region surrounded by confident water should tip over). That's structured, non-linear inference over the region adjacency graph: CRFs are the classical answer, but exploiting the non-linearity really wants a learned model, e.g. a graph network over the adjacency structure. A funny place for an AI-free pipeline to land: the classical stack carries cleaning, segmentation and classification, and the first place a network earns its keep is the last mile of spatial reasoning. For this project, argmax plus robust smoothing was the right trade; for production, that last mile is where I'd spend the learning budget.
06Ideas to push it further
- Scattering features instead of moments. The fifteen moments are blind to texture, the known ceiling of the current features. The invariant scattering transform (Mallat, 2012; Bruna and Mallat, 2013) gives translation-invariant, deformation-stable texture descriptors with no training data. The single change I'd make first; with more data, a SIFT-CNN hybrid on scattering is the natural next rung.
- A principled geodesic weight. The color-delta edge weight is a hand-tuned proxy for "structure at the scale I care about." Deriving weights from scale-space theory (Lindeberg, 1994) or renormalization-group coarse-graining would make them a function of actual scale rather than an eyeballed knob, and likely speed up the diagram computation too.
- Retinex as a soft cleaner, early in the chain. Retinex (Land and McCann, 1971) and its multi-scale variant (Jobson et al., 1997) separate illumination from reflectance, precisely the satellite problem: sun angle and haze modulate the same ground truth differently across tiles, leaking into the color moments as inter-tile drift.
- Marchenko-Pastur as a spectral denoiser. Random matrix theory says pure-noise covariance eigenvalues fill a known Marchenko-Pastur bulk; anything escaping is signal. Shrinking in-bulk eigenvalues of patch covariances, cutoff via Gavish and Donoho's optimal threshold, could replace the first two cleaning stages with something nearly parameter-free.
- Conditional Random Field for smoothing+reconstitution. The current pain point is argmax compression and its translation back to images, which only works well with a low-error forecaster. A CRF jointly reasons about the entire label field instead of predicting each superpixel independently, letting neighboring distributions influence each other without leaking noise. CRF example on satellite images.
07Closing thoughts
Can we do multi-class segmentation without a DNN? On this problem, with 72 labels and compressed RGB: yes, comfortably. Every piece of structure injected without learning it, restoration physics, geometric priors, robust regional statistics, is a few hundred labels not needed. The pipeline bets that forty years of classical vision literature is worth a few thousand annotations, and the bet paid.
It also clarified where deep learning belongs: not the encoder, which the classical stack covers well at small scale, but the structured, non-linear spatial reasoning after prediction, where for a toy project a robust smoothing filter got close enough. Julia made experimentation fast, C++ made the final version fast, neither asked for a GPU.
08Further reading
- Zeng et al., Structure-Sensitive Superpixels via Geodesic Distance, ICCV 2011, the geodesic Voronoi formulation this pipeline adapts.
- Lebrun, Buades and Morel, A Nonlocal Bayesian Image Denoising Algorithm, 2013, NL-Bayes, with a runnable online demo at IPOL.
- Romano and Elad, Boosting of Image Denoising Algorithms, 2015, SOS boosting, works on top of any denoiser.
- Achanta et al., SLIC Superpixels Compared to State-of-the-Art, PAMI 2012, the standard superpixel baseline and its evaluation metrics.
- Chen and Guestrin, XGBoost: A Scalable Tree Boosting System, KDD 2016.
- Bruna and Mallat, Invariant Scattering Convolution Networks, PAMI 2013, training-free texture features, the most promising upgrade.
- Gavish and Donoho, The Optimal Hard Threshold for Singular Values is 4/sqrt(3), 2014, the practical entry point to Marchenko-Pastur denoising.