Updated light management and creation, refactored module. Started work on rendering and integrators, fixed up maths functions

src/core/options.rs

@@ -0,0 +1,89 @@
+use std::ops::Deref;
+use std::sync::OnceLock;
+
+#[derive(Debug, Clone, Copy, PartialEq, Eq)]
+pub enum RenderingCoordinateSystem {
+    Camera,
+    CameraWorld,
+    World,
+}
+
+#[derive(Debug, Clone, Copy)]
+pub struct BasicPBRTOptions {
+    pub seed: i32,
+    pub quiet: bool,
+    pub disable_pixel_jitter: bool,
+    pub disable_wavelength_jitter: bool,
+    pub disable_texture_filtering: bool,
+    pub force_diffuse: bool,
+    pub use_gpu: bool,
+    pub wavefront: bool,
+    pub rendering_space: RenderingCoordinateSystem,
+}
+
+impl Default for BasicPBRTOptions {
+    fn default() -> Self {
+        Self {
+            seed: 0,
+            quiet: false,
+            disable_pixel_jitter: false,
+            disable_wavelength_jitter: false,
+            disable_texture_filtering: false,
+            force_diffuse: false,
+            use_gpu: false,
+            wavefront: false,
+            rendering_space: RenderingCoordinateSystem::CameraWorld,
+        }
+    }
+}
+
+#[derive(Debug, Clone)]
+pub struct PBRTOptions {
+    pub basic: BasicPBRTOptions,
+
+    pub n_threads: usize,
+    pub log_level: String,
+    pub image_file: String,
+}
+
+impl Default for PBRTOptions {
+    fn default() -> Self {
+        Self {
+            basic: BasicPBRTOptions::default(),
+            n_threads: 0, // 0 usually implies "autodetect"
+            log_level: "info".to_string(),
+            image_file: "output.exr".to_string(),
+        }
+    }
+}
+
+impl Deref for PBRTOptions {
+    type Target = BasicPBRTOptions;
+
+    fn deref(&self) -> &Self::Target {
+        &self.basic
+    }
+}
+
+// -----------------------------------------------------------------------
+// Global State Management
+// -----------------------------------------------------------------------
+
+static OPTIONS: OnceLock<PBRTOptions> = OnceLock::new();
+
+pub fn init_pbrt(options: PBRTOptions) {
+    OPTIONS
+        .set(options)
+        .expect("PBRT has already been initialized!");
+}
+
+pub fn cleanup_pbrt() {
+    todo!()
+}
+
+pub fn get_options() -> &'static PBRTOptions {
+    OPTIONS.get().unwrap_or_else(|| {
+        static DEFAULT: OnceLock<PBRTOptions> = OnceLock::new();
+        DEFAULT.get_or_init(PBRTOptions::default)
+    })
+}

src/integrators/mod.rs

@@ -0,0 +1,520 @@
+mod pipeline;
+
+use pipeline::*;
+
+use crate::camera::{Camera, CameraTrait};
+use crate::core::bxdf::{BSDF, BxDFFlags, BxDFTrait, FArgs, TransportMode};
+use crate::core::film::{FilmTrait, VisibleSurface};
+use crate::core::interaction::{Interaction, InteractionTrait, SurfaceInteraction};
+use crate::core::options::get_options;
+use crate::core::pbrt::{Float, SHADOW_EPSILON};
+use crate::core::primitive::{Primitive, PrimitiveTrait};
+use crate::core::sampler::{CameraSample, Sampler, SamplerTrait};
+use crate::geometry::{Bounds2i, Point2f, Point2i, Point3fi, Ray, Vector3f, VectorLike};
+use crate::lights::sampler::{LightSampler, UniformLightSampler};
+use crate::lights::{Light, LightSampleContext, LightTrait, sampler::LightSamplerTrait};
+use crate::shapes::ShapeIntersection;
+use crate::spectra::{SampledSpectrum, SampledWavelengths};
+use crate::utils::math::square;
+use crate::utils::sampling::{
+    power_heuristic, sample_uniform_hemisphere, sample_uniform_sphere, uniform_hemisphere_pdf,
+    uniform_sphere_pdf,
+};
+use bumpalo::Bump;
+use rayon::prelude::*;
+
+use std::sync::Arc;
+
+#[derive(Clone, Debug)]
+pub struct IntegratorBase {
+    pub aggregate: Arc<Primitive>,
+    pub lights: Vec<Arc<Light>>,
+    pub infinite_lights: Vec<Arc<Light>>,
+}
+
+impl IntegratorBase {
+    pub fn new(aggregate: Arc<Primitive>, lights: Vec<Arc<Light>>) -> Self {
+        let infinite_lights = lights
+            .iter()
+            .filter(|light| light.light_type().is_infinite())
+            .cloned()
+            .collect();
+
+        Self {
+            aggregate,
+            lights,
+            infinite_lights,
+        }
+    }
+
+    pub fn intersect(&self, ray: &Ray, t_max: Option<Float>) -> Option<ShapeIntersection> {
+        self.aggregate.intersect(ray, t_max)
+    }
+
+    pub fn intersect_p(&self, ray: &Ray, t_max: Option<Float>) -> bool {
+        self.aggregate.intersect_p(ray, t_max)
+    }
+
+    pub fn unoccluded(&self, p0: &Interaction, p1: &Interaction) -> bool {
+        !self.intersect_p(&p0.spawn_ray_to_interaction(p1), Some(1. - SHADOW_EPSILON))
+    }
+}
+
+pub trait IntegratorTrait {
+    fn render(&self);
+}
+
+pub trait RayIntegratorTrait: Send + Sync + std::fmt::Debug {
+    fn evaluate_pixel_sample(
+        &self,
+        p_pixel: Point2i,
+        sample_ind: usize,
+        sampler: &mut Sampler,
+        scratch: &Bump,
+    );
+
+    fn li(
+        &self,
+        ray: Ray,
+        lambda: &SampledWavelengths,
+        sampler: &mut Sampler,
+        scratch: &Bump,
+        visible_surface: bool,
+    ) -> (SampledSpectrum, Option<VisibleSurface>);
+}
+
+struct NextRaySample {
+    wi: Vector3f,
+    f: SampledSpectrum,
+    pdf: Float,
+    is_specular: bool,
+}
+
+#[derive(Clone, Debug)]
+pub struct SimplePathIntegrator {
+    base: IntegratorBase,
+    camera: Arc<Camera>,
+    light_sampler: UniformLightSampler,
+    sample_lights: bool,
+    sample_bsdf: bool,
+    max_depth: usize,
+}
+
+impl SimplePathIntegrator {
+    fn sample_direct_lighting(
+        &self,
+        isect: &SurfaceInteraction,
+        bsdf: &BSDF,
+        sampler: &mut Sampler,
+        lambda: &SampledWavelengths,
+    ) -> SampledSpectrum {
+        let Some(sampled_light) = self.light_sampler.sample(sampler.get1d()) else {
+            return SampledSpectrum::new(0.0);
+        };
+        let u_light = sampler.get2d();
+        let p0 = LightSampleContext::from(isect);
+
+        let Some(ls) = sampled_light.light.sample_li(&p0, u_light, lambda, false) else {
+            return SampledSpectrum::new(0.0);
+        };
+
+        if ls.pdf == 0.0 || ls.l.is_black() {
+            return SampledSpectrum::new(0.0);
+        }
+
+        let wi = ls.wi;
+        let wo = -isect.wo();
+
+        let Some(mut f) = bsdf.f(wo, wi, TransportMode::Radiance) else {
+            return SampledSpectrum::new(0.);
+        };
+
+        f *= wi.abs_dot(isect.shading.n.into());
+
+        if f.is_black() {
+            return SampledSpectrum::new(0.0);
+        }
+
+        // Visibility check
+        let p0_surf = Interaction::Surface(isect.clone());
+        if !self.base.unoccluded(&p0_surf, ls.p_light.as_ref()) {
+            return SampledSpectrum::new(0.0);
+        }
+
+        f * ls.l / (sampled_light.p * ls.pdf)
+    }
+
+    fn sample_bsdf_direction(
+        &self,
+        wo: Vector3f,
+        bsdf: &BSDF,
+        sampler: &mut Sampler,
+        _isect: &SurfaceInteraction,
+    ) -> Option<NextRaySample> {
+        let u = sampler.get1d();
+        let bs = bsdf.sample_f(wo, u, sampler.get2d(), FArgs::default())?;
+
+        Some(NextRaySample {
+            wi: bs.wi,
+            f: bs.f,
+            pdf: bs.pdf,
+            is_specular: bs.is_specular(),
+        })
+    }
+
+    fn sample_uniform_direction(
+        &self,
+        wo: Vector3f,
+        bsdf: &BSDF,
+        sampler: &mut Sampler,
+        isect: &SurfaceInteraction,
+    ) -> Option<NextRaySample> {
+        let flags = bsdf.flags();
+        let mut wi;
+        let pdf;
+
+        if BxDFFlags::is_reflective(&flags) && BxDFFlags::is_transmissive(&flags) {
+            wi = sample_uniform_sphere(sampler.get2d());
+            pdf = uniform_sphere_pdf();
+        } else {
+            wi = sample_uniform_hemisphere(sampler.get2d());
+            pdf = uniform_hemisphere_pdf();
+
+            let same_hemi = wo.dot(isect.n().into()) * wi.dot(isect.n().into()) > 0.0;
+            if BxDFFlags::is_reflective(&flags) && !same_hemi
+                || BxDFFlags::is_transmissive(&flags) && same_hemi
+            {
+                wi = -wi;
+            }
+        }
+
+        let f = bsdf.f(wo, wi, TransportMode::Radiance)?;
+
+        Some(NextRaySample {
+            wi,
+            f,
+            pdf,
+            is_specular: false,
+        })
+    }
+}
+
+impl RayIntegratorTrait for SimplePathIntegrator {
+    fn evaluate_pixel_sample(
+        &self,
+        p_pixel: Point2i,
+        sample_ind: usize,
+        sampler: &mut Sampler,
+        scratch: &Bump,
+    ) {
+        pipeline::evaluate_pixel_sample(
+            self,
+            self.camera.as_ref(),
+            sampler,
+            p_pixel,
+            sample_ind,
+            scratch,
+        );
+    }
+
+    fn li(
+        &self,
+        mut ray: Ray,
+        lambda: &SampledWavelengths,
+        sampler: &mut Sampler,
+        scratch: &Bump,
+        _visible_surface: bool,
+    ) -> (SampledSpectrum, Option<VisibleSurface>) {
+        let mut l = SampledSpectrum::new(0.0);
+        let mut beta = SampledSpectrum::new(1.0);
+        let mut specular_bounce = true;
+        let mut depth = 0;
+        while !beta.is_black() {
+            let Some(mut si) = self.base.intersect(&ray, None) else {
+                if !self.sample_lights || specular_bounce {
+                    for light in &self.base.infinite_lights {
+                        l += light.le(&ray, lambda);
+                    }
+                }
+                break;
+            };
+
+            let t_hit = si.t_hit();
+            let isect = &mut si.intr;
+            if !self.sample_lights || specular_bounce {
+                l += beta * isect.le(-ray.d, lambda);
+            }
+            depth += 1;
+            if depth == self.max_depth {
+                break;
+            }
+
+            let Some(bsdf) = isect.get_bsdf(&ray, lambda, self.camera.as_ref(), scratch, sampler)
+            else {
+                specular_bounce = false;
+                isect.skip_intersection(&mut ray, t_hit);
+                continue;
+            };
+
+            let wo = -ray.d;
+            if self.sample_lights {
+                l += beta * self.sample_direct_lighting(isect, &bsdf, sampler, lambda);
+            }
+
+            let sample = if self.sample_bsdf {
+                self.sample_bsdf_direction(wo, &bsdf, sampler, isect)
+            } else {
+                self.sample_uniform_direction(wo, &bsdf, sampler, isect)
+            };
+
+            let Some(bs) = sample else {
+                break;
+            };
+
+            beta *= bs.f * bs.wi.abs_dot(isect.shading.n.into()) / bs.pdf;
+            specular_bounce = bs.is_specular;
+            ray = isect.spawn_ray(bs.wi);
+        }
+        assert!(beta.y(lambda) > 0.);
+        debug_assert!(beta.y(lambda).is_finite());
+        (l, None)
+    }
+}
+
+#[derive(Clone, Debug)]
+pub struct PathIntegrator {
+    base: IntegratorBase,
+    camera: Arc<Camera>,
+    sampler_prototype: Sampler,
+    light_sampler: LightSampler,
+    regularize: bool,
+    max_depth: usize,
+}
+
+impl PathIntegrator {
+    fn sample_ld(
+        &self,
+        intr: &SurfaceInteraction,
+        bsdf: &BSDF,
+        lambda: &SampledWavelengths,
+        sampler: &mut Sampler,
+    ) -> Option<SampledSpectrum> {
+        let mut ctx: LightSampleContext = intr.into();
+        let flags = bsdf.flags();
+        if BxDFFlags::is_reflective(&flags) ^ !BxDFFlags::is_transmissive(&flags) {
+            ctx.pi = Point3fi::new_from_point(intr.offset_ray_vector(-intr.wo()));
+        }
+
+        let u = sampler.get1d();
+        let sampled_light = self.light_sampler.sample_with_context(&ctx, u)?;
+
+        let u_light = sampler.get2d();
+        let light = sampled_light.light;
+        let ls = light.sample_li(&ctx, u_light, lambda, true)?;
+
+        if ls.l.is_black() || ls.pdf == 0.0 {
+            return None;
+        }
+
+        let wo = intr.wo();
+        let wi = ls.wi;
+        let f = bsdf.f(wo, wi, TransportMode::Radiance)? * wi.abs_dot(intr.shading.n.into());
+
+        if !self
+            .base
+            .unoccluded(&Interaction::Surface(intr.clone()), &ls.p_light)
+        {
+            return None;
+        }
+
+        let pl = sampled_light.p * ls.pdf;
+        if light.light_type().is_delta_light() {
+            Some(ls.l * f / pl)
+        } else {
+            let pb = bsdf.pdf(wo, wi, FArgs::default());
+            let wl = power_heuristic(1, pl, 1, pb);
+            Some(wl * ls.l * f / pl)
+        }
+    }
+}
+
+impl RayIntegratorTrait for PathIntegrator {
+    fn evaluate_pixel_sample(
+        &self,
+        p_pixel: Point2i,
+        sample_ind: usize,
+        sampler: &mut Sampler,
+        scratch: &Bump,
+    ) {
+        pipeline::evaluate_pixel_sample(
+            self,
+            self.camera.as_ref(),
+            sampler,
+            p_pixel,
+            sample_ind,
+            scratch,
+        );
+    }
+
+    fn li(
+        &self,
+        mut ray: Ray,
+        lambda: &SampledWavelengths,
+        sampler: &mut Sampler,
+        scratch: &Bump,
+        visible_surface: bool,
+    ) -> (SampledSpectrum, Option<VisibleSurface>) {
+        let mut l = SampledSpectrum::new(0.0);
+        let mut beta = SampledSpectrum::new(1.0);
+        let mut visible: Option<VisibleSurface> = None;
+        let mut depth = 0;
+        let mut pb = 1.;
+        let mut eta_scale = 1.;
+        let mut specular_bounce = false;
+        let mut any_non_specular_bounces = false;
+
+        let mut prev_int_ctx = LightSampleContext::default();
+        loop {
+            let Some(mut si) = self.base.intersect(&ray, None) else {
+                for light in &self.base.infinite_lights {
+                    let le = light.le(&ray, lambda);
+                    if depth == 0 || specular_bounce {
+                        l += beta * le;
+                    } else {
+                        let pl = self.light_sampler.pmf_with_context(&prev_int_ctx, light)
+                            * light.pdf_li(&prev_int_ctx, ray.d, true);
+                        let w_b = power_heuristic(1, pb, 1, pl);
+                        l += beta * w_b * le;
+                    }
+                }
+                break;
+            };
+
+            let t_hit = si.t_hit();
+            let isect = &mut si.intr;
+            let le = isect.le(-ray.d, lambda);
+            if depth == 0 || specular_bounce {
+                l += beta * le;
+            } else if let Some(area_light) = &isect.area_light {
+                let pl = self
+                    .light_sampler
+                    .pmf_with_context(&prev_int_ctx, area_light)
+                    + area_light.pdf_li(&prev_int_ctx, ray.d, true);
+                let wl = power_heuristic(1, pb, 1, pl);
+                l += beta * wl * le;
+            }
+
+            let Some(mut bsdf) =
+                isect.get_bsdf(&ray, lambda, self.camera.as_ref(), scratch, sampler)
+            else {
+                specular_bounce = true;
+                isect.skip_intersection(&mut ray, t_hit);
+                continue;
+            };
+
+            if depth == 0 && visible_surface {
+                const N_RHO_SAMPLES: usize = 16;
+                const UC_RHO: [Float; N_RHO_SAMPLES] = [
+                    0.75741637,
+                    0.37870818,
+                    0.7083487,
+                    0.18935409,
+                    0.9149363,
+                    0.35417435,
+                    0.5990858,
+                    0.09467703,
+                    0.8578725,
+                    0.45746812,
+                    0.686759,
+                    0.17708716,
+                    0.9674518,
+                    0.2995429,
+                    0.5083201,
+                    0.047338516,
+                ];
+                let u_rho: [Point2f; N_RHO_SAMPLES] = [
+                    Point2f::new(0.855985, 0.570367),
+                    Point2f::new(0.381823, 0.851844),
+                    Point2f::new(0.285328, 0.764262),
+                    Point2f::new(0.733380, 0.114073),
+                    Point2f::new(0.542663, 0.344465),
+                    Point2f::new(0.127274, 0.414848),
+                    Point2f::new(0.964700, 0.947162),
+                    Point2f::new(0.594089, 0.643463),
+                    Point2f::new(0.095109, 0.170369),
+                    Point2f::new(0.825444, 0.263359),
+                    Point2f::new(0.429467, 0.454469),
+                    Point2f::new(0.244460, 0.816459),
+                    Point2f::new(0.756135, 0.731258),
+                    Point2f::new(0.516165, 0.152852),
+                    Point2f::new(0.180888, 0.214174),
+                    Point2f::new(0.898579, 0.503897),
+                ];
+
+                let albedo = bsdf.rho_wo(isect.wo(), &UC_RHO, &u_rho);
+                visible = Some(VisibleSurface::new(isect, &albedo, lambda));
+            }
+
+            if self.regularize && any_non_specular_bounces {
+                bsdf.regularize();
+            }
+
+            depth += 1;
+
+            if depth == self.max_depth {
+                break;
+            }
+
+            if BxDFFlags::is_non_specular(&bsdf.flags()) {
+                if let Some(ld) = self.sample_ld(isect, &bsdf, lambda, sampler) {
+                    l += beta * ld;
+                }
+            }
+
+            let wo = -ray.d;
+            let u = sampler.get1d();
+            let Some(bs) = bsdf.sample_f(wo, u, sampler.get2d(), FArgs::default()) else {
+                break;
+            };
+
+            beta *= bs.f * bs.wi.abs_dot(isect.shading.n.into()) / bs.pdf;
+            pb = if bs.pdf_is_proportional {
+                bsdf.pdf(wo, bs.wi, FArgs::default())
+            } else {
+                bs.pdf
+            };
+            specular_bounce = bs.is_specular();
+            any_non_specular_bounces |= !bs.is_specular();
+            if bs.is_transmissive() {
+                eta_scale *= square(bs.eta);
+            }
+            prev_int_ctx = LightSampleContext::from(&si.intr);
+            ray = si
+                .intr
+                .spawn_ray_with_differentials(&ray, bs.wi, bs.flags, bs.eta);
+
+            let rr_beta = beta * eta_scale;
+            if rr_beta.max_component_value() < 1. && depth > 1 {
+                let q = (1. - rr_beta.max_component_value()).max(0.);
+                if sampler.get1d() < q {
+                    break;
+                }
+                beta /= 1. - q;
+            }
+        }
+
+        (l, visible)
+    }
+}
+
+pub enum Integrator {
+    BPDT(BPDTIntegrator),
+    MLT(MLTIntegrator),
+    SPPM(SPPMIntegrator),
+    Sampler(SamplerIntegrator),
+}
+
+pub struct BPDTIntegrator;
+pub struct MLTIntegrator;
+pub struct SPPMIntegrator;
+pub struct SamplerIntegrator;

src/integrators/pipeline.rs

@@ -0,0 +1,264 @@
+use crate::core::{options::PBRTOptions, sampler::get_camera_sample};
+use crate::image::{Image, ImageMetadata};
+use indicatif::{ProgressBar, ProgressStyle};
+use std::io::Write;
+use std::path::Path;
+
+use super::*;
+
+struct PbrtProgress {
+    bar: ProgressBar,
+}
+
+impl PbrtProgress {
+    fn new(total_work: u64, description: &str, quiet: bool) -> Self {
+        if quiet {
+            return Self {
+                bar: ProgressBar::hidden(),
+            };
+        }
+
+        let bar = ProgressBar::new(total_work);
+
+        bar.set_style(
+            ProgressStyle::default_bar()
+                .template("[{elapsed_precise}] {bar:40.cyan/blue} {pos:>7}/{len:7} {msg}")
+                .unwrap()
+                .progress_chars("=>-"),
+        );
+
+        bar.set_message(description.to_string());
+
+        Self { bar }
+    }
+
+    fn update(&self, amount: u64) {
+        self.bar.inc(amount);
+    }
+
+    fn done(&self) {
+        self.bar.finish_with_message("Done");
+    }
+
+    fn elapsed_seconds(&self) -> f32 {
+        self.bar.elapsed().as_secs_f32()
+    }
+}
+
+fn generate_tiles(bounds: Bounds2i) -> Vec<Bounds2i> {
+    let mut tiles = Vec::new();
+    const TILE_SIZE: i32 = 16;
+    for y in (bounds.p_min.y()..bounds.p_max.y()).step_by(TILE_SIZE as usize) {
+        for x in (bounds.p_min.x()..bounds.p_max.x()).step_by(TILE_SIZE as usize) {
+            let p_min = Point2i::new(x, y);
+            let p_max = Point2i::new(
+                (x + TILE_SIZE).min(bounds.p_max.x()),
+                (y + TILE_SIZE).min(bounds.p_max.y()),
+            );
+            tiles.push(Bounds2i::from_points(p_min, p_max));
+        }
+    }
+    tiles
+}
+
+pub fn render<T>(
+    integrator: &T,
+    _base: &IntegratorBase,
+    camera: &Camera,
+    sampler_prototype: &Sampler,
+) where
+    T: RayIntegratorTrait,
+{
+    let options = get_options();
+    if let Some((p_pixel, sample_index)) = options.debug_start {
+        let s_index = sample_index as usize;
+        let scratch = Bump::new();
+        let mut tile_sampler = sampler_prototype.clone();
+
+        tile_sampler.start_pixel_sample(p_pixel, s_index, None);
+
+        evaluate_pixel_sample(
+            integrator,
+            camera,
+            &mut tile_sampler,
+            p_pixel,
+            s_index,
+            &scratch,
+        );
+        return;
+    }
+
+    let pixel_bounds = camera.get_film().pixel_bounds();
+    let spp = sampler_prototype.samples_per_pixel();
+    let total_work = (pixel_bounds.area() as u64) * (spp as u64);
+    let progress = PbrtProgress::new(total_work, "Rendering", options.quiet);
+    let mut wave_start = 0;
+    let mut wave_end = 1;
+    let mut next_wave_size = 1;
+
+    let mut reference_image: Option<Image> = None;
+    let mut mse_out_file: Option<std::fs::File> = None;
+
+    if let Some(ref_path) = &options.mse_reference_image {
+        let image_and_metadata =
+            Image::read(Path::new(&ref_path), None).expect("Could not load image");
+        let image = image_and_metadata.image;
+        let metadata = image_and_metadata.metadata;
+        let resolution = image.resolution();
+        // reference_image = Some(image);
+        let mse_pixel_bounds = metadata
+            .pixel_bounds
+            .unwrap_or_else(|| Bounds2i::from_points(Point2i::new(0, 0), resolution));
+
+        if !mse_pixel_bounds.overlaps(&pixel_bounds) {
+            panic!("Output pixel bounds dont fit inside the reference image");
+        }
+
+        let crop_p_min = Point2i::from(pixel_bounds.p_min - mse_pixel_bounds.p_min);
+        let crop_p_max = Point2i::from(pixel_bounds.p_max - mse_pixel_bounds.p_min);
+
+        let crop_bounds = Bounds2i::from_points(crop_p_min, crop_p_max);
+
+        let cropped_image = image.crop(crop_bounds).clone();
+        let cropped_resolution = cropped_image.resolution();
+
+        let expected_res = Point2i::new(
+            pixel_bounds.p_max.x() - pixel_bounds.p_min.x(),
+            pixel_bounds.p_max.y() - pixel_bounds.p_min.y(),
+        );
+
+        reference_image = Some(cropped_image);
+
+        assert_eq!(
+            cropped_resolution, expected_res,
+            "Cropped reference image resolution mismatch"
+        );
+
+        if let Some(out_path) = &options.mse_reference_output {
+            mse_out_file = Some(
+                std::fs::File::create(out_path)
+                    .expect(&format!("Failed to create MSE output file: {}", out_path)),
+            );
+        }
+    }
+
+    let tiles = generate_tiles(pixel_bounds);
+    while wave_start < spp {
+        tiles.par_iter().for_each(|tile_bounds| {
+            let mut arena = Bump::with_capacity(65 * 1024);
+            let mut sampler = sampler_prototype.clone();
+
+            for p_pixel in tile_bounds {
+                for sample_index in wave_start..wave_end {
+                    sampler.start_pixel_sample(*p_pixel, sample_index, None);
+
+                    evaluate_pixel_sample(
+                        integrator,
+                        camera,
+                        &mut sampler,
+                        *p_pixel,
+                        sample_index,
+                        &arena,
+                    );
+
+                    arena.reset();
+                }
+            }
+
+            let work_done = (tile_bounds.area() as u64) * ((wave_end - wave_start) as u64);
+            progress.update(work_done);
+        });
+
+        wave_start = wave_end;
+        wave_end = (wave_end + next_wave_size).min(spp);
+
+        if reference_image.is_none() {
+            next_wave_size = (2 * next_wave_size).min(64);
+        }
+
+        if wave_start == spp {
+            progress.done();
+        }
+
+        if wave_start == spp || options.write_partial_images || reference_image.is_some() {
+            let mut metadata = ImageMetadata {
+                render_time_seconds: Some(progress.elapsed_seconds()),
+                samples_per_pixel: Some(wave_start as i32),
+                ..Default::default()
+            };
+
+            if wave_start == spp || options.write_partial_images {
+                camera.init_metadata(&mut metadata);
+                camera
+                    .get_film()
+                    .write_image(&metadata, 1.0 / wave_start as Float);
+            }
+
+            if let Some(ref_img) = &reference_image {
+                let splat_scale = 1.0 / (wave_start as Float);
+
+                let film_metadata = ImageMetadata::default();
+                let film_image = camera.get_film().get_image(&film_metadata, splat_scale);
+
+                let (mse_values, _mse_debug_img) =
+                    film_image.mse(film_image.all_channels_desc(), ref_img, false);
+                let mse_avg = mse_values.average();
+
+                if let Some(file) = &mut mse_out_file {
+                    writeln!(file, "{}, {:.9}", wave_start, mse_avg).ok();
+                    file.flush().ok();
+                }
+
+                metadata.mse = Some(mse_avg);
+            }
+        }
+    }
+}
+
+pub fn evaluate_pixel_sample<T: RayIntegratorTrait>(
+    integrator: &T,
+    camera: &Camera,
+    sampler: &mut Sampler,
+    pixel: Point2i,
+    _sample_index: usize,
+    scratch: &Bump,
+) {
+    let mut lu = sampler.get1d();
+    if get_options().disable_wavelength_jitter {
+        lu = 0.5;
+    }
+    let lambda = camera.get_film().sample_wavelengths(lu);
+    let mut film = camera.get_film();
+    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() < 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);
+        }
+
+        let initialize_visible_surface = film.uses_visible_surface();
+        let (mut l, visible_surface) = integrator.li(
+            camera_ray.ray,
+            &lambda,
+            sampler,
+            scratch,
+            initialize_visible_surface,
+        );
+        l *= camera_ray.weight;
+
+        if l.has_nans() || l.y(&lambda).is_infinite() {
+            l = SampledSpectrum::new(0.);
+        }
+
+        film.add_sample(
+            pixel,
+            l,
+            &lambda,
+            visible_surface.as_ref(),
+            camera_sample.filter_weight,
+        );
+    }
+}

src/lights/diffuse.rs

@@ -1,19 +1,255 @@
+use super::{
+    DenselySampledSpectrum, LightBase, LightBounds, LightLiSample, LightSampleContext, LightTrait,
+    LightType, RGBIlluminantSpectrum, SampledSpectrum, SampledWavelengths, Spectrum, SpectrumTrait,
+};
+use crate::core::interaction::{
+    Interaction, InteractionTrait, SimpleInteraction, SurfaceInteraction,
+};
 use crate::core::medium::MediumInterface;
-use crate::core::pbrt::Float;
+use crate::core::pbrt::{Float, PI};
-use crate::shapes::Shape;
+use crate::core::texture::{
-use crate::spectra::Spectrum;
+    FloatTexture, FloatTextureTrait, TextureEvalContext, TextureEvaluator,
+    UniversalTextureEvaluator,
+};
+use crate::geometry::{
+    Bounds3f, Normal3f, Point2f, Point2fi, Point2i, Point3f, Point3fi, Ray, Vector3f, VectorLike,
+};
+use crate::image::Image;
+use crate::shapes::{Shape, ShapeSample, ShapeSampleContext, ShapeTrait};
+use crate::utils::color::RGB;
+use crate::utils::colorspace::RGBColorSpace;
+use crate::utils::hash::hash_float;
 use crate::utils::transform::Transform;
 
 use std::sync::Arc;
 
+#[derive(Clone, Debug)]
 pub struct DiffuseAreaLight {
-    pub l_emit: Spectrum,
+    base: LightBase,
-    pub shape: Arc<Shape>,
+    shape: Shape,
-    pub two_sided: bool,
+    alpha: Option<FloatTexture>,
-    pub area: Float,
+    area: Float,
-    pub flags: u8,
+    two_sided: bool,
-    pub n_samples: i32,
+    lemit: Arc<DenselySampledSpectrum>,
-    pub medium_interface: MediumInterface,
+    scale: Float,
-    light_to_world: Transform<Float>,
+    image: Option<Image>,
-    world_to_light: Transform<Float>,
+    image_color_space: Option<Arc<RGBColorSpace>>,
+}
+
+impl DiffuseAreaLight {
+    #[allow(clippy::too_many_arguments)]
+    pub fn new(
+        render_from_light: Transform<Float>,
+        medium_interface: MediumInterface,
+        le: Spectrum,
+        scale: Float,
+        shape: Shape,
+        alpha: FloatTexture,
+        image: Option<Image>,
+        image_color_space: Option<Arc<RGBColorSpace>>,
+        two_sided: bool,
+    ) -> Self {
+        let is_constant_zero = match &alpha {
+            FloatTexture::FloatConstant(tex) => tex.evaluate(&TextureEvalContext::default()) == 0.0,
+            _ => false,
+        };
+
+        let (light_type, stored_alpha) = if is_constant_zero {
+            (LightType::DeltaPosition, None)
+        } else {
+            (LightType::Area, Some(alpha))
+        };
+
+        let base = LightBase::new(light_type, &render_from_light, &medium_interface);
+
+        let lemit = LightBase::lookup_spectrum(&le);
+        if let Some(im) = &image {
+            let desc = im
+                .get_channel_desc(&["R", "G", "B"])
+                .expect("Image used for DiffuseAreaLight doesn't have R, G, B channels");
+
+            assert_eq!(3, desc.size(), "Image channel description size mismatch");
+            assert!(
+                desc.is_identity(),
+                "Image channel description is not identity"
+            );
+
+            assert!(
+                image_color_space.is_some(),
+                "Image provided but ColorSpace is missing"
+            );
+        }
+        let is_triangle_or_bilinear = matches!(shape, Shape::Triangle(_) | Shape::BilinearPatch(_));
+        if render_from_light.has_scale(None) && !is_triangle_or_bilinear {
+            println!(
+                "Scaling detected in rendering to light space transformation! \
+                 The system has numerous assumptions, implicit and explicit, \
+                 that this transform will have no scale factors in it. \
+                 Proceed at your own risk; your image may have errors."
+            );
+        }
+
+        Self {
+            base,
+            area: shape.area(),
+            shape,
+            alpha: stored_alpha,
+            two_sided,
+            lemit,
+            scale,
+            image,
+            image_color_space,
+        }
+    }
+
+    fn alpha_masked(&self, intr: &Interaction) -> bool {
+        let Some(alpha_tex) = &self.alpha else {
+            return false;
+        };
+        let ctx = TextureEvalContext::from(intr);
+        let a = UniversalTextureEvaluator.evaluate_float(alpha_tex, &ctx);
+        if a >= 1.0 {
+            return false;
+        }
+        if a <= 0.0 {
+            return true;
+        }
+        hash_float(&intr.p()) > a
+    }
+}
+
+impl LightTrait for DiffuseAreaLight {
+    fn base(&self) -> &LightBase {
+        &self.base
+    }
+
+    fn phi(&self, lambda: SampledWavelengths) -> SampledSpectrum {
+        let mut l = SampledSpectrum::new(0.);
+        if let Some(image) = &self.image {
+            for y in 0..image.resolution().y() {
+                for x in 0..image.resolution().x() {
+                    let mut rgb = RGB::default();
+                    for c in 0..3 {
+                        rgb[c] = image.get_channel(Point2i::new(x, y), c);
+                    }
+                    l += RGBIlluminantSpectrum::new(
+                        self.image_color_space.as_ref().unwrap(),
+                        RGB::clamp_zero(rgb),
+                    )
+                    .sample(&lambda);
+                }
+            }
+            l *= self.scale / (image.resolution().x() * image.resolution().y()) as Float;
+        } else {
+            l = self.lemit.sample(&lambda) * self.scale;
+        }
+        let two_side = if self.two_sided { 2. } else { 1. };
+        PI * two_side * self.area * l
+    }
+
+    fn sample_li(
+        &self,
+        ctx: &LightSampleContext,
+        u: Point2f,
+        lambda: &SampledWavelengths,
+        _allow_incomplete_pdf: bool,
+    ) -> Option<LightLiSample> {
+        let shape_ctx = ShapeSampleContext::new(ctx.pi, ctx.n, ctx.ns, 0.0);
+        let ss = self.shape.sample_from_context(&shape_ctx, u)?;
+        let mut intr: SurfaceInteraction = ss.intr.as_ref().clone();
+
+        intr.common.medium_interface = Some(self.base.medium_interface.clone());
+        let p = intr.p();
+        let n = intr.n();
+        let uv = intr.uv;
+
+        let generic_intr = Interaction::Surface(intr);
+        if self.alpha_masked(&generic_intr) {
+            return None;
+        }
+
+        let wi = (p - ctx.p()).normalize();
+        let le = self.l(p, n, uv, -wi, &lambda);
+
+        if le.is_black() {
+            return None;
+        }
+
+        Some(LightLiSample::new(le, wi, ss.pdf, generic_intr))
+    }
+
+    fn pdf_li(&self, ctx: &LightSampleContext, wi: Vector3f, _allow_incomplete_pdf: bool) -> Float {
+        let shape_ctx = ShapeSampleContext::new(ctx.pi, ctx.n, ctx.ns, 0.);
+        self.shape.pdf_from_context(&shape_ctx, wi)
+    }
+
+    fn l(
+        &self,
+        p: Point3f,
+        n: Normal3f,
+        mut uv: Point2f,
+        w: Vector3f,
+        lambda: &SampledWavelengths,
+    ) -> SampledSpectrum {
+        if self.two_sided && n.dot(w.into()) < 0. {
+            return SampledSpectrum::new(0.);
+        }
+        let intr = Interaction::Surface(SurfaceInteraction::new_minimal(
+            Point3fi::new_from_point(p),
+            uv,
+        ));
+
+        if self.alpha_masked(&intr) {
+            return SampledSpectrum::new(0.);
+        }
+        if let Some(image) = &self.image {
+            let mut rgb = RGB::default();
+            uv[1] = 1. - uv[1];
+            for c in 0..3 {
+                rgb[c] = image.bilerp_channel(uv, c);
+            }
+
+            let spec = RGBIlluminantSpectrum::new(
+                self.image_color_space.as_ref().unwrap(),
+                RGB::clamp_zero(rgb),
+            );
+
+            self.scale * spec.sample(lambda)
+        } else {
+            self.scale * self.lemit.sample(lambda)
+        }
+    }
+
+    fn le(&self, _ray: &Ray, _lambda: &SampledWavelengths) -> SampledSpectrum {
+        todo!()
+    }
+
+    fn preprocess(&mut self, _scene_bounds: &Bounds3f) {
+        unimplemented!()
+    }
+
+    fn bounds(&self) -> Option<LightBounds> {
+        let mut phi = 0.;
+        if let Some(image) = &self.image {
+            for y in 0..image.resolution.y() {
+                for x in 0..image.resolution.x() {
+                    for c in 0..3 {
+                        phi += image.get_channel(Point2i::new(x, y), c);
+                    }
+                }
+            }
+        } else {
+            phi = self.lemit.max_value();
+        }
+
+        let nb = self.shape.normal_bounds();
+        Some(LightBounds::new(
+            &self.shape.bounds(),
+            nb.w,
+            phi,
+            nb.cos_theta,
+            (PI / 2.).cos(),
+            self.two_sided,
+        ))
+    }
 }

src/lights/infinite.rs

@@ -0,0 +1,593 @@
+use crate::{
+    core::medium::Medium,
+    geometry::Frame,
+    spectra::RGBIlluminantSpectrum,
+    utils::{
+        color::RGB,
+        colorspace::RGBColorSpace,
+        math::equal_area_square_to_sphere,
+        sampling::{
+            AliasTable, PiecewiseConstant2D, WindowedPiecewiseConstant2D, sample_uniform_sphere,
+            uniform_sphere_pdf,
+        },
+    },
+};
+
+use rayon::prelude::*;
+
+use crate::core::pbrt::clamp_t;
+use crate::image::{PixelFormat, WrapMode};
+
+use std::path::Path;
+
+use super::{
+    Arc, Bounds2f, Bounds3f, DenselySampledSpectrum, Float, Image, Interaction, LightBase,
+    LightBounds, LightLiSample, LightSampleContext, LightTrait, LightType, MediumInterface,
+    Normal3f, PI, Point2f, Point2i, Point3f, Ray, SampledSpectrum, SampledWavelengths,
+    SimpleInteraction, Spectrum, Transform, Vector3f, VectorLike, equal_area_sphere_to_square,
+    square,
+};
+
+#[derive(Debug, Clone)]
+pub struct InfiniteUniformLight {
+    base: LightBase,
+    lemit: Arc<DenselySampledSpectrum>,
+    scale: Float,
+    scene_center: Point3f,
+    scene_radius: Float,
+}
+
+impl InfiniteUniformLight {
+    pub fn new(render_from_light: Transform<Float>, le: Spectrum, scale: Float) -> Self {
+        let base = LightBase::new(
+            LightType::Infinite,
+            &render_from_light,
+            &MediumInterface::default(),
+        );
+        let lemit = LightBase::lookup_spectrum(&le);
+        Self {
+            base,
+            lemit,
+            scale,
+            scene_center: Point3f::default(),
+            scene_radius: 0.,
+        }
+    }
+}
+
+impl LightTrait for InfiniteUniformLight {
+    fn base(&self) -> &LightBase {
+        &self.base
+    }
+    fn sample_li(
+        &self,
+        ctx: &LightSampleContext,
+        u: Point2f,
+        lambda: &SampledWavelengths,
+        allow_incomplete_pdf: bool,
+    ) -> Option<LightLiSample> {
+        if allow_incomplete_pdf {
+            return None;
+        }
+        let wi = sample_uniform_sphere(u);
+        let pdf = uniform_sphere_pdf();
+        let intr_simple = SimpleInteraction::new_interface(
+            ctx.p() + wi * (2. * self.scene_radius),
+            Some(MediumInterface::default()),
+        );
+
+        let intr = Interaction::Simple(intr_simple);
+        Some(LightLiSample::new(
+            self.scale * self.lemit.sample(&lambda),
+            wi,
+            pdf,
+            intr,
+        ))
+    }
+
+    fn phi(&self, lambda: SampledWavelengths) -> SampledSpectrum {
+        4. * PI * PI * square(self.scene_radius) * self.scale * self.lemit.sample(&lambda)
+    }
+
+    fn pdf_li(
+        &self,
+        _ctx: &LightSampleContext,
+        _wi: Vector3f,
+        allow_incomplete_pdf: bool,
+    ) -> Float {
+        if allow_incomplete_pdf {
+            return 0.;
+        }
+        uniform_sphere_pdf()
+    }
+
+    fn l(
+        &self,
+        _p: Point3f,
+        _n: Normal3f,
+        _uv: Point2f,
+        _w: Vector3f,
+        _lambda: &SampledWavelengths,
+    ) -> SampledSpectrum {
+        todo!()
+    }
+
+    fn le(&self, _ray: &Ray, lambda: &SampledWavelengths) -> SampledSpectrum {
+        self.scale * self.lemit.sample(lambda)
+    }
+    fn preprocess(&mut self, _scene_bounds: &Bounds3f) {
+        todo!()
+    }
+    fn bounds(&self) -> Option<LightBounds> {
+        todo!()
+    }
+}
+
+#[derive(Clone, Debug)]
+pub struct InfiniteImageLight {
+    base: LightBase,
+    image: Image,
+    image_color_space: Arc<RGBColorSpace>,
+    scale: Float,
+    scene_radius: Float,
+    scene_center: Point3f,
+    distrib: PiecewiseConstant2D,
+    compensated_distrib: PiecewiseConstant2D,
+}
+
+impl InfiniteImageLight {
+    pub fn new(
+        render_from_light: Transform<Float>,
+        image: Image,
+        image_color_space: Arc<RGBColorSpace>,
+        scale: Float,
+        filename: String,
+    ) -> Self {
+        let base = LightBase::new(
+            LightType::Infinite,
+            &render_from_light,
+            &MediumInterface::default(),
+        );
+
+        let desc = image
+            .get_channel_desc(&["R", "G", "B"])
+            .expect("Image used for DiffuseAreaLight doesn't have R, G, B channels");
+
+        assert_eq!(3, desc.size());
+        assert!(desc.is_identity());
+        if image.resolution().x() != image.resolution().y() {
+            panic!(
+                "{}: image resolution ({}, {}) is non-square. It's unlikely this is an equal area environment map.",
+                filename,
+                image.resolution.x(),
+                image.resolution.y()
+            );
+        }
+        let mut d = image.get_sampling_distribution_uniform();
+        let domain = Bounds2f::from_points(Point2f::new(0., 0.), Point2f::new(1., 1.));
+        let distrib = PiecewiseConstant2D::new_with_bounds(&d, domain);
+        let slice = &mut d.values; // or d.as_slice_mut()
+        let count = slice.len() as Float;
+        let sum: Float = slice.iter().sum();
+        let average = sum / count;
+
+        for v in slice.iter_mut() {
+            *v = (*v - average).max(0.0);
+        }
+
+        let all_zero = slice.iter().all(|&v| v == 0.0);
+        if all_zero {
+            for v in slice.iter_mut() {
+                *v = 1.0;
+            }
+        }
+
+        let compensated_distrib = PiecewiseConstant2D::new_with_bounds(&d, domain);
+
+        Self {
+            base,
+            image,
+            image_color_space,
+            scene_center: Point3f::default(),
+            scene_radius: 0.,
+            scale,
+            distrib,
+            compensated_distrib,
+        }
+    }
+
+    fn image_le(&self, uv: Point2f, lambda: &SampledWavelengths) -> SampledSpectrum {
+        let mut rgb = RGB::default();
+        for c in 0..3 {
+            rgb[c] = self.image.lookup_nearest_channel_with_wrap(
+                uv,
+                c,
+                WrapMode::OctahedralSphere.into(),
+            );
+        }
+        let spec =
+            RGBIlluminantSpectrum::new(self.image_color_space.as_ref(), RGB::clamp_zero(rgb));
+        self.scale * spec.sample(lambda)
+    }
+}
+
+impl LightTrait for InfiniteImageLight {
+    fn base(&self) -> &LightBase {
+        &self.base
+    }
+    fn sample_li(
+        &self,
+        ctx: &LightSampleContext,
+        u: Point2f,
+        lambda: &SampledWavelengths,
+        allow_incomplete_pdf: bool,
+    ) -> Option<LightLiSample> {
+        let (uv, map_pdf, _) = if allow_incomplete_pdf {
+            self.compensated_distrib.sample(u)
+        } else {
+            self.distrib.sample(u)
+        };
+
+        if map_pdf == 0. {
+            return None;
+        }
+        // Convert infinite light sample point to direction
+        let w_light = equal_area_square_to_sphere(uv);
+        let wi = self.base.render_from_light.apply_to_vector(w_light);
+        let pdf = map_pdf / (4. * PI);
+
+        // Return radiance value for infinite light direction
+        let mut simple_intr = SimpleInteraction::new_interface(
+            ctx.p() + wi * (2. * self.scene_radius),
+            Some(MediumInterface::default()),
+        );
+
+        simple_intr.common.medium_interface = Some(self.base.medium_interface.clone());
+        let intr = Interaction::Simple(simple_intr);
+        Some(LightLiSample::new(
+            self.image_le(uv, &lambda),
+            wi,
+            pdf,
+            intr,
+        ))
+    }
+
+    fn phi(&self, lambda: SampledWavelengths) -> SampledSpectrum {
+        let mut sum_l = SampledSpectrum::new(0.);
+        let width = self.image.resolution.x();
+        let height = self.image.resolution.y();
+        for v in 0..height {
+            for u in 0..width {
+                let mut rgb = RGB::default();
+                for c in 0..3 {
+                    rgb[c] = self.image.get_channel_with_wrap(
+                        Point2i::new(u, v),
+                        c,
+                        WrapMode::OctahedralSphere.into(),
+                    );
+                }
+                sum_l += RGBIlluminantSpectrum::new(
+                    self.image_color_space.as_ref(),
+                    RGB::clamp_zero(rgb),
+                )
+                .sample(&lambda);
+            }
+        }
+        4. * PI * PI * square(self.scene_radius) * self.scale * sum_l / (width * height) as Float
+    }
+
+    fn pdf_li(&self, _ctx: &LightSampleContext, wi: Vector3f, allow_incomplete_pdf: bool) -> Float {
+        let w_light = self.base.render_from_light.apply_inverse_vector(wi);
+        let uv = equal_area_sphere_to_square(w_light);
+        let pdf = if allow_incomplete_pdf {
+            self.compensated_distrib.pdf(uv)
+        } else {
+            self.distrib.pdf(uv)
+        };
+        pdf / (4. * PI)
+    }
+
+    fn l(
+        &self,
+        _p: Point3f,
+        _n: Normal3f,
+        _uv: Point2f,
+        _w: Vector3f,
+        _lambda: &SampledWavelengths,
+    ) -> SampledSpectrum {
+        todo!()
+    }
+
+    fn le(&self, ray: &Ray, lambda: &SampledWavelengths) -> SampledSpectrum {
+        let w_light = self
+            .base
+            .render_from_light
+            .apply_inverse_vector(ray.d)
+            .normalize();
+        let uv = equal_area_sphere_to_square(w_light);
+        self.image_le(uv, lambda)
+    }
+
+    fn preprocess(&mut self, scene_bounds: &Bounds3f) {
+        let (scene_center, scene_radius) = scene_bounds.bounding_sphere();
+        self.scene_center = scene_center;
+        self.scene_radius = scene_radius;
+    }
+
+    fn bounds(&self) -> Option<LightBounds> {
+        None
+    }
+}
+
+#[derive(Debug, Clone)]
+pub struct InfinitePortalLight {
+    pub base: LightBase,
+    pub image: Image,
+    pub image_color_space: Arc<RGBColorSpace>,
+    pub scale: Float,
+    pub filename: String,
+    pub portal: [Point3f; 4],
+    pub portal_frame: Frame,
+    pub distribution: WindowedPiecewiseConstant2D,
+    pub scene_center: Point3f,
+    pub scene_radius: Float,
+}
+
+impl InfinitePortalLight {
+    fn base(&self) -> &LightBase {
+        &self.base
+    }
+    pub fn new(
+        render_from_light: Transform<Float>,
+        equal_area_image: &Image,
+        image_color_space: Arc<RGBColorSpace>,
+        scale: Float,
+        filename: String,
+        points: Vec<Point3f>,
+    ) -> Self {
+        let base = LightBase::new(
+            LightType::Infinite,
+            &render_from_light,
+            &MediumInterface::default(),
+        );
+
+        let desc = equal_area_image
+            .get_channel_desc(&["R", "G", "B"])
+            .unwrap_or_else(|_| {
+                panic!(
+                    "{}: image used for PortalImageInfiniteLight doesn't have R, G, B channels.",
+                    filename
+                )
+            });
+
+        assert_eq!(3, desc.offset.len());
+        let src_res = equal_area_image.resolution;
+        if src_res.x() != src_res.y() {
+            panic!(
+                "{}: image resolution ({}, {}) is non-square. It's unlikely this is an equal area environment map.",
+                filename,
+                src_res.x(),
+                src_res.y()
+            );
+        }
+
+        if points.len() != 4 {
+            panic!(
+                "Expected 4 vertices for infinite light portal but given {}",
+                points.len()
+            );
+        }
+
+        let portal: [Point3f; 4] = [points[0], points[1], points[2], points[3]];
+
+        let p01 = (portal[1] - portal[0]).normalize();
+        let p12 = (portal[2] - portal[1]).normalize();
+        let p32 = (portal[2] - portal[3]).normalize();
+        let p03 = (portal[3] - portal[0]).normalize();
+
+        if (p01.dot(p32) - 1.0).abs() > 0.001 || (p12.dot(p03) - 1.0).abs() > 0.001 {
+            panic!("Infinite light portal isn't a planar quadrilateral (opposite edges)");
+        }
+
+        if p01.dot(p12).abs() > 0.001
+            || p12.dot(p32).abs() > 0.001
+            || p32.dot(p03).abs() > 0.001
+            || p03.dot(p01).abs() > 0.001
+        {
+            panic!("Infinite light portal isn't a planar quadrilateral (perpendicular edges)");
+        }
+
+        let portal_frame = Frame::from_xy(p03, p01);
+
+        let width = src_res.x();
+        let height = src_res.y();
+
+        let mut new_pixels = vec![0.0 as Float; (width * height * 3) as usize];
+
+        new_pixels
+            .par_chunks_mut((width * 3) as usize)
+            .enumerate()
+            .for_each(|(y, row_pixels)| {
+                let y = y as i32;
+
+                for x in 0..width {
+                    let uv = Point2f::new(
+                        (x as Float + 0.5) / width as Float,
+                        (y as Float + 0.5) / height as Float,
+                    );
+
+                    let (w_world, _) = Self::render_from_image(portal_frame, uv);
+                    let w_local = render_from_light.apply_inverse_vector(w_world).normalize();
+                    let uv_equi = equal_area_sphere_to_square(w_local);
+
+                    let pixel_idx = (x * 3) as usize;
+
+                    for c in 0..3 {
+                        let val = equal_area_image.bilerp_channel_with_wrap(
+                            uv_equi,
+                            c,
+                            WrapMode::OctahedralSphere.into(),
+                        );
+                        row_pixels[pixel_idx + c] = val;
+                    }
+                }
+            });
+
+        let image = Image::new(
+            PixelFormat::F32,
+            src_res,
+            &["R", "G", "B"],
+            equal_area_image.encoding,
+        );
+
+        let duv_dw_closure = |p: Point2f| -> Float {
+            let (_, jacobian) = Self::render_from_image(portal_frame, p);
+            jacobian
+        };
+
+        let d = image.get_sampling_distribution(
+            duv_dw_closure,
+            Bounds2f::from_points(Point2f::new(0., 0.), Point2f::new(1., 1.)),
+        );
+
+        let distribution = WindowedPiecewiseConstant2D::new(d);
+
+        Self {
+            base,
+            image,
+            image_color_space,
+            scale,
+            scene_center: Point3f::default(),
+            scene_radius: 0.,
+            filename,
+            portal,
+            portal_frame,
+            distribution,
+        }
+    }
+
+    pub fn image_lookup(&self, uv: Point2f, lambda: &SampledWavelengths) -> SampledSpectrum {
+        let mut rgb = RGB::default();
+        for c in 0..3 {
+            rgb[c] = self.image.lookup_nearest_channel(uv, c)
+        }
+        let spec =
+            RGBIlluminantSpectrum::new(self.image_color_space.as_ref(), RGB::clamp_zero(rgb));
+        self.scale * spec.sample(lambda)
+    }
+
+    pub fn image_from_render(&self, w_render: Vector3f) -> Option<(Point2f, Float)> {
+        let w = self.portal_frame.to_local(w_render);
+
+        if w.z() <= 0.0 {
+            return None;
+        }
+
+        let alpha = w.x().atan2(w.z());
+        let beta = w.y().atan2(w.z());
+
+        let duv_dw = square(PI) * (1. - square(w.x())) * (1. - square(w.y())) / w.z();
+
+        Some((
+            Point2f::new(
+                clamp_t((alpha + PI / 2.0) / PI, 0.0, 1.0),
+                clamp_t((beta + PI / 2.0) / PI, 0.0, 1.0),
+            ),
+            duv_dw,
+        ))
+    }
+
+    pub fn image_bounds(&self, p: Point3f) -> Option<Bounds2f> {
+        let (p0, _) = self.image_from_render((self.portal[0] - p).normalize())?;
+        let (p1, _) = self.image_from_render((self.portal[2] - p).normalize())?;
+
+        Some(Bounds2f::from_points(p0, p1))
+    }
+
+    pub fn area(&self) -> Float {
+        (self.portal[1] - self.portal[0]).norm() * (self.portal[3] - self.portal[0]).norm()
+    }
+
+    pub fn render_from_image(portal_frame: Frame, uv: Point2f) -> (Vector3f, Float) {
+        let alpha = -PI / 2.0 + uv.x() * PI;
+        let beta = -PI / 2.0 + uv.y() * PI;
+
+        let x = alpha.tan();
+        let y = beta.tan();
+
+        let w = Vector3f::new(x, y, 1.0).normalize();
+
+        let duv_dw = square(PI) * (1.0 - square(w.x())) * (1.0 - square(w.y())) / w.z();
+
+        (portal_frame.from_local(w), duv_dw)
+    }
+}
+
+impl LightTrait for InfinitePortalLight {
+    fn base(&self) -> &LightBase {
+        &self.base
+    }
+    fn sample_li(
+        &self,
+        ctx: &LightSampleContext,
+        u: Point2f,
+        lambda: &SampledWavelengths,
+        _allow_incomplete_pdf: bool,
+    ) -> Option<LightLiSample> {
+        let b = self.image_bounds(ctx.p())?;
+        let (uv, map_pdf) = self.distribution.sample(u, b)?;
+        let (wi, duv_dw) = Self::render_from_image(self.portal_frame, uv);
+        if duv_dw == 0. {
+            return None;
+        }
+        let pdf = map_pdf / duv_dw;
+        let l = self.image_lookup(uv, &lambda);
+        let pl = ctx.p() + 2. * self.scene_radius * wi;
+        let sintr = SimpleInteraction::new_interface(pl, Some(self.base.medium_interface.clone()));
+        let intr = Interaction::Simple(sintr);
+        Some(LightLiSample::new(l, wi, pdf, intr))
+    }
+
+    fn phi(&self, _lambda: SampledWavelengths) -> SampledSpectrum {
+        todo!()
+    }
+
+    fn pdf_li(&self, ctx: &LightSampleContext, wi: Vector3f, _allow_incomplete_pdf: bool) -> Float {
+        let Some((uv, duv_dw)) = self.image_from_render(wi) else {
+            return 0.;
+        };
+        let Some(b) = self.image_bounds(ctx.p()) else {
+            return 0.;
+        };
+        let pdf = self.distribution.pdf(uv, b);
+        pdf / duv_dw
+    }
+
+    fn l(
+        &self,
+        _p: Point3f,
+        _n: Normal3f,
+        _uv: Point2f,
+        _w: Vector3f,
+        _lambda: &SampledWavelengths,
+    ) -> SampledSpectrum {
+        todo!()
+    }
+
+    fn le(&self, ray: &Ray, lambda: &SampledWavelengths) -> SampledSpectrum {
+        let uv = self.image_from_render(ray.d.normalize());
+        let b = self.image_bounds(ray.o);
+        match (uv, b) {
+            (Some((p, duv_dw)), Some(bounds)) if bounds.contains(p) => self.image_lookup(p, lambda),
+            _ => SampledSpectrum::new(0.0),
+        }
+    }
+
+    fn preprocess(&mut self, _scene_bounds: &Bounds3f) {
+        todo!()
+    }
+
+    fn bounds(&self) -> Option<LightBounds> {
+        None
+    }
+}

src/lights/mod.rs

@@ -1,27 +1,918 @@
 pub mod diffuse;
+pub mod infinite;
+pub mod sampler;
 
-#[derive(Debug, Clone)]
+use crate::core::interaction::{
-pub struct DiffuseAreaLight;
+    Interaction, InteractionTrait, MediumInteraction, SimpleInteraction, SurfaceInteraction,
-#[derive(Debug, Clone)]
+};
-pub struct DistantLight;
+use crate::core::medium::MediumInterface;
-#[derive(Debug, Clone)]
+use crate::core::pbrt::{Float, InternCache, PI};
-pub struct GonioPhotometricLight;
+use crate::geometry::{
-#[derive(Debug, Clone)]
+    Bounds2f, Bounds3f, DirectionCone, Normal3f, Point2f, Point2i, Point3f, Point3fi, Ray,
-pub struct InfiniteAreaLight;
+    Vector3f, VectorLike, cos_theta,
-#[derive(Debug, Clone)]
+};
-pub struct PointLight;
+use crate::image::Image;
-#[derive(Debug, Clone)]
+use crate::spectra::{
-pub struct ProjectionLight;
+    DenselySampledSpectrum, LAMBDA_MAX, LAMBDA_MIN, RGBIlluminantSpectrum, SampledSpectrum,
-#[derive(Debug, Clone)]
+    SampledWavelengths, Spectrum, SpectrumTrait,
-pub struct SpotLight;
+};
+use crate::utils::color::RGB;
+use crate::utils::colorspace::RGBColorSpace;
+use crate::utils::math::{equal_area_sphere_to_square, radians, safe_sqrt, smooth_step, square};
+use crate::utils::sampling::PiecewiseConstant2D;
+use crate::utils::transform::Transform;
+use bitflags::bitflags;
 
-#[derive(Debug, Clone)]
+use enum_dispatch::enum_dispatch;
-pub enum Light {
+use std::sync::{Arc, OnceLock};
-    DiffuseArea(Box<DiffuseAreaLight>),
+
-    Distant(Box<DistantLight>),
+use diffuse::DiffuseAreaLight;
-    GonioPhotometric(Box<GonioPhotometricLight>),
+use infinite::{InfiniteImageLight, InfinitePortalLight, InfiniteUniformLight};
-    InfiniteArea(Box<InfiniteAreaLight>),
+
-    Point(Box<PointLight>),
+static SPECTRUM_CACHE: OnceLock<InternCache<DenselySampledSpectrum>> = OnceLock::new();
-    Projection(Box<ProjectionLight>),
+
-    Spot(Box<SpotLight>),
+fn get_spectrum_cache() -> &'static InternCache<DenselySampledSpectrum> {
+    SPECTRUM_CACHE.get_or_init(InternCache::new)
+}
+
+bitflags! {
+    #[derive(Clone, Copy, Debug, PartialEq, Eq, Hash)]
+    pub struct LightType: u32 {
+        const DeltaPosition  = 1;
+        const DeltaDirection = 2;
+        const Area           = 4;
+        const Infinite       = 8;
+    }
+}
+
+impl LightType {
+    pub fn is_infinite(&self) -> bool {
+        self.contains(LightType::Infinite)
+    }
+
+    pub fn is_delta_light(&self) -> bool {
+        self.contains(LightType::DeltaPosition) || self.contains(LightType::DeltaDirection)
+    }
+}
+
+#[derive(Debug, Clone)]
+pub struct LightLeSample {
+    l: SampledSpectrum,
+    ray: Ray,
+    intr: Option<Interaction>,
+    pdf_pos: Float,
+    pdf_dir: Float,
+}
+
+#[derive(Debug, Default, Clone)]
+pub struct LightSampleContext {
+    pub pi: Point3fi,
+    pub n: Normal3f,
+    pub ns: Normal3f,
+}
+
+impl LightSampleContext {
+    pub fn new(pi: Point3fi, n: Normal3f, ns: Normal3f) -> Self {
+        Self { pi, n, ns }
+    }
+
+    pub fn p(&self) -> Point3f {
+        self.pi.into()
+    }
+}
+
+impl From<&SurfaceInteraction> for LightSampleContext {
+    fn from(si: &SurfaceInteraction) -> Self {
+        Self {
+            pi: si.common.pi,
+            n: si.common.n,
+            ns: si.shading.n,
+        }
+    }
+}
+
+impl From<&MediumInteraction> for LightSampleContext {
+    fn from(mi: &MediumInteraction) -> Self {
+        Self {
+            pi: mi.common.pi,
+            n: Normal3f::default(),
+            ns: Normal3f::default(),
+        }
+    }
+}
+
+impl From<&Interaction> for LightSampleContext {
+    fn from(intr: &Interaction) -> Self {
+        match intr {
+            Interaction::Surface(si) => Self {
+                pi: si.common.pi,
+                n: si.common.n,
+                ns: si.shading.n,
+            },
+
+            Interaction::Medium(mi) => Self {
+                pi: mi.common.pi,
+                n: mi.common.n,
+                ns: mi.common.n,
+            },
+
+            Interaction::Simple(sim) => Self {
+                pi: sim.common.pi,
+                n: sim.common.n,
+                ns: sim.common.n,
+            },
+        }
+    }
+}
+
+#[derive(Debug, Clone)]
+pub struct LightLiSample {
+    pub l: SampledSpectrum,
+    pub wi: Vector3f,
+    pub pdf: Float,
+    pub p_light: Arc<Interaction>,
+}
+
+impl LightLiSample {
+    pub fn new(l: SampledSpectrum, wi: Vector3f, pdf: Float, p_light: Interaction) -> Self {
+        Self {
+            l,
+            wi,
+            pdf,
+            p_light: Arc::new(p_light),
+        }
+    }
+}
+
+#[derive(Debug, Clone)]
+pub struct LightBase {
+    pub light_type: LightType,
+    pub render_from_light: Transform<Float>,
+    pub medium_interface: MediumInterface,
+}
+
+impl LightBase {
+    pub fn new(
+        light_type: LightType,
+        render_from_light: &Transform<Float>,
+        medium_interface: &MediumInterface,
+    ) -> Self {
+        Self {
+            light_type,
+            render_from_light: *render_from_light,
+            medium_interface: medium_interface.clone(),
+        }
+    }
+
+    fn l(
+        &self,
+        _p: Point3f,
+        _n: Normal3f,
+        _uv: Point2f,
+        _w: Vector3f,
+        _lambda: &SampledWavelengths,
+    ) -> SampledSpectrum {
+        SampledSpectrum::default()
+    }
+
+    pub fn light_type(&self) -> LightType {
+        self.light_type
+    }
+
+    pub fn le(&self, _ray: &Ray, _lambda: &SampledWavelengths) -> SampledSpectrum {
+        SampledSpectrum::default()
+    }
+
+    pub fn lookup_spectrum(s: &Spectrum) -> Arc<DenselySampledSpectrum> {
+        let cache = SPECTRUM_CACHE.get_or_init(InternCache::new);
+        let dense_spectrum = DenselySampledSpectrum::from_spectrum(s, LAMBDA_MIN, LAMBDA_MAX);
+        cache.lookup(dense_spectrum)
+    }
+}
+
+#[derive(Debug, Clone)]
+pub struct LightBounds {
+    bounds: Bounds3f,
+    phi: Float,
+    w: Vector3f,
+    cos_theta_o: Float,
+    cos_theta_e: Float,
+    two_sided: bool,
+}
+
+impl LightBounds {
+    pub fn new(
+        bounds: &Bounds3f,
+        w: Vector3f,
+        phi: Float,
+        cos_theta_o: Float,
+        cos_theta_e: Float,
+        two_sided: bool,
+    ) -> Self {
+        Self {
+            bounds: *bounds,
+            phi,
+            w,
+            cos_theta_o,
+            cos_theta_e,
+            two_sided,
+        }
+    }
+
+    pub fn centroid(&self) -> Point3f {
+        self.bounds.p_min + Vector3f::from(self.bounds.p_max) / 2.
+    }
+
+    pub fn importance(&self, p: Point3f, n: Normal3f) -> Float {
+        // Compute clamped squared distance to reference point
+        let pc = self.centroid();
+        let d2_raw = p.distance_squared(pc);
+        let d2 = d2_raw.max(self.bounds.diagonal().norm()) / 2.;
+        let cos_sub_clamped = |sin_theta_a: Float,
+                               cos_theta_a: Float,
+                               sin_theta_b: Float,
+                               cos_theta_b: Float|
+         -> Float {
+            if cos_theta_a > cos_theta_b {
+                return 1.;
+            }
+            cos_theta_a * cos_theta_b + sin_theta_a * sin_theta_b
+        };
+
+        let sin_sub_clamped = |sin_theta_a: Float,
+                               cos_theta_a: Float,
+                               sin_theta_b: Float,
+                               cos_theta_b: Float|
+         -> Float {
+            if cos_theta_a > cos_theta_b {
+                return 1.;
+            }
+            sin_theta_a * cos_theta_b - cos_theta_a * sin_theta_b
+        };
+
+        let wi = (p - pc).normalize();
+        let mut cos_theta_w = self.w.dot(wi);
+        if self.two_sided {
+            cos_theta_w = cos_theta_w.abs();
+        }
+        let sin_theta_w = safe_sqrt(1. - square(cos_theta_w));
+        let cos_theta_b = DirectionCone::bound_subtended_directions(&self.bounds, p).cos_theta;
+        let sin_theta_b = safe_sqrt(1. - square(cos_theta_b));
+        let sin_theta_o = safe_sqrt(1. - square(self.cos_theta_o));
+        let cos_theta_x = cos_sub_clamped(sin_theta_w, cos_theta_w, sin_theta_o, self.cos_theta_o);
+        let sin_theta_x = sin_sub_clamped(sin_theta_w, cos_theta_w, sin_theta_o, self.cos_theta_o);
+        let cos_theta_p = cos_sub_clamped(sin_theta_x, cos_theta_x, sin_theta_b, cos_theta_b);
+
+        if cos_theta_p <= self.cos_theta_e {
+            return 0.;
+        }
+
+        let mut importance = self.phi * cos_theta_p / d2;
+        if n != Normal3f::new(0., 0., 0.) {
+            let cos_theta_i = wi.abs_dot(n.into());
+            let sin_theta_i = safe_sqrt(1. - square(cos_theta_i));
+            let cos_thetap_i = cos_sub_clamped(sin_theta_i, cos_theta_i, sin_theta_b, cos_theta_b);
+            importance *= cos_thetap_i;
+        }
+        importance
+    }
+
+    pub fn union(a: &Self, b: &Self) -> Self {
+        if a.phi == 0. {
+            return a.clone();
+        }
+        if b.phi == 0. {
+            return b.clone();
+        }
+
+        let a_cone = DirectionCone::new(a.w, a.cos_theta_o);
+        let b_cone = DirectionCone::new(b.w, b.cos_theta_o);
+        let cone = DirectionCone::union(&a_cone, &b_cone);
+        let cos_theta_o = cone.cos_theta;
+        let cos_theta_e = a.cos_theta_e.min(b.cos_theta_e);
+        LightBounds::new(
+            &a.bounds.union(b.bounds),
+            cone.w,
+            a.phi + b.phi,
+            cos_theta_o,
+            cos_theta_e,
+            a.two_sided || b.two_sided,
+        )
+    }
+}
+
+#[enum_dispatch]
+pub trait LightTrait: Send + Sync + std::fmt::Debug {
+    fn base(&self) -> &LightBase;
+    fn phi(&self, lambda: SampledWavelengths) -> SampledSpectrum;
+    fn sample_li(
+        &self,
+        ctx: &LightSampleContext,
+        u: Point2f,
+        lambda: &SampledWavelengths,
+        allow_incomplete_pdf: bool,
+    ) -> Option<LightLiSample>;
+    fn pdf_li(&self, ctx: &LightSampleContext, wi: Vector3f, allow_incomplete_pdf: bool) -> Float;
+    fn l(
+        &self,
+        p: Point3f,
+        n: Normal3f,
+        uv: Point2f,
+        w: Vector3f,
+        lambda: &SampledWavelengths,
+    ) -> SampledSpectrum;
+    fn le(&self, ray: &Ray, lambda: &SampledWavelengths) -> SampledSpectrum;
+    fn preprocess(&mut self, scene_bounds: &Bounds3f);
+    fn bounds(&self) -> Option<LightBounds>;
+    fn light_type(&self) -> LightType {
+        self.base().light_type()
+    }
+}
+
+#[derive(Debug, Clone)]
+#[enum_dispatch(LightTrait)]
+#[allow(clippy::large_enum_variant)]
+pub enum Light {
+    DiffuseArea(DiffuseAreaLight),
+    Distant(DistantLight),
+    Goniometric(GoniometricLight),
+    InfiniteUniform(InfiniteUniformLight),
+    InfiniteImage(InfiniteImageLight),
+    InfinitePortal(InfinitePortalLight),
+    Point(PointLight),
+    Projection(ProjectionLight),
+    Spot(SpotLight),
+}
+
+#[derive(Debug, Clone)]
+pub struct DistantLight {
+    pub base: LightBase,
+    lemit: Arc<DenselySampledSpectrum>,
+    scale: Float,
+    scene_center: Point3f,
+    scene_radius: Float,
+}
+
+impl DistantLight {
+    pub fn new(render_from_light: &Transform<Float>, lemit: Spectrum, scale: Float) -> Self {
+        let base = LightBase::new(
+            LightType::DeltaDirection,
+            render_from_light,
+            &MediumInterface::default(),
+        );
+
+        let l_interned = LightBase::lookup_spectrum(&lemit);
+        Self {
+            base,
+            lemit: l_interned,
+            scale,
+            scene_center: Point3f::default(),
+            scene_radius: 0.,
+        }
+    }
+}
+
+impl LightTrait for DistantLight {
+    fn base(&self) -> &LightBase {
+        &self.base
+    }
+
+    fn phi(&self, lambda: SampledWavelengths) -> SampledSpectrum {
+        self.scale * self.lemit.sample(&lambda) * PI * self.scene_radius.sqrt()
+    }
+
+    fn sample_li(
+        &self,
+        ctx: &LightSampleContext,
+        _u: Point2f,
+        lambda: &SampledWavelengths,
+        _allow_incomplete_pdf: bool,
+    ) -> Option<LightLiSample> {
+        let wi = self
+            .base
+            .render_from_light
+            .apply_to_vector(Vector3f::new(0., 0., 1.))
+            .normalize();
+        let p_outside = ctx.p() + wi * 2. * self.scene_radius;
+
+        let li = self.scale * self.lemit.sample(&lambda);
+        let intr = SimpleInteraction::new(
+            Point3fi::new_from_point(p_outside),
+            0.0,
+            Some(self.base.medium_interface.clone()),
+        );
+
+        Some(LightLiSample::new(li, wi, 1., Interaction::Simple(intr)))
+    }
+
+    fn pdf_li(
+        &self,
+        _ctx: &LightSampleContext,
+        _wi: Vector3f,
+        _allow_incomplete_pdf: bool,
+    ) -> Float {
+        0.
+    }
+
+    fn l(
+        &self,
+        _p: Point3f,
+        _n: Normal3f,
+        _uv: Point2f,
+        _w: Vector3f,
+        _lambda: &SampledWavelengths,
+    ) -> SampledSpectrum {
+        todo!()
+    }
+    fn le(&self, _ray: &Ray, _lambda: &SampledWavelengths) -> SampledSpectrum {
+        todo!()
+    }
+
+    fn preprocess(&mut self, scene_bounds: &Bounds3f) {
+        let (center, radius) = scene_bounds.bounding_sphere();
+        self.scene_center = center;
+        self.scene_radius = radius;
+    }
+
+    fn bounds(&self) -> Option<LightBounds> {
+        None
+    }
+}
+
+#[derive(Debug, Clone)]
+pub struct GoniometricLight {
+    pub base: LightBase,
+    iemit: Arc<DenselySampledSpectrum>,
+    scale: Float,
+    image: Image,
+    distrib: PiecewiseConstant2D,
+}
+
+impl GoniometricLight {
+    pub fn new(
+        render_from_light: &Transform<Float>,
+        medium_interface: &MediumInterface,
+        iemit: Spectrum,
+        scale: Float,
+        image: Image,
+    ) -> Self {
+        let base = LightBase::new(
+            LightType::DeltaPosition,
+            render_from_light,
+            medium_interface,
+        );
+
+        let i_interned = LightBase::lookup_spectrum(&iemit);
+        let d = image.get_sampling_distribution_uniform();
+        let distrib = PiecewiseConstant2D::new_with_data(&d);
+        Self {
+            base,
+            iemit: i_interned,
+            scale,
+            image,
+            distrib,
+        }
+    }
+
+    pub fn i(&self, w: Vector3f, lambda: &SampledWavelengths) -> SampledSpectrum {
+        let uv = equal_area_sphere_to_square(w);
+        self.scale * self.iemit.sample(&lambda) * self.image.lookup_nearest_channel(uv, 0)
+    }
+}
+
+impl LightTrait for GoniometricLight {
+    fn base(&self) -> &LightBase {
+        &self.base
+    }
+
+    fn phi(&self, lambda: SampledWavelengths) -> SampledSpectrum {
+        let mut sum_y = 0.;
+        for y in 0..self.image.resolution.y() {
+            for x in 0..self.image.resolution.x() {
+                sum_y += self.image.get_channel(Point2i::new(x, y), 0);
+            }
+        }
+        self.scale * self.iemit.sample(&lambda) * 4. * PI * sum_y
+            / (self.image.resolution.x() * self.image.resolution.y()) as Float
+    }
+
+    fn sample_li(
+        &self,
+        _ctx: &LightSampleContext,
+        _u: Point2f,
+        _lambda: &SampledWavelengths,
+        _allow_incomplete_pdf: bool,
+    ) -> Option<LightLiSample> {
+        todo!()
+    }
+
+    fn pdf_li(
+        &self,
+        _ctx: &LightSampleContext,
+        _wi: Vector3f,
+        _allow_incomplete_pdf: bool,
+    ) -> Float {
+        0.
+    }
+
+    fn l(
+        &self,
+        _p: Point3f,
+        _n: Normal3f,
+        _uv: Point2f,
+        _w: Vector3f,
+        _lambda: &SampledWavelengths,
+    ) -> SampledSpectrum {
+        todo!()
+    }
+    fn le(&self, _ray: &Ray, _lambda: &SampledWavelengths) -> SampledSpectrum {
+        todo!()
+    }
+    fn preprocess(&mut self, _scene_bounds: &Bounds3f) {
+        todo!()
+    }
+    fn bounds(&self) -> Option<LightBounds> {
+        todo!()
+    }
+}
+
+#[derive(Debug, Clone)]
+pub struct PointLight {
+    base: LightBase,
+    i: Arc<DenselySampledSpectrum>,
+    scale: Float,
+}
+
+impl PointLight {
+    pub fn new(
+        render_from_light: Transform<Float>,
+        medium_interface: MediumInterface,
+        i: &Spectrum,
+        scale: Float,
+    ) -> Self {
+        let base = LightBase::new(
+            LightType::DeltaPosition,
+            &render_from_light,
+            &medium_interface,
+        );
+
+        let i_interned = LightBase::lookup_spectrum(i);
+
+        Self {
+            base,
+            i: i_interned,
+            scale,
+        }
+    }
+}
+
+impl LightTrait for PointLight {
+    fn base(&self) -> &LightBase {
+        &self.base
+    }
+
+    fn sample_li(
+        &self,
+        ctx: &LightSampleContext,
+        _u: Point2f,
+        lambda: &SampledWavelengths,
+        _allow_incomplete_pdf: bool,
+    ) -> Option<LightLiSample> {
+        let pi = self
+            .base
+            .render_from_light
+            .apply_to_interval(&Point3fi::default());
+        let p: Point3f = pi.into();
+        let wi = (p - ctx.p()).normalize();
+        let li = self.scale * self.i.sample(&lambda) / p.distance_squared(ctx.p());
+        let intr = SimpleInteraction::new(pi, 0.0, Some(self.base.medium_interface.clone()));
+        Some(LightLiSample::new(li, wi, 1., Interaction::Simple(intr)))
+    }
+
+    fn phi(&self, lambda: SampledWavelengths) -> SampledSpectrum {
+        4. * PI * self.scale * self.i.sample(&lambda)
+    }
+
+    fn pdf_li(
+        &self,
+        _ctx: &LightSampleContext,
+        _wi: Vector3f,
+        _allow_incomplete_pdf: bool,
+    ) -> Float {
+        0.
+    }
+
+    fn l(
+        &self,
+        _p: Point3f,
+        _n: Normal3f,
+        _uv: Point2f,
+        _w: Vector3f,
+        _lambda: &SampledWavelengths,
+    ) -> SampledSpectrum {
+        todo!()
+    }
+    fn le(&self, _ray: &Ray, _lambda: &SampledWavelengths) -> SampledSpectrum {
+        todo!()
+    }
+    fn preprocess(&mut self, _scene_bounds: &Bounds3f) {
+        todo!()
+    }
+    fn bounds(&self) -> Option<LightBounds> {
+        let p = self
+            .base
+            .render_from_light
+            .apply_to_point(Point3f::new(0., 0., 0.));
+        let phi = 4. * PI * self.scale * self.i.max_value();
+        Some(LightBounds::new(
+            &Bounds3f::from_points(p, p),
+            Vector3f::new(0., 0., 1.),
+            phi,
+            PI.cos(),
+            (PI / 2.).cos(),
+            false,
+        ))
+    }
+}
+
+#[derive(Debug, Clone)]
+pub struct ProjectionLight {
+    base: LightBase,
+    image: Image,
+    image_color_space: Arc<RGBColorSpace>,
+    scale: Float,
+    screen_bounds: Bounds2f,
+    hither: Float,
+    screen_from_light: Transform<Float>,
+    light_from_screen: Transform<Float>,
+    a: Float,
+    distrib: PiecewiseConstant2D,
+}
+
+impl ProjectionLight {
+    pub fn new(
+        render_from_light: Transform<Float>,
+        medium_interface: MediumInterface,
+        image: Image,
+        image_color_space: Arc<RGBColorSpace>,
+        scale: Float,
+        fov: Float,
+    ) -> Self {
+        let base = LightBase::new(
+            LightType::DeltaPosition,
+            &render_from_light,
+            &medium_interface,
+        );
+        let aspect = image.resolution().x() as Float / image.resolution().y() as Float;
+        let screen_bounds = if aspect > 1. {
+            Bounds2f::from_points(Point2f::new(-aspect, -1.), Point2f::new(aspect, 1.))
+        } else {
+            Bounds2f::from_points(
+                Point2f::new(-1., 1. / aspect),
+                Point2f::new(1., 1. / aspect),
+            )
+        };
+
+        let hither = 1e-3;
+        let screen_from_light = Transform::perspective(fov, hither, 1e30).unwrap();
+        let light_from_screen = screen_from_light.inverse();
+        let opposite = (radians(fov) / 2.).tan();
+        let aspect_ratio = if aspect > 1. { aspect } else { 1. / aspect };
+        let a = 4. * square(opposite) * aspect_ratio;
+        let dwda = |p: Point2f| {
+            let w =
+                Vector3f::from(light_from_screen.apply_to_point(Point3f::new(p.x(), p.y(), 0.)));
+            cos_theta(w.normalize()).powi(3)
+        };
+
+        let d = image.get_sampling_distribution(dwda, screen_bounds);
+        let distrib = PiecewiseConstant2D::new_with_bounds(&d, screen_bounds);
+
+        Self {
+            base,
+            image,
+            image_color_space,
+            screen_bounds,
+            screen_from_light,
+            light_from_screen,
+            scale,
+            hither,
+            a,
+            distrib,
+        }
+    }
+
+    pub fn i(&self, w: Vector3f, lambda: SampledWavelengths) -> SampledSpectrum {
+        if w.z() < self.hither {
+            return SampledSpectrum::new(0.);
+        }
+        let ps = self.screen_from_light.apply_to_point(w.into());
+        if !self.screen_bounds.contains(Point2f::new(ps.x(), ps.y())) {
+            return SampledSpectrum::new(0.);
+        }
+        let uv = Point2f::from(self.screen_bounds.offset(&Point2f::new(ps.x(), ps.y())));
+        let mut rgb = RGB::default();
+        for c in 0..3 {
+            rgb[c] = self.image.lookup_nearest_channel(uv, c);
+        }
+        let s = RGBIlluminantSpectrum::new(self.image_color_space.as_ref(), RGB::clamp_zero(rgb));
+        self.scale * s.sample(&lambda)
+    }
+}
+
+impl LightTrait for ProjectionLight {
+    fn base(&self) -> &LightBase {
+        &self.base
+    }
+    fn phi(&self, lambda: SampledWavelengths) -> SampledSpectrum {
+        let mut sum = SampledSpectrum::new(0.);
+        for y in 0..self.image.resolution.y() {
+            for x in 0..self.image.resolution.x() {
+                let ps = self.screen_bounds.lerp(Point2f::new(
+                    (x as Float + 0.5) / self.image.resolution.x() as Float,
+                    (y as Float + 0.5) / self.image.resolution.y() as Float,
+                ));
+                let w_raw = Vector3f::from(self.light_from_screen.apply_to_point(Point3f::new(
+                    ps.x(),
+                    ps.y(),
+                    0.,
+                )));
+                let w = w_raw.normalize();
+                let dwda = cos_theta(w).powi(3);
+                let mut rgb = RGB::default();
+                for c in 0..3 {
+                    rgb[c] = self.image.get_channel(Point2i::new(x, y), c).into();
+                }
+
+                let s = RGBIlluminantSpectrum::new(
+                    self.image_color_space.as_ref(),
+                    RGB::clamp_zero(rgb),
+                );
+                sum += s.sample(&lambda) * dwda;
+            }
+        }
+        self.scale * self.a * sum / (self.image.resolution.x() * self.image.resolution.y()) as Float
+    }
+
+    fn sample_li(
+        &self,
+        _ctx: &LightSampleContext,
+        _u: Point2f,
+        _lambda: &SampledWavelengths,
+        _allow_incomplete_pdf: bool,
+    ) -> Option<LightLiSample> {
+        todo!()
+    }
+    fn pdf_li(
+        &self,
+        _ctx: &LightSampleContext,
+        _wi: Vector3f,
+        _allow_incomplete_pdf: bool,
+    ) -> Float {
+        todo!()
+    }
+    fn l(
+        &self,
+        _p: Point3f,
+        _n: Normal3f,
+        _uv: Point2f,
+        _w: Vector3f,
+        _lambda: &SampledWavelengths,
+    ) -> SampledSpectrum {
+        todo!()
+    }
+    fn le(&self, _ray: &Ray, _lambda: &SampledWavelengths) -> SampledSpectrum {
+        todo!()
+    }
+    fn preprocess(&mut self, _scene_bounds: &Bounds3f) {
+        todo!()
+    }
+    fn bounds(&self) -> Option<LightBounds> {
+        todo!()
+    }
+}
+
+#[derive(Debug, Clone)]
+pub struct SpotLight {
+    base: LightBase,
+    iemit: Arc<DenselySampledSpectrum>,
+    scale: Float,
+    cos_fallof_start: Float,
+    cos_fallof_end: Float,
+}
+
+impl SpotLight {
+    pub fn new(
+        render_from_light: &Transform<Float>,
+        medium_interface: &MediumInterface,
+        iemit: Spectrum,
+        scale: Float,
+        total_width: Float,
+        fallof_start: Float,
+    ) -> Self {
+        let base = LightBase::new(
+            LightType::DeltaPosition,
+            render_from_light,
+            medium_interface,
+        );
+
+        let i_interned = LightBase::lookup_spectrum(&iemit);
+        let cos_fallof_end = radians(total_width).cos();
+        let cos_fallof_start = radians(fallof_start).cos();
+        assert!(fallof_start < total_width);
+        Self {
+            base,
+            iemit: i_interned,
+            scale,
+            cos_fallof_start,
+            cos_fallof_end,
+        }
+    }
+
+    pub fn i(&self, w: Vector3f, lambda: SampledWavelengths) -> SampledSpectrum {
+        smooth_step(cos_theta(w), self.cos_fallof_end, self.cos_fallof_start)
+            * self.scale
+            * self.iemit.sample(&lambda)
+    }
+}
+
+impl LightTrait for SpotLight {
+    fn base(&self) -> &LightBase {
+        &self.base
+    }
+    fn sample_li(
+        &self,
+        ctx: &LightSampleContext,
+        _u: Point2f,
+        lambda: &SampledWavelengths,
+        _allow_incomplete_pdf: bool,
+    ) -> Option<LightLiSample> {
+        let pi = self
+            .base
+            .render_from_light
+            .apply_to_interval(&Point3fi::default());
+        let p: Point3f = pi.into();
+        let wi = (p - ctx.p()).normalize();
+        let w_light = self.base.render_from_light.apply_inverse_vector(-wi);
+        let li = self.i(w_light, *lambda) / p.distance_squared(ctx.p());
+
+        let intr = SimpleInteraction::new(pi, 0.0, Some(self.base.medium_interface.clone()));
+        Some(LightLiSample::new(li, wi, 1., Interaction::Simple(intr)))
+    }
+
+    fn phi(&self, lambda: SampledWavelengths) -> SampledSpectrum {
+        self.scale
+            * self.iemit.sample(&lambda)
+            * 2.
+            * PI
+            * ((1. - self.cos_fallof_start) + (self.cos_fallof_start - self.cos_fallof_end) / 2.)
+    }
+
+    fn pdf_li(
+        &self,
+        _ctx: &LightSampleContext,
+        _wi: Vector3f,
+        _allow_incomplete_pdf: bool,
+    ) -> Float {
+        0.
+    }
+
+    fn l(
+        &self,
+        _p: Point3f,
+        _n: Normal3f,
+        _uv: Point2f,
+        _w: Vector3f,
+        _lambda: &SampledWavelengths,
+    ) -> SampledSpectrum {
+        todo!()
+    }
+    fn le(&self, _ray: &Ray, _lambda: &SampledWavelengths) -> SampledSpectrum {
+        todo!()
+    }
+    fn preprocess(&mut self, _scene_bounds: &Bounds3f) {
+        todo!()
+    }
+
+    fn bounds(&self) -> Option<LightBounds> {
+        let p = self
+            .base
+            .render_from_light
+            .apply_to_point(Point3f::default());
+        let w = self
+            .base
+            .render_from_light
+            .apply_to_vector(Vector3f::new(0., 0., 1.))
+            .normalize();
+        let phi = self.scale * self.iemit.max_value() * 4. * PI;
+        let cos_theta_e = (self.cos_fallof_end.acos() - self.cos_fallof_start.acos()).cos();
+        Some(LightBounds::new(
+            &Bounds3f::from_points(p, p),
+            w,
+            phi,
+            self.cos_fallof_start,
+            cos_theta_e,
+            false,
+        ))
+    }
 }

src/lights/sampler.rs

@@ -0,0 +1,521 @@
+use super::{
+    Bounds3f, Float, Light, LightBounds, LightSampleContext, LightTrait, Normal3f, PI, Point3f,
+    SampledSpectrum, SampledWavelengths, Vector3f, VectorLike, safe_sqrt, square,
+};
+use crate::geometry::primitives::OctahedralVector;
+use crate::geometry::{DirectionCone, Normal};
+use crate::utils::math::sample_discrete;
+use std::collections::HashMap;
+use std::sync::Arc;
+
+use crate::core::pbrt::{ONE_MINUS_EPSILON, clamp_t, lerp};
+use crate::utils::sampling::AliasTable;
+use enum_dispatch::enum_dispatch;
+
+#[derive(Clone, Copy, Debug, Default)]
+#[repr(C)]
+pub struct CompactLightBounds {
+    pub w: OctahedralVector,
+    pub phi: Float,
+
+    // [0..15]  = qCosTheta_o
+    // [15..30] = qCosTheta_e
+    // [30..31] = twoSided
+    // [31..32] = Unused/Padding
+    packed_info: u32,
+
+    pub qb: [[u16; 3]; 2],
+}
+
+const _: () = assert!(std::mem::size_of::<CompactLightBounds>() == 24);
+
+impl CompactLightBounds {
+    pub fn new(lb: &LightBounds, all_b: &Bounds3f) -> Self {
+        let q_cos_o = Self::quantize_cos(lb.cos_theta_o);
+        let q_cos_e = Self::quantize_cos(lb.cos_theta_e);
+        let two_sided = if lb.two_sided { 1 } else { 0 };
+
+        // | - twoSided (1) - | - qCosTheta_e (15) - | - qCosTheta_o (15) - |
+        let packed_info = (q_cos_o & 0x7FFF) | ((q_cos_e & 0x7FFF) << 15) | (two_sided << 30);
+
+        let mut qb = [[0u16; 3]; 2];
+        for i in 0..3 {
+            qb[0][i] = Self::quantize_bounds(lb.bounds.p_min[i], all_b.p_min[i], all_b.p_max[i])
+                .floor() as u16;
+            qb[1][i] = Self::quantize_bounds(lb.bounds.p_max[i], all_b.p_min[i], all_b.p_max[i])
+                .ceil() as u16;
+        }
+
+        Self {
+            w: OctahedralVector::new(lb.w.normalize()),
+            phi: lb.phi,
+            packed_info,
+            qb,
+        }
+    }
+
+    #[inline(always)]
+    pub fn two_sided(&self) -> bool {
+        (self.packed_info >> 30) & 1 == 1
+    }
+
+    #[inline(always)]
+    fn q_cos_theta_o(&self) -> u32 {
+        self.packed_info & 0x7FFF
+    }
+
+    #[inline(always)]
+    fn q_cos_theta_e(&self) -> u32 {
+        (self.packed_info >> 15) & 0x7FFF
+    }
+
+    #[inline]
+    pub fn cos_theta_o(&self) -> Float {
+        2.0 * (self.q_cos_theta_o() as Float / 32767.0) - 1.0
+    }
+
+    #[inline]
+    pub fn cos_theta_e(&self) -> Float {
+        2.0 * (self.q_cos_theta_e() as Float / 32767.0) - 1.0
+    }
+
+    pub fn importance(&self, p: Point3f, n: Normal3f, all_b: &Bounds3f) -> Float {
+        let bounds = self.bounds(all_b);
+        let cos_o = self.cos_theta_o();
+        let cos_e = self.cos_theta_e();
+
+        let pc = bounds.centroid();
+        let d2 = p.distance_squared(pc).max(bounds.diagonal().norm() * 0.5);
+
+        let cos_sub_clamped = |sin_a: Float, cos_a: Float, sin_b: Float, cos_b: Float| {
+            if cos_a > cos_b {
+                1.0
+            } else {
+                cos_a * cos_b + sin_a * sin_b
+            }
+        };
+        let sin_sub_clamped = |sin_a: Float, cos_a: Float, sin_b: Float, cos_b: Float| {
+            if cos_a > cos_b {
+                0.0
+            } else {
+                sin_a * cos_b - cos_a * sin_b
+            }
+        };
+
+        let wi = (p - pc).normalize();
+        let w_vec = self.w.to_vector();
+        let mut cos_w = w_vec.dot(wi);
+        if self.two_sided() {
+            cos_w = cos_w.abs();
+        }
+        let sin_w = safe_sqrt(1.0 - square(cos_w));
+        let cos_b = DirectionCone::bound_subtended_directions(&bounds, p).cos_theta;
+        let sin_b = safe_sqrt(1. - square(cos_b));
+        let sin_o = safe_sqrt(1. - square(cos_o));
+        let cos_x = cos_sub_clamped(sin_w, cos_w, sin_o, cos_o);
+        let sin_x = sin_sub_clamped(sin_w, cos_w, sin_o, cos_o);
+        let cos_p = cos_sub_clamped(sin_x, cos_x, sin_b, cos_b);
+        if cos_p <= cos_e {
+            return 0.;
+        }
+        let mut importance = self.phi * cos_p / d2;
+        if n != Normal3f::zero() {
+            let cos_i = wi.abs_dot(n.into());
+            let sin_i = safe_sqrt(1. - square(cos_i));
+            let cos_pi = cos_sub_clamped(sin_i, cos_i, sin_b, cos_b);
+            importance *= cos_pi;
+        }
+        importance
+    }
+
+    pub fn bounds(&self, all_b: &Bounds3f) -> Bounds3f {
+        let mut p_min = Point3f::default();
+        let mut p_max = Point3f::default();
+
+        for i in 0..3 {
+            let t_min = self.qb[0][i] as Float / 65535.0;
+            let t_max = self.qb[1][i] as Float / 65535.0;
+            p_min[i] = lerp(t_min, all_b.p_min[i], all_b.p_max[i]);
+            p_max[i] = lerp(t_max, all_b.p_min[i], all_b.p_max[i]);
+        }
+        Bounds3f::from_points(p_min, p_max)
+    }
+
+    fn quantize_cos(c: Float) -> u32 {
+        (32767.0 * ((c + 1.0) * 0.5)).floor() as u32
+    }
+
+    fn quantize_bounds(c: Float, min: Float, max: Float) -> Float {
+        if min == max {
+            return 0.0;
+        }
+        65535.0 * clamp_t((c - min) / (max - min), 0.0, 1.0)
+    }
+}
+
+#[derive(Debug, Clone)]
+pub struct SampledLight {
+    pub light: Arc<Light>,
+    pub p: Float,
+}
+
+impl SampledLight {
+    pub fn new(light: Arc<Light>, p: Float) -> Self {
+        Self { light, p }
+    }
+}
+
+#[enum_dispatch]
+pub trait LightSamplerTrait: Send + Sync + std::fmt::Debug {
+    fn sample_with_context(&self, ctx: &LightSampleContext, u: Float) -> Option<SampledLight>;
+    fn pmf_with_context(&self, ctx: &LightSampleContext, light: &Arc<Light>) -> Float;
+    fn sample(&self, u: Float) -> Option<SampledLight>;
+    fn pmf(&self, light: &Arc<Light>) -> Float;
+}
+
+#[derive(Clone, Debug)]
+#[enum_dispatch(LightSamplerTrait)]
+pub enum LightSampler {
+    Uniform(UniformLightSampler),
+    Power(PowerLightSampler),
+    BVH(BVHLightSampler),
+}
+
+#[derive(Clone, Debug)]
+pub struct UniformLightSampler {
+    lights: Vec<Arc<Light>>,
+}
+
+impl UniformLightSampler {
+    pub fn new(lights: &[Arc<Light>]) -> Self {
+        Self {
+            lights: lights.to_vec(),
+        }
+    }
+}
+
+impl LightSamplerTrait for UniformLightSampler {
+    fn sample_with_context(&self, _ctx: &LightSampleContext, u: Float) -> Option<SampledLight> {
+        self.sample(u)
+    }
+    fn pmf_with_context(&self, _ctx: &LightSampleContext, light: &Arc<Light>) -> Float {
+        self.pmf(light)
+    }
+    fn sample(&self, u: Float) -> Option<SampledLight> {
+        if self.lights.is_empty() {
+            return None;
+        }
+
+        let light_index = (u as usize * self.lights.len()).min(self.lights.len() - 1);
+        Some(SampledLight {
+            light: self.lights[light_index].clone(),
+            p: 1. / self.lights.len() as Float,
+        })
+    }
+    fn pmf(&self, _light: &Arc<Light>) -> Float {
+        if self.lights.is_empty() {
+            return 0.;
+        }
+        1. / self.lights.len() as Float
+    }
+}
+
+#[derive(Clone, Debug)]
+pub struct PowerLightSampler {
+    lights: Vec<Arc<Light>>,
+    light_to_index: HashMap<usize, usize>,
+    alias_table: AliasTable,
+}
+
+impl PowerLightSampler {
+    pub fn new(lights: &[Arc<Light>]) -> Self {
+        if lights.is_empty() {
+            return Self {
+                lights: Vec::new(),
+                light_to_index: HashMap::new(),
+                alias_table: AliasTable::new(&[]),
+            };
+        }
+
+        let mut lights_vec = Vec::with_capacity(lights.len());
+        let mut light_to_index = HashMap::with_capacity(lights.len());
+        let mut light_power = Vec::with_capacity(lights.len());
+
+        let lambda = SampledWavelengths::sample_visible(0.5);
+
+        for (i, light) in lights.iter().enumerate() {
+            lights_vec.push(light.clone());
+
+            let ptr = Arc::as_ptr(light) as usize;
+            light_to_index.insert(ptr, i);
+
+            let phi = SampledSpectrum::safe_div(&light.phi(lambda), &lambda.pdf());
+            light_power.push(phi.average());
+        }
+
+        let alias_table = AliasTable::new(&light_power);
+
+        Self {
+            lights: lights_vec,
+            light_to_index,
+            alias_table,
+        }
+    }
+}
+
+impl LightSamplerTrait for PowerLightSampler {
+    fn sample_with_context(&self, _ctx: &LightSampleContext, u: Float) -> Option<SampledLight> {
+        self.sample(u)
+    }
+
+    fn pmf_with_context(&self, _ctx: &LightSampleContext, light: &Arc<Light>) -> Float {
+        self.pmf(light)
+    }
+
+    fn sample(&self, u: Float) -> Option<SampledLight> {
+        if self.alias_table.size() == 0 {
+            return None;
+        }
+
+        let (light_index, pmf, _) = self.alias_table.sample(u);
+
+        Some(SampledLight {
+            light: self.lights[light_index].clone(),
+            p: pmf,
+        })
+    }
+
+    fn pmf(&self, light: &Arc<Light>) -> Float {
+        if self.alias_table.size() == 0 {
+            return 0.;
+        }
+
+        let ptr = Arc::as_ptr(light) as usize;
+
+        if let Some(&index) = self.light_to_index.get(&ptr) {
+            self.alias_table.pmf(index)
+        } else {
+            0.0
+        }
+    }
+}
+
+#[derive(Clone, Copy, Debug, Default)]
+#[repr(C, align(32))]
+pub struct LightBVHNode {
+    pub light_bounds: CompactLightBounds,
+
+    // Bit 31 (MSB) : isLeaf (1 bit)
+    // Bits 0..31   : childOrLightIndex (31 bits)
+    packed_data: u32,
+}
+
+const _: () = assert!(std::mem::size_of::<LightBVHNode>() == 32);
+
+impl LightBVHNode {
+    /// Mask to isolate the Leaf Flag (Bit 31)
+    const LEAF_MASK: u32 = 0x8000_0000;
+    /// Mask to isolate the Index (Bits 0-30)
+    const INDEX_MASK: u32 = 0x7FFF_FFFF;
+
+    pub fn make_leaf(light_index: u32, cb: CompactLightBounds) -> Self {
+        debug_assert!(
+            (light_index & Self::LEAF_MASK) == 0,
+            "Light index too large"
+        );
+
+        Self {
+            light_bounds: cb,
+            // Set index and flip the MSB to 1
+            packed_data: light_index | Self::LEAF_MASK,
+        }
+    }
+
+    pub fn make_interior(child_index: u32, cb: CompactLightBounds) -> Self {
+        debug_assert!(
+            (child_index & Self::LEAF_MASK) == 0,
+            "Child index too large"
+        );
+
+        Self {
+            light_bounds: cb,
+            // Set index, MSB remains 0
+            packed_data: child_index,
+        }
+    }
+
+    #[inline(always)]
+    pub fn is_leaf(&self) -> bool {
+        (self.packed_data & Self::LEAF_MASK) != 0
+    }
+
+    #[inline(always)]
+    pub fn light_index(&self) -> u32 {
+        debug_assert!(self.is_leaf());
+        self.packed_data & Self::INDEX_MASK
+    }
+
+    #[inline(always)]
+    pub fn child_index(&self) -> u32 {
+        debug_assert!(!self.is_leaf());
+        self.packed_data & Self::INDEX_MASK
+    }
+
+    #[inline(always)]
+    pub fn child_or_light_index(&self) -> u32 {
+        self.packed_data & Self::INDEX_MASK
+    }
+
+    pub fn sample(&self, _ctx: &LightSampleContext, _u: Float) -> Option<SampledLight> {
+        todo!("Implement LightBVHNode::Sample logic")
+    }
+}
+
+#[derive(Clone, Debug)]
+pub struct BVHLightSampler {
+    lights: Vec<Arc<Light>>,
+    infinite_lights: Vec<Arc<Light>>,
+    all_light_bounds: Bounds3f,
+    nodes: Vec<LightBVHNode>,
+    light_to_bit_trail: HashMap<usize, usize>,
+}
+
+impl BVHLightSampler {
+    fn evaluate_cost(&self, b: &LightBounds, bounds: &Bounds3f, dim: usize) -> Float {
+        let theta_o = b.cos_theta_o.acos();
+        let theta_e = b.cos_theta_e.acos();
+        let theta_w = (theta_o + theta_e).min(PI);
+        let sin_o = safe_sqrt(1. - square(b.cos_theta_o));
+        let m_omega = 2. * PI * (1. - b.cos_theta_o)
+            + PI / 2.
+                * (2. * theta_w - (theta_o - 2. * theta_w).cos() - 2. * theta_o * sin_o
+                    + b.cos_theta_o);
+        let kr = bounds.diagonal().max_component_value() / bounds.diagonal()[dim];
+        b.phi * m_omega * kr * b.bounds.surface_area()
+    }
+}
+
+impl LightSamplerTrait for BVHLightSampler {
+    fn sample_with_context(&self, ctx: &LightSampleContext, mut u: Float) -> Option<SampledLight> {
+        let empty_nodes = if self.nodes.is_empty() { 0. } else { 1. };
+        let inf_size = self.infinite_lights.len() as Float;
+        let light_size = self.lights.len() as Float;
+
+        let p_inf = inf_size / (inf_size + empty_nodes);
+
+        if u < p_inf {
+            u /= p_inf;
+            let ind = (u * light_size).min(light_size - 1.) as usize;
+            let pmf = p_inf / inf_size;
+            Some(SampledLight::new(self.infinite_lights[ind].clone(), pmf))
+        } else {
+            if self.nodes.is_empty() {
+                return None;
+            }
+            let p = ctx.p();
+            let n = ctx.ns;
+            u = ((u - p_inf) / (1. - p_inf)).min(ONE_MINUS_EPSILON);
+            let mut node_ind = 0;
+            let mut pmf = 1. - p_inf;
+
+            loop {
+                let node = self.nodes[node_ind];
+                if !node.is_leaf() {
+                    let children: [LightBVHNode; 2] = [
+                        self.nodes[node_ind + 1],
+                        self.nodes[node.child_or_light_index() as usize],
+                    ];
+                    let ci: [Float; 2] = [
+                        children[0]
+                            .light_bounds
+                            .importance(p, n, &self.all_light_bounds),
+                        children[1]
+                            .light_bounds
+                            .importance(p, n, &self.all_light_bounds),
+                    ];
+
+                    if ci[0] == 0. && ci[1] == 0. {
+                        return None;
+                    }
+
+                    let mut node_pmf: Float = 0.;
+                    let child = sample_discrete(&ci, u, Some(&mut node_pmf), Some(&mut u))?;
+                    pmf *= node_pmf;
+                    node_ind = if child == 0 {
+                        node_ind + 1
+                    } else {
+                        node.child_or_light_index() as usize
+                    };
+                } else {
+                    if node_ind > 0
+                        || node.light_bounds.importance(p, n, &self.all_light_bounds) > 0.
+                    {
+                        return Some(SampledLight::new(
+                            self.lights[node.child_or_light_index() as usize].clone(),
+                            pmf,
+                        ));
+                    }
+                    return None;
+                }
+            }
+        }
+    }
+
+    fn pmf_with_context(&self, ctx: &LightSampleContext, light: &Arc<Light>) -> Float {
+        let ptr = Arc::as_ptr(light) as usize;
+        let empty_nodes = if self.nodes.is_empty() { 0. } else { 1. };
+        if self.light_to_bit_trail.contains_key(&ptr) {
+            return 1. / (self.infinite_lights.len() as Float + empty_nodes);
+        }
+
+        let mut bit_trail = self.light_to_bit_trail[&ptr];
+        let p = ctx.p();
+        let n = ctx.ns;
+        let p_inf = self.infinite_lights.len() as Float
+            / (self.infinite_lights.len() as Float + empty_nodes);
+        let mut pmf = 1. - p_inf;
+        let mut node_ind = 0;
+
+        loop {
+            let node = self.nodes[node_ind];
+            if node.is_leaf() {
+                return pmf;
+            }
+            let child0 = self.nodes[node_ind + 1];
+            let child1 = self.nodes[node.child_or_light_index() as usize];
+            let ci = [
+                child0.light_bounds.importance(p, n, &self.all_light_bounds),
+                child1.light_bounds.importance(p, n, &self.all_light_bounds),
+            ];
+            pmf *= ci[bit_trail & 1] / (ci[0] + ci[1]);
+            node_ind = if (bit_trail & 1) != 0 {
+                node.child_or_light_index() as usize
+            } else {
+                node_ind + 1
+            };
+            bit_trail >>= 1;
+        }
+    }
+
+    fn sample(&self, u: Float) -> Option<SampledLight> {
+        if self.lights.is_empty() {
+            return None;
+        }
+
+        let light_ind =
+            (u * self.lights.len() as Float).min(self.lights.len() as Float - 1.) as usize;
+
+        Some(SampledLight::new(
+            self.lights[light_ind].clone(),
+            1. / self.lights.len() as Float,
+        ))
+    }
+
+    fn pmf(&self, _light: &Arc<Light>) -> Float {
+        if self.lights.is_empty() {
+            return 0.;
+        }
+
+        1. / self.lights.len() as Float
+    }
+}