Fixing logic bugs in Sampling, texture loading. Moving on to issues with ImageTexture, the long dreaded dive

.gitignore

@@ -15,3 +15,7 @@ tests/
 scenes/
 compile.sh
 output/
+*.md
+!README.md
+!INSTALL.md
+docs/

shared/src/bxdfs/diffuse.rs

@@ -41,7 +41,7 @@ impl BxDFTrait for DiffuseBxDF {
             return None;
         }
         let mut wi = sample_cosine_hemisphere(u);
-        if wo.z() == 0. {
+        if wo.z() < 0. {
             wi[2] *= -1.;
         }
         let pdf = cosine_hemisphere_pdf(abs_cos_theta(wi));

shared/src/cameras/orthographic.rs

@@ -114,7 +114,7 @@ impl CameraTrait for OrthographicCamera {
 
         Some(CameraRay {
             ray: camera_ray,
-            weight: SampledSpectrum::default(),
+            weight: SampledSpectrum::new(1.),
         })
     }
 

shared/src/cameras/spherical.rs

@@ -58,7 +58,7 @@ impl CameraTrait for SphericalCamera {
         );
         Some(CameraRay {
             ray: self.render_from_camera(&ray, &mut None),
-            weight: SampledSpectrum::default(),
+            weight: SampledSpectrum::new(1.),
         })
     }
 }

shared/src/core/bssrdf.rs

@@ -86,7 +86,7 @@ impl From<&SubsurfaceInteraction> for SurfaceInteraction {
             dvdx: 0.,
             dudy: 0.,
             dvdy: 0.,
-            shape: Ptr::from(&Shape::default()),
+            shape: Ptr::null(),
         }
     }
 }

shared/src/core/color.rs

@@ -8,6 +8,7 @@ use core::fmt;
 use core::ops::{
     Add, AddAssign, Div, DivAssign, Index, IndexMut, Mul, MulAssign, Neg, Sub, SubAssign,
 };
+use anyhow::{Result, bail};
 use enum_dispatch::enum_dispatch;
 use num_traits::Float as NumFloat;
 
@@ -680,6 +681,16 @@ pub enum ColorEncoding {
     SRGB(SRGBEncoding),
 }
 
+impl ColorEncoding {
+    pub fn from_name(name: &str) -> Result<Self> {
+        match name {
+            "sRGB" | "srgb" => Ok(ColorEncoding::SRGB(SRGBEncoding)),
+            "linear" => Ok(ColorEncoding::Linear(LinearEncoding)),
+            _ => bail!("Unknown color encoding: {}", name),
+        }
+    }
+}
+
 impl fmt::Display for ColorEncoding {
     fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
         write!(f, "Encoding")

shared/src/core/geometry/bounds.rs

@@ -5,6 +5,7 @@ use crate::core::geometry::{max, min};
 use crate::utils::gpu_array_from_fn;
 use crate::utils::interval::Interval;
 use crate::utils::math::lerp;
+use crate::{gamma, gamma_t};
 use core::mem;
 use core::ops::{Add, Div, DivAssign, Mul, Sub};
 use num_traits::{Bounded, Num};
@@ -220,7 +221,7 @@ where
         (center, radius)
     }
 
-    pub fn insersect(&self, o: Point3<T>, d: Vector3<T>, t_max: T) -> Option<(T, T)> {
+    pub fn intersect(&self, o: Point3<T>, d: Vector3<T>, t_max: T) -> Option<(T, T)> {
         let mut t0 = T::zero();
         let mut t1 = t_max;
 
@@ -231,6 +232,8 @@ where
             if t_near > t_far {
                 mem::swap(&mut t_near, &mut t_far);
             }
+
+            t_far = t_far * (T::one() + (T::one() + T::one()) * gamma_t::<T>(3));
             t0 = if t_near > t0 { t_near } else { t0 };
             t1 = if t_far < t1 { t_far } else { t1 };
             if t0 > t1 {
@@ -274,7 +277,10 @@ impl Bounds3f {
 
         // Check Y
         let ty_min = (bounds[dir_is_neg[1]].y() - o.y()) * inv_dir.y();
-        let ty_max = (bounds[1 - dir_is_neg[1]].y() - o.y()) * inv_dir.y();
+        let mut ty_max = (bounds[1 - dir_is_neg[1]].y() - o.y()) * inv_dir.y();
+
+        t_max = t_max * (1. + 2. * gamma(3));
+        ty_max = ty_max * (1. + 2. * gamma(3));
 
         if t_min > ty_max || ty_min > t_max {
             return None;
@@ -288,7 +294,8 @@ impl Bounds3f {
 
         // Check Z
         let tz_min = (bounds[dir_is_neg[2]].z() - o.z()) * inv_dir.z();
-        let tz_max = (bounds[1 - dir_is_neg[2]].z() - o.z()) * inv_dir.z();
+        let mut tz_max = (bounds[1 - dir_is_neg[2]].z() - o.z()) * inv_dir.z();
+        tz_max = tz_max * (1. + 2. * gamma(3));
 
         if t_min > tz_max || tz_min > t_max {
             return None;
@@ -321,7 +328,11 @@ impl Bounds3f {
         let mut t_min = (bounds[dir_is_neg[0]].x() - o.x()) * inv_dir.x();
         let mut t_max = (bounds[1 - dir_is_neg[0]].x() - o.x()) * inv_dir.x();
         let ty_min = (bounds[dir_is_neg[1]].y() - o.y()) * inv_dir.y();
-        let ty_max = (bounds[1 - dir_is_neg[1]].y() - o.y()) * inv_dir.y();
+        let mut ty_max = (bounds[1 - dir_is_neg[1]].y() - o.y()) * inv_dir.y();
+
+
+        t_max = t_max * (1. + 2. * gamma(3));
+        ty_max = ty_max * (1. + 2. * gamma(3));
 
         if t_min > ty_max || ty_min > t_max {
             return false;
@@ -334,7 +345,8 @@ impl Bounds3f {
         }
 
         let tz_min = (bounds[dir_is_neg[2]].z() - o.z()) * inv_dir.z();
-        let tz_max = (bounds[1 - dir_is_neg[2]].z() - o.z()) * inv_dir.z();
+        let mut tz_max = (bounds[1 - dir_is_neg[2]].z() - o.z()) * inv_dir.z();
+        tz_max = tz_max * (1. + 2. * gamma(3));
 
         if t_min > tz_max || tz_min > t_max {
             return false;

shared/src/core/geometry/cone.rs

@@ -52,13 +52,13 @@ impl DirectionCone {
                 * Vector3f::new(
                     w.x()
                         * (wp.y() * w.y() + wp.z() * w.z()
-                            - wp.x() * (square(w.y() + square(w.z())))),
+                            - wp.x() * (square(w.y()) + square(w.z()))),
                     w.y()
                         * (wp.x() * w.x() + wp.z() * w.z()
-                            - wp.y() * (square(w.x() + square(w.z())))),
+                            - wp.y() * (square(w.x()) + square(w.z()))),
                     w.z()
                         * (wp.x() * w.x() + wp.y() * w.y()
-                            - wp.z() * (square(w.x() + square(w.y())))),
+                            - wp.z() * (square(w.x()) + square(w.y()))),
                 )
     }
 
@@ -91,10 +91,10 @@ impl DirectionCone {
         let theta_b = safe_acos(b.cos_theta);
         let theta_d = a.w.angle_between(b.w);
 
-        if (theta_d + theta_b).min(PI) <= theta_b {
+        if (theta_d + theta_b).min(PI) <= theta_a {
             return a.clone();
         }
-        if (theta_d + theta_a).min(PI) <= theta_a {
+        if (theta_d + theta_a).min(PI) <= theta_b {
             return b.clone();
         }
 
@@ -107,7 +107,7 @@ impl DirectionCone {
         // Find the merged cone's axis and return cone union
         let theta_r = theta_o - theta_a;
         let wr = a.w.cross(b.w);
-        if wr.norm_squared() >= 0. {
+        if wr.norm_squared() == 0. {
             return DirectionCone::entire_sphere();
         }
 

shared/src/core/geometry/primitives.rs

@@ -2,13 +2,13 @@ use super::traits::{SqrtExt, Tuple, VectorLike};
 use super::{Float, NumFloat, PI};
 use crate::utils::interval::Interval;
 use crate::utils::math::{clamp, difference_of_products, quadratic, safe_asin};
+use core::fmt;
 use core::hash::{Hash, Hasher};
 use core::iter::Sum;
 use core::ops::{
     Add, AddAssign, Div, DivAssign, Index, IndexMut, Mul, MulAssign, Neg, Sub, SubAssign,
 };
 use num_traits::{AsPrimitive, FloatConst, Num, Signed, Zero};
-use core::fmt;
 
 pub trait MulAdd<M = Self, A = Self> {
     type Output;
@@ -23,6 +23,14 @@ impl MulAdd<Float, Float> for Float {
     }
 }
 
+impl MulAdd<f64, f64> for f64 {
+    type Output = f64;
+    #[inline(always)]
+    fn mul_add(self, multiplier: f64, addend: f64) -> Self::Output {
+        num_traits::Float::mul_add(self, multiplier, addend)
+    }
+}
+
 // N-dimensional displacement
 #[repr(C)]
 #[derive(Debug, Copy, Clone, PartialEq, Eq)]
@@ -217,6 +225,27 @@ macro_rules! impl_tuple_core {
     };
 }
 
+#[macro_export]
+macro_rules! impl_num_zero {
+    ($Struct:ident) => {
+        impl<T, const N: usize> num_traits::Zero for $Struct<T, N>
+        where
+            T: num_traits::Zero + Copy + PartialEq,
+        {
+            #[inline]
+            fn zero() -> Self {
+                Self([T::zero(); N])
+            }
+            #[inline]
+            fn is_zero(&self) -> bool {
+                self.0.iter().all(|c| c.is_zero())
+            }
+        }
+    };
+}
+impl_num_zero!(Vector);
+impl_num_zero!(Normal);
+
 #[macro_export]
 macro_rules! impl_scalar_ops {
     ($Struct:ident) => {
@@ -607,33 +636,33 @@ impl<T: Copy> Vector4<T> {
 // Vector operations
 impl<T> Vector3<T>
 where
-    T: Num + Copy + Neg<Output = T>,
+    T: Num + Copy + Neg<Output = T> + Zero + MulAdd<T, T, Output = T>,
 {
     pub fn cross(self, rhs: Self) -> Self {
         Self([
-            self[1] * rhs[2] - self[2] * rhs[1],
-            self[2] * rhs[0] - self[0] * rhs[2],
-            self[0] * rhs[1] - self[1] * rhs[0],
+            difference_of_products(self[1], rhs[2], self[2], rhs[1]),
+            difference_of_products(self[2], rhs[0], self[0], rhs[2]),
+            difference_of_products(self[0], rhs[1], self[1], rhs[0]),
         ])
     }
 }
 
 impl<T> Normal3<T>
 where
-    T: Num + Copy + Neg<Output = T>,
+    T: Num + Copy + Neg<Output = T> + Zero + MulAdd<T, T, Output = T>,
 {
     pub fn cross(self, rhs: Self) -> Self {
         Self([
-            self[1] * rhs[2] - self[2] * rhs[1],
-            self[2] * rhs[0] - self[0] * rhs[2],
-            self[0] * rhs[1] - self[1] * rhs[0],
+            difference_of_products(self[1], rhs[2], self[2], rhs[1]),
+            difference_of_products(self[2], rhs[0], self[0], rhs[2]),
+            difference_of_products(self[0], rhs[1], self[1], rhs[0]),
         ])
     }
 }
 
 impl<T> Vector3<T>
 where
-    T: Num + NumFloat + Copy + Neg<Output = T>,
+    T: Num + NumFloat + Copy + Neg<Output = T> + Zero + MulAdd<T, T, Output = T>,
 {
     pub fn coordinate_system(&self) -> (Self, Self)
     where
@@ -663,7 +692,7 @@ where
 
 impl<T> Normal3<T>
 where
-    T: Num + NumFloat + Copy + Neg<Output = T>,
+    T: Num + NumFloat + Copy + Neg<Output = T> + Zero + MulAdd<T, T, Output = T>,
 {
     pub fn coordinate_system(&self) -> (Self, Self)
     where

shared/src/core/geometry/ray.rs

@@ -43,25 +43,21 @@ impl Ray {
         self.o + self.d * t
     }
 
-    pub fn offset_origin(p: &Point3fi, n: &Normal3f, w: &Vector3f) -> Point3f {
-        let d: Float = Vector3f::from(n.abs()).dot(p.error());
-        let normal: Vector3f = Vector3f::from(*n);
-
-        let mut offset = p.midpoint();
-        if w.dot(normal) < 0.0 {
-            offset -= normal * d;
-        } else {
-            offset += normal * d;
+    pub fn offset_origin(pi: &Point3fi, n: &Normal3f, w: &Vector3f) -> Point3f {
+        let d: Float = Vector3f::from(n.abs()).dot(pi.error());
+        let mut disp: Vector3f = Vector3f::from(*n) * d;
+        if w.dot(Vector3f::from(*n)) < 0.0 {
+            disp = -disp;
         }
-
+        let mut po = pi.midpoint() + disp;
         for i in 0..3 {
-            if n[i] > 0.0 {
-                offset[i] = next_float_up(offset[i]);
-            } else if n[i] < 0.0 {
-                offset[i] = next_float_down(offset[i]);
+            if disp[i] > 0.0 {
+                po[i] = next_float_up(po[i]);
+            } else if disp[i] < 0.0 {
+                po[i] = next_float_down(po[i]);
             }
         }
-        offset
+        po
     }
 
     pub fn spawn(pi: &Point3fi, n: &Normal3f, time: Float, d: Vector3f) -> Ray {

shared/src/core/image.rs

@@ -2,10 +2,12 @@ use crate::core::color::{ColorEncoding, ColorEncodingTrait, LINEAR};
 use crate::core::geometry::{Bounds2f, Point2f, Point2fi, Point2i};
 use crate::utils::math::{f16_to_f32_software, lerp, square};
 use crate::{gvec_with_capacity, Float, GVec};
+use anyhow::{Result, bail};
 use core::hash;
 use core::ops::{Deref, DerefMut};
 use num_traits::Float as NumFloat;
 
+#[repr(C)]
 #[derive(Debug, Clone, Copy, PartialEq, Eq, Hash)]
 pub enum WrapMode {
     Black,
@@ -14,6 +16,18 @@ pub enum WrapMode {
     OctahedralSphere,
 }
 
+impl WrapMode {
+    pub fn parse(name: &str) -> Result<WrapMode> {
+        match name {
+            "clamp" => Ok(WrapMode::Clamp),
+            "black" => Ok(WrapMode::Black),
+            "repeat" => Ok(WrapMode::Repeat),
+            "octahedralsphere" => Ok(WrapMode::OctahedralSphere),
+            _ => bail!("{:?}: wrap mode unknown", name)
+        }
+    }
+}
+
 #[derive(Debug, Clone, Copy, PartialEq, Eq, Hash)]
 pub struct WrapMode2D {
     pub uv: [WrapMode; 2],
@@ -212,6 +226,7 @@ pub struct ImageBase {
     pub n_channels: i32,
 }
 
+#[repr(C)]
 #[derive(Clone, Debug)]
 pub struct Image {
     pub format: PixelFormat,

shared/src/core/interaction.rs

@@ -561,6 +561,8 @@ impl SurfaceInteraction {
         self.shading.n = ns;
         if orientation {
             self.common.n = self.n().face_forward(self.shading.n);
+        } else {
+            self.shading.n = self.shading.n.face_forward(self.common.n);
         }
         self.shading.dpdu = dpdus;
         self.shading.dpdv = dpdvs;

shared/src/core/medium.rs

@@ -11,7 +11,7 @@ use crate::utils::containers::SampledGrid;
 use crate::utils::math::{clamp, square};
 use crate::utils::rng::Rng;
 use crate::utils::transform::Transform;
-use crate::{gvec_with_capacity, GVec, Ptr};
+use crate::{gvec_with_capacity, leak, GVec, Ptr};
 use enum_dispatch::enum_dispatch;
 use num_traits::Float as NumFloat;
 
@@ -170,7 +170,7 @@ pub struct RayMajorantSegment {
 }
 
 #[repr(C)]
-#[derive(Clone, Copy, Debug)]
+#[derive(Clone, Debug)]
 pub enum RayMajorantIterator {
     Homogeneous(HomogeneousMajorantIterator),
     DDA(DDAMajorantIterator),
@@ -191,7 +191,7 @@ impl Iterator for RayMajorantIterator {
 }
 
 #[repr(C)]
-#[derive(Clone, Copy, Debug)]
+#[derive(Clone, Debug)]
 pub struct HomogeneousMajorantIterator {
     called: bool,
     seg: RayMajorantSegment,
@@ -224,7 +224,7 @@ impl Iterator for HomogeneousMajorantIterator {
 }
 
 #[repr(C)]
-#[derive(Debug, Clone, Copy)]
+#[derive(Debug, Clone)]
 pub struct DDAMajorantIterator {
     sigma_t: SampledSpectrum,
     t_min: Float,
@@ -727,33 +727,13 @@ impl Default for MediumInterface {
     }
 }
 
-impl From<Medium> for MediumInterface {
-    fn from(medium: Medium) -> Self {
-        Self {
-            inside: Ptr::from(&medium),
-            outside: Ptr::from(&medium),
-        }
-    }
-}
-
-impl From<&Medium> for MediumInterface {
-    fn from(medium: &Medium) -> Self {
-        Self::from(medium.clone())
-    }
-}
-
 impl MediumInterface {
-    pub fn new(inside: &Medium, outside: &Medium) -> Self {
-        Self {
-            inside: Ptr::from(inside),
-            outside: Ptr::from(outside),
+    pub fn new(inside: Ptr<Medium>, outside: Ptr<Medium>) -> Self {
+        Self { inside, outside }
     }
-    }
-
     pub fn empty() -> Self {
         Self::default()
     }
-
     pub fn is_medium_transition(&self) -> bool {
         self.inside != self.outside
     }

shared/src/core/pbrt.rs

@@ -1,6 +1,6 @@
 use crate::core::geometry::Lerp;
 use core::ops::{Add, Mul};
-use num_traits::{Num, PrimInt};
+use num_traits::{Float as NumFloat, Num, NumCast, PrimInt};
 
 use crate::core::light::LightTrait;
 use crate::core::shape::Shape;
@@ -104,9 +104,16 @@ pub const PI_OVER_2: Float = core::f32::consts::FRAC_PI_2;
 pub const PI_OVER_4: Float = core::f32::consts::FRAC_PI_4;
 pub const SQRT_2: Float = core::f32::consts::SQRT_2;
 
+#[inline]
+pub fn gamma_t<T: NumFloat + NumCast>(n: i32) -> T {
+    let n = T::from(n).unwrap();
+    let eps = T::epsilon() / (T::one() + T::one());
+    n * eps / (T::one() - n * eps)
+}
+
 #[inline]
 pub fn gamma(n: i32) -> Float {
-    n as Float * MACHINE_EPSILON / (1. - n as Float * MACHINE_EPSILON)
+    gamma_t::<Float>(n)
 }
 
 #[cfg(feature = "cpu_debug")]

shared/src/core/sampler.rs

@@ -5,8 +5,7 @@ use crate::utils::math::{
     clamp, encode_morton_2, inverse_radical_inverse, lerp, log2_int,
     owen_scrambled_radical_inverse, permutation_element, radical_inverse, round_up_pow2,
     scrambled_radical_inverse, sobol_interval_to_index, sobol_sample, BinaryPermuteScrambler,
-    DigitPermutation, FastOwenScrambler, NoRandomizer, OwenScrambler, Scrambler,
-    PRIME_TABLE_SIZE,
+    DigitPermutation, FastOwenScrambler, NoRandomizer, OwenScrambler, Scrambler, PRIME_TABLE_SIZE,
 };
 use crate::utils::rng::Rng;
 use crate::utils::sobol::N_SOBOL_DIMENSIONS;
@@ -193,9 +192,10 @@ impl SamplerTrait for HaltonSampler {
     }
 
     fn get1d(&mut self) -> Float {
-        if self.dim >= PRIME_TABLE_SIZE as u32 {
+        if self.dim + 1 >= PRIME_TABLE_SIZE as u32 {
             self.dim = 2;
         }
+        self.dim += 1;
         self.sample_dimension(self.dim)
     }
 
@@ -211,7 +211,7 @@ impl SamplerTrait for HaltonSampler {
     fn get_pixel2d(&mut self) -> Point2f {
         Point2f::new(
             radical_inverse(0, self.halton_index >> self.base_exponents[0]),
-            radical_inverse(1, self.halton_index >> self.base_exponents[1]),
+            radical_inverse(1, self.halton_index / self.base_scales[1]),
         )
     }
 }

shared/src/shapes/cylinder.rs

@@ -9,8 +9,8 @@ use crate::core::shape::{
 use crate::utils::splines::{
     bound_cubic_bezier, cubic_bezier_control_points, evaluate_cubic_bezier, subdivide_cubic_bezier,
 };
-use crate::utils::transform::{Transform, look_at};
-use crate::{Float, PI, gamma};
+use crate::utils::transform::{look_at, Transform};
+use crate::{gamma, Float, PI};
 
 use crate::core::geometry::{SqrtExt, Tuple};
 use crate::utils::interval::Interval;
@@ -61,11 +61,10 @@ impl CylinderShape {
         let di = self
             .object_from_render
             .apply_to_vector_interval(&Vector3fi::new_from_vector(r.d));
-        // Solve quadratic equation to find cylinder t0 and t1 values>>
-        let a: Interval = square(di.x()) + square(di.y()) + square(di.z());
-        let b: Interval = 2. * (di.x() * oi.x() + di.y() * oi.y() + di.z() * oi.z());
-        let c: Interval =
-            square(oi.x()) + square(oi.y()) + square(oi.z()) - square(Interval::new(self.radius));
+        // Solve quadratic equation to find cylinder t0 and t1 values
+        let a: Interval = square(di.x()) + square(di.y());
+        let b: Interval = 2. * (di.x() * oi.x() + di.y() * oi.y());
+        let c: Interval = square(oi.x()) + square(oi.y()) - square(Interval::new(self.radius));
         let f = b / (2. * a);
         let vx: Interval = oi.x() - f * di.x();
         let vy: Interval = oi.y() - f * di.y();

shared/src/shapes/disk.rs

@@ -6,9 +6,10 @@ use crate::core::interaction::{Interaction, InteractionTrait, SurfaceInteraction
 use crate::core::shape::{
     QuadricIntersection, ShapeIntersection, ShapeSample, ShapeSampleContext, ShapeTrait,
 };
-use crate::utils::Transform;
+use crate::utils::interval::Interval;
 use crate::utils::math::square;
 use crate::utils::sampling::sample_uniform_disk_concentric;
+use crate::utils::Transform;
 use crate::{Float, PI};
 use num_traits::Float as NumFloat;
 
@@ -48,25 +49,30 @@ impl DiskShape {
     }
 
     fn basic_intersect(&self, r: &Ray, t_max: Float) -> Option<QuadricIntersection> {
-        let oi = self.object_from_render.apply_to_point(r.o);
-        let di = self.object_from_render.apply_to_vector(r.d);
-        // Reject disk intersections for rays parallel to the disk’s plane
-        if di.z() == 0. {
+        let oi = self
+            .object_from_render
+            .apply_to_interval(&Point3fi::new_from_point(r.o));
+        let di = self
+            .object_from_render
+            .apply_to_vector_interval(&Vector3fi::new_from_vector(r.d));
+
+        if Float::from(di.z()) == 0. {
+            return None;
+        }
+        let t_shape_hit: Interval = (self.height - oi.z()) / di.z();
+        if t_shape_hit.high <= 0. || t_shape_hit.low >= t_max {
             return None;
         }
 
-        let t_shape_hit = (self.height - oi.z()) / di.z();
-        if t_shape_hit == 0. || t_shape_hit >= t_max {
-            return None;
-        }
+        let oi_f = Point3f::from(oi);
+        let di_f = Vector3f::from(di);
+        let t = Float::from(t_shape_hit);
+        let p_hit: Point3f = oi_f + di_f * t;
 
-        // See if hit point is inside disk radii and phi_max
-        let p_hit: Point3f = oi + t_shape_hit * di;
         let dist2 = square(p_hit.x()) + square(p_hit.y());
         if dist2 > square(self.radius) || dist2 < square(self.inner_radius) {
             return None;
         }
-
         let mut phi = p_hit.y().atan2(p_hit.x());
         if phi < 0. {
             phi += 2. * PI;
@@ -76,7 +82,7 @@ impl DiskShape {
         }
 
         Some(QuadricIntersection {
-            t_hit: t_shape_hit,
+            t_hit: t,
             p_obj: p_hit,
             phi,
         })
@@ -105,7 +111,7 @@ impl DiskShape {
         let p_error = Vector3f::zero();
         let flip_normal = self.reverse_orientation ^ self.transform_swap_handedness;
         let wo_object = self.object_from_render.apply_to_vector(wo);
-        SurfaceInteraction::new(
+        let intr = SurfaceInteraction::new(
             Point3fi::new_with_error(p_hit, p_error),
             Point2f::new(u, v),
             wo_object,
@@ -115,7 +121,15 @@ impl DiskShape {
             dndv,
             time,
             flip_normal,
-        )
+        );
+
+        match self
+            .render_from_object
+            .apply_to_interaction(&Interaction::Surface(intr))
+        {
+            Interaction::Surface(si) => si,
+            _ => unreachable!("Only surfaces need apply"),
+        }
     }
 }
 
@@ -197,7 +211,7 @@ impl ShapeTrait for DiskShape {
         }
         wi = wi.normalize();
 
-        ss.pdf = Vector3f::from(ss.intr.n()).dot(-wi).abs() / ctx.p().distance_squared(ss.intr.p());
+        ss.pdf /= Vector3f::from(ss.intr.n()).abs_dot(-wi) / ctx.p().distance_squared(ss.intr.p());
         if ss.pdf.is_infinite() {
             return None;
         }

shared/src/shapes/sphere.rs

@@ -1,17 +1,17 @@
+use crate::core::geometry::{spherical_direction, Frame, SqrtExt};
 use crate::core::geometry::{
     Bounds3f, DirectionCone, Normal3f, Point2f, Point3f, Point3fi, Ray, Vector2f, Vector3f,
     Vector3fi, VectorLike,
 };
-use crate::core::geometry::{Frame, SqrtExt, spherical_direction};
 use crate::core::interaction::{Interaction, InteractionTrait, SurfaceInteraction};
 use crate::core::pbrt::gamma;
 use crate::core::shape::{
     QuadricIntersection, ShapeIntersection, ShapeSample, ShapeSampleContext, ShapeTrait,
 };
-use crate::utils::Transform;
 use crate::utils::interval::Interval;
 use crate::utils::math::{clamp, difference_of_products, radians, safe_acos, safe_sqrt, square};
 use crate::utils::sampling::sample_uniform_sphere;
+use crate::utils::Transform;
 use crate::{Float, PI};
 use num_traits::Float as NumFloat;
 
@@ -86,7 +86,7 @@ impl SphereShape {
         let c: Interval =
             square(oi.x()) + square(oi.y()) + square(oi.z()) - square(Interval::new(self.radius));
 
-        let v: Vector3fi = (oi - b / Vector3fi::from((2. * a) * di)).into();
+        let v: Vector3fi = (oi - b / 2. * a * di).into();
         let length: Interval = v.norm();
         let discrim =
             4. * a * (Interval::new(self.radius) + length) * (Interval::new(self.radius) - length);
@@ -108,7 +108,7 @@ impl SphereShape {
             mem::swap(&mut t0, &mut t1);
         }
 
-        if t0.high > t_max || t1.low < 0. {
+        if t0.high >= t_max || t1.low <= 0. {
             return None;
         }
         let mut t_shape_hit = t0;
@@ -120,6 +120,9 @@ impl SphereShape {
         }
 
         let mut p_hit = Point3f::from(oi) + Float::from(t_shape_hit) * Vector3f::from(di);
+        let scale = self.radius / p_hit.distance(Point3f::new(0., 0., 0.));
+        p_hit = Point3f::from(Vector3f::from(p_hit) * scale);
+
         if p_hit.x() == 0. && p_hit.y() == 0. {
             p_hit[0] = 1e-5 * self.radius;
         }

shared/src/shapes/triangle.rs

@@ -81,7 +81,7 @@ impl TriangleShape {
 
     fn get_uvs(&self) -> Option<[Point2f; 3]> {
         let mesh = self.mesh();
-        if mesh.s.is_empty() {
+        if mesh.uv.is_empty() {
             return None;
         }
         let [v0, v1, v2] = self.get_vertex_indices();

shared/src/textures/image.rs

@@ -1,5 +1,5 @@
-use crate::Float;
 use crate::core::color::{RGB, XYZ};
+use crate::core::image::Image;
 use crate::core::spectrum::SpectrumTrait;
 use crate::core::texture::{SpectrumType, TextureEvalContext, TextureMapping2D};
 use crate::spectra::{
@@ -7,6 +7,7 @@ use crate::spectra::{
     SampledWavelengths,
 };
 use crate::utils::Ptr;
+use crate::Float;
 
 /* GPU heavy code, dont know if this will ever work the way Im doing things.
 * Leaving it here isolated, for careful handling */
@@ -15,6 +16,7 @@ use crate::utils::Ptr;
 pub struct GPUSpectrumImageTexture {
     pub mapping: TextureMapping2D,
     pub tex_obj: u64,
+    pub image: Ptr<Image>,
     pub scale: Float,
     pub invert: bool,
     pub is_single_channel: bool,

shared/src/textures/mod.rs

@@ -6,7 +6,6 @@ pub mod fbm;
 pub mod image;
 pub mod marble;
 pub mod mix;
-pub mod ptex;
 pub mod scaled;
 pub mod windy;
 pub mod wrinkled;
@@ -19,7 +18,6 @@ pub use fbm::*;
 pub use image::*;
 pub use marble::*;
 pub use mix::*;
-pub use ptex::*;
 pub use scaled::*;
 pub use windy::*;
 pub use wrinkled::*;

shared/src/textures/ptex.rs

@@ -1,51 +0,0 @@
-use crate::Float;
-use crate::core::color::{ColorEncoding, RGB};
-use crate::core::spectrum::{SpectrumTrait, StandardSpectra};
-use crate::core::texture::{SpectrumType, TextureEvalContext};
-use crate::spectra::{
-    DeviceStandardColorSpaces, RGBAlbedoSpectrum, RGBColorSpace, RGBIlluminantSpectrum,
-    RGBUnboundedSpectrum, SampledSpectrum, SampledWavelengths,
-};
-use crate::utils::Ptr;
-
-#[repr(C)]
-#[derive(Debug, Clone, Copy)]
-pub struct GPUFloatPtexTexture {
-    pub face_values: *const Float,
-}
-
-impl GPUFloatPtexTexture {
-    pub fn evaluate(&self, ctx: &TextureEvalContext) -> Float {
-        unsafe { *self.face_values.add(ctx.face_index as usize) }
-    }
-}
-
-#[repr(C)]
-#[derive(Clone, Debug, Copy)]
-pub struct GPUSpectrumPtexTexture {
-    pub face_values: Ptr<RGB>,
-    pub n_faces: u32,
-    pub spectrum_type: SpectrumType,
-    pub colorspaces: DeviceStandardColorSpaces,
-}
-
-impl GPUSpectrumPtexTexture {
-    pub fn evaluate(
-        &self,
-        ctx: &TextureEvalContext,
-        lambda: &SampledWavelengths,
-    ) -> SampledSpectrum {
-        let index = (ctx.face_index as u32).clamp(0, self.n_faces.saturating_sub(1));
-        let rgb = unsafe { &*self.face_values.add(index as usize) };
-        let s_rgb = self.colorspaces.srgb;
-
-        match self.spectrum_type {
-            SpectrumType::Unbounded => RGBUnboundedSpectrum::new(&s_rgb, *rgb).sample(lambda),
-            SpectrumType::Albedo => {
-                let clamped_rgb = rgb.clamp(0.0, 1.0);
-                RGBAlbedoSpectrum::new(&s_rgb, clamped_rgb).sample(lambda)
-            }
-            SpectrumType::Illuminant => RGBIlluminantSpectrum::new(&s_rgb, *rgb).sample(lambda),
-        }
-    }
-}

shared/src/utils/math.rs

@@ -71,15 +71,16 @@ pub fn evaluate_polynomial(t: Float, coeffs: &[Float]) -> Float {
     result
 }
 
-pub fn difference_of_products<T>(a: Float, b: T, c: Float, d: T) -> T
+pub fn difference_of_products<T, U, V>(a: T, b: U, c: T, d: U) -> V
 where
-    T: Copy + Neg<Output = T> + Mul<Float, Output = T> + Add<Output = T>,
-    T: MulAdd<Float, T, Output = T>,
+    T: Copy + Neg<Output = T>,
+    U: Copy + MulAdd<T, V, Output = V>,
+    V: Copy + Zero + Neg<Output = V> + Add<Output = V>,
 {
-    let cd = d * c;
-    let diff = b.mul_add(a, -cd);
-    let error = d.mul_add(-c, cd);
-    diff + error
+    let cd: V = d.mul_add(c, V::zero());
+    let diff: V = b.mul_add(a, -cd);
+    let err:  V = d.mul_add(-c, cd);
+    diff + err
 }
 
 #[inline]

shared/src/utils/transform.rs

@@ -94,46 +94,44 @@ impl<T: NumFloat + Sum + Product> Default for TransformGeneric<T> {
     }
 }
 
+impl<T: NumFloat> TransformGeneric<T> {
+    #[inline]
+    pub fn apply_to_point(&self, p: Point<T, 3>) -> Point<T, 3> {
+        let x = self.m[0][0]*p.x() + self.m[0][1]*p.y() + self.m[0][2]*p.z() + self.m[0][3];
+        let y = self.m[1][0]*p.x() + self.m[1][1]*p.y() + self.m[1][2]*p.z() + self.m[1][3];
+        let z = self.m[2][0]*p.x() + self.m[2][1]*p.y() + self.m[2][2]*p.z() + self.m[2][3];
+        let w = self.m[3][0]*p.x() + self.m[3][1]*p.y() + self.m[3][2]*p.z() + self.m[3][3];
+        if w == T::one() {
+            Point([x, y, z])
+        } else {
+            Point([x / w, y / w, z / w])
+        }
+    }
+
+    #[inline]
+    pub fn apply_to_vector(&self, v: Vector<T, 3>) -> Vector<T, 3> {
+        // directions: linear part only, no translation, no w-divide
+        Vector([
+            self.m[0][0]*v.x() + self.m[0][1]*v.y() + self.m[0][2]*v.z(),
+            self.m[1][0]*v.x() + self.m[1][1]*v.y() + self.m[1][2]*v.z(),
+            self.m[2][0]*v.x() + self.m[2][1]*v.y() + self.m[2][2]*v.z(),
+        ])
+    }
+
+    #[inline]
+    pub fn apply_to_normal(&self, n: Normal<T, 3>) -> Normal<T, 3> {
+        // normals: inverse-transpose, no translation, no w-divide
+        Normal([
+            self.m_inv[0][0]*n.x() + self.m_inv[1][0]*n.y() + self.m_inv[2][0]*n.z(),
+            self.m_inv[0][1]*n.x() + self.m_inv[1][1]*n.y() + self.m_inv[2][1]*n.z(),
+            self.m_inv[0][2]*n.x() + self.m_inv[1][2]*n.y() + self.m_inv[2][2]*n.z(),
+        ])
+    }
+}
+
 pub type Transform = TransformGeneric<Float>;
 
-impl TransformGeneric<Float> {
-    pub fn apply_to_point(&self, p: Point3f) -> Point3f {
-        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 + self.m[0][3];
-        let yp = self.m[1][0] * x + self.m[1][1] * y + self.m[1][2] * z + self.m[1][3];
-        let zp = self.m[2][0] * x + self.m[2][1] * y + self.m[2][2] * z + self.m[2][3];
-        let wp = self.m[3][0] * x + self.m[3][1] * y + self.m[3][2] * z + self.m[3][3];
-
-        if wp == 1. {
-            Point3f::new(xp, yp, zp)
-        } else {
-            Point3f::new(xp / wp, yp / wp, zp / wp)
-        }
-    }
-
-    pub fn apply_to_vector(&self, v: Vector3f) -> Vector3f {
-        let x = v.x();
-        let y = v.y();
-        let z = v.z();
-        let xv = self.m[0][0] * x + self.m[0][1] * y + self.m[0][2] * z;
-        let yv = self.m[1][0] * x + self.m[1][1] * y + self.m[1][2] * z;
-        let zv = self.m[2][0] * x + self.m[2][1] * y + self.m[2][2] * z;
-        Vector3f::new(xv, yv, zv)
-
-    }
-
-    pub fn apply_to_normal(&self, p: Normal3f) -> Normal3f {
-        let x = p.x();
-        let y = p.y();
-        let z = p.z();
-        let xn = self.m_inv[0][0] * x + self.m_inv[1][1] * y + self.m_inv[2][0] * z;
-        let yn = self.m_inv[0][1] * x + self.m_inv[1][1] * y + self.m_inv[2][1] * z;
-        let zn = self.m_inv[0][2] * x + self.m_inv[1][2] * y + self.m_inv[2][2] * z;
-        Normal3f::new(xn, yn, zn)
-    }
-
+impl Transform {
     pub fn apply_to_bounds(&self, b: Bounds3f) -> Bounds3f {
         let mut new_bounds = Bounds3f::default();
         for i in 0..8 {
@@ -669,61 +667,49 @@ impl Mul<TransformGeneric<Float>> for Float {
     }
 }
 
-impl<T> Mul<Point<T, 3>> for TransformGeneric<T>
-where
-    T: NumFloat,
-{
+impl<T: NumFloat> Mul<Point<T, 3>> for TransformGeneric<T> {
     type Output = Point<T, 3>;
-    fn mul(self, p: Point<T, 3>) -> Self::Output {
-        let xp = self.m[0][0] * p.x() + self.m[0][1] * p.y() + self.m[0][2] * p.z() + self.m[0][3];
-        let yp = self.m[1][0] * p.x() + self.m[1][1] * p.y() + self.m[1][2] * p.z() + self.m[1][3];
-        let zp = self.m[2][0] * p.x() + self.m[2][1] * p.y() + self.m[2][2] * p.z() + self.m[2][3];
-        let wp = self.m[3][0] * p.x() + self.m[3][1] * p.y() + self.m[3][2] * p.z() + self.m[3][3];
-
-        if wp == T::one() {
-            Point([xp, yp, zp])
-        } else {
-            Point([xp / wp, yp / wp, zp / wp])
+    #[inline]
+    fn mul(self, p: Point<T, 3>) -> Point<T, 3> {
+        self.apply_to_point(p)
     }
 }
-}
-
-impl<T> Mul<Vector<T, 3>> for TransformGeneric<T>
-where
-    T: NumFloat,
-{
+impl<T: NumFloat> Mul<Vector<T, 3>> for TransformGeneric<T> {
     type Output = Vector<T, 3>;
-    fn mul(self, p: Vector<T, 3>) -> Self::Output {
-        let xp = self.m[0][0] * p.x() + self.m[0][1] * p.y() + self.m[0][2] * p.z() + self.m[0][3];
-        let yp = self.m[1][0] * p.x() + self.m[1][1] * p.y() + self.m[1][2] * p.z() + self.m[1][3];
-        let zp = self.m[2][0] * p.x() + self.m[2][1] * p.y() + self.m[2][2] * p.z() + self.m[2][3];
-        let wp = self.m[3][0] * p.x() + self.m[3][1] * p.y() + self.m[3][2] * p.z() + self.m[3][3];
-
-        if wp == T::one() {
-            Vector([xp, yp, zp])
-        } else {
-            Vector([xp / wp, yp / wp, zp / wp])
+    #[inline]
+    fn mul(self, v: Vector<T, 3>) -> Vector<T, 3> {
+        self.apply_to_vector(v)
     }
 }
-}
-
-impl<T> Mul<Normal<T, 3>> for TransformGeneric<T>
-where
-    T: NumFloat,
-{
+impl<T: NumFloat> Mul<Normal<T, 3>> for TransformGeneric<T> {
     type Output = Normal<T, 3>;
-    fn mul(self, p: Normal<T, 3>) -> Self::Output {
-        let xp = self.m[0][0] * p.x() + self.m[0][1] * p.y() + self.m[0][2] * p.z() + self.m[0][3];
-        let yp = self.m[1][0] * p.x() + self.m[1][1] * p.y() + self.m[1][2] * p.z() + self.m[1][3];
-        let zp = self.m[2][0] * p.x() + self.m[2][1] * p.y() + self.m[2][2] * p.z() + self.m[2][3];
-        let wp = self.m[3][0] * p.x() + self.m[3][1] * p.y() + self.m[3][2] * p.z() + self.m[3][3];
-
-        if wp == T::one() {
-            Normal([xp, yp, zp])
-        } else {
-            Normal([xp / wp, yp / wp, zp / wp])
+    #[inline]
+    fn mul(self, n: Normal<T, 3>) -> Normal<T, 3> {
+        self.apply_to_normal(n)
     }
 }
+
+// reference versions, so `&t * v` works without consuming the transform
+impl<T: NumFloat> Mul<Point<T, 3>> for &TransformGeneric<T> {
+    type Output = Point<T, 3>;
+    #[inline]
+    fn mul(self, p: Point<T, 3>) -> Point<T, 3> {
+        self.apply_to_point(p)
+    }
+}
+impl<T: NumFloat> Mul<Vector<T, 3>> for &TransformGeneric<T> {
+    type Output = Vector<T, 3>;
+    #[inline]
+    fn mul(self, v: Vector<T, 3>) -> Vector<T, 3> {
+        self.apply_to_vector(v)
+    }
+}
+impl<T: NumFloat> Mul<Normal<T, 3>> for &TransformGeneric<T> {
+    type Output = Normal<T, 3>;
+    #[inline]
+    fn mul(self, n: Normal<T, 3>) -> Normal<T, 3> {
+        self.apply_to_normal(n)
+    }
 }
 
 impl From<Frame> for TransformGeneric<Float> {

shared/src/wavefront/workitems.rs

@@ -1,4 +1,5 @@
 use crate::core::bxdf::BxDFFlags;
+use crate::core::film::VisibleSurface;
 use crate::core::geometry::{
     Normal3f, Point2f, Point2i, Point3f, Point3fi, Ray, RayDifferential, Vector3f,
 };
@@ -15,22 +16,13 @@ use crate::{Float, Ptr};
 #[repr(C)]
 #[derive(Clone, Copy)]
 pub struct PixelSampleState {
-    pub filter_weight: SoABuffer<Float>,
-    pub p_film: SoABuffer<Point2f>,
+    pub pixel: SoABuffer<Point2i>,
     pub l: SoABuffer<SampledSpectrum>,
     pub lambda: SoABuffer<SampledWavelengths>,
-    pub r_u: SoABuffer<SampledSpectrum>,
-    pub r_l: SoABuffer<SampledSpectrum>,
-    pub prev_intr_ctx: SoABuffer<LightSampleContext>,
-    pub beta: SoABuffer<SampledSpectrum>,
-    pub depth: SoABuffer<u32>,
-    pub specular_bounce: SoABuffer<u8>,
-    pub any_non_specular_bounces: SoABuffer<u8>,
-    pub eta_scale: SoABuffer<Float>,
+    pub filter_weight: SoABuffer<Float>,
     pub camera_ray_weight: SoABuffer<SampledSpectrum>,
-    pub visible_surface_idx: SoABuffer<u32>,
+    pub visible_surface: SoABuffer<VisibleSurface>,
     pub samples: SoABuffer<RaySamples>,
-    pub p_pixel: SoABuffer<Point2i>,
 }
 
 impl SoA for PixelSampleState {
@@ -38,22 +30,13 @@ impl SoA for PixelSampleState {
 
     fn allocate(n: u32, alloc: &dyn SoAAllocator) -> Self {
         Self {
-            filter_weight: alloc_soa_buffer(n, alloc),
-            p_film: alloc_soa_buffer(n, alloc),
+            pixel: alloc_soa_buffer(n, alloc),
             l: alloc_soa_buffer(n, alloc),
             lambda: alloc_soa_buffer(n, alloc),
-            r_u: alloc_soa_buffer(n, alloc),
-            r_l: alloc_soa_buffer(n, alloc),
-            prev_intr_ctx: alloc_soa_buffer(n, alloc),
-            beta: alloc_soa_buffer(n, alloc),
-            depth: alloc_soa_buffer(n, alloc),
-            specular_bounce: alloc_soa_buffer(n, alloc),
-            any_non_specular_bounces: alloc_soa_buffer(n, alloc),
-            eta_scale: alloc_soa_buffer(n, alloc),
+            filter_weight: alloc_soa_buffer(n, alloc),
             camera_ray_weight: alloc_soa_buffer(n, alloc),
-            visible_surface_idx: alloc_soa_buffer(n, alloc),
+            visible_surface: alloc_soa_buffer(n, alloc),
             samples: alloc_soa_buffer(n, alloc),
-            p_pixel: alloc_soa_buffer(n, alloc),
         }
     }
 
@@ -73,8 +56,8 @@ pub struct RayWorkItem {
     pub r_l: SampledSpectrum,
     pub prev_intr_ctx: LightSampleContext,
     pub eta_scale: Float,
-    pub specular_bounce: u8,
-    pub any_non_specular_bounces: u8,
+    pub specular_bounce: bool,
+    pub any_non_specular_bounces: bool,
 }
 
 #[repr(C)]
@@ -89,8 +72,8 @@ pub struct RayWorkItemSoA {
     pub r_l: SoABuffer<SampledSpectrum>,
     pub prev_intr_ctx: SoABuffer<LightSampleContext>,
     pub eta_scale: SoABuffer<Float>,
-    pub specular_bounce: SoABuffer<u8>,
-    pub any_non_specular_bounces: SoABuffer<u8>,
+    pub specular_bounce: SoABuffer<bool>,
+    pub any_non_specular_bounces: SoABuffer<bool>,
 }
 
 impl SoA for RayWorkItemSoA {
@@ -332,14 +315,13 @@ pub struct MaterialEvalWorkItem {
     pub lambda: SampledWavelengths,
     pub beta: SampledSpectrum,
     pub r_u: SampledSpectrum,
-    // pub r_l: SampledSpectrum,
     pub any_non_specular_bounces: bool,
     pub depth: u32,
     pub eta_scale: Float,
     pub dpdus: Vector3f,
     pub dpdvs: Vector3f,
     pub dndus: Normal3f,
-    pub dndvs: Normal3f
+    pub dndvs: Normal3f,
 }
 
 #[repr(C)]
@@ -361,15 +343,13 @@ pub struct MaterialEvalWorkItemSoA {
     pub lambda: SoABuffer<SampledWavelengths>,
     pub beta: SoABuffer<SampledSpectrum>,
     pub r_u: SoABuffer<SampledSpectrum>,
-    // pub r_l: SoABuffer<SampledSpectrum>,
     pub any_non_specular_bounces: SoABuffer<u8>,
     pub depth: SoABuffer<u32>,
     pub eta_scale: SoABuffer<Float>,
     pub dpdus: SoABuffer<Vector3f>,
     pub dpdvs: SoABuffer<Vector3f>,
     pub dndus: SoABuffer<Normal3f>,
-    pub dndvs: SoABuffer<Normal3f>
-
+    pub dndvs: SoABuffer<Normal3f>,
 }
 
 impl SoA for MaterialEvalWorkItemSoA {
@@ -393,7 +373,6 @@ impl SoA for MaterialEvalWorkItemSoA {
             lambda: alloc_soa_buffer(n, alloc),
             beta: alloc_soa_buffer(n, alloc),
             r_u: alloc_soa_buffer(n, alloc),
-            // r_l: alloc_soa_buffer(n, alloc),
             any_non_specular_bounces: alloc_soa_buffer(n, alloc),
             depth: alloc_soa_buffer(n, alloc),
             eta_scale: alloc_soa_buffer(n, alloc),
@@ -422,7 +401,6 @@ impl SoA for MaterialEvalWorkItemSoA {
             lambda: self.lambda.get(i),
             beta: self.beta.get(i),
             r_u: self.r_u.get(i),
-            // r_l: self.r_l.get(i),
             any_non_specular_bounces: self.any_non_specular_bounces.get(i) != 0,
             depth: self.depth.get(i),
             eta_scale: self.eta_scale.get(i),
@@ -430,7 +408,6 @@ impl SoA for MaterialEvalWorkItemSoA {
             dpdvs: self.dpdvs.get(i),
             dndus: self.dndus.get(i),
             dndvs: self.dndvs.get(i),
-
         }
     }
 
@@ -451,7 +428,6 @@ impl SoA for MaterialEvalWorkItemSoA {
         self.lambda.set(i, v.lambda);
         self.beta.set(i, v.beta);
         self.r_u.set(i, v.r_u);
-        // self.r_l.set(i, v.r_l);
         self.any_non_specular_bounces
             .set(i, v.any_non_specular_bounces as u8);
         self.depth.set(i, v.depth);

src/core/film.rs

@@ -15,7 +15,7 @@ use shared::spectra::{
     cie::SWATCHES_RAW, DenselySampledSpectrum, PiecewiseLinearSpectrum, RGBColorSpace,
 };
 use shared::utils::math::{linear_least_squares, SquareMatrix};
-use shared::{Float, Ptr, leak};
+use shared::{leak, Float, Ptr};
 use std::sync::atomic::{AtomicUsize, Ordering};
 use std::sync::{Arc, LazyLock};
 
@@ -100,7 +100,7 @@ impl CreatePixelSensor for PixelSensor {
                 output_colorspace.as_ref(),
                 sensor_illum.as_deref(),
                 imaging_ratio,
-                arena
+                arena,
             ))
         } else {
             let r_opt = get_named_spectrum(&format!("{}_r", sensor_name));
@@ -129,7 +129,7 @@ impl CreatePixelSensor for PixelSensor {
                         .expect("Sensor must have illuminant"),
                 ),
                 imaging_ratio,
-                arena
+                arena,
             ))
         }
     }
@@ -205,7 +205,7 @@ impl CreatePixelSensor for PixelSensor {
         output_colorspace: &RGBColorSpace,
         sensor_illum: Option<&Spectrum>,
         imaging_ratio: Float,
-        arena: &Arena
+        arena: &Arena,
     ) -> Self {
         let spectra = get_spectra_context();
         let r_bar = CIE_X_DATA.clone();
@@ -229,7 +229,6 @@ impl CreatePixelSensor for PixelSensor {
             imaging_ratio,
         }
     }
-
 }
 
 pub trait CreateFilmBase {

src/textures/image.rs

@@ -3,16 +3,16 @@ use crate::core::texture::{
     CreateFloatTexture, CreateSpectrumTexture, FloatTexture, FloatTextureTrait, SpectrumTexture,
     SpectrumTextureTrait,
 };
-use crate::utils::mipmap::{MIPMap, MIPMapFilterOptions};
-use crate::utils::{FileLoc, TextureParameterDictionary};
-use crate::Arena;
+use crate::utils::mipmap::{FilterFunction, MIPMap, MIPMapFilterOptions};
+use crate::utils::{FileLoc, TextureParameterDictionary, resolve_filename};
+use crate::{Arena};
 use anyhow::Result;
 use shared::core::color::RGB;
 use shared::core::color::{ColorEncoding, SRGBEncoding};
 use shared::core::geometry::Vector2f;
 use shared::core::image::WrapMode;
 use shared::core::spectrum::SpectrumTrait;
-use shared::core::texture::{SpectrumType, TextureEvalContext, TextureMapping2D};
+use shared::core::texture::{SpectrumType, TexCoord2D, TextureEvalContext, TextureMapping2D};
 use shared::spectra::{
     RGBAlbedoSpectrum, RGBIlluminantSpectrum, RGBUnboundedSpectrum, SampledSpectrum,
     SampledWavelengths,
@@ -168,6 +168,7 @@ impl CreateSpectrumTexture for SpectrumImageTexture {
 
         let filter_options = MIPMapFilterOptions::default();
         let wrap_str = parameters.get_one_string("wrap", "repeat")?;
+        // let wrap_mode = WrapMode::parse(wrap_str)?;
         let wrap_mode = match wrap_str.as_str() {
             "repeat" => WrapMode::Repeat,
             "clamp" => WrapMode::Clamp,
@@ -222,18 +223,65 @@ impl FloatImageTexture {
 }
 
 impl FloatTextureTrait for FloatImageTexture {
-    fn evaluate(&self, _ctx: &TextureEvalContext) -> Float {
-        todo!()
+    fn evaluate(&self, ctx: &TextureEvalContext) -> Float {
+        let mut c: TexCoord2D = self.base.mapping.map(ctx);
+        c.st[1] = 1. - c.st[1];
+        let v: Float = self.base.scale
+            * self.base.mipmap.filter::<Float>(
+                c.st,
+                Vector2f::new(c.dsdx, c.dtdx),
+                Vector2f::new(c.dsdy, c.dtdy),
+            );
+
+        if self.base.invert {
+            (1. - v).max(0.)
+        } else {
+            v
+        }
     }
 }
 
 impl CreateFloatTexture for FloatImageTexture {
     fn create(
-        _render_from_texture: Transform,
-        _parameters: TextureParameterDictionary,
-        _loc: FileLoc,
+        render_from_texture: Transform,
+        parameters: TextureParameterDictionary,
+        loc: FileLoc,
         _arena: &Arena,
     ) -> Result<FloatTexture> {
-        todo!()
+        let mapping = TextureMapping2D::create(&parameters, &render_from_texture, &loc)?;
+        let max_aniso = parameters.get_one_float("maxanisotropy", 8.)?;
+        let filter = parameters.get_one_string("filter", "bilinear")?;
+        let mut filter_options = MIPMapFilterOptions::default();
+        filter_options.max_anisotropy = max_aniso;
+
+        let ff = FilterFunction::parse(&filter)?;
+        filter_options.filter = ff;
+        let wrap_string = parameters.get_one_string("wrap", "repeat")?;
+        let wrap_mode = WrapMode::parse(&wrap_string)?;
+        let scale = parameters.get_one_float("scale", 1.)?;
+        let invert = parameters.get_one_bool("invert", false)?;
+        let filename = resolve_filename(&parameters.get_one_string("filename", "")?);
+        let default_encoding = if Path::new(&filename)
+            .extension()
+            .map_or(false, |ext| ext == "png")
+        {
+            "sRGB"
+        } else {
+            "linear"
+        };
+        let encoding_str = parameters.get_one_string("encoding", default_encoding)?;
+        let encoding = ColorEncoding::from_name(&encoding_str)?;
+
+        let tex = FloatImageTexture::new(
+            mapping,
+            filename,
+            filter_options,
+            wrap_mode,
+            scale,
+            invert,
+            encoding,
+        );
+
+        Ok(FloatTexture::Image(tex))
     }
 }

src/utils/mipmap.rs

@@ -3,6 +3,7 @@ use shared::core::color::{ColorEncoding, RGB};
 use shared::core::geometry::{Point2f, Point2i, Vector2f, VectorLike};
 use shared::core::image::{WrapMode, WrapMode2D};
 use shared::spectra::RGBColorSpace;
+use anyhow::{bail, Result};
 use shared::utils::math::{lerp, safe_sqrt, square};
 use shared::Float;
 use std::hash::{Hash, Hasher};
@@ -21,6 +22,19 @@ pub enum FilterFunction {
     Ewa,
 }
 
+impl FilterFunction {
+    pub fn parse(name: &str) -> Result<FilterFunction> {
+        match name {
+            "ewa" | "EWA" => Ok(FilterFunction::Ewa),
+            "trilinear" => Ok(FilterFunction::Trilinear),
+            "bilinear" => Ok(FilterFunction::Bilinear),
+            "point" => Ok(FilterFunction::Point),
+            _ => bail!("Filter function unknown")
+        }
+    }
+}
+
+
 impl std::fmt::Display for FilterFunction {
     fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
         let s = match self {

src/wavefront/aggregate.rs

@@ -1,4 +1,5 @@
 use crate::globals::get_options;
+use log::debug;
 use rayon::prelude::*;
 use shared::core::geometry::{Bounds3f, Ray, VectorLike};
 use shared::core::interaction::{InteractionTrait, SurfaceInteraction};
@@ -60,7 +61,7 @@ impl WavefrontAggregate for CpuAggregate {
                     r_u: r.r_u,
                     r_l: r.r_l,
                     depth: r.depth,
-                    specular_bounce: r.specular_bounce != 0,
+                    specular_bounce: r.specular_bounce,
                     prev_intr_ctx: r.prev_intr_ctx,
                 });
                 return;
@@ -83,14 +84,14 @@ impl WavefrontAggregate for CpuAggregate {
                     p: intr.p(),
                     n: intr.n(),
                     uv: intr.common.uv,
-                    wo: -r.ray.d,
+                    wo: intr.wo(),
                     lambda: r.lambda,
                     pixel_index: r.pixel_index,
                     beta: r.beta,
                     r_u: r.r_u,
                     r_l: r.r_l,
                     depth: r.depth,
-                    specular_bounce: r.specular_bounce != 0,
+                    specular_bounce: r.specular_bounce,
                     prev_intr_ctx: r.prev_intr_ctx,
                 });
             }
@@ -103,12 +104,20 @@ impl WavefrontAggregate for CpuAggregate {
                 universal_eval_mtl_q
             };
 
+            if material.is_conductor() {
+                debug!("shading frame: {:?}", intr.shading.dpdu);
+                debug!(
+                    "dot product: {:?}",
+                    intr.shading.dpdu.normalize().dot(intr.n().into())
+                );
+            }
+
             eval_q.push(MaterialEvalWorkItem {
                 p: intr.pi(),
                 n: intr.n(),
                 ns: intr.shading.n,
-                dpdu: intr.shading.dpdu,
-                dpdv: intr.shading.dpdv,
+                dpdu: intr.dpdu,
+                dpdv: intr.dpdv,
                 uv: intr.common.uv,
                 wo: intr.wo(),
                 time: r.ray.time,
@@ -120,7 +129,7 @@ impl WavefrontAggregate for CpuAggregate {
                 lambda: r.lambda,
                 beta: r.beta,
                 r_u: r.r_u,
-                any_non_specular_bounces: r.any_non_specular_bounces != 0,
+                any_non_specular_bounces: r.any_non_specular_bounces,
                 depth: r.depth,
                 eta_scale: r.eta_scale,
                 dpdus: intr.shading.dpdu,

src/wavefront/integrator.rs

@@ -75,8 +75,6 @@ where
         let scanlines_per_pass = (max_samples / res_x).max(1);
         let max_queue_size = res_x * scanlines_per_pass;
 
-        eprintln!("wavefront got {} lights", lights.len());
-
         let mut infinite_lights = gvec();
         for light in &lights {
             if light.light_type().is_infinite() {
@@ -84,10 +82,6 @@ where
             }
         }
 
-        eprintln!("infinite_lights len = {}", infinite_lights.len());
-
-        // for light in
-
         let cpu_aggregate = CpuAggregate::new(*aggregate);
 
         CpuWavefrontRenderer(WavefrontPathIntegrator {
@@ -241,7 +235,7 @@ impl CpuWavefrontRenderer {
                     pixel_bounds.p_min.x() + (pixel_index as i32 % x_resolution),
                     y0 + (pixel_index as i32 / x_resolution),
                 );
-                pixel_sample_state.p_pixel.set(pixel_index, p_pixel);
+                pixel_sample_state.pixel.set(pixel_index, p_pixel);
 
                 // Skipped pixels contribute nothing; their slots are simply never
                 // populated, and update_film filters them by the same bounds test.
@@ -263,9 +257,6 @@ impl CpuWavefrontRenderer {
                 pixel_sample_state
                     .filter_weight
                     .set(pixel_index, camera_sample.filter_weight);
-                pixel_sample_state
-                    .p_film
-                    .set(pixel_index, camera_sample.p_film);
 
                 let Some(camera_ray) = camera.generate_ray(camera_sample, &lambda) else {
                     pixel_sample_state
@@ -288,8 +279,8 @@ impl CpuWavefrontRenderer {
                     r_l: SampledSpectrum::new(1.0),
                     prev_intr_ctx: LightSampleContext::default(),
                     eta_scale: 1.0,
-                    specular_bounce: 0,
-                    any_non_specular_bounces: 0,
+                    specular_bounce: false,
+                    any_non_specular_bounces: false,
                 });
             });
     }
@@ -318,7 +309,7 @@ impl CpuWavefrontRenderer {
                 if w.depth == 0 || w.specular_bounce {
                     l_contrib += w.beta * le / w.r_u.average();
                 } else {
-                    // MIS: combine BSDF and light sampling weights via ratio tracking
+                    // Compute MIS-weighted radiance contribution from infinite light
                     let ctx = w.prev_intr_ctx;
                     let light_choice_pdf = light_sampler.pmf_with_context(&ctx, light);
                     let r_l = w.r_l * light_choice_pdf * light.pdf_li(&ctx, w.ray_d, true);
@@ -353,13 +344,14 @@ impl CpuWavefrontRenderer {
             let l_contrib = if w.depth == 0 || w.specular_bounce {
                 w.beta * le / w.r_u.average()
             } else {
+                let wi = -w.wo;
                 let ctx = w.prev_intr_ctx;
                 let light_choice_pdf = light_sampler.pmf_with_context(&ctx, light);
                 // wi from previous interaction to this light hit
-                let wi = (w.p - Point3f::from(ctx.pi)).normalize();
                 let light_pdf = light_choice_pdf * light.pdf_li(&ctx, wi, true);
+                let r_u = w.r_u;
                 let r_l = w.r_l * light_pdf;
-                w.beta * le / (w.r_u + r_l).average()
+                w.beta * le / (r_u + r_l).average()
             };
 
             if !l_contrib.is_black() {
@@ -382,6 +374,9 @@ impl CpuWavefrontRenderer {
         } else {
             &self.basic_eval_material_queue
         };
+        if self.regularize {
+            eprintln!("regularize=true");
+        }
 
         let n = queue.size();
         let next = ((depth + 1) % 2) as usize;
@@ -459,8 +454,14 @@ impl CpuWavefrontRenderer {
                     eta_scale *= square(bs.eta);
                 }
 
+                if material.is_conductor() {
+                    if bs.is_specular() {
+                        eprintln!("NON SPECULAR");
+                    }
+                }
+
                 let rr_beta = (beta * eta_scale / r_u.average()).max_component_value();
-                if rr_beta < 1.0 && w.depth > 1 {
+                if rr_beta < 1.0 && w.depth >= 1 {
                     let q = (1.0 - rr_beta).max(0.0_f32);
                     if rs.indirect.rr < q {
                         beta = SampledSpectrum::new(0.0);
@@ -489,8 +490,8 @@ impl CpuWavefrontRenderer {
                         r_l,
                         prev_intr_ctx: ctx,
                         eta_scale,
-                        specular_bounce: bs.is_specular() as u8,
-                        any_non_specular_bounces: any_non_specular as u8,
+                        specular_bounce: bs.is_specular(),
+                        any_non_specular_bounces: any_non_specular,
                     });
                 }
             }
@@ -505,7 +506,7 @@ impl CpuWavefrontRenderer {
                 };
 
                 if flags.is_reflective() && !flags.is_transmissive() {
-                    light_ctx.pi = Point3fi::new_from_point(Ray::offset_origin(&w.p, &w.n, &(-wo)));
+                    light_ctx.pi = Point3fi::new_from_point(Ray::offset_origin(&w.p, &w.n, &wo));
                 } else if flags.is_transmissive() && flags.is_reflective() {
                     light_ctx.pi = Point3fi::new_from_point(Ray::offset_origin(&w.p, &w.n, &(-wo)));
                 }
@@ -573,7 +574,7 @@ impl CpuWavefrontRenderer {
         (0..self.max_queue_size as usize)
             .into_par_iter()
             .for_each(|pixel_index| {
-                let p_pixel = self.pixel_sample_state.p_pixel.get(pixel_index);
+                let p_pixel = self.pixel_sample_state.pixel.get(pixel_index);
                 if !pixel_bounds.contains_exclusive(p_pixel) {
                     return;
                 }
@@ -600,7 +601,7 @@ impl CpuWavefrontRenderer {
             let w = unsafe { ray_queue.storage.get(i) };
             let dimension = 6 + 7 * w.depth;
             let pi = w.pixel_index as usize;
-            let p_pixel = pixel_sample_state.p_pixel.get(pi);
+            let p_pixel = pixel_sample_state.pixel.get(pi);
 
             let mut sampler = sampler_proto.clone();
             sampler.start_pixel_sample(p_pixel, sample_index as i32, Some(dimension));