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