Add next event estimation

ray-tracing-core/src/math/mod.rs

@@ -2,6 +2,7 @@ pub mod dir3;
 pub mod frame;
 pub mod mat3;
 pub mod pos3;
+pub mod sampling;
 
 pub use dir3::Dir3;
 pub use frame::Frame;

ray-tracing-core/src/math/pos3.rs

@@ -31,6 +31,20 @@ impl Pos3 {
             z: 0.0,
         }
     }
+
+    pub fn from_barycentric(vertices: [Self; 3], coordinates: [Float; 3]) -> Self {
+        Self {
+            x: vertices[0].x() * coordinates[0]
+                + vertices[1].x() * coordinates[1]
+                + vertices[2].x() * coordinates[2],
+            y: vertices[0].y() * coordinates[0]
+                + vertices[1].y() * coordinates[1]
+                + vertices[2].y() * coordinates[2],
+            z: vertices[0].z() * coordinates[0]
+                + vertices[1].z() * coordinates[1]
+                + vertices[2].z() * coordinates[2],
+        }
+    }
 }
 
 impl Sub for Pos3 {

ray-tracing-core/src/math/sampling.rs

@@ -0,0 +1,8 @@
+use crate::prelude::*;
+
+/// Sample a triangle
+pub fn sample_triangle(uv: [Float; 2]) -> [Float; 3] {
+    let s = 1.0 - Float::sqrt(1.0 - uv[0]);
+    let t = (1.0 - s) * uv[1];
+    [s, t, 1.0 - s - t]
+}

ray-tracing-core/src/scene.rs

@@ -2,6 +2,13 @@ use crate::prelude::*;
 
 pub trait Scene<R: Rng> {
     fn intersect(&self, ray: Ray, min: Float, max: Float) -> Option<Intersection<'_, R>>;
+
+    fn sample_light(
+        &self,
+        w_in: Dir3,
+        intersection: &Intersection<'_, R>,
+        rng: &mut R,
+    ) -> Option<LightSample<'_, R>>;
 }
 
 pub struct Intersection<'sc, R: Rng> {
@@ -46,3 +53,41 @@ impl<'sc, R: Rng> Intersection<'sc, R> {
         Frame::from_normal(self.normal)
     }
 }
+
+pub struct LightSample<'sc, R: Rng> {
+    pos: Pos3,
+    pdf: Float,
+    normal: Dir3,
+    light: &'sc dyn Light<R>,
+}
+
+impl<'sc, R: Rng> LightSample<'sc, R> {
+    pub fn new(pos: Pos3, pdf: Float, normal: Dir3, light: &'sc dyn Light<R>) -> Self {
+        Self {
+            pos,
+            pdf,
+            normal,
+            light,
+        }
+    }
+
+    pub fn pdf(&self) -> Float {
+        self.pdf
+    }
+
+    pub fn pos(&self) -> Pos3 {
+        self.pos
+    }
+
+    pub fn light(&self) -> &'sc dyn Light<R> {
+        self.light
+    }
+
+    pub fn normal(&self) -> Dir3 {
+        self.normal
+    }
+
+    pub fn tangent_frame(&self) -> Frame {
+        Frame::from_normal(self.normal)
+    }
+}

ray-tracing-material/src/lib.rs

@@ -1,4 +1,4 @@
 pub mod default;
 pub mod diffuse;
-
+pub mod microfacet;
 pub mod mirror;

ray-tracing-material/src/microfacet.rs

@@ -0,0 +1,64 @@
+use ray_tracing_core::{material::Material, prelude::*};
+use std::fmt::Debug;
+
+#[derive(Debug)]
+struct Microfacet<D: MicrofacetDistribution + Debug> {
+    dist: D,
+    color: Color,
+}
+
+fn fresnel(w_in: Dir3, n: Dir3, nu_rel: Float) -> Float {
+    todo!()
+}
+
+impl<R: Rng, D: MicrofacetDistribution + Debug + Sync> Material<R> for Microfacet<D> {
+    fn eval(
+        &self,
+        w_in: ray_tracing_core::prelude::Dir3,
+        w_out: ray_tracing_core::prelude::Dir3,
+        _rng: &mut R,
+    ) -> ray_tracing_core::prelude::Color {
+        let w_h = (w_in + w_out).normalize();
+
+        let g = self.dist.g1(w_in, w_h) * self.dist.g1(w_out, w_h);
+
+        self.color * fresnel(w_in, w_h, 1.0) * g * self.dist.d(w_h) / (4.0 * w_in.y() * w_out.y())
+    }
+}
+
+struct BeckmannDistribution {
+    alpha: Float,
+}
+
+impl MicrofacetDistribution for BeckmannDistribution {
+    fn g1(&self, w: Dir3, w_h: Dir3) -> Float {
+        if w.y() > 0.0 && Dir3::dot(w, w_h) > 0.0 {
+            let cos_theta = w.y();
+            let a = self.alpha * (Float::sqrt(1.0 - cos_theta * cos_theta) / cos_theta);
+            (3.535 * a + 2.181 * a * a) / (1.0 + 2.276 * a + 2.577 * a * a)
+        } else {
+            1.0
+        }
+    }
+
+    fn d(&self, w_h: Dir3) -> Float {
+        if w_h.y() > 0.0 {
+            let cos_theta_squared = w_h.y() * w_h.y();
+            1.0 / (FloatConsts::PI
+                * self.alpha
+                * self.alpha
+                * cos_theta_squared
+                * cos_theta_squared)
+        } else {
+            0.0
+        }
+    }
+}
+
+trait MicrofacetDistribution {
+    fn g1(&self, w: Dir3, w_h: Dir3) -> Float;
+
+    fn d(&self, w_h: Dir3) -> Float;
+
+    // fn sample_d(&self) -> Dir3;
+}

ray-tracing-renderer/src/lib.rs

@@ -1,3 +1,4 @@
 pub mod depth_renderer;
+pub mod next_event_estimation;
 pub mod path_tracer;
 pub mod path_tracer_importance;

ray-tracing-renderer/src/next_event_estimation.rs

@@ -0,0 +1,102 @@
+use ray_tracing_core::{camera::Camera, prelude::*, renderer::ClassicalRenderer, scene::Scene};
+use std::marker::PhantomData;
+
+pub struct NextEventEstimation<R>
+where
+    R: Rng,
+{
+    width: u32,
+    height: u32,
+    rng: PhantomData<R>,
+}
+
+impl<R> NextEventEstimation<R>
+where
+    R: Rng,
+{
+    pub fn new(width: u32, height: u32) -> Self {
+        Self {
+            width,
+            height,
+            rng: PhantomData {},
+        }
+    }
+}
+
+impl<S, C, R> ClassicalRenderer<R, S, C> for NextEventEstimation<R>
+where
+    S: Scene<R>,
+    C: Camera<R>,
+    R: Rng,
+{
+    fn render_pixel(&self, scene: &S, camera: &C, x: u32, y: u32, rng: &mut R) -> Color {
+        let mut sum = Color::black();
+        let mut alpha = Color::white();
+
+        let mut r = camera.forward(x, y, rng);
+
+        let mut count = 0;
+
+        while let Some(i) = scene.intersect(r, 0.001, Float::INFINITY) {
+            let frame = i.tangent_frame();
+
+            let intersect_pos = r.at(i.t());
+
+            let w_in = frame.to_frame(-r.dir());
+
+            if let Some(light) = i.light() {
+                if count == 0 {
+                    sum += alpha * light.emit(w_in, rng);
+                }
+            }
+
+            let w_out = if let Some(material) = i.material() {
+                if let Some(l) = scene.sample_light(w_in, &i, rng) {
+                    let light_frame = l.tangent_frame();
+
+                    let light_dir = l.pos() - intersect_pos;
+
+                    let light_ray = Ray::new(intersect_pos, light_dir.normalize(), r.time());
+
+                    if scene
+                        .intersect(light_ray, 0.001, light_dir.length() - 0.001)
+                        .is_none()
+                    {
+                        let g = (Dir3::dot(i.normal(), light_dir.normalize())
+                            * Dir3::dot(l.normal(), -light_dir.normalize()))
+                            / light_dir.length_squared();
+
+                        sum += alpha
+                            * g
+                            * material.eval(w_in, frame.to_frame(light_dir.normalize()), rng)
+                            * l.light()
+                                .emit(light_frame.to_frame(-light_dir.normalize()), rng)
+                            / l.pdf();
+                    }
+                }
+
+                let sample_result = material.sample(w_in, rng);
+                alpha *= sample_result.color();
+                sample_result.w_out()
+            } else {
+                return sum;
+            };
+
+            r = Ray::new(r.at(i.t()), frame.to_world(w_out), r.time());
+            count += 1;
+            if count > 4 {
+                break;
+            }
+        }
+
+        sum
+    }
+
+    fn width(&self) -> u32 {
+        self.width
+    }
+
+    fn height(&self) -> u32 {
+        self.height
+    }
+}

ray-tracing-scene/Cargo.toml

@@ -6,5 +6,6 @@ edition = "2021"
 [dependencies]
 miette = { version = "7.2.0", features = ["fancy"] }
 nom = "7.1.3"
+rand = "0.8.5"
 ray-tracing-core = { path = "../ray-tracing-core" }
 ray-tracing-material = { path = "../ray-tracing-material" }

ray-tracing-scene/src/basic_scene.rs

@@ -52,4 +52,13 @@ impl<R: Rng> Scene<R> for BasicScene {
 
         intersection
     }
+
+    fn sample_light(
+        &self,
+        _w_in: Dir3,
+        _intersection: &Intersection<'_, R>,
+        _rng: &mut R,
+    ) -> Option<ray_tracing_core::scene::LightSample<'_, R>> {
+        None
+    }
 }

ray-tracing-scene/src/examples.rs

@@ -44,7 +44,7 @@ pub fn cornell2<R: Rng + Debug>() -> ExampleScene<R> {
         })
         .collect();
 
-    let scene = dbg!(TriangleBVH::new(obj.vertices, triangles, materials));
+    let scene = TriangleBVH::new(obj.vertices, triangles, materials);
 
     ExampleScene {
         scene,

ray-tracing-scene/src/triangle_bvh.rs

@@ -1,7 +1,9 @@
+use rand::seq::SliceRandom;
 use ray_tracing_core::{
     prelude::*,
-    scene::{Intersection, Scene},
+    scene::{Intersection, LightSample, Scene},
 };
+use sampling::sample_triangle;
 
 type Index = u32;
 
@@ -11,6 +13,7 @@ pub struct TriangleBVH<R: Rng> {
     triangles: Vec<Triangle>,
     materials: Vec<BVHMaterial<R>>,
     bvh: Vec<Node>,
+    lights: Vec<Triangle>,
 }
 
 #[derive(Debug)]
@@ -211,20 +214,29 @@ impl<R: Rng> TriangleBVH<R> {
         }];
         let aabb = calculate_aabb(&vertices, &triangles);
         build_bvh(&vertices, &mut triangles, &mut bvh, 0, aabb);
+        let lights: Vec<_> = triangles
+            .iter()
+            .filter(|f| materials[f.material as usize].light.is_some())
+            .map(Triangle::clone)
+            .collect();
         Self {
             vertices,
             bvh,
             triangles,
             materials,
+            lights,
         }
     }
 
-    fn get_vertices(&self, triangle: Index) -> [Pos3; 3] {
+    fn get_vertices_from_index(&self, triangle: Index) -> [Pos3; 3] {
-        let t = self.triangles[triangle as usize];
+        self.get_vertices(&self.triangles[triangle as usize])
+    }
+
+    fn get_vertices(&self, triangle: &Triangle) -> [Pos3; 3] {
         [
-            self.vertices[t.vertices[0] as usize],
+            self.vertices[triangle.vertices[0] as usize],
-            self.vertices[t.vertices[1] as usize],
+            self.vertices[triangle.vertices[1] as usize],
-            self.vertices[t.vertices[2] as usize],
+            self.vertices[triangle.vertices[2] as usize],
         ]
     }
 
@@ -278,7 +290,7 @@ impl<R: Rng> TriangleBVH<R> {
             Node::Leaf { start, count } => {
                 let mut intersection = None;
                 for i in start..(start + count) {
-                    if let Some(t) = triangle_intersection(ray, self.get_vertices(i)) {
+                    if let Some(t) = triangle_intersection(ray, self.get_vertices_from_index(i)) {
                         if min <= t && t <= max && !intersection.is_some_and(|(_, old_t)| t > old_t)
                         {
                             intersection.replace((i, t));
@@ -304,7 +316,7 @@ impl<R: Rng> Scene<R> for TriangleBVH<R> {
         let triangle = self.triangles[i as usize];
         let material = &self.materials[triangle.material as usize];
 
-        let n = triangle_normal(self.get_vertices(i)).normalize();
+        let n = triangle_normal(self.get_vertices_from_index(i)).normalize();
 
         Some(Intersection::new(
             t,
@@ -313,4 +325,34 @@ impl<R: Rng> Scene<R> for TriangleBVH<R> {
             material.light.as_deref(),
         ))
     }
+
+    fn sample_light(
+        &self,
+        _w_in: Dir3,
+        _intersection: &Intersection<'_, R>,
+        rng: &mut R,
+    ) -> Option<ray_tracing_core::scene::LightSample<'_, R>> {
+        let t = self.lights.choose(rng);
+
+        if let Some(t) = t {
+            let light = self.materials[t.material as usize]
+                .light
+                .as_ref()
+                .unwrap()
+                .as_ref();
+            let b = sample_triangle(rng.gen());
+
+            let n = triangle_normal(self.get_vertices(t));
+            let area = n.length() * 0.5;
+
+            Some(LightSample::new(
+                Pos3::from_barycentric(self.get_vertices(t), b),
+                1.0 / ((self.lights.len() as Float) * area),
+                n.normalize(),
+                light,
+            ))
+        } else {
+            None
+        }
+    }
 }

ray-tracing-tev/src/main.rs

@@ -7,8 +7,8 @@ use ray_tracing_core::{
     scene::Scene,
 };
 use ray_tracing_renderer::{
-    depth_renderer::DepthRenderer, path_tracer::PathTracer,
+    depth_renderer::DepthRenderer, next_event_estimation::NextEventEstimation,
-    path_tracer_importance::PathTracerImportance,
+    path_tracer::PathTracer, path_tracer_importance::PathTracerImportance,
 };
 use ray_tracing_scene::examples::example_scenes;
 use rayon::prelude::*;
@@ -17,14 +17,14 @@ use tev_client::{PacketCreateImage, PacketUpdateImage, TevClient, TevError};
 
 #[derive(Parser)]
 struct Args {
-    #[arg(default_value = "127.0.0.1:14158")]
+    #[arg(long, default_value = "127.0.0.1:14158")]
     tev: String,
     scenes: Vec<String>,
-    #[arg(default_value_t = 400)]
+    #[arg(long, default_value_t = 400)]
     width: u32,
-    #[arg(default_value_t = 400)]
+    #[arg(long, default_value_t = 400)]
     height: u32,
-    #[arg(default_value_t = 1024)]
+    #[arg(long, default_value_t = 1024)]
     samples_per_pixel: usize,
 }
 
@@ -115,20 +115,18 @@ fn main() {
             &r,
             &s,
             &c,
-            128,
+            args.samples_per_pixel,
             &mut client,
         )
         .unwrap();
 
-        let samples_per_pixel = 1024;
-
         let r = PathTracer::new(args.width, args.height);
         render_image(
             format!("{scene} - path tracer"),
             &r,
             &s,
             &c,
-            samples_per_pixel,
+            args.samples_per_pixel,
             &mut client,
         )
         .unwrap();
@@ -139,7 +137,18 @@ fn main() {
             &r,
             &s,
             &c,
-            samples_per_pixel,
+            args.samples_per_pixel,
+            &mut client,
+        )
+        .unwrap();
+
+        let r = NextEventEstimation::new(args.width, args.height);
+        render_image(
+            format!("{scene} - next event estimation"),
+            &r,
+            &s,
+            &c,
+            args.samples_per_pixel,
             &mut client,
         )
         .unwrap();