3D and Materials: twisted walls and bevels are drawn back to front

At 180° and 360° of twist, Extrude & Bevel and Inflate drew some wall
strips in the wrong order, and a far strip showed through a near one.
Each strip ran the whole length of an outline edge and was sorted by
the depth of its center, while a twisted wall can cross it far from
that center. The caps were sorted the same way.

Twisted walls and bevels are now drawn in pieces about as long as a
wall slice is wide (fewer near the 32,768-facet limit), and twisted
bevel steps turn at most 6° each, as wall slices do. A cap, or
Inflate's flat back, that faces the viewer is drawn after the walls,
and edge-on pieces of a twisted wall keep their own depth. Pieces that
face the viewer are drawn where their strip as a whole faces away, so
a twisted solid seen from behind shows its walls' corners around the
back cap. Untwisted solids give the same facets as on main.

Over 80 twisted squares, the facet seen over a point lay up to 79.7 pt
behind a nearer facet on main and now lies at most 5.0 pt behind.
"Twisted" in 72 pt type, extruded with a 360° twist, takes 12,311
facets instead of 7,622 and 3.1 ms instead of 1.9 ms.

Docs: docs/three-d.md.

crates/three-d/src/solid.rs

@@ -211,6 +211,9 @@ pub struct Facet {
     pub contours: Vec<Vec<Point>>,
     pub depth: f64,
     pub brightness: f64,
+    /// The contours in view space, before projection, for tests that measure depths.
+    #[cfg(test)]
+    view: Vec<Vec<V3>>,
 }
 
 impl Facet {
@@ -276,6 +279,8 @@ impl Inflate {
 const MAX_SOLID_POINTS: usize = 2048;
 /// Outline points Rotate samples at most: it only turns the outline, so an outlined paragraph fits.
 const MAX_PLANE_POINTS: usize = 65536;
+/// Pieces a twisted wall's strip is cut into along one edge at most.
+const MAX_PIECES: usize = 64;
 
 /// How many straight pieces follow `c` within `tolerance` (Wang's bound on a cubic's distance from
 /// its chords), 1 to 32.
@@ -413,7 +418,8 @@ struct Camera {
     origin: Point,
     /// The eye's distance in front of the centre; `None` without perspective.
     eye: Option<f64>,
-    /// Deeper than any part of the object: faces drawn first go this far behind it.
+    /// Deeper than any part of the object: faces drawn first go this far behind it, and faces
+    /// drawn last go this far in front of it.
     far: f64,
     light: V3,
 }
@@ -456,9 +462,23 @@ struct Builder<'a> {
     facets: Vec<Facet>,
 }
 
+/// Where a face goes in the order the facets are drawn, back to front.
+#[derive(Clone, Copy, PartialEq, Eq)]
+enum Order {
+    /// The face is sorted by the depth of its center, except that a face seen almost edge-on goes
+    /// first, under everything.
+    Depth,
+    /// The face is sorted by the depth of its center even when it is seen almost edge-on, so that
+    /// a piece of a twisted wall at the solid's outline covers the walls behind it.
+    Twisted,
+    /// The face is drawn after every other face when it is turned toward the viewer. A flat end of
+    /// the solid has the whole solid behind its plane, so it covers whatever it overlaps.
+    Last,
+}
+
 impl Builder<'_> {
-    /// Add a face of view-space contours with outward normal `n` (view space).
-    fn face(&mut self, contours: Vec<Vec<V3>>, n: V3) -> Result<(), String> {
+    /// Add a face of view-space contours with outward normal `n` (view space), drawn in `order`.
+    fn face(&mut self, contours: Vec<Vec<V3>>, n: V3, order: Order) -> Result<(), String> {
         let count: usize = contours.iter().map(Vec::len).sum();
         if count < 3 {
             return Ok(());
@@ -467,9 +487,10 @@ impl Builder<'_> {
         let centre = V3(centre.0 / count as f64, centre.1 / count as f64, centre.2 / count as f64);
         let n = n.unit();
         let facing = self.camera.facing(centre, n);
-        // A closed solid hides the faces turned away. One seen almost exactly edge-on is kept but
-        // drawn first, under everything: on a curved surface rounding alone turns such faces
-        // either way, and dropping every other one would leave gaps.
+        // A closed solid hides the faces turned away. One seen almost exactly edge-on is kept and,
+        // unless it is a piece of a twisted wall, drawn first, under everything: on a curved
+        // surface rounding alone turns such faces either way, and dropping every other one would
+        // leave gaps.
         let (n, behind) = if facing > 0.0 {
             (n, false)
         } else if self.closed && facing > -0.02 {
@@ -479,6 +500,8 @@ impl Builder<'_> {
         } else {
             (V3(-n.0, -n.1, -n.2), false)
         };
+        #[cfg(test)]
+        let view = contours.clone();
         let contours: Vec<Vec<Point>> = contours.iter().map(|c| c.iter().map(|v| self.camera.project(*v)).collect()).collect();
         if contours.iter().flatten().any(|p| !p.x.is_finite() || !p.y.is_finite()) {
             return Err("The 3D projection is not finite".into());
@@ -486,31 +509,59 @@ impl Builder<'_> {
         if contours.iter().map(|c| signed_area(c).abs()).sum::<f64>() < 1e-10 {
             return Ok(());
         }
-        let depth = if behind { centre.2 - self.camera.far } else { centre.2 };
-        self.facets.push(Facet { contours, depth, brightness: self.camera.brightness(n) });
+        let depth = if behind && order != Order::Twisted {
+            centre.2 - self.camera.far
+        } else if order == Order::Last {
+            centre.2 + self.camera.far
+        } else {
+            centre.2
+        };
+        self.facets.push(Facet {
+            contours,
+            depth,
+            brightness: self.camera.brightness(n),
+            #[cfg(test)]
+            view,
+        });
         if self.facets.len() > MAX_FACES {
             return Err("The 3D surface is too complex; simplify the shape".into());
         }
         Ok(())
     }
 
-    /// The band between two partner rings (same length), as one quad per edge. `outward` is the
+    /// The band between two partner rings (same length), as one quad per edge, or for a twisted
+    /// wall with `cuts` (see [`cuts`]) as `cuts[i]` quads along edge `i`. `outward` is the
     /// direction (before rotation) the band faces away from the solid, for orienting normals.
-    fn band(&mut self, a: &[V3], b: &[V3], outward: impl Fn(usize) -> V3) -> Result<(), String> {
+    fn band(&mut self, a: &[V3], b: &[V3], cuts: &[usize], outward: impl Fn(usize) -> V3) -> Result<(), String> {
         let n = a.len().min(b.len());
+        // At 0 and 1, the corners themselves.
+        let lerp = |p: V3, q: V3, s: f64| {
+            if s <= 0.0 {
+                p
+            } else if s >= 1.0 {
+                q
+            } else {
+                V3(p.0 + (q.0 - p.0) * s, p.1 + (q.1 - p.1) * s, p.2 + (q.2 - p.2) * s)
+            }
+        };
+        let order = if cuts.is_empty() { Order::Depth } else { Order::Twisted };
         for i in 0..n {
             let j = (i + 1) % n;
-            let quad = vec![a[i], a[j], b[j], b[i]];
-            // A twist bends the strip between two slices into a saddle: its normal is the whole
-            // quad's average (Newell's), not one corner's triangle, which leans far to one side.
-            let mut normal = newell(std::slice::from_ref(&quad));
-            if normal.dot(normal) < 1e-18 {
-                continue;
+            let pieces = cuts.get(i).copied().unwrap_or(1).max(1);
+            for k in 0..pieces {
+                let (s0, s1) = (k as f64 / pieces as f64, (k + 1) as f64 / pieces as f64);
+                let quad = vec![lerp(a[i], a[j], s0), lerp(a[i], a[j], s1), lerp(b[i], b[j], s1), lerp(b[i], b[j], s0)];
+                // A twist bends the strip between two slices into a saddle: its normal is the whole
+                // quad's average (Newell's), not one corner's triangle, which leans far to one side.
+                let mut normal = newell(std::slice::from_ref(&quad));
+                if normal.dot(normal) < 1e-18 {
+                    continue;
+                }
+                if normal.dot(self.camera.rotate(outward(i))) < 0.0 {
+                    normal = V3(-normal.0, -normal.1, -normal.2);
+                }
+                self.face(vec![quad], normal, order)?;
             }
-            if normal.dot(self.camera.rotate(outward(i))) < 0.0 {
-                normal = V3(-normal.0, -normal.1, -normal.2);
-            }
-            self.face(vec![quad], normal)?;
         }
         Ok(())
     }
@@ -521,6 +572,38 @@ impl Builder<'_> {
     }
 }
 
+/// How many pieces each edge of a twisted wall is drawn in, from two neighboring slices of the
+/// wall, `a` and `b` (partner rings): as many as make each piece about as long along the edge as
+/// the band between the slices is wide there, 1 to [`MAX_PIECES`]. Twisted walls pass in front of
+/// one another, and a strip the length of its edge has no one depth where another crosses it.
+fn cuts(a: &[V3], b: &[V3]) -> Vec<usize> {
+    let n = a.len().min(b.len());
+    let length = |p: V3, q: V3| q.sub(p).dot(q.sub(p)).sqrt();
+    (0..n)
+        .map(|i| {
+            let j = (i + 1) % n;
+            let (along, across) = (length(a[i], a[j]).max(length(b[i], b[j])), length(a[i], b[i]).max(length(a[j], b[j])));
+            let count = (along / across.max(1e-9)).ceil();
+            if count.is_finite() { (count as usize).clamp(1, MAX_PIECES) } else { MAX_PIECES }
+        })
+        .collect()
+}
+
+/// Lowers the [`cuts`] of every ring's edges (one list per ring) so that `bands` bands of each
+/// make no more than `room` quads in all, or one per edge and band when that is more: the pieces
+/// past each edge's first share the room in proportion.
+fn fit(cuts: &mut [Vec<usize>], bands: usize, room: usize) {
+    let edges = cuts.iter().map(Vec::len).sum::<usize>().saturating_mul(bands);
+    let total = cuts.iter().flatten().sum::<usize>().saturating_mul(bands);
+    let room = room.max(edges);
+    if total > room {
+        let share = (room - edges) as f64 / (total - edges) as f64;
+        for c in cuts.iter_mut().flatten() {
+            *c = 1 + ((c.saturating_sub(1)) as f64 * share).floor() as usize;
+        }
+    }
+}
+
 /// The area-weighted normal of (possibly non-planar) contours, by Newell's method: for contours
 /// winding with their filled side on the left as seen from in front, it points to the front.
 fn newell(contours: &[Vec<V3>]) -> V3 {
@@ -622,8 +705,28 @@ fn steps(o: &Extrude) -> Vec<Step> {
         Some(b) if !b.inside => -b.width,
         _ => 0.0,
     };
+    // A twist turns a bevel's steps too: each step between two samples is cut into steps that
+    // turn no further than 6°, as the walls' slices do at most.
+    let refine = |pts: Vec<(f64, f64)>| -> Vec<(f64, f64)> {
+        if o.twist == 0.0 {
+            return pts;
+        }
+        let mut out = Vec::with_capacity(pts.len());
+        for (k, p) in pts.iter().enumerate() {
+            if let Some(q) = k.checked_sub(1).and_then(|k| pts.get(k)) {
+                let turn = (p.1 - q.1).abs() * height / d * o.twist.abs();
+                let n = ((turn / 6.0).ceil() as usize).clamp(1, 64);
+                for m in 1..n {
+                    let s = m as f64 / n as f64;
+                    out.push((q.0 + (p.0 - q.0) * s, q.1 + (p.1 - q.1) * s));
+                }
+            }
+            out.push(*p);
+        }
+        out
+    };
     if let Some(b) = bevel {
-        for (x, y) in profile(&b) {
+        for (x, y) in refine(profile(&b)) {
             let inset = if b.inside { b.width * (1.0 - x) } else { -b.width * x };
             out.push(Step { inset, t: y * height / d, face: 1.0 });
         }
@@ -639,7 +742,7 @@ fn steps(o: &Extrude) -> Vec<Step> {
         out.push(Step { inset: wall, t, face: 0.0 });
     }
     if let Some(b) = bevel.filter(|b| b.both_sides) {
-        for (x, y) in profile(&b).into_iter().rev() {
+        for (x, y) in refine(profile(&b)).into_iter().rev() {
             let inset = if b.inside { b.width * (1.0 - x) } else { -b.width * x };
             out.push(Step { inset, t: 1.0 - y * height / d, face: -1.0 });
         }
@@ -664,7 +767,22 @@ pub fn extrude(path: &PathData, rule: FillRule, bounds: Rect, options: Extrude,
     let z = |t: f64| o.depth / 2.0 - o.depth * t;
     let ring_at =
         |ring: &[Point], s: &Step| -> Vec<V3> { section(&inset(ring, s.inset), bounds, &o, s.t).iter().map(|p| camera.place(*p, z(s.t))).collect() };
-    for ring in &rings {
+    // Every band of a twisted outline, bevels included, is cut as its walls' first slice is (with
+    // bevels and no walls, as its shallowest band is), within the facet limit.
+    let depth_of = |s: &[Step]| s.get(1).zip(s.first()).map_or(0.0, |(b, a)| (b.t - a.t).abs());
+    let wall = steps
+        .windows(2)
+        .find(|s| s[0].face == 0.0 && s[1].face == 0.0)
+        .or_else(|| steps.windows(2).filter(|s| depth_of(s) > 1e-12).min_by(|a, c| depth_of(a).total_cmp(&depth_of(c))));
+    let mut pieces: Vec<Vec<usize>> = rings
+        .iter()
+        .map(|ring| match wall {
+            Some([s0, s1]) if o.twist != 0.0 => cuts(&ring_at(ring, s0), &ring_at(ring, s1)),
+            _ => vec![],
+        })
+        .collect();
+    fit(&mut pieces, steps.len().saturating_sub(1), MAX_FACES);
+    for (ring, pieces) in rings.iter().zip(&pieces) {
         let placed: Vec<Vec<V3>> = steps.iter().map(|s| ring_at(ring, s)).collect();
         for (k, pair) in placed.windows(2).enumerate() {
             let ([a, c], Some(s0), Some(s1)) = (pair, steps.get(k), steps.get(k + 1)) else { continue };
@@ -677,13 +795,16 @@ pub fn extrude(path: &PathData, rule: FillRule, bounds: Rect, options: Extrude,
             } else {
                 0.0
             };
-            b.band(a, c, |i| outward(ring, i, face))?;
+            b.band(a, c, pieces, |i| outward(ring, i, face))?;
         }
     }
     if o.cap {
+        // A cap turned toward the viewer covers the twisted walls it overlaps, even those nearer
+        // than its center.
+        let ends = if o.twist != 0.0 { Order::Last } else { Order::Depth };
         for (s, n) in [(steps.first(), 1.0), (steps.last(), -1.0)] {
             let Some(s) = s else { continue };
-            b.face(rings.iter().map(|r| ring_at(r, s)).collect(), camera.rotate(V3(0.0, 0.0, n)))?;
+            b.face(rings.iter().map(|r| ring_at(r, s)).collect(), camera.rotate(V3(0.0, 0.0, n)), ends)?;
         }
     }
     Ok(b.finish())
@@ -697,7 +818,7 @@ pub fn plane(path: &PathData, rule: FillRule, bounds: Rect, view: View) -> Resul
     let camera = Camera::new(view, bounds, 0.0);
     let mut b = Builder { camera: &camera, closed: false, facets: vec![] };
     let contours: Vec<Vec<V3>> = rings.iter().map(|r| r.iter().map(|p| camera.place(*p, 0.0)).collect()).collect();
-    b.face(contours, camera.rotate(V3(0.0, 0.0, 1.0)))?;
+    b.face(contours, camera.rotate(V3(0.0, 0.0, 1.0)), Order::Depth)?;
     Ok(b.finish())
 }
 
@@ -1175,7 +1296,7 @@ pub fn inflate(path: &PathData, rule: FillRule, bounds: Rect, options: Inflate,
                     let front: Vec<Vec<V3>> =
                         pieces.iter().zip(&lifts).map(|(p, h)| p.iter().zip(h).map(|(q, h)| camera.place(*q, half + h)).collect()).collect();
                     let n = newell(&front);
-                    b.face(front, n)?;
+                    b.face(front, n, Order::Depth)?;
                     if o.both_sides && o.cap {
                         // The back winds the same way round as the front, so its outward normal is
                         // reversed.
@@ -1185,7 +1306,7 @@ pub fn inflate(path: &PathData, rule: FillRule, bounds: Rect, options: Inflate,
                             .map(|(p, h)| p.iter().zip(h).map(|(q, h)| camera.place(back_turn.apply(*q), -half - h * back_turn.scale)).collect())
                             .collect();
                         let n = newell(&back);
-                        b.face(back, V3(-n.0, -n.1, -n.2))?;
+                        b.face(back, V3(-n.0, -n.1, -n.2), Order::Depth)?;
                     }
                 }
             }
@@ -1193,19 +1314,23 @@ pub fn inflate(path: &PathData, rule: FillRule, bounds: Rect, options: Inflate,
     }
     if !o.both_sides && o.cap {
         let back: Vec<Vec<V3>> = rings.iter().map(|r| r.iter().map(|p| camera.place(back_turn.apply(*p), -half)).collect()).collect();
-        b.face(back, camera.rotate(V3(0.0, 0.0, -1.0)))?;
+        // The back turned toward the viewer covers the twisted walls it overlaps, even those
+        // nearer than its center.
+        b.face(back, camera.rotate(V3(0.0, 0.0, -1.0)), if o.twist != 0.0 { Order::Last } else { Order::Depth })?;
     }
     if o.depth > 0.0 {
         // The sides, cut into slices that follow a twist or taper.
         let slices = if o.twist != 0.0 || o.taper != 1.0 { ((o.twist.abs() / 6.0).ceil() as usize).clamp(8, 64) } else { 1 };
-        for ring in &rings {
-            let at = |k: usize| -> Vec<V3> {
-                let t = k as f64 / slices as f64;
-                let turn = Turn::new(bounds, o.twist, o.taper, t);
-                ring.iter().map(|p| camera.place(turn.apply(*p), half - o.depth * t)).collect()
-            };
+        let at = |ring: &[Point], k: usize| -> Vec<V3> {
+            let t = k as f64 / slices as f64;
+            let turn = Turn::new(bounds, o.twist, o.taper, t);
+            ring.iter().map(|p| camera.place(turn.apply(*p), half - o.depth * t)).collect()
+        };
+        let mut pieces: Vec<Vec<usize>> = rings.iter().map(|r| if o.twist != 0.0 { cuts(&at(r, 0), &at(r, 1)) } else { vec![] }).collect();
+        fit(&mut pieces, slices, MAX_FACES.saturating_sub(b.facets.len()));
+        for (ring, pieces) in rings.iter().zip(&pieces) {
             for k in 0..slices {
-                b.band(&at(k), &at(k + 1), |i| outward(ring, i, 0.0))?;
+                b.band(&at(ring, k), &at(ring, k + 1), pieces, |i| outward(ring, i, 0.0))?;
             }
         }
     }

crates/three-d/src/tests_solid.rs

@@ -91,11 +91,15 @@ fn taper_shrinks_the_back_and_twist_turns_it() {
     let (flo, fhi) = at_front.iter().fold((f64::INFINITY, f64::NEG_INFINITY), |(lo, hi), q| (lo.min(q.y), hi.max(q.y)));
     let (back_h, front_h) = (hi - lo, fhi - flo);
     assert!(close(front_h, 100.0, 1e-6) && close(back_h, 50.0, 1e-6), "front {front_h}, back {back_h}");
-    // Seen from behind, a 45° twist shows the back cap turned into a diamond.
+    // Seen from behind, a 45° twist shows the back cap turned into a diamond, drawn last, and the
+    // corners of the walls nearer the front reaching past it.
     let back = View { rotation_y: 180.0, ..front() };
     let twisted = extrude(&p, FillRule::EvenOdd, b, Extrude { twist: 45.0, cap: true, ..o }, back).unwrap();
-    assert_eq!(twisted.len(), 1, "only the back cap faces the viewer");
+    let cap = twisted.last().unwrap();
+    assert_eq!(cap.contours[0].len(), 4);
+    assert!(close(extent(std::slice::from_ref(cap)).width(), 100.0 * 2f64.sqrt(), 1e-6), "{:?}", extent(&twisted));
     assert!(close(extent(&twisted).width(), 100.0 * 2f64.sqrt(), 1e-6), "{:?}", extent(&twisted));
+    assert!(depth_at(cap, Point::new(94.0, 94.0)).is_none() && twisted.iter().any(|f| depth_at(f, Point::new(94.0, 94.0)).is_some()));
     // Its walls are cut into slices that follow the turn.
     let open = extrude(&p, FillRule::EvenOdd, b, Extrude { twist: 45.0, ..o }, View { rotation_x: -20.0, rotation_y: 30.0, ..front() }).unwrap();
     assert!(open.len() > 2 * 8, "{}", open.len());
@@ -369,6 +373,16 @@ fn small_curves_take_few_points_and_rotate_takes_large_outlines() {
     assert!(ring.len() > 4 * 16, "a 100 pt circle within 0.05 pt: {}", ring.len());
 }
 
+/// The pieces twisted walls and bevels are cut into stay within the facet limit: 200 bumps with
+/// a rounded bevel, turned 360°, give no more facets than the limit.
+#[test]
+fn twisted_pieces_stay_within_the_facet_limit() {
+    let (p, r) = bumps(200);
+    let o = Extrude { twist: 360.0, bevel: Some(Bevel { shape: BevelShape::Rounded, ..Bevel::default() }), ..Extrude::default() };
+    let f = extrude(&p, FillRule::EvenOdd, r, o, View { rotation_x: -18.0, rotation_y: -26.0, rotation_z: 8.0, ..front() });
+    assert!(f.is_ok_and(|f| f.len() <= MAX_FACES));
+}
+
 /// The bucketed nearest-edge search agrees with checking every edge, inside, outside and far
 /// beyond the shape.
 #[test]
@@ -410,10 +424,19 @@ fn twist_holds_between_0_and_360_degrees() {
 
 /// A twisted extrusion shows every wall facing the viewer: a twist bends each strip between two
 /// slices into a saddle, and its normal from one corner's triangle leaned so far back that from
-/// 60° on every wall was dropped. Checked against each strip's average normal, oriented away from
-/// the square's axis.
+/// 60° on every wall was dropped. Checked against the average normal of each piece the strip is
+/// cut into along its edge ([`cuts`]), oriented away from the square's axis.
 #[test]
 fn twisted_walls_face_the_way_they_turn() {
+    let lerp = |p: V3, q: V3, s: f64| {
+        if s <= 0.0 {
+            p
+        } else if s >= 1.0 {
+            q
+        } else {
+            V3(p.0 + (q.0 - p.0) * s, p.1 + (q.1 - p.1) * s, p.2 + (q.2 - p.2) * s)
+        }
+    };
     let (p, b) = shape(&square(0.0, 0.0, 100.0));
     let c = b.center();
     for twist in [60.0, 120.0, 180.0, 300.0] {
@@ -427,6 +450,8 @@ fn twisted_walls_face_the_way_they_turn() {
             let ring = &rings(&p, FillRule::EvenOdd, b.width().max(b.height()), MAX_SOLID_POINTS).unwrap()[0];
             let slices = ((twist / 6.0).ceil() as usize).clamp(8, 64);
             let z = |t: f64| o.depth / 2.0 - o.depth * t;
+            let placed = |t: f64| -> Vec<V3> { section(ring, b, &o, t).iter().map(|p| camera.place(*p, z(t))).collect() };
+            let pieces = cuts(&placed(0.0), &placed(1.0 / slices as f64));
             let mut expected = 0;
             for k in 0..slices {
                 let (t0, t1) = (k as f64 / slices as f64, (k + 1) as f64 / slices as f64);
@@ -434,17 +459,24 @@ fn twisted_walls_face_the_way_they_turn() {
                 for i in 0..ring.len() {
                     let j = (i + 1) % ring.len();
                     let q = [camera.place(r0[i], z(t0)), camera.place(r0[j], z(t0)), camera.place(r1[j], z(t1)), camera.place(r1[i], z(t1))];
-                    let mut n = newell(&[q.to_vec()]).unit();
-                    let mid = r0[i].midpoint(r1[j]);
-                    let away = camera.rotate(V3(mid.x - c.x, mid.y - c.y, 0.0));
-                    if n.dot(away) < 0.0 {
-                        n = V3(-n.0, -n.1, -n.2);
-                    }
-                    let centre = V3((q[0].0 + q[2].0) / 2.0, (q[0].1 + q[2].1) / 2.0, (q[0].2 + q[2].2) / 2.0);
-                    // Kept as the builder keeps them: facing, or seen almost edge-on, with an area.
-                    let projected: Vec<Point> = q.iter().map(|v| camera.project(*v)).collect();
-                    if camera.facing(centre, n) > -0.02 && signed_area(&projected).abs() >= 1e-10 {
-                        expected += 1;
+                    let pieces = pieces[i];
+                    for k in 0..pieces {
+                        let (s0, s1) = (k as f64 / pieces as f64, (k + 1) as f64 / pieces as f64);
+                        let q = [lerp(q[0], q[1], s0), lerp(q[0], q[1], s1), lerp(q[3], q[2], s1), lerp(q[3], q[2], s0)];
+                        let mut n = newell(&[q.to_vec()]).unit();
+                        let s = (s0 + s1) / 2.0;
+                        let mid = r0[i].lerp(r0[j], s).midpoint(r1[i].lerp(r1[j], s));
+                        let away = camera.rotate(V3(mid.x - c.x, mid.y - c.y, 0.0));
+                        if n.dot(away) < 0.0 {
+                            n = V3(-n.0, -n.1, -n.2);
+                        }
+                        let sum = q.iter().fold(V3::default(), |s, v| V3(s.0 + v.0, s.1 + v.1, s.2 + v.2));
+                        let centre = V3(sum.0 / 4.0, sum.1 / 4.0, sum.2 / 4.0);
+                        // Kept as the builder keeps them: facing, or seen almost edge-on, with an area.
+                        let projected: Vec<Point> = q.iter().map(|v| camera.project(*v)).collect();
+                        if camera.facing(centre, n) > -0.02 && signed_area(&projected).abs() >= 1e-10 {
+                            expected += 1;
+                        }
                     }
                 }
             }
@@ -459,6 +491,142 @@ fn twisted_walls_face_the_way_they_turn() {
     }
 }
 
+/// The depth (larger is nearer the eye) of `f`'s own surface over the page point `p`, seen without
+/// perspective, or `None` where `f` doesn't cover `p`. A flat facet is its plane, a four-cornered
+/// one that isn't flat is the bilinear patch between its corners, and any other is a fan of
+/// triangles from its center.
+fn depth_at(f: &Facet, p: Point) -> Option<f64> {
+    let winding: i32 = f
+        .contours
+        .iter()
+        .flat_map(|c| c.iter().zip(c.iter().cycle().skip(1)).take(c.len()))
+        .map(|(a, b)| {
+            let side = (b.x - a.x) * (p.y - a.y) - (p.x - a.x) * (b.y - a.y);
+            if a.y <= p.y && b.y > p.y && side > 0.0 {
+                1
+            } else if a.y > p.y && b.y <= p.y && side < 0.0 {
+                -1
+            } else {
+                0
+            }
+        })
+        .sum();
+    if winding == 0 {
+        return None;
+    }
+    // Without perspective a projected point is its view-space point moved by one offset.
+    let pairs: Vec<Vec<(Point, f64)>> = f.contours.iter().zip(&f.view).map(|(c, v)| c.iter().zip(v).map(|(q, v)| (*q, v.2)).collect()).collect();
+    let n = newell(&f.view);
+    let (anchor, z0) = pairs[0][0];
+    let flat = n.2.abs() > 1e-12
+        && pairs.iter().flatten().all(|(q, z)| (n.0 * (q.x - anchor.x) + n.1 * (q.y - anchor.y) + n.2 * (z - z0)).abs() < 1e-9 * n.dot(n).sqrt());
+    if flat {
+        return Some(z0 - (n.0 * (p.x - anchor.x) + n.1 * (p.y - anchor.y)) / n.2);
+    }
+    let c = &pairs[0];
+    if pairs.len() == 1 && c.len() == 4 {
+        let at = |u: f64, v: f64| {
+            let w = [(1.0 - u) * (1.0 - v), u * (1.0 - v), u * v, (1.0 - u) * v];
+            c.iter().zip(w).fold((0.0, 0.0, 0.0), |s, ((q, z), w)| (s.0 + w * q.x, s.1 + w * q.y, s.2 + w * z))
+        };
+        // Newton's method for the patch point over `p`.
+        let (mut u, mut v) = (0.5, 0.5);
+        for _ in 0..50 {
+            let (x, y, _) = at(u, v);
+            let h = 1e-7;
+            let (xu, yu, _) = at(u + h, v);
+            let (xv, yv, _) = at(u, v + h);
+            let (a, b, cc, d) = ((xu - x) / h, (xv - x) / h, (yu - y) / h, (yv - y) / h);
+            let det = a * d - b * cc;
+            if det.abs() < 1e-18 {
+                return None;
+            }
+            let (ex, ey) = (p.x - x, p.y - y);
+            u += (d * ex - b * ey) / det;
+            v += (a * ey - cc * ex) / det;
+        }
+        return (u > -1e-6 && u < 1.0 + 1e-6 && v > -1e-6 && v < 1.0 + 1e-6).then(|| at(u, v).2);
+    }
+    let m = c.iter().fold((0.0, 0.0, 0.0), |s, (q, z)| (s.0 + q.x, s.1 + q.y, s.2 + z));
+    let k = c.len() as f64;
+    let (mx, my, mz) = (m.0 / k, m.1 / k, m.2 / k);
+    for ((a, za), (b, zb)) in c.iter().zip(c.iter().cycle().skip(1)).take(c.len()) {
+        let det = (a.x - mx) * (b.y - my) - (b.x - mx) * (a.y - my);
+        if det.abs() < 1e-18 {
+            continue;
+        }
+        let s = ((p.x - mx) * (b.y - my) - (b.x - mx) * (p.y - my)) / det;
+        let t = ((a.x - mx) * (p.y - my) - (p.x - mx) * (a.y - my)) / det;
+        if s >= -1e-9 && t >= -1e-9 && s + t <= 1.0 + 1e-9 {
+            return Some(mz + s * (za - mz) + t * (zb - mz));
+        }
+    }
+    None
+}
+
+/// Over a grid of `samples` × `samples` points across the facets, seen without perspective: how
+/// far, at most, the facet drawn last over a point (the one seen) is behind a facet drawn before
+/// it there, counting only facets it lies wholly behind the plane of (which must be drawn after
+/// it), and at which point.
+fn worst_misorder(facets: &[Facet], samples: usize) -> (f64, Point) {
+    let e = extent(facets);
+    let boxes: Vec<Rect> = facets.iter().map(|f| extent(std::slice::from_ref(f))).collect();
+    let behind = |f: &Facet, other: &Facet| {
+        // The other facet's plane through its corners' average, its normal toward the eye (+z).
+        let n = newell(&other.view);
+        let n = if n.2 < 0.0 { V3(-n.0, -n.1, -n.2) } else { n }.unit();
+        let count = other.view.iter().map(Vec::len).sum::<usize>() as f64;
+        let c = other.view.iter().flatten().fold(V3::default(), |s, v| V3(s.0 + v.0 / count, s.1 + v.1 / count, s.2 + v.2 / count));
+        n.2 > 1e-9 && f.view.iter().flatten().all(|v| v.sub(c).dot(n) < -1e-3)
+    };
+    let mut worst = (0.0, e.origin());
+    for j in 0..samples {
+        for i in 0..samples {
+            let p = Point::new(e.x0 + (i as f64 + 0.5) / samples as f64 * e.width(), e.y0 + (j as f64 + 0.5) / samples as f64 * e.height());
+            let covering: Vec<(&Facet, f64)> =
+                facets.iter().zip(&boxes).filter(|(_, r)| r.contains(p)).filter_map(|(f, _)| depth_at(f, p).map(|d| (f, d))).collect();
+            let Some(((seen, depth), under)) = covering.split_last() else { continue };
+            for (f, d) in under {
+                if d - depth > worst.0 && behind(seen, f) {
+                    worst = (d - depth, p);
+                }
+            }
+        }
+    }
+    worst
+}
+
+/// A twist of 180° or 360° turns a square's walls into a twisted bar, whose ridges pass in front
+/// of the walls between them, in front of its bevels and in front of its caps. Over every point
+/// the facet drawn last is less than 6 pt behind any facet drawn before it whose plane it lies
+/// wholly behind. Inflate's walls twist the same way (at 0 % volume its front is flat).
+#[test]
+fn twisted_walls_are_drawn_back_to_front() {
+    let (p, b) = shape(&square(0.0, 0.0, 100.0));
+    let rounded = Bevel { shape: BevelShape::Rounded, width: 8.0, height: 8.0, both_sides: true, ..Bevel::default() };
+    for (x, y, z) in [(-18.0, -26.0, 8.0), (-60.0, 20.0, 0.0), (75.0, 10.0, 30.0), (-20.0, 150.0, 10.0), (30.0, 200.0, 0.0)] {
+        let v = View { rotation_x: x, rotation_y: y, rotation_z: z, ..front() };
+        for twist in [180.0, 360.0] {
+            for depth in [50.0, 200.0] {
+                let solids = [
+                    ("Extrude", extrude(&p, FillRule::EvenOdd, b, Extrude { depth, twist, ..Extrude::default() }, v)),
+                    ("beveled", extrude(&p, FillRule::EvenOdd, b, Extrude { depth, twist, bevel: Some(rounded), ..Extrude::default() }, v)),
+                    ("Inflate", inflate(&p, FillRule::EvenOdd, b, Inflate { depth, twist, volume: 0.0, ..Inflate::default() }, v)),
+                    // 10 pt deep, the bevels take all of it; 40 pt deep, they leave a wall.
+                    (
+                        "bevels alone",
+                        extrude(&p, FillRule::EvenOdd, b, Extrude { depth: depth / 5.0, twist, bevel: Some(rounded), ..Extrude::default() }, v),
+                    ),
+                ];
+                for (kind, f) in solids {
+                    let (gap, at) = worst_misorder(&f.unwrap(), 120);
+                    assert!(gap < 6.0, "{kind}, view ({x}, {y}, {z}), twist {twist}, depth {depth}: {gap} pt behind at {at:?}");
+                }
+            }
+        }
+    }
+}
+
 /// Inflate takes Twist, Taper and Cap like the other types: the back turned and scaled (its puff
 /// too), the sides following, and without a cap no back at all.
 #[test]

docs/three-d.md

@@ -121,6 +121,9 @@ out, perspective). No Adobe software, assets, presets or bevel shapes were used.
 - The bevel shapes and their names are our own; custom bevel paths are not supported.
 - Visibility is a painter's algorithm: faces of a closed solid turned away are left out, the rest drawn back to
   front. Concave solids seen at steep angles can overlap wrongly; open (hollow) solids show both sides.
+  Twisted walls and bevels are drawn in pieces about as long as a wall slice is wide (fewer near the limit of
+  32,768 facets), and a cap or flat back turned toward the viewer is drawn over them. The piece drawn on top can
+  still lie a few points behind another nearby piece.
 - Shading is flat per facet: an inflated rim seen almost edge-on shows faint banding where the facets get finer.
 - Paths and compound paths only: convert type, images and groups to outlines first. One 3D effect per object.
 - Curves are followed within 1/2000 of the shape's size (at least 0.05 pt), so a small curve takes a few points