Historia-Urbis/road_parceling/tests/degenerate.rs

//! Degenerate-case suite from spec Table `tab:degenerate` (§5).
//!
//! Every named test exists. Tests that require milestone-0.2
//! features (full deformation pipeline, curved-road support, sliver
//! merging) are gated with `#[ignore]` and a clear TODO so they
//! show up red in `cargo test -- --ignored` and green in the default
//! run.

use glam::DVec2;
use road_parceling::{
    apply_road_edit, subdivide_all, EdgeKind, RoadEdit, RoadGraph, SubdivisionError,
    SubdivisionParams,
};

// ---------- helpers ----------

fn rectangle_graph(w: f64, h: f64) -> RoadGraph {
    let mut g = RoadGraph::new();
    let a = g.add_node(DVec2::new(0.0, 0.0));
    let b = g.add_node(DVec2::new(w, 0.0));
    let c = g.add_node(DVec2::new(w, h));
    let d = g.add_node(DVec2::new(0.0, h));
    g.add_road(a, b).unwrap();
    g.add_road(b, c).unwrap();
    g.add_road(c, d).unwrap();
    g.add_road(d, a).unwrap();
    g.rebuild_topology().unwrap();
    g
}

fn assert_invariants_i1_i3(parcels: &road_parceling::ParcelSet) {
    // I1 — every parcel polygon is simple, CCW, non-degenerate. The
    // `Polygon` constructor enforces this on insertion, so iterating
    // and re-constructing each polygon is sufficient.
    for (_, p) in parcels.iter() {
        let v: Vec<DVec2> = p.vertices().to_vec();
        let _ = road_parceling::geometry::Polygon::new(v).expect("I1: polygon must validate");
    }
    // I2 — each parcel has exactly one Frontage edge.
    for (_, p) in parcels.iter() {
        let n_frontage = p
            .edge_kinds()
            .iter()
            .filter(|k| matches!(k, EdgeKind::Frontage))
            .count();
        assert_eq!(n_frontage, 1, "I2 violation: {n_frontage} frontages");
    }
    // I3 — non-overlap. Coarse cross-parcel check is deferred to
    // milestone 0.2 (requires polygon-polygon intersection); the
    // per-block area sum check lives in subdivide.rs unit tests.
}

// ----------------------------------------------------------------------
//  Tests #1 – #3 (acute / colinear)
// ----------------------------------------------------------------------

#[test]
fn acute_intersection_15deg() {
    // Two roads sharing a node meet at 15°. With the bisector-clip
    // at acute corners (milestone 0.3), no proper sliver-merge yet,
    // but I1–I3 must hold.
    let mut g = RoadGraph::new();
    let apex = g.add_node(DVec2::new(0.0, 0.0));
    let p1 = g.add_node(DVec2::new(100.0, 0.0));
    let angle = 15.0_f64.to_radians();
    let p2 = g.add_node(DVec2::new(100.0 * angle.cos(), 100.0 * angle.sin()));
    g.add_road(apex, p1).unwrap();
    g.add_road(apex, p2).unwrap();
    g.add_road(p1, p2).unwrap();
    g.rebuild_topology().unwrap();
    let params = SubdivisionParams::default();
    let parcels = subdivide_all(&g, &params).unwrap();
    // We don't insist parcels exist (a 15° wedge may be too narrow
    // for any parcel to satisfy `min_area`), but if any do, they
    // must satisfy I1–I3.
    assert_invariants_i1_i3(&parcels);
}

#[test]
fn acute_intersection_5deg() {
    // Knife-edge 5° angle. Library must not panic.
    let mut g = RoadGraph::new();
    let apex = g.add_node(DVec2::new(0.0, 0.0));
    let p1 = g.add_node(DVec2::new(100.0, 0.0));
    let angle = 5.0_f64.to_radians();
    let p2 = g.add_node(DVec2::new(100.0 * angle.cos(), 100.0 * angle.sin()));
    g.add_road(apex, p1).unwrap();
    g.add_road(apex, p2).unwrap();
    g.add_road(p1, p2).unwrap();
    g.rebuild_topology().unwrap();
    let params = SubdivisionParams::default();
    // Either succeeds with valid parcels or returns a typed error;
    // both are acceptable, but it must not panic.
    if let Ok(parcels) = subdivide_all(&g, &params) {
        assert_invariants_i1_i3(&parcels);
    }
}

#[test]
fn colinear_roads() {
    // Two end-to-end segments with zero turn at the shared node.
    // Currently we treat them as two short frontages; the spec
    // ultimately wants them merged into one continuous frontage,
    // but invariants must hold either way (milestone 0.1 baseline).
    let mut g = RoadGraph::new();
    let a = g.add_node(DVec2::new(0.0, 0.0));
    let b = g.add_node(DVec2::new(100.0, 0.0));
    let c = g.add_node(DVec2::new(200.0, 0.0));
    let d = g.add_node(DVec2::new(200.0, 100.0));
    let e = g.add_node(DVec2::new(0.0, 100.0));
    g.add_road(a, b).unwrap();
    g.add_road(b, c).unwrap(); // colinear with a-b
    g.add_road(c, d).unwrap();
    g.add_road(d, e).unwrap();
    g.add_road(e, a).unwrap();
    g.rebuild_topology().unwrap();
    let params = SubdivisionParams::default();
    let parcels = subdivide_all(&g, &params).unwrap();
    assert!(
        !parcels.is_empty(),
        "should subdivide despite colinear edges"
    );
    assert_invariants_i1_i3(&parcels);
}

// ----------------------------------------------------------------------
//  Tests #4 – #6 (graph degeneracies)
// ----------------------------------------------------------------------

#[test]
fn zero_length_segment() {
    // Two coincident-position nodes, road between them.
    let mut g = RoadGraph::new();
    let a = g.add_node(DVec2::new(0.0, 0.0));
    let b = g.add_node(DVec2::new(0.0, 0.0));
    let _ = g.add_road(a, b); // currently accepted (same id check only)
    let params = SubdivisionParams::default();
    // Either typed error or skip — both acceptable.
    let _ = subdivide_all(&g, &params);
}

#[test]
fn near_duplicate_nodes() {
    let mut g = RoadGraph::new();
    let a = g.add_node(DVec2::new(0.0, 0.0));
    let b = g.add_node(DVec2::new(1e-9, 0.0)); // within EPS_GEOM
    let c = g.add_node(DVec2::new(1.0, 0.0));
    let _ = g.add_road(a, b);
    let _ = g.add_road(b, c);
    let params = SubdivisionParams::default();
    let _ = subdivide_all(&g, &params); // must not panic
}

#[test]
fn self_intersecting_graph() {
    // Two roads cross with no node at the crossing.
    let mut g = RoadGraph::new();
    let a = g.add_node(DVec2::new(0.0, 0.0));
    let b = g.add_node(DVec2::new(2.0, 2.0));
    let c = g.add_node(DVec2::new(0.0, 2.0));
    let d = g.add_node(DVec2::new(2.0, 0.0));
    g.add_road(a, b).unwrap();
    g.add_road(c, d).unwrap();
    let err = g.rebuild_topology().unwrap_err();
    assert!(matches!(err, SubdivisionError::NonPlanarGraph(_)));
}

// ----------------------------------------------------------------------
//  Tests #7 – #9 (intersection topologies)
// ----------------------------------------------------------------------

#[test]
fn cul_de_sac() {
    // A road approaches a polygonal bulb (12-segment approximation
    // of a circle). The bulb's interior face becomes a block; its
    // boundary is the circle. Parcels in the bulb interior are
    // small rectangles facing each segment — not yet true pie
    // slices (that requires straight-skeleton subdivision,
    // milestone 0.4) — but I1–I3 must hold.
    use std::f64::consts::TAU;
    let mut g = RoadGraph::new();
    // Approach road: from (0, -200) up to (0, 0).
    let approach_start = g.add_node(DVec2::new(0.0, -200.0));
    let bulb_radius = 60.0;
    let bulb_segments = 12;
    let mut bulb_nodes = Vec::new();
    for i in 0..bulb_segments {
        let a = (i as f64 / bulb_segments as f64) * TAU;
        let p = DVec2::new(bulb_radius * a.cos(), bulb_radius * a.sin());
        bulb_nodes.push(g.add_node(p));
    }
    // Approach connects to the bulb at the bottom node.
    let bottom_idx = bulb_segments * 3 / 4; // angle 270°
    g.add_road(approach_start, bulb_nodes[bottom_idx]).unwrap();
    for i in 0..bulb_segments {
        let next = (i + 1) % bulb_segments;
        g.add_road(bulb_nodes[i], bulb_nodes[next]).unwrap();
    }
    g.rebuild_topology().unwrap();
    let params = SubdivisionParams::default();
    let parcels = subdivide_all(&g, &params).unwrap();
    // Bulb interior gets parcels along its inside; we don't insist
    // on a specific count, just that invariants hold.
    assert_invariants_i1_i3(&parcels);
}

#[test]
fn t_intersection() {
    //  D ---- F           y=h  ---->
    //  |      |
    //  A      |     vertical road b--c bisecting top edge
    //  |      |
    //  C ---- E           y=0
    //  Three blocks: left, right, top.
    //  But actually a real T has two cells. Build: bottom edge plus
    //  a vertical road meeting it.
    //  Layout:
    //    H───G───F
    //    │   │   │
    //    │   M   │
    //    │   │   │
    //    A───B───C
    //  with a vertical road B-M-G (M is intermediate). Two blocks.
    let mut g = RoadGraph::new();
    let a = g.add_node(DVec2::new(0.0, 0.0));
    let b = g.add_node(DVec2::new(100.0, 0.0));
    let c = g.add_node(DVec2::new(200.0, 0.0));
    let d = g.add_node(DVec2::new(200.0, 100.0));
    let e = g.add_node(DVec2::new(100.0, 100.0));
    let f = g.add_node(DVec2::new(0.0, 100.0));
    g.add_road(a, b).unwrap();
    g.add_road(b, c).unwrap();
    g.add_road(c, d).unwrap();
    g.add_road(d, e).unwrap();
    g.add_road(e, f).unwrap();
    g.add_road(f, a).unwrap();
    g.add_road(b, e).unwrap(); // T stem
    g.rebuild_topology().unwrap();
    let params = SubdivisionParams::default();
    let parcels = subdivide_all(&g, &params).unwrap();
    assert!(parcels.len() > 4);
    assert_invariants_i1_i3(&parcels);
}

#[test]
fn y_intersection() {
    // Three roads at 120°, meeting at the origin, with a triangular
    // outer boundary for the blocks to be bounded.
    let mut g = RoadGraph::new();
    let center = g.add_node(DVec2::new(0.0, 0.0));
    let r = 100.0;
    let third = std::f64::consts::TAU / 3.0;
    let p1 = g.add_node(DVec2::new(r * 0.0_f64.cos(), r * 0.0_f64.sin()));
    let p2 = g.add_node(DVec2::new(r * third.cos(), r * third.sin()));
    let p3 = g.add_node(DVec2::new(r * (2.0 * third).cos(), r * (2.0 * third).sin()));
    g.add_road(center, p1).unwrap();
    g.add_road(center, p2).unwrap();
    g.add_road(center, p3).unwrap();
    g.add_road(p1, p2).unwrap();
    g.add_road(p2, p3).unwrap();
    g.add_road(p3, p1).unwrap();
    g.rebuild_topology().unwrap();
    let params = SubdivisionParams::default();
    let parcels = subdivide_all(&g, &params).unwrap();
    assert!(!parcels.is_empty());
    assert_invariants_i1_i3(&parcels);
}

// ----------------------------------------------------------------------
//  Tests #10 – #12 (size extremes / curvature)
// ----------------------------------------------------------------------

#[test]
fn tiny_block() {
    // Perimeter < 4 * w_min. Should produce 0 parcels without
    // panicking.
    let g = rectangle_graph(5.0, 4.0);
    let params = SubdivisionParams::default(); // min_frontage = 6
    let parcels = subdivide_all(&g, &params).unwrap();
    assert!(parcels.len() <= 1);
    assert_invariants_i1_i3(&parcels);
}

#[test]
fn huge_block() {
    let g = rectangle_graph(1000.0, 1000.0);
    let params = SubdivisionParams::default();
    let parcels = subdivide_all(&g, &params).unwrap();
    // Sanity: with frontage_width=20 and perimeter=4000, expect
    // hundreds of parcels — but bounded.
    assert!(parcels.len() > 50);
    assert!(parcels.len() < 4000);
    assert_invariants_i1_i3(&parcels);
}

#[test]
fn curved_road_high_curv() {
    // A polyline approximating a tight curve: radius 20m, less than
    // the 30m default depth. Per-edge depth caps (computed by ray
    // cast across the block from each edge midpoint) should clamp
    // depth so parcels don't self-intersect through the centerline.
    use std::f64::consts::TAU;
    let mut g = RoadGraph::new();
    let radius = 20.0_f64; // tight: < params.depth=30
    let segments = 16;
    let mut inner = Vec::new();
    for i in 0..segments {
        let a = (i as f64 / segments as f64) * TAU;
        let p = DVec2::new(radius * a.cos(), radius * a.sin());
        inner.push(g.add_node(p));
    }
    let outer_radius = 100.0;
    let mut outer = Vec::new();
    for i in 0..segments {
        let a = (i as f64 / segments as f64) * TAU;
        let p = DVec2::new(outer_radius * a.cos(), outer_radius * a.sin());
        outer.push(g.add_node(p));
    }
    for i in 0..segments {
        let n = (i + 1) % segments;
        g.add_road(inner[i], inner[n]).unwrap();
        g.add_road(outer[i], outer[n]).unwrap();
    }
    // Spoke from inner[0] to outer[0] to bound the annulus into a
    // single block.
    g.add_road(inner[0], outer[0]).unwrap();
    g.rebuild_topology().unwrap();
    let params = SubdivisionParams::default();
    let parcels = subdivide_all(&g, &params).unwrap();
    // Library must not panic on tight curvature; invariants must
    // hold. (The inner ring of parcels has tighter depth caps from
    // the ray-cast per-edge midpoint computation.)
    assert_invariants_i1_i3(&parcels);
}

// ----------------------------------------------------------------------
//  Tests #13 – #18 (edits)
// ----------------------------------------------------------------------

#[test]
fn road_edit_micro_move() {
    // A 0.01 m move; with a true deform pipeline every parcel
    // would survive as `Deformed`. Milestone 0.1 regenerates, but
    // the post-edit set must still be valid.
    let mut g = rectangle_graph(200.0, 100.0);
    let params = SubdivisionParams::default();
    let mut parcels = subdivide_all(&g, &params).unwrap();
    let some_node = g.road_endpoints().next().unwrap().1;
    let pos = g.node_position(some_node).unwrap();
    let edit = RoadEdit::MoveNode {
        node: some_node,
        to: pos + DVec2::new(0.01, 0.0),
    };
    let _outcome = apply_road_edit(&mut parcels, &mut g, edit, &params).unwrap();
    assert!(!parcels.is_empty(), "post-edit set is non-empty");
    assert_invariants_i1_i3(&parcels);
}

#[test]
fn road_edit_large_move() {
    let mut g = rectangle_graph(200.0, 100.0);
    let params = SubdivisionParams::default();
    let mut parcels = subdivide_all(&g, &params).unwrap();
    let some_node = g.road_endpoints().next().unwrap().1;
    let pos = g.node_position(some_node).unwrap();
    let edit = RoadEdit::MoveNode {
        node: some_node,
        to: pos + DVec2::new(50.0, 0.0),
    };
    let _outcome = apply_road_edit(&mut parcels, &mut g, edit, &params).unwrap();
    assert_invariants_i1_i3(&parcels);
}

#[test]
fn road_edit_inverse_restores() {
    // Apply edit then its inverse — surviving parcels stay valid and
    // drift by no more than the magnitude of the edit delta. The
    // minimum-change deformation path (D14) skips parcels whose
    // frontage was unaffected by the line change, so strict
    // vertex-by-vertex restoration is not the contract; centroid-
    // bounded drift is.
    let mut g = rectangle_graph(200.0, 100.0);
    let params = SubdivisionParams::default();
    let mut parcels = subdivide_all(&g, &params).unwrap();

    let initial: std::collections::HashMap<_, DVec2> = parcels
        .iter()
        .map(|(id, p)| (id, p.polygon().centroid()))
        .collect();

    let some_node = g.road_endpoints().next().unwrap().1;
    let pos = g.node_position(some_node).unwrap();
    let delta = DVec2::new(0.5, 0.0);
    let _ = apply_road_edit(
        &mut parcels,
        &mut g,
        RoadEdit::MoveNode {
            node: some_node,
            to: pos + delta,
        },
        &params,
    )
    .unwrap();
    let _ = apply_road_edit(
        &mut parcels,
        &mut g,
        RoadEdit::MoveNode {
            node: some_node,
            to: pos,
        },
        &params,
    )
    .unwrap();

    for (id, p) in parcels.iter() {
        let Some(orig) = initial.get(&id) else {
            continue;
        };
        let drift = (p.polygon().centroid() - *orig).length();
        assert!(
            drift <= delta.length() + 1e-6,
            "centroid drift {drift} > edit delta {} for parcel {id:?}",
            delta.length()
        );
    }
    assert_invariants_i1_i3(&parcels);
}

#[test]
fn road_delete_condemns() {
    let mut g = rectangle_graph(200.0, 100.0);
    let params = SubdivisionParams::default();
    let mut parcels = subdivide_all(&g, &params).unwrap();
    let road = g.road_endpoints().next().unwrap().0;
    let n_on_road = parcels.parcels_on_road(road).count();
    let outcome = apply_road_edit(
        &mut parcels,
        &mut g,
        RoadEdit::DeleteSegment { road },
        &params,
    )
    .unwrap();
    assert!(!outcome.condemned.is_empty());
    assert!(outcome.condemned.len() >= n_on_road);
}

#[test]
fn road_split_preserves() {
    // Splitting a segment should preserve the parcels on it. Each
    // pre-split parcel either stays (frontage rebinds to one of the
    // two new roads) or is split into two parcels along the
    // perpendicular through the split point.
    let mut g = rectangle_graph(200.0, 100.0);
    let params = SubdivisionParams::default();
    let mut parcels = subdivide_all(&g, &params).unwrap();
    let pre_count = parcels.len();

    // Pick the bottom road (longest, has many parcels) and split at
    // x = 100, y = 0.
    let bottom_road = g
        .road_endpoints()
        .find(|&(_, _, b)| {
            let pos = g.node_position(b).unwrap();
            pos.y.abs() < 1e-6 && pos.x > 100.0
        })
        .map(|(rid, _, _)| rid)
        .unwrap();

    let outcome = apply_road_edit(
        &mut parcels,
        &mut g,
        RoadEdit::SplitSegment {
            road: bottom_road,
            at: DVec2::new(100.0, 0.0),
        },
        &params,
    )
    .unwrap();
    // Most parcels should be deformed (frontage_road rebound, no
    // geometric change). At most one is condemned (the one whose
    // frontage spans the split point — and it gets two replacements).
    assert!(
        outcome.regenerated.is_empty(),
        "split path must not trigger block regenerate; got {} regenerated",
        outcome.regenerated.len()
    );
    assert!(
        !outcome.deformed.is_empty(),
        "split path should preserve at least some parcels"
    );
    // Net parcel count should be ≥ pre_count (a split can ADD
    // parcels but never lose them outright).
    assert!(parcels.len() >= pre_count - 1);
    assert_invariants_i1_i3(&parcels);
}

#[test]
fn building_footprint_persists() {
    // A stub building's `fits_in` returns true iff the parcel area
    // stays above a threshold. After a tiny edit it stays above, so
    // the building survives.
    use road_parceling::{BuildingFitCheck, BuildingHandle, Parcel};

    struct StubBuilding {
        min_area: f64,
    }
    impl BuildingFitCheck for StubBuilding {
        fn fits_in(&self, parcel: &Parcel) -> bool {
            parcel.area() >= self.min_area
        }
    }

    let mut g = rectangle_graph(200.0, 100.0);
    let params = SubdivisionParams::default();
    let mut parcels = subdivide_all(&g, &params).unwrap();
    let pid = parcels.iter().map(|(id, _)| id).next().unwrap();
    let attach_area = parcels.get(pid).unwrap().area();
    parcels
        .get_mut(pid)
        .unwrap()
        .attach_building(BuildingHandle::new(StubBuilding {
            min_area: attach_area * 0.5,
        }));
    assert!(parcels.get(pid).unwrap().has_building());
    let some_node = g.road_endpoints().next().unwrap().1;
    let pos = g.node_position(some_node).unwrap();
    let _ = apply_road_edit(
        &mut parcels,
        &mut g,
        RoadEdit::MoveNode {
            node: some_node,
            to: pos + DVec2::new(0.1, 0.0),
        },
        &params,
    )
    .unwrap();
    if let Some(p) = parcels.get(pid) {
        assert!(p.has_building(), "building must persist on tiny edit");
    }
}

// ----------------------------------------------------------------------
//  Tests #19 – #21 (graph oddities)
// ----------------------------------------------------------------------

#[test]
fn degenerate_isolated_node() {
    let mut g = RoadGraph::new();
    let _isolated = g.add_node(DVec2::new(50.0, 50.0));
    // Plus a normal block for there to be something to subdivide.
    let a = g.add_node(DVec2::new(0.0, 0.0));
    let b = g.add_node(DVec2::new(100.0, 0.0));
    let c = g.add_node(DVec2::new(100.0, 100.0));
    let d = g.add_node(DVec2::new(0.0, 100.0));
    g.add_road(a, b).unwrap();
    g.add_road(b, c).unwrap();
    g.add_road(c, d).unwrap();
    g.add_road(d, a).unwrap();
    g.rebuild_topology().unwrap();
    let params = SubdivisionParams::default();
    let parcels = subdivide_all(&g, &params).unwrap();
    assert!(!parcels.is_empty());
    assert_invariants_i1_i3(&parcels);
}

#[test]
fn disconnected_graph() {
    // Two disjoint rectangles. Each must subdivide independently.
    let mut g = RoadGraph::new();
    // Component 1
    let a = g.add_node(DVec2::new(0.0, 0.0));
    let b = g.add_node(DVec2::new(100.0, 0.0));
    let c = g.add_node(DVec2::new(100.0, 100.0));
    let d = g.add_node(DVec2::new(0.0, 100.0));
    g.add_road(a, b).unwrap();
    g.add_road(b, c).unwrap();
    g.add_road(c, d).unwrap();
    g.add_road(d, a).unwrap();
    // Component 2
    let e = g.add_node(DVec2::new(500.0, 500.0));
    let f = g.add_node(DVec2::new(600.0, 500.0));
    let h = g.add_node(DVec2::new(600.0, 600.0));
    let i = g.add_node(DVec2::new(500.0, 600.0));
    g.add_road(e, f).unwrap();
    g.add_road(f, h).unwrap();
    g.add_road(h, i).unwrap();
    g.add_road(i, e).unwrap();
    g.rebuild_topology().unwrap();
    let params = SubdivisionParams::default();
    let parcels = subdivide_all(&g, &params).unwrap();
    assert!(!parcels.is_empty());
    // Roughly twice as many as a single block.
    let single = subdivide_all(&rectangle_graph(100.0, 100.0), &params).unwrap();
    assert!(parcels.len() >= 2 * single.len() - 4);
    assert_invariants_i1_i3(&parcels);
}

#[test]
fn numerical_precision_stress() {
    // Coords near 1e6 (well inside f64 mantissa). Spec says 1e20 but
    // anything past ~1e15 loses meter-level precision; we test a
    // city-scale-plus offset that f64 can still handle within
    // tolerances.
    let offset = 1.0e6;
    let mut g = RoadGraph::new();
    let a = g.add_node(DVec2::new(offset, offset));
    let b = g.add_node(DVec2::new(offset + 200.0, offset));
    let c = g.add_node(DVec2::new(offset + 200.0, offset + 100.0));
    let d = g.add_node(DVec2::new(offset, offset + 100.0));
    g.add_road(a, b).unwrap();
    g.add_road(b, c).unwrap();
    g.add_road(c, d).unwrap();
    g.add_road(d, a).unwrap();
    g.rebuild_topology().unwrap();
    let params = SubdivisionParams::default();
    let parcels = subdivide_all(&g, &params).unwrap();
    assert!(!parcels.is_empty());
    assert_invariants_i1_i3(&parcels);
}

#[test]
fn y_intersection_no_overlaps() {
    // Programmatic I3 check: no parcel's centroid is contained inside
    // any *other* parcel's polygon.
    let mut g = RoadGraph::new();
    let center = g.add_node(DVec2::new(0.0, 0.0));
    let r = 100.0;
    let third = std::f64::consts::TAU / 3.0;
    let p1 = g.add_node(DVec2::new(r * 0.0_f64.cos(), r * 0.0_f64.sin()));
    let p2 = g.add_node(DVec2::new(r * third.cos(), r * third.sin()));
    let p3 = g.add_node(DVec2::new(r * (2.0 * third).cos(), r * (2.0 * third).sin()));
    g.add_road(center, p1).unwrap();
    g.add_road(center, p2).unwrap();
    g.add_road(center, p3).unwrap();
    g.add_road(p1, p2).unwrap();
    g.add_road(p2, p3).unwrap();
    g.add_road(p3, p1).unwrap();
    g.rebuild_topology().unwrap();
    let params = SubdivisionParams::default();
    let parcels = subdivide_all(&g, &params).unwrap();

    let parcels_vec: Vec<_> = parcels.iter().collect();
    for (i, (_, pi)) in parcels_vec.iter().enumerate() {
        let centroid_i = pi.polygon().centroid();
        for (j, (_, pj_pair)) in parcels_vec.iter().enumerate() {
            if i == j {
                continue;
            }
            let pj = pj_pair.polygon();
            assert!(
                !pj.contains(centroid_i),
                "I3 violation: centroid of parcel {} ({:?}) is inside parcel {}",
                i,
                centroid_i,
                j
            );
        }
    }
}