/**
 * Zero-dependency marching-squares contour tracer.
 *
 * Runs unchanged in the browser and in Node: the module has no imports at all,
 * never touches the DOM and never writes to the console.
 *
 * Coordinate system
 * -----------------
 * Pixel (x, y) is the unit square whose top-left corner sits at (x, y), so the
 * centre of pixel (x, y) is at (x + 0.5, y + 0.5) and y grows downwards
 * (row 0 is the top row of the input).
 *
 * Contours
 * --------
 * Every returned contour is a *closed* polyline of sub-pixel points; the first
 * point is not repeated at the end. Contours are oriented so that the inside of
 * the shape (the region where the sample is at or above the threshold) lies to
 * the left of the direction of travel. See `contourArea` for the sign this
 * implies.
 */

const EDGE_TOP = 0;
const EDGE_RIGHT = 1;
const EDGE_BOTTOM = 2;
const EDGE_LEFT = 3;

/**
 * Directed segments emitted by each unambiguous marching-squares case.
 *
 * The case index is built from the cell corners as
 * `tl | tr << 1 | br << 2 | bl << 3`. Entries are flat `[from, to, from, to]`
 * pairs of edge ids; the direction keeps the inside region on the left.
 */
const CELL_SEGMENTS = [
  null, //  0 - nothing inside
  [EDGE_LEFT, EDGE_TOP, -1, -1], //  1 - top-left only
  [EDGE_TOP, EDGE_RIGHT, -1, -1], //  2 - top-right only
  [EDGE_LEFT, EDGE_RIGHT, -1, -1], //  3 - top row
  [EDGE_RIGHT, EDGE_BOTTOM, -1, -1], //  4 - bottom-right only
  null, //  5 - saddle (top-left + bottom-right)
  [EDGE_TOP, EDGE_BOTTOM, -1, -1], //  6 - right column
  [EDGE_LEFT, EDGE_BOTTOM, -1, -1], //  7 - everything but bottom-left
  [EDGE_BOTTOM, EDGE_LEFT, -1, -1], //  8 - bottom-left only
  [EDGE_BOTTOM, EDGE_TOP, -1, -1], //  9 - left column
  null, // 10 - saddle (top-right + bottom-left)
  [EDGE_BOTTOM, EDGE_RIGHT, -1, -1], // 11 - everything but bottom-right
  [EDGE_RIGHT, EDGE_LEFT, -1, -1], // 12 - bottom row
  [EDGE_RIGHT, EDGE_TOP, -1, -1], // 13 - everything but top-right
  [EDGE_TOP, EDGE_LEFT, -1, -1], // 14 - everything but top-left
  null, // 15 - everything inside
];

// Case 5 (top-left + bottom-right inside). When the centre of the cell is
// inside, the two inside corners are joined through the middle and the two
// outside corners are separated; otherwise the inside corners are separated.
const SADDLE_5_CONNECTED = [EDGE_LEFT, EDGE_BOTTOM, EDGE_RIGHT, EDGE_TOP];
const SADDLE_5_SPLIT = [EDGE_LEFT, EDGE_TOP, EDGE_RIGHT, EDGE_BOTTOM];

// Case 10 (top-right + bottom-left inside): the mirror image of case 5.
const SADDLE_10_CONNECTED = [EDGE_TOP, EDGE_LEFT, EDGE_BOTTOM, EDGE_RIGHT];
const SADDLE_10_SPLIT = [EDGE_TOP, EDGE_RIGHT, EDGE_BOTTOM, EDGE_LEFT];

const DUPLICATE_EPS = 1e-9;

/**
 * Signed area of a closed polyline, via the shoelace formula. The polyline is
 * implicitly closed (the last point is joined back to the first), so an open
 * ring is fine.
 *
 * Sign convention: this is the plain shoelace sum `Σ (x_i·y_{i+1} − x_{i+1}·y_i) / 2`
 * evaluated in the tracer's y-down pixel coordinates. A ring that runs
 * clockwise *as seen on screen* is therefore positive, and a counter-clockwise
 * one is negative. Because `traceAlphaContours` keeps the inside on the left,
 * the outer boundary of a filled region comes out negative and a hole in it
 * comes out positive.
 *
 * @param {Array<{x: number, y: number}>} points
 * @returns {number} signed area in square pixels
 */
export function contourArea(points) {
  const n = points.length;
  if (!points || n < 3) return 0;
  let sum = 0;
  for (let i = 0; i < n; i++) {
    const a = points[i];
    const b = i + 1 === n ? points[0] : points[i + 1];
    sum += a.x * b.y - b.x * a.y;
  }
  return sum / 2;
}

/**
 * Build an SVG path `d` attribute with one closed subpath per contour.
 *
 * `mapPoint(x, y)` returns the `[X, Y]` pair written to the output, which is
 * what makes it possible to flip the y axis or apply a scale without touching
 * the tracer. Coordinates are rounded to `decimals` places.
 *
 * @param {Array<Array<{x: number, y: number}>>} contours
 * @param {(x: number, y: number) => [number, number]} mapPoint
 * @param {number} [decimals]
 * @returns {string}
 */
export function contoursToPathData(contours, mapPoint, decimals = 2) {
  const places = Math.max(0, Math.floor(decimals));
  const factor = Math.pow(10, places);
  const parts = [];
  for (let c = 0; c < contours.length; c++) {
    const points = contours[c];
    if (!points || points.length < 2) continue;
    for (let i = 0; i < points.length; i++) {
      const mapped = mapPoint(points[i].x, points[i].y);
      // Rounding before formatting keeps `-0.00` out of the output.
      const rx = Math.round(mapped[0] * factor) / factor;
      const ry = Math.round(mapped[1] * factor) / factor;
      parts.push((i === 0 ? 'M' : 'L') + rx.toFixed(places) + ' ' + ry.toFixed(places));
    }
    parts.push('Z');
  }
  return parts.join(' ');
}

// --- closed-ring simplification helpers ------------------------------------

/** Drop points that repeat their predecessor (including across the wrap). */
function removeConsecutiveDuplicates(points) {
  const out = [];
  for (let i = 0; i < points.length; i++) {
    const p = points[i];
    const last = out[out.length - 1];
    if (last && Math.abs(last.x - p.x) <= DUPLICATE_EPS && Math.abs(last.y - p.y) <= DUPLICATE_EPS) {
      continue;
    }
    out.push(p);
  }
  while (out.length > 1) {
    const first = out[0];
    const last = out[out.length - 1];
    if (Math.abs(first.x - last.x) <= DUPLICATE_EPS && Math.abs(first.y - last.y) <= DUPLICATE_EPS) {
      out.pop();
    } else {
      break;
    }
  }
  return out;
}

/** Drop points that sit on the straight segment between their two neighbours. */
function removeCollinear(points) {
  let list = points;
  let changed = true;
  while (changed && list.length > 3) {
    changed = false;
    const n = list.length;
    const out = [];
    for (let i = 0; i < n; i++) {
      const a = list[i === 0 ? n - 1 : i - 1];
      const b = list[i];
      const c = list[i + 1 === n ? 0 : i + 1];
      const abx = b.x - a.x;
      const aby = b.y - a.y;
      const bcx = c.x - b.x;
      const bcy = c.y - b.y;
      const cross = abx * bcy - aby * bcx;
      // |cross| / (|ab| * |bc|) is sin(turn angle); a small value means a
      // straight-through point, which carries no shape information.
      const scale = Math.sqrt((abx * abx + aby * aby) * (bcx * bcx + bcy * bcy));
      const straight = scale <= DUPLICATE_EPS || Math.abs(cross) <= 1e-9 * scale;
      if (straight && abx * bcx + aby * bcy >= 0) {
        changed = true;
      } else {
        out.push(b);
      }
    }
    list = out;
  }
  return list;
}

/**
 * Iterative Douglas–Peucker for an *open* polyline. The two end points are
 * always kept; the recursion uses an explicit stack so long contours cannot
 * overflow the call stack.
 */
function douglasPeucker(points, tolerance) {
  const n = points.length;
  if (n <= 2) return points.slice();
  const keep = new Uint8Array(n);
  keep[0] = 1;
  keep[n - 1] = 1;
  const stack = [0, n - 1];
  while (stack.length > 0) {
    const i1 = stack.pop();
    const i0 = stack.pop();
    if (i1 <= i0 + 1) continue;
    const a = points[i0];
    const b = points[i1];
    const dx = b.x - a.x;
    const dy = b.y - a.y;
    const len = Math.sqrt(dx * dx + dy * dy);
    let maxDistance = -1;
    let maxIndex = -1;
    if (len <= DUPLICATE_EPS) {
      // Degenerate segment: fall back to the distance from the anchor point.
      for (let i = i0 + 1; i < i1; i++) {
        const px = points[i].x - a.x;
        const py = points[i].y - a.y;
        const d = Math.sqrt(px * px + py * py);
        if (d > maxDistance) {
          maxDistance = d;
          maxIndex = i;
        }
      }
    } else {
      for (let i = i0 + 1; i < i1; i++) {
        const p = points[i];
        const d = Math.abs(dy * (p.x - a.x) - dx * (p.y - a.y)) / len;
        if (d > maxDistance) {
          maxDistance = d;
          maxIndex = i;
        }
      }
    }
    if (maxDistance > tolerance && maxIndex > i0) {
      keep[maxIndex] = 1;
      stack.push(i0, maxIndex, maxIndex, i1);
    }
  }
  const out = [];
  for (let i = 0; i < n; i++) {
    if (keep[i]) out.push(points[i]);
  }
  return out;
}

/**
 * Simplify a closed ring, wrap-around segment included.
 *
 * The ring is cut at the point farthest from `points[0]`, which gives two open
 * polylines that together cover every segment of the loop exactly once; each
 * half is then simplified with Douglas–Peucker and the halves are stitched
 * back together (without duplicating the shared anchors).
 */
function simplifyClosedRing(points, tolerance) {
  const n = points.length;
  if (n <= 3) return points.slice();
  let far = 0;
  let farDistance = -1;
  const first = points[0];
  for (let i = 1; i < n; i++) {
    const dx = points[i].x - first.x;
    const dy = points[i].y - first.y;
    const d = dx * dx + dy * dy;
    if (d > farDistance) {
      farDistance = d;
      far = i;
    }
  }
  // A ring whose points all coincide carries no shape; leave it to minArea.
  if (far <= 0 || farDistance <= DUPLICATE_EPS) return points.slice();

  const head = douglasPeucker(points.slice(0, far + 1), tolerance);
  const tail = douglasPeucker(points.slice(far).concat([first]), tolerance);

  // `head` ends and `tail` starts on the same anchor, and both end on
  // `points[0]`; drop the duplicated join so every point appears once.
  return head.slice(0, -1).concat(tail.slice(0, -1));
}

/**
 * Trace the iso-contour of a scalar field at `options.threshold`.
 *
 * @param {ArrayLike<number>} alpha   width*height samples, row-major, row 0 on top
 * @param {number} width
 * @param {number} height
 * @param {object} [options]
 * @param {number} [options.threshold=0.5]        inside when `alpha/255 >= threshold`
 * @param {number} [options.simplifyTolerance=0.35] Douglas-Peucker tolerance, px
 * @param {number} [options.minArea=2]            drop rings smaller than this, px²
 * @returns {Array<Array<{x: number, y: number}>>}
 */
export function traceAlphaContours(alpha, width, height, options = {}) {
  const threshold = options.threshold ?? 0.5;
  const tolerance = options.simplifyTolerance ?? 0.35;
  const minArea = options.minArea ?? 2;

  const w = Math.floor(width);
  const h = Math.floor(height);
  if (!alpha || w < 1 || h < 1 || alpha.length < w * h) return [];
  if (w < 2 && h < 2) return [];

  // Work on a signed field, padded with a 1px outside border: `s >= 0` is
  // inside. The padding guarantees every crossing is strictly interior, so the
  // marching always yields closed rings and never touches the array edges.
  // The border holds the same value a fully transparent pixel maps to, which
  // makes shapes that run off the image close exactly on the image edge.
  const pw = w + 2;
  const ph = h + 2;
  const s = new Float32Array(pw * ph);
  s.fill(-threshold);
  for (let y = 0; y < h; y++) {
    const src = y * w;
    const dst = (y + 1) * pw + 1;
    for (let x = 0; x < w; x++) s[dst + x] = alpha[src + x] / 255 - threshold;
  }
  const sample = (x, y) => s[y * pw + x];

  // Crossing points, keyed by the grid edge they sit on. Both cells sharing an
  // edge call these with the same sample pair in the same order (top→bottom,
  // left→right), so the coordinates come out bit-identical and can be matched
  // by key alone.
  const points = new Map();
  function crossing(key, a, b, x0, y0, dx, dy) {
    let p = points.get(key);
    if (p === undefined) {
      const t = a / (a - b);
      p = { x: x0 + dx * t, y: y0 + dy * t };
      points.set(key, p);
    }
    return p;
  }
  // Horizontal edge of the padded grid at row `py`, spanning columns px..px+1.
  const hKey = (px, py) => `h:${px}:${py}`;
  const hPoint = (px, py) => crossing(hKey(px, py), sample(px, py), sample(px + 1, py), px - 0.5, py - 0.5, 1, 0);
  // Vertical edge of the padded grid at column `px`, spanning rows py..py+1.
  const vKey = (px, py) => `v:${px}:${py}`;
  const vPoint = (px, py) => crossing(vKey(px, py), sample(px, py), sample(px, py + 1), px - 0.5, py - 0.5, 0, 1);

  const edgeKey = [
    (px, py) => hKey(px, py), // EDGE_TOP
    (px, py) => vKey(px + 1, py), // EDGE_RIGHT
    (px, py) => hKey(px, py + 1), // EDGE_BOTTOM
    (px, py) => vKey(px, py), // EDGE_LEFT
  ];
  const edgePoint = [
    (px, py) => hPoint(px, py), // EDGE_TOP
    (px, py) => vPoint(px + 1, py), // EDGE_RIGHT
    (px, py) => hPoint(px, py + 1), // EDGE_BOTTOM
    (px, py) => vPoint(px, py), // EDGE_LEFT
  ];

  // Directed segments: `from` -> `to`, inside region on the left of travel.
  const fromKeys = [];
  const toKeys = [];
  const fromPoints = [];

  for (let py = 0; py < ph - 1; py++) {
    for (let px = 0; px < pw - 1; px++) {
      const tl = sample(px, py);
      const tr = sample(px + 1, py);
      const br = sample(px + 1, py + 1);
      const bl = sample(px, py + 1);
      const inside = (v) => (v >= 0 ? 1 : 0);
      const code = inside(tl) | (inside(tr) << 1) | (inside(br) << 2) | (inside(bl) << 3);
      let segments = CELL_SEGMENTS[code];
      if (code === 5) {
        // With the cell centre inside, the two inside corners join through the
        // middle; otherwise each is cut off on its own.
        segments = (tl + tr + br + bl) / 4 >= 0 ? SADDLE_5_CONNECTED : SADDLE_5_SPLIT;
      } else if (code === 10) {
        segments = (tl + tr + br + bl) / 4 >= 0 ? SADDLE_10_CONNECTED : SADDLE_10_SPLIT;
      }
      if (!segments) continue;
      for (let i = 0; i < segments.length; i += 2) {
        const a = segments[i];
        const b = segments[i + 1];
        if (a < 0 || b < 0) continue; // padding of the single-segment cases
        fromKeys.push(edgeKey[a](px, py));
        fromPoints.push(edgePoint[a](px, py));
        toKeys.push(edgeKey[b](px, py));
        edgePoint[b](px, py); // make sure the shared crossing exists
      }
    }
  }

  // Every crossing has exactly one incoming and one outgoing segment, so the
  // segments can be walked into closed rings without any ambiguity.
  const outgoing = new Map();
  for (let i = 0; i < fromKeys.length; i++) outgoing.set(fromKeys[i], i);

  const segmentCount = fromKeys.length;
  const used = new Uint8Array(segmentCount);
  const contours = [];

  for (let start = 0; start < segmentCount; start++) {
    if (used[start]) continue;
    const ring = [];
    let current = start;
    for (let guard = 0; guard <= segmentCount; guard++) {
      used[current] = 1;
      ring.push(fromPoints[current]);
      const next = outgoing.get(toKeys[current]);
      if (next === undefined || next === start) break;
      if (used[next]) break;
      current = next;
    }
    if (ring.length >= 3) contours.push(ring);
  }

  const result = [];
  for (let i = 0; i < contours.length; i++) {
    let ring = removeConsecutiveDuplicates(contours[i]);
    ring = simplifyClosedRing(ring, tolerance);
    ring = removeCollinear(ring);
    if (ring.length < 3) continue;
    if (Math.abs(contourArea(ring)) < minArea) continue;
    result.push(ring);
  }
  return result;
}