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/cameras/orthographic.rs

@@ -96,7 +96,7 @@ impl CameraTrait for OrthographicCamera {
             p_camera,
             Vector3f::new(0., 0., 1.),
             Some(self.sample_time(sample.time)),
-            &self.base().medium,
+            self.base().medium,
         );
         if self.lens_radius > 0. {
             let p_lens_vec =

shared/src/cameras/perspective.rs

@@ -98,7 +98,7 @@ impl CameraTrait for PerspectiveCamera {
             Point3f::new(0., 0., 0.),
             p_vector.normalize(),
             Some(self.sample_time(sample.time)),
-            &*self.base().medium,
+            self.base().medium,
         );
         // Modify ray for depth of field
         if self.lens_radius > 0. {

shared/src/cameras/realistic.rs

@@ -132,7 +132,7 @@ impl RealisticCamera {
             Point3f::new(0., x, self.lens_front_z() + 1.),
             Vector3f::new(0., 0., -1.),
             None,
-            &Ptr::null(),
+            Ptr::null(),
         );
         let Some((_, r_film)) = self.trace_lenses_from_film(&r_scene) else {
             panic!(
@@ -144,7 +144,7 @@ impl RealisticCamera {
             Point3f::new(x, 0., self.lens_rear_z() - 1.),
             Vector3f::new(0., 0., 1.),
             None,
-            &Ptr::null(),
+            Ptr::null(),
         );
         let Some((_, r_scene)) = self.trace_lenses_from_film(&r_film) else {
             panic!(
@@ -209,7 +209,7 @@ impl RealisticCamera {
             // Expand pupil bounds if ray makes it through the lens system
             if !pupil_bounds.contains(Point2f::new(p_rear.x(), p_rear.y()))
                 && trace_lenses_from_film(
-                    Ray::new(p_film, p_rear - p_film, None, &Ptr::null()),
+                    Ray::new(p_film, p_rear - p_film, None, Ptr::null()),
                     None,
                 )
             {
@@ -270,7 +270,7 @@ impl RealisticCamera {
             Point3f::new(r_camera.o.x(), r_camera.o.y(), -r_camera.o.z()),
             Vector3f::new(r_camera.d.x(), r_camera.d.y(), -r_camera.d.z()),
             Some(r_camera.time),
-            &Ptr::null(),
+            Ptr::null(),
         );
 
         for i in (0..self.n_elements - 1).rev() {
@@ -337,7 +337,7 @@ impl RealisticCamera {
             Point3f::new(r_lens.o.x(), r_lens.o.y(), -r_lens.o.z()),
             Vector3f::new(r_lens.d.x(), r_lens.d.y(), -r_lens.d.z()),
             Some(r_lens.time),
-            &Ptr::null(),
+            Ptr::null(),
         );
 
         Some((weight, r_out))
@@ -415,7 +415,7 @@ impl CameraTrait for RealisticCamera {
         let eps = self.sample_exit_pupil(Point2f::new(p_film.x(), p_film.y()), sample.p_lens)?;
 
         let p_pupil = Point3f::new(0., 0., 0.);
-        let r_film = Ray::new(p_film, p_pupil - p_film, None, &Ptr::null());
+        let r_film = Ray::new(p_film, p_pupil - p_film, None, Ptr::null());
         let (weight, mut ray) = self.trace_lenses_from_film(&r_film)?;
         if weight == 0. {
             return None;

shared/src/cameras/spherical.rs

@@ -54,7 +54,7 @@ impl CameraTrait for SphericalCamera {
             Point3f::new(0., 0., 0.),
             dir,
             Some(self.sample_time(sample.time)),
-            &self.base().medium,
+            self.base().medium,
         );
         Some(CameraRay {
             ray: self.render_from_camera(&ray, &mut None),

shared/src/core/aggregates.rs

@@ -1,10 +1,19 @@
-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)]
+#[derive(Debug, Clone, Copy, PartialEq, Eq)]
+pub enum SplitMethod {
+    SAH,
+    Hlbvh,
+    Middle,
+    EqualCounts,
+}
+
+#[repr(C)]
+#[derive(Default, Debug, Clone, Copy)]
 pub struct LinearBVHNode {
     pub bounds: Bounds3f,
     pub primitives_offset: usize,
@@ -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 stack_ptr = 0;
-        let mut node_idx = 0usize;
+        let mut to_visit_offset = 0;
+        let mut current_node_index = 0;
+        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 {
-                    break;
-                }
-                stack_ptr -= 1;
-                node_idx = stack[stack_ptr];
-                continue;
-            }
+                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 node.n_primitives > 0 {
-                // Leaf: test all primitives
-                for i in 0..node.n_primitives {
-                    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();
-                        best_si = Some(si);
+                        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;
                     }
                 }
-
-                if stack_ptr == 0 {
+            } else {
+                // The ray missed the AABB of this node. Pop stack to try the next node.
+                if to_visit_offset == 0 {
                     break;
                 }
-                stack_ptr -= 1;
-                node_idx = stack[stack_ptr];
-            } 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
-                }
+                to_visit_offset -= 1;
+                current_node_index = nodes_to_visit[to_visit_offset];
             }
         }
 
@@ -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 stack_ptr = 0;
-        let mut node_idx = 0usize;
+        let mut to_visit_offset = 0;
+        let mut current_node_index = 0;
+        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 {
-                    break;
-                }
-                stack_ptr -= 1;
-                node_idx = stack[stack_ptr];
-                continue;
-            }
+                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 node.n_primitives > 0 {
-                for i in 0..node.n_primitives {
-                    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;
+                        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;
                     }
                 }
-
-                if stack_ptr == 0 {
+            } else {
+                if to_visit_offset == 0 {
                     break;
                 }
-                stack_ptr -= 1;
-                node_idx = stack[stack_ptr];
-            } 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;
-                }
+                to_visit_offset -= 1;
+                current_node_index = nodes_to_visit[to_visit_offset];
             }
         }
-
         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::utils::Ptr;
+use crate::spectra::{SampledSpectrum, N_SPECTRUM_SAMPLES};
 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 radius_samples: Ptr<Float>,
-    pub profile: Ptr<Float>,
-    pub rho_eff: Ptr<Float>,
-    pub profile_cdf: Ptr<Float>,
+    pub rho_samples: GVec<Float>,
+    pub radius_samples: GVec<Float>,
+    pub profile: GVec<Float>,
+    pub rho_eff: GVec<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;
-        unsafe { core::slice::from_raw_parts(self.profile.as_ref(), n_profile) }
+        // let n_profile = (self.n_rho * self.n_radius) as usize;
+        &self.profile
     }
 
     pub fn get_cdf(&self) -> &[Float] {
-        let n_profile = (self.n_rho * self.n_radius) as usize;
-        unsafe { core::slice::from_raw_parts(self.profile_cdf.as_ref(), n_profile) }
+        // let n_profile = (self.n_rho * self.n_radius) as usize;
+        &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/core/camera.rs

@@ -227,13 +227,13 @@ pub trait CameraTrait {
             Point3f::new(0., 0., 0.) + self.base().min_pos_differential_x,
             Vector3f::new(0., 0., 1.) + self.base().min_dir_differential_x,
             None,
-            &Ptr::default(),
+            Ptr::default(),
         );
         let y_ray = Ray::new(
             Point3f::new(0., 0., 0.) + self.base().min_pos_differential_y,
             Vector3f::new(0., 0., 1.) + self.base().min_dir_differential_y,
             None,
-            &Ptr::default(),
+            Ptr::default(),
         );
         let n_down = Vector3f::from(n_down_z);
         let tx = -(n_down.dot(y_ray.o.into())) / n_down.dot(x_ray.d);

shared/src/core/geometry/ray.rs

@@ -30,12 +30,12 @@ impl Default for Ray {
 }
 
 impl Ray {
-    pub fn new(o: Point3f, d: Vector3f, time: Option<Float>, medium: &Medium) -> Self {
+    pub fn new(o: Point3f, d: Vector3f, time: Option<Float>, medium: Ptr<Medium>) -> Self {
         Self {
             o,
             d,
             time: time.unwrap_or_else(|| Self::default().time),
-            medium: Ptr::from(medium),
+            medium,
             ..Self::default()
         }
     }

shared/src/core/image.rs

@@ -167,6 +167,17 @@ impl Pixels {
         }
     }
 
+    pub fn as_u8(&self) -> &[u8] {
+        &self.data
+    }
+
+    pub fn as_f16(&mut self) -> &[u16] {
+        assert_eq!(self.format, PixelFormat::F16);
+        unsafe {
+            core::slice::from_raw_parts(self.data.as_ptr() as *const u16, self.data.len() / 2)
+        }
+    }
+
     pub fn as_u8_mut(&mut self) -> &mut [u8] {
         &mut self.data
     }

shared/src/core/primitive.rs

@@ -1,5 +1,5 @@
-use crate::core::geometry::{Bounds3f, Ray};
 use crate::core::aggregates::BVHAggregate;
+use crate::core::geometry::{Bounds3f, Ray};
 use crate::core::interaction::{Interaction, InteractionTrait, SurfaceInteraction};
 use crate::core::light::Light;
 use crate::core::material::Material;
@@ -104,7 +104,12 @@ impl PrimitiveTrait for SimplePrimitive {
 
     fn intersect(&self, r: &Ray, t_max: Option<Float>) -> Option<ShapeIntersection> {
         let mut si = self.shape.intersect(r, t_max)?;
-        si.set_intersection_properties(self.material, Ptr::null(), MediumInterface::default(), r.medium);
+        si.set_intersection_properties(
+            self.material,
+            Ptr::null(),
+            MediumInterface::default(),
+            r.medium,
+        );
         Some(si)
     }
 
@@ -183,7 +188,7 @@ impl PrimitiveTrait for AnimatedPrimitive {
 }
 
 #[repr(C)]
-#[derive(Debug, Clone, Copy)]
+#[derive(Default, Debug, Clone, Copy)]
 pub struct LinearBVHNode {
     bounds: Bounds3f,
 }
@@ -217,7 +222,6 @@ pub enum Primitive {
     KdTree(KdTreeAggregate),
 }
 
-
 impl<T: PrimitiveTrait> PrimitiveTrait for Ptr<T> {
     fn bounds(&self) -> Bounds3f {
         unsafe { self.as_ref().bounds() }
@@ -231,5 +235,3 @@ impl<T: PrimitiveTrait> PrimitiveTrait for Ptr<T> {
         unsafe { self.as_ref().intersect_p(r, t_max) }
     }
 }
-
-

shared/src/core/shape.rs

@@ -118,7 +118,7 @@ impl ShapeSampleContext {
     }
 
     pub fn spawn_ray(&self, w: Vector3f) -> Ray {
-        Ray::new(self.offset_ray_origin(w), w, Some(self.time), &Ptr::null())
+        Ray::new(self.offset_ray_origin(w), w, Some(self.time), Ptr::null())
     }
 }
 

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/shapes/mod.rs

@@ -12,4 +12,4 @@ pub use cylinder::*;
 pub use disk::*;
 pub use sphere::*;
 pub use triangle::*;
-pub use mesh::*;
+pub use mesh::{TriangleMesh, BilinearPatchMesh};

shared/src/shapes/sphere.rs

@@ -57,7 +57,7 @@ impl SphereShape {
         phi_max: Float,
     ) -> Self {
         let theta_z_min = clamp(z_min.min(z_max) / radius, -1., 1.).acos();
-        let theta_z_max = clamp(z_max.min(z_max) / radius, -1., 1.).acos();
+        let theta_z_max = clamp(z_min.max(z_max) / radius, -1., 1.).acos();
         let phi_max = radians(clamp(phi_max, 0., 360.0));
         Self {
             render_from_object: render_from_object.clone(),

shared/src/shapes/triangle.rs

@@ -78,7 +78,6 @@ impl TriangleShape {
         Some([mesh.s[v0], mesh.s[v1], mesh.s[v2]])
     }
 
-
     fn get_uvs(&self) -> Option<[Point2f; 3]> {
         let mesh = self.mesh();
         if mesh.s.is_empty() {
@@ -107,13 +106,125 @@ impl TriangleShape {
 
     fn intersect_triangle(
         &self,
-        _ray: &Ray,
-        _t_max: Float,
-        _p0: Point3f,
-        _p1: Point3f,
-        _p2: Point3f,
+        ray: &Ray,
+        t_max: Float,
+        p0: Point3f,
+        p1: Point3f,
+        p2: Point3f,
     ) -> Option<TriangleIntersection> {
-        todo!()
+        if (p2 - p0).cross(p1 - p0).norm_squared() == 0. {
+            return None;
+        }
+
+        // Transform triangle vertices to ray coordinate space
+        // Translate vertices based on ray origin
+        let mut p0t = p0 - Vector3f::from(ray.o);
+        let mut p1t = p1 - Vector3f::from(ray.o);
+        let mut p2t = p2 - Vector3f::from(ray.o);
+
+        // Permute components of triangle vertices and ray direction
+        let kz = ray.d.abs().max_component_index();
+        let mut kx = kz + 1;
+        if kx == 3 {
+            kx = 0;
+        }
+        let mut ky = kx + 1;
+        if ky == 3 {
+            ky = 0;
+        }
+        let d = ray.d.permute([kx, ky, kz]);
+        p0t = p0t.permute([kx, ky, kz]);
+        p1t = p1t.permute([kx, ky, kz]);
+        p2t = p2t.permute([kx, ky, kz]);
+
+        // Apply shear transformation to translated vertex positions
+        let sx = -d.x() / d.z();
+        let sy = -d.y() / d.z();
+        let sz = 1. / d.z();
+        p0t[0] += sx * p0t[2];
+        p0t[1] += sy * p0t[2];
+        p1t[0] += sx * p1t[2];
+        p1t[1] += sy * p1t[2];
+        p2t[0] += sx * p2t[2];
+        p2t[1] += sy * p2t[2];
+
+        // Compute edge function coefficients _e0_, _e1_, and _e2_
+        let e0 = difference_of_products(p1t.x(), p2t.y(), p1t.y(), p2t.x());
+        let e1 = difference_of_products(p2t.x(), p0t.y(), p2t.y(), p0t.x());
+        let e2 = difference_of_products(p0t.x(), p1t.y(), p0t.y(), p1t.x());
+
+        // Fall back to double-precision test at triangle edges
+        // if sizeof(Float) == sizeof(float) && (e0 == 0.0f || e1 == 0.0f || e2 == 0.0f)) {
+        //     double p2txp1ty = (double)p2t.x * (double)p1t.y;
+        //     double p2typ1tx = (double)p2t.y * (double)p1t.x;
+        //     e0 = (float)(p2typ1tx - p2txp1ty);
+        //     double p0txp2ty = (double)p0t.x * (double)p2t.y;
+        //     double p0typ2tx = (double)p0t.y * (double)p2t.x;
+        //     e1 = (float)(p0typ2tx - p0txp2ty);
+        //     double p1txp0ty = (double)p1t.x * (double)p0t.y;
+        //     double p1typ0tx = (double)p1t.y * (double)p0t.x;
+        //     e2 = (float)(p1typ0tx - p1txp0ty);
+        // }
+
+        // Perform triangle edge and determinant tests
+        if (e0 < 0. || e1 < 0. || e2 < 0.) && (e0 > 0. || e1 > 0. || e2 > 0.) {
+            return None;
+        }
+        let det = e0 + e1 + e2;
+        if det == 0. {
+            return None;
+        }
+
+        // Compute scaled hit distance to triangle and test against ray $t$ range
+        p0t[2] *= sz;
+        p1t[2] *= sz;
+        p2t[2] *= sz;
+        let t_scaled = e0 * p0t.z() + e1 * p1t.z() + e2 * p2t.z();
+        if det < 0. && (t_scaled >= 0. || t_scaled < t_max * det) {
+            return None;
+        } else if det > 0. && (t_scaled <= 0. || t_scaled > t_max * det) {
+            return None;
+        }
+
+        // Compute barycentric coordinates and $t$ value for triangle intersection
+        let inv_det = 1. / det;
+        let b0 = e0 * inv_det;
+        let b1 = e1 * inv_det;
+        let b2 = e2 * inv_det;
+        let t = t_scaled * inv_det;
+
+        debug_assert!(t.is_finite());
+
+        // Ensure that computed triangle $t$ is conservatively greater than zero
+        // Compute $\delta_z$ term for triangle $t$ error bounds
+        let maxZt = Vector3f::new(p0t.z(), p1t.z(), p2t.z())
+            .abs()
+            .max_component_value();
+        let deltaZ = gamma(3) * maxZt;
+
+        // Compute $\delta_x$ and $\delta_y$ terms for triangle $t$ error bounds
+        let maxXt = Vector3f::new(p0t.x(), p1t.x(), p2t.x())
+            .abs()
+            .max_component_value();
+        let maxYt = Vector3f::new(p0t.y(), p1t.y(), p2t.y())
+            .abs()
+            .max_component_value();
+        let deltaX = gamma(5) * (maxXt + maxZt);
+        let deltaY = gamma(5) * (maxYt + maxZt);
+
+        // Compute $\delta_e$ term for triangle $t$ error bounds
+        let deltaE = 2. * (gamma(2) * maxXt * maxYt + deltaY * maxXt + deltaX * maxYt);
+
+        // Compute $\delta_t$ term for triangle $t$ error bounds and check _t_
+        let maxE = Vector3f::new(e0, e1, e2).abs().max_component_value();
+        let deltaT =
+            3. * (gamma(3) * maxE * maxZt + deltaE * maxZt + deltaZ * maxE) * inv_det.abs();
+        if t <= deltaT {
+            return None;
+        }
+
+        // Return _TriangleIntersection_ for intersection
+        return Some(TriangleIntersection { b0, b1, b2, t });
     }
 
     fn interaction_from_intersection(

shared/src/spectra/simple.rs

@@ -3,7 +3,8 @@ use super::sampled::{LAMBDA_MAX, LAMBDA_MIN};
 use crate::core::spectrum::{Spectrum, SpectrumTrait};
 use crate::spectra::{SampledSpectrum, SampledWavelengths, N_SPECTRUM_SAMPLES};
 use crate::utils::find_interval;
-use crate::{gvec, gvec_with_capacity, Float, GVec, Ptr};
+use crate::utils::math::square;
+use crate::{gvec, gvec_from_slice, gvec_with_capacity, Float, GVec, Ptr};
 use core::hash::{Hash, Hasher};
 use num_traits::Float as NumFloat;
 
@@ -88,6 +89,41 @@ impl DenselySampledSpectrum {
         }
     }
 
+    pub fn generate_cie_d(temperature: Float) -> Self {
+        let cct = temperature * 1.4388 / 1.4380;
+
+        if cct < 4000.0 {
+            return Self::from_function(
+                |lambda| BlackbodySpectrum::new(cct).evaluate(lambda),
+                LAMBDA_MIN,
+                LAMBDA_MAX,
+            );
+        }
+
+        let x = if cct < 7000. {
+            -4.607 * 1e9 / cct.powi(3) + 2.9678 * 1e6 / square(cct) + 0.09911 * 1e3 / cct + 0.244063
+        } else {
+            -2.0064 * 1e9 / cct.powi(3) + 1.9018 * 1e6 / square(cct) + 0.24748 * 1e3 / cct + 0.23704
+        };
+        let y = -3. * x + 2.87 * x - 0.275;
+        let m = 0.0241 + 0.2562 * x - 0.7341 * y;
+        let m1 = (-1.3515 - 1.7703 * x + 5.9114 * y) / m;
+        let m2 = (0.0300 - 31.4424 * x + 30.0717 * y) / m;
+
+        let mut coarse_values = gvec_with_capacity(N_CIES);
+        for i in 0..N_CIES {
+            coarse_values.push((CIE_S0[i] + CIE_S1[i] * m1 + CIE_S2[i] * m2) * 0.01);
+        }
+
+        let temp_pls = PiecewiseLinearSpectrum {
+            lambdas: gvec_from_slice(&CIE_S_LAMBDA),
+            values: gvec_from_slice(&coarse_values),
+            count: N_CIES as u32,
+        };
+
+        Self::from_function(|lambda| temp_pls.evaluate(lambda), LAMBDA_MIN, LAMBDA_MAX)
+    }
+
     pub fn scale(&mut self, s: Float) {
         for v in &mut self.values {
             *v *= s;
@@ -183,7 +219,6 @@ impl PiecewiseLinearSpectrum {
         self.count.try_into().unwrap()
     }
 
-
     #[inline(always)]
     pub fn lambda(&self, idx: u32) -> Float {
         unsafe { *self.lambdas.as_ptr().add(idx as usize) }

shared/src/utils/math.rs

@@ -113,7 +113,7 @@ pub fn safe_asin<T: NumFloat>(x: T) -> T {
 #[inline]
 pub fn safe_acos(x: Float) -> Float {
     if (-1.001..1.001).contains(&x) {
-        clamp(x, -1., 1.).asin()
+        clamp(x, -1., 1.).acos()
     } else {
         panic!("Not valid value for acos")
     }
@@ -542,9 +542,9 @@ pub fn next_float_up(v: Float) -> Float {
 
     let mut ui = float_to_bits(v);
     if v >= 0.0 {
-        ui += 1;
+        ui = ui.wrapping_add(1);
     } else {
-        ui -= 1;
+        ui = ui.wrapping_sub(1);
     }
     bits_to_float(ui)
 }
@@ -558,9 +558,9 @@ pub fn next_float_down(v: Float) -> Float {
 
     let mut ui = float_to_bits(v);
     if v > 0.0 {
-        ui -= 1;
+        ui = ui.wrapping_sub(1);
     } else {
-        ui += 1;
+        ui = ui.wrapping_add(1);
     }
     bits_to_float(ui)
 }
@@ -798,7 +798,11 @@ impl DigitPermutation {
             inv_base_m *= inv_base;
         }
 
-        let mut permutations = gvec_with_capacity(n_digits as usize * base as usize);
+        let mut permutations = {
+            let mut v = gvec_with_capacity(n_digits as usize * base as usize);
+            v.resize(n_digits as usize * base as usize, 0u16);
+            v
+        };
 
         for digit_index in 0..n_digits {
             let hash_input = [base as u64, digit_index as u64, seed];
@@ -829,7 +833,11 @@ impl DigitPermutation {
 
 pub fn compute_radical_inverse_permutations(seed: u64) -> GVec<DigitPermutation> {
     let mut result = gvec();
-    result.extend(PRIMES.iter().map(|&base| DigitPermutation::new(base as i32, seed)));
+    result.extend(
+        PRIMES
+            .iter()
+            .map(|&base| DigitPermutation::new(base as i32, seed)),
+    );
     result
 }
 

shared/src/utils/sampling.rs

@@ -774,6 +774,9 @@ impl PiecewiseConstant1D {
 
     pub fn find_interval(&self, u: Float) -> usize {
         let n = self.func.len();
+        if n == 0 {
+            return 0
+        }
         let mut size = n;
         let mut first = 0usize;
         while size > 0 {

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)]
@@ -165,7 +166,7 @@ impl TransformGeneric<Float> {
                 *t -= dt;
             }
         }
-        Ray::new(o.into(), r.d, Some(r.time), &*r.medium)
+        Ray::new(o.into(), r.d, Some(r.time), r.medium)
     }
 
     pub fn apply_to_interval(&self, pi: &Point3fi) -> Point3fi {
@@ -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 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;
-        let wp = self.m[3][0] * x + self.m[3][1] * y + self.m[3][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 + self.m[1][3];
+        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 + 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. {
@@ -366,7 +391,7 @@ impl TransformGeneric<Float> {
                 t = t_max - dt;
             }
         }
-        (Ray::new(Point3f::from(o), d, Some(r.time), &*r.medium), t)
+        (Ray::new(Point3f::from(o), d, Some(r.time), r.medium), t)
     }
 
     pub fn to_quaternion(self) -> Quaternion {
@@ -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()?;
-    Some(TransformGeneric::new(camera_from_world, world_from_camera))
+    let camera_from_world = world_from_camera.inverse().context("Failed to inverse viewing matrix")?;
+    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, 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 crate::Arena;
-use shared::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 {
@@ -26,7 +17,6 @@ struct BVHSplitBucket {
     pub bounds: Bounds3f,
 }
 
-
 #[derive(Debug, Clone, Copy, Default)]
 struct MortonPrimitive {
     primitive_index: usize,
@@ -113,51 +103,27 @@ impl BVHBuildNode {
     }
 }
 
-pub struct SharedPrimitiveBuffer<'a, P> {
-    ptr: *mut P,
-    pub offset: &'a AtomicUsize,
-    _marker: std::marker::PhantomData<&'a mut [P]>,
+pub trait CreateBVH {
+    fn new(primitives: Vec<Primitive>, max_prims_in_node: usize, split_method: SplitMethod) -> Self;
+    fn build_hlbvh(
+        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> {}
-unsafe impl<'a, P> Send for SharedPrimitiveBuffer<'a, P> {}
-
-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>,
+impl CreateBVH for BVHAggregate {
+    fn new(
+        mut primitives: Vec<Primitive>,
         max_prims_in_node: usize,
         split_method: SplitMethod,
     ) -> Self {
@@ -165,10 +131,11 @@ impl<P: PrimitiveTrait + Clone + Send + Sync> BVHAggregate<P> {
 
         if primitives.is_empty() {
             return Self {
-                max_prims_in_node,
-                primitives,
+                max_prims_in_node: max_prims_in_node.try_into().unwrap(),
+                node_count: 0,
+                primitives: gvec_from_slice(&primitives),
                 split_method,
-                nodes: Vec::new(),
+                nodes: gvec(),
             };
         }
 
@@ -178,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 {
@@ -190,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());
-        Self::leaf_order(&root, &mut leaf_order);
-        Self::reorder(&mut primitives, &leaf_order);
-        drop(leaf_order);
+        let mut leaf_vec = Vec::with_capacity(primitives.len());
+        leaf_order(&root, &mut leaf_vec);
+        reorder(&mut primitives, &leaf_vec);
+        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,
-            primitives,
+            node_count: node_count.try_into().unwrap(),
+            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]) {
-        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(
+    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
@@ -587,329 +476,252 @@ impl<P: PrimitiveTrait + Clone + Send + Sync> BVHAggregate<P> {
             ))
         }
     }
+}
 
-    fn build_recursive(
-        bvh_primitives: &mut [BVHPrimitiveInfo],
-        total_nodes: &AtomicUsize,
-        original_primitives: &[P],
-        max_prims_in_node: usize,
-        split_method: SplitMethod,
-    ) -> Box<BVHBuildNode> {
-        total_nodes.fetch_add(1, AtomicOrdering::Relaxed);
-        let bounds = bvh_primitives
-            .iter()
-            .fold(Bounds3f::default(), |b, p| b.union(p.bounds));
+fn build_recursive(
+    bvh_primitives: &mut [BVHPrimitiveInfo],
+    total_nodes: &AtomicUsize,
+    original_primitives: &[Primitive],
+    max_prims_in_node: usize,
+    split_method: SplitMethod,
+) -> Box<BVHBuildNode> {
+    total_nodes.fetch_add(1, AtomicOrdering::Relaxed);
+    let bounds = bvh_primitives
+        .iter()
+        .fold(Bounds3f::default(), |b, p| b.union(p.bounds));
 
-        let n_primitives = bvh_primitives.len();
-        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();
-            return Box::new(BVHBuildNode::new_leaf(n_primitives, bounds, indices));
-        }
-        let centroid_bounds = bvh_primitives.iter().fold(Bounds3f::default(), |b, p| {
-            b.union_point(p.bounds.centroid())
-        });
+    let n_primitives = bvh_primitives.len();
+    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();
+        return Box::new(BVHBuildNode::new_leaf(n_primitives, bounds, indices));
+    }
+    let centroid_bounds = bvh_primitives.iter().fold(Bounds3f::default(), |b, p| {
+        b.union_point(p.bounds.centroid())
+    });
 
-        let dim = centroid_bounds.max_dimension();
-        if centroid_bounds.p_max[dim] == centroid_bounds.p_min[dim] {
-            let indices: Vec<usize> = bvh_primitives.iter().map(|p| p.primitive_number).collect();
-            return Box::new(BVHBuildNode::new_leaf(n_primitives, bounds, indices));
-        }
+    let dim = centroid_bounds.max_dimension();
+    if centroid_bounds.p_max[dim] == centroid_bounds.p_min[dim] {
+        let indices: Vec<usize> = bvh_primitives.iter().map(|p| p.primitive_number).collect();
+        return Box::new(BVHBuildNode::new_leaf(n_primitives, bounds, indices));
+    }
 
-        let mut mid: usize;
-        match split_method {
-            SplitMethod::Middle => {
-                let pmid = (centroid_bounds.p_min[dim] + centroid_bounds.p_max[dim]) / 2.;
-                mid = partition_slice(bvh_primitives, |p| p.centroid[dim] < pmid);
+    let mut mid: usize;
+    match split_method {
+        SplitMethod::Middle => {
+            let pmid = (centroid_bounds.p_min[dim] + centroid_bounds.p_max[dim]) / 2.;
+            mid = partition_slice(bvh_primitives, |p| p.centroid[dim] < pmid);
 
-                if mid != 0 && mid != n_primitives {
-                } 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 => {
+            if mid != 0 && mid != n_primitives {
+            } else {
                 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 {
-                    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 {
-                    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])
+        }
+        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 {
+                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 {
+                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;
-                        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;
-                            }
-                            b <= min_cost_split_bucket
+                        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)
                         });
-                        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 build_leaf = |prims: &mut [BVHPrimitiveInfo]| -> Box<BVHBuildNode> {
-            Self::build_recursive(
-                prims,
-                total_nodes,
-                original_primitives,
-                max_prims_in_node,
-                split_method,
-            )
-        };
+    let (left_prims, right_prims) = bvh_primitives.split_at_mut(mid);
+    let build_leaf = |prims: &mut [BVHPrimitiveInfo]| -> Box<BVHBuildNode> {
+        build_recursive(
+            prims,
+            total_nodes,
+            original_primitives,
+            max_prims_in_node,
+            split_method,
+        )
+    };
 
-        let (child0, child1) = if n_primitives > 128 * 1024 {
-            rayon::join(|| build_leaf(left_prims), || build_leaf(right_prims))
-        } else {
-            (build_leaf(left_prims), build_leaf(right_prims))
-        };
+    let (child0, child1) = if n_primitives > 128 * 1024 {
+        rayon::join(|| build_leaf(left_prims), || build_leaf(right_prims))
+    } else {
+        (build_leaf(left_prims), build_leaf(right_prims))
+    };
 
-        let axis = dim as u8;
-        Box::new(BVHBuildNode::new_interior(axis, child0, child1))
+    let axis = dim as u8;
+    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> {
-        if self.nodes.is_empty() {
-            return None;
+    local_offset
+}
+
+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 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 },
-            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];
-            }
+        let mut prev = i;
+        let mut curr = order[i];
+        while curr != i {
+            primitives.swap(prev, curr);
+            done[prev] = true;
+            prev = curr;
+            curr = order[prev];
         }
-
-        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
+        done[prev] = true;
     }
 }
 
-impl BVHAggregate<Primitive> {
-    pub fn to_device(&self, arena: &mut Arena) -> DeviceBVHAggregate {
-        let (prims_ptr, _) = arena.alloc_slice(&self.primitives);
-        
-        let shared_nodes: Vec<shared::core::aggregates::LinearBVHNode> = self.nodes
-            .iter()
-            .map(|n| shared::core::aggregates::LinearBVHNode {
-                bounds: n.bounds,
-                primitives_offset: n.primitives_offset,
-                n_primitives: n.n_primitives,
-                axis: n.axis,
-                pad: 0,
-            })
-            .collect();
-        let (nodes_ptr, _) = arena.alloc_slice(&shared_nodes);
-
-        DeviceBVHAggregate {
-            max_prims_in_node: self.max_prims_in_node as u32,
-            primitives: prims_ptr,
-            primitive_count: self.primitives.len() as u32,
-            nodes: nodes_ptr,
-            node_count: self.nodes.len() as u32,
+fn leaf_order(node: &BVHBuildNode, out: &mut Vec<usize>) {
+    match node {
+        BVHBuildNode::Leaf {
+            primitive_indices, ..
+        } => {
+            out.extend_from_slice(primitive_indices);
+        }
+        BVHBuildNode::Interior { children, .. } => {
+            leaf_order(&children[0], out);
+            leaf_order(&children[1], out);
         }
     }
 }
+
+// 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/camera.rs

@@ -1,8 +1,7 @@
-use crate::core::image::ImageMetadata;
-use crate::core::image::{Image, ImageIO};
+use crate::core::image::{HostImage, ImageIO, ImageMetadata};
 use crate::globals::get_options;
 use crate::utils::read_float_file;
-use crate::{Arena, FileLoc, ParameterDictionary};
+use crate::{Arena, FileLoc, ParameterDictionary, ArenaUpload};
 use anyhow::{anyhow, Result};
 use shared::cameras::*;
 use shared::core::camera::{Camera, CameraBase, CameraTrait, CameraTransform};
@@ -232,8 +231,8 @@ impl CameraFactory for Camera {
                 }
 
                 let builtin_res = 256;
-                let rasterize = |vert: &[Point2f]| -> Image {
-                    let mut image = Image::new(
+                let rasterize = |vert: &[Point2f]| -> HostImage {
+                    let mut image = HostImage::new(
                         PixelFormat::F32,
                         Point2i::new(builtin_res, builtin_res),
                         &["Y"],
@@ -279,12 +278,12 @@ impl CameraFactory for Camera {
                 };
 
                 let aperture_name = params.get_one_string("aperture", "")?;
-                let mut aperture_image: Option<Image> = None;
+                let mut aperture_image: Option<HostImage> = None;
 
                 if !aperture_name.is_empty() {
                     match aperture_name.as_str() {
                         "gaussian" => {
-                            let mut img = Image::new(
+                            let mut img = HostImage::new(
                                 PixelFormat::F32,
                                 Point2i::new(builtin_res, builtin_res),
                                 &["Y"],
@@ -306,7 +305,7 @@ impl CameraFactory for Camera {
                             aperture_image = Some(img);
                         }
                         "square" => {
-                            let mut img = Image::new(
+                            let mut img = HostImage::new(
                                 PixelFormat::F32,
                                 Point2i::new(builtin_res, builtin_res),
                                 &["Y"],
@@ -316,7 +315,7 @@ impl CameraFactory for Camera {
                             let high = (0.75 * builtin_res as Float) as i32;
                             for y in low..high {
                                 for x in low..high {
-                                    img.set_channel(Point2i::new(x, y), 0, 4.0);
+                                    img.inner.set_channel(Point2i::new(x, y), 0, 4.0);
                                 }
                             }
                             aperture_image = Some(img);
@@ -353,9 +352,9 @@ impl CameraFactory for Camera {
                             aperture_image = Some(rasterize(&vert));
                         }
                         _ => {
-                            if let Ok(im) = Image::read(Path::new(&aperture_name), None) {
+                            if let Ok(im) = HostImage::read(Path::new(&aperture_name), None) {
                                 if im.image.n_channels() > 1 {
-                                    let mut mono = Image::new(
+                                    let mut mono = HostImage::new(
                                         PixelFormat::F32,
                                         im.image.resolution(),
                                         &["Y"],
@@ -366,7 +365,7 @@ impl CameraFactory for Camera {
                                         for x in 0..res.x() {
                                             let avg =
                                                 im.image.get_channels(Point2i::new(x, y)).average();
-                                            mono.set_channel(Point2i::new(x, y), 0, avg);
+                                            mono.inner.set_channel(Point2i::new(x, y), 0, avg);
                                         }
                                     }
                                     aperture_image = Some(mono);
@@ -383,10 +382,11 @@ impl CameraFactory for Camera {
                     &lens_params,
                     focal_distance,
                     aperture_diameter,
-                    Ptr::from(&*aperture_image.unwrap()),
+                    arena.upload(aperture_image)
+                    
                 );
 
-                arena.alloc(camera);
+                arena.alloc(camera.clone());
                 Ok(Camera::Realistic(camera))
             }
             "spherical" => {

src/core/color.rs

@@ -15,6 +15,7 @@ impl CreateRGBToSpectrumTable for RGBToSpectrumTable {
         assert_eq!(z_nodes.len(), RES as usize);
         assert_eq!(coeffs.len(), (RES * RES * RES) as usize * 3 * 3);
         Self {
+            n_nodes: z_nodes.len().try_into().unwrap(),
             z_nodes: gvec_from_slice(z_nodes),
             coeffs: gvec_from_slice(unsafe {
                 core::slice::from_raw_parts(

src/core/film.rs

@@ -15,7 +15,7 @@ use shared::spectra::{
     cie::SWATCHES_RAW, DenselySampledSpectrum, PiecewiseLinearSpectrum, RGBColorSpace,
 };
 use shared::utils::math::{linear_least_squares, SquareMatrix};
-use shared::{Float, Ptr};
+use shared::{Float, Ptr, leak};
 use std::sync::atomic::{AtomicUsize, Ordering};
 use std::sync::{Arc, LazyLock};
 
@@ -27,16 +27,21 @@ const SWATCH_REFLECTANCES: LazyLock<[Spectrum; N_SWATCH_REFLECTANCES]> = LazyLoc
     std::array::from_fn(|i| {
         let raw_data = SWATCHES_RAW[i];
         let pls = PiecewiseLinearSpectrum::from_interleaved(raw_data, false);
-        Spectrum::Piecewise(pls)
+        Spectrum::Piecewise(leak(pls))
     })
 });
 
+pub fn get_swatches() -> Arc<[Spectrum; N_SWATCH_REFLECTANCES]> {
+    Arc::new(*SWATCH_REFLECTANCES)
+}
+
 pub trait CreatePixelSensor: Sized {
     fn create(
         params: &ParameterDictionary,
         output_colorspace: Arc<RGBColorSpace>,
         exposure_time: Float,
         loc: &FileLoc,
+        arena: &Arena,
     ) -> Result<Self>;
 
     fn new(
@@ -46,12 +51,14 @@ pub trait CreatePixelSensor: Sized {
         output_colorspace: Arc<RGBColorSpace>,
         sensor_illum: Option<&Spectrum>,
         imaging_ratio: Float,
+        arena: &Arena,
     ) -> Self;
 
     fn new_with_white_balance(
         output_colorspace: &RGBColorSpace,
         sensor_illum: Option<&Spectrum>,
         imaging_ratio: Float,
+        arena: &Arena,
     ) -> Self;
 }
 
@@ -61,6 +68,7 @@ impl CreatePixelSensor for PixelSensor {
         output_colorspace: Arc<RGBColorSpace>,
         exposure_time: Float,
         loc: &FileLoc,
+        arena: &Arena,
     ) -> Result<Self>
     where
         Self: Sized,
@@ -79,8 +87,10 @@ impl CreatePixelSensor for PixelSensor {
             DenselySampledSpectrum::generate_cie_d(white_balance_temp)
         };
 
+        let d_ptr = arena.alloc(d_illum);
+
         let sensor_illum: Option<Arc<Spectrum>> = if white_balance_temp != 0. {
-            Some(Spectrum::Dense(d_illum.device()).into())
+            Some(Spectrum::Dense(d_ptr).into())
         } else {
             None
         };
@@ -90,6 +100,7 @@ impl CreatePixelSensor for PixelSensor {
                 output_colorspace.as_ref(),
                 sensor_illum.as_deref(),
                 imaging_ratio,
+                arena
             ));
         } else {
             let r_opt = get_named_spectrum(&format!("{}_r", sensor_name));
@@ -118,6 +129,7 @@ impl CreatePixelSensor for PixelSensor {
                         .expect("Sensor must have illuminant"),
                 ),
                 imaging_ratio,
+                arena
             ));
         }
     }
@@ -129,6 +141,7 @@ impl CreatePixelSensor for PixelSensor {
         output_colorspace: Arc<RGBColorSpace>,
         sensor_illum: Option<&Spectrum>,
         imaging_ratio: Float,
+        arena: &Arena,
     ) -> Self {
         let illum: &Spectrum = match sensor_illum {
             Some(arc_illum) => arc_illum,
@@ -138,17 +151,20 @@ impl CreatePixelSensor for PixelSensor {
         let r_bar = DenselySampledSpectrum::from_spectrum(r);
         let g_bar = DenselySampledSpectrum::from_spectrum(g);
         let b_bar = DenselySampledSpectrum::from_spectrum(b);
+        let r_ptr = arena.alloc(r_bar);
+        let g_ptr = arena.alloc(g_bar);
+        let b_ptr = arena.alloc(b_bar);
         let mut rgb_camera = [[0.; 3]; N_SWATCH_REFLECTANCES];
 
-        let swatches = Self::get_swatches();
+        let swatches = get_swatches();
 
         for i in 0..N_SWATCH_REFLECTANCES {
             let rgb = PixelSensor::project_reflectance::<RGB>(
                 &swatches[i],
                 illum,
-                &Spectrum::Dense(r_bar),
-                &Spectrum::Dense(g_bar),
-                &Spectrum::Dense(b_bar),
+                &Spectrum::Dense(r_ptr),
+                &Spectrum::Dense(g_ptr),
+                &Spectrum::Dense(b_ptr),
             );
             for c in 0..3 {
                 rgb_camera[i][c] = rgb[c];
@@ -157,7 +173,7 @@ impl CreatePixelSensor for PixelSensor {
 
         let mut xyz_output = [[0.; 3]; N_SWATCH_REFLECTANCES];
         let spectra = get_spectra_context();
-        let sensor_white_g = illum.inner_product(&Spectrum::Dense(g_bar));
+        let sensor_white_g = illum.inner_product(&Spectrum::Dense(g_ptr));
         let sensor_white_y = illum.inner_product(&Spectrum::Dense(spectra.y));
         for i in 0..N_SWATCH_REFLECTANCES {
             let s = swatches[i].clone();
@@ -177,9 +193,9 @@ impl CreatePixelSensor for PixelSensor {
             .expect("Could not convert sensor illuminance to XYZ space");
 
         PixelSensor {
-            r_bar: r_bar.clone(),
-            g_bar: g_bar.clone(),
-            b_bar: b_bar.clone(),
+            r_bar: r_ptr,
+            g_bar: g_ptr,
+            b_bar: b_ptr,
             imaging_ratio,
             xyz_from_sensor_rgb,
         }
@@ -189,6 +205,7 @@ impl CreatePixelSensor for PixelSensor {
         output_colorspace: &RGBColorSpace,
         sensor_illum: Option<&Spectrum>,
         imaging_ratio: Float,
+        arena: &Arena
     ) -> Self {
         let spectra = get_spectra_context();
         let r_bar = CIE_X_DATA.clone();
@@ -205,24 +222,21 @@ impl CreatePixelSensor for PixelSensor {
         }
 
         PixelSensor {
-            r_bar: r_bar.clone(),
-            g_bar: g_bar.clone(),
-            b_bar: b_bar.clone(),
+            r_bar: arena.alloc(r_bar),
+            g_bar: arena.alloc(g_bar),
+            b_bar: arena.alloc(b_bar),
             xyz_from_sensor_rgb,
             imaging_ratio,
         }
     }
 
-    fn get_swatches() -> Arc<[Spectrum; N_SWATCH_REFLECTANCES]> {
-        Arc::new(*SWATCH_REFLECTANCES)
-    }
 }
 
 pub trait CreateFilmBase {
     fn create(
         params: &ParameterDictionary,
         filter: Filter,
-        sensor: Option<&PixelSensor>,
+        sensor: Ptr<PixelSensor>,
         loc: &FileLoc,
     ) -> Result<Self>
     where
@@ -233,7 +247,7 @@ impl CreateFilmBase for FilmBase {
     fn create(
         params: &ParameterDictionary,
         filter: Filter,
-        sensor: Option<&PixelSensor>,
+        sensor: Ptr<PixelSensor>,
         loc: &FileLoc,
     ) -> Result<Self>
     where
@@ -287,7 +301,7 @@ impl CreateFilmBase for FilmBase {
             pixel_bounds,
             filter,
             diagonal: diagonal_mm * 0.001,
-            sensor: Ptr::from(sensor.unwrap()),
+            sensor,
         })
     }
 }

src/core/filter.rs

@@ -1,10 +1,10 @@
-use crate::utils::sampling::PiecewiseConstant2D;
-use crate::utils::{FileLoc, ParameterDictionary};
 use crate::Arena;
-use anyhow::{anyhow, Result};
+use crate::{FileLoc, ParameterDictionary};
+use anyhow::{bail, Result};
 use shared::core::filter::Filter;
 use shared::core::geometry::{Bounds2f, Point2f, Vector2f};
 use shared::filters::*;
+use shared::utils::sampling::PiecewiseConstant2D;
 use shared::{Array2D, Float};
 
 pub trait FilterFactory {
@@ -28,13 +28,14 @@ impl FilterFactory for Filter {
                 let xw = params.get_one_float("xradius", 0.5)?;
                 let yw = params.get_one_float("yradius", 0.5)?;
                 let filter = BoxFilter::new(Vector2f::new(xw, yw));
-                Ok(Filter::Box(filter))
+                Ok(Filter::Box(arena.alloc(filter)))
             }
             "gaussian" => {
                 let xw = params.get_one_float("xradius", 1.5)?;
                 let yw = params.get_one_float("yradius", 1.5)?;
                 let sigma = params.get_one_float("sigma", 0.5)?;
                 let filter = GaussianFilter::new(Vector2f::new(xw, yw), sigma);
+                Ok(Filter::Gaussian(arena.alloc(filter)))
             }
             "mitchell" => {
                 let xw = params.get_one_float("xradius", 2.)?;
@@ -42,22 +43,22 @@ impl FilterFactory for Filter {
                 let b = params.get_one_float("B", 1. / 3.)?;
                 let c = params.get_one_float("C", 1. / 3.)?;
                 let filter = MitchellFilter::new(Vector2f::new(xw, yw), b, c);
-                Ok(Filter::Mitchell(filter))
+                Ok(Filter::Mitchell(arena.alloc(filter)))
             }
             "sinc" => {
                 let xw = params.get_one_float("xradius", 4.)?;
                 let yw = params.get_one_float("yradius", 4.)?;
                 let tau = params.get_one_float("tau", 3.)?;
                 let filter = LanczosSincFilter::new(Vector2f::new(xw, yw), tau);
-                Ok(Filter::LanczosSinc(filter))
+                Ok(Filter::LanczosSinc(arena.alloc(filter)))
             }
             "triangle" => {
                 let xw = params.get_one_float("xradius", 2.)?;
                 let yw = params.get_one_float("yradius", 2.)?;
                 let filter = TriangleFilter::new(Vector2f::new(xw, yw));
-                Ok(Filter::Triangle(filter))
+                Ok(Filter::Triangle(arena.alloc(filter)))
             }
-            _ => Err(anyhow!("Film type '{}' unknown at {}", name, loc)),
+            _ => bail!("Film type '{}' unknown at {}", name, loc),
         }
     }
 }

src/core/image/io.rs

@@ -1,15 +1,14 @@
-use super::{Image, ImageAndMetadata, ImageMetadata};
-use crate::core::image::{PixelStorage, WrapMode};
-use crate::utils::error::ImageError;
+use super::{HostImage, ImageAndMetadata, ImageMetadata};
+use crate::core::image::WrapMode;
 use anyhow::{Context, Result, bail};
 use exr::prelude::{read_first_rgba_layer_from_file, write_rgba_file};
-use image_rs::{DynamicImage, ImageReader};
 use shared::Float;
 use shared::core::color::{ColorEncoding, LINEAR, SRGB};
 use shared::core::geometry::Point2i;
 use shared::core::image::PixelFormat;
 use std::fs::File;
 use std::io::{BufRead, BufReader, BufWriter, Read, Write};
+use image_rs::{DynamicImage, ImageReader};
 use std::path::Path;
 
 pub trait ImageIO {
@@ -22,7 +21,7 @@ pub trait ImageIO {
     fn to_u8_buffer(&self) -> Vec<u8>;
 }
 
-impl ImageIO for Image {
+impl ImageIO for HostImage {
     fn read(path: &Path, encoding: Option<ColorEncoding>) -> Result<ImageAndMetadata> {
         let ext = path
             .extension()
@@ -40,55 +39,39 @@ impl ImageIO for Image {
     fn write(&self, filename: &str, metadata: &ImageMetadata) -> Result<()> {
         let path = Path::new(filename);
         let ext = path.extension().and_then(|s| s.to_str()).unwrap_or("");
-        let res = match ext.to_lowercase().as_str() {
+        match ext.to_lowercase().as_str() {
             "exr" => self.write_exr(path, metadata),
             "png" => self.write_png(path),
             "pfm" => self.write_pfm(path),
             "qoi" => self.write_qoi(path),
             _ => Err(anyhow::anyhow!("Unsupported write format: {}", ext)),
-        };
-        res.map_err(|e| ImageError::Io(std::io::Error::other(e)))?;
-        Ok(())
+        }
     }
 
     fn write_png(&self, path: &Path) -> Result<()> {
-        let w = self.resolution().x() as u32;
-        let h = self.resolution().y() as u32;
-
-        // Convert whatever we have to u8 [0..255]
+        let w = self.inner.resolution().x() as u32;
+        let h = self.inner.resolution().y() as u32;
         let data = self.to_u8_buffer();
-        let channels = self.n_channels();
+        let channels = self.inner.n_channels();
 
         match channels {
             1 => {
-                // Luma
                 image_rs::save_buffer_with_format(
-                    path,
-                    &data,
-                    w,
-                    h,
+                    path, &data, w, h,
                     image_rs::ColorType::L8,
                     image_rs::ImageFormat::Png,
                 )?;
             }
             3 => {
-                // RGB
                 image_rs::save_buffer_with_format(
-                    path,
-                    &data,
-                    w,
-                    h,
+                    path, &data, w, h,
                     image_rs::ColorType::Rgb8,
                     image_rs::ImageFormat::Png,
                 )?;
             }
             4 => {
-                // RGBA
                 image_rs::save_buffer_with_format(
-                    path,
-                    &data,
-                    w,
-                    h,
+                    path, &data, w, h,
                     image_rs::ColorType::Rgba8,
                     image_rs::ImageFormat::Png,
                 )?;
@@ -99,37 +82,30 @@ impl ImageIO for Image {
     }
 
     fn write_qoi(&self, path: &Path) -> Result<()> {
-        let w = self.resolution().x() as u32;
-        let h = self.resolution().y() as u32;
+        let w = self.inner.resolution().x() as u32;
+        let h = self.inner.resolution().y() as u32;
         let data = self.to_u8_buffer();
 
-        let color_type = match self.n_channels() {
+        let color_type = match self.inner.n_channels() {
             3 => image_rs::ColorType::Rgb8,
             4 => image_rs::ColorType::Rgba8,
             _ => bail!("QOI only supports 3 or 4 channels"),
         };
 
         image_rs::save_buffer_with_format(
-            path,
-            &data,
-            w,
-            h,
-            color_type,
-            image_rs::ImageFormat::Qoi,
+            path, &data, w, h, color_type, image_rs::ImageFormat::Qoi,
         )?;
         Ok(())
     }
 
     fn write_exr(&self, path: &Path, _metadata: &ImageMetadata) -> Result<()> {
-        // EXR requires F32
-        let w = self.resolution().x() as usize;
-        let h = self.resolution().y() as usize;
-        let c = self.n_channels();
+        let w = self.inner.resolution().x() as usize;
+        let h = self.inner.resolution().y() as usize;
+        let c = self.inner.n_channels();
 
         write_rgba_file(path, w, h, |x, y| {
-            // Helper to get float value regardless of internal storage
             let get = |ch| {
-                self.get_channel_with_wrap(
+                self.inner.get_channel_with_wrap(
                     Point2i::new(x as i32, y as i32),
                     ch,
                     WrapMode::Clamp.into(),
@@ -154,27 +130,22 @@ impl ImageIO for Image {
         let file = File::create(path)?;
         let mut writer = BufWriter::new(file);
 
-        if self.n_channels() != 3 {
+        if self.inner.n_channels() != 3 {
             bail!("PFM writing currently only supports 3 channels (RGB)");
         }
 
-        // Header
-        let res = self.resolution();
+        let res = self.inner.resolution();
         writeln!(writer, "PF")?;
         writeln!(writer, "{} {}", res.x(), res.y())?;
-        let scale = if cfg!(target_endian = "little") {
-            -1.0
-        } else {
-            1.0
-        };
+        let scale = if cfg!(target_endian = "little") { -1.0 } else { 1.0 };
         writeln!(writer, "{}", scale)?;
 
-        // PBRT stores top-to-bottom.
         for y in (0..res.y()).rev() {
             for x in 0..res.x() {
                 for c in 0..3 {
-                    let val =
-                        self.get_channel_with_wrap(Point2i::new(x, y), c, WrapMode::Clamp.into());
+                    let val = self.inner.get_channel_with_wrap(
+                        Point2i::new(x, y), c, WrapMode::Clamp.into(),
+                    );
                     writer.write_all(&val.to_le_bytes())?;
                 }
             }
@@ -183,20 +154,18 @@ impl ImageIO for Image {
         Ok(())
     }
 
-    // TODO: Change Image to use Vec for data, always. Only convert to Device types on
-    // constructors/creation
     fn to_u8_buffer(&self) -> Vec<u8> {
-        match &self.pixels {
-            PixelStorage::U8(data) => data.to_vec(),
-            PixelStorage::F16(data) => data
-                .iter()
-                .map(|v| (v.to_f32().clamp(0.0, 1.0) * 255.0 + 0.5) as u8)
-                .collect(),
-            PixelStorage::F32(data) => data
-                .iter()
-                .map(|v| (v.clamp(0.0, 1.0) * 255.0 + 0.5) as u8)
-                .collect(),
+        let res = self.inner.resolution();
+        let n_pixels = (res.x() * res.y()) as usize;
+        let nc = self.inner.n_channels() as usize;
+        let total = n_pixels * nc;
+
+        let mut buf = Vec::with_capacity(total);
+        for i in 0..total {
+            let val = unsafe { self.inner.pixels.read(i, &self.inner.encoding) };
+            buf.push((val.clamp(0.0, 1.0) * 255.0 + 0.5) as u8);
         }
+        buf
     }
 }
 
@@ -209,21 +178,24 @@ fn read_generic(path: &Path, encoding: Option<ColorEncoding>) -> Result<ImageAnd
     let h = dyn_img.height() as i32;
     let res = Point2i::new(w, h);
 
-    // Check if it was loaded as high precision or standard
-    let rgb_names = vec!["R", "G", "B"];
-    let rgba_names = vec!["R", "G", "B", "A"];
+    let rgb_names: &[&str] = &["R", "G", "B"];
+    let rgba_names: &[&str] = &["R", "G", "B", "A"];
+
     let image = match dyn_img {
-        DynamicImage::ImageRgb32F(buf) => Image::from_f32(buf.into_raw(), res, &rgb_names),
-        DynamicImage::ImageRgba32F(buf) => Image::from_f32(buf.into_raw(), res, &rgba_names),
+        DynamicImage::ImageRgb32F(buf) => {
+            HostImage::from_f32(&buf.into_raw(), res, rgb_names)
+        }
+        DynamicImage::ImageRgba32F(buf) => {
+            HostImage::from_f32(&buf.into_raw(), res, rgba_names)
+        }
         _ => {
-            // Default to RGB8 for everything else
             let enc = encoding.unwrap_or(SRGB);
             if dyn_img.color().has_alpha() {
                 let buf = dyn_img.to_rgba8();
-                Image::from_u8(buf.into_raw(), res, &rgba_names, enc)
+                HostImage::from_u8(&buf.into_raw(), res, rgba_names, enc)
             } else {
                 let buf = dyn_img.to_rgb8();
-                Image::from_u8(buf.into_raw(), res, &rgb_names, enc)
+                HostImage::from_u8(&buf.into_raw(), res, rgb_names, enc)
             }
         }
     };
@@ -242,7 +214,6 @@ fn read_exr(path: &Path) -> Result<ImageAndMetadata> {
         |buffer, position, pixel| {
             let width = position.width();
             let idx = (position.y() * width + position.x()) * 4;
-            // Map exr pixel struct to our buffer
             buffer[idx] = pixel.0;
             buffer[idx + 1] = pixel.1;
             buffer[idx + 2] = pixel.2;
@@ -254,10 +225,11 @@ fn read_exr(path: &Path) -> Result<ImageAndMetadata> {
     let w = image.layer_data.size.width() as i32;
     let h = image.layer_data.size.height() as i32;
 
-    let image = Image::from_f32(
-        image.layer_data.channel_data.pixels,
+    let rgba_names: &[&str] = &["R", "G", "B", "A"];
+    let image = HostImage::from_f32(
+        &image.layer_data.channel_data.pixels,
         Point2i::new(w, h),
-        &vec!["R", "G", "B", "A"],
+        rgba_names,
     );
 
     let metadata = ImageMetadata::default();
@@ -268,7 +240,6 @@ fn read_pfm(path: &Path) -> Result<ImageAndMetadata> {
     let file = File::open(path)?;
     let mut reader = BufReader::new(file);
 
-    // PFM Headers are: "PF\nwidth height\nscale\n" (or Pf for grayscale)
     let mut header_word = String::new();
     reader.read_line(&mut header_word)?;
     let header_word = header_word.trim();
@@ -309,9 +280,7 @@ fn read_pfm(path: &Path) -> Result<ImageAndMetadata> {
 
     let mut pixels = vec![0.0 as Float; (w * h * channels) as usize];
 
-    // PFM is Bottom-to-Top
     for y in 0..h {
-        // Flippety-do
         let src_y = h - 1 - y;
         for x in 0..w {
             for c in 0..channels {
@@ -331,13 +300,9 @@ fn read_pfm(path: &Path) -> Result<ImageAndMetadata> {
         }
     }
 
-    let names = if channels == 1 {
-        vec!["Y"]
-    } else {
-        vec!["R", "G", "B"]
-    };
+    let names: &[&str] = if channels == 1 { &["Y"] } else { &["R", "G", "B"] };
 
-    let image = Image::new(PixelFormat::F32, Point2i::new(w, h), &names, LINEAR.into());
+    let image = HostImage::new(PixelFormat::F32, Point2i::new(w, h), names, LINEAR.into());
     let metadata = ImageMetadata::default();
 
     Ok(ImageAndMetadata { image, metadata })

src/core/image/mod.rs

@@ -1,14 +1,12 @@
-use crate::utils::containers::Array2D;
-use anyhow::{Result, anyhow};
+use anyhow::{anyhow, Result};
 use half::f16;
 use rayon::prelude::{IndexedParallelIterator, ParallelIterator, ParallelSliceMut};
-use shared::Float;
-use shared::Ptr;
 use shared::core::color::{ColorEncoding, ColorEncodingTrait, LINEAR};
 use shared::core::geometry::{Bounds2f, Point2f, Point2i};
 use shared::core::image::{Image, ImageBase, PixelFormat, Pixels, WrapMode, WrapMode2D};
 use shared::utils::math::square;
-use smallvec::{SmallVec, smallvec};
+use shared::{Array2D, Float, Ptr};
+use smallvec::{smallvec, SmallVec};
 use std::ops::{Deref, DerefMut};
 use std::sync::Arc;
 
@@ -90,31 +88,32 @@ impl HostImage {
         let n_channels = channel_names.len() as i32;
         Self {
             inner: Image::from_u8(data, resolution, n_channels, encoding),
-            channel_names: channel_names.iter().map(|s| s.as_ref().to_string()).collect(),
+            channel_names: channel_names
+                .iter()
+                .map(|s| s.as_ref().to_string())
+                .collect(),
         }
     }
 
-    pub fn from_f32(
-        data: &[f32],
-        resolution: Point2i,
-        channel_names: &[impl AsRef<str>],
-    ) -> Self {
+    pub fn from_f32(data: &[f32], resolution: Point2i, channel_names: &[impl AsRef<str>]) -> Self {
         let n_channels = channel_names.len() as i32;
         Self {
             inner: Image::from_f32(data, resolution, n_channels),
-            channel_names: channel_names.iter().map(|s| s.as_ref().to_string()).collect(),
+            channel_names: channel_names
+                .iter()
+                .map(|s| s.as_ref().to_string())
+                .collect(),
         }
     }
 
-    pub fn from_f16(
-        data: &[u16],
-        resolution: Point2i,
-        channel_names: &[impl AsRef<str>],
-    ) -> Self {
+    pub fn from_f16(data: &[u16], resolution: Point2i, channel_names: &[impl AsRef<str>]) -> Self {
         let n_channels = channel_names.len() as i32;
         Self {
             inner: Image::from_f16(data, resolution, n_channels),
-            channel_names: channel_names.iter().map(|s| s.as_ref().to_string()).collect(),
+            channel_names: channel_names
+                .iter()
+                .map(|s| s.as_ref().to_string())
+                .collect(),
         }
     }
 
@@ -127,7 +126,10 @@ impl HostImage {
         let n_channels = channel_names.len() as i32;
         Self {
             inner: Image::new(format, resolution, n_channels, encoding),
-            channel_names: channel_names.iter().map(|s| s.as_ref().to_string()).collect(),
+            channel_names: channel_names
+                .iter()
+                .map(|s| s.as_ref().to_string())
+                .collect(),
         }
     }
 
@@ -137,7 +139,11 @@ impl HostImage {
         values: &[f32],
     ) -> Self {
         let n_channels = channel_names.len();
-        assert_eq!(values.len(), n_channels, "values length must match channel count");
+        assert_eq!(
+            values.len(),
+            n_channels,
+            "values length must match channel count"
+        );
 
         let n_pixels = (resolution.x() * resolution.y()) as usize;
         let mut data = Vec::with_capacity(n_pixels * n_channels);
@@ -179,14 +185,37 @@ impl HostImage {
     }
 
     pub fn channel_names_from_desc(&self, desc: &ImageChannelDesc) -> Vec<&str> {
-        desc.offset.iter().map(|&i| self.channel_names[i].as_str()).collect()
+        desc.offset
+            .iter()
+            .map(|&i| self.channel_names[i].as_str())
+            .collect()
+    }
+
+    pub fn bilerp_channel_with_wrap(&self, p: Point2f, c: i32, wrap: WrapMode2D) -> Float {
+        self.inner.bilerp_channel_with_wrap(p, c, wrap)
+    }
+
+    pub fn bilerp_channel(&self, p: Point2f, c: i32) -> Float {
+        self.bilerp_channel_with_wrap(p, c, WrapMode::Clamp.into())
+    }
+
+    pub fn get_channel(&self, p: Point2i, c: i32) -> Float {
+        self.get_channel_with_wrap(p, c, WrapMode::Clamp.into())
+    }
+
+    pub fn get_channel_with_wrap(&self, p: Point2i, c: i32, wrap_mode: WrapMode2D) -> Float {
+        self.inner.get_channel_with_wrap(p, c, wrap_mode)
     }
 
     pub fn get_channels(&self, p: Point2i) -> ImageChannelValues {
         self.get_channels_with_wrap(p, WrapMode::Clamp.into())
     }
 
-    pub fn get_channels_with_wrap(&self, mut p: Point2i, wrap_mode: WrapMode2D) -> ImageChannelValues {
+    pub fn get_channels_with_wrap(
+        &self,
+        mut p: Point2i,
+        wrap_mode: WrapMode2D,
+    ) -> ImageChannelValues {
         if !self.inner.remap_pixel_coords(&mut p, wrap_mode) {
             return ImageChannelValues(SmallVec::from_elem(0.0, self.inner.n_channels() as usize));
         }
@@ -212,7 +241,11 @@ impl HostImage {
         let pixel_offset = self.inner.pixel_offset(pp);
         let mut values = SmallVec::with_capacity(desc.offset.len());
         for &c in &desc.offset {
-            values.push(unsafe { self.inner.pixels.read(pixel_offset + c, &self.inner.encoding) });
+            values.push(unsafe {
+                self.inner
+                    .pixels
+                    .read(pixel_offset + c, &self.inner.encoding)
+            });
         }
         ImageChannelValues(values)
     }
@@ -223,20 +256,29 @@ impl HostImage {
         }
     }
 
+    pub fn set_channel(&mut self, p: Point2i, c: i32, val: Float) {
+        self.inner.set_channel(p, c, val);
+    }
+
     pub fn select_channels(&self, desc: &ImageChannelDesc) -> Self {
-        let new_names: Vec<String> = desc.offset.iter().map(|&i| self.channel_names[i].clone()).collect();
+        let new_names: Vec<String> = desc
+            .offset
+            .iter()
+            .map(|&i| self.channel_names[i].clone())
+            .collect();
         let res = self.inner.resolution();
         let pixel_count = (res.x() * res.y()) as usize;
         let src_nc = self.inner.n_channels() as usize;
         let dst_nc = desc.offset.len();
 
-        // Always produce f32 output for simplicity
         let mut dst = vec![0.0f32; pixel_count * dst_nc];
         for i in 0..pixel_count {
             let src_offset = i * src_nc;
             for (out_idx, &in_c) in desc.offset.iter().enumerate() {
                 dst[i * dst_nc + out_idx] = unsafe {
-                    self.inner.pixels.read(src_offset + in_c, &self.inner.encoding)
+                    self.inner
+                        .pixels
+                        .read(src_offset + in_c, &self.inner.encoding)
                 };
             }
         }
@@ -256,7 +298,7 @@ impl HostImage {
         dist.as_mut_slice()
             .par_chunks_mut(width as usize)
             .enumerate()
-            .for_each(|(y, row)| {
+            .for_each(|(y, row): (usize, &mut [Float])| {
                 let y = y as i32;
                 for (x, out_val) in row.iter_mut().enumerate() {
                     let x = x as i32;
@@ -310,15 +352,18 @@ impl HostImage {
 
         for y in 0..res.y() {
             for x in 0..res.x() {
-                let v = self.get_channels_with_desc(
-                    Point2i::new(x, y), desc, WrapMode::Clamp.into(),
-                );
+                let v =
+                    self.get_channels_with_desc(Point2i::new(x, y), desc, WrapMode::Clamp.into());
                 let v_ref = ref_img.get_channels_with_desc(
-                    Point2i::new(x, y), &ref_desc, WrapMode::Clamp.into(),
+                    Point2i::new(x, y),
+                    &ref_desc,
+                    WrapMode::Clamp.into(),
                 );
                 for c in 0..n_channels {
                     let se = square(v[c] as f64 - v_ref[c] as f64);
-                    if se.is_infinite() { continue; }
+                    if se.is_infinite() {
+                        continue;
+                    }
                     sum_se[c] += se;
                     if generate_mse_image {
                         let idx = (y as usize * width + x as usize) * n_channels + c;

src/core/image/ops.rs

@@ -1,6 +1,4 @@
 use super::HostImage;
-use crate::core::image::pixel::PixelStorageTrait;
-use crate::core::image::PixelStorage;
 use rayon::prelude::*;
 use shared::core::color::ColorEncoding;
 use shared::core::geometry::{Bounds2i, Point2i};
@@ -18,98 +16,84 @@ pub struct ResampleWeight {
 impl HostImage {
     pub fn flip_y(&mut self) {
         let res = self.inner.resolution;
-        let nc = self.inner.n_channels as usize;
+        let nc = self.inner.n_channels;
 
-        match self.inner.format {
-            PixelFormat::U8 => flip_y_kernel(self.inner.pixels.as_u8_mut(), res, nc),
-            PixelFormat::F16 => flip_y_kernel(self.inner.pixels.as_f16_mut(), res, nc),
-            PixelFormat::F32 => flip_y_kernel(self.inner.pixels.as_f32_slice_mut(), res, nc),
+        for y in 0..res.y() / 2 {
+            let y2 = res.y() - 1 - y;
+            for x in 0..res.x() {
+                for c in 0..nc {
+                    let a = self.inner.get_channel(Point2i::new(x, y), c);
+                    let b = self.inner.get_channel(Point2i::new(x, y2), c);
+                    self.inner.set_channel(Point2i::new(x, y), c, b);
+                    self.inner.set_channel(Point2i::new(x, y2), c, a);
+                }
+            }
         }
     }
 
     pub fn crop(&self, bounds: Bounds2i) -> HostImage {
-        let n_channels = self.inner.n_channels as usize;
+        let nc = self.inner.n_channels;
         let new_res = Point2i::new(
             bounds.p_max.x() - bounds.p_min.x(),
             bounds.p_max.y() - bounds.p_min.y(),
         );
 
-        let new_names = self.channel_names.clone();
         let mut new_image =
-            HostImage::new(self.inner.format, new_res, &new_names, self.inner.encoding);
+            HostImage::new(self.inner.format, new_res, &self.channel_names, self.inner.encoding);
 
-        match self.inner.format {
-            PixelFormat::U8 => crop_kernel(
-                self.inner.pixels.as_u8(),
-                new_image.inner.pixels.as_u8_mut(),
-                self.inner.resolution,
-                bounds,
-                n_channels,
-            ),
-            PixelFormat::F16 => crop_kernel(
-                self.inner.pixels.as_f16(),
-                new_image.inner.pixels.as_f16_mut(),
-                self.inner.resolution,
-                bounds,
-                n_channels,
-            ),
-            PixelFormat::F32 => crop_kernel(
-                self.inner.pixels.as_f32_slice(),
-                new_image.inner.pixels.as_f32_slice_mut(),
-                self.inner.resolution,
-                bounds,
-                n_channels,
-            ),
+        for y in 0..new_res.y() {
+            for x in 0..new_res.x() {
+                let src = Point2i::new(bounds.p_min.x() + x, bounds.p_min.y() + y);
+                let dst = Point2i::new(x, y);
+                for c in 0..nc {
+                    let val = self.inner.get_channel(src, c);
+                    new_image.inner.set_channel(dst, c, val);
+                }
+            }
         }
 
         new_image
     }
 
     pub fn copy_rect_out(&self, extent: Bounds2i, buf: &mut [Float], wrap: WrapMode2D) {
-        match self.inner.format {
-            PixelFormat::U8 => {
-                copy_rect_out_kernel(self.inner.pixels.as_u8(), self, extent, buf, wrap)
-            }
-            PixelFormat::F16 => {
-                copy_rect_out_kernel(self.inner.pixels.as_f16(), self, extent, buf, wrap)
-            }
-            PixelFormat::F32 => {
-                copy_rect_out_kernel(self.inner.pixels.as_f32_slice(), self, extent, buf, wrap)
-            }
-        }
+        let w = (extent.p_max.x() - extent.p_min.x()) as usize;
+        let channels = self.inner.n_channels as usize;
+
+        buf.par_chunks_mut(w * channels)
+            .enumerate()
+            .for_each(|(y_rel, row_buf)| {
+                let y = extent.p_min.y() + y_rel as i32;
+                for x_rel in 0..w {
+                    let x = extent.p_min.x() + x_rel as i32;
+                    let p = Point2i::new(x, y);
+                    for c in 0..channels {
+                        row_buf[x_rel * channels + c] =
+                            self.inner.get_channel_with_wrap(p, c as i32, wrap);
+                    }
+                }
+            });
     }
 
     pub fn copy_rect_in(&mut self, extent: Bounds2i, buf: &[Float]) {
-        let res = self.inner.resolution;
-        let n_channels = self.inner.n_channels as usize;
-        let encoding = self.inner.encoding;
-        let format = self.inner.format;
+        let w = (extent.p_max.x() - extent.p_min.x()) as usize;
+        let channels = self.inner.n_channels as usize;
 
-        match format {
-            PixelFormat::U8 => copy_rect_in_kernel(
-                self.inner.pixels.as_u8_mut(),
-                res,
-                n_channels,
-                encoding,
-                extent,
-                buf,
-            ),
-            PixelFormat::F16 => copy_rect_in_kernel(
-                self.inner.pixels.as_f16_mut(),
-                res,
-                n_channels,
-                encoding,
-                extent,
-                buf,
-            ),
-            PixelFormat::F32 => copy_rect_in_kernel(
-                self.inner.pixels.as_f32_slice_mut(),
-                res,
-                n_channels,
-                encoding,
-                extent,
-                buf,
-            ),
+        for (y_rel, row) in buf.chunks(w * channels).enumerate() {
+            let y = extent.p_min.y() + y_rel as i32;
+            if y < 0 || y >= self.inner.resolution.y() {
+                continue;
+            }
+
+            for x_rel in 0..w {
+                let x = extent.p_min.x() + x_rel as i32;
+                if x < 0 || x >= self.inner.resolution.x() {
+                    continue;
+                }
+                for c in 0..channels {
+                    let val = row[x_rel * channels + c];
+                    self.inner.set_channel(Point2i::new(x, y), c as i32, val);
+                }
+            }
         }
     }
 
@@ -180,6 +164,7 @@ impl HostImage {
             }
 
             let new_res = Point2i::new((old.x() / 2).max(1), (old.y() / 2).max(1));
+            let nc = prev.n_channels();
             let mut next = HostImage::new(
                 prev.inner.format,
                 new_res,
@@ -187,177 +172,37 @@ impl HostImage {
                 prev.inner.encoding,
             );
 
-            match next.inner.format {
-                PixelFormat::U8 => downsample_kernel(
-                    next.inner.pixels.as_u8_mut(),
-                    new_res,
-                    &prev.inner,
-                    internal_wrap,
-                ),
-                PixelFormat::F16 => downsample_kernel(
-                    next.inner.pixels.as_f16_mut(),
-                    new_res,
-                    &prev.inner,
-                    internal_wrap,
-                ),
-                PixelFormat::F32 => downsample_kernel(
-                    next.inner.pixels.as_f32_slice_mut(),
-                    new_res,
-                    &prev.inner,
-                    internal_wrap,
-                ),
+            for y in 0..new_res.y() {
+                for x in 0..new_res.x() {
+                    let src_x = x * 2;
+                    let src_y = y * 2;
+                    for c in 0..nc {
+                        let mut sum = 0.0;
+                        let mut count = 0.0;
+                        for dy in 0..2i32 {
+                            for dx in 0..2i32 {
+                                let sx = src_x + dx;
+                                let sy = src_y + dy;
+                                if sx < old.x() && sy < old.y() {
+                                    sum += prev.inner.get_channel_with_wrap(
+                                        Point2i::new(sx, sy), c, internal_wrap,
+                                    );
+                                    count += 1.0;
+                                }
+                            }
+                        }
+                        let avg = if count > 0.0 { sum / count } else { 0.0 };
+                        next.inner.set_channel(Point2i::new(x, y), c, avg);
+                    }
+                }
             }
+
             levels.push(next);
         }
         levels
     }
 }
 
-fn flip_y_kernel<T: PixelStorageTrait>(pixels: &mut [T], res: Point2i, channels: usize) {
-    let w = res.x() as usize;
-    let h = res.y() as usize;
-    let stride = w * channels;
-
-    for y in 0..(h / 2) {
-        let bot = h - 1 - y;
-        for i in 0..stride {
-            pixels.swap(y * stride + i, bot * stride + i);
-        }
-    }
-}
-
-fn crop_kernel<T: PixelStorageTrait>(
-    src: &[T],
-    dst: &mut [T],
-    src_res: Point2i,
-    bounds: Bounds2i,
-    channels: usize,
-) {
-    let dst_w = (bounds.p_max.x() - bounds.p_min.x()) as usize;
-    // let dst_h = (bounds.p_max.y() - bounds.p_min.y()) as usize;
-
-    dst.par_chunks_mut(dst_w * channels)
-        .enumerate()
-        .for_each(|(dy, dst_row)| {
-            let sy = bounds.p_min.y() as usize + dy;
-            let sx_start = bounds.p_min.x() as usize;
-
-            let src_offset = (sy * src_res.x() as usize + sx_start) * channels;
-            let count = dst_w * channels;
-
-            dst_row.copy_from_slice(&src[src_offset..src_offset + count]);
-        });
-}
-
-fn copy_rect_out_kernel<T: PixelStorageTrait>(
-    src: &[T],
-    image: &HostImage,
-    extent: Bounds2i,
-    buf: &mut [Float],
-    wrap: WrapMode2D,
-) {
-    let w = (extent.p_max.x() - extent.p_min.x()) as usize;
-    let channels = image.n_channels() as usize;
-    let enc = image.encoding();
-    let res = image.resolution();
-
-    buf.par_chunks_mut(w * channels)
-        .enumerate()
-        .for_each(|(y_rel, row_buf)| {
-            let y = extent.p_min.y() + y_rel as i32;
-            for x_rel in 0..w {
-                let x = extent.p_min.x() + x_rel as i32;
-
-                if x >= 0 && x < res.x() && y >= 0 && y < res.y() {
-                    let offset = (y as usize * res.x() as usize + x as usize) * channels;
-
-                    for c in 0..channels {
-                        row_buf[x_rel * channels + c] = T::to_linear(src[offset + c], enc);
-                    }
-                } else {
-                    // We fall back to get_channel which handles the wrapping math.
-                    let p = Point2i::new(x, y);
-                    for c in 0..channels {
-                        row_buf[x_rel * channels + c] =
-                            image.get_channel_with_wrap(p, c.try_into().unwrap(), wrap);
-                    }
-                }
-            }
-        });
-}
-
-fn copy_rect_in_kernel<T: PixelStorageTrait>(
-    dst: &mut [T],
-    res: Point2i,
-    channels: usize,
-    enc: ColorEncoding,
-    extent: Bounds2i,
-    buf: &[Float],
-) {
-    let w = (extent.p_max.x() - extent.p_min.x()) as usize;
-    let res_x = res.x() as usize;
-
-    let rows = buf.chunks(w * channels);
-    for (y_rel, row) in rows.enumerate() {
-        let y = extent.p_min.y() + y_rel as i32;
-        if y < 0 || y >= res.y() {
-            continue;
-        }
-
-        let dst_row_start = (y as usize * res_x) * channels;
-
-        for (x_rel, &val) in row.iter().enumerate() {
-            let c = x_rel % channels;
-            let x_pixel = x_rel / channels;
-            let x = extent.p_min.x() + x_pixel as i32;
-
-            if x >= 0 && x < res.x() {
-                let idx = dst_row_start + (x as usize * channels) + c;
-                dst[idx] = T::from_linear(val, enc);
-            }
-        }
-    }
-}
-
-fn downsample_kernel<T: PixelStorageTrait>(
-    dst: &mut [T],
-    dst_res: Point2i,
-    prev: &Image,
-    wrap: WrapMode2D,
-) {
-    let w = dst_res.x() as usize;
-    let channels = prev.n_channels();
-    let enc = prev.encoding();
-    let old_res = prev.resolution();
-
-    dst.par_chunks_mut(w * channels as usize)
-        .enumerate()
-        .for_each(|(y, row)| {
-            let src_y = y * 2;
-            for x in 0..w {
-                let src_x = x * 2;
-                for c in 0..channels {
-                    let mut sum = 0.0;
-                    let mut count = 0.0;
-
-                    for dy in 0..2 {
-                        for dx in 0..2 {
-                            let sx = src_x as i32 + dx;
-                            let sy = src_y as i32 + dy;
-                            if sx < old_res.x() && sy < old_res.y() {
-                                sum += prev.get_channel_with_wrap(Point2i::new(sx, sy), c, wrap);
-                                count += 1.0;
-                            }
-                        }
-                    }
-
-                    let avg = if count > 0.0 { sum / count } else { 0.0 };
-                    row[x * channels as usize + c as usize] = T::from_linear(avg, enc);
-                }
-            }
-        });
-}
-
 fn resample_weights(old_res: usize, new_res: usize) -> Vec<ResampleWeight> {
     let filter_radius = 2.0;
     let tau = 2.0;
@@ -416,7 +261,6 @@ fn compute_resize_tile(
 
     let mut x_buf = vec![0.0; n_channels * ny_in * nx_out];
 
-    // Resize X
     for y in 0..ny_in {
         for x in 0..nx_out {
             let x_global = out_extent.p_min.x() + x as i32;
@@ -438,7 +282,6 @@ fn compute_resize_tile(
 
     let mut out_buf = vec![0.0; n_channels * nx_out * ny_out];
 
-    // Resize Y
     for x in 0..nx_out {
         for y in 0..ny_out {
             let y_global = out_extent.p_min.y() + y as i32;

src/core/light.rs

@@ -6,9 +6,8 @@ use shared::core::camera::CameraTransform;
 use shared::core::light::Light;
 use shared::core::medium::Medium;
 use shared::core::shape::Shape;
-use shared::spectra::DenselySampledSpectrum;
 use shared::core::spectrum::Spectrum;
-use shared::spectra::RGBColorSpace;
+use shared::spectra::{DenselySampledSpectrum, RGBColorSpace};
 use shared::Transform;
 use std::sync::Arc;
 
@@ -18,9 +17,6 @@ pub fn lookup_spectrum(s: &Spectrum) -> Arc<DenselySampledSpectrum> {
     cache.lookup(dense_spectrum).into()
 }
 
-// Placeholders for non-area lights that never inspect these arguments.
-// TODO: refactor each light to only take what it actually needs,
-// then delete this bullshit
 fn dummy_shape() -> Shape {
     Shape::default()
 }
@@ -29,6 +25,8 @@ fn dummy_alpha() -> FloatTexture {
     FloatTexture::default()
 }
 
+/// Create a non-area light. Returns a Light value — the caller decides
+/// whether to wrap it in Arc (host ownership) or arena.alloc (GPU).
 pub fn create_light(
     name: &str,
     render_from_light: Transform,
@@ -36,79 +34,47 @@ pub fn create_light(
     parameters: &ParameterDictionary,
     loc: &FileLoc,
     camera_transform: CameraTransform,
-    arena: &mut Arena,
+    arena: &Arena,
 ) -> Result<Light> {
     let shape = dummy_shape();
     let alpha = dummy_alpha();
-
     match name {
         "point" => crate::lights::point::create(
-            render_from_light,
-            medium,
-            parameters,
-            loc,
-            &shape,
-            &alpha,
-            None,
-            arena,
+            render_from_light, medium, parameters, loc,
+            &shape, &alpha, None, arena,
         ),
         "spot" => crate::lights::spot::create(
-            render_from_light,
-            medium,
-            parameters,
-            loc,
-            &shape,
-            &alpha,
-            None,
-            arena,
+            render_from_light, medium, parameters, loc,
+            &shape, &alpha, None, arena,
         ),
         "distant" => crate::lights::distant::create(
-            render_from_light,
-            medium,
-            parameters,
-            loc,
-            &shape,
-            &alpha,
-            None,
-            arena,
+            render_from_light, medium, parameters, loc,
+            &shape, &alpha, None, arena,
         ),
         "goniometric" => crate::lights::goniometric::create(
-            render_from_light,
-            medium,
-            parameters,
-            loc,
-            &shape,
-            &alpha,
-            None,
-            arena,
+            render_from_light, medium, parameters, loc,
+            &shape, &alpha, None, arena,
         ),
         "projection" => crate::lights::projection::create(
-            render_from_light,
-            medium,
-            parameters,
-            loc,
-            &shape,
-            &alpha,
-            None,
-            arena,
+            render_from_light, medium, parameters, loc,
+            &shape, &alpha, None, arena,
         ),
         "infinite" => crate::lights::infinite::create(
-            render_from_light,
-            medium.into(),
-            camera_transform,
-            parameters,
-            None,
-            loc,
-            arena,
+            render_from_light, medium.into(), camera_transform,
+            parameters, None, loc, arena,
         ),
         "diffuse" => Err(anyhow!(
-            "{}: \"diffuse\" is an area light. Use create_area_light with a shape",
+            "{}: \"diffuse\" is an area light; use create_area_light with a shape",
             loc
         )),
         _ => Err(anyhow!("{}: unknown light type \"{}\"", loc, name)),
     }
 }
 
+/// Create a diffuse area light bound to a specific shape.
+/// Returns a Light value. The individual light constructor still uses
+/// the arena internally for sub-allocations (shape, image, spectra),
+/// but the Light itself is returned as a value for the caller to place.
 pub fn create_area_light(
     render_from_light: Transform,
     medium: Option<Medium>,
@@ -117,11 +83,10 @@ pub fn create_area_light(
     shape: &Shape,
     alpha_tex: &FloatTexture,
     colorspace: Option<&RGBColorSpace>,
-    arena: &mut Arena,
-) -> Result<Ptr<Light>> {
-    let light = crate::lights::diffuse::create(
+    arena: &Arena,
+) -> Result<Light> {
+    crate::lights::diffuse::create(
         render_from_light, medium, parameters, loc,
         shape, alpha_tex, colorspace, arena,
-    )?;
-    Ok(arena.alloc(light))
+    )
 }

src/core/material.rs

@@ -1,5 +1,5 @@
 use crate::Arena;
-use crate::core::image::Image;
+use crate::core::image::HostImage;
 use crate::utils::TextureParameterDictionary;
 use crate::utils::error::FileLoc;
 use anyhow::{Result, anyhow};
@@ -11,7 +11,7 @@ use std::sync::Arc;
 pub trait CreateMaterial: Sized {
     fn create(
         parameters: &TextureParameterDictionary,
-        normal_map: Option<Arc<Image>>,
+        normal_map: Option<Arc<HostImage>>,
         named_materials: &HashMap<String, Material>,
         loc: &FileLoc,
         arena: &Arena,
@@ -22,7 +22,7 @@ pub trait MaterialFactory {
     fn create(
         name: &str,
         params: &TextureParameterDictionary,
-        normal_map: Option<Arc<Image>>,
+        normal_map: Option<Arc<HostImage>>,
         named_materials: &HashMap<String, Material>,
         loc: FileLoc,
         arena: &Arena,
@@ -35,7 +35,7 @@ impl MaterialFactory for Material {
     fn create(
         name: &str,
         parameters: &TextureParameterDictionary,
-        normal_map: Option<Arc<Image>>,
+        normal_map: Option<Arc<HostImage>>,
         named_materials: &HashMap<String, Material>,
         loc: FileLoc,
         arena: &Arena,

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));
-        match result {
-            Some(t) => {
-                self.for_active_transforms(|cur| cur * &t);
-            }
-            None => {
-                eprintln!("Error: Could not invert transform at {}", loc);
-            }
-        }
+        let t = transform::look_at((ex, ey, ez), (lx, ly, lz), (ux, uy, uz))
+            .with_context(|| format!("at {}", loc))?;
+        self.for_active_transforms(|cur| cur * &t);
         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,20 +1,23 @@
 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;
-use crate::core::image::{io::ImageIO, Image};
+use crate::core::image::{HostImage, ImageIO};
 use crate::core::material::MaterialFactory;
 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, Upload};
+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;
@@ -25,9 +28,10 @@ use shared::core::medium::{Medium, MediumInterface};
 use shared::core::primitive::{AnimatedPrimitive, GeometricPrimitive, Primitive, SimplePrimitive};
 use shared::core::sampler::Sampler;
 use shared::core::shape::Shape;
-use shared::core::texture::SpectrumType;
+use shared::core::texture::{GPUFloatTexture, SpectrumType};
+use shared::lights::sampler::LightSampler;
 use shared::spectra::RGBColorSpace;
-use shared::utils::Ptr;
+use shared::{Ptr, Transform};
 use std::collections::HashMap;
 use std::sync::Arc;
 
@@ -130,8 +134,9 @@ impl BasicScene {
             *self.film_colorspace.lock() = Some(Arc::clone(cs));
         }
 
-        let filter = Filter::create(&filter.name, &filter.parameters, &filter.loc, &arena)
+        let filter = Filter::create(&filter.name, &filter.parameters, &filter.loc, arena)
             .map_err(|e| anyhow!("Failed to create filter: {}", e))?;
+
         let shutter_close = camera.base.parameters.get_one_float("shutterclose", 1.)?;
         let shutter_open = camera.base.parameters.get_one_float("shutteropen", 0.)?;
         let exposure_time = shutter_close - shutter_open;
@@ -144,7 +149,7 @@ impl BasicScene {
                 filter,
                 Some(camera.camera_transform.clone()),
                 &film.loc,
-                &arena,
+                arena,
             )
             .map_err(|e| anyhow!("Failed to create film: {}", e))?,
         );
@@ -154,38 +159,33 @@ impl BasicScene {
             job: None,
         };
 
-        let arena_sampler = Arc::clone(&arena);
-        let sampler_film = Arc::clone(&film_instance);
-        let sampler_job = run_async(move || {
-            let res = sampler_film.as_ref().base().full_resolution;
-            Sampler::create(
-                &sampler.name,
-                &sampler.parameters,
-                res,
-                &sampler.loc,
-                &arena_sampler,
-            )
-            .map_err(|e| anyhow!("Failed to create sampler: {}", e))
-        });
-        self.sampler_state.lock().job = Some(sampler_job);
+        let res = film_instance.as_ref().base().full_resolution;
+        let sampler_result =
+            Sampler::create(&sampler.name, &sampler.parameters, res, &sampler.loc, arena)
+                .map_err(|e| anyhow!("Failed to create sampler: {}", e))?;
+
+        *self.sampler_state.lock() = SingletonState {
+            result: Some(Arc::new(sampler_result)),
+            job: None,
+        };
+
+        let medium = self.get_medium(&camera.medium, &camera.base.loc);
+        let camera_result = Camera::create(
+            &camera.base.name,
+            &camera.base.parameters,
+            &camera.camera_transform,
+            medium,
+            Arc::clone(&film_instance),
+            &camera.base.loc,
+            arena,
+        )
+        .map_err(|e| anyhow!("Failed to create camera: {}", e))?;
+
+        *self.camera_state.lock() = SingletonState {
+            result: Some(Arc::new(camera_result)),
+            job: None,
+        };
 
-        let arena_camera = Arc::clone(&arena);
-        let camera_film = Arc::clone(&film_instance);
-        let scene_ptr = Arc::clone(self);
-        let camera_job = run_async(move || {
-            let medium = scene_ptr.get_medium(&camera.medium, &camera.base.loc);
-            Camera::create(
-                &camera.base.name,
-                &camera.base.parameters,
-                &camera.camera_transform,
-                medium,
-                camera_film,
-                &camera.base.loc,
-                &arena_camera,
-            )
-            .map_err(|e| anyhow!("Failed to create camera: {}", e))
-        });
-        self.camera_state.lock().job = Some(camera_job);
         Ok(())
     }
 
@@ -351,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();
@@ -412,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();
 
@@ -508,11 +508,7 @@ impl BasicScene {
         Ok((named_materials, materials))
     }
 
-    pub fn create_lights(
-        &self,
-        camera_transform: &CameraTransform,
-        arena: &mut Arena,
-    ) -> Vec<Light> {
+    pub fn create_lights(&self, camera_transform: &CameraTransform, arena: &Arena) -> Vec<Light> {
         let state = self.light_state.lock();
 
         state
@@ -583,7 +579,6 @@ impl BasicScene {
             let default_alpha = Arc::new(FloatTexture::default());
             let alpha_ref = alpha_tex.as_ref().unwrap_or(&default_alpha);
 
-            // Use the film colorspace as fallback for area light emission
             let film_cs = self.film_colorspace.lock();
             let colorspace_ref = al_entity
                 .parameters
@@ -628,8 +623,8 @@ impl BasicScene {
         textures: &NamedTextures,
         named_materials: &HashMap<String, Material>,
         materials: &[Material],
-        arena: &mut Arena,
-    ) -> (Vec<Primitive>, Vec<Arc<Light>>) {
+        arena: &Arena,
+    ) -> (Arc<Primitive>, Vec<Arc<Light>>) {
         let entities = self.shapes.lock();
         let animated = self.animated_shapes.lock();
         let light_state = self.light_state.lock();
@@ -669,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(
@@ -677,10 +701,10 @@ impl BasicScene {
         mtl: Material,
         alpha_tex: &Option<Arc<FloatTexture>>,
         mi: MediumInterface,
-        al_params: Option<&AreaLightEntity>,
+        al_params: Option<&SceneEntity>,
         render_from_light: Transform,
         film_cs: Option<&RGBColorSpace>,
-        arena: &mut Arena,
+        arena: &Arena,
         area_lights: &mut Vec<Arc<Light>>,
     ) -> Vec<(Ptr<Shape>, Ptr<Light>, Ptr<GPUFloatTexture>)> {
         shapes
@@ -691,7 +715,7 @@ impl BasicScene {
                         let cs = al_entity.parameters.color_space.as_deref().or(film_cs);
                         let default_alpha = Arc::new(FloatTexture::default());
                         let alpha_ref = alpha_tex.as_ref().unwrap_or(&default_alpha);
-                        crate::core::light::create_area_light(
+                        match crate::core::light::create_area_light(
                             render_from_light,
                             None,
                             &al_entity.parameters,
@@ -700,15 +724,21 @@ impl BasicScene {
                             alpha_ref,
                             cs,
                             arena,
-                        )
-                        .ok()
+                        ) {
+                            Ok(light) => {
+                                // Keep an Arc copy for the host-side light list
+                                area_lights.push(Arc::new(light));
+                                // Alloc into arena for the GPU primitive
+                                Some(arena.alloc(light))
+                            }
+                            Err(e) => {
+                                log::error!("Failed to create area light: {}", e);
+                                None
+                            }
+                        }
                     })
                     .unwrap_or(Ptr::null());
 
-                if !area_light_ptr.is_null() {
-                    area_lights.push(Arc::new(unsafe { &*area_light_ptr.as_raw() }.clone()));
-                }
-
                 let alpha_ptr = alpha_tex
                     .as_ref()
                     .map(|t| arena.upload(t.as_ref()))
@@ -727,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>>,
     ) {
@@ -810,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>>,
     ) {
@@ -934,23 +964,28 @@ impl BasicScene {
                     .unwrap_or(Ptr::null()),
             };
 
-            let shape_lights_opt = lookup
+            // Light is &Light from the HashMap — alloc into arena for the Ptr
+            let light_ptr = lookup
                 .shape_lights
                 .get(&shape_ctx.entity_index)
-                .and_then(|lights| lights.get(shape_ctx.shape_index));
+                .and_then(|lights| lights.get(shape_ctx.shape_index))
+                .map(|l| arena.alloc(*l))
+                .unwrap_or(Ptr::null());
 
-            let light_ptr = arena.alloc_opt(shape_lights_opt);
             let shape_ptr = shape_ctx.shape;
-            let prim = if area_light.is_null() && !mi.is_medium_transition() && alpha_tex.is_none()
-            {
+
+            let prim = if light_ptr.is_null() && !mi.is_medium_transition() && alpha_tex.is_none() {
                 Primitive::Simple(SimplePrimitive::new(shape_ptr, Ptr::from(&mtl)))
             } else {
                 Primitive::Geometric(GeometricPrimitive::new(
                     shape_ptr,
                     arena.alloc(mtl),
-                    area_light,
+                    light_ptr,
                     mi.clone(),
-                    arena.upload(alpha_tex),
+                    alpha_tex
+                        .as_ref()
+                        .map(|t| arena.upload(t.as_ref()))
+                        .unwrap_or(Ptr::null()),
                 ))
             };
 
@@ -971,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>> {
@@ -1027,7 +1076,7 @@ impl BasicScene {
         let filename_clone = filename.clone();
         let job = run_async(move || {
             let path = std::path::Path::new(&filename_clone);
-            let immeta = Image::read(path, Some(LINEAR)).expect(&format!(
+            let immeta = HostImage::read(path, Some(LINEAR)).expect(&format!(
                 "{}: normal map must contain R, G, B channels",
                 filename_clone
             ));
@@ -1051,7 +1100,7 @@ impl BasicScene {
         &self,
         state: &MaterialState,
         params: &ParameterDictionary,
-    ) -> Result<Option<Arc<Image>>> {
+    ) -> Result<Option<Arc<HostImage>>> {
         let filename = resolve_filename(&params.get_one_string("normalmap", "")?);
         if filename.is_empty() {
             return Ok(None);

src/core/scene/state.rs

@@ -1,5 +1,5 @@
 use super::{LightSceneEntity, SceneEntity, TextureSceneEntity};
-use crate::core::image::Image;
+use crate::core::image::HostImage;
 use crate::core::texture::{FloatTexture, SpectrumTexture};
 use crate::utils::parallel::AsyncJob;
 use anyhow::Result;
@@ -22,8 +22,8 @@ pub struct TextureState {
 pub struct MaterialState {
     pub named_materials: Vec<(String, SceneEntity)>,
     pub materials: Vec<SceneEntity>,
-    pub normal_map_jobs: HashMap<String, AsyncJob<Arc<Image>>>,
-    pub normal_maps: HashMap<String, Arc<Image>>,
+    pub normal_map_jobs: HashMap<String, AsyncJob<Arc<HostImage>>>,
+    pub normal_maps: HashMap<String, Arc<HostImage>>,
 }
 
 #[derive(Debug, Default)]

src/core/shape.rs

@@ -116,8 +116,8 @@ impl ShapeFactory for Shape {
                 global_store.push(host_arc.clone());
                 drop(global_store);
 
-                let n_tris = host_arc..n_triangles;
-                let mesh_ptr = Ptr::from(&host_arc);
+                let n_tris = host_arc.n_triangles;
+                let mesh_ptr = arena.alloc_arc(host_arc);
                 let shapes: Vec<Ptr<Shape>> = (0..n_tris)
                     .map(|i| {
                         let tri_shape = Shape::Triangle(TriangleShape {

src/core/texture.rs

@@ -1,5 +1,5 @@
 use crate::textures::*;
-use crate::utils::TextureParameterDictionary;
+use crate::utils::{MIPMap, MIPMapFilterOptions, TextureParameterDictionary};
 use crate::{Arena, FileLoc};
 use anyhow::{anyhow, Result};
 use enum_dispatch::enum_dispatch;
@@ -8,8 +8,8 @@ use shared::core::geometry::Vector3f;
 use shared::core::image::WrapMode;
 use shared::core::texture::SpectrumType;
 use shared::core::texture::{
-    CylindricalMapping, GPUFloatTexture, GPUSpectrumTexture, PlanarMapping, SphericalMapping,
-    TextureEvalContext, TextureMapping2D, UVMapping,
+    CylindricalMapping, PlanarMapping, SphericalMapping, TextureEvalContext, TextureMapping2D,
+    UVMapping,
 };
 use shared::spectra::{SampledSpectrum, SampledWavelengths};
 use shared::textures::*;
@@ -133,7 +133,10 @@ impl SpectrumTexture {
             invert: inner.base.invert,
             is_single_channel: inner.base.mipmap.is_single_channel(),
             color_space: arena.alloc(
-                inner.base.mipmap.color_space
+                inner
+                    .base
+                    .mipmap
+                    .color_space
                     .clone()
                     .unwrap_or_else(crate::spectra::default_colorspace),
             ),
@@ -262,4 +265,3 @@ pub struct TexInfo {
     pub wrap_mode: WrapMode,
     pub encoding: ColorEncoding,
 }
-

src/films/gbuffer.rs

@@ -1,9 +1,8 @@
 use super::*;
-use crate::core::film::{CreateFilmBase, PixelSensor};
-use crate::utils::containers::Array2D;
+use crate::core::film::{CreateFilmBase, CreatePixelSensor};
+use shared::core::film::PixelSensor;
 use anyhow::{Result, anyhow};
 use shared::core::film::{FilmBase, GBufferFilm};
-use shared::core::filter::FilterTrait;
 use shared::spectra::RGBColorSpace;
 use shared::utils::AnimatedTransform;
 use std::path::Path;
@@ -16,13 +15,14 @@ impl CreateFilm for GBufferFilm {
         filter: Filter,
         camera_transform: Option<CameraTransform>,
         loc: &FileLoc,
-        _arena: &Arena,
+        arena: &Arena,
     ) -> Result<Film> {
         let colorspace = params.color_space.as_ref().unwrap();
         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)?;
-        let film_base = FilmBase::create(params, filter, Some(&sensor.device()), loc)?;
+        let sensor = PixelSensor::create(params, colorspace.clone(), exposure_time, loc, arena)?;
+        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

@@ -1,11 +1,9 @@
 use super::*;
-use crate::core::film::{CreateFilmBase, PixelSensor};
-use crate::utils::containers::Array2D;
+use crate::core::film::{CreateFilmBase, CreatePixelSensor};
 use crate::Arena;
 use anyhow::Result;
 use shared::core::camera::CameraTransform;
-use shared::core::film::{Film, FilmBase, RGBFilm, RGBPixel};
-use shared::core::filter::FilterTrait;
+use shared::core::film::{Film, FilmBase, RGBFilm, RGBPixel, PixelSensor};
 use shared::spectra::RGBColorSpace;
 
 impl CreateFilm for RGBFilm {
@@ -15,13 +13,14 @@ impl CreateFilm for RGBFilm {
         filter: Filter,
         _camera_transform: Option<CameraTransform>,
         loc: &FileLoc,
-        _arena: &Arena,
+        arena: &Arena,
     ) -> Result<Film> {
         let colorspace = params.color_space.as_ref().cloned().unwrap_or_else(crate::spectra::default_colorspace_arc);
         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)?;
-        let film_base = FilmBase::create(params, filter, Some(&sensor.device()), loc)?;
+        let sensor = PixelSensor::create(params, colorspace.clone(), exposure_time, loc, arena)?;
+        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

@@ -1,16 +1,14 @@
 use super::*;
-use crate::core::film::{CreateFilmBase, PixelSensor};
+use crate::core::film::{CreateFilmBase, CreatePixelSensor};
 use crate::{Arena, FileLoc, ParameterDictionary};
-use anyhow::{Result, anyhow};
-use shared::Float;
+use anyhow::{anyhow, Result};
 use shared::core::camera::CameraTransform;
-use shared::core::film::{FilmBase, SpectralFilm};
-use shared::core::filter::FilterTrait;
+use shared::core::film::{FilmBase, PixelSensor, SpectralFilm};
 use shared::spectra::{LAMBDA_MAX, LAMBDA_MIN};
 use shared::utils::math::SquareMatrix;
+use shared::Float;
 use std::path::Path;
 
-
 impl CreateFilm for SpectralFilm {
     fn create(
         params: &ParameterDictionary,
@@ -18,14 +16,15 @@ impl CreateFilm for SpectralFilm {
         filter: Filter,
         _camera_transform: Option<CameraTransform>,
         loc: &FileLoc,
-        _arena: &Arena,
+        arena: &Arena,
     ) -> Result<Film> {
         // Missing default illuminant, use srgb
         let colorspace = params.color_space.as_ref().unwrap();
         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)?;
-        let film_base = FilmBase::create(params, filter, Some(&sensor.device()), loc)?;
+        let sensor = PixelSensor::create(params, colorspace.clone(), exposure_time, loc, arena)?;
+        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

@@ -2,7 +2,7 @@ use super::base::IntegratorBase;
 use super::RayIntegratorTrait;
 use crate::core::camera::InitMetadata;
 use crate::core::film::FilmTrait;
-use crate::core::image::{Image, ImageIO, ImageMetadata};
+use crate::core::image::{HostImage, ImageIO, ImageMetadata};
 use crate::globals::get_options;
 use crate::spectra::get_spectra_context;
 use crate::Arena;
@@ -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);
@@ -108,12 +114,12 @@ pub fn render<T>(
     let mut wave_end = 1;
     let mut next_wave_size = 1;
 
-    let mut reference_image: Option<Image> = None;
+    let mut reference_image: Option<HostImage> = None;
     let mut mse_out_file: Option<std::fs::File> = None;
 
     if let Some(ref_path) = &options.mse_reference_image {
         let image_and_metadata =
-            Image::read(Path::new(&ref_path), None).expect("Could not load image");
+            HostImage::read(Path::new(&ref_path), None).expect("Could not load image");
         let image = image_and_metadata.image;
         let metadata = image_and_metadata.metadata;
         let resolution = image.resolution();
@@ -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/lib.rs

@@ -2,7 +2,6 @@
 #[allow(dead_code)]
 pub mod core;
 pub mod films;
-pub mod filters;
 pub mod globals;
 pub mod integrators;
 pub mod lights;

src/lights/diffuse.rs

@@ -1,6 +1,4 @@
-use std::path::Path;
-
-use crate::core::image::{Image, ImageIO};
+use crate::core::image::{HostImage, ImageIO};
 use crate::core::light::lookup_spectrum;
 use crate::core::spectrum::spectrum_to_photometric;
 use crate::core::texture::FloatTexture;
@@ -13,11 +11,12 @@ use shared::core::light::{Light, LightBase, LightType};
 use shared::core::medium::{Medium, MediumInterface};
 use shared::core::shape::{Shape, ShapeTrait};
 use shared::core::spectrum::Spectrum;
-use shared::core::texture::{SpectrumType, TextureEvalContext};
+use shared::core::texture::{GPUFloatTexture, SpectrumType, TextureEvalContext};
 use shared::lights::DiffuseAreaLight;
 use shared::spectra::RGBColorSpace;
 use shared::utils::Transform;
 use shared::{Float, PI};
+use std::path::Path;
 
 pub fn create(
     render_from_light: Transform,
@@ -37,13 +36,13 @@ pub fn create(
     let two_sided = params.get_one_bool("twosided", false)?;
 
     let filename = resolve_filename(&params.get_one_string("filename", "")?);
-    let (image, image_color_space): (Option<Image>, Option<RGBColorSpace>) =
+    let (image, image_color_space): (Option<HostImage>, Option<RGBColorSpace>) =
         if !filename.is_empty() {
             if l.is_some() {
                 return Err(anyhow!("{}: both \"L\" and \"filename\" specified", loc));
             }
 
-            let im = Image::read(Path::new(&filename), None)?;
+            let im = HostImage::read(Path::new(&filename), None)?;
 
             if im.image.has_any_infinite_pixels() {
                 return Err(anyhow!("{}: image has infinite pixel values", loc));
@@ -100,12 +99,13 @@ pub fn create(
         scale *= phi_v / k_e;
     }
 
-    // Upload alpha texture to GPU and check for constant-zero
+    // Upload alpha texture to GPU and check for null texture
     let alpha_ptr = arena.upload(alpha);
-    let light_type = if alpha_ptr.is_constant_zero() {
-        LightType::DeltaPosition
-    } else {
-        LightType::Area
+    let light_type = match unsafe { alpha_ptr.as_ref() } {
+        GPUFloatTexture::Constant(t) if t.evaluate(&TextureEvalContext::default()) == 0.0 => {
+            LightType::DeltaPosition
+        }
+        _ => LightType::Area,
     };
 
     let mi = match medium {
@@ -146,7 +146,7 @@ pub fn create(
         area: shape.area(),
         shape: arena.alloc(*shape),
         alpha: alpha_ptr,
-        image: arena.alloc(image),
+        image: arena.upload(image),
         colorspace: arena.alloc_opt(image_color_space),
         lemit: arena.alloc((*lookup_spectrum(l_for_scale)).clone()),
         two_sided,

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/goniometric.rs

@@ -1,13 +1,10 @@
-use std::path::Path;
-
-use crate::core::image::{Image, ImageIO};
+use crate::core::image::{HostImage, ImageIO};
 use crate::core::light::lookup_spectrum;
 use crate::core::spectrum::spectrum_to_photometric;
 use crate::core::texture::FloatTexture;
-use crate::utils::sampling::PiecewiseConstant2D;
 use crate::utils::resolve_filename;
 use crate::{Arena, FileLoc, ParameterDictionary};
-use anyhow::{Result, anyhow};
+use anyhow::{anyhow, Result};
 use shared::core::geometry::Point2i;
 use shared::core::light::{Light, LightBase, LightType};
 use shared::core::medium::{Medium, MediumInterface};
@@ -16,8 +13,10 @@ use shared::core::spectrum::Spectrum;
 use shared::core::texture::SpectrumType;
 use shared::lights::GoniometricLight;
 use shared::spectra::RGBColorSpace;
+use shared::utils::sampling::PiecewiseConstant2D;
 use shared::utils::{Ptr, Transform};
 use shared::{Float, PI};
+use std::path::Path;
 
 pub fn create(
     render_from_light: Transform,
@@ -25,11 +24,10 @@ pub fn create(
     params: &ParameterDictionary,
     loc: &FileLoc,
     _shape: &Shape,
-    _alpha_text: &FloatTexture,
+    _alpha_tex: &FloatTexture,
     colorspace: Option<&RGBColorSpace>,
     arena: &Arena,
 ) -> Result<Light> {
-
     let default_cs = crate::spectra::default_colorspace();
     let cs = colorspace.unwrap_or(&default_cs);
     let i = params
@@ -39,14 +37,13 @@ pub fn create(
             SpectrumType::Illuminant,
         )
         .expect("Could not retrieve spectrum");
+
     let mut scale = params.get_one_float("scale", 1.)?;
     let filename = resolve_filename(&params.get_one_string("filename", "")?);
-    let image: Ptr<Image> = if filename.is_empty() {
-        Ptr::null()
-    } else {
-        let im = Image::read(Path::new(&filename), None)
-            .map_err(|e| anyhow!("could not load image '{}': {}", filename, e))?;
 
+    let host_image = if !filename.is_empty() {
+        let im = HostImage::read(Path::new(&filename), None)
+            .map_err(|e| anyhow!("could not load image '{}': {}", filename, e))?;
         let loaded = im.image;
         let res = loaded.resolution();
 
@@ -56,7 +53,6 @@ pub fn create(
                 filename
             ));
         }
-
         if res.x() != res.y() {
             return Err(anyhow!(
                 "image resolution ({}, {}) is non-square; unlikely to be an equal-area map",
@@ -65,15 +61,19 @@ pub fn create(
             ));
         }
 
-        Ptr::from(&convert_to_luminance_image(&loaded, &filename, loc)?)
+        Some(convert_to_luminance_image(&loaded, &filename, loc)?)
+    } else {
+        None
     };
 
     scale /= spectrum_to_photometric(i);
-    let phi_v = params.get_one_float("power", -1.0)?;
 
+    let phi_v = params.get_one_float("power", -1.0)?;
     if phi_v > 0.0 {
-        let k_e = compute_emissive_power(&image);
-        scale *= phi_v / k_e;
+        if let Some(ref img) = host_image {
+            let k_e = compute_emissive_power(img);
+            scale *= phi_v / k_e;
+        }
     }
 
     let swap_yz: [Float; 16] = [
@@ -81,7 +81,7 @@ pub fn create(
     ];
     let t =
         Transform::from_flat(&swap_yz).expect("Could not create transform for GoniometricLight");
-    let _final_render_from_light = render_from_light * t;
+    let final_render_from_light = render_from_light * t;
 
     let mi = match medium {
         Some(m) => {
@@ -93,31 +93,35 @@ pub fn create(
         }
         None => MediumInterface::default(),
     };
-    let base = LightBase::new(LightType::DeltaPosition, render_from_light, mi);
 
+    let base = LightBase::new(LightType::DeltaPosition, final_render_from_light, mi);
     let iemit = lookup_spectrum(&i);
 
-    let image_ptr = if !image.is_null() {
-        let distrib = PiecewiseConstant2D::from_image(&image);
-        let distrib_ptr = arena.alloc(distrib);
-        let img_ptr = arena.alloc(image);
-        (img_ptr, distrib_ptr)
-    } else {
-        (Ptr::null(), Ptr::null())
+    // Build distribution from host image, then upload both to arena
+    let (image_ptr, distrib_ptr) = match host_image {
+        Some(img) => {
+            let distrib = PiecewiseConstant2D::from_image(&img.inner);
+            (arena.alloc(img.inner), arena.alloc(distrib))
+        }
+        None => (Ptr::null(), Ptr::null()),
     };
 
     let specific = GoniometricLight {
         base,
-        iemit: arena.alloc(*iemit),
+        iemit: arena.alloc((*iemit).clone()),
         scale,
-        image: image_ptr.0,
-        distrib: image_ptr.1,
+        image: image_ptr,
+        distrib: distrib_ptr,
     };
 
     Ok(Light::Goniometric(specific))
 }
 
-fn convert_to_luminance_image(image: &Image, filename: &str, loc: &FileLoc) -> Result<Image> {
+fn convert_to_luminance_image(
+    image: &HostImage,
+    filename: &str,
+    loc: &FileLoc,
+) -> Result<HostImage> {
     let res = image.resolution();
     let rgb_desc = image.get_channel_desc(&["R", "G", "B"]);
     let y_desc = image.get_channel_desc(&["Y"]);
@@ -142,7 +146,7 @@ fn convert_to_luminance_image(image: &Image, filename: &str, loc: &FileLoc) -> R
                 }
             }
 
-            Ok(Image::from_f32(y_pixels, res, &["Y"].to_vec()))
+            Ok(HostImage::from_f32(&y_pixels, res, &["Y"].to_vec()))
         }
 
         (Err(_), Ok(_)) => {
@@ -158,7 +162,7 @@ fn convert_to_luminance_image(image: &Image, filename: &str, loc: &FileLoc) -> R
     }
 }
 
-fn compute_emissive_power(image: &Image) -> Float {
+fn compute_emissive_power(image: &HostImage) -> Float {
     let res = image.resolution();
     let mut sum_y = 0.0;
 

src/lights/infinite.rs

@@ -3,19 +3,18 @@ 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::utils::sampling::{PiecewiseConstant2D, WindowedPiecewiseConstant2D};
-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;
 use shared::core::geometry::{cos_theta, Bounds2f, Frame, Point2f, Point2i, Point3f, VectorLike};
 use shared::core::image::WrapMode;
 use shared::core::light::{Light, LightBase, LightType};
-use shared::core::medium::{Medium, MediumInterface};
+use shared::core::medium::Medium;
 use shared::core::spectrum::Spectrum;
 use shared::core::texture::SpectrumType;
 use shared::lights::{ImageInfiniteLight, PortalInfiniteLight, UniformInfiniteLight};
-use shared::spectra::{DenselySampledSpectrum, RGBColorSpace};
+use shared::spectra::RGBColorSpace;
 use shared::utils::math::{equal_area_sphere_to_square, equal_area_square_to_sphere};
 use shared::utils::sampling::{PiecewiseConstant2D, WindowedPiecewiseConstant2D};
 use shared::{Float, Ptr, Transform, PI};
@@ -63,7 +62,7 @@ pub fn create(
         }
 
         let lemit = lookup_spectrum(&spectrum);
-        let light = UniformInfiniteLight::new(render_from_light, scale, arena.alloc(*lemit));
+        let light = UniformInfiniteLight::new(render_from_light, scale, arena.alloc_arc(lemit));
         return Ok(Light::InfiniteUniform(light));
     }
 
@@ -96,7 +95,7 @@ pub fn create(
 fn create_image_light(
     render_from_light: Transform,
     scale: Float,
-    image: Image,
+    image: HostImage,
     image_cs: RGBColorSpace,
     arena: &Arena,
 ) -> Result<Light> {
@@ -110,7 +109,6 @@ fn create_image_light(
     let (n_u, n_v) = (res.x() as usize, res.y() as usize);
 
     // Extract luminance data
-    let image_ptr = image.upload(arena);
     let value = &image;
     let mut data: Vec<Float> = (0..n_v)
         .flat_map(|v| {
@@ -139,7 +137,7 @@ fn create_image_light(
     let light = ImageInfiniteLight::new(
         render_from_light,
         scale,
-        image_ptr,
+        arena.upload(image),
         arena.alloc(image_cs),
         arena.alloc(distrib),
         arena.alloc(compensated_distrib),
@@ -151,7 +149,7 @@ fn create_image_light(
 fn create_portal_light(
     render_from_light: Transform,
     scale: Float,
-    image: Image,
+    image: HostImage,
     image_cs: RGBColorSpace,
     portal_points: &[Point3f],
     camera_transform: CameraTransform,
@@ -198,7 +196,7 @@ fn create_portal_light(
     let light = PortalInfiniteLight::new(
         render_from_light,
         scale,
-        arena.alloc(remapped),
+        arena.upload(remapped),
         arena.alloc(image_cs),
         portal,
         portal_frame,
@@ -230,10 +228,10 @@ fn validate_and_build_portal_frame(portal: &[Point3f; 4], loc: &FileLoc) -> Resu
 }
 
 fn remap_image_through_portal(
-    image: &Image,
+    image: &HostImage,
     render_from_light: &Transform,
     portal_frame: &Frame,
-) -> Image {
+) -> HostImage {
     let res = image.resolution();
     let (width, height) = (res.x() as usize, res.y() as usize);
 
@@ -263,7 +261,7 @@ fn remap_image_through_portal(
             }
         });
 
-    Image::from_f32(pixels, res, &["R", "G", "B"])
+    HostImage::from_f32(&pixels, res, &["R", "G", "B"])
 }
 
 fn load_image(
@@ -271,16 +269,16 @@ fn load_image(
     l: &[Spectrum],
     colorspace: &RGBColorSpace,
     loc: &FileLoc,
-) -> Result<(Image, RGBColorSpace)> {
+) -> Result<(HostImage, RGBColorSpace)> {
     if filename.is_empty() {
         let stdspec = get_spectra_context();
         let rgb = l[0].to_rgb(colorspace, &stdspec);
         let image =
-            Image::new_constant(Point2i::new(1, 1), &["R", "G", "B"], &[rgb.r, rgb.g, rgb.b]);
+            HostImage::new_constant(Point2i::new(1, 1), &["R", "G", "B"], &[rgb.r, rgb.g, rgb.b]);
         return Ok((image, colorspace.clone()));
     }
 
-    let im = Image::read(Path::new(filename), None)
+    let im = HostImage::read(Path::new(filename), None)
         .map_err(|e| anyhow!("failed to load '{}': {}", filename, e))?;
 
     if im.image.has_any_infinite_pixels() || im.image.has_any_nan_pixels() {
@@ -296,7 +294,7 @@ fn load_image(
     Ok((im.image.select_channels(&desc), cs))
 }
 
-fn compute_hemisphere_illuminance(image: &Image, cs: &RGBColorSpace) -> Float {
+fn compute_hemisphere_illuminance(image: &HostImage, cs: &RGBColorSpace) -> Float {
     let lum = cs.luminance_vector();
     let res = image.resolution();
 

src/lights/projection.rs

@@ -1,13 +1,11 @@
-use crate::core::image::{Image, ImageIO};
+use crate::core::image::{HostImage, ImageIO};
 use crate::core::spectrum::spectrum_to_photometric;
 use crate::core::texture::FloatTexture;
-use crate::utils::sampling::PiecewiseConstant2D;
 use crate::utils::resolve_filename;
-use crate::{Arena, FileLoc, ParameterDictionary};
-use anyhow::{Result, anyhow};
-use shared::Float;
+use crate::{Arena, ArenaUpload, FileLoc, ParameterDictionary};
+use anyhow::{anyhow, Result};
 use shared::core::geometry::{
-    Bounds2f, Point2f, Point2i, Point3f, Vector3f, VectorLike, cos_theta,
+    cos_theta, Bounds2f, Point2f, Point2i, Point3f, Vector3f, VectorLike,
 };
 use shared::core::light::{Light, LightBase, LightType};
 use shared::core::medium::{Medium, MediumInterface};
@@ -15,8 +13,10 @@ use shared::core::shape::Shape;
 use shared::core::spectrum::Spectrum;
 use shared::lights::ProjectionLight;
 use shared::spectra::RGBColorSpace;
-use shared::utils::Transform;
 use shared::utils::math::{radians, square};
+use shared::utils::sampling::PiecewiseConstant2D;
+use shared::utils::Transform;
+use shared::Float;
 use std::path::Path;
 
 pub fn create(
@@ -41,7 +41,7 @@ pub fn create(
         ));
     }
 
-    let im = Image::read(Path::new(&filename), None)
+    let im = HostImage::read(Path::new(&filename), None)
         .map_err(|e| anyhow!("{}: could not load image '{}': {}", loc, filename, e))?;
 
     if im.image.has_any_infinite_pixels() {
@@ -125,9 +125,9 @@ pub fn create(
 
     let specific = ProjectionLight {
         base,
-        image: image.upload(arena),
-        image_color_space: colorspace.upload(arena),
-        distrib: distrib.upload(arena),
+        image: arena.upload(image),
+        image_color_space: arena.alloc(colorspace),
+        distrib: arena.alloc(distrib),
         screen_bounds,
         screen_from_light,
         light_from_screen,
@@ -150,7 +150,7 @@ fn compute_screen_bounds(aspect: Float) -> Bounds2f {
     }
 }
 
-fn compute_emissive_power(image: &Image, colorspace: &RGBColorSpace, fov: Float) -> Float {
+fn compute_emissive_power(image: &HostImage, colorspace: &RGBColorSpace, fov: Float) -> Float {
     let res = image.resolution();
     let aspect = res.x() as Float / res.y() as Float;
     let screen_bounds = compute_screen_bounds(aspect);

src/lights/sampler.rs

@@ -1,60 +1,83 @@
 use crate::Arena;
 use shared::core::light::{Light, LightTrait};
-use shared::lights::sampler::PowerLightSampler;
+use shared::lights::sampler::{
+    BVHLightSampler, LightSampler, PowerLightSampler, UniformLightSampler,
+};
 use shared::utils::sampling::AliasTable;
 use shared::spectra::{SampledSpectrum, SampledWavelengths};
-use std::collections::HashMap;
+use shared::utils::Ptr;
+use shared::Float;
 use std::sync::Arc;
 
-pub struct PowerSamplerHost {
-    pub lights: Vec<Arc<Light>>,
-    pub light_to_index: HashMap<usize, usize>,
-    pub alias_table: AliasTable,
-}
-
-impl PowerSamplerHost {
-    pub fn new(lights: &[Arc<Light>]) -> Self {
-        if lights.is_empty() {
-            return Self {
-                lights: Vec::new(),
-                light_to_index: HashMap::new(),
-                alias_table: AliasTable::new(&[]),
-            };
+pub fn create_light_sampler(
+    name: &str,
+    lights: &[Arc<Light>],
+    arena: &Arena,
+) -> LightSampler {
+    let device_lights = lights_to_slice(lights, arena);
+    match name {
+        "uniform" => LightSampler::Uniform(create_uniform(device_lights, lights.len())),
+        "power" => LightSampler::Power(create_power(lights, device_lights, arena)),
+        "bvh" => {
+            log::warn!("BVH light sampler not yet implemented, falling back to power");
+            LightSampler::Power(create_power(lights, device_lights, arena))
         }
-
-        let mut lights_vec = Vec::with_capacity(lights.len());
-        let mut light_to_index = HashMap::with_capacity(lights.len());
-        let mut light_power = Vec::with_capacity(lights.len());
-
-        let lambda = SampledWavelengths::sample_visible(0.5);
-
-        for (i, light) in lights.iter().enumerate() {
-            lights_vec.push(light.clone());
-
-            let ptr = Arc::as_ptr(light) as usize;
-            light_to_index.insert(ptr, i);
-
-            let phi = SampledSpectrum::safe_div(&light.phi(lambda), &lambda.pdf());
-            light_power.push(phi.average());
-        }
-
-        let alias_table = AliasTable::new(&light_power);
-
-        Self {
-            lights: lights_vec,
-            light_to_index,
-            alias_table,
+        _ => {
+            log::error!("Unknown light sampler \"{}\", using power", name);
+            LightSampler::Power(create_power(lights, device_lights, arena))
         }
     }
-    // pub fn to_device(&self, arena: &Arena) -> PowerLightSampler {
-    //     let device_lights: Vec<Light> = self.lights.iter().map(|l| (**l).clone()).collect();
-    //     let (lights_ptr, _) = arena.alloc_slice(&device_lights);
-    //     let alias_device = self.alias_table.to_device(arena);
-    //
-    //     PowerLightSampler {
-    //         lights: lights_ptr,
-    //         lights_len: self.lights.len() as u32,
-    //         alias_table: alias_device,
-    //     }
-    // }
+}
+
+fn lights_to_slice(lights: &[Arc<Light>], arena: &Arena) -> (Ptr<Light>, u32) {
+    if lights.is_empty() {
+        return (Ptr::null(), 0);
+    }
+    let vals: Vec<Light> = lights.iter().map(|l| **l).collect();
+    let (ptr, _) = arena.alloc_slice(&vals);
+    (ptr, lights.len() as u32)
+}
+
+fn create_uniform(
+    (lights, lights_len): (Ptr<Light>, u32),
+    _count: usize,
+) -> UniformLightSampler {
+    UniformLightSampler::new(lights, lights_len)
+}
+
+fn create_power(
+    host_lights: &[Arc<Light>],
+    (lights, lights_len): (Ptr<Light>, u32),
+    arena: &Arena,
+) -> PowerLightSampler {
+    if host_lights.is_empty() {
+        return PowerLightSampler {
+            lights: Ptr::null(),
+            lights_len: 0,
+            alias_table: Ptr::null(),
+        };
+    }
+
+    let lambda = SampledWavelengths::sample_visible(0.5);
+    let mut light_power: Vec<Float> = host_lights
+        .iter()
+        .map(|l| {
+            let phi = SampledSpectrum::safe_div(&l.phi(lambda), &lambda.pdf());
+            phi.average()
+        })
+        .collect();
+
+    // If all lights have zero power, treat as uniform
+    if light_power.iter().sum::<Float>() == 0.0 {
+        light_power.fill(1.0);
+    }
+
+    let alias_table = AliasTable::new(&light_power);
+    let alias_ptr = arena.alloc(alias_table);
+
+    PowerLightSampler {
+        lights,
+        lights_len,
+        alias_table: alias_ptr,
+    }
 }

src/lights/spot.rs

@@ -1,7 +1,7 @@
 use crate::core::light::lookup_spectrum;
 use crate::core::spectrum::spectrum_to_photometric;
 use crate::core::texture::FloatTexture;
-use crate::utils::{Arena, FileLoc, ParameterDictionary};
+use crate::{Arena, FileLoc, ParameterDictionary};
 use anyhow::Result;
 use shared::core::geometry::{Frame, Point3f, VectorLike};
 use shared::core::light::{Light, LightBase, LightType};
@@ -12,7 +12,6 @@ use shared::core::texture::SpectrumType;
 use shared::lights::SpotLight;
 use shared::spectra::RGBColorSpace;
 use shared::utils::math::radians;
-use shared::utils::{Ptr, Transform};
 use shared::{Float, Ptr, Transform, PI};
 
 trait CreateSpotLight {
@@ -23,6 +22,7 @@ trait CreateSpotLight {
         scale: Float,
         cos_falloff_start: Float,
         total_width: Float,
+        arena: &Arena
     ) -> Self;
 }
 
@@ -34,6 +34,7 @@ impl CreateSpotLight for SpotLight {
         scale: Float,
         cos_falloff_start: Float,
         total_width: Float,
+        arena: &Arena
     ) -> Self {
         let base = LightBase::new(
             LightType::DeltaPosition,
@@ -42,7 +43,7 @@ impl CreateSpotLight for SpotLight {
         );
 
         let i = lookup_spectrum(&le);
-        let iemit = arena.alloc(i);
+        let iemit = arena.alloc_arc(i);
         Self {
             base,
             iemit,
@@ -101,7 +102,7 @@ pub fn create(
         None => MediumInterface::default(),
     };
 
-    let specific = SpotLight::new(final_render, mi, i, scale, coneangle, coneangle - conedelta);
+    let specific = SpotLight::new(final_render, mi, i, scale, coneangle, coneangle - conedelta, arena);
     arena.alloc(specific);
     Ok(Light::Spot(specific))
 }

src/materials/coated.rs

@@ -1,10 +1,10 @@
-use crate::core::image::Image;
+use crate::core::image::HostImage;
 use crate::core::material::CreateMaterial;
 use crate::core::texture::SpectrumTexture;
 use crate::globals::get_options;
 use crate::spectra::data::get_named_spectrum;
 use crate::utils::TextureParameterDictionary;
-use crate::{Arena, FileLoc, Upload, ArenaUpload};
+use crate::{Arena, FileLoc, ArenaUpload};
 use anyhow::{bail, Result};
 use shared::core::material::Material;
 use shared::core::spectrum::Spectrum;
@@ -18,7 +18,7 @@ use std::sync::Arc;
 impl CreateMaterial for CoatedDiffuseMaterial {
     fn create(
         parameters: &TextureParameterDictionary,
-        normal_map: Option<Arc<Image>>,
+        normal_map: Option<Arc<HostImage>>,
         _named_materials: &HashMap<String, Material>,
         _loc: &FileLoc,
         arena: &Arena,
@@ -64,7 +64,7 @@ impl CreateMaterial for CoatedDiffuseMaterial {
             arena.upload(g),
             arena.upload(displacement),
             arena.alloc(eta),
-            arena.alloc(normal_map),
+            arena.upload(normal_map),
             remap_roughness,
             max_depth as u32,
             n_samples as u32,
@@ -78,7 +78,7 @@ impl CreateMaterial for CoatedDiffuseMaterial {
 impl CreateMaterial for CoatedConductorMaterial {
     fn create(
         parameters: &TextureParameterDictionary,
-        normal_map: Option<Arc<Image>>,
+        normal_map: Option<Arc<HostImage>>,
         _named_materials: &HashMap<String, Material>,
         loc: &FileLoc,
         arena: &Arena,
@@ -154,7 +154,7 @@ impl CreateMaterial for CoatedConductorMaterial {
         let remap_roughness = parameters.get_one_bool("remaproughness", true)?;
 
         let material = Self::new(
-            arena.upload(displacement)
+            arena.upload(displacement),
             arena.upload(interface_u_roughness),
             arena.upload(interface_v_roughness),
             arena.upload(thickness),
@@ -165,7 +165,7 @@ impl CreateMaterial for CoatedConductorMaterial {
             arena.upload(conductor_eta),
             arena.upload(k),
             arena.upload(reflectance),
-            arena.alloc(normal_map),
+            arena.upload(normal_map),
             arena.alloc(interface_eta),
             max_depth as u32,
             n_samples as u32,

src/materials/complex.rs

@@ -1,9 +1,9 @@
-use crate::core::image::Image;
+use crate::core::image::HostImage;
 use crate::core::material::CreateMaterial;
 use crate::core::texture::SpectrumTexture;
 use crate::spectra::get_colorspace_device;
 use crate::utils::TextureParameterDictionary;
-use crate::{Arena, FileLoc, Upload, ArenaUpload};
+use crate::{Arena, ArenaUpload, FileLoc};
 use anyhow::Result;
 use shared::bxdfs::HairBxDF;
 use shared::core::material::Material;
@@ -17,7 +17,7 @@ use std::sync::Arc;
 impl CreateMaterial for HairMaterial {
     fn create(
         parameters: &TextureParameterDictionary,
-        _normal_map: Option<Arc<Image>>,
+        _normal_map: Option<Arc<HostImage>>,
         _named_materials: &HashMap<String, Material>,
         _loc: &FileLoc,
         arena: &Arena,
@@ -60,11 +60,10 @@ impl CreateMaterial for HairMaterial {
     }
 }
 
-
 impl CreateMaterial for SubsurfaceMaterial {
     fn create(
         _parameters: &TextureParameterDictionary,
-        _normal_map: Option<Arc<Image>>,
+        _normal_map: Option<Arc<HostImage>>,
         _named_materials: &HashMap<String, Material>,
         _loc: &FileLoc,
         _arena: &Arena,
@@ -76,7 +75,7 @@ impl CreateMaterial for SubsurfaceMaterial {
 impl CreateMaterial for MeasuredMaterial {
     fn create(
         _parameters: &TextureParameterDictionary,
-        _normal_map: Option<Arc<Image>>,
+        _normal_map: Option<Arc<HostImage>>,
         _named_materials: &HashMap<String, Material>,
         _loc: &FileLoc,
         _arena: &Arena,

src/materials/conductor.rs

@@ -1,4 +1,4 @@
-use crate::core::image::Image;
+use crate::core::image::HostImage;
 use crate::core::material::CreateMaterial;
 // use crate::core::scattering::TrowbridgeReitzDistribution;
 use crate::utils::TextureParameterDictionary;
@@ -10,7 +10,7 @@ use std::sync::Arc;
 impl CreateMaterial for ConductorMaterial {
     fn create(
         _parameters: &TextureParameterDictionary,
-        _normal_map: Option<Arc<Image>>,
+        _normal_map: Option<Arc<HostImage>>,
         _named_materials: &std::collections::HashMap<String, shared::core::material::Material>,
         _loc: &crate::utils::FileLoc,
         _arena: &crate::Arena,

src/materials/dielectric.rs

@@ -1,5 +1,5 @@
 use crate::Arena;
-use crate::core::image::Image;
+use crate::core::image::HostImage;
 use crate::core::material::CreateMaterial;
 use crate::utils::{FileLoc, TextureParameterDictionary};
 use anyhow::Result;
@@ -11,7 +11,7 @@ use std::sync::Arc;
 impl CreateMaterial for DielectricMaterial {
     fn create(
         _parameters: &TextureParameterDictionary,
-        _normal_map: Option<Arc<Image>>,
+        _normal_map: Option<Arc<HostImage>>,
         _named_materials: &HashMap<String, Material>,
         _loc: &FileLoc,
         _arena: &Arena,
@@ -23,7 +23,7 @@ impl CreateMaterial for DielectricMaterial {
 impl CreateMaterial for ThinDielectricMaterial {
     fn create(
         _parameters: &TextureParameterDictionary,
-        _normal_map: Option<Arc<Image>>,
+        _normal_map: Option<Arc<HostImage>>,
         _named_materials: &HashMap<String, Material>,
         _loc: &FileLoc,
         _arena: &Arena,

src/materials/diffuse.rs

@@ -1,29 +1,49 @@
-use crate::Arena;
-use crate::core::image::Image;
+use crate::core::image::HostImage;
 use crate::core::material::CreateMaterial;
+use crate::core::texture::SpectrumTexture;
+use crate::utils::upload::ArenaUpload;
 use crate::utils::{FileLoc, TextureParameterDictionary};
+use crate::Arena;
 use anyhow::Result;
 use shared::core::material::Material;
+use shared::core::spectrum::Spectrum;
+use shared::core::texture::SpectrumType;
 use shared::materials::{DiffuseMaterial, DiffuseTransmissionMaterial};
+use shared::spectra::ConstantSpectrum;
+use shared::textures::SpectrumConstantTexture;
 use std::collections::HashMap;
 use std::sync::Arc;
 
 impl CreateMaterial for DiffuseMaterial {
     fn create(
-        _parameters: &TextureParameterDictionary,
-        _normal_map: Option<Arc<Image>>,
+        parameters: &TextureParameterDictionary,
+        normal_map: Option<Arc<HostImage>>,
         _named_materials: &HashMap<String, Material>,
         _loc: &FileLoc,
-        _arena: &Arena,
+        arena: &Arena,
     ) -> Result<Material> {
-        todo!()
+        let reflectance = parameters
+            .get_spectrum_texture("reflectance", None, SpectrumType::Albedo)
+            .unwrap_or_else(|| {
+                Arc::new(SpectrumTexture::Constant(
+                    SpectrumConstantTexture::new(Spectrum::Constant(ConstantSpectrum::new(0.5))),
+                ))
+            });
+        let displacement = parameters.get_float_texture_or_null("displacement")?;
+
+        let specific = DiffuseMaterial {
+            reflectance: arena.upload(reflectance),
+            displacement: arena.upload(displacement),
+            normal_map: arena.upload(normal_map),
+        };
+        Ok(Material::Diffuse(specific))
     }
 }
 
 impl CreateMaterial for DiffuseTransmissionMaterial {
     fn create(
         _parameters: &TextureParameterDictionary,
-        _normal_map: Option<Arc<Image>>,
+        _normal_map: Option<Arc<HostImage>>,
         _named_materials: &HashMap<String, Material>,
         _loc: &FileLoc,
         _arena: &Arena,

src/materials/mix.rs

@@ -1,4 +1,4 @@
-use crate::core::image::Image;
+use crate::core::image::HostImage;
 use crate::core::material::CreateMaterial;
 use crate::utils::{Arena, FileLoc, TextureParameterDictionary};
 use anyhow::Result;
@@ -10,7 +10,7 @@ use std::sync::Arc;
 impl CreateMaterial for MixMaterial {
     fn create(
         _parameters: &TextureParameterDictionary,
-        _normal_map: Option<Arc<Image>>,
+        _normal_map: Option<Arc<HostImage>>,
         _named_materials: &HashMap<String, Material>,
         _loc: &FileLoc,
         _arena: &Arena,

src/shapes/bilinear.rs

@@ -1,13 +1,12 @@
-use crate::core::image::{Image, ImageIO};
+use crate::core::image::{HostImage, ImageIO};
 use crate::core::shape::{CreateShape, ALL_BILINEAR_MESHES};
 use crate::core::texture::FloatTexture;
-use crate::shapes::mesh::BilinearPatchMesh;
-use crate::utils::sampling::PiecewiseConstant2D;
 use crate::{Arena, FileLoc, ParameterDictionary};
 use anyhow::{anyhow, Result};
 use log::warn;
 use shared::core::shape::Shape;
-use shared::shapes::BilinearPatchShape;
+use shared::shapes::{BilinearPatchMesh, BilinearPatchShape};
+use shared::utils::sampling::PiecewiseConstant2D;
 use shared::{Ptr, Transform};
 use std::collections::HashMap;
 use std::path::Path;
@@ -93,10 +92,10 @@ impl CreateShape for BilinearPatchShape {
                 );
                 None
             } else {
-                let im = Image::read(Path::new(&filename), None)?;
+                let im = HostImage::read(Path::new(&filename), None)?;
                 let mut img = im.image;
                 img.flip_y();
-                Some(PiecewiseConstant2D::from_image(&img))
+                Some(PiecewiseConstant2D::from_image(&img.inner))
             }
         } else {
             None
@@ -105,11 +104,11 @@ impl CreateShape for BilinearPatchShape {
         let host = BilinearPatchMesh::new(
             &render_from_object,
             reverse_orientation,
-            vertex_indices,
-            p,
-            n,
-            uv,
-            image_dist,
+            &vertex_indices,
+            &p,
+            &n,
+            &uv,
+            image_dist.as_ref(),
         );
 
         let host_arc = Arc::new(host);
@@ -117,8 +116,8 @@ impl CreateShape for BilinearPatchShape {
         // let mesh_index = global_store.len() as u32;
         global_store.push(host_arc.clone());
         drop(global_store);
-        let n_patches = host_arc.device.n_patches;
-        let mesh_ptr = Ptr::from(&host_arc.device);
+        let n_patches = host_arc.n_patches;
+        let mesh_ptr = arena.alloc_arc(host_arc);
         let mut shapes = Vec::with_capacity(n_patches as usize);
         for i in 0..n_patches as i32 {
             shapes.push(arena.alloc(Shape::BilinearPatch(BilinearPatchShape {

src/shapes/mesh.rs

@@ -1,11 +1,10 @@
-use crate::utils::sampling::PiecewiseConstant2D;
-use crate::Arena;
 use anyhow::{bail, Context, Result as AnyResult};
 use ply_rs::parser::Parser;
 use ply_rs::ply::{DefaultElement, Property};
-use shared::core::geometry::{Normal3f, Point2f, Point3f, Vector3f, VectorLike};
-use shared::shapes::mesh::{BilinearPatchMesh, TriangleMesh};
-use shared::{Ptr, Transform};
+use shared::core::geometry::{Normal3f, Point2f, Point3f, VectorLike};
+use shared::shapes::mesh::TriangleMesh;
+use shared::utils::sampling::PiecewiseConstant2D;
+use shared::Transform;
 use std::fs::File;
 use std::path::Path;
 
@@ -269,11 +268,13 @@ impl TriQuadMesh {
 }
 
 pub trait ReadTriangleMesh {
-    pub fn from_ply<P: AsRef<Path>>(
+    fn from_ply<P: AsRef<Path>>(
         filename: P,
         render_from_object: &Transform,
         reverse_orientation: bool,
-    ) -> AnyResult<Self>;
+    ) -> AnyResult<Self>
+    where
+        Self: Sized;
 }
 
 impl ReadTriangleMesh for TriangleMesh {

src/shapes/triangle.rs

@@ -1,11 +1,10 @@
 use crate::core::shape::{CreateShape, ALL_TRIANGLE_MESHES};
 use crate::core::texture::FloatTexture;
-use crate::shapes::mesh::TriangleMesh;
 use crate::{Arena, FileLoc, ParameterDictionary};
 use anyhow::{bail, Result};
 use log::warn;
 use shared::core::shape::Shape;
-use shared::shapes::TriangleShape;
+use shared::shapes::{TriangleMesh, TriangleShape};
 use shared::{Ptr, Transform};
 use std::collections::HashMap;
 use std::sync::Arc;
@@ -87,12 +86,12 @@ impl CreateShape for TriangleShape {
         let host = TriangleMesh::new(
             &render_from_object,
             reverse_orientation,
-            vertex_indices,
-            p,
-            n,
-            s,
-            uvs,
-            face_indices,
+            &vertex_indices,
+            &p,
+            &n,
+            &s,
+            &uvs,
+            &face_indices,
         );
 
         let host_arc = Arc::new(host);
@@ -100,8 +99,8 @@ impl CreateShape for TriangleShape {
         // let mesh_index = global_store.len() as u32;
         global_store.push(host_arc.clone());
         drop(global_store);
-        let n_patches = host_arc.device.n_triangles;
-        let mesh_ptr = Ptr::from(&host_arc.device);
+        let n_patches = host_arc.n_triangles;
+        let mesh_ptr = arena.alloc_arc(host_arc);
         let mut shapes = Vec::with_capacity(n_patches as usize);
         for i in 0..n_patches {
             shapes.push(arena.alloc(Shape::Triangle(TriangleShape {

src/spectra/colorspace.rs

@@ -12,8 +12,8 @@ pub trait CreateRGBColorSpace {
         r: Point2f,
         g: Point2f,
         b: Point2f,
-        illuminant: &DenselySampledSpectrum,
-        rgb_to_spectrum_table: &RGBToSpectrumTable,
+        illuminant: Ptr<DenselySampledSpectrum>,
+        rgb_to_spectrum_table: Ptr<RGBToSpectrumTable>,
     ) -> Self;
 }
 
@@ -22,19 +22,16 @@ impl CreateRGBColorSpace for RGBColorSpace {
         r: Point2f,
         g: Point2f,
         b: Point2f,
-        illuminant: &DenselySampledSpectrum,
-        rgb_to_spectrum_table: &RGBToSpectrumTable,
+        illuminant: Ptr<DenselySampledSpectrum>,
+        rgb_to_spectrum_table: Ptr<RGBToSpectrumTable>,
     ) -> Self {
         let stdspec = get_spectra_context();
-        let illum_ptr = Ptr::from(illuminant);
-        let illum_spectrum = Spectrum::Dense(illum_ptr);
+        let illum_spectrum = Spectrum::Dense(illuminant);
         let w_xyz: XYZ = illum_spectrum.to_xyz(&stdspec);
         let w = w_xyz.xy();
-
         let r_xyz = XYZ::from_xyy(r, Some(1.0));
         let g_xyz = XYZ::from_xyy(g, Some(1.0));
         let b_xyz = XYZ::from_xyy(b, Some(1.0));
-
         let rgb_values = [
             [r_xyz.x(), g_xyz.x(), b_xyz.x()],
             [r_xyz.y(), g_xyz.y(), b_xyz.y()],
@@ -44,14 +41,13 @@ impl CreateRGBColorSpace for RGBColorSpace {
         let c: RGB = rgb.inverse().unwrap() * w_xyz;
         let xyz_from_rgb = rgb * SquareMatrix::diag(&[c.r, c.g, c.b]);
         let rgb_from_xyz = xyz_from_rgb.inverse().expect("singular");
-
         RGBColorSpace {
             r,
             g,
             b,
             w,
-            illuminant: illum_ptr,
-            rgb_to_spectrum_table: Ptr::from(rgb_to_spectrum_table),
+            illuminant,
+            rgb_to_spectrum_table,
             xyz_from_rgb,
             rgb_from_xyz,
         }

src/spectra/data.rs

@@ -1,7 +1,7 @@
 use shared::core::spectrum::Spectrum;
 use shared::spectra::cie::*;
 use shared::spectra::{DenselySampledSpectrum, PiecewiseLinearSpectrum};
-use shared::{gbox, leak, Float, Ptr};
+use shared::{gbox, gvec_from_slice, leak, Float, Ptr};
 use std::collections::HashMap;
 use std::sync::LazyLock;
 
@@ -11,13 +11,11 @@ pub fn create_cie(data: &[Float]) -> DenselySampledSpectrum {
         95 => (300.0, 5.0),
         n => panic!("Unexpected CIE data length: {}", n),
     };
-
     let lambdas: Vec<Float> = (0..data.len())
         .map(|i| start_lambda + i as Float * step)
         .collect();
-
-    let buffer = PiecewiseLinearSpectrum::new(lambdas, data.to_vec());
-    let spec = Spectrum::Piecewise(Ptr::from(&*buffer));
+    let buffer = PiecewiseLinearSpectrum::new(gvec_from_slice(&lambdas), gvec_from_slice(data));
+    let spec = Spectrum::Piecewise(leak(buffer));
     DenselySampledSpectrum::from_spectrum(&spec)
 }
 

src/spectra/mod.rs

@@ -6,8 +6,7 @@ use shared::core::spectrum::{Spectrum, StandardSpectra};
 use shared::spectra::cie::{CIE_D65, CIE_X, CIE_Y, CIE_Z};
 use shared::spectra::{DenselySampledSpectrum, DeviceStandardColorSpaces, RGBColorSpace};
 use shared::Ptr;
-use std::sync::Arc;
-use std::sync::LazyLock;
+use std::sync::{Arc, LazyLock, OnceLock};
 
 pub mod colorspace;
 pub mod data;
@@ -22,10 +21,6 @@ pub static CIE_Z_DATA: LazyLock<DenselySampledSpectrum> =
 pub static CIE_D65_DATA: LazyLock<DenselySampledSpectrum> =
     LazyLock::new(|| data::create_cie(&CIE_D65));
 
-fn get_d65_illuminant_buffer() -> Arc<DenselySampledSpectrum> {
-    Arc::new(CIE_D65_DATA.clone())
-}
-
 pub fn cie_x() -> Spectrum {
     Spectrum::Dense(Ptr::from(&*CIE_X_DATA))
 }
@@ -48,50 +43,54 @@ pub fn get_spectra_context() -> StandardSpectra {
     }
 }
 
+static D65_ILLUMINANT: LazyLock<DenselySampledSpectrum> = LazyLock::new(|| CIE_D65_DATA.clone());
 
 pub static SRGB: LazyLock<Arc<RGBColorSpace>> = LazyLock::new(|| {
-    let illum = get_d65_illuminant_buffer();
     let r = Point2f::new(0.64, 0.33);
     let g = Point2f::new(0.3, 0.6);
     let b = Point2f::new(0.15, 0.06);
     let table_ptr = Ptr::from(&*SRGB_TABLE);
+    let illum_ptr = Ptr::from(&*D65_ILLUMINANT);
 
-    Arc::new(RGBColorSpace::new(r, g, b, &illum, &table_ptr))
+    Arc::new(RGBColorSpace::new(r, g, b, illum_ptr, table_ptr))
 });
 
 pub static DCI_P3: LazyLock<Arc<RGBColorSpace>> = LazyLock::new(|| {
-    let illum = get_d65_illuminant_buffer();
     let r = Point2f::new(0.680, 0.320);
     let g = Point2f::new(0.265, 0.690);
     let b = Point2f::new(0.150, 0.060);
     let table_ptr = Ptr::from(&*DCI_P3_TABLE);
-    Arc::new(RGBColorSpace::new(r, g, b, &illum, &table_ptr))
+    let illum_ptr = Ptr::from(&*D65_ILLUMINANT);
+
+    Arc::new(RGBColorSpace::new(r, g, b, illum_ptr, table_ptr))
 });
 
-pub static REC2020: LazyLock<Arc<RGBColorSpaceData>> = LazyLock::new(|| {
-    let illum = get_d65_illuminant_buffer();
+pub static REC2020: LazyLock<Arc<RGBColorSpace>> = LazyLock::new(|| {
     let r = Point2f::new(0.708, 0.292);
     let g = Point2f::new(0.170, 0.797);
     let b = Point2f::new(0.131, 0.046);
     let table_ptr = Ptr::from(&*REC2020_TABLE);
-    Arc::new(RGBColorSpace::new(r, g, b, &illum, &table_ptr))
+    let illum_ptr = Ptr::from(&*D65_ILLUMINANT);
+
+    Arc::new(RGBColorSpace::new(r, g, b, illum_ptr, table_ptr))
 });
 
-pub static ACES: LazyLock<Arc<RGBColorSpaceData>> = LazyLock::new(|| {
-    let illum = get_d65_illuminant_buffer();
+pub static ACES: LazyLock<Arc<RGBColorSpace>> = LazyLock::new(|| {
     let r = Point2f::new(0.7347, 0.2653);
     let g = Point2f::new(0.0000, 1.0000);
     let b = Point2f::new(0.0001, -0.0770);
-    let table_ptr = Ptr::from(&ACES_TABLE);
-    Arc::new(RGBColorSpace::new(r, g, b, &illum, &table_ptr))
+    let table_ptr = Ptr::from(&*ACES_TABLE);
+    let illum_ptr = Ptr::from(&*D65_ILLUMINANT);
+
+    Arc::new(RGBColorSpace::new(r, g, b, illum_ptr, table_ptr))
 });
 
 #[derive(Debug, Clone)]
 pub struct StandardColorSpaces {
-    pub srgb: Arc<RGBColorSpaceData>,
-    pub dci_p3: Arc<RGBColorSpaceData>,
-    pub rec2020: Arc<RGBColorSpaceData>,
-    pub aces2065_1: Arc<RGBColorSpaceData>,
+    pub srgb: Arc<RGBColorSpace>,
+    pub dci_p3: Arc<RGBColorSpace>,
+    pub rec2020: Arc<RGBColorSpace>,
+    pub aces2065_1: Arc<RGBColorSpace>,
 }
 
 impl StandardColorSpaces {
@@ -135,7 +134,8 @@ pub fn default_colorspace_arc() -> Arc<RGBColorSpace> {
 
 pub fn default_colorspace_ref() -> &'static RGBColorSpace {
     static CS: OnceLock<RGBColorSpace> = OnceLock::new();
-    CS.get_or_init(|| stdcs.srgb)
+    let stdcs = get_colorspace_device();
+    CS.get_or_init(|| *stdcs.srgb)
 }
 
 pub fn default_illuminant() -> Spectrum {

src/textures/mix.rs

@@ -7,10 +7,6 @@ use anyhow::Result;
 use shared::core::geometry::{Vector3f, VectorLike};
 use shared::core::texture::{SpectrumType, TextureEvalContext};
 use shared::spectra::{SampledSpectrum, SampledWavelengths};
-use shared::textures::{
-    GPUFloatDirectionMixTexture, GPUFloatMixTexture, GPUSpectrumDirectionMixTexture,
-    GPUSpectrumMixTexture,
-};
 use shared::utils::Transform;
 use shared::Float;
 use std::sync::Arc;

src/textures/scaled.rs

@@ -5,7 +5,6 @@ use crate::Arena;
 use anyhow::Result;
 use shared::core::texture::{SpectrumType, TextureEvalContext};
 use shared::spectra::{SampledSpectrum, SampledWavelengths};
-use shared::textures::{GPUFloatScaledTexture, GPUSpectrumScaledTexture};
 use shared::utils::Transform;
 use shared::Float;
 use std::sync::Arc;

src/utils/arena.rs

@@ -148,6 +148,13 @@ impl<A: GpuAllocator> Arena<A> {
         }
     }
 
+    pub fn alloc_arc<T: Clone>(&self, value: Arc<T>) -> Ptr<T> {
+        match Arc::try_unwrap(value) {
+            Ok(inner) => self.alloc(inner),
+            Err(arc) => self.alloc((*arc).clone()),
+        }
+    }
+
     pub fn alloc_slice<T: Copy>(&self, values: &[T]) -> (Ptr<T>, usize) {
         let mut bump = self.bump.lock();
         let (ptr, len) = bump.alloc_slice(values);

src/utils/mipmap.rs

@@ -1,10 +1,10 @@
-use crate::core::image::{Image, ImageIO};
-use shared::Float;
+use crate::core::image::{HostImage, ImageIO};
 use shared::core::color::{ColorEncoding, RGB};
 use shared::core::geometry::{Point2f, Point2i, Vector2f, VectorLike};
 use shared::core::image::{WrapMode, WrapMode2D};
 use shared::spectra::RGBColorSpace;
 use shared::utils::math::{lerp, safe_sqrt, square};
+use shared::Float;
 use std::hash::{Hash, Hasher};
 use std::ops::{Add, Mul, Sub};
 use std::path::Path;
@@ -69,18 +69,18 @@ pub trait MIPMapSample:
     Copy + Add<Output = Self> + Sub<Output = Self> + Mul<Float, Output = Self> + std::fmt::Debug
 {
     fn zero() -> Self;
-    fn sample_bilerp(image: &Image, st: Point2f, wrap: WrapMode2D) -> Self;
-    fn sample_texel(image: &Image, st: Point2i, wrap: WrapMode2D) -> Self;
+    fn sample_bilerp(image: &HostImage, st: Point2f, wrap: WrapMode2D) -> Self;
+    fn sample_texel(image: &HostImage, st: Point2i, wrap: WrapMode2D) -> Self;
 }
 
 impl MIPMapSample for Float {
     fn zero() -> Self {
         0.
     }
-    fn sample_bilerp(image: &Image, st: Point2f, wrap: WrapMode2D) -> Self {
+    fn sample_bilerp(image: &HostImage, st: Point2f, wrap: WrapMode2D) -> Self {
         image.bilerp_channel_with_wrap(st, 0, wrap)
     }
-    fn sample_texel(image: &Image, st: Point2i, wrap: WrapMode2D) -> Self {
+    fn sample_texel(image: &HostImage, st: Point2i, wrap: WrapMode2D) -> Self {
         image.get_channel_with_wrap(st, 0, wrap)
     }
 }
@@ -89,7 +89,7 @@ impl MIPMapSample for RGB {
     fn zero() -> Self {
         RGB::new(0., 0., 0.)
     }
-    fn sample_bilerp(image: &Image, st: Point2f, wrap: WrapMode2D) -> Self {
+    fn sample_bilerp(image: &HostImage, st: Point2f, wrap: WrapMode2D) -> Self {
         let nc = image.n_channels();
         if nc >= 3 {
             let r = image.bilerp_channel_with_wrap(st, 0, wrap);
@@ -101,7 +101,7 @@ impl MIPMapSample for RGB {
             RGB::new(v, v, v)
         }
     }
-    fn sample_texel(image: &Image, st: Point2i, wrap: WrapMode2D) -> Self {
+    fn sample_texel(image: &HostImage, st: Point2i, wrap: WrapMode2D) -> Self {
         let nc = image.n_channels();
         if nc >= 3 {
             let r = image.get_channel_with_wrap(st, 0, wrap);
@@ -117,7 +117,7 @@ impl MIPMapSample for RGB {
 
 #[derive(Clone, Debug)]
 pub struct MIPMap {
-    pub pyramid: Vec<Image>,
+    pub pyramid: Vec<HostImage>,
     pub color_space: Option<RGBColorSpace>,
     pub wrap_mode: WrapMode,
     pub options: MIPMapFilterOptions,
@@ -127,12 +127,12 @@ pub struct MIPMap {
 
 impl MIPMap {
     pub fn new(
-        image: Image,
+        image: HostImage,
         color_space: Option<RGBColorSpace>,
         wrap_mode: WrapMode,
         options: MIPMapFilterOptions,
     ) -> Self {
-        let pyramid = Image::generate_pyramid(image, wrap_mode);
+        let pyramid = HostImage::generate_pyramid(image, wrap_mode);
         Self {
             pyramid,
             color_space,
@@ -160,11 +160,11 @@ impl MIPMap {
         self.color_space.clone()
     }
 
-    pub fn get_level(&self, level: usize) -> &Image {
+    pub fn get_level(&self, level: usize) -> &HostImage {
         &self.pyramid[level]
     }
 
-    pub fn base_image(&self) -> &Image {
+    pub fn base_image(&self) -> &HostImage {
         &self.pyramid[0]
     }
 
@@ -321,7 +321,7 @@ impl MIPMap {
         wrap_mode: WrapMode,
         encoding: ColorEncoding,
     ) -> Result<MIPMap, ()> {
-        let image_and_metadata = Image::read(filename, Some(encoding)).unwrap();
+        let image_and_metadata = HostImage::read(filename, Some(encoding)).unwrap();
         let image = image_and_metadata.image;
         Ok(MIPMap::new(
             image,
@@ -345,7 +345,7 @@ impl MIPMap {
 }
 
 #[cfg(feature = "cuda")]
-fn create_cuda_texture(pyramid: &[Image], wrap_mode: WrapMode) -> u64 {
+fn create_cuda_texture(pyramid: &[HostImage], wrap_mode: WrapMode) -> u64 {
     use cuda_runtime_sys::*;
 
     let base = &pyramid[0];

src/utils/mod.rs

@@ -8,7 +8,6 @@ pub mod mipmap;
 pub mod parallel;
 pub mod parameters;
 pub mod parser;
-pub mod sampling;
 pub mod strings;
 pub mod upload;
 
@@ -18,6 +17,7 @@ pub use parameters::{
     ParameterDictionary, ParsedParameter, ParsedParameterVector, TextureParameterDictionary,
 };
 pub use strings::*;
+pub use mipmap::{MIPMap, MIPMapFilterOptions};
 pub use upload::{Upload, ArenaUpload};
 
 #[cfg(feature = "vulkan")]

src/utils/parameters.rs

@@ -12,7 +12,7 @@ use shared::spectra::{
     PiecewiseLinearSpectrum, RGBAlbedoSpectrum, RGBColorSpace, RGBIlluminantSpectrum,
     RGBUnboundedSpectrum,
 };
-use shared::Float;
+use shared::{gvec_from_slice, leak, Float};
 
 use std::collections::HashMap;
 use std::sync::{
@@ -682,8 +682,8 @@ impl ParameterDictionary {
             .unzip();
 
         vec![Spectrum::Piecewise(leak(PiecewiseLinearSpectrum {
-            lambdas: gvec_from_slice(lambdas),
-            values: gvec_from_slice(values),
+            lambdas: gvec_from_slice(&lambdas),
+            values: gvec_from_slice(&values),
             count: lambdas.len() as u32,
         }))]
     }
@@ -915,7 +915,12 @@ impl TextureParameterDictionary {
                         panic!("[{:?}] Negative RGB values for '{}'", p.loc, name);
                     }
 
-                    let cs = self.dict.color_space.as_ref().unwrap();
+                    let cs = self
+                        .dict
+                        .color_space
+                        .as_ref()
+                        .map(|arc| arc.as_ref())
+                        .unwrap_or_else(|| crate::spectra::default_colorspace_ref());
                     let s: Spectrum = match stype {
                         SpectrumType::Illuminant => {
                             Spectrum::RGBIlluminant(RGBIlluminantSpectrum::new(cs, rgb))

src/utils/upload.rs

@@ -1,4 +1,6 @@
 use crate::core::texture::{FloatTexture, SpectrumTexture};
+use crate::core::image::HostImage;
+use shared::core::image::Image;
 use crate::Arena;
 use shared::core::texture::{GPUFloatTexture, GPUSpectrumTexture};
 use shared::textures::*;
@@ -149,6 +151,27 @@ impl Upload for Option<Arc<SpectrumTexture>> {
     }
 }
 
+impl Upload for HostImage {
+    type Target = Ptr<Image>;
+    fn upload(self, arena: &Arena) -> Ptr<Image> {
+        arena.alloc(self.inner)
+    }
+}
+
+impl Upload for Option<HostImage> {
+    type Target = Ptr<Image>;
+    fn upload(self, arena: &Arena) -> Ptr<Image> {
+        arena.alloc_opt(self.map(|h| h.inner))
+    }
+}
+
+impl Upload for Option<Arc<HostImage>> {
+    type Target = Ptr<Image>;
+    fn upload(self, arena: &Arena) -> Ptr<Image> {
+        arena.alloc_opt(self.map(|h| h.as_ref().inner.clone()))
+    }
+}
+
 pub trait ArenaUpload {
     fn upload<T: Upload>(&self, value: T) -> T::Target;
 }