From 8459f3110f25d09034a57aacc63f2683f4e2a531 Mon Sep 17 00:00:00 2001
From: Matthew Gordon <matthew.scott.gordon@gmail.com>
Date: Mon, 4 Oct 2021 09:01:58 -0400
Subject: [PATCH] Add PDF to RandomDistribution

---
 src/random_distributions/mod.rs               |  1 +
 .../uniform_hemisphere.rs                     | 18 ++++++++++++++++++
 src/random_distributions/uniform_square.rs    | 19 +++++++++++++++++++
 src/random_distributions/unit_disc.rs         | 16 ++++++++++++++++
 4 files changed, 54 insertions(+)

diff --git a/src/random_distributions/mod.rs b/src/random_distributions/mod.rs
index 3e49b86..3686fe3 100644
--- a/src/random_distributions/mod.rs
+++ b/src/random_distributions/mod.rs
@@ -9,4 +9,5 @@ pub use uniform_hemisphere::UniformHemisphere;
 
 pub trait RandomDistribution<T> {
     fn value(&self) -> T;
+    fn pdf(&self, value: T) -> f64;
 }
diff --git a/src/random_distributions/uniform_hemisphere.rs b/src/random_distributions/uniform_hemisphere.rs
index 506f4f5..8fc7572 100644
--- a/src/random_distributions/uniform_hemisphere.rs
+++ b/src/random_distributions/uniform_hemisphere.rs
@@ -1,3 +1,5 @@
+use std::f64::consts::PI;
+
 use rand::distributions::{Open01, OpenClosed01};
 use rand::{thread_rng, Rng};
 
@@ -30,6 +32,10 @@ impl RandomDistribution<Vec3> for UniformHemisphere {
         }
         result.normalize()
     }
+
+    fn pdf(&self, _: Vec3) -> f64 {
+        1.0 / (2.0 * PI)
+    }
 }
 
 #[cfg(test)]
@@ -45,4 +51,16 @@ mod tests {
             println!("{}, {}, {}", value.x(), value.y(), value.z());
         }
     }
+
+    #[test]
+    #[ignore]
+    fn integral_is_near_area() {
+        let target = UniformHemisphere::new();
+        let integral = (0..1000)
+            .map(|_| target.value())
+            .map(|value| 1.0 / target.pdf(value))
+            .sum::<f64>()
+            / 1000.0;
+        println!("Area: {}\nIntegral: {}", 2.0 * PI, integral);
+    }
 }
diff --git a/src/random_distributions/uniform_square.rs b/src/random_distributions/uniform_square.rs
index 004f490..23d7ac9 100644
--- a/src/random_distributions/uniform_square.rs
+++ b/src/random_distributions/uniform_square.rs
@@ -23,6 +23,10 @@ impl RandomDistribution<Vec2> for UniformSquare {
         self.corner
             + Vec2::new(rng.sample::<f64, _>(Open01), rng.sample::<f64, _>(Open01)) * self.size
     }
+
+    fn pdf(&self, _value: Vec2) -> f64 {
+        1.0 / (self.size * self.size)
+    }
 }
 
 #[cfg(test)]
@@ -41,4 +45,19 @@ mod tests {
             println!("{}, {}", value.x(), value.y());
         }
     }
+
+    #[test]
+    #[ignore]
+    fn integral_is_near_area() {
+        let target = UniformSquare {
+            corner: Vec2::new(1.5, -2.5),
+            size: 3.0,
+        };
+        let integral = (0..1000)
+            .map(|_| target.value())
+            .map(|value| 1.0 / target.pdf(value))
+            .sum::<f64>()
+            / 1000.0;
+        println!("Area: {}\nIntegral: {}", 3.0 * 3.0, integral);
+    }
 }
diff --git a/src/random_distributions/unit_disc.rs b/src/random_distributions/unit_disc.rs
index f9a5057..7513b1c 100644
--- a/src/random_distributions/unit_disc.rs
+++ b/src/random_distributions/unit_disc.rs
@@ -32,6 +32,10 @@ impl RandomDistribution<Vec2> for UnitDisc {
             Vec2::new(angle.cos(), angle.sin()) * radius
         }
     }
+
+    fn pdf(&self, _: Vec2) -> f64 {
+        1.0 / PI
+    }
 }
 
 #[cfg(test)]
@@ -47,4 +51,16 @@ mod tests {
             println!("{}, {}", value.x(), value.y());
         }
     }
+
+    #[test]
+    #[ignore]
+    fn integral_is_near_area() {
+        let target = UnitDisc::new();
+        let integral = (0..1000)
+            .map(|_| target.value())
+            .map(|value| 1.0 / target.pdf(value))
+            .sum::<f64>()
+            / 1000.0;
+        println!("Area: {}\nIntegral: {}", PI, integral);
+    }
 }