Loft: pair section vertices by best fit (#881) (#956)

Each section now starts at the rotation whose vertices line up best with the lead section's, vertex for vertex, instead of at the one vertex closest in direction to the lead's first. Fixes the twisted loft of matching non-convex sections; regression test loft_pairs_the_vertices_of_matching_non_convex_sections.

crates/kernel/src/build.rs

@@ -1075,17 +1075,17 @@ pub(crate) fn loft_wires(sections: &[(Plane, Loop2)]) -> Result<(Vec<mt::Wire>,
         // Break angles: the lead section's vertex directions seen in this plane.
         let angles: Vec<f64> = lead_dirs.iter().map(|d| Vec2::new(d.dot(p.x), d.dot(p.y)).angle()).collect();
         let mut r = l.resample(n, &angles);
-        // Start at the vertex closest in direction to the lead's first vertex.
-        if let Some(d0) = lead_dirs.first() {
-            let target = Vec2::new(d0.dot(p.x), d0.dot(p.y)).angle();
-            let ang = |s: &Seg2| {
-                let d = ((s.start() - c).angle() - target).rem_euclid(std::f64::consts::TAU);
-                d.min(std::f64::consts::TAU - d)
-            };
-            let best = r.lp.segs.iter().enumerate().min_by(|a, b| ang(a.1).total_cmp(&ang(b.1))).map(|(k, _)| k);
-            if let Some(k) = best {
-                r.rotate_left(k);
-            }
+        // Start where every vertex lines up best with the lead's vertex of the same number (the
+        // largest sum of dot products of their directions from the centroids). The one vertex
+        // closest in direction to the lead's first vertex twisted non-convex sections (#881): a
+        // zig-zag's centroid sits off its outline, so its first vertex can point away from the
+        // lead's.
+        let dirs: Vec<Vec2> = lead_dirs.iter().map(|d| Vec2::new(d.dot(p.x), d.dot(p.y))).collect();
+        let m = r.lp.segs.len();
+        let fit =
+            |k: usize| -> f64 { dirs.iter().enumerate().filter_map(|(i, d)| r.lp.segs.get((i + k) % m).map(|s| (s.start() - c).dot(*d))).sum() };
+        if let Some(k) = (0..m).min_by(|a, b| fit(*b).total_cmp(&fit(*a))) {
+            r.rotate_left(k);
         }
         wires.push(uniform_wire(p, &r.lp)?);
         edge_sources.push(r.from);

crates/kernel/src/tests.rs

@@ -2326,3 +2326,62 @@ fn retained_bytes_matches_the_iterator_counts() {
         assert_eq!(body.retained_bytes(), walked);
     }
 }
+
+/// #881: a pleated band lofted to a nearly flat one, both 18 vertices from (-26, 3), counter-
+/// clockwise and mirror-symmetric about x = 0. The zig-zag's centroid sits well above its first
+/// vertex, so that vertex pointed nowhere near the lead's: the loft started three vertices round,
+/// twisted (14 820 mm³, centre of mass at x = 5.3). Vertex to vertex, it is a prismoid.
+#[test]
+fn loft_pairs_the_vertices_of_matching_non_convex_sections() {
+    let bottom = [
+        (-26.0, 3.0),
+        (-24.76384, 3.0),
+        (-16.0568, 14.97217),
+        (-8.0, 4.90117),
+        (0.0, 14.90117),
+        (8.0, 4.90117),
+        (16.0568, 14.97217),
+        (24.76384, 3.0),
+        (26.0, 3.0),
+        (26.0, 0.0),
+        (23.23616, 0.0),
+        (15.9432, 10.02783),
+        (8.0, 0.09883),
+        (0.0, 10.09883),
+        (-8.0, 0.09883),
+        (-15.9432, 10.02783),
+        (-23.23616, 0.0),
+        (-26.0, 0.0),
+    ];
+    let top = [
+        (-26.0, 3.0),
+        (-24.02812, 3.0),
+        (-16.0, 3.30105),
+        (-8.0, 3.00105),
+        (0.0, 3.30105),
+        (8.0, 3.00105),
+        (16.0, 3.30105),
+        (24.02812, 3.0),
+        (26.0, 3.0),
+        (26.0, 0.0),
+        (23.97188, 0.0),
+        (16.0, 0.29895),
+        (8.0, -0.00105),
+        (0.0, 0.29895),
+        (-8.0, -0.00105),
+        (-16.0, 0.29895),
+        (-23.97188, 0.0),
+        (-26.0, 0.0),
+    ];
+    let v = |p: &[(f64, f64)]| p.iter().map(|&(x, y)| Vec2::new(x, y)).collect::<Vec<_>>();
+    let mid: Vec<Vec2> = v(&bottom).into_iter().zip(v(&top)).map(|(a, b)| (a + b) * 0.5).collect();
+    // Prismoid: h/6 · (A0 + 4 Amid + A1) = 17 914.43.
+    let area = |p: &[Vec2]| Loop2::polygon(p).signed_area().abs();
+    let want = 90.0 / 6.0 * (area(&v(&bottom)) + 4.0 * area(&mid) + area(&v(&top)));
+    let secs = [(Plane::XY, Loop2::polygon(&v(&bottom))), (Plane::XY.offset(90.0), Loop2::polygon(&v(&top)))];
+    let m = measure(&loft(&secs).unwrap()).unwrap();
+    // Twisted ruled faces measure a little under the prismoid (17 884; a square lofted to a
+    // skewed quad, one twisted face, 0.3 % under), so the bounds leave room for that.
+    assert!(rel(m.volume, want) < 5e-3, "{} vs {want}", m.volume);
+    assert!(m.centroid.x.abs() < 0.05, "centre of mass x {}", m.centroid.x);
+}