s2json_core/geometry/
vector_point.rs

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