// Copyright 2018 the Kurbo Authors

// SPDX-License-Identifier: Apache-2.0 OR MIT

//! Quadratic Bézier segments.

use core::ops::{Mul, Range};

use arrayvec::ArrayVec;

use crate::common::solve_cubic;

use crate::MAX_EXTREMA;

use crate::{

    Affine, CubicBez, Line, Nearest, ParamCurve, ParamCurveArclen, ParamCurveArea,

    ParamCurveCurvature, ParamCurveDeriv, ParamCurveExtrema, ParamCurveNearest, PathEl, Point,

    Rect, Shape,

};

#[cfg(not(feature = "std"))]

use crate::common::FloatFuncs;

/// A single quadratic Bézier segment.

#[derive(Clone, Copy, Debug, PartialEq)]

#[cfg_attr(feature = "schemars", derive(schemars::JsonSchema))]

#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]

#[allow(missing_docs)]

pub struct QuadBez {

    pub p0: Point,

    pub p1: Point,

    pub p2: Point,

}

impl QuadBez {

    /// Create a new quadratic Bézier segment.

    #[inline]

    pub fn new<V: Into<Point>>(p0: V, p1: V, p2: V) -> QuadBez {

        QuadBez {

            p0: p0.into(),

            p1: p1.into(),

            p2: p2.into(),

        }

    }

    /// Raise the order by 1.

    ///

    /// Returns a cubic Bézier segment that exactly represents this quadratic.

    #[inline]

    pub fn raise(&self) -> CubicBez {

        CubicBez::new(

            self.p0,

            self.p0 + (2.0 / 3.0) * (self.p1 - self.p0),

            self.p2 + (2.0 / 3.0) * (self.p1 - self.p2),

            self.p2,

        )

    }

    /// Estimate the number of subdivisions for flattening.

    pub(crate) fn estimate_subdiv(&self, sqrt_tol: f64) -> FlattenParams {

        // Determine transformation to $y = x^2$ parabola.

        let d01 = self.p1 - self.p0;

        let d12 = self.p2 - self.p1;

        let dd = d01 - d12;

        let cross = (self.p2 - self.p0).cross(dd);

        let x0 = d01.dot(dd) * cross.recip();

        let x2 = d12.dot(dd) * cross.recip();

        let scale = (cross / (dd.hypot() * (x2 - x0))).abs();

        // Compute number of subdivisions needed.

        let a0 = approx_parabola_integral(x0);

        let a2 = approx_parabola_integral(x2);

        let val = if scale.is_finite() {

            let da = (a2 - a0).abs();

            let sqrt_scale = scale.sqrt();

            if x0.signum() == x2.signum() {

                da * sqrt_scale

            } else {

                // Handle cusp case (segment contains curvature maximum)

                let xmin = sqrt_tol / sqrt_scale;

                sqrt_tol * da / approx_parabola_integral(xmin)

            }

        } else {

            0.0

        };

        let u0 = approx_parabola_inv_integral(a0);

        let u2 = approx_parabola_inv_integral(a2);

        let uscale = (u2 - u0).recip();

        FlattenParams {

            a0,

            a2,

            u0,

            uscale,

            val,

        }

    }

    // Maps a value from 0..1 to 0..1.

    pub(crate) fn determine_subdiv_t(&self, params: &FlattenParams, x: f64) -> f64 {

        let a = params.a0 + (params.a2 - params.a0) * x;

        let u = approx_parabola_inv_integral(a);

        (u - params.u0) * params.uscale

    }

    /// Is this quadratic Bezier curve finite?

    #[inline]

    pub const fn is_finite(&self) -> bool {

        self.p0.is_finite() && self.p1.is_finite() && self.p2.is_finite()

    }

    /// Is this quadratic Bezier curve NaN?

    #[inline]

    pub const fn is_nan(&self) -> bool {

        self.p0.is_nan() || self.p1.is_nan() || self.p2.is_nan()

    }

}

/// An iterator for quadratic beziers.

pub struct QuadBezIter {

    quad: QuadBez,

    ix: usize,

}

impl Shape for QuadBez {

    type PathElementsIter<'iter> = QuadBezIter;

    #[inline]

    fn path_elements(&self, _tolerance: f64) -> QuadBezIter {

        QuadBezIter { quad: *self, ix: 0 }

    }

    fn area(&self) -> f64 {

        0.0

    }

    #[inline]

    fn perimeter(&self, accuracy: f64) -> f64 {

        self.arclen(accuracy)

    }

    fn winding(&self, _pt: Point) -> i32 {

        0

    }

    #[inline]

    fn bounding_box(&self) -> Rect {

        ParamCurveExtrema::bounding_box(self)

    }

}

impl Iterator for QuadBezIter {

    type Item = PathEl;

    fn next(&mut self) -> Option<PathEl> {

        self.ix += 1;

        match self.ix {

            1 => Some(PathEl::MoveTo(self.quad.p0)),

            2 => Some(PathEl::QuadTo(self.quad.p1, self.quad.p2)),

            _ => None,

        }

    }

}

pub(crate) struct FlattenParams {

    a0: f64,

    a2: f64,

    u0: f64,

    uscale: f64,

    /// The number of `subdivisions * 2 * sqrt_tol`.

    pub(crate) val: f64,

}

/// An approximation to $\int (1 + 4x^2) ^ -0.25 dx$

///

/// This is used for flattening curves.

fn approx_parabola_integral(x: f64) -> f64 {

    const D: f64 = 0.67;

    x / (1.0 - D + (D.powi(4) + 0.25 * x * x).sqrt().sqrt())

}

/// An approximation to the inverse parabola integral.

fn approx_parabola_inv_integral(x: f64) -> f64 {

    const B: f64 = 0.39;

    x * (1.0 - B + (B * B + 0.25 * x * x).sqrt())

}

impl ParamCurve for QuadBez {

    #[inline]

    fn eval(&self, t: f64) -> Point {

        let mt = 1.0 - t;

        (self.p0.to_vec2() * (mt * mt)

            + (self.p1.to_vec2() * (mt * 2.0) + self.p2.to_vec2() * t) * t)

            .to_point()

    }

    fn subsegment(&self, range: Range<f64>) -> QuadBez {

        let (t0, t1) = (range.start, range.end);

        let p0 = self.eval(t0);

        let p2 = self.eval(t1);

        let p1 = p0 + (self.p1 - self.p0).lerp(self.p2 - self.p1, t0) * (t1 - t0);

        QuadBez { p0, p1, p2 }

    }

    /// Subdivide into halves, using de Casteljau.

    #[inline]

    fn subdivide(&self) -> (QuadBez, QuadBez) {

        let pm = self.eval(0.5);

        (

            QuadBez::new(self.p0, self.p0.midpoint(self.p1), pm),

            QuadBez::new(pm, self.p1.midpoint(self.p2), self.p2),

        )

    }

    #[inline]

    fn start(&self) -> Point {

        self.p0

    }

    #[inline]

    fn end(&self) -> Point {

        self.p2

    }

}

impl ParamCurveDeriv for QuadBez {

    type DerivResult = Line;

    #[inline]

    fn deriv(&self) -> Line {

        Line::new(

            (2.0 * (self.p1.to_vec2() - self.p0.to_vec2())).to_point(),

            (2.0 * (self.p2.to_vec2() - self.p1.to_vec2())).to_point(),

        )

    }

}

impl ParamCurveArclen for QuadBez {

    /// Arclength of a quadratic Bézier segment.

    ///

    /// This computation is based on an analytical formula. Since that formula suffers

    /// from numerical instability when the curve is very close to a straight line, we

    /// detect that case and fall back to Legendre-Gauss quadrature.

    ///

    /// Accuracy should be better than 1e-13 over the entire range.

    ///

    /// Adapted from <http://www.malczak.linuxpl.com/blog/quadratic-bezier-curve-length/>

    /// with permission.

    fn arclen(&self, _accuracy: f64) -> f64 {

        let d2 = self.p0.to_vec2() - 2.0 * self.p1.to_vec2() + self.p2.to_vec2();

        let a = d2.hypot2();

        let d1 = self.p1 - self.p0;

        let c = d1.hypot2();

        if a < 5e-4 * c {

            // This case happens for nearly straight Béziers.

            //

            // Calculate arclength using Legendre-Gauss quadrature using formula from Behdad

            // in https://github.com/Pomax/BezierInfo-2/issues/77

            let v0 = (-0.492943519233745 * self.p0.to_vec2()

                + 0.430331482911935 * self.p1.to_vec2()

                + 0.0626120363218102 * self.p2.to_vec2())

            .hypot();

            let v1 = ((self.p2 - self.p0) * 0.4444444444444444).hypot();

            let v2 = (-0.0626120363218102 * self.p0.to_vec2()

                - 0.430331482911935 * self.p1.to_vec2()

                + 0.492943519233745 * self.p2.to_vec2())

            .hypot();

            return v0 + v1 + v2;

        }

        let b = 2.0 * d2.dot(d1);

        let sabc = (a + b + c).sqrt();

        let a2 = a.powf(-0.5);

        let a32 = a2.powi(3);

        let c2 = 2.0 * c.sqrt();

        let ba_c2 = b * a2 + c2;

        let v0 = 0.25 * a2 * a2 * b * (2.0 * sabc - c2) + sabc;

        // TODO: justify and fine-tune this exact constant.

        // The factor of a2 here is a little arbitrary: we really want

        // to test whether ba_c2 is small, but it's also important for

        // this comparison to be scale-invariant. We chose a2 (instead of,

        // for example, c2 on the rhs) because it's unchanged under

        // reversing the parametrization.

        if ba_c2 * a2 < 1e-13 {

            // This case happens for Béziers with a sharp kink.

            v0

        } else {

            v0 + 0.25

                * a32

                * (4.0 * c * a - b * b)

                * (((2.0 * a + b) * a2 + 2.0 * sabc) / ba_c2).ln()

        }

    }

}

impl ParamCurveArea for QuadBez {

    #[inline]

    fn signed_area(&self) -> f64 {

        (self.p0.x * (2.0 * self.p1.y + self.p2.y) + 2.0 * self.p1.x * (self.p2.y - self.p0.y)

            - self.p2.x * (self.p0.y + 2.0 * self.p1.y))

            * (1.0 / 6.0)

    }

}

impl ParamCurveNearest for QuadBez {

    /// Find the nearest point, using analytical algorithm based on cubic root finding.

    fn nearest(&self, p: Point, _accuracy: f64) -> Nearest {

        fn eval_t(p: Point, t_best: &mut f64, r_best: &mut Option<f64>, t: f64, p0: Point) {

            let r = (p0 - p).hypot2();

            if r_best.map(|r_best| r < r_best).unwrap_or(true) {

                *r_best = Some(r);

                *t_best = t;

            }

        }

        fn try_t(

            q: &QuadBez,

            p: Point,

            t_best: &mut f64,

            r_best: &mut Option<f64>,

            t: f64,

        ) -> bool {

            if !(0.0..=1.0).contains(&t) {

                return true;

            }

            eval_t(p, t_best, r_best, t, q.eval(t));

            false

        }

        let d0 = self.p1 - self.p0;

        let d1 = self.p0.to_vec2() + self.p2.to_vec2() - 2.0 * self.p1.to_vec2();

        let d = self.p0 - p;

        let c0 = d.dot(d0);

        let c1 = 2.0 * d0.hypot2() + d.dot(d1);

        let c2 = 3.0 * d1.dot(d0);

        let c3 = d1.hypot2();

        let roots = solve_cubic(c0, c1, c2, c3);

        let mut r_best = None;

        let mut t_best = 0.0;

        let mut need_ends = false;

        if roots.is_empty() {

            need_ends = true;

        }

        for &t in &roots {

            need_ends |= try_t(self, p, &mut t_best, &mut r_best, t);

        }

        if need_ends {

            eval_t(p, &mut t_best, &mut r_best, 0.0, self.p0);

            eval_t(p, &mut t_best, &mut r_best, 1.0, self.p2);

        }

        Nearest {

            t: t_best,

            distance_sq: r_best.unwrap(),

        }

    }

}

impl ParamCurveCurvature for QuadBez {}

impl ParamCurveExtrema for QuadBez {

    fn extrema(&self) -> ArrayVec<f64, MAX_EXTREMA> {

        let mut result = ArrayVec::new();

        let d0 = self.p1 - self.p0;

        let d1 = self.p2 - self.p1;

        let dd = d1 - d0;

        if dd.x != 0.0 {

            let t = -d0.x / dd.x;

            if t > 0.0 && t < 1.0 {

                result.push(t);

            }

        }

        if dd.y != 0.0 {

            let t = -d0.y / dd.y;

            if t > 0.0 && t < 1.0 {

                result.push(t);

                if result.len() == 2 && result[0] > t {

                    result.swap(0, 1);

                }

            }

        }

        result

    }

}

impl Mul<QuadBez> for Affine {

    type Output = QuadBez;

    #[inline]

    fn mul(self, other: QuadBez) -> QuadBez {

        QuadBez {

            p0: self * other.p0,

            p1: self * other.p1,

            p2: self * other.p2,

        }

    }

}

#[cfg(test)]

mod tests {

    use crate::{

        Affine, Nearest, ParamCurve, ParamCurveArclen, ParamCurveArea, ParamCurveDeriv,

        ParamCurveExtrema, ParamCurveNearest, Point, QuadBez,

    };

    fn assert_near(p0: Point, p1: Point, epsilon: f64) {

        assert!((p1 - p0).hypot() < epsilon, "{p0:?} != {p1:?}");

    }

    #[test]

    fn quadbez_deriv() {

        let q = QuadBez::new((0.0, 0.0), (0.0, 0.5), (1.0, 1.0));

        let deriv = q.deriv();

        let n = 10;

        for i in 0..=n {

            let t = (i as f64) * (n as f64).recip();

            let delta = 1e-6;

            let p = q.eval(t);

            let p1 = q.eval(t + delta);

            let d_approx = (p1 - p) * delta.recip();

            let d = deriv.eval(t).to_vec2();

            assert!((d - d_approx).hypot() < delta * 2.0);

        }

    }

    #[test]

    fn quadbez_arclen() {

        let q = QuadBez::new((0.0, 0.0), (0.0, 0.5), (1.0, 1.0));

        let true_arclen = 0.5 * 5.0f64.sqrt() + 0.25 * (2.0 + 5.0f64.sqrt()).ln();

        for i in 0..12 {

            let accuracy = 0.1f64.powi(i);

            let est = q.arclen(accuracy);

            let error = est - true_arclen;

            assert!(error.abs() < accuracy, "{est} != {true_arclen}");

        }

    }

    #[test]

    fn quadbez_arclen_pathological() {

        let q = QuadBez::new((-1.0, 0.0), (1.03, 0.0), (1.0, 0.0));

        let true_arclen = 2.0008737864167325; // A rough empirical calculation

        let accuracy = 1e-11;

        let est = q.arclen(accuracy);

        assert!(

            (est - true_arclen).abs() < accuracy,

            "{est} != {true_arclen}"

        );

    }

    #[test]

    fn quadbez_subsegment() {

        let q = QuadBez::new((3.1, 4.1), (5.9, 2.6), (5.3, 5.8));

        let t0 = 0.1;

        let t1 = 0.8;

        let qs = q.subsegment(t0..t1);

        let epsilon = 1e-12;

        let n = 10;

        for i in 0..=n {

            let t = (i as f64) * (n as f64).recip();

            let ts = t0 + t * (t1 - t0);

            assert_near(q.eval(ts), qs.eval(t), epsilon);

        }

    }

    #[test]

    fn quadbez_raise() {

        let q = QuadBez::new((3.1, 4.1), (5.9, 2.6), (5.3, 5.8));

        let c = q.raise();

        let qd = q.deriv();

        let cd = c.deriv();

        let epsilon = 1e-12;

        let n = 10;

        for i in 0..=n {

            let t = (i as f64) * (n as f64).recip();

            assert_near(q.eval(t), c.eval(t), epsilon);

            assert_near(qd.eval(t), cd.eval(t), epsilon);

        }

    }

    #[test]

    fn quadbez_signed_area() {

        // y = 1 - x^2

        let q = QuadBez::new((1.0, 0.0), (0.5, 1.0), (0.0, 1.0));

        let epsilon = 1e-12;

        assert!((q.signed_area() - 2.0 / 3.0).abs() < epsilon);

        assert!(((Affine::rotate(0.5) * q).signed_area() - 2.0 / 3.0).abs() < epsilon);

        assert!(((Affine::translate((0.0, 1.0)) * q).signed_area() - 3.5 / 3.0).abs() < epsilon);

        assert!(((Affine::translate((1.0, 0.0)) * q).signed_area() - 3.5 / 3.0).abs() < epsilon);

    }

    fn verify(result: Nearest, expected: f64) {

        assert!(

            (result.t - expected).abs() < 1e-6,

            "got {result:?} expected {expected}"

        );

    }

    #[test]

    fn quadbez_nearest() {

        // y = x^2

        let q = QuadBez::new((-1.0, 1.0), (0.0, -1.0), (1.0, 1.0));

        verify(q.nearest((0.0, 0.0).into(), 1e-3), 0.5);

        verify(q.nearest((0.0, 0.1).into(), 1e-3), 0.5);

        verify(q.nearest((0.0, -0.1).into(), 1e-3), 0.5);

        verify(q.nearest((0.5, 0.25).into(), 1e-3), 0.75);

        verify(q.nearest((1.0, 1.0).into(), 1e-3), 1.0);

        verify(q.nearest((1.1, 1.1).into(), 1e-3), 1.0);

        verify(q.nearest((-1.1, 1.1).into(), 1e-3), 0.0);

        let a = Affine::rotate(0.5);

        verify((a * q).nearest(a * Point::new(0.5, 0.25), 1e-3), 0.75);

    }

    // This test exposes a degenerate case in the solver used internally

    // by the "nearest" calculation - the cubic term is zero.

    #[test]

    fn quadbez_nearest_low_order() {

        let q = QuadBez::new((-1.0, 0.0), (0.0, 0.0), (1.0, 0.0));

        verify(q.nearest((0.0, 0.0).into(), 1e-3), 0.5);

        verify(q.nearest((0.0, 1.0).into(), 1e-3), 0.5);

    }

    #[test]

    fn quadbez_nearest_rounding_panic() {

        let quad = QuadBez::new(

            (-1.0394736842105263, 0.0),

            (0.8210526315789474, -1.511111111111111),

            (0.0, 1.9333333333333333),

        );

        let test = Point::new(-1.7976931348623157e308, 0.8571428571428571);

        // accuracy ignored

        let _res = quad.nearest(test, 1e-6);

        // if we got here then we didn't panic

    }

    #[test]

    fn quadbez_extrema() {

        // y = x^2

        let q = QuadBez::new((-1.0, 1.0), (0.0, -1.0), (1.0, 1.0));

        let extrema = q.extrema();

        assert_eq!(extrema.len(), 1);

        assert!((extrema[0] - 0.5).abs() < 1e-6);

        let q = QuadBez::new((0.0, 0.5), (1.0, 1.0), (0.5, 0.0));

        let extrema = q.extrema();

        assert_eq!(extrema.len(), 2);

        assert!((extrema[0] - 1.0 / 3.0).abs() < 1e-6);

        assert!((extrema[1] - 2.0 / 3.0).abs() < 1e-6);

        // Reverse direction

        let q = QuadBez::new((0.5, 0.0), (1.0, 1.0), (0.0, 0.5));

        let extrema = q.extrema();

        assert_eq!(extrema.len(), 2);

        assert!((extrema[0] - 1.0 / 3.0).abs() < 1e-6);

        assert!((extrema[1] - 2.0 / 3.0).abs() < 1e-6);

    }

    // A regression test for #477: the approximate-linearity test for

    // using the analytic solution needs to be scale-invariant.

    #[test]

    fn perimeter_not_nan() {

        let q = QuadBez::new((2685., -1251.), (2253., -1303.), (2253., -1303.));

        let len = q.arclen(crate::DEFAULT_ACCURACY);

        assert!(len.is_finite());

    }

}