//! Face and edge provenance found by geometry: where each face (and edge) of an operation's
//! result came from, for operations that do not yet report exact [`crate::history`].
//!
//! [`trace`] works for any operation: a result face that lies on an input face (on its surface,
//! inside its trimmed region) comes from that face; one that lies on none is generated by the
//! operation. Distinct input faces on one surface (a target's and a tool's coplanar tops) are
//! told apart by their trimmed regions, so a merged face names both and a split face's pieces
//! name the same one. Operations that know more (which profile segment a side face sweeps,
//! which edge a blend replaces) fill in generated faces' origins themselves.

use std::collections::BTreeSet;

use solidcraft_geom::{Mesh, Plane, Region2, Seg2, Vec2, Vec3};

use crate::Result;
use crate::body::Body;
use crate::history::{Change, EdgeOrigin, EdgeProv, FaceOrigin, FaceProv, History};

/// Up to this many sample points per result face.
const SAMPLES: usize = 24;

/// A face's triangles with their unit normals, and its box.
struct FaceTris {
    tris: Vec<([Vec3; 3], Vec3)>,
    lo: Vec3,
    hi: Vec3,
    area: f64,
}

fn face_tris(m: &Mesh) -> Vec<FaceTris> {
    let n = m.tri_face.iter().copied().max().map(|x| x as usize + 1).unwrap_or(0);
    let mut out: Vec<FaceTris> = (0..n)
        .map(|_| FaceTris {
            tris: Vec::new(),
            lo: Vec3::new(f64::INFINITY, f64::INFINITY, f64::INFINITY),
            hi: Vec3::new(f64::NEG_INFINITY, f64::NEG_INFINITY, f64::NEG_INFINITY),
            area: 0.0,
        })
        .collect();
    for (t, f) in m.triangles.iter().zip(&m.tri_face) {
        let (Some(slot), Some(p)) = (out.get_mut(*f as usize), m.tri(t)) else { continue };
        let c = (p[1] - p[0]).cross(p[2] - p[0]);
        let Some(nrm) = c.normalized() else { continue };
        slot.area += c.len() * 0.5;
        for q in p {
            slot.lo = slot.lo.min(q);
            slot.hi = slot.hi.max(q);
        }
        slot.tris.push((p, nrm));
    }
    out
}

fn closest_on_tri(p: Vec3, [a, b, c]: [Vec3; 3]) -> Vec3 {
    let (ab, ac, ap) = (b - a, c - a, p - a);
    let (d1, d2) = (ab.dot(ap), ac.dot(ap));
    if d1 <= 0.0 && d2 <= 0.0 {
        return a;
    }
    let bp = p - b;
    let (d3, d4) = (ab.dot(bp), ac.dot(bp));
    if d3 >= 0.0 && d4 <= d3 {
        return b;
    }
    let vc = d1 * d4 - d3 * d2;
    if vc <= 0.0 && d1 >= 0.0 && d3 <= 0.0 {
        return a + ab * (d1 / (d1 - d3));
    }
    let cp = p - c;
    let (d5, d6) = (ab.dot(cp), ac.dot(cp));
    if d6 >= 0.0 && d5 <= d6 {
        return c;
    }
    let vb = d5 * d2 - d1 * d6;
    if vb <= 0.0 && d2 >= 0.0 && d6 <= 0.0 {
        return a + ac * (d2 / (d2 - d6));
    }
    let va = d3 * d6 - d5 * d4;
    if va <= 0.0 && (d4 - d3) >= 0.0 && (d5 - d6) >= 0.0 {
        return b + (c - b) * ((d4 - d3) / ((d4 - d3) + (d5 - d6)));
    }
    let den = 1.0 / (va + vb + vc);
    a + ab * (vb * den) + ac * (vc * den)
}

/// Does point `p` with normal `n` lie on the face (within `tol`, normals along each other)?
fn on_face(f: &FaceTris, p: Vec3, n: Vec3, tol: f64) -> bool {
    if p.x < f.lo.x - tol || p.y < f.lo.y - tol || p.z < f.lo.z - tol || p.x > f.hi.x + tol || p.y > f.hi.y + tol || p.z > f.hi.z + tol {
        return false;
    }
    // Either side: a cut's tool faces come out turned inside out.
    f.tris.iter().any(|(t, tn)| tn.dot(n).abs() > 0.95 && closest_on_tri(p, *t).dist(p) <= tol)
}

/// Sample points (triangle centroids, spread out) with normals of a result face.
fn face_samples(f: &FaceTris) -> Vec<(Vec3, Vec3)> {
    let all: Vec<(Vec3, Vec3)> = f.tris.iter().map(|(t, n)| ((t[0] + t[1] + t[2]) * (1.0 / 3.0), *n)).collect();
    if all.len() <= SAMPLES {
        return all;
    }
    let step = all.len() as f64 / SAMPLES as f64;
    (0..SAMPLES).filter_map(|i| all.get((i as f64 * step) as usize).copied()).collect()
}

/// Provenance of `result` from the bodies it was made of (in operand order), by geometry.
pub fn trace(result: &Body, inputs: &[&Body]) -> Result<History> {
    let size = inputs.iter().map(|b| b.size()).fold(result.size(), f64::max).max(1e-9);
    let mesh_tol = size * 2e-4;
    let tol = mesh_tol * 2.0;
    let rm = result.tessellate(mesh_tol)?;
    let rfaces = face_tris(&rm);
    let mut imeshes = Vec::with_capacity(inputs.len());
    for b in inputs {
        imeshes.push(b.tessellate(mesh_tol)?);
    }
    let ifaces: Vec<Vec<FaceTris>> = imeshes.iter().map(face_tris).collect();
    // Each result face: the input faces most of its samples lie on.
    let mut origins: Vec<Vec<FaceOrigin>> = Vec::with_capacity(rfaces.len());
    for rf in &rfaces {
        let samples = face_samples(rf);
        let mut found: BTreeSet<FaceOrigin> = BTreeSet::new();
        for (op, faces) in ifaces.iter().enumerate() {
            for (fi, f) in faces.iter().enumerate() {
                let hits = samples.iter().filter(|(p, n)| on_face(f, *p, *n, tol)).count();
                // A face on the input face's surface and region: most of its samples land on it
                // (a merged face's samples split between its origins).
                if hits * 4 >= samples.len().max(1) || (hits > 0 && samples.len() <= 2) {
                    found.insert(FaceOrigin::Face { operand: op, face: fi });
                }
            }
        }
        origins.push(found.into_iter().collect());
    }
    // How each relates: several result faces on one origin are pieces of a split; several
    // origins are a merge; one-to-one with the same area is kept.
    let mut uses: std::collections::HashMap<FaceOrigin, usize> = std::collections::HashMap::new();
    for o in origins.iter().flatten() {
        *uses.entry(*o).or_insert(0) += 1;
    }
    let mut faces = Vec::with_capacity(origins.len());
    for (ri, os) in origins.into_iter().enumerate() {
        let change = match os.as_slice() {
            [] => Change::Generated,
            [o] if uses.get(o).copied().unwrap_or(0) > 1 => Change::Split,
            [FaceOrigin::Face { operand, face }] => {
                let a = rfaces.get(ri).map(|f| f.area).unwrap_or(0.0);
                let b = ifaces.get(*operand).and_then(|v| v.get(*face)).map(|f| f.area).unwrap_or(0.0);
                if (a - b).abs() <= 1e-4 * a.max(b).max(1e-12) { Change::Kept } else { Change::Modified }
            }
            [_] => Change::Modified,
            _ => Change::Merged,
        };
        faces.push(FaceProv { origins: os, change });
    }
    // Edges: the faces they separate, and the input edge they lie on.
    let mut edges = Vec::with_capacity(rm.edges.len());
    for (ei, pts) in rm.edges.iter().enumerate() {
        let fs: Vec<usize> = rm.edge_faces.get(ei).map(|v| v.iter().map(|x| *x as usize).collect()).unwrap_or_default();
        let probe: Vec<Vec3> = if pts.len() >= 2 {
            let k = pts.len();
            vec![pts[k / 2], (pts[0] + pts[1]) * 0.5, (pts[k - 2] + pts[k - 1]) * 0.5]
        } else {
            pts.clone()
        };
        let mut origin = None;
        'search: for (op, m) in imeshes.iter().enumerate() {
            for (ie, ipts) in m.edges.iter().enumerate() {
                if probe.iter().all(|p| ipts.windows(2).any(|w| p.dist_to_segment(w[0], w[1]) <= tol)) {
                    origin = Some((op, ie));
                    break 'search;
                }
            }
        }
        edges.push(EdgeProv { faces: fs, origins: origin.map(|(operand, edge)| EdgeOrigin::Edge { operand, edge }).into_iter().collect() });
    }
    Ok(History { operands: inputs.len(), faces, edges })
}

/// The profile segments of the loops, in [`FaceOrigin::Profile`] order (outer, then holes).
fn segments(regions: &[Region2]) -> Vec<Seg2> {
    crate::numbered::profile_loops(regions).flat_map(|l| l.segs.iter().copied()).collect()
}

/// The segment nearest the 2D points (by mean distance).
fn nearest_segment(segs: &[Seg2], pts: &[Vec2]) -> Option<usize> {
    let dist = |s: &Seg2, q: Vec2| {
        let t = s.closest_param(q);
        s.point_at(t).dist(q)
    };
    (0..segs.len()).min_by(|a, b| {
        let m = |i: usize| segs.get(i).map(|s| pts.iter().map(|q| dist(s, *q)).sum::<f64>()).unwrap_or(f64::INFINITY);
        m(*a).total_cmp(&m(*b))
    })
}

/// Origins of an extrusion's faces (`tool` from [`crate::extrude`] of `regions` on `plane`):
/// caps on the start and end planes, sides by the profile segment they sweep.
pub fn extrude_origins(tool: &Body, plane: &Plane, regions: &[Region2]) -> Result<Vec<FaceOrigin>> {
    let size = tool.size().max(1e-9);
    let m = tool.tessellate(size * 2e-4)?;
    let faces = face_tris(&m);
    let n = plane.normal();
    let segs = segments(regions);
    let heights: Vec<f64> = m.positions.iter().map(|p| (*p - plane.origin).dot(n)).collect();
    let (lo, hi) = heights.iter().fold((f64::INFINITY, f64::NEG_INFINITY), |(a, b), h| (a.min(*h), b.max(*h)));
    let mut out = Vec::with_capacity(faces.len());
    for f in &faces {
        let samples = face_samples(f);
        let cap = samples.iter().all(|(_, fn_)| fn_.dot(n).abs() > 1.0 - 1e-6);
        if cap {
            let h = samples.first().map(|(p, _)| (*p - plane.origin).dot(n)).unwrap_or(lo);
            // The start cap is the one at the low end along the plane's normal.
            out.push(FaceOrigin::Cap { end: (h - lo).abs() > (h - hi).abs() });
            continue;
        }
        let pts: Vec<Vec2> = samples.iter().map(|(p, _)| plane.to_local(*p - n * (*p - plane.origin).dot(n))).collect();
        match nearest_segment(&segs, &pts) {
            Some(k) => out.push(FaceOrigin::Profile { segment: k }),
            None => out.push(FaceOrigin::Cap { end: false }),
        }
    }
    Ok(out)
}

/// Origins of a revolution's faces (`tool` from [`crate::revolve`] of `regions` about the axis
/// through `axis_origin` along `axis_dir`, in the plane): the angle caps, and each swept face
/// by the profile segment it turns.
pub fn revolve_origins(tool: &Body, plane: &Plane, regions: &[Region2], axis_origin: Vec2, axis_dir: Vec2) -> Result<Vec<FaceOrigin>> {
    let size = tool.size().max(1e-9);
    let m = tool.tessellate(size * 2e-4)?;
    let faces = face_tris(&m);
    let segs = segments(regions);
    let o = plane.to_world(axis_origin);
    let a = plane.dir_to_world(axis_dir).normalized().ok_or_else(|| crate::KernelError::Invalid("axis".into()))?;
    // The profile's side of the axis, in the plane.
    let side = {
        let c = regions.first().map(|r| plane.to_world(r.outer.centroid())).unwrap_or(o);
        let w = c - o;
        (w - a * w.dot(a)).normalized().unwrap_or_else(|| plane.normal().cross(a))
    };
    let mut out = Vec::with_capacity(faces.len());
    for f in &faces {
        let samples = face_samples(f);
        // An angle cap is planar and contains the axis.
        let planar = samples.windows(2).all(|w| w[0].1.dot(w[1].1) > 1.0 - 1e-9);
        let n0 = samples.first().map(|x| x.1).unwrap_or(Vec3::Z);
        if planar && n0.dot(a).abs() < 1e-6 && samples.iter().all(|(p, _)| (*p - o).dot(n0).abs() < size * 1e-6) {
            // The start cap lies in the profile's plane.
            let start = samples.iter().all(|(p, _)| plane.height(*p).abs() < size * 1e-6);
            out.push(FaceOrigin::Cap { end: !start });
            continue;
        }
        // Turn each sample back into the profile's half-plane.
        let pts: Vec<Vec2> = samples
            .iter()
            .map(|(p, _)| {
                let w = *p - o;
                let h = w.dot(a);
                let r = (w - a * h).len();
                plane.to_local(o + a * h + side * r)
            })
            .collect();
        match nearest_segment(&segs, &pts) {
            Some(k) => out.push(FaceOrigin::Profile { segment: k }),
            None => out.push(FaceOrigin::Cap { end: false }),
        }
    }
    Ok(out)
}

/// Fill in generated faces of a fillet or chamfer (`result` from `input` with `edges` blended,
/// edge indices in `input`'s mesh order): each blend face comes from the edge it lies along.
pub fn blend_origins(result: &Body, input: &Body, edges: &[usize], prov: &mut History) -> Result<()> {
    let size = input.size().max(1e-9);
    let rm = result.tessellate(size * 2e-4)?;
    let im = input.tessellate(size * 2e-4)?;
    let rfaces = face_tris(&rm);
    for (fi, fp) in prov.faces.iter_mut().enumerate() {
        if fp.change != Change::Generated {
            continue;
        }
        let Some(rf) = rfaces.get(fi) else { continue };
        let samples = face_samples(rf);
        let mean = |e: usize| {
            let Some(pts) = im.edges.get(e) else { return f64::INFINITY };
            samples.iter().map(|(p, _)| pts.windows(2).map(|w| p.dist_to_segment(w[0], w[1])).fold(f64::INFINITY, f64::min)).sum::<f64>()
                / samples.len().max(1) as f64
        };
        if let Some(e) = edges.iter().copied().min_by(|a, b| mean(*a).total_cmp(&mean(*b))) {
            fp.origins = vec![FaceOrigin::Edge { operand: 0, edge: e }];
        }
    }
    Ok(())
}

/// Fill in generated faces of a shell or offset (`result` from `input`, walls `thickness`
/// inside): an inner face whose points, moved back out by the thickness, land on an input face
/// comes from that face ([`FaceOrigin::Offset`]).
pub fn offset_origins(result: &Body, input: &Body, thickness: f64, prov: &mut History) -> Result<()> {
    let size = input.size().max(1e-9);
    let tol = size * 4e-4;
    let rm = result.tessellate(size * 2e-4)?;
    let im = input.tessellate(size * 2e-4)?;
    let rfaces = face_tris(&rm);
    let ifaces = face_tris(&im);
    for (fi, fp) in prov.faces.iter_mut().enumerate() {
        if fp.change != Change::Generated {
            continue;
        }
        let Some(rf) = rfaces.get(fi) else { continue };
        let samples = face_samples(rf);
        // Inner walls face into the cavity: their normal points back toward the input face.
        let found = ifaces.iter().position(|f| samples.iter().all(|(p, n)| on_face(f, *p - *n * thickness, *n, tol)));
        if let Some(face) = found {
            fp.origins = vec![FaceOrigin::Offset { operand: 0, face }];
        }
    }
    Ok(())
}

/// Origins of any swept tool's faces (sweep, loft, pipe; extrude and revolve too) from where
/// they start: the face lying in the profile's plane is the start cap, a face with an edge in
/// that plane sweeps the profile segment under the edge, and the rest are end caps.
pub fn swept_origins(tool: &Body, plane: &Plane, regions: &[Region2]) -> Result<Vec<FaceOrigin>> {
    let size = tool.size().max(1e-9);
    let m = tool.tessellate(size * 2e-4)?;
    let faces = face_tris(&m);
    let segs = segments(regions);
    let flat = |p: Vec3| plane.height(p).abs() <= size * 1e-6;
    let mut out = Vec::with_capacity(faces.len());
    for (fi, f) in faces.iter().enumerate() {
        if f.tris.iter().all(|(t, _)| t.iter().all(|p| flat(*p))) {
            out.push(FaceOrigin::Cap { end: false });
            continue;
        }
        // Its edges in the profile's plane, as 2D points.
        let pts: Vec<Vec2> = m
            .edges
            .iter()
            .zip(&m.edge_faces)
            .filter(|(_, fs)| fs.contains(&(fi as u32)))
            .filter(|(e, _)| e.len() >= 2 && e.iter().all(|p| flat(*p)))
            .flat_map(|(e, _)| e.windows(2).map(|w| plane.to_local((w[0] + w[1]) * 0.5)).collect::<Vec<_>>())
            .collect();
        match (pts.is_empty(), nearest_segment(&segs, &pts)) {
            (false, Some(k)) => out.push(FaceOrigin::Profile { segment: k }),
            _ => out.push(FaceOrigin::Cap { end: true }),
        }
    }
    Ok(out)
}

/// Origins of a [`crate::sweep_path`] tool's faces. Its walls come in profile order (once for
/// smooth walls, once per station for ruled ones) followed by the start and end caps; checked
/// against the start cap lying in the profile's plane, else found as [`swept_origins`] does.
pub fn sweep_origins(tool: &Body, plane: &Plane, regions: &[Region2]) -> Result<Vec<FaceOrigin>> {
    let n = tool.face_count();
    let k = segments(regions).len();
    let size = tool.size().max(1e-9);
    if k > 0 && n >= 2 && (n - 2).is_multiple_of(k) {
        let m = tool.tessellate(size * 2e-4)?;
        let faces = face_tris(&m);
        let start_flat = faces.get(n - 2).is_some_and(|f| f.tris.iter().all(|(t, _)| t.iter().all(|p| plane.height(*p).abs() <= size * 1e-6)));
        if start_flat {
            return Ok((0..n)
                .map(|i| {
                    if i + 2 == n {
                        FaceOrigin::Cap { end: false }
                    } else if i + 1 == n {
                        FaceOrigin::Cap { end: true }
                    } else {
                        FaceOrigin::Profile { segment: i % k }
                    }
                })
                .collect());
        }
    }
    swept_origins(tool, plane, regions)
}

/// Replace origins on operand `operand` (a tool made by the operation) with that tool's own
/// origins (`tool[k]` for its face k): a result face from the tool's side face k comes from the
/// profile segment that swept it.
pub fn compose(prov: &mut History, operand: usize, tool: &[FaceOrigin]) {
    for fp in &mut prov.faces {
        for o in &mut fp.origins {
            if let FaceOrigin::Face { operand: op, face } = *o
                && op == operand
                && let Some(t) = tool.get(face)
            {
                *o = *t;
            }
        }
    }
}

#[cfg(test)]
mod tests {
    use solidcraft_geom::Vec3;

    use super::{Change, FaceOrigin, trace};
    use crate::{BoolOp, boolean, box_solid, cylinder};

    fn face_at(b: &crate::Body, p: Vec3) -> usize {
        let m = b.tessellate(1e-3).unwrap();
        let mut best = (f64::INFINITY, 0);
        for (t, f) in m.triangles.iter().zip(&m.tri_face) {
            let q = m.tri(t).unwrap();
            let d = super::closest_on_tri(p, q).dist(p);
            if d < best.0 {
                best = (d, *f as usize);
            }
        }
        best.1
    }

    #[test]
    fn provenance_through_booleans() {
        let a = box_solid(Vec3::ZERO, Vec3::new(40.0, 30.0, 10.0)).unwrap();
        // A through hole: the plate's top and bottom are modified, the sides kept, the hole's
        // wall comes from the tool's side.
        let pin = cylinder(Vec3::new(20.0, 15.0, -1.0), Vec3::Z, 4.0, 12.0).unwrap();
        let cut = boolean(&a, &pin, BoolOp::Cut).unwrap().unwrap();
        let p = trace(&cut, &[&a, &pin]).unwrap();
        let top = face_at(&cut, Vec3::new(5.0, 5.0, 10.0));
        let fp = &p.faces[top];
        assert_eq!(fp.origins, vec![FaceOrigin::Face { operand: 0, face: face_at(&a, Vec3::new(5.0, 5.0, 10.0)) }]);
        assert_eq!(fp.change, Change::Modified);
        let side = face_at(&cut, Vec3::new(0.0, 15.0, 5.0));
        assert_eq!(p.faces[side].change, Change::Kept);
        let wall = face_at(&cut, Vec3::new(24.0, 15.0, 5.0));
        assert!(matches!(p.faces[wall].origins[..], [FaceOrigin::Face { operand: 1, .. }]), "{:?}", p.faces[wall]);
        // A slot cut across the plate splits its top in two pieces of one origin.
        let slot = box_solid(Vec3::new(18.0, -1.0, 5.0), Vec3::new(22.0, 31.0, 11.0)).unwrap();
        let split = boolean(&a, &slot, BoolOp::Cut).unwrap().unwrap();
        let p = trace(&split, &[&a, &slot]).unwrap();
        let (l, r) = (face_at(&split, Vec3::new(5.0, 5.0, 10.0)), face_at(&split, Vec3::new(35.0, 5.0, 10.0)));
        assert_ne!(l, r);
        assert_eq!(p.faces[l].origins, p.faces[r].origins);
        assert_eq!(p.faces[l].change, Change::Split);
        // Two boxes side by side, flush on top: the union's top merges both tops.
        let b = box_solid(Vec3::new(40.0, 0.0, 0.0), Vec3::new(60.0, 30.0, 10.0)).unwrap();
        let u = boolean(&a, &b, BoolOp::Union).unwrap().unwrap();
        let p = trace(&u, &[&a, &b]).unwrap();
        let top = face_at(&u, Vec3::new(50.0, 15.0, 10.0));
        if u.face_count() == 6 {
            assert_eq!(p.faces[top].change, Change::Merged, "{:?}", p.faces[top]);
            assert_eq!(p.faces[top].origins.len(), 2);
        }
        // Edges: the plate's untouched vertical edges come from the plate's edges.
        let m = cut.tessellate(1e-3).unwrap();
        let p = trace(&cut, &[&a, &pin]).unwrap();
        let corner = (0..m.edges.len()).find(|i| m.edges[*i].iter().all(|q| q.x.abs() < 1e-9 && q.y.abs() < 1e-9)).unwrap();
        assert!(matches!(p.edges[corner].origins[..], [super::EdgeOrigin::Edge { operand: 0, .. }]));
        assert_eq!(p.edges[corner].faces.len(), 2);
    }

    #[test]
    fn provenance_of_sweeps_blends_and_shells() {
        use solidcraft_geom::{Loop2, Plane, Region2, Vec2};

        use super::{History, blend_origins, compose, extrude_origins, offset_origins, revolve_origins};
        let plane = Plane::named("XY").unwrap();
        // A 40 × 30 rectangle (segments 0: bottom, 1: right, 2: top, 3: left) with a hole.
        let rect = Loop2::polygon(&[Vec2::new(0.0, 0.0), Vec2::new(40.0, 0.0), Vec2::new(40.0, 30.0), Vec2::new(0.0, 30.0)]);
        let region = Region2 { outer: rect, holes: vec![Loop2::circle(Vec2::new(20.0, 15.0), 5.0).reversed()] };
        let tool = crate::extrude(&plane, std::slice::from_ref(&region), 0.0, 10.0).unwrap().remove(0);
        let o = extrude_origins(&tool, &plane, std::slice::from_ref(&region)).unwrap();
        assert_eq!(o[face_at(&tool, Vec3::new(40.0, 15.0, 5.0))], FaceOrigin::Profile { segment: 1 });
        assert_eq!(o[face_at(&tool, Vec3::new(20.0, 0.0, 5.0))], FaceOrigin::Profile { segment: 0 });
        assert!(matches!(o[face_at(&tool, Vec3::new(25.0, 15.0, 5.0))], FaceOrigin::Profile { segment: 4 | 5 }));
        assert_eq!(o[face_at(&tool, Vec3::new(5.0, 5.0, 0.0))], FaceOrigin::Cap { end: false });
        assert_eq!(o[face_at(&tool, Vec3::new(5.0, 5.0, 10.0))], FaceOrigin::Cap { end: true });
        // Joined to a base plate: the tool's faces keep their profile origins through the union.
        let base = box_solid(Vec3::new(-10.0, -10.0, -5.0), Vec3::new(60.0, 50.0, 0.0)).unwrap();
        let joined = boolean(&base, &tool, BoolOp::Union).unwrap().unwrap();
        let mut p = trace(&joined, &[&base, &tool]).unwrap();
        compose(&mut p, 1, &o);
        let right = face_at(&joined, Vec3::new(40.0, 15.0, 5.0));
        assert_eq!(p.faces[right].origins, vec![FaceOrigin::Profile { segment: 1 }]);
        // A revolve about the Y axis: the rectangle's right side sweeps the outer wall.
        let rv = Loop2::polygon(&[Vec2::new(10.0, 0.0), Vec2::new(20.0, 0.0), Vec2::new(20.0, 30.0), Vec2::new(10.0, 30.0)]);
        let rr = Region2 { outer: rv, holes: vec![] };
        let rt = crate::revolve(&plane, std::slice::from_ref(&rr), Vec2::ZERO, Vec2::new(0.0, 1.0), std::f64::consts::FRAC_PI_2).unwrap().remove(0);
        let ro = revolve_origins(&rt, &plane, std::slice::from_ref(&rr), Vec2::ZERO, Vec2::new(0.0, 1.0)).unwrap();
        let p45 = |r: f64, y: f64| Vec3::new(r * std::f64::consts::FRAC_1_SQRT_2, y, r * std::f64::consts::FRAC_1_SQRT_2);
        let outer = face_at(&rt, p45(20.0, 15.0));
        assert_eq!(ro[outer], FaceOrigin::Profile { segment: 1 }, "{ro:?}");
        assert_eq!(ro[face_at(&rt, p45(15.0, 30.0))], FaceOrigin::Profile { segment: 2 });
        assert_eq!(ro[face_at(&rt, Vec3::new(15.0, 15.0, 0.0))], FaceOrigin::Cap { end: false });
        // A fillet: the blend face comes from the edge it replaced.
        let b = box_solid(Vec3::ZERO, Vec3::new(40.0, 30.0, 20.0)).unwrap();
        let bm = b.tessellate(1e-3).unwrap();
        let e = (0..bm.edges.len()).find(|i| bm.edges[*i].iter().all(|q| (q.y - 30.0).abs() < 1e-9 && (q.z - 20.0).abs() < 1e-9)).unwrap();
        let f = crate::fillet(&b, &[Vec3::new(20.0, 30.0, 20.0)], 5.0).unwrap();
        let mut p: History = trace(&f, &[&b]).unwrap();
        blend_origins(&f, &b, &[e], &mut p).unwrap();
        let blend = face_at(&f, Vec3::new(20.0, 25.0 + 5.0 * std::f64::consts::FRAC_1_SQRT_2, 15.0 + 5.0 * std::f64::consts::FRAC_1_SQRT_2));
        assert_eq!(p.faces[blend].origins, vec![FaceOrigin::Edge { operand: 0, edge: e }]);
        assert_eq!(p.faces[blend].change, Change::Generated);
        // A sweep along a 3D path: each wall from the segment it starts on.
        let path: Vec<Vec3> = (0..=20).map(|i| Vec3::new(0.0, 0.0, i as f64) + Vec3::new((i as f64 * 0.2).sin() * 10.0, 0.0, 0.0)).collect();
        let sq = Loop2::polygon(&[Vec2::new(-2.0, -2.0), Vec2::new(2.0, -2.0), Vec2::new(2.0, 2.0), Vec2::new(-2.0, 2.0)]);
        let sr = Region2 { outer: sq, holes: vec![] };
        let sw = crate::sweep_path(&plane, &sr, &path).unwrap();
        let so = super::sweep_origins(&sw, &plane, std::slice::from_ref(&sr)).unwrap();
        let walls: std::collections::BTreeSet<_> =
            so.iter().filter_map(|o| if let FaceOrigin::Profile { segment } = o { Some(*segment) } else { None }).collect();
        assert_eq!(walls.len(), 4, "{so:?}");
        // The first ring of walls agrees with what their start edges say.
        let geo = super::swept_origins(&sw, &plane, std::slice::from_ref(&sr)).unwrap();
        for i in 0..4 {
            assert_eq!(so[i], geo[i], "wall {i}");
        }
        assert_eq!(so.iter().filter(|o| **o == FaceOrigin::Cap { end: false }).count(), 1);
        assert_eq!(so.iter().filter(|o| **o == FaceOrigin::Cap { end: true }).count(), 1);
        // The same for an extrusion.
        assert_eq!(super::swept_origins(&tool, &plane, std::slice::from_ref(&region)).unwrap(), o);
        // A shell: the inner walls are offsets of the outer ones.
        let sh = crate::shell(&b, &[Vec3::new(20.0, 15.0, 20.0)], 2.0).unwrap();
        let mut p = trace(&sh, &[&b]).unwrap();
        offset_origins(&sh, &b, 2.0, &mut p).unwrap();
        let inner = face_at(&sh, Vec3::new(2.0, 15.0, 10.0));
        let outer = face_at(&b, Vec3::new(0.0, 15.0, 10.0));
        assert_eq!(p.faces[inner].origins, vec![FaceOrigin::Offset { operand: 0, face: outer }]);
    }
}