gistools/geometry/s2/
point.rs

1use super::{ST_TO_UV, face_uv_to_xyz};
2use crate::geometry::{LonLat, S2CellId, xyz_to_face_st, xyz_to_face_uv};
3use core::{
4    cmp::Ordering,
5    fmt::Debug,
6    ops::{Add, AddAssign, Div, DivAssign, Mul, MulAssign, Neg, Rem, RemAssign, Sub, SubAssign},
7};
8use libm::{atan2, fabs, sqrt};
9use s2json::{GetXY, GetZ, NewXY, NewXYZ, SetXY, SetZ, VectorPoint};
10use serde::{Deserialize, Serialize};
11
12/// An S2Point represents a point on the unit sphere as a 3D vector. Usually
13/// points are normalized to be unit length, but some methods do not require
14/// this.  See util/math/vector.h for the methods available.  Among other
15/// things, there are overloaded operators that make it convenient to write
16/// arithmetic expressions (e.g. p1 + p2).
17///
18/// NOTE: asumes only f64 or greater is used.
19///
20/// Implements the [`GetXY`] and [`GetZ`] traits
21///
22/// ## Usage
23///
24/// Methods that are available:
25/// - [`S2Point::new`]: Create a new S2Point
26/// - [`S2Point::is_empty`]: Check if the S2Point is empty
27/// - [`S2Point::face`]: Returns the S2 face assocated with this point
28/// - [`S2Point::angle`]: Returns the angle between this point and another
29/// - [`S2Point::cross`]: Get the cross product of two XYZ Points
30/// - [`S2Point::to_face_st`]: Convert an S2Point to an S2Point in normalized vector coordinates
31/// - [`S2Point::get_face`]: Returns the S2 face assocated with this point
32/// - [`S2Point::dot`]: Get the dot product of two XYZ Points
33/// - [`S2Point::abs`]: Returns the absolute value of the point
34/// - [`S2Point::invert`]: Inverts the point
35/// - [`S2Point::len`]: Returns the length of the point
36/// - [`S2Point::norm`]: Returns the vector's squared norm.
37/// - [`S2Point::norm2`]: The dot product of the point with itself
38/// - [`S2Point::normalize`]: Normalizes the point
39/// - [`S2Point::distance`]: return the distance from this point to the other point
40/// - [`S2Point::largest_abs_component`]: Returns the largest absolute component of the point
41/// - [`S2Point::intermediate`]: Returns the intermediate point between this and the other point
42/// - [`S2Point::from_face_uv`]: Convert an Face-U-V coordinate to an S2Point
43/// - [`S2Point::from_face_st`]: Convert an Face-S-T coordinate to an S2Point
44#[derive(Debug, Copy, Clone, Default, PartialEq, Serialize, Deserialize)]
45#[repr(C)]
46pub struct S2Point {
47    /// The x component.
48    pub x: f64,
49    /// The y component.
50    pub y: f64,
51    /// The z component.
52    pub z: f64,
53}
54impl GetXY for S2Point {
55    fn x(&self) -> f64 {
56        self.x
57    }
58    fn y(&self) -> f64 {
59        self.y
60    }
61}
62impl GetZ for S2Point {
63    fn z(&self) -> Option<f64> {
64        Some(self.z)
65    }
66}
67impl NewXY for S2Point {
68    fn new_xy(x: f64, y: f64) -> Self {
69        S2Point { x, y, z: 0.0 }
70    }
71}
72impl NewXYZ for S2Point {
73    fn new_xyz(x: f64, y: f64, z: f64) -> Self {
74        S2Point { x, y, z }
75    }
76}
77impl SetXY for S2Point {
78    fn set_x(&mut self, x: f64) {
79        self.x = x;
80    }
81    fn set_y(&mut self, y: f64) {
82        self.y = y;
83    }
84}
85impl SetZ for S2Point {
86    fn set_z(&mut self, z: f64) {
87        self.z = z;
88    }
89}
90impl S2Point {
91    /// Creates a new S2Point.
92    pub fn new(x: f64, y: f64, z: f64) -> Self {
93        S2Point { x, y, z }
94    }
95
96    /// Returns true if the point is the zero vector.
97    pub fn is_empty(&self) -> bool {
98        let zero = f64::default();
99        self.x == zero && self.y == zero && self.z == zero
100    }
101
102    /// Returns the S2 face assocated with this point
103    pub fn face(&self, f: u8) -> f64 {
104        if f == 0 {
105            self.x
106        } else if f == 1 {
107            self.y
108        } else {
109            self.z
110        }
111    }
112
113    /// Returns the angle between "this" and v in radians, in the range [0, pi]. If
114    /// either vector is zero-length, or nearly zero-length, the result will be
115    /// zero, regardless of the other value.
116    pub fn angle(&self, b: &Self) -> f64 {
117        atan2(self.cross(b).norm(), self.dot(b))
118    }
119
120    /// Get the cross product of two XYZ Points
121    pub fn cross(&self, b: &Self) -> Self {
122        Self::new(
123            self.y * b.z - self.z * b.y,
124            self.z * b.x - self.x * b.z,
125            self.x * b.y - self.y * b.x,
126        )
127    }
128
129    /// Returns a Face-ST representation of this point
130    pub fn to_face_st(&self) -> (u8, f64, f64) {
131        xyz_to_face_st(self)
132    }
133
134    /// Returns the S2 face assocated with this point
135    pub fn get_face(&self) -> u8 {
136        xyz_to_face_uv(self).0
137    }
138
139    /// dot returns the standard dot product of v and ov.
140    pub fn dot(&self, b: &Self) -> f64 {
141        self.x * b.x + self.y * b.y + self.z * b.z
142    }
143
144    /// Returns the absolute value of the point.
145    pub fn abs(&self) -> Self {
146        Self::new(fabs(self.x), fabs(self.y), fabs(self.z))
147    }
148
149    /// Returns the inverse of the point
150    pub fn invert(&self) -> Self {
151        Self::new(-self.x, -self.y, -self.z)
152    }
153
154    /// Returns the length of the point.
155    pub fn len(&self) -> f64 {
156        self.norm()
157    }
158
159    /// Returns the vector's squared norm.
160    pub fn norm(&self) -> f64 {
161        sqrt(self.norm2())
162    }
163
164    /// The dot product of the point with itself
165    pub fn norm2(&self) -> f64 {
166        self.dot(self)
167    }
168
169    /// Normalize this point to unit length.
170    pub fn normalize(&mut self) {
171        let len = self.len();
172        if len > 0.0 {
173            self.x /= len;
174            self.y /= len;
175            self.z /= len;
176        }
177    }
178
179    /// return the distance from this point to the other point in radians
180    pub fn distance(&self, b: &Self) -> f64 {
181        let d = *self - *b;
182        d.len()
183    }
184
185    /// Returns the largest absolute component of the point.
186    pub fn largest_abs_component(&self) -> u8 {
187        let tmp = self.abs();
188        if tmp.x > tmp.y {
189            if tmp.x > tmp.z { 0 } else { 2 }
190        } else if tmp.y > tmp.z {
191            1
192        } else {
193            2
194        }
195    }
196
197    /// Returns the intermediate point between this and the other point.
198    pub fn intermediate(&self, b: &Self, t: f64) -> Self {
199        Self::new(
200            self.x + ((b.x - self.x) * (1.0 - t)),
201            self.y + ((b.y - self.y) * (1.0 - t)),
202            self.z + ((b.z - self.z) * (1.0 - t)),
203        )
204    }
205
206    /// Convert a u-v coordinate to an XYZ Point.
207    pub fn from_face_uv(face: u8, u: f64, v: f64) -> Self {
208        let mut p = face_uv_to_xyz(face, u, v);
209        p.normalize();
210        p
211    }
212
213    /// Convert an s-t coordinate to an XYZ Point.
214    pub fn from_face_st(face: u8, s: f64, t: f64) -> Self {
215        let u = ST_TO_UV(s);
216        let v = ST_TO_UV(t);
217        Self::from_face_uv(face, u, v)
218    }
219}
220impl<M: Clone + Default> From<&LonLat<M>> for S2Point {
221    fn from(lonlat: &LonLat<M>) -> Self {
222        lonlat.to_point()
223    }
224}
225impl<M: Clone + Default> From<&VectorPoint<M>> for S2Point {
226    fn from(v: &VectorPoint<M>) -> Self {
227        Self { x: v.x, y: v.y, z: v.z.unwrap_or(0.0) }
228    }
229}
230impl<M: Clone + Default> From<&mut VectorPoint<M>> for S2Point {
231    fn from(v: &mut VectorPoint<M>) -> Self {
232        Self { x: v.x, y: v.y, z: v.z.unwrap_or(0.0) }
233    }
234}
235impl From<S2CellId> for S2Point {
236    fn from(cellid: S2CellId) -> Self {
237        cellid.to_point()
238    }
239}
240// Implementing the Add trait for S2Point
241impl Add<S2Point> for S2Point {
242    type Output = Self;
243    fn add(self, other: Self) -> Self::Output {
244        S2Point { x: self.x + other.x, y: self.y + other.y, z: self.z + other.z }
245    }
246}
247impl AddAssign<S2Point> for S2Point {
248    fn add_assign(&mut self, other: S2Point) {
249        self.x += other.x;
250        self.y += other.y;
251        self.z += other.z;
252    }
253}
254impl Add<f64> for S2Point {
255    type Output = Self;
256    fn add(self, other: f64) -> Self::Output {
257        S2Point { x: self.x + other, y: self.y + other, z: self.z + other }
258    }
259}
260impl AddAssign<f64> for S2Point {
261    fn add_assign(&mut self, other: f64) {
262        self.x += other;
263        self.y += other;
264        self.z += other;
265    }
266}
267// Implementing the Sub trait for S2Point
268impl Sub<S2Point> for S2Point {
269    type Output = Self;
270    fn sub(self, other: Self) -> Self::Output {
271        S2Point { x: self.x - other.x, y: self.y - other.y, z: self.z - other.z }
272    }
273}
274impl SubAssign<S2Point> for S2Point {
275    fn sub_assign(&mut self, other: S2Point) {
276        self.x -= other.x;
277        self.y -= other.y;
278        self.z -= other.z;
279    }
280}
281impl Sub<f64> for S2Point {
282    type Output = Self;
283    fn sub(self, other: f64) -> Self::Output {
284        S2Point { x: self.x - other, y: self.y - other, z: self.z - other }
285    }
286}
287impl SubAssign<f64> for S2Point {
288    fn sub_assign(&mut self, other: f64) {
289        self.x -= other;
290        self.y -= other;
291        self.z -= other;
292    }
293}
294// Implementing the Neg trait for S2Point
295impl Neg for S2Point {
296    type Output = Self;
297    fn neg(self) -> Self::Output {
298        S2Point { x: -self.x, y: -self.y, z: -self.z }
299    }
300}
301// Implementing the Div trait for S2Point
302impl Div<S2Point> for S2Point {
303    type Output = Self;
304    fn div(self, other: Self) -> Self::Output {
305        S2Point { x: self.x / other.x, y: self.y / other.y, z: self.z / other.z }
306    }
307}
308impl DivAssign<S2Point> for S2Point {
309    fn div_assign(&mut self, other: S2Point) {
310        self.x /= other.x;
311        self.y /= other.y;
312        self.z /= other.z;
313    }
314}
315impl Div<f64> for S2Point {
316    type Output = Self;
317    fn div(self, other: f64) -> Self::Output {
318        S2Point { x: self.x / other, y: self.y / other, z: self.z / other }
319    }
320}
321impl DivAssign<f64> for S2Point {
322    fn div_assign(&mut self, other: f64) {
323        self.x /= other;
324        self.y /= other;
325        self.z /= other;
326    }
327}
328// Implementing the Mul trait for S2Point
329impl Mul<S2Point> for S2Point {
330    type Output = Self;
331    fn mul(self, other: Self) -> Self::Output {
332        S2Point { x: self.x * other.x, y: self.y * other.y, z: self.z * other.z }
333    }
334}
335impl MulAssign<S2Point> for S2Point {
336    fn mul_assign(&mut self, other: S2Point) {
337        self.x *= other.x;
338        self.y *= other.y;
339        self.z *= other.z;
340    }
341}
342impl Mul<f64> for S2Point {
343    type Output = Self;
344    fn mul(self, other: f64) -> Self::Output {
345        S2Point { x: self.x * other, y: self.y * other, z: self.z * other }
346    }
347}
348impl MulAssign<f64> for S2Point {
349    fn mul_assign(&mut self, other: f64) {
350        self.x *= other;
351        self.y *= other;
352        self.z *= other;
353    }
354}
355impl Rem<f64> for S2Point {
356    type Output = Self;
357    fn rem(self, other: f64) -> Self::Output {
358        S2Point { x: self.x % other, y: self.y % other, z: self.z % other }
359    }
360}
361impl RemAssign<f64> for S2Point {
362    fn rem_assign(&mut self, other: f64) {
363        self.x %= other;
364        self.y %= other;
365        self.z %= other;
366    }
367}
368impl Eq for S2Point {}
369impl Ord for S2Point {
370    fn cmp(&self, other: &Self) -> Ordering {
371        match self.x.partial_cmp(&other.x) {
372            Some(Ordering::Equal) => {}
373            other => return other.unwrap(), // Handle cases where `x` comparison is not equal
374        }
375        match self.y.partial_cmp(&other.y) {
376            Some(Ordering::Equal) => {}
377            other => return other.unwrap(), // Handle cases where `y` comparison is not equal
378        }
379        match self.z.partial_cmp(&other.z) {
380            Some(order) => order,
381            None => Ordering::Equal, // This handles the NaN case safely
382        }
383    }
384}
385impl PartialOrd for S2Point {
386    fn partial_cmp(&self, other: &Self) -> Option<Ordering> {
387        Some(self.cmp(other))
388    }
389}
gistools/geometry/s2/point.rs

gistools/geometry/s2/
point.rs
