Added GGX glossy material, and simplified surface closure API.

psychoblend/__init__.py

@@ -88,7 +88,7 @@ class PsychopathMesh(bpy.types.PropertyGroup):
 class PsychopathMaterial(bpy.types.PropertyGroup):
     surface_shader_type = EnumProperty(
         name="Surface Shader Type", description="",
-        items=[('Emit', 'Emit', ""), ('Lambert', 'Lambert', ""), ('GTR', 'GTR', "")],
+        items=[('Emit', 'Emit', ""), ('Lambert', 'Lambert', ""), ('GTR', 'GTR', ""), ('GGX', 'GGX', "")],
         default="Lambert"
         )
 

psychoblend/assembly.py

@@ -327,6 +327,12 @@ class Material:
             w.write("Roughness [%f]\n" % self.mat.psychopath.roughness)
             w.write("TailShape [%f]\n" % self.mat.psychopath.tail_shape)
             w.write("Fresnel [%f]\n" % self.mat.psychopath.fresnel)
+        elif self.mat.psychopath.surface_shader_type == 'GGX':
+            w.write("Type [GGX]\n")
+            color = self.mat.psychopath.color
+            w.write("Color [%f %f %f]\n" % (color[0], color[1], color[2]))
+            w.write("Roughness [%f]\n" % self.mat.psychopath.roughness)
+            w.write("Fresnel [%f]\n" % self.mat.psychopath.fresnel)
         else:
             raise "Unsupported surface shader type '%s'" % self.mat.psychopath.surface_shader_type
         w.unindent()

psychoblend/ui.py

@@ -249,6 +249,10 @@ class MATERIAL_PT_psychopath_surface(PsychopathPanel, bpy.types.Panel):
             layout.prop(mat.psychopath, "tail_shape")
             layout.prop(mat.psychopath, "fresnel")
 
+        if mat.psychopath.surface_shader_type == 'GGX':
+            layout.prop(mat.psychopath, "roughness")
+            layout.prop(mat.psychopath, "fresnel")
+
 
 def register():
     bpy.utils.register_class(RENDER_PT_psychopath_render_settings)

src/parse/psy_surface_shader.rs

@@ -160,6 +160,69 @@ pub fn parse_surface_shader<'a>(
                 fresnel: fresnel,
             })
         }
+
+        "GGX" => {
+            // Color
+            let color = if let Some((_, contents, byte_offset)) =
+                tree.iter_leaf_children_with_type("Color").nth(0)
+            {
+                if let IResult::Done(_, color) =
+                    closure!(tuple!(ws_f32, ws_f32, ws_f32))(contents.as_bytes())
+                {
+                    // TODO: handle color space conversions properly.
+                    // Probably will need a special color type with its
+                    // own parser...?
+                    XYZ::from_tuple(rec709_e_to_xyz(color))
+                } else {
+                    // Found color, but its contents is not in the right format
+                    return Err(PsyParseError::UnknownError(byte_offset));
+                }
+            } else {
+                return Err(PsyParseError::MissingNode(
+                    tree.byte_offset(),
+                    "Expected a Color field in GTR SurfaceShader.",
+                ));
+            };
+
+            // Roughness
+            let roughness = if let Some((_, contents, byte_offset)) =
+                tree.iter_leaf_children_with_type("Roughness").nth(0)
+            {
+                if let IResult::Done(_, roughness) = ws_f32(contents.as_bytes()) {
+                    roughness
+                } else {
+                    return Err(PsyParseError::UnknownError(byte_offset));
+                }
+            } else {
+                return Err(PsyParseError::MissingNode(
+                    tree.byte_offset(),
+                    "Expected a Roughness field in GTR SurfaceShader.",
+                ));
+            };
+
+            // Fresnel
+            let fresnel = if let Some((_, contents, byte_offset)) =
+                tree.iter_leaf_children_with_type("Fresnel").nth(0)
+            {
+                if let IResult::Done(_, fresnel) = ws_f32(contents.as_bytes()) {
+                    fresnel
+                } else {
+                    return Err(PsyParseError::UnknownError(byte_offset));
+                }
+            } else {
+                return Err(PsyParseError::MissingNode(
+                    tree.byte_offset(),
+                    "Expected a Fresnel field in GTR SurfaceShader.",
+                ));
+            };
+
+            arena.alloc(SimpleSurfaceShader::GGX {
+                color: color,
+                roughness: roughness,
+                fresnel: fresnel,
+            })
+        }
+
         _ => unimplemented!(),
     };
 

src/renderer.rs

@@ -492,10 +492,8 @@ impl LightPath {
 
                             // Distant light
                             SceneLightSample::Distant { direction, .. } => {
-                                let attenuation =
+                                let (attenuation, closure_pdf) =
                                     material.evaluate(ray.dir, direction, idata.nor, idata.nor_g);
-                                let closure_pdf =
-                                    material.sample_pdf(ray.dir, direction, idata.nor, idata.nor_g);
                                 let mut shadow_ray = {
                                     // Calculate the shadow ray for testing if the light is
                                     // in shadow or not.
@@ -520,10 +518,8 @@ impl LightPath {
                             // Surface light
                             SceneLightSample::Surface { sample_geo, .. } => {
                                 let dir = sample_geo.0 - idata.pos;
-                                let attenuation =
+                                let (attenuation, closure_pdf) =
                                     material.evaluate(ray.dir, dir, idata.nor, idata.nor_g);
-                                let closure_pdf =
-                                    material.sample_pdf(ray.dir, dir, idata.nor, idata.nor_g);
                                 let shadow_ray = {
                                     // Calculate the shadow ray for testing if the light is
                                     // in shadow or not.

src/shading/mod.rs

@@ -2,7 +2,7 @@ pub mod surface_closure;
 
 use std::fmt::Debug;
 
-use self::surface_closure::{EmitClosure, GTRClosure, LambertClosure, SurfaceClosureUnion};
+use self::surface_closure::{EmitClosure, GGXClosure, GTRClosure, LambertClosure, SurfaceClosureUnion};
 use color::{Color, XYZ};
 use surface::SurfaceIntersectionData;
 
@@ -43,6 +43,11 @@ pub enum SimpleSurfaceShader {
         tail_shape: f32,
         fresnel: f32,
     },
+    GGX {
+        color: XYZ,
+        roughness: f32,
+        fresnel: f32,
+    },
 }
 
 impl SurfaceShader for SimpleSurfaceShader {
@@ -72,6 +77,15 @@ impl SurfaceShader for SimpleSurfaceShader {
                 tail_shape,
                 fresnel,
             )),
+            SimpleSurfaceShader::GGX {
+                color,
+                roughness,
+                fresnel,
+            } => SurfaceClosureUnion::GGXClosure(GGXClosure::new(
+                color.to_spectral_sample(wavelength),
+                roughness,
+                fresnel,
+            )),
         }
     }
 }

src/shading/surface_closure.rs

@@ -15,6 +15,7 @@ pub enum SurfaceClosureUnion {
     EmitClosure(EmitClosure),
     LambertClosure(LambertClosure),
     GTRClosure(GTRClosure),
+    GGXClosure(GGXClosure),
 }
 
 impl SurfaceClosureUnion {
@@ -23,6 +24,7 @@ impl SurfaceClosureUnion {
             SurfaceClosureUnion::EmitClosure(ref closure) => closure as &SurfaceClosure,
             SurfaceClosureUnion::LambertClosure(ref closure) => closure as &SurfaceClosure,
             SurfaceClosureUnion::GTRClosure(ref closure) => closure as &SurfaceClosure,
+            SurfaceClosureUnion::GGXClosure(ref closure) => closure as &SurfaceClosure,
         }
     }
 }
@@ -32,6 +34,10 @@ impl SurfaceClosureUnion {
 /// Note: each surface closure is assumed to be bound to a particular hero
 /// wavelength.  This is implicit in the `sample`, `evaluate`, and `sample_pdf`
 /// functions below.
+///
+/// Also important is that _both_ the color filter and pdf returned from
+/// `sample()` and `evaluate()` should be identical for the same parameters
+/// and outgoing light direction.
 pub trait SurfaceClosure {
     /// Returns whether the closure has a delta distribution or not.
     fn is_delta(&self) -> bool;
@@ -59,19 +65,10 @@ pub trait SurfaceClosure {
     /// out:   The outgoing light direction.
     /// nor:   The shading surface normal at the surface point.
     /// nor_g: The geometric surface normal at the surface point.
-    /// wavelength: The wavelength of light to evaluate for.
     ///
-    /// Returns the resulting filter color.
+    /// Returns the resulting filter color and pdf of if this had been generated
-    fn evaluate(&self, inc: Vector, out: Vector, nor: Normal, nor_g: Normal) -> SpectralSample;
+    /// by `sample()`.
-
+    fn evaluate(&self, inc: Vector, out: Vector, nor: Normal, nor_g: Normal) -> (SpectralSample, f32);
-    /// Returns the pdf for the given 'in' direction producing the given 'out'
-    /// direction with the given differential geometry.
-    ///
-    /// inc: The incoming light direction.
-    /// out: The outgoing light direction.
-    /// nor:   The shading surface normal at the surface point.
-    /// nor_g: The geometric surface normal at the surface point.
-    fn sample_pdf(&self, inc: Vector, out: Vector, nor: Normal, nor_g: Normal) -> f32;
 
     /// Returns an estimate of the sum total energy that evaluate() would return
     /// when integrated over a spherical light source with a center at relative
@@ -194,16 +191,10 @@ impl SurfaceClosure for EmitClosure {
         (Vector::new(0.0, 0.0, 0.0), self.col, 1.0)
     }
 
-    fn evaluate(&self, inc: Vector, out: Vector, nor: Normal, nor_g: Normal) -> SpectralSample {
+    fn evaluate(&self, inc: Vector, out: Vector, nor: Normal, nor_g: Normal) -> (SpectralSample, f32) {
         let _ = (inc, out, nor, nor_g); // Not using these, silence warning
 
-        self.col
+        (self.col, 1.0)
-    }
-
-    fn sample_pdf(&self, inc: Vector, out: Vector, nor: Normal, nor_g: Normal) -> f32 {
-        let _ = (inc, out, nor, nor_g); // Not using these, silence warning
-
-        1.0
     }
 
     fn estimate_eval_over_sphere_light(
@@ -260,14 +251,13 @@ impl SurfaceClosure for LambertClosure {
 
         // Make sure it's not on the wrong side of the geometric normal.
         if dot(flipped_nor_g, out) >= 0.0 {
-            let filter = self.evaluate(inc, out, nor, nor_g);
+            (out, self.col * pdf, pdf)
-            (out, filter, pdf)
         } else {
             (out, SpectralSample::new(0.0), 0.0)
         }
     }
 
-    fn evaluate(&self, inc: Vector, out: Vector, nor: Normal, nor_g: Normal) -> SpectralSample {
+    fn evaluate(&self, inc: Vector, out: Vector, nor: Normal, nor_g: Normal) -> (SpectralSample, f32) {
         let (nn, flipped_nor_g) = if dot(nor_g.into_vector(), inc) <= 0.0 {
             (nor.normalized().into_vector(), nor_g.into_vector())
         } else {
@@ -276,23 +266,9 @@ impl SurfaceClosure for LambertClosure {
 
         if dot(flipped_nor_g, out) >= 0.0 {
             let fac = dot(nn, out.normalized()).max(0.0) * INV_PI;
-            self.col * fac
+            (self.col * fac, fac)
         } else {
-            SpectralSample::new(0.0)
+            (SpectralSample::new(0.0), 0.0)
-        }
-    }
-
-    fn sample_pdf(&self, inc: Vector, out: Vector, nor: Normal, nor_g: Normal) -> f32 {
-        let (nn, flipped_nor_g) = if dot(nor_g.into_vector(), inc) <= 0.0 {
-            (nor.normalized().into_vector(), nor_g.into_vector())
-        } else {
-            (-nor.normalized().into_vector(), -nor_g.into_vector())
-        };
-
-        if dot(flipped_nor_g, out) >= 0.0 {
-            dot(nn, out.normalized()).max(0.0) * INV_PI
-        } else {
-            0.0
         }
     }
 
@@ -356,6 +332,17 @@ impl SurfaceClosure for LambertClosure {
 
             let cos_nv = dot(nn, v).max(-1.0).min(1.0);
 
+            // Alt implementation from the SPI paper.
+            // Worse sampling, but here for reference.
+            // {
+            //     let nl_ang = cos_nv.acos();
+            //     let rad_ang = cos_theta_max.acos();
+            //     let min_ang = (nl_ang - rad_ang).max(0.0);
+            //     let lamb = min_ang.cos().max(0.0);
+
+            //     return lamb / dist2;
+            // }
+
             return sphere_lambert(cos_nv, cos_theta_max);
         }
     }
@@ -487,15 +474,14 @@ impl SurfaceClosure for GTRClosure {
 
         // Make sure it's not on the wrong side of the geometric normal.
         if dot(flipped_nor_g, out) >= 0.0 {
-            let filter = self.evaluate(inc, out, nor, nor_g);
+            let (filter, pdf) = self.evaluate(inc, out, nor, nor_g);
-            let pdf = self.sample_pdf(inc, out, nor, nor_g);
             (out, filter, pdf)
         } else {
             (out, SpectralSample::new(0.0), 0.0)
         }
     }
 
-    fn evaluate(&self, inc: Vector, out: Vector, nor: Normal, nor_g: Normal) -> SpectralSample {
+    fn evaluate(&self, inc: Vector, out: Vector, nor: Normal, nor_g: Normal) -> (SpectralSample, f32) {
         // Calculate needed vectors, normalized
         let aa = -inc.normalized(); // Vector pointing to where "in" came from
         let bb = out.normalized(); // Out
@@ -510,7 +496,7 @@ impl SurfaceClosure for GTRClosure {
 
         // Make sure everything's on the correct side of the surface
         if dot(nn, aa) < 0.0 || dot(nn, bb) < 0.0 || dot(flipped_nor_g, bb) < 0.0 {
-            return SpectralSample::new(0.0);
+            return (SpectralSample::new(0.0), 0.0);
         }
 
         // Calculate needed dot products
@@ -559,15 +545,10 @@ impl SurfaceClosure for GTRClosure {
         // Calculate everything else
         if self.roughness == 0.0 {
             // If sharp mirror, just return col * fresnel factor
-            return col_f;
+            return (col_f, 0.0);
         } else {
             // Calculate D - Distribution
-            let dist = if nh > 0.0 {
+            let dist = self.dist(nh, self.roughness);
-                let nh2 = nh * nh;
-                self.normalization_factor / (1.0 + ((roughness2 - 1.0) * nh2)).powf(self.tail_shape)
-            } else {
-                0.0
-            };
 
             // Calculate G1 - Geometric microfacet shadowing
             let g1 = {
@@ -590,34 +571,10 @@ impl SurfaceClosure for GTRClosure {
             };
 
             // Final result
-            col_f * (dist * g1 * g2) * INV_PI
+            (col_f * (dist * g1 * g2) * INV_PI, dist * INV_PI)
         }
     }
 
-    fn sample_pdf(&self, inc: Vector, out: Vector, nor: Normal, nor_g: Normal) -> f32 {
-        // Calculate needed vectors, normalized
-        let aa = -inc.normalized(); // Vector pointing to where "in" came from
-        let bb = out.normalized(); // Out
-        let hh = (aa + bb).normalized(); // Half-way between aa and bb
-
-        // Surface normal
-        let (nn, flipped_nor_g) = if dot(nor_g.into_vector(), inc) <= 0.0 {
-            (nor.normalized().into_vector(), nor_g.into_vector())
-        } else {
-            (-nor.normalized().into_vector(), -nor_g.into_vector())
-        };
-
-        // Make sure everything's on the correct side of the surface
-        if dot(nn, aa) < 0.0 || dot(nn, bb) < 0.0 || dot(flipped_nor_g, bb) < 0.0 {
-            return 0.0;
-        }
-
-        // Calculate needed dot products
-        let nh = clamp(dot(nn, hh), -1.0, 1.0);
-
-        self.dist(nh, self.roughness) * INV_PI
-    }
-
     fn estimate_eval_over_sphere_light(
         &self,
         inc: Vector,
@@ -676,3 +633,243 @@ impl SurfaceClosure for GTRClosure {
         fac * (1.0f32).min(1.0 - cos_theta_max) * INV_PI
     }
 }
+
+
+/// The GGX microfacet BRDF.
+#[derive(Debug, Copy, Clone)]
+pub struct GGXClosure {
+    col: SpectralSample,
+    roughness: f32,
+    fresnel: f32, // [0.0, 1.0] determines how much fresnel reflection comes into play
+}
+
+impl GGXClosure {
+    pub fn new(col: SpectralSample, roughness: f32, fresnel: f32) -> GGXClosure {
+        let mut closure = GGXClosure {
+            col: col,
+            roughness: roughness,
+            fresnel: fresnel,
+        };
+
+        closure.validate();
+
+        closure
+    }
+
+    // Makes sure values are in a valid range
+    fn validate(&mut self) {
+        debug_assert!(self.fresnel >= 0.0 && self.fresnel <= 1.0);
+        debug_assert!(self.roughness >= 0.0 && self.roughness <= 1.0);
+    }
+
+    // Returns the cosine of the half-angle that should be sampled, given
+    // a random variable in [0,1]
+    fn half_theta_sample(u: f32, rough: f32) -> f32 {
+        let rough2 = rough * rough;
+
+        // Calculate top half of equation
+        let top = 1.0 - u;
+
+        // Calculate bottom half of equation
+        let bottom = 1.0 + ((rough2 - 1.0) * u);
+
+        (top / bottom).sqrt()
+    }
+
+    /// The GGX microfacet distribution function.
+    ///
+    /// nh: cosine of the angle between the surface normal and the microfacet normal.
+    fn ggx_d(nh: f32, rough: f32) -> f32 {
+        if nh <= 0.0 {
+            return 0.0;
+        }
+
+        let rough2 = rough * rough;
+        let tmp = 1.0 + ((rough2 - 1.0) * (nh * nh));
+        rough2 / (PI_32 * tmp * tmp)
+    }
+
+    /// The GGX Smith shadow-masking function.
+    ///
+    /// vh: cosine of the angle between the view vector and the microfacet normal.
+    /// vn: cosine of the angle between the view vector and surface normal.
+    fn ggx_g(vh: f32, vn: f32, rough: f32) -> f32 {
+        if (vh * vn) <= 0.0 {
+            0.0
+        } else {
+            2.0 / (1.0 + (1.0 + rough * rough * (1.0 - vn * vn) / (vn * vn)).sqrt())
+        }
+    }
+}
+
+impl SurfaceClosure for GGXClosure {
+    fn is_delta(&self) -> bool {
+        self.roughness == 0.0
+    }
+
+    fn sample(
+        &self,
+        inc: Vector,
+        nor: Normal,
+        nor_g: Normal,
+        uv: (f32, f32),
+    ) -> (Vector, SpectralSample, f32) {
+        // Get normalized surface normal
+        let (nn, flipped_nor_g) = if dot(nor_g.into_vector(), inc) <= 0.0 {
+            (nor.normalized().into_vector(), nor_g.into_vector())
+        } else {
+            (-nor.normalized().into_vector(), -nor_g.into_vector())
+        };
+
+        // Generate a random ray direction in the hemisphere
+        // of the surface.
+        let theta_cos = Self::half_theta_sample(uv.0, self.roughness);
+        let theta_sin = (1.0 - (theta_cos * theta_cos)).sqrt();
+        let angle = uv.1 * PI_32 * 2.0;
+        let mut half_dir = Vector::new(angle.cos() * theta_sin, angle.sin() * theta_sin, theta_cos);
+        half_dir = zup_to_vec(half_dir, nn).normalized();
+
+        let out = inc - (half_dir * 2.0 * dot(inc, half_dir));
+
+        // Make sure it's not on the wrong side of the geometric normal.
+        if dot(flipped_nor_g, out) >= 0.0 {
+            let (filter, pdf) = self.evaluate(inc, out, nor, nor_g);
+            (out, filter, pdf)
+        } else {
+            (out, SpectralSample::new(0.0), 0.0)
+        }
+    }
+
+    fn evaluate(&self, inc: Vector, out: Vector, nor: Normal, nor_g: Normal) -> (SpectralSample, f32) {
+        // Calculate needed vectors, normalized
+        let aa = -inc.normalized(); // Vector pointing to where "in" came from
+        let bb = out.normalized(); // Out
+        let hh = (aa + bb).normalized(); // Half-way between aa and bb
+
+        // Surface normal
+        let (nn, flipped_nor_g) = if dot(nor_g.into_vector(), inc) <= 0.0 {
+            (nor.normalized().into_vector(), nor_g.into_vector())
+        } else {
+            (-nor.normalized().into_vector(), -nor_g.into_vector())
+        };
+
+        // Make sure everything's on the correct side of the surface
+        if dot(nn, aa) < 0.0 || dot(nn, bb) < 0.0 || dot(flipped_nor_g, bb) < 0.0 {
+            return (SpectralSample::new(0.0), 0.0);
+        }
+
+        // Calculate needed dot products
+        let na = clamp(dot(nn, aa), -1.0, 1.0);
+        let nb = clamp(dot(nn, bb), -1.0, 1.0);
+        let ha = clamp(dot(hh, aa), -1.0, 1.0);
+        let hb = clamp(dot(hh, bb), -1.0, 1.0);
+        let nh = clamp(dot(nn, hh), -1.0, 1.0);
+
+        // Calculate F - Fresnel
+        let col_f = {
+            let rev_fresnel = 1.0 - self.fresnel;
+            let c0 = lerp(
+                schlick_fresnel_from_fac(self.col.e.get_0(), hb),
+                self.col.e.get_0(),
+                rev_fresnel,
+            );
+            let c1 = lerp(
+                schlick_fresnel_from_fac(self.col.e.get_1(), hb),
+                self.col.e.get_1(),
+                rev_fresnel,
+            );
+            let c2 = lerp(
+                schlick_fresnel_from_fac(self.col.e.get_2(), hb),
+                self.col.e.get_2(),
+                rev_fresnel,
+            );
+            let c3 = lerp(
+                schlick_fresnel_from_fac(self.col.e.get_3(), hb),
+                self.col.e.get_3(),
+                rev_fresnel,
+            );
+
+            let mut col_f = self.col;
+            col_f.e.set_0(c0);
+            col_f.e.set_1(c1);
+            col_f.e.set_2(c2);
+            col_f.e.set_3(c3);
+
+            col_f
+        };
+
+        // Calculate everything else
+        if self.roughness == 0.0 {
+            // If sharp mirror, just return col * fresnel factor
+            return (col_f, 0.0);
+        } else {
+            // Calculate D - Distribution
+            let dist = Self::ggx_d(nh, self.roughness);
+
+            // Calculate G1 and G2- Geometric microfacet shadowing
+            let g1 = Self::ggx_g(ha, na, self.roughness);
+            let g2 = Self::ggx_g(hb, nb, self.roughness);
+
+            // Final result
+            (col_f * (dist * g1 * g2) * INV_PI, dist * INV_PI)
+        }
+    }
+
+    fn estimate_eval_over_sphere_light(
+        &self,
+        inc: Vector,
+        to_light_center: Vector,
+        light_radius_squared: f32,
+        nor: Normal,
+        nor_g: Normal,
+    ) -> f32 {
+        // TODO: all of the stuff in this function is horribly hacky.
+        // Find a proper way to approximate the light contribution from a
+        // solid angle.
+
+        let _ = nor_g; // Not using this, silence warning
+
+        let dist2 = to_light_center.length2();
+        let sin_theta_max2 = (light_radius_squared / dist2).min(1.0);
+        let cos_theta_max = (1.0 - sin_theta_max2).sqrt();
+
+        assert!(cos_theta_max >= -1.0);
+        assert!(cos_theta_max <= 1.0);
+
+        // Surface normal
+        let nn = if dot(nor.into_vector(), inc) < 0.0 {
+            nor.normalized()
+        } else {
+            -nor.normalized() // If back-facing, flip normal
+        }.into_vector();
+
+        let aa = -inc.normalized(); // Vector pointing to where "in" came from
+        let bb = to_light_center.normalized(); // Out
+
+        // Brute-force method
+        //let mut fac = 0.0;
+        //const N: usize = 256;
+        //for i in 0..N {
+        //    let uu = Halton::sample(0, i);
+        //    let vv = Halton::sample(1, i);
+        //    let mut samp = uniform_sample_cone(uu, vv, cos_theta_max);
+        //    samp = zup_to_vec(samp, bb).normalized();
+        //    if dot(nn, samp) > 0.0 {
+        //        let hh = (aa+samp).normalized();
+        //        fac += self.dist(dot(nn, hh), roughness);
+        //    }
+        //}
+        //fac /= N * N;
+
+        // Approximate method
+        let theta = cos_theta_max.acos();
+        let hh = (aa + bb).normalized();
+        let nh = clamp(dot(nn, hh), -1.0, 1.0);
+        let fac = Self::ggx_d(
+            nh,
+            (1.0f32).min(self.roughness.sqrt() + (2.0 * theta / PI_32)),
+        );
+
+        fac * (1.0f32).min(1.0 - cos_theta_max) * INV_PI
+    }
+}