Skip to main content

brep_kernel/intersect/
arrangement.rs

1use serde::{Deserialize, Serialize};
2
3#[derive(Clone, Copy, Debug, Deserialize, Serialize)]
4pub struct Vec2 {
5    pub x: f64,
6    pub y: f64,
7}
8
9impl Vec2 {
10    pub(crate) fn sub(self, other: Self) -> Self {
11        Self {
12            x: self.x - other.x,
13            y: self.y - other.y,
14        }
15    }
16    pub(crate) fn add(self, other: Self) -> Self {
17        Self {
18            x: self.x + other.x,
19            y: self.y + other.y,
20        }
21    }
22    pub(crate) fn scale(self, factor: f64) -> Self {
23        Self {
24            x: self.x * factor,
25            y: self.y * factor,
26        }
27    }
28    fn distance(self, other: Self) -> f64 {
29        self.sub(other).length()
30    }
31    pub(crate) fn length(self) -> f64 {
32        (self.x * self.x + self.y * self.y).sqrt()
33    }
34    fn cross(self, other: Self) -> f64 {
35        self.x * other.y - self.y * other.x
36    }
37    pub(crate) fn dot(self, other: Self) -> f64 {
38        self.x * other.x + self.y * other.y
39    }
40}
41
42#[derive(Clone, Debug, Deserialize, Serialize)]
43pub struct Segment2 {
44    pub a: Vec2,
45    pub b: Vec2,
46    #[serde(default)]
47    pub tag: serde_json::Value,
48}
49
50#[derive(Clone, Debug, Serialize)]
51pub struct ArrangementPiece {
52    pub a: Vec2,
53    pub b: Vec2,
54    pub parent: Segment2,
55}
56
57#[derive(Clone, Debug, Serialize)]
58pub struct CycleUse {
59    pub piece: ArrangementPiece,
60    pub forward: bool,
61}
62
63#[derive(Clone, Debug, Serialize)]
64pub struct ArrangementRegion {
65    pub outer: Vec<CycleUse>,
66    pub holes: Vec<Vec<CycleUse>>,
67    pub area: f64,
68}
69
70#[derive(Clone)]
71struct Node {
72    point: Vec2,
73    outgoing: Vec<usize>,
74}
75
76fn tolerance_parameter(segment: &Segment2, tolerance: f64) -> f64 {
77    let length = segment.a.distance(segment.b);
78    if length <= tolerance {
79        0.5
80    } else {
81        tolerance / length
82    }
83}
84
85fn interpolate(a: Vec2, b: Vec2, parameter: f64) -> Vec2 {
86    Vec2 {
87        x: a.x + (b.x - a.x) * parameter,
88        y: a.y + (b.y - a.y) * parameter,
89    }
90}
91
92pub fn segment_intersection(
93    a1: Vec2,
94    b1: Vec2,
95    a2: Vec2,
96    b2: Vec2,
97    tolerance: f64,
98) -> Option<[f64; 2]> {
99    let d1 = b1.sub(a1);
100    let d2 = b2.sub(a2);
101    let denominator = d1.cross(d2);
102    let length1 = d1.length();
103    let length2 = d2.length();
104    if denominator.abs() <= 1e-12 * length1 * length2 {
105        return None;
106    }
107    let offset = a2.sub(a1);
108    let t1 = offset.cross(d2) / denominator;
109    let t2 = offset.cross(d1) / denominator;
110    let epsilon1 = tolerance / length1.max(tolerance);
111    let epsilon2 = tolerance / length2.max(tolerance);
112    if t1 < -epsilon1 || t1 > 1.0 + epsilon1 || t2 < -epsilon2 || t2 > 1.0 + epsilon2 {
113        None
114    } else {
115        Some([t1.clamp(0.0, 1.0), t2.clamp(0.0, 1.0)])
116    }
117}
118
119#[derive(Clone, Copy, PartialEq)]
120enum PolygonClass {
121    Inside,
122    Outside,
123    Boundary,
124}
125
126fn point_segment_distance(point: Vec2, a: Vec2, b: Vec2) -> f64 {
127    let segment = b.sub(a);
128    let length_squared = segment.dot(segment);
129    if length_squared <= 1e-300 {
130        return point.distance(a);
131    }
132    let parameter = (point.sub(a).dot(segment) / length_squared).clamp(0.0, 1.0);
133    point.distance(Vec2 {
134        x: a.x + segment.x * parameter,
135        y: a.y + segment.y * parameter,
136    })
137}
138
139fn point_in_polygon_class(point: Vec2, polygon: &[Vec2], tolerance: f64) -> PolygonClass {
140    for index in 0..polygon.len() {
141        if point_segment_distance(point, polygon[index], polygon[(index + 1) % polygon.len()])
142            <= tolerance
143        {
144            return PolygonClass::Boundary;
145        }
146    }
147    let mut inside = false;
148    for index in 0..polygon.len() {
149        let a = polygon[index];
150        let b = polygon[(index + 1) % polygon.len()];
151        if (a.y > point.y) != (b.y > point.y) {
152            let crossing = a.x + (point.y - a.y) / (b.y - a.y) * (b.x - a.x);
153            if crossing > point.x {
154                inside = !inside;
155            }
156        }
157    }
158    if inside {
159        PolygonClass::Inside
160    } else {
161        PolygonClass::Outside
162    }
163}
164
165pub fn point_in_polygon(point: Vec2, polygon: &[Vec2], tolerance: f64) -> &'static str {
166    match point_in_polygon_class(point, polygon, tolerance) {
167        PolygonClass::Inside => "in",
168        PolygonClass::Outside => "out",
169        PolygonClass::Boundary => "boundary",
170    }
171}
172
173#[derive(Clone)]
174struct Cycle {
175    uses: Vec<CycleUse>,
176    area: f64,
177    node_indices: Vec<usize>,
178}
179
180pub fn arrange_segments(
181    segments: &[Segment2],
182    tolerance: f64,
183) -> Result<Vec<ArrangementRegion>, String> {
184    let mut cuts = vec![Vec::<f64>::new(); segments.len()];
185    for first in 0..segments.len() {
186        for second in first + 1..segments.len() {
187            let Some([first_parameter, second_parameter]) = segment_intersection(
188                segments[first].a,
189                segments[first].b,
190                segments[second].a,
191                segments[second].b,
192                tolerance,
193            ) else {
194                continue;
195            };
196            let first_tolerance = tolerance_parameter(&segments[first], tolerance);
197            let second_tolerance = tolerance_parameter(&segments[second], tolerance);
198            if first_parameter > first_tolerance && first_parameter < 1.0 - first_tolerance {
199                cuts[first].push(first_parameter);
200            }
201            if second_parameter > second_tolerance && second_parameter < 1.0 - second_tolerance {
202                cuts[second].push(second_parameter);
203            }
204        }
205    }
206
207    let mut pieces = Vec::new();
208    for (index, segment) in segments.iter().enumerate() {
209        let mut parameters = vec![0.0];
210        cuts[index].sort_by(f64::total_cmp);
211        parameters.extend(cuts[index].iter().copied());
212        parameters.push(1.0);
213        let parameter_tolerance = tolerance_parameter(segment, tolerance);
214        let mut deduplicated: Vec<f64> = Vec::new();
215        for parameter in parameters {
216            if deduplicated
217                .last()
218                .is_none_or(|previous| parameter - *previous > parameter_tolerance)
219            {
220                deduplicated.push(parameter);
221            } else if parameter == 1.0 {
222                *deduplicated.last_mut().unwrap() = 1.0;
223            }
224        }
225        for pair in deduplicated.windows(2) {
226            let a = interpolate(segment.a, segment.b, pair[0]);
227            let b = interpolate(segment.a, segment.b, pair[1]);
228            if a.distance(b) <= tolerance {
229                continue;
230            }
231            pieces.push(ArrangementPiece {
232                a,
233                b,
234                parent: segment.clone(),
235            });
236        }
237    }
238
239    let mut nodes: Vec<Node> = Vec::new();
240    let find_node = |point: Vec2, nodes: &mut Vec<Node>| {
241        if let Some(index) = nodes
242            .iter()
243            .position(|node| node.point.distance(point) <= tolerance)
244        {
245            index
246        } else {
247            nodes.push(Node {
248                point,
249                outgoing: Vec::new(),
250            });
251            nodes.len() - 1
252        }
253    };
254    let mut piece_start = Vec::with_capacity(pieces.len());
255    let mut piece_end = Vec::with_capacity(pieces.len());
256    for piece in &pieces {
257        piece_start.push(find_node(piece.a, &mut nodes));
258        piece_end.push(find_node(piece.b, &mut nodes));
259    }
260
261    let mut alive = vec![true; pieces.len()];
262    loop {
263        let mut pruned = false;
264        let mut degree = vec![0usize; nodes.len()];
265        for index in 0..pieces.len() {
266            if !alive[index] {
267                continue;
268            }
269            if piece_start[index] == piece_end[index] {
270                alive[index] = false;
271                pruned = true;
272                continue;
273            }
274            degree[piece_start[index]] += 1;
275            degree[piece_end[index]] += 1;
276        }
277        for index in 0..pieces.len() {
278            if alive[index] && (degree[piece_start[index]] == 1 || degree[piece_end[index]] == 1) {
279                alive[index] = false;
280                pruned = true;
281            }
282        }
283        if !pruned {
284            break;
285        }
286    }
287
288    let half_tail = |half_edge: usize| {
289        if half_edge % 2 == 0 {
290            piece_start[half_edge >> 1]
291        } else {
292            piece_end[half_edge >> 1]
293        }
294    };
295    let half_head = |half_edge: usize| {
296        if half_edge % 2 == 0 {
297            piece_end[half_edge >> 1]
298        } else {
299            piece_start[half_edge >> 1]
300        }
301    };
302    for node in &mut nodes {
303        node.outgoing.clear();
304    }
305    for index in 0..pieces.len() {
306        if !alive[index] {
307            continue;
308        }
309        nodes[piece_start[index]].outgoing.push(index * 2);
310        nodes[piece_end[index]].outgoing.push(index * 2 + 1);
311    }
312    let angles: Vec<f64> = (0..pieces.len() * 2)
313        .map(|half_edge| {
314            let tail = nodes[half_tail(half_edge)].point;
315            let head = nodes[half_head(half_edge)].point;
316            (head.y - tail.y).atan2(head.x - tail.x)
317        })
318        .collect();
319    for node in &mut nodes {
320        node.outgoing
321            .sort_by(|a, b| angles[*a].total_cmp(&angles[*b]));
322    }
323
324    let mut visited = rustc_hash::FxHashSet::default();
325    let mut positive = Vec::new();
326    let mut negative = Vec::new();
327    for index in 0..pieces.len() {
328        if !alive[index] {
329            continue;
330        }
331        for start in [index * 2, index * 2 + 1] {
332            if visited.contains(&start) {
333                continue;
334            }
335            let mut uses = Vec::new();
336            let mut node_indices = Vec::new();
337            let mut area = 0.0;
338            let mut half_edge = start;
339            let mut guard = 0;
340            loop {
341                visited.insert(half_edge);
342                let piece_index = half_edge >> 1;
343                uses.push(CycleUse {
344                    piece: pieces[piece_index].clone(),
345                    forward: half_edge % 2 == 0,
346                });
347                let tail_index = half_tail(half_edge);
348                let head_index = half_head(half_edge);
349                let tail = nodes[tail_index].point;
350                let head = nodes[head_index].point;
351                area += 0.5 * (tail.x * head.y - head.x * tail.y);
352                node_indices.push(tail_index);
353                let reverse = half_edge ^ 1;
354                let outgoing = &nodes[head_index].outgoing;
355                let reverse_index = outgoing
356                    .iter()
357                    .position(|candidate| *candidate == reverse)
358                    .ok_or_else(|| "arrangeSegments: reverse half-edge missing".to_string())?;
359                half_edge = outgoing[(reverse_index + outgoing.len() - 1) % outgoing.len()];
360                guard += 1;
361                if guard > pieces.len() * 4 + 8 {
362                    return Err("arrangeSegments: face walk did not terminate".into());
363                }
364                if half_edge == start {
365                    break;
366                }
367            }
368            let cycle = Cycle {
369                uses,
370                area,
371                node_indices,
372            };
373            if area > tolerance * tolerance {
374                positive.push(cycle);
375            } else if area < -tolerance * tolerance {
376                negative.push(cycle);
377            }
378        }
379    }
380
381    positive.sort_by(|a, b| a.area.total_cmp(&b.area));
382    let mut regions: Vec<(ArrangementRegion, Vec<Vec2>)> = positive
383        .into_iter()
384        .map(|cycle| {
385            let polygon = cycle
386                .node_indices
387                .iter()
388                .map(|index| nodes[*index].point)
389                .collect();
390            (
391                ArrangementRegion {
392                    outer: cycle.uses,
393                    holes: Vec::new(),
394                    area: cycle.area,
395                },
396                polygon,
397            )
398        })
399        .collect();
400    for cycle in negative {
401        for (region, polygon) in &mut regions {
402            let mut inside = false;
403            let mut on_boundary = false;
404            for node_index in &cycle.node_indices {
405                match point_in_polygon_class(nodes[*node_index].point, polygon, tolerance) {
406                    PolygonClass::Boundary => {
407                        on_boundary = true;
408                        continue;
409                    }
410                    PolygonClass::Inside => {
411                        inside = true;
412                        on_boundary = false;
413                        break;
414                    }
415                    PolygonClass::Outside => {
416                        inside = false;
417                        on_boundary = false;
418                        break;
419                    }
420                }
421            }
422            if on_boundary {
423                continue;
424            }
425            if inside {
426                region.holes.push(cycle.uses.clone());
427                region.area += cycle.area;
428                break;
429            }
430        }
431    }
432    Ok(regions.into_iter().map(|(region, _)| region).collect())
433}
434
435#[cfg(test)]
436mod tests {
437    use super::*;
438    use serde_json::json;
439
440    fn segment(a: [f64; 2], b: [f64; 2], tag: &str) -> Segment2 {
441        Segment2 {
442            a: Vec2 { x: a[0], y: a[1] },
443            b: Vec2 { x: b[0], y: b[1] },
444            tag: json!(tag),
445        }
446    }
447
448    #[test]
449    fn crossing_cut_splits_square_into_two_regions() {
450        let segments = vec![
451            segment([0.0, 0.0], [4.0, 0.0], "bottom"),
452            segment([4.0, 0.0], [4.0, 3.0], "right"),
453            segment([4.0, 3.0], [0.0, 3.0], "top"),
454            segment([0.0, 3.0], [0.0, 0.0], "left"),
455            segment([2.0, 0.0], [2.0, 3.0], "cut"),
456        ];
457        let regions = arrange_segments(&segments, 1e-7).unwrap();
458        assert_eq!(regions.len(), 2);
459        assert!(regions
460            .iter()
461            .all(|region| (region.area - 6.0).abs() < 1e-9));
462    }
463
464    #[test]
465    fn nested_cycles_assign_hole() {
466        let segments = vec![
467            segment([0.0, 0.0], [6.0, 0.0], "outer"),
468            segment([6.0, 0.0], [6.0, 6.0], "outer"),
469            segment([6.0, 6.0], [0.0, 6.0], "outer"),
470            segment([0.0, 6.0], [0.0, 0.0], "outer"),
471            segment([2.0, 2.0], [2.0, 4.0], "inner"),
472            segment([2.0, 4.0], [4.0, 4.0], "inner"),
473            segment([4.0, 4.0], [4.0, 2.0], "inner"),
474            segment([4.0, 2.0], [2.0, 2.0], "inner"),
475        ];
476        let regions = arrange_segments(&segments, 1e-7).unwrap();
477        assert!(regions
478            .iter()
479            .any(|region| { region.holes.len() == 1 && (region.area - 32.0).abs() < 1e-9 }));
480    }
481}