axiolid_construct/

revolve.rs

//! Revolution of a profile about an axis.
//!
//! The profile is flattened to rings in its own plane, then each ring point
//! is rotated about the axis in steps chosen by the chord budget. A full
//! turn closes on itself; a partial turn is capped by the profile at each
//! end.
//!
//! Pappus gives the oracle: a full revolution has volume `2*pi*R*A`, where
//! `R` is the centroid's distance from the axis and `A` the profile area.
//!
//! # How far the walls sit from the surface (#231)
//!
//! [`revolve`] bounds the distance from every point of the surface swept by
//! the GIVEN rings to its triangles by `tolerance.linear()`. Rings that
//! chord a curved profile add their own deviation on top: a caller that
//! wants the exact surface of revolution within a budget `c` flattens the
//! profile to `c_p` and revolves to `c - c_p`. The reference compiler takes
//! half each, which minimises the triangle count for a fixed sum: the steps
//! round the profile and round the axis each grow as the inverse square
//! root of their share, so their product is smallest at an even split. The
//! two shares add because a rotation is an isometry: a point of the exact
//! surface is the rotation of a profile point within `c_p` of the ring
//! polygon, so it is within `c_p` of the same rotation of the polygon.
//!
//! Each ring segment `q_a -> q_b` sweeps, between the stations at angles
//! `t` and `t + h`, the surface `R(s) q(x)` with `q(x) = (1-x) q_a + x q_b`.
//! For fixed `x` the point `R(s) q(x)` runs along an arc of radius
//! `rho(x)`, its distance from the axis, and the chord between the arc's
//! ends is the row at `x` of the bilinear patch over the quad
//! `R(t) q_a, R(t) q_b, R(t+h) q_a, R(t+h) q_b` (a rotation is linear in
//! the point). An arc of angle `h` stays within its sagitta
//! `rho (1 - cos(h/2))` of its chord, and `rho(x)` is the norm of an affine
//! function, so convex in `x` and largest at a ring vertex: every point of
//! the swept surface is within `rho_max (1 - cos(h/2))` of the patch, with
//! `rho_max` the largest vertex distance from the axis over every ring --
//! the outer equator, for a torus, never the tube centre. The patch is in
//! turn within `loft::quad_deviation` of the two triangles the
//! loft cuts it into. When the axis lies in the profile's plane every quad
//! is a planar trapezoid (its two chords are images of parallel radial
//! vectors under the same `R(t+h) - R(t)`), convex while the segment stays
//! on one side of the axis, and that term is zero: the bound is the sagitta
//! alone, whatever the sign of the Gaussian curvature (the inner side of a
//! torus is saddle-shaped, but its trapezoids are still flat). An axis off
//! the profile's plane sweeps hyperboloidal quads; their measured twist is
//! added, and the step shrinks until the sum fits.
//!
//! Every span of a revolution is the first span rotated, so the twist is
//! measured once.
//!
//! # Vertices on the axis (#363)
//!
//! A ring vertex exactly on an in-plane axis (the exact `side_of_axis` sign
//! is zero) is fixed under the rotation, so all its stations coincide. It
//! gets one mesh vertex, shared by every station and by both caps of a
//! partial turn. Its quads collapse: the one with a single axis end keeps the
//! triangle that is not degenerate, and the one with both ends on the axis (a
//! profile edge lying along it) sweeps nothing and emits nothing; an axis
//! edge of a partial turn is then an edge between the two caps. The bound
//! above is untouched. A point at radius 0 traces no circle, so the step
//! count, which the largest radius governs, is unchanged. A collapsed quad is
//! the triangle `(a, b, d)` (for `b = d`) or `(a, d, c)` (for `a = c`), which
//! the patch bounds of `loft::quad_deviation` cover exactly: `|c - a|` is
//! zero for `a = c`, and `2|d - b|` is zero for `b = d`, where the patch
//! degenerates into the triangle `(a, d, c)`. The merge changes neither. Merging moves no
//! vertex: the shared one is the profile point itself, which every station's
//! rotated copy equals up to rounding.
//!
//! In a full turn only a vertex that ends a ring edge lying along the axis is
//! accepted (sphere, cylinder and cone profiles). One whose both ring
//! neighbours are off the axis is refused (`UnsupportedInput`): two nappes of
//! the surface meet at that point, the boundary is not a two-manifold there,
//! and its mesh would be non-manifold at the vertex whatever the edge audit
//! says, as for the horn torus of #248. In a partial turn the two caps close
//! the point into a sector point, as in the exact path, and it is meshed.
//!
//! The weld is exact: it needs the axis origin and direction to have z
//! exactly zero. With z noise within the tolerance, which the in-plane check
//! forgives, a 2D-on-axis vertex is off the real axis and its orbit is a real
//! circle, so a touching profile is refused (`UnsupportedInput`) instead.
//!
//! # Vertices near the axis (#367)
//!
//! A vertex that touches the axis need not be on it to the bit: the joint of
//! two arcs is trigonometric (x of about 6e-17 on the Y axis), and almost no
//! vertex is on a tilted axis. For the exact in-plane axis a vertex within
//! `tolerance.linear()` of the axis line, measured to the line, is treated
//! as on it: it is projected onto the line and welded as above, and the
//! pinch refusal applies to it in a full turn. Moving a ring vertex by `d`
//! moves the swept surface by at most `d`, so the largest projection is
//! taken from the chord budget the steps must meet (`BudgetExceeded` when it
//! leaves none), and the bound stays the budget. The distance is a float
//! evaluation with a rounding bound; a distance within that bound of the
//! tolerance is undecidable and `UnsupportedInput`, never welded. Under z
//! noise nothing is projected, as nothing is welded. A ring left with no area
//! by the projection (every vertex in the band) is `UnsupportedInput`: it
//! would vanish from the solid or leave it empty.

use axiolid_contracts::{GeomError, GeomResult, Operation};
use axiolid_core::{Point2, Point3, Scalar, Tolerance, Vec3};
use axiolid_guarantees::Sign;
use axiolid_mesh::TriMesh;
use axiolid_predicates::arithmetic::{expansion_sign, expansion_sum};
use axiolid_predicates::two_product;

use crate::profile::Rings;

/// The most steps a revolution or a circle takes before it refuses the
/// budget, rather than meet it with a coarser mesh that breaks the bound
/// the caller asked for.
pub(crate) const MAX_STEPS: usize = 4096;

/// Rotation of `p` about the axis through `origin` along unit `dir`.
///
/// Rodrigues' formula. Written out rather than pulled from a matrix type
/// so the axis stays arbitrary: a revolution is not a Z-up operation.
pub(crate) fn rotate(p: Point3, origin: Point3, dir: Vec3, angle: Scalar) -> Point3 {
    let v = p - origin;
    let (s, c) = angle.sin_cos();
    origin + v * c + dir.cross(v) * s + dir * (dir.dot(v) * (1.0 - c))
}

/// Angular steps so the swept arc meets the chord budget.
///
/// The widest point of the profile governs: a ring point at distance `r`
/// from the axis traces a circle of that radius, and its sagitta is
/// `r(1 - cos(dtheta/2))`. Using the maximum radius means every other
/// point is sampled at least as finely. A budget that needs more than
/// [`MAX_STEPS`] steps is refused, never met by fewer.
pub(crate) fn steps(max_radius: Scalar, angle: Scalar, tol: Scalar) -> GeomResult<usize> {
    if !(tol.is_finite() && tol > 0.0) {
        return Err(GeomError::InvalidInput(format!(
            "chord budget must be positive and finite, got {tol}"
        )));
    }
    if !(max_radius.is_finite() && max_radius > 0.0) {
        return Ok(8);
    }
    let ratio = (1.0 - (tol / max_radius).min(1.0)).clamp(-1.0, 1.0);
    let per = 2.0 * ratio.acos();
    let n = (angle.abs() / per).ceil();
    // NaN (a zero step) counts as too many.
    if n.is_nan() || n > MAX_STEPS as Scalar {
        return Err(too_many_steps());
    }
    Ok((n as usize).max(2))
}

/// The refusal for a budget that needs more than [`MAX_STEPS`] steps.
pub(crate) fn too_many_steps() -> GeomError {
    GeomError::BudgetExceeded {
        resource: "revolution angular steps",
    }
}

/// The next step count to try when `n` misses the budget: a quarter more,
/// so a near miss costs little and a far one converges in a few rounds.
pub(crate) fn grow(n: usize) -> GeomResult<usize> {
    let next = n + n / 4 + 1;
    if next > MAX_STEPS {
        return Err(too_many_steps());
    }
    Ok(next)
}

/// Steps for a revolution of `rings` through `angle` whose walls stay
/// within `budget` of the surface the rings sweep (see the module notes).
fn bounded_steps(
    rings: &Rings,
    axis_origin: Point3,
    dir: Vec3,
    angle: Scalar,
    max_r: Scalar,
    budget: Scalar,
) -> GeomResult<usize> {
    let at = |t: Scalar| {
        crate::loft::place(rings, |p| {
            rotate(Point3::new(p.x, p.y, 0.0), axis_origin, dir, t)
        })
    };
    let mut n = steps(max_r, angle, budget)?;
    loop {
        let h = angle / n as Scalar;
        let sagitta = max_r * (1.0 - (0.5 * h).cos());
        // The loft walks the stations from the far end back, so each span
        // runs from angle `t + h` to `t`; every span is this one rotated.
        let twist = crate::loft::span_deviation(&at(h), &at(0.0));
        if sagitta + twist <= budget {
            return Ok(n);
        }
        n = grow(n)?;
    }
}

/// Exact sign of `d x (p - a)`, the side of `p` of the line through `a`
/// along `d`: `d.x*p.y - d.x*a.y - d.y*p.x + d.y*a.x`, summed as expansions.
///
/// Not `orient2d(a, a + d, p)`: `a + d` rounds, which tilts the line it
/// decides against.
pub(crate) fn side_of_axis(a: Point2, d: Vec3, p: Point2) -> Sign {
    let mut sum = expansion_sum(&term(d.x, p.y, false), &term(d.x, a.y, true));
    sum = expansion_sum(&sum, &term(d.y, p.x, true));
    sum = expansion_sum(&sum, &term(d.y, a.x, false));
    expansion_sign(&sum)
}

/// Exact sign of `d . (to - from)`, how far `to` lies along the axis
/// direction `d` from `from`: `d.x*to.x - d.x*from.x + d.y*to.y - d.y*from.y`,
/// summed as expansions. Zero exactly when the segment is perpendicular to
/// the axis (#450).
pub(crate) fn along_axis(d: Vec3, from: Point2, to: Point2) -> Sign {
    let mut sum = expansion_sum(&term(d.x, to.x, false), &term(d.x, from.x, true));
    sum = expansion_sum(&sum, &term(d.y, to.y, false));
    sum = expansion_sum(&sum, &term(d.y, from.y, true));
    expansion_sign(&sum)
}

/// `x * y` (or its negation) as the two-term expansion of its exact value.
fn term(x: Scalar, y: Scalar, negate: bool) -> [Scalar; 2] {
    let (product, error) = two_product(x, y);
    if negate {
        [-error, -product]
    } else {
        [error, product]
    }
}

/// The exact side of every ring vertex of the line through `a` along the
/// `x`, `y` of `d`, in flat ring order (outer ring, then each hole).
fn ring_sides(rings: &Rings, a: Point2, d: Vec3) -> GeomResult<Vec<Sign>> {
    let mut sides = Vec::new();
    for p in rings.outer.iter().chain(rings.holes.iter().flatten()) {
        if !(p.x.is_finite() && p.y.is_finite()) {
            return Err(GeomError::InvalidInput(
                "revolved profile must have finite vertices".to_owned(),
            ));
        }
        sides.push(side_of_axis(a, d, *p));
    }
    Ok(sides)
}

/// Judge an axis leaning out of the profile plane (`axis_direction.z`
/// beyond the tolerance) (#362).
///
/// A point `p` of the plane has the invariants `h`, its position along the
/// axis, and `rho`, its distance from it, and two points sweep the same
/// orbit exactly when they agree in both. `h` is linear in the plane, so
/// each level set is a line perpendicular to the axis's projection `N` into
/// the plane; along it `rho^2 = s^2 + c` is even about `N`. Two points
/// share an orbit exactly when they are mirror images across `N`, so the
/// sweep is free of self-intersection exactly when no two points of the
/// profile are, and for a region that holds no mirror pair, that is when it
/// lies on one closed side of `N`. The exact side test of the in-plane axis,
/// applied to `N` (the line through the origin's `x`, `y` along the
/// direction's), decides that, with no tolerance band: rings on both
/// strict sides are `InvalidInput`. This subsumes an axis piercing the
/// interior of the profile, which puts rings on both sides.
///
/// The axis meets the plane at one point `X` on `N`. A profile touching `N`
/// is fine, its points have `rho > 0`, unless `X` is on the profile: that
/// point is fixed under the rotation and the solid pinches there. `X` is not
/// rational in the inputs without division, so it is decided with a margin:
/// the profile's contact with `N` (its vertices on `N`, and the edges
/// between two of them) must keep clear of `X` by more than
/// `tolerance.linear()` plus a rounding allowance, measured along `N` (the
/// distance between the two points of the plane, not their orbit radius
/// about the axis), and is otherwise `UnsupportedInput`, a case that is not
/// certified. A profile strictly to one side of `N` has no contact and needs
/// no margin. An axis exactly perpendicular to the plane has no projection:
/// every copy of the profile lies in the plane and they coincide,
/// `InvalidInput`.
///
/// A tapered revolution blends `start` into `end` linearly, so every
/// station is judged: a vertex on one side of `N` at one end and on the
/// other at the other end crosses `N` between them, which the pooled side
/// test refuses; a vertex on `N` at both ends slides along it, and its
/// whole closed span of offsets, and an edge on `N` at both ends the span of
/// its four, must keep clear of `X`. A vertex or edge on `N` at one end only
/// touches it there alone. A plain revolution passes its rings as both.
///
/// No vertex is welded: a vertex on `N` is off the axis. A partial turn is
/// judged alike, as an in-plane axis is.
fn check_leaning_axis(
    start: &Rings,
    end: &Rings,
    axis_origin: Point3,
    axis_direction: Vec3,
    tolerance: Tolerance,
) -> GeomResult<()> {
    if axis_direction.x == 0.0 && axis_direction.y == 0.0 {
        return Err(GeomError::InvalidInput(
            "revolved profile lies in a plane perpendicular to the axis: its copies coincide"
                .to_owned(),
        ));
    }
    let structure = |r: &Rings| -> Vec<usize> {
        core::iter::once(r.outer.len())
            .chain(r.holes.iter().map(Vec::len))
            .collect()
    };
    let shape = structure(start);
    if shape != structure(end) {
        return Err(GeomError::InvalidInput(
            "the end rings of a revolution must share their ring structure".to_owned(),
        ));
    }
    let a = Point2::new(axis_origin.x, axis_origin.y);
    let sides_start = ring_sides(start, a, axis_direction)?;
    let sides_end = ring_sides(end, a, axis_direction)?;
    let both = |sign: Sign| sides_start.contains(&sign) || sides_end.contains(&sign);
    if both(Sign::Positive) && both(Sign::Negative) {
        return Err(GeomError::InvalidInput(
            "revolved profile crosses the axis: its rings lie on both sides of the axis's \
             projection into the profile plane"
                .to_owned(),
        ));
    }
    if !both(Sign::Zero) {
        return Ok(());
    }
    // Positions along N measured from where the axis meets the plane.
    let along = axis_direction.x.hypot(axis_direction.y);
    let t = -axis_origin.z / axis_direction.z;
    let uncertified = || GeomError::UnsupportedInput {
        backend: crate::BACKEND_ID,
        operation: Operation::Sweep,
        input: "revolution about an axis leaning out of the profile plane that may meet the \
                profile",
    };
    if !t.is_finite() || !along.is_finite() {
        return Err(uncertified());
    }
    let unit = (axis_direction.x / along, axis_direction.y / along);
    let reach = t * along;
    let offset_of = |p: Point2| (p.x - a.x) * unit.0 + (p.y - a.y) * unit.1 - reach;
    let slack = |p: Point2| {
        let scale =
            p.x.abs()
                .max(p.y.abs())
                .max(a.x.abs())
                .max(a.y.abs())
                .max(reach.abs());
        tolerance.linear() + 64.0 * Scalar::EPSILON * scale
    };
    let flat = |r: &Rings| -> Vec<Point2> {
        r.outer
            .iter()
            .chain(r.holes.iter().flatten())
            .copied()
            .collect()
    };
    let (from, to) = (flat(start), flat(end));
    // The offset along N of ring point `i` at one end, if it is on N.
    let at = |sides: &[Sign], points: &[Point2], i: usize| {
        (sides[i] == Sign::Zero).then(|| offset_of(points[i]))
    };
    let mut first = 0;
    for len in shape {
        for k in 0..len {
            let (f, g) = (first + k, first + (k + 1) % len);
            let (s0, s1) = (at(&sides_start, &from, f), at(&sides_start, &from, g));
            let (e0, e1) = (at(&sides_end, &to, f), at(&sides_end, &to, g));
            let margin = slack(from[f]).max(slack(to[f]));
            // The closed spans of offsets the vertex and the edge cover: at
            // each end alone, and sliding along N between the two.
            for group in [
                vec![s0],
                vec![e0],
                vec![s0, e0],
                vec![s0, s1],
                vec![e0, e1],
                vec![s0, s1, e0, e1],
            ] {
                // A group counts only when every member is on N.
                if group.iter().any(Option::is_none) {
                    continue;
                }
                let offsets: Vec<Scalar> = group.into_iter().flatten().collect();
                let low = offsets.iter().copied().fold(Scalar::INFINITY, Scalar::min);
                let high = offsets
                    .iter()
                    .copied()
                    .fold(Scalar::NEG_INFINITY, Scalar::max);
                if !(low.is_finite() && high.is_finite())
                    || (low - margin <= 0.0 && 0.0 <= high + margin)
                {
                    return Err(uncertified());
                }
            }
        }
        first += len;
    }
    Ok(())
}

/// [`check_axis_side`] for a tapered revolution from `start` to `end` (#364,
/// #362): the flags are over `[start.outer, end.outer, start.holes..,
/// end.holes..]`, the order its welding reads.
///
/// An axis leaning out of the plane is judged over the whole blend (see
/// [`check_leaning_axis`]) and flags nothing. Otherwise the rings are judged
/// together, each as a ring of its own, so a section that starts on one
/// side of the axis and ends on the other is refused (#351).
pub(crate) fn check_blend_axis_side(
    start: &Rings,
    end: &Rings,
    axis_origin: Point3,
    axis_direction: Vec3,
    tolerance: Tolerance,
    full: bool,
) -> GeomResult<BlendContact> {
    let ends = Rings {
        outer: start.outer.clone(),
        holes: core::iter::once(&end.outer)
            .chain(&start.holes)
            .chain(&end.holes)
            .cloned()
            .collect(),
    };
    if axis_direction.z.abs() > tolerance.linear() * axis_direction.length() {
        check_leaning_axis(start, end, axis_origin, axis_direction, tolerance)?;
        let count = ends.outer.len() + ends.holes.iter().map(Vec::len).sum::<usize>();
        return Ok(BlendContact {
            flags: vec![false; count],
            start: start.clone(),
            end: end.clone(),
            projection: 0.0,
        });
    }
    let contact = check_axis_side(&ends, axis_origin, axis_direction, tolerance, full)?;
    let mut holes = contact.rings.holes.into_iter();
    let end_outer = holes.next().unwrap_or_default();
    let start_holes: Vec<Vec<Point2>> = holes.by_ref().take(start.holes.len()).collect();
    Ok(BlendContact {
        flags: contact.on_axis,
        start: Rings {
            outer: contact.rings.outer,
            holes: start_holes,
        },
        end: Rings {
            outer: end_outer,
            holes: holes.collect(),
        },
        projection: contact.projection,
    })
}

/// What [`check_axis_side`] found: the flags of the vertices on the axis, the
/// rings with those vertices projected onto it, and the farthest a vertex
/// moved (#367).
pub(crate) struct AxisContact {
    pub(crate) on_axis: Vec<bool>,
    pub(crate) rings: Rings,
    pub(crate) projection: Scalar,
}

/// [`AxisContact`] of a tapered revolution: the flags over `[start.outer,
/// end.outer, start.holes.., end.holes..]` and the two profiles projected.
pub(crate) struct BlendContact {
    pub(crate) flags: Vec<bool>,
    pub(crate) start: Rings,
    pub(crate) end: Rings,
    pub(crate) projection: Scalar,
}

/// The chord budget left once the walls have moved by `projection` (#367):
/// moving a ring vertex by `d` moves the swept surface by at most `d`, so
/// the projection adds to the deviation bound. A projection that uses the
/// whole budget leaves none, which is `BudgetExceeded`.
pub(crate) fn budget_after_projection(
    tolerance: Tolerance,
    projection: Scalar,
) -> GeomResult<Scalar> {
    let left = tolerance.linear() - projection;
    if left > 0.0 {
        Ok(left)
    } else {
        Err(too_many_steps())
    }
}

/// Whether `p` is within `tolerance.linear()` of the in-plane axis line
/// through `a` along `d` (#367): `Some(distance bound)` when it is, `None`
/// when it is not. The distance is a float evaluation of
/// `|d x (p - a)| / |d|`, whose rounding error is bounded by eight machine
/// epsilons of the operand magnitudes; a distance within that error of the
/// tolerance cannot be decided and is `UnsupportedInput`.
fn near_axis(a: Point2, d: Vec3, p: Point2, tolerance: Tolerance) -> GeomResult<Option<Scalar>> {
    let tol = tolerance.linear();
    if tol <= 0.0 {
        return Ok(None);
    }
    let norm = d.x.hypot(d.y);
    let cross = d.x * (p.y - a.y) - d.y * (p.x - a.x);
    let distance = cross.abs() / norm;
    let error = 8.0
        * Scalar::EPSILON
        * (d.x.abs() * (p.y.abs() + a.y.abs()) + d.y.abs() * (p.x.abs() + a.x.abs()))
        / norm;
    if distance + error <= tol {
        Ok(Some(distance + error))
    } else if distance - error > tol {
        Ok(None)
    } else {
        Err(touching_unsupported())
    }
}

/// How far apart `p` and `q`, both on the in-plane axis, are along it:
/// `Some(distance)` when it is at most `tolerance.linear()` with room for the
/// rounding of the two evaluations, else `None` (a slide, or undecidable
/// within rounding, both refused by the caller) (#367).
pub(crate) fn slide_along_axis(
    axis_origin: Point3,
    d: Vec3,
    p: Point2,
    q: Point2,
    tolerance: Tolerance,
) -> Option<Scalar> {
    if p == q {
        return Some(0.0);
    }
    let norm = d.x.hypot(d.y);
    let along = |r: Point2| (d.x * (r.x - axis_origin.x) + d.y * (r.y - axis_origin.y)) / norm;
    let error = |r: Point2| {
        8.0 * Scalar::EPSILON
            * (d.x.abs() * (r.x.abs() + axis_origin.x.abs())
                + d.y.abs() * (r.y.abs() + axis_origin.y.abs()))
            / norm
    };
    let shift = (along(p) - along(q)).abs() + error(p) + error(q);
    (shift <= tolerance.linear()).then_some(shift)
}

/// `p` projected onto the in-plane axis line through `a` along `d`.
fn project_onto_axis(a: Point2, d: Vec3, p: Point2) -> Point2 {
    let t = (d.x * (p.x - a.x) + d.y * (p.y - a.y)) / (d.x * d.x + d.y * d.y);
    Point2::new(a.x + d.x * t, a.y + d.y * t)
}

/// Refuse rings that lie on both sides of the axis, and report which touch
/// it (#351, #363).
///
/// Decided when the axis lies in the profile's plane, where "side" is
/// meaningful: each vertex is classified by the exact sign of its cross
/// product with the axis direction, with no tolerance band. Vertices
/// strictly on both sides are a crossing (the swept solid passes through
/// itself) and `InvalidInput`. A vertex on the axis that ends an edge lying
/// on it is not an error: the result holds `true` at its index in the flat
/// ring order (outer ring, then each hole) and [`revolve`] meshes it once. A
/// vertex on the axis with both ring neighbours off it pinches the solid to
/// a point and, in a `full` turn, is `UnsupportedInput`; a partial turn's
/// caps close it into a sector point. An edge cannot meet the axis between
/// two vertices without one of those cases. A vertex on the projected axis
/// is `UnsupportedInput` too when the axis origin or direction has any
/// nonzero z: the vertex is then off the real axis by the noise, and welding
/// it would collapse a real orbit. For the exact axis a vertex within
/// `tolerance.linear()` of the line counts as on it and is projected onto it
/// (#367, see the module notes); the result carries the moved rings.
///
/// "In the plane" follows the exact path: the axis origin and direction
/// leave it only beyond `tolerance.linear()` in z, so placement noise does
/// not switch the check off. An axis parallel to the plane but offset from
/// it is checked for crossing only: no ring point reaches the axis, since
/// its radius stays at least the offset, so nothing is reported. An axis
/// leaning out of the plane is judged by [`check_leaning_axis`] (#362).
pub(crate) fn check_axis_side(
    rings: &Rings,
    axis_origin: Point3,
    axis_direction: Vec3,
    tolerance: Tolerance,
    full: bool,
) -> GeomResult<AxisContact> {
    let count = rings.outer.len() + rings.holes.iter().map(Vec::len).sum::<usize>();
    if axis_direction.z.abs() > tolerance.linear() * axis_direction.length() {
        check_leaning_axis(rings, rings, axis_origin, axis_direction, tolerance)?;
        return Ok(AxisContact {
            on_axis: vec![false; count],
            rings: rings.clone(),
            projection: 0.0,
        });
    }
    let offset = axis_origin.z.abs() > tolerance.linear();
    let exact = axis_origin.z == 0.0 && axis_direction.z == 0.0;
    let a = Point2::new(axis_origin.x, axis_origin.y);
    let (mut positive, mut negative, mut noisy_touch) = (false, false, false);
    let mut on_axis = Vec::with_capacity(count);
    let mut snapped = Vec::with_capacity(count);
    let mut moved = Vec::with_capacity(count);
    let mut projection: Scalar = 0.0;
    let points = rings.outer.iter().chain(rings.holes.iter().flatten());
    for (side, p) in ring_sides(rings, a, axis_direction)?
        .into_iter()
        .zip(points)
    {
        let (side, p, near) = if side != Sign::Zero && exact && !offset {
            match near_axis(a, axis_direction, *p, tolerance)? {
                Some(distance) => {
                    projection = projection.max(distance);
                    (Sign::Zero, project_onto_axis(a, axis_direction, *p), true)
                }
                None => (side, *p, false),
            }
        } else {
            (side, *p, false)
        };
        snapped.push(p);
        moved.push(near);
        match side {
            Sign::Positive => positive = true,
            Sign::Negative => negative = true,
            _ => {}
        }
        noisy_touch |= side == Sign::Zero && !offset && !exact;
        on_axis.push(side == Sign::Zero && !offset);
    }
    if positive && negative {
        return Err(GeomError::InvalidInput(
            "revolved profile crosses the axis: its rings lie on both sides".to_owned(),
        ));
    }
    if noisy_touch {
        return Err(touching_unsupported());
    }
    // In a full turn a vertex on the axis whose ring neighbours are both off
    // it is a pinch: two nappes of the swept surface meet at one point, so
    // the boundary is not a two-manifold there (the horn torus of #248 is
    // refused alike). A partial turn's caps close it into a sector point.
    let rings = {
        let mut flat = snapped.into_iter();
        let mut take =
            |ring: &Vec<Point2>| -> Vec<Point2> { flat.by_ref().take(ring.len()).collect() };
        Rings {
            outer: take(&rings.outer),
            holes: rings.holes.iter().map(&mut take).collect(),
        }
    };
    // A ring that the projection flattens onto the axis sweeps nothing: it
    // would vanish from the solid (a hole) or leave it empty (an outer ring).
    let mut first = 0;
    for ring in core::iter::once(&rings.outer).chain(&rings.holes) {
        let end = first + ring.len();
        if moved[first..end].contains(&true) {
            let area: Scalar = (0..ring.len())
                .map(|k| {
                    let (p, q) = (ring[k], ring[(k + 1) % ring.len()]);
                    p.x * q.y - q.x * p.y
                })
                .sum();
            if on_axis[first..end].iter().all(|&on| on) || area == 0.0 {
                return Err(flattened_unsupported());
            }
        }
        first = end;
    }
    let contact = |on_axis| AxisContact {
        on_axis,
        rings: rings.clone(),
        projection,
    };
    if !full {
        return Ok(contact(on_axis));
    }
    let mut start = 0;
    for len in core::iter::once(rings.outer.len()).chain(rings.holes.iter().map(Vec::len)) {
        let ring = &on_axis[start..start + len];
        for k in 0..len {
            if ring[k] && !ring[(k + len - 1) % len] && !ring[(k + 1) % len] {
                return Err(GeomError::UnsupportedInput {
                    backend: crate::BACKEND_ID,
                    operation: Operation::Sweep,
                    input: "revolution of a profile that touches the axis at an isolated \
                            point: the solid pinches there",
                });
            }
        }
        start += len;
    }
    Ok(contact(on_axis))
}

/// The refusal for a ring that lies within the tolerance of the axis (#367).
fn flattened_unsupported() -> GeomError {
    GeomError::UnsupportedInput {
        backend: crate::BACKEND_ID,
        operation: Operation::Sweep,
        input: "revolution of a ring within the tolerance of the axis: it sweeps no area",
    }
}

/// The refusal for a family that cannot yet mesh a profile touching the
/// axis (#363).
pub(crate) fn touching_unsupported() -> GeomError {
    GeomError::UnsupportedInput {
        backend: crate::BACKEND_ID,
        operation: Operation::Sweep,
        input: "revolution of a profile touching the axis",
    }
}

/// Revolve a profile about an axis into a closed solid.
///
/// A full turn wraps its rings by index, so the seam shares vertices by
/// construction rather than by two samplings agreeing numerically. That is
/// the same structure `tessellate_primitive` uses for a cylinder, and it
/// avoids the trim-based path recorded as broken in issue #2.
///
/// Every point of the surface the rings sweep lies within
/// `tolerance.linear()` of the mesh (see the module notes for the proof);
/// the rings' own distance from a curved profile is the caller's to add.
/// A budget that needs more than 4096 steps round the axis is refused with
/// `BudgetExceeded`. A profile with vertices on both sides of an in-plane
/// axis is `InvalidInput`. A vertex on the axis is meshed once, the point
/// being fixed under the rotation (see the module notes).
pub fn revolve(
    rings: &Rings,
    axis_origin: Point3,
    axis_direction: Vec3,
    angle: Scalar,
    tolerance: Tolerance,
) -> GeomResult<TriMesh> {
    if !angle.is_finite() || angle == 0.0 {
        return Err(GeomError::InvalidInput(format!(
            "revolution angle must be finite and non-zero, got {angle}"
        )));
    }
    let len = axis_direction.length();
    if !len.is_finite() || len <= 0.0 {
        return Err(GeomError::InvalidInput(
            "revolution axis must be a finite non-zero direction".to_owned(),
        ));
    }
    let full = (angle.abs() - core::f64::consts::TAU).abs() <= 1e-9;
    let contact = check_axis_side(rings, axis_origin, axis_direction, tolerance, full)?;
    let (on_axis, rings, projection) = (contact.on_axis, &contact.rings, contact.projection);
    let dir = axis_direction / len;
    // Widest distance from the axis governs the angular step: a point at
    // radius r traces a circle of that radius.
    let mut max_r: Scalar = 0.0;
    for p in rings.outer.iter().chain(rings.holes.iter().flatten()) {
        let v = Point3::new(p.x, p.y, 0.0) - axis_origin;
        max_r = max_r.max((v - dir * dir.dot(v)).length());
    }
    let budget = budget_after_projection(tolerance, projection)?;
    let n = bounded_steps(rings, axis_origin, dir, angle, max_r, budget)?;
    // A full turn emits n stations and wraps onto the first; a partial turn
    // emits n + 1 so both ends exist to be capped.
    let count = if full { n } else { n + 1 };
    let stations: Vec<crate::loft::Station> = (0..count)
        .rev()
        .map(|s| {
            let t = angle * (s as Scalar) / (n as Scalar);
            crate::loft::place(rings, |p| {
                rotate(Point3::new(p.x, p.y, 0.0), axis_origin, dir, t)
            })
        })
        .collect();
    let mut mesh = crate::loft::loft_welded(rings, &stations, full, &on_axis)?;
    // The station order above (end of the sweep first) winds the walls
    // outward for a positive angle about `dir`. A negative angle sweeps the
    // other way round, so the same order winds them inward instead (#221):
    // `angle` and `dir` only matter through their product (the axis-angle
    // rotation vector), and `(dir, angle)` and `(-dir, -angle)` describe the
    // identical rotation, so the fix cannot key off `angle`'s stored sign
    // alone without also depending on `dir`'s. Settling it from the built
    // mesh's own orientation sidesteps that: whichever combination produced
    // it, a closed solid's sign is unambiguous.
    if matches!(crate::extrude::outward_orientation(&mesh), Some(false)) {
        for triangle in mesh.indices.as_chunks_mut::<3>().0 {
            triangle.swap(1, 2);
        }
    }
    Ok(mesh)
}