Measuring and mass properties in 3D

/** Centre, radius and plane normal of a circle edge of an exact regen (edges[i]: a, b, at; mesh segs give the normal). */
export function circleOf(e, build = null, index = -1) {
  const full = V.len(V.sub(e.a, e.b)) < 1e-9 * Math.max(1, V.len(e.a));
  let c;
  if (full) c = V.mul(V.add(e.a, e.at), 0.5);
  else {
    // circumcentre of a, at, b
    const A = e.a, B = e.at, C = e.b;
    const ab = V.sub(B, A), ac = V.sub(C, A), n = V.cross(ab, ac), nn = V.dot(n, n);
    if (nn < 1e-24) return null;
    c = V.add(A, V.mul(V.add(V.mul(V.cross(n, ab), V.dot(ac, ac)), V.mul(V.cross(ac, n), V.dot(ab, ab))), 1 / (2 * nn)));
  }
  const r = full ? V.len(V.sub(e.a, e.at)) / 2 : V.len(V.sub(e.a, c));
  let n = full ? null : V.norm(V.cross(V.sub(e.a, c), V.sub(e.b, c)));
  if (full || V.len(n) < 0.5) {
    // a third point off the a–at diameter from OCCT's polyline of this edge
    n = null;
    const segs = build?.mesh?.edgeSegs, se = build?.mesh?.segEdge;
    if (segs && se) {
      let best = 0;
      for (let k = 0; k < se.length; k++) {
        if (se[k] !== index) continue;
        const q = [segs[k * 6], segs[k * 6 + 1], segs[k * 6 + 2]];
        const x = V.cross(V.sub(e.a, c), V.sub(q, c));
        if (V.len(x) > best) { best = V.len(x); n = V.norm(x); }
      }
    }
  }
  return { c, r, n };
}

export function pointTri(p, [a, b, c]) {
  // closest point on a triangle (Ericson, Real-Time Collision Detection 5.1.5)
  const ab = V.sub(b, a), ac = V.sub(c, a), ap = V.sub(p, a);
  const d1 = V.dot(ab, ap), d2 = V.dot(ac, ap);
  if (d1 <= 0 && d2 <= 0) return a;
  const bp = V.sub(p, b), d3 = V.dot(ab, bp), d4 = V.dot(ac, bp);
  if (d3 >= 0 && d4 <= d3) return b;
  const vc = d1 * d4 - d3 * d2;
  if (vc <= 0 && d1 >= 0 && d3 <= 0) return V.add(a, V.mul(ab, d1 / (d1 - d3)));
  const cp = V.sub(p, c), d5 = V.dot(ab, cp), d6 = V.dot(ac, cp);
  if (d6 >= 0 && d5 <= d6) return c;
  const vb = d5 * d2 - d1 * d6;
  if (vb <= 0 && d2 >= 0 && d6 <= 0) return V.add(a, V.mul(ac, d2 / (d2 - d6)));
  const va = d3 * d6 - d5 * d4;
  if (va <= 0 && d4 - d3 >= 0 && d5 - d6 >= 0) return V.add(b, V.mul(V.sub(c, b), (d4 - d3) / ((d4 - d3) + (d5 - d6))));
  const den = 1 / (va + vb + vc);
  return V.add(a, V.add(V.mul(ab, vb * den), V.mul(ac, vc * den)));
}

export function segSeg([p1, q1], [p2, q2]) {
  // closest points of two segments (Ericson 5.1.9)
  const d1 = V.sub(q1, p1), d2 = V.sub(q2, p2), r = V.sub(p1, p2);
  const a = V.dot(d1, d1), e = V.dot(d2, d2), f = V.dot(d2, r);
  let s, t;
  const clamp = (x) => Math.max(0, Math.min(1, x));
  if (a <= 1e-18 && e <= 1e-18) return [p1, p2];
  if (a <= 1e-18) { s = 0; t = clamp(f / e); }
  else {
    const c = V.dot(d1, r);
    if (e <= 1e-18) { t = 0; s = clamp(-c / a); }
    else {
      const b = V.dot(d1, d2), den = a * e - b * b;
      s = den > 1e-18 ? clamp((b * f - c * e) / den) : 0;
      t = (b * s + f) / e;
      if (t < 0) { t = 0; s = clamp(-c / a); } else if (t > 1) { t = 1; s = clamp((b - c) / a); }
    }
  }
  return [V.add(p1, V.mul(d1, s)), V.add(p2, V.mul(d2, t))];
}

export function pair(A, B, ia, ib, exactDistance = null) {
  const out = {};
  if (isPointPick(A) && isPointPick(B)) Object.assign(out, { distance: V.len(V.sub(ib.p, ia.p)), p1: ia.p, p2: ib.p, method: 'exact' });
  else if (exactDistance) Object.assign(out, { distance: exactDistance.distance, p1: exactDistance.p1, p2: exactDistance.p2, method: 'exact', inside: exactDistance.inside });
  else if (A.geometry && B.geometry) {
    const d = meshDistance(A.geometry, B.geometry);
    Object.assign(out, d, { method: 'facets' });
  }
  if (out.p1) out.delta = V.sub(out.p2, out.p1);
  const da = dirOf(ia), db = dirOf(ib);
  if (da && db) {
    const c = V.dot(da.d, db.d);
    if (da.kind === 'plane' && db.kind === 'plane') out.angle = deg(acosC(c));
    else if (da.kind === 'line' && db.kind === 'line') out.angle = deg(acosC(Math.abs(c)));
    else out.angle = 90 - deg(acosC(Math.abs(c)));
  }
  // parallel planes: their gap; circles/holes: centre to centre (and axis to axis when parallel)
  if (da?.kind === 'plane' && db?.kind === 'plane' && Math.abs(Math.abs(V.dot(da.d, db.d)) - 1) < 1e-9) out.planeGap = Math.abs(V.dot(V.sub(ib.o, ia.o), da.d));
  if (ia.centre && ib.centre) out.centre = V.len(V.sub(ib.centre, ia.centre));
  const la = centreLine(ia), lb = centreLine(ib);
  if (la && lb && Math.abs(Math.abs(V.dot(la.d, lb.d)) - 1) < 1e-9) {
    const w = V.sub(lb.o, la.o);
    out.axisGap = V.len(V.sub(w, V.mul(la.d, V.dot(w, la.d))));
  }
  return out;
}

/** Mass-weighted properties in world axes: { mass (g), volume, area, centre, tensor (g·mm²) }. */
export function placeMass(props, matrix, density) {
  const M = matrix || IDENT;
  const R = [[M[0], M[4], M[8]], [M[1], M[5], M[9]], [M[2], M[6], M[10]]];
  const k = density / 1000;
  const T = props.tensor;
  const RT = [0, 1, 2].map((i) => [0, 1, 2].map((j) => R[i][0] * T[0][j] + R[i][1] * T[1][j] + R[i][2] * T[2][j]));
  const tensor = [0, 1, 2].map((i) => [0, 1, 2].map((j) => k * (RT[i][0] * R[j][0] + RT[i][1] * R[j][1] + RT[i][2] * R[j][2])));
  return { mass: massOf(props.volume, density), volume: props.volume, area: props.area, centre: V.applyMat(M, props.centre), tensor, exact: !!props.exact };
}

export function combineMass(list) {
  const mass = list.reduce((s, x) => s + x.mass, 0);
  const centre = [0, 1, 2].map((k) => list.reduce((s, x) => s + x.mass * x.centre[k], 0) / (mass || 1));
  const tensor = [[0, 0, 0], [0, 0, 0], [0, 0, 0]];
  for (const x of list) {
    const d = V.sub(x.centre, centre), dd = V.dot(d, d);
    for (let i = 0; i < 3; i++) for (let j = 0; j < 3; j++) tensor[i][j] += x.tensor[i][j] + x.mass * ((i === j ? dd : 0) - d[i] * d[j]);
  }
  return { mass, volume: list.reduce((s, x) => s + x.volume, 0), area: list.reduce((s, x) => s + x.area, 0), centre, tensor, exact: list.every((x) => x.exact) };
}

/** Eigenvalues (ascending) and unit eigenvectors of a symmetric 3×3 (cyclic Jacobi). */
export function principal(T) {
  const a = T.map((r) => [...r]);
  const v = [[1, 0, 0], [0, 1, 0], [0, 0, 1]];
  for (let sweep = 0; sweep < 50; sweep++) {
    const off = Math.abs(a[0][1]) + Math.abs(a[0][2]) + Math.abs(a[1][2]);
    if (off < 1e-15 * (Math.abs(a[0][0]) + Math.abs(a[1][1]) + Math.abs(a[2][2]) + 1e-300)) break;
    for (const [p, q] of [[0, 1], [0, 2], [1, 2]]) {
      if (Math.abs(a[p][q]) < 1e-300) continue;
      const th = (a[q][q] - a[p][p]) / (2 * a[p][q]);
      const t = Math.sign(th || 1) / (Math.abs(th) + Math.sqrt(th * th + 1));
      const c = 1 / Math.sqrt(t * t + 1), s = t * c;
      for (let k = 0; k < 3; k++) { const akp = a[k][p], akq = a[k][q]; a[k][p] = c * akp - s * akq; a[k][q] = s * akp + c * akq; }
      for (let k = 0; k < 3; k++) { const apk = a[p][k], aqk = a[q][k]; a[p][k] = c * apk - s * aqk; a[q][k] = s * apk + c * aqk; }
      for (let k = 0; k < 3; k++) { const vkp = v[k][p], vkq = v[k][q]; v[k][p] = c * vkp - s * vkq; v[k][q] = s * vkp + c * vkq; }
    }
  }
  const order = [0, 1, 2].sort((i, j) => a[i][i] - a[j][j]);
  return { moments: order.map((i) => a[i][i]), axes: order.map((i) => [v[0][i], v[1][i], v[2][i]]) };
}