3DM surfaces into triangles

/** A mesh record's faces as triangles (a quad whose last two corners match is a triangle). */
export function meshTriangles(m) {
  const nv = m.nv;
  const positions = new Float64Array(nv * 3);
  for (let i = 0; i < nv * 3; i++) positions[i] = m.V[i];
  const idx = [];
  for (let f = 0; f < m.nf; f++) {
    const a = m.F[f * 4], b = m.F[f * 4 + 1], c = m.F[f * 4 + 2], d = m.F[f * 4 + 3];
    if (c === d) { if (a !== b && b !== c && a !== c) idx.push(a, b, c); }
    else { if (a !== b && b !== c && a !== c) idx.push(a, b, c); if (a !== c && c !== d && a !== d) idx.push(a, c, d); }
  }
  const out = { positions, indices: Uint32Array.from(idx) };
  if (m.N) out.normals = Float64Array.from(m.N);
  if (m.T) out.uvs = Float64Array.from(m.T);
  if (m.C) { const c = new Float32Array(nv * 4); for (let i = 0; i < nv * 4; i++) c[i] = m.C[i] / 255; for (let i = 0; i < nv; i++) c[i * 4 + 3] = 1; out.colors = c; }
  return out;
}

  done() {
    const P = this.P, n = P.length / 3;
    let lo = [Infinity, Infinity, Infinity], hi = [-Infinity, -Infinity, -Infinity];
    for (let i = 0; i < n; i++) for (let k = 0; k < 3; k++) { lo[k] = Math.min(lo[k], P[i * 3 + k]); hi[k] = Math.max(hi[k], P[i * 3 + k]); }
    const q = Math.max(Math.hypot(hi[0] - lo[0], hi[1] - lo[1], hi[2] - lo[2]) * 1e-9, 1e-12);
    const key = new Map(), map = new Int32Array(n), pos = [], nor = [];
    for (let i = 0; i < n; i++) {
      const k = Math.round(P[i * 3] / q) + ',' + Math.round(P[i * 3 + 1] / q) + ',' + Math.round(P[i * 3 + 2] / q);
      let j = key.get(k);
      if (j === undefined) { j = pos.length / 3; key.set(k, j); pos.push(P[i * 3], P[i * 3 + 1], P[i * 3 + 2]); nor.push(this.N[i * 3], this.N[i * 3 + 1], this.N[i * 3 + 2]); }
      map[i] = j;
    }
    const I = [], TF = [];
    for (let t = 0; t < this.I.length; t += 3) {
      const a = map[this.I[t]], b = map[this.I[t + 1]], c = map[this.I[t + 2]];
      if (a === b || b === c || a === c) continue;
      I.push(a, b, c); TF.push(this.TF[t / 3]);
    }
    return { positions: Float64Array.from(pos), normals: Float32Array.from(nor), indices: Uint32Array.from(I), triFace: Int32Array.from(TF) };
  }
}

/** An edge's samples: { s: fractions along the edge's own direction, pts } on its 3D curve, or null when it has none we read. */
function edgeSamples(b, ei, tol, cache) {
  if (cache.has(ei)) return cache.get(ei);
  let r = null;
  const e = b.E[ei];
  const c = e && e.c3i >= 0 ? b.c3[e.c3i] : null;
  if (c && curveOk(c)) {
    const [p0, p1] = e.proxy;
    const { pts, ts } = sampleCurve(c, tol, [p0, p1]);
    const span = p1 - p0 || 1;
    let s = ts.map((t) => (t - p0) / span), P = pts;
    if (e.rev) { s = s.map((x) => 1 - x).reverse(); P = P.slice().reverse(); }
    r = { s, pts: P };
  }
  cache.set(ei, r);
  return r;
}

// triangulate in (u, v), scaled so a step in u and one in v are about as long on the surface
const du = dom[0][1] - dom[0][0] || 1, dv = dom[1][1] - dom[1][0] || 1;
const along = (fn) => { let L = 0, prev = null; for (let k = 0; k <= 8; k++) { const p = fn(k / 8); if (prev) L += len(sub(p, prev)); prev = p; } return L; };
const vm = dom[1][0] + dv / 2, um = dom[0][0] + du / 2;
const Lu = Math.max(along((s) => surfaceAt(S, dom[0][0] + s * du, vm)), 1e-12), Lv = Math.max(along((s) => surfaceAt(S, um, dom[1][0] + s * dv)), 1e-12);
const ku = Lu / du, kv = Lv / dv;

const half = (a, c) => {
  const key = a < c ? a + ',' + c : c + ',' + a;
  if (mids.has(key)) return mids.get(key);
  let m = -1;
  if (free || !(a < nb && c < nb && boundary.has(Math.min(a, c) * nb + Math.max(a, c)))) {
    const pa = pos(a), pc = pos(c);
    if (len(sub(pa, pc)) > minLen) {
      const u = (vu[a] + vu[c]) / 2, v = (vv[a] + vv[c]) / 2;
      const pm = surfaceAt(S, u, v);
      const chord = [(pa[0] + pc[0]) / 2, (pa[1] + pc[1]) / 2, (pa[2] + pc[2]) / 2];
      const na = vn[a], nc = vn[c];
      if (len(sub(pm, chord)) > tol || (na && nc && dot(na, nc) < cosMax)) {
        if (out.I.length / 3 >= maxTriangles) capped = true;
        else { m = out.add(pm, nrm(u, v)) - base; vu.push(u); vv.push(v); vn.push(nrm(u, v)); }
      }
    }
  }
  mids.set(key, m);
  return m;
};

/**
 * The grid step (in the scaled (u, v) plane, about surface length) for a face: from the most the surface bends, so a chord
 * stays within `tol` (sagitta L²κ/8) and turns at most ~20°; null for a flat face.
 */
function gridStep(S, dom, tol, Lu, Lv) {
  const n = 8, N = [], Pp = [];
  for (let i = 0; i <= n; i++) for (let j = 0; j <= n; j++) {
    const u = dom[0][0] + (dom[0][1] - dom[0][0]) * i / n, v = dom[1][0] + (dom[1][1] - dom[1][0]) * j / n;
    N.push(normalAt(S, u, v, dom)); Pp.push(surfaceAt(S, u, v));
  }
  const ks = [];
  const at = (i, j) => i * (n + 1) + j;
  const near = Math.min(Lu, Lv) / n / 4;   // pairs closer than this sit at a pole or a pinch: their normals say nothing
  for (let i = 0; i <= n; i++) for (let j = 0; j <= n; j++) for (const [a, b] of [[i + 1, j], [i, j + 1]]) {
    if (a > n || b > n) continue;
    const n0 = N[at(i, j)], n1 = N[at(a, b)];
    const d = len(sub(Pp[at(i, j)], Pp[at(a, b)]));
    if (!n0 || !n1 || d <= near) continue;
    const ang = Math.acos(Math.max(-1, Math.min(1, dot(n0, n1))));
    ks.push(ang / d);
  }
  // the grid follows how the face bends over most of it; the edge halving after it handles the tight places (an apex)
  ks.sort((a, b) => a - b);
  const k = ks.length ? ks[Math.floor(ks.length * 0.75)] : 0;
  if (k < 1e-9) return null;
  let h = Math.min(Math.sqrt(8 * tol / k), 0.349 / k);
  const big = Math.max(Lu, Lv);
  h = Math.max(h, big / 400, Math.sqrt((Lu * Lv) / MAX_GRID));
  return h < big ? h : null;
}

/** Grid points inside the face's loops (even-odd), at least half a step from every side. */
function insidePoints(loops, ku, kv, h, dom) {
  const segs = [];
  for (const l of loops) for (let k = 0; k < l.uv.length; k++) { const a = l.uv[k], b = l.uv[(k + 1) % l.uv.length]; segs.push([a[0] * ku, a[1] * kv, b[0] * ku, b[1] * kv]); }
  const x0 = dom[0][0] * ku, x1 = dom[0][1] * ku, y0 = dom[1][0] * kv, y1 = dom[1][1] * kv;
  const lo = [Math.min(x0, x1), Math.min(y0, y1)], hi = [Math.max(x0, x1), Math.max(y0, y1)];
  const nx = Math.floor((hi[0] - lo[0]) / h), ny = Math.floor((hi[1] - lo[1]) / h);
  if (nx < 1 && ny < 1) return [];
  const gx = (hi[0] - lo[0]) / (nx + 1), gy = (hi[1] - lo[1]) / (ny + 1);
  // the sides in buckets by row, for the distance test
  const rows = Math.max(1, ny + 1), bucket = Array.from({ length: rows + 2 }, () => []);
  const rowOf = (y) => Math.max(0, Math.min(rows + 1, Math.floor((y - lo[1]) / gy)));
  segs.forEach((sg, k) => { const a = rowOf(Math.min(sg[1], sg[3])), b = rowOf(Math.max(sg[1], sg[3])); for (let r = a; r <= b; r++) bucket[r].push(k); });
  const out = [];
  for (let j = 1; j <= ny; j++) {
    const y = lo[1] + j * gy;
    const r = rowOf(y);
    const near = new Set([...(bucket[r - 1] || []), ...bucket[r], ...(bucket[r + 1] || [])]);
    // crossings of this row by every side (even-odd)
    const xs = [];
    for (const sg of segs) { const [ax, ay, bx, by] = sg; if ((ay > y) !== (by > y)) xs.push(ax + (y - ay) * (bx - ax) / (by - ay)); }
    xs.sort((a, b) => a - b);
    for (let i = 1; i <= nx; i++) {
      const x = lo[0] + i * gx;
      let c = 0;
      while (c < xs.length && xs[c] < x) c++;
      if (c % 2 === 0) continue;
      let ok = true;
      for (const k of near) {
        const [ax, ay, bx, by] = segs[k];
        const dx = bx - ax, dy = by - ay, L = dx * dx + dy * dy;
        const t = L ? Math.max(0, Math.min(1, ((x - ax) * dx + (y - ay) * dy) / L)) : 0;
        if (Math.hypot(x - ax - t * dx, y - ay - t * dy) < 0.5 * Math.min(gx, gy)) { ok = false; break; }
      }
      if (ok) out.push([x, y]);
    }
  }
  return out;
}

/** An extrusion: walls along the path from each profile, caps from the closed ones (outer with its holes). */
export function tessellateExtrusion(e, { tol = 0.01, mode = 'cached-first', stats = newStats() } = {}) {
  if (mode !== 'ours' && e.mesh) { const out = new Builder(); out.mesh(e.mesh, 0); stats.cached++; return out.done(); }
  const profiles = e.count > 1 && e.profile && e.profile.cls === 'polyCurve' ? e.profile.segs : e.profile ? [e.profile] : [];
  if (!profiles.length || !profiles.every((p) => curveOk(p))) { bump(stats, 'notRead', 'extrusion'); return null; }
  const T = unit(sub(e.path[1], e.path[0])), Uu = unit(e.up);
  if (!T || !Uu) { bump(stats, 'notRead', 'extrusion'); return null; }
  const X = unit(cross(Uu, T));
  if (!X) { bump(stats, 'notRead', 'extrusion'); return null; }
  const Y = cross(T, X);
  const at = (s) => [0, 1, 2].map((k) => e.path[0][k] + (e.path[1][k] - e.path[0][k]) * s);
  const P0 = at(e.t[0]), P1 = at(e.t[1]);
  const place = (O, q) => [O[0] + q[0] * X[0] + q[1] * Y[0], O[1] + q[0] * X[1] + q[1] * Y[1], O[2] + q[0] * X[2] + q[1] * Y[2]];
  if (e.mitered) stats.mitered++;
  const out = new Builder();
  const rings = [];
  profiles.forEach((pc, k) => {
    const { pts } = sampleCurve(pc, tol);
    const closed = pts.length > 2 && Math.hypot(pts[0][0] - pts[pts.length - 1][0], pts[0][1] - pts[pts.length - 1][1]) <= Math.max(tol, 1e-9);
    let ring = closed ? pts.slice(0, -1) : pts.slice();
    const area = ring.reduce((A, p, i) => { const q = ring[(i + 1) % ring.length]; return A + p[0] * q[1] - q[0] * p[1]; }, 0);
    if (closed && ((k === 0) !== (area > 0))) ring = ring.reverse();   // outer turns left, holes right: walls face out
    rings.push({ ring, closed });
  });
  let faceId = 0;
  for (const { ring, closed } of rings) {
    const n = ring.length, segs = closed ? n : n - 1;
    // a new face where the profile turns sharply (a corner of a rectangle), one face along a smooth run
    let prevDir = null;
    for (let i = 0; i < segs; i++) {
      const a = ring[i], c = ring[(i + 1) % n];
      const dir = unit([c[0] - a[0], c[1] - a[1], 0]);
      if (!dir) continue;
      if (!prevDir || dot(dir, prevDir) < Math.cos(30 * Math.PI / 180)) faceId++;
      prevDir = dir;
      const nn = [dir[1] * X[0] - dir[0] * Y[0], dir[1] * X[1] - dir[0] * Y[1], dir[1] * X[2] - dir[0] * Y[2]];
      const i0 = out.add(place(P0, a), nn), i1 = out.add(place(P0, c), nn), i2 = out.add(place(P1, c), nn), i3 = out.add(place(P1, a), nn);
      out.tri(i0, i1, i2, faceId); out.tri(i0, i2, i3, faceId);
    }
  }
  const capRings = rings.filter((r) => r.closed);
  const capped = e.oldCaps ? rings[0] && rings[0].closed : true;
  if (capRings.length && rings[0].closed && capped) {
    const flat = [], holes = [], pts = [];
    capRings.forEach((r, k) => { if (k) holes.push(pts.length); for (const p of r.ring) { flat.push(p[0], p[1]); pts.push(p); } });
    const tris = earcut(flat, holes);
    for (const [end, O, nn] of [[0, P0, T.map((x) => -x)], [1, P1, T]]) {
      if (!e.cap[end] && !e.oldCaps) continue;
      faceId++;
      const base = out.nv;
      for (const p of pts) out.add(place(O, p), nn);
      for (let k = 0; k < tris.length; k += 3) {
        const [a, c, d] = [tris[k], tris[k + 1], tris[k + 2]];
        const ar = (pts[c][0] - pts[a][0]) * (pts[d][1] - pts[a][1]) - (pts[d][0] - pts[a][0]) * (pts[c][1] - pts[a][1]);
        const up = (ar > 0) === (end === 1);
        if (up) out.tri(base + a, base + c, base + d, faceId); else out.tri(base + a, base + d, base + c, faceId);
      }
    }
  }
  bump(stats, 'ours', 'extrusion');
  return out.done();
}

const baseTol = tol || model.tolerance || 0.01;
const tolFor = (bbox) => {
  if (!bbox || !(bbox[0][0] <= bbox[1][0] && bbox[0][1] <= bbox[1][1] && bbox[0][2] <= bbox[1][2])) return baseTol;
  const d = Math.hypot(bbox[1][0] - bbox[0][0], bbox[1][1] - bbox[0][1], bbox[1][2] - bbox[0][2]);
  return Number.isFinite(d) && d > 0 ? Math.max(baseTol, d * 1.5e-3) : baseTol;
};

/**
 * A box around a brep or extrusion from what is cheap to reach (its vertices, its curves' break points, control points),
 * for the tessellation tolerance only. The box a file stores is not used: some programs write an empty one.
 */
function roughBox(g) {
  const lo = [Infinity, Infinity, Infinity], hi = [-Infinity, -Infinity, -Infinity];
  const add = (p) => { if (p && p.every(Number.isFinite)) for (let k = 0; k < 3; k++) { lo[k] = Math.min(lo[k], p[k]); hi[k] = Math.max(hi[k], p[k]); } };
  const curve = (c) => { if (curveOk(c)) for (const p of sampleCurve(c, Infinity).pts) add(p); };
  if (g.cls === 'brep') {
    g.V.forEach(add);
    g.c3.forEach(curve);
    for (const s of g.S) {
      if (!s) continue;
      if (s.cls === 'nurbsSurface') { const w = s.dim + (s.rat ? 1 : 0); for (let i = 0; i < s.cvs.length; i += w) add(s.rat ? [s.cvs[i] / s.cvs[i + s.dim], s.cvs[i + 1] / s.cvs[i + s.dim], s.cvs[i + 2] / s.cvs[i + s.dim]] : [s.cvs[i], s.cvs[i + 1], s.dim > 2 ? s.cvs[i + 2] : 0]); }
      if (s.cls === 'revSurface' && curveOk(s.curve)) {
        const pts = sampleCurve(s.curve, Infinity).pts, [o, e] = s.axis, d = [e[0] - o[0], e[1] - o[1], e[2] - o[2]], L = Math.hypot(...d) || 1, z = d.map((x) => x / L);
        for (const p of pts) { const v = [p[0] - o[0], p[1] - o[1], p[2] - o[2]], t = v[0] * z[0] + v[1] * z[1] + v[2] * z[2], r = Math.hypot(v[0] - t * z[0], v[1] - t * z[1], v[2] - t * z[2]); add([o[0] + t * z[0] - r, o[1] + t * z[1] - r, o[2] + t * z[2] - r]); add([o[0] + t * z[0] + r, o[1] + t * z[1] + r, o[2] + t * z[2] + r]); }
      }
    }
  } else if (g.cls === 'extrusion') {
    add(g.path[0]); add(g.path[1]);
    if (curveOk(g.profile)) { const pts = sampleCurve(g.profile, Infinity).pts; let r = 0; for (const p of pts) r = Math.max(r, Math.hypot(p[0], p[1])); add(g.path[0].map((x) => x - r)); add(g.path[1].map((x) => x + r)); }
  }
  return lo[0] <= hi[0] ? [lo, hi] : null;
}