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