open-standards/map/map-rs/crates/map-index/src/projection.rs

// S2 face projection math — quadratic projection (same as Google S2).

use crate::point::S2Point;

/// Maximum valid level for S2 cells.
pub const MAX_LEVEL: u8 = 30;

/// Number of cells along one axis at a given level.
#[inline]
pub const fn grid_size(level: u8) -> u32 {
    1u32 << level
}

/// Determine which face of the cube a point projects onto, and compute the UV coordinates.
/// Returns (face, u, v).
#[inline]
pub fn xyz_to_face_uv(p: &S2Point) -> (u8, f64, f64) {
    let face = largest_abs_component(p);
    let (u, v) = valid_face_xyz_to_uv(face, p);
    (face, u, v)
}

/// Convert face + UV coordinates back to an S2Point (not normalized).
#[inline]
pub fn face_uv_to_xyz(face: u8, u: f64, v: f64) -> S2Point {
    match face {
        0 => S2Point::new(1.0, u, v),
        1 => S2Point::new(-u, 1.0, v),
        2 => S2Point::new(-u, -v, 1.0),
        3 => S2Point::new(-1.0, -v, -u),
        4 => S2Point::new(v, -1.0, -u),
        5 => S2Point::new(v, u, -1.0),
        _ => S2Point::new(0.0, 0.0, 0.0),
    }
}

/// Quadratic ST-to-UV projection (S2 standard).
#[inline]
pub fn st_to_uv(s: f64) -> f64 {
    if s >= 0.5 {
        (1.0 / 3.0) * (4.0 * s * s - 1.0)
    } else {
        (1.0 / 3.0) * (1.0 - 4.0 * (1.0 - s) * (1.0 - s))
    }
}

/// Quadratic UV-to-ST projection (inverse of st_to_uv).
#[inline]
pub fn uv_to_st(u: f64) -> f64 {
    if u >= 0.0 {
        0.5 * (1.0 + 3.0 * u).sqrt()
    } else {
        1.0 - 0.5 * (1.0 - 3.0 * u).sqrt()
    }
}

/// Convert an ST coordinate to a discrete grid coordinate (I or J) at the given level.
/// The grid has 2^level cells along each axis, so I/J range is [0, 2^level - 1].
/// ST is in [0, 1].
#[inline]
pub fn st_to_ij(s: f64, level: u8) -> u32 {
    let size = grid_size(level) as f64;
    let ij = (s * size).floor() as i64;
    ij.clamp(0, (grid_size(level) as i64) - 1) as u32
}

/// Convert a discrete grid coordinate back to the ST coordinate of the cell center.
#[inline]
pub fn ij_to_st(i: u32, level: u8) -> f64 {
    let size = grid_size(level) as f64;
    (i as f64 + 0.5) / size
}

/// Determine which face the point projects onto (face with largest absolute component).
/// Face 0 = +X, 1 = +Y, 2 = +Z, 3 = -X, 4 = -Y, 5 = -Z.
#[inline]
fn largest_abs_component(p: &S2Point) -> u8 {
    let ax = p.x.abs();
    let ay = p.y.abs();
    let az = p.z.abs();

    if ax > ay {
        if ax > az {
            if p.x > 0.0 { 0 } else { 3 }
        } else {
            if p.z > 0.0 { 2 } else { 5 }
        }
    } else {
        if ay > az {
            if p.y > 0.0 { 1 } else { 4 }
        } else {
            if p.z > 0.0 { 2 } else { 5 }
        }
    }
}

/// Given a point and its face, compute UV coordinates.
#[inline]
fn valid_face_xyz_to_uv(face: u8, p: &S2Point) -> (f64, f64) {
    match face {
        0 => (p.y / p.x, p.z / p.x),
        1 => (-p.x / p.y, p.z / p.y),
        2 => (-p.x / p.z, -p.y / p.z),
        3 => (p.z / p.x, p.y / p.x),
        4 => (p.z / p.y, -p.x / p.y),
        5 => (-p.y / p.z, p.x / p.z),
        _ => (0.0, 0.0),
    }
}

/// Convert face + UV to face + ST (applying the quadratic projection).
#[inline]
pub fn face_uv_to_st(u: f64, v: f64) -> (f64, f64) {
    (uv_to_st(u), uv_to_st(v))
}

/// Convert face + ST to face + UV (inverse quadratic projection).
#[inline]
pub fn face_st_to_uv(s: f64, t: f64) -> (f64, f64) {
    (st_to_uv(s), st_to_uv(t))
}

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

    #[test]
    fn st_uv_round_trip() {
        for &u in &[-1.0, -0.5, 0.0, 0.25, 0.5, 0.75, 1.0] {
            let s = uv_to_st(u);
            let u2 = st_to_uv(s);
            assert!(
                (u - u2).abs() < 1e-14,
                "round trip failed for u={u}: got {u2}"
            );
        }
    }

    #[test]
    fn face_detection() {
        // +X face
        let p = S2Point::new(1.0, 0.1, 0.2);
        let (face, _, _) = xyz_to_face_uv(&p);
        assert_eq!(face, 0);

        // +Z face
        let p = S2Point::new(0.1, 0.2, 1.0);
        let (face, _, _) = xyz_to_face_uv(&p);
        assert_eq!(face, 2);

        // -Y face
        let p = S2Point::new(0.1, -1.0, 0.2);
        let (face, _, _) = xyz_to_face_uv(&p);
        assert_eq!(face, 4);
    }

    #[test]
    fn face_uv_xyz_round_trip() {
        let p = S2Point::from_lat_lng(0.3, 0.7).normalize();
        let (face, u, v) = xyz_to_face_uv(&p);
        let p2 = face_uv_to_xyz(face, u, v).normalize();
        assert!((p.x - p2.x).abs() < 1e-14);
        assert!((p.y - p2.y).abs() < 1e-14);
        assert!((p.z - p2.z).abs() < 1e-14);
    }

    #[test]
    fn st_to_ij_boundaries() {
        // At level 1, grid_size = 2, so valid IJ = [0, 1]
        assert_eq!(st_to_ij(0.0, 1), 0);
        assert_eq!(st_to_ij(0.25, 1), 0);
        assert_eq!(st_to_ij(0.5, 1), 1);
        assert_eq!(st_to_ij(0.99, 1), 1);
    }

    #[test]
    fn ij_to_st_center() {
        // At level 2, grid = 4. Cell 0 center = 0.125, cell 1 center = 0.375, etc.
        assert!((ij_to_st(0, 2) - 0.125).abs() < 1e-14);
        assert!((ij_to_st(1, 2) - 0.375).abs() < 1e-14);
        assert!((ij_to_st(2, 2) - 0.625).abs() < 1e-14);
        assert!((ij_to_st(3, 2) - 0.875).abs() < 1e-14);
    }

    #[test]
    fn all_six_faces_covered() {
        let points = [
            S2Point::new(1.0, 0.0, 0.0),  // face 0
            S2Point::new(0.0, 1.0, 0.0),  // face 1
            S2Point::new(0.0, 0.0, 1.0),  // face 2
            S2Point::new(-1.0, 0.0, 0.0), // face 3
            S2Point::new(0.0, -1.0, 0.0), // face 4
            S2Point::new(0.0, 0.0, -1.0), // face 5
        ];
        for (expected, p) in points.iter().enumerate() {
            let (face, _, _) = xyz_to_face_uv(p);
            assert_eq!(face, expected as u8, "expected face {expected}");
        }
    }

    #[test]
    fn lat0_lon0_is_face0() {
        let p = S2Point::from_lat_lng(0.0, 0.0);
        let (face, u, v) = xyz_to_face_uv(&p);
        assert_eq!(face, 0);
        assert!(u.abs() < 1e-14);
        assert!(v.abs() < 1e-14);
    }
}