Arcs, circles, smooth curves, divide and offset

/** A regular polygon (or a circle with `sides` segments) centred at c, radius r, normal n, first corner along u. */
export function ngonPts(c, n, u, r, sides) {
  const v = cross(norm(n), u);
  const out = [];
  for (let i = 0; i < sides; i++) {
    const a = (i / sides) * Math.PI * 2;
    out.push(add(c, add(mul(u, r * Math.cos(a)), mul(v, r * Math.sin(a)))));
  }
  return out;
}

/** Points of an arc: centre c, unit normal n, first direction u (from c to the start), radius r, sweep (radians), segs pieces. */
export function arcPts(c, n, u, r, sweep, segs) {
  const v = cross(norm(n), u);
  const out = [];
  for (let i = 0; i <= segs; i++) { const a = (sweep * i) / segs; out.push(add(c, add(mul(u, r * Math.cos(a)), mul(v, r * Math.sin(a))))); }
  return out;
}

/**
 * What a curve is: {id, kind, edges, verts (in order along it), closed, length, segs, round (all on one circle), center,
 * normal, radius, flat (no face it bounds leaves its plane: its radius can change)} — or null when it has no edges.
 */
export function curveInfo(def, cid) {
  const eids = curveEdges(def, cid);
  if (!eids.length) return null;
  const adj = new Map();
  let length = 0;
  for (const id of eids) {
    const e = def.edges.get(id);
    for (const [x, y] of [[e.a, e.b], [e.b, e.a]]) { if (!adj.has(x)) adj.set(x, []); adj.get(x).push(y); }
    length += dist(def.verts.get(e.a), def.verts.get(e.b));
  }
  const ends = [...adj].filter(([, l]) => l.length === 1).map(([v]) => v);
  const closed = !ends.length && [...adj.values()].every((l) => l.length === 2);
  const chain = [];
  let cur = ends[0] ?? [...adj.keys()][0], prev = null;
  while (cur !== undefined && chain.length <= adj.size) {
    chain.push(cur);
    const nx = (adj.get(cur) || []).find((x) => x !== prev && !chain.includes(x));
    prev = cur; cur = nx;
  }
  const P = chain.map((v) => def.verts.get(v));
  const out = { id: cid, kind: cid.split(':')[0], edges: eids, verts: chain, closed, length, segs: eids.length, round: false, center: null, normal: null, radius: 0, flat: false, sweep: 0 };
  if (P.length >= 3 && chain.length === adj.size) {
    const a = P[0], b = P[Math.floor(P.length / 3)], c = P[Math.floor((2 * P.length) / 3)];
    const ab = sub(b, a), ac = sub(c, a), n = cross(ab, ac), nn = dot(n, n);
    if (nn > 1e-24) {
      const cc = add(a, mul(add(mul(cross(n, ab), dot(ac, ac)), mul(cross(ac, n), dot(ab, ab))), 1 / (2 * nn)));
      const r = dist(cc, a), nrm = norm(n), tol = Math.max(EPS, r * 1e-7);
      if (P.every((p) => Math.abs(dist(p, cc) - r) <= tol && Math.abs(dot(sub(p, cc), nrm)) <= tol)) {
        Object.assign(out, { round: true, center: cc, normal: nrm, radius: r });
        const ef = edgeFaces(def);
        out.flat = eids.every((id) => { const e = def.edges.get(id); return (ef.get(key(e.a, e.b)) || []).every((fid) => Math.abs(dot(faceNormal(def, def.faces.get(fid)), nrm)) > 1 - 1e-9); });
        let sw = 0;
        for (let i = 0; i + 1 < P.length + (closed ? 1 : 0); i++) {
          const u1 = norm(sub(P[i], cc)), u2 = norm(sub(P[(i + 1) % P.length], cc));
          sw += Math.atan2(dot(cross(u1, u2), nrm), dot(u1, u2));
        }
        out.sweep = Math.abs(sw);
      }
    }
  }
  return out;
}

/** Give a round, flat curve a new radius about its centre (its points move; faces in its plane follow). */
export function setCurveRadius(def, cid, r) {
  const ci = curveInfo(def, cid);
  if (!ci || !ci.round) throw new Error(t('m3d.model.only_arcs_circles_and_polygons_have'));
  if (!ci.flat) throw new Error(t('m3d.model.this_curve_bounds_a_face_that'));
  if (!(r > EPS)) throw new Error(t('m3d.model.type_a_radius_above_zero'));
  const k = r / ci.radius, c = ci.center;
  for (const v of ci.verts) def.verts.set(v, add(c, mul(sub(def.verts.get(v), c), k)));
  return curveInfo(def, cid);
}

/**
 * Where n equal parts of an edge or curve would meet (world-free: the definition's coordinates), for previews and for
 * the divide itself. `what` = {edge: id} or {curve: id, edge?: the edge picked on it}. → {points (the n − 1 new places, or for an arc or circle that is
 * re-drawn: every corner of the new curve), resegment (true: an arc or circle whose segments are replaced), length} |
 * null when it is gone.
 */
export function dividePlan(def, what, n) {
  if (!(Number.isInteger(n) && n >= 2 && n <= DIVIDE_MAX)) throw new Error(t('m3d.model.divide_count', { max: DIVIDE_MAX }));
  if (what.curve === undefined) {
    const e = def.edges.get(what.edge);
    if (!e) return null;
    const a = def.verts.get(e.a), b = def.verts.get(e.b);
    return { points: Array.from({ length: n - 1 }, (_, i) => add(a, mul(sub(b, a), (i + 1) / n))), resegment: false, length: dist(a, b) };
  }
  let ci = curveInfo(def, what.curve);
  if (!ci) return null;
  // one curve id can name several separate runs (push/pull gives a pulled circle's top the bottom's id): the run with
  // the edge picked (else the first) is divided, along its own length
  const run = curveRun(def, what.curve, what.edge ?? ci.edges[0]);
  if (run.edges.length !== ci.edges.length) ci = { ...ci, ...run, round: false };
  // an arc or circle whose own corners nothing else uses, lying flat in its faces' planes: re-drawn with n segments
  if (ci.round && ci.flat && (ci.kind === 'arc' || ci.kind === 'circle') && ci.verts.length === ci.segs + (ci.closed ? 0 : 1)) {
    const ve = vertEdges(def), inner = ci.closed ? ci.verts : ci.verts.slice(1, -1);
    if (inner.every((v) => (ve.get(v) || []).length === 2)) {
      const p0 = def.verts.get(ci.verts[0]), p1 = def.verts.get(ci.verts[1]), u = norm(sub(p0, ci.center));
      // which way round the chain runs about the normal
      const sgn = dot(cross(sub(p0, ci.center), sub(p1, ci.center)), ci.normal) >= 0 ? 1 : -1;
      const pts = arcPts(ci.center, mul(ci.normal, sgn), u, ci.radius, ci.closed ? 2 * Math.PI : ci.sweep, n);
      if (ci.closed) pts.pop();
      else pts[pts.length - 1] = def.verts.get(ci.verts[ci.verts.length - 1]).slice();   // the far end stays exactly where it was
      pts[0] = p0.slice();
      return { points: pts, resegment: true, length: ci.radius * (ci.closed ? 2 * Math.PI : ci.sweep), curve: ci };
    }
  }
  // any other curve: equal lengths along its segments
  const chain = ci.closed ? [...ci.verts, ci.verts[0]] : ci.verts;
  const P = chain.map((v) => def.verts.get(v));
  const L = ci.length, out = [];
  let s0 = 0, i = 0;
  for (let k = 1; k < n; k++) {
    const want = (k * L) / n;
    while (i < P.length - 2 && s0 + dist(P[i], P[i + 1]) < want) { s0 += dist(P[i], P[i + 1]); i++; }
    const seg = dist(P[i], P[i + 1]);
    out.push(add(P[i], mul(sub(P[i + 1], P[i]), seg ? Math.min(1, (want - s0) / seg) : 0)));
  }
  return { points: out, resegment: false, length: L, curve: ci };
}

/**
 * Divide an edge (`{edge: id}`) or a curve (`{curve: id}`) of `def` into n equal parts (0.7 MC-DIVIDE). A straight edge
 * gets n − 1 corners; an arc or circle that only its own faces touch is re-drawn with n segments on the true arc (as a
 * curve: it still selects as one and keeps its radius); any other curve gets corners at equal lengths along it. Faces
 * keep their place. → {verts (the corners added), edges (every piece, in order)}.
 */
export function divideEdges(def, what, n) {
  const plan = dividePlan(def, what, n);
  if (!plan) throw new Error(t('m3d.model.divide_gone'));
  if (plan.resegment) return resegmentCurve(def, plan.curve, plan.points);
  let chain;
  if (what.curve === undefined) { const e = def.edges.get(what.edge); chain = [e.a, e.b]; }
  else chain = plan.curve.closed ? [...plan.curve.verts, plan.curve.verts[0]] : plan.curve.verts.slice();
  // walk the chain, splitting the segment each new corner lands on (a corner already there is kept, not doubled)
  const made = [], P = chain.map((v) => def.verts.get(v));
  let seg = 0, from = chain[0];
  const cum = [0];
  for (let k = 0; k + 1 < P.length; k++) cum.push(cum[k] + dist(P[k], P[k + 1]));
  for (let k = 1; k < n; k++) {
    const want = (k * cum[cum.length - 1]) / n;
    while (seg < P.length - 2 && cum[seg + 1] < want - EPS) { seg++; from = chain[seg]; }
    if (Math.abs(cum[seg + 1] - want) <= EPS) continue;          // lands on an existing corner
    if (Math.abs(cum[seg] - want) <= EPS) continue;
    const span = cum[seg + 1] - cum[seg], p = add(P[seg], mul(sub(P[seg + 1], P[seg]), (want - cum[seg]) / span));
    const eid = edgeBetween(def, from, chain[seg + 1]);
    if (eid === null) continue;
    const v = vAdd(def, p);
    splitEdge(def, eid, v);
    made.push(v);
    from = v;
  }
  const edges = [];
  const all = [];
  for (let k = 0; k + 1 < chain.length; k++) {
    all.push(chain[k]);
    const a = def.verts.get(chain[k]), b = def.verts.get(chain[k + 1]);
    all.push(...made.filter((v) => { const p = def.verts.get(v); return onSegment(p, a, b); }).sort((x, y) => dist(def.verts.get(x), a) - dist(def.verts.get(y), a)));
  }
  all.push(chain[chain.length - 1]);
  for (let k = 0; k + 1 < all.length; k++) { const id = edgeBetween(def, all[k], all[k + 1]); if (id !== null) edges.push(id); }
  return { verts: made, edges };
}

/** Points along a centripetal Catmull–Rom spline through P (3D), segs per span; closed wraps round. */
export function smoothPts(P, closed = false, segs = SPLINE_SEGS) {
  const n = P.length;
  if (n < 3) return P.map((p) => p.slice());
  const at = (i) => (closed ? P[(i + n) % n] : i < 0 ? sub(mul(P[0], 2), P[1]) : i >= n ? sub(mul(P[n - 1], 2), P[n - 2]) : P[i]);
  const out = [];
  const spans = closed ? n : n - 1;
  const lerp = (a, b, ta, tb, t) => add(mul(a, (tb - t) / (tb - ta)), mul(b, (t - ta) / (tb - ta)));
  for (let i = 0; i < spans; i++) {
    const p0 = at(i - 1), p1 = at(i), p2 = at(i + 1), p3 = at(i + 2);
    const tj = (a, b) => Math.max(1e-9, Math.sqrt(dist(a, b)));
    const t1 = tj(p0, p1), t2 = t1 + tj(p1, p2), t3 = t2 + tj(p2, p3);
    for (let k = 0; k < segs; k++) {
      const t = t1 + ((t2 - t1) * k) / segs;
      const A1 = lerp(p0, p1, 0, t1, t), A2 = lerp(p1, p2, t1, t2, t), A3 = lerp(p2, p3, t2, t3, t);
      const B1 = lerp(A1, A2, 0, t2, t), B2 = lerp(A2, A3, t1, t3, t);
      out.push(lerp(B1, B2, t1, t2, t));
    }
  }
  if (!closed) out.push(P[n - 1].slice());
  return out;
}

/**
 * Offset a 2D polyline: every segment's line moves `d` to its left, corners where neighbouring lines meet (closed:
 * the loop wraps). Returns the new points, or null when a segment would flip (the offset is larger than the shape allows).
 */
export function offset2(P, d, closed) {
  const n = P.length;
  if (n < 2) return null;
  const segs = [];
  for (let i = 0; i < (closed ? n : n - 1); i++) {
    const a = P[i], b = P[(i + 1) % n], dx = b[0] - a[0], dy = b[1] - a[1], L = Math.hypot(dx, dy);
    if (!(L > 1e-12)) return null;
    const nx = -dy / L, ny = dx / L;
    segs.push({ a: [a[0] + nx * d, a[1] + ny * d], u: [dx / L, dy / L] });
  }
  const meet = (s, t) => {
    const den = s.u[0] * t.u[1] - s.u[1] * t.u[0];
    if (Math.abs(den) < 1e-9) return [t.a[0], t.a[1]];   // in line: the shifted point itself
    const k = ((t.a[0] - s.a[0]) * t.u[1] - (t.a[1] - s.a[1]) * t.u[0]) / den;
    return [s.a[0] + s.u[0] * k, s.a[1] + s.u[1] * k];
  };
  const Q = [];
  for (let i = 0; i < n; i++) {
    if (!closed && i === 0) { Q.push(segs[0].a); continue; }
    if (!closed && i === n - 1) { const s = segs[n - 2]; const L = Math.hypot(P[n - 1][0] - P[n - 2][0], P[n - 1][1] - P[n - 2][1]); Q.push([s.a[0] + s.u[0] * L, s.a[1] + s.u[1] * L]); continue; }
    Q.push(meet(segs[(i - 1 + segs.length) % segs.length], segs[i]));
  }
  for (let i = 0; i < segs.length; i++) {
    const a = Q[i], b = Q[(i + 1) % n];
    if ((b[0] - a[0]) * segs[i].u[0] + (b[1] - a[1]) * segs[i].u[1] <= 1e-9) return null;
  }
  return Q;
}

function segsCross2(A, B, C, D) {
  const o = (p, q, r) => (q[0] - p[0]) * (r[1] - p[1]) - (q[1] - p[1]) * (r[0] - p[0]);
  const d1 = o(C, D, A), d2 = o(C, D, B), d3 = o(A, B, C), d4 = o(A, B, D);
  return ((d1 > 1e-12 && d2 < -1e-12) || (d1 < -1e-12 && d2 > 1e-12)) && ((d3 > 1e-12 && d4 < -1e-12) || (d3 < -1e-12 && d4 > 1e-12));
}

/**
 * The loops of face fid offset by d (positive = into the face, negative = out of it), in the def's coordinates:
 * [[p…]…] in the face's loop order, or throws with the reason.
 */
export function offsetFacePts(def, fid, d) {
  const f = def.faces.get(fid);
  if (!f) throw new Error(t('m3d.model.there_is_no_such_face'));
  const n = faceNormal(def, f), fr = frame(n), o = def.verts.get(f.loops[0][0]);
  const h = dot(o, n);
  const loops2 = f.loops.map((l) => l.map((v) => to2(fr, def.verts.get(v))));
  const out2 = [];
  for (const L of loops2) {
    const Q = offset2(L, d, true);
    if (!Q) throw new Error(t('m3d.model.that_is_further_than_this_face'));
    out2.push(Q);
  }
  if (loopsCross(out2)) throw new Error(t('m3d.model.that_is_further_than_this_face_2'));
  return out2.map((Q) => Q.map(([x, y]) => add(add(mul(fr.u, x), mul(fr.v, y)), mul(n, h))));
}

/**
 * Offset a chain of connected edges (all in one plane with normal n) by d to the left of its direction (seen from n's tip).
 * Returns the new points (def coordinates), or throws with the reason.
 */
export function offsetChainPts(def, eids, d, n) {
  const adj = new Map();
  for (const id of eids) {
    const e = def.edges.get(id);
    if (!e) continue;
    for (const [x, y] of [[e.a, e.b], [e.b, e.a]]) { if (!adj.has(x)) adj.set(x, []); adj.get(x).push(y); }
  }
  if ([...adj.values()].some((l) => l.length > 2)) throw new Error(t('m3d.model.offset_works_on_one_chain_of'));
  const ends = [...adj].filter(([, l]) => l.length === 1).map(([v]) => v);
  const closed = !ends.length;
  const chain = [];
  let cur = ends[0] ?? [...adj.keys()][0], prev = null;
  while (cur !== undefined && !chain.includes(cur)) { chain.push(cur); const nx = adj.get(cur).find((x) => x !== prev && !chain.includes(x)); prev = cur; cur = nx; }
  if (chain.length !== adj.size) throw new Error(t('m3d.model.offset_works_on_one_connected_chain'));
  const P3 = chain.map((v) => def.verts.get(v));
  const fr = frame(norm(n)), h = dot(P3[0], fr.n);
  if (P3.some((p) => Math.abs(dot(p, fr.n) - h) > Math.max(EPS, 1e-7 * dist(p, P3[0])))) throw new Error(t('m3d.model.these_edges_are_not_all_in'));
  const Q = offset2(P3.map((p) => to2(fr, p)), d, closed);
  if (!Q || (closed && loopsCross([Q]))) throw new Error(t('m3d.model.that_is_further_than_these_edges'));
  return { pts: Q.map(([x, y]) => add(add(mul(fr.u, x), mul(fr.v, y)), mul(fr.n, h))), closed, chain };
}