geop_core_topology/contains/

shell.rs

//! Point-in-solid containment via ray casting against a shell's vertices,
//! edges, and faces in 3-D — the 3-D analog of `super::face::face_contains`.
//!
//! Coincidence with the boundary is checked first (vertex, then edge, then
//! face — each a full pass, so a higher-priority coincidence is never
//! shadowed by iteration order). Otherwise a ray is cast from the query
//! point in a random direction (drawn from a seeded PRNG, see
//! [`super::rng::Rng`]), retried with a fresh direction whenever it grazes a
//! vertex or an edge (both ambiguous to count reliably — an edge is shared
//! by two faces, a vertex by several edges), until one is found whose only
//! crossings are clean face-interior hits. The even/odd parity of that
//! crossing count then gives inside/outside — this needs no consistently
//! oriented face normal, only a clean ray.

use std::collections::HashSet;

use crate::{Body, CoedgeGeometry, EdgeId, FaceId, Model, ShellId, SolidId, VertexId};
use geop_core_geometry::{
    contains::{curve::curve_could_contain, surface::surface_could_contain},
    intersection::{
        curve_curve_intersect, curve_surface_intersect, refine_crossing,
        refine_curve_curve_crossing,
    },
    nurb_curve::NurbCurve3D,
};
use geop_core_math::{
    geop_error::{GeopError, GeopResult},
    scalars::Scalar,
    solvers::probe::{Probe, attempts},
    vector::{Vector2, Vector3},
};

use super::{
    face::{PointClassification as FaceClassification, face_contains},
    rng::Rng,
};

/// Result of classifying a query point against a shell's boundary.
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
pub enum PointClassification {
    /// The query point coincides with a vertex.
    OnVertex,
    /// The query point lies on an edge, away from its endpoints.
    OnEdge,
    /// The query point lies on a face's interior surface.
    OnFace,
    /// The query point is strictly inside the shell.
    Inside,
    /// The query point is strictly outside the shell.
    Outside,
}

const MAX_RAY_ATTEMPTS: usize = 64;

/// A shell of a model as rays are cast against it: the distinct vertices
/// and edges its faces reference, found once for every ray from a point.
pub(crate) struct ShellBoundary<'m> {
    pub(crate) model: &'m Model,
    pub(crate) shell_id: ShellId,
    pub(crate) vertices: Vec<VertexId>,
    pub(crate) edges: Vec<EdgeId>,
}

impl<'m> ShellBoundary<'m> {
    /// The boundary of `shell_id`: the distinct vertices and edges its
    /// faces reference.
    pub(crate) fn of(model: &'m Model, shell_id: ShellId) -> Self {
        let (vertices, edges) = shell_vertices_and_edges(model, shell_id);
        Self {
            model,
            shell_id,
            vertices,
            edges,
        }
    }

    /// A ray length that clears the shell from `point` (see
    /// [`ray_length_for`]).
    pub(crate) fn ray_length(&self, point: &Vector3) -> GeopResult<Scalar> {
        ray_length_for(self.model, &self.vertices, point)
    }
}

/// The distinct vertices and edges referenced by `shell_id`'s faces.
fn shell_vertices_and_edges(model: &Model, shell_id: ShellId) -> (Vec<VertexId>, Vec<EdgeId>) {
    let shell = &model.shells[&shell_id];
    let mut vertex_ids = Vec::new();
    let mut edge_ids = Vec::new();
    let mut seen_v = HashSet::new();
    let mut seen_e = HashSet::new();
    for &face_id in &shell.faces {
        for coedge_id in model.iterate_face_coedges(face_id) {
            let coedge = &model.coedges[&coedge_id];
            match coedge.geometry {
                CoedgeGeometry::Edge(edge_id) => {
                    if seen_e.insert(edge_id) {
                        edge_ids.push(edge_id);
                    }
                    let edge = &model.edges[&edge_id];
                    if seen_v.insert(edge.start_vertex) {
                        vertex_ids.push(edge.start_vertex);
                    }
                    if seen_v.insert(edge.end_vertex) {
                        vertex_ids.push(edge.end_vertex);
                    }
                }
                CoedgeGeometry::Vertex(vertex_id) => {
                    if seen_v.insert(vertex_id) {
                        vertex_ids.push(vertex_id);
                    }
                }
            }
        }
    }
    (vertex_ids, edge_ids)
}

/// A ray length guaranteed to clear `vertex_ids`' whole spatial footprint
/// from `point` — long enough that a ray this long finding zero crossings
/// unambiguously means "missed the shell entirely" (i.e. outside), not
/// "wasn't cast far enough".
fn ray_length_for(model: &Model, vertex_ids: &[VertexId], point: &Vector3) -> GeopResult<Scalar> {
    let mut max_dist_sq = Scalar::ONE;
    for &vertex_id in vertex_ids {
        let d = model.vertices[&vertex_id].point.sub(point).norm_sq();
        if d.could_be_greater(max_dist_sq) {
            max_dist_sq = d;
        }
    }
    Ok(max_dist_sq
        .sqrt()?
        .mul(Scalar::from_f64(3.0))
        .add(Scalar::ONE))
}

/// Casts a single ray from `point` in `direction` against `boundary`'s
/// vertices/edges/faces (`ray_length` as [`ShellBoundary::ray_length`]
/// gives it) and
/// returns what it crossed ([`Ray`]) — or `Ok(Err(why))` if
/// this particular direction is degenerate (grazes a vertex or edge, or
/// hits a face right on its own trim boundary) and the caller should retry
/// with a different one. Factored out of [`shell_contains`] so other
/// callers (see `validation`'s manifold check) can cast several different
/// directions from the same point without redoing its per-point setup.
///
/// `leaving` is a face of the shell that `point` lies inside of, for a ray
/// cast from a point on a face (see [`solid_contains_beside`]): the ray's
/// start on that face is where it leaves it, not a crossing, and every
/// crossing past it counts.
///
/// `seed` seeds the `face_contains` sub-calls used to classify face-interior
/// hits; callers casting multiple rays from one call site should vary it
/// per ray (e.g. by direction or attempt index) so those sub-calls don't all
/// retry along the identical degenerate path.
#[allow(clippy::too_many_arguments)]
pub(crate) fn cast_ray(
    boundary: &ShellBoundary,
    point: Vector3,
    direction: Vector3,
    ray_length: Scalar,
    leaving: Option<FaceId>,
    max_nodes: usize,
    epsilon: Scalar,
    seed: u64,
) -> GeopResult<Result<Ray, String>> {
    let ShellBoundary {
        model,
        shell_id,
        vertices: vertex_ids,
        edges: edge_ids,
    } = boundary;
    let (model, shell_id) = (*model, *shell_id);
    let shell = &model.shells[&shell_id];
    let far = point.add(&direction.prod_scalar(ray_length));
    let ray = NurbCurve3D::line(point, far)?;

    // A ray grazing a vertex is ambiguous (shared by several edges).
    for &vertex_id in vertex_ids {
        let vp = model.vertices[&vertex_id].point;
        if curve_could_contain(&ray, &vp, max_nodes, epsilon)?.is_some() {
            return Ok(Err(format!(
                "it could pass through vertex {vertex_id} at {vp:?}"
            )));
        }
    }

    // A ray crossing an edge is ambiguous (shared by two faces), and so is a
    // ray lying along one. An `Err` (the search exhausted its node budget) is
    // no answer at all, so it too asks the caller to retry. One point per
    // hit is all this needs, so `max_solutions` is 1: a coincident ray is
    // rejected either way, and a larger cap would only buy samples of it.
    for &edge_id in edge_ids {
        let full_curve = &model.edges[&edge_id].curve;
        let hits = match curve_curve_intersect(&ray, full_curve, 1, max_nodes, epsilon) {
            Ok(hits) => hits,
            Err(e) => {
                return Ok(Err(format!(
                    "its search against edge {edge_id} failed: {e}"
                )));
            }
        };
        if hits.is_coincident() {
            return Ok(Err(format!("it runs along edge {edge_id}")));
        }
        for (t_hit, s_hit) in hits.into_vec() {
            // Only a crossing *ahead* of the query point is a graze. The point
            // is not on this edge (checked by the caller), so a hit the
            // search cannot place beyond it — refined, in case it is just
            // close — is its own resolution around the query, not a crossing;
            // and were it one, the faces meeting along the edge would be hit
            // right on their trim, which they reject below.
            let (t_hit, _) = if t_hit.definitely_greater(Scalar::ZERO) {
                (t_hit, s_hit)
            } else {
                refine_curve_curve_crossing(&ray, full_curve, t_hit, s_hit)
            };
            if t_hit.definitely_greater(Scalar::ZERO) {
                return Ok(Err(format!("it could cross edge {edge_id} at t={t_hit:?}")));
            }
        }
    }

    let mut crossings = Vec::new();
    for &face_id in &shell.faces {
        let surface = &model.faces[&face_id].surface;
        // A ray lying in a face is ambiguous, and an `Err` (node budget
        // exhausted) is no answer; both ask the caller to retry. Every
        // crossing is needed for the parity count, so the cap is `max_nodes`
        // — the most crossings a search of that budget could ever report.
        // But a ray lying in the surface is rejected whatever its points, and
        // asked for `max_nodes` of them it gets that many samples of the
        // overlap: so it is asked for one first, and for every crossing only
        // if it found one and lies in no surface.
        let search = |max_solutions: usize| {
            curve_surface_intersect(&ray, surface, max_solutions, max_nodes, epsilon)
                .map_err(|e| format!("its search against face {face_id} failed: {e}"))
        };
        let hits = match search(1) {
            Ok(hits) if hits.is_coincident() => {
                return Ok(Err(format!("it lies in face {face_id}'s surface")));
            }
            Ok(hits) if hits.len() == 1 => search(max_nodes),
            other => other,
        };
        let hits = match hits {
            Ok(hits) => hits,
            Err(why) => return Ok(Err(why)),
        };
        for (t_hit, uv) in hits.into_vec() {
            // As in `face::loops_contain`: a crossing next to the query point
            // is still a crossing, so one the search cannot place beyond it
            // is refined, not dropped.
            let (t_hit, uv) = if t_hit.definitely_greater(Scalar::ZERO) {
                (t_hit, uv)
            } else {
                refine_crossing(&ray, surface, t_hit, uv)
            };
            // Sharpen — `curve_surface_intersect` honestly returns the
            // whole surviving span of its converged leaf, not an
            // arbitrarily narrowed midpoint.
            let (u, v) = (uv[0].midpoint(), uv[1].midpoint());
            let face_seed = seed ^ face_id.0;
            // A crossing that could be at the ray's far end could be on
            // either side of it. A ray cast past the shell never meets one
            // there; one ending at a point already classified (see
            // `shell_contains_from`) meets one only if that point is all but
            // on the boundary.
            if !t_hit.definitely_less(Scalar::ONE) {
                return Ok(Err(format!(
                    "its hit on face {face_id} at t={t_hit:?} could be at its far end"
                )));
            }
            match face_contains(model, face_id, u, v, max_nodes, epsilon, face_seed) {
                Ok(FaceClassification::Inside) if t_hit.definitely_greater(Scalar::ZERO) => {
                    crossings.push((face_id, t_hit, uv))
                }
                // The ray's own start, on the face it leaves.
                Ok(FaceClassification::Inside) if leaving == Some(face_id) => {}
                // A hit inside the trim that still cannot be told from the
                // query point: the query is on this face to within what the
                // numbers can tell, though the sharp query passed the caller's
                // own check (see `face::loops_contain`).
                Ok(FaceClassification::Inside) => return Ok(Ok(Ray::OnFace)),
                // Outside the trim, wherever along the ray: not a crossing.
                Ok(FaceClassification::Outside) => {}
                // A hit right on this face's own trim boundary, or an
                // outright error resolving it, is the same ambiguity as a
                // vertex/edge graze one dimension down.
                other => {
                    return Ok(Err(format!(
                        "its hit on face {face_id} at uv=({u:?}, {v:?}) classified as {other:?}"
                    )));
                }
            }
        }
    }
    Ok(Ok(Ray::Crossings(crossings)))
}

/// What a clean ray found (see [`cast_ray`]).
#[derive(Clone, Debug)]
pub(crate) enum Ray {
    /// Its start is on a face of the shell.
    OnFace,
    /// The faces it crosses, each at its parameter along the ray, in the
    /// shell's order of faces.
    Crossings(Vec<(FaceId, Scalar, Vector2)>),
}

impl Ray {
    /// The classification of the ray's start: inside the shell if it
    /// crossed the boundary an odd number of times.
    ///
    /// Unlike `face_contains` (whose ray is guaranteed to cross the outer
    /// trim loop for any query point within the surface's own bounded
    /// parameter domain), a shell occupies a bounded region of otherwise
    /// unbounded 3-D space: a ray past the shell that finds zero crossings
    /// has, by construction, missed the shell's spatial footprint entirely,
    /// which can only happen when the query point is outside — zero is a
    /// legitimate (even) count, not degenerate.
    pub(crate) fn classification(&self) -> PointClassification {
        match self {
            Ray::OnFace => PointClassification::OnFace,
            Ray::Crossings(crossings) if crossings.len() % 2 == 1 => PointClassification::Inside,
            Ray::Crossings(_) => PointClassification::Outside,
        }
    }
}

/// Classify `point` against `shell_id`'s boundary: [`PointClassification::OnVertex`] /
/// [`PointClassification::OnEdge`] / [`PointClassification::OnFace`] if the
/// query point itself coincides with the boundary, else
/// [`PointClassification::Inside`]/[`PointClassification::Outside`] via ray
/// casting.
///
/// The ray direction is drawn from a seeded PRNG (see [`Rng`]) and retried
/// (up to a bounded number of attempts, see [`cast_ray`]) until one resolves
/// cleanly. The even/odd parity of that crossing count then determines
/// inside/outside, with no dependence on any face's normal orientation
/// (faces in this crate aren't guaranteed to wind consistently).
///
/// `max_nodes` bounds both the BFS containment searches and the DFS
/// intersection searches; `epsilon` is the shared geometric tolerance;
/// `seed` seeds the direction PRNG.
pub fn shell_contains(
    model: &Model,
    shell_id: ShellId,
    point: Vector3,
    max_nodes: usize,
    epsilon: Scalar,
    seed: u64,
) -> GeopResult<PointClassification> {
    classify(model, shell_id, point, None, max_nodes, epsilon, seed)
}

/// Like [`shell_contains`], from `from`: a point already classified against
/// `shell_id`, strictly inside it (`from_inside`) or strictly outside.
///
/// The first path tried ends there, and `point` is then on the same side
/// as `from` exactly when that path crosses the boundary an even number of
/// times — the same parity argument as a ray cast past the whole shell,
/// with the answer at the far end known instead of "outside". A point near
/// `from` — another point of the same face, say — makes that path short,
/// and a short path rules out nearly every face of the shell before any
/// search runs (see `curve_could_meet_aabb`), where a ray cast past the
/// shell is tested against all of them. A degenerate one (grazing a vertex
/// or an edge, hitting a face at its trim or at an end of the path) falls
/// back to the rays of [`shell_contains`].
#[allow(clippy::too_many_arguments)]
pub fn shell_contains_from(
    model: &Model,
    shell_id: ShellId,
    point: Vector3,
    from: Vector3,
    from_inside: bool,
    max_nodes: usize,
    epsilon: Scalar,
    seed: u64,
) -> GeopResult<PointClassification> {
    classify(
        model,
        shell_id,
        point,
        Some((from, from_inside)),
        max_nodes,
        epsilon,
        seed,
    )
}

/// [`shell_contains`] and [`shell_contains_from`]: the coincidence checks,
/// then the ray to `from` if there is one, then random rays.
fn classify(
    model: &Model,
    shell_id: ShellId,
    point: Vector3,
    from: Option<(Vector3, bool)>,
    max_nodes: usize,
    epsilon: Scalar,
    seed: u64,
) -> GeopResult<PointClassification> {
    let boundary = ShellBoundary::of(model, shell_id);
    let (vertex_ids, edge_ids) = (&boundary.vertices, &boundary.edges);

    // Coincidence pre-check: vertex, then edge, then face — each a full
    // pass, so higher-priority coincidences are deterministic regardless of
    // iteration order.
    for &vertex_id in vertex_ids {
        if model.vertices[&vertex_id].point.could_be_equal(&point) {
            return Ok(PointClassification::OnVertex);
        }
    }
    for &edge_id in edge_ids {
        if curve_could_contain(&model.edges[&edge_id].curve, &point, max_nodes, epsilon)?.is_some()
        {
            return Ok(PointClassification::OnEdge);
        }
    }
    if shell_face_at(model, shell_id, &point, max_nodes, epsilon, seed)?.is_some() {
        return Ok(PointClassification::OnFace);
    }

    let mut last_rejection = String::new();
    if let Some((from, from_inside)) = from {
        // `point` is inside if it is on `from`'s side and `from` is, or on
        // the other side and `from` is not.
        let inside = |other_side: bool| match other_side != from_inside {
            true => PointClassification::Inside,
            false => PointClassification::Outside,
        };
        // Two boxes that overlap, neither touching the boundary, are one
        // connected region clear of it: the same side.
        if point.could_be_equal(&from) {
            return Ok(inside(false));
        }
        // Not straight to `from`, but by way of a point `via` off to a
        // random side, as far from the middle as the two are apart: a
        // straight segment between two points of one face lies in that
        // face's plane, and in any face of the shell flush with it; and
        // points that structured, the corners of a symmetric face, say, can
        // as well graze a cylinder, a crossing that counts once and is
        // none. A random detour is as clear of these as a random ray is.
        let mut rng = Rng::new(seed);
        let gap = from.sub(&point).norm().div(Scalar::TWO)?;
        let via = point
            .add(&from)
            .prod_scalar(Scalar::ONE.div(Scalar::TWO)?)
            .add(&rng.next_direction3().prod_scalar(gap))
            .sharpen();
        let leg = |start: Vector3| {
            cast_ray(
                &boundary,
                start,
                via.sub(&start),
                Scalar::ONE,
                None,
                max_nodes,
                epsilon,
                seed,
            )
        };
        // The leg from `point` may find it on the boundary; the one from
        // `from`, already classified clear of it, finding it there is
        // ambiguous.
        match (
            leg(point)?.map(|ray| ray.classification()),
            leg(from)?.map(|ray| ray.classification()),
        ) {
            (
                Ok(
                    on @ (PointClassification::OnFace
                    | PointClassification::OnEdge
                    | PointClassification::OnVertex),
                ),
                _,
            ) => {
                return Ok(on);
            }
            (
                Ok(first),
                Ok(second @ (PointClassification::Inside | PointClassification::Outside)),
            ) => {
                let odd = |c: PointClassification| c == PointClassification::Inside;
                return Ok(inside(odd(first) != odd(second)));
            }
            (Err(reason), _) | (_, Err(reason)) => {
                last_rejection = format!("the detour to {from:?} was rejected: {reason}");
            }
            (_, Ok(on)) => {
                last_rejection =
                    format!("the detour to {from:?} found that point {on:?}, though classified");
            }
        }
    }

    let ray_length = boundary.ray_length(&point)?;

    let mut rng = Rng::new(seed);
    // Directions tried in turn until one is clear.
    let found = attempts(MAX_RAY_ATTEMPTS, |&attempt| -> GeopResult<_> {
        let direction = rng.next_direction3();
        let attempt_seed = seed ^ (attempt as u64).wrapping_mul(0x9E3779B97F4A7C15);
        match cast_ray(
            &boundary,
            point,
            direction,
            ray_length,
            None,
            max_nodes,
            epsilon,
            attempt_seed,
        )? {
            Ok(ray) => Ok(Probe::Found(ray.classification())),
            Err(reason) => {
                last_rejection = format!("the ray along {direction:?} was rejected: {reason}");
                Ok(Probe::Next)
            }
        }
    })?;
    if let Some(classification) = found {
        return Ok(classification);
    }
    Err(GeopError::new(format!(
        "shell_contains: could not find a ray direction clear of every vertex and edge after many \
         attempts; the last one was rejected because {last_rejection}"
    )))
}

/// What [`solid_contains_beside`] found: whether the solid is on the
/// normal's side of the face, and the ray that decided it, for the record.
#[derive(Clone, Debug)]
pub struct Beside {
    pub inside: bool,
    pub direction: Vector3,
    /// The faces the ray crossed past its start, each at its parameter
    /// along the ray.
    pub crossings: Vec<(FaceId, Scalar, Vector2)>,
}

/// Classifies the points of the solid `solid_id` just off its face
/// `face_id` (of its shell `shell_id`) on the side `normal` points to, at
/// `point`, a point inside the face: `inside` if the solid is on that side
/// of the face there. Other arguments as for [`shell_contains`].
///
/// Decided by a ray from `point` in a random direction into the normal's
/// half space: every point of the ray just past the face is on that side,
/// and the parity of the ray's crossings past the face, over every shell of
/// the solid, classifies them all at once (see [`cast_ray`], which does not
/// count the ray's start on the face it leaves). No point a chosen distance
/// off the face is asked about: that asks about a point that far away, and
/// on a non-convex solid a distance that looked small reaches through the
/// local material into another part of it. A bore 18 from the side of a
/// block, probed at 29 both ways, has the block on the normal's side and
/// the air past the block behind, and its wall reads as inverted.
#[allow(clippy::too_many_arguments)]
pub fn solid_contains_beside(
    model: &Model,
    solid_id: SolidId,
    shell_id: ShellId,
    face_id: FaceId,
    point: Vector3,
    normal: Vector3,
    max_nodes: usize,
    epsilon: Scalar,
    seed: u64,
) -> GeopResult<Beside> {
    let boundaries: Vec<ShellBoundary> = model
        .get_solid(solid_id)?
        .shells
        .iter()
        .map(|&s| ShellBoundary::of(model, s))
        .collect();
    let mut rng = Rng::new(seed);
    let mut last_rejection = String::new();
    let found = attempts(MAX_RAY_ATTEMPTS, |&attempt| -> GeopResult<_> {
        let drawn = rng.next_direction3();
        let along = drawn.prod_dot(&normal);
        // Turned into the normal's half space. One that could run along
        // the face leaves it to neither side.
        let direction = if along.definitely_greater(Scalar::ZERO) {
            drawn
        } else if along.definitely_less(Scalar::ZERO) {
            drawn.neg()
        } else {
            last_rejection = format!("the direction {drawn:?} could run along the face");
            return Ok(Probe::Next);
        };
        let attempt_seed = seed ^ (attempt as u64).wrapping_mul(0x9E3779B97F4A7C15);
        let mut crossings = Vec::new();
        for boundary in &boundaries {
            let leaving = (boundary.shell_id == shell_id).then_some(face_id);
            match cast_ray(
                boundary,
                point,
                direction,
                boundary.ray_length(&point)?,
                leaving,
                max_nodes,
                epsilon,
                attempt_seed,
            )? {
                Ok(Ray::Crossings(crossed)) => crossings.extend(crossed),
                Ok(Ray::OnFace) => {
                    last_rejection = format!(
                        "the ray along {direction:?} found the point on a face of shell {}",
                        boundary.shell_id
                    );
                    return Ok(Probe::Next);
                }
                Err(reason) => {
                    last_rejection = format!("the ray along {direction:?} was rejected: {reason}");
                    return Ok(Probe::Next);
                }
            }
        }
        // Inside an odd number of shells: the shells of one solid never
        // cross, so that is inside the outer shell and no void.
        Ok(Probe::Found(Beside {
            inside: crossings.len() % 2 == 1,
            direction,
            crossings,
        }))
    })?;
    found.ok_or_else(|| {
        GeopError::new(format!(
            "solid_contains_beside: no ray from {point:?} off face {face_id} of solid {solid_id} \
             was clear of every vertex and edge after {MAX_RAY_ATTEMPTS} attempts; the last one \
             was rejected because {last_rejection}"
        ))
    })
}

/// Whether `point`, a point inside face `face_id` of the solid `solid_id`,
/// is on the solid's boundary anywhere else as well: on a vertex, on an
/// edge, or on another face, of any of its shells. A body touching the face
/// at a point leaves a vertex of its own in the face's interior with no
/// topology between them: an inscribed sphere's pole on the cube's face.
/// Other arguments as for [`shell_contains`].
pub fn on_solid_boundary_besides(
    model: &Model,
    solid_id: SolidId,
    face_id: FaceId,
    point: &Vector3,
    max_nodes: usize,
    epsilon: Scalar,
    seed: u64,
) -> GeopResult<bool> {
    for &shell_id in &model.get_solid(solid_id)?.shells {
        let boundary = ShellBoundary::of(model, shell_id);
        for &vertex_id in &boundary.vertices {
            if model.get_vertex(vertex_id)?.point.could_be_equal(point) {
                return Ok(true);
            }
        }
        for &edge_id in &boundary.edges {
            if curve_could_contain(&model.get_edge(edge_id)?.curve, point, max_nodes, epsilon)?
                .is_some()
            {
                return Ok(true);
            }
        }
        for &other in &model.get_shell(shell_id)?.faces {
            if other == face_id {
                continue;
            }
            let Some((u, v)) =
                surface_could_contain(&model.get_face(other)?.surface, point, max_nodes, epsilon)?
            else {
                continue;
            };
            if !matches!(
                face_contains(model, other, u, v, max_nodes, epsilon, seed ^ other.0)?,
                FaceClassification::Outside
            ) {
                return Ok(true);
            }
        }
    }
    Ok(false)
}

/// The face of `shell_id` that `point` lies on, if any, and where on its
/// surface: the first, in the shell's order.
///
/// Proximity to a face's *untrimmed* surface is not membership (see
/// AGENTS.md): a face's surface generally extends well past its own trim
/// loop (two faces split from one larger one still share that surface), so
/// a point can sit squarely on the surface nowhere near the face. The face
/// has to contain the point within its trim, or on it.
pub fn shell_face_at(
    model: &Model,
    shell_id: ShellId,
    point: &Vector3,
    max_nodes: usize,
    epsilon: Scalar,
    seed: u64,
) -> GeopResult<Option<(FaceId, Scalar, Scalar)>> {
    for &face_id in &model.get_shell(shell_id)?.faces {
        let Some((u, v)) =
            surface_could_contain(&model.faces[&face_id].surface, point, max_nodes, epsilon)?
        else {
            continue;
        };
        if !matches!(
            face_contains(model, face_id, u, v, max_nodes, epsilon, seed ^ face_id.0)?,
            FaceClassification::Outside
        ) {
            return Ok(Some((face_id, u, v)));
        }
    }
    Ok(None)
}

/// Points already classified against shells, each strictly inside or
/// outside: what [`solid_contains`] classifies later points from, along a
/// short path from the nearest (see [`shell_contains_from`]) instead of a
/// ray across the whole shell.
#[derive(Clone, Debug, Default)]
pub struct Classified {
    points: Vec<(ShellId, Vector3, bool)>,
}

impl Classified {
    /// The point classified against `shell` nearest to `point`, and whether
    /// it is inside: which one is a free choice, so it is chosen by the
    /// midpoints.
    fn nearest(&self, shell: ShellId, point: &Vector3) -> Option<(Vector3, bool)> {
        let distance = |p: &Vector3| p.sub(point).norm_sq().to_f64();
        self.points
            .iter()
            .filter(|(s, _, _)| *s == shell)
            .min_by(|a, b| distance(&a.1).total_cmp(&distance(&b.1)))
            .map(|&(_, p, inside)| (p, inside))
    }
}

/// The lumps of the solid `solid_id`: each shell bounding material from
/// outside, with the voids in it. A shell is a void when inside an odd
/// number of the solid's other shells, and a void belongs to the
/// innermost shell around it. Asked of each shell at one of its vertices
/// off the other shell: shells of one solid never cross, so where one
/// vertex is, all of it is. A solid of one shell is one lump without
/// voids. Lumps are in the order of their shells in the solid. Other
/// arguments as for [`shell_contains`].
pub fn solid_lumps(
    model: &Model,
    solid_id: SolidId,
    max_nodes: usize,
    epsilon: Scalar,
    seed: u64,
) -> GeopResult<Vec<(ShellId, Vec<ShellId>)>> {
    let shells = model.get_solid(solid_id)?.shells.clone();
    if shells.len() == 1 {
        return Ok(vec![(shells[0], Vec::new())]);
    }
    // Which other shells hold each shell.
    let mut around: Vec<Vec<usize>> = vec![Vec::new(); shells.len()];
    for (i, &shell) in shells.iter().enumerate() {
        let vertices: Vec<VertexId> = model.iter_body_vertices(Body::Sheet(shell))?.collect();
        for (j, &other) in shells.iter().enumerate() {
            if i == j {
                continue;
            }
            let mut decided = None;
            for &vertex in &vertices {
                let point = model.get_vertex(vertex)?.point;
                match shell_contains(model, other, point, max_nodes, epsilon, seed)? {
                    PointClassification::Inside => decided = Some(true),
                    PointClassification::Outside => decided = Some(false),
                    _ => continue,
                }
                break;
            }
            match decided {
                Some(true) => around[i].push(j),
                Some(false) => {}
                None => {
                    return Err(GeopError::new(format!(
                        "solid_lumps: every vertex of shell {shell} of solid {solid_id} lies on shell {other}: the two touch everywhere, and neither is inside the other"
                    )));
                }
            }
        }
    }
    let depth: Vec<usize> = around.iter().map(Vec::len).collect();
    let mut lumps: Vec<(ShellId, Vec<ShellId>)> = Vec::new();
    let mut lump_of: Vec<Option<usize>> = vec![None; shells.len()];
    for i in 0..shells.len() {
        if depth[i].is_multiple_of(2) {
            lump_of[i] = Some(lumps.len());
            lumps.push((shells[i], Vec::new()));
        }
    }
    for i in 0..shells.len() {
        if !depth[i].is_multiple_of(2) {
            let parent = around[i]
                .iter()
                .copied()
                .find(|&j| depth[j] + 1 == depth[i])
                .and_then(|j| lump_of[j])
                .ok_or_else(|| {
                    GeopError::new(format!(
                        "solid_lumps: shell {} of solid {solid_id} is a void inside no lump",
                        shells[i]
                    ))
                })?;
            lumps[parent].1.push(shells[i]);
        }
    }
    Ok(lumps)
}

/// Classify `point` against the solid `solid_id`: on its boundary if on any
/// of its shells', else inside exactly when inside an odd number of them —
/// inside the outer shell and no void, say — for the shells of one solid
/// never cross. Asking each shell alone would put a point in a void inside
/// the solid. Each shell is asked from the nearest point `classified`
/// already holds for it, if any (see [`shell_contains_from`]), and a point
/// found strictly inside or outside a shell is added to it. Other arguments
/// as for [`shell_contains`].
pub fn solid_contains(
    model: &Model,
    solid_id: SolidId,
    point: Vector3,
    max_nodes: usize,
    epsilon: Scalar,
    seed: u64,
    classified: &mut Classified,
) -> GeopResult<PointClassification> {
    let mut inside = false;
    for &shell_id in &model.get_solid(solid_id)?.shells {
        let from = classified.nearest(shell_id, &point);
        let class = classify(model, shell_id, point, from, max_nodes, epsilon, seed)?;
        match class {
            PointClassification::Inside | PointClassification::Outside => {
                let is_inside = class == PointClassification::Inside;
                classified.points.push((shell_id, point, is_inside));
                inside ^= is_inside;
            }
            on => return Ok(on),
        }
    }
    Ok(if inside {
        PointClassification::Inside
    } else {
        PointClassification::Outside
    })
}

/// The face of the solid `solid_id` that `point` lies on, if any, and where
/// on its surface (see [`shell_face_at`]).
pub fn solid_face_at(
    model: &Model,
    solid_id: SolidId,
    point: &Vector3,
    max_nodes: usize,
    epsilon: Scalar,
    seed: u64,
) -> GeopResult<Option<(FaceId, Scalar, Scalar)>> {
    for &shell_id in &model.get_solid(solid_id)?.shells {
        if let Some(at) = shell_face_at(model, shell_id, point, max_nodes, epsilon, seed)? {
            return Ok(Some(at));
        }
    }
    Ok(None)
}

#[cfg(test)]
mod tests {
    use super::{PointClassification, shell_contains, shell_contains_from};
    use crate::{
        Coedge, CoedgeGeometry, CoedgeId, Edge, EdgeId, Face, FaceId, Model, Sense, Shell, ShellId,
        Vertex, VertexId, boundary::BoundaryType,
    };
    use geop_core_geometry::{
        nurb_curve::{NurbCurve, NurbCurve2D, NurbCurve3D},
        nurb_surface::NurbSurface3D,
    };
    use geop_core_math::{
        scalars::Scalar,
        vector::{Vector3, Vector4},
    };

    const MAX: usize = 200;
    const MIN_SUBDIVISION_SIZE: f64 = 1e-3;
    const SEED: u64 = 424_242;

    type P3 = (f64, f64, f64);

    /// A single planar quad face, corners given in trim-loop order
    /// (`p00 -> p10 -> p11 -> p01`), each face independently vertexed/edged
    /// (geometrically watertight, topologically not shared — fine for pure
    /// containment testing).
    fn quad_face(model: &mut Model, shell: ShellId, p00: P3, p10: P3, p11: P3, p01: P3) -> FaceId {
        let p4 = |p: P3| {
            Vector4::from_array([
                Scalar::from_f64(p.0),
                Scalar::from_f64(p.1),
                Scalar::from_f64(p.2),
                Scalar::ONE,
            ])
        };
        let surface = NurbSurface3D::try_new(
            1,
            1,
            vec![p4(p00), p4(p01), p4(p10), p4(p11)],
            vec![Scalar::ZERO, Scalar::ZERO, Scalar::ONE, Scalar::ONE],
            vec![Scalar::ZERO, Scalar::ZERO, Scalar::ONE, Scalar::ONE],
        )
        .unwrap();
        let face_id = model.insert_face(Face {
            surface,
            outer: BoundaryType::Vertex(crate::VertexId(0)),
            holes: Vec::new(),
            shell,
        });

        let corners = [p00, p10, p11, p01];
        let verts: Vec<VertexId> = corners
            .iter()
            .map(|&(x, y, z)| {
                model.insert_vertex(Vertex {
                    point: Vector3::from_array([
                        Scalar::from_f64(x),
                        Scalar::from_f64(y),
                        Scalar::from_f64(z),
                    ]),
                })
            })
            .collect();
        let edges: Vec<EdgeId> = (0..4)
            .map(|i| {
                model.insert_edge(Edge {
                    curve: NurbCurve3D::try_new(
                        1,
                        vec![p4(corners[i]), p4(corners[(i + 1) % 4])],
                        vec![Scalar::ZERO, Scalar::ZERO, Scalar::ONE, Scalar::ONE],
                    )
                    .unwrap(),
                    start_vertex: verts[i],
                    end_vertex: verts[(i + 1) % 4],
                })
            })
            .collect();
        let uv = [(0.0, 0.0), (1.0, 0.0), (1.0, 1.0), (0.0, 1.0)];
        let pc = |a: (f64, f64), b: (f64, f64)| -> NurbCurve2D {
            NurbCurve::try_new(
                1,
                vec![
                    Vector3::from_array([
                        Scalar::from_f64(a.0),
                        Scalar::from_f64(a.1),
                        Scalar::ONE,
                    ]),
                    Vector3::from_array([
                        Scalar::from_f64(b.0),
                        Scalar::from_f64(b.1),
                        Scalar::ONE,
                    ]),
                ],
                vec![Scalar::ZERO, Scalar::ZERO, Scalar::ONE, Scalar::ONE],
            )
            .unwrap()
        };
        let coedges: Vec<CoedgeId> = (0..4)
            .map(|i| {
                model.insert_coedge(Coedge {
                    geometry: CoedgeGeometry::Edge(edges[i]),
                    sense: Sense::Forward,
                    pcurve: pc(uv[i], uv[(i + 1) % 4]),
                    next: CoedgeId(0),
                    prev: CoedgeId(0),
                    face: face_id,
                })
            })
            .collect();
        for i in 0..4 {
            model.coedges.get_mut(&coedges[i]).unwrap().next = coedges[(i + 1) % 4];
            model.coedges.get_mut(&coedges[i]).unwrap().prev = coedges[(i + 3) % 4];
        }
        model.faces.get_mut(&face_id).unwrap().outer = BoundaryType::Loop(coedges[0]);

        face_id
    }

    /// The unit cube `[0,1]^3`, as six independently-vertexed quad faces
    /// (mixed winding across faces on purpose — containment must not depend
    /// on consistent orientation).
    fn unit_cube(model: &mut Model) -> ShellId {
        let shell_id = model.insert_shell(Shell {
            faces: vec![],
            solid: None,
        });
        let faces = vec![
            quad_face(
                model,
                shell_id,
                (0., 0., 1.),
                (1., 0., 1.),
                (1., 1., 1.),
                (0., 1., 1.),
            ), // +Z
            quad_face(
                model,
                shell_id,
                (0., 0., 0.),
                (0., 1., 0.),
                (1., 1., 0.),
                (1., 0., 0.),
            ), // -Z (reversed winding)
            quad_face(
                model,
                shell_id,
                (1., 0., 0.),
                (1., 1., 0.),
                (1., 1., 1.),
                (1., 0., 1.),
            ), // +X
            quad_face(
                model,
                shell_id,
                (0., 0., 0.),
                (0., 0., 1.),
                (0., 1., 1.),
                (0., 1., 0.),
            ), // -X (reversed winding)
            quad_face(
                model,
                shell_id,
                (0., 1., 0.),
                (1., 1., 0.),
                (1., 1., 1.),
                (0., 1., 1.),
            ), // +Y
            quad_face(
                model,
                shell_id,
                (0., 0., 0.),
                (1., 0., 0.),
                (1., 0., 1.),
                (0., 0., 1.),
            ), // -Y
        ];
        model.shells.get_mut(&shell_id).unwrap().faces = faces;
        shell_id
    }

    #[test]
    fn center_is_inside() {
        let mut model = Model::new();
        let shell_id = unit_cube(&mut model);
        let p = Vector3::from_array([Scalar::from_f64(0.5); 3]);
        assert_eq!(
            shell_contains(
                &model,
                shell_id,
                p,
                MAX,
                Scalar::from_f64(MIN_SUBDIVISION_SIZE),
                SEED
            )
            .unwrap(),
            PointClassification::Inside
        );
    }

    #[test]
    fn far_point_is_outside() {
        let mut model = Model::new();
        let shell_id = unit_cube(&mut model);
        let p = Vector3::from_array([
            Scalar::from_f64(-5.0),
            Scalar::from_f64(0.5),
            Scalar::from_f64(0.5),
        ]);
        assert_eq!(
            shell_contains(
                &model,
                shell_id,
                p,
                MAX,
                Scalar::from_f64(MIN_SUBDIVISION_SIZE),
                SEED
            )
            .unwrap(),
            PointClassification::Outside
        );
    }

    #[test]
    fn point_just_outside_face_is_outside() {
        let mut model = Model::new();
        let shell_id = unit_cube(&mut model);
        let p = Vector3::from_array([
            Scalar::from_f64(-0.1),
            Scalar::from_f64(0.5),
            Scalar::from_f64(0.5),
        ]);
        assert_eq!(
            shell_contains(
                &model,
                shell_id,
                p,
                MAX,
                Scalar::from_f64(MIN_SUBDIVISION_SIZE),
                SEED
            )
            .unwrap(),
            PointClassification::Outside
        );
    }

    #[test]
    fn face_point_is_on_face() {
        let mut model = Model::new();
        let shell_id = unit_cube(&mut model);
        let p = Vector3::from_array([Scalar::ZERO, Scalar::from_f64(0.5), Scalar::from_f64(0.5)]);
        assert_eq!(
            shell_contains(
                &model,
                shell_id,
                p,
                MAX,
                Scalar::from_f64(MIN_SUBDIVISION_SIZE),
                SEED
            )
            .unwrap(),
            PointClassification::OnFace
        );
    }

    #[test]
    fn edge_point_is_on_edge() {
        let mut model = Model::new();
        let shell_id = unit_cube(&mut model);
        let p = Vector3::from_array([Scalar::ZERO, Scalar::ZERO, Scalar::from_f64(0.5)]);
        assert_eq!(
            shell_contains(
                &model,
                shell_id,
                p,
                MAX,
                Scalar::from_f64(MIN_SUBDIVISION_SIZE),
                SEED
            )
            .unwrap(),
            PointClassification::OnEdge
        );
    }

    #[test]
    fn vertex_point_is_on_vertex() {
        let mut model = Model::new();
        let shell_id = unit_cube(&mut model);
        let p = Vector3::from_array([Scalar::ZERO, Scalar::ZERO, Scalar::ZERO]);
        assert_eq!(
            shell_contains(
                &model,
                shell_id,
                p,
                MAX,
                Scalar::from_f64(MIN_SUBDIVISION_SIZE),
                SEED
            )
            .unwrap(),
            PointClassification::OnVertex
        );
    }

    /// Classified from a point already classified, inside or out, every
    /// point gets what `shell_contains` gives it — points in the cube's
    /// plane of symmetry and on its faces' planes among them, whose
    /// straight segments to each other lie in a face's plane or run
    /// through an edge.
    #[test]
    fn classified_from_a_known_point() {
        let mut model = Model::new();
        let shell_id = unit_cube(&mut model);
        let p = |x: f64, y: f64, z: f64| Vector3::from_array([x, y, z].map(Scalar::from_f64));
        let points = [
            p(0.5, 0.5, 0.5),
            p(0.25, 0.5, 0.5),
            p(-0.5, 0.5, 0.5),
            p(1.5, 0.5, 0.5),
            p(-0.5, -0.5, 0.5),
            p(0.5, 0.5, 2.0),
            p(-0.5, 0.0, 0.0),
            p(0.0, -0.5, 0.0),
            p(0.0, 0.0, 0.0),
            p(0.0, 0.5, 0.5),
        ];
        let eps = Scalar::from_f64(MIN_SUBDIVISION_SIZE);
        for from in &points[..6] {
            let from_class = shell_contains(&model, shell_id, *from, MAX, eps, SEED).unwrap();
            let from_inside = from_class == PointClassification::Inside;
            for (k, point) in points.iter().enumerate() {
                let want = shell_contains(&model, shell_id, *point, MAX, eps, SEED).unwrap();
                let got = shell_contains_from(
                    &model,
                    shell_id,
                    *point,
                    *from,
                    from_inside,
                    MAX,
                    eps,
                    SEED ^ k as u64,
                )
                .unwrap();
                assert_eq!(got, want, "{point:?} from {from:?}");
            }
        }
    }
}