brepkit_geometry/sampling/
deflection.rs

1//! Adaptive deflection-based curve sampling via recursive midpoint subdivision.
2
3use brepkit_math::traits::ParametricCurve;
4use brepkit_math::vec::{Point3, Vec3};
5
6/// Maximum recursion depth to guard against degenerate curves.
7const MAX_DEPTH: u32 = 20;
8
9/// Compute the perpendicular distance from point `p` to the chord `a → b`.
10///
11/// Uses the formula `|(p - a) × (b - a)| / |b - a|`.
12/// Returns `0.0` when `a` and `b` coincide (degenerate chord).
13fn chord_deviation(p: Point3, a: Point3, b: Point3) -> f64 {
14    let ab: Vec3 = b - a;
15    let ab_len = ab.length();
16    if ab_len < f64::EPSILON {
17        return 0.0;
18    }
19    let ap: Vec3 = p - a;
20    ap.cross(ab).length() / ab_len
21}
22
23/// Recursively subdivide the interval `[t_a, t_b]` until the chord deviation
24/// at the midpoint is below `max_deflection`.
25///
26/// `p_a` and `p_b` are the already-evaluated curve points at `t_a` and `t_b`.
27/// New points are appended to `out` (excluding `p_a`; `p_b` is added by the
28/// outermost caller after all recursion completes).
29#[allow(clippy::too_many_arguments)]
30fn subdivide<C: ParametricCurve>(
31    curve: &C,
32    t_a: f64,
33    p_a: Point3,
34    t_b: f64,
35    p_b: Point3,
36    max_deflection: f64,
37    depth: u32,
38    out: &mut Vec<(f64, Point3)>,
39) {
40    if depth >= MAX_DEPTH {
41        return;
42    }
43    let t_m = 0.5 * (t_a + t_b);
44    let p_m = curve.evaluate(t_m);
45
46    if chord_deviation(p_m, p_a, p_b) <= max_deflection {
47        // Chord is within tolerance — no need to subdivide further.
48        return;
49    }
50
51    subdivide(curve, t_a, p_a, t_m, p_m, max_deflection, depth + 1, out);
52    out.push((t_m, p_m));
53    subdivide(curve, t_m, p_m, t_b, p_b, max_deflection, depth + 1, out);
54}
55
56/// Adaptively sample a curve so that every chord's midpoint deviation is
57/// below `max_deflection`.
58///
59/// Returns `(t, Point3)` pairs sorted by increasing `t`, always including
60/// the endpoints `t_start` and `t_end`.
61///
62/// If `max_deflection` is non-positive, the function returns only the two
63/// endpoints (no subdivision).
64#[must_use]
65pub fn sample_deflection<C: ParametricCurve>(
66    curve: &C,
67    t_start: f64,
68    t_end: f64,
69    max_deflection: f64,
70) -> Vec<(f64, Point3)> {
71    let p_start = curve.evaluate(t_start);
72    let p_end = curve.evaluate(t_end);
73
74    let mut out = Vec::new();
75    out.push((t_start, p_start));
76
77    if max_deflection > 0.0 {
78        subdivide(
79            curve,
80            t_start,
81            p_start,
82            t_end,
83            p_end,
84            max_deflection,
85            0,
86            &mut out,
87        );
88    }
89
90    out.push((t_end, p_end));
91    out
92}
93
94#[cfg(test)]
95mod tests {
96    #![allow(clippy::unwrap_used, clippy::expect_used)]
97
98    use super::*;
99    use brepkit_math::curves::Circle3D;
100    use brepkit_math::vec::{Point3, Vec3};
101    use std::f64::consts::TAU;
102
103    fn circle(radius: f64) -> Circle3D {
104        Circle3D::new(Point3::new(0.0, 0.0, 0.0), Vec3::new(0.0, 0.0, 1.0), radius).unwrap()
105    }
106
107    #[test]
108    fn deflection_circle_r10_many_points() {
109        let c = circle(10.0);
110        let max_dev = 0.01;
111        let pairs = sample_deflection(&c, 0.0, TAU, max_dev);
112
113        // Should produce many points to satisfy the tight tolerance on a large circle.
114        assert!(
115            pairs.len() > 10,
116            "expected many points, got {}",
117            pairs.len()
118        );
119    }
120
121    #[test]
122    fn deflection_every_midpoint_within_tolerance() {
123        let c = circle(10.0);
124        let max_dev = 0.01;
125        let pairs = sample_deflection(&c, 0.0, TAU, max_dev);
126
127        // For every consecutive pair, verify the midpoint chord deviation is ≤ max_dev.
128        for window in pairs.windows(2) {
129            let (t_a, p_a) = window[0];
130            let (t_b, p_b) = window[1];
131            let t_m = 0.5 * (t_a + t_b);
132            let p_m = c.evaluate(t_m);
133            let dev = chord_deviation(p_m, p_a, p_b);
134            assert!(
135                dev <= max_dev + 1e-12,
136                "chord deviation {dev} exceeds max {max_dev} between t={t_a} and t={t_b}"
137            );
138        }
139    }
140
141    #[test]
142    fn deflection_endpoints_always_included() {
143        let c = circle(1.0);
144        let pairs = sample_deflection(&c, 0.0, TAU, 0.1);
145        assert!(!pairs.is_empty());
146        assert!((pairs.first().unwrap().0 - 0.0).abs() < 1e-12);
147        assert!((pairs.last().unwrap().0 - TAU).abs() < 1e-12);
148    }
149
150    #[test]
151    fn non_positive_deflection_returns_two_endpoints() {
152        let c = circle(1.0);
153        let pairs = sample_deflection(&c, 0.0, TAU, 0.0);
154        assert_eq!(pairs.len(), 2);
155        let pairs_neg = sample_deflection(&c, 0.0, TAU, -1.0);
156        assert_eq!(pairs_neg.len(), 2);
157    }
158
159    #[test]
160    fn points_sorted_by_parameter() {
161        let c = circle(5.0);
162        let pairs = sample_deflection(&c, 0.0, TAU, 0.05);
163        for w in pairs.windows(2) {
164            assert!(
165                w[0].0 < w[1].0,
166                "parameters not sorted: {} >= {}",
167                w[0].0,
168                w[1].0
169            );
170        }
171    }
172}
brepkit_geometry/sampling/deflection.rs

brepkit_geometry/sampling/
deflection.rs
