A lot of fixes, rendering is somewhat working. CAREFUL: REMOVED as_ref, need to sub by get but got lazy

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.radius()
     }
     fn integral(&self) -> Float {
-        unsafe { self.as_ref().integral() }
+        self.as_ref().integral()
     }
  
     fn evaluate(&self, p: Point2f) -> Float {
-        unsafe { self.as_ref().evaluate(p) }
+        self.evaluate(p)
     }
 
     fn sample(&self, p: Point2f) -> FilterSample {
-        unsafe { self.as_ref().sample(p) }
+        self.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

@@ -590,13 +590,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 {

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.evaluate(lambda)
     }
     fn max_value(&self) -> Float {
-        unsafe { self.as_ref().max_value() }
+        self.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/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,3 +1,4 @@
+use core::hash::{Hash, Hasher};
 use core::marker::PhantomData;
 use core::ops::Index;
 
@@ -25,6 +26,12 @@ 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 +51,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 +60,15 @@ impl<T> Ptr<T> {
         }
     }
 
+    #[inline(always)]
+    pub unsafe 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 +108,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(
@@ -61,7 +61,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)
     }

src/core/render.rs

@@ -13,29 +13,36 @@ 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)?;
+    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 +82,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() } {
                             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,30 @@ 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 +139,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,
@@ -343,7 +380,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 +395,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 +563,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 +721,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" {
@@ -739,12 +776,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 +823,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 +873,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/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;
@@ -227,6 +227,12 @@ impl RayIntegratorTrait for PathIntegrator {
             let t_hit = si.t_hit();
             let isect = &mut si.intr;
 
+            if state.depth == 0 {
+                let n = isect.n();
+                let v = (n.x() + 1.0) * 0.5;
+                return (SampledSpectrum::from_array(&[v, v, v, v]), None);
+            }
+
             // Emission from hit surface
             let le = isect.le(-ray.d, lambda);
             if !le.is_black() {
@@ -308,6 +314,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);

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/utils/loopsubdiv.rs

@@ -0,0 +1,534 @@
+use crate::utils::backend::GpuAllocator;
+use shared::core::geometry::{Normal3f, Point3f, Vector3f};
+use shared::{Float, Ptr};
+use std::collections::HashMap;
+
+pub struct SDVertex {
+    pub p: Point3f,
+    pub start_face: Ptr<SDFace>,
+    pub child: Ptr<SDVertex>,
+    pub regular: bool,
+    pub boundary: bool,
+}
+
+impl SDVertex {
+    pub fn new(p: Point3f) -> Self {
+        Self {
+            p,
+            start_face: Ptr::null(),
+            child: Ptr::null(),
+            regular: false,
+            boundary: false,
+        }
+    }
+}
+
+pub struct SDFace {
+    pub v: [Ptr<SDVertex>; 3],
+    pub f: [Ptr<SDFace>; 3],
+    pub children: [Ptr<SDFace>; 4],
+}
+
+impl SDFace {
+    pub fn new(v0: Ptr<SDVertex>, v1: Ptr<SDVertex>, v2: Ptr<SDVertex>) -> Self {
+        Self {
+            v: [v0, v1, v2],
+            f: [Ptr::null(); 3],
+            children: [Ptr::null(); 4],
+        }
+    }
+
+    pub fn vnum(&self, vert: Ptr<SDVertex>) -> usize {
+        for i in 0..3 {
+            if self.v[i] == vert {
+                return i;
+            }
+        }
+        panic!("Basic logic error in SDFace::vnum()");
+    }
+
+    pub fn next_face(&self, vert: Ptr<SDVertex>) -> Ptr<SDFace> {
+        self.f[self.vnum(vert)]
+    }
+
+    pub fn prev_face(&self, vert: Ptr<SDVertex>) -> Ptr<SDFace> {
+        self.f[prev(self.vnum(vert))]
+    }
+
+    pub fn next_vert(&self, vert: Ptr<SDVertex>) -> Ptr<SDVertex> {
+        self.v[next(self.vnum(vert))]
+    }
+
+    pub fn prev_vert(&self, vert: Ptr<SDVertex>) -> Ptr<SDVertex> {
+        self.v[prev(self.vnum(vert))]
+    }
+
+    pub fn other_vert(&self, v0: Ptr<SDVertex>, v1: Ptr<SDVertex>) -> Ptr<SDVertex> {
+        for i in 0..3 {
+            if self.v[i] != v0 && self.v[i] != v1 {
+                return self.v[i];
+            }
+        }
+        panic!("Basic logic error in SDFace::other_vert()");
+    }
+}
+
+const fn next(i: usize) -> usize {
+    (i + 1) % 3
+}
+
+const fn prev(i: usize) -> usize {
+    (i + 2) % 3
+}
+
+fn canonical_key<T>(v0: Ptr<T>, v1: Ptr<T>) -> (usize, usize) {
+    let a = v0.as_ptr() as usize;
+    let b = v1.as_ptr() as usize;
+    if a <= b {
+        (a, b)
+    } else {
+        (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(v: Ptr<SDVertex>) -> usize {
+    let vertex = v;
+    let start = vertex.start_face;
+    if start.is_null() {
+        return 0;
+    }
+
+    if !vertex.boundary {
+        let mut nf = 1;
+        let mut f = start;
+        loop {
+            f = f.next_face(v);
+            if f.is_null() {
+                panic!("interior vertex with broken face loop");
+            }
+            if f == start {
+                break;
+            }
+            nf += 1;
+        }
+        nf
+    } else {
+        let mut nf = 1;
+        let mut f = start;
+        while !f.next_face(v).is_null() {
+            f = f.next_face(v);
+            nf += 1;
+        }
+        f = start;
+        while !f.prev_face(v).is_null() {
+            f = f.prev_face(v);
+            nf += 1;
+        }
+        nf + 1
+    }
+}
+
+fn one_ring(v: Ptr<SDVertex>, out: &mut [Point3f]) {
+    let vertex = v;
+    if !vertex.boundary {
+        let start = vertex.start_face;
+        let mut face = start;
+        let mut i = 0;
+        loop {
+            out[i] = face.next_vert(v).p;
+            i += 1;
+            face = face.next_face(v);
+            if face == start {
+                break;
+            }
+        }
+    } else {
+        let start = vertex.start_face;
+        let mut face = start;
+        while !face.next_face(v).is_null() {
+            face = face.next_face(v);
+        }
+        out[0] = face.next_vert(v).p;
+        let mut i = 1;
+        loop {
+            out[i] = face.prev_vert(v).p;
+            i += 1;
+            let pf = face.prev_face(v);
+            if pf.is_null() {
+                break;
+            }
+            face = pf;
+        }
+    }
+}
+
+fn weight_one_ring(v: Ptr<SDVertex>, beta: Float) -> Point3f {
+    let valence = valence(v);
+    let mut p_ring = vec![Point3f::default(); valence];
+    one_ring(v, &mut p_ring);
+
+    let vertex = v;
+    let mut p = vertex.p * (1.0 - valence as Float * beta);
+    for i in 0..valence {
+        p += p_ring[i] * beta;
+    }
+    p
+}
+
+fn weight_boundary(v: Ptr<SDVertex>, beta: Float) -> Point3f {
+    let valence = valence(v);
+    let mut p_ring = vec![Point3f::default(); valence];
+    one_ring(v, &mut p_ring);
+
+    let vertex = v;
+    let mut p = vertex.p * (1.0 - 2.0 * beta);
+    p += p_ring[0] * beta;
+    p += p_ring[valence - 1] * beta;
+    p
+}
+
+pub fn loop_subdivide<A: GpuAllocator>(
+    arena: &Arena<A>,
+    render_from_object: &Transform,
+    reverse_orientation: bool,
+    n_levels: usize,
+    vertex_indices: &[usize],
+    p: &[Point3f],
+) -> TriangleMesh {
+    let mut vertices: Vec<Ptr<SDVertex>> = Vec::with_capacity(p.len());
+    for &pos in p {
+        vertices.push(arena.alloc(SDVertex::new(pos)));
+    }
+
+    let n_faces = vertex_indices.len() / 3;
+    let mut faces: Vec<Ptr<SDFace>> = Vec::with_capacity(n_faces);
+    for i in 0..n_faces {
+        let v0 = vertices[vertex_indices[3 * i]];
+        let v1 = vertices[vertex_indices[3 * i + 1]];
+        let v2 = vertices[vertex_indices[3 * i + 2]];
+        faces.push(arena.alloc(SDFace::new(v0, v1, v2)));
+    }
+
+    for &f in &faces {
+        let face = f;
+        for j in 0..3 {
+            face.v[j].get_mut().unwrap().start_face = f;
+        }
+    }
+
+    let mut edges: HashMap<(usize, usize), (Ptr<SDFace>, usize)> = HashMap::new();
+    for &f in &faces {
+        let face = f;
+        for edge_num in 0..3 {
+            let v0 = face.v[edge_num];
+            let v1 = face.v[next(edge_num)];
+            let key = canonical_key(v0, v1);
+
+            if let Some(&(first_face, first_edge_num)) = edges.get(&key) {
+                first_face.get_mut().unwrap().f[first_edge_num] = f;
+                f.get_mut().f[edge_num] = first_face;
+                edges.remove(&key);
+            } else {
+                edges.insert(key, (f, edge_num));
+            }
+        }
+    }
+
+    for &v in &vertices {
+        let vertex = v.get_mut().unwrap();
+        let start = vertex.start_face;
+        let mut f = start;
+        let mut is_boundary = false;
+        loop {
+            let nf = f.next_face(v);
+            if nf.is_null() {
+                is_boundary = true;
+                break;
+            }
+            f = nf;
+            if f == start {
+                break;
+            }
+        }
+
+        let val = valence(v);
+        vertex.boundary = is_boundary;
+        vertex.regular = if !is_boundary && val == 6 {
+            true
+        } else if is_boundary && val == 4 {
+            true
+        } else {
+            false
+        };
+    }
+
+    // -- Subdivision levels ---------------------------------------------------
+    let mut f = faces;
+    let mut v = vertices;
+
+    for _ in 0..n_levels {
+        let mut new_faces: Vec<Ptr<SDFace>> = Vec::new();
+        let mut new_vertices: Vec<Ptr<SDVertex>> = Vec::new();
+
+        // Allocate vertex children
+        for &vertex_ptr in &v {
+            let vertex = vertex_ptr;
+            let child = arena.alloc(SDVertex {
+                p: Point3f::default(),
+                start_face: Ptr::null(),
+                child: Ptr::null(),
+                regular: vertex.regular,
+                boundary: vertex.boundary,
+            });
+            vertex_ptr.get_mut().unwrap().child = child;
+            new_vertices.push(child);
+        }
+
+        for &face_ptr in &f {
+            for k in 0..4 {
+                let child = arena.alloc(SDFace::new(Ptr::null(), Ptr::null(), Ptr::null()));
+                face_ptr.get_mut().unwrap().children[k] = child;
+                new_faces.push(child);
+            }
+        }
+
+        for &vertex_ptr in &v {
+            let vertex = vertex_ptr;
+            let child = vertex.child;
+            let new_p = if !vertex.boundary {
+                if vertex.regular {
+                    weight_one_ring(vertex_ptr, 1.0 / 16.0)
+                } else {
+                    let val = valence(vertex_ptr);
+                    weight_one_ring(vertex_ptr, beta(val))
+                }
+            } else {
+                weight_boundary(vertex_ptr, 1.0 / 8.0)
+            };
+            child.get_mut().unwrap().p = new_p;
+        }
+
+        let mut edge_verts: HashMap<(usize, usize), Ptr<SDVertex>> = HashMap::new();
+        for &face_ptr in &f {
+            let face = face_ptr;
+            for k in 0..3 {
+                let v0 = face.v[k];
+                let v1 = face.v[next(k)];
+                let key = canonical_key(v0, v1);
+
+                if !edge_verts.contains_key(&key) {
+                    let is_boundary = face.f[k].is_null();
+                    let start_face = face.children[3];
+
+                    let new_p = if is_boundary {
+                        let p0 = v0.p;
+                        let p1 = v1.p;
+                        p0 * 0.5 + p1 * 0.5
+                    } else {
+                        let p0 = v0.p;
+                        let p1 = v1.p;
+                        let p2 = face.other_vert(v0, v1).p;
+                        let f2 = face.f[k];
+                        let p3 = f2.other_vert(v0, v1).p;
+                        p0 * (3.0 / 8.0) + p1 * (3.0 / 8.0) + p2 * (1.0 / 8.0) + p3 * (1.0 / 8.0)
+                    };
+
+                    let vert = arena.alloc(SDVertex {
+                        p: new_p,
+                        start_face,
+                        child: Ptr::null(),
+                        regular: true,
+                        boundary: is_boundary,
+                    });
+                    edge_verts.insert(key, vert);
+                    new_vertices.push(vert);
+                }
+            }
+        }
+
+        for &vertex_ptr in &v {
+            let vertex = vertex_ptr;
+            let start_face = vertex.start_face;
+            let vert_num = start_face.vnum(vertex_ptr);
+            let child = vertex.child;
+            let child_start_face = start_face.children[vert_num];
+            child.get_mut().unwrap().start_face = child_start_face;
+        }
+
+        for &face_ptr in &f {
+            let face = face_ptr;
+            // Copy the data we need so we don't hold a borrow into the parent
+            // while we mutably borrow its children.
+            let face_children = face.children;
+            let face_f = face.f;
+            let face_v = face.v;
+
+            for j in 0..3 {
+                let c3 = face_children[3];
+                let c_j = face_children[j];
+                let c_next_j = face_children[next(j)];
+
+                // Sibling links
+                if !c3.is_null() {
+                    c3.get_mut().unwrap().f[j] = c_next_j;
+                }
+                if !c_j.is_null() {
+                    c_j.get_mut().unwrap().f[next(j)] = c3;
+                }
+
+                // Neighbor links
+                let f2 = face_f[j];
+                if !f2.is_null() {
+                    let f2_vnum = f2.vnum(face_v[j]);
+                    let f2_child = f2.children[f2_vnum];
+                    if !c_j.is_null() {
+                        c_j.get_mut().unwrap().f[j] = f2_child;
+                    }
+                } else if !c_j.is_null() {
+                    c_j.get_mut().unwrap().f[j] = Ptr::null();
+                }
+
+                let f2_prev = face_f[prev(j)];
+                if !f2_prev.is_null() {
+                    let f2_vnum = f2_prev.vnum(face_v[j]);
+                    let f2_child = f2_prev.children[f2_vnum];
+                    if !c_j.is_null() {
+                        c_j.get_mut().unwrap().f[prev(j)] = f2_child;
+                    }
+                } else if !c_j.is_null() {
+                    c_j.get_mut().unwrap().f[prev(j)] = Ptr::null();
+                }
+            }
+        }
+
+        // Update face vertex pointers
+        for &face_ptr in &f {
+            let face = face_ptr;
+            let face_v = face.v;
+            let face_children = face.children;
+
+            for j in 0..3 {
+                let v_j_child = face_v[j].child;
+                let c_j = face_children[j];
+                if !c_j.is_null() {
+                    c_j.get_mut().unwrap().v[j] = v_j_child;
+                }
+
+                let key = canonical_key(face_v[j], face_v[next(j)]);
+                let vert = edge_verts[&key];
+
+                let c_next_j = face_children[next(j)];
+                let c3 = face_children[3];
+
+                if !c_j.is_null() {
+                    c_j.get_mut().unwrap().v[next(j)] = vert;
+                }
+                if !c_next_j.is_null() {
+                    c_next_j.get_mut().unwrap().v[j] = vert;
+                }
+                if !c3.is_null() {
+                    c3.get_mut().unwrap().v[j] = vert;
+                }
+            }
+        }
+
+        // Prepare for next level
+        f = new_faces;
+        v = new_vertices;
+    }
+
+    let mut p_limit: Vec<Point3f> = Vec::with_capacity(v.len());
+    for &vertex_ptr in &v {
+        let vertex = vertex_ptr;
+        if vertex.boundary {
+            p_limit.push(weight_boundary(vertex_ptr, 1.0 / 5.0));
+        } else {
+            let val = valence(vertex_ptr);
+            p_limit.push(weight_one_ring(vertex_ptr, loop_gamma(val)));
+        }
+    }
+    for (i, &vertex_ptr) in v.iter().enumerate() {
+        vertex_ptr.get_mut().unwrap().p = p_limit[i];
+    }
+
+    let mut ns: Vec<Normal3f> = Vec::with_capacity(v.len());
+    let mut p_ring: Vec<Point3f> = Vec::with_capacity(16);
+
+    for &vertex_ptr in &v {
+        let vertex = vertex_ptr;
+        let valence = valence(vertex_ptr);
+        p_ring.resize(valence, Point3f::default());
+        one_ring(vertex_ptr, &mut p_ring[..valence]);
+
+        let mut s = Vector3f::default();
+        let mut t = Vector3f::default();
+
+        if !vertex.boundary {
+            for j in 0..valence {
+                let angle = 2.0 * std::f64::consts::PI * j as f64 / valence as f64;
+                s += Vector3f::from(p_ring[j]) * angle.cos() as Float;
+                t += Vector3f::from(p_ring[j]) * angle.sin() as Float;
+            }
+        } else {
+            s = Vector3f::from(p_ring[valence - 1] - p_ring[0]);
+            if valence == 2 {
+                t = Vector3f::from(p_ring[0] + p_ring[1] - vertex_ptr.p * 2.0);
+            } else if valence == 3 {
+                t = Vector3f::from(p_ring[1] - vertex_ptr.p);
+            } else if valence == 4 {
+                t = Vector3f::from(
+                    p_ring[0] * -1.0
+                        + p_ring[1] * 2.0
+                        + p_ring[2] * 2.0
+                        + p_ring[3] * -1.0
+                        + vertex_ptr.p * -2.0,
+                );
+            } else {
+                let theta = std::f64::consts::PI / (valence - 1) as f64;
+                t = Vector3f::from(p_ring[0] + p_ring[valence - 1].into()) * theta.sin() as Float;
+                for k in 1..(valence - 1) {
+                    let wt = (2.0 * theta.cos() - 2.0) * ((k as f64) * theta).sin();
+                    t += Vector3f::from(p_ring[k]) * wt as Float;
+                }
+                t = -t;
+            }
+        }
+        ns.push(Normal3f::from(Vector3f::cross(s, t)));
+    }
+
+    let ntris = f.len();
+    let mut verts: Vec<usize> = Vec::with_capacity(3 * ntris);
+    let mut used_verts: HashMap<Ptr<SDVertex>, usize> = HashMap::new();
+    for (i, &vertex_ptr) in v.iter().enumerate() {
+        used_verts.insert(vertex_ptr, i);
+    }
+    for &face_ptr in &f {
+        let face = face_ptr;
+        for j in 0..3 {
+            verts.push(used_verts[&face.v[j]]);
+        }
+    }
+
+    TriangleMesh::new(
+        render_from_object,
+        reverse_orientation,
+        &verts,
+        &p_limit,
+        &[],
+        &ns,
+        &[],
+        &[],
+    )
+}

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

@@ -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(())
@@ -976,8 +976,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 +1263,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(