diff --git a/.github/workflows/pages.yml b/.github/workflows/pages.yml index a3ca058..2b9c83d 100644 --- a/.github/workflows/pages.yml +++ b/.github/workflows/pages.yml @@ -1,36 +1,3 @@ -# on: -# push: -# branches: -# - master - -# jobs: -# generate-docs: -# runs-on: ubuntu-latest - -# steps: -# - name: Checkout code -# uses: actions/checkout@v4 - -# - name: Setup Pages -# id: pages -# uses: actions/configure-pages@v4 - -# - name: Set up Bun -# uses: oven-sh/setup-bun@v1 -# with: -# bun-version: 1.1.29 - -# - name: Install dependencies -# run: bun install - -# - name: Build the docs -# run: bun run docs - -# - name: Upload artifacts -# uses: actions/upload-artifact@v4 -# with: -# path: ./docs # Path to the docs folder - # Simple workflow for deploying static content to GitHub Pages name: Deploy static content to Pages diff --git a/rust/convert.rs b/rust/convert.rs deleted file mode 100644 index bc8c6e3..0000000 --- a/rust/convert.rs +++ /dev/null @@ -1,83 +0,0 @@ -// import { toWM } from './s2'; -// import { toS2, toUnitScale, toVector } from './wm'; - -use alloc::vec; -use alloc::vec::Vec; - -use crate::{Feature, JSONCollection, Projection, VectorFeature, WMFeature}; - -/// Given an input data, convert it to a vector of VectorFeature -pub fn convert( - projection: Projection, - data: &JSONCollection, - tolerance: Option, - maxzoom: Option, - build_bbox: Option, -) -> Vec> { - let mut res: Vec> = vec![]; - - match data { - JSONCollection::FeatureCollection(feature_collection) => { - for feature in &feature_collection.features { - match &feature { - WMFeature::Feature(feature) => { - res.extend(convert_feature( - projection, feature, tolerance, maxzoom, build_bbox, - )); - } - WMFeature::VectorFeature(feature) => { - res.extend(convert_vector_feature(projection, feature, tolerance, maxzoom)) - } - } - } - } - JSONCollection::S2FeatureCollection(feature_collection) => { - for feature in &feature_collection.features { - res.extend(convert_vector_feature(projection, feature, tolerance, maxzoom)); - } - } - JSONCollection::Feature(feature) => { - res.extend(convert_feature(projection, feature, tolerance, maxzoom, build_bbox)); - } - JSONCollection::VectorFeature(feature) => { - res.extend(convert_vector_feature(projection, feature, tolerance, maxzoom)); - } - } - - res -} - -/// Convert a GeoJSON Feature to the appropriate VectorFeature -fn convert_feature( - projection: Projection, - data: &Feature, - tolerance: Option, - maxzoom: Option, - build_bbox: Option, -) -> Vec> { - let mut vf: VectorFeature = Feature::::to_vector(data, build_bbox); - match projection { - Projection::S2 => vf.to_s2(tolerance, maxzoom), - Projection::WM => { - vf.to_unit_scale(tolerance, maxzoom); - vec![vf] - } - } -} - -/// Convert a GeoJSON VectorFeature to the appropriate VectorFeature -fn convert_vector_feature( - projection: Projection, - data: &VectorFeature, - tolerance: Option, - maxzoom: Option, -) -> Vec> { - match projection { - Projection::S2 => data.to_s2(tolerance, maxzoom), - Projection::WM => { - let mut vf = data.to_wm(); - vf.to_unit_scale(tolerance, maxzoom); - vec![vf] - } - } -} diff --git a/rust/lib.rs b/rust/lib.rs index 7dc2929..3783620 100644 --- a/rust/lib.rs +++ b/rust/lib.rs @@ -8,31 +8,13 @@ extern crate alloc; -/// Global conversion tool -pub mod convert; /// All geometry types and structs pub mod geometry; -/// All S2 tooling -pub mod s2; -/// All simplify tooling -pub mod simplify; -/// Tile Structure -pub mod tile; -/// All utility tools -pub mod util; /// All values types and structs pub mod values; -/// All WM tooling -pub mod wm; -pub use convert::*; pub use geometry::*; -pub use s2::*; -pub use simplify::*; -pub use tile::*; -pub use util::*; pub use values::*; -pub use wm::*; use serde::{Deserialize, Serialize}; diff --git a/rust/s2/cellid.rs b/rust/s2/cellid.rs deleted file mode 100644 index 8b6d8bf..0000000 --- a/rust/s2/cellid.rs +++ /dev/null @@ -1,750 +0,0 @@ -extern crate alloc; - -use libm::{fmax, fmin}; - -use alloc::{string::String, vec::Vec}; - -use crate::{ - face_si_ti_to_xyz, face_uv_to_xyz, ij_to_st, si_ti_to_st, st_to_ij, xyz_to_face_uv, BBox, - LonLat, S2Point, K_INVERT_MASK, K_MAX_CELL_LEVEL, K_SWAP_MASK, LOOKUP_POS, ST_TO_UV, UV_TO_ST, -}; - -use super::LOOKUP_IJ; - -/// Cell ID works with both S2 and WM with a common interface -pub type CellId = S2CellId; - -// The following lookup tables are used to convert efficiently between an -// (i,j) cell index and the corresponding position along the Hilbert curve. -// "lookupPos" maps 4 bits of "i", 4 bits of "j", and 2 bits representing the -// orientation of the current cell into 8 bits representing the order in which -// that subcell is visited by the Hilbert curve, plus 2 bits indicating the -// new orientation of the Hilbert curve within that subcell. (Cell -// orientations are represented as combination of K_SWAP_MASK and kInvertMask.) -// -// "lookupIJ" is an inverted table used for mapping in the opposite -// direction. -// -// We also experimented with looking up 16 bits at a time (14 bits of position -// plus 2 of orientation) but found that smaller lookup tables gave better -// performance. (2KB fits easily in the primary cache.) - -// Values for these constants are *declared* in the *.h file. Even though -// the declaration specifies a value for the constant, that declaration -// is not a *definition* of storage for the value. Because the values are -// supplied in the declaration, we don't need the values here. Failing to -// define storage causes link errors for any code that tries to take the -// address of one of these values. - -// Although only 60 bits are needed to represent the index of a leaf cell, the -// extra position bit lets us encode each cell as its Hilbert curve position -// at the cell center, which is halfway along the portion of the Hilbert curve -// that fills that cell. - -/// The number of bits used to encode the face of the cell. -pub const K_FACE_BITS: u8 = 3; -/// The number of faces in the S2 cell projection -pub const K_NUM_FACES: u8 = 6; -/// The maximum level in the S2 cell decomposition -pub const K_MAX_LEVEL: u64 = K_MAX_CELL_LEVEL as u64; // Valid levels: 0..K_MAX_LEVEL -/// The number of bits used to encode the position along the Hilbert curve -pub const K_POS_BITS: u64 = 2 * K_MAX_LEVEL + 1; -/// The maximum number of cells in the S2 cell decomposition -pub const K_MAX_SIZE: u32 = 1 << K_MAX_LEVEL; -/// The number of bits used to encode the orientation of the Hilbert curve -pub const K_LOOKUP_BITS: u8 = 4; - -/// This is the offset required to wrap around from the beginning of the -/// Hilbert curve to the end or vice versa; see nextWrap() and prevWrap(). -/// SWIG doesn't understand uint64{}, so use static_cast. -pub const K_WRAP_OFFSET: u64 = (K_NUM_FACES as u64) << K_POS_BITS; - -/// An S2CellId is a 64-bit unsigned integer that uniquely identifies a -/// cell in the S2 cell decomposition. It has the following format: -/// -/// id = [face][face_pos] -/// -/// face: a 3-bit number (range 0..5) encoding the cube face. -/// -/// face_pos: a 61-bit number encoding the position of the center of this -/// cell along the Hilbert curve over this face (see the Wiki -/// pages for details). -/// -/// Sequentially increasing cell ids follow a continuous space-filling curve -/// over the entire sphere. They have the following properties: -/// -/// - The id of a cell at level k consists of a 3-bit face number followed -/// by k bit pairs that recursively select one of the four children of -/// each cell. The next bit is always 1, and all other bits are 0. -/// Therefore, the level of a cell is determined by the position of its -/// lowest-numbered bit that is turned on (for a cell at level k, this -/// position is 2 * (K_MAX_LEVEL - k).) -/// -/// - The id of a parent cell is at the midpoint of the range of ids spanned -/// by its children (or by its descendants at any level). -/// -/// Leaf cells are often used to represent points on the unit sphere, and -/// this class provides methods for converting directly between these two -/// representations. For cells that represent 2D regions rather than -/// discrete point, it is better to use the S2Cell class. -/// -/// This class is intended to be copied by value as desired. It uses -/// the default copy constructor and assignment operator. -#[derive(Debug, Copy, Clone, PartialEq, PartialOrd, Ord, Eq, Hash)] -#[repr(C)] -pub struct S2CellId { - /// the id contains the face, s, and t components - pub id: u64, -} -impl S2CellId { - /// Construct a cell id from the given 64-bit integer. - pub fn new(id: u64) -> Self { - S2CellId { id } - } - - /// Returns an empty cell id. - pub fn none() -> Self { - 0_u64.into() - } - - /// Returns an invalid cell id guaranteed to be larger than any - /// valid cell id. Useful for creating indexes. - pub fn sentinel() -> Self { - (!0_u64).into() - } - - /// Return the cell corresponding to a given S2 cube face. - pub fn from_face(face: u8) -> S2CellId { - (((face as u64) << K_POS_BITS) + lsb_for_level(0)).into() - } - - /// Construct a leaf cell containing the given normalized S2LatLng. - pub fn from_lon_lat(ll: &LonLat) -> S2CellId { - S2CellId::from_s2_point(&ll.to_point()) - } - - /// Construct a leaf cell containing the given point "p". Usually there is - /// exactly one such cell, but for points along the edge of a cell, any - /// adjacent cell may be (deterministically) chosen. This is because - /// S2CellIds are considered to be closed sets. The returned cell will - /// always contain the given point, i.e. - /// - /// S2Cell(S2CellId(p)).contains(p) - /// - /// is always true. The point "p" does not need to be normalized. - /// - /// If instead you want every point to be contained by exactly one S2Cell, - /// you will need to convert the S2CellIds to S2Loops (which implement point - /// containment this way). - pub fn from_s2_point(p: &S2Point) -> Self { - let (face, u, v) = xyz_to_face_uv(p); - let i: u32 = st_to_ij(UV_TO_ST(u)); - let j: u32 = st_to_ij(UV_TO_ST(v)); - Self::from_face_ij(face, i, j, None) - } - - /// Construct a leaf cell given its face and (u,v) coordinates. - pub fn from_face_uv(face: u8, u: f64, v: f64) -> Self { - S2CellId::from_face_st(face, UV_TO_ST(u), UV_TO_ST(v)) - } - - /// Construct a leaf cell given its face and (s,t) coordinates. - pub fn from_face_st(face: u8, s: f64, t: f64) -> Self { - let i: u32 = st_to_ij(s); - let j: u32 = st_to_ij(t); - Self::from_face_ij(face, i, j, None) - } - - /// Return a leaf cell given its cube face (range 0..5) and - /// i- and j-coordinates (see s2coords.h). - pub fn from_face_ij(face: u8, i: u32, j: u32, level: Option) -> S2CellId { - let mut i = i as u64; - let mut j = j as u64; - if let Some(level) = level { - i <<= K_MAX_LEVEL - level as u64; - j <<= K_MAX_LEVEL - level as u64; - } - // Optimization notes: - // - Non-overlapping bit fields can be combined with either "+" or "|". - // Generally "+" seems to produce better code, but not always. - - // Note that this value gets shifted one bit to the left at the end - // of the function. - let mut n: u64 = (face as u64) << (K_POS_BITS - 1); - - // Alternating faces have opposite Hilbert curve orientations; this - // is necessary in order for all faces to have a right-handed - // coordinate system. - let mut bits: u64 = (face & K_SWAP_MASK) as u64; - - // Each iteration maps 4 bits of "i" and "j" into 8 bits of the Hilbert - // curve position. The lookup table transforms a 10-bit key of the form - // "iiiijjjjoo" to a 10-bit value of the form "ppppppppoo", where the - // letters [ijpo] denote bits of "i", "j", Hilbert curve position, and - // Hilbert curve orientation respectively. - let mut k: u8 = 7; - loop { - let mask: u64 = (1 << 4) - 1; - bits += ((i >> (k * 4)) & mask) << (4 + 2); - bits += ((j >> (k * 4)) & mask) << 2; - bits = LOOKUP_POS[bits as usize] as u64; - n |= (bits >> 2) << (k * 2 * 4); - bits &= (K_SWAP_MASK | K_INVERT_MASK) as u64; - if k == 0 { - break; - } - k -= 1; - } - - let id: S2CellId = ((n << 1) + 1).into(); - - if let Some(level) = level { - id.parent(Some(level)) - } else { - id - } - } - - /// Given a distance and optional zoom level, construct a cell ID at that distance - pub fn from_distance(distance: u64, level: Option) -> S2CellId { - let mut level: u64 = level.unwrap_or(K_MAX_LEVEL as u8) as u64; - level = 2 * (K_MAX_LEVEL - level); - ((distance << (level + 1)) + (1 << level)).into() - } - - /// convert an id to a zoom-i-j after setting it to a new parent zoom - pub fn to_zoom_ij(&self, level: Option) -> (u8, u32, u32) { - let (face, i, j, _or) = self.to_face_ij_orientation(level); - (face, i, j) - } - - /// Return the (face, i, j) coordinates for the leaf cell corresponding to - /// this cell id. Since cells are represented by the Hilbert curve position - /// at the center of the cell, the returned (i,j) for non-leaf cells will be - /// a leaf cell adjacent to the cell center. If "orientation" is non-nullptr, - /// also return the Hilbert curve orientation for the current cell. - /// Returns (face, i, j, orientation) - pub fn to_face_ij_orientation(&self, level: Option) -> (u8, u32, u32, u8) { - let mut i: u32 = 0; - let mut j: u32 = 0; - let face: u8 = (self.id >> K_POS_BITS) as u8; - let mut bits: u64 = (face & 1) as u64; - - // Each iteration maps 8 bits of the Hilbert curve position into - // 4 bits of "i" and "j". The lookup table transforms a key of the - // form "ppppppppoo" to a value of the form "iiiijjjjoo", where the - // letters [ijpo] represents bits of "i", "j", the Hilbert curve - // position, and the Hilbert curve orientation respectively. - // - // On the first iteration we need to be careful to clear out the bits - // representing the cube face. - let mut k = 7; - while k >= 0 { - let nbits: u64 = if k == 7 { 2 } else { 4 }; - let kk: u64 = k as u64; - bits += ((self.id >> (kk * 8 + 1)) & ((1 << (2 * nbits)) - 1)) << 2; - bits = LOOKUP_IJ[bits as usize] as u64; - i += (bits as u32 >> K_NUM_FACES as u32) << (kk * 4); - j += (((bits >> 2) & 15) << (kk * 4)) as u32; - bits &= K_FACE_BITS as u64; - k -= 1; - } - - // adjust bits to the orientation - let lsb = self.id & (!self.id + 1); - if (lsb & 1229782938247303424) != 0 { - bits ^= 1; - } - - if let Some(level) = level { - i >>= K_MAX_LEVEL as u32 - level as u32; - j >>= K_MAX_LEVEL as u32 - level as u32; - } - - (face, i, j, bits as u8) - } - - /// Return the (face, si, ti) coordinates of the center of the cell. Note - /// that although (si,ti) coordinates span the range [0,2**31] in general, - /// the cell center coordinates are always in the range [1,2**31-1] and - /// therefore can be represented using a signed 32-bit integer. - /// Returns (face, si, ti) - pub fn get_center_si_ti(&self) -> (u8, u32, u32) { - // First we compute the discrete (i,j) coordinates of a leaf cell contained - // within the given cell. Given that cells are represented by the Hilbert - // curve position corresponding at their center, it turns out that the cell - // returned by to_face_ij_orientation is always one of two leaf cells closest - // to the center of the cell (unless the given cell is a leaf cell itself, - // in which case there is only one possibility). - // - // Given a cell of size s >= 2 (i.e. not a leaf cell), and letting (imin, - // jmin) be the coordinates of its lower left-hand corner, the leaf cell - // returned by to_face_ij_orientation() is either (imin + s/2, jmin + s/2) - // (imin + s/2 - 1, jmin + s/2 - 1). The first case is the one we want. - // We can distinguish these two cases by looking at the low bit of "i" or - // "j". In the second case the low bit is one, unless s == 2 (i.e. the - // level just above leaf cells) in which case the low bit is zero. - // - // In the code below, the expression ((i ^ (int(id_) >> 2)) & 1) is true - // if we are in the second case described above. - let (face, i, j, _or) = self.to_face_ij_orientation(None); - let mut delta: u32 = 0; - if self.is_leaf() { - delta = 1; - } else if ((i ^ (self.id >> 2) as u32) & 1) != 0 { - delta = 2; - } - - // Note that (2 * {i,j} + delta) will never overflow a 32-bit integer. - let psi = 2 * i + delta; - let pti = 2 * j + delta; - - (face, psi, pti) - } - - /// Creates a human readable debug string. Used for << and available for - /// direct usage as well. The format is "f/dd..d" where "f" is a digit in - /// the range [0-5] representing the S2CellId face, and "dd..d" is a string - /// of digits in the range [0-3] representing each child's position with - /// respect to its parent. (Note that the latter string may be empty.) - /// - /// For example "4/" represents S2CellId::from_face(4), and "3/02" represents - /// S2CellId::from_face(3).child(0).child(2). - pub fn display_name(&self) -> String { - let mut out = String::new(); - if !self.is_valid() { - out.push_str("Invalid"); - return out; - } - out.push((self.face() + 48) as char); - out += "/"; - let mut cur_level: u8 = 1; - while cur_level <= self.level() { - out.push((self.child_position(cur_level) + 48) as char); - cur_level += 1; - } - - out - } - - /// Return the child position (0..3) of this cell's ancestor at the given - /// level within its parent. For example, childPosition(1) returns the - /// position of this cell's level-1 ancestor within its top-level face cell. - /// REQUIRES: 1 <= level <= this->level(). - pub fn child_position(&self, level: u8) -> u8 { - if !self.is_valid() || level > self.level() { - unreachable!(); - } - // return @as(u3, @intCast(self.id >> (2 * (K_MAX_LEVEL - @as(u6, level)) + 1) & 3)); - ((self.id >> (2 * (K_MAX_LEVEL - level as u64) + 1)) & 3) as u8 - } - - /// Converts a string in the format returned by ToString() to an S2CellId. - /// Returns S2CellId.None() if the string could not be parsed. - /// - /// The method name includes "Debug" in order to avoid possible confusion - /// with FromToken() above. - pub fn from_string(val: &str) -> S2CellId { - // This function is reasonably efficient, but is only intended for use in tests. - let val_bytes = val.as_bytes(); - let level = (val.len() - 2) as u64; - if !(0..=K_MAX_LEVEL).contains(&level) { - return S2CellId::none(); - } - let face: u8 = val_bytes[0] - b'0'; - if !(0..=5).contains(&face) || val_bytes[1] != b'/' { - return S2CellId::none(); - } - let mut id = S2CellId::from_face(face); - let mut i = 2; - while i < val.len() { - let child_pos = val_bytes[i] - b'0'; - if !(0..=3).contains(&child_pos) { - return S2CellId::none(); - } - id = id.child(child_pos); - i += 1; - } - - id - } - - /// Return the lowest-numbered bit that is on for this cell id, which is - /// equal to (uint64{1} << (2 * (K_MAX_LEVEL - level))). So for example, - /// a.lsb() <= b.lsb() if and only if a.level() >= b.level(), but the - /// first test is more efficient. - pub fn lsb(&self) -> u64 { - if self.id == 0 { - return 0; - } - self.id & (!self.id + 1_u64) - } - - /// Return the immediate child of this cell at the given traversal order - /// position (in the range 0 to 3). This cell must not be a leaf cell. - pub fn child(&self, position: u8) -> S2CellId { - if !self.is_valid() || self.is_leaf() || position > 3 { - unreachable!(); - } - // To change the level, we need to move the least-significant bit two - // positions downward. We do this by subtracting (4 * new_lsb) and adding - // new_lsb. Then to advance to the given child cell, we add - // (2 * position * new_lsb). - let new_lsb = self.lsb() >> 2; - (self.id - (3 * new_lsb) + (2 * (position as u64) * new_lsb)).into() - } - - /// Which cube face this cell belongs to, in the range 0..5. - pub fn face(&self) -> u8 { - (self.id >> K_POS_BITS) as u8 - } - - /// The position of the cell center along the Hilbert curve over this face, - /// in the range 0..(2**K_POS_BITS-1). - pub fn pos(&self) -> u64 { - self.id & (0u64 >> K_FACE_BITS) - } - - /// Return the subdivision level of the cell (range 0..K_MAX_LEVEL). - pub fn level(&self) -> u8 { - // We can't just S2_DCHECK(isValid()) because we want level() to be - // defined for end-iterators, i.e. S2CellId::End(kLevel). However there is - // no good way to define S2CellId::None().level(), so we do prohibit that. - if self.id == 0 { - unreachable!(); - } - - K_MAX_LEVEL as u8 - (self.id.trailing_zeros() as u8 >> 1u8) - } - - /// Return true if id represents a valid cell. - /// - /// All methods require isValid() to be true unless otherwise specified - /// (although not all methods enforce this). - pub fn is_valid(&self) -> bool { - self.face() < K_NUM_FACES && (self.lsb() & 0x1555555555555555) != 0_u64 - } - - /// Return true if this is a leaf cell (more efficient than checking - /// whether level() == K_MAX_LEVEL). - pub fn is_leaf(&self) -> bool { - (self.id & 1u64) != 0 - } - - /// Convert an S2CellID to an S2Point - pub fn to_point_raw(&self) -> S2Point { - let (face, si, ti) = self.get_center_si_ti(); - face_si_ti_to_xyz(face, si, ti) - } - - /// Convert an S2CellID to an S2Point in normalized vector coordinates - pub fn to_point(&self) -> S2Point { - let mut p = self.to_point_raw(); - p.normalize(); - p - } - - /// Convert an S2CellID to an Face-S-T coordinate - /// Returns (face, s, t) - pub fn to_st(&self) -> (u8, f64, f64) { - let (face, i, j, _or) = self.to_face_ij_orientation(None); - let s = ij_to_st(i); - let t = ij_to_st(j); - - (face, s, t) - } - - /// Convert an S2CellID to an Face-U-V coordinate - /// Returns (face, u, v) - pub fn to_uv(&self) -> (u8, f64, f64) { - let (face, s, t) = self.to_st(); - let u = ST_TO_UV(s); - let v = ST_TO_UV(t); - (face, u, v) - } - - /// Given an S2CellID, check if it is a Face Cell. - pub fn is_face(&self) -> bool { - (self.id & ((1 << 60) - 1)) == 0 - } - - /// Given an S2CellID, get the distance it spans (or length it covers) - pub fn distance(&self, lev: Option) -> f64 { - let level = lev.unwrap_or(self.level()); - (self.id >> (2 * (30 - level) + 1)) as f64 - } - - /// Given an S2CellID, get all the quad children tiles - pub fn children(&self, orientation: Option) -> [S2CellId; 4] { - let mut childs = [self.child(0), self.child(3), self.child(2), self.child(1)]; - - if let Some(orientation) = orientation { - if orientation == 0 { - childs.swap(1, 3); - } - } - - childs - } - - /// Given an S2CellID, get the quad children tiles using a face-zoom-ij input - pub fn children_ij(&self, face: u8, level: u8, i: u32, j: u32) -> [Self; 4] { - let i = i << 1; - let j = j << 1; - let level = level + 1; - - [ - S2CellId::from_face_ij(face, i, j, Some(level)), - S2CellId::from_face_ij(face, i + 1, j, Some(level)), - S2CellId::from_face_ij(face, i, j + 1, Some(level)), - S2CellId::from_face_ij(face, i + 1, j + 1, Some(level)), - ] - } - - /// Given an S2CellID, get the parent quad tile - pub fn parent(&self, level: Option) -> S2CellId { - let new_lsb = match level { - Some(level) => 1 << (2 * (K_MAX_LEVEL - level as u64)), - None => (self.id & (!self.id + 1)) << 2, - }; - - ((self.id & (!new_lsb + 1)) | new_lsb).into() - } - - /// Given an S2CellID, get the hilbert range it spans - /// returns (min, max) - pub fn range(&self) -> (Self, Self) { - let id = self.id; - let lsb = id & (!id + 1); - - ((id - (lsb - 1)).into(), (id + (lsb - 1)).into()) - } - - /// Check if the first S2CellID contains the second. - pub fn contains(&self, other: &S2CellId) -> bool { - let (min, max) = self.range(); - *other >= min && *other <= max - } - - /// Check if an S2CellID intersects another. This includes edges touching. - pub fn intersects(&self, other: &S2CellId) -> bool { - let (min_self, max_self) = self.range(); - let (min_other, max_other) = other.range(); - min_other <= max_self && max_other >= min_self - } - - /// Get the next S2CellID in the hilbert space - pub fn next(&self) -> S2CellId { - let id = self.id; - let n = id + ((id & (!id + 1)) << 1); - if n < K_WRAP_OFFSET { - n.into() - } else { - (n - K_WRAP_OFFSET).into() - } - } - - /// Get the previous S2CellID in the hilbert space - pub fn prev(&self) -> S2CellId { - let id = self.id; - let p = id - ((id & (!id + 1)) << 1); - if p < K_WRAP_OFFSET { - p.into() - } else { - (p + K_WRAP_OFFSET).into() - } - } - - /// Given an S2CellID and level (zoom), get the center point of that cell in S-T space - /// returns (face, s, t) - pub fn center_st(&self) -> (u8, f64, f64) { - let id = self.id; - let (face, i, j, _or) = self.to_face_ij_orientation(None); - let delta = if (id & 1) != 0 { - 1 - } else if ((i as u64 ^ (id >> 2)) & 1) != 0 { - 2 - } else { - 0 - }; - // Note that (2 * {i,j} + delta) will never overflow a 32-bit integer. - let si = 2 * i + delta; - let ti = 2 * j + delta; - - (face, si_ti_to_st(si), si_ti_to_st(ti)) - } - - /// Given an S2CellID and level (zoom), get the S-T bounding range of that cell - /// returns (s_min, t_min, s_max, t_max) - pub fn bounds_st(&self, level: Option) -> BBox { - let level = level.unwrap_or(self.level()); - - let (_face, s, t) = self.center_st(); - let half_size = size_st(level) * 0.5; - BBox::new(s - half_size, t - half_size, s + half_size, t + half_size) - } - - /// Given an S2CellID, find the neighboring S2CellIDs - /// returns neighbors: [up, right, down, left] - pub fn neighbors(&self) -> [S2CellId; 4] { - let level = self.level(); - let size = size_ij(level); - let (face, i, j, _or) = self.to_face_ij_orientation(None); - - [ - S2CellId::from_ij_same(face, i, j - size, j as i32 - size as i32 >= 0) - .parent(Some(level)), - S2CellId::from_ij_same(face, i + size, j, i + size < K_MAX_SIZE).parent(Some(level)), - S2CellId::from_ij_same(face, i, j + size, j + size < K_MAX_SIZE).parent(Some(level)), - S2CellId::from_ij_same(face, i - size, j, i as i32 - size as i32 >= 0) - .parent(Some(level)), - ] - } - - /// Given a Face-I-J and a desired level (zoom), find the neighboring S2CellIDs - /// return neighbors: [down, right, up, left] - pub fn neighbors_ij(&self, face: u8, i: u32, j: u32, level: u8) -> [S2CellId; 4] { - let size = size_ij(level); - - [ - S2CellId::from_ij_same(face, i, j - size, j as i32 - size as i32 >= 0) - .parent(Some(level)), - S2CellId::from_ij_same(face, i + size, j, i + size < K_MAX_SIZE).parent(Some(level)), - S2CellId::from_ij_same(face, i, j + size, j + size < K_MAX_SIZE).parent(Some(level)), - S2CellId::from_ij_same(face, i - size, j, i as i32 - size as i32 >= 0) - .parent(Some(level)), - ] - } - - /// Build an S2CellID given a Face-I-J, but ensure the face is the same if desired - pub fn from_ij_same(face: u8, i: u32, j: u32, same_face: bool) -> S2CellId { - if same_face { - S2CellId::from_face_ij(face, i, j, None) - } else { - S2CellId::from_face_ij_wrap(face, i, j) - } - } - - /// Build an S2CellID given a Face-I-J, but ensure it's a legal value, otherwise wrap before creation - pub fn from_face_ij_wrap(face: u8, i: u32, j: u32) -> S2CellId { - // Convert i and j to the coordinates of a leaf cell just beyond the - // boundary of this face. This prevents 32-bit overflow in the case - // of finding the neighbors of a face cell. - let i = i.clamp(0, K_MAX_SIZE); - let j = j.clamp(0, K_MAX_SIZE); - - // We want to wrap these coordinates onto the appropriate adjacent face. - // The easiest way to do this is to convert the (i,j) coordinates to (x,y,z) - // (which yields a point outside the normal face boundary), and then call - // S2::XYZtoFaceUV() to project back onto the correct face. - // - // The code below converts (i,j) to (si,ti), and then (si,ti) to (u,v) using - // the linear projection (u=2*s-1 and v=2*t-1). (The code further below - // converts back using the inverse projection, s=0.5*(u+1) and t=0.5*(v+1). - // Any projection would work here, so we use the simplest.) We also clamp - // the (u,v) coordinates so that the point is barely outside the - // [-1,1]x[-1,1] face rectangle, since otherwise the reprojection step - // (which divides by the new z coordinate) might change the other - // coordinates enough so that we end up in the wrong leaf cell. - let k_scale = 1. / K_MAX_SIZE as f64; - let k_limit = 1. + 2.220_446_049_250_313e-16; - let u = fmax( - -k_limit, - fmin(k_limit, k_scale * (2. * (i as f64 - (K_MAX_SIZE as f64) / 2.) + 1.)), - ); - let v = fmax( - -k_limit, - fmin(k_limit, k_scale * (2. * (j as f64 - (K_MAX_SIZE as f64) / 2.) + 1.)), - ); - - // Find the leaf cell coordinates on the adjacent face, and convert - // them to a cell id at the appropriate level. - let (n_face, nu, nv) = xyz_to_face_uv(&face_uv_to_xyz(face, u, v)); - S2CellId::from_face_ij(n_face, st_to_ij(0.5 * (nu + 1.)), st_to_ij(0.5 * (nv + 1.)), None) - } - - /// Given an S2CellID, find it's nearest neighbors associated with it - pub fn vertex_neighbors(&self, level: Option) -> Vec { - let level = level.unwrap_or(self.level()); - let mut neighbors: Vec = Vec::new(); - - let (face, i, j, _or) = self.to_face_ij_orientation(None); - - // Determine the i- and j-offsets to the closest neighboring cell in each - // direction. This involves looking at the next bit of "i" and "j" to - // determine which quadrant of this->parent(level) this cell lies in. - let halfsize = size_ij(level + 1); - let size = halfsize << 1; - let isame; - let jsame; - let ioffset; - let joffset; - - if (i & halfsize) != 0 { - ioffset = size; - isame = i + size < K_MAX_SIZE; - } else { - ioffset = -(size as i32) as u32; - isame = i as i32 - size as i32 >= 0; - } - if (j & halfsize) != 0 { - joffset = size; - jsame = j + size < K_MAX_SIZE; - } else { - joffset = -(size as i32) as u32; - jsame = j as i32 - size as i32 >= 0; - } - - neighbors.push(self.parent(Some(level))); - neighbors.push(S2CellId::from_ij_same(face, i + ioffset, j, isame).parent(Some(level))); - neighbors.push(S2CellId::from_ij_same(face, i, j + joffset, jsame).parent(Some(level))); - if isame || jsame { - neighbors.push( - S2CellId::from_ij_same(face, i + ioffset, j + joffset, isame && jsame) - .parent(Some(level)), - ); - } - - neighbors - } -} -impl From for S2CellId { - fn from(value: u64) -> Self { - S2CellId::new(value) - } -} -impl From for S2CellId { - fn from(value: LonLat) -> Self { - S2CellId::from_lon_lat(&value) - } -} -impl From for S2CellId { - fn from(value: S2Point) -> Self { - S2CellId::from_s2_point(&value) - } -} -impl From for S2CellId { - fn from(s: String) -> Self { - S2CellId::from_string(&s) - } -} -impl From<&str> for S2CellId { - fn from(s: &str) -> Self { - S2CellId::from_string(s) - } -} - -/// Return the lowest-numbered bit that is on for cells at the given level. -pub fn lsb_for_level(level: u8) -> u64 { - 1_u64 << (2 * (K_MAX_LEVEL - (level as u64))) -} - -/// Return the range maximum of a level (zoom) in S-T space -pub fn size_st(level: u8) -> f64 { - ij_to_st(size_ij(level)) -} - -/// Return the range maximum of a level (zoom) in I-J space -pub fn size_ij(level: u8) -> u32 { - 1 << (30 - level) -} diff --git a/rust/s2/convert.rs b/rust/s2/convert.rs deleted file mode 100644 index 4b53b14..0000000 --- a/rust/s2/convert.rs +++ /dev/null @@ -1,46 +0,0 @@ -use crate::{Face, LonLat, S2CellId, VectorFeature, VectorGeometry, VectorPoint}; - -impl VectorFeature { - /// Convert an S2 Feature to a GeoJSON Vector Feature - pub fn to_wm(&self) -> Self { - if self._type == "VectorFeature" { - return self.clone(); - } - let mut geometry = self.geometry.clone(); - convert_geometry(self.face, &mut geometry); - VectorFeature::::new_wm( - self.id, - self.properties.clone(), - geometry, - self.metadata.clone(), - ) - } -} - -/// Underlying conversion mechanic to move S2Geometry to GeoJSON Geometry -fn convert_geometry(face: Face, geometry: &mut VectorGeometry) { - match geometry { - VectorGeometry::Point(point) => convert_geometry_point(face, &mut point.coordinates), - VectorGeometry::LineString(points) | VectorGeometry::MultiPoint(points) => { - points.coordinates.iter_mut().for_each(|point| convert_geometry_point(face, point)) - } - VectorGeometry::Polygon(lines) | VectorGeometry::MultiLineString(lines) => lines - .coordinates - .iter_mut() - .for_each(|line| line.iter_mut().for_each(|point| convert_geometry_point(face, point))), - VectorGeometry::MultiPolygon(polygons) => { - polygons.coordinates.iter_mut().for_each(|polygon| { - polygon.iter_mut().for_each(|line| { - line.iter_mut().for_each(|point| convert_geometry_point(face, point)) - }) - }) - } - } -} - -/// Mutate an S2 Point to a GeoJSON Point -fn convert_geometry_point(face: Face, point: &mut VectorPoint) { - let LonLat { lon, lat } = (&S2CellId::from_face_st(face.into(), point.x, point.y)).into(); - point.x = lon; - point.y = lat; -} diff --git a/rust/s2/coords.rs b/rust/s2/coords.rs deleted file mode 100644 index b27946d..0000000 --- a/rust/s2/coords.rs +++ /dev/null @@ -1,517 +0,0 @@ -use core::f64::consts::PI; - -use libm::{atan, floor, round, sqrt, tan}; - -use crate::{ - s2::{S2Point, K_FACE_UVW_AXES, K_FACE_UVW_FACES}, - Point, -}; - -// This file contains documentation of the various coordinate systems used -// throughout the library. Most importantly, S2 defines a framework for -// decomposing the unit sphere into a hierarchy of "cells". Each cell is a -// quadrilateral bounded by four geodesics. The top level of the hierarchy is -// obtained by projecting the six faces of a cube onto the unit sphere, and -// lower levels are obtained by subdividing each cell into four children -// recursively. Cells are numbered such that sequentially increasing cells -// follow a continuous space-filling curve over the entire sphere. The -// transformation is designed to make the cells at each level fairly uniform -// in size. -// -// -////////////////////////// S2Cell Decomposition ///////////////////////// -// -// The following methods define the cube-to-sphere projection used by -// the S2Cell decomposition. -// -// In the process of converting a latitude-longitude pair to a 64-bit cell -// id, the following coordinate systems are used: -// -// (id) -// An S2CellId is a 64-bit encoding of a face and a Hilbert curve position -// on that face. The Hilbert curve position implicitly encodes both the -// position of a cell and its subdivision level (see s2cell_id.h). -// -// (face, i, j) -// Leaf-cell coordinates. "i" and "j" are integers in the range -// [0,(2**30)-1] that identify a particular leaf cell on the given face. -// The (i, j) coordinate system is right-handed on each face, and the -// faces are oriented such that Hilbert curves connect continuously from -// one face to the next. -// -// (face, s, t) -// Cell-space coordinates. "s" and "t" are real numbers in the range -// [0,1] that identify a point on the given face. For example, the point -// (s, t) = (0.5, 0.5) corresponds to the center of the top-level face -// cell. This point is also a vertex of exactly four cells at each -// subdivision level greater than zero. -// -// (face, si, ti) -// Discrete cell-space coordinates. These are obtained by multiplying -// "s" and "t" by 2**31 and rounding to the nearest u32eger. -// Discrete coordinates lie in the range [0,2**31]. This coordinate -// system can represent the edge and center positions of all cells with -// no loss of precision (including non-leaf cells). In binary, each -// coordinate of a level-k cell center ends with a 1 followed by -// (30 - k) 0s. The coordinates of its edges end with (at least) -// (31 - k) 0s. -// -// (face, u, v) -// Cube-space coordinates in the range [-1,1]. To make the cells at each -// level more uniform in size after they are projected onto the sphere, -// we apply a nonlinear transformation of the form u=f(s), v=f(t). -// The (u, v) coordinates after this transformation give the actual -// coordinates on the cube face (modulo some 90 degree rotations) before -// it is projected onto the unit sphere. -// -// (face, u, v, w) -// Per-face coordinate frame. This is an extension of the (face, u, v) -// cube-space coordinates that adds a third axis "w" in the direction of -// the face normal. It is always a right-handed 3D coordinate system. -// Cube-space coordinates can be converted to this frame by setting w=1, -// while (u,v,w) coordinates can be projected onto the cube face by -// dividing by w, i.e. (face, u/w, v/w). -// -// (x, y, z) -// Direction vector (S2Point). Direction vectors are not necessarily unit -// length, and are often chosen to be points on the biunit cube -// [-1,+1]x[-1,+1]x[-1,+1]. They can be be normalized to obtain the -// corresponding point on the unit sphere. -// -// (lon, lat) -// Latitude and longitude (lonlat). Latitudes must be between -90 and -// 90 degrees inclusive, and longitudes must be between -180 and 180 -// degrees inclusive. -// -// Note that the (i, j), (s, t), (si, ti), and (u, v) coordinate systems are -// right-handed on all six faces. - -// S2 is a namespace for constants and simple utility functions that are used -// throughout the S2 library. The name "S2" is derived from the mathematical -// symbol for the two-dimensional unit sphere (note that the "2" refers to the -// dimension of the surface, not the space it is embedded in). - -/// This is the number of levels needed to specify a leaf cell. This -/// constant is defined here so that the S2::Metric class and the conversion -/// functions below can be implemented without including s2cell_id.h. Please -/// see s2cell_id.h for other useful constants and conversion functions. -pub const K_MAX_CELL_LEVEL: u8 = 30; - -/// The maximum index of a valid leaf cell plus one. The range of valid leaf -/// cell indices is [0..kLimitIJ-1]. -pub const K_LIMIT_IJ: u32 = 1 << K_MAX_CELL_LEVEL; // == S2CellId::kMaxSize - -/// The maximum value of an si- or ti-coordinate. The range of valid (si,ti) -/// values is [0..kMaxSiTi]. -pub const K_MAX_SI_TI: u32 = 1 << (K_MAX_CELL_LEVEL + 1); - -/// We have implemented three different projections from cell-space (s,t) to -/// cube-space (u,v): linear, quadratic, and tangent. They have the following -/// tradeoffs: -/// -/// Linear - This is the fastest transformation, but also produces the least -/// uniform cell sizes. Cell areas vary by a factor of about 5.2, with the -/// largest cells at the center of each face and the smallest cells in -/// the corners. -/// -/// Tangent - Transforming the coordinates via atan() makes the cell sizes -/// more uniform. The areas vary by a maximum ratio of 1.4 as opposed to a -/// maximum ratio of 5.2. However, each call to atan() is about as expensive -/// as all of the other calculations combined when converting from points to -/// cell ids, i.e. it reduces performance by a factor of 3. -/// -/// Quadratic - This is an approximation of the tangent projection that -/// is much faster and produces cells that are almost as uniform in size. -/// It is about 3 times faster than the tangent projection for converting -/// cell ids to points or vice versa. Cell areas vary by a maximum ratio of -/// about 2.1. -/// -/// Here is a table comparing the cell uniformity using each projection. "Area -/// ratio" is the maximum ratio over all subdivision levels of the largest cell -/// area to the smallest cell area at that level, "edge ratio" is the maximum -/// ratio of the longest edge of any cell to the shortest edge of any cell at -/// the same level, and "diag ratio" is the ratio of the longest diagonal of -/// any cell to the shortest diagonal of any cell at the same level. "ToPoint" -/// and "FromPoint" are the times in microseconds required to convert cell ids -/// to and from points (unit vectors) respectively. "ToPointRaw" is the time -/// to convert to a non-unit-length vector, which is all that is needed for -/// some purposes. -/// Area Edge Diag ToPointRaw ToPoint FromPoint -/// Ratio Ratio Ratio (microseconds) -/// ------------------------------------------------------------------- -/// Linear: 5.200 2.117 2.959 0.020 0.087 0.085 -/// Tangent: 1.414 1.414 1.704 0.237 0.299 0.258 -/// Quadratic: 2.082 1.802 1.932 0.033 0.096 0.108 -/// -/// The worst-case cell aspect ratios are about the same with all three -/// projections. The maximum ratio of the longest edge to the shortest edge -/// within the same cell is about 1.4 and the maximum ratio of the diagonals -/// within the same cell is about 1.7. -/// -/// NOTE: Currently Tan only has 1e-12 accuracy while Quadratic is within 1e-15. -#[derive(Default)] -pub enum S2Projection { - /// Linear projection - S2LinearProjection, - /// Tangent projection - S2TanProjection, - /// Quadratic projection - #[default] - S2QuadraticProjection, -} - -/// Convert an s- or t-value to the corresponding u- or v-value. This is -/// a non-linear transformation from [0,1] to [-1,1] that attempts to -/// make the cell sizes more uniform. -pub fn st_to_uvlinear(s: f64) -> f64 { - 2. * s - 1. -} -/// Convert an s- or t-value to the corresponding u- or v-value. This is -/// a non-linear transformation from [0,1] to [-1,1] that attempts to -/// make the cell sizes more uniform. -pub fn st_to_uvquadratic(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)) - } -} -/// Convert an s- or t-value to the corresponding u- or v-value. This is -/// a non-linear transformation from [0,1] to [-1,1] that attempts to -/// make the cell sizes more uniform. -pub fn st_to_uvtan(s_: f64) -> f64 { - use core::f64::consts::PI; - // Unfortunately, tan(M_PI_4) is slightly less than 1.0. This isn't due to - // a flaw in the implementation of tan(), it's because the derivative of - // tan(x) at x=pi/4 is 2, and it happens that the two adjacent floating - // point numbers on either side of the infinite-precision value of pi/4 have - // tangents that are slightly below and slightly above 1.0 when rounded to - // the nearest double-precision result. - let s = tan(PI / 2.0 * s_ - PI / 4.0); - s + (1.0 / (1_u64 << 53) as f64) * s -} - -/// The default projection is quadratic -#[cfg(feature = "quadratic")] -pub const ST_TO_UV: fn(f64) -> f64 = st_to_uvquadratic; -/// If settings are updated you can use the tangent projection -#[cfg(all(not(feature = "quadratic"), feature = "tan"))] -pub const ST_TO_UV: fn(f64) -> f64 = st_to_uvtan; -/// If settings are updated you can use the linear projection -#[cfg(all(not(feature = "quadratic"), not(feature = "tan")))] -pub const ST_TO_UV: fn(f64) -> f64 = st_to_uvlinear; - -/// The inverse of the STtoUV transformation. Note that it is not always -/// true that UV_TO_ST(STtoUV(x)) == x due to numerical errors. -pub fn uv_to_stlinear(u: f64) -> f64 { - 0.5 * (u + 1.0) -} -/// The inverse of the STtoUV transformation. Note that it is not always -/// true that UV_TO_ST(STtoUV(x)) == x due to numerical errors. -pub fn uv_to_st_quadratic(u: f64) -> f64 { - if u >= 0. { - 0.5 * sqrt(1.0 + 3.0 * u) - } else { - 1.0 - 0.5 * sqrt(1.0 - 3.0 * u) - } -} -/// The inverse of the STtoUV transformation. Note that it is not always -/// true that UV_TO_ST(STtoUV(x)) == x due to numerical errors. -pub fn uv_to_st_tan(u: f64) -> f64 { - let a: f64 = atan(u); - (2.0 * (1.0 / PI)) * (a + (PI / 4.0)) -} - -/// The default projection is quadratic -#[cfg(feature = "quadratic")] -pub const UV_TO_ST: fn(f64) -> f64 = uv_to_st_quadratic; -/// If settings are updated you can use the tangent projection -#[cfg(all(not(feature = "quadratic"), feature = "tan"))] -pub const UV_TO_ST: fn(f64) -> f64 = uv_to_st_tan; -/// If settings are updated you can use the linear projection -#[cfg(all(not(feature = "quadratic"), not(feature = "tan")))] -pub const UV_TO_ST: fn(f64) -> f64 = uv_to_stlinear; - -/// Convert the i- or j-index of a leaf cell to the minimum corresponding s- -/// or t-value contained by that cell. The argument must be in the range -/// [0..2**30], i.e. up to one position beyond the normal range of valid leaf -/// cell indices. -pub fn ij_to_st(i: u32) -> f64 { - if !(0..=K_LIMIT_IJ).contains(&i) { - unreachable!(); - } - i as f64 / K_LIMIT_IJ as f64 -} - -/// Return the i- or j-index of the leaf cell containing the given -/// s- or t-value. If the argument is outside the range spanned by valid -/// leaf cell indices, return the index of the closest valid leaf cell (i.e., -/// return values are clamped to the range of valid leaf cell indices). -pub fn st_to_ij(s: f64) -> u32 { - (round(K_LIMIT_IJ as f64 * s - 0.5) as u32).clamp(0, K_LIMIT_IJ - 1) -} - -/// Convert an si- or ti-value to the corresponding s- or t-value. -pub fn si_ti_to_st(si: u32) -> f64 { - if si > K_MAX_SI_TI { - unreachable!(); - } - (1.0 / K_LIMIT_IJ as f64) * si as f64 -} - -/// Return the si- or ti-coordinate that is nearest to the given s- or -/// t-value. The result may be outside the range of valid (si,ti)-values. -pub fn st_to_si_ti(s: f64) -> u32 { - round(s * K_MAX_SI_TI as f64) as u32 -} - -/// Convert a direction vector (not necessarily unit length) to an (s,t) point. -pub struct ST { - /// the s coordinate - pub s: f64, - /// the t coordinate - pub t: f64, -} - -/// Convert an S2Point to an (s,t) point. -pub fn to_face_st(p: &S2Point, face: u8) -> ST { - let uv = to_face_uv(p, face); - ST { s: UV_TO_ST(uv.u), t: UV_TO_ST(uv.v) } -} - -/// A U-V coordinate pair. -pub struct UV { - /// the u coordinate - pub u: f64, - /// the v coordinate - pub v: f64, -} - -/// Convert an S2Point to an (u,v) point. -pub fn to_face_uv(p: &S2Point, face: u8) -> UV { - let (valid, u, v) = face_xyz_to_uv(face, p); - debug_assert!(valid, "face_xyz_to_uv failed"); - UV { u, v } -} - -/// Convert (face, u, v) coordinates to a direction vector (not -/// necessarily unit length). -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), - _ => S2Point::new(v, u, -1.0), - } -} - -/// Given a *valid* face for the given point p (meaning that dot product -/// of p with the face normal is positive), return the corresponding -/// u and v values (which may lie outside the range [-1,1]). -/// Returns (pu, pv). -pub fn valid_face_xyz_to_uv(face: u8, p: &S2Point) -> (f64, f64) { - if p.dot(&get_norm(face)) <= 0.0 { - unreachable!(); - } - 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), - _ => (-p.y / p.z, -p.x / p.z), - } -} - -/// Return the face containing the given direction vector. (For points on -/// the boundary between faces, the result is arbitrary but repeatable.) -pub fn get_face(p: &S2Point) -> u8 { - let mut face: u8 = p.largest_abs_component(); - if p.face(face) < 0.0 { - face += 3; - } - face -} - -/// Convert a direction vector (not necessarily unit length) to -/// (face, u, v) coordinates. -/// Returns (face, pu, pv) -pub fn xyz_to_face_uv(p: &S2Point) -> (u8, f64, f64) { - let face: u8 = get_face(p); - let (pu, pv) = valid_face_xyz_to_uv(face, p); - (face, pu, pv) -} - -/// Convert a direction vector (not necessarily unit length) to -/// (face, u, v) coordinates. -/// Returns (face, ps, pt) -pub fn xyz_to_face_st(p: &S2Point) -> (u8, f64, f64) { - let face: u8 = get_face(p); - let (pu, pv) = valid_face_xyz_to_uv(face, p); - (face, UV_TO_ST(pu), UV_TO_ST(pv)) -} - -/// If the dot product of p with the given face normal is positive, -/// set the corresponding u and v values (which may lie outside the range -/// [-1,1]) and return true. Otherwise return false. -pub fn face_xyz_to_uv(face: u8, p: &S2Point) -> (bool, f64, f64) { - if face < 3 { - if p.face(face) <= 0.0 { - return (false, 0.0, 0.0); - } - } else if p.face(face - 3) >= 0.0 { - return (false, 0.0, 0.0); - } - let (pu, pv) = valid_face_xyz_to_uv(face, p); - (true, pu, pv) -} - -/// Transform the given point P to the (u,v,w) coordinate frame of the given -/// face (where the w-axis represents the face normal). -pub fn face_xyz_to_uvw(face: u8, p: &S2Point) -> S2Point { - // The result coordinates are simply the dot products of P with the (u,v,w) - // axes for the given face (see kFaceUVWAxes). - match face { - 0 => S2Point::new(p.y, p.z, p.x), - 1 => S2Point::new(-p.x, p.z, p.y), - 2 => S2Point::new(-p.x, -p.y, p.z), - 3 => S2Point::new(-p.z, -p.y, -p.x), - 4 => S2Point::new(-p.z, p.x, -p.y), - _ => S2Point::new(p.y, p.x, -p.z), - } -} - -/// Return the right-handed normal (not necessarily unit length) for an -/// edge in the direction of the positive v-axis at the given u-value on -/// the given face. (This vector is perpendicular to the plane through -/// the sphere origin that contains the given edge.) -pub fn get_u_norm(face: u8, u: f64) -> S2Point { - match face { - 0 => S2Point::new(u, -1.0, 0.0), - 1 => S2Point::new(1.0, u, 0.0), - 2 => S2Point::new(1.0, 0.0, u), - 3 => S2Point::new(-u, 0.0, 1.0), - 4 => S2Point::new(0.0, -u, 1.0), - _ => S2Point::new(0.0, -1., -u), - } -} - -/// Return the right-handed normal (not necessarily unit length) for an -/// edge in the direction of the positive u-axis at the given v-value on -/// the given face. -pub fn get_v_norm(face: u8, v: f64) -> S2Point { - match face { - 0 => S2Point::new(-v, 0.0, 1.0), - 1 => S2Point::new(0.0, -v, 1.0), - 2 => S2Point::new(0.0, -1.0, -v), - 3 => S2Point::new(v, -1.0, 0.0), - 4 => S2Point::new(1.0, v, 0.0), - _ => S2Point::new(1.0, 0.0, v), - } -} - -/// Return the unit-length normal for the given face. -pub fn get_norm(face: u8) -> S2Point { - get_uvw_axis(face, 2) -} -/// Return the u-axis for the given face. -pub fn get_u_axis(face: u8) -> S2Point { - get_uvw_axis(face, 0) -} -/// Return the v-axis for the given face. -pub fn get_v_axis(face: u8) -> S2Point { - get_uvw_axis(face, 1) -} - -/// Return the given axis of the given face (u=0, v=1, w=2). -pub fn get_uvw_axis(face: u8, axis: usize) -> S2Point { - let p = K_FACE_UVW_AXES[face as usize][axis]; - S2Point::new(p[0], p[1], p[2]) -} - -/// With respect to the (u,v,w) coordinate system of a given face, return the -/// face that lies in the given direction (negative=0, positive=1) of the -/// given axis (u=0, v=1, w=2). For example, GetUVWFace(4, 0, 1) returns the -/// face that is adjacent to face 4 in the positive u-axis direction. -pub fn get_uvw_face(face: u8, axis: usize, direction: usize) -> i32 { - if !(0..=5).contains(&face) { - unreachable!(); - } - if !(0..=2).contains(&axis) { - unreachable!(); - } - if !(0..=1).contains(&direction) { - unreachable!(); - } - K_FACE_UVW_FACES[face as usize][axis][direction] -} - -/// Convert a direction vector (not necessarily unit length) to -/// (face, si, ti) coordinates and, if p is exactly equal to the center of a -/// cell, return the level of this cell (31 otherwise as its outside the bounds of levels). -/// Return (face, zoom, si, ti). -pub fn xyz_to_face_si_ti(p: &S2Point) -> (u8, u8, u32, u32) { - let (face, u, v) = xyz_to_face_uv(p); - let si = st_to_si_ti(UV_TO_ST(u)); - let ti = st_to_si_ti(UV_TO_ST(v)); - // If the levels corresponding to si,ti are not equal, then p is not a cell - // center. The si,ti values 0 and kMaxSiTi need to be handled specially - // because they do not correspond to cell centers at any valid level; they - // are mapped to level -1 by the code below. - // let level = K_MAX_CELL_LEVEL - (si.* | K_MAX_SI_TI); - let level: i32 = K_MAX_CELL_LEVEL as i32 - ((si | K_MAX_SI_TI).trailing_zeros() as i32); - // if (level < 0 or level != K_MAX_CELL_LEVEL - @ctz(ti.* | K_MAX_SI_TI)) return 31; - if level < 0 || level != K_MAX_CELL_LEVEL as i32 - ((ti | K_MAX_SI_TI).trailing_zeros() as i32) - { - return (face, 31, si, ti); - } - // In infinite precision, this test could be changed to ST == SiTi. However, - // due to rounding errors, UV_TO_ST(xyz_to_face_uv(face_uv_to_xyz(STtoUV(...)))) is - // not idempotent. On the other hand, center_raw is computed exactly the same - // way p was originally computed (if it is indeed the center of an S2Cell): - // the comparison can be exact. - // let center: &S2Point = FaceSiTitoXYZ(face, si, ti).Normalize(); - let mut center = face_si_ti_to_xyz(face, si, ti); - center.normalize(); - if *p == center { - (face, level as u8, si, ti) - } else { - (face, 31, si, ti) - } -} - -/// Convert (face, si, ti) coordinates to a direction vector (not necessarily -/// unit length). -pub fn face_si_ti_to_xyz(face: u8, si: u32, ti: u32) -> S2Point { - let u = ST_TO_UV(si_ti_to_st(si)); - let v = ST_TO_UV(si_ti_to_st(ti)); - face_uv_to_xyz(face, u, v) -} - -/// Convert an u-v-zoom coordinate to a tile coordinate -/// returns the tile X-Y coordinate -pub fn tile_xy_from_uv_zoom(u: f64, v: f64, zoom: u8) -> Point { - tile_xy_from_st_zoom(UV_TO_ST(u), UV_TO_ST(v), zoom) -} - -/// Convert an s-t-zoom coordinate to a tile coordinate -/// returns the tile X-Y coordinate -pub fn tile_xy_from_st_zoom(s: f64, t: f64, zoom: u8) -> Point { - let division_factor = (2 / (1 << zoom)) as f64 * 0.5; - - (floor(s / division_factor), floor(t / division_factor)) -} - -// **** TEST FUNCTIONS **** - -/// swap the i and j axes -pub fn swap_axes(ij: usize) -> usize { - ((ij >> 1) & 1) + ((ij & 1) << 1) -} - -/// invert the i and j axes -pub fn invert_bits(ij: usize) -> usize { - ij ^ 3 -} diff --git a/rust/s2/coords_internal.rs b/rust/s2/coords_internal.rs deleted file mode 100644 index c0fd709..0000000 --- a/rust/s2/coords_internal.rs +++ /dev/null @@ -1,246 +0,0 @@ -// The canonical Hilbert traversal order looks like an inverted 'U': -// the subcells are visited in the order (0,0), (0,1), (1,1), (1,0). -// The following tables encode the traversal order for various -// orientations of the Hilbert curve (axes swapped and/or directions -// of the axes reversed). - -/// Together these flags define a cell orientation. If 'kSwapMask' -/// is true, then canonical traversal order is flipped around the -/// diagonal (i.e. i and j are swapped with each other). If -pub const K_SWAP_MASK: u8 = 0x01; -/// 'kInvertMask' is true, then the traversal order is rotated by 180 -/// degrees (i.e. the bits of i and j are inverted, or equivalently, -/// the axis directions are reversed). -pub const K_INVERT_MASK: u8 = 0x02; - -/// kIJtoPos[orientation][ij] -> pos -/// -/// Given a cell orientation and the (i,j)-index of a subcell (0=(0,0), -/// 1=(0,1), 2=(1,0), 3=(1,1)), return the order in which this subcell is -/// visited by the Hilbert curve (a position in the range [0..3]). -pub const K_IJTO_POS: [[usize; 4]; 4] = [ - // (0,0) (0,1) (1,0) (1,1) - [0, 1, 3, 2], // canonical order - [0, 3, 1, 2], // axes swapped - [2, 3, 1, 0], // bits inverted - [2, 1, 3, 0], // swapped & inverted -]; - -/// kPosToIJ[orientation][pos] -> ij -/// -/// Return the (i,j) index of the subcell at the given position 'pos' in the -/// Hilbert curve traversal order with the given orientation. This is the -/// inverse of the previous table: -/// -/// kPosToIJ[r][kIJtoPos[r][ij]] == ij -pub const K_POS_TO_IJ: [[usize; 4]; 4] = [ - // 0 1 2 3 - [0, 1, 3, 2], // canonical order: (0,0), (0,1), (1,1), (1,0) - [0, 2, 3, 1], // axes swapped: (0,0), (1,0), (1,1), (0,1) - [3, 2, 0, 1], // bits inverted: (1,1), (1,0), (0,0), (0,1) - [3, 1, 0, 2], // swapped & inverted: (1,1), (0,1), (0,0), (1,0) -]; - -/// kPosToOrientation[pos] -> orientation_modifier -/// -/// Return a modifier indicating how the orientation of the child subcell -/// with the given traversal position [0..3] is related to the orientation -/// of the parent cell. The modifier should be XOR-ed with the parent -/// orientation to obtain the curve orientation in the child. -pub const K_POS_TO_ORIENTATION: [u8; 4] = [K_SWAP_MASK, 0, 0, K_INVERT_MASK + K_SWAP_MASK]; - -/// The U,V,W axes for each face. -pub const K_FACE_UVW_FACES: [[[i32; 2]; 3]; 6] = [ - // Face 0 - [ - [4, 1], // U - [5, 2], // V - [3, 0], // W - ], - // Face 1 - [ - [0, 3], // U - [5, 2], // V - [4, 1], // W - ], - // Face 2 - [ - [0, 3], // U - [1, 4], // V - [5, 2], // W - ], - // Face 3 - [ - [2, 5], // U - [1, 4], // V - [0, 3], // W - ], - // Face 4 - [ - [2, 5], // U - [3, 0], // V - [1, 4], // W - ], - // Face 5 - [ - [4, 1], // U - [3, 0], // V - [2, 5], // W - ], -]; - -/// The precomputed neighbors of each face (see GetUVWFace). -pub const K_FACE_UVW_AXES: [[[f64; 3]; 3]; 6] = [ - // Face 0 - [ - [0.0, 1.0, 0.0], // U - [0.0, 0.0, 1.0], // V - [1.0, 0.0, 0.0], // W - ], - // Face 1 - [ - [-1.0, 0.0, 0.0], // U - [0.0, 0.0, 1.0], // V - [0.0, 1.0, 0.0], // W - ], - // Face 2 - [ - [-1.0, 0.0, 0.0], // U - [0.0, -1.0, 0.0], // V - [0.0, 0.0, 1.0], // W - ], - // Face 3 - [ - [0.0, 0.0, -1.0], // U - [0.0, -1.0, 0.0], // V - [-1.0, 0.0, 0.0], // W - ], - // Face 4 - [ - [0.0, 0.0, -1.0], // U - [1.0, 0.0, 0.0], // V - [0.0, -1.0, 0.0], // W - ], - // Face 5 - [ - [0.0, 1.0, 0.0], // U - [1.0, 0.0, 0.0], // V - [0.0, 0.0, -1.0], // W - ], -]; - -/// [1 << (2 * kLookupBits + 2)]u16; -pub const LOOKUP_POS: [u16; 1024] = [ - 0, 1, 682, 683, 14, 4, 685, 679, 17, 58, 688, 667, 20, 61, 702, 663, 234, 64, 705, 619, 237, - 78, 708, 615, 240, 81, 762, 603, 254, 84, 765, 599, 257, 938, 768, 427, 260, 941, 782, 423, - 314, 944, 785, 411, 317, 958, 788, 407, 320, 961, 1002, 363, 334, 964, 1005, 359, 337, 1018, - 1008, 347, 340, 1021, 1022, 343, 5, 15, 678, 684, 9, 8, 675, 674, 31, 54, 693, 668, 24, 51, - 697, 658, 230, 69, 719, 620, 227, 73, 712, 610, 245, 95, 758, 604, 249, 88, 755, 594, 271, 934, - 773, 428, 264, 931, 777, 418, 310, 949, 799, 412, 307, 953, 792, 402, 325, 975, 998, 364, 329, - 968, 995, 354, 351, 1014, 1013, 348, 344, 1011, 1017, 338, 59, 16, 666, 689, 55, 30, 669, 692, - 33, 32, 651, 650, 36, 46, 647, 653, 218, 123, 720, 625, 221, 119, 734, 628, 203, 97, 736, 586, - 199, 100, 750, 589, 272, 922, 827, 433, 286, 925, 823, 436, 288, 907, 801, 394, 302, 903, 804, - 397, 379, 976, 986, 369, 375, 990, 989, 372, 353, 992, 971, 330, 356, 1006, 967, 333, 60, 21, - 662, 703, 50, 25, 659, 696, 47, 37, 652, 646, 40, 41, 642, 643, 214, 124, 725, 639, 211, 114, - 729, 632, 204, 111, 741, 582, 194, 104, 745, 579, 277, 918, 828, 447, 281, 915, 818, 440, 293, - 908, 815, 390, 297, 898, 808, 387, 380, 981, 982, 383, 370, 985, 979, 376, 367, 997, 972, 326, - 360, 1001, 962, 323, 65, 235, 618, 704, 68, 231, 621, 718, 122, 219, 624, 721, 125, 215, 638, - 724, 129, 128, 555, 554, 132, 142, 551, 557, 186, 145, 539, 560, 189, 148, 535, 574, 491, 874, - 833, 448, 487, 877, 836, 462, 475, 880, 890, 465, 471, 894, 893, 468, 384, 811, 897, 298, 398, - 807, 900, 301, 401, 795, 954, 304, 404, 791, 957, 318, 79, 236, 614, 709, 72, 226, 611, 713, - 118, 220, 629, 735, 115, 210, 633, 728, 143, 133, 556, 550, 136, 137, 546, 547, 182, 159, 540, - 565, 179, 152, 530, 569, 492, 870, 847, 453, 482, 867, 840, 457, 476, 885, 886, 479, 466, 889, - 883, 472, 389, 812, 911, 294, 393, 802, 904, 291, 415, 796, 950, 309, 408, 786, 947, 313, 80, - 241, 602, 763, 94, 244, 605, 759, 96, 202, 587, 737, 110, 205, 583, 740, 144, 187, 561, 538, - 158, 183, 564, 541, 160, 161, 522, 523, 174, 164, 525, 519, 497, 858, 848, 507, 500, 861, 862, - 503, 458, 843, 864, 481, 461, 839, 878, 484, 443, 817, 912, 282, 439, 820, 926, 285, 417, 778, - 928, 267, 420, 781, 942, 263, 85, 255, 598, 764, 89, 248, 595, 754, 101, 198, 588, 751, 105, - 195, 578, 744, 149, 188, 575, 534, 153, 178, 568, 531, 165, 175, 518, 524, 169, 168, 515, 514, - 511, 854, 853, 508, 504, 851, 857, 498, 454, 844, 869, 495, 451, 834, 873, 488, 444, 831, 917, - 278, 434, 824, 921, 275, 431, 774, 933, 268, 424, 771, 937, 258, 939, 256, 426, 769, 935, 270, - 429, 772, 923, 273, 432, 826, 919, 276, 446, 829, 875, 490, 449, 832, 871, 493, 452, 846, 859, - 496, 506, 849, 855, 510, 509, 852, 513, 512, 171, 170, 516, 526, 167, 173, 570, 529, 155, 176, - 573, 532, 151, 190, 576, 746, 107, 193, 590, 749, 103, 196, 593, 752, 91, 250, 596, 766, 87, - 253, 940, 261, 422, 783, 930, 265, 419, 776, 924, 287, 437, 822, 914, 280, 441, 819, 876, 486, - 463, 837, 866, 483, 456, 841, 860, 501, 502, 863, 850, 505, 499, 856, 527, 517, 172, 166, 520, - 521, 162, 163, 566, 543, 156, 181, 563, 536, 146, 185, 581, 742, 108, 207, 585, 739, 98, 200, - 607, 757, 92, 246, 600, 761, 82, 243, 945, 315, 410, 784, 948, 311, 413, 798, 906, 289, 395, - 800, 909, 292, 391, 814, 881, 474, 464, 891, 884, 477, 478, 887, 842, 459, 480, 865, 845, 455, - 494, 868, 528, 571, 177, 154, 542, 567, 180, 157, 544, 545, 138, 139, 558, 548, 141, 135, 635, - 730, 113, 208, 631, 733, 116, 222, 609, 715, 74, 224, 612, 711, 77, 238, 959, 316, 406, 789, - 952, 306, 403, 793, 902, 303, 396, 805, 899, 296, 386, 809, 895, 470, 469, 892, 888, 467, 473, - 882, 838, 460, 485, 879, 835, 450, 489, 872, 533, 572, 191, 150, 537, 562, 184, 147, 549, 559, - 134, 140, 553, 552, 131, 130, 636, 726, 127, 213, 626, 723, 120, 217, 623, 716, 70, 229, 616, - 706, 67, 233, 960, 321, 362, 1003, 974, 324, 365, 999, 977, 378, 368, 987, 980, 381, 382, 983, - 810, 385, 299, 896, 813, 388, 295, 910, 816, 442, 283, 913, 830, 445, 279, 916, 747, 577, 192, - 106, 743, 580, 206, 109, 731, 634, 209, 112, 727, 637, 212, 126, 640, 641, 42, 43, 654, 644, - 45, 39, 657, 698, 48, 27, 660, 701, 62, 23, 965, 335, 358, 1004, 969, 328, 355, 994, 991, 374, - 373, 988, 984, 371, 377, 978, 806, 399, 300, 901, 803, 392, 290, 905, 821, 438, 284, 927, 825, - 435, 274, 920, 748, 591, 197, 102, 738, 584, 201, 99, 732, 630, 223, 117, 722, 627, 216, 121, - 645, 655, 38, 44, 649, 648, 35, 34, 671, 694, 53, 28, 664, 691, 57, 18, 1019, 336, 346, 1009, - 1015, 350, 349, 1012, 993, 352, 331, 970, 996, 366, 327, 973, 794, 400, 305, 955, 797, 414, - 308, 951, 779, 416, 266, 929, 775, 430, 269, 932, 753, 592, 251, 90, 756, 606, 247, 93, 714, - 608, 225, 75, 717, 622, 228, 71, 699, 656, 26, 49, 695, 670, 29, 52, 673, 672, 11, 10, 676, - 686, 7, 13, 1020, 341, 342, 1023, 1010, 345, 339, 1016, 1007, 357, 332, 966, 1000, 361, 322, - 963, 790, 405, 319, 956, 787, 409, 312, 946, 780, 421, 262, 943, 770, 425, 259, 936, 767, 597, - 252, 86, 760, 601, 242, 83, 710, 613, 239, 76, 707, 617, 232, 66, 700, 661, 22, 63, 690, 665, - 19, 56, 687, 677, 12, 6, 680, 681, 2, 3, -]; - -/// [1 << (2 * kLookupBits + 2)]u16; -pub const LOOKUP_IJ: [u16; 1024] = [ - 0, 1, 1022, 1023, 65, 4, 959, 1018, 69, 68, 955, 954, 6, 67, 1016, 957, 9, 128, 1015, 894, 12, - 193, 1010, 831, 76, 197, 946, 827, 75, 134, 949, 888, 137, 136, 887, 886, 140, 201, 882, 823, - 204, 205, 818, 819, 203, 142, 821, 880, 198, 79, 824, 945, 135, 74, 889, 948, 131, 10, 893, - 1012, 192, 13, 830, 1011, 257, 16, 767, 1006, 260, 81, 762, 943, 324, 85, 698, 939, 323, 22, - 701, 1000, 384, 25, 638, 999, 449, 28, 575, 994, 453, 92, 571, 930, 390, 91, 632, 933, 392, - 153, 630, 871, 457, 156, 567, 866, 461, 220, 563, 802, 398, 219, 624, 805, 335, 214, 689, 808, - 330, 151, 692, 873, 266, 147, 756, 877, 269, 208, 755, 814, 273, 272, 751, 750, 276, 337, 746, - 687, 340, 341, 682, 683, 339, 278, 685, 744, 400, 281, 622, 743, 465, 284, 559, 738, 469, 348, - 555, 674, 406, 347, 616, 677, 408, 409, 614, 615, 473, 412, 551, 610, 477, 476, 547, 546, 414, - 475, 608, 549, 351, 470, 673, 552, 346, 407, 676, 617, 282, 403, 740, 621, 285, 464, 739, 558, - 222, 463, 800, 561, 159, 458, 865, 564, 155, 394, 869, 628, 216, 397, 806, 627, 215, 334, 809, - 688, 210, 271, 812, 753, 146, 267, 876, 757, 149, 328, 875, 694, 87, 326, 937, 696, 82, 263, - 940, 761, 18, 259, 1004, 765, 21, 320, 1003, 702, 24, 385, 998, 639, 89, 388, 935, 634, 93, - 452, 931, 570, 30, 451, 992, 573, 33, 512, 991, 510, 36, 577, 986, 447, 100, 581, 922, 443, 99, - 518, 925, 504, 160, 521, 862, 503, 225, 524, 799, 498, 229, 588, 795, 434, 166, 587, 856, 437, - 168, 649, 854, 375, 233, 652, 791, 370, 237, 716, 787, 306, 174, 715, 848, 309, 111, 710, 913, - 312, 106, 647, 916, 377, 42, 643, 980, 381, 45, 704, 979, 318, 48, 769, 974, 255, 113, 772, - 911, 250, 117, 836, 907, 186, 54, 835, 968, 189, 57, 896, 967, 126, 60, 961, 962, 63, 124, 965, - 898, 59, 123, 902, 901, 120, 185, 904, 839, 118, 188, 969, 834, 55, 252, 973, 770, 51, 251, - 910, 773, 112, 246, 847, 776, 177, 183, 842, 841, 180, 179, 778, 845, 244, 240, 781, 782, 243, - 304, 785, 718, 239, 369, 788, 655, 234, 373, 852, 651, 170, 310, 851, 712, 173, 313, 912, 711, - 110, 316, 977, 706, 47, 380, 981, 642, 43, 379, 918, 645, 104, 441, 920, 583, 102, 444, 985, - 578, 39, 508, 989, 514, 35, 507, 926, 517, 96, 502, 863, 520, 161, 439, 858, 585, 164, 435, - 794, 589, 228, 496, 797, 526, 227, 495, 734, 529, 288, 490, 671, 532, 353, 426, 667, 596, 357, - 429, 728, 595, 294, 366, 727, 656, 297, 303, 722, 721, 300, 299, 658, 725, 364, 360, 661, 662, - 363, 358, 599, 664, 425, 295, 594, 729, 428, 291, 530, 733, 492, 352, 533, 670, 491, 417, 536, - 607, 486, 420, 601, 602, 423, 484, 605, 538, 419, 483, 542, 541, 480, 545, 544, 479, 478, 548, - 609, 474, 415, 612, 613, 410, 411, 611, 550, 413, 472, 672, 553, 350, 471, 737, 556, 287, 466, - 741, 620, 283, 402, 678, 619, 344, 405, 680, 681, 342, 343, 745, 684, 279, 338, 749, 748, 275, - 274, 686, 747, 336, 277, 623, 742, 401, 280, 618, 679, 404, 345, 554, 675, 468, 349, 557, 736, - 467, 286, 560, 801, 462, 223, 625, 804, 399, 218, 629, 868, 395, 154, 566, 867, 456, 157, 569, - 928, 455, 94, 572, 993, 450, 31, 636, 997, 386, 27, 635, 934, 389, 88, 697, 936, 327, 86, 700, - 1001, 322, 23, 764, 1005, 258, 19, 763, 942, 261, 80, 758, 879, 264, 145, 695, 874, 329, 148, - 691, 810, 333, 212, 752, 813, 270, 211, 816, 817, 206, 207, 881, 820, 143, 202, 885, 884, 139, - 138, 822, 883, 200, 141, 825, 944, 199, 78, 828, 1009, 194, 15, 892, 1013, 130, 11, 891, 950, - 133, 72, 953, 952, 71, 70, 956, 1017, 66, 7, 1020, 1021, 2, 3, 1019, 958, 5, 64, 1014, 895, 8, - 129, 951, 890, 73, 132, 947, 826, 77, 196, 1008, 829, 14, 195, 1007, 766, 17, 256, 1002, 703, - 20, 321, 938, 699, 84, 325, 941, 760, 83, 262, 878, 759, 144, 265, 815, 754, 209, 268, 811, - 690, 213, 332, 872, 693, 150, 331, 870, 631, 152, 393, 807, 626, 217, 396, 803, 562, 221, 460, - 864, 565, 158, 459, 929, 568, 95, 454, 932, 633, 90, 391, 996, 637, 26, 387, 995, 574, 29, 448, - 990, 511, 32, 513, 927, 506, 97, 516, 923, 442, 101, 580, 984, 445, 38, 579, 983, 382, 41, 640, - 978, 319, 44, 705, 914, 315, 108, 709, 917, 376, 107, 646, 855, 374, 169, 648, 850, 311, 172, - 713, 786, 307, 236, 717, 789, 368, 235, 654, 792, 433, 230, 591, 857, 436, 167, 586, 861, 500, - 163, 522, 798, 499, 224, 525, 735, 494, 289, 528, 730, 431, 292, 593, 666, 427, 356, 597, 669, - 488, 355, 534, 606, 487, 416, 537, 543, 482, 481, 540, 539, 418, 485, 604, 600, 421, 422, 603, - 598, 359, 424, 665, 535, 354, 489, 668, 531, 290, 493, 732, 592, 293, 430, 731, 657, 296, 367, - 726, 660, 361, 362, 663, 724, 365, 298, 659, 723, 302, 301, 720, 719, 238, 305, 784, 714, 175, - 308, 849, 650, 171, 372, 853, 653, 232, 371, 790, 590, 231, 432, 793, 527, 226, 497, 796, 523, - 162, 501, 860, 584, 165, 438, 859, 582, 103, 440, 921, 519, 98, 505, 924, 515, 34, 509, 988, - 576, 37, 446, 987, 641, 40, 383, 982, 644, 105, 378, 919, 708, 109, 314, 915, 707, 46, 317, - 976, 768, 49, 254, 975, 833, 52, 191, 970, 837, 116, 187, 906, 774, 115, 248, 909, 777, 176, - 247, 846, 780, 241, 242, 783, 844, 245, 178, 779, 843, 182, 181, 840, 905, 184, 119, 838, 908, - 249, 114, 775, 972, 253, 50, 771, 971, 190, 53, 832, 966, 127, 56, 897, 903, 122, 121, 900, - 899, 58, 125, 964, 960, 61, 62, 963, -]; diff --git a/rust/s2/mod.rs b/rust/s2/mod.rs deleted file mode 100644 index 54c11a2..0000000 --- a/rust/s2/mod.rs +++ /dev/null @@ -1,10 +0,0 @@ -mod cellid; -mod convert; -mod coords; -mod coords_internal; -mod point; - -pub use cellid::*; -pub use coords::*; -pub use coords_internal::*; -pub use point::*; diff --git a/rust/s2/point.rs b/rust/s2/point.rs deleted file mode 100644 index b7712c9..0000000 --- a/rust/s2/point.rs +++ /dev/null @@ -1,242 +0,0 @@ -use core::cmp::Ordering; -use core::fmt::Debug; -use core::ops::{Add, Div, Mul, Neg, Sub}; - -use libm::{fabs, sqrt}; - -use crate::{s2::xyz_to_face_uv, wm::LonLat, xyz_to_face_st}; - -use super::S2CellId; - -/// An S2Point represents a point on the unit sphere as a 3D vector. Usually -/// points are normalized to be unit length, but some methods do not require -/// this. See util/math/vector.h for the methods available. Among other -/// things, there are overloaded operators that make it convenient to write -/// arithmetic expressions (e.g. (1-x)*p1 + x*p2). -/// NOTE: asumes only f64 or greater is used. -#[derive(Debug, Copy, Clone, Default, PartialEq)] -pub struct S2Point { - /// The x component. - pub x: f64, - /// The y component. - pub y: f64, - /// The z component. - pub z: f64, -} -impl S2Point { - /// Creates a new S2Point. - pub fn new(x: f64, y: f64, z: f64) -> Self { - S2Point { x, y, z } - } - - /// Returns true if the point is the zero vector. - pub fn is_empty(&self) -> bool { - let zero = f64::default(); - self.x == zero && self.y == zero && self.z == zero - } - - /// Returns the S2 face assocated with this point - pub fn face(&self, f: u8) -> f64 { - if f == 0 { - self.x - } else if f == 1 { - self.y - } else { - self.z - } - } - - /// Returns a Face-ST representation of this point - pub fn to_face_st(&self) -> (u8, f64, f64) { - xyz_to_face_st(self) - } - - /// Returns the S2 face assocated with this point - pub fn get_face(&self) -> u8 { - xyz_to_face_uv(self).0 - } - - /// dot returns the standard dot product of v and ov. - pub fn dot(&self, b: &Self) -> f64 { - self.x * b.x + self.y * b.y + self.z * b.z - } - - /// Returns the absolute value of the point. - pub fn abs(&self) -> Self { - Self::new(fabs(self.x), fabs(self.y), fabs(self.z)) - } - - /// Returns the length of the point. - pub fn len(&self) -> f64 { - sqrt(self.x * self.x + self.y * self.y + self.z * self.z) - } - - /// return the distance from this point to the other point - pub fn distance(&self, b: &Self) -> f64 { - let d = *self - *b; - d.len() - } - - /// Returns the largest absolute component of the point. - pub fn largest_abs_component(&self) -> u8 { - let tmp = self.abs(); - if tmp.x > tmp.y { - if tmp.x > tmp.z { - 0 - } else { - 2 - } - } else if tmp.y > tmp.z { - 1 - } else { - 2 - } - } - - /// Normalize this point to unit length. - pub fn normalize(&mut self) { - let mut len = self.x * self.x + self.y * self.y + self.z * self.z; - if len > 0.0 { - len = sqrt(len); - self.x /= len; - self.y /= len; - self.z /= len; - } - } -} -impl From<&LonLat> for S2Point { - fn from(lonlat: &LonLat) -> Self { - lonlat.to_point() - } -} -impl From<&S2CellId> for S2Point { - fn from(cellid: &S2CellId) -> Self { - cellid.to_point() - } -} -// Implementing the Add trait for S2Point -impl Add for S2Point { - type Output = Self; - fn add(self, other: Self) -> Self::Output { - S2Point { x: self.x + other.x, y: self.y + other.y, z: self.z + other.z } - } -} -impl Add for S2Point { - type Output = Self; - fn add(self, other: f64) -> Self::Output { - S2Point { x: self.x + other, y: self.y + other, z: self.z + other } - } -} -// Implementing the Sub trait for S2Point -impl Sub for S2Point { - type Output = Self; - fn sub(self, other: Self) -> Self::Output { - S2Point { x: self.x - other.x, y: self.y - other.y, z: self.z - other.z } - } -} -impl Sub for S2Point { - type Output = Self; - fn sub(self, other: f64) -> Self::Output { - S2Point { x: self.x - other, y: self.y - other, z: self.z - other } - } -} -// Implementing the Neg trait for S2Point -impl Neg for S2Point { - type Output = Self; - fn neg(self) -> Self::Output { - S2Point { x: -self.x, y: -self.y, z: -self.z } - } -} -// Implementing the Div trait for S2Point -impl Div for S2Point { - type Output = Self; - fn div(self, other: Self) -> Self::Output { - S2Point { x: self.x / other.x, y: self.y / other.y, z: self.z / other.z } - } -} -impl Div for S2Point { - type Output = Self; - fn div(self, other: f64) -> Self::Output { - S2Point { x: self.x / other, y: self.y / other, z: self.z / other } - } -} -// Implementing the Mul trait for S2Point -impl Mul for S2Point { - type Output = Self; - fn mul(self, other: Self) -> Self::Output { - S2Point { x: self.x * other.x, y: self.y * other.y, z: self.z * other.z } - } -} -impl Mul for S2Point { - type Output = Self; - fn mul(self, other: f64) -> Self::Output { - S2Point { x: self.x * other, y: self.y * other, z: self.z * other } - } -} -impl Eq for S2Point {} -impl Ord for S2Point { - fn cmp(&self, other: &Self) -> Ordering { - match self.x.partial_cmp(&other.x) { - Some(Ordering::Equal) => {} - other => return other.unwrap(), // Handle cases where `x` comparison is not equal - } - match self.y.partial_cmp(&other.y) { - Some(Ordering::Equal) => {} - other => return other.unwrap(), // Handle cases where `y` comparison is not equal - } - match self.z.partial_cmp(&other.z) { - Some(order) => order, - None => Ordering::Equal, // This handles the NaN case safely - } - } -} -impl PartialOrd for S2Point { - fn partial_cmp(&self, other: &Self) -> Option { - Some(self.cmp(other)) - } -} - -// /// norm returns the vector's norm. -// pub fn norm(&self) T { -// return @sqrt(self.norm2()); -// } - -// /// norm2 returns the square of the norm. -// pub fn norm2(&self) T { -// return self.dot(self); -// } - -// pub fn normalize(&self) S2PointType(T) { -// const len = self.norm(); -// return Init(self.x / len, self.y / len, self.z / len); -// } - -// pub fn distance(a: *const Self, b: *const S2PointType(T)) T { -// return @sqrt(pow(@abs(b.x - a.x), 2) + pow(@abs(b.y - a.y), 2) + pow(@abs(b.z - a.z), 2)); -// } - -// pub fn cross(&self, b: *const S2PointType(T)) S2PointType(T) { -// return Init( -// self.y * b.z - self.z * b.y, -// self.z * b.x - self.x * b.z, -// self.x * b.y - self.y * b.x, -// ); -// } - -// pub fn intermediate(&self, b: *const S2PointType(T), t: T) S2PointType(T) { -// var c = .{ -// .x = (self.x) + ((b.x - self.x) * (1 - t)), -// .y = (self.y) + ((b.y - self.y) * (1 - t)), -// .z = (self.z) + ((b.z - self.z) * (1 - t)), -// }; -// return c.normalize(); -// } - -// /// Returns the angle between "this" and v in radians, in the range [0, pi]. If -// /// either vector is zero-length, or nearly zero-length, the result will be -// /// zero, regardless of the other value. -// pub fn angle(&self, b: *const S2PointType(T)) T { -// return atan2(T, self.cross(b).norm(), self.dot(b)); -// } -// }; -// } diff --git a/rust/simplify.rs b/rust/simplify.rs deleted file mode 100644 index b9c9624..0000000 --- a/rust/simplify.rs +++ /dev/null @@ -1,194 +0,0 @@ -use crate::{VectorGeometry, VectorLineString, VectorPoint}; - -use libm::pow; - -use alloc::vec; - -impl VectorGeometry { - /// Build sqdistances for a vector geometry - pub fn build_sq_dists(&mut self, tolerance: f64, maxzoom: Option) { - build_sq_dists(self, tolerance, maxzoom); - } - - /// Simplify the geometry to have a tolerance which will be relative to the tile's zoom level. - pub fn simplify(&mut self, tolerance: f64, zoom: u8, maxzoom: Option) { - simplify(self, tolerance, zoom, maxzoom); - } -} - -/// build sqdistances for line vector data -pub fn build_sq_dists(geometry: &mut VectorGeometry, tolerance: f64, maxzoom: Option) { - let maxzoom = maxzoom.unwrap_or(16); - let tol = pow(tolerance / ((1 << maxzoom) as f64 * 4_096.), 2.); - - match geometry { - VectorGeometry::LineString(geo) => { - let coords = &mut geo.coordinates; - build_sq_dist(coords, 0, coords.len() - 1, tol); - } - VectorGeometry::Polygon(geo) | VectorGeometry::MultiLineString(geo) => { - let coords = &mut geo.coordinates; - coords.iter_mut().for_each(|line| build_sq_dist(line, 0, line.len() - 1, tol)); - } - VectorGeometry::MultiPolygon(geo) => { - let coords = &mut geo.coordinates; - coords.iter_mut().for_each(|polygon| { - polygon.iter_mut().for_each(|line| build_sq_dist(line, 0, line.len() - 1, tol)) - }); - } - _ => {} - } -} - -/// calculate simplification of line vector data using -/// optimized Douglas-Peucker algorithm -fn build_sq_dist(coords: &mut VectorLineString, first: usize, last: usize, sq_tolerance: f64) { - coords[first].t = Some(1.); - _build_sq_dist(coords, first, last, sq_tolerance); - coords[last].t = Some(1.); -} - -/// calculate simplification of line vector data using -/// optimized Douglas-Peucker algorithm -fn _build_sq_dist(coords: &mut VectorLineString, first: usize, last: usize, sq_tolerance: f64) { - let mid = (last - first) >> 1; - let mut max_sq_dist = sq_tolerance; - let mut min_pos_to_mid = last - first; - let mut index: Option = None; - - let VectorPoint { x: ax, y: ay, .. } = coords[first]; - let VectorPoint { x: bx, y: by, .. } = coords[last]; - - let mut i = first; - while i < last { - let VectorPoint { x, y, .. } = coords[i]; - let d = get_sq_seg_dist(x, y, ax, ay, bx, by); - - if d > max_sq_dist { - index = Some(i); - max_sq_dist = d; - } else if d == max_sq_dist { - // a workaround to ensure we choose a pivot close to the middle of the list, - // reducing recursion depth, for certain degenerate inputs - let pos_to_mid = isize::abs(i as isize - mid as isize) as usize; - if pos_to_mid < min_pos_to_mid { - index = Some(i); - min_pos_to_mid = pos_to_mid; - } - } - - i += 1; - } - - if max_sq_dist > sq_tolerance { - if let Some(index) = index { - if index - first > 1 { - _build_sq_dist(coords, first, index, sq_tolerance); - } - coords[index].t = Some(max_sq_dist); - if last - index > 1 { - _build_sq_dist(coords, index, last, sq_tolerance); - } - } - } -} - -/// square distance from a point to a segment -fn get_sq_seg_dist(ps: f64, pt: f64, s: f64, t: f64, bs: f64, bt: f64) -> f64 { - let mut s = s; - let mut t = t; - let mut ds = bs - s; - let mut dt = bt - t; - - if ds != 0. || dt != 0. { - let m = ((ps - s) * ds + (pt - t) * dt) / (ds * ds + dt * dt); - - if m > 1. { - s = bs; - t = bt; - } else if m > 0. { - s += ds * m; - t += dt * m; - } - } - - ds = ps - s; - dt = pt - t; - - ds * ds + dt * dt -} - -/// Simplify a vector geometry -pub fn simplify(geometry: &mut VectorGeometry, tolerance: f64, zoom: u8, maxzoom: Option) { - let maxzoom = maxzoom.unwrap_or(16); - let zoom_tol = if zoom >= maxzoom { 0. } else { tolerance / ((1 << zoom) as f64 * 4_096.) }; - match geometry { - VectorGeometry::LineString(geo) => { - geo.coordinates = simplify_line(&geo.coordinates, zoom_tol, false, false); - } - VectorGeometry::Polygon(geo) | VectorGeometry::MultiLineString(geo) => { - geo.coordinates = geo - .coordinates - .iter() - .map(|line| simplify_line(line, zoom_tol, false, false)) - .collect(); - } - VectorGeometry::MultiPolygon(geo) => { - geo.coordinates = geo - .coordinates - .iter() - .map(|polygon| { - polygon.iter().map(|line| simplify_line(line, zoom_tol, true, true)).collect() - }) - .collect(); - } - _ => (), - } -} - -/// simplified a vector line -fn simplify_line( - line: &VectorLineString, - tolerance: f64, - is_polygon: bool, - is_outer: bool, -) -> VectorLineString { - let sq_tolerance = tolerance * tolerance; - let size = line.len() as f64; - if tolerance > 0. && size < (if is_polygon { sq_tolerance } else { tolerance }) { - return line.clone(); - } - - let mut ring: VectorLineString = vec![]; - for point in line { - if tolerance == 0. || (point.t.unwrap_or(0.0)) > sq_tolerance { - ring.push(point.clone()); - } - } - if is_polygon { - rewind(&mut ring, is_outer); - } - - ring -} - -/// In place adjust the ring if necessary -pub fn rewind(ring: &mut VectorLineString, clockwise: bool) { - let len = ring.len(); - let mut area: f64 = 0.; - let mut i = 0; - let mut j = len - 2; - while i < len { - area += (ring[i].x - ring[j].x) * (ring[i].y + ring[j].y); - j = i; - i += 2; - } - i = 0; - if (area > 0.) == clockwise { - let len_half = len / 2; - while i < len_half { - ring.swap(i, len - i - 1); - i += 2; - } - } -} diff --git a/rust/tile.rs b/rust/tile.rs deleted file mode 100644 index 32624d9..0000000 --- a/rust/tile.rs +++ /dev/null @@ -1,287 +0,0 @@ -use alloc::collections::{BTreeMap, BTreeSet}; -use alloc::string::String; -use alloc::string::ToString; -use alloc::vec; -use alloc::vec::Vec; - -use libm::round; - -use crate::{ - convert, CellId, Face, JSONCollection, Projection, TileChildren, VectorFeature, VectorGeometry, - VectorPoint, -}; - -/// If a user creates metadata for a VectorFeature, it needs to define a get_layer function -pub trait HasLayer { - /// Get the layer from metadata if it exists - fn get_layer(&self) -> Option; -} -impl HasLayer for () { - fn get_layer(&self) -> Option { - None - } -} - -/// Tile Class to contain the tile information for splitting or simplifying -pub struct Tile { - /// the tile id - pub id: CellId, - /// the tile's layers - pub layers: BTreeMap>, - /// whether the tile feature geometry has been transformed - pub transformed: bool, -} -impl Tile { - /// Create a new Tile - pub fn new(id: CellId) -> Self { - Self { id, layers: BTreeMap::new(), transformed: false } - } - - /// Returns true if the tile is empty of features - pub fn is_empty(&self) -> bool { - for layer in self.layers.values() { - if !layer.features.is_empty() { - return false; - } - } - - true - } - - /// Add a feature to the tile - pub fn add_feature(&mut self, feature: VectorFeature, layer: Option) { - let layer_name = feature - .metadata - .as_ref() - .and_then(|meta| meta.get_layer()) // Get the layer from metadata if it exists - .or(layer) // Fall back to the provided layer - .unwrap_or_else(|| "default".to_string()); // Fall back to "default" if none found - - if !self.layers.contains_key(&layer_name) { - self.layers.insert(layer_name.clone(), Layer::new(layer_name.clone())); - } - self.layers.get_mut(&layer_name).unwrap().features.push(feature); - } - - /// Simplify the geometry to have a tolerance which will be relative to the tile's zoom level. - /// NOTE: This should be called after the tile has been split into children if that functionality - /// is needed. - pub fn transform(&mut self, tolerance: f64, maxzoom: Option) { - if self.transformed { - return; - } - let (zoom, i, j) = self.id.to_zoom_ij(None); - - for layer in self.layers.values_mut() { - for feature in layer.features.iter_mut() { - feature.geometry.simplify(tolerance, zoom, maxzoom); - feature.geometry.transform(zoom, i as f64, j as f64) - } - } - - self.transformed = true; - } -} - -/// Layer Class to contain the layer information for splitting or simplifying -pub struct Layer { - /// the layer name - pub name: String, - /// the layer's features - pub features: Vec>, -} -impl Layer { - /// Create a new Layer - pub fn new(name: String) -> Self { - Self { name, features: vec![] } - } -} - -/// Options for creating a TileStore -pub struct TileStoreOptions { - /// manually set the projection, otherwise it defaults to whatever the data type is - pub projection: Option, - /// min zoom to generate data on - pub minzoom: Option, - /// max zoom level to cluster the points on - pub maxzoom: Option, - /// tile buffer on each side in pixels - pub index_maxzoom: Option, - /// simplification tolerance (higher means simpler) - pub tolerance: Option, - /// tile buffer on each side so lines and polygons don't get clipped - pub buffer: Option, -} - -/// TileStore Class is a tile-lookup system that splits and simplifies as needed for each tile request */ -pub struct TileStore { - minzoom: u8, // min zoom to preserve detail on - maxzoom: u8, // max zoom to preserve detail on - faces: BTreeSet, // store which faces are active. 0 face could be entire WM projection - index_maxzoom: u8, // max zoom in the tile index - tolerance: f64, // simplification tolerance (higher means simpler) - buffer: f64, // tile buffer for lines and polygons - tiles: BTreeMap>, // stores both WM and S2 tiles - projection: Projection, // projection to build tiles for -} -impl Default for TileStore { - fn default() -> Self { - Self { - minzoom: 0, - maxzoom: 18, - faces: BTreeSet::::new(), - index_maxzoom: 4, - tolerance: 3., - buffer: 0.0625, - tiles: BTreeMap::>::new(), - projection: Projection::S2, - } - } -} -impl TileStore { - /// Create a new TileStore - pub fn new(data: JSONCollection, options: TileStoreOptions) -> Self { - let mut tile_store = Self { - minzoom: options.minzoom.unwrap_or(0), - maxzoom: options.maxzoom.unwrap_or(20), - faces: BTreeSet::::new(), - index_maxzoom: options.index_maxzoom.unwrap_or(4), - tolerance: options.tolerance.unwrap_or(3.), - buffer: options.buffer.unwrap_or(64.), - tiles: BTreeMap::>::new(), - projection: options.projection.unwrap_or(Projection::S2), - }; - // sanity check - debug_assert!( - tile_store.minzoom <= tile_store.maxzoom - && tile_store.maxzoom > 0 - && tile_store.maxzoom <= 20, - "maxzoom should be in the 0-20 range" - ); - // convert features - let features: Vec> = convert( - tile_store.projection, - &data, - Some(tile_store.tolerance), - Some(tile_store.maxzoom), - None, - ); - features.into_iter().for_each(|feature| tile_store.add_feature(feature)); - for i in 0..6 { - tile_store.split_tile(CellId::from_face(i), None, None); - } - - tile_store - } - - /// Add a feature to the tile store - pub fn add_feature(&mut self, feature: VectorFeature) { - let face: u8 = feature.face.into(); - let tile = self.tiles.entry(CellId::from_face(face)).or_insert_with(|| { - self.faces.insert(feature.face); - Tile::new(CellId::from_face(face)) - }); - - tile.add_feature(feature, None); - } - - /// Split tiles given a range - fn split_tile(&mut self, start_id: CellId, end_id: Option, end_zoom: Option) { - let TileStore { buffer, tiles, tolerance, maxzoom, index_maxzoom, .. } = self; - let end_zoom = end_zoom.unwrap_or(*maxzoom); - let mut stack: Vec = vec![start_id]; - // avoid recursion by using a processing queue - while !stack.is_empty() { - // find our next tile to split - let stack_id = stack.pop(); - if stack_id.is_none() { - break; - } - let tile = tiles.get_mut(&stack_id.unwrap()); - // if the tile we need does not exist, is empty, or already transformed, skip it - if tile.is_none() { - continue; - } - let tile = tile.unwrap(); - if tile.is_empty() || tile.transformed { - continue; - } - let tile_zoom = tile.id.level(); - // 1: stop tiling if we reached a defined end - // 2: stop tiling if it's the first-pass tiling, and we reached max zoom for indexing - // 3: stop at currently needed maxzoom OR current tile does not include child - if tile_zoom >= *maxzoom || // 1 - (end_id.is_none() && tile_zoom >= *index_maxzoom) || // 2 - (end_id.is_some() && (tile_zoom > end_zoom || !tile.id.contains(&end_id.unwrap()))) - { - continue; - } - - // split the tile - let TileChildren { - bottom_left: bl_id, - bottom_right: br_id, - top_left: tl_id, - top_right: tr_id, - } = tile.split(Some(*buffer)); - // now that the tile has been split, we can transform it - tile.transform(*tolerance, Some(*maxzoom)); - // push the new features to the stack - stack.extend(vec![bl_id.id, br_id.id, tl_id.id, tr_id.id]); - // store the children - tiles.insert(bl_id.id, bl_id); - tiles.insert(br_id.id, br_id); - tiles.insert(tl_id.id, tl_id); - tiles.insert(tr_id.id, tr_id); - } - } - - /// Get a tile - pub fn get_tile(&mut self, id: CellId) -> Option<&Tile> { - let zoom = id.level(); - let face = id.face(); - // If the zoom is out of bounds, return nothing - if !(0..=20).contains(&zoom) || !self.faces.contains(&face.into()) { - return None; - } - - // we want to find the closest tile to the data. - let mut p_id = id; - while !self.tiles.contains_key(&p_id) && !p_id.is_face() { - p_id = p_id.parent(None); - } - // split as necessary, the algorithm will know if the tile is already split - self.split_tile(p_id, Some(id), Some(zoom)); - - // grab the tile and split it if necessary - self.tiles.get(&id) - } -} - -impl VectorGeometry { - /// Transform the geometry from the 0->1 coordinate system to a tile coordinate system - pub fn transform(&mut self, zoom: u8, ti: f64, tj: f64) { - let zoom = (1 << (zoom as u64)) as f64; - match self { - VectorGeometry::Point(p) => p.coordinates.transform(zoom, ti, tj), - VectorGeometry::LineString(l) | VectorGeometry::MultiPoint(l) => { - l.coordinates.iter_mut().for_each(|p| p.transform(zoom, ti, tj)) - } - VectorGeometry::MultiLineString(l) | VectorGeometry::Polygon(l) => l - .coordinates - .iter_mut() - .for_each(|l| l.iter_mut().for_each(|p| p.transform(zoom, ti, tj))), - VectorGeometry::MultiPolygon(l) => l.coordinates.iter_mut().for_each(|p| { - p.iter_mut().for_each(|l| l.iter_mut().for_each(|p| p.transform(zoom, ti, tj))) - }), - } - } -} - -impl VectorPoint { - /// Transform the point from the 0->1 coordinate system to a tile coordinate system - pub fn transform(&mut self, zoom: f64, ti: f64, tj: f64) { - self.x = round(self.x * zoom - ti); - self.y = round(self.y * zoom - tj); - } -} diff --git a/rust/util.rs b/rust/util.rs deleted file mode 100644 index 3881513..0000000 --- a/rust/util.rs +++ /dev/null @@ -1,22 +0,0 @@ -/// Earth's radius in meters -pub const EARTH_RADIUS: f64 = 6_371_008.8; -/// Earth's equatorial radius in meters -pub const EARTH_RADIUS_EQUATORIAL: f64 = 6_378_137.0; -/// Earth's polar radius in meters -pub const EARTH_RADIUS_POLAR: f64 = 6_356_752.3; -/// The average circumference of the world in meters -pub const EARTH_CIRCUMFERENCE: f64 = 2.0 * core::f64::consts::PI * EARTH_RADIUS; - -/// Mars' radius in meters -pub const MARS_RADIUS: f64 = 3_389_500.0; -/// Mars' equatorial radius in meters -pub const MARS_RADIUS_EQUATORIAL: f64 = 3_396_200.0; -/// Mars' polar radius in meters -pub const MARS_RADIUS_POLAR: f64 = 3_376_200.0; - -/// 900913 (Web Mercator) constant -pub const A: f64 = 6_378_137.0; -/// 900913 (Web Mercator) max extent -pub const MAXEXTENT: f64 = 20_037_508.342789244; -/// 900913 (Web Mercator) maximum latitude -pub const MAXLAT: f64 = 85.0511287798; diff --git a/rust/wm/clip.rs b/rust/wm/clip.rs deleted file mode 100644 index d9e2cef..0000000 --- a/rust/wm/clip.rs +++ /dev/null @@ -1,553 +0,0 @@ -use alloc::string::ToString; -use alloc::vec; -use alloc::vec::Vec; - -use crate::{ - Axis, BBox3D, HasLayer, MValue, Tile, VectorFeature, VectorGeometry, VectorGeometryType, - VectorLineString, VectorLineStringGeometry, VectorMultiLineOffset, VectorMultiLineString, - VectorMultiLineStringGeometry, VectorMultiPointGeometry, VectorMultiPolygon, - VectorMultiPolygonGeometry, VectorMultiPolygonOffset, VectorPoint, VectorPointGeometry, - VectorPolygonGeometry, -}; - -// TODO: Cases of `to_vec` clones large swathes of data. Can we optimize this? - -/// The children of a tile -pub struct TileChildren { - /// The bottom left child tile - pub bottom_left: Tile, - /// The bottom right child tile - pub bottom_right: Tile, - /// The top left child tile - pub top_left: Tile, - /// The top right child tile - pub top_right: Tile, -} - -impl Tile { - /// split tile into 4 children - pub fn split(&mut self, buffer: Option) -> TileChildren { - split_tile(self, buffer) - } -} - -/** - * @param tile - the tile to split - * @param buffer - the buffer around the tile for lines and polygons - * @returns - the tile's children split into 4 sub-tiles - */ -pub fn split_tile(tile: &mut Tile, buffer: Option) -> TileChildren { - let buffer = buffer.unwrap_or(0.0625); - let face = tile.id.face(); - let (zoom, i, j) = tile.id.to_zoom_ij(None); - let [bl_id, br_id, tl_id, tr_id] = tile.id.children_ij(face, zoom, i, j); - let mut children = TileChildren { - bottom_left: Tile::new(bl_id), - bottom_right: Tile::new(br_id), - top_left: Tile::new(tl_id), - top_right: Tile::new(tr_id), - }; - let scale = (1 << zoom) as f64; - let k1 = 0.; - let k2 = 0.5; - let k3 = 0.5; - let k4 = 1.; - let i = i as f64; - let j = j as f64; - - let mut tl: Option>>; - let mut bl: Option>>; - let mut tr: Option>>; - let mut br: Option>>; - - for (name, layer) in tile.layers.iter_mut() { - let features = &layer.features; - let left = _clip(features, scale, i - k1, i + k3, Axis::X, buffer); - let right = _clip(features, scale, i + k2, i + k4, Axis::X, buffer); - - if let Some(left) = left { - bl = _clip(&left, scale, j - k1, j + k3, Axis::Y, buffer); - tl = _clip(&left, scale, j + k2, j + k4, Axis::Y, buffer); - if let Some(bl) = bl { - for d in bl { - children.bottom_left.add_feature(d, Some(name.to_string())); - } - } - if let Some(tl) = tl { - for d in tl { - children.top_left.add_feature(d, Some(name.to_string())); - } - } - } - - if let Some(right) = right { - br = _clip(&right, scale, j - k1, j + k3, Axis::Y, buffer); - tr = _clip(&right, scale, j + k2, j + k4, Axis::Y, buffer); - if let Some(br) = br { - for d in br { - children.bottom_right.add_feature(d, Some(name.to_string())); - } - } - if let Some(tr) = tr { - for d in tr { - children.top_right.add_feature(d, Some(name.to_string())); - } - } - } - } - - children -} - -/// Internal clip function for a collection of VectorFeatures -fn _clip( - features: &[VectorFeature], - scale: f64, - k1: f64, - k2: f64, - axis: Axis, - base_buffer: f64, -) -> Option>> -where - M: Clone, -{ - // scale - let k1 = k1 / scale; - let k2 = k2 / scale; - // prep buffer and result container - let buffer = base_buffer / scale; - let k1b = k1 - buffer; - let k2b = k2 + buffer; - let mut clipped: Vec> = vec![]; - let axis_x = axis == Axis::X; - - for feature in features { - let geometry = &feature.geometry; - // trivial accept and reject cases - if let Some(vec_bbox) = geometry.vec_bbox() { - let min = if axis_x { vec_bbox.left } else { vec_bbox.bottom }; - let max = if axis_x { vec_bbox.right } else { vec_bbox.top }; - if min >= k1 && max < k2 { - clipped.push(feature.clone()); - continue; - } else if max < k1 || min >= k2 { - continue; - } - } - // build the new clipped geometry - let new_geometry: Option = match geometry { - VectorGeometry::Point(geo) => clip_point(geo, axis, k1, k2), - VectorGeometry::MultiPoint(geo) => clip_multi_point(geo, axis, k1, k2), - VectorGeometry::LineString(geo) => clip_line_string(geo, axis, k1b, k2b), - VectorGeometry::MultiLineString(geo) => { - clip_multi_line_string(geo, axis, k1b, k2b, false) - } - VectorGeometry::Polygon(geo) => clip_polygon(geo, axis, k1b, k2b), - VectorGeometry::MultiPolygon(geo) => clip_multi_polygon(geo, axis, k1b, k2b), - }; - // store if the geometry was inside the range - if let Some(new_geometry) = new_geometry { - clipped.push(VectorFeature::from_vector_feature(feature, Some(new_geometry))); - } - } - - if clipped.is_empty() { - None - } else { - Some(clipped) - } -} - -/// Clip a point to an axis and range -fn clip_point( - geometry: &VectorPointGeometry, - axis: Axis, - k1: f64, - k2: f64, -) -> Option { - let coords = &geometry.coordinates; - let value = if axis == Axis::X { coords.x } else { coords.y }; - if value >= k1 && value < k2 { - Some(VectorGeometry::Point(geometry.clone())) - } else { - None - } -} - -/// Clip a MultiPoint to an axis and range -fn clip_multi_point( - geometry: &VectorMultiPointGeometry, - axis: Axis, - k1: f64, - k2: f64, -) -> Option { - let mut new_geo = geometry.clone(); - new_geo.coordinates = geometry - .coordinates - .iter() - .filter(|point| { - let value = if axis == Axis::X { point.x } else { point.y }; - value >= k1 && value < k2 - }) - .cloned() - .collect(); - - if new_geo.coordinates.is_empty() { - None - } else { - Some(VectorGeometry::MultiPoint(new_geo)) - } -} - -/// Clip a LineString to an axis and range -fn clip_line_string( - geometry: &VectorLineStringGeometry, - axis: Axis, - k1: f64, - k2: f64, -) -> Option { - let VectorLineStringGeometry { is_3d, coordinates: line, bbox, vec_bbox, .. } = geometry; - let init_o = geometry.offset.unwrap_or(0.); - let mut new_offsets: VectorMultiLineOffset = vec![]; - let mut new_lines: VectorMultiLineString = vec![]; - for clip in - _clip_line(ClipLineResult { line: line.to_vec(), offset: init_o }, k1, k2, axis, false) - { - new_offsets.push(clip.offset); - new_lines.push(clip.line); - } - if new_lines.is_empty() { - None - } else { - Some(VectorGeometry::MultiLineString(VectorMultiLineStringGeometry { - _type: VectorGeometryType::MultiLineString, - is_3d: *is_3d, - coordinates: new_lines, - bbox: *bbox, - offset: Some(new_offsets), - vec_bbox: Some(vec_bbox.unwrap_or_default().clip(axis, k1, k2)), - ..Default::default() - })) - } -} - -/// Clip a MultiLineString geometry to an axis and range -fn clip_multi_line_string( - geometry: &VectorMultiLineStringGeometry, - axis: Axis, - k1: f64, - k2: f64, - is_polygon: bool, -) -> Option { - let VectorMultiLineStringGeometry { is_3d, coordinates, bbox, vec_bbox, .. } = geometry; - let init_o = - geometry.offset.clone().unwrap_or_else(|| coordinates.iter().map(|_| 0.).collect()); - let mut new_offsets: VectorMultiLineOffset = vec![]; - let mut new_lines: VectorMultiLineString = vec![]; - let vec_bbox = vec_bbox.unwrap_or_default().clip(axis, k1, k2); - coordinates.iter().enumerate().for_each(|(i, line)| { - for clip in _clip_line( - ClipLineResult { line: line.to_vec(), offset: init_o[i] }, - k1, - k2, - axis, - is_polygon, - ) { - new_offsets.push(clip.offset); - new_lines.push(clip.line); - } - }); - if new_lines.is_empty() || (is_polygon && new_lines[0].len() < 4) { - None - } else if !is_polygon { - Some(VectorGeometry::MultiLineString(VectorMultiLineStringGeometry { - _type: VectorGeometryType::MultiLineString, - is_3d: *is_3d, - coordinates: new_lines, - bbox: *bbox, - offset: Some(new_offsets), - vec_bbox: Some(vec_bbox), - ..Default::default() - })) - } else { - Some(VectorGeometry::Polygon(VectorPolygonGeometry { - _type: VectorGeometryType::Polygon, - is_3d: *is_3d, - coordinates: new_lines, - bbox: *bbox, - offset: Some(new_offsets), - vec_bbox: Some(vec_bbox), - ..Default::default() - })) - } -} - -/// Clip a Polygon geometry to an axis and range -fn clip_polygon( - geometry: &VectorPolygonGeometry, - axis: Axis, - k1: f64, - k2: f64, -) -> Option { - clip_multi_line_string(geometry, axis, k1, k2, true) -} - -/// Clip a MultiPolygon geometry to an axis and range -fn clip_multi_polygon( - geometry: &VectorMultiPolygonGeometry, - axis: Axis, - k1: f64, - k2: f64, -) -> Option { - let VectorMultiPolygonGeometry { is_3d, coordinates, bbox, vec_bbox, .. } = geometry; - let init_o = geometry - .offset - .clone() - .unwrap_or_else(|| coordinates.iter().map(|l| l.iter().map(|_| 0.).collect()).collect()); - let mut new_coordinates: VectorMultiPolygon = vec![]; - let mut new_offsets: VectorMultiPolygonOffset = vec![]; - coordinates.iter().enumerate().for_each(|(p, polygon)| { - let new_polygon = clip_polygon( - &VectorPolygonGeometry { - _type: VectorGeometryType::Polygon, - is_3d: *is_3d, - coordinates: polygon.to_vec(), - offset: Some(init_o[p].clone()), - ..Default::default() - }, - axis, - k1, - k2, - ); - if let Some(VectorGeometry::Polygon(new_polygon)) = new_polygon { - new_coordinates.push(new_polygon.coordinates); - if let Some(new_offset) = new_polygon.offset { - new_offsets.push(new_offset); - } - } - }); - - if new_coordinates.is_empty() { - None - } else { - Some(VectorGeometry::MultiPolygon(VectorMultiPolygonGeometry { - _type: VectorGeometryType::MultiPolygon, - is_3d: *is_3d, - coordinates: new_coordinates, - bbox: *bbox, - offset: Some(new_offsets), - vec_bbox: Some(vec_bbox.unwrap_or_default().clip(axis, k1, k2)), - ..Default::default() - })) - } -} - -/// After clipping a line, return the altered line, -/// the offset the new line starts at, -/// and if the line is ccw -pub struct ClipLineResult { - /// The clipped line - pub line: VectorLineString, - /// The offset the new line starts at - pub offset: f64, -} -/// Ensuring `vec_bbox` exists -pub struct ClipLineResultWithBBox { - /// The clipped line - pub line: VectorLineString, - /// The offset the new line starts at - pub offset: f64, - /// The new vector bounding box - pub vec_bbox: BBox3D, -} - -/// clip an input line to a bounding box -/// Data should always be in a 0->1 coordinate system to use this clip function -pub fn clip_line( - geom: &VectorLineString, - bbox: BBox3D, - is_polygon: bool, - offset: Option, - buffer: Option, // default for a full size tile. Assuming 1024 extent and a 64 point buffer -) -> Vec { - let offset = offset.unwrap_or(0.); - let buffer = buffer.unwrap_or(0.0625); - let mut res: Vec = vec![]; - let BBox3D { left, bottom, right, top, .. } = bbox; - // clip horizontally - let horizontal_clips = _clip_line( - ClipLineResult { line: geom.clone(), offset }, - left - buffer, - right + buffer, - Axis::X, - is_polygon, - ); - for clip in horizontal_clips { - // clip vertically - let mut vertical_clips = - _clip_line(clip, bottom - buffer, top + buffer, Axis::Y, is_polygon); - res.append(&mut vertical_clips); - } - res.iter_mut() - .map(|clip| { - let mut vec_bbox: Option = None; - for p in clip.line.iter() { - match &mut vec_bbox { - Some(bbox) => bbox.extend_from_point(p), - None => vec_bbox = Some(BBox3D::from_point(p)), - } - } - ClipLineResultWithBBox { - line: core::mem::take(&mut clip.line), - offset: clip.offset, - vec_bbox: vec_bbox.unwrap(), - } - }) - .collect() -} - -/// Interal clip tool -fn _clip_line( - input: ClipLineResult, - k1: f64, - k2: f64, - axis: Axis, - is_polygon: bool, -) -> Vec { - // let { line: geom, offset: startOffset } = input; - let geom = &input.line; - let start_offset = input.offset; - let mut new_geom: Vec = vec![]; - let mut slice: VectorLineString = vec![]; - let mut last = geom.len() - 1; - let intersect = if axis == Axis::X { intersect_x } else { intersect_y }; - - let mut cur_offset = start_offset; - let mut acc_offset = start_offset; - let mut prev_p = &geom[0]; - let mut first_enter = false; - - let mut i = 0; - while i < last { - let VectorPoint { x: ax, y: ay, z: az, m: am, .. } = &geom[i]; - let VectorPoint { x: bx, y: by, z: bz, m: bm, .. } = &geom[i + 1]; - let a: f64 = if axis == Axis::X { *ax } else { *ay }; - let b: f64 = if axis == Axis::X { *bx } else { *by }; - let z: Option = match (az, bz) { - (Some(az), Some(bz)) => Some((az + bz) / 2.0), - (Some(az), None) => Some(*az), - (None, Some(bz)) => Some(*bz), - _ => None, - }; - let mut entered = false; - let mut exited = false; - let mut int_p: Option = None; - - // ENTER OR CONTINUE CASES - if a < k1 { - // ---|--> | (line enters the clip region from the left) - if b > k1 { - int_p = Some(intersect(*ax, *ay, *bx, *by, k1, z, bm)); - slice.push(int_p.clone().unwrap()); - entered = true; - } - } else if a > k2 { - // | <--|--- (line enters the clip region from the right) - if b < k2 { - int_p = Some(intersect(*ax, *ay, *bx, *by, k2, z, bm)); - slice.push(int_p.clone().unwrap()); - entered = true; - } - } else { - int_p = Some(VectorPoint { x: *ax, y: *ay, z: *az, m: am.clone(), t: None }); - slice.push(int_p.clone().unwrap()); - } - - // Update the intersection point and offset if the int_p exists - if let Some(int_p) = int_p.as_ref() { - // our first enter will change the offset for the line - if entered && !first_enter { - cur_offset = acc_offset + prev_p.distance(int_p); - first_enter = true; - } - } - - // EXIT CASES - if b < k1 && a >= k1 { - // <--|--- | or <--|-----|--- (line exits the clip region on the left) - int_p = Some(intersect(*ax, *ay, *bx, *by, k1, z, if bm.is_some() { bm } else { am })); - slice.push(int_p.unwrap()); - exited = true; - } - if b > k2 && a <= k2 { - // | ---|--> or ---|-----|--> (line exits the clip region on the right) - int_p = Some(intersect(*ax, *ay, *bx, *by, k2, z, if bm.is_some() { bm } else { am })); - slice.push(int_p.unwrap()); - exited = true; - } - - // update the offset - acc_offset += prev_p.distance(&geom[i + 1]); - prev_p = &geom[i + 1]; - - // If not a polygon, we can cut it into parts, otherwise we just keep tracking the edges - if !is_polygon && exited { - new_geom.push(ClipLineResult { line: slice, offset: cur_offset }); - slice = vec![]; - first_enter = false; - } - - i += 1; - } - - // add the last point if inside the clip - let last_point = geom[last].clone(); - let a = if axis == Axis::X { last_point.x } else { last_point.y }; - if a >= k1 && a <= k2 { - slice.push(last_point.clone()); - } - - // close the polygon if its endpoints are not the same after clipping - if !slice.is_empty() && is_polygon { - last = slice.len() - 1; - let first_p = &slice[0]; - if last >= 1 && (slice[last].x != first_p.x || slice[last].y != first_p.y) { - slice.push(first_p.clone()); - } - } - - // add the final slice - if !slice.is_empty() { - new_geom.push(ClipLineResult { line: slice, offset: cur_offset }); - } - - new_geom -} - -/// Find the intersection of two points on the X axis -fn intersect_x( - ax: f64, - ay: f64, - bx: f64, - by: f64, - x: f64, - z: Option, - m: &Option, -) -> VectorPoint { - let t = (x - ax) / (bx - ax); - VectorPoint { x, y: ay + (by - ay) * t, z, m: m.clone(), t: Some(1.) } -} - -/// Find the intersection of two points on the Y axis -fn intersect_y( - ax: f64, - ay: f64, - bx: f64, - by: f64, - y: f64, - z: Option, - m: &Option, -) -> VectorPoint { - let t = (y - ay) / (by - ay); - VectorPoint { x: ax + (bx - ax) * t, y, z, m: m.clone(), t: Some(1.) } -} diff --git a/rust/wm/convert.rs b/rust/wm/convert.rs deleted file mode 100644 index 04ce757..0000000 --- a/rust/wm/convert.rs +++ /dev/null @@ -1,703 +0,0 @@ -use crate::{ - build_sq_dists, clip_line, BBox3D, ClipLineResultWithBBox, Face, Feature, Geometry, LonLat, - S2Feature, S2Point, STPoint, VectorFeature, VectorGeometry, VectorGeometryType, - VectorLineString, VectorLineStringGeometry, VectorMultiLineStringGeometry, - VectorMultiPointGeometry, VectorMultiPolygonGeometry, VectorPoint, VectorPointGeometry, - VectorPolygon, VectorPolygonGeometry, -}; - -use alloc::collections::BTreeSet; -use alloc::vec; -use alloc::vec::Vec; - -// TODO: We are cloning geometry twice at times. Let's optimize (check "to_vec" and "clone" cases) - -impl Feature { - /// Convert a GeoJSON Feature to a GeoJSON Vector Feature - pub fn to_vector(data: &Feature, build_bbox: Option) -> VectorFeature { - let build_bbox = build_bbox.unwrap_or(false); - let Feature { id, properties, metadata, geometry, .. } = data; - let vector_geo = convert_geometry(geometry, build_bbox); - VectorFeature::new_wm(*id, properties.clone(), vector_geo, metadata.clone()) - } -} - -impl VectorFeature { - /// Reproject GeoJSON geometry coordinates from lon-lat to a 0->1 coordinate system in place - pub fn to_unit_scale(&mut self, tolerance: Option, maxzoom: Option) { - let mut bbox: Option = None; - match &mut self.geometry { - VectorGeometry::Point(geo) => { - geo.coordinates.project(&mut bbox); - geo.vec_bbox = bbox; - } - VectorGeometry::LineString(geo) | VectorGeometry::MultiPoint(geo) => { - geo.coordinates.iter_mut().for_each(|p| p.project(&mut bbox)); - geo.vec_bbox = bbox; - } - VectorGeometry::Polygon(geo) | VectorGeometry::MultiLineString(geo) => { - geo.coordinates - .iter_mut() - .for_each(|p| p.iter_mut().for_each(|p| p.project(&mut bbox))); - geo.vec_bbox = bbox; - } - VectorGeometry::MultiPolygon(geo) => { - geo.coordinates.iter_mut().for_each(|p| { - p.iter_mut().for_each(|p| p.iter_mut().for_each(|p| p.project(&mut bbox))) - }); - geo.vec_bbox = bbox; - } - } - - if let Some(tolerance) = tolerance { - build_sq_dists(&mut self.geometry, tolerance, maxzoom); - } - } - - /// Reproject GeoJSON geometry coordinates from lon-lat to a 0->1 coordinate system in place - pub fn to_ll(&mut self) { - match &mut self.geometry { - VectorGeometry::Point(geo) => { - geo.coordinates.unproject(); - } - VectorGeometry::LineString(geo) | VectorGeometry::MultiPoint(geo) => { - geo.coordinates.iter_mut().for_each(|p| p.unproject()); - } - VectorGeometry::Polygon(geo) | VectorGeometry::MultiLineString(geo) => { - geo.coordinates.iter_mut().for_each(|p| p.iter_mut().for_each(|p| p.unproject())); - } - VectorGeometry::MultiPolygon(geo) => { - geo.coordinates.iter_mut().for_each(|p| { - p.iter_mut().for_each(|p| p.iter_mut().for_each(|p| p.unproject())) - }); - } - } - } - - /// Convet a GeoJSON Feature to an S2Feature - pub fn to_s2(&self, tolerance: Option, maxzoom: Option) -> Vec> { - let VectorFeature { _type, id, properties, metadata, geometry, .. } = self; - let mut res: Vec> = vec![]; - - if _type == "S2Feature" { - res.push(self.clone()); - } else { - let vector_geo = convert_geometry_wm_to_s2(geometry, tolerance, maxzoom); - for ConvertedGeometry { geometry, face } in vector_geo { - res.push(S2Feature::::new_s2( - *id, - face, - properties.clone(), - geometry, - metadata.clone(), - )); - } - } - - res - } -} - -/// Convert a GeoJSON Geometry to an Vector Geometry -fn convert_geometry(geometry: &Geometry, _build_bbox: bool) -> VectorGeometry { - // TODO: build a bbox if user wants it - match geometry { - Geometry::Point(geo) => { - let mut coordinates: VectorPoint = (&geo.coordinates).into(); - coordinates.m = geo.m_values.clone(); - VectorGeometry::Point(VectorPointGeometry { - _type: VectorGeometryType::Point, - is_3d: false, - coordinates, - bbox: geo.bbox.as_ref().map(|bbox| bbox.into()), - offset: None, - vec_bbox: None, - indices: None, - tesselation: None, - }) - } - Geometry::Point3D(geo) => { - let mut coordinates: VectorPoint = (&geo.coordinates).into(); - coordinates.m = geo.m_values.clone(); - VectorGeometry::Point(VectorPointGeometry { - _type: VectorGeometryType::Point, - is_3d: true, - coordinates, - bbox: geo.bbox, - offset: None, - vec_bbox: None, - indices: None, - tesselation: None, - }) - } - Geometry::MultiPoint(geo) => VectorGeometry::MultiPoint(VectorMultiPointGeometry { - _type: VectorGeometryType::MultiPoint, - is_3d: false, - coordinates: geo - .coordinates - .iter() - .enumerate() - .map(|(i, p)| { - let mut vp = VectorPoint::from(p); - if let Some(m) = &geo.m_values { - vp.m = Some(m[i].clone()); - } - vp - }) - .collect(), - bbox: geo.bbox.as_ref().map(|bbox| bbox.into()), - ..Default::default() - }), - Geometry::MultiPoint3D(geo) => VectorGeometry::MultiPoint(VectorMultiPointGeometry { - _type: VectorGeometryType::MultiPoint, - is_3d: true, - coordinates: geo - .coordinates - .iter() - .enumerate() - .map(|(i, p)| { - let mut vp = VectorPoint::from(p); - if let Some(m) = &geo.m_values { - vp.m = Some(m[i].clone()); - } - vp - }) - .collect(), - bbox: geo.bbox, - ..Default::default() - }), - Geometry::LineString(geo) => VectorGeometry::LineString(VectorLineStringGeometry { - _type: VectorGeometryType::LineString, - is_3d: false, - coordinates: geo - .coordinates - .iter() - .enumerate() - .map(|(i, p)| { - let mut vp = VectorPoint::from(p); - if let Some(m) = &geo.m_values { - vp.m = Some(m[i].clone()); - } - vp - }) - .collect(), - bbox: geo.bbox.as_ref().map(|bbox| bbox.into()), - ..Default::default() - }), - Geometry::LineString3D(geo) => VectorGeometry::LineString(VectorLineStringGeometry { - _type: VectorGeometryType::LineString, - is_3d: true, - coordinates: geo - .coordinates - .iter() - .enumerate() - .map(|(i, p)| { - let mut vp = VectorPoint::from(p); - if let Some(m) = &geo.m_values { - vp.m = Some(m[i].clone()); - } - vp - }) - .collect(), - bbox: geo.bbox, - ..Default::default() - }), - Geometry::MultiLineString(geo) => { - VectorGeometry::MultiLineString(VectorMultiLineStringGeometry { - _type: VectorGeometryType::MultiLineString, - is_3d: false, - coordinates: geo - .coordinates - .iter() - .enumerate() - .map(|(i, line)| { - line.iter() - .enumerate() - .map(|(j, p)| { - let mut vp = VectorPoint::from(p); - if let Some(m) = &geo.m_values { - vp.m = Some(m[i][j].clone()); - } - vp - }) - .collect() - }) - .collect(), - bbox: geo.bbox.as_ref().map(|bbox| bbox.into()), - ..Default::default() - }) - } - Geometry::MultiLineString3D(geo) => { - VectorGeometry::MultiLineString(VectorMultiLineStringGeometry { - _type: VectorGeometryType::MultiLineString, - is_3d: true, - coordinates: geo - .coordinates - .iter() - .enumerate() - .map(|(i, line)| { - line.iter() - .enumerate() - .map(|(j, p)| { - let mut vp = VectorPoint::from(p); - if let Some(m) = &geo.m_values { - vp.m = Some(m[i][j].clone()); - } - vp - }) - .collect() - }) - .collect(), - bbox: geo.bbox, - ..Default::default() - }) - } - Geometry::Polygon(geo) => VectorGeometry::Polygon(VectorPolygonGeometry { - _type: VectorGeometryType::Polygon, - is_3d: false, - coordinates: geo - .coordinates - .iter() - .enumerate() - .map(|(i, line)| { - line.iter() - .enumerate() - .map(|(j, p)| { - let mut vp = VectorPoint::from(p); - if let Some(m) = &geo.m_values { - vp.m = Some(m[i][j].clone()); - } - vp - }) - .collect() - }) - .collect(), - bbox: geo.bbox.as_ref().map(|bbox| bbox.into()), - ..Default::default() - }), - Geometry::Polygon3D(geo) => VectorGeometry::Polygon(VectorPolygonGeometry { - _type: VectorGeometryType::Polygon, - is_3d: true, - coordinates: geo - .coordinates - .iter() - .enumerate() - .map(|(i, line)| { - line.iter() - .enumerate() - .map(|(j, p)| { - let mut vp = VectorPoint::from(p); - if let Some(m) = &geo.m_values { - vp.m = Some(m[i][j].clone()); - } - vp - }) - .collect() - }) - .collect(), - bbox: geo.bbox, - ..Default::default() - }), - Geometry::MultiPolygon(geo) => VectorGeometry::MultiPolygon(VectorMultiPolygonGeometry { - _type: VectorGeometryType::MultiPolygon, - is_3d: false, - coordinates: geo - .coordinates - .iter() - .enumerate() - .map(|(i, polygon)| { - polygon - .iter() - .enumerate() - .map(|(j, line)| { - line.iter() - .enumerate() - .map(|(k, p)| { - let mut vp = VectorPoint::from(p); - if let Some(m) = &geo.m_values { - vp.m = Some(m[i][j][k].clone()); - } - vp - }) - .collect() - }) - .collect() - }) - .collect(), - bbox: geo.bbox.as_ref().map(|bbox| bbox.into()), - ..Default::default() - }), - Geometry::MultiPolygon3D(geo) => VectorGeometry::MultiPolygon(VectorMultiPolygonGeometry { - _type: VectorGeometryType::MultiPolygon, - is_3d: true, - coordinates: geo - .coordinates - .iter() - .enumerate() - .map(|(i, polygon)| { - polygon - .iter() - .enumerate() - .map(|(j, line)| { - line.iter() - .enumerate() - .map(|(k, p)| { - let mut vp = VectorPoint::from(p); - if let Some(m) = &geo.m_values { - vp.m = Some(m[i][j][k].clone()); - } - vp - }) - .collect() - }) - .collect() - }) - .collect(), - bbox: geo.bbox, - ..Default::default() - }), - } -} - -/// The resultant geometry after conversion -pub struct ConvertedGeometry { - pub geometry: VectorGeometry, - pub face: Face, -} -pub type ConvertedGeometryList = Vec; - -/// Underlying conversion mechanic to move GeoJSON Geometry to S2Geometry -fn convert_geometry_wm_to_s2( - geometry: &VectorGeometry, - tolerance: Option, - maxzoom: Option, -) -> ConvertedGeometryList { - let mut res: ConvertedGeometryList = vec![]; - - match geometry { - VectorGeometry::Point(geo) => { - res.extend(convert_geometry_point(geo)); - } - VectorGeometry::MultiPoint(geo) => { - res.extend(convert_geometry_multipoint(geo)); - } - VectorGeometry::LineString(geo) => { - res.extend(convert_geometry_linestring(geo)); - } - VectorGeometry::MultiLineString(geo) => { - res.extend(convert_geometry_multilinestring(geo)); - } - VectorGeometry::Polygon(geo) => { - res.extend(convert_geometry_polygon(geo)); - } - VectorGeometry::MultiPolygon(geo) => { - res.extend(convert_geometry_multipolygon(geo)); - } - } - - if let Some(tolerance) = tolerance { - for c_geo in &mut res { - build_sq_dists(&mut c_geo.geometry, tolerance, maxzoom); - } - } - - res -} - -// /** -// * @param geometry - GeoJSON PointGeometry -// * @returns - S2 PointGeometry -// */ -/// Convert a GeoJSON PointGeometry to a S2 PointGeometry -fn convert_geometry_point(geometry: &VectorPointGeometry) -> ConvertedGeometryList { - let VectorPointGeometry { _type, is_3d, coordinates, bbox, .. } = geometry; - let mut new_point = coordinates.clone(); - let ll: S2Point = (&LonLat::new(new_point.x, new_point.y)).into(); - let (face, s, t) = ll.to_face_st(); - new_point.x = s; - new_point.y = t; - let vec_bbox = Some(BBox3D::from_point(&new_point)); - vec![ConvertedGeometry { - face: face.into(), - geometry: VectorGeometry::Point(VectorPointGeometry { - _type: VectorGeometryType::Point, - coordinates: new_point, - is_3d: *is_3d, - bbox: *bbox, - vec_bbox, - offset: None, - indices: None, - tesselation: None, - }), - }] -} - -// /** -// * @param geometry - GeoJSON PointGeometry -// * @returns - S2 PointGeometry -// */ -fn convert_geometry_multipoint(geometry: &VectorMultiPointGeometry) -> ConvertedGeometryList { - let VectorMultiPointGeometry { is_3d, coordinates, bbox, .. } = geometry; - coordinates - .iter() - .flat_map(|coordinates| { - convert_geometry_point(&VectorPointGeometry { - _type: VectorGeometryType::Point, - is_3d: *is_3d, - coordinates: coordinates.clone(), - bbox: *bbox, - offset: None, - vec_bbox: None, - indices: None, - tesselation: None, - }) - }) - .collect() -} - -/// Convert a GeoJSON LineStringGeometry to S2 LineStringGeometry -fn convert_geometry_linestring(geometry: &VectorLineStringGeometry) -> ConvertedGeometryList { - let VectorLineStringGeometry { _type, is_3d, coordinates, bbox, .. } = geometry; - - convert_line_string(coordinates, false) - .into_iter() - .map(|cline| { - let ConvertedLineString { face, line, offset, vec_bbox } = cline; - ConvertedGeometry { - face, - geometry: VectorGeometry::LineString(VectorLineStringGeometry { - _type: VectorGeometryType::LineString, - is_3d: *is_3d, - coordinates: line.to_vec(), - bbox: *bbox, - offset: Some(offset), - vec_bbox: Some(vec_bbox), - ..Default::default() - }), - } - }) - .collect() -} - -/// Convert a GeoJSON MultiLineStringGeometry to S2 MultiLineStringGeometry -fn convert_geometry_multilinestring( - geometry: &VectorMultiLineStringGeometry, -) -> ConvertedGeometryList { - let VectorMultiLineStringGeometry { is_3d, coordinates, bbox, .. } = geometry; - - coordinates - .iter() - .flat_map(|line| convert_line_string(line, false)) - .map(|ConvertedLineString { face, line, offset, vec_bbox }| ConvertedGeometry { - face, - geometry: VectorGeometry::LineString(VectorLineStringGeometry { - _type: VectorGeometryType::LineString, - is_3d: *is_3d, - coordinates: line, - bbox: *bbox, - offset: Some(offset), - vec_bbox: Some(vec_bbox), - ..Default::default() - }), - }) - .collect() -} - -/// Convert a GeoJSON PolygonGeometry to S2 PolygonGeometry -fn convert_geometry_polygon(geometry: &VectorPolygonGeometry) -> ConvertedGeometryList { - let VectorPolygonGeometry { _type, is_3d, coordinates, bbox, .. } = geometry; - let mut res: ConvertedGeometryList = vec![]; - - // conver all lines - let mut outer_ring = convert_line_string(&coordinates[0], true); - let mut inner_rings = coordinates[1..].iter().flat_map(|line| convert_line_string(line, true)); - - // for each face, build a new polygon - for ConvertedLineString { face, line, offset, vec_bbox: poly_bbox } in &mut outer_ring { - let mut polygon: VectorPolygon = vec![line.to_vec()]; - let mut polygon_offsets = vec![*offset]; - let mut poly_bbox = *poly_bbox; - for ConvertedLineString { - face: inner_face, - line: inner_line, - offset: inner_offset, - vec_bbox, - } in &mut inner_rings - { - if inner_face == *face { - polygon.push(inner_line); - polygon_offsets.push(inner_offset); - poly_bbox = poly_bbox.merge(&vec_bbox); - } - } - - res.push(ConvertedGeometry { - face: *face, - geometry: VectorGeometry::Polygon(VectorPolygonGeometry { - _type: VectorGeometryType::Polygon, - is_3d: *is_3d, - coordinates: polygon, - bbox: *bbox, - offset: Some(polygon_offsets), - vec_bbox: Some(poly_bbox), - ..Default::default() - }), - }); - } - - res -} - -/// Convert a GeoJSON MultiPolygonGeometry to S2 MultiPolygonGeometry -fn convert_geometry_multipolygon(geometry: &VectorMultiPolygonGeometry) -> ConvertedGeometryList { - let VectorMultiPolygonGeometry { is_3d, coordinates, bbox, offset, .. } = geometry; - coordinates - .iter() - .enumerate() - .flat_map(|(i, polygon)| { - let offset: Option> = offset.as_ref().map(|offset| offset[i].clone()); - convert_geometry_polygon(&VectorPolygonGeometry { - _type: VectorGeometryType::Polygon, - is_3d: *is_3d, - coordinates: polygon.to_vec(), - bbox: *bbox, - offset, - ..Default::default() - }) - }) - .collect() -} - -/// LineString converted from WM to S2 -pub struct ConvertedLineString { - face: Face, - line: VectorLineString, - offset: f64, - vec_bbox: BBox3D, -} - -/// Convert WM LineString to S2 -fn convert_line_string(line: &VectorLineString, is_polygon: bool) -> Vec { - let mut res: Vec = vec![]; - // first re-project all the coordinates to S2 - let mut new_geometry: Vec = vec![]; - for VectorPoint { x: lon, y: lat, z, m, .. } in line { - let ll: S2Point = (&LonLat::new(*lon, *lat)).into(); - let (face, s, t) = ll.to_face_st(); - new_geometry.push(STPoint { face: face.into(), s, t, z: *z, m: m.clone() }); - } - // find all the faces that exist in the line - let mut faces = BTreeSet::::new(); - new_geometry.iter().for_each(|stpoint| { - faces.insert(stpoint.face); - }); - // for each face, build a line - for face in faces { - let mut line: VectorLineString = vec![]; - for st_point in &new_geometry { - line.push(st_point_to_face(face, st_point)); - } - let clipped_lines = clip_line(&line, BBox3D::default(), is_polygon, None, None); - for ClipLineResultWithBBox { line, offset, vec_bbox } in clipped_lines { - res.push(ConvertedLineString { face, line, offset, vec_bbox }); - } - } - - res -} - -/// Given a face, rotate the point into it's 0->1 coordinate system -fn st_point_to_face(target_face: Face, stp: &STPoint) -> VectorPoint { - let cur_face = stp.face; - if target_face == cur_face { - return VectorPoint { x: stp.s, y: stp.t, z: stp.z, m: stp.m.clone(), t: None }; - } - - let (rot, x, y) = &FACE_RULE_SET[target_face as usize][cur_face as usize]; - let (new_s, new_t) = rotate(*rot, stp.s, stp.t); - - VectorPoint { x: new_s + *x as f64, y: new_t + *y as f64, z: stp.z, m: stp.m.clone(), t: None } -} - -/** - * @param rot - rotation - * @param s - input s - * @param t - input t - * @returns - new [s, t] after rotating - */ -fn rotate(rot: Rotation, s: f64, t: f64) -> (f64, f64) { - match rot { - Rotation::_0 => (s, t), - Rotation::_90 => (t, 1. - s), - Rotation::_Neg90 => (1. - t, s), - } -} - -#[derive(Debug, PartialEq, Copy, Clone)] -pub enum Rotation { - _0, - _90, - _Neg90, -} - -/// Ruleset for converting an S2Point from a face to another. -/// While this this set includes opposite side faces, without axis mirroring, -/// it is not technically accurate and shouldn't be used. Instead, data should let two points travel -/// further than a full face width. -/// FACE_RULE_SET[target_face][currentFace] = [rot, x, y] -pub const FACE_RULE_SET: [[(Rotation, i8, i8); 6]; 6] = [ - // Target Face 0 - [ - (Rotation::_0, 0, 0), // Current Face 0 - (Rotation::_0, 1, 0), // Current Face 1 - (Rotation::_90, 0, 1), // Current Face 2 - (Rotation::_Neg90, 2, 0), // Current Face 3 - (Rotation::_Neg90, -1, 0), // Current Face 4 - (Rotation::_0, 0, -1), // Current Face 5 - ], - // Target Face 1 - [ - (Rotation::_0, -1, 0), // Current Face 0 - (Rotation::_0, 0, 0), // Current Face 1 - (Rotation::_0, 0, 1), // Current Face 2 - (Rotation::_Neg90, 1, 0), // Current Face 3 - (Rotation::_Neg90, 2, 0), // Current Face 4 - (Rotation::_90, 0, -1), // Current Face 5 - ], - // Target Face 2 - [ - (Rotation::_Neg90, -1, 0), // Current Face 0 - (Rotation::_0, 0, -1), // Current Face 1 - (Rotation::_0, 0, 0), // Current Face 2 - (Rotation::_0, 1, 0), // Current Face 3 - (Rotation::_90, 0, 1), // Current Face 4 - (Rotation::_Neg90, 2, 0), // Current Face 5 - ], - // Target Face 3 - [ - (Rotation::_Neg90, 2, 0), // Current Face 0 - (Rotation::_90, 0, -1), // Current Face 1 - (Rotation::_0, -1, 0), // Current Face 2 - (Rotation::_0, 0, 0), // Current Face 3 - (Rotation::_0, 0, 1), // Current Face 4 - (Rotation::_Neg90, 1, 0), // Current Face 5 - ], - // Target Face 4 - [ - (Rotation::_90, 0, 1), // Current Face 0 - (Rotation::_Neg90, 2, 0), // Current Face 1 - (Rotation::_Neg90, -1, 0), // Current Face 2 - (Rotation::_0, 0, -1), // Current Face 3 - (Rotation::_0, 0, 0), // Current Face 4 - (Rotation::_0, 1, 0), // Current Face 5 - ], - // Target Face 5 - [ - (Rotation::_0, 0, 1), // Current Face 0 - (Rotation::_Neg90, 1, 0), // Current Face 1 - (Rotation::_Neg90, 2, 0), // Current Face 2 - (Rotation::_90, 0, -1), // Current Face 3 - (Rotation::_0, -1, 0), // Current Face 4 - (Rotation::_0, 0, 0), // Current Face 5 - ], -]; diff --git a/rust/wm/coords.rs b/rust/wm/coords.rs deleted file mode 100644 index dc9ef55..0000000 --- a/rust/wm/coords.rs +++ /dev/null @@ -1,275 +0,0 @@ -use libm::{atan, cos, exp, floor, log, pow, sin, tan}; - -use crate::{A, EARTH_CIRCUMFERENCE, MAXEXTENT}; - -use core::f64::consts::PI; - -/// The source of the coordinate inputs -#[derive(Debug, Clone, PartialEq)] -pub enum Source { - /// The WGS84 projection - WGS84, - /// The Google (900913) projection - Google, -} - -/// Given a zoom and tilesize, build mercator positional attributes -pub fn get_zoom_size(zoom: u8, tile_size: f64) -> (f64, f64, f64, f64) { - let size = tile_size * pow(2., zoom as f64); - (size / 360., size / (2. * PI), size / 2., size) -} - -/// Convert Longitude and Latitude to a mercator pixel coordinate -/// Return the mercator pixel coordinate -pub fn ll_to_px( - lonlat: (f64, f64), - zoom: u8, - anti_meridian: Option, - tile_size: Option, -) -> (f64, f64) { - let anti_meridian = anti_meridian.unwrap_or(false); - let tile_size = tile_size.unwrap_or(512.); - - let (bc, cc, zc, ac) = get_zoom_size(zoom, tile_size); - let expansion = if anti_meridian { 2. } else { 1. }; - let d = zc; - let f = sin((lonlat.1).to_radians()).clamp(-0.999999999999, 0.999999999999); - let mut x = d + lonlat.0 * bc; - let mut y = d + 0.5 * log((1. + f) / (1. - f)) * -cc; - if x > ac * expansion { - x = ac * expansion; - } - if y > ac { - y = ac; - } - - (x, y) -} - -/// Convert mercator pixel coordinates to Longitude and Latitude -/// Return the Longitude and Latitude -pub fn px_to_ll(xy: (f64, f64), zoom: u8, tile_size: Option) -> (f64, f64) { - let tile_size = tile_size.unwrap_or(512.); - let (bc, cc, zc, _) = get_zoom_size(zoom, tile_size); - let g = (xy.1 - zc) / -cc; - let lon = (xy.0 - zc) / bc; - let lat = 2. * atan(exp(g)).to_degrees() - 0.5 * PI; - - (lon, lat) -} - -/// Convert Longitude and Latitude to a mercator x-y coordinates -/// Return the mercator x-y coordinates -pub fn ll_to_merc(lonlan: (f64, f64)) -> (f64, f64) { - let mut x = (A * lonlan.0).to_radians(); - let mut y = A * log(tan(PI * 0.25 + (0.5 * lonlan.1).to_radians())); - // if xy value is beyond maxextent (e.g. poles), return maxextent. - x = x.clamp(-MAXEXTENT, MAXEXTENT); - y = y.clamp(-MAXEXTENT, MAXEXTENT); - - (x, y) -} - -/// Convert mercator x-y coordinates to Longitude and Latitude -/// Return the Longitude and Latitude -pub fn merc_to_ll(merc: (f64, f64)) -> (f64, f64) { - let x = (merc.0 / A).to_degrees(); - let y = (0.5 * PI - 2. * atan(exp(-merc.1 / A))).to_degrees(); - (x, y) -} - -/// Convert a pixel coordinate to a tile x-y coordinate -/// Return the tile x-y -pub fn px_to_tile(px: (f64, f64), tile_size: Option) -> (u32, u32) { - let tile_size = tile_size.unwrap_or(512.); - (floor(px.0 / tile_size) as u32, floor(px.1 / tile_size) as u32) -} - -/// Convert a tile x-y-z to a bbox of the form `[w, s, e, n]` -/// Return the bbox -pub fn tile_to_bbox(tile: (u8, u32, u32), tile_size: Option) -> (f64, f64, f64, f64) { - let tile_size = tile_size.unwrap_or(512.); - let (_zoom, x, y) = tile; - let min_x = x as f64 * tile_size; - let min_y = y as f64 * tile_size; - let max_x = min_x + tile_size; - let max_y = min_y + tile_size; - - (min_x, min_y, max_x, max_y) -} - -/// Convert a lat-lon and zoom to the tile's x-y coordinates -/// Return the tile x-y -pub fn ll_to_tile(lonlat: (f64, f64), zoom: u8, tile_size: Option) -> (u32, u32) { - let px = ll_to_px(lonlat, zoom, Some(false), tile_size); - px_to_tile(px, tile_size) -} - -/// given a lon-lat and tile, find the offset in pixels -/// return the tile xy pixel -pub fn ll_to_tile_px( - lonlat: (f64, f64), - tile: (u8, u32, u32), - tile_size: Option, -) -> (f64, f64) { - let (zoom, x, y) = tile; - let tile_size = tile_size.unwrap_or(512.); - let px = ll_to_px(lonlat, zoom, Some(false), Some(tile_size)); - let tile_x_start = x as f64 * tile_size; - let tile_y_start = y as f64 * tile_size; - - ((px.0 - tile_x_start) / tile_size, (px.1 - tile_y_start) / tile_size) -} - -/// Convert a bbox of the form `[w, s, e, n]` to a bbox of the form `[w, s, e, n]` -/// The result can be in lon-lat (WGS84) or WebMercator (900913) -pub fn convert_bbox(bbox: &(f64, f64, f64, f64), source: Source) -> (f64, f64, f64, f64) { - let low: (f64, f64); - let high: (f64, f64); - match source { - Source::WGS84 => { - low = merc_to_ll((bbox.0, bbox.1)); - high = merc_to_ll((bbox.2, bbox.3)); - } - Source::Google => { - low = ll_to_merc((bbox.0, bbox.1)); - high = ll_to_merc((bbox.2, bbox.3)); - } - }; - (low.0, low.1, high.0, high.1) -} - -/// Convert a tile x-y-z to a bbox of the form `[w, s, e, n]` -/// The result can be in lon-lat (WGS84) or WebMercator (900913) -/// The default result is in WebMercator (900913) -pub fn xyz_to_bbox( - x: f64, - y: f64, - zoom: u8, - tms_style: Option, - source: Option, - tile_size: Option, -) -> (f64, f64, f64, f64) { - let tms_style = tms_style.unwrap_or(true); - let source = source.unwrap_or(Source::Google); - let tile_size = tile_size.unwrap_or(512.0); - let mut y = y; - // Convert xyz into bbox with srs WGS84 - // if tmsStyle, the y is inverted - if tms_style { - y = pow(2., zoom as f64) - 1. - y; - } - // Use +y to make sure it's a number to avoid inadvertent concatenation. - let bl: (f64, f64) = (x * tile_size, (y + 1.) * tile_size); - // Use +x to make sure it's a number to avoid inadvertent concatenation. - let tr: (f64, f64) = ((x + 1.) * tile_size, y * tile_size); - // to pixel-coordinates - let px_bl = px_to_ll(bl, zoom, Some(tile_size)); - let px_tr = px_to_ll(tr, zoom, Some(tile_size)); - - match source { - Source::Google => { - let ll_bl = ll_to_merc(px_bl); - let ll_tr = ll_to_merc(px_tr); - (ll_bl.0, ll_bl.1, ll_tr.0, ll_tr.1) - } - _ => (px_bl.0, px_bl.1, px_tr.0, px_tr.1), - } -} - -/// Convert a bbox of the form `[w, s, e, n]` to a tile's bounding box -/// in the form of { minX, maxX, minY, maxY } -/// The bbox can be in lon-lat (WGS84) or WebMercator (900913) -/// The default expectation is in WebMercator (900913) -/// returns the tile's bounding box -pub fn bbox_to_xyz_bounds( - bbox: (f64, f64, f64, f64), - zoom: u8, - tms_style: Option, - source: Option, - tile_size: Option, -) -> (f64, f64, f64, f64) { - let tms_style = tms_style.unwrap_or(true); - let source = source.unwrap_or(Source::Google); - let tile_size = tile_size.unwrap_or(512.0); - - let mut bl = (bbox.0, bbox.1); // bottom left - let mut tr = (bbox.2, bbox.3); // top right - - if source == Source::Google { - bl = ll_to_merc(bl); - tr = ll_to_merc(tr); - } - let px_bl = ll_to_px(bl, zoom, Some(false), Some(tile_size)); - let px_tr = ll_to_px(tr, zoom, Some(false), Some(tile_size)); - let x = (floor(px_bl.0) / tile_size, floor((px_tr.0 - 1.0) / tile_size)); - let y = (floor(px_tr.1) / tile_size, floor((px_bl.1 - 1.0) / tile_size)); - - let mut bounds = - (f64::min(x.0, x.1), f64::min(y.0, y.1), f64::max(x.0, x.1), f64::max(y.0, y.1)); - if bounds.0 < 0. { - bounds.0 = 0. - } - if bounds.1 < 0. { - bounds.1 = 0. - } - - if tms_style { - let zoom_diff = pow(2., zoom as f64) - 1.; - bounds.1 = zoom_diff - bounds.3; - bounds.3 = zoom_diff - bounds.1; - } - - bounds -} - -/// The circumference at a line of latitude in meters. -pub fn circumference_at_latitude(latitude: f64) -> f64 { - EARTH_CIRCUMFERENCE * cos(latitude.to_radians()) -} - -/// Convert longitude to mercator projection X-Value -/// returns the X-Value -pub fn lng_to_mercator_x(lng: f64) -> f64 { - (180.0 + lng) / 360.0 -} - -/// Convert latitude to mercator projection Y-Value -/// returns the Y-Value -pub fn lat_to_mercator_y(lat: f64) -> f64 { - (180. - (180. / PI) * log(tan(PI / 4. + (lat * PI) / 360.))) / 360. -} - -/// Convert altitude to mercator projection Z-Value -/// returns the Z-Value -pub fn altitude_to_mercator_z(altitude: f64, lat: f64) -> f64 { - altitude / circumference_at_latitude(lat) -} - -/// Convert mercator projection's X-Value to longitude -/// returns the longitude -pub fn lng_from_mercator_x(x: f64) -> f64 { - x * 360. - 180. -} - -/// Convert mercator projection's Y-Value to latitude -/// returns the latitude -pub fn lat_from_mercator_y(y: f64) -> f64 { - let y2 = 180. - y * 360.; - (360. / PI) * atan(exp((y2 * PI) / 180.)) - 90. -} - -/// Convert mercator projection's Z-Value to altitude -/// returns the altitude -pub fn altitude_from_mercator_z(z: f64, y: f64) -> f64 { - z * circumference_at_latitude(lat_from_mercator_y(y)) -} - -/// Determine the Mercator scale factor for a given latitude, see -/// https://en.wikipedia.org/wiki/Mercator_projection#Scale_factor -/// -/// At the equator the scale factor will be 1, which increases at higher latitudes. -/// returns the scale factor -pub fn mercator_lat_scale(lat: f64) -> f64 { - 1. / cos((lat * PI) / 180.) -} diff --git a/rust/wm/lonlat.rs b/rust/wm/lonlat.rs deleted file mode 100644 index ec9b4db..0000000 --- a/rust/wm/lonlat.rs +++ /dev/null @@ -1,238 +0,0 @@ -use libm::{asin, atan2, cos, fabs, sin, sqrt}; - -use core::cmp::Ordering; -use core::f64::consts::PI; -use core::ops::{Add, Div, Mul, Neg, Sub}; - -use crate::s2::{S2CellId, S2Point}; - -/// This class represents a point on the unit sphere as a pair -/// of latitude-longitude coordinates. Like the rest of the "geometry" -/// package, the intent is to represent spherical geometry as a mathematical -/// abstraction, so functions that are specifically related to the Earth's -/// geometry (e.g. easting/northing conversions) should be put elsewhere. -#[derive(Clone, Copy, Default, PartialEq, Debug)] -pub struct LonLat { - /// longitude in degrees - pub lon: f64, - /// latitude in degrees - pub lat: f64, -} -impl LonLat { - /// The default constructor sets the latitude and longitude to zero. This is - /// mainly useful when declaring arrays, STL containers, etc. - pub fn new(lon: f64, lat: f64) -> Self { - LonLat { lon, lat } - } - - /// Convert a S2CellId to an LonLat - pub fn from_s2cellid(cellid: &S2CellId) -> LonLat { - let p = cellid.to_point(); - LonLat::from_s2_point(&p) - } - - /// Convert a direction vector (not necessarily unit length) to an LonLat. - pub fn from_s2_point(p: &S2Point) -> LonLat { - let lon_rad = atan2(p.y + 0.0, p.x + 0.0); - let lat_rad = atan2(p.z, sqrt(p.x * p.x + p.y * p.y)); - let ll = LonLat::new((lon_rad).to_degrees(), (lat_rad).to_degrees()); - if !ll.is_valid() { - unreachable!(); - } - ll - } - - /// return the value of the axis - pub fn from_axis(&self, axis: u8) -> f64 { - if axis == 0 { - self.lon - } else { - self.lat - } - } - - /// Normalize the coordinates to the range [-180, 180] and [-90, 90] deg. - pub fn normalize(&mut self) { - let mut lon = self.lon; - let mut lat = self.lat; - while lon < -180. { - lon += 360. - } - while lon > 180. { - lon -= 360. - } - while lat < -90. { - lat += 180. - } - while lat > 90. { - lat -= 180. - } - } - - /// Return the latitude or longitude coordinates in degrees. - pub fn coords(self) -> (f64, f64) { - (self.lon, self.lat) - } - - /// Return true if the latitude is between -90 and 90 degrees inclusive - /// and the longitude is between -180 and 180 degrees inclusive. - pub fn is_valid(&self) -> bool { - fabs(self.lat) <= (PI / 2.0) && fabs(self.lon) <= PI - } - - // // Clamps the latitude to the range [-90, 90] degrees, and adds or subtracts - // // a multiple of 360 degrees to the longitude if necessary to reduce it to - // // the range [-180, 180]. - // LonLat Normalized() const; - - /// Converts an LonLat to the equivalent unit-length vector. Unnormalized - /// values (see Normalize()) are wrapped around the sphere as would be expected - /// based on their definition as spherical angles. So for example the - /// following pairs yield equivalent points (modulo numerical error): - /// (90.5, 10) =~ (89.5, -170) - /// (a, b) =~ (a + 360 * n, b) - /// The maximum error in the result is 1.5 * DBL_EPSILON. (This does not - /// include the error of converting degrees, E5, E6, or E7 to radians.) - /// - /// Can be used just like an S2Point constructor. For example: - /// S2Cap cap; - /// cap.AddPoint(S2Point(latlon)); - pub fn to_point(&self) -> S2Point { - if !self.is_valid() { - unreachable!(); - } - let lon: f64 = (self.lon).to_degrees(); - let lat: f64 = (self.lat).to_degrees(); - S2Point::new( - cos(lat) * cos(lon), // x - cos(lat) * sin(lon), // y - sin(lat), // z - ) - } - - /// An alternative to to_point() that returns a GPU compatible vector. - pub fn to_point_gl(&self) -> S2Point { - let lon: f64 = (self.lon).to_degrees(); - let lat: f64 = (self.lat).to_degrees(); - S2Point::new( - cos(lat) * sin(lon), // y - sin(lat), // z - cos(lat) * cos(lon), // x - ) - } - - /// Returns the distance (measured along the surface of the sphere) to the - /// given LonLat, implemented using the Haversine formula. This is - /// equivalent to - /// - /// S1Angle(ToPoint(), o.ToPoint()) - /// - /// except that this function is slightly faster, and is also somewhat less - /// accurate for distances approaching 180 degrees (see s1angle.h for - /// details). Both LngLats must be normalized. - pub fn get_distance(&self, b: &LonLat) -> f64 { - // This implements the Haversine formula, which is numerically stable for - // small distances but only gets about 8 digits of precision for very large - // distances (e.g. antipodal points). Note that 8 digits is still accurate - // to within about 10cm for a sphere the size of the Earth. - // - // This could be fixed with another sin() and cos() below, but at that point - // you might as well just convert both arguments to S2Points and compute the - // distance that way (which gives about 15 digits of accuracy for all - // distances). - if !self.is_valid() || !b.is_valid() { - unreachable!(); - } - - let lat1 = self.lat; - let lat2 = b.lat; - let lon1 = self.lon; - let lon2 = b.lon; - let dlat = sin(0.5 * (lat2 - lat1)); - let dlon = sin(0.5 * (lon2 - lon1)); - let x = dlat * dlat + dlon * dlon * cos(lat1) * cos(lat2); - 2. * asin(sqrt(f64::min(1., x))) - } -} -impl From<&S2CellId> for LonLat { - fn from(c: &S2CellId) -> Self { - LonLat::from_s2cellid(c) - } -} -impl From<&S2Point> for LonLat { - fn from(p: &S2Point) -> Self { - LonLat::from_s2_point(p) - } -} -impl Add for LonLat { - type Output = LonLat; - fn add(self, rhs: f64) -> Self::Output { - LonLat::new(self.lat + rhs, self.lon + rhs) - } -} -impl Add for LonLat { - type Output = LonLat; - fn add(self, rhs: LonLat) -> Self::Output { - LonLat::new(self.lat + rhs.lat, self.lon + rhs.lon) - } -} -impl Sub for LonLat { - type Output = LonLat; - fn sub(self, rhs: LonLat) -> Self::Output { - LonLat::new(self.lat - rhs.lat, self.lon - rhs.lon) - } -} -impl Sub for LonLat { - type Output = LonLat; - fn sub(self, rhs: f64) -> Self::Output { - LonLat::new(self.lat - rhs, self.lon - rhs) - } -} -impl Mul for LonLat { - type Output = LonLat; - fn mul(self, rhs: f64) -> Self::Output { - LonLat::new(self.lat * rhs, self.lon * rhs) - } -} -impl Mul for f64 { - type Output = LonLat; - fn mul(self, rhs: LonLat) -> Self::Output { - LonLat::new(self * rhs.lat, self * rhs.lon) - } -} -impl Div for LonLat { - type Output = LonLat; - fn div(self, rhs: f64) -> Self::Output { - LonLat::new(self.lat / rhs, self.lon / rhs) - } -} -impl Div for LonLat { - type Output = LonLat; - fn div(self, rhs: LonLat) -> Self::Output { - LonLat::new(self.lat / rhs.lat, self.lon / rhs.lon) - } -} -impl Neg for LonLat { - type Output = LonLat; - fn neg(self) -> Self::Output { - LonLat::new(-self.lat, -self.lon) - } -} -impl Eq for LonLat {} -impl Ord for LonLat { - fn cmp(&self, other: &Self) -> Ordering { - match self.lon.partial_cmp(&other.lon) { - Some(Ordering::Equal) => {} - other => return other.unwrap(), // Handle cases where `lon` comparison is not equal - } - match self.lat.partial_cmp(&other.lat) { - Some(order) => order, - None => Ordering::Equal, // This handles the NaN case safely - } - } -} -impl PartialOrd for LonLat { - fn partial_cmp(&self, other: &Self) -> Option { - Some(self.cmp(other)) - } -} diff --git a/rust/wm/mod.rs b/rust/wm/mod.rs deleted file mode 100644 index 9fbd1cb..0000000 --- a/rust/wm/mod.rs +++ /dev/null @@ -1,8 +0,0 @@ -mod clip; -mod convert; -mod coords; -mod lonlat; - -pub use clip::*; -pub use coords::*; -pub use lonlat::*;