//! Frontage-first block subdivision (spec §3.1, milestone 0.2).
//!
//! For each block face we:
//!
//! 1. Identify "real" corners — block-boundary vertices whose
//!    underlying graph node has degree ≥3, or degree 2 with a bend
//!    angle below 150°. Smooth continuations (collinear-roads case)
//!    are not real corners.
//! 2. At each real corner, build a corner parcel: a parallelogram
//!    extending `frontage_width` along each adjacent block edge and
//!    closing diagonally inward. The longer adjacent edge becomes the
//!    parcel's frontage; the other road-adjacent edge is reclassified
//!    as `Side`. (See §11 D9–D11.)
//! 3. For each block boundary edge, walk the *remaining* frontage
//!    (between the corner footprints, if any) at arc-length intervals
//!    of approximately `params.frontage_width`, extruding perpendicular
//!    to the block interior. No more bisector clipping — the corners
//!    take care of corner geometry up front.
//!
//! Determinism (invariant **I6**) is preserved: random jitter is
//! drawn from a `ChaCha8Rng` whose seed is mixed from `params.seed`
//! and the road-segment id; corner parcels are pure functions of
//! geometry.

use std::time::{Duration, Instant};

use glam::DVec2;
use rand::{Rng, SeedableRng};
use rand_chacha::ChaCha8Rng;
use slotmap::Key;

use crate::config::SubdivisionParams;
use crate::error::SubdivisionError;
use crate::geometry::{Polygon, EPS_GEOM};
use crate::network::blocks::{extract_blocks, Block};
use crate::network::graph::FaceId;
use crate::network::{NodeId, RoadGraph, RoadId};

use super::classify::classify_edges;
use super::regularize::regularize_parcel;
use super::{Parcel, ParcelSet};

/// Per-phase timing collected by [`subdivide_all_with_stats`].
///
/// Phases are non-overlapping wall-clock spans. `total` is measured
/// independently and may be slightly larger than the sum of the
/// phases due to overhead between phases.
#[derive(Debug, Clone, Default)]
pub struct SubdivisionStats {
    /// Number of bounded faces (blocks) extracted from the graph.
    pub block_count: usize,
    /// Total parcels produced across all blocks.
    pub parcel_count: usize,
    /// Time to (re)build the half-edge / DCEL topology if needed.
    pub topology_rebuild: Duration,
    /// Time to extract bounded faces into [`Block`] objects.
    pub block_extraction: Duration,
    /// Cumulative time spent in [`subdivide_block`] across all blocks.
    pub block_subdivide_total: Duration,
    /// Wall-clock total of [`subdivide_all_with_stats`].
    pub total: Duration,
}

impl SubdivisionStats {
    /// Average parcels generated per second over the whole call.
    #[must_use]
    pub fn parcels_per_second(&self) -> f64 {
        let secs = self.total.as_secs_f64();
        if secs > 0.0 {
            self.parcel_count as f64 / secs
        } else {
            0.0
        }
    }

    /// Average wall-clock time per parcel, in microseconds.
    #[must_use]
    pub fn time_per_parcel_us(&self) -> f64 {
        if self.parcel_count > 0 {
            (self.total.as_nanos() as f64 / self.parcel_count as f64) / 1_000.0
        } else {
            0.0
        }
    }
}

/// Subdivide every bounded face of `graph` into parcels.
///
/// # Errors
///
/// Returns a [`SubdivisionError`] if `params` is invalid, the graph
/// is non-planar, or an internal geometric failure occurs.
pub fn subdivide_all(
    graph: &RoadGraph,
    params: &SubdivisionParams,
) -> Result<ParcelSet, SubdivisionError> {
    subdivide_all_with_stats(graph, params).map(|(p, _)| p)
}

/// Subdivide every bounded face and return per-phase performance
/// statistics alongside the parcel set.
///
/// # Errors
///
/// Same as [`subdivide_all`].
pub fn subdivide_all_with_stats(
    graph: &RoadGraph,
    params: &SubdivisionParams,
) -> Result<(ParcelSet, SubdivisionStats), SubdivisionError> {
    params.validate()?;
    let mut stats = SubdivisionStats::default();
    let total_start = Instant::now();

    let topo_start = Instant::now();
    let mut working: RoadGraph;
    let g_ref = if graph.topology_valid {
        graph
    } else {
        working = graph.clone();
        working.rebuild_topology()?;
        &working
    };
    stats.topology_rebuild = topo_start.elapsed();

    let blocks_start = Instant::now();
    let blocks = extract_blocks(g_ref);
    stats.block_extraction = blocks_start.elapsed();
    stats.block_count = blocks.len();

    let mut parcels = ParcelSet::default();
    for (block_idx, block) in blocks.iter().enumerate() {
        let block_start = Instant::now();
        let new = subdivide_block(g_ref, block, params, block_idx)?;
        stats.block_subdivide_total += block_start.elapsed();
        for parcel in new {
            parcels.insert(parcel);
        }
    }
    stats.parcel_count = parcels.len();
    stats.total = total_start.elapsed();
    Ok((parcels, stats))
}

/// Subdivide a single block. Internal entry used both by
/// `subdivide_all` and by `apply_road_edit`'s regenerate path.
pub(crate) fn subdivide_block(
    graph: &RoadGraph,
    block: &Block,
    params: &SubdivisionParams,
    _block_idx: usize,
) -> Result<Vec<Parcel>, SubdivisionError> {
    let n = block.polygon.len();
    if n < 3 {
        return Ok(Vec::new());
    }
    let verts: Vec<DVec2> = block.polygon.vertices().to_vec();
    let face = block.face_id();

    // Edge lengths (parallel to verts).
    let edge_lens: Vec<f64> = (0..n)
        .map(|i| (verts[(i + 1) % n] - verts[i]).length())
        .collect();

    // Real-corner classification per vertex (D11). Acute corners
    // (interior angle < 60°) are skipped here — sliver-merge is
    // milestone 0.3 — so parcels there fall back to whatever the
    // regular walk produces (usually a malformed shape that triggers
    // I1 rejection during validation, but never a panic).
    // Per-vertex interior angle.
    let interior_angles: Vec<f64> = (0..n)
        .map(|i| {
            let v_prev = block.polygon.vertices()[(i + n - 1) % n];
            let v_curr = block.polygon.vertices()[i];
            let v_next = block.polygon.vertices()[(i + 1) % n];
            let t_in = (v_prev - v_curr).normalize_or_zero();
            let t_out = (v_next - v_curr).normalize_or_zero();
            let t_arrive = -t_in;
            let cross = t_arrive.x * t_out.y - t_arrive.y * t_out.x;
            let dot = t_arrive.dot(t_out);
            let turn = cross.atan2(dot);
            std::f64::consts::PI - turn
        })
        .collect();

    let real_corner: Vec<bool> = (0..n)
        .map(|i| {
            is_real_corner(graph, block, i, params) && interior_angles[i] > 60.0_f64.to_radians()
        })
        .collect();
    // Acute corners: real corners (degree-3 or sharp degree-2) with
    // interior < 60°. We don't build a parcel here, but adjacent
    // regulars get bisector-clipped so their territories stay
    // separated.
    let acute_corner: Vec<bool> = (0..n)
        .map(|i| {
            is_real_corner(graph, block, i, params) && interior_angles[i] <= 60.0_f64.to_radians()
        })
        .collect();

    // Per-edge depth cap (block half-depth bound).
    let depth_caps = edge_depth_caps(&block.polygon);

    let r_target = params.frontage_width;

    // First pass: decide which adjacent road wins frontage at each
    // real corner. This drives both the corner geometry below and the
    // per-edge walk bounds.
    let mut frontage_on_next: Vec<bool> = vec![false; n];
    for i in 0..n {
        if !real_corner[i] {
            continue;
        }
        let prev_len = edge_lens[(i + n - 1) % n];
        let next_len = edge_lens[i];
        let prev_road = block.roads[(i + n - 1) % n];
        let next_road = block.roads[i];
        frontage_on_next[i] = if (next_len - prev_len).abs() < EPS_GEOM {
            // tie → deterministic by RoadId
            next_road.data().as_ffi() <= prev_road.data().as_ffi()
        } else {
            next_len > prev_len
        };
    }

    let mut out: Vec<Parcel> = Vec::new();

    // Build corner parcels.
    for i in 0..n {
        if !real_corner[i] {
            continue;
        }
        let v_curr = verts[i];
        let v_prev = verts[(i + n - 1) % n];
        let v_next = verts[(i + 1) % n];

        let prev_road = block.roads[(i + n - 1) % n];
        let next_road = block.roads[i];

        let t_in = (v_prev - v_curr).normalize_or_zero();
        let t_out = (v_next - v_curr).normalize_or_zero();
        if t_in.length_squared() < EPS_GEOM * EPS_GEOM
            || t_out.length_squared() < EPS_GEOM * EPS_GEOM
        {
            continue;
        }

        let prev_len = edge_lens[(i + n - 1) % n];
        let next_len = edge_lens[i];
        let r = r_target
            .min(prev_len * 0.5)
            .min(next_len * 0.5)
            .max(params.min_frontage * 0.5);
        if r < params.min_frontage * 0.5 {
            continue;
        }

        // Inward normals.
        let n_out = DVec2::new(-t_out.y, t_out.x);
        let n_in = DVec2::new(t_in.y, -t_in.x);

        // Two flavors of corner parcel, both with frontage = R on the
        // *winning* road and side2 = depth along the other:
        //
        //   A. frontage on NEXT road: v0 → v0+t_out·R is frontage,
        //      v0+t_in·depth → v0 is the side along the prev road.
        //   B. frontage on PREV road: v0+t_in·R → v0 is frontage,
        //      v0 → v0+t_out·depth is the side along the next road.
        let (frontage_road, polygon, frontage_idx, consume_next, consume_prev) =
            if frontage_on_next[i] {
                let perp_depth = params.depth.min(prev_len * 0.5);
                let v0 = v_curr;
                let v1 = v_curr + t_out * r;
                let v2 = v_curr + t_out * r + n_out * perp_depth;
                let v3 = v_curr + t_in * perp_depth;
                let Ok(p) = Polygon::new(vec![v0, v1, v2, v3]) else {
                    continue;
                };
                (next_road, p, 0_usize, r, perp_depth)
            } else {
                let perp_depth = params.depth.min(next_len * 0.5);
                let v0 = v_curr;
                let v1 = v_curr + t_out * perp_depth;
                let v2 = v_curr + t_out * perp_depth + n_out * r;
                let v3 = v_curr + t_in * r;
                let Ok(p) = Polygon::new(vec![v0, v1, v2, v3]) else {
                    continue;
                };
                // Frontage edge in CCW order is v3 → v0, i.e. edge 3.
                (prev_road, p, 3_usize, perp_depth, r)
            };
        // The variables `consume_next` and `consume_prev` propagate
        // out via the per-vertex helpers below.
        let _ = (consume_next, consume_prev, n_in);

        // Clip to the block boundary so the parcel can never extend
        // past it. Required for non-rectangular blocks where the
        // R×depth rectangle would otherwise cross another block edge.
        let Some(polygon) = clip_polygon_to_block(&polygon, &block.polygon) else {
            continue;
        };
        if polygon.area() < params.min_area {
            continue;
        }
        let edge_kinds = classify_edges(polygon.len(), frontage_idx);

        let mut parcel = Parcel {
            polygon,
            edge_kinds,
            vertex_ids: Vec::new(),
            frontage_road,
            frontage_edge_index: frontage_idx,
            block: face,
            building: None,
        };
        regularize_parcel(&mut parcel, params);
        out.push(parcel);
    }

    // For each edge, walk between the corner footprints (if any) and
    // emit regular frontage parcels.
    for i in 0..n {
        let p = verts[i];
        let q = verts[(i + 1) % n];
        let road = block.roads[i];
        let edge_vec = q - p;
        let edge_len = edge_vec.length();
        if edge_len < params.min_frontage {
            continue;
        }
        let edge_dir = edge_vec / edge_len;
        let inward = DVec2::new(-edge_dir.y, edge_dir.x);

        // Walk-bounds. The corner at vertex i (start of this edge)
        // consumes the corner's `r` along this edge (it's the corner's
        // *next* side). The corner at vertex i+1 (end of this edge)
        // consumes the corner's `perp_depth` along this edge (it's
        // *that* corner's *prev* side, i.e., its Side2 runs along
        // here for a depth-equivalent distance).
        let consume_at_start = if real_corner[i] {
            // Edge i is the *next* edge of corner at vertex i.
            // - frontage_on_next[i] true  → consume R (frontage).
            // - frontage_on_next[i] false → consume depth (side2).
            if frontage_on_next[i] {
                corner_consume_short(&edge_lens, i, r_target, params)
            } else {
                corner_consume_long(&edge_lens, i, params)
            }
        } else {
            0.0
        };
        let next_v = (i + 1) % n;
        let consume_at_end = if real_corner[next_v] {
            // Edge i is the *prev* edge of corner at vertex (i+1).
            // - frontage_on_next[next_v] true  → consume depth (side2 along prev).
            // - frontage_on_next[next_v] false → consume R (frontage on prev).
            if frontage_on_next[next_v] {
                corner_consume_long(&edge_lens, next_v, params)
            } else {
                corner_consume_short(&edge_lens, next_v, r_target, params)
            }
        } else {
            0.0
        };
        let t_start = consume_at_start;
        let t_end = edge_len - consume_at_end;
        if t_end - t_start < params.min_frontage {
            continue;
        }

        let mut rng = rng_for_road(params.seed, road);
        let max_depth = depth_caps[i].min(params.depth + params.depth_variance.abs() + EPS_GEOM);
        let setback = params.setback.max(0.0);

        // Place split positions between t_start and t_end.
        let splits = split_positions_between(t_start, t_end, params, &mut rng);

        for k in 0..splits.len() - 1 {
            let t_a = splits[k];
            let t_b = splits[k + 1];
            if t_b - t_a < params.min_frontage - EPS_GEOM {
                continue;
            }
            let p_a = p + edge_dir * t_a;
            let p_b = p + edge_dir * t_b;
            let depth_jitter = params.depth_variance * (rng.gen::<f64>() * 2.0 - 1.0);
            let raw_depth = (params.depth + depth_jitter).max(params.min_frontage);
            let total_depth = (setback + raw_depth)
                .min(max_depth)
                .max(params.min_frontage * 0.5);
            if total_depth < params.min_frontage * 0.25 {
                continue;
            }
            let q_b = p_b + inward * total_depth;
            let q_a = p_a + inward * total_depth;

            let Ok(raw_poly) = Polygon::new(vec![p_a, p_b, q_b, q_a]) else {
                continue;
            };
            // At acute corners, bisector-clip so adjacent edges'
            // parcels can't extend past the corner bisector into each
            // other's territory. Obtuse corners are handled by their
            // corner parcel instead.
            let frontage_mid = (p_a + p_b) * 0.5;
            let mut working = raw_poly;
            if acute_corner[i] {
                if let Some(b) = corner_bisector(&verts, i) {
                    if let Some(clipped) = clip_with_bisector(&working, verts[i], b, frontage_mid) {
                        working = clipped;
                    } else {
                        continue;
                    }
                }
            }
            if acute_corner[(i + 1) % n] {
                let v_idx = (i + 1) % n;
                if let Some(b) = corner_bisector(&verts, v_idx) {
                    if let Some(clipped) =
                        clip_with_bisector(&working, verts[v_idx], b, frontage_mid)
                    {
                        working = clipped;
                    } else {
                        continue;
                    }
                }
            }
            // Clip to block boundary so parcels can never escape the
            // block (matters at acute corners and on non-rectangular
            // blocks).
            let Some(polygon) = clip_polygon_to_block(&working, &block.polygon) else {
                continue;
            };
            if polygon.area() < params.min_area {
                continue;
            }
            // After clipping, find the frontage edge — it's the one
            // that lies on the road segment within tolerance.
            let frontage_idx = match find_frontage_edge_after_clip(&polygon, p, q) {
                Some(idx) => idx,
                None => continue,
            };
            let frontage_len = {
                let v = polygon.vertices();
                (v[(frontage_idx + 1) % v.len()] - v[frontage_idx]).length()
            };
            if frontage_len < params.min_frontage {
                continue;
            }
            let edge_kinds = classify_edges(polygon.len(), frontage_idx);
            let mut parcel = Parcel {
                polygon,
                edge_kinds,
                vertex_ids: Vec::new(),
                frontage_road: road,
                frontage_edge_index: frontage_idx,
                block: face,
                building: None,
            };
            regularize_parcel(&mut parcel, params);
            out.push(parcel);
        }
    }

    let out = cleanup_block_parcel_overlaps(graph, out, params);
    Ok(out)
}

/// Polygon-difference cleanup: iterate generated parcels in placement
/// order, subtract previously-claimed territory from each. Guarantees
/// strict pairwise non-overlap (I3) regardless of upstream depth-cap
/// quirks. Order: corner parcels are generated before regular parcels,
/// so they win territory ties.
///
/// `geo`'s boolean ops occasionally panic in the internal sweep
/// algorithm on near-coincident edges. We snap inputs to a 1\,mm
/// grid (well below parcel-relevant precision) and wrap each
/// difference/union in `catch_unwind` so a failed cleanup falls
/// back to keeping the parcel as-is rather than crashing the
/// subdivision.
fn cleanup_block_parcel_overlaps(
    graph: &RoadGraph,
    parcels: Vec<Parcel>,
    params: &SubdivisionParams,
) -> Vec<Parcel> {
    use geo::Area;
    use geo::BooleanOps;
    if parcels.len() < 2 {
        return parcels;
    }
    let mut result: Vec<Parcel> = Vec::with_capacity(parcels.len());
    let mut claimed: Option<geo::MultiPolygon<f64>> = None;
    for parcel in parcels {
        let pgon = to_geo_polygon(&parcel.polygon);
        let pgon_mp = geo::MultiPolygon::new(vec![pgon]);
        let remaining = match claimed.as_ref() {
            Some(c) => {
                let safe_diff = std::panic::catch_unwind(std::panic::AssertUnwindSafe(|| {
                    pgon_mp.difference(c)
                }));
                match safe_diff {
                    Ok(r) => r,
                    Err(_) => {
                        // Boolean op crashed; we can't trust the
                        // claimed-territory tracking for this parcel,
                        // so drop it rather than risk introducing an
                        // I3 violation in subsequent parcels.
                        continue;
                    }
                }
            }
            None => pgon_mp,
        };
        let mut largest_idx: Option<usize> = None;
        let mut largest_area = 0.0;
        for (i, p) in remaining.0.iter().enumerate() {
            let a = p.unsigned_area();
            if a > largest_area {
                largest_area = a;
                largest_idx = Some(i);
            }
        }
        let Some(idx) = largest_idx else {
            continue;
        };
        if largest_area < params.min_area {
            continue;
        }
        let geo_poly = remaining.0[idx].clone();
        let Some(new_poly) = polygon_from_geo(&geo_poly) else {
            // Boolean output is too degenerate to reconstruct as a
            // valid Polygon. Drop rather than push the original
            // (which would introduce overlap with downstream parcels).
            continue;
        };
        let road = parcel.frontage_road;
        let Some((a, b)) = graph.road_nodes(road) else {
            continue;
        };
        let (Some(pos_a), Some(pos_b)) = (graph.node_position(a), graph.node_position(b)) else {
            continue;
        };
        // Cleanup pass output has been geo-snapped to a 1\,mm grid, so
        // give the matcher a couple-mm tolerance instead of the strict
        // 0.1\,mm used for unsnapped polygons.
        let frontage_idx = match find_frontage_edge_with_tol(&new_poly, pos_a, pos_b, 2.0e-3) {
            Some(i) => i,
            None => continue,
        };
        let frontage_len = {
            let v = new_poly.vertices();
            (v[(frontage_idx + 1) % v.len()] - v[frontage_idx]).length()
        };
        if frontage_len < params.min_frontage {
            continue;
        }
        let edge_kinds = classify_edges(new_poly.len(), frontage_idx);
        let updated = Parcel {
            polygon: new_poly,
            vertex_ids: Vec::new(),
            edge_kinds,
            frontage_road: parcel.frontage_road,
            frontage_edge_index: frontage_idx,
            block: parcel.block,
            building: parcel.building,
        };
        let kept_mp = geo::MultiPolygon::new(vec![geo_poly]);
        let safe_union = std::panic::catch_unwind(std::panic::AssertUnwindSafe(|| {
            match claimed.as_ref() {
                Some(c) => c.union(&kept_mp),
                None => kept_mp,
            }
        }));
        match safe_union {
            Ok(u) => {
                claimed = Some(u);
                result.push(updated);
            }
            Err(_) => {
                // Union panicked: we can't track this parcel's claim,
                // so drop it to keep downstream parcels overlap-free.
            }
        }
    }
    result
}

fn to_geo_polygon(p: &Polygon) -> geo::Polygon<f64> {
    // Snap to a 1\,mm grid before handing to `geo`. This is well
    // below any parcel-relevant precision and helps the sweep-line
    // boolean algorithm avoid near-coincident edges that trip its
    // internal invariants.
    let scale = 1000.0_f64;
    let snapped: Vec<DVec2> = p
        .vertices()
        .iter()
        .map(|v| DVec2::new((v.x * scale).round() / scale, (v.y * scale).round() / scale))
        .collect();
    let cleaned = strip_collinear_and_short(snapped);
    let mut coords: Vec<geo::Coord<f64>> = cleaned
        .iter()
        .map(|v| geo::Coord { x: v.x, y: v.y })
        .collect();
    if let Some(first) = coords.first().copied() {
        coords.push(first);
    }
    geo::Polygon::new(geo::LineString::from(coords), vec![])
}

fn polygon_from_geo(p: &geo::Polygon<f64>) -> Option<Polygon> {
    let exterior = p.exterior();
    let mut verts: Vec<DVec2> = exterior
        .points()
        .map(|pt| DVec2::new(pt.x(), pt.y()))
        .collect();
    if verts.len() >= 2 {
        let first = verts[0];
        let last = *verts.last().expect("len >= 2");
        if (first - last).length_squared() < EPS_GEOM * EPS_GEOM {
            verts.pop();
        }
    }
    let cleaned = strip_collinear_and_short(verts);
    if cleaned.len() < 3 {
        return None;
    }
    // Use strict `Polygon::new` so the polygon satisfies I1.
    Polygon::new(cleaned).ok()
}

/// Remove collinear-triple vertices and near-zero-length edges from
/// a polygon ring. Boolean-op output frequently includes both kinds
/// of artifact: tiny fp-noise edges along a straight cut, and
/// collinear interpolation points where two pieces of the same line
/// were stitched together. Strict [`Polygon::new`] (I1) rejects
/// both, so we clean them up before reconstructing the polygon.
fn strip_collinear_and_short(verts: Vec<DVec2>) -> Vec<DVec2> {
    if verts.len() < 3 {
        return verts;
    }
    // Pass 1: drop near-zero-length edges.
    let mut pass1: Vec<DVec2> = Vec::with_capacity(verts.len());
    for v in verts {
        if let Some(&last) = pass1.last() {
            if (v - last).length_squared() < EPS_GEOM * EPS_GEOM {
                continue;
            }
        }
        pass1.push(v);
    }
    if pass1.len() >= 2 {
        let first = pass1[0];
        let last = *pass1.last().expect("len >= 2");
        if (first - last).length_squared() < EPS_GEOM * EPS_GEOM {
            pass1.pop();
        }
    }
    if pass1.len() < 3 {
        return pass1;
    }
    // Pass 2: drop collinear triples (in unit-tangent cross sense).
    let mut pass2: Vec<DVec2> = Vec::with_capacity(pass1.len());
    let n = pass1.len();
    for i in 0..n {
        let prev = pass1[(i + n - 1) % n];
        let curr = pass1[i];
        let next = pass1[(i + 1) % n];
        let in_dir = (curr - prev).normalize_or_zero();
        let out_dir = (next - curr).normalize_or_zero();
        if in_dir.length_squared() < EPS_GEOM * EPS_GEOM
            || out_dir.length_squared() < EPS_GEOM * EPS_GEOM
        {
            continue;
        }
        let cross = in_dir.x * out_dir.y - in_dir.y * out_dir.x;
        // Collinear iff cross ~ 0 and they point the same way.
        if cross.abs() < 1e-3 && in_dir.dot(out_dir) > 0.0 {
            continue;
        }
        pass2.push(curr);
    }
    pass2
}

/// True iff vertex `i` of the block is a "real corner" — graph degree
/// ≥3, or degree 2 with a bend sharper than 150°. (D11.)
fn is_real_corner(graph: &RoadGraph, block: &Block, i: usize, _params: &SubdivisionParams) -> bool {
    let n = block.polygon.len();
    let half_edge_id = block.boundary_edges[i];
    let Some(he) = graph.half_edges.get(half_edge_id) else {
        return false;
    };
    let node_id = he.origin;
    let Some(node) = graph.nodes.get(node_id) else {
        return false;
    };
    if node.outgoing.len() >= 3 {
        return true;
    }
    // Degree-2: bend test.
    let v_prev = block.polygon.vertices()[(i + n - 1) % n];
    let v_curr = block.polygon.vertices()[i];
    let v_next = block.polygon.vertices()[(i + 1) % n];
    let in_dir = (v_curr - v_prev).normalize_or_zero();
    let out_dir = (v_next - v_curr).normalize_or_zero();
    if in_dir.length_squared() < EPS_GEOM * EPS_GEOM
        || out_dir.length_squared() < EPS_GEOM * EPS_GEOM
    {
        return false;
    }
    let cos_turn = in_dir.dot(out_dir).clamp(-1.0, 1.0);
    let turn_angle = cos_turn.acos();
    // turn > 30° means interior < 150° = "real corner".
    turn_angle > 30.0_f64.to_radians()
}

/// Inward bisector direction at vertex `i` of a CCW polygon.
fn corner_bisector(verts: &[DVec2], i: usize) -> Option<DVec2> {
    let n = verts.len();
    let v_prev = verts[(i + n - 1) % n];
    let v_curr = verts[i];
    let v_next = verts[(i + 1) % n];
    let in_tan = (v_curr - v_prev).normalize_or_zero();
    let out_tan = (v_next - v_curr).normalize_or_zero();
    let in_normal = DVec2::new(-in_tan.y, in_tan.x);
    let out_normal = DVec2::new(-out_tan.y, out_tan.x);
    let d = (in_normal + out_normal).normalize_or_zero();
    if d.length_squared() < EPS_GEOM * EPS_GEOM {
        return None;
    }
    Some(d)
}

/// Clip `polygon` against the half-plane defined by the bisector at
/// `corner_v` with direction `bisector_dir`, keeping the side that
/// contains `reference_inside`.
fn clip_with_bisector(
    polygon: &Polygon,
    corner_v: DVec2,
    bisector_dir: DVec2,
    reference_inside: DVec2,
) -> Option<Polygon> {
    let n_perp = DVec2::new(-bisector_dir.y, bisector_dir.x);
    let s = (reference_inside - corner_v).dot(n_perp);
    let inward_normal = if s >= 0.0 { n_perp } else { -n_perp };
    polygon.clip_half_plane(corner_v, inward_normal)
}

/// Iteratively clip `parcel` by each inward half-plane of `block`.
/// Returns `None` if the result is empty or invalid. Correct for
/// convex blocks; may over-clip mild concavities (acceptable for
/// milestone 0.2).
fn clip_polygon_to_block(parcel: &Polygon, block: &Polygon) -> Option<Polygon> {
    let mut current = parcel.clone();
    let block_verts = block.vertices();
    let n = block_verts.len();
    for i in 0..n {
        let a = block_verts[i];
        let b = block_verts[(i + 1) % n];
        let edge = b - a;
        let len = edge.length();
        if len < EPS_GEOM {
            continue;
        }
        let dir = edge / len;
        let inward = DVec2::new(-dir.y, dir.x);
        current = current.clip_half_plane(a, inward)?;
    }
    Some(current)
}

/// Locate the polygon edge that lies on the road segment
/// `(road_a, road_b)` within tolerance.
fn find_frontage_edge_after_clip(polygon: &Polygon, road_a: DVec2, road_b: DVec2) -> Option<usize> {
    find_frontage_edge_with_tol(polygon, road_a, road_b, 100.0 * EPS_GEOM)
}

fn find_frontage_edge_with_tol(
    polygon: &Polygon,
    road_a: DVec2,
    road_b: DVec2,
    tol: f64,
) -> Option<usize> {
    let road_dir = (road_b - road_a).normalize_or_zero();
    if road_dir.length_squared() < EPS_GEOM * EPS_GEOM {
        return None;
    }
    let road_normal = DVec2::new(-road_dir.y, road_dir.x);
    let v = polygon.vertices();
    let n = v.len();
    let mut best: Option<(usize, f64)> = None;
    for i in 0..n {
        let p = v[i];
        let q = v[(i + 1) % n];
        let dp = (p - road_a).dot(road_normal).abs();
        let dq = (q - road_a).dot(road_normal).abs();
        if dp < tol && dq < tol {
            let edge_dir = (q - p).normalize_or_zero();
            if edge_dir.dot(road_dir).abs() > 0.5 {
                let len = (q - p).length();
                if best.is_none_or(|(_, blen)| len > blen) {
                    best = Some((i, len));
                }
            }
        }
    }
    best.map(|(i, _)| i)
}

/// "Short" corner consumption — equal to the corner radius `R`,
/// capped at half each adjacent edge.
fn corner_consume_short(
    edge_lens: &[f64],
    i: usize,
    r_target: f64,
    params: &SubdivisionParams,
) -> f64 {
    let n = edge_lens.len();
    let prev_len = edge_lens[(i + n - 1) % n];
    let next_len = edge_lens[i];
    r_target
        .min(prev_len * 0.5)
        .min(next_len * 0.5)
        .max(params.min_frontage * 0.5)
}

/// "Long" corner consumption — equal to the parcel depth, capped at
/// half each adjacent edge.
fn corner_consume_long(edge_lens: &[f64], i: usize, params: &SubdivisionParams) -> f64 {
    let n = edge_lens.len();
    let prev_len = edge_lens[(i + n - 1) % n];
    let next_len = edge_lens[i];
    params
        .depth
        .min(prev_len * 0.5)
        .min(next_len * 0.5)
        .max(params.min_frontage * 0.5)
}

/// Per-edge depth cap (half the perpendicular distance from edge
/// midpoint to the next polygon edge along the inward normal).
fn edge_depth_caps(poly: &Polygon) -> Vec<f64> {
    let verts = poly.vertices();
    let n = verts.len();
    let mut caps = Vec::with_capacity(n);
    for i in 0..n {
        let p = verts[i];
        let q = verts[(i + 1) % n];
        let edge = q - p;
        let len = edge.length();
        if len < EPS_GEOM {
            caps.push(0.0);
            continue;
        }
        let dir = edge / len;
        let inward = DVec2::new(-dir.y, dir.x);
        let mid = (p + q) * 0.5;
        let mut min_t = f64::INFINITY;
        for j in 0..n {
            if j == i {
                continue;
            }
            let a = verts[j];
            let b = verts[(j + 1) % n];
            if let Some(t) = ray_segment_distance(mid, inward, a, b) {
                if t > EPS_GEOM && t < min_t {
                    min_t = t;
                }
            }
        }
        caps.push(min_t * 0.5);
    }
    caps
}

fn ray_segment_distance(origin: DVec2, dir: DVec2, a: DVec2, b: DVec2) -> Option<f64> {
    let v1 = origin - a;
    let v2 = b - a;
    let v3 = DVec2::new(-dir.y, dir.x);
    let denom = v2.dot(v3);
    if denom.abs() < EPS_GEOM {
        return None;
    }
    let t = v2.perp_dot(v1) / denom;
    let s = v1.dot(v3) / denom;
    if t >= 0.0 && (-EPS_GEOM..=1.0 + EPS_GEOM).contains(&s) {
        Some(t)
    } else {
        None
    }
}

/// Split positions between [`start`, `end`] along an edge. Begins
/// with `start` and ends with `end`.
fn split_positions_between(
    start: f64,
    end: f64,
    params: &SubdivisionParams,
    rng: &mut ChaCha8Rng,
) -> Vec<f64> {
    let mut splits: Vec<f64> = vec![start];
    let mut t = start;
    loop {
        let jitter = params.frontage_variance * (rng.gen::<f64>() * 2.0 - 1.0);
        let dt = (params.frontage_width + jitter).max(params.min_frontage);
        t += dt;
        if t >= end - params.min_frontage * 0.5 {
            break;
        }
        splits.push(t);
    }
    splits.push(end);
    splits
}

fn rng_for_road(seed: u64, road: RoadId) -> ChaCha8Rng {
    let road_bits = road.data().as_ffi();
    let mixed = seed
        .wrapping_mul(0x9E37_79B9_7F4A_7C15)
        .wrapping_add(road_bits)
        .wrapping_mul(0xBF58_476D_1CE4_E5B9);
    ChaCha8Rng::seed_from_u64(mixed)
}

#[allow(dead_code)]
fn _unused_node_id() -> NodeId {
    NodeId::null()
}

#[allow(dead_code)]
fn _unused_face_id() -> FaceId {
    FaceId::null()
}

#[cfg(test)]
mod tests {
    use super::*;

    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
    }

    #[test]
    fn rectangle_subdivides_into_many_parcels() {
        let g = rectangle_graph(200.0, 100.0);
        let params = SubdivisionParams::default();
        let set = subdivide_all(&g, &params).unwrap();
        assert!(set.len() > 8, "got {} parcels", set.len());
    }

    #[test]
    fn rectangle_parcels_have_one_frontage_each() {
        let g = rectangle_graph(200.0, 100.0);
        let params = SubdivisionParams::default();
        let set = subdivide_all(&g, &params).unwrap();
        for (_, p) in set.iter() {
            let frontage_count = p
                .edge_kinds()
                .iter()
                .filter(|k| matches!(k, super::super::EdgeKind::Frontage))
                .count();
            assert_eq!(frontage_count, 1, "I2 violation");
        }
    }

    #[test]
    fn rectangle_parcels_do_not_overlap() {
        let g = rectangle_graph(200.0, 100.0);
        let params = SubdivisionParams::default();
        let set = subdivide_all(&g, &params).unwrap();
        let total_area: f64 = set.iter().map(|(_, p)| p.area()).sum();
        let block_area = 200.0 * 100.0;
        assert!(
            total_area <= block_area + 1e-3,
            "parcel areas {total_area} > block area {block_area}",
        );
    }

    #[test]
    fn deterministic_with_fixed_seed() {
        let g = rectangle_graph(200.0, 100.0);
        let params = SubdivisionParams::default();
        let a = subdivide_all(&g, &params).unwrap();
        let b = subdivide_all(&g, &params).unwrap();
        assert_eq!(a.len(), b.len());
        for ((_, p1), (_, p2)) in a.iter().zip(b.iter()) {
            assert_eq!(p1.vertices().len(), p2.vertices().len());
            for (v1, v2) in p1.vertices().iter().zip(p2.vertices().iter()) {
                assert!((v1.x - v2.x).abs() < 1e-12);
                assert!((v1.y - v2.y).abs() < 1e-12);
            }
        }
    }

    #[test]
    fn invalid_params_rejected() {
        let g = rectangle_graph(50.0, 50.0);
        let params = SubdivisionParams {
            frontage_width: -1.0,
            ..SubdivisionParams::default()
        };
        let err = subdivide_all(&g, &params).unwrap_err();
        assert!(matches!(err, SubdivisionError::InvalidParams(_)));
    }

    #[test]
    fn rectangle_has_corner_parcels() {
        let g = rectangle_graph(200.0, 100.0);
        let params = SubdivisionParams::default();
        let set = subdivide_all(&g, &params).unwrap();
        // Parcel at the (0, 0) corner: its polygon should contain the
        // origin (or a vertex right at it) and have 4 vertices.
        let near_origin = set
            .iter()
            .find(|(_, p)| p.vertices().iter().any(|v| v.length_squared() < 1.0))
            .map(|(_, p)| p.vertices().len());
        assert_eq!(
            near_origin,
            Some(4),
            "expected a 4-vertex corner parcel near the origin"
        );
    }

    #[test]
    fn stats_report_nonzero_phases() {
        let g = rectangle_graph(200.0, 100.0);
        let params = SubdivisionParams::default();
        let (set, stats) = subdivide_all_with_stats(&g, &params).unwrap();
        assert_eq!(stats.parcel_count, set.len());
        assert!(stats.block_count >= 1);
        assert!(stats.total > Duration::ZERO);
        assert!(stats.parcels_per_second() > 0.0);
    }
}
