axiolid_nurbs/
curve_projection.rs

1//! Bounded inverse queries for NURBS curves.
2
3use crate::axis::active_spans;
4use crate::projection::{CurveProjection2, CurveProjection3, ProjectionOptions, ProjectionStatus};
5use axiolid_contracts::{GeomError, GeomResult};
6use axiolid_core::{Point2, Point3, Scalar};
7use axiolid_curve::{BSplineCurve2, BSplineCurve3};
8use axiolid_evaluate::curve::{bspline_jet2, bspline_jet3};
9
10/// Find the best planar-curve projection candidate within explicit budgets.
11///
12/// Every active knot span and both domain endpoints are seeded. The result is
13/// not a certified global minimum; inspect its status and distance.
14pub fn project_curve2(
15    curve: &BSplineCurve2,
16    target: Point2,
17    options: ProjectionOptions,
18) -> GeomResult<CurveProjection2> {
19    let spans = active_spans(
20        &curve.knots,
21        &curve.multiplicities,
22        curve.degree,
23        curve.control_points.len(),
24    )?;
25    let (t, p, d, it, status, on_boundary) = project([target.x, target.y], &spans, options, |t| {
26        let j = bspline_jet2(curve, t)?;
27        Ok((
28            [j.point.x, j.point.y],
29            [j.first.x, j.first.y],
30            [j.second.x, j.second.y],
31        ))
32    })?;
33    Ok(CurveProjection2 {
34        parameter: t,
35        point: Point2::new(p[0], p[1]),
36        distance: d,
37        iterations: it,
38        on_boundary,
39        status,
40    })
41}
42
43/// Find the best spatial-curve projection candidate within explicit budgets.
44///
45/// Every active knot span and both domain endpoints are seeded. The result is
46/// not a certified global minimum; inspect its status and distance.
47pub fn project_curve3(
48    curve: &BSplineCurve3,
49    target: Point3,
50    options: ProjectionOptions,
51) -> GeomResult<CurveProjection3> {
52    let spans = active_spans(
53        &curve.knots,
54        &curve.multiplicities,
55        curve.degree,
56        curve.control_points.len(),
57    )?;
58    let (t, p, d, it, status, on_boundary) =
59        project([target.x, target.y, target.z], &spans, options, |t| {
60            let j = bspline_jet3(curve, t)?;
61            Ok((
62                [j.point.x, j.point.y, j.point.z],
63                [j.first.x, j.first.y, j.first.z],
64                [j.second.x, j.second.y, j.second.z],
65            ))
66        })?;
67    Ok(CurveProjection3 {
68        parameter: t,
69        point: Point3::new(p[0], p[1], p[2]),
70        distance: d,
71        iterations: it,
72        on_boundary,
73        status,
74    })
75}
76
77type Candidate<const N: usize> = (Scalar, [Scalar; N], Scalar, u16, ProjectionStatus, bool);
78
79fn project<const N: usize>(
80    target: [Scalar; N],
81    spans: &[(Scalar, Scalar)],
82    options: ProjectionOptions,
83    jet: impl Fn(Scalar) -> GeomResult<([Scalar; N], [Scalar; N], [Scalar; N])>,
84) -> GeomResult<Candidate<N>> {
85    if target.iter().any(|x| !x.is_finite()) {
86        return Err(GeomError::InvalidInput(
87            "projection target must be finite".to_owned(),
88        ));
89    }
90    let lo = spans[0].0;
91    let hi = spans[spans.len() - 1].1;
92    let mut best: Option<Candidate<N>> = None;
93    let mut starts = 0_u32;
94    for &(a, b) in spans {
95        for sample in 0..=options.samples_per_span() {
96            starts = starts.checked_add(1).ok_or(GeomError::BudgetExceeded {
97                resource: "projection starts",
98            })?;
99            if starts > options.max_starts() {
100                return Err(GeomError::BudgetExceeded {
101                    resource: "projection starts",
102                });
103            }
104            let start =
105                a + (b - a) * Scalar::from(sample) / Scalar::from(options.samples_per_span());
106            let candidate = refine(target, start, lo, hi, options, &jet)?;
107            if best.as_ref().is_none_or(|current| candidate.2 < current.2) {
108                best = Some(candidate);
109            }
110        }
111    }
112    best.ok_or_else(|| GeomError::Degenerate("projection produced no candidate".to_owned()))
113}
114
115fn refine<const N: usize>(
116    target: [Scalar; N],
117    start: Scalar,
118    lo: Scalar,
119    hi: Scalar,
120    options: ProjectionOptions,
121    jet: &impl Fn(Scalar) -> GeomResult<([Scalar; N], [Scalar; N], [Scalar; N])>,
122) -> GeomResult<Candidate<N>> {
123    let mut t = start;
124    let mut iterations = 0;
125    let mut status = ProjectionStatus::BudgetExhausted;
126    for iteration in 0..options.max_iterations() {
127        iterations = iteration + 1;
128        let (p, d1, d2) = jet(t)?;
129        let r = sub(p, target);
130        let speed = dot(d1, d1).sqrt();
131        let gradient = dot(r, d1);
132        if gradient.abs() <= options.tolerance().linear() * speed.max(1.0) {
133            status = ProjectionStatus::Converged;
134            break;
135        }
136        let hessian = dot(d1, d1) + dot(r, d2);
137        if !hessian.is_finite() || hessian.abs() <= Scalar::EPSILON * dot(d1, d1).max(1.0) {
138            break;
139        }
140        let next = (t - gradient / hessian).clamp(lo, hi);
141        if (next - t).abs() * speed <= options.tolerance().linear() {
142            t = next;
143            status = ProjectionStatus::Converged;
144            break;
145        }
146        t = next;
147    }
148    let (point, _, _) = jet(t)?;
149    let distance = dot(sub(point, target), sub(point, target)).sqrt();
150    if !distance.is_finite() {
151        return Err(GeomError::Degenerate(
152            "projection distance is non-finite".to_owned(),
153        ));
154    }
155    Ok((t, point, distance, iterations, status, t == lo || t == hi))
156}
157
158fn sub<const N: usize>(a: [Scalar; N], b: [Scalar; N]) -> [Scalar; N] {
159    core::array::from_fn(|i| a[i] - b[i])
160}
161fn dot<const N: usize>(a: [Scalar; N], b: [Scalar; N]) -> Scalar {
162    (0..N).map(|i| a[i] * b[i]).sum()
163}
axiolid_nurbs/curve_projection.rs

axiolid_nurbs/
curve_projection.rs
