pwaLUPMIS2/src/geom/polygonDivide.js

/**
 * Pure geometry functions for dividing a polygon into N equal-area pieces.
 *
 * No OpenLayers dependency — operates on raw coordinate arrays.
 *
 * The algorithm finds the polygon's longest edge, then places N-1 cutting
 * lines perpendicular to that edge.  Each cutting-line position is found
 * via binary search so that the piece it cuts off has exactly 1/N of the
 * remaining area.  The actual cut is delegated to `splitPolygonByLine()`.
 */

import { splitPolygonByLine } from './polygonSplit.js';

// ── Utility helpers (self-contained) ─────────────────────────────────────────

function dist2(a, b) {
  return (a[0] - b[0]) ** 2 + (a[1] - b[1]) ** 2;
}

/**
 * Signed area of a ring (shoelace formula).
 */
function signedArea(ring) {
  let area = 0;
  for (let i = 0, n = ring.length; i < n - 1; i++) {
    area += ring[i][0] * ring[i + 1][1] - ring[i + 1][0] * ring[i][1];
  }
  return area / 2;
}

/**
 * Absolute polygon area, accounting for holes.
 */
function polygonArea(coords) {
  let area = Math.abs(signedArea(coords[0]));
  for (let i = 1; i < coords.length; i++) {
    area -= Math.abs(signedArea(coords[i]));
  }
  return area;
}

/**
 * Find the longest edge of a ring and return direction vectors.
 *
 * @param {number[][]} ring  Closed ring
 * @returns {{ p0: number[], p1: number[], along: number[], perp: number[] }}
 *   `along` = unit vector along the longest edge,
 *   `perp`  = unit vector perpendicular to `along` (rotated 90° CCW)
 */
function longestEdge(ring) {
  const n = ring.length - 1; // unique vertices
  let bestLen = -1;
  let bestI = 0;

  for (let i = 0; i < n; i++) {
    const d = dist2(ring[i], ring[i + 1]);
    if (d > bestLen) {
      bestLen = d;
      bestI = i;
    }
  }

  const p0 = ring[bestI];
  const p1 = ring[bestI + 1];
  const len = Math.sqrt(bestLen);
  const along = [(p1[0] - p0[0]) / len, (p1[1] - p0[1]) / len];
  // Perpendicular: rotate 90° CCW
  const perp = [-along[1], along[0]];

  return { p0, p1, along, perp };
}

/**
 * Build a cutting line perpendicular to `along` at parameter `t`.
 *
 * The line passes through `origin + t * along` and extends `extent` units
 * in both `perp` directions — long enough to fully cross the polygon.
 */
function makeCuttingLine(origin, along, perp, t, extent) {
  const cx = origin[0] + t * along[0];
  const cy = origin[1] + t * along[1];
  return [
    [cx - extent * perp[0], cy - extent * perp[1]],
    [cx + extent * perp[0], cy + extent * perp[1]],
  ];
}

/**
 * Project the centroid of a polygon's exterior ring onto the `along` axis.
 * Returns the scalar parameter `t` relative to `origin`.
 */
function centroidT(coords, origin, along) {
  const ring = coords[0];
  const n = ring.length - 1;
  let sx = 0, sy = 0;
  for (let i = 0; i < n; i++) {
    sx += ring[i][0];
    sy += ring[i][1];
  }
  const cx = sx / n - origin[0];
  const cy = sy / n - origin[1];
  return cx * along[0] + cy * along[1];
}

// ── Main export ──────────────────────────────────────────────────────────────

/**
 * Divide a polygon into N equal-area pieces by parallel cuts perpendicular
 * to a user-selected edge.
 *
 * @param {number[][][]} polygonCoords  Polygon coordinates [ring, ...holes]
 * @param {number} n  Number of pieces (must be >= 1)
 * @param {number[][]} edgeCoords  The selected edge `[p0, p1]` — cuts will be
 *   perpendicular to this edge direction.
 * @returns {{ pieces: number[][][][], error?: undefined } | { pieces: null, error: string }}
 */
export function dividePolygon(polygonCoords, n, edgeCoords) {
  if (!Number.isInteger(n) || n < 1) {
    return { pieces: null, error: 'Number of divisions must be a positive integer.' };
  }
  if (n === 1) {
    return { pieces: [polygonCoords] };
  }

  const ring = polygonCoords[0];
  const totalArea = polygonArea(polygonCoords);

  if (totalArea < 1e-6) {
    return { pieces: null, error: 'Polygon has no measurable area.' };
  }

  // 1. Determine cutting direction from the selected edge
  let p0, along, perp;
  if (edgeCoords && edgeCoords.length === 2) {
    p0 = edgeCoords[0];
    const dx = edgeCoords[1][0] - edgeCoords[0][0];
    const dy = edgeCoords[1][1] - edgeCoords[0][1];
    const len = Math.sqrt(dx * dx + dy * dy);
    if (len < 1e-10) {
      return { pieces: null, error: 'Selected edge has zero length.' };
    }
    along = [dx / len, dy / len];
    perp = [-along[1], along[0]];
  } else {
    // Fallback: use longest edge
    const edge = longestEdge(ring);
    p0 = edge.p0;
    along = edge.along;
    perp = edge.perp;
  }
  const origin = p0;

  // 2. Project all vertices onto the `along` axis to find extent
  const nVerts = ring.length - 1;
  let tMin = Infinity, tMax = -Infinity;
  for (let i = 0; i < nVerts; i++) {
    const dx = ring[i][0] - origin[0];
    const dy = ring[i][1] - origin[1];
    const t = dx * along[0] + dy * along[1];
    if (t < tMin) tMin = t;
    if (t > tMax) tMax = t;
  }

  // Cutting line extent: enough to cross the polygon in the `perp` direction
  let perpMin = Infinity, perpMax = -Infinity;
  for (let i = 0; i < nVerts; i++) {
    const dx = ring[i][0] - origin[0];
    const dy = ring[i][1] - origin[1];
    const p = dx * perp[0] + dy * perp[1];
    if (p < perpMin) perpMin = p;
    if (p > perpMax) perpMax = p;
  }
  const extent = (perpMax - perpMin) * 1.5; // generous overshoot

  // 3. Iteratively cut pieces
  const pieces = [];
  let remaining = polygonCoords;
  let remainingCount = n;

  for (let i = 0; i < n - 1; i++) {
    const remainingArea = polygonArea(remaining);
    const targetArea = remainingArea / remainingCount;

    // Re-project the remaining polygon to get its current t-range
    const remRing = remaining[0];
    const remN = remRing.length - 1;
    let rMin = Infinity, rMax = -Infinity;
    for (let j = 0; j < remN; j++) {
      const dx = remRing[j][0] - origin[0];
      const dy = remRing[j][1] - origin[1];
      const t = dx * along[0] + dy * along[1];
      if (t < rMin) rMin = t;
      if (t > rMax) rMax = t;
    }

    // Binary search for the cutting position
    let lo = rMin;
    let hi = rMax;
    let bestT = (lo + hi) / 2;
    let bestPiece = null;
    let bestRemaining = null;
    let bestError = Infinity;

    for (let iter = 0; iter < 40; iter++) {
      const mid = (lo + hi) / 2;
      const line = makeCuttingLine(origin, along, perp, mid, extent);
      const result = splitPolygonByLine(remaining, line);

      if (!result) {
        // Cutting line didn't produce a valid split — nudge and retry
        // Try slightly shifted positions
        const nudge = (hi - lo) * 0.01;
        const lineA = makeCuttingLine(origin, along, perp, mid + nudge, extent);
        const resultA = splitPolygonByLine(remaining, lineA);
        if (resultA) {
          const [halfA, halfB] = resultA;
          const tA = centroidT(halfA, origin, along);
          const tB = centroidT(halfB, origin, along);
          const nearPiece = tA < tB ? halfA : halfB;
          const farPiece = tA < tB ? halfB : halfA;
          const nearArea = polygonArea(nearPiece);
          const err = Math.abs(nearArea - targetArea);
          if (err < bestError) {
            bestError = err;
            bestT = mid + nudge;
            bestPiece = nearPiece;
            bestRemaining = farPiece;
          }
        }
        // Try the other direction
        const lineB = makeCuttingLine(origin, along, perp, mid - nudge, extent);
        const resultB = splitPolygonByLine(remaining, lineB);
        if (resultB) {
          const [halfA, halfB] = resultB;
          const tA = centroidT(halfA, origin, along);
          const tB = centroidT(halfB, origin, along);
          const nearPiece = tA < tB ? halfA : halfB;
          const farPiece = tA < tB ? halfB : halfA;
          const nearArea = polygonArea(nearPiece);
          const err = Math.abs(nearArea - targetArea);
          if (err < bestError) {
            bestError = err;
            bestT = mid - nudge;
            bestPiece = nearPiece;
            bestRemaining = farPiece;
          }
        }
        // Bisect anyway to keep converging
        lo = mid;
        continue;
      }

      const [halfA, halfB] = result;
      const tA = centroidT(halfA, origin, along);
      const tB = centroidT(halfB, origin, along);
      const nearPiece = tA < tB ? halfA : halfB;
      const farPiece = tA < tB ? halfB : halfA;
      const nearArea = polygonArea(nearPiece);

      const err = Math.abs(nearArea - targetArea);
      if (err < bestError) {
        bestError = err;
        bestT = mid;
        bestPiece = nearPiece;
        bestRemaining = farPiece;
      }

      // Converged?
      if (err / remainingArea < 0.001) break;

      // Adjust search range
      if (nearArea < targetArea) {
        lo = mid; // need to cut farther out
      } else {
        hi = mid; // need to cut closer
      }
    }

    if (!bestPiece || !bestRemaining) {
      return {
        pieces: null,
        error: `Could not find a valid cut for piece ${i + 1} of ${n}. The polygon shape may be too irregular for equal division.`,
      };
    }

    pieces.push(bestPiece);
    remaining = bestRemaining;
    remainingCount--;
  }

  // The last remaining piece is the Nth piece
  pieces.push(remaining);

  return { pieces };
}