smlib: surface/surface intersection wrapper (smlib_surface_*, smlib_ssi)

ffi/ssi.rs declares the C API's SSI block; ssi.rs wraps it: Surface
(plane, cylinder, cone, sphere, torus, NURBS, with the SSI corpus's
fields, angles in degrees), Method (Global, Advanced, Basic) and
intersect(), which returns every curve's type, UV flags, deviation,
interval and samples (xyz, uv on each surface, each surface at its uv)
and the isolated points. Without the native feature every call returns
NotBuilt.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

crates/smlib/src/ffi.rs

@@ -184,6 +184,7 @@ pub mod flat;
 pub mod history;
 pub mod measure;
 pub mod ops;
+pub mod ssi;
 
 /// Where [`set_log`] sends SMLib's log lines.
 #[cfg(feature = "native")]

crates/smlib/src/ffi/ssi.rs

@@ -0,0 +1,361 @@
+//! The surface/surface intersection block of `capi/include/smlib.h` (`smlib_surface_*`,
+//! `smlib_ssi*`), declared by hand. Part of the FFI boundary, like its parent: called only
+//! through [`crate::ssi`].
+#![allow(unsafe_code)]
+
+/// `smlib_surface`, opaque.
+#[repr(C)]
+pub struct RawSurface {
+    _private: [u8; 0],
+}
+
+/// `smlib_ssi_result`, opaque.
+#[repr(C)]
+pub struct RawResult {
+    _private: [u8; 0],
+}
+
+#[cfg(feature = "native")]
+unsafe extern "C" {
+    fn smlib_surface_plane(origin: *const f64, xaxis: *const f64, yaxis: *const f64, u: *const f64, v: *const f64, out: *mut *mut RawSurface) -> i32;
+    fn smlib_surface_cylinder(
+        origin: *const f64,
+        xaxis: *const f64,
+        yaxis: *const f64,
+        radius: f64,
+        angle_deg: *const f64,
+        height: *const f64,
+        out: *mut *mut RawSurface,
+    ) -> i32;
+    fn smlib_surface_cone(
+        apex: *const f64,
+        xaxis: *const f64,
+        yaxis: *const f64,
+        half_angle_deg: f64,
+        height: *const f64,
+        out: *mut *mut RawSurface,
+    ) -> i32;
+    fn smlib_surface_sphere(
+        center: *const f64,
+        xaxis: *const f64,
+        yaxis: *const f64,
+        radius: f64,
+        lon_deg: *const f64,
+        lat_deg: *const f64,
+        out: *mut *mut RawSurface,
+    ) -> i32;
+    fn smlib_surface_torus(
+        center: *const f64,
+        xaxis: *const f64,
+        yaxis: *const f64,
+        major_radius: f64,
+        minor_radius: f64,
+        major_deg: *const f64,
+        minor_deg: *const f64,
+        out: *mut *mut RawSurface,
+    ) -> i32;
+    fn smlib_surface_nurbs(
+        degree_u: i32,
+        degree_v: i32,
+        nu: i32,
+        nv: i32,
+        knots_u: *const f64,
+        n_knots_u: i64,
+        knots_v: *const f64,
+        n_knots_v: i64,
+        points: *const f64,
+        weights: *const f64,
+        out: *mut *mut RawSurface,
+    ) -> i32;
+    fn smlib_surface_free(s: *mut RawSurface);
+    fn smlib_surface_domain(s: *const RawSurface, out: *mut f64) -> i32;
+    fn smlib_surface_evaluate(s: *const RawSurface, u: f64, v: f64, out: *mut f64) -> i32;
+    fn smlib_ssi(a: *const RawSurface, b: *const RawSurface, method: i32, approx_tol: f64, samples_per_curve: i32, out: *mut *mut RawResult) -> i32;
+    fn smlib_ssi_result_free(r: *mut RawResult);
+    fn smlib_ssi_curve_count(r: *const RawResult, out: *mut i32) -> i32;
+    fn smlib_ssi_curve_info(
+        r: *const RawResult,
+        i: i32,
+        ty: *mut i32,
+        flags: *mut i32,
+        samples: *mut i32,
+        deviation: *mut f64,
+        interval: *mut f64,
+    ) -> i32;
+    fn smlib_ssi_curve_samples(
+        r: *const RawResult,
+        i: i32,
+        cap: i64,
+        xyz: *mut f64,
+        uv_a: *mut f64,
+        uv_b: *mut f64,
+        a_at_uv: *mut f64,
+        b_at_uv: *mut f64,
+    ) -> i32;
+    fn smlib_ssi_point_count(r: *const RawResult, out: *mut i32) -> i32;
+    fn smlib_ssi_point(r: *const RawResult, i: i32, xyz: *mut f64, uv_a: *mut f64, uv_b: *mut f64) -> i32;
+}
+
+/// An owned `smlib_surface`, freed once on drop.
+#[cfg_attr(not(feature = "native"), allow(dead_code))]
+pub struct SurfaceHandle(std::ptr::NonNull<RawSurface>);
+
+/// An owned `smlib_ssi_result`, freed once on drop.
+#[cfg_attr(not(feature = "native"), allow(dead_code))]
+pub struct ResultHandle(std::ptr::NonNull<RawResult>);
+
+// Only ever used through the C API, which serialises every call.
+unsafe impl Send for SurfaceHandle {}
+unsafe impl Sync for SurfaceHandle {}
+unsafe impl Send for ResultHandle {}
+unsafe impl Sync for ResultHandle {}
+
+impl Drop for SurfaceHandle {
+    fn drop(&mut self) {
+        #[cfg(feature = "native")]
+        unsafe {
+            smlib_surface_free(self.0.as_ptr())
+        };
+    }
+}
+
+impl Drop for ResultHandle {
+    fn drop(&mut self) {
+        #[cfg(feature = "native")]
+        unsafe {
+            smlib_ssi_result_free(self.0.as_ptr())
+        };
+    }
+}
+
+/// One curve's `smlib_ssi_curve_info`.
+#[derive(Clone, Copy, Debug)]
+pub struct CurveInfo {
+    pub ty: i32,
+    pub flags: i32,
+    pub samples: i32,
+    pub deviation: f64,
+    pub interval: [f64; 2],
+}
+
+/// One curve's samples, flat: 3 (xyz, a_at_uv, b_at_uv) or 2 (uv_a, uv_b) values per sample.
+#[derive(Clone, Debug, Default)]
+pub struct Samples {
+    pub xyz: Vec<f64>,
+    pub uv_a: Vec<f64>,
+    pub uv_b: Vec<f64>,
+    pub a_at_uv: Vec<f64>,
+    pub b_at_uv: Vec<f64>,
+}
+
+#[cfg(feature = "native")]
+mod calls {
+    use super::*;
+    use crate::ffi::{ERR_ARG, ERR_EMPTY, OK};
+
+    fn new_surface(f: impl FnOnce(*mut *mut RawSurface) -> i32) -> Result<SurfaceHandle, i32> {
+        let mut out: *mut RawSurface = std::ptr::null_mut();
+        match (f(&mut out), std::ptr::NonNull::new(out)) {
+            (OK, Some(p)) => Ok(SurfaceHandle(p)),
+            (OK, None) => Err(ERR_EMPTY),
+            (s, _) => Err(s),
+        }
+    }
+
+    pub fn plane(origin: &[f64; 3], xaxis: &[f64; 3], yaxis: &[f64; 3], u: &[f64; 2], v: &[f64; 2]) -> Result<SurfaceHandle, i32> {
+        new_surface(|out| unsafe { smlib_surface_plane(origin.as_ptr(), xaxis.as_ptr(), yaxis.as_ptr(), u.as_ptr(), v.as_ptr(), out) })
+    }
+    pub fn cylinder(
+        origin: &[f64; 3],
+        xaxis: &[f64; 3],
+        yaxis: &[f64; 3],
+        radius: f64,
+        angle_deg: &[f64; 2],
+        height: &[f64; 2],
+    ) -> Result<SurfaceHandle, i32> {
+        new_surface(|out| unsafe {
+            smlib_surface_cylinder(origin.as_ptr(), xaxis.as_ptr(), yaxis.as_ptr(), radius, angle_deg.as_ptr(), height.as_ptr(), out)
+        })
+    }
+    pub fn cone(apex: &[f64; 3], xaxis: &[f64; 3], yaxis: &[f64; 3], half_angle_deg: f64, height: &[f64; 2]) -> Result<SurfaceHandle, i32> {
+        new_surface(|out| unsafe { smlib_surface_cone(apex.as_ptr(), xaxis.as_ptr(), yaxis.as_ptr(), half_angle_deg, height.as_ptr(), out) })
+    }
+    pub fn sphere(
+        center: &[f64; 3],
+        xaxis: &[f64; 3],
+        yaxis: &[f64; 3],
+        radius: f64,
+        lon_deg: &[f64; 2],
+        lat_deg: &[f64; 2],
+    ) -> Result<SurfaceHandle, i32> {
+        new_surface(|out| unsafe {
+            smlib_surface_sphere(center.as_ptr(), xaxis.as_ptr(), yaxis.as_ptr(), radius, lon_deg.as_ptr(), lat_deg.as_ptr(), out)
+        })
+    }
+    pub fn torus(
+        center: &[f64; 3],
+        xaxis: &[f64; 3],
+        yaxis: &[f64; 3],
+        major_radius: f64,
+        minor_radius: f64,
+        major_deg: &[f64; 2],
+        minor_deg: &[f64; 2],
+    ) -> Result<SurfaceHandle, i32> {
+        new_surface(|out| unsafe {
+            smlib_surface_torus(
+                center.as_ptr(),
+                xaxis.as_ptr(),
+                yaxis.as_ptr(),
+                major_radius,
+                minor_radius,
+                major_deg.as_ptr(),
+                minor_deg.as_ptr(),
+                out,
+            )
+        })
+    }
+    /// `points` holds `nu * nv` points (3 values each), `weights` `nu * nv` or none; the C side
+    /// checks every count against the degrees and sizes.
+    pub fn nurbs(
+        degree: [i32; 2],
+        nu: i32,
+        nv: i32,
+        knots_u: &[f64],
+        knots_v: &[f64],
+        points: &[f64],
+        weights: Option<&[f64]>,
+    ) -> Result<SurfaceHandle, i32> {
+        let n = i64::from(nu).checked_mul(i64::from(nv)).ok_or(ERR_ARG)?;
+        let fits = |len: usize, per: i64| n.checked_mul(per).is_some_and(|need| i64::try_from(len).is_ok_and(|l| l == need));
+        if !fits(points.len(), 3) || weights.is_some_and(|w| !fits(w.len(), 1)) {
+            return Err(ERR_ARG);
+        }
+        let ku = i64::try_from(knots_u.len()).map_err(|_| ERR_ARG)?;
+        let kv = i64::try_from(knots_v.len()).map_err(|_| ERR_ARG)?;
+        let w = weights.map_or(std::ptr::null(), <[f64]>::as_ptr);
+        new_surface(|out| unsafe {
+            smlib_surface_nurbs(degree[0], degree[1], nu, nv, knots_u.as_ptr(), ku, knots_v.as_ptr(), kv, points.as_ptr(), w, out)
+        })
+    }
+    pub fn domain(s: &SurfaceHandle) -> Result<[f64; 4], i32> {
+        let mut d = [0.0; 4];
+        match unsafe { smlib_surface_domain(s.0.as_ptr(), d.as_mut_ptr()) } {
+            OK => Ok(d),
+            st => Err(st),
+        }
+    }
+    pub fn evaluate(s: &SurfaceHandle, u: f64, v: f64) -> Result<[f64; 3], i32> {
+        let mut p = [0.0; 3];
+        match unsafe { smlib_surface_evaluate(s.0.as_ptr(), u, v, p.as_mut_ptr()) } {
+            OK => Ok(p),
+            st => Err(st),
+        }
+    }
+    pub fn ssi(a: &SurfaceHandle, b: &SurfaceHandle, method: i32, approx_tol: f64, samples: i32) -> Result<ResultHandle, i32> {
+        let mut out: *mut RawResult = std::ptr::null_mut();
+        match (unsafe { smlib_ssi(a.0.as_ptr(), b.0.as_ptr(), method, approx_tol, samples, &mut out) }, std::ptr::NonNull::new(out)) {
+            (OK, Some(p)) => Ok(ResultHandle(p)),
+            (OK, None) => Err(ERR_EMPTY),
+            (s, _) => Err(s),
+        }
+    }
+    pub fn curve_count(r: &ResultHandle) -> Result<i32, i32> {
+        let mut n = 0;
+        match unsafe { smlib_ssi_curve_count(r.0.as_ptr(), &mut n) } {
+            OK => Ok(n),
+            s => Err(s),
+        }
+    }
+    pub fn curve_info(r: &ResultHandle, i: i32) -> Result<CurveInfo, i32> {
+        let mut c = CurveInfo { ty: 0, flags: 0, samples: 0, deviation: 0.0, interval: [0.0; 2] };
+        match unsafe { smlib_ssi_curve_info(r.0.as_ptr(), i, &mut c.ty, &mut c.flags, &mut c.samples, &mut c.deviation, c.interval.as_mut_ptr()) } {
+            OK => Ok(c),
+            s => Err(s),
+        }
+    }
+    pub fn curve_samples(r: &ResultHandle, i: i32, n: usize) -> Result<Samples, i32> {
+        let cap = i64::try_from(n).map_err(|_| ERR_ARG)?;
+        let three = n.checked_mul(3).ok_or(ERR_ARG)?;
+        let mut s =
+            Samples { xyz: vec![0.0; three], uv_a: vec![0.0; 2 * n], uv_b: vec![0.0; 2 * n], a_at_uv: vec![0.0; three], b_at_uv: vec![0.0; three] };
+        let st = unsafe {
+            smlib_ssi_curve_samples(
+                r.0.as_ptr(),
+                i,
+                cap,
+                s.xyz.as_mut_ptr(),
+                s.uv_a.as_mut_ptr(),
+                s.uv_b.as_mut_ptr(),
+                s.a_at_uv.as_mut_ptr(),
+                s.b_at_uv.as_mut_ptr(),
+            )
+        };
+        if st == OK { Ok(s) } else { Err(st) }
+    }
+    pub fn point_count(r: &ResultHandle) -> Result<i32, i32> {
+        let mut n = 0;
+        match unsafe { smlib_ssi_point_count(r.0.as_ptr(), &mut n) } {
+            OK => Ok(n),
+            s => Err(s),
+        }
+    }
+    pub fn point(r: &ResultHandle, i: i32) -> Result<([f64; 3], [f64; 2], [f64; 2]), i32> {
+        let (mut p, mut a, mut b) = ([0.0; 3], [0.0; 2], [0.0; 2]);
+        match unsafe { smlib_ssi_point(r.0.as_ptr(), i, p.as_mut_ptr(), a.as_mut_ptr(), b.as_mut_ptr()) } {
+            OK => Ok((p, a, b)),
+            s => Err(s),
+        }
+    }
+}
+
+/// The same calls when SMLib is not built: each one fails with `NOT_BUILT`.
+#[cfg(not(feature = "native"))]
+mod calls {
+    use super::*;
+    use crate::ffi::NOT_BUILT;
+
+    pub fn plane(_: &[f64; 3], _: &[f64; 3], _: &[f64; 3], _: &[f64; 2], _: &[f64; 2]) -> Result<SurfaceHandle, i32> {
+        Err(NOT_BUILT)
+    }
+    pub fn cylinder(_: &[f64; 3], _: &[f64; 3], _: &[f64; 3], _: f64, _: &[f64; 2], _: &[f64; 2]) -> Result<SurfaceHandle, i32> {
+        Err(NOT_BUILT)
+    }
+    pub fn cone(_: &[f64; 3], _: &[f64; 3], _: &[f64; 3], _: f64, _: &[f64; 2]) -> Result<SurfaceHandle, i32> {
+        Err(NOT_BUILT)
+    }
+    pub fn sphere(_: &[f64; 3], _: &[f64; 3], _: &[f64; 3], _: f64, _: &[f64; 2], _: &[f64; 2]) -> Result<SurfaceHandle, i32> {
+        Err(NOT_BUILT)
+    }
+    pub fn torus(_: &[f64; 3], _: &[f64; 3], _: &[f64; 3], _: f64, _: f64, _: &[f64; 2], _: &[f64; 2]) -> Result<SurfaceHandle, i32> {
+        Err(NOT_BUILT)
+    }
+    pub fn nurbs(_: [i32; 2], _: i32, _: i32, _: &[f64], _: &[f64], _: &[f64], _: Option<&[f64]>) -> Result<SurfaceHandle, i32> {
+        Err(NOT_BUILT)
+    }
+    pub fn domain(_: &SurfaceHandle) -> Result<[f64; 4], i32> {
+        Err(NOT_BUILT)
+    }
+    pub fn evaluate(_: &SurfaceHandle, _: f64, _: f64) -> Result<[f64; 3], i32> {
+        Err(NOT_BUILT)
+    }
+    pub fn ssi(_: &SurfaceHandle, _: &SurfaceHandle, _: i32, _: f64, _: i32) -> Result<ResultHandle, i32> {
+        Err(NOT_BUILT)
+    }
+    pub fn curve_count(_: &ResultHandle) -> Result<i32, i32> {
+        Err(NOT_BUILT)
+    }
+    pub fn curve_info(_: &ResultHandle, _: i32) -> Result<CurveInfo, i32> {
+        Err(NOT_BUILT)
+    }
+    pub fn curve_samples(_: &ResultHandle, _: i32, _: usize) -> Result<Samples, i32> {
+        Err(NOT_BUILT)
+    }
+    pub fn point_count(_: &ResultHandle) -> Result<i32, i32> {
+        Err(NOT_BUILT)
+    }
+    pub fn point(_: &ResultHandle, _: i32) -> Result<([f64; 3], [f64; 2], [f64; 2]), i32> {
+        Err(NOT_BUILT)
+    }
+}
+
+pub use calls::*;

crates/smlib/src/lib.rs

@@ -23,6 +23,7 @@ pub mod history;
 pub mod import;
 mod measure;
 mod ops;
+pub mod ssi;
 
 pub use check::{CheckLevel, CheckReport};
 pub use measure::{MEASURE_ACCURACY, MassProperties};

crates/smlib/src/ssi.rs

@@ -0,0 +1,272 @@
+//! Surface/surface intersection (`capi/include/smlib.h`, "Surface/surface intersection"): two
+//! surfaces described in plain numbers, intersected by one of SMLib's SSI entry points, each
+//! intersection curve sampled uniformly in its parameters.
+//!
+//! The constructors take exactly the fields of an SSI corpus case (`tests/ssi-corpus`), angles in
+//! degrees. Each surface's UV domain is SMLib's natural one ([`Surface::domain`]): (u, v) itself
+//! for planes, the parameterisation of SMLib's internal NURBS form for the other analytic kinds.
+
+use crate::{Error, Result, error, ffi, init};
+
+/// `SMLIB_SSI_MAX_SAMPLES`: most samples per curve.
+pub const MAX_SAMPLES: usize = 100_000;
+
+/// A surface with its UV domain, owned on the C side.
+pub struct Surface(ffi::ssi::SurfaceHandle);
+
+impl std::fmt::Debug for Surface {
+    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
+        f.write_str("smlib::ssi::Surface")
+    }
+}
+
+fn int(n: usize) -> Result<i32> {
+    i32::try_from(n).map_err(|_| Error::Argument)
+}
+
+impl Surface {
+    /// `origin + u*x + v*y` (x, y the axes made unit), for u in `u`, v in `v`.
+    pub fn plane(origin: [f64; 3], xaxis: [f64; 3], yaxis: [f64; 3], u: [f64; 2], v: [f64; 2]) -> Result<Surface> {
+        init()?;
+        ffi::ssi::plane(&origin, &xaxis, &yaxis, &u, &v).map(Surface).map_err(error)
+    }
+
+    /// A cylinder of `radius` about the axis through `origin` along `xaxis × yaxis`, angle
+    /// `angle_deg` from `xaxis` towards `yaxis`, for heights `height` measured from `origin`.
+    pub fn cylinder(origin: [f64; 3], xaxis: [f64; 3], yaxis: [f64; 3], radius: f64, angle_deg: [f64; 2], height: [f64; 2]) -> Result<Surface> {
+        init()?;
+        ffi::ssi::cylinder(&origin, &xaxis, &yaxis, radius, &angle_deg, &height).map(Surface).map_err(error)
+    }
+
+    /// A full cone with its apex at `apex` along `xaxis × yaxis`, radius `h * tan(half_angle_deg)`
+    /// at distance h along the axis, for h in `height` (from 0: through the apex).
+    pub fn cone(apex: [f64; 3], xaxis: [f64; 3], yaxis: [f64; 3], half_angle_deg: f64, height: [f64; 2]) -> Result<Surface> {
+        init()?;
+        ffi::ssi::cone(&apex, &xaxis, &yaxis, half_angle_deg, &height).map(Surface).map_err(error)
+    }
+
+    /// A sphere of `radius` at `center`, over longitudes `lon_deg` and latitudes `lat_deg`.
+    pub fn sphere(center: [f64; 3], xaxis: [f64; 3], yaxis: [f64; 3], radius: f64, lon_deg: [f64; 2], lat_deg: [f64; 2]) -> Result<Surface> {
+        init()?;
+        ffi::ssi::sphere(&center, &xaxis, &yaxis, radius, &lon_deg, &lat_deg).map(Surface).map_err(error)
+    }
+
+    /// A ring torus at `center` about `xaxis × yaxis`, `major_deg` round the axis and `minor_deg`
+    /// round the tube.
+    pub fn torus(
+        center: [f64; 3],
+        xaxis: [f64; 3],
+        yaxis: [f64; 3],
+        major_radius: f64,
+        minor_radius: f64,
+        major_deg: [f64; 2],
+        minor_deg: [f64; 2],
+    ) -> Result<Surface> {
+        init()?;
+        ffi::ssi::torus(&center, &xaxis, &yaxis, major_radius, minor_radius, &major_deg, &minor_deg).map(Surface).map_err(error)
+    }
+
+    /// A (rational) B-spline surface: `nu × nv` control points, point (i, j) (i along u) at
+    /// `points[i * nv + j]`; `weights` likewise, or `None` for all 1.
+    pub fn nurbs(
+        degree: [usize; 2],
+        nu: usize,
+        nv: usize,
+        knots_u: &[f64],
+        knots_v: &[f64],
+        points: &[[f64; 3]],
+        weights: Option<&[f64]>,
+    ) -> Result<Surface> {
+        init()?;
+        let flat: Vec<f64> = points.iter().flatten().copied().collect();
+        ffi::ssi::nurbs([int(degree[0])?, int(degree[1])?], int(nu)?, int(nv)?, knots_u, knots_v, &flat, weights).map(Surface).map_err(error)
+    }
+
+    /// The UV domain the intersection runs within, as `[u0, u1, v0, v1]`.
+    pub fn domain(&self) -> Result<[f64; 4]> {
+        ffi::ssi::domain(&self.0).map_err(error)
+    }
+
+    /// The surface's point at (u, v).
+    pub fn evaluate(&self, u: f64, v: f64) -> Result<[f64; 3]> {
+        ffi::ssi::evaluate(&self.0, u, v).map_err(error)
+    }
+}
+
+/// Which of SMLib's intersectors runs.
+#[derive(Clone, Copy, Debug, PartialEq, Eq)]
+pub enum Method {
+    /// `SmSurface::GlobalSurfaceIntersect`: analytic special cases, else [`Method::Advanced`].
+    Global,
+    /// `SmAdvSurfaceIntersector`, the general NURBS path.
+    Advanced,
+    /// `SmSurfaceIntersector`, the basic intersector.
+    Basic,
+}
+
+impl Method {
+    fn code(self) -> i32 {
+        match self {
+            Method::Global => 0,
+            Method::Advanced => 1,
+            Method::Basic => 2,
+        }
+    }
+}
+
+/// A curve's type, SMLib's `SmTsectCurveType`.
+#[derive(Clone, Copy, Debug, PartialEq, Eq)]
+pub enum CurveType {
+    /// The method reported no type.
+    None,
+    Unknown,
+    /// Degenerate: one point where the surfaces touch.
+    Touching,
+    /// Transversal.
+    Crossing,
+    /// The surfaces are tangent along it.
+    Tangent,
+    /// Where the surfaces butt up against each other.
+    Coincident,
+    /// A small angle of intersection.
+    NearTangent,
+    /// Bounds a region where the surfaces coincide.
+    RegionBoundary,
+    /// A value this crate does not know.
+    Other(i32),
+}
+
+impl CurveType {
+    fn from_code(c: i32) -> CurveType {
+        match c {
+            -1 => CurveType::None,
+            0 => CurveType::Unknown,
+            1 => CurveType::Touching,
+            2 => CurveType::Crossing,
+            3 => CurveType::Tangent,
+            4 => CurveType::Coincident,
+            5 => CurveType::NearTangent,
+            6 => CurveType::RegionBoundary,
+            c => CurveType::Other(c),
+        }
+    }
+}
+
+/// One intersection curve, sampled at the same count of parameters on the 3D curve and on each
+/// UV curve. A value SMLib could not evaluate is NaN, and so is every UV value (with the surface
+/// at it) of a UV curve the method did not return ([`SsiCurve::has_uv_a`]).
+#[derive(Clone, Debug)]
+pub struct SsiCurve {
+    pub curve_type: CurveType,
+    pub has_uv_a: bool,
+    pub has_uv_b: bool,
+    /// SMLib's deviation estimate (NaN: none).
+    pub deviation: f64,
+    /// The 3D curve's natural parameter interval, which `xyz` samples uniformly.
+    pub interval: [f64; 2],
+    pub xyz: Vec<[f64; 3]>,
+    pub uv_a: Vec<[f64; 2]>,
+    pub uv_b: Vec<[f64; 2]>,
+    /// Surface a at `uv_a`, b at `uv_b`.
+    pub a_at_uv: Vec<[f64; 3]>,
+    pub b_at_uv: Vec<[f64; 3]>,
+}
+
+/// An isolated point: a [`CurveType::Touching`] curve's first sample.
+#[derive(Clone, Copy, Debug)]
+pub struct SsiPoint {
+    pub xyz: [f64; 3],
+    pub uv_a: [f64; 2],
+    pub uv_b: [f64; 2],
+}
+
+/// What an intersection found.
+#[derive(Clone, Debug, Default)]
+pub struct SsiResult {
+    pub curves: Vec<SsiCurve>,
+    pub points: Vec<SsiPoint>,
+}
+
+fn triples(v: &[f64]) -> Vec<[f64; 3]> {
+    v.as_chunks::<3>().0.to_vec()
+}
+
+fn pairs(v: &[f64]) -> Vec<[f64; 2]> {
+    v.as_chunks::<2>().0.to_vec()
+}
+
+/// Intersects `a` and `b` within their domains with `method` at approximation tolerance
+/// `approx_tol` (<= 0: SMLib's default, 1/1000 of the surfaces' size), sampling each curve at
+/// `samples` parameters (2..=[`MAX_SAMPLES`]). An intersector failure is SMLib's status.
+pub fn intersect(a: &Surface, b: &Surface, method: Method, approx_tol: f64, samples: usize) -> Result<SsiResult> {
+    if !(2..=MAX_SAMPLES).contains(&samples) {
+        return Err(Error::Argument);
+    }
+    let r = ffi::ssi::ssi(&a.0, &b.0, method.code(), approx_tol, int(samples)?).map_err(error)?;
+    let n = ffi::ssi::curve_count(&r).map_err(error)?;
+    let mut curves = Vec::new();
+    for i in 0..n {
+        let info = ffi::ssi::curve_info(&r, i).map_err(error)?;
+        let count = usize::try_from(info.samples).map_err(|_| Error::Argument)?;
+        let s = ffi::ssi::curve_samples(&r, i, count).map_err(error)?;
+        curves.push(SsiCurve {
+            curve_type: CurveType::from_code(info.ty),
+            has_uv_a: info.flags & 1 != 0,
+            has_uv_b: info.flags & 2 != 0,
+            deviation: info.deviation,
+            interval: info.interval,
+            xyz: triples(&s.xyz),
+            uv_a: pairs(&s.uv_a),
+            uv_b: pairs(&s.uv_b),
+            a_at_uv: triples(&s.a_at_uv),
+            b_at_uv: triples(&s.b_at_uv),
+        });
+    }
+    let n = ffi::ssi::point_count(&r).map_err(error)?;
+    let mut points = Vec::new();
+    for i in 0..n {
+        let (xyz, uv_a, uv_b) = ffi::ssi::point(&r, i).map_err(error)?;
+        points.push(SsiPoint { xyz, uv_a, uv_b });
+    }
+    Ok(SsiResult { curves, points })
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    /// A plane through a sphere's centre: one great circle.
+    #[test]
+    fn plane_through_sphere_is_a_circle() {
+        let z = [0.0, 0.0, 1.0];
+        let (x, y) = ([1.0, 0.0, 0.0], [0.0, 1.0, 0.0]);
+        let plane = Surface::plane([0.0; 3], x, y, [-2.0, 2.0], [-2.0, 2.0]);
+        let sphere = Surface::sphere([0.0; 3], x, z, 1.0, [0.0, 360.0], [-90.0, 90.0]);
+        let (plane, sphere) = match (plane, sphere) {
+            (Err(Error::NotBuilt), _) | (_, Err(Error::NotBuilt)) => return,
+            (Ok(p), Ok(s)) => (p, s),
+            (p, s) => panic!("{p:?} {s:?}"),
+        };
+        let r = intersect(&plane, &sphere, Method::Global, 1e-6, 101).expect("intersect");
+        assert!(!r.curves.is_empty());
+        for c in &r.curves {
+            assert_eq!(c.xyz.len(), 101);
+            for p in &c.xyz {
+                assert!((p[0].hypot(p[1]) - 1.0).abs() < 1e-5 && p[2].abs() < 1e-5, "{p:?}");
+            }
+        }
+        assert_eq!(intersect(&plane, &sphere, Method::Global, 1e-6, 1).err(), Some(Error::Argument));
+    }
+
+    #[test]
+    fn bad_input_is_an_argument_error() {
+        match Surface::plane([0.0; 3], [1.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0], [0.0, 1.0]) {
+            Err(Error::NotBuilt) => {}
+            r => assert_eq!(r.err(), Some(Error::Argument)),
+        }
+        match Surface::nurbs([1, 1], 2, 2, &[0.0, 0.0, 1.0, 1.0], &[0.0, 0.0, 1.0, 1.0], &[[0.0; 3]; 3], None) {
+            Err(Error::NotBuilt) => {}
+            r => assert_eq!(r.err(), Some(Error::Argument)),
+        }
+    }
+}