rustial_engine/

tessellator.rs

// ---------------------------------------------------------------------------
//! # Vector tessellation -- polygon fill and line stroke
//!
//! This module converts geographic vector geometries into GPU-ready
//! triangle meshes.  It provides the following public functions:
//!
//! - [`triangulate_polygon`] -- fill a simple polygon with triangles
//!   using **ear-clipping** (handles concave polygons correctly).
//! - [`triangulate_polygon_with_holes`] -- fill a polygon with interior
//!   rings (holes) using ear-clipping.
//! - [`stroke_line`] -- expand a polyline into a thick ribbon of
//!   triangles with configurable half-width.
//!
//! Both operate in the **degree plane** (lon as X, lat as Y).  The
//! caller ([`VectorLayer::tessellate`](crate::layers::VectorLayer::tessellate))
//! projects the resulting vertices into Web Mercator world space before
//! GPU upload.
//!
//! ## Coordinate space
//!
//! Positions are in raw (lon, lat) degrees.  At the equator 1 degree
//! ~ 111 km in both axes, so the geometry is approximately isotropic.
//! For high-latitude features the caller should consider projecting
//! first, but for typical map-display use cases the distortion is
//! masked by the Mercator projection.
// ---------------------------------------------------------------------------

use crate::layers::{LineCap, LineJoin};
use rustial_math::GeoCoord;

// ---------------------------------------------------------------------------
// Internal helpers
// ---------------------------------------------------------------------------

/// A lightweight 2D vector used for tangent / normal calculations.
///
/// Intentionally kept private -- callers work with `GeoCoord` and the
/// `[f64; 2]` output of [`stroke_line`].
#[derive(Debug, Clone, Copy)]
struct Vec2 {
    x: f64,
    y: f64,
}

/// Normalize a 2D vector to unit length.
///
/// Returns the zero vector when the input length is below `1e-15`
/// (sub-nanometre in degree space), avoiding division by zero for
/// coincident or nearly-coincident points.
#[inline]
fn normalize(v: Vec2) -> Vec2 {
    let len = (v.x * v.x + v.y * v.y).sqrt();
    if len < 1e-15 {
        Vec2 { x: 0.0, y: 0.0 }
    } else {
        Vec2 {
            x: v.x / len,
            y: v.y / len,
        }
    }
}

// ---------------------------------------------------------------------------
// Public API -- polygon triangulation (ear-clipping)
// ---------------------------------------------------------------------------

/// Strip the closing duplicate vertex from a ring if present.
///
/// Returns the effective vertex count after removing a closing duplicate
/// that matches the first vertex within 1e-12 degrees.  The GeoJSON
/// convention repeats the first vertex at the end; this function detects
/// that and excludes it from triangulation.
fn effective_ring_len(coords: &[GeoCoord]) -> usize {
    let n = coords.len();
    if n > 3
        && (coords[0].lat - coords[n - 1].lat).abs() < 1e-12
        && (coords[0].lon - coords[n - 1].lon).abs() < 1e-12
    {
        n - 1
    } else {
        n
    }
}

/// Signed area of a ring (lon as X, lat as Y) using the shoelace formula.
/// Positive for counter-clockwise, negative for clockwise.
fn _signed_ring_area(coords: &[GeoCoord], len: usize) -> f64 {
    let mut area = 0.0;
    let mut j = len - 1;
    for i in 0..len {
        area += (coords[j].lon - coords[i].lon) * (coords[j].lat + coords[i].lat);
        j = i;
    }
    area * 0.5
}

/// Test whether point `p` lies strictly inside triangle `(a, b, c)`.
///
/// Uses the cross-product sign test.  Points exactly on an edge return
/// `false` to avoid degenerate ear removal.
fn point_in_triangle(
    px: f64, py: f64,
    ax: f64, ay: f64,
    bx: f64, by: f64,
    cx: f64, cy: f64,
) -> bool {
    let d1 = (px - bx) * (ay - by) - (ax - bx) * (py - by);
    let d2 = (px - cx) * (by - cy) - (bx - cx) * (py - cy);
    let d3 = (px - ax) * (cy - ay) - (cx - ax) * (py - ay);
    let has_neg = (d1 < 0.0) || (d2 < 0.0) || (d3 < 0.0);
    let has_pos = (d1 > 0.0) || (d2 > 0.0) || (d3 > 0.0);
    !(has_neg && has_pos)
}

/// Cross product of vectors (b - a) and (c - a), returning just the Z component.
fn cross_z(ax: f64, ay: f64, bx: f64, by: f64, cx: f64, cy: f64) -> f64 {
    (bx - ax) * (cy - ay) - (by - ay) * (cx - ax)
}

/// Ear-clipping triangulation of a flat coordinate array with hole bridging.
///
/// `flat` contains `[lon, lat, lon, lat, ...]` for the exterior ring
/// followed by each hole ring.  `hole_starts` contains the index into
/// `flat` (counting pairs, i.e. vertex indices) where each hole begins.
///
/// Returns triangle indices into the flat vertex array (vertex index,
/// not byte index).
fn earcut_flat(flat: &[f64], hole_starts: &[usize]) -> Vec<u32> {
    let total_verts = flat.len() / 2;
    if total_verts < 3 {
        return Vec::new();
    }

    // Build the initial doubly-linked list of vertex indices.
    let mut prev = vec![0usize; total_verts];
    let mut next = vec![0usize; total_verts];

    // Link the exterior ring (vertices 0..hole_starts[0] or 0..total_verts).
    let outer_end = if hole_starts.is_empty() {
        total_verts
    } else {
        hole_starts[0]
    };

    if outer_end < 3 {
        return Vec::new();
    }

    // Ensure CCW winding for the outer ring.
    let outer_ccw = ring_area_flat(flat, 0, outer_end) > 0.0;
    link_ring(&mut prev, &mut next, 0, outer_end, outer_ccw);

    // Process holes: ensure CW winding, then bridge each hole into the
    // outer polygon.
    let mut outer_start = 0usize;
    for (hi, &h_start) in hole_starts.iter().enumerate() {
        let h_end = if hi + 1 < hole_starts.len() {
            hole_starts[hi + 1]
        } else {
            total_verts
        };
        if h_end - h_start < 3 {
            continue;
        }
        let hole_ccw = ring_area_flat(flat, h_start, h_end) > 0.0;
        // Holes should be CW (negative area in our convention).
        link_ring(&mut prev, &mut next, h_start, h_end, !hole_ccw);

        // Bridge: find the rightmost vertex of the hole, then find a
        // visible vertex on the outer boundary to bridge to.
        let bridge_hole = rightmost_vertex(flat, h_start, h_end);
        outer_start = bridge_hole_to_outer(
            flat, &mut prev, &mut next, outer_start, outer_end, bridge_hole,
        );
    }

    // Run ear-clipping on the linked list.
    let mut indices = Vec::with_capacity((total_verts - 2) * 3);
    ear_clip(flat, &prev, &mut next.clone(), outer_start, &mut indices, total_verts);
    indices
}

/// Signed area of a ring in a flat `[lon, lat, ...]` array.
fn ring_area_flat(flat: &[f64], start: usize, end: usize) -> f64 {
    let mut area = 0.0;
    let mut j = end - 1;
    for i in start..end {
        let jx = flat[j * 2];
        let jy = flat[j * 2 + 1];
        let ix = flat[i * 2];
        let iy = flat[i * 2 + 1];
        area += (jx - ix) * (jy + iy);
        j = i;
    }
    area * 0.5
}

/// Link a ring into a circular doubly-linked list.
/// If `forward` is true, link in index order; otherwise reverse.
fn link_ring(prev: &mut [usize], next: &mut [usize], start: usize, end: usize, forward: bool) {
    if forward {
        for i in start..end {
            let p = if i == start { end - 1 } else { i - 1 };
            let n = if i == end - 1 { start } else { i + 1 };
            prev[i] = p;
            next[i] = n;
        }
    } else {
        for i in start..end {
            let p = if i == end - 1 { start } else { i + 1 };
            let n = if i == start { end - 1 } else { i - 1 };
            prev[i] = p;
            next[i] = n;
        }
    }
}

/// Find the vertex with the largest X (lon) in a ring.
fn rightmost_vertex(flat: &[f64], start: usize, end: usize) -> usize {
    let mut best = start;
    for i in (start + 1)..end {
        if flat[i * 2] > flat[best * 2]
            || (flat[i * 2] == flat[best * 2] && flat[i * 2 + 1] < flat[best * 2 + 1])
        {
            best = i;
        }
    }
    best
}

/// Bridge a hole vertex into the outer polygon by inserting mutual links.
fn bridge_hole_to_outer(
    flat: &[f64],
    prev: &mut [usize],
    next: &mut [usize],
    outer_start: usize,
    _outer_end: usize,
    hole_vertex: usize,
) -> usize {
    // Find the outer vertex closest to `hole_vertex` that is visible.
    let hx = flat[hole_vertex * 2];
    let hy = flat[hole_vertex * 2 + 1];
    let mut best = outer_start;
    let mut best_dist = f64::INFINITY;
    let mut i = outer_start;
    loop {
        let ix = flat[i * 2];
        let iy = flat[i * 2 + 1];
        let d = (ix - hx) * (ix - hx) + (iy - hy) * (iy - hy);
        if d < best_dist {
            best_dist = d;
            best = i;
        }
        i = next[i];
        if i == outer_start {
            break;
        }
    }

    // Splice the hole into the outer ring at `best`.
    // After bridging: ... -> best -> hole_vertex -> ... hole ring ... -> hole_vertex -> best -> ...
    // We create two copies conceptually but reuse the same indices by
    // re-linking:  best.next = hole_vertex, and at the end of the hole
    // ring we link back to best.
    let best_next = next[best];
    let hole_prev = prev[hole_vertex];

    next[best] = hole_vertex;
    prev[hole_vertex] = best;
    next[hole_prev] = best_next;
    prev[best_next] = hole_prev;

    best
}

/// Core ear-clipping loop.
fn ear_clip(
    flat: &[f64],
    orig_prev: &[usize],
    next: &mut [usize],
    start: usize,
    indices: &mut Vec<u32>,
    _total_verts: usize,
) {
    // We need a mutable prev as well.
    let mut prev = orig_prev.to_vec();
    let mut remaining = {
        // Count vertices in the linked list starting from `start`.
        let mut count = 1usize;
        let mut i = next[start];
        while i != start {
            count += 1;
            i = next[i];
        }
        count
    };

    let mut ear = start;
    let mut _stop = ear;
    let mut pass = 0;

    while remaining > 2 {
        let a = prev[ear];
        let b = ear;
        let c = next[ear];

        if is_ear(flat, &prev, next, a, b, c) {
            indices.push(a as u32);
            indices.push(b as u32);
            indices.push(c as u32);

            // Remove vertex b.
            next[a] = c;
            prev[c] = a;

            remaining -= 1;
            _stop = c;
            ear = c;
            pass = 0;
        } else {
            ear = next[ear];
            pass += 1;
            if pass >= remaining {
                // No ear found in a full pass -- degenerate polygon.
                // Fall back to fan triangulation of remaining vertices.
                fan_remaining(next, start, ear, remaining, indices);
                break;
            }
        }
    }
}

/// Check if vertex `b` is an ear (the triangle a-b-c contains no other vertices).
fn is_ear(
    flat: &[f64],
    prev: &[usize],
    next: &[usize],
    a: usize,
    b: usize,
    c: usize,
) -> bool {
    let ax = flat[a * 2];
    let ay = flat[a * 2 + 1];
    let bx = flat[b * 2];
    let by = flat[b * 2 + 1];
    let cx = flat[c * 2];
    let cy = flat[c * 2 + 1];

    // The ear triangle must be convex (positive cross product for CCW).
    if cross_z(ax, ay, bx, by, cx, cy) <= 0.0 {
        return false;
    }

    // Check that no other vertex lies inside the triangle.
    let mut p = next[c];
    while p != a {
        let px = flat[p * 2];
        let py = flat[p * 2 + 1];
        if point_in_triangle(px, py, ax, ay, bx, by, cx, cy)
            && cross_z(
                flat[prev[p] * 2], flat[prev[p] * 2 + 1],
                px, py,
                flat[next[p] * 2], flat[next[p] * 2 + 1],
            ) >= 0.0
        {
            return false;
        }
        p = next[p];
    }
    true
}

/// Fan-triangulate the remaining linked-list vertices as a last resort.
fn fan_remaining(
    next: &[usize],
    _start: usize,
    first: usize,
    remaining: usize,
    indices: &mut Vec<u32>,
) {
    if remaining < 3 {
        return;
    }
    let anchor = first;
    let mut b = next[anchor];
    for _ in 0..remaining - 2 {
        let c = next[b];
        indices.push(anchor as u32);
        indices.push(b as u32);
        indices.push(c as u32);
        b = c;
    }
}

/// Triangulate a simple polygon (no holes) into a list of triangle indices.
///
/// Uses **ear-clipping** which handles both convex and concave polygons
/// correctly.  For polygons with holes, use
/// [`triangulate_polygon_with_holes`] instead.
///
/// # Returns
///
/// Indices into the input `coords` slice, grouped in triples.
///
/// # Closing vertex
///
/// If the last coordinate duplicates the first (within 1e-12 degrees),
/// it is treated as a ring-closing sentinel and excluded from
/// triangulation.  This matches the GeoJSON convention where polygon
/// rings repeat the first vertex.
///
/// # Edge cases
///
/// | Input | Result |
/// |-------|--------|
/// | Fewer than 3 coordinates | empty `Vec` |
/// | Exactly 3 unique coordinates | 1 triangle (3 indices) |
/// | Closing duplicate reducing count below 3 | empty `Vec` |
///
/// # Example
///
/// ```
/// use rustial_engine::triangulate_polygon;
/// use rustial_engine::GeoCoord;
///
/// let square = vec![
///     GeoCoord::from_lat_lon(0.0, 0.0),
///     GeoCoord::from_lat_lon(0.0, 1.0),
///     GeoCoord::from_lat_lon(1.0, 1.0),
///     GeoCoord::from_lat_lon(1.0, 0.0),
/// ];
/// let indices = triangulate_polygon(&square);
/// assert_eq!(indices.len(), 6); // 2 triangles
/// ```
pub fn triangulate_polygon(coords: &[GeoCoord]) -> Vec<u32> {
    let n = coords.len();
    if n < 3 {
        return Vec::new();
    }

    let effective_len = effective_ring_len(coords);
    if effective_len < 3 {
        return Vec::new();
    }

    // Build flat coordinate array [lon, lat, lon, lat, ...].
    let mut flat = Vec::with_capacity(effective_len * 2);
    for c in &coords[..effective_len] {
        flat.push(c.lon);
        flat.push(c.lat);
    }

    earcut_flat(&flat, &[])
}

/// Triangulate a polygon with interior rings (holes).
///
/// The exterior ring and each hole are provided as separate slices.
/// Closing duplicate vertices are stripped automatically.  Winding
/// order is normalised internally (exterior CCW, holes CW).
///
/// Returns indices into a combined vertex array where the exterior
/// ring vertices come first (`0..exterior.effective_len`), followed
/// by hole vertices in order.
///
/// # Example
///
/// ```
/// use rustial_engine::triangulate_polygon_with_holes;
/// use rustial_engine::GeoCoord;
///
/// let exterior = vec![
///     GeoCoord::from_lat_lon(0.0, 0.0),
///     GeoCoord::from_lat_lon(0.0, 4.0),
///     GeoCoord::from_lat_lon(4.0, 4.0),
///     GeoCoord::from_lat_lon(4.0, 0.0),
/// ];
/// let hole = vec![
///     GeoCoord::from_lat_lon(1.0, 1.0),
///     GeoCoord::from_lat_lon(1.0, 3.0),
///     GeoCoord::from_lat_lon(3.0, 3.0),
///     GeoCoord::from_lat_lon(3.0, 1.0),
/// ];
/// let indices = triangulate_polygon_with_holes(&exterior, &[&hole]);
/// assert!(!indices.is_empty());
/// ```
pub fn triangulate_polygon_with_holes(
    exterior: &[GeoCoord],
    holes: &[&[GeoCoord]],
) -> Vec<u32> {
    let ext_len = effective_ring_len(exterior);
    if ext_len < 3 {
        return Vec::new();
    }

    let hole_lens: Vec<usize> = holes.iter().map(|h| effective_ring_len(h)).collect();
    let total_verts: usize = ext_len + hole_lens.iter().sum::<usize>();
    let mut flat = Vec::with_capacity(total_verts * 2);

    // Exterior ring.
    for c in &exterior[..ext_len] {
        flat.push(c.lon);
        flat.push(c.lat);
    }

    // Holes.
    let mut hole_starts = Vec::with_capacity(holes.len());
    let mut offset = ext_len;
    for (i, hole) in holes.iter().enumerate() {
        let hl = hole_lens[i];
        if hl < 3 {
            continue;
        }
        hole_starts.push(offset);
        for c in &hole[..hl] {
            flat.push(c.lon);
            flat.push(c.lat);
        }
        offset += hl;
    }

    earcut_flat(&flat, &hole_starts)
}

// ---------------------------------------------------------------------------
// Public API -- line stroke
// ---------------------------------------------------------------------------

/// Expand a polyline into a thick triangle-strip ribbon.
///
/// Each input vertex is extruded along its local normal (perpendicular
/// to the tangent) by `+/- half_width`, producing two output vertices.
/// Adjacent quads are then connected with two triangles each.
///
/// # Tangent computation
///
/// | Vertex position | Tangent source |
/// |-----------------|----------------|
/// | First | Forward difference to next vertex |
/// | Interior | Central difference (average of prev->curr and curr->next) |
/// | Last | Backward difference from previous vertex |
///
/// At sharp bends the averaged tangent may cause the ribbon to narrow
/// or self-intersect.  A miter-limit strategy is planned for a future
/// release.
///
/// # Returns
///
/// `(positions, indices)` where:
///
/// - `positions` -- `[lon, lat]` pairs in degree space, two per input
///   vertex (left and right of the centreline).
/// - `indices` -- triangle indices into `positions` (6 per segment).
///
/// # Edge cases
///
/// | Input | Result |
/// |-------|--------|
/// | Fewer than 2 coordinates | `(empty, empty)` |
/// | `half_width <= 0.0` | Degenerate zero-area ribbon (positions collapse onto the centreline) |
/// | Coincident consecutive vertices | Zero-length tangent yields zero normal -- the two extruded vertices coincide, producing degenerate triangles that are invisible on screen |
///
/// # Example
///
/// ```
/// use rustial_engine::stroke_line;
/// use rustial_engine::GeoCoord;
///
/// let line = vec![
///     GeoCoord::from_lat_lon(0.0, 0.0),
///     GeoCoord::from_lat_lon(0.0, 1.0),
/// ];
/// let (positions, indices) = stroke_line(&line, 0.01);
/// assert_eq!(positions.len(), 4); // 2 vertices * 2 sides
/// assert_eq!(indices.len(), 6);   // 1 segment * 6 indices
/// ```
pub fn stroke_line(coords: &[GeoCoord], half_width: f64) -> (Vec<[f64; 2]>, Vec<u32>) {
    if coords.len() < 2 {
        return (Vec::new(), Vec::new());
    }

    let vertex_count = coords.len() * 2;
    let segment_count = coords.len() - 1;
    let mut positions = Vec::with_capacity(vertex_count);
    let mut indices = Vec::with_capacity(segment_count * 6);

    for i in 0..coords.len() {
        let curr = Vec2 {
            x: coords[i].lon,
            y: coords[i].lat,
        };

        // Compute the tangent direction at this vertex.
        let tangent = if i == 0 {
            // First vertex: forward difference.
            let next = Vec2 {
                x: coords[1].lon,
                y: coords[1].lat,
            };
            normalize(Vec2 {
                x: next.x - curr.x,
                y: next.y - curr.y,
            })
        } else if i == coords.len() - 1 {
            // Last vertex: backward difference.
            let prev = Vec2 {
                x: coords[i - 1].lon,
                y: coords[i - 1].lat,
            };
            normalize(Vec2 {
                x: curr.x - prev.x,
                y: curr.y - prev.y,
            })
        } else {
            // Interior vertex: central difference (averaged tangent).
            let prev = Vec2 {
                x: coords[i - 1].lon,
                y: coords[i - 1].lat,
            };
            let next = Vec2 {
                x: coords[i + 1].lon,
                y: coords[i + 1].lat,
            };
            normalize(Vec2 {
                x: next.x - prev.x,
                y: next.y - prev.y,
            })
        };

        // The extrusion normal is the 90-degree rotation of the tangent.
        let normal = Vec2 {
            x: -tangent.y,
            y: tangent.x,
        };

        // Left and right extruded positions.
        positions.push([
            curr.x + normal.x * half_width,
            curr.y + normal.y * half_width,
        ]);
        positions.push([
            curr.x - normal.x * half_width,
            curr.y - normal.y * half_width,
        ]);
    }

    // Connect each pair of extruded vertices into a quad (2 triangles).
    //
    //   left[i]  ---- left[i+1]        Indices per quad:
    //      |    \         |              (base, base+1, base+2)
    //      |     \        |              (base+1, base+3, base+2)
    //   right[i] --- right[i+1]
    //
    for i in 0..segment_count as u32 {
        let base = i * 2;
        indices.push(base);
        indices.push(base + 1);
        indices.push(base + 2);
        indices.push(base + 1);
        indices.push(base + 3);
        indices.push(base + 2);
    }

    (positions, indices)
}

// ---------------------------------------------------------------------------
// stroke_line_styled — full line tessellation with cap, join, and per-vertex
//                      line normals and cumulative distances
// ---------------------------------------------------------------------------

/// Output of [`stroke_line_styled`].
#[derive(Debug, Clone)]
pub struct StrokeLineResult {
    /// Vertex positions `[lon, lat]` in degree space.
    pub positions: Vec<[f64; 2]>,
    /// Triangle indices into `positions`.
    pub indices: Vec<u32>,
    /// Per-vertex extrusion normal `[nx, ny]` (unit-length, perpendicular to
    /// centreline).  Same length as `positions`.
    pub normals: Vec<[f64; 2]>,
    /// Per-vertex cumulative distance along the polyline centreline (degrees).
    /// Same length as `positions`.
    pub distances: Vec<f64>,
    /// Per-vertex cap/join flag.  `1.0` for vertices belonging to round
    /// cap or round join fan geometry (center and perimeter), `0.0` for
    /// ribbon body, bevel, miter, and square cap vertices.
    ///
    /// The GPU shader uses this flag to switch from linear edge AA
    /// (ribbon body) to SDF circle-based AA (round caps/joins).
    pub cap_join: Vec<f32>,
}

/// Number of triangular fan segments used for round caps and round joins.
const ROUND_SEGMENTS: u32 = 8;

/// Circumscription scale factor: fan perimeter vertices are placed at
/// `half_width * CIRCUMSCRIBE` so the polygon circumscribes the ideal
/// circle.  The fragment shader then clips at SDF distance = 1.0,
/// producing a pixel-perfect round edge regardless of segment count.
const CIRCUMSCRIBE: f64 = {
    // 1 / cos(π / (2 * ROUND_SEGMENTS))
    // For 8 segments: 1 / cos(π/16) ≈ 1.01959
    // Precomputed as a const because f64::cos isn't const-fn.
    1.01959
};

/// Expand a polyline into a thick triangle ribbon with proper line caps,
/// line joins, per-vertex extrusion normals, and cumulative distances.
///
/// This is the richer counterpart of [`stroke_line`].  The additional
/// outputs enable the GPU line shader to evaluate dash patterns (via
/// distances) and perform antialiasing (via normals), while the
/// tessellation itself produces correct cap and join geometry.
///
/// # Parameters
///
/// - `coords` — polyline vertices in `(lon, lat)` degree space.
/// - `half_width` — extrusion half-width in the same degree-space units.
/// - `cap` — line cap style applied at both endpoints.
/// - `join` — line join style applied at interior vertices.
/// - `miter_limit` — ratio threshold for automatic miter → bevel fallback.
pub fn stroke_line_styled(
    coords: &[GeoCoord],
    half_width: f64,
    cap: LineCap,
    join: LineJoin,
    miter_limit: f32,
) -> StrokeLineResult {
    let empty = StrokeLineResult {
        positions: Vec::new(),
        indices: Vec::new(),
        normals: Vec::new(),
        distances: Vec::new(),
        cap_join: Vec::new(),
    };
    if coords.len() < 2 {
        return empty;
    }

    // Pre-compute per-segment tangent, normal, and length.
    let seg_count = coords.len() - 1;
    let mut seg_tangents = Vec::with_capacity(seg_count);
    let mut seg_normals = Vec::with_capacity(seg_count);
    let mut seg_lengths = Vec::with_capacity(seg_count);

    for i in 0..seg_count {
        let dx = coords[i + 1].lon - coords[i].lon;
        let dy = coords[i + 1].lat - coords[i].lat;
        let len = (dx * dx + dy * dy).sqrt();
        let t = if len < 1e-15 {
            Vec2 { x: 1.0, y: 0.0 }
        } else {
            Vec2 { x: dx / len, y: dy / len }
        };
        let n = Vec2 { x: -t.y, y: t.x };
        seg_tangents.push(t);
        seg_normals.push(n);
        seg_lengths.push(len);
    }

    // Cumulative distance at each vertex along the centreline.
    let mut cum_dist = Vec::with_capacity(coords.len());
    cum_dist.push(0.0);
    for i in 0..seg_count {
        cum_dist.push(cum_dist[i] + seg_lengths[i]);
    }

    // --- Build output arrays ---
    // Conservative capacity: 2 verts per original vertex + cap + join extras.
    let cap_extra = match cap { LineCap::Round => ROUND_SEGMENTS as usize * 2, _ => 2 };
    let join_extra = match join { LineJoin::Round => ROUND_SEGMENTS as usize, _ => 2 };
    let est_verts = coords.len() * 2 + cap_extra * 2 + join_extra * seg_count;
    let est_indices = seg_count * 6 + cap_extra * 6 + join_extra * 3 * seg_count;
    let mut positions = Vec::with_capacity(est_verts);
    let mut normals = Vec::with_capacity(est_verts);
    let mut distances = Vec::with_capacity(est_verts);
    let mut cap_join_flags = Vec::with_capacity(est_verts);
    let mut indices = Vec::with_capacity(est_indices);

    // Helper: push a vertex and return its index.
    #[inline]
    fn push_vert(
        positions: &mut Vec<[f64; 2]>,
        normals: &mut Vec<[f64; 2]>,
        distances: &mut Vec<f64>,
        cap_join_flags: &mut Vec<f32>,
        pos: [f64; 2],
        nrm: [f64; 2],
        dist: f64,
        cj: f32,
    ) -> u32 {
        let idx = positions.len() as u32;
        positions.push(pos);
        normals.push(nrm);
        distances.push(dist);
        cap_join_flags.push(cj);
        idx
    }

    // --- Start cap ---
    let first_n = seg_normals[0];
    let first_t = seg_tangents[0];
    let cx = coords[0].lon;
    let cy = coords[0].lat;

    match cap {
        LineCap::Butt => {
            // Two vertices at the start, perpendicular to the first segment.
            let l = push_vert(&mut positions, &mut normals, &mut distances, &mut cap_join_flags,
                [cx + first_n.x * half_width, cy + first_n.y * half_width],
                [first_n.x, first_n.y], 0.0, 0.0);
            let r = push_vert(&mut positions, &mut normals, &mut distances, &mut cap_join_flags,
                [cx - first_n.x * half_width, cy - first_n.y * half_width],
                [-first_n.x, -first_n.y], 0.0, 0.0);
            let _ = (l, r); // stored in the buffer for the first segment
        }
        LineCap::Square => {
            // Extend half_width backwards along the tangent.
            let bx = cx - first_t.x * half_width;
            let by = cy - first_t.y * half_width;
            push_vert(&mut positions, &mut normals, &mut distances, &mut cap_join_flags,
                [bx + first_n.x * half_width, by + first_n.y * half_width],
                [first_n.x, first_n.y], 0.0, 0.0);
            push_vert(&mut positions, &mut normals, &mut distances, &mut cap_join_flags,
                [bx - first_n.x * half_width, by - first_n.y * half_width],
                [-first_n.x, -first_n.y], 0.0, 0.0);
        }
        LineCap::Round => {
            // Semicircular fan at the start, sweeping from left-normal
            // backwards through -tangent to right-normal.
            let center_idx = push_vert(&mut positions, &mut normals, &mut distances, &mut cap_join_flags,
                [cx, cy], [0.0, 0.0], 0.0, 1.0);
            // Start angle: direction of left normal.
            let start_angle = first_n.y.atan2(first_n.x);
            // We sweep π radians (semicircle) on the back side.
            // Fan perimeter vertices are placed at the circumscribed radius
            // so the polygon extends slightly beyond the ideal circle.
            // The shader's SDF clip at normal_length = 1.0 produces a
            // pixel-perfect round edge.
            let hw_circ = half_width * CIRCUMSCRIBE;
            let mut fan_verts = Vec::with_capacity(ROUND_SEGMENTS as usize + 1);
            for k in 0..=ROUND_SEGMENTS {
                let a = start_angle + std::f64::consts::PI * k as f64 / ROUND_SEGMENTS as f64;
                let nx = a.cos();
                let ny = a.sin();
                let v = push_vert(&mut positions, &mut normals, &mut distances, &mut cap_join_flags,
                    [cx + nx * hw_circ, cy + ny * hw_circ],
                    [nx * CIRCUMSCRIBE, ny * CIRCUMSCRIBE], 0.0, 1.0);
                fan_verts.push(v);
            }
            for k in 0..ROUND_SEGMENTS {
                indices.push(center_idx);
                indices.push(fan_verts[k as usize + 1]);
                indices.push(fan_verts[k as usize]);
            }
            // The first and last fan vertices become the left/right of the
            // ribbon at vertex 0.  We need two "current" vertices for the
            // quad strip — reuse the first and last fan vertices.
            // Push the ribbon-start pair referencing the same positions.
            push_vert(&mut positions, &mut normals, &mut distances, &mut cap_join_flags,
                [cx + first_n.x * half_width, cy + first_n.y * half_width],
                [first_n.x, first_n.y], 0.0, 0.0);
            push_vert(&mut positions, &mut normals, &mut distances, &mut cap_join_flags,
                [cx - first_n.x * half_width, cy - first_n.y * half_width],
                [-first_n.x, -first_n.y], 0.0, 0.0);
        }
    }

    // The last two vertices in the buffer are always the current left/right
    // pair that the next quad strip segment connects from.
    // After the start cap they are at positions.len()-2 and positions.len()-1.

    // --- Interior segments with joins ---
    for i in 0..seg_count {
        let n = seg_normals[i];
        let dist = cum_dist[i + 1];
        let vx = coords[i + 1].lon;
        let vy = coords[i + 1].lat;

        // The previous left/right pair.
        let prev_left = positions.len() as u32 - 2;
        let prev_right = positions.len() as u32 - 1;

        if i < seg_count - 1 {
            // Interior vertex — need a join.
            let n_next = seg_normals[i + 1];

            // Determine the miter direction (bisector normal).
            let bx = n.x + n_next.x;
            let by = n.y + n_next.y;
            let blen = (bx * bx + by * by).sqrt();

            // Cross product of segment tangents to determine turn direction.
            let cross = seg_tangents[i].x * seg_tangents[i + 1].y
                      - seg_tangents[i].y * seg_tangents[i + 1].x;

            if blen < 1e-12 {
                // Near-reversal: segments are anti-parallel.  Use bevel.
                let l = push_vert(&mut positions, &mut normals, &mut distances, &mut cap_join_flags,
                    [vx + n.x * half_width, vy + n.y * half_width],
                    [n.x, n.y], dist, 0.0);
                let r = push_vert(&mut positions, &mut normals, &mut distances, &mut cap_join_flags,
                    [vx - n.x * half_width, vy - n.y * half_width],
                    [-n.x, -n.y], dist, 0.0);
                // Quad from previous pair to current.
                indices.extend_from_slice(&[prev_left, prev_right, l, prev_right, r, l]);
                // Bevel triangle to next-segment normal.
                let l2 = push_vert(&mut positions, &mut normals, &mut distances, &mut cap_join_flags,
                    [vx + n_next.x * half_width, vy + n_next.y * half_width],
                    [n_next.x, n_next.y], dist, 0.0);
                let r2 = push_vert(&mut positions, &mut normals, &mut distances, &mut cap_join_flags,
                    [vx - n_next.x * half_width, vy - n_next.y * half_width],
                    [-n_next.x, -n_next.y], dist, 0.0);
                indices.extend_from_slice(&[l, r, l2, r, r2, l2]);
            } else {
                let miter_nx = bx / blen;
                let miter_ny = by / blen;
                // Miter length ratio = 1 / sin(half-angle).
                let dot = n.x * miter_nx + n.y * miter_ny;
                let miter_len = if dot.abs() > 1e-12 { 1.0 / dot } else { miter_limit as f64 + 1.0 };

                let use_miter = matches!(join, LineJoin::Miter) && miter_len <= miter_limit as f64;

                if use_miter {
                    // Miter join: single vertex pair along the bisector.
                    let hw_m = half_width * miter_len;
                    let l = push_vert(&mut positions, &mut normals, &mut distances, &mut cap_join_flags,
                        [vx + miter_nx * hw_m, vy + miter_ny * hw_m],
                        [miter_nx, miter_ny], dist, 0.0);
                    let r = push_vert(&mut positions, &mut normals, &mut distances, &mut cap_join_flags,
                        [vx - miter_nx * hw_m, vy - miter_ny * hw_m],
                        [-miter_nx, -miter_ny], dist, 0.0);
                    indices.extend_from_slice(&[prev_left, prev_right, l, prev_right, r, l]);
                } else if matches!(join, LineJoin::Round) {
                    // Close the current segment with the incoming normal.
                    let l_in = push_vert(&mut positions, &mut normals, &mut distances, &mut cap_join_flags,
                        [vx + n.x * half_width, vy + n.y * half_width],
                        [n.x, n.y], dist, 0.0);
                    let r_in = push_vert(&mut positions, &mut normals, &mut distances, &mut cap_join_flags,
                        [vx - n.x * half_width, vy - n.y * half_width],
                        [-n.x, -n.y], dist, 0.0);
                    indices.extend_from_slice(&[prev_left, prev_right, l_in, prev_right, r_in, l_in]);

                    // Round join fan on the outer side.
                    // Fan perimeter vertices use circumscribed radius for SDF AA.
                    let hw_circ = half_width * CIRCUMSCRIBE;
                    let center_idx = push_vert(&mut positions, &mut normals, &mut distances, &mut cap_join_flags,
                        [vx, vy], [0.0, 0.0], dist, 1.0);

                    if cross > 0.0 {
                        // Left turn: outer side is the left (positive normal).
                        let a0 = n.y.atan2(n.x);
                        let a1 = n_next.y.atan2(n_next.x);
                        let mut da = a1 - a0;
                        if da < 0.0 { da += 2.0 * std::f64::consts::PI; }
                        let steps = ROUND_SEGMENTS;
                        let mut prev_v = l_in;
                        for k in 1..=steps {
                            let a = a0 + da * k as f64 / steps as f64;
                            let nx = a.cos();
                            let ny = a.sin();
                            let v = push_vert(&mut positions, &mut normals, &mut distances, &mut cap_join_flags,
                                [vx + nx * hw_circ, vy + ny * hw_circ],
                                [nx * CIRCUMSCRIBE, ny * CIRCUMSCRIBE], dist, 1.0);
                            indices.extend_from_slice(&[center_idx, prev_v, v]);
                            prev_v = v;
                        }
                        // The new left/right pair for the next segment.
                        push_vert(&mut positions, &mut normals, &mut distances, &mut cap_join_flags,
                            [vx + n_next.x * half_width, vy + n_next.y * half_width],
                            [n_next.x, n_next.y], dist, 0.0);
                        push_vert(&mut positions, &mut normals, &mut distances, &mut cap_join_flags,
                            [vx - n_next.x * half_width, vy - n_next.y * half_width],
                            [-n_next.x, -n_next.y], dist, 0.0);
                    } else {
                        // Right turn: outer side is the right (negative normal).
                        let a0 = (-n.y).atan2(-n.x);
                        let a1 = (-n_next.y).atan2(-n_next.x);
                        let mut da = a1 - a0;
                        if da > 0.0 { da -= 2.0 * std::f64::consts::PI; }
                        let steps = ROUND_SEGMENTS;
                        let mut prev_v = r_in;
                        for k in 1..=steps {
                            let a = a0 + da * k as f64 / steps as f64;
                            let nx = a.cos();
                            let ny = a.sin();
                            let v = push_vert(&mut positions, &mut normals, &mut distances, &mut cap_join_flags,
                                [vx + nx * hw_circ, vy + ny * hw_circ],
                                [nx * CIRCUMSCRIBE, ny * CIRCUMSCRIBE], dist, 1.0);
                            indices.extend_from_slice(&[center_idx, v, prev_v]);
                            prev_v = v;
                        }
                        push_vert(&mut positions, &mut normals, &mut distances, &mut cap_join_flags,
                            [vx + n_next.x * half_width, vy + n_next.y * half_width],
                            [n_next.x, n_next.y], dist, 0.0);
                        push_vert(&mut positions, &mut normals, &mut distances, &mut cap_join_flags,
                            [vx - n_next.x * half_width, vy - n_next.y * half_width],
                            [-n_next.x, -n_next.y], dist, 0.0);
                    }
                } else {
                    // Bevel join (also miter fallback).
                    let l_in = push_vert(&mut positions, &mut normals, &mut distances, &mut cap_join_flags,
                        [vx + n.x * half_width, vy + n.y * half_width],
                        [n.x, n.y], dist, 0.0);
                    let r_in = push_vert(&mut positions, &mut normals, &mut distances, &mut cap_join_flags,
                        [vx - n.x * half_width, vy - n.y * half_width],
                        [-n.x, -n.y], dist, 0.0);
                    indices.extend_from_slice(&[prev_left, prev_right, l_in, prev_right, r_in, l_in]);

                    // Bevel triangle to the next-segment normal.
                    let l2 = push_vert(&mut positions, &mut normals, &mut distances, &mut cap_join_flags,
                        [vx + n_next.x * half_width, vy + n_next.y * half_width],
                        [n_next.x, n_next.y], dist, 0.0);
                    let r2 = push_vert(&mut positions, &mut normals, &mut distances, &mut cap_join_flags,
                        [vx - n_next.x * half_width, vy - n_next.y * half_width],
                        [-n_next.x, -n_next.y], dist, 0.0);
                    indices.extend_from_slice(&[l_in, r_in, l2, r_in, r2, l2]);
                }
            }
        } else {
            // Last vertex — no join needed.  Emit the endpoint pair.
            let l = push_vert(&mut positions, &mut normals, &mut distances, &mut cap_join_flags,
                [vx + n.x * half_width, vy + n.y * half_width],
                [n.x, n.y], dist, 0.0);
            let r = push_vert(&mut positions, &mut normals, &mut distances, &mut cap_join_flags,
                [vx - n.x * half_width, vy - n.y * half_width],
                [-n.x, -n.y], dist, 0.0);
            indices.extend_from_slice(&[prev_left, prev_right, l, prev_right, r, l]);
        }
    }

    // --- End cap ---
    let last_n = seg_normals[seg_count - 1];
    let last_t = seg_tangents[seg_count - 1];
    let ex = coords[coords.len() - 1].lon;
    let ey = coords[coords.len() - 1].lat;
    let end_dist = cum_dist[coords.len() - 1];

    match cap {
        LineCap::Butt => { /* already terminated by the last quad */ }
        LineCap::Square => {
            let prev_left = positions.len() as u32 - 2;
            let prev_right = positions.len() as u32 - 1;
            let fx = ex + last_t.x * half_width;
            let fy = ey + last_t.y * half_width;
            let l = push_vert(&mut positions, &mut normals, &mut distances, &mut cap_join_flags,
                [fx + last_n.x * half_width, fy + last_n.y * half_width],
                [last_n.x, last_n.y], end_dist, 0.0);
            let r = push_vert(&mut positions, &mut normals, &mut distances, &mut cap_join_flags,
                [fx - last_n.x * half_width, fy - last_n.y * half_width],
                [-last_n.x, -last_n.y], end_dist, 0.0);
            indices.extend_from_slice(&[prev_left, prev_right, l, prev_right, r, l]);
        }
        LineCap::Round => {
            let center_idx = push_vert(&mut positions, &mut normals, &mut distances, &mut cap_join_flags,
                [ex, ey], [0.0, 0.0], end_dist, 1.0);
            // Semicircular fan on the forward side with circumscribed radius.
            let hw_circ = half_width * CIRCUMSCRIBE;
            let start_angle = last_n.y.atan2(last_n.x);
            let mut fan_verts = Vec::with_capacity(ROUND_SEGMENTS as usize + 1);
            for k in 0..=ROUND_SEGMENTS {
                let a = start_angle - std::f64::consts::PI * k as f64 / ROUND_SEGMENTS as f64;
                let nx = a.cos();
                let ny = a.sin();
                let v = push_vert(&mut positions, &mut normals, &mut distances, &mut cap_join_flags,
                    [ex + nx * hw_circ, ey + ny * hw_circ],
                    [nx * CIRCUMSCRIBE, ny * CIRCUMSCRIBE], end_dist, 1.0);
                fan_verts.push(v);
            }
            // Connect the fan to the last ribbon pair.
            let prev_left = center_idx - 2;
            let prev_right = center_idx - 1;
            indices.extend_from_slice(&[prev_left, fan_verts[0], center_idx]);
            indices.extend_from_slice(&[prev_right, center_idx, fan_verts[ROUND_SEGMENTS as usize]]);
            for k in 0..ROUND_SEGMENTS {
                indices.push(center_idx);
                indices.push(fan_verts[k as usize]);
                indices.push(fan_verts[k as usize + 1]);
            }
        }
    }

    StrokeLineResult { positions, indices, normals, distances, cap_join: cap_join_flags }
}

// ---------------------------------------------------------------------------
// Tests
// ---------------------------------------------------------------------------

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

    // =====================================================================
    // triangulate_polygon
    // =====================================================================

    #[test]
    fn triangulate_empty() {
        assert!(triangulate_polygon(&[]).is_empty());
    }

    #[test]
    fn triangulate_single_point() {
        assert!(triangulate_polygon(&[GeoCoord::from_lat_lon(0.0, 0.0)]).is_empty());
    }

    #[test]
    fn triangulate_two_points() {
        let coords = vec![
            GeoCoord::from_lat_lon(0.0, 0.0),
            GeoCoord::from_lat_lon(1.0, 1.0),
        ];
        assert!(triangulate_polygon(&coords).is_empty());
    }

    #[test]
    fn triangulate_triangle() {
        let coords = vec![
            GeoCoord::from_lat_lon(0.0, 0.0),
            GeoCoord::from_lat_lon(0.0, 1.0),
            GeoCoord::from_lat_lon(1.0, 0.0),
        ];
        let indices = triangulate_polygon(&coords);
        assert_eq!(indices.len(), 3);
        // All three vertices must appear exactly once.
        let mut sorted = indices.clone();
        sorted.sort();
        assert_eq!(sorted, vec![0, 1, 2]);
    }

    #[test]
    fn triangulate_square() {
        let coords = vec![
            GeoCoord::from_lat_lon(0.0, 0.0),
            GeoCoord::from_lat_lon(0.0, 1.0),
            GeoCoord::from_lat_lon(1.0, 1.0),
            GeoCoord::from_lat_lon(1.0, 0.0),
        ];
        let indices = triangulate_polygon(&coords);
        assert_eq!(indices.len(), 6); // 2 triangles
        assert!(indices.iter().all(|&i| i < 4));
    }

    #[test]
    fn triangulate_closed_polygon() {
        // Closing vertex (duplicate of first) should be stripped.
        let coords = vec![
            GeoCoord::from_lat_lon(0.0, 0.0),
            GeoCoord::from_lat_lon(0.0, 1.0),
            GeoCoord::from_lat_lon(1.0, 1.0),
            GeoCoord::from_lat_lon(1.0, 0.0),
            GeoCoord::from_lat_lon(0.0, 0.0),
        ];
        let indices = triangulate_polygon(&coords);
        assert_eq!(indices.len(), 6); // 2 triangles, closing vertex stripped
        // All indices must reference the first 4 vertices (not the 5th).
        assert!(indices.iter().all(|&i| i < 4));
    }

    #[test]
    fn triangulate_closing_vertex_reduces_below_three() {
        // 3 vertices where last == first -> effective_len = 2 -> degenerate.
        let coords = vec![
            GeoCoord::from_lat_lon(0.0, 0.0),
            GeoCoord::from_lat_lon(1.0, 1.0),
            GeoCoord::from_lat_lon(0.0, 0.0),
        ];
        // n=3, closing detected (but only when n > 3), so no stripping.
        // effective_len stays 3. This produces 1 degenerate triangle --
        // the function does not reject degenerate triangles, which is
        // fine (they render as nothing).
        let indices = triangulate_polygon(&coords);
        assert_eq!(indices.len(), 3);
    }

    #[test]
    fn triangulate_pentagon() {
        let coords = vec![
            GeoCoord::from_lat_lon(0.0, 0.0),
            GeoCoord::from_lat_lon(0.0, 2.0),
            GeoCoord::from_lat_lon(1.0, 3.0),
            GeoCoord::from_lat_lon(2.0, 2.0),
            GeoCoord::from_lat_lon(2.0, 0.0),
        ];
        let indices = triangulate_polygon(&coords);
        // 5 vertices -> 3 triangles -> 9 indices.
        assert_eq!(indices.len(), 9);
        assert!(indices.iter().all(|&i| i < 5));
    }

    #[test]
    fn triangulate_indices_are_valid() {
        // For any N-gon, every index must be in [0, effective_len).
        let coords: Vec<GeoCoord> = (0..20)
            .map(|i| {
                let angle = 2.0 * std::f64::consts::PI * i as f64 / 20.0;
                GeoCoord::from_lat_lon(angle.sin() * 10.0, angle.cos() * 10.0)
            })
            .collect();
        let indices = triangulate_polygon(&coords);
        assert_eq!(indices.len(), 18 * 3); // 20 - 2 = 18 triangles
        assert!(indices.iter().all(|&i| (i as usize) < coords.len()));
    }

    // =====================================================================
    // stroke_line
    // =====================================================================

    #[test]
    fn stroke_empty() {
        let (p, i) = stroke_line(&[], 1.0);
        assert!(p.is_empty());
        assert!(i.is_empty());
    }

    #[test]
    fn stroke_single_point() {
        let (p, i) = stroke_line(&[GeoCoord::from_lat_lon(0.0, 0.0)], 1.0);
        assert!(p.is_empty());
        assert!(i.is_empty());
    }

    #[test]
    fn stroke_line_two_points() {
        let coords = vec![
            GeoCoord::from_lat_lon(0.0, 0.0),
            GeoCoord::from_lat_lon(0.0, 1.0),
        ];
        let (positions, indices) = stroke_line(&coords, 0.01);
        assert_eq!(positions.len(), 4); // 2 vertices * 2 sides
        assert_eq!(indices.len(), 6);   // 1 segment * 6 indices
    }

    #[test]
    fn stroke_line_three_points() {
        let coords = vec![
            GeoCoord::from_lat_lon(0.0, 0.0),
            GeoCoord::from_lat_lon(0.0, 1.0),
            GeoCoord::from_lat_lon(0.0, 2.0),
        ];
        let (positions, indices) = stroke_line(&coords, 0.01);
        assert_eq!(positions.len(), 6); // 3 vertices * 2 sides
        assert_eq!(indices.len(), 12);  // 2 segments * 6 indices
    }

    #[test]
    fn stroke_ribbon_has_nonzero_width() {
        // A horizontal east-west line should be extruded north-south.
        let coords = vec![
            GeoCoord::from_lat_lon(0.0, 0.0),
            GeoCoord::from_lat_lon(0.0, 1.0),
        ];
        let hw = 0.5;
        let (positions, _) = stroke_line(&coords, hw);

        // Vertex 0 (left) and vertex 1 (right) at the first point.
        let left = positions[0];
        let right = positions[1];

        // The extrusion should be in the Y (lat) direction because the
        // tangent is east-west (X = lon direction), so the normal is
        // north-south.
        let dy = (left[1] - right[1]).abs();
        assert!(
            (dy - 2.0 * hw).abs() < 1e-12,
            "expected ribbon width {}, got {dy}",
            2.0 * hw
        );
    }

    #[test]
    fn stroke_indices_are_valid() {
        let coords: Vec<GeoCoord> = (0..10)
            .map(|i| GeoCoord::from_lat_lon(0.0, i as f64))
            .collect();
        let (positions, indices) = stroke_line(&coords, 0.01);
        let max_idx = positions.len() as u32;
        assert!(
            indices.iter().all(|&i| i < max_idx),
            "all indices must be within the position buffer"
        );
    }

    #[test]
    fn stroke_zero_width() {
        // Zero half-width should produce degenerate (zero-area) quads
        // without panicking.
        let coords = vec![
            GeoCoord::from_lat_lon(0.0, 0.0),
            GeoCoord::from_lat_lon(0.0, 1.0),
        ];
        let (positions, indices) = stroke_line(&coords, 0.0);
        assert_eq!(positions.len(), 4);
        assert_eq!(indices.len(), 6);
        // Left and right should coincide.
        assert_eq!(positions[0], positions[1]);
    }

    #[test]
    fn stroke_coincident_points() {
        // Two identical points: tangent is zero-length, normal is zero.
        // Should not panic.
        let coords = vec![
            GeoCoord::from_lat_lon(5.0, 10.0),
            GeoCoord::from_lat_lon(5.0, 10.0),
        ];
        let (positions, indices) = stroke_line(&coords, 0.01);
        assert_eq!(positions.len(), 4);
        assert_eq!(indices.len(), 6);
    }

    // =====================================================================
    // normalize
    // =====================================================================

    #[test]
    fn normalize_unit_vector() {
        let v = normalize(Vec2 { x: 3.0, y: 4.0 });
        let len = (v.x * v.x + v.y * v.y).sqrt();
        assert!((len - 1.0).abs() < 1e-12);
    }

    #[test]
    fn normalize_zero_vector() {
        let v = normalize(Vec2 { x: 0.0, y: 0.0 });
        assert_eq!(v.x, 0.0);
        assert_eq!(v.y, 0.0);
    }

    #[test]
    fn normalize_tiny_vector() {
        // Below the 1e-15 threshold.
        let v = normalize(Vec2 { x: 1e-16, y: 0.0 });
        assert_eq!(v.x, 0.0);
        assert_eq!(v.y, 0.0);
    }

    // =====================================================================
    // earcut -- concave polygons
    // =====================================================================

    #[test]
    fn triangulate_concave_l_shape() {
        // An L-shaped polygon that is concave.
        //
        //  (0,2)---(1,2)
        //    |       |
        //  (0,1)---(1,1)---(2,1)
        //                    |
        //  (0,0)-----------(2,0)
        let coords = vec![
            GeoCoord::from_lat_lon(0.0, 0.0),
            GeoCoord::from_lat_lon(0.0, 2.0),
            GeoCoord::from_lat_lon(1.0, 2.0),
            GeoCoord::from_lat_lon(1.0, 1.0),
            GeoCoord::from_lat_lon(2.0, 1.0),
            GeoCoord::from_lat_lon(2.0, 0.0),
        ];
        let indices = triangulate_polygon(&coords);
        // 6 vertices -> 4 triangles -> 12 indices.
        assert_eq!(indices.len(), 12);
        assert!(indices.iter().all(|&i| i < 6));

        // Verify total signed area matches the expected L-shape area:
        // Full 2x2 square minus 1x1 corner = 3.
        let area: f64 = indices.chunks(3).map(|tri| {
            let a = &coords[tri[0] as usize];
            let b = &coords[tri[1] as usize];
            let c = &coords[tri[2] as usize];
            // Signed area of triangle (positive for CCW, negative for CW).
            ((b.lon - a.lon) * (c.lat - a.lat) - (c.lon - a.lon) * (b.lat - a.lat)) * 0.5
        }).sum::<f64>().abs();
        assert!((area - 3.0).abs() < 1e-6, "L-shape area should be 3.0, got {area}");
    }

    // =====================================================================
    // earcut -- polygon with holes
    // =====================================================================

    #[test]
    fn triangulate_with_hole() {
        let exterior = vec![
            GeoCoord::from_lat_lon(0.0, 0.0),
            GeoCoord::from_lat_lon(0.0, 4.0),
            GeoCoord::from_lat_lon(4.0, 4.0),
            GeoCoord::from_lat_lon(4.0, 0.0),
        ];
        let hole = vec![
            GeoCoord::from_lat_lon(1.0, 1.0),
            GeoCoord::from_lat_lon(1.0, 3.0),
            GeoCoord::from_lat_lon(3.0, 3.0),
            GeoCoord::from_lat_lon(3.0, 1.0),
        ];
        let indices = triangulate_polygon_with_holes(&exterior, &[&hole]);
        // 4 exterior + 4 hole = 8 vertices -> 6 triangles -> 18 indices.
        assert_eq!(indices.len() % 3, 0);
        assert!(!indices.is_empty());
        // All indices must reference valid vertices.
        assert!(indices.iter().all(|&i| (i as usize) < 8));
    }

    #[test]
    fn triangulate_with_hole_empty_exterior() {
        let exterior: Vec<GeoCoord> = Vec::new();
        let hole = vec![
            GeoCoord::from_lat_lon(1.0, 1.0),
            GeoCoord::from_lat_lon(1.0, 3.0),
            GeoCoord::from_lat_lon(3.0, 3.0),
        ];
        assert!(triangulate_polygon_with_holes(&exterior, &[&hole]).is_empty());
    }

    // =====================================================================
    // stroke_line_styled
    // =====================================================================

    #[test]
    fn styled_stroke_empty() {
        let r = stroke_line_styled(&[], 1.0, LineCap::Butt, LineJoin::Miter, 2.0);
        assert!(r.positions.is_empty());
        assert!(r.indices.is_empty());
        assert!(r.normals.is_empty());
        assert!(r.distances.is_empty());
    }

    #[test]
    fn styled_stroke_single_point() {
        let r = stroke_line_styled(
            &[GeoCoord::from_lat_lon(0.0, 0.0)],
            1.0, LineCap::Butt, LineJoin::Miter, 2.0,
        );
        assert!(r.positions.is_empty());
    }

    #[test]
    fn styled_stroke_butt_cap_two_points() {
        let coords = vec![
            GeoCoord::from_lat_lon(0.0, 0.0),
            GeoCoord::from_lat_lon(0.0, 1.0),
        ];
        let r = stroke_line_styled(&coords, 0.01, LineCap::Butt, LineJoin::Miter, 2.0);
        // 2 start-cap verts + 2 end-segment verts = 4
        assert_eq!(r.positions.len(), 4);
        assert_eq!(r.normals.len(), 4);
        assert_eq!(r.distances.len(), 4);
        assert!(!r.indices.is_empty());
        // All indices valid.
        let max_idx = r.positions.len() as u32;
        assert!(r.indices.iter().all(|&i| i < max_idx));
    }

    #[test]
    fn styled_stroke_square_cap_extends_beyond_endpoints() {
        let coords = vec![
            GeoCoord::from_lat_lon(0.0, 0.0),
            GeoCoord::from_lat_lon(0.0, 1.0),
        ];
        let hw = 0.1;
        let r = stroke_line_styled(&coords, hw, LineCap::Square, LineJoin::Miter, 2.0);
        // Square caps extend half_width backwards, so min lon should be < 0.
        let min_lon: f64 = r.positions.iter().map(|p| p[0]).fold(f64::INFINITY, f64::min);
        assert!(min_lon < -0.05, "square cap should extend before the first vertex, got {min_lon}");
    }

    #[test]
    fn styled_stroke_round_cap_produces_more_vertices() {
        let coords = vec![
            GeoCoord::from_lat_lon(0.0, 0.0),
            GeoCoord::from_lat_lon(0.0, 1.0),
        ];
        let butt = stroke_line_styled(&coords, 0.01, LineCap::Butt, LineJoin::Miter, 2.0);
        let round = stroke_line_styled(&coords, 0.01, LineCap::Round, LineJoin::Miter, 2.0);
        assert!(
            round.positions.len() > butt.positions.len(),
            "round cap should produce more vertices than butt: {} vs {}",
            round.positions.len(), butt.positions.len(),
        );
    }

    #[test]
    fn styled_stroke_distances_are_monotonic() {
        let coords = vec![
            GeoCoord::from_lat_lon(0.0, 0.0),
            GeoCoord::from_lat_lon(0.0, 1.0),
            GeoCoord::from_lat_lon(0.0, 2.0),
        ];
        let r = stroke_line_styled(&coords, 0.01, LineCap::Butt, LineJoin::Miter, 2.0);
        // Distances should be non-negative and the max should be > 0.
        assert!(r.distances.iter().all(|&d| d >= 0.0));
        let max_dist = r.distances.iter().cloned().fold(0.0_f64, f64::max);
        assert!(max_dist > 0.0, "max distance should be positive");
    }

    #[test]
    fn styled_stroke_normals_are_unit_length_or_zero() {
        let coords = vec![
            GeoCoord::from_lat_lon(0.0, 0.0),
            GeoCoord::from_lat_lon(0.0, 1.0),
            GeoCoord::from_lat_lon(1.0, 1.0),
        ];
        let r = stroke_line_styled(&coords, 0.01, LineCap::Round, LineJoin::Round, 2.0);
        for (i, nrm) in r.normals.iter().enumerate() {
            let len = (nrm[0] * nrm[0] + nrm[1] * nrm[1]).sqrt();
            if r.cap_join[i] > 0.5 {
                // Circumscribed cap/join fan vertices: magnitude ≈ CIRCUMSCRIBE or 0 (center).
                assert!(
                    len < CIRCUMSCRIBE + 1e-10 || len < 1e-10,
                    "cap/join normal length should be ≤{CIRCUMSCRIBE} or ~0, got {len}"
                );
            } else {
                // Ribbon body vertices: unit length.
                assert!(
                    len < 1.0 + 1e-10 || len < 1e-10,
                    "body normal length should be ≤1 or ~0, got {len}"
                );
            }
        }
    }

    #[test]
    fn styled_stroke_cap_join_flags_correct_for_round_caps() {
        let coords = vec![
            GeoCoord::from_lat_lon(0.0, 0.0),
            GeoCoord::from_lat_lon(0.0, 1.0),
        ];
        let r = stroke_line_styled(&coords, 0.01, LineCap::Round, LineJoin::Round, 2.0);
        assert_eq!(r.cap_join.len(), r.positions.len());
        // At least some vertices should be flagged as cap/join.
        let cap_join_count = r.cap_join.iter().filter(|&&f| f > 0.5).count();
        assert!(cap_join_count > 0, "round cap should flag some vertices");
        // Ribbon body vertices should exist and be flagged 0.0.
        let body_count = r.cap_join.iter().filter(|&&f| f < 0.5).count();
        assert!(body_count > 0, "ribbon body vertices should exist");
    }

    #[test]
    fn styled_stroke_cap_join_flags_zero_for_butt_caps() {
        let coords = vec![
            GeoCoord::from_lat_lon(0.0, 0.0),
            GeoCoord::from_lat_lon(0.0, 1.0),
        ];
        let r = stroke_line_styled(&coords, 0.01, LineCap::Butt, LineJoin::Bevel, 2.0);
        assert_eq!(r.cap_join.len(), r.positions.len());
        // No cap/join vertices for butt caps + bevel joins.
        for &flag in &r.cap_join {
            assert!(flag < 0.5, "butt cap + bevel join should have no cap_join flags");
        }
    }

    #[test]
    fn styled_stroke_circumscribed_fan_extends_beyond_unit() {
        // Round cap fan perimeter normals should have magnitude > 1.0
        // (circumscribed) for SDF clipping.
        let coords = vec![
            GeoCoord::from_lat_lon(0.0, 0.0),
            GeoCoord::from_lat_lon(0.0, 1.0),
        ];
        let r = stroke_line_styled(&coords, 0.01, LineCap::Round, LineJoin::Round, 2.0);
        let mut found_circumscribed = false;
        for (i, nrm) in r.normals.iter().enumerate() {
            if r.cap_join[i] > 0.5 {
                let len = (nrm[0] * nrm[0] + nrm[1] * nrm[1]).sqrt();
                if len > 1.0 + 1e-10 {
                    found_circumscribed = true;
                    assert!(
                        (len - CIRCUMSCRIBE).abs() < 1e-6,
                        "circumscribed normal should be ~{CIRCUMSCRIBE}, got {len}"
                    );
                }
            }
        }
        assert!(found_circumscribed, "should find circumscribed fan normals with length > 1");
    }

    #[test]
    fn styled_stroke_bevel_join_at_sharp_angle() {
        // 90-degree turn: miter_limit 1.0 forces bevel.
        let coords = vec![
            GeoCoord::from_lat_lon(0.0, 0.0),
            GeoCoord::from_lat_lon(0.0, 1.0),
            GeoCoord::from_lat_lon(1.0, 1.0),
        ];
        let r = stroke_line_styled(&coords, 0.01, LineCap::Butt, LineJoin::Bevel, 2.0);
        // Bevel should produce more vertices than a straight line.
        let straight = stroke_line_styled(
            &[GeoCoord::from_lat_lon(0.0, 0.0), GeoCoord::from_lat_lon(0.0, 2.0)],
            0.01, LineCap::Butt, LineJoin::Bevel, 2.0,
        );
        assert!(
            r.positions.len() > straight.positions.len(),
            "bevel join should add extra vertices: {} vs {}",
            r.positions.len(), straight.positions.len(),
        );
    }

    #[test]
    fn styled_stroke_round_join_produces_fan_vertices() {
        let coords = vec![
            GeoCoord::from_lat_lon(0.0, 0.0),
            GeoCoord::from_lat_lon(0.0, 1.0),
            GeoCoord::from_lat_lon(1.0, 1.0),
        ];
        let bevel = stroke_line_styled(&coords, 0.01, LineCap::Butt, LineJoin::Bevel, 2.0);
        let round = stroke_line_styled(&coords, 0.01, LineCap::Butt, LineJoin::Round, 2.0);
        assert!(
            round.positions.len() > bevel.positions.len(),
            "round join should produce more vertices than bevel: {} vs {}",
            round.positions.len(), bevel.positions.len(),
        );
    }

    #[test]
    fn styled_stroke_miter_join_collinear() {
        // Collinear points: miter should be simple (no extra vertices).
        let coords = vec![
            GeoCoord::from_lat_lon(0.0, 0.0),
            GeoCoord::from_lat_lon(0.0, 1.0),
            GeoCoord::from_lat_lon(0.0, 2.0),
        ];
        let r = stroke_line_styled(&coords, 0.01, LineCap::Butt, LineJoin::Miter, 2.0);
        // All indices valid.
        let max_idx = r.positions.len() as u32;
        assert!(r.indices.iter().all(|&i| i < max_idx));
        assert!(!r.indices.is_empty());
    }

    #[test]
    fn styled_stroke_indices_valid_for_complex_line() {
        // 10-point zigzag.
        let coords: Vec<GeoCoord> = (0..10)
            .map(|i| {
                let lat = if i % 2 == 0 { 0.0 } else { 1.0 };
                GeoCoord::from_lat_lon(lat, i as f64)
            })
            .collect();
        for cap in [LineCap::Butt, LineCap::Round, LineCap::Square] {
            for join in [LineJoin::Miter, LineJoin::Bevel, LineJoin::Round] {
                let r = stroke_line_styled(&coords, 0.01, cap, join, 2.0);
                let max_idx = r.positions.len() as u32;
                assert!(
                    r.indices.iter().all(|&i| i < max_idx),
                    "invalid index for cap={cap:?} join={join:?}: max valid={max_idx}",
                );
                assert_eq!(r.positions.len(), r.normals.len());
                assert_eq!(r.positions.len(), r.distances.len());
            }
        }
    }
}