ogeom_geom/

traits.rs

//! The adaptor traits: what every algorithm sees geometry through.
//!
//! Intersection, projection, extrema and tessellation are written against
//! [`Curve3d`], [`Curve2d`] and [`Surface`], never against a concrete type.
//! Nothing downstream needs to know whether it holds an analytic cylinder or a
//! NURBS patch, which is what keeps one implementation of each algorithm rather
//! than one per surface type.
//!
//! # Fast paths, opted into
//!
//! [`Curve3d::kind`] and [`Surface::kind`] let an algorithm *ask*. Plane/plane
//! intersection is two lines of algebra and should never go through a marching
//! intersector. A caller that wants that shortcut matches on the kind and takes
//! it. The general path never has to know what it is looking at, so adding a
//! surface type does not break existing algorithms. It only forgoes a
//! shortcut until someone writes one.

use ogeom_core::{OgeomResult, Tolerances};
use ogeom_math::{Direction, Point, Point2, Transform, Vector, Vector2};

/// How far past its ends a curve's domain still answers: the parametric
/// confusion, or a few ulps of the ends themselves where those are large.
/// A line's unbounded domain ends at a billion, where one ulp is a tenth of
/// a micron, and a step to the end computed as `a + (b - a) * t` lands an
/// ulp or so past it.
pub(crate) fn domain_slack(a: f64, b: f64, tol: Tolerances) -> f64 {
    tol.parametric()
        .max(a.abs().max(b.abs()) * 4.0 * f64::EPSILON)
}

/// Whether a curve evaluated outside its domain says where it was asked
/// from. Read once: the refusal is routine inside booleans, and reading the
/// environment takes a lock.
static DEBUG_DOMAIN: std::sync::LazyLock<bool> =
    std::sync::LazyLock::new(|| std::env::var_os("OGEOM_DEBUG_DOMAIN").is_some());

/// How smooth a curve or surface is.
///
/// Ordered from least to most smooth, so `>=` is a meaningful test.
#[derive(Debug, Clone, Copy, PartialEq, Eq, PartialOrd, Ord)]
pub enum Continuity {
    /// Positions agree. Tangents may not. A corner.
    C0,
    /// Tangent *directions* agree, but their magnitudes need not. Enough for a
    /// visually smooth join, and the usual requirement for a wire.
    G1,
    /// First derivatives agree exactly.
    C1,
    /// Curvature agrees.
    G2,
    /// Second derivatives agree exactly.
    C2,
    /// Differentiable to any order: an analytic surface away from its
    /// degeneracies.
    CInfinity,
}

/// What kind of curve this is, for analytic fast paths.
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
#[non_exhaustive]
pub enum CurveKind {
    /// A straight line.
    Line,
    /// A circle or circular arc.
    Circle,
    /// An ellipse or elliptical arc.
    Ellipse,
    /// One branch of a hyperbola.
    Hyperbola,
    /// A parabola.
    Parabola,
    /// A polynomial or rational Bézier.
    Bezier,
    /// A polynomial or rational B-spline.
    BSpline,
    /// A helix about an axis.
    Helix,
    /// A restriction of another curve to a sub-interval.
    Trimmed,
    /// A curve offset from another.
    Offset,
    /// A pcurve composed with the surface it is drawn on.
    OnSurface,
    /// An affine-plus-trigonometric chart curve.
    Trig,
}

/// What kind of surface this is, for analytic fast paths.
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
#[non_exhaustive]
pub enum SurfaceKind {
    /// A plane.
    Plane,
    /// A circular cylinder.
    Cylinder,
    /// A circular cone.
    Cone,
    /// A sphere.
    Sphere,
    /// A torus.
    Torus,
    /// A polynomial or rational Bézier patch.
    Bezier,
    /// A polynomial or rational B-spline patch.
    BSpline,
    /// A surface of revolution.
    Revolution,
    /// A surface swept by translating a curve.
    Extrusion,
    /// A restriction of another surface to a sub-rectangle.
    Trimmed,
    /// A surface offset from another.
    Offset,
}

impl SurfaceKind {
    /// Whether this surface is a quadric: a plane, cylinder, cone or sphere.
    ///
    /// Quadric pairs have closed-form intersections, which is worth a great
    /// deal: it is the difference between an exact conic and a marched
    /// approximation of one.
    #[must_use]
    pub const fn is_quadric(self) -> bool {
        matches!(
            self,
            Self::Plane | Self::Cylinder | Self::Cone | Self::Sphere
        )
    }

    /// Whether this surface has an exact analytic form, as opposed to being
    /// defined by control points.
    #[must_use]
    pub const fn is_analytic(self) -> bool {
        matches!(
            self,
            Self::Plane | Self::Cylinder | Self::Cone | Self::Sphere | Self::Torus
        )
    }
}

/// A parametric curve in space.
///
/// Implementors must guarantee:
///
/// - [`Curve3d::domain`] returns `(a, b)` with `a < b`, both finite;
/// - [`Curve3d::point_at`] and the derivative methods agree: `d1_at` is the
///   derivative of `point_at`, and so on;
/// - a periodic curve's period is exactly its domain width.
pub trait Curve3d {
    /// The parameter interval over which the curve is defined.
    fn domain(&self) -> (f64, f64);

    /// The point at `u`.
    ///
    /// # Errors
    ///
    /// [`OgeomError::Domain`](ogeom_core::OgeomError::Domain) if `u` is outside the
    /// domain and the curve is not periodic.
    fn point_at(&self, u: f64, tol: Tolerances) -> OgeomResult<Point>;

    /// The first derivative at `u`.
    ///
    /// # Errors
    ///
    /// As [`Curve3d::point_at`].
    fn d1_at(&self, u: f64, tol: Tolerances) -> OgeomResult<Vector>;

    /// Derivatives up to order `n`, with `result[0]` the point itself.
    ///
    /// # Errors
    ///
    /// As [`Curve3d::point_at`].
    fn derivatives_at(&self, u: f64, n: usize, tol: Tolerances) -> OgeomResult<Vec<Vector>>;

    /// What kind of curve this is.
    fn kind(&self) -> CurveKind;

    /// How smooth the curve is across its domain.
    fn continuity(&self) -> Continuity;

    /// Whether the curve's ends meet.
    fn is_closed(&self, tol: Tolerances) -> bool;

    /// Whether the curve continues past its domain by repeating.
    ///
    /// Distinct from being closed: a full circle is both, a closed B-spline
    /// whose ends merely coincide is closed but not periodic, and evaluating
    /// the latter outside its domain is an error.
    fn is_periodic(&self) -> bool;

    /// The unit tangent at `u`.
    ///
    /// # Errors
    ///
    /// [`OgeomError::Construction`](ogeom_core::OgeomError::Construction) at a cusp,
    /// where the derivative vanishes.
    fn tangent_at(&self, u: f64, tol: Tolerances) -> OgeomResult<Direction> {
        Direction::new(self.d1_at(u, tol)?, tol)
    }

    /// The curvature at `u`.
    ///
    /// # Errors
    ///
    /// As [`Curve3d::point_at`].
    fn curvature_at(&self, u: f64, tol: Tolerances) -> OgeomResult<f64> {
        let d = self.derivatives_at(u, 2, tol)?;
        let speed = d[1].magnitude();
        if speed == 0.0 {
            return Ok(0.0);
        }
        Ok(d[1].cross(d[2]).magnitude() / (speed * speed * speed))
    }

    /// The start point.
    ///
    /// # Errors
    ///
    /// As [`Curve3d::point_at`].
    fn start(&self, tol: Tolerances) -> OgeomResult<Point> {
        self.point_at(self.domain().0, tol)
    }

    /// The end point.
    ///
    /// # Errors
    ///
    /// As [`Curve3d::point_at`].
    fn end(&self, tol: Tolerances) -> OgeomResult<Point> {
        self.point_at(self.domain().1, tol)
    }

    /// Bring `u` into the domain, wrapping if the curve is periodic.
    ///
    /// # Errors
    ///
    /// [`OgeomError::Domain`](ogeom_core::OgeomError::Domain) if `u` is outside the
    /// domain of a non-periodic curve by more than `tol.parametric()`.
    fn normalize_parameter(&self, u: f64, tol: Tolerances) -> OgeomResult<f64> {
        let (a, b) = self.domain();
        if self.is_periodic() && u.is_finite() {
            if u >= a && u < b {
                return Ok(u);
            }
            let period = b - a;
            return Ok(a + (u - a).rem_euclid(period));
        }
        let slack = domain_slack(a, b, tol);
        if !u.is_finite() || u < a - slack || u > b + slack {
            if *DEBUG_DOMAIN {
                eprintln!(
                    "DOMAIN {u} outside [{a}, {b}]:\n{}",
                    std::backtrace::Backtrace::force_capture()
                );
            }
            return Err(ogeom_core::ogeom_err!(
                Domain,
                "parameter {u} outside curve domain [{a}, {b}]"
            ));
        }
        Ok(u.clamp(a, b))
    }
}

/// A parametric curve in the plane.
///
/// The same contract as [`Curve3d`], one dimension down. Kept as a separate
/// trait rather than a generic parameter because pcurves (curves in a
/// surface's parameter space) are used differently enough from spatial curves
/// that conflating them invites mistakes.
pub trait Curve2d {
    /// The parameter interval over which the curve is defined.
    fn domain(&self) -> (f64, f64);

    /// The point at `u`.
    ///
    /// # Errors
    ///
    /// [`OgeomError::Domain`](ogeom_core::OgeomError::Domain) if `u` is outside the
    /// domain and the curve is not periodic.
    fn point_at(&self, u: f64, tol: Tolerances) -> OgeomResult<Point2>;

    /// The first derivative at `u`.
    ///
    /// # Errors
    ///
    /// As [`Curve2d::point_at`].
    fn d1_at(&self, u: f64, tol: Tolerances) -> OgeomResult<Vector2>;

    /// Derivatives up to order `n`, with `result[0]` the point itself.
    ///
    /// # Errors
    ///
    /// As [`Curve2d::point_at`].
    fn derivatives_at(&self, u: f64, n: usize, tol: Tolerances) -> OgeomResult<Vec<Vector2>>;

    /// What kind of curve this is.
    fn kind(&self) -> CurveKind;

    /// Whether the curve's ends meet.
    fn is_closed(&self, tol: Tolerances) -> bool;

    /// Whether the curve continues past its domain by repeating.
    fn is_periodic(&self) -> bool;

    /// The start point.
    ///
    /// # Errors
    ///
    /// As [`Curve2d::point_at`].
    fn start(&self, tol: Tolerances) -> OgeomResult<Point2> {
        self.point_at(self.domain().0, tol)
    }

    /// The end point.
    ///
    /// # Errors
    ///
    /// As [`Curve2d::point_at`].
    fn end(&self, tol: Tolerances) -> OgeomResult<Point2> {
        self.point_at(self.domain().1, tol)
    }
}

/// A surface's point and its derivatives through second order, at one place.
///
/// What a foot-point solve needs in one go: see [`Surface::jet_at`].
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct SurfaceJet {
    /// The point at the parameters asked for.
    pub point: Point,
    /// The first derivative along `u`.
    pub du: Vector,
    /// The first derivative along `v`.
    pub dv: Vector,
    /// The second derivative along `u`.
    pub d2u: Vector,
    /// The mixed second derivative.
    pub duv: Vector,
    /// The second derivative along `v`.
    pub d2v: Vector,
}

/// A surface's curvature at one place: the principal curvatures and the
/// directions they are taken in, with the mean and Gaussian curvatures
/// they combine to.
///
/// What curvature display, zebra analysis and a fillet's seat all ask.
/// The principal directions are unit tangents in space, perpendicular to
/// each other. At an umbilic (where every direction curves the same, as
/// everywhere on a sphere or a plane) they are any perpendicular pair
/// in the tangent plane, and [`SurfaceCurvature::is_umbilic`] says so.
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct SurfaceCurvature {
    /// The largest normal curvature.
    pub max: f64,
    /// The smallest normal curvature.
    pub min: f64,
    /// The tangent direction of the largest.
    pub max_direction: Direction,
    /// The tangent direction of the smallest.
    pub min_direction: Direction,
    /// The unit normal the curvatures are signed against: positive where
    /// the surface curves towards it.
    pub normal: Direction,
}

impl SurfaceCurvature {
    /// The mean curvature, `(max + min) / 2`.
    #[must_use]
    pub fn mean(&self) -> f64 {
        f64::midpoint(self.max, self.min)
    }

    /// The Gaussian curvature, `max · min`.
    #[must_use]
    pub fn gaussian(&self) -> f64 {
        self.max * self.min
    }

    /// Whether the two principal curvatures are equal to within `tol`'s
    /// angular resolution of their size, so no direction is principal
    /// above another.
    #[must_use]
    pub fn is_umbilic(&self, tol: Tolerances) -> bool {
        (self.max - self.min).abs() <= tol.angular() * self.max.abs().max(self.min.abs()).max(1.0)
    }

    /// The normal curvature along the tangent direction `along`, signed as
    /// the principal ones are.
    ///
    /// Euler's formula: `max cos²θ + min sin²θ`, with `θ` the angle from the
    /// largest's direction. The part of `along` that leaves the tangent plane
    /// is dropped first. `None` if nothing is left of it, as for a vector
    /// along the normal or a zero one.
    #[must_use]
    pub fn normal_curvature(&self, along: Vector) -> Option<f64> {
        let c = along.dot(self.max_direction.vector());
        let s = along.dot(self.min_direction.vector());
        let size = c.mul_add(c, s * s);
        // Relative to the vector's own length, so a short tangent is still
        // a tangent and only one lying along the normal is refused.
        let length = along.dot(along);
        if length == 0.0 || !size.is_finite() || size <= 1e-24 * length {
            return None;
        }
        Some(self.max.mul_add(c * c, self.min * s * s) / size)
    }
}

/// A parametric surface.
///
/// Implementors must guarantee:
///
/// - both domain intervals are non-empty and finite;
/// - the derivative methods agree with [`Surface::point_at`];
/// - `du` and `dv` are the derivatives along the first and second parameters
///   respectively, in that order, for every surface.
pub trait Surface {
    /// The `u` and `v` parameter intervals.
    fn domain(&self) -> ((f64, f64), (f64, f64));

    /// The point at `(u, v)`.
    ///
    /// # Errors
    ///
    /// [`OgeomError::Domain`](ogeom_core::OgeomError::Domain) if the parameters lie
    /// outside the domain in a direction that is not periodic.
    fn point_at(&self, u: f64, v: f64, tol: Tolerances) -> OgeomResult<Point>;

    /// The two first derivatives at `(u, v)`.
    ///
    /// # Errors
    ///
    /// As [`Surface::point_at`].
    fn d1_at(&self, u: f64, v: f64, tol: Tolerances) -> OgeomResult<(Vector, Vector)>;

    /// The three second derivatives at `(u, v)`: `d2u`, `duv`, `d2v`.
    ///
    /// # Errors
    ///
    /// As [`Surface::point_at`].
    fn d2_at(&self, u: f64, v: f64, tol: Tolerances) -> OgeomResult<(Vector, Vector, Vector)>;

    /// The point and first derivatives together: `(point, du, dv)`.
    ///
    /// For a caller that needs both at every sample (an integral over the
    /// surface): a patch has the point in the same derivative table, and
    /// asking for it separately sums the control grid twice. Agreement with
    /// the separate accessors is to rounding, as for [`Surface::jet_at`].
    ///
    /// # Errors
    ///
    /// As [`Surface::point_at`].
    fn point_d1_at(&self, u: f64, v: f64, tol: Tolerances) -> OgeomResult<(Point, Vector, Vector)> {
        let point = self.point_at(u, v, tol)?;
        let (du, dv) = self.d1_at(u, v, tol)?;
        Ok((point, du, dv))
    }

    /// The point and every derivative through second order, together.
    ///
    /// The default asks the three accessors above, which is right for a
    /// surface carrying its answer in closed form: a plane or a cylinder costs
    /// the same either way.
    ///
    /// It is not right for a tensor-product patch, where each accessor
    /// re-locates the knot spans, rebuilds the basis functions and sums the
    /// control grid again, three times, for numbers the order-two table
    /// already holds. Anything walking a surface and asking for all six at a
    /// point, as a foot-point solve does at every step, pays that over and
    /// over. Such surfaces override this.
    ///
    /// **The contract is agreement to rounding, not to the bit.** A patch sums
    /// its point by de Boor and its derivatives by basis functions, and those
    /// reassociate differently. The combined answer may differ from the
    /// separate accessors' in the last ulp. Callers needing one consistent
    /// jet (every value from the same evaluation) should use this and not
    /// mix it with the accessors at the same parameters.
    ///
    /// # Errors
    ///
    /// As [`Surface::point_at`].
    fn jet_at(&self, u: f64, v: f64, tol: Tolerances) -> OgeomResult<SurfaceJet> {
        let point = self.point_at(u, v, tol)?;
        let (du, dv) = self.d1_at(u, v, tol)?;
        let (d2u, duv, d2v) = self.d2_at(u, v, tol)?;
        Ok(SurfaceJet {
            point,
            du,
            dv,
            d2u,
            duv,
            d2v,
        })
    }

    /// What kind of surface this is.
    fn kind(&self) -> SurfaceKind;

    /// How smooth the surface is across its domain.
    fn continuity(&self) -> Continuity;

    /// Whether the surface closes on itself along `u`.
    fn is_closed_u(&self, tol: Tolerances) -> bool;

    /// Whether the surface closes on itself along `v`.
    fn is_closed_v(&self, tol: Tolerances) -> bool;

    /// Whether `u` repeats past the domain.
    fn is_periodic_u(&self) -> bool;

    /// Whether `v` repeats past the domain.
    fn is_periodic_v(&self) -> bool;

    /// The unit normal at `(u, v)`, following `du x dv`.
    ///
    /// # Errors
    ///
    /// [`OgeomError::Construction`](ogeom_core::OgeomError::Construction) at a
    /// degeneracy (a pole or an apex) where the tangents determine no normal.
    fn normal_at(&self, u: f64, v: f64, tol: Tolerances) -> OgeomResult<Direction> {
        let (du, dv) = self.d1_at(u, v, tol)?;
        // Compared against the square of the larger tangent: at a pole one
        // tangent vanishes, and a test relative to the *product* would find the
        // cross product large by comparison and call the point healthy.
        let scale = du.magnitude().max(dv.magnitude());
        if du.cross(dv).magnitude() <= tol.angular() * scale * scale {
            return Err(ogeom_core::ogeom_err!(
                Construction,
                "surface is degenerate at ({u}, {v}); no normal is determined"
            ));
        }
        Direction::new(du.cross(dv), tol)
    }

    /// The principal curvatures and directions at `(u, v)`.
    ///
    /// From the two fundamental forms of the jet: the principal curvatures
    /// are the eigenvalues of the shape operator `I⁻¹ II`, the directions
    /// its eigenvectors carried into space. Signed against the surface's
    /// own normal, [`Surface::normal_at`], so a sphere of radius `r` seen
    /// from outside reads `-1/r` twice and a cylinder `-1/r` round and `0`
    /// along.
    ///
    /// # Errors
    ///
    /// As [`Surface::normal_at`]: at a degenerate point there is no normal
    /// to sign against.
    fn curvature_at(&self, u: f64, v: f64, tol: Tolerances) -> OgeomResult<SurfaceCurvature> {
        let jet = self.jet_at(u, v, tol)?;
        let normal = self.normal_at(u, v, tol)?;
        let n = normal.vector();
        let (e, f, g) = (jet.du.dot(jet.du), jet.du.dot(jet.dv), jet.dv.dot(jet.dv));
        let (l, m, nn) = (jet.d2u.dot(n), jet.duv.dot(n), jet.d2v.dot(n));
        let det = e * g - f * f;
        let gaussian = (l * nn - m * m) / det;
        let mean = (e * nn - 2.0 * f * m + g * l) / (2.0 * det);
        let spread = (mean * mean - gaussian).max(0.0).sqrt();
        let (max, min) = (mean + spread, mean - spread);
        // A principal direction `a du + b dv` solves `(II - k I)(a, b) = 0`;
        // either row of that matrix gives it, and the row with the larger
        // entries is the one to trust. At an umbilic both rows vanish and
        // any tangent will do: `du` and the normal turned against it.
        let direction_for = |k: f64| -> OgeomResult<Vector> {
            let row1 = (l - k * e, m - k * f);
            let row2 = (m - k * f, nn - k * g);
            let (a, b) = if row1.0.hypot(row1.1) >= row2.0.hypot(row2.1) {
                (row1.1, -row1.0)
            } else {
                (row2.1, -row2.0)
            };
            let along = jet.du * a + jet.dv * b;
            let scale = jet.du.magnitude().max(jet.dv.magnitude());
            if along.magnitude() <= tol.angular() * scale * (a.abs() + b.abs()).max(1.0)
                || a.abs() + b.abs() <= tol.angular() * (l.abs() + m.abs() + nn.abs()).max(1.0)
            {
                return Ok(jet.du);
            }
            Ok(along)
        };
        let max_direction = Direction::new(direction_for(max)?, tol)?;
        let min_direction = if spread <= tol.angular() * max.abs().max(min.abs()).max(1.0) {
            Direction::new(n.cross(max_direction.vector()), tol)?
        } else {
            Direction::new(direction_for(min)?, tol)?
        };
        Ok(SurfaceCurvature {
            max,
            min,
            max_direction,
            min_direction,
            normal,
        })
    }

    /// Whether the surface degenerates at `(u, v)`.
    ///
    /// # Errors
    ///
    /// As [`Surface::point_at`].
    fn is_degenerate_at(&self, u: f64, v: f64, tol: Tolerances) -> OgeomResult<bool> {
        let (du, dv) = self.d1_at(u, v, tol)?;
        let scale = du.magnitude().max(dv.magnitude());
        Ok(du.cross(dv).magnitude() <= tol.angular() * scale * scale)
    }

    /// Bring `(u, v)` into the domain, wrapping in whichever directions are
    /// periodic, or closed.
    ///
    /// A surface that merely closes on itself (a clamped B-spline tube
    /// whose first and last control columns coincide) has the same points
    /// at both ends of its domain exactly as a periodic one does, and a
    /// parameter a whole period past the end names a point it has. A face
    /// whose trim runs right round such a tube has a ring that straddles
    /// the join whichever way it is slid, so somewhere it is asked past the
    /// end, and a refusal there leaves the face unmeshable. Closure is
    /// consulted only once a parameter is actually outside, so the common
    /// case pays nothing for it.
    ///
    /// # Errors
    ///
    /// [`OgeomError::Domain`](ogeom_core::OgeomError::Domain) if a parameter is outside
    /// a direction's domain by more than `tol.parametric()` and that
    /// direction neither repeats nor closes.
    fn normalize_parameters(&self, u: f64, v: f64, tol: Tolerances) -> OgeomResult<(f64, f64)> {
        // Every evaluation of every surface passes through here, so the
        // in-range path is two comparisons, and everything rarer (wrapping
        // a periodic direction, wrapping a closed one, refusing) lives out
        // of line. An indirect call for closure on the common path costs a
        // measurable share of a large mesh's time, for nothing in-range.
        let ((ua, ub), (va, vb)) = self.domain();
        let slack = tol.parametric();
        let u = if u >= ua - slack && u <= ub + slack && !self.is_periodic_u() {
            u.clamp(ua, ub)
        } else if u >= ua && u < ub {
            // Periodic and already in its first period: nothing to wrap.
            u
        } else {
            self.parameter_outside(u, ua, ub, self.is_periodic_u(), true, tol)?
        };
        let v = if v >= va - slack && v <= vb + slack && !self.is_periodic_v() {
            v.clamp(va, vb)
        } else if v >= va && v < vb {
            v
        } else {
            self.parameter_outside(v, va, vb, self.is_periodic_v(), false, tol)?
        };
        Ok((u, v))
    }

    /// A parameter outside its domain, or on a periodic direction: wrapped
    /// where the direction repeats or closes on itself, refused otherwise.
    ///
    /// A surface that merely closes on itself has the same points at both
    /// ends of its domain exactly as a periodic one does, so a parameter a
    /// whole period past the end names a point it has and is wrapped to it,
    /// for the reason given on `normalize_parameters`. Closure is consulted
    /// only here, once a parameter is actually outside.
    ///
    /// # Errors
    ///
    /// [`OgeomError::Domain`](ogeom_core::OgeomError::Domain) if the parameter is
    /// outside by more than `tol.parametric()` and the direction neither
    /// repeats nor closes.
    #[cold]
    fn parameter_outside(
        &self,
        t: f64,
        a: f64,
        b: f64,
        periodic: bool,
        across: bool,
        tol: Tolerances,
    ) -> OgeomResult<f64> {
        if periodic && t.is_finite() {
            return Ok(a + (t - a).rem_euclid(b - a));
        }
        let closed = if across {
            self.is_closed_u(tol)
        } else {
            self.is_closed_v(tol)
        };
        if t.is_finite() && b > a && closed {
            return Ok(a + (t - a).rem_euclid(b - a));
        }
        Err(ogeom_core::ogeom_err!(
            Domain,
            "{} parameter {t} outside [{a}, {b}]",
            if across { "u" } else { "v" }
        ))
    }
}

/// Geometry that can be moved by a similarity.
///
/// Separate from the evaluation traits because a *view* of geometry (a
/// borrowed adaptor over someone else's data) can be evaluated but not moved.
pub trait Transformable: Sized {
    /// This geometry moved by `t`.
    ///
    /// # Errors
    ///
    /// [`OgeomError::Construction`](ogeom_core::OgeomError::Construction) if the
    /// transformed geometry would be degenerate.
    fn transformed(&self, t: &Transform, tol: Tolerances) -> OgeomResult<Self>;
}

/// Reversible geometry: the same point set, traversed the other way.
pub trait Reversible: Sized {
    /// This geometry with its parameter direction reversed.
    ///
    /// The domain is preserved, so a curve reversed still runs over the same
    /// interval. Only the direction of travel changes. Preserving the domain
    /// matters because trimming ranges elsewhere refer to it.
    #[must_use]
    fn reversed(&self) -> Self;
}

#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn continuity_is_ordered_from_least_to_most_smooth() {
        assert!(Continuity::C0 < Continuity::G1);
        assert!(Continuity::G1 < Continuity::C1);
        assert!(Continuity::C1 < Continuity::G2);
        assert!(Continuity::G2 < Continuity::C2);
        assert!(Continuity::C2 < Continuity::CInfinity);
        // The ordering is what makes a requirement expressible as a comparison.
        assert!(Continuity::CInfinity >= Continuity::C1);
    }

    #[test]
    fn quadrics_and_analytics_are_classified_correctly() {
        for k in [
            SurfaceKind::Plane,
            SurfaceKind::Cylinder,
            SurfaceKind::Cone,
            SurfaceKind::Sphere,
        ] {
            assert!(k.is_quadric(), "{k:?}");
            assert!(k.is_analytic(), "{k:?}");
        }
        // A torus is analytic but quartic, not quadric, a distinction that
        // matters, since quadric pairs have closed-form intersections and
        // torus pairs do not.
        assert!(!SurfaceKind::Torus.is_quadric());
        assert!(SurfaceKind::Torus.is_analytic());

        for k in [
            SurfaceKind::BSpline,
            SurfaceKind::Bezier,
            SurfaceKind::Revolution,
            SurfaceKind::Extrusion,
            SurfaceKind::Trimmed,
            SurfaceKind::Offset,
        ] {
            assert!(!k.is_quadric(), "{k:?}");
            assert!(!k.is_analytic(), "{k:?}");
        }
    }

    /// A line's unbounded domain ends at a billion, where an ulp is a tenth
    /// of a micron: a parameter an ulp or two past the end, which a step
    /// computed toward the end lands on, is the end. A short segment's end
    /// keeps the parametric slack, and a micron past either is refused.
    #[test]
    #[allow(clippy::unwrap_used)]
    fn a_domain_end_answers_an_ulp_past_it_at_any_magnitude() {
        let tol = Tolerances::millimetres();
        let axis = ogeom_math::Axis::new(Point::ORIGIN, Direction::X);
        let line = crate::LineCurve::new(axis);
        let (a, b) = line.domain();
        for u in [b, b.next_up(), b.next_up().next_up(), a.next_down()] {
            assert!(line.point_at(u, tol).is_ok(), "{u} on [{a}, {b}]");
        }
        assert!(line.point_at(b + 1e-3, tol).is_err());
        let segment =
            crate::LineCurve::segment(Point::ORIGIN, Point::new(1.0, 0.0, 0.0), tol).unwrap();
        assert!(segment.point_at(1.0 + 1e-10, tol).is_ok());
        assert!(segment.point_at(1.0 + 1e-6, tol).is_err());
        let flat = crate::Line2d::over(
            ogeom_math::Axis2::new(Point2::ORIGIN, ogeom_math::Direction2::X),
            -crate::curve::LINE_EXTENT,
            crate::curve::LINE_EXTENT,
        )
        .unwrap();
        assert!(Curve2d::point_at(&flat, crate::curve::LINE_EXTENT.next_up(), tol).is_ok());
        assert!(Curve2d::point_at(&flat, crate::curve::LINE_EXTENT + 1e-3, tol).is_err());
    }
}