import type { Feature, MultiLineString, Position } from 'geojson';

/**
 * Uniform-grid spatial index over the PI boundary rings.
 *
 * Why this exists: every boundary question we ask during drawing ("is this
 * point inside PI?", "where does this stroke meet the boundary?", "trace the
 * boundary from A to B") used to be answered by scanning the whole ring with
 * turf's `nearestPointOnLine` -- O(n) per query over tens of thousands of
 * vertices, allocating a Feature per call. On the finish path that added up to
 * ~2.3 s of blocked main thread, which the browser reported to the user as an
 * "unresponsive page" dialog.
 *
 * The grid buckets each boundary SEGMENT into every cell its bbox touches.
 * A nearest-point query then only tests segments in the cells near the query
 * point, expanding the search ring until the best candidate is provably closer
 * than the next unexplored ring.
 *
 * Distances use an equirectangular approximation (lng scaled by cos(lat)),
 * which is accurate to well under a metre across the PI latitude band and far
 * cheaper than great-circle maths.
 */

/** Result of snapping an arbitrary coord onto the boundary. */
export interface BoundaryHit {
  /** Distance from the query point to the boundary, in kilometres. */
  distKm: number;
  /** The closest point ON the boundary. */
  snapped: Position;
  /** Which ring of the boundary the hit landed on. */
  ringIdx: number;
  /** Index of the segment within that ring (segment i spans vertex i..i+1). */
  segIdx: number;
  /** Parametric position along that segment, 0..1. */
  t: number;
}

const CELL_DEG = 0.02; // ~2 km at PI latitudes
const M_PER_DEG_LAT = 111_320;
/** cos(31.5 deg) -- PI sits in a narrow latitude band, so one factor is enough. */
const LNG_SCALE = Math.cos((31.5 * Math.PI) / 180);

function cellKey(cx: number, cy: number): string {
  return `${cx}|${cy}`;
}

/** Squared distance in "scaled degrees", where x is pre-multiplied by LNG_SCALE. */
function distSqScaled(ax: number, ay: number, bx: number, by: number): number {
  const dx = (ax - bx) * LNG_SCALE;
  const dy = ay - by;
  return dx * dx + dy * dy;
}

export class BoundaryIndex {
  /** The boundary rings, in original (undensified) resolution. */
  readonly rings: Position[][];
  private readonly grid = new Map<string, number[]>();
  /** Packed segment records: ringIdx, segIdx per entry index. */
  private readonly segRing: number[] = [];
  private readonly segIdx: number[] = [];

  constructor(rings: Position[][]) {
    this.rings = rings;
    for (let r = 0; r < rings.length; r++) {
      const ring = rings[r];
      for (let i = 0; i + 1 < ring.length; i++) {
        const [ax, ay] = ring[i];
        const [bx, by] = ring[i + 1];
        const entry = this.segRing.length;
        this.segRing.push(r);
        this.segIdx.push(i);

        const minX = Math.min(ax, bx);
        const maxX = Math.max(ax, bx);
        const minY = Math.min(ay, by);
        const maxY = Math.max(ay, by);
        const cx0 = Math.floor(minX / CELL_DEG);
        const cx1 = Math.floor(maxX / CELL_DEG);
        const cy0 = Math.floor(minY / CELL_DEG);
        const cy1 = Math.floor(maxY / CELL_DEG);
        for (let cx = cx0; cx <= cx1; cx++) {
          for (let cy = cy0; cy <= cy1; cy++) {
            const k = cellKey(cx, cy);
            const bucket = this.grid.get(k);
            if (bucket) bucket.push(entry);
            else this.grid.set(k, [entry]);
          }
        }
      }
    }
  }

  /** Total segment count -- useful for diagnostics. */
  get segmentCount(): number {
    return this.segRing.length;
  }

  /**
   * Closest point on the boundary to `coord`.
   *
   * Expands the searched cell ring outward until the best hit found so far is
   * closer than the nearest possible point in the next unexplored ring, so the
   * result is exact, not approximate. Returns null only for an empty index.
   */
  nearest(coord: Position): BoundaryHit | null {
    if (this.segRing.length === 0) return null;
    const [qx, qy] = coord;
    const cx = Math.floor(qx / CELL_DEG);
    const cy = Math.floor(qy / CELL_DEG);

    let bestSq = Infinity;
    let best: BoundaryHit | null = null;
    const seen = new Set<number>();

    for (let radius = 0; radius < 4096; radius++) {
      // Once the best hit is closer than the inner edge of this ring of cells,
      // nothing further out can beat it.
      if (best) {
        const ringInnerDeg = (radius - 1) * CELL_DEG;
        if (ringInnerDeg > 0 && bestSq < ringInnerDeg * ringInnerDeg * LNG_SCALE * LNG_SCALE) {
          break;
        }
      }
      let anyCell = false;
      for (let dx = -radius; dx <= radius; dx++) {
        for (let dy = -radius; dy <= radius; dy++) {
          // Only the perimeter of the square is new on each iteration.
          if (radius > 0 && Math.abs(dx) !== radius && Math.abs(dy) !== radius) continue;
          const bucket = this.grid.get(cellKey(cx + dx, cy + dy));
          if (!bucket) continue;
          anyCell = true;
          for (const entry of bucket) {
            if (seen.has(entry)) continue;
            seen.add(entry);
            const r = this.segRing[entry];
            const i = this.segIdx[entry];
            const ring = this.rings[r];
            const hit = pointToSegment(qx, qy, ring[i], ring[i + 1]);
            if (hit.distSq < bestSq) {
              bestSq = hit.distSq;
              best = {
                distKm: 0, // filled in below
                snapped: hit.snapped,
                ringIdx: r,
                segIdx: i,
                t: hit.t,
              };
            }
          }
        }
      }
      // Keep expanding past empty regions until we have at least one candidate.
      if (!anyCell && !best && radius > 512) break;
    }

    if (!best) return null;
    best.distKm = (Math.sqrt(bestSq) * M_PER_DEG_LAT) / 1000;
    return best;
  }

  /** True when `coord` lies within `km` of the boundary. */
  isWithinKm(coord: Position, km: number): boolean {
    const hit = this.nearest(coord);
    return hit !== null && hit.distKm <= km;
  }

  /**
   * First point at which segment a->b crosses the boundary, or null.
   *
   * Only segments in the grid cells the query segment passes through are
   * tested, so this replaces a full-boundary `turf.lineIntersect` that ran on
   * every mousemove while the cursor sat outside PI.
   *
   * When several crossings exist the one nearest `b` is returned -- that is the
   * last place the stroke left the region.
   */
  intersectSegment(a: Position, b: Position): Position | null {
    const minX = Math.min(a[0], b[0]);
    const maxX = Math.max(a[0], b[0]);
    const minY = Math.min(a[1], b[1]);
    const maxY = Math.max(a[1], b[1]);
    const cx0 = Math.floor(minX / CELL_DEG);
    const cx1 = Math.floor(maxX / CELL_DEG);
    const cy0 = Math.floor(minY / CELL_DEG);
    const cy1 = Math.floor(maxY / CELL_DEG);

    let best: Position | null = null;
    let bestSq = Infinity;
    const seen = new Set<number>();
    for (let cx = cx0; cx <= cx1; cx++) {
      for (let cy = cy0; cy <= cy1; cy++) {
        const bucket = this.grid.get(cellKey(cx, cy));
        if (!bucket) continue;
        for (const entry of bucket) {
          if (seen.has(entry)) continue;
          seen.add(entry);
          const ring = this.rings[this.segRing[entry]];
          const i = this.segIdx[entry];
          const hit = segmentIntersection(a, b, ring[i], ring[i + 1]);
          if (!hit) continue;
          const d = distSqScaled(hit[0], hit[1], b[0], b[1]);
          if (d < bestSq) {
            bestSq = d;
            best = hit;
          }
        }
      }
    }
    return best;
  }

  /**
   * Trace the boundary ring from one hit to another, returning the vertex path
   * INCLUSIVE of both snapped endpoints.
   *
   * Walks the ring by vertex index in both directions and returns the shorter
   * path. This is the part turf's `lineSlice` cannot do: a closed ring stored
   * as an array has a seam at index 0, and `lineSlice` will not walk across it,
   * so for any pair of points straddling the seam it silently returns the long
   * way round (up to the entire 1200 km perimeter).
   */
  traceBetween(from: BoundaryHit, to: BoundaryHit): Position[] | null {
    const [first] = this.traceCandidates(from, to);
    return first ?? null;
  }

  /**
   * Both ways round the ring, shorter path first.
   *
   * `traceBetween` commits to the shorter path, which is right most of the
   * time but not always: the caller also has to satisfy an area sanity check,
   * and when the short way encloses far more ground than the user implied, the
   * whole auto-complete used to be abandoned and the shape closed with a
   * straight line -- without ever trying the other direction, which is often
   * the one the user meant. Handing back both lets the caller choose on area
   * rather than on length alone.
   */
  traceCandidates(from: BoundaryHit, to: BoundaryHit): Position[][] {
    if (from.ringIdx !== to.ringIdx) return [];
    const ring = this.rings[from.ringIdx];
    const segCount = ring.length - 1;
    if (segCount <= 0) return [];

    const forward = walk(ring, segCount, from, to, +1);
    const backward = walk(ring, segCount, from, to, -1);
    const both = [forward, backward].filter((p): p is Position[] => p !== null);
    return both.sort((a, b) => pathLengthDeg(a) - pathLengthDeg(b));
  }
}

/** Build an index from the MultiLineString the app already carries around. */
export function buildBoundaryIndex(
  boundary: Feature<MultiLineString>,
): BoundaryIndex {
  return new BoundaryIndex(boundary.geometry.coordinates as Position[][]);
}

function lerp(a: Position, b: Position, t: number): Position {
  return [a[0] + (b[0] - a[0]) * t, a[1] + (b[1] - a[1]) * t];
}

function pointToSegment(
  px: number,
  py: number,
  a: Position,
  b: Position,
): { distSq: number; snapped: Position; t: number } {
  const ax = a[0] * LNG_SCALE;
  const ay = a[1];
  const bx = b[0] * LNG_SCALE;
  const by = b[1];
  const qx = px * LNG_SCALE;
  const qy = py;
  const vx = bx - ax;
  const vy = by - ay;
  const lenSq = vx * vx + vy * vy;
  let t = lenSq === 0 ? 0 : ((qx - ax) * vx + (qy - ay) * vy) / lenSq;
  t = t < 0 ? 0 : t > 1 ? 1 : t;
  const sx = ax + vx * t;
  const sy = ay + vy * t;
  const dx = qx - sx;
  const dy = qy - sy;
  return { distSq: dx * dx + dy * dy, snapped: lerp(a, b, t), t };
}

/**
 * Walk the ring from `from` to `to` in the given direction, wrapping across the
 * array seam. Returns the inclusive vertex path, or null if the walk exceeds
 * one full lap (which would mean the endpoints were not on this ring).
 */
function walk(
  ring: Position[],
  segCount: number,
  from: BoundaryHit,
  to: BoundaryHit,
  dir: 1 | -1,
): Position[] | null {
  const out: Position[] = [from.snapped];

  // Same segment and the walk direction agrees with the parametric order:
  // the path is just the two snapped points.
  if (from.segIdx === to.segIdx) {
    if (dir === 1 && to.t >= from.t) {
      out.push(to.snapped);
      return out;
    }
    if (dir === -1 && to.t <= from.t) {
      out.push(to.snapped);
      return out;
    }
    // Otherwise we have to go all the way around the ring and come back.
  }

  let seg = from.segIdx;
  for (let steps = 0; steps <= segCount; steps++) {
    // Moving forward we leave segment `seg` via its end vertex (seg+1);
    // moving backward we leave via its start vertex (seg).
    const vertex = dir === 1 ? seg + 1 : seg;
    out.push(ring[vertex]);
    seg = (seg + dir + segCount) % segCount;
    if (seg === to.segIdx) {
      out.push(to.snapped);
      return out;
    }
  }
  return null;
}

/** Proper intersection point of segments p1-p2 and p3-p4, or null. */
function segmentIntersection(
  p1: Position,
  p2: Position,
  p3: Position,
  p4: Position,
): Position | null {
  const x1 = p1[0], y1 = p1[1];
  const x2 = p2[0], y2 = p2[1];
  const x3 = p3[0], y3 = p3[1];
  const x4 = p4[0], y4 = p4[1];
  const denom = (y4 - y3) * (x2 - x1) - (x4 - x3) * (y2 - y1);
  if (denom === 0) return null; // parallel or coincident
  const ua = ((x4 - x3) * (y1 - y3) - (y4 - y3) * (x1 - x3)) / denom;
  const ub = ((x2 - x1) * (y1 - y3) - (y2 - y1) * (x1 - x3)) / denom;
  if (ua < 0 || ua > 1 || ub < 0 || ub > 1) return null;
  return [x1 + ua * (x2 - x1), y1 + ua * (y2 - y1)];
}

/** Path length in scaled degrees -- only used to compare two candidate paths. */
function pathLengthDeg(path: Position[]): number {
  let total = 0;
  for (let i = 1; i < path.length; i++) {
    total += Math.sqrt(
      distSqScaled(path[i - 1][0], path[i - 1][1], path[i][0], path[i][1]),
    );
  }
  return total;
}