SKP import: placements, texture coordinates and face loops

/** SketchUp's 13 doubles → column-major 4×4. The first nine are the 3×3 part ROW by row (x' = m0·x + m1·y +
 * m2·z + m9, as OpenSKP reads it; checked on gondola_v20.skp, whose rotated components meet their neighbours'
 * vertices only this way round), then the translation, then w. w is 1 in every file we have; we treat it as
 * the homogeneous divisor. */
function matrixOf(xf) {
  const w = okNum(xf[12]) && Math.abs(xf[12]) > 1e-12 ? xf[12] : 1;
  const m = [xf[0], xf[3], xf[6], 0, xf[1], xf[4], xf[7], 0, xf[2], xf[5], xf[8], 0, xf[9], xf[10], xf[11], w];
  if (!m.every(okNum)) return null;
  if (w !== 1) { for (let i = 0; i < 15; i++) m[i] /= w; m[15] = 1; }
  return m;
}

/** Texture coordinates for a face's loops (OpenSKP's recipe): plane basis xr = Z × n, yr = n × xr; untouched
 * textures project straight onto it, positioned ones go through the inverse of their 3×3 mapping. */
function faceUv(loops, n, ftc, tileW, tileH) {
  const [nx, ny, nz] = n;
  let xr, yr;
  const cl = Math.hypot(nx, ny);
  if (cl < 1e-9) { xr = [1, 0, 0]; yr = [0, nz >= 0 ? 1 : -1, 0]; }
  else { xr = [-ny / cl, nx / cl, 0]; yr = [-nz * xr[1], nz * xr[0], nx * xr[1] - ny * xr[0]]; }
  let inv = null;
  if (ftc) {
    const [a, b, c, d, e, f, g, h, i] = ftc.front;
    const det = a * (e * i - f * h) - b * (d * i - f * g) + c * (d * h - e * g);
    if (Number.isFinite(det) && Math.abs(det) > 1e-12) {
      const k = 1 / det;
      inv = [(e * i - f * h) * k, (c * h - b * i) * k, (b * f - c * e) * k, (f * g - d * i) * k, (a * i - c * g) * k,
        (c * d - a * f) * k, (d * h - e * g) * k, (b * g - a * h) * k, (a * e - b * d) * k];
    }
  }
  const tw = okNum(tileW) && tileW > 1e-9 ? tileW : 1, th = okNum(tileH) && tileH > 1e-9 ? tileH : 1;
  return loops.map((pts) => pts.map((p) => {
    const px = p[0] * xr[0] + p[1] * xr[1] + p[2] * xr[2];
    const py = p[0] * yr[0] + p[1] * yr[1] + p[2] * yr[2];
    if (!inv) return [px / tw, py / th];
    let q = px * inv[2] + py * inv[5] + inv[8];
    if (Math.abs(q) < 1e-12) q = 1;
    return [(px * inv[0] + py * inv[3] + inv[6]) / q / tw, (px * inv[1] + py * inv[4] + inv[7]) / q / th];
  }));
}

/** Newell normal of a loop (for checking orientation against the stored plane). */
function newell(pts) {
  let x = 0, y = 0, z = 0;
  for (let i = 0; i < pts.length; i++) {
    const a = pts[i], b = pts[(i + 1) % pts.length];
    x += (a[1] - b[1]) * (a[2] + b[2]); y += (a[2] - b[2]) * (a[0] + b[0]); z += (a[0] - b[0]) * (a[1] + b[1]);
  }
  return [x, y, z];
}

const w = m13[12];
if (!(Number.isFinite(w) && Math.abs(w) > 1e-12)) throw damaged('ccskp.damaged.placement_scale');
const k = 1 / w;
for (let i = 0; i < 12; i++) if (!okNum(m13[i] * k)) throw damaged('ccskp.damaged.placement_numbers');
// rows r0 = m13[0..2], r1 = m13[3..5], r2 = m13[6..8]; column-major out: column j = (r0[j], r1[j], r2[j]).
m = [m13[0] * k, m13[3] * k, m13[6] * k, 0,
  m13[1] * k, m13[4] * k, m13[7] * k, 0,
  m13[2] * k, m13[5] * k, m13[8] * k, 0,
  m13[9] * k, m13[10] * k, m13[11] * k, 1];

/** Texture coordinates of p on a face with normal n (SketchUp's default projection, or the face's placement matrix). */
function faceUv(p, n, place, tw, th) {
  // The face's plane axes: x = normalize(Z × n), y = n × x (vertical normals: x = X, y = ±Y).
  let xr, yr;
  const cl = Math.hypot(n[0], n[1]);
  if (cl < 1e-9) { xr = [1, 0, 0]; yr = [0, n[2] >= 0 ? 1 : -1, 0]; } else {
    xr = [-n[1] / cl, n[0] / cl, 0];
    yr = [n[1] * xr[2] - n[2] * xr[1], n[2] * xr[0] - n[0] * xr[2], n[0] * xr[1] - n[1] * xr[0]];
  }
  const px = p[0] * xr[0] + p[1] * xr[1] + p[2] * xr[2];
  const py = p[0] * yr[0] + p[1] * yr[1] + p[2] * yr[2];
  if (!place) return [px / tw, py / th];
  // place = 3×3 row-major, texture → face plane, row vectors: [u v 1]·M = [px py 1]; invert it.
  const [a, b, c, d, e, f, g, h, i] = place;
  const det = a * (e * i - f * h) - b * (d * i - f * g) + c * (d * h - e * g);
  if (!(Math.abs(det) > 1e-12)) return [px / tw, py / th];
  const k = 1 / det;
  const inv = [(e * i - f * h) * k, (c * h - b * i) * k, (b * f - c * e) * k,
    (f * g - d * i) * k, (a * i - c * g) * k, (c * d - a * f) * k,
    (d * h - e * g) * k, (b * g - a * h) * k, (a * e - b * d) * k];
  const u = px * inv[0] + py * inv[3] + inv[6];
  const v = px * inv[1] + py * inv[4] + inv[7];
  let q = px * inv[2] + py * inv[5] + inv[8];
  if (Math.abs(q) < 1e-12) q = 1;
  return [u / q / tw, v / q / th];
}

const n = f.normal && Number.isFinite(f.normal[0] + f.normal[1] + f.normal[2]) ? f.normal : newell(loops[0]);
const nl = Math.hypot(n[0], n[1], n[2]) || 1;
const unit = [n[0] / nl, n[1] / nl, n[2] / nl];
const area = (l) => { const w = newell(l); return (w[0] * unit[0] + w[1] * unit[1] + w[2] * unit[2]) / 2; };
let outer = 0;
for (let i = 1; i < loops.length; i++) if (Math.abs(area(loops[i])) > Math.abs(area(loops[outer]))) outer = i;
if (outer !== 0) { const [o] = loops.splice(outer, 1); loops.unshift(o); }
if (area(loops[0]) < 0) { for (const l of loops) l.reverse(); flipped++; }
for (let i = 1; i < loops.length; i++) if (area(loops[i]) > 0) loops[i].reverse();