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 { 0 } else { 2 }
224        } else if tmp.y > tmp_z {
225            1
226        } else {
227            2
228        }
229    }
230
231    /// Returns the intermediate point between this and the other point.
232    pub fn intermediate<M2: Clone>(&self, b: &VectorPoint<M2>, t: f64) -> Self {
233        if let (Some(z), Some(bz)) = (self.z, b.z) {
234            Self::new_xyz(
235                self.x + ((b.x - self.x) * (1.0 - t)),
236                self.y + ((b.y - self.y) * (1.0 - t)),
237                z + ((bz - z) * (1.0 - t)),
238                None,
239            )
240        } else {
241            Self::new_xy(
242                self.x + ((b.x - self.x) * (1.0 - t)),
243                self.y + ((b.y - self.y) * (1.0 - t)),
244                None,
245            )
246        }
247    }
248
249    /// Returns the perpendicular vector
250    pub fn perpendicular(&self) -> Self {
251        if let Some(z) = self.z {
252            let ref_point = if fabs(self.x) > fabs(z) {
253                Self::new_xyz(0., 0., 1., None)
254            } else {
255                Self::new_xyz(1., 0., 0., None)
256            };
257            let cross = self.cross(&ref_point);
258            Self::new_xyz(-self.y, self.x, cross.z.unwrap(), None)
259        } else {
260            Self::new_xy(-self.y, self.x, None)
261        }
262    }
263}
264impl<M: Clone> From<Point> for VectorPoint<M> {
265    fn from(p: Point) -> Self {
266        Self { x: p.0, y: p.1, z: None, m: None, t: None }
267    }
268}
269impl<M: Clone> From<&Point> for VectorPoint<M> {
270    fn from(p: &Point) -> Self {
271        Self { x: p.0, y: p.1, z: None, m: None, t: None }
272    }
273}
274impl<M: Clone> From<Point3D> for VectorPoint<M> {
275    fn from(p: Point3D) -> Self {
276        Self { x: p.0, y: p.1, z: Some(p.2), m: None, t: None }
277    }
278}
279impl<M: Clone> From<&Point3D> for VectorPoint<M> {
280    fn from(p: &Point3D) -> Self {
281        Self { x: p.0, y: p.1, z: Some(p.2), m: None, t: None }
282    }
283}
284impl<M1: Clone, M2: Clone> Add<&VectorPoint<M2>> for &VectorPoint<M1> {
285    type Output = VectorPoint<M1>;
286    fn add(self, other: &VectorPoint<M2>) -> Self::Output {
287        VectorPoint {
288            x: self.x + other.x,
289            y: self.y + other.y,
290            // Only add `z` if both `self.z` and `other.z` are `Some`
291            z: match (self.z, other.z) {
292                (Some(z1), Some(z2)) => Some(z1 + z2),
293                _ => None, // If either `z` is None, the result is `None`
294            },
295            m: None, // Combine m values
296            t: None, // Handle `t` as necessary
297        }
298    }
299}
300impl<M1: Clone, M2: Clone> AddAssign<&VectorPoint<M2>> for VectorPoint<M1> {
301    fn add_assign(&mut self, other: &VectorPoint<M2>) {
302        self.x += other.x;
303        self.y += other.y;
304        if let (Some(z), Some(other_z)) = (self.z, other.z) {
305            self.z = Some(z + other_z);
306        }
307    }
308}
309impl<M: Clone> Add<f64> for &VectorPoint<M> {
310    type Output = VectorPoint<M>;
311    fn add(self, other: f64) -> Self::Output {
312        VectorPoint {
313            x: self.x + other,
314            y: self.y + other,
315            z: self.z.map(|z| z + other),
316            m: None,
317            t: None,
318        }
319    }
320}
321impl<M: Clone> AddAssign<f64> for VectorPoint<M> {
322    fn add_assign(&mut self, other: f64) {
323        self.x += other;
324        self.y += other;
325        if let Some(z) = self.z {
326            self.z = Some(z + other);
327        }
328    }
329}
330// Implementing the Sub trait for VectorPoint
331impl<M1: Clone, M2: Clone> Sub<&VectorPoint<M2>> for &VectorPoint<M1> {
332    type Output = VectorPoint<M1>;
333    fn sub(self, other: &VectorPoint<M2>) -> Self::Output {
334        VectorPoint {
335            x: self.x - other.x,
336            y: self.y - other.y,
337            z: match (self.z, other.z) {
338                (Some(z1), Some(z2)) => Some(z1 - z2),
339                _ => None, // If either `z` is None, the result is `None`
340            },
341            m: None, // Combine m values
342            t: None, // Handle `t` as necessary
343        }
344    }
345}
346impl<M1: Clone, M2: Clone> SubAssign<&VectorPoint<M2>> for VectorPoint<M1> {
347    fn sub_assign(&mut self, other: &VectorPoint<M2>) {
348        self.x -= other.x;
349        self.y -= other.y;
350        if let (Some(z), Some(other_z)) = (self.z, other.z) {
351            self.z = Some(z - other_z);
352        }
353    }
354}
355impl<M: Clone> Sub<f64> for &VectorPoint<M> {
356    type Output = VectorPoint<M>;
357    fn sub(self, other: f64) -> Self::Output {
358        VectorPoint {
359            x: self.x - other,
360            y: self.y - other,
361            z: self.z.map(|z| z - other),
362            m: None,
363            t: None,
364        }
365    }
366}
367impl<M: Clone> SubAssign<f64> for VectorPoint<M> {
368    fn sub_assign(&mut self, other: f64) {
369        self.x -= other;
370        self.y -= other;
371        if let Some(z) = self.z {
372            self.z = Some(z - other);
373        }
374    }
375}
376// Implementing the Neg trait for VectorPoint
377impl<M: Clone> Neg for &VectorPoint<M> {
378    type Output = VectorPoint<M>;
379    fn neg(self) -> Self::Output {
380        VectorPoint { x: -self.x, y: -self.y, z: self.z.map(|z| -z), m: None, t: None }
381    }
382}
383// Implementing the Div trait for VectorPoint
384impl<M1: Clone, M2: Clone> Div<&VectorPoint<M2>> for &VectorPoint<M1> {
385    type Output = VectorPoint<M1>;
386    fn div(self, other: &VectorPoint<M2>) -> Self::Output {
387        VectorPoint {
388            x: self.x / other.x,
389            y: self.y / other.y,
390            z: match (self.z, other.z) {
391                (Some(z1), Some(z2)) => Some(z1 / z2),
392                _ => None, // If either `z` is None, the result is `None`
393            },
394            m: None, // Combine m values
395            t: None, // Handle `t` as necessary
396        }
397    }
398}
399impl<M1: Clone, M2: Clone> DivAssign<&VectorPoint<M2>> for VectorPoint<M1> {
400    fn div_assign(&mut self, other: &VectorPoint<M2>) {
401        self.x /= other.x;
402        self.y /= other.y;
403        if let (Some(z), Some(other_z)) = (self.z, other.z) {
404            self.z = Some(z / other_z);
405        }
406    }
407}
408impl<M: Clone> Div<f64> for &VectorPoint<M> {
409    type Output = VectorPoint<M>;
410    fn div(self, other: f64) -> Self::Output {
411        VectorPoint {
412            x: self.x / other,
413            y: self.y / other,
414            z: self.z.map(|z| z / other),
415            m: None,
416            t: None,
417        }
418    }
419}
420impl<M: Clone> DivAssign<f64> for VectorPoint<M> {
421    fn div_assign(&mut self, other: f64) {
422        self.x /= other;
423        self.y /= other;
424        if let Some(z) = self.z {
425            self.z = Some(z / other);
426        }
427    }
428}
429// Implementing the Mul trait for VectorPoint
430impl<M1: Clone, M2: Clone> Mul<&VectorPoint<M2>> for &VectorPoint<M1> {
431    type Output = VectorPoint<M1>;
432    fn mul(self, other: &VectorPoint<M2>) -> Self::Output {
433        VectorPoint {
434            x: self.x * other.x,
435            y: self.y * other.y,
436            z: match (self.z, other.z) {
437                (Some(z1), Some(z2)) => Some(z1 * z2),
438                _ => None, // If either `z` is None, the result is `None`
439            },
440            m: None, // Combine m values
441            t: None, // Handle `t` as necessary
442        }
443    }
444}
445impl<M1: Clone, M2: Clone> MulAssign<&VectorPoint<M2>> for VectorPoint<M1> {
446    fn mul_assign(&mut self, other: &VectorPoint<M2>) {
447        self.x *= other.x;
448        self.y *= other.y;
449        if let (Some(z), Some(other_z)) = (self.z, other.z) {
450            self.z = Some(z * other_z);
451        }
452    }
453}
454impl<M: Clone> Mul<f64> for &VectorPoint<M> {
455    type Output = VectorPoint<M>;
456    fn mul(self, other: f64) -> Self::Output {
457        VectorPoint {
458            x: self.x * other,
459            y: self.y * other,
460            z: self.z.map(|z| z * other),
461            m: None,
462            t: None,
463        }
464    }
465}
466impl<M: Clone> MulAssign<f64> for VectorPoint<M> {
467    fn mul_assign(&mut self, other: f64) {
468        self.x *= other;
469        self.y *= other;
470        if let Some(z) = self.z {
471            self.z = Some(z * other);
472        }
473    }
474}
475impl<M: Clone> Rem<f64> for VectorPoint<M> {
476    type Output = VectorPoint<M>;
477
478    fn rem(self, modulus: f64) -> Self::Output {
479        self.modulo(modulus)
480    }
481}
482impl<M: Clone> RemAssign<f64> for VectorPoint<M> {
483    fn rem_assign(&mut self, modulus: f64) {
484        let modulus = fabs(modulus);
485        self.x = fmod(self.x, modulus);
486        self.y = fmod(self.y, modulus);
487        if let Some(z) = self.z {
488            self.z = Some(fmod(z, modulus));
489        }
490    }
491}
492impl<M: Clone> Eq for VectorPoint<M> {}
493impl<M: Clone> Ord for VectorPoint<M> {
494    fn cmp(&self, other: &VectorPoint<M>) -> Ordering {
495        match self.x.partial_cmp(&other.x) {
496            Some(Ordering::Equal) => {}
497            other => return other.unwrap_or(Ordering::Greater), /* Handle cases where `x` comparison is not equal */
498        }
499        match self.y.partial_cmp(&other.y) {
500            Some(Ordering::Equal) => {}
501            other => return other.unwrap_or(Ordering::Greater), /* Handle cases where `y` comparison is not equal */
502        }
503        match self.z.partial_cmp(&other.z) {
504            Some(order) => order,
505            None => Ordering::Equal, // This handles the NaN case safely
506        }
507    }
508}
509impl<M: Clone> PartialEq for VectorPoint<M> {
510    fn eq(&self, other: &VectorPoint<M>) -> bool {
511        self.x == other.x && self.y == other.y && self.z == other.z
512    }
513}
514impl<M: Clone> PartialOrd for VectorPoint<M> {
515    fn partial_cmp(&self, other: &VectorPoint<M>) -> Option<Ordering> {
516        Some(self.cmp(other))
517    }
518}