s2json_core/geometry/
vector_point.rs

1use crate::*;
2use core::{
3    cmp::Ordering,
4    f64::consts::PI,
5    fmt::Debug,
6    ops::{Add, AddAssign, Div, DivAssign, Mul, MulAssign, Neg, Rem, RemAssign, Sub, SubAssign},
7};
8use libm::{atan, atan2, fabs, fmod, log, sin, sinh, sqrt};
9use serde::{Deserialize, Serialize};
10
11/// # Vector Points
12///
13/// ## Description
14/// A Vector Point uses a structure for 2D or 3D points. Useful for geometry that also has m-values
15///
16/// Implements the ops [`Add`], [`AddAssign`], [`Sub`], [`SubAssign`], [`Mul`], [`MulAssign`], [`Div`], [`DivAssign`], [`Rem`], [`RemAssign`], [`Neg`], [`PartialEq`], and [`Eq`].
17///
18/// Also implements local traits [`GetXY`], [`GetZ`], [`GetM`], [`SetXY`], [`SetZ`], [`SetM`], [`NewXY`], [`NewXYM`], [`NewXYZ`], [`NewXYZM`], [`MValueCompatible`], [`MValue`], and [`Interpolate`].
19///
20/// ## Fields
21/// - `x`: X coordinate
22/// - `y`: Y coordinate
23/// - `z`: Z coordinate or "altitude". May be None
24/// - `m`: M-Value
25/// - `t`: T for tolerance. A tmp value often used for simplification if needed.
26///
27/// ## Usage
28///
29/// - [`VectorPoint::new`]: Creates a new VectorPoint
30/// - [`VectorPoint::to_m_value`]: Convert to an MValue VectorPoint using to_m_value
31/// - [`VectorPoint::project`]: Project the point into the 0->1 coordinate system
32/// - [`VectorPoint::unproject`]: Unproject the point from the 0->1 coordinate system
33/// - [`VectorPoint::is_empty`]: Returns true if the point is the zero vector
34/// - [`VectorPoint::face`]: Returns the S2 face assocated with this point
35/// - [`VectorPoint::modulo]: Apply modular arithmetic
36/// - [`VectorPoint::angle`]: 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
37/// - [`VectorPoint::cross`]: Get the cross product of two Vector Points
38/// - [`VectorPoint::dot`]: dot returns the standard dot product of v and ov.
39/// - [`VectorPoint::abs`]: Returns the absolute value of the point.
40/// - [`VectorPoint::invert`]: Returns the inverse of the point
41/// - [`VectorPoint::len`]: Returns the length of the point
42/// - [`VectorPoint::norm`]: norm returns the vector's norm. (sqrt of norm2)
43/// - [`VectorPoint::norm2`]: norm2 returns the vector's squared norm (self.dot(self))
44/// - [`VectorPoint::normalize`]: Normalize this point to unit length
45/// - [`VectorPoint::distance`]: return the distance from this point to the other point
46/// - [`VectorPoint::largest_abs_component`]: Returns the largest absolute component of the point
47/// - [`VectorPoint::interpolate`]: Interpolate between two points with a given t value
48#[derive(Serialize, Deserialize, Debug, Clone, Default)]
49#[repr(C)]
50pub struct VectorPoint<M: Clone = MValue> {
51    /// X coordinate
52    pub x: f64,
53    /// Y coordinate
54    pub y: f64,
55    /// Z coordinate or "altitude". May be None
56    #[serde(skip_serializing_if = "Option::is_none")]
57    pub z: Option<f64>,
58    /// M-Value
59    #[serde(skip_serializing_if = "Option::is_none")]
60    pub m: Option<M>,
61    /// T for tolerance. A tmp value used for simplification
62    #[serde(skip)]
63    pub t: Option<f64>,
64}
65// GET
66impl<M: Clone> GetXY for VectorPoint<M> {
67    fn x(&self) -> f64 {
68        self.x
69    }
70    fn y(&self) -> f64 {
71        self.y
72    }
73}
74impl<M: Clone> GetZ for VectorPoint<M> {
75    fn z(&self) -> Option<f64> {
76        self.z
77    }
78}
79impl<M: Clone> GetM<M> for VectorPoint<M> {
80    fn m(&self) -> Option<&M> {
81        self.m.as_ref()
82    }
83}
84// SET
85impl<M: Clone> SetXY for VectorPoint<M> {
86    fn set_xy(&mut self, x: f64, y: f64) {
87        self.x = x;
88        self.y = y;
89    }
90
91    fn set_x(&mut self, x: f64) {
92        self.x = x;
93    }
94
95    fn set_y(&mut self, y: f64) {
96        self.y = y;
97    }
98}
99impl<M: Clone> SetZ for VectorPoint<M> {
100    fn set_z(&mut self, z: f64) {
101        self.z = Some(z);
102    }
103}
104impl<M: Clone> SetM<M> for VectorPoint<M> {
105    fn set_m(&mut self, m: M) {
106        self.m = Some(m);
107    }
108}
109// NEW
110impl<M: Clone> NewXY for VectorPoint<M> {
111    fn new_xy(x: f64, y: f64) -> Self {
112        Self { x, y, z: None, m: None, t: None }
113    }
114}
115impl<M: Clone> NewXYM<M> for VectorPoint<M> {
116    fn new_xym(x: f64, y: f64, m: M) -> Self {
117        Self { x, y, z: None, m: Some(m), t: None }
118    }
119}
120impl<M: Clone> NewXYZ for VectorPoint<M> {
121    fn new_xyz(x: f64, y: f64, z: f64) -> Self {
122        Self { x, y, z: Some(z), m: None, t: None }
123    }
124}
125impl<M: Clone> NewXYZM<M> for VectorPoint<M> {
126    fn new_xyzm(x: f64, y: f64, z: f64, m: M) -> Self {
127        Self { x, y, z: Some(z), m: Some(m), t: None }
128    }
129}
130// MValue METHODS
131impl VectorPoint<MValue> {
132    /// Helper function for tests. Create a new point quickly from an xy coordinate
133    pub fn from_xy(x: f64, y: f64) -> Self {
134        Self { x, y, z: None, m: None, t: None }
135    }
136
137    /// Helper function for tests. Create a new point quickly from an xyz coordinate
138    pub fn from_xyz(x: f64, y: f64, z: f64) -> Self {
139        Self { x, y, z: Some(z), m: None, t: None }
140    }
141}
142// M METHODS
143impl<M: Clone> VectorPoint<M> {
144    /// Create a new point
145    pub fn new(x: f64, y: f64, z: Option<f64>, m: Option<M>) -> Self {
146        Self { x, y, z, m, t: None }
147    }
148
149    /// Create a new point with xy
150    pub fn new_xy(x: f64, y: f64, m: Option<M>) -> Self {
151        Self { x, y, z: None, m, t: None }
152    }
153
154    /// Create a new point with xyz
155    pub fn new_xyz(x: f64, y: f64, z: f64, m: Option<M>) -> Self {
156        Self { x, y, z: Some(z), m, t: None }
157    }
158
159    /// Convert to an MValue VectorPoint
160    pub fn to_m_value(&self) -> VectorPoint<MValue>
161    where
162        M: Into<MValue>,
163    {
164        VectorPoint {
165            x: self.x,
166            y: self.y,
167            z: self.z,
168            m: self.m.clone().map(|m| m.into()),
169            t: None,
170        }
171    }
172
173    /// Project the point into the 0->1 coordinate system
174    pub fn project(&mut self, bbox: Option<&mut BBox3D>) {
175        let y = self.y;
176        let x = self.x;
177        let sin = sin((y * PI) / 180.);
178        let y2 = 0.5 - (0.25 * log((1. + sin) / (1. - sin))) / PI;
179        self.x = x / 360. + 0.5;
180        self.y = y2.clamp(0., 1.);
181
182        if let Some(bbox) = bbox {
183            bbox.extend_from_point(self)
184        };
185    }
186
187    /// Unproject the point from the 0->1 coordinate system back to a lon-lat coordinate
188    pub fn unproject(&mut self) {
189        let lon = (self.x - 0.5) * 360.;
190        let y2 = 0.5 - self.y;
191        let lat = atan(sinh(PI * (y2 * 2.))).to_degrees();
192
193        self.x = lon;
194        self.y = lat;
195    }
196
197    /// Returns true if the point is the zero vector.
198    pub fn is_empty(&self) -> bool {
199        let zero = f64::default();
200        self.x == zero && self.y == zero && (self.z.is_none() || self.z.unwrap() == zero)
201    }
202
203    /// Returns the S2 face assocated with this point
204    pub fn face(&self, f: u8) -> f64 {
205        if f == 0 {
206            self.x
207        } else if f == 1 {
208            self.y
209        } else {
210            self.z.unwrap_or(0.)
211        }
212    }
213
214    /// Apply modular arithmetic to x, y, and z using `modulus`
215    pub fn modulo(self, modulus: f64) -> Self {
216        let modulus = fabs(modulus); // Ensure positive modulus
217        Self {
218            x: fmod(self.x, modulus),
219            y: fmod(self.y, modulus),
220            z: self.z.map(|z| fmod(z, modulus)),
221            m: self.m,
222            t: None,
223        }
224    }
225
226    /// Returns the angle between "this" and v in radians, in the range [0, pi]. If
227    /// either vector is zero-length, or nearly zero-length, the result will be
228    /// zero, regardless of the other value.
229    pub fn angle<M2: Clone>(&self, b: &VectorPoint<M2>) -> f64 {
230        atan2(self.cross(b).norm(), self.dot(b))
231    }
232
233    /// Get the cross product of two Vector Points
234    pub fn cross<M2: Clone>(&self, b: &VectorPoint<M2>) -> Self {
235        if let (Some(z), Some(bz)) = (self.z, b.z) {
236            Self::new_xyz(
237                self.y * bz - z * b.y,
238                z * b.x - self.x * bz,
239                self.x * b.y - self.y * b.x,
240                None,
241            )
242        } else {
243            Self::new_xy(self.x * b.y - self.y * b.x, self.y * b.x - self.x * b.y, None)
244        }
245    }
246
247    /// dot returns the standard dot product of v and ov.
248    pub fn dot<M2: Clone>(&self, b: &VectorPoint<M2>) -> f64 {
249        if let (Some(z), Some(bz)) = (self.z, b.z) {
250            self.x * b.x + self.y * b.y + z * bz
251        } else {
252            self.x * b.x + self.y * b.y
253        }
254    }
255
256    /// Returns the absolute value of the point.
257    pub fn abs(&self) -> Self {
258        Self::new(fabs(self.x), fabs(self.y), self.z.map(fabs), None)
259    }
260
261    /// Returns the inverse of the point
262    pub fn invert(&self) -> Self {
263        -self
264    }
265
266    /// Returns the length of the point.
267    pub fn len(&self) -> f64 {
268        self.norm()
269    }
270
271    /// norm returns the vector's norm.
272    pub fn norm(&self) -> f64 {
273        sqrt(self.norm2())
274    }
275
276    /// norm2 returns the vector's squared norm.
277    pub fn norm2(&self) -> f64 {
278        self.dot(self)
279    }
280
281    /// Normalize this point to unit length.
282    pub fn normalize(&mut self) {
283        let len = self.len();
284        if len > 0.0 {
285            self.x /= len;
286            self.y /= len;
287            if let Some(z) = self.z {
288                self.z = Some(z / len);
289            }
290        }
291    }
292
293    /// return the distance from this point to the other point in radians
294    pub fn distance<M2: Clone>(&self, b: &VectorPoint<M2>) -> f64 {
295        let d: VectorPoint<M> = self - b;
296        d.len()
297    }
298
299    /// Returns the largest absolute component of the point.
300    pub fn largest_abs_component(&self) -> u8 {
301        let tmp = self.abs();
302        let tmp_z = tmp.z.unwrap_or(-2.);
303        if tmp.x > tmp.y {
304            if tmp.x > tmp_z { 0 } else { 2 }
305        } else if tmp.y > tmp_z {
306            1
307        } else {
308            2
309        }
310    }
311
312    /// Returns the intermediate point between this and the other point.
313    pub fn intermediate<M2: Clone>(&self, b: &VectorPoint<M2>, t: f64) -> Self {
314        if let (Some(z), Some(bz)) = (self.z, b.z) {
315            Self::new_xyz(
316                self.x + ((b.x - self.x) * (1.0 - t)),
317                self.y + ((b.y - self.y) * (1.0 - t)),
318                z + ((bz - z) * (1.0 - t)),
319                None,
320            )
321        } else {
322            Self::new_xy(
323                self.x + ((b.x - self.x) * (1.0 - t)),
324                self.y + ((b.y - self.y) * (1.0 - t)),
325                None,
326            )
327        }
328    }
329
330    /// Returns the perpendicular vector
331    pub fn perpendicular(&self) -> Self {
332        if let Some(z) = self.z {
333            let ref_point = if fabs(self.x) > fabs(z) {
334                Self::new_xyz(0., 0., 1., None)
335            } else {
336                Self::new_xyz(1., 0., 0., None)
337            };
338            let cross = self.cross(&ref_point);
339            Self::new_xyz(-self.y, self.x, cross.z.unwrap(), None)
340        } else {
341            Self::new_xy(-self.y, self.x, None)
342        }
343    }
344}
345impl<M: Clone, P: GetXYZ> From<&P> for VectorPoint<M> {
346    fn from(p: &P) -> Self {
347        Self { x: p.x(), y: p.y(), z: p.z(), m: None, t: None }
348    }
349}
350impl<M: Clone> From<Point> for VectorPoint<M> {
351    fn from(p: Point) -> Self {
352        Self { x: p.0, y: p.1, z: None, m: None, t: None }
353    }
354}
355impl<M: Clone> From<Point3D> for VectorPoint<M> {
356    fn from(p: Point3D) -> Self {
357        Self { x: p.0, y: p.1, z: Some(p.2), m: None, t: None }
358    }
359}
360impl<M1: Clone, M2: Clone> Add<&VectorPoint<M2>> for &VectorPoint<M1> {
361    type Output = VectorPoint<M1>;
362    fn add(self, other: &VectorPoint<M2>) -> Self::Output {
363        VectorPoint {
364            x: self.x + other.x,
365            y: self.y + other.y,
366            // Only add `z` if both `self.z` and `other.z` are `Some`
367            z: match (self.z, other.z) {
368                (Some(z1), Some(z2)) => Some(z1 + z2),
369                _ => None, // If either `z` is None, the result is `None`
370            },
371            m: None, // Combine m values
372            t: None, // Handle `t` as necessary
373        }
374    }
375}
376impl<M1: Clone, M2: Clone> AddAssign<&VectorPoint<M2>> for VectorPoint<M1> {
377    fn add_assign(&mut self, other: &VectorPoint<M2>) {
378        self.x += other.x;
379        self.y += other.y;
380        if let (Some(z), Some(other_z)) = (self.z, other.z) {
381            self.z = Some(z + other_z);
382        }
383    }
384}
385impl<M: Clone> Add<f64> for &VectorPoint<M> {
386    type Output = VectorPoint<M>;
387    fn add(self, other: f64) -> Self::Output {
388        VectorPoint {
389            x: self.x + other,
390            y: self.y + other,
391            z: self.z.map(|z| z + other),
392            m: None,
393            t: None,
394        }
395    }
396}
397impl<M: Clone> AddAssign<f64> for VectorPoint<M> {
398    fn add_assign(&mut self, other: f64) {
399        self.x += other;
400        self.y += other;
401        if let Some(z) = self.z {
402            self.z = Some(z + other);
403        }
404    }
405}
406// Implementing the Sub trait for VectorPoint
407impl<M1: Clone, M2: Clone> Sub<&VectorPoint<M2>> for &VectorPoint<M1> {
408    type Output = VectorPoint<M1>;
409    fn sub(self, other: &VectorPoint<M2>) -> Self::Output {
410        VectorPoint {
411            x: self.x - other.x,
412            y: self.y - other.y,
413            z: match (self.z, other.z) {
414                (Some(z1), Some(z2)) => Some(z1 - z2),
415                _ => None, // If either `z` is None, the result is `None`
416            },
417            m: None, // Combine m values
418            t: None, // Handle `t` as necessary
419        }
420    }
421}
422impl<M1: Clone, M2: Clone> SubAssign<&VectorPoint<M2>> for VectorPoint<M1> {
423    fn sub_assign(&mut self, other: &VectorPoint<M2>) {
424        self.x -= other.x;
425        self.y -= other.y;
426        if let (Some(z), Some(other_z)) = (self.z, other.z) {
427            self.z = Some(z - other_z);
428        }
429    }
430}
431impl<M: Clone> Sub<f64> for &VectorPoint<M> {
432    type Output = VectorPoint<M>;
433    fn sub(self, other: f64) -> Self::Output {
434        VectorPoint {
435            x: self.x - other,
436            y: self.y - other,
437            z: self.z.map(|z| z - other),
438            m: None,
439            t: None,
440        }
441    }
442}
443impl<M: Clone> SubAssign<f64> for VectorPoint<M> {
444    fn sub_assign(&mut self, other: f64) {
445        self.x -= other;
446        self.y -= other;
447        if let Some(z) = self.z {
448            self.z = Some(z - other);
449        }
450    }
451}
452// Implementing the Neg trait for VectorPoint
453impl<M: Clone> Neg for &VectorPoint<M> {
454    type Output = VectorPoint<M>;
455    fn neg(self) -> Self::Output {
456        VectorPoint { x: -self.x, y: -self.y, z: self.z.map(|z| -z), m: None, t: None }
457    }
458}
459// Implementing the Div trait for VectorPoint
460impl<M1: Clone, M2: Clone> Div<&VectorPoint<M2>> for &VectorPoint<M1> {
461    type Output = VectorPoint<M1>;
462    fn div(self, other: &VectorPoint<M2>) -> Self::Output {
463        VectorPoint {
464            x: self.x / other.x,
465            y: self.y / other.y,
466            z: match (self.z, other.z) {
467                (Some(z1), Some(z2)) => Some(z1 / z2),
468                _ => None, // If either `z` is None, the result is `None`
469            },
470            m: None, // Combine m values
471            t: None, // Handle `t` as necessary
472        }
473    }
474}
475impl<M1: Clone, M2: Clone> DivAssign<&VectorPoint<M2>> for VectorPoint<M1> {
476    fn div_assign(&mut self, other: &VectorPoint<M2>) {
477        self.x /= other.x;
478        self.y /= other.y;
479        if let (Some(z), Some(other_z)) = (self.z, other.z) {
480            self.z = Some(z / other_z);
481        }
482    }
483}
484impl<M: Clone> Div<f64> for &VectorPoint<M> {
485    type Output = VectorPoint<M>;
486    fn div(self, other: f64) -> Self::Output {
487        VectorPoint {
488            x: self.x / other,
489            y: self.y / other,
490            z: self.z.map(|z| z / other),
491            m: None,
492            t: None,
493        }
494    }
495}
496impl<M: Clone> DivAssign<f64> for VectorPoint<M> {
497    fn div_assign(&mut self, other: f64) {
498        self.x /= other;
499        self.y /= other;
500        if let Some(z) = self.z {
501            self.z = Some(z / other);
502        }
503    }
504}
505// Implementing the Mul trait for VectorPoint
506impl<M1: Clone, M2: Clone> Mul<&VectorPoint<M2>> for &VectorPoint<M1> {
507    type Output = VectorPoint<M1>;
508    fn mul(self, other: &VectorPoint<M2>) -> Self::Output {
509        VectorPoint {
510            x: self.x * other.x,
511            y: self.y * other.y,
512            z: match (self.z, other.z) {
513                (Some(z1), Some(z2)) => Some(z1 * z2),
514                _ => None, // If either `z` is None, the result is `None`
515            },
516            m: None, // Combine m values
517            t: None, // Handle `t` as necessary
518        }
519    }
520}
521impl<M1: Clone, M2: Clone> MulAssign<&VectorPoint<M2>> for VectorPoint<M1> {
522    fn mul_assign(&mut self, other: &VectorPoint<M2>) {
523        self.x *= other.x;
524        self.y *= other.y;
525        if let (Some(z), Some(other_z)) = (self.z, other.z) {
526            self.z = Some(z * other_z);
527        }
528    }
529}
530impl<M: Clone> Mul<f64> for &VectorPoint<M> {
531    type Output = VectorPoint<M>;
532    fn mul(self, other: f64) -> Self::Output {
533        VectorPoint {
534            x: self.x * other,
535            y: self.y * other,
536            z: self.z.map(|z| z * other),
537            m: None,
538            t: None,
539        }
540    }
541}
542impl<M: Clone> MulAssign<f64> for VectorPoint<M> {
543    fn mul_assign(&mut self, other: f64) {
544        self.x *= other;
545        self.y *= other;
546        if let Some(z) = self.z {
547            self.z = Some(z * other);
548        }
549    }
550}
551impl<M: Clone> Rem<f64> for VectorPoint<M> {
552    type Output = VectorPoint<M>;
553
554    fn rem(self, modulus: f64) -> Self::Output {
555        self.modulo(modulus)
556    }
557}
558impl<M: Clone> RemAssign<f64> for VectorPoint<M> {
559    fn rem_assign(&mut self, modulus: f64) {
560        let modulus = fabs(modulus);
561        self.x = fmod(self.x, modulus);
562        self.y = fmod(self.y, modulus);
563        if let Some(z) = self.z {
564            self.z = Some(fmod(z, modulus));
565        }
566    }
567}
568impl<M: Clone> Eq for VectorPoint<M> {}
569impl<M: Clone> Ord for VectorPoint<M> {
570    fn cmp(&self, other: &VectorPoint<M>) -> Ordering {
571        match self.x.partial_cmp(&other.x) {
572            Some(Ordering::Equal) => {}
573            other => return other.unwrap_or(Ordering::Greater), /* Handle cases where `x` comparison is not equal */
574        }
575        match self.y.partial_cmp(&other.y) {
576            Some(Ordering::Equal) => {}
577            other => return other.unwrap_or(Ordering::Greater), /* Handle cases where `y` comparison is not equal */
578        }
579        match self.z.partial_cmp(&other.z) {
580            Some(order) => order,
581            None => Ordering::Equal, // This handles the NaN case safely
582        }
583    }
584}
585impl<M: Clone> PartialEq for VectorPoint<M> {
586    fn eq(&self, other: &VectorPoint<M>) -> bool {
587        self.x == other.x && self.y == other.y && self.z == other.z
588    }
589}
590impl<M: Clone> PartialOrd for VectorPoint<M> {
591    fn partial_cmp(&self, other: &VectorPoint<M>) -> Option<Ordering> {
592        Some(self.cmp(other))
593    }
594}
s2json_core/geometry/vector_point.rs

s2json_core/geometry/
vector_point.rs
