Smooths a mesh by Catmull-Clark subdivision, with sharp and semi-sharp creases, and finds where its vertices lie on the limit surface.
How the precision is kept
New positions are weighted means of old ones, computed in 64-bit and stored in a Float64Array. The per-edge sharpness that blends the smooth and the sharp rule is held in a Float32Array, and the crease layer it starts from is a 32-bit layer of the mesh. For a semi-sharp crease that 32-bit weight multiplies a 64-bit difference of positions.
Open item. Keep the per-level edge sharpness in a Float64Array, so a semi-sharp crease blends positions with a 64-bit weight.
Calculations under this policy
6 calculations in 1 source file. The code is read from the released app when this page loads. It is the calculation people run, not a copy kept for this page.
mesh/subdiv.jsModelCraft 0.8.1 32-bit in part6 calculationsReviewed
Tests tools/modelcraft-test/mesh.mjs
Face points
64-bit
The new point in the middle of each face.
F = Σ corners ÷ n
The corners are added in order round the face and divided once by their count. All values are 64-bit.
mesh/subdiv.js
lines 43 to 53
ModelCraft 0.8.1
// face points
const FP = new Float64Array(nf * 3);
for (let f = 0; f < nf; f++) {
if (fIdx[f] < 0) continue;
let x = 0, y = 0, z = 0, n = 0, h = m.fHe[f];
const h0 = h;
do { const v = hV[h] * 3; x += P[v]; y += P[v + 1]; z += P[v + 2]; n++; h = nx[h]; } while (h !== h0 && n < 1e6);
FP[f * 3] = x / n; FP[f * 3 + 1] = y / n; FP[f * 3 + 2] = z / n;
const o = (nv + ne + fIdx[f]) * 3;
out[o] = FP[f * 3]; out[o + 1] = FP[f * 3 + 1]; out[o + 2] = FP[f * 3 + 2];
}
Edge points
32-bit in part
The new point on each edge.
m = (a + b) ÷ 2; e = (a + b + F₁ + F₂) ÷ 4; point = e + (m − e) × s for 0 ≤ s < 1
The weight s is read from a Float32Array. With s = 0 the smooth point is used unchanged and with s ≥ 1 the midpoint, so plain and fully sharp edges are 64-bit throughout. For a semi-sharp edge (0 < s < 1) the offset (m − e) × s carries the weight's 32-bit rounding, about 7 digits of the offset.
smooth = (F̄ + 2 R̄ + (n − 3) V) ÷ n; crease = (6 V + a + b) ÷ 8; k = min(1, (sₐ + s_b) ÷ 2); point = smooth + (crease − smooth) × k
The smooth and crease positions are 64-bit means. The blend factor k is the mean of two sharpness values read from the Float32Array, so a vertex on a semi-sharp crease carries that 32-bit weight. A vertex with three or more sharp edges, or a border vertex with two edges, keeps its position exactly.
mesh/subdiv.js
lines 65 to 97
ModelCraft 0.8.1
// vertex points
for (let v = 0; v < nv; v++) {
const o = v * 3;
if (m.vertDead(v)) continue;
out[o] = P[o]; out[o + 1] = P[o + 1]; out[o + 2] = P[o + 2];
const ring = m.outgoing(v);
if (!ring.length) continue;
const border = tw[ring[ring.length - 1]] < 0;
// the edges round v: outgoing half-edges, and for a border vertex the incoming open edge too
const edges = ring.map((h) => ({ h, w: m.head(h), s: S[h] }));
if (border) { const p = m.hPv[ring[0]]; edges.push({ h: p, w: hV[p], s: S[p] }); }
const sharpE = edges.filter((e) => e.s > 0);
const n = edges.length;
// smooth position
let fx = 0, fy = 0, fz = 0;
for (const h of ring) { const f = hF[h] * 3; fx += FP[f]; fy += FP[f + 1]; fz += FP[f + 2]; }
fx /= ring.length; fy /= ring.length; fz /= ring.length;
let rx = 0, ry = 0, rz = 0;
for (const e of edges) { const w = e.w * 3; rx += (P[o] + P[w]) / 2; ry += (P[o + 1] + P[w + 1]) / 2; rz += (P[o + 2] + P[w + 2]) / 2; }
rx /= n; ry /= n; rz /= n;
const smooth = [(fx + 2 * rx + (n - 3) * P[o]) / n, (fy + 2 * ry + (n - 3) * P[o + 1]) / n, (fz + 2 * rz + (n - 3) * P[o + 2]) / n];
const hard = sharpE.filter((e) => e.s >= 1).length, soft = sharpE.length;
let res = smooth;
if (soft >= 3 || (border && n === 2)) res = [P[o], P[o + 1], P[o + 2]]; // a corner stays
else if (soft === 2) {
// along a crease or the border: the curve rule, blended in by the creases' average sharpness (borders count as sharp)
const [a, b] = sharpE, A = a.w * 3, B = b.w * 3;
const crease = [(6 * P[o] + P[A] + P[B]) / 8, (6 * P[o + 1] + P[A + 1] + P[B + 1]) / 8, (6 * P[o + 2] + P[A + 2] + P[B + 2]) / 8];
const k = Math.min(1, (Math.min(a.s, 1e9) + Math.min(b.s, 1e9)) / 2);
res = hard === 2 || k >= 1 ? crease : smooth.map((x, i) => x + (crease[i] - x) * k);
}
out[o] = res[0]; out[o + 1] = res[1]; out[o + 2] = res[2];
}
Edge sharpness per level
32-bit in part
The sharpness each new edge inherits.
s′ = max(0, s − 1); open edge: ∞; crease layer = min(1, s′ ÷ 4)
Sharpness starts as the crease layer times 4 (creaseOf) and goes down by one per level. It is stored in a Float32Array at every level, and the starting crease layer is itself 32-bit, so a crease typed as 0.3 enters as the nearest 32-bit value.
Interpolates per-corner values onto the new quads.
corner: own value; edge: mean of two; face: mean of the face's corners
The means are taken in 64-bit and written back to the layer's own array, which is 32-bit by the mesh's layer type. These layers hold texture coordinates, colors and weights, not positions.
mesh/subdiv.js
lines 129 to 138
ModelCraft 0.8.1
// corners: v_i = h's own value, e_i = mean(h, next), face = mean of the face, e_{i−1} = mean(prev, h)
const cs = m.faceCorners(f, []);
for (const L of cornerL) {
const A = r.attrArray(L.name), s = L.size, src = L.col.a;
for (let j = 0; j < s; j++) {
let fm = 0; for (const c of cs) fm += src[c * s + j]; fm /= cs.length;
const vals = [src[h * s + j], (src[h * s + j] + src[nxh * s + j]) / 2, fm, (src[p * s + j] + src[h * s + j]) / 2];
for (let c = 0; c < 4; c++) A[keep.cornerOf[k * 4 + c] * s + j] = vals[c];
}
}
Limit positions
64-bit
Where each vertex of an all-quad mesh lies on the limit surface.
smooth: (n² V + 4 Σ E + Σ F) ÷ (n (n + 5)); crease or border: (4 V + a + b) ÷ 6
64-bit sums over the ring, divided once. Here an edge counts as sharp only when its sharpness is 1 or more, so the 32-bit weight is compared, never multiplied in. Vertices with a face that is not a quad keep their position.
mesh/subdiv.js
lines 173 to 203
ModelCraft 0.8.1
/**
* Where each vertex of a subdivided (all-quad) mesh lies on the limit surface. Interior, smooth: (n²V + 4ΣE + ΣF) / (n(n + 5))
* with E its edge neighbours and F the far corners of its faces; on a crease or border: (4V + a + b) / 6; a corner stays.
*/
export function limitPositions(m) {
const P = m.P, out = new Float64Array(m.nv * 3), crease = m.attrArray('crease');
for (let v = 0; v < m.nv; v++) {
const o = v * 3;
out[o] = P[o]; out[o + 1] = P[o + 1]; out[o + 2] = P[o + 2];
if (m.vertDead(v) || m.vHe[v] < 0) continue;
const ring = m.outgoing(v), border = m.hTw[ring[ring.length - 1]] < 0;
const edges = ring.map((h) => ({ w: m.head(h), s: creaseOf(m, h, crease) }));
if (border) { const p = m.hPv[ring[0]]; edges.push({ w: m.hV[p], s: Infinity }); }
const sharp = edges.filter((e) => e.s >= 1);
if (sharp.length >= 3 || (border && edges.length === 2)) continue;
if (sharp.length === 2) {
const A = sharp[0].w * 3, B = sharp[1].w * 3;
for (let k = 0; k < 3; k++) out[o + k] = (4 * P[o + k] + P[A + k] + P[B + k]) / 6;
continue;
}
if (ring.some((h) => m.faceSize(m.hF[h]) !== 4)) continue; // only all-quad neighbourhoods (after one level every one is)
const n = ring.length;
let e = [0, 0, 0], f = [0, 0, 0];
for (const h of ring) {
const w = m.head(h) * 3, far = m.hV[m.hNx[m.hNx[h]]] * 3;
for (let k = 0; k < 3; k++) { e[k] += P[w + k]; f[k] += P[far + k]; }
}
for (let k = 0; k < 3; k++) out[o + k] = (n * n * P[o + k] + 4 * e[k] + f[k]) / (n * (n + 5));
}
return out;
}
In detail
Method
Catmull and Clark 1978 for the face, edge and vertex rules; crease rules after DeRose, Kass and Truong 1998; limit positions by the standard limit mask. The file's header cites all three.
Number types
Float64Array for positions; Float32Array for edge sharpness and for attribute layers (texture coordinates, colors, weights); Int32Array for topology.
Limits
At most 8 levels. Limit positions are worked out only where every face round a vertex is a quad.
Approach: 32-bit step, flagged
One step stores or exchanges 32-bit values that reach a number a person reads. It is listed here with its open item, and the step is either moved to 64-bit or the app will warn where the number is read.