diff --git a/CLAUDE.md b/CLAUDE.md index 9810392..6f5dc26 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -35,10 +35,25 @@ Logging uses `tracing`; set verbosity with `-l debug|info|warn|error|trace` (def ## Workspace layout -Four crates under `crates/`: +Five crates under `crates/`: -- **`crust-core`** — the whole engine as a library (`crust_core`): renderer, integrator, - materials, primitives, lights, BVH, USD import. Everything of substance lives here. +- **`crust-rt`** (lib name `crust_rt`) — the intersection kernel, factored out the way + `openqmc-rs` was, behind a deliberately **Embree-shaped API**: `Geometry` values + (triangle meshes with optional per-vertex shading normals, analytic spheres, round + curve segments, single-level `Instance`s with transform motion blur) attach to a + `SceneBuilder` with per-geometry visibility masks; `commit()` builds the acceleration + structure; `Scene::intersect`/`Scene::occluded` mirror `rtcIntersect1`/`rtcOccluded1`. + Hits are plain `Copy` `RayHit`s carrying `geom_id`/`prim_id` — the kernel never sees + materials. Instanced hits report the *instance's* top-level `geom_id` with the inner + `prim_id`. Internals: watertight Woop-2013 triangles, rounded-cone curves, and the + parallel deterministic SBVH build collapsed to BVH4 (details below). Depends only on + glam + rayon; deliberately swappable for Embree bindings behind the same seam. + +- **`crust-core`** — the engine as a library (`crust_core`): renderer, integrator, + materials, lights, volumes, path guiding, USD import — everything above the + intersection layer, which it consumes from `crust-rt` through `rt_world.rs` + (`WorldBuilder`/`World`: kernel geometries paired with a `geom_id`-indexed + material table). UI-free by design: no progress-bar or image-encoding dependencies; progress is reported through a `ProgressCallback`, and fallible entry points return `crust_core::Error` instead of exiting. @@ -156,14 +171,30 @@ material types, `simple_scene`, `get_settings`). Prefer importing from `crust_co ## Core traits (extension points) -- **`Hittable`** (`hittable.rs`) — `hit(ray, t_min, t_max) -> Option` + `bounding_box()`. - `HitRecord` is `Copy` geometry only (point, normal, `t`, `front_face`); `Hit` pairs it with - a **borrowed** `&dyn Material`, so traversal never touches an `Arc` refcount. Implemented by - `Sphere`, `Triangle`, `SmoothTriangle`, `Mesh`, `HittableList`, `Bvh`. Rendering uses a - **two-level BVH**: `Renderer::new` builds a top-level `Bvh` (`bvh.rs`) over the scene's - `HittableList`, and each imported mesh carries its own nested `Bvh` over triangles - (`usd_import.rs`). `Bvh` is a flat node array with binned-SAH splits and iterative - traversal — deterministic for a given scene. +- **Geometry & intersection** — there is no `Hittable` trait anymore: all intersection + lives in the **`crust-rt`** kernel crate (see the workspace layout), and crust-core + talks to it through `rt_world.rs`: a `WorldBuilder` pairs every attached + `rt::Geometry` with its `Arc` (`attach(...) -> geom_id`), and the + committed `World` resolves kernel hits back to materials **by `geom_id`** — + `World::intersect` returns a `WorldHit { rec: HitRecord, mat, geom_id, prim_id }`, + `World::occluded` is the shadow-ray early-exit query (NEE never uses closest-hit). + `HitRecord` (`hittable.rs`) remains the material-facing `Copy` hit geometry (point, + ray-facing normal, `t`, `front_face`). Inside the kernel: `Sphere`, + triangles with one shared **watertight** intersector (Woop et al. 2013 — dominant-axis + shear, 2D edge functions with f64 fallback on exact-zero ties; no pinholes along + shared edges), `RoundCurves` (sphere-swept cones for hair), and `Instance` (a + committed inner `rt::Scene` placed by a transform — rays transform into local space + with unnormalized direction so `t` carries over, normals map back by + inverse-transpose; an optional end-of-shutter transform lerps per-ray for motion + blur). Per-geometry masks gate intersection on `ray.mask` (`MASK_*` consts, + re-exported from the kernel). The build is reference-based **SBVH** (binned object + SAH + spatial splits gated by the α-overlap test, references clipped via exact + Sutherland-Hodgman for triangles and duplicated across children), runs subtrees in + parallel via `rayon::join` above 4096 refs, is **deterministic** (input-only + decisions, pinned by a build-twice test), and is then **collapsed to BVH4**: 4-wide + SoA nodes whose slab tests run on `Vec4` lanes (`safe_inv` keeps zero-direction + components NaN-free; closest-hit traversal orders lanes near-to-far, occlusion + traversal early-exits). - **`Material`** (`material/material.rs`) — `scatter_importance(r_in, rec) -> Option` used by the integrator (`ScatterSample.delta` marks singular lobes like transmission: never mixed with a @@ -183,10 +214,10 @@ material types, `simple_scene`, `get_settings`). Prefer importing from `crust_co - **`Light`** (`light.rs`) — `sample_point`/`pdf`/`emission`/`material`. The one implementation is **`AreaLight`**: a `LightShape` (pure emitting geometry — `SphereShape`, `RectShape`) paired with the `Arc` its scene geometry carries. - Lights are stored in a `LightList` and their surfaces are also added to `world` as + Lights are stored in a `LightList` and their surfaces are also attached to `world` as emissive geometry (Cornell-box semantics: a light is both light and visible object) — - sharing one `Emissive` Arc, which is how the integrator attributes a bounce-hit - emissive surface to its light (`LightList::find_by_material`, address identity). + the `AreaLight` records the geometry's `geom_id`, which is how the integrator + attributes a bounce-hit emissive surface to its light (`LightList::find_by_geom`). **NEE samples one light per vertex** (uniform pick), so the light strategy's MIS density is `light.pdf / n_lights` — the bounce side evaluates the exact same expression for the light it hit; keep the two sides identical or emission is @@ -212,7 +243,19 @@ The only scene format. `load_scene` opens the stage, imports `RenderSettings` fi camera needs the aspect ratio), then traverses prims with an explicit stack that bakes the Xform hierarchy into world matrices. Schema mapping: -- `UsdGeomMesh` → world-baked triangles wrapped in a `BVHNode`; `UsdGeomSphere` → analytic `Sphere`. +- `UsdGeomMesh` → a *local-space* committed `rt::Scene` (one `TriangleMesh` geometry) + placed by an `rt::Geometry::Instance`; prims with identical points/topology/material + share one kernel scene (content hash + memoized material Arcs, so binding paths + compare by pointer). Non-invertible transforms fall back to world-space baking. + `UsdGeomSphere` → analytic `Sphere` geometry. +- `UsdGeomBasisCurves` → an instanced `rt::Geometry::RoundCurves` batch: `linear` curves + directly, `cubic` (bezier | bspline | catmullRom) flattened at 8 samples per span; widths + (USD diameters) resolve per-vertex / per-curve / constant by array length. +- Any geometry prim may author `crust:rayMask` (int; bit 0 camera, bit 1 shadow, bit 2 + indirect — default all) to hide from ray categories, and `crust:motion:translate` + (float3, world-space) to streak through that translation over the shutter (transform + motion blur; primary rays draw a `K_TIME` shutter sample and every secondary/shadow ray + inherits the path's time). Sample scenes: `samples/motionblur.usda`, `samples/curves.usda`. - Materials resolve via `MaterialBindingAPI`, dispatched on the bound shader's `info:id`: - `UsdPreviewSurface` → mapped into `OpenPBR` (portable; `diffuseColor→baseColor`, `metallic→baseMetalness`, `roughness→specularRoughness`, etc.). @@ -247,6 +290,16 @@ feature flag. ## Known incomplete work +- **Geometry/acceleration caveats.** Motion blur is transform-only and lerps the *matrix* + linearly (no deformation blur, no quaternion motion — a large shutter rotation bows + slightly, but the union-of-endpoints bbox stays conservative). Curve import flattens + cubic spans to polylines (no exact cubic intersector) and lerps widths across a span in + parameter; the rounded-cone can report an interior sphere surface for rays *starting + inside* the hull (irrelevant for opaque hair). Mesh-BVH sharing needs identical + points/topology *and* material binding — `UsdGeomPointInstancer` and `instanceable` + composition arcs are not consumed. Emissive curves/instances are not light-list entries + (BSDF-sampled only, like emissive volumes). + - **Upstream `openusd` xformOp bug, worked around locally.** `openusd` 0.5.0 (latest as of 2026-06) composes multi-op `xformOpOrder` stacks in the wrong order (the authored translate comes back multiplied by the scale), which used to make diff --git a/README.md b/README.md index 7bd5424..ea62d3c 100644 --- a/README.md +++ b/README.md @@ -66,6 +66,29 @@ cargo run --release -- -i samples/cornellbox.usda -o cornell.exr cargo run --release ``` +### 📐 Geometry & acceleration + +All intersection lives in the **`crust-rt`** kernel crate, behind an +**Embree-shaped API** (`Geometry` → `SceneBuilder` → `commit()` → +`intersect`/`occluded`, ID-based hits — swappable for Embree bindings behind +the same seam, in 100 % safe Rust). Meshes triangulate with a **watertight** +ray/triangle test (Woop et al. 2013 — no pinholes along shared edges) and +build a **BVH4**: a parallel, deterministic, reference-based SAH build with +SBVH **spatial splits** (Stich et al. 2009), collapsed into 4-wide SIMD nodes +whose slab tests run on `Vec4` lanes. Shadow rays use a dedicated early-exit +**occlusion query** (the `rtcIntersect`/`rtcOccluded` split). Each mesh's +kernel scene is built in local space and **instanced**: prims sharing +points/topology/material share one triangle BVH under an instance transform. +`UsdGeomBasisCurves` import as **round curve segments** (sphere-swept cones; +cubic bezier/bspline/catmullRom spans flatten to polylines) — see +`samples/curves.usda`. Two per-prim extras: + +- `crust:motion:translate = (x, y, z)` — **transform motion blur**: the prim + streaks through that world-space translation over the shutter + (`samples/motionblur.usda`). +- `crust:rayMask = ` — **ray visibility mask** (bit 0 camera, bit 1 + shadow, bit 2 indirect): e.g. `6` is a shadow-caster hidden from the camera. + ### 🎨 Materials Every `UsdGeomMesh` / `UsdGeomSphere` binds a `UsdShadeMaterial` via diff --git a/crates/crust-core/Cargo.toml b/crates/crust-core/Cargo.toml index c01763a..49e8833 100644 --- a/crates/crust-core/Cargo.toml +++ b/crates/crust-core/Cargo.toml @@ -11,6 +11,7 @@ path = "src/lib.rs" [dependencies] utils = { path = "../utils" } openqmc = { path = "../openqmc-rs", package = "openqmc-rs" } +crust-rt = { path = "../crust-rt" } rayon = "1.10.0" glam.workspace = true tracing.workspace = true diff --git a/crates/crust-core/src/aabb.rs b/crates/crust-core/src/aabb.rs index 5007181..4c0a1e7 100644 --- a/crates/crust-core/src/aabb.rs +++ b/crates/crust-core/src/aabb.rs @@ -1,87 +1,5 @@ -use crate::ray::Ray; -use glam::Vec3A; +//! The axis-aligned bounding box now lives in the `crust-rt` kernel +//! crate; re-exported here for the modules (guiding, volumes) that use it +//! as a plain geometric type. -#[derive(Debug, Clone, Copy)] -pub struct AABB { - pub minimum: Vec3A, - pub maximum: Vec3A, -} - -impl AABB { - pub fn new(minimum: Vec3A, maximum: Vec3A) -> Self { - AABB { minimum, maximum } - } - - pub fn surrounding_box(box0: AABB, box1: AABB) -> AABB { - let small = Vec3A::new( - box0.minimum.x.min(box1.minimum.x), - box0.minimum.y.min(box1.minimum.y), - box0.minimum.z.min(box1.minimum.z), - ); - let big = Vec3A::new( - box0.maximum.x.max(box1.maximum.x), - box0.maximum.y.max(box1.maximum.y), - box0.maximum.z.max(box1.maximum.z), - ); - AABB { - minimum: small, - maximum: big, - } - } - - pub fn hit(&self, ray: &Ray, mut t_min: f32, mut t_max: f32) -> bool { - for a in 0..3 { - let inv_d = 1.0 / ray.direction()[a]; - let mut t0 = (self.minimum[a] - ray.origin()[a]) * inv_d; - let mut t1 = (self.maximum[a] - ray.origin()[a]) * inv_d; - - if inv_d < 0.0 { - std::mem::swap(&mut t0, &mut t1); - } - - t_min = t_min.max(t0); - t_max = t_max.min(t1); - - if t_max <= t_min { - return false; - } - } - true - } - - pub fn compare_x(a: AABB, b: AABB) -> std::cmp::Ordering { - a.minimum.x.partial_cmp(&b.minimum.x).unwrap() - } - - pub fn compare_y(a: AABB, b: AABB) -> std::cmp::Ordering { - a.minimum.y.partial_cmp(&b.minimum.y).unwrap() - } - - pub fn compare_z(a: AABB, b: AABB) -> std::cmp::Ordering { - a.minimum.z.partial_cmp(&b.minimum.z).unwrap() - } -} - -pub fn triangle_aabb(v0: Vec3A, v1: Vec3A, v2: Vec3A) -> AABB { - let mut min = Vec3A::new( - v0[0].min(v1[0]).min(v2[0]), - v0[1].min(v1[1]).min(v2[1]), - v0[2].min(v1[2]).min(v2[2]), - ); - let mut max = Vec3A::new( - v0[0].max(v1[0]).max(v2[0]), - v0[1].max(v1[1]).max(v2[1]), - v0[2].max(v1[2]).max(v2[2]), - ); - // An axis-aligned triangle has zero extent along one axis, and the slab - // test in `AABB::hit` rejects zero-thickness boxes (`t_max <= t_min`). - // Pad degenerate axes so flat geometry stays hittable. - const PAD: f32 = 1e-4; - for a in 0..3 { - if max[a] - min[a] < PAD { - min[a] -= PAD; - max[a] += PAD; - } - } - AABB::new(min, max) -} +pub use crust_rt::AABB; diff --git a/crates/crust-core/src/bvh.rs b/crates/crust-core/src/bvh.rs deleted file mode 100644 index 01a4973..0000000 --- a/crates/crust-core/src/bvh.rs +++ /dev/null @@ -1,389 +0,0 @@ -//! Flattened bounding-volume hierarchy. -//! -//! One structure serves both levels of the scene: a per-mesh BVH over -//! triangles (built at USD import) and the top-level BVH over all scene -//! objects (built by `Renderer::new`). Nodes live in a contiguous array in -//! depth-first order and are traversed iteratively with a fixed stack — -//! no per-node allocation, no `Arc` pointer chasing. Splits use binned -//! surface-area-heuristic partitioning on the longest centroid axis, which -//! is deterministic: the same input always builds the same tree. - -use crate::aabb::AABB; -use crate::hittable::{Hit, Hittable}; -use crate::ray::Ray; -use glam::Vec3A; - -/// Leaves are forced at this depth so the traversal stack can never -/// overflow; SAH partitions are otherwise free to be arbitrarily uneven. -const MAX_DEPTH: usize = 60; -/// Ranges at or below this size always become a leaf. -const MIN_LEAF: usize = 2; -/// A range no larger than this may stay a leaf when splitting is not -/// worth it by SAH cost; larger ranges are always split. -const MAX_LEAF: usize = 8; -/// Number of candidate split planes tested per axis. -const BINS: usize = 12; - -struct Node { - bbox: AABB, - /// Leaf (`count > 0`): index of the first primitive in `prims`. - /// Internal (`count == 0`): index of the right child — the left child - /// immediately follows the node itself in depth-first order. - first_or_right: u32, - count: u32, - /// Split axis of an internal node, for front-to-back child ordering. - axis: u8, -} - -pub struct Bvh { - nodes: Vec, - /// Primitives permuted into leaf order, so each leaf owns a contiguous run. - prims: Vec>, - /// Objects without a bounding box cannot enter the tree and are tested - /// linearly on every ray. - unbounded: Vec>, -} - -impl Bvh { - pub fn new(objects: Vec>) -> Self { - let mut bounded = Vec::with_capacity(objects.len()); - let mut bboxes = Vec::with_capacity(objects.len()); - let mut unbounded = Vec::new(); - for obj in objects { - match obj.bounding_box() { - Some(b) => { - bboxes.push(b); - bounded.push(obj); - } - None => unbounded.push(obj), - } - } - - let centroids: Vec = bboxes - .iter() - .map(|b| 0.5 * (b.minimum + b.maximum)) - .collect(); - let mut order: Vec = (0..bounded.len() as u32).collect(); - let mut nodes = Vec::new(); - if !bounded.is_empty() { - nodes.reserve(2 * bounded.len()); - build_range(&mut nodes, &mut order, 0, bounded.len(), &bboxes, ¢roids, 0); - } - - // Permute the primitives so leaf ranges are contiguous. - let mut slots: Vec>> = bounded.into_iter().map(Some).collect(); - let prims = order - .iter() - .map(|&i| slots[i as usize].take().expect("permutation visits each slot once")) - .collect(); - - Bvh { - nodes, - prims, - unbounded, - } - } - - pub fn count(&self) -> usize { - self.prims.len() + self.unbounded.len() - } -} - -impl Hittable for Bvh { - fn hit(&self, ray: &Ray, t_min: f32, t_max: f32) -> Option> { - let mut closest = t_max; - let mut best: Option = None; - - for obj in &self.unbounded { - if let Some(hit) = obj.hit(ray, t_min, closest) { - closest = hit.rec.t; - best = Some(hit); - } - } - - if self.nodes.is_empty() { - return best; - } - - // Depth is bounded by MAX_DEPTH and each traversal level holds at - // most one deferred sibling on the stack. - let mut stack = [0u32; MAX_DEPTH + 4]; - stack[0] = 0; - let mut sp = 1usize; - - while sp > 0 { - sp -= 1; - let idx = stack[sp]; - let node = &self.nodes[idx as usize]; - if !node.bbox.hit(ray, t_min, closest) { - continue; - } - if node.count > 0 { - let first = node.first_or_right as usize; - for prim in &self.prims[first..first + node.count as usize] { - if let Some(hit) = prim.hit(ray, t_min, closest) { - closest = hit.rec.t; - best = Some(hit); - } - } - } else { - // Visit the child on the ray's near side first so its hits - // shrink `closest` before the far child is tested. - let left = idx + 1; - let right = node.first_or_right; - let (near, far) = if ray.direction()[node.axis as usize] < 0.0 { - (right, left) - } else { - (left, right) - }; - stack[sp] = far; - stack[sp + 1] = near; - sp += 2; - } - } - - best - } - - fn bounding_box(&self) -> Option { - if !self.unbounded.is_empty() { - return None; - } - self.nodes.first().map(|n| n.bbox) - } -} - -fn surface_area(b: &AABB) -> f32 { - let d = b.maximum - b.minimum; - 2.0 * (d.x * d.y + d.y * d.z + d.z * d.x) -} - -/// Builds the node for `order[start..end]` (indices into the bbox/centroid -/// arrays), appending it and its subtree to `nodes` in depth-first order. -/// Returns the node's index. Leaf `first` indices refer to positions in -/// `order`, which after the build become positions in the permuted `prims`. -fn build_range( - nodes: &mut Vec, - order: &mut [u32], - start: usize, - end: usize, - bboxes: &[AABB], - centroids: &[Vec3A], - depth: usize, -) -> u32 { - let count = end - start; - let mut bbox = bboxes[order[start] as usize]; - for &i in &order[start + 1..end] { - bbox = AABB::surrounding_box(bbox, bboxes[i as usize]); - } - - let idx = nodes.len() as u32; - let leaf = |nodes: &mut Vec| { - nodes.push(Node { - bbox, - first_or_right: start as u32, - count: count as u32, - axis: 0, - }); - idx - }; - - if count <= MIN_LEAF || depth >= MAX_DEPTH { - return leaf(nodes); - } - - // Split along the longest axis of the centroid bounds. - let mut cmin = centroids[order[start] as usize]; - let mut cmax = cmin; - for &i in &order[start + 1..end] { - cmin = cmin.min(centroids[i as usize]); - cmax = cmax.max(centroids[i as usize]); - } - let extent = cmax - cmin; - let axis = if extent.x >= extent.y && extent.x >= extent.z { - 0 - } else if extent.y >= extent.z { - 1 - } else { - 2 - }; - if extent[axis] <= 1e-6 { - // All centroids coincide — nothing to partition on. - return leaf(nodes); - } - - // Binned SAH: histogram the centroids, then score every split plane - // between adjacent bins by `area · count` on each side. - let bin_of = |i: u32| -> usize { - let t = (centroids[i as usize][axis] - cmin[axis]) / extent[axis]; - ((t * BINS as f32) as usize).min(BINS - 1) - }; - let mut bin_counts = [0usize; BINS]; - let mut bin_bounds: [Option; BINS] = [None; BINS]; - for &i in &order[start..end] { - let b = bin_of(i); - bin_counts[b] += 1; - bin_bounds[b] = Some(match bin_bounds[b] { - Some(existing) => AABB::surrounding_box(existing, bboxes[i as usize]), - None => bboxes[i as usize], - }); - } - - let side_cost = |bounds: Option, count: usize| -> f32 { - match bounds { - Some(b) if count > 0 => surface_area(&b) * count as f32, - _ => 0.0, - } - }; - let mut best: Option<(usize, f32)> = None; // (last bin of the left side, cost) - for split in 0..BINS - 1 { - let mut lb = None; - let mut lc = 0usize; - for b in 0..=split { - lc += bin_counts[b]; - lb = match (lb, bin_bounds[b]) { - (Some(x), Some(y)) => Some(AABB::surrounding_box(x, y)), - (x, y) => x.or(y), - }; - } - let mut rb = None; - let mut rc = 0usize; - for b in split + 1..BINS { - rc += bin_counts[b]; - rb = match (rb, bin_bounds[b]) { - (Some(x), Some(y)) => Some(AABB::surrounding_box(x, y)), - (x, y) => x.or(y), - }; - } - if lc == 0 || rc == 0 { - continue; - } - let cost = side_cost(lb, lc) + side_cost(rb, rc); - if best.is_none_or(|(_, c)| cost < c) { - best = Some((split, cost)); - } - } - - let mid = match best { - // Small ranges may stay a leaf when the best split costs more than - // intersecting every primitive directly. - Some((_, cost)) if count <= MAX_LEAF && cost >= surface_area(&bbox) * count as f32 => { - return leaf(nodes); - } - Some((split, _)) => { - // Unstable in-place partition by bin — deterministic for a given - // input order. - let mut mid = start; - for i in start..end { - if bin_of(order[i]) <= split { - order.swap(i, mid); - mid += 1; - } - } - mid - } - // Every populated bin on one side (can happen with extreme float - // distributions): fall back to a median split. - None => { - order[start..end] - .sort_unstable_by(|&a, &b| { - centroids[a as usize][axis].total_cmp(¢roids[b as usize][axis]) - }); - start + count / 2 - } - }; - - nodes.push(Node { - bbox, - first_or_right: 0, // patched once the left subtree is laid out - count: 0, - axis: axis as u8, - }); - build_range(nodes, order, start, mid, bboxes, centroids, depth + 1); - let right = build_range(nodes, order, mid, end, bboxes, centroids, depth + 1); - nodes[idx as usize].first_or_right = right; - idx -} - -#[cfg(test)] -mod tests { - use super::*; - use crate::material::OpenPBR; - use crate::primitives::Sphere; - use std::sync::Arc; - - fn sphere_grid(n: i32) -> Vec> { - let mat = Arc::new(OpenPBR::diffuse(Vec3A::splat(0.5))); - let mut out: Vec> = Vec::new(); - for x in 0..n { - for y in 0..n { - for z in 0..n { - out.push(Box::new(Sphere::new( - Vec3A::new(x as f32, y as f32, z as f32) * 3.0, - 0.5, - mat.clone(), - ))); - } - } - } - out - } - - /// The BVH must find exactly the hits a linear scan finds. - #[test] - fn matches_linear_scan() { - let bvh = Bvh::new(sphere_grid(4)); - let mut list = crate::hittable_list::HittableList::new(); - for obj in sphere_grid(4) { - list.add(obj); - } - - let origins = [ - Vec3A::new(-5.0, 4.5, 4.5), - Vec3A::new(20.0, 3.0, 3.0), - Vec3A::new(4.5, -5.0, 4.5), - Vec3A::new(0.0, 0.0, -10.0), - ]; - let dirs = [ - Vec3A::new(1.0, 0.0, 0.0), - Vec3A::new(-1.0, 0.05, 0.02).normalize(), - Vec3A::new(0.0, 1.0, 0.0), - Vec3A::new(0.3, 0.3, 1.0).normalize(), - Vec3A::new(0.0, 0.0, -1.0), - ]; - for o in origins { - for d in dirs { - let ray = Ray::new(o, d); - let a = bvh.hit(&ray, 0.001, f32::INFINITY); - let b = list.hit(&ray, 0.001, f32::INFINITY); - match (a, b) { - (Some(x), Some(y)) => { - assert!((x.rec.t - y.rec.t).abs() < 1e-4, "t mismatch for {o:?} {d:?}"); - } - (None, None) => {} - (x, y) => panic!( - "hit disagreement for {o:?} {d:?}: bvh={} linear={}", - x.is_some(), - y.is_some() - ), - } - } - } - } - - #[test] - fn empty_bvh_misses() { - let bvh = Bvh::new(Vec::new()); - let ray = Ray::new(Vec3A::ZERO, Vec3A::X); - assert!(bvh.hit(&ray, 0.001, f32::INFINITY).is_none()); - assert!(bvh.bounding_box().is_none()); - } - - #[test] - fn bounding_box_covers_all_prims() { - let bvh = Bvh::new(sphere_grid(3)); - let bbox = bvh.bounding_box().expect("grid is fully bounded"); - assert!(bbox.minimum.cmple(Vec3A::splat(-0.5)).all()); - assert!(bbox.maximum.cmpge(Vec3A::splat(6.5)).all()); - } -} diff --git a/crates/crust-core/src/camera.rs b/crates/crust-core/src/camera.rs index ab40b97..87e1f4b 100644 --- a/crates/crust-core/src/camera.rs +++ b/crates/crust-core/src/camera.rs @@ -66,7 +66,9 @@ impl Camera { /// - `lens_uv`: A 2D uniform sample used to sample the lens for depth of /// field. Ignored when the camera has zero aperture. Callers pass a QMC /// sample so this dimension is decorrelated from the pixel jitter. - pub fn get_ray(&self, s: f32, t: f32, lens_uv: [f32; 2]) -> Ray { + /// - `time`: Shutter time in `[0, 1)`, carried on the ray for motion + /// blur (moving instances interpolate their transform at this time). + pub fn get_ray(&self, s: f32, t: f32, lens_uv: [f32; 2], time: f32) -> Ray { let offset = if self.lens_radius > 0.0 { let rd = self.lens_radius * concentric_disk(lens_uv); self.u * rd.x + self.v * rd.y @@ -77,5 +79,7 @@ impl Camera { self.origin + offset, self.lower_left_corner + s * self.horizontal + t * self.vertical - self.origin - offset, ) + .with_time(time) + .with_mask(crate::ray::MASK_CAMERA) } } diff --git a/crates/crust-core/src/hittable.rs b/crates/crust-core/src/hittable.rs index f2537a5..7a4dd04 100644 --- a/crates/crust-core/src/hittable.rs +++ b/crates/crust-core/src/hittable.rs @@ -1,13 +1,12 @@ use crate::ray::Ray; use glam::Vec3A; -use crate::aabb::AABB; -use crate::material::Material; - -/// The `HitRecord` struct stores the geometry of a ray-object intersection: -/// the intersection point, surface normal, ray parameter, and facing. -/// It is plain `Copy` data — the material of the hit surface travels -/// alongside it in [`Hit`], borrowed from the primitive that was hit. +/// The `HitRecord` struct stores the geometry of a ray-surface +/// intersection as the *materials* consume it: the intersection point, +/// the ray-facing surface normal, the ray parameter, and the facing flag. +/// Intersection itself happens in the `crust-rt` kernel (which reports +/// ID-based [`crust_rt::RayHit`]s); [`crate::World`] converts kernel hits +/// into `HitRecord`s and looks up the material by `geom_id`. #[derive(Clone, Copy, Default)] pub struct HitRecord { /// The point of intersection. @@ -22,9 +21,6 @@ pub struct HitRecord { impl HitRecord { /// Creates a new, default `HitRecord`. - /// - /// # Returns - /// - A new instance of `HitRecord` with default values. pub fn new() -> HitRecord { Default::default() } @@ -46,31 +42,3 @@ impl HitRecord { }; } } - -/// A successful ray-object intersection: the geometric [`HitRecord`] plus the -/// material at the hit point, borrowed from the scene for the duration of the -/// trace. Borrowing (instead of the former `Option>` field) -/// keeps intersection records `Copy` and removes an atomic refcount bump per -/// candidate hit in the traversal hot loop. -#[derive(Clone, Copy)] -pub struct Hit<'a> { - pub rec: HitRecord, - pub mat: &'a dyn Material, -} - -/// The `Hittable` trait defines objects that can be intersected by rays. -/// Implementing this trait allows objects to participate in ray tracing. -pub trait Hittable: Send + Sync { - /// Determines if a ray intersects the object. - /// - /// # Parameters - /// - `ray`: The ray to test for intersection. - /// - `t_min`: The minimum value of the parameter `t` to consider. - /// - `t_max`: The maximum value of the parameter `t` to consider. - /// - /// # Returns - /// - `Some(Hit)` describing the closest intersection in `(t_min, t_max)`, - /// or `None` if the ray misses the object. - fn hit(&self, ray: &Ray, t_min: f32, t_max: f32) -> Option>; - fn bounding_box(&self) -> Option; -} diff --git a/crates/crust-core/src/hittable_list.rs b/crates/crust-core/src/hittable_list.rs deleted file mode 100644 index 9c6c176..0000000 --- a/crates/crust-core/src/hittable_list.rs +++ /dev/null @@ -1,83 +0,0 @@ -use crate::aabb::AABB; -use crate::hittable::{Hit, Hittable}; -use crate::ray::Ray; - -/// The `HittableList` struct represents a collection of objects that can be intersected by rays. -/// It allows for managing multiple `Hittable` objects and testing for ray intersections with all of them. -/// -/// This is the scene-construction container; for rendering, `Renderer::new` -/// converts it into a top-level [`crate::Bvh`] via [`HittableList::into_objects`]. -#[derive(Default)] -pub struct HittableList { - /// A vector of objects implementing the `Hittable` trait. - objects: Vec>, -} - -impl HittableList { - /// Creates a new, empty `HittableList`. - /// - /// # Returns - /// - A new instance of `HittableList`. - pub fn new() -> HittableList { - Default::default() - } - - /// Adds a `Hittable` object to the list. - /// - /// # Parameters - /// - `object`: A boxed object implementing the `Hittable` trait. - pub fn add(&mut self, object: Box) { - self.objects.push(object); - } - /// Returns the count of objects in the list. - /// - /// # Returns - /// - The number of objects in the list. - pub fn count(&self) -> usize { - self.objects.len() - } - - /// Consumes the list and returns its objects, e.g. to build an - /// acceleration structure over them. - pub fn into_objects(self) -> Vec> { - self.objects - } -} - -impl Hittable for HittableList { - /// Finds the closest intersection of the ray with any object in the list - /// by linear scan. - fn hit(&self, ray: &Ray, t_min: f32, t_max: f32) -> Option> { - let mut closest_so_far = t_max; - let mut best: Option = None; - - for object in &self.objects { - if let Some(hit) = object.hit(ray, t_min, closest_so_far) { - closest_so_far = hit.rec.t; - best = Some(hit); - } - } - - best - } - fn bounding_box(&self) -> Option { - if self.objects.is_empty() { - return None; - } - - let mut temp_box: Option = None; - - for object in &self.objects { - if let Some(bbox) = object.bounding_box() { - temp_box = Some(match temp_box { - Some(existing) => AABB::surrounding_box(existing, bbox), - None => bbox, - }); - } else { - return None; // fail if any object doesn't have a bounding box - } - } - - temp_box - } -} diff --git a/crates/crust-core/src/lib.rs b/crates/crust-core/src/lib.rs index 374800a..8154c75 100644 --- a/crates/crust-core/src/lib.rs +++ b/crates/crust-core/src/lib.rs @@ -1,16 +1,14 @@ mod aabb; mod buffer; -mod bvh; mod camera; mod error; mod guiding; mod hittable; -mod hittable_list; mod light; mod material; mod medium; -mod primitives; mod ray; +mod rt_world; mod scene; mod tracer; mod volume; @@ -21,20 +19,22 @@ mod world; /// edit can swap in another OpenQMC sampler (e.g. `SobolBnSampler`). pub type PathSampler = openqmc::SobolSampler; +/// The intersection kernel (Embree-shaped scene/geometry API), re-exported +/// so applications can build [`rt::Geometry`] values for [`WorldBuilder`]. +pub use crust_rt as rt; + pub use aabb::AABB; pub use buffer::Buffer; -pub use bvh::Bvh; pub use camera::Camera; pub use error::Error; pub use glam::{Mat4, Vec3A}; pub use guiding::{GuidingConfig, GuidingField, SampleData}; -pub use hittable::{Hit, HitRecord, Hittable}; -pub use hittable_list::HittableList; +pub use hittable::HitRecord; pub use light::{AreaLight, Light, LightList, LightShape, RectShape, SphereShape}; pub use material::*; pub use medium::Medium; -pub use primitives::{SmoothTriangle, Sphere, Triangle}; -pub use ray::Ray; +pub use ray::{MASK_ALL, MASK_CAMERA, MASK_INDIRECT, MASK_SHADOW, Ray}; +pub use rt_world::{World, WorldBuilder, WorldHit}; pub use scene::Scene; pub use tracer::{ProgressCallback, RenderSettings, Renderer, SamplingStrategy, ray_color}; pub use volume::{DensityField, PhaseMix, VolumeEvent, VolumeRegion, Volumes}; diff --git a/crates/crust-core/src/light.rs b/crates/crust-core/src/light.rs index 7b980f5..69edd3a 100644 --- a/crates/crust-core/src/light.rs +++ b/crates/crust-core/src/light.rs @@ -82,8 +82,8 @@ impl LightShape for RectShape { /// The `Light` trait is what the integrator's light-sampling strategy (NEE) /// needs from a light: a surface point to aim a shadow ray at, the /// solid-angle density of that choice for MIS, the emitted radiance, and — -/// for lights with scene geometry — the material identity that lets a -/// bounce ray recognize the light it hit. +/// for lights with scene geometry — the geometry id that lets a bounce +/// ray recognize the light it hit. pub trait Light: Send + Sync { /// Samples a point on the light source, uniform by area. /// @@ -99,10 +99,10 @@ pub trait Light: Send + Sync { /// material bound to the light's scene geometry. fn emission(&self) -> Vec3A; - /// The emissive material bound to this light's scene geometry, used to - /// recognize the light when a bounce ray hits it (by address identity). - /// `None` for lights with no geometry in the world. - fn material(&self) -> Option<&dyn Material>; + /// The `geom_id` of this light's scene geometry in the world, used to + /// recognize the light when a bounce ray hits it. `None` for lights + /// with no geometry in the world. + fn geom_id(&self) -> Option; } /// A geometric area light: any [`LightShape`] paired with the [`Emissive`] @@ -111,11 +111,18 @@ pub trait Light: Send + Sync { pub struct AreaLight { shape: Box, material: Arc, + /// The world `geom_id` of the emissive geometry this light shares its + /// surface with — how bounce hits are attributed back to the light. + geom_id: u32, } impl AreaLight { - pub fn new(shape: Box, material: Arc) -> Self { - Self { shape, material } + pub fn new(shape: Box, material: Arc, geom_id: u32) -> Self { + Self { + shape, + material, + geom_id, + } } } @@ -144,8 +151,8 @@ impl Light for AreaLight { self.material.emitted() } - fn material(&self) -> Option<&dyn Material> { - Some(self.material.as_ref()) + fn geom_id(&self) -> Option { + Some(self.geom_id) } } @@ -190,15 +197,14 @@ impl LightList { } } - /// Finds the light whose scene geometry carries `mat`, by address - /// identity of the material. Used by the integrator to attribute a - /// bounce-hit emissive surface to its light for MIS; emissive geometry - /// with no light-list entry returns `None`. - pub fn find_by_material(&self, mat: &dyn Material) -> Option<&Arc> { - self.lights.iter().find(|l| { - l.material() - .is_some_and(|m| std::ptr::addr_eq(m as *const dyn Material, mat)) - }) + /// Finds the light whose scene geometry has world id `geom_id`. Used + /// by the integrator to attribute a bounce-hit emissive surface to its + /// light for MIS; emissive geometry with no light-list entry returns + /// `None`. + pub fn find_by_geom(&self, geom_id: u32) -> Option<&Arc> { + self.lights + .iter() + .find(|l| l.geom_id() == Some(geom_id)) } /// Returns the number of lights in the `LightList`. @@ -248,6 +254,7 @@ mod tests { radius: 1.0, }), Arc::new(Emissive::new(Vec3A::splat(10.0))), + 0, ); // Nearest point on the sphere as seen from below. let pdf = light.pdf(Vec3A::ZERO, Vec3A::new(0.0, 4.0, 0.0)); @@ -256,20 +263,19 @@ mod tests { } #[test] - fn find_by_material_matches_by_address() { - let mat_a = Arc::new(Emissive::new(Vec3A::splat(1.0))); - let mat_b = Arc::new(Emissive::new(Vec3A::splat(1.0))); + fn find_by_geom_matches_by_id() { + let mat = Arc::new(Emissive::new(Vec3A::splat(1.0))); let mut lights = LightList::new(); lights.add(Arc::new(AreaLight::new( Box::new(SphereShape { center: Vec3A::ZERO, radius: 1.0, }), - mat_a.clone(), + mat, + 7, ))); - assert!(lights.find_by_material(mat_a.as_ref()).is_some()); - // Equal color but a different allocation must not match. - assert!(lights.find_by_material(mat_b.as_ref()).is_none()); + assert!(lights.find_by_geom(7).is_some()); + assert!(lights.find_by_geom(8).is_none()); } } diff --git a/crates/crust-core/src/primitives/mesh.rs b/crates/crust-core/src/primitives/mesh.rs deleted file mode 100644 index ae0f860..0000000 --- a/crates/crust-core/src/primitives/mesh.rs +++ /dev/null @@ -1,85 +0,0 @@ -use crate::aabb::AABB; -use crate::hittable::{Hit, HitRecord, Hittable}; -use crate::material::Material; -use crate::primitives::triangle::triangle_hit; -use crate::ray::Ray; -use glam::Vec3A; -use std::sync::Arc; - -pub struct Mesh { - vertices: Vec, - indices: Vec, - pub material: Arc, -} -// `Mesh` is not re-exported at the crate root yet, so its constructors and -// accessors count as dead code from inside the crate. -#[allow(dead_code)] -impl Mesh { - pub fn new(vertices: Vec, indices: Vec, material: Arc) -> Self { - Self { - vertices, - indices, - material, - } - } - - pub fn get_vertices(&self) -> &Vec { - &self.vertices - } - - pub fn get_indices(&self) -> &Vec { - &self.indices - } -} - -impl Hittable for Mesh { - fn bounding_box(&self) -> Option { - let mut verts = self.vertices.iter(); - let first = *verts.next()?; - let (mut min, mut max) = (first, first); - for &v in verts { - min = min.min(v); - max = max.max(v); - } - // Pad degenerate axes so flat meshes survive the slab test, matching - // `triangle_aabb`. - const PAD: f32 = 1e-4; - for a in 0..3 { - if max[a] - min[a] < PAD { - min[a] -= PAD; - max[a] += PAD; - } - } - Some(AABB::new(min, max)) - } - fn hit(&self, r: &Ray, t_min: f32, t_max: f32) -> Option> { - indexed_mesh_hit(r, &self.vertices, &self.indices, t_min, t_max).map(|rec| Hit { - rec, - mat: self.material.as_ref(), - }) - } -} - -fn indexed_mesh_hit( - ray: &Ray, - vertices: &[Vec3A], - indices: &[u32], - t_min: f32, - t_max: f32, -) -> Option { - let mut best: Option = None; - let mut closest_so_far = t_max; - - for tri in indices.chunks_exact(3) { - let v0 = vertices[tri[0] as usize]; - let v1 = vertices[tri[1] as usize]; - let v2 = vertices[tri[2] as usize]; - - if let Some(rec) = triangle_hit(ray, v0, v1, v2, t_min, closest_so_far) { - closest_so_far = rec.t; - best = Some(rec); - } - } - - best -} diff --git a/crates/crust-core/src/primitives/mod.rs b/crates/crust-core/src/primitives/mod.rs deleted file mode 100644 index 5b25878..0000000 --- a/crates/crust-core/src/primitives/mod.rs +++ /dev/null @@ -1,9 +0,0 @@ -mod sphere; -pub use sphere::Sphere; -mod triangle; -pub use triangle::Triangle; -mod mesh; -#[allow(unused_imports)] -pub use mesh::Mesh; -mod smooth_triangle; -pub use smooth_triangle::SmoothTriangle; diff --git a/crates/crust-core/src/primitives/smooth_triangle.rs b/crates/crust-core/src/primitives/smooth_triangle.rs deleted file mode 100644 index b59c778..0000000 --- a/crates/crust-core/src/primitives/smooth_triangle.rs +++ /dev/null @@ -1,91 +0,0 @@ -use crate::aabb::{AABB, triangle_aabb}; -use crate::hittable::{Hit, HitRecord, Hittable}; -use crate::material::Material; -use crate::ray::Ray; -use glam::Vec3A; -use std::sync::Arc; - -pub struct SmoothTriangle { - pub v0: Vec3A, - pub v1: Vec3A, - pub v2: Vec3A, - pub n0: Vec3A, - pub n1: Vec3A, - pub n2: Vec3A, - pub material: Arc, -} - -impl SmoothTriangle { - pub fn new( - v0: Vec3A, - v1: Vec3A, - v2: Vec3A, - n0: Vec3A, - n1: Vec3A, - n2: Vec3A, - material: Arc, - ) -> Self { - Self { - v0, - v1, - v2, - n0, - n1, - n2, - material, - } - } -} - -impl Hittable for SmoothTriangle { - fn hit(&self, ray: &Ray, t_min: f32, t_max: f32) -> Option> { - let edge1 = self.v1 - self.v0; - let edge2 = self.v2 - self.v0; - let h = ray.direction().cross(edge2); - let a = edge1.dot(h); - - // Use a relative epsilon based on triangle size - let triangle_size = (edge1.length() + edge2.length()) * 0.5; - let epsilon = triangle_size * 1e-6; - - if a.abs() < epsilon { - return None; - } - - let f = 1.0 / a; - let s = ray.origin() - self.v0; - let u = f * s.dot(h); - if !(0.0..=1.0).contains(&u) { - return None; - } - - let q = s.cross(edge1); - let v = f * ray.direction().dot(q); - if v < 0.0 || u + v > 1.0 { - return None; - } - - let t = f * edge2.dot(q); - if t < t_min || t > t_max { - return None; - } - - let mut rec = HitRecord::new(); - rec.t = t; - rec.p = ray.at(t); - - // Interpolate normal using barycentric weights - let w = 1.0 - u - v; - let interpolated_normal = (self.n0 * w + self.n1 * u + self.n2 * v).normalize(); - rec.set_face_normal(ray, interpolated_normal); - - Some(Hit { - rec, - mat: self.material.as_ref(), - }) - } - - fn bounding_box(&self) -> Option { - Some(triangle_aabb(self.v0, self.v1, self.v2)) - } -} diff --git a/crates/crust-core/src/primitives/sphere.rs b/crates/crust-core/src/primitives/sphere.rs deleted file mode 100644 index 9f0f27a..0000000 --- a/crates/crust-core/src/primitives/sphere.rs +++ /dev/null @@ -1,60 +0,0 @@ -use crate::aabb::AABB; -use crate::hittable::{Hit, HitRecord, Hittable}; -use crate::material::Material; -use crate::ray::Ray; -use glam::Vec3A; -use std::sync::Arc; - -pub struct Sphere { - center: Vec3A, - radius: f32, - material: Arc, -} -impl Sphere { - pub fn new(center: Vec3A, radius: f32, material: Arc) -> Self { - Self { - center, - radius, - material, - } - } -} -impl Hittable for Sphere { - fn bounding_box(&self) -> Option { - Some(AABB::new( - self.center - Vec3A::new(self.radius, self.radius, self.radius), - self.center + Vec3A::new(self.radius, self.radius, self.radius), - )) - } - fn hit(&self, r: &Ray, t_min: f32, t_max: f32) -> Option> { - let oc = r.origin() - self.center; - let a = r.direction().length_squared(); - let half_b = oc.dot(r.direction()); - let c = oc.length_squared() - self.radius * self.radius; - let discriminant = half_b * half_b - a * c; - - if discriminant < 0.0 { - return None; - } - - let sqrt_d = discriminant.sqrt(); - let mut root = (-half_b - sqrt_d) / a; - - if root <= t_min || root >= t_max { - root = (-half_b + sqrt_d) / a; - if root <= t_min || root >= t_max { - return None; - } - } - - let mut rec = HitRecord::new(); - rec.t = root; - rec.p = r.at(root); - let outward_normal = (rec.p - self.center) / self.radius; - rec.set_face_normal(r, outward_normal); - Some(Hit { - rec, - mat: self.material.as_ref(), - }) - } -} diff --git a/crates/crust-core/src/primitives/triangle.rs b/crates/crust-core/src/primitives/triangle.rs deleted file mode 100644 index 5675752..0000000 --- a/crates/crust-core/src/primitives/triangle.rs +++ /dev/null @@ -1,82 +0,0 @@ -use crate::aabb::{AABB, triangle_aabb}; -use crate::hittable::{Hit, HitRecord, Hittable}; -use crate::material::Material; -use crate::ray::Ray; -use glam::Vec3A; -use std::sync::Arc; - -pub struct Triangle { - v0: Vec3A, - v1: Vec3A, - v2: Vec3A, - material: Arc, -} -impl Triangle { - pub fn new(v0: Vec3A, v1: Vec3A, v2: Vec3A, material: Arc) -> Self { - Self { - v0, - v1, - v2, - material, - } - } -} - -impl Hittable for Triangle { - fn bounding_box(&self) -> Option { - Some(triangle_aabb(self.v0, self.v1, self.v2)) - } - fn hit(&self, r: &Ray, t_min: f32, t_max: f32) -> Option> { - triangle_hit(r, self.v0, self.v1, self.v2, t_min, t_max).map(|rec| Hit { - rec, - mat: self.material.as_ref(), - }) - } -} - -pub(crate) fn triangle_hit( - ray: &Ray, - v0: Vec3A, - v1: Vec3A, - v2: Vec3A, - t_min: f32, - t_max: f32, -) -> Option { - let edge1 = v1 - v0; - let edge2 = v2 - v0; - let h = ray.direction().cross(edge2); - let a = edge1.dot(h); - - // Use a relative epsilon based on triangle size - let triangle_size = (edge1.length() + edge2.length()) * 0.5; - let epsilon = triangle_size * 1e-6; - - if a.abs() < epsilon { - return None; - } - - let f = 1.0 / a; - let s = ray.origin() - v0; - let u = f * s.dot(h); - if !(0.0..=1.0).contains(&u) { - return None; - } - - let q = s.cross(edge1); - let v = f * ray.direction().dot(q); - if v < 0.0 || u + v > 1.0 { - return None; - } - - let t = f * edge2.dot(q); - if t < t_min || t > t_max { - return None; - } - - let mut rec = HitRecord::new(); - rec.t = t; - rec.p = ray.at(t); - let normal = edge1.cross(edge2).normalize(); - rec.set_face_normal(ray, normal); - Some(rec) -} diff --git a/crates/crust-core/src/ray.rs b/crates/crust-core/src/ray.rs index d2e3661..f0c5cdf 100644 --- a/crates/crust-core/src/ray.rs +++ b/crates/crust-core/src/ray.rs @@ -2,26 +2,26 @@ use crate::medium::Medium; use glam::Vec3A; use std::sync::Arc; -/// The `Ray` struct represents a ray in 3D space, defined by an origin and a direction. -/// Rays are used in ray tracing to determine intersections with objects in the scene. -/// -/// An optional `medium` describes the participating medium the ray is -/// currently travelling through — used by transmissive OpenPBR materials so -/// the tracer can apply Beer-Lambert attenuation between surface hits. +pub use crust_rt::{MASK_ALL, MASK_CAMERA, MASK_INDIRECT, MASK_SHADOW}; + +/// The renderer's ray: the kernel ray (origin, direction, shutter time, +/// visibility mask — see [`crust_rt::Ray`]) plus the renderer-side state +/// the kernel deliberately does not know about: the participating +/// `medium` the ray is currently travelling through, used by transmissive +/// OpenPBR materials so the tracer can apply Beer-Lambert attenuation +/// between surface hits. #[derive(Default, Clone)] pub struct Ray { - orig: Vec3A, - dir: Vec3A, + rt: crust_rt::Ray, medium: Option>, } impl Ray { /// Creates a new `Ray` with the specified origin and direction, in - /// vacuum (no medium). + /// vacuum (no medium), at shutter time 0, visible to all geometry. pub fn new(origin: Vec3A, direction: Vec3A) -> Ray { Ray { - orig: origin, - dir: direction, + rt: crust_rt::Ray::new(origin, direction), medium: None, } } @@ -30,25 +30,49 @@ impl Ray { /// a refraction that enters a transmissive volume. pub fn new_in_medium(origin: Vec3A, direction: Vec3A, medium: Arc) -> Ray { Ray { - orig: origin, - dir: direction, + rt: crust_rt::Ray::new(origin, direction), medium: Some(medium), } } + /// Same ray with the shutter time replaced. + pub fn with_time(mut self, time: f32) -> Ray { + self.rt.time = time; + self + } + + /// Same ray with the visibility mask replaced. + pub fn with_mask(mut self, mask: u32) -> Ray { + self.rt.mask = mask; + self + } + + /// The kernel view of this ray — what `crust_rt` queries take. + pub fn rt(&self) -> &crust_rt::Ray { + &self.rt + } + pub fn origin(&self) -> Vec3A { - self.orig + self.rt.origin } pub fn direction(&self) -> Vec3A { - self.dir + self.rt.dir } pub fn medium(&self) -> Option<&Arc> { self.medium.as_ref() } + pub fn time(&self) -> f32 { + self.rt.time + } + + pub fn mask(&self) -> u32 { + self.rt.mask + } + pub fn at(&self, t: f32) -> Vec3A { - self.orig + t * self.dir + self.rt.at(t) } } diff --git a/crates/crust-core/src/rt_world.rs b/crates/crust-core/src/rt_world.rs new file mode 100644 index 0000000..202b403 --- /dev/null +++ b/crates/crust-core/src/rt_world.rs @@ -0,0 +1,112 @@ +//! The bridge between the `crust-rt` kernel and the renderer: a +//! [`WorldBuilder`] pairs every attached kernel [`Geometry`] with its +//! [`Material`], and the committed [`World`] resolves kernel hits +//! (`geom_id`) back to materials — the ID-based split that keeps shading +//! state out of the intersection kernel. + +use crate::hittable::HitRecord; +use crate::material::Material; +use crate::ray::Ray; +use crust_rt::{AABB, Geometry, MASK_ALL, SceneBuilder}; +use std::sync::Arc; + +/// Scene-construction container: kernel geometries plus the material +/// table indexed by `geom_id`. Commit into a [`World`] to render. +#[derive(Default)] +pub struct WorldBuilder { + rt: SceneBuilder, + materials: Vec>, +} + +impl WorldBuilder { + pub fn new() -> Self { + Default::default() + } + + /// Attaches a geometry with its material; returns the `geom_id` that + /// hits on this geometry will carry. + pub fn attach(&mut self, geometry: Geometry, material: Arc) -> u32 { + self.attach_masked(geometry, material, MASK_ALL) + } + + /// Attaches a geometry visible only to the ray categories in `mask` + /// (see the `MASK_*` constants). + pub fn attach_masked( + &mut self, + geometry: Geometry, + material: Arc, + mask: u32, + ) -> u32 { + let id = self.rt.attach_masked(geometry, mask); + self.materials.push(material); + debug_assert_eq!(id as usize + 1, self.materials.len()); + id + } + + /// Number of geometries attached so far. + pub fn count(&self) -> usize { + self.materials.len() + } + + /// Builds the acceleration structure (parallel, deterministic). + pub fn commit(self) -> World { + World { + scene: self.rt.commit(), + materials: self.materials, + } + } +} + +/// A successful world intersection: the material-facing [`HitRecord`], +/// the material looked up from the hit's `geom_id`, and the IDs +/// themselves (the integrator attributes bounce-hit lights by `geom_id`). +pub struct WorldHit<'a> { + pub rec: HitRecord, + pub mat: &'a dyn Material, + pub geom_id: u32, + pub prim_id: u32, +} + +/// The committed scene geometry the renderer traces against. +pub struct World { + scene: crust_rt::Scene, + materials: Vec>, +} + +impl World { + /// Closest hit in `(t_min, t_max)` with its material resolved. + pub fn intersect(&self, ray: &Ray, t_min: f32, t_max: f32) -> Option> { + let h = self.scene.intersect(ray.rt(), t_min, t_max)?; + Some(WorldHit { + rec: HitRecord { + p: ray.at(h.t), + normal: h.normal, + t: h.t, + front_face: h.front_face, + }, + mat: self.materials[h.geom_id as usize].as_ref(), + geom_id: h.geom_id, + prim_id: h.prim_id, + }) + } + + /// Early-exit occlusion query — the shadow-ray fast path. + pub fn occluded(&self, ray: &Ray, t_min: f32, t_max: f32) -> bool { + self.scene.occluded(ray.rt(), t_min, t_max) + } + + /// World bounds of all geometry; `None` for an empty world. + pub fn bounds(&self) -> Option { + self.scene.bounds() + } + + /// Number of geometries (`geom_id`s are `0..count`). + pub fn count(&self) -> usize { + self.materials.len() + } + + /// The material bound to a geometry. + pub fn material(&self, geom_id: u32) -> &dyn Material { + self.materials[geom_id as usize].as_ref() + } +} diff --git a/crates/crust-core/src/scene.rs b/crates/crust-core/src/scene.rs index 777cd5f..d0cb67c 100644 --- a/crates/crust-core/src/scene.rs +++ b/crates/crust-core/src/scene.rs @@ -1,16 +1,17 @@ use crate::camera::Camera; -use crate::hittable_list::HittableList; use crate::light::LightList; +use crate::rt_world::World; use crate::tracer::RenderSettings; use crate::volume::VolumeRegion; /// The renderer's runtime scene, produced from a USD stage /// (`Scene::from_usd`) or assembled by hand (`Scene::new`, e.g. from the /// procedural `world::simple_scene`). `Renderer::new` consumes this -/// directly. +/// directly. `world` is a committed [`World`] — kernel scene plus the +/// per-geometry material table. pub struct Scene { pub camera: Camera, - pub world: HittableList, + pub world: World, pub lights: LightList, pub settings: RenderSettings, /// Participating-media regions (smoke, fog, …), kept outside `world` @@ -21,7 +22,7 @@ pub struct Scene { impl Scene { pub fn new( camera: Camera, - world: HittableList, + world: World, lights: LightList, settings: RenderSettings, ) -> Self { diff --git a/crates/crust-core/src/scene/usd_import.rs b/crates/crust-core/src/scene/usd_import.rs index d810b6d..4c6bc38 100644 --- a/crates/crust-core/src/scene/usd_import.rs +++ b/crates/crust-core/src/scene/usd_import.rs @@ -1,27 +1,29 @@ //! USD scene import: opens a stage and produces a runtime `Scene` //! (camera, world, lights, render settings). See `Scene::from_usd`. +use std::collections::HashMap; +use std::hash::{Hash, Hasher}; use std::path::Path; use std::sync::Arc; use glam::Mat4 as GMat4; use tracing::{debug, info, warn}; -use crate::bvh::Bvh; use crate::camera::Camera; -use crate::hittable::Hittable; -use crate::hittable_list::HittableList; use crate::light::{AreaLight, LightList, RectShape, SphereShape}; use crate::material::{Emissive, Material, OpenPBR}; -use crate::primitives::{Sphere as CrustSphere, Triangle}; +use crate::ray::MASK_ALL; +use crate::rt_world::WorldBuilder; use crate::scene::Scene; +use crust_rt::{CurveSegment, Geometry, Scene as RtScene, SceneBuilder as RtSceneBuilder}; use crate::tracer::{RenderSettings, SamplingStrategy}; use crate::volume::{DensityField, VolumeRegion}; -use glam::{Vec3, Vec3A}; +use glam::{Affine3A, Vec3, Vec3A}; use openusd::gf::{Matrix4d, Vec3f}; use openusd::schemas::geom::{ - Camera as UsdCamera, Mesh as UsdMesh, PointBased, Sphere as UsdSphere, Xform, Xformable, + BasisCurves as UsdBasisCurves, Camera as UsdCamera, Curves as UsdCurves, Mesh as UsdMesh, + PointBased, Sphere as UsdSphere, Xform, Xformable, }; use openusd::schemas::lux::{ CylinderLight, DiskLight, DistantLight, DomeLight, Light as UsdLight, RectLight, SphereLight, @@ -54,10 +56,16 @@ pub(crate) fn load_scene(path: &Path) -> Result { // Render settings come first — the camera importer needs the aspect ratio. let settings = import_render_settings(&stage); - let mut world = HittableList::new(); + let mut world = WorldBuilder::new(); let mut lights = LightList::new(); let mut volumes: Vec = Vec::new(); let mut camera_candidate: Option = None; + // Prims binding the same material path share one Arc, and prims with + // identical local geometry + material share one committed kernel + // scene through instancing — N placements of a mesh cost one copy of + // its triangles. + let mut materials = MaterialCache::default(); + let mut meshes: HashMap> = HashMap::new(); let mut stack: Vec<(Prim, GMat4)> = vec![(stage.prim_at(sdf::Path::abs_root()), GMat4::IDENTITY)]; @@ -75,11 +83,14 @@ pub(crate) fn load_scene(path: &Path) -> Result { if custom_token(&prim, "crust:volume:type").is_some() { emit_volume(&prim, this_world, &mut volumes); } else if let Ok(Some(mesh)) = UsdMesh::get(&stage, prim.path().clone()) { - let mat = resolve_material(&stage, &prim); - emit_mesh(&mut world, &prim, &mesh, this_world, mat); + let mat = resolve_material(&stage, &prim, &mut materials); + emit_mesh(&mut world, &prim, &mesh, this_world, mat, &mut meshes); } else if let Ok(Some(sphere)) = UsdSphere::get(&stage, prim.path().clone()) { - let mat = resolve_material(&stage, &prim); + let mat = resolve_material(&stage, &prim, &mut materials); emit_sphere(&mut world, &prim, &sphere, this_world, mat); + } else if let Ok(Some(curves)) = UsdBasisCurves::get(&stage, prim.path().clone()) { + let mat = resolve_material(&stage, &prim, &mut materials); + emit_curves(&mut world, &prim, &curves, this_world, mat); } else if UsdCamera::get(&stage, prim.path().clone()) .ok() .flatten() @@ -116,7 +127,7 @@ pub(crate) fn load_scene(path: &Path) -> Result { crate::world::get_settings().0 }); - Ok(Scene::new(camera, world, lights, settings).with_volumes(volumes)) + Ok(Scene::new(camera, world.commit(), lights, settings).with_volumes(volumes)) } // ----------------------------------------------------------------------- @@ -421,16 +432,71 @@ fn resets_xform_stack_at(stage: &Stage, prim: &Prim) -> bool { false } +// ----------------------------------------------------------------------- +// Per-prim geometry attributes (visibility mask, motion) +// ----------------------------------------------------------------------- + +/// `crust:rayMask` — which ray categories see this geometry (bit 0 camera, +/// bit 1 shadow, bit 2 indirect; default: all). E.g. `crust:rayMask = 6` +/// makes a light-blocker invisible to the camera. +fn prim_ray_mask(prim: &Prim) -> u32 { + custom_i32(prim, "crust:rayMask") + .map(|m| m as u32) + .unwrap_or(MASK_ALL) +} + +/// `crust:motion:translate` — a world-space translation the prim moves +/// through over the shutter interval (transform motion blur). +fn prim_motion_translate(prim: &Prim) -> Option { + custom_color3(prim, "crust:motion:translate").map(|v| Vec3::new(v.x, v.y, v.z)) +} + // ----------------------------------------------------------------------- // Mesh // ----------------------------------------------------------------------- +/// Identity of an imported mesh's shared geometry: a content hash of the +/// authored points/counts/indices plus the (memoized, so pointer-comparable) +/// material. Prims agreeing on all of it share one local-space triangle BVH. +#[derive(PartialEq, Eq, Hash)] +struct MeshKey { + geo_hash: u64, + n_points: usize, + n_indices: usize, + material: usize, +} + +impl MeshKey { + fn new( + points: &[Vec3f], + counts: &[i32], + indices: &[i32], + material: &Arc, + ) -> Self { + let mut h = std::collections::hash_map::DefaultHasher::new(); + for p in points { + p.x.to_bits().hash(&mut h); + p.y.to_bits().hash(&mut h); + p.z.to_bits().hash(&mut h); + } + counts.hash(&mut h); + indices.hash(&mut h); + MeshKey { + geo_hash: h.finish(), + n_points: points.len(), + n_indices: indices.len(), + material: Arc::as_ptr(material) as *const u8 as usize, + } + } +} + fn emit_mesh( - world: &mut HittableList, + world: &mut WorldBuilder, prim: &Prim, mesh: &UsdMesh, world_xf: GMat4, material: Arc, + meshes: &mut HashMap>, ) { let points: Option> = mesh .points_attr() @@ -471,45 +537,112 @@ fn emit_mesh( } }; - let verts: Vec = points - .iter() - .map(|p| { - let v = world_xf.transform_point3(Vec3::new(p.x, p.y, p.z)); - Vec3A::new(v.x, v.y, v.z) - }) - .collect(); + let mask = prim_ray_mask(prim); + let motion = prim_motion_translate(prim); + + // Non-invertible placements (a zero scale axis) cannot be instanced — + // bake the degenerate transform into world-space triangles as before. + if world_xf.determinant().abs() < 1e-12 { + warn!( + "Mesh at {} has a non-invertible transform — baking instead of instancing", + prim.path() + ); + if motion.is_some() { + warn!( + "Mesh at {}: crust:motion:translate is ignored on baked (non-invertible) geometry", + prim.path() + ); + } + let verts: Vec = points + .iter() + .map(|p| { + let v = world_xf.transform_point3(Vec3::new(p.x, p.y, p.z)); + Vec3A::new(v.x, v.y, v.z) + }) + .collect(); + match triangulate(&counts, &indices, verts.len()) { + Some(tris) => { + world.attach_masked( + Geometry::TriangleMesh { + vertices: verts, + indices: tris, + normals: None, + }, + material, + mask, + ); + } + None => debug!("Mesh at {} produced no triangles", prim.path()), + } + return; + } + + // The mesh's kernel scene is built in the prim's *local* space and + // shared by every prim with identical geometry + material; the + // Instance geometry carries the placement. + let key = MeshKey::new(&points, &counts, &indices, &material); + let inner = match meshes.get(&key) { + Some(shared) => { + debug!("Mesh at {} shares geometry with an earlier prim", prim.path()); + shared.clone() + } + None => { + let verts: Vec = points.iter().map(|p| Vec3A::new(p.x, p.y, p.z)).collect(); + let Some(tris) = triangulate(&counts, &indices, verts.len()) else { + debug!("Mesh at {} produced no triangles", prim.path()); + return; + }; + let mut b = RtSceneBuilder::new(); + b.attach(Geometry::TriangleMesh { + vertices: verts, + indices: tris, + normals: None, + }); + let scene = Arc::new(b.commit()); + meshes.insert(key, scene.clone()); + scene + } + }; + + let l2w = Affine3A::from_mat4(world_xf); + world.attach_masked( + Geometry::Instance { + scene: inner, + transform: l2w, + transform_end: motion.map(|v| Affine3A::from_translation(v) * l2w), + }, + material, + mask, + ); +} - let mut tris: Vec> = Vec::new(); +/// Fan-triangulates the faces into an index-triple list; `None` if +/// nothing survives. +fn triangulate(counts: &[i32], indices: &[i32], n_verts: usize) -> Option> { + let mut tris: Vec<[u32; 3]> = Vec::new(); let mut offset = 0usize; - for &fc in &counts { + for &fc in counts { let fc = fc as usize; if fc < 3 || offset + fc > indices.len() { offset += fc; continue; } for k in 1..(fc - 1) { - let i0 = indices[offset] as usize; - let i1 = indices[offset + k] as usize; - let i2 = indices[offset + k + 1] as usize; - if i0 >= verts.len() || i1 >= verts.len() || i2 >= verts.len() { + let i0 = indices[offset]; + let i1 = indices[offset + k]; + let i2 = indices[offset + k + 1]; + if i0 < 0 || i1 < 0 || i2 < 0 { continue; } - tris.push(Box::new(Triangle::new( - verts[i0], - verts[i1], - verts[i2], - material.clone(), - ))); + let (i0, i1, i2) = (i0 as u32, i1 as u32, i2 as u32); + if i0 as usize >= n_verts || i1 as usize >= n_verts || i2 as usize >= n_verts { + continue; + } + tris.push([i0, i1, i2]); } offset += fc; } - - if tris.is_empty() { - debug!("Mesh at {} produced no triangles", prim.path()); - return; - } - - world.add(Box::new(Bvh::new(tris))); + if tris.is_empty() { None } else { Some(tris) } } // ----------------------------------------------------------------------- @@ -517,7 +650,7 @@ fn emit_mesh( // ----------------------------------------------------------------------- fn emit_sphere( - world: &mut HittableList, + world: &mut WorldBuilder, prim: &Prim, sphere: &UsdSphere, world_xf: GMat4, @@ -542,7 +675,239 @@ fn emit_sphere( radius, center ); - world.add(Box::new(CrustSphere::new(center, radius, material))); + let mask = prim_ray_mask(prim); + match prim_motion_translate(prim) { + // A moving sphere rides an identity-placed instance whose end + // transform is the shutter translation. + Some(v) => { + let mut b = RtSceneBuilder::new(); + b.attach(Geometry::Sphere { center, radius }); + world.attach_masked( + Geometry::Instance { + scene: Arc::new(b.commit()), + transform: Affine3A::IDENTITY, + transform_end: Some(Affine3A::from_translation(v)), + }, + material, + mask, + ); + } + None => { + world.attach_masked(Geometry::Sphere { center, radius }, material, mask); + } + } +} + +// ----------------------------------------------------------------------- +// Curves +// ----------------------------------------------------------------------- + +/// Cubic basis matrices (row-major, `[t³ t² t 1] · M · [P0 P1 P2 P3]ᵀ`). +const BEZIER_M: [[f32; 4]; 4] = [ + [-1.0, 3.0, -3.0, 1.0], + [3.0, -6.0, 3.0, 0.0], + [-3.0, 3.0, 0.0, 0.0], + [1.0, 0.0, 0.0, 0.0], +]; +const BSPLINE_M: [[f32; 4]; 4] = [ + [-1.0 / 6.0, 3.0 / 6.0, -3.0 / 6.0, 1.0 / 6.0], + [3.0 / 6.0, -6.0 / 6.0, 3.0 / 6.0, 0.0], + [-3.0 / 6.0, 0.0, 3.0 / 6.0, 0.0], + [1.0 / 6.0, 4.0 / 6.0, 1.0 / 6.0, 0.0], +]; +const CATMULL_ROM_M: [[f32; 4]; 4] = [ + [-0.5, 1.5, -1.5, 0.5], + [1.0, -2.5, 2.0, -0.5], + [-0.5, 0.0, 0.5, 0.0], + [0.0, 1.0, 0.0, 0.0], +]; + +fn eval_cubic(m: &[[f32; 4]; 4], cp: &[Vec3A; 4], t: f32) -> Vec3A { + let pow = [t * t * t, t * t, t, 1.0]; + let mut p = Vec3A::ZERO; + for (row, &w) in m.iter().zip(&pow) { + for (c, &coeff) in row.iter().enumerate() { + p += cp[c] * (coeff * w); + } + } + p +} + +/// Segments each cubic span is flattened into before intersection. +const CURVE_FLATTEN_SEGS: usize = 8; + +/// Import a `UsdGeomBasisCurves` batch as round (sphere-swept) curve +/// segments: `linear` curves directly, `cubic` curves (bezier / bspline / +/// catmullRom) flattened to a polyline at `CURVE_FLATTEN_SEGS` samples per +/// span. Widths (diameters, per USD) may be authored per point (`vertex`), +/// per curve, or constant; anything else falls back to the first value. +/// The segments live in local space under an `Instance`, like meshes. +fn emit_curves( + world: &mut WorldBuilder, + prim: &Prim, + curves: &UsdBasisCurves, + world_xf: GMat4, + material: Arc, +) { + let points: Option> = curves + .points_attr() + .get::() + .ok() + .flatten() + .and_then(|v| match v { + sdf::Value::Vec3fVec(v) => Some(v), + _ => None, + }); + let counts: Option> = curves + .curve_vertex_counts_attr() + .get::() + .ok() + .flatten() + .and_then(|v| match v { + sdf::Value::IntVec(v) => Some(v), + _ => None, + }); + let (points, counts) = match (points, counts) { + (Some(p), Some(c)) => (p, c), + _ => { + debug!( + "BasisCurves at {} missing points / curveVertexCounts — skipped", + prim.path() + ); + return; + } + }; + let pts: Vec = points.iter().map(|p| Vec3A::new(p.x, p.y, p.z)).collect(); + + let widths: Vec = curves + .widths_attr() + .get::() + .ok() + .flatten() + .and_then(|v| match v { + sdf::Value::FloatVec(v) => Some(v), + _ => None, + }) + .unwrap_or_else(|| vec![1.0]); + + // USD defaults: type = cubic, basis = bezier. + let ty = custom_token(prim, "type").unwrap_or_else(|| "cubic".to_string()); + let basis_name = custom_token(prim, "basis").unwrap_or_else(|| "bezier".to_string()); + let (basis, vstep) = match basis_name.as_str() { + "bezier" => (&BEZIER_M, 3usize), + "bspline" => (&BSPLINE_M, 1), + "catmullRom" => (&CATMULL_ROM_M, 1), + other => { + warn!( + "BasisCurves at {}: basis \"{}\" is not supported (bezier | bspline | catmullRom) — skipped", + prim.path(), + other + ); + return; + } + }; + + // Width (diameter) of control point `global_idx` under the authored + // interpolation, resolved structurally from the array length. + let n_points = pts.len(); + let n_curves = counts.len(); + let width_of = |global_idx: usize, curve_idx: usize| -> f32 { + if widths.len() == n_points { + widths[global_idx] // vertex + } else if widths.len() == n_curves { + widths[curve_idx] // uniform (per curve) + } else { + widths[0] // constant / fallback + } + }; + + let mut segments: Vec = Vec::new(); + let mut offset = 0usize; + for (curve_idx, &cnt) in counts.iter().enumerate() { + let cnt = cnt as usize; + if offset + cnt > pts.len() { + warn!( + "BasisCurves at {}: curveVertexCounts overruns points — remaining curves skipped", + prim.path() + ); + break; + } + let cp = &pts[offset..offset + cnt]; + let radius = + |k: usize| 0.5 * width_of(offset + k, curve_idx).max(1e-6); + + if ty == "linear" { + for k in 0..cnt.saturating_sub(1) { + segments.push(CurveSegment { + p0: cp[k], + p1: cp[k + 1], + r0: radius(k), + r1: radius(k + 1), + }); + } + } else { + // Cubic: flatten each span. Span k uses control points + // [k·vstep .. k·vstep+3]; widths interpolate linearly over the + // curve parameter between the span's end control points. + if cnt < 4 { + offset += cnt; + continue; + } + let n_spans = (cnt - 4) / vstep + 1; + for s in 0..n_spans { + let base = s * vstep; + let ctrl = [cp[base], cp[base + 1], cp[base + 2], cp[base + 3]]; + let (w0, w1) = (radius(base), radius(base + 3)); + let mut prev_p = eval_cubic(basis, &ctrl, 0.0); + let mut prev_r = w0; + for i in 1..=CURVE_FLATTEN_SEGS { + let t = i as f32 / CURVE_FLATTEN_SEGS as f32; + let p = eval_cubic(basis, &ctrl, t); + let r = w0 + (w1 - w0) * t; + segments.push(CurveSegment { + p0: prev_p, + p1: p, + r0: prev_r, + r1: r, + }); + prev_p = p; + prev_r = r; + } + } + } + offset += cnt; + } + + if segments.is_empty() { + debug!("BasisCurves at {} produced no segments", prim.path()); + return; + } + info!( + "Imported BasisCurves at {} ({} {} curves, {} segments)", + prim.path(), + counts.len(), + ty, + segments.len() + ); + + if world_xf.determinant().abs() < 1e-12 { + warn!( + "BasisCurves at {} has a non-invertible transform — skipped", + prim.path() + ); + return; + } + let mut b = RtSceneBuilder::new(); + b.attach(Geometry::RoundCurves { segments }); + world.attach_masked( + Geometry::Instance { + scene: Arc::new(b.commit()), + transform: Affine3A::from_mat4(world_xf), + transform_end: None, + }, + material, + prim_ray_mask(prim), + ); } // ----------------------------------------------------------------------- @@ -635,7 +1000,7 @@ fn lux_emission(light: &impl UsdLight) -> Vec3A { } fn emit_sphere_light( - world: &mut HittableList, + world: &mut WorldBuilder, lights: &mut LightList, light: &SphereLight, world_xf: GMat4, @@ -645,17 +1010,24 @@ fn emit_sphere_light( let pos_v = world_xf.transform_point3(Vec3::ZERO); let position = Vec3A::new(pos_v.x, pos_v.y, pos_v.z); - // One Emissive material shared between the visible sphere geometry and - // the AreaLight — the integrator identifies the light by this material - // when a bounce ray hits it. + // The visible sphere geometry and the AreaLight share one surface; + // the integrator attributes a bounce hit to the light by the + // geometry id the attach returns. let material = Arc::new(Emissive::new(effective)); - world.add(Box::new(CrustSphere::new(position, radius, material.clone()))); + let geom_id = world.attach( + Geometry::Sphere { + center: position, + radius, + }, + material.clone(), + ); lights.add(Arc::new(AreaLight::new( Box::new(SphereShape { center: position, radius, }), material, + geom_id, ))); debug!( "SphereLight: pos={:?} radius={} effective_color={:?}", @@ -664,7 +1036,7 @@ fn emit_sphere_light( } fn emit_rect_light( - world: &mut HittableList, + world: &mut WorldBuilder, lights: &mut LightList, light: &RectLight, world_xf: GMat4, @@ -684,9 +1056,9 @@ fn emit_rect_light( let edge_v = Vec3A::new(ev.x, ev.y, ev.z); let normal = Vec3A::new(nz.x, nz.y, nz.z); - // One Emissive material shared between the visible geometry (two - // triangles spanning the rectangle) and the AreaLight, so the - // integrator can attribute bounce hits to this light by identity. + // The visible geometry (one mesh: two triangles spanning the + // rectangle) and the AreaLight share one surface; bounce hits are + // attributed to the light by the geometry id. let material = Arc::new(Emissive::new(effective)); let (c00, c10, c11, c01) = ( origin, @@ -694,11 +1066,18 @@ fn emit_rect_light( origin + edge_u + edge_v, origin + edge_v, ); - world.add(Box::new(Triangle::new(c00, c10, c11, material.clone()))); - world.add(Box::new(Triangle::new(c00, c11, c01, material.clone()))); + let geom_id = world.attach( + Geometry::TriangleMesh { + vertices: vec![c00, c10, c11, c01], + indices: vec![[0, 1, 2], [0, 2, 3]], + normals: None, + }, + material.clone(), + ); lights.add(Arc::new(AreaLight::new( Box::new(RectShape::new(origin, edge_u, edge_v, normal)), material, + geom_id, ))); debug!( "RectLight: origin={:?} edge_u={:?} edge_v={:?} effective_color={:?}", @@ -745,16 +1124,47 @@ fn warn_unsupported_light(stage: &Stage, prim: &Prim) { // Materials // ----------------------------------------------------------------------- -fn resolve_material(stage: &Stage, prim: &Prim) -> Arc { +/// Memoizes resolved materials by binding path (and shares one default), +/// so prims bound to the same USD material get pointer-identical Arcs — +/// which is what lets `MeshKey` recognize shared mesh geometry. +#[derive(Default)] +struct MaterialCache { + by_path: HashMap>, + default: Option>, +} + +impl MaterialCache { + fn default_material(&mut self) -> Arc { + self.default.get_or_insert_with(default_material).clone() + } +} + +fn resolve_material( + stage: &Stage, + prim: &Prim, + cache: &mut MaterialCache, +) -> Arc { let mat_path = MaterialBindingAPI::get(stage, prim.path().clone()) .ok() .flatten() .and_then(|b| b.direct_binding("").ok().flatten()); let Some(mat_path) = mat_path else { - return default_material(); + return cache.default_material(); }; + if let Some(hit) = cache.by_path.get(mat_path.as_str()) { + return hit.clone(); + } + let resolved = resolve_material_uncached(stage, &mat_path); + cache + .by_path + .insert(mat_path.as_str().to_string(), resolved.clone()); + resolved +} + +fn resolve_material_uncached(stage: &Stage, mat_path: &sdf::Path) -> Arc { + let mat = match UsdMaterial::get(stage, mat_path.clone()) { Ok(Some(m)) => m, _ => { diff --git a/crates/crust-core/src/tracer.rs b/crates/crust-core/src/tracer.rs index 5370b2a..4f71d19 100644 --- a/crates/crust-core/src/tracer.rs +++ b/crates/crust-core/src/tracer.rs @@ -1,12 +1,12 @@ use crate::buffer::Buffer; -use crate::bvh::Bvh; use crate::guiding::{GuidingConfig, GuidingField, SampleData, luminance}; -use crate::hittable::{HitRecord, Hittable}; +use crate::hittable::HitRecord; use crate::material::{Material, ScatterSample}; use crate::medium::sample_henyey_greenstein; use crate::ray::Ray; +use crate::rt_world::{World, WorldHit}; use crate::volume::{PhaseMix, VolumeEvent, Volumes}; -use crate::{LightList, PathSampler, camera::Camera, hittable_list::HittableList}; +use crate::{LightList, PathSampler, camera::Camera}; use glam::Vec3A; use rayon::prelude::*; use std::sync::atomic::{AtomicU64, Ordering}; @@ -17,6 +17,7 @@ use tracing::{info, warn}; // current vertex domain. Distinct keys give independent 4D sub-patterns. const K_CAMERA: i32 = 0; // off root: jitter (0,1) + lens uv (2,3) const K_PATH: i32 = 1; // off root: the bounce subtree +const K_TIME: i32 = 2; // off root: shutter time for motion blur const K_NEE: i32 = 0; // off vertex: light pick (0) + area uv (1,2) const K_NEE_SHADOW: i32 = 1; // off vertex: shadow-ray volume transmittance const K_BSDF: i32 = 2; // off vertex: material scatter block @@ -131,9 +132,9 @@ struct PassStats { pub struct Renderer { pub camera: Camera, - /// Top-level BVH over every scene object (meshes carry their own - /// nested BVH over triangles, built at import). - pub world: Bvh, + /// The committed world: the `crust-rt` kernel scene plus the material + /// table hits resolve through (by `geom_id`). + pub world: World, pub lights: LightList, pub settings: RenderSettings, /// Participating-media regions, sampled outside the BVH (see @@ -145,13 +146,11 @@ pub struct Renderer { impl Renderer { pub fn new( camera: Camera, - world: HittableList, + world: World, lights: LightList, settings: RenderSettings, ) -> Self { - let object_count = world.count(); - let world = Bvh::new(world.into_objects()); - info!("built top-level BVH over {} scene objects", object_count); + info!("world holds {} geometries", world.count()); Renderer { camera, world, @@ -221,7 +220,7 @@ impl Renderer { /// removes here, and the final pass renders unguided instead (the /// training passes still blend in — they are unbiased either way). fn render_guided(&self, tiled: bool, progress: Option) -> Buffer { - let bounds = match self.world.bounding_box() { + let bounds = match self.world.bounds() { Some(b) => b, None => { warn!("path guiding enabled but the scene has no bounding box; rendering unguided"); @@ -474,7 +473,8 @@ impl Renderer { let cam = root.new_domain(K_CAMERA).draw_sample_f32::<4>(); let u = ((i as f32) + cam[0]) / (self.settings.width - 1) as f32; let v = ((j as f32) + cam[1]) / (self.settings.height - 1) as f32; - let r = self.camera.get_ray(u, v, [cam[2], cam[3]]); + let time = root.new_domain(K_TIME).draw_sample_f32::<1>()[0]; + let r = self.camera.get_ray(u, v, [cam[2], cam[3]], time); let color = trace_path( &r, &self.world, @@ -595,7 +595,7 @@ impl RenderSettings { pub fn ray_color( r: &Ray, - world: &dyn Hittable, + world: &World, lights: &LightList, volumes: &Volumes, depth: i32, @@ -758,7 +758,7 @@ enum PrevVertex<'a> { fn bounce_emission_weight( prev: &PrevVertex, lights: &LightList, - hit: &crate::hittable::Hit, + hit: &WorldHit, strategy: SamplingStrategy, ) -> f32 { let (from, bounce_pdf) = match prev { @@ -770,7 +770,7 @@ fn bounce_emission_weight( } PrevVertex::Phase { pos, pdf } => (*pos, *pdf), }; - match lights.find_by_material(hit.mat) { + match lights.find_by_geom(hit.geom_id) { Some(light) => { let light_pdf = (light.pdf(from, hit.rec.p) / lights.count() as f32).max(1e-6); strategy.bounce_weight(bounce_pdf, light_pdf) @@ -785,13 +785,15 @@ fn bounce_emission_weight( /// exact for homogeneous ones). MIS weights are unaffected — transmittance /// is part of the integrand on both strategies, not of either pdf. fn shadow_transmittance( - world: &dyn Hittable, + world: &World, volumes: &Volumes, shadow_ray: &Ray, distance: f32, vertex: PathSampler, ) -> Vec3A { - if world.hit(shadow_ray, 0.001, distance - 0.001).is_some() { + // Dedicated occlusion query: any hit in range means full shadow, so the + // early-exit traversal beats searching for the closest hit. + if world.occluded(shadow_ray, 0.001, distance - 0.001) { return Vec3A::ZERO; } if volumes.is_empty() { @@ -810,11 +812,12 @@ fn volume_nee( p: Vec3A, wi: Vec3A, phase: &PhaseMix, - world: &dyn Hittable, + world: &World, volumes: &Volumes, lights: &LightList, strategy: SamplingStrategy, vertex: PathSampler, + time: f32, ) -> Vec3A { if !strategy.samples_lights() { return Vec3A::ZERO; @@ -828,7 +831,9 @@ fn volume_nee( let light_dir = light_point - p; let light_distance = light_dir.length(); let light_dir_unit = light_dir / light_distance.max(1e-6); - let shadow_ray = Ray::new(p, light_dir_unit); + let shadow_ray = Ray::new(p, light_dir_unit) + .with_time(time) + .with_mask(crate::ray::MASK_SHADOW); let tr = shadow_transmittance(world, volumes, &shadow_ray, light_distance, vertex); if tr == Vec3A::ZERO { return Vec3A::ZERO; @@ -849,7 +854,7 @@ fn volume_nee( /// path and therefore cannot be computed forward. fn trace_path( r: &Ray, - world: &dyn Hittable, + world: &World, lights: &LightList, volumes: &Volumes, depth: i32, @@ -885,7 +890,7 @@ fn trace_path( // the ray itself) but never the background — reproduce both, // attenuating through any media the final segment crosses. if let Some(p) = &prev { - if let Some(hit) = world.hit(&ray, 0.001, f32::INFINITY) { + if let Some(hit) = world.intersect(&ray, 0.001, f32::INFINITY) { let cos_o = ray.direction().normalize().dot(hit.rec.normal).abs(); let mut emitted = hit.mat.emitted_directional(cos_o); if emitted.length_squared() > 0.0 { @@ -905,7 +910,7 @@ fn trace_path( break; } - let hit_opt = world.hit(&ray, 0.001, f32::INFINITY); + let hit_opt = world.intersect(&ray, 0.001, f32::INFINITY); let t_surf = hit_opt.as_ref().map_or(f32::INFINITY, |h| h.rec.t); // Free-flight candidate in the carried homogeneous medium @@ -950,7 +955,8 @@ fn trace_path( let ps = v.new_domain(K_PHASE).draw_sample_f32::<4>(); let dir = phase.sample(wi, ps[0], [ps[1], ps[2]]); let phase_pdf = phase.pdf(wi.dot(dir)).max(1e-6); - let nee = volume_nee(p, wi, &phase, world, volumes, lights, strategy, v); + let nee = + volume_nee(p, wi, &phase, world, volumes, lights, strategy, v, ray.time()); // The walk weight goes into `atten` (it multiplies NEE and // everything beyond); the continuation factor is ONE @@ -992,7 +998,9 @@ fn trace_path( ray = match ray.medium() { Some(m) => Ray::new_in_medium(p, dir, m.clone()), None => Ray::new(p, dir), - }; + } + .with_time(ray.time()) + .with_mask(crate::ray::MASK_INDIRECT); remaining -= 1; continue; } @@ -1055,7 +1063,9 @@ fn trace_path( break; } records.push(vrec); - ray = Ray::new_in_medium(pos, dir, medium); + ray = Ray::new_in_medium(pos, dir, medium) + .with_time(ray.time()) + .with_mask(crate::ray::MASK_INDIRECT); remaining -= 1; prev = None; continue; @@ -1135,7 +1145,9 @@ fn trace_path( let light_distance = light_dir.length(); let light_dir_unit = light_dir.normalize(); - let shadow_ray = Ray::new(rec.p, light_dir_unit); + let shadow_ray = Ray::new(rec.p, light_dir_unit) + .with_time(ray.time()) + .with_mask(crate::ray::MASK_SHADOW); let shadow_tr = shadow_transmittance(world, volumes, &shadow_ray, light_distance, v); @@ -1236,7 +1248,12 @@ fn trace_path( delta: sample.delta, })); records.push(vrec); - ray = sample.ray; + // Materials build the scattered ray without path context; + // stamp the path's shutter time and the indirect category. + ray = sample + .ray + .with_time(ray.time()) + .with_mask(crate::ray::MASK_INDIRECT); remaining -= 1; continue; } diff --git a/crates/crust-core/src/volume.rs b/crates/crust-core/src/volume.rs index 5b9f094..5a8d610 100644 --- a/crates/crust-core/src/volume.rs +++ b/crates/crust-core/src/volume.rs @@ -388,7 +388,7 @@ impl Volumes { if region.majorant_sigma_t <= 0.0 { continue; } - if !region.world_aabb.hit(ray, t_eps, t_max) { + if !region.world_aabb.hit(ray.rt(), t_eps, t_max) { continue; } if let Some((t0, t1)) = region.intersect(ray) { diff --git a/crates/crust-core/src/world.rs b/crates/crust-core/src/world.rs index 6c34063..b63916f 100644 --- a/crates/crust-core/src/world.rs +++ b/crates/crust-core/src/world.rs @@ -1,42 +1,44 @@ use crate::RenderSettings; -use crate::Sphere; use crate::camera::Camera; -use crate::hittable_list::HittableList; use crate::light::{AreaLight, LightList, SphereShape}; use crate::material::{Emissive, Material, OpenPBR}; +use crate::rt_world::{World, WorldBuilder}; +use crust_rt::Geometry; use glam::Vec3A; use std::sync::Arc; use utils::{random_range3, random3}; +/// Attaches an analytic sphere with its material. +fn add_sphere(world: &mut WorldBuilder, center: Vec3A, radius: f32, material: Arc) { + world.attach(Geometry::Sphere { center, radius }, material); +} + /// Adds a sphere light to the scene: emissive sphere geometry in `world` /// plus an `AreaLight` over the same surface in `lights`, tied together by -/// sharing one `Emissive` material (Cornell-box semantics). +/// the geometry id (Cornell-box semantics — one surface, both roles). fn add_sphere_light( - world: &mut HittableList, + world: &mut WorldBuilder, lights: &mut LightList, color: Vec3A, center: Vec3A, radius: f32, ) { let material = Arc::new(Emissive::new(color)); - world.add(Box::new(Sphere::new(center, radius, material.clone()))); + let geom_id = world.attach(Geometry::Sphere { center, radius }, material.clone()); lights.add(Arc::new(AreaLight::new( Box::new(SphereShape { center, radius }), material, + geom_id, ))); } #[allow(dead_code)] -pub fn random_scene() -> (HittableList, LightList) { - let mut world = HittableList::new(); +pub fn random_scene() -> (World, LightList) { + let mut world = WorldBuilder::new(); let mut lights = LightList::new(); let ground_material = Arc::new(OpenPBR::diffuse(Vec3A::new(0.5, 0.5, 0.5))); - world.add(Box::new(Sphere::new( - Vec3A::new(0.0, -1000.0, 0.0), - 1000.0, - ground_material, - ))); + add_sphere(&mut world, Vec3A::new(0.0, -1000.0, 0.0), 1000.0, ground_material); for a in -11..11 { for b in -11..11 { @@ -68,31 +70,19 @@ pub fn random_scene() -> (HittableList, LightList) { // Glass Arc::new(OpenPBR::glass(1.5)) }; - world.add(Box::new(Sphere::new(center, 0.2, sphere_material))); + add_sphere(&mut world, center, 0.2, sphere_material); } } } let material1 = Arc::new(OpenPBR::glass(1.5)); - world.add(Box::new(Sphere::new( - Vec3A::new(0.0, 1.0, 0.0), - 1.0, - material1, - ))); + add_sphere(&mut world, Vec3A::new(0.0, 1.0, 0.0), 1.0, material1); let material2 = Arc::new(OpenPBR::diffuse(Vec3A::new(0.4, 0.2, 0.1))); - world.add(Box::new(Sphere::new( - Vec3A::new(-4.0, 1.0, 0.0), - 1.0, - material2, - ))); + add_sphere(&mut world, Vec3A::new(-4.0, 1.0, 0.0), 1.0, material2); let material3 = Arc::new(OpenPBR::metal(Vec3A::new(0.7, 0.6, 0.5), 0.0)); - world.add(Box::new(Sphere::new( - Vec3A::new(4.0, 1.0, 0.0), - 1.0, - material3, - ))); + add_sphere(&mut world, Vec3A::new(4.0, 1.0, 0.0), 1.0, material3); add_sphere_light( &mut world, @@ -109,19 +99,15 @@ pub fn random_scene() -> (HittableList, LightList) { 1.0, ); - (world, lights) + (world.commit(), lights) } -pub fn simple_scene() -> (HittableList, LightList) { - let mut world = HittableList::new(); +pub fn simple_scene() -> (World, LightList) { + let mut world = WorldBuilder::new(); let mut lights = LightList::new(); let ground_material = Arc::new(OpenPBR::diffuse(Vec3A::new(0.8, 0.5, 0.5))); - world.add(Box::new(Sphere::new( - Vec3A::new(0.0, -1000.0, 0.0), - 1000.0, - ground_material, - ))); + add_sphere(&mut world, Vec3A::new(0.0, -1000.0, 0.0), 1000.0, ground_material); // Deterministic grid of spheres with preset materials for a in -2..3 { @@ -134,38 +120,22 @@ pub fn simple_scene() -> (HittableList, LightList) { _ => Arc::new(OpenPBR::glossy(Vec3A::new(0.9, 0.9, 0.9), 0.2, 0.5)), }; - world.add(Box::new(Sphere::new(center, 0.2, material))); + add_sphere(&mut world, center, 0.2, material); } } // Center spheres let material1 = Arc::new(OpenPBR::glass(1.5)); - world.add(Box::new(Sphere::new( - Vec3A::new(0.0, 1.0, 0.0), - 1.0, - material1, - ))); + add_sphere(&mut world, Vec3A::new(0.0, 1.0, 0.0), 1.0, material1); let material2 = Arc::new(OpenPBR::diffuse(Vec3A::new(0.4, 0.2, 0.1))); - world.add(Box::new(Sphere::new( - Vec3A::new(-4.0, 1.0, 0.0), - 1.0, - material2, - ))); + add_sphere(&mut world, Vec3A::new(-4.0, 1.0, 0.0), 1.0, material2); let material3 = Arc::new(OpenPBR::metal(Vec3A::new(0.7, 0.6, 0.5), 0.0)); - world.add(Box::new(Sphere::new( - Vec3A::new(4.0, 1.0, 0.0), - 1.0, - material3, - ))); + add_sphere(&mut world, Vec3A::new(4.0, 1.0, 0.0), 1.0, material3); let material4 = Arc::new(OpenPBR::glossy(Vec3A::new(0.5, 0.5, 0.5), 0.2, 0.0)); - world.add(Box::new(Sphere::new( - Vec3A::new(0.0, 1.0, 4.0), - 1.0, - material4, - ))); + add_sphere(&mut world, Vec3A::new(0.0, 1.0, 4.0), 1.0, material4); let material5 = Arc::new(OpenPBR { base_color: Vec3A::new(0.5, 0.5, 0.5), @@ -173,11 +143,7 @@ pub fn simple_scene() -> (HittableList, LightList) { coat_weight: 1.0, ..OpenPBR::default() }); - world.add(Box::new(Sphere::new( - Vec3A::new(0.0, 1.0, -4.0), - 1.0, - material5, - ))); + add_sphere(&mut world, Vec3A::new(0.0, 1.0, -4.0), 1.0, material5); // Lights add_sphere_light( @@ -195,7 +161,7 @@ pub fn simple_scene() -> (HittableList, LightList) { 1.0, ); - (world, lights) + (world.commit(), lights) } pub fn get_settings() -> (Camera, RenderSettings) { diff --git a/crates/crust-core/tests/usd_scene.rs b/crates/crust-core/tests/usd_scene.rs index 2a62a79..c390368 100644 --- a/crates/crust-core/tests/usd_scene.rs +++ b/crates/crust-core/tests/usd_scene.rs @@ -1,6 +1,6 @@ use std::path::PathBuf; -use crust_core::{Hittable, Scene}; +use crust_core::Scene; use openusd::schemas::shade::{Material as UsdMaterial, MaterialBindingAPI}; use openusd::sdf; use openusd::usd::{PrimPredicate, Stage}; @@ -81,7 +81,7 @@ fn cornellbox_transforms_compose_correctly() { .expect("failed to open cornellbox.usda"); let bbox = scene .world - .bounding_box() + .bounds() .expect("cornellbox world must be bounded"); let tol = 0.1; @@ -110,8 +110,8 @@ fn loads_rectlight_usda() { // Ball sphere + floor mesh BVH + two triangles of rect-light geometry. assert_eq!( scene.world.count(), - 4, - "expected 4 hittables (sphere, floor, 2 light triangles), got {}", + 3, + "expected 3 geometries (sphere, floor, rect-light mesh), got {}", scene.world.count() ); // The RectLight must import as a real light, not warn-and-skip. @@ -160,8 +160,8 @@ fn loads_fog_usda() { // must import as a volume region, NOT as geometry. assert_eq!( scene.world.count(), - 4, - "expected 4 hittables (room, ball, 2 light triangles), got {}", + 3, + "expected 3 geometries (room, ball, rect-light mesh), got {}", scene.world.count() ); assert_eq!(scene.lights.count(), 1); @@ -194,8 +194,8 @@ fn loads_smoke_usda() { // as regions, not geometry. assert_eq!( scene.world.count(), - 3, - "expected 3 hittables (room, 2 light triangles), got {}", + 2, + "expected 2 geometries (room, rect-light mesh), got {}", scene.world.count() ); assert_eq!(scene.lights.count(), 1); @@ -285,3 +285,98 @@ fn openpbr_showcase_materials_all_decode() { bound ); } + +#[test] +fn loads_curves_usda() { + let scene = Scene::from_usd(&sample("curves.usda")).expect("failed to open curves.usda"); + + // Tuft instance + Tripod instance + floor instance + 2 light triangles. + assert_eq!( + scene.world.count(), + 4, + "expected 4 geometries (2 curve batches, floor, rect-light mesh), got {}", + scene.world.count() + ); + + // The linear Tripod strand's first segment rises from (1.6, 0, -0.5) + // to (1.6, 0.8, -0.3) (after the prim's translate); a -Z ray at its + // mid-height must hit it, and slightly to the side must miss it. + let on_axis = crust_core::Ray::new( + crust_core::Vec3A::new(1.6, 0.4, 5.0), + -crust_core::Vec3A::Z, + ); + let hit = scene + .world + .intersect(&on_axis, 0.001, f32::INFINITY) + .expect("ray through the strand must hit"); + assert!( + (hit.rec.t - 5.4).abs() < 0.1, + "strand hit at t={} (expected ~5.4)", + hit.rec.t + ); + let wide = crust_core::Ray::new( + crust_core::Vec3A::new(2.6, 0.4, 5.0), + -crust_core::Vec3A::Z, + ); + // The wide ray flies past every strand and over the floor edge... but the + // floor extends to z=-8, so aim slightly upward to clear it entirely. + let wide_up = crust_core::Ray::new( + crust_core::Vec3A::new(2.6, 0.4, 5.0), + (crust_core::Vec3A::new(2.6, 3.0, -8.0) - crust_core::Vec3A::new(2.6, 0.4, 5.0)).normalize(), + ); + assert!(scene.world.intersect(&wide, 0.001, 4.0).is_none()); + assert!(scene.world.intersect(&wide_up, 0.001, f32::INFINITY).is_none()); +} + +#[test] +fn loads_motionblur_usda() { + let scene = + Scene::from_usd(&sample("motionblur.usda")).expect("failed to open motionblur.usda"); + + // Mover sphere, Riser cube, floor, shadow card, 2 light triangles. + assert_eq!( + scene.world.count(), + 5, + "expected 5 geometries, got {}", + scene.world.count() + ); + + // The sphere starts at (-1.5, 0.6, 0) and streaks +1 in x over the + // shutter: a time-0 ray down its start position hits, a time-1 ray at + // the same spot misses, and a time-1 ray at the end position hits. + let at = |x: f32, time: f32| { + crust_core::Ray::new( + crust_core::Vec3A::new(x, 0.6, 6.0), + -crust_core::Vec3A::Z, + ) + .with_time(time) + }; + assert!(scene.world.intersect(&at(-1.5, 0.0), 0.001, 5.9).is_some()); + assert!(scene.world.intersect(&at(-1.5, 1.0), 0.001, 5.9).is_none()); + assert!(scene.world.intersect(&at(-0.5, 1.0), 0.001, 5.9).is_some()); + + // The shadow card (crust:rayMask = 6) is invisible to camera rays but + // opaque to shadow rays. + let down = crust_core::Ray::new( + crust_core::Vec3A::new(0.0, 5.0, 0.0), + -crust_core::Vec3A::Y, + ); + let cam_hit = scene + .world + .intersect(&down.clone().with_mask(crust_core::MASK_CAMERA), 0.001, f32::INFINITY) + .expect("camera ray passes the card and hits the floor"); + assert!( + (cam_hit.rec.t - 5.0).abs() < 1e-3, + "camera ray should reach the floor at t=5, got {}", + cam_hit.rec.t + ); + let shadow_hit = scene + .world + .intersect(&down.clone().with_mask(crust_core::MASK_SHADOW), 0.001, f32::INFINITY) + .expect("shadow ray must be blocked by the card"); + assert!( + (shadow_hit.rec.t - 3.0).abs() < 1e-3, + "shadow ray should stop at the card at t=3, got {}", + shadow_hit.rec.t + ); +} diff --git a/crates/crust-rt/Cargo.toml b/crates/crust-rt/Cargo.toml new file mode 100644 index 0000000..45a837b --- /dev/null +++ b/crates/crust-rt/Cargo.toml @@ -0,0 +1,14 @@ +[package] +name = "crust-rt" +version = "0.1.0" +edition = { workspace = true } +license-file = { workspace = true } +description = "Crust Render's ray tracing kernel: Embree-shaped scene/geometry API over a parallel SBVH build collapsed to BVH4, in safe Rust" + +[lib] +name = "crust_rt" +path = "src/lib.rs" + +[dependencies] +glam.workspace = true +rayon = "1.10.0" diff --git a/crates/crust-rt/src/aabb.rs b/crates/crust-rt/src/aabb.rs new file mode 100644 index 0000000..493a676 --- /dev/null +++ b/crates/crust-rt/src/aabb.rs @@ -0,0 +1,58 @@ +use crate::ray::Ray; +use glam::Vec3A; + +#[derive(Debug, Clone, Copy)] +pub struct AABB { + pub minimum: Vec3A, + pub maximum: Vec3A, +} + +impl AABB { + pub fn new(minimum: Vec3A, maximum: Vec3A) -> Self { + AABB { minimum, maximum } + } + + pub fn surrounding_box(box0: AABB, box1: AABB) -> AABB { + AABB { + minimum: box0.minimum.min(box1.minimum), + maximum: box0.maximum.max(box1.maximum), + } + } + + /// Scalar slab test. The BVH traversal uses its own 4-wide variant; + /// this is for callers working with a single box. + pub fn hit(&self, ray: &Ray, mut t_min: f32, mut t_max: f32) -> bool { + for a in 0..3 { + let inv_d = 1.0 / ray.dir[a]; + let mut t0 = (self.minimum[a] - ray.origin[a]) * inv_d; + let mut t1 = (self.maximum[a] - ray.origin[a]) * inv_d; + + if inv_d < 0.0 { + std::mem::swap(&mut t0, &mut t1); + } + + t_min = t_min.max(t0); + t_max = t_max.min(t1); + + if t_max <= t_min { + return false; + } + } + true + } +} + +/// Bounds of a triangle, padded on degenerate axes so axis-aligned +/// geometry survives the slab test (which rejects zero-thickness boxes). +pub fn triangle_aabb(v0: Vec3A, v1: Vec3A, v2: Vec3A) -> AABB { + let mut min = v0.min(v1).min(v2); + let mut max = v0.max(v1).max(v2); + const PAD: f32 = 1e-4; + for a in 0..3 { + if max[a] - min[a] < PAD { + min[a] -= PAD; + max[a] += PAD; + } + } + AABB::new(min, max) +} diff --git a/crates/crust-rt/src/bvh.rs b/crates/crust-rt/src/bvh.rs new file mode 100644 index 0000000..d5c3fa6 --- /dev/null +++ b/crates/crust-rt/src/bvh.rs @@ -0,0 +1,972 @@ +//! Bounding-volume hierarchy over the internal primitives. +//! +//! The build is reference-based, in the SBVH mold (Stich et al. 2009): +//! every subtree owns a list of `PrimRef`s (bbox + primitive index). Each +//! node first evaluates a binned surface-area-heuristic *object* split on +//! the longest centroid axis; when the object split's children overlap by +//! more than `SBVH_ALPHA` of the root surface area it also evaluates a +//! binned *spatial* split, which chops straddling references in two at the +//! split plane (bounds via `Prim::clipped_aabb` — exact polygon clipping +//! for triangles) and duplicates them into both children. The cheaper +//! split wins. Leaves therefore index into a shared `indices` array (a +//! primitive can appear in several leaves) while the primitives +//! themselves are stored exactly once. +//! +//! The binary tree is then collapsed into 4-wide SoA nodes whose slab +//! tests run on `Vec4` lanes. Large subtrees build in parallel via +//! `rayon::join`; every split decision depends only on the input, so the +//! tree is deterministic — threads only change *when* subtrees are built, +//! never *what*. + +use crate::aabb::AABB; +use crate::prim::{Prim, PrimHit}; +use crate::ray::Ray; +use glam::{Vec3A, Vec4}; + +/// Leaves are forced at this depth so the traversal stack can never +/// overflow; SAH partitions are otherwise free to be arbitrarily uneven. +const MAX_DEPTH: usize = 60; +/// Ranges at or below this size always become a leaf. +const MIN_LEAF: usize = 2; +/// A range no larger than this may stay a leaf when splitting is not +/// worth it by SAH cost; larger ranges are always split. +const MAX_LEAF: usize = 8; +/// Number of candidate split planes tested per axis. +const BINS: usize = 12; +/// Subtrees larger than this build their children on parallel rayon tasks. +const PARALLEL_THRESHOLD: usize = 4096; +/// Spatial splits are only *considered* when the object split's children +/// overlap by more than this fraction of the root surface area (the SBVH +/// α of Stich et al. 2009) — the natural brake on reference duplication. +const SBVH_ALPHA: f32 = 1e-5; +/// And never below this depth: chopping tiny deep subtrees duplicates +/// references for negligible traversal gain. +const SBVH_MAX_DEPTH: usize = 32; + +/// Binary build node — an intermediate: the finished tree is the 4-wide +/// [`WideNode`] array produced by collapsing these. +struct Node { + bbox: AABB, + /// Leaf (`count > 0`): offset of the first entry in `indices`. + /// Internal (`count == 0`): index of the right child — the left child + /// immediately follows the node itself in depth-first order. + first_or_right: u32, + count: u32, +} + +/// Marks an unused lane of a [`WideNode`] (together with `count == 0`). +const EMPTY_LANE: u32 = u32::MAX; + +/// A 4-wide BVH node in SoA layout: lane `k` of each `Vec4` holds child +/// `k`'s slab bounds, so one round of vector min/max tests all four child +/// boxes against the ray at once (Embree's BVH4 idea). +struct WideNode { + bmin_x: Vec4, + bmin_y: Vec4, + bmin_z: Vec4, + bmax_x: Vec4, + bmax_y: Vec4, + bmax_z: Vec4, + /// Leaf lane (`count > 0`): first offset into `indices`. Internal lane + /// (`count == 0`): wide-node index, or [`EMPTY_LANE`]. + child: [u32; 4], + count: [u32; 4], +} + +impl WideNode { + fn empty() -> Self { + WideNode { + // +INF/+INF bounds: the swapped slab test can never pass them + // (an *inverted* box would — min/max swapping un-inverts it). + bmin_x: Vec4::INFINITY, + bmin_y: Vec4::INFINITY, + bmin_z: Vec4::INFINITY, + bmax_x: Vec4::INFINITY, + bmax_y: Vec4::INFINITY, + bmax_z: Vec4::INFINITY, + child: [EMPTY_LANE; 4], + count: [0; 4], + } + } + + fn set_lane_bounds(&mut self, lane: usize, b: &AABB) { + self.bmin_x[lane] = b.minimum.x; + self.bmin_y[lane] = b.minimum.y; + self.bmin_z[lane] = b.minimum.z; + self.bmax_x[lane] = b.maximum.x; + self.bmax_y[lane] = b.maximum.y; + self.bmax_z[lane] = b.maximum.z; + } +} + +pub(crate) struct Bvh { + wide: Vec, + /// Leaf ranges index into this; spatial splits may list a primitive in + /// more than one leaf. + indices: Vec, + /// The primitives, stored once each, in input order. + prims: Vec>, + /// Bounds of the whole tree (the binary root's, kept through collapse). + root_bbox: Option, +} + +/// One build reference: conservative bounds of (a fragment of) primitive +/// `idx`. Spatial splits shrink the bounds and duplicate the reference. +#[derive(Clone, Copy)] +struct PrimRef { + bbox: AABB, + idx: u32, +} + +impl PrimRef { + fn centroid(&self) -> Vec3A { + 0.5 * (self.bbox.minimum + self.bbox.maximum) + } +} + +/// A built subtree with node indices and leaf offsets local to itself; +/// `merge` splices children under a parent, offsetting as it goes. +struct Subtree { + nodes: Vec, + indices: Vec, +} + +impl Bvh { + pub(crate) fn new(prims: Vec>) -> Self { + let refs: Vec = prims + .iter() + .enumerate() + .map(|(i, p)| PrimRef { + bbox: p.bbox(), + idx: i as u32, + }) + .collect(); + + let (wide, indices, root_bbox) = if refs.is_empty() { + (Vec::new(), Vec::new(), None) + } else { + let root_bbox = union_all(&refs); + let subtree = build_subtree(&prims, refs, 0, surface_area(&root_bbox)); + (collapse(&subtree.nodes), subtree.indices, Some(root_bbox)) + }; + + Bvh { + wide, + indices, + prims, + root_bbox, + } + } + + pub(crate) fn prim_count(&self) -> usize { + self.prims.len() + } + + pub(crate) fn bounds(&self) -> Option { + self.root_bbox + } + + /// Closest hit in `(t_min, t_max)`. + pub(crate) fn hit(&self, ray: &Ray, t_min: f32, t_max: f32) -> Option { + if self.wide.is_empty() { + return None; + } + let mut closest = t_max; + let mut best: Option = None; + + let o = ray.origin; + let inv = Vec3A::new(safe_inv(ray.dir.x), safe_inv(ray.dir.y), safe_inv(ray.dir.z)); + + let mut stack = [0u32; WIDE_STACK]; + stack[0] = 0; + let mut sp = 1usize; + + while sp > 0 { + sp -= 1; + let node = &self.wide[stack[sp] as usize]; + let (tnear, tfar) = slab4(node, o, inv, t_min, closest); + + // Hit lanes, insertion-sorted near-to-far (≤ 4 entries). + let mut order = [(0f32, 0usize); 4]; + let mut n_hit = 0; + for l in 0..4 { + if node.count[l] == 0 && node.child[l] == EMPTY_LANE { + continue; + } + if tnear[l] <= tfar[l] { + let mut i = n_hit; + while i > 0 && order[i - 1].0 > tnear[l] { + order[i] = order[i - 1]; + i -= 1; + } + order[i] = (tnear[l], l); + n_hit += 1; + } + } + + // Leaf lanes intersect immediately (near first, shrinking + // `closest`); internal lanes are pushed far-to-near so the + // nearest pops first. + for i in 0..n_hit { + let l = order[i].1; + if node.count[l] > 0 { + let first = node.child[l] as usize; + for &pi in &self.indices[first..first + node.count[l] as usize] { + if let Some(hit) = self.prims[pi as usize].hit(ray, t_min, closest) { + closest = hit.t; + best = Some(hit); + } + } + } + } + for i in (0..n_hit).rev() { + let l = order[i].1; + if node.count[l] == 0 { + stack[sp] = node.child[l]; + sp += 1; + } + } + } + + best + } + + /// Early-exit occlusion traversal: no ordering, returns on the first + /// confirmed hit anywhere in `(t_min, t_max)`. + pub(crate) fn hit_any(&self, ray: &Ray, t_min: f32, t_max: f32) -> bool { + if self.wide.is_empty() { + return false; + } + let o = ray.origin; + let inv = Vec3A::new(safe_inv(ray.dir.x), safe_inv(ray.dir.y), safe_inv(ray.dir.z)); + + let mut stack = [0u32; WIDE_STACK]; + stack[0] = 0; + let mut sp = 1usize; + + while sp > 0 { + sp -= 1; + let node = &self.wide[stack[sp] as usize]; + let (tnear, tfar) = slab4(node, o, inv, t_min, t_max); + for l in 0..4 { + if node.count[l] == 0 && node.child[l] == EMPTY_LANE { + continue; + } + if tnear[l] > tfar[l] { + continue; + } + if node.count[l] > 0 { + let first = node.child[l] as usize; + for &pi in &self.indices[first..first + node.count[l] as usize] { + if self.prims[pi as usize].hit_any(ray, t_min, t_max) { + return true; + } + } + } else { + stack[sp] = node.child[l]; + sp += 1; + } + } + } + false + } +} + +/// Finite reciprocal of a direction component: zero (and denormal-tiny) +/// components become a huge same-signed value instead of ±∞, so the slab +/// arithmetic can never produce the 0·∞ = NaN that poisons vector min/max. +#[inline] +fn safe_inv(x: f32) -> f32 { + if x.abs() < 1e-20 { + 1e20f32.copysign(x) + } else { + 1.0 / x + } +} + +/// Traversal stack: wide depth ≤ binary `MAX_DEPTH`, and each visited node +/// leaves at most 3 deferred lanes behind, so 3·MAX_DEPTH + 4 slots hold +/// the worst case. +const WIDE_STACK: usize = 3 * MAX_DEPTH + 4; + +/// The 4-lane slab test: entry/exit distances for all four child boxes of +/// `node` at once. A lane hits iff `tnear[l] <= tfar[l]`. +#[inline] +fn slab4(node: &WideNode, o: Vec3A, inv: Vec3A, t_min: f32, t_max: f32) -> (Vec4, Vec4) { + let t0x = (node.bmin_x - Vec4::splat(o.x)) * Vec4::splat(inv.x); + let t1x = (node.bmax_x - Vec4::splat(o.x)) * Vec4::splat(inv.x); + let t0y = (node.bmin_y - Vec4::splat(o.y)) * Vec4::splat(inv.y); + let t1y = (node.bmax_y - Vec4::splat(o.y)) * Vec4::splat(inv.y); + let t0z = (node.bmin_z - Vec4::splat(o.z)) * Vec4::splat(inv.z); + let t1z = (node.bmax_z - Vec4::splat(o.z)) * Vec4::splat(inv.z); + let tnear = t0x + .min(t1x) + .max(t0y.min(t1y)) + .max(t0z.min(t1z)) + .max(Vec4::splat(t_min)); + let tfar = t0x + .max(t1x) + .min(t0y.max(t1y)) + .min(t0z.max(t1z)) + .min(Vec4::splat(t_max)); + (tnear, tfar) +} + +fn surface_area(b: &AABB) -> f32 { + let d = b.maximum - b.minimum; + 2.0 * (d.x * d.y + d.y * d.z + d.z * d.x) +} + +fn union_all(refs: &[PrimRef]) -> AABB { + refs.iter() + .skip(1) + .fold(refs[0].bbox, |acc, r| AABB::surrounding_box(acc, r.bbox)) +} + +/// Component-wise intersection; `None` when the boxes do not overlap. +fn intersect_aabb(a: &AABB, b: &AABB) -> Option { + let lo = a.minimum.max(b.minimum); + let hi = a.maximum.min(b.maximum); + if lo.cmple(hi).all() { + Some(AABB::new(lo, hi)) + } else { + None + } +} + +fn leaf(bbox: AABB, refs: &[PrimRef]) -> Subtree { + Subtree { + nodes: vec![Node { + bbox, + first_or_right: 0, + count: refs.len() as u32, + }], + indices: refs.iter().map(|r| r.idx).collect(), + } +} + +/// Splices `left`/`right` under a fresh internal node, rewriting the +/// children's local node indices and leaf offsets into the merged frame. +fn merge(bbox: AABB, left: Subtree, right: Subtree) -> Subtree { + let mut nodes = Vec::with_capacity(1 + left.nodes.len() + right.nodes.len()); + let right_node_offset = 1 + left.nodes.len() as u32; + nodes.push(Node { + bbox, + first_or_right: right_node_offset, + count: 0, + }); + for mut n in left.nodes { + if n.count == 0 { + n.first_or_right += 1; + } + nodes.push(n); + } + let leaf_offset = left.indices.len() as u32; + for mut n in right.nodes { + if n.count == 0 { + n.first_or_right += right_node_offset; + } else { + n.first_or_right += leaf_offset; + } + nodes.push(n); + } + let mut indices = left.indices; + indices.extend(right.indices); + Subtree { nodes, indices } +} + +/// The winning object split: partition predicate parameters plus the SAH +/// cost and the (unclipped) child bounds the overlap test needs. +struct ObjSplit { + axis: usize, + /// Centroid-to-bin mapping: `bin = ((c - cmin) * scale) as usize`. + cmin: f32, + scale: f32, + split_bin: usize, + cost: f32, + left_bbox: AABB, + right_bbox: AABB, +} + +/// The winning spatial split: chop plane on `axis` at `pos`. +struct SpatSplit { + axis: usize, + pos: f32, + cost: f32, +} + +fn best_object_split(refs: &[PrimRef]) -> Option { + // Split along the longest axis of the centroid bounds. + let mut cmin = refs[0].centroid(); + let mut cmax = cmin; + for r in &refs[1..] { + let c = r.centroid(); + cmin = cmin.min(c); + cmax = cmax.max(c); + } + let extent = cmax - cmin; + let axis = if extent.x >= extent.y && extent.x >= extent.z { + 0 + } else if extent.y >= extent.z { + 1 + } else { + 2 + }; + if extent[axis] <= 1e-6 { + // All centroids coincide — nothing to partition on. + return None; + } + let scale = BINS as f32 / extent[axis]; + let bin_of = |r: &PrimRef| (((r.centroid()[axis] - cmin[axis]) * scale) as usize).min(BINS - 1); + + // Binned SAH: histogram the centroids, then score every split plane + // between adjacent bins by `area · count` on each side. + let mut bin_counts = [0usize; BINS]; + let mut bin_bounds: [Option; BINS] = [None; BINS]; + for r in refs { + let b = bin_of(r); + bin_counts[b] += 1; + bin_bounds[b] = Some(match bin_bounds[b] { + Some(existing) => AABB::surrounding_box(existing, r.bbox), + None => r.bbox, + }); + } + + let mut best: Option = None; + for split in 0..BINS - 1 { + let mut lb: Option = None; + let mut lc = 0usize; + for b in 0..=split { + lc += bin_counts[b]; + lb = match (lb, bin_bounds[b]) { + (Some(x), Some(y)) => Some(AABB::surrounding_box(x, y)), + (x, y) => x.or(y), + }; + } + let mut rb: Option = None; + let mut rc = 0usize; + for b in split + 1..BINS { + rc += bin_counts[b]; + rb = match (rb, bin_bounds[b]) { + (Some(x), Some(y)) => Some(AABB::surrounding_box(x, y)), + (x, y) => x.or(y), + }; + } + if lc == 0 || rc == 0 { + continue; + } + let (lb, rb) = (lb.expect("lc > 0"), rb.expect("rc > 0")); + let cost = surface_area(&lb) * lc as f32 + surface_area(&rb) * rc as f32; + if best.as_ref().is_none_or(|b| cost < b.cost) { + best = Some(ObjSplit { + axis, + cmin: cmin[axis], + scale, + split_bin: split, + cost, + left_bbox: lb, + right_bbox: rb, + }); + } + } + best +} + +/// Binned spatial split on the node's longest axis: references straddling +/// a bin contribute their *clipped* bounds to it (entry/exit counting), so +/// the candidate children reflect what duplication would actually produce. +fn best_spatial_split(prims: &[Box], refs: &[PrimRef], bbox: &AABB) -> Option { + let extent = bbox.maximum - bbox.minimum; + let axis = if extent.x >= extent.y && extent.x >= extent.z { + 0 + } else if extent.y >= extent.z { + 1 + } else { + 2 + }; + if extent[axis] <= 1e-6 { + return None; + } + let lo = bbox.minimum[axis]; + let width = extent[axis] / BINS as f32; + let bin_of = |x: f32| (((x - lo) / width) as usize).clamp(0, BINS - 1); + + let mut entry = [0usize; BINS]; + let mut exit = [0usize; BINS]; + let mut bounds: [Option; BINS] = [None; BINS]; + let mut add = |b: usize, aabb: AABB| { + bounds[b] = Some(match bounds[b] { + Some(existing) => AABB::surrounding_box(existing, aabb), + None => aabb, + }); + }; + + for r in refs { + let b0 = bin_of(r.bbox.minimum[axis]); + let b1 = bin_of(r.bbox.maximum[axis]); + entry[b0] += 1; + exit[b1] += 1; + if b0 == b1 { + add(b0, r.bbox); + continue; + } + for b in b0..=b1 { + let (bin_lo, bin_hi) = (lo + b as f32 * width, lo + (b + 1) as f32 * width); + if let Some(c) = prims[r.idx as usize] + .clipped_aabb(axis, bin_lo, bin_hi) + .and_then(|c| intersect_aabb(&c, &r.bbox)) + { + add(b, c); + } + } + } + + let mut best: Option = None; + for split in 0..BINS - 1 { + let mut lb: Option = None; + let mut lc = 0usize; + for b in 0..=split { + lc += entry[b]; + lb = match (lb, bounds[b]) { + (Some(x), Some(y)) => Some(AABB::surrounding_box(x, y)), + (x, y) => x.or(y), + }; + } + let mut rb: Option = None; + let mut rc = 0usize; + for b in split + 1..BINS { + rc += exit[b]; + rb = match (rb, bounds[b]) { + (Some(x), Some(y)) => Some(AABB::surrounding_box(x, y)), + (x, y) => x.or(y), + }; + } + if lc == 0 || rc == 0 { + continue; + } + let cost = surface_area(&lb.expect("lc > 0")) * lc as f32 + + surface_area(&rb.expect("rc > 0")) * rc as f32; + if best.as_ref().is_none_or(|b| cost < b.cost) { + best = Some(SpatSplit { + axis, + pos: lo + (split + 1) as f32 * width, + cost, + }); + } + } + best +} + +/// Builds the subtree for `refs`, appending nodes in depth-first order +/// (locally indexed — `merge` rebases children). `root_area` normalizes +/// the SBVH overlap test. +fn build_subtree( + prims: &[Box], + mut refs: Vec, + depth: usize, + root_area: f32, +) -> Subtree { + let bbox = union_all(&refs); + let count = refs.len(); + if count <= MIN_LEAF || depth >= MAX_DEPTH { + return leaf(bbox, &refs); + } + + let object = best_object_split(&refs); + + // Consider a spatial split only when the object split's children + // overlap enough for chopping to pay for the duplicated references. + let spatial = match &object { + Some(o) if depth < SBVH_MAX_DEPTH => { + let overlap = + intersect_aabb(&o.left_bbox, &o.right_bbox).map_or(0.0, |b| surface_area(&b)); + if overlap / root_area > SBVH_ALPHA { + best_spatial_split(prims, &refs, &bbox).filter(|s| s.cost < o.cost) + } else { + None + } + } + _ => None, + }; + + let (left_refs, right_refs) = if let Some(s) = spatial { + // Chop: straddling references are clipped into both children. + let mut left = Vec::with_capacity(count); + let mut right = Vec::with_capacity(count); + for r in refs { + if r.bbox.maximum[s.axis] <= s.pos { + left.push(r); + } else if r.bbox.minimum[s.axis] >= s.pos { + right.push(r); + } else { + let prim = &prims[r.idx as usize]; + if let Some(c) = prim + .clipped_aabb(s.axis, f32::NEG_INFINITY, s.pos) + .and_then(|c| intersect_aabb(&c, &r.bbox)) + { + left.push(PrimRef { bbox: c, idx: r.idx }); + } + if let Some(c) = prim + .clipped_aabb(s.axis, s.pos, f32::INFINITY) + .and_then(|c| intersect_aabb(&c, &r.bbox)) + { + right.push(PrimRef { bbox: c, idx: r.idx }); + } + } + } + // Degenerate chop (numeric edge): fall back to a leaf-or-object + // path rather than recursing on an empty side. + if left.is_empty() || right.is_empty() { + let mut refs: Vec = left; + refs.extend(right); + return object_partition_or_leaf(prims, refs, bbox, object, depth, root_area); + } + (left, right) + } else { + match object { + Some(o) => { + // Leaf when splitting costs more than intersecting through. + if count <= MAX_LEAF && o.cost >= surface_area(&bbox) * count as f32 { + return leaf(bbox, &refs); + } + partition_by_bin(&mut refs, &o) + } + // Every centroid coincides: median split by input order. + None => { + let mid = count / 2; + let right = refs.split_off(mid); + (refs, right) + } + } + }; + + let parallel = left_refs.len().max(right_refs.len()) > PARALLEL_THRESHOLD; + let (l, r) = if parallel { + rayon::join( + || build_subtree(prims, left_refs, depth + 1, root_area), + || build_subtree(prims, right_refs, depth + 1, root_area), + ) + } else { + ( + build_subtree(prims, left_refs, depth + 1, root_area), + build_subtree(prims, right_refs, depth + 1, root_area), + ) + }; + merge(bbox, l, r) +} + +/// The non-spatial tail of `build_subtree`, reused by the degenerate-chop +/// fallback: object-partition when possible, else leaf. +fn object_partition_or_leaf( + prims: &[Box], + mut refs: Vec, + bbox: AABB, + object: Option, + depth: usize, + root_area: f32, +) -> Subtree { + let count = refs.len(); + match object { + Some(o) if count > MIN_LEAF => { + let (l, r) = partition_by_bin(&mut refs, &o); + if l.is_empty() || r.is_empty() { + let mut all = l; + all.extend(r); + return leaf(bbox, &all); + } + let left = build_subtree(prims, l, depth + 1, root_area); + let right = build_subtree(prims, r, depth + 1, root_area); + merge(bbox, left, right) + } + _ => leaf(bbox, &refs), + } +} + +/// Order-preserving partition of `refs` by the object split's centroid +/// bin — deterministic for a given input order. +fn partition_by_bin(refs: &mut Vec, o: &ObjSplit) -> (Vec, Vec) { + let mut left = Vec::with_capacity(refs.len()); + let mut right = Vec::with_capacity(refs.len()); + for r in refs.drain(..) { + let bin = (((r.centroid()[o.axis] - o.cmin) * o.scale) as usize).min(BINS - 1); + if bin <= o.split_bin { + left.push(r); + } else { + right.push(r); + } + } + (left, right) +} + +/// Collapses the binary tree into 4-wide nodes: each wide node adopts its +/// binary node's two children, then repeatedly replaces the largest-area +/// internal child with that child's own two children until four lanes are +/// filled (or only leaves remain). Purely input-driven, so determinism is +/// preserved. +fn collapse(binary: &[Node]) -> Vec { + let mut out = Vec::with_capacity(binary.len() / 2 + 1); + if binary.is_empty() { + return out; + } + if binary[0].count > 0 { + // Single-leaf tree. + let mut w = WideNode::empty(); + w.set_lane_bounds(0, &binary[0].bbox); + w.child[0] = binary[0].first_or_right; + w.count[0] = binary[0].count; + out.push(w); + return out; + } + collapse_node(binary, 0, &mut out); + out +} + +fn collapse_node(binary: &[Node], b_idx: u32, out: &mut Vec) -> u32 { + let slot = out.len(); + out.push(WideNode::empty()); + + let mut kids = [0u32; 4]; + kids[0] = b_idx + 1; + kids[1] = binary[b_idx as usize].first_or_right; + let mut n_kids = 2; + while n_kids < 4 { + // Expand the internal kid with the largest surface area; ties + // resolve to the first (deterministic). + let mut best: Option<(usize, f32)> = None; + for (i, &k) in kids.iter().enumerate().take(n_kids) { + if binary[k as usize].count == 0 { + let a = surface_area(&binary[k as usize].bbox); + if best.is_none_or(|(_, ba)| a > ba) { + best = Some((i, a)); + } + } + } + let Some((i, _)) = best else { break }; + let k = kids[i]; + kids[i] = k + 1; + kids[n_kids] = binary[k as usize].first_or_right; + n_kids += 1; + } + + for lane in 0..n_kids { + let k = kids[lane] as usize; + let bounds = binary[k].bbox; + out[slot].set_lane_bounds(lane, &bounds); + if binary[k].count > 0 { + out[slot].child[lane] = binary[k].first_or_right; + out[slot].count[lane] = binary[k].count; + } else { + let ci = collapse_node(binary, kids[lane], out); + out[slot].child[lane] = ci; + } + } + slot as u32 +} + +#[cfg(test)] +mod tests { + use super::*; + use crate::prim::{SpherePrim, TrianglePrim}; + use crate::ray::MASK_ALL; + + fn sphere_grid(n: i32) -> Vec> { + let mut out: Vec> = Vec::new(); + for x in 0..n { + for y in 0..n { + for z in 0..n { + out.push(Box::new(SpherePrim { + center: Vec3A::new(x as f32, y as f32, z as f32) * 3.0, + radius: 0.5, + geom_id: (x * n * n + y * n + z) as u32, + mask: MASK_ALL, + })); + } + } + } + out + } + + /// Long thin diagonal triangles — the geometry spatial splits exist + /// for. Built so object splits alone leave heavily overlapping + /// children. + fn diagonal_shards(n: i32) -> Vec> { + let mut out: Vec> = Vec::new(); + for i in 0..n { + let o = i as f32 * 0.35; + out.push(Box::new(TrianglePrim { + v0: Vec3A::new(o, o, o), + v1: Vec3A::new(o + 10.0, o + 10.0, o + 10.2), + v2: Vec3A::new(o + 10.0, o + 10.3, o + 10.0), + normals: None, + geom_id: 0, + prim_id: i as u32, + mask: MASK_ALL, + })); + } + out + } + + fn linear_scan(prims: &[Box], ray: &Ray, t_min: f32, t_max: f32) -> Option { + let mut closest = t_max; + let mut best = None; + for p in prims { + if let Some(h) = p.hit(ray, t_min, closest) { + closest = h.t; + best = Some(h); + } + } + best + } + + fn assert_matches_linear(objects: impl Fn() -> Vec>) { + let bvh = Bvh::new(objects()); + let reference = objects(); + + let origins = [ + Vec3A::new(-5.0, 4.5, 4.5), + Vec3A::new(20.0, 3.0, 3.0), + Vec3A::new(4.5, -5.0, 4.5), + Vec3A::new(0.0, 0.0, -10.0), + Vec3A::new(5.0, 5.2, -3.0), + ]; + let dirs = [ + Vec3A::new(1.0, 0.0, 0.0), + Vec3A::new(-1.0, 0.05, 0.02).normalize(), + Vec3A::new(0.0, 1.0, 0.0), + Vec3A::new(0.3, 0.3, 1.0).normalize(), + Vec3A::new(0.0, 0.0, -1.0), + Vec3A::new(0.577, 0.577, 0.577), + ]; + for o in origins { + for d in dirs { + let ray = Ray::new(o, d); + let a = bvh.hit(&ray, 0.001, f32::INFINITY); + let b = linear_scan(&reference, &ray, 0.001, f32::INFINITY); + match (a, b) { + (Some(x), Some(y)) => { + assert!((x.t - y.t).abs() < 1e-4, "t mismatch for {o:?} {d:?}"); + assert_eq!(x.geom_id, y.geom_id, "id mismatch for {o:?} {d:?}"); + } + (None, None) => {} + (x, y) => panic!( + "hit disagreement for {o:?} {d:?}: bvh={} linear={}", + x.is_some(), + y.is_some() + ), + } + assert_eq!( + bvh.hit_any(&ray, 0.001, f32::INFINITY), + b.is_some(), + "occlusion disagreement for {o:?} {d:?}" + ); + } + } + } + + /// The BVH must find exactly the hits a linear scan finds. + #[test] + fn matches_linear_scan() { + assert_matches_linear(|| sphere_grid(4)); + } + + /// Same, on geometry that triggers spatial splits (verified below). + #[test] + fn spatial_splits_match_linear_scan() { + assert_matches_linear(|| diagonal_shards(64)); + } + + /// Diagonal shards must actually produce duplicated references — + /// otherwise the spatial-split path is dead code. + #[test] + fn spatial_splits_duplicate_references() { + let bvh = Bvh::new(diagonal_shards(64)); + assert!( + bvh.indices.len() > bvh.prims.len(), + "no reference duplication: {} indices for {} prims", + bvh.indices.len(), + bvh.prims.len() + ); + } + + /// `hit_any` must agree with `hit(..).is_some()` for every ray and range. + #[test] + fn hit_any_matches_hit() { + let bvh = Bvh::new(sphere_grid(4)); + let origins = [ + Vec3A::new(-5.0, 4.5, 4.5), + Vec3A::new(20.0, 3.0, 3.0), + Vec3A::new(4.5, 4.5, 4.5), + ]; + let dirs = [ + Vec3A::new(1.0, 0.0, 0.0), + Vec3A::new(-1.0, 0.05, 0.02).normalize(), + Vec3A::new(0.3, 0.3, 1.0).normalize(), + ]; + for o in origins { + for d in dirs { + let ray = Ray::new(o, d); + for t_max in [0.5, 3.0, f32::INFINITY] { + assert_eq!( + bvh.hit_any(&ray, 0.001, t_max), + bvh.hit(&ray, 0.001, t_max).is_some(), + "occlusion disagreement for {o:?} {d:?} t_max={t_max}" + ); + } + } + } + } + + /// Parallel subtree builds must not change the tree: the same input + /// always produces byte-identical topology. + #[test] + fn build_is_deterministic() { + let a = Bvh::new(sphere_grid(6)); + let b = Bvh::new(sphere_grid(6)); + assert_eq!(a.wide.len(), b.wide.len()); + assert_eq!(a.indices, b.indices); + for (x, y) in a.wide.iter().zip(&b.wide) { + assert_eq!(x.child, y.child); + assert_eq!(x.count, y.count); + assert_eq!(x.bmin_x, y.bmin_x); + assert_eq!(x.bmax_z, y.bmax_z); + } + } + + /// The collapse must actually widen: a big tree ends up with far + /// fewer wide nodes than a binary tree would need. + #[test] + fn collapse_widens_the_tree() { + let bvh = Bvh::new(sphere_grid(6)); // 216 prims + let n_leaf_slots: usize = bvh + .wide + .iter() + .flat_map(|w| w.count.iter()) + .filter(|&&c| c > 0) + .count(); + assert!(n_leaf_slots > 0); + // A binary tree over L leaves has L-1 internal nodes; BVH4 should + // need roughly a third of that. + assert!( + bvh.wide.len() * 2 < n_leaf_slots.max(2) * 2 - 1, + "{} wide nodes for {} leaves", + bvh.wide.len(), + n_leaf_slots + ); + } + + #[test] + fn empty_bvh_misses() { + let bvh = Bvh::new(Vec::new()); + let ray = Ray::new(Vec3A::ZERO, Vec3A::X); + assert!(bvh.hit(&ray, 0.001, f32::INFINITY).is_none()); + assert!(bvh.bounds().is_none()); + } + + #[test] + fn bounds_cover_all_prims() { + let bvh = Bvh::new(sphere_grid(3)); + let bbox = bvh.bounds().expect("grid is fully bounded"); + assert!(bbox.minimum.cmple(Vec3A::splat(-0.5)).all()); + assert!(bbox.maximum.cmpge(Vec3A::splat(6.5)).all()); + } +} diff --git a/crates/crust-rt/src/curve.rs b/crates/crust-rt/src/curve.rs new file mode 100644 index 0000000..fa5bdba --- /dev/null +++ b/crates/crust-rt/src/curve.rs @@ -0,0 +1,181 @@ +//! Rounded-cone intersection — the round (sphere-swept) curve segment. + +use crate::ray::Ray; +use glam::Vec3A; + +/// Nearest boundary hit of the rounded cone (sphere-swept segment) in +/// `[t_min, t_max]`, as `(t, outward_normal)`. The solid is the union of +/// the two cap spheres and the tangent cone body between them; taking the +/// minimum valid `t` over all three surfaces yields the hull entry point +/// for rays starting outside (rays starting *inside* may report an +/// interior sphere surface — irrelevant for opaque hair). +/// +/// Following Quilez's rounded-cone intersector, evaluated with a +/// normalized direction and rescaled back to the caller's parameter. +pub(crate) fn rounded_cone_intersect( + ray: &Ray, + p0: Vec3A, + p1: Vec3A, + r0: f32, + r1: f32, + t_min: f32, + t_max: f32, +) -> Option<(f32, Vec3A)> { + let r0 = r0.max(1e-6); + let r1 = r1.max(1e-6); + let len = ray.dir.length(); + if len < 1e-20 { + return None; + } + let rd = ray.dir / len; + let o = ray.origin; + // Candidate ts below are in normalized units; the range test happens + // there too, then the winner converts back. + let (n_min, n_max) = (t_min * len, t_max * len); + + let ba = p1 - p0; + let oa = o - p0; + let ob = o - p1; + let rr = r0 - r1; + let m0 = ba.dot(ba); + let m1 = ba.dot(oa); + let m2 = ba.dot(rd); + let m3 = rd.dot(oa); + let m5 = oa.dot(oa); + let m6 = ob.dot(rd); + let m7 = ob.dot(ob); + + let mut best: Option<(f32, Vec3A)> = None; + let mut consider = |t: f32, n: Vec3A| { + if t >= n_min && t <= n_max && best.map_or(true, |(bt, _)| t < bt) { + best = Some((t, n)); + } + }; + + // Cone body between the tangency circles (0 < y < d2), when neither + // sphere swallows the other. + let d2 = m0 - rr * rr; + if d2 > 0.0 { + let k2 = d2 - m2 * m2; + let k1 = d2 * m3 - m1 * m2 + m2 * rr * r0; + let k0 = d2 * m5 - m1 * m1 + m1 * rr * r0 * 2.0 - m0 * r0 * r0; + let h = k1 * k1 - k0 * k2; + if h >= 0.0 && k2.abs() > 1e-12 { + let sq = h.sqrt(); + for t in [(-k1 - sq) / k2, (-k1 + sq) / k2] { + let y = m1 - r0 * rr + t * m2; + if y > 0.0 && y < d2 { + consider(t, (d2 * (oa + t * rd) - ba * y).normalize()); + } + } + } + } + + // Cap spheres. + let h0 = m3 * m3 - m5 + r0 * r0; + if h0 >= 0.0 { + let sq = h0.sqrt(); + for t in [-m3 - sq, -m3 + sq] { + consider(t, (oa + t * rd) / r0); + } + } + let h1 = m6 * m6 - m7 + r1 * r1; + if h1 >= 0.0 { + let sq = h1.sqrt(); + for t in [-m6 - sq, -m6 + sq] { + consider(t, (ob + t * rd) / r1); + } + } + + best.map(|(t, n)| (t / len, n)) +} + +#[cfg(test)] +mod tests { + use super::*; + + fn hit(o: Vec3A, d: Vec3A, p0: Vec3A, p1: Vec3A, r0: f32, r1: f32) -> Option<(f32, Vec3A)> { + rounded_cone_intersect(&Ray::new(o, d), p0, p1, r0, r1, 0.001, f32::INFINITY) + } + + #[test] + fn capsule_axial_hit_at_cap() { + // Capsule from (0,0,0) to (4,0,0), radius 0.5; axial ray from -X + // hits the p0 cap sphere at x = -0.5. + let (t, n) = hit( + Vec3A::new(-3.0, 0.0, 0.0), + Vec3A::X, + Vec3A::ZERO, + Vec3A::new(4.0, 0.0, 0.0), + 0.5, + 0.5, + ) + .expect("axial hit"); + assert!((t - 2.5).abs() < 1e-4, "t = {t}"); + assert!(n.abs_diff_eq(-Vec3A::X, 1e-4)); + } + + #[test] + fn capsule_perpendicular_hit_at_radius() { + let (t, n) = hit( + Vec3A::new(2.0, 3.0, 0.0), + -Vec3A::Y, + Vec3A::ZERO, + Vec3A::new(4.0, 0.0, 0.0), + 0.5, + 0.5, + ) + .expect("side hit"); + assert!((t - 2.5).abs() < 1e-3, "t = {t}"); + assert!(n.abs_diff_eq(Vec3A::Y, 1e-3)); + } + + #[test] + fn cone_radius_shrinks_along_axis() { + // r goes 0.5 -> 0.1 over x in [0,4]; a perpendicular ray near the + // thin end must graze closer to the axis than one near the thick end. + let (p0, p1) = (Vec3A::ZERO, Vec3A::new(4.0, 0.0, 0.0)); + let t_thick = hit(Vec3A::new(0.5, 3.0, 0.0), -Vec3A::Y, p0, p1, 0.5, 0.1) + .expect("thick hit") + .0; + let t_thin = hit(Vec3A::new(3.5, 3.0, 0.0), -Vec3A::Y, p0, p1, 0.5, 0.1) + .expect("thin hit") + .0; + let surf_thick = 3.0 - t_thick; // height of surface above axis + let surf_thin = 3.0 - t_thin; + assert!( + surf_thick > surf_thin + 0.2, + "cone does not taper: {surf_thick} vs {surf_thin}" + ); + assert!(surf_thick <= 0.5 + 1e-3 && surf_thin >= 0.1 - 1e-3); + } + + #[test] + fn respects_t_range_and_misses() { + let (p0, p1) = (Vec3A::ZERO, Vec3A::new(4.0, 0.0, 0.0)); + let o = Vec3A::new(2.0, 3.0, 0.0); + let r = Ray::new(o, -Vec3A::Y); + assert!(rounded_cone_intersect(&r, p0, p1, 0.5, 0.5, 0.001, 2.0).is_none()); + // Ray passing wide of the capsule. + assert!(hit(Vec3A::new(2.0, 3.0, 2.0), -Vec3A::Y, p0, p1, 0.5, 0.5).is_none()); + // Unnormalized direction: same surface point, halved t. + let (t, _) = hit(o, -Vec3A::Y * 2.0, p0, p1, 0.5, 0.5).expect("hit"); + assert!((t - 1.25).abs() < 1e-3); + } + + #[test] + fn degenerate_swallowed_sphere_still_hits() { + // Segment shorter than the radius difference: the big sphere + // swallows the small one; hits must come from the big sphere. + let (t, _) = hit( + Vec3A::new(0.0, 0.0, -5.0), + Vec3A::Z, + Vec3A::ZERO, + Vec3A::new(0.1, 0.0, 0.0), + 1.0, + 0.05, + ) + .expect("hit"); + assert!((t - 4.0).abs() < 1e-3); + } +} diff --git a/crates/crust-rt/src/lib.rs b/crates/crust-rt/src/lib.rs new file mode 100644 index 0000000..480d025 --- /dev/null +++ b/crates/crust-rt/src/lib.rs @@ -0,0 +1,47 @@ +//! Crust Render's ray tracing kernel — the intersection layer only, +//! factored out of the renderer the way `openqmc-rs` factored out +//! sampling. The API is deliberately **Embree-shaped** (geometry objects +//! attached to a scene, a commit step, `intersect`/`occluded` queries, +//! ID-based hits) so a renderer written against it could swap in Embree +//! bindings behind the same seam — while this native implementation stays +//! 100% safe Rust. +//! +//! ``` +//! use crust_rt::{Geometry, Ray, SceneBuilder}; +//! use glam::Vec3A; +//! +//! let mut builder = SceneBuilder::new(); +//! let ball = builder.attach(Geometry::Sphere { center: Vec3A::ZERO, radius: 1.0 }); +//! let scene = builder.commit(); +//! +//! let hit = scene +//! .intersect(&Ray::new(Vec3A::new(0.0, 0.0, -5.0), Vec3A::Z), 0.001, f32::INFINITY) +//! .expect("hit"); +//! assert_eq!(hit.geom_id, ball); +//! assert!((hit.t - 4.0).abs() < 1e-4); +//! assert!(!scene.occluded(&Ray::new(Vec3A::new(0.0, 0.0, -5.0), Vec3A::Z), 0.001, 3.9)); +//! ``` +//! +//! What Embree calls `rtcIntersect1` / `rtcOccluded1` are +//! [`Scene::intersect`] / [`Scene::occluded`]; geometry types map to the +//! [`Geometry`] enum (triangle meshes with optional per-vertex shading +//! normals, analytic spheres, round curve segments, and single-level +//! instances with optional transform motion blur); per-geometry visibility +//! masks test against the ray's category bit exactly like Embree's. Hits +//! carry `geom_id`/`prim_id` — the application owns the mapping from IDs +//! to materials or anything else; this crate never sees shading data. + +mod aabb; +mod bvh; +mod curve; +mod prim; +mod ray; +mod scene; +mod triangle; + +pub use aabb::AABB; +pub use ray::{MASK_ALL, MASK_CAMERA, MASK_INDIRECT, MASK_SHADOW, Ray}; +pub use scene::{CurveSegment, Geometry, RayHit, Scene, SceneBuilder}; + +/// The `geom_id`/`prim_id` value that never names a real geometry. +pub const INVALID_ID: u32 = u32::MAX; diff --git a/crates/crust-rt/src/prim.rs b/crates/crust-rt/src/prim.rs new file mode 100644 index 0000000..888a526 --- /dev/null +++ b/crates/crust-rt/src/prim.rs @@ -0,0 +1,308 @@ +//! Internal primitives the BVH is built over. Each carries the IDs and +//! the visibility mask of the geometry it came from; a hit is plain +//! `Copy` data (no lifetimes, no shading state). + +use crate::aabb::{AABB, triangle_aabb}; +use crate::curve::rounded_cone_intersect; +use crate::ray::Ray; +use crate::scene::Scene; +use crate::triangle::{clip_triangle_aabb, triangle_intersect}; +use glam::{Affine3A, Mat3A, Vec3A}; +use std::sync::Arc; + +/// An intersection as the primitives report it: `outward` is the +/// *geometric outward* normal (not yet oriented against the ray) — the +/// public API flips it and derives `front_face` at the query edge, so +/// instance transforms can map it without bookkeeping. +#[derive(Clone, Copy)] +pub(crate) struct PrimHit { + pub t: f32, + pub outward: Vec3A, + pub u: f32, + pub v: f32, + pub geom_id: u32, + pub prim_id: u32, +} + +pub(crate) trait Prim: Send + Sync { + fn hit(&self, ray: &Ray, t_min: f32, t_max: f32) -> Option; + + /// Boolean occlusion variant; overridden where cheaper than `hit`. + fn hit_any(&self, ray: &Ray, t_min: f32, t_max: f32) -> bool { + self.hit(ray, t_min, t_max).is_some() + } + + fn bbox(&self) -> AABB; + + /// Conservative bounds of the part inside the axis slab, for the + /// BVH's spatial splits. Default: bbox clipped to the slab. + fn clipped_aabb(&self, axis: usize, min: f32, max: f32) -> Option { + let b = self.bbox(); + if b.minimum[axis] > max || b.maximum[axis] < min { + return None; + } + let mut c = b; + c.minimum[axis] = c.minimum[axis].max(min); + c.maximum[axis] = c.maximum[axis].min(max); + Some(c) + } +} + +#[inline] +fn masked_out(ray: &Ray, mask: u32) -> bool { + ray.mask & mask == 0 +} + +// --------------------------------------------------------------------- +// Triangle (optionally with per-vertex shading normals) +// --------------------------------------------------------------------- + +pub(crate) struct TrianglePrim { + pub v0: Vec3A, + pub v1: Vec3A, + pub v2: Vec3A, + /// Per-vertex shading normals; the reported normal interpolates them + /// by the hit barycentrics when present. + pub normals: Option<[Vec3A; 3]>, + pub geom_id: u32, + pub prim_id: u32, + pub mask: u32, +} + +impl Prim for TrianglePrim { + fn hit(&self, ray: &Ray, t_min: f32, t_max: f32) -> Option { + if masked_out(ray, self.mask) { + return None; + } + let (t, u, v) = triangle_intersect(ray, self.v0, self.v1, self.v2, t_min, t_max)?; + let outward = match &self.normals { + Some([n0, n1, n2]) => (*n0 * (1.0 - u - v) + *n1 * u + *n2 * v).normalize(), + None => { + let n = (self.v1 - self.v0).cross(self.v2 - self.v0); + if n == Vec3A::ZERO { + return None; // degenerate sliver + } + n.normalize() + } + }; + Some(PrimHit { + t, + outward, + u, + v, + geom_id: self.geom_id, + prim_id: self.prim_id, + }) + } + + fn hit_any(&self, ray: &Ray, t_min: f32, t_max: f32) -> bool { + !masked_out(ray, self.mask) + && triangle_intersect(ray, self.v0, self.v1, self.v2, t_min, t_max).is_some() + } + + fn bbox(&self) -> AABB { + triangle_aabb(self.v0, self.v1, self.v2) + } + + fn clipped_aabb(&self, axis: usize, min: f32, max: f32) -> Option { + clip_triangle_aabb(self.v0, self.v1, self.v2, axis, min, max) + } +} + +// --------------------------------------------------------------------- +// Sphere +// --------------------------------------------------------------------- + +pub(crate) struct SpherePrim { + pub center: Vec3A, + pub radius: f32, + pub geom_id: u32, + pub mask: u32, +} + +impl Prim for SpherePrim { + fn hit(&self, ray: &Ray, t_min: f32, t_max: f32) -> Option { + if masked_out(ray, self.mask) { + return None; + } + let oc = ray.origin - self.center; + let a = ray.dir.length_squared(); + let half_b = oc.dot(ray.dir); + let c = oc.length_squared() - self.radius * self.radius; + let discriminant = half_b * half_b - a * c; + if discriminant < 0.0 { + return None; + } + let sqrt_d = discriminant.sqrt(); + let mut root = (-half_b - sqrt_d) / a; + if root <= t_min || root >= t_max { + root = (-half_b + sqrt_d) / a; + if root <= t_min || root >= t_max { + return None; + } + } + Some(PrimHit { + t: root, + outward: (ray.at(root) - self.center) / self.radius, + u: 0.0, + v: 0.0, + geom_id: self.geom_id, + prim_id: 0, + }) + } + + fn bbox(&self) -> AABB { + AABB::new( + self.center - Vec3A::splat(self.radius), + self.center + Vec3A::splat(self.radius), + ) + } +} + +// --------------------------------------------------------------------- +// Round curve segment (sphere-swept cone) +// --------------------------------------------------------------------- + +pub(crate) struct CurvePrim { + pub p0: Vec3A, + pub p1: Vec3A, + pub r0: f32, + pub r1: f32, + pub geom_id: u32, + pub prim_id: u32, + pub mask: u32, +} + +impl Prim for CurvePrim { + fn hit(&self, ray: &Ray, t_min: f32, t_max: f32) -> Option { + if masked_out(ray, self.mask) { + return None; + } + let (t, outward) = + rounded_cone_intersect(ray, self.p0, self.p1, self.r0, self.r1, t_min, t_max)?; + Some(PrimHit { + t, + outward, + u: 0.0, + v: 0.0, + geom_id: self.geom_id, + prim_id: self.prim_id, + }) + } + + fn bbox(&self) -> AABB { + let a = AABB::new(self.p0 - Vec3A::splat(self.r0), self.p0 + Vec3A::splat(self.r0)); + let b = AABB::new(self.p1 - Vec3A::splat(self.r1), self.p1 + Vec3A::splat(self.r1)); + AABB::surrounding_box(a, b) + } +} + +// --------------------------------------------------------------------- +// Instance: a committed scene placed by a transform, with optional +// transform motion blur (linear matrix interpolation at the ray's time). +// --------------------------------------------------------------------- + +pub(crate) struct InstancePrim { + pub scene: Arc, + /// Local-to-world at shutter time 0. + pub l2w: Affine3A, + /// Cached world-to-local and inverse-transpose at time 0. + pub w2l: Affine3A, + pub normal_mat: Mat3A, + /// Local-to-world at shutter time 1, when the instance moves. + pub l2w_end: Option, + pub bounds: AABB, + pub geom_id: u32, + pub mask: u32, +} + +/// Element-wise linear interpolation of two affine transforms. Every +/// interpolated point stays inside the convex hull of its endpoint +/// positions, so a union-of-endpoints bound is conservative. +pub(crate) fn lerp_affine(a: &Affine3A, b: &Affine3A, t: f32) -> Affine3A { + Affine3A { + matrix3: Mat3A::from_cols( + a.matrix3.x_axis.lerp(b.matrix3.x_axis, t), + a.matrix3.y_axis.lerp(b.matrix3.y_axis, t), + a.matrix3.z_axis.lerp(b.matrix3.z_axis, t), + ), + translation: a.translation.lerp(b.translation, t), + } +} + +/// World-space box of `local` under `m`: transformed corner bounds, +/// padded on degenerate axes like `triangle_aabb`. +pub(crate) fn transformed_aabb(local: &AABB, m: &Affine3A) -> AABB { + let mut min = Vec3A::splat(f32::INFINITY); + let mut max = Vec3A::splat(f32::NEG_INFINITY); + for i in 0..8 { + let corner = Vec3A::new( + if i & 1 == 0 { local.minimum.x } else { local.maximum.x }, + if i & 2 == 0 { local.minimum.y } else { local.maximum.y }, + if i & 4 == 0 { local.minimum.z } else { local.maximum.z }, + ); + let p = m.transform_point3a(corner); + min = min.min(p); + max = max.max(p); + } + const PAD: f32 = 1e-4; + for a in 0..3 { + if max[a] - min[a] < PAD { + min[a] -= PAD; + max[a] += PAD; + } + } + AABB::new(min, max) +} + +impl InstancePrim { + /// World-to-local and normal transform at the ray's shutter time. + fn transforms_at(&self, time: f32) -> (Affine3A, Mat3A) { + match &self.l2w_end { + Some(end) if time > 0.0 => { + let w2l = lerp_affine(&self.l2w, end, time).inverse(); + (w2l, w2l.matrix3.transpose()) + } + _ => (self.w2l, self.normal_mat), + } + } + + fn to_local(&self, ray: &Ray, w2l: &Affine3A) -> Ray { + Ray { + origin: w2l.transform_point3a(ray.origin), + // Unnormalized on purpose: local t == world t. + dir: w2l.transform_vector3a(ray.dir), + time: ray.time, + mask: ray.mask, + } + } +} + +impl Prim for InstancePrim { + fn hit(&self, ray: &Ray, t_min: f32, t_max: f32) -> Option { + if masked_out(ray, self.mask) { + return None; + } + let (w2l, normal_mat) = self.transforms_at(ray.time); + let local = self.to_local(ray, &w2l); + let mut hit = self.scene.intersect_outward(&local, t_min, t_max)?; + hit.outward = (normal_mat * hit.outward).normalize(); + // The hit is attributed to the *instance's* geometry id: the + // application maps materials per top-level geometry. The inner + // primitive index is kept. + hit.geom_id = self.geom_id; + Some(hit) + } + + fn hit_any(&self, ray: &Ray, t_min: f32, t_max: f32) -> bool { + if masked_out(ray, self.mask) { + return false; + } + let (w2l, _) = self.transforms_at(ray.time); + self.scene.occluded(&self.to_local(ray, &w2l), t_min, t_max) + } + + fn bbox(&self) -> AABB { + self.bounds + } +} diff --git a/crates/crust-rt/src/ray.rs b/crates/crust-rt/src/ray.rs new file mode 100644 index 0000000..f72a8e7 --- /dev/null +++ b/crates/crust-rt/src/ray.rs @@ -0,0 +1,55 @@ +use glam::Vec3A; + +/// Ray visibility categories, Embree-mask style: a ray carries the bit of +/// the category it belongs to, geometry carries a mask of the categories +/// it is visible to, and an intersection only counts when +/// `geometry_mask & ray_mask != 0`. Geometry and rays both default to +/// [`MASK_ALL`], so masking is zero-cost until something opts in. +pub const MASK_CAMERA: u32 = 1 << 0; +pub const MASK_SHADOW: u32 = 1 << 1; +pub const MASK_INDIRECT: u32 = 1 << 2; +pub const MASK_ALL: u32 = u32::MAX; + +/// A ray: origin, (unnormalized) direction, shutter `time` in `[0, 1)` for +/// motion blur, and the visibility category `mask`. Plain `Copy` data — +/// renderer-side state like the participating medium a path is inside +/// belongs to the caller, not the kernel. +#[derive(Clone, Copy, Debug)] +pub struct Ray { + pub origin: Vec3A, + pub dir: Vec3A, + pub time: f32, + pub mask: u32, +} + +impl Default for Ray { + fn default() -> Self { + Ray::new(Vec3A::ZERO, Vec3A::ZERO) + } +} + +impl Ray { + /// A ray at shutter time 0, visible to all geometry. + pub fn new(origin: Vec3A, dir: Vec3A) -> Ray { + Ray { + origin, + dir, + time: 0.0, + mask: MASK_ALL, + } + } + + pub fn with_time(mut self, time: f32) -> Ray { + self.time = time; + self + } + + pub fn with_mask(mut self, mask: u32) -> Ray { + self.mask = mask; + self + } + + pub fn at(&self, t: f32) -> Vec3A { + self.origin + t * self.dir + } +} diff --git a/crates/crust-rt/src/scene.rs b/crates/crust-rt/src/scene.rs new file mode 100644 index 0000000..6996239 --- /dev/null +++ b/crates/crust-rt/src/scene.rs @@ -0,0 +1,494 @@ +//! The Embree-shaped public API: attach [`Geometry`] objects to a +//! [`SceneBuilder`], `commit()` into an immutable [`Scene`], query with +//! `intersect` / `occluded`. + +use crate::aabb::AABB; +use crate::bvh::Bvh; +use crate::prim::{ + CurvePrim, InstancePrim, Prim, PrimHit, SpherePrim, TrianglePrim, transformed_aabb, +}; +use crate::ray::{MASK_ALL, Ray}; +use glam::{Affine3A, Vec3A}; +use std::sync::Arc; + +/// One round (sphere-swept) curve segment: a cone frustum tangent to the +/// spheres `(p0, r0)` and `(p1, r1)`, with spherical caps. +#[derive(Clone, Copy, Debug)] +pub struct CurveSegment { + pub p0: Vec3A, + pub p1: Vec3A, + pub r0: f32, + pub r1: f32, +} + +/// A geometry to attach to a scene. The variants mirror Embree's geometry +/// types (the subset crust needs): triangle meshes, analytic spheres, +/// round curves, and single-level instances of another committed scene. +pub enum Geometry { + TriangleMesh { + vertices: Vec, + indices: Vec<[u32; 3]>, + /// Optional per-vertex shading normals; hits interpolate them by + /// the barycentrics (`SmoothTriangle` semantics). + normals: Option>, + }, + Sphere { + center: Vec3A, + radius: f32, + }, + RoundCurves { + segments: Vec, + }, + /// Another committed scene placed by `transform` (local-to-world). + /// With `transform_end`, the placement interpolates linearly (per + /// matrix element) at the ray's shutter time — transform motion blur. + /// `transform` must be invertible. + Instance { + scene: Arc, + transform: Affine3A, + transform_end: Option, + }, +} + +/// An intersection: distance, ray-facing normal (flipped to oppose the +/// ray, with `front_face` recording the original orientation), the hit +/// barycentrics where meaningful, and the IDs that let the application +/// map the hit back to its own data (materials, lights, …). +/// +/// For hits inside an [`Geometry::Instance`], `geom_id` is the +/// *instance's* id in the queried scene and `prim_id` the primitive index +/// within the instanced scene — the application maps per top-level +/// geometry. +#[derive(Clone, Copy, Debug)] +pub struct RayHit { + pub t: f32, + pub normal: Vec3A, + pub front_face: bool, + pub u: f32, + pub v: f32, + pub geom_id: u32, + pub prim_id: u32, +} + +/// Accumulates geometries, then builds the acceleration structure once in +/// [`SceneBuilder::commit`] (Embree's `rtcCommitScene`). +#[derive(Default)] +pub struct SceneBuilder { + geoms: Vec<(Geometry, u32)>, +} + +impl SceneBuilder { + pub fn new() -> Self { + Self::default() + } + + /// Attaches a geometry visible to every ray category; returns its + /// `geom_id` (dense, starting at 0 — usable as a table index). + pub fn attach(&mut self, geometry: Geometry) -> u32 { + self.attach_masked(geometry, MASK_ALL) + } + + /// Attaches a geometry visible only to ray categories in `mask`. + pub fn attach_masked(&mut self, geometry: Geometry, mask: u32) -> u32 { + self.geoms.push((geometry, mask)); + (self.geoms.len() - 1) as u32 + } + + /// Number of geometries attached so far. + pub fn count(&self) -> usize { + self.geoms.len() + } + + /// Expands every geometry into primitives and builds the BVH. + pub fn commit(self) -> Scene { + let n_geoms = self.geoms.len() as u32; + let mut prims: Vec> = Vec::new(); + for (geom_id, (geom, mask)) in self.geoms.into_iter().enumerate() { + let geom_id = geom_id as u32; + match geom { + Geometry::TriangleMesh { + vertices, + indices, + normals, + } => { + for (prim_id, tri) in indices.into_iter().enumerate() { + let [i0, i1, i2] = tri; + let (i0, i1, i2) = (i0 as usize, i1 as usize, i2 as usize); + if i0 >= vertices.len() || i1 >= vertices.len() || i2 >= vertices.len() { + continue; + } + let tri_normals = normals.as_ref().and_then(|ns| { + (i0 < ns.len() && i1 < ns.len() && i2 < ns.len()) + .then(|| [ns[i0], ns[i1], ns[i2]]) + }); + prims.push(Box::new(TrianglePrim { + v0: vertices[i0], + v1: vertices[i1], + v2: vertices[i2], + normals: tri_normals, + geom_id, + prim_id: prim_id as u32, + mask, + })); + } + } + Geometry::Sphere { center, radius } => { + prims.push(Box::new(SpherePrim { + center, + radius, + geom_id, + mask, + })); + } + Geometry::RoundCurves { segments } => { + for (prim_id, s) in segments.into_iter().enumerate() { + prims.push(Box::new(CurvePrim { + p0: s.p0, + p1: s.p1, + r0: s.r0, + r1: s.r1, + geom_id, + prim_id: prim_id as u32, + mask, + })); + } + } + Geometry::Instance { + scene, + transform, + transform_end, + } => { + let Some(inner_bounds) = scene.bounds() else { + continue; // empty instanced scene + }; + let w2l = transform.inverse(); + let bounds = match &transform_end { + Some(end) => AABB::surrounding_box( + transformed_aabb(&inner_bounds, &transform), + transformed_aabb(&inner_bounds, end), + ), + None => transformed_aabb(&inner_bounds, &transform), + }; + prims.push(Box::new(InstancePrim { + scene, + l2w: transform, + normal_mat: w2l.matrix3.transpose(), + w2l, + l2w_end: transform_end, + bounds, + geom_id, + mask, + })); + } + } + } + Scene { + bvh: Bvh::new(prims), + n_geoms, + } + } +} + +/// A committed, immutable scene. Queries are `&self` and thread-safe. +pub struct Scene { + bvh: Bvh, + n_geoms: u32, +} + +impl Scene { + /// Closest hit in `(t_min, t_max)`, or `None` (Embree's + /// `rtcIntersect1`). + pub fn intersect(&self, ray: &Ray, t_min: f32, t_max: f32) -> Option { + let hit = self.bvh.hit(ray, t_min, t_max)?; + let front_face = ray.dir.dot(hit.outward) < 0.0; + Some(RayHit { + t: hit.t, + normal: if front_face { hit.outward } else { -hit.outward }, + front_face, + u: hit.u, + v: hit.v, + geom_id: hit.geom_id, + prim_id: hit.prim_id, + }) + } + + /// Does the ray hit *anything* in `(t_min, t_max)`? Early-exit + /// traversal — the shadow-ray fast path (Embree's `rtcOccluded1`). + pub fn occluded(&self, ray: &Ray, t_min: f32, t_max: f32) -> bool { + self.bvh.hit_any(ray, t_min, t_max) + } + + /// World bounds of everything in the scene; `None` when empty. + pub fn bounds(&self) -> Option { + self.bvh.bounds() + } + + /// Number of attached geometries (`geom_id`s are `0..count`). + pub fn geometry_count(&self) -> u32 { + self.n_geoms + } + + /// Number of primitives the geometries expanded into. + pub fn primitive_count(&self) -> usize { + self.bvh.prim_count() + } + + /// Internal closest-hit that keeps the *outward* (unoriented) normal, + /// so instance transforms can map it without re-deriving orientation. + pub(crate) fn intersect_outward(&self, ray: &Ray, t_min: f32, t_max: f32) -> Option { + self.bvh.hit(ray, t_min, t_max) + } +} + +#[cfg(test)] +mod tests { + use super::*; + use crate::ray::{MASK_CAMERA, MASK_SHADOW}; + use glam::Mat4; + + fn unit_sphere_scene() -> Arc { + let mut b = SceneBuilder::new(); + b.attach(Geometry::Sphere { + center: Vec3A::ZERO, + radius: 1.0, + }); + Arc::new(b.commit()) + } + + #[test] + fn ids_map_back_to_geometries() { + let mut b = SceneBuilder::new(); + let ball = b.attach(Geometry::Sphere { + center: Vec3A::new(-3.0, 0.0, 0.0), + radius: 1.0, + }); + let quad = b.attach(Geometry::TriangleMesh { + vertices: vec![ + Vec3A::new(2.0, -1.0, -1.0), + Vec3A::new(2.0, -1.0, 1.0), + Vec3A::new(2.0, 1.0, 1.0), + Vec3A::new(2.0, 1.0, -1.0), + ], + indices: vec![[0, 1, 2], [0, 2, 3]], + normals: None, + }); + let scene = b.commit(); + assert_eq!(scene.geometry_count(), 2); + assert_eq!(scene.primitive_count(), 3); + + let hit_ball = scene + .intersect( + &Ray::new(Vec3A::new(-3.0, 0.0, -5.0), Vec3A::Z), + 0.001, + f32::INFINITY, + ) + .expect("ball hit"); + assert_eq!(hit_ball.geom_id, ball); + + // Aim at the second triangle of the quad (upper-left half). + let hit_quad = scene + .intersect( + &Ray::new(Vec3A::new(0.0, 0.5, -0.5), Vec3A::X), + 0.001, + f32::INFINITY, + ) + .expect("quad hit"); + assert_eq!(hit_quad.geom_id, quad); + assert_eq!(hit_quad.prim_id, 1); + } + + #[test] + fn front_face_semantics_match_ray_side() { + let scene = unit_sphere_scene(); + let outside = scene + .intersect(&Ray::new(Vec3A::new(0.0, 0.0, -5.0), Vec3A::Z), 0.001, 100.0) + .expect("outside hit"); + assert!(outside.front_face); + assert!(outside.normal.abs_diff_eq(-Vec3A::Z, 1e-4)); + + // From inside the sphere the normal flips toward the origin. + let inside = scene + .intersect(&Ray::new(Vec3A::ZERO, Vec3A::Z), 0.001, 100.0) + .expect("inside hit"); + assert!(!inside.front_face); + assert!(inside.normal.abs_diff_eq(-Vec3A::Z, 1e-4)); + } + + #[test] + fn masks_filter_by_ray_category() { + let mut b = SceneBuilder::new(); + b.attach_masked( + Geometry::Sphere { + center: Vec3A::ZERO, + radius: 1.0, + }, + MASK_SHADOW, + ); + let scene = b.commit(); + let ray = Ray::new(Vec3A::new(0.0, 0.0, -5.0), Vec3A::Z); + assert!(scene.intersect(&ray.with_mask(MASK_CAMERA), 0.001, 100.0).is_none()); + assert!(scene.intersect(&ray.with_mask(MASK_SHADOW), 0.001, 100.0).is_some()); + assert!(!scene.occluded(&ray.with_mask(MASK_CAMERA), 0.001, 100.0)); + assert!(scene.occluded(&ray.with_mask(MASK_SHADOW), 0.001, 100.0)); + } + + #[test] + fn smooth_normals_interpolate() { + // One triangle with vertex normals fanned outward; the hit normal + // at an interior point must be a blend, not the face normal. + let mut b = SceneBuilder::new(); + b.attach(Geometry::TriangleMesh { + vertices: vec![ + Vec3A::new(-1.0, -1.0, 2.0), + Vec3A::new(1.0, -1.0, 2.0), + Vec3A::new(0.0, 1.0, 2.0), + ], + indices: vec![[0, 1, 2]], + normals: Some(vec![ + Vec3A::new(-0.5, 0.0, -1.0).normalize(), + Vec3A::new(0.5, 0.0, -1.0).normalize(), + Vec3A::new(0.0, 0.5, -1.0).normalize(), + ]), + }); + let scene = b.commit(); + // Straight at the v1 corner region: x > 0 → normal tilts +x. + let hit = scene + .intersect(&Ray::new(Vec3A::new(0.6, -0.7, 0.0), Vec3A::Z), 0.001, 100.0) + .expect("hit"); + assert!(hit.normal.x > 0.1, "normal not interpolated: {:?}", hit.normal); + assert!(hit.normal.z < 0.0); + } + + #[test] + fn translated_instance_matches_baked() { + let mut b = SceneBuilder::new(); + b.attach(Geometry::Instance { + scene: unit_sphere_scene(), + transform: Affine3A::from_translation(glam::Vec3::new(3.0, 0.0, 0.0)), + transform_end: None, + }); + let scene = b.commit(); + let ray = Ray::new(Vec3A::new(3.0, 0.0, -5.0), Vec3A::Z); + let hit = scene.intersect(&ray, 0.001, f32::INFINITY).expect("hit"); + assert!((hit.t - 4.0).abs() < 1e-4); + assert!(hit.normal.abs_diff_eq(-Vec3A::Z, 1e-4)); + assert!(scene.occluded(&ray, 0.001, f32::INFINITY)); + assert!(!scene.occluded(&ray, 0.001, 3.9)); + } + + #[test] + fn nonuniform_scale_transforms_normals_correctly() { + // Sphere scaled 2x in X: probe an oblique point where the naive + // (non inverse-transpose) normal mapping would be wrong. + let mut b = SceneBuilder::new(); + b.attach(Geometry::Instance { + scene: unit_sphere_scene(), + transform: Affine3A::from_scale(glam::Vec3::new(2.0, 1.0, 1.0)), + transform_end: None, + }); + let scene = b.commit(); + // Hit the ellipsoid straight down above x=1 (local x=0.5). + let hit = scene + .intersect(&Ray::new(Vec3A::new(1.0, 5.0, 0.0), -Vec3A::Y), 0.001, f32::INFINITY) + .expect("hit"); + // Implicit ellipsoid (x/2)^2 + y^2 + z^2 = 1: gradient at + // (1, sqrt(3)/2, 0) is proportional to (0.5, sqrt(3), 0). + let expected = Vec3A::new(0.5, 3.0f32.sqrt(), 0.0).normalize(); + assert!( + hit.normal.abs_diff_eq(expected, 1e-3), + "normal {:?} != expected {:?}", + hit.normal, + expected + ); + } + + #[test] + fn rotated_instance_hits_where_baked_triangle_would() { + let mut inner = SceneBuilder::new(); + inner.attach(Geometry::TriangleMesh { + vertices: vec![ + Vec3A::new(-1.0, -1.0, 0.0), + Vec3A::new(1.0, -1.0, 0.0), + Vec3A::new(0.0, 1.0, 0.0), + ], + indices: vec![[0, 1, 2]], + normals: None, + }); + let mut b = SceneBuilder::new(); + b.attach(Geometry::Instance { + scene: Arc::new(inner.commit()), + transform: Affine3A::from_mat4( + Mat4::from_rotation_y(std::f32::consts::FRAC_PI_2) + * Mat4::from_translation(glam::Vec3::new(0.0, 0.0, 2.0)), + ), + transform_end: None, + }); + let scene = b.commit(); + // Local (0,0,2) maps to world (2,0,0); triangle now faces +X. + let hit = scene + .intersect(&Ray::new(Vec3A::new(5.0, 0.0, 0.0), -Vec3A::X), 0.001, f32::INFINITY) + .expect("hit"); + assert!((hit.t - 3.0).abs() < 1e-4); + assert!(hit.normal.abs_diff_eq(Vec3A::X, 1e-4)); + } + + #[test] + fn motion_blur_interpolates_position() { + let mut b = SceneBuilder::new(); + b.attach(Geometry::Instance { + scene: unit_sphere_scene(), + transform: Affine3A::IDENTITY, + transform_end: Some(Affine3A::from_translation(glam::Vec3::new(4.0, 0.0, 0.0))), + }); + let scene = b.commit(); + // At time 0 the sphere is at the origin... + let r0 = Ray::new(Vec3A::new(0.0, 0.0, -5.0), Vec3A::Z).with_time(0.0); + assert!(scene.intersect(&r0, 0.001, f32::INFINITY).is_some()); + // ...at time 1 it has moved to x=4... + let r1 = Ray::new(Vec3A::new(4.0, 0.0, -5.0), Vec3A::Z).with_time(1.0); + assert!(scene.intersect(&r1, 0.001, f32::INFINITY).is_some()); + let r1_origin = Ray::new(Vec3A::new(0.0, 0.0, -5.0), Vec3A::Z).with_time(1.0); + assert!(scene.intersect(&r1_origin, 0.001, f32::INFINITY).is_none()); + // ...and at time 0.5 it is halfway. + let rh = Ray::new(Vec3A::new(2.0, 0.0, -5.0), Vec3A::Z).with_time(0.5); + let hit = scene.intersect(&rh, 0.001, f32::INFINITY).expect("halfway hit"); + assert!((hit.t - 4.0).abs() < 1e-4); + // The shutter-union bounding box covers both endpoints. + let bb = scene.bounds().unwrap(); + assert!(bb.minimum.x <= -1.0 && bb.maximum.x >= 5.0); + } + + #[test] + fn instance_hits_report_instance_geom_id_and_inner_prim_id() { + let mut inner = SceneBuilder::new(); + inner.attach(Geometry::TriangleMesh { + vertices: vec![ + Vec3A::new(-1.0, -1.0, 0.0), + Vec3A::new(1.0, -1.0, 0.0), + Vec3A::new(1.0, 1.0, 0.0), + Vec3A::new(-1.0, 1.0, 0.0), + ], + indices: vec![[0, 1, 2], [0, 2, 3]], + normals: None, + }); + let inner = Arc::new(inner.commit()); + + let mut b = SceneBuilder::new(); + let _floor = b.attach(Geometry::Sphere { + center: Vec3A::new(0.0, -100.0, 0.0), + radius: 1.0, + }); + let inst = b.attach(Geometry::Instance { + scene: inner, + transform: Affine3A::from_translation(glam::Vec3::new(0.0, 0.0, 5.0)), + transform_end: None, + }); + let scene = b.commit(); + // Upper-left region → second triangle of the instanced mesh. + let hit = scene + .intersect(&Ray::new(Vec3A::new(-0.5, 0.5, 0.0), Vec3A::Z), 0.001, f32::INFINITY) + .expect("hit"); + assert_eq!(hit.geom_id, inst); + assert_eq!(hit.prim_id, 1); + } +} diff --git a/crates/crust-rt/src/triangle.rs b/crates/crust-rt/src/triangle.rs new file mode 100644 index 0000000..3e2b905 --- /dev/null +++ b/crates/crust-rt/src/triangle.rs @@ -0,0 +1,279 @@ +//! Watertight ray/triangle intersection and exact triangle clipping. + +use crate::aabb::AABB; +use crate::ray::Ray; +use glam::Vec3A; + +/// Watertight ray/triangle intersection (Woop, Benthin & Wald 2013, +/// "Watertight Ray/Triangle Intersection"). Returns `(t, u, v)` where `u` +/// is the barycentric weight of `v1` and `v` of `v2` (`w = 1 - u - v` of +/// `v0`), or `None` on a miss. +/// +/// The ray is transformed so its dominant axis becomes +Z (a winding- +/// preserving permutation plus a shear), the triangle is projected onto +/// z = 0, and the hit test becomes three 2D signed edge functions. Because +/// adjacent triangles share exact vertex coordinates, the shared edge's +/// function evaluates to the *same* value (with opposite sign convention) +/// for both triangles, so a ray crossing the edge can never miss both — +/// the guarantee epsilon-based Möller-Trumbore lacks. Edge functions that +/// come out exactly 0.0 in f32 are recomputed in f64, which resolves the +/// on-edge ties consistently. +pub(crate) fn triangle_intersect( + ray: &Ray, + v0: Vec3A, + v1: Vec3A, + v2: Vec3A, + t_min: f32, + t_max: f32, +) -> Option<(f32, f32, f32)> { + let d = ray.dir; + + // Permute so the dominant direction axis is Z; swap X/Y when d.z < 0 so + // the winding (and edge-function signs) are preserved. + let ad = d.abs(); + let kz = if ad.x > ad.y { + if ad.x > ad.z { 0 } else { 2 } + } else if ad.y > ad.z { + 1 + } else { + 2 + }; + let mut kx = (kz + 1) % 3; + let mut ky = (kz + 2) % 3; + if d[kz] < 0.0 { + std::mem::swap(&mut kx, &mut ky); + } + + // Shear constants mapping the ray direction onto +Z. + let sx = d[kx] / d[kz]; + let sy = d[ky] / d[kz]; + let sz = 1.0 / d[kz]; + + // Vertices relative to the ray origin, sheared into ray space. + let a = v0 - ray.origin; + let b = v1 - ray.origin; + let c = v2 - ray.origin; + let ax = a[kx] - sx * a[kz]; + let ay = a[ky] - sy * a[kz]; + let bx = b[kx] - sx * b[kz]; + let by = b[ky] - sy * b[kz]; + let cx = c[kx] - sx * c[kz]; + let cy = c[ky] - sy * c[kz]; + + // Signed 2D edge functions; e0 is opposite v0, etc. + let mut e0 = bx * cy - by * cx; + let mut e1 = cx * ay - cy * ax; + let mut e2 = ax * by - ay * bx; + + // An exact 0.0 means the ray passes through an edge or vertex in f32; + // re-evaluate in f64 so both triangles sharing the edge agree on which + // side the ray falls. + if e0 == 0.0 || e1 == 0.0 || e2 == 0.0 { + e0 = (bx as f64 * cy as f64 - by as f64 * cx as f64) as f32; + e1 = (cx as f64 * ay as f64 - cy as f64 * ax as f64) as f32; + e2 = (ax as f64 * by as f64 - ay as f64 * bx as f64) as f32; + } + + // Inside iff all edge functions share a sign (0.0 counts as both, + // keeping shared edges hittable from either triangle). + if (e0 < 0.0 || e1 < 0.0 || e2 < 0.0) && (e0 > 0.0 || e1 > 0.0 || e2 > 0.0) { + return None; + } + let det = e0 + e1 + e2; + if det == 0.0 { + return None; + } + + // Scaled hit distance; range-tested against t_min/t_max without a + // division, minding det's sign. + let az = sz * a[kz]; + let bz = sz * b[kz]; + let cz = sz * c[kz]; + let t_scaled = e0 * az + e1 * bz + e2 * cz; + if det < 0.0 && (t_scaled > t_min * det || t_scaled < t_max * det) { + return None; + } + if det > 0.0 && (t_scaled < t_min * det || t_scaled > t_max * det) { + return None; + } + + let inv_det = 1.0 / det; + Some((t_scaled * inv_det, e1 * inv_det, e2 * inv_det)) +} + +/// Exact bounds of a triangle clipped to the axis slab `[min, max]`: +/// Sutherland-Hodgman against the slab's two planes, then the bounds of +/// the surviving polygon (padded on degenerate axes like `triangle_aabb`). +/// This is what makes spatial splits effective — a long diagonal +/// triangle's *clipped* bounds shrink on every axis, where the default +/// bbox∩slab only shrinks along the split axis. +pub(crate) fn clip_triangle_aabb( + v0: Vec3A, + v1: Vec3A, + v2: Vec3A, + axis: usize, + min: f32, + max: f32, +) -> Option { + // Clip the polygon against x[axis] >= min, then x[axis] <= max. + let mut poly = [Vec3A::ZERO; 8]; + let mut n = 3; + poly[0] = v0; + poly[1] = v1; + poly[2] = v2; + + for &(bound, keep_ge) in &[(min, true), (max, false)] { + let mut out = [Vec3A::ZERO; 8]; + let mut m = 0; + for i in 0..n { + let a = poly[i]; + let b = poly[(i + 1) % n]; + let (da, db) = if keep_ge { + (a[axis] - bound, b[axis] - bound) + } else { + (bound - a[axis], bound - b[axis]) + }; + if da >= 0.0 { + out[m] = a; + m += 1; + } + if (da > 0.0) != (db > 0.0) && da != db { + // Edge crosses the plane. + let t = da / (da - db); + out[m] = a + (b - a) * t; + m += 1; + } + } + poly = out; + n = m; + if n == 0 { + return None; + } + } + + let mut lo = poly[0]; + let mut hi = poly[0]; + for &p in &poly[1..n] { + lo = lo.min(p); + hi = hi.max(p); + } + // Snap the split axis exactly to the slab (interpolation rounding can + // stick out) and pad flat axes so the slab test keeps working. + lo[axis] = lo[axis].max(min); + hi[axis] = hi[axis].min(max); + const PAD: f32 = 1e-4; + for a in 0..3 { + if hi[a] - lo[a] < PAD { + lo[a] -= PAD; + hi[a] += PAD; + } + } + Some(AABB::new(lo, hi)) +} + +#[cfg(test)] +mod tests { + use super::*; + + fn ray(o: Vec3A, d: Vec3A) -> Ray { + Ray::new(o, d) + } + + #[test] + fn interior_hit_has_expected_t_and_barycentrics() { + let (v0, v1, v2) = ( + Vec3A::new(0.0, 0.0, 5.0), + Vec3A::new(4.0, 0.0, 5.0), + Vec3A::new(0.0, 4.0, 5.0), + ); + // Aim at the barycentric point w=0.25, u=0.5, v=0.25 -> (2.0, 1.0, 5). + let r = ray(Vec3A::new(2.0, 1.0, 0.0), Vec3A::Z); + let (t, u, v) = + triangle_intersect(&r, v0, v1, v2, 0.001, f32::INFINITY).expect("interior hit"); + assert!((t - 5.0).abs() < 1e-5); + assert!((u - 0.5).abs() < 1e-5, "u={u}"); + assert!((v - 0.25).abs() < 1e-5, "v={v}"); + } + + #[test] + fn respects_t_range() { + let (v0, v1, v2) = ( + Vec3A::new(-1.0, -1.0, 5.0), + Vec3A::new(1.0, -1.0, 5.0), + Vec3A::new(0.0, 1.0, 5.0), + ); + let r = ray(Vec3A::ZERO, Vec3A::Z); + assert!(triangle_intersect(&r, v0, v1, v2, 0.001, 4.9).is_none()); + assert!(triangle_intersect(&r, v0, v1, v2, 5.1, 100.0).is_none()); + assert!(triangle_intersect(&r, v0, v1, v2, 0.001, f32::INFINITY).is_some()); + // Negative-t (behind the origin) must never report. + let back = ray(Vec3A::ZERO, -Vec3A::Z); + assert!(triangle_intersect(&back, v0, v1, v2, 0.001, f32::INFINITY).is_none()); + } + + /// The watertightness guarantee: rays crossing the shared edge of two + /// triangles forming a quad must hit at least one of them — no pinholes. + /// The diagonal is deliberately irrational-ish so the edge points do not + /// land on exact float grid values. + #[test] + fn shared_edge_is_watertight() { + let p00 = Vec3A::new(-1.3371, -0.7713, 3.7); + let p10 = Vec3A::new(1.9241, -1.1157, 4.3); + let p11 = Vec3A::new(1.6083, 1.4127, 3.9); + let p01 = Vec3A::new(-0.9743, 0.8291, 4.1); + // Shared edge p00 -> p11. + let tris = [(p00, p10, p11), (p00, p11, p01)]; + + for i in 0..=10_000 { + let s = i as f32 / 10_000.0; + let target = p00 + s * (p11 - p00); + let origin = Vec3A::new(0.1731, -0.0913, 0.0); + let r = ray(origin, target - origin); + let hits = tris + .iter() + .filter(|(a, b, c)| { + triangle_intersect(&r, *a, *b, *c, 0.001, f32::INFINITY).is_some() + }) + .count(); + assert!(hits >= 1, "pinhole at s={s} (target {target:?})"); + } + } + + /// Vertices shared by several triangles must also be covered. + #[test] + fn shared_vertex_is_covered() { + // Fan of 4 triangles around a hub vertex. + let hub = Vec3A::new(0.2137, 0.5391, 5.1); + let rim = [ + Vec3A::new(1.3, 0.4, 5.0), + Vec3A::new(0.3, 1.7, 5.3), + Vec3A::new(-1.1, 0.6, 4.9), + Vec3A::new(-0.2, -1.2, 5.2), + ]; + let r = ray(Vec3A::new(0.0, 0.0, 0.0), hub); + let hits = (0..4) + .filter(|&k| { + triangle_intersect(&r, hub, rim[k], rim[(k + 1) % 4], 0.001, f32::INFINITY) + .is_some() + }) + .count(); + assert!(hits >= 1, "ray through shared vertex missed the whole fan"); + } + + /// Axis-aligned rays (a zero direction component on the dominant-axis + /// permutation path) must still intersect correctly. + #[test] + fn axis_aligned_rays_hit() { + let (v0, v1, v2) = ( + Vec3A::new(-1.0, -1.0, 2.0), + Vec3A::new(1.0, -1.0, 2.0), + Vec3A::new(0.0, 1.0, 2.0), + ); + for d in [Vec3A::Z, -Vec3A::Z] { + let o = Vec3A::new(0.0, 0.0, if d.z > 0.0 { 0.0 } else { 4.0 }); + let hit = triangle_intersect(&ray(o, d), v0, v1, v2, 0.001, f32::INFINITY); + assert!(hit.is_some(), "axis-aligned ray {d:?} missed"); + assert!((hit.unwrap().0 - 2.0).abs() < 1e-5); + } + } +} diff --git a/docs/embree_comparison.md b/docs/embree_comparison.md new file mode 100644 index 0000000..ce6e935 --- /dev/null +++ b/docs/embree_comparison.md @@ -0,0 +1,286 @@ +# Embree vs. Crust Render — a deep comparison + +*Written against Embree 4.4.x ([RenderKit/embree](https://github.com/RenderKit/embree), +latest release 4.4.0, master documents 4.4.1) and the current state of this +repository (2026-07).* + +The single most important thing to understand before comparing the two: **they are +different kinds of software**. Embree is a *ray tracing kernel library* — it builds +acceleration structures and answers ray queries, and does nothing else. Crust is a +*complete renderer* — integrator, materials, lights, media, sampling, scene I/O — of +which the part Embree covers (primitives + BVH + traversal) is roughly 600 lines +(`bvh.rs`, `aabb.rs`, `primitives/`, `hittable.rs`). So the comparison splits into +three questions: + +1. On the ground Embree actually occupies (intersection kernels), how does crust's + implementation compare? +2. What does crust do that Embree deliberately does not? +3. What Embree ideas are worth stealing (or binding to) for crust? + +--- + +## 1. What Embree is + +Intel's open-source (Apache-2.0) collection of high-performance ray tracing kernels, +targeting professional renderers — it is the intersection backend of Blender Cycles, +Autodesk Arnold, Chaos V-Ray/Corona, DreamWorks MoonRay, and many others. It exposes a +C99 API (`rtcNewDevice` → `rtcNewScene` → `rtcNewGeometry` → `rtcCommitScene` → +`rtcIntersect*`/`rtcOccluded*`) plus a SYCL path for Intel GPUs. It ships **no +shading, no lights, no sampler, no image output** — the application owns the +integrator; Embree owns "given this ray, what does it hit, fast." + +### Geometry types + +| Type | Details | +|---|---| +| Triangle meshes | Indexed; watertight intersection (Woop et al. 2013) | +| Quad meshes | Stored as one primitive, intersected as two triangles | +| Grid meshes | Regular grids, compressed in-memory representation for heavy displaced geometry | +| Subdivision surfaces | Catmull-Clark, with edge creases, corners, holes, face-varying interpolation | +| Curves | Linear / Bézier / B-spline / Hermite / Catmull-Rom bases; each as *flat* ribbons (distant hair), *round* swept tubes (close-up), or *normal-oriented* ribbons | +| Points | Spheres, ray-oriented discs, normal-oriented discs (particles, point clouds) | +| User geometry | App-provided bounds + intersect/occluded callbacks (analytic prims, out-of-core geometry) | +| Instances | **Multi-level** instancing and instance arrays (low-overhead massive instancing) | + +### Motion blur + +Multi-segment motion blur with **2–129 time steps** per geometry, on *every* geometry +type: vertex (deformation) blur, affine instance transform blur, and **quaternion +motion blur** for correctly interpolated rotations (propellers, wheels). Rays carry a +`time` in [0,1]; the BVH stores time-bounded bounds (MBlur BVH with 4D +spatial-temporal splits). + +### Query API + +- `rtcIntersect1/4/8/16` — closest-hit, single rays or SIMD packets of 4/8/16. +- `rtcOccluded1/4/8/16` — dedicated boolean occlusion queries that **terminate on the + first hit found** — the shadow-ray fast path. +- `rtcForwardIntersect/Occluded` — continuing rays through instance boundaries in SYCL. +- **Filter functions** — per-geometry or per-ray-query callbacks invoked on *every* + candidate hit, able to accept or reject it: transparency-mapped shadows, alpha + cutouts, self-intersection avoidance, collecting *all* hits along a ray. +- **Ray masks** — per-geometry/per-ray 32-bit visibility masks (e.g. camera-only or + shadow-exempt geometry). +- **Point queries** (`rtcPointQuery`) — closest-point-on-surface within a radius + (photon lookups, signed distance). +- **Collision detection** (`rtcCollide`) — BVH-vs-BVH scene intersection tests (cloth). + +### Acceleration structures and builders + +- Internally **wide BVHs** — BVH4 / BVH8 chosen per ISA, with SIMD-parallel + child-box tests during traversal; quantized/compressed node variants + (`EMBREE_COMPACT_POLYS`, compact scene flag) for memory-constrained scenes. +- Three build qualities per geometry/scene: + - `RTC_BUILD_QUALITY_LOW` — Morton-code builder, near-linear time, for + per-frame rebuilds of dynamic geometry; + - `RTC_BUILD_QUALITY_MEDIUM` — binned SAH (the default); + - `RTC_BUILD_QUALITY_HIGH` — **spatial-split BVH (SBVH)**, presplitting large + triangles for markedly better traversal on scenes with non-uniform triangle + sizes; + - `RTC_BUILD_QUALITY_REFIT` — keep topology, refit bounds for deformations. +- Builds are **fully parallel** (TBB by default; PPL or an internal task system as + alternatives) and scale to many-core machines; 4.4 further improved build + performance when applications oversubscribe threads. +- The builder is also exposed *standalone* (`rtcBuildBVH`) so applications can build + their own BVHs over custom primitives with Embree's high-quality builders. + +### Hardware targets + +- **x86**: SSE2 → SSE4.2 → AVX → AVX2 → AVX-512, runtime ISA dispatch in one binary. +- **ARM**: NEON on Apple Silicon and AArch64 Linux (Windows-on-ARM experimental). +- **GPU (SYCL)**: Intel Xe HPG (Arc) and Xe HPC (Data Center Flex/Max), using + hardware ray tracing units; out of beta as of 4.4. CPU and GPU share the same + API and feature set (minus a few CPU-only features like collision detection). + +--- + +## 2. Crust's intersection layer, on Embree's terms + +Crust's equivalent surface is deliberately small and lives in safe Rust: + +- **Primitives**: analytic `Sphere`, `Triangle`/`SmoothTriangle`, and indexed `Mesh` + (`primitives/`). Triangle intersection is Möller–Trumbore with a + size-relative epsilon (`triangle.rs:37`) — *not* the watertight scheme; rays can + slip through shared edges/vertices in pathological cases, which Embree's default + kernels rule out by construction. +- **BVH** (`bvh.rs`): a single flat-array **binary** BVH used at both levels — a + top-level tree over scene objects (built in `Renderer::new`) and one nested tree + per imported mesh (built at USD import). Construction is binned SAH — 12 bins on + the longest centroid axis, leaf ranges of 2–8 primitives, SAH leaf-termination + test, median-split fallback — i.e. exactly Embree's `MEDIUM` quality class, + minus spatial splits. The build is **serial and recursive**; determinism is an + explicit design goal (same scene → same tree), which Embree's parallel builders + do not promise by default. +- **Traversal**: iterative with a fixed 64-slot stack, front-to-back child ordering + by split axis and ray direction sign. One ray at a time; each primitive test is a + virtual call through `Box`. SIMD exists only *within* vector math + (`glam::Vec3A` is 4-wide) — there is no packet/stream traversal and no wide-node + SIMD box test. +- **Occlusion**: there is no dedicated any-hit query. NEE shadow rays reuse the + closest-hit traversal (`shadow_transmittance`, `tracer.rs:787`), i.e. they keep + searching for the *nearest* hit when *any* hit would do. Distance-bounded + (`t_max = d − ε`) but without first-hit early-out. +- **Dynamism**: none. Transforms are baked into world-space triangles at USD import; + there is no instancing (N copies of a mesh = N copies of its triangles and its + BVH), no refit, no motion blur (rays carry no time; the camera has no shutter). +- **Parallelism**: rendering is Rayon-parallel over pixels/tiles, which keeps every + core busy at the *renderer* level even though each individual ray query is scalar. + This is the same architecture Embree-based renderers use — Embree parallelizes + *within* a query only insofar as packets do; wavefront parallelism is always the + application's job. + +### Head-to-head on the kernel ground + +*(Table updated after the adoption pass — the ✅⁺ rows were ❌ or weaker +when this document was first written; see the addendum in §4.)* + +| | Embree 4.4 | Crust | +|---|---|---| +| Language / safety | C++ (+ ISPC, SYCL), unsafe by nature | 100 % safe Rust | +| Triangles | ✅ watertight | ✅⁺ watertight (Woop 2013, f64 tie fallback) | +| Spheres / points | ✅ spheres + 2 disc types | ✅ analytic sphere (also a `LightShape`) | +| Quads / grids / subdivision | ✅ | ❌ (USD import triangulates) | +| Curves / hair | ✅ 5 bases × 3 modes | ✅⁺ round linear segments; cubic bezier/bspline/catmullRom flattened | +| User geometry | ✅ callbacks | ✅-ish (`Hittable` trait — same idea, idiomatic Rust) | +| Instancing | ✅ multi-level + arrays | ✅⁺ single-level `Instance`, shared local-space BLAS with content dedup | +| Motion blur | ✅ 2–129 steps, quaternion | ✅⁺ 2-step transform blur (linear matrix lerp) | +| BVH arity | 4–8 wide, SIMD node tests | ✅⁺ 4-wide SoA nodes, `Vec4` slab tests | +| Build quality tiers | low / medium / high(SBVH) / refit | ✅⁺ medium + SBVH spatial splits (α-gated) | +| Parallel build | ✅ TBB, many-core | ✅⁺ rayon subtree tasks, deterministic | +| Two-level structure | ✅ (scene of instances) | ✅⁺ scene of instances sharing BLASes | +| Closest-hit query | ✅ 1/4/8/16 | ✅ single ray | +| Occlusion query | ✅ early-exit `rtcOccluded` | ✅⁺ early-exit `hit_any` on every shadow ray | +| Filter / any-hit callbacks | ✅ | ❌ (no alpha-cutout shadows) | +| Ray masks | ✅ | ✅⁺ camera/shadow/indirect bits via `crust:rayMask` | +| Point queries / collision | ✅ | ❌ | +| ISA dispatch / packets | SSE2→AVX-512, NEON, SYCL GPUs | scalar rays + 4-wide node tests via glam | +| Memory options | compact/quantized nodes | one wide-node layout | + +The honest performance summary: for the scenes crust targets (a few meshes, up to a +few hundred thousand triangles), a well-built binary SAH BVH with ordered traversal is +within a small factor of Embree. The gaps that would actually show up in a profile, +in order: shadow rays doing closest-hit work (every NEE vertex), per-primitive +dynamic dispatch in leaves, binary vs. wide nodes, and serial build time on heavy +meshes. The gaps that show up in *capability* rather than speed: instancing, motion +blur, curves, and watertightness. + +--- + +## 3. What crust has that Embree does not (and never will) + +Everything above the query line — this is the part of crust that would *survive* +adopting Embree, not be replaced by it: + +- The **iterative two-pass path integrator** with MIS (selectable power/balance + heuristics + diagnostic single-strategy modes), Russian roulette, and + emission-ownership bookkeeping (`tracer.rs`). +- **OpenPBR übershader** aligned against the MaterialX/Adobe references (aniso GGX + VNDF, F82/Schlick, EON diffuse, sheen, thin-film, per-channel Cauchy dispersion, + thin-walled and thick microfacet transmission). +- **Volumetrics**: carried interior media (weighted analog free-flight with chromatic + correction) and free-standing `VolumeRegion`s (homogeneous / fBm noise / voxel + grids) with delta tracking, ratio-tracked shadow transmittance, phase MIS. Embree + has no concept of participating media at all. +- **Path guiding** (Practical Path Guiding SD-tree with variance-weighted pass + blending and a guiding-efficiency gate) — renderer research territory Embree never + touches. +- **QMC sampling** via the from-scratch `openqmc-rs` port (bit-for-bit against + upstream), domain-tree threaded through the integrator by value. +- **Adaptive sampling**, area-light NEE, **USD scene import** (with a local xformOp + composer working around an upstream `openusd` bug), EXR/PNG output. + +Embree's tutorials contain a toy pathtracer, but it is demo code; the library's +contract stops at "what did the ray hit." Everything in this list is crust's actual +value; Embree is (by design) a component such a renderer would sit on top of. + +--- + +## 4. Takeaways — what is worth adopting + +Two distinct routes, not mutually exclusive: + +> **Addendum 2 (kernel extraction).** The "openqmc move" has since been made: +> the whole intersection layer now lives in a dedicated **`crust-rt`** crate +> behind an Embree-shaped API — `Geometry` objects attached to a +> `SceneBuilder`, `commit()`, `Scene::intersect`/`Scene::occluded` +> (`rtcIntersect1`/`rtcOccluded1`), ID-based `RayHit`s (`geom_id`/`prim_id`), +> per-geometry masks. The renderer maps hits to materials through a +> `geom_id`-indexed table (`crust_core::World`), and light attribution moved +> from `Arc`-address identity to `geom_id`. This is exactly the seam Route A +> below would need: an optional Embree backend is now a matter of +> implementing the same four query/build entry points over `embree4-rs`. + +### Route A: bind to Embree + +Rust bindings exist (`embree4-sys`, `embree4-rs`). Cycles-class traversal speed, +instancing, curves, and motion blur essentially for free — at the cost of the +project's two founding constraints: **safe Rust** (an FFI kernel is unsafe C++ at the +bottom of every ray) and self-containedness (a large native dependency, ISA/build +matrix, no wasm story). Given that crust is explicitly a toy/learning renderer in +safe Rust, this route is coherent only as an optional, feature-gated backend behind +the existing `Hittable` seam — `trait` object in, `Hit` out — which the codebase's +narrow query interface (`hit(ray, t_min, t_max)`) would make straightforward to slot +in and A/B against the native BVH. + +### Route B: adopt Embree's ideas natively (recommended order) + +> **Addendum (implemented).** All seven items below have since been adopted +> natively, in safe Rust: +> 1. `Hittable::hit_any` early-exit occlusion traversal, used by every NEE +> shadow ray. +> 2. Watertight Woop-2013 triangle intersection shared by +> `Triangle`/`SmoothTriangle`, pinned by 10k-sample shared-edge tests. +> 3. `Instance` (shared local-space mesh BVHs, content-hash dedup at USD +> import). +> 4. Parallel deterministic BVH build (`rayon::join` subtrees). +> 5. BVH4: post-build collapse into 4-wide SoA nodes with `Vec4` slab tests. +> 6. SBVH spatial splits with exact triangle clipping (`clipped_aabb`). +> 7. Ray masks (`crust:rayMask`), transform motion blur +> (`crust:motion:translate` + shutter time on rays), and round curve +> primitives (`UsdGeomBasisCurves` → sphere-swept cones). +> +> The original recommendation text is kept below as the design rationale. + +1. **Dedicated occlusion query** — add `hit_any(ray, t_min, t_max) -> bool` to + `Hittable` (default-implemented via `hit`) with a real early-exit BVH traversal + (no front-to-back ordering needed, accept first confirmed hit) and use it in + `shadow_transmittance`. Cheapest change with the broadest payoff: one shadow ray + per NEE vertex on every path. +2. **Watertight triangle test** — adopt Woop-style watertight intersection (or at + minimum a shared-edge-consistent formulation) in `triangle_hit`; crust's + epsilon-based Möller–Trumbore is the classic source of pinhole leaks on shared + edges of thin geometry. +3. **Instancing** — an `Instance` `Hittable` that stores an `Arc` (sharing the + nested BVH) plus a world↔local transform, intersecting in local space. Removes the + N× memory cost of repeated USD prims and is the prerequisite for `UsdGeomPointInstancer`. +4. **Parallel BVH build** — the per-mesh builds are embarrassingly parallel across + meshes today (a Rayon `par_iter` at import time); a parallel top-down build of a + single large tree (subtree tasks below a size threshold) is the second step. + Determinism can be preserved: parallelism changes *when* subtrees are built, not + *what* is built, if splits stay deterministic. +5. **Wide BVH (BVH4)** — collapse the binary tree post-build and test 4 child boxes + with `glam`-friendly SoA layout (4×6 f32 slabs). This is the standard CPU answer + to per-node cost and works in safe Rust; measure against the added leaf/dispatch + cost before committing. +6. **Spatial splits (SBVH)** at import for meshes with high triangle-size variance — + Embree's `HIGH` quality tier; helps architectural scenes far more than the + showcase scenes currently in `samples/`. +7. Longer-term, only with a driving use case: motion blur (needs shutter-time in + `Ray` and time-bounded AABBs), curve primitives for hair, and ray masks + (trivially expressible as a bitset on `Hit`/geometry once needed). + +Packet/stream traversal is deliberately *not* on the list: crust's Rayon-over-pixels +parallelism already saturates cores, packets pay off mainly for coherent primary +rays, and incoherent path-traced bounces defeat them — Embree itself steers users +toward single-ray queries for incoherent workloads. + +--- + +## Sources + +- [RenderKit/embree](https://github.com/RenderKit/embree) — README (master, v4.4.1) +- [Embree releases](https://github.com/RenderKit/embree/releases/) — 4.4.0 latest release +- [Embree CHANGELOG](https://github.com/RenderKit/embree/blob/master/CHANGELOG.md) +- This repository: `crates/crust-core/src/{bvh,hittable,tracer,volume}.rs`, + `crates/crust-core/src/primitives/`, `CLAUDE.md`, + `docs/openpbr_reference_alignment.md` diff --git a/samples/curves.usda b/samples/curves.usda new file mode 100644 index 0000000..0422452 --- /dev/null +++ b/samples/curves.usda @@ -0,0 +1,77 @@ +#usda 1.0 +( + doc = "BasisCurves import: a tuft of cubic bezier hair strands and a linear tripod over a floor, lit by a RectLight. Exercises the RoundCurveSegment primitive and cubic flattening." + defaultPrim = "World" + upAxis = "Y" +) + +def Xform "World" +{ + def Camera "Cam" + { + float focalLength = 24 + float horizontalAperture = 20.955 + double3 xformOp:translate = (0, 1.2, 5) + uniform token[] xformOpOrder = ["xformOp:translate"] + } + + # Three cubic bezier strands curling up from the floor (vertex widths: + # thick at the root, thin at the tip). + def BasisCurves "Tuft" + { + uniform token type = "cubic" + uniform token basis = "bezier" + int[] curveVertexCounts = [4, 4, 4] + point3f[] points = [ + (-0.6, 0, 0), (-0.7, 0.6, 0.1), (-0.4, 1.2, -0.1), (-0.8, 1.7, 0), + (0.0, 0, 0.1), (0.1, 0.7, 0.0), (-0.1, 1.3, 0.2), (0.3, 1.9, 0.1), + (0.6, 0, -0.1), (0.7, 0.5, 0.0), (0.5, 1.1, 0.1), (0.9, 1.6, -0.2) + ] + float[] widths = [0.12, 0.10, 0.06, 0.03, 0.12, 0.10, 0.06, 0.03, 0.12, 0.10, 0.06, 0.03] ( + interpolation = "vertex" + ) + } + + # A linear three-segment strand, constant width. + def BasisCurves "Tripod" + { + uniform token type = "linear" + int[] curveVertexCounts = [4] + point3f[] points = [(1.6, 0, 0), (1.6, 0.8, 0.2), (1.3, 1.4, 0), (1.7, 1.8, -0.1)] + float[] widths = [0.08] ( + interpolation = "constant" + ) + double3 xformOp:translate = (0, 0, -0.5) + uniform token[] xformOpOrder = ["xformOp:translate"] + } + + def Mesh "Floor" + { + int[] faceVertexCounts = [4] + int[] faceVertexIndices = [0, 1, 2, 3] + point3f[] points = [(-8, 0, -8), (8, 0, -8), (8, 0, 8), (-8, 0, 8)] + } + + def RectLight "KeyLight" + { + float inputs:width = 2 + float inputs:height = 1.5 + color3f inputs:color = (8, 8, 8) + float inputs:intensity = 1.0 + double3 xformOp:translate = (0, 2.5, 4) + uniform token[] xformOpOrder = ["xformOp:translate"] + } +} + +def Scope "Render" +{ + def RenderSettings "settings" + { + int2 resolution = (64, 64) + int crust:samplesPerPixel = 8 + int crust:maxDepth = 4 + int crust:minSamplesPerPixel = 4 + float crust:varianceThreshold = 0.05 + int crust:frame = 0 + } +} diff --git a/samples/motionblur.usda b/samples/motionblur.usda new file mode 100644 index 0000000..d2c978b --- /dev/null +++ b/samples/motionblur.usda @@ -0,0 +1,85 @@ +#usda 1.0 +( + doc = "Transform motion blur and ray masks: a sphere and a cube streak sideways over the shutter (crust:motion:translate), and a floating blocker card is visible to shadow/indirect rays but hidden from the camera (crust:rayMask = 6)." + defaultPrim = "World" + upAxis = "Y" +) + +def Xform "World" +{ + def Camera "Cam" + { + float focalLength = 24 + float horizontalAperture = 20.955 + double3 xformOp:translate = (0, 1.5, 7) + uniform token[] xformOpOrder = ["xformOp:translate"] + } + + # Streaks one unit to the right over the shutter interval. + def Sphere "Mover" + { + double radius = 0.6 + float3 crust:motion:translate = (1.0, 0, 0) + double3 xformOp:translate = (-1.5, 0.6, 0) + uniform token[] xformOpOrder = ["xformOp:translate"] + } + + # A cube mesh rising over the shutter. + def Mesh "Riser" + { + int[] faceVertexCounts = [4, 4, 4, 4, 4, 4] + int[] faceVertexIndices = [ + 0, 1, 3, 2, 4, 6, 7, 5, 0, 4, 5, 1, 2, 3, 7, 6, 0, 2, 6, 4, 1, 5, 7, 3 + ] + point3f[] points = [ + (-0.4, -0.4, -0.4), (-0.4, -0.4, 0.4), (-0.4, 0.4, -0.4), (-0.4, 0.4, 0.4), + (0.4, -0.4, -0.4), (0.4, -0.4, 0.4), (0.4, 0.4, -0.4), (0.4, 0.4, 0.4) + ] + float3 crust:motion:translate = (0, 0.8, 0) + double3 xformOp:translate = (1.3, 0.4, 0) + uniform token[] xformOpOrder = ["xformOp:translate"] + } + + def Mesh "Floor" + { + int[] faceVertexCounts = [4] + int[] faceVertexIndices = [0, 1, 2, 3] + point3f[] points = [(-8, 0, -8), (8, 0, -8), (8, 0, 8), (-8, 0, 8)] + } + + # Visible to shadow (2) + indirect (4) rays, hidden from camera rays: + # it darkens the floor without appearing in the image. + def Mesh "ShadowCard" + { + int[] faceVertexCounts = [4] + int[] faceVertexIndices = [0, 1, 2, 3] + point3f[] points = [(-1, 0, -1), (1, 0, -1), (1, 0, 1), (-1, 0, 1)] + int crust:rayMask = 6 + double3 xformOp:translate = (0, 2.0, 0) + uniform token[] xformOpOrder = ["xformOp:translate"] + } + + def RectLight "KeyLight" + { + float inputs:width = 3 + float inputs:height = 2 + color3f inputs:color = (10, 10, 10) + float inputs:intensity = 1.0 + double3 xformOp:translate = (0, 4.5, 3) + float xformOp:rotateX = -50 + uniform token[] xformOpOrder = ["xformOp:translate", "xformOp:rotateX"] + } +} + +def Scope "Render" +{ + def RenderSettings "settings" + { + int2 resolution = (96, 54) + int crust:samplesPerPixel = 16 + int crust:maxDepth = 4 + int crust:minSamplesPerPixel = 8 + float crust:varianceThreshold = 0.05 + int crust:frame = 0 + } +} diff --git a/scripts/gen_stress_scene.py b/scripts/gen_stress_scene.py new file mode 100644 index 0000000..8e26b57 --- /dev/null +++ b/scripts/gen_stress_scene.py @@ -0,0 +1,154 @@ +#!/usr/bin/env python3 +"""Generate a traversal-heavy stress scene for benchmarking the BVH. + +Layout (all sizes in scene units, Y up): + - a 6x6 grid of identical ~9.8k-triangle UV spheres, every prim authoring + the same points/indices and binding the same material, so the importer's + instancing dedup collapses them to one shared BVH (~353k triangles when + world-baked instead); + - one 30k-triangle field of long, thin, randomly-oriented shards -- the + overlapping-bounds geometry SBVH spatial splits exist for; + - a floor, a RectLight overhead, and a thin occluder slab between the + light and the center of the field so a large fraction of NEE shadow + rays are occluded (exercising the early-exit occlusion query). + +Usage: gen_stress_scene.py [out.usda] +The scene is deliberately NOT checked in -- it is ~12 MB of generated text. +""" + +import math +import random +import sys + +RINGS = 50 # UV sphere rings +SEGS = 100 # UV sphere segments -> 9800 tris, 4902 points +GRID = 6 # GRID x GRID sphere instances +SPACING = 2.5 +SHARDS = 30000 # thin-triangle count +FIELD = 16.0 # shard field half-extent + + +def uv_sphere(radius): + pts = [(0.0, radius, 0.0)] + for r in range(1, RINGS): + phi = math.pi * r / RINGS + y = radius * math.cos(phi) + s = radius * math.sin(phi) + for k in range(SEGS): + th = 2.0 * math.pi * k / SEGS + pts.append((s * math.cos(th), y, s * math.sin(th))) + pts.append((0.0, -radius, 0.0)) + last = len(pts) - 1 + + def ring(r, k): + return 1 + (r - 1) * SEGS + (k % SEGS) + + idx = [] + for k in range(SEGS): # top fan + idx += [0, ring(1, k + 1), ring(1, k)] + for r in range(1, RINGS - 1): # quads as triangle pairs + for k in range(SEGS): + a, b = ring(r, k), ring(r, k + 1) + c, d = ring(r + 1, k), ring(r + 1, k + 1) + idx += [a, b, d, a, d, c] + for k in range(SEGS): # bottom fan + idx += [last, ring(RINGS - 1, k), ring(RINGS - 1, k + 1)] + return pts, idx + + +def fmt_pts(pts): + return ", ".join(f"({x:.5g}, {y:.5g}, {z:.5g})" for x, y, z in pts) + + +def main(): + out = sys.argv[1] if len(sys.argv) > 1 else "stress.usda" + rng = random.Random(42) + pts, idx = uv_sphere(0.8) + counts = ", ".join(["3"] * (len(idx) // 3)) + tris_per_sphere = len(idx) // 3 + + shard_pts, shard_idx = [], [] + for _ in range(SHARDS): + x = rng.uniform(-FIELD, FIELD) + z = rng.uniform(-FIELD, FIELD) + ln = rng.uniform(1.0, 3.0) + dx, dz = rng.uniform(-1, 1), rng.uniform(-1, 1) + n = math.hypot(dx, dz) or 1.0 + dx, dz = dx / n * ln, dz / n * ln + h = rng.uniform(0.3, 1.2) + b = len(shard_pts) + shard_pts += [(x, 0.0, z), (x + dx, h, z + dz), (x + dx * 0.9 + 0.03, h, z + dz * 0.9)] + shard_idx += [b, b + 1, b + 2] + shard_counts = ", ".join(["3"] * SHARDS) + + with open(out, "w") as f: + f.write('#usda 1.0\n(\n doc = "Generated BVH stress scene -- see scripts/gen_stress_scene.py"\n' + ' defaultPrim = "World"\n upAxis = "Y"\n)\n\n') + f.write('def Xform "World"\n{\n') + f.write(' def Camera "Cam"\n {\n' + ' float focalLength = 24\n' + ' float horizontalAperture = 20.955\n' + ' double3 xformOp:translate = (0, 9, 20)\n' + ' float xformOp:rotateX = -22\n' + ' uniform token[] xformOpOrder = ["xformOp:translate", "xformOp:rotateX"]\n }\n\n') + + f.write(' def Scope "Looks"\n {\n' + ' def Material "Grey"\n {\n' + ' token outputs:surface.connect = \n' + ' def Shader "Shader"\n {\n' + ' uniform token info:id = "UsdPreviewSurface"\n' + ' color3f inputs:diffuseColor = (0.55, 0.55, 0.6)\n' + ' float inputs:roughness = 0.4\n' + ' token outputs:surface\n }\n }\n }\n\n') + + half = (GRID - 1) * SPACING * 0.5 + for i in range(GRID): + for j in range(GRID): + x, z = i * SPACING - half, j * SPACING - half + f.write(f' def Mesh "Ball_{i}_{j}" (prepend apiSchemas = ["MaterialBindingAPI"])\n {{\n') + f.write(' rel material:binding = \n') + f.write(f' int[] faceVertexCounts = [{counts}]\n') + f.write(f' int[] faceVertexIndices = [{", ".join(map(str, idx))}]\n') + f.write(f' point3f[] points = [{fmt_pts(pts)}]\n') + f.write(f' double3 xformOp:translate = ({x}, 0.8, {z})\n') + f.write(' uniform token[] xformOpOrder = ["xformOp:translate"]\n }\n\n') + + f.write(' def Mesh "Shards"\n {\n') + f.write(f' int[] faceVertexCounts = [{shard_counts}]\n') + f.write(f' int[] faceVertexIndices = [{", ".join(map(str, shard_idx))}]\n') + f.write(f' point3f[] points = [{fmt_pts(shard_pts)}]\n }}\n\n') + + f.write(' def Mesh "Floor"\n {\n' + ' int[] faceVertexCounts = [4]\n' + ' int[] faceVertexIndices = [0, 1, 2, 3]\n' + ' point3f[] points = [(-40, 0, -40), (40, 0, -40), (40, 0, 40), (-40, 0, 40)]\n }\n\n') + + f.write(' def Mesh "OccluderSlab"\n {\n' + ' int[] faceVertexCounts = [4]\n' + ' int[] faceVertexIndices = [0, 1, 2, 3]\n' + ' point3f[] points = [(-8, 6, -8), (8, 6, -8), (8, 6, 8), (-8, 6, 8)]\n }\n\n') + + f.write(' def RectLight "Key"\n {\n' + ' float inputs:width = 12\n' + ' float inputs:height = 12\n' + ' color3f inputs:color = (6, 6, 6)\n' + ' float inputs:intensity = 1\n' + ' double3 xformOp:translate = (0, 12, 0)\n' + ' float xformOp:rotateX = -90\n' + ' uniform token[] xformOpOrder = ["xformOp:translate", "xformOp:rotateX"]\n }\n}\n\n') + + f.write('def Scope "Render"\n{\n def RenderSettings "settings"\n {\n' + ' int2 resolution = (640, 360)\n' + ' int crust:samplesPerPixel = 16\n' + ' int crust:maxDepth = 6\n' + ' int crust:minSamplesPerPixel = 16\n' + ' float crust:varianceThreshold = 0\n' + ' int crust:frame = 0\n }\n}\n') + + total = GRID * GRID * tris_per_sphere + SHARDS + 4 + print(f"{out}: {GRID * GRID} spheres x {tris_per_sphere} tris + {SHARDS} shards" + f" = {total} triangles (world-baked)") + + +if __name__ == "__main__": + main()