Illustrate: projection, hidden lines and vector pictures

export function projector(cam, aspect = 1) {
  const dir = norm(sub(cam.target, cam.pos));
  let right = norm(cross(dir, cam.up));
  if (len(right) < 1e-9) right = norm(cross(dir, Math.abs(dir[1]) < 0.9 ? [0, 1, 0] : [1, 0, 0]));
  const up = cross(right, dir);
  const dist = len(sub(cam.target, cam.pos));
  if (cam.type === 'persp') {
    const tanH = Math.tan((cam.fov || 35) * Math.PI / 360);   // half the vertical field
    return {
      type: 'persp', aspect, right, up, dir, mmPerUnit: 2 * dist * tanH,
      project(p) {
        const q = sub(p, cam.pos), d = dot(q, dir);
        const k = 1 / (2 * tanH * Math.max(d, 1e-9));
        return [aspect / 2 + dot(q, right) * k, 0.5 + dot(q, up) * k, d];
      },
      ray(u, v) {
        const x = (u - aspect / 2) * 2 * tanH, y = (v - 0.5) * 2 * tanH;
        return { o: [...cam.pos], d: norm(add(dir, add(mul(right, x), mul(up, y)))) };
      },
    };
  }
  const h = cam.height || 100;
  return {
    type: 'ortho', aspect, right, up, dir, mmPerUnit: 2 * h,
    project(p) {
      const q = sub(p, cam.target);
      return [aspect / 2 + dot(q, right) / (2 * h), 0.5 + dot(q, up) / (2 * h), dot(sub(p, cam.pos), dir)];
    },
    ray(u, v) {
      const o = add(cam.target, add(mul(right, (u - aspect / 2) * 2 * h), mul(up, (v - 0.5) * 2 * h)));
      return { o: sub(o, mul(dir, dist + 1e6)), d: dir };
    },
  };
}

function occluded(u, v, d, self = -1) {
  const i = Math.min(G - 1, Math.max(0, Math.floor((u - u0) / cw))), j = Math.min(G - 1, Math.max(0, Math.floor((v - v0) / ch)));
  for (const k of grid[j * G + i]) {
    const p = parts[occ[k * 2]], tt = occ[k * 2 + 1], I = p.mesh.indices, Q = p.Q;
    const a = I[tt * 3] * 3, b = I[tt * 3 + 1] * 3, c = I[tt * 3 + 2] * 3;
    const x1 = Q[a], y1 = Q[a + 1], x2 = Q[b], y2 = Q[b + 1], x3 = Q[c], y3 = Q[c + 1];
    const den = (y2 - y3) * (x1 - x3) + (x3 - x2) * (y1 - y3);
    if (Math.abs(den) < 1e-18) continue;
    const l1 = ((y2 - y3) * (u - x3) + (x3 - x2) * (v - y3)) / den, l2 = ((y3 - y1) * (u - x3) + (x1 - x3) * (v - y3)) / den, l3 = 1 - l1 - l2;
    if (l1 < margin || l2 < margin || l3 < margin) continue;
    let w1 = l1, w2 = l2, w3 = l3, depth;
    if (persp) { w1 = l1 / Q[a + 2]; w2 = l2 / Q[b + 2]; w3 = l3 / Q[c + 2]; const s = w1 + w2 + w3; depth = 1 / s; w1 /= s; w2 /= s; w3 /= s; } else depth = l1 * Q[a + 2] + l2 * Q[b + 2] + l3 * Q[c + 2];
    if (depth >= d - eps * 10) continue;

export function featureSegments(mesh, { angle = 20, topo = null } = {}) {
  const { edges } = topo || topology(mesh);
  const P = mesh.positions, lim = Math.cos(angle * Math.PI / 180);
  const out = [];
  for (const e of edges) {
    if (e.t2 >= 0 && e.cos >= lim) continue;
    out.push(P[e.a * 3], P[e.a * 3 + 1], P[e.a * 3 + 2], P[e.b * 3], P[e.b * 3 + 1], P[e.b * 3 + 2]);
  }
  return Float32Array.from(out);
}

const num = (v) => { const r = Math.abs(v) < 1e-9 ? 0 : v; return String(Number(r.toFixed(3))); };