diff --git a/Cargo.lock b/Cargo.lock index d44b6fbc0..0b9ce1154 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -7685,6 +7685,7 @@ dependencies = [ "js-sys", "pollster", "rayon", + "smallvec", "tang", "thiserror 2.0.18", "vcad-ir", @@ -7951,12 +7952,15 @@ version = "0.9.3" dependencies = [ "clap", "image", + "pollster", + "rayon", "serde_json", "vcad-ecad-pcb", "vcad-ecad-symbols", "vcad-eval", "vcad-ir", "vcad-kernel", + "vcad-kernel-gpu", "vcad-kernel-raytrace", "vcad-loon", ] diff --git a/crates/vcad-kernel-raytrace/Cargo.toml b/crates/vcad-kernel-raytrace/Cargo.toml index 77a20bd0f..020e2446c 100644 --- a/crates/vcad-kernel-raytrace/Cargo.toml +++ b/crates/vcad-kernel-raytrace/Cargo.toml @@ -8,6 +8,10 @@ repository.workspace = true [dependencies] rayon = "1" +# Ray-surface intersection returns at most a handful of hits (four, on a +# torus) and is called millions of times per frame; a heap allocation per +# test is pure overhead. +smallvec = "1" vcad-kernel-math = { workspace = true } vcad-ir = { workspace = true } vcad-kernel-topo = { workspace = true } diff --git a/crates/vcad-kernel-raytrace/src/bvh.rs b/crates/vcad-kernel-raytrace/src/bvh.rs index e3b835090..5d17da099 100644 --- a/crates/vcad-kernel-raytrace/src/bvh.rs +++ b/crates/vcad-kernel-raytrace/src/bvh.rs @@ -10,13 +10,61 @@ use vcad_kernel_tessellate::TriangleMesh; use vcad_kernel_topo::FaceId; use crate::intersect::{intersect_surface, intersect_triangle, surface_tangent}; -use crate::trim::{face_normal, point_in_face}; +use crate::trim::{face_normal, FaceTrim}; use crate::{Ray, RayHit}; /// A flattened BVH node tuple for GPU upload. /// Contains: (AABB, is_leaf, left_or_first, right_or_count) pub type FlatBvhNode = (Aabb3, bool, u32, u32); +/// One triangle, resolved out of the mesh's index/vertex arrays and narrowed +/// to `f32` for GPU upload. +/// +/// The BVH stores triangles as indices into shared vertex arrays, which is how +/// you want them in memory but not how the shader reads them: the WGSL tracer +/// has no spare storage-buffer binding for a vertex array (the browser cap of +/// ten is already spent), so each triangle travels self-contained inside a +/// `GpuSurface`'s parameter block. De-indexing happens here, once, at flatten +/// time. +#[derive(Debug, Clone, Copy, PartialEq)] +pub struct FlatTriangle { + /// The three corner positions, in winding order. + pub positions: [[f32; 3]; 3], + /// Per-corner shading normals, or `None` when the source mesh carried + /// none — in which case consumers fall back to the geometric normal, the + /// same fallback [`Bvh::trace`] applies on the CPU. + pub normals: Option<[[f32; 3]; 3]>, +} + +/// What a flattened BVH's leaves index into. +/// +/// A leaf's `(first, count)` range addresses this list whichever arm it is. +/// The two arms are the two things a [`Bvh`] can be built over, so a consumer +/// learns from the value itself whether it is holding trimmed analytic faces +/// or triangles, rather than having to ask the BVH again. +#[derive(Debug, Clone)] +pub enum FlatPrims { + /// Trimmed analytic faces, in leaf order. From a BRep-backed BVH. + Faces(Vec), + /// De-indexed triangles, in leaf order. From a mesh-backed BVH. + Triangles(Vec), +} + +impl FlatPrims { + /// Number of primitives, whichever kind they are. + pub fn len(&self) -> usize { + match self { + Self::Faces(f) => f.len(), + Self::Triangles(t) => t.len(), + } + } + + /// Whether the BVH flattened to no primitives at all. + pub fn is_empty(&self) -> bool { + self.len() == 0 + } +} + /// A BVH node - either a leaf containing primitives or an internal node with /// children. #[derive(Debug, Clone)] @@ -95,6 +143,28 @@ impl MeshGeom { tri, )) } + + /// De-index one triangle into the self-contained, `f32` form the GPU + /// wants. Normals come along only when the mesh actually carried them; + /// see [`FlatTriangle::normals`]. + fn flat(&self, tri: u32) -> FlatTriangle { + let [i0, i1, i2] = self.tris[tri as usize]; + let pos = |i: u32| { + let p = self.positions[i as usize]; + [p.x as f32, p.y as f32, p.z as f32] + }; + let normals = (!self.normals.is_empty()).then(|| { + let nrm = |i: u32| { + let n = self.normals[i as usize]; + [n.x as f32, n.y as f32, n.z as f32] + }; + [nrm(i0), nrm(i1), nrm(i2)] + }); + FlatTriangle { + positions: [pos(i0), pos(i1), pos(i2)], + normals, + } + } } /// The geometry a [`Bvh`] was built over. @@ -106,6 +176,15 @@ enum BvhGeom { brep: Arc, /// Primitive index -> face ID. faces: Vec, + /// Primitive index -> that face's trim boundary, projected into UV + /// once at build time. + /// + /// Without this, every ray-face hit test reprojected the face's whole + /// trim loop — Newton-iterating each vertex back onto the surface and + /// allocating three `Vec`s — to answer one point-in-polygon query. + /// The work does not depend on the query point, and a frame asks it + /// millions of times. + trims: Vec, }, /// Triangles of a mesh-only solid. Mesh(Arc), @@ -135,6 +214,10 @@ impl Bvh { pub fn build_shared(brep: Arc) -> Self { // Collect all faces with their AABBs let faces: Vec = brep.topology.faces.iter().map(|(id, _)| id).collect(); + let trims: Vec = faces + .iter() + .map(|&face_id| FaceTrim::build(&brep, face_id)) + .collect(); let mut prim_data: Vec = faces .iter() .enumerate() @@ -152,7 +235,7 @@ impl Bvh { Self { root, - geom: BvhGeom::BRep { brep, faces }, + geom: BvhGeom::BRep { brep, faces, trims }, } } @@ -337,15 +420,14 @@ impl Bvh { /// of an any-hit query. fn prim_occludes(&self, ray: &Ray, prim: u32, t_min: f64, t_max: f64) -> bool { match &self.geom { - BvhGeom::BRep { brep, faces } => { + BvhGeom::BRep { brep, faces, trims } => { let face_id = faces[prim as usize]; + let trim = &trims[prim as usize]; let face = &brep.topology.faces[face_id]; let surface = &brep.geometry.surfaces[face.surface_index]; intersect_surface(ray, surface.as_ref()) .into_iter() - .any(|hit| { - hit.t > t_min && hit.t < t_max && point_in_face(brep, face_id, hit.uv) - }) + .any(|hit| hit.t > t_min && hit.t < t_max && trim.contains(hit.uv)) } // A triangle is convex: at most one hit, so there is nothing to // short-circuit past. @@ -438,8 +520,9 @@ impl Bvh { /// Test a ray against a single primitive, appending every hit. fn test_prim(&self, ray: &Ray, prim: u32, hits: &mut Vec) { match &self.geom { - BvhGeom::BRep { brep, faces } => { + BvhGeom::BRep { brep, faces, trims } => { let face_id = faces[prim as usize]; + let trim = &trims[prim as usize]; let face = &brep.topology.faces[face_id]; let surface = &brep.geometry.surfaces[face.surface_index]; @@ -447,7 +530,7 @@ impl Bvh { for hit in surface_hits { // Check if the hit is within the face's trim boundaries - if point_in_face(brep, face_id, hit.uv) { + if trim.contains(hit.uv) { let point = ray.at(hit.t); let normal = face_normal(brep, face_id, hit.uv); let tangent = surface_tangent(surface.as_ref(), hit.uv); @@ -470,8 +553,9 @@ impl Bvh { /// strictly past `t_min`. fn test_prim_single(&self, ray: &Ray, prim: u32, t_min: f64) -> Option { match &self.geom { - BvhGeom::BRep { brep, faces } => { + BvhGeom::BRep { brep, faces, trims } => { let face_id = faces[prim as usize]; + let trim = &trims[prim as usize]; let face = &brep.topology.faces[face_id]; let surface = &brep.geometry.surfaces[face.surface_index]; @@ -481,7 +565,7 @@ impl Bvh { for hit in surface_hits { if hit.t > t_min - && point_in_face(brep, face_id, hit.uv) + && trim.contains(hit.uv) && (closest.is_none() || hit.t < closest.as_ref().unwrap().t) { let point = ray.at(hit.t); @@ -526,26 +610,30 @@ impl Bvh { /// - For internal nodes: left_or_first = left child index, right_or_count = right child index /// - For leaf nodes: left_or_first = start face index in faces array, right_or_count = face count /// - /// Also returns the list of face IDs in leaf order. - /// - /// BRep-backed BVHs only — the GPU pipeline traces analytic surfaces. - /// A mesh-backed BVH flattens to nothing. - pub fn flatten(&self) -> (Vec, Vec) { + /// Also returns the primitives in leaf order: face IDs for a BRep-backed + /// BVH, de-indexed [`FlatTriangle`]s for a mesh-backed one. Both kinds + /// upload; the caller matches on [`FlatPrims`] to learn which it got. + pub fn flatten(&self) -> (Vec, FlatPrims) { let mut nodes = Vec::new(); - let mut faces = Vec::new(); - let BvhGeom::BRep { - faces: face_ids, .. - } = &self.geom - else { - return (nodes, faces); + let prims = match &self.geom { + BvhGeom::BRep { faces, .. } => { + let mut out = Vec::new(); + if let Some(root) = &self.root { + flatten_node(root, &mut nodes, &mut out, |p| faces[p as usize]); + } + FlatPrims::Faces(out) + } + BvhGeom::Mesh(mesh) => { + let mut out = Vec::new(); + if let Some(root) = &self.root { + flatten_node(root, &mut nodes, &mut out, |p| mesh.flat(p)); + } + FlatPrims::Triangles(out) + } }; - if let Some(root) = &self.root { - flatten_node(root, face_ids, &mut nodes, &mut faces); - } - - (nodes, faces) + (nodes, prims) } } @@ -565,19 +653,23 @@ fn get_aabb(node: &BvhNode) -> Aabb3 { } /// Recursively flatten a BVH node into a vector. -fn flatten_node( +/// +/// Generic over the primitive `resolve` produces from a build-order primitive +/// index, so the BRep (face ID) and mesh (triangle) arms share one traversal +/// rather than keeping two copies of the index bookkeeping in step. +fn flatten_node

( node: &BvhNode, - face_ids: &[FaceId], nodes: &mut Vec, - faces: &mut Vec, + prims_out: &mut Vec

, + resolve: impl Fn(u32) -> P + Copy, ) -> usize { let idx = nodes.len(); match node { BvhNode::Leaf { aabb, prims } => { - let start = faces.len() as u32; + let start = prims_out.len() as u32; let count = prims.len() as u32; - faces.extend(prims.iter().map(|&p| face_ids[p as usize])); + prims_out.extend(prims.iter().map(|&p| resolve(p))); nodes.push((*aabb, true, start, count)); } BvhNode::Internal { aabb, left, right } => { @@ -585,8 +677,8 @@ fn flatten_node( nodes.push((*aabb, false, 0, 0)); // Recursively flatten children - let left_idx = flatten_node(left, face_ids, nodes, faces); - let right_idx = flatten_node(right, face_ids, nodes, faces); + let left_idx = flatten_node(left, nodes, prims_out, resolve); + let right_idx = flatten_node(right, nodes, prims_out, resolve); // Update this node with child indices nodes[idx].2 = left_idx as u32; @@ -1022,20 +1114,110 @@ mod tests { } #[test] - fn mesh_bvh_flattens_to_nothing() { - // The GPU pipeline traces analytic surfaces only; a mesh BVH must - // not hand it face IDs it doesn't have. + fn mesh_bvh_flattens_to_triangles_in_leaf_order() { let bvh = Bvh::build_mesh(&cube_mesh()); - let (nodes, faces) = bvh.flatten(); - assert!(nodes.is_empty() && faces.is_empty()); + let (nodes, prims) = bvh.flatten(); + assert!(!nodes.is_empty()); + + let FlatPrims::Triangles(tris) = prims else { + panic!("a mesh-backed BVH must flatten to triangles"); + }; + // A cube mesh is 12 triangles and the BVH drops none of them. + assert_eq!(tris.len(), 12); + + // Every leaf's (first, count) range must address the triangle list, + // and between them the leaves must cover it exactly once — that is + // the invariant the GPU leaf walk relies on. + let mut covered = vec![0u32; tris.len()]; + for (_, is_leaf, first, count) in &nodes { + if !is_leaf { + continue; + } + for i in *first..*first + *count { + covered[i as usize] += 1; + } + } + assert!( + covered.iter().all(|&c| c == 1), + "leaf ranges must partition the triangle list, got {covered:?}" + ); + + // De-indexing must reproduce real geometry, not zeros: the cube + // spans 0..10 on every axis, so each corner coordinate is 0 or 10. + for t in &tris { + for p in &t.positions { + assert!( + p.iter().all(|c| c.abs() < 1e-4 || (c - 10.0).abs() < 1e-4), + "corner {p:?} is not on the cube" + ); + } + // Non-degenerate: the build filter already rejected zero-area + // triangles, so a collapsed one here means de-indexing is wrong. + let [a, b, c] = t.positions; + let e1 = [b[0] - a[0], b[1] - a[1], b[2] - a[2]]; + let e2 = [c[0] - a[0], c[1] - a[1], c[2] - a[2]]; + let cr = [ + e1[1] * e2[2] - e1[2] * e2[1], + e1[2] * e2[0] - e1[0] * e2[2], + e1[0] * e2[1] - e1[1] * e2[0], + ]; + assert!(cr.iter().any(|c| c.abs() > 1e-6), "degenerate {t:?}"); + } + } + + #[test] + fn mesh_without_normals_flattens_without_them() { + // A mesh with no normal array must flatten to triangles that SAY so, + // rather than to zeroed normals — the shader cannot tell an absent + // normal from a zero one, and would normalize the latter into NaN. + let mut mesh = cube_mesh(); + mesh.normals.clear(); + let bvh = Bvh::build_mesh(&mesh); + let (_, prims) = bvh.flatten(); + let FlatPrims::Triangles(tris) = prims else { + unreachable!() + }; + assert!(tris.iter().all(|t| t.normals.is_none())); + } + + #[test] + fn mesh_with_normals_carries_them_through_flatten() { + let mut mesh = cube_mesh(); + // One normal per vertex, parallel to `vertices` — the only shape + // `build_mesh` trusts. Point them all +Z; the values are arbitrary, + // what matters is that flatten preserves them per corner. + mesh.normals = (0..mesh.vertices.len() / 3) + .flat_map(|_| [0.0f32, 0.0, 1.0]) + .collect(); + let bvh = Bvh::build_mesh(&mesh); + let (_, prims) = bvh.flatten(); + let FlatPrims::Triangles(tris) = prims else { + unreachable!() + }; + assert!(!tris.is_empty()); + for t in &tris { + let ns = t.normals.expect("normals survive flatten"); + assert!(ns.iter().all(|n| *n == [0.0, 0.0, 1.0]), "{ns:?}"); + } + } + + #[test] + fn empty_mesh_bvh_flattens_to_nothing() { + let bvh = Bvh::build_mesh(&TriangleMesh::new()); + let (nodes, prims) = bvh.flatten(); + assert!(nodes.is_empty()); + assert!(prims.is_empty()); } #[test] fn brep_flatten_still_round_trips_face_ids() { let cube = make_cube(10.0, 10.0, 10.0); let bvh = Bvh::build(&cube); - let (nodes, faces) = bvh.flatten(); + let (nodes, prims) = bvh.flatten(); assert!(!nodes.is_empty()); + let FlatPrims::Faces(faces) = prims else { + panic!("a BRep-backed BVH must flatten to face IDs"); + }; assert_eq!(faces.len(), cube.topology.faces.len()); assert!(faces.iter().all(|&f| cube.topology.faces.contains_key(f))); } diff --git a/crates/vcad-kernel-raytrace/src/gpu/buffers.rs b/crates/vcad-kernel-raytrace/src/gpu/buffers.rs index f33e6518b..68cb2b7b2 100644 --- a/crates/vcad-kernel-raytrace/src/gpu/buffers.rs +++ b/crates/vcad-kernel-raytrace/src/gpu/buffers.rs @@ -4,6 +4,7 @@ use bytemuck::{Pod, Zeroable}; use vcad_kernel_booleans::bbox::face_aabb; use vcad_kernel_geom::{Surface, SurfaceKind}; use vcad_kernel_primitives::BRepSolid; +use vcad_kernel_tessellate::TriangleMesh; use vcad_kernel_topo::FaceId; use crate::bvh::Bvh; @@ -16,11 +17,42 @@ pub const MAX_SURFACES: usize = 1024; pub const MAX_FACES: usize = 4096; /// Maximum BVH nodes. -pub const MAX_BVH_NODES: usize = 8192; +/// +/// Raised from the original 8192 when `vcad-render --photoreal --gpu` started +/// merging one mesh BLAS per solid into a single tree: a real assembly is +/// hundreds of thousands of triangles and its BVH has roughly one node per +/// two of them, so 8192 refused every scene bigger than a bracket. Nothing in +/// the shader is sized by this — the node array is a storage buffer sized +/// from the data — so the ceiling that actually matters is the device's +/// `max_storage_buffer_binding_size`, which [`GpuScene::validate`] checks +/// separately. This is the "obviously absurd" guard, not the real limit. +pub const MAX_BVH_NODES: usize = 4_000_000; + +/// Deepest root-to-leaf path the WGSL tracer can walk. +/// +/// `trace_bvh` holds its traversal stack in a fixed `array` and +/// *silently drops* a push that would overflow it — geometry simply +/// disappears from the render. [`GpuScene::validate`] measures the packed +/// tree against this so an over-deep scene is a message instead of a hole. +pub const MAX_TRAVERSAL_DEPTH: usize = 64; /// Maximum trim loop vertices. pub const MAX_TRIM_VERTS: usize = 32768; +/// Maximum triangles in a single mesh-backed scene. +/// +/// Deliberately far above [`MAX_FACES`]: those caps are sized for a BRep, +/// where a "face" is a whole trimmed surface and a few thousand is a large +/// part. A tessellated mesh counts in the hundreds of thousands for the same +/// object, so the mesh path gets its own ceiling. At +/// [`size_of::`](GpuSurface) = 144 bytes per triangle this bounds +/// the surface buffer at ~144 MB, which is above the 128 MB that +/// `maxStorageBufferBindingSize` defaults to — so a mesh near this cap can +/// still be refused by the driver. The limit is here to turn an absurd mesh +/// into a message rather than an OOM, not to certify that everything under it +/// uploads. +pub const MAX_MESH_TRIANGLES: usize = 1_000_000; + /// GPU-compatible surface representation. /// /// Each surface type is packed into 32 floats: @@ -29,7 +61,8 @@ pub const MAX_TRIM_VERTS: usize = 32768; #[repr(C)] #[derive(Clone, Copy, Debug, Pod, Zeroable)] pub struct GpuSurface { - /// Surface type: 0=Plane, 1=Cylinder, 2=Sphere, 3=Cone, 4=Torus, 5=Bilinear + /// Surface type: 0=Plane, 1=Cylinder, 2=Sphere, 3=Cone, 4=Torus, + /// 5=Bilinear, 6=BSpline, 7=Triangle. See [`GpuSurface::type_name`]. pub surface_type: u32, /// Padding for alignment pub _pad: [u32; 3], @@ -171,8 +204,93 @@ impl GpuSurface { params, } } + + /// Pack one mesh triangle into a surface record. + /// + /// A triangle is not a parametric surface, but the `params` block is 32 + /// floats of otherwise-idle space and a triangle needs 19 of them, so it + /// rides in the surface array rather than in a vertex buffer of its own. + /// That choice is what keeps the mesh path inside the browser's cap of + /// ten storage-buffer bindings — see the binding block at the top of + /// `shaders/raytrace.wgsl`. The cost is that shared vertices are stored + /// once per incident triangle: 144 bytes per triangle flat, so a 500k-tri + /// mesh is ~72 MB of surface buffer. On unified memory that is a fair + /// trade for not forking the shader; a native-only indexed path would be + /// roughly 4x leaner and is the obvious follow-up if it ever bites. + /// + /// Layout, mirrored by `intersect_triangle` and `compute_normal` in the + /// WGSL: + /// + /// ```text + /// [0..3) v0 [9..12) n0 + /// [3..6) v1 [12..15) n1 + /// [6..9) v2 [15..18) n2 + /// [18] 1.0 when the normals above are real, 0.0 otherwise + /// ``` + pub fn triangle(tri: &crate::bvh::FlatTriangle) -> Self { + let mut params = [0.0f32; 32]; + for (i, p) in tri.positions.iter().enumerate() { + params[i * 3..i * 3 + 3].copy_from_slice(p); + } + // Absent normals stay zero AND are flagged, because zero is a + // legitimate-looking vector the shader would happily normalize into + // NaN. The flag is what makes the shader take the geometric-normal + // fallback instead — the same fallback `MeshGeom::test` takes on the + // CPU, so a normal-less mesh shades identically in both renderers. + if let Some(normals) = &tri.normals { + for (i, n) in normals.iter().enumerate() { + params[9 + i * 3..12 + i * 3].copy_from_slice(n); + } + params[18] = 1.0; + } + + Self { + surface_type: SURFACE_TYPE_TRIANGLE, + _pad: [0; 3], + params, + } + } + + /// Whether the WGSL `intersect_surface` switch has a case for this + /// surface type. + /// + /// Types 0-4 (plane, cylinder, sphere, cone, torus) are traced + /// analytically and type 7 (triangle) by Möller-Trumbore. Bilinear (5) + /// and B-spline (6) are packed by [`Self::from_surface`] but fall into + /// the shader's `default` arm, which returns a miss — such a face would + /// silently vanish from the render, so [`GpuScene::from_brep`] rejects a + /// scene containing one. + pub fn is_gpu_traceable(&self) -> bool { + self.surface_type <= SURFACE_TYPE_TORUS || self.surface_type == SURFACE_TYPE_TRIANGLE + } + + /// Human-readable name of a packed surface type code. + pub fn type_name(surface_type: u32) -> &'static str { + match surface_type { + 0 => "Plane", + 1 => "Cylinder", + 2 => "Sphere", + 3 => "Cone", + 4 => "Torus", + 5 => "Bilinear", + 6 => "BSpline", + 7 => "Triangle", + _ => "Unknown", + } + } } +/// Highest surface type code the WGSL `intersect_surface` switch handles +/// analytically. +const SURFACE_TYPE_TORUS: u32 = 4; + +/// Surface type code for a mesh triangle. +/// +/// Deliberately past B-spline (6) rather than reusing a hole: the codes 0-6 +/// are a one-to-one image of [`SurfaceKind`] and a triangle is not one of +/// those, so it gets its own code and the mapping stays honest. +pub const SURFACE_TYPE_TRIANGLE: u32 = 7; + /// GPU-compatible material representation (PBR). /// /// Mirrors [`crate::pathtrace::Pbr`] field-for-field so the GPU path tracer and @@ -384,16 +502,39 @@ pub struct GpuCamera { pub target: [f32; 4], /// Up vector. pub up: [f32; 4], + /// World direction mapping to screen +x. Only read when + /// [`basis_mode`](Self::basis_mode) is [`CAMERA_BASIS_EXPLICIT`]. + pub right: [f32; 4], /// Field of view in radians. pub fov: f32, /// Image width. pub width: u32, /// Image height. pub height: u32, - /// Padding. - pub _pad: u32, + /// How the shader builds the screen basis. + /// + /// * [`CAMERA_BASIS_DERIVED`] — build it right-handedly from `position`, + /// `target` and `up` (`right = forward × up`). What the viewport has + /// always done, and what [`GpuCamera::new`] still sets. + /// * [`CAMERA_BASIS_EXPLICIT`] — use [`right`](Self::right) and + /// [`up`](Self::up) *verbatim*, with `forward = normalize(target - + /// position)`. + /// + /// The explicit mode exists because [`crate::pathtrace::Camera`] can + /// carry a **mirrored** (left-handed) screen basis — `View::Isometric` + /// and the named CAD views in `vcad-render` all do — and no + /// `look_at`-plus-up-hint construction can reproduce one. Rebuilding such + /// a view right-handedly flips the image left-for-right. + pub basis_mode: u32, } +/// [`GpuCamera::basis_mode`]: derive the screen basis from the up hint. +pub const CAMERA_BASIS_DERIVED: u32 = 0; + +/// [`GpuCamera::basis_mode`]: use the supplied `right`/`up` verbatim, so a +/// mirrored basis survives the trip to the shader. +pub const CAMERA_BASIS_EXPLICIT: u32 = 1; + /// Render state for progressive rendering. /// /// Layout (128 bytes, 16-byte aligned — matches `RenderState` in raytrace.wgsl): @@ -473,10 +614,50 @@ pub struct GpuRenderState { pub env_rotation: f32, /// Normaliser for the environment's uv-space PDF. pub env_marg_int: f32, + /// Extra decorrelation term folded into the WGSL per-pixel RNG seed. + /// + /// The shader's hash is `pixel.x*1973 + pixel.y*9277 + sample*26699 + + /// frame_index*12345 + seed*2654435761 + 1`. Zero — the value the + /// viewport uses and the value every existing constructor sets — + /// reproduces the pre-seed behaviour bit for bit, so the browser path is + /// unchanged. An offline render sets it to get a different but + /// *reproducible* sample sequence for the same frame indices. + /// + /// Occupies what used to be the first padding word, so the uniform is + /// still 128 bytes. + pub seed: u32, + /// What a camera ray that hits nothing returns. + /// + /// * [`BACKGROUND_SKY`] (0) — `sky_color`, the *themed viewport backdrop*. + /// The historical behaviour, and what every viewport constructor sets. + /// * [`BACKGROUND_ENVIRONMENT`] (1) — `env_radiance`, the same sky the + /// integrator lights with. This is what `vcad-render --photoreal` + /// shows behind the subject, so an offline render asks for it. + /// + /// It has to be a shader-side choice rather than a CPU composite: a pixel + /// on the subject's silhouette averages background and surface samples + /// together, and once that mean exists the two contributions cannot be + /// separated again. + /// + /// Occupies what used to be the first word of `_pad3`, so the uniform is + /// still 128 bytes. + pub background_mode: u32, /// Padding to a 16-byte multiple (required for uniform buffers). - pub _pad3: [u32; 3], + pub _pad3: u32, } +/// [`GpuRenderState::background_mode`]: draw the themed viewport backdrop. +pub const BACKGROUND_SKY: u32 = 0; + +/// [`GpuRenderState::background_mode`]: draw the lighting environment, as the +/// CPU renderer does with `PathTraceOptions::show_background`. +pub const BACKGROUND_ENVIRONMENT: u32 = 1; + +/// [`GpuRenderState::background_mode`]: leave the backdrop black, matching +/// the CPU renderer with `show_background` off. Paired with the film's +/// coverage alpha this is what makes a transparent PNG. +pub const BACKGROUND_BLACK: u32 = 2; + /// Default silhouette line color: near-black, slightly cool. const DEFAULT_SILHOUETTE_COLOR: [f32; 4] = [0.08, 0.08, 0.10, 1.0]; /// Default crease line color: slightly lighter than silhouette. @@ -550,7 +731,9 @@ impl GpuRenderState { env_height: 0, env_rotation: 0.0, env_marg_int: 0.0, - _pad3: [0; 3], + seed: 0, + background_mode: BACKGROUND_SKY, + _pad3: 0, } } @@ -664,7 +847,9 @@ impl GpuRenderState { env_height: 0, env_rotation: 0.0, env_marg_int: 0.0, - _pad3: [0; 3], + seed: 0, + background_mode: BACKGROUND_SKY, + _pad3: 0, } } @@ -690,6 +875,15 @@ impl GpuRenderState { } } +/// Sub-pixel jitter for one accumulation frame, in `[-0.5, 0.5]`. +/// +/// The same low-discrepancy offsets [`GpuRenderState::new`] bakes in, exposed +/// so the offline sample loop can advance the jitter without rebuilding the +/// whole render state each sample. +pub fn halton_jitter(frame_index: u32) -> (f32, f32) { + halton_2_3(frame_index) +} + /// Generate Halton sequence sample for bases 2 and 3. /// Returns values in range [-0.5, 0.5] for sub-pixel jittering. fn halton_2_3(index: u32) -> (f32, f32) { @@ -723,10 +917,49 @@ impl GpuCamera { position: [position[0], position[1], position[2], 1.0], target: [target[0], target[1], target[2], 1.0], up: [up[0], up[1], up[2], 0.0], + // Unread in derived mode; zero rather than a made-up axis so a + // stale value can never be mistaken for a real basis. + right: [0.0; 4], fov, width, height, - _pad: 0, + basis_mode: CAMERA_BASIS_DERIVED, + } + } + + /// Create a camera from an explicit — possibly mirrored — screen basis. + /// + /// `forward`, `right` and `up` are used as given (the shader normalises + /// them but does not re-orthogonalise), so a left-handed CAD view reaches + /// the GPU unflipped. `focus_dist` only positions the `target` point the + /// shader derives `forward` from; it does not focus anything, since the + /// GPU tracer is a pinhole. + #[allow(clippy::too_many_arguments)] + pub fn from_basis( + position: [f32; 3], + forward: [f32; 3], + right: [f32; 3], + up: [f32; 3], + fov: f32, + focus_dist: f32, + width: u32, + height: u32, + ) -> Self { + let d = focus_dist.max(1.0); + Self { + position: [position[0], position[1], position[2], 1.0], + target: [ + position[0] + forward[0] * d, + position[1] + forward[1] * d, + position[2] + forward[2] * d, + 1.0, + ], + up: [up[0], up[1], up[2], 0.0], + right: [right[0], right[1], right[2], 0.0], + fov, + width, + height, + basis_mode: CAMERA_BASIS_EXPLICIT, } } } @@ -773,6 +1006,34 @@ pub enum GpuSceneError { TooManyBvhNodes(usize), /// Too many trim vertices. TooManyTrimVerts(usize), + /// A surface kind the WGSL tracer has no intersection case for. + /// + /// Carries the surface's index in `brep.geometry.surfaces` and the packed + /// type name, so the caller can say *which* geometry it cannot render + /// instead of handing back a blank frame. + UnsupportedSurface { + /// Index into `brep.geometry.surfaces`. + index: usize, + /// Packed surface type code. + surface_type: u32, + /// Human-readable name of that type. + name: &'static str, + }, + /// Too many mesh triangles (exceeds [`MAX_MESH_TRIANGLES`]). + TooManyTriangles(usize), + /// [`GpuScene::from_mesh_bvh`] was handed a BRep-backed BVH. + NotAMeshBvh, + /// The packed BVH is deeper than the shader's traversal stack + /// ([`MAX_TRAVERSAL_DEPTH`]). + BvhTooDeep(usize), + /// The largest storage buffer this scene needs is bigger than the + /// device's `max_storage_buffer_binding_size`. + ExceedsDeviceBinding { + /// Bytes the largest single binding would need. + bytes: u64, + /// The device's limit. + cap: u64, + }, } impl std::fmt::Display for GpuSceneError { @@ -788,6 +1049,43 @@ impl std::fmt::Display for GpuSceneError { Self::TooManyTrimVerts(n) => { write!(f, "too many trim vertices: {} (max {})", n, MAX_TRIM_VERTS) } + Self::UnsupportedSurface { + index, + surface_type, + name, + } => write!( + f, + "surface {index} is a {name} (type {surface_type}), which the GPU tracer \ + cannot intersect; faces on it would render as empty space" + ), + Self::TooManyTriangles(n) => write!( + f, + "too many mesh triangles: {n} (max {MAX_MESH_TRIANGLES}) -- \ + each triangle costs one {tri_bytes}-byte surface record, so \ + this mesh would need {mb} MB of surface buffer alone", + tri_bytes = std::mem::size_of::(), + mb = n * std::mem::size_of::() / (1024 * 1024), + ), + Self::BvhTooDeep(d) => write!( + f, + "merged BVH is {d} levels deep (max {MAX_TRAVERSAL_DEPTH}) -- the \ + GPU tracer's traversal stack cannot hold it and would silently \ + drop geometry. Render fewer parts per pass, or fall back to the \ + CPU tracer (drop --gpu)" + ), + Self::ExceedsDeviceBinding { bytes, cap } => write!( + f, + "scene needs a {} MB storage binding but this adapter caps one at \ + {} MB -- split the render or drop --gpu to trace on the CPU", + bytes / (1024 * 1024), + cap / (1024 * 1024), + ), + Self::NotAMeshBvh => write!( + f, + "from_mesh_bvh needs a BVH built by Bvh::build_mesh; this one is \ + BRep-backed. Use GpuScene::from_brep, which traces its surfaces \ + analytically rather than a tessellation of them" + ), } } } @@ -799,6 +1097,92 @@ impl std::error::Error for GpuSceneError {} /// Delegates to [`crate::pathtrace::studio_rig`] — the SAME function /// `vcad-render --photoreal` calls — so the viewport and the CPU renderer are /// lit by an identical rig rather than by two hand-tuned approximations. +/// Convert a flattened BVH into the shader's node layout. +/// +/// An empty tree still yields one (zeroed) node: WebGPU rejects a zero-sized +/// storage buffer, and a zero AABB is missed by every ray, so the empty scene +/// renders as pure background rather than failing to bind. +fn gpu_bvh_nodes(flat_nodes: &[crate::bvh::FlatBvhNode]) -> Vec { + if flat_nodes.is_empty() { + return vec![GpuBvhNode::zeroed()]; + } + flat_nodes + .iter() + .map( + |(aabb, is_leaf, left_or_first, right_or_count)| GpuBvhNode { + aabb_min: [aabb.min.x as f32, aabb.min.y as f32, aabb.min.z as f32, 0.0], + aabb_max: [aabb.max.x as f32, aabb.max.y as f32, aabb.max.z as f32, 0.0], + // For leaves `left_or_first` is a start index into `faces`, which + // is built in BVH leaf order, so it maps across directly. + left_or_first: *left_or_first, + right_or_count: *right_or_count, + is_leaf: u32::from(*is_leaf), + _pad: 0, + }, + ) + .collect() +} + +/// A triangle moved into world space by `t`. +/// +/// Positions go through the full matrix; shading normals through +/// `apply_normal`, which is the inverse-transpose — under a non-uniform scale +/// a normal does *not* transform like the surface it sits on, and using the +/// plain vector transform would light a squashed part as if it were not. +/// Zero-length results are left alone so the shader's +/// "normals are real" flag keeps meaning what it says. +fn place_triangle( + tri: &crate::bvh::FlatTriangle, + t: &vcad_kernel_math::Transform, +) -> crate::bvh::FlatTriangle { + use vcad_kernel_math::{Point3, Vec3}; + + let mut out = *tri; + for p in out.positions.iter_mut() { + let w = t.apply_point(&Point3::new(p[0] as f64, p[1] as f64, p[2] as f64)); + *p = [w.x as f32, w.y as f32, w.z as f32]; + } + if let Some(normals) = out.normals.as_mut() { + for n in normals.iter_mut() { + let w = t.apply_normal(&Vec3::new(n[0] as f64, n[1] as f64, n[2] as f64)); + let len = w.norm(); + if len > 1e-12 { + *n = [(w.x / len) as f32, (w.y / len) as f32, (w.z / len) as f32]; + } + } + } + out +} + +/// Re-fit a packed BVH node's AABB around the eight transformed corners of +/// its old one. Conservative under rotation — the new box is axis-aligned +/// around a rotated box — which costs some traversal and can never miss a +/// primitive the old box contained. +fn place_aabb(node: &mut GpuBvhNode, t: &vcad_kernel_math::Transform) { + use vcad_kernel_math::Point3; + + let (lo, hi) = (node.aabb_min, node.aabb_max); + if !lo[..3].iter().chain(&hi[..3]).all(|v| v.is_finite()) { + return; + } + let mut min = [f32::INFINITY; 3]; + let mut max = [f32::NEG_INFINITY; 3]; + for i in 0..8 { + let c = Point3::new( + if i & 1 == 0 { lo[0] } else { hi[0] } as f64, + if i & 2 == 0 { lo[1] } else { hi[1] } as f64, + if i & 4 == 0 { lo[2] } else { hi[2] } as f64, + ); + let w = t.apply_point(&c); + for (a, v) in [w.x, w.y, w.z].into_iter().enumerate() { + min[a] = min[a].min(v as f32); + max[a] = max[a].max(v as f32); + } + } + node.aabb_min = [min[0], min[1], min[2], 0.0]; + node.aabb_max = [max[0], max[1], max[2], 0.0]; +} + fn studio_lights_for_bvh(bvh_nodes: &[GpuBvhNode]) -> Vec { let Some(root) = bvh_nodes.first() else { return Vec::new(); @@ -877,9 +1261,27 @@ impl GpuScene { return Err(GpuSceneError::TooManySurfaces(surfaces.len())); } + // Fail closed on surface kinds the WGSL tracer has no case for. + // Packing them and uploading anyway makes those faces disappear from + // the image with no diagnostic at all. + if let Some((index, s)) = surfaces + .iter() + .enumerate() + .find(|(_, s)| !s.is_gpu_traceable()) + { + return Err(GpuSceneError::UnsupportedSurface { + index, + surface_type: s.surface_type, + name: GpuSurface::type_name(s.surface_type), + }); + } + // Build BVH first to get the face ordering let bvh = Bvh::build(brep); - let (flat_nodes, bvh_faces) = bvh.flatten(); + let (flat_nodes, prims) = bvh.flatten(); + let crate::bvh::FlatPrims::Faces(bvh_faces) = prims else { + unreachable!("Bvh::build always produces a BRep-backed BVH") + }; // Build face list in BVH traversal order (so BVH leaf indices are contiguous) let mut faces = Vec::with_capacity(bvh_faces.len()); @@ -1012,35 +1414,7 @@ impl GpuScene { // Convert flattened BVH to GPU format // Faces are now in BVH order, so leaf indices map directly - let mut bvh_nodes = Vec::with_capacity(flat_nodes.len().max(1)); - - if flat_nodes.is_empty() { - // Empty BVH - add a dummy node - bvh_nodes.push(GpuBvhNode::zeroed()); - } else { - for (aabb, is_leaf, left_or_first, right_or_count) in &flat_nodes { - if *is_leaf { - // For leaves: left_or_first is start index in faces array (which is now BVH-ordered) - bvh_nodes.push(GpuBvhNode { - aabb_min: [aabb.min.x as f32, aabb.min.y as f32, aabb.min.z as f32, 0.0], - aabb_max: [aabb.max.x as f32, aabb.max.y as f32, aabb.max.z as f32, 0.0], - left_or_first: *left_or_first, - right_or_count: *right_or_count, - is_leaf: 1, - _pad: 0, - }); - } else { - bvh_nodes.push(GpuBvhNode { - aabb_min: [aabb.min.x as f32, aabb.min.y as f32, aabb.min.z as f32, 0.0], - aabb_max: [aabb.max.x as f32, aabb.max.y as f32, aabb.max.z as f32, 0.0], - left_or_first: *left_or_first, - right_or_count: *right_or_count, - is_leaf: 0, - _pad: 0, - }); - } - } - } + let bvh_nodes = gpu_bvh_nodes(&flat_nodes); if bvh_nodes.len() > MAX_BVH_NODES { return Err(GpuSceneError::TooManyBvhNodes(bvh_nodes.len())); @@ -1069,6 +1443,148 @@ impl GpuScene { }) } + /// Build GPU scene data from a triangle mesh. + /// + /// The counterpart of [`Self::from_brep`] for geometry that has no + /// analytic surfaces: frozen `topology_optimize` results, imported + /// STL/GLB parts, and the cached tessellations `--photoreal` traces by + /// default. Builds the BVH with [`Bvh::build_mesh`] — the *same* BLAS the + /// CPU path tracer walks — so the two renderers agree on geometry and + /// differ only in arithmetic precision. + pub fn from_mesh(mesh: &TriangleMesh) -> Result { + Self::from_mesh_bvh(&Bvh::build_mesh(mesh)) + } + + /// Build GPU scene data from an already-built mesh BVH. + /// + /// Split out from [`Self::from_mesh`] so a caller that already traces the + /// BVH on the CPU can upload that exact tree rather than rebuilding a + /// second one that might partition differently. + /// + /// Rejects a BRep-backed BVH with [`GpuSceneError::NotAMeshBvh`] rather + /// than silently producing an empty scene. + /// + /// Every face gets material 0, the GPU default grey. Use + /// [`Self::from_mesh_bvh_placed`] to give the part its own material and + /// world placement, which is what a multi-part scene needs. + pub fn from_mesh_bvh(bvh: &Bvh) -> Result { + Self::from_mesh_bvh_placed(bvh, GpuMaterial::default(), None) + } + + /// [`Self::from_mesh_bvh`] with the part's own material and world + /// transform. + /// + /// **The transform is baked into the packed vertices**, not carried as an + /// instance: the WGSL tracer walks one flat node array with no instancing + /// layer, so there is nowhere to put a per-object matrix. Positions go + /// through `apply_point`, shading normals through `apply_normal` (the + /// inverse-transpose — a non-uniform scale rotates a normal differently + /// from the surface it belongs to, and using `apply_vec` here would shade + /// a squashed part wrong), and the BVH node AABBs are re-fitted around + /// their transformed corners. Re-fitting corners is conservative rather + /// than tight under rotation, which costs a little traversal and cannot + /// lose a hit. + /// + /// Baking means the result is a *static* snapshot: a new pose needs a new + /// scene. That is why `--animate` stays on the CPU. + pub fn from_mesh_bvh_placed( + bvh: &Bvh, + material: GpuMaterial, + transform: Option<&vcad_kernel_math::Transform>, + ) -> Result { + let mut scene = Self::pack_mesh_bvh(bvh, transform)?; + scene.materials = vec![material]; + Ok(scene) + } + + fn pack_mesh_bvh( + bvh: &Bvh, + transform: Option<&vcad_kernel_math::Transform>, + ) -> Result { + let (flat_nodes, prims) = bvh.flatten(); + let crate::bvh::FlatPrims::Triangles(tris) = prims else { + return Err(GpuSceneError::NotAMeshBvh); + }; + + if tris.len() > MAX_MESH_TRIANGLES { + return Err(GpuSceneError::TooManyTriangles(tris.len())); + } + + // One surface and one face per triangle, in BVH leaf order, so a + // leaf's (first, count) range indexes `faces` directly — exactly as + // in the BRep path. The face carries no trim loops: a triangle's + // Möller-Trumbore test already answers the containment question that + // trimming answers for an analytic surface, and `point_in_face` + // short-circuits on the triangle type code rather than running a + // winding test over an empty polygon. + let mut surfaces = Vec::with_capacity(tris.len()); + let mut faces = Vec::with_capacity(tris.len()); + for (i, tri) in tris.iter().enumerate() { + let placed; + let tri = match transform { + None => tri, + Some(t) => { + placed = place_triangle(tri, t); + &placed + } + }; + surfaces.push(GpuSurface::triangle(tri)); + + let mut lo = tri.positions[0]; + let mut hi = tri.positions[0]; + for p in &tri.positions[1..] { + for a in 0..3 { + lo[a] = lo[a].min(p[a]); + hi[a] = hi[a].max(p[a]); + } + } + + faces.push(GpuFace { + surface_idx: i as u32, + // Triangle normals come from the mesh's own winding and + // vertex normals; there is no topological orientation to + // apply on top, and flipping here would invert the shading + // relative to the CPU tracer. + orientation: 0, + trim_start: 0, + trim_count: 0, + aabb_min: [lo[0], lo[1], lo[2], 0.0], + aabb_max: [hi[0], hi[1], hi[2], 0.0], + inner_start: 0, + inner_count: 0, + inner_loop_count: 0, + inner_desc_start: 0, + material_idx: 0, + _pad2: [0; 3], + }); + } + + let mut bvh_nodes = gpu_bvh_nodes(&flat_nodes); + if let Some(t) = transform { + for n in bvh_nodes.iter_mut() { + place_aabb(n, t); + } + } + let lights = studio_lights_for_bvh(&bvh_nodes); + + Ok(Self { + surfaces, + faces, + materials: vec![GpuMaterial::default()], + bvh_nodes, + // WebGPU refuses a zero-sized storage buffer, and a mesh scene + // has nothing to put in either of these. One dummy element each, + // referenced by nothing (every face has trim_count 0). + trim_verts: vec![GpuVec2 { x: 0.0, y: 0.0 }], + inner_loop_descs: vec![0], + // Face IDs are a BRep concept; a triangle has none, so nothing + // here can be keyed by one. + face_index_map: std::collections::HashMap::new(), + environment: None, + lights, + }) + } + /// Merge another GpuScene into this one. Combines surfaces, faces, /// materials, trim verts, inner-loop descriptors, and BVH nodes. /// @@ -1162,6 +1678,12 @@ impl GpuScene { merged_bvh.extend(self.bvh_nodes); merged_bvh.extend(adjusted_nodes); + // Re-derive the studio rig from the combined bounds. Keeping self's + // rig would light the merged scene as if only self's half of it + // existed — with a mesh part merged alongside a BRep one, the softbox + // distances come out wrong for whichever half did not set them. + self.lights = studio_lights_for_bvh(&merged_bvh); + self.surfaces.extend(other.surfaces); self.faces.extend(adjusted_faces); self.materials.extend(other.materials); @@ -1172,6 +1694,122 @@ impl GpuScene { self } + /// Fold many scenes into one with a **balanced** tree of merges. + /// + /// [`Self::merge`] adds exactly one level of depth per call, so folding N + /// parts linearly (`a.merge(b).merge(c)…`) costs N-1 levels of traversal + /// stack on top of the deepest part's own tree. A pairwise fold costs + /// `ceil(log2(N))` instead — for the 60-odd parts of a real assembly that + /// is 6 levels rather than 59, which is the difference between fitting in + /// the shader's stack and silently dropping geometry. + /// + /// Returns `None` for an empty input: a scene with nothing in it has no + /// root AABB, and every downstream consumer would rather be told than + /// handed a zeroed tree. + pub fn merge_all(mut scenes: Vec) -> Option { + if scenes.is_empty() { + return None; + } + while scenes.len() > 1 { + let mut next = Vec::with_capacity(scenes.len().div_ceil(2)); + let mut it = scenes.into_iter(); + while let Some(a) = it.next() { + match it.next() { + Some(b) => next.push(a.merge(b)), + None => next.push(a), + } + } + scenes = next; + } + scenes.pop() + } + + /// Depth of the packed BVH, in nodes from root to deepest leaf. + /// + /// This is what the shader's traversal stack has to hold. Computed + /// iteratively — a merged assembly tree is deep enough that a recursive + /// walk is a real stack-overflow risk on the host too — and defensive + /// against a malformed tree: a node index that repeats on the current + /// path, or points past the array, terminates that branch rather than + /// looping forever. + pub fn bvh_depth(&self) -> usize { + if self.bvh_nodes.is_empty() { + return 0; + } + let mut best = 0usize; + // (node index, depth). Depth is 1 at the root. + let mut stack = vec![(0u32, 1usize)]; + let mut visited = vec![false; self.bvh_nodes.len()]; + while let Some((idx, depth)) = stack.pop() { + let Some(node) = self.bvh_nodes.get(idx as usize) else { + continue; + }; + if std::mem::replace(&mut visited[idx as usize], true) { + continue; + } + best = best.max(depth); + if node.is_leaf == 0 { + stack.push((node.left_or_first, depth + 1)); + stack.push((node.right_or_count, depth + 1)); + } + } + best + } + + /// Check the packed scene against every limit that would otherwise fail + /// silently or as a driver error, *before* anything is uploaded. + /// + /// `max_binding_bytes` is the device's `max_storage_buffer_binding_size` + /// (`ctx.device.limits()`); pass `None` to skip that check. + /// + /// [`Self::merge`] deliberately does not validate — it is a building + /// block, and checking N times while folding N parts would report the + /// wrong totals. This is the gate to call once, on the finished scene. + pub fn validate(&self, max_binding_bytes: Option) -> Result<(), GpuSceneError> { + let mesh = self.is_mesh_scene(); + if mesh { + if self.surfaces.len() > MAX_MESH_TRIANGLES { + return Err(GpuSceneError::TooManyTriangles(self.surfaces.len())); + } + } else { + if self.surfaces.len() > MAX_SURFACES { + return Err(GpuSceneError::TooManySurfaces(self.surfaces.len())); + } + if self.faces.len() > MAX_FACES { + return Err(GpuSceneError::TooManyFaces(self.faces.len())); + } + } + if self.bvh_nodes.len() > MAX_BVH_NODES { + return Err(GpuSceneError::TooManyBvhNodes(self.bvh_nodes.len())); + } + if self.trim_verts.len() > MAX_TRIM_VERTS { + return Err(GpuSceneError::TooManyTrimVerts(self.trim_verts.len())); + } + let depth = self.bvh_depth(); + if depth > MAX_TRAVERSAL_DEPTH { + return Err(GpuSceneError::BvhTooDeep(depth)); + } + if let Some(cap) = max_binding_bytes { + let bytes = (self.surfaces.len() * std::mem::size_of::()) + .max(self.faces.len() * std::mem::size_of::()) + .max(self.bvh_nodes.len() * std::mem::size_of::()) + as u64; + if bytes > cap { + return Err(GpuSceneError::ExceedsDeviceBinding { bytes, cap }); + } + } + Ok(()) + } + + /// Whether this scene's geometry is triangles rather than trimmed + /// analytic surfaces. Mesh scenes count in the hundreds of thousands and + /// are held to [`MAX_MESH_TRIANGLES`], not to the BRep-scale caps. + fn is_mesh_scene(&self) -> bool { + self.surfaces + .iter() + .all(|s| s.surface_type == SURFACE_TYPE_TRIANGLE) + } + /// Set the material for all faces in the scene. /// /// This replaces the default gray material with the specified color. @@ -1214,3 +1852,229 @@ impl GpuScene { GpuMaterial::from_pbr(crate::pathtrace::Pbr::from_material_def(mat, tint)); } } + +#[cfg(test)] +mod tests { + use super::*; + use vcad_kernel_geom::BilinearSurface; + use vcad_kernel_math::Point3; + use vcad_kernel_primitives::make_cube; + + #[test] + fn analytic_surface_types_are_traceable() { + for t in 0..=4u32 { + let s = GpuSurface { + surface_type: t, + _pad: [0; 3], + params: [0.0; 32], + }; + assert!(s.is_gpu_traceable(), "{} should be traceable", t); + } + } + + #[test] + fn bilinear_and_bspline_are_not_traceable() { + // Both are packed by `from_surface` but hit the WGSL `default` arm, + // which reports a miss — so they must never reach the GPU. + for t in [5u32, 6] { + let s = GpuSurface { + surface_type: t, + _pad: [0; 3], + params: [0.0; 32], + }; + assert!(!s.is_gpu_traceable(), "{} must be rejected", t); + } + } + + #[test] + fn cube_builds_a_gpu_scene() { + let scene = GpuScene::from_brep(&make_cube(10.0, 10.0, 10.0)).expect("cube is analytic"); + assert_eq!(scene.faces.len(), 6); + } + + #[test] + fn unsupported_surface_names_the_offending_geometry() { + // Swap one of the cube's planes for a bilinear patch: the shader has + // no case for it, so building the scene must fail rather than drop + // the face from the image. + let mut cube = make_cube(10.0, 10.0, 10.0); + cube.geometry.surfaces[2] = Box::new(BilinearSurface::new( + Point3::new(0.0, 0.0, 0.0), + Point3::new(1.0, 0.0, 0.0), + Point3::new(0.0, 1.0, 0.0), + Point3::new(1.0, 1.0, 0.0), + )); + + let Err(err) = GpuScene::from_brep(&cube) else { + panic!("bilinear surface must be rejected, not silently dropped"); + }; + match err { + GpuSceneError::UnsupportedSurface { index, name, .. } => { + assert_eq!(index, 2); + assert_eq!(name, "Bilinear"); + } + other => panic!("wrong error: {other:?}"), + } + assert!(err.to_string().contains("Bilinear"), "{err}"); + } + + /// A triangle mesh, tessellated from a cube. + fn cube_mesh() -> TriangleMesh { + vcad_kernel_tessellate::tessellate_brep(&make_cube(10.0, 10.0, 10.0), 16) + } + + #[test] + fn triangles_are_traceable() { + // The counterpart of `bilinear_and_bspline_are_not_traceable`: the + // shader DOES have a case for type 7, and `from_mesh` would reject + // its own output if this said otherwise. + let s = GpuSurface { + surface_type: SURFACE_TYPE_TRIANGLE, + _pad: [0; 3], + params: [0.0; 32], + }; + assert!(s.is_gpu_traceable()); + assert_eq!(GpuSurface::type_name(SURFACE_TYPE_TRIANGLE), "Triangle"); + } + + #[test] + fn triangle_packing_round_trips_positions_and_normals() { + let tri = crate::bvh::FlatTriangle { + positions: [[1.0, 2.0, 3.0], [4.0, 5.0, 6.0], [7.0, 8.0, 9.0]], + normals: Some([[0.0, 0.0, 1.0], [0.0, 1.0, 0.0], [1.0, 0.0, 0.0]]), + }; + let s = GpuSurface::triangle(&tri); + + assert_eq!(s.surface_type, SURFACE_TYPE_TRIANGLE); + // Positions in slots 0..9, normals in 9..18, flag in 18. The WGSL + // reads these by literal index, so the layout is load-bearing. + assert_eq!( + &s.params[0..9], + &[1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0] + ); + assert_eq!( + &s.params[9..18], + &[0.0, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0, 0.0, 0.0] + ); + assert_eq!(s.params[18], 1.0); + // Nothing may spill past the documented block. + assert!(s.params[19..].iter().all(|&v| v == 0.0)); + } + + #[test] + fn triangle_without_normals_clears_the_flag() { + // The shader cannot distinguish an absent normal from a zero one, so + // the flag is the only thing standing between a normal-less mesh and + // normalize(vec3(0)) — i.e. NaN across the whole surface. + let tri = crate::bvh::FlatTriangle { + positions: [[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0]], + normals: None, + }; + let s = GpuSurface::triangle(&tri); + assert_eq!(s.params[18], 0.0); + assert!(s.params[9..18].iter().all(|&v| v == 0.0)); + } + + #[test] + fn mesh_builds_a_gpu_scene_of_triangles() { + let mesh = cube_mesh(); + let scene = GpuScene::from_mesh(&mesh).expect("mesh scene builds"); + + let tris = mesh.indices.len() / 3; + assert_eq!(scene.faces.len(), tris); + assert_eq!(scene.surfaces.len(), tris); + assert!(scene + .surfaces + .iter() + .all(|s| s.surface_type == SURFACE_TYPE_TRIANGLE)); + + // One surface per face, in step. A face pointing at the wrong + // surface would render the wrong triangle in the right BVH slot, + // which is exactly the kind of bug the image tests see only as noise. + assert!(scene + .faces + .iter() + .enumerate() + .all(|(i, f)| f.surface_idx == i as u32)); + + // Triangles carry no trim loops; the shader short-circuits on the + // type code instead. If a nonzero count ever appeared here, every + // mesh hit would be rejected by the winding test. + assert!(scene.faces.iter().all(|f| f.trim_count == 0)); + assert!(scene.faces.iter().all(|f| f.inner_loop_count == 0)); + + // Every BVH leaf must address the face array it was built alongside. + for n in &scene.bvh_nodes { + if n.is_leaf == 1 { + assert!( + (n.left_or_first + n.right_or_count) as usize <= scene.faces.len(), + "leaf range runs past the end of the face array" + ); + } + } + + // WebGPU rejects zero-sized storage buffers, so even the unused + // trim arrays must carry a dummy element. + assert!(!scene.trim_verts.is_empty()); + assert!(!scene.inner_loop_descs.is_empty()); + + // The rig has to be sized to the mesh, or a mesh-only scene renders + // unlit. + assert!(!scene.lights.is_empty()); + } + + #[test] + fn from_mesh_bvh_rejects_a_brep_bvh() { + // A BRep BVH flattens to face IDs, not triangles. Building a mesh + // scene from it would produce an empty one; say so instead. + let bvh = Bvh::build(&make_cube(10.0, 10.0, 10.0)); + let Err(err) = GpuScene::from_mesh_bvh(&bvh) else { + panic!("a BRep-backed BVH is not a mesh"); + }; + assert!(matches!(err, GpuSceneError::NotAMeshBvh)); + assert!(err.to_string().contains("from_brep"), "{err}"); + } + + #[test] + fn merging_a_mesh_into_an_analytic_scene_rebases_every_index() { + let brep = GpuScene::from_brep(&make_cube(10.0, 10.0, 10.0)).expect("analytic half"); + let mesh = GpuScene::from_mesh(&cube_mesh()).expect("mesh half"); + let (brep_faces, brep_surfaces) = (brep.faces.len(), brep.surfaces.len()); + let (mesh_faces, mesh_surfaces) = (mesh.faces.len(), mesh.surfaces.len()); + let mesh_nodes = mesh.bvh_nodes.len(); + + let merged = brep.merge(mesh); + + assert_eq!(merged.faces.len(), brep_faces + mesh_faces); + assert_eq!(merged.surfaces.len(), brep_surfaces + mesh_surfaces); + + // The merged tree gains a new root spanning both halves. + assert_eq!(merged.bvh_nodes[0].is_leaf, 0); + + // No index may dangle after rebasing — this is the whole risk of the + // merge, and a dangling one reads out of bounds in the shader. + for f in &merged.faces { + assert!((f.surface_idx as usize) < merged.surfaces.len()); + assert!((f.material_idx as usize) < merged.materials.len()); + } + for n in &merged.bvh_nodes { + if n.is_leaf == 1 { + assert!((n.left_or_first + n.right_or_count) as usize <= merged.faces.len()); + } else { + assert!((n.left_or_first as usize) < merged.bvh_nodes.len()); + assert!((n.right_or_count as usize) < merged.bvh_nodes.len()); + } + } + + // The mesh half's triangles must still be triangles, and must still + // be reachable from the faces that were rebased onto them. + let mesh_tris = merged.faces[brep_faces..] + .iter() + .filter(|f| { + merged.surfaces[f.surface_idx as usize].surface_type == SURFACE_TYPE_TRIANGLE + }) + .count(); + assert_eq!(mesh_tris, mesh_faces, "merge lost track of the triangles"); + assert!(mesh_nodes > 0); + } +} diff --git a/crates/vcad-kernel-raytrace/src/gpu/mod.rs b/crates/vcad-kernel-raytrace/src/gpu/mod.rs index 1375b1f82..09ac7b838 100644 --- a/crates/vcad-kernel-raytrace/src/gpu/mod.rs +++ b/crates/vcad-kernel-raytrace/src/gpu/mod.rs @@ -8,8 +8,16 @@ mod pipeline; pub mod shaders; pub use buffers::{ - depth_for_frame, GpuAreaLight, GpuBvhNode, GpuCamera, GpuFace, GpuMaterial, GpuRenderState, - GpuScene, GpuSceneError, GpuSurface, GpuVec2, DEFAULT_ENV_INTENSITY, DEFAULT_FIREFLY_CLAMP, - DEFAULT_MAX_DEPTH, DEFAULT_RR_START, + depth_for_frame, halton_jitter, GpuAreaLight, GpuBvhNode, GpuCamera, GpuFace, GpuMaterial, + GpuRenderState, GpuScene, GpuSceneError, GpuSurface, GpuVec2, BACKGROUND_BLACK, + BACKGROUND_ENVIRONMENT, BACKGROUND_SKY, CAMERA_BASIS_DERIVED, CAMERA_BASIS_EXPLICIT, + DEFAULT_ENV_INTENSITY, DEFAULT_FIREFLY_CLAMP, DEFAULT_MAX_DEPTH, DEFAULT_RR_START, + MAX_TRAVERSAL_DEPTH, }; pub use pipeline::RayTracePipeline; + +/// Offline (non-interactive) rendering: upload the scene once, accumulate N +/// samples, read back linear HDR radiance once. Native only — it blocks on +/// the device, which would deadlock the browser event loop. +#[cfg(all(feature = "gpu", not(target_arch = "wasm32")))] +pub use pipeline::{OfflineOptions, OfflineResult}; diff --git a/crates/vcad-kernel-raytrace/src/gpu/pipeline.rs b/crates/vcad-kernel-raytrace/src/gpu/pipeline.rs index bdc1da094..f9df2a10b 100644 --- a/crates/vcad-kernel-raytrace/src/gpu/pipeline.rs +++ b/crates/vcad-kernel-raytrace/src/gpu/pipeline.rs @@ -918,6 +918,580 @@ impl RayTracePipeline { } } +/// Settings for an offline (non-interactive) GPU render. +/// +/// Deliberately narrow compared with [`GpuRenderState`]: the viewport-only +/// knobs (edge overlay, theme, debug modes, adaptive refinement) are all +/// forced off, because an offline render wants the integrator's estimate and +/// nothing painted on top of it. +#[cfg(all(feature = "gpu", not(target_arch = "wasm32")))] +#[derive(Debug, Clone, Copy)] +pub struct OfflineOptions { + /// Output width in pixels. + pub width: u32, + /// Output height in pixels. + pub height: u32, + /// Samples per pixel. Each is one dispatch of the main kernel. + pub spp: u32, + /// Maximum path length. Held constant across all samples — the + /// viewport's `depth_for_frame` escalation exists to make the *first* + /// frame land fast, which an offline render does not care about, and it + /// would bias the running mean toward the shallow early samples. + pub max_depth: u32, + /// Depth at which Russian roulette begins. + pub rr_start: u32, + /// Clamp on indirect radiance to kill fireflies (0 disables). + pub firefly_clamp: f32, + /// Overall multiplier on the studio environment. + pub env_intensity: f32, + /// Whether the implicit ground plane participates in the path trace. + pub ground_enabled: bool, + /// The exact counterpart of `PathTraceOptions::show_background`. + /// + /// `true` puts the *lighting environment* behind the subject; `false` + /// leaves it black, which paired with the film's coverage alpha gives a + /// transparent RGBA render. Either way the viewport's themed `sky_color` + /// backdrop — a UI choice unrelated to the sky the integrator samples — + /// is out of the picture. + /// + /// This has to be a shader-side switch rather than a composite after the + /// fact: a pixel on the silhouette has already averaged its background + /// and surface samples together. + /// + /// Defaults to `true` — an offline render is not a viewport. + pub show_background: bool, + /// RNG decorrelation seed. The same seed and the same scene give the + /// same image, every run, on the same adapter. + pub seed: u32, +} + +#[cfg(all(feature = "gpu", not(target_arch = "wasm32")))] +impl Default for OfflineOptions { + fn default() -> Self { + Self { + width: 512, + height: 512, + spp: 64, + max_depth: super::buffers::DEFAULT_MAX_DEPTH, + rr_start: super::buffers::DEFAULT_RR_START, + firefly_clamp: super::buffers::DEFAULT_FIREFLY_CLAMP, + env_intensity: super::buffers::DEFAULT_ENV_INTENSITY, + ground_enabled: true, + show_background: true, + seed: 0, + } + } +} + +/// The HDR result of an offline render. +/// +/// This is the accumulation buffer verbatim: **linear radiance**, not +/// tonemapped and not gamma-encoded. Exposure, ACES and sRGB encoding belong +/// to the caller — feed [`Self::to_film`] into +/// [`crate::pathtrace::Film::to_srgb8`] to get exactly the CPU renderer's +/// output transform rather than a second, subtly different one. +#[cfg(all(feature = "gpu", not(target_arch = "wasm32")))] +#[derive(Debug, Clone)] +pub struct OfflineResult { + /// Width in pixels. + pub width: u32, + /// Height in pixels. + pub height: u32, + /// Samples per pixel that were actually accumulated. + pub spp: u32, + /// Linear RGB radiance plus coverage, 4 floats per pixel, row-major + /// top-to-bottom. The alpha channel is the shader's sample-count marker + /// (1.0 from the main pass), not opacity. + pub rgba: Vec, +} + +#[cfg(all(feature = "gpu", not(target_arch = "wasm32")))] +impl OfflineResult { + /// Linear RGB at one pixel. + pub fn pixel(&self, x: u32, y: u32) -> [f32; 3] { + let i = ((y * self.width + x) * 4) as usize; + [self.rgba[i], self.rgba[i + 1], self.rgba[i + 2]] + } + + /// Mean relative luminance over the whole image. + pub fn mean_luminance(&self) -> f32 { + if self.rgba.is_empty() { + return 0.0; + } + let sum: f64 = self + .rgba + .as_chunks::<4>() + .0 + .iter() + .map(|p| (0.2126 * p[0] + 0.7152 * p[1] + 0.0722 * p[2]) as f64) + .sum(); + (sum / (self.rgba.len() / 4) as f64) as f32 + } + + /// Repackage as a [`crate::pathtrace::Film`] so the CPU renderer's + /// exposure/ACES/sRGB path can be reused unchanged. + /// + /// The normal, depth, albedo and variance guide buffers are left zeroed — + /// the GPU does not read them back — so the result must not be fed to + /// [`crate::pathtrace::denoise`], which needs them. + pub fn to_film(&self) -> crate::pathtrace::Film { + let n = (self.width * self.height) as usize; + let mut rgb = Vec::with_capacity(n * 3); + let mut alpha = Vec::with_capacity(n); + for p in self.rgba.as_chunks::<4>().0 { + rgb.extend_from_slice(&p[..3]); + alpha.push(p[3]); + } + crate::pathtrace::Film { + width: self.width, + height: self.height, + rgb, + alpha, + normal: vec![0.0; n * 3], + depth: vec![0.0; n], + albedo: vec![0.0; n * 3], + variance: vec![0.0; n], + } + } +} + +#[cfg(all(feature = "gpu", not(target_arch = "wasm32")))] +impl RayTracePipeline { + /// Render `spp` samples per pixel offline and read back the HDR + /// accumulation buffer once. + /// + /// This exists because [`Self::render_with_render_state`] is shaped for a + /// 1-spp progressive viewport: every call recreates every scene buffer, + /// rebuilds the bind group, and reads back the tonemapped `Rgba8Unorm` + /// texture. At 1 frame per user gesture that is free; at 512 spp it is + /// 512 scene uploads and 512 GPU→CPU round trips, and the round trips + /// alone dominate the render. + /// + /// Here the scene is uploaded once, the bind group is built once, and the + /// per-sample loop rewrites only the 128-byte `RenderState` uniform + /// before dispatching the main kernel. The viewport's refine and denoise + /// passes are skipped: both are interactivity aids that trade bias for + /// speed at low sample counts, which is the wrong trade when the whole + /// point is to converge. + /// + /// The readback is the f32 accumulation buffer (binding 8), not the + /// tonemapped texture, so the caller gets linear HDR radiance and can + /// apply exposure and ACES itself. + /// + /// Native only — it blocks on `device.poll(Maintain::Wait)`, which + /// deadlocks the browser's single-threaded event loop. + pub fn render_offline( + &self, + ctx: &GpuContext, + scene: &GpuScene, + camera: &GpuCamera, + opts: &OfflineOptions, + ) -> Result { + use wgpu::util::DeviceExt; + + let (width, height) = (opts.width.max(1), opts.height.max(1)); + let spp = opts.spp.max(1); + + // ── one-time uploads ────────────────────────────────────────────── + let camera_buffer = ctx + .device + .create_buffer_init(&wgpu::util::BufferInitDescriptor { + label: Some("Offline Camera Buffer"), + contents: bytemuck::bytes_of(camera), + usage: wgpu::BufferUsages::UNIFORM, + }); + + // WGSL cannot bind a zero-length storage array, so every empty scene + // buffer still gets one dummy element; the counts in the scene + // structs are what the shader actually loops over. + let storage = |label: &str, bytes: &[u8]| { + ctx.device + .create_buffer_init(&wgpu::util::BufferInitDescriptor { + label: Some(label), + contents: bytes, + usage: wgpu::BufferUsages::STORAGE, + }) + }; + + let surfaces = if scene.surfaces.is_empty() { + vec![super::buffers::GpuSurface::zeroed()] + } else { + scene.surfaces.clone() + }; + let surfaces_buffer = storage("Offline Surfaces", bytemuck::cast_slice(&surfaces)); + + let faces = if scene.faces.is_empty() { + vec![super::buffers::GpuFace::zeroed()] + } else { + scene.faces.clone() + }; + let faces_buffer = storage("Offline Faces", bytemuck::cast_slice(&faces)); + + let bvh_nodes = if scene.bvh_nodes.is_empty() { + vec![super::buffers::GpuBvhNode::zeroed()] + } else { + scene.bvh_nodes.clone() + }; + let bvh_buffer = storage("Offline BVH", bytemuck::cast_slice(&bvh_nodes)); + + let trim_verts = if scene.trim_verts.is_empty() { + vec![super::buffers::GpuVec2 { x: 0.0, y: 0.0 }] + } else { + scene.trim_verts.clone() + }; + let trim_buffer = storage("Offline Trim", bytemuck::cast_slice(&trim_verts)); + + let inner_loop_descs = if scene.inner_loop_descs.is_empty() { + vec![0u32] + } else { + scene.inner_loop_descs.clone() + }; + let inner_loop_descs_buffer = storage( + "Offline Inner Loops", + bytemuck::cast_slice(&inner_loop_descs), + ); + + let materials = if scene.materials.is_empty() { + vec![super::buffers::GpuMaterial::default()] + } else { + scene.materials.clone() + }; + let materials_buffer = storage("Offline Materials", bytemuck::cast_slice(&materials)); + + let lights: Vec = if scene.lights.is_empty() { + vec![super::buffers::GpuAreaLight::default()] + } else { + scene.lights.clone() + }; + let light_buf = storage("Offline Area Lights", bytemuck::cast_slice(&lights)); + + // The output texture is never read back — the HDR accumulation + // buffer is the deliverable — but binding 6 is not optional, and + // every dispatch writes the tonemapped pixel there. + let output_texture = ctx.device.create_texture(&wgpu::TextureDescriptor { + label: Some("Offline Output Texture"), + size: wgpu::Extent3d { + width, + height, + depth_or_array_layers: 1, + }, + mip_level_count: 1, + sample_count: 1, + dimension: wgpu::TextureDimension::D2, + format: wgpu::TextureFormat::Rgba8Unorm, + usage: wgpu::TextureUsages::STORAGE_BINDING, + view_formats: &[], + }); + let output_view = output_texture.create_view(&Default::default()); + + let px_buf_size = (width as u64) * (height as u64) * 16; + let accum = ctx.device.create_buffer(&wgpu::BufferDescriptor { + label: Some("Offline Accumulation Buffer"), + size: px_buf_size, + usage: wgpu::BufferUsages::STORAGE | wgpu::BufferUsages::COPY_SRC, + mapped_at_creation: false, + }); + let depth_normal_buffer = ctx.device.create_buffer(&wgpu::BufferDescriptor { + label: Some("Offline Depth Normal Buffer"), + size: px_buf_size, + usage: wgpu::BufferUsages::STORAGE, + mapped_at_creation: false, + }); + let feature_id_buffer = ctx.device.create_buffer(&wgpu::BufferDescriptor { + label: Some("Offline Feature ID Buffer"), + size: (width as u64) * (height as u64) * 4, + usage: wgpu::BufferUsages::STORAGE, + mapped_at_creation: false, + }); + + let (env_pixels_view, env_cdf_view) = self.offline_env_views(ctx, scene); + + // The one buffer that changes between samples. COPY_DST so the loop + // can `write_buffer` a new frame index instead of allocating. + let mut render_state = self.offline_render_state(scene, opts, 1); + let render_state_buffer = + ctx.device + .create_buffer_init(&wgpu::util::BufferInitDescriptor { + label: Some("Offline Render State Buffer"), + contents: bytemuck::bytes_of(&render_state), + usage: wgpu::BufferUsages::UNIFORM | wgpu::BufferUsages::COPY_DST, + }); + + let bind_group = ctx.device.create_bind_group(&wgpu::BindGroupDescriptor { + label: Some("Offline Ray Trace Bind Group"), + layout: &self.bind_group_layout, + entries: &[ + wgpu::BindGroupEntry { + binding: 0, + resource: camera_buffer.as_entire_binding(), + }, + wgpu::BindGroupEntry { + binding: 1, + resource: surfaces_buffer.as_entire_binding(), + }, + wgpu::BindGroupEntry { + binding: 2, + resource: faces_buffer.as_entire_binding(), + }, + wgpu::BindGroupEntry { + binding: 3, + resource: bvh_buffer.as_entire_binding(), + }, + wgpu::BindGroupEntry { + binding: 4, + resource: trim_buffer.as_entire_binding(), + }, + wgpu::BindGroupEntry { + binding: 5, + resource: inner_loop_descs_buffer.as_entire_binding(), + }, + wgpu::BindGroupEntry { + binding: 6, + resource: wgpu::BindingResource::TextureView(&output_view), + }, + wgpu::BindGroupEntry { + binding: 7, + resource: render_state_buffer.as_entire_binding(), + }, + wgpu::BindGroupEntry { + binding: 8, + resource: accum.as_entire_binding(), + }, + wgpu::BindGroupEntry { + binding: 9, + resource: materials_buffer.as_entire_binding(), + }, + wgpu::BindGroupEntry { + binding: 10, + resource: depth_normal_buffer.as_entire_binding(), + }, + wgpu::BindGroupEntry { + binding: 11, + resource: light_buf.as_entire_binding(), + }, + wgpu::BindGroupEntry { + binding: 12, + resource: feature_id_buffer.as_entire_binding(), + }, + wgpu::BindGroupEntry { + binding: 13, + resource: wgpu::BindingResource::TextureView(&env_pixels_view), + }, + wgpu::BindGroupEntry { + binding: 14, + resource: wgpu::BindingResource::TextureView(&env_cdf_view), + }, + ], + }); + + // ── the sample loop ─────────────────────────────────────────────── + // One submit per sample. `write_buffer` is staged and applied at the + // next submit, so the frame index cannot be updated for several + // dispatches inside a single encoder — they would all read the same + // uniform and collapse the running mean. + let (groups_x, groups_y) = (width.div_ceil(8), height.div_ceil(8)); + for frame_index in 1..=spp { + render_state.frame_index = frame_index; + let (jx, jy) = super::buffers::halton_jitter(frame_index); + render_state.jitter_x = jx; + render_state.jitter_y = jy; + ctx.queue + .write_buffer(&render_state_buffer, 0, bytemuck::bytes_of(&render_state)); + + let mut encoder = ctx + .device + .create_command_encoder(&wgpu::CommandEncoderDescriptor { + label: Some("Offline Ray Trace Encoder"), + }); + { + let mut pass = encoder.begin_compute_pass(&wgpu::ComputePassDescriptor { + label: Some("Offline Ray Trace Pass"), + timestamp_writes: None, + }); + pass.set_pipeline(&self.pipeline); + pass.set_bind_group(0, &bind_group, &[]); + pass.dispatch_workgroups(groups_x, groups_y, 1); + } + ctx.queue.submit(Some(encoder.finish())); + } + + // ── one readback of the HDR accumulation buffer ─────────────────── + let readback = ctx.device.create_buffer(&wgpu::BufferDescriptor { + label: Some("Offline HDR Readback Buffer"), + size: px_buf_size, + usage: wgpu::BufferUsages::COPY_DST | wgpu::BufferUsages::MAP_READ, + mapped_at_creation: false, + }); + let mut encoder = ctx + .device + .create_command_encoder(&wgpu::CommandEncoderDescriptor { + label: Some("Offline Readback Encoder"), + }); + encoder.copy_buffer_to_buffer(&accum, 0, &readback, 0, px_buf_size); + ctx.queue.submit(Some(encoder.finish())); + + let slice = readback.slice(..); + let mapped = { + use std::sync::atomic::{AtomicBool, Ordering}; + use std::sync::Arc; + + let ok = Arc::new(AtomicBool::new(false)); + let ok_clone = ok.clone(); + slice.map_async(wgpu::MapMode::Read, move |r| { + if r.is_ok() { + ok_clone.store(true, Ordering::SeqCst); + } + }); + ctx.device.poll(wgpu::Maintain::Wait); + ok.load(Ordering::SeqCst) + }; + if !mapped { + return Err(GpuError::BufferMapping); + } + + let data = slice.get_mapped_range(); + let rgba: Vec = bytemuck::cast_slice::(&data).to_vec(); + drop(data); + readback.unmap(); + + Ok(OfflineResult { + width, + height, + spp, + rgba, + }) + } + + /// The render state shared by every sample in an offline render. + /// + /// Everything the viewport uses to stay responsive is off: no edge + /// overlay, no stylisation, no debug mode, no adaptive refinement, and a + /// fixed `max_depth` instead of the per-frame escalation. + fn offline_render_state( + &self, + scene: &GpuScene, + opts: &OfflineOptions, + frame_index: u32, + ) -> GpuRenderState { + let mut s = GpuRenderState::new(frame_index); + s.enable_edges = 0; + s.stylize = 0; + s.debug_mode = 0; + s.refine_sample_count = 0; + s.max_depth = opts.max_depth.max(1); + s.rr_start = opts.rr_start; + s.firefly_clamp = opts.firefly_clamp; + s.env_intensity = opts.env_intensity; + s.ground_enabled = u32::from(opts.ground_enabled); + s.seed = opts.seed; + s.background_mode = if opts.show_background { + super::buffers::BACKGROUND_ENVIRONMENT + } else { + super::buffers::BACKGROUND_BLACK + }; + s.light_count = scene.lights.len() as u32; + match &scene.environment { + Some(e) => { + s.env_mode = 1; + s.env_width = e.width; + s.env_height = e.height; + s.env_intensity = e.intensity; + s.env_rotation = e.rotation; + s.env_marg_int = e.marg_int; + } + None => s.env_mode = 0, + } + s + } + + /// Environment texture views for the offline bind group. + /// + /// A gradient-lit scene still needs 1x1 dummies: a texture binding cannot + /// be null, and `env_mode` is what the shader branches on. + fn offline_env_views( + &self, + ctx: &GpuContext, + scene: &GpuScene, + ) -> (wgpu::TextureView, wgpu::TextureView) { + let mk = |label: &str, w: u32, h: u32, fmt: wgpu::TextureFormat, data: &[f32]| { + let tex = ctx.device.create_texture(&wgpu::TextureDescriptor { + label: Some(label), + size: wgpu::Extent3d { + width: w, + height: h, + depth_or_array_layers: 1, + }, + mip_level_count: 1, + sample_count: 1, + dimension: wgpu::TextureDimension::D2, + format: fmt, + usage: wgpu::TextureUsages::TEXTURE_BINDING | wgpu::TextureUsages::COPY_DST, + view_formats: &[], + }); + let bpp = if fmt == wgpu::TextureFormat::Rgba32Float { + 16 + } else { + 4 + }; + ctx.queue.write_texture( + wgpu::ImageCopyTexture { + texture: &tex, + mip_level: 0, + origin: wgpu::Origin3d::ZERO, + aspect: wgpu::TextureAspect::All, + }, + bytemuck::cast_slice(data), + wgpu::ImageDataLayout { + offset: 0, + bytes_per_row: Some(w * bpp), + rows_per_image: Some(h), + }, + wgpu::Extent3d { + width: w, + height: h, + depth_or_array_layers: 1, + }, + ); + tex.create_view(&wgpu::TextureViewDescriptor::default()) + }; + + match &scene.environment { + Some(e) if e.width > 0 && e.height > 0 => ( + mk( + "Offline Environment Pixels", + e.width, + e.height, + wgpu::TextureFormat::Rgba32Float, + &e.pixels, + ), + mk( + "Offline Environment CDF", + e.width + 1, + e.height + 1, + wgpu::TextureFormat::R32Float, + &e.cdf, + ), + ), + _ => ( + mk( + "Offline Environment Pixels (unused)", + 1, + 1, + wgpu::TextureFormat::Rgba32Float, + &[0.0, 0.0, 0.0, 1.0], + ), + mk( + "Offline Environment CDF (unused)", + 1, + 1, + wgpu::TextureFormat::R32Float, + &[0.0], + ), + ), + } + } +} + /// Stub for when GPU feature is not enabled. #[cfg(not(feature = "gpu"))] pub struct RayTracePipeline; diff --git a/crates/vcad-kernel-raytrace/src/gpu/shaders/raytrace.wgsl b/crates/vcad-kernel-raytrace/src/gpu/shaders/raytrace.wgsl index a4e035932..d1420b60c 100644 --- a/crates/vcad-kernel-raytrace/src/gpu/shaders/raytrace.wgsl +++ b/crates/vcad-kernel-raytrace/src/gpu/shaders/raytrace.wgsl @@ -37,10 +37,15 @@ struct Camera { position: vec4, look_at: vec4, up: vec4, + // Screen +x. Only read when basis_mode == 1. + right: vec4, fov: f32, width: u32, height: u32, - _pad: u32, + // 0 = derive the basis right-handedly from `up` (viewport default), + // 1 = use `right`/`up` verbatim, so a MIRRORED CAD view survives the + // trip and `vcad-render --photoreal --gpu` does not flip isometrics. + basis_mode: u32, } struct RenderState { @@ -84,8 +89,13 @@ struct RenderState { env_height: u32, env_rotation: f32, env_marg_int: f32, - _pad3: u32, - _pad4: u32, + // Extra RNG decorrelation term; 0 reproduces the pre-seed noise exactly. + seed: u32, + // What a camera ray that hits nothing returns. + // 0 = `sky_color`, the themed viewport backdrop (default). + // 1 = `env_radiance`, the same sky the integrator lights with — which + // is what the CPU renderer shows behind the subject. + background_mode: u32, _pad5: u32, } @@ -147,6 +157,26 @@ fn pixel_index_i32(coord: vec2) -> u32 { // Utility functions +// The camera's screen basis, as columns (right, up, forward). +// +// In derived mode (basis_mode == 0) `right` is reconstructed right-handedly +// from the up hint, which is what the viewport wants. In explicit mode the +// supplied axes are used as given: a CAD projection basis can be MIRRORED, +// and re-deriving it would flip the render left-for-right against every +// other output style. +fn camera_basis() -> mat3x3 { + let forward = normalize(camera.look_at.xyz - camera.position.xyz); + if camera.basis_mode == 1u { + return mat3x3( + normalize(camera.right.xyz), + normalize(camera.up.xyz), + forward, + ); + } + let r = normalize(cross(forward, camera.up.xyz)); + return mat3x3(r, cross(r, forward), forward); +} + // Core ray generation with an explicit sub-pixel offset. // offset is in pixels, typically in [-0.5, 0.5]. fn ray_origin_and_direction_offset(pixel: vec2, offset: vec2) -> mat2x3 { @@ -160,9 +190,10 @@ fn ray_origin_and_direction_offset(pixel: vec2, offset: vec2) -> mat2x ); // Build camera coordinate system - let forward = normalize(camera.look_at.xyz - camera.position.xyz); - let right = normalize(cross(forward, camera.up.xyz)); - let up = cross(right, forward); + let basis = camera_basis(); + let right = basis[0]; + let up = basis[1]; + let forward = basis[2]; // Compute ray direction let dir = normalize( @@ -613,10 +644,74 @@ fn intersect_torus(origin: vec3, dir: vec3, params: array) -> return hit; } +// Möller-Trumbore ray/triangle test. +// +// Ports `intersect::intersect_triangle`, the CPU tracer's intersector, with +// two deliberate differences forced by f32: +// +// * the degeneracy guard is 1e-8 relative rather than 1e-12 — at f32 +// precision 1e-12 is below the noise floor of the determinant itself, so +// it would admit near-parallel rays whose barycentrics are pure rounding +// error; +// * the barycentric slack is 1e-6 rather than 1e-12, for the same reason. +// Slack matters here: a shared edge must be inclusive from both sides, or +// a mesh grows a lace of single-pixel holes along every triangle border. +// +// Returns barycentrics in `uv`, which `compute_normal` then interpolates the +// vertex normals with. `hit.uv` is (u, v); the third weight is 1 - u - v. +fn intersect_triangle(origin: vec3, dir: vec3, params: array) -> RayHit { + var hit: RayHit; + hit.t = MAX_T; + hit.face_idx = 0xFFFFFFFFu; + + let v0 = vec3(params[0], params[1], params[2]); + let v1 = vec3(params[3], params[4], params[5]); + let v2 = vec3(params[6], params[7], params[8]); + + let e1 = v1 - v0; + let e2 = v2 - v0; + + let pvec = cross(dir, e2); + let det = dot(e1, pvec); + + // Size-relative: an absolute epsilon is meaningless when the model might + // be dimensioned in millimetres or in metres. + let scale = length(e1) * length(e2); + if abs(det) <= 1e-8 * max(scale, 1.0) { + return hit; + } + + let inv_det = 1.0 / det; + let tvec = origin - v0; + + let u = dot(tvec, pvec) * inv_det; + if u < -1e-6 || u > 1.0 + 1e-6 { + return hit; + } + + let qvec = cross(tvec, e1); + let v = dot(dir, qvec) * inv_det; + if v < -1e-6 || u + v > 1.0 + 1e-6 { + return hit; + } + + let t = dot(e2, qvec) * inv_det; + if t <= 0.0 { + return hit; + } + + hit.t = t; + hit.uv = vec2(u, v); + return hit; +} + fn intersect_surface(origin: vec3, dir: vec3, surface_idx: u32) -> RayHit { let surface = surfaces[surface_idx]; switch surface.surface_type { + case SURFACE_TRIANGLE: { + return intersect_triangle(origin, dir, surface.params); + } case SURFACE_PLANE: { return intersect_plane(origin, dir, surface.params); } @@ -715,6 +810,14 @@ fn uv_in_trim_bounds(uv: vec2, start: u32, count: u32) -> bool { fn point_in_face(uv: vec2, face_idx: u32) -> bool { let face = faces[face_idx]; + // A triangle carries no trim loops: `intersect_triangle` already answered + // the containment question that trimming answers for an analytic surface, + // and `uv` here holds barycentrics, not surface parameters. Falling into + // the winding test below would reject every mesh hit (trim_count is 0). + if surfaces[face.surface_idx].surface_type == SURFACE_TRIANGLE { + return true; + } + // Check outer loop - point must be inside if face.trim_count < 3u { // For faces with < 3 trim vertices (e.g., full cylinder walls), @@ -798,8 +901,14 @@ fn trace_bvh(origin: vec3, dir: vec3) -> RayHit { let inv_dir = 1.0 / dir; - // Stack-based traversal - var stack: array; + // Stack-based traversal. + // + // 64 entries, not 32: a merged offline scene (`vcad-render --photoreal + // --gpu` folds one BLAS per solid into a single tree) is deeper than any + // viewport scene, and overflowing here DROPS geometry silently. The host + // validates the packed tree's depth against `MAX_TRAVERSAL_DEPTH` before + // upload, so this bound is checked rather than hoped for. + var stack: array; var stack_ptr = 0; stack[0] = 0u; // Root node stack_ptr = 1; @@ -832,11 +941,11 @@ fn trace_bvh(origin: vec3, dir: vec3) -> RayHit { } } else { // Internal node: push children - if stack_ptr < 31 { + if stack_ptr < 63 { stack[stack_ptr] = node.left_or_first; stack_ptr++; } - if stack_ptr < 31 { + if stack_ptr < 63 { stack[stack_ptr] = node.right_or_count; stack_ptr++; } @@ -988,8 +1097,13 @@ fn in_shadow(p: vec3, light_dir: vec3, max_t: f32) -> bool { // PCG hash → uniform [0, 1) noise. Per-pixel + per-frame seed so the noise // decorrelates across pixels (prevents banding) and animates per frame // (so progressive accumulation averages out). +// +// `render_state.seed` decorrelates whole renders from one another. The +// viewport leaves it at 0, which reproduces the original hash exactly; an +// offline render sets it so a re-run under a different seed is an +// independent — but still reproducible — estimate of the same image. fn rand_uniform(pixel: vec2, sample_idx: u32) -> f32 { - var state = pixel.x * 1973u + pixel.y * 9277u + sample_idx * 26699u + render_state.frame_index * 12345u + 1u; + var state = pixel.x * 1973u + pixel.y * 9277u + sample_idx * 26699u + render_state.frame_index * 12345u + render_state.seed * 2654435761u + 1u; state = state * 747796405u + 2891336453u; let word = ((state >> ((state >> 28u) + 4u)) ^ state) * 277803737u; let r = (word >> 22u) ^ word; @@ -1054,9 +1168,10 @@ fn world_pos_from_depth(pixel: vec2, t: f32) -> vec3 { (f32(pixel.x) + 0.5) / f32(camera.width) * 2.0 - 1.0, 1.0 - (f32(pixel.y) + 0.5) / f32(camera.height) * 2.0 ); - let forward = normalize(camera.look_at.xyz - camera.position.xyz); - let right = normalize(cross(forward, camera.up.xyz)); - let up_cam = cross(right, forward); + let basis = camera_basis(); + let right = basis[0]; + let up_cam = basis[1]; + let forward = basis[2]; let dir = normalize(forward + right * ndc.x * fov_tan * aspect + up_cam * ndc.y * fov_tan); return camera.position.xyz + dir * t; } @@ -1064,9 +1179,10 @@ fn world_pos_from_depth(pixel: vec2, t: f32) -> vec3 { // Project a world-space point onto the screen. Returns pixel coords, or // (-1, -1) when the point is behind the camera or outside the viewport. fn world_to_screen_coords(world_pos: vec3) -> vec2 { - let forward = normalize(camera.look_at.xyz - camera.position.xyz); - let right = normalize(cross(forward, camera.up.xyz)); - let up_cam = cross(right, forward); + let basis = camera_basis(); + let right = basis[0]; + let up_cam = basis[1]; + let forward = basis[2]; let fov_tan = tan(camera.fov * 0.5); let aspect = f32(camera.width) / f32(camera.height); let p = world_pos - camera.position.xyz; @@ -1087,6 +1203,33 @@ fn compute_normal(hit: RayHit) -> vec3 { var normal: vec3; switch surface.surface_type { + case SURFACE_TRIANGLE: { + let v0 = vec3(surface.params[0], surface.params[1], surface.params[2]); + let v1 = vec3(surface.params[3], surface.params[4], surface.params[5]); + let v2 = vec3(surface.params[6], surface.params[7], surface.params[8]); + + // Geometric normal. Non-zero by construction: `Bvh::build_mesh` + // drops zero-area triangles, so this is always a usable fallback. + let geometric = cross(v1 - v0, v2 - v0); + + // Smooth shading: barycentric blend of the vertex normals, so a + // mesh part doesn't read as faceted next to an analytic one. + // Mirrors `MeshGeom::test` on the CPU, including both of its + // fallbacks — no normal array (params[18] == 0), and a blend that + // cancels (opposed vertex normals across a degenerate crease). + normal = geometric; + if surface.params[18] > 0.5 { + let n0 = vec3(surface.params[9], surface.params[10], surface.params[11]); + let n1 = vec3(surface.params[12], surface.params[13], surface.params[14]); + let n2 = vec3(surface.params[15], surface.params[16], surface.params[17]); + let u = hit.uv.x; + let v = hit.uv.y; + let blended = n0 * (1.0 - u - v) + n1 * u + n2 * v; + if length(blended) > 1e-6 { + normal = blended; + } + } + } case SURFACE_PLANE: { normal = vec3(surface.params[9], surface.params[10], surface.params[11]); } @@ -1601,6 +1744,20 @@ fn shade(hit: RayHit, origin: vec3, dir: vec3, pixel: vec2) -> ve // than the lighting environment — the backdrop is a viewport choice, // and `vcad-render` composites its own. The area lights are skipped // here for the same reason `path_trace` skips them at depth 0. + // + // An offline render asks for `env_radiance` instead: `vcad-render + // --photoreal` shows the lighting environment behind the subject, and + // compositing a different backdrop in afterwards is impossible once + // the two have been averaged together inside a partially-covered + // edge pixel. + if render_state.background_mode == 1u { + return vec4(env_radiance(dir), 0.0); + } + // Mode 2 is the CPU's `show_background == false`: leave the backdrop + // black so an RGBA render composites onto any page. + if render_state.background_mode == 2u { + return vec4(0.0, 0.0, 0.0, 0.0); + } return vec4(sky_color(dir), 0.0); } @@ -1892,11 +2049,18 @@ fn main(@builtin(global_invocation_id) global_id: vec3) { // Trace ray using BVH acceleration, then test the implicit ground // plane and pick whichever is closer. var hit = trace_bvh(origin, dir); - let ground = intersect_ground(origin, dir); - if ground.t < hit.t { - hit.t = ground.t; - hit.face_idx = FACE_IDX_GROUND; - hit.uv = vec2(ground.fade, 0.0); + // Gated, like every other `intersect_ground` call. It used not to be, + // which meant `ground_enabled = 0` still drew the implicit floor to + // camera rays — invisible in the viewport, which always enables it, and + // very visible to an offline render that supplies its own floor geometry + // and got a second one underneath it. + if render_state.ground_enabled != 0u { + let ground = intersect_ground(origin, dir); + if ground.t < hit.t { + hit.t = ground.t; + hit.face_idx = FACE_IDX_GROUND; + hit.uv = vec2(ground.fade, 0.0); + } } let new_color = shade(hit, origin, dir, pixel); @@ -2021,9 +2185,16 @@ fn main(@builtin(global_invocation_id) global_id: vec3) { final_color = vec4(heat_color(0.0), 1.0); } - // Store accumulated color with sample count in alpha (1.0 from main pass). - // The refine pass may update this for edge pixels. - accumulated.a = 1.0; + // Alpha carries the refine pass's sample count: 1.0 from the main pass, + // overwritten per edge pixel by `refine` with the number of rays it + // actually fired. Written ONLY when refinement is enabled — with it off + // nothing ever reads the marker, and clobbering alpha throws away the + // integrator's coverage estimate, which is what an offline RGBA render + // turns into its transparency. (`--backdrop none` was coming back fully + // opaque for exactly this reason.) + if render_state.refine_sample_count > 0u { + accumulated.a = 1.0; + } // Store to accumulation buffer and output accum_buffer[pixel_index_i32(pixel_coord)] = accumulated; @@ -2075,11 +2246,13 @@ fn refine(@builtin(global_invocation_id) global_id: vec3) { let dir = ray[1]; var hit = trace_bvh(origin, dir); - let ground = intersect_ground(origin, dir); - if ground.t < hit.t { - hit.t = ground.t; - hit.face_idx = FACE_IDX_GROUND; - hit.uv = vec2(ground.fade, 0.0); + if render_state.ground_enabled != 0u { + let ground = intersect_ground(origin, dir); + if ground.t < hit.t { + hit.t = ground.t; + hit.face_idx = FACE_IDX_GROUND; + hit.uv = vec2(ground.fade, 0.0); + } } color_sum += shade(hit, origin, dir, pixel).rgb; diff --git a/crates/vcad-kernel-raytrace/src/gpu/shaders/surface.wgsl b/crates/vcad-kernel-raytrace/src/gpu/shaders/surface.wgsl index d31ea4865..363ef7aee 100644 --- a/crates/vcad-kernel-raytrace/src/gpu/shaders/surface.wgsl +++ b/crates/vcad-kernel-raytrace/src/gpu/shaders/surface.wgsl @@ -10,6 +10,12 @@ const SURFACE_SPHERE: u32 = 2u; const SURFACE_CONE: u32 = 3u; const SURFACE_TORUS: u32 = 4u; const SURFACE_BILINEAR: u32 = 5u; +// A mesh triangle, packed into a GpuSurface's params by +// `GpuSurface::triangle` in buffers.rs. Not a parametric surface: it has no +// dP/du, so `surface_dpdu` deliberately leaves it in the default arm and the +// shading frame falls back to an arbitrary orthonormal basis — which is what +// an isotropic BSDF wants anyway. +const SURFACE_TRIANGLE: u32 = 7u; const MAX_T: f32 = 1e10; const EPSILON: f32 = 1e-6; diff --git a/crates/vcad-kernel-raytrace/src/intersect/bilinear.rs b/crates/vcad-kernel-raytrace/src/intersect/bilinear.rs index 217ebb971..30fce6325 100644 --- a/crates/vcad-kernel-raytrace/src/intersect/bilinear.rs +++ b/crates/vcad-kernel-raytrace/src/intersect/bilinear.rs @@ -1,6 +1,6 @@ //! Ray-bilinear surface intersection (Newton iteration). -use super::SurfaceHit; +use super::{SurfaceHit, SurfaceHits}; use crate::Ray; use vcad_kernel_geom::BilinearSurface; use vcad_kernel_math::Point2; @@ -14,7 +14,7 @@ const TOLERANCE: f64 = 1e-10; /// /// Uses Newton iteration to find intersections. Returns all valid intersections /// with t >= 0 and (u, v) within [0, 1]. -pub fn intersect_bilinear(ray: &Ray, surface: &BilinearSurface) -> Vec { +pub fn intersect_bilinear(ray: &Ray, surface: &BilinearSurface) -> SurfaceHits { // For planar bilinear patches, use the simpler plane intersection if surface.is_planar() { return intersect_planar_quad(ray, surface); @@ -29,7 +29,7 @@ pub fn intersect_bilinear(ray: &Ray, surface: &BilinearSurface) -> Vec Opti } /// Intersect ray with a planar quad (degenerate bilinear surface). -fn intersect_planar_quad(ray: &Ray, surface: &BilinearSurface) -> Vec { +fn intersect_planar_quad(ray: &Ray, surface: &BilinearSurface) -> SurfaceHits { // Compute plane from first three corners let e1 = surface.p10 - surface.p00; let e2 = surface.p01 - surface.p00; @@ -147,7 +147,7 @@ fn intersect_planar_quad(ray: &Ray, surface: &BilinearSurface) -> Vec Vec Vec Vec { +pub fn intersect_bspline(ray: &Ray, surface: &dyn Surface) -> SurfaceHits { let ((u_min, u_max), (v_min, v_max)) = surface.domain(); // Use subdivision to find initial guesses, then refine with Newton - let mut hits = Vec::new(); + let mut hits = SurfaceHits::new(); subdivide_and_intersect(ray, surface, u_min, u_max, v_min, v_max, 0, &mut hits); // Remove duplicates @@ -40,7 +40,7 @@ fn subdivide_and_intersect( v_min: f64, v_max: f64, depth: usize, - hits: &mut Vec, + hits: &mut SurfaceHits, ) { // Check if the ray might intersect this patch by testing corner bounding box let corners = [ diff --git a/crates/vcad-kernel-raytrace/src/intersect/cone.rs b/crates/vcad-kernel-raytrace/src/intersect/cone.rs index dfc104cce..9574c6d33 100644 --- a/crates/vcad-kernel-raytrace/src/intersect/cone.rs +++ b/crates/vcad-kernel-raytrace/src/intersect/cone.rs @@ -1,6 +1,6 @@ //! Ray-cone intersection (quadratic equation). -use super::SurfaceHit; +use super::{SurfaceHit, SurfaceHits}; use crate::Ray; use std::f64::consts::PI; use vcad_kernel_geom::ConeSurface; @@ -10,7 +10,7 @@ use vcad_kernel_math::Point2; /// /// Returns up to 2 intersections, sorted by t. /// Only intersections with t >= 0 and v >= 0 (on the cone, not the nappes) are returned. -pub fn intersect_cone(ray: &Ray, cone: &ConeSurface) -> Vec { +pub fn intersect_cone(ray: &Ray, cone: &ConeSurface) -> SurfaceHits { let axis = cone.axis.as_ref(); let d = ray.direction.as_ref(); let co = ray.origin - cone.apex; @@ -33,7 +33,7 @@ pub fn intersect_cone(ray: &Ray, cone: &ConeSurface) -> Vec { let b = 2.0 * (d_dot_a * co_dot_a - cos2 * d.dot(co)); let c = co_dot_a * co_dot_a - cos2 * co.dot(co); - let mut hits = Vec::new(); + let mut hits = SurfaceHits::new(); if a.abs() < 1e-12 { // Linear case (ray direction makes exactly the cone half-angle with axis) diff --git a/crates/vcad-kernel-raytrace/src/intersect/cylinder.rs b/crates/vcad-kernel-raytrace/src/intersect/cylinder.rs index 5c42af0f0..730a3928a 100644 --- a/crates/vcad-kernel-raytrace/src/intersect/cylinder.rs +++ b/crates/vcad-kernel-raytrace/src/intersect/cylinder.rs @@ -1,6 +1,6 @@ //! Ray-cylinder intersection (quadratic equation). -use super::SurfaceHit; +use super::{SurfaceHit, SurfaceHits}; use crate::Ray; use vcad_kernel_geom::CylinderSurface; use vcad_kernel_math::Point2; @@ -9,7 +9,7 @@ use vcad_kernel_math::Point2; /// /// Returns up to 2 intersections (entry and exit points), sorted by t. /// Only intersections with t >= 0 are returned. -pub fn intersect_cylinder(ray: &Ray, cylinder: &CylinderSurface) -> Vec { +pub fn intersect_cylinder(ray: &Ray, cylinder: &CylinderSurface) -> SurfaceHits { let axis = cylinder.axis.as_ref(); let d = ray.direction.as_ref(); let oc = ray.origin - cylinder.center; @@ -27,19 +27,19 @@ pub fn intersect_cylinder(ray: &Ray, cylinder: &CylinderSurface) -> Vec; + /// Result of a ray-surface intersection (before trim testing). #[derive(Debug, Clone, Copy)] pub struct SurfaceHit { @@ -73,13 +82,13 @@ pub fn surface_tangent(surface: &dyn Surface, uv: Point2) -> Option { /// Intersect a ray with a surface, returning all intersections sorted by t. /// /// This dispatches to the appropriate intersector based on surface type. -pub fn intersect_surface(ray: &Ray, surface: &dyn Surface) -> Vec { +pub fn intersect_surface(ray: &Ray, surface: &dyn Surface) -> SurfaceHits { match surface.surface_type() { SurfaceKind::Plane => { if let Some(plane) = surface.as_any().downcast_ref::() { intersect_plane(ray, plane).into_iter().collect() } else { - Vec::new() + SurfaceHits::new() } } SurfaceKind::Cylinder => { @@ -89,7 +98,7 @@ pub fn intersect_surface(ray: &Ray, surface: &dyn Surface) -> Vec { { intersect_cylinder(ray, cyl) } else { - Vec::new() + SurfaceHits::new() } } SurfaceKind::Sphere => { @@ -99,7 +108,7 @@ pub fn intersect_surface(ray: &Ray, surface: &dyn Surface) -> Vec { { intersect_sphere(ray, sph) } else { - Vec::new() + SurfaceHits::new() } } SurfaceKind::Cone => { @@ -109,7 +118,7 @@ pub fn intersect_surface(ray: &Ray, surface: &dyn Surface) -> Vec { { intersect_cone(ray, cone) } else { - Vec::new() + SurfaceHits::new() } } SurfaceKind::Torus => { @@ -119,7 +128,7 @@ pub fn intersect_surface(ray: &Ray, surface: &dyn Surface) -> Vec { { intersect_torus(ray, torus) } else { - Vec::new() + SurfaceHits::new() } } SurfaceKind::Bilinear => { @@ -129,7 +138,7 @@ pub fn intersect_surface(ray: &Ray, surface: &dyn Surface) -> Vec { { intersect_bilinear(ray, bil) } else { - Vec::new() + SurfaceHits::new() } } SurfaceKind::BSpline => { diff --git a/crates/vcad-kernel-raytrace/src/intersect/sphere.rs b/crates/vcad-kernel-raytrace/src/intersect/sphere.rs index f8ee6095b..d0f4314ee 100644 --- a/crates/vcad-kernel-raytrace/src/intersect/sphere.rs +++ b/crates/vcad-kernel-raytrace/src/intersect/sphere.rs @@ -1,6 +1,6 @@ //! Ray-sphere intersection (quadratic equation). -use super::SurfaceHit; +use super::{SurfaceHit, SurfaceHits}; use crate::Ray; use std::f64::consts::PI; use vcad_kernel_geom::SphereSurface; @@ -10,7 +10,7 @@ use vcad_kernel_math::Point2; /// /// Returns up to 2 intersections (entry and exit points), sorted by t. /// Only intersections with t >= 0 are returned. -pub fn intersect_sphere(ray: &Ray, sphere: &SphereSurface) -> Vec { +pub fn intersect_sphere(ray: &Ray, sphere: &SphereSurface) -> SurfaceHits { let oc = ray.origin - sphere.center; let d = ray.direction.as_ref(); @@ -21,14 +21,14 @@ pub fn intersect_sphere(ray: &Ray, sphere: &SphereSurface) -> Vec { let discriminant = b * b - 4.0 * a * c; if discriminant < 0.0 { - return Vec::new(); + return SurfaceHits::new(); } let sqrt_disc = discriminant.sqrt(); let t1 = (-b - sqrt_disc) / (2.0 * a); let t2 = (-b + sqrt_disc) / (2.0 * a); - let mut hits = Vec::new(); + let mut hits = SurfaceHits::new(); for t in [t1, t2] { if t < 0.0 { diff --git a/crates/vcad-kernel-raytrace/src/intersect/torus.rs b/crates/vcad-kernel-raytrace/src/intersect/torus.rs index e44e5aa82..d53fd3a20 100644 --- a/crates/vcad-kernel-raytrace/src/intersect/torus.rs +++ b/crates/vcad-kernel-raytrace/src/intersect/torus.rs @@ -2,7 +2,7 @@ //! //! Uses Ferrari's method to solve the quartic polynomial analytically. -use super::SurfaceHit; +use super::{SurfaceHit, SurfaceHits}; use crate::Ray; use std::f64::consts::PI; use vcad_kernel_geom::TorusSurface; @@ -12,7 +12,7 @@ use vcad_kernel_math::Point2; /// /// Returns up to 4 intersections, sorted by t. /// Only intersections with t >= 0 are returned. -pub fn intersect_torus(ray: &Ray, torus: &TorusSurface) -> Vec { +pub fn intersect_torus(ray: &Ray, torus: &TorusSurface) -> SurfaceHits { let r = torus.major_radius; let r2 = r * r; let a = torus.minor_radius; @@ -48,7 +48,7 @@ pub fn intersect_torus(ray: &Ray, torus: &TorusSurface) -> Vec { // Solve the quartic let roots = solve_quartic(c4, c3, c2, c1, c0); - let mut hits: Vec = roots + let mut hits: SurfaceHits = roots .into_iter() .filter(|&t| t >= 0.0) .map(|t| { diff --git a/crates/vcad-kernel-raytrace/src/lib.rs b/crates/vcad-kernel-raytrace/src/lib.rs index 2fec7965b..08e24495e 100644 --- a/crates/vcad-kernel-raytrace/src/lib.rs +++ b/crates/vcad-kernel-raytrace/src/lib.rs @@ -44,7 +44,7 @@ pub mod trim; #[cfg(feature = "gpu")] pub mod gpu; -pub use bvh::Bvh; +pub use bvh::{Bvh, FlatPrims, FlatTriangle}; pub use cpu::{render_scene, render_scene_samples, CpuRenderer}; pub use pathtrace::{ studio_rig, AreaLight, Camera, Environment, Film, Ground, Object, PathTraceOptions, Pbr, Scene, diff --git a/crates/vcad-kernel-raytrace/src/pathtrace.rs b/crates/vcad-kernel-raytrace/src/pathtrace.rs index e3cc8ee27..f3ca04a79 100644 --- a/crates/vcad-kernel-raytrace/src/pathtrace.rs +++ b/crates/vcad-kernel-raytrace/src/pathtrace.rs @@ -923,6 +923,17 @@ pub struct PathTraceOptions { pub show_background: bool, /// Random seed. pub seed: u64, + /// Stop sampling a pixel early once its own variance estimate says the + /// remaining budget cannot move it visibly. + /// + /// [`spp`](Self::spp) becomes a *ceiling* rather than a fixed count. Every + /// pixel still gets at least a floor of samples, and the decision is made + /// from the pixel's own running sums, so the film stays deterministic and + /// independent of how the frame was tiled. + /// + /// Set `false` for a reference render, where a uniform sample count is + /// the point. + pub adaptive: bool, /// Run the edge-aware à-trous denoiser over the film before returning. /// /// This is a pure post-process on the accumulated radiance — it consumes @@ -953,6 +964,7 @@ impl Default for PathTraceOptions { firefly_clamp: Some(12.0), show_background: true, seed: 0x5eed_1234, + adaptive: true, denoise: true, denoise_iters: 5, sigma_normal: 0.35, @@ -997,6 +1009,45 @@ impl Rng { } } +// ─── low-discrepancy sampling ───────────────────────────────────────────── + +/// Van der Corput radical inverse of `i` in `BASE`. +/// +/// Reflects `i`'s digits in `BASE` about the radix point, which spreads +/// consecutive indices as far apart as the base allows. Successive prime +/// bases give the Halton sequence; pairing base 2 with `i / N` gives +/// Hammersley, which is what the camera dimensions use. +#[inline] +fn radical_inverse(mut i: u64) -> f64 { + let inv_base = 1.0 / BASE as f64; + let mut inv_bn = 1.0; + let mut acc = 0u64; + // Accumulate the reversed digits as an integer, then scale once: doing + // the division per digit accumulates rounding error over ~50 digits. + while i > 0 { + let digit = i % BASE as u64; + acc = acc * BASE as u64 + digit; + i /= BASE as u64; + inv_bn *= inv_base; + } + (acc as f64 * inv_bn).min(1.0 - f64::EPSILON) +} + +/// Cranley-Patterson rotation: shift `x` by `offset` on the unit torus. +/// +/// Preserves the point set's discrepancy while randomising its absolute +/// placement, which is what lets every pixel share one low-discrepancy set +/// without the shared structure showing up as a visible pattern. +#[inline] +fn cp_rotate(x: f64, offset: f64) -> f64 { + let v = x + offset; + if v >= 1.0 { + v - 1.0 + } else { + v + } +} + // ─── small math helpers ─────────────────────────────────────────────────── #[inline] @@ -1794,11 +1845,204 @@ pub struct Film { pub variance: Vec, } +/// Side of a square render tile, in pixels. +/// +/// 16 is small enough that a slow corner of the image cannot monopolise a +/// core for long, and large enough that per-tile bookkeeping disappears next +/// to 256 pixels of tracing. +const TILE: usize = 16; + +/// Samples traced between convergence checks. +/// +/// The check needs a sample variance to be worth anything, so it cannot run +/// after every sample; 16 gives a usable estimate and is fine enough that a +/// converged pixel wastes at most 15 samples past the line. +const ADAPTIVE_BATCH: u32 = 16; + +/// Minimum samples every pixel gets, whatever the variance estimate says. +/// +/// A pixel that happens to draw several near-equal samples early reports a +/// tiny variance and would quit while genuinely unconverged — the classic +/// adaptive-sampling failure, and it shows up as blotching in exactly the +/// smooth regions adaptivity was meant to speed up. +const ADAPTIVE_FLOOR: u32 = 32; + +/// Relative tolerance on the 95% confidence half-width of pixel luminance. +/// +/// Tuned against the PSNR gate in `scripts/photoreal-quality.sh`: see the +/// sweep in this change's commit message. +const ADAPTIVE_TOL: f32 = 0.10; + +/// Absolute luminance added to the mean before applying [`ADAPTIVE_TOL`]. +/// +/// Pure relative error never converges in shadow, where the mean approaches +/// zero; pure absolute error over-samples highlights. Adding the two is the +/// usual compromise. +const ADAPTIVE_LUM_FLOOR: f32 = 0.02; + +/// Everything one pixel contributes to the [`Film`]. +/// +/// Traced into a tile-local buffer and blitted into the film afterwards, so +/// the parallel pass needs no aliasing tricks over the shared buffers. +struct PixelResult { + rgb: [f32; 3], + alpha: f32, + normal: [f32; 3], + depth: f32, + albedo: [f32; 3], + variance: f32, +} + +/// One finished tile, in row-major order within its own `w` x `h` extent. +struct TileResult { + x: usize, + y: usize, + w: usize, + h: usize, + pixels: Vec, +} + +/// Trace one pixel's full sample budget. +/// +/// The RNG is seeded from the pixel's coordinates and the option seed alone, +/// never from its position in a work queue — that is what makes the film +/// independent of how [`render`] chose to divide the frame up. +#[allow(clippy::too_many_arguments)] +fn trace_pixel( + scene: &Scene, + accel: &SceneAccel, + cam: &Camera, + opts: &PathTraceOptions, + width: u32, + height: u32, + aspect: f64, + spp: u32, + px: usize, + py: usize, +) -> PixelResult { + let mut rng = + Rng::new(opts.seed ^ ((py as u64) << 32) ^ (px as u64).wrapping_mul(0x9E37_79B9_7F4A_7C15)); + // Cranley-Patterson rotations for the four camera dimensions, drawn once + // per pixel. The low-discrepancy point set below is the *same* for every + // pixel; rotating it by a per-pixel random offset keeps each pixel's + // stratification intact while decorrelating neighbours, so the residual + // error looks like noise rather than a repeating pattern locked to the + // pixel grid. Drawing them from the existing PCG is what keeps --seed + // determinism: no global state, no thread-dependent order. + let rot = [rng.f64(), rng.f64(), rng.f64(), rng.f64()]; + let mut acc = [0.0f32; 3]; + let mut cov = 0.0f32; + // Running sums for the estimator's own variance. + let mut lsum = 0.0f32; + let mut lsum2 = 0.0f32; + let mut out = PixelResult { + rgb: [0.0; 3], + alpha: 0.0, + normal: [0.0; 3], + depth: 0.0, + albedo: [0.0; 3], + variance: 0.0, + }; + + // Sample in batches so the estimator can be asked, between batches, + // whether it has already resolved this pixel. `traced` is the count + // actually spent, which is <= spp under adaptive sampling. + let mut traced = 0u32; + while traced < spp { + let batch = ADAPTIVE_BATCH.min(spp - traced); + for k in 0..batch { + let s = traced + k; + // Pixel jitter and lens position come from a 4D Halton set + // rotated into this pixel's frame, not from four fresh uniforms. + // Four independent uniforms can clump — at 32spp a purely random + // jitter leaves visibly uneven coverage of the pixel footprint, + // and that shows up as extra aliasing on every silhouette. A + // low-discrepancy set covers the square evenly by construction. + // + // Halton rather than Hammersley: Hammersley's first dimension is + // `s / N`, which needs the final sample count up front. Adaptive + // sampling does not know it, and a set that changes shape when + // the loop stops early is worse than a slightly weaker set that + // is correct at every prefix. + let (jx, jy) = ( + cp_rotate(radical_inverse::<2>(s as u64), rot[0]), + cp_rotate(radical_inverse::<3>(s as u64), rot[1]), + ); + let sx = 2.0 * ((px as f64 + jx) / width as f64) - 1.0; + let sy = 1.0 - 2.0 * ((py as f64 + jy) / height as f64); + let (lu, lv) = concentric_disc( + cp_rotate(radical_inverse::<5>(s as u64), rot[2]), + cp_rotate(radical_inverse::<7>(s as u64), rot[3]), + ); + + let ray = cam.ray(sx, sy, aspect, lu, lv); + let (l, primary) = radiance(scene, accel, opts, ray, &mut rng); + acc = add3(acc, l); + let ls = luminance(l); + lsum += ls; + lsum2 += ls * ls; + if primary.hit { + cov += 1.0; + } + if s == 0 { + // Guide buffers come from one primary ray, not an average: + // averaging normals and depths across samples would soften + // exactly the silhouettes the edge-stopping weights exist to + // protect. + out.normal = primary.normal; + out.depth = primary.depth; + out.albedo = primary.albedo; + } + } + traced += batch; + + // Stop once the estimator's own error bar says the remaining samples + // cannot move this pixel by anything a viewer could see. The floor is + // non-negotiable: a pixel that happened to draw several near-equal + // samples early would otherwise report a tiny variance and quit while + // genuinely unconverged. + if opts.adaptive && traced >= ADAPTIVE_FLOOR.min(spp) && traced < spp { + let n = traced as f32; + let mean = lsum / n; + // The clamp is load-bearing, not defensive: once the samples agree + // closely, `lsum2 / n` and `mean * mean` cancel to within f32 + // rounding and can land just below zero, which would put a NaN + // through the sqrt below — and a NaN compares false, so the pixel + // would never converge. + let sample_var = (lsum2 / n - mean * mean).max(0.0) * n / (n - 1.0); + // Half-width of the 95% confidence interval on the mean. + let ci = 1.96 * (sample_var / n).sqrt(); + // Relative tolerance with an absolute floor: pure relative error + // never converges in shadow, where the mean approaches zero, and + // pure absolute error over-samples highlights. + if ci <= ADAPTIVE_TOL * (mean + ADAPTIVE_LUM_FLOOR) { + break; + } + } + } + + let inv = 1.0 / traced as f32; + out.rgb = [acc[0] * inv, acc[1] * inv, acc[2] * inv]; + out.alpha = cov * inv; + // Variance of the *mean*: sample variance / n. A single sample carries no + // information about its own spread, so fall back to the estimate itself + // as a scale. + out.variance = if traced > 1 { + let n = traced as f32; + let mean = lsum * inv; + let sample_var = (lsum2 * inv - mean * mean).max(0.0) * n / (n - 1.0); + sample_var / n + } else { + lsum * lsum + }; + out +} + /// Render `scene` from `cam` into a linear-space [`Film`]. /// -/// Scanlines are traced in parallel. Each pixel's RNG is seeded from its +/// Square tiles are traced in parallel. Each pixel's RNG is seeded from its /// coordinates and the option seed, so output is deterministic and -/// independent of thread scheduling. +/// independent of thread scheduling *and* of the tiling. /// /// When [`PathTraceOptions::denoise`] is set (the default), the film is run /// through [`denoise`] before returning. Pass `denoise: false` for a @@ -1826,79 +2070,59 @@ pub fn render( let mut albedo = vec![0.0f32; (width * height * 3) as usize]; let mut variance = vec![0.0f32; (width * height) as usize]; - let w3 = width as usize * 3; - let w1 = width as usize; - rgb.par_chunks_mut(w3) - .zip(alpha.par_chunks_mut(w1)) - .zip(normal.par_chunks_mut(w3)) - .zip(depth.par_chunks_mut(w1)) - .zip(albedo.par_chunks_mut(w3)) - .zip(variance.par_chunks_mut(w1)) - .enumerate() - .for_each(|(py, (((((row, arow), nrow), drow), brow), vrow))| { - for px in 0..width as usize { - let mut rng = Rng::new( - opts.seed - ^ ((py as u64) << 32) - ^ (px as u64).wrapping_mul(0x9E37_79B9_7F4A_7C15), - ); - let mut acc = [0.0f32; 3]; - let mut cov = 0.0f32; - // Running sums for the estimator's own variance. - let mut lsum = 0.0f32; - let mut lsum2 = 0.0f32; - - for s in 0..spp { - // Jittered pixel position. - let jx = rng.f64(); - let jy = rng.f64(); - let sx = 2.0 * ((px as f64 + jx) / width as f64) - 1.0; - let sy = 1.0 - 2.0 * ((py as f64 + jy) / height as f64); - let (lu, lv) = concentric_disc(rng.f64(), rng.f64()); - - let ray = cam.ray(sx, sy, aspect, lu, lv); - let (l, primary) = radiance(scene, &accel, opts, ray, &mut rng); - acc = add3(acc, l); - let ls = luminance(l); - lsum += ls; - lsum2 += ls * ls; - if primary.hit { - cov += 1.0; - } - if s == 0 { - // Guide buffers come from one primary ray, not an - // average: averaging normals and depths across - // samples would soften exactly the silhouettes the - // edge-stopping weights exist to protect. - nrow[px * 3] = primary.normal[0]; - nrow[px * 3 + 1] = primary.normal[1]; - nrow[px * 3 + 2] = primary.normal[2]; - drow[px] = primary.depth; - brow[px * 3] = primary.albedo[0]; - brow[px * 3 + 1] = primary.albedo[1]; - brow[px * 3 + 2] = primary.albedo[2]; - } + // Traced in tiles rather than scanlines. A scanline is a poor unit of + // work for a path tracer: cost per pixel varies enormously (a pixel that + // escapes to the environment is nearly free, one deep inside a fillet is + // not), and a row spans the whole image, so every row averages to roughly + // the same cost and rayon has nothing left to steal with. Square tiles + // concentrate the cheap and the expensive regions into *different* tasks, + // which is exactly the imbalance work-stealing exists to fix. Tiles are + // also cache-friendlier: a 16x16 neighbourhood of rays hits the same + // corner of the BVH. + // + // Seeding stays per-pixel, so the decomposition is invisible in the + // output: this produces byte-identical films to the scanline version. + let tiles_x = (width as usize).div_ceil(TILE); + let tiles_y = (height as usize).div_ceil(TILE); + let tiles: Vec = (0..tiles_x * tiles_y) + .into_par_iter() + .map(|ti| { + let tx = (ti % tiles_x) * TILE; + let ty = (ti / tiles_x) * TILE; + let tw = TILE.min(width as usize - tx); + let th = TILE.min(height as usize - ty); + let mut pixels = Vec::with_capacity(tw * th); + for py in ty..ty + th { + for px in tx..tx + tw { + pixels.push(trace_pixel( + scene, &accel, cam, opts, width, height, aspect, spp, px, py, + )); } - - let inv = 1.0 / spp as f32; - row[px * 3] = acc[0] * inv; - row[px * 3 + 1] = acc[1] * inv; - row[px * 3 + 2] = acc[2] * inv; - arow[px] = cov * inv; - // Variance of the *mean*: sample variance / spp. A single - // sample carries no information about its own spread, so fall - // back to the estimate itself as a scale. - vrow[px] = if spp > 1 { - let mean = lsum * inv; - let sample_var = - (lsum2 * inv - mean * mean).max(0.0) * spp as f32 / (spp - 1) as f32; - sample_var / spp as f32 - } else { - let mean = lsum; - mean * mean - }; } - }); + TileResult { + x: tx, + y: ty, + w: tw, + h: th, + pixels, + } + }) + .collect(); + + for tile in &tiles { + for ly in 0..tile.h { + for lx in 0..tile.w { + let p = &tile.pixels[ly * tile.w + lx]; + let i = (tile.y + ly) * width as usize + tile.x + lx; + rgb[i * 3..i * 3 + 3].copy_from_slice(&p.rgb); + normal[i * 3..i * 3 + 3].copy_from_slice(&p.normal); + albedo[i * 3..i * 3 + 3].copy_from_slice(&p.albedo); + alpha[i] = p.alpha; + depth[i] = p.depth; + variance[i] = p.variance; + } + } + } let mut film = Film { width, @@ -2373,6 +2597,199 @@ mod tests { assert_eq!(a.rgb, b.rgb, "render must be seed-deterministic"); } + #[test] + fn radical_inverse_matches_hand_computed_values() { + // Base 2: 1 -> 0.1b = 1/2, 2 -> 0.01b = 1/4, 3 -> 0.11b = 3/4. + assert_eq!(radical_inverse::<2>(0), 0.0); + assert!((radical_inverse::<2>(1) - 0.5).abs() < 1e-12); + assert!((radical_inverse::<2>(2) - 0.25).abs() < 1e-12); + assert!((radical_inverse::<2>(3) - 0.75).abs() < 1e-12); + // Base 3: 1 -> 1/3, 2 -> 2/3, 4 = 11_3 -> 0.11_3 = 4/9. + assert!((radical_inverse::<3>(1) - 1.0 / 3.0).abs() < 1e-12); + assert!((radical_inverse::<3>(2) - 2.0 / 3.0).abs() < 1e-12); + assert!((radical_inverse::<3>(4) - 4.0 / 9.0).abs() < 1e-12); + } + + /// The whole point of the point set: no gaps and no clumps. A purely + /// random 2D sample would routinely leave a stratum empty at these + /// counts, which is the aliasing this replaced. + /// + /// Unrotated, `(i/N, phi_2(i))` for `N = 2^m` is a (0, m)-net in base 2: + /// every 8x8 stratum of a 64-point set holds exactly one sample. A + /// Cranley-Patterson rotation shifts the set on the torus and so is no + /// longer a net, but it stays far more even than random — the strongest + /// clump it can produce is two, and no stratum empties by more than that + /// allows. + #[test] + fn camera_point_set_covers_every_stratum() { + let n = 64u32; + let sample = |s: u32, ox: f64, oy: f64| { + ( + cp_rotate((s as f64 + 0.5) / n as f64, ox), + cp_rotate(radical_inverse::<2>(s as u64), oy), + ) + }; + + let mut hits = [[0u32; 8]; 8]; + for s in 0..n { + let (x, y) = sample(s, 0.0, 0.0); + hits[(y * 8.0) as usize][(x * 8.0) as usize] += 1; + } + for (row, counts) in hits.iter().enumerate() { + for (col, &c) in counts.iter().enumerate() { + assert_eq!(c, 1, "unrotated stratum ({col}, {row}) got {c}, want 1"); + } + } + + for &(ox, oy) in &[(0.317, 0.61), (0.94, 0.02)] { + let mut hits = [[0u32; 8]; 8]; + for s in 0..n { + let (x, y) = sample(s, ox, oy); + assert!((0.0..1.0).contains(&x) && (0.0..1.0).contains(&y)); + hits[(y * 8.0) as usize][(x * 8.0) as usize] += 1; + } + let worst = hits.iter().flatten().copied().max().unwrap(); + assert!( + worst <= 2, + "rotated set clumped {worst} samples in a stratum" + ); + } + } + + /// Below the floor, adaptive sampling must be a no-op — not "almost" a + /// no-op. A low-spp render is exactly where an early stop would do the + /// most damage, so the floor has to be a hard gate, not a heuristic. + #[test] + fn adaptive_is_inert_below_the_sample_floor() { + let scene = test_scene(); + let cam = test_camera(); + let base = PathTraceOptions { + spp: ADAPTIVE_FLOOR, + denoise: false, + ..Default::default() + }; + let fixed = render( + &scene, + &cam, + 24, + 24, + &PathTraceOptions { + adaptive: false, + ..base + }, + ); + let adaptive = render( + &scene, + &cam, + 24, + 24, + &PathTraceOptions { + adaptive: true, + ..base + }, + ); + assert_eq!( + fixed.rgb, adaptive.rgb, + "adaptive sampling fired at or below the floor" + ); + } + + /// Adaptive sampling trades samples for time, and the trade is only + /// honest if the picture barely moves. Measured against a converged + /// reference, the adaptive film must land close to the fixed-count film + /// of the same budget — not merely "differ from it". + #[test] + fn adaptive_tracks_the_fixed_count_estimate() { + let scene = test_scene(); + let cam = test_camera(); + let (w, h) = (64, 64); + let reference = render( + &scene, + &cam, + w, + h, + &PathTraceOptions { + spp: 1024, + denoise: false, + adaptive: false, + ..Default::default() + }, + ); + let budget = PathTraceOptions { + spp: 128, + denoise: false, + ..Default::default() + }; + let fixed = render( + &scene, + &cam, + w, + h, + &PathTraceOptions { + adaptive: false, + ..budget + }, + ); + let adaptive = render( + &scene, + &cam, + w, + h, + &PathTraceOptions { + adaptive: true, + ..budget + }, + ); + + let e_fixed = rmse(&fixed, &reference); + let e_adaptive = rmse(&adaptive, &reference); + eprintln!("RMSE vs 1024spp: fixed {e_fixed:.5}, adaptive {e_adaptive:.5}"); + // Adaptive spends fewer samples, so it must be *somewhat* worse; the + // bound is on how much. Measured 1.59x on this scene, which is the + // pessimistic end — it is small, uniformly lit, and gives adaptivity + // almost nothing to skip, where the harness scenes lose only + // ~0.3 dB. 2.0x pins the trade without being brittle. + assert!( + e_adaptive < e_fixed * 2.0, + "adaptive sampling lost too much accuracy: RMSE {e_fixed} -> {e_adaptive}" + ); + } + + /// A frame whose dimensions are not multiples of [`TILE`] must still be + /// covered completely: the edge tiles are clipped, and an off-by-one in + /// the blit would leave an unwritten seam that reads as a black stripe. + #[test] + fn ragged_frame_leaves_no_untraced_seam() { + let scene = test_scene(); + let cam = test_camera(); + // 37x23: both axes straddle the tile grid, neither is a multiple. + let (w, h) = (37u32, 23u32); + let film = render( + &scene, + &cam, + w, + h, + &PathTraceOptions { + spp: 4, + denoise: false, + ..Default::default() + }, + ); + // Every pixel either lands on the subject or looks at the lit + // environment, so none may hold the all-zero value a skipped blit + // would leave behind. + for i in 0..(w * h) as usize { + let lit = + film.rgb[i * 3] != 0.0 || film.rgb[i * 3 + 1] != 0.0 || film.rgb[i * 3 + 2] != 0.0; + assert!( + lit || film.alpha[i] > 0.0, + "pixel ({}, {}) was never written", + i % w as usize, + i / w as usize + ); + } + } + /// The BSDF sampling PDF must match the analytic PDF used by MIS, or /// light sampling and BSDF sampling silently disagree and the image is /// energy-wrong in a way that is hard to see by eye. diff --git a/crates/vcad-kernel-raytrace/src/trim.rs b/crates/vcad-kernel-raytrace/src/trim.rs index 8bb7afe1e..f395312e5 100644 --- a/crates/vcad-kernel-raytrace/src/trim.rs +++ b/crates/vcad-kernel-raytrace/src/trim.rs @@ -9,68 +9,118 @@ use vcad_kernel_math::Point2; use vcad_kernel_primitives::BRepSolid; use vcad_kernel_topo::{FaceId, Orientation}; -/// Test if a UV point is inside a face's trim boundaries. +/// A face's trim boundary, projected into UV once and reusable for every +/// subsequent point test. /// -/// Returns `true` if the point is inside the outer loop and outside all inner loops (holes). -pub fn point_in_face(brep: &BRepSolid, face_id: FaceId, uv: Point2) -> bool { - let topo = &brep.topology; - let face = &topo.faces[face_id]; - let surface = &brep.geometry.surfaces[face.surface_index]; +/// Building this is the expensive half of [`point_in_face`]: every loop +/// vertex is inverse-projected onto the surface (Newton iteration for +/// B-spline and bilinear faces), pole vertices are repaired, and degenerate +/// caps are re-synthesised from the adjacent surface. None of it depends on +/// the query point, so a ray tracer that tests millions of hits against the +/// same face should do it once — see [`FaceTrim::build`] and +/// [`FaceTrim::contains`]. +#[derive(Debug, Clone)] +pub struct FaceTrim { + /// The face spans its whole surface: the outer loop is only a seam. + untrimmed: bool, + /// Outer boundary in UV. Meaningless when `untrimmed`. + outer: Vec, + /// For an untrimmed cylinder or cone, the extent of the loop along the + /// unbounded `v` parameter. + v_range: Option<(f64, f64)>, + /// Hole boundaries in UV. + inners: Vec>, +} - // Get UV coordinates of the outer loop vertices - let raw_uvs = loop_uv_coords(brep, face.outer_loop, surface.as_ref()); - - // A closed surface covering the whole primitive (a full sphere or - // torus) is bounded only by its seam: the outer loop projects to a - // zero-area polygon in UV, which would reject every hit. Treat a - // degenerate outer loop as "untrimmed" — the face spans the entire - // surface — and still honour inner loops (holes) below. - // - // The degeneracy verdict is taken on the *raw* projection, before any - // pole repair: a full sphere's seam loop passes through both poles, and - // repairing those would hand it a non-zero area and reject every hit. - let mut untrimmed = polygon_area(&raw_uvs).abs() < 1e-9; - let mut outer_uvs = if untrimmed { - raw_uvs - } else { - repair_pole_vertices(surface.as_ref(), &raw_uvs) - }; +impl FaceTrim { + /// Project a face's trim loops into UV. + pub fn build(brep: &BRepSolid, face_id: FaceId) -> Self { + let topo = &brep.topology; + let face = &topo.faces[face_id]; + let surface = &brep.geometry.surfaces[face.surface_index]; + + // Get UV coordinates of the outer loop vertices + let raw_uvs = loop_uv_coords(brep, face.outer_loop, surface.as_ref()); + + // A closed surface covering the whole primitive (a full sphere or + // torus) is bounded only by its seam: the outer loop projects to a + // zero-area polygon in UV, which would reject every hit. Treat a + // degenerate outer loop as "untrimmed" — the face spans the entire + // surface — and still honour inner loops (holes) below. + // + // The degeneracy verdict is taken on the *raw* projection, before any + // pole repair: a full sphere's seam loop passes through both poles, + // and repairing those would hand it a non-zero area and reject every + // hit. + let mut untrimmed = polygon_area(&raw_uvs).abs() < 1e-9; + let mut outer = if untrimmed { + raw_uvs + } else { + repair_pole_vertices(surface.as_ref(), &raw_uvs) + }; - // A planar cap bounded by a single closed circle edge (cylinder/cone - // caps) projects to a degenerate UV polygon too — but "untrimmed" on a - // plane means the infinite plane. Rebuild the circle from the adjacent - // surface instead. - if untrimmed { - if let Some(poly) = synthesize_planar_cap_polygon(brep, face_id) { - outer_uvs = poly; - untrimmed = false; + // A planar cap bounded by a single closed circle edge (cylinder/cone + // caps) projects to a degenerate UV polygon too — but "untrimmed" on + // a plane means the infinite plane. Rebuild the circle from the + // adjacent surface instead. + if untrimmed { + if let Some(poly) = synthesize_planar_cap_polygon(brep, face_id) { + outer = poly; + untrimmed = false; + } } - } - if untrimmed { // On a cylinder or cone, v is an unbounded length parameter, so a - // seam-degenerate loop (e.g. a full cylinder wall: only seam - // vertices survive projection, the rim circles collapse) must still - // clamp v to the loop's extent — otherwise the wall traces as an - // infinite cylinder. u legitimately wraps the full turn. - if let Some((v_min, v_max)) = unbounded_v_range(surface.as_ref(), &outer_uvs) { - if uv.y < v_min || uv.y > v_max { - return false; - } + // seam-degenerate loop (e.g. a full cylinder wall: only seam vertices + // survive projection, the rim circles collapse) must still clamp v to + // the loop's extent — otherwise the wall traces as an infinite + // cylinder. u legitimately wraps the full turn. + let v_range = if untrimmed { + unbounded_v_range(surface.as_ref(), &outer) + } else { + None + }; + + let inners = face + .inner_loops + .iter() + .map(|&inner_loop| loop_uv_coords(brep, inner_loop, surface.as_ref())) + .collect(); + + Self { + untrimmed, + outer, + v_range, + inners, } - } else if !point_in_polygon(&wrap_u_into_polygon(uv, &outer_uvs), &outer_uvs) { - return false; } - // Check if point is inside any hole (should be outside all holes) - for &inner_loop in &face.inner_loops { - let inner_uvs = loop_uv_coords(brep, inner_loop, surface.as_ref()); - if point_in_polygon(&uv, &inner_uvs) { - return false; // Inside a hole + /// Is `uv` inside the outer boundary and outside every hole? + pub fn contains(&self, uv: Point2) -> bool { + if self.untrimmed { + if let Some((v_min, v_max)) = self.v_range { + if uv.y < v_min || uv.y > v_max { + return false; + } + } + } else if !point_in_polygon(&wrap_u_into_polygon(uv, &self.outer), &self.outer) { + return false; } + + // Inside a hole is outside the face. + !self.inners.iter().any(|inner| point_in_polygon(&uv, inner)) } +} - true +/// Test if a UV point is inside a face's trim boundaries. +/// +/// Returns `true` if the point is inside the outer loop and outside all inner loops (holes). +/// +/// This reprojects the face's loops on every call. Callers testing the same +/// face repeatedly — the ray tracer, above all — should hold a [`FaceTrim`] +/// instead. +pub fn point_in_face(brep: &BRepSolid, face_id: FaceId, uv: Point2) -> bool { + FaceTrim::build(brep, face_id).contains(uv) } /// Latitude magnitude above which a spherical loop vertex counts as sitting diff --git a/crates/vcad-kernel-raytrace/tests/gpu_mesh.rs b/crates/vcad-kernel-raytrace/tests/gpu_mesh.rs new file mode 100644 index 000000000..a0befe67f --- /dev/null +++ b/crates/vcad-kernel-raytrace/tests/gpu_mesh.rs @@ -0,0 +1,529 @@ +//! Triangle meshes in the wgpu path tracer. +//! +//! `--photoreal` traces cached triangle meshes by default, so until the GPU +//! could intersect a triangle the fast path was unavailable for exactly the +//! geometry the renderer actually feeds it. These tests cover the seam: +//! a mesh-only scene must render, a mixed BRep+mesh scene must render *both* +//! halves, and both must land on top of what the CPU path tracer produces +//! from the same BVH. +//! +//! Parity is measured against the CPU tracer driving the **same** +//! `Bvh::build_mesh` tree, not against the analytic solid the mesh was +//! tessellated from. That isolates what is under test — the WGSL triangle +//! intersector and normal interpolation — from tessellation error, which is a +//! property of the mesh both renderers share. +//! +//! These tests need a real adapter and are `#[ignore]`-tagged like +//! `gpu_smoke.rs` and `gpu_offline.rs`. Run locally with: +//! +//! ```text +//! cargo test -p vcad-kernel-raytrace --features gpu --test gpu_mesh -- --ignored --nocapture +//! ``` + +#![cfg(all(feature = "gpu", not(target_arch = "wasm32")))] + +use std::sync::Arc; + +use vcad_kernel_gpu::{GpuContext, GpuError}; +use vcad_kernel_math::{Point3, Vec3}; +use vcad_kernel_primitives::{make_cube, make_sphere}; +use vcad_kernel_raytrace::gpu::{GpuCamera, GpuScene, OfflineOptions, OfflineResult}; +use vcad_kernel_raytrace::pathtrace::{ + self, Camera, Environment, Object, PathTraceOptions, Pbr, Scene, +}; +use vcad_kernel_raytrace::Bvh; +use vcad_kernel_tessellate::TriangleMesh; + +use vcad_kernel_raytrace::gpu::RayTracePipeline; + +/// Skip with a clear message when no adapter is available. +fn ctx_or_skip(test_name: &str) -> Option<&'static GpuContext> { + match pollster::block_on(GpuContext::init()) { + Ok(ctx) => Some(ctx), + Err(GpuError::NoAdapter) => { + eprintln!("[{test_name}] skipped: no compatible GPU adapter"); + None + } + Err(e) => panic!("GPU init failed unexpectedly: {e}"), + } +} + +/// Subject radius and the camera distance that makes it overfill the frame. +/// Same framing as `gpu_offline.rs`, and for the same reason: every pixel +/// lands on the subject, so a whole-image comparison never straddles the +/// backdrop — where the GPU's themed `sky_color` and the CPU's `env_radiance` +/// legitimately differ. +const R: f64 = 5.0; +const EYE_Z: f64 = 14.0; +const FOV_DEG: f64 = 24.0; + +/// The GPU's default material, spelled as a CPU `Pbr`. `GpuMaterial::default` +/// is 0.7 grey at roughness 0.5, which is *not* `Pbr::default`. +fn gpu_default_material() -> Pbr { + Pbr { + base_color: [0.7, 0.7, 0.7], + metallic: 0.0, + roughness: 0.5, + anisotropy: 0.0, + clearcoat: 0.0, + clearcoat_roughness: 0.1, + ior: 1.5, + emissive: [0.0; 3], + } +} + +/// A tessellated sphere: the mesh-heavy subject the parity tests trace. +fn mesh_sphere(segments: u32) -> TriangleMesh { + vcad_kernel_tessellate::tessellate(&make_sphere(R, segments), segments) +} + +/// The CPU counterpart of the scene the GPU builders produce: the given +/// objects, the studio rig sized to the union of their BVH bounds, the shared +/// gradient environment, no ground. +/// +/// The rig derivation mirrors `studio_lights_for_bvh` in `gpu/buffers.rs`, +/// including how `GpuScene::merge` re-derives it from the merged root — so a +/// mixed scene is lit identically on both sides. +fn cpu_scene(bvhs: Vec>) -> Scene { + let mut min = [f64::MAX; 3]; + let mut max = [f64::MIN; 3]; + for bvh in &bvhs { + let aabb = match bvh.root().expect("BVH has a root") { + vcad_kernel_raytrace::bvh::BvhNode::Leaf { aabb, .. } + | vcad_kernel_raytrace::bvh::BvhNode::Internal { aabb, .. } => *aabb, + }; + let extents = [ + (aabb.min.x, aabb.max.x), + (aabb.min.y, aabb.max.y), + (aabb.min.z, aabb.max.z), + ]; + for (a, (lo, hi)) in extents.into_iter().enumerate() { + min[a] = min[a].min(lo); + max[a] = max[a].max(hi); + } + } + + let center = Point3::new( + (min[0] + max[0]) * 0.5, + (min[1] + max[1]) * 0.5, + (min[2] + max[2]) * 0.5, + ); + let radius = 0.5 + * ((max[0] - min[0]).powi(2) + (max[1] - min[1]).powi(2) + (max[2] - min[2]).powi(2)) + .sqrt(); + + Scene { + objects: bvhs + .into_iter() + .map(|b| Object::new(b, gpu_default_material())) + .collect(), + lights: pathtrace::studio_rig(center, radius), + env: Environment::default(), + ground: None, + } +} + +fn gpu_camera(w: u32, h: u32) -> GpuCamera { + GpuCamera::new( + [0.0, 0.0, EYE_Z as f32], + [0.0, 0.0, 0.0], + [0.0, 1.0, 0.0], + (FOV_DEG as f32).to_radians(), + w, + h, + ) +} + +fn cpu_camera() -> Camera { + Camera::look_at( + Point3::new(0.0, 0.0, EYE_Z), + Point3::new(0.0, 0.0, 0.0), + Vec3::new(0.0, 1.0, 0.0), + FOV_DEG, + ) +} + +fn offline_opts(w: u32, h: u32, spp: u32) -> OfflineOptions { + OfflineOptions { + width: w, + height: h, + spp, + ground_enabled: false, + seed: 7, + ..Default::default() + } +} + +fn cpu_opts(spp: u32) -> PathTraceOptions { + PathTraceOptions { + spp, + adaptive: false, + denoise: false, + show_background: true, + seed: 7, + ..Default::default() + } +} + +/// Peak signal-to-noise ratio between a GPU HDR readback and a CPU render, +/// in dB, over the whole frame. +/// +/// Both images are clamped into [0, 1] first — the peak the ratio is taken +/// against. Radiance above 1 is a specular highlight whose *absolute* error +/// scales with its magnitude, so leaving it unclamped would let a single hot +/// pixel dominate a metric meant to describe the whole image. Clamping is +/// also what a display does, which makes the number mean "how different do +/// these look" rather than "how different are these float buffers". +fn psnr(gpu: &OfflineResult, cpu: &pathtrace::Film) -> f64 { + psnr_masked(gpu, cpu, false) +} + +/// As [`psnr`], but optionally restricted to pixels where the CPU render hit +/// geometry. +/// +/// Masking is required whenever the subject does not fill the frame. The two +/// renderers deliberately draw *different backdrops* — the GPU's `sky_color` +/// is a themed UI choice, the CPU's `env_radiance` is the shared lighting +/// environment — so background pixels compare two things that were never +/// meant to match, and on a mostly-empty frame they swamp the metric. Where +/// the subject fills the frame (`mesh_render_matches_cpu_reference`) the mask +/// is a no-op and the unmasked form is used to keep the comparison honest. +fn psnr_masked(gpu: &OfflineResult, cpu: &pathtrace::Film, subject_only: bool) -> f64 { + let mut mse = 0.0f64; + let mut counted = 0usize; + let n = gpu.rgba.len() / 4; + assert_eq!(n, cpu.rgb.len() / 3, "image sizes differ"); + for i in 0..n { + if subject_only && cpu.depth[i] <= 0.0 { + continue; + } + counted += 1; + for c in 0..3 { + let a = (gpu.rgba[i * 4 + c] as f64).clamp(0.0, 1.0); + let b = (cpu.rgb[i * 3 + c] as f64).clamp(0.0, 1.0); + mse += (a - b) * (a - b); + } + } + assert!(counted > 0, "nothing to compare -- the mask kept no pixels"); + mse /= (counted * 3) as f64; + if mse <= 0.0 { + return f64::INFINITY; + } + 10.0 * (1.0 / mse).log10() +} + +/// The core parity test: a mesh-only GPU render must land on top of the CPU +/// render of the same mesh BVH, at the same spp and seed. +/// +/// The threshold is 25 dB, which is loose on purpose. The two integrators +/// share a BSDF, an environment and a light rig, but not their sampling: +/// different RNG streams, different MIS bookkeeping, f32 against f64. At 64 +/// spp each image still carries visible Monte Carlo noise, and two +/// *independently* noisy estimates of the same signal differ by roughly the +/// sum of their variances — which on this subject sits in the low 30s dB even +/// when both are correct. 25 dB leaves headroom for that while remaining far +/// below what a real defect produces: a dropped shading normal (flat-faceted +/// mesh) lands near 20 dB, and a broken intersector renders background and +/// scores in the single digits. +#[test] +#[ignore = "requires GPU"] +fn mesh_render_matches_cpu_reference() { + let Some(ctx) = ctx_or_skip("mesh_render_matches_cpu_reference") else { + return; + }; + let pipeline = RayTracePipeline::new(ctx).expect("pipeline creation"); + + let mesh = mesh_sphere(48); + let bvh = Arc::new(Bvh::build_mesh(&mesh)); + let scene = GpuScene::from_mesh(&mesh).expect("mesh scene builds"); + + // Every triangle became one surface and one face. + assert_eq!( + scene.surfaces.len(), + scene.faces.len(), + "mesh scenes carry one surface per face" + ); + assert!( + scene.surfaces.iter().all(|s| s.is_gpu_traceable()), + "a packed triangle must be traceable, or every mesh hit is a miss" + ); + + let (w, h) = (64u32, 64u32); + let spp = 64; + let gpu = pipeline + .render_offline(ctx, &scene, &gpu_camera(w, h), &offline_opts(w, h, spp)) + .expect("offline render"); + + assert!( + gpu.rgba.iter().all(|v| v.is_finite()), + "mesh render produced NaN or infinity -- a zero shading normal \ + normalized into garbage is the usual cause" + ); + + let mean = gpu.mean_luminance(); + assert!( + mean > 1e-4, + "mesh render is black (mean luminance {mean}) -- the triangles were \ + uploaded but the shader never hit one" + ); + + let cpu = pathtrace::render(&cpu_scene(vec![bvh]), &cpu_camera(), w, h, &cpu_opts(spp)); + + // Assert the framing premise the whole-image comparison rests on. + let covered = cpu.depth.iter().filter(|&&d| d > 0.0).count(); + let total = (w * h) as usize; + assert!( + covered * 100 >= total * 98, + "the subject does not fill the frame ({covered}/{total} pixels hit)" + ); + + let cpu_mean = cpu + .rgb + .as_chunks::<3>() + .0 + .iter() + .map(|p| (0.2126 * p[0] + 0.7152 * p[1] + 0.0722 * p[2]) as f64) + .sum::() + / total as f64; + let ratio = mean as f64 / cpu_mean; + let db = psnr(&gpu, &cpu); + eprintln!("[mesh parity] GPU mean {mean:.4}, CPU mean {cpu_mean:.4}, ratio {ratio:.3}, PSNR {db:.1} dB"); + + assert!( + (0.5..=2.0).contains(&ratio), + "GPU mean luminance {mean} vs CPU {cpu_mean} (ratio {ratio:.3}) -- \ + the two paths disagree about lighting, not about noise" + ); + assert!( + db >= 25.0, + "GPU/CPU mesh renders differ by more than Monte Carlo noise \ + (PSNR {db:.1} dB, threshold 25)" + ); +} + +/// Shading normals must actually be interpolated. A tessellated sphere with +/// vertex normals reads smooth; the same mesh stripped of them reads +/// faceted. If the shader ignored `params[18]` and the packed normals, the +/// two would render identically. +#[test] +#[ignore = "requires GPU"] +fn interpolated_normals_change_the_shading() { + let Some(ctx) = ctx_or_skip("interpolated_normals_change_the_shading") else { + return; + }; + let pipeline = RayTracePipeline::new(ctx).expect("pipeline creation"); + + // Coarse on purpose: at high tessellation smooth and faceted converge, + // and the difference this test looks for would vanish into the noise. + let smooth = mesh_sphere(12); + let mut faceted = smooth.clone(); + faceted.normals.clear(); + + let (w, h) = (64u32, 64u32); + let cam = gpu_camera(w, h); + let opts = offline_opts(w, h, 32); + + let a = pipeline + .render_offline( + ctx, + &GpuScene::from_mesh(&smooth).expect("smooth scene"), + &cam, + &opts, + ) + .expect("smooth render"); + let b = pipeline + .render_offline( + ctx, + &GpuScene::from_mesh(&faceted).expect("faceted scene"), + &cam, + &opts, + ) + .expect("faceted render"); + + // Same geometry, same seed, same sample count: any difference is the + // shading normal and nothing else. + let diff: f64 = a + .rgba + .iter() + .zip(&b.rgba) + .map(|(x, y)| (x - y).abs() as f64) + .sum::() + / a.rgba.len() as f64; + eprintln!("[normals] mean |smooth - faceted| = {diff:.5}"); + assert!( + diff > 1e-3, + "dropping the vertex normals changed nothing (mean diff {diff:.6}) -- \ + the shader is not reading the packed shading normals" + ); + + // ...and the faceted one must still be a valid image, not NaN soup from + // normalizing the zeroed normal slots. + assert!(b.rgba.iter().all(|v| v.is_finite())); + assert!(b.mean_luminance() > 1e-4); +} + +/// A merged BRep + mesh scene must render both halves. +/// +/// `GpuScene::merge` rebases every cross-buffer index; a triangle's +/// `surface_idx` has to survive that just as an analytic face's does. The two +/// subjects are placed side by side so each owns a known half of the frame, +/// which lets the test check coverage per region rather than trusting a +/// whole-image statistic to notice one of them missing. +#[test] +#[ignore = "requires GPU"] +fn merged_brep_and_mesh_scene_shows_both() { + let Some(ctx) = ctx_or_skip("merged_brep_and_mesh_scene_shows_both") else { + return; + }; + let pipeline = RayTracePipeline::new(ctx).expect("pipeline creation"); + + // Analytic sphere at the origin, mesh cube well clear of it to the +x + // side. The cube is displaced by editing its mesh vertices rather than + // by transforming the solid: this test is about the GPU consuming + // triangles, and a mesh translation is exact, so it introduces nothing + // the parity check would then have to tolerate. + let sphere = make_sphere(2.0, 32); + let mut cube_mesh = vcad_kernel_tessellate::tessellate(&make_cube(3.0, 3.0, 3.0), 16); + // `make_cube` spans 0..3 on each axis; centre it on y and z, then push + // it out to x in 5..8. + for v in cube_mesh.vertices.as_chunks_mut::<3>().0 { + v[0] += 5.0; + v[1] -= 1.5; + v[2] -= 1.5; + } + + let scene = GpuScene::from_brep(&sphere) + .expect("analytic half builds") + .merge(GpuScene::from_mesh(&cube_mesh).expect("mesh half builds")); + + // The merge must have kept every surface traceable and every face + // pointing at a surface that exists. + assert!(scene.surfaces.iter().all(|s| s.is_gpu_traceable())); + assert!( + scene + .faces + .iter() + .all(|f| (f.surface_idx as usize) < scene.surfaces.len()), + "merge left a face pointing past the end of the surface array" + ); + + let (w, h) = (96u32, 96u32); + let cam = GpuCamera::new( + [0.0, 0.0, 20.0], + [0.0, 0.0, 0.0], + [0.0, 1.0, 0.0], + 45f32.to_radians(), + w, + h, + ); + let out = pipeline + .render_offline(ctx, &scene, &cam, &offline_opts(w, h, 32)) + .expect("merged render"); + + assert!(out.rgba.iter().all(|v| v.is_finite())); + + // Depth is the honest coverage signal: the background has its own + // radiance, so a bright pixel does not prove a subject is there. The + // offline result carries no depth, so use the CPU render of the same + // scene for the geometry check and the GPU render for the parity one. + let cpu = pathtrace::render( + &cpu_scene(vec![ + Arc::new(Bvh::build(&sphere)), + Arc::new(Bvh::build_mesh(&cube_mesh)), + ]), + &Camera::look_at( + Point3::new(0.0, 0.0, 20.0), + Point3::new(0.0, 0.0, 0.0), + Vec3::new(0.0, 1.0, 0.0), + 45.0, + ), + w, + h, + &cpu_opts(32), + ); + + // Column ranges each subject projects into. At z=20 with a 45deg vertical + // FOV the frame half-width at z=0 is 20*tan(22.5deg) = 8.28, so world x + // maps to column 48 + 48*x/8.28: the sphere (x in -2..2) lands in 36..60 + // and the cube (x in 5..8) in 77..94. The bands below are those, trimmed + // in a little so a pixel of antialiasing at the silhouette cannot leak + // one subject's coverage into the other's count. + let hits_in = |x0: u32, x1: u32| -> usize { + (0..h) + .flat_map(|y| (x0..x1).map(move |x| (x, y))) + .filter(|(x, y)| cpu.depth[(y * w + x) as usize] > 0.0) + .count() + }; + let (left, right) = (hits_in(38, 58), hits_in(79, 92)); + assert!(left > 100, "the analytic sphere is missing ({left} px)"); + assert!(right > 100, "the mesh cube is missing ({right} px)"); + + // Both halves must also agree between the renderers. This is the check + // that would catch a merge rebasing triangles onto the wrong surfaces: + // the image would still be non-empty on both sides, but wrong. + let db = psnr_masked(&out, &cpu, true); + eprintln!("[merged] left {left} px, right {right} px, subject PSNR {db:.1} dB"); + assert!( + db >= 20.0, + "merged BRep+mesh render disagrees with the CPU reference \ + (PSNR {db:.1} dB) -- index rebasing across the merge is suspect" + ); +} + +/// Measurement, not an assertion: mesh-heavy render, GPU against CPU. +/// Run with `--ignored --nocapture` to see the numbers. +#[test] +#[ignore = "benchmark; requires GPU"] +fn bench_mesh_gpu_vs_cpu() { + let Some(ctx) = ctx_or_skip("bench_mesh_gpu_vs_cpu") else { + return; + }; + let pipeline = RayTracePipeline::new(ctx).expect("pipeline creation"); + + let mesh = mesh_sphere(300); + let tris = mesh.indices.len() / 3; + let (w, h) = (512u32, 512u32); + let spp = 128u32; + + let t_build = std::time::Instant::now(); + let scene = GpuScene::from_mesh(&mesh).expect("mesh scene builds"); + let build = t_build.elapsed(); + + let surface_mb = + scene.surfaces.len() * std::mem::size_of::(); + eprintln!( + "[bench] {tris} triangles -> {} surfaces / {} faces / {} BVH nodes, \ + {:.1} MB of surface buffer, scene build {:?}", + scene.surfaces.len(), + scene.faces.len(), + scene.bvh_nodes.len(), + surface_mb as f64 / (1024.0 * 1024.0), + build, + ); + + // Warm up shader/pipeline caches so the timed run is not paying for them. + let cam = gpu_camera(w, h); + let _ = pipeline + .render_offline(ctx, &scene, &cam, &offline_opts(w, h, 2)) + .expect("warmup"); + + let t0 = std::time::Instant::now(); + let out = pipeline + .render_offline(ctx, &scene, &cam, &offline_opts(w, h, spp)) + .expect("gpu render"); + let gpu = t0.elapsed(); + assert!(out.mean_luminance() > 0.0); + + let bvh = Arc::new(Bvh::build_mesh(&mesh)); + let cpu_sc = cpu_scene(vec![bvh]); + let t1 = std::time::Instant::now(); + let _ = pathtrace::render(&cpu_sc, &cpu_camera(), w, h, &cpu_opts(spp)); + let cpu = t1.elapsed(); + + eprintln!( + "[bench] {w}x{h} @ {spp} spp, {tris} tris: GPU {:?}, CPU {:?} ({:.1}x)", + gpu, + cpu, + cpu.as_secs_f64() / gpu.as_secs_f64(), + ); +} diff --git a/crates/vcad-kernel-raytrace/tests/gpu_offline.rs b/crates/vcad-kernel-raytrace/tests/gpu_offline.rs new file mode 100644 index 000000000..15f945168 --- /dev/null +++ b/crates/vcad-kernel-raytrace/tests/gpu_offline.rs @@ -0,0 +1,382 @@ +//! Offline GPU render: persistent buffers, N-spp accumulation, one HDR +//! readback. +//! +//! `RayTracePipeline::render_offline` is the entry point `vcad-render` will +//! eventually sit on. The viewport's `render_with_render_state` cannot serve +//! that role: it rebuilds every scene buffer and reads back a tonemapped +//! `Rgba8Unorm` texture on *every* call, which at 512 spp means 512 scene +//! uploads and 512 GPU->CPU round trips. +//! +//! These tests need a real adapter and are `#[ignore]`-tagged like +//! `gpu_smoke.rs`. Run locally with: +//! +//! ```text +//! cargo test -p vcad-kernel-raytrace --features gpu --test gpu_offline -- --ignored --nocapture +//! ``` + +#![cfg(all(feature = "gpu", not(target_arch = "wasm32")))] + +use std::sync::Arc; + +use vcad_kernel_gpu::{GpuContext, GpuError}; +use vcad_kernel_math::{Point3, Vec3}; +use vcad_kernel_primitives::make_sphere; +use vcad_kernel_raytrace::gpu::{GpuCamera, GpuScene, OfflineOptions, RayTracePipeline}; +use vcad_kernel_raytrace::pathtrace::{ + self, Camera, Environment, Object, PathTraceOptions, Pbr, Scene, +}; +use vcad_kernel_raytrace::Bvh; + +/// Skip with a clear message when no adapter is available. +fn ctx_or_skip(test_name: &str) -> Option<&'static GpuContext> { + match pollster::block_on(GpuContext::init()) { + Ok(ctx) => Some(ctx), + Err(GpuError::NoAdapter) => { + eprintln!("[{test_name}] skipped: no compatible GPU adapter"); + None + } + Err(e) => panic!("GPU init failed unexpectedly: {e}"), + } +} + +/// Sphere radius, and the camera distance that makes it overfill the frame. +const R: f64 = 5.0; +const EYE_Z: f64 = 14.0; +/// Vertical FOV. The sphere's silhouette half-angle is asin(5/14) = 20.9deg; +/// a 24deg vertical FOV puts the frame *corner* at 16.7deg, comfortably +/// inside it. Every pixel therefore hits the sphere, which is what lets the +/// CPU/GPU comparison below use a whole-image mean: the two renderers agree +/// on lighting but deliberately differ on the *visible backdrop* (the GPU's +/// `sky_color` is a themed UI choice, `env_radiance` is the shared one), so a +/// frame with visible background would compare two different things. +const FOV_DEG: f64 = 24.0; + +/// The GPU's default material, spelled as a CPU `Pbr`. `GpuMaterial::default` +/// is 0.7 grey at roughness 0.5, which is *not* `Pbr::default`. +fn gpu_default_material() -> Pbr { + Pbr { + base_color: [0.7, 0.7, 0.7], + metallic: 0.0, + roughness: 0.5, + anisotropy: 0.0, + clearcoat: 0.0, + clearcoat_roughness: 0.1, + ior: 1.5, + emissive: [0.0; 3], + } +} + +/// The CPU counterpart of the scene `GpuScene::from_brep` builds: same solid, +/// same material, same studio rig derived from the same BVH root bounds, same +/// gradient environment, no ground. +fn cpu_scene(bvh: Arc) -> Scene { + let root = bvh.root().expect("sphere BVH has a root"); + let aabb = match root { + vcad_kernel_raytrace::bvh::BvhNode::Leaf { aabb, .. } + | vcad_kernel_raytrace::bvh::BvhNode::Internal { aabb, .. } => *aabb, + }; + let center = Point3::new( + (aabb.min.x + aabb.max.x) * 0.5, + (aabb.min.y + aabb.max.y) * 0.5, + (aabb.min.z + aabb.max.z) * 0.5, + ); + let radius = 0.5 + * ((aabb.max.x - aabb.min.x).powi(2) + + (aabb.max.y - aabb.min.y).powi(2) + + (aabb.max.z - aabb.min.z).powi(2)) + .sqrt(); + + Scene { + objects: vec![Object::new(bvh, gpu_default_material())], + lights: pathtrace::studio_rig(center, radius), + env: Environment::default(), + ground: None, + } +} + +fn gpu_camera(w: u32, h: u32) -> GpuCamera { + GpuCamera::new( + [0.0, 0.0, EYE_Z as f32], + [0.0, 0.0, 0.0], + [0.0, 1.0, 0.0], + (FOV_DEG as f32).to_radians(), + w, + h, + ) +} + +fn cpu_camera() -> Camera { + Camera::look_at( + Point3::new(0.0, 0.0, EYE_Z), + Point3::new(0.0, 0.0, 0.0), + Vec3::new(0.0, 1.0, 0.0), + FOV_DEG, + ) +} + +fn offline_opts(w: u32, h: u32, spp: u32) -> OfflineOptions { + OfflineOptions { + width: w, + height: h, + spp, + ground_enabled: false, + seed: 7, + ..Default::default() + } +} + +/// The core test: a 64-spp offline render must return finite, non-negative, +/// non-constant HDR radiance whose mean luminance lands near an equivalent +/// CPU render. +#[test] +#[ignore = "requires GPU"] +fn offline_hdr_matches_cpu_mean_luminance() { + let Some(ctx) = ctx_or_skip("offline_hdr_matches_cpu_mean_luminance") else { + return; + }; + let pipeline = RayTracePipeline::new(ctx).expect("pipeline creation"); + + let sphere = make_sphere(R, 32); + let scene = GpuScene::from_brep(&sphere).expect("scene builds"); + + let (w, h) = (64u32, 64u32); + let spp = 64; + let out = pipeline + .render_offline(ctx, &scene, &gpu_camera(w, h), &offline_opts(w, h, spp)) + .expect("offline render"); + + assert_eq!(out.width, w); + assert_eq!(out.height, h); + assert_eq!(out.spp, spp); + assert_eq!( + out.rgba.len() as u32, + w * h * 4, + "HDR buffer is the wrong size" + ); + + // The readback is raw f32 from a compute shader: NaN or a negative + // radiance means the integrator is broken, and no tonemap will save it. + assert!( + out.rgba.iter().all(|v| v.is_finite()), + "HDR buffer contains NaN or infinity" + ); + assert!( + out.rgba + .as_chunks::<4>() + .0 + .iter() + .all(|p| p[..3].iter().all(|&c| c >= 0.0)), + "HDR buffer contains negative radiance" + ); + + let mean = out.mean_luminance(); + assert!( + mean > 1e-4, + "offline render is black (mean luminance {mean}) -- the dispatch \ + loop ran but wrote nothing to the accumulation buffer" + ); + + // A constant image would mean the accumulation buffer was never actually + // shaded (e.g. every dispatch no-op'd and we read back a cleared buffer). + let lum: Vec = out + .rgba + .as_chunks::<4>() + .0 + .iter() + .map(|p| 0.2126 * p[0] + 0.7152 * p[1] + 0.0722 * p[2]) + .collect(); + let (lo, hi) = lum + .iter() + .fold((f32::MAX, f32::MIN), |(lo, hi), &v| (lo.min(v), hi.max(v))); + assert!( + hi - lo > 1e-3, + "offline render is a constant image (luminance {lo}..{hi})" + ); + + // CPU reference over the same scene, camera and sample count. + let bvh = Arc::new(Bvh::build(&sphere)); + let cpu = pathtrace::render( + &cpu_scene(bvh), + &cpu_camera(), + w, + h, + &PathTraceOptions { + spp, + adaptive: false, + denoise: false, + show_background: true, + seed: 7, + ..Default::default() + }, + ); + + // Every pixel must be on the sphere, or the two renderers are being + // compared over different backdrops and the mean is meaningless. + let covered = cpu.depth.iter().filter(|&&d| d > 0.0).count(); + let total = (w * h) as usize; + assert!( + covered * 100 >= total * 98, + "the subject does not fill the frame ({covered}/{total} pixels hit) -- \ + the CPU/GPU mean-luminance comparison assumes it does" + ); + + let cpu_mean = cpu + .rgb + .as_chunks::<3>() + .0 + .iter() + .map(|p| (0.2126 * p[0] + 0.7152 * p[1] + 0.0722 * p[2]) as f64) + .sum::() + / total as f64; + let ratio = mean as f64 / cpu_mean; + + // Deliberately loose. The two integrators share a BSDF and an + // environment but not their sampling: different RNG, different MIS + // bookkeeping, f32 on the GPU against f64 on the CPU. A factor-of-two + // band still catches every failure that matters -- a dead light rig, a + // dropped bounce, a tonemap sneaking into the accumulation buffer -- and + // does not fire on Monte Carlo disagreement. + assert!( + (0.5..=2.0).contains(&ratio), + "GPU mean luminance {mean} vs CPU {cpu_mean} (ratio {ratio:.3}) -- \ + these should agree to within a factor of two; a large gap means the \ + two paths disagree about lighting, not about noise" + ); +} + +/// The same seed must give the same image, twice. Without this the offline +/// path is unusable for regression baselines. +#[test] +#[ignore = "requires GPU"] +fn offline_render_is_deterministic_for_a_fixed_seed() { + let Some(ctx) = ctx_or_skip("offline_render_is_deterministic_for_a_fixed_seed") else { + return; + }; + let pipeline = RayTracePipeline::new(ctx).expect("pipeline creation"); + let sphere = make_sphere(R, 32); + let scene = GpuScene::from_brep(&sphere).expect("scene builds"); + + let (w, h) = (32u32, 32u32); + let opts = offline_opts(w, h, 8); + let cam = gpu_camera(w, h); + + let a = pipeline + .render_offline(ctx, &scene, &cam, &opts) + .expect("render a"); + let b = pipeline + .render_offline(ctx, &scene, &cam, &opts) + .expect("render b"); + assert_eq!(a.rgba, b.rgba, "same seed produced a different image"); + + // ...and a different seed must actually change the noise, or the seed + // plumbing into RenderState is not reaching the shader's RNG. + let c = pipeline + .render_offline(ctx, &scene, &cam, &OfflineOptions { seed: 99, ..opts }) + .expect("render c"); + assert_ne!( + a.rgba, c.rgba, + "changing the seed did not change the image -- render_state.seed is \ + not reaching rand_uniform" + ); +} + +/// More samples must reduce the estimator's noise. This is the cheapest +/// end-to-end proof that the accumulation loop is actually *accumulating* +/// rather than overwriting the buffer each dispatch. +#[test] +#[ignore = "requires GPU"] +fn more_samples_reduce_noise() { + let Some(ctx) = ctx_or_skip("more_samples_reduce_noise") else { + return; + }; + let pipeline = RayTracePipeline::new(ctx).expect("pipeline creation"); + let sphere = make_sphere(R, 32); + let scene = GpuScene::from_brep(&sphere).expect("scene builds"); + + let (w, h) = (48u32, 48u32); + let cam = gpu_camera(w, h); + + // Local neighbour-difference energy: a noisy image has large pixel-to- + // pixel jumps on a smooth sphere, a converged one does not. + let roughness = |img: &vcad_kernel_raytrace::gpu::OfflineResult| -> f64 { + let mut acc = 0.0f64; + for y in 0..h { + for x in 1..w { + let a = img.pixel(x, y); + let b = img.pixel(x - 1, y); + acc += (0..3).map(|c| (a[c] - b[c]).abs() as f64).sum::(); + } + } + acc / ((w - 1) * h) as f64 + }; + + let low = pipeline + .render_offline(ctx, &scene, &cam, &offline_opts(w, h, 1)) + .expect("1 spp"); + let high = pipeline + .render_offline(ctx, &scene, &cam, &offline_opts(w, h, 64)) + .expect("64 spp"); + + let (r_low, r_high) = (roughness(&low), roughness(&high)); + assert!( + r_high < r_low * 0.8, + "64 spp is no smoother than 1 spp ({r_high:.5} vs {r_low:.5}) -- the \ + sample loop is overwriting the accumulation buffer instead of \ + averaging into it" + ); +} + +/// Measurement, not an assertion: per-spp cost of the persistent-buffer loop +/// against calling the viewport entry point in a loop, same scene and spp. +/// Run with `--ignored --nocapture` to see the numbers. +#[test] +#[ignore = "benchmark; requires GPU"] +fn bench_offline_vs_viewport_loop() { + let Some(ctx) = ctx_or_skip("bench_offline_vs_viewport_loop") else { + return; + }; + let pipeline = RayTracePipeline::new(ctx).expect("pipeline creation"); + let sphere = make_sphere(R, 32); + let scene = GpuScene::from_brep(&sphere).expect("scene builds"); + + let (w, h) = (512u32, 512u32); + let spp = 128u32; + let cam = gpu_camera(w, h); + + // Warm up shader/pipeline caches so the first timed run is not paying + // for them. + let _ = pipeline + .render_offline(ctx, &scene, &cam, &offline_opts(w, h, 2)) + .expect("warmup"); + + let t0 = std::time::Instant::now(); + let out = pipeline + .render_offline(ctx, &scene, &cam, &offline_opts(w, h, spp)) + .expect("offline render"); + let offline = t0.elapsed(); + assert!(out.mean_luminance() > 0.0); + + let t1 = std::time::Instant::now(); + let mut accum = None; + for frame in 1..=spp { + let state = vcad_kernel_raytrace::gpu::GpuRenderState::new(frame); + let (_pixels, buf) = pollster::block_on( + pipeline.render_with_render_state(ctx, &scene, &cam, w, h, accum, state), + ) + .expect("viewport render"); + accum = Some(buf); + } + let viewport = t1.elapsed(); + + eprintln!( + "\n{w}x{h} @ {spp} spp\n \ + render_offline: {:>9.3?} ({:>8.3?} / spp)\n \ + render_with_render_state x N: {:>9.3?} ({:>8.3?} / spp)\n \ + speedup: {:.2}x\n", + offline, + offline / spp, + viewport, + viewport / spp, + viewport.as_secs_f64() / offline.as_secs_f64(), + ); +} diff --git a/crates/vcad-render/Cargo.toml b/crates/vcad-render/Cargo.toml index 6c320f740..27e024581 100644 --- a/crates/vcad-render/Cargo.toml +++ b/crates/vcad-render/Cargo.toml @@ -15,7 +15,20 @@ raster = ["dep:image"] cli = ["dep:clap", "dep:vcad-loon"] # Direct BRep ray tracing (`--raytrace`): pixel-perfect raster output via # vcad-kernel-raytrace. Off for the WASM build so it doesn't grow. -raytrace = ["raster", "dep:vcad-kernel-raytrace"] +# vcad-kernel-raytrace already pulls rayon in for the tracer itself; the +# photoreal scene build uses it directly to fan out per-solid BVH construction. +raytrace = ["raster", "dep:vcad-kernel-raytrace", "dep:rayon"] +# `--photoreal --gpu`: path-trace the same scene on a wgpu compute pipeline +# instead of on rayon. Off by default — it drags in wgpu and its whole backend +# stack, which is heavy for a CLI most people run on the CPU. A build without +# it still *has* the flag, and answers with "not compiled in" rather than +# "unrecognised argument", so the diagnostic points at the build. +photoreal-gpu = [ + "raytrace", + "vcad-kernel-raytrace/gpu", + "dep:vcad-kernel-gpu", + "dep:pollster", +] [dependencies] vcad-eval = { workspace = true } @@ -27,6 +40,9 @@ clap = { version = "4", features = ["derive"], optional = true } image = { version = "0.25", default-features = false, features = ["jpeg", "png", "hdr"], optional = true } vcad-kernel-raytrace = { workspace = true, optional = true } vcad-loon = { workspace = true, optional = true } +rayon = { version = "1", optional = true } +vcad-kernel-gpu = { workspace = true, optional = true } +pollster = { workspace = true, optional = true } [lib] name = "vcad_render" @@ -42,3 +58,7 @@ required-features = ["cli"] # need it unconditionally. vcad-loon = { workspace = true } vcad-ecad-symbols = { workspace = true } +# Integration tests decode the PNGs the renderer emits. `image` is an +# *optional* normal dependency (behind `raster`), which an integration test +# cannot reach through the library, so it is named again here. +image = { version = "0.25", default-features = false, features = ["jpeg", "png"] } diff --git a/crates/vcad-render/examples/psnr.rs b/crates/vcad-render/examples/psnr.rs new file mode 100644 index 000000000..889663fa2 --- /dev/null +++ b/crates/vcad-render/examples/psnr.rs @@ -0,0 +1,225 @@ +//! Image-quality harness for photoreal render changes. +//! +//! Compares a candidate PNG against a reference PNG and reports PSNR (dB) and +//! grayscale SSIM. Used to gate sampling changes in the path tracer: a +//! candidate render at the normal sample count must stay above a PSNR floor +//! relative to a high-spp reference of the same scene. +//! +//! ```text +//! cargo run --release -p vcad-render --example psnr -- ref.png test.png [--min-psnr 35] +//! ``` +//! +//! Exits non-zero when `--min-psnr` is given and the candidate falls below it, +//! so it can drive a shell gate directly. +//! +//! Both metrics run on the images' 8-bit sRGB samples, which is what a viewer +//! actually sees; PSNR uses all three channels, SSIM the Rec. 709 luma. + +use std::process::ExitCode; + +fn main() -> ExitCode { + let args: Vec = std::env::args().skip(1).collect(); + let mut paths: Vec<&str> = Vec::new(); + let mut min_psnr: Option = None; + let mut i = 0; + while i < args.len() { + match args[i].as_str() { + "--min-psnr" => { + i += 1; + match args.get(i).and_then(|s| s.parse::().ok()) { + Some(v) => min_psnr = Some(v), + None => { + eprintln!("--min-psnr expects a number"); + return ExitCode::from(2); + } + } + } + "-h" | "--help" => { + println!("usage: psnr [--min-psnr ]"); + return ExitCode::SUCCESS; + } + other => paths.push(other), + } + i += 1; + } + if paths.len() != 2 { + eprintln!("usage: psnr [--min-psnr ]"); + return ExitCode::from(2); + } + + let reference = match load_rgb(paths[0]) { + Ok(v) => v, + Err(e) => { + eprintln!("{}: {e}", paths[0]); + return ExitCode::from(2); + } + }; + let candidate = match load_rgb(paths[1]) { + Ok(v) => v, + Err(e) => { + eprintln!("{}: {e}", paths[1]); + return ExitCode::from(2); + } + }; + if reference.width != candidate.width || reference.height != candidate.height { + eprintln!( + "size mismatch: {}x{} vs {}x{}", + reference.width, reference.height, candidate.width, candidate.height + ); + return ExitCode::from(2); + } + + let psnr = psnr(&reference, &candidate); + let ssim = ssim(&reference, &candidate); + let psnr_str = if psnr.is_infinite() { + "inf (identical)".to_string() + } else { + format!("{psnr:.2} dB") + }; + println!("psnr {psnr_str}"); + println!("ssim {ssim:.5}"); + + match min_psnr { + Some(floor) if psnr < floor => { + eprintln!("FAIL: psnr {psnr:.2} dB below floor {floor:.2} dB"); + ExitCode::FAILURE + } + _ => ExitCode::SUCCESS, + } +} + +/// An 8-bit RGB image, three bytes per pixel. +struct Image { + width: u32, + height: u32, + rgb: Vec, +} + +fn load_rgb(path: &str) -> Result { + let img = image::open(path).map_err(|e| e.to_string())?.to_rgb8(); + Ok(Image { + width: img.width(), + height: img.height(), + rgb: img.into_raw(), + }) +} + +/// Peak signal-to-noise ratio over all three channels, 255 peak. +/// +/// Infinite for byte-identical images. +fn psnr(a: &Image, b: &Image) -> f64 { + let mut sse = 0.0f64; + for (x, y) in a.rgb.iter().zip(b.rgb.iter()) { + let d = *x as f64 - *y as f64; + sse += d * d; + } + let mse = sse / a.rgb.len() as f64; + if mse == 0.0 { + return f64::INFINITY; + } + 10.0 * (255.0 * 255.0 / mse).log10() +} + +/// Rec. 709 luma, 0..255. +fn luma(rgb: &[u8], i: usize) -> f64 { + 0.2126 * rgb[i * 3] as f64 + 0.7152 * rgb[i * 3 + 1] as f64 + 0.0722 * rgb[i * 3 + 2] as f64 +} + +/// Mean grayscale SSIM over 8x8 windows (stride 4). +/// +/// The standard 11x11 Gaussian window is overkill for a gate; a uniform +/// window over a dense stride tracks the same structural differences and is +/// far simpler to read. +fn ssim(a: &Image, b: &Image) -> f64 { + const WIN: usize = 8; + const STRIDE: usize = 4; + // Stabilising constants from Wang et al. 2004, for L = 255. + let c1 = (0.01 * 255.0f64).powi(2); + let c2 = (0.03 * 255.0f64).powi(2); + + let w = a.width as usize; + let h = a.height as usize; + if w < WIN || h < WIN { + return f64::NAN; + } + + let mut total = 0.0f64; + let mut windows = 0usize; + let mut y0 = 0; + while y0 + WIN <= h { + let mut x0 = 0; + while x0 + WIN <= w { + let (mut sa, mut sb, mut saa, mut sbb, mut sab) = (0.0, 0.0, 0.0, 0.0, 0.0); + for dy in 0..WIN { + for dx in 0..WIN { + let i = (y0 + dy) * w + x0 + dx; + let va = luma(&a.rgb, i); + let vb = luma(&b.rgb, i); + sa += va; + sb += vb; + saa += va * va; + sbb += vb * vb; + sab += va * vb; + } + } + let n = (WIN * WIN) as f64; + let ma = sa / n; + let mb = sb / n; + // Unbiased (n-1) variance/covariance, as in the reference paper. + let va = (saa - sa * ma) / (n - 1.0); + let vb = (sbb - sb * mb) / (n - 1.0); + let cov = (sab - sa * mb) / (n - 1.0); + let num = (2.0 * ma * mb + c1) * (2.0 * cov + c2); + let den = (ma * ma + mb * mb + c1) * (va + vb + c2); + total += num / den; + windows += 1; + x0 += STRIDE; + } + y0 += STRIDE; + } + total / windows as f64 +} + +#[cfg(test)] +mod tests { + use super::*; + + fn solid(width: u32, height: u32, v: u8) -> Image { + Image { + width, + height, + rgb: vec![v; (width * height * 3) as usize], + } + } + + #[test] + fn identical_images_score_perfectly() { + let a = solid(16, 16, 128); + let b = solid(16, 16, 128); + assert!(psnr(&a, &b).is_infinite()); + assert!((ssim(&a, &b) - 1.0).abs() < 1e-9); + } + + #[test] + fn known_offset_matches_closed_form() { + // Every sample off by 1 => MSE 1 => PSNR = 20*log10(255). + let a = solid(16, 16, 100); + let b = solid(16, 16, 101); + let expected = 20.0 * 255.0f64.log10(); + assert!((psnr(&a, &b) - expected).abs() < 1e-9); + } + + #[test] + fn structural_difference_drops_ssim() { + let a = solid(32, 32, 10); + let mut b = solid(32, 32, 10); + for i in 0..(32 * 32) { + if (i / 32) % 2 == 0 { + b.rgb[i * 3] = 240; + b.rgb[i * 3 + 1] = 240; + b.rgb[i * 3 + 2] = 240; + } + } + assert!(ssim(&a, &b) < 0.1); + } +} diff --git a/crates/vcad-render/src/animate.rs b/crates/vcad-render/src/animate.rs index c877574de..e596a4f1a 100644 --- a/crates/vcad-render/src/animate.rs +++ b/crates/vcad-render/src/animate.rs @@ -237,7 +237,7 @@ pub fn render_photoreal_animation( // Static scene roots first (their transforms never change), then the // assembly instances, whose object→world transform is rewritten per // frame. - let mut objects = build_objects(&ev.statics).unwrap_or_default(); + let mut objects = build_objects(&ev.statics, pr.mesh_segments()).unwrap_or_default(); let static_count = objects.len(); let part_solids: Vec = std::mem::take(&mut ev.parts) @@ -245,7 +245,7 @@ pub fn render_photoreal_animation( .map(crate::part_as_local_scene_solid) .collect(); let instance_ids: Vec = part_solids.iter().map(|s| s.id.clone()).collect(); - objects.extend(build_objects(&part_solids)?); + objects.extend(build_objects(&part_solids, pr.mesh_segments())?); let articulated = objects.len() - static_count; // build_objects fails closed on untraceable parts, so counts can only // agree here; keep the invariant visible without implying a live branch. diff --git a/crates/vcad-render/src/lib.rs b/crates/vcad-render/src/lib.rs index fe05a49e1..a673e4553 100644 --- a/crates/vcad-render/src/lib.rs +++ b/crates/vcad-render/src/lib.rs @@ -47,6 +47,8 @@ pub mod materials; pub mod pcb; #[cfg(feature = "raytrace")] pub mod photoreal; +#[cfg(feature = "photoreal-gpu")] +pub mod photoreal_gpu; /// First PCB in a raw `.vcad` document, if any: `PcbBoard` nodes are /// checked in node-id order, then the legacy top-level `pcb` field. @@ -620,9 +622,20 @@ pub fn tessellation_segments(raster_size_px: Option) -> u32 { /// `segments` (see [`tessellation_segments`]). /// /// A cache hit yields a mesh-backed `Solid` (no BRep), which the tessellated -/// raster and SVG paths render identically to a freshly evaluated root. The -/// ray-traced, photoreal and `--section` paths need analytic surfaces, so -/// callers must not wrap those in a cache scope. +/// raster and SVG paths render identically to a freshly evaluated root. +/// +/// The photoreal path is cacheable too, and by default is: it traces +/// triangles at `photoreal::MESH_SEGMENTS`, so a cached mesh is all the +/// geometry it wants, and it tessellates on a miss as well so cold and warm +/// renders agree. Wrap it in a scope built with those segments — +/// [`tessellation_segments`] answers for the raster and SVG paths, not for +/// this one. +/// +/// `--raytrace`, `--section` and photoreal's `--exact` +/// (`PhotorealOptions::exact`) do need analytic BRep surfaces, +/// which a cached mesh cannot supply. Callers must not wrap those in a cache +/// scope: they would pay to populate a cache they can never read, and a hit +/// would silently downgrade them to triangles. pub fn with_root_cache( cache: std::rc::Rc, segments: u32, @@ -3008,7 +3021,9 @@ mod raster { /// Segment count for canvases at or above [`HIRES_THRESHOLD_PX`]. /// Adjacent facets differ by 2.8°, still well under the ~10° coplanar /// tolerance, so no facet stripes appear. - const RASTER_SEGMENTS_HIRES: u32 = 128; + /// Also the photoreal path's fixed count — see + /// `photoreal::MESH_SEGMENTS`. + pub(crate) const RASTER_SEGMENTS_HIRES: u32 = 128; /// The segment count a raster canvas of `size_px` tessellates at. pub(super) fn segments_for(size_px: u32) -> u32 { diff --git a/crates/vcad-render/src/main.rs b/crates/vcad-render/src/main.rs index efb306180..3312c6e66 100644 --- a/crates/vcad-render/src/main.rs +++ b/crates/vcad-render/src/main.rs @@ -291,9 +291,9 @@ struct Cli { /// `~/.cache/vcad`; `VCAD_CACHE=0` has the same effect). The cache is /// keyed on each root's resolved expression plus the kernel build, so /// it never serves geometry from a different kernel or an edited root; - /// this flag is for timing the kernel itself. `--raytrace`, - /// `--photoreal` and `--section` need BRep surfaces and bypass the - /// cache regardless. + /// this flag is for timing the kernel itself. `--raytrace`, `--section` + /// and `--photoreal --exact` need analytic BRep surfaces, which a cached + /// mesh cannot supply, and bypass the cache regardless. #[arg(long)] no_cache: bool, @@ -304,6 +304,51 @@ struct Cli { #[arg(long, conflicts_with = "raytrace")] photoreal: bool, + /// Trace analytic BRep surface intersection (`--photoreal`): slower, + /// sharper curved silhouettes. + /// + /// By default the path tracer runs against a 128-segment tessellation, + /// which the on-disk root-mesh cache can serve — so a re-render skips + /// re-evaluating the whole document and pays only for the tracing. A + /// BRep cannot be cached (it has no serialized form), so `--exact` + /// re-evaluates the kernel on every run — and, on a document of any size, + /// takes far longer than the tracing does. Reach for it on hero renders + /// containing small-radius cylinders or fillets, where the tessellation's + /// facets show along the silhouette. Note it also reframes very slightly + /// (the camera fits bounds that a tessellation inscribes), so `--exact` + /// and default output never overlay pixel-for-pixel. + #[arg(long, requires = "photoreal")] + exact: bool, + + /// Path-trace on the GPU instead of the CPU (`--photoreal`). + /// + /// Same scene, same studio rig, same exposure and tonemap — a wgpu + /// compute pipeline integrates it instead of rayon. + /// + /// HONOURED: --spp, --max-depth, --exposure, --fov, --seed, --size, + /// --fill, --auto-aspect, --view/--azimuth/--elevation, --env, + /// --env-rotation, per-part materials, and --backdrop studio|none. + /// + /// REFUSED, with a message rather than a wrong image: --exact, + /// --aperture, --ortho, --backdrop shadow-catcher, --animate. + /// + /// IGNORED: --no-adaptive (the GPU has no adaptive sampler, so --spp is + /// an exact count rather than a ceiling), and denoising — it needs + /// normal/depth/albedo guide buffers the GPU tracer does not read back, + /// so --gpu turns it off and says so on stderr. + /// + /// Two knowing differences from the CPU image: the studio floor is a + /// large quad rather than an infinite plane, and the integrator runs in + /// f32 rather than f64. The GPU render converges, but to a slightly + /// different picture — around 30 dB PSNR from the CPU render on rose-pro, + /// and no closer at 1024 spp than at 256. Use it for fast looks and + /// sweeps; render the final hero on the CPU. + /// + /// An explicit --gpu with no usable adapter is an error, never a silent + /// fall back to the CPU. + #[arg(long, requires = "photoreal")] + gpu: bool, + /// Samples per pixel for `--photoreal`. 32 for a quick look, 512+ for a /// clean hero render. #[arg(long, default_value_t = 128, requires = "photoreal")] @@ -359,6 +404,17 @@ struct Cli { #[arg(long, requires = "photoreal")] no_denoise: bool, + /// Sample every pixel to the full `--spp` count (`--photoreal`), instead + /// of stopping pixels whose own variance estimate says the remaining + /// samples cannot move them visibly. + /// + /// Adaptive sampling is on by default and typically cuts render time + /// substantially at unchanged quality; `--spp` is a ceiling under it. + /// Turn it off for reference renders, or to compare two images at a + /// genuinely equal sample count. + #[arg(long, requires = "photoreal")] + no_adaptive: bool, + /// Render a jointed assembly over time (`--photoreal`): one PNG per /// timeline sample into the directory given by `-o`, evaluating the /// document's geometry exactly once for the whole sequence. @@ -507,9 +563,65 @@ fn photoreal_options(cli: &Cli) -> vcad_render::photoreal::PhotorealOptions { }, seed: cli.seed, denoise: !cli.no_denoise, + adaptive: !cli.no_adaptive, + exact: cli.exact, } } +/// `--photoreal --gpu`: the same scene, path-traced on a wgpu compute +/// pipeline. +/// +/// Denoising is turned off here rather than in `photoreal_options`, so the +/// note fires exactly once and only when the user actually asked for both +/// `--gpu` and a denoised render. The denoiser is guided by the per-pixel +/// normal/depth/albedo buffers the CPU integrator records, and the GPU reads +/// back only radiance; run against the zeroed guides it would not produce a +/// worse image, it would silently produce *no* filtering at all, which is the +/// version of this the user must not be left guessing about. +#[cfg(all(feature = "raytrace", feature = "photoreal-gpu"))] +fn render_photoreal_gpu( + raw: &str, + opts: &vcad_render::RasterOptions, + pr: &vcad_render::photoreal::PhotorealOptions, + png: bool, +) -> Result, String> { + use vcad_render::photoreal_gpu; + + // Refusals first, so a run that is about to be rejected outright does not + // print a note about a denoiser it will never reach. + photoreal_gpu::check_supported(opts, pr)?; + + let mut pr = pr.clone(); + if pr.denoise { + eprintln!("{}", photoreal_gpu::denoise_is_unavailable()); + pr.denoise = false; + } + if png { + photoreal_gpu::render_photoreal_gpu_png_str(raw, opts, &pr) + } else { + photoreal_gpu::render_photoreal_gpu_jpeg_str(raw, opts, &pr) + } +} + +/// `--gpu` in a build compiled without the `photoreal-gpu` feature. +/// +/// Deliberately an error naming the feature rather than a silent CPU render: +/// a user who asked for the GPU and got the CPU would take every timing +/// afterwards as a GPU number. +#[cfg(all(feature = "raytrace", not(feature = "photoreal-gpu")))] +fn render_photoreal_gpu( + _raw: &str, + _opts: &vcad_render::RasterOptions, + _pr: &vcad_render::photoreal::PhotorealOptions, + _png: bool, +) -> Result, String> { + Err("--gpu: this build of vcad-render was compiled without the \ + `photoreal-gpu` feature. Rebuild with \ + `cargo build -p vcad-render --features photoreal-gpu`, or drop --gpu \ + to path-trace on the CPU." + .to_string()) +} + /// Render an animation: one PNG per timeline sample into the `-o` directory, /// then (unless told not to) mux them into an mp4. /// @@ -522,6 +634,21 @@ fn run_animation(input: &Path, cli: &Cli) -> Result<(), String> { assemble_mp4, parse_timeline_spec, render_photoreal_animation, AnimateOptions, }; + // `--gpu` bakes each part's transform into the uploaded vertices — the + // WGSL tracer has no instancing layer to put a per-object matrix in — so + // every frame of a sequence would mean repacking and re-uploading the + // whole scene. That is precisely the cost `--animate` exists to pay once, + // so the two do not compose and saying so beats quietly re-uploading. + if cli.gpu { + return Err( + "--animate does not compose with --gpu: the GPU scene bakes \ + each part's pose into its uploaded vertices, so every frame \ + would repack and re-upload the whole assembly — which is the \ + one cost --animate exists to avoid. Drop --gpu." + .to_string(), + ); + } + let spec = cli.animate.as_ref().expect("caller checked --animate"); let Some(out_dir) = cli.output.clone() else { return Err("--animate writes a directory of frames: pass -o

".to_string()); @@ -601,8 +728,9 @@ fn render_raster(raw: &str, cli: &Cli, format: Format) -> Result, String let png = format == Format::Png; let opts = raster_opts(cli, png); if cli.photoreal { - // Same constraint as --raytrace: the path tracer needs analytic BRep - // surfaces, and the overlays are drawn by the projected 2D path. + // The overlays are drawn by the projected 2D path, which has no + // counterpart here; --section additionally leaves mesh-backed solids + // that --exact could not trace analytically anyway. if cli.section.is_some() || cli.annotations().any() { return Err( "--photoreal does not compose with --section/--axes/--labels/--dims; \ @@ -613,6 +741,9 @@ fn render_raster(raw: &str, cli: &Cli, format: Format) -> Result, String #[cfg(feature = "raytrace")] { let pr = photoreal_options(cli); + if cli.gpu { + return render_photoreal_gpu(raw, &opts, &pr, png); + } return if png { vcad_render::photoreal::render_photoreal_png_str(raw, &opts, &pr) } else { @@ -796,13 +927,35 @@ fn is_loon(path: &Path) -> bool { /// The root-mesh cache for this invocation, if any: off by flag or /// environment, and off for the paths that need analytic BRep surfaces /// (which a cached mesh can't supply). +/// +/// Plain `--photoreal` is *not* one of those paths — it traces a +/// tessellation, and re-evaluating every BRep to produce one on each run is +/// the single largest fixed cost of a photoreal render. Only `--exact` opts +/// back into analytic surfaces, and therefore out of the cache. fn root_cache(cli: &Cli) -> Option> { - if cli.no_cache || cli.raytrace || cli.photoreal || cli.section.is_some() { + if cli.no_cache || cli.raytrace || cli.section.is_some() { + return None; + } + if cli.photoreal && cli.exact { return None; } vcad_eval::cache::DiskMeshCache::from_env().map(std::rc::Rc::new) } +/// The segment count the cached mesh must be tessellated at for this +/// invocation to be able to use it — the facet count the render would +/// otherwise produce itself. +/// +/// The photoreal path fixes its own count regardless of canvas size (see +/// `photoreal::MESH_SEGMENTS`); every other path scales with the output. +fn cache_segments(cli: &Cli, format: Format) -> u32 { + #[cfg(feature = "raytrace")] + if cli.photoreal { + return vcad_render::photoreal::MESH_SEGMENTS; + } + vcad_render::tessellation_segments(raster_size_px(cli, format)) +} + /// The raster canvas size this render will use, or `None` for SVG. #[cfg(feature = "raster")] fn raster_size_px(cli: &Cli, format: Format) -> Option { @@ -822,7 +975,7 @@ fn render_one(input: &Path, dest: Option<&Path>, format: Format, cli: &Cli) -> R Some(cache) => { // The cached mesh must carry the facet count this output would // tessellate at, or a hit renders coarser than a miss. - let segments = vcad_render::tessellation_segments(raster_size_px(cli, format)); + let segments = cache_segments(cli, format); let r = vcad_render::with_root_cache(cache.clone(), segments, || { render_one_uncached(input, dest, format, cli) }); diff --git a/crates/vcad-render/src/photoreal.rs b/crates/vcad-render/src/photoreal.rs index 4de82b778..7be79e697 100644 --- a/crates/vcad-render/src/photoreal.rs +++ b/crates/vcad-render/src/photoreal.rs @@ -3,13 +3,16 @@ //! Where the drafting path projects tessellated triangles onto a tonal ramp //! and the `--raytrace` path swaps in analytic intersection but keeps the //! same ramp, this path solves the rendering equation properly: a -//! physically-based path tracer over the untessellated BRep, lit by a -//! three-softbox studio rig and an analytic sky, viewed through a camera with -//! a real focal length and aperture. +//! physically-based path tracer lit by a three-softbox studio rig and an +//! analytic sky, viewed through a camera with a real focal length and +//! aperture. //! -//! The geometry advantage carries over — silhouettes and specular highlights -//! on fillets come from analytic ray–surface intersection, so they are exact -//! at any resolution with no facet banding. +//! **Geometry comes from a tessellation by default** ([`MESH_SEGMENTS`]), so +//! the content-addressed root-mesh cache can serve it and a re-render pays +//! only for path tracing rather than re-evaluating every BRep. `--exact` +//! ([`PhotorealOptions::exact`]) swaps in analytic ray–surface intersection, +//! which is exact at any resolution with no facet banding on curved +//! silhouettes, at the price of a full kernel evaluation every run. use std::sync::Arc; @@ -71,6 +74,74 @@ pub struct PhotorealOptions { /// On by default: it is worth far more per second of render time than the /// equivalent extra samples. Turn it off for reference renders. pub denoise: bool, + /// Stop sampling a pixel early once its own variance says the rest of the + /// budget cannot move it visibly (`--no-adaptive` turns this off). + /// + /// On by default. `spp` becomes a ceiling rather than a fixed count; the + /// stopping decision is per-pixel and made from that pixel's own running + /// sums, so a fixed `--seed` still gives a byte-stable image. + pub adaptive: bool, + /// Intersect the analytic BRep surfaces instead of a tessellation of them + /// (`--exact`). + /// + /// **Off by default.** A `vcad_kernel::Solid` has no serialized form, so + /// a BRep can only ever come from re-evaluating the document — which for + /// a real assembly is most of the wall time of a moderate-`--spp` render, + /// and is paid again on every invocation. Tracing triangles at + /// [`MESH_SEGMENTS`] instead lets the content-addressed root-mesh cache + /// (`vcad_eval::cache`) supply the geometry, leaving only the path + /// tracing to redo. + /// + /// The cost is curved silhouettes. At 128 segments adjacent facets differ + /// by 2.8°, which is invisible on flat sheet and on large-radius bosses, + /// but reads as a stair-step along the rim of a small cylinder at hero + /// resolutions — on rose-pro at `--size 1200` it is plain on the D55 + /// actuator cans. That is what `exact` is for. + /// + /// Two further differences are worth knowing before diffing two renders: + /// + /// * **Framing shifts slightly.** The camera fits the objects' BVH root + /// AABBs, and a tessellation is inscribed in the surface it approximates, + /// so the subject's bounds are a hair smaller in mesh mode and the fit + /// lands at a marginally different scale. Nothing is clipped; the two + /// images simply do not overlay pixel-for-pixel anywhere, not just on + /// curves. + /// * **Shading is not uniformly worse.** Mesh mode goes through + /// `render_bake`'s crease-aware normals, which resolve plate edges more + /// cleanly than analytic evaluation of the same BRep does; `exact` in + /// exchange shows faint banding across some large planar faces. Mesh + /// mode loses on curves and wins on creases. + /// + /// Note this is a *pixel* switch, not only a speed one, and it is + /// deliberately not conditioned on whether a cache happens to be + /// installed: a render must not change its output depending on the state + /// of an accelerator. Mesh mode therefore tessellates on a cache miss + /// too, so a cold render and a warm one are identical. + pub exact: bool, +} + +/// Segment count the default (non-[`exact`](PhotorealOptions::exact)) +/// photoreal path tessellates curved faces at. +/// +/// This is the raster path's hi-res count, but — unlike the raster path — it +/// is used at *every* canvas size rather than scaling with it. Two reasons: +/// the path tracer resolves a silhouette far more sharply than the flat +/// raster shader (a specular highlight tracks facet normals; a tonal ramp +/// mostly doesn't), so the coarser 64-segment count shows at sizes the raster +/// path gets away with; and the extra triangles cost almost nothing here, +/// since BVH traversal is logarithmic in triangle count while the +/// tessellation itself is cached. +/// +/// A size-independent count also means one cache entry per root serves every +/// `--size`, instead of one entry per size bucket. +pub const MESH_SEGMENTS: u32 = crate::raster::RASTER_SEGMENTS_HIRES; + +impl PhotorealOptions { + /// Segment count to tessellate BRep-backed solids at before tracing, or + /// `None` when [`Self::exact`] asks for analytic surfaces. + pub(crate) fn mesh_segments(&self) -> Option { + (!self.exact).then_some(MESH_SEGMENTS) + } } impl Default for PhotorealOptions { @@ -87,6 +158,8 @@ impl Default for PhotorealOptions { backdrop: Backdrop::Studio, seed: 0x5eed_1234, denoise: true, + adaptive: true, + exact: false, } } } @@ -115,23 +188,50 @@ pub fn render_photoreal_png_str( /// Build one BVH per solid, world-placed at the identity. /// -/// BRep-backed solids trace analytically; mesh-only parts (frozen -/// topology-optimization results, imported STL/GLB) trace as crease-baked -/// triangles, same as the `--raytrace` path. Fails closed when any part has -/// no traceable geometry — a silently missing part reads as a design that -/// doesn't have it. -pub(crate) fn build_objects(solids: &[crate::SceneSolid]) -> Result, String> { - let mut objects: Vec = Vec::new(); - let mut untraceable: Vec = Vec::new(); - for s in solids { - let bvh = match s.solid.as_brep() { - Some(brep) => Bvh::build(brep), - None => { - let mut mesh = s.solid.to_mesh(0); +/// `mesh_segments` picks the geometry representation, and comes from +/// [`PhotorealOptions::mesh_segments`]: +/// +/// * `Some(n)` — the default. Every solid traces as crease-baked triangles, +/// BRep-backed ones tessellated at `n` segments. A solid that arrived +/// mesh-backed (a root-mesh cache hit, a frozen topology-optimization +/// result, an imported STL/GLB) already carries its triangles and `n` is +/// moot for it, which is exactly why cold and warm renders agree. +/// * `None` — `--exact`. BRep-backed solids intersect analytically; mesh-only +/// parts still trace as triangles, since there is nothing else to trace. +/// +/// Fails closed when any part has no traceable geometry — a silently missing +/// part reads as a design that doesn't have it. +/// +/// BVH construction is the one genuinely parallel-friendly stage of scene +/// setup (each solid is independent, and for a many-root assembly it is the +/// bulk of the time between evaluation and the first traced ray), so the +/// builds fan out over rayon. Results are re-zipped with `solids` afterwards, +/// so object order — and therefore the untraceable-part diagnostics — is +/// unchanged from the serial version. +pub(crate) fn build_objects( + solids: &[crate::SceneSolid], + mesh_segments: Option, +) -> Result, String> { + use rayon::prelude::*; + + let bvhs: Vec = solids + .par_iter() + .map(|s| match (mesh_segments, s.solid.as_brep()) { + (None, Some(brep)) => Bvh::build(brep), + (segments, _) => { + // `to_mesh` ignores the segment count for a mesh-backed + // solid, so `unwrap_or(0)` keeps --exact's mesh-only branch + // byte-identical to what it did before this parameter existed. + let mut mesh = s.solid.to_mesh(segments.unwrap_or(0)); vcad_kernel::vcad_kernel_tessellate::render_bake_default(&mut mesh); Bvh::build_mesh(&mesh) } - }; + }) + .collect(); + + let mut objects: Vec = Vec::new(); + let mut untraceable: Vec = Vec::new(); + for (s, bvh) in solids.iter().zip(bvhs) { if bvh.root().is_none() { untraceable.push(s.name.clone().unwrap_or_else(|| s.id.clone())); continue; @@ -352,6 +452,7 @@ pub(crate) fn trace_options(pr: &PhotorealOptions, png: bool) -> PathTraceOption show_background: !png || pr.backdrop == Backdrop::Studio, seed: pr.seed, denoise: pr.denoise, + adaptive: pr.adaptive, ..PathTraceOptions::default() } } @@ -397,7 +498,7 @@ fn rasterize( if solids.is_empty() { return Err("no solids produced".to_string()); } - let objects = build_objects(&solids)?; + let objects = build_objects(&solids, pr.mesh_segments())?; let corners: Vec<[f64; 3]> = objects.iter().flat_map(object_corners).collect(); let framing = frame_view(&corners, opts, pr)?; let scene = dress_scene(objects, &framing, pr)?; diff --git a/crates/vcad-render/src/photoreal_gpu.rs b/crates/vcad-render/src/photoreal_gpu.rs new file mode 100644 index 000000000..fcad92b03 --- /dev/null +++ b/crates/vcad-render/src/photoreal_gpu.rs @@ -0,0 +1,498 @@ +//! `--photoreal --gpu`: the photoreal path traced on a wgpu compute pipeline. +//! +//! This is the *same scene* as [`crate::photoreal`], not a second renderer +//! with its own opinions. Geometry comes from `photoreal::build_objects` (one +//! cached triangle BLAS per solid), framing from `photoreal::frame_view`, and +//! lights/environment/floor from `photoreal::dress_scene`; the film goes back +//! through `pathtrace::Film::to_srgb8`, so exposure, ACES and sRGB encoding +//! are byte-for-byte the CPU path's. What changes is only who integrates the +//! rendering equation. +//! +//! # What `--gpu` honours, and what it refuses +//! +//! The GPU tracer is a narrower renderer than the CPU one, and the honest +//! thing is to say so in the error message rather than quietly render +//! something else. Refused outright, with a message naming the flag: +//! +//! * `--exact` — the GPU's BRep path ([`GpuScene::from_brep`]) caps at a +//! thousand analytic surfaces, which is a bracket, not an assembly. `--gpu` +//! traces the cached tessellation, exactly as CPU `--photoreal` does by +//! default. +//! * `--aperture` — the WGSL camera is a pinhole; there is no lens to sample. +//! * `--ortho` — the WGSL camera is projective only. +//! * `--backdrop shadow-catcher` — the shadow catcher is a CPU integrator +//! feature (a surface that contributes occlusion to alpha and nothing to +//! colour); there is no shader counterpart. +//! * `--animate` — transforms are *baked* into the uploaded vertices (see +//! [`GpuScene::from_mesh_bvh_placed`]), so a new pose means a new upload. +//! A sequence is exactly the case where that is the wrong trade. +//! +//! Honoured, and matching the CPU path: `--spp`, `--max-depth`, `--exposure`, +//! `--fov`, `--seed`, `--size`, `--fill`, `--auto-aspect`, `--view` / +//! `--azimuth` / `--elevation` (including the mirrored isometric basis), +//! `--env` and `--env-rotation` (gradient and HDRI both), per-part materials, +//! and `--backdrop studio` / `--backdrop none`. +//! +//! Two knowing divergences, both documented where they happen: +//! +//! * **The studio floor is a large quad, not an infinite plane.** The CPU +//! `Ground` is analytic and unbounded; the shader has no such primitive at +//! an arbitrary height, so the floor is uploaded as two triangles +//! [`GROUND_EXTENT`] scene-radii across. Within the frame of any sane +//! product shot the two agree; a camera aimed at the horizon would see +//! where the quad stops. +//! * **No denoiser.** See [`denoise_is_unavailable`]. +//! +//! `--no-adaptive` is accepted and ignored: adaptive sampling is a per-pixel +//! early-out the shader has no equivalent for, so `--spp` on the GPU is an +//! exact count rather than a ceiling. +//! +//! # How close the two get +//! +//! Close, and not identical, and the gap does not close with samples. On +//! rose-pro at 800px the GPU render sits ~29 dB PSNR from a 1024-spp CPU +//! reference at 64 spp and ~30 dB at 1024 spp — it converges, but to a +//! slightly different picture, because the integrator is f32 where the CPU's +//! is f64 and because of the floor above. A CPU render at the same 64 spp +//! reaches ~39 dB. So `--gpu` is the right tool for a fast look or a large +//! sweep, and the CPU path is still the one to render a final hero on. +//! +//! Setting `VCAD_GPU_DEBUG` prints the packed scene's part count, triangle +//! count, node count, BVH depth, light count and environment mode to stderr — +//! the numbers every "why was this refused" question turns out to be about. + +use vcad_kernel::vcad_kernel_math::{Point3, Transform, Vec3}; +use vcad_kernel::vcad_kernel_tessellate::TriangleMesh; +use vcad_kernel_gpu::{GpuContext, GpuError}; +use vcad_kernel_raytrace::gpu::{ + GpuAreaLight, GpuCamera, GpuMaterial, GpuScene, OfflineOptions, RayTracePipeline, +}; +use vcad_kernel_raytrace::pathtrace::{Environment, Ground, Object, Scene}; +use vcad_kernel_raytrace::Bvh; + +use super::photoreal::{self, Backdrop, Framing, PhotorealOptions}; +use super::raster::{encode_jpeg, encode_png, Frame}; +use super::{evaluate_vcad, RasterOptions}; + +/// Half-width of the studio floor quad, in scene radii. +/// +/// The CPU renderer's floor is an unbounded analytic plane. The shader has no +/// plane-at-arbitrary-height primitive (its built-in ground is pinned to +/// z = 0 and faded by distance, which is a viewport look, not this one), so +/// the floor ships as two triangles instead. +/// +/// The value is a measured compromise, not a guess. Too small and the quad's +/// edge walks into frame. Too large and f32 loses the plane: at 600 radii a +/// rose-pro render came back with visible wavy banding across the floor — +/// self-shadowing acne from ray/triangle arithmetic on coordinates six +/// hundred times the size of the millimetre geometry beside them. 50 radii +/// puts the edge comfortably outside the frame of any camera actually aimed +/// at the subject, and the banding is gone. +pub const GROUND_EXTENT: f64 = 50.0; + +/// Why the CPU denoiser is not run on a GPU film. +/// +/// `pathtrace::denoise` is guided by the per-pixel normal, depth and albedo +/// the CPU integrator records from each primary ray. `render_offline` reads +/// back the radiance accumulation buffer and nothing else, so +/// `OfflineResult::to_film` leaves those guides zeroed — and the denoiser +/// treats `depth == 0` as "this pixel escaped to the background, pass it +/// through untouched". Run against a zeroed film it is therefore not +/// *wrong*, it is a silent no-op: the user asks for denoising, waits for it, +/// and gets the raw film back with no indication. +/// +/// So `--gpu` disables it and says so on stderr, once. Writing the guide +/// buffers from the shader is the obvious follow-up; it is a real slice of +/// work (three more storage buffers and a first-sample-only write path), not +/// something to fake here. +pub const fn denoise_is_unavailable() -> &'static str { + "--photoreal --gpu: the denoiser needs the normal/depth/albedo guide \ + buffers, which the GPU tracer does not read back; rendering without it \ + (raise --spp to compensate, or drop --gpu)" +} + +/// Render raw `.vcad` document JSON to a photorealistic JPEG on the GPU. +pub fn render_photoreal_gpu_jpeg_str( + raw_vcad: &str, + opts: &RasterOptions, + pr: &PhotorealOptions, +) -> Result, String> { + encode_jpeg(rasterize(raw_vcad, opts, pr, false)?, opts) +} + +/// Render raw `.vcad` document JSON to a photorealistic RGBA PNG on the GPU. +pub fn render_photoreal_gpu_png_str( + raw_vcad: &str, + opts: &RasterOptions, + pr: &PhotorealOptions, +) -> Result, String> { + encode_png(rasterize(raw_vcad, opts, pr, true)?, opts) +} + +/// Reject the option combinations the GPU tracer cannot render *faithfully*. +/// +/// Every one of these would otherwise produce a plausible-looking image that +/// is not the render the user asked for, which is the worst failure mode a +/// renderer has. Called before any GPU work, so the diagnostic arrives +/// immediately rather than after a scene build. +pub fn check_supported(opts: &RasterOptions, pr: &PhotorealOptions) -> Result<(), String> { + let _ = opts; + if pr.exact { + return Err("--gpu does not support --exact: the GPU tracer's analytic \ + BRep path is capped at ~1k surfaces, far below a real \ + assembly. Drop --exact to trace the cached tessellation \ + (what --photoreal does by default), or drop --gpu." + .to_string()); + } + if pr.aperture_frac > 0.0 { + return Err(format!( + "--gpu does not support --aperture (got {}): the GPU camera is a \ + pinhole, so depth of field would silently render sharp. Drop \ + --aperture, or drop --gpu.", + pr.aperture_frac + )); + } + if pr.orthographic { + return Err("--gpu does not support --ortho: the GPU camera is \ + projective only, and would render a perspective image \ + under an orthographic flag. Drop --ortho, or drop --gpu." + .to_string()); + } + if pr.backdrop == Backdrop::ShadowCatcher { + return Err("--gpu does not support --backdrop shadow-catcher: the \ + shadow catcher is a CPU integrator feature (a surface \ + that contributes to alpha but not to colour) with no \ + shader counterpart. Use --backdrop studio or none, or \ + drop --gpu." + .to_string()); + } + Ok(()) +} + +/// Acquire the GPU, or explain why not. +/// +/// `--gpu` is an explicit request, so an unavailable adapter is a hard error +/// rather than a quiet fallback: a user who asked for the GPU and silently +/// got a CPU render would draw the wrong conclusion from every timing they +/// took afterwards. The message names the fallback rather than taking it. +fn context() -> Result<&'static GpuContext, String> { + match pollster::block_on(GpuContext::init()) { + Ok(ctx) => Ok(ctx), + Err(GpuError::NoAdapter) => Err("--gpu: no compatible GPU adapter found. \ + Drop --gpu to path-trace on the CPU." + .to_string()), + Err(e) => Err(format!( + "--gpu: could not initialise the GPU ({e}). \ + Drop --gpu to path-trace on the CPU." + )), + } +} + +/// Two triangles at `z`, centred on `center`, spanning `half` in x and y. +/// +/// Wound counter-clockwise seen from +z and given explicit up-normals, so it +/// shades as a floor rather than taking the geometric-normal fallback. +fn ground_mesh(center: Point3, z: f64, half: f64) -> TriangleMesh { + let (x0, x1) = ((center.x - half) as f32, (center.x + half) as f32); + let (y0, y1) = ((center.y - half) as f32, (center.y + half) as f32); + let z = z as f32; + TriangleMesh { + vertices: vec![x0, y0, z, x1, y0, z, x1, y1, z, x0, y1, z], + indices: vec![0, 1, 2, 0, 2, 3], + normals: vec![0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 0.0, 1.0], + face_kinds: Vec::new(), + } +} + +/// Pack one CPU [`Object`] as its own single-material GPU scene. +fn object_scene(obj: &Object) -> Result { + // An identity transform is the overwhelmingly common case (the static + // photoreal path never places anything); skipping the bake there keeps + // the packed vertices bit-identical to the un-transformed ones rather + // than round-tripping every coordinate through an f64 matrix multiply. + let transform = (!is_identity(&obj.transform)).then_some(&obj.transform); + GpuScene::from_mesh_bvh_placed(&obj.bvh, GpuMaterial::from_pbr(obj.material), transform) + .map_err(|e| format!("--gpu: cannot upload part geometry: {e}")) +} + +/// Whether `t` moves nothing. +/// +/// Probed rather than read out of the matrix: an affine map is pinned by +/// where it sends the origin and the three basis points, so this is exact +/// without depending on the matrix's storage order. +fn is_identity(t: &Transform) -> bool { + const EPS: f64 = 1e-12; + [ + Point3::new(0.0, 0.0, 0.0), + Point3::new(1.0, 0.0, 0.0), + Point3::new(0.0, 1.0, 0.0), + Point3::new(0.0, 0.0, 1.0), + ] + .iter() + .all(|p| { + let q = t.apply_point(p); + (q.x - p.x).abs() < EPS && (q.y - p.y).abs() < EPS && (q.z - p.z).abs() < EPS + }) +} + +/// Turn the CPU-side [`Scene`] into a single validated [`GpuScene`]. +/// +/// The order here matters. `GpuScene::merge` re-derives the studio rig from +/// the combined bounds, which is right for the subject and very wrong once the +/// floor quad (hundreds of radii across) joins the scene — so the rig is +/// written once, after the last merge. The CPU path sizes its rig the same way: +/// `dress_scene` works from the subject's bounds, not the floor's. +fn build_scene(scene: &Scene, framing: &Framing) -> Result { + let parts: Vec = scene + .objects + .iter() + .map(object_scene) + .collect::>()?; + let mut merged = + GpuScene::merge_all(parts).ok_or_else(|| "--gpu: scene has no geometry".to_string())?; + + if let Some(ground) = &scene.ground { + merged = merged.merge(ground_scene(ground, framing)?); + } + + // After the last merge, never before — see above. The lights themselves + // come from `dress_scene`, which is where the CPU renderer gets them, + // including the empty rig an HDRI environment implies: an image + // environment already carries its own lighting. + merged.lights = scene + .lights + .iter() + .map(GpuAreaLight::from_area_light) + .collect(); + merged.set_environment(match &scene.env { + Environment::Image(map) => Some(map), + // The shader's analytic gradient is a transcription of + // `GradientEnv::default()`; passing None selects it. + Environment::Gradient(_) => None, + }); + Ok(merged) +} + +fn ground_scene(ground: &Ground, framing: &Framing) -> Result { + let mesh = ground_mesh(framing.center, ground.z, framing.radius * GROUND_EXTENT); + GpuScene::from_mesh_bvh_placed( + &Bvh::build_mesh(&mesh), + GpuMaterial::from_pbr(ground.material), + None, + ) + .map_err(|e| format!("--gpu: cannot upload the studio floor: {e}")) +} + +/// The GPU camera for a CPU [`Framing`]. +/// +/// Hands the screen basis over verbatim. `View::Isometric` and the named CAD +/// views carry a *mirrored* basis (see `View::right`), and the shader's +/// default reconstruction — `right = forward × up` — cannot represent one, so +/// rebuilding it there flips the render left-for-right against every other +/// output style. `CAMERA_BASIS_EXPLICIT` is exactly this fix. +fn camera_for(framing: &Framing, w: u32, h: u32) -> GpuCamera { + let c = &framing.camera; + let f32v = |v: Vec3| [v.x as f32, v.y as f32, v.z as f32]; + GpuCamera::from_basis( + [c.eye.x as f32, c.eye.y as f32, c.eye.z as f32], + f32v(c.forward), + f32v(c.right), + f32v(c.up), + (c.fov_deg as f32).to_radians(), + c.focus_dist as f32, + w, + h, + ) +} + +fn rasterize( + raw_vcad: &str, + opts: &RasterOptions, + pr: &PhotorealOptions, + png: bool, +) -> Result { + photoreal::check_raster_opts(opts)?; + check_supported(opts, pr)?; + + let ctx = context()?; + let pipeline = + RayTracePipeline::new(ctx).map_err(|e| format!("--gpu: pipeline creation failed: {e}"))?; + + let solids = evaluate_vcad(raw_vcad)?; + if solids.is_empty() { + return Err("no solids produced".to_string()); + } + // `mesh_segments()` is Some here: --exact is refused above. + let objects = photoreal::build_objects(&solids, pr.mesh_segments())?; + let corners: Vec<[f64; 3]> = objects.iter().flat_map(photoreal::object_corners).collect(); + let framing = photoreal::frame_view(&corners, opts, pr)?; + let scene = photoreal::dress_scene(objects, &framing, pr)?; + + let gpu_scene = build_scene(&scene, &framing)?; + if std::env::var_os("VCAD_GPU_DEBUG").is_some() { + eprintln!( + "--gpu scene: {} parts, {} triangles, {} bvh nodes, depth {}, {} lights, env {}", + scene.objects.len(), + gpu_scene.surfaces.len(), + gpu_scene.bvh_nodes.len(), + gpu_scene.bvh_depth(), + gpu_scene.lights.len(), + if gpu_scene.environment.is_some() { + "image" + } else { + "gradient" + }, + ); + } + gpu_scene + .validate(Some( + ctx.device.limits().max_storage_buffer_binding_size as u64, + )) + .map_err(|e| format!("--gpu: {e}"))?; + + let canvas = framing.canvas; + let (w, h) = (canvas.w as u32, canvas.h as u32); + let cpu_opts = photoreal::trace_options(pr, png); + let out = pipeline + .render_offline( + ctx, + &gpu_scene, + &camera_for(&framing, w, h), + &OfflineOptions { + width: w, + height: h, + spp: cpu_opts.spp, + max_depth: cpu_opts.max_depth, + rr_start: cpu_opts.rr_start, + // Only read for the gradient: an image environment carries + // its own intensity, which `offline_render_state` takes from + // the packed map. + env_intensity: match &scene.env { + Environment::Gradient(g) => g.intensity, + Environment::Image(_) => 1.0, + }, + firefly_clamp: cpu_opts.firefly_clamp.unwrap_or(0.0), + // The floor is real geometry here, so the shader's own + // implicit z=0 ground must stay out of the way — it would + // otherwise add a second, differently-placed floor. + ground_enabled: false, + show_background: cpu_opts.show_background, + // The CPU seed is 64-bit and only ever mixed into a per-pixel + // hash; the shader's is 32. Fold rather than truncate so two + // seeds that differ only in the high word still differ here. + seed: (pr.seed as u32) ^ ((pr.seed >> 32) as u32), + }, + ) + .map_err(|e| format!("--gpu: render failed: {e}"))?; + + let film = out.to_film(); + let rgba = film.to_srgb8(pr.exposure, png && pr.backdrop != Backdrop::Studio); + + let n = canvas.len(); + let mut rgb = vec![0u8; n * 3]; + let mut mask = vec![0u8; n]; + for i in 0..n { + rgb[i * 3] = rgba[i * 4]; + rgb[i * 3 + 1] = rgba[i * 4 + 1]; + rgb[i * 3 + 2] = rgba[i * 4 + 2]; + mask[i] = rgba[i * 4 + 3]; + } + Ok(Frame { rgb, mask, canvas }) +} + +#[cfg(test)] +mod tests { + use super::*; + + fn opts() -> RasterOptions { + RasterOptions { + size_px: 64, + ..Default::default() + } + } + + #[test] + fn exact_is_refused() { + let pr = PhotorealOptions { + exact: true, + ..Default::default() + }; + let err = check_supported(&opts(), &pr).expect_err("should refuse"); + assert!(err.contains("--exact"), "unhelpful: {err}"); + } + + #[test] + fn aperture_is_refused() { + let pr = PhotorealOptions { + aperture_frac: 0.03, + ..Default::default() + }; + let err = check_supported(&opts(), &pr).expect_err("should refuse"); + assert!(err.contains("--aperture"), "unhelpful: {err}"); + } + + #[test] + fn ortho_is_refused() { + let pr = PhotorealOptions { + orthographic: true, + ..Default::default() + }; + let err = check_supported(&opts(), &pr).expect_err("should refuse"); + assert!(err.contains("--ortho"), "unhelpful: {err}"); + } + + #[test] + fn shadow_catcher_is_refused() { + let pr = PhotorealOptions { + backdrop: Backdrop::ShadowCatcher, + ..Default::default() + }; + let err = check_supported(&opts(), &pr).expect_err("should refuse"); + assert!(err.contains("shadow-catcher"), "unhelpful: {err}"); + } + + #[test] + fn the_defaults_are_supported() { + check_supported(&opts(), &PhotorealOptions::default()).expect("defaults must work on GPU"); + } + + /// The floor quad must actually be under the subject and wider than it, + /// or the "studio sweep" is a stripe. + #[test] + fn ground_quad_spans_the_subject() { + let mesh = ground_mesh(Point3::new(5.0, -2.0, 10.0), 1.5, 100.0); + assert_eq!(mesh.indices.len(), 6); + let xs: Vec = mesh + .vertices + .as_chunks::<3>() + .0 + .iter() + .map(|v| v[0]) + .collect(); + let zs: Vec = mesh + .vertices + .as_chunks::<3>() + .0 + .iter() + .map(|v| v[2]) + .collect(); + assert!( + zs.iter().all(|&z| (z - 1.5).abs() < 1e-6), + "floor is not flat" + ); + assert!( + xs.iter().cloned().fold(f32::MIN, f32::max) >= 105.0, + "floor does not reach past the subject" + ); + assert!( + mesh.normals.as_chunks::<3>().0.iter().all(|n| n[2] > 0.9), + "floor normals must point up" + ); + } + + #[test] + fn identity_transform_is_detected() { + assert!(is_identity(&Transform::identity())); + assert!(!is_identity(&Transform::translation(0.0, 0.0, 1e-3))); + } +} diff --git a/crates/vcad-render/tests/photoreal_gpu.rs b/crates/vcad-render/tests/photoreal_gpu.rs new file mode 100644 index 000000000..d22410753 --- /dev/null +++ b/crates/vcad-render/tests/photoreal_gpu.rs @@ -0,0 +1,234 @@ +//! `--photoreal --gpu` against `--photoreal`, through the full `vcad-render` +//! entry points. +//! +//! The unit tests inside `photoreal_gpu.rs` cover the refusals; this is the +//! one that costs a GPU and answers the only question that matters: does the +//! GPU path render *the same picture* the CPU path does? Same document, same +//! seed, same sample count, same framing — compared as PSNR over the encoded +//! sRGB pixels, which is what a viewer actually sees. +//! +//! Run with: +//! +//! ```text +//! cargo test -p vcad-render --features photoreal-gpu --test photoreal_gpu +//! ``` +//! +//! Skipped, not failed, when there is no adapter: CI machines without a GPU +//! must not turn red over a feature they cannot exercise. + +#![cfg(feature = "photoreal-gpu")] + +use vcad_render::photoreal::{Backdrop, PhotorealOptions}; +use vcad_render::photoreal_gpu::render_photoreal_gpu_png_str; +use vcad_render::RasterOptions; + +/// A stepped block: flat faces at several angles, a through hole, and enough +/// self-shadowing that a lighting disagreement between the two integrators +/// would show up rather than averaging out. +fn doc() -> String { + r#"{ + "version": "0.1", + "nodes": { + "1": { "id": 1, "name": "base", + "op": { "type": "Cube", "size": { "x": 30, "y": 20, "z": 6 } } }, + "2": { "id": 2, "name": "post", + "op": { "type": "Cylinder", "radius": 5, "height": 18, "segments": 64 } } + }, + "materials": { + "aluminum": { "name": "aluminum", "color": [0.91, 0.92, 0.93], + "metallic": 1.0, "roughness": 0.4 }, + "abs": { "name": "abs", "color": [0.15, 0.16, 0.18], + "metallic": 0.0, "roughness": 0.55 } + }, + "part_materials": {}, + "roots": [{ "root": 1, "material": "aluminum" }, + { "root": 2, "material": "abs" }] +}"# + .to_string() +} + +fn raster_opts(size: u32) -> RasterOptions { + RasterOptions { + size_px: size, + ..Default::default() + } +} + +/// Both renderers get an *identical* brief: the same sample count spent on +/// every pixel (no adaptive early-out, which the GPU has no equivalent for), +/// and no denoiser (which the GPU cannot run at all). Comparing a denoised +/// CPU image against a raw GPU one would measure the denoiser, not the +/// integrators. +fn photoreal_opts(spp: u32) -> PhotorealOptions { + PhotorealOptions { + spp, + seed: 7, + denoise: false, + adaptive: false, + backdrop: Backdrop::Studio, + ..Default::default() + } +} + +struct Image { + w: u32, + h: u32, + rgb: Vec, +} + +fn decode(png: &[u8]) -> Image { + let img = image::load_from_memory(png).expect("valid PNG").to_rgb8(); + Image { + w: img.width(), + h: img.height(), + rgb: img.into_raw(), + } +} + +/// Peak signal-to-noise ratio over the 8-bit sRGB samples, in dB. Infinite +/// for identical images, which no two integrators ever are. +fn psnr(a: &Image, b: &Image) -> f64 { + assert_eq!((a.w, a.h), (b.w, b.h), "size mismatch"); + let mse: f64 = a + .rgb + .iter() + .zip(&b.rgb) + .map(|(&x, &y)| { + let d = x as f64 - y as f64; + d * d + }) + .sum::() + / a.rgb.len() as f64; + if mse <= 0.0 { + return f64::INFINITY; + } + 10.0 * (255.0f64 * 255.0 / mse).log10() +} + +/// Render on the GPU, or `None` when this machine has no adapter. +/// +/// The skip is keyed on the specific "no adapter" message the `--gpu` path +/// produces. Every *other* failure is a real failure and is allowed to +/// propagate: a test that skipped on any error at all would go green on a +/// broken shader. +fn gpu_png(opts: &RasterOptions, pr: &PhotorealOptions) -> Option> { + match render_photoreal_gpu_png_str(&doc(), opts, pr) { + Ok(png) => Some(png), + Err(e) if e.contains("no compatible GPU adapter") => { + eprintln!("skipped: {e}"); + None + } + Err(e) => panic!("GPU render failed: {e}"), + } +} + +/// The gate. 24 dB is deliberately below the ~28-30 dB the two paths actually +/// hit on real scenes: at a test-sized sample count both images are still +/// visibly noisy, and the two integrators' noise is independent, so a chunk +/// of the measured error here is Monte Carlo disagreement rather than a +/// difference of opinion about the scene. Anything that gets the geometry, +/// the framing, the lighting or the tonemap wrong lands far below this — a +/// mirrored camera basis scores in the low teens, a dropped part or a missing +/// floor worse still. +const MIN_PSNR: f64 = 24.0; + +#[test] +fn gpu_render_matches_the_cpu_render() { + let opts = raster_opts(128); + let pr = photoreal_opts(96); + + let Some(gpu) = gpu_png(&opts, &pr) else { + return; + }; + let cpu = + vcad_render::photoreal::render_photoreal_png_str(&doc(), &opts, &pr).expect("CPU render"); + + let (gpu, cpu) = (decode(&gpu), decode(&cpu)); + let db = psnr(&gpu, &cpu); + eprintln!("CPU vs GPU: {db:.2} dB"); + assert!( + db >= MIN_PSNR, + "GPU render is {db:.2} dB from the CPU render (floor {MIN_PSNR} dB) -- \ + the two paths disagree about the scene, not merely about noise" + ); +} + +/// The subject must land in the same place, not merely look similar on +/// average. A mirrored camera basis — which is exactly what the GPU's default +/// right-handed reconstruction does to `View::Isometric` — leaves the overall +/// brightness untouched and moves every pixel, so a whole-image metric is a +/// weak detector for it and a column profile is a strong one. +#[test] +fn the_isometric_view_is_not_mirrored() { + let opts = RasterOptions { + size_px: 128, + view: vcad_render::View::Isometric, + ..Default::default() + }; + let pr = photoreal_opts(48); + + let Some(gpu) = gpu_png(&opts, &pr) else { + return; + }; + let cpu = + vcad_render::photoreal::render_photoreal_png_str(&doc(), &opts, &pr).expect("CPU render"); + let (gpu, cpu) = (decode(&gpu), decode(&cpu)); + + // Per-column mean luminance: a silhouette signature that survives noise + // but not a left-right flip. + let profile = |img: &Image| -> Vec { + (0..img.w) + .map(|x| { + (0..img.h) + .map(|y| { + let i = ((y * img.w + x) * 3) as usize; + 0.2126 * img.rgb[i] as f64 + + 0.7152 * img.rgb[i + 1] as f64 + + 0.0722 * img.rgb[i + 2] as f64 + }) + .sum::() + / img.h as f64 + }) + .collect() + }; + let (pg, pc) = (profile(&gpu), profile(&cpu)); + let err = |a: &[f64], b: &[f64]| -> f64 { + a.iter().zip(b).map(|(x, y)| (x - y).abs()).sum::() / a.len() as f64 + }; + + let mut flipped = pg.clone(); + flipped.reverse(); + let (straight, mirrored) = (err(&pg, &pc), err(&flipped, &pc)); + assert!( + straight < mirrored, + "the GPU isometric render matches the CPU one better MIRRORED \ + ({mirrored:.3}) than straight ({straight:.3}) -- the camera basis is \ + being rebuilt right-handedly instead of handed over verbatim" + ); +} + +/// `--backdrop none` must leave the background transparent, exactly as the +/// CPU path does. This is the shader's `background_mode` reaching the +/// coverage alpha; if it did not, the escaped rays would paint the viewport's +/// themed sky and the PNG would be opaque. +#[test] +fn backdrop_none_leaves_the_background_transparent() { + let opts = raster_opts(96); + let pr = PhotorealOptions { + backdrop: Backdrop::None, + ..photoreal_opts(16) + }; + let Some(png) = gpu_png(&opts, &pr) else { + return; + }; + let img = image::load_from_memory(&png).expect("valid PNG").to_rgba8(); + + let clear = img.pixels().filter(|p| p.0[3] < 8).count(); + let opaque = img.pixels().filter(|p| p.0[3] > 200).count(); + assert!( + clear > 0 && opaque > 0, + "expected a mix of covered and transparent pixels, got \ + {clear} clear / {opaque} opaque of {}", + img.pixels().len() + ); +} diff --git a/scripts/photoreal-quality.sh b/scripts/photoreal-quality.sh new file mode 100755 index 000000000..b42ed508c --- /dev/null +++ b/scripts/photoreal-quality.sh @@ -0,0 +1,117 @@ +#!/usr/bin/env bash +# Photoreal image-quality gate. +# +# Renders each scene twice — a high-spp reference and a candidate at the +# normal sample count — and compares them with the `psnr` example. Sampling +# changes in crates/vcad-kernel-raytrace/src/pathtrace.rs must keep every +# scene above MIN_PSNR. +# +# References are large PNGs and are deliberately kept OUT of the repo, under +# $REF_DIR (default /tmp/vcad-photoreal-ref). They are regenerated only when +# missing, or when --regen is passed, so a candidate sweep costs one cheap +# render per scene. +# +# Usage: +# scripts/photoreal-quality.sh # gate current build +# scripts/photoreal-quality.sh --regen # rebuild references first +# REF_SPP=1024 CAND_SPP=32 scripts/photoreal-quality.sh +# +# The candidate render command can be extended with EXTRA_ARGS, which is how +# a flag under test (e.g. --no-adaptive) gets swept against the same +# references. +set -euo pipefail + +ROOT="$(cd "$(dirname "${BASH_SOURCE[0]}")/.." && pwd)" +cd "$ROOT" + +REF_DIR="${REF_DIR:-/tmp/vcad-photoreal-ref}" +OUT_DIR="${OUT_DIR:-/tmp/vcad-photoreal-cand}" +REF_SPP="${REF_SPP:-1024}" +CAND_SPP="${CAND_SPP:-32}" +SIZE="${SIZE:-800}" +SEED="${SEED:-7}" +MIN_PSNR="${MIN_PSNR:-35}" +EXTRA_ARGS="${EXTRA_ARGS:-}" + +REGEN=0 +[ "${1:-}" = "--regen" ] && REGEN=1 + +# BIN can point at a binary built from another revision, which is how a +# candidate is A/B'd against a baseline over one shared set of references. +# Setting it also skips the build, so the comparison binary is not clobbered. +BIN="${BIN:-}" +PSNR="$ROOT/target/release/examples/psnr" + +if [ -z "$BIN" ]; then + echo "building..." + cargo build --release -p vcad-render --bin vcad-render --example psnr >/dev/null 2>&1 + BIN="$ROOT/target/release/vcad-render" +fi + +mkdir -p "$REF_DIR" "$OUT_DIR" +LOG="$OUT_DIR/render.log" +: >"$LOG" + +# scene name : path : size +SCENES=( + "rose-pro:hardware/rose-pro/rose-pro.loon:$SIZE" + "plate:examples/parametric-plate.vcad:400" +) + +status=0 +for entry in "${SCENES[@]}"; do + name="${entry%%:*}" + rest="${entry#*:}" + path="${rest%%:*}" + size="${rest#*:}" + + # References come in two flavours. The un-denoised one is the honest + # ground truth for the integrator; the denoised one is what a user + # actually looks at, and is the flavour the gate scores, because a + # sampling change that the denoiser papers over is not a regression the + # user can see. + for variant in raw denoised; do + ref="$REF_DIR/$name-$variant.png" + # `dn=()` plus `set -u` is an unbound expansion on bash 3.2 (macOS), so + # keep the array non-empty by carrying --seed inside it. + if [ "$variant" = raw ]; then + dn=(--no-denoise --seed "$SEED") + else + dn=(--seed "$SEED") + fi + + if [ $REGEN -eq 1 ] || [ ! -f "$ref" ]; then + echo "reference: $name/$variant @ ${REF_SPP}spp ${size}px" + # --no-adaptive: a reference is only ground truth if every pixel + # actually received the full budget. Adaptive sampling would stop the + # easy pixels early and quietly make the reference the thing under test. + "$BIN" "$path" --photoreal --spp "$REF_SPP" --size "$size" \ + --no-adaptive "${dn[@]}" -o "$ref" >>"$LOG" 2>&1 + fi + + cand="$OUT_DIR/$name-$variant.png" + # shellcheck disable=SC2086 + "$BIN" "$path" --photoreal --spp "$CAND_SPP" --size "$size" \ + "${dn[@]}" $EXTRA_ARGS -o "$cand" >>"$LOG" 2>&1 + + # The raw film at CAND_SPP is pure Monte Carlo noise measured against a + # 1024spp reference; report it for information, gate only on `denoised`. + if [ "$variant" = denoised ]; then + floor=(--min-psnr "$MIN_PSNR") + else + floor=(--min-psnr 0) + fi + # Capture rather than pipe: a pipeline's exit status is the *last* + # command's, which would swallow the gate's failure. + if report="$("$PSNR" "$ref" "$cand" "${floor[@]}")"; then + verdict=ok + else + verdict=FAIL + status=1 + fi + printf '%-12s %-9s %-32s %s\n' "$name" "$variant" \ + "$(echo "$report" | tr '\n' ' ')" "$verdict" + done +done + +exit $status