axiolid_nurbs/
certified_curve_projection.rs

1//! Exhaustive closest-point bounds for piecewise rational Bézier curves.
2
3use crate::certified_bezier::{
4    distance_to_point_interval_upper, next_up, representative_distance, Cell, HomogeneousPoint,
5};
6use crate::certified_projection::{
7    CertifiedProjectionOptions, CurveProjectionCertificate2, CurveProjectionCertificate3,
8    ParameterInterval,
9};
10use crate::certified_refinement::{piecewise_bezier_cells, RefinementBudget};
11use axiolid_contracts::{GeomError, GeomResult};
12use axiolid_core::{Point2, Point3, Scalar};
13use axiolid_curve::{BSplineCurve2, BSplineCurve3};
14use axiolid_evaluate::curve::{bspline_jet2, bspline_jet3};
15use core::cmp::Ordering;
16use std::collections::BinaryHeap;
17
18/// Globally bound the closest point on a clamped planar B-spline curve.
19///
20/// Positive rational weights make every refined segment lie in the convex hull
21/// of its Euclidean control points. Interval-aware homogeneous knot insertion
22/// and outward-rounded subdivision tighten those hulls until the global
23/// distance gap meets the requested tolerance.
24pub fn project_curve2_certified(
25    curve: &BSplineCurve2,
26    target: Point2,
27    options: CertifiedProjectionOptions,
28) -> GeomResult<CurveProjectionCertificate2> {
29    let mut refinement_budget =
30        RefinementBudget::new(options.max_nodes(), "certified projection nodes");
31    let roots = piecewise_bezier_cells(
32        curve,
33        |point| [point.x, point.y, 0.0],
34        &mut refinement_budget,
35    )?;
36    let result = project(roots, [target.x, target.y, 0.0], 2, options, |parameter| {
37        let point = bspline_jet2(curve, parameter)?.point;
38        Ok([point.x, point.y, 0.0])
39    })?;
40    Ok(CurveProjectionCertificate2 {
41        parameter: result.parameter,
42        point: Point2::new(result.point[0], result.point[1]),
43        distance: result.distance,
44        distance_lower_bound: result.lower,
45        distance_upper_bound: result.upper,
46        possible_minimizer_intervals: result.intervals,
47        visited_nodes: result.nodes,
48    })
49}
50
51/// Globally bound the closest point on a clamped spatial B-spline curve.
52///
53/// See [`project_curve2_certified`] for the certification and current input
54/// scope. A budget failure never returns a partial certificate.
55pub fn project_curve3_certified(
56    curve: &BSplineCurve3,
57    target: Point3,
58    options: CertifiedProjectionOptions,
59) -> GeomResult<CurveProjectionCertificate3> {
60    let mut refinement_budget =
61        RefinementBudget::new(options.max_nodes(), "certified projection nodes");
62    let roots = piecewise_bezier_cells(
63        curve,
64        |point| [point.x, point.y, point.z],
65        &mut refinement_budget,
66    )?;
67    let result = project(roots, target.to_array(), 3, options, |parameter| {
68        let point = bspline_jet3(curve, parameter)?.point;
69        Ok(point.to_array())
70    })?;
71    Ok(CurveProjectionCertificate3 {
72        parameter: result.parameter,
73        point: Point3::from_array(result.point),
74        distance: result.distance,
75        distance_lower_bound: result.lower,
76        distance_upper_bound: result.upper,
77        possible_minimizer_intervals: result.intervals,
78        visited_nodes: result.nodes,
79    })
80}
81
82#[derive(Debug)]
83struct QueueCell {
84    lower: Scalar,
85    serial: u64,
86    cell: Cell,
87}
88
89impl PartialEq for QueueCell {
90    fn eq(&self, other: &Self) -> bool {
91        self.serial == other.serial
92    }
93}
94impl Eq for QueueCell {}
95impl PartialOrd for QueueCell {
96    fn partial_cmp(&self, other: &Self) -> Option<Ordering> {
97        Some(self.cmp(other))
98    }
99}
100impl Ord for QueueCell {
101    fn cmp(&self, other: &Self) -> Ordering {
102        other
103            .lower
104            .total_cmp(&self.lower)
105            .then_with(|| other.serial.cmp(&self.serial))
106    }
107}
108
109#[derive(Debug, Clone)]
110struct Candidate {
111    parameter: Scalar,
112    point: [Scalar; 3],
113    distance: Scalar,
114    upper: Scalar,
115}
116
117struct ProjectionCore {
118    parameter: Scalar,
119    point: [Scalar; 3],
120    distance: Scalar,
121    lower: Scalar,
122    upper: Scalar,
123    intervals: Vec<ParameterInterval>,
124    nodes: u32,
125}
126
127fn project(
128    roots: Vec<Cell>,
129    target: [Scalar; 3],
130    dimensions: usize,
131    options: CertifiedProjectionOptions,
132    evaluate: impl Fn(Scalar) -> GeomResult<[Scalar; 3]>,
133) -> GeomResult<ProjectionCore> {
134    if target[..dimensions].iter().any(|value| !value.is_finite()) {
135        return Err(GeomError::InvalidInput(
136            "projection target must be finite".to_owned(),
137        ));
138    }
139    let mut nodes = u32::try_from(roots.len()).map_err(|_| GeomError::BudgetExceeded {
140        resource: "certified projection nodes",
141    })?;
142    if nodes > options.max_nodes() {
143        return Err(GeomError::BudgetExceeded {
144            resource: "certified projection nodes",
145        });
146    }
147
148    let mut best: Option<Candidate> = None;
149    let mut heap = BinaryHeap::new();
150    let mut serial = 0_u64;
151    for cell in roots {
152        consider_parameter(
153            &mut best,
154            cell.start,
155            &cell.controls[0],
156            target,
157            dimensions,
158            &evaluate,
159        )?;
160        consider_parameter(
161            &mut best,
162            cell.end,
163            &cell.controls[cell.controls.len() - 1],
164            target,
165            dimensions,
166            &evaluate,
167        )?;
168        let midpoint = cell.start * 0.5 + cell.end * 0.5;
169        consider_parameter(
170            &mut best,
171            midpoint,
172            &cell.midpoint_point()?,
173            target,
174            dimensions,
175            &evaluate,
176        )?;
177        let lower = cell.lower_bound(target, dimensions)?;
178        heap.push(QueueCell {
179            lower,
180            serial,
181            cell,
182        });
183        serial += 1;
184    }
185    let mut best = best.ok_or_else(|| {
186        GeomError::InvalidInput("certified projection has no curve segments".to_owned())
187    })?;
188
189    loop {
190        while heap.peek().is_some_and(|entry| entry.lower > best.upper) {
191            heap.pop();
192        }
193        let global_lower = heap
194            .peek()
195            .map_or(best.upper, |entry| entry.lower.min(best.upper));
196        if best.upper - global_lower <= options.tolerance().linear() {
197            let mut intervals: Vec<_> = heap
198                .iter()
199                .filter(|entry| entry.lower <= best.upper)
200                .map(|entry| ParameterInterval {
201                    start: entry.cell.start,
202                    end: entry.cell.end,
203                })
204                .collect();
205            intervals.push(ParameterInterval {
206                start: best.parameter,
207                end: best.parameter,
208            });
209            intervals.sort_by(|left, right| left.start.total_cmp(&right.start));
210            return Ok(ProjectionCore {
211                parameter: best.parameter,
212                point: best.point,
213                distance: best.distance,
214                lower: global_lower,
215                upper: best.upper,
216                intervals,
217                nodes,
218            });
219        }
220
221        let current = heap.pop().ok_or_else(|| {
222            GeomError::Degenerate("certified projection queue became empty".to_owned())
223        })?;
224        if current.cell.depth >= options.max_depth() {
225            return Err(GeomError::BudgetExceeded {
226                resource: "certified projection depth",
227            });
228        }
229        let next_nodes = nodes.checked_add(2).ok_or(GeomError::BudgetExceeded {
230            resource: "certified projection nodes",
231        })?;
232        if next_nodes > options.max_nodes() {
233            return Err(GeomError::BudgetExceeded {
234                resource: "certified projection nodes",
235            });
236        }
237        nodes = next_nodes;
238        let (left, right) = current.cell.split()?;
239        for child in [left, right] {
240            let midpoint = child.start * 0.5 + child.end * 0.5;
241            consider_parameter(
242                &mut best,
243                midpoint,
244                &child.midpoint_point()?,
245                target,
246                dimensions,
247                &evaluate,
248            )?;
249            let lower = child.lower_bound(target, dimensions)?;
250            if lower <= best.upper {
251                heap.push(QueueCell {
252                    lower,
253                    serial,
254                    cell: child,
255                });
256                serial = serial.checked_add(1).ok_or(GeomError::BudgetExceeded {
257                    resource: "certified projection serials",
258                })?;
259            }
260        }
261    }
262}
263
264fn consider_parameter(
265    best: &mut impl CandidateSlot,
266    parameter: Scalar,
267    enclosure: &HomogeneousPoint,
268    target: [Scalar; 3],
269    dimensions: usize,
270    evaluate: &impl Fn(Scalar) -> GeomResult<[Scalar; 3]>,
271) -> GeomResult<()> {
272    let point = evaluate(parameter)?;
273    if point[..dimensions].iter().any(|value| !value.is_finite()) {
274        return Err(GeomError::Degenerate(
275            "scalar projection evaluation is non-finite".to_owned(),
276        ));
277    }
278    let distance = representative_distance(point, target, dimensions)?;
279    let upper = distance_to_point_interval_upper(target, enclosure.euclidean()?, dimensions)?
280        .max(next_up(distance));
281    let candidate = Candidate {
282        parameter,
283        point,
284        distance,
285        upper,
286    };
287    best.consider(candidate);
288    Ok(())
289}
290
291trait CandidateSlot {
292    fn consider(&mut self, candidate: Candidate);
293}
294
295impl CandidateSlot for Option<Candidate> {
296    fn consider(&mut self, candidate: Candidate) {
297        if self.as_ref().is_none_or(|current| {
298            candidate.upper < current.upper
299                || (candidate.upper == current.upper && candidate.parameter < current.parameter)
300        }) {
301            *self = Some(candidate);
302        }
303    }
304}
305
306impl CandidateSlot for Candidate {
307    fn consider(&mut self, candidate: Candidate) {
308        if candidate.upper < self.upper
309            || (candidate.upper == self.upper && candidate.parameter < self.parameter)
310        {
311            *self = candidate;
312        }
313    }
314}
axiolid_nurbs/certified_curve_projection.rs

axiolid_nurbs/
certified_curve_projection.rs
