ogeom_math/
point.rs

1//! Positions in 2D and 3D.
2//!
3//! A point is affected by translation; a [`Vector`] is not. Keeping them
4//! distinct types means the compiler rejects the classic errors (translating a
5//! normal, adding two positions) rather than letting them produce plausible
6//! nonsense.
7//!
8//! The algebra that results is the affine one: point − point is a vector, point
9//! + vector is a point, and point + point is not defined.
10
11use core::ops::{Add, AddAssign, Sub, SubAssign};
12
13use ogeom_core::{OgeomResult, Tolerances, ogeom_bail};
14
15use crate::{Vector, Vector2};
16
17/// A position in space.
18#[derive(Debug, Clone, Copy, PartialEq, Default)]
19pub struct Point {
20    /// X coordinate.
21    pub x: f64,
22    /// Y coordinate.
23    pub y: f64,
24    /// Z coordinate.
25    pub z: f64,
26}
27
28/// A position in the plane.
29#[derive(Debug, Clone, Copy, PartialEq, Default)]
30pub struct Point2 {
31    /// X coordinate.
32    pub x: f64,
33    /// Y coordinate.
34    pub y: f64,
35}
36
37impl Point {
38    /// The origin.
39    pub const ORIGIN: Self = Self::new(0.0, 0.0, 0.0);
40
41    /// A point from coordinates.
42    #[must_use]
43    pub const fn new(x: f64, y: f64, z: f64) -> Self {
44        Self { x, y, z }
45    }
46
47    /// Coordinates as an array.
48    #[must_use]
49    pub const fn to_array(self) -> [f64; 3] {
50        [self.x, self.y, self.z]
51    }
52
53    /// A point from an array.
54    #[must_use]
55    pub const fn from_array([x, y, z]: [f64; 3]) -> Self {
56        Self::new(x, y, z)
57    }
58
59    /// The position vector from the origin.
60    #[must_use]
61    pub const fn to_vector(self) -> Vector {
62        Vector::new(self.x, self.y, self.z)
63    }
64
65    /// The point at the tip of `v` placed at the origin.
66    #[must_use]
67    pub const fn from_vector(v: Vector) -> Self {
68        Self::new(v.x, v.y, v.z)
69    }
70
71    /// Coordinate by index, `0..3`.
72    ///
73    /// # Errors
74    ///
75    /// [`OgeomError::Range`](ogeom_core::OgeomError::Range) if `index >= 3`.
76    pub fn coord(self, index: usize) -> OgeomResult<f64> {
77        match index {
78            0 => Ok(self.x),
79            1 => Ok(self.y),
80            2 => Ok(self.z),
81            _ => ogeom_bail!(Range, "point coordinate {index} of 3"),
82        }
83    }
84
85    /// Squared distance to `other`. Prefer this when comparing distances.
86    #[must_use]
87    pub fn square_distance(self, other: Self) -> f64 {
88        (other - self).square_magnitude()
89    }
90
91    /// Distance to `other`.
92    #[must_use]
93    pub fn distance(self, other: Self) -> f64 {
94        self.square_distance(other).sqrt()
95    }
96
97    /// Whether two points coincide within `tol.confusion()`.
98    #[must_use]
99    pub fn is_equal(self, other: Self, tol: Tolerances) -> bool {
100        self.square_distance(other) <= tol.confusion() * tol.confusion()
101    }
102
103    /// Whether two points coincide within an explicit distance.
104    #[must_use]
105    pub fn is_within(self, other: Self, distance: f64) -> bool {
106        self.square_distance(other) <= distance * distance
107    }
108
109    /// Whether every coordinate is finite.
110    #[must_use]
111    pub fn is_finite(self) -> bool {
112        self.x.is_finite() && self.y.is_finite() && self.z.is_finite()
113    }
114
115    /// The midpoint of `self` and `other`.
116    #[must_use]
117    pub fn midpoint(self, other: Self) -> Self {
118        Self::new(
119            f64::midpoint(self.x, other.x),
120            f64::midpoint(self.y, other.y),
121            f64::midpoint(self.z, other.z),
122        )
123    }
124
125    /// Linear interpolation, `t = 0` giving `self`.
126    #[must_use]
127    pub fn lerp(self, other: Self, t: f64) -> Self {
128        self + (other - self) * t
129    }
130
131    /// The centroid of a set of points.
132    ///
133    /// # Errors
134    ///
135    /// [`OgeomError::Construction`](ogeom_core::OgeomError::Construction) if the set is
136    /// empty.
137    pub fn centroid(points: &[Self]) -> OgeomResult<Self> {
138        let Some((first, rest)) = points.split_first() else {
139            ogeom_bail!(Construction, "centroid of an empty point set");
140        };
141        // Accumulate offsets from the first point rather than absolute
142        // coordinates: for a cluster far from the origin the absolute sum loses
143        // precision to cancellation, while the offsets stay small.
144        let mut sum = Vector::ZERO;
145        for p in rest {
146            sum += *p - *first;
147        }
148        #[allow(clippy::cast_precision_loss)]
149        Ok(*first + sum / points.len() as f64)
150    }
151
152    /// Component-wise minimum.
153    #[must_use]
154    pub fn min(self, other: Self) -> Self {
155        Self::new(
156            self.x.min(other.x),
157            self.y.min(other.y),
158            self.z.min(other.z),
159        )
160    }
161
162    /// Component-wise maximum.
163    #[must_use]
164    pub fn max(self, other: Self) -> Self {
165        Self::new(
166            self.x.max(other.x),
167            self.y.max(other.y),
168            self.z.max(other.z),
169        )
170    }
171
172    /// This point with the Z coordinate dropped.
173    #[must_use]
174    pub const fn xy(self) -> Point2 {
175        Point2::new(self.x, self.y)
176    }
177}
178
179impl Point2 {
180    /// The origin.
181    pub const ORIGIN: Self = Self::new(0.0, 0.0);
182
183    /// A point from coordinates.
184    #[must_use]
185    pub const fn new(x: f64, y: f64) -> Self {
186        Self { x, y }
187    }
188
189    /// Coordinates as an array.
190    #[must_use]
191    pub const fn to_array(self) -> [f64; 2] {
192        [self.x, self.y]
193    }
194
195    /// A point from an array.
196    #[must_use]
197    pub const fn from_array([x, y]: [f64; 2]) -> Self {
198        Self::new(x, y)
199    }
200
201    /// The position vector from the origin.
202    #[must_use]
203    pub const fn to_vector(self) -> Vector2 {
204        Vector2::new(self.x, self.y)
205    }
206
207    /// The point at the tip of `v` placed at the origin.
208    #[must_use]
209    pub const fn from_vector(v: Vector2) -> Self {
210        Self::new(v.x, v.y)
211    }
212
213    /// Squared distance to `other`.
214    #[must_use]
215    pub fn square_distance(self, other: Self) -> f64 {
216        (other - self).square_magnitude()
217    }
218
219    /// Distance to `other`.
220    #[must_use]
221    pub fn distance(self, other: Self) -> f64 {
222        self.square_distance(other).sqrt()
223    }
224
225    /// Whether two points coincide within `tol.confusion()`.
226    #[must_use]
227    pub fn is_equal(self, other: Self, tol: Tolerances) -> bool {
228        self.square_distance(other) <= tol.confusion() * tol.confusion()
229    }
230
231    /// Whether both coordinates are finite.
232    #[must_use]
233    pub fn is_finite(self) -> bool {
234        self.x.is_finite() && self.y.is_finite()
235    }
236
237    /// The midpoint of `self` and `other`.
238    #[must_use]
239    pub fn midpoint(self, other: Self) -> Self {
240        Self::new(
241            f64::midpoint(self.x, other.x),
242            f64::midpoint(self.y, other.y),
243        )
244    }
245
246    /// Linear interpolation, `t = 0` giving `self`.
247    #[must_use]
248    pub fn lerp(self, other: Self, t: f64) -> Self {
249        self + (other - self) * t
250    }
251
252    /// This point in the XY plane of space.
253    #[must_use]
254    pub const fn to_3d(self) -> Point {
255        Point::new(self.x, self.y, 0.0)
256    }
257}
258
259impl Add<Vector> for Point {
260    type Output = Self;
261    fn add(self, v: Vector) -> Self {
262        Self::new(self.x + v.x, self.y + v.y, self.z + v.z)
263    }
264}
265
266impl Sub<Vector> for Point {
267    type Output = Self;
268    fn sub(self, v: Vector) -> Self {
269        Self::new(self.x - v.x, self.y - v.y, self.z - v.z)
270    }
271}
272
273impl Sub for Point {
274    type Output = Vector;
275    /// The displacement from `other` to `self`.
276    fn sub(self, other: Self) -> Vector {
277        Vector::new(self.x - other.x, self.y - other.y, self.z - other.z)
278    }
279}
280
281impl AddAssign<Vector> for Point {
282    fn add_assign(&mut self, v: Vector) {
283        *self = *self + v;
284    }
285}
286
287impl SubAssign<Vector> for Point {
288    fn sub_assign(&mut self, v: Vector) {
289        *self = *self - v;
290    }
291}
292
293impl Add<Vector2> for Point2 {
294    type Output = Self;
295    fn add(self, v: Vector2) -> Self {
296        Self::new(self.x + v.x, self.y + v.y)
297    }
298}
299
300impl Sub<Vector2> for Point2 {
301    type Output = Self;
302    fn sub(self, v: Vector2) -> Self {
303        Self::new(self.x - v.x, self.y - v.y)
304    }
305}
306
307impl Sub for Point2 {
308    type Output = Vector2;
309    fn sub(self, other: Self) -> Vector2 {
310        Vector2::new(self.x - other.x, self.y - other.y)
311    }
312}
313
314impl AddAssign<Vector2> for Point2 {
315    fn add_assign(&mut self, v: Vector2) {
316        *self = *self + v;
317    }
318}
319
320impl SubAssign<Vector2> for Point2 {
321    fn sub_assign(&mut self, v: Vector2) {
322        *self = *self - v;
323    }
324}
325
326#[cfg(test)]
327#[allow(clippy::unwrap_used)]
328mod tests {
329    use super::*;
330    use approx::assert_relative_eq;
331
332    const T: Tolerances = Tolerances::millimetres();
333
334    #[test]
335    fn affine_algebra() {
336        let a = Point::new(1.0, 2.0, 3.0);
337        let b = Point::new(4.0, 6.0, 3.0);
338        let d: Vector = b - a;
339        assert_eq!(d, Vector::new(3.0, 4.0, 0.0));
340        assert_eq!(a + d, b);
341        assert_eq!(b - d, a);
342        assert_relative_eq!(a.distance(b), 5.0);
343        assert_relative_eq!(a.square_distance(b), 25.0);
344    }
345
346    #[test]
347    fn midpoint_does_not_overflow_for_extreme_coordinates() {
348        // The naive (a + b) / 2 overflows to infinity here; f64::midpoint does
349        // not. Coordinates this large are pathological, but a kernel that
350        // produces infinities on them is worse than one that does not.
351        let a = Point::new(f64::MAX, 0.0, 0.0);
352        let b = Point::new(f64::MAX, 0.0, 0.0);
353        assert!(a.midpoint(b).is_finite());
354        assert_eq!(a.midpoint(b).x, f64::MAX);
355    }
356
357    #[test]
358    fn midpoint_and_lerp_agree() {
359        let a = Point::new(-1.0, 5.0, 2.0);
360        let b = Point::new(3.0, 1.0, -4.0);
361        assert_eq!(a.midpoint(b), a.lerp(b, 0.5));
362        assert_eq!(a.lerp(b, 0.0), a);
363        assert_eq!(a.lerp(b, 1.0), b);
364    }
365
366    #[test]
367    fn centroid_keeps_precision_far_from_the_origin() {
368        // Summing absolute coordinates around 1e9 and dividing loses the
369        // millimetre-scale detail we care about. Summing offsets does not.
370        let base = 1.0e9;
371        let pts = [
372            Point::new(base, base, base),
373            Point::new(base + 2.0, base, base),
374            Point::new(base + 1.0, base + 3.0, base),
375        ];
376        let c = Point::centroid(&pts).unwrap();
377        assert_relative_eq!(c.x, base + 1.0, epsilon = 1e-9);
378        assert_relative_eq!(c.y, base + 1.0, epsilon = 1e-9);
379        assert_relative_eq!(c.z, base, epsilon = 1e-9);
380    }
381
382    #[test]
383    fn centroid_of_nothing_is_an_error_not_the_origin() {
384        assert!(Point::centroid(&[]).is_err());
385        let single = Point::new(1.0, 2.0, 3.0);
386        assert_eq!(Point::centroid(&[single]).unwrap(), single);
387    }
388
389    #[test]
390    fn equality_uses_tolerance() {
391        let a = Point::new(1.0, 1.0, 1.0);
392        let near = Point::new(1.0 + 1e-9, 1.0, 1.0);
393        let far = Point::new(1.0 + 1e-3, 1.0, 1.0);
394        assert!(a.is_equal(near, T));
395        assert!(!a.is_equal(far, T));
396        assert!(a.is_within(far, 1e-2));
397    }
398
399    #[test]
400    fn coordinate_access_is_bounds_checked() {
401        let p = Point::new(7.0, 8.0, 9.0);
402        assert_eq!(p.coord(1).unwrap(), 8.0);
403        assert!(p.coord(3).is_err());
404    }
405
406    #[test]
407    fn dimension_round_trip() {
408        let p = Point::new(1.0, 2.0, 3.0);
409        assert_eq!(p.xy(), Point2::new(1.0, 2.0));
410        assert_eq!(p.xy().to_3d(), Point::new(1.0, 2.0, 0.0));
411    }
412
413    #[test]
414    fn point2_affine_algebra() {
415        let a = Point2::new(1.0, 2.0);
416        let b = Point2::new(4.0, 6.0);
417        assert_eq!(b - a, Vector2::new(3.0, 4.0));
418        assert_relative_eq!(a.distance(b), 5.0);
419        assert_eq!(a + (b - a), b);
420    }
421}
ogeom_math/point.rs

ogeom_math/
point.rs
