//! Ellipsoid geometry, shared by the synthetic scene (replay), the
//! occlusion proxy (splatgen), and tests.

use glam::Vec3;

/// An axis-aligned ellipsoid: all points `p` with
/// `((p - center) / radii)·length() = 1`.
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct Ellipsoid {
    pub center: Vec3,
    /// Per-axis radii, metres. All positive.
    pub radii: Vec3,
}

impl Ellipsoid {
    pub fn new(center: Vec3, radii: Vec3) -> Self {
        Self { center, radii }
    }

    /// Does the ellipsoid contain the point (strictly inside)?
    pub fn contains(&self, p: Vec3) -> bool {
        let s = (p - self.center) / self.radii;
        s.dot(s) < 1.0
    }

    /// The outward unit normal at a surface point.
    pub fn normal_at(&self, p: Vec3) -> Vec3 {
        let grad = (p - self.center) / (self.radii * self.radii);
        grad.normalize_or_zero()
    }

    /// The first intersection of the ray `origin + t * dir` with the
    /// ellipsoid, at `t > 0`, with the surface normal at the hit.
    ///
    /// `dir` need not be normalized; the returned `t` is in the ray's
    /// own parameter space. Rays starting inside exit through the far
    /// side.
    pub fn intersect_ray(&self, origin: Vec3, dir: Vec3) -> Option<(f32, Vec3)> {
        // Scale space so the ellipsoid is the unit sphere.
        let o = (origin - self.center) / self.radii;
        let d = dir / self.radii;
        let a = d.dot(d);
        if a <= f32::EPSILON {
            return None;
        }
        let b = 2.0 * o.dot(d);
        let c = o.dot(o) - 1.0;
        let disc = b * b - 4.0 * a * c;
        if disc < 0.0 {
            return None;
        }
        let sqrt_disc = disc.sqrt();
        let mut t = (-b - sqrt_disc) / (2.0 * a);
        if t <= f32::EPSILON {
            t = (-b + sqrt_disc) / (2.0 * a);
        }
        if t <= f32::EPSILON {
            return None;
        }
        let position = origin + dir * t;
        Some((t, self.normal_at(position)))
    }

    /// An approximate bounding check: does the ray ever come within
    /// `margin` of the ellipsoid's surface? Used to skip work cheaply.
    #[allow(dead_code)]
    pub fn near_ray(&self, origin: Vec3, dir: Vec3, margin: f32) -> bool {
        let o = (origin - self.center) / self.radii;
        let d = dir / self.radii;
        let a = d.dot(d);
        if a <= f32::EPSILON {
            return o.dot(o) < (1.0 + margin).powi(2);
        }
        let t = -o.dot(d) / a;
        let closest = if t > 0.0 { o + d * t } else { o };
        closest.dot(closest) < (1.0 + margin).powi(2)
    }
}

#[cfg(test)]
mod tests {
    use super::*;
    use glam::vec3;

    fn unit() -> Ellipsoid {
        Ellipsoid::new(vec3(0.0, 0.0, 2.0), vec3(1.0, 1.0, 1.0))
    }

    #[test]
    fn membership() {
        let e = unit();
        assert!(e.contains(vec3(0.0, 0.0, 2.0)));
        assert!(e.contains(vec3(0.5, 0.0, 2.5)));
        assert!(!e.contains(vec3(1.0, 0.0, 2.0)));
        assert!(!e.contains(vec3(0.0, 0.0, 0.0)));
    }

    #[test]
    fn ray_hits_sphere() {
        let e = unit();
        // Straight on: hits the front pole, normal faces the camera.
        let (t, n) = e
            .intersect_ray(Vec3::ZERO, vec3(0.0, 0.0, 1.0))
            .expect("hit");
        assert!((t - 1.0).abs() < 1e-5);
        assert!((n - vec3(0.0, 0.0, -1.0)).length() < 1e-5);

        // Miss.
        assert!(
            e.intersect_ray(vec3(5.0, 0.0, 0.0), vec3(0.0, 0.0, 1.0))
                .is_none()
        );
        // Behind the ray.
        assert!(
            e.intersect_ray(vec3(0.0, 0.0, 5.0), vec3(0.0, 0.0, 1.0))
                .is_none()
        );
        // From inside: exits at the far side.
        let (t, _) = e
            .intersect_ray(vec3(0.0, 0.0, 2.0), vec3(0.0, 0.0, 1.0))
            .expect("exits");
        assert!((t - 1.0).abs() < 1e-5);
    }

    #[test]
    fn ray_hits_scaled_ellipsoid() {
        let e = Ellipsoid::new(vec3(0.0, 0.18, 2.0), vec3(0.09, 0.11, 0.10));
        // Straight at the front: t is in unnormalized ray units.
        let (t, n) = e
            .intersect_ray(Vec3::ZERO, vec3(0.0, 0.18, 2.0))
            .expect("hit");
        let p = Vec3::ZERO + vec3(0.0, 0.18, 2.0) * t;
        // The hit is on the surface and the normal points back at us.
        let s = (p - e.center) / e.radii;
        assert!((s.length() - 1.0).abs() < 1e-4);
        assert!(n.z < 0.0);

        // Side rays hit with tilted normals.
        let (t, n) = e
            .intersect_ray(vec3(-1.0, 0.18, 2.0), vec3(1.0, 0.0, 0.0))
            .expect("side hit");
        let p = vec3(-1.0, 0.18, 2.0) + vec3(1.0, 0.0, 0.0) * t;
        assert!((p.x - (-0.09)).abs() < 1e-4, "x = {}", p.x);
        assert!(n.x < 0.0 && n.z.abs() < 0.01);
    }
}