Added missing geometry. Implementing missing material methods

shared/src/bxdfs/diffuse.rs

@@ -7,7 +7,7 @@ use crate::{Float, INV_PI};
 use core::any::Any;
 
 #[repr(C)]
-#[derive(Debug, Copy, Clone)]
+#[derive(Debug, Copy, Clone, Default)]
 pub struct DiffuseBxDF {
     pub r: SampledSpectrum,
 }

shared/src/cameras/perspective.rs

@@ -115,10 +115,11 @@ impl CameraTrait for PerspectiveCamera {
             r.d = (p_focus - r.o).normalize();
         }
 
-        let ray = self.render_from_camera(&r, &mut None);
+        let mut ray = self.render_from_camera(&r, &mut None);
+        ray.d = ray.d.normalize();
         Some(CameraRay {
             ray,
-            weight: SampledSpectrum::default(),
+            weight: SampledSpectrum::new(1.),
         })
     }
 }

shared/src/core/bsdf.rs

@@ -1,18 +1,18 @@
-use crate::Float;
 use crate::core::bxdf::{BxDF, BxDFFlags, BxDFTrait, FArgs, TransportMode};
 use crate::core::geometry::{Frame, Normal3f, Point2f, Vector3f, VectorLike};
 use crate::spectra::SampledSpectrum;
 use crate::utils::Ptr;
+use crate::Float;
 
 #[repr(C)]
 #[derive(Copy, Clone, Debug, Default)]
 pub struct BSDF {
-    bxdf: Ptr<BxDF>,
+    bxdf: BxDF,
     shading_frame: Frame,
 }
 
 impl BSDF {
-    pub fn new(ns: Normal3f, dpdus: Vector3f, bxdf: Ptr<BxDF>) -> Self {
+    pub fn new(ns: Normal3f, dpdus: Vector3f, bxdf: BxDF) -> Self {
         Self {
             bxdf,
             shading_frame: Frame::new(dpdus.normalize(), Vector3f::from(ns)),
@@ -20,11 +20,11 @@ impl BSDF {
     }
 
     pub fn is_valid(&self) -> bool {
-        !self.bxdf.is_null()
+        !self.bxdf.flags().is_empty()
     }
 
     pub fn flags(&self) -> BxDFFlags {
-        if self.bxdf.is_null() {
+        if !self.is_valid() {
             // Either this, or transmissive for seethrough
             return BxDFFlags::empty();
         }
@@ -45,7 +45,7 @@ impl BSDF {
         wi_render: Vector3f,
         mode: TransportMode,
     ) -> Option<SampledSpectrum> {
-        if self.bxdf.is_null() {
+        if !self.is_valid() {
             return None;
         }
 
@@ -66,7 +66,7 @@ impl BSDF {
         u2: Point2f,
         f_args: FArgs,
     ) -> Option<BSDFSample> {
-        let bxdf = unsafe { self.bxdf.as_ref() };
+        let bxdf = self.bxdf;
 
         let sampling_flags = BxDFFlags::from_bits_truncate(f_args.sample_flags.bits());
         let wo = self.render_to_local(wo_render);
@@ -85,7 +85,7 @@ impl BSDF {
     }
 
     pub fn pdf(&self, wo_render: Vector3f, wi_render: Vector3f, f_args: FArgs) -> Float {
-        if self.bxdf.is_null() {
+        if !self.is_valid() {
             return 0.0;
         }
         let sample_flags = BxDFFlags::from_bits_truncate(f_args.sample_flags.bits());
@@ -101,7 +101,7 @@ impl BSDF {
     }
 
     pub fn rho_u(&self, u1: &[Point2f], uc: &[Float], u2: &[Point2f]) -> SampledSpectrum {
-        if self.bxdf.is_null() {
+        if !self.is_valid() {
             return SampledSpectrum::default();
         }
 
@@ -109,7 +109,7 @@ impl BSDF {
     }
 
     pub fn rho_wo(&self, wo_render: Vector3f, uc: &[Float], u: &[Point2f]) -> SampledSpectrum {
-        if self.bxdf.is_null() {
+        if !self.is_valid() {
             return SampledSpectrum::default();
         }
 
@@ -118,8 +118,8 @@ impl BSDF {
     }
 
     pub fn regularize(&mut self) {
-        if !self.bxdf.is_null() {
+        if !self.is_valid() {
-            let mut bxdf = unsafe { *self.bxdf.as_raw() };
+            let bxdf = &mut self.bxdf;
             bxdf.regularize();
         }
     }

shared/src/core/bxdf.rs

@@ -155,3 +155,9 @@ pub enum BxDF {
     CoatedConductor(CoatedConductorBxDF),
     NormalizedFresnel(NormalizedFresnelBxDF),
 }
+
+impl Default for BxDF {
+    fn default() -> Self {
+        BxDF::Diffuse(DiffuseBxDF::default())
+    }
+}

shared/src/core/camera.rs

@@ -41,6 +41,10 @@ pub struct CameraTransform {
 }
 
 impl CameraTransform {
+    pub fn render_from_world(&self) -> Transform {
+        self.world_from_render.inverse()
+    }
+
     pub fn from_world(
         world_from_camera: AnimatedTransform,
         rendering_space: RenderingCoordinateSystem,

shared/src/core/film.rs

@@ -115,6 +115,9 @@ impl RGBFilm {
         _vi: Option<&VisibleSurface>,
         weight: Float,
     ) {
+        if !self.base.pixel_bounds.contains_exclusive(p_film) {
+            return;
+        }
         let sensor = self.get_sensor();
         let mut rgb = sensor.to_sensor_rgb(l, lambda);
         let m = rgb.into_iter().copied().fold(f32::NEG_INFINITY, f32::max);

shared/src/core/filter.rs

@@ -75,17 +75,17 @@ pub enum Filter {
 
 impl<T: FilterTrait> FilterTrait for Ptr<T> {
     fn radius(&self) -> Vector2f {
-        unsafe { self.as_ref().radius() }
+        self.get().unwrap().radius()
     }
     fn integral(&self) -> Float {
-        unsafe { self.as_ref().integral() }
+        self.get().unwrap().integral()
     }
  
     fn evaluate(&self, p: Point2f) -> Float {
-        unsafe { self.as_ref().evaluate(p) }
+        self.get().unwrap().evaluate(p)
     }
 
     fn sample(&self, p: Point2f) -> FilterSample {
-        unsafe { self.as_ref().sample(p) }
+        self.get().unwrap().sample(p)
     }
 }

shared/src/core/geometry/bounds.rs

@@ -1,5 +1,5 @@
 use super::{Float, NumFloat};
-use super::{Point, Point2f, Point3, Point3f, Vector, Vector2, Vector2f, Vector3, Vector3f};
+use super::{Point,  Point2i, Point2f, Point3, Point3f, Vector, Vector2, Vector2f, Vector3, Vector3f};
 use crate::core::geometry::traits::{SqrtExt, VectorLike};
 use crate::core::geometry::{max, min};
 use crate::utils::gpu_array_from_fn;
@@ -349,3 +349,35 @@ impl Bounds3f {
         (t_min < ray_t_max) && (t_max > 0.0)
     }
 }
+
+pub struct BoundsPixelIterator {
+    bounds: Bounds2i,
+    current: Point2i,
+}
+
+impl Iterator for BoundsPixelIterator {
+    type Item = Point2i;
+    fn next(&mut self) -> Option<Point2i> {
+        if self.current.y() >= self.bounds.p_max.y() {
+            return None;
+        }
+        let result = self.current;
+        let mut x = self.current.x() + 1;
+        let mut y = self.current.y();
+        if x >= self.bounds.p_max.x() {
+            x = self.bounds.p_min.x();
+            y += 1;
+        }
+        self.current = Point2i::new(x, y);
+        Some(result)
+    }
+}
+
+impl Bounds2i {
+    pub fn pixels(&self) -> BoundsPixelIterator {
+        BoundsPixelIterator {
+            bounds: *self,
+            current: self.p_min,
+        }
+    }
+}

shared/src/core/interaction.rs

@@ -354,8 +354,7 @@ impl SurfaceInteraction {
         let (dpdu, dpdv) = if !displacement.is_null() {
             bump_map(tex_eval, &displacement, &ctx)
         } else if !normal_image.is_null() {
-            let map = unsafe { normal_image.as_ref() };
+            normal_map(&normal_image, &ctx)
-            normal_map(map, &ctx)
         } else {
             (self.shading.dpdu, self.shading.dpdv)
         };
@@ -590,13 +589,13 @@ impl SurfaceInteraction {
     #[cfg(not(target_os = "cuda"))]
     pub fn set_intersection_properties(
         &mut self,
-        mtl: &Material,
+        mtl: Ptr<Material>,
-        area: &Light,
+        area: Ptr<Light>,
-        ray_medium: &Medium,
+        ray_medium: Ptr<Medium>,
         prim_medium_interface: MediumInterface,
     ) {
-        self.material = Ptr::from(mtl);
+        self.material = mtl;
-        self.area_light = Ptr::from(area);
+        self.area_light = area;
 
         if prim_medium_interface.is_medium_transition() {
             self.common.medium_interface = prim_medium_interface;

shared/src/core/material.rs

@@ -2,7 +2,6 @@ use crate::materials::*;
 use core::ops::Deref;
 use enum_dispatch::enum_dispatch;
 
-use crate::Float;
 use crate::bxdfs::{
     CoatedConductorBxDF, CoatedDiffuseBxDF, ConductorBxDF, DielectricBxDF, DiffuseBxDF,
 };
@@ -19,9 +18,10 @@ use crate::core::texture::{
 };
 use crate::materials::*;
 use crate::spectra::{SampledSpectrum, SampledWavelengths};
-use crate::utils::Ptr;
 use crate::utils::hash::hash_float;
 use crate::utils::math::clamp;
+use crate::utils::Ptr;
+use crate::Float;
 
 #[repr(C)]
 #[derive(Clone, Debug, Copy)]
@@ -194,7 +194,6 @@ pub enum Material {
     Mix(MixMaterial),
 }
 
-
 // TODO: THIS IS A HACK JUST FOR TESTING
 impl PartialEq for Material {
     fn eq(&self, other: &Self) -> bool {

shared/src/core/primitive.rs

@@ -43,7 +43,7 @@ impl PrimitiveTrait for GeometricPrimitive {
     fn intersect(&self, r: &Ray, t_max: Option<Float>) -> Option<ShapeIntersection> {
         let mut si = self.shape.intersect(r, t_max)?;
         if !self.alpha.is_null() {
-            let alpha = unsafe { &self.alpha.as_ref() };
+            let alpha = &self.alpha.get().unwrap();
             let ctx = TextureEvalContext::from(&si.intr);
             let a = alpha.evaluate(&ctx);
             if a < 1.0 {
@@ -224,14 +224,14 @@ pub enum Primitive {
 
 impl<T: PrimitiveTrait> PrimitiveTrait for Ptr<T> {
     fn bounds(&self) -> Bounds3f {
-        unsafe { self.as_ref().bounds() }
+        self.get().unwrap().bounds()
     }
 
     fn intersect(&self, r: &Ray, t_max: Option<Float>) -> Option<ShapeIntersection> {
-        unsafe { self.as_ref().intersect(r, t_max) }
+        self.get().unwrap().intersect(r, t_max)
     }
 
     fn intersect_p(&self, r: &Ray, t_max: Option<Float>) -> bool {
-        unsafe { self.as_ref().intersect_p(r, t_max) }
+        self.get().unwrap().intersect_p(r, t_max)
     }
 }

shared/src/core/shape.rs

@@ -43,7 +43,7 @@ impl ShapeIntersection {
         ray_medium: Ptr<Medium>,
     ) {
         self.intr
-            .set_intersection_properties(&mtl, &area, &ray_medium, prim_medium_interface);
+            .set_intersection_properties(mtl, area, ray_medium, prim_medium_interface);
     }
 }
 

shared/src/core/spectrum.rs

@@ -40,10 +40,10 @@ pub enum Spectrum {
 
 impl<T: SpectrumTrait> SpectrumTrait for Ptr<T> {
     fn evaluate(&self, lambda: Float) -> Float {
-        unsafe { self.as_ref().evaluate(lambda) }
+        self.get().unwrap().evaluate(lambda)
     }
     fn max_value(&self) -> Float {
-        unsafe { self.as_ref().max_value() }
+        self.get().unwrap().max_value()
     }
 }
 

shared/src/lights/diffuse.rs

@@ -126,8 +126,7 @@ impl LightTrait for DiffuseAreaLight {
                 rgb[c] = self.image.bilerp_channel(uv, c as i32);
             }
 
-            let cs_ref = unsafe { self.colorspace.as_ref() };
+            let spec = RGBIlluminantSpectrum::new(&self.colorspace, rgb.clamp_zero());
-            let spec = RGBIlluminantSpectrum::new(cs_ref, rgb.clamp_zero());
 
             self.scale * spec.sample(lambda)
         } else {
@@ -150,8 +149,7 @@ impl LightTrait for DiffuseAreaLight {
                         rgb[c] = self.image.get_channel(Point2i::new(x, y), c as i32);
                     }
 
-                    let cs_ref = unsafe { self.colorspace.as_ref() };
+                    l += RGBIlluminantSpectrum::new(&self.colorspace, rgb.clamp_zero()).sample(&lambda);
-                    l += RGBIlluminantSpectrum::new(cs_ref, rgb.clamp_zero()).sample(&lambda);
                 }
             }
             l *= self.scale / (self.image.resolution().x() * self.image.resolution().y()) as Float;

shared/src/materials/coated.rs

@@ -118,7 +118,7 @@ impl MaterialTrait for CoatedDiffuseMaterial {
             self.seed,
         ));
 
-        BSDF::new(ctx.ns, ctx.dpdus, Ptr::from(&bxdf))
+        BSDF::new(ctx.ns, ctx.dpdus, bxdf)
     }
 
     fn get_bssrdf<T>(
@@ -245,7 +245,7 @@ impl MaterialTrait for CoatedConductorMaterial {
         let (mut ce, mut ck) = if !self.conductor_eta.is_null() {
             let k_tex = self.k;
             let ce = tex_eval.evaluate_spectrum(&self.conductor_eta, ctx, lambda);
-            let ck = tex_eval.evaluate_spectrum(unsafe { k_tex.as_ref() }, ctx, lambda);
+            let ck = tex_eval.evaluate_spectrum(k_tex.get().unwrap(), ctx, lambda);
             (ce, ck)
         } else {
             let r = SampledSpectrum::clamp(
@@ -288,7 +288,7 @@ impl MaterialTrait for CoatedConductorMaterial {
             self.n_samples,
             self.seed,
         ));
-        BSDF::new(ctx.ns, ctx.dpdus, Ptr::from(&bxdf))
+        BSDF::new(ctx.ns, ctx.dpdus, bxdf)
     }
 
     fn get_bssrdf<T>(

shared/src/materials/conductor.rs

@@ -10,34 +10,85 @@ use crate::core::scattering::TrowbridgeReitzDistribution;
 use crate::core::spectrum::{Spectrum, SpectrumTrait};
 use crate::core::texture::{GPUFloatTexture, GPUSpectrumTexture, TextureEvaluator};
 use crate::spectra::{SampledSpectrum, SampledWavelengths};
-use crate::utils::Ptr;
 use crate::utils::math::clamp;
+use crate::utils::Ptr;
 
 #[repr(C)]
 #[derive(Clone, Copy, Debug)]
 pub struct ConductorMaterial {
-    pub displacement: Ptr<GPUFloatTexture>,
+    pub normal_map: Ptr<Image>,
+    pub reflectance: Ptr<GPUSpectrumTexture>,
     pub eta: Ptr<GPUSpectrumTexture>,
     pub k: Ptr<GPUSpectrumTexture>,
-    pub reflectance: Ptr<GPUSpectrumTexture>,
     pub u_roughness: Ptr<GPUFloatTexture>,
     pub v_roughness: Ptr<GPUFloatTexture>,
+    pub displacement: Ptr<GPUFloatTexture>,
     pub remap_roughness: bool,
-    pub normal_map: Ptr<Image>,
+}
+
+impl ConductorMaterial {
+    pub fn new(
+        normal_map: Ptr<Image>,
+        reflectance: Ptr<GPUSpectrumTexture>,
+        eta: Ptr<GPUSpectrumTexture>,
+        k: Ptr<GPUSpectrumTexture>,
+        u_roughness: Ptr<GPUFloatTexture>,
+        v_roughness: Ptr<GPUFloatTexture>,
+        displacement: Ptr<GPUFloatTexture>,
+        remap_roughness: bool,
+    ) -> Self {
+        Self {
+            normal_map,
+            reflectance,
+            eta,
+            k,
+            u_roughness,
+            v_roughness,
+            displacement,
+            remap_roughness,
+        }
+    }
 }
 
 impl MaterialTrait for ConductorMaterial {
     fn get_bsdf<T: TextureEvaluator>(
         &self,
-        _tex_eval: &T,
+        tex_eval: &T,
-        _ctx: &MaterialEvalContext,
+        ctx: &MaterialEvalContext,
-        _lambda: &SampledWavelengths,
+        lambda: &SampledWavelengths,
     ) -> BSDF {
-        todo!()
+        let mut u_rough = tex_eval.evaluate_float(&self.u_roughness, ctx);
+        let mut v_rough = tex_eval.evaluate_float(&self.v_roughness, ctx);
+        if self.remap_roughness {
+            u_rough = TrowbridgeReitzDistribution::roughness_to_alpha(u_rough);
+            v_rough = TrowbridgeReitzDistribution::roughness_to_alpha(v_rough);
+        }
+        let (etas, ks) = if !self.eta.is_null() {
+            (
+                tex_eval.evaluate_spectrum(&self.eta, ctx, lambda),
+                tex_eval.evaluate_spectrum(&self.k, ctx, lambda),
+            )
+        } else {
+            let r = SampledSpectrum::clamp(
+                &tex_eval.evaluate_spectrum(&self.reflectance, ctx, lambda),
+                0.,
+                0.9999,
+            );
+            let one_minus_r = SampledSpectrum::new(1.) - r;
+            (
+                SampledSpectrum::new(1.),
+                2. * r.sqrt() / SampledSpectrum::clamp_zero(&one_minus_r).sqrt(),
+            )
+        };
+
+        let distrib = TrowbridgeReitzDistribution::new(u_rough, v_rough);
+        let bxdf = BxDF::Conductor(ConductorBxDF::new(&distrib, etas, ks));
+        BSDF::new(ctx.ns, ctx.dpdus, bxdf)
     }
+
     fn get_bssrdf<T>(
         &self,
-        _tex_eval: &T,
+        tex_eval: &T,
         _ctx: &MaterialEvalContext,
         _lambda: &SampledWavelengths,
     ) -> Option<BSSRDF> {
@@ -51,14 +102,14 @@ impl MaterialTrait for ConductorMaterial {
     }
 
     fn get_normal_map(&self) -> Option<&Image> {
-        todo!()
+        self.normal_map.get()
     }
 
     fn get_displacement(&self) -> Ptr<GPUFloatTexture> {
-        todo!()
+        self.displacement
     }
 
     fn has_subsurface_scattering(&self) -> bool {
-        todo!()
+        false
     }
 }

shared/src/materials/dielectric.rs

@@ -10,8 +10,8 @@ use crate::core::scattering::TrowbridgeReitzDistribution;
 use crate::core::spectrum::{Spectrum, SpectrumTrait};
 use crate::core::texture::{GPUFloatTexture, GPUSpectrumTexture, TextureEvaluator};
 use crate::spectra::{SampledSpectrum, SampledWavelengths};
-use crate::utils::Ptr;
 use crate::utils::math::clamp;
+use crate::Ptr;
 
 #[repr(C)]
 #[derive(Clone, Copy, Debug)]
@@ -51,7 +51,7 @@ impl MaterialTrait for DielectricMaterial {
         let distrib = TrowbridgeReitzDistribution::new(u_rough, v_rough);
         let bxdf = BxDF::Dielectric(DielectricBxDF::new(sampled_eta, distrib));
 
-        BSDF::new(ctx.ns, ctx.dpdus, Ptr::from(&bxdf))
+        BSDF::new(ctx.ns, ctx.dpdus, bxdf)
     }
 
     fn get_bssrdf<T>(

shared/src/materials/diffuse.rs

@@ -11,7 +11,7 @@ use crate::core::scattering::TrowbridgeReitzDistribution;
 use crate::core::spectrum::{Spectrum, SpectrumTrait};
 use crate::core::texture::{GPUFloatTexture, GPUSpectrumTexture, TextureEvaluator};
 use crate::spectra::{SampledSpectrum, SampledWavelengths};
-use crate::utils::Ptr;
+use crate::Ptr;
 use crate::utils::math::clamp;
 
 #[repr(C)]
@@ -31,7 +31,7 @@ impl MaterialTrait for DiffuseMaterial {
     ) -> BSDF {
         let r = tex_eval.evaluate_spectrum(&self.reflectance, ctx, lambda);
         let bxdf = BxDF::Diffuse(DiffuseBxDF::new(r));
-        BSDF::new(ctx.ns, ctx.dpdus, Ptr::from(&bxdf))
+        BSDF::new(ctx.ns, ctx.dpdus, bxdf)
     }
 
     fn get_bssrdf<T>(
@@ -48,7 +48,7 @@ impl MaterialTrait for DiffuseMaterial {
     }
 
     fn get_normal_map(&self) -> Option<&Image> {
-        Some(&*self.normal_map)
+        self.normal_map.get()
     }
 
     fn get_displacement(&self) -> Ptr<GPUFloatTexture> {
@@ -93,7 +93,7 @@ impl MaterialTrait for DiffuseTransmissionMaterial {
     }
 
     fn get_normal_map(&self) -> Option<&Image> {
-        Some(&*self.image)
+        self.image.get()
     }
 
     fn get_displacement(&self) -> Ptr<GPUFloatTexture> {

shared/src/spectra/rgb.rs

@@ -1,8 +1,8 @@
 use super::{
-    DenselySampledSpectrum, LAMBDA_MAX, LAMBDA_MIN, N_SPECTRUM_SAMPLES, RGBColorSpace,
+    DenselySampledSpectrum, RGBColorSpace, SampledSpectrum, SampledWavelengths, LAMBDA_MAX,
-    SampledSpectrum, SampledWavelengths,
+    LAMBDA_MIN, N_SPECTRUM_SAMPLES,
 };
-use crate::core::color::{RGB, RGBSigmoidPolynomial, XYZ};
+use crate::core::color::{RGBSigmoidPolynomial, RGB, XYZ};
 use crate::core::spectrum::SpectrumTrait;
 use crate::utils::Ptr;
 
@@ -77,7 +77,7 @@ impl RGBIlluminantSpectrum {
         let illuminant = cs.illuminant;
         let m = rgb.max_component_value();
         let scale = 2. * m;
-        let rsp = cs.to_rgb_coeffs(if scale == 1. {
+        let rsp = cs.to_rgb_coeffs(if scale != 0. {
             rgb / scale
         } else {
             RGB::new(0., 0., 0.)

shared/src/utils/options.rs

@@ -78,7 +78,7 @@ impl Default for PBRTOptions {
             image_file: "output.exr",
             mse_reference_image: None,
             mse_reference_output: None,
-            debug_start: Some((Point2i::default(), 0)),
+            debug_start: None,
             display_server: "",
             crop_window: None,
             pixel_bounds: None,

shared/src/utils/ptr.rs

@@ -1,5 +1,7 @@
+use core::hash::{Hash, Hasher};
 use core::marker::PhantomData;
 use core::ops::Index;
+use core::cmp;
 
 #[repr(transparent)]
 #[derive(Debug)]
@@ -20,11 +22,29 @@ impl<T: ?Sized> PartialEq for Ptr<T> {
     }
 }
 
+impl<T: ?Sized> PartialOrd for Ptr<T> {
+    fn partial_cmp(&self, other: &Self) -> Option<cmp::Ordering> {
+        Some(self.ptr.cmp(&other.ptr))
+    }
+}
+
+impl<T: ?Sized> Ord for Ptr<T> {
+    fn cmp(&self, other: &Self) -> cmp::Ordering {
+        self.ptr.cmp(&other.ptr)
+    }
+}
+
 impl<T: ?Sized> Eq for Ptr<T> {}
 
 unsafe impl<T: ?Sized + Send> Send for Ptr<T> {}
 unsafe impl<T: ?Sized + Sync> Sync for Ptr<T> {}
 
+impl<T> Hash for Ptr<T> {
+    fn hash<H: Hasher>(&self, state: &mut H) {
+        self.ptr.hash(state);
+    }
+}
+
 impl<T> Ptr<T> {
     pub const fn null() -> Self {
         Self {
@@ -44,12 +64,6 @@ impl<T> Ptr<T> {
         self.ptr
     }
 
-    #[inline(always)]
-    pub unsafe fn as_ref<'a>(self) -> &'a T {
-        debug_assert!(!self.is_null(), "null Ptr dereference");
-        unsafe { &*self.ptr }
-    }
-
     #[inline(always)]
     pub fn get<'a>(self) -> Option<&'a T> {
         if self.is_null() {
@@ -59,6 +73,15 @@ impl<T> Ptr<T> {
         }
     }
 
+    #[inline(always)]
+    pub fn get_mut<'a>(self) -> Option<&'a mut T> {
+        if self.is_null() {
+            None
+        } else {
+            Some(unsafe { &mut *(self.ptr as *mut T) })
+        }
+    }
+
     #[inline(always)]
     pub unsafe fn at<'a>(self, index: usize) -> &'a T {
         debug_assert!(!self.is_null(), "null Ptr array access");
@@ -98,7 +121,11 @@ impl<T> Ptr<T> {
     #[inline(always)]
     pub fn get_slice<'a>(self, len: usize) -> Option<&'a [T]> {
         if self.is_null() {
-            if len == 0 { Some(&[]) } else { None }
+            if len == 0 {
+                Some(&[])
+            } else {
+                None
+            }
         } else {
             Some(unsafe { core::slice::from_raw_parts(self.ptr, len) })
         }

shared/src/utils/transform.rs

@@ -12,9 +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::{gamma, Float};
+use anyhow::{bail, Context, Result};
 
 #[repr(C)]
 #[derive(Debug, Copy, Clone)]
@@ -113,36 +113,25 @@ impl TransformGeneric<Float> {
         }
     }
 
-    pub fn apply_to_vector(&self, p: Vector3f) -> Vector3f {
+    pub fn apply_to_vector(&self, v: Vector3f) -> Vector3f {
-        let x = p.x();
+        let x = v.x();
-        let y = p.y();
+        let y = v.y();
-        let z = p.z();
+        let z = v.z();
-        let xp = self.m[0][0] * x + self.m[0][1] * y + self.m[0][2] * z;
+        let xv = 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 yv = 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 zv = 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;
+        Vector3f::new(xv, yv, zv)
 
-        if wp == 1. {
-            Vector3f::new(xp, yp, zp)
-        } else {
-            Vector3f::new(xp / wp, yp / wp, zp / wp)
-        }
     }
 
     pub fn apply_to_normal(&self, p: Normal3f) -> Normal3f {
         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 xn = self.m_inv[0][0] * x + self.m_inv[1][1] * y + self.m_inv[2][0] * z;
-        let yp = self.m[1][0] * x + self.m[1][1] * y + self.m[1][2] * z;
+        let yn = self.m_inv[0][1] * x + self.m_inv[1][1] * y + self.m_inv[2][1] * z;
-        let zp = self.m[2][0] * x + self.m[2][1] * y + self.m[2][2] * z;
+        let zn = self.m_inv[0][2] * x + self.m_inv[1][2] * y + self.m_inv[2][2] * z;
-        let wp = self.m[3][0] * x + self.m[3][1] * y + self.m[3][2] * z;
+        Normal3f::new(xn, yn, zn)
-
-        if wp == 1. {
-            Normal3f::new(xp, yp, zp)
-        } else {
-            Normal3f::new(xp / wp, yp / wp, zp / wp)
-        }
     }
 
     pub fn apply_to_bounds(&self, b: Bounds3f) -> Bounds3f {
@@ -155,18 +144,20 @@ impl TransformGeneric<Float> {
     }
 
     pub fn apply_to_ray(&self, r: &Ray, t_max: &mut Option<Float>) -> Ray {
-        let norm_squared = r.d.norm_squared();
+        let pfi = Point3fi::new_from_point(r.o);
-        let mut o = Point3fi::new_from_point(r.o);
+        let mut o = self.apply_to_interval(&pfi);
+        let d = self.apply_to_vector(r.d);
+        let norm_squared = d.norm_squared();
         if norm_squared > 0. {
-            let dt = r.d.abs().dot(o.error()) / norm_squared;
+            let dt = d.abs().dot(o.error()) / norm_squared;
-            let offset = Vector3fi::new_from_vector(r.d * dt);
+            let offset = Vector3fi::new_from_vector(d * dt);
             o = o + offset;
 
             if let Some(t) = t_max.as_mut() {
                 *t -= dt;
             }
         }
-        Ray::new(o.into(), r.d, Some(r.time), r.medium)
+        Ray::new(o.into(), d, Some(r.time), r.medium)
     }
 
     pub fn apply_to_interval(&self, pi: &Point3fi) -> Point3fi {
@@ -493,7 +484,7 @@ impl TransformGeneric<Float> {
     pub fn rotate(sin_theta: Float, cos_theta: Float, axis: impl Into<Vector3f>) -> Self {
         let vec_axis: Vector3f = axis.into();
         let a = vec_axis.normalize();
-        let mut m: SquareMatrix<Float, 4> = SquareMatrix::default();
+        let mut m: SquareMatrix<Float, 4> = SquareMatrix::identity();
         m[0][0] = a.x() * a.x() + (1. - a.x() * a.x()) * cos_theta;
         m[0][1] = a.x() * a.y() * (1. - cos_theta) - a.z() * sin_theta;
         m[0][2] = a.x() * a.z() * (1. - cos_theta) + a.y() * sin_theta;
@@ -527,7 +518,7 @@ impl TransformGeneric<Float> {
         let uu = u.dot(u);
         let vv = v.dot(v);
         let uv = u.dot(v);
-        let mut r: SquareMatrix<Float, 4> = SquareMatrix::default();
+        let mut r: SquareMatrix<Float, 4> = SquareMatrix::identity();
         for i in 0..3 {
             for j in 0..3 {
                 let k = if i == j { 1. } else { 0. };
@@ -770,27 +761,27 @@ impl From<Quaternion> for TransformGeneric<Float> {
         let wy = q.v.y() * q.w;
         let wz = q.v.z() * q.w;
 
-        let mut m = [[0.0; 4]; 4];
+        let mut m_inv = [[0.0; 4]; 4];
 
-        m[0][0] = 1. - 2. * (yy + zz);
+        m_inv[0][0] = 1. - 2. * (yy + zz);
-        m[0][1] = 2. * (xy - wz);
+        m_inv[0][1] = 2. * (xy + wz);
-        m[0][2] = 2. * (xz + wy);
+        m_inv[0][2] = 2. * (xz - wy);
 
-        m[1][0] = 2. * (xy + wz);
+        m_inv[1][0] = 2. * (xy - wz);
-        m[1][1] = 1. - 2. * (xx + zz);
+        m_inv[1][1] = 1. - 2. * (xx + zz);
-        m[1][2] = 2. * (yz - wx);
+        m_inv[1][2] = 2. * (yz + wx);
 
-        m[2][0] = 2. * (xz - wy);
+        m_inv[2][0] = 2. * (xz + wy);
-        m[2][1] = 2. * (yz + wx);
+        m_inv[2][1] = 2. * (yz - wx);
-        m[2][2] = 1. - 2. * (xx + yy);
+        m_inv[2][2] = 1. - 2. * (xx + yy);
 
-        m[3][3] = 1.;
+        m_inv[3][3] = 1.;
 
-        let m_sq = SquareMatrix::new(m);
+        let m_inv_mat = SquareMatrix::new(m_inv);
         // For a pure rotation, the inverse is the transpose.
-        let m_inv_sq = m_sq.transpose();
+        let m = m_inv_mat.transpose();
 
-        TransformGeneric::new(m_sq, m_inv_sq)
+        TransformGeneric::new(m, m_inv_mat)
     }
 }
 
@@ -872,7 +863,22 @@ impl AnimatedTransform {
         let actually_animated = start_transform != end_transform;
 
         if !actually_animated {
-            return Self::default();
+            return Self {
+                start_transform: *start_transform,
+                end_transform: *end_transform,
+                start_time,
+                end_time,
+                actually_animated: false,
+                t: gpu_array_from_fn(|_| Vector3f::default()),
+                r: gpu_array_from_fn(|_| Quaternion::default()),
+                s: gpu_array_from_fn(|_| SquareMatrix::default()),
+                has_rotation: false,
+                c1: gpu_array_from_fn(|_| DerivativeTerm::default()),
+                c2: gpu_array_from_fn(|_| DerivativeTerm::default()),
+                c3: gpu_array_from_fn(|_| DerivativeTerm::default()),
+                c4: gpu_array_from_fn(|_| DerivativeTerm::default()),
+                c5: gpu_array_from_fn(|_| DerivativeTerm::default()),
+            };
         }
 
         let (t0, r_temp, s0) = start_transform.decompose();
@@ -2120,7 +2126,7 @@ pub fn look_at(
     look: impl Into<Point3f>,
     up: impl Into<Point3f>,
 ) -> Result<TransformGeneric<Float>> {
-    let mut world_from_camera: SquareMatrix<Float, 4> = SquareMatrix::default();
+    let mut world_from_camera: SquareMatrix<Float, 4> = SquareMatrix::identity();
     // Initialize fourth column of viewing matrix
     let pos: Point3f = pos.into();
     let look: Point3f = look.into();
@@ -2158,6 +2164,8 @@ 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().context("Failed to inverse viewing matrix")?;
+    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/film.rs

@@ -218,7 +218,7 @@ impl CreatePixelSensor for PixelSensor {
             let target_white = output_colorspace.w;
             xyz_from_sensor_rgb = white_balance(source_white, target_white);
         } else {
-            xyz_from_sensor_rgb = SquareMatrix::<Float, 3>::default();
+            xyz_from_sensor_rgb = SquareMatrix::<Float, 3>::identity();
         }
 
         PixelSensor {

src/core/filter.rs

@@ -35,7 +35,9 @@ impl FilterFactory for Filter {
                 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)))
+                let ptr = arena.alloc(filter);
+                log::warn!("GaussianFilter allocated at {:p}", ptr.as_raw());
+                Ok(Filter::Gaussian(ptr))
             }
             "mitchell" => {
                 let xw = params.get_one_float("xradius", 2.)?;

src/core/interaction.rs

@@ -1,5 +1,4 @@
 use crate::globals::get_options;
-use shared::Ptr;
 use shared::bxdfs::DiffuseBxDF;
 use shared::core::bsdf::BSDF;
 use shared::core::bssrdf::BSSRDF;
@@ -11,6 +10,7 @@ use shared::core::material::{Material, MaterialEvalContext, MaterialTrait};
 use shared::core::sampler::{Sampler, SamplerTrait};
 use shared::core::texture::UniversalTextureEvaluator;
 use shared::spectra::SampledWavelengths;
+use shared::Ptr;
 
 pub trait InteractionGetter {
     fn get_bsdf(
@@ -39,7 +39,9 @@ impl InteractionGetter for SurfaceInteraction {
     ) -> Option<BSDF> {
         self.compute_differentials(r, camera, sampler.samples_per_pixel() as i32);
         let material = {
-            let mut active_mat = unsafe { self.material.as_ref() };
+            let Some(mut active_mat) = self.material.get() else {
+                return None;
+            };
             let tex_eval = UniversalTextureEvaluator;
             while let Material::Mix(mix) = active_mat {
                 let ctx = MaterialEvalContext::from(&*self);
@@ -61,7 +63,7 @@ impl InteractionGetter for SurfaceInteraction {
         if get_options().force_diffuse {
             let rho = bsdf.rho_wo(self.common.wo, &[sampler.get1d()], &[sampler.get2d()]);
             let diff_bxdf = BxDF::Diffuse(DiffuseBxDF::new(rho));
-            bsdf = BSDF::new(self.shading.n, self.shading.dpdu, Ptr::from(&diff_bxdf));
+            bsdf = BSDF::new(self.shading.n, self.shading.dpdu, diff_bxdf);
         }
         Some(bsdf)
     }
@@ -72,7 +74,9 @@ impl InteractionGetter for SurfaceInteraction {
         lambda: &SampledWavelengths,
         _camera: &Camera,
     ) -> Option<BSSRDF> {
-        let mut active_mat = unsafe { self.material.as_ref() };
+        let Some(mut active_mat) = self.material.get() else {
+            return None;
+        };
         let tex_eval = UniversalTextureEvaluator;
         while let Material::Mix(mix) = active_mat {
             let ctx = MaterialEvalContext::from(self);

src/core/render.rs

@@ -13,29 +13,43 @@ use shared::spectra::{SampledWavelengths, LAMBDA_MAX, LAMBDA_MIN};
 use shared::Float;
 use std::sync::Arc;
 
-fn render_scene(scene: &BasicScene, arena: &Arena) -> Result<()> {
+pub fn render_scene(scene: &BasicScene, arena: &Arena) -> Result<()> {
     let media = scene.create_media();
     let textures = scene.create_textures(arena);
-    let (lights, _) = scene.create_lights(&textures, arena);
     let (named_materials, materials) = scene.create_materials(&textures, arena)?;
+    for (i, m) in materials.iter().enumerate() {
+        eprintln!("materials[{}]: {:?}", i, std::mem::discriminant(m));
+    }
+    let lights = scene.create_lights(&textures, arena);
+
+    let have_scattering = {
+        let shapes = scene.shapes.lock();
+        let animated = scene.animated_shapes.lock();
+        shapes
+            .iter()
+            .any(|sh| !sh.inside_medium.is_empty() || !sh.outside_medium.is_empty())
+            || animated
+                .iter()
+                .any(|sh| !sh.inside_medium.is_empty() || !sh.outside_medium.is_empty())
+    };
+
     let (aggregate, area_lights) =
-        scene.create_aggregate(&textures, &named_materials, &materials, arena);
+        scene.create_aggregate(&textures, &named_materials, &materials, &media, arena);
+
+    let mut all_lights = lights;
+    all_lights.extend(area_lights);
+
     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.clone(), sampler.clone(), aggregate.clone(), lights, arena);
+    let integrator = scene.create_integrator(
-    let mut have_scattering = false;
+        camera.clone(),
-    for sh in scene.shapes.lock().iter() {
+        sampler.clone(),
-        if !sh.inside_medium.is_empty() || !sh.outside_medium.is_empty() {
+        aggregate.clone(),
-            have_scattering = true;
+        all_lights,
-        }
+        arena,
-    }
+    );
-    for sh in scene.animated_shapes.lock().iter() {
-        if !sh.inside_medium.is_empty() || !sh.outside_medium.is_empty() {
-            have_scattering = true;
-        }
-    }
 
     if get_options().pixel_material.is_some() {
         let lambda =
@@ -75,20 +89,13 @@ fn render_scene(scene: &BasicScene, arena: &Arena) -> Result<()> {
                     );
                     log::debug!("Distance from camera: {}\n", intr.p().distance(cr.ray.o));
 
-                    let mut is_named = false;
                     for (name, mtl) in &named_materials {
-                        if *mtl == unsafe { *intr.material.as_ref() } {
+                        if *mtl == *intr.material.get().unwrap() {
                             log::debug!("Named material: {}\n\n", name);
-                            is_named = true;
                             break;
                         }
                     }
 
-                    // if !is_named {
-                    //     log::warn!("{}\n\n", intr.material.as_ref().to_str());
-                    // }
-                    //
-
                     depth += 1;
                     ray = intr.spawn_ray(ray.d);
                 }

src/core/scene/builder.rs

@@ -61,7 +61,7 @@ struct PendingAreaLight {
     loc: FileLoc,
 }
 
-#[derive(Default, Debug, Clone)]
+#[derive(Debug, Clone)]
 struct GraphicsState {
     pub current_inside_medium: String,
     pub current_outside_medium: String,
@@ -81,6 +81,29 @@ struct GraphicsState {
     pub transform_end_time: Float,
 }
 
+impl Default for GraphicsState {
+    fn default() -> Self {
+        Self {
+            color_space: Some(crate::spectra::default_colorspace_arc()),
+            active_transform_bits: BasicSceneBuilder::ALL_TRANSFORM_BITS,
+            current_inside_medium: String::new(),
+            current_outside_medium: String::new(),
+            current_material_name: String::new(),
+            current_material_index: None,
+            pending_area_light: None,
+            shape_attributes: Vec::new(),
+            light_attributes: Vec::new(),
+            material_attributes: Vec::new(),
+            medium_attributes: Vec::new(),
+            texture_attributes: Vec::new(),
+            reverse_orientation: false,
+            ctm: TransformSet::default(),
+            transform_start_time: 0.0,
+            transform_end_time: 1.0,
+        }
+    }
+}
+
 #[derive(PartialEq, Eq)]
 enum BlockState {
     OptionsBlock,
@@ -115,6 +138,19 @@ impl BasicSceneBuilder {
     pub const END_TRANSFORM_BITS: u32 = 1 << 1;
     pub const ALL_TRANSFORM_BITS: u32 = (1 << MAX_TRANSFORMS) - 1;
 
+    fn render_from_object_at(&self, index: usize) -> Transform {
+        self.render_from_world * self.graphics_state.ctm[index]
+    }
+
+    fn render_from_object(&self) -> AnimatedTransform {
+        AnimatedTransform::new(
+            &self.render_from_object_at(0),
+            self.graphics_state.transform_start_time,
+            &self.render_from_object_at(1),
+            self.graphics_state.transform_end_time,
+        )
+    }
+
     pub fn new(scene: Arc<BasicScene>) -> Self {
         Self {
             scene,
@@ -226,7 +262,7 @@ impl ParserTarget for BasicSceneBuilder {
         let stdcs = get_colorspace_device();
         let _ = match stdcs.get_named(name) {
             Some(cs) => {
-                self.graphics_state.color_space = unsafe { Some(Arc::new(*cs.as_ref())) };
+                self.graphics_state.color_space = Some(Arc::new(*cs.get().unwrap()));
             }
             None => {
                 eprintln!("Error: Color space '{}' unknown at {}", name, loc);
@@ -343,7 +379,6 @@ impl ParserTarget for BasicSceneBuilder {
         self.verify_options("Camera", &loc)?;
 
         let camera_from_world = self.graphics_state.ctm;
-
         let world_from_camera = camera_from_world.inverse();
 
         self.named_coordinate_systems
@@ -359,7 +394,8 @@ impl ParserTarget for BasicSceneBuilder {
         let rendering_space = RenderingCoordinateSystem::CameraWorld;
         let camera_transform =
             CameraTransform::from_world(animated_world_from_cam, rendering_space);
-        self.render_from_world = camera_from_world.t[0];
+
+        self.render_from_world = camera_transform.render_from_world();
 
         let parameters = self.make_params(params, &loc)?;
 
@@ -526,7 +562,7 @@ impl ParserTarget for BasicSceneBuilder {
         Ok(())
     }
 
-    fn world_begin(&mut self, loc: FileLoc, arena: Arc<Arena>) -> Result<(), ParserError> {
+    fn world_begin(&mut self, loc: FileLoc, arena: &Arena) -> Result<(), ParserError> {
         self.verify_options("WorldBegin", &loc)?;
         self.current_block = BlockState::WorldBlock;
         for i in 0..MAX_TRANSFORMS {
@@ -684,7 +720,7 @@ impl ParserTarget for BasicSceneBuilder {
         };
         let entity = TextureSceneEntity {
             base,
-            render_from_object: AnimatedTransform::from_transform(&self.graphics_state.ctm[0]),
+            render_from_object: self.render_from_object(),
         };
 
         if type_name == "float" {
@@ -710,9 +746,12 @@ impl ParserTarget for BasicSceneBuilder {
             loc,
             parameters: ParameterDictionary::new(params.clone(), None).unwrap(),
         };
-        self.graphics_state.current_material_name = self.scene.add_material(entity).to_string();
+        let idx = self.scene.add_material(entity);
+        self.graphics_state.current_material_index = Some(idx);
+        self.graphics_state.current_material_name = String::new();
         Ok(())
     }
+
     fn make_named_material(
         &mut self,
         _name: &str,
@@ -739,12 +778,7 @@ impl ParserTarget for BasicSceneBuilder {
             self.graphics_state.color_space.clone(),
         )?;
 
-        let render_from_light = AnimatedTransform::new(
+        let render_from_light = self.render_from_object();
-            &self.graphics_state.ctm.t[0],
-            self.graphics_state.transform_start_time,
-            &self.graphics_state.ctm.t[1],
-            self.graphics_state.transform_end_time,
-        );
 
         let entity = LightSceneEntity {
             transformed_base: TransformedSceneEntity {
@@ -791,7 +825,7 @@ impl ParserTarget for BasicSceneBuilder {
             self.graphics_state.color_space.clone(),
         )?;
 
-        let render_from_object = self.graphics_state.ctm[0];
+        let render_from_object = self.render_from_object_at(0);
         let object_from_render = render_from_object.inverse();
 
         let light_index = if let Some(ref al) = self.graphics_state.pending_area_light {
@@ -841,6 +875,7 @@ impl ParserTarget for BasicSceneBuilder {
 
         Ok(())
     }
+
     fn object_begin(&mut self, _name: &str, _loc: FileLoc) -> Result<(), ParserError> {
         Ok(())
     }

src/core/scene/mod.rs

@@ -5,7 +5,5 @@ pub mod state;
 
 pub use builder::BasicSceneBuilder;
 pub use entities::*;
-pub use scene::{BasicScene, SceneLookup};
+pub use scene::BasicScene;
 pub use state::*;
-
-

src/core/shape.rs

@@ -1,12 +1,13 @@
 use crate::core::texture::FloatTexture;
 use crate::shapes::*;
+use crate::utils::loop_subdivide;
 use crate::utils::resolve_filename;
 use crate::{Arena, FileLoc, ParameterDictionary};
 use anyhow::{anyhow, bail, Result};
 use parking_lot::Mutex;
 use shared::core::shape::*;
 use shared::shapes::*;
-use shared::{Transform, Ptr};
+use shared::{Ptr, Transform};
 use std::collections::HashMap;
 use std::sync::Arc;
 
@@ -130,6 +131,79 @@ impl ShapeFactory for Shape {
 
                 Ok(shapes)
             }
+            "loopsubdiv" => {
+                let n_levels = parameters.get_one_int("levels", 3)? as usize;
+
+                let Ok(vertex_indices_i32) = parameters.get_int_array("indices") else {
+                    bail!("Vertex indices \"indices\" not provided for \"LoopSubdiv\" shape.");
+                };
+                let Ok(p) = parameters.get_point3f_array("P") else {
+                    bail!("Vertex positions \"P\" not provided for LoopSubdiv shape.");
+                };
+
+                if vertex_indices_i32.is_empty() {
+                    bail!("Vertex indices \"indices\" must be provided with loopsubdiv.");
+                }
+                if vertex_indices_i32.len() % 3 != 0 {
+                    bail!(
+                        "Number of vertex indices {} not a multiple of 3.",
+                        vertex_indices_i32.len()
+                    );
+                }
+                if p.is_empty() {
+                    bail!("Vertex positions \"P\" must be provided with loopsubdiv.");
+                }
+
+                let vertex_indices: Vec<usize> = vertex_indices_i32
+                    .iter()
+                    .map(|&i| {
+                        if i < 0 {
+                            bail!("Negative vertex index {} in loopsubdiv.", i);
+                        }
+                        let ui = i as usize;
+                        if ui >= p.len() {
+                            bail!(
+                                "Vertex index {} out of bounds ({} \"P\" values given).",
+                                i,
+                                p.len()
+                            );
+                        }
+                        Ok(ui)
+                    })
+                    .collect::<Result<Vec<_>>>()?;
+
+                let scheme = parameters.get_one_string("scheme", "loop")?;
+                if scheme != "loop" {
+                    bail!("Unsupported subdivision scheme \"{}\".", scheme);
+                }
+
+                let mesh = loop_subdivide(
+                    arena,
+                    &render_from_object,
+                    reverse_orientation,
+                    n_levels,
+                    &vertex_indices,
+                    &p,
+                );
+
+                let host_arc = Arc::new(mesh);
+                let mut global_store = ALL_TRIANGLE_MESHES.lock();
+                global_store.push(host_arc.clone());
+                drop(global_store);
+
+                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| {
+                        arena.alloc(Shape::Triangle(TriangleShape {
+                            mesh: mesh_ptr,
+                            tri_index: i as i32,
+                        }))
+                    })
+                    .collect();
+
+                Ok(shapes)
+            }
             _ => Err(anyhow!("Unknown shape name")),
         }
     }

src/integrators/path.rs

@@ -1,10 +1,10 @@
-use super::RayIntegratorTrait;
 use super::base::IntegratorBase;
 use super::constants::*;
 use super::state::PathState;
-use crate::Arena;
+use super::RayIntegratorTrait;
 use crate::core::interaction::InteractionGetter;
-use shared::core::bsdf::{BSDF, BSDFSample};
+use crate::Arena;
+use shared::core::bsdf::{BSDFSample, BSDF};
 use shared::core::bxdf::{BxDFFlags, FArgs, TransportMode};
 use shared::core::camera::Camera;
 use shared::core::film::VisibleSurface;
@@ -308,6 +308,8 @@ impl RayIntegratorTrait for PathIntegrator {
         sampler: &mut Sampler,
         arena: &Arena,
     ) {
+        use shared::core::primitive::PrimitiveTrait;
+        eprintln!("BVH bounds: {:?}", &self.base.aggregate.bounds());
         crate::integrators::pipeline::evaluate_pixel_sample(
             self,
             self.camera.as_ref(),

src/integrators/pipeline.rs

@@ -101,7 +101,11 @@ pub fn render<T>(
         return;
     }
 
-    let pixel_bounds = camera.get_film().pixel_bounds();
+    let sample_bounds = camera.get_film().sample_bounds();
+    let pixel_bounds = Bounds2i::from_points(
+    Point2i::new(sample_bounds.p_min.x().floor() as i32, sample_bounds.p_min.y().floor() as i32),
+    Point2i::new(sample_bounds.p_max.x().ceil() as i32, sample_bounds.p_max.y().ceil() as i32),
+);
     println!(
         "pixel_bounds: {:?}, area: {}",
         pixel_bounds,
@@ -165,15 +169,14 @@ pub fn render<T>(
         tiles.par_iter().for_each(|tile_bounds| {
             let mut sampler = sampler_prototype.clone();
 
-            for p_pixel in tile_bounds {
+            for p_pixel in tile_bounds.pixels() {
                 for sample_index in wave_start..wave_end {
-                    sampler.start_pixel_sample(*p_pixel, sample_index, None);
+                    sampler.start_pixel_sample(p_pixel, sample_index, None);
-                    println!("Evaluating pixel {:?} sample {}", p_pixel, sample_index);
                     evaluate_pixel_sample(
                         integrator,
                         camera,
                         &mut sampler,
-                        *p_pixel,
+                        p_pixel,
                         sample_index.try_into().unwrap(),
                         &arena,
                     );
@@ -252,8 +255,8 @@ pub fn evaluate_pixel_sample<T: RayIntegratorTrait>(
     let filter = film.get_filter();
     let camera_sample = get_camera_sample(sampler, pixel, filter);
     if let Some(mut camera_ray) = camera.generate_ray_differential(camera_sample, &lambda) {
-        debug_assert!(camera_ray.ray.d.norm() > 0.999);
+        // debug_assert!(camera_ray.ray.d.norm() > 0.999);
-        debug_assert!(camera_ray.ray.d.norm() < 1.001);
+        // debug_assert!(camera_ray.ray.d.norm() < 1.001);
         let ray_diff_scale = (sampler.samples_per_pixel() as Float).sqrt().max(0.125);
         if get_options().disable_pixel_jitter {
             camera_ray.ray.scale_differentials(ray_diff_scale);
@@ -275,10 +278,10 @@ pub fn evaluate_pixel_sample<T: RayIntegratorTrait>(
         }
 
         if pixel.x() == 352 && pixel.y() == 352 && sample_index == 0 {
-            println!("Center pixel: L = {:?}", l);
+            eprintln!("Center pixel: L = {:?}", l);
-            println!("  ray origin: {:?}", camera_ray.ray.o);
+            eprintln!("  ray origin: {:?}", camera_ray.ray.o);
-            println!("  ray dir: {:?}", camera_ray.ray.d);
+            eprintln!("  ray dir: {:?}", camera_ray.ray.d);
-            println!("  camera_sample.p_film: {:?}", camera_sample.p_film);
+            eprintln!("  camera_sample.p_film: {:?}", camera_sample.p_film);
         }
 
         film.add_sample(

src/lights/diffuse.rs

@@ -28,7 +28,7 @@ pub fn create(
     colorspace: Option<&RGBColorSpace>,
     arena: &Arena,
 ) -> Result<Light> {
-    let mut l = params.get_one_spectrum("l", None, SpectrumType::Illuminant);
+    let mut l = params.get_one_spectrum("L", None, SpectrumType::Illuminant);
     let default_cs = crate::spectra::default_colorspace();
     let cs = colorspace.unwrap_or(&default_cs);
     let illum_spec = Spectrum::Dense(cs.illuminant);
@@ -101,7 +101,7 @@ pub fn create(
 
     // Upload alpha texture to GPU and check for null texture
     let alpha_ptr = arena.upload(alpha);
-    let light_type = match unsafe { alpha_ptr.as_ref() } {
+    let light_type = match alpha_ptr.get().unwrap() {
         GPUFloatTexture::Constant(t) if t.evaluate(&TextureEvalContext::default()) == 0.0 => {
             LightType::DeltaPosition
         }

src/materials/conductor.rs

@@ -1,20 +1,68 @@
 use crate::core::image::HostImage;
 use crate::core::material::CreateMaterial;
-// use crate::core::scattering::TrowbridgeReitzDistribution;
+use crate::core::texture::SpectrumTexture;
+use crate::spectra::data::get_named_spectrum;
 use crate::utils::TextureParameterDictionary;
-use anyhow::Result;
+use crate::{Arena, ArenaUpload, FileLoc};
+use anyhow::{bail, Result};
+use shared::core::texture::SpectrumType;
 use shared::core::material::Material;
 use shared::materials::ConductorMaterial;
+use shared::textures::SpectrumConstantTexture;
+use std::collections::HashMap;
 use std::sync::Arc;
 
 impl CreateMaterial for ConductorMaterial {
     fn create(
-        _parameters: &TextureParameterDictionary,
+        parameters: &TextureParameterDictionary,
-        _normal_map: Option<Arc<HostImage>>,
+        normal_map: Option<Arc<HostImage>>,
-        _named_materials: &std::collections::HashMap<String, shared::core::material::Material>,
+        _named_materials: &HashMap<String, Material>,
-        _loc: &crate::utils::FileLoc,
+        loc: &FileLoc,
-        _arena: &crate::Arena,
+        arena: &Arena,
     ) -> Result<Material> {
-        todo!()
+        let eta = parameters.get_spectrum_texture_or_null("eta", SpectrumType::Unbounded);
+        let k = parameters.get_spectrum_texture_or_null("k", SpectrumType::Unbounded);
+        let reflectance =
+            parameters.get_spectrum_texture_or_null("reflectance", SpectrumType::Albedo);
+        let (eta, k) = if reflectance.is_some() {
+            if eta.is_some() || k.is_some() {
+                bail!(
+                    "{}: For the coated conductor material, both \"reflectance\" \
+         and \"eta\"/\"k\" cannot be provided simultaneously.",
+                    loc
+                );
+            }
+            (None, None)
+        } else {
+            let eta = eta.unwrap_or_else(|| {
+                let s = get_named_spectrum("metal-Cu-eta").expect("Missing copper spectrum");
+                Arc::new(SpectrumTexture::Constant(SpectrumConstantTexture::new(s)))
+            });
+            let k = k.unwrap_or_else(|| {
+                let s = get_named_spectrum("metal-Cu-k").expect("Missing copper spectrum");
+                Arc::new(SpectrumTexture::Constant(SpectrumConstantTexture::new(s)))
+            });
+            (Some(eta), Some(k))
+        };
+
+        let u_roughness = parameters.get_float_texture("roughness", 0.)?;
+        let v_roughness = parameters.get_float_texture("roughness", 0.)?;
+        let displacement = parameters.get_float_texture_or_null("displacement")?;
+        let remap_roughness = parameters.get_one_bool("remaproughness", true)?;
+        eprintln!("ConductorMaterial::create eta={:?} k={:?} reflectance={:?}", eta.is_some(), k.is_some(), reflectance.is_some());
+
+        let material = Self::new(
+            arena.upload(normal_map),
+            arena.upload(reflectance),
+            arena.upload(eta),
+            arena.upload(k),
+            arena.upload(u_roughness),
+            arena.upload(v_roughness),
+            arena.upload(displacement),
+            remap_roughness,
+        );
+
+        arena.alloc(material);
+        Ok(Material::Conductor(material))
     }
 }

src/samplers/halton.rs

@@ -1,9 +1,11 @@
 use super::*;
 use crate::globals::get_options;
-use shared::utils::math::compute_radical_inverse_permutations;
+use crate::Arena;
-use anyhow::{Result, anyhow};
+use anyhow::{anyhow, Result};
 use shared::core::geometry::Point2i;
-use shared::core::sampler::{HaltonSampler, MAX_HALTON_RESOLUTION, RandomizeStrategy};
+use shared::core::sampler::{HaltonSampler, RandomizeStrategy, MAX_HALTON_RESOLUTION};
+use shared::utils::math::compute_radical_inverse_permutations;
+use shared::{gbox, Ptr};
 
 pub trait CreateHaltonSampler {
     fn new(
@@ -11,6 +13,7 @@ pub trait CreateHaltonSampler {
         full_res: Point2i,
         randomize: RandomizeStrategy,
         seed: u64,
+        arena: &Arena,
     ) -> Self;
 }
 
@@ -20,8 +23,11 @@ impl CreateHaltonSampler for HaltonSampler {
         full_res: Point2i,
         randomize: RandomizeStrategy,
         seed: u64,
+        arena: &Arena,
     ) -> Self {
         let digit_permutations = compute_radical_inverse_permutations(seed);
+        let leaked = Box::leak(gbox(digit_permutations));
+        let perm_ptr = Ptr::from(leaked.as_slice());
         let mut base_scales = [0u64; 2];
         let mut base_exponents = [0u64; 2];
         let bases = [2, 3];
@@ -46,14 +52,14 @@ impl CreateHaltonSampler for HaltonSampler {
         let mut mult_inverse = [0u64; 2];
 
         mult_inverse[0] =
-            Self::multiplicative_inverse(base_scales[0] as i64, base_scales[0] as i64);
+            Self::multiplicative_inverse(base_scales[1] as i64, base_scales[0] as i64);
         mult_inverse[1] =
-            Self::multiplicative_inverse(base_scales[1] as i64, base_scales[1] as i64);
+            Self::multiplicative_inverse(base_scales[0] as i64, base_scales[1] as i64);
 
         Self {
             samples_per_pixel,
             randomize,
-            digit_permutations: digit_permutations.as_ptr().into(),
+            digit_permutations: perm_ptr,
             base_scales,
             base_exponents,
             mult_inverse,
@@ -68,7 +74,7 @@ impl CreateSampler for HaltonSampler {
         params: &ParameterDictionary,
         full_res: Point2i,
         loc: &FileLoc,
-        _arena: &Arena,
+        arena: &Arena,
     ) -> Result<Sampler> {
         let options = get_options();
         let nsamp = if let Some(n) = options.quick_render.then_some(1).or(options.pixel_samples) {
@@ -94,8 +100,7 @@ impl CreateSampler for HaltonSampler {
             }
         };
 
-        let sampler = HaltonSampler::new(nsamp, full_res, s, seed as u64);
+        let sampler = HaltonSampler::new(nsamp, full_res, s, seed as u64, arena);
-        // arena.alloc(sampler);
         Ok(Sampler::Halton(sampler))
     }
 }

src/textures/scaled.rs

@@ -3,8 +3,9 @@ use crate::core::texture::{FloatTextureTrait, SpectrumTextureTrait};
 use crate::utils::{FileLoc, TextureParameterDictionary};
 use crate::Arena;
 use anyhow::Result;
+use shared::core::spectrum::Spectrum;
 use shared::core::texture::{SpectrumType, TextureEvalContext};
-use shared::spectra::{SampledSpectrum, SampledWavelengths};
+use shared::spectra::{SampledSpectrum, SampledWavelengths, ConstantSpectrum};
 use shared::utils::Transform;
 use shared::Float;
 use std::sync::Arc;
@@ -70,11 +71,32 @@ pub struct SpectrumScaledTexture {
 impl CreateSpectrumTexture for SpectrumScaledTexture {
     fn create(
         _render_from_texture: Transform,
-        _parameters: TextureParameterDictionary,
+        parameters: TextureParameterDictionary,
-        _spectrum_type: SpectrumType,
+        spectrum_type: SpectrumType,
         _loc: FileLoc,
     ) -> Result<SpectrumTexture> {
-        todo!()
+        let one = Spectrum::Constant(ConstantSpectrum::new(1.0));
+        let tex = parameters
+            .get_spectrum_texture("tex", Some(one), spectrum_type)
+            .ok_or_else(|| anyhow::anyhow!("{}: missing \"tex\" parameter", _loc))?;
+        let scale = parameters.get_float_texture("scale", 1.0)?;
+
+        if let FloatTexture::Constant(ref cscale) = *scale {
+            let cs = cscale.value;
+            if cs == 1.0 {
+                return Ok((*tex).clone());
+            }
+            if let SpectrumTexture::Image(ref image) = *tex {
+                let mut image_copy = image.clone();
+                image_copy.base.multiply_scale(cs);
+                return Ok(SpectrumTexture::Image(image_copy));
+            }
+        }
+
+        Ok(SpectrumTexture::Scaled(SpectrumScaledTexture {
+            tex,
+            scale,
+        }))
     }
 }
 

src/utils/loopsubdiv.rs

@@ -0,0 +1,473 @@
+use crate::Arena;
+use shared::core::geometry::{Normal3f, Point3f, Vector3f};
+use shared::shapes::TriangleMesh;
+use shared::{Float, Transform};
+use std::collections::HashMap;
+
+const NULL: u32 = u32::MAX;
+
+#[repr(C)]
+struct SDVertex {
+    p: Point3f,
+    start_face: u32,
+    child: u32,
+    regular: bool,
+    boundary: bool,
+}
+
+#[repr(C)]
+struct SDFace {
+    v: [u32; 3],
+    f: [u32; 3],
+    children: [u32; 4],
+}
+
+impl SDFace {
+    fn vnum(&self, vert: u32) -> usize {
+        for i in 0..3 {
+            if self.v[i] == vert {
+                return i;
+            }
+        }
+        panic!("SDFace::vnum: vertex {} not in face", vert);
+    }
+
+    fn next_face(&self, vert: u32) -> u32 { self.f[self.vnum(vert)] }
+    fn prev_face(&self, vert: u32) -> u32 { self.f[prev(self.vnum(vert))] }
+    fn next_vert(&self, vert: u32) -> u32 { self.v[next(self.vnum(vert))] }
+    fn prev_vert(&self, vert: u32) -> u32 { self.v[prev(self.vnum(vert))] }
+
+    fn other_vert(&self, v0: u32, v1: u32) -> u32 {
+        for i in 0..3 {
+            if self.v[i] != v0 && self.v[i] != v1 {
+                return self.v[i];
+            }
+        }
+        panic!("SDFace::other_vert: face doesn't contain both vertices");
+    }
+}
+
+#[inline(always)]
+const fn next(i: usize) -> usize { (i + 1) % 3 }
+
+#[inline(always)]
+const fn prev(i: usize) -> usize { (i + 2) % 3 }
+
+#[derive(Hash, Eq, PartialEq, Copy, Clone)]
+struct EdgeKey(u32, u32);
+
+impl EdgeKey {
+    fn new(a: u32, b: u32) -> Self {
+        if a <= b { Self(a, b) } else { Self(b, a) }
+    }
+}
+
+fn beta(valence: usize) -> Float {
+    if valence == 3 { 3.0 / 16.0 } else { 3.0 / (8.0 * valence as Float) }
+}
+
+fn loop_gamma(valence: usize) -> Float {
+    1.0 / (valence as Float + 3.0 / (8.0 * beta(valence)))
+}
+
+fn valence(verts: &[SDVertex], faces: &[SDFace], vi: u32) -> usize {
+    let v = &verts[vi as usize];
+    let start = v.start_face;
+    debug_assert_ne!(start, NULL);
+
+    if !v.boundary {
+        let mut nf = 1usize;
+        let mut fi = faces[start as usize].next_face(vi);
+        while fi != start {
+            assert_ne!(fi, NULL, "interior vertex with broken face loop");
+            nf += 1;
+            fi = faces[fi as usize].next_face(vi);
+        }
+        nf
+    } else {
+        let mut nf = 1usize;
+        let mut fi = faces[start as usize].next_face(vi);
+        while fi != NULL {
+            nf += 1;
+            fi = faces[fi as usize].next_face(vi);
+        }
+        fi = faces[start as usize].prev_face(vi);
+        while fi != NULL {
+            nf += 1;
+            fi = faces[fi as usize].prev_face(vi);
+        }
+        nf + 1
+    }
+}
+
+fn one_ring(verts: &[SDVertex], faces: &[SDFace], vi: u32) -> Vec<Point3f> {
+    let v = &verts[vi as usize];
+    let mut ring = Vec::new();
+
+    if !v.boundary {
+        let start = v.start_face;
+        let mut fi = start;
+        loop {
+            ring.push(verts[faces[fi as usize].next_vert(vi) as usize].p);
+            fi = faces[fi as usize].next_face(vi);
+            if fi == start { break; }
+        }
+    } else {
+        // Walk forward to the last face in the fan
+        let mut fi = v.start_face;
+        loop {
+            let nf = faces[fi as usize].next_face(vi);
+            if nf == NULL { break; }
+            fi = nf;
+        }
+        ring.push(verts[faces[fi as usize].next_vert(vi) as usize].p);
+        loop {
+            ring.push(verts[faces[fi as usize].prev_vert(vi) as usize].p);
+            let pf = faces[fi as usize].prev_face(vi);
+            if pf == NULL { break; }
+            fi = pf;
+        }
+    }
+    ring
+}
+
+fn weight_one_ring(verts: &[SDVertex], faces: &[SDFace], vi: u32, beta: Float) -> Point3f {
+    let ring = one_ring(verts, faces, vi);
+    let val = ring.len();
+    let mut acc = (1.0 - val as Float * beta) * Vector3f::from(verts[vi as usize].p);
+    for r in &ring {
+        acc += beta * Vector3f::from(*r);
+    }
+    Point3f::from(acc)
+}
+
+fn weight_boundary(verts: &[SDVertex], faces: &[SDFace], vi: u32, beta: Float) -> Point3f {
+    let ring = one_ring(verts, faces, vi);
+    let val = ring.len();
+    let acc = (1.0 - 2.0 * beta) * Vector3f::from(verts[vi as usize].p)
+        + beta * Vector3f::from(ring[0])
+        + beta * Vector3f::from(ring[val - 1]);
+    Point3f::from(acc)
+}
+
+pub fn loop_subdivide(
+    _arena: &Arena,
+    render_from_object: &Transform,
+    reverse_orientation: bool,
+    n_levels: usize,
+    vertex_indices: &[usize],
+    points: &[Point3f],
+) -> TriangleMesh {
+    let n_faces_init = vertex_indices.len() / 3;
+
+    let mut verts: Vec<SDVertex> = points
+        .iter()
+        .map(|&p| SDVertex {
+            p,
+            start_face: NULL,
+            child: NULL,
+            regular: false,
+            boundary: false,
+        })
+        .collect();
+
+    let mut faces: Vec<SDFace> = (0..n_faces_init)
+        .map(|i| {
+            let base = i * 3;
+            SDFace {
+                v: [
+                    vertex_indices[base] as u32,
+                    vertex_indices[base + 1] as u32,
+                    vertex_indices[base + 2] as u32,
+                ],
+                f: [NULL; 3],
+                children: [NULL; 4],
+            }
+        })
+        .collect();
+
+    for fi in 0..n_faces_init {
+        for j in 0..3 {
+            verts[faces[fi].v[j] as usize].start_face = fi as u32;
+        }
+    }
+
+    {
+        let mut edge_map: HashMap<EdgeKey, (u32, usize)> =
+            HashMap::with_capacity(n_faces_init * 3 / 2);
+
+        for fi in 0..n_faces_init {
+            for edge_num in 0..3 {
+                let v0 = faces[fi].v[edge_num];
+                let v1 = faces[fi].v[next(edge_num)];
+                let key = EdgeKey::new(v0, v1);
+
+                if let Some(&(other_fi, other_edge)) = edge_map.get(&key) {
+                    faces[other_fi as usize].f[other_edge] = fi as u32;
+                    faces[fi].f[edge_num] = other_fi;
+                    edge_map.remove(&key);
+                } else {
+                    edge_map.insert(key, (fi as u32, edge_num));
+                }
+            }
+        }
+    }
+
+    // Classify vertices as boundary/interior and regular/extraordinary
+    for vi in 0..verts.len() {
+        let start = verts[vi].start_face;
+        let mut fi = start;
+        let is_boundary = loop {
+            let nf = faces[fi as usize].next_face(vi as u32);
+            if nf == NULL { break true; }
+            fi = nf;
+            if fi == start { break false; }
+        };
+        verts[vi].boundary = is_boundary;
+        let val = valence(&verts, &faces, vi as u32);
+        verts[vi].regular = if !is_boundary { val == 6 } else { val == 4 };
+    }
+
+    let mut live_verts: Vec<u32> = (0..verts.len() as u32).collect();
+    let mut live_faces: Vec<u32> = (0..faces.len() as u32).collect();
+
+    for _level in 0..n_levels {
+        let old_verts = live_verts.clone();
+        let old_faces = live_faces.clone();
+        let mut new_verts: Vec<u32> = Vec::with_capacity(old_verts.len() + old_faces.len() * 3);
+        let mut new_faces: Vec<u32> = Vec::with_capacity(old_faces.len() * 4);
+
+        // Allocate child (even) vertices — one per existing vertex
+        for &vi in &old_verts {
+            let child_id = verts.len() as u32;
+            verts.push(SDVertex {
+                p: Point3f::default(),
+                start_face: NULL,
+                child: NULL,
+                regular: verts[vi as usize].regular,
+                boundary: verts[vi as usize].boundary,
+            });
+            verts[vi as usize].child = child_id;
+            new_verts.push(child_id);
+        }
+
+        // Allocate 4 child faces per existing face
+        for &fi in &old_faces {
+            let base = faces.len() as u32;
+            for _ in 0..4 {
+                faces.push(SDFace {
+                    v: [NULL; 3],
+                    f: [NULL; 3],
+                    children: [NULL; 4],
+                });
+            }
+            faces[fi as usize].children = [base, base + 1, base + 2, base + 3];
+            for k in 0..4u32 {
+                new_faces.push(base + k);
+            }
+        }
+
+        // Reposition even vertices using the Loop stencil
+        for &vi in &old_verts {
+            let child = verts[vi as usize].child;
+            let new_p = if !verts[vi as usize].boundary {
+                if verts[vi as usize].regular {
+                    weight_one_ring(&verts, &faces, vi, 1.0 / 16.0)
+                } else {
+                    let val = valence(&verts, &faces, vi);
+                    weight_one_ring(&verts, &faces, vi, beta(val))
+                }
+            } else {
+                weight_boundary(&verts, &faces, vi, 1.0 / 8.0)
+            };
+            verts[child as usize].p = new_p;
+        }
+
+        // Create odd (edge) vertices
+        let mut edge_verts: HashMap<EdgeKey, u32> = HashMap::new();
+
+        for &fi in &old_faces {
+            for k in 0..3 {
+                let v0 = faces[fi as usize].v[k];
+                let v1 = faces[fi as usize].v[next(k)];
+                let key = EdgeKey::new(v0, v1);
+
+                if edge_verts.contains_key(&key) { continue; }
+
+                let is_boundary = faces[fi as usize].f[k] == NULL;
+                let new_vi = verts.len() as u32;
+
+                let new_p = if is_boundary {
+                    Point3f::from(
+                        0.5 * Vector3f::from(verts[v0 as usize].p)
+                            + 0.5 * Vector3f::from(verts[v1 as usize].p),
+                    )
+                } else {
+                    let neighbor_fi = faces[fi as usize].f[k];
+                    let other_a = faces[fi as usize].other_vert(v0, v1);
+                    let other_b = faces[neighbor_fi as usize].other_vert(v0, v1);
+                    Point3f::from(
+                        (3.0 / 8.0) * Vector3f::from(verts[v0 as usize].p)
+                            + (3.0 / 8.0) * Vector3f::from(verts[v1 as usize].p)
+                            + (1.0 / 8.0) * Vector3f::from(verts[other_a as usize].p)
+                            + (1.0 / 8.0) * Vector3f::from(verts[other_b as usize].p),
+                    )
+                };
+
+                verts.push(SDVertex {
+                    p: new_p,
+                    start_face: faces[fi as usize].children[3],
+                    child: NULL,
+                    regular: true,
+                    boundary: is_boundary,
+                });
+                new_verts.push(new_vi);
+                edge_verts.insert(key, new_vi);
+            }
+        }
+
+        // Snapshot parent arrays before writing into children to avoid
+        // reading stale data from the same Vec.
+
+        // Even vertex start_face pointers
+        for &vi in &old_verts {
+            let sf = verts[vi as usize].start_face;
+            let vnum = faces[sf as usize].vnum(vi);
+            let child = verts[vi as usize].child;
+            verts[child as usize].start_face = faces[sf as usize].children[vnum];
+        }
+
+        // Face neighbor pointers
+        for &fi in &old_faces {
+            let children = faces[fi as usize].children;
+            let face_f = faces[fi as usize].f;
+            let face_v = faces[fi as usize].v;
+
+            for j in 0..3usize {
+                // Center child (children[3]) neighbors the three corner children
+                faces[children[3] as usize].f[j] = children[next(j)];
+                faces[children[j] as usize].f[next(j)] = children[3];
+
+                // Corner child j's f[j] = neighbor face's child at the shared vertex
+                let f2 = face_f[j];
+                faces[children[j] as usize].f[j] = if f2 != NULL {
+                    let vnum2 = faces[f2 as usize].vnum(face_v[j]);
+                    faces[f2 as usize].children[vnum2]
+                } else {
+                    NULL
+                };
+
+                let f2 = face_f[prev(j)];
+                faces[children[j] as usize].f[prev(j)] = if f2 != NULL {
+                    let vnum2 = faces[f2 as usize].vnum(face_v[j]);
+                    faces[f2 as usize].children[vnum2]
+                } else {
+                    NULL
+                };
+            }
+        }
+
+        // Face vertex pointers
+        for &fi in &old_faces {
+            let children = faces[fi as usize].children;
+            let face_v = faces[fi as usize].v;
+
+            for j in 0..3usize {
+                let child_v = verts[face_v[j] as usize].child;
+                let edge_v = edge_verts[&EdgeKey::new(face_v[j], face_v[next(j)])];
+
+                faces[children[j] as usize].v[j] = child_v;
+                faces[children[j] as usize].v[next(j)] = edge_v;
+                faces[children[next(j)] as usize].v[j] = edge_v;
+                faces[children[3] as usize].v[j] = edge_v;
+            }
+        }
+
+        live_verts = new_verts;
+        live_faces = new_faces;
+    }
+
+    let n_final = live_verts.len();
+    let mut p_limit: Vec<Point3f> = Vec::with_capacity(n_final);
+
+    for &vi in &live_verts {
+        if verts[vi as usize].boundary {
+            p_limit.push(weight_boundary(&verts, &faces, vi, 1.0 / 5.0));
+        } else {
+            let val = valence(&verts, &faces, vi);
+            p_limit.push(weight_one_ring(&verts, &faces, vi, loop_gamma(val)));
+        }
+    }
+    for (i, &vi) in live_verts.iter().enumerate() {
+        verts[vi as usize].p = p_limit[i];
+    }
+
+    let pi = std::f64::consts::PI as Float;
+    let mut normals: Vec<Normal3f> = Vec::with_capacity(n_final);
+
+    for &vi in &live_verts {
+        let ring = one_ring(&verts, &faces, vi);
+        let val = ring.len();
+        let vert_p = verts[vi as usize].p;
+
+        let (s, t) = if !verts[vi as usize].boundary {
+            let mut s = Vector3f::default();
+            let mut t = Vector3f::default();
+            for j in 0..val {
+                let angle = 2.0 * pi * j as Float / val as Float;
+                s += angle.cos() * Vector3f::from(ring[j]);
+                t += angle.sin() * Vector3f::from(ring[j]);
+            }
+            (s, t)
+        } else {
+            let s = Vector3f::from(ring[val - 1]) - Vector3f::from(ring[0]);
+            let t = if val == 2 {
+                Vector3f::from(ring[0]) + Vector3f::from(ring[1])
+                    - 2.0 * Vector3f::from(vert_p)
+            } else if val == 3 {
+                Vector3f::from(ring[1]) - Vector3f::from(vert_p)
+            } else if val == 4 {
+                -1.0 * Vector3f::from(ring[0])
+                    + 2.0 * Vector3f::from(ring[1])
+                    + 2.0 * Vector3f::from(ring[2])
+                    - 1.0 * Vector3f::from(ring[3])
+                    - 2.0 * Vector3f::from(vert_p)
+            } else {
+                let theta = pi / (val - 1) as Float;
+                let mut t = theta.sin()
+                    * (Vector3f::from(ring[0]) + Vector3f::from(ring[val - 1]));
+                for k in 1..(val - 1) {
+                    let wt = (2.0 * theta.cos() - 2.0) * (k as Float * theta).sin();
+                    t += wt * Vector3f::from(ring[k]);
+                }
+                -t
+            };
+            (s, t)
+        };
+
+        normals.push(Normal3f::from(s.cross(t)));
+    }
+
+    let mut vert_remap: HashMap<u32, i32> = HashMap::with_capacity(n_final);
+    for (i, &vi) in live_verts.iter().enumerate() {
+        vert_remap.insert(vi, i as i32);
+    }
+
+    let mut indices: Vec<i32> = Vec::with_capacity(live_faces.len() * 3);
+    for &fi in &live_faces {
+        for j in 0..3 {
+            indices.push(vert_remap[&faces[fi as usize].v[j]]);
+        }
+    }
+
+    TriangleMesh::new(
+        render_from_object,
+        reverse_orientation,
+        &indices,
+        &p_limit,
+        &normals,
+        &[],
+        &[],
+        &[],
+    )
+}

src/utils/mod.rs

@@ -4,6 +4,7 @@ pub mod containers;
 pub mod error;
 pub mod file;
 pub mod io;
+pub mod loopsubdiv;
 pub mod mipmap;
 pub mod parallel;
 pub mod parameters;
@@ -12,11 +13,12 @@ pub mod upload;
 
 pub use error::FileLoc;
 pub use file::{read_float_file, resolve_filename};
+pub use loopsubdiv::*;
+pub use mipmap::{MIPMap, MIPMapFilterOptions};
 pub use parameters::{
     ParameterDictionary, ParsedParameter, ParsedParameterVector, TextureParameterDictionary,
 };
-pub use mipmap::{MIPMap, MIPMapFilterOptions};
+pub use upload::{ArenaUpload, Upload};
-pub use upload::{Upload, ArenaUpload};
 
 #[cfg(feature = "vulkan")]
 pub type Arena = arena::Arena<backend::vulkan::VulkanAllocator>;

src/utils/parameters.rs

@@ -2,7 +2,7 @@ use crate::core::spectrum::SPECTRUM_FILE_CACHE;
 use crate::core::texture::{FloatTexture, SpectrumTexture};
 use crate::spectra::data::get_named_spectrum;
 use crate::spectra::piecewise::ReadFromFile;
-use crate::utils::FileLoc;
+use crate::utils::{FileLoc, resolve_filename};
 use anyhow::{bail, Result};
 use shared::core::color::RGB;
 use shared::core::geometry::{Normal3f, Point2f, Point3f, Vector2f, Vector3f};
@@ -597,7 +597,13 @@ impl ParameterDictionary {
         match param.type_name.as_str() {
             "rgb" | "color" => self.extract_rgb_spectrum(param, spectrum_type),
             "blackbody" => self.extract_file_spectrum(param),
-            "spectrum" => self.extract_sampled_spectrum(param),
+            "spectrum" => {
+                if !param.strings.is_empty() {
+                    self.extract_file_spectrum(param)
+                } else {
+                    self.extract_sampled_spectrum(param)
+                }
+            }
             _ => Vec::new(),
         }
     }
@@ -695,7 +701,10 @@ impl ParameterDictionary {
             .map(|s| {
                 get_named_spectrum(s)
                     .ok_or(())
-                    .or_else(|_| read_spectrum_from_file(s).map_err(|_| ()))
+                    .or_else(|_| {
+                        let resolved = resolve_filename(s);
+                        read_spectrum_from_file(&resolved).map_err(|_| ())
+                    })
                     .unwrap_or_else(|_| panic!("{}: {}: unable to read spectrum", param.loc, s))
             })
             .collect()

src/utils/parser.rs

@@ -7,9 +7,9 @@ use std::io::{self, Read};
 use std::path::{Path, PathBuf};
 use std::sync::Arc;
 
-use crate::Arena;
 use crate::utils::error::FileLoc;
 use crate::utils::parameters::{ParameterDictionary, ParsedParameter, ParsedParameterVector};
+use crate::Arena;
 use shared::Float;
 
 pub trait ParserTarget {
@@ -112,7 +112,7 @@ pub trait ParserTarget {
         loc: FileLoc,
     ) -> Result<(), ParserError>;
 
-    fn world_begin(&mut self, loc: FileLoc, arena: Arc<Arena>) -> Result<(), ParserError>;
+    fn world_begin(&mut self, loc: FileLoc, arena: &Arena) -> Result<(), ParserError>;
     fn attribute_begin(&mut self, loc: FileLoc) -> Result<(), ParserError>;
     fn attribute_end(&mut self, loc: FileLoc) -> Result<(), ParserError>;
     fn attribute(
@@ -576,7 +576,7 @@ impl ParserTarget for FormattingParserTarget {
         Ok(())
     }
 
-    fn world_begin(&mut self, _loc: FileLoc, _arena: Arc<Arena>) -> Result<(), ParserError> {
+    fn world_begin(&mut self, _loc: FileLoc, _arena: &Arena) -> Result<(), ParserError> {
         println!("{}WorldBegin", self.indent(0));
         self.cat_indent_count += 4;
         Ok(())
@@ -886,9 +886,11 @@ impl<'a> SceneParser<'a> {
             let mut val_type = match type_name {
                 "integer" => ValType::Int,
                 "bool" => ValType::Bool,
-                "float" | "point" | "vector" | "normal" | "color" | "spectrum" | "rgb"
+                "float" | "point" | "vector" | "normal" | "color" | "rgb" | "blackbody" => {
-                | "blackbody" => ValType::Float,
+                    ValType::Float
+                }
                 "string" | "texture" => ValType::String,
+                "spectrum" => ValType::Unknown,
                 _ => ValType::Unknown,
             };
 
@@ -976,8 +978,7 @@ impl<'a> SceneParser<'a> {
         }
     }
 
-    pub fn run(&mut self) -> Result<(), ParserError> {
+    pub fn run(&mut self, arena: Arc<Arena>) -> Result<(), ParserError> {
-        let arena = Arc::new(Arena::default());
         loop {
             let token = match self.next_token()? {
                 Some(t) => t,
@@ -1264,7 +1265,7 @@ impl<'a> SceneParser<'a> {
                 },
 
                 'W' => match token.text.as_str() {
-                    "WorldBegin" => self.target.world_begin(token.loc, arena.clone())?,
+                    "WorldBegin" => self.target.world_begin(token.loc, &arena)?,
                     "WorldEnd" => {}
                     _ => {
                         return Err(ParserError::Generic(