emd-ring

Exact Earth Mover's Distance on a 9-bin ring, about 9x faster than OpenCV.

What it is

A tiny, zero-dependency Rust solver for the exact circular EMD between two 9-bin histograms. It is a drop-in match for what OpenCV's cvCalcEMD2(DIST_USER, ring-penalty) returns on the same inputs - same semantics, including unequal total masses (Rubner partial matching, cost normalized by flow) - but it runs at 149 ns/column vs OpenCV's 1,356 ns/column on a reference machine, with no allocation and no FFI.

When to use it

Use it when you compute EMD over 9-bin histograms on a ring The ground distance is the circular one, min(|i-j|, 9-|i-j|).

In practice that is the shape you get from HOG orientation histograms (9 unsigned orientation bins): comparing gradient-orientation distributions between two image patches or two consecutive video frames - motion and similarity scoring, background subtraction, and the like.

If you are calling cvCalcEMD2 in a hot loop on 9-bin circular signatures, this is a direct replacement.

use emd_ring::ring_emd_exact;

// two 9-bin histograms (non-negative bins; total mass must be > 0)
let a = [0.0, 1.0, 3.0, 2.0, 0.0, 0.0, 1.0, 0.0, 0.0];
let b = [0.0, 0.0, 2.0, 3.0, 1.0, 0.0, 0.0, 0.0, 1.0];

let cost = ring_emd_exact(&a, &b); // == cvCalcEMD2(DIST_USER, ring) for these inputs

When NOT to use it

Reach for OpenCV or a general optimal-transport solver instead when:

• The bin count is not 9. This is specialized for n = 9 and takes &[f32; 9].
• The ground distance is not the ring min(|i-j|, 9-|i-j|) - 2-D signatures, arbitrary cost matrices, and so on.
• Your reference is a different EMD. This matches one specific oracle: OpenCV's cvCalcEMD2 with DIST_USER and the ring penalty, including its Rubner unequal-mass behavior. Against a textbook balanced EMD or another library the numbers will not line up.

Preconditions (not an OpenCV fallback)

Bins must be non-negative and both sides must have total mass > ~1e-6. On input that violates this the result is unspecified - no panic, but no meaningful value. Identical inputs are fine (the answer is 0), but because that answer is trivial, callers in a hot loop usually skip both empty and identical columns before calling any solver, with a w1 > 1e-6 && w2 > 1e-6 guard, so emd-ring only ever sees columns that pass it.

Benchmarks

Median ns per column pair, criterion, on a reference x86-64 machine, fed real (non-synthetic) histogram-column pairs:

solver	ns/column	note
OpenCV cvCalcEMD2 (full path, FFI+Mat)	1,356	the thing being replaced
ring_emd_exact (this crate)	149	exact, ~9x faster

The 149 ns figure is the end of an optimization chain - 1,130 (DP reference) -> 224 (PAVA + window scan) -> 152 (SSE2 scan) -> 149 - landing in the same ballpark as a lossy O(N) approximation (~95 ns) while staying exact. Reproduce with cargo bench.

How it was validated

The oracle is OpenCV itself: cvCalcEMD2 is treated as a black box and matched, not re-derived from papers.

• Ground-truth corpus: 21.3M labeled histogram-column pairs, split into train and holdout groups so the holdout is genuinely unseen.
• 100% agreement across all 21.3M records (19.7M train + 1.6M holdout) at tolerance: relative error <= 1e-5 or absolute <= 1e-7 on the column contribution max(w1,w2) * cost. The tolerance is not tighter because cvCalcEMD2 is itself only exact to ~1.8e-4 * W - the oracle drifts more than the solver does.
• Committed 50k-record fixture - a stratified sample across mass-ratio deciles and edge cases - so the agreement test runs anywhere with no corpus present.
• Property tests (quickcheck): symmetry, identity, scale and rotation invariance, ring-diameter bound; the fast PAVA path is pinned bit-identical to a slow, obviously-correct DP reference.
• Drop-in check: used as a replacement for the OpenCV cvCalcEMD2 call, it produced byte-identical end-to-end results.

Running tests and benchmarks

# unit + property tests + committed-fixture agreement (no corpus needed)
cargo test --release

# full 21.3M-record corpus sweep (the corpus is multi-GB and not in git)
EMD_RING_CORPUS=/path/to/corpus cargo test --release

# benchmarks
cargo bench

The public API is small: ring_emd_exact(&[f32; 9], &[f32; 9]) -> f32, plus NUM_BINS and ring_distance. Zero runtime dependencies.

How it works

The speed comes from not treating this as a general optimal-transport problem. Because the bins sit on a ring, the answer is almost a formula.

Picture each histogram as piles of dirt, one per bin. EMD is the least work to reshape the first set of piles into the second, where work = (dirt moved) x (distance). On a ring, bin 0 and bin 8 are neighbors, so dirt can travel either way around the loop.

On a line, the plan is forced. Between each pair of neighbouring bins is a "gap." The dirt that has to cross a given gap is just the running total of a - b up to that point: whatever surplus has piled up on the left has nowhere to go but across that gap to the right. Add up the sizes of those gap-flows and you have the exact cost - and distance takes care of itself, because dirt that travels far is counted at every gap it passes.

   a:          3    0    0    0
   b:          0    1    1    1
   a - b:     +3   -1   -1   -1      per-bin surplus / deficit
   gap flow:   3    2    1    0      running total = dirt crossing each gap

(Three scoops leave bin 0; one drops off at each of the next bins, so 3 cross the first gap, then 2, then 1.)

The ring adds exactly one free choice: how much dirt circulates all the way around the loop. That is a single constant subtracted from every gap-flow, and the value that makes the total smallest is the median of the gap-flows - the classic "point with the least total distance to a set of numbers." So the whole equal-mass solve is: one pass for the gap-flows, take their median, subtract it, sum the sizes. No searching, no trying options.

Unequal totals - the case that matters most in practice - add a second freedom. Only the smaller pile's worth of dirt has to move; the surplus stays put, for free, wherever is cheapest to abandon it. That turns "subtract a flat median" into "subtract the best non-decreasing staircase, allowed to climb by the total surplus": each step spends a bit of surplus to drop the line where it saves the most crossings. ring_emd_exact finds that staircase with PAVA (pool-adjacent-violators isotonic regression) plus a small scan over the ~18 circulation offsets. Equal masses are just the special case where the staircase can't climb at all, so it collapses back to the flat median.

It is all fixed-size arithmetic on [f32; 9] stack arrays - no heap, no solver iterations, no FFI - with the inner scan hand-vectorized. OpenCV, by contrast, hands the same nine numbers to a general transportation-simplex LP that never learns the bins form a ring, and pays for that generality on every call. That is the ~9x.