s2rst 0.4.0

// SPDX-License-Identifier: Apache-2.0
// SPDX-FileCopyrightText: 2026 Torgeir Børresen <tb@starkad.no>
// Rust port of Google's S2 Geometry library — a derivative work, modified from
// the upstream Apache-2.0 source(s) below (Copyright Google Inc.). See LICENSE.
//   - C++:  google/s2geometry
//   - Java: google/s2-geometry-library-java

//! [`S2Fractal`] generates fractal loops on the unit sphere.
//!
//! Used primarily as a testing utility to create complex, realistic
//! geometry (e.g., Koch-snowflake-like coastlines) at configurable
//! subdivision levels and fractal dimensions.
//!
//! Corresponds to C++ `s2fractal.h/cc`.

#![expect(
    clippy::cast_sign_loss,
    reason = "level (i32) cast to u32 for pow — always non-negative"
)]
#![expect(
    clippy::cast_possible_truncation,
    reason = "level (i32->u32) for pow — always non-negative"
)]
use crate::r3::Matrix3x3;
use crate::s1::Angle;
use crate::s2::Loop;
use crate::s2::point::{Point, from_frame, get_frame};

/// A simple xorshift64 pseudo-random number generator.
///
/// Used internally by [`S2Fractal`] for deterministic fractal generation.
#[derive(Debug)]
struct Rng {
    state: u64,
}

impl Rng {
    fn new(seed: u64) -> Self {
        Rng {
            state: if seed == 0 {
                0x_dead_beef_cafe_babe
            } else {
                seed
            },
        }
    }

    /// Returns a pseudo-random `u64`.
    fn next_u64(&mut self) -> u64 {
        let mut x = self.state;
        x ^= x << 13;
        x ^= x >> 7;
        x ^= x << 17;
        self.state = x;
        x
    }

    /// Returns a pseudo-random `f64` in `[0, 1)`.
    fn next_f64(&mut self) -> f64 {
        (self.next_u64() >> 11) as f64 / (1u64 << 53) as f64
    }

    /// Returns `true` with probability `p`.
    fn bernoulli(&mut self, p: f64) -> bool {
        self.next_f64() < p
    }
}

/// A 2D point in the tangent plane, used during fractal generation.
#[derive(Clone, Copy)]
struct R2Point {
    x: f64,
    y: f64,
}

impl R2Point {
    fn new(x: f64, y: f64) -> Self {
        R2Point { x, y }
    }

    /// Returns the perpendicular vector (rotated 90 degrees counter-clockwise).
    fn ortho(self) -> R2Point {
        R2Point::new(-self.y, self.x)
    }
}

impl std::ops::Add for R2Point {
    type Output = R2Point;
    fn add(self, rhs: R2Point) -> R2Point {
        R2Point::new(self.x + rhs.x, self.y + rhs.y)
    }
}

impl std::ops::Sub for R2Point {
    type Output = R2Point;
    fn sub(self, rhs: R2Point) -> R2Point {
        R2Point::new(self.x - rhs.x, self.y - rhs.y)
    }
}

impl std::ops::Mul<f64> for R2Point {
    type Output = R2Point;
    fn mul(self, s: f64) -> R2Point {
        R2Point::new(self.x * s, self.y * s)
    }
}

impl std::ops::Mul<R2Point> for f64 {
    type Output = R2Point;
    fn mul(self, p: R2Point) -> R2Point {
        R2Point::new(self * p.x, self * p.y)
    }
}

/// Generates fractal loops on the unit sphere.
///
/// The fractal is a Koch-snowflake-like curve generated by recursively
/// subdividing edges of an equilateral triangle. The fractal dimension
/// controls how "jagged" the result is: 1.0 produces a simple triangle,
/// values approaching 2.0 produce nearly space-filling curves.
///
/// # Example
///
/// ```ignore
/// use s2rst::s2::fractal::S2Fractal;
/// use s2rst::s1::Angle;
/// use s2rst::s2::point::get_frame;
/// use s2rst::s2::Point;
///
/// let mut fractal = S2Fractal::new(42);
/// fractal.set_max_level(5);
/// fractal.set_fractal_dimension(1.5);
/// let frame = get_frame(Point::from_coords(0.0, 0.0, 1.0));
/// let loop_ = fractal.make_loop(&frame, Angle::from_degrees(10.0));
/// ```
#[derive(Debug)]
pub struct S2Fractal {
    rng: Rng,
    max_level: i32,
    min_level_arg: i32,
    min_level: i32,
    dimension: f64,
    edge_fraction: f64,
    offset_fraction: f64,
}

impl S2Fractal {
    /// Creates a new fractal generator with the given random seed.
    pub fn new(seed: u64) -> Self {
        let dimension = (4.0_f64).ln() / (3.0_f64).ln(); // ~1.2619
        let edge_fraction = 4.0_f64.powf(-1.0 / dimension);
        let offset_fraction = (edge_fraction - 0.25).sqrt();
        S2Fractal {
            rng: Rng::new(seed),
            max_level: -1,
            min_level_arg: -1,
            min_level: -1,
            dimension,
            edge_fraction,
            offset_fraction,
        }
    }

    /// Sets the maximum subdivision level.
    ///
    /// The number of vertices at a given level is `3 * 4^level`.
    /// Must be called before `make_loop()`.
    ///
    /// # Panics
    ///
    /// Panics if `max_level` is negative.
    pub fn set_max_level(&mut self, max_level: i32) {
        assert!(max_level >= 0);
        self.max_level = max_level;
        self.compute_min_level();
    }

    /// Returns the maximum subdivision level.
    pub fn max_level(&self) -> i32 {
        self.max_level
    }

    /// Sets the minimum subdivision level.
    ///
    /// If set to -1 (default), the minimum level equals the maximum level
    /// (uniform subdivision). If set to a value less than `max_level`, the
    /// fractal will have variable detail, with random early termination
    /// between `min_level` and `max_level`.
    ///
    /// # Panics
    ///
    /// Panics if `min_level_arg` is less than -1.
    pub fn set_min_level(&mut self, min_level_arg: i32) {
        assert!(min_level_arg >= -1);
        self.min_level_arg = min_level_arg;
        self.compute_min_level();
    }

    /// Returns the minimum subdivision level.
    pub fn min_level(&self) -> i32 {
        self.min_level
    }

    /// Sets the level such that the approximate number of edges is at least
    /// `min_edges`. The actual level is `round(0.5 * log2(min_edges / 3))`.
    pub fn level_for_approx_min_edges(&mut self, min_edges: i32) {
        let level = (0.5 * (f64::from(min_edges) / 3.0).log2()).round() as i32;
        self.set_min_level(level.max(0));
    }

    /// Sets the level such that the approximate number of edges is at most
    /// `max_edges`. The actual level is `round(0.5 * log2(max_edges / 3))`.
    pub fn level_for_approx_max_edges(&mut self, max_edges: i32) {
        let level = (0.5 * (f64::from(max_edges) / 3.0).log2()).round() as i32;
        self.set_max_level(level.max(0));
    }

    /// Sets the fractal dimension in range `[1.0, 2.0)`.
    ///
    /// - 1.0: simple triangle (no fractal subdivision)
    /// - log(4)/log(3) ≈ 1.26: Koch snowflake
    /// - approaching 2.0: nearly space-filling
    ///
    /// # Panics
    ///
    /// Panics if `dimension` is outside `[1.0, 2.0)`.
    pub fn set_fractal_dimension(&mut self, dimension: f64) {
        assert!((1.0..2.0).contains(&dimension));
        self.dimension = dimension;
        self.compute_offsets();
    }

    /// Returns the fractal dimension.
    pub fn fractal_dimension(&self) -> f64 {
        self.dimension
    }

    /// Returns the minimum distance from the fractal center to any boundary
    /// point, as a factor of the nominal radius.
    pub fn min_radius_factor(&self) -> f64 {
        // The minimum radius is attained at one of the vertices created by the
        // first subdivision step as long as the dimension is not too small (at
        // least kMinDimensionForMinRadiusAtLevel1, see below).  Otherwise we
        // fall back on the incircle radius of the original triangle, which is
        // always a lower bound (and is attained when dimension = 1).
        //
        // The value below is equal to -log(4)/log((2 + cbrt(2) - cbrt(4))/6).
        const MIN_DIMENSION_FOR_MIN_RADIUS_AT_LEVEL1: f64 = 1.0852230903040407;
        if self.dimension >= MIN_DIMENSION_FOR_MIN_RADIUS_AT_LEVEL1 {
            (1.0 + 3.0 * self.edge_fraction * (self.edge_fraction - 1.0)).sqrt()
        } else {
            0.5
        }
    }

    /// Returns the maximum distance from the fractal center to any boundary
    /// point, as a factor of the nominal radius.
    pub fn max_radius_factor(&self) -> f64 {
        1.0_f64.max(self.offset_fraction * 3.0_f64.sqrt() + 0.5)
    }

    /// Generates a fractal loop on the unit sphere.
    ///
    /// The loop is centered on the z-axis of `frame` with the given
    /// `nominal_radius`. The actual radius varies by
    /// `min_radius_factor()` to `max_radius_factor()`.
    pub fn make_loop(&mut self, frame: &Matrix3x3, nominal_radius: Angle) -> Loop {
        let r2_vertices = self.get_r2_vertices();
        let r = nominal_radius.radians();
        let mut vertices = Vec::with_capacity(r2_vertices.len());
        for v in &r2_vertices {
            let p = Point(crate::r3::Vector {
                x: v.x * r,
                y: v.y * r,
                z: 1.0,
            });
            vertices.push(from_frame(frame, p).normalize());
        }
        Loop::new(vertices)
    }

    /// Generates a fractal loop centered at the given point with the given radius.
    pub fn make_loop_at(&mut self, center: Point, nominal_radius: Angle) -> Loop {
        let frame = get_frame(center);
        self.make_loop(&frame, nominal_radius)
    }

    fn compute_min_level(&mut self) {
        if self.min_level_arg >= 0 && self.min_level_arg <= self.max_level {
            self.min_level = self.min_level_arg;
        } else {
            self.min_level = self.max_level;
        }
    }

    fn compute_offsets(&mut self) {
        self.edge_fraction = 4.0_f64.powf(-1.0 / self.dimension);
        self.offset_fraction = (self.edge_fraction - 0.25).sqrt();
    }

    fn get_r2_vertices(&mut self) -> Vec<R2Point> {
        let mut vertices = Vec::new();
        let sqrt3_2 = 3.0_f64.sqrt() / 2.0;
        let v0 = R2Point::new(1.0, 0.0);
        let v1 = R2Point::new(-0.5, sqrt3_2);
        let v2 = R2Point::new(-0.5, -sqrt3_2);

        self.get_r2_vertices_helper(v0, v1, 0, &mut vertices);
        self.get_r2_vertices_helper(v1, v2, 0, &mut vertices);
        self.get_r2_vertices_helper(v2, v0, 0, &mut vertices);

        vertices
    }

    fn get_r2_vertices_helper(
        &mut self,
        v0: R2Point,
        v4: R2Point,
        level: i32,
        vertices: &mut Vec<R2Point>,
    ) {
        if level >= self.min_level {
            let levels_remaining = self.max_level - level + 1;
            if self.rng.bernoulli(1.0 / f64::from(levels_remaining)) {
                vertices.push(v0);
                return;
            }
        }

        let dir = v4 - v0;
        let v1 = v0 + dir * self.edge_fraction;
        let v2 = 0.5 * (v0 + v4) - dir.ortho() * self.offset_fraction;
        let v3 = v4 - dir * self.edge_fraction;

        self.get_r2_vertices_helper(v0, v1, level + 1, vertices);
        self.get_r2_vertices_helper(v1, v2, level + 1, vertices);
        self.get_r2_vertices_helper(v2, v3, level + 1, vertices);
        self.get_r2_vertices_helper(v3, v4, level + 1, vertices);
    }
}

/// Returns the number of vertices at a given subdivision level.
///
/// At level `n`, the fractal has `3 * 4^n` edges (and vertices).
pub fn num_vertices_at_level(level: i32) -> usize {
    3 * (1usize << (2 * level as u32))
}

#[cfg(test)]
mod tests {
    use super::*;
    use crate::s1::Angle;
    use crate::s2::point::get_frame;

    #[test]
    fn test_num_vertices_at_level() {
        assert_eq!(num_vertices_at_level(0), 3);
        assert_eq!(num_vertices_at_level(1), 12);
        assert_eq!(num_vertices_at_level(2), 48);
        assert_eq!(num_vertices_at_level(3), 192);
        assert_eq!(num_vertices_at_level(4), 768);
    }

    #[test]
    fn test_default_dimension() {
        let f = S2Fractal::new(1);
        let expected = (4.0_f64).ln() / (3.0_f64).ln();
        assert!((f.fractal_dimension() - expected).abs() < 1e-14);
    }

    fn test_fractal(min_level: i32, max_level: i32, dimension: f64) {
        let mut fractal = S2Fractal::new(42);
        fractal.set_fractal_dimension(dimension);
        fractal.set_max_level(max_level);
        if min_level >= 0 {
            fractal.set_min_level(min_level);
        }

        let min_rf = fractal.min_radius_factor();
        let max_rf = fractal.max_radius_factor();

        let frame = get_frame(Point::from_coords(0.0, 0.0, 1.0));
        let radius = Angle::from_degrees(10.0);
        let lp = fractal.make_loop(&frame, radius);
        let nv = lp.num_vertices();

        // Verify vertex count is in expected range.
        let effective_min = if min_level >= 0 { min_level } else { max_level };
        let min_verts = num_vertices_at_level(effective_min);
        let max_verts = num_vertices_at_level(max_level);
        assert!(
            nv >= min_verts && nv <= max_verts,
            "vertex count {nv} not in range [{min_verts}, {max_verts}] \
             for levels [{effective_min}, {max_level}], dim={dimension}",
        );

        // Verify radius factors are positive and sensible.
        assert!(min_rf > 0.0, "min_radius_factor = {min_rf}");
        assert!(max_rf >= min_rf, "max_rf {max_rf} < min_rf {min_rf}");

        // Verify the loop is valid.
        assert!(lp.validate().is_ok(), "loop validation failed");
    }

    #[test]
    fn test_triangle_fractal() {
        // dimension=1.0 → simple triangle, no fractal subdivision
        test_fractal(-1, 5, 1.0);
    }

    #[test]
    fn test_triangle_multi_fractal() {
        test_fractal(2, 6, 1.0);
    }

    #[test]
    fn test_koch_curve_fractal() {
        // Standard Koch snowflake dimension.
        let dim = (4.0_f64).ln() / (3.0_f64).ln();
        test_fractal(-1, 5, dim);
    }

    #[test]
    fn test_koch_curve_multi_fractal() {
        let dim = (4.0_f64).ln() / (3.0_f64).ln();
        test_fractal(3, 6, dim);
    }

    #[test]
    fn test_space_filling_fractal() {
        // Nearly space-filling (dimension close to 2.0).
        test_fractal(-1, 3, 1.999);
    }

    #[test]
    fn test_cesaro_fractal() {
        test_fractal(-1, 5, 1.8);
    }

    #[test]
    fn test_cesaro_multi_fractal() {
        test_fractal(2, 5, 1.8);
    }

    #[test]
    fn test_make_loop_at() {
        let mut fractal = S2Fractal::new(123);
        fractal.set_max_level(3);
        let center = Point::from_coords(1.0, 0.0, 0.0);
        let lp = fractal.make_loop_at(center, Angle::from_degrees(5.0));

        assert!(lp.num_vertices() >= 3);
        assert!(lp.validate().is_ok());
    }

    #[test]
    fn test_level_for_approx_edges() {
        let mut fractal = S2Fractal::new(1);
        fractal.level_for_approx_max_edges(100);
        // 3 * 4^n ~ 100 → n ~ 2.5 → round to 3 → 192 edges max
        assert!(fractal.max_level() >= 2 && fractal.max_level() <= 4);

        fractal.level_for_approx_min_edges(50);
        assert!(fractal.min_level() >= 1 && fractal.min_level() <= 3);
    }

    #[test]
    fn test_radius_factors() {
        // For Koch snowflake, radius factors should be predictable.
        let dim = (4.0_f64).ln() / (3.0_f64).ln();
        let mut fractal = S2Fractal::new(1);
        fractal.set_fractal_dimension(dim);
        fractal.set_max_level(5);

        let min_rf = fractal.min_radius_factor();
        let max_rf = fractal.max_radius_factor();

        // Koch snowflake: min_rf ≈ 0.866 (√3/2), max_rf ≈ 1.115
        assert!(
            min_rf > 0.5 && min_rf < 1.0,
            "min_rf = {min_rf}, expected ~0.866"
        );
        assert!(
            (1.0..1.5).contains(&max_rf),
            "max_rf = {max_rf}, expected ~1.115"
        );
    }

    #[test]
    fn test_deterministic_with_same_seed() {
        let mut f1 = S2Fractal::new(42);
        f1.set_max_level(3);
        let frame = get_frame(Point::from_coords(0.0, 0.0, 1.0));
        let l1 = f1.make_loop(&frame, Angle::from_degrees(5.0));

        let mut f2 = S2Fractal::new(42);
        f2.set_max_level(3);
        let l2 = f2.make_loop(&frame, Angle::from_degrees(5.0));

        assert_eq!(l1.num_vertices(), l2.num_vertices());
        for i in 0..l1.num_vertices() {
            assert_eq!(l1.vertex(i), l2.vertex(i));
        }
    }

    #[test]
    fn test_different_seeds_differ() {
        let frame = get_frame(Point::from_coords(0.0, 0.0, 1.0));

        let mut f1 = S2Fractal::new(1);
        f1.set_max_level(4);
        f1.set_min_level(2);
        let l1 = f1.make_loop(&frame, Angle::from_degrees(5.0));

        let mut f2 = S2Fractal::new(999);
        f2.set_max_level(4);
        f2.set_min_level(2);
        let l2 = f2.make_loop(&frame, Angle::from_degrees(5.0));

        // Different seeds should (almost certainly) produce different loops.
        let same = l1.num_vertices() == l2.num_vertices()
            && (0..l1.num_vertices()).all(|i| l1.vertex(i) == l2.vertex(i));
        assert!(!same, "different seeds should produce different loops");
    }
}