ogeom_fillet/analyse.rs
1//! What a blend actually achieved, measured rather than asserted.
2//!
3//! A blend claims two things: that it meets each of its supports along a
4//! curve, and that it meets them *smoothly*: the two surfaces sharing a
5//! normal there. Both are claims about geometry that construction can get
6//! subtly wrong: a fitted section drifts, a marched spine carries its
7//! chord budget, a rebuilt face lands on a neighbour a hair off. This
8//! module reports the two numbers instead of trusting them.
9//!
10//! The report is per shared edge, because that is where the claim lives.
11
12use ogeom_core::{OgeomResult, Tolerances, ogeom_bail};
13use ogeom_geom::{Curve3d as _, Surface as _};
14use ogeom_topo::{EdgeRepr, Model, NodeData, Shape, ShapeType, explore_unique};
15
16/// How well a blend face meets one neighbour along one edge.
17#[derive(Debug, Clone, PartialEq)]
18pub struct BlendContact {
19 /// The shared edge.
20 pub edge: Shape,
21 /// The face on the other side of it.
22 pub neighbour: Shape,
23 /// The largest angle, in radians, between the two surfaces' normals at
24 /// the sampled stations: zero for a tangent join.
25 pub tangency_error: f64,
26 /// The largest distance, in model units, between the edge's curve and
27 /// the two surfaces it is supposed to lie on.
28 pub gap: f64,
29 /// How many stations were sampled along the edge.
30 pub stations: usize,
31}
32
33/// Measure a blend face against every face it shares an edge with.
34///
35/// Both numbers are sampled: `stations` points spread over each shared
36/// edge's range, the surfaces read at the parameters the edge's own pcurves
37/// give, so nothing is inverted and nothing is guessed. A blend that is
38/// tangent everywhere except between two stations reports zero, which is
39/// the honest limit of sampling and the reason the count is in the report.
40///
41/// An edge whose pcurve on either face is missing cannot be measured this
42/// way at all: it is reported with its `gap` and `tangency_error` set to
43/// infinity rather than quietly skipped, because a blend nobody can measure
44/// is not a blend anybody should trust.
45///
46/// # Errors
47///
48/// [`OgeomError::Construction`](ogeom_core::OgeomError::Construction) if
49/// `blend` is not a face of `shape`, or `stations` is less than two.
50pub fn analyse_blend(
51 model: &Model,
52 shape: &Shape,
53 blend: &Shape,
54 stations: usize,
55 tol: Tolerances,
56) -> OgeomResult<Vec<BlendContact>> {
57 if stations < 2 {
58 ogeom_bail!(
59 Construction,
60 "a blend measured at {stations} stations is not measured"
61 );
62 }
63 let faces = explore_unique(model, shape, ShapeType::Face)?;
64 if !faces.iter().any(|f| f.node() == blend.node()) {
65 ogeom_bail!(Construction, "that face is not part of this shape");
66 }
67 let surface_of =
68 |face: &Shape| -> Option<(ogeom_topo::SurfaceId, ogeom_geom::SurfaceGeometry)> {
69 let NodeData::Face(data) = model.node(face)?.data() else {
70 return None;
71 };
72 let geometry = model.geometry().surface(data.surface)?.clone();
73 Some((data.surface, geometry))
74 };
75 let Some((blend_id, blend_surface)) = surface_of(blend) else {
76 ogeom_bail!(Construction, "the blend face carries no surface");
77 };
78 let blend_edges = explore_unique(model, blend, ShapeType::Edge)?;
79
80 let mut out = Vec::new();
81 for neighbour in &faces {
82 if neighbour.node() == blend.node() {
83 continue;
84 }
85 let Some((other_id, other_surface)) = surface_of(neighbour) else {
86 continue;
87 };
88 for edge in explore_unique(model, neighbour, ShapeType::Edge)? {
89 if !blend_edges.iter().any(|e| e.node() == edge.node()) {
90 continue;
91 }
92 let Some(data) = model.node(&edge).and_then(|n| n.data().as_edge()) else {
93 continue;
94 };
95 let Some(EdgeRepr::Curve3d { curve, range, .. }) = data.curve3d() else {
96 continue;
97 };
98 let Some(geometry) = model.geometry().curve(*curve).cloned() else {
99 continue;
100 };
101 let chart =
102 |id: ogeom_topo::SurfaceId| -> Option<(ogeom_geom::PlanarCurve, (f64, f64))> {
103 match data.pcurve_for(id, edge.location())? {
104 EdgeRepr::PCurve { curve, range, .. } => {
105 Some((model.geometry().pcurve(*curve)?.clone(), *range))
106 }
107 EdgeRepr::Seam { forward, range, .. } => {
108 Some((model.geometry().pcurve(*forward)?.clone(), *range))
109 }
110 _ => None,
111 }
112 };
113 let (Some((pc_blend, pr_blend)), Some((pc_other, pr_other))) =
114 (chart(blend_id), chart(other_id))
115 else {
116 out.push(BlendContact {
117 edge: edge.clone(),
118 neighbour: neighbour.clone(),
119 tangency_error: f64::INFINITY,
120 gap: f64::INFINITY,
121 stations: 0,
122 });
123 continue;
124 };
125 let (mut worst_angle, mut worst_gap) = (0.0f64, 0.0f64);
126 for k in 0..stations {
127 #[expect(
128 clippy::cast_precision_loss,
129 reason = "a station index, far below the mantissa"
130 )]
131 let f = k as f64 / (stations - 1) as f64;
132 let t = (range.1 - range.0).mul_add(f, range.0);
133 let on_curve = geometry.point_at(t, tol)?;
134 let mut normals = Vec::with_capacity(2);
135 for ((pcurve, prange), surface) in [
136 ((&pc_blend, pr_blend), &blend_surface),
137 ((&pc_other, pr_other), &other_surface),
138 ] {
139 let pt = (prange.1 - prange.0).mul_add(f, prange.0);
140 let uv = ogeom_geom::Curve2d::point_at(pcurve, pt, tol)?;
141 // A pcurve ending on a chart's pole may stand a rounding
142 // past the chart's window there; the station is read at
143 // the window's edge, and the gap says what that costs.
144 let uv = {
145 use ogeom_geom::Surface as _;
146 let ((u0, u1), (v0, v1)) = surface.domain();
147 ogeom_math::Point2::new(
148 if surface.is_periodic_u() {
149 uv.x
150 } else {
151 uv.x.clamp(u0, u1)
152 },
153 if surface.is_periodic_v() {
154 uv.y
155 } else {
156 uv.y.clamp(v0, v1)
157 },
158 )
159 };
160 worst_gap =
161 worst_gap.max(surface.point_at(uv.x, uv.y, tol)?.distance(on_curve));
162 // A station on a chart's pole (a corner patch's own
163 // corner sits on the ball's pole by construction) has
164 // no normal from the chart; the stations beside it say
165 // what the join does there.
166 if let Ok(normal) = surface.normal_at(uv.x, uv.y, tol) {
167 normals.push(normal.vector());
168 }
169 }
170 if normals.len() < 2 {
171 continue;
172 }
173 // Orientation is the topology's business, not the join's:
174 // two faces meeting smoothly may still be wound opposite
175 // ways, so the angle is taken to the nearer sense.
176 let dot = normals[0].dot(normals[1]).abs().min(1.0);
177 worst_angle = worst_angle.max(dot.acos());
178 }
179 out.push(BlendContact {
180 edge: edge.clone(),
181 neighbour: neighbour.clone(),
182 tangency_error: worst_angle,
183 gap: worst_gap,
184 stations,
185 });
186 }
187 }
188 Ok(out)
189}