Drawing an edge: welding, splitting and re-facing

/** The vertex at p: an existing one, or a new one (splitting the edge it lands on). */
export function vertexAt(def, p) {
  const have = findVert(def, p);
  if (have !== null) return have;
  const id = vAdd(def, p);
  for (const [eid, e] of [...def.edges]) {
    const A = def.verts.get(e.a), B = def.verts.get(e.b), u = sub(B, A), L2 = dot(u, u);
    const s = dot(sub(p, A), u) / L2;
    if (s <= 0 || s >= 1) continue;
    if (dist(add(A, mul(u, s)), p) <= EPS) { splitEdge(def, eid, id); break; }
  }
  return id;
}

const A = def.verts.get(a), B = def.verts.get(b), u = sub(B, A), L2 = dot(u, u);
const stops = [{ s: 0, v: a }, { s: 1, v: b }];
for (const [vid, P] of def.verts) {
  if (vid === a || vid === b) continue;
  const s = dot(sub(P, A), u) / L2;
  if (s <= 0 || s >= 1) continue;
  if (dist(add(A, mul(u, s)), P) <= EPS) stops.push({ s, v: vid });
}
for (const [eid, e] of [...def.edges]) {
  if (e.a === a || e.a === b || e.b === a || e.b === b) continue;
  if (stops.some((st) => st.v === e.a || st.v === e.b)) continue;
  const C = def.verts.get(e.a), D = def.verts.get(e.b), w = sub(D, C);
  const r = segSeg(A, u, C, w);
  if (!r) continue;
  const v = vAdd(def, r.p);
  splitEdge(def, eid, v);
  stops.push({ s: r.s, v });
}
stops.sort((x, y) => x.s - y.s);

/** Proper crossing of segments A + s·u and C + t·w (both strictly inside), coplanar within EPS. */
function segSeg(A, u, C, w) {
  const n = cross(u, w), nn = dot(n, n);
  if (nn < 1e-18 * dot(u, u) * dot(w, w)) return null;
  const AC = sub(C, A);
  const s = dot(cross(AC, w), n) / nn, t = dot(cross(AC, u), n) / nn;
  const tolS = EPS / Math.sqrt(dot(u, u)), tolT = EPS / Math.sqrt(dot(w, w));
  if (s <= tolS || s >= 1 - tolS || t <= tolT || t >= 1 - tolT) return null;
  const p1 = add(A, mul(u, s)), p2 = add(C, mul(w, t));
  if (dist(p1, p2) > EPS) return null;
  return { s, t, p: p1 };
}

/** Planes an edge could close a face on: planes of faces it lies on, and of each neighbour edge. */
function candidatePlanes(def, eids) {
  const planes = [];
  const put = (pl) => { if (!planes.some((q) => samePlane(q, pl))) planes.push(pl); };
  const ve = vertEdges(def);
  for (const eid of eids) {
    const e = def.edges.get(eid);
    if (!e) continue;
    const A = def.verts.get(e.a), B = def.verts.get(e.b), u = sub(B, A);
    for (const f of def.faces.values()) {
      const pl = facePlane(def, f);
      if (len(pl.n) && onPlane(pl, A) && onPlane(pl, B)) put(planeFrom(pl.n, A));
    }
    for (const end of [e.a, e.b]) for (const id2 of ve.get(end) || []) {
      if (id2 === eid) continue;
      const e2 = def.edges.get(id2), C = def.verts.get(e2.a === end ? e2.b : e2.a);
      const n = cross(u, sub(C, A));
      if (len(n) > 1e-9 * len(u) * dist(C, A)) put(planeFrom(n, A));
    }
  }
  return planes;
}

/**
 * Rebuild the faces of one plane from its edge arrangement (see the header). `newEdges` may close
 * new faces; `hint` (a direction) picks the side a new face faces.
 */
export function reface(def, plane, { newEdges = new Set(), hint = null } = {}) {
  const fr = frame(plane.n);
  const onP = new Set();
  for (const [id, p] of def.verts) if (onPlane(plane, p)) onP.add(id);
  const old = [];
  for (const [fid, f] of def.faces) {
    if (!f.loops.every((l) => l.every((v) => onP.has(v)))) continue;
    const n = faceNormal(def, f);
    if (!len(n) || Math.abs(dot(n, plane.n)) < 1 - 1e-6) continue;
    old.push({ fid, f, s: Math.sign(dot(n, plane.n)), poly: f.loops.map((l) => l.map((v) => to2(fr, def.verts.get(v)))) });
  }
  const E = [];
  for (const [eid, e] of def.edges) if (onP.has(e.a) && onP.has(e.b)) E.push(eid);
  if (!old.length && !E.some((id) => newEdges.has(id))) return false;
  const P2 = new Map();
  for (const v of onP) P2.set(v, to2(fr, def.verts.get(v)));

  const live = new Set(E);
  let cycles = [];
  for (let guard = 0; guard < 200; guard++) {
    pruneDangling(def, live);
    cycles = walkCycles(def, live, P2);
    const bridges = new Set();
    for (const c of cycles) { const seen = new Set(); for (const e of c.edges) { if (seen.has(e)) bridges.add(e); seen.add(e); } }
    if (!bridges.size) break;
    for (const b of bridges) live.delete(b);
  }
  const pos = cycles.filter((c) => c.area > EPS2), negs = cycles.filter((c) => c.area < -EPS2);
  const holesOf = new Map(pos.map((c) => [c, []]));
  for (const h of negs) {
    let best = null;
    for (const c of pos) {
      const vs = new Set(c.verts);
      const probe = h.verts.find((v) => !vs.has(v));
      if (probe === undefined || !inPoly(P2.get(probe), c.poly)) continue;
      if (!best || c.area < best.area) best = c;
    }
    if (best) holesOf.get(best).push(h);
  }
  const made = [];
  const reused = new Set();
  for (const c of pos) {
    const holes = holesOf.get(c);
    const pt = samplePoint([c.poly, ...holes.map((h) => h.poly)], old.map((o) => o.poly));
    if (!pt) continue;
    let inherit = null, sgn = 0, mixed = false;
    for (const o of old) {
      if (!inRegion(pt, o.poly)) continue;
      if (!inherit) inherit = o;
      if (sgn && sgn !== o.s) mixed = true;
      sgn = o.s;
    }
    if (mixed) continue;
    if (!inherit) {
      if (!c.edges.some((e) => newEdges.has(e))) continue;
      sgn = hint && dot(hint, plane.n) < 0 ? -1 : 1;
    }
    let loops = [c.verts.slice(), ...holes.map((h) => h.verts.slice())];
    if (sgn < 0) loops = loops.map((l) => l.reverse());
    const fid = inherit && !reused.has(inherit.fid) ? inherit.fid : null;
    if (fid !== null) reused.add(fid);
    made.push({ fid, loops, color: inherit?.f.color, hidden: inherit?.f.hidden });
  }
  for (const o of old) def.faces.delete(o.fid);
  for (const m of made) fAdd(def, m.loops, m.color, m.fid ?? undefined, m.hidden);
  return true;
}

/** Faces of a planar graph: each cycle keeps its region on the left (CCW positive area). */
function walkCycles(def, live, P2) {
  const out = new Map();   // vertex → [{to, e, ang}] sorted by angle
  for (const id of live) {
    const e = def.edges.get(id);
    for (const [u, v] of [[e.a, e.b], [e.b, e.a]]) {
      const p = P2.get(u), q = P2.get(v);
      if (!out.has(u)) out.set(u, []);
      out.get(u).push({ from: u, to: v, e: id, ang: Math.atan2(q[1] - p[1], q[0] - p[0]) });
    }
  }
  for (const l of out.values()) l.sort((x, y) => x.ang - y.ang);
  const seen = new Set();
  const hk = (h) => h.from + '>' + h.to;
  const cycles = [];
  for (const l of out.values()) for (const h0 of l) {
    if (seen.has(hk(h0))) continue;
    const verts = [], edges = [];
    let h = h0, guard = 0;
    while (!seen.has(hk(h)) && guard++ < 100000) {
      seen.add(hk(h));
      verts.push(h.from);
      edges.push(h.e);
      const list = out.get(h.to);
      const i = list.findIndex((x) => x.to === h.from);
      h = list[(i - 1 + list.length) % list.length];
    }
    const poly = verts.map((v) => P2.get(v));
    cycles.push({ verts, edges, poly, area: area2(poly) });
  }
  return cycles;
}

function inPoly(pt, poly) {
  let inside = false;
  for (let i = 0, j = poly.length - 1; i < poly.length; j = i++) {
    const a = poly[i], b = poly[j];
    if ((a[1] > pt[1]) !== (b[1] > pt[1]) && pt[0] < (b[0] - a[0]) * (pt[1] - a[1]) / (b[1] - a[1]) + a[0]) inside = !inside;
  }
  return inside;
}

/**
 * A point well inside a region (outer + holes): the centroid of its largest triangle. With `avoid` (polygons: lists of
 * loops of 2D points) the point is also kept off their edges: reface asks which old faces hold the point, and a point ON
 * an old face's edge (the edge being merged away, say) is in both faces or in neither depending on the last bit of
 * rounding, so a region could inherit nothing and be dropped (kernel-print, 2026-10-23: a 2 cm part lost its sides,
 * because both centroids of a square cut along the other diagonal lie on the old one). Several points of the largest
 * triangles are tried and the one farthest from every avoided edge wins (the first clear by 10 µm, or the best).
 */
function samplePoint(loops, avoid = null) {
  const flat = [], holes = [], pts = [];
  loops.forEach((l, i) => { if (i) holes.push(pts.length); for (const p of l) { flat.push(p[0], p[1]); pts.push(p); } });
  const t = earcut(flat, holes);
  const tris = [];
  for (let i = 0; i < t.length; i += 3) {
    const a = pts[t[i]], b = pts[t[i + 1]], c = pts[t[i + 2]];
    tris.push({ a, b, c, ar: Math.abs((b[0] - a[0]) * (c[1] - a[1]) - (c[0] - a[0]) * (b[1] - a[1])) });
  }
  if (!tris.length) return null;
  tris.sort((x, y) => y.ar - x.ar);
  if (!(tris[0].ar > 1e-14)) return null;
  const at = (q, wa, wb, wc) => [q.a[0] * wa + q.b[0] * wb + q.c[0] * wc, q.a[1] * wa + q.b[1] * wb + q.c[1] * wc];
  const first = at(tris[0], 1 / 3, 1 / 3, 1 / 3);
  if (!avoid || !avoid.length) return first;
  const clear = (p) => {
    let m = Infinity;
    for (const poly of avoid) for (const l of poly) for (let i = 0; i < l.length; i++) { const d = segDist2(p, l[i], l[(i + 1) % l.length]); if (d < m) m = d; }
    return m;
  };
  const W = [[1 / 3, 1 / 3, 1 / 3], [0.5, 0.25, 0.25], [0.25, 0.5, 0.25], [0.25, 0.25, 0.5]];
  let best = first, bd = -1;
  for (const q of tris.slice(0, 4)) {
    if (!(q.ar > 1e-14)) break;
    for (const w of W) {
      const p = at(q, w[0], w[1], w[2]), d = clear(p);
      if (d > bd) { bd = d; best = p; }
      if (d > 10 * EPS) return best;
    }
  }
  return best;
}

export function cleanup(def) {
  pruneSplines(def);
  const ids = [...def.verts.keys()].sort((a, b) => def.verts.get(a)[0] - def.verts.get(b)[0]);
  const to = new Map();
  for (let i = 0; i < ids.length; i++) {
    if (to.has(ids[i])) continue;
    const p = def.verts.get(ids[i]);
    for (let j = i + 1; j < ids.length; j++) {
      const q = def.verts.get(ids[j]);
      if (q[0] - p[0] > EPS) break;
      if (!to.has(ids[j]) && near(p, q)) to.set(ids[j], ids[i]);
    }
  }
  let changed = to.size > 0;
  if (to.size) {
    const m = (v) => to.get(v) ?? v;
    const edges = [...def.edges];
    def.edges.clear(); def._ek = null;
    for (const [id, e] of edges) {
      const a = m(e.a), b = m(e.b);
      if (a === b || edgeBetween(def, a, b) !== null) continue;
      def.edges.set(id, { ...e, a, b });
      ek(def).set(key(a, b), id);
    }
    for (const f of def.faces.values()) f.loops = f.loops.map((l) => l.map(m));
    for (const v of to.keys()) def.verts.delete(v);
  }
  for (const [fid, f] of [...def.faces]) {
    const loops = [];
    f.loops.forEach((l, i) => {
      const c = l.filter((v, k) => v !== l[(k + 1) % l.length]);
      if (c.length >= 3 && polyArea3(c.map((v) => def.verts.get(v))) > EPS2) loops.push(c);
      else if (i === 0) loops.push(null);
    });
    if (!loops.length || loops[0] === null) { def.faces.delete(fid); changed = true; continue; }
    if (loops.length !== f.loops.length || loops.some((l, i) => l.length !== f.loops[i].length)) { f.loops = loops; changed = true; }
    for (const l of f.loops) for (let i = 0; i < l.length; i++) eAdd(def, l[i], l[(i + 1) % l.length]);
  }
  removeLooseVerts(def);
  return changed;
}

/** A vertex between two collinear edges (and nothing else) disappears; the edges become one. */
export function healCollinear(def, vids) {
  for (const v of vids) {
    if (!def.verts.has(v)) continue;
    const ve = vertEdges(def).get(v) || [];
    if (ve.length !== 2) continue;
    const e1 = def.edges.get(ve[0]), e2 = def.edges.get(ve[1]);
    const a = e1.a === v ? e1.b : e1.a, b = e2.a === v ? e2.b : e2.a;
    if (a === b) continue;
    const V = def.verts.get(v), A = def.verts.get(a), B = def.verts.get(b);
    const u1 = norm(sub(A, V)), u2 = norm(sub(B, V));
    if (dot(u1, u2) > -1 + 1e-9) continue;
    if (edgeBetween(def, a, b) !== null) continue;
    const ef = edgeFaces(def);
    const f1 = (ef.get(key(v, a)) || []).slice().sort().join(), f2 = (ef.get(key(v, b)) || []).slice().sort().join();
    if (f1 !== f2) continue;
    for (const f of def.faces.values()) for (const l of f.loops) { const i = l.indexOf(v); if (i >= 0) l.splice(i, 1); }
    eRemove(def, ve[0]); eRemove(def, ve[1]);
    eAdd(def, a, b, e1.soft && e2.soft, e1.curve === e2.curve ? e1.curve : undefined, e1.hidden && e2.hidden);
    def.verts.delete(v);
  }
}

/**
 * Remove the edges of `eids` that lie between two coplanar faces facing the same way, merging those faces. The edges are
 * taken out a plane at a time and each plane is re-faced once (kernel-print, 2026-10-23: one reface per edge made a
 * boolean result of a few thousand triangles take minutes); merging only joins faces of one plane, so an edge that
 * qualifies at the start still does after any other merge on that plane.
 */
export function mergeCoplanar(def, eids) {
  const ef = edgeFaces(def), planes = [];
  for (const eid of eids) {
    const e = def.edges.get(eid);
    if (!e) continue;
    const fs = ef.get(key(e.a, e.b)) || [];
    if (fs.length !== 2) continue;
    const f1 = def.faces.get(fs[0]), f2 = def.faces.get(fs[1]);
    const p1 = facePlane(def, f1), p2 = facePlane(def, f2);
    if (dot(p1.n, p2.n) < 1 - 1e-9 || Math.abs(p1.d - p2.d) > EPS) continue;
    const pl = planeFrom(p1.n, def.verts.get(e.a));
    let g = planes.find((q) => samePlane(q.pl, pl));
    if (!g) planes.push((g = { pl, eids: [], ends: [] }));
    g.eids.push(eid); g.ends.push(e.a, e.b);
  }
  for (const g of planes) {
    for (const eid of g.eids) eRemove(def, eid);
    reface(def, g.pl);
  }
  if (!planes.length) return;
  removeLooseVerts(def);
  healCollinear(def, [...new Set(planes.flatMap((g) => g.ends))]);
}