i_curve 0.1.2

use crate::collections::stack_vec::StackVec;
use crate::int::CurveInt;
use crate::kernel::int::curve::cubic::CubicSegment;
use crate::kernel::int::curve::param::SegmentParam;
use crate::kernel::int::normalization::monotone::decomposition::MonotoneDecompositionDirection;
#[cfg(test)]
use crate::kernel::int::normalization::monotone::decomposition::{DecomposeIntoMonotone, roots_to_segments};
use crate::kernel::int::normalization::unit_quadratic::solve_unit_quadratic;
use i_overlay::i_float::int::number::wide_int::WideIntNumber;
use i_overlay::i_shape::int::IntPoint;

#[cfg(test)]
impl<I: CurveInt> DecomposeIntoMonotone for CubicSegment<I> {
    type Output = StackVec<CubicSegment<I>, 5>;

    fn decompose_into_monotone(&self) -> Self::Output {
        roots_to_segments(self, self.monotone_roots())
    }
}

impl<I: CurveInt> CubicSegment<I> {
    pub(crate) fn monotone_roots(&self) -> StackVec<SegmentParam<I>, 4> {
        let x_roots = self.monotone_roots_by_direction(MonotoneDecompositionDirection::X);
        let y_roots = self.monotone_roots_by_direction(MonotoneDecompositionDirection::Y);

        let mut roots = StackVec::new();
        roots.extend_from_slice(x_roots.as_slice());
        roots.extend_from_slice(y_roots.as_slice());
        roots.as_mut_slice().sort_unstable_by_key(|root| root.value());
        roots.dedup();
        roots
    }

    fn monotone_roots_by_direction(
        &self,
        direction: MonotoneDecompositionDirection,
    ) -> StackVec<SegmentParam<I>, 2> {
        let [a, b, c, _] = Self::axis_abcd(self.control_points, direction);

        // Cubic Bezier in Bernstein form:
        //
        //   q(t) = (1 - t)^3 * k0
        //        + 3 * (1 - t)^2 * t * k1
        //        + 3 * (1 - t) * t^2 * k2
        //        + t^3 * k3
        //
        // In power form:
        //
        //   q(t) = a * t^3 + b * t^2 + c * t + d
        //
        // A cubic is monotone on intervals that do not contain an extremum.
        // Extrema are roots of the derivative:
        //
        //   q'(t) = 3 * a * t^2 + 2 * b * t + c = 0
        //
        // The derivative is quadratic, so it can contribute up to two split
        // parameters per axis.
        solve_unit_quadratic::<I>(I::Wide::from_u32(3) * a, I::Wide::TWO * b, c)
    }

    #[inline]
    fn axis_abcd(
        control_points: [IntPoint<I>; 4],
        direction: MonotoneDecompositionDirection,
    ) -> [I::Wide; 4] {
        let [k0, k1, k2, k3] = match direction {
            MonotoneDecompositionDirection::X => control_points.map(|p| p.x.to_wide()),
            MonotoneDecompositionDirection::Y => control_points.map(|p| p.y.to_wide()),
        };

        // Expanding the Bernstein cubic gives:
        //
        //   q(t) = a * t^3 + b * t^2 + c * t + d
        //
        // where:
        //
        //   a = -k0 + 3 * k1 - 3 * k2 + k3
        //   b =  3 * k0 - 6 * k1 + 3 * k2
        //   c = -3 * k0 + 3 * k1
        //   d =  k0
        let three = I::Wide::from_u32(3);
        let six = I::Wide::from_u32(6);

        [
            k3 - three * k2 + three * k1 - k0,
            three * k0 - six * k1 + three * k2,
            three * (k1 - k0),
            k0,
        ]
    }
}

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

    #[test]
    fn splits_cubic_by_distinct_axis_extrema() {
        let cubic = CubicSegment::<i32> {
            control_points: [
                IntPoint::new(0, 0),
                IntPoint::new(4, 3),
                IntPoint::new(-1, -3),
                IntPoint::new(3, 0),
            ],
        };

        let parts = cubic.decompose_into_monotone();
        let parts = parts.as_slice();

        assert_eq!(parts.len(), 3);
        assert!(
            parts
                .iter()
                .all(|part| part.control_points[0] != part.control_points[3])
        );
        assert_eq!(parts[0].control_points[0], cubic.control_points[0]);
        assert_eq!(parts[2].control_points[3], cubic.control_points[3]);
    }

    #[test]
    fn keeps_fractional_derivative_roots_in_fixed_scale() {
        let cubic = CubicSegment::<i32> {
            control_points: [
                IntPoint::new(0, 0),
                IntPoint::new(0, 4),
                IntPoint::new(0, -1),
                IntPoint::new(0, 3),
            ],
        };

        let roots = cubic.monotone_roots_by_direction(MonotoneDecompositionDirection::Y);
        let roots = roots.as_slice();

        assert_eq!(roots.len(), 2);
        assert!((roots[0].value() - SegmentParam::<i32>::from_int(1, 3).value()).abs() <= 1i64);
        assert!((roots[1].value() - SegmentParam::<i32>::from_int(2, 3).value()).abs() <= 1i64);
    }

    #[test]
    fn keeps_irrational_derivative_roots_in_fixed_scale() {
        let cubic = CubicSegment::<i32> {
            control_points: [
                IntPoint::new(0, 0),
                IntPoint::new(0, 3),
                IntPoint::new(0, -3),
                IntPoint::new(0, 0),
            ],
        };

        let roots = cubic.monotone_roots_by_direction(MonotoneDecompositionDirection::Y);
        let roots = roots.as_slice();

        let scale = SegmentParam::<i32>::DENOMINATOR as f64;
        let offset = (1.0_f64 / 3.0).sqrt();
        let expected_0 = (scale * (1.0 - offset) * 0.5).round() as i64;
        let expected_1 = (scale * (1.0 + offset) * 0.5).round() as i64;

        assert_eq!(roots.len(), 2);
        assert!((roots[0].value() - expected_0).abs() <= 1i64);
        assert!((roots[1].value() - expected_1).abs() <= 1i64);
    }
}