nodedb_spatial/predicates/
edge.rs

1// SPDX-License-Identifier: Apache-2.0
2
3//! Low-level segment and edge geometry primitives.
4//!
5//! Used by ST_Contains, ST_Intersects, and ST_Distance. All operations
6//! are planar (coordinate-space), not spherical. This is correct for the
7//! small-area geometries typical of spatial predicates (city blocks, not
8//! hemispheres). For large-area work, haversine handles the globe.
9
10/// Orientation of three points (collinear, clockwise, counter-clockwise).
11/// Uses the cross product of vectors (p→q) and (p→r).
12#[derive(Debug, Clone, Copy, PartialEq, Eq)]
13pub enum Orientation {
14    Collinear,
15    Clockwise,
16    CounterClockwise,
17}
18
19/// Compute orientation of ordered triplet (p, q, r).
20pub fn orientation(p: [f64; 2], q: [f64; 2], r: [f64; 2]) -> Orientation {
21    let val = (q[1] - p[1]) * (r[0] - q[0]) - (q[0] - p[0]) * (r[1] - q[1]);
22    // Use epsilon for floating-point tolerance.
23    if val.abs() < 1e-12 {
24        Orientation::Collinear
25    } else if val > 0.0 {
26        Orientation::Clockwise
27    } else {
28        Orientation::CounterClockwise
29    }
30}
31
32/// Whether point q lies on segment p-r (given that p, q, r are collinear).
33pub fn on_segment(p: [f64; 2], q: [f64; 2], r: [f64; 2]) -> bool {
34    q[0] <= p[0].max(r[0])
35        && q[0] >= p[0].min(r[0])
36        && q[1] <= p[1].max(r[1])
37        && q[1] >= p[1].min(r[1])
38}
39
40/// Whether two segments (p1-q1) and (p2-q2) intersect.
41///
42/// Uses the standard orientation-based algorithm. Handles collinear
43/// overlapping segments correctly.
44pub fn segments_intersect(p1: [f64; 2], q1: [f64; 2], p2: [f64; 2], q2: [f64; 2]) -> bool {
45    let o1 = orientation(p1, q1, p2);
46    let o2 = orientation(p1, q1, q2);
47    let o3 = orientation(p2, q2, p1);
48    let o4 = orientation(p2, q2, q1);
49
50    // General case: different orientations mean crossing.
51    if o1 != o2 && o3 != o4 {
52        return true;
53    }
54
55    // Collinear special cases: check if endpoints lie on the other segment.
56    if o1 == Orientation::Collinear && on_segment(p1, p2, q1) {
57        return true;
58    }
59    if o2 == Orientation::Collinear && on_segment(p1, q2, q1) {
60        return true;
61    }
62    if o3 == Orientation::Collinear && on_segment(p2, p1, q2) {
63        return true;
64    }
65    if o4 == Orientation::Collinear && on_segment(p2, q1, q2) {
66        return true;
67    }
68
69    false
70}
71
72/// Whether a point lies exactly on a line segment (within epsilon tolerance).
73pub fn point_on_segment(pt: [f64; 2], seg_a: [f64; 2], seg_b: [f64; 2]) -> bool {
74    // Check collinearity via cross product.
75    let cross =
76        (pt[0] - seg_a[0]) * (seg_b[1] - seg_a[1]) - (pt[1] - seg_a[1]) * (seg_b[0] - seg_a[0]);
77    if cross.abs() > 1e-10 {
78        return false;
79    }
80    // Check that pt is within the segment's bounding box.
81    pt[0] >= seg_a[0].min(seg_b[0]) - 1e-10
82        && pt[0] <= seg_a[0].max(seg_b[0]) + 1e-10
83        && pt[1] >= seg_a[1].min(seg_b[1]) - 1e-10
84        && pt[1] <= seg_a[1].max(seg_b[1]) + 1e-10
85}
86
87/// Whether a point lies on any edge of a polygon ring.
88pub fn point_on_ring_boundary(pt: [f64; 2], ring: &[[f64; 2]]) -> bool {
89    if ring.len() < 2 {
90        return false;
91    }
92    for i in 0..ring.len() - 1 {
93        if point_on_segment(pt, ring[i], ring[i + 1]) {
94            return true;
95        }
96    }
97    // Check closing segment if ring isn't explicitly closed.
98    if ring.first() != ring.last()
99        && let (Some(&first), Some(&last)) = (ring.first(), ring.last())
100        && point_on_segment(pt, last, first)
101    {
102        return true;
103    }
104    false
105}
106
107/// Minimum squared distance from a point to a line segment.
108///
109/// Returns the squared distance (avoid sqrt for comparison purposes).
110pub fn point_to_segment_dist_sq(pt: [f64; 2], seg_a: [f64; 2], seg_b: [f64; 2]) -> f64 {
111    let dx = seg_b[0] - seg_a[0];
112    let dy = seg_b[1] - seg_a[1];
113    let len_sq = dx * dx + dy * dy;
114
115    if len_sq < 1e-20 {
116        // Degenerate segment (zero length) — distance to point.
117        let ex = pt[0] - seg_a[0];
118        let ey = pt[1] - seg_a[1];
119        return ex * ex + ey * ey;
120    }
121
122    // Project pt onto the line, clamped to [0, 1].
123    let t = ((pt[0] - seg_a[0]) * dx + (pt[1] - seg_a[1]) * dy) / len_sq;
124    let t = t.clamp(0.0, 1.0);
125
126    let proj_x = seg_a[0] + t * dx;
127    let proj_y = seg_a[1] + t * dy;
128
129    let ex = pt[0] - proj_x;
130    let ey = pt[1] - proj_y;
131    ex * ex + ey * ey
132}
133
134/// Minimum squared distance between two line segments.
135pub fn segment_to_segment_dist_sq(a1: [f64; 2], a2: [f64; 2], b1: [f64; 2], b2: [f64; 2]) -> f64 {
136    if segments_intersect(a1, a2, b1, b2) {
137        return 0.0;
138    }
139    // Min of all endpoint-to-segment distances.
140    let d1 = point_to_segment_dist_sq(a1, b1, b2);
141    let d2 = point_to_segment_dist_sq(a2, b1, b2);
142    let d3 = point_to_segment_dist_sq(b1, a1, a2);
143    let d4 = point_to_segment_dist_sq(b2, a1, a2);
144    d1.min(d2).min(d3).min(d4)
145}
146
147/// Extract all edges from a polygon ring as segment pairs.
148pub fn ring_edges(ring: &[[f64; 2]]) -> Vec<([f64; 2], [f64; 2])> {
149    if ring.len() < 2 {
150        return Vec::new();
151    }
152    // no-governor: hot-path ring edge extraction; ring.len() = vertex count, predicate path
153    let mut edges = Vec::with_capacity(ring.len());
154    for i in 0..ring.len() - 1 {
155        edges.push((ring[i], ring[i + 1]));
156    }
157    // Close the ring if not explicitly closed.
158    if ring.first() != ring.last()
159        && let (Some(&first), Some(&last)) = (ring.first(), ring.last())
160    {
161        edges.push((last, first));
162    }
163    edges
164}
165
166#[cfg(test)]
167mod tests {
168    use super::*;
169
170    #[test]
171    fn segments_cross() {
172        assert!(segments_intersect(
173            [0.0, 0.0],
174            [10.0, 10.0],
175            [0.0, 10.0],
176            [10.0, 0.0]
177        ));
178    }
179
180    #[test]
181    fn segments_parallel_no_cross() {
182        assert!(!segments_intersect(
183            [0.0, 0.0],
184            [10.0, 0.0],
185            [0.0, 1.0],
186            [10.0, 1.0]
187        ));
188    }
189
190    #[test]
191    fn segments_collinear_overlap() {
192        assert!(segments_intersect(
193            [0.0, 0.0],
194            [5.0, 0.0],
195            [3.0, 0.0],
196            [8.0, 0.0]
197        ));
198    }
199
200    #[test]
201    fn segments_collinear_no_overlap() {
202        assert!(!segments_intersect(
203            [0.0, 0.0],
204            [2.0, 0.0],
205            [3.0, 0.0],
206            [5.0, 0.0]
207        ));
208    }
209
210    #[test]
211    fn segments_endpoint_touch() {
212        assert!(segments_intersect(
213            [0.0, 0.0],
214            [5.0, 5.0],
215            [5.0, 5.0],
216            [10.0, 0.0]
217        ));
218    }
219
220    #[test]
221    fn point_on_segment_middle() {
222        assert!(point_on_segment([5.0, 5.0], [0.0, 0.0], [10.0, 10.0]));
223    }
224
225    #[test]
226    fn point_on_segment_endpoint() {
227        assert!(point_on_segment([0.0, 0.0], [0.0, 0.0], [10.0, 10.0]));
228    }
229
230    #[test]
231    fn point_off_segment() {
232        assert!(!point_on_segment([5.0, 6.0], [0.0, 0.0], [10.0, 10.0]));
233    }
234
235    #[test]
236    fn point_on_ring_edge() {
237        let ring = vec![
238            [0.0, 0.0],
239            [10.0, 0.0],
240            [10.0, 10.0],
241            [0.0, 10.0],
242            [0.0, 0.0],
243        ];
244        assert!(point_on_ring_boundary([5.0, 0.0], &ring)); // bottom edge
245        assert!(point_on_ring_boundary([0.0, 5.0], &ring)); // left edge
246        assert!(!point_on_ring_boundary([5.0, 5.0], &ring)); // interior
247    }
248
249    #[test]
250    fn point_on_ring_vertex() {
251        let ring = vec![
252            [0.0, 0.0],
253            [10.0, 0.0],
254            [10.0, 10.0],
255            [0.0, 10.0],
256            [0.0, 0.0],
257        ];
258        assert!(point_on_ring_boundary([0.0, 0.0], &ring));
259        assert!(point_on_ring_boundary([10.0, 10.0], &ring));
260    }
261
262    #[test]
263    fn point_to_segment_perpendicular() {
264        let d = point_to_segment_dist_sq([5.0, 1.0], [0.0, 0.0], [10.0, 0.0]);
265        assert!((d - 1.0).abs() < 1e-10); // 1 unit away
266    }
267
268    #[test]
269    fn point_to_segment_endpoint() {
270        let d = point_to_segment_dist_sq([12.0, 0.0], [0.0, 0.0], [10.0, 0.0]);
271        assert!((d - 4.0).abs() < 1e-10); // 2 units past endpoint
272    }
273
274    #[test]
275    fn segment_to_segment_parallel() {
276        let d = segment_to_segment_dist_sq([0.0, 0.0], [10.0, 0.0], [0.0, 3.0], [10.0, 3.0]);
277        assert!((d - 9.0).abs() < 1e-10);
278    }
279
280    #[test]
281    fn segment_to_segment_crossing() {
282        let d = segment_to_segment_dist_sq([0.0, 0.0], [10.0, 10.0], [0.0, 10.0], [10.0, 0.0]);
283        assert!(d < 1e-10);
284    }
285
286    #[test]
287    fn ring_edges_closed() {
288        let ring = vec![
289            [0.0, 0.0],
290            [10.0, 0.0],
291            [10.0, 10.0],
292            [0.0, 10.0],
293            [0.0, 0.0],
294        ];
295        let edges = ring_edges(&ring);
296        assert_eq!(edges.len(), 4);
297    }
298}
nodedb_spatial/predicates/edge.rs

nodedb_spatial/predicates/
edge.rs
