axiolid_exact_arrangement/
points.rs

1//! Exact points of the plane in homogeneous coordinates.
2//!
3//! A point is `(x, y, w)` with `w > 0` and stands for `(x / w, y / w)`. The
4//! coordinates are [`Dyadic`] numbers, so a point can be the crossing of two
5//! segments or of an edge with a plane without being a binary64 position.
6//!
7//! The predicates are polynomials in these values. `orient` is the 3 x 3
8//! determinant and `compare` is `c_a w_b - c_b w_a` for a coordinate `c`. When
9//! every point has `w = 1` and binary64 coordinates (a "flat" point) they are
10//! the exact `orient2d` and a plain comparison. Otherwise each is evaluated
11//! twice at most: in `Interval`s of the stored enclosures, which decides
12//! nearly every call, and, when that straddles zero, in `Dyadic`s of the
13//! stored exact values. `Dyadic` has no range limit, so products of tiny
14//! coordinates are no reason to refuse.
15//!
16//! The boxes in [`ExactPoints::bounds`] are sound enclosures of the points.
17//! They only prune candidates and never decide.
18
19use std::cell::Cell;
20use std::cmp::Ordering;
21
22use axiolid_contracts::{GeomError, GeomResult, Sign};
23use axiolid_core::Point2;
24use axiolid_exact::{Arith, Dyadic, Interval};
25use axiolid_predicates::orient2d;
26
27use crate::ratio::round_ratio;
28
29/// Homogeneous coordinates `(x, y, w)`.
30type Homogeneous<T> = [T; 3];
31
32/// A point of an [`ExactPoints`] store.
33#[derive(Debug, Clone, Copy, PartialEq, Eq, PartialOrd, Ord, Hash)]
34pub struct PointId(pub usize);
35
36#[derive(Debug, Clone)]
37struct ExactPoint {
38    exact: Homogeneous<Dyadic>,
39    enclosure: Homogeneous<Interval>,
40    /// `[x_lo, y_lo, x_hi, y_hi]` of the point.
41    bounds: [f64; 4],
42    /// The coordinates when `w = 1` and they are binary64.
43    flat: Option<Point2>,
44}
45
46/// Where a point lies relative to a ring.
47#[derive(Debug, Clone, Copy, PartialEq, Eq)]
48pub enum Location {
49    /// In the interior.
50    Inside,
51    /// Outside the closed ring.
52    Outside,
53    /// On the ring.
54    Boundary,
55}
56
57/// The extent of a ring, for pruning candidate pairs. Never decides.
58#[derive(Debug, Clone, Copy)]
59pub struct RingBox {
60    /// `[min_x, min_y, max_x, max_y]` of the exact points, with every vertex
61    /// box at its outer limit: contains the ring.
62    bounds: [f64; 4],
63    /// The same with every vertex box at its inner limit: the extremes of the
64    /// ring reach at least this far.
65    core: [f64; 4],
66}
67
68impl RingBox {
69    /// Whether this ring can lie in `other` by their boxes: its smallest
70    /// coordinate cannot be below the box of `other`, nor its largest above.
71    pub fn possibly_within(&self, other: &Self) -> bool {
72        self.core[0] >= other.bounds[0]
73            && self.core[1] >= other.bounds[1]
74            && self.core[2] <= other.bounds[2]
75            && self.core[3] <= other.bounds[3]
76    }
77}
78
79/// A store of exact points and the exact predicates on them.
80///
81/// Orientation is counter-clockwise in `(x, y)` times the frame sign given to
82/// [`ExactPoints::with_frame_sign`]; a caller whose `(x, y)` are the
83/// coordinates of a mirrored frame passes `Sign::Negative` and reads
84/// "counter-clockwise" in its own frame.
85#[derive(Debug, Clone)]
86pub struct ExactPoints {
87    frame_sign: Sign,
88    points: Vec<ExactPoint>,
89    escalations: Cell<usize>,
90}
91
92impl Default for ExactPoints {
93    fn default() -> Self {
94        Self::new()
95    }
96}
97
98fn not_finite() -> GeomError {
99    GeomError::Degenerate("an exact point is not finite".into())
100}
101
102pub(crate) fn undecided() -> GeomError {
103    GeomError::Degenerate("an exact sign could not be decided".into())
104}
105
106pub(crate) fn inconsistent() -> GeomError {
107    GeomError::Degenerate("the exact arrangement state is inconsistent".into())
108}
109
110fn exact(value: f64) -> GeomResult<Dyadic> {
111    Dyadic::try_from_f64(value).ok_or_else(not_finite)
112}
113
114/// `p_i q_j - q_i p_j`.
115fn minor<T: Arith>(p: &Homogeneous<T>, q: &Homogeneous<T>, i: usize, j: usize) -> T {
116    p[i].mul(&q[j]).sub(&q[i].mul(&p[j]))
117}
118
119/// `p x q`: the line through two points, or the point where two lines meet.
120fn cross<T: Arith>(p: &Homogeneous<T>, q: &Homogeneous<T>) -> Homogeneous<T> {
121    [minor(p, q, 1, 2), minor(p, q, 2, 0), minor(p, q, 0, 1)]
122}
123
124/// `det[a; b; c]`: positive when the points turn counter-clockwise (their
125/// weights are positive).
126fn determinant<T: Arith>(a: &Homogeneous<T>, b: &Homogeneous<T>, c: &Homogeneous<T>) -> T {
127    a[0].mul(&minor(b, c, 1, 2))
128        .sub(&a[1].mul(&minor(b, c, 0, 2)))
129        .add(&a[2].mul(&minor(b, c, 0, 1)))
130}
131
132/// `c_a w_b - c_b w_a` for coordinate `axis` (0 = x, 1 = y).
133fn difference<T: Arith>(a: &Homogeneous<T>, b: &Homogeneous<T>, axis: usize) -> T {
134    a[axis].mul(&b[2]).sub(&b[axis].mul(&a[2]))
135}
136
137pub(crate) fn ordering(sign: Sign) -> Ordering {
138    match sign {
139        Sign::Positive => Ordering::Greater,
140        Sign::Negative => Ordering::Less,
141        _ => Ordering::Equal,
142    }
143}
144
145fn coordinate(point: Point2, axis: usize) -> f64 {
146    if axis == 0 {
147        point.x
148    } else {
149        point.y
150    }
151}
152
153/// Product of two signs.
154fn sign_times(left: Sign, right: Sign) -> Sign {
155    if left == Sign::Zero || right == Sign::Zero {
156        Sign::Zero
157    } else if left == right {
158        Sign::Positive
159    } else {
160        Sign::Negative
161    }
162}
163
164/// `numerator / weight` as a binary64 number, when that is exact. `hint` is a
165/// candidate, such as an end of an edge that does not move in this coordinate.
166fn exact_quotient(numerator: &Dyadic, weight: &Dyadic, hint: Option<f64>) -> Option<f64> {
167    let candidate = match hint {
168        Some(hint) => hint,
169        // Both convert exactly up to 53 bits, so the division is then
170        // correctly rounded and equals the quotient whenever that is
171        // representable.
172        None if numerator.bits() <= 53 && weight.bits() <= 53 => {
173            numerator.to_f64() / weight.to_f64()
174        }
175        None => return None,
176    };
177    let exact_candidate = Dyadic::try_from_f64(candidate)?;
178    (exact_candidate.mul(weight) == *numerator).then_some(candidate)
179}
180
181impl ExactPoints {
182    /// An empty store whose orientation is counter-clockwise in `(x, y)`.
183    pub fn new() -> Self {
184        Self::with_frame_sign(Sign::Positive)
185    }
186
187    /// An empty store whose orientation is `frame_sign` times the sense of
188    /// `(x, y)`. `frame_sign` must not be `Sign::Zero`.
189    pub fn with_frame_sign(frame_sign: Sign) -> Self {
190        Self {
191            frame_sign,
192            points: Vec::new(),
193            escalations: Cell::new(0),
194        }
195    }
196
197    /// The frame sign the store was made with.
198    pub fn frame_sign(&self) -> Sign {
199        self.frame_sign
200    }
201
202    /// The number of points.
203    pub fn len(&self) -> usize {
204        self.points.len()
205    }
206
207    /// Whether the store has no point.
208    pub fn is_empty(&self) -> bool {
209        self.points.is_empty()
210    }
211
212    /// A binary64 position.
213    pub fn push_flat(&mut self, flat: Point2) -> GeomResult<PointId> {
214        let point = [exact(flat.x)?, exact(flat.y)?, exact(1.0)?];
215        let enclosure = [
216            Interval::point(flat.x),
217            Interval::point(flat.y),
218            Interval::point(1.0),
219        ];
220        self.points.push(ExactPoint {
221            exact: point,
222            enclosure,
223            bounds: [flat.x, flat.y, flat.x, flat.y],
224            flat: Some(flat),
225        });
226        Ok(PointId(self.points.len() - 1))
227    }
228
229    /// The point `(x / w, y / w)` of homogeneous coordinates with `w != 0`.
230    /// A negative `w` is normalised. A point that is a binary64 position is
231    /// stored as one, so it takes the fast path of the predicates.
232    pub fn push(&mut self, homogeneous: [Dyadic; 3]) -> GeomResult<PointId> {
233        self.push_hinted(homogeneous, [None, None])
234    }
235
236    /// Like [`ExactPoints::push`], with a candidate value for each
237    /// coordinate to try first, such as the coordinate of an edge that does
238    /// not move in it.
239    pub fn push_hinted(
240        &mut self,
241        homogeneous: [Dyadic; 3],
242        hints: [Option<f64>; 2],
243    ) -> GeomResult<PointId> {
244        let mut point = homogeneous;
245        match point[2].sign() {
246            Some(Sign::Positive) => {}
247            Some(Sign::Negative) => point = point.map(|value| value.neg()),
248            _ => return Err(inconsistent()),
249        }
250        if let (Some(x), Some(y)) = (
251            exact_quotient(&point[0], &point[2], hints[0]),
252            exact_quotient(&point[1], &point[2], hints[1]),
253        ) {
254            return self.push_flat(Point2::new(x, y));
255        }
256        let enclosure = [
257            point[0].enclosure(),
258            point[1].enclosure(),
259            point[2].enclosure(),
260        ];
261        let (x, y) = (
262            enclosure[0].quotient(enclosure[2]),
263            enclosure[1].quotient(enclosure[2]),
264        );
265        self.points.push(ExactPoint {
266            exact: point,
267            enclosure,
268            bounds: [x.lo(), y.lo(), x.hi(), y.hi()],
269            flat: None,
270        });
271        Ok(PointId(self.points.len() - 1))
272    }
273
274    /// The coordinates of a point that is a binary64 position (`w = 1`).
275    pub fn flat(&self, point: PointId) -> Option<Point2> {
276        self.points[point.0].flat
277    }
278
279    /// The exact homogeneous coordinates of a point, `w > 0`.
280    pub fn exact(&self, point: PointId) -> &[Dyadic; 3] {
281        &self.points[point.0].exact
282    }
283
284    /// The point where the segments `ends[0]`-`ends[1]` and
285    /// `ends[2]`-`ends[3]` cross properly.
286    ///
287    /// The lines `L = a x b` and the meeting point `X = L1 x L2` have degrees
288    /// 2 and (3, 3, 2) in the coordinates; `w` is made positive. When the
289    /// point is a binary64 position it takes the fast path. Parallel segments
290    /// have no crossing.
291    pub fn segment_crossing(&mut self, ends: [PointId; 4]) -> GeomResult<PointId> {
292        let [a, b, c, d] = ends.map(|end| &self.points[end.0]);
293        let point = cross(&cross(&a.exact, &b.exact), &cross(&c.exact, &d.exact));
294        // An edge that does not move in a coordinate crosses at that value.
295        let fixed = |axis: usize| -> Option<f64> {
296            [(a, b), (c, d)].into_iter().find_map(|(from, to)| {
297                let (from, to) = (coordinate(from.flat?, axis), coordinate(to.flat?, axis));
298                (from == to).then_some(from)
299            })
300        };
301        let hints = [fixed(0), fixed(1)];
302        self.push_hinted(point, hints)
303    }
304
305    /// The parameter `t` along the first of two crossing segments (`a` to `b`
306    /// and `c` to `d`, as in [`ExactPoints::segment_crossing`]) at which they
307    /// cross: `a + t (b - a)`. For output only.
308    pub fn segment_parameter(&self, ends: [PointId; 4]) -> GeomResult<f64> {
309        let [a, b, c, d] = ends.map(|end| &self.points[end.0].exact);
310        let difference = |to: &Homogeneous<Dyadic>, from: &Homogeneous<Dyadic>| {
311            [to[0].sub(&from[0]), to[1].sub(&from[1])]
312        };
313        let turn = |p: [Dyadic; 2], q: [Dyadic; 2]| p[0].mul(&q[1]).sub(&p[1].mul(&q[0]));
314        let along = difference(d, c);
315        let numerator = turn(difference(c, a), along.clone());
316        let denominator = turn(difference(b, a), along);
317        if denominator.sign() == Some(Sign::Zero) {
318            return Err(inconsistent());
319        }
320        let (m_n, e_n) = numerator.approx_parts();
321        let (m_d, e_d) = denominator.approx_parts();
322        let shift = (e_n - e_d).clamp(-1100, 1100) as i32;
323        let t = m_n / m_d * 2f64.powi(shift / 2) * 2f64.powi(shift - shift / 2);
324        if t.is_finite() {
325            Ok(t)
326        } else {
327            Err(not_finite())
328        }
329    }
330
331    /// Exact orientation of `a`, `b`, `c`: positive when they turn
332    /// counter-clockwise in the store's sense (see [`ExactPoints`]).
333    pub fn orient(&self, a: PointId, b: PointId, c: PointId) -> GeomResult<Sign> {
334        let (a, b, c) = (&self.points[a.0], &self.points[b.0], &self.points[c.0]);
335        if let (Some(fa), Some(fb), Some(fc)) = (a.flat, b.flat, c.flat) {
336            let sign = orient2d(fa, fb, fc).sign().ok_or_else(undecided)?;
337            return Ok(sign_times(sign, self.frame_sign));
338        }
339        let sign = match determinant(&a.enclosure, &b.enclosure, &c.enclosure).sign() {
340            Some(sign) => sign,
341            None => {
342                self.escalations.set(self.escalations.get() + 1);
343                determinant(&a.exact, &b.exact, &c.exact)
344                    .sign()
345                    .ok_or_else(undecided)?
346            }
347        };
348        Ok(sign_times(sign, self.frame_sign))
349    }
350
351    /// Order of the coordinates of `a` and `b` along `axis` (0 = x, 1 = y).
352    /// The stored axes, not the frame's.
353    pub fn compare(&self, a: PointId, b: PointId, axis: usize) -> GeomResult<Ordering> {
354        let (a, b) = (&self.points[a.0], &self.points[b.0]);
355        if let (Some(fa), Some(fb)) = (a.flat, b.flat) {
356            return coordinate(fa, axis)
357                .partial_cmp(&coordinate(fb, axis))
358                .ok_or_else(not_finite);
359        }
360        let sign = match difference(&a.enclosure, &b.enclosure, axis).sign() {
361            Some(sign) => sign,
362            None => {
363                self.escalations.set(self.escalations.get() + 1);
364                difference(&a.exact, &b.exact, axis)
365                    .sign()
366                    .ok_or_else(undecided)?
367            }
368        };
369        Ok(ordering(sign))
370    }
371
372    /// Whether `a` and `b` are the same point.
373    pub fn same(&self, a: PointId, b: PointId) -> GeomResult<bool> {
374        Ok(a == b
375            || (self.compare(a, b, 0)? == Ordering::Equal
376                && self.compare(a, b, 1)? == Ordering::Equal))
377    }
378
379    /// Order by `x`, then `y`.
380    pub fn lexicographic(&self, a: PointId, b: PointId) -> GeomResult<Ordering> {
381        match self.compare(a, b, 0)? {
382            Ordering::Equal => self.compare(a, b, 1),
383            order => Ok(order),
384        }
385    }
386
387    /// `[x_lo, y_lo, x_hi, y_hi]`, a sound box of the point.
388    pub fn bounds(&self, point: PointId) -> [f64; 4] {
389        self.points[point.0].bounds
390    }
391
392    /// The axis along which two distinct points of a line differ.
393    pub fn dominant(&self, a: PointId, b: PointId) -> GeomResult<usize> {
394        Ok(if self.compare(a, b, 0)? == Ordering::Equal {
395            1
396        } else {
397            0
398        })
399    }
400
401    /// `(low, high)` of two points along `axis`.
402    pub fn ordered(&self, a: PointId, b: PointId, axis: usize) -> GeomResult<(PointId, PointId)> {
403        Ok(if self.compare(a, b, axis)? == Ordering::Greater {
404            (b, a)
405        } else {
406            (a, b)
407        })
408    }
409
410    /// Whether collinear `at` lies strictly between `before` and `after`.
411    pub fn strictly_between(
412        &self,
413        before: PointId,
414        at: PointId,
415        after: PointId,
416    ) -> GeomResult<bool> {
417        let axis = self.dominant(before, after)?;
418        let (low, high) = self.ordered(before, after, axis)?;
419        Ok(self.compare(low, at, axis)? == Ordering::Less
420            && self.compare(at, high, axis)? == Ordering::Less)
421    }
422
423    /// Whether collinear `point` lies on the closed segment `start`-`end`.
424    pub fn within_segment(&self, point: PointId, start: PointId, end: PointId) -> GeomResult<bool> {
425        let axis = self.dominant(start, end)?;
426        let (low, high) = self.ordered(start, end, axis)?;
427        Ok(self.compare(low, point, axis)? != Ordering::Greater
428            && self.compare(point, high, axis)? != Ordering::Greater)
429    }
430
431    /// Where `point` lies relative to the closed `ring`: a crossing-number
432    /// walk with exact orientation tests.
433    ///
434    /// `orient` is in the frame's sense and the `y` comparisons are in the
435    /// stored axes. When the frame mirrors them the walk counts the crossings
436    /// of the opposite ray, which is the winding number negated; only zero or
437    /// nonzero is used.
438    pub fn locate(&self, point: PointId, ring: &[PointId]) -> GeomResult<Location> {
439        let mut winding = 0i32;
440        for index in 0..ring.len() {
441            let (a, b) = (ring[index], ring[(index + 1) % ring.len()]);
442            let side = self.orient(a, b, point)?;
443            if side == Sign::Zero && self.within_segment(point, a, b)? {
444                return Ok(Location::Boundary);
445            }
446            if self.compare(a, point, 1)? != Ordering::Greater {
447                if self.compare(b, point, 1)? == Ordering::Greater && side == Sign::Positive {
448                    winding += 1;
449                }
450            } else if self.compare(b, point, 1)? != Ordering::Greater && side == Sign::Negative {
451                winding -= 1;
452            }
453        }
454        Ok(if winding == 0 {
455            Location::Outside
456        } else {
457            Location::Inside
458        })
459    }
460
461    /// The boxes of a ring of exact points.
462    pub fn ring_box(&self, ring: &[PointId]) -> RingBox {
463        let mut bounds = [
464            f64::INFINITY,
465            f64::INFINITY,
466            f64::NEG_INFINITY,
467            f64::NEG_INFINITY,
468        ];
469        let mut core = bounds;
470        for &node in ring {
471            let [x_lo, y_lo, x_hi, y_hi] = self.bounds(node);
472            bounds = [
473                bounds[0].min(x_lo),
474                bounds[1].min(y_lo),
475                bounds[2].max(x_hi),
476                bounds[3].max(y_hi),
477            ];
478            core = [
479                core[0].min(x_hi),
480                core[1].min(y_hi),
481                core[2].max(x_lo),
482                core[3].max(y_lo),
483            ];
484        }
485        RingBox { bounds, core }
486    }
487
488    /// Each coordinate `x / w` and `y / w` rounded once to the nearest
489    /// binary64 number, ties to even. `Degenerate` when one is beyond the
490    /// binary64 range.
491    pub fn round(&self, point: PointId) -> GeomResult<Point2> {
492        let point = &self.points[point.0];
493        if let Some(flat) = point.flat {
494            return Ok(flat);
495        }
496        let [x, y, w] = &point.exact;
497        match (round_ratio(x, w), round_ratio(y, w)) {
498            (Some(x), Some(y)) => Ok(Point2::new(x, y)),
499            _ => Err(not_finite()),
500        }
501    }
502
503    /// The number of predicate calls that left the interval filter for exact
504    /// arithmetic. Not counted by the flat fast path.
505    pub fn escalations(&self) -> usize {
506        self.escalations.get()
507    }
508}
509
510#[cfg(test)]
511mod tests {
512    use super::*;
513
514    fn flats(points: &mut ExactPoints, corners: &[(f64, f64)]) -> Vec<PointId> {
515        corners
516            .iter()
517            .map(|&(x, y)| points.push_flat(Point2::new(x, y)).expect("point"))
518            .collect()
519    }
520
521    #[test]
522    fn the_diagonals_of_the_unit_square_meet_at_its_centre_on_the_fast_path() {
523        let mut points = ExactPoints::new();
524        let at = flats(
525            &mut points,
526            &[(0.0, 0.0), (1.0, 1.0), (1.0, 0.0), (0.0, 1.0), (0.5, 0.5)],
527        );
528        let crossing = points
529            .segment_crossing([at[0], at[1], at[2], at[3]])
530            .expect("crossing");
531        assert_eq!(points.flat(crossing), Some(Point2::new(0.5, 0.5)));
532        assert!(points.same(crossing, at[4]).expect("same"));
533    }
534
535    #[test]
536    fn a_segment_crossing_does_not_depend_on_the_direction_of_its_segments() {
537        let mut points = ExactPoints::new();
538        let at = flats(
539            &mut points,
540            &[(0.0, 0.0), (0.5, 1.0), (1.0, 0.0), (0.0, 1.0)],
541        );
542        let forward = points
543            .segment_crossing([at[0], at[1], at[2], at[3]])
544            .expect("crossing");
545        let swapped = points
546            .segment_crossing([at[2], at[3], at[0], at[1]])
547            .expect("crossing");
548        let reversed = points
549            .segment_crossing([at[1], at[0], at[3], at[2]])
550            .expect("crossing");
551        assert_eq!(points.flat(forward), None);
552        let exact = points.exact(forward);
553        // (1/3, 2/3): the numerators differ by a factor of two.
554        assert_eq!(exact[0].add(&exact[0]), exact[1]);
555        assert_eq!(exact[2].sign(), Some(Sign::Positive));
556        for other in [swapped, reversed] {
557            assert!(points.same(forward, other).expect("same"));
558            assert_eq!(points.exact(forward), points.exact(other));
559        }
560    }
561
562    #[test]
563    fn crossings_on_one_segment_are_collinear_in_either_frame_sense() {
564        for (frame_sign, turn) in [
565            (Sign::Positive, Sign::Positive),
566            (Sign::Negative, Sign::Negative),
567        ] {
568            let mut points = ExactPoints::with_frame_sign(frame_sign);
569            // The segment 0-1 from (0.1, 0.1) to (3.7, 2.9), crossed by three
570            // segments whose crossing points are not binary64.
571            let mut corners = vec![(0.1, 0.1), (3.7, 2.9), (2.0, 0.0)];
572            for x in [0.7, 1.3, 2.9] {
573                corners.extend([(x, -1.0), (x + 0.3, 4.0)]);
574            }
575            let at = flats(&mut points, &corners);
576            let crossing: Vec<PointId> = (0..3)
577                .map(|k| {
578                    points
579                        .segment_crossing([at[0], at[1], at[3 + 2 * k], at[4 + 2 * k]])
580                        .expect("crossing")
581                })
582                .collect();
583            assert!(crossing.iter().all(|&p| points.flat(p).is_none()));
584            assert_eq!(
585                points
586                    .orient(crossing[0], crossing[1], crossing[2])
587                    .expect("orient"),
588                Sign::Zero
589            );
590            // (2, 0) lies below the segment, to its right going from the
591            // first crossing to the second.
592            let side = points
593                .orient(crossing[0], crossing[1], at[2])
594                .expect("orient");
595            assert_eq!(side, sign_times(Sign::Negative, turn));
596        }
597    }
598
599    #[test]
600    fn a_mirrored_frame_reverses_the_orientation() {
601        let corners = [(0.0, 0.0), (1.0, 0.0), (0.0, 1.0)];
602        let mut counter = ExactPoints::new();
603        let mut mirrored = ExactPoints::with_frame_sign(Sign::Negative);
604        let (a, b) = (
605            flats(&mut counter, &corners),
606            flats(&mut mirrored, &corners),
607        );
608        assert_eq!(counter.orient(a[0], a[1], a[2]).unwrap(), Sign::Positive);
609        assert_eq!(mirrored.orient(b[0], b[1], b[2]).unwrap(), Sign::Negative);
610    }
611
612    #[test]
613    fn the_parameter_of_a_segment_crossing_is_its_position_along_the_first_segment() {
614        let mut points = ExactPoints::new();
615        let at = flats(
616            &mut points,
617            &[(0.0, 0.0), (0.5, 1.0), (1.0, 0.0), (0.0, 1.0)],
618        );
619        let t = points
620            .segment_parameter([at[0], at[1], at[2], at[3]])
621            .expect("parameter");
622        assert!((t - 2.0 / 3.0).abs() < 1e-15, "{t}");
623        let back = points
624            .segment_parameter([at[2], at[3], at[0], at[1]])
625            .expect("parameter");
626        assert!((back - 2.0 / 3.0).abs() < 1e-15, "{back}");
627        assert!(points
628            .segment_parameter([at[0], at[1], at[0], at[1]])
629            .is_err());
630    }
631
632    #[test]
633    fn parallel_segments_have_no_crossing() {
634        let mut points = ExactPoints::new();
635        let at = flats(
636            &mut points,
637            &[(0.0, 0.0), (1.0, 0.0), (0.0, 1.0), (1.0, 1.0)],
638        );
639        assert!(points
640            .segment_crossing([at[0], at[1], at[2], at[3]])
641            .is_err());
642    }
643
644    #[test]
645    fn a_homogeneous_point_with_a_binary64_position_takes_the_fast_path() {
646        let mut points = ExactPoints::new();
647        let three = Dyadic::try_from_f64(3.0).unwrap();
648        let (x, y) = (
649            Dyadic::try_from_f64(1.5).unwrap(),
650            Dyadic::try_from_f64(0.75).unwrap(),
651        );
652        // (3 * 0.5, 3 * 0.25, 3) is (0.5, 0.25); a negative weight is the same.
653        let half = points
654            .push([
655                x.mul(&Dyadic::try_from_f64(1.0).unwrap()),
656                y.mul(&Dyadic::try_from_f64(1.0).unwrap()),
657                three.clone(),
658            ])
659            .expect("point");
660        let negated = points.push([x.neg(), y.neg(), three.neg()]).expect("point");
661        assert_eq!(points.flat(half), Some(Point2::new(0.5, 0.25)));
662        assert!(points.same(half, negated).expect("same"));
663        assert!(points.push([x, y, Dyadic::zero()]).is_err());
664    }
665
666    #[test]
667    fn a_crossing_rounds_once_to_the_nearest_binary64() {
668        let mut points = ExactPoints::new();
669        let at = flats(
670            &mut points,
671            &[(0.0, 0.0), (0.5, 1.0), (1.0, 0.0), (0.0, 1.0)],
672        );
673        let crossing = points
674            .segment_crossing([at[0], at[1], at[2], at[3]])
675            .expect("crossing");
676        let rounded = points.round(crossing).expect("rounded");
677        assert_eq!(rounded.x.to_bits(), (1.0f64 / 3.0).to_bits());
678        assert_eq!(rounded.y.to_bits(), (2.0f64 / 3.0).to_bits());
679        assert_eq!(points.round(at[1]).unwrap(), Point2::new(0.5, 1.0));
680    }
681
682    #[test]
683    fn locate_tells_inside_outside_and_boundary() {
684        let mut points = ExactPoints::new();
685        let at = flats(
686            &mut points,
687            &[
688                (0.0, 0.0),
689                (2.0, 0.0),
690                (2.0, 2.0),
691                (0.0, 2.0),
692                (1.0, 1.0),
693                (3.0, 1.0),
694                (1.0, 0.0),
695            ],
696        );
697        let ring = &at[..4];
698        assert_eq!(points.locate(at[4], ring).unwrap(), Location::Inside);
699        assert_eq!(points.locate(at[5], ring).unwrap(), Location::Outside);
700        assert_eq!(points.locate(at[6], ring).unwrap(), Location::Boundary);
701        assert!(points
702            .ring_box(ring)
703            .possibly_within(&points.ring_box(ring)));
704    }
705
706    #[test]
707    fn escalations_count_the_calls_the_interval_filter_could_not_decide() {
708        let mut points = ExactPoints::new();
709        let at = flats(
710            &mut points,
711            &[(0.0, 0.0), (0.5, 1.0), (1.0, 0.0), (0.0, 1.0)],
712        );
713        let crossing = points
714            .segment_crossing([at[0], at[1], at[2], at[3]])
715            .expect("crossing");
716        assert_eq!(points.escalations(), 0);
717        // The crossing lies exactly on its own segment: the enclosure
718        // straddles zero and the exact value decides.
719        assert_eq!(
720            points.orient(at[0], crossing, at[1]).expect("orient"),
721            Sign::Zero
722        );
723        assert!(points.escalations() >= 1);
724    }
725}
axiolid_exact_arrangement/points.rs

axiolid_exact_arrangement/
points.rs
