yawn/server/web/assets/sculpt/topology.js

import { buildMesh, sphere } from "./mesh.js";

export const MAX_VERTICES = 65536;
export const MAX_TRIANGLES = 131072;

function checked(positions, indices, level = 0, masks) {
  if (positions.length / 3 > MAX_VERTICES || indices.length / 3 > MAX_TRIANGLES)
    throw new Error("Mesh exceeds topology limits.");
  return buildMesh(positions, indices, level, masks);
}

function gridSurface(slices, stacks, point) {
  const p = [],
    i = [];
  for (let y = 0; y < stacks; y++)
    for (let x = 0; x < slices; x++) p.push(...point(x / slices, y / stacks));
  const at = (x, y) => y * slices + ((x + slices) % slices);
  for (let y = 0; y < stacks; y++)
    for (let x = 0; x < slices; x++) {
      const a = at(x, y),
        b = at(x + 1, y),
        c = at(x + 1, (y + 1) % stacks),
        d = at(x, (y + 1) % stacks);
      i.push(a, b, c, a, c, d);
    }
  return [p, i];
}

function cubeMesh(n) {
  const p = [],
    faces = [],
    map = new Map();
  const vertex = (v) => {
    const key = v.join(",");
    if (!map.has(key)) {
      map.set(key, p.length / 3);
      p.push(...v.map((x) => x / n));
    }
    return map.get(key);
  };
  const sides = [
    [
      [n, -n, -n],
      [0, 2 * n, 0],
      [0, 0, 2 * n],
    ],
    [
      [-n, -n, n],
      [0, 2 * n, 0],
      [0, 0, -2 * n],
    ],
    [
      [-n, n, -n],
      [0, 0, 2 * n],
      [2 * n, 0, 0],
    ],
    [
      [-n, -n, n],
      [0, 0, -2 * n],
      [2 * n, 0, 0],
    ],
    [
      [-n, -n, -n],
      [0, 2 * n, 0],
      [2 * n, 0, 0],
    ],
    [
      [n, -n, n],
      [0, 2 * n, 0],
      [-2 * n, 0, 0],
    ],
  ];
  for (const [o, u, v] of sides)
    for (let y = 0; y < n; y++)
      for (let x = 0; x < n; x++) {
        const q = (xx, yy) =>
          vertex(o.map((z, k) => z + (u[k] * xx) / n + (v[k] * yy) / n));
        const a = q(x, y),
          b = q(x + 1, y),
          c = q(x + 1, y + 1),
          d = q(x, y + 1);
        faces.push(a, b, c, a, c, d);
      }
  return [p, faces];
}

function radial(kind, slices, rings) {
  const p = [],
    idx = [];
  const cone = kind === "cone";
  for (let y = 0; y <= (cone ? rings - 1 : rings); y++) {
    const t = y / rings,
      r = cone ? 1 - t : 1,
      z = -1 + 2 * t;
    for (let x = 0; x < slices; x++) {
      const a = (2 * Math.PI * x) / slices;
      p.push(r * Math.cos(a), r * Math.sin(a), z);
    }
  }
  for (let y = 0; y < rings - (cone ? 1 : 0); y++)
    for (let x = 0; x < slices; x++) {
      const a = y * slices + x,
        b = y * slices + ((x + 1) % slices),
        c = (y + 1) * slices + ((x + 1) % slices),
        d = (y + 1) * slices + x;
      idx.push(a, b, c, a, c, d);
    }
  if (cone) {
    const apex = p.length / 3,
      o = (rings - 1) * slices;
    p.push(0, 0, 1);
    for (let x = 0; x < slices; x++)
      idx.push(o + x, o + ((x + 1) % slices), apex);
  }
  const bottom = p.length / 3;
  p.push(0, 0, -1);
  for (let x = 0; x < slices; x++) idx.push(bottom, (x + 1) % slices, x);
  if (!cone) {
    const top = p.length / 3;
    p.push(0, 0, 1);
    const o = rings * slices;
    for (let x = 0; x < slices; x++)
      idx.push(top, o + x, o + ((x + 1) % slices));
  }
  return [p, idx];
}

export function preset(name, level = 4) {
  if (!Number.isInteger(level) || level < 1 || level > 6)
    throw new Error("Choose a supported resolution.");
  if (name === "sphere") return sphere(level);
  let data;
  const scale = 2 ** level;
  if (name === "cube") data = cubeMesh(scale);
  else if (name === "cylinder" || name === "cone")
    data = radial(name, scale * 4, scale + 8);
  else if (name === "torus")
    data = gridSurface(scale * 4, scale * 2, (u, v) => {
      const a = 2 * Math.PI * u,
        b = 2 * Math.PI * v,
        r = 0.32,
        R = 0.68;
      return [
        (R + r * Math.cos(b)) * Math.cos(a),
        (R + r * Math.cos(b)) * Math.sin(a),
        r * Math.sin(b),
      ];
    });
  else throw new Error(`Unknown topology preset: ${name}`);
  return checked(data[0], data[1], level);
}

export function subdivide(mesh) {
  const tri = mesh.indices.length / 3,
    vertices = mesh.positions.length / 3;
  const edges = new Map();
  for (let k = 0; k < mesh.indices.length; k += 3)
    for (const [a, b] of [
      [mesh.indices[k], mesh.indices[k + 1]],
      [mesh.indices[k + 1], mesh.indices[k + 2]],
      [mesh.indices[k + 2], mesh.indices[k]],
    ])
      edges.set(a < b ? `${a},${b}` : `${b},${a}`, 0);
  if (vertices + edges.size > MAX_VERTICES || tri * 4 > MAX_TRIANGLES)
    throw new Error("Subdivision exceeds topology limits.");
  const p = Array.from(mesh.positions),
    masks = Array.from(mesh.masks || new Float32Array(vertices));
  for (const key of edges.keys()) {
    const [a, b] = key.split(",").map(Number),
      v = p.length / 3;
    edges.set(key, v);
    for (let j = 0; j < 3; j++) p.push((p[a * 3 + j] + p[b * 3 + j]) / 2);
    masks.push((masks[a] + masks[b]) / 2);
  }
  const mid = (a, b) => edges.get(a < b ? `${a},${b}` : `${b},${a}`),
    out = [];
  for (let k = 0; k < mesh.indices.length; k += 3) {
    const a = mesh.indices[k],
      b = mesh.indices[k + 1],
      c = mesh.indices[k + 2],
      ab = mid(a, b),
      bc = mid(b, c),
      ca = mid(c, a);
    out.push(a, ab, ca, b, bc, ab, c, ca, bc, ab, bc, ca);
  }
  return checked(p, out, mesh.level + 1, masks);
}

function topology(indices, count) {
  const neighbors = Array.from({ length: count }, () => new Set()),
    incident = Array.from({ length: count }, () => []),
    edgeFaces = new Map(),
    faces = new Set();
  let valid = true;
  for (let k = 0; k < indices.length; k += 3) {
    const face = [indices[k], indices[k + 1], indices[k + 2]],
      faceKey = [...face].sort((a, b) => a - b).join(",");
    if (new Set(face).size !== 3 || faces.has(faceKey)) valid = false;
    faces.add(faceKey);
    for (const v of face) incident[v].push(k / 3);
    for (const [a, b] of [
      [indices[k], indices[k + 1]],
      [indices[k + 1], indices[k + 2]],
      [indices[k + 2], indices[k]],
    ]) {
      neighbors[a].add(b);
      neighbors[b].add(a);
      const key = a < b ? `${a},${b}` : `${b},${a}`;
      const list = edgeFaces.get(key) || [];
      list.push(k / 3);
      edgeFaces.set(key, list);
      if (list.length > 2) valid = false;
    }
  }

  // A manifold vertex link is one cycle (interior) or one path (boundary).
  // Edge incidence alone does not detect two components touching at a vertex.
  for (let v = 0; valid && v < count; v++) {
    if (!incident[v].length) continue;
    const link = new Map();
    for (const f of incident[v]) {
      const other = indices.slice(f * 3, f * 3 + 3).filter((x) => x !== v);
      if (other.length !== 2) {
        valid = false;
        break;
      }
      if (!link.has(other[0])) link.set(other[0], new Set());
      if (!link.has(other[1])) link.set(other[1], new Set());
      link.get(other[0]).add(other[1]);
      link.get(other[1]).add(other[0]);
    }
    if (!valid) break;
    const boundary = [...neighbors[v]].filter((n) => {
      const key = v < n ? `${v},${n}` : `${n},${v}`;
      return edgeFaces.get(key).length === 1;
    });
    const degrees = [...link.values()].map((x) => x.size);
    if (
      degrees.some((d) => d < 1 || d > 2) ||
      (boundary.length === 0
        ? degrees.some((d) => d !== 2)
        : boundary.length !== 2 || degrees.filter((d) => d === 1).length !== 2)
    ) {
      valid = false;
      break;
    }
    const reached = new Set(),
      pending = [link.keys().next().value];
    while (pending.length) {
      const x = pending.pop();
      if (reached.has(x)) continue;
      reached.add(x);
      for (const n of link.get(x)) pending.push(n);
    }
    if (reached.size !== link.size) valid = false;
  }
  return { neighbors, incident, edgeFaces, valid };
}

export function decimate(mesh, ratio = 0.5) {
  if (!Number.isFinite(ratio) || ratio <= 0 || ratio > 1)
    throw new Error("Ratio must be in (0, 1].");
  let p = Array.from(mesh.positions),
    idx = Array.from(mesh.indices),
    masks = Array.from(mesh.masks || new Float32Array(p.length / 3));
  const target = Math.max(4, Math.floor((idx.length / 3) * ratio)),
    deadline = Date.now() + 2500;
  let accepted = 0,
    work = 0;
  const MAX_PASSES = 64,
    MAX_WORK = 1000000;
  for (
    let pass = 0;
    pass < MAX_PASSES && idx.length / 3 > target && Date.now() < deadline;
    pass++
  ) {
    const { neighbors, incident, edgeFaces, valid } = topology(
        idx,
        p.length / 3,
      ),
      protectedVertices = new Set();
    if (!valid) break;
    for (const [key, faces] of edgeFaces)
      if (faces.length === 1)
        for (const v of key.split(",").map(Number)) protectedVertices.add(v);
    masks.forEach((m, v) => {
      if (m !== 0) {
        protectedVertices.add(v);
        for (const n of neighbors[v]) protectedVertices.add(n);
      }
    });
    const candidates = [];
    for (const [key, edgeIncident] of edgeFaces) {
      if (Date.now() >= deadline || work++ >= MAX_WORK) break;
      if (edgeIncident.length === 2) {
        const [a, b] = key.split(",").map(Number);
        if (protectedVertices.has(a) || protectedVertices.has(b)) continue;
        const opposite = edgeIncident.map((f) =>
            idx.slice(f * 3, f * 3 + 3).find((v) => v !== a && v !== b),
          ),
          common = [...neighbors[a]].filter((v) => neighbors[b].has(v));
        if (
          new Set(opposite).size !== 2 ||
          common.length !== 2 ||
          common.some((v) => !opposite.includes(v))
        )
          continue;
        let d = 0;
        for (let j = 0; j < 3; j++) d += (p[a * 3 + j] - p[b * 3 + j]) ** 2;
        candidates.push([d, a, b]);
      }
    }
    candidates.sort((a, b) => a[0] - b[0] || a[1] - b[1] || a[2] - b[2]);
    const replace = new Map(),
      locked = new Set();
    let removed = 0;
    for (const [, a, b] of candidates) {
      if (Date.now() >= deadline || work++ >= MAX_WORK) break;
      if (idx.length / 3 - removed * 2 <= target) break;
      if (locked.has(a) || locked.has(b)) continue;
      const pos = [0, 1, 2].map((j) =>
          Math.fround((p[a * 3 + j] + p[b * 3 + j]) / 2),
        ),
        localFaces = [...new Set([...incident[a], ...incident[b]])],
        edgeKey = a < b ? `${a},${b}` : `${b},${a}`,
        omitted = new Set(edgeFaces.get(edgeKey)),
        seen = new Set();
      let candidateValid = true;
      for (const f of localFaces) {
        if (Date.now() >= deadline || work++ >= MAX_WORK) {
          candidateValid = false;
          break;
        }
        const k = f * 3;
        if (omitted.has(f)) continue;
        const vs = [idx[k], idx[k + 1], idx[k + 2]].map((v) =>
          v === b ? a : v,
        );
        const triangleKey = [...vs].sort((x, y) => x - y).join(",");
        if (new Set(vs).size < 3 || seen.has(triangleKey)) {
          candidateValid = false;
          break;
        }
        seen.add(triangleKey);
        const q = vs.map((v) => (v === a ? pos : p.slice(v * 3, v * 3 + 3)));
        const ux = q[1][0] - q[0][0],
          uy = q[1][1] - q[0][1],
          uz = q[1][2] - q[0][2],
          vx = q[2][0] - q[0][0],
          vy = q[2][1] - q[0][1],
          vz = q[2][2] - q[0][2];
        const nx = uy * vz - uz * vy,
          ny = uz * vx - ux * vz,
          nz = ux * vy - uy * vx;
        const ox = meshLikeNormal(p, idx, k);
        if (
          nx * nx + ny * ny + nz * nz < 1e-16 ||
          nx * ox[0] + ny * ox[1] + nz * ox[2] <= 0
        ) {
          candidateValid = false;
          break;
        }
      }
      if (!candidateValid) continue;
      replace.set(b, a);
      p.splice(a * 3, 3, ...pos);
      locked.add(a);
      locked.add(b);
      for (const n of neighbors[a]) locked.add(n);
      for (const n of neighbors[b]) locked.add(n);
      removed++;
      accepted++;
    }
    if (!replace.size) break;
    idx = idx.map((v) => replace.get(v) ?? v);
    const seen = new Set();
    idx = idx.reduce((o, v, k) => {
      if (k % 3) return o;
      const f = idx.slice(k, k + 3);
      if (new Set(f).size < 3) return o;
      const key = [...f].sort((a, b) => a - b).join(",");
      // Candidate validation guarantees this. Never use duplicate removal to
      // repair a collapse which violated the link condition.
      if (seen.has(key))
        throw new Error("Decimation produced duplicate faces.");
      seen.add(key);
      o.push(...f);
      return o;
    }, []);
  }
  if (!accepted)
    return checked(mesh.positions, mesh.indices, mesh.level, mesh.masks);
  const used = [...new Set(idx)].sort((a, b) => a - b),
    remap = new Map(used.map((v, i) => [v, i])),
    outp = [],
    outm = [];
  for (const v of used) {
    outp.push(...p.slice(v * 3, v * 3 + 3));
    outm.push(masks[v]);
  }
  return checked(
    outp,
    idx.map((v) => remap.get(v)),
    mesh.level,
    outm,
  );
}

function meshLikeNormal(p, idx, k) {
  const a = idx[k],
    b = idx[k + 1],
    c = idx[k + 2],
    ux = p[b * 3] - p[a * 3],
    uy = p[b * 3 + 1] - p[a * 3 + 1],
    uz = p[b * 3 + 2] - p[a * 3 + 2],
    vx = p[c * 3] - p[a * 3],
    vy = p[c * 3 + 1] - p[a * 3 + 1],
    vz = p[c * 3 + 2] - p[a * 3 + 2];
  return [uy * vz - uz * vy, uz * vx - ux * vz, ux * vy - uy * vx];
}