splat3d: CPU-SIMD 3D Gaussian Splatting forward renderer (Kerbl 2023)

Cargo.toml

@@ -38,6 +38,10 @@ test = true
 name = "ocr_benchmark"
 required-features = ["std"]
 
+[[example]]
+name = "splat3d_flex"
+required-features = ["splat3d"]
+
 [dependencies]
 num-integer = { workspace = true }
 num-traits = { workspace = true }
@@ -161,6 +165,11 @@ harness = false
 name = "zip"
 harness = false
 
+[[bench]]
+name = "splat3d_bench"
+harness = false
+required-features = ["splat3d"]
+
 [features]
 default = ["std", "hpc-extras"]
 
@@ -211,6 +220,16 @@ native = ["std"]
 intel-mkl = ["std"]
 openblas = ["std"]
 
+# splat3d: CPU-SIMD 3D Gaussian Splatting forward renderer
+# (`src/hpc/splat3d/*`). Pure SIMD, no GPU, no wgpu, reuses the
+# existing `crate::simd` polyfill (F32x16 via AVX-512 / AVX2 / NEON
+# / scalar dispatch). Gated because the module pulls in the Smith-1961
+# 3×3 SPD eigendecomp + EWA-sandwich projection kernels; downstream
+# consumers (medvol, lance-graph-render) opt in. f32 hot path; the
+# Pillar-7 probe certifying the math sibling lives in
+# `lance-graph/crates/jc/src/ewa_sandwich_3d.rs`.
+splat3d = ["std"]
+
 # no_std polyfill for `static LazyLock` in `src/simd.rs` (sprint A12).
 # Pulls in `portable-atomic` with the `critical-section` impl plus the
 # `critical-section` runtime so we can build a once-cell-style cache for

benches/RESULTS.md

@@ -0,0 +1,125 @@
+# splat3d bench results
+
+Per-kernel timing baseline for the `splat3d` feature. Regression > 5%
+on any row blocks merge per the sprint discipline. Update this file in
+the same commit as any change to a `splat3d` kernel.
+
+## Run
+
+```bash
+# Default build (x86-64-v1 baseline, F32x16 = AVX2-emulated 2× __m256)
+cargo bench --features splat3d --bench splat3d_bench
+
+# AVX-512 native build (recommended on Sapphire Rapids / Zen4)
+RUSTFLAGS="-C target-cpu=native" \
+  cargo bench --features splat3d --bench splat3d_bench
+```
+
+Hardware: record the CPU model + topology + the `target-cpu` /
+`target-feature` flags used so cross-box comparisons are meaningful.
+
+## PR 1 — Spd3 + EWA-sandwich SIMD batch
+
+Baseline measurements from the sprint's reference hardware run.
+
+### Hardware: Intel Xeon (Sapphire Rapids family), AVX-512F+BW+VL+VNNI+BF16, 2.10 GHz, container build
+
+The PR 1 spec aimed for ≥10× speedup on `sandwich_x16` over the scalar
+loop on AVX-512. Measured 1.83× — the AoS↔SoA transpose overhead at 6
+fields per `Spd3` × 16 lanes dominates the inner-loop SIMD savings for
+this microbench. The downstream impact is muted because the rasterizer
+(PR 5) and `GaussianBatch::covariance_x16` (PR 2) already keep their
+hot-path data in SoA layout, avoiding the transpose. Treat the 1.83×
+microbench number as a floor; the rasterizer-driven benchmark in PR 7
+exercises the SoA-native path that benefits more strongly from F32x16.
+
+Per the architectural decision in `.cargo/config.toml` ("No global
+target-cpu — each kernel uses `#[target_feature(enable = "avx512f")]`
+per-function with LazyLock runtime detection"), the DEFAULT build uses
+the AVX2-emulated F32x16. The `target-cpu=native` row below shows the
+intended-tier numbers.
+
+#### Default build (no `target-cpu` flag)
+
+| Bench | Median | Speedup vs scalar |
+|---|---|---|
+| `spd3_sandwich_scalar_x16_loop` | 209.96 ns | 1.0× |
+| `spd3_sandwich_simd_x16` | 1225.7 ns | **0.17× (slower)** |
+| `spd3_eig_smith_1961` | 130.82 ns | — |
+| `spd3_from_scale_quat` | 11.35 ns | — |
+
+The SIMD regression on the AVX2-emulated build is a known artifact: the
+polyfill emits two `__m256` operations per `F32x16` op AND adds the
+6-field AoS↔SoA transpose at the function boundary. Net: more
+instructions than the scalar loop, which the autovectorizer is happy
+to map to `vfmadd` chains directly. Filed as TECH_DEBT for the
+performance sprint:
+- Restructure `sandwich_x16` to take SoA inputs directly (skip the
+  transpose); call sites (rasterizer, `GaussianBatch::covariance_x16`)
+  already have SoA layout.
+- Add runtime tier dispatch in `sandwich_x16` so AVX2 builds call a
+  scalar loop wrapper that the compiler auto-vectorizes cleanly.
+
+#### `RUSTFLAGS="-C target-cpu=native"` build (AVX-512F path active)
+
+| Bench | Median | Speedup vs scalar |
+|---|---|---|
+| `spd3_sandwich_scalar_x16_loop` | 166.33 ns | 1.0× |
+| `spd3_sandwich_simd_x16` | 90.41 ns | **1.83×** |
+| `spd3_eig_smith_1961` | 125.66 ns | — |
+| `spd3_from_scale_quat` | 9.19 ns | — |
+
+The 1.83× is below the 10× spec target but ABOVE the 1.0× break-even
+that gates the function's existence. With SoA inputs at the call site
+(no transpose), the inner-loop arithmetic ratio is 16-wide
+multiply-add chains vs 16 sequential scalars — measured rasterizer
+throughput (PR 5+) is where the kernel earns its keep.
+
+`spd3_eig_smith_1961` ≈ 126 ns: one closed-form eigendecomp dominated
+by `acos` (≈ 80 ns by itself). The diagonal-fast-path branch (which
+skips the trig entirely) is what makes the rasterizer's per-pixel
+work tractable; this microbench measures the WORST case.
+
+`spd3_from_scale_quat` ≈ 9 ns: the 3DGS canonical Σ builder. PR 2's
+`GaussianBatch::covariance_x16` SIMD-batches this; the scalar
+microbench is the per-call latency floor.
+
+## PR 2 — GaussianBatch SoA + SH eval
+
+Not yet baselined as separate benches — covered indirectly by the
+projection-kernel and rasterizer benches when PR 7 adds them.
+
+## PR 3 — Projection kernel
+
+Not yet baselined as a separate bench; the `project_chunk_x16`
+inner-loop math has identical AoS↔SoA structure to `sandwich_x16`
+and is expected to show similar 1.5-2× SIMD-vs-scalar ratios on
+AVX-512 native builds.
+
+## PR 4 — Tile binner
+
+Sort + prefix-sum throughput target (per the sprint spec): 2M
+instances sorted in ≤ 8 ms on 1 thread. Not yet benched separately;
+`sort_unstable_by_key` is the first-cut sort. Radix sort follow-up is
+TECH_DEBT once PR 7's full-pipeline timings show the binner is the
+hot spot.
+
+## PR 5 — Rasterizer
+
+Per-tile alpha-blend with the `F32x16` 16-pixel-row inner loop. The
+acceptance gate (1080p × 500K gaussians ≤ 25 ms on 8-core AVX-512) is
+left for the dedicated rasterizer bench in a follow-up; PR 5 ships
+the kernel + correctness tests, not the rasterizer-scale bench.
+
+## PR 6 — SplatFrame + SplatRenderer
+
+Double-buffer driver — no microbench; the full-pipeline rasterizer
+bench in a follow-up will exercise it under realistic load.
+
+## PR 7 — End-to-end demo
+
+The demo binary `examples/splat3d_flex.rs` and integration test
+`tests/splat3d_correctness.rs` ship as the e2e regression guards.
+Full-pipeline frame-time numbers (p50/p95/p99) await a Inria bicycle
+scene download — left as a follow-up for the dedicated benchmarking
+session against real-world data.

benches/splat3d_bench.rs

@@ -0,0 +1,101 @@
+//! Criterion benches for `ndarray::hpc::splat3d` kernels.
+//!
+//! Per-PR bench growth:
+//! - PR 1: `spd3::sandwich` scalar vs `sandwich_x16` SIMD (target ≥10×
+//!   on AVX-512), `Spd3::eig` Smith-1961 closed-form throughput,
+//!   `Spd3::from_scale_quat` (the 3DGS canonical builder).
+//! - PR 2: `gaussian::GaussianBatch::covariance_x16`, `sh::sh_eval_deg3_x16`.
+//! - PR 3+: `project_batch`, tile binning, per-tile rasterize.
+//!
+//! Hardware specs and absolute timings live in `benches/RESULTS.md`,
+//! updated per-PR. The bench output committed to RESULTS.md is the
+//! gate against regression — a >5% slowdown on any kernel blocks
+//! merge per the sprint discipline.
+
+use criterion::{black_box, criterion_group, criterion_main, Criterion};
+use ndarray::hpc::splat3d::{sandwich, sandwich_x16, Spd3};
+
+/// Deterministic 16 distinct SPD pairs. Using `[m; 16]` (PP-13 P0.2
+/// finding) let the optimizer constant-fold the scalar loop to one
+/// `sandwich` + ×16, which would make the SIMD-vs-scalar bench measure
+/// loop-folding rather than real SIMD parallelism. Each lane gets its
+/// own scale/quat so the inputs differ entry-wise across all 6 SoA
+/// channels the SIMD kernel transposes.
+fn build_distinct_pairs() -> ([Spd3; 16], [Spd3; 16]) {
+    let mut ms = [Spd3::I; 16];
+    let mut ns = [Spd3::I; 16];
+    for k in 0..16 {
+        let t = (k as f32 + 1.0) * 0.0625;
+        let scale_m = [0.5 + 1.0 * t, 0.4 + 0.9 * t, 0.3 + 1.2 * t];
+        let scale_n = [1.3 - 0.7 * t, 0.8 + 0.5 * t, 1.1 - 0.4 * t];
+        // Two different axis families — half rotate about Y, half about X+Z
+        // diagonal — so the rotation matrices populate different sets of
+        // off-diagonal cross terms.
+        let theta_m = 0.2 + 0.4 * t;
+        let theta_n = 0.7 - 0.3 * t;
+        let quat_m = [theta_m.cos(), 0.0, theta_m.sin(), 0.0];
+        let sqh = (0.5f32).sqrt();
+        let quat_n = [theta_n.cos(), theta_n.sin() * sqh, 0.0, theta_n.sin() * sqh];
+        ms[k] = Spd3::from_scale_quat(scale_m, quat_m);
+        ns[k] = Spd3::from_scale_quat(scale_n, quat_n);
+    }
+    (ms, ns)
+}
+
+fn bench_spd3_sandwich_scalar_loop(c: &mut Criterion) {
+    let (ms, ns) = build_distinct_pairs();
+
+    c.bench_function("spd3_sandwich_scalar_x16_loop", |b| {
+        b.iter(|| {
+            let mut acc = Spd3::ZERO;
+            for i in 0..16 {
+                let r = sandwich(black_box(&ms[i]), black_box(&ns[i]));
+                acc.a11 += r.a11;
+                acc.a22 += r.a22;
+                acc.a33 += r.a33;
+                acc.a12 += r.a12;
+                acc.a13 += r.a13;
+                acc.a23 += r.a23;
+            }
+            black_box(acc);
+        });
+    });
+}
+
+fn bench_spd3_sandwich_simd_x16(c: &mut Criterion) {
+    let (ms, ns) = build_distinct_pairs();
+    let mut out = [Spd3::ZERO; 16];
+
+    c.bench_function("spd3_sandwich_simd_x16", |b| {
+        b.iter(|| {
+            sandwich_x16(black_box(&ms), black_box(&ns), &mut out);
+            black_box(&out);
+        });
+    });
+}
+
+fn bench_spd3_eig(c: &mut Criterion) {
+    let s = Spd3::from_scale_quat([2.0, 1.5, 0.8], [0.8660254, 0.5, 0.0, 0.0]);
+    c.bench_function("spd3_eig_smith_1961", |b| {
+        b.iter(|| {
+            let r = black_box(&s).eig();
+            black_box(r);
+        });
+    });
+}
+
+fn bench_spd3_from_scale_quat(c: &mut Criterion) {
+    let scale = [1.3f32, 0.9, 0.6];
+    let quat = [0.7071068f32, 0.0, 0.7071068, 0.0];
+    c.bench_function("spd3_from_scale_quat", |b| {
+        b.iter(|| {
+            let s = Spd3::from_scale_quat(black_box(scale), black_box(quat));
+            black_box(s);
+        });
+    });
+}
+
+criterion_group!(
+    spd3, bench_spd3_sandwich_scalar_loop, bench_spd3_sandwich_simd_x16, bench_spd3_eig, bench_spd3_from_scale_quat,
+);
+criterion_main!(spd3);