klyff_msdf 0.1.3

#![allow(dead_code)] //why does expect produce warning here? there are dead code

//! Math utilities for line segments.
//!
//! Based on chapter 2 of chlumsky's thesis.

use std::f64::consts::PI;

use super::math_curve::cross;
use glam::{Vec2, vec2};

#[derive(Clone, Copy, PartialEq)]
pub(crate) struct Segment {
    pub(crate) from: Vec2,
    pub(crate) to: Vec2,
    pub(crate) color: SegmentColor,
    pub(crate) prev_angle: Option<i16>,

    /// Quadratic bezier mid point.
    /// All bezier edges in this context are flattened to very short and almost straight lines.
    /// This straight line is used in most calculation (closest segment, polygon clipping, etc.).
    /// The final MSDF calculation step for each pixel uses this quad bezier midpoint for
    /// additional accuracy.
    pub(crate) mid: Vec2,
}

#[derive(Clone, Copy, PartialEq)]
pub(crate) enum SegmentColor {
    White,
    RedGreen,
    GreenBlue,
    BlueRed,
}

#[derive(Clone, Copy, PartialEq)]
pub enum SegmentValidity {
    Unknown,
    Valid,
    Invalid,
}

impl Segment {
    pub(crate) fn has_color(&self, channel: usize) -> bool {
        match self.color {
            SegmentColor::White => true,
            SegmentColor::RedGreen => channel != 2,
            SegmentColor::GreenBlue => channel != 0,
            SegmentColor::BlueRed => channel != 1,
        }
    }

    pub(crate) fn angle(&self) -> i16 {
        angle_of_segment(self.to - self.from)
    }

    pub(crate) fn split_at(&self, pos: Vec2) -> (Segment, Segment) {
        fn interpolate(a: Vec2, b: Vec2, t: f32) -> Vec2 {
            a + t * (b - a)
        }

        let seg1 = Segment {
            to: pos,
            mid: self.from.midpoint(pos),
            ..*self
        };
        let seg2 = Segment {
            from: pos,
            mid: pos.midpoint(self.to),
            prev_angle: Some(seg1.angle()),
            ..*self
        };
        (seg1, seg2)
    }
}

pub(crate) fn angle_of_segment(dir: Vec2) -> i16 {
    let atan = f32::atan2(dir.y, dir.x);
    (atan / std::f32::consts::PI * (i16::MAX as f32)) as i16
}

impl std::fmt::Debug for Segment {
    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
        write!(
            f,
            "Segment({:.2},{:.2} : {:.2},{:.2} : {:.2},{:.2} : {:?}) ",
            self.from.x, self.from.y, self.mid.x, self.mid.y, self.to.x, self.to.y, self.color,
        )
    }
}
impl std::fmt::Debug for SegmentColor {
    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
        match self {
            SegmentColor::White => write!(f, "rgb"),
            SegmentColor::RedGreen => write!(f, "rg."),
            SegmentColor::GreenBlue => write!(f, ".gb"),
            SegmentColor::BlueRed => write!(f, "r.b"),
        }
    }
}

/// A bounding box.
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct Aabb {
    pub min: Vec2,
    pub max: Vec2,
}

impl Aabb {
    pub fn new(min_x: f32, min_y: f32, max_x: f32, max_y: f32) -> Self {
        Self {
            min: vec2(min_x, min_y),
            max: vec2(max_x, max_y),
        }
    }
    /// Creates from a line segment.
    pub(crate) fn from_segment(seg: Segment) -> Self {
        Self {
            min: vec2(seg.from.x.min(seg.to.x), seg.from.y.min(seg.to.y)),
            max: vec2(seg.from.x.max(seg.to.x), seg.from.y.max(seg.to.y)),
        }
    }

    /// Check if a bounding box is wholely contained inside this bounding box.
    pub fn contains(&self, other: &Self) -> bool {
        let xr = self.min.x..=self.max.x;
        let yr = self.min.y..=self.max.y;
        xr.contains(&other.min.x)
            && xr.contains(&other.max.x)
            && yr.contains(&other.min.y)
            && yr.contains(&other.max.y)
    }

    /// Check if a point belong to this bounding box.
    pub fn contains_pos(&self, pos: Vec2) -> bool {
        let xr = self.min.x..=self.max.x;
        let yr = self.min.y..=self.max.y;
        xr.contains(&pos.x) && yr.contains(&pos.y)
    }

    /// Corner with maximum x and minimum y component.
    pub fn maxx_miny(&self) -> Vec2 {
        vec2(self.max.x, self.min.y)
    }

    /// Corner with minimum x and maximum y component.
    pub fn minx_maxy(&self) -> Vec2 {
        vec2(self.min.x, self.max.y)
    }

    pub(crate) fn overlap(&self, other: &Aabb) -> bool {
        fn intersect(afrom: f32, ato: f32, bfrom: f32, bto: f32) -> bool {
            ato >= bfrom && afrom <= bto
        }
        intersect(self.min.x, self.max.x, other.min.x, other.max.x)
            && intersect(self.min.y, self.max.y, other.min.y, other.max.y)
    }
}

const EPS: f32 = 1e-6;
// Cardano's method
fn solve_cubic_polynomial(
    a: f64,
    b: f64,
    c: f64,
    d: f64,
) -> (Option<f64>, Option<f64>, Option<f64>) {
    // Normalize: x³ + ax² + bx + c
    let (a, b, c) = (b / a, c / a, d / a);

    let a2 = a * a;
    let q = 1. / 9. * (a2 - 3. * b);
    let r = 1. / 54. * (a * (2. * a2 - 9. * b) + 27. * c);
    let r2 = r * r;
    let q3 = q * q * q;
    let a = a / 3.;
    if r2 < q3 {
        let t = (r / q3.sqrt()).clamp(-1., 1.);
        let t = t.acos();
        let q = -2. * q.sqrt();
        (
            Some(q * (t / 3.).cos() - a),
            Some(q * ((t + 2. * PI) / 3.).cos() - a),
            Some(q * ((t - 2. * PI) / 3.).cos() - a),
        )
    } else {
        let u = (r.abs() + (r2 - q3).sqrt()).cbrt() * if r < 0. { 1. } else { -1. };
        let v = if u == 0. { 0. } else { q / u };
        let x0 = u + v - a;
        if u == v || (u - v).abs() < (u + v).abs() * 1e-12 {
            let x1 = (u + v) * -0.5 - a;
            (Some(x0), Some(x1), None)
        } else {
            (Some(x0), None, None)
        }
    }
}

/// Evaluates a position on a quad bezier.
fn quad_bezier(p0: Vec2, p1: Vec2, p2: Vec2, t: f32) -> Vec2 {
    p0 + 2. * t * (p1 - p0) + t * t * (p2 - 2. * p1 + p0)
}

/// Evaluates a tangent on a quad bezier.
fn d_quad_bezier(p0: Vec2, p1: Vec2, p2: Vec2, t: f32) -> Vec2 {
    2. * (p1 - p0) + 2. * t * (p2 - 2. * p1 + p0)
}

/// True signed distance from a point to a segment.
///
/// Chlumsky's thesis, Formula 2.41
pub(crate) fn signed_dist_point_to_segment(pos: Vec2, from: Vec2, mid: Vec2, to: Vec2) -> f32 {
    let p = pos - from;
    let p1 = mid - from;
    let p2 = to - 2. * mid + from;
    if p2.dot(p2) <= EPS {
        let tangent = to - from;
        if tangent.length_squared() <= 1e-6 {
            let perpendicular = pos - from.midpoint(to);
            return cross(tangent, perpendicular).signum() * (pos - from).length();
        }
        // Can NOT use clamped t for sign calculation here, since the cross product might flip sign
        // when the point is outside the [0, 1] range.
        let t = (pos - from).dot(tangent) / tangent.length_squared();
        let t_clamp = t.clamp(0., 1.);
        let projected = from + (to - from) * t;
        let projected_clamp = from + (to - from) * t_clamp;
        let perpendicular = projected - pos;
        let shortest = projected_clamp - pos;
        cross(tangent, perpendicular).signum() * shortest.length()
    } else {
        let (t1, t2, t3) = solve_cubic_polynomial(
            dot_f64(p2, p2),
            3. * dot_f64(p1, p2),
            2. * dot_f64(p1, p1) - dot_f64(p2, p),
            -dot_f64(p1, p),
        );

        let mut min_dist = f32::MAX;
        for t in [t1, t2, t3, Some(0.), Some(1.)].into_iter().flatten() {
            if !(0. ..=1.).contains(&t) {
                continue;
            }

            let perpendicular = quad_bezier(from, mid, to, t as f32) - pos;
            let dist = perpendicular.length();
            if dist < min_dist.abs() {
                min_dist =
                    dist * cross(d_quad_bezier(from, mid, to, t as f32), perpendicular).signum();
            }
        }

        min_dist
    }
}

fn dot_f64(a: Vec2, b: Vec2) -> f64 {
    a.x as f64 * b.x as f64 + a.y as f64 * b.y as f64
}

/// Signed pseudo distance from a point to a line segment.
///
/// This distance is only meaningful if taken to the closest line segment.
/// See Chlumsky's thesis, Formula 2.28 and section 2.4
pub(crate) fn pseudo_signed_dist_point_to_segment(
    pos: Vec2,
    from: Vec2,
    mid: Vec2,
    to: Vec2,
) -> f32 {
    let p = pos - from;
    let p1 = mid - from;
    let p2 = to - 2. * mid + from;
    if p2.dot(p2) <= EPS {
        let tangent = to - from;
        if tangent.length_squared() <= 1e-6 {
            let perpendicular = pos - from.midpoint(to);
            return cross(tangent, perpendicular).signum() * (pos - from).length();
        }
        let t = (pos - from).dot(tangent) / tangent.length_squared();
        let projected = from + (to - from) * t;
        let perpendicular = projected - pos;
        cross(tangent, perpendicular).signum() * perpendicular.length()
    } else {
        let (t1, t2, t3) = solve_cubic_polynomial(
            dot_f64(p2, p2),
            3. * dot_f64(p1, p2),
            2. * dot_f64(p1, p1) - dot_f64(p2, p),
            -dot_f64(p1, p),
        );

        let mut min_dist = f32::MAX;
        for t in [t1, t2, t3].into_iter().flatten() {
            if t < 0. {
                let tangent = mid - from;
                let perpendicular = if tangent.length_squared() <= 1e-6 {
                    pos - from.midpoint(mid)
                } else {
                    let t = (pos - from).dot(tangent) / tangent.length_squared();
                    let projected = from + tangent * t;
                    projected - pos
                };
                let dist = perpendicular.length();
                if dist < min_dist.abs() {
                    min_dist = dist * cross(tangent, perpendicular).signum();
                }
            } else if t > 1. {
                let tangent = to - mid;
                let perpendicular = if tangent.length_squared() <= 1e-6 {
                    pos - mid.midpoint(to)
                } else {
                    let t = (pos - mid).dot(tangent) / tangent.length_squared();
                    let projected = mid + tangent * t;
                    projected - pos
                };
                let dist = perpendicular.length();
                if dist < min_dist.abs() {
                    min_dist = dist * cross(tangent, perpendicular).signum();
                }
            } else {
                let tangent = d_quad_bezier(from, mid, to, t as f32);
                let perpendicular = quad_bezier(from, mid, to, t as f32) - pos;
                debug_assert!(
                    tangent.dot(perpendicular).abs() < 1.,
                    "{}",
                    tangent.dot(perpendicular).abs()
                );
                let dist = perpendicular.length();
                if dist < min_dist.abs() {
                    min_dist = dist * cross(tangent, perpendicular).signum();
                }
            }
        }

        min_dist
    }
}

/// Unsigned, true distance from a point to a line segment, squared.
pub(crate) fn unsigned_dist_point_to_segment_sqr(pos: Vec2, from: Vec2, to: Vec2) -> f32 {
    let tangent = to - from;
    if tangent.length_squared() <= 1e-6 {
        return (pos - from).length_squared();
    }
    let t = (pos - from).dot(tangent) / tangent.length_squared();
    let t = t.clamp(0., 1.);
    let projected = from + (to - from) * t;
    let shortest = projected - pos;
    shortest.length_squared()
}

/// Orthogonality measure from a point to a segment.
///
/// Used to disambiguate between segments with the same distance to a point.
/// See Chlumsky's thesis, Formula 2.43
pub(crate) fn orthogonality(pos: Vec2, from: Vec2, to: Vec2) -> f32 {
    let tangent = to - from;
    let t = (pos - from).dot(tangent) / tangent.length_squared();
    let t_clamp = t.clamp(0., 1.);
    let projected = from + (to - from) * t_clamp;
    cross(tangent.normalize(), (pos - projected).normalize()).abs()
}

/// Distance from a point to an infinite straight line.
pub(crate) fn distance_point_to_line(pos: Vec2, from: Vec2, to: Vec2) -> f32 {
    let v1 = to - from;
    let v2 = pos - from;
    let perp_dot = (v1.perp_dot(v2)).abs();
    perp_dot / v1.length()
}

/// Unsigned distance from a point to an AABB, squared.
pub(crate) fn dist_point_to_aabb_sqr(pos: Vec2, aabb: Aabb) -> f32 {
    // See RafiRome's answer:
    // https://stackoverflow.com/questions/5254838/calculating-distance-between-a-point-and-a-rectangular-box-nearest-point
    let dx = (aabb.min.x - pos.x).max(0.).max(pos.x - aabb.max.x);
    let dy = (aabb.min.y - pos.y).max(0.).max(pos.y - aabb.max.y);
    dx * dx + dy * dy
}

/// Calculate the AABB bounding box of some number of segments.
pub(crate) fn bound_of_segments(segs: &[Segment]) -> Aabb {
    let mut iter = segs.iter();
    let Some(first_seg) = iter.next() else {
        return Aabb {
            min: Vec2::ZERO,
            max: Vec2::ZERO,
        };
    };

    let mut aabb = Aabb {
        min: Vec2 {
            x: first_seg.from.x.min(first_seg.to.x),
            y: first_seg.from.y.min(first_seg.to.y),
        },
        max: Vec2 {
            x: first_seg.from.x.max(first_seg.to.x),
            y: first_seg.from.y.max(first_seg.to.y),
        },
    };

    for seg in iter {
        aabb.min.x = aabb.min.x.min(f32::min(seg.from.x, seg.to.x));
        aabb.max.x = aabb.max.x.max(f32::max(seg.from.x, seg.to.x));
        aabb.min.y = aabb.min.y.min(f32::min(seg.from.y, seg.to.y));
        aabb.max.y = aabb.max.y.max(f32::max(seg.from.y, seg.to.y));
    }

    aabb
}

#[cfg(test)]
use proptest::prelude::*;
#[cfg(test)]
impl Segment {
    pub(crate) fn arbitrary_poly(
        reverse_probability: f64,
        edge_count_range: std::ops::Range<usize>,
    ) -> impl Strategy<Value = Vec<Self>> {
        let coord_range = -100.0..=100.0f32;
        let radius_range = 10.0..=200.0f32;
        (
            proptest::bool::weighted(reverse_probability),
            coord_range.clone(),
            coord_range,
            0.0f32..=1.0,
            0.0f32..=1.0,
            radius_range,
            edge_count_range,
        )
            .prop_flat_map(
                |(reverse, center_x, center_y, irregularity, spikiness, radius, edge_count)| {
                    let irregularity = irregularity.clamp(0.2, 0.8) * 2. * std::f32::consts::PI
                        / (edge_count as f32);
                    let spikiness = radius * spikiness.clamp(0.2, 0.8);
                    let angle_step = 2. * std::f32::consts::PI / (edge_count as f32);
                    proptest::collection::vec(
                        (
                            -irregularity..=irregularity,
                            radius - spikiness..=radius + spikiness,
                        ),
                        edge_count..=edge_count,
                    )
                    .prop_map(move |offsets| {
                        offsets
                            .into_iter()
                            .enumerate()
                            .map(|(i, (angle_offset, radius_offset))| {
                                let r = radius + radius_offset;
                                let a = if reverse {
                                    (angle_step * i as f32) + angle_offset
                                } else {
                                    (-angle_step * i as f32) + angle_offset
                                };
                                vec2(center_x + r * f32::cos(a), center_y + r * f32::sin(a))
                            })
                            .collect::<Vec<_>>()
                    })
                },
            )
            .prop_map(|verts| {
                use crate::{math_curve::is_corner, segment::cycle_color};

                let mut segs: Vec<Segment> = vec![];
                let wrap = [*verts.last().unwrap(), *verts.first().unwrap()];
                let mut c = SegmentColor::BlueRed;
                let mut angle = Segment {
                    from: wrap[0],
                    to: wrap[1],
                    color: c,
                    prev_angle: None,
                    mid: wrap[0].midpoint(wrap[1]),
                }
                .angle();
                for seg in verts.windows(2).chain(std::iter::once(wrap.as_slice())) {
                    let from = seg[0];
                    let to = seg[1];

                    let prev_angle = if let Some(last) = segs.last()
                        && is_corner(last.to - last.from, to - from)
                    {
                        cycle_color(&mut c);
                        None
                    } else {
                        Some(angle)
                    };
                    let seg = Segment {
                        from,
                        mid: from.midpoint(to),
                        to,
                        color: c,
                        prev_angle,
                    };
                    segs.push(seg);
                    angle = seg.angle();
                }
                segs
            })
    }

    pub(crate) fn arbitrary_poly_multi(
        reverse_probability: f64,
        shape_count_range: std::ops::Range<usize>,
        edge_count_range: std::ops::Range<usize>,
    ) -> impl Strategy<Value = Vec<Segment>> {
        proptest::collection::vec(
            Segment::arbitrary_poly(reverse_probability, edge_count_range),
            shape_count_range,
        )
        .prop_map(|polys| {
            let mut combined = vec![];
            for mut poly in polys {
                combined.append(&mut poly);
            }
            combined
        })
    }
}