TheriapolisV3/Theriapolis.Core/World/Polylines/PolylineBuilder.cs

using Theriapolis.Core.Util;
using Theriapolis.Core.World;

namespace Theriapolis.Core.World.Polylines;

/// <summary>
/// Static helpers for building, smoothing, and rasterizing polylines in world-pixel space.
/// </summary>
public static class PolylineBuilder
{
    // ── Catmull-Rom spline smoothing ─────────────────────────────────────────

    /// <summary>
    /// Apply Catmull-Rom spline smoothing to a list of control points.
    /// Produces <see cref="C.SPLINE_SUBDIVISIONS"/> intermediate points per segment.
    /// </summary>
    public static List<Vec2> CatmullRomSmooth(IReadOnlyList<Vec2> points, int subdivisions = C.SPLINE_SUBDIVISIONS)
    {
        var result = new List<Vec2>(points.Count * subdivisions);
        if (points.Count < 2) { result.AddRange(points); return result; }

        for (int i = 0; i < points.Count - 1; i++)
        {
            var p0 = points[Math.Max(0, i - 1)];
            var p1 = points[i];
            var p2 = points[i + 1];
            var p3 = points[Math.Min(points.Count - 1, i + 2)];

            for (int s = 0; s < subdivisions; s++)
            {
                float t = s / (float)subdivisions;
                result.Add(CatmullRomPoint(p0, p1, p2, p3, t));
            }
        }
        result.Add(points[^1]);
        return result;
    }

    private static Vec2 CatmullRomPoint(Vec2 p0, Vec2 p1, Vec2 p2, Vec2 p3, float t)
    {
        float t2 = t * t;
        float t3 = t2 * t;
        float x = 0.5f * ((2f * p1.X)
                       + (-p0.X + p2.X) * t
                       + (2f * p0.X - 5f * p1.X + 4f * p2.X - p3.X) * t2
                       + (-p0.X + 3f * p1.X - 3f * p2.X + p3.X) * t3);
        float y = 0.5f * ((2f * p1.Y)
                       + (-p0.Y + p2.Y) * t
                       + (2f * p0.Y - 5f * p1.Y + 4f * p2.Y - p3.Y) * t2
                       + (-p0.Y + 3f * p1.Y - 3f * p2.Y + p3.Y) * t3);
        return new Vec2(x, y);
    }

    // ── Perpendicular meander noise ──────────────────────────────────────────

    /// <summary>
    /// Apply a perpendicular noise offset to each point along a smoothed polyline.
    /// Uses FastNoiseLite for deterministic, seeded noise.
    /// </summary>
    public static void ApplyMeanderNoise(
        List<Vec2> points,
        float amplitude,
        float frequency,
        ulong noiseSeed)
    {
        if (points.Count < 2) return;

        // Build a FastNoiseLite instance seeded for this polyline
        var noise = new FastNoiseLite
        {
            Seed      = (int)(noiseSeed & 0x7FFFFFFFu),
            Noise     = FastNoiseLite.NoiseType.OpenSimplex2,
            Frequency = frequency,
        };

        for (int i = 0; i < points.Count; i++)
        {
            // Perpendicular direction at this point
            Vec2 tangent;
            if (i == 0)
                tangent = (points[1] - points[0]).Normalized;
            else if (i == points.Count - 1)
                tangent = (points[^1] - points[^2]).Normalized;
            else
                tangent = (points[i + 1] - points[i - 1]).Normalized;

            Vec2 perp = tangent.Perp;

            // Sample noise at world-pixel position
            float offset = noise.GetNoise(points[i].X * 0.01f, points[i].Y * 0.01f) * amplitude;
            points[i] = points[i] + perp * offset;
        }
    }

    // ── Ramer-Douglas-Peucker simplification ─────────────────────────────────

    /// <summary>
    /// Ramer-Douglas-Peucker polyline simplification for LOD rendering.
    /// Returns a simplified point list with perpendicular distance tolerance in world pixels.
    /// </summary>
    public static List<Vec2> RDPSimplify(IReadOnlyList<Vec2> points, float tolerance = C.RDP_TOLERANCE)
    {
        if (points.Count <= 2) return new List<Vec2>(points);

        var result = new List<Vec2>();
        var stack  = new Stack<(int start, int end)>();
        var keep   = new bool[points.Count];
        keep[0] = keep[^1] = true;

        stack.Push((0, points.Count - 1));
        while (stack.Count > 0)
        {
            var (start, end) = stack.Pop();
            float maxDist  = 0f;
            int   maxIndex = start;

            for (int i = start + 1; i < end; i++)
            {
                float d = PerpendicularDistance(points[i], points[start], points[end]);
                if (d > maxDist) { maxDist = d; maxIndex = i; }
            }

            if (maxDist > tolerance)
            {
                keep[maxIndex] = true;
                stack.Push((start, maxIndex));
                stack.Push((maxIndex, end));
            }
        }

        for (int i = 0; i < points.Count; i++)
            if (keep[i]) result.Add(points[i]);

        return result;
    }

    private static float PerpendicularDistance(Vec2 pt, Vec2 lineStart, Vec2 lineEnd)
    {
        Vec2 d = lineEnd - lineStart;
        float len = d.Length;
        if (len < 1e-6f) return Vec2.Dist(pt, lineStart);
        return MathF.Abs(Vec2.Dot(d.Perp, pt - lineStart)) / len;
    }

    // ── Tile flag rasterization ──────────────────────────────────────────────

    /// <summary>
    /// Rasterizes a polyline (in world-pixel space) onto the world tile grid,
    /// setting the appropriate FeatureFlags and direction on each touched tile and its neighbors.
    /// </summary>
    public static void RasterizeToTileFlags(
        WorldState  world,
        Polyline    polyline,
        bool        setAdjacency = true)
    {
        var pts = polyline.Points;
        if (pts.Count < 2) return;

        int px = C.WORLD_TILE_PIXELS;

        for (int seg = 0; seg < pts.Count - 1; seg++)
        {
            int x0 = (int)(pts[seg].X / px);
            int y0 = (int)(pts[seg].Y / px);
            int x1 = (int)(pts[seg + 1].X / px);
            int y1 = (int)(pts[seg + 1].Y / px);

            // Direction of this segment
            int ddx = Math.Sign(x1 - x0);
            int ddy = Math.Sign(y1 - y0);
            byte dir = Dir.FromDelta(ddx == 0 && ddy == 0 ? 0 : ddx, ddy);

            // Bresenham line in tile space
            BresenhamLine(x0, y0, x1, y1, (tx, ty) =>
            {
                if ((uint)tx >= C.WORLD_WIDTH_TILES || (uint)ty >= C.WORLD_HEIGHT_TILES) return;
                if (world.TileAt(tx, ty).Biome == BiomeId.Ocean) return;
                ref var tile = ref world.TileAt(tx, ty);

                switch (polyline.Type)
                {
                    case PolylineType.River:
                        tile.Features |= FeatureFlags.HasRiver;
                        if (tile.RiverFlowDir == Dir.None) tile.RiverFlowDir = dir;
                        break;
                    case PolylineType.Road:
                        // Allow HasRoad on river tiles (bridge crossings) so that
                        // SplitByExistingFeature treats shared crossings as existing road
                        // and doesn't create redundant short segments.
                        // Still skip rail tiles (Addendum A §2).
                        if ((tile.Features & FeatureFlags.HasRail) == 0 ||
                            (tile.Features & FeatureFlags.IsSettlement) != 0)
                            tile.Features |= FeatureFlags.HasRoad;
                        break;
                    case PolylineType.Rail:
                        // Bridges: don't mark a river tile as also having rail.
                        // The crossing is implied by the route; co-location would be a violation.
                        if ((tile.Features & FeatureFlags.HasRiver) == 0 ||
                            (tile.Features & FeatureFlags.IsSettlement) != 0)
                        {
                            tile.Features |= FeatureFlags.HasRail;
                            if (tile.RailDir == Dir.None) tile.RailDir = dir;
                        }
                        break;
                }

                if (!setAdjacency) return;

                // Mark 8 neighbors as adjacent
                for (int ny = ty - 1; ny <= ty + 1; ny++)
                for (int nx = tx - 1; nx <= tx + 1; nx++)
                {
                    if (nx == tx && ny == ty) continue;
                    if ((uint)nx >= C.WORLD_WIDTH_TILES || (uint)ny >= C.WORLD_HEIGHT_TILES) continue;
                    if (world.TileAt(nx, ny).Biome == BiomeId.Ocean) continue;
                    ref var neighbor = ref world.TileAt(nx, ny);

                    switch (polyline.Type)
                    {
                        case PolylineType.River:
                            neighbor.Features |= FeatureFlags.RiverAdjacent;
                            if (neighbor.RiverFlowDir == Dir.None) neighbor.RiverFlowDir = dir;
                            break;
                        case PolylineType.Rail:
                            neighbor.Features |= FeatureFlags.RailroadAdjacent;
                            if (neighbor.RailDir == Dir.None) neighbor.RailDir = dir;
                            break;
                    }
                }
            });
        }
    }

    private static void BresenhamLine(int x0, int y0, int x1, int y1, Action<int, int> visit)
    {
        int dx = Math.Abs(x1 - x0);
        int dy = Math.Abs(y1 - y0);
        int sx = x0 < x1 ? 1 : -1;
        int sy = y0 < y1 ? 1 : -1;
        int err = dx - dy;

        while (true)
        {
            visit(x0, y0);
            if (x0 == x1 && y0 == y1) break;
            int e2 = 2 * err;
            if (e2 > -dy) { err -= dy; x0 += sx; }
            if (e2 <  dx) { err += dx; y0 += sy; }
        }
    }

    // ── Staircase smoothing ───────────────────────────────────────────────────

    /// <summary>
    /// Post-processes an A* tile path to reduce staircase artefacts in the
    /// final smoothed polyline. A* tiebreaking can produce sequences like
    /// <c>SE, E, SE, SE, SE</c> where a cardinal move lands in the middle of
    /// runs of diagonal moves — Catmull-Rom smoothing then develops a visible
    /// kink.
    ///
    /// This walks the path looking for (cardinal, diagonal) move pairs and
    /// swaps them to (diagonal, cardinal) when the alternative intermediate
    /// tile is free of preserve-flag features and remains passable under
    /// <paramref name="costFn"/>. Leaves total movement cost unchanged (the
    /// two legs just get reordered) so A*'s optimality choice is respected.
    /// Bubbling is forward-only, so any remaining cardinals end up at the
    /// end of diagonal runs.
    /// </summary>
    public static void SmoothStaircases(
        WorldState world,
        List<(int X, int Y)> path,
        FeatureFlags preserveMask,
        Func<int, int, int, int, byte, float> costFn)
    {
        if (path.Count < 3) return;

        bool changed = true;
        int guard = path.Count; // bound passes to prevent oscillation
        while (changed && guard-- > 0)
        {
            changed = false;
            for (int i = 0; i < path.Count - 2; i++)
            {
                var a = path[i];
                var b = path[i + 1];
                var c = path[i + 2];

                // Preserve tiles with load-bearing features (existing road/rail/
                // river bridges, settlement footprint). Those are deliberate A*
                // choices and junction/endpoint anchors downstream rely on them.
                var bFeatures = world.Tiles[b.X, b.Y].Features;
                if ((bFeatures & preserveMask) != 0) continue;

                int dx1 = b.X - a.X, dy1 = b.Y - a.Y;
                int dx2 = c.X - b.X, dy2 = c.Y - b.Y;

                bool isCard1 = (dx1 == 0) ^ (dy1 == 0);
                bool isDiag2 = (dx2 != 0) && (dy2 != 0);
                if (!isCard1 || !isDiag2) continue;

                // Alternative intermediate: do the diagonal first, then cardinal.
                int nbx = a.X + dx2, nby = a.Y + dy2;
                if ((uint)nbx >= C.WORLD_WIDTH_TILES || (uint)nby >= C.WORLD_HEIGHT_TILES) continue;

                // Don't swap into a tile whose features A* would have weighed
                // differently — stays equivalent in character to the old b.
                if ((world.Tiles[nbx, nby].Features & preserveMask) != 0) continue;

                byte dirToNb = Dir.FromDelta(dx2, dy2);
                byte dirToC  = Dir.FromDelta(dx1, dy1);
                if (float.IsPositiveInfinity(costFn(a.X, a.Y, nbx, nby, dirToNb))) continue;
                if (float.IsPositiveInfinity(costFn(nbx, nby, c.X, c.Y, dirToC)))  continue;

                path[i + 1] = (nbx, nby);
                changed = true;
            }
        }
    }

    // ── Turn-angle limiter ────────────────────────────────────────────────────

    /// <summary>
    /// Caps the deflection angle at every internal vertex of an A* tile path.
    /// A vertex whose turn (angle between incoming and outgoing edges) exceeds
    /// <paramref name="maxTurnDegrees"/> is elided when the shortcut from its
    /// predecessor to successor is a single passable A* step under
    /// <paramref name="costFn"/>. Vertices carrying any flag in
    /// <paramref name="preserveMask"/> are never removed — those are load-bearing
    /// (existing rail/road junctions, river-bridge crossings, settlement footprint),
    /// matching the contract used by <see cref="SmoothStaircases"/>.
    ///
    /// Used by rail generation to avoid 90°/135° bends that look unrealistic
    /// for heavy rail traffic. The grid produces turns in 45° steps, so a cap
    /// of 75° permits 0°/45° turns and forces 90°/135° corners to be smoothed
    /// into two consecutive diagonals.
    /// </summary>
    public static void LimitTurnAngle(
        WorldState world,
        List<(int X, int Y)> path,
        float maxTurnDegrees,
        FeatureFlags preserveMask,
        Func<int, int, int, int, byte, float> costFn)
    {
        if (path.Count < 3) return;
        float minCosAllowed = MathF.Cos(maxTurnDegrees * MathF.PI / 180f);

        bool changed = true;
        int guard = path.Count * 2;
        while (changed && guard-- > 0)
        {
            changed = false;
            for (int i = 1; i < path.Count - 1; i++)
            {
                var a = path[i - 1];
                var b = path[i];
                var c = path[i + 1];

                float vx1 = b.X - a.X, vy1 = b.Y - a.Y;
                float vx2 = c.X - b.X, vy2 = c.Y - b.Y;
                float l1 = MathF.Sqrt(vx1 * vx1 + vy1 * vy1);
                float l2 = MathF.Sqrt(vx2 * vx2 + vy2 * vy2);
                if (l1 < 1e-6f || l2 < 1e-6f) continue;
                float cosTurn = (vx1 * vx2 + vy1 * vy2) / (l1 * l2);
                if (cosTurn >= minCosAllowed) continue; // already gentle enough

                // Keep load-bearing vertices — junctions, bridges, endpoints.
                if ((world.Tiles[b.X, b.Y].Features & preserveMask) != 0) continue;

                // Shortcut a→c must be a single-step move (Chebyshev 1) and
                // remain passable. Larger gaps would span multi-tile stretches
                // the pathfinder never evaluated.
                int dx = c.X - a.X, dy = c.Y - a.Y;
                if (dx == 0 && dy == 0) continue;
                if (Math.Max(Math.Abs(dx), Math.Abs(dy)) != 1) continue;
                byte dir = Dir.FromDelta(Math.Sign(dx), Math.Sign(dy));
                if (float.IsPositiveInfinity(costFn(a.X, a.Y, c.X, c.Y, dir))) continue;

                path.RemoveAt(i);
                changed = true;
                break; // restart — indices shifted
            }
        }
    }

    // ── Path splitting for shared infrastructure ──────────────────────────────

    /// <summary>
    /// Splits a tile path into sub-paths covering only "new construction" — runs of tiles
    /// that do NOT already have the given feature flag. Each sub-path includes one overlap
    /// tile at each junction so the smoothed polyline visually connects to the existing feature.
    /// Returns an empty list if the entire path is already covered.
    ///
    /// </summary>
    public static List<List<(int X, int Y)>> SplitByExistingFeature(
        WorldState world,
        List<(int X, int Y)> path,
        FeatureFlags feature)
    {
        var segments = new List<List<(int X, int Y)>>();
        List<(int X, int Y)>? current = null;

        for (int i = 0; i < path.Count; i++)
        {
            var (x, y) = path[i];
            var tileFeatures = world.Tiles[x, y].Features;
            bool exists = (tileFeatures & feature) != 0;

            if (!exists)
            {
                if (current == null)
                {
                    current = new List<(int X, int Y)>();
                    // Include the last "existing" tile as a junction anchor
                    if (i > 0) current.Add(path[i - 1]);
                }
                current.Add(path[i]);
            }
            else if (current != null)
            {
                // Ending a new-construction run — include this existing tile as endpoint anchor
                current.Add(path[i]);
                if (current.Count >= 2) segments.Add(current);
                current = null;
            }
        }

        // Trailing new-construction segment (path ends on new tiles)
        if (current != null && current.Count >= 2)
            segments.Add(current);

        return segments;
    }

    // ── Control point from tile coords ───────────────────────────────────────

    /// <summary>Convert world tile coordinate to world-pixel center.</summary>
    public static Vec2 TileToWorldPixel(int tileX, int tileY)
    {
        float px = C.WORLD_TILE_PIXELS;
        return new Vec2(tileX * px + px * 0.5f, tileY * px + px * 0.5f);
    }
}