axiolid_construct/profile/
paired.rs

1//! Profiles that change between stations, flattened together (#361).
2//!
3//! [`super::profile_rings`] flattens each profile on its own, so two
4//! different curved profiles get unrelated vertex sets and no
5//! correspondence between them: a run lofted from one to the other has no
6//! bound on its profile share. Here the profiles whose rings pair segment
7//! for segment are flattened *together*: segment `j` of ring `r` is cut
8//! into the same number of equal steps of its own parameter in every
9//! profile, and vertex `m` of one is matched to vertex `m` of the others.
10//! The loft between two sections interpolates those matched vertices, which
11//! is the chord polygon of the interpolated exact section sampled at the
12//! same fractions, so its distance from that section is a bound that holds
13//! for every interpolated section, not only the stations.
14//!
15//! # Lowering
16//!
17//! A profile lowers to segments: a rectangle, a section or a contour
18//! through its exact contour; a circle to one full-turn circle (its bore,
19//! if any, another); an ellipse to one full-turn ellipse; a derived profile
20//! to its basis plus the transforms applied to the flattened points. Each
21//! ring is oriented by the sign of its exact area (outer counter-clockwise,
22//! holes clockwise) by reversing its segments. A ring reversed in one
23//! profile and not in another (a mirrored profile against an unmirrored
24//! one) does not pair: its start point stays mirrored, so the vertices would
25//! match points of the opposite side. Rings and segments are then matched
26//! by index. Anything else (a polyline or B-spline segment, a centre line, a composite,
27//! another ring or segment count) does not pair, and the caller keeps the
28//! unbounded result it names.
29//!
30//! # The bound
31//!
32//! With `m2` the supremum of `|C''|` of a segment over its parameter (zero
33//! for a line, `r * hypot(|x|, |y|)` for a circle, `hypot(a |x|, b |y|)`
34//! for an ellipse, times the largest stretch of a derived profile's
35//! transforms) and `dt` its parameter span, a segment cut into `n` steps
36//! lies within `(dt / n)^2 * m2 / 8` of its chords, and so does every
37//! convex combination of two such segments. `n` is the smallest count that
38//! keeps that within the chord error, taken over all profiles. Joints are
39//! snapped as [`super::profile_rings`] does, and the distance a joint moved
40//! (times the stretch) is added.
41//!
42//! A run whose every profile is a circle starts each ring half a step off
43//! the axes, as the circle flattener does (#194).
44
45use axiolid_contracts::{GeomError, GeomResult};
46use axiolid_core::{Frame2, Interval, Point2, Scalar, Tolerance, Transform2, Vec2};
47use axiolid_curve::{Circle2, Curve2, Ellipse2};
48use axiolid_profile::{Contour, ContourProfile, Profile, ProfileSegment};
49use axiolid_reference::curve::evaluate2;
50
51use super::{near2, snap_joints, stretch2, Rings};
52
53/// The most vertices the pairing produces over all rings of one profile.
54const MAX_VERTICES: Scalar = (1u32 << 20) as Scalar;
55
56/// Rings flattened together, and how far the interpolated sections lie from
57/// them.
58#[non_exhaustive]
59#[derive(Debug, Clone, PartialEq)]
60pub struct PairedRings {
61    /// One [`Rings`] per profile, in the order given, matching vertex for
62    /// vertex.
63    pub rings: Vec<Rings>,
64    /// Every point of the exact section interpolated between two of the
65    /// profiles, at any fraction, is within this of the polygon the
66    /// interpolated rings form, and the reverse.
67    pub deviation: Scalar,
68}
69
70/// A profile lowered to segments.
71struct Lowered {
72    /// Outer ring first, then the holes.
73    rings: Vec<Vec<ProfileSegment>>,
74    /// Applied to the flattened points, innermost first.
75    transforms: Vec<Transform2>,
76    /// Product of the transforms' largest stretches.
77    stretch: Scalar,
78    /// Whether [`orient`] reversed each ring.
79    reversed: Vec<bool>,
80}
81
82/// Flatten `profiles` together under `chord_error`, or `None` when they do
83/// not pair segment for segment (see the module docs).
84///
85/// # Errors
86///
87/// `InvalidInput` for a chord error that is not positive and finite or a
88/// profile that does not validate; as the exact lowering of the profile
89/// otherwise.
90pub fn paired_profile_rings(
91    profiles: &[&Profile],
92    chord_error: Scalar,
93    tolerance: Tolerance,
94) -> GeomResult<Option<PairedRings>> {
95    if !chord_error.is_finite() || chord_error <= 0.0 {
96        return Err(GeomError::InvalidInput(format!(
97            "paired profile chord error must be positive and finite, got {chord_error}"
98        )));
99    }
100    if profiles.is_empty() {
101        return Ok(None);
102    }
103    let mut lowered = Vec::with_capacity(profiles.len());
104    for profile in profiles {
105        let Some(mut one) = lower(profile)? else {
106            return Ok(None);
107        };
108        if !orient(&mut one)? {
109            return Ok(None);
110        }
111        lowered.push(one);
112    }
113    // A ring reversed in one profile and not in another would pair each
114    // vertex with a point of the opposite side (a mirrored start point
115    // stays mirrored when the ring is reversed).
116    if lowered.iter().any(|l| l.reversed != lowered[0].reversed) {
117        return Ok(None);
118    }
119    let first = &lowered[0];
120    if lowered.iter().any(|l| {
121        l.rings.len() != first.rings.len()
122            || l.rings
123                .iter()
124                .zip(&first.rings)
125                .any(|(a, b)| a.len() != b.len())
126    }) {
127        return Ok(None);
128    }
129    // The common count of every segment, and the vertices it makes.
130    let mut counts: Vec<Vec<usize>> = Vec::with_capacity(first.rings.len());
131    let mut total: Scalar = 0.0;
132    for (r, ring) in first.rings.iter().enumerate() {
133        let mut row = Vec::with_capacity(ring.len());
134        for j in 0..ring.len() {
135            let mut n: Scalar = 1.0;
136            for l in &lowered {
137                let segment = &l.rings[r][j];
138                let Some(m2) = second_bound(&segment.curve) else {
139                    return Ok(None);
140                };
141                let steps = span(segment)? * (m2 * l.stretch / (8.0 * chord_error)).sqrt();
142                n = n.max(steps.ceil());
143            }
144            total += n;
145            if total.is_nan() || total > MAX_VERTICES {
146                return Ok(None);
147            }
148            row.push(n as usize);
149        }
150        counts.push(row);
151    }
152    // Every profile a circle: a half step off the axes (#194).
153    let circles = profiles.iter().all(|p| matches!(p, Profile::Circle(_)));
154    let mut rings = Vec::with_capacity(lowered.len());
155    let mut deviation: Scalar = 0.0;
156    for l in &lowered {
157        let mut moved: Scalar = 0.0;
158        let mut chord: Scalar = 0.0;
159        let mut flat: Vec<Vec<Point2>> = Vec::with_capacity(l.rings.len());
160        for (ring, row) in l.rings.iter().zip(&counts) {
161            let mut pieces: Vec<Vec<Point2>> = Vec::with_capacity(ring.len());
162            for (segment, &n) in ring.iter().zip(row) {
163                let step = span(segment)? / n as Scalar;
164                let m2 = second_bound(&segment.curve).unwrap_or(0.0) * l.stretch;
165                chord = chord.max(step * step * m2 / 8.0);
166                let shift = if circles {
167                    core::f64::consts::PI / n as Scalar
168                } else {
169                    0.0
170                };
171                pieces.push(sample(segment, n, shift)?);
172            }
173            let contour = Contour::new(ring.clone());
174            snap_joints(&contour, &mut pieces, tolerance, &mut moved, true)?;
175            let mut points = Vec::new();
176            for mut piece in pieces {
177                piece.pop();
178                points.extend(piece);
179            }
180            transform(&mut points, &l.transforms);
181            flat.push(points);
182        }
183        deviation = deviation.max(chord + moved * l.stretch);
184        if flat.iter().any(|ring| !distinct(ring, tolerance)) {
185            return Ok(None);
186        }
187        let mut it = flat.into_iter();
188        let outer = it.next().unwrap_or_default();
189        rings.push(Rings {
190            outer,
191            holes: it.collect(),
192        });
193    }
194    Ok(Some(PairedRings {
195        rings,
196        deviation: deviation * (1.0 + 1e-12),
197    }))
198}
199
200fn transform(points: &mut [Point2], transforms: &[Transform2]) {
201    for p in points {
202        for t in transforms {
203            *p = t.transform_point2(*p);
204        }
205    }
206}
207
208/// Whether a closed ring has at least three vertices, none repeating the
209/// one before it (nor the last the first) within the linear tolerance.
210fn distinct(ring: &[Point2], tolerance: Tolerance) -> bool {
211    ring.len() >= 3
212        && (0..ring.len()).all(|i| !near2(ring[i], ring[(i + 1) % ring.len()], tolerance.linear()))
213}
214
215/// A segment's parameter span.
216fn span(segment: &ProfileSegment) -> GeomResult<Scalar> {
217    let span = (segment.domain.end - segment.domain.start).abs();
218    if span.is_finite() {
219        Ok(span)
220    } else {
221        Err(GeomError::InvalidInput(
222            "a profile segment's parameter domain must be finite".into(),
223        ))
224    }
225}
226
227/// `sup |C''|` over a segment's parameter: a line has none, a circle
228/// `r * hypot(|x|, |y|)`, an ellipse `hypot(a |x|, b |y|)`; other families
229/// do not pair.
230fn second_bound(curve: &Curve2) -> Option<Scalar> {
231    match curve {
232        Curve2::Line(_) => Some(0.0),
233        Curve2::Circle(c) => Some(c.radius * c.frame.x.length().hypot(c.frame.y.length())),
234        Curve2::Ellipse(e) => {
235            Some((e.semi_axis_x * e.frame.x.length()).hypot(e.semi_axis_y * e.frame.y.length()))
236        }
237        _ => None,
238    }
239}
240
241/// `count + 1` points of a segment in contour order, at equal fractions of
242/// its parameter span from its start in that order, the domain moved by
243/// `shift` first.
244fn sample(segment: &ProfileSegment, count: usize, shift: Scalar) -> GeomResult<Vec<Point2>> {
245    let (lo, hi) = (segment.domain.start + shift, segment.domain.end + shift);
246    let (from, to) = if segment.same_sense {
247        (lo, hi)
248    } else {
249        (hi, lo)
250    };
251    (0..=count)
252        .map(|m| {
253            let t = if m == 0 {
254                from
255            } else if m == count {
256                to
257            } else {
258                from + (to - from) * (m as Scalar / count as Scalar)
259            };
260            evaluate2(&segment.curve, t)
261        })
262        .collect()
263}
264
265/// Reverse every ring whose area has the wrong sign: the outer
266/// counter-clockwise, holes clockwise. `false` when a ring's area is zero or
267/// not finite, so its orientation is not known.
268fn orient(lowered: &mut Lowered) -> GeomResult<bool> {
269    let det: Scalar = lowered
270        .transforms
271        .iter()
272        .map(|t| t.matrix2.determinant())
273        .product();
274    for (r, ring) in lowered.rings.iter_mut().enumerate() {
275        let mut area = 0.0;
276        for segment in ring.iter() {
277            area += segment_area(segment)?;
278        }
279        area *= det;
280        if area == 0.0 || !area.is_finite() {
281            return Ok(false);
282        }
283        let reverse = (r == 0 && area < 0.0) || (r > 0 && area > 0.0);
284        if reverse {
285            ring.reverse();
286            for segment in ring.iter_mut() {
287                segment.same_sense = !segment.same_sense;
288            }
289        }
290        lowered.reversed.push(reverse);
291    }
292    Ok(true)
293}
294
295/// `1/2 * integral of (x dy - y dx)` along a segment in contour order: the
296/// exact signed area it contributes to a closed ring. A line is the cross
297/// product of its end points; a circle or ellipse `O + A cos t + B sin t`
298/// gives `(A x B)(t1 - t0) + O x B (sin t1 - sin t0) + O x A (cos t1 - cos t0)`.
299fn segment_area(segment: &ProfileSegment) -> GeomResult<Scalar> {
300    let (t0, t1) = (segment.domain.start, segment.domain.end);
301    let cross = |a: Vec2, b: Vec2| a.x * b.y - a.y * b.x;
302    let conic = |o: Point2, a: Vec2, b: Vec2| {
303        let o = Vec2::new(o.x, o.y);
304        0.5 * (cross(a, b) * (t1 - t0)
305            + cross(o, b) * (t1.sin() - t0.sin())
306            + cross(o, a) * (t1.cos() - t0.cos()))
307    };
308    let forward = match &segment.curve {
309        Curve2::Line(_) => {
310            let p = evaluate2(&segment.curve, t0)?;
311            let q = evaluate2(&segment.curve, t1)?;
312            0.5 * (p.x * q.y - q.x * p.y)
313        }
314        Curve2::Circle(c) => conic(c.frame.origin, c.frame.x * c.radius, c.frame.y * c.radius),
315        Curve2::Ellipse(e) => conic(
316            e.frame.origin,
317            e.frame.x * e.semi_axis_x,
318            e.frame.y * e.semi_axis_y,
319        ),
320        _ => return Ok(Scalar::NAN),
321    };
322    Ok(if segment.same_sense {
323        forward
324    } else {
325        -forward
326    })
327}
328
329/// A profile as segments, or `None` for a family that does not pair.
330fn lower(profile: &Profile) -> GeomResult<Option<Lowered>> {
331    let plain = |rings: Vec<Vec<ProfileSegment>>| Lowered {
332        rings,
333        transforms: Vec::new(),
334        stretch: 1.0,
335        reversed: Vec::new(),
336    };
337    let contour = |c: &ContourProfile| {
338        plain(
339            core::iter::once(&c.outer)
340                .chain(&c.holes)
341                .map(|contour| contour.segments.clone())
342                .collect(),
343        )
344    };
345    Ok(match profile {
346        Profile::Rectangle(r) => Some(contour(&crate::section_lower::rectangle_contour(r)?)),
347        Profile::Section(s) => Some(contour(&crate::section_lower::section_contour(s)?)),
348        Profile::Contour(c) => Some(contour(c)),
349        Profile::Circle(c) => {
350            if !c.radius.is_finite() || c.radius <= 0.0 {
351                return Err(GeomError::InvalidInput(format!(
352                    "circle profile radius must be positive and finite, got {}",
353                    c.radius
354                )));
355            }
356            let circle = |radius| {
357                full_turn(Curve2::Circle(Circle2 {
358                    frame: unit_frame(),
359                    radius,
360                }))
361            };
362            let mut rings = vec![vec![circle(c.radius)]];
363            if let Some(t) = c.thickness {
364                if !(t > 0.0 && t < c.radius) {
365                    return Err(GeomError::InvalidInput(format!(
366                        "annulus wall thickness {t} does not fit inside radius {}",
367                        c.radius
368                    )));
369                }
370                rings.push(vec![circle(c.radius - t)]);
371            }
372            Some(plain(rings))
373        }
374        Profile::Ellipse(e) => {
375            if !(e.semi_axis_x.is_finite() && e.semi_axis_x > 0.0)
376                || !(e.semi_axis_y.is_finite() && e.semi_axis_y > 0.0)
377            {
378                return Err(GeomError::InvalidInput(format!(
379                    "ellipse semi-axes must be positive and finite, got {} and {}",
380                    e.semi_axis_x, e.semi_axis_y
381                )));
382            }
383            Some(plain(vec![vec![full_turn(Curve2::Ellipse(Ellipse2 {
384                frame: unit_frame(),
385                semi_axis_x: e.semi_axis_x,
386                semi_axis_y: e.semi_axis_y,
387            }))]]))
388        }
389        Profile::Derived { basis, transform } => lower(basis)?.map(|mut inner| {
390            inner.transforms.push(*transform);
391            inner.stretch *= stretch2(transform);
392            inner
393        }),
394        _ => None,
395    })
396}
397
398fn unit_frame() -> Frame2 {
399    Frame2 {
400        origin: Point2::new(0.0, 0.0),
401        x: Vec2::new(1.0, 0.0),
402        y: Vec2::new(0.0, 1.0),
403    }
404}
405
406fn full_turn(curve: Curve2) -> ProfileSegment {
407    ProfileSegment {
408        curve,
409        domain: Interval {
410            start: 0.0,
411            end: core::f64::consts::TAU,
412        },
413        same_sense: true,
414    }
415}
axiolid_construct/profile/paired.rs

axiolid_construct/profile/
paired.rs
