axiolid_nurbs/
surface_projection.rs

1//! Bounded inverse queries for tensor-product NURBS surfaces.
2
3use crate::axis::active_spans;
4use crate::projection::{ProjectionOptions, ProjectionStatus, SurfaceProjection};
5use axiolid_contracts::{GeomError, GeomResult};
6use axiolid_core::{Point3, Scalar};
7use axiolid_evaluate::surface::bspline_jet;
8use axiolid_surface::BSplineSurface;
9
10/// Find the best tensor-product surface projection candidate within budgets.
11///
12/// The Cartesian product of active U/V knot-span seeds is refined with the
13/// exact squared-distance Hessian. The result is not a certified global minimum.
14pub fn project_surface(
15    surface: &BSplineSurface,
16    target: Point3,
17    options: ProjectionOptions,
18) -> GeomResult<SurfaceProjection> {
19    if !target.is_finite() {
20        return Err(GeomError::InvalidInput(
21            "projection target must be finite".to_owned(),
22        ));
23    }
24    bspline_jet(surface, 0.0, 0.0)?;
25    let us = active_spans(
26        &surface.u_knots,
27        &surface.u_multiplicities,
28        surface.u_degree,
29        surface.control_points.len(),
30    )?;
31    let vc = surface.control_points.first().map_or(0, Vec::len);
32    let vs = active_spans(
33        &surface.v_knots,
34        &surface.v_multiplicities,
35        surface.v_degree,
36        vc,
37    )?;
38    let bounds = (us[0].0, us[us.len() - 1].1, vs[0].0, vs[vs.len() - 1].1);
39    let mut best: Option<SurfaceProjection> = None;
40    let mut starts = 0_u32;
41    for &(ua, ub) in &us {
42        for &(va, vb) in &vs {
43            for iu in 0..=options.samples_per_span() {
44                for iv in 0..=options.samples_per_span() {
45                    starts = starts.checked_add(1).ok_or(GeomError::BudgetExceeded {
46                        resource: "projection starts",
47                    })?;
48                    if starts > options.max_starts() {
49                        return Err(GeomError::BudgetExceeded {
50                            resource: "projection starts",
51                        });
52                    }
53                    let u = ua
54                        + (ub - ua) * Scalar::from(iu) / Scalar::from(options.samples_per_span());
55                    let v = va
56                        + (vb - va) * Scalar::from(iv) / Scalar::from(options.samples_per_span());
57                    let candidate = refine(surface, target, u, v, bounds, options)?;
58                    if best
59                        .as_ref()
60                        .is_none_or(|current| candidate.distance < current.distance)
61                    {
62                        best = Some(candidate);
63                    }
64                }
65            }
66        }
67    }
68    best.ok_or_else(|| GeomError::Degenerate("projection produced no candidate".to_owned()))
69}
70
71fn refine(
72    surface: &BSplineSurface,
73    target: Point3,
74    mut u: Scalar,
75    mut v: Scalar,
76    bounds: (Scalar, Scalar, Scalar, Scalar),
77    options: ProjectionOptions,
78) -> GeomResult<SurfaceProjection> {
79    let (ulo, uhi, vlo, vhi) = bounds;
80    let mut iterations = 0;
81    let mut status = ProjectionStatus::BudgetExhausted;
82    for iteration in 0..options.max_iterations() {
83        iterations = iteration + 1;
84        let j = bspline_jet(surface, u, v)?;
85        let r = j.point - target;
86        let gu = r.dot(j.du);
87        let gv = r.dot(j.dv);
88        let scale = j.du.length().max(j.dv.length()).max(1.0);
89        if gu.hypot(gv) <= options.tolerance().linear() * scale {
90            status = ProjectionStatus::Converged;
91            break;
92        }
93        let huu = j.du.dot(j.du) + r.dot(j.duu);
94        let huv = j.du.dot(j.dv) + r.dot(j.duv);
95        let hvv = j.dv.dot(j.dv) + r.dot(j.dvv);
96        let det = huu * hvv - huv * huv;
97        let hscale = huu.abs().max(huv.abs()).max(hvv.abs()).max(1.0);
98        if !det.is_finite() || det.abs() <= Scalar::EPSILON * hscale * hscale {
99            break;
100        }
101        let du = (-gu * hvv + huv * gv) / det;
102        let dv = (huv * gu - huu * gv) / det;
103        let next_u = (u + du).clamp(ulo, uhi);
104        let next_v = (v + dv).clamp(vlo, vhi);
105        let movement = (j.du * (next_u - u) + j.dv * (next_v - v)).length();
106        u = next_u;
107        v = next_v;
108        if movement <= options.tolerance().linear() {
109            status = ProjectionStatus::Converged;
110            break;
111        }
112    }
113    let point = bspline_jet(surface, u, v)?.point;
114    let distance = point.distance(target);
115    if !distance.is_finite() {
116        return Err(GeomError::Degenerate(
117            "projection distance is non-finite".to_owned(),
118        ));
119    }
120    Ok(SurfaceProjection {
121        u,
122        v,
123        point,
124        distance,
125        iterations,
126        on_boundary: boundary(u, v, bounds),
127        status,
128    })
129}
130
131fn boundary(u: Scalar, v: Scalar, bounds: (Scalar, Scalar, Scalar, Scalar)) -> bool {
132    u == bounds.0 || u == bounds.1 || v == bounds.2 || v == bounds.3
133}
axiolid_nurbs/surface_projection.rs

axiolid_nurbs/
surface_projection.rs
