diff --git a/CHANGELOG.md b/CHANGELOG.md index 158a096..c6ec41c 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -3,6 +3,10 @@ This changelog summarizes meaningful work on Joy Engine, Joy Editor, FORM LAB, PARTICLE LAB and the engine examples. Entries are grouped by day, newest first, using America/Detroit dates. Work before 10-08-2026 happened inside [JoyGames](https://github.com/trentpolack/JoyGames); its [CHANGELOG](https://github.com/trentpolack/JoyGames/blob/develop/CHANGELOG.md) is the historical record for that period. ## 10-16-2026 +### SDF Collisions +- Added signed distance field colliders to `physics` as a prototype: `SdfShape` (analytical primitives, or any triangle mesh baked into a voxel distance grid that tolerates open or overlapping procedural meshes) and `PhysicsWorld::create_sdf_body`. Contacts follow Giles and Andrews' caching and importance-sampling method (Pacific Graphics 2026) through a chained narrow-phase dispatcher, so worlds without SDF bodies simulate exactly as before, and SDF simulations stay deterministic natively and in WebAssembly. See [SDF colliders](docs/rust-simulation.md#sdf-colliders). +- Added the SDF COLLISIONS example: FORM props piling up in an SDF bin beside the same props as bounding boxes, with live costs and a `--bench`/`?bench=1` benchmark. Web builds' `THIRD-PARTY-NOTICES.txt` now also carries the reference implementation's MIT license. + ### Rust Runtime Default - Retired the deprecated engine, editor, lab, and example runtime implementation and its renderer dependencies. Rust is the shared native/browser implementation; canonical npm scripts wrap Cargo and browser builds. - Replaced current runtime and contribution documentation with Rust guidance and generated rustdoc, retaining useful Node tooling, asset schemas, capture orchestration, and recorded compatibility fixtures. diff --git a/Cargo.toml b/Cargo.toml index fdb4736..4736a81 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -47,6 +47,10 @@ path = "examples/form-building/main.rs" name = "form-destruction" path = "examples/form-destruction/main.rs" +[[example]] +name = "sdf-collisions" +path = "examples/sdf-collisions/main.rs" + [[test]] name = "engine" path = "test/rust/main.rs" @@ -58,7 +62,10 @@ bytemuck = { version = "1.25.2", features = ["derive"] } serde = { version = "1.0.229", features = ["derive"] } # float_roundtrip keeps authored documents byte-stable through a parse/serialize round trip. -serde_json = { version = "1.0.154", features = ["preserve_order", "float_roundtrip"] } +serde_json = { version = "1.0.154", features = [ + "preserve_order", + "float_roundtrip", +] } indexmap = { version = "2.14.2", features = ["serde"] } thiserror = "2.0.21" log = "0.4.34" @@ -82,14 +89,25 @@ image-webp = "0.2.4" rapier3d = { version = "0.36.1", features = ["enhanced-determinism"] } # Rendering, windows, and UI (rendering and platform layers). These match joy-apps' pins. -wgpu = { version = "30.0.1", default-features = false, features = ["std", "wgsl", "webgpu", "webgl", "vulkan", "metal", "dx12", "gles"] } +wgpu = { version = "30.0.1", default-features = false, features = [ + "std", + "wgsl", + "webgpu", + "webgl", + "vulkan", + "metal", + "dx12", + "gles", +] } winit = "0.30.13" egui = "0.36.2" egui-wgpu = { version = "0.36.2", default-features = false } [target.'cfg(not(target_arch = "wasm32"))'.dependencies] # egui-winit's pinned release does not build for WebAssembly; browsers use `browser::egui_input` instead. -egui-winit = { version = "0.36.2", default-features = false, features = ["clipboard"] } +egui-winit = { version = "0.36.2", default-features = false, features = [ + "clipboard", +] } pollster = "0.4.0" env_logger = "0.11.11" image = { version = "0.25.10", default-features = false, features = ["png"] } @@ -101,7 +119,9 @@ cpal = "0.18.2" # Native file dialogs for authoring tools (platform layer, behind development::FileDialogs); the version joy-apps uses. # The desktop portal backend loads libdbus at run time (falling back to zenity), so Linux builds need no system # packages; Wayland parent-window support is left out. -rfd = { version = "0.17.2", default-features = false, features = ["xdg-portal"] } +rfd = { version = "0.17.2", default-features = false, features = [ + "xdg-portal", +] } [target.'cfg(not(target_arch = "wasm32"))'.dev-dependencies] @@ -113,10 +133,40 @@ wasm-bindgen = "0.2.106" wasm-bindgen-futures = "0.4.56" js-sys = "0.3.83" web-sys = { version = "0.3.83", features = [ - "AudioBuffer", "AudioBufferSourceNode", "AudioContext", "AudioContextOptions", "AudioContextState", "AudioDestinationNode", - "AudioNode", "AudioParam", "AudioScheduledSourceNode", "BaseAudioContext", "Blob", "BlobPropertyBag", "Clipboard", "ClipboardEvent", "CssStyleDeclaration", - "DataTransfer", "Document", "Element", "EventTarget", "File", "FileList", "GainNode", "Gamepad", "GamepadButton", "GamepadMappingType", "HtmlCanvasElement", "HtmlElement", "HtmlInputElement", "Location", - "Navigator", "Storage", "Url", "UrlSearchParams", "Window", + "AudioBuffer", + "AudioBufferSourceNode", + "AudioContext", + "AudioContextOptions", + "AudioContextState", + "AudioDestinationNode", + "AudioNode", + "AudioParam", + "AudioScheduledSourceNode", + "BaseAudioContext", + "Blob", + "BlobPropertyBag", + "Clipboard", + "ClipboardEvent", + "CssStyleDeclaration", + "DataTransfer", + "Document", + "Element", + "EventTarget", + "File", + "FileList", + "GainNode", + "Gamepad", + "GamepadButton", + "GamepadMappingType", + "HtmlCanvasElement", + "HtmlElement", + "HtmlInputElement", + "Location", + "Navigator", + "Storage", + "Url", + "UrlSearchParams", + "Window", ] } console_error_panic_hook = "0.1.7" console_log = "1.0.0" diff --git a/docs/gameplay.md b/docs/gameplay.md index d718724..4632484 100644 --- a/docs/gameplay.md +++ b/docs/gameplay.md @@ -22,6 +22,8 @@ Unknown authored components remain data; only registered runtime behavior runs. ## Physics and Destruction Rapier provides rigid bodies and configurable raycast vehicles through engine-owned types. Vehicle tires are suspension rays rather than solid rolling colliders. Keep barriers thick enough and speeds appropriate for the discrete solver. See [simulation and physics](rust-simulation.md) for determinism, ownership, and vehicle limits. +Bodies are boxes, spheres, or planes by default. `SdfShape` colliders give non-convex procedural geometry (FORM meshes baked to distance grids, or analytical shapes) its true silhouette in collision; see [SDF colliders](rust-simulation.md#sdf-colliders). + `ChunkStructure` batches damage and support evaluation into a single committed revision; the renderer and physics backend synchronize from that revision. Project-specific collapse, damage, and debris rules stay in the game. ## Authoring diff --git a/docs/rust-simulation.md b/docs/rust-simulation.md index a3c1670..3b09eec 100644 --- a/docs/rust-simulation.md +++ b/docs/rust-simulation.md @@ -6,6 +6,7 @@ The headless simulation layer of Joy Engine's Rust crate (`src/environment`, `sr | --- | --- | | `environment` | `EnvironmentOptions` (serde), `EnvironmentSystem`, `EnvironmentFrame` | | `physics` | `PhysicsWorld`, `BodyOptions`/`BODY_DEFAULTS`, `VehicleOptions`/`WheelOptions`/`VEHICLE_DEFAULTS`/`WHEEL_DEFAULTS`, `BodyHandle`, `VehicleHandle`, `RaycastHit` | +| `physics` (SDF colliders) | `SdfShape`, `SdfPrimitive`, `SdfBakeOptions`/`SDF_BAKE_DEFAULTS`, `SdfCollisionStatistics`, `PhysicsWorld::create_sdf_body` | | `physics` (backend contract) | `PhysicsBackend` (implemented by `PhysicsWorld`), `PhysicsResidency`, `PhysicsCollisionShape`, `PhysicsMutationBatch`, `PhysicsQueryHit`, `create_removal_mutation_batch` | | `physics::destruction` | `ChunkStructure`, `DestructibleChunkDefinition`, `RadialDamage`, `ChunkChangeSet`, `ChunkDamageState`, `ChunkRaycastHit`, `ChunkSupportOptions`, `ChunkSupportState`, `create_chunk_mesh_data`, `DestructionError` | | `gameplay` (documents) | `EntityDefinition`, `EntityTransform`, `LevelData`, `parse_level`/`serialize_level`, `LevelDocument`, `ObjectDocument`, `resolve_level_objects`, `asset_path` | @@ -40,6 +41,19 @@ Physics uses rapier3d 0.36 with `enhanced-determinism`; identical worlds stepped - Rolling/brake friction is divided by grounded-wheel count before clamping. Opposite side impulses are not applied to dynamic ground bodies. - Wheel ray queries see new/moved colliders after the next solver step. Discrete solver behavior and authored tuning can produce wheelies or tunneling; test driving interactively. +## SDF Colliders +`SdfShape` describes a solid by its signed distance: `SdfShape::primitive` takes an analytical box, sphere, capsule, cylinder, or torus, and `SdfShape::from_triangle_mesh` bakes any triangle mesh (FORM output included) into a voxel grid of distances and gradients, with no convex decomposition. Baking signs the grid by winding number along X, Y, and Z grid lines with a majority vote, so open seams, flipped winding, and overlapping parts still enclose a solid. A shape precomputes its bounds, unit-density mass properties, and high-curvature seed points once; it is an immutable shared handle, so every body built from it reuses the grid. `PhysicsWorld::create_sdf_body` takes the shape plus ordinary `BodyOptions` (mass, pose, materials, and damping apply; the primitive fields are ignored), and the shape's frame is the body frame. + +Contacts follow Giles and Andrews, *Real-Time Collision Handling for Signed Distance Fields via Caching and Importance Sampling* (Pacific Graphics 2026), adapted from its MIT-licensed reference implementation (see `src/physics/sdf/LICENSE-sdf-dcd-demo`; web builds include it in `THIRD-PARTY-NOTICES.txt`). For each pair with overlapping bounds, Nelder-Mead minimizes a smooth maximum of both distances, minus a repulsion from points already found, from a few seeds: last step's points (a temporal cache that also carries the solver's warm-start impulses) and a random high-curvature seed. A chained narrow-phase dispatcher handles pairs involving an SDF collider (against SDFs, boxes, spheres, and planes) and passes every other pair to rapier's default dispatcher unchanged, so worlds without SDF bodies simulate exactly as before. Each manifold point becomes its own rapier manifold, which rapier's contact clustering merges by normal. `PhysicsWorld::sdf_collision_statistics` reports the latest fixed step's pairs, contacts, iterations, and distance queries. + +The prototype departs from the paper where the engine needed it: +- Normals and depths are read one-sidedly at the optimized point's projections onto both surfaces, keeping the deeper reading, because the paper's "deeper field at the point" rule tilts resting contacts wherever an edge or corner meets a face. +- Cached points live in both bodies' frames and ride their poses, refine to a quarter of the tolerance, keep their position unless the refinement is deeper, and skip optimization when the bodies barely moved relative to each other. +- Pairs with SDF colliders opt out of rapier's contact recycling (through a no-op contact filter hook), which would otherwise freeze the cache's step-by-step refinement while bodies settle. +- Sampling uses fixed per-pair seeded streams and the solver's scalar math, so identical worlds produce identical results, natively and in WebAssembly. + +Limits: collision stays discrete, so fast or thin bodies can tunnel; SDF colliders ignore meshes, compounds, sensors, shape casts, and continuous collision detection; and several contact points spread across a small curved-on-flat overlap, which acts as rolling resistance (an SDF sphere launched at 4 m/s on a plane rolls at about 1.85 m/s after three seconds, where rapier's ball keeps 2.57 m/s). Cost is CPU work per overlapping pair. In the [SDF collisions example](../examples/sdf-collisions/README.md) benchmark, 50 non-convex props cost about 4 ms per fixed step on average and 8 ms once the pile is complete natively, and about twice that in WebAssembly, against 0.1 to 0.2 ms for the same props as boxes. Baked grids take about 0.4 to 2.6 MB at 40 to 48 cells per side. + ## Destruction `ChunkStructure` is solver-independent and uses `f64` arithmetic and portable distance helpers, preserving recorded deterministic results. Damage is queued and applied atomically at `commit()`: overlapping hits accumulate, chunks at zero health leave first, then support is evaluated once and unsupported chunks detach in the same revision (the support diagnostics describe that evaluation). Support is a geometric masonry rule: chunks touching the ground plane anchor, stacked contacts and side-face contacts (never corners) carry support, and horizontal offsets spend the cantilever allowance; a shortest-path search with generation-order tie-breaking resolves reach. Revisions increment only when chunks leave and on every `reset()`. `create_chunk_mesh_data` returns closed boxes as plain `MeshAssetData` (positions, indices, per-vertex material tints, and `chunkCount`/`closedInteriorSurfaces` metadata) for the rendering layer to upload. diff --git a/examples/README.md b/examples/README.md index 1524f5c..1d11b7e 100644 --- a/examples/README.md +++ b/examples/README.md @@ -7,6 +7,7 @@ Small feature demonstrations live beside the engine and consume its public `joy- | [Scene Effects](scene-effects/README.md) | Lighting, shadows, and post-processing inspection | `/examples/scene-effects/` | | [FORM Building](form-building/README.md) | Editable procedural recipes and baked building assets | `/examples/form-building/` | | [FORM Destruction](form-destruction/README.md) | Surface-projected impacts and synchronized chunk removal | `/examples/form-destruction/` | +| [SDF Collisions](sdf-collisions/README.md) | Signed distance colliders for FORM props beside bounding boxes, with a native and browser benchmark | Local only (no `gameSite` metadata yet) | ## Run From the JoyGames root, replace `form-destruction` with any example folder name: diff --git a/examples/sdf-collisions/README.md b/examples/sdf-collisions/README.md new file mode 100644 index 0000000..aaddc8a --- /dev/null +++ b/examples/sdf-collisions/README.md @@ -0,0 +1,23 @@ +# SDF Collisions +A Rust example that compares signed distance field (SDF) colliders with bounding boxes. Two bins fill with the same FORM props, dropped in the same poses on the same schedule: the left bin collides through each prop's SDF, the right through the box the engine could offer before. Chairs hook legs, jacks interlock, and staples nest on the left; on the right the same props stack as crates. + +```bash +npm run dev --workspace @joy-engine/example-sdf-collisions +npm run build --workspace @joy-engine/example-sdf-collisions +npm run test --workspace @joy-engine/example-sdf-collisions +npm run bench --workspace @joy-engine/example-sdf-collisions +npm run capture --workspace @joy-engine/example-sdf-collisions -- --built +``` + +Drag to orbit, scroll to zoom, Space or DROP MORE to add ten props per bin, P to pause, and R to reset. The sidebar shows each bin's CPU physics time per fixed step, the SDF contact work of the latest step, and each collider's build cost and memory. + +## Props and Colliders +The chair, jack, staple, and I-beam are FORM scripts in `assets/` built from axis-aligned boxes. Each evaluates once, draws through the engine's PBR path, and bakes into a 40-cell voxel distance grid with `SdfShape::from_triangle_mesh`; every body shares its prop's grid. The ring uses the analytical `SdfPrimitive::Torus`, and its FORM mesh only draws it. The box bin uses each collider's local bounds, and both bins give a prop the same mass (its SDF volume times a pine-like density). + +Box-only props keep the baked colliders identical on every target: FORM's spheres and rotated instances use platform trigonometry, which can differ in the last bit between native and WebAssembly builds. + +## Benchmark +`--bench` (native) runs a headless scenario, prints one line, and exits; `?bench=1` runs it in browsers before the demo starts, logs it to the console, and shows it in the sidebar. Both bins fill with 50 props over 1,200 fixed steps. The report gives CPU milliseconds per fixed step for the whole run and for the final, settled two seconds, plus a hash of every body's final position. Matching SDF hashes between a native run and a browser run show that SDF contact generation is deterministic across targets. + +## Verification +`test/rust/example_sdf_collisions.rs` checks that the props bake into colliders that fit their meshes, and that the benchmark piles repeat bit for bit with no SDF prop leaving its bin. `test/rust/physics.rs` covers SDF resting contact, stacking, raycasts, validation, mesh baking, and determinism in the engine itself. Browser captures add canvas startup and image evidence. CPU timing does not measure GPU execution. diff --git a/examples/sdf-collisions/assets/beam.form b/examples/sdf-collisions/assets/beam.form new file mode 100644 index 0000000..ca51262 --- /dev/null +++ b/examples/sdf-collisions/assets/beam.form @@ -0,0 +1,18 @@ +# I-BEAM: two flanges joined by a web. Flanges nest into each other's channels when beams stack. +project meta(name="I-beam", info="SDF COLLISIONS / PROP") + +param length = 1.1 meta(min=0.4, max=2, step=0.05) +param flange = 0.42 meta(min=0.2, max=0.8, step=0.02) +param depth = 0.42 meta(min=0.2, max=0.8, step=0.02) +param plate = 0.07 meta(min=0.03, max=0.15, step=0.01) + +# White vertices, so the material's base color shows as authored. +color(1, 1, 1) +let unit = meshBox(1, 1, 1) +let finish = materialPbr(rgba(0.5, 0.64, 0.9, 1), 0.15, 0.45) + +repeat 2 as side { + let y = (side*2 - 1)*(depth/2 - plate/2) + instance(unit, vector(0, y, 0), vector(0, 0, 0), vector(length, plate, flange), finish) +} +instance(unit, vector(0, 0, 0), vector(0, 0, 0), vector(length, depth - plate*2, plate), finish) diff --git a/examples/sdf-collisions/assets/chair.form b/examples/sdf-collisions/assets/chair.form new file mode 100644 index 0000000..e9ea0a6 --- /dev/null +++ b/examples/sdf-collisions/assets/chair.form @@ -0,0 +1,25 @@ +# CHAIR: a seat, four legs, and a back rest. Its legs and back hook into neighbors, the classic case for +# convex decomposition and the paper's "pile of chairs". Axis-aligned boxes only, so the baked collider is +# identical on every platform. +project meta(name="Chair", info="SDF COLLISIONS / PROP") + +param seat = 0.7 meta(min=0.4, max=1.2, step=0.05) +param leg = 0.5 meta(min=0.2, max=1, step=0.05) +param back = 0.6 meta(min=0.2, max=1, step=0.05) +param thickness = 0.08 meta(min=0.04, max=0.2, step=0.01) + +# White vertices, so the material's base color shows as authored. +color(1, 1, 1) +let unit = meshBox(1, 1, 1) +let finish = materialPbr(rgba(0.86, 0.5, 0.26, 1), 0, 0.62) +let inset = seat/2 - thickness/2 + +instance(unit, vector(0, 0, 0), vector(0, 0, 0), vector(seat, thickness, seat), finish) +repeat 2 as column { + repeat 2 as row { + let x = (column*2 - 1)*inset + let z = (row*2 - 1)*inset + instance(unit, vector(x, -(leg + thickness)/2, z), vector(0, 0, 0), vector(thickness, leg, thickness), finish) + } +} +instance(unit, vector(0, (back + thickness)/2, -inset), vector(0, 0, 0), vector(seat, back, thickness), finish) diff --git a/examples/sdf-collisions/assets/jack.form b/examples/sdf-collisions/assets/jack.form new file mode 100644 index 0000000..1ff3bd0 --- /dev/null +++ b/examples/sdf-collisions/assets/jack.form @@ -0,0 +1,22 @@ +# JACK: three crossed bars with blocky end caps. Every orientation lands on a different set of arms. +project meta(name="Jack", info="SDF COLLISIONS / PROP") + +param span = 1 meta(min=0.4, max=1.6, step=0.05) +param bar = 0.14 meta(min=0.06, max=0.3, step=0.01) +param cap = 0.24 meta(min=0.1, max=0.4, step=0.01) + +# White vertices, so the material's base color shows as authored. +color(1, 1, 1) +let unit = meshBox(1, 1, 1) +let finish = materialPbr(rgba(0.3, 0.86, 0.82, 1), 0.05, 0.45) +let reach = span/2 - cap/2 + +instance(unit, vector(0, 0, 0), vector(0, 0, 0), vector(span - cap, bar, bar), finish) +instance(unit, vector(0, 0, 0), vector(0, 0, 0), vector(bar, span - cap, bar), finish) +instance(unit, vector(0, 0, 0), vector(0, 0, 0), vector(bar, bar, span - cap), finish) +repeat 2 as side { + let sign = side*2 - 1 + instance(unit, vector(sign*reach, 0, 0), vector(0, 0, 0), vector(cap, cap, cap), finish) + instance(unit, vector(0, sign*reach, 0), vector(0, 0, 0), vector(cap, cap, cap), finish) + instance(unit, vector(0, 0, sign*reach), vector(0, 0, 0), vector(cap, cap, cap), finish) +} diff --git a/examples/sdf-collisions/assets/ring.form b/examples/sdf-collisions/assets/ring.form new file mode 100644 index 0000000..79b0d55 --- /dev/null +++ b/examples/sdf-collisions/assets/ring.form @@ -0,0 +1,16 @@ +# RING: the render mesh of the analytical torus collider (major radius 0.36, minor radius 0.1). Its collider is +# an exact distance function, so this mesh never feeds physics. +project meta(name="Ring", info="SDF COLLISIONS / PROP") + +param major = 0.36 meta(min=0.2, max=0.8, step=0.01) +param minor = 0.1 meta(min=0.03, max=0.3, step=0.01) +param detail = 40 meta(min=12, max=64, step=4) + +let profile = geometry() +repeat 17 as i { + let angle = i/16*6.283185307179586 + addPoint(profile, vector(major + cos(angle)*minor, sin(angle)*minor, 0)) +} +let ring = revolve(profile, detail) +color(ring, rgba(0.86, 0.3, 0.44, 1)) +output("Ring", smoothNormals(ring)) diff --git a/examples/sdf-collisions/assets/staple.form b/examples/sdf-collisions/assets/staple.form new file mode 100644 index 0000000..5d36df3 --- /dev/null +++ b/examples/sdf-collisions/assets/staple.form @@ -0,0 +1,17 @@ +# STAPLE: two posts joined by a bridge. Open shapes like this interlock and wedge into each other. +project meta(name="Staple", info="SDF COLLISIONS / PROP") + +param width = 0.8 meta(min=0.3, max=1.4, step=0.05) +param height = 0.8 meta(min=0.3, max=1.4, step=0.05) +param bar = 0.14 meta(min=0.06, max=0.3, step=0.01) + +# White vertices, so the material's base color shows as authored. +color(1, 1, 1) +let unit = meshBox(1, 1, 1) +let finish = materialPbr(rgba(0.93, 0.76, 0.23, 1), 0.1, 0.5) +let post = width/2 - bar/2 + +repeat 2 as side { + instance(unit, vector((side*2 - 1)*post, -bar/2, 0), vector(0, 0, 0), vector(bar, height - bar, bar), finish) +} +instance(unit, vector(0, height/2 - bar/2, 0), vector(0, 0, 0), vector(width, bar, bar), finish) diff --git a/examples/sdf-collisions/main.rs b/examples/sdf-collisions/main.rs new file mode 100644 index 0000000..37aaa4a --- /dev/null +++ b/examples/sdf-collisions/main.rs @@ -0,0 +1,57 @@ +// Copyright (c) 2026 Trent Polack. All Rights Reserved. +// Licensed under the MIT License. + +//! SDF COLLISIONS: the same FORM props dropped into two bins, one colliding through signed distance fields and one +//! through bounding boxes, to compare how the piles settle and what each costs. +//! +//! cargo run --example sdf-collisions +//! cargo run --release --example sdf-collisions -- --bench +//! node scripts/rust/build-web.mjs --example sdf-collisions +//! +//! Drag to orbit, scroll to zoom, Space to drop more props, P to pause, and R to reset. Natively, `--bench` prints the +//! headless benchmark and exits; in browsers, `?bench=1` logs it to the console and shows it in the sidebar. + +#[path = "src/mod.rs"] +mod example; + +use example::sdf_collisions_demo::{SdfCollisionsDemo, input_map}; +use joy_engine::prelude::*; + +/// Shell settings: default renderer with drag-to-orbit and the demo's keys. +fn options() -> ShellOptions { + ShellOptions { title: "SDF COLLISIONS".into(), input: input_map(), ..Default::default() } +} + +fn start() { + let result = SdfCollisionsDemo::new() + .map_err(|error| ShellError::Game(error.to_string())) + .and_then(|demo| run(options(), demo)); + if let Err(error) = result { + joy_engine::log::error!("{error}"); + } +} + +/// Print the headless benchmark without opening a window. +#[cfg(not(target_arch = "wasm32"))] +fn print_benchmark() { + use example::prop_library::PropLibrary; + use example::sdf_collisions_demo::benchmark; + + let report = PropLibrary::load().and_then(|library| benchmark(&library)); + match report { + Ok(report) => println!("{report}"), + Err(error) => eprintln!("{error}"), + } +} + +joy_engine::web_main!(start); + +fn main() { + // Browsers start through `web_main!` instead, and run the benchmark inside the demo. + #[cfg(not(target_arch = "wasm32"))] + if example::sdf_collisions_demo::benchmark_requested() { + print_benchmark(); + } else { + start(); + } +} diff --git a/examples/sdf-collisions/package.json b/examples/sdf-collisions/package.json new file mode 100644 index 0000000..dda34f6 --- /dev/null +++ b/examples/sdf-collisions/package.json @@ -0,0 +1,20 @@ +{ + "name": "@joy-engine/example-sdf-collisions", + "version": "0.1.0", + "private": true, + "type": "module", + "author": "Trent Polack", + "license": "MIT", + "scripts": { + "dev": "cargo run --package joy-engine --example sdf-collisions", + "build": "node ../../scripts/rust/build-web.mjs --example sdf-collisions --release --title \"SDF COLLISIONS · Joy Engine\" --out dist", + "preview": "node ../../scripts/preview-web.mjs", + "check": "cargo check --package joy-engine --example sdf-collisions", + "test": "cargo test --package joy-engine example_sdf_collisions", + "bench": "cargo run --release --package joy-engine --example sdf-collisions -- --bench", + "capture": "node ../../scripts/capture-web.mjs" + }, + "dependencies": { + "joy-engine": "file:../.." + } +} diff --git a/examples/sdf-collisions/src/benchmark.rs b/examples/sdf-collisions/src/benchmark.rs new file mode 100644 index 0000000..67edc7a --- /dev/null +++ b/examples/sdf-collisions/src/benchmark.rs @@ -0,0 +1,101 @@ +// Copyright (c) 2026 Trent Polack. All Rights Reserved. +// Licensed under the MIT License. + +//! The headless benchmark: both bins fill with the same props on the same schedule, without a renderer, and report +//! their CPU physics time per fixed step and a hash of the final pile. +//! +//! Natively `--bench` prints the report and exits; in browsers `?bench=1` logs it to the console and shows it in the +//! HUD. Comparing the two hashes across targets checks that SDF contact generation is deterministic there too. + +use std::fmt; + +use joy_engine::physics::{PhysicsError, SdfCollisionStatistics}; +use joy_engine::runtime::Stopwatch; + +use super::config::sdf_collisions_config::PHYSICS; +use super::pile_arena::{ArenaSnapshot, ColliderMode, PileArena}; +use super::prop_library::PropLibrary; + +/// Final steps averaged for the settled-pile timing. +const SETTLED_STEPS: u32 = 120; + +/// One bin's benchmark result. +#[derive(Clone, Debug, PartialEq)] +pub struct BinReport { + /// CPU milliseconds per fixed step, averaged over the whole run (drops, impacts, and settling). + pub milliseconds_per_step: f64, + /// CPU milliseconds per fixed step over the final steps, once the pile has (mostly) come to rest. + pub settled_milliseconds_per_step: f64, + /// Position hash of the final pile. + pub hash: u64, + /// The final pile's diagnostics. + pub snapshot: ArenaSnapshot, +} + +/// Both bins' results for the same scenario. +#[derive(Clone, Debug, PartialEq)] +pub struct BenchmarkReport { + /// Fixed steps simulated per bin. + pub steps: u32, + /// Props dropped into each bin. + pub props: usize, + /// The SDF bin and the bounding-box bin. + pub sdf: BinReport, + pub boxes: BinReport, +} + +/// Fill both bins with `props` props over `steps` fixed steps and time them separately. +pub fn run_benchmark(library: &PropLibrary, steps: u32, props: usize) -> Result { + Ok(BenchmarkReport { + steps, + props, + sdf: run_bin(ColliderMode::Sdf, library, steps, props)?, + boxes: run_bin(ColliderMode::BoundingBox, library, steps, props)?, + }) +} + +fn run_bin(mode: ColliderMode, library: &PropLibrary, steps: u32, props: usize) -> Result { + let mut arena = PileArena::new(mode, 0.0, None)?; + arena.set_target(props); + + let settled_from = steps.saturating_sub(SETTLED_STEPS); + let mut total_milliseconds = 0.0; + let mut settled_milliseconds = 0.0; + for step in 0..steps { + let stopwatch = Stopwatch::start(); + arena.update(PHYSICS.fixed_step, library, None)?; + let elapsed = stopwatch.elapsed_milliseconds(); + total_milliseconds += elapsed; + if step >= settled_from { + settled_milliseconds += elapsed; + } + } + Ok(BinReport { + milliseconds_per_step: total_milliseconds / f64::from(steps.max(1)), + settled_milliseconds_per_step: settled_milliseconds / f64::from((steps - settled_from).max(1)), + hash: arena.position_hash(), + snapshot: arena.snapshot(), + }) +} + +impl fmt::Display for BenchmarkReport { + fn fmt(&self, formatter: &mut fmt::Formatter<'_>) -> fmt::Result { + let SdfCollisionStatistics { pairs, contacts, iterations, distance_queries } = self.sdf.snapshot.statistics; + write!( + formatter, + "SDF benchmark: {} props, {} steps | SDF {:.3} ms/step ({:.3} settled), hash {:016x}, {} left the bin | \ + boxes {:.3} ms/step ({:.3} settled), hash {:016x}, {} left the bin | last step: {pairs} pairs, \ + {contacts} contacts, {iterations} iterations, {distance_queries} distance queries", + self.props, + self.steps, + self.sdf.milliseconds_per_step, + self.sdf.settled_milliseconds_per_step, + self.sdf.hash, + self.sdf.snapshot.escaped, + self.boxes.milliseconds_per_step, + self.boxes.settled_milliseconds_per_step, + self.boxes.hash, + self.boxes.snapshot.escaped, + ) + } +} diff --git a/examples/sdf-collisions/src/config/sdf_collisions_config.rs b/examples/sdf-collisions/src/config/sdf_collisions_config.rs new file mode 100644 index 0000000..7d0643b --- /dev/null +++ b/examples/sdf-collisions/src/config/sdf_collisions_config.rs @@ -0,0 +1,173 @@ +// Copyright (c) 2026 Trent Polack. All Rights Reserved. +// Licensed under the MIT License. + +//! Named tuning for the SDF collisions example: arena layout, the shared spawn sequence, physics materials, the bake +//! resolution, the benchmark scenario, the camera, and presentation. Lengths are meters, times seconds, masses +//! kilograms, and angles radians. + +use joy_engine::physics::SdfBakeOptions; +use joy_engine::prelude::*; + +/// Two identical bins side by side: SDF colliders on the left, the props' bounding boxes on the right. +pub struct ArenaSettings { + /// Distance between the two bins' centers along X. + pub separation: f32, + /// Half the inner width of a bin along X and Z. + pub inner_half_extent: f32, + /// Wall height above the bin floor. + pub wall_height: f32, + /// Wall and floor slab thickness. + pub wall_thickness: f32, + /// Bodies farther than this below the floor, or outside the bin by this much, count as escaped (tunneled). + pub escape_margin: f32, +} + +/// The spawn sequence. Both bins draw it from identically seeded streams, so they receive the same props in the same +/// poses at the same times. +pub struct SpawnSettings { + /// Seed of both bins' spawn streams. + pub seed: i64, + /// Seconds between drops. + pub interval_seconds: f32, + /// Drop height above the bin floor. + pub height: f32, + /// Half the width of the square drop area, centered on the bin. + pub spread: f32, + /// Props per bin before the pile is complete. + pub initial_count: usize, + /// Props added by DROP MORE. + pub burst: usize, + /// Most props a bin holds. + pub maximum_count: usize, +} + +/// Rigid-body materials and stepping. +pub struct PhysicsSettings { + /// Fixed simulation step in seconds (the engine default). + pub fixed_step: f32, + /// Most fixed steps per frame. + pub max_sub_steps: u32, + /// Prop density in kilograms per cubic meter (about pine); mass is the SDF volume times this in both bins. + pub density: f32, + /// Friction of props, floors, and walls. + pub friction: f32, + /// Restitution of props. + pub restitution: f32, + /// Per-second linear and angular velocity loss fractions (engine semantics). + pub linear_damping: f32, + pub angular_damping: f32, + /// Distance-grid resolution for meshes baked from FORM. + pub bake: SdfBakeOptions, +} + +/// The headless benchmark behind `--bench` (native) and `?bench=1` (browsers). +pub struct BenchmarkSettings { + /// Launch parameter that runs it. + pub parameter: &'static str, + /// Fixed steps simulated per bin. + pub steps: u32, + /// Props dropped into each bin. + pub props: usize, +} + +/// Orbit camera framing and controls. +pub struct ViewSettings { + /// Orbit center, between the two bins. + pub target: Vec3, + /// Starting distance from the target. + pub distance: f32, + /// Starting heading and elevation. + pub yaw: f32, + pub pitch: f32, + /// Near and far clip distances. + pub clips: (f32, f32), + /// Elevation and zoom limits. + pub pitch_limits: (f32, f32), + pub distance_limits: (f32, f32), + /// Radians of orbit per logical pixel dragged, and the exponential zoom rate per scroll unit. + pub orbit_sensitivity: f32, + pub zoom_sensitivity: f32, +} + +/// Lights, colors, and timing. +pub struct PresentationSettings { + /// Ambient light strength and the background color. + pub ambient: f32, + pub clear_color: [f32; 4], + /// A warm shadowed key light above the front-left, and a cool fill from behind. + pub key_light: PointLight, + pub fill_light: PointLight, + /// Bin floor tints: SDF bin, then box bin. The HUD legend uses the same colors. + pub sdf_tint: Rgb, + pub box_tint: Rgb, + pub wall_tint: Rgb, + /// Timing samples kept per distribution. + pub timing_samples: usize, + /// Longest frame the simulation advances, in seconds; hidden tabs and stalls never fast-forward the piles. + pub maximum_frame_seconds: f32, +} + +pub const ARENA: ArenaSettings = ArenaSettings { + separation: 6.2, + inner_half_extent: 1.9, + wall_height: 1.0, + wall_thickness: 0.2, + escape_margin: 1.0, +}; + +pub const SPAWN: SpawnSettings = SpawnSettings { + seed: 2026, + interval_seconds: 0.3, + height: 4.2, + spread: 1.2, + initial_count: 30, + burst: 10, + maximum_count: 90, +}; + +pub const PHYSICS: PhysicsSettings = PhysicsSettings { + fixed_step: 1.0 / 60.0, + max_sub_steps: 4, + density: 450.0, + friction: 0.6, + restitution: 0.05, + linear_damping: 0.05, + angular_damping: 0.1, + bake: SdfBakeOptions { resolution: 40 }, +}; + +pub const BENCHMARK: BenchmarkSettings = BenchmarkSettings { parameter: "bench", steps: 1200, props: 50 }; + +pub const VIEW: ViewSettings = ViewSettings { + target: Vec3::new(0.0, 0.6, 0.0), + distance: 12.5, + yaw: 0.0, + pitch: 0.8, + clips: (0.1, 200.0), + pitch_limits: (0.05, 1.45), + distance_limits: (3.0, 60.0), + orbit_sensitivity: 0.005, + zoom_sensitivity: 0.001, +}; + +pub const PRESENTATION: PresentationSettings = PresentationSettings { + ambient: 0.75, + clear_color: [0.055, 0.075, 0.095, 1.0], + key_light: PointLight { + position: Vec3::new(-4.0, 11.0, 7.0), + color: [8.0, 7.4, 6.6], + range: 40.0, + casts_shadow: true, + }, + fill_light: PointLight { + position: Vec3::new(6.0, 7.0, -6.0), + color: [2.0, 2.4, 3.0], + range: 30.0, + casts_shadow: false, + }, + sdf_tint: [0.27, 0.62, 0.58], + box_tint: [0.74, 0.42, 0.27], + wall_tint: [0.32, 0.35, 0.38], + timing_samples: 120, + maximum_frame_seconds: 0.1, +}; diff --git a/examples/sdf-collisions/src/hud.rs b/examples/sdf-collisions/src/hud.rs new file mode 100644 index 0000000..f4eaf17 --- /dev/null +++ b/examples/sdf-collisions/src/hud.rs @@ -0,0 +1,198 @@ +// Copyright (c) 2026 Trent Polack. All Rights Reserved. +// Licensed under the MIT License. + +//! The example's sidebar: what the two bins compare, live pile and contact diagnostics per bin, the colliders' build +//! costs, the benchmark result, and DROP MORE, PAUSE, and RESET. Drawing returns the requested changes instead of +//! mutating the demo. + +use joy_engine::development::dock; +use joy_engine::egui::{Align2, Color32, FontId, Id, Margin, RichText, Stroke, Vec2 as UiVec2}; +use joy_engine::prelude::*; + +use super::config::sdf_collisions_config::PRESENTATION; +use super::pile_arena::ArenaSnapshot; +use super::prop_library::Prop; + +/// Sidebar width in logical pixels, and on viewports at most [`NARROW_VIEWPORT`] wide. +pub const SIDEBAR_WIDTH: f32 = 330.0; +const SIDEBAR_WIDTH_NARROW: f32 = 260.0; +const NARROW_VIEWPORT: f32 = 760.0; + +// The sidebar palette. +const FILL: Color32 = Color32::from_rgba_premultiplied(0x0d, 0x14, 0x1a, 0xf0); +const TITLE: Color32 = Color32::from_rgb(0x8f, 0xe3, 0xd2); +const TEXT: Color32 = Color32::from_rgb(0xe3, 0xea, 0xee); +const MUTED: Color32 = Color32::from_rgb(0x9d, 0xab, 0xb5); +const RULE: Color32 = Color32::from_rgb(0x2b, 0x36, 0x40); +const BUTTON: Color32 = Color32::from_rgb(0x1a, 0x25, 0x2e); +const BUTTON_LINE: Color32 = Color32::from_rgb(0x46, 0x56, 0x63); +const WARNING: Color32 = Color32::from_rgb(0xff, 0xb0, 0x8a); + +/// Fields smaller than this are analytical: they store dimensions, not a grid. +const ANALYTICAL_FIELD_BYTES: usize = 4096; + +/// A change requested from the sidebar. +#[derive(Clone, Copy, Debug, PartialEq, Eq)] +pub enum HudCommand { + /// Grow both piles by a burst of props. + DropMore, + /// Pause or resume both simulations. + TogglePause, + /// Empty both bins and start the drop sequence again. + Reset, +} + +/// Everything the sidebar shows, borrowed from the demo for one frame. +pub struct HudView<'a> { + pub sdf: &'a ArenaSnapshot, + pub boxes: &'a ArenaSnapshot, + pub props: &'a [Prop], + pub field_bytes: usize, + pub paused: bool, + pub backend: &'a str, + pub render_milliseconds: f64, + pub benchmark: Option<&'a str>, +} + +/// Sidebar width for a viewport width in logical pixels. +pub fn sidebar_width(viewport_width: f32) -> f32 { + if viewport_width <= NARROW_VIEWPORT { SIDEBAR_WIDTH_NARROW } else { SIDEBAR_WIDTH } +} + +/// Draw the sidebar on the right edge and return the commands it requested. +pub fn draw(context: &egui::Context, view: &HudView<'_>) -> Vec { + let mut commands = Vec::new(); + let screen = context.content_rect(); + let width = sidebar_width(screen.width()); + let margin: i8 = if width == SIDEBAR_WIDTH { 20 } else { 12 }; + let reserved = dock::reserved_height(context); + egui::Area::new(Id::new("sdf-collisions-sidebar")).anchor(Align2::RIGHT_TOP, UiVec2::ZERO).constrain(false).show( + context, + |ui| { + egui::Frame::new().fill(FILL).inner_margin(Margin::same(margin)).show(ui, |ui| { + let inner = UiVec2::new(width, screen.height()) - UiVec2::splat(f32::from(margin) * 2.0); + ui.set_width(inner.x); + ui.set_height(inner.y); + style(ui); + egui::ScrollArea::vertical().auto_shrink([false, false]).show(ui, |ui| { + contents(ui, view, &mut commands); + ui.add_space(reserved); + }); + }); + }, + ); + commands +} + +fn contents(ui: &mut egui::Ui, view: &HudView<'_>, commands: &mut Vec) { + ui.label(RichText::new("SDF COLLISIONS").font(FontId::proportional(20.0)).strong().color(TITLE)); + ui.label( + RichText::new("The same FORM props, dropped in the same poses. Left: signed distance colliders. Right: the boxes the engine could use before.") + .color(MUTED), + ); + ui.label( + RichText::new("Drag to orbit · scroll to zoom · Space drops more · P pauses · R resets").small().color(MUTED), + ); + + section(ui, "SDF COLLIDERS", PRESENTATION.sdf_tint); + bin_lines(ui, view.sdf); + let statistics = view.sdf.statistics; + ui.label(format!("{} pairs · {} contacts", statistics.pairs, statistics.contacts)); + ui.label(format!("{} iterations · {} distance queries", statistics.iterations, statistics.distance_queries)); + + section(ui, "BOUNDING-BOX COLLIDERS", PRESENTATION.box_tint); + bin_lines(ui, view.boxes); + + section(ui, "COLLIDERS", PRESENTATION.wall_tint); + ui.label( + RichText::new(format!("Distance grids share {:.1} MB across every body.", megabytes(view.field_bytes))) + .small() + .color(MUTED), + ); + for prop in view.props { + let field = if prop.field_bytes >= ANALYTICAL_FIELD_BYTES { + format!("{} triangles · {:.1} MB grid", prop.triangles, megabytes(prop.field_bytes)) + } else { + "analytical torus".to_owned() + }; + ui.label(format!( + "{}: {field} · {:.0} ms · {} seeds", + prop.name, + prop.collider_milliseconds, + prop.shape.seed_count() + )); + } + + ui.add_space(6.0); + let full = |ui: &mut egui::Ui, label: &str| { + ui.add(egui::Button::new(label).min_size(UiVec2::new(ui.available_width(), 28.0))).clicked() + }; + if full(ui, "DROP MORE") { + commands.push(HudCommand::DropMore); + } + if full(ui, if view.paused { "RESUME" } else { "PAUSE" }) { + commands.push(HudCommand::TogglePause); + } + if full(ui, "RESET") { + commands.push(HudCommand::Reset); + } + + ui.add_space(6.0); + ui.label( + RichText::new(format!("{} · render encode {:.2} ms (CPU)", view.backend, view.render_milliseconds)) + .small() + .color(MUTED), + ); + if let Some(benchmark) = view.benchmark { + ui.label(RichText::new(benchmark).small().monospace().color(TEXT)); + } + ui.label( + RichText::new("Contacts: Giles and Andrews, Real-Time Collision Handling for Signed Distance Fields via Caching and Importance Sampling (Pacific Graphics 2026).") + .small() + .color(MUTED), + ); +} + +/// Bytes to megabytes for display. +fn megabytes(bytes: usize) -> f64 { + bytes as f64 / (1024.0 * 1024.0) +} + +/// A section heading with a color swatch matching the bin floor. +fn section(ui: &mut egui::Ui, title: &str, tint: Rgb) { + ui.add_space(4.0); + ui.separator(); + ui.horizontal(|ui| { + let (rect, _) = ui.allocate_exact_size(UiVec2::splat(12.0), egui::Sense::hover()); + let [red, green, blue] = tint.map(|channel| (channel.clamp(0.0, 1.0) * 255.0) as u8); + ui.painter().rect_filled(rect, 2.0, Color32::from_rgb(red, green, blue)); + ui.label(RichText::new(title).strong().color(TEXT)); + }); +} + +/// Pile size, escapes, and physics time for one bin. +fn bin_lines(ui: &mut egui::Ui, snapshot: &ArenaSnapshot) { + ui.label(format!("{} / {} props", snapshot.bodies, snapshot.target)); + ui.label(format!( + "physics {:.2} ms/step · worst {:.2} ms", + snapshot.mean_step_milliseconds, snapshot.worst_step_milliseconds + )); + if snapshot.escaped > 0 { + ui.label(RichText::new(format!("{} escaped the bin", snapshot.escaped)).color(WARNING)); + } +} + +/// The sidebar's text, rule, and button look. +fn style(ui: &mut egui::Ui) { + let visuals = ui.visuals_mut(); + visuals.override_text_color = Some(TEXT); + visuals.widgets.noninteractive.bg_stroke = Stroke::new(1.0, RULE); + for widget in [&mut visuals.widgets.inactive, &mut visuals.widgets.hovered, &mut visuals.widgets.active] { + widget.weak_bg_fill = BUTTON; + widget.bg_fill = BUTTON; + widget.bg_stroke = Stroke::new(1.0, BUTTON_LINE); + widget.corner_radius = 4.0.into(); + } + visuals.widgets.hovered.bg_stroke = Stroke::new(1.0, TITLE); + ui.spacing_mut().item_spacing.y = 6.0; +} diff --git a/examples/sdf-collisions/src/mod.rs b/examples/sdf-collisions/src/mod.rs new file mode 100644 index 0000000..b9f8408 --- /dev/null +++ b/examples/sdf-collisions/src/mod.rs @@ -0,0 +1,16 @@ +// Copyright (c) 2026 Trent Polack. All Rights Reserved. +// Licensed under the MIT License. + +//! The SDF collisions example's modules: tuning (`config/`), the FORM props and their colliders, one bin and its pile, +//! the headless benchmark, the sidebar, and the demo itself. `main.rs` and the integration tests in `test/rust/` both +//! include this tree. + +/// Named tuning values. +pub mod config { + pub mod sdf_collisions_config; +} +pub mod benchmark; +pub mod hud; +pub mod pile_arena; +pub mod prop_library; +pub mod sdf_collisions_demo; diff --git a/examples/sdf-collisions/src/pile_arena.rs b/examples/sdf-collisions/src/pile_arena.rs new file mode 100644 index 0000000..c95d143 --- /dev/null +++ b/examples/sdf-collisions/src/pile_arena.rs @@ -0,0 +1,326 @@ +// Copyright (c) 2026 Trent Polack. All Rights Reserved. +// Licensed under the MIT License. + +//! One bin and its growing pile: an independent physics world with a static floor and walls, props dropped on a +//! seeded schedule, and (when a scene is supplied) one mesh per body that follows it. +//! +//! The SDF bin gives each prop its signed distance collider. The box bin gives each prop its bounding box, which is +//! what the engine could offer before SDF colliders; boxes are centered on the bounds, so their meshes are offset to +//! match. Both bins draw the same seeded spawn stream and step the same fixed time, so their piles differ only in +//! collision. Worlds use meters, seconds, and Y-up space. Scenes are borrowed per call; `None` runs headless. + +use joy_engine::core::PrimitiveShape; +use joy_engine::core::js_math; +use joy_engine::physics::{BodyHandle, BodyOptions, PhysicsError, SdfCollisionStatistics}; +use joy_engine::prelude::*; +use joy_engine::runtime::Stopwatch; + +use super::config::sdf_collisions_config::{ARENA, PHYSICS, PRESENTATION, SPAWN}; +use super::prop_library::PropLibrary; + +/// Which collider each prop receives. +#[derive(Clone, Copy, Debug, PartialEq, Eq)] +pub enum ColliderMode { + /// The prop's signed distance field. + Sdf, + /// The prop's local bounding box. + BoundingBox, +} + +/// One dropped prop: its body, its prop index, and its mesh when a scene is attached. +struct Piece { + body: BodyHandle, + prop: usize, + mesh: Option, +} + +/// Copied bin diagnostics; timings are CPU milliseconds per fixed step. +#[derive(Clone, Debug, PartialEq)] +pub struct ArenaSnapshot { + /// Props dropped so far. + pub bodies: usize, + /// Props the pile grows to. + pub target: usize, + /// Bodies that left the bin (fell through the floor or out over the walls). + pub escaped: usize, + /// Mean and worst recent physics time per fixed step. + pub mean_step_milliseconds: f64, + pub worst_step_milliseconds: f64, + /// SDF contact work of the latest fixed step (all zero in the box bin). + pub statistics: SdfCollisionStatistics, +} + +/// A bin, its world, and its pile. +pub struct PileArena { + // Configuration. + mode: ColliderMode, + center_x: f32, + + // Simulation. Pieces are in drop order, which is also the order the hash and meshes visit them. + world: PhysicsWorld, + random: SeededRandom, + pieces: Vec, + target: usize, + spawn_clock: f32, + + // Presentation; empty when headless. + bin_meshes: Vec, + + // CPU timing samples in milliseconds per fixed step. + step_samples: Vec, +} + +impl PileArena { + /// An empty bin centered at `center_x`, with its floor top at y = 0. Adds the bin's meshes to `scene` if given. + pub fn new(mode: ColliderMode, center_x: f32, scene: Option<&mut RenderScene>) -> Result { + let mut world = PhysicsWorld::new(PhysicsWorldOptions { + fixed_step: PHYSICS.fixed_step, + max_sub_steps: PHYSICS.max_sub_steps, + ..PhysicsWorldOptions::default() + })?; + let slabs = bin_slabs(center_x); + for (center, half_extents) in &slabs { + world.create_body(BodyOptions { + position: *center, + half_extents: *half_extents, + friction: PHYSICS.friction, + ..BodyOptions::new(PrimitiveShape::Box) + })?; + } + + let mut bin_meshes = Vec::new(); + if let Some(scene) = scene { + let floor_tint = if mode == ColliderMode::Sdf { PRESENTATION.sdf_tint } else { PRESENTATION.box_tint }; + for (index, (center, half_extents)) in slabs.iter().enumerate() { + let tint = if index == 0 { floor_tint } else { PRESENTATION.wall_tint }; + let transform = Transform::new(*center, Vec3::ZERO, *half_extents * 2.0); + bin_meshes + .push(scene.insert(MeshInstance::new(primitive_box()).with_transform(transform).with_tint(tint))); + } + } + + Ok(Self { + mode, + center_x, + world, + random: SeededRandom::new(SPAWN.seed), + pieces: Vec::new(), + target: SPAWN.initial_count, + spawn_clock: 0.0, + bin_meshes, + step_samples: Vec::new(), + }) + } + + /// Drop props that are due, advance physics by `elapsed_seconds` in fixed steps, and move meshes to their bodies. + pub fn update( + &mut self, + elapsed_seconds: f32, + library: &PropLibrary, + mut scene: Option<&mut RenderScene>, + ) -> Result<(), PhysicsError> { + self.spawn_clock += elapsed_seconds; + while self.spawn_clock >= SPAWN.interval_seconds && self.pieces.len() < self.target { + self.spawn_clock -= SPAWN.interval_seconds; + self.spawn(library, scene.as_deref_mut())?; + } + if self.pieces.len() >= self.target { + self.spawn_clock = 0.0; + } + + let stopwatch = Stopwatch::start(); + let steps = self.world.step(elapsed_seconds)?; + if steps > 0 { + retain(&mut self.step_samples, stopwatch.elapsed_milliseconds() / f64::from(steps)); + } + + if let Some(scene) = scene { + for piece in &self.pieces { + let (Some(pose), Some(mesh)) = (self.world.body_transform(piece.body), piece.mesh) else { + continue; + }; + if let Some(instance) = scene.get_mut(mesh) { + instance.transform.position = pose.position - self.mesh_offset(library, piece.prop, pose.rotation); + instance.transform.rotation = pose.rotation; + } + } + } + Ok(()) + } + + /// Grow the pile by `count` props, up to the configured maximum. + pub fn raise_target(&mut self, count: usize) { + self.set_target(self.target + count); + } + + /// Set how many props the pile grows to, up to the configured maximum; props already dropped stay. + pub fn set_target(&mut self, count: usize) { + self.target = count.min(SPAWN.maximum_count); + } + + /// Remove every prop and its mesh, restart the spawn stream, and restore the initial pile size. + pub fn reset(&mut self, mut scene: Option<&mut RenderScene>) { + for piece in self.pieces.drain(..) { + self.world.remove_body(piece.body); + if let (Some(scene), Some(mesh)) = (scene.as_deref_mut(), piece.mesh) { + scene.remove(mesh); + } + } + self.random = SeededRandom::new(SPAWN.seed); + self.target = SPAWN.initial_count; + self.spawn_clock = 0.0; + self.step_samples.clear(); + } + + /// Copied diagnostics. + pub fn snapshot(&self) -> ArenaSnapshot { + let worst = self.step_samples.iter().copied().fold(0.0, f64::max); + ArenaSnapshot { + bodies: self.pieces.len(), + target: self.target, + escaped: self.escaped(), + mean_step_milliseconds: mean(&self.step_samples), + worst_step_milliseconds: worst, + statistics: self.world.sdf_collision_statistics(), + } + } + + /// FNV-1a hash of every body's position bits in drop order; equal hashes mean identical piles. + /// + /// Positions come straight from the solver. Rotations are left out: their Euler decomposition uses platform math, + /// and a pile whose rotations diverged would not keep identical positions anyway. + pub fn position_hash(&self) -> u64 { + const OFFSET: u64 = 0xcbf2_9ce4_8422_2325; + const PRIME: u64 = 0x0000_0100_0000_01b3; + let mut hash = OFFSET; + for piece in &self.pieces { + let Some(pose) = self.world.body_transform(piece.body) else { + continue; + }; + for value in pose.position.to_array() { + for byte in value.to_bits().to_le_bytes() { + hash = (hash ^ u64::from(byte)).wrapping_mul(PRIME); + } + } + } + hash + } + + /// Release the bin's meshes and every prop mesh from `scene`; the world drops with `self`. + #[allow(dead_code, reason = "part of the test API that test/rust drives")] + pub fn dispose(mut self, scene: &mut RenderScene) { + self.reset(Some(scene)); + for mesh in self.bin_meshes.drain(..) { + scene.remove(mesh); + } + } + + /// Drop the next prop of the seeded stream above the bin. + fn spawn(&mut self, library: &PropLibrary, scene: Option<&mut RenderScene>) -> Result<(), PhysicsError> { + // Every drop draws the same six values in the same order, so both bins stay in lockstep. + let prop_index = ((unit(&mut self.random) * library.len() as f32) as usize).min(library.len() - 1); + let offset_x = (unit(&mut self.random) * 2.0 - 1.0) * SPAWN.spread; + let offset_z = (unit(&mut self.random) * 2.0 - 1.0) * SPAWN.spread; + let rotation = Vec3::new( + unit(&mut self.random) * std::f32::consts::TAU, + unit(&mut self.random) * std::f32::consts::TAU, + unit(&mut self.random) * std::f32::consts::TAU, + ); + let prop = library.prop(prop_index); + let origin = Vec3::new(self.center_x + offset_x, SPAWN.height, offset_z); + + let options = BodyOptions { + mass: prop.mass, + rotation, + friction: PHYSICS.friction, + restitution: PHYSICS.restitution, + linear_damping: PHYSICS.linear_damping, + angular_damping: PHYSICS.angular_damping, + ..BodyOptions::default() + }; + let body = match self.mode { + ColliderMode::Sdf => { + self.world.create_sdf_body(&prop.shape, BodyOptions { position: origin, ..options })? + } + ColliderMode::BoundingBox => { + let position = origin + rotate_portably(prop.box_center, rotation); + let half_extents = prop.box_half_extents; + self.world.create_body(BodyOptions { position, half_extents, ..options })? + } + }; + + let mesh = scene.map(|scene| { + let transform = Transform::new(origin, rotation, Vec3::ONE); + scene.insert(MeshInstance::new(prop.asset.clone()).with_transform(transform).with_shading(Shading::Pbr)) + }); + self.pieces.push(Piece { body, prop: prop_index, mesh }); + Ok(()) + } + + /// World offset from a body's origin to its mesh's origin: none for SDF bodies, the rotated box center for boxes. + fn mesh_offset(&self, library: &PropLibrary, prop: usize, rotation: Vec3) -> Vec3 { + match self.mode { + ColliderMode::Sdf => Vec3::ZERO, + ColliderMode::BoundingBox => rotate_portably(library.prop(prop).box_center, rotation), + } + } + + /// Bodies below the floor or outside the walls by more than the escape margin. + fn escaped(&self) -> usize { + let limit = ARENA.inner_half_extent + ARENA.wall_thickness + ARENA.escape_margin; + self.pieces + .iter() + .filter_map(|piece| self.world.body_transform(piece.body)) + .filter(|pose| { + let outside = (pose.position.x - self.center_x).abs() > limit || pose.position.z.abs() > limit; + pose.position.y < -ARENA.escape_margin || outside + }) + .count() + } +} + +/// The bin's floor slab followed by its four walls, as (center, half extents) in world meters. +fn bin_slabs(center_x: f32) -> [(Vec3, Vec3); 5] { + let inner = ARENA.inner_half_extent; + let thickness = ARENA.wall_thickness; + let outer = inner + thickness; + let wall_y = ARENA.wall_height * 0.5; + let half_height = ARENA.wall_height * 0.5; + [ + (Vec3::new(center_x, -thickness * 0.5, 0.0), Vec3::new(outer, thickness * 0.5, outer)), + (Vec3::new(center_x - inner - (thickness * 0.5), wall_y, 0.0), Vec3::new(thickness * 0.5, half_height, outer)), + (Vec3::new(center_x + inner + (thickness * 0.5), wall_y, 0.0), Vec3::new(thickness * 0.5, half_height, outer)), + (Vec3::new(center_x, wall_y, -inner - (thickness * 0.5)), Vec3::new(inner, half_height, thickness * 0.5)), + (Vec3::new(center_x, wall_y, inner + (thickness * 0.5)), Vec3::new(inner, half_height, thickness * 0.5)), + ] +} + +/// Rotate `offset` by XYZ Euler radians applied X, then Y, then Z, with portable trigonometry, so box bodies start at +/// the same place on every target (the engine's `quat_from_euler` uses platform math, which is fine for drawing). +fn rotate_portably(offset: Vec3, rotation: Vec3) -> Vec3 { + let [x, y, z] = offset.to_array().map(f64::from); + let [(sine_x, cosine_x), (sine_y, cosine_y), (sine_z, cosine_z)] = + rotation.to_array().map(|angle| (js_math::sin(f64::from(angle)), js_math::cos(f64::from(angle)))); + let (y1, z1) = ((y * cosine_x) - (z * sine_x), (y * sine_x) + (z * cosine_x)); + let (x2, z2) = ((x * cosine_y) + (z1 * sine_y), (z1 * cosine_y) - (x * sine_y)); + let (x3, y3) = ((x2 * cosine_z) - (y1 * sine_z), (x2 * sine_z) + (y1 * cosine_z)); + Vec3::new(x3 as f32, y3 as f32, z2 as f32) +} + +/// Next value of the spawn stream in `[0, 1)`. +fn unit(random: &mut SeededRandom) -> f32 { + random.next_unit() as f32 +} + +/// Keep the newest samples of a bounded timing distribution. +fn retain(samples: &mut Vec, value: f64) { + samples.push(value); + if samples.len() > PRESENTATION.timing_samples { + samples.remove(0); + } +} + +/// Mean of a timing distribution, or zero without samples. +fn mean(samples: &[f64]) -> f64 { + if samples.is_empty() { 0.0 } else { samples.iter().sum::() / samples.len() as f64 } +} diff --git a/examples/sdf-collisions/src/prop_library.rs b/examples/sdf-collisions/src/prop_library.rs new file mode 100644 index 0000000..836249e --- /dev/null +++ b/examples/sdf-collisions/src/prop_library.rs @@ -0,0 +1,140 @@ +// Copyright (c) 2026 Trent Polack. All Rights Reserved. +// Licensed under the MIT License. + +//! The props both bins drop: FORM sources evaluated once into a render mesh, a signed distance collider, and the +//! bounding box the comparison bin uses as its collider. +//! +//! Four props bake their FORM triangles into voxel distance grids; the ring uses the analytical torus field, and its +//! FORM mesh only draws it. Every prop is evaluated once at startup and shared by every body that uses it. + +use joy_engine::form::{CompiledScript, GeometryData}; +use joy_engine::physics::{SdfPrimitive, SdfShape}; +use joy_engine::prelude::*; +use joy_engine::rendering::form_mesh_asset; +use joy_engine::runtime::Stopwatch; +use joy_engine::serde_json::{Map, Value}; + +use super::config::sdf_collisions_config::PHYSICS; + +/// How a prop's collider is built. +#[derive(Clone, Copy, Debug, PartialEq)] +pub enum ColliderSource { + /// Bake the FORM triangles into a voxel distance grid. + BakedMesh, + /// An analytical torus (meters); the FORM mesh only draws it. + Torus { major_radius: f32, minor_radius: f32 }, +} + +/// One authored prop. +pub struct PropDefinition { + /// Display name. + pub name: &'static str, + /// FORM source. + pub source: &'static str, + /// Collider construction. + pub collider: ColliderSource, +} + +/// The authored props, in spawn-table order. +pub const PROPS: [PropDefinition; 5] = [ + PropDefinition { name: "Chair", source: include_str!("../assets/chair.form"), collider: ColliderSource::BakedMesh }, + PropDefinition { name: "Jack", source: include_str!("../assets/jack.form"), collider: ColliderSource::BakedMesh }, + PropDefinition { + name: "Staple", + source: include_str!("../assets/staple.form"), + collider: ColliderSource::BakedMesh, + }, + PropDefinition { name: "I-beam", source: include_str!("../assets/beam.form"), collider: ColliderSource::BakedMesh }, + PropDefinition { + name: "Ring", + source: include_str!("../assets/ring.form"), + collider: ColliderSource::Torus { major_radius: 0.36, minor_radius: 0.1 }, + }, +]; + +/// A prop ready to spawn: its render mesh, SDF collider, bounding-box collider, and mass. +pub struct Prop { + /// Display name. + pub name: &'static str, + /// Render mesh in the prop's local frame (the body frame of the SDF bin). + pub asset: MeshAsset, + /// Signed distance collider in the same local frame. + pub shape: SdfShape, + /// Center and half extents of the local bounding box, used as the comparison bin's box collider. + pub box_center: Vec3, + pub box_half_extents: Vec3, + /// Mass in kilograms: the field's enclosed volume times the configured density, used in both bins. + pub mass: f32, + /// Source triangles, CPU milliseconds spent building the collider, and the collider's field memory in bytes. + pub triangles: usize, + pub collider_milliseconds: f64, + pub field_bytes: usize, +} + +/// Every prop, evaluated and baked once; owners borrow it. +pub struct PropLibrary { + props: Vec, +} + +impl PropLibrary { + /// Evaluate and bake every prop. Fails if a bundled FORM source or its geometry is invalid. + pub fn load() -> GameResult { + let props = PROPS.iter().map(build_prop).collect::>>()?; + Ok(Self { props }) + } + + /// The props in spawn-table order. + pub fn props(&self) -> &[Prop] { + &self.props + } + + /// The prop at a spawn-table index. + pub fn prop(&self, index: usize) -> &Prop { + &self.props[index] + } + + /// Number of props. + pub fn len(&self) -> usize { + self.props.len() + } + + /// Field memory of every collider in bytes; bodies share their prop's field. + pub fn field_bytes(&self) -> usize { + self.props.iter().map(|prop| prop.field_bytes).sum() + } +} + +/// Evaluate one prop's FORM source and build its colliders. +fn build_prop(definition: &PropDefinition) -> GameResult { + let geometry = CompiledScript::parse(definition.source)?.evaluate(&Value::Object(Map::new()), None)?; + let asset = form_mesh_asset(&geometry, Shading::Pbr)?; + + let stopwatch = Stopwatch::start(); + let shape = match definition.collider { + ColliderSource::BakedMesh => { + SdfShape::from_triangle_mesh(&positions(&geometry), &geometry.triangles, PHYSICS.bake)? + } + ColliderSource::Torus { major_radius, minor_radius } => { + SdfShape::primitive(SdfPrimitive::Torus { major_radius, minor_radius })? + } + }; + let collider_milliseconds = stopwatch.elapsed_milliseconds(); + + let (minimum, maximum) = shape.local_bounds(); + Ok(Prop { + name: definition.name, + asset, + box_center: (minimum + maximum) * 0.5, + box_half_extents: (maximum - minimum) * 0.5, + mass: shape.volume() * PHYSICS.density, + triangles: geometry.triangles.len() / 3, + collider_milliseconds, + field_bytes: shape.field_bytes(), + shape, + }) +} + +/// FORM positions (script units, which these sources author in meters) as vectors. +fn positions(geometry: &GeometryData) -> Vec { + geometry.positions.chunks_exact(3).map(|xyz| Vec3::new(xyz[0] as f32, xyz[1] as f32, xyz[2] as f32)).collect() +} diff --git a/examples/sdf-collisions/src/sdf_collisions_demo.rs b/examples/sdf-collisions/src/sdf_collisions_demo.rs new file mode 100644 index 0000000..bc60c98 --- /dev/null +++ b/examples/sdf-collisions/src/sdf_collisions_demo.rs @@ -0,0 +1,292 @@ +// Copyright (c) 2026 Trent Polack. All Rights Reserved. +// Licensed under the MIT License. + +//! The SDF collisions demo: two bins side by side fill with the same FORM props, dropped in the same poses on the same +//! schedule. The left bin collides through signed distance fields; the right through the props' bounding boxes. Both +//! simulations step every frame, so the difference in how the piles settle, and what each costs, is visible at once. +//! +//! Drag orbits, the wheel zooms, Space drops more props, P pauses, and R empties both bins and starts over. + +use joy_engine::prelude::*; +use joy_engine::rendering::ShiftedCamera; +use joy_engine::runtime::{ + Stopwatch, + launch_parameter +}; + +use super::benchmark::{ + BenchmarkReport, + run_benchmark +}; +use super::config::sdf_collisions_config::{ + ARENA, + BENCHMARK, + PRESENTATION, + SPAWN, + VIEW +}; +use super::hud::{ + self, + HudCommand, + HudView, + sidebar_width +}; +use super::pile_arena::{ + ArenaSnapshot, + ColliderMode, + PileArena +}; +use super::prop_library::PropLibrary; + +/// Input action held to orbit (primary mouse button or a touch). +pub const ORBIT_ACTION: &str = "orbit"; +/// Input actions for DROP MORE, PAUSE, and RESET. +pub const DROP_ACTION: &str = "drop"; +pub const PAUSE_ACTION: &str = "pause"; +pub const RESET_ACTION: &str = "reset"; + +/// The demo's input: drag with the primary button (or a touch), Space, P, and R. +pub fn input_map() -> InputMap { + InputMap::new() + .action(ORBIT_ACTION, [InputBinding::MouseButton(0)]) + .action(DROP_ACTION, [InputBinding::key(key::SPACE)]) + .action(PAUSE_ACTION, [InputBinding::key(key::KEY_P)]) + .action(RESET_ACTION, [InputBinding::key(key::KEY_R)]) +} + +/// Whether the launch parameters ask for the benchmark (`--bench` natively, `?bench=1` in browsers). +pub fn benchmark_requested() -> bool { + launch_parameter(BENCHMARK.parameter).is_some() +} + +/// Run the configured benchmark and return its one-line report. +pub fn benchmark(library: &PropLibrary) -> GameResult { + Ok(run_benchmark(library, BENCHMARK.steps, BENCHMARK.props)?) +} + +/// Copied diagnostics of both bins. +#[derive(Clone, Debug, PartialEq)] +#[allow(dead_code, reason = "part of the capture and test API that test/rust drives")] +pub struct DemoSnapshot { + /// The SDF bin. + pub sdf: ArenaSnapshot, + /// The bounding-box bin. + pub boxes: ArenaSnapshot, + /// Whether both simulations are paused. + pub paused: bool, + /// Frames rendered. + pub frames: u64, +} + +/// The demo: owns the scene, the prop library, both bins, the camera, and HUD state. +pub struct SdfCollisionsDemo { + // Content and simulation. Bins borrow the library's shapes and meshes per call. + scene: RenderScene, + library: PropLibrary, + sdf_bin: PileArena, + box_bin: PileArena, + paused: bool, + + // Presentation and input. + camera: OrbitCamera, + drag: Option, + reset_history: bool, + + // Diagnostics: CPU render-encoding samples in milliseconds, the backend label, and the benchmark line. + frame_count: u64, + render_samples: Vec, + backend: &'static str, + benchmark: Option, +} + +impl SdfCollisionsDemo { + /// Both empty bins and the evaluated props. + pub fn new() -> GameResult { + let library = PropLibrary::load()?; + let mut scene = RenderScene::new(); + let half_separation = ARENA.separation * 0.5; + let sdf_bin = PileArena::new(ColliderMode::Sdf, -half_separation, Some(&mut scene))?; + let box_bin = PileArena::new(ColliderMode::BoundingBox, half_separation, Some(&mut scene))?; + let camera = OrbitCamera { + target: VIEW.target, + distance: VIEW.distance, + yaw: VIEW.yaw, + pitch: VIEW.pitch, + near_clip: Some(VIEW.clips.0), + far_clip: Some(VIEW.clips.1), + ..OrbitCamera::default() + }; + Ok(Self { + scene, + library, + sdf_bin, + box_bin, + paused: false, + camera, + drag: None, + reset_history: false, + frame_count: 0, + render_samples: Vec::new(), + backend: "starting", + benchmark: None, + }) + } + + /// Advance both bins by `elapsed_seconds` (bounded per frame) unless paused. + pub fn advance(&mut self, elapsed_seconds: f32) -> GameResult { + if self.paused { + return Ok(()); + } + let elapsed = elapsed_seconds.min(PRESENTATION.maximum_frame_seconds); + self.sdf_bin.update(elapsed, &self.library, Some(&mut self.scene))?; + self.box_bin.update(elapsed, &self.library, Some(&mut self.scene))?; + Ok(()) + } + + /// Grow both piles by a burst of props. + pub fn drop_more(&mut self) { + self.sdf_bin.raise_target(SPAWN.burst); + self.box_bin.raise_target(SPAWN.burst); + } + + /// Empty both bins and restart the shared drop sequence. + pub fn reset(&mut self) { + self.sdf_bin.reset(Some(&mut self.scene)); + self.box_bin.reset(Some(&mut self.scene)); + self.reset_history = true; + } + + /// Pause or resume both simulations. + pub fn toggle_pause(&mut self) { + self.paused = !self.paused; + } + + /// Copied diagnostics. + #[allow(dead_code, reason = "part of the capture and test API that test/rust drives")] + pub fn snapshot(&self) -> DemoSnapshot { + DemoSnapshot { + sdf: self.sdf_bin.snapshot(), + boxes: self.box_bin.snapshot(), + paused: self.paused, + frames: self.frame_count, + } + } + + fn apply(&mut self, command: HudCommand) { + match command { + HudCommand::DropMore => self.drop_more(), + HudCommand::TogglePause => self.toggle_pause(), + HudCommand::Reset => self.reset(), + } + } + + /// The orbit camera with its image centered beside the sidebar. + fn view_camera(&self, viewport: Vec2) -> ShiftedCamera<'_, OrbitCamera> { + ShiftedCamera::between_insets(&self.camera, viewport.x, 0.0, sidebar_width(viewport.x)) + } +} + +impl Game for SdfCollisionsDemo { + /// Record the backend, and run the benchmark when the launch parameters request it. It runs here rather than in + /// [`SdfCollisionsDemo::new`] because the shell routes logs to the console only once it starts. + fn start(&mut self, context: &mut StartContext<'_>) -> GameResult { + self.backend = context.capabilities.backend.label(); + if benchmark_requested() { + let report = benchmark(&self.library)?; + joy_engine::log::info!("{report}"); + self.benchmark = Some(report.to_string()); + } + Ok(()) + } + + fn input(&mut self, frame: &mut Frame<'_>) -> GameResult { + if frame.input.pressed(DROP_ACTION) { + self.drop_more(); + } + if frame.input.pressed(PAUSE_ACTION) { + self.toggle_pause(); + } + if frame.input.pressed(RESET_ACTION) { + self.reset(); + } + + let pointer = frame.input.pointer(); + if frame.input.action(ORBIT_ACTION).down && pointer.active { + if let Some(previous) = self.drag { + let delta = pointer.position - previous; + self.camera.yaw -= delta.x * VIEW.orbit_sensitivity; + self.camera.pitch = (self.camera.pitch + (delta.y * VIEW.orbit_sensitivity)) + .clamp(VIEW.pitch_limits.0, VIEW.pitch_limits.1); + } + self.drag = Some(pointer.position); + } else { + self.drag = None; + } + Ok(()) + } + + fn simulation(&mut self, frame: &mut Frame<'_>) -> GameResult { + self.advance(frame.time.elapsed_seconds) + } + + fn render(&mut self, frame: &mut Frame<'_>, target: &mut RenderTarget<'_>) -> GameResult { + if std::mem::take(&mut self.reset_history) { + target.renderer().reset_history(); + } + let lights = [PRESENTATION.key_light, PRESENTATION.fill_light]; + let inputs = FrameInputs { + lights: &lights, + ambient: Some(PRESENTATION.ambient), + clear_color: PRESENTATION.clear_color, + time_seconds: frame.time.time_seconds as f32, + ..FrameInputs::default() + }; + let camera = self.view_camera(frame.viewport.logical_size); + let stopwatch = Stopwatch::start(); + target.draw(&self.scene, &camera, &inputs)?; + self.render_samples.push(stopwatch.elapsed_milliseconds()); + if self.render_samples.len() > PRESENTATION.timing_samples { + self.render_samples.remove(0); + } + self.frame_count += 1; + + let sdf = self.sdf_bin.snapshot(); + let boxes = self.box_bin.snapshot(); + frame.dev_tools.set("SDF bin", format!("{} props · {:.2} ms/step", sdf.bodies, sdf.mean_step_milliseconds)); + frame.dev_tools.set("box bin", format!("{} props · {:.2} ms/step", boxes.bodies, boxes.mean_step_milliseconds)); + frame + .dev_tools + .set("SDF contacts", format!("{} pairs · {} contacts", sdf.statistics.pairs, sdf.statistics.contacts)); + Ok(()) + } + + fn ui(&mut self, context: &egui::Context, _frame: &mut Frame<'_>) { + let sdf = self.sdf_bin.snapshot(); + let boxes = self.box_bin.snapshot(); + let render_milliseconds = if self.render_samples.is_empty() { + 0.0 + } else { + self.render_samples.iter().sum::() / self.render_samples.len() as f64 + }; + let view = HudView { + sdf: &sdf, + boxes: &boxes, + props: self.library.props(), + field_bytes: self.library.field_bytes(), + paused: self.paused, + backend: self.backend, + render_milliseconds, + benchmark: self.benchmark.as_deref(), + }; + for command in hud::draw(context, &view) { + self.apply(command); + } + + // Wheel zoom belongs to the canvas, not to the scrolling sidebar. + let scroll = context.input(|input| input.smooth_scroll_delta.y); + if scroll != 0.0 && !context.is_pointer_over_egui() { + let distance = self.camera.distance * (-scroll * VIEW.zoom_sensitivity).exp(); + self.camera.distance = distance.clamp(VIEW.distance_limits.0, VIEW.distance_limits.1); + } + } +} diff --git a/package-lock.json b/package-lock.json index a788a7f..2e0e842 100644 --- a/package-lock.json +++ b/package-lock.json @@ -72,6 +72,14 @@ "joy-engine": "file:../.." } }, + "examples/sdf-collisions": { + "name": "@joy-engine/example-sdf-collisions", + "version": "0.1.0", + "license": "MIT", + "dependencies": { + "joy-engine": "file:../.." + } + }, "examples/terrain": { "name": "@joy-engine/example-terrain", "version": "0.1.0", @@ -92,6 +100,10 @@ "resolved": "examples/scene-effects", "link": true }, + "node_modules/@joy-engine/example-sdf-collisions": { + "resolved": "examples/sdf-collisions", + "link": true + }, "node_modules/@joy-engine/example-terrain": { "resolved": "examples/terrain", "link": true diff --git a/rustfmt.toml b/rustfmt.toml index 6be60a8..0799122 100644 --- a/rustfmt.toml +++ b/rustfmt.toml @@ -1,6 +1,3 @@ -# Copyright (c) 2026 Trent Polack. All Rights Reserved. -# Licensed under the MIT License. - -# Matches joy-apps so engine, games, and tools share one Rust layout. +# Matches `joy-apps` so engine, games, and tools share one Rust layout. max_width = 120 use_small_heuristics = "Max" diff --git a/src/physics/mod.rs b/src/physics/mod.rs index 61f524a..d223eb6 100644 --- a/src/physics/mod.rs +++ b/src/physics/mod.rs @@ -16,6 +16,9 @@ //! //! [`destruction`] holds solver-independent destructible chunk structures, and [`PhysicsBackend`] is the //! shape-residency contract they synchronize through; `PhysicsWorld` implements it as the CPU backend. +//! +//! [`SdfShape`] colliders describe solids by signed distance (analytical primitives or meshes baked into voxel grids) +//! and collide through optimization-based contact search; see [`PhysicsWorld::create_sdf_body`]. mod backend_contract; mod convert; @@ -23,6 +26,7 @@ pub mod destruction; mod physics_body; mod physics_vehicle; mod physics_world; +mod sdf; pub use backend_contract::{ PhysicsBackend, PhysicsCollisionShape, PhysicsMutationBatch, PhysicsQueryHit, PhysicsResidency, @@ -37,6 +41,7 @@ pub use physics_vehicle::{ VEHICLE_DEFAULTS, VehicleInput, VehicleOptions, VehicleTuning, WHEEL_DEFAULTS, WheelOptions, }; pub use physics_world::{PhysicsWorld, PhysicsWorldOptions}; +pub use sdf::{SDF_BAKE_DEFAULTS, SdfBakeOptions, SdfCollisionStatistics, SdfPrimitive, SdfShape}; use glam::Vec3; @@ -132,4 +137,7 @@ pub enum PhysicsError { /// A batch added a shape whose center was not finite or whose half extents were not finite and positive. #[error("Physics shape {0} needs a finite center and positive half extents.")] InvalidShape(String), + /// An SDF primitive or mesh was rejected; the message names the problem. + #[error("Invalid SDF shape: {0}.")] + InvalidSdfShape(&'static str), } diff --git a/src/physics/physics_body.rs b/src/physics/physics_body.rs index 7eba333..4bdd592 100644 --- a/src/physics/physics_body.rs +++ b/src/physics/physics_body.rs @@ -17,13 +17,15 @@ use glam::Vec3; use rapier3d::dynamics::{CoefficientCombineRule, RigidBodyBuilder}; use rapier3d::geometry::{ColliderBuilder, SharedShape}; use rapier3d::math::{Pose, Vector}; +use rapier3d::pipeline::ActiveHooks; use super::convert::{ deterministic_ln_1p, euler_from_solver_rotation, finite_number, finite_vector, from_solver_vector, solver_rotation_from_euler, to_solver_vector, }; use super::physics_world::BodyRecord; -use super::{BodyHandle, PhysicsError, PhysicsTransform, PhysicsWorld}; +use super::sdf::SdfCollider; +use super::{BodyHandle, PhysicsError, PhysicsTransform, PhysicsWorld, SdfShape}; /// Smallest accepted sphere radius (JavaScript `Number.EPSILON`), compared in `f64`. const MIN_SPHERE_RADIUS: f64 = f64::EPSILON; @@ -91,43 +93,28 @@ impl PhysicsWorld { /// Options are validated in the TypeScript engine's order (mass, friction, restitution, position, rotation, /// damping, shape dimensions, static planes); invalid configuration leaves the world unchanged. pub fn create_body(&mut self, options: BodyOptions) -> Result { - let mass = finite_number(options.mass, "mass", 0.0, f64::INFINITY)?; - let friction = finite_number(options.friction, "friction", 0.0, f64::INFINITY)?; - let restitution = finite_number(options.restitution, "restitution", 0.0, 1.0)?; - let position = finite_vector(options.position, "position")?; - let rotation = finite_vector(options.rotation, "rotation")?; - let linear_damping = finite_number(options.linear_damping, "linearDamping", 0.0, 1.0)?; - let angular_damping = finite_number(options.angular_damping, "angularDamping", 0.0, 1.0)?; + let validated = ValidatedBody::new(&options)?; let shape = create_shape(&options)?; - if options.shape == PrimitiveShape::Plane && mass != 0.0 { + if options.shape == PrimitiveShape::Plane && validated.mass != 0.0 { return Err(PhysicsError::DynamicPlane); } + Ok(self.insert_body(validated, shape, ActiveHooks::empty())) + } - // Mass 0 is static, exactly like cannon-es. - let dynamic = mass > 0.0; - let builder = if dynamic { RigidBodyBuilder::dynamic() } else { RigidBodyBuilder::fixed() }; - let body = builder - .pose(Pose::from_parts(to_solver_vector(position), solver_rotation_from_euler(rotation))) - .linear_damping(damping_coefficient(linear_damping)) - .angular_damping(damping_coefficient(angular_damping)) - .can_sleep(false) - .build(); - - let collider = ColliderBuilder::new(shape) - .friction(friction) - .friction_combine_rule(CoefficientCombineRule::Multiply) - .restitution(restitution) - .restitution_combine_rule(CoefficientCombineRule::Multiply); - let collider = if dynamic { collider.mass(mass) } else { collider.density(0.0) }; + /// Create a world-owned rigid body whose collider is a signed distance shape, and return its handle. + /// + /// The shape replaces the primitive collider: `options.shape`, `half_extents`, and `radius` are ignored, while mass, + /// pose, friction, restitution, and damping are validated and applied exactly as in [`PhysicsWorld::create_body`]. + /// The shape's local frame is the body frame, and its center of mass and inertia come from the field. The body + /// shares `shape`'s data; dropping the last body and handle releases it. Invalid options leave the world unchanged. + pub fn create_sdf_body(&mut self, shape: &SdfShape, options: BodyOptions) -> Result { + let validated = ValidatedBody::new(&options)?; - let raw = self.solver.bodies.insert(body); - let collider = self.solver.colliders.insert_with_parent(collider.build(), raw, &mut self.solver.bodies); - if let Some(body) = self.solver.bodies.get_mut(raw) { - // Mass is queried (and used by vehicles) before the first step, so resolve it immediately. - body.recompute_mass_properties_from_colliders(&self.solver.colliders); - } - self.bodies.insert(raw, BodyRecord { collider, dynamic, vehicle: None }); - Ok(BodyHandle { world: self.world_id, raw }) + // Rapier recycles a pair's previous contacts while its relative pose barely moves, but SDF manifolds are sparse + // and refine themselves step by step (each new sample can displace a worse cached point), so stale ones let + // bodies tip. Pairs with active hooks are never recycled; the world's default hooks keep every solver flag. + let collider = SharedShape::new(SdfCollider { shape: shape.clone() }); + Ok(self.insert_body(validated, collider, ActiveHooks::FILTER_CONTACT_PAIRS)) } /// Remove a body and any vehicle using it as a chassis. Returns whether a live body was removed. @@ -244,6 +231,66 @@ impl PhysicsWorld { } } +impl PhysicsWorld { + /// Insert a validated body with its single collider; static bodies (mass zero) are fixed, exactly like cannon-es. + fn insert_body(&mut self, validated: ValidatedBody, shape: SharedShape, hooks: ActiveHooks) -> BodyHandle { + let dynamic = validated.mass > 0.0; + let builder = if dynamic { RigidBodyBuilder::dynamic() } else { RigidBodyBuilder::fixed() }; + let body = builder + .pose(Pose::from_parts( + to_solver_vector(validated.position), + solver_rotation_from_euler(validated.rotation), + )) + .linear_damping(damping_coefficient(validated.linear_damping)) + .angular_damping(damping_coefficient(validated.angular_damping)) + .can_sleep(false) + .build(); + + let collider = ColliderBuilder::new(shape) + .friction(validated.friction) + .friction_combine_rule(CoefficientCombineRule::Multiply) + .restitution(validated.restitution) + .restitution_combine_rule(CoefficientCombineRule::Multiply) + .active_hooks(hooks); + let collider = if dynamic { collider.mass(validated.mass) } else { collider.density(0.0) }; + + let raw = self.solver.bodies.insert(body); + let collider = self.solver.colliders.insert_with_parent(collider.build(), raw, &mut self.solver.bodies); + if let Some(body) = self.solver.bodies.get_mut(raw) { + // Mass is queried (and used by vehicles) before the first step, so resolve it immediately. + body.recompute_mass_properties_from_colliders(&self.solver.colliders); + } + self.bodies.insert(raw, BodyRecord { collider, dynamic, vehicle: None }); + BodyHandle { world: self.world_id, raw } + } +} + +/// Body options that passed validation, independent of the collider shape. +struct ValidatedBody { + mass: f32, + friction: f32, + restitution: f32, + position: Vec3, + rotation: Vec3, + linear_damping: f32, + angular_damping: f32, +} + +impl ValidatedBody { + /// Validate in the TypeScript engine's order: mass, friction, restitution, position, rotation, then damping. + fn new(options: &BodyOptions) -> Result { + Ok(Self { + mass: finite_number(options.mass, "mass", 0.0, f64::INFINITY)?, + friction: finite_number(options.friction, "friction", 0.0, f64::INFINITY)?, + restitution: finite_number(options.restitution, "restitution", 0.0, 1.0)?, + position: finite_vector(options.position, "position")?, + rotation: finite_vector(options.rotation, "rotation")?, + linear_damping: finite_number(options.linear_damping, "linearDamping", 0.0, 1.0)?, + angular_damping: finite_number(options.angular_damping, "angularDamping", 0.0, 1.0)?, + }) + } +} + /// Validate collider dimensions and build the solver shape. fn create_shape(options: &BodyOptions) -> Result { match options.shape { diff --git a/src/physics/physics_world.rs b/src/physics/physics_world.rs index a639650..5df5686 100644 --- a/src/physics/physics_world.rs +++ b/src/physics/physics_world.rs @@ -4,17 +4,20 @@ //! The physics world owner: fixed-step advancement with bounded catch-up, body and vehicle bookkeeping, and raycasts. use std::fmt; +use std::sync::Arc; use std::sync::atomic::{AtomicU64, Ordering}; use glam::Vec3; use indexmap::IndexMap; use rapier3d::control::DynamicRayCastVehicleController; use rapier3d::dynamics::RigidBodyHandle; -use rapier3d::geometry::{ColliderHandle, Ray}; +use rapier3d::geometry::{ColliderHandle, NarrowPhase, Ray}; +use rapier3d::parry::query::{DefaultQueryDispatcher, QueryDispatcher}; use rapier3d::pipeline::QueryFilter; use super::convert::{finite_number, finite_vector, from_solver_vector, to_solver_vector}; -use super::{BodyHandle, PhysicsError, RaycastHit, WheelOptions}; +use super::sdf::{SdfCollisionCounters, SdfQueryDispatcher}; +use super::{BodyHandle, PhysicsError, RaycastHit, SdfCollisionStatistics, WheelOptions}; // --------------------------------------------------------------------------------------------------------------------- // Constants @@ -122,6 +125,9 @@ pub struct PhysicsWorld { // mutation batch sequence. pub(super) resident_shapes: IndexMap, pub(super) last_mutation_sequence: Option, + + // SDF contact-search diagnostics for the latest fixed step, shared with the solver's narrow-phase dispatcher. + sdf_counters: Arc, } impl fmt::Debug for PhysicsWorld { @@ -160,6 +166,11 @@ impl PhysicsWorld { solver.gravity = to_solver_vector(gravity); solver.integration_parameters.dt = fixed_step; + // SDF pairs are handled first; every other pair falls through to the solver's default dispatcher unchanged. + let sdf_counters = Arc::new(SdfCollisionCounters::default()); + let dispatcher = SdfQueryDispatcher::new(Arc::clone(&sdf_counters)).chain(DefaultQueryDispatcher); + solver.narrow_phase = NarrowPhase::with_query_dispatcher(dispatcher); + Ok(Self { world_id: NEXT_WORLD_ID.fetch_add(1, Ordering::Relaxed), fixed_step, @@ -172,6 +183,7 @@ impl PhysicsWorld { forced_bodies: Vec::new(), resident_shapes: IndexMap::new(), last_mutation_sequence: None, + sdf_counters, }) } @@ -229,6 +241,12 @@ impl PhysicsWorld { from_solver_vector(self.solver.gravity) } + /// Copied SDF contact-search counts from the most recent fixed step; all zero before the first step or without + /// SDF bodies. + pub fn sdf_collision_statistics(&self) -> SdfCollisionStatistics { + self.sdf_counters.snapshot() + } + // ----------------------------------------------------------------------------------------------------------------- // Internal helpers shared by the body and vehicle modules // ----------------------------------------------------------------------------------------------------------------- @@ -312,6 +330,7 @@ impl PhysicsWorld { record.controller.update_vehicle(dt, queries); } + self.sdf_counters.reset(); self.solver.step(); for handle in self.forced_bodies.drain(..) { diff --git a/src/physics/sdf/LICENSE-sdf-dcd-demo b/src/physics/sdf/LICENSE-sdf-dcd-demo new file mode 100644 index 0000000..e82a25f --- /dev/null +++ b/src/physics/sdf/LICENSE-sdf-dcd-demo @@ -0,0 +1,21 @@ +MIT License + +Copyright (c) 2026 Chris Giles + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE. diff --git a/src/physics/sdf/mod.rs b/src/physics/sdf/mod.rs new file mode 100644 index 0000000..c44504a --- /dev/null +++ b/src/physics/sdf/mod.rs @@ -0,0 +1,58 @@ +// Copyright (c) 2026 Trent Polack. All Rights Reserved. +// Licensed under the MIT License. + +//! Signed distance field (SDF) colliders with optimization-based contact generation. +//! +//! An [`SdfShape`] describes a solid by its signed distance function: analytical primitives, or a voxel grid baked from +//! any triangle mesh (FORM output included) without convex decomposition. Bodies built with +//! [`PhysicsWorld::create_sdf_body`](super::PhysicsWorld::create_sdf_body) collide with each other and with the +//! world's boxes, spheres, and planes; everything else about them (solver, materials, damping, raycasts) works like +//! any other body. +//! +//! Contacts follow Giles and Andrews, *Real-Time Collision Handling for Signed Distance Fields via Caching and +//! Importance Sampling* (Pacific Graphics 2026), adapted from their MIT-licensed reference implementation +//! (, Copyright (c) 2026 Chris Giles; see `LICENSE-sdf-dcd-demo`). For each +//! pair with overlapping bounds, a few seed points are optimized toward the deepest overlap with Nelder-Mead: last +//! step's points (a temporal cache that also warm starts the solver) and a random high-curvature seed. The resulting +//! sparse manifold goes to rapier's solver. +//! +//! This is a prototype. Differences from the paper: +//! - Contact normals and depths are read one-sidedly at the optimized point's projections onto both surfaces (the +//! deeper reading wins) instead of from the deeper field at the point itself (Eq. 5), which tilts resting contacts +//! wherever an edge or corner meets a face. +//! - Cached points are stored in both bodies' frames and carried by their poses rather than advected by velocity +//! (Eq. 6). They refine to a tighter tolerance and keep their carried position unless the refinement is deeper, so +//! points on flat faces stop trading places between steps. +//! - Pairs involving SDF colliders opt out of rapier's contact recycling, which would freeze the cache's step-by-step +//! refinement while bodies settle. +//! - Contact sampling uses fixed per-pair seeded streams, so identical worlds produce identical results on every +//! target; the solver's scalar math keeps native and WebAssembly builds in agreement. +//! +//! Limits: collision remains discrete, so fast or thin bodies can tunnel. SDF bodies ignore meshes, compounds, sensors, +//! shape casts, and continuous collision detection, and only the Nelder-Mead optimizer is implemented. Manifold points +//! spread across small curved-on-flat overlaps, which acts like rolling resistance on rolling SDF bodies. + +mod sdf_collider; +mod sdf_contact; +mod sdf_dispatcher; +mod sdf_primitives; +mod sdf_shape; +mod sdf_voxel; + +pub use sdf_shape::{SDF_BAKE_DEFAULTS, SdfBakeOptions, SdfPrimitive, SdfShape}; + +pub(crate) use sdf_collider::SdfCollider; +pub(crate) use sdf_dispatcher::{SdfCollisionCounters, SdfQueryDispatcher}; + +/// SDF contact search work during the most recent fixed step, for diagnostics. Counts never affect simulation. +#[derive(Clone, Copy, Debug, Default, PartialEq, Eq)] +pub struct SdfCollisionStatistics { + /// Collider pairs involving an SDF whose bounds overlapped. + pub pairs: u32, + /// Manifold points within the solver's contact distance. + pub contacts: u32, + /// Nelder-Mead iterations across all pairs. + pub iterations: u32, + /// Distance field evaluations across all pairs (two per objective evaluation). + pub distance_queries: u32, +} diff --git a/src/physics/sdf/sdf_collider.rs b/src/physics/sdf/sdf_collider.rs new file mode 100644 index 0000000..c54092a --- /dev/null +++ b/src/physics/sdf/sdf_collider.rs @@ -0,0 +1,146 @@ +// Copyright (c) 2026 Trent Polack. All Rights Reserved. +// Licensed under the MIT License. + +//! The solver-side collider for an [`SdfShape`]: a custom parry shape with bounds, mass properties, sphere-traced +//! ray casts, and Newton point projection. Contacts between it and other shapes come from +//! [`SdfQueryDispatcher`](super::sdf_dispatcher::SdfQueryDispatcher); the solver itself treats it like any collider. + +use rapier3d::math::{Real, Vector}; +use rapier3d::parry::bounding_volume::{Aabb, BoundingSphere, BoundingVolume}; +use rapier3d::parry::mass_properties::MassProperties; +use rapier3d::parry::query::{PointProjection, PointQuery, Ray, RayCast, RayIntersection}; +use rapier3d::parry::shape::{FeatureId, Shape, ShapeType, TypedShape}; + +use super::sdf_shape::{SdfShape, SdfShapeData}; + +/// Sphere-tracing steps before a ray gives up. +const RAY_MAX_STEPS: usize = 256; +/// A ray hits once it is closer to the surface than this fraction of the shape's size. +const RAY_HIT_FRACTION: Real = 1.0e-4; +/// Fraction of the distance a ray advances per step; baked grids interpolate distances, which can slightly overestimate +/// them near sharp features. +const RAY_STEP_FRACTION: Real = 0.9; +/// Angle (radians) the shape may rotate before contacts can change; conservative for sharp, non-convex fields. +const CCD_ANGULAR_THICKNESS: Real = std::f32::consts::FRAC_PI_4; + +/// A shared [`SdfShape`] attached to a solver collider. +#[derive(Clone)] +pub(crate) struct SdfCollider { + pub(crate) shape: SdfShape, +} + +impl SdfCollider { + pub(crate) fn data(&self) -> &SdfShapeData { + &self.shape.data + } +} + +impl Shape for SdfCollider { + fn compute_local_aabb(&self) -> Aabb { + self.data().bounds + } + + fn compute_local_bounding_sphere(&self) -> BoundingSphere { + self.data().bounds.bounding_sphere() + } + + fn clone_dyn(&self) -> Box { + Box::new(self.clone()) + } + + fn scale_dyn(&self, _scale: Vector, _num_subdivisions: u32) -> Option> { + // Scaling would need a rebaked field and new seeds; build a new SdfShape instead. + None + } + + fn mass_properties(&self, density: Real) -> MassProperties { + let data = self.data(); + MassProperties::with_inertia_matrix(data.center_of_mass, data.volume * density, data.inertia * density) + } + + fn shape_type(&self) -> ShapeType { + ShapeType::Custom + } + + fn as_typed_shape(&self) -> TypedShape<'_> { + TypedShape::Custom(self) + } + + fn ccd_thickness(&self) -> Real { + self.data().bounds.half_extents().min_element() + } + + fn ccd_angular_thickness(&self) -> Real { + CCD_ANGULAR_THICKNESS + } + + fn is_convex(&self) -> bool { + self.data().convex + } +} + +impl RayCast for SdfCollider { + /// Sphere-trace the field from where the ray enters the bounds. A ray starting inside reports a zero time of + /// impact when `solid`, and otherwise the exit point with an inward normal. + fn cast_local_ray_and_get_normal( + &self, + ray: &Ray, + max_time_of_impact: Real, + solid: bool, + ) -> Option { + let data = self.data(); + let speed = ray.dir.length(); + if speed <= 0.0 || !speed.is_finite() { + return None; + } + let hit_distance = RAY_HIT_FRACTION * data.size; + let bounds = data.bounds.loosened(2.0 * hit_distance); + let mut time = bounds.cast_local_ray(ray, max_time_of_impact, true)?; + + // Inside the shape at the start: solid casts stop immediately, hollow casts march to the exit. + let mut sample = data.sample(ray.point_at(time)); + if sample.distance < 0.0 { + if solid { + return Some(RayIntersection::new(time, Vector::ZERO, FeatureId::Unknown)); + } + for _ in 0..RAY_MAX_STEPS { + sample = data.sample(ray.point_at(time)); + if sample.distance >= 0.0 { + return Some(RayIntersection::new(time, -sample.gradient.normalize_or_zero(), FeatureId::Unknown)); + } + time += (-sample.distance).max(hit_distance) / speed; + if time > max_time_of_impact || !bounds.contains_local_point(ray.point_at(time)) { + return None; + } + } + return None; + } + + for _ in 0..RAY_MAX_STEPS { + sample = data.sample(ray.point_at(time)); + if sample.distance < hit_distance { + return Some(RayIntersection::new(time, sample.gradient.normalize_or_zero(), FeatureId::Unknown)); + } + time += (RAY_STEP_FRACTION * sample.distance) / speed; + if time > max_time_of_impact || !bounds.contains_local_point(ray.point_at(time)) { + return None; + } + } + None + } +} + +impl PointQuery for SdfCollider { + fn project_local_point(&self, point: Vector, solid: bool) -> PointProjection { + let data = self.data(); + let inside = data.sample(point).distance < 0.0; + if inside && solid { + return PointProjection::new(true, point); + } + PointProjection::new(inside, data.project(point)) + } + + fn project_local_point_and_get_feature(&self, point: Vector) -> (PointProjection, FeatureId) { + (self.project_local_point(point, false), FeatureId::Unknown) + } +} diff --git a/src/physics/sdf/sdf_contact.rs b/src/physics/sdf/sdf_contact.rs new file mode 100644 index 0000000..53ef48c --- /dev/null +++ b/src/physics/sdf/sdf_contact.rs @@ -0,0 +1,673 @@ +// Copyright (c) 2026 Trent Polack. All Rights Reserved. +// Licensed under the MIT License. + +//! Contact search between two distance fields: Algorithm 1 of Giles and Andrews 2026. +//! +//! Points inside both shapes are where both fields are negative, so the deepest overlap minimizes +//! `max(φ1, φ2)`. Each step minimizes a smooth upper bound of that maximum (Eq. 2) minus a repulsion from the points +//! already found (Eq. 3), starting from a handful of seeds: last step's points carried along with the bodies (the +//! temporal cache, Sec. 3.5.1) and random high-curvature seeds inside the search domain (Sec. 3.5.2). The result is a +//! sparse, well-spread manifold of at most a few points per pair. +//! +//! Everything is expressed in the first shape's local frame; `pose12` places the second shape in that frame. Lengths +//! are meters, and the paper's unit-size tolerances scale with the smaller shape's size. + +use rapier3d::math::{Pose, Real, Vector}; +use rapier3d::parry::bounding_volume::{Aabb, BoundingVolume}; + +use super::sdf_primitives::{FieldSample, box_field, half_space_field, sphere_field}; +use super::sdf_shape::SdfShapeData; +use super::sdf_voxel::unit; +use crate::core::SeededRandom; + +// --------------------------------------------------------------------------------------------------------------------- +// Settings +// --------------------------------------------------------------------------------------------------------------------- + +/// Search parameters; the defaults are the paper's (Sec. 4) plus the reference implementation's stability additions. +struct ContactSearchSettings { + /// Points carried to the next step (`n_cache`). + cache_size: usize, + /// New stochastic samples per pair and step (`n_rand`). + random_samples: usize, + /// Skip random samples for two convex shapes whose full cache is in contact: they overlap in one convex region, + /// which the cached points already track. Roughly halves the cost of resting convex stacks. + skip_convex_samples: bool, + /// Softmax smoothing `ε` of Eq. 2 for unit-size shapes; scales with the squared size. + softmax_smoothing: Real, + /// Repulsion weight `α` of Eq. 3. + repulsion: Real, + /// Optimizer tolerance `τ` for unit-size shapes, also the duplicate distance; scales with size. + tolerance: Real, + /// Improvement of the objective a new sample needs to displace a cached point, for unit-size shapes. Keeps equally + /// deep points on flat faces from trading places every step, which would defeat solver warm starting. + replacement_threshold: Real, + /// Optimizer iteration limit. + max_iterations: usize, + /// Initial simplex size of new samples, as a fraction of the search domain. + simplex_fraction: Real, + /// Smallest initial simplex of cached points, in tolerances; cached points only need to follow the bodies. + cached_simplex_tolerances: Real, + /// Convergence tolerance of cached points, as a fraction of `tolerance`. + cached_tolerance_fraction: Real, + /// Depth gain (in tolerances) a cached point needs before it leaves its carried position. + cached_improvement: Real, + /// Cached points carried less than this fraction of the tolerance are re-measured without optimizing. + settled_fraction: Real, +} + +const CONTACT_SEARCH: ContactSearchSettings = ContactSearchSettings { + cache_size: 4, + random_samples: 1, + skip_convex_samples: true, + softmax_smoothing: 0.1, + repulsion: 0.01, + tolerance: 0.01, + replacement_threshold: 0.01, + max_iterations: 50, + simplex_fraction: 0.1, + cached_simplex_tolerances: 5.0, + cached_tolerance_fraction: 0.25, + cached_improvement: 0.01, + settled_fraction: 0.1, +}; + +/// Smallest squared gradient length used as a Newton denominator. +const MIN_GRADIENT_SQUARED: Real = 1.0e-8; +/// Newton steps when projecting an optimized point onto a surface, and the precision (in tolerances) that ends them +/// early. Exact fields land in one step; baked grids take a few. +const PROJECTION_STEPS: usize = 4; +const PROJECTION_PRECISION: Real = 0.01; +/// Normals shorter than this (degenerate gradients) cannot orient a contact. +const MIN_NORMAL_LENGTH_SQUARED: Real = 0.25; +/// Nelder-Mead also stops when the objective varies less than this fraction of the tolerance across the simplex. +const FLAT_SIMPLEX_FRACTION: Real = 0.05; +/// Only half-spaces this close to an axis clip the search domain. +const AXIS_ALIGNED_COSINE: Real = 0.999; + +// --------------------------------------------------------------------------------------------------------------------- +// Records +// --------------------------------------------------------------------------------------------------------------------- + +/// One side of a pair, as a distance field in its own local frame. +pub(crate) enum PairShape<'a> { + /// A signed distance shape. + Sdf(&'a SdfShapeData), + /// A solver box, measured by its exact field. + Cuboid { half_extents: Vector }, + /// A solver sphere. + Ball { radius: Real }, + /// A solver half-space through the origin with an outward unit normal; unbounded. + HalfSpace { normal: Vector }, +} + +/// A manifold point carried to the next step, stored in both shapes' frames so it moves with both bodies. +#[derive(Clone, Copy, Debug)] +pub(crate) struct CachedContact { + /// The optimized point in the first shape's frame. + pub(crate) position1: Vector, + /// The same point in the second shape's frame. + pub(crate) position2: Vector, + /// Contact normal in the first shape's frame. + pub(crate) normal1: Vector, + /// Whether the point was a contact (within the prediction distance). + pub(crate) active: bool, +} + +/// One manifold point in the first shape's frame. +#[derive(Clone, Copy, Debug)] +pub(crate) struct ContactPoint { + /// Optimized point `x*` (Eq. 4). + pub(crate) position: Vector, + /// Unit normal from the first shape toward the second. + pub(crate) normal: Vector, + /// Signed gap along the normal in meters; negative when the shapes overlap. + pub(crate) separation: Real, + /// Contact anchor on the first shape's surface. + pub(crate) anchor1: Vector, + /// Contact anchor on the second shape's surface; `anchor2 - anchor1 = normal * separation`. + pub(crate) anchor2: Vector, + /// Cache slot this point was warm started from, if any. + pub(crate) cache_slot: Option, + /// Value of Eq. 2 at the point. + softmax: Real, + /// Eq. 3 against the final manifold, used for ordering. + rank: Real, +} + +/// Work done by a search, for diagnostics. +#[derive(Clone, Copy, Debug, Default)] +pub(crate) struct SearchWork { + /// Optimizer iterations. + pub(crate) iterations: u32, + /// Objective evaluations; each queries both fields. + pub(crate) evaluations: u32, +} + +/// One-sided contact measurement: a surface point of one shape read against the other shape's field. +#[derive(Clone, Copy)] +struct SurfaceReading { + normal: Vector, + separation: Real, + anchor1: Vector, + anchor2: Vector, +} + +/// A cached point re-optimized this step: where it came from and how far the bodies carried it. +#[derive(Clone, Copy)] +struct CacheSeed { + slot: usize, + normal: Vector, + advected: Real, +} + +// --------------------------------------------------------------------------------------------------------------------- +// Search +// --------------------------------------------------------------------------------------------------------------------- + +/// Run Algorithm 1 for one pair and return the manifold sorted best first; its leading points form the next cache. +/// +/// `prediction` is the solver's contact distance in meters: shapes closer than it get contacts, and the search domain +/// grows by half of it on each side. Returns no points when the bounds do not overlap. +pub(crate) fn search_contacts( + shape1: &PairShape<'_>, + shape2: &PairShape<'_>, + pose12: &Pose, + prediction: Real, + cache: &[CachedContact], + random: &mut SeededRandom, + work: &mut SearchWork, +) -> Vec { + let Some(domain) = search_domain(shape1, shape2, pose12, prediction) else { + return Vec::new(); + }; + let Some(size) = pair_size(shape1, shape2) else { + return Vec::new(); + }; + let mut search = ContactSearch { + shape1, + shape2, + pose12, + domain, + tolerance: CONTACT_SEARCH.tolerance * size, + smoothing: CONTACT_SEARCH.softmax_smoothing * (size * size), + points: Vec::new(), + work, + }; + search.work.evaluations += 1; + if !search.may_touch(prediction) { + return Vec::new(); + } + + // Alg. 1, lines 4-10: warm start from the cache, carried with the bodies. Each cached point is stored in both + // frames, so the midpoint moves with the average motion of the two bodies (Eq. 6) without needing velocities. + let cached = cache.len().min(CONTACT_SEARCH.cache_size); + for (slot, entry) in cache[..cached].iter().enumerate() { + let carried2 = pose12.transform_point(entry.position2); + let seed = (entry.position1 + carried2) * 0.5; + let advected = (seed - entry.position1).length(); + search.add_point(seed, Some(CacheSeed { slot, normal: entry.normal1, advected })); + } + + // Alg. 1, lines 12-18: stochastic samples find new contact regions and top the cache back up. + let cached_contacts = cache[..cached].iter().filter(|entry| entry.active).count(); + let convex_pair_tracked = CONTACT_SEARCH.skip_convex_samples + && shape1.convex() + && shape2.convex() + && cached_contacts == CONTACT_SEARCH.cache_size; + let random_samples = if convex_pair_tracked { 0 } else { CONTACT_SEARCH.random_samples }; + let samples = random_samples + CONTACT_SEARCH.cache_size.saturating_sub(cached); + for _ in 0..samples { + let seed = search.stochastic_sample(random); + search.add_point(seed, None); + } + + search.finish(size) +} + +/// Working state of one pair's search. +struct ContactSearch<'a> { + shape1: &'a PairShape<'a>, + shape2: &'a PairShape<'a>, + pose12: &'a Pose, + domain: Aabb, + tolerance: Real, + smoothing: Real, + points: Vec, + work: &'a mut SearchWork, +} + +impl ContactSearch<'_> { + /// Minimize the objective from `seed`, measure the contact there (Sec. 3.4), and add it to the manifold. + fn add_point(&mut self, seed: Vector, cached: Option) { + let point = match cached { + // Bodies that barely moved relative to each other leave a cached point where it converged; re-optimizing + // it would only let the repulsion shuffle it. Settled piles skip most of their optimizer work this way. + Some(cached) if cached.advected < (CONTACT_SEARCH.settled_fraction * self.tolerance) => { + seed.clamp(self.domain.mins, self.domain.maxs) + } + _ => self.optimize(seed, cached), + }; + + // The optimizer converges where both distances are about equal, so the paper's rule (Eq. 5, the more deeply + // penetrating field at x* gives the normal) is ambiguous wherever an edge or corner meets a face, and tilts + // resting contacts sideways. Instead, measure one-sidedly from both surfaces near x*, as SDF-mesh methods do: + // project x* onto each shape's surface and read the other field there. The deeper reading wins; near-ties on + // cached points keep the normal closest to the previous one, so contacts do not flip between steps. + let from_first = self.measure_from_first(point); + let from_second = self.measure_from_second(point); + let reading = match (from_first, from_second) { + (Some(first), Some(second)) => { + let tied = (first.separation - second.separation).abs() < self.tolerance; + match cached { + Some(cached) if tied => { + if first.normal.dot(cached.normal) >= second.normal.dot(cached.normal) { + first + } else { + second + } + } + _ => { + if first.separation <= second.separation { + first + } else { + second + } + } + } + } + (Some(reading), None) | (None, Some(reading)) => reading, + (None, None) => return, + }; + + self.points.push(ContactPoint { + position: point, + normal: reading.normal, + separation: reading.separation, + anchor1: reading.anchor1, + anchor2: reading.anchor2, + cache_slot: cached.map(|cached| cached.slot), + softmax: self.depth(point), + rank: 0.0, + }); + } + + /// Minimize the objective from `seed` with Nelder-Mead and return the point to measure. + /// + /// The initial simplex covers a tenth of the domain. A cached point only has to follow the bodies, so its simplex + /// covers how far it was carried (otherwise it would wander across flat faces), and it converges more tightly: a + /// lagging contact on a rolling body tilts its normal and acts as drag. It also keeps its carried position unless + /// the optimizer found a deeper one: on flat faces the depth is level and only the repulsion moves points, whose + /// balance changes with the order points are found in, so without this they trade places from step to step. + fn optimize(&mut self, seed: Vector, cached: Option) -> Vector { + let mut step = + (self.domain.extents() * CONTACT_SEARCH.simplex_fraction).max(Vector::splat(0.5 * self.tolerance)); + let mut tolerance = self.tolerance; + if let Some(cached) = cached { + let carried = cached.advected.max(CONTACT_SEARCH.cached_simplex_tolerances * self.tolerance); + step = step.min(Vector::splat(carried)); + tolerance *= CONTACT_SEARCH.cached_tolerance_fraction; + } + let minimum = nelder_mead(|point| self.objective(point), seed, &self.domain, step, tolerance); + self.work.iterations += minimum.iterations; + self.work.evaluations += minimum.evaluations; + + if cached.is_none() { + return minimum.point; + } + self.work.evaluations += 2; + let carried_depth = self.depth(seed); + if self.depth(minimum.point) > carried_depth - (CONTACT_SEARCH.cached_improvement * self.tolerance) { + return seed.clamp(self.domain.mins, self.domain.maxs); + } + minimum.point + } + + /// Project `point` onto the first shape's surface and read the second field there: the normal is the second + /// shape's inward surface direction and the separation its signed distance. + fn measure_from_first(&self, point: Vector) -> Option { + let surface1 = self.shape1.project(point, self.tolerance); + let sample = self.sample2(surface1); + let normal = -sample.gradient.normalize_or_zero(); + if normal.length_squared() < MIN_NORMAL_LENGTH_SQUARED { + return None; + } + Some(SurfaceReading { + normal, + separation: sample.distance, + anchor1: surface1, + anchor2: surface1 + (normal * sample.distance), + }) + } + + /// Project `point` onto the second shape's surface and read the first field there. + fn measure_from_second(&self, point: Vector) -> Option { + let local = self.pose12.inverse_transform_point(point); + let surface2 = self.pose12.transform_point(self.shape2.project(local, self.tolerance)); + let sample = self.shape1.sample(surface2); + let normal = sample.gradient.normalize_or_zero(); + if normal.length_squared() < MIN_NORMAL_LENGTH_SQUARED { + return None; + } + Some(SurfaceReading { + normal, + separation: sample.distance, + anchor1: surface2 - (normal * sample.distance), + anchor2: surface2, + }) + } + + /// Eq. 2 alone: how deep a point lies in both shapes, without repulsion. + fn depth(&self, point: Vector) -> Real { + softmax(self.shape1.distance(point), self.distance2(point), self.smoothing) + } + + /// The second shape's distance at a point of the first shape's frame. + fn distance2(&self, point: Vector) -> Real { + self.shape2.distance(self.pose12.inverse_transform_point(point)) + } + + /// Whether some point of the domain could lie within `prediction` of both surfaces. Distance fields change at a + /// bounded rate, so a center farther from either surface than that rate times the domain's half-diagonal (plus + /// the prediction distance) proves the domain holds no contact. Overlapping bounds of non-convex shapes often + /// fail this test, which saves the whole search. + fn may_touch(&self, prediction: Real) -> bool { + let center = self.domain.center(); + let radius = self.domain.half_extents().length(); + let reach1 = (self.shape1.lipschitz_bound() * radius) + prediction; + let reach2 = (self.shape2.lipschitz_bound() * radius) + prediction; + self.shape1.distance(center) <= reach1 && self.distance2(center) <= reach2 + } + + /// Eq. 3: the smooth maximum of both distances minus the repulsion from points already in the manifold. + fn objective(&self, point: Vector) -> Real { + let repulsion: Real = self.points.iter().map(|found| (point - found.position).length()).sum(); + self.depth(point) - (CONTACT_SEARCH.repulsion * repulsion) + } + + /// The second shape's field at a point of the first shape's frame, with its gradient in that frame. + fn sample2(&self, point: Vector) -> FieldSample { + let sample = self.shape2.sample(self.pose12.inverse_transform_point(point)); + FieldSample { distance: sample.distance, gradient: self.pose12.transform_vector(sample.gradient) } + } + + /// STOCHASTICSAMPLE (Sec. 3.5.2): a uniformly chosen high-curvature seed of either shape inside the domain, or a + /// uniform point of the domain when no seed lies inside it. + fn stochastic_sample(&self, random: &mut SeededRandom) -> Vector { + let mut sample = + self.domain.mins + (Vector::new(unit(random), unit(random), unit(random)) * self.domain.extents()); + let mut candidates = 0_u32; + let mut consider = |seed: Vector, random: &mut SeededRandom| { + if !self.domain.contains_local_point(seed) { + return; + } + // Reservoir sampling keeps each inside seed with equal probability without collecting them. + candidates += 1; + if ((unit(random) * candidates as Real) as u32).min(candidates - 1) == 0 { + sample = seed; + } + }; + self.shape1.for_each_seed(|seed| consider(seed, random)); + + // The second shape's seeds are culled against the domain's bounds in their own frame first, so only the few + // near the overlap are transformed (the reference implementation culls both shapes this way). + let local_domain = self.domain.transform_by(&self.pose12.inverse()); + self.shape2.for_each_seed(|seed| { + if local_domain.contains_local_point(seed) { + consider(self.pose12.transform_point(seed), random); + } + }); + sample + } + + /// Alg. 1, lines 20-28: drop points that converged within the tolerance of an earlier point (cached points come + /// first and survive), then sort by Eq. 3 against the final manifold, preferring deep points far from the others. + fn finish(self, size: Real) -> Vec { + let duplicate_squared = self.tolerance * self.tolerance; + let mut unique: Vec = Vec::with_capacity(self.points.len()); + for point in &self.points { + if unique.iter().all(|kept| (kept.position - point.position).length_squared() >= duplicate_squared) { + unique.push(*point); + } + } + + let replacement = CONTACT_SEARCH.replacement_threshold * size; + let positions: Vec = unique.iter().map(|point| point.position).collect(); + for point in &mut unique { + let spread: Real = positions.iter().map(|other| (point.position - *other).length()).sum(); + let novelty = if point.cache_slot.is_some() { 0.0 } else { replacement }; + point.rank = point.softmax + novelty - (CONTACT_SEARCH.repulsion * spread); + } + unique.sort_by(|left, right| left.rank.total_cmp(&right.rank)); + unique + } +} + +impl PairShape<'_> { + /// Signed distance and gradient at a point of this shape's frame. + fn sample(&self, point: Vector) -> FieldSample { + match self { + Self::Sdf(data) => data.sample(point), + Self::Cuboid { half_extents } => box_field(point, *half_extents), + Self::Ball { radius } => sphere_field(point, *radius), + Self::HalfSpace { normal } => half_space_field(point, *normal), + } + } + + /// Signed distance alone at a point of this shape's frame. + fn distance(&self, point: Vector) -> Real { + match self { + Self::Sdf(data) => data.distance(point), + _ => self.sample(point).distance, + } + } + + /// Upper bound on the field's rate of change in meters per meter. + fn lipschitz_bound(&self) -> Real { + match self { + Self::Sdf(data) => data.lipschitz_bound(), + Self::Cuboid { .. } | Self::Ball { .. } | Self::HalfSpace { .. } => 1.0, + } + } + + /// Move a point of this shape's frame onto its surface with Newton steps, which also suit inexact fields. + fn project(&self, mut point: Vector, tolerance: Real) -> Vector { + let precision = tolerance * PROJECTION_PRECISION; + for _ in 0..PROJECTION_STEPS { + let sample = self.sample(point); + if sample.distance.abs() < precision { + break; + } + point -= sample.gradient * (sample.distance / gradient_squared(sample)); + } + point + } + + /// Local bounds, or `None` for unbounded half-spaces. + fn local_bounds(&self) -> Option { + match self { + Self::Sdf(data) => Some(data.bounds), + Self::Cuboid { half_extents } => Some(Aabb::new(-*half_extents, *half_extents)), + Self::Ball { radius } => Some(Aabb::new(Vector::splat(-*radius), Vector::splat(*radius))), + Self::HalfSpace { .. } => None, + } + } + + /// Largest side of the bounds, or `None` for half-spaces. + fn size(&self) -> Option { + self.local_bounds().map(|bounds| bounds.extents().max_element()) + } + + fn convex(&self) -> bool { + match self { + Self::Sdf(data) => data.convex, + Self::Cuboid { .. } | Self::Ball { .. } | Self::HalfSpace { .. } => true, + } + } + + /// Visit the high-curvature seeds in this shape's frame. A solver box's corners are its sharpest points. + fn for_each_seed(&self, mut visit: impl FnMut(Vector)) { + match self { + Self::Sdf(data) => data.seeds.iter().for_each(|seed| visit(*seed)), + Self::Cuboid { half_extents } => { + for corner in 0..8 { + let signs = Vector::new( + if corner & 1 == 0 { -1.0 } else { 1.0 }, + if corner & 2 == 0 { -1.0 } else { 1.0 }, + if corner & 4 == 0 { -1.0 } else { 1.0 }, + ); + visit(*half_extents * signs); + } + } + Self::Ball { .. } | Self::HalfSpace { .. } => {} + } + } +} + +// --------------------------------------------------------------------------------------------------------------------- +// Helpers +// --------------------------------------------------------------------------------------------------------------------- + +/// Search domain V (Sec. 3.1): the overlap of both shapes' bounds, grown by half the prediction distance, in the first +/// shape's frame. A half-space contributes no bounds, but clips the domain to its side when it is axis aligned. +fn search_domain(shape1: &PairShape<'_>, shape2: &PairShape<'_>, pose12: &Pose, prediction: Real) -> Option { + let margin = 0.5 * prediction; + let bounds1 = shape1.local_bounds().map(|bounds| bounds.loosened(margin)); + let bounds2 = shape2.local_bounds().map(|bounds| bounds.transform_by(pose12).loosened(margin)); + let mut domain = match (bounds1, bounds2) { + (Some(first), Some(second)) => first.intersection(&second)?, + (Some(bounds), None) | (None, Some(bounds)) => bounds, + (None, None) => return None, + }; + + // A half-space keeps the points with `normal · x <= offset`; the margin lets near contacts in. + let plane = match (shape1, shape2) { + (PairShape::HalfSpace { normal }, _) => Some((*normal, margin)), + (_, PairShape::HalfSpace { normal }) => { + let normal = pose12.transform_vector(*normal); + Some((normal, normal.dot(pose12.translation) + margin)) + } + _ => None, + }; + if let Some((normal, offset)) = plane { + for axis in 0..3 { + if normal[axis] > AXIS_ALIGNED_COSINE { + domain.maxs[axis] = domain.maxs[axis].min(offset / normal[axis]); + } else if normal[axis] < -AXIS_ALIGNED_COSINE { + domain.mins[axis] = domain.mins[axis].max(offset / normal[axis]); + } + } + } + + let valid = domain.mins.is_finite() && domain.maxs.is_finite() && domain.mins.cmple(domain.maxs).all(); + valid.then_some(domain) +} + +/// The smaller bounded shape's size; tolerances scale with it. +fn pair_size(shape1: &PairShape<'_>, shape2: &PairShape<'_>) -> Option { + match (shape1.size(), shape2.size()) { + (Some(first), Some(second)) => Some(first.min(second)), + (Some(size), None) | (None, Some(size)) => Some(size), + (None, None) => None, + } +} + +/// Eq. 2: a smooth upper bound of `max(first, second)`, exact when `smoothing` is zero. +fn softmax(first: Real, second: Real, smoothing: Real) -> Real { + let difference = first - second; + 0.5 * (first + second + ((difference * difference) + smoothing).sqrt()) +} + +fn gradient_squared(sample: FieldSample) -> Real { + sample.gradient.length_squared().max(MIN_GRADIENT_SQUARED) +} + +/// Result of a Nelder-Mead run. +struct Minimum { + point: Vector, + iterations: u32, + evaluations: u32, +} + +/// Derivative-free Nelder-Mead simplex minimization (Sec. 3.3) inside `domain`, with the standard coefficients +/// (reflection 1, expansion 2, contraction and shrink 1/2). Trial points are projected onto the domain. Stops when the +/// simplex fits within `tolerance` of its best vertex, when the objective is flat across it, or after the iteration +/// limit. Only distances are queried, which makes it robust to discontinuous or inexact gradients. +fn nelder_mead( + objective: impl Fn(Vector) -> Real, + start: Vector, + domain: &Aabb, + step: Vector, + tolerance: Real, +) -> Minimum { + let project = |point: Vector| point.clamp(domain.mins, domain.maxs); + let mut evaluations = 0; + let mut evaluate = |point: Vector| { + evaluations += 1; + objective(point) + }; + + // Initial simplex: the start plus one step along each axis, pointing into the domain. + let start = project(start); + let mut vertices = [start; 4]; + for axis in 0..3 { + let toward_inside = if start[axis] + step[axis] <= domain.maxs[axis] { step[axis] } else { -step[axis] }; + vertices[axis + 1][axis] += toward_inside; + vertices[axis + 1] = project(vertices[axis + 1]); + } + let mut values = vertices.map(&mut evaluate); + + let mut iterations: u32 = 0; + while (iterations as usize) < CONTACT_SEARCH.max_iterations { + // Order the vertices from best to worst. + for i in 1..4 { + let mut j = i; + while j > 0 && values[j] < values[j - 1] { + values.swap(j, j - 1); + vertices.swap(j, j - 1); + j -= 1; + } + } + + let extent = (1..4).map(|i| (vertices[i] - vertices[0]).length()).fold(0.0, Real::max); + if extent < tolerance || (values[3] - values[0]) < (FLAT_SIMPLEX_FRACTION * tolerance) { + break; + } + iterations += 1; + + let centroid = (vertices[0] + vertices[1] + vertices[2]) / 3.0; + let reflected = project(centroid + (centroid - vertices[3])); + let reflected_value = evaluate(reflected); + if reflected_value < values[0] { + let expanded = project(centroid + ((centroid - vertices[3]) * 2.0)); + let expanded_value = evaluate(expanded); + if expanded_value < reflected_value { + vertices[3] = expanded; + values[3] = expanded_value; + } else { + vertices[3] = reflected; + values[3] = reflected_value; + } + } else if reflected_value < values[2] { + vertices[3] = reflected; + values[3] = reflected_value; + } else { + // Contract toward the better of the reflected and worst points; shrink toward the best if that fails. + let toward = if reflected_value < values[3] { reflected } else { vertices[3] }; + let contracted = centroid + ((toward - centroid) * 0.5); + let contracted_value = evaluate(contracted); + if contracted_value < reflected_value.min(values[3]) { + vertices[3] = contracted; + values[3] = contracted_value; + } else { + for i in 1..4 { + vertices[i] = vertices[0] + ((vertices[i] - vertices[0]) * 0.5); + values[i] = evaluate(vertices[i]); + } + } + } + } + + let best = (1..4).fold(0, |best, i| if values[i] < values[best] { i } else { best }); + Minimum { point: vertices[best], iterations, evaluations } +} diff --git a/src/physics/sdf/sdf_dispatcher.rs b/src/physics/sdf/sdf_dispatcher.rs new file mode 100644 index 0000000..9f350c5 --- /dev/null +++ b/src/physics/sdf/sdf_dispatcher.rs @@ -0,0 +1,266 @@ +// Copyright (c) 2026 Trent Polack. All Rights Reserved. +// Licensed under the MIT License. + +//! Narrow-phase bridge between the solver and SDF contact search. +//! +//! [`SdfQueryDispatcher`] runs first in the solver's dispatcher chain. Pairs without an SDF collider return +//! `Unsupported`, so the solver's default dispatcher handles them exactly as before. Pairs with one produce one contact +//! manifold per search point: manifolds carry a single normal, and SDF contacts can face many ways at once. The solver +//! merges manifolds with nearly parallel normals before solving. +//! +//! Each pair keeps a workspace with its deterministic random stream and the cached points that warm start the next +//! step. Cached points align with the pair's manifolds, so a carried point also carries its solver impulses. + +use std::sync::Arc; +use std::sync::atomic::{AtomicU32, Ordering}; + +use rapier3d::math::{Pose, Real, Vector}; +use rapier3d::parry::query::details::NormalConstraints; +use rapier3d::parry::query::{ + ClosestPoints, Contact, ContactManifold, ContactManifoldsWorkspace, NonlinearRigidMotion, + PersistentQueryDispatcher, QueryDispatcher, ShapeCastHit, ShapeCastOptions, ShapeDistance, ShapeIntersection, + TrackedContact, TypedWorkspaceData, Unsupported, WorkspaceData, +}; +use rapier3d::parry::shape::{Ball, Cuboid, HalfSpace, PackedFeatureId, Shape}; + +use super::SdfCollisionStatistics; +use super::sdf_collider::SdfCollider; +use super::sdf_contact::{CachedContact, PairShape, SearchWork, search_contacts}; +use crate::core::SeededRandom; + +/// Seed of every pair's sampling stream. Pairs start identically, so identical worlds stay identical. +const PAIR_RANDOM_SEED: i64 = 0x5DF_C0DE; + +/// Collision work of the current fixed step, shared between a world and its dispatcher. +/// +/// These are diagnostics only: counters never feed back into the simulation. +#[derive(Debug, Default)] +pub(crate) struct SdfCollisionCounters { + pairs: AtomicU32, + contacts: AtomicU32, + iterations: AtomicU32, + evaluations: AtomicU32, +} + +impl SdfCollisionCounters { + /// Start counting a new fixed step. + pub(crate) fn reset(&self) { + for counter in [&self.pairs, &self.contacts, &self.iterations, &self.evaluations] { + counter.store(0, Ordering::Relaxed); + } + } + + /// Copy the counts. + pub(crate) fn snapshot(&self) -> SdfCollisionStatistics { + SdfCollisionStatistics { + pairs: self.pairs.load(Ordering::Relaxed), + contacts: self.contacts.load(Ordering::Relaxed), + iterations: self.iterations.load(Ordering::Relaxed), + distance_queries: self.evaluations.load(Ordering::Relaxed).saturating_mul(2), + } + } + + fn record(&self, contacts: u32, work: SearchWork) { + self.pairs.fetch_add(1, Ordering::Relaxed); + self.contacts.fetch_add(contacts, Ordering::Relaxed); + self.iterations.fetch_add(work.iterations, Ordering::Relaxed); + self.evaluations.fetch_add(work.evaluations, Ordering::Relaxed); + } +} + +/// Contact generation for pairs involving an [`SdfCollider`]; chain it before the solver's default dispatcher. +#[derive(Clone)] +pub(crate) struct SdfQueryDispatcher { + counters: Arc, +} + +impl SdfQueryDispatcher { + pub(crate) fn new(counters: Arc) -> Self { + Self { counters } + } +} + +/// Per-pair state kept by the solver between steps. +#[derive(Clone)] +struct SdfPairWorkspace { + random: SeededRandom, + /// Last step's points in manifold order; the search reads its leading entries as the cache. + cache: Vec, +} + +impl WorkspaceData for SdfPairWorkspace { + fn as_typed_workspace_data(&self) -> TypedWorkspaceData<'_> { + TypedWorkspaceData::Custom + } + + fn clone_dyn(&self) -> Box { + Box::new(self.clone()) + } +} + +impl PersistentQueryDispatcher for SdfQueryDispatcher +where + ManifoldData: Default + Clone + Send + Sync, + ContactData: Default + Copy + Send + Sync, +{ + /// Replace the pair's manifolds with this step's SDF contacts; `pos12` places the second shape in the first's + /// frame and `prediction` is the contact distance in meters. + fn contact_manifolds( + &self, + pos12: &Pose, + g1: &dyn Shape, + g2: &dyn Shape, + prediction: Real, + manifolds: &mut Vec>, + workspace: &mut Option, + ) -> Result<(), Unsupported> { + if !is_sdf(g1) && !is_sdf(g2) { + return Err(Unsupported); + } + let (Some(shape1), Some(shape2)) = (pair_shape(g1), pair_shape(g2)) else { + // Shapes without a distance field (meshes, compounds) never touch SDF colliders in this prototype. + manifolds.clear(); + return Ok(()); + }; + + let state = pair_workspace(workspace); + if state.cache.len() != manifolds.len() { + // Something else rewrote the manifolds; their order no longer matches the cache. + state.cache.clear(); + } + let mut work = SearchWork::default(); + let points = search_contacts(&shape1, &shape2, pos12, prediction, &state.cache, &mut state.random, &mut work); + + // One manifold per point. Points carried from the cache reuse their manifold and its solver impulses. + let mut previous: Vec>> = + std::mem::take(manifolds).into_iter().map(Some).collect(); + state.cache.clear(); + let mut contacts = 0; + for point in &points { + let reused = point.cache_slot.and_then(|slot| previous.get_mut(slot).and_then(Option::take)); + let mut manifold = reused.unwrap_or_else(|| ContactManifold::with_data(0, 0, ManifoldData::default())); + let warm_start = manifold.points.first().map(|contact| contact.data).unwrap_or_default(); + + manifold.points.clear(); + manifold.local_n1 = point.normal; + manifold.local_n2 = pos12.inverse_transform_vector(-point.normal); + manifold.subshape1 = 0; + manifold.subshape2 = 0; + manifold.subshape_poses = None; + let mut contact = TrackedContact::new( + point.anchor1, + pos12.inverse_transform_point(point.anchor2), + PackedFeatureId::UNKNOWN, + PackedFeatureId::UNKNOWN, + point.separation, + ); + contact.data = warm_start; + manifold.points.push(contact); + manifolds.push(manifold); + + let active = point.separation < prediction; + contacts += u32::from(active); + state.cache.push(CachedContact { + position1: point.position, + position2: pos12.inverse_transform_point(point.position), + normal1: point.normal, + active, + }); + } + self.counters.record(contacts, work); + Ok(()) + } + + fn contact_manifold_convex_convex( + &self, + _pos12: &Pose, + _g1: &dyn Shape, + _g2: &dyn Shape, + _normal_constraints1: Option<&dyn NormalConstraints>, + _normal_constraints2: Option<&dyn NormalConstraints>, + _prediction: Real, + _manifold: &mut ContactManifold, + ) -> Result<(), Unsupported> { + // Only composite shapes call this for their parts; SDF colliders are never parts. + Err(Unsupported) + } +} + +/// Single-shot queries (sensors, shape casts, continuous collision) are not implemented for SDF colliders; they fall +/// through to the default dispatcher, which reports them unsupported, so SDF bodies are skipped by those queries. +impl QueryDispatcher for SdfQueryDispatcher { + fn intersection_test(&self, _: &Pose, _: &dyn Shape, _: &dyn Shape) -> Result { + Err(Unsupported) + } + + fn distance(&self, _: &Pose, _: &dyn Shape, _: &dyn Shape) -> Result { + Err(Unsupported) + } + + fn contact(&self, _: &Pose, _: &dyn Shape, _: &dyn Shape, _: Real) -> Result, Unsupported> { + Err(Unsupported) + } + + fn closest_points(&self, _: &Pose, _: &dyn Shape, _: &dyn Shape, _: Real) -> Result { + Err(Unsupported) + } + + fn cast_shapes( + &self, + _: &Pose, + _: Vector, + _: &dyn Shape, + _: &dyn Shape, + _: ShapeCastOptions, + ) -> Result, Unsupported> { + Err(Unsupported) + } + + fn cast_shapes_nonlinear( + &self, + _: &NonlinearRigidMotion, + _: &dyn Shape, + _: &NonlinearRigidMotion, + _: &dyn Shape, + _: Real, + _: Real, + _: bool, + ) -> Result, Unsupported> { + Err(Unsupported) + } +} + +fn is_sdf(shape: &dyn Shape) -> bool { + (shape as &dyn std::any::Any).is::() +} + +/// The distance field of a collider, when it has one the search can evaluate. +fn pair_shape(shape: &dyn Shape) -> Option> { + let any = shape as &dyn std::any::Any; + if let Some(sdf) = any.downcast_ref::() { + return Some(PairShape::Sdf(sdf.data())); + } + if let Some(cuboid) = any.downcast_ref::() { + return Some(PairShape::Cuboid { half_extents: cuboid.half_extents }); + } + if let Some(ball) = any.downcast_ref::() { + return Some(PairShape::Ball { radius: ball.radius }); + } + if let Some(half_space) = any.downcast_ref::() { + return Some(PairShape::HalfSpace { normal: half_space.normal }); + } + None +} + +/// This pair's workspace, created on first contact or when the solver dropped it (for example after a shape change). +fn pair_workspace(workspace: &mut Option) -> &mut SdfPairWorkspace { + let ours = workspace.as_ref().is_some_and(|existing| existing.0.is::()); + if !ours { + let fresh = SdfPairWorkspace { random: SeededRandom::new(PAIR_RANDOM_SEED), cache: Vec::new() }; + *workspace = Some(ContactManifoldsWorkspace(Box::new(fresh))); + } + workspace + .as_mut() + .and_then(|existing| existing.0.downcast_mut::()) + .expect("the pair workspace was just ensured to be an SdfPairWorkspace") +} diff --git a/src/physics/sdf/sdf_primitives.rs b/src/physics/sdf/sdf_primitives.rs new file mode 100644 index 0000000..c9b2e95 --- /dev/null +++ b/src/physics/sdf/sdf_primitives.rs @@ -0,0 +1,121 @@ +// Copyright (c) 2026 Trent Polack. All Rights Reserved. +// Licensed under the MIT License. + +//! Analytical signed distance functions with gradients, evaluated in a shape's local frame. +//! +//! Every function returns the signed distance (negative inside) and its gradient, a unit vector pointing away from the +//! nearest surface. Symmetric primitives use +Y as their axis; the torus lies in the XZ plane. Values are meters in the +//! solver's scalar math, so native and WebAssembly builds evaluate them bit for bit alike. + +use rapier3d::math::{Real, Vector}; + +/// Distances below this are treated as zero when normalizing a gradient direction. +const DIRECTION_EPSILON: Real = 1.0e-9; + +/// One evaluation of a field: signed distance in meters and the gradient of that distance. +#[derive(Clone, Copy, Debug, PartialEq)] +pub(crate) struct FieldSample { + /// Signed distance; negative inside the shape. + pub distance: Real, + /// Gradient of the distance. Unit length for exact fields; approximately unit length for baked grids. + pub gradient: Vector, +} + +/// Axis-aligned box centered on the origin. +pub(crate) fn box_field(point: Vector, half_extents: Vector) -> FieldSample { + let offset = point.abs() - half_extents; + let outside = offset.max(Vector::ZERO); + let outside_length = outside.length(); + let inside_distance = offset.max_element().min(0.0); + + // Outside, the gradient points from the closest point on the box; inside, along the nearest face's normal. + let gradient = if outside_length > 0.0 { + (outside / outside_length) * signs(point) + } else { + let axis = largest_axis(offset); + let mut direction = Vector::ZERO; + direction[axis] = sign(point[axis]); + direction + }; + FieldSample { distance: outside_length + inside_distance, gradient } +} + +/// Sphere centered on the origin. +pub(crate) fn sphere_field(point: Vector, radius: Real) -> FieldSample { + let length = point.length(); + let gradient = if length > DIRECTION_EPSILON { point / length } else { Vector::Y }; + FieldSample { distance: length - radius, gradient } +} + +/// Capsule along Y: a segment from `-half_length` to `half_length` swept by `radius`. +pub(crate) fn capsule_field(point: Vector, radius: Real, half_length: Real) -> FieldSample { + let along_axis = point.y.clamp(-half_length, half_length); + sphere_field(Vector::new(point.x, point.y - along_axis, point.z), radius) +} + +/// Solid cylinder along Y with flat caps at `±half_height`. +pub(crate) fn cylinder_field(point: Vector, radius: Real, half_height: Real) -> FieldSample { + let radial = ((point.x * point.x) + (point.z * point.z)).sqrt(); + let radial_offset = radial - radius; + let height_offset = point.y.abs() - half_height; + let radial_outside = radial_offset.max(0.0); + let height_outside = height_offset.max(0.0); + let outside = ((radial_outside * radial_outside) + (height_outside * height_outside)).sqrt(); + + let radial_direction = radial_direction(point, radial); + let gradient = if outside > 0.0 { + (radial_direction * (radial_outside / outside)) + + Vector::new(0.0, sign(point.y) * (height_outside / outside), 0.0) + } else if radial_offset > height_offset { + radial_direction + } else { + Vector::new(0.0, sign(point.y), 0.0) + }; + FieldSample { distance: outside + radial_offset.max(height_offset).min(0.0), gradient } +} + +/// Torus around Y: a circle of `major_radius` in the XZ plane swept by `minor_radius`. +pub(crate) fn torus_field(point: Vector, major_radius: Real, minor_radius: Real) -> FieldSample { + let radial = ((point.x * point.x) + (point.z * point.z)).sqrt(); + let ring_offset = radial - major_radius; + let length = ((ring_offset * ring_offset) + (point.y * point.y)).sqrt(); + + let radial_direction = radial_direction(point, radial); + let gradient = if length > DIRECTION_EPSILON { + (radial_direction * (ring_offset / length)) + Vector::new(0.0, point.y / length, 0.0) + } else { + radial_direction + }; + FieldSample { distance: length - minor_radius, gradient } +} + +/// Half-space whose boundary plane passes through the origin with the unit outward `normal`. +pub(crate) fn half_space_field(point: Vector, normal: Vector) -> FieldSample { + FieldSample { distance: point.dot(normal), gradient: normal } +} + +/// `-1` for negative values and `+1` otherwise, so points on a symmetry plane still receive a direction. +fn sign(value: Real) -> Real { + if value < 0.0 { -1.0 } else { 1.0 } +} + +/// Component-wise [`sign`]. +fn signs(value: Vector) -> Vector { + Vector::new(sign(value.x), sign(value.y), sign(value.z)) +} + +/// Index of the largest component; ties prefer the later axis, matching the reference implementation. +fn largest_axis(value: Vector) -> usize { + if value.x > value.y { + if value.x > value.z { 0 } else { 2 } + } else if value.y > value.z { + 1 + } else { + 2 + } +} + +/// Unit direction away from the Y axis in the XZ plane, or +X on the axis itself. +fn radial_direction(point: Vector, radial: Real) -> Vector { + if radial > DIRECTION_EPSILON { Vector::new(point.x / radial, 0.0, point.z / radial) } else { Vector::X } +} diff --git a/src/physics/sdf/sdf_shape.rs b/src/physics/sdf/sdf_shape.rs new file mode 100644 index 0000000..45266f9 --- /dev/null +++ b/src/physics/sdf/sdf_shape.rs @@ -0,0 +1,461 @@ +// Copyright (c) 2026 Trent Polack. All Rights Reserved. +// Licensed under the MIT License. + +//! Signed distance collision shapes: analytical primitives and voxel grids baked from triangle meshes. +//! +//! An [`SdfShape`] is an immutable, cheaply cloned handle. Construction validates the input, then precomputes +//! everything collision needs from the field alone: tight local bounds, unit-density mass properties, and the +//! high-curvature seed points that steer stochastic contact sampling (Eq. 7 of Giles and Andrews 2026). Bodies created +//! with [`PhysicsWorld::create_sdf_body`](super::super::PhysicsWorld::create_sdf_body) share the handle's data, so one +//! baked grid serves any number of instances. + +use std::fmt; +use std::sync::Arc; + +use glam::Vec3; +use rapier3d::math::{Matrix, Real, Vector}; +use rapier3d::parry::bounding_volume::{Aabb, BoundingVolume}; + +use super::super::PhysicsError; +use super::super::convert::{from_solver_vector, to_solver_vector}; +use super::sdf_primitives::{FieldSample, box_field, capsule_field, cylinder_field, sphere_field, torus_field}; +use super::sdf_voxel::{TriangleMesh, VoxelField, unit}; +use crate::core::SeededRandom; + +// --------------------------------------------------------------------------------------------------------------------- +// Constants +// --------------------------------------------------------------------------------------------------------------------- + +/// Defaults for [`SdfShape::from_triangle_mesh`]: 48 cells along the longest side (about 1.7 MB for a cube-like mesh). +pub const SDF_BAKE_DEFAULTS: SdfBakeOptions = SdfBakeOptions { resolution: 48 }; + +/// Smallest and largest accepted bake resolutions, in cells along the mesh's longest side. +const MIN_BAKE_RESOLUTION: u32 = 8; +const MAX_BAKE_RESOLUTION: u32 = 128; + +/// Integration cells along the longest side when computing mass properties. +const MASS_CELLS: Real = 48.0; + +/// Shapes enclosing less volume than this fraction of their bounding box cannot carry mass. +const MIN_VOLUME_FRACTION: f64 = 1.0e-6; + +/// Surface candidates tested for curvature seeds, and the fixed stream that draws them (seeds are part of the shape, +/// so the same input always yields the same seeds). +const SEED_CANDIDATES: usize = 8000; +const SEED_RANDOM_SEED: i64 = 0x5DF_5EED; + +/// Curvature threshold for a unit-size shape (`|Δφ| > κ`, Eq. 7); it scales with `1 / size`. +const SEED_CURVATURE_THRESHOLD: Real = 4.0; +/// Finite-difference step of the Laplacian, as a fraction of the shape's size. +const SEED_LAPLACIAN_STEP: Real = 0.001; +/// Minimum spacing between seeds, as a fraction of the shape's size, so seeds spread over every sharp region. +const SEED_SPACING: Real = 0.08; +/// Analytical candidates start this far (as a fraction of size) outside the bounds and are projected onto the surface. +const SEED_CANDIDATE_MARGIN: Real = 0.1; +/// Candidates farther than this fraction of size from the surface after projection are discarded. +const SEED_SURFACE_TOLERANCE: Real = 0.001; + +/// Largest rate of change of a baked grid's distance, in meters per meter (see `SdfShapeData::lipschitz_bound`). +const VOXEL_LIPSCHITZ_BOUND: Real = 2.0; + +/// Newton projection limits for moving a point onto the zero level set. +const PROJECTION_ITERATIONS: usize = 32; +const PROJECTION_TOLERANCE: Real = 1.0e-5; +/// Smallest squared gradient length used as a Newton denominator. +const MIN_GRADIENT_SQUARED: Real = 1.0e-8; + +// --------------------------------------------------------------------------------------------------------------------- +// Public records +// --------------------------------------------------------------------------------------------------------------------- + +/// Analytical shapes for [`SdfShape::primitive`]. Dimensions are meters in the shape's local frame; symmetric shapes +/// use +Y as their axis and are centered on the origin. +#[derive(Clone, Copy, Debug, PartialEq)] +pub enum SdfPrimitive { + /// Box with positive half extents. + Box { + /// Half extents along X, Y, and Z. + half_extents: Vec3, + }, + /// Sphere. + Sphere { + /// Radius. + radius: f32, + }, + /// Capsule: a Y segment from `-half_length` to `half_length` swept by `radius`. + Capsule { + /// Radius of the swept sphere. + radius: f32, + /// Half the length of the core segment. + half_length: f32, + }, + /// Solid cylinder along Y with flat caps. + Cylinder { + /// Radius. + radius: f32, + /// Half the height. + half_height: f32, + }, + /// Ring torus around Y in the XZ plane; `minor_radius` must be smaller than `major_radius`. + Torus { + /// Distance from the center to the middle of the tube. + major_radius: f32, + /// Tube radius. + minor_radius: f32, + }, +} + +/// Grid settings for [`SdfShape::from_triangle_mesh`]. +#[derive(Clone, Copy, Debug, PartialEq)] +pub struct SdfBakeOptions { + /// Cells along the mesh's longest side, in `[8, 128]`. Memory grows with the cube: 48 cells take about 1.7 MB for + /// a cube-like mesh, 100 cells (the paper's setting) about 17 MB. + pub resolution: u32, +} + +impl Default for SdfBakeOptions { + fn default() -> Self { + SDF_BAKE_DEFAULTS + } +} + +/// Immutable signed distance collision shape, shared by every body created from it. +/// +/// Cloning copies a reference; the field, bounds, mass properties, and seeds are computed once at construction and +/// released when the last clone and the last body using it are dropped. Coordinates are meters in the shape's local +/// frame, which is the frame of the body that uses it. +#[derive(Clone)] +pub struct SdfShape { + pub(crate) data: Arc, +} + +// --------------------------------------------------------------------------------------------------------------------- +// Shape data +// --------------------------------------------------------------------------------------------------------------------- + +/// The distance field behind a shape. +pub(crate) enum SdfField { + Box { half_extents: Vector }, + Sphere { radius: Real }, + Capsule { radius: Real, half_length: Real }, + Cylinder { radius: Real, half_height: Real }, + Torus { major_radius: Real, minor_radius: Real }, + Voxel(VoxelField), +} + +/// Precomputed collision data of one shape, in its local frame. +pub(crate) struct SdfShapeData { + field: SdfField, + /// Tight local bounds of the zero level set. + pub(crate) bounds: Aabb, + /// Largest side of the bounds; scales tolerances so shapes of any size behave alike. + pub(crate) size: Real, + /// High-curvature surface points for stochastic contact sampling (Sec. 3.5.2). + pub(crate) seeds: Vec, + /// Two convex shapes overlap in one convex region, which lets contact search skip random samples. + pub(crate) convex: bool, + /// Unit-density mass properties: volume (cubic meters), center of mass, and inertia about the center of mass. + pub(crate) volume: Real, + pub(crate) center_of_mass: Vector, + pub(crate) inertia: Matrix, +} + +impl SdfShape { + /// Create an analytical shape. Fails unless every dimension is finite and positive, and a torus's minor radius is + /// smaller than its major radius. + pub fn primitive(primitive: SdfPrimitive) -> Result { + let (field, half_extents, convex) = match primitive { + SdfPrimitive::Box { half_extents } => { + positive(&[half_extents.x, half_extents.y, half_extents.z])?; + let half_extents = to_solver_vector(half_extents); + (SdfField::Box { half_extents }, half_extents, true) + } + SdfPrimitive::Sphere { radius } => { + positive(&[radius])?; + (SdfField::Sphere { radius }, Vector::splat(radius), true) + } + SdfPrimitive::Capsule { radius, half_length } => { + positive(&[radius, half_length])?; + (SdfField::Capsule { radius, half_length }, Vector::new(radius, radius + half_length, radius), true) + } + SdfPrimitive::Cylinder { radius, half_height } => { + positive(&[radius, half_height])?; + (SdfField::Cylinder { radius, half_height }, Vector::new(radius, half_height, radius), true) + } + SdfPrimitive::Torus { major_radius, minor_radius } => { + positive(&[major_radius, minor_radius])?; + if minor_radius >= major_radius { + return Err(PhysicsError::InvalidSdfShape("torus minor radius must be below its major radius")); + } + let outer = major_radius + minor_radius; + (SdfField::Torus { major_radius, minor_radius }, Vector::new(outer, minor_radius, outer), false) + } + }; + let bounds = Aabb::new(-half_extents, half_extents); + let data = SdfShapeData::finalize(field, bounds, convex, SurfaceSource::Projection)?; + Ok(Self { data: Arc::new(data) }) + } + + /// Bake a closed (or nearly closed) triangle mesh into a voxel distance grid. + /// + /// `positions` are meters in the mesh's frame, which becomes the shape's local frame, and `indices` lists three + /// position indices per triangle. Winding may face either way; small holes and overlapping parts are tolerated. + /// Fails on non-finite positions, malformed indices, a resolution outside `[8, 128]`, or a mesh that encloses no + /// volume. Baking is CPU work proportional to the grid size; build shapes once and share them. + pub fn from_triangle_mesh( + positions: &[Vec3], + indices: &[u32], + options: SdfBakeOptions, + ) -> Result { + if !(MIN_BAKE_RESOLUTION..=MAX_BAKE_RESOLUTION).contains(&options.resolution) { + return Err(PhysicsError::OutOfRange { + name: "resolution", + minimum: f64::from(MIN_BAKE_RESOLUTION), + maximum: f64::from(MAX_BAKE_RESOLUTION), + }); + } + let positions: Vec = positions.iter().map(|position| to_solver_vector(*position)).collect(); + let mesh = TriangleMesh::new(&positions, indices).map_err(PhysicsError::InvalidSdfShape)?; + let field = SdfField::Voxel(VoxelField::bake(&mesh, options.resolution)); + let data = SdfShapeData::finalize(field, mesh.bounds(), false, SurfaceSource::Mesh(&mesh))?; + Ok(Self { data: Arc::new(data) }) + } + + /// Signed distance in meters from a local point to the surface; negative inside. + pub fn distance(&self, point: Vec3) -> f32 { + self.data.sample(to_solver_vector(point)).distance + } + + /// Tight local bounds of the surface as `(minimum, maximum)` corners in meters. + pub fn local_bounds(&self) -> (Vec3, Vec3) { + (from_solver_vector(self.data.bounds.mins), from_solver_vector(self.data.bounds.maxs)) + } + + /// Enclosed volume in cubic meters, integrated from the field. + pub fn volume(&self) -> f32 { + self.data.volume + } + + /// Number of precomputed high-curvature seed points. + pub fn seed_count(&self) -> usize { + self.data.seeds.len() + } + + /// Whether the shape is a convex primitive (boxes, spheres, capsules, and cylinders). + pub fn is_convex(&self) -> bool { + self.data.convex + } + + /// Bytes held by the field (the voxel grid for baked meshes; analytical shapes store only their dimensions). + pub fn field_bytes(&self) -> usize { + match &self.data.field { + SdfField::Voxel(voxel) => voxel.memory_bytes(), + _ => std::mem::size_of::(), + } + } +} + +impl fmt::Debug for SdfShape { + fn fmt(&self, formatter: &mut fmt::Formatter<'_>) -> fmt::Result { + formatter + .debug_struct("SdfShape") + .field("bounds", &self.local_bounds()) + .field("volume", &self.data.volume) + .field("seeds", &self.data.seeds.len()) + .field("convex", &self.data.convex) + .finish() + } +} + +/// Where curvature seed candidates come from. +enum SurfaceSource<'a> { + /// Random points around the bounds, projected onto the surface along the gradient. + Projection, + /// Area-weighted points on the source triangles. + Mesh(&'a TriangleMesh), +} + +impl SdfShapeData { + /// Signed distance and gradient at a local point. + pub(crate) fn sample(&self, point: Vector) -> FieldSample { + match &self.field { + SdfField::Box { half_extents } => box_field(point, *half_extents), + SdfField::Sphere { radius } => sphere_field(point, *radius), + SdfField::Capsule { radius, half_length } => capsule_field(point, *radius, *half_length), + SdfField::Cylinder { radius, half_height } => cylinder_field(point, *radius, *half_height), + SdfField::Torus { major_radius, minor_radius } => torus_field(point, *major_radius, *minor_radius), + SdfField::Voxel(voxel) => voxel.sample(point), + } + } + + /// Signed distance alone at a local point; cheaper than [`Self::sample`] for baked grids. + pub(crate) fn distance(&self, point: Vector) -> Real { + match &self.field { + SdfField::Voxel(voxel) => voxel.distance(point), + _ => self.sample(point).distance, + } + } + + /// Upper bound on how fast the distance can change per meter. Exact fields change at most one meter per meter; + /// trilinear interpolation and extrapolation of a baked grid can exceed that, by at most a factor of two. + pub(crate) fn lipschitz_bound(&self) -> Real { + match &self.field { + SdfField::Voxel(_) => VOXEL_LIPSCHITZ_BOUND, + _ => 1.0, + } + } + + /// Move a local point onto the zero level set with Newton steps along the gradient (valid for inexact fields too). + pub(crate) fn project(&self, mut point: Vector) -> Vector { + let tolerance = PROJECTION_TOLERANCE * self.size; + for _ in 0..PROJECTION_ITERATIONS { + let sample = self.sample(point); + if sample.distance.abs() < tolerance { + break; + } + point -= sample.gradient * (sample.distance / sample.gradient.length_squared().max(MIN_GRADIENT_SQUARED)); + } + point + } + + /// Compute mass properties and curvature seeds for a field. + fn finalize(field: SdfField, bounds: Aabb, convex: bool, surface: SurfaceSource<'_>) -> Result { + let size = bounds.extents().max_element(); + let mut data = Self { + field, + bounds, + size, + seeds: Vec::new(), + convex, + volume: 0.0, + center_of_mass: bounds.center(), + inertia: Matrix::IDENTITY, + }; + data.integrate_mass_properties()?; + data.seeds = data.curvature_seeds(surface); + Ok(data) + } + + /// Unit-density volume, center of mass, and inertia by integrating over a voxelization of the field. + /// + /// Moments accumulate in `f64` so large grids do not lose small contributions; each inside cell adds its own + /// second moment (`h² / 12` per axis) so coarse grids keep the right inertia for boxes. + fn integrate_mass_properties(&mut self) -> Result<(), PhysicsError> { + let extents = self.bounds.extents(); + let spacing = self.size / MASS_CELLS; + let counts = [0, 1, 2].map(|axis| ((extents[axis] / spacing).ceil() as usize).max(1)); + let cell = + Vector::new(extents.x / counts[0] as Real, extents.y / counts[1] as Real, extents.z / counts[2] as Real); + let cell_volume = f64::from(cell.x * cell.y * cell.z); + + let mut mass = 0.0_f64; + let mut first = [0.0_f64; 3]; + let mut second = [[0.0_f64; 3]; 3]; + for z in 0..counts[2] { + for y in 0..counts[1] { + for x in 0..counts[0] { + let center = + self.bounds.mins + (cell * Vector::new(x as Real + 0.5, y as Real + 0.5, z as Real + 0.5)); + if self.sample(center).distance >= 0.0 { + continue; + } + let center = [center.x, center.y, center.z].map(f64::from); + mass += cell_volume; + for i in 0..3 { + first[i] += center[i] * cell_volume; + for j in 0..3 { + let own = if i == j { f64::from(cell[i] * cell[i]) / 12.0 } else { 0.0 }; + second[i][j] += ((center[i] * center[j]) + own) * cell_volume; + } + } + } + } + } + + let bounds_volume = f64::from(extents.x) * f64::from(extents.y) * f64::from(extents.z); + if mass <= MIN_VOLUME_FRACTION * bounds_volume { + return Err(PhysicsError::InvalidSdfShape("shape must enclose a volume")); + } + let center = first.map(|moment| moment / mass); + + // Second moment about the center of mass S, then inertia I = tr(S) Id - S. + let mut moment = [[0.0_f64; 3]; 3]; + for i in 0..3 { + for j in 0..3 { + moment[i][j] = second[i][j] - (mass * center[i] * center[j]); + } + } + let trace = moment[0][0] + moment[1][1] + moment[2][2]; + let inertia = |i: usize, j: usize| (if i == j { trace } else { 0.0 } - moment[i][j]) as Real; + + self.volume = mass as Real; + self.center_of_mass = Vector::new(center[0] as Real, center[1] as Real, center[2] as Real); + self.inertia = Matrix::from_cols( + Vector::new(inertia(0, 0), inertia(1, 0), inertia(2, 0)), + Vector::new(inertia(0, 1), inertia(1, 1), inertia(2, 1)), + Vector::new(inertia(0, 2), inertia(1, 2), inertia(2, 2)), + ); + Ok(()) + } + + /// Surface points whose Laplacian magnitude exceeds the curvature threshold, thinned to a minimum spacing. + fn curvature_seeds(&self, surface: SurfaceSource<'_>) -> Vec { + let mut random = SeededRandom::new(SEED_RANDOM_SEED); + let candidates = match surface { + SurfaceSource::Mesh(mesh) => mesh.surface_samples(SEED_CANDIDATES, &mut random), + SurfaceSource::Projection => self.projected_surface_samples(&mut random), + }; + + let threshold = SEED_CURVATURE_THRESHOLD / self.size; + let step = SEED_LAPLACIAN_STEP * self.size; + let spacing_squared = (SEED_SPACING * self.size) * (SEED_SPACING * self.size); + let mut seeds: Vec = Vec::new(); + for candidate in candidates { + if self.laplacian(candidate, step).abs() <= threshold { + continue; + } + // Poisson-disk rejection keeps seeds spread over every high-curvature region instead of piling up. + if seeds.iter().all(|seed| (*seed - candidate).length_squared() >= spacing_squared) { + seeds.push(candidate); + } + } + seeds + } + + /// Random points around the bounds, projected onto the surface; points that fail to converge are dropped. + fn projected_surface_samples(&self, random: &mut SeededRandom) -> Vec { + let region = self.bounds.loosened(SEED_CANDIDATE_MARGIN * self.size); + let extents = region.extents(); + let tolerance = SEED_SURFACE_TOLERANCE * self.size; + let mut samples = Vec::with_capacity(SEED_CANDIDATES); + for _ in 0..SEED_CANDIDATES { + let offset = Vector::new(unit(random), unit(random), unit(random)) * extents; + let point = self.project(region.mins + offset); + if self.sample(point).distance.abs() < tolerance { + samples.push(point); + } + } + samples + } + + /// Laplacian of the field by central differences of its gradient (Eq. 7). + fn laplacian(&self, point: Vector, step: Real) -> Real { + let mut sum = 0.0; + for axis in 0..3 { + let mut offset = Vector::ZERO; + offset[axis] = step; + let ahead = self.sample(point + offset).gradient[axis]; + let behind = self.sample(point - offset).gradient[axis]; + sum += (ahead - behind) / (2.0 * step); + } + sum + } +} + +/// Require finite, positive dimensions. +fn positive(values: &[f32]) -> Result<(), PhysicsError> { + if values.iter().all(|value| value.is_finite() && *value > 0.0) { + return Ok(()); + } + Err(PhysicsError::InvalidSdfShape("dimensions must be finite and positive")) +} diff --git a/src/physics/sdf/sdf_voxel.rs b/src/physics/sdf/sdf_voxel.rs new file mode 100644 index 0000000..0abb187 --- /dev/null +++ b/src/physics/sdf/sdf_voxel.rs @@ -0,0 +1,479 @@ +// Copyright (c) 2026 Trent Polack. All Rights Reserved. +// Licensed under the MIT License. + +//! Voxel signed distance fields baked from triangle meshes (Sec. 3.6 of Giles and Andrews 2026). +//! +//! A regular grid around the mesh stores the signed distance and its gradient at every node; queries interpolate them +//! trilinearly and extrapolate outside the grid by adding the distance to it. Distances come from exact closest points +//! on the triangles. The sign comes from the winding number along grid lines in X, Y, and Z with a majority vote, so +//! procedural meshes with small holes, open seams, or overlapping parts still classify their interior. +//! +//! Positions and distances are meters in the mesh's own frame, which becomes the shape's local frame. + +use rapier3d::math::{Real, Vector}; +use rapier3d::parry::bounding_volume::Aabb; + +use super::sdf_primitives::FieldSample; +use crate::core::{RandomSource, SeededRandom}; + +/// Empty grid cells around the mesh bounds on every side, so near-surface queries stay inside the grid. +const GRID_PADDING_CELLS: usize = 3; + +/// Offsets of the classification lines from the grid nodes, as fractions of a cell. Irrational-looking offsets keep the +/// lines off shared vertices and edges, where a crossing could be counted twice or missed. +const LINE_JITTER: [Real; 2] = [0.123_456_7, 0.345_678_9]; + +/// Triangles whose projected area is below this fraction of the squared mesh size never cross a grid line. +const DEGENERATE_AREA_FRACTION: Real = 1.0e-12; + +/// Three axis directions vote on the sign; a node is inside when at least this many agree. +const INSIDE_VOTES: u8 = 2; + +/// Triangles per bounding volume hierarchy leaf. +const BVH_LEAF_TRIANGLES: usize = 4; + +/// Validated triangle mesh with outward-facing winding, used while baking and for surface sampling. +pub(crate) struct TriangleMesh { + positions: Vec, + triangles: Vec<[usize; 3]>, + bounds: Aabb, +} + +/// A baked distance grid: node `(x, y, z)` sits at `origin + (x, y, z) * cell`. +pub(crate) struct VoxelField { + origin: Vector, + cell: Real, + counts: [usize; 3], + // Per-node signed distances and gradients in X-major order. Distances are stored apart because contact search + // evaluates them far more often than gradients, and a compact array keeps those lookups in cache. + distances: Vec, + gradients: Vec, +} + +/// Where a query falls in the grid: the clamped point, its cell's first node, and its fraction across the cell. +struct GridLookup { + clamped: Vector, + base: usize, + fraction: Vector, +} + +impl TriangleMesh { + /// Validate positions and triangle indices and orient the triangles outward (positive enclosed volume). + /// + /// Fails on non-finite positions, an index count that is not a positive multiple of three, out-of-range indices, + /// or a mesh without spatial extent. Degenerate triangles are kept; they contribute nothing to the bake. + pub(crate) fn new(positions: &[Vector], indices: &[u32]) -> Result { + if positions.iter().any(|position| !position.is_finite()) { + return Err("mesh positions must be finite"); + } + if indices.is_empty() || !indices.len().is_multiple_of(3) { + return Err("mesh indices must describe at least one triangle"); + } + if indices.iter().any(|index| *index as usize >= positions.len()) { + return Err("mesh indices must reference existing positions"); + } + + let mut triangles: Vec<[usize; 3]> = indices + .chunks_exact(3) + .map(|triangle| [triangle[0] as usize, triangle[1] as usize, triangle[2] as usize]) + .collect(); + let bounds = Aabb::from_points(triangles.iter().flatten().map(|index| positions[*index])); + if bounds.extents().max_element() <= 0.0 { + return Err("mesh must have spatial extent"); + } + + // Winding decides which side is inside, so triangles must face outward; flip them all when the enclosed signed + // volume is negative. + let signed_volume: Real = + triangles.iter().map(|[a, b, c]| positions[*a].dot(positions[*b].cross(positions[*c])) / 6.0).sum(); + if signed_volume < 0.0 { + for triangle in &mut triangles { + triangle.swap(1, 2); + } + } + Ok(Self { positions: positions.to_vec(), triangles, bounds }) + } + + /// Tight local bounds of the referenced vertices. + pub(crate) fn bounds(&self) -> Aabb { + self.bounds + } + + /// Area-weighted points on the triangles, drawn from `random`. + pub(crate) fn surface_samples(&self, count: usize, random: &mut SeededRandom) -> Vec { + let mut cumulative_area = Vec::with_capacity(self.triangles.len()); + let mut total_area = 0.0; + for triangle in &self.triangles { + let [a, b, c] = self.vertices(*triangle); + total_area += 0.5 * (b - a).cross(c - a).length(); + cumulative_area.push(total_area); + } + if total_area <= 0.0 { + return Vec::new(); + } + + let mut samples = Vec::with_capacity(count); + for _ in 0..count { + let target = unit(random) * total_area; + let index = cumulative_area.partition_point(|area| *area < target).min(self.triangles.len() - 1); + let [a, b, c] = self.vertices(self.triangles[index]); + + // Uniform barycentric sampling: the square root keeps the density even across the triangle. + let first = unit(random).sqrt(); + let second = unit(random); + samples.push((a * (1.0 - first)) + (b * (first * (1.0 - second))) + (c * (first * second))); + } + samples + } + + fn vertices(&self, triangle: [usize; 3]) -> [Vector; 3] { + triangle.map(|index| self.positions[index]) + } +} + +impl VoxelField { + /// Bake `mesh` into a grid with `resolution` cells along its longest side, plus padding. + pub(crate) fn bake(mesh: &TriangleMesh, resolution: u32) -> Self { + let bounds = mesh.bounds(); + let extents = bounds.extents(); + let cell = extents.max_element() / resolution as Real; + let origin = bounds.mins - (Vector::ONE * (GRID_PADDING_CELLS as Real * cell)); + let counts = [0, 1, 2].map(|axis| (extents[axis] / cell).ceil() as usize + 1 + (2 * GRID_PADDING_CELLS)); + let mut field = Self { origin, cell, counts, distances: Vec::new(), gradients: Vec::new() }; + + let votes = field.classify_inside(mesh); + (field.distances, field.gradients) = field.closest_point_samples(mesh, &votes); + field + } + + /// Signed distance at a local point: trilinear inside the grid, plus the distance to the grid outside it. + pub(crate) fn distance(&self, point: Vector) -> Real { + let lookup = self.lookup(point); + self.interpolate(&lookup, |index| self.distances[index]) + (point - lookup.clamped).length() + } + + /// Signed distance and gradient at a local point. Beyond the grid the gradient leans toward the query point. + pub(crate) fn sample(&self, point: Vector) -> FieldSample { + let lookup = self.lookup(point); + let distance = self.interpolate(&lookup, |index| self.distances[index]); + let mut gradient = self.interpolate(&lookup, |index| self.gradients[index]); + let outside = point - lookup.clamped; + let outside_length = outside.length(); + if outside_length > 0.0 { + gradient = (gradient + (outside / outside_length)).normalize_or_zero(); + } + FieldSample { distance: distance + outside_length, gradient } + } + + /// Bytes held by the grid. + pub(crate) fn memory_bytes(&self) -> usize { + (self.distances.len() * std::mem::size_of::()) + (self.gradients.len() * std::mem::size_of::()) + } + + /// Locate a point's cell, clamping it into the grid first. + fn lookup(&self, point: Vector) -> GridLookup { + let last = self.counts.map(|count| (count - 1) as Real); + let highest = self.origin + (Vector::new(last[0], last[1], last[2]) * self.cell); + let clamped = point.clamp(self.origin, highest); + let grid = (clamped - self.origin) / self.cell; + let cell = [0, 1, 2].map(|axis| (grid[axis] as usize).min(self.counts[axis] - 2)); + let fraction = grid - Vector::new(cell[0] as Real, cell[1] as Real, cell[2] as Real); + GridLookup { clamped, base: self.index(cell), fraction } + } + + /// Trilinear interpolation of per-node values around a lookup, as nested lerps along X, then Y, then Z. + fn interpolate(&self, lookup: &GridLookup, value: impl Fn(usize) -> T) -> T + where + T: Copy + std::ops::Add + std::ops::Sub + std::ops::Mul, + { + let lerp = |from: T, to: T, amount: Real| from + ((to - from) * amount); + let row = self.counts[0]; + let slice = self.counts[0] * self.counts[1]; + let Vector { x, y, z } = lookup.fraction; + let at = |offset: usize| value(lookup.base + offset); + let near = lerp(lerp(at(0), at(1), x), lerp(at(row), at(row + 1), x), y); + let far = lerp(lerp(at(slice), at(slice + 1), x), lerp(at(slice + row), at(slice + row + 1), x), y); + lerp(near, far, z) + } + + /// Count, for every node, how many of the three axis directions place it inside the mesh by winding number. + fn classify_inside(&self, mesh: &TriangleMesh) -> Vec { + let mut votes = vec![0_u8; self.node_count()]; + let mesh_size = mesh.bounds().extents().max_element(); + let degenerate_area = DEGENERATE_AREA_FRACTION * (mesh_size * mesh_size); + for axis in 0..3 { + let across = [(axis + 1) % 3, (axis + 2) % 3]; + let jitter = LINE_JITTER.map(|fraction| fraction * self.cell); + let line_counts = [self.counts[across[0]], self.counts[across[1]]]; + let mut crossings: Vec> = vec![Vec::new(); line_counts[0] * line_counts[1]]; + + // Record where every grid line along `axis` passes through a triangle, and in which direction. + for triangle in &mesh.triangles { + let [v0, v1, v2] = mesh.vertices(*triangle); + let [b, c] = across; + let determinant = ((v1[b] - v0[b]) * (v2[c] - v0[c])) - ((v2[b] - v0[b]) * (v1[c] - v0[c])); + if determinant.abs() < degenerate_area { + continue; + } + let line_range = |dimension: usize, slot: usize| { + let low = v0[dimension].min(v1[dimension]).min(v2[dimension]); + let high = v0[dimension].max(v1[dimension]).max(v2[dimension]); + let first = ((low - self.origin[dimension] - jitter[slot]) / self.cell).ceil().max(0.0) as usize; + let last = ((high - self.origin[dimension] - jitter[slot]) / self.cell).floor(); + (first, if last < 0.0 { None } else { Some((last as usize).min(line_counts[slot] - 1)) }) + }; + let (first_b, Some(last_b)) = line_range(b, 0) else { + continue; + }; + let (first_c, Some(last_c)) = line_range(c, 1) else { + continue; + }; + for k in first_c..=last_c { + for j in first_b..=last_b { + let along_b = self.origin[b] + (j as Real * self.cell) + jitter[0] - v0[b]; + let along_c = self.origin[c] + (k as Real * self.cell) + jitter[1] - v0[c]; + let w1 = ((along_b * (v2[c] - v0[c])) - ((v2[b] - v0[b]) * along_c)) / determinant; + let w2 = (((v1[b] - v0[b]) * along_c) - (along_b * (v1[c] - v0[c]))) / determinant; + if w1 < 0.0 || w2 < 0.0 || (w1 + w2) > 1.0 { + continue; + } + // The determinant's sign is the normal's sign along the line, so negative values enter. + let position = v0[axis] + (w1 * (v1[axis] - v0[axis])) + (w2 * (v2[axis] - v0[axis])); + let direction = if determinant < 0.0 { 1 } else { -1 }; + crossings[j + (k * line_counts[0])].push((position, direction)); + } + } + } + + // Sweep each line in order, accumulating the winding number at every node it passes. + for k in 0..line_counts[1] { + for j in 0..line_counts[0] { + let line = &mut crossings[j + (k * line_counts[0])]; + line.sort_by(|left, right| left.0.total_cmp(&right.0).then(left.1.cmp(&right.1))); + let mut passed = 0; + let mut winding = 0; + for i in 0..self.counts[axis] { + let position = self.origin[axis] + (i as Real * self.cell); + while passed < line.len() && line[passed].0 < position { + winding += line[passed].1; + passed += 1; + } + if winding > 0 { + let mut node = [0; 3]; + node[axis] = i; + node[across[0]] = j; + node[across[1]] = k; + votes[self.index(node)] += 1; + } + } + } + } + } + votes + } + + /// Exact closest-point distances and gradients for every node, signed by the inside votes. + fn closest_point_samples(&self, mesh: &TriangleMesh, votes: &[u8]) -> (Vec, Vec) { + let hierarchy = TriangleHierarchy::new(mesh); + let mut distances = vec![0.0; self.node_count()]; + let mut gradients = vec![Vector::ZERO; self.node_count()]; + let direction_epsilon = self.cell * 1.0e-4; + for z in 0..self.counts[2] { + for y in 0..self.counts[1] { + // The previous node's closest triangle usually stays closest and bounds the next search. + let mut previous_triangle: Option = None; + for x in 0..self.counts[0] { + let point = self.origin + (Vector::new(x as Real, y as Real, z as Real) * self.cell); + let mut best = ClosestPoint { distance_squared: Real::INFINITY, point, triangle: 0 }; + if let Some(triangle) = previous_triangle { + let [a, b, c] = mesh.vertices(mesh.triangles[triangle]); + let closest = closest_point_on_triangle(point, a, b, c); + // Slightly loosened so the same triangle is found again by the full search. + let distance_squared = ((closest - point).length_squared() * 1.0001) + 1.0e-12; + best = ClosestPoint { distance_squared, point: closest, triangle }; + } + hierarchy.closest(mesh, point, &mut best); + previous_triangle = Some(best.triangle); + + let index = self.index([x, y, z]); + let inside = votes[index] >= INSIDE_VOTES; + let distance = (point - best.point).length(); + let gradient = if distance > direction_epsilon { + let away_from_surface = (point - best.point) / distance; + if inside { -away_from_surface } else { away_from_surface } + } else { + let [a, b, c] = mesh.vertices(mesh.triangles[best.triangle]); + (b - a).cross(c - a).normalize_or_zero() + }; + distances[index] = if inside { -distance } else { distance }; + gradients[index] = gradient; + } + } + } + (distances, gradients) + } + + fn node_count(&self) -> usize { + self.counts[0] * self.counts[1] * self.counts[2] + } + + fn index(&self, node: [usize; 3]) -> usize { + node[0] + (self.counts[0] * (node[1] + (self.counts[1] * node[2]))) + } +} + +// --------------------------------------------------------------------------------------------------------------------- +// Closest points +// --------------------------------------------------------------------------------------------------------------------- + +/// Best closest point found so far; searches only accept points nearer than `distance_squared`. +struct ClosestPoint { + distance_squared: Real, + point: Vector, + triangle: usize, +} + +/// One bounding volume hierarchy node. Leaves own `count` triangles starting at `first` in the ordering. +struct HierarchyNode { + bounds: Aabb, + children: Option<[usize; 2]>, + first: usize, + count: usize, +} + +/// Bounding volume hierarchy over a mesh's triangles, split at centroid medians, for closest-point queries. +struct TriangleHierarchy { + nodes: Vec, + order: Vec, +} + +impl TriangleHierarchy { + fn new(mesh: &TriangleMesh) -> Self { + let centroids: Vec = mesh + .triangles + .iter() + .map(|triangle| { + let [a, b, c] = mesh.vertices(*triangle); + (a + b + c) / 3.0 + }) + .collect(); + let mut hierarchy = Self { nodes: Vec::new(), order: (0..mesh.triangles.len()).collect() }; + hierarchy.build(mesh, ¢roids, 0, mesh.triangles.len()); + hierarchy + } + + fn build(&mut self, mesh: &TriangleMesh, centroids: &[Vector], first: usize, count: usize) -> usize { + let range = first..(first + count); + let bounds = Aabb::from_points( + self.order[range.clone()].iter().flat_map(|triangle| mesh.vertices(mesh.triangles[*triangle])), + ); + let index = self.nodes.len(); + self.nodes.push(HierarchyNode { bounds, children: None, first, count }); + if count <= BVH_LEAF_TRIANGLES { + return index; + } + + // Median split along the longest axis of the centroids. + let centers = Aabb::from_points(self.order[range.clone()].iter().map(|triangle| centroids[*triangle])); + let extents = centers.extents(); + let axis = if extents.x > extents.y { + if extents.x > extents.z { 0 } else { 2 } + } else if extents.y > extents.z { + 1 + } else { + 2 + }; + let half = count / 2; + self.order[range].select_nth_unstable_by(half, |left, right| { + centroids[*left][axis].total_cmp(¢roids[*right][axis]).then(left.cmp(right)) + }); + let left = self.build(mesh, centroids, first, half); + let right = self.build(mesh, centroids, first + half, count - half); + self.nodes[index].children = Some([left, right]); + index + } + + /// Improve `best` with any triangle closer to `point`. + fn closest(&self, mesh: &TriangleMesh, point: Vector, best: &mut ClosestPoint) { + let mut stack = vec![0]; + while let Some(node_index) = stack.pop() { + let node = &self.nodes[node_index]; + if distance_squared_to_box(&node.bounds, point) >= best.distance_squared { + continue; + } + let Some([left, right]) = node.children else { + for triangle in &self.order[node.first..(node.first + node.count)] { + let [a, b, c] = mesh.vertices(mesh.triangles[*triangle]); + let closest = closest_point_on_triangle(point, a, b, c); + let distance_squared = (closest - point).length_squared(); + if distance_squared < best.distance_squared { + *best = ClosestPoint { distance_squared, point: closest, triangle: *triangle }; + } + } + continue; + }; + + // Visit the nearer child first so it tightens the bound for the farther one. + let left_first = distance_squared_to_box(&self.nodes[left].bounds, point) + < distance_squared_to_box(&self.nodes[right].bounds, point); + if left_first { + stack.extend([right, left]); + } else { + stack.extend([left, right]); + } + } + } +} + +/// Closest point on triangle `abc` to `point` (Ericson, *Real-Time Collision Detection*, 5.1.5). +fn closest_point_on_triangle(point: Vector, a: Vector, b: Vector, c: Vector) -> Vector { + let ab = b - a; + let ac = c - a; + let ap = point - a; + let d1 = ab.dot(ap); + let d2 = ac.dot(ap); + if d1 <= 0.0 && d2 <= 0.0 { + return a; + } + + let bp = point - b; + let d3 = ab.dot(bp); + let d4 = ac.dot(bp); + if d3 >= 0.0 && d4 <= d3 { + return b; + } + + let vc = (d1 * d4) - (d3 * d2); + if vc <= 0.0 && d1 >= 0.0 && d3 <= 0.0 { + return a + (ab * (d1 / (d1 - d3))); + } + + let cp = point - c; + let d5 = ab.dot(cp); + let d6 = ac.dot(cp); + if d6 >= 0.0 && d5 <= d6 { + return c; + } + + let vb = (d5 * d2) - (d1 * d6); + if vb <= 0.0 && d2 >= 0.0 && d6 <= 0.0 { + return a + (ac * (d2 / (d2 - d6))); + } + + let va = (d3 * d6) - (d5 * d4); + if va <= 0.0 && (d4 - d3) >= 0.0 && (d5 - d6) >= 0.0 { + return b + ((c - b) * ((d4 - d3) / ((d4 - d3) + (d5 - d6)))); + } + + let denominator = 1.0 / (va + vb + vc); + a + (ab * (vb * denominator)) + (ac * (vc * denominator)) +} + +fn distance_squared_to_box(bounds: &Aabb, point: Vector) -> Real { + (point - point.clamp(bounds.mins, bounds.maxs)).length_squared() +} + +/// Next value of a seeded stream as a solver scalar in `[0, 1)`. +pub(crate) fn unit(random: &mut SeededRandom) -> Real { + random.next_unit() as Real +} diff --git a/test/rust/example_sdf_collisions.rs b/test/rust/example_sdf_collisions.rs new file mode 100644 index 0000000..0e59d78 --- /dev/null +++ b/test/rust/example_sdf_collisions.rs @@ -0,0 +1,103 @@ +// Copyright (c) 2026 Trent Polack. All Rights Reserved. +// Licensed under the MIT License. + +//! The SDF COLLISIONS example (`examples/sdf-collisions`): its FORM props bake into SDF colliders that fit their +//! meshes, the headless benchmark fills both bins identically on every run with the SDF pile kept inside its bin, and +//! the demo drops, pauses, adds, and resets props through real key input. +//! +//! The demo test needs a GPU adapter (llvmpipe or SwiftShader is enough); without one it prints a note and returns. +//! Set `JOY_RENDER_CAPTURE_DIR` to keep a PNG capture. Native only. +#![cfg(not(target_arch = "wasm32"))] + +#[path = "../../examples/sdf-collisions/src/mod.rs"] +mod example; + +use example::benchmark::run_benchmark; +use example::config::sdf_collisions_config::SPAWN; +use example::prop_library::{PROPS, PropLibrary}; +use example::sdf_collisions_demo::{SdfCollisionsDemo, input_map}; +use joy_engine::runtime::key; +use joy_engine::{HeadlessOptions, HeadlessShell}; + +#[test] +fn props_bake_into_colliders_that_fit_their_meshes() { + let library = PropLibrary::load().unwrap(); + assert_eq!(library.len(), PROPS.len()); + for prop in library.props() { + // Each collider fits its render mesh's bounds and carries mass and seeds for stochastic sampling. + let mesh_bounds = prop.asset.bounds(); + let (minimum, maximum) = prop.shape.local_bounds(); + assert!( + (minimum - mesh_bounds.min).abs().max_element() < 0.02, + "{}: {minimum} vs {}", + prop.name, + mesh_bounds.min + ); + assert!( + (maximum - mesh_bounds.max).abs().max_element() < 0.02, + "{}: {maximum} vs {}", + prop.name, + mesh_bounds.max + ); + assert!(prop.mass > 1.0 && prop.shape.seed_count() > 0, "{}", prop.name); + assert!(prop.box_half_extents.min_element() > 0.0); + } +} + +#[test] +fn benchmark_piles_repeat_exactly_and_stay_in_their_bins() { + let library = PropLibrary::load().unwrap(); + let first = run_benchmark(&library, 240, 8).unwrap(); + let second = run_benchmark(&library, 240, 8).unwrap(); + assert_eq!(first.sdf.hash, second.sdf.hash, "SDF contact sampling is seeded, so piles repeat bit for bit"); + assert_eq!(first.boxes.hash, second.boxes.hash); + assert_ne!(first.sdf.hash, first.boxes.hash, "the bins differ only in collision, and that shows"); + assert_eq!(first.sdf.snapshot.bodies, 8); + assert_eq!(first.sdf.snapshot.escaped, 0, "no SDF prop falls through the floor or walls"); + assert!(first.sdf.snapshot.statistics.contacts > 0, "{:?}", first.sdf.snapshot.statistics); +} + +#[test] +fn the_demo_drops_pauses_adds_and_resets_props_through_key_input() { + let options = HeadlessOptions { width: 640, height: 400, input: input_map(), ..Default::default() }; + let mut shell = match HeadlessShell::new(options) { + Ok(shell) => shell, + Err(error) => { + eprintln!("skipping GPU test: {error}"); + return; + } + }; + let mut demo = SdfCollisionsDemo::new().unwrap(); + shell.start(&mut demo).unwrap(); + // Both bins receive the same drops on the same schedule (one every 0.3 seconds at 60 frames per second). + shell.run_frames(&mut demo, 60).unwrap(); + let dropped = demo.snapshot(); + assert_eq!(dropped.sdf.bodies, dropped.boxes.bodies); + assert!(dropped.sdf.bodies >= 3, "{dropped:?}"); + if let Ok(directory) = std::env::var("JOY_RENDER_CAPTURE_DIR") { + std::fs::create_dir_all(&directory).unwrap(); + shell.save_png(&std::path::Path::new(&directory).join("example-sdf-collisions.png")).unwrap(); + } + + // P freezes both piles; Space raises both targets; R empties both bins and restores the initial pile size. + press(&mut shell, &mut demo, key::KEY_P); + let paused = demo.snapshot(); + shell.run_frames(&mut demo, 30).unwrap(); + assert!(paused.paused && demo.snapshot().sdf.bodies == paused.sdf.bodies); + press(&mut shell, &mut demo, key::KEY_P); + press(&mut shell, &mut demo, key::SPACE); + assert_eq!(demo.snapshot().sdf.target, SPAWN.initial_count + SPAWN.burst); + assert_eq!(demo.snapshot().boxes.target, SPAWN.initial_count + SPAWN.burst); + press(&mut shell, &mut demo, key::KEY_R); + let reset = demo.snapshot(); + assert_eq!((reset.sdf.bodies, reset.boxes.bodies), (0, 0)); + assert_eq!(reset.sdf.target, SPAWN.initial_count); +} + +/// Press and release a key over two frames. +fn press(shell: &mut HeadlessShell, demo: &mut SdfCollisionsDemo, code: &str) { + shell.input().handle_key(code, true); + shell.run_frames(demo, 1).unwrap(); + shell.input().handle_key(code, false); + shell.run_frames(demo, 1).unwrap(); +} diff --git a/test/rust/main.rs b/test/rust/main.rs index b71a475..bee020c 100644 --- a/test/rust/main.rs +++ b/test/rust/main.rs @@ -14,6 +14,7 @@ mod environment_system; mod example_form_building; mod example_form_destruction; mod example_scene_effects; +mod example_sdf_collisions; mod form_evaluate_cli; mod form_evaluator; mod form_evaluator_parity; diff --git a/test/rust/physics.rs b/test/rust/physics.rs index 8e41020..a3a969b 100644 --- a/test/rust/physics.rs +++ b/test/rust/physics.rs @@ -2,13 +2,13 @@ // Licensed under the MIT License. //! Physics behavior the TypeScript suite protects: resting contact, bounded fixed stepping, ownership, raycast -//! vehicles, transform conventions, and determinism. +//! vehicles, transform conventions, and determinism, plus signed distance field colliders. use glam::Vec3; use joy_engine::core::{PrimitiveShape, quat_from_euler}; use joy_engine::physics::{ - BodyHandle, BodyOptions, PhysicsError, PhysicsWorld, PhysicsWorldOptions, VehicleHandle, VehicleInput, - VehicleOptions, VehicleTuning, WheelOptions, + BodyHandle, BodyOptions, PhysicsError, PhysicsWorld, PhysicsWorldOptions, SdfBakeOptions, SdfPrimitive, SdfShape, + VehicleHandle, VehicleInput, VehicleOptions, VehicleTuning, WheelOptions, }; /// Advance in 60 Hz frames for a number of seconds. @@ -59,6 +59,25 @@ fn make_vehicle(world: &mut PhysicsWorld) -> (BodyHandle, VehicleHandle) { (chassis, vehicle) } +/// An axis-aligned box mesh of `half` extents centered on the origin, optionally without its top face (an open crate). +fn box_mesh(half: Vec3, open_top: bool) -> (Vec, Vec) { + let corner = |index: u32| { + let sign = |bit: u32| if index & bit == 0 { -1.0 } else { 1.0 }; + Vec3::new(sign(1), sign(2), sign(4)) * half + }; + let positions = (0..8).map(corner).collect(); + // Faces as quads (outward winding): -X, +X, -Y, +Y, -Z, +Z. + let faces: [[u32; 4]; 6] = [[0, 4, 6, 2], [1, 3, 7, 5], [0, 1, 5, 4], [2, 6, 7, 3], [0, 2, 3, 1], [4, 5, 7, 6]]; + let mut indices = Vec::new(); + for (face, [a, b, c, d]) in faces.into_iter().enumerate() { + if open_top && face == 3 { + continue; + } + indices.extend([a, b, c, a, c, d]); + } + (positions, indices) +} + #[test] fn dynamic_box_rests_on_static_ground_and_queries_copy_state() { let mut world = PhysicsWorld::new(PhysicsWorldOptions::default()).unwrap(); @@ -313,6 +332,101 @@ fn colliders_and_wheels_follow_the_x_then_y_then_z_rotation_convention() { } } +#[test] +fn sdf_bodies_rest_stack_and_answer_raycasts() { + let mut world = PhysicsWorld::new(PhysicsWorldOptions::default()).unwrap(); + ground(&mut world); + let crate_shape = SdfShape::primitive(SdfPrimitive::Box { half_extents: Vec3::splat(0.5) }).unwrap(); + let ring_shape = SdfShape::primitive(SdfPrimitive::Torus { major_radius: 0.5, minor_radius: 0.15 }).unwrap(); + let (positions, indices) = box_mesh(Vec3::new(0.6, 0.3, 0.4), true); + let baked = SdfShape::from_triangle_mesh(&positions, &indices, SdfBakeOptions { resolution: 16 }).unwrap(); + + // A three-high SDF stack, a non-convex ring, and a baked mesh, all dropped onto a primitive box ground. + let stack: Vec = (0..3) + .map(|level| { + let position = Vec3::new(-3.0, 0.55 + (level as f32 * 1.02), 0.0); + world.create_sdf_body(&crate_shape, BodyOptions { mass: 1.0, position, ..BodyOptions::default() }).unwrap() + }) + .collect(); + let ring = world + .create_sdf_body( + &ring_shape, + BodyOptions { + mass: 1.0, + position: Vec3::new(0.0, 1.0, 0.0), + rotation: Vec3::new(0.3, 0.0, 0.2), + ..BodyOptions::default() + }, + ) + .unwrap(); + let tray = world + .create_sdf_body( + &baked, + BodyOptions { mass: 1.0, position: Vec3::new(3.0, 1.0, 0.0), ..BodyOptions::default() }, + ) + .unwrap(); + simulate(&mut world, 4.0); + + // Everything settles at its resting height without tipping, sinking, or drifting. + for (level, body) in stack.iter().enumerate() { + let transform = world.body_transform(*body).unwrap(); + assert!((transform.position.y - (0.5 + level as f32)).abs() < 0.02, "level {level}: {transform:?}"); + assert!((transform.position.x + 3.0).abs() < 0.02 && transform.rotation.abs().max_element() < 0.02); + assert!(world.body_velocity(*body).unwrap().length() < 0.02); + } + let ring_pose = world.body_transform(ring).unwrap(); + assert!((ring_pose.position.y - 0.15).abs() < 0.01, "{ring_pose:?}"); + assert!(ring_pose.rotation.x.abs() < 0.02 && ring_pose.rotation.z.abs() < 0.02); + assert!((world.body_transform(tray).unwrap().position.y - 0.3).abs() < 0.02); + let statistics = world.sdf_collision_statistics(); + assert!(statistics.pairs >= 5 && statistics.contacts >= 5, "{statistics:?}"); + + // Rays stop at SDF surfaces, pass through the ring's hole to the ground, and report outward normals. + let top = world.raycast(Vec3::new(-3.0, 5.0, 0.0), Vec3::new(-3.0, -1.0, 0.0)).unwrap(); + assert_eq!(top.body, stack[2]); + assert!((top.point.y - 3.0).abs() < 0.03 && top.normal.y > 0.95); + let tube = world.raycast(Vec3::new(0.5, 2.0, 0.0), Vec3::new(0.5, -1.0, 0.0)).unwrap(); + assert_eq!(tube.body, ring); + assert!((tube.point.y - 0.3).abs() < 0.03); + let hole = world.raycast(Vec3::new(0.0, 2.0, 0.0), Vec3::new(0.0, -1.0, 0.0)).unwrap(); + assert_ne!(hole.body, ring); + assert!(hole.point.y.abs() < 0.01); +} + +#[test] +fn sdf_shapes_validate_input_and_bake_open_meshes() { + let invalid = [ + SdfPrimitive::Box { half_extents: Vec3::new(0.5, 0.0, 0.5) }, + SdfPrimitive::Sphere { radius: f32::NAN }, + SdfPrimitive::Torus { major_radius: 0.2, minor_radius: 0.3 }, + ]; + for primitive in invalid { + assert!(matches!(SdfShape::primitive(primitive), Err(PhysicsError::InvalidSdfShape(_))), "{primitive:?}"); + } + let (positions, indices) = box_mesh(Vec3::splat(0.5), false); + let coarse = SdfBakeOptions { resolution: 4 }; + assert!(matches!(SdfShape::from_triangle_mesh(&positions, &indices, coarse), Err(PhysicsError::OutOfRange { .. }))); + let options = SdfBakeOptions { resolution: 16 }; + assert!(SdfShape::from_triangle_mesh(&positions, &indices[..4], options).is_err()); + assert!(SdfShape::from_triangle_mesh(&positions, &[0, 1, 99], options).is_err()); + let flat = [Vec3::ZERO, Vec3::X, Vec3::Z]; + assert!(matches!(SdfShape::from_triangle_mesh(&flat, &[0, 1, 2], options), Err(PhysicsError::InvalidSdfShape(_)))); + + // Procedural meshes are often open or inconsistently wound: an open-top crate with flipped winding still bakes a + // solid whose inside, outside, and volume match the closed box. + let (positions, mut indices) = box_mesh(Vec3::splat(0.5), true); + for triangle in indices.chunks_exact_mut(3) { + triangle.swap(1, 2); + } + let crate_shape = SdfShape::from_triangle_mesh(&positions, &indices, options).unwrap(); + assert!((crate_shape.distance(Vec3::ZERO) + 0.5).abs() < 0.02); + assert!((crate_shape.distance(Vec3::new(1.0, 0.0, 0.0)) - 0.5).abs() < 0.02); + assert!((crate_shape.distance(Vec3::new(0.0, -1.0, 0.0)) - 0.5).abs() < 0.02); + assert!((crate_shape.volume() - 1.0).abs() < 0.05, "{crate_shape:?}"); + assert_eq!(crate_shape.local_bounds(), (Vec3::splat(-0.5), Vec3::splat(0.5))); + assert!(crate_shape.seed_count() > 0 && !crate_shape.is_convex()); +} + #[test] fn identical_worlds_produce_bit_identical_results() { let run = || { @@ -327,12 +441,41 @@ fn identical_worlds_produce_bit_identical_results() { ..BodyOptions::new(PrimitiveShape::Box) }) .unwrap(); + + // SDF contacts sample seeded random streams; tumbling non-convex bodies exercise them every step. + let ring = SdfShape::primitive(SdfPrimitive::Torus { major_radius: 0.4, minor_radius: 0.12 }).unwrap(); + let (positions, indices) = box_mesh(Vec3::new(0.3, 0.2, 0.25), true); + let tray = SdfShape::from_triangle_mesh(&positions, &indices, SdfBakeOptions { resolution: 12 }).unwrap(); + let tumbling = [ + world + .create_sdf_body( + &ring, + BodyOptions { + mass: 1.0, + position: Vec3::new(-0.3, 4.5, 4.0), + rotation: Vec3::new(0.9, 0.2, 0.4), + ..BodyOptions::default() + }, + ) + .unwrap(), + world + .create_sdf_body( + &tray, + BodyOptions { + mass: 1.0, + position: Vec3::new(0.4, 5.5, 4.2), + rotation: Vec3::new(-0.5, 0.7, 0.1), + ..BodyOptions::default() + }, + ) + .unwrap(), + ]; let mut trace = Vec::new(); for frame in 0..240 { let input = VehicleInput { throttle: 1.0, steer: if frame > 90 { 0.6 } else { 0.0 }, brake: 0.0 }; world.set_vehicle_input(vehicle, input).unwrap(); world.step(1.0 / 60.0).unwrap(); - for body in [chassis, crate_box] { + for body in [chassis, crate_box, tumbling[0], tumbling[1]] { let transform = world.body_transform(body).unwrap(); trace.extend(transform.position.to_array().map(f32::to_bits)); trace.extend(transform.rotation.to_array().map(f32::to_bits));