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 satellite samples, 6 cartography classes and no realistic path to the thousands of labels a network wants, fighting the data budget was pointless, so I unbundled what a segmentation network does implicitly (denoise, propose regions, classify regions) and implemented each stage explicitly. Prototyping happened in Julia, where iterating on image pipelines is pleasant; the final pipeline was rewritten in C++.
The whole pipeline looks like this:
(All other figures in this post are illustrative reconstructions on a synthetic toy scene, made to show each mechanism clearly, not screenshots of the actual project data.)
01The data, and why it was worse than it sounds
Satellite imagery is usually multispectral or hyperspectral: a dozen or more bands, including near-infrared channels where water and vegetation practically classify themselves (this is why indices like NDVI exist). I had none of that. My source was plain RGB, and on top of losing the spectral bands, the images carried heavy compression artifacts: block boundaries, ringing around edges, quantization banding in smooth areas.
That single fact shaped the whole design. Compression artifacts are structured noise: they create fake edges on an 8x8 grid and destroy the subtle texture that distinguishes, say, bare soil from a construction site. Any region-growing or edge-aware algorithm downstream will happily latch onto block boundaries if you let it. This is why the cleaning stack in this pipeline is unusually heavy for a segmentation project. It is not perfectionism, it is damage control.
The 72 annotated samples came in varying dimensions, which is awkward on its own. Looking at them, it was clear they were fragments of a smaller number of scenes, so the first step was to reconstruct the 8 original satellite tiles they were cut from. From those, an adaptive tiling pass generated about 12k patches of 128x128. The grid is chosen per tile so patches land at or near native resolution, which minimizes bilinear downsampling. Bilinear resampling is a low-pass filter, and after compression has already eaten part of the high frequencies, you really do not want to smear what is left. The fixed 128x128 size also caps per-patch compute, which kept experimentation fast.
02Cleaning: four algorithms, four different problems
No single filter fixes compressed RGB. The restoration chain runs four stages, each aimed at one specific degradation, and the order is not negotiable.
1. Non-Local Bayes + SOS boosting, against stochastic noise. NL-Bayes (Lebrun, Buades and Morel, 2013) groups similar patches across the image and denoises each group with a Bayesian estimate under a Gaussian patch model. Satellite scenes are extremely self-similar, sand, water and rooftops repeat everywhere, so patch-based methods have a lot of material to work with. On top of it, SOS boosting (Romano and Elad, 2015), for Strengthen-Operate-Subtract, feeds the denoiser its own previous estimate added back to the input, then subtracts it. It recovers detail the first pass over-smoothed. Denoising must come first, because everything after this point, deconvolution above all, amplifies whatever noise it receives.
2. Short-term wavelet shrinkage, against residual grain. Patch models leave behind a floor of low-amplitude, spatially uncorrelated residue. A light soft-threshold on the finest wavelet scales (Donoho and Johnstone, 1994) removes it cheaply without touching structural scales. This stage is deliberately conservative; its only job is to hand a clean floor to the deconvolution.
3. Richardson-Lucy deconvolution, against optical and atmospheric blur. The imaging chain convolves the ground truth with a blur kernel: optics, atmosphere, platform motion. Richardson (1972) and Lucy (1974) give the classic iterative maximum-likelihood deblur under Poisson statistics, which matches the photon-counting nature of the sensor. RL has two famous vices: it amplifies noise, which is why stages 1 and 2 exist, and it rings at image borders, which is one more reason to deconvolve reconstructed tiles rather than small fragments. Iterations are capped early; past a handful, RL drifts toward a noise-amplifying fixed point.
4. CLAHE, for the segmenter, not for the eye. Contrast Limited Adaptive Histogram Equalization (Zuiderveld, 1994) equalizes histograms in local windows with a clip limit so flat regions do not blow up into noise. This last stage is not restoration, it is conditioning: the superpixel step ahead is driven by local color gradients, and CLAHE makes weak-but-real boundaries, a dirt road against sand, a shallow water edge, produce gradients strong enough to stop a geodesic front. It runs last so its nonlinear local remapping does not violate the statistical assumptions of the model-based stages before it.
03Superpixels: a geodesic Voronoi diagram with a shortcut
Superpixels are this pipeline's replacement for a learned encoder. They turn 16k pixels into a few hundred perceptually coherent regions, and everything downstream reasons about regions, not pixels.
The segmentation is a geodesic Voronoi diagram, following the idea of Zeng et al., Structure-Sensitive Superpixels via Geodesic Distance (ICCV 2011): each pixel belongs to the seed with the smallest geodesic distance, where the path cost accumulates local color variation along the way. Crossing an edge is expensive, so region borders snap onto image structure. A plain Euclidean Voronoi diagram, or a spatially dominated clustering like SLIC, will cut straight through a shoreline when the spatial term wins. The geodesic one will not.
The expensive part of the original paper is the loop: place seeds, compute the partition, relocate centroids, repeat until stable. Most of those iterations exist to recover from a bad start. So instead of iterating more, I initialized better:
- Quantize the cleaned patch to 16 colors, dynamically. This detail matters: a mass-weighted quantizer like median-cut allocates bins proportionally to pixel count, so a scene that is 70% sand burns ten shades of beige and leaves roads, rooftops and water edges fighting over the leftovers. The quantization here allocates bins by coverage of the color space instead, so every distinct material gets a bin regardless of how much area it occupies, and the dominant terrain collapses into a few large ones.
- Flood-fill the quantized map into connected regions of near-constant color.
- Lay a grid of seeds inside each region, with density proportional to its area.
The seeds are already partially color-aware before the first geodesic pass runs: none of them straddles a strong edge, big homogeneous regions like water get sparse coverage, busy industrial areas get dense coverage. From that starting point the partition needs dramatically fewer iterations than the uniform-grid initialization in the paper, close to a single pass in practice, and boundary adherence came out better, because iterative relocation can drift centroids next to edges and produce ragged borders.
04From regions to labels
Each superpixel is summarized by five statistics per RGB channel: min, max, mean, standard deviation and skewness. Fifteen features total. That is deliberately spartan. With 6 classes and a training set derived from 72 annotations, feature count is a regularization decision, and these five moments already separate the cartography classes well at this resolution: water is dark, low-variance and negatively skewed; sand is bright and flat; built-up areas have high variance and heavy tails from roof and shadow mixtures.
The classifier is XGBoost (Chen and Guestrin, 2016), taking the argmax over 6 classes per superpixel. Gradient-boosted trees fit this setting unusually well: they are strong on small tabular datasets, indifferent to feature scaling (the moments live on very different ranges), and their axis-aligned splits match the thresholded nature of the problem. "Mean blue below X and green variance below Y" simply is water.
05The hardest part: turning predictions back into a map
Here is the honest confession: everything up to this point is the part people expect to be hard, and it was not. The genuinely hard part starts after the classifier has spoken.
The classifier gives one prediction per superpixel, and I made that a design decision by construction: a max() over the segmentation, one label per superpixel, take the class scores, keep the argmax, done. This simplification is what makes the rest of the job tractable. Regions are internally consistent by definition, so the reconstruction problem becomes "assemble a few hundred labeled puzzle pieces per patch" instead of "reconcile 16k independent pixel opinions."
But assembling those pieces is where the real difficulty lives, and it splits into three sub-problems:
Fusing superpixels back to tile scale. Predictions are made on 128x128 patches; the deliverable is a labeled tile. Superpixel masks from adjacent patches have to be merged, adjacent regions with the same label fused into single cartographic units, and the geometry realigned from patch coordinates back to the original tile frame by inverting the adaptive tiling transform. None of this is conceptually deep, all of it is where the bugs concentrate: off-by-one alignments at patch borders show up as one-pixel label seams that look exactly like a classifier error and are not.
Smoothing the label field. The classifier is right far more often than not; mispredictions are sparse, isolated superpixels rather than large wrong areas. That sparsity is exactly what makes them fixable: classification errors are close to spatially uncorrelated, while true land cover is heavily autocorrelated. A lone "urban" region surrounded by water is almost certainly wrong. A robust spatial smoothing pass, majority reassignment of isolated regions from their neighborhood, exploits that mismatch and converts many isolated mistakes into correct labels essentially for free. It is the cheapest accuracy gain in the whole pipeline.
The road not taken: keeping the full distribution. The argmax throws information away. The alternative is to keep the entire class-probability vector per superpixel and predict all of them jointly, letting neighboring distributions negotiate: a 55/45 water-versus-sand region surrounded by confident water should tip over, a confident island of urban should survive smoothing. That negotiation is a structured, non-linear inference problem over the region adjacency graph. Classical tools exist, CRFs over superpixel graphs are the traditional answer, but to actually exploit the non-linearity of these interactions you end up wanting a learned model, in practice a DNN, something like a graph network over the superpixel adjacency structure. Which is a funny place for an AI-free pipeline to land: the classical stack carries you all the way through cleaning, segmentation and classification, and the first place a neural network genuinely earns its keep is the last mile of spatial reasoning. For this project, argmax plus robust smoothing was the right trade. For a production system, that last mile is where I would spend the learning budget.
06Ideas to push it further
A few directions I would explore next, in rough order of expected value:
- Scattering features instead of moments. The fifteen distribution moments are blind to texture, which is the known ceiling of the current features. The invariant scattering transform (Mallat, 2012; Bruna and Mallat, 2013) provides translation-invariant, deformation-stable texture descriptors with no training data, exactly the constraint here. Mallat's texture classification results with scattering networks are impressive enough that this is the single change I would make first. And if the goal shifted to production with more data, a SIFT-CNN style hybrid built on scattering is the natural next rung.
- A principled geodesic weight. The color-delta edge weight in the Voronoi diagram is a hand-tuned proxy for "structure at the scale I care about." Deriving the weights from scale-space theory (Lindeberg, 1994), or more speculatively from a renormalization-group style coarse-graining, would make them a function of the scale at which structure actually exists, instead of a knob tuned by eye, and would likely speed up the diagram computation as well.
- Retinex as a soft cleaner, early in the chain. Retinex (Land and McCann, 1971), and its multi-scale variant (Jobson et al., 1997), separates illumination from reflectance. That is precisely the satellite problem: sun angle and haze modulate the same ground truth differently across tiles, and that variance currently leaks into the color moments as inter-tile feature drift.
- Marchenko-Pastur as a spectral denoiser. Random matrix theory says the eigenvalues of a pure-noise covariance matrix fill a known bulk, the Marchenko-Pastur distribution, and anything escaping the bulk is signal. Shrinking or discarding in-bulk eigenvalues of patch covariances, with the cutoff set by Gavish and Donoho's optimal threshold, could replace the first two cleaning stages with something essentially parameter-free.
- Conditional Random Field as superpixel smoothing+reconstitution the current painpoint is argmax superpixel segment compression and its translation back to images. As shown on previous images the current method is only good with low-error rate forecast model. A CRF is a structured prediction model that lets you encode dependencies between the output variables, in this case, the class labels of adjacent superpixels. Instead of predicting each superpixel independently, a CRF jointly reasons about the entire label field, letting neighboring distributions influence each other which create an information bridge between superpixel without leaking noise. CRF example of Satellite images
07Closing thoughts
So, can we do multi-class segmentation without a DNN? On this problem, with 72 labels and compressed RGB: yes, and comfortably. Every piece of structure you inject without learning it, restoration physics, geometric priors, robust regional statistics, is a few hundred labels you did not need. The pipeline is a bet that forty years of classical vision literature is worth a few thousand annotations, and the bet paid.
It also clarified where deep learning actually belongs in this problem. Not in the encoder, the classical stack covers that surprisingly well at small scale. The one place a network genuinely earned an invitation was the structured, non-linear spatial reasoning after prediction, and for a toy project, a robust smoothing filter got close enough. Julia made the experimentation loop fast, C++ made the final version fast, and neither of them 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.