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

s2json_core/geometry/
vector_point.rs
