axiolid_exact_arrangement/
arrangement.rs

1//! The arrangement of straight segments on exact points.
2//!
3//! The segments may touch, cross, abut and overlap. They are cut at every
4//! point where they meet into pieces that no other piece crosses, and every
5//! decision is an exact sign on [`ExactPoints`]:
6//!
7//! 1. Equal points: the end points of the segments are joined where they are
8//!    the same point, whatever their [`PointId`].
9//! 2. Cuts: a sweep over the segments' boxes proposes pairs. A proper
10//!    crossing creates a point ([`ExactPoints::segment_crossing`]); an end of
11//!    one segment lying strictly inside the other, which covers T-junctions
12//!    and collinear overlaps, splits it.
13//! 3. Equal points again: along each segment its points are ordered exactly
14//!    and equal neighbours are joined, so three segments through one point
15//!    give one vertex.
16//! 4. Pieces: the stretch between two consecutive distinct vertices of a
17//!    segment, once however many segments run over it.
18//!
19//! Vertices are numbered in a deterministic order: those with input points
20//! first, by their smallest [`PointId`], then the crossings, by their
21//! smallest pair of segment indices. Pieces are numbered by their ends.
22
23use std::collections::{BTreeMap, BTreeSet};
24
25use axiolid_contracts::{ExecutionOptions, GeomError, GeomResult, Sign};
26
27use crate::points::{inconsistent, ExactPoints, PointId};
28
29/// Bytes charged for one crossing point: the Dyadic coordinates, the
30/// enclosures and the bookkeeping of a segment pair.
31pub const CROSSING_SCRATCH_BYTES: usize = 512;
32
33/// Bytes charged for one end of a segment lying inside another one (a
34/// T-junction split): an entry in two cut lists and its share of the sort.
35pub const SPLIT_SCRATCH_BYTES: usize = 64;
36
37/// A vertex of an [`Arrangement`].
38#[derive(Debug, Clone, Copy, PartialEq, Eq, PartialOrd, Ord, Hash)]
39pub struct VertexId(pub usize);
40
41/// A piece of an [`Arrangement`].
42#[derive(Debug, Clone, Copy, PartialEq, Eq, PartialOrd, Ord, Hash)]
43pub struct PieceId(pub usize);
44
45/// A count of candidate pairs against a limit, shared by the steps of one
46/// request.
47#[derive(Debug, Clone)]
48pub struct PairBudget {
49    /// Pairs counted so far.
50    pub used: usize,
51    /// Pairs allowed.
52    pub max: usize,
53    /// The resource named by `GeomError::BudgetExceeded` when exceeded.
54    pub resource: &'static str,
55}
56
57impl PairBudget {
58    /// A budget of `max` pairs, none used.
59    pub const fn new(max: usize, resource: &'static str) -> Self {
60        Self {
61            used: 0,
62            max,
63            resource,
64        }
65    }
66
67    /// Count one candidate pair.
68    pub fn count(&mut self) -> GeomResult<()> {
69        self.used += 1;
70        if self.used > self.max {
71            return Err(GeomError::BudgetExceeded {
72                resource: self.resource,
73            });
74        }
75        Ok(())
76    }
77}
78
79/// A segment between two points of the store. The caller keeps its own data
80/// per segment index.
81#[derive(Debug, Clone, Copy, PartialEq, Eq)]
82pub struct Segment {
83    /// The ends.
84    pub ends: [PointId; 2],
85}
86
87/// How the point of a vertex without input points came about.
88#[derive(Debug, Clone, Copy, PartialEq, Eq)]
89pub struct Crossing {
90    /// The smallest pair of segments crossing here, ascending.
91    pub segments: [usize; 2],
92    /// The ends of those segments, in that order.
93    pub ends: [PointId; 4],
94}
95
96/// A point of the arrangement.
97#[derive(Debug, Clone)]
98pub struct Vertex {
99    /// An exact point at this position.
100    pub point: PointId,
101    /// The input points joined here, ascending; empty for a crossing.
102    pub inputs: Vec<PointId>,
103    /// For a vertex without input points, the crossing that created it.
104    pub crossing: Option<Crossing>,
105}
106
107/// A stretch of segments no other segment crosses or touches inside.
108#[derive(Debug, Clone)]
109pub struct Piece {
110    /// The vertices, `ends[0] < ends[1]`.
111    pub ends: [VertexId; 2],
112    /// The segments that run over the piece, ascending.
113    pub segments: Vec<usize>,
114}
115
116/// The pieces of a set of segments.
117#[derive(Debug, Clone)]
118pub struct Arrangement {
119    /// The vertices.
120    pub vertices: Vec<Vertex>,
121    /// The pieces, ordered by their ends.
122    pub pieces: Vec<Piece>,
123    /// Per input segment, its pieces in order from its first end to its
124    /// second, each with whether that direction runs from `ends[0]` to
125    /// `ends[1]` of the piece.
126    pub paths: Vec<Vec<(PieceId, bool)>>,
127    /// The number of vertices without input points.
128    pub crossings: usize,
129    /// For each vertex, the pieces at it and whether each leaves it at
130    /// `ends[0]`.
131    incident: Vec<Vec<(PieceId, bool)>>,
132}
133
134/// A point while the arrangement is built: its exact point and its label.
135#[derive(Clone, Copy)]
136enum Label {
137    Input,
138    Crossing([usize; 2], [PointId; 4]),
139}
140
141struct Cut {
142    ends: [usize; 2],
143    /// `[min_x, min_y, max_x, max_y]`.
144    bounds: [f64; 4],
145    /// Local points inside the segment, from the pairs examined.
146    splits: Vec<usize>,
147    /// All local points from the first end to the second, in order.
148    path: Vec<usize>,
149}
150
151/// Points while building: exact point, label and union-find parent.
152struct Local {
153    exact: Vec<PointId>,
154    labels: Vec<Label>,
155    parent: Vec<usize>,
156}
157
158impl Local {
159    fn push(&mut self, exact: PointId, label: Label) -> usize {
160        self.exact.push(exact);
161        self.labels.push(label);
162        self.parent.push(self.parent.len());
163        self.parent.len() - 1
164    }
165
166    fn find(&mut self, mut point: usize) -> usize {
167        while self.parent[point] != point {
168            self.parent[point] = self.parent[self.parent[point]];
169            point = self.parent[point];
170        }
171        point
172    }
173
174    /// Join two points; the root is the smaller number.
175    fn union(&mut self, left: usize, right: usize) {
176        let (left, right) = (self.find(left), self.find(right));
177        self.parent[left.max(right)] = left.min(right);
178    }
179}
180
181/// Whether two `[min_x, min_y, max_x, max_y]` boxes overlap, edges included.
182pub fn boxes_overlap(first: &[f64; 4], second: &[f64; 4]) -> bool {
183    first[0] <= second[2] && second[0] <= first[2] && first[1] <= second[3] && second[1] <= first[3]
184}
185
186/// The smallest box holding all `boxes`.
187pub fn union_box(boxes: &[[f64; 4]]) -> [f64; 4] {
188    boxes.iter().fold(
189        [
190            f64::INFINITY,
191            f64::INFINITY,
192            f64::NEG_INFINITY,
193            f64::NEG_INFINITY,
194        ],
195        |all, one| {
196            [
197                all[0].min(one[0]),
198                all[1].min(one[1]),
199                all[2].max(one[2]),
200                all[3].max(one[3]),
201            ]
202        },
203    )
204}
205
206/// Sort `items` by an exact comparison that can fail.
207pub(crate) fn sort_exact<T: Copy>(
208    items: &mut [T],
209    mut compare: impl FnMut(T, T) -> GeomResult<std::cmp::Ordering>,
210) -> GeomResult<()> {
211    let mut failure = None;
212    items.sort_by(|&left, &right| match compare(left, right) {
213        Ok(order) => order,
214        Err(error) => {
215            failure.get_or_insert(error);
216            std::cmp::Ordering::Equal
217        }
218    });
219    failure.map_or(Ok(()), Err)
220}
221
222/// Cut `segments` into pieces.
223///
224/// Every examined pair of segments counts against `budget`. Every crossing
225/// and every T-junction split is charged: refused with `BudgetExceeded` for
226/// `memory` when `declared_bytes`, plus `CROSSING_SCRATCH_BYTES` for each
227/// crossing (the next one included), plus `SPLIT_SCRATCH_BYTES` for each
228/// split, exceeds `options.memory_budget_bytes()`. `declared_bytes` is what the caller has declared per element of its own
229/// input. Cancellation is checked per segment.
230///
231/// Input contract: the ends are points of `points`, a segment may repeat no
232/// end and may not be empty, and segments may cross, touch, share an end or
233/// overlap in any way. Two segments that share a `PointId` and lie on one
234/// line are cut where one's other end lies strictly inside the other. Two
235/// identical segments yield one piece with both as carriers.
236pub fn arrange(
237    points: &mut ExactPoints,
238    segments: &[Segment],
239    budget: &mut PairBudget,
240    declared_bytes: usize,
241    options: &ExecutionOptions,
242) -> GeomResult<Arrangement> {
243    // Local points: the distinct end points ascending, then the crossings.
244    let mut inputs: Vec<PointId> = segments.iter().flat_map(|segment| segment.ends).collect();
245    inputs.sort_unstable();
246    inputs.dedup();
247    let mut local = Local {
248        exact: Vec::new(),
249        labels: Vec::new(),
250        parent: Vec::new(),
251    };
252    for &point in &inputs {
253        local.push(point, Label::Input);
254    }
255    // Points at one position share a vertex.
256    let mut by_position: Vec<usize> = (0..inputs.len()).collect();
257    sort_exact(&mut by_position, |left, right| {
258        points.lexicographic(local.exact[left], local.exact[right])
259    })?;
260    for pair in by_position.windows(2) {
261        if points.same(local.exact[pair[0]], local.exact[pair[1]])? {
262            local.union(pair[0], pair[1]);
263        }
264    }
265
266    let local_of = |point: PointId| inputs.binary_search(&point).map_err(|_| inconsistent());
267    let mut cuts = Vec::with_capacity(segments.len());
268    for segment in segments {
269        let ends = [local_of(segment.ends[0])?, local_of(segment.ends[1])?];
270        let (low, high) = (
271            points.bounds(segment.ends[0]),
272            points.bounds(segment.ends[1]),
273        );
274        cuts.push(Cut {
275            ends,
276            bounds: [
277                low[0].min(high[0]),
278                low[1].min(high[1]),
279                low[2].max(high[2]),
280                low[3].max(high[3]),
281            ],
282            splits: Vec::new(),
283            path: Vec::new(),
284        });
285    }
286
287    cut_segments(
288        segments,
289        &mut cuts,
290        points,
291        &mut local,
292        (budget, declared_bytes),
293        options,
294    )?;
295    order_cuts(&mut cuts, points, &mut local, options)?;
296
297    let (vertices, vertex_of) = number_vertices(&mut local);
298    let crossings = vertices
299        .iter()
300        .filter(|vertex| vertex.inputs.is_empty())
301        .count();
302
303    // Pieces, once each however many segments run over them.
304    let mut by_ends = BTreeMap::<(usize, usize), usize>::new();
305    let mut carriers = Vec::<BTreeSet<usize>>::new();
306    let mut on_segments = Vec::<(usize, usize, bool)>::new();
307    for (index, cut) in cuts.iter().enumerate() {
308        options.check_cancelled()?;
309        let mut previous: Option<usize> = None;
310        for &point in &cut.path {
311            let vertex = vertex_of[&local.find(point)];
312            if let Some(from) = previous {
313                if from != vertex {
314                    let ends = (from.min(vertex), from.max(vertex));
315                    let next = carriers.len();
316                    let piece = *by_ends.entry(ends).or_insert(next);
317                    if piece == next {
318                        carriers.push(BTreeSet::new());
319                    }
320                    carriers[piece].insert(index);
321                    on_segments.push((index, piece, from < vertex));
322                }
323            }
324            previous = Some(vertex);
325        }
326    }
327    // Piece numbers in order of their ends.
328    let mut ends_of = vec![(0, 0); carriers.len()];
329    for (&ends, &piece) in &by_ends {
330        ends_of[piece] = ends;
331    }
332    let mut renumber = vec![0usize; carriers.len()];
333    let mut order: Vec<usize> = (0..carriers.len()).collect();
334    order.sort_by_key(|&piece| ends_of[piece]);
335    for (rank, &piece) in order.iter().enumerate() {
336        renumber[piece] = rank;
337    }
338    let pieces: Vec<Piece> = order
339        .iter()
340        .map(|&piece| Piece {
341            ends: [VertexId(ends_of[piece].0), VertexId(ends_of[piece].1)],
342            segments: carriers[piece].iter().copied().collect(),
343        })
344        .collect();
345    let mut paths = vec![Vec::new(); segments.len()];
346    for (segment, piece, forward) in on_segments {
347        paths[segment].push((PieceId(renumber[piece]), forward));
348    }
349    let mut incident = vec![Vec::new(); vertices.len()];
350    for (index, piece) in pieces.iter().enumerate() {
351        incident[piece.ends[0].0].push((PieceId(index), true));
352        incident[piece.ends[1].0].push((PieceId(index), false));
353    }
354    Ok(Arrangement {
355        vertices,
356        pieces,
357        paths,
358        crossings,
359        incident,
360    })
361}
362
363/// Examine the pairs of segments whose boxes overlap: join the proper
364/// crossings and the ends that lie inside the other segment.
365fn cut_segments(
366    segments: &[Segment],
367    cuts: &mut [Cut],
368    points: &mut ExactPoints,
369    local: &mut Local,
370    (budget, declared_bytes): (&mut PairBudget, usize),
371    options: &ExecutionOptions,
372) -> GeomResult<()> {
373    let mut order: Vec<usize> = (0..cuts.len()).collect();
374    order.sort_by(|&left, &right| {
375        cuts[left].bounds[0]
376            .total_cmp(&cuts[right].bounds[0])
377            .then(left.cmp(&right))
378    });
379    // The boxes in sweep order, to scan them without looking segments up.
380    let boxes: Vec<[f64; 4]> = order.iter().map(|&index| cuts[index].bounds).collect();
381    let mut charged = 0usize;
382    for (rank, &first) in order.iter().enumerate() {
383        options.check_cancelled()?;
384        for (&second, held) in order[rank + 1..].iter().zip(&boxes[rank + 1..]) {
385            if held[0] > boxes[rank][2] {
386                break;
387            }
388            if !boxes_overlap(&boxes[rank], held) {
389                continue;
390            }
391            options.check_cancelled()?;
392            budget.count()?;
393            let (low, high) = (first.min(second), first.max(second));
394            let [a, b] = segments[low].ends;
395            let [c, d] = segments[high].ends;
396            // The same segment twice is one piece and cuts nothing.
397            if [a, b] == [c, d] || [a, b] == [d, c] {
398                continue;
399            }
400            // Segments of one point meet only there, unless they lie on one
401            // line: then the nearer of the two other ends is inside the
402            // longer segment and splits it.
403            if let Some((shared, near, far)) = shared_end([a, b], [c, d]) {
404                let (low_ends, high_ends) = (cuts[low].ends, cuts[high].ends);
405                let local = |end: PointId| {
406                    let (own, cut) = if [a, b].contains(&end) {
407                        ([a, b], low_ends)
408                    } else {
409                        ([c, d], high_ends)
410                    };
411                    cut[usize::from(own[1] == end)]
412                };
413                if points.orient(shared, near, far)? == Sign::Zero {
414                    // `near` and `far` belong to different segments.
415                    if points.strictly_between(shared, near, far)? {
416                        let longer = if [a, b].contains(&far) { low } else { high };
417                        charge(options, declared_bytes, &mut charged, SPLIT_SCRATCH_BYTES)?;
418                        cuts[longer].splits.push(local(near));
419                    } else if points.strictly_between(shared, far, near)? {
420                        let longer = if [a, b].contains(&near) { low } else { high };
421                        charge(options, declared_bytes, &mut charged, SPLIT_SCRATCH_BYTES)?;
422                        cuts[longer].splits.push(local(far));
423                    }
424                }
425                continue;
426            }
427            // Both ends of one segment strictly on one side of the other's
428            // line: they neither cross nor touch, and most pairs end here.
429            let (past_c, past_d) = (points.orient(a, b, c)?, points.orient(a, b, d)?);
430            if past_c == past_d && past_c != Sign::Zero {
431                continue;
432            }
433            let (past_a, past_b) = (points.orient(c, d, a)?, points.orient(c, d, b)?);
434            if past_a == past_b && past_a != Sign::Zero {
435                continue;
436            }
437            let around = [past_c, past_d, past_a, past_b];
438            let proper = around.iter().all(|&sign| sign != Sign::Zero)
439                && around[0] != around[1]
440                && around[2] != around[3];
441            if proper {
442                charge(
443                    options,
444                    declared_bytes,
445                    &mut charged,
446                    CROSSING_SCRATCH_BYTES,
447                )?;
448                let ends = [a, b, c, d];
449                let point = points.segment_crossing(ends)?;
450                let new = local.push(point, Label::Crossing([low, high], ends));
451                cuts[low].splits.push(new);
452                cuts[high].splits.push(new);
453            } else {
454                let (low_ends, high_ends) = (cuts[low].ends, cuts[high].ends);
455                if around[0] == Sign::Zero && points.strictly_between(a, c, b)? {
456                    charge(options, declared_bytes, &mut charged, SPLIT_SCRATCH_BYTES)?;
457                    cuts[low].splits.push(high_ends[0]);
458                }
459                if around[1] == Sign::Zero && points.strictly_between(a, d, b)? {
460                    charge(options, declared_bytes, &mut charged, SPLIT_SCRATCH_BYTES)?;
461                    cuts[low].splits.push(high_ends[1]);
462                }
463                if around[2] == Sign::Zero && points.strictly_between(c, a, d)? {
464                    charge(options, declared_bytes, &mut charged, SPLIT_SCRATCH_BYTES)?;
465                    cuts[high].splits.push(low_ends[0]);
466                }
467                if around[3] == Sign::Zero && points.strictly_between(c, b, d)? {
468                    charge(options, declared_bytes, &mut charged, SPLIT_SCRATCH_BYTES)?;
469                    cuts[high].splits.push(low_ends[1]);
470                }
471            }
472        }
473    }
474    Ok(())
475}
476
477/// The end two different segments share, with the other end of each (first
478/// segment's first), or `None` when they share none.
479fn shared_end(first: [PointId; 2], second: [PointId; 2]) -> Option<(PointId, PointId, PointId)> {
480    for (i, &end) in first.iter().enumerate() {
481        if let Some(j) = second.iter().position(|&other| other == end) {
482            return Some((end, first[1 - i], second[1 - j]));
483        }
484    }
485    None
486}
487
488/// Refuse the next charge when the scratch declared per element and the
489/// crossings and splits charged so far no longer fit the memory budget.
490pub(crate) fn charge(
491    options: &ExecutionOptions,
492    declared_bytes: usize,
493    charged: &mut usize,
494    cost: usize,
495) -> GeomResult<()> {
496    let Some(budget) = options.memory_budget_bytes() else {
497        return Ok(());
498    };
499    let needed = charged
500        .checked_add(cost)
501        .and_then(|total| declared_bytes.checked_add(total));
502    match needed {
503        Some(needed) if needed <= budget => {
504            *charged += cost;
505            Ok(())
506        }
507        _ => Err(GeomError::BudgetExceeded { resource: "memory" }),
508    }
509}
510
511/// Put the points of every segment in order along it, and join equal ones.
512fn order_cuts(
513    cuts: &mut [Cut],
514    points: &ExactPoints,
515    local: &mut Local,
516    options: &ExecutionOptions,
517) -> GeomResult<()> {
518    for cut in cuts {
519        options.check_cancelled()?;
520        let [first, last] = cut.ends.map(|end| local.exact[end]);
521        let axis = points.dominant(first, last)?;
522        let descending = points.compare(first, last, axis)? == std::cmp::Ordering::Greater;
523        let mut inner = std::mem::take(&mut cut.splits);
524        inner.sort_unstable();
525        inner.dedup();
526        sort_exact(&mut inner, |left, right| {
527            let order = points.compare(local.exact[left], local.exact[right], axis)?;
528            Ok(if descending { order.reverse() } else { order })
529        })?;
530        cut.path = std::iter::once(cut.ends[0])
531            .chain(inner)
532            .chain(std::iter::once(cut.ends[1]))
533            .collect();
534        for pair in cut.path.windows(2) {
535            if points.same(local.exact[pair[0]], local.exact[pair[1]])? {
536                local.union(pair[0], pair[1]);
537            }
538        }
539    }
540    Ok(())
541}
542
543/// Number the groups of equal points: those with input points first, by their
544/// smallest one, then the crossings, by their smallest pair of segments.
545fn number_vertices(local: &mut Local) -> (Vec<Vertex>, BTreeMap<usize, usize>) {
546    struct Group {
547        inputs: Vec<PointId>,
548        /// The smallest crossing label of the group, with its point.
549        crossing: Option<(Crossing, usize)>,
550        /// Any point of the group.
551        any: usize,
552    }
553    let mut groups = BTreeMap::<usize, Group>::new();
554    for point in 0..local.parent.len() {
555        let root = local.find(point);
556        let group = groups.entry(root).or_insert(Group {
557            inputs: Vec::new(),
558            crossing: None,
559            any: point,
560        });
561        match local.labels[point] {
562            Label::Input => group.inputs.push(local.exact[point]),
563            Label::Crossing(segments, ends) => {
564                if group
565                    .crossing
566                    .is_none_or(|(known, _)| segments < known.segments)
567                {
568                    group.crossing = Some((Crossing { segments, ends }, point));
569                }
570            }
571        }
572    }
573    let mut ordered: Vec<(usize, Group)> = groups.into_iter().collect();
574    for (_, group) in &mut ordered {
575        group.inputs.sort_unstable();
576    }
577    ordered.sort_by(
578        |(_, left), (_, right)| match (left.inputs.first(), right.inputs.first()) {
579            (Some(a), Some(b)) => a.cmp(b),
580            (Some(_), None) => std::cmp::Ordering::Less,
581            (None, Some(_)) => std::cmp::Ordering::Greater,
582            (None, None) => left
583                .crossing
584                .map(|c| c.0.segments)
585                .cmp(&right.crossing.map(|c| c.0.segments)),
586        },
587    );
588    let mut vertex_of = BTreeMap::new();
589    let mut vertices = Vec::with_capacity(ordered.len());
590    for (number, (root, group)) in ordered.into_iter().enumerate() {
591        vertex_of.insert(root, number);
592        let (crossing, point) = match (group.inputs.is_empty(), group.crossing) {
593            (true, Some((crossing, point))) => (Some(crossing), point),
594            _ => (None, group.any),
595        };
596        vertices.push(Vertex {
597            point: local.exact[point],
598            inputs: group.inputs,
599            crossing,
600        });
601    }
602    (vertices, vertex_of)
603}
604
605impl Arrangement {
606    /// The pieces at `vertex` in counter-clockwise order in the points'
607    /// frame, each with whether it leaves the vertex at `ends[0]`.
608    ///
609    /// Both `orient` and the half-plane split are in the stored axes up to
610    /// the frame sign, which reverses the order within a half but keeps the
611    /// cyclic order. Two pieces that leave a vertex the same way overlap,
612    /// which `arrange` has cut apart, so they are an inconsistency.
613    pub fn fan(&self, points: &ExactPoints, vertex: VertexId) -> GeomResult<Vec<(PieceId, bool)>> {
614        let here = self.vertices[vertex.0].point;
615        let mut spokes = Vec::with_capacity(self.incident[vertex.0].len());
616        for &(piece, forward) in &self.incident[vertex.0] {
617            let ends = self.pieces[piece.0].ends;
618            let to = self.vertices[ends[usize::from(forward)].0].point;
619            let upper = points.compare(to, here, 1)? == std::cmp::Ordering::Greater
620                || (points.compare(to, here, 1)? == std::cmp::Ordering::Equal
621                    && points.compare(to, here, 0)? == std::cmp::Ordering::Greater);
622            spokes.push((piece, forward, to, u8::from(!upper)));
623        }
624        sort_exact(&mut spokes, |left, right| {
625            if left.3 != right.3 {
626                return Ok(left.3.cmp(&right.3));
627            }
628            match points.orient(here, left.2, right.2)? {
629                Sign::Positive => Ok(std::cmp::Ordering::Less),
630                Sign::Negative => Ok(std::cmp::Ordering::Greater),
631                _ => Err(inconsistent()),
632            }
633        })?;
634        Ok(spokes
635            .into_iter()
636            .map(|(piece, forward, _, _)| (piece, forward))
637            .collect())
638    }
639}
axiolid_exact_arrangement/arrangement.rs

axiolid_exact_arrangement/
arrangement.rs
