axiolid_exact_arrangement/
distance.rs

1//! Exact distance predicates on [`ExactPoints`].
2//!
3//! Distances are Euclidean in the stored coordinates and never leave exact
4//! arithmetic: a square root never appears, because every question is a
5//! comparison of squares. For homogeneous points `P = (px, py, pw)` with
6//! `pw > 0`, put `D(P, Q) = (qx pw - px qw, qy pw - py qw)`, which is
7//! `pw qw (Q - P)`. For a point `V` and a segment `A`-`B`, let `u = D(A, B)`,
8//! `v = D(A, V)`, `dot = u . v`, `uu = u . u` and `cross = ux vy - uy vx`.
9//! Every weight is positive, so:
10//!
11//! - `V` projects at or before `A` when `dot <= 0`, and at or beyond `B` when
12//!   `dot wb >= uu wv`;
13//! - the foot of the perpendicular is
14//!   `(ax uu wv + dot ux, ay uu wv + dot uy, wa uu wv)`;
15//! - the distance from `V` to the line `AB` is compared with `g` by the sign
16//!   of `cross^2 - g^2 wa^2 wv^2 uu`;
17//! - the distance between two points is compared with `g` by the sign of
18//!   `v . v - g^2 (wa wv)^2` with `v = D(a, b)`.
19//!
20//! Each sign is the interval tier first, then `Dyadic`, through `certify`.
21
22use std::cmp::Ordering;
23
24use axiolid_contracts::{GeomError, GeomResult, Sign};
25use axiolid_exact::{certify, Arith, Dyadic, SignExpr};
26
27use crate::points::{ordering, undecided, ExactPoints, PointId};
28
29/// Where on a segment the point nearest to a given point lies.
30#[derive(Debug, Clone, Copy, PartialEq, Eq)]
31pub enum Nearest {
32    /// At the first end, which includes the projection falling on it.
33    Start,
34    /// Strictly between the ends: the foot of the perpendicular.
35    Interior,
36    /// At the second end, which includes the projection falling on it.
37    End,
38}
39
40type Homogeneous = [Dyadic; 3];
41
42fn lift<T: Arith>(point: &Homogeneous) -> [T; 3] {
43    [
44        T::from_dyadic(&point[0]),
45        T::from_dyadic(&point[1]),
46        T::from_dyadic(&point[2]),
47    ]
48}
49
50/// `D(p, q) = (qx pw - px qw, qy pw - py qw)`.
51fn delta<T: Arith>(p: &[T; 3], q: &[T; 3]) -> [T; 2] {
52    [
53        q[0].mul(&p[2]).sub(&p[0].mul(&q[2])),
54        q[1].mul(&p[2]).sub(&p[1].mul(&q[2])),
55    ]
56}
57
58fn dot<T: Arith>(u: &[T; 2], v: &[T; 2]) -> T {
59    u[0].mul(&v[0]).add(&u[1].mul(&v[1]))
60}
61
62#[derive(Clone, Copy)]
63enum Question {
64    /// The sign of `dot`.
65    Behind,
66    /// The sign of `dot wb - uu wv`.
67    Beyond,
68    /// The sign of `cross^2 - g^2 wa^2 wv^2 uu`.
69    Line,
70}
71
72struct Projection<'a> {
73    question: Question,
74    start: &'a Homogeneous,
75    end: &'a Homogeneous,
76    point: &'a Homogeneous,
77    length: f64,
78}
79
80impl SignExpr for Projection<'_> {
81    fn sign_in<T: Arith>(&self) -> Option<Sign> {
82        let (a, b, v) = (
83            lift::<T>(self.start),
84            lift::<T>(self.end),
85            lift::<T>(self.point),
86        );
87        let (u, w) = (delta(&a, &b), delta(&a, &v));
88        let along = dot(&u, &w);
89        let uu = dot(&u, &u);
90        match self.question {
91            Question::Behind => along.sign(),
92            Question::Beyond => along.mul(&b[2]).sub(&uu.mul(&v[2])).sign(),
93            Question::Line => {
94                let cross = u[0].mul(&w[1]).sub(&u[1].mul(&w[0]));
95                let g = T::from_f64(self.length);
96                let scale = g.mul(&g).mul(&a[2]).mul(&a[2]).mul(&v[2]).mul(&v[2]);
97                cross.mul(&cross).sub(&scale.mul(&uu)).sign()
98            }
99        }
100    }
101}
102
103struct PointPair<'a> {
104    first: &'a Homogeneous,
105    second: &'a Homogeneous,
106    length: f64,
107}
108
109impl SignExpr for PointPair<'_> {
110    fn sign_in<T: Arith>(&self) -> Option<Sign> {
111        let (a, b) = (lift::<T>(self.first), lift::<T>(self.second));
112        let v = delta(&a, &b);
113        let g = T::from_f64(self.length);
114        let weights = a[2].mul(&b[2]);
115        dot(&v, &v)
116            .sub(&g.mul(&g).mul(&weights).mul(&weights))
117            .sign()
118    }
119}
120
121fn certified<E: SignExpr>(expr: &E) -> GeomResult<Sign> {
122    certify(expr).map_err(|_| undecided())
123}
124
125fn check_length(length: f64) -> GeomResult<()> {
126    if length.is_finite() && length >= 0.0 {
127        Ok(())
128    } else {
129        Err(GeomError::InvalidInput(
130            "a distance bound is not finite and non-negative".into(),
131        ))
132    }
133}
134
135impl ExactPoints {
136    fn check_segment(&self, start: PointId, end: PointId) -> GeomResult<()> {
137        if self.same(start, end)? {
138            return Err(GeomError::InvalidInput("a segment has equal ends".into()));
139        }
140        Ok(())
141    }
142
143    /// Where the point of the segment `start`-`end` nearest to `point` lies:
144    /// an end when the projection of `point` falls on or beyond it, else the
145    /// foot of the perpendicular. `InvalidInput` when the ends are equal.
146    pub fn nearest_on_segment(
147        &self,
148        point: PointId,
149        start: PointId,
150        end: PointId,
151    ) -> GeomResult<Nearest> {
152        self.check_segment(start, end)?;
153        let ask = |question| {
154            certified(&Projection {
155                question,
156                start: self.exact(start),
157                end: self.exact(end),
158                point: self.exact(point),
159                length: 0.0,
160            })
161        };
162        if ask(Question::Behind)? != Sign::Positive {
163            return Ok(Nearest::Start);
164        }
165        if ask(Question::Beyond)? != Sign::Negative {
166            return Ok(Nearest::End);
167        }
168        Ok(Nearest::Interior)
169    }
170
171    /// The foot of the perpendicular from `point` to the line through `start`
172    /// and `end`, an exact rational point. It lies on the segment exactly
173    /// when [`ExactPoints::nearest_on_segment`] says `Interior`.
174    /// `InvalidInput` when the ends are equal.
175    pub fn foot(&mut self, point: PointId, start: PointId, end: PointId) -> GeomResult<PointId> {
176        self.check_segment(start, end)?;
177        let [a, b, v] = [start, end, point].map(|id| self.exact(id).clone());
178        let (u, w) = (delta(&a, &b), delta(&a, &v));
179        let (along, uu) = (dot(&u, &w), dot(&u, &u));
180        let scale = uu.mul(&v[2]);
181        let foot = [
182            a[0].mul(&scale).add(&along.mul(&u[0])),
183            a[1].mul(&scale).add(&along.mul(&u[1])),
184            a[2].mul(&scale),
185        ];
186        self.push(foot)
187    }
188
189    /// The order of the distance between `a` and `b` against `length`:
190    /// `Less` when they are closer, `Equal` at exactly `length`.
191    /// `InvalidInput` unless `length` is finite and not negative.
192    pub fn compare_distance(&self, a: PointId, b: PointId, length: f64) -> GeomResult<Ordering> {
193        check_length(length)?;
194        certified(&PointPair {
195            first: self.exact(a),
196            second: self.exact(b),
197            length,
198        })
199        .map(ordering)
200    }
201
202    /// The order of the distance from `point` to the line through `start` and
203    /// `end` against `length`: `Less` when it is closer, `Equal` at exactly
204    /// `length`. `InvalidInput` unless `length` is finite and not negative,
205    /// or when the ends are equal.
206    pub fn compare_line_distance(
207        &self,
208        point: PointId,
209        start: PointId,
210        end: PointId,
211        length: f64,
212    ) -> GeomResult<Ordering> {
213        check_length(length)?;
214        self.check_segment(start, end)?;
215        certified(&Projection {
216            question: Question::Line,
217            start: self.exact(start),
218            end: self.exact(end),
219            point: self.exact(point),
220            length,
221        })
222        .map(ordering)
223    }
224}
axiolid_exact_arrangement/distance.rs

axiolid_exact_arrangement/
distance.rs
