use bone_types::{
    Aabb3, AngleTolerance, ChordHeightTolerance, Parameter, Point3, Tolerance, UnitVec3, Vec3,
};

use crate::KernelError;
use crate::mesh::{MeshVertex, TriMesh};
use crate::nurbs::{
    Degree, KnotVector, RationalControlPoint, evaluate_homogeneous, project, project_derivative,
};
use crate::surface3::Surface3;

const MIN_SEGMENTS: u32 = 4;
const MAX_SEGMENTS: u32 = 256;
const SEED_SAMPLES: u32 = 8;
const PROJECTION_STEPS: u32 = 24;
const PROBE_LINES: [f64; 5] = [0.0, 0.25, 0.5, 0.75, 1.0];
const NORMAL_TOLERANCE: Tolerance = Tolerance::new(1.0e-12);

#[derive(Clone, Debug, PartialEq)]
pub struct NurbsSurface3 {
    degree_u: Degree,
    degree_v: Degree,
    knots_u: KnotVector,
    knots_v: KnotVector,
    rows: Vec<Vec<RationalControlPoint>>,
    homogeneous: Vec<Vec<[f64; 4]>>,
}

impl NurbsSurface3 {
    pub fn new(
        degree_u: Degree,
        degree_v: Degree,
        knots_u: KnotVector,
        knots_v: KnotVector,
        rows: Vec<Vec<RationalControlPoint>>,
    ) -> Result<Self, KernelError> {
        let columns = rows.first().map_or(0, Vec::len);
        let rectangular = rows.iter().all(|row| row.len() == columns);
        let u_enough = rows.len() > degree_u.get();
        let v_enough = columns > degree_v.get();
        let u_knots_ok = knots_u.values().len() == rows.len() + degree_u.get() + 1
            && knots_u.is_clamped(degree_u);
        let v_knots_ok =
            knots_v.values().len() == columns + degree_v.get() + 1 && knots_v.is_clamped(degree_v);
        if !rectangular || !u_enough || !v_enough || !u_knots_ok || !v_knots_ok {
            return Err(KernelError::DegenerateNurbs);
        }
        let homogeneous = rows
            .iter()
            .map(|row| row.iter().map(|point| point.homogeneous()).collect())
            .collect();
        Ok(Self {
            degree_u,
            degree_v,
            knots_u,
            knots_v,
            rows,
            homogeneous,
        })
    }

    #[must_use]
    pub const fn degree_u(&self) -> Degree {
        self.degree_u
    }

    #[must_use]
    pub const fn degree_v(&self) -> Degree {
        self.degree_v
    }

    #[must_use]
    fn param_u(&self, t: f64) -> f64 {
        let (lo, hi) = self.knots_u.domain(self.degree_u);
        lo + (hi - lo) * t
    }

    #[must_use]
    fn param_v(&self, t: f64) -> f64 {
        let (lo, hi) = self.knots_v.domain(self.degree_v);
        lo + (hi - lo) * t
    }

    #[must_use]
    fn span_u(&self) -> f64 {
        let (lo, hi) = self.knots_u.domain(self.degree_u);
        hi - lo
    }

    #[must_use]
    fn span_v(&self) -> f64 {
        let (lo, hi) = self.knots_v.domain(self.degree_v);
        hi - lo
    }

    #[must_use]
    fn homogeneous_at(&self, u: f64, v: f64) -> ([f64; 4], [f64; 4], [f64; 4]) {
        let pu = self.param_u(u);
        let pv = self.param_v(v);
        let (row_value, row_dv): (Vec<[f64; 4]>, Vec<[f64; 4]>) = self
            .homogeneous
            .iter()
            .map(|row| evaluate_homogeneous(self.degree_v.get(), self.knots_v.values(), row, pv))
            .unzip();
        let (value, du) =
            evaluate_homogeneous(self.degree_u.get(), self.knots_u.values(), &row_value, pu);
        let (dv, _) = evaluate_homogeneous(self.degree_u.get(), self.knots_u.values(), &row_dv, pu);
        (value, du, dv)
    }

    #[must_use]
    fn point_at(&self, u: f64, v: f64) -> Point3 {
        project(self.homogeneous_at(u, v).0)
    }

    #[must_use]
    fn normal_at(&self, u: f64, v: f64) -> UnitVec3 {
        let (du, dv) = self.partials_at(u, v);
        du.cross(dv)
            .try_normalize(NORMAL_TOLERANCE)
            .unwrap_or_else(|_| UnitVec3::z_axis())
    }

    #[must_use]
    fn partials_at(&self, u: f64, v: f64) -> (Vec3, Vec3) {
        let (value, du, dv) = self.homogeneous_at(u, v);
        (
            project_derivative(value, du) * self.span_u(),
            project_derivative(value, dv) * self.span_v(),
        )
    }

    #[must_use]
    fn sample_along(&self, t: f64, fixed: f64, along_u: bool) -> (Point3, Vec3) {
        let (u, v) = if along_u { (t, fixed) } else { (fixed, t) };
        let (du, dv) = self.partials_at(u, v);
        (self.point_at(u, v), if along_u { du } else { dv })
    }

    #[must_use]
    fn segment_within(
        &self,
        t0: f64,
        t1: f64,
        fixed: f64,
        along_u: bool,
        chord: f64,
        angle: f64,
    ) -> bool {
        let (p0, d0) = self.sample_along(t0, fixed, along_u);
        let (p1, d1) = self.sample_along(t1, fixed, along_u);
        let (pm, _) = self.sample_along(f64::midpoint(t0, t1), fixed, along_u);
        let chord_ok = chord <= 0.0 || (pm - midpoint(p0, p1)).norm_mm() <= chord;
        let angle_ok = angle <= 0.0 || turn(d0, d1) <= angle;
        chord_ok && angle_ok
    }

    #[must_use]
    fn direction_resolved(&self, n: u32, along_u: bool, chord: f64, angle: f64) -> bool {
        let inv = 1.0 / f64::from(n);
        PROBE_LINES.iter().all(|&fixed| {
            (0..n).all(|i| {
                self.segment_within(
                    f64::from(i) * inv,
                    f64::from(i + 1) * inv,
                    fixed,
                    along_u,
                    chord,
                    angle,
                )
            })
        })
    }

    #[must_use]
    fn direction_samples(
        &self,
        chord: ChordHeightTolerance,
        angle: AngleTolerance,
        along_u: bool,
    ) -> u32 {
        let chord_mm = chord.millimeters();
        let angle_rad = angle.radians();
        core::iter::successors(Some(MIN_SEGMENTS), |&n| {
            (n < MAX_SEGMENTS).then(|| (2 * n).min(MAX_SEGMENTS))
        })
        .find(|&n| self.direction_resolved(n, along_u, chord_mm, angle_rad))
        .unwrap_or(MAX_SEGMENTS)
    }

    #[must_use]
    fn seed_uv(&self, p: Point3) -> (f64, f64) {
        let inv = 1.0 / f64::from(SEED_SAMPLES);
        (0..=SEED_SAMPLES)
            .flat_map(|i| (0..=SEED_SAMPLES).map(move |j| (i, j)))
            .map(|(i, j)| {
                let (u, v) = (f64::from(i) * inv, f64::from(j) * inv);
                (u, v, (self.point_at(u, v) - p).norm_mm())
            })
            .min_by(|a, b| a.2.total_cmp(&b.2))
            .map_or((0.0, 0.0), |(u, v, _)| (u, v))
    }

    #[must_use]
    fn project(&self, p: Point3, seed: (f64, f64), tolerance: Tolerance) -> (f64, f64) {
        let stepped = (0..PROJECTION_STEPS).try_fold(seed, |(u, v), _| {
            let point = self.point_at(u, v);
            let (su, sv) = self.partials_at(u, v);
            let r = point - p;
            let a11 = su.dot_mm2(su);
            let a12 = su.dot_mm2(sv);
            let a22 = sv.dot_mm2(sv);
            let det = a11 * a22 - a12 * a12;
            if det.abs() <= f64::EPSILON {
                return Err((u, v));
            }
            let b1 = r.dot_mm2(su);
            let b2 = r.dot_mm2(sv);
            let next = (
                (u + (b2 * a12 - b1 * a22) / det).clamp(0.0, 1.0),
                (v + (b1 * a12 - b2 * a11) / det).clamp(0.0, 1.0),
            );
            if (self.point_at(next.0, next.1) - point).norm_mm() <= tolerance.value() {
                Err(next)
            } else {
                Ok(next)
            }
        });
        match stepped {
            Ok(uv) | Err(uv) => uv,
        }
    }
}

fn midpoint(a: Point3, b: Point3) -> Point3 {
    let (ax, ay, az) = a.coords_mm();
    let (bx, by, bz) = b.coords_mm();
    Point3::from_mm(
        f64::midpoint(ax, bx),
        f64::midpoint(ay, by),
        f64::midpoint(az, bz),
    )
}

fn turn(a: Vec3, b: Vec3) -> f64 {
    a.cross(b).norm_mm().atan2(a.dot_mm2(b))
}

impl Surface3 for NurbsSurface3 {
    fn evaluate(&self, u: Parameter, v: Parameter) -> Point3 {
        self.point_at(u.value(), v.value())
    }

    fn partials(&self, u: Parameter, v: Parameter) -> (Vec3, Vec3) {
        self.partials_at(u.value(), v.value())
    }

    fn normal(&self, u: Parameter, v: Parameter) -> UnitVec3 {
        self.normal_at(u.value(), v.value())
    }

    fn bounding_box(&self) -> Aabb3 {
        Aabb3::from_points(
            self.rows
                .iter()
                .flat_map(|row| row.iter().map(RationalControlPoint::position)),
        )
        .unwrap_or_else(|| Aabb3::from_corners(Point3::origin(), Point3::origin()))
    }

    fn contains_point(&self, p: Point3, tolerance: Tolerance) -> bool {
        let (u, v) = self.project(p, self.seed_uv(p), tolerance);
        (self.point_at(u, v) - p).norm_mm() <= tolerance.value()
    }

    fn tessellate(&self, chord: ChordHeightTolerance, angle: AngleTolerance) -> TriMesh {
        let cols = self.direction_samples(chord, angle, true);
        let rows = self.direction_samples(chord, angle, false);
        let inv_u = 1.0 / f64::from(cols);
        let inv_v = 1.0 / f64::from(rows);
        TriMesh::from_grid(cols + 1, rows + 1, |i, j| {
            let u = f64::from(i) * inv_u;
            let v = f64::from(j) * inv_v;
            MeshVertex::new(self.point_at(u, v), self.normal_at(u, v))
        })
    }
}

impl core::fmt::Display for NurbsSurface3 {
    fn fmt(&self, f: &mut core::fmt::Formatter<'_>) -> core::fmt::Result {
        write!(
            f,
            "nurbs_surface3{{ degree=({}, {}), control_net={}x{} }}",
            self.degree_u.get(),
            self.degree_v.get(),
            self.rows.len(),
            self.rows.first().map_or(0, Vec::len),
        )
    }
}
