diff --git a/.claude/rules/agent-cargo-hygiene.md b/.claude/rules/agent-cargo-hygiene.md new file mode 100644 index 00000000..35c1bc4f --- /dev/null +++ b/.claude/rules/agent-cargo-hygiene.md @@ -0,0 +1,33 @@ +# Agent Cargo Hygiene — one target dir, no 12× build residue + +## The problem + +When the orchestrator fans out a fleet of subagents (the Sonnet build/fix +agents), each agent that runs a full `cargo build`/`check`/`test` in its **own** +isolated working copy materialises its own `target/`. This workspace's +`target/` is ~7 GB. Twelve agents in twelve worktrees = ~84 GB of duplicated +build residue and twelve cold compiles competing for the same cores. + +## The rule + +- **Opus (orchestrator + Opus agents): run cargo freely.** No restriction. +- **Sonnet fleet agents: do NOT each run a full compile.** They edit code and + reason; they must not spawn isolated worktrees or trigger their own cold + `cargo build`/`check`/`test` that each grow a separate 7 GB `target/`. + - "tests yes, compile no": a targeted `cargo test`/`clippy` against the + **shared** workspace `target/` is fine; a bare compile-only + (`cargo check`/`build`) is wasted residue — clippy already compiles. +- **Verification is centralised.** The orchestrator (Opus) runs + `cargo fmt` + `cargo clippy` + `cargo test` **once**, in the single shared + `target/`, after the fleet's edits land. One build, not twelve. + +## How the orchestrator fans out work + +- Spawn the fleet **without** `isolation: "worktree"` so all agents share the + one repo checkout and one `target/`. +- Tell each agent explicitly: *edit only; do not run `cargo build`/`check`; do + not create a worktree; the orchestrator compiles and lints centrally.* +- After edits, the orchestrator runs the gates (`cargo fmt -p `, + `cargo clippy -p `, `cargo test -p `) — keeping the tree + `cargo clippy -- -D warnings`-clean (see `CLAUDE.md` Hard Rules) and + `cargo fmt`-clean on the pinned toolchain, with no residue blowup. diff --git a/CLAUDE.md b/CLAUDE.md index 08d80267..bb7490e4 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -21,6 +21,7 @@ This project uses specialized agents in `.claude/agents/`. Follow these rules: - Feature prioritization, gap analysis, phase planning → `l3-strategist` 4. When encountering `unsafe` code, **always** delegate to sentinel-qa for audit 5. Write decisions to the blackboard, not just to chat +6. **Cargo build residue** — fan out the Sonnet fleet in the *shared* checkout (no per-agent worktrees), edit-only; the Opus orchestrator compiles/lints/tests **once** in the single 7 GB `target/`. Opus may run cargo freely. See `.claude/rules/agent-cargo-hygiene.md`. ## Hard Rules - OpenBLAS and MKL are **mutually exclusive** feature gates. Never both. diff --git a/Cargo.lock b/Cargo.lock index e63d2546..05525ca1 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -218,6 +218,10 @@ dependencies = [ "shlex", ] +[[package]] +name = "cesium" +version = "0.0.0" + [[package]] name = "cfg-if" version = "1.0.4" diff --git a/crates/cesium/Cargo.toml b/crates/cesium/Cargo.toml new file mode 100644 index 00000000..afa31c2a --- /dev/null +++ b/crates/cesium/Cargo.toml @@ -0,0 +1,21 @@ +[package] +name = "cesium" +version = "0.0.0" +edition = "2021" +publish = false +description = "OPTIONAL Cesium/ArcGIS reference + parity oracle — NOT a production dependency. Reverse-engineers 3D Tiles / KHR_gaussian_splatting / ArcGIS-PBF into ndarray CAM SoA for parity diffing against splat3d. Dev/reference only; relocates to lance-graph/crates/cesium when flawless." + +# Workspace member via `crates/*`, deliberately NOT in `default-members` → +# optional/dev by construction (a member, never a default build target). +# RULE 3: no serde, no json crates. Prefer ArcGIS PBF (binary) over f=json. +# JSON, if unavoidable at the cold import boundary, never reaches the hotpath. + +[dependencies] +# Wired when modules go live (see src/lib.rs). The skeleton is dependency-free +# so it compiles trivially and CI stays green. +# ndarray = { workspace = true } # splat3d::GaussianBatch CAM SoA target + cam_pq +# prost / quick-protobuf = ... # ArcGIS PBF decode (binary, no-JSON) — module `arcgis_pbf` + +[features] +default = [] +# oracle = [] # gates module `oracle` (reference-render diff harness) once a real reference renderer is wired diff --git a/crates/cesium/src/arcgis_pbf.rs b/crates/cesium/src/arcgis_pbf.rs new file mode 100644 index 00000000..ece13226 --- /dev/null +++ b/crates/cesium/src/arcgis_pbf.rs @@ -0,0 +1,428 @@ +//! `arcgis_pbf` — ArcGIS FeatureService **PBF** (binary protobuf, `f=pbf`) ingestion. +//! +//! # No-JSON Mandate +//! This module is THE no-JSON path. JSON (`f=json`) is **forbidden in the hotpath** and +//! must never be introduced here. The ArcGIS FeatureService binary response (`f=pbf`) uses +//! Protocol Buffer encoding defined by Esri's open `arcgis-pbf` repository +//! (github.com/Esri/arcgis-pbf) and decodes ≥ 4× faster than the equivalent JSON at wire +//! size ≈ 30-50% smaller. The cold import boundary (URL construction, HTTP) may tolerate a +//! small JSON config blob; the frame-decode loop must not. +//! +//! # Protobuf Schema — `FeatureCollectionPBuffer` +//! Source: `github.com/Esri/arcgis-pbf/proto/FeatureCollection/FeatureCollection.proto` +//! (confirmed 2026-05-26 via raw fetch). All field numbers below are verbatim from that file. +//! +//! ```text +//! // package esriPBuffer; +//! +//! // ── Enumerations ──────────────────────────────────────────────────────────── +//! enum GeometryType { +//! esriGeometryTypePoint = 0; +//! esriGeometryTypeMultipoint = 1; +//! esriGeometryTypePolyline = 2; +//! esriGeometryTypePolygon = 3; +//! esriGeometryTypeMultipatch = 4; +//! esriGeometryTypeEnvelope = 5; +//! esriGeometryTypeNone = 127; +//! } +//! +//! enum QuantizeOriginPostion { // sic: "Postion" in the real proto +//! upperLeft = 0; +//! lowerLeft = 1; +//! } +//! +//! // ── Spatial reference ─────────────────────────────────────────────────────── +//! message SpatialReference { +//! uint32 wkid = 1; // well-known ID (horizontal) +//! uint32 lastestWkid = 2; // sic: "lastest" in the real proto (latest effective) +//! uint32 vcsWkid = 3; // vertical coordinate system WKID +//! uint32 latestVcsWkid = 4; +//! string wkt = 5; +//! string wkt2 = 6; +//! } +//! +//! // ── Quantised-coordinate transform ────────────────────────────────────────── +//! // Applied to Geometry.coords — multiply + translate to recover real-world coords. +//! message Scale { double xScale = 1; double yScale = 2; double mScale = 3; double zScale = 4; } +//! message Translate { double xTranslate = 1; double yTranslate = 2; double mTranslate = 3; double zTranslate = 4; } +//! message Transform { +//! QuantizeOriginPostion quantizeOriginPostion = 1; +//! Scale scale = 2; +//! Translate translate = 3; +//! } +//! +//! // ── Field descriptors ─────────────────────────────────────────────────────── +//! message Field { +//! string name = 1; +//! FieldType fieldType = 2; +//! string alias = 3; +//! SQLType sqlType = 4; +//! string domain = 5; +//! string defaultValue = 6; +//! } +//! +//! // ── Geometry ──────────────────────────────────────────────────────────────── +//! // `coords` is a delta-encoded, zigzag-signed sequence of quantised integers. +//! // `lengths` gives the ring/part lengths for multi-part geometries. +//! message Geometry { +//! GeometryType geometryType = 1; +//! repeated uint32 lengths = 2; // packed +//! repeated sint64 coords = 3; // packed — delta-zigzag quantised +//! repeated uint64 ids = 4; // packed — OIDs (multipoint) +//! } +//! +//! // ── Attribute value (oneof) ────────────────────────────────────────────────── +//! message Value { +//! oneof value_type { +//! string string_value = 1; +//! float float_value = 2; +//! double double_value = 3; +//! sint32 sint_value = 4; +//! uint32 uint_value = 5; +//! int64 int64_value = 6; +//! uint64 uint64_value = 7; +//! sint64 sint64_value = 8; +//! bool bool_value = 9; +//! bool null_value = 10; // true = NULL attribute +//! } +//! uint32 index = 11; // index into FeatureResult.fields for this value's column +//! } +//! +//! // ── Feature ───────────────────────────────────────────────────────────────── +//! message Feature { +//! repeated Value attributes = 1; +//! oneof compressed_geometry { +//! Geometry geometry = 2; +//! esriShapeBuffer shapeBuffer = 3; +//! CurveGeometry curveGeometry = 7; +//! } +//! Geometry centroid = 4; +//! repeated Geometry aggregateGeometries = 5; +//! Envelope envelope = 6; +//! } +//! +//! // ── Top-level result wrapper ───────────────────────────────────────────────── +//! message FeatureResult { +//! string objectIdFieldName = 1; +//! UniqueIdField uniqueIdField = 2; +//! string globalIdFieldName = 3; +//! string geohashFieldName = 4; +//! GeometryProperties geometryProperties = 5; +//! ServerGens serverGens = 6; +//! GeometryType geometryType = 7; +//! SpatialReference spatialReference = 8; +//! bool exceededTransferLimit = 9; +//! bool hasZ = 10; +//! bool hasM = 11; +//! Transform transform = 12; +//! repeated Field fields = 13; +//! repeated Value values = 14; // column-ordered attribute pool +//! repeated Feature features = 15; +//! repeated GeometryField geometryFields = 16; +//! } +//! +//! message QueryResult { +//! oneof results { +//! FeatureResult featureResult = 1; +//! CountResult countResult = 2; +//! ObjectIdsResult idsResult = 3; +//! ExtentCountResult extentCountResult = 4; +//! } +//! } +//! +//! message FeatureCollectionPBuffer { +//! string version = 1; // e.g. "2.1" +//! QueryResult queryResult = 2; +//! } +//! ``` +//! +//! # Coordinate dequantisation +//! `Geometry.coords` is a **delta-encoded, zigzag-signed** sint64 packed field. +//! The cumulative sum of decoded zigzag values gives quantised integers Q. +//! Real-world coordinates are recovered by: +//! ```text +//! x_real = (Q_x * transform.scale.xScale) + transform.translate.xTranslate +//! y_real = (Q_y * transform.scale.yScale) + transform.translate.yTranslate +//! z_real = (Q_z * transform.scale.zScale) + transform.translate.zTranslate // if hasZ +//! ``` +//! Origin convention (`QuantizeOriginPostion`): upperLeft ⟹ Y increases downward; +//! lowerLeft ⟹ standard math convention. +//! +//! # CAM SoA target +//! After decode the intermediate structs (`PbfFeatureCollection`, `PbfFeature`, `PbfGeometry`) +//! must be consumed by `to_cam_soa` (agent C) to produce `GaussianBatch` + `cam_pq`. +//! These structs do NOT survive past the import boundary. +//! +//! # Request construction (cold boundary — JSON allowed here, not in decode loop) +//! ```text +//! GET /FeatureServer/0/query +//! ?f=pbf +//! &outFields=* +//! &returnGeometry=true +//! &returnZ=true +//! &quantizationParameters={"tolerance":0.001,"extent":{...},"mode":"view"} +//! ``` +//! `quantizationParameters` is JSON but lives in the URL query string, not in the +//! binary response payload. It is constructed once at cold import, not per frame. +//! +//! # Grounding +//! - Proto field numbers: `github.com/Esri/arcgis-pbf` (fetched 2026-05-26, confirmed) +//! - R-ArcGIS/arcpbf — R reference implementation (dequant, delta decode) +//! - dfridkin/arcGIS-Rust-runtime — Rust approach patterns +//! - EsriDevEvents/locup — Location Platform Rust usage +//! - AdaWorldAPI/arcgisutils — UNVERIFIED (private / unreachable); scaffold only + +// ══════════════════════════════════════════════════════════════════════════════ +// Intermediate structs — compile-only (no live impl until Opus review) +// ══════════════════════════════════════════════════════════════════════════════ + +/// Decoded `SpatialReference` from the PBF `FeatureResult` header. +/// Passed directly to `esri_crs::EsriCrs` for reprojection before CAM SoA hand-off. +/// +/// ```rust +/// // pub struct PbfSpatialReference { +/// // /// Field 1 — original WKID (horizontal). May be legacy (e.g. 102100). +/// // pub wkid: u32, +/// // /// Field 2 — latest WKID (sic "lastest" in proto); use when non-zero. +/// // pub latest_wkid: u32, +/// // /// Field 3 — vertical CRS WKID (0 = none / assumed ellipsoidal). +/// // pub vcs_wkid: u32, +/// // /// Field 4 — latest vertical WKID. +/// // pub latest_vcs_wkid: u32, +/// // /// Field 5/6 — WKT fallback when WKID is 0; expensive, avoid in hotpath. +/// // pub wkt: Option, +/// // } +/// ``` +pub struct PbfSpatialReference; + +/// Quantised-coordinate transform decoded from `FeatureResult.transform` (field 12). +/// Must be applied to every `Geometry.coords` delta-sum before the result is meaningful. +/// +/// ```rust +/// // pub struct PbfTransform { +/// // /// upperLeft (0) or lowerLeft (1) — affects Y-axis direction. +/// // pub origin: u8, +/// // /// Multipliers: xScale, yScale, zScale (mScale unused for splats). +/// // pub scale: [f64; 4], // [x, y, m, z] +/// // /// Additive offsets after scaling. +/// // pub translate: [f64; 4], // [x, y, m, z] +/// // } +/// // +/// // impl PbfTransform { +/// // /// Apply to a zigzag-decoded, delta-accumulated quantised coordinate triple. +/// // /// `q` is [Q_x, Q_y, Q_z] after cumulative sum. +/// // #[inline] +/// // pub fn dequantise(&self, q: [i64; 3]) -> [f64; 3] { +/// // [ +/// // q[0] as f64 * self.scale[0] + self.translate[0], +/// // q[1] as f64 * self.scale[1] + self.translate[1], +/// // q[2] as f64 * self.scale[3] + self.translate[3], // z uses index 3 +/// // ] +/// // } +/// // } +/// ``` +pub struct PbfTransform; + +/// A single decoded feature geometry — positions only; attributes are stripped +/// at the import boundary (we only need xyz for CAM SoA splat seeding). +/// +/// `coords_raw` contains the raw zigzag-decoded delta values **before** dequantisation +/// so that `PbfTransform::dequantise` can be applied in one pass without per-coord +/// allocation. +/// +/// ```rust +/// // pub struct PbfGeometry { +/// // /// Geometry type — we only process Point / Multipoint / Multipatch for splats. +/// // pub geometry_type: u32, // raw proto enum value +/// // /// Delta-accumulated zigzag sint64 coords, stride = 2 (xy) or 3 (xyz, hasZ). +/// // pub coords_raw: Vec, +/// // /// Part/ring lengths for multi-part geometries. +/// // pub lengths: Vec, +/// // /// True when Z is present (from `FeatureResult.hasZ` field 10). +/// // pub has_z: bool, +/// // } +/// ``` +pub struct PbfGeometry; + +/// One decoded feature (geometry + object ID only — other attributes dropped). +/// +/// ```rust +/// // pub struct PbfFeature { +/// // pub object_id: u64, +/// // pub geometry: Option, +/// // /// Centroid geometry (field 4), present for polygon/polyline features. +/// // pub centroid: Option, +/// // } +/// ``` +pub struct PbfFeature; + +/// Top-level decoded PBF response — the unit consumed by `to_cam_soa`. +/// This struct must NOT outlive the import boundary; it is transmuted into CAM SoA +/// by agent C's `to_cam_soa::pbf_to_gaussian_batch`. +/// +/// ```rust +/// // pub struct PbfFeatureCollection { +/// // /// Protocol version string from `FeatureCollectionPBuffer.version` (field 1). +/// // pub version: String, +/// // /// Spatial reference for reprojection via `esri_crs`. +/// // pub spatial_ref: PbfSpatialReference, +/// // /// Coordinate dequantisation parameters. +/// // pub transform: PbfTransform, +/// // /// Geometry type from `FeatureResult.geometryType` (field 7). +/// // pub geometry_type: u32, +/// // /// Decoded features — geometry only, attributes stripped. +/// // pub features: Vec, +/// // /// True if the server truncated the result set. +/// // pub exceeded_transfer_limit: bool, +/// // } +/// ``` +pub struct PbfFeatureCollection; + +// ══════════════════════════════════════════════════════════════════════════════ +// Decode pipeline — commented until protobuf dep is wired +// ══════════════════════════════════════════════════════════════════════════════ + +/// Decode a raw `f=pbf` HTTP response body into a [`PbfFeatureCollection`]. +/// +/// # No-JSON guarantee +/// This function operates entirely on binary Protocol Buffer bytes. No JSON +/// parsing occurs here or in any function it calls. JSON is only permitted in +/// the URL query string at request construction time (cold boundary). +/// +/// # Wire format +/// The bytes are a length-delimited protobuf encoding of `FeatureCollectionPBuffer` +/// (package `esriPBuffer`). Decode with `prost` or `quick-protobuf` — either is +/// fine; the generated code is a compile-time artefact, not a runtime allocation. +/// +/// # Decode steps +/// 1. Parse outer `FeatureCollectionPBuffer` → extract `version` + `QueryResult`. +/// 2. Assert `QueryResult` is `featureResult` variant (field 1); others (count/ids/extent) +/// carry no geometry and should return `Err(PbfError::NoGeometry)`. +/// 3. Read `FeatureResult` fields: +/// - field 8 → `PbfSpatialReference` +/// - field 10 (`hasZ`) + field 11 (`hasM`) → stride hint +/// - field 12 → `PbfTransform` (scale + translate + origin) +/// - field 9 (`exceededTransferLimit`) +/// 4. For each `Feature` (field 15): decode `Geometry` (field 2) or `centroid` (field 4). +/// - Delta-accumulate `Geometry.coords` (packed sint64, zigzag). +/// - Apply `PbfTransform::dequantise` → `[f64; 3]` per vertex. +/// - Collect into `PbfGeometry.coords_raw` (pre-dequant) to defer allocation. +/// 5. Drop all `Value` attribute fields (field 14 / feature field 1) — not needed for splats. +/// +/// ```rust +/// // pub fn decode_pbf(bytes: &[u8]) -> Result { +/// // use prost::Message; +/// // // Step 1 — outer wrapper +/// // let pbf = esri_pb::FeatureCollectionPBuffer::decode(bytes) +/// // .map_err(PbfError::Decode)?; +/// // // Step 2 — must be featureResult +/// // let fr = pbf.query_result +/// // .and_then(|qr| qr.results) +/// // .and_then(|r| if let Results::FeatureResult(fr) = r { Some(fr) } else { None }) +/// // .ok_or(PbfError::NoGeometry)?; +/// // // Step 3 — header fields +/// // let spatial_ref = decode_spatial_ref(fr.spatial_reference); +/// // let transform = decode_transform(fr.transform); +/// // let has_z = fr.has_z; +/// // let has_m = fr.has_m; +/// // // Step 4 — features +/// // let features: Vec = fr.features.into_iter() +/// // .map(|f| decode_feature(f, has_z, has_m)) +/// // .collect::>()?; +/// // Ok(PbfFeatureCollection { +/// // version: pbf.version, +/// // spatial_ref, +/// // transform, +/// // geometry_type: fr.geometry_type as u32, +/// // features, +/// // exceeded_transfer_limit: fr.exceeded_transfer_limit, +/// // }) +/// // } +/// ``` +pub fn decode_pbf(_bytes: &[u8]) -> Result { + Err(PbfError::NotImplemented("scaffold only — not yet implemented")) +} + +/// Zigzag-decode + delta-accumulate a packed sint64 coordinate sequence from a +/// `Geometry.coords` field. +/// +/// The protobuf packed sint64 encoding stores each value as its zigzag-mapped uint64: +/// `zigzag(n) = (n << 1) ^ (n >> 63)`. After decoding each zigzag value the running +/// cumulative sum gives the actual quantised integer. +/// +/// ```rust +/// // fn decode_delta_coords(raw_zigzag: &[u64]) -> Vec { +/// // let mut acc = 0i64; +/// // raw_zigzag.iter().map(|&z| { +/// // let delta = ((z >> 1) as i64) ^ -((z & 1) as i64); // zigzag decode +/// // acc = acc.wrapping_add(delta); +/// // acc +/// // }).collect() +/// // } +/// ``` +/// +/// # Stride +/// After delta accumulation the flat `Vec` is strided: +/// - `hasZ = false`: stride 2, layout `[x0,y0, x1,y1, …]` +/// - `hasZ = true`: stride 3, layout `[x0,y0,z0, x1,y1,z1, …]` +/// +/// Dequantise with [`PbfTransform`] before passing to `to_cam_soa`. +pub(crate) fn decode_delta_coords(_raw: &[u64], _has_z: bool) -> Vec { + unimplemented!("scaffold only") +} + +// ══════════════════════════════════════════════════════════════════════════════ +// Error type +// ══════════════════════════════════════════════════════════════════════════════ + +/// Errors produced by the PBF decode pipeline. +/// +/// ```rust +/// // #[derive(Debug)] +/// // pub enum PbfError { +/// // /// Protobuf wire-format decode failure (wraps prost::DecodeError). +/// // Decode(prost::DecodeError), +/// // /// The QueryResult did not contain a FeatureResult with geometry. +/// // NoGeometry, +/// // /// coord count is not divisible by stride (corrupted payload). +/// // CoordStrideMismatch { count: usize, stride: usize }, +/// // /// Transform scale fields are zero — division by zero guard. +/// // DegenerateTransform, +/// // } +/// ``` +#[derive(Debug)] +pub enum PbfError { + /// Function is a scaffold stub and has not been implemented yet. + NotImplemented(&'static str), +} + +impl std::fmt::Display for PbfError { + fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { + match self { + PbfError::NotImplemented(msg) => write!(f, "not implemented: {}", msg), + } + } +} + +impl std::error::Error for PbfError {} + +// ══════════════════════════════════════════════════════════════════════════════ +// Internal helpers — commented stubs +// ══════════════════════════════════════════════════════════════════════════════ + +// fn decode_spatial_ref(sr: Option) -> PbfSpatialReference { ... } +// fn decode_transform(t: Option) -> PbfTransform { ... } +// fn decode_feature(f: esri_pb::Feature, has_z: bool, has_m: bool) +// -> Result { ... } +// +// NOTE(grounding): `esri_pb` refers to the prost-generated module from +// `crates/cesium/proto/FeatureCollection.proto` (copy of Esri/arcgis-pbf). +// Do NOT hand-write the proto structs — generate them with `prost-build` in build.rs. +// The build.rs compileProtos call is also commented until dep is wired: +// +// fn main() { +// prost_build::compile_protos( +// &["proto/FeatureCollection/FeatureCollection.proto"], +// &["proto/"], +// ).unwrap(); +// } diff --git a/crates/cesium/src/esri_crs.rs b/crates/cesium/src/esri_crs.rs new file mode 100644 index 00000000..86464662 --- /dev/null +++ b/crates/cesium/src/esri_crs.rs @@ -0,0 +1,358 @@ +//! `esri_crs` — `ESRI_crs` extension: WKID (horizontal + vertical) → WGS84 / local-origin. +//! +//! # What is ESRI_crs? +//! `ESRI_crs` is an Esri-proprietary 3D Tiles extension that attaches a full +//! coordinate reference system definition to a tileset or glTF asset via Well-Known IDs +//! (WKIDs). It is required for Gaussian splat layers produced by ArcGIS Pro / ArcGIS +//! Online — the tileset will not reproject correctly without it. The extension appears +//! in `tileset.json` under `extensions.ESRI_crs` and/or in per-tile glTF under +//! `extensions.ESRI_crs` at the asset level. +//! +//! # JSON schema (UNVERIFIED — grounded against I3S spec + ArcGIS Pro docs; no +//! dedicated ESRI_crs glTF schema URL was publicly reachable as of 2026-05-26) +//! ```json +//! { +//! "extensions": { +//! "ESRI_crs": { +//! "wkid": 103142, // original horizontal WKID +//! "latestWkid": 6565, // current WKID for same CRS (may differ from wkid) +//! "vcsWkid": 105703, // vertical CRS WKID (0 or absent = ellipsoidal) +//! "latestVcsWkid": 6360 // current vertical WKID +//! } +//! } +//! } +//! ``` +//! Sources: Esri I3S-spec `spatialReference.cmn.md` (confirmed); ArcGIS Pro +//! "Work with Gaussian splat layers" page (confirmed ESRI_crs is required); +//! Web-map specification `spatialReference` object (confirmed wkid/latestWkid/vcsWkid). +//! +//! # WKID semantics +//! - `wkid` — original ID assigned when the CRS was created. May be a deprecated alias +//! (e.g. 102100 = old Web Mercator; latestWkid = 3857). Always prefer `latestWkid` +//! when non-zero. +//! - `latestWkid` — current EPSG/Esri code. Use for projection library lookups. +//! - `vcsWkid` — vertical CRS (0 = not specified / use ellipsoidal height). +//! - `latestVcsWkid` — current vertical WKID. +//! +//! # Common WKIDs for geospatial splats +//! | WKID | latestWkid | Meaning | +//! |-------|------------|----------------------------------| +//! | 4326 | 4326 | WGS84 geographic (lat/lon/ellht) | +//! | 102100| 3857 | Web Mercator (metres, XY only) | +//! | 4978 | 4978 | WGS84 geocentric ECEF | +//! | 4979 | 4979 | WGS84 geographic 3D | +//! +//! # Reprojection pipeline +//! This module only carries the parsed CRS + coordinate-transform math stubs. +//! Full PROJ/GDAL integration is a live-impl decision deferred to Opus review. +//! The intermediate output `EsriCrsReprojected` feeds `to_cam_soa` (agent C). +//! +//! ## WGS84 target +//! All coordinates must ultimately reach WGS84 geographic (EPSG:4979) or local-ENU +//! (east-north-up) before CAM SoA hand-off. Steps: +//! 1. Parse `ESRI_crs` → `EsriCrs`. +//! 2. If `wkid == 4326` or `latestWkid == 4326` (or 4979): identity pass. +//! 3. If `wkid == 102100 / 3857` (Web Mercator): +//! apply inverse Mercator projection → WGS84 lat/lon. +//! 4. Otherwise: invoke PROJ pipeline `EPSG:{latestWkid} → EPSG:4979`. +//! (UNVERIFIED: exact PROJ string; mark as PROJ-required at runtime) +//! 5. Apply vertical transform if `vcsWkid != 0`: +//! `EPSG:{latestVcsWkid} → ellipsoidal height`. +//! (UNVERIFIED: geoid separation lookup not implemented in scaffold) +//! 6. Optionally rebase to local ENU origin (provided externally by caller). +//! +//! ## Web-Mercator inverse (inline, no PROJ dependency) +//! For WKID 3857 / 102100 this can be done with closed-form math (no external dep): +//! ```text +//! λ (lon, rad) = x / R_EARTH +//! φ (lat, rad) = 2·atan(exp(y / R_EARTH)) - π/2 +//! R_EARTH = 6_378_137.0 // WGS84 semi-major axis, metres +//! ``` +//! +//! # Grounding +//! - I3S-spec `spatialReference.cmn.md` (github.com/Esri/i3s-spec, confirmed 2026-05-26) +//! - ArcGIS Web Map spec `spatialReference` object (confirmed wkid/latestWkid/vcsWkid) +//! - ArcGIS PBF proto `SpatialReference` message (field numbers confirmed 2026-05-26) +//! - "ESRI_crs required for Gaussian splat layers": ArcGIS Pro docs (confirmed intent; +//! exact JSON schema not independently fetched — UNVERIFIED for glTF-specific form) + +// ══════════════════════════════════════════════════════════════════════════════ +// Constants +// ══════════════════════════════════════════════════════════════════════════════ + +/// WGS84 semi-major axis in metres — used in Web Mercator inverse projection. +/// +/// ```rust +/// // pub const WGS84_A: f64 = 6_378_137.0; +/// ``` +pub const WGS84_A: f64 = 6_378_137.0; + +/// WKID for WGS84 geographic 2D (EPSG:4326). +pub const WKID_WGS84_GEO2D: u32 = 4326; + +/// WKID for WGS84 geographic 3D / lat-lon-ellht (EPSG:4979). +pub const WKID_WGS84_GEO3D: u32 = 4979; + +/// WKID for WGS84 geocentric ECEF (EPSG:4978). +pub const WKID_WGS84_ECEF: u32 = 4978; + +/// Legacy Web Mercator WKID (Esri alias for EPSG:3857). +pub const WKID_WEB_MERCATOR_LEGACY: u32 = 102100; + +/// Current Web Mercator WKID (EPSG:3857). +pub const WKID_WEB_MERCATOR: u32 = 3857; + +// ══════════════════════════════════════════════════════════════════════════════ +// Parsed ESRI_crs +// ══════════════════════════════════════════════════════════════════════════════ + +/// Parsed `ESRI_crs` extension object. +/// +/// Consumed from the cold JSON import boundary (tileset.json or glTF asset extension). +/// JSON parsing is allowed here — this is the cold boundary, not the hotpath. +/// +/// ```rust +/// // #[derive(Debug, Clone)] +/// // pub struct EsriCrs { +/// // /// Original horizontal WKID (may be legacy alias). +/// // pub wkid: u32, +/// // /// Latest/current horizontal WKID; use for projection if non-zero. +/// // pub latest_wkid: u32, +/// // /// Vertical CRS WKID; 0 = ellipsoidal height (no conversion needed). +/// // pub vcs_wkid: u32, +/// // /// Latest vertical CRS WKID. +/// // pub latest_vcs_wkid: u32, +/// // } +/// // +/// // impl EsriCrs { +/// // /// Return the effective horizontal WKID (prefer latestWkid over wkid). +/// // pub fn effective_h_wkid(&self) -> u32 { +/// // if self.latest_wkid != 0 { self.latest_wkid } else { self.wkid } +/// // } +/// // +/// // /// Return the effective vertical WKID (0 = none). +/// // pub fn effective_v_wkid(&self) -> u32 { +/// // if self.latest_vcs_wkid != 0 { self.latest_vcs_wkid } else { self.vcs_wkid } +/// // } +/// // +/// // /// True when the horizontal CRS is already WGS84 (identity pass). +/// // pub fn is_wgs84(&self) -> bool { +/// // let h = self.effective_h_wkid(); +/// // matches!(h, WKID_WGS84_GEO2D | WKID_WGS84_GEO3D | WKID_WGS84_ECEF) +/// // } +/// // +/// // /// True when this is Web Mercator (fast closed-form inverse available). +/// // pub fn is_web_mercator(&self) -> bool { +/// // let h = self.effective_h_wkid(); +/// // matches!(h, WKID_WEB_MERCATOR | WKID_WEB_MERCATOR_LEGACY) +/// // } +/// // } +/// ``` +pub struct EsriCrs; + +/// Construct an `EsriCrs` from the four WKID fields that appear in both the +/// ArcGIS PBF `SpatialReference` proto message and the `ESRI_crs` JSON extension. +/// +/// This is the bridge between the PBF decode path and the CRS handling path — +/// `PbfSpatialReference` maps 1:1 onto this struct. +/// +/// ```rust +/// // pub fn from_pbf_spatial_ref( +/// // wkid: u32, +/// // latest_wkid: u32, +/// // vcs_wkid: u32, +/// // latest_vcs_wkid: u32, +/// // ) -> EsriCrs { +/// // EsriCrs { wkid, latest_wkid, vcs_wkid, latest_vcs_wkid } +/// // } +/// ``` +pub(crate) fn from_pbf_spatial_ref(_wkid: u32, _latest_wkid: u32, _vcs_wkid: u32, _latest_vcs_wkid: u32) -> EsriCrs { + unimplemented!("scaffold only") +} + +/// Parse `ESRI_crs` from a raw JSON value at the cold import boundary. +/// JSON is permitted here — this is NOT in the hotpath. +/// +/// Expected JSON shape (see module-level doc for schema): +/// ```json +/// { "wkid": 4326, "latestWkid": 4326, "vcsWkid": 0, "latestVcsWkid": 0 } +/// ``` +/// +/// ```rust +/// // pub fn from_json_value(v: &serde_json::Value) -> Result { +/// // // NOTE: serde_json import is ONLY in this cold-boundary function. +/// // // Never call this from the frame-decode loop. +/// // Ok(EsriCrs { +/// // wkid: v["wkid"].as_u64().unwrap_or(0) as u32, +/// // latest_wkid: v["latestWkid"].as_u64().unwrap_or(0) as u32, +/// // vcs_wkid: v["vcsWkid"].as_u64().unwrap_or(0) as u32, +/// // latest_vcs_wkid: v["latestVcsWkid"].as_u64().unwrap_or(0) as u32, +/// // }) +/// // } +/// ``` +// UNVERIFIED: field-name casing in the glTF ESRI_crs extension JSON may differ from +// the PBF proto (proto uses "lastestWkid" — sic — vs JSON uses "latestWkid"). +// Confirm against an actual tileset.json produced by ArcGIS Pro before wiring. +pub fn from_json_value(_v: &()) -> Result { + Err(CrsError::NotImplemented("scaffold only — not yet implemented")) +} + +// ══════════════════════════════════════════════════════════════════════════════ +// Reprojection output +// ══════════════════════════════════════════════════════════════════════════════ + +/// Coordinates reprojected to WGS84 geographic 3D (EPSG:4979): lon, lat, ellht. +/// This is the intermediate consumed by `to_cam_soa` before local-ENU rebasing. +/// +/// ```rust +/// // pub struct Wgs84Coord { +/// // pub lon_deg: f64, // longitude, degrees +/// // pub lat_deg: f64, // latitude, degrees +/// // pub ellht_m: f64, // ellipsoidal height, metres +/// // } +/// ``` +pub struct Wgs84Coord; + +/// Local-ENU rebased coordinate (east, north, up in metres relative to origin). +/// Produced by `rebase_to_local_enu` after WGS84 reprojection. +/// Consumed by `to_cam_soa` as the final position passed into `GaussianBatch`. +/// +/// ```rust +/// // pub struct EnuCoord { +/// // pub east_m: f32, +/// // pub north_m: f32, +/// // pub up_m: f32, +/// // } +/// ``` +pub struct EnuCoord; + +// ══════════════════════════════════════════════════════════════════════════════ +// Reprojection functions +// ══════════════════════════════════════════════════════════════════════════════ + +/// Reproject a batch of XYZ coordinates from the CRS described by `crs` to WGS84 (EPSG:4979). +/// +/// `coords` is a flat `[x0,y0,z0, x1,y1,z1, …]` slice in source-CRS units. +/// Returns a `Vec` of length `coords.len() / 3`. +/// +/// Dispatch logic: +/// 1. `crs.is_wgs84()` → identity (copy to `Wgs84Coord` with lon=x, lat=y, h=z) +/// 2. `crs.is_web_mercator()` → closed-form inverse Mercator (no external dep) +/// 3. Otherwise → UNVERIFIED: requires PROJ runtime; return `Err(CrsError::ProjRequired)` +/// +/// ```rust +/// // pub fn reproject_to_wgs84( +/// // crs: &EsriCrs, +/// // coords: &[f64], +/// // ) -> Result, CrsError> { +/// // assert_eq!(coords.len() % 3, 0, "coords must be xyz triples"); +/// // let n = coords.len() / 3; +/// // let mut out = Vec::with_capacity(n); +/// // if crs.is_wgs84() { +/// // for i in 0..n { +/// // out.push(Wgs84Coord { +/// // lon_deg: coords[i*3], +/// // lat_deg: coords[i*3+1], +/// // ellht_m: coords[i*3+2], +/// // }); +/// // } +/// // } else if crs.is_web_mercator() { +/// // for i in 0..n { +/// // let (lon, lat) = inverse_mercator(coords[i*3], coords[i*3+1]); +/// // out.push(Wgs84Coord { lon_deg: lon, lat_deg: lat, ellht_m: coords[i*3+2] }); +/// // } +/// // } else { +/// // return Err(CrsError::ProjRequired(crs.effective_h_wkid())); +/// // } +/// // Ok(out) +/// // } +/// ``` +pub fn reproject_to_wgs84(_crs: &EsriCrs, _coords: &[f64]) -> Result, CrsError> { + Err(CrsError::NotImplemented("scaffold only — not yet implemented")) +} + +/// Closed-form inverse Web Mercator: (x_m, y_m) → (lon_deg, lat_deg). +/// Valid for EPSG:3857 / WKID 102100. No external dependency. +/// +/// ```rust +/// // #[inline] +/// // fn inverse_mercator(x: f64, y: f64) -> (f64, f64) { +/// // let lon = x / WGS84_A; +/// // let lat = 2.0 * (y / WGS84_A).exp().atan() - std::f64::consts::FRAC_PI_2; +/// // (lon.to_degrees(), lat.to_degrees()) +/// // } +/// ``` +pub(crate) fn inverse_mercator(_x: f64, _y: f64) -> (f64, f64) { + unimplemented!("scaffold only") +} + +/// Rebase a batch of WGS84 coordinates to a local ENU frame centred on `origin`. +/// +/// `origin` is the WGS84 reference point (lon_deg, lat_deg, ellht_m) that maps to (0,0,0). +/// Uses the standard WGS84 → ECEF → ENU rotation (3×3 matrix, closed-form). +/// +/// ```rust +/// // pub fn rebase_to_local_enu( +/// // coords: &[Wgs84Coord], +/// // origin: &Wgs84Coord, +/// // ) -> Vec { +/// // // 1. Convert origin to ECEF +/// // // 2. Build ENU rotation matrix from origin lat/lon +/// // // 3. For each coord: ECEF → subtract origin ECEF → rotate → ENU +/// // // UNVERIFIED: geoid separation for vertical not applied here; +/// // // caller must handle vcsWkid != 0 separately. +/// // todo!() +/// // } +/// // +/// // fn wgs84_to_ecef(lon_deg: f64, lat_deg: f64, h_m: f64) -> [f64; 3] { +/// // let a = WGS84_A; +/// // let f = 1.0 / 298.257_223_563; // WGS84 flattening +/// // let e2 = 2.0*f - f*f; +/// // let lat = lat_deg.to_radians(); +/// // let lon = lon_deg.to_radians(); +/// // let n = a / (1.0 - e2 * lat.sin().powi(2)).sqrt(); +/// // [(n + h_m)*lat.cos()*lon.cos(), +/// // (n + h_m)*lat.cos()*lon.sin(), +/// // (n*(1.0-e2) + h_m)*lat.sin()] +/// // } +/// ``` +pub(crate) fn rebase_to_local_enu(_coords: &[Wgs84Coord], _origin: &Wgs84Coord) -> Vec { + unimplemented!("scaffold only") +} + +// ══════════════════════════════════════════════════════════════════════════════ +// Error type +// ══════════════════════════════════════════════════════════════════════════════ + +/// Errors produced by the CRS reprojection pipeline. +/// +/// ```rust +/// // #[derive(Debug)] +/// // pub enum CrsError { +/// // /// JSON parsing failed at the cold import boundary. +/// // JsonParse(String), +/// // /// WKID requires PROJ runtime (not implemented in scaffold). +/// // ProjRequired(u32), +/// // /// WKID is unknown / unrecognised. +/// // UnknownWkid(u32), +/// // /// Input coordinate slice is not a multiple of 3. +/// // BadCoordStride(usize), +/// // /// Vertical CRS conversion requires geoid grid (not yet wired). +/// // VerticalCrsRequired(u32), +/// // } +/// ``` +#[derive(Debug)] +pub enum CrsError { + /// Function is a scaffold stub and has not been implemented yet. + NotImplemented(&'static str), +} + +impl std::fmt::Display for CrsError { + fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { + match self { + CrsError::NotImplemented(msg) => write!(f, "not implemented: {}", msg), + } + } +} + +impl std::error::Error for CrsError {} diff --git a/crates/cesium/src/fixtures.rs b/crates/cesium/src/fixtures.rs new file mode 100644 index 00000000..2929ef9d --- /dev/null +++ b/crates/cesium/src/fixtures.rs @@ -0,0 +1,476 @@ +//! `fixtures` (group D) — golden scenes + parity gates. +//! +//! Defines the canonical test scenes used to validate `splat3d` parity against +//! reference renders, and the SSIM / PSNR thresholds that constitute a "pass". +//! +//! # Philosophy +//! +//! The oracle **measures** parity; this module **gates** it. A `ParityGate` +//! bundles a scene descriptor with its acceptable metric bounds. Whether to +//! treat a failed gate as a CI hard-failure is the caller's decision. +//! +//! # First reference scene +//! +//! The first reference is an Inria `.ply` scene, loaded via +//! `splat3d::ply::read_ply` (confirmed: `src/hpc/splat3d/ply.rs` line 110). +//! The `staehlli/3D-settlement-development` scene is the named fixture source +//! from `lib.rs` grounding references. +//! +//! # Threshold rationale +//! +//! | Gate | SSIM | PSNR (dB) | Rationale | +//! |---------------------|-------|-----------|------------------------------------------------| +//! | `SYNTHETIC_TIGHT` | 0.99 | 40.0 | Synthetic scene, controlled camera, no sort risk | +//! | `INRIA_PLY_NOMINAL` | 0.95 | 30.0 | Real scene; sort-order risk may cost ~3–5 dB | +//! | `SMOKE_TEST` | 0.80 | 20.0 | Minimal bar: confirms pipeline runs at all | +//! +//! # Sort-order risk (see `oracle` module doc) +//! +//! All thresholds above are set conservatively to account for the front-to-back +//! vs back-to-front accumulation delta documented in `oracle.rs`. The +//! `INRIA_PLY_NOMINAL` gate in particular leaves 5–10 dB of headroom that a +//! sort-order-aware comparison could reclaim. +//! +//! # Implementation status +//! +//! **Commented scaffold only.** The `Fixture` and `ParityGate` structs are +//! live `std`-only. All references to `splat3d` symbols are commented out +//! pending the ndarray dep re-enable + Opus + CodeRabbit review. +//! +//! # Source grounding +//! +//! Symbols cited (confirmed from local source): +//! - `read_ply` — `src/hpc/splat3d/ply.rs:110` +//! - Framebuffer layout — `raster.rs` doc: interleaved RGB, len = `3·W·H` +//! - `GaussianBatch::len` — `gaussian.rs:89` (active gaussian count) + +// ───────────────────────────────────────────────────────────────────────────── +// Parity threshold constants (live) +// ───────────────────────────────────────────────────────────────────────────── + +/// SSIM / PSNR pair for a parity gate. +/// +/// # Examples +/// +/// ``` +/// use cesium::fixtures::Thresholds; +/// +/// let t = Thresholds::SYNTHETIC_TIGHT; +/// assert_eq!(t.min_ssim, 0.99); +/// assert_eq!(t.min_psnr_db, 40.0); +/// +/// // Thresholds are ordered: SMOKE < INRIA < SYNTHETIC +/// assert!(Thresholds::SMOKE_TEST.min_ssim < Thresholds::INRIA_PLY_NOMINAL.min_ssim); +/// assert!(Thresholds::INRIA_PLY_NOMINAL.min_ssim < Thresholds::SYNTHETIC_TIGHT.min_ssim); +/// ``` +#[derive(Debug, Clone, Copy, PartialEq)] +pub struct Thresholds { + /// Minimum acceptable SSIM in `[0, 1]`. + pub min_ssim: f32, + /// Minimum acceptable PSNR in dB. + pub min_psnr_db: f32, +} + +impl Thresholds { + /// Synthetic scene, controlled camera, minimal sort-order risk. + pub const SYNTHETIC_TIGHT: Self = Self { + min_ssim: 0.99, + min_psnr_db: 40.0, + }; + + /// Real Inria PLY scene; sort-order risk may cost ~3–5 dB of PSNR. + /// Threshold is set conservatively to account for front-to-back vs + /// back-to-front accumulation divergence (see `oracle` module doc). + pub const INRIA_PLY_NOMINAL: Self = Self { + min_ssim: 0.95, + min_psnr_db: 30.0, + }; + + /// Smoke test: confirms the pipeline runs and produces a plausible image. + pub const SMOKE_TEST: Self = Self { + min_ssim: 0.80, + min_psnr_db: 20.0, + }; +} + +// ───────────────────────────────────────────────────────────────────────────── +// Fixture descriptor (live — std-only) +// ───────────────────────────────────────────────────────────────────────────── + +/// Descriptor for a golden test scene. +/// +/// A `Fixture` names a scene and records its provenance. It does NOT hold the +/// actual scene data (too large for in-memory constants); the oracle pipeline +/// loads data from disk at test time. +/// +/// # Examples +/// +/// ``` +/// use cesium::fixtures::scenes; +/// +/// let f = scenes::SYNTHETIC_UNIT.clone(); +/// assert_eq!(f.id, "synthetic_unit"); +/// assert_eq!(f.expected_gaussian_count, Some(4)); +/// assert_eq!(f.width, 64); +/// assert_eq!(f.height, 64); +/// ``` +#[derive(Debug, Clone)] +pub struct Fixture { + /// Short identifier used in test output (e.g. `"settlement_dev"`). + pub id: &'static str, + /// Human-readable description of the scene and its provenance. + pub description: &'static str, + /// Expected number of Gaussians in the loaded batch (informational). + /// `None` means unchecked. + pub expected_gaussian_count: Option, + /// Image width used for the reference render, in pixels. + pub width: u32, + /// Image height used for the reference render, in pixels. + pub height: u32, + /// Background colour `[R, G, B]` in `[0, 1]` used for the reference render. + pub background: [f32; 3], +} + +/// Named fixtures. Path resolution is left to the test harness so the crate +/// has no compile-time dependency on the filesystem layout. +pub mod scenes { + use super::Fixture; + + /// Inria 3DGS garden scene (small variant, ~300K Gaussians). + /// Source: Kerbl et al. 2023 supplemental material. + // UNVERIFIED: exact gaussian count — 300K is an estimate from the paper. + pub const INRIA_GARDEN_SMALL: Fixture = Fixture { + id: "inria_garden_small", + description: "Inria 3DGS garden scene (small variant), Kerbl 2023. \ + Loaded via splat3d::ply::read_ply. \ + // UNVERIFIED: ~300K gaussians.", + expected_gaussian_count: None, // UNVERIFIED + width: 1920, + height: 1080, + background: [0.0, 0.0, 0.0], + }; + + /// staehlli/3D-settlement-development scene (named grounding ref in lib.rs). + /// This is the primary parity fixture for the cesium crate. + // UNVERIFIED: gaussian count and exact image dimensions for this scene. + pub const SETTLEMENT_DEV: Fixture = Fixture { + id: "settlement_dev", + description: "staehlli/3D-settlement-development 3DGS scene. \ + Grounding reference named in cesium/src/lib.rs. \ + Loaded via splat3d::ply::read_ply. \ + // UNVERIFIED: gaussian count and image dims.", + expected_gaussian_count: None, // UNVERIFIED + width: 1920, + height: 1080, + background: [0.0, 0.0, 0.0], + }; + + /// Synthetic unit scene: 4 Gaussians at known positions, controlled camera. + /// Used for `SYNTHETIC_TIGHT` gate — no sort-order ambiguity. + pub const SYNTHETIC_UNIT: Fixture = Fixture { + id: "synthetic_unit", + description: "4 Gaussians at axis-aligned positions, identity camera. \ + Fully deterministic: no sort-order risk, SSIM ≥ 0.99 expected.", + expected_gaussian_count: Some(4), + width: 64, + height: 64, + background: [0.0, 0.0, 0.0], + }; +} + +// ───────────────────────────────────────────────────────────────────────────── +// ParityGate (live — std-only) +// ───────────────────────────────────────────────────────────────────────────── + +/// A parity gate bundles a [`Fixture`] with its [`Thresholds`]. +/// +/// Use [`ParityGate::check`] to evaluate a [`crate::oracle::ParityMetrics`] +/// against the gate. The gate reports pass/fail but does not panic — the +/// caller decides whether a failed gate is a CI hard failure. +/// +/// # Examples +/// +/// ``` +/// use cesium::fixtures::gates; +/// +/// let gate = gates::synthetic_unit_tight(); +/// assert_eq!(gate.fixture.id, "synthetic_unit"); +/// assert_eq!(gate.thresholds.min_ssim, 0.99); +/// ``` +#[derive(Debug, Clone)] +pub struct ParityGate { + /// The scene this gate applies to. + pub fixture: Fixture, + /// SSIM / PSNR thresholds. + pub thresholds: Thresholds, +} + +/// Result of evaluating a parity gate. +/// +/// # Examples +/// +/// ``` +/// use cesium::fixtures::{gates, GateResult}; +/// use cesium::oracle::ParityMetrics; +/// +/// let fb: Vec = vec![0.5, 0.3, 0.8, 0.1, 0.9, 0.2]; +/// let metrics = ParityMetrics::compute(&fb, &fb).unwrap(); +/// let gate = gates::settlement_dev_smoke(); +/// let result: GateResult = gate.check(&metrics); +/// assert!(result.passed); +/// assert!(result.ssim_margin >= 0.0); +/// assert!(result.psnr_margin_db.is_infinite()); +/// ``` +#[derive(Debug, Clone, PartialEq)] +pub struct GateResult { + /// Whether the gate passed. + pub passed: bool, + /// Observed SSIM. + pub ssim: f32, + /// Observed PSNR in dB. + pub psnr: f32, + /// SSIM margin (observed − minimum; negative = failing). + pub ssim_margin: f32, + /// PSNR margin in dB (observed − minimum; negative = failing). + pub psnr_margin_db: f32, +} + +impl ParityGate { + /// Evaluate the gate against observed metrics. + /// + /// Returns a [`GateResult`] with `passed = true` when both `ssim` and + /// `psnr` meet the gate's [`Thresholds`]. Identical buffers produce + /// SSIM = 1.0 and infinite PSNR, which always passes. + /// + /// # Examples + /// + /// ``` + /// use cesium::fixtures::gates; + /// use cesium::oracle::ParityMetrics; + /// + /// // Two identical small framebuffers (3 pixels × RGB = 9 samples). + /// let fb: Vec = vec![0.2, 0.5, 0.8, 0.1, 0.4, 0.9, 0.3, 0.6, 0.7]; + /// let metrics = ParityMetrics::compute(&fb, &fb).unwrap(); + /// + /// // Identical buffers → SSIM 1.0, infinite PSNR → passes even the tight gate. + /// let gate = cesium::fixtures::gates::synthetic_unit_tight(); + /// let result = gate.check(&metrics); + /// assert!(result.passed); + /// assert_eq!(result.ssim, 1.0); + /// assert!(result.psnr.is_infinite()); + /// ``` + pub fn check(&self, metrics: &crate::oracle::ParityMetrics) -> GateResult { + let ssim_margin = metrics.ssim - self.thresholds.min_ssim; + let psnr_margin_db = if metrics.psnr.is_infinite() { + f32::INFINITY + } else { + metrics.psnr - self.thresholds.min_psnr_db + }; + let passed = metrics.passes(self.thresholds.min_ssim, self.thresholds.min_psnr_db); + GateResult { + passed, + ssim: metrics.ssim, + psnr: metrics.psnr, + ssim_margin, + psnr_margin_db, + } + } +} + +/// Pre-built gates for the named scenes. +pub mod gates { + use super::{scenes, ParityGate, Thresholds}; + + /// Smoke test gate for the settlement_dev scene. + /// + /// Uses [`Thresholds::SMOKE_TEST`] (SSIM ≥ 0.80, PSNR ≥ 20 dB). + /// + /// # Examples + /// + /// ``` + /// use cesium::fixtures::gates; + /// + /// let gate = gates::settlement_dev_smoke(); + /// assert_eq!(gate.fixture.id, "settlement_dev"); + /// assert_eq!(gate.thresholds.min_ssim, 0.80); + /// assert_eq!(gate.thresholds.min_psnr_db, 20.0); + /// ``` + pub fn settlement_dev_smoke() -> ParityGate { + ParityGate { + fixture: scenes::SETTLEMENT_DEV.clone(), + thresholds: Thresholds::SMOKE_TEST, + } + } + + /// Nominal gate for the settlement_dev scene (accounts for sort-order risk). + /// + /// Uses [`Thresholds::INRIA_PLY_NOMINAL`] (SSIM ≥ 0.95, PSNR ≥ 30 dB). + /// + /// # Examples + /// + /// ``` + /// use cesium::fixtures::gates; + /// + /// let gate = gates::settlement_dev_nominal(); + /// assert_eq!(gate.fixture.id, "settlement_dev"); + /// assert_eq!(gate.thresholds.min_ssim, 0.95); + /// assert_eq!(gate.thresholds.min_psnr_db, 30.0); + /// ``` + pub fn settlement_dev_nominal() -> ParityGate { + ParityGate { + fixture: scenes::SETTLEMENT_DEV.clone(), + thresholds: Thresholds::INRIA_PLY_NOMINAL, + } + } + + /// Tight gate for the synthetic unit scene. + /// + /// Uses [`Thresholds::SYNTHETIC_TIGHT`] (SSIM ≥ 0.99, PSNR ≥ 40 dB). + /// + /// # Examples + /// + /// ``` + /// use cesium::fixtures::gates; + /// use cesium::oracle::ParityMetrics; + /// + /// // Identical framebuffers → SSIM 1.0, infinite PSNR → passes the tight gate. + /// let fb: Vec = vec![0.2, 0.5, 0.8, 0.1, 0.4, 0.9, 0.3, 0.6, 0.7]; + /// let metrics = ParityMetrics::compute(&fb, &fb).unwrap(); + /// + /// let gate = gates::synthetic_unit_tight(); + /// assert_eq!(gate.fixture.id, "synthetic_unit"); + /// assert_eq!(gate.thresholds.min_ssim, 0.99); + /// assert_eq!(gate.thresholds.min_psnr_db, 40.0); + /// + /// let result = gate.check(&metrics); + /// assert!(result.passed); + /// ``` + pub fn synthetic_unit_tight() -> ParityGate { + ParityGate { + fixture: scenes::SYNTHETIC_UNIT.clone(), + thresholds: Thresholds::SYNTHETIC_TIGHT, + } + } +} + +// ───────────────────────────────────────────────────────────────────────────── +// Commented scaffold — live oracle integration blocked on ndarray dep +// ───────────────────────────────────────────────────────────────────────────── + +// DEFERRED: When ndarray dep is re-enabled, wire the full fixture runner: +// +// use crate::oracle::{load_ply_batch, render_batch, ParityMetrics, OracleError}; +// use ndarray_hpc::hpc::splat3d::project::Camera; +// +// /// Run the `settlement_dev` smoke test fixture end-to-end. +// /// +// /// Loads `ply_path`, renders at `SETTLEMENT_DEV` resolution, diffs against +// /// `reference_fb`, and evaluates against `settlement_dev_smoke()`. +// /// +// /// # Sort-order risk +// /// +// /// The `SMOKE_TEST` threshold (SSIM ≥ 0.80, PSNR ≥ 20 dB) is intentionally +// /// loose to tolerate front-to-back vs back-to-front divergence. See `oracle` +// /// module doc for the named risk. +// pub fn run_settlement_dev_smoke( +// ply_path: &std::path::Path, +// reference_fb: &[f32], +// ) -> Result { +// let fixture = &scenes::SETTLEMENT_DEV; +// // UNVERIFIED: Camera::identity_at_origin signature — confirm in project.rs +// let camera = Camera::identity_at_origin(fixture.width, fixture.height); +// let metrics = crate::oracle::run_ply_oracle( +// ply_path, +// reference_fb, +// &camera, +// fixture.background, +// )?; +// Ok(gates::settlement_dev_smoke().check(&metrics)) +// } + +// ───────────────────────────────────────────────────────────────────────────── +// Unit tests (live — std-only) +// ───────────────────────────────────────────────────────────────────────────── + +#[cfg(test)] +mod tests { + use super::*; + use crate::oracle::ParityMetrics; + + fn make_metrics(ssim: f32, psnr: f32) -> ParityMetrics { + ParityMetrics { + ssim, + psnr, + mse: 0.0, + sample_count: 48, + } + } + + #[test] + fn thresholds_constants_are_ordered() { + // SMOKE < INRIA < SYNTHETIC + assert!(Thresholds::SMOKE_TEST.min_ssim < Thresholds::INRIA_PLY_NOMINAL.min_ssim); + assert!(Thresholds::INRIA_PLY_NOMINAL.min_ssim < Thresholds::SYNTHETIC_TIGHT.min_ssim); + assert!(Thresholds::SMOKE_TEST.min_psnr_db < Thresholds::INRIA_PLY_NOMINAL.min_psnr_db); + assert!(Thresholds::INRIA_PLY_NOMINAL.min_psnr_db < Thresholds::SYNTHETIC_TIGHT.min_psnr_db); + } + + #[test] + fn gate_check_passes_above_threshold() { + let gate = gates::synthetic_unit_tight(); + let m = make_metrics(0.995, 45.0); + let r = gate.check(&m); + assert!(r.passed, "metrics above threshold must pass"); + assert!(r.ssim_margin > 0.0, "ssim_margin must be positive"); + assert!(r.psnr_margin_db > 0.0, "psnr_margin_db must be positive"); + } + + #[test] + fn gate_check_fails_below_ssim() { + let gate = gates::synthetic_unit_tight(); + let m = make_metrics(0.97, 45.0); // SSIM below 0.99 + let r = gate.check(&m); + assert!(!r.passed, "SSIM below threshold must fail"); + assert!(r.ssim_margin < 0.0, "ssim_margin must be negative"); + } + + #[test] + fn gate_check_fails_below_psnr() { + let gate = gates::synthetic_unit_tight(); + let m = make_metrics(0.995, 35.0); // PSNR below 40 dB + let r = gate.check(&m); + assert!(!r.passed, "PSNR below threshold must fail"); + assert!(r.psnr_margin_db < 0.0, "psnr_margin_db must be negative"); + } + + #[test] + fn gate_check_infinite_psnr_passes() { + let gate = gates::settlement_dev_nominal(); + let m = make_metrics(0.97, f32::INFINITY); + let r = gate.check(&m); + assert!(r.passed, "infinite PSNR (identical buffers) must pass any gate"); + assert!(r.psnr_margin_db.is_infinite(), "psnr_margin_db must be infinite"); + } + + #[test] + fn smoke_gate_is_looser_than_nominal() { + let smoke = gates::settlement_dev_smoke(); + let nominal = gates::settlement_dev_nominal(); + assert!(smoke.thresholds.min_ssim <= nominal.thresholds.min_ssim); + assert!(smoke.thresholds.min_psnr_db <= nominal.thresholds.min_psnr_db); + } + + #[test] + fn fixture_ids_are_unique() { + let fixtures = [scenes::INRIA_GARDEN_SMALL.id, scenes::SETTLEMENT_DEV.id, scenes::SYNTHETIC_UNIT.id]; + let mut seen = std::collections::HashSet::new(); + for id in &fixtures { + assert!(seen.insert(*id), "duplicate fixture id: {id}"); + } + } + + #[test] + fn synthetic_unit_fixture_has_known_gaussian_count() { + assert_eq!(scenes::SYNTHETIC_UNIT.expected_gaussian_count, Some(4)); + } +} diff --git a/crates/cesium/src/hlod.rs b/crates/cesium/src/hlod.rs new file mode 100644 index 00000000..5a5bc662 --- /dev/null +++ b/crates/cesium/src/hlod.rs @@ -0,0 +1,665 @@ +//! `hlod` (group C) — HLOD (Hierarchical LOD) tree traversal with refine +//! `ADD` vs `REPLACE` semantics, gated by SSE. +//! +//! Implements the OGC 3D Tiles 1.1 tile refinement logic used by Cesium to +//! walk the tile tree and decide which tiles to render at any camera position. +//! +//! # Refinement modes (grounded against OGC 3D Tiles 1.1 §7.5) +//! +//! | Mode | Meaning | +//! |-----------|---------------------------------------------------------------| +//! | `ADD` | Child tiles are *added to* the parent's geometry. The parent | +//! | | is always rendered; children supplement it when refined. | +//! | `REPLACE` | Child tiles *replace* the parent. When refined, the parent | +//! | | is NOT rendered; children take over completely. | +//! +//! # Traversal algorithm (depth-first, greedy) +//! +//! ```text +//! fn traverse(tile, camera, frustum, max_sse, result): +//! sse = sse_for_tile(tile.geometric_error, dist(camera, tile), frustum) +//! if sse <= max_sse OR tile is leaf: +//! // This tile's LOD is sufficient — add to render list. +//! result.push(tile.content_uri) +//! return +//! // SSE too large — must refine. +//! match tile.refine: +//! REPLACE: +//! // Do NOT render the parent — recurse into children only. +//! for child in tile.children: traverse(child, ...) +//! ADD: +//! // Render the parent AND recurse into children. +//! result.push(tile.content_uri) +//! for child in tile.children: traverse(child, ...) +//! ``` +//! +//! # UNVERIFIED items +//! - Whether Cesium performs frustum culling before or after SSE testing. +//! (Likely before, but the exact pass order in `Cesium3DTileset._update` is +//! not confirmed here.) +//! - Whether `geometricError == 0` on leaf tiles is guaranteed by the spec or +//! just a convention. +//! - Exact behaviour when `tile.children` is empty but `sse > max_sse` +//! (spec says treat as leaf, but Cesium may log a warning). +//! +//! # Commented scaffold only +//! All impl is `//`-commented until reviewed live by Opus + CodeRabbit. +//! Compilable stub types and unit tests are live (std-only). + +// ───────────────────────────────────────────────────────────────────────────── +// Compilable stub types (std-only) +// ───────────────────────────────────────────────────────────────────────────── + +/// Maximum recursion depth for [`traverse_tile`]. +/// +/// 3D-Tiles tile trees are typically very shallow (< 20 levels in practice). +/// This cap prevents infinite loops on cyclic child-index graphs produced by +/// malformed or adversarial input. +const MAX_HLOD_DEPTH: usize = 64; + +/// Tile refinement strategy (OGC 3D Tiles 1.1 §7.5). +/// +/// # Examples +/// +/// ``` +/// use cesium::hlod::Refine; +/// +/// let r = Refine::Replace; +/// assert_ne!(r, Refine::Add); +/// +/// // Refine is Copy — can be used after a move. +/// let r2 = r; +/// assert_eq!(r, r2); +/// ``` +#[derive(Clone, Copy, Debug, PartialEq, Eq)] +pub enum Refine { + /// Children ADD to parent — parent is always rendered when visible. + Add, + /// Children REPLACE parent — parent is NOT rendered when refined. + Replace, +} + +/// A node in the HLOD tile tree (stub representation). +/// +/// In a live implementation this would reference the `tileset.json` content +/// and bounding-volume structures; here all geometry is reduced to the minimum +/// information needed for traversal decisions. +/// +/// # Examples +/// +/// ``` +/// use cesium::hlod::{HlodTile, Refine}; +/// +/// let tile = HlodTile { +/// geometric_error: 100.0, +/// bounding_sphere_center: [0.0, 0.0, 0.0], +/// bounding_sphere_radius: 10.0, +/// refine: Refine::Replace, +/// content_uri: Some("root.glb".into()), +/// children: vec![], +/// }; +/// assert!(tile.is_leaf()); +/// ``` +#[derive(Clone, Debug)] +pub struct HlodTile { + /// Tile's `geometricError` from `tileset.json` (world units). + pub geometric_error: f32, + /// World-space centre of the tile's bounding sphere. + pub bounding_sphere_center: [f32; 3], + /// Radius of the tile's bounding sphere (world units). + pub bounding_sphere_radius: f32, + /// Refinement strategy for this tile. + pub refine: Refine, + /// URI of this tile's content (e.g. `"0/0/0.glb"`). + /// `None` for implicit/empty tiles that have no renderable content. + pub content_uri: Option, + /// Child tile indices (into the same flat `HlodTree::tiles` Vec). + pub children: Vec, +} + +impl HlodTile { + /// True if this tile has no children (it is a leaf by structure). + /// + /// # Examples + /// + /// ``` + /// use cesium::hlod::{HlodTile, Refine}; + /// + /// let leaf = HlodTile { + /// geometric_error: 0.0, + /// bounding_sphere_center: [1.0, 2.0, 3.0], + /// bounding_sphere_radius: 5.0, + /// refine: Refine::Replace, + /// content_uri: Some("leaf.glb".into()), + /// children: vec![], + /// }; + /// assert!(leaf.is_leaf()); + /// + /// let parent = HlodTile { children: vec![1], ..leaf.clone() }; + /// assert!(!parent.is_leaf()); + /// ``` + pub fn is_leaf(&self) -> bool { + self.children.is_empty() + } + + /// Approximate distance from `camera` to the bounding-sphere surface. + /// + /// Returns 0.0 if the camera is inside the sphere. + /// + /// # Examples + /// + /// ``` + /// use cesium::hlod::{HlodTile, Refine}; + /// + /// let tile = HlodTile { + /// geometric_error: 0.0, + /// bounding_sphere_center: [0.0, 0.0, 0.0], + /// bounding_sphere_radius: 1.0, + /// refine: Refine::Replace, + /// content_uri: None, + /// children: vec![], + /// }; + /// + /// // Camera on the surface — distance is 0. + /// assert_eq!(tile.distance_to_camera([1.0, 0.0, 0.0]), 0.0); + /// + /// // Camera 5 units away from centre, radius 1 → surface distance = 4. + /// let d = tile.distance_to_camera([5.0, 0.0, 0.0]); + /// assert!((d - 4.0).abs() < 1e-5); + /// ``` + pub fn distance_to_camera(&self, camera: [f32; 3]) -> f32 { + let dx = camera[0] - self.bounding_sphere_center[0]; + let dy = camera[1] - self.bounding_sphere_center[1]; + let dz = camera[2] - self.bounding_sphere_center[2]; + let dist_center = (dx * dx + dy * dy + dz * dz).sqrt(); + (dist_center - self.bounding_sphere_radius).max(0.0) + } +} + +/// A flat HLOD tile tree (arena-allocated for cache-friendly traversal). +/// +/// # Arena layout and root invariant +/// +/// Tiles are stored in a plain `Vec` and referenced by index. **Index 0 is +/// the root by construction.** [`push`](HlodTree::push) returns the arena +/// index of each newly inserted tile; callers MUST push the root tile first +/// (before any child tiles) so that the root occupies index 0. Violating +/// this invariant causes [`traverse_hlod`] to start at the wrong tile. +/// +/// # Examples +/// +/// ``` +/// use cesium::hlod::{HlodTile, HlodTree, Refine}; +/// +/// let mut tree = HlodTree::new(); +/// +/// // Push the root first — it gets index 0. +/// let root_idx = tree.push(HlodTile { +/// geometric_error: 50.0, +/// bounding_sphere_center: [0.0, 0.0, 0.0], +/// bounding_sphere_radius: 20.0, +/// refine: Refine::Replace, +/// content_uri: Some("root.glb".into()), +/// children: vec![], // no children yet +/// }); +/// assert_eq!(root_idx, 0); +/// assert_eq!(tree.tiles.len(), 1); +/// ``` +#[derive(Clone, Debug, Default)] +pub struct HlodTree { + /// All tiles in the tree. Root is at index 0. + pub tiles: Vec, +} + +impl HlodTree { + /// Construct an empty tree. + /// + /// # Examples + /// + /// ``` + /// use cesium::hlod::HlodTree; + /// + /// let tree = HlodTree::new(); + /// assert!(tree.tiles.is_empty()); + /// ``` + pub fn new() -> Self { + Self { tiles: Vec::new() } + } + + /// Push a tile and return its index. + /// + /// **Callers must push the root tile first** (index 0) before any child + /// tiles; see the [`HlodTree`] type documentation for the root invariant. + /// + /// # Examples + /// + /// ``` + /// use cesium::hlod::{HlodTile, HlodTree, Refine}; + /// + /// let mut tree = HlodTree::new(); + /// let idx = tree.push(HlodTile { + /// geometric_error: 0.0, + /// bounding_sphere_center: [0.0, 0.0, 0.0], + /// bounding_sphere_radius: 1.0, + /// refine: Refine::Replace, + /// content_uri: Some("leaf.glb".into()), + /// children: vec![], + /// }); + /// assert_eq!(idx, 0); + /// ``` + pub fn push(&mut self, tile: HlodTile) -> usize { + let idx = self.tiles.len(); + self.tiles.push(tile); + idx + } +} + +/// Result of one HLOD traversal: the set of content URIs to render. +/// +/// # Examples +/// +/// ``` +/// use cesium::hlod::TraversalResult; +/// +/// let mut r = TraversalResult::default(); +/// assert!(r.render_list.is_empty()); +/// r.render_list.push("tile.glb".into()); +/// assert_eq!(r.render_list.len(), 1); +/// ``` +#[derive(Clone, Debug, Default)] +pub struct TraversalResult { + /// Content URIs of tiles selected for rendering (no duplicates guaranteed + /// only if the tree structure itself has no cycles). + pub render_list: Vec, +} + +// ───────────────────────────────────────────────────────────────────────────── +// Live traversal stub (calls the commented-out recursive impl via shim) +// ───────────────────────────────────────────────────────────────────────────── + +/// Traverse the HLOD tree from the root (index 0) and collect content URIs +/// to render given a camera position and SSE parameters. +/// +/// **Stub.** The real recursive body is `//`-commented below; this entry point +/// delegates to the stub shim so callers and tests compile. +/// +/// # Root invariant +/// +/// Traversal always starts at index 0, which is the root **by construction**. +/// Callers must build the tree by pushing the root tile first (index 0); see +/// [`HlodTree`] for the full invariant. If `tree.tiles` is empty this +/// function returns an empty [`TraversalResult`] immediately. +/// +/// # Defensive safety +/// +/// `traverse_tile` skips any child index that is out of range for the current +/// `tiles` slice (silent ignore — scaffold behaviour). It also caps recursion +/// at [`MAX_HLOD_DEPTH`] (64) to prevent infinite loops on cyclic child-index +/// graphs produced by malformed input. These guards do not affect well-formed +/// trees. +/// +/// # Arguments +/// - `tree` — the HLOD tile tree (root at index 0). +/// - `camera` — world-space camera position. +/// - `sse_denominator` — pre-computed value from +/// [`crate::sse::SseFrustum::sse_denominator`]. +/// - `maximum_sse` — accept threshold (Cesium default 16.0 pixels). +/// +/// # Returns +/// [`TraversalResult`] containing the render list. +/// +/// # Examples +/// +/// ``` +/// use cesium::hlod::{HlodTile, HlodTree, Refine, traverse_hlod}; +/// use cesium::sse::SseFrustum; +/// +/// // Build a tiny two-tile tree: root (Replace) → one leaf child. +/// let mut tree = HlodTree::new(); +/// +/// // Push root first (index 0). +/// let _root = tree.push(HlodTile { +/// geometric_error: 100.0, +/// bounding_sphere_center: [0.0, 0.0, 0.0], +/// bounding_sphere_radius: 10.0, +/// refine: Refine::Replace, +/// content_uri: Some("root.glb".into()), +/// children: vec![1], +/// }); +/// // Push child (index 1). +/// let _child = tree.push(HlodTile { +/// geometric_error: 0.0, +/// bounding_sphere_center: [0.0, 0.0, 0.0], +/// bounding_sphere_radius: 1.0, +/// refine: Refine::Replace, +/// content_uri: Some("child.glb".into()), +/// children: vec![], +/// }); +/// +/// let denom = SseFrustum { +/// viewport_height_px: 1080.0, +/// fovy_rad: std::f32::consts::FRAC_PI_3, +/// } +/// .sse_denominator(); +/// +/// // Camera very far away → root SSE tiny → root accepted, no child. +/// let result = traverse_hlod(&tree, [0.0, 0.0, 1_000_000.0], denom, 16.0); +/// assert!(result.render_list.contains(&"root.glb".to_string())); +/// assert!(!result.render_list.contains(&"child.glb".to_string())); +/// +/// // Camera at origin → root SSE enormous → Replace refines to child. +/// let result2 = traverse_hlod(&tree, [0.0, 0.0, 0.0], denom, 16.0); +/// assert!(!result2.render_list.contains(&"root.glb".to_string())); +/// assert!(result2.render_list.contains(&"child.glb".to_string())); +/// ``` +pub fn traverse_hlod(tree: &HlodTree, camera: [f32; 3], sse_denominator: f32, maximum_sse: f32) -> TraversalResult { + let mut result = TraversalResult::default(); + if tree.tiles.is_empty() { + return result; + } + traverse_tile(tree, 0, camera, sse_denominator, maximum_sse, &mut result, 0); + result +} + +/// Recursive inner traversal (live stub — ADD/REPLACE logic is real). +/// +/// `depth` tracks the current recursion depth. When `depth >= MAX_HLOD_DEPTH` +/// the function returns immediately, preventing infinite loops on cyclic input. +/// Any child index that is out of range for `tree.tiles` is silently skipped. +fn traverse_tile( + tree: &HlodTree, idx: usize, camera: [f32; 3], sse_denominator: f32, maximum_sse: f32, + result: &mut TraversalResult, depth: usize, +) { + // Depth guard — prevents infinite recursion on cyclic child-index graphs. + if depth >= MAX_HLOD_DEPTH { + return; + } + + let tile = &tree.tiles[idx]; + let distance = tile.distance_to_camera(camera); + let sse = crate::sse::sse_for_tile(tile.geometric_error, distance, sse_denominator); + let refine_needed = sse > maximum_sse && !tile.is_leaf(); + + match tile.refine { + Refine::Replace => { + if refine_needed { + // REPLACE + must refine → skip parent, recurse into children. + // Clone children indices to avoid borrow conflict on `tree`. + let children: Vec = tile.children.clone(); + for &child_idx in &children { + // Bounds guard — skip out-of-range indices silently. + if child_idx >= tree.tiles.len() { + continue; + } + traverse_tile(tree, child_idx, camera, sse_denominator, maximum_sse, result, depth + 1); + } + } else { + // REPLACE + SSE satisfied (or leaf) → render this tile. + if let Some(uri) = &tile.content_uri { + result.render_list.push(uri.clone()); + } + } + } + Refine::Add => { + // ADD → always render this tile (if it has content). + if let Some(uri) = &tile.content_uri { + result.render_list.push(uri.clone()); + } + if refine_needed { + // ADD + must refine → also recurse into children. + let children: Vec = tile.children.clone(); + for &child_idx in &children { + // Bounds guard — skip out-of-range indices silently. + if child_idx >= tree.tiles.len() { + continue; + } + traverse_tile(tree, child_idx, camera, sse_denominator, maximum_sse, result, depth + 1); + } + } + } + } +} + +// ───────────────────────────────────────────────────────────────────────────── +// COMMENTED-OUT full implementation notes +// ───────────────────────────────────────────────────────────────────────────── +// +// When wiring to the live tileset parser (`crate::tileset`) and the SSE module: +// +// ```rust +// // Integration sketch (not yet live): +// // +// // use crate::sse::SseFrustum; +// // use crate::tileset::Tileset; +// // +// // pub fn render_list_from_tileset( +// // tileset: &Tileset, +// // camera: [f32; 3], +// // frustum: SseFrustum, +// // maximum_sse: f32, +// // ) -> TraversalResult { +// // let tree = hlod_tree_from_tileset(tileset); // UNVERIFIED: exact mapping +// // let denom = frustum.sse_denominator(); +// // traverse_hlod(&tree, camera, denom, maximum_sse) +// // } +// // +// // Frustum-culling pass (UNVERIFIED: exact Cesium pass order): +// // Before calling traverse_tile, test whether tile bounding sphere +// // intersects the view frustum. If outside frustum → skip entirely +// // (neither render nor recurse). +// // +// // UNVERIFIED: Cesium also applies `tileset.dynamicScreenSpaceError` for +// // tiles that are both far away and have high density — a progressive +// // density-based demotion. Not implemented here. +// ``` + +// ───────────────────────────────────────────────────────────────────────────── +// Unit tests (std-only, live) +// ───────────────────────────────────────────────────────────────────────────── + +#[cfg(test)] +mod tests { + use super::*; + use crate::sse::SseFrustum; + + fn default_frustum() -> f32 { + SseFrustum { + viewport_height_px: 1080.0, + fovy_rad: std::f32::consts::FRAC_PI_3, + } + .sse_denominator() + } + + /// Build a minimal two-level tree: root → one child. + fn two_level_tree(root_refine: Refine, ge_root: f32, ge_child: f32) -> HlodTree { + let mut tree = HlodTree::new(); + let child_idx = tree.push(HlodTile { + geometric_error: ge_child, + bounding_sphere_center: [0.0, 0.0, 0.0], + bounding_sphere_radius: 1.0, + refine: Refine::Replace, + content_uri: Some("child.glb".into()), + children: vec![], + }); + tree.push(HlodTile { + geometric_error: ge_root, + bounding_sphere_center: [0.0, 0.0, 0.0], + bounding_sphere_radius: 10.0, + refine: root_refine, + content_uri: Some("root.glb".into()), + children: vec![child_idx], + }); + // Swap so root is at index 0 (last push is not root; rebuild correctly) + // Actually we must re-build: root must be index 0. + // Rebuild with correct order: + let mut t2 = HlodTree::new(); + // root first (idx 0), child will be idx 1 + let _root = t2.push(HlodTile { + geometric_error: ge_root, + bounding_sphere_center: [0.0, 0.0, 0.0], + bounding_sphere_radius: 10.0, + refine: root_refine, + content_uri: Some("root.glb".into()), + children: vec![1], // child is index 1 + }); + let _child = t2.push(HlodTile { + geometric_error: ge_child, + bounding_sphere_center: [0.0, 0.0, 0.0], + bounding_sphere_radius: 1.0, + refine: Refine::Replace, + content_uri: Some("child.glb".into()), + children: vec![], + }); + let _ = tree; // discard the mis-built one + t2 + } + + #[test] + fn empty_tree_returns_empty_list() { + let tree = HlodTree::new(); + let result = traverse_hlod(&tree, [0.0, 0.0, 100.0], default_frustum(), 16.0); + assert!(result.render_list.is_empty()); + } + + #[test] + fn single_leaf_always_rendered() { + let mut tree = HlodTree::new(); + tree.push(HlodTile { + geometric_error: 0.0, + bounding_sphere_center: [0.0, 0.0, 0.0], + bounding_sphere_radius: 1.0, + refine: Refine::Replace, + content_uri: Some("leaf.glb".into()), + children: vec![], + }); + // Camera very far away — SSE is 0 (ge=0) → accepted + let r = traverse_hlod(&tree, [0.0, 0.0, 1_000_000.0], default_frustum(), 16.0); + assert_eq!(r.render_list, vec!["leaf.glb"]); + } + + #[test] + fn replace_refine_renders_child_when_sse_too_large() { + // Root has huge geometric error close up → must refine. + // Camera at origin, root bounding sphere at origin r=1 → distance ≈ 0 → enormous SSE. + let tree = two_level_tree(Refine::Replace, 10_000.0, 0.0); + let r = traverse_hlod(&tree, [0.0, 0.0, 0.0], default_frustum(), 16.0); + // REPLACE + refine → parent NOT rendered, child rendered. + assert!(!r.render_list.contains(&"root.glb".to_string()), "root must not appear on REPLACE refine"); + assert!(r.render_list.contains(&"child.glb".to_string()), "child must appear"); + } + + #[test] + fn replace_refine_renders_root_when_sse_small() { + // Camera very far away → SSE tiny → root accepted, no refinement. + let tree = two_level_tree(Refine::Replace, 1.0, 0.0); + let r = traverse_hlod(&tree, [0.0, 0.0, 1_000_000.0], default_frustum(), 16.0); + assert!(r.render_list.contains(&"root.glb".to_string())); + assert!(!r.render_list.contains(&"child.glb".to_string()), "child must NOT appear when root accepted"); + } + + #[test] + fn add_refine_renders_both_when_sse_too_large() { + // ADD: parent always rendered; children also rendered when SSE demands it. + let tree = two_level_tree(Refine::Add, 10_000.0, 0.0); + let r = traverse_hlod(&tree, [0.0, 0.0, 0.0], default_frustum(), 16.0); + assert!(r.render_list.contains(&"root.glb".to_string()), "ADD root must always render"); + assert!(r.render_list.contains(&"child.glb".to_string()), "ADD child must render when refining"); + } + + #[test] + fn add_refine_renders_only_root_when_sse_small() { + let tree = two_level_tree(Refine::Add, 1.0, 0.0); + let r = traverse_hlod(&tree, [0.0, 0.0, 1_000_000.0], default_frustum(), 16.0); + assert!(r.render_list.contains(&"root.glb".to_string())); + // SSE satisfied → no refinement → child not included for ADD either + assert!(!r.render_list.contains(&"child.glb".to_string())); + } + + #[test] + fn tile_without_content_uri_does_not_pollute_render_list() { + let mut tree = HlodTree::new(); + // Root with no content_uri but two children with content. + tree.push(HlodTile { + geometric_error: 100.0, + bounding_sphere_center: [0.0, 0.0, 0.0], + bounding_sphere_radius: 10.0, + refine: Refine::Replace, + content_uri: None, // no content + children: vec![1, 2], + }); + tree.push(HlodTile { + geometric_error: 0.0, + bounding_sphere_center: [0.0, 0.0, 0.0], + bounding_sphere_radius: 1.0, + refine: Refine::Replace, + content_uri: Some("a.glb".into()), + children: vec![], + }); + tree.push(HlodTile { + geometric_error: 0.0, + bounding_sphere_center: [0.0, 0.0, 0.0], + bounding_sphere_radius: 1.0, + refine: Refine::Replace, + content_uri: Some("b.glb".into()), + children: vec![], + }); + // Camera close → root SSE enormous → refine → children rendered. + let r = traverse_hlod(&tree, [0.0, 0.0, 0.0], default_frustum(), 16.0); + assert!(!r.render_list.iter().any(|s| s.is_empty()), "no empty URI strings"); + assert!(r.render_list.contains(&"a.glb".to_string())); + assert!(r.render_list.contains(&"b.glb".to_string())); + assert_eq!(r.render_list.len(), 2); + } + + #[test] + fn refine_enum_semantics_documented() { + assert_ne!(Refine::Add, Refine::Replace); + // Confirm Copy semantics + let r = Refine::Add; + let r2 = r; // Copy + assert_eq!(r, r2); + } + + #[test] + fn out_of_range_child_index_is_skipped() { + // Root references a child index (99) that does not exist. + let mut tree = HlodTree::new(); + tree.push(HlodTile { + geometric_error: 10_000.0, + bounding_sphere_center: [0.0, 0.0, 0.0], + bounding_sphere_radius: 10.0, + refine: Refine::Replace, + content_uri: Some("root.glb".into()), + children: vec![99], // out of range + }); + // Must not panic; render list is empty (REPLACE + refine, but no valid child). + let r = traverse_hlod(&tree, [0.0, 0.0, 0.0], default_frustum(), 16.0); + assert!(r.render_list.is_empty()); + } + + #[test] + fn depth_guard_prevents_cycle_infinite_loop() { + // Tile 0 points to tile 1 which points back to tile 0 — a cycle. + let mut tree = HlodTree::new(); + tree.push(HlodTile { + geometric_error: 10_000.0, + bounding_sphere_center: [0.0, 0.0, 0.0], + bounding_sphere_radius: 10.0, + refine: Refine::Add, + content_uri: None, + children: vec![1], + }); + tree.push(HlodTile { + geometric_error: 10_000.0, + bounding_sphere_center: [0.0, 0.0, 0.0], + bounding_sphere_radius: 10.0, + refine: Refine::Add, + content_uri: None, + children: vec![0], // cycle back to root + }); + // Must terminate (depth guard fires at MAX_HLOD_DEPTH). + let r = traverse_hlod(&tree, [0.0, 0.0, 0.0], default_frustum(), 16.0); + // Both tiles have no content_uri — render list empty, but no panic/hang. + let _ = r; + } +} diff --git a/crates/cesium/src/implicit_tiling.rs b/crates/cesium/src/implicit_tiling.rs new file mode 100644 index 00000000..11d9b227 --- /dev/null +++ b/crates/cesium/src/implicit_tiling.rs @@ -0,0 +1,442 @@ +//! `implicit_tiling` (group A) — OGC 3D Tiles 1.1 implicit tiling: +//! subtree availability bitstreams + Morton/Z-order tile indexing. +//! +//! # Grounding +//! - OGC 3D Tiles 1.1 §9 "Implicit Tiling" (CesiumGS/3d-tiles, 22-025r4) +//! - `tile.implicitTiling.schema.json`, `availability.schema.json` (CesiumGS/3d-tiles) +//! - Subtree binary file format: magic `0x74627573` ("subt" LE), 24-byte header +//! +//! # Design contract +//! - Subtree binary files are decoded once (cold); availability bitstreams are +//! stored as compact `Vec` bit-arrays after decoding — no JSON lives in +//! any struct past the cold boundary. +//! - Morton index computation is `const fn`-compatible and branchless. +//! - All code is `//`-commented scaffold; no live implementation yet. +//! Reviewed by Opus + CodeRabbit before any impl is uncommented. +//! +//! # Availability model +//! A **subtree** covers `subtreeLevels` levels of a QUADTREE or OCTREE. +//! Within one subtree three boolean arrays express what exists: +//! +//! | Array | Source field | Always present? | +//! |---|---|---| +//! | Tile availability | `tileAvailability` | yes | +//! | Content availability | `contentAvailability` | if tile has content | +//! | Child-subtree avail | `childSubtreeAvailability` | yes | +//! +//! Each array is either a **bitstream** (index into `bufferViews`) or a +//! **constant** (0 = all absent, 1 = all present). +//! +//! # Morton indexing +//! Within each level tiles are enumerated in Morton Z-order. +//! The linear Morton index for a tile at `(x, y)` in a QUADTREE level is +//! computed by bit-interleaving `x` and `y`. For an OCTREE tile `(x, y, z)`, +//! bits of all three coordinates are interleaved. +//! +//! The **subtree-local** linear index of a tile at absolute `(level, x, y[, z])` +//! is computed as: +//! ```text +//! local_level = level mod subtree_levels +//! offset = sum_{i=0}^{local_level-1} branching_factor^i +//! morton = morton_index(x mod 2^local_level, y mod 2^local_level [, z ...]) +//! index = offset + morton +//! ``` +//! +//! Child-subtree availability uses Morton index in the leaf-below layer. + +// ───────────────────────────────────────────────────────────────────────────── +// Subtree binary header — verified against OGC 22-025r4 §9 + CesiumGS README +// ───────────────────────────────────────────────────────────────────────────── +// +// Binary layout (all fields little-endian): +// +// Offset Size Type Field +// ────── ──── ────── ────────────────────────────────────────────────────── +// 0 4 u32 magic = 0x74627573 (ASCII "subt", LE) +// 4 4 u32 version = 1 +// 8 8 u64 jsonByteLength (bytes, padded to 8-byte boundary) +// 16 8 u64 binaryByteLength (bytes, 0 if no binary chunk) +// ────── ──── ────── ────────────────────────────────────────────────────── +// 24 jsonByteLength JSON chunk (UTF-8, padded with 0x20 to 8b align) +// 24+jsonByteLength binaryByteLength Binary chunk (padded with 0x00) +// +// ```rust +// /// Magic number for subtree binary files ("subt" in ASCII, little-endian u32). +// const SUBTREE_MAGIC: u32 = 0x74627573; +// /// Only version currently defined by the spec. +// const SUBTREE_VERSION: u32 = 1; +// /// Total header size in bytes. +// const SUBTREE_HEADER_BYTES: usize = 24; +// ``` + +// ───────────────────────────────────────────────────────────────────────────── +// Planned types (all COMMENTED OUT) +// ───────────────────────────────────────────────────────────────────────────── + +// ┌─────────────────────────────────────────────────────────────────────────┐ +// │ Buffer / BufferView (mirrors subtree JSON; cold-parse only) │ +// │ │ +// │ Verified field names from subtree JSON spec: │ +// │ buffer: uri (optional string), byteLength (u64 required) │ +// │ bufferView: buffer (u32 index), byteOffset (u64), byteLength (u64) │ +// │ │ +// │ The FIRST buffer may omit `uri`, meaning it refers to the binary chunk │ +// │ embedded in the same .subtree file ("internal buffer"). │ +// └─────────────────────────────────────────────────────────────────────────┘ +// +// ```rust +// /// Decoded subtree buffer descriptor. +// pub struct SubtreeBuffer { +// /// External URI, or `None` if this is the internal binary chunk. +// pub uri: Option, +// pub byte_length: u64, +// } +// +// /// Decoded subtree buffer-view (a contiguous range inside a buffer). +// pub struct SubtreeBufferView { +// /// Index into the `buffers` array. +// pub buffer: u32, +// pub byte_offset: u64, +// pub byte_length: u64, +// } +// ``` + +// ┌─────────────────────────────────────────────────────────────────────────┐ +// │ Availability │ +// │ │ +// │ Verified field names from availability.schema.json: │ +// │ bitstream — u32 index into bufferViews (mutually excl. w/ constant)│ +// │ constant — 0 (all unavailable) or 1 (all available) │ +// │ availableCount — u64, number of 1-bits (informational, may be absent) │ +// │ │ +// │ After cold decode we materialise the bitstream into owned Vec. │ +// └─────────────────────────────────────────────────────────────────────────┘ +// +// ```rust +// /// Decoded availability array for one availability category. +// pub enum Availability { +// /// All bits unpacked from a bufferView; index `i` → bit `i`. +// /// Stored LSB-first within each byte (spec §9.5.2). +// Bitstream { +// /// Raw bytes copied from bufferView; length = ceil(total_tiles / 8). +// bits: Vec, +// /// Optional pre-counted number of available (=1) entries. +// available_count: Option, +// }, +// /// All entries share the same value (0 = absent, 1 = present). +// Constant(bool), +// } +// +// impl Availability { +// /// Test whether tile/content/child-subtree at linear index `i` is available. +// /// O(1) for both variants. +// pub fn get(&self, i: usize) -> bool { +// match self { +// Availability::Bitstream { bits, .. } => { +// (bits[i / 8] >> (i % 8)) & 1 == 1 +// } +// Availability::Constant(v) => *v, +// } +// } +// } +// ``` + +// ┌─────────────────────────────────────────────────────────────────────────┐ +// │ Decoded subtree │ +// │ │ +// │ One .subtree file spans `subtree_levels` levels. │ +// │ Verified JSON root fields (subtree JSON body after binary header): │ +// │ tileAvailability — availability object, required │ +// │ contentAvailability — array of availability objects, optional │ +// │ childSubtreeAvailability — availability object, required │ +// │ buffers — array of buffer objects, optional │ +// │ bufferViews — array of bufferView objects, optional │ +// └─────────────────────────────────────────────────────────────────────────┘ +// +// ```rust +// /// Fully decoded subtree: availability bitstreams materialised, buffers resolved. +// /// Constructed once from a .subtree binary file; no JSON values inside. +// pub struct Subtree { +// /// Which tiles exist within this subtree's level range. +// pub tile_availability: Availability, +// /// Which tiles have loadable content (one entry per content layer). +// /// Empty if the implicit root tile had no `content`/`contents`. +// pub content_availability: Vec, +// /// Which leaf-+1 cells have a child subtree file. +// pub child_subtree_availability: Availability, +// /// Number of levels this subtree covers (copy of ImplicitTilingRef::subtree_levels). +// pub subtree_levels: u32, +// /// Subdivision scheme (needed for Morton index math). +// pub scheme: crate::tileset::SubdivisionScheme, +// } +// ``` + +// ┌─────────────────────────────────────────────────────────────────────────┐ +// │ Tile coordinate │ +// └─────────────────────────────────────────────────────────────────────────┘ +// +// ```rust +// /// Absolute tile coordinate in the implicit tileset. +// /// `z` is always 0 for QUADTREE tiles. +// pub struct TileCoord { +// pub level: u32, +// pub x: u64, +// pub y: u64, +// /// z-coordinate for OCTREE only; 0 for QUADTREE. +// pub z: u64, +// } +// ``` + +// ───────────────────────────────────────────────────────────────────────────── +// Planned functions (all COMMENTED OUT) +// ───────────────────────────────────────────────────────────────────────────── + +// ┌─────────────────────────────────────────────────────────────────────────┐ +// │ Subtree binary decode │ +// └─────────────────────────────────────────────────────────────────────────┘ +// +// ```rust +// /// Decode a `.subtree` binary file into a [`Subtree`]. +// /// +// /// Validates the 24-byte header (`magic`, `version`), extracts the JSON +// /// chunk and optional binary chunk, then parses availability. +// /// +// /// The binary chunk (if present) is the "internal buffer" — buffer index 0 +// /// with no `uri`. External buffers referenced by URI are NOT fetched here; +// /// callers must supply a resolver callback. +// /// +// /// # Errors +// /// Returns [`SubtreeError`] on header mismatch, truncated data, or JSON parse +// /// failure. Unused padding bits in bitstreams (the spec requires them to be +// /// 0x00) are validated in debug builds, silently accepted in release builds. +// pub fn decode_subtree( +// bytes: &[u8], +// subtree_levels: u32, +// scheme: crate::tileset::SubdivisionScheme, +// ) -> Result { +// // 1. Check len ≥ SUBTREE_HEADER_BYTES. +// // 2. Read magic (LE u32) — must be SUBTREE_MAGIC. +// // 3. Read version (LE u32) — must be SUBTREE_VERSION. +// // 4. Read json_byte_len (LE u64) and bin_byte_len (LE u64). +// // 5. Slice JSON chunk: bytes[24 .. 24 + json_byte_len]. +// // 6. Slice binary chunk: bytes[24 + json_byte_len ..]. +// // 7. Parse JSON → SubtreeJson (internal repr, immediately discarded). +// // 8. Materialise availability bitstreams from buffer/bufferView indices. +// // 9. Return Subtree. +// todo!() +// } +// +// /// Resolve an availability JSON object into an [`Availability`] value. +// /// `buffer_views` and `binary_chunk` are passed to handle both internal +// /// and external buffers. +// fn resolve_availability( +// // json availability object fields: bitstream index or constant +// bitstream_idx: Option, +// constant: Option, +// available_count: Option, +// buffer_views: &[SubtreeBufferView], +// buffers: &[SubtreeBuffer], +// binary_chunk: &[u8], +// ) -> Result { +// // Exactly one of bitstream_idx or constant must be Some (spec invariant). +// // If bitstream: locate bufferView → locate buffer → slice bytes → copy. +// // If constant: return Availability::Constant(constant == 1). +// todo!() +// } +// ``` + +// ┌─────────────────────────────────────────────────────────────────────────┐ +// │ Morton index (Z-order curve) │ +// │ │ +// │ Standard bit-interleave algorithm. No external dep. │ +// │ For QUADTREE: interleave bits of (x, y) → 2-bit groups. │ +// │ For OCTREE: interleave bits of (x, y, z) → 3-bit groups. │ +// │ │ +// │ Reference: any standard Morton encoding — technique is well-known and │ +// │ unambiguous; no UNVERIFIED marks needed for the algorithm itself. │ +// └─────────────────────────────────────────────────────────────────────────┘ +// +// ```rust +// /// Compute the Morton Z-order index for a QUADTREE tile at `(x, y)`. +// /// +// /// Bits of `x` and `y` are interleaved: result bit 2k = x bit k, +// /// result bit 2k+1 = y bit k. Both coordinates must be < 2^32. +// /// +// /// This is the tile's position within its level's linear enumeration, +// /// as required by OGC 3D Tiles 1.1 §9. +// pub const fn morton2(x: u64, y: u64) -> u64 { +// // Spread bits of x into even positions, y into odd positions. +// // Classic "magic number" spread technique (no loop, O(log bits) shifts): +// // spread(n) via masks 0x5555_5555_5555_5555 etc. +// // let sx = spread_bits(x); +// // let sy = spread_bits(y); +// // sx | (sy << 1) +// todo!() +// } +// +// /// Compute the Morton Z-order index for an OCTREE tile at `(x, y, z)`. +// /// +// /// Bits of x, y, z are interleaved in groups of 3: +// /// result bit 3k = x bit k, 3k+1 = y bit k, 3k+2 = z bit k. +// pub const fn morton3(x: u64, y: u64, z: u64) -> u64 { +// // Spread each coordinate into positions 0,3,6,... / 1,4,7,... / 2,5,8,... +// // let sx = spread_bits3(x); +// // let sy = spread_bits3(y); +// // let sz = spread_bits3(z); +// // sx | (sy << 1) | (sz << 2) +// todo!() +// } +// +// /// Spread the 21 low bits of `n` into every-other bit (for morton2). +// /// Bits land at positions 0,2,4,...,40; upper bits of result are 0. +// const fn spread_bits2(mut n: u64) -> u64 { +// // n = (n | (n << 16)) & 0x0000_FFFF_0000_FFFF; +// // n = (n | (n << 8)) & 0x00FF_00FF_00FF_00FF; +// // n = (n | (n << 4)) & 0x0F0F_0F0F_0F0F_0F0F; +// // n = (n | (n << 2)) & 0x3333_3333_3333_3333; +// // n = (n | (n << 1)) & 0x5555_5555_5555_5555; +// // n +// todo!() +// } +// +// /// Spread the 21 low bits of `n` into every-third bit (for morton3). +// const fn spread_bits3(mut n: u64) -> u64 { +// // Standard 3D interleave masks for 21 bits. +// // n &= 0x1F_FFFF; +// // n = (n | n << 32) & 0x001F_0000_0000_FFFF; +// // n = (n | n << 16) & 0x001F_0000_FF00_00FF; +// // n = (n | n << 8) & 0x100F_00F0_0F00_F00F; +// // n = (n | n << 4) & 0x10C3_0C30_C30C_30C3; +// // n = (n | n << 2) & 0x1249_2492_4924_9249; +// // n +// todo!() +// } +// ``` + +// ┌─────────────────────────────────────────────────────────────────────────┐ +// │ Subtree-local linear index │ +// │ │ +// │ Converts absolute (level, x, y[, z]) → subtree-local linear index. │ +// │ Used to look up Availability::get(index). │ +// └─────────────────────────────────────────────────────────────────────────┘ +// +// ```rust +// /// Compute the subtree-local linear index for a tile at absolute coordinate. +// /// +// /// For QUADTREE the formula is: +// /// local_level = coord.level % subtree_levels +// /// level_offset = (4^local_level - 1) / 3 [geometric series sum] +// /// index = level_offset + morton2(x_local, y_local) +// /// where x_local = coord.x % 2^local_level, y_local = coord.y % 2^local_level +// /// +// /// For OCTREE replace 4→8, /3→/7, morton2→morton3 with z_local. +// pub fn subtree_local_index( +// coord: &TileCoord, +// subtree_levels: u32, +// scheme: &crate::tileset::SubdivisionScheme, +// ) -> usize { +// todo!() +// } +// +// /// Compute the child-subtree availability index for a tile at the leaf level +// /// of a subtree (the level immediately below the last level in the subtree). +// /// This is simply morton2/morton3 of (x % branching, y % branching[, z %]). +// pub fn child_subtree_index( +// coord: &TileCoord, +// subtree_levels: u32, +// scheme: &crate::tileset::SubdivisionScheme, +// ) -> usize { +// todo!() +// } +// ``` + +// ┌─────────────────────────────────────────────────────────────────────────┐ +// │ URI template expansion │ +// │ │ +// │ Template variable substitution for subtrees.uri patterns. │ +// │ QUADTREE: {level}, {x}, {y} │ +// │ OCTREE: {level}, {x}, {y}, {z} │ +// └─────────────────────────────────────────────────────────────────────────┘ +// +// ```rust +// /// Expand a subtrees URI template for the given tile coordinate. +// /// +// /// Template variables `{level}`, `{x}`, `{y}` (and `{z}` for octrees) are +// /// substituted with their decimal representations. The result is an +// /// ASCII-only relative URI suitable for I/O resolution. +// /// +// /// # Example (QUADTREE) +// /// template = `"subtrees/{level}/{x}/{y}.subtree"` +// /// coord = TileCoord { level: 2, x: 1, y: 3, z: 0 } +// /// result = `"subtrees/2/1/3.subtree"` +// pub fn expand_subtree_uri(template: &str, coord: &TileCoord) -> String { +// todo!() +// } +// ``` + +// ┌─────────────────────────────────────────────────────────────────────────┐ +// │ Implicit tile tree expansion │ +// │ │ +// │ Converts an ImplicitTilingRef + loaded Subtrees into concrete │ +// │ TileNode entries that the rest of the crate can traverse uniformly. │ +// │ This is the bridge back to `tileset::TileNode`. │ +// └─────────────────────────────────────────────────────────────────────────┘ +// +// ```rust +// /// Expand an implicit tiling root into a flat list of available tile coordinates. +// /// +// /// Iterates through all subtrees referenced by `implicit_ref`, decodes each +// /// .subtree file via `fetch_subtree` (caller-supplied I/O), checks tile +// /// availability, and emits [`TileCoord`]s for every available tile. +// /// +// /// Content URIs are constructed by the caller using the tile's `content.uri` +// /// template from the parent `Tileset`, which is out of scope here. +// /// +// /// # Design note +// /// This function does no I/O itself — it drives I/O through a callback so +// /// the crate stays dependency-free at this level. +// pub fn expand_implicit( +// implicit_ref: &crate::tileset::ImplicitTilingRef, +// fetch_subtree: F, +// ) -> Result, SubtreeError> +// where +// F: Fn(&str) -> Result, SubtreeError>, +// { +// // 1. Iterate root subtree coord (level=0, x=0, y=0, z=0). +// // 2. Expand subtree_uri template → call fetch_subtree. +// // 3. decode_subtree() → Subtree. +// // 4. Walk all tiles within subtree: for each available tile, emit coord. +// // 5. For each available child-subtree, recurse (BFS or DFS). +// // 6. Stop when level >= implicit_ref.available_levels. +// todo!() +// } +// ``` + +// ───────────────────────────────────────────────────────────────────────────── +// Error type (COMMENTED OUT) +// ───────────────────────────────────────────────────────────────────────────── +// +// ```rust +// pub enum SubtreeError { +// /// File too short to contain the 24-byte header. +// TruncatedHeader, +// /// Magic bytes do not match 0x74627573. +// BadMagic(u32), +// /// Version is not 1. +// UnsupportedVersion(u32), +// /// JSON chunk is not valid UTF-8. +// JsonUtf8Error, +// /// JSON syntax or missing required field. +// JsonParse(String), +// /// bufferView index out of range. +// BadBufferViewIndex(u32), +// /// buffer index out of range. +// BadBufferIndex(u32), +// /// Buffer slice out of bounds. +// BufferSliceBounds, +// /// Availability had neither `bitstream` nor `constant` (spec violation). +// MalformedAvailability, +// /// External buffer fetch callback returned an error. +// FetchError(String), +// } +// ``` diff --git a/crates/cesium/src/khr_gs.rs b/crates/cesium/src/khr_gs.rs new file mode 100644 index 00000000..ecd6966e --- /dev/null +++ b/crates/cesium/src/khr_gs.rs @@ -0,0 +1,472 @@ +//! `khr_gs` (group A) — glTF `KHR_gaussian_splatting` parse +//! (attributes POSITION, ROTATION quat, SCALE vec3, OPACITY scalar, +//! SH_DEGREE_n_COEF_i coefficients; COLOR_0 fallback; cold import → intermediate structs) +//! +//! # Grounding +//! - Khronos `KHR_gaussian_splatting` extension README +//! (KhronosGroup/glTF, extensions/2.0/Khronos/KHR_gaussian_splatting) +//! - glTF 2.0 core spec (accessor componentType / type encoding, primitive modes) +//! +//! # Design contract +//! - **Cold import only.** glTF JSON is parsed once; no serde/JSON values +//! survive into [`GsSplats`]. Accessor binary data is read once and stored +//! as typed intermediate arrays. +//! - **Intermediate structs → `to_cam_soa`.** [`GsSplats`] maps 1-to-1 to +//! `splat3d::gaussian::GaussianBatch` fields. The actual conversion is owned +//! by [`crate::to_cam_soa`]; this module only parses. +//! - **No live implementation.** All planned code is `//`-commented scaffold; +//! reviewed by Opus + CodeRabbit before any impl is uncommented. +//! +//! # Extension location (verified) +//! The extension object lives at: +//! ```text +//! glTF → meshes[i] → primitives[j] → extensions → KHR_gaussian_splatting +//! ``` +//! The primitive `mode` MUST be `0` (POINTS). +//! Splat attributes live in `primitive.attributes`, NOT inside the extension object. +//! +//! # Extension object fields (verified) +//! | Field | Type | Required | Default | +//! |---|---|---|---| +//! | `kernel` | string | yes | — | +//! | `colorSpace` | string | yes | — | +//! | `projection` | string | no | `"perspective"` | +//! | `sortingMethod`| string | no | `"cameraDistance"`| +//! +//! `kernel` is currently always `"ellipse"`. +//! `colorSpace` is one of `"srgb_rec709_display"` or `"lin_rec709_display"`. +//! +//! # Attribute semantics (verified) +//! All attributes are at `primitive.attributes`; the extension-namespaced ones +//! use the prefix `KHR_gaussian_splatting:`. +//! +//! | Attribute semantic | Accessor type | Component types | Required | +//! |---|---|---|---| +//! | `POSITION` | VEC3 | FLOAT | yes | +//! | `KHR_gaussian_splatting:ROTATION` | VEC4 | FLOAT, BYTE_NORM, SHORT_NORM | yes | +//! | `KHR_gaussian_splatting:SCALE` | VEC3 | FLOAT, UBYTE, UBYTE_NORM, USHORT, USHORT_NORM | yes | +//! | `KHR_gaussian_splatting:OPACITY` | SCALAR | FLOAT, UBYTE_NORM, USHORT_NORM | yes | +//! | `KHR_gaussian_splatting:SH_DEGREE_0_COEF_0` | VEC3 | FLOAT | yes (if SH used) | +//! | `KHR_gaussian_splatting:SH_DEGREE_1_COEF_[0-2]` | VEC3 | FLOAT | optional | +//! | `KHR_gaussian_splatting:SH_DEGREE_2_COEF_[0-4]` | VEC3 | FLOAT | optional | +//! | `KHR_gaussian_splatting:SH_DEGREE_3_COEF_[0-6]` | VEC3 | FLOAT | optional | +//! | `COLOR_0` | VEC4 | per glTF core spec | optional fallback | +//! +//! ROTATION quaternion component order: (x, y, z, w) — standard glTF VEC4 order. +//! Quantized ROTATION: signed byte/short normalized to [-1, 1]. +//! Quantized SCALE: unsigned byte/short (normalized or raw — see below). +//! OPACITY: linear value in [0.0, 1.0]. +//! COLOR_0: fallback diffuse ≈ clamp(SH_0_0 * 0.282095 + 0.5, 0, 1). +//! +//! SH coefficients must be provided in ascending degree order; +//! lower degrees must be fully populated before higher ones are used. +//! +//! # UNVERIFIED items +//! - Whether quantized SCALE uses UBYTE (un-normalized) or UBYTE_NORM — the +//! search result listed both; confirm against the normative schema. +//! - Exact glTF componentType integer codes for BYTE_NORM, SHORT_NORM, etc. +//! (glTF core: BYTE=5120, UBYTE=5121, SHORT=5122, USHORT=5123, FLOAT=5126). + +// ───────────────────────────────────────────────────────────────────────────── +// glTF accessor constants — verified against glTF 2.0 core spec §3.6.2 +// ───────────────────────────────────────────────────────────────────────────── +// +// ```rust +// /// glTF accessor component type codes. +// mod gltf_component_type { +// pub const BYTE: u32 = 5120; +// pub const UNSIGNED_BYTE: u32 = 5121; +// pub const SHORT: u32 = 5122; +// pub const UNSIGNED_SHORT: u32 = 5123; +// pub const UNSIGNED_INT: u32 = 5125; +// pub const FLOAT: u32 = 5126; +// } +// +// /// glTF accessor type strings (from `accessor.type`). +// mod gltf_accessor_type { +// pub const SCALAR: &str = "SCALAR"; +// pub const VEC2: &str = "VEC2"; +// pub const VEC3: &str = "VEC3"; +// pub const VEC4: &str = "VEC4"; +// pub const MAT2: &str = "MAT2"; +// pub const MAT3: &str = "MAT3"; +// pub const MAT4: &str = "MAT4"; +// } +// +// /// glTF primitive mode — POINTS is the only valid mode for KHR_gaussian_splatting. +// const PRIMITIVE_MODE_POINTS: u32 = 0; +// ``` + +// ───────────────────────────────────────────────────────────────────────────── +// Attribute name constants (verified against KHR_gaussian_splatting README) +// ───────────────────────────────────────────────────────────────────────────── +// +// ```rust +// /// Attribute semantic names as they appear in `primitive.attributes`. +// mod attr { +// pub const POSITION: &str = "POSITION"; +// pub const COLOR_0: &str = "COLOR_0"; +// pub const ROTATION: &str = "KHR_gaussian_splatting:ROTATION"; +// pub const SCALE: &str = "KHR_gaussian_splatting:SCALE"; +// pub const OPACITY: &str = "KHR_gaussian_splatting:OPACITY"; +// // Degree-0 (DC) spherical harmonic — required when SH is used. +// pub const SH_0_0: &str = "KHR_gaussian_splatting:SH_DEGREE_0_COEF_0"; +// // Degree-1 (3 coefficients: indices 0-2). +// pub const SH_1: [&str; 3] = [ +// "KHR_gaussian_splatting:SH_DEGREE_1_COEF_0", +// "KHR_gaussian_splatting:SH_DEGREE_1_COEF_1", +// "KHR_gaussian_splatting:SH_DEGREE_1_COEF_2", +// ]; +// // Degree-2 (5 coefficients: indices 0-4). +// pub const SH_2: [&str; 5] = [ +// "KHR_gaussian_splatting:SH_DEGREE_2_COEF_0", +// "KHR_gaussian_splatting:SH_DEGREE_2_COEF_1", +// "KHR_gaussian_splatting:SH_DEGREE_2_COEF_2", +// "KHR_gaussian_splatting:SH_DEGREE_2_COEF_3", +// "KHR_gaussian_splatting:SH_DEGREE_2_COEF_4", +// ]; +// // Degree-3 (7 coefficients: indices 0-6). +// pub const SH_3: [&str; 7] = [ +// "KHR_gaussian_splatting:SH_DEGREE_3_COEF_0", +// "KHR_gaussian_splatting:SH_DEGREE_3_COEF_1", +// "KHR_gaussian_splatting:SH_DEGREE_3_COEF_2", +// "KHR_gaussian_splatting:SH_DEGREE_3_COEF_3", +// "KHR_gaussian_splatting:SH_DEGREE_3_COEF_4", +// "KHR_gaussian_splatting:SH_DEGREE_3_COEF_5", +// "KHR_gaussian_splatting:SH_DEGREE_3_COEF_6", +// ]; +// } +// ``` + +// ───────────────────────────────────────────────────────────────────────────── +// Planned types (all COMMENTED OUT) +// ───────────────────────────────────────────────────────────────────────────── + +// ┌─────────────────────────────────────────────────────────────────────────┐ +// │ Parsed KHR_gaussian_splatting extension object │ +// │ │ +// │ The extension object at mesh.primitive.extensions.KHR_gaussian_splatting│ +// │ contains metadata only — attributes are parsed separately. │ +// └─────────────────────────────────────────────────────────────────────────┘ +// +// ```rust +// /// Decoded `KHR_gaussian_splatting` extension descriptor (metadata only). +// pub struct KhrGsExtension { +// /// The Gaussian kernel shape. Currently only `"ellipse"` is defined. +// pub kernel: KernelType, +// /// Color space of reconstructed splat color. +// pub color_space: ColorSpace, +// /// Projection type. Default: `"perspective"`. +// pub projection: ProjectionType, +// /// Sorting method for alpha compositing. Default: `"cameraDistance"`. +// pub sorting_method: SortingMethod, +// } +// +// /// Kernel types defined by the extension (extensible in future revisions). +// pub enum KernelType { +// /// 2D ellipse projecting an ellipsoid (the only currently defined value). +// Ellipse, +// /// UNVERIFIED: forward-compat catch-all for unknown kernel strings. +// Unknown(String), +// } +// +// /// Color space of the reconstructed splat color output. +// pub enum ColorSpace { +// /// `"srgb_rec709_display"` — BT.709 sRGB display-referred. +// SrgbRec709Display, +// /// `"lin_rec709_display"` — BT.709 linear display-referred. +// LinRec709Display, +// /// UNVERIFIED: forward-compat for additional color spaces. +// Unknown(String), +// } +// +// /// Projection type for Gaussian rasterisation. +// pub enum ProjectionType { +// Perspective, +// /// UNVERIFIED: whether orthographic is defined in the extension. +// Unknown(String), +// } +// +// /// Sorting method for depth-ordered alpha compositing. +// pub enum SortingMethod { +// CameraDistance, +// /// UNVERIFIED: additional sorting methods may be defined. +// Unknown(String), +// } +// ``` + +// ┌─────────────────────────────────────────────────────────────────────────┐ +// │ Intermediate parsed splat batch │ +// │ │ +// │ All attribute accessors decoded into flat f32 arrays (SoA orientation). │ +// │ Layout mirrors splat3d::gaussian::GaussianBatch so to_cam_soa can │ +// │ copy/transform with minimal indirection. │ +// │ │ +// │ SH layout: coefficients packed contiguously per splat: │ +// │ [sh_deg0 (3 f32)] [sh_deg1_0..2 (9 f32)] [sh_deg2_0..4 (15 f32)] │ +// │ [sh_deg3_0..6 (21 f32)] (only present up to the degree loaded) │ +// │ Total SH floats per splat: 3 + 9 + 15 + 21 = 48 (degree-3 full). │ +// └─────────────────────────────────────────────────────────────────────────┘ +// +// ```rust +// /// Intermediate representation of all KHR_gaussian_splatting splats from +// /// one glTF mesh primitive. Owned; no references into original bytes. +// /// +// /// All decoded values are in f32. Quantized source data is dequantised +// /// during decode (not deferred to to_cam_soa). +// /// +// /// Field lengths: every non-empty Vec has length `splat_count` (for +// /// per-splat scalars) or `splat_count * N` (for per-splat VecN). +// pub struct GsSplats { +// /// Number of Gaussian splats in this primitive. +// pub splat_count: usize, +// /// Extension metadata (kernel, color space, etc.). +// pub extension: KhrGsExtension, +// /// Splat centre positions in local space. Layout: [x0,y0,z0, x1,y1,z1, ...]. +// /// Source attribute: `POSITION` (VEC3 FLOAT, dequant n/a). +// pub position: Vec, // len = splat_count * 3 +// /// Unit quaternions (x,y,z,w). Layout: [x0,y0,z0,w0, x1,...]. +// /// Source attribute: `KHR_gaussian_splatting:ROTATION` (VEC4). +// /// Quantized byte/short inputs are dequantized to f32 here. +// pub rotation_xyzw: Vec, // len = splat_count * 4 +// /// Scale along each principal axis. Layout: [sx0,sy0,sz0, sx1,...]. +// /// Source attribute: `KHR_gaussian_splatting:SCALE` (VEC3). +// /// UNVERIFIED: whether raw UBYTE means literal scale or needs a decode table. +// pub scale: Vec, // len = splat_count * 3 +// /// Linear opacity in [0.0, 1.0]. +// /// Source attribute: `KHR_gaussian_splatting:OPACITY` (SCALAR). +// pub opacity: Vec, // len = splat_count +// /// Spherical harmonics coefficients (all degrees, contiguous per splat). +// /// Absent (empty) if the primitive carries no SH attributes. +// /// Source attributes: SH_DEGREE_n_COEF_i (VEC3 FLOAT). +// pub sh_coefficients: Vec, // len = splat_count * (3 * (degree+1)^2), or 0 +// /// Highest SH degree loaded (0-3); 0 = only DC term. +// /// None = no SH attributes present (use color_0_rgba fallback). +// pub sh_degree: Option, +// /// Optional fallback diffuse RGBA (from `COLOR_0`). +// /// Present only when COLOR_0 attribute exists in the primitive. +// /// Layout: [r0,g0,b0,a0, r1,...] in [0.0, 1.0]. +// pub color_0_rgba: Option>, // len = splat_count * 4, or None +// } +// ``` + +// ┌─────────────────────────────────────────────────────────────────────────┐ +// │ Quantization dequantise helpers │ +// │ │ +// │ These are called during cold decode to produce uniform f32 arrays. │ +// └─────────────────────────────────────────────────────────────────────────┘ +// +// ```rust +// /// Dequantise a signed byte normalized value to f32 in [-1.0, 1.0]. +// /// glTF core spec §3.6.2.3: max(value / 127.0, -1.0). +// #[inline] +// fn dequant_i8_norm(v: i8) -> f32 { +// ((v as f32) / 127.0_f32).max(-1.0_f32) +// } +// +// /// Dequantise a signed short normalized value to f32 in [-1.0, 1.0]. +// /// glTF core spec §3.6.2.3: max(value / 32767.0, -1.0). +// #[inline] +// fn dequant_i16_norm(v: i16) -> f32 { +// ((v as f32) / 32767.0_f32).max(-1.0_f32) +// } +// +// /// Dequantise an unsigned byte normalized value to f32 in [0.0, 1.0]. +// /// glTF core spec §3.6.2.3: value / 255.0. +// #[inline] +// fn dequant_u8_norm(v: u8) -> f32 { +// (v as f32) / 255.0_f32 +// } +// +// /// Dequantise an unsigned short normalized value to f32 in [0.0, 1.0]. +// /// glTF core spec §3.6.2.3: value / 65535.0. +// #[inline] +// fn dequant_u16_norm(v: u16) -> f32 { +// (v as f32) / 65535.0_f32 +// } +// ``` + +// ───────────────────────────────────────────────────────────────────────────── +// Planned functions (all COMMENTED OUT) +// ───────────────────────────────────────────────────────────────────────────── + +// ┌─────────────────────────────────────────────────────────────────────────┐ +// │ Top-level parse entry point │ +// └─────────────────────────────────────────────────────────────────────────┘ +// +// ```rust +// /// Parse all KHR_gaussian_splatting splat primitives from a glTF binary +// /// (GLB) or JSON+BIN pair. +// /// +// /// Returns one [`GsSplats`] per qualifying mesh primitive (mode=POINTS, +// /// extension present). A file may contain zero, one, or many. +// /// +// /// This is the ONLY function that reads JSON and binary accessor data. +// /// After it returns, no JSON/serde values survive. +// /// +// /// # Arguments +// /// - `gltf_json`: raw bytes of the glTF JSON (the JSON chunk of a GLB, or +// /// a standalone `.gltf` file). +// /// - `bin_chunk`: the BIN chunk from a GLB, or an external `.bin` buffer. +// /// Pass `None` if the file has no binary data. +// /// +// /// # Errors +// /// Returns [`KhrGsError`] on JSON parse failure, missing required attributes, +// /// unsupported accessor types, or buffer bounds violations. +// pub fn parse_khr_gaussian_splatting( +// gltf_json: &[u8], +// bin_chunk: Option<&[u8]>, +// ) -> Result, KhrGsError> { +// // 1. Decode glTF JSON → minimal internal repr (no serde dep yet): +// // - meshes[*].primitives[*] where mode == 0 (POINTS) +// // - filter primitives that have extensions.KHR_gaussian_splatting +// // 2. For each qualifying primitive, call parse_gs_primitive(). +// // 3. Return Vec. +// todo!() +// } +// +// /// Parse one qualifying mesh primitive into [`GsSplats`]. +// fn parse_gs_primitive( +// // primitive attributes map (semantic → accessor index) +// // extension object fields +// // accessors array +// // bufferViews array +// // buffers array +// // bin_chunk +// ) -> Result { +// // 1. Parse extension object → KhrGsExtension. +// // 2. Verify mode == PRIMITIVE_MODE_POINTS. +// // 3. Extract POSITION accessor → decode_vec3_f32() → position vec. +// // 4. Extract ROTATION accessor → decode_rotation() (handles quantized types). +// // 5. Extract SCALE accessor → decode_scale() (handles quantized types). +// // 6. Extract OPACITY accessor → decode_opacity() (handles quantized types). +// // 7. Determine SH degree: probe SH_DEGREE_3_* first, step down. +// // 8. Extract SH accessors for each present degree → decode_sh_vec3_f32(). +// // 9. Optionally extract COLOR_0 accessor. +// // 10. Validate all lengths == splat_count. +// // 11. Return GsSplats. +// todo!() +// } +// ``` + +// ┌─────────────────────────────────────────────────────────────────────────┐ +// │ Accessor decode helpers │ +// │ │ +// │ Each helper resolves accessor → bufferView → buffer slice and returns │ +// │ a decoded Vec. All dequantisation happens here (cold path). │ +// └─────────────────────────────────────────────────────────────────────────┘ +// +// ```rust +// /// Decode a VEC3 FLOAT accessor into a flat Vec. +// /// Used for POSITION and SH coefficients. +// fn decode_vec3_f32( +// // accessor index, accessors, buffer_views, buffers, bin_chunk +// ) -> Result, KhrGsError> { +// // Validate componentType == FLOAT, type == "VEC3". +// // Slice buffer bytes, cast to &[f32] (checking alignment). +// // Copy to owned Vec (respect byteStride if nonzero). +// todo!() +// } +// +// /// Decode ROTATION VEC4 accessor; dequantise BYTE_NORM or SHORT_NORM if needed. +// /// Output: flat Vec with [x,y,z,w] per splat, all normalised. +// fn decode_rotation( +// // accessor index + context +// ) -> Result, KhrGsError> { +// // Match componentType: +// // FLOAT (5126): direct copy (4 f32 per element). +// // BYTE (5120) normalized: dequant_i8_norm for each of 4 components. +// // SHORT (5122) normalized: dequant_i16_norm for each of 4 components. +// // Other: KhrGsError::UnsupportedComponentType. +// todo!() +// } +// +// /// Decode SCALE VEC3 accessor; dequantise if quantized. +// /// UNVERIFIED: exact dequant formula for UBYTE (non-normalized) variant. +// fn decode_scale( +// // accessor index + context +// ) -> Result, KhrGsError> { +// // Match componentType: +// // FLOAT (5126): direct copy. +// // UNSIGNED_BYTE (5121) normalized: dequant_u8_norm. +// // UNSIGNED_BYTE (5121) not normalized: UNVERIFIED — cast as f32 directly? +// // UNSIGNED_SHORT (5123) normalized: dequant_u16_norm. +// // UNSIGNED_SHORT (5123) not normalized: UNVERIFIED. +// // Other: KhrGsError::UnsupportedComponentType. +// todo!() +// } +// +// /// Decode OPACITY SCALAR accessor; dequantise normalized types. +// /// Output: flat Vec in [0.0, 1.0]. +// fn decode_opacity( +// // accessor index + context +// ) -> Result, KhrGsError> { +// // Match componentType: +// // FLOAT (5126): direct copy; clamp to [0,1] in debug builds. +// // UNSIGNED_BYTE (5121) normalized: dequant_u8_norm. +// // UNSIGNED_SHORT (5123) normalized: dequant_u16_norm. +// // Other: KhrGsError::UnsupportedComponentType. +// todo!() +// } +// +// /// Detect the highest SH degree present by probing attribute names. +// /// Returns 0 if only SH_DEGREE_0_COEF_0 is present, 3 if all degrees present. +// /// Returns None if SH_DEGREE_0_COEF_0 is absent (no SH data at all). +// fn detect_sh_degree( +// // attribute name set +// ) -> Option { +// // Probe SH_3[6] present → return Some(3). +// // Probe SH_2[4] present → return Some(2). +// // Probe SH_1[2] present → return Some(1). +// // Probe SH_0_0 present → return Some(0). +// // None. +// todo!() +// } +// +// /// Decode all SH coefficient accessors (all degrees up to `max_degree`) +// /// into a single contiguous Vec per splat. +// /// +// /// Layout per splat: [DC(3), deg1_coef0(3),..,deg1_coef2(3), +// /// deg2_coef0(3),..,deg2_coef4(3), +// /// deg3_coef0(3),..,deg3_coef6(3)] +// /// +// /// All coefficients are FLOAT (VEC3); no dequantisation needed. +// fn decode_sh_coefficients( +// // attribute map, max_degree, accessor context +// ) -> Result, KhrGsError> { +// todo!() +// } +// ``` + +// ───────────────────────────────────────────────────────────────────────────── +// Error type (COMMENTED OUT) +// ───────────────────────────────────────────────────────────────────────────── +// +// ```rust +// pub enum KhrGsError { +// /// glTF JSON is invalid UTF-8. +// JsonUtf8, +// /// glTF JSON syntax error. +// JsonSyntax(String), +// /// Required primitive attribute is missing. +// MissingAttribute(&'static str), +// /// Accessor componentType is not supported for this attribute. +// UnsupportedComponentType { attribute: &'static str, component_type: u32 }, +// /// Accessor type string mismatch (e.g. expected VEC3, got SCALAR). +// AccessorTypeMismatch { attribute: &'static str, expected: &'static str, got: String }, +// /// Primitive mode is not POINTS (0). +// NotPointsPrimitive(u32), +// /// Buffer or bufferView index out of range. +// BufferIndexOutOfRange { index: u32, max: usize }, +// /// Byte slice is out of bounds for the given accessor offset/count/stride. +// AccessorBoundsError, +// /// SH coefficient accessor count mismatch (lower degrees absent). +// ShDegreeGap(u8), +// /// Splat count inconsistency across attribute accessors. +// SplatCountMismatch { attribute: &'static str, expected: usize, got: usize }, +// /// Extension object missing required field. +// MissingExtensionField(&'static str), +// /// UNVERIFIED: extension field had an unrecognised string value (forward-compat warn). +// UnrecognizedExtensionValue { field: &'static str, value: String }, +// } +// ``` diff --git a/crates/cesium/src/lib.rs b/crates/cesium/src/lib.rs new file mode 100644 index 00000000..bd3e649b --- /dev/null +++ b/crates/cesium/src/lib.rs @@ -0,0 +1,58 @@ +#![allow(dead_code)] +//! # cesium — optional Cesium / ArcGIS reference + parity oracle +//! +//! **NOT a production dependency.** This crate is the *runnable legacy pillar* +//! (see `.claude/knowledge/crate-registry.md`, modesty clause): it +//! reverse-engineers external geospatial-splat formats into ndarray CAM SoA and +//! diffs our `splat3d` render against a reference, to keep "parity / better" +//! honest by **measurement**. It is a workspace member but **NOT** in +//! `default-members`, so it is optional/dev by construction. It relocates to +//! `lance-graph/crates/cesium` once flawless. +//! +//! ## Hard rules (from crate-registry.md) +//! - **Reverse-engineer ONLY** — extract data, rebuild as CAM SoA; never depend +//! on, round-trip, or emit the source format. Refactor > reverse-engineer. +//! - **CAM SoA target** — output is ndarray's splat SoA +//! (`splat3d::gaussian::GaussianBatch`) + `cam_pq` codebook indices. No AoS, +//! no source-native structs survive past the import boundary. +//! - **NO JSON IN THE HOTPATH, EVER** — JSON only at the cold import boundary. +//! Prefer the **ArcGIS PBF (protobuf, binary)** path over `f=json`. +//! +//! ## Grounding references (reverse-engineer against these — do NOT fabricate; +//! mark anything you cannot confirm against a source as `// UNVERIFIED:`) +//! - `CesiumGS/cesium` — 3D Tiles / HLOD / KHR_gaussian_splatting behaviour +//! - Khronos `KHR_gaussian_splatting` + `KHR_gaussian_splatting_compression_spz` +//! - `R-ArcGIS/arcpbf`, `AdaWorldAPI/arcgisutils`, `dfridkin/arcGIS-Rust-runtime`, +//! `EsriDevEvents/locup` — ArcGIS FeatureService **PBF** (binary, no-JSON) + `ESRI_crs` +//! - `staehlli/3D-settlement-development` — a real 3DGS scene for fixtures +//! - Inria `.ply` via ndarray `splat3d::ply` — first reference render for the oracle +//! +//! ## Module map (12 slots — `///` commented scaffold ONLY until reviewed live +//! by Opus + CodeRabbit/Codex) +//! | module | group | purpose | +//! |---|---|---| +//! | [`tileset`] | A | 3D Tiles `tileset.json` / `.3tz` parse (cold import) | +//! | [`implicit_tiling`] | A | subtree availability bitstreams + Morton/Z-order | +//! | [`khr_gs`] | A | KHR_gaussian_splatting glTF parse (POSITION/ROTATION/SCALE/COLOR_0/SH_DEGREE_n) | +//! | [`arcgis_pbf`] | B | ArcGIS FeatureService **PBF** ingestion (binary, the no-JSON path) | +//! | [`spz`] | B | `KHR_gaussian_splatting_compression_spz` decode (quantize+gzip) | +//! | [`esri_crs`] | B | `ESRI_crs` WKID (h+v) reproject → WGS84 / local-origin | +//! | [`to_cam_soa`] | C | reverse-engineer parsed splats → CAM SoA (`GaussianBatch` + `cam_pq`) | +//! | [`sse`] | C | screen-space-error LOD selection (`geometricError·viewportH/(dist·2tan(fovy/2))`) | +//! | [`hlod`] | C | HLOD tree + refine ADD·REPLACE | +//! | [`point_fallback`] | D | point-cloud fallback (KHR base extension = point primitives) | +//! | [`oracle`] | D | reference-render invocation + SSIM/PSNR diff vs `splat3d` | +//! | [`fixtures`] | D | golden scenes + parity gates | + +pub mod tileset; +pub mod implicit_tiling; +pub mod khr_gs; +pub mod arcgis_pbf; +pub mod spz; +pub mod esri_crs; +pub mod to_cam_soa; +pub mod sse; +pub mod hlod; +pub mod point_fallback; +pub mod oracle; +pub mod fixtures; diff --git a/crates/cesium/src/oracle.rs b/crates/cesium/src/oracle.rs new file mode 100644 index 00000000..568a38cd --- /dev/null +++ b/crates/cesium/src/oracle.rs @@ -0,0 +1,517 @@ +//! `oracle` (group D) — parity harness: load a reference, render via `splat3d`, +//! diff (SSIM / PSNR). +//! +//! # Purpose +//! +//! This is a **measurement instrument**, not an assertion engine. It quantifies +//! how close our `splat3d` render is to a reference image — it never asserts +//! "must be above N dB" at runtime. Threshold gates live in [`crate::fixtures`]. +//! +//! # Reference pipeline (first reference = Inria `.ply`) +//! +//! ```text +//! Inria .ply file +//! │ read_ply::>() ← splat3d::ply::read_ply +//! ▼ +//! GaussianBatch (SoA) ← splat3d::gaussian::GaussianBatch +//! │ (project → tile → rasterize) ← splat3d::raster::rasterize_frame +//! ▼ +//! framebuffer: Vec (interleaved RGB, 3·W·H) +//! │ ParityMetrics::compute(&ref_fb, &our_fb) +//! ▼ +//! ParityMetrics { ssim, psnr } +//! ``` +//! +//! # Named risk: sort-order divergence +//! +//! `splat3d::tile::TileBinning::from_projected` sorts instances by packed u64 +//! key `(tile_id << 32) | depth_bits` using `sort_unstable_by_key`, giving +//! **front-to-back** depth order within each tile. The rasterizer in +//! `raster::rasterize_tile` then performs a **front-to-back alpha-blend** with +//! early-out at `T < T_SATURATION_EPS = 1e-4`. +//! +//! The reference renderers (brush, Inria `diff-gaussian-rasterization`) may +//! sort **back-to-front** (painter's-algorithm) instead. When two or more +//! Gaussians have nearly equal depth, floating-point accumulation order is +//! reversed between the two approaches, producing per-pixel colour differences +//! of order `O(alpha²)`. For high-opacity scenes this can be several dB of +//! PSNR. Any SSIM/PSNR delta unexplained by numerical noise should be +//! investigated for sort-order as the primary suspect. +//! +//! # Reference renderer — same binary, NO FFI +//! +//! There is **no Cesium/cesium-native FFI** and there never will be. The +//! oracle's "reference" is one of the existing in-scope 3DGS renderers wired +//! as an ordinary workspace dependency in the **same binary** (A-to-Z, no +//! `extern "C"`, no browser/Node, no C++ ABI): +//! +//! - **`brush`** (`AdaWorldAPI/brush`, Rust + `wgpu`) — runnable in-binary on +//! CPU/any-GPU with no CUDA toolchain; the default reference path. +//! - **Inria `gaussian-splatting`** + **`diff-gaussian-rasterization`** +//! (`AdaWorldAPI/*`, CUDA) — the ground-truth rasterizer when a CUDA device +//! is present; same binary, linked as a normal Rust crate, not via FFI. +//! +//! All three speak our CAM SoA (`splat3d::gaussian::GaussianBatch`) at the +//! boundary, so the oracle compares like-for-like: load → render(ours) vs +//! render(reference) → SSIM/PSNR. None of these is a foreign binary; the +//! wiring is `[dependencies]` + a direct call, gated by `cfg`/feature, once +//! `ndarray` is re-enabled in `Cargo.toml`. +//! +//! # Implementation status +//! +//! **Commented scaffold only.** `ndarray` is commented out in `Cargo.toml`; +//! the stub types (`ParityMetrics`, `OracleError`) are `std`-only and compile +//! clean. All impl blocks that touch `splat3d` symbols are commented out. +//! +//! # Source grounding +//! +//! Symbols cited from the local codebase (confirmed by reading before writing): +//! - `read_ply` — `src/hpc/splat3d/ply.rs` line 110: `pub fn read_ply(reader: R) -> Result` +//! - `rasterize_frame` — `src/hpc/splat3d/raster.rs` line 213: `pub fn rasterize_frame(binning: &TileBinning, projected: &ProjectedBatch, framebuffer: &mut [f32], width: u32, height: u32, background: [f32; 3])` +//! - `TileBinning::from_projected` — `src/hpc/splat3d/tile.rs` line 105: sort by `(tile_id << 32) | depth_bits`, front-to-back +//! - `GaussianBatch` — `src/hpc/splat3d/gaussian.rs`: SoA fields `mean_x`, `mean_y`, `mean_z`, `len` +//! - Framebuffer layout — `raster.rs` doc: interleaved RGB `[R0,G0,B0,R1,…]`, len = `3·width·height` + +// ───────────────────────────────────────────────────────────────────────────── +// Stub types (live — std-only, no ndarray dep) +// ───────────────────────────────────────────────────────────────────────────── + +/// Errors the oracle can return. +/// +/// # Examples +/// +/// ``` +/// use cesium::oracle::OracleError; +/// +/// let e = OracleError::DimensionMismatch { ref_len: 12, candidate_len: 9 }; +/// assert!(format!("{e}").contains("12")); +/// assert!(format!("{e}").contains("9")); +/// +/// let e2 = OracleError::EmptyFramebuffer; +/// assert!(format!("{e2}").contains("empty")); +/// ``` +#[derive(Debug)] +pub enum OracleError { + /// I/O error reading a reference file. + Io(std::io::Error), + /// PLY parse failed. Payload is a human-readable description. + PlyParse(String), + /// Framebuffer dimensions mismatch between reference and candidate. + DimensionMismatch { ref_len: usize, candidate_len: usize }, + /// The render pipeline returned an empty framebuffer (width or height = 0). + EmptyFramebuffer, + // DEFERRED: ReferenceRender(String) — reserved for failures from the + // in-binary reference renderer (brush / Inria), wired as a workspace + // dependency once `ndarray` is re-enabled in Cargo.toml. NOT an FFI path. + // ReferenceRender(String), +} + +impl std::fmt::Display for OracleError { + fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { + match self { + OracleError::Io(e) => write!(f, "oracle I/O error: {e}"), + OracleError::PlyParse(msg) => write!(f, "oracle PLY parse error: {msg}"), + OracleError::DimensionMismatch { ref_len, candidate_len } => { + write!(f, "oracle dimension mismatch: ref={ref_len} candidate={candidate_len}") + } + OracleError::EmptyFramebuffer => write!(f, "oracle: empty framebuffer (width or height = 0)"), + } + } +} + +impl From for OracleError { + fn from(e: std::io::Error) -> Self { + OracleError::Io(e) + } +} + +/// SSIM / PSNR parity metrics between a reference and a candidate framebuffer. +/// +/// Both metrics are computed over the full `3·W·H` interleaved RGB buffer +/// (not per-channel). Channel weighting and luminance-only variants are +/// `// DEFERRED:` until the fixture suite validates a specific convention. +/// +/// # SSIM +/// +/// Structural Similarity Index Measure in `[−1, 1]`; 1.0 = identical. +/// A coarse approximation (global mean/variance, not 11×11 Gaussian window) +/// is acceptable for the parity gate; the fixture thresholds in +/// [`crate::fixtures`] are set conservatively to hide the approximation error. +/// +/// # PSNR +/// +/// Peak Signal-to-Noise Ratio in dB: `10·log10(MAX² / MSE)` where `MAX=1.0` +/// (framebuffer in `[0,1]`). `f32::INFINITY` when the buffers are identical +/// (MSE = 0). Practically ≥ 30 dB = visually acceptable; ≥ 40 dB = excellent. +/// +/// # Sort-order risk +/// +/// See module-level doc. Any SSIM < 0.98 or PSNR < 35 dB for a scene with +/// depth-ordered overlapping Gaussians should be attributed to sort-order +/// divergence before tuning thresholds. +/// +/// # Examples +/// +/// ``` +/// use cesium::oracle::ParityMetrics; +/// +/// // Construct metrics directly to inspect fields. +/// let m = ParityMetrics { ssim: 0.97, psnr: 38.5, mse: 0.00014, sample_count: 48 }; +/// assert!(m.ssim > 0.95); +/// assert!(m.psnr.is_finite()); +/// ``` +#[derive(Debug, Clone, PartialEq)] +pub struct ParityMetrics { + /// SSIM in `[−1, 1]`. 1.0 = identical. + pub ssim: f32, + /// PSNR in dB. `f32::INFINITY` when buffers are identical. + pub psnr: f32, + /// Mean squared error (raw, before log). + pub mse: f32, + /// Number of f32 samples compared (`3 · width · height`). + pub sample_count: usize, +} + +impl ParityMetrics { + /// Compute SSIM and PSNR between two equal-length interleaved RGB + /// framebuffers. Both must have the same `len = 3 · width · height`. + /// + /// Returns `Err(OracleError::DimensionMismatch)` if lengths differ, or + /// `Err(OracleError::EmptyFramebuffer)` if either is empty. + /// + /// # Algorithm (coarse approximation — acceptable for parity gating) + /// + /// ```text + /// mu_r = mean(ref), mu_c = mean(candidate) + /// var_r = var(ref), var_c = var(candidate) + /// cov = covariance(ref, candidate) + /// C1 = (0.01)², C2 = (0.03)² [standard luminance/contrast constants] + /// SSIM = (2·mu_r·mu_c + C1)(2·cov + C2) + /// / ((mu_r² + mu_c² + C1)(var_r + var_c + C2)) + /// MSE = mean((ref[i] - cand[i])²) + /// PSNR = 10·log10(1.0 / MSE) (MAX=1.0) + /// ``` + /// + /// # Examples + /// + /// ``` + /// fn main() -> Result<(), cesium::oracle::OracleError> { + /// use cesium::oracle::ParityMetrics; + /// + /// // Identical buffers → PSNR = ∞, SSIM = 1.0. + /// let fb = vec![0.5f32, 0.3, 0.8, 0.1, 0.9, 0.2]; + /// let m = ParityMetrics::compute(&fb, &fb)?; + /// assert!(m.psnr.is_infinite()); + /// assert!((m.ssim - 1.0).abs() < 1e-5); + /// + /// // Non-equal buffers → finite PSNR. + /// let ref_fb = vec![1.0f32; 6]; + /// let cand_fb = vec![0.0f32; 6]; + /// let m2 = ParityMetrics::compute(&ref_fb, &cand_fb)?; + /// assert!(m2.psnr.is_finite()); + /// + /// Ok(()) + /// } + /// ``` + pub fn compute(reference: &[f32], candidate: &[f32]) -> Result { + let n = reference.len(); + if n == 0 || candidate.is_empty() { + return Err(OracleError::EmptyFramebuffer); + } + if n != candidate.len() { + return Err(OracleError::DimensionMismatch { + ref_len: n, + candidate_len: candidate.len(), + }); + } + + let nf = n as f64; + + // Means + let mu_r: f64 = reference.iter().map(|&v| v as f64).sum::() / nf; + let mu_c: f64 = candidate.iter().map(|&v| v as f64).sum::() / nf; + + // Variances and covariance (single-pass after mean) + let mut var_r = 0.0f64; + let mut var_c = 0.0f64; + let mut cov = 0.0f64; + let mut mse = 0.0f64; + for (&r, &c) in reference.iter().zip(candidate.iter()) { + let dr = r as f64 - mu_r; + let dc = c as f64 - mu_c; + var_r += dr * dr; + var_c += dc * dc; + cov += dr * dc; + let diff = r as f64 - c as f64; + mse += diff * diff; + } + var_r /= nf; + var_c /= nf; + cov /= nf; + mse /= nf; + + // SSIM constants (standard Wang et al. 2004 values for [0,1] range) + const C1: f64 = 0.01 * 0.01; + const C2: f64 = 0.03 * 0.03; + let num = (2.0 * mu_r * mu_c + C1) * (2.0 * cov + C2); + let den = (mu_r * mu_r + mu_c * mu_c + C1) * (var_r + var_c + C2); + let ssim = if den.abs() < 1e-30 { 1.0f32 } else { (num / den) as f32 }; + + let psnr = if mse < 1e-30 { + f32::INFINITY + } else { + (10.0 * (1.0 / mse).log10()) as f32 + }; + + Ok(ParityMetrics { + ssim, + psnr, + mse: mse as f32, + sample_count: n, + }) + } + + /// Returns `true` if this result passes the given SSIM and PSNR thresholds. + /// + /// The oracle measures parity — it does not assert it. Callers in + /// [`crate::fixtures`] decide whether to gate on these results. + /// + /// # Examples + /// + /// ``` + /// fn main() -> Result<(), cesium::oracle::OracleError> { + /// use cesium::oracle::ParityMetrics; + /// + /// // Identical buffers pass any reasonable threshold. + /// let fb = vec![0.5f32, 0.3, 0.8, 0.1, 0.9, 0.2]; + /// let m = ParityMetrics::compute(&fb, &fb)?; + /// assert!(m.passes(0.99, 40.0)); + /// + /// // Significantly different buffers fail a tight threshold. + /// let m2 = ParityMetrics::compute(&vec![1.0f32; 6], &vec![0.0f32; 6])?; + /// assert!(!m2.passes(0.99, 40.0)); + /// + /// Ok(()) + /// } + /// ``` + pub fn passes(&self, min_ssim: f32, min_psnr_db: f32) -> bool { + self.ssim >= min_ssim && (self.psnr >= min_psnr_db || self.psnr == f32::INFINITY) + } +} + +// ───────────────────────────────────────────────────────────────────────────── +// Commented scaffold — live impl blocked on ndarray dep +// ───────────────────────────────────────────────────────────────────────────── + +// DEFERRED: wire the full pipeline once `ndarray = { workspace = true }` is +// re-enabled and the module passes Opus + CodeRabbit review. +// +// use std::fs::File; +// use std::io::BufReader; +// +// // ── splat3d symbols (confirmed from local source, not fabricated) ──────────── +// // read_ply — ply.rs:110 pub fn read_ply(reader: R) -> Result +// // rasterize_frame — raster.rs:213 pub fn rasterize_frame(binning: &TileBinning, projected: &ProjectedBatch, framebuffer: &mut [f32], width: u32, height: u32, background: [f32; 3]) +// // TileBinning::from_projected — tile.rs:105 sorts front-to-back by (tile_id<<32)|depth_bits +// // GaussianBatch — gaussian.rs: SoA with len, mean_x/y/z, scale_x/y/z, quat_w/x/y/z, opacity, sh +// // Framebuffer: interleaved RGB [R0,G0,B0,R1,…], len = 3·width·height (raster.rs doc) +// // +// use ndarray_hpc::hpc::splat3d::ply::{read_ply, PlyError}; +// use ndarray_hpc::hpc::splat3d::gaussian::GaussianBatch; +// use ndarray_hpc::hpc::splat3d::raster::rasterize_frame; +// use ndarray_hpc::hpc::splat3d::tile::TileBinning; +// use ndarray_hpc::hpc::splat3d::project::{Camera, ProjectedBatch, project_batch}; +// +// impl From for OracleError { +// fn from(e: PlyError) -> Self { +// OracleError::PlyParse(format!("{e:?}")) +// } +// } +// +// /// Load an Inria `.ply` file and return a `GaussianBatch`. +// /// +// /// Entry point for the first parity reference. Applies the canonical +// /// activation transforms inline (sigmoid(opacity), exp(scale), normalize(quat)) +// /// as documented in `ply.rs`. +// /// +// /// # Sort-order risk (see module doc) +// /// +// /// The batch returned here is unsorted. `TileBinning::from_projected` +// /// will sort instances front-to-back via the packed u64 key. If the +// /// reference render was produced with back-to-front order, the PSNR +// /// comparison will reflect the sort-order delta. +// pub fn load_ply_batch(path: &std::path::Path) -> Result { +// let file = File::open(path)?; +// let reader = BufReader::new(file); +// let batch = read_ply(reader)?; // ply.rs::read_ply +// Ok(batch) +// } +// +// /// Render a `GaussianBatch` via the `splat3d` pipeline and return the +// /// framebuffer as `Vec` (interleaved RGB, len = `3·width·height`). +// /// +// /// # Named risk: sort-order +// /// +// /// `rasterize_frame` (raster.rs:213) consumes instances in the order +// /// produced by `TileBinning::from_projected` (tile.rs:105), which is +// /// **front-to-back** (ascending depth). If the reference was rendered +// /// back-to-front, every pixel with overlapping Gaussians will differ. +// /// +// /// # UNVERIFIED: camera parameters +// /// +// /// The `Camera` struct (project.rs) takes `(width, height, fovy_rad, +// /// near, far)`. The identity camera `Camera::identity_at_origin` is used +// /// here as a placeholder; real parity comparisons need matched intrinsics. +// // UNVERIFIED: Camera::identity_at_origin exists in project.rs — confirm +// // when wiring this function. +// pub fn render_batch( +// batch: &GaussianBatch, +// camera: &Camera, +// background: [f32; 3], +// ) -> Result, OracleError> { +// let w = camera.width; +// let h = camera.height; +// if w == 0 || h == 0 { +// return Err(OracleError::EmptyFramebuffer); +// } +// let mut projected = ProjectedBatch::with_capacity(batch.len); +// project_batch(batch, camera, &mut projected); // UNVERIFIED: project_batch fn name +// let binning = TileBinning::from_projected(&projected, camera); // tile.rs:105 +// let mut framebuffer = vec![0.0f32; 3 * w as usize * h as usize]; +// rasterize_frame(&binning, &projected, &mut framebuffer, w, h, background); // raster.rs:213 +// Ok(framebuffer) +// } +// +// /// Full oracle run: load `.ply`, render, compare to a reference framebuffer. +// /// +// /// Returns `ParityMetrics`; caller decides whether to gate on the result. +// pub fn run_ply_oracle( +// ply_path: &std::path::Path, +// reference_fb: &[f32], +// camera: &Camera, +// background: [f32; 3], +// ) -> Result { +// let batch = load_ply_batch(ply_path)?; +// let candidate_fb = render_batch(&batch, camera, background)?; +// ParityMetrics::compute(reference_fb, &candidate_fb) +// } + +// ───────────────────────────────────────────────────────────────────────────── +// DEFERRED: in-binary reference renderer (brush / Inria) — NO FFI +// ───────────────────────────────────────────────────────────────────────────── + +// The cross-renderer reference is an existing in-scope 3DGS renderer wired as +// an ordinary workspace dependency in the SAME binary — no `extern "C"`, no +// `#[link]`, no browser/Node, no foreign ABI. It consumes the same CAM SoA +// `GaussianBatch` we feed `splat3d`, renders to the same interleaved-RGB +// framebuffer layout, and we diff the two with `ParityMetrics::compute`. +// +// DEFERRED until `ndarray` (and the chosen reference crate) are re-enabled in +// Cargo.toml and the module passes Opus + CodeRabbit review: +// +// // Default reference: brush (Rust + wgpu, runnable without CUDA). +// // UNVERIFIED: brush's public render entry point + camera/SoA adapter shape; +// // confirm against AdaWorldAPI/brush before wiring. +// fn render_reference_brush( +// batch: &GaussianBatch, +// camera: &Camera, +// background: [f32; 3], +// ) -> Result, OracleError> { +// // brush::render(...) → Vec interleaved RGB, len = 3·W·H +// todo!("wire brush as same-binary workspace dep") +// } +// +// // Ground truth (when a CUDA device is present): Inria gaussian-splatting / +// // diff-gaussian-rasterization, linked as a normal Rust crate (NOT FFI), +// // gated behind a `cuda` cfg/feature. +// // UNVERIFIED: the Rust-facing entry point exposed by the Inria crates. +// #[cfg(feature = "cuda")] +// fn render_reference_inria( +// batch: &GaussianBatch, +// camera: &Camera, +// background: [f32; 3], +// ) -> Result, OracleError> { +// todo!("wire diff-gaussian-rasterization as same-binary workspace dep") +// } + +// ───────────────────────────────────────────────────────────────────────────── +// Unit tests (live — std-only) +// ───────────────────────────────────────────────────────────────────────────── + +#[cfg(test)] +mod tests { + use super::*; + + fn make_fb(values: &[f32]) -> Vec { + values.to_vec() + } + + #[test] + fn parity_metrics_identical_buffers_gives_inf_psnr_and_ssim_one() { + let fb = make_fb(&[0.5f32, 0.3, 0.8, 0.1, 0.9, 0.2]); + let m = ParityMetrics::compute(&fb, &fb).unwrap(); + assert!(m.psnr.is_infinite(), "identical buffers must give infinite PSNR"); + assert!((m.ssim - 1.0).abs() < 1e-5, "identical buffers must give SSIM = 1.0, got {}", m.ssim); + assert!(m.mse < 1e-10, "identical buffers MSE must be 0, got {}", m.mse); + } + + #[test] + fn parity_metrics_zero_vs_one_gives_finite_psnr() { + let ref_fb = vec![1.0f32; 48]; + let cand_fb = vec![0.0f32; 48]; + let m = ParityMetrics::compute(&ref_fb, &cand_fb).unwrap(); + assert!(m.psnr.is_finite(), "PSNR must be finite for non-identical buffers, got {}", m.psnr); + // MSE = 1.0, PSNR = 10*log10(1/1) = 0 dB + assert!((m.psnr - 0.0).abs() < 0.01, "PSNR for all-zero vs all-one should be 0 dB, got {}", m.psnr); + assert!((m.mse - 1.0).abs() < 1e-5, "MSE must be 1.0, got {}", m.mse); + } + + #[test] + fn parity_metrics_passes_respects_thresholds() { + let ref_fb = vec![0.5f32; 12]; + let cand_fb = ref_fb.clone(); + let m = ParityMetrics::compute(&ref_fb, &cand_fb).unwrap(); + // Identical → passes any threshold + assert!(m.passes(0.99, 40.0), "identical buffers must pass any threshold"); + // Fabricate a bad metrics to test failure path + let bad = ParityMetrics { + ssim: 0.5, + psnr: 20.0, + mse: 0.01, + sample_count: 12, + }; + assert!(!bad.passes(0.99, 40.0), "bad metrics must fail tight threshold"); + assert!(bad.passes(0.4, 15.0), "bad metrics must pass loose threshold"); + } + + #[test] + fn parity_metrics_dimension_mismatch_returns_error() { + let a = vec![0.0f32; 12]; + let b = vec![0.0f32; 9]; + match ParityMetrics::compute(&a, &b) { + Err(OracleError::DimensionMismatch { + ref_len: 12, + candidate_len: 9, + }) => {} + other => panic!("expected DimensionMismatch, got {other:?}"), + } + } + + #[test] + fn parity_metrics_empty_returns_error() { + match ParityMetrics::compute(&[], &[]) { + Err(OracleError::EmptyFramebuffer) => {} + other => panic!("expected EmptyFramebuffer, got {other:?}"), + } + } + + #[test] + fn oracle_error_display_is_human_readable() { + let e = OracleError::DimensionMismatch { + ref_len: 100, + candidate_len: 99, + }; + let s = format!("{e}"); + assert!(s.contains("100") && s.contains("99"), "display must include lens: {s}"); + } +} diff --git a/crates/cesium/src/point_fallback.rs b/crates/cesium/src/point_fallback.rs new file mode 100644 index 00000000..0e67053f --- /dev/null +++ b/crates/cesium/src/point_fallback.rs @@ -0,0 +1,273 @@ +//! `point_fallback` (group D) — point-cloud fallback renderer. +//! +//! When a [`GaussianBatch`] cannot be rendered as full 3DGS splats (e.g. the +//! receiver only supports the **KHR_gaussian_splatting base extension**, which +//! mandates point-primitive output for coarse-tier HHTL preview), this module +//! extracts the mean positions `(x, y, z)` from the SoA buffers and emits a +//! flat `[f32]` point cloud suitable for upload as a KHR point-primitive mesh. +//! +//! # KHR_gaussian_splatting base extension +//! +//! The Khronos base extension (as of 2024-Q3 draft) requires the runtime to +//! fall back to rendering splats as point primitives when the implementation +//! does not support the full covariance / SH pipeline. The HHTL coarse-tier +//! preview uses exactly this path: positions only, no colour or opacity. +//! +//! # What this module does NOT do +//! +//! - Project positions onto screen (that's `splat3d::project`). +//! - Emit colour / opacity / SH — point primitives carry `POSITION` only in +//! the base extension. A colour channel could be synthesised from the DC +//! SH term (`sh[0]` per channel) but is left `// DEFERRED:` until the +//! extension spec stabilises. +//! +//! # Implementation status +//! +//! **Commented scaffold only.** No live impl that touches ndarray exists here; +//! `ndarray` is commented out in `Cargo.toml`. The stub structs and unit +//! tests below are `std`-only and compile clean. +//! +//! # Source grounding +//! +//! Symbols cited from the local codebase (read before writing, not fabricated): +//! - `GaussianBatch` — `src/hpc/splat3d/gaussian.rs` (fields `mean_x`, +//! `mean_y`, `mean_z`, `len`, `capacity`). +//! - `read_ply` — `src/hpc/splat3d/ply.rs` (loader entry point that produces +//! `GaussianBatch` from an Inria binary-PLY reader). + +// ───────────────────────────────────────────────────────────────────────────── +// Stub types (live — std-only, no ndarray dep) +// ───────────────────────────────────────────────────────────────────────────── + +/// Output of the point-cloud fallback extraction. +/// +/// Flat interleaved XYZ: `[x0, y0, z0, x1, y1, z1, …]`, length `3 * count`. +/// Suitable for upload as a KHR_gaussian_splatting point-primitive `POSITION` +/// accessor (component type `FLOAT`, count = `point_count`). +/// +/// # Invariant +/// +/// The safe constructor [`PointCloud::new`] guarantees `xyz.len() == 3 * +/// point_count`. The public fields are kept for backwards-compatible literal +/// construction, but callers are strongly encouraged to use `new` to avoid +/// silent desync between `xyz` and `point_count`. +/// +/// # Examples +/// +/// ``` +/// let pc = cesium::point_fallback::PointCloud::new(vec![1.0, 2.0, 3.0, 4.0, 5.0, 6.0]); +/// assert_eq!(pc.point_count, 2); +/// assert_eq!(pc.position(0), [1.0, 2.0, 3.0]); +/// assert!(pc.try_position(5).is_none()); +/// ``` +#[derive(Debug, Clone, PartialEq)] +pub struct PointCloud { + /// Flat XYZ buffer: length `3 * point_count`. + pub xyz: Vec, + /// Number of points (= `GaussianBatch::len` of the source batch). + pub point_count: usize, +} + +impl PointCloud { + /// Construct a [`PointCloud`] from a flat XYZ buffer, deriving + /// `point_count` as `xyz.len() / 3`. + /// + /// This is the preferred constructor because it guarantees the internal + /// invariant `xyz.len() == 3 * point_count` by construction. If `xyz` + /// has 1 or 2 trailing floats that do not form a complete triplet they are + /// ignored by integer division (e.g. 7 floats → `point_count = 2`, the + /// 7th float is kept in the buffer but unreachable through `position` / + /// `try_position`). + /// + /// # Examples + /// + /// ``` + /// let pc = cesium::point_fallback::PointCloud::new(vec![1.0, 2.0, 3.0, 4.0, 5.0, 6.0]); + /// assert_eq!(pc.point_count, 2); + /// assert_eq!(pc.xyz.len(), 6); + /// + /// // Stray trailing float is ignored by integer division. + /// let pc2 = cesium::point_fallback::PointCloud::new(vec![1.0, 2.0, 3.0, 4.0]); + /// assert_eq!(pc2.point_count, 1); + /// ``` + pub fn new(xyz: Vec) -> Self { + let point_count = xyz.len() / 3; + Self { xyz, point_count } + } + + /// Return the position of point `i` as `[x, y, z]`. + /// + /// # Panics + /// + /// Panics if `i >= point_count` or if `xyz` has been manually desynced + /// from `point_count` so that the required floats are missing. Use + /// [`PointCloud::try_position`] for a non-panicking alternative, and + /// [`PointCloud::new`] to construct values that cannot desync. + /// + /// # Examples + /// + /// ``` + /// let pc = cesium::point_fallback::PointCloud::new(vec![1.0, 2.0, 3.0, 4.0, 5.0, 6.0]); + /// assert_eq!(pc.position(0), [1.0, 2.0, 3.0]); + /// assert_eq!(pc.position(1), [4.0, 5.0, 6.0]); + /// ``` + /// + /// ```should_panic + /// let pc = cesium::point_fallback::PointCloud::new(vec![1.0, 2.0, 3.0]); + /// let _ = pc.position(1); // index out of range — panics + /// ``` + pub fn position(&self, i: usize) -> [f32; 3] { + self.try_position(i) + .expect("PointCloud::position: index out of range or xyz/point_count desynced") + } + + /// Return the position of point `i` as `Some([x, y, z])`, or `None` if + /// `i` is out of range or the struct has been manually desynced so that + /// the required floats are absent from `xyz`. + /// + /// This is the non-panicking counterpart to [`PointCloud::position`]. + /// + /// # Examples + /// + /// ``` + /// let pc = cesium::point_fallback::PointCloud::new(vec![1.0, 2.0, 3.0, 4.0, 5.0, 6.0]); + /// assert_eq!(pc.try_position(0), Some([1.0, 2.0, 3.0])); + /// assert_eq!(pc.try_position(1), Some([4.0, 5.0, 6.0])); + /// assert!(pc.try_position(2).is_none()); + /// assert!(pc.try_position(5).is_none()); + /// ``` + pub fn try_position(&self, i: usize) -> Option<[f32; 3]> { + if i >= self.point_count { + return None; + } + let base = i * 3; + if base + 2 >= self.xyz.len() { + return None; + } + Some([self.xyz[base], self.xyz[base + 1], self.xyz[base + 2]]) + } + + /// Byte length of the XYZ buffer (for KHR accessor `byteLength`). + /// + /// # Examples + /// + /// ``` + /// let pc = cesium::point_fallback::PointCloud::new(vec![0.0, 0.0, 0.0, 1.0, 1.0, 1.0]); + /// // 2 points × 3 floats × 4 bytes = 24 + /// assert_eq!(pc.byte_length(), 24); + /// ``` + pub fn byte_length(&self) -> usize { + self.xyz.len() * 4 + } +} + +// ───────────────────────────────────────────────────────────────────────────── +// Commented scaffold — live impl blocked on ndarray dep +// ───────────────────────────────────────────────────────────────────────────── + +// DEFERRED: uncomment when `ndarray = { workspace = true }` is re-enabled in +// Cargo.toml and the module goes live post Opus + CodeRabbit review. +// +// use ndarray_hpc::hpc::splat3d::gaussian::GaussianBatch; +// +// /// Extract mean positions from a [`GaussianBatch`] into a [`PointCloud`]. +// /// +// /// Reads `batch.mean_x`, `batch.mean_y`, `batch.mean_z` (SoA, len = +// /// `batch.len`) and interleaves them into a flat `[x, y, z, …]` buffer. +// /// +// /// # Complexity +// /// +// /// O(n) time, O(3n) extra allocation. No SIMD — this is a cold-path +// /// extraction, not a hot-path accumulation. +// /// +// /// # KHR accessor layout +// /// +// /// The output xyz buffer maps directly to a glTF `POSITION` accessor: +// /// ```text +// /// accessors[k] = { +// /// bufferView: , +// /// componentType: 5126, // FLOAT +// /// type: "VEC3", +// /// count: point_count, +// /// } +// /// ``` +// pub fn extract_point_cloud(batch: &GaussianBatch) -> PointCloud { +// let n = batch.len; // GaussianBatch::len — active gaussian count +// let mut xyz = Vec::with_capacity(3 * n); +// for i in 0..n { +// // SoA fields: batch.mean_x[i], batch.mean_y[i], batch.mean_z[i] +// // (confirmed from gaussian.rs — pub Vec per axis) +// xyz.push(batch.mean_x[i]); +// xyz.push(batch.mean_y[i]); +// xyz.push(batch.mean_z[i]); +// } +// PointCloud { xyz, point_count: n } +// } +// +// // DEFERRED: synthesise colour from DC SH term once KHR base-ext colour +// // semantics stabilise. +// // +// // The DC SH coefficient for channel R is `batch.sh[i * 48 + 0]`, G at +// // `i * 48 + 16`, B at `i * 48 + 32` (SH_COEFFS_PER_CHANNEL = 16 from +// // gaussian.rs). The Inria colour offset (+0.5) must be applied before +// // clamping to [0, 1]: +// // color_r = (SH_C0 * sh_r + 0.5).clamp(0.0, 1.0) +// // where SH_C0 ≈ 0.282095 (from sh.rs::SH_C0). +// // +// // UNVERIFIED: KHR point-primitive extension colour accessor name — may be +// // COLOR_0 (glTF standard) or a custom attribute. + +// ───────────────────────────────────────────────────────────────────────────── +// Unit tests (live — std-only) +// ───────────────────────────────────────────────────────────────────────────── + +#[cfg(test)] +mod tests { + use super::*; + + fn make_point_cloud(xyzs: &[[f32; 3]]) -> PointCloud { + let mut xyz = Vec::with_capacity(xyzs.len() * 3); + for p in xyzs { + xyz.extend_from_slice(p); + } + PointCloud { + xyz, + point_count: xyzs.len(), + } + } + + #[test] + fn point_cloud_position_round_trips() { + let pts = [[1.0f32, 2.0, 3.0], [4.0, 5.0, 6.0], [-1.0, 0.5, 0.0]]; + let pc = make_point_cloud(&pts); + assert_eq!(pc.point_count, 3); + for (i, expected) in pts.iter().enumerate() { + let got = pc.position(i); + assert_eq!(got, *expected, "position({i}) mismatch"); + } + } + + #[test] + fn point_cloud_byte_length() { + let pc = make_point_cloud(&[[0.0, 0.0, 0.0], [1.0, 1.0, 1.0]]); + // 2 points × 3 floats × 4 bytes = 24 + assert_eq!(pc.byte_length(), 24); + } + + #[test] + fn point_cloud_empty() { + let pc = PointCloud { + xyz: vec![], + point_count: 0, + }; + assert_eq!(pc.byte_length(), 0); + assert_eq!(pc.point_count, 0); + } + + #[test] + #[should_panic] + fn point_cloud_position_out_of_bounds_panics() { + let pc = make_point_cloud(&[[1.0, 2.0, 3.0]]); + let _ = pc.position(1); // must panic + } +} diff --git a/crates/cesium/src/spz.rs b/crates/cesium/src/spz.rs new file mode 100644 index 00000000..d05712ee --- /dev/null +++ b/crates/cesium/src/spz.rs @@ -0,0 +1,462 @@ +//! `spz` — `KHR_gaussian_splatting_compression_spz` decode. +//! +//! Niantic SPZ is the binary compressed format ratified as +//! `KHR_gaussian_splatting_compression_spz` inside the Khronos glTF Gaussian Splatting +//! extension family (release-candidate as of 2026-02). It achieves ≈90% size reduction +//! vs equivalent PLY through per-attribute quantisation + per-stream ZSTD compression +//! (v4; earlier versions used gzip). +//! +//! # Format overview +//! Source: `github.com/nianticlabs/spz` (open source, fetched 2026-05-26). +//! All struct fields and constants below are grounded against `load-spz.cc` / `spz.h`. +//! +//! ```text +//! ┌─────────────────────────────────────────────────┐ +//! │ NgspFileHeader (32 bytes, plaintext, LE) │ +//! ├─────────────────────────────────────────────────┤ +//! │ optional extension blocks (tocByteOffset − 32)│ +//! ├─────────────────────────────────────────────────┤ +//! │ TOC entries (numStreams × 16 bytes) │ +//! │ each: [compressedSize: u64, uncompSize: u64] │ +//! ├─────────────────────────────────────────────────┤ +//! │ compressed stream 0: positions │ +//! │ compressed stream 1: alphas │ +//! │ compressed stream 2: colors (DC SH coeff) │ +//! │ compressed stream 3: scales │ +//! │ compressed stream 4: rotations │ +//! │ compressed stream 5: SH (degree 1–4 rest) │ +//! └─────────────────────────────────────────────────┘ +//! ``` +//! +//! # NgspFileHeader layout (32 bytes, little-endian) +//! ```text +//! offset size field +//! 0 4 magic — 0x5053474e ("NGSP" ASCII LE) +//! 4 4 version — currently 4 +//! 8 4 numPoints — Gaussian count +//! 12 1 shDegree — 0–4 +//! 13 1 fractionalBits — fixed-point precision (default 12) +//! 14 1 flags — FlagAntialiased=0x1, FlagHasExtensions=0x2 +//! 15 1 numStreams — typically 6 +//! 16 4 tocByteOffset — byte offset of TOC from file start +//! 20 12 reserved — must be zero +//! ``` +//! +//! # Quantisation details (grounded from `load-spz.cc`) +//! +//! ## Positions — stream 0 +//! - **9 bytes/point** (3 coords × 3 bytes each) +//! - 24-bit signed fixed-point, little-endian, sign-extended from bit 23 +//! - `real = fixed32 / (1 << fractionalBits)` where `fractionalBits=12` → ~0.25 mm +//! - Coordinate flip may be applied post-dequant (implementation-defined axis convention) +//! +//! ## Alphas — stream 1 +//! - **1 byte/point** (uint8) +//! - `opacity_logit = invSigmoid(u8 / 255.0)` where `invSigmoid(x) = ln(x / (1-x))` +//! - Maps linear opacity to logit space for GaussianBatch compatibility +//! +//! ## Colors (DC) — stream 2 +//! - **3 bytes/point** (uint8 × 3, RGB) +//! - `sh0[c] = (u8 / 255.0 - 0.5) / colorScale` where `colorScale = 0.15` +//! - These are SH degree-0 coefficients, not sRGB display values +//! +//! ## Scales — stream 3 +//! - **3 bytes/point** (uint8 × 3, log-encoded) +//! - `scale[j] = (u8 / 16.0) - 10.0` (log-space; exponentiate for Gaussian σ) +//! +//! ## Rotations — stream 4 +//! - **v3 and earlier (3 bytes/point):** xyz of unit quaternion as int8, derive w +//! - **v3+ / "smallest-three" (4 bytes/point):** +//! 2-bit index of the omitted (largest) component + three 10-bit signed magnitudes +//! packed into 32 bits. `MIN_SMALLEST_THREE_QUATERNIONS_VERSION = 3`. +//! +//! Decode sketch: +//! ```text +//! packed u32 → largest_idx (bits 30-31) + 3 × 10-bit signed components (bits 0-29) +//! reconstruct omitted component = sqrt(1 - dot(smallest_three, smallest_three)) +//! renormalise → unit quaternion [w, x, y, z] +//! ``` +//! // UNVERIFIED: exact bit-field layout within the 32-bit word; ground against +//! // nianticlabs/spz `unpackSmallestThreeQuaternion` before implementing. +//! +//! ## Spherical Harmonics — stream 5 +//! - Only present if `shDegree > 0` +//! - **Bit widths:** degree 1 → 5 bits/coeff (default), degrees 2–4 → 4 bits/coeff +//! (configurable via `sh1Bits`/`shRestBits` in range 1–8) +//! - In the stream as-stored all SH coefficients are repacked to **8 bits**: +//! `sh = (u8 - 128.0) / 128.0` +//! - Layout: for each SH coefficient index across all points, then across R/G/B channels +//! (column-major per coefficient, not interleaved per point) +//! - `shDim` = number of SH coefficients per channel: +//! degree 0 → 1, degree 1 → 4, degree 2 → 9, degree 3 → 16 +//! // UNVERIFIED: confirm shDim formula against spz source (degree 0 is in colors stream) +//! +//! # KHR_gaussian_splatting_compression_spz glTF binding +//! The SPZ blob appears as a glTF `bufferView` referenced from the mesh primitive's +//! extension object under `KHR_gaussian_splatting_compression_spz`. The extension +//! replaces the per-attribute accessor layout defined by `KHR_gaussian_splatting`; +//! the splat renderer reconstructs all attributes from the SPZ streams. +//! +//! # CAM SoA target +//! `spz_to_gaussian_batch` (agent C, `to_cam_soa`) consumes `SpzGaussianCloud` and +//! produces `GaussianBatch` + `cam_pq` indices. This module stops at the intermediate. +//! +//! # Grounding +//! - `github.com/nianticlabs/spz` (MIT, fetched 2026-05-26) +//! - Khronos `KHR_gaussian_splatting_compression_spz` PR #2490 (release-candidate 2026-02) +//! - Polyvia3D SPZ format overview (secondary confirmation) + +// ══════════════════════════════════════════════════════════════════════════════ +// Constants +// ══════════════════════════════════════════════════════════════════════════════ + +/// SPZ file magic: bytes `NGSP` in ASCII little-endian → `0x5053474e`. +/// +/// ```rust +/// // pub const NGSP_MAGIC: u32 = 0x5053_474e; +/// ``` +pub const NGSP_MAGIC: u32 = 0x5053_474e; + +/// Minimum SPZ version that uses "smallest-three" quaternion encoding. +/// Versions < 3 use raw int8 xyz + derived w (3 bytes/point). +/// +/// ```rust +/// // pub const MIN_SMALLEST_THREE_VERSION: u32 = 3; +/// ``` +pub const MIN_SMALLEST_THREE_VERSION: u32 = 3; + +/// Default fixed-point fractional bit count (yields ~0.25 mm resolution). +/// +/// ```rust +/// // pub const DEFAULT_FRACTIONAL_BITS: u8 = 12; +/// ``` +pub const DEFAULT_FRACTIONAL_BITS: u8 = 12; + +/// DC colour scaling factor applied when dequantising SH degree-0 coefficients. +/// `sh0 = (u8/255 - 0.5) / COLOR_SCALE`. +/// +/// ```rust +/// // pub const COLOR_SCALE: f32 = 0.15; +/// ``` +pub const COLOR_SCALE: f32 = 0.15; + +// ══════════════════════════════════════════════════════════════════════════════ +// Header +// ══════════════════════════════════════════════════════════════════════════════ + +/// Parsed SPZ file header (32 bytes, little-endian plaintext). +/// Grounded against `NgspFileHeader` in `nianticlabs/spz`. +/// +/// ```rust +/// // #[repr(C)] +/// // pub struct SpzHeader { +/// // pub magic: u32, // must == NGSP_MAGIC +/// // pub version: u32, // must be 1–4; 4 = ZSTD streams +/// // pub num_points: u32, // Gaussian count +/// // pub sh_degree: u8, // 0–4 +/// // pub fractional_bits: u8, // fixed-point precision +/// // pub flags: u8, // bit 0 = antialiased; bit 1 = has_extensions +/// // pub num_streams: u8, // typically 6 +/// // pub toc_byte_offset: u32, // TOC start offset from file start +/// // pub reserved: [u8; 12], // must be zero +/// // } +/// // +/// // impl SpzHeader { +/// // pub const SIZE: usize = 32; +/// // +/// // pub fn from_bytes(b: &[u8; 32]) -> Result { +/// // let magic = u32::from_le_bytes(b[0..4].try_into().unwrap()); +/// // if magic != NGSP_MAGIC { return Err(SpzError::BadMagic(magic)); } +/// // let version = u32::from_le_bytes(b[4..8].try_into().unwrap()); +/// // Ok(Self { +/// // magic, +/// // version, +/// // num_points: u32::from_le_bytes(b[8..12].try_into().unwrap()), +/// // sh_degree: b[12], +/// // fractional_bits: b[13], +/// // flags: b[14], +/// // num_streams: b[15], +/// // toc_byte_offset: u32::from_le_bytes(b[16..20].try_into().unwrap()), +/// // reserved: b[20..32].try_into().unwrap(), +/// // }) +/// // } +/// // +/// // pub fn has_extensions(&self) -> bool { self.flags & 0x2 != 0 } +/// // pub fn is_antialiased(&self) -> bool { self.flags & 0x1 != 0 } +/// // } +/// ``` +pub struct SpzHeader; + +/// One entry in the stream Table of Contents (16 bytes each). +/// Located at `header.toc_byte_offset`, with `header.num_streams` entries. +/// +/// ```rust +/// // pub struct SpzTocEntry { +/// // pub compressed_size: u64, +/// // pub uncompressed_size: u64, +/// // } +/// ``` +pub struct SpzTocEntry; + +// ══════════════════════════════════════════════════════════════════════════════ +// Decoded intermediate +// ══════════════════════════════════════════════════════════════════════════════ + +/// Fully decoded SPZ Gaussian cloud — intermediate struct consumed by `to_cam_soa`. +/// All values are in their natural floating-point representation post-dequantisation. +/// This struct must NOT outlive the import boundary. +/// +/// Layout is **SoA** (structure-of-arrays) to match CAM SoA target and enable SIMD +/// conversion passes in `to_cam_soa`. +/// +/// ```rust +/// // pub struct SpzGaussianCloud { +/// // pub num_points: usize, +/// // pub sh_degree: u8, +/// // pub antialiased: bool, +/// // +/// // /// Positions [x,y,z] — f32, world-space after fixed-point decode. +/// // /// Length: num_points * 3. +/// // pub positions: Vec, +/// // +/// // /// Opacity logits (invSigmoid of linear alpha) — f32. +/// // /// Length: num_points. +/// // pub opacity_logits: Vec, +/// // +/// // /// DC SH coefficients / base colour, logit-scaled — f32. +/// // /// Length: num_points * 3 (RGB interleaved: r0,g0,b0, r1,g1,b1, …). +/// // pub sh_dc: Vec, +/// // +/// // /// Log-space scale [sx, sy, sz] per Gaussian — f32. +/// // /// Length: num_points * 3. +/// // pub log_scales: Vec, +/// // +/// // /// Unit quaternions [w, x, y, z] per Gaussian — f32. +/// // /// Length: num_points * 4. +/// // pub rotations: Vec, +/// // +/// // /// Higher-degree SH coefficients (degrees 1–shDegree) — f32. +/// // /// Length: num_points * (shDim-1) * 3 where shDim = (shDegree+1)^2. +/// // /// UNVERIFIED: confirm exact shDim formula and channel-interleave order. +/// // pub sh_rest: Vec, +/// // } +/// ``` +pub struct SpzGaussianCloud; + +// ══════════════════════════════════════════════════════════════════════════════ +// Decode pipeline +// ══════════════════════════════════════════════════════════════════════════════ + +/// Top-level entry point: decode a raw SPZ blob into [`SpzGaussianCloud`]. +/// +/// # Steps +/// 1. Parse [`SpzHeader`] from bytes[0..32]. +/// 2. Parse TOC from `bytes[header.toc_byte_offset..]`. +/// 3. Locate stream data start = `toc_byte_offset + num_streams * 16`. +/// 4. For each stream: slice compressed bytes, decompress with ZSTD (v4) or gzip (v1–3). +/// 5. Decode each attribute stream (see below). +/// 6. Return `SpzGaussianCloud`. +/// +/// # Version compatibility +/// - v1–2: gzip compression, 3-byte/point rotation (int8 xyz) +/// - v3: gzip + smallest-three quaternion (4 bytes/point) +/// - v4: ZSTD per-stream + smallest-three quaternion +/// +/// ```rust +/// // pub fn decode_spz(bytes: &[u8]) -> Result { +/// // let hdr_bytes: &[u8; 32] = bytes[..32].try_into().map_err(|_| SpzError::TooShort)?; +/// // let hdr = SpzHeader::from_bytes(hdr_bytes)?; +/// // let toc = parse_toc(bytes, &hdr)?; +/// // let stream_data_start = hdr.toc_byte_offset as usize + hdr.num_streams as usize * 16; +/// // let streams = decompress_streams(bytes, &hdr, &toc, stream_data_start)?; +/// // Ok(SpzGaussianCloud { +/// // num_points: hdr.num_points as usize, +/// // sh_degree: hdr.sh_degree, +/// // antialiased: hdr.is_antialiased(), +/// // positions: decode_positions(&streams[0], &hdr)?, +/// // opacity_logits:decode_alphas(&streams[1], hdr.num_points as usize)?, +/// // sh_dc: decode_colors(&streams[2], hdr.num_points as usize)?, +/// // log_scales: decode_scales(&streams[3], hdr.num_points as usize)?, +/// // rotations: decode_rotations(&streams[4], &hdr)?, +/// // sh_rest: decode_sh_rest(&streams[5], &hdr)?, +/// // }) +/// // } +/// ``` +pub fn decode_spz(_bytes: &[u8]) -> Result { + Err(SpzError::NotImplemented("scaffold only — not yet implemented")) +} + +/// Decode stream 0: 24-bit fixed-point positions → f32 world coordinates. +/// +/// ```rust +/// // fn decode_positions(stream: &[u8], hdr: &SpzHeader) -> Result, SpzError> { +/// // let n = hdr.num_points as usize; +/// // assert_eq!(stream.len(), n * 9, "positions stream size mismatch"); +/// // let scale = 1.0_f32 / (1u32 << hdr.fractional_bits) as f32; +/// // let mut out = Vec::with_capacity(n * 3); +/// // for i in 0..n { +/// // for j in 0..3 { +/// // let base = i * 9 + j * 3; +/// // let mut v = stream[base] as u32 +/// // | ((stream[base+1] as u32) << 8) +/// // | ((stream[base+2] as u32) << 16); +/// // if v & 0x80_0000 != 0 { v |= 0xFF00_0000; } // sign extend +/// // out.push((v as i32) as f32 * scale); +/// // } +/// // } +/// // Ok(out) +/// // } +/// ``` +pub(crate) fn decode_positions(_stream: &[u8], _num_points: usize, _fractional_bits: u8) -> Vec { + unimplemented!("scaffold only") +} + +/// Decode stream 1: uint8 alpha → opacity logit (invSigmoid). +/// +/// ```rust +/// // fn decode_alphas(stream: &[u8], n: usize) -> Result, SpzError> { +/// // assert_eq!(stream.len(), n); +/// // Ok(stream.iter().map(|&a| { +/// // let p = a as f32 / 255.0; +/// // let p = p.clamp(1e-6, 1.0 - 1e-6); // guard against log(0) +/// // (p / (1.0 - p)).ln() // invSigmoid / logit +/// // }).collect()) +/// // } +/// ``` +pub(crate) fn decode_alphas(_stream: &[u8]) -> Vec { + unimplemented!("scaffold only") +} + +/// Decode stream 2: uint8 RGB → DC SH coefficients. +/// `sh0 = (u8 / 255.0 - 0.5) / COLOR_SCALE`. +/// +/// ```rust +/// // fn decode_colors(stream: &[u8], n: usize) -> Result, SpzError> { +/// // assert_eq!(stream.len(), n * 3); +/// // Ok(stream.iter().map(|&c| { +/// // (c as f32 / 255.0 - 0.5) / COLOR_SCALE +/// // }).collect()) +/// // } +/// ``` +pub(crate) fn decode_colors(_stream: &[u8]) -> Vec { + unimplemented!("scaffold only") +} + +/// Decode stream 3: uint8 log-scale → f32 log-scale values. +/// `log_scale = u8 / 16.0 - 10.0`. +/// +/// ```rust +/// // fn decode_scales(stream: &[u8], n: usize) -> Result, SpzError> { +/// // assert_eq!(stream.len(), n * 3); +/// // Ok(stream.iter().map(|&s| s as f32 / 16.0 - 10.0).collect()) +/// // } +/// ``` +pub(crate) fn decode_scales(_stream: &[u8]) -> Vec { + unimplemented!("scaffold only") +} + +/// Decode stream 4: quaternion bytes → [w, x, y, z] f32 unit quaternions. +/// +/// Two sub-cases depending on `header.version`: +/// - version < `MIN_SMALLEST_THREE_VERSION`: int8 xyz, `w = sqrt(1 - x²-y²-z²)` +/// - version >= `MIN_SMALLEST_THREE_VERSION`: smallest-three packed u32 +/// +/// Smallest-three decode (UNVERIFIED — confirm bit layout against spz source): +/// ```text +/// packed u32: +/// bits 31-30: index of omitted (largest magnitude) component +/// bits 29-20: component A (10-bit signed) +/// bits 19-10: component B (10-bit signed) +/// bits 9-0: component C (10-bit signed) +/// 10-bit signed scale factor: UNVERIFIED — ground against unpackSmallestThreeQuaternion +/// omitted = sqrt(1 - A²-B²-C²), sign = positive by convention +/// reconstruct full quaternion, renormalise +/// ``` +/// +/// ```rust +/// // fn decode_rotations(stream: &[u8], hdr: &SpzHeader) -> Result, SpzError> { +/// // let n = hdr.num_points as usize; +/// // if hdr.version < MIN_SMALLEST_THREE_VERSION { +/// // // legacy: 3 bytes/point, int8 xyz +/// // assert_eq!(stream.len(), n * 3); +/// // // ... decode xyz, derive w ... +/// // } else { +/// // // smallest-three: 4 bytes/point +/// // assert_eq!(stream.len(), n * 4); +/// // // ... unpack 2-bit index + 3×10-bit signed ... +/// // } +/// // } +/// ``` +pub(crate) fn decode_rotations(_stream: &[u8], _version: u32, _num_points: usize) -> Vec { + unimplemented!("scaffold only") +} + +/// Decode stream 5: higher-degree SH coefficients (degrees 1–shDegree). +/// `sh = (u8 - 128.0) / 128.0`. +/// +/// Stream size: `num_points × (shDim - 1) × 3` bytes +/// where `shDim = (shDegree + 1)^2` (UNVERIFIED: degree-0 already in colors stream). +/// +/// ```rust +/// // fn decode_sh_rest(stream: &[u8], hdr: &SpzHeader) -> Result, SpzError> { +/// // let sh_dim = ((hdr.sh_degree as usize + 1).pow(2)).saturating_sub(1); // UNVERIFIED +/// // let expected = hdr.num_points as usize * sh_dim * 3; +/// // assert_eq!(stream.len(), expected, "SH rest stream size mismatch"); +/// // Ok(stream.iter().map(|&v| (v as f32 - 128.0) / 128.0).collect()) +/// // } +/// ``` +pub(crate) fn decode_sh_rest(_stream: &[u8], _sh_degree: u8, _num_points: usize) -> Vec { + unimplemented!("scaffold only") +} + +// ══════════════════════════════════════════════════════════════════════════════ +// Error type +// ══════════════════════════════════════════════════════════════════════════════ + +/// Errors produced by the SPZ decode pipeline. +/// +/// ```rust +/// // #[derive(Debug)] +/// // pub enum SpzError { +/// // /// File is shorter than the 32-byte header. +/// // TooShort, +/// // /// Magic bytes do not match NGSP_MAGIC. +/// // BadMagic(u32), +/// // /// Unsupported version (>4). +/// // UnsupportedVersion(u32), +/// // /// ZSTD or gzip decompression failure. +/// // DecompressFailed(String), +/// // /// Stream length does not match expected size for attribute. +/// // StreamSizeMismatch { stream: usize, expected: usize, actual: usize }, +/// // } +/// ``` +#[derive(Debug)] +pub enum SpzError { + /// Function is a scaffold stub and has not been implemented yet. + NotImplemented(&'static str), +} + +impl std::fmt::Display for SpzError { + fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { + match self { + SpzError::NotImplemented(msg) => write!(f, "not implemented: {}", msg), + } + } +} + +impl std::error::Error for SpzError {} + +// ══════════════════════════════════════════════════════════════════════════════ +// Internal helpers — commented stubs +// ══════════════════════════════════════════════════════════════════════════════ + +// fn parse_toc(bytes: &[u8], hdr: &SpzHeader) -> Result, SpzError> { ... } +// fn decompress_streams(bytes: &[u8], hdr: &SpzHeader, toc: &[SpzTocEntry], start: usize) +// -> Result>, SpzError> { +// // v4: zstd::bulk::decompress per stream +// // v1-3: flate2::read::GzDecoder +// ... +// } +// +// NOTE: do NOT use gzip for v4 — the version field is the authoritative switch. +// NOTE: streams may be shorter than num_streams if shDegree=0 (stream 5 is absent/empty). diff --git a/crates/cesium/src/sse.rs b/crates/cesium/src/sse.rs new file mode 100644 index 00000000..5e484744 --- /dev/null +++ b/crates/cesium/src/sse.rs @@ -0,0 +1,303 @@ +//! `sse` (group C) — screen-space-error (SSE) LOD selection. +//! +//! Implements the OGC 3D Tiles / Cesium SSE formula used to decide whether a +//! tile's geometric error is small enough at the current camera to be accepted +//! as a leaf (no further refinement needed): +//! +//! ```text +//! sse = geometricError * viewportHeight / (distance * 2 * tan(fovy / 2)) +//! ``` +//! +//! A tile (or HLOD cluster) is accepted when `sse ≤ maximumScreenSpaceError` +//! (Cesium default: 16.0 pixels). +//! +//! # Grounding +//! - OGC 3D Tiles 1.1 spec §7.4 "Screen-Space Error" +//! - CesiumGS/cesium `Cesium3DTileset.maximumScreenSpaceError` (default 16) +//! - CesiumGS/cesium `Cesium3DTile._screenSpaceError` computation +//! +//! # UNVERIFIED items +//! - Whether Cesium additionally applies a `sseDenominatorFactor` or a +//! per-tile content-error multiplier in newer 3D Tiles 1.1 revisions. +//! - Exact behaviour when `distance ≈ 0` (camera inside tile bounding volume); +//! CesiumGS clamps to some minimum distance — exact value unconfirmed. +//! +//! # Commented scaffold only +//! All impl is `//`-commented until reviewed live by Opus + CodeRabbit. +//! The compilable stub types and their unit tests are live (std-only). + +// ───────────────────────────────────────────────────────────────────────────── +// Compilable stub types (std-only) +// ───────────────────────────────────────────────────────────────────────────── + +/// Camera frustum parameters needed for SSE computation. +/// +/// # Fields +/// - `viewport_height_px` — render target height in pixels (e.g. 1080.0). +/// - `fovy_rad` — vertical field-of-view in radians (e.g. π/3 ≈ 1.047 for 60°). +/// +/// # Pre-computation +/// The hot-path helper [`SseFrustum::sse_denominator`] pre-computes +/// `2 * tan(fovy/2)` once per frame so the per-tile call is a +/// multiply-divide only. +/// +/// # Examples +/// ``` +/// use cesium::sse::{SseFrustum, sse_for_tile, tile_meets_sse}; +/// +/// let frustum = SseFrustum { +/// viewport_height_px: 1080.0, +/// fovy_rad: std::f32::consts::FRAC_PI_3, +/// }; +/// let denom = frustum.sse_denominator(); // ≈ 935.307 +/// let sse = sse_for_tile(20.0, 500.0, denom); // ≈ 37.41 px +/// assert!((sse - 37.412).abs() < 0.1); +/// assert!(tile_meets_sse(1.0, 1000.0, denom, 16.0)); // tiny error — accepted +/// assert!(!tile_meets_sse(10_000.0, 1.0, denom, 16.0)); // huge error — rejected +/// ``` +#[derive(Clone, Copy, Debug, PartialEq)] +pub struct SseFrustum { + /// Render target height in pixels. + pub viewport_height_px: f32, + /// Vertical field-of-view in radians. + pub fovy_rad: f32, +} + +impl SseFrustum { + /// Pre-compute the constant denominator factor `viewportHeight / (2 * tan(fovy/2))`. + /// + /// Cache this value for the duration of a frame; call [`sse_for_tile`] with + /// the result to avoid recomputing `tan` per tile. + /// + /// # Examples + /// ``` + /// use cesium::sse::SseFrustum; + /// + /// let frustum = SseFrustum { + /// viewport_height_px: 1080.0, + /// fovy_rad: std::f32::consts::FRAC_PI_3, + /// }; + /// let denom = frustum.sse_denominator(); + /// // 1080 / (2 * tan(π/6)) ≈ 935.307 + /// assert!((denom - 935.307).abs() < 0.1); + /// ``` + #[inline] + pub fn sse_denominator(self) -> f32 { + // 2 * tan(fovy / 2) — the angular subtension of the viewport + let two_tan_half_fov = 2.0 * (self.fovy_rad * 0.5).tan(); + // UNVERIFIED: guard against zero fovy (degenerate frustum) + if two_tan_half_fov.abs() < f32::EPSILON { + return f32::MAX; + } + self.viewport_height_px / two_tan_half_fov + } +} + +/// Compute the screen-space error for one tile. +/// +/// # Arguments +/// - `geometric_error` — tile's `geometricError` from `tileset.json` (world units). +/// - `distance` — distance from camera to tile bounding volume (world units, +/// must be positive). Values ≤ 0 are clamped to a small epsilon to avoid +/// division by zero. +/// - `sse_denominator` — pre-computed `viewportHeight / (2 * tan(fovy/2))`; +/// obtain once per frame from [`SseFrustum::sse_denominator`]. +/// +/// # Returns +/// Screen-space error in pixels. Compare against +/// `maximumScreenSpaceError` (Cesium default 16.0) to decide refinement. +/// +/// # Formula +/// ```text +/// sse = geometricError * (viewportHeight / (distance * 2 * tan(fovy/2))) +/// = geometricError * sse_denominator / distance +/// ``` +/// +/// # Examples +/// ``` +/// use cesium::sse::{SseFrustum, sse_for_tile}; +/// +/// let frustum = SseFrustum { +/// viewport_height_px: 1080.0, +/// fovy_rad: std::f32::consts::FRAC_PI_3, +/// }; +/// let denom = frustum.sse_denominator(); // ≈ 935.307 +/// let sse = sse_for_tile(20.0, 500.0, denom); // 20 * 935.307 / 500 ≈ 37.41 +/// assert!((sse - 37.412).abs() < 0.1); +/// // zero geometric error → zero SSE regardless of distance +/// assert_eq!(sse_for_tile(0.0, 500.0, denom), 0.0); +/// ``` +#[inline] +pub fn sse_for_tile(geometric_error: f32, distance: f32, sse_denominator: f32) -> f32 { + // Clamp distance to avoid division by zero; UNVERIFIED: Cesium's exact + // minimum-distance clamp value. + let d = distance.max(1e-6_f32); + geometric_error * sse_denominator / d +} + +/// Returns `true` when the tile meets the LOD acceptance criterion. +/// +/// Equivalent to `sse_for_tile(...) <= maximum_sse`. +/// +/// # Examples +/// ``` +/// use cesium::sse::{SseFrustum, tile_meets_sse}; +/// +/// let frustum = SseFrustum { +/// viewport_height_px: 1080.0, +/// fovy_rad: std::f32::consts::FRAC_PI_3, +/// }; +/// let denom = frustum.sse_denominator(); +/// // small geometric error far away — SSE ≈ 0.935 px ≤ 16.0 → accepted +/// assert!(tile_meets_sse(1.0, 1000.0, denom, 16.0)); +/// // large geometric error up close — SSE enormous → rejected +/// assert!(!tile_meets_sse(10_000.0, 1.0, denom, 16.0)); +/// ``` +#[inline] +pub fn tile_meets_sse(geometric_error: f32, distance: f32, sse_denominator: f32, maximum_sse: f32) -> bool { + sse_for_tile(geometric_error, distance, sse_denominator) <= maximum_sse +} + +// ───────────────────────────────────────────────────────────────────────────── +// COMMENTED-OUT extended helpers (no external deps required, but kept commented +// until the full HLOD traversal is reviewed live) +// ───────────────────────────────────────────────────────────────────────────── +// +// ```rust +// /// Approximate distance from a camera position to a tile's bounding-sphere +// /// centre, minus the sphere radius (so the distance to the *surface*). +// /// +// /// Returns 0.0 if the camera is inside the sphere. +// /// +// /// UNVERIFIED: Cesium uses oriented bounding box (OBB) or bounding sphere +// /// depending on the tile's `boundingVolume` type. Sphere variant shown here. +// #[inline] +// pub fn distance_to_bounding_sphere( +// camera: [f32; 3], +// sphere_center: [f32; 3], +// sphere_radius: f32, +// ) -> f32 { +// let dx = camera[0] - sphere_center[0]; +// let dy = camera[1] - sphere_center[1]; +// let dz = camera[2] - sphere_center[2]; +// let dist_to_center = (dx*dx + dy*dy + dz*dz).sqrt(); +// (dist_to_center - sphere_radius).max(0.0) +// } +// ``` + +// ───────────────────────────────────────────────────────────────────────────── +// Unit tests (std-only, live) +// ───────────────────────────────────────────────────────────────────────────── + +#[cfg(test)] +mod tests { + use super::*; + + fn approx(a: f32, b: f32, tol: f32) -> bool { + (a - b).abs() <= tol + } + + /// Reference computation using the full formula directly (no pre-computation). + fn sse_reference(geometric_error: f32, viewport_height: f32, distance: f32, fovy_rad: f32) -> f32 { + let d = distance.max(1e-6); + geometric_error * viewport_height / (d * 2.0 * (fovy_rad * 0.5).tan()) + } + + #[test] + fn sse_denominator_matches_reference() { + let frustum = SseFrustum { + viewport_height_px: 1080.0, + fovy_rad: std::f32::consts::FRAC_PI_3, + }; + let denom = frustum.sse_denominator(); + // reference: 1080 / (2 * tan(π/6)) = 1080 / (2 * 0.57735) ≈ 935.3 + let expected = 1080.0_f32 / (2.0 * (std::f32::consts::FRAC_PI_3 / 2.0).tan()); + assert!(approx(denom, expected, 0.1), "denom={denom} expected={expected}"); + } + + #[test] + fn sse_for_tile_matches_reference() { + let ge = 20.0_f32; + let vh = 1080.0_f32; + let fovy = std::f32::consts::FRAC_PI_3; + let dist = 500.0_f32; + + let frustum = SseFrustum { + viewport_height_px: vh, + fovy_rad: fovy, + }; + let denom = frustum.sse_denominator(); + let got = sse_for_tile(ge, dist, denom); + let expected = sse_reference(ge, vh, dist, fovy); + + assert!(approx(got, expected, 1e-3), "got={got} expected={expected}"); + } + + #[test] + fn tile_meets_sse_accepts_small_error() { + // geometricError=1.0, distance=1000.0 → very small SSE → accepted + let frustum = SseFrustum { + viewport_height_px: 1080.0, + fovy_rad: std::f32::consts::FRAC_PI_3, + }; + let denom = frustum.sse_denominator(); + assert!(tile_meets_sse(1.0, 1000.0, denom, 16.0)); + } + + #[test] + fn tile_meets_sse_rejects_large_error() { + // geometricError=10000.0, distance=1.0 → enormous SSE → rejected + let frustum = SseFrustum { + viewport_height_px: 1080.0, + fovy_rad: std::f32::consts::FRAC_PI_3, + }; + let denom = frustum.sse_denominator(); + assert!(!tile_meets_sse(10000.0, 1.0, denom, 16.0)); + } + + #[test] + fn zero_geometric_error_always_accepted() { + let frustum = SseFrustum { + viewport_height_px: 1080.0, + fovy_rad: std::f32::consts::FRAC_PI_3, + }; + let denom = frustum.sse_denominator(); + // geometricError=0 → SSE=0 → always meets any positive threshold + assert!(tile_meets_sse(0.0, 1.0, denom, 16.0)); + } + + #[test] + fn very_small_distance_does_not_panic() { + let frustum = SseFrustum { + viewport_height_px: 1080.0, + fovy_rad: std::f32::consts::FRAC_PI_3, + }; + let denom = frustum.sse_denominator(); + let sse = sse_for_tile(10.0, 0.0, denom); // distance=0 → clamped + assert!(sse.is_finite(), "SSE must be finite even at zero distance"); + } + + #[test] + fn sse_scales_linearly_with_geometric_error() { + let frustum = SseFrustum { + viewport_height_px: 1080.0, + fovy_rad: std::f32::consts::FRAC_PI_3, + }; + let denom = frustum.sse_denominator(); + let s1 = sse_for_tile(10.0, 100.0, denom); + let s2 = sse_for_tile(20.0, 100.0, denom); + assert!(approx(s2, 2.0 * s1, 1e-4), "s1={s1} s2={s2}"); + } + + #[test] + fn sse_inversely_proportional_to_distance() { + let frustum = SseFrustum { + viewport_height_px: 1080.0, + fovy_rad: std::f32::consts::FRAC_PI_3, + }; + let denom = frustum.sse_denominator(); + let s1 = sse_for_tile(10.0, 100.0, denom); + let s2 = sse_for_tile(10.0, 200.0, denom); + assert!(approx(s1, 2.0 * s2, 1e-4), "s1={s1} s2={s2}"); + } +} diff --git a/crates/cesium/src/tileset.rs b/crates/cesium/src/tileset.rs new file mode 100644 index 00000000..c4664ba1 --- /dev/null +++ b/crates/cesium/src/tileset.rs @@ -0,0 +1,353 @@ +//! `tileset` (group A) — OGC 3D Tiles `tileset.json` / `.3tz` parse (cold import → tile-tree model) +//! +//! # Grounding +//! - OGC 3D Tiles 1.1 specification (CesiumGS/3d-tiles, 22-025r4) +//! - Schema: `tile.schema.json`, `tileset.schema.json`, `content.schema.json`, +//! `boundingVolume.schema.json`, `tile.implicitTiling.schema.json` +//! +//! # Design contract +//! - **COLD IMPORT ONLY.** `tileset.json` is parsed exactly once at load time. +//! No JSON value ever reaches the hotpath; all JSON-derived data is immediately +//! converted to owned intermediate structs defined below. +//! - **Intermediate structs only.** No source-native types survive past the import +//! boundary. The parsed tree maps directly to [`to_cam_soa`] bridge inputs. +//! - **No live implementation.** All planned code is `//`-commented scaffold; +//! nothing here is `unsafe` or `pub fn` yet. Reviewed by Opus + CodeRabbit +//! before any impl is uncommented. +//! +//! # Module responsibilities +//! 1. Parse `tileset.json` (or `.3tz` ZIP envelope) cold → [`Tileset`] + [`TileNode`] tree. +//! 2. Resolve `implicitTiling` back-references so implicit tiles become concrete [`TileNode`] +//! entries (deferred to [`implicit_tiling`] — only the hook lives here). +//! 3. Expose a flat tile iterator (`iter_leaves`) that the SSE module uses for LOD +//! selection without any JSON involvement. +//! +//! [`to_cam_soa`]: crate::to_cam_soa +//! [`implicit_tiling`]: crate::implicit_tiling + +// ───────────────────────────────────────────────────────────────────────────── +// Planned public types (all COMMENTED OUT — not yet live) +// ───────────────────────────────────────────────────────────────────────────── + +// ┌─────────────────────────────────────────────────────────────────────────┐ +// │ Top-level tileset document │ +// │ │ +// │ Fields verified against tileset.schema.json (CesiumGS/3d-tiles main): │ +// │ asset — object, required │ +// │ geometricError — f64 ≥ 0, required │ +// │ root — TileNode, required │ +// │ schema — optional metadata class map │ +// │ schemaUri — optional URI string │ +// │ statistics — optional │ +// │ groups — optional array │ +// │ metadata — optional metadataEntity │ +// │ extensionsUsed — Vec │ +// │ extensionsRequired — Vec │ +// └─────────────────────────────────────────────────────────────────────────┘ +// +// ```rust +// /// Top-level parsed representation of a `tileset.json` document. +// /// Constructed once at cold-import time; never re-parsed. +// pub struct Tileset { +// /// Asset metadata (version string, tilesetVersion, etc.). +// pub asset: TilesetAsset, +// /// Maximum geometric error for the entire tileset (meters). +// /// Source field: `geometricError` (f64 ≥ 0). +// pub geometric_error: f64, +// /// Root of the tile tree. Children are stored inline (recursive). +// pub root: TileNode, +// /// Extensions declared in `extensionsUsed`. +// pub extensions_used: Vec, +// /// Extensions declared in `extensionsRequired`. +// pub extensions_required: Vec, +// } +// ``` + +// ┌─────────────────────────────────────────────────────────────────────────┐ +// │ Asset metadata │ +// │ │ +// │ Fields from asset.schema.json: │ +// │ version — string, required ("1.1") │ +// │ tilesetVersion — optional string (dataset version, not spec version) │ +// └─────────────────────────────────────────────────────────────────────────┘ +// +// ```rust +// pub struct TilesetAsset { +// /// Spec version string, e.g. `"1.1"`. +// pub version: String, +// /// Optional dataset/revision version tag. +// pub tileset_version: Option, +// } +// ``` + +// ┌─────────────────────────────────────────────────────────────────────────┐ +// │ Bounding volume │ +// │ │ +// │ Verified against boundingVolume.schema.json — exactly one of: │ +// │ box — [f64; 12]: center (3) + half-axes (3×3) column-major │ +// │ region — [f64; 6]: [west, south, east, north] radians + [min,max] m │ +// │ sphere — [f64; 4]: center (3) + radius │ +// │ Only one key may be present per object (spec: anyOf). │ +// └─────────────────────────────────────────────────────────────────────────┘ +// +// ```rust +// /// Parsed bounding volume — one of three mutually exclusive shapes. +// pub enum BoundingVolume { +// /// `boundingVolume.box`: [cx,cy,cz, hx0,hx1,hx2, hy0,hy1,hy2, hz0,hz1,hz2] +// Box([f64; 12]), +// /// `boundingVolume.region`: [west,south,east,north] (rad) + [minH,maxH] (m) +// Region([f64; 6]), +// /// `boundingVolume.sphere`: [cx,cy,cz,radius] +// Sphere([f64; 4]), +// } +// ``` + +// ┌─────────────────────────────────────────────────────────────────────────┐ +// │ Refine mode │ +// │ │ +// │ Verified: tile.schema.json `refine` enum — only "ADD" and "REPLACE" │ +// └─────────────────────────────────────────────────────────────────────────┘ +// +// ```rust +// /// Refinement strategy for this tile's children. +// pub enum Refine { +// /// ADD — render both parent and children simultaneously. +// Add, +// /// REPLACE — render children instead of parent when refined. +// Replace, +// } +// ``` + +// ┌─────────────────────────────────────────────────────────────────────────┐ +// │ Tile content │ +// │ │ +// │ Verified against content.schema.json: │ +// │ uri — string, required (was `url` in 1.0, renamed in 1.1) │ +// │ boundingVolume — optional, same shape as tile bounding volume │ +// │ metadata — optional metadataEntity │ +// │ │ +// │ A tile may have EITHER `content` (singular) OR `contents` (array ≥ 1) │ +// │ but NOT both (spec: mutually exclusive). │ +// └─────────────────────────────────────────────────────────────────────────┘ +// +// ```rust +// /// A single content reference within a tile. +// pub struct TileContent { +// /// Relative or absolute URI to the tile content payload +// /// (e.g. `.glb`, `.b3dm`, `.pnts`, `.cmpt`). +// /// Source field: `uri` (renamed from `url` in 3D Tiles 1.1). +// pub uri: String, +// /// Optional tighter bounding volume for this content within the tile. +// pub bounding_volume: Option, +// } +// ``` + +// ┌─────────────────────────────────────────────────────────────────────────┐ +// │ Implicit tiling descriptor (cold-parse only) │ +// │ │ +// │ Verified against tile.implicitTiling.schema.json — all four required: │ +// │ subdivisionScheme — enum "QUADTREE" | "OCTREE" │ +// │ subtreeLevels — u32 ≥ 1 (levels per subtree binary file) │ +// │ availableLevels — u32 ≥ 1 (total levels with available tiles) │ +// │ subtrees — object with `uri` template string │ +// │ QUADTREE template vars: {level} {x} {y} │ +// │ OCTREE template vars: {level} {x} {y} {z} │ +// └─────────────────────────────────────────────────────────────────────────┘ +// +// ```rust +// /// Implicit tiling configuration attached to a tile node. +// /// Expanding implicit tiles is delegated to [`crate::implicit_tiling`]. +// pub struct ImplicitTilingRef { +// pub scheme: SubdivisionScheme, +// /// Number of levels stored in a single subtree binary file (≥ 1). +// pub subtree_levels: u32, +// /// Total number of available levels across all subtrees (≥ 1). +// pub available_levels: u32, +// /// URI template for subtree files. +// /// QUADTREE: `{level}/{x}/{y}.subtree` +// /// OCTREE: `{level}/{x}/{y}/{z}.subtree` +// pub subtrees_uri_template: String, +// } +// +// pub enum SubdivisionScheme { +// Quadtree, +// Octree, +// } +// ``` + +// ┌─────────────────────────────────────────────────────────────────────────┐ +// │ Tile node (recursive tree) │ +// │ │ +// │ Verified fields from tile.schema.json: │ +// │ boundingVolume — required │ +// │ geometricError — f64 ≥ 0, required │ +// │ viewerRequestVolume — optional │ +// │ refine — optional enum ADD|REPLACE (inherited if absent) │ +// │ content — optional (mutually exclusive with `contents`) │ +// │ contents — optional array ≥ 1 (multi-content) │ +// │ children — optional array of child TileNodes ≥ 1 │ +// │ transform — optional [f64; 16] column-major 4×4 affine │ +// │ implicitTiling — optional (marks implicit root tile) │ +// └─────────────────────────────────────────────────────────────────────────┘ +// +// ```rust +// /// A single node in the parsed tile tree. +// /// All JSON-derived strings are owned; no serde values survive here. +// pub struct TileNode { +// pub bounding_volume: BoundingVolume, +// /// Geometric error in meters (used by SSE LOD selection). +// pub geometric_error: f64, +// /// Optional tighter volume within which a viewer can request this tile. +// pub viewer_request_volume: Option, +// /// Refinement mode (ADD or REPLACE); None means "inherit from parent". +// pub refine: Option, +// /// Column-major 4×4 affine transform (local → parent). None = identity. +// /// Source field: `transform` ([f64; 16]). +// pub transform: Option<[f64; 16]>, +// /// Content payloads — zero, one, or many. +// pub contents: Vec, +// /// Explicit children (empty for leaf tiles and implicit roots). +// pub children: Vec, +// /// Populated for implicit root tiles; expanded by [`crate::implicit_tiling`]. +// pub implicit_tiling: Option, +// } +// ``` + +// ───────────────────────────────────────────────────────────────────────────── +// Planned functions (all COMMENTED OUT) +// ───────────────────────────────────────────────────────────────────────────── + +// ┌─────────────────────────────────────────────────────────────────────────┐ +// │ Cold parse entry point │ +// │ │ +// │ Accepts raw bytes of `tileset.json` or a `.3tz` ZIP envelope. │ +// │ JSON is decoded here and ONLY here — no JSON value escapes this fn. │ +// │ │ +// │ .3tz format: a ZIP64 archive containing `tileset.json` at the root │ +// │ plus relative tile payloads. The archive is a "3D Tiles Archive" (3TZ).│ +// │ UNVERIFIED: exact .3tz magic bytes and central-directory layout; │ +// │ ground against OGC 3D Tiles 1.1 §13 before uncommenting. │ +// └─────────────────────────────────────────────────────────────────────────┘ +// +// ```rust +// /// Parse a `tileset.json` byte slice (UTF-8) into a [`Tileset`] tree. +// /// +// /// This is the ONLY entry point where JSON is parsed. After this call no +// /// JSON/serde values remain alive. Cheap to call once; do not call in +// /// a loop or on the hotpath. +// /// +// /// # Errors +// /// Returns [`ParseError`] on malformed JSON, missing required fields, or +// /// unknown enum values (fail-fast — do not silently default). +// pub fn parse_tileset(bytes: &[u8]) -> Result { +// // 1. Decode UTF-8. +// // 2. Pass to minimal JSON tokeniser (no external dep until uncommented). +// // 3. Walk object keys, populate Tileset fields. +// // 4. Recursively parse `root` → TileNode via parse_tile_node(). +// // 5. Drop all intermediate JSON state before returning. +// todo!() +// } +// +// /// Parse a `.3tz` ZIP64 archive: locate `tileset.json` entry, extract it, +// /// then delegate to [`parse_tileset`]. +// /// +// /// UNVERIFIED: .3tz central-directory magic and entry naming convention; +// /// ground against OGC 22-025r4 §13 / CesiumGS/3d-tiles before enabling. +// pub fn parse_3tz(archive_bytes: &[u8]) -> Result { +// // 1. Locate ZIP64 end-of-central-directory record. +// // 2. Find entry named "tileset.json" (case-sensitive per spec). +// // 3. Decompress (DEFLATE or STORE). +// // 4. Call parse_tileset() on decompressed bytes. +// todo!() +// } +// +// /// Recursively parse a tile JSON object into a [`TileNode`]. +// /// Internal — called by parse_tileset() for `root` and every element of +// /// `children`. Not pub: callers always enter via [`parse_tileset`]. +// fn parse_tile_node(/* json object value */) -> Result { +// // Fields parsed in order: +// // boundingVolume (required) +// // geometricError (required, f64 ≥ 0) +// // viewerRequestVolume (optional) +// // refine (optional "ADD"|"REPLACE") +// // transform (optional [f64;16]) +// // content (optional, mutually exclusive with contents) +// // contents (optional array) +// // implicitTiling (optional → ImplicitTilingRef) +// // children (optional array, recursive) +// todo!() +// } +// +// /// Parse a `boundingVolume` JSON object into [`BoundingVolume`]. +// /// Exactly one of `box`, `region`, `sphere` must be present. +// fn parse_bounding_volume(/* json object */) -> Result { +// todo!() +// } +// +// /// Parse a `content` JSON object into [`TileContent`]. +// fn parse_content(/* json object */) -> Result { +// // `uri` field required (renamed from `url` in 3D Tiles 1.1 — accept +// // both spellings for 1.0 compat, prefer `uri`). +// todo!() +// } +// ``` + +// ┌─────────────────────────────────────────────────────────────────────────┐ +// │ Tile tree traversal │ +// │ │ +// │ Flat iterator over leaf TileNodes (explicit) or implicit-root stubs. │ +// │ Used by `sse` for LOD selection — no JSON involved. │ +// └─────────────────────────────────────────────────────────────────────────┘ +// +// ```rust +// impl Tileset { +// /// Iterator over all leaf [`TileNode`]s (nodes with no explicit children +// /// and no `implicitTiling`). Implicit roots are yielded as-is; callers +// /// must expand them via [`crate::implicit_tiling::expand_implicit`]. +// /// +// /// Traversal order: pre-order depth-first (matches CesiumJS traversal). +// pub fn iter_leaves(&self) -> impl Iterator { +// // Iterative DFS over self.root using a stack of &TileNode. +// // Yield node if node.children.is_empty(). +// todo!() +// } +// +// /// Walk the tile tree and collect every content URI reachable from +// /// explicit (non-implicit) nodes. Implicit subtree URIs are NOT +// /// expanded here — use [`crate::implicit_tiling`] for that. +// pub fn collect_content_uris(&self) -> Vec<&str> { +// todo!() +// } +// } +// ``` + +// ───────────────────────────────────────────────────────────────────────────── +// Error type (COMMENTED OUT) +// ───────────────────────────────────────────────────────────────────────────── +// +// ```rust +// /// Cold-parse errors for tileset.json / .3tz. +// /// Detailed enough for CLI diagnostics; not used in any hotpath. +// pub enum ParseError { +// /// Input is not valid UTF-8. +// InvalidUtf8, +// /// JSON tokenisation failed at byte offset. +// JsonSyntax(usize), +// /// Required field absent. +// MissingField(&'static str), +// /// Field value had wrong type. +// WrongType { field: &'static str, expected: &'static str }, +// /// `refine` value was not "ADD" or "REPLACE". +// UnknownRefine(String), +// /// `subdivisionScheme` was not "QUADTREE" or "OCTREE". +// UnknownSubdivisionScheme(String), +// /// Both `content` and `contents` are present (spec violation). +// ConflictingContent, +// /// .3tz archive entry `tileset.json` not found. +// /// UNVERIFIED: whether the spec mandates this exact filename or +// /// allows arbitrary roots inside the archive. +// ArchiveEntryNotFound, +// /// Decompression failure inside .3tz. +// DecompressError, +// } +// ``` diff --git a/crates/cesium/src/to_cam_soa.rs b/crates/cesium/src/to_cam_soa.rs new file mode 100644 index 00000000..62ddadf0 --- /dev/null +++ b/crates/cesium/src/to_cam_soa.rs @@ -0,0 +1,530 @@ +//! `to_cam_soa` (group C) — Rule-2 bridge: reverse-engineer parsed KHR splats +//! into ndarray CAM SoA (`splat3d::gaussian::GaussianBatch` + `cam_pq` codebook +//! indices). This is the canonical import target — the only module allowed to +//! produce `GaussianBatch` from external data. +//! +//! # Design contract +//! - **Reverse-engineer, never copy.** The parsed intermediate ([`khr_gs::GsSplats`]) +//! is restructured into SoA; the original attribute layout is abandoned at the +//! boundary. +//! - **No JSON in the hotpath** (rule 3). All JSON parsing happened in +//! [`crate::khr_gs`]; by the time we touch data here it is already `Vec`. +//! - **Commented scaffold only** — all impl touching `ndarray` is `//`-commented +//! until the `ndarray` dep in `Cargo.toml` is uncommented and reviewed live by +//! Opus + CodeRabbit. Compilable stub structs (std-only) and their unit tests +//! are live. +//! +//! # GaussianBatch field mapping (grounded against gaussian.rs) +//! +//! ```text +//! GsSplats field GaussianBatch field Notes +//! ─────────────────────── ────────────────────────── ────────────────────────────────── +//! position[i*3+0] → mean_x[i] +//! position[i*3+1] → mean_y[i] +//! position[i*3+2] → mean_z[i] +//! scale[i*3+0] → scale_x[i] +//! scale[i*3+1] → scale_y[i] +//! scale[i*3+2] → scale_z[i] +//! rotation_xyzw[i*4+3] → quat_w[i] *** glTF xyzw → splat3d wxyz *** +//! rotation_xyzw[i*4+0] → quat_x[i] x stays x +//! rotation_xyzw[i*4+1] → quat_y[i] y stays y +//! rotation_xyzw[i*4+2] → quat_z[i] z stays z +//! opacity[i] → opacity[i] +//! sh_coefficients[..] → sh[i*48 + ch*16 + basis] see SH repack below +//! ``` +//! +//! # SH repack (khr_gs layout → GaussianBatch layout) +//! +//! `GsSplats::sh_coefficients` stores coefficients contiguously per splat in +//! *attribute order* (one VEC3 per coefficient: [r,g,b] interleaved): +//! +//! ```text +//! khr_gs per-splat layout (degree-3 = 16 coefs × 3 channels = 48 f32): +//! [deg0_c0_r, deg0_c0_g, deg0_c0_b, ← basis k=0, channels rgb +//! deg1_c0_r, deg1_c0_g, deg1_c0_b, ← basis k=1 +//! deg1_c1_r, ..., +//! ... +//! deg3_c6_r, deg3_c6_g, deg3_c6_b] ← basis k=15 +//! ``` +//! +//! `GaussianBatch::sh` stores coefficients in *gaussian-major, channel-major*: +//! +//! ```text +//! sh[i * 48 + ch * 16 + k] ch ∈ {0=R, 1=G, 2=B}, k ∈ 0..16 +//! ``` +//! +//! The repack loop is therefore: +//! `batch.sh[i*48 + 0*16 + k] = gs_sh[k*3 + 0]` (red) +//! `batch.sh[i*48 + 1*16 + k] = gs_sh[k*3 + 1]` (green) +//! `batch.sh[i*48 + 2*16 + k] = gs_sh[k*3 + 2]` (blue) +//! where `gs_sh` is the 48-f32 slice for splat i. +//! +//! # CAM-PQ fingerprint (grounded against cam_pq.rs) +//! +//! After the SoA is built, the bridge optionally encodes splat position vectors +//! through `CamCodebook::encode` to produce a [`CamFingerprint`] (`[u8;6]`) +//! per splat. Encoding requires a pre-trained codebook; callers that have no +//! codebook receive `None` for the fingerprint slice. +//! +//! # UNVERIFIED items +//! - Whether Cesium-exported KHR files always normalise the ROTATION quaternion +//! before writing, or whether we must renormalise after dequantisation. +//! - Whether COLOR_0 is already linearised or is sRGB-encoded in practice. + +// ───────────────────────────────────────────────────────────────────────────── +// Compilable stub types (std-only, no ndarray dep) +// ───────────────────────────────────────────────────────────────────────────── + +/// Outcome of [`convert_to_cam_soa`] (stub — real fields gated behind ndarray). +/// +/// When the ndarray dep is live, `batch` becomes +/// `splat3d::gaussian::GaussianBatch` and `cam_fingerprints` becomes +/// `Vec`. +/// +/// # Examples +/// +/// ``` +/// use cesium::to_cam_soa::{convert_to_cam_soa, CamSoaResult}; +/// +/// let result: CamSoaResult = convert_to_cam_soa( +/// 1, +/// &[0.0, 0.0, 0.0], +/// &[0.0, 0.0, 0.0, 1.0], +/// &[1.0, 1.0, 1.0], +/// &[0.5], +/// None, +/// ) +/// .unwrap(); +/// assert_eq!(result.splat_count, 1); +/// ``` +#[derive(Debug)] +pub struct CamSoaResult { + /// Number of Gaussian splats converted. + pub splat_count: usize, + /// UNVERIFIED: placeholder for the live GaussianBatch field once ndarray dep + /// is uncommented. True type: `splat3d::gaussian::GaussianBatch`. + pub _batch_placeholder: (), + /// UNVERIFIED: placeholder for per-splat CAM fingerprints once cam_pq dep + /// is uncommented. True type: `Vec` = `Vec<[u8;6]>`. + pub _cam_fingerprints_placeholder: Option<()>, +} + +/// Errors emitted by the conversion bridge. +/// +/// # Examples +/// +/// ``` +/// use cesium::to_cam_soa::{convert_to_cam_soa, ToCamSoaError}; +/// +/// // Empty splat count → EmptySplats +/// let err = convert_to_cam_soa(0, &[], &[], &[], &[], None).unwrap_err(); +/// assert_eq!(err, ToCamSoaError::EmptySplats); +/// +/// // Wrong position buffer length → BufferLengthMismatch +/// let err = convert_to_cam_soa( +/// 1, +/// &[0.0], // too short: expected 3 (= 1*3) +/// &[0.0, 0.0, 0.0, 1.0], +/// &[1.0, 1.0, 1.0], +/// &[0.5], +/// None, +/// ) +/// .unwrap_err(); +/// assert!(matches!( +/// err, +/// ToCamSoaError::BufferLengthMismatch { field: "position", expected: 3, actual: 1 } +/// )); +/// ``` +#[derive(Debug, PartialEq, Eq)] +pub enum ToCamSoaError { + /// The source splat count is zero — nothing to convert. + EmptySplats, + /// A per-splat input buffer has the wrong length. + /// + /// `expected` is `splat_count * stride` (e.g. 3 for position/scale, 4 for + /// rotation, 1 for opacity). `actual` is `slice.len()`. + BufferLengthMismatch { + /// Name of the offending field (`"position"`, `"rotation_xyzw"`, + /// `"scale"`, or `"opacity"`). + field: &'static str, + /// Expected length (`splat_count * per-splat stride`). + expected: usize, + /// Actual length of the slice that was passed in. + actual: usize, + }, + /// SH coefficient array length is inconsistent with `splat_count * sh_stride`. + ShLengthMismatch { + splat_count: usize, + sh_len: usize, + expected_stride: usize, + }, + /// UNVERIFIED: quaternion norm was below threshold after dequantisation. + DegenerateRotation { splat_index: usize }, +} + +impl core::fmt::Display for ToCamSoaError { + fn fmt(&self, f: &mut core::fmt::Formatter<'_>) -> core::fmt::Result { + match self { + Self::EmptySplats => write!(f, "to_cam_soa: source splat count is zero"), + Self::BufferLengthMismatch { + field, + expected, + actual, + } => write!(f, "to_cam_soa: buffer `{field}` length {actual} != expected {expected}",), + Self::ShLengthMismatch { + splat_count, + sh_len, + expected_stride, + } => write!( + f, + "to_cam_soa: sh_coefficients length {} != splat_count {} * stride {}", + sh_len, splat_count, expected_stride, + ), + Self::DegenerateRotation { splat_index } => { + write!(f, "to_cam_soa: degenerate (near-zero-norm) quaternion at splat {}", splat_index,) + } + } + } +} + +// ───────────────────────────────────────────────────────────────────────────── +// Live stub function (returns a placeholder result; real body commented below) +// ───────────────────────────────────────────────────────────────────────────── + +/// Convert parsed KHR gaussian splats into ndarray CAM SoA. +/// +/// **Stub.** Returns `Ok(CamSoaResult { splat_count, .. })` so callers +/// and tests compile. Real body is `//`-commented until ndarray dep is live. +/// +/// # Arguments +/// - `splat_count` — number of splats in the source arrays (for validation). +/// - `position` — flat f32 slice, layout `[x0,y0,z0, x1,y1,z1, ...]`; must have length `splat_count * 3`. +/// - `rotation_xyzw` — flat f32 slice, layout `[x0,y0,z0,w0, x1,...]`; must have length `splat_count * 4`. +/// **Note:** glTF stores (x,y,z,w); the bridge reorders to GaussianBatch's (w,x,y,z). +/// - `scale` — flat f32 slice, layout `[sx0,sy0,sz0, ...]`; must have length `splat_count * 3`. +/// - `opacity` — flat f32 slice, one value per splat; must have length `splat_count`. +/// - `sh_coefficients` — optional flat f32 slice; `None` means no SH. +/// Layout per splat: `[k=0 rgb, k=1 rgb, ..., k=15 rgb]` (48 f32 for degree-3); +/// when `Some`, must have length `splat_count * 48`. +/// +/// # Returns +/// `Ok(CamSoaResult)` on success; `Err(ToCamSoaError)` on validation failure. +/// +/// # Errors +/// - [`ToCamSoaError::EmptySplats`] — `splat_count == 0`. +/// - [`ToCamSoaError::BufferLengthMismatch`] — any per-splat buffer has the wrong length. +/// - [`ToCamSoaError::ShLengthMismatch`] — `sh_coefficients` is `Some` but has the wrong length. +/// +/// # Examples +/// +/// ``` +/// use cesium::to_cam_soa::{convert_to_cam_soa, ToCamSoaError}; +/// +/// // Happy path: one splat, all buffers correctly sized. +/// let result = convert_to_cam_soa( +/// 1, +/// &[0.0, 0.0, 0.0], // position: 1*3 = 3 elements +/// &[0.0, 0.0, 0.0, 1.0], // rotation_xyzw: 1*4 = 4 elements (identity quat) +/// &[1.0, 1.0, 1.0], // scale: 1*3 = 3 elements +/// &[0.5], // opacity: 1 element +/// None, +/// ) +/// .unwrap(); +/// assert_eq!(result.splat_count, 1); +/// +/// // Error path: position buffer too short. +/// let err = convert_to_cam_soa( +/// 2, +/// &[0.0, 0.0, 0.0], // wrong: expected 6 (= 2*3) +/// &[0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0], +/// &[1.0, 1.0, 1.0, 1.0, 1.0, 1.0], +/// &[0.5, 0.5], +/// None, +/// ) +/// .unwrap_err(); +/// assert!(matches!( +/// err, +/// ToCamSoaError::BufferLengthMismatch { field: "position", expected: 6, actual: 3 } +/// )); +/// ``` +pub fn convert_to_cam_soa( + splat_count: usize, position: &[f32], rotation_xyzw: &[f32], scale: &[f32], opacity: &[f32], + sh_coefficients: Option<&[f32]>, +) -> Result { + if splat_count == 0 { + return Err(ToCamSoaError::EmptySplats); + } + // Validate all per-splat buffer lengths before proceeding. + let checks: &[(&'static str, usize, usize)] = &[ + ("position", splat_count * 3, position.len()), + ("rotation_xyzw", splat_count * 4, rotation_xyzw.len()), + ("scale", splat_count * 3, scale.len()), + ("opacity", splat_count, opacity.len()), + ]; + for &(field, expected, actual) in checks { + if actual != expected { + return Err(ToCamSoaError::BufferLengthMismatch { + field, + expected, + actual, + }); + } + } + // Validate SH length when present. + if let Some(sh) = sh_coefficients { + const SH_STRIDE: usize = 48; // 3 channels × 16 basis functions + let expected = splat_count * SH_STRIDE; + if sh.len() != expected { + return Err(ToCamSoaError::ShLengthMismatch { + splat_count, + sh_len: sh.len(), + expected_stride: SH_STRIDE, + }); + } + } + Ok(CamSoaResult { + splat_count, + _batch_placeholder: (), + _cam_fingerprints_placeholder: None, + }) +} + +// ───────────────────────────────────────────────────────────────────────────── +// COMMENTED-OUT full implementation (ndarray dep required) +// ───────────────────────────────────────────────────────────────────────────── +// +// Uncomment when `ndarray = { workspace = true }` is re-enabled in Cargo.toml +// and this module has been reviewed live by Opus + CodeRabbit. +// +// ```rust +// use ndarray_crate_alias::hpc::splat3d::gaussian::{GaussianBatch, SH_COEFFS_PER_GAUSSIAN}; +// use ndarray_crate_alias::hpc::cam_pq::{CamCodebook, CamFingerprint}; +// +// /// Full conversion: GsSplats → GaussianBatch + optional CAM fingerprints. +// /// +// /// # Quaternion reorder (CRITICAL) +// /// glTF / KHR_gaussian_splatting stores ROTATION as VEC4 (x,y,z,w). +// /// GaussianBatch stores quaternions as (w,x,y,z) — the `quat` field of +// /// `Gaussian3D` is documented as `[w, x, y, z]` (gaussian.rs line 53). +// /// The mapping is therefore: +// /// rotation_xyzw[i*4 + 3] → quat_w[i] ← w is at index 3 in glTF +// /// rotation_xyzw[i*4 + 0] → quat_x[i] +// /// rotation_xyzw[i*4 + 1] → quat_y[i] +// /// rotation_xyzw[i*4 + 2] → quat_z[i] +// /// +// /// # SH repack +// /// khr_gs layout: per-splat block of (n_coefs × 3) f32s, interleaved rgb per basis: +// /// gs_sh_block[k * 3 + ch] where k ∈ 0..16, ch ∈ {0=R,1=G,2=B} +// /// GaussianBatch layout: sh[i*48 + ch*16 + k] +// /// Repack: batch.sh[i*48 + ch*16 + k] = gs_sh_block[k*3 + ch] +// pub fn convert_to_cam_soa_live( +// splat_count: usize, +// position: &[f32], +// rotation_xyzw: &[f32], +// scale: &[f32], +// opacity: &[f32], +// sh_coefficients: Option<&[f32]>, +// codebook: Option<&CamCodebook>, +// ) -> Result<(GaussianBatch, Option>), ToCamSoaError> { +// if splat_count == 0 { +// return Err(ToCamSoaError::EmptySplats); +// } +// +// let mut batch = GaussianBatch::with_capacity(splat_count); +// +// // ── 1. Scatter SoA fields ────────────────────────────────────────────── +// for i in 0..splat_count { +// // Position → mean_x / mean_y / mean_z +// batch.mean_x[i] = position[i * 3 + 0]; +// batch.mean_y[i] = position[i * 3 + 1]; +// batch.mean_z[i] = position[i * 3 + 2]; +// +// // Scale → scale_x / scale_y / scale_z +// batch.scale_x[i] = scale[i * 3 + 0]; +// batch.scale_y[i] = scale[i * 3 + 1]; +// batch.scale_z[i] = scale[i * 3 + 2]; +// +// // Rotation reorder: glTF (x,y,z,w) → GaussianBatch (w,x,y,z) +// batch.quat_w[i] = rotation_xyzw[i * 4 + 3]; // *** w is at [3] in glTF *** +// batch.quat_x[i] = rotation_xyzw[i * 4 + 0]; +// batch.quat_y[i] = rotation_xyzw[i * 4 + 1]; +// batch.quat_z[i] = rotation_xyzw[i * 4 + 2]; +// +// // Opacity +// batch.opacity[i] = opacity[i]; +// } +// batch.len = splat_count; +// +// // ── 2. SH repack (optional) ──────────────────────────────────────────── +// if let Some(sh) = sh_coefficients { +// const SH_STRIDE: usize = SH_COEFFS_PER_GAUSSIAN; // 48 +// let expected = splat_count * SH_STRIDE; +// if sh.len() != expected { +// return Err(ToCamSoaError::ShLengthMismatch { +// splat_count, +// sh_len: sh.len(), +// expected_stride: SH_STRIDE, +// }); +// } +// for i in 0..splat_count { +// let gs_base = i * SH_STRIDE; // start of splat i's khr_gs block +// let batch_base = i * SH_STRIDE; // start of splat i's batch block +// // khr_gs: gs_sh_block[k * 3 + ch] +// // batch: sh[batch_base + ch * 16 + k] +// for k in 0..16_usize { +// for ch in 0..3_usize { +// batch.sh[batch_base + ch * 16 + k] = sh[gs_base + k * 3 + ch]; +// } +// } +// } +// } +// +// // ── 3. CAM-PQ encode (optional) ─────────────────────────────────────── +// // Encode position (3D) padded to codebook's total_dim with zeros. +// // UNVERIFIED: whether position alone or a richer feature vector should be +// // encoded. Using position as a spatial index is the conservative choice. +// let cam_fingerprints = if let Some(cb) = codebook { +// let mut fps: Vec = Vec::with_capacity(splat_count); +// let mut padded = vec![0.0f32; cb.total_dim]; +// for i in 0..splat_count { +// padded[0] = batch.mean_x[i]; +// padded[1] = batch.mean_y[i]; +// if cb.total_dim > 2 { padded[2] = batch.mean_z[i]; } +// fps.push(cb.encode(&padded)); +// } +// Some(fps) +// } else { +// None +// }; +// +// Ok((batch, cam_fingerprints)) +// } +// ``` + +// ───────────────────────────────────────────────────────────────────────────── +// Unit tests (std-only, live) +// ───────────────────────────────────────────────────────────────────────────── + +#[cfg(test)] +mod tests { + use super::*; + + fn make_flat_pos(n: usize) -> Vec { + (0..n * 3).map(|i| i as f32 * 0.1).collect() + } + fn make_flat_rot_xyzw(n: usize) -> Vec { + // Identity quaternion in glTF xyzw = (0,0,0,1) + let mut v = vec![0.0f32; n * 4]; + for i in 0..n { + v[i * 4 + 3] = 1.0; // w at index 3 + } + v + } + fn make_flat_scale(n: usize) -> Vec { + vec![1.0f32; n * 3] + } + fn make_flat_opacity(n: usize) -> Vec { + vec![1.0f32; n] + } + + #[test] + fn zero_splats_returns_error() { + let err = convert_to_cam_soa(0, &[], &[], &[], &[], None).unwrap_err(); + assert_eq!(err, ToCamSoaError::EmptySplats); + } + + #[test] + fn stub_returns_splat_count() { + let n = 4; + let result = convert_to_cam_soa( + n, + &make_flat_pos(n), + &make_flat_rot_xyzw(n), + &make_flat_scale(n), + &make_flat_opacity(n), + None, + ) + .unwrap(); + assert_eq!(result.splat_count, n); + } + + #[test] + fn sh_wrong_length_returns_error() { + let n = 2; + // correct length = 2 * 48 = 96; pass 47 instead + let bad_sh = vec![0.0f32; 47]; + let err = convert_to_cam_soa( + n, + &make_flat_pos(n), + &make_flat_rot_xyzw(n), + &make_flat_scale(n), + &make_flat_opacity(n), + Some(&bad_sh), + ) + .unwrap_err(); + assert!(matches!(err, ToCamSoaError::ShLengthMismatch { .. })); + } + + #[test] + fn sh_correct_length_ok() { + let n = 3; + let sh = vec![0.0f32; n * 48]; // 3 × 48 = 144 + let result = convert_to_cam_soa( + n, + &make_flat_pos(n), + &make_flat_rot_xyzw(n), + &make_flat_scale(n), + &make_flat_opacity(n), + Some(&sh), + ) + .unwrap(); + assert_eq!(result.splat_count, n); + } + + /// Document the quaternion reorder contract so it is machine-checkable once + /// the live impl is uncommented. + /// + /// glTF ROTATION VEC4 = (x, y, z, w); index mapping to GaussianBatch: + /// rotation_xyzw[i*4 + 3] → quat_w[i] (w is at offset 3 in glTF) + /// rotation_xyzw[i*4 + 0] → quat_x[i] + /// rotation_xyzw[i*4 + 1] → quat_y[i] + /// rotation_xyzw[i*4 + 2] → quat_z[i] + /// + /// This test records the offset map and will fail if anyone changes it. + #[test] + fn quat_reorder_offset_map_is_documented() { + // Encode the expected mapping as a const array: [(gltf_offset, field_name_tag)] + // We use a simple u8 tag: 0=w 1=x 2=y 3=z + const GLTF_TO_BATCH: [(usize, u8); 4] = [ + (3, 0), // gltf[3] → w + (0, 1), // gltf[0] → x + (1, 2), // gltf[1] → y + (2, 3), // gltf[2] → z + ]; + assert_eq!(GLTF_TO_BATCH[0], (3, 0), "w must come from gltf offset 3"); + assert_eq!(GLTF_TO_BATCH[1], (0, 1), "x must come from gltf offset 0"); + assert_eq!(GLTF_TO_BATCH[2], (1, 2), "y must come from gltf offset 1"); + assert_eq!(GLTF_TO_BATCH[3], (2, 3), "z must come from gltf offset 2"); + } + + /// Document SH repack formula. + /// khr_gs: gs_block[k * 3 + ch] → batch: sh[i*48 + ch*16 + k] + #[test] + fn sh_repack_index_formula() { + let n_splats = 1_usize; + let sh_stride = 48_usize; // SH_COEFFS_PER_GAUSSIAN + let i = 0_usize; + // For k=0, ch=0 (R DC term): + let gs_idx = 0 * 3 + 0; // k=0, ch=0 + let batch_idx = i * sh_stride + 0 * 16 + 0; // ch=0, k=0 + assert_eq!(gs_idx, 0); + assert_eq!(batch_idx, 0); + // For k=15, ch=2 (B last basis): + let gs_idx2 = 15 * 3 + 2; + let batch_idx2 = i * sh_stride + 2 * 16 + 15; + assert_eq!(gs_idx2, 47); + assert_eq!(batch_idx2, 47); + // Confirm n_splats is used (suppress unused warning) + let _ = n_splats; + } +}