qualia_core_db/specialized_libs/computational_geometry/
hull.rs

1use core::cmp::Ordering;
2
3use crate::tensor::Tensor10D;
4
5use super::kernel::{FilteredF64Kernel, GeometryKernel};
6use super::primitives::{Orientation, Point2};
7
8#[derive(Debug, Clone, Copy, PartialEq, Eq)]
9pub enum HullError {
10    TooManyPoints,
11    ScratchTooSmall { required: usize },
12    OutputTooSmall { required: usize },
13    NonFiniteCoordinate { index: usize },
14}
15
16#[inline]
17fn point_cmp(a: (f64, f64), b: (f64, f64)) -> Ordering {
18    a.0.total_cmp(&b.0).then_with(|| a.1.total_cmp(&b.1))
19}
20
21fn hull_indices_by<K, F>(
22    kernel: &K,
23    count: usize,
24    xy: F,
25    scratch: &mut [u32],
26    out: &mut [u32],
27) -> Result<usize, HullError>
28where
29    K: GeometryKernel,
30    F: Fn(usize) -> (f64, f64) + Copy,
31{
32    if count > u32::MAX as usize {
33        return Err(HullError::TooManyPoints);
34    }
35    let scratch_required = count.saturating_mul(3);
36    if scratch.len() < scratch_required {
37        return Err(HullError::ScratchTooSmall {
38            required: scratch_required,
39        });
40    }
41    if out.len() < count {
42        return Err(HullError::OutputTooSmall { required: count });
43    }
44    if count == 0 {
45        return Ok(0);
46    }
47
48    // Split: order (n for sorting) + stack (2n for Andrew's monotone chain).
49    // The stack needs 2n because the upper hull phase can temporarily hold
50    // up to 2n-1 entries (n from the lower hull + n-1 from the upper hull
51    // before the cross-product check pops duplicates). With only n stack
52    // slots, a convex polygon input (all points on hull) overflows.
53    let (order, stack) = scratch[..scratch_required].split_at_mut(count);
54    for (i, slot) in order.iter_mut().enumerate() {
55        let p = xy(i);
56        if !p.0.is_finite() || !p.1.is_finite() {
57            return Err(HullError::NonFiniteCoordinate { index: i });
58        }
59        *slot = i as u32;
60    }
61    order.sort_unstable_by(|a, b| point_cmp(xy(*a as usize), xy(*b as usize)));
62
63    // Deduplicate coincident input points in place.
64    let mut unique = 1usize;
65    for read in 1..count {
66        if xy(order[read] as usize) != xy(order[unique - 1] as usize) {
67            order[unique] = order[read];
68            unique += 1;
69        }
70    }
71    if unique == 1 {
72        out[0] = order[0];
73        return Ok(1);
74    }
75
76    let turn = |ia: u32, ib: u32, ic: u32| {
77        let a = xy(ia as usize);
78        let b = xy(ib as usize);
79        let c = xy(ic as usize);
80        kernel.orientation_2(
81            Point2::new(a.0, a.1),
82            Point2::new(b.0, b.1),
83            Point2::new(c.0, c.1),
84        )
85    };
86
87    let mut len = 0usize;
88    for &index in order[..unique].iter() {
89        while len >= 2
90            && turn(stack[len - 2], stack[len - 1], index) != Orientation::CounterClockwise
91        {
92            len -= 1;
93        }
94        stack[len] = index;
95        len += 1;
96    }
97
98    let lower_len = len;
99    for &index in order[..unique - 1].iter().rev() {
100        while len > lower_len
101            && turn(stack[len - 2], stack[len - 1], index) != Orientation::CounterClockwise
102        {
103            len -= 1;
104        }
105        if index == order[0] {
106            break;
107        }
108        stack[len] = index;
109        len += 1;
110    }
111    out[..len].copy_from_slice(&stack[..len]);
112    Ok(len)
113}
114
115/// Compute the CCW convex-hull vertex indices with no heap allocation, using
116/// the default [`FilteredF64Kernel`].
117///
118/// `scratch` requires `3 * points.len()` entries and `out` requires
119/// `points.len()` entries. Collinear interior points and duplicates are omitted.
120pub fn convex_hull_indices_2(
121    points: &[Point2],
122    scratch: &mut [u32],
123    out: &mut [u32],
124) -> Result<usize, HullError> {
125    convex_hull_indices_2_with_kernel(&FilteredF64Kernel::default(), points, scratch, out)
126}
127
128/// Kernel-generic variant of [`convex_hull_indices_2`] — the algorithm runs
129/// unchanged over any [`GeometryKernel`] (filtered `f64` today, exact
130/// arithmetic in P1.7). This is the seam where the kernel is swapped.
131pub fn convex_hull_indices_2_with_kernel<K: GeometryKernel>(
132    kernel: &K,
133    points: &[Point2],
134    scratch: &mut [u32],
135    out: &mut [u32],
136) -> Result<usize, HullError> {
137    hull_indices_by(
138        kernel,
139        points.len(),
140        |i| (points[i].x, points[i].y),
141        scratch,
142        out,
143    )
144}
145
146/// Compute CCW convex-hull points into a caller-owned output slice, using
147/// the default [`FilteredF64Kernel`].
148pub fn convex_hull_2(
149    points: &[Point2],
150    scratch: &mut [u32],
151    out: &mut [Point2],
152) -> Result<usize, HullError> {
153    convex_hull_2_with_kernel(&FilteredF64Kernel::default(), points, scratch, out)
154}
155
156/// Kernel-generic variant of [`convex_hull_2`].
157pub fn convex_hull_2_with_kernel<K: GeometryKernel>(
158    kernel: &K,
159    points: &[Point2],
160    scratch: &mut [u32],
161    out: &mut [Point2],
162) -> Result<usize, HullError> {
163    if out.len() < points.len() {
164        return Err(HullError::OutputTooSmall {
165            required: points.len(),
166        });
167    }
168    let required = points.len().saturating_mul(3);
169    if scratch.len() < required {
170        return Err(HullError::ScratchTooSmall { required });
171    }
172    // Split: order (n for sorting) + hull (2n for Andrew's monotone chain
173    // stack — see hull_indices_by for the 2n requirement).
174    let (order, hull) = scratch[..required].split_at_mut(points.len());
175    let count = hull_indices_by_local(kernel, points, order, hull)?;
176    for i in 0..count {
177        out[i] = points[hull[i] as usize];
178    }
179    Ok(count)
180}
181
182fn hull_indices_by_local<K: GeometryKernel>(
183    kernel: &K,
184    points: &[Point2],
185    order: &mut [u32],
186    stack: &mut [u32],
187) -> Result<usize, HullError> {
188    let count = points.len();
189    if count == 0 {
190        return Ok(0);
191    }
192    for (i, slot) in order.iter_mut().enumerate() {
193        let p = points[i];
194        if !p.x.is_finite() || !p.y.is_finite() {
195            return Err(HullError::NonFiniteCoordinate { index: i });
196        }
197        *slot = i as u32;
198    }
199    order.sort_unstable_by(|a, b| {
200        point_cmp(
201            (points[*a as usize].x, points[*a as usize].y),
202            (points[*b as usize].x, points[*b as usize].y),
203        )
204    });
205    let mut unique = 1usize;
206    for read in 1..count {
207        if points[order[read] as usize] != points[order[unique - 1] as usize] {
208            order[unique] = order[read];
209            unique += 1;
210        }
211    }
212    if unique == 1 {
213        stack[0] = order[0];
214        return Ok(1);
215    }
216    let mut len = 0usize;
217    for &index in order[..unique].iter() {
218        while len >= 2
219            && kernel.orientation_2(
220                points[stack[len - 2] as usize],
221                points[stack[len - 1] as usize],
222                points[index as usize],
223            ) != Orientation::CounterClockwise
224        {
225            len -= 1;
226        }
227        stack[len] = index;
228        len += 1;
229    }
230    let lower_len = len;
231    for &index in order[..unique - 1].iter().rev() {
232        while len > lower_len
233            && kernel.orientation_2(
234                points[stack[len - 2] as usize],
235                points[stack[len - 1] as usize],
236                points[index as usize],
237            ) != Orientation::CounterClockwise
238        {
239            len -= 1;
240        }
241        if index == order[0] {
242            break;
243        }
244        stack[len] = index;
245        len += 1;
246    }
247    Ok(len)
248}
249
250/// Convex hull of the spatial `(x,y)` projection of 10D manifold nodes, using
251/// the default [`FilteredF64Kernel`].
252///
253/// The returned indices continue to address the original tensors, preserving
254/// their q/v/w/t/spectral coordinates for graph and reasoning consumers.
255pub fn convex_hull_tensor_xy(
256    points: &[Tensor10D],
257    scratch: &mut [u32],
258    out: &mut [u32],
259) -> Result<usize, HullError> {
260    convex_hull_tensor_xy_with_kernel(&FilteredF64Kernel::default(), points, scratch, out)
261}
262
263/// Kernel-generic variant of [`convex_hull_tensor_xy`].
264pub fn convex_hull_tensor_xy_with_kernel<K: GeometryKernel>(
265    kernel: &K,
266    points: &[Tensor10D],
267    scratch: &mut [u32],
268    out: &mut [u32],
269) -> Result<usize, HullError> {
270    hull_indices_by(
271        kernel,
272        points.len(),
273        |i| (points[i].x as f64, points[i].y as f64),
274        scratch,
275        out,
276    )
277}
278
279/// Check that a polygon is CCW and strongly convex, using the default
280/// [`FilteredF64Kernel`].
281pub fn is_ccw_strongly_convex_2(points: &[Point2]) -> bool {
282    is_ccw_strongly_convex_2_with_kernel(&FilteredF64Kernel::default(), points)
283}
284
285/// Kernel-generic variant of [`is_ccw_strongly_convex_2`].
286pub fn is_ccw_strongly_convex_2_with_kernel<K: GeometryKernel>(
287    kernel: &K,
288    points: &[Point2],
289) -> bool {
290    if points.len() < 3 {
291        return false;
292    }
293    for i in 0..points.len() {
294        if kernel.orientation_2(
295            points[i],
296            points[(i + 1) % points.len()],
297            points[(i + 2) % points.len()],
298        ) != Orientation::CounterClockwise
299        {
300            return false;
301        }
302    }
303    true
304}
305
306#[cfg(test)]
307mod tests {
308    use super::*;
309
310    #[test]
311    fn square_hull_omits_duplicate_and_interior_points() {
312        let points = [
313            Point2::new(0.0, 0.0),
314            Point2::new(1.0, 0.0),
315            Point2::new(0.5, 0.5),
316            Point2::new(1.0, 1.0),
317            Point2::new(0.0, 1.0),
318            Point2::new(0.0, 0.0),
319        ];
320        let mut scratch = [0u32; 18];
321        let mut out = [0u32; 6];
322        let n = convex_hull_indices_2(&points, &mut scratch, &mut out).unwrap();
323        assert_eq!(n, 4);
324        let hull: Vec<Point2> = out[..n]
325            .iter()
326            .map(|&index| points[index as usize])
327            .collect();
328        assert_eq!(
329            hull,
330            vec![
331                Point2::new(0.0, 0.0),
332                Point2::new(1.0, 0.0),
333                Point2::new(1.0, 1.0),
334                Point2::new(0.0, 1.0),
335            ]
336        );
337        assert!(is_ccw_strongly_convex_2(&hull));
338    }
339
340    #[test]
341    fn collinear_hull_contains_only_extremes() {
342        let points = [
343            Point2::new(1.0, 0.0),
344            Point2::new(0.0, 0.0),
345            Point2::new(2.0, 0.0),
346        ];
347        let mut scratch = [0u32; 9];
348        let mut out = [0u32; 3];
349        let n = convex_hull_indices_2(&points, &mut scratch, &mut out).unwrap();
350        assert_eq!(n, 2);
351        assert_eq!(points[out[0] as usize], Point2::new(0.0, 0.0));
352        assert_eq!(points[out[1] as usize], Point2::new(2.0, 0.0));
353    }
354
355    #[test]
356    fn ten_dimensional_hull_preserves_source_indices() {
357        let mut points = [Tensor10D::default(); 5];
358        for (point, xy) in
359            points
360                .iter_mut()
361                .zip([[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0], [0.5, 0.5]])
362        {
363            point.x = xy[0];
364            point.y = xy[1];
365        }
366        points[2].q = 7.0;
367        let mut scratch = [0u32; 15];
368        let mut out = [0u32; 5];
369        let n = convex_hull_tensor_xy(&points, &mut scratch, &mut out).unwrap();
370        assert_eq!(&out[..n], &[0, 1, 2, 3]);
371        assert_eq!(points[out[2] as usize].q, 7.0);
372    }
373
374    #[test]
375    fn reports_caller_buffer_requirements() {
376        let points = [Point2::new(0.0, 0.0); 4];
377        let mut scratch = [0u32; 11];
378        let mut out = [0u32; 4];
379        assert_eq!(
380            convex_hull_indices_2(&points, &mut scratch, &mut out),
381            Err(HullError::ScratchTooSmall { required: 12 })
382        );
383    }
384
385    #[test]
386    fn kernel_generic_path_matches_default_path() {
387        // The P1.2 contract: the same algorithm over any GeometryKernel
388        // produces identical output. The default (FilteredF64Kernel) and the
389        // explicit kernel-generic call must agree byte-for-byte on the hull
390        // indices.
391        use super::FilteredF64Kernel;
392
393        let points = [
394            Point2::new(0.0, 0.0),
395            Point2::new(1.0, 0.0),
396            Point2::new(0.5, 0.5),
397            Point2::new(1.0, 1.0),
398            Point2::new(0.0, 1.0),
399            Point2::new(0.25, 0.75),
400        ];
401        let mut scratch_a = [0u32; 18];
402        let mut out_a = [0u32; 6];
403        let n_a = convex_hull_indices_2(&points, &mut scratch_a, &mut out_a).unwrap();
404
405        let mut scratch_b = [0u32; 18];
406        let mut out_b = [0u32; 6];
407        let n_b = convex_hull_indices_2_with_kernel(
408            &FilteredF64Kernel::default(),
409            &points,
410            &mut scratch_b,
411            &mut out_b,
412        )
413        .unwrap();
414
415        assert_eq!(n_a, n_b);
416        assert_eq!(&out_a[..n_a], &out_b[..n_b]);
417    }
418
419    #[test]
420    fn strongly_convex_check_matches_through_kernel() {
421        use super::FilteredF64Kernel;
422
423        let hull = [
424            Point2::new(0.0, 0.0),
425            Point2::new(1.0, 0.0),
426            Point2::new(1.0, 1.0),
427            Point2::new(0.0, 1.0),
428        ];
429        assert!(is_ccw_strongly_convex_2(&hull));
430        assert!(is_ccw_strongly_convex_2_with_kernel(
431            &FilteredF64Kernel::default(),
432            &hull
433        ));
434
435        let non_convex = [
436            Point2::new(0.0, 0.0),
437            Point2::new(1.0, 0.0),
438            Point2::new(0.5, 0.5),
439            Point2::new(1.0, 1.0),
440        ];
441        assert!(!is_ccw_strongly_convex_2(&non_convex));
442        assert!(!is_ccw_strongly_convex_2_with_kernel(
443            &FilteredF64Kernel::default(),
444            &non_convex
445        ));
446    }
447
448    /// Regression test: all points on the convex hull must not overflow the
449    /// stack buffer. Before the fix, the upper hull phase of Andrew's
450    /// monotone chain could temporarily hold 2n entries, overflowing the
451    /// n-element stack. This hexagon (all 6 vertices on the hull) reproduces
452    /// the original panic.
453    #[test]
454    fn convex_polygon_all_points_on_hull_no_overflow() {
455        let points = [
456            Point2::new(0.0, 0.0),
457            Point2::new(2.0, 0.0),
458            Point2::new(3.0, 2.0),
459            Point2::new(2.0, 4.0),
460            Point2::new(0.0, 4.0),
461            Point2::new(-1.0, 2.0),
462        ];
463        let mut scratch = [0u32; 18]; // 3 * 6
464        let mut out = [0u32; 6];
465        let n = convex_hull_indices_2(&points, &mut scratch, &mut out).unwrap();
466        assert_eq!(n, 6, "all 6 vertices should be on the hull");
467        let hull: Vec<Point2> = out[..n].iter().map(|&i| points[i as usize]).collect();
468        assert!(is_ccw_strongly_convex_2(&hull));
469    }
470}
qualia_core_db/specialized_libs/computational_geometry/hull.rs

qualia_core_db/specialized_libs/computational_geometry/
hull.rs
