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

s2json_core/geometry/
vector_point.rs
