Measure the best-achievable speedups of NDVI operations with our fully-fused method, under single-threaded and multi-threaded query workloads.
NDVI is a two-band operation. I modified the single-strean codec to decode both bands in lockstep: for each 256-element sub-block, OutReg j is extracted from band A into OutReg_A and from band B into OutReg_B, then op(OutReg_A, OutReg_B) is computed immediately - all in vector registers, never materialising either band to memory.
I fuse NDVI into standalone bitpacking only (no FoR at the logical level).
The single-band sum benchmarks show standalone bitpacking gives the best speedups on Landsat & Sentinel’s NDVI-relevant bands - better than FoR+bitpacking, PFor, and FoR+PFor. Fully fusing NDVI into the more complex codecs would be significant implementation effort for diminishing returns on these datasets.
I implemented SUM(NDVI) and COUNT(NDVI > 0.3) - both common geospatial aggregations (SUM is essentially the mean; the scalar divide for mean at the end is negligible). I chose aggregations over stores because they fit naturally into the streaming framework and are cleanly end-to-end comparable.
SUM(NDVI): the divide inside NDVI cannot be factored out, so the operation must do a floating-point divide (AVX2 has no integer divide; other approaches like using libdivide requires a constant denominator).
This is a worst-case for instruction cost. Three implementations:
_mm256_div_ps - standard, high-accuracy, FP-unit bound_mm256_rcp_ps - approximate reciprocal_mm256_rcp_ps + 1 Newton-Raphson step - approximate reciprocal with accuracy correctionCOUNT(NDVI > 0.3): the divide can be eliminated algebraically:
NDVI > 0.3 ⇒ (N - R) / (N + R) > 0.3 ⇒ 0.7N - 1.3R > 0
Computable with fixed-point multiplies and a comparison mask - no float conversion or divide. Best-case for instruction cost.
As a comparison baseline, I also implemented SUM(NIR + Red) - a simple integer widen-sum representing the minimum possible instruction overhead for a two-band operation.
These span a range of instruction latencies, letting us explore how operation cost interacts with the compression speedup.
Setup is the same as the single-band case, with two extensions: (1) two grids in memory, one per band; (2) a multi-threaded benchmark (below). For fairness, the baseline (uncompressed) and experiment (fused compressed) use the same SIMD instructions for computing the operation.
Geospatial computations are commonly parallelised. As thread count increases, RAM bandwidth saturates and the benefit of a smaller compressed working set amplifies - so single-threaded is the worst case for our method.
For concurrency level X (number of reader threads), we compute the operation over X independently equally-spaced regions. Each thread reads its blocks and performs fused decode+compute independently (no concurrent writes). We measure wall-clock time for all threads to complete, aggregated appropriately.