fix: fix polygonize dropping body area for polygons with touching holes

.gitignore

@@ -1,2 +1,3 @@
 /target
-/.vscode
+/.vscode
+.idea/

fixtures/polygonizer/input/taiwan_bug_minimal.geojson

@@ -0,0 +1,108 @@
+{
+  "type": "FeatureCollection",
+  "features": [
+    {
+      "type": "Feature",
+      "geometry": {
+        "type": "Polygon",
+        "coordinates": [
+          [
+            [
+              80.16730278417222,
+              11.513844408273066
+            ],
+            [
+              127.68336218282334,
+              3.464753784417475
+            ],
+            [
+              138.98936479970567,
+              47.93827024865087
+            ],
+            [
+              46.77538031026475,
+              30.091086782693232
+            ],
+            [
+              80.16730278417222,
+              11.513844408273066
+            ]
+          ],
+          [
+            [
+              112.20145815297715,
+              22.908672250985468
+            ],
+            [
+              79.86003225728624,
+              15.676307358025872
+            ],
+            [
+              71.10004442617051,
+              21.868254102944697
+            ],
+            [
+              112.20145815297715,
+              22.908672250985468
+            ]
+          ]
+        ]
+      },
+      "properties": null
+    },
+    {
+      "type": "Feature",
+      "geometry": {
+        "type": "Polygon",
+        "coordinates": [
+          [
+            [
+              108.0560302734375,
+              21.533888876052067
+            ],
+            [
+              122.62939453125,
+              41.409775793166325
+            ],
+            [
+              108.10546875,
+              24.846565305800603
+            ],
+            [
+              108.0560302734375,
+              21.533888876052067
+            ]
+          ]
+        ]
+      },
+      "properties": null
+    },
+    {
+      "type": "Feature",
+      "geometry": {
+        "type": "Polygon",
+        "coordinates": [
+          [
+            [
+              124.30344892254274,
+              41.20182861410828
+            ],
+            [
+              117.32219020774764,
+              41.20182861410828
+            ],
+            [
+              117.32219020774764,
+              36.76893152856083
+            ],
+            [
+              124.30344892254274,
+              41.20182861410828
+            ]
+          ]
+        ]
+      },
+      "properties": null
+    }
+  ]
+}

src/graph/topology_cleanup.rs

@@ -369,7 +369,23 @@ fn split_touching_boundary_polygon<T: GeoFloat + rstar::RTreeNum>(
     let deduplicated_faces = deduplicate_faces_by_exterior(face_candidates);
     let split_faces = prune_container_faces(&deduplicated_faces);
 
-    if split_faces.len() >= 2 {
+    // Guard against dropping body area on holed polygons: if the polygon
+    // has holes, the body (exterior minus holes) is multiply-connected
+    // and generally not representable as minimal edge rings. The
+    // extracted sub-faces may only cover "pockets" between touching
+    // holes, silently discarding the bulk of the body. Require that a
+    // representative body point remain covered by the split faces.
+    let body_interior_covered = polygon.interiors().is_empty()
+        || polygon
+            .interior_point()
+            .map(|interior_point| {
+                split_faces
+                    .iter()
+                    .any(|face| face.contains(&interior_point))
+            })
+            .unwrap_or(false);
+
+    if split_faces.len() >= 2 && body_interior_covered {
         split_faces
     } else if polygon.interiors().is_empty() {
         split_no_hole_polygon_on_repeated_vertex(&polygon).unwrap_or_else(|| vec![polygon])

src/tests.rs

@@ -138,6 +138,52 @@ fn load_input_lines(name: &str) -> Vec<Line<f64>> {
         .collect()
 }
 
+fn load_input_polygons(name: &str) -> Vec<Polygon<f64>> {
+    let collection = load_feature_collection("input", name);
+    let mut polygons = Vec::new();
+    for feature in collection.features {
+        let geometry = feature
+            .geometry
+            .unwrap_or_else(|| panic!("input fixture {name} has a feature with no geometry"));
+        match geometry.value {
+            GeometryValue::Polygon { .. } => {
+                let polygon: Polygon<f64> = geometry.try_into().unwrap();
+                polygons.push(polygon);
+            }
+            GeometryValue::MultiPolygon { .. } => {
+                let multi: MultiPolygon<f64> = geometry.try_into().unwrap();
+                polygons.extend(multi.0);
+            }
+            other => {
+                panic!("input fixture {name} uses unsupported polygon geometry type: {other:?}")
+            }
+        }
+    }
+    polygons
+}
+
+fn multipolygon_to_geojson_string(multi: &MultiPolygon<f64>) -> String {
+    let features = multi
+        .0
+        .iter()
+        .map(|polygon| geojson::Feature {
+            bbox: None,
+            geometry: Some(geojson::Geometry::new(GeometryValue::from(polygon))),
+            id: None,
+            properties: None,
+            foreign_members: None,
+        })
+        .collect();
+
+    let collection = geojson::FeatureCollection {
+        bbox: None,
+        features,
+        foreign_members: None,
+    };
+
+    serde_json::to_string_pretty(&collection).expect("serialize multipolygon to geojson")
+}
+
 fn load_expected_polygons_from(output_kind: &str, name: &str) -> MultiPolygon<f64> {
     let collection = load_feature_collection(output_kind, name);
     let mut polygons = Vec::new();
@@ -923,6 +969,82 @@ fn polygonize_nested_shell_overlap_minimal_matches_fixture() {
     assert_polygonize_fixture("nested_shell_overlap_minimal");
 }
 
+// ---------------------------------------------------------------------------
+// Large real-world bug reproducers
+// ---------------------------------------------------------------------------
+
+/// Coverage preservation through `nodify(1e-6) + polygonize`: any point
+/// that lies in the interior of some input polygon must still lie in the
+/// interior of exactly one output polygon.
+///
+/// The fixture is a 3-feature, 17-vertex reduction of a real-world
+/// MultiPolygon whose holes touch each other inside a single shell. On
+/// that shape, `split_touching_boundary_polygons` used to replace the
+/// shell with sub-face "pocket" rings extracted between touching holes,
+/// discarding the bulk of the body; this regression test pins the
+/// invariant at probe point (121.0, 23.5), which sits in the body.
+///
+/// Depends on pass 2 (split_touching).
+#[test]
+fn polygonize_preserves_enclosed_taiwan_polygon_minimal_matches_fixture() {
+    let taiwan_probe = point! { x: 121.0_f64, y: 23.5_f64 };
+
+    let input_polygons = load_input_polygons("taiwan_bug_minimal");
+    let input_owners: Vec<usize> = input_polygons
+        .iter()
+        .enumerate()
+        .filter_map(|(index, polygon)| {
+            if polygon.contains(&taiwan_probe) {
+                Some(index)
+            } else {
+                None
+            }
+        })
+        .collect();
+    assert!(
+        !input_owners.is_empty(),
+        "pre-condition: Taiwan probe (121.0, 23.5) must be contained by at least one \
+         input polygon in the taiwan_bug_minimal fixture, got none"
+    );
+
+    let input_lines: Vec<_> = input_polygons
+        .iter()
+        .flat_map(|polygon| polygon.lines_iter())
+        .collect();
+    let noded_lines = nodify_lines(input_lines.into_iter(), 1e-6_f64);
+    let polygons = polygonize(noded_lines.into_iter());
+
+    eprintln!(
+        "polygonize output for taiwan_bug_minimal ({} polygons):\n{}",
+        polygons.0.len(),
+        multipolygon_to_geojson_string(&polygons)
+    );
+
+    let containing_polygons: Vec<usize> = polygons
+        .0
+        .iter()
+        .enumerate()
+        .filter_map(|(index, polygon)| {
+            if polygon.contains(&taiwan_probe) {
+                Some(index)
+            } else {
+                None
+            }
+        })
+        .collect();
+
+    assert_eq!(
+        containing_polygons.len(),
+        1,
+        "Taiwan probe (121.0, 23.5) is contained by {} input polygon(s) (indices {:?}) \
+         but by {} output polygon(s) (indices {:?})",
+        input_owners.len(),
+        input_owners,
+        containing_polygons.len(),
+        containing_polygons
+    );
+}
+
 // ===========================================================================
 // Standalone nodify fixture tests
 // ===========================================================================