ogeom_algo/
project_plane.rs

1//! An edge projected orthogonally onto a plane, as an exact curve in the
2//! plane's own coordinates.
3//!
4//! Orthogonal projection onto a plane is affine, so it keeps every curve
5//! whose family is closed under affine maps: a line stays a line (or
6//! collapses to a point), a circle or ellipse becomes an ellipse (a circle
7//! when its plane is parallel, a segment when it stands square), and a
8//! B-spline becomes the B-spline on its projected control points, weights
9//! unchanged. Any other curve is fitted, and the fit's error is returned
10//! with it.
11
12use ogeom_core::{OgeomResult, Tolerances, ogeom_bail};
13use ogeom_geom::{BSpline2d, BSplineCurve, Curve, Curve3d as _, Transformable as _};
14use ogeom_math::{KnotVector, Plane, Point, Point2, Vector, Vector2, Weighted};
15use ogeom_topo::{EdgeRepr, Model, NodeData, Shape};
16
17/// An edge's orthogonal projection onto a plane, in the plane frame's `x`
18/// and `y` coordinates.
19#[derive(Debug, Clone, PartialEq)]
20pub enum ProjectedCurve {
21    /// The edge projects to one point: a line along the plane's normal.
22    Point(Point2),
23    /// A segment.
24    Line {
25        /// One end.
26        start: Point2,
27        /// The other end.
28        end: Point2,
29    },
30    /// A circle or circular arc, `centre + radius (cos t, sin t)`, turning
31    /// counter-clockwise about the plane's normal from `range.0` to
32    /// `range.1`; a full turn where the range spans one.
33    Circle {
34        /// The centre.
35        centre: Point2,
36        /// The radius.
37        radius: f64,
38        /// The arc's angles, `range.0 < range.1`.
39        range: (f64, f64),
40    },
41    /// An ellipse or elliptical arc, `centre + major cos t + minor sin t`
42    /// with `minor` the major axis turned a quarter counter-clockwise and
43    /// scaled by `ratio`, from `range.0` to `range.1`.
44    Ellipse {
45        /// The centre.
46        centre: Point2,
47        /// The major semi-axis, as a vector.
48        major: Vector2,
49        /// The minor semi-axis over the major, in `(0, 1)`.
50        ratio: f64,
51        /// The arc's eccentric angles, `range.0 < range.1`.
52        range: (f64, f64),
53    },
54    /// A B-spline: exact for a B-spline edge, fitted for any other curve,
55    /// with the fit's largest miss in `fit_error`.
56    BSpline {
57        /// The curve.
58        curve: BSpline2d,
59        /// `None` where the projection is exact.
60        fit_error: Option<f64>,
61    },
62}
63
64/// The orthogonal projection of an edge onto a plane, in the plane's own
65/// 2D coordinates (its frame's `x` and `y` axes). The edge's placement is
66/// applied first.
67///
68/// # Errors
69///
70/// [`OgeomError::Construction`](ogeom_core::OgeomError::Construction) if
71/// `edge` is not an edge with a curve of its own, or if a curve with no
72/// closed-form projection cannot be fitted; plus whatever evaluating the
73/// curve refuses.
74pub fn project_edge_onto_plane(
75    model: &Model,
76    edge: &Shape,
77    plane: &Plane,
78    tol: Tolerances,
79) -> OgeomResult<ProjectedCurve> {
80    let Some(NodeData::Edge(data)) = model.node(edge).map(|n| n.data()) else {
81        ogeom_bail!(Construction, "only an edge projects onto a plane");
82    };
83    let Some(EdgeRepr::Curve3d { curve, range, .. }) = data.curve3d() else {
84        ogeom_bail!(
85            Construction,
86            "the edge has no curve in space (a degenerate edge), so there is nothing to project"
87        );
88    };
89    let Some(curve) = model.geometry().curve(*curve) else {
90        ogeom_bail!(Dangling, "the edge's curve is not in this model");
91    };
92    let curve = curve
93        .clone()
94        .transformed(&edge.transform(model.datums())?, tol)?;
95    project_curve(&curve, *range, &Onto::new(plane), tol)
96}
97
98/// The plane, as the map it is.
99struct Onto {
100    origin: Point,
101    x: Vector,
102    y: Vector,
103}
104
105impl Onto {
106    fn new(plane: &Plane) -> Self {
107        let frame = plane.frame();
108        Self {
109            origin: frame.origin(),
110            x: frame.x().vector(),
111            y: frame.y().vector(),
112        }
113    }
114
115    fn point(&self, p: Point) -> Point2 {
116        let d = p - self.origin;
117        Point2::new(d.dot(self.x), d.dot(self.y))
118    }
119
120    fn vector(&self, v: Vector) -> Vector2 {
121        Vector2::new(v.dot(self.x), v.dot(self.y))
122    }
123}
124
125fn project_curve(
126    curve: &Curve,
127    range: (f64, f64),
128    onto: &Onto,
129    tol: Tolerances,
130) -> OgeomResult<ProjectedCurve> {
131    match curve {
132        Curve::Line(_) => {
133            let start = onto.point(curve.point_at(range.0, tol)?);
134            let end = onto.point(curve.point_at(range.1, tol)?);
135            Ok(if start.distance(end) <= tol.confusion() {
136                ProjectedCurve::Point(start.midpoint(end))
137            } else {
138                ProjectedCurve::Line { start, end }
139            })
140        }
141        Curve::Circle(c) => {
142            let circle = c.circle();
143            let frame = circle.frame();
144            let r = circle.radius();
145            conic(
146                onto.point(frame.origin()),
147                onto.vector(frame.x().vector() * r),
148                onto.vector(frame.y().vector() * r),
149                range,
150                tol,
151            )
152        }
153        Curve::Ellipse(e) => {
154            let ellipse = e.ellipse();
155            let frame = ellipse.frame();
156            conic(
157                onto.point(frame.origin()),
158                onto.vector(frame.x().vector() * ellipse.major_radius()),
159                onto.vector(frame.y().vector() * ellipse.minor_radius()),
160                range,
161                tol,
162            )
163        }
164        Curve::BSpline(spline) => {
165            let piece = restricted(spline, range, tol)?;
166            Ok(ProjectedCurve::BSpline {
167                curve: projected_spline(&piece, onto, tol)?,
168                fit_error: None,
169            })
170        }
171        Curve::Trimmed(trimmed) => {
172            // A trimmed curve's parameter is its basis's, run backwards
173            // where it is reversed; the projection takes no side.
174            let (s, e) = trimmed.domain();
175            let on_basis = if trimmed.is_reversed() {
176                (s + e - range.1, s + e - range.0)
177            } else {
178                range
179            };
180            project_curve(trimmed.basis(), on_basis, onto, tol)
181        }
182        _ => fitted(curve, range, onto, tol),
183    }
184}
185
186/// The projection of `centre + u cos t + v sin t` over `range`: an
187/// ellipse in general, a circle where `u` and `v` stay square and equal,
188/// a segment where the ellipse closes to no width.
189fn conic(
190    centre: Point2,
191    u: Vector2,
192    v: Vector2,
193    range: (f64, f64),
194    tol: Tolerances,
195) -> OgeomResult<ProjectedCurve> {
196    // Principal axes of the conjugate pair: at `t0` the point is furthest
197    // from the centre, and the pair turned by `t0` is square.
198    let t0 = 0.5 * (2.0 * u.dot(v)).atan2(u.dot(u) - v.dot(v));
199    let major = u * t0.cos() + v * t0.sin();
200    let minor = v * t0.cos() - u * t0.sin();
201    let (major, minor, t0) = if minor.magnitude() > major.magnitude() {
202        (minor, -major, t0 + core::f64::consts::FRAC_PI_2)
203    } else {
204        (major, minor, t0)
205    };
206    // The point is `centre + major cos(t - t0) + minor sin(t - t0)`.
207    let length = major.magnitude();
208    if length <= tol.confusion() {
209        return Ok(ProjectedCurve::Point(centre));
210    }
211    let width = minor.magnitude();
212    if width <= tol.confusion() {
213        return Ok(segment(centre, major, (range.0 - t0, range.1 - t0)));
214    }
215    // Counter-clockwise where the minor axis leads the major by a quarter
216    // turn; otherwise the parameter runs the other way round.
217    let (from, to) = if major.cross(minor) > 0.0 {
218        (range.0 - t0, range.1 - t0)
219    } else {
220        (t0 - range.1, t0 - range.0)
221    };
222    let ratio = width / length;
223    if (1.0 - ratio) * length <= tol.confusion() {
224        let angle = major.y.atan2(major.x);
225        return Ok(ProjectedCurve::Circle {
226            centre,
227            radius: length,
228            range: normalized(from + angle, to + angle),
229        });
230    }
231    Ok(ProjectedCurve::Ellipse {
232        centre,
233        major,
234        ratio,
235        range: normalized(from, to),
236    })
237}
238
239/// A range moved by whole turns to start in `[0, 2 pi)`.
240fn normalized(from: f64, to: f64) -> (f64, f64) {
241    let start = from.rem_euclid(core::f64::consts::TAU);
242    (start, start + (to - from))
243}
244
245/// The segment `centre + major cos s` covers for `s` over `range`: its
246/// ends are the extremes of `cos s`, which an arc through a multiple of pi
247/// reaches inside the range.
248fn segment(centre: Point2, major: Vector2, range: (f64, f64)) -> ProjectedCurve {
249    let (lo, hi) = (range.0.min(range.1), range.0.max(range.1));
250    let mut values = vec![lo.cos(), hi.cos()];
251    let pi = core::f64::consts::PI;
252    let mut k = (lo / pi).ceil();
253    while k * pi <= hi && values.len() < 4 {
254        values.push((k * pi).cos());
255        k += 1.0;
256    }
257    let least = values.iter().copied().fold(f64::INFINITY, f64::min);
258    let most = values.iter().copied().fold(f64::NEG_INFINITY, f64::max);
259    ProjectedCurve::Line {
260        start: centre + major * least,
261        end: centre + major * most,
262    }
263}
264
265/// The part of a spline an edge covers, split out where the edge covers
266/// less than the whole.
267fn restricted(
268    spline: &BSplineCurve,
269    range: (f64, f64),
270    tol: Tolerances,
271) -> OgeomResult<BSplineCurve> {
272    let (lo, hi) = spline.domain();
273    let reach = tol.parametric();
274    let mut piece = spline.clone();
275    if range.0 > lo + reach {
276        piece = piece.split_at(range.0, tol)?.1;
277    }
278    if range.1 < hi - reach {
279        piece = piece.split_at(range.1, tol)?.0;
280    }
281    Ok(piece)
282}
283
284/// A spline's projection: the same knots, each control point projected,
285/// each weight kept.
286fn projected_spline(spline: &BSplineCurve, onto: &Onto, tol: Tolerances) -> OgeomResult<BSpline2d> {
287    let knots: KnotVector = spline.knots().clone();
288    if spline.is_rational() {
289        // The homogeneous point `w p` maps to `w` times the projected point,
290        // which is the projection of `w p` less `w` times the origin's.
291        let control = spline
292            .control_points()
293            .iter()
294            .map(|c| {
295                let d = c.scaled.to_vector() - onto.origin.to_vector() * c.weight;
296                Weighted {
297                    scaled: Point2::new(d.dot(onto.x), d.dot(onto.y)),
298                    weight: c.weight,
299                }
300            })
301            .collect();
302        BSpline2d::rational(knots, control)
303    } else {
304        let control = spline
305            .control_points()
306            .iter()
307            .map(|c| onto.point(c.scaled))
308            .collect();
309        BSpline2d::new(knots, control, tol)
310    }
311}
312
313/// A curve with no closed-form projection, fitted through its projected
314/// points to a hundred times the confusion distance.
315fn fitted(
316    curve: &Curve,
317    range: (f64, f64),
318    onto: &Onto,
319    tol: Tolerances,
320) -> OgeomResult<ProjectedCurve> {
321    const SAMPLES: u32 = 200;
322    let mut points = Vec::with_capacity(SAMPLES as usize + 1);
323    for k in 0..=SAMPLES {
324        let t = range.0 + (range.1 - range.0) * f64::from(k) / f64::from(SAMPLES);
325        let p = onto.point(curve.point_at(t, tol)?);
326        points.push(Point::new(p.x, p.y, 0.0));
327    }
328    let fit = ogeom_geom::fit::fit_points(&points, 3, tol.confusion() * 100.0, tol)?;
329    if !fit.met {
330        ogeom_bail!(
331            Construction,
332            "the projected curve could not be fitted within {}; its best fit misses by {}",
333            tol.confusion() * 100.0,
334            fit.error
335        );
336    }
337    let flat = BSplineCurve::new(
338        fit.curve.knots().clone(),
339        fit.curve
340            .control_points()
341            .iter()
342            .map(|c| c.scaled)
343            .collect(),
344        tol,
345    )?;
346    let onto_xy = Onto {
347        origin: Point::ORIGIN,
348        x: Vector::X,
349        y: Vector::Y,
350    };
351    Ok(ProjectedCurve::BSpline {
352        curve: projected_spline(&flat, &onto_xy, tol)?,
353        fit_error: Some(fit.error),
354    })
355}
ogeom_algo/project_plane.rs

ogeom_algo/
project_plane.rs
