Merge pull request #108 from CypherNaught-0x/fix/basin-rim-blends

Fix fillets and chamfers around sloped pocket rims

crates/kernel/src/blend.rs

@@ -134,7 +134,7 @@ fn setback(r: f64, phi: f64, shape: Shape) -> f64 {
 /// An extended corner edge must stay in the trimmed side face, not just its infinite
 /// supporting plane. Reject crossings and contact with any boundary after the old corner,
 /// including holes. Curved boundaries are sampled more finely than the clearance tolerance.
-fn chamfer_extension_fits(face: &mt::Face, p: Vec3, q: Vec3, tol: f64) -> bool {
+pub(crate) fn extension_fits(face: &mt::Face, p: Vec3, q: Vec3, tol: f64) -> bool {
     use mt::{BoundedCurve, ParameterDivision1D};
     // truck requires a curve sampling tolerance of at least 1e-6. Keep the
     // boundary clearance four times larger than the chord approximation error.
@@ -260,7 +260,7 @@ fn blend_geometry(
             // vertex. This is valid only within the trimmed side face.
             if (q - p - dir * along).len() > tol * 100.0 + 1e-7
                 || along >= len - tol
-                || (along < 0.0 && !faces.get(fi).is_some_and(|f| chamfer_extension_fits(f, p, q, tol)))
+                || (along < 0.0 && !faces.get(fi).is_some_and(|f| extension_fits(f, p, q, tol)))
             {
                 return Err(unsupported("the blend is larger than the neighbouring faces"));
             }
@@ -785,6 +785,14 @@ pub fn chamfer(body: &Body, edges: &[Vec3], distance: f64) -> Result<Body> {
             return Ok(r);
         }
     }
+    // At a pocket rim the cutters cannot extend beyond the selected edges. Rebuild
+    // the corners locally, allowing adjacent chamfers to reach different depths
+    // on their common wall seam.
+    if edges.len() > 1
+        && let Ok(result) = crate::blend_corner::chamfer(body, edges, distance)
+    {
+        return Ok(result);
+    }
     // Coplanar neighbours as one face (as after a join), so corners look as they should.
     let mut cur = Body::new(crate::heal::heal_keep(body.deep_copy(), body.size(), body.split_keep()))?;
     if let Some(groups) = loop_groups(&cur, edges)? {
@@ -961,7 +969,7 @@ mod tests {
         for plane in [Plane::XY, Plane::XZ] {
             let face = crate::build::face(&plane, &region).unwrap();
             let p = plane.to_world(Vec2::new(4.0, 4.0));
-            let fits = |y| chamfer_extension_fits(&face, p, plane.to_world(Vec2::new(4.0, y)), 1e-6);
+            let fits = |y| extension_fits(&face, p, plane.to_world(Vec2::new(4.0, y)), 1e-6);
             assert!(fits(2.0));
             assert!(!fits(1.5)); // endpoint in the circular hole
             assert!(!fits(0.5)); // endpoint in material, but the extension crossed the hole

crates/kernel/src/blend_corner.rs

@@ -1,7 +1,7 @@
 //! Fillets with a radius per edge, and the corners where filleted edges meet.
 //!
 //! The edges are straight and convex, and every face at their end vertices is planar with three
-//! edges meeting at each such vertex (box-like corners). The faces there are rebuilt around the
+//! edges meeting at each such vertex, including reflex corners of pocket rims. The faces there are rebuilt around the
 //! blends, faces away from them are kept. At a vertex:
 //! - one filleted edge ends on the third face: the blend ends in that face's section of the
 //!   cylinder (an exact arc or ellipse);
@@ -38,7 +38,7 @@ pub fn fillet_radii(body: &Body, edges: &[(Vec3, f64)]) -> Result<Body> {
     let mut first: Option<KernelError> = None;
     // Earlier rounds at the edges' ends come off and are blended again with the new edges.
     for _ in 0..4 {
-        let err = match build(&cur, &sel) {
+        let err = match build(&cur, &sel, true) {
             Ok(b) => return Ok(b),
             Err(e) => e,
         };
@@ -54,6 +54,14 @@ pub fn fillet_radii(body: &Body, edges: &[(Vec3, f64)]) -> Result<Body> {
     Err(first.unwrap_or_else(|| unsupported("these fillets")))
 }
 
+/// Local equal-distance chamfers, including pocket rims where adjacent walls have
+/// different slopes. The planar strips trim each other at each selected corner.
+pub(crate) fn chamfer(body: &Body, edges: &[Vec3], distance: f64) -> Result<Body> {
+    let cur = Body::new(crate::heal::heal_keep(body.deep_copy(), body.size(), body.split_keep()))?;
+    let selected: Vec<_> = edges.iter().map(|p| (*p, distance)).collect();
+    build(&cur, &selected, false)
+}
+
 /// A face of the shell being rebuilt, or a new face.
 #[derive(Clone, Copy, PartialEq)]
 enum Owner {
@@ -112,7 +120,7 @@ fn edge(a: &mt::Vertex, b: &mt::Vertex, c: mt::Curve) -> Result<mt::Edge> {
 /// along `u`.
 fn hit(v: Vec3, u: Vec3, t: Vec3, s: f64) -> Result<Vec3> {
     let k = u.dot(t);
-    if k < 1e-6 {
+    if k.abs() < 1e-6 {
         return Err(unsupported("a fillet that runs along its neighbouring edge"));
     }
     Ok(v + u * (s / k))
@@ -122,8 +130,8 @@ fn hit(v: Vec3, u: Vec3, t: Vec3, s: f64) -> Result<Vec3> {
 /// directions `t1`, `t2`, setbacks `s1`, `s2`) cross.
 fn meet(v: Vec3, (u1, t1, s1): (Vec3, Vec3, f64), (u2, t2, s2): (Vec3, Vec3, f64)) -> Result<Vec3> {
     let (k1, k2) = (u1.dot(t2), u2.dot(t1));
-    if k1 < 1e-6 || k2 < 1e-6 {
-        return Err(unsupported("fillets of edges that meet at a reflex corner"));
+    if k1.abs() < 1e-6 || k2.abs() < 1e-6 {
+        return Err(unsupported("fillets of nearly parallel edges at a corner"));
     }
     Ok(v + u1 * (s2 / k1) + u2 * (s1 / k2))
 }
@@ -174,7 +182,7 @@ fn connect(from: &mt::Vertex, to: &mt::Vertex, pool: &[(mt::Edge, Vec<Owner>)],
     Err(KernelError::Failed("fillet corner: a boundary does not close".into()))
 }
 
-fn build(cur: &Body, sel_pts: &[(Vec3, f64)]) -> Result<Body> {
+fn build(cur: &Body, sel_pts: &[(Vec3, f64)], round: bool) -> Result<Body> {
     let size = cur.size();
     let tol = (size * 2e-3).max(1e-3);
     let solid = cur.deep_copy();
@@ -252,16 +260,20 @@ fn build(cur: &Body, sel_pts: &[(Vec3, f64)]) -> Result<Body> {
             return Err(unsupported("the faces at the edge are tangent"));
         }
         let r = radius.get(id).copied().unwrap_or(1.0);
-        let s = r / (alpha / 2.0).tan();
+        let s = if round { r / (alpha / 2.0).tan() } else { r };
         let w = (alpha / 2.0).sin();
         // The cylinder, well past both ends (sections on oblique end faces reach beyond them).
         let m = 4.0 * (s + r);
         let (q0, q1) = (p0 - d * m, p1 + d * m);
         let ctrl = vec![vec![h(q0 + ta * s, 1.0), h(q1 + ta * s, 1.0)], vec![h(q0, w), h(q1, w)], vec![h(q0 + tb * s, 1.0), h(q1 + tb * s, 1.0)]];
-        let surf = mt::Surface::NurbsSurface(mt::NurbsSurface::new(mt::BSplineSurface::new(
-            (mt::KnotVec::bezier_knot(2), mt::KnotVec::bezier_knot(1)),
-            ctrl,
-        )));
+        let surf = if round {
+            mt::Surface::NurbsSurface(mt::NurbsSurface::new(mt::BSplineSurface::new(
+                (mt::KnotVec::bezier_knot(2), mt::KnotVec::bezier_knot(1)),
+                ctrl,
+            )))
+        } else {
+            mt::Surface::Plane(mt::Plane::new(p3(q0 + ta * s), p3(q1 + ta * s), p3(q0 + tb * s)))
+        };
         let mid = (p0 + p1) * 0.5;
         let axis = mid + ta * s - na * r;
         let out = ((ta + tb) * -1.0).normalized().unwrap_or(na);
@@ -311,20 +323,33 @@ fn build(cur: &Body, sel_pts: &[(Vec3, f64)]) -> Result<Body> {
             let rest: Vec<(mt::EdgeID, Vec3)> = ar.edges.iter().filter(|(id, _)| !sels.contains_key(id)).copied().collect();
             chosen.sort_by_key(|(id, _)| order.iter().position(|x| x == id));
             let sel = |id: mt::EdgeID| sels.get(&id).ok_or_else(|| KernelError::Failed("fillet".into()));
+            // At an inner rim a contact line can extend the next edge backwards.
+            // Its new endpoint must remain in the trimmed face, including its holes.
+            let check_extension = |q: Vec3, fi: usize| -> Result<Vec3> {
+                let face = faces.get(fi).ok_or_else(|| KernelError::Failed("face".into()))?;
+                if !crate::blend::extension_fits(face, v, q, size * 1e-7) {
+                    return Err(unsupported("a fillet larger than the faces next to it"));
+                }
+                Ok(q)
+            };
+            let hit = |u: Vec3, t: Vec3, s: f64, fi: usize| -> Result<Vec3> {
+                let q = hit(v, u, t, s)?;
+                if u.dot(t) < 0.0 { check_extension(q, fi) } else { Ok(q) }
+            };
             let list = pool.entry(*vid).or_default();
             match (&chosen[..], &rest[..]) {
                 ([(e1, _)], [(e2, u2), (e3, u3)]) => {
                     let s1 = sel(*e1)?;
                     let (f12, f31) = (face_with(*e1, *e2)?, face_with(*e1, *e3)?);
                     let f23 = face_with(*e2, *e3)?;
-                    let p2 = hit(v, *u2, s1.t(f12), s1.s)?;
-                    let p3_ = hit(v, *u3, s1.t(f31), s1.s)?;
+                    let p2 = hit(*u2, s1.t(f12), s1.s, f12)?;
+                    let p3_ = hit(*u3, s1.t(f31), s1.s, f31)?;
                     let (a, b) = (builder::vertex(p3(p2)), builder::vertex(p3(p3_)));
                     endpt.insert((*e2, *vid, None), a.clone());
                     endpt.insert((*e1, *vid, Some(f12)), a.clone());
                     endpt.insert((*e3, *vid, None), b.clone());
                     endpt.insert((*e1, *vid, Some(f31)), b.clone());
-                    let c = edge(&a, &b, mt::Curve::NurbsCurve(conic(p2, v, p3_, s1.w)))?;
+                    let c = if round { edge(&a, &b, mt::Curve::NurbsCurve(conic(p2, v, p3_, s1.w)))? } else { builder::line(&a, &b) };
                     list.push((c, vec![Owner::Face(f23), Owner::Cyl(*e1)]));
                 }
                 ([(ea, ua), (eb, ub)], [(e3, u3)]) => {
@@ -332,14 +357,19 @@ fn build(cur: &Body, sel_pts: &[(Vec3, f64)]) -> Result<Body> {
                     let fab = face_with(*ea, *eb)?;
                     let (fb3, fa3) = (face_with(*eb, *e3)?, face_with(*ea, *e3)?);
                     let q = meet(v, (*ua, sa.t(fab), sa.s), (*ub, sb.t(fab), sb.s))?;
+                    if ua.dot(sb.t(fab)) < 0.0 || ub.dot(sa.t(fab)) < 0.0 {
+                        check_extension(q, fab)?;
+                    }
                     // Where each cylinder touches the third edge.
-                    let xb = hit(v, *u3, sb.t(fb3), sb.s)?;
-                    let xa = hit(v, *u3, sa.t(fa3), sa.s)?;
+                    let xb = hit(*u3, sb.t(fb3), sb.s, fb3)?;
+                    let xa = hit(*u3, sa.t(fa3), sa.s, fa3)?;
                     let qv = builder::vertex(p3(q));
                     endpt.insert((*ea, *vid, Some(fab)), qv.clone());
                     endpt.insert((*eb, *vid, Some(fab)), qv.clone());
                     let (la, lb) = ((xa - v).len(), (xb - v).len());
-                    let crease_to = |end: Vec3| -> Result<mt::Curve> { crease(sa, sb, q, end, size) };
+                    let crease_to = |end: Vec3| -> Result<mt::Curve> {
+                        if round { crease(sa, sb, q, end, size) } else { Ok(mt::Curve::Line(mt::Line(p3(q), p3(end)))) }
+                    };
                     if (la - lb).abs() <= size * 1e-9 {
                         let xv = builder::vertex(p3(xb));
                         endpt.insert((*e3, *vid, None), xv.clone());
@@ -354,16 +384,22 @@ fn build(cur: &Body, sel_pts: &[(Vec3, f64)]) -> Result<Body> {
                             if lb > la { (sb, sa, *ua, fb3, fa3, *eb, *ea, xb) } else { (sa, sb, *ub, fa3, fb3, *ea, *eb, xa) };
                         // The big cylinder's section by the small one's end face: from the third
                         // edge to the small edge's line, where it crosses the small contact line.
-                        let z = hit(v, us, big.t(fab), big.s)?;
+                        let z = hit(us, big.t(fab), big.s, fab)?;
                         let mut sec = conic(x, v, z, big.w);
                         let tt = small.t(fsmall3);
                         let c0 = (x - v).dot(tt) - small.s;
                         let c1 = big.w * (0.0 - small.s);
                         let c2 = (z - v).dot(tt) - small.s;
-                        let t = bernstein_root(c0, c1, c2).ok_or_else(|| unsupported("fillets this different in size at one corner"))?;
-                        let y = {
+                        // A chamfer section is a line, not the cylinder's conic.
+                        // Where it meets the other strip's contact line is also
+                        // where the crease ends; the remainder trims that wall.
+                        let t = if round { bernstein_root(c0, c1, c2) } else { Some(c0 / (c0 - c2)).filter(|t| *t > 0.0 && *t < 1.0) }
+                            .ok_or_else(|| unsupported("blends this different in size at one corner"))?;
+                        let y = if round {
                             use mt::ParametricCurve;
                             crate::body::from_p3(sec.subs(t))
+                        } else {
+                            x + (z - x) * t
                         };
                         {
                             use mt::Cut;
@@ -374,13 +410,16 @@ fn build(cur: &Body, sel_pts: &[(Vec3, f64)]) -> Result<Body> {
                         endpt.insert((*e3, *vid, None), xv.clone());
                         endpt.insert((eb_, *vid, Some(fbig3)), xv.clone());
                         endpt.insert((es, *vid, Some(fsmall3)), yv.clone());
-                        let piece = edge(&xv, &yv, mt::Curve::NurbsCurve(sec))?;
+                        let piece = if round { edge(&xv, &yv, mt::Curve::NurbsCurve(sec))? } else { builder::line(&xv, &yv) };
                         list.push((piece, vec![Owner::Face(fsmall3), Owner::Cyl(eb_)]));
                         let c = edge(&qv, &yv, crease_to(y)?)?;
                         list.push((c, vec![Owner::Cyl(*ea), Owner::Cyl(*eb)]));
                     }
                 }
                 ([(e1, u1), (e2, u2), (e3, u3)], []) => {
+                    if !round {
+                        return Err(unsupported("three chamfers at one corner in the local blend"));
+                    }
                     let (s1, s2, s3) = (sel(*e1)?, sel(*e2)?, sel(*e3)?);
                     let (f12, f23, f31) = (face_with(*e1, *e2)?, face_with(*e2, *e3)?, face_with(*e3, *e1)?);
                     let q12 = meet(v, (*u1, s1.t(f12), s1.s), (*u2, s2.t(f12), s2.s))?;
@@ -1086,4 +1125,107 @@ mod tests {
             .unwrap();
         assert!(fillet_radii(&l, &[(Vec3::new(10.0, 10.0, 10.0), 1.0), (c, 2.0)]).is_err());
     }
+
+    #[test]
+    fn sloped_basin_rim_fillets() {
+        use solvecraft_geom::{Loop2, Plane, Vec2};
+        let rectangle = |x0, y0, x1, y1| Loop2::polygon(&[Vec2::new(x0, y0), Vec2::new(x1, y0), Vec2::new(x1, y1), Vec2::new(x0, y1)]);
+        let slab = box_solid(Vec3::ZERO, Vec3::new(100.0, 60.0, 30.0)).unwrap();
+        let pocket =
+            crate::loft(&[(Plane::XY.offset(30.0), rectangle(20.0, 4.0, 80.0, 50.0)), (Plane::XY.offset(10.0), rectangle(30.0, 20.0, 70.0, 40.0))])
+                .unwrap();
+        let basin = crate::boolean(&slab, &pocket, crate::BoolOp::Cut).unwrap().unwrap();
+        let picks = [Vec3::new(50.0, 4.0, 30.0), Vec3::new(80.0, 27.0, 30.0), Vec3::new(50.0, 50.0, 30.0), Vec3::new(20.0, 27.0, 30.0)];
+        for count in 1..=4 {
+            let result = fillet(&basin, &picks[..count], 2.0).unwrap();
+            assert!(result.validity().is_empty());
+            assert_eq!(result.unmeshed_faces(), 0);
+            assert!(vol(&result) < vol(&basin));
+            let mesh = result.tessellate(0.01).unwrap();
+            assert!(mesh.contains(Vec3::new(50.0, 2.0, 27.0)));
+            assert!(!mesh.contains(Vec3::new(50.0, 3.9, 29.99)));
+            // The small blend removes only a strip at the rim, and never escapes
+            // the slab. Topology checks alone do not catch crossed faces.
+            assert!(vol(&basin) - vol(&result) < 2.0_f64.powi(2) * (60.0 + 46.0) * 2.0);
+            for v in result.solid.vertex_iter() {
+                let p = vtx(&v);
+                assert!(p.x >= 0.0 && p.x <= 100.0 && p.y >= 0.0 && p.y <= 60.0 && p.z >= 0.0 && p.z <= 30.0 + 1e-8);
+            }
+        }
+    }
+
+    #[test]
+    fn sloped_rim_fillets_cannot_extend_outside_the_material() {
+        use solvecraft_geom::{Loop2, Plane, Vec2};
+        for margin in [0.5, 1.0] {
+            let rectangle = |x0, y0, x1, y1| Loop2::polygon(&[Vec2::new(x0, y0), Vec2::new(x1, y0), Vec2::new(x1, y1), Vec2::new(x0, y1)]);
+            let slab = box_solid(Vec3::ZERO, Vec3::new(100.0, 60.0, 30.0)).unwrap();
+            let pocket = crate::loft(&[
+                (Plane::XY.offset(30.0), rectangle(20.0, margin, 80.0, 50.0)),
+                (Plane::XY.offset(10.0), rectangle(30.0, 20.0, 70.0, 40.0)),
+            ])
+            .unwrap();
+            let basin = crate::boolean(&slab, &pocket, crate::BoolOp::Cut).unwrap().unwrap();
+            let picks = [Vec3::new(50.0, margin, 30.0), Vec3::new(80.0, 27.0, 30.0), Vec3::new(50.0, 50.0, 30.0), Vec3::new(20.0, 27.0, 30.0)];
+            for count in [1, 2, 4] {
+                assert!(fillet(&basin, &picks[..count], 4.0).is_err());
+                assert!(fillet(&basin, &picks[..count], margin * 0.5).is_ok());
+            }
+        }
+    }
+
+    #[test]
+    fn chamfers_meet_on_unequal_sloped_rim_walls() {
+        use solvecraft_geom::{Loop2, Plane, Vec2};
+        for margin in [1.0, 4.0] {
+            let rectangle = |x0, y0, x1, y1| Loop2::polygon(&[Vec2::new(x0, y0), Vec2::new(x1, y0), Vec2::new(x1, y1), Vec2::new(x0, y1)]);
+            let slab = box_solid(Vec3::ZERO, Vec3::new(100.0, 60.0, 30.0)).unwrap();
+            let pocket = crate::loft(&[
+                (Plane::XY.offset(30.0), rectangle(20.0, margin, 80.0, 50.0)),
+                (Plane::XY.offset(10.0), rectangle(30.0, 20.0, 70.0, 40.0)),
+            ])
+            .unwrap();
+            let basin = crate::boolean(&slab, &pocket, crate::BoolOp::Cut).unwrap().unwrap();
+            let picks = [Vec3::new(50.0, margin, 30.0), Vec3::new(80.0, 27.0, 30.0), Vec3::new(50.0, 50.0, 30.0), Vec3::new(20.0, 27.0, 30.0)];
+            for count in [2, 3, 4] {
+                if margin == 1.0 {
+                    assert!(crate::chamfer(&basin, &picks[..count], 2.0).is_err());
+                    assert!(crate::chamfer(&basin, &picks[..count], 1.0).is_err());
+                }
+                let distance = if margin == 1.0 { 0.5 } else { 2.0 };
+                let result = crate::chamfer(&basin, &picks[..count], distance).unwrap();
+                assert!(result.validity().is_empty());
+                assert_eq!(result.unmeshed_faces(), 0);
+                assert!(vol(&result) < vol(&basin));
+                assert!(vol(&basin) - vol(&result) < distance.powi(2) * (60.0 + 50.0 - margin));
+                let reverse: Vec<_> = picks[..count].iter().copied().rev().collect();
+                let reversed = crate::chamfer(&basin, &reverse, distance).unwrap();
+                assert!((vol(&result) - vol(&reversed)).abs() < 1e-6);
+                // These are exact planes, with every boundary on its owning face.
+                for f in result.solid.face_iter() {
+                    let Some(n) = plane_normal(f) else { continue };
+                    let origin = vtx(&f.vertex_iter().next().unwrap());
+                    for v in f.vertex_iter() {
+                        assert!((vtx(&v) - origin).dot(n).abs() < 1e-7);
+                    }
+                }
+            }
+        }
+    }
+
+    #[test]
+    fn equal_slope_rim_loops_cannot_escape_the_face() {
+        use solvecraft_geom::{Loop2, Plane, Vec2};
+        let rectangle = |x0, y0, x1, y1| Loop2::polygon(&[Vec2::new(x0, y0), Vec2::new(x1, y0), Vec2::new(x1, y1), Vec2::new(x0, y1)]);
+        let slab = box_solid(Vec3::ZERO, Vec3::new(100.0, 60.0, 30.0)).unwrap();
+        let pocket =
+            crate::loft(&[(Plane::XY.offset(30.0), rectangle(20.0, 1.0, 80.0, 50.0)), (Plane::XY.offset(10.0), rectangle(30.0, 11.0, 70.0, 40.0))])
+                .unwrap();
+        let basin = crate::boolean(&slab, &pocket, crate::BoolOp::Cut).unwrap().unwrap();
+        let picks = [Vec3::new(50.0, 1.0, 30.0), Vec3::new(80.0, 27.0, 30.0), Vec3::new(50.0, 50.0, 30.0), Vec3::new(20.0, 27.0, 30.0)];
+        assert!(crate::chamfer(&basin, &picks, 2.0).is_err());
+        assert!(fillet(&basin, &picks, 4.0).is_err());
+        assert!(crate::chamfer(&basin, &picks, 0.5).is_ok());
+        assert!(fillet(&basin, &picks, 0.5).is_ok());
+    }
 }

crates/kernel/src/loopblend.rs

@@ -270,7 +270,7 @@ fn build(
         if !(phi > 1e-3 && phi < std::f64::consts::PI - 1e-3) {
             return None;
         }
-        let sb = r / (phi / 2.0).tan();
+        let sb = if round { r / (phi / 2.0).tan() } else { r };
         Some((m, t, sb, r / (phi / 2.0).sin()))
     };
     // The straight seam between the walls at a sharp joint (direction away from the face).
@@ -328,6 +328,19 @@ fn build(
         }
         pb.push(p + e * (sb / h1));
     }
+    // Sharp joints between straight rim edges can extend beyond the old corner.
+    // Equal wall slopes reach this path instead of the corner builder, and must
+    // also reject an extension through the outer wall or another hole.
+    for (i, q) in pa.iter().enumerate() {
+        // Smooth arc transitions and open-run ends retain their existing checks.
+        if !mitre.get(i).copied().unwrap_or(false) {
+            continue;
+        }
+        let p = joint_point(i).ok_or_else(|| unsupported("degenerate loop joint"))?;
+        if !crate::blend::extension_fits(f, p, *q, size * 1e-7) {
+            return Err(unsupported("the blend is larger than the neighbouring faces"));
+        }
+    }
     // The offset arcs must not turn inside out; an arc of exactly the blend's radius shrinks
     // to a point (a rounded corner rounded again: the blend there is a sphere).
     let mut collapsed = vec![false; k.max(nj)];