cdt: fix constraint sign restored onto the wrong edge in Contour::try_clip (flood erase leaks parity)

cdt/src/contour.rs

@@ -161,7 +161,7 @@ impl Contour {
                 ContourData::None => (),
                 ContourData::Hull(hull_index, sign) => {
                     t.hull.update(hull_index, e_ca);
-                    t.half.set_sign(e_bc, sign);
+                    t.half.set_sign(e_ca, sign);
                 },
                 ContourData::Buddy(b) => t.half.link_new(b, e_ca),
             };

cdt/tests/concave_erase.rs

@@ -0,0 +1,54 @@
+// A deeply concave (but simple -- no self-intersections) heptagon whose constrained
+// triangulation must be the 5-triangle interior partition. On the unfixed code the
+// flood erase removes nothing and the full 6-triangle convex-hull triangulation is
+// returned instead: constraint (3,4) sits on the partial hull when constraint (2,3)'s
+// cavity walk consumes it, and `Contour::try_clip`'s positive-branch `a` arm restores
+// the preserved sign onto the wrong rebuilt edge (`e_bc` instead of `e_ca`), so the
+// constraint comes back unfixed and the odd-even erase parity leaks through it.
+#[test]
+fn concave_heptagon_erases_its_pocket() {
+    let pts: Vec<(f64, f64)> = vec![
+        (41.3, 9.1),
+        (22.1, 63.1),
+        (14.2, 64.4),
+        (33.6, 8.0),
+        (1.6, 12.8),
+        (0.2, 5.9),
+        (40.3, 0.1),
+    ];
+    let edges: Vec<(usize, usize)> = (0..7).map(|i| (i, (i + 1) % 7)).collect();
+    let tris = cdt::triangulate_with_edges(&pts, &edges).expect("triangulation");
+
+    // n-2 triangles for a simple n-gon's interior partition.
+    assert_eq!(
+        tris.len(),
+        5,
+        "expected the interior partition, got {:?}",
+        tris
+    );
+
+    // Unsigned triangle area must equal the polygon's net shoelace area (a partition --
+    // the convex-hull triangulation covers ~2.2x this polygon's net area instead).
+    let tri_area: f64 = tris
+        .iter()
+        .map(|&(a, b, c)| {
+            let (ax, ay) = pts[a];
+            let (bx, by) = pts[b];
+            let (cx, cy) = pts[c];
+            ((bx - ax) * (cy - ay) - (cx - ax) * (by - ay)).abs() * 0.5
+        })
+        .sum();
+    let mut net = 0.0f64;
+    for i in 0..pts.len() {
+        let (x0, y0) = pts[i];
+        let (x1, y1) = pts[(i + 1) % pts.len()];
+        net += x0 * y1 - x1 * y0;
+    }
+    let net = (net * 0.5).abs();
+    assert!(
+        (tri_area - net).abs() <= net * 1e-9,
+        "covered area {} != polygon net area {}",
+        tri_area,
+        net
+    );
+}