Fix some Microfacet bugs

Cargo.lock

@@ -1382,6 +1382,12 @@ dependencies = [
  "windows-targets 0.52.6",
 ]
 
+[[package]]
+name = "libm"
+version = "0.2.11"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "8355be11b20d696c8f18f6cc018c4e372165b1fa8126cef092399c9951984ffa"
+
 [[package]]
 name = "libredox"
 version = "0.1.3"
@@ -1662,6 +1668,7 @@ source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "071dfc062690e90b734c0b2273ce72ad0ffa95f0c74596bc250dcfd960262841"
 dependencies = [
  "autocfg",
+ "libm",
 ]
 
 [[package]]
@@ -2105,6 +2112,16 @@ dependencies = [
  "getrandom",
 ]
 
+[[package]]
+name = "rand_distr"
+version = "0.4.3"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "32cb0b9bc82b0a0876c2dd994a7e7a2683d3e7390ca40e6886785ef0c7e3ee31"
+dependencies = [
+ "num-traits",
+ "rand",
+]
+
 [[package]]
 name = "rav1e"
 version = "0.7.1"
@@ -2201,6 +2218,7 @@ dependencies = [
 name = "ray-tracing-material"
 version = "0.1.0"
 dependencies = [
+ "rand_distr",
  "ray-tracing-core",
 ]
 

ray-tracing-material-visualizer/src/bin/fresnel.rs

@@ -0,0 +1,27 @@
+use plotters::prelude::*;
+use ray_tracing_core::prelude::*;
+use ray_tracing_material::microfacet::fresnel_real;
+
+fn main() -> Result<(), Box<dyn std::error::Error>> {
+    let root = BitMapBackend::new("fresnel.png", (640, 480)).into_drawing_area();
+    root.fill(&WHITE)?;
+    let mut chart = ChartBuilder::on(&root)
+        .caption("Fresnel", ("sans-serif", 50).into_font())
+        .margin(5)
+        .x_label_area_size(30)
+        .y_label_area_size(30)
+        .build_cartesian_2d(0.0..FloatConsts::FRAC_PI_2, 0.0..1.0f32)?;
+
+    chart.configure_mesh().draw()?;
+
+    chart.draw_series(LineSeries::new(
+        (0..=300)
+            .map(|x| x as Float * FloatConsts::FRAC_PI_2 / 300.0)
+            .map(|x| (x, fresnel_real(Float::cos(x), 1.0, 1.5))),
+        &RED,
+    ))?;
+
+    root.present()?;
+
+    Ok(())
+}

ray-tracing-material-visualizer/src/main.rs

@@ -1,3 +1,5 @@
+use core::f64;
+
 use plotters::{
     chart::ChartBuilder,
     coord::Shift,
@@ -19,17 +21,17 @@ fn main() -> Result<(), Box<dyn std::error::Error>> {
 
     let color = Color::new(1.0, 1.0, 1.0);
 
-    let m = Mirror::new(color);
+    // let m = Mirror::new(color);
-    generate_chart("mirror.png", &m, 2, w_in)?;
+    // generate_chart("mirror.png", &m, 2, w_in)?;
 
-    let m = DiffuseMaterial::new(color);
+    // let m = DiffuseMaterial::new(color);
-    generate_chart("diffuse.png", &m, 100, w_in)?;
+    // generate_chart("diffuse.png", &m, 100, w_in)?;
 
-    let m = Microfacet::new(BeckmannDistribution::new(0.5), color);
+    let m = Microfacet::new(BeckmannDistribution::new(0.1), color);
-    generate_chart("microfacet.png", &m, 100, w_in)?;
+    generate_chart("microfacet.png", &m, 400, w_in)?;
 
-    let m = OrenNayar::new(0.5, color);
+    // let m = OrenNayar::new(0.5, color);
-    generate_chart("oren-nayar.png", &m, 100, w_in)?;
+    // generate_chart("oren-nayar.png", &m, 100, w_in)?;
 
     Ok(())
 }
@@ -52,7 +54,10 @@ fn generate_chart<M: Material<SmallRng>>(
     let eval_histogram = create_histogram_evaled(m, 400, 400, executions, w_in);
     let sample_histogram = create_historgram_sampled(m, 400, 400, executions, w_in);
 
-    let max = Float::max(eval_histogram.max(), sample_histogram.max());
+    let max = Float::max(
+        Float::max(eval_histogram.max(), sample_histogram.max()),
+        Float::MIN_POSITIVE,
+    );
 
     let area = areas[0].titled("Evaled", ("sans-serif", 30))?;
     plot_material(&area, eval_histogram, max)?;
@@ -142,7 +147,7 @@ where
                 _ => unreachable!(),
             };
 
-            assert!(f >= 0.0 && f <= max && f.is_finite(), "{f}");
+            assert!(f >= 0.0 && f <= max && f.is_finite() && max != 0.0, "{f}");
             let color = plotters::prelude::ViridisRGB::get_color_normalized(f, 0.0, max);
             if dir.y() < 0.0 {
                 lower_area.draw_pixel((dir.x() as f64, dir.z() as f64), &color)?;
@@ -160,7 +165,7 @@ struct Histogram {
     inserted: usize,
     width: usize,
     height: usize,
-    max: f32,
+    max: f64,
 }
 
 impl Histogram {
@@ -194,14 +199,15 @@ impl Histogram {
     fn insert(&mut self, pos: Dir3, value: Color) {
         let index = self.discretize(pos);
         self.data[index] += value;
-        self.max = f32::max(self.max, self.data[index].r());
+        self.max = f64::max(self.max, self.data[index].r() as f64);
-        self.max = f32::max(self.max, self.data[index].g());
+        self.max = f64::max(self.max, self.data[index].g() as f64);
-        self.max = f32::max(self.max, self.data[index].b());
+        self.max = f64::max(self.max, self.data[index].b() as f64);
         self.inserted += 1;
     }
 
     fn max(&self) -> f32 {
-        self.max / (self.inserted as f32)
+        dbg!(self.max, self.inserted, self.inserted as f64);
+        (self.max / (self.inserted as f64)) as f32
     }
 
     fn iter(&self) -> impl Iterator<Item = (Dir3, Color)> + use<'_> {
@@ -280,8 +286,11 @@ fn create_historgram_sampled<M: Material<SmallRng>>(
                         && sample.color().g() >= 0.0
                         && sample.color().g().is_finite()
                         && sample.color().b() >= 0.0
-                        && sample.color().b().is_finite(),
+                        && sample.color().b().is_finite()
-                    "w_in: {w_in:?}; w_out: {:?}; material: {m:?}; color: {:?}",
+                        && (-1.0..=1.0).contains(&sample.w_out().x())
+                        && (-1.0..=1.0).contains(&sample.w_out().y())
+                        && (-1.0..=1.0).contains(&sample.w_out().z()),
+                    "invalid_sampling: w_in: {w_in:?}; w_out: {:?}; material: {m:?}; color: {:?}",
                     sample.w_out(),
                     sample.color()
                 );
@@ -317,7 +326,7 @@ fn create_histogram_evaled<M: Material<SmallRng>>(
                         && color.g().is_finite()
                         && color.b() >= 0.0
                         && color.b().is_finite(),
-                    "w_in: {w_in:?}; w_out: {w_out:?}; material: {m:?}; color: {color:?}"
+                    "invalid_eval: w_in: {w_in:?}; w_out: {w_out:?}; material: {m:?}; color: {color:?}"
                 );
 
                 (w_out, color)

ray-tracing-material/Cargo.toml

@@ -4,4 +4,5 @@ version = "0.1.0"
 edition = "2021"
 
 [dependencies]
+rand_distr = "0.4.3"
 ray-tracing-core = { path = "../ray-tracing-core" }

ray-tracing-material/src/microfacet.rs

@@ -1,3 +1,4 @@
+use rand_distr::{Distribution, Normal};
 use ray_tracing_core::{material::Material, prelude::*};
 use std::fmt::Debug;
 
@@ -13,7 +14,7 @@ impl<D: MicrofacetDistribution + Debug> Microfacet<D> {
     }
 }
 
-fn fresnel_real(cos_theta_in: Float, nu1: Float, nu2: Float) -> Float {
+pub fn fresnel_real(cos_theta_in: Float, nu1: Float, nu2: Float) -> Float {
     let nu_rel = nu1 / nu2;
     let cos_theta_tr = Float::sqrt(1.0 - nu_rel * nu_rel * (1.0 - cos_theta_in * cos_theta_in));
 
@@ -34,12 +35,28 @@ impl<R: Rng, D: MicrofacetDistribution + Debug + Sync> Material<R> for Microface
 
             let g = self.dist.g1(w_in, w_h) * self.dist.g1(w_out, w_h);
 
-            self.color * fresnel_real(Dir3::dot(w_in, w_h), 1.0, 1.3) * g * self.dist.d(w_h)
+            self.color * fresnel_real(Dir3::dot(w_in, w_h), 1.0, 1.5) * g * self.dist.d(w_h)
                 / (4.0 * w_in.y() * w_out.y()).max(0.0)
         } else {
             Color::black()
         }
     }
+
+    fn sample(&self, w_in: Dir3, rng: &mut R) -> ray_tracing_core::material::SampleResult {
+        let w_h = self.dist.sample_d(rng);
+
+        let w_out = (2.0 * (Dir3::dot(w_in, w_h) * w_h) - w_in).normalize();
+
+        if w_out.y() > 0.0 {
+            let g = self.dist.g1(w_in, w_h) * self.dist.g1(w_out, w_h);
+            ray_tracing_core::material::SampleResult::new(
+                w_out,
+                self.color * fresnel_real(Dir3::dot(w_in, w_h), 1.0, 1.5) * g,
+            )
+        } else {
+            ray_tracing_core::material::SampleResult::new(w_out, Color::black())
+        }
+    }
 }
 
 #[derive(Debug)]
@@ -57,25 +74,43 @@ 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);
+            let a = 1.0 / (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)
+            if a < 1.6 {
+                (3.535 * a + 2.181 * a * a) / (1.0 + 2.276 * a + 2.577 * a * a)
+            } else {
+                1.0
+            }
         } else {
-            1.0
+            0.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
+            let tan_theta_squared =
-                * self.alpha
+                w_h.x() * w_h.x() / cos_theta_squared + w_h.z() * w_h.z() / cos_theta_squared;
-                * self.alpha
+            Float::exp(-tan_theta_squared / (self.alpha * self.alpha))
-                * cos_theta_squared
+                / (FloatConsts::PI
-                * cos_theta_squared)
+                    * self.alpha
+                    * self.alpha
+                    * cos_theta_squared
+                    * cos_theta_squared)
         } else {
             0.0
         }
     }
+
+    fn sample_d<R: Rng>(&self, rng: &mut R) -> Dir3 {
+        let n = Normal::new(0.0, self.alpha).unwrap();
+
+        let mx = n.sample(rng);
+        let mz = n.sample(rng);
+
+        let y = Float::sqrt(1.0 / (mx * mx + 1.0 + mz * mz));
+
+        Dir3::new(mx * y, y, mz * y).normalize()
+    }
 }
 
 pub trait MicrofacetDistribution {
@@ -83,5 +118,5 @@ pub trait MicrofacetDistribution {
 
     fn d(&self, w_h: Dir3) -> Float;
 
-    // fn sample_d(&self) -> Dir3;
+    fn sample_d<R: Rng>(&self, rng: &mut R) -> Dir3;
 }