Add iridescent material

ray-tracing-material/src/iridescent.rs

@@ -0,0 +1,207 @@
+use core::f32;
+
+use ray_tracing_core::prelude::*;
+
+#[derive(Debug)]
+pub struct Iridescent {
+    d1: f32,
+    d2: f32,
+    n1: f32,
+    n2: f32,
+    n: f32,
+}
+
+impl Iridescent {
+    pub fn new(d1: f32, d2: f32, n1: f32, n2: f32, n: u32) -> Self {
+        Self {
+            d1,
+            d2,
+            n1,
+            n2,
+            n: n as f32,
+        }
+    }
+
+    fn fresnel_reflection(&self, theta1: f32, theta2: f32, s_polaized: bool) -> f32 {
+        if s_polaized {
+            (self.n1 * theta1.cos() - self.n2 * theta2.cos())
+                / (self.n1 * theta1.cos() + self.n2 * theta2.cos())
+        } else {
+            (self.n2 * theta1.cos() - self.n1 * theta2.cos())
+                / (self.n2 * theta1.cos() + self.n1 * theta2.cos())
+        }
+    }
+
+    fn abs_r_1_squared(&self, l: f32, theta1: f32, theta2: f32, s_polaized: bool) -> f32 {
+        let r12 = self.fresnel_reflection(theta1, theta2, s_polaized);
+
+        let phi = 2.0 * f32::consts::PI * self.n2 * self.d2 * theta2.cos() / l;
+
+        let num_real = r12 * (1.0 - f32::cos(-2.0 * phi));
+        let num_imag = r12 * (1.0 - f32::sin(-2.0 * phi));
+
+        let denom_real = 1.0 - r12 * r12 * f32::cos(-2.0 * phi);
+        let denom_imag = 1.0 - r12 * r12 * f32::sin(-2.0 * phi);
+
+        let real = (num_real * denom_real + num_imag * denom_imag)
+            / (denom_real * denom_real + denom_imag * denom_imag);
+
+        let imag = (num_imag * denom_real - num_real * denom_imag)
+            / (denom_real * denom_real + denom_imag * denom_imag);
+
+        real * real + imag * imag
+    }
+
+    fn abs_C_squared(&self, l: f32, theta1: f32, theta2: f32, s_polaized: bool) -> f32 {
+        let local_abs_r_1_squared = self.abs_r_1_squared(l, theta1, theta2, s_polaized);
+
+        local_abs_r_1_squared / (local_abs_r_1_squared + 1.0)
+    }
+
+    fn cos_K_Delta(&self, l: f32, theta1: f32, theta2: f32, s_polaized: bool) -> f32 {
+        let k1z = 2.0 * f32::consts::PI * self.n1 * theta1.cos() / l;
+        let k2z = 2.0 * f32::consts::PI * self.n2 * theta2.cos() / l;
+
+        let omega = if s_polaized {
+            k2z / k1z + k1z / k2z
+        } else {
+            (self.n1 * self.n1 * k2z) / (self.n2 * self.n2 * k1z)
+                + (self.n2 * self.n2 * k1z) / (self.n1 * self.n1 * k2z)
+        };
+
+        f32::cos(k1z * self.d1) * f32::cos(k2z * self.d2)
+            - 0.5 * omega * f32::sin(k1z * self.d1) * f32::sin(k2z * self.d2)
+    }
+
+    fn reflectance(&self, l: f32, theta1: f32, s_polaized: bool) -> f32 {
+        let theta2 = f32::asin(f32::sin(theta1) * self.n1 / self.n2);
+
+        let local_cos_k_delta = self.cos_K_Delta(l, theta1, theta2, s_polaized);
+
+        let local_c_squred = self.abs_C_squared(l, theta1, theta2, s_polaized);
+
+        if !local_c_squred.is_normal() || !local_cos_k_delta.is_normal() {
+            dbg!((
+                l,
+                theta1,
+                theta2,
+                local_cos_k_delta,
+                local_c_squred,
+                s_polaized
+            ));
+        }
+
+        if local_cos_k_delta.abs() <= 1.0 {
+            let k_delta = f32::acos(local_cos_k_delta);
+
+            if !k_delta.is_normal() {
+                dbg!((
+                    l,
+                    theta1,
+                    theta2,
+                    local_cos_k_delta,
+                    local_c_squred,
+                    s_polaized,
+                    k_delta,
+                ));
+            }
+            let u = f32::sin(k_delta) / f32::sin(self.n * k_delta);
+
+            local_c_squred / (local_c_squred + u * u)
+        } else {
+            let imk_delta = -f32::ln(f32::abs(
+                local_cos_k_delta - f32::sqrt(local_cos_k_delta * local_cos_k_delta - 1.0),
+            ));
+
+            if !imk_delta.is_normal() {
+                dbg!((
+                    l,
+                    theta1,
+                    theta2,
+                    local_cos_k_delta,
+                    local_c_squred,
+                    s_polaized,
+                    imk_delta,
+                ));
+            }
+            let u = f32::sinh(imk_delta) / f32::sinh(self.n * imk_delta);
+
+            local_c_squred / (local_c_squred + u * u)
+        }
+    }
+
+    fn averaged_reflectance(&self, l: f32, theta1: f32) -> f32 {
+        0.5 * (self.reflectance(l, theta1, true) + self.reflectance(l, theta1, false))
+    }
+
+    fn piecewise_gaussian(l: f32, mu: f32, tau1: f32, tau2: f32) -> f32 {
+        if l < mu {
+            let s = tau1 * (l - mu);
+            f32::exp(-0.5 * s * s)
+        } else {
+            let s = tau2 * (l - mu);
+            f32::exp(-0.5 * s * s)
+        }
+    }
+
+    fn color_matching_function(l: f32) -> Color {
+        Color::new(
+            1.056 * Self::piecewise_gaussian(l, 599.8, 0.0264, 0.0323)
+                + 0.362 * Self::piecewise_gaussian(l, 442.0, 0.0624, 0.0374)
+                - 0.065 * Self::piecewise_gaussian(l, 501.1, 0.0490, 0.0382),
+            0.821 * Self::piecewise_gaussian(l, 568.8, 0.0213, 0.0247)
+                + 0.286 * Self::piecewise_gaussian(l, 530.9, 0.0613, 0.0322),
+            1.217 * Self::piecewise_gaussian(l, 437.0, 0.0845, 0.0278)
+                + 0.681 * Self::piecewise_gaussian(l, 459.0, 0.0385, 0.0725),
+        )
+    }
+
+    fn monte_carlo_reflectance<R: Rng>(&self, theta1: f32, rng: &mut R) -> Color {
+        let range = 400.0..800.0;
+        let size = range.end - range.start;
+
+        let l = rng.gen_range(range);
+
+        let r = self.averaged_reflectance(l, theta1);
+
+        if !r.is_normal() {
+            dbg!((theta1, l, r));
+        }
+
+        let color = Self::color_matching_function(l);
+
+        if !color.r().is_normal() || !color.g().is_normal() || !color.b().is_normal() {
+            dbg!(color);
+        }
+
+        color * r
+    }
+}
+
+impl<R: Rng> Material<R> for Iridescent {
+    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_in;
+        let _ = w_out;
+        let _ = rng;
+        Color::black()
+    }
+
+    fn sample(&self, w_in: Dir3, rng: &mut R) -> ray_tracing_core::material::SampleResult {
+        let w_out = (2.0 * Dir3::up() * w_in.y()) - w_in;
+
+        let color = self.monte_carlo_reflectance(f32::acos(w_out.z()), rng);
+
+        ray_tracing_core::material::SampleResult::new(w_out, color, true)
+    }
+
+    fn pdf(&self, w_in: Dir3, w_out: Dir3) -> Float {
+        let _ = w_out;
+        let _ = w_in;
+        0.0
+    }
+}

ray-tracing-material/src/lib.rs

@@ -1,5 +1,6 @@
 pub mod default;
 pub mod diffuse;
+pub mod iridescent;
 pub mod microfacet;
 pub mod mirror;
 pub mod oren_nayar;