Add sphere type and ray-sphere intersection test

Author: jdahlstrom

core/src/geom/prim.rs

@@ -3,13 +3,12 @@
 //! Includes vertices, polygons, planes, rays, and more.
 
 use alloc::vec::Vec;
-use core::fmt::Debug;
+use core::fmt::{self, Debug, Formatter};
 
 use crate::math::{
     Affine, Lerp, Linear, Mat4, Parametric, Point, Point2, Point3, Vec2, Vec3,
-    Vector, space::Real, vec2, vec3,
+    Vector, pt3, space::Real, vec2, vec3,
 };
-
 use crate::render::Model;
 
 /// Vertex with a position and arbitrary other attributes.
@@ -42,6 +41,8 @@ pub type Plane3<B = ()> = Plane<Vector<[f32; 4], Real<3, B>>>;
 #[derive(Copy, Clone, Debug, Eq, PartialEq)]
 pub struct Ray<T: Affine>(pub T, pub T::Diff);
 
+pub type Ray3<B = ()> = Ray<Point3<B>>;
+
 /// A curve composed of a chain of line segments.
 ///
 /// The polyline is represented as a list of points, or vertices, with each
@@ -60,6 +61,9 @@ pub struct Polygon<T>(pub Vec<T>);
 #[derive(Copy, Clone, Debug, Default, Eq, PartialEq)]
 pub struct Edge<T>(pub T, pub T);
 
+#[derive(Copy, Clone, PartialEq)]
+pub struct Sphere<B = ()>(pub Point3<B>, pub f32);
+
 /// A surface normal in 3D.
 // TODO Use distinct type rather than alias
 pub type Normal3 = Vec3;
@@ -748,6 +752,22 @@ impl<B> Default for Plane3<B> {
     }
 }
 
+impl<B: Debug + Default> Debug for Sphere<B> {
+    fn fmt(&self, f: &mut Formatter<'_>) -> fmt::Result {
+        f.debug_tuple("Sphere")
+            .field(&self.0)
+            .field(&self.1)
+            .finish()
+    }
+}
+
+impl<B> Default for Sphere<B> {
+    /// Returns a unit sphere, with the center at the origin and radius 1.
+    fn default() -> Self {
+        Self(pt3(0.0, 0.0, 0.0), 1.0)
+    }
+}
+
 #[cfg(test)]
 mod tests {
     use crate::assert_approx_eq;

geom/src/isect.rs

@@ -1,7 +1,9 @@
 use core::fmt::Debug;
 
+#[cfg(feature = "std")] // TODO separate fp feature for geom
+use retrofire_core::geom::Sphere;
 use retrofire_core::{
-    geom::{Plane3, Ray},
+    geom::{Plane3, Ray, Ray3},
     math::{ApproxEq, Point3, vec3},
     render::scene::BBox,
 };
@@ -18,8 +20,10 @@ pub trait Intersect<T> {
     fn intersect(&self, other: &T) -> Self::Result;
 }
 
-impl<B> Intersect<Plane3<B>> for Ray<Point3<B>> {
-    type Result = Option<(f32, Point3<B>)>;
+type RayIntersect3<B> = Option<(f32, Point3<B>)>;
+
+impl<B> Intersect<Plane3<B>> for Ray3<B> {
+    type Result = RayIntersect3<B>;
 
     /// Returns the unique intersection point of `self` and a plane,
     /// or `None` if they do not intersect.
@@ -52,8 +56,8 @@ impl<B> Intersect<Plane3<B>> for Ray<Point3<B>> {
     }
 }
 
-impl<B: Debug + Default> Intersect<BBox<B>> for Ray<Point3<B>> {
-    type Result = Option<(f32, Point3<B>)>; // Only closest for now
+impl<B: Debug + Default> Intersect<BBox<B>> for Ray3<B> {
+    type Result = RayIntersect3<B>; // Only closest for now
 
     /// Returns the nearest intersection point of `self` and a box,
     /// or `None` if they do not intersect.
@@ -120,6 +124,93 @@ impl<B: Debug + Default> Intersect<BBox<B>> for Ray<Point3<B>> {
     }
 }
 
+#[cfg(feature = "std")]
+impl<B> Intersect<Sphere<B>> for Ray3<B> {
+    type Result = RayIntersect3<B>; // Only closest for now
+
+    /// Returns the intersection point of `self` and a sphere closest to the
+    /// origin of `self`, or `None` if they do not intersect.
+    ///
+    /// # Examples
+    /// ```
+    /// ```
+    fn intersect(&self, &Sphere(center, r): &Sphere<B>) -> Self::Result {
+        let &Ray(orig, dir) = self;
+
+        //             > r, no intersection
+        // If |C - C'| = r, one    -"-
+        //             < r, two    -"-
+        //           _______
+        //          /       \
+        //         /         \
+        //        |     C--r--|
+        //         \    |    /
+        //          \___|___/
+        //             _|
+        //   O--------+-C'----> d
+        //
+        // Find point P = (x, y, z) given sphere (C, r) and ray (O, d)
+        //
+        // Sphere equation:
+        //   (x - c_x)² + (y - c_y)² + (z - c_z)² = r²
+        // or in vector form
+        //   (P - C) · (P - C) = r²
+        //
+        // Ray equation:
+        //   P = O + t·d
+        //
+        // Intersection:
+        //
+        //   Substitute ray equation to sphere equation:
+        //   (o + t·d - c) · (o + t·d - c) = r²
+        //
+        //   Multiply out
+        //   o (o + td - c) + t·d (o + td - c) - c (o + td - c) = r²
+        //
+        //   Distribute
+        //   o·o + o·td - o·c + o·td + td·td - td·c - c·o - c·td + c·c = r²
+        //
+        //   Reorder
+        //   td·td + o·td + o·td - td·c + o·o - o·c - c·o + c·c - r² = 0
+        //
+        //   Factor out t's
+        //   (d·d) t² + (o·d + o·d - c·d) t + o·o - 2o·c + c·c - r² = 0
+        //
+        // Solve quadratic equation:
+        //   (d·d) t²  +  2(o - c)·d t  +  (o - c)² - r²  =  0
+        //
+        //   t = (-b ± √(b² - 4ac)) / 2a
+
+        let c_to_o = orig - center;
+        let a = dir.len_sqr(); // >= 0
+        let b = 2.0 * c_to_o.dot(&dir);
+        let c = c_to_o.len_sqr() - r * r;
+
+        let discriminant = b * b - 4.0 * a * c;
+        if discriminant < 0.0 {
+            // the line of the ray does not hit the sphere
+            return None;
+        }
+
+        use retrofire_core::math::float::f32;
+        let sqrt = f32::sqrt(discriminant);
+        // sqrt >= 0.0, thus t0 <= t1 always
+        let (t0, t1) = (-b - sqrt, -b + sqrt);
+        let t = if t0 >= 0.0 {
+            // ray hits both points
+            t0
+        } else if t1 >= 0.0 {
+            // ray origin is inside sphere
+            t1
+        } else {
+            // sphere is behind ray
+            return None;
+        };
+        let t = t / (2.0 * a);
+        Some((t, orig + t * dir))
+    }
+}
+
 #[cfg(test)]
 mod tests {
     use retrofire_core::math::{Linear, Vec3, pt3};
@@ -299,4 +390,46 @@ mod tests {
             assert_eq!(ray.intersect(&empty), None);
         }
     }
+
+    // TODO until sqrt has a fallback
+    #[cfg(feature = "std")]
+    mod ray_sphere {
+        use super::*;
+
+        const SPHERE: Sphere = Sphere(pt3(0.0, 0.0, 1.0), 2.0);
+
+        #[test]
+        fn ray_passes_through_sphere() {
+            let ray: Ray3 = Ray(pt3(0.0, 0.0, -3.0), vec3(0.0, 0.0, 2.0));
+            assert_eq!(
+                ray.intersect(&SPHERE),
+                Some((1.0, pt3(0.0, 0.0, -1.0)))
+            );
+        }
+        #[test]
+        fn ray_tangent_to_sphere() {
+            let ray: Ray3 = Ray(pt3(0.0, 2.0, -3.0), vec3(0.0, 0.0, 2.0));
+            assert_eq!(ray.intersect(&SPHERE), Some((2.0, pt3(0.0, 2.0, 1.0))));
+        }
+
+        #[test]
+        fn ray_origin_inside_sphere() {
+            let ray: Ray3 = Ray(pt3(0.0, 0.0, 0.0), vec3(0.0, 0.0, 2.0));
+            assert_eq!(ray.intersect(&SPHERE), Some((1.5, pt3(0.0, 0.0, 3.0))));
+
+            let ray: Ray3 = Ray(pt3(0.0, 0.0, 2.0), vec3(0.0, 0.0, 2.0));
+            assert_eq!(ray.intersect(&SPHERE), Some((0.5, pt3(0.0, 0.0, 3.0))));
+        }
+
+        fn sphere_behind_ray() {
+            let ray: Ray3 = Ray(pt3(0.0, 0.0, -3.0), vec3(0.0, 0.0, -1.0));
+            assert_eq!(ray.intersect(&SPHERE), None);
+        }
+
+        #[test]
+        fn ray_misses_sphere() {
+            let ray: Ray3 = Ray(pt3(0.0, 0.0, -3.0), vec3(0.0, 1.0, 1.0));
+            assert_eq!(ray.intersect(&SPHERE), None);
+        }
+    }
 }