Push/pull, move, stretch, follow me and intersect

/**
 * Push or pull face `fid` by `distance` along its normal (positive = out of its front).
 * Returns {ok, mode: 'move'|'extrude', distance (after clamping), face (the moved/new face id)}.
 */
export function pushPull(def, fid, distance) {
  const f = def.faces.get(fid);
  if (!f || !(Math.abs(distance) > EPS)) return { ok: false };
  const n = faceNormal(def, f);
  if (!len(n)) return { ok: false };
  const ef = edgeFaces(def);
  const vids = new Set(f.loops.flat());
  const segs = loopSegs(f);
  const contains = (g) => Math.abs(dot(faceNormal(def, g), n)) < 1e-6;

  // Move mode: every boundary edge has a neighbour face that runs along n, and so does every
  // other face touching the face's corners.
  let move = segs.every(([a, b]) => {
    const others = (ef.get(key(a, b)) || []).filter((x) => x !== fid);
    return others.length > 0 && others.every((o) => contains(def.faces.get(o)));
  });
  if (move) {
    const vf = vertFaces(def);
    for (const v of vids) for (const g of vf.get(v) || []) if (g !== fid && !contains(def.faces.get(g))) move = false;
  }
  if (move) {
    // the face can come back no further than its shortest side
    let limit = Infinity;
    const ve = vertEdges(def);
    for (const v of vids) for (const eid of ve.get(v) || []) {
      const e = def.edges.get(eid), w = e.a === v ? e.b : e.a;
      if (vids.has(w)) continue;
      const along = dot(sub(def.verts.get(w), def.verts.get(v)), n);
      if (along < -EPS) limit = Math.min(limit, -along);
    }
    if (distance < 0 && -distance >= limit - EPS) distance = -limit;
    for (const v of vids) def.verts.set(v, add(def.verts.get(v), mul(n, distance)));
    const anchor = add(def.verts.get(f.loops[0][0]), [0, 0, 0]);
    cleanup(def);
    reface(def, planeFrom(n, anchor));
    cleanup(def);
    return { ok: true, mode: 'move', distance, face: def.faces.has(fid) ? fid : null };
  }

  const solid = backed(def, fid, ef);
  if (solid && distance < 0) {
    const t = depthBehind(def, fid);
    if (-distance > t - EPS) distance = -t;
  }
  const map = new Map();
  for (const v of vids) map.set(v, vAdd(def, add(def.verts.get(v), mul(n, distance))));
  const orig = f.loops.map((l) => l.slice());
  const top = orig.map((l) => l.map((v) => map.get(v)));
  const lone = !solid;
  const flip = lone && distance < 0;
  if (lone && distance > 0) f.loops = orig.map((l) => l.slice().reverse());
  else if (!lone) def.faces.delete(fid);
  const topId = fAdd(def, flip ? top.map((l) => l.slice().reverse()) : top, f.color);
  const sideEdges = [];
  for (const l of orig) {
    for (let i = 0; i < l.length; i++) {
      const a = l[i], b = l[(i + 1) % l.length], a2 = map.get(a), b2 = map.get(b);
      const quad = [a, b, b2, a2];
      fAdd(def, [flip ? quad.reverse() : quad], f.color);
      sideEdges.push(edgeBetween(def, a, a2), edgeBetween(def, a, b));
    }
  }
  cleanup(def);
  const topFace = def.faces.get(topId);
  if (topFace) reface(def, planeFrom(n, def.verts.get(topFace.loops[0][0])));
  cleanup(def);
  mergeCoplanar(def, sideEdges.filter((x) => x !== null));
  // the sides of a pulled circle read as one smooth surface: soften edges between side faces < 25° apart
  const ef2 = edgeFaces(def);
  for (const l of orig) for (const a of l) {
    const eid = edgeBetween(def, a, map.get(a));
    if (eid === null) continue;
    const e = def.edges.get(eid), fs = ef2.get(key(e.a, e.b)) || [];
    if (fs.length === 2 && dot(faceNormal(def, def.faces.get(fs[0])), faceNormal(def, def.faces.get(fs[1]))) > Math.cos(25 * Math.PI / 180)) e.soft = true;
  }
  return { ok: true, mode: 'extrude', distance, face: def.faces.has(topId) ? topId : null };
}

/** How far a face can be pushed in before it meets the first face behind it (Infinity if none). */
export function depthBehind(def, fid) {
  const f = def.faces.get(fid), n = faceNormal(def, f);
  const fr = frame(n);
  const pt2 = samplePoint(f.loops.map((l) => l.map((v) => to2(fr, def.verts.get(v)))));
  if (!pt2) return Infinity;
  const P0 = def.verts.get(f.loops[0][0]);
  const o = add(add(mul(fr.u, pt2[0]), mul(fr.v, pt2[1])), mul(n, dot(n, P0)));
  const d = mul(n, -1);
  let best = Infinity;
  for (const [gid, g] of def.faces) {
    if (gid === fid) continue;
    for (const [a, b, c] of faceTris(def, g)) {
      const t = rayTri(o, d, def.verts.get(a), def.verts.get(b), def.verts.get(c));
      if (t !== null && t > EPS && t < best) best = t;
    }
  }
  return best;
}

/** Move vertices by matrix m; faces bent out of plane fold into triangles. */
export function transformVerts(def, vs, m) {
  if (vs.size) {
    const vf = vertFaces(def);
    const touched = new Set();
    for (const v of vs) for (const f of vf.get(v) || []) touched.add(f);
    const before = new Map();
    for (const fid of touched) before.set(fid, faceTris(def, def.faces.get(fid)));
    for (const v of vs) def.verts.set(v, apply(m, def.verts.get(v)));
    // a mirror (Flip, a negative scale) turns faces inside out: those it moved whole are reversed so they face out again
    if (det3(m) < 0) for (const f of def.faces.values()) if (f.loops.every((l) => l.every((v) => vs.has(v)))) f.loops = f.loops.map((l) => l.slice().reverse());
    // a curve moved whole keeps its points (Move, Rotate, Scale, Flip); stretched in part, it is no longer editable
    for (const [cid, sp] of def.splines) {
      const own = [...def.edges.values()].filter((e) => e.curve === cid);
      const moved = own.filter((e) => vs.has(e.a) && vs.has(e.b)).length;
      if (own.length && moved === own.length) sp.pts = sp.pts.map((p) => apply(m, p));
      else if (moved || own.some((e) => vs.has(e.a) || vs.has(e.b))) def.splines.delete(cid);
    }
    for (const fid of touched) {
      const f = def.faces.get(fid);
      if (planar(def, f)) continue;
      def.faces.delete(fid);
      for (const tri of before.get(fid)) fAdd(def, [tri], f.color, undefined, f.hidden);
    }
    cleanup(def);
  }
}

function planar(def, f) {
  const pl = facePlane(def, f);
  if (!len(pl.n)) return false;
  return f.loops.every((l) => l.every((v) => Math.abs(dot(pl.n, def.verts.get(v)) - pl.d) <= EPS * 10));
}

/**
 * Stretch (resize after the fact): move every vertex of `def` at or beyond `cut` along unit direction `d` by
 * `delta` — only the vertices in `only`, when given — and every instance standing wholly beyond the cut on that
 * geometry. Faces between the moved and the still vertices stretch. Throws (changing nothing) when a face would
 * bend out of its plane or turn inside out. `boxOf(defId)` gives a definition's own box (to move instances).
 * Returns the number of vertices moved.
 */
export function stretch(def, d, cut, delta, { only = null, boxOf = () => null } = {}) {
  const tol = EPS * 10;
  const vs = new Set();
  for (const [id, p] of def.verts) if ((!only || only.has(id)) && dot(p, d) >= cut - tol) vs.add(id);
  if (!vs.size && !def.insts.size) throw new Error(t('m3d.model.there_is_nothing_on_that_side'));
  const at = new Map([...vs].map((v) => [v, add(def.verts.get(v), mul(d, delta))]));
  const pos = (v) => at.get(v) || def.verts.get(v);
  for (const f of def.faces.values()) {
    const all = f.loops.flat();
    const n = all.filter((v) => vs.has(v)).length;
    if (!n || n === all.length) continue;
    const n0 = faceNormal(def, f), n1 = newell(f.loops[0].map(pos));
    if (len(n1) <= EPS2 || dot(norm(n1), n0) <= 0) throw new Error(t('m3d.model.a_face_would_shrink_to_nothing'));
    const k = norm(n1), c = dot(k, pos(f.loops[0][0]));
    if (!all.every((v) => Math.abs(dot(k, pos(v)) - c) <= tol)) throw new Error(t('m3d.model.a_face_would_bend_out_of'));
  }
  let bb = null;
  for (const v of only || []) bb = bboxAdd(bb, def.verts.get(v));
  for (const [v, p] of at) def.verts.set(v, p);
  for (const inst of def.insts.values()) {
    const local = boxOf(inst.def);
    if (!local) continue;
    // the instance's corners in def's coordinates: it moves when it stands wholly beyond the cut, on this piece
    const corners = [];
    for (let i = 0; i < 8; i++) corners.push(apply(inst.m, [i & 1 ? local.max[0] : local.min[0], i & 2 ? local.max[1] : local.min[1], i & 4 ? local.max[2] : local.min[2]]));
    if (corners.some((p) => dot(p, d) < cut - tol)) continue;
    let ib = null;
    for (const p of corners) ib = bboxAdd(ib, p);
    if (bb && [0, 1, 2].some((i) => ib.min[i] > bb.max[i] + tol || ib.max[i] < bb.min[i] - tol)) continue;
    inst.m = mmul(translation(mul(d, delta)), inst.m);
  }
  cleanup(def);
  return vs.size;
}

for (let k = 1; k <= nSeg; k++) {
  const d = dirs[k - 1];
  let n, o;
  if (closed && k === nSeg) { n = np; o = p0; }   // back to the profile's own plane: the ring closes on it
  else if (!closed && k === nSeg) { n = d; o = P[k]; }
  else {
    const dn = dirs[k % nSeg];
    if (dot(d, dn) < -0.999) throw new Error(t('m3d.model.the_path_turns_straight_back_on'));
    n = norm(add(d, dn)); o = P[k];
  }
  const nd = dot(n, d);
  if (Math.abs(nd) < 1e-9) throw new Error(t('m3d.model.the_path_turns_too_sharply_here'));
  ring = ring.map((L) => L.map((q) => add(q, mul(d, dot(n, sub(o, q)) / nd))));
  if (closed && k === nSeg) rings.push(rings[0]);
  else rings.push(ring.map((L) => L.map((q) => vAdd(def, q))));
}

const spans = (F, x0, u) => {
  const w = cross(F.n, u), ts = [];
  for (const loop of F.loops) for (let i = 0; i < loop.length; i++) {
    const P = loop[i], Q = loop[(i + 1) % loop.length];
    const sp = dot(sub(P, x0), w), sq = dot(sub(Q, x0), w);
    if ((sp > 0) === (sq > 0)) continue;
    const c = add(P, mul(sub(Q, P), sp / (sp - sq)));
    ts.push(dot(sub(c, x0), u));
  }
  ts.sort((a, b) => a - b);
  const out = [];
  for (let i = 0; i + 1 < ts.length; i += 2) if (ts[i + 1] - ts[i] > EPS) out.push([ts[i], ts[i + 1]]);
  return out;
};

  const dir = cross(a.n, b.n), L = len(dir);
  if (L < 1e-9) continue;   // parallel or coplanar
  const u = mul(dir, 1 / L);
  const x0 = mul(cross(sub(mul(b.n, a.d), mul(a.n, b.d)), dir), 1 / (L * L));
  const sa = spans(a, x0, u), sb = spans(b, x0, u);
  for (const [a0, a1] of sa) for (const [b0, b1] of sb) {
    const t0 = Math.max(a0, b0), t1 = Math.min(a1, b1);
    if (t1 - t0 > EPS) segs.push([add(x0, mul(u, t0)), add(x0, mul(u, t1))]);
  }
}

/** The mirror through the plane at point c with unit normal n (an affine 4×4, column-major). */
export function mirrorMatrix(c, n) {
  const [x, y, z] = norm(n);
  const L = [1 - 2 * x * x, -2 * x * y, -2 * x * z, 0, -2 * x * y, 1 - 2 * y * y, -2 * y * z, 0, -2 * x * z, -2 * y * z, 1 - 2 * z * z, 0, 0, 0, 0, 1];
  return mmul(translation(c), mmul(L, translation(mul(c, -1))));
}