modify nalgebra sphere-ray intersection

del-geo-core/src/vec3.rs

@@ -255,6 +255,15 @@ where
     ]
 }
 
+pub fn mirror_reflection(v: &[f32;3], nrm: &[f32;3]) -> [f32;3] {
+    let a = dot(nrm, v);
+    [
+        v[0] - nrm[0] * 2. * a,
+        v[1] - nrm[1] * 2. * a,
+        v[2] - nrm[2] * 2. * a,
+    ]
+}
+
 // ------------------------------------------
 pub struct XYZ<'a, Real> {
     pub p: &'a [Real; 3],

del-geo-nalgebra/src/sphere.rs

@@ -1,22 +1,28 @@
-pub fn intersection_ray(
-    center: &nalgebra::Vector3<f32>,
-    rad: f32,
-    ray_src: &nalgebra::Vector3<f32>,
-    ray_dir: &nalgebra::Vector3<f32>,
-) -> Option<(nalgebra::Vector3<f32>, nalgebra::Vector3<f32>, f32)> {
-    let depth0 = (center - ray_src).dot(ray_dir);
-    if depth0 < 0f32 {
-        return None;
+
+/// ray_dir is not always a unit vector
+pub fn intersection_ray<T>(
+    center: &nalgebra::Vector3<T>,
+    rad: T,
+    ray_src: &nalgebra::Vector3<T>,
+    ray_dir: &nalgebra::Vector3<T>,
+) -> Option<T>
+where T: nalgebra::RealField + Copy
+{
+    // Solve t^2*d.d + 2*t*(o-p).d + (o-p).(o-p)-R^2 = 0
+    let op = ray_src - center;
+    let a = ray_dir.dot(&ray_dir);
+    let b = op.dot(&ray_dir);
+    let c = op.dot(&op) - rad * rad;
+    let det: T = b * b - c * a;
+    if det < T::zero() {
+        None
+    } else {
+        let det = det.sqrt();
+        if -b - det >= T::zero() {
+            return Some((-b - det)/a);
+        } else if -b + det >= T::zero() {
+            return Some((-b + det)/a);
+        }
+        None
     }
-    let sqdist = (ray_src + depth0 * ray_dir - center).norm_squared();
-    if rad * rad - sqdist < 0f32 {
-        return None;
-    }
-    let depth1 = depth0 - (rad * rad - sqdist).sqrt();
-    if depth1 < 0f32 {
-        return None;
-    }
-    let hit_pos = ray_src + depth1 * ray_dir;
-    let hit_normal = (hit_pos - center).normalize();
-    Some((hit_pos, hit_normal, depth1))
 }