use bone_types::{AngleTolerance, ChordHeightTolerance, Point3, Tolerance, UnitVec3, Vec3};
use robust::{Coord3D, orient3d};

use super::super::edges::EdgeCurve3;
use super::super::{BrepError, BrepSolid};
use super::as_vec;
use super::model::{Body, Face, PlanarFace, Tube, dot};
use crate::curve3::Curve3;

const RAY_SPAN_MM: f64 = 1.0e6;
const SOUP_CHORD_MM: f64 = 0.02;
const SOUP_ANGLE_RAD: f64 = 0.08;

fn ray_directions() -> [Vec3; 4] {
    [
        Vec3::from_mm(0.732_191, 0.451_77, 0.509_63),
        Vec3::from_mm(-0.611_3, 0.701_1, 0.367_9),
        Vec3::from_mm(0.289_4, -0.534_2, 0.795_1),
        Vec3::from_mm(0.456_8, 0.241_3, -0.856_2),
    ]
}

#[derive(Copy, Clone, Debug, PartialEq, Eq)]
pub enum FaceClass {
    In,
    Out,
    OnSame,
    OnOpposite,
}

pub(crate) struct SolidSoup {
    tris: Vec<[Point3; 3]>,
}

impl SolidSoup {
    pub(crate) fn new(solid: &BrepSolid, _tolerance: Tolerance) -> Result<Self, BrepError> {
        let mesh = solid
            .tessellate(
                ChordHeightTolerance::from_mm(SOUP_CHORD_MM),
                AngleTolerance::from_radians(SOUP_ANGLE_RAD),
            )
            .map_err(|_| BrepError::BooleanIndeterminate)?;
        let tris = mesh
            .faces()
            .iter()
            .flat_map(|slab| {
                let positions = slab.positions().to_vec();
                slab.triangles().iter().map(move |triangle| {
                    [
                        positions[triangle[0] as usize],
                        positions[triangle[1] as usize],
                        positions[triangle[2] as usize],
                    ]
                })
            })
            .collect();
        Ok(Self { tris })
    }

    pub(crate) fn inside(&self, point: Point3) -> Option<bool> {
        ray_directions()
            .into_iter()
            .find_map(|dir| self.cast(point, dir))
    }

    fn cast(&self, point: Point3, dir: Vec3) -> Option<bool> {
        let far = point + dir * RAY_SPAN_MM;
        let crossings = self.tris.iter().try_fold(0u32, |count, tri| {
            match segment_triangle(point, far, tri[0], tri[1], tri[2]) {
                Crossing::Ambiguous => Err(()),
                Crossing::Hit => Ok(count + 1),
                Crossing::Miss => Ok(count),
            }
        });
        crossings.ok().map(|count| count % 2 == 1)
    }
}

fn coord(p: Point3) -> Coord3D<f64> {
    let (x, y, z) = p.coords_mm();
    Coord3D { x, y, z }
}

enum Crossing {
    Hit,
    Miss,
    Ambiguous,
}

fn segment_triangle(c: Point3, q: Point3, p0: Point3, p1: Point3, p2: Point3) -> Crossing {
    let (cc, cq) = (coord(c), coord(q));
    let (a0, a1, a2) = (coord(p0), coord(p1), coord(p2));
    let side_c = orient3d(a0, a1, a2, cc);
    let side_q = orient3d(a0, a1, a2, cq);
    if side_c == 0.0 {
        return Crossing::Miss;
    }
    if side_q == 0.0 {
        return Crossing::Ambiguous;
    }
    if (side_c > 0.0) == (side_q > 0.0) {
        return Crossing::Miss;
    }
    let s0 = orient3d(cc, cq, a0, a1);
    let s1 = orient3d(cc, cq, a1, a2);
    let s2 = orient3d(cc, cq, a2, a0);
    if s0 == 0.0 || s1 == 0.0 || s2 == 0.0 {
        return Crossing::Ambiguous;
    }
    if (s0 > 0.0) == (s1 > 0.0) && (s1 > 0.0) == (s2 > 0.0) {
        Crossing::Hit
    } else {
        Crossing::Miss
    }
}

fn cross2(o: (f64, f64), p: (f64, f64), q: (f64, f64)) -> f64 {
    (p.0 - o.0) * (q.1 - o.1) - (p.1 - o.1) * (q.0 - o.0)
}

fn point_in_triangle(a: (f64, f64), b: (f64, f64), c: (f64, f64), p: (f64, f64)) -> bool {
    let d0 = cross2(a, b, p);
    let d1 = cross2(b, c, p);
    let d2 = cross2(c, a, p);
    let any_neg = d0 < 0.0 || d1 < 0.0 || d2 < 0.0;
    let any_pos = d0 > 0.0 || d1 > 0.0 || d2 > 0.0;
    !(any_neg && any_pos)
}

fn distance_to_line(a: (f64, f64), b: (f64, f64), p: (f64, f64)) -> f64 {
    let length = (b.0 - a.0).hypot(b.1 - a.1);
    if length <= 0.0 {
        (p.0 - a.0).hypot(p.1 - a.1)
    } else {
        cross2(a, b, p).abs() / length
    }
}

pub(super) fn local(face: &PlanarFace, p: Point3) -> (f64, f64) {
    let offset = p - face.plane.origin();
    (
        offset.dot_mm2(as_vec(face.plane.x_axis())),
        offset.dot_mm2(as_vec(face.plane.y_axis())),
    )
}

fn polygon_locals(face: &PlanarFace) -> Vec<(f64, f64)> {
    face.polygon
        .iter()
        .map(|v| local(face, v.position))
        .collect()
}

fn interior_point(face: &PlanarFace) -> Point3 {
    let locals = polygon_locals(face);
    let count = locals.len();
    let lowest = (0..count)
        .min_by(|&i, &j| {
            (locals[i].1, locals[i].0)
                .partial_cmp(&(locals[j].1, locals[j].0))
                .unwrap_or(core::cmp::Ordering::Equal)
        })
        .unwrap_or(0);
    let before = (lowest + count - 1) % count;
    let after = (lowest + 1) % count;
    let (a, ear, c) = (locals[before], locals[lowest], locals[after]);
    let contained: Vec<(f64, f64)> = (0..count)
        .filter(|&i| i != before && i != lowest && i != after)
        .map(|i| locals[i])
        .filter(|&p| point_in_triangle(a, ear, c, p))
        .collect();
    let inside = contained
        .into_iter()
        .max_by(|p, q| {
            distance_to_line(a, c, *p)
                .partial_cmp(&distance_to_line(a, c, *q))
                .unwrap_or(core::cmp::Ordering::Equal)
        })
        .map_or(
            ((a.0 + ear.0 + c.0) / 3.0, (a.1 + ear.1 + c.1) / 3.0),
            |deepest| {
                (
                    f64::midpoint(ear.0, deepest.0),
                    f64::midpoint(ear.1, deepest.1),
                )
            },
        );
    face.plane.point_at_local(inside.0, inside.1, 0.0)
}

pub(super) fn on_plane(face: &PlanarFace, p: Point3, tolerance: Tolerance) -> bool {
    (p - face.plane.origin())
        .dot_mm2(as_vec(face.plane.normal()))
        .abs()
        <= tolerance.value()
}

fn on_polygon_boundary(locals: &[(f64, f64)], p: (f64, f64), tolerance: f64) -> bool {
    let count = locals.len();
    (0..count).any(|index| near_segment(p, locals[index], locals[(index + 1) % count], tolerance))
}

fn winding_inside(locals: &[(f64, f64)], p: (f64, f64)) -> bool {
    let (px, py) = p;
    let count = locals.len();
    (0..count).fold(false, |inside, index| {
        let (xi, yi) = locals[index];
        let (xj, yj) = locals[(index + count - 1) % count];
        let crosses = (yi > py) != (yj > py) && px < (xj - xi) * (py - yi) / (yj - yi) + xi;
        inside ^ crosses
    })
}

fn near_segment(p: (f64, f64), a: (f64, f64), b: (f64, f64), tolerance: f64) -> bool {
    let (dx, dy) = (b.0 - a.0, b.1 - a.1);
    let length_sq = dx * dx + dy * dy;
    let t = if length_sq <= 0.0 {
        0.0
    } else {
        (((p.0 - a.0) * dx + (p.1 - a.1) * dy) / length_sq).clamp(0.0, 1.0)
    };
    let foot = (a.0 + t * dx, a.1 + t * dy);
    (p.0 - foot.0).hypot(p.1 - foot.1) <= tolerance
}

fn tessellate_loop(verts: &[super::model::PlanarVertex], curves: &[EdgeCurve3]) -> Vec<Point3> {
    if curves.is_empty() {
        return verts.iter().map(|v| v.position).collect();
    }
    curves
        .iter()
        .flat_map(|curve| curve.tessellate(ChordHeightTolerance::from_mm(0.05)))
        .collect()
}

fn is_line(curve: &EdgeCurve3) -> bool {
    matches!(curve, EdgeCurve3::Line(_))
}

fn is_straight_polygon(face: &PlanarFace) -> bool {
    face.holes.is_empty() && face.polygon.len() >= 3 && face.curves.iter().all(is_line)
}

fn average(points: &[Point3]) -> Option<Point3> {
    if points.is_empty() {
        return None;
    }
    let (count, sx, sy, sz) = points
        .iter()
        .fold((0.0, 0.0, 0.0, 0.0), |(n, sx, sy, sz), p| {
            let (x, y, z) = p.coords_mm();
            (n + 1.0, sx + x, sy + y, sz + z)
        });
    Some(Point3::from_mm(sx / count, sy / count, sz / count))
}

fn in_circle_hole(hole: &super::model::PlanarLoop, p: (f64, f64), face: &PlanarFace) -> bool {
    if let Some(EdgeCurve3::Circle(circle)) = hole.curves.first() {
        let c = local(face, circle.center());
        (p.0 - c.0).hypot(p.1 - c.1) < circle.radius_mm()
    } else {
        let locals: Vec<(f64, f64)> = hole.verts.iter().map(|v| local(face, v.position)).collect();
        winding_inside(&locals, p)
    }
}

fn in_any_hole(face: &PlanarFace, p: (f64, f64)) -> bool {
    face.holes.iter().any(|hole| in_circle_hole(hole, p, face))
}

fn planar_center(face: &PlanarFace) -> Option<Point3> {
    const FRACTIONS: [f64; 5] = [0.5, 0.7, 0.85, 0.95, 0.99];
    if is_straight_polygon(face) {
        return Some(interior_point(face));
    }
    let outer = tessellate_loop(&face.polygon, &face.curves);
    let centroid = average(&outer)?;
    if face.holes.is_empty() {
        return Some(centroid);
    }
    let centroid_local = local(face, centroid);
    if !in_any_hole(face, centroid_local) {
        return Some(centroid);
    }
    outer
        .iter()
        .flat_map(|rim| {
            FRACTIONS.iter().map(move |&t| {
                let (cx, cy, cz) = centroid.coords_mm();
                let (rx, ry, rz) = rim.coords_mm();
                Point3::from_mm(cx + (rx - cx) * t, cy + (ry - cy) * t, cz + (rz - cz) * t)
            })
        })
        .find(|cand| !in_any_hole(face, local(face, *cand)))
        .or(Some(centroid))
}

fn face_center(face: &Face) -> Option<Point3> {
    match face {
        Face::Tube(tube) => Some(tube.center()),
        Face::Planar(planar) => planar_center(planar),
    }
}

fn face_outward_normal(face: &Face, tolerance: Tolerance) -> Option<UnitVec3> {
    match face {
        Face::Planar(planar) => Some(planar.plane.normal()),
        Face::Tube(tube) => tube.surface_normal(tolerance),
    }
}

#[derive(Copy, Clone, PartialEq, Eq)]
enum Containment {
    Outside,
    Boundary,
    Interior,
}

pub(super) fn region_contains(face: &PlanarFace, p: Point3, tolerance: Tolerance) -> bool {
    !matches!(containment(face, p, tolerance), Containment::Outside)
}

pub(super) fn region_interior(face: &PlanarFace, p: Point3, tolerance: Tolerance) -> bool {
    matches!(containment(face, p, tolerance), Containment::Interior)
}

fn containment(face: &PlanarFace, p: Point3, tolerance: Tolerance) -> Containment {
    if !on_plane(face, p, tolerance) || in_any_hole(face, local(face, p)) {
        return Containment::Outside;
    }
    let l = local(face, p);
    let eps = tolerance.value();
    if let Some(EdgeCurve3::Circle(circle)) = face.curves.first()
        && face.polygon.len() == 1
    {
        let c = local(face, circle.center());
        let distance = (l.0 - c.0).hypot(l.1 - c.1);
        let radius = circle.radius_mm();
        return if (distance - radius).abs() <= eps {
            Containment::Boundary
        } else if distance < radius {
            Containment::Interior
        } else {
            Containment::Outside
        };
    }
    let locals = if is_straight_polygon(face) {
        polygon_locals(face)
    } else {
        tessellate_loop(&face.polygon, &face.curves)
            .iter()
            .map(|q| local(face, *q))
            .collect()
    };
    if on_polygon_boundary(&locals, l, eps) {
        Containment::Boundary
    } else if winding_inside(&locals, l) {
        Containment::Interior
    } else {
        Containment::Outside
    }
}

fn tube_contains(tube: &Tube, p: Point3, tolerance: Tolerance) -> bool {
    let radial = super::model::radial_distance(tube.cyl.axis_origin, tube.cyl.axis, p);
    if (radial - tube.cyl.radius_mm()).abs() > tolerance.value() {
        return false;
    }
    let z = super::model::project_onto_axis(tube.cyl.axis_origin, tube.cyl.axis, p);
    if z < tube.z0 - tolerance.value() || z > tube.z1 + tolerance.value() {
        return false;
    }
    let (theta, _) = tube.cyl.theta_z(p);
    let span = tube.theta1 - tube.theta0;
    if (span - core::f64::consts::TAU).abs() <= tolerance.value() {
        return true;
    }
    let rel = (theta - tube.theta0).rem_euclid(core::f64::consts::TAU);
    rel <= span + tolerance.value()
}

fn supporting_normal(other: &Body, p: Point3, tolerance: Tolerance) -> Option<UnitVec3> {
    other.iter().find_map(|face| match face {
        Face::Planar(planar) => {
            region_contains(planar, p, tolerance).then(|| planar.plane.normal())
        }
        Face::Tube(tube) => {
            if tube_contains(tube, p, tolerance) {
                let (theta, _) = tube.cyl.theta_z(p);
                let radial = tube.cyl.radial(theta, tolerance)?;
                Some(if tube.outward {
                    radial
                } else {
                    radial.reversed()
                })
            } else {
                None
            }
        }
    })
}

pub(crate) fn classify_face(
    face: &Face,
    other: &Body,
    soup: &SolidSoup,
    tolerance: Tolerance,
) -> Option<FaceClass> {
    let center = face_center(face)?;
    if let Some(support_normal) = supporting_normal(other, center, tolerance) {
        let this_normal = face_outward_normal(face, tolerance)?;
        let aligned = dot(this_normal, support_normal);
        return Some(if aligned >= 0.0 {
            FaceClass::OnSame
        } else {
            FaceClass::OnOpposite
        });
    }
    match soup.inside(center) {
        Some(true) => Some(FaceClass::In),
        Some(false) => Some(FaceClass::Out),
        None => None,
    }
}
