Fixing rendering issues, unifying rendering pipeline

.gitignore

@@ -14,3 +14,4 @@ tests/
 *.txt
 scenes/
 compile.sh
+output/

shared/Cargo.toml

@@ -4,6 +4,7 @@ version = "0.1.0"
 edition = "2024"
 
 [dependencies]
+anyhow = "1.0.100"
 bitflags = "2.10.0"
 half = "2.7.1"
 bytemuck = { version = "1.24.0", features = ["derive"] }

shared/src/core/aggregates.rs

@@ -1,7 +1,16 @@
-use crate::core::geometry::{Bounds3f, Ray, Vector3f};
+use crate::core::geometry::{Bounds3f, Point3f, Ray, Vector3f};
 use crate::core::primitive::{Primitive, PrimitiveTrait};
 use crate::core::shape::ShapeIntersection;
-use crate::{Float, Ptr, GVec, gvec};
+use crate::{gvec, Float, GVec, Ptr};
+
+#[repr(C)]
+#[derive(Debug, Clone, Copy, PartialEq, Eq)]
+pub enum SplitMethod {
+    SAH,
+    Hlbvh,
+    Middle,
+    EqualCounts,
+}
 
 #[repr(C)]
 #[derive(Default, Debug, Clone, Copy)]
@@ -18,7 +27,7 @@ pub struct LinearBVHNode {
 pub struct BVHAggregate {
     pub node_count: u32,
     pub max_prims_in_node: u32,
-    pub primitive_count: u32,
+    pub split_method: SplitMethod,
     pub primitives: GVec<Primitive>,
     pub nodes: GVec<LinearBVHNode>,
 }
@@ -26,11 +35,11 @@ pub struct BVHAggregate {
 impl BVHAggregate {
     pub fn empty() -> Self {
         Self {
-            max_prims_in_node: 0,
-            primitives: gvec(),
-            primitive_count: 0,
-            nodes: gvec(),
             node_count: 0,
+            max_prims_in_node: 0,
+            split_method: SplitMethod::SAH,
+            primitives: gvec(),
+            nodes: gvec(),
         }
     }
 
@@ -59,9 +68,10 @@ impl PrimitiveTrait for BVHAggregate {
             return None;
         }
 
-        let mut hit_t = t_max.unwrap_or(Float::INFINITY);
         let mut best_si: Option<ShapeIntersection> = None;
 
+        let mut hit_t = t_max.unwrap_or(Float::INFINITY);
+
         let inv_dir = Vector3f::new(1.0 / r.d.x(), 1.0 / r.d.y(), 1.0 / r.d.z());
         let dir_is_neg = [
             if inv_dir.x() < 0.0 { 1 } else { 0 },
@@ -69,53 +79,65 @@ impl PrimitiveTrait for BVHAggregate {
             if inv_dir.z() < 0.0 { 1 } else { 0 },
         ];
 
-        let mut stack = [0usize; 64];
+        let mut to_visit_offset = 0;
-        let mut stack_ptr = 0;
+        let mut current_node_index = 0;
-        let mut node_idx = 0usize;
+        let mut nodes_to_visit = [0usize; 64];
 
         loop {
-            let node = self.node(node_idx);
+            let node = &self.nodes[current_node_index];
 
+            // Check ray against BVH node bounds using the current closest hit_t
             if node
                 .bounds
                 .intersect_p(r.o, hit_t, inv_dir, &dir_is_neg)
-                .is_none()
+                .is_some()
             {
-                if stack_ptr == 0 {
+                if node.n_primitives > 0 {
-                    break;
+                    // Intersect ray with all primitives in this leaf
-                }
+                    for i in 0..node.n_primitives {
-                stack_ptr -= 1;
+                        let prim_idx = node.primitives_offset + i as usize;
-                node_idx = stack[stack_ptr];
+                        let prim = &self.primitives[prim_idx];
-                continue;
-            }
 
-            if node.n_primitives > 0 {
+                        if let Some(si) = prim.intersect(r, Some(hit_t)) {
-                // Leaf: test all primitives
+                            hit_t = si.t_hit();
-                for i in 0..node.n_primitives {
+                            best_si = Some(si);
-                    let prim_idx = node.primitives_offset + i as usize;
+                        }
-                    let prim = self.primitive(prim_idx);
+                    }
-                    if let Some(si) = prim.intersect(r, Some(hit_t)) {
+
-                        hit_t = si.t_hit();
+                    if to_visit_offset == 0 {
-                        best_si = Some(si);
+                        break;
+                    }
+                    to_visit_offset -= 1;
+                    current_node_index = nodes_to_visit[to_visit_offset];
+                } else {
+                    // Check the sign of the ray direction against the split axis
+                    if dir_is_neg[node.axis as usize] == 1 {
+                        // Ray is negative (Right -> Left).
+                        // Near child is Second Child (stored in primitives_offset).
+                        // Far child is First Child (current + 1).
+
+                        // Push Far
+                        nodes_to_visit[to_visit_offset] = current_node_index + 1;
+                        to_visit_offset += 1;
+
+                        // Visit Near immediately
+                        current_node_index = node.primitives_offset;
+                    } else {
+                        // Ray is positive (Left -> Right).
+                        // Push Far
+                        nodes_to_visit[to_visit_offset] = node.primitives_offset;
+                        to_visit_offset += 1;
+
+                        current_node_index += 1;
                     }
                 }
-
+            } else {
-                if stack_ptr == 0 {
+                // The ray missed the AABB of this node. Pop stack to try the next node.
+                if to_visit_offset == 0 {
                     break;
                 }
-                stack_ptr -= 1;
+                to_visit_offset -= 1;
-                node_idx = stack[stack_ptr];
+                current_node_index = nodes_to_visit[to_visit_offset];
-            } else {
-                // Interior: push far, visit near
-                if dir_is_neg[node.axis as usize] == 1 {
-                    stack[stack_ptr] = node_idx + 1;
-                    stack_ptr += 1;
-                    node_idx = node.primitives_offset; // second child
-                } else {
-                    stack[stack_ptr] = node.primitives_offset;
-                    stack_ptr += 1;
-                    node_idx += 1; // first child
-                }
             }
         }
 
@@ -123,7 +145,7 @@ impl PrimitiveTrait for BVHAggregate {
     }
 
     fn intersect_p(&self, r: &Ray, t_max: Option<Float>) -> bool {
-        if self.nodes.is_empty() || self.node_count == 0 {
+        if self.nodes.is_empty() {
             return false;
         }
 
@@ -136,53 +158,57 @@ impl PrimitiveTrait for BVHAggregate {
             if inv_dir.z() < 0.0 { 1 } else { 0 },
         ];
 
-        let mut stack = [0usize; 64];
+        let mut to_visit_offset = 0;
-        let mut stack_ptr = 0;
+        let mut current_node_index = 0;
-        let mut node_idx = 0usize;
+        let mut nodes_to_visit = [0usize; 64];
 
         loop {
-            let node = self.node(node_idx);
+            let node = &self.nodes[current_node_index];
 
+            // Check AABB
             if node
                 .bounds
                 .intersect_p(r.o, t_max, inv_dir, &dir_is_neg)
-                .is_none()
+                .is_some()
             {
-                if stack_ptr == 0 {
+                if node.n_primitives > 0 {
-                    break;
+                    for i in 0..node.n_primitives {
-                }
+                        let prim_idx = node.primitives_offset + i as usize;
-                stack_ptr -= 1;
+                        let prim = &self.primitives[prim_idx];
-                node_idx = stack[stack_ptr];
-                continue;
-            }
 
-            if node.n_primitives > 0 {
+                        if prim.intersect_p(r, Some(t_max)) {
-                for i in 0..node.n_primitives {
+                            return true;
-                    let prim_idx = node.primitives_offset + i as usize;
+                        }
-                    let prim = self.primitive(prim_idx);
+                    }
-                    if prim.intersect_p(r, Some(t_max)) {
+
-                        return true;
+                    // No intersection in this leaf, try next node in stack
+                    if to_visit_offset == 0 {
+                        break;
+                    }
+                    to_visit_offset -= 1;
+                    current_node_index = nodes_to_visit[to_visit_offset];
+                } else {
+                    // Standard front-to-back traversal order helps find an occlusion
+                    // closer to the origin faster, potentially saving work.
+
+                    if dir_is_neg[node.axis as usize] == 1 {
+                        nodes_to_visit[to_visit_offset] = current_node_index + 1;
+                        to_visit_offset += 1;
+                        current_node_index = node.primitives_offset;
+                    } else {
+                        nodes_to_visit[to_visit_offset] = node.primitives_offset;
+                        to_visit_offset += 1;
+                        current_node_index += 1;
                     }
                 }
-
+            } else {
-                if stack_ptr == 0 {
+                if to_visit_offset == 0 {
                     break;
                 }
-                stack_ptr -= 1;
+                to_visit_offset -= 1;
-                node_idx = stack[stack_ptr];
+                current_node_index = nodes_to_visit[to_visit_offset];
-            } else {
-                if dir_is_neg[node.axis as usize] == 1 {
-                    stack[stack_ptr] = node_idx + 1;
-                    stack_ptr += 1;
-                    node_idx = node.primitives_offset;
-                } else {
-                    stack[stack_ptr] = node.primitives_offset;
-                    stack_ptr += 1;
-                    node_idx += 1;
-                }
             }
         }
-
         false
     }
 }

shared/src/core/bssrdf.rs

@@ -3,11 +3,11 @@ use crate::core::bsdf::BSDF;
 use crate::core::geometry::{Frame, Normal3f, Point2f, Point3f, Point3fi, Vector3f};
 use crate::core::interaction::{InteractionBase, ShadingGeom, SurfaceInteraction};
 use crate::core::shape::Shape;
-use crate::spectra::{N_SPECTRUM_SAMPLES, SampledSpectrum};
+use crate::spectra::{SampledSpectrum, N_SPECTRUM_SAMPLES};
-use crate::utils::Ptr;
 use crate::utils::math::{catmull_rom_weights, square};
 use crate::utils::sampling::sample_catmull_rom_2d;
-use crate::{Float, PI};
+use crate::utils::Ptr;
+use crate::{gvec_with_capacity, Float, GVec, PI};
 use enum_dispatch::enum_dispatch;
 use num_traits::Float as NumFloat;
 
@@ -92,41 +92,59 @@ impl From<&SubsurfaceInteraction> for SurfaceInteraction {
 }
 
 #[repr(C)]
-#[derive(Clone, Copy, Debug)]
+#[derive(Clone, Debug)]
 pub struct BSSRDFTable {
     pub n_rho: u32,
     pub n_radius: u32,
-    pub rho_samples: Ptr<Float>,
+    pub rho_samples: GVec<Float>,
-    pub radius_samples: Ptr<Float>,
+    pub radius_samples: GVec<Float>,
-    pub profile: Ptr<Float>,
+    pub profile: GVec<Float>,
-    pub rho_eff: Ptr<Float>,
+    pub rho_eff: GVec<Float>,
-    pub profile_cdf: Ptr<Float>,
+    pub profile_cdf: GVec<Float>,
 }
 
 impl BSSRDFTable {
+    pub fn new(n_rho: usize, n_radius: usize) -> Self {
+        let rho_samples: GVec<Float> = gvec_with_capacity(n_rho);
+        let radius_samples: GVec<Float> = gvec_with_capacity(n_radius);
+        let profile: GVec<Float> = gvec_with_capacity(n_radius * n_rho);
+        let rho_eff: GVec<Float> = gvec_with_capacity(n_rho);
+        let profile_cdf: GVec<Float> = gvec_with_capacity(n_radius * n_rho);
+        Self {
+            n_rho: n_rho.try_into().unwrap(),
+            n_radius: n_radius.try_into().unwrap(),
+            rho_samples,
+            radius_samples,
+            profile,
+            rho_eff,
+            profile_cdf,
+
+        }
+    }
+
     pub fn get_rho(&self) -> &[Float] {
-        unsafe { core::slice::from_raw_parts(self.rho_samples.as_ref(), self.n_rho as usize) }
+        &self.rho_samples
     }
 
     pub fn get_radius(&self) -> &[Float] {
-        unsafe { core::slice::from_raw_parts(self.radius_samples.as_ref(), self.n_radius as usize) }
+        &self.radius_samples
     }
 
     pub fn get_profile(&self) -> &[Float] {
-        let n_profile = (self.n_rho * self.n_radius) as usize;
+        // let n_profile = (self.n_rho * self.n_radius) as usize;
-        unsafe { core::slice::from_raw_parts(self.profile.as_ref(), n_profile) }
+        &self.profile
     }
 
     pub fn get_cdf(&self) -> &[Float] {
-        let n_profile = (self.n_rho * self.n_radius) as usize;
+        // let n_profile = (self.n_rho * self.n_radius) as usize;
-        unsafe { core::slice::from_raw_parts(self.profile_cdf.as_ref(), n_profile) }
+        &self.profile_cdf
     }
 
     pub fn eval_profile(&self, rho_index: u32, radius_index: u32) -> Float {
         debug_assert!(rho_index < self.n_rho);
         debug_assert!(radius_index < self.n_radius);
         let idx = (rho_index * self.n_radius + radius_index) as usize;
-        unsafe { *self.profile.add(idx) }
+        unsafe { *self.profile.as_ptr().add(idx) }
     }
 }
 

shared/src/lights/diffuse.rs

@@ -164,7 +164,7 @@ impl LightTrait for DiffuseAreaLight {
 
     #[cfg(not(target_os = "cuda"))]
     fn preprocess(&mut self, _scene_bounds: &Bounds3f) {
-        unimplemented!()
+        return
     }
 
     #[cfg(not(target_os = "cuda"))]

shared/src/materials/complex.rs

@@ -149,7 +149,7 @@ pub struct SubsurfaceMaterial {
     pub u_roughness: Ptr<GPUFloatTexture>,
     pub v_roughness: Ptr<GPUFloatTexture>,
     pub remap_roughness: bool,
-    pub table: BSSRDFTable,
+    pub table: Ptr<BSSRDFTable>,
 }
 
 impl MaterialTrait for SubsurfaceMaterial {

shared/src/utils/transform.rs

@@ -2,7 +2,7 @@ use core::iter::{Product, Sum};
 use core::ops::{Add, Div, Index, IndexMut, Mul};
 use num_traits::Float as NumFloat;
 
-use super::math::{SquareMatrix, radians, safe_acos};
+use super::math::{radians, safe_acos, SquareMatrix};
 use super::quaternion::Quaternion;
 use crate::core::color::{RGB, XYZ};
 use crate::core::geometry::{
@@ -12,8 +12,9 @@ use crate::core::geometry::{
 use crate::core::interaction::{
     Interaction, InteractionBase, InteractionTrait, MediumInteraction, SurfaceInteraction,
 };
+use anyhow::{bail, Context, Result};
 use crate::utils::gpu_array_from_fn;
-use crate::{Float, gamma};
+use crate::{gamma, Float};
 
 #[repr(C)]
 #[derive(Debug, Copy, Clone)]
@@ -173,27 +174,27 @@ impl TransformGeneric<Float> {
         let x = p.x();
         let y = p.y();
         let z = p.z();
-        let xp = self.m[0][0] * x + self.m[0][1] * y + self.m[0][2] * z;
+        let xp = self.m[0][0] * x + self.m[0][1] * y + self.m[0][2] * z + self.m[0][3];
-        let yp = self.m[1][0] * x + self.m[1][1] * y + self.m[1][2] * z;
+        let yp = self.m[1][0] * x + self.m[1][1] * y + self.m[1][2] * z + self.m[1][3];
-        let zp = self.m[2][0] * x + self.m[2][1] * y + self.m[2][2] * z;
+        let zp = self.m[2][0] * x + self.m[2][1] * y + self.m[2][2] * z + self.m[2][3];
-        let wp = self.m[3][0] * x + self.m[3][1] * y + self.m[3][2] * z;
+        let wp = self.m[3][0] * x + self.m[3][1] * y + self.m[3][2] * z + self.m[3][3];
         let mut p_error = Vector3f::default();
 
         if pi.is_exact() {
             p_error[0] = gamma(3)
                 * ((self.m[0][0] * x).abs()
                     + (self.m[0][1] * y).abs()
-                    + (self.m[0][2] + z).abs()
+                    + (self.m[0][2] * z).abs()
                     + self.m[0][3].abs());
             p_error[1] = gamma(3)
                 * ((self.m[1][0] * x).abs()
                     + (self.m[1][1] * y).abs()
-                    + (self.m[1][2] + z).abs()
+                    + (self.m[1][2] * z).abs()
                     + self.m[1][3].abs());
             p_error[2] = gamma(3)
                 * ((self.m[2][0] * x).abs()
                     + (self.m[2][1] * y).abs()
-                    + (self.m[2][2] + z).abs()
+                    + (self.m[2][2] * z).abs()
                     + self.m[2][3].abs());
         }
 
@@ -210,7 +211,6 @@ impl TransformGeneric<Float> {
         let y = v.y();
         let z = v.z();
 
-        // Transform the midpoint of the vector interval
         let xp = self.m[0][0] * x + self.m[0][1] * y + self.m[0][2] * z;
         let yp = self.m[1][0] * x + self.m[1][1] * y + self.m[1][2] * z;
         let zp = self.m[2][0] * x + self.m[2][1] * y + self.m[2][2] * z;
@@ -218,7 +218,6 @@ impl TransformGeneric<Float> {
 
         let mut v_error = Vector3f::default();
 
-        // Propagate the error, ignoring the translational part of the matrix
         if vi.is_exact() {
             v_error[0] = gamma(3)
                 * ((self.m[0][0] * x).abs() + (self.m[0][1] * y).abs() + (self.m[0][2] * z).abs());
@@ -226,6 +225,32 @@ impl TransformGeneric<Float> {
                 * ((self.m[1][0] * x).abs() + (self.m[1][1] * y).abs() + (self.m[1][2] * z).abs());
             v_error[2] = gamma(3)
                 * ((self.m[2][0] * x).abs() + (self.m[2][1] * y).abs() + (self.m[2][2] * z).abs());
+        } else {
+            let vin_error = vi.error();
+            v_error[0] = (gamma(3) + 1.)
+                * ((self.m[0][0] * x).abs() * vin_error.x()
+                    + (self.m[0][1] * y).abs() * vin_error.y()
+                    + (self.m[0][2] * z).abs() * vin_error.z())
+                + gamma(3)
+                    * ((self.m[0][0] * x).abs()
+                        + (self.m[0][1] * y).abs()
+                        + (self.m[0][2] * z).abs());
+            v_error[1] = (gamma(3) + 1.)
+                * ((self.m[1][0] * x).abs() * vin_error.x()
+                    + (self.m[1][1] * y).abs() * vin_error.y()
+                    + (self.m[1][2] * z).abs() * vin_error.z())
+                + gamma(3)
+                    * ((self.m[1][0] * x).abs()
+                        + (self.m[1][1] * y).abs()
+                        + (self.m[1][2] * z).abs());
+            v_error[2] = (gamma(3) + 1.)
+                * ((self.m[2][0] * x).abs() * vin_error.x()
+                    + (self.m[2][1] * y).abs() * vin_error.y()
+                    + (self.m[2][2] * z).abs() * vin_error.z())
+                + gamma(3)
+                    * ((self.m[2][0] * x).abs()
+                        + (self.m[2][1] * y).abs()
+                        + (self.m[2][2] * z).abs());
         }
 
         if wp == 1. {
@@ -2094,7 +2119,7 @@ pub fn look_at(
     pos: impl Into<Point3f>,
     look: impl Into<Point3f>,
     up: impl Into<Point3f>,
-) -> Option<TransformGeneric<Float>> {
+) -> Result<TransformGeneric<Float>> {
     let mut world_from_camera: SquareMatrix<Float, 4> = SquareMatrix::default();
     // Initialize fourth column of viewing matrix
     let pos: Point3f = pos.into();
@@ -2108,7 +2133,7 @@ pub fn look_at(
     // Initialize first three columns of viewing matrix
     let dir = (look - pos).normalize();
     if Vector3f::from(up).normalize().cross(dir).norm() == 0. {
-        panic!(
+        bail!(
             "LookAt: \"up\" vector ({}, {}, {}) and viewing direction ({}, {}, {}) passed to LookAt are pointing in the same direction.",
             up.x(),
             up.y(),
@@ -2133,6 +2158,6 @@ pub fn look_at(
     world_from_camera[2][2] = dir.z();
     world_from_camera[3][2] = 0.;
 
-    let camera_from_world = world_from_camera.inverse()?;
+    let camera_from_world = world_from_camera.inverse().context("Failed to inverse viewing matrix")?;
-    Some(TransformGeneric::new(camera_from_world, world_from_camera))
+    Ok(TransformGeneric::new(camera_from_world, world_from_camera))
 }

src/core/aggregates.rs

@@ -1,24 +1,15 @@
 use crate::Arena;
 use rayon::prelude::*;
-use shared::core::aggregates::{BVHAggregate as DeviceBVHAggregate, LinearBVHNode};
+use shared::core::aggregates::{BVHAggregate, LinearBVHNode, SplitMethod};
 use shared::core::geometry::{Bounds3f, Point3f, Ray, Vector3f};
 use shared::core::primitive::{Primitive, PrimitiveTrait};
 use shared::core::shape::ShapeIntersection;
 use shared::utils::math::encode_morton_3;
 use shared::utils::{find_interval, partition_slice};
-use shared::{gvec_from_slice, Float};
+use shared::{gvec, gvec_from_slice, Float};
 use std::cmp::Ordering;
 use std::sync::atomic::{AtomicUsize, Ordering as AtomicOrdering};
 
-#[repr(C)]
-#[derive(Debug, Clone, Copy, PartialEq, Eq)]
-pub enum SplitMethod {
-    SAH,
-    Hlbvh,
-    Middle,
-    EqualCounts,
-}
-
 #[repr(C)]
 #[derive(Debug, Default, Clone, Copy, PartialEq)]
 struct BVHSplitBucket {
@@ -112,51 +103,27 @@ impl BVHBuildNode {
     }
 }
 
-pub struct SharedPrimitiveBuffer<'a, P> {
+pub trait CreateBVH {
-    ptr: *mut P,
+    fn new(primitives: Vec<Primitive>, max_prims_in_node: usize, split_method: SplitMethod) -> Self;
-    pub offset: &'a AtomicUsize,
+    fn build_hlbvh(
-    _marker: std::marker::PhantomData<&'a mut [P]>,
+        bvh_primitives: &[BVHPrimitiveInfo],
+        total_nodes: &AtomicUsize,
+        _original_primitives: &[Primitive],
+        max_prims_in_node: usize,
+    ) -> Box<BVHBuildNode>;
+    fn emit_lbvh(
+        bvh_primitives: &[BVHPrimitiveInfo],
+        morton_prims: &[MortonPrimitive],
+        total_nodes: &mut usize,
+        bit_index: i32,
+        max_prims_in_node: usize,
+    ) -> Box<BVHBuildNode>;
+    fn build_upper_sah(nodes: &mut [BVHBuildNode], total_nodes: &AtomicUsize) -> Box<BVHBuildNode>;
 }
 
-unsafe impl<'a, P> Sync for SharedPrimitiveBuffer<'a, P> {}
+impl CreateBVH for BVHAggregate {
-unsafe impl<'a, P> Send for SharedPrimitiveBuffer<'a, P> {}
+    fn new(
-
+        mut primitives: Vec<Primitive>,
-impl<'a, P> SharedPrimitiveBuffer<'a, P> {
-    pub fn new(slice: &'a mut [P], offset: &'a AtomicUsize) -> Self {
-        Self {
-            ptr: slice.as_mut_ptr(),
-            offset,
-            _marker: std::marker::PhantomData,
-        }
-    }
-
-    pub fn append(&self, primitives: &[P], indices: &[BVHPrimitiveInfo]) -> usize
-    where
-        P: Clone,
-    {
-        let count = indices.len();
-        let start_index = self.offset.fetch_add(count, AtomicOrdering::Relaxed);
-
-        unsafe {
-            for (i, info) in indices.iter().enumerate() {
-                let target_ptr = self.ptr.add(start_index + i);
-                std::ptr::write(target_ptr, primitives[info.primitive_number].clone());
-            }
-        }
-        start_index
-    }
-}
-
-pub struct BVHAggregate<P: PrimitiveTrait + Clone + Send + Sync> {
-    pub max_prims_in_node: usize,
-    pub primitives: Vec<P>,
-    pub split_method: SplitMethod,
-    pub nodes: Vec<LinearBVHNode>,
-}
-
-impl<P: PrimitiveTrait + Clone + Send + Sync> BVHAggregate<P> {
-    pub fn new(
-        mut primitives: Vec<P>,
         max_prims_in_node: usize,
         split_method: SplitMethod,
     ) -> Self {
@@ -164,10 +131,11 @@ impl<P: PrimitiveTrait + Clone + Send + Sync> BVHAggregate<P> {
 
         if primitives.is_empty() {
             return Self {
-                max_prims_in_node,
+                max_prims_in_node: max_prims_in_node.try_into().unwrap(),
-                primitives,
+                node_count: 0,
+                primitives: gvec_from_slice(&primitives),
                 split_method,
-                nodes: Vec::new(),
+                nodes: gvec(),
             };
         }
 
@@ -177,7 +145,7 @@ impl<P: PrimitiveTrait + Clone + Send + Sync> BVHAggregate<P> {
             .map(|(i, p)| BVHPrimitiveInfo::new(i, p.bounds()))
             .collect();
 
-        let total_nodes_count: usize;
+        let node_count: usize;
         let root: Box<BVHBuildNode>;
 
         match split_method {
@@ -189,124 +157,46 @@ impl<P: PrimitiveTrait + Clone + Send + Sync> BVHAggregate<P> {
                     &primitives,
                     max_prims_in_node,
                 );
-                total_nodes_count = nodes_counter.load(AtomicOrdering::Relaxed);
+                node_count = nodes_counter.load(AtomicOrdering::Relaxed);
             }
             _ => {
                 let nodes_counter = AtomicUsize::new(0);
-                root = Self::build_recursive(
+                root = build_recursive(
                     &mut primitive_info,
                     &nodes_counter,
                     &primitives,
                     max_prims_in_node,
                     split_method,
                 );
-                total_nodes_count = nodes_counter.load(AtomicOrdering::Relaxed);
+                node_count = nodes_counter.load(AtomicOrdering::Relaxed);
             }
         };
 
         // Walk the tree and collect primitive indices in the exact order
         // the linear layout will visit them (left-to-right, depth-first)
-        let mut leaf_order = Vec::with_capacity(primitives.len());
+        let mut leaf_vec = Vec::with_capacity(primitives.len());
-        Self::leaf_order(&root, &mut leaf_order);
+        leaf_order(&root, &mut leaf_vec);
-        Self::reorder(&mut primitives, &leaf_order);
+        reorder(&mut primitives, &leaf_vec);
-        drop(leaf_order);
+        drop(leaf_vec);
 
-        let mut nodes = vec![LinearBVHNode::default(); total_nodes_count];
+        let mut nodes = vec![LinearBVHNode::default(); node_count];
         let mut offset = 0;
         let mut prim_offset = 0;
-        Self::flatten(&root, &mut nodes, &mut offset, &mut prim_offset);
+        flatten(&root, &mut nodes, &mut offset, &mut prim_offset);
 
         Self {
-            max_prims_in_node,
+            node_count: node_count.try_into().unwrap(),
-            primitives,
+            max_prims_in_node: max_prims_in_node.try_into().unwrap(),
             split_method,
-            nodes,
+            primitives: gvec_from_slice(&primitives),
+            nodes: gvec_from_slice(&nodes),
         }
     }
 
-    fn reorder(primitives: &mut [P], order: &[usize]) {
+    fn build_hlbvh(
-        let n = primitives.len();
-        assert_eq!(n, order.len());
-
-        let mut done = vec![false; n];
-        for i in 0..n {
-            if done[i] || order[i] == i {
-                done[i] = true;
-                continue;
-            }
-
-            let mut prev = i;
-            let mut curr = order[i];
-            while curr != i {
-                primitives.swap(prev, curr);
-                done[prev] = true;
-                prev = curr;
-                curr = order[prev];
-            }
-            done[prev] = true;
-        }
-    }
-
-    fn leaf_order(node: &BVHBuildNode, out: &mut Vec<usize>) {
-        match node {
-            BVHBuildNode::Leaf {
-                primitive_indices, ..
-            } => {
-                out.extend_from_slice(primitive_indices);
-            }
-            BVHBuildNode::Interior { children, .. } => {
-                Self::leaf_order(&children[0], out);
-                Self::leaf_order(&children[1], out);
-            }
-        }
-    }
-
-    fn flatten(
-        node: &BVHBuildNode,
-        nodes: &mut [LinearBVHNode],
-        offset: &mut usize,
-        prim_offset: &mut usize,
-    ) -> usize {
-        let local_offset = *offset;
-        *offset += 1;
-
-        match node {
-            BVHBuildNode::Leaf {
-                n_primitives,
-                bounds,
-                ..
-            } => {
-                let n = *n_primitives;
-                let linear_node = &mut nodes[local_offset];
-                linear_node.bounds = *bounds;
-                linear_node.n_primitives = n as u16;
-                linear_node.primitives_offset = *prim_offset;
-                linear_node.axis = 0; // Irrelevant for leaves
-                *prim_offset += n;
-            }
-
-            BVHBuildNode::Interior {
-                split_axis,
-                children,
-                bounds,
-            } => {
-                nodes[local_offset].bounds = *bounds;
-                nodes[local_offset].axis = *split_axis;
-                nodes[local_offset].n_primitives = 0;
-
-                Self::flatten(&children[0], nodes, offset, prim_offset);
-                let second_child_offset = Self::flatten(&children[1], nodes, offset, prim_offset);
-                nodes[local_offset].primitives_offset = second_child_offset;
-            }
-        }
-
-        local_offset
-    }
-
-    pub fn build_hlbvh(
         bvh_primitives: &[BVHPrimitiveInfo],
         total_nodes: &AtomicUsize,
-        _original_primitives: &[P],
+        _original_primitives: &[Primitive],
         max_prims_in_node: usize,
     ) -> Box<BVHBuildNode> {
         let bounds = bvh_primitives
@@ -586,327 +476,252 @@ impl<P: PrimitiveTrait + Clone + Send + Sync> BVHAggregate<P> {
             ))
         }
     }
+}
 
-    fn build_recursive(
+fn build_recursive(
-        bvh_primitives: &mut [BVHPrimitiveInfo],
+    bvh_primitives: &mut [BVHPrimitiveInfo],
-        total_nodes: &AtomicUsize,
+    total_nodes: &AtomicUsize,
-        original_primitives: &[P],
+    original_primitives: &[Primitive],
-        max_prims_in_node: usize,
+    max_prims_in_node: usize,
-        split_method: SplitMethod,
+    split_method: SplitMethod,
-    ) -> Box<BVHBuildNode> {
+) -> Box<BVHBuildNode> {
-        total_nodes.fetch_add(1, AtomicOrdering::Relaxed);
+    total_nodes.fetch_add(1, AtomicOrdering::Relaxed);
-        let bounds = bvh_primitives
+    let bounds = bvh_primitives
-            .iter()
+        .iter()
-            .fold(Bounds3f::default(), |b, p| b.union(p.bounds));
+        .fold(Bounds3f::default(), |b, p| b.union(p.bounds));
 
-        let n_primitives = bvh_primitives.len();
+    let n_primitives = bvh_primitives.len();
-        if bounds.surface_area() == 0.0 || n_primitives == 1 || n_primitives <= max_prims_in_node {
+    if bounds.surface_area() == 0.0 || n_primitives == 1 || n_primitives <= max_prims_in_node {
-            let indices: Vec<usize> = bvh_primitives.iter().map(|p| p.primitive_number).collect();
+        let indices: Vec<usize> = bvh_primitives.iter().map(|p| p.primitive_number).collect();
-            return Box::new(BVHBuildNode::new_leaf(n_primitives, bounds, indices));
+        return Box::new(BVHBuildNode::new_leaf(n_primitives, bounds, indices));
-        }
+    }
-        let centroid_bounds = bvh_primitives.iter().fold(Bounds3f::default(), |b, p| {
+    let centroid_bounds = bvh_primitives.iter().fold(Bounds3f::default(), |b, p| {
-            b.union_point(p.bounds.centroid())
+        b.union_point(p.bounds.centroid())
-        });
+    });
 
-        let dim = centroid_bounds.max_dimension();
+    let dim = centroid_bounds.max_dimension();
-        if centroid_bounds.p_max[dim] == centroid_bounds.p_min[dim] {
+    if centroid_bounds.p_max[dim] == centroid_bounds.p_min[dim] {
-            let indices: Vec<usize> = bvh_primitives.iter().map(|p| p.primitive_number).collect();
+        let indices: Vec<usize> = bvh_primitives.iter().map(|p| p.primitive_number).collect();
-            return Box::new(BVHBuildNode::new_leaf(n_primitives, bounds, indices));
+        return Box::new(BVHBuildNode::new_leaf(n_primitives, bounds, indices));
-        }
+    }
 
-        let mut mid: usize;
+    let mut mid: usize;
-        match split_method {
+    match split_method {
-            SplitMethod::Middle => {
+        SplitMethod::Middle => {
-                let pmid = (centroid_bounds.p_min[dim] + centroid_bounds.p_max[dim]) / 2.;
+            let pmid = (centroid_bounds.p_min[dim] + centroid_bounds.p_max[dim]) / 2.;
-                mid = partition_slice(bvh_primitives, |p| p.centroid[dim] < pmid);
+            mid = partition_slice(bvh_primitives, |p| p.centroid[dim] < pmid);
 
-                if mid != 0 && mid != n_primitives {
+            if mid != 0 && mid != n_primitives {
-                } else {
+            } else {
-                    mid = n_primitives / 2;
-                    bvh_primitives.select_nth_unstable_by(mid, |a, b| {
-                        a.centroid[dim].partial_cmp(&b.centroid[dim]).unwrap()
-                    });
-                }
-            }
-            SplitMethod::EqualCounts => {
                 mid = n_primitives / 2;
                 bvh_primitives.select_nth_unstable_by(mid, |a, b| {
                     a.centroid[dim].partial_cmp(&b.centroid[dim]).unwrap()
                 });
             }
-            SplitMethod::SAH | _ => {
+        }
-                if n_primitives < 2 {
+        SplitMethod::EqualCounts => {
-                    mid = n_primitives / 2;
+            mid = n_primitives / 2;
-                    bvh_primitives.select_nth_unstable_by(mid, |a, b| {
+            bvh_primitives.select_nth_unstable_by(mid, |a, b| {
-                        a.centroid[dim]
+                a.centroid[dim].partial_cmp(&b.centroid[dim]).unwrap()
-                            .partial_cmp(&b.centroid[dim])
+            });
-                            .unwrap_or(Ordering::Equal)
+        }
-                    });
+        SplitMethod::SAH | _ => {
-                } else {
+            if n_primitives < 2 {
-                    const N_BUCKETS: usize = 12;
+                mid = n_primitives / 2;
-                    let mut buckets = [BVHSplitBucket::default(); N_BUCKETS];
+                bvh_primitives.select_nth_unstable_by(mid, |a, b| {
-                    for prim in bvh_primitives.iter() {
+                    a.centroid[dim]
-                        let mut b = (N_BUCKETS as Float
+                        .partial_cmp(&b.centroid[dim])
-                            * centroid_bounds.offset(&prim.centroid)[dim])
+                        .unwrap_or(Ordering::Equal)
+                });
+            } else {
+                const N_BUCKETS: usize = 12;
+                let mut buckets = [BVHSplitBucket::default(); N_BUCKETS];
+                for prim in bvh_primitives.iter() {
+                    let mut b =
+                        (N_BUCKETS as Float * centroid_bounds.offset(&prim.centroid)[dim]) as usize;
+                    if b == N_BUCKETS {
+                        b = N_BUCKETS - 1;
+                    }
+                    buckets[b].count += 1;
+                    buckets[b].bounds = buckets[b].bounds.union(prim.bounds);
+                }
+                // Compute costs for splitting after each bucket>
+                const N_SPLITS: usize = N_BUCKETS - 1;
+                let mut costs = [0.0 as Float; N_SPLITS];
+                let mut count_below = 0;
+                let mut bound_below = Bounds3f::default();
+                for i in 0..N_SPLITS {
+                    bound_below = bound_below.union(buckets[i].bounds);
+                    count_below += buckets[i].count;
+                    costs[i] += count_below as Float * bound_below.surface_area();
+                }
+                // Finish initializing costs using a backward scan over splits
+                let mut count_above = 0;
+                let mut bound_above = Bounds3f::default();
+                for i in (0..N_SPLITS).rev() {
+                    bound_above = bound_above.union(buckets[i + 1].bounds);
+                    count_above += buckets[i + 1].count;
+                    costs[i] += count_above as Float * bound_above.surface_area();
+                }
+
+                // Find bucket to split at that minimizes SAH metric>
+                let mut min_cost = Float::INFINITY;
+                let mut min_cost_split_bucket = 0;
+                for (i, &cost) in costs.iter().enumerate().take(N_SPLITS) {
+                    if cost < min_cost {
+                        min_cost = cost;
+                        min_cost_split_bucket = i;
+                    }
+                }
+
+                // Compute leaf cost and SAH split cost for chosen split
+                let leaf_cost = n_primitives as Float;
+                min_cost = 0.5 + min_cost / bounds.surface_area();
+
+                // Either create leaf or split primitives at selected SAH bucket>
+                if n_primitives > max_prims_in_node || min_cost < leaf_cost {
+                    mid = partition_slice(bvh_primitives, |bp| {
+                        let mut b = (N_BUCKETS as Float * centroid_bounds.offset(&bp.centroid)[dim])
                             as usize;
                         if b == N_BUCKETS {
                             b = N_BUCKETS - 1;
                         }
-                        buckets[b].count += 1;
+                        b <= min_cost_split_bucket
-                        buckets[b].bounds = buckets[b].bounds.union(prim.bounds);
+                    });
-                    }
+                    if mid == 0 || mid == n_primitives {
-                    // Compute costs for splitting after each bucket>
+                        mid = n_primitives / 2;
-                    const N_SPLITS: usize = N_BUCKETS - 1;
+                        bvh_primitives.select_nth_unstable_by(mid, |a, b| {
-                    let mut costs = [0.0 as Float; N_SPLITS];
+                            a.centroid[dim]
-                    let mut count_below = 0;
+                                .partial_cmp(&b.centroid[dim])
-                    let mut bound_below = Bounds3f::default();
+                                .unwrap_or(Ordering::Equal)
-                    for i in 0..N_SPLITS {
-                        bound_below = bound_below.union(buckets[i].bounds);
-                        count_below += buckets[i].count;
-                        costs[i] += count_below as Float * bound_below.surface_area();
-                    }
-                    // Finish initializing costs using a backward scan over splits
-                    let mut count_above = 0;
-                    let mut bound_above = Bounds3f::default();
-                    for i in (0..N_SPLITS).rev() {
-                        bound_above = bound_above.union(buckets[i + 1].bounds);
-                        count_above += buckets[i + 1].count;
-                        costs[i] += count_above as Float * bound_above.surface_area();
-                    }
-
-                    // Find bucket to split at that minimizes SAH metric>
-                    let mut min_cost = Float::INFINITY;
-                    let mut min_cost_split_bucket = 0;
-                    for (i, &cost) in costs.iter().enumerate().take(N_SPLITS) {
-                        if cost < min_cost {
-                            min_cost = cost;
-                            min_cost_split_bucket = i;
-                        }
-                    }
-
-                    // Compute leaf cost and SAH split cost for chosen split
-                    let leaf_cost = n_primitives as Float;
-                    min_cost = 0.5 + min_cost / bounds.surface_area();
-
-                    // Either create leaf or split primitives at selected SAH bucket>
-                    if n_primitives > max_prims_in_node || min_cost < leaf_cost {
-                        mid = partition_slice(bvh_primitives, |bp| {
-                            let mut b = (N_BUCKETS as Float
-                                * centroid_bounds.offset(&bp.centroid)[dim])
-                                as usize;
-                            if b == N_BUCKETS {
-                                b = N_BUCKETS - 1;
-                            }
-                            b <= min_cost_split_bucket
                         });
-                        if mid == 0 || mid == n_primitives {
-                            mid = n_primitives / 2;
-                            bvh_primitives.select_nth_unstable_by(mid, |a, b| {
-                                a.centroid[dim]
-                                    .partial_cmp(&b.centroid[dim])
-                                    .unwrap_or(Ordering::Equal)
-                            });
-                        }
-                    } else {
-                        let indices: Vec<usize> =
-                            bvh_primitives.iter().map(|p| p.primitive_number).collect();
-                        return Box::new(BVHBuildNode::new_leaf(n_primitives, bounds, indices));
                     }
+                } else {
+                    let indices: Vec<usize> =
+                        bvh_primitives.iter().map(|p| p.primitive_number).collect();
+                    return Box::new(BVHBuildNode::new_leaf(n_primitives, bounds, indices));
                 }
             }
-        };
+        }
+    };
 
-        let (left_prims, right_prims) = bvh_primitives.split_at_mut(mid);
+    let (left_prims, right_prims) = bvh_primitives.split_at_mut(mid);
-        let build_leaf = |prims: &mut [BVHPrimitiveInfo]| -> Box<BVHBuildNode> {
+    let build_leaf = |prims: &mut [BVHPrimitiveInfo]| -> Box<BVHBuildNode> {
-            Self::build_recursive(
+        build_recursive(
-                prims,
+            prims,
-                total_nodes,
+            total_nodes,
-                original_primitives,
+            original_primitives,
-                max_prims_in_node,
+            max_prims_in_node,
-                split_method,
+            split_method,
-            )
+        )
-        };
+    };
 
-        let (child0, child1) = if n_primitives > 128 * 1024 {
+    let (child0, child1) = if n_primitives > 128 * 1024 {
-            rayon::join(|| build_leaf(left_prims), || build_leaf(right_prims))
+        rayon::join(|| build_leaf(left_prims), || build_leaf(right_prims))
-        } else {
+    } else {
-            (build_leaf(left_prims), build_leaf(right_prims))
+        (build_leaf(left_prims), build_leaf(right_prims))
-        };
+    };
 
-        let axis = dim as u8;
+    let axis = dim as u8;
-        Box::new(BVHBuildNode::new_interior(axis, child0, child1))
+    Box::new(BVHBuildNode::new_interior(axis, child0, child1))
+}
+
+fn flatten(
+    node: &BVHBuildNode,
+    nodes: &mut [LinearBVHNode],
+    offset: &mut usize,
+    prim_offset: &mut usize,
+) -> usize {
+    let local_offset = *offset;
+    *offset += 1;
+
+    match node {
+        BVHBuildNode::Leaf {
+            n_primitives,
+            bounds,
+            ..
+        } => {
+            let n = *n_primitives;
+            let linear_node = &mut nodes[local_offset];
+            linear_node.bounds = *bounds;
+            linear_node.n_primitives = n as u16;
+            linear_node.primitives_offset = *prim_offset;
+            linear_node.axis = 0; // Irrelevant for leaves
+            *prim_offset += n;
+        }
+
+        BVHBuildNode::Interior {
+            split_axis,
+            children,
+            bounds,
+        } => {
+            nodes[local_offset].bounds = *bounds;
+            nodes[local_offset].axis = *split_axis;
+            nodes[local_offset].n_primitives = 0;
+
+            flatten(&children[0], nodes, offset, prim_offset);
+            let second_child_offset = flatten(&children[1], nodes, offset, prim_offset);
+            nodes[local_offset].primitives_offset = second_child_offset;
+        }
     }
 
-    pub fn intersect(&self, r: &Ray, t_max: Option<Float>) -> Option<ShapeIntersection> {
+    local_offset
-        if self.nodes.is_empty() {
+}
-            return None;
+
+fn reorder(primitives: &mut [Primitive], order: &[usize]) {
+    let n = primitives.len();
+    assert_eq!(n, order.len());
+
+    let mut done = vec![false; n];
+    for i in 0..n {
+        if done[i] || order[i] == i {
+            done[i] = true;
+            continue;
         }
 
-        let mut best_si: Option<ShapeIntersection> = None;
+        let mut prev = i;
-
+        let mut curr = order[i];
-        let mut hit_t = t_max.unwrap_or(Float::INFINITY);
+        while curr != i {
-
+            primitives.swap(prev, curr);
-        let inv_dir = Vector3f::new(1.0 / r.d.x(), 1.0 / r.d.y(), 1.0 / r.d.z());
+            done[prev] = true;
-        let dir_is_neg = [
+            prev = curr;
-            if inv_dir.x() < 0.0 { 1 } else { 0 },
+            curr = order[prev];
-            if inv_dir.y() < 0.0 { 1 } else { 0 },
-            if inv_dir.z() < 0.0 { 1 } else { 0 },
-        ];
-
-        let mut to_visit_offset = 0;
-        let mut current_node_index = 0;
-        let mut nodes_to_visit = [0usize; 64];
-
-        loop {
-            let node = &self.nodes[current_node_index];
-
-            // Check ray against BVH node bounds using the current closest hit_t
-            if node
-                .bounds
-                .intersect_p(r.o, hit_t, inv_dir, &dir_is_neg)
-                .is_some()
-            {
-                if node.n_primitives > 0 {
-                    // Intersect ray with all primitives in this leaf
-                    for i in 0..node.n_primitives {
-                        let prim_idx = node.primitives_offset + i as usize;
-                        let prim = &self.primitives[prim_idx];
-
-                        if let Some(si) = prim.intersect(r, Some(hit_t)) {
-                            hit_t = si.t_hit();
-                            best_si = Some(si);
-                        }
-                    }
-
-                    if to_visit_offset == 0 {
-                        break;
-                    }
-                    to_visit_offset -= 1;
-                    current_node_index = nodes_to_visit[to_visit_offset];
-                } else {
-                    // Check the sign of the ray direction against the split axis
-                    if dir_is_neg[node.axis as usize] == 1 {
-                        // Ray is negative (Right -> Left).
-                        // Near child is Second Child (stored in primitives_offset).
-                        // Far child is First Child (current + 1).
-
-                        // Push Far
-                        nodes_to_visit[to_visit_offset] = current_node_index + 1;
-                        to_visit_offset += 1;
-
-                        // Visit Near immediately
-                        current_node_index = node.primitives_offset;
-                    } else {
-                        // Ray is positive (Left -> Right).
-                        // Push Far
-                        nodes_to_visit[to_visit_offset] = node.primitives_offset;
-                        to_visit_offset += 1;
-
-                        current_node_index += 1;
-                    }
-                }
-            } else {
-                // The ray missed the AABB of this node. Pop stack to try the next node.
-                if to_visit_offset == 0 {
-                    break;
-                }
-                to_visit_offset -= 1;
-                current_node_index = nodes_to_visit[to_visit_offset];
-            }
         }
-
+        done[prev] = true;
-        best_si
-    }
-
-    fn intersect_p(&self, r: &Ray, t_max: Option<Float>) -> bool {
-        if self.nodes.is_empty() {
-            return false;
-        }
-
-        let t_max = t_max.unwrap_or(Float::INFINITY);
-
-        let inv_dir = Vector3f::new(1.0 / r.d.x(), 1.0 / r.d.y(), 1.0 / r.d.z());
-        let dir_is_neg = [
-            if inv_dir.x() < 0.0 { 1 } else { 0 },
-            if inv_dir.y() < 0.0 { 1 } else { 0 },
-            if inv_dir.z() < 0.0 { 1 } else { 0 },
-        ];
-
-        let mut to_visit_offset = 0;
-        let mut current_node_index = 0;
-        let mut nodes_to_visit = [0usize; 64];
-
-        loop {
-            let node = &self.nodes[current_node_index];
-
-            // Check AABB
-            if node
-                .bounds
-                .intersect_p(r.o, t_max, inv_dir, &dir_is_neg)
-                .is_some()
-            {
-                if node.n_primitives > 0 {
-                    for i in 0..node.n_primitives {
-                        let prim_idx = node.primitives_offset + i as usize;
-                        let prim = &self.primitives[prim_idx];
-
-                        if prim.intersect_p(r, Some(t_max)) {
-                            return true;
-                        }
-                    }
-
-                    // No intersection in this leaf, try next node in stack
-                    if to_visit_offset == 0 {
-                        break;
-                    }
-                    to_visit_offset -= 1;
-                    current_node_index = nodes_to_visit[to_visit_offset];
-                } else {
-                    // Standard front-to-back traversal order helps find an occlusion
-                    // closer to the origin faster, potentially saving work.
-
-                    if dir_is_neg[node.axis as usize] == 1 {
-                        nodes_to_visit[to_visit_offset] = current_node_index + 1;
-                        to_visit_offset += 1;
-                        current_node_index = node.primitives_offset;
-                    } else {
-                        nodes_to_visit[to_visit_offset] = node.primitives_offset;
-                        to_visit_offset += 1;
-                        current_node_index += 1;
-                    }
-                }
-            } else {
-                if to_visit_offset == 0 {
-                    break;
-                }
-                to_visit_offset -= 1;
-                current_node_index = nodes_to_visit[to_visit_offset];
-            }
-        }
-        false
     }
 }
 
-impl BVHAggregate<Primitive> {
+fn leaf_order(node: &BVHBuildNode, out: &mut Vec<usize>) {
-    pub fn to_device(&self, arena: &mut Arena) -> DeviceBVHAggregate {
+    match node {
-        let shared_nodes: Vec<shared::core::aggregates::LinearBVHNode> = self
+        BVHBuildNode::Leaf {
-            .nodes
+            primitive_indices, ..
-            .iter()
+        } => {
-            .map(|n| shared::core::aggregates::LinearBVHNode {
+            out.extend_from_slice(primitive_indices);
-                bounds: n.bounds,
+        }
-                primitives_offset: n.primitives_offset,
+        BVHBuildNode::Interior { children, .. } => {
-                n_primitives: n.n_primitives,
+            leaf_order(&children[0], out);
-                axis: n.axis,
+            leaf_order(&children[1], out);
-                pad: 0,
-            })
-            .collect();
-
-        DeviceBVHAggregate {
-            max_prims_in_node: self.max_prims_in_node as u32,
-            primitives: gvec_from_slice(&self.primitives),
-            primitive_count: self.primitives.len() as u32,
-            nodes: gvec_from_slice(&shared_nodes),
-            node_count: self.nodes.len() as u32,
         }
     }
 }
+
+// BVHAggregate *BVHAggregate::Create(std::vector<Primitive> prims,
+//                                    const ParameterDictionary &parameters) {
+//     std::string splitMethodName = parameters.GetOneString("splitmethod", "sah");
+//     BVHAggregate::SplitMethod splitMethod;
+//     if (splitMethodName == "sah")
+//         splitMethod = BVHAggregate::SplitMethod::SAH;
+//     else if (splitMethodName == "hlbvh")
+//         splitMethod = BVHAggregate::SplitMethod::HLBVH;
+//     else if (splitMethodName == "middle")
+//         splitMethod = BVHAggregate::SplitMethod::Middle;
+//     else if (splitMethodName == "equal")
+//         splitMethod = BVHAggregate::SplitMethod::EqualCounts;
+//     else {
+//         Warning(R"(BVH split method "%s" unknown.  Using "sah".)", splitMethodName);
+//         splitMethod = BVHAggregate::SplitMethod::SAH;
+//     }
+//
+//     int maxPrimsInNode = parameters.GetOneInt("maxnodeprims", 4);
+//     return new BVHAggregate(std::move(prims), maxPrimsInNode, splitMethod);
+// }

src/core/bssrdf.rs

@@ -1,40 +0,0 @@
-use shared::Float;
-use shared::Ptr;
-use shared::core::bssrdf::BSSRDFTable;
-
-pub struct BSSRDFTableData {
-    pub rho_samples: Vec<Float>,
-    pub radius_samples: Vec<Float>,
-    pub profile: Vec<Float>,
-    pub rho_eff: Vec<Float>,
-    pub profile_cdf: Vec<Float>,
-}
-
-impl BSSRDFTableData {
-    pub fn new(n_rho_samples: usize, n_radius_samples: usize) -> Self {
-        let rho_samples: Vec<Float> = Vec::with_capacity(n_rho_samples);
-        let radius_samples: Vec<Float> = Vec::with_capacity(n_radius_samples);
-        let profile: Vec<Float> = Vec::with_capacity(n_radius_samples * n_rho_samples);
-        let rho_eff: Vec<Float> = Vec::with_capacity(n_rho_samples);
-        let profile_cdf: Vec<Float> = Vec::with_capacity(n_radius_samples * n_rho_samples);
-        Self {
-            rho_samples,
-            radius_samples,
-            profile,
-            rho_eff,
-            profile_cdf,
-        }
-    }
-
-    pub fn view(&self, rho_ptr: *const Float, radius_ptr: *const Float) -> BSSRDFTable {
-        BSSRDFTable {
-            rho_samples: rho_ptr.into(),
-            n_rho: self.rho_samples.len() as u32,
-            radius_samples: radius_ptr.into(),
-            n_radius: self.radius_samples.len() as u32,
-            profile: Ptr::from(self.profile.as_ptr()),
-            profile_cdf: Ptr::from(self.profile_cdf.as_ptr()),
-            rho_eff: Ptr::from(self.rho_eff.as_ptr()),
-        }
-    }
-}

src/core/film.rs

@@ -236,7 +236,7 @@ pub trait CreateFilmBase {
     fn create(
         params: &ParameterDictionary,
         filter: Filter,
-        sensor: Option<&PixelSensor>,
+        sensor: Ptr<PixelSensor>,
         loc: &FileLoc,
     ) -> Result<Self>
     where
@@ -247,7 +247,7 @@ impl CreateFilmBase for FilmBase {
     fn create(
         params: &ParameterDictionary,
         filter: Filter,
-        sensor: Option<&PixelSensor>,
+        sensor: Ptr<PixelSensor>,
         loc: &FileLoc,
     ) -> Result<Self>
     where
@@ -301,7 +301,7 @@ impl CreateFilmBase for FilmBase {
             pixel_bounds,
             filter,
             diagonal: diagonal_mm * 0.001,
-            sensor: Ptr::from(sensor.unwrap()),
+            sensor,
         })
     }
 }

src/core/mod.rs

@@ -1,5 +1,5 @@
 pub mod aggregates;
-pub mod bssrdf;
+// pub mod bssrdf;
 pub mod camera;
 pub mod color;
 pub mod film;
@@ -15,3 +15,4 @@ pub mod scene;
 pub mod shape;
 pub mod spectrum;
 pub mod texture;
+pub mod render;

src/core/render.rs

@@ -0,0 +1,20 @@
+use crate::core::scene::BasicScene;
+use crate::Arena;
+use anyhow::Result;
+use log::warn;
+use shared::core::camera::CameraTrait;
+
+fn render_scene(scene: &BasicScene, arena: &Arena) -> Result<()> {
+    let media = scene.create_media();
+    let textures = scene.create_textures(arena);
+    let lights = scene.create_lights()
+    let (named_materials, materials) = scene.create_materials(&textures, arena)?;
+    let (aggregate, area_lights) =
+        scene.create_aggregate(&textures, &named_materials, &materials, arena);
+    let camera = scene.get_camera().unwrap();
+    let film = camera.get_film();
+    warn!("Creating integrator");
+    let sampler = scene.get_sampler()?;
+    let integrator = scene.create_integrator(camera, sampler, aggregate, lights, arena);
+    Ok(())
+}

src/core/scene/builder.rs

@@ -1,18 +1,19 @@
-use super::BasicScene;
 use super::entities::*;
-use crate::Arena;
+use super::BasicScene;
 use crate::spectra::get_colorspace_device;
 use crate::utils::error::FileLoc;
 use crate::utils::normalize_utf8;
 use crate::utils::parameters::{ParameterDictionary, ParsedParameterVector};
 use crate::utils::parser::{ParserError, ParserTarget};
-use shared::Float;
+use crate::Arena;
+use anyhow::Context;
 use shared::core::camera::CameraTransform;
 use shared::core::geometry::Vector3f;
 use shared::spectra::RGBColorSpace;
 use shared::utils::options::RenderingCoordinateSystem;
 use shared::utils::transform;
 use shared::utils::transform::{AnimatedTransform, Transform};
+use shared::Float;
 use std::collections::{HashMap, HashSet};
 use std::ops::{Index, IndexMut};
 use std::sync::Arc;
@@ -279,15 +280,9 @@ impl ParserTarget for BasicSceneBuilder {
         uz: Float,
         loc: FileLoc,
     ) -> Result<(), ParserError> {
-        let result = transform::look_at((ex, ey, ez), (lx, ly, lz), (ux, uy, uz));
+        let t = transform::look_at((ex, ey, ez), (lx, ly, lz), (ux, uy, uz))
-        match result {
+            .with_context(|| format!("at {}", loc))?;
-            Some(t) => {
+        self.for_active_transforms(|cur| cur * &t);
-                self.for_active_transforms(|cur| cur * &t);
-            }
-            None => {
-                eprintln!("Error: Could not invert transform at {}", loc);
-            }
-        }
         Ok(())
     }
 
@@ -640,7 +635,10 @@ impl ParserTarget for BasicSceneBuilder {
         loc: FileLoc,
         arena: Arc<Arena>,
     ) -> Result<(), ParserError> {
-        eprintln!("TEXTURE: name='{}' type='{}' tex='{}'", orig_name, type_name, tex_name);
+        eprintln!(
+            "TEXTURE: name='{}' type='{}' tex='{}'",
+            orig_name, type_name, tex_name
+        );
         let name = normalize_utf8(orig_name);
         self.verify_world("Texture", &loc)?;
         let dict = ParameterDictionary::from_array(

src/core/scene/scene.rs

@@ -1,5 +1,6 @@
 use super::entities::*;
 use super::state::*;
+use crate::core::aggregates::CreateBVH;
 use crate::core::camera::CameraFactory;
 use crate::core::film::FilmFactory;
 use crate::core::filter::FilterFactory;
@@ -9,12 +10,14 @@ use crate::core::primitive::{CreateGeometricPrimitive, CreateSimplePrimitive};
 use crate::core::sampler::SamplerFactory;
 use crate::core::shape::{ShapeFactory, ShapeWithContext};
 use crate::core::texture::{FloatTexture, SpectrumTexture};
+use crate::integrators::{CreateIntegrator, PathConfig, PathIntegrator};
 use crate::utils::parallel::{run_async, AsyncJob};
 use crate::utils::parameters::{NamedTextures, ParameterDictionary, TextureParameterDictionary};
 use crate::utils::resolve_filename;
 use crate::{Arena, ArenaUpload, FileLoc};
 use anyhow::{anyhow, Result};
 use parking_lot::Mutex;
+use shared::core::aggregates::{BVHAggregate, SplitMethod};
 use shared::core::camera::{Camera, CameraTransform};
 use shared::core::color::LINEAR;
 use shared::core::film::Film;
@@ -26,6 +29,7 @@ use shared::core::primitive::{AnimatedPrimitive, GeometricPrimitive, Primitive,
 use shared::core::sampler::Sampler;
 use shared::core::shape::Shape;
 use shared::core::texture::{GPUFloatTexture, SpectrumType};
+use shared::lights::sampler::LightSampler;
 use shared::spectra::RGBColorSpace;
 use shared::{Ptr, Transform};
 use std::collections::HashMap;
@@ -347,7 +351,7 @@ impl BasicScene {
         self.instances.lock().extend(uses);
     }
 
-    pub fn create_textures(&self, arena: &mut Arena) -> NamedTextures {
+    pub fn create_textures(&self, arena: &Arena) -> NamedTextures {
         let mut state = self.texture_state.lock();
 
         let mut float_textures: HashMap<String, Arc<FloatTexture>> = HashMap::new();
@@ -408,7 +412,7 @@ impl BasicScene {
     pub fn create_materials(
         &self,
         textures: &NamedTextures,
-        arena: &mut Arena,
+        arena: &Arena,
     ) -> Result<(HashMap<String, Material>, Vec<Material>)> {
         let mut state = self.material_state.lock();
 
@@ -619,8 +623,8 @@ impl BasicScene {
         textures: &NamedTextures,
         named_materials: &HashMap<String, Material>,
         materials: &[Material],
-        arena: &mut Arena,
+        arena: &Arena,
-    ) -> (Vec<Primitive>, Vec<Arc<Light>>) {
+    ) -> (Arc<Primitive>, Vec<Arc<Light>>) {
         let entities = self.shapes.lock();
         let animated = self.animated_shapes.lock();
         let light_state = self.light_state.lock();
@@ -660,7 +664,36 @@ impl BasicScene {
             );
         }
 
-        (primitives, area_lights)
+        let aggregate = if !primitives.is_empty() {
+            BVHAggregate::new(primitives.clone(), 4, SplitMethod::SAH)
+        } else {
+            BVHAggregate::empty()
+        };
+
+        let agg_ptr = arena.alloc(aggregate);
+
+        (Arc::new(Primitive::BVH(agg_ptr)), area_lights)
+    }
+
+    pub fn create_integrator(
+        &self,
+        camera: Arc<Camera>,
+        sampler: Arc<Sampler>,
+        aggregate: Arc<Primitive>,
+        lights: Vec<Arc<Light>>,
+        arena: &Arena
+    ) -> PathIntegrator {
+        let integrator = &self.integrator.lock().clone().unwrap();
+        PathIntegrator::create(
+            integrator.parameters.clone(),
+            camera,
+            sampler,
+            aggregate,
+            lights,
+            PathConfig::SIMPLE,
+            arena
+        )
+        .expect("Integrator creation has failed")
     }
 
     fn build_primitives_inner(
@@ -724,7 +757,7 @@ impl BasicScene {
         light_state: &LightState,
         media: &MediaState,
         film_cs: Option<&RGBColorSpace>,
-        arena: &mut Arena,
+        arena: &Arena,
         primitives: &mut Vec<Primitive>,
         area_lights: &mut Vec<Arc<Light>>,
     ) {
@@ -807,7 +840,7 @@ impl BasicScene {
         light_state: &LightState,
         media: &MediaState,
         film_cs: Option<&RGBColorSpace>,
-        arena: &mut Arena,
+        arena: &Arena,
         primitives: &mut Vec<Primitive>,
         area_lights: &mut Vec<Arc<Light>>,
     ) {
@@ -973,6 +1006,20 @@ impl BasicScene {
         Vec::new()
     }
 
+    pub fn create_media(&self) -> HashMap<String, Arc<Medium>> {
+        let mut state = self.media_state.lock();
+
+        if !state.jobs.is_empty() {
+            let jobs: Vec<(String, AsyncJob<Medium>)> = state.jobs.drain().collect();
+            for (name, job) in jobs {
+                let medium = Arc::new(job.wait());
+                state.map.insert(name, medium);
+            }
+        }
+
+        state.map.clone()
+    }
+
     // Getters
 
     pub fn get_camera(&self) -> Result<Arc<Camera>> {

src/films/gbuffer.rs

@@ -21,7 +21,8 @@ impl CreateFilm for GBufferFilm {
         let max_component_value = params.get_one_float("maxcomponentvalue", Float::INFINITY)?;
         let write_fp16 = params.get_one_bool("savefp16", true)?;
         let sensor = PixelSensor::create(params, colorspace.clone(), exposure_time, loc, arena)?;
-        let film_base = FilmBase::create(params, filter, Some(&sensor), loc)?;
+        let sensor_ptr = arena.alloc(sensor);
+        let film_base = FilmBase::create(params, filter, sensor_ptr, loc)?;
 
         let filename = params.get_one_string("filename", "pbrt.exr")?;
         if Path::new(&filename).extension() != Some("exr".as_ref()) {

src/films/rgb.rs

@@ -19,7 +19,8 @@ impl CreateFilm for RGBFilm {
         let max_component_value = params.get_one_float("maxcomponentvalue", Float::INFINITY)?;
         let write_fp16 = params.get_one_bool("savefp16", true)?;
         let sensor = PixelSensor::create(params, colorspace.clone(), exposure_time, loc, arena)?;
-        let film_base = FilmBase::create(params, filter, Some(&sensor), loc)?;
+        let sensor_ptr = arena.alloc(sensor);
+        let film_base = FilmBase::create(params, filter, sensor_ptr, loc)?;
         let film = RGBFilm::new(film_base, &colorspace, max_component_value, write_fp16);
         Ok(Film::RGB(film))
     }

src/films/spectral.rs

@@ -23,7 +23,8 @@ impl CreateFilm for SpectralFilm {
         let max_component_value = params.get_one_float("maxcomponentvalue", Float::INFINITY)?;
         let write_fp16 = params.get_one_bool("savefp16", true)?;
         let sensor = PixelSensor::create(params, colorspace.clone(), exposure_time, loc, arena)?;
-        let film_base = FilmBase::create(params, filter, Some(&sensor), loc)?;
+        let sensor_ptr = arena.alloc(sensor);
+        let film_base = FilmBase::create(params, filter, sensor_ptr, loc)?;
 
         let filename = params.get_one_string("filename", "pbrt.exr")?;
         if Path::new(&filename).extension() != Some("exr".as_ref()) {

src/integrators/base.rs

@@ -18,7 +18,24 @@ pub struct IntegratorBase {
 }
 
 impl IntegratorBase {
-    pub fn new(aggregate: Arc<Primitive>, lights: Vec<Arc<Light>>) -> Self {
+    pub fn new(aggregate: Arc<Primitive>, mut lights: Vec<Arc<Light>>) -> Self {
+        let scene_bounds = aggregate.bounds();
+
+        for light in &mut lights {
+            Arc::get_mut(light)
+                .expect("Light has multiple owners during setup")
+                .preprocess(&scene_bounds);
+        }
+
+        println!(
+            "IntegratorBase: {} lights, scene_bounds: {:?}",
+            lights.len(),
+            scene_bounds
+        );
+        for (i, l) in lights.iter().enumerate() {
+            println!("  light[{}]: type={:?}", i, l.light_type());
+        }
+
         let infinite_lights = lights
             .iter()
             .filter(|light| light.light_type().is_infinite())

src/integrators/mod.rs

@@ -4,13 +4,20 @@ pub mod path;
 pub mod pipeline;
 pub mod state;
 
-pub use path::PathIntegrator;
+pub use path::{PathConfig, PathIntegrator};
 
-use crate::Arena;
+use crate::lights::sampler::create_light_sampler;
+use crate::{Arena, ParameterDictionary};
+use anyhow::Result;
+use shared::core::camera::Camera;
 use shared::core::film::VisibleSurface;
 use shared::core::geometry::{Point2i, Ray};
+use shared::core::light::Light;
+use shared::core::primitive::Primitive;
 use shared::core::sampler::Sampler;
+use shared::lights::sampler::LightSampler;
 use shared::spectra::{SampledSpectrum, SampledWavelengths};
+use std::sync::Arc;
 
 pub trait IntegratorTrait {
     fn render(&self);
@@ -34,3 +41,33 @@ pub trait RayIntegratorTrait {
         arena: &Arena,
     ) -> (SampledSpectrum, Option<VisibleSurface>);
 }
+
+pub trait CreateIntegrator {
+    fn create(
+        parameters: ParameterDictionary,
+        camera: Arc<Camera>,
+        sampler: Arc<Sampler>,
+        aggregate: Arc<Primitive>,
+        lights: Vec<Arc<Light>>,
+        config: PathConfig,
+        arena: &Arena,
+    ) -> Result<PathIntegrator>;
+}
+
+impl CreateIntegrator for PathIntegrator {
+    fn create(
+        parameters: ParameterDictionary,
+        camera: Arc<Camera>,
+        sampler: Arc<Sampler>,
+        aggregate: Arc<Primitive>,
+        lights: Vec<Arc<Light>>,
+        config: PathConfig,
+        arena: &Arena,
+    ) -> Result<PathIntegrator> {
+        let _max_depth = parameters.get_one_int("maxdepth", 5)?;
+        let _regularize = parameters.get_one_bool("regularize", false)?;
+        let light_sampler = create_light_sampler("bvh", &lights, arena);
+        let integrator = PathIntegrator::new(aggregate, lights, camera, light_sampler, config);
+        Ok(integrator)
+    }
+}

src/integrators/path.rs

@@ -81,7 +81,7 @@ impl PathIntegrator {
         sampler: LightSampler,
         config: PathConfig,
     ) -> Self {
-        let base = IntegratorBase::new(aggregate, lights.clone());
+        let base = IntegratorBase::new(aggregate, lights);
         Self {
             base,
             camera,

src/integrators/pipeline.rs

@@ -82,6 +82,7 @@ pub fn render<T>(
 ) where
     T: RayIntegratorTrait + Sync,
 {
+    println!("RENDER CALLED");
     let options = get_options();
     if let Some((p_pixel, sample_index)) = options.debug_start {
         let s_index = sample_index as usize;
@@ -101,6 +102,11 @@ pub fn render<T>(
     }
 
     let pixel_bounds = camera.get_film().pixel_bounds();
+    println!(
+        "pixel_bounds: {:?}, area: {}",
+        pixel_bounds,
+        pixel_bounds.area()
+    );
     let spp = sampler_prototype.samples_per_pixel();
     let total_work = (pixel_bounds.area() as u64) * (spp as u64);
     let progress = PbrtProgress::new(total_work, "Rendering", options.quiet);
@@ -162,6 +168,7 @@ pub fn render<T>(
             for p_pixel in tile_bounds {
                 for sample_index in wave_start..wave_end {
                     sampler.start_pixel_sample(*p_pixel, sample_index, None);
+                    println!("Evaluating pixel {:?} sample {}", p_pixel, sample_index);
                     evaluate_pixel_sample(
                         integrator,
                         camera,
@@ -232,7 +239,7 @@ pub fn evaluate_pixel_sample<T: RayIntegratorTrait>(
     camera: &Camera,
     sampler: &mut Sampler,
     pixel: Point2i,
-    _sample_index: usize,
+    sample_index: usize,
     arena: &Arena,
 ) {
     let mut lu = sampler.get1d();
@@ -267,6 +274,13 @@ pub fn evaluate_pixel_sample<T: RayIntegratorTrait>(
             l = SampledSpectrum::new(0.);
         }
 
+        if pixel.x() == 352 && pixel.y() == 352 && sample_index == 0 {
+            println!("Center pixel: L = {:?}", l);
+            println!("  ray origin: {:?}", camera_ray.ray.o);
+            println!("  ray dir: {:?}", camera_ray.ray.d);
+            println!("  camera_sample.p_film: {:?}", camera_sample.p_film);
+        }
+
         film.add_sample(
             pixel,
             l,

src/lights/distant.rs

@@ -15,11 +15,11 @@ use shared::utils::{Ptr, Transform};
 use shared::Float;
 
 trait CreateDistantLight {
-    fn new(render_from_light: Transform, le: Spectrum, scale: Float) -> Self;
+    fn new(render_from_light: Transform, le: Spectrum, scale: Float, arena: &Arena) -> Self;
 }
 
 impl CreateDistantLight for DistantLight {
-    fn new(render_from_light: Transform, le: Spectrum, scale: Float) -> Self {
+    fn new(render_from_light: Transform, le: Spectrum, scale: Float, arena: &Arena) -> Self {
         let base = LightBase::new(
             LightType::DeltaDirection,
             render_from_light,
@@ -28,7 +28,7 @@ impl CreateDistantLight for DistantLight {
         let lemit = lookup_spectrum(&le);
         Self {
             base,
-            lemit: Ptr::from(&*lemit),
+            lemit: arena.alloc_arc(lemit),
             scale,
             scene_center: Point3f::default(),
             scene_radius: 0.,
@@ -44,7 +44,7 @@ pub fn create(
     _shape: &Shape,
     _alpha_text: &FloatTexture,
     colorspace: Option<&RGBColorSpace>,
-    _arena: &Arena,
+    arena: &Arena,
 ) -> Result<Light> {
     let default_cs = crate::spectra::default_colorspace();
     let cs = colorspace.unwrap_or(&default_cs);
@@ -88,7 +88,7 @@ pub fn create(
         scale *= e_v;
     }
 
-    let specific = DistantLight::new(final_render, l, scale);
+    let specific = DistantLight::new(final_render, l, scale, arena);
 
     Ok(Light::Distant(specific))
 }

src/lights/infinite.rs

@@ -3,7 +3,7 @@ use crate::core::light::lookup_spectrum;
 use crate::core::spectrum::spectrum_to_photometric;
 use crate::spectra::get_spectra_context;
 use crate::utils::resolve_filename;
-use crate::{Arena, FileLoc, ParameterDictionary, ArenaUpload, Upload};
+use crate::{Arena, FileLoc, ParameterDictionary, ArenaUpload};
 use anyhow::{anyhow, Result};
 use rayon::prelude::*;
 use shared::core::camera::CameraTransform;

src/lights/sampler.rs

@@ -9,7 +9,6 @@ use shared::utils::Ptr;
 use shared::Float;
 use std::sync::Arc;
 
-/// Top-level dispatcher matching the C++ LightSampler::Create.
 pub fn create_light_sampler(
     name: &str,
     lights: &[Arc<Light>],