Flattens a piece of a mesh into the UV plane with as little distortion of angles as the method allows.
How the precision is kept
All arithmetic is 64-bit in Float64Arrays, and the answers come from iterations that stop at stated tolerances or iteration caps. The accuracy of the flattened shape is therefore set by those stops and by the conditioning of the island, not by the number type.
Calculations under this policy
8 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.
Finds the unknowns that minimize |Ax − b| for a sparse set of rows.
columns scaled by 1 ÷ |column|; stop when |Aᵀr|² ≤ tol²·max(1, |Aᵀr₀|²), tol = 1e-12
Conjugate gradients on the normal equations without forming them, in Float64Arrays; dot products are plain sums in index order. It stops at the tolerance, when q·q is not positive, or at the iteration cap, and returns the current iterate in each case.
mesh/uv.js
lines 106 to 137
ModelCraft 0.8.1
/** Rows: [{c: Int32Array cols, v: Float64Array vals, b}]; n unknowns; x0 start. Minimises |Ax − b|. */
function cgls(rows, n, x0, { iters = 4000, tol = 1e-12 } = {}) {
const colN = new Float64Array(n);
for (const r of rows) for (let k = 0; k < r.c.length; k++) colN[r.c[k]] += r.v[k] * r.v[k];
const D = new Float64Array(n);
for (let j = 0; j < n; j++) D[j] = colN[j] > 0 ? 1 / Math.sqrt(colN[j]) : 0;
const Av = (y, out) => { for (let i = 0; i < rows.length; i++) { const r = rows[i]; let s = 0; for (let k = 0; k < r.c.length; k++) s += r.v[k] * D[r.c[k]] * y[r.c[k]]; out[i] = s; } return out; };
const ATv = (z, out) => { out.fill(0); for (let i = 0; i < rows.length; i++) { const r = rows[i], zi = z[i]; if (!zi) continue; for (let k = 0; k < r.c.length; k++) out[r.c[k]] += r.v[k] * D[r.c[k]] * zi; } return out; };
const y = new Float64Array(n);
for (let j = 0; j < n; j++) y[j] = D[j] > 0 ? x0[j] / D[j] : 0;
const R = rows.length, r = new Float64Array(R), q = new Float64Array(R), s = new Float64Array(n), p = new Float64Array(n);
Av(y, r);
for (let i = 0; i < R; i++) r[i] = rows[i].b - r[i];
ATv(r, s); p.set(s);
let gamma = dot(s, s);
const g0 = gamma;
for (let it = 0; it < iters && gamma > tol * tol * Math.max(1, g0); it++) {
Av(p, q);
const qq = dot(q, q);
if (!(qq > 0)) break;
const a = gamma / qq;
for (let j = 0; j < n; j++) y[j] += a * p[j];
for (let i = 0; i < R; i++) r[i] -= a * q[i];
ATv(r, s);
const g2 = dot(s, s), b = g2 / gamma;
gamma = g2;
for (let j = 0; j < n; j++) p[j] = s[j] + b * p[j];
}
const x = new Float64Array(n);
for (let j = 0; j < n; j++) x[j] = y[j] * D[j];
return x;
}
Conformal map (LSCM)
64-bit
UV positions for the corners of an island from the 2D shape of each of its triangles, with two points pinned.
per triangle, with d = twice its signed area and s = 1 ÷ √|d|: Σⱼ Wⱼ·(uⱼ + i·vⱼ) = 0, W = s × the opposite edge vectors
Each triangle gives two rows scaled by 1 ÷ √|d|; a triangle with |d| not above 1e-30 is left out. Pinned points move to the right-hand side. The iteration cap passed to CGLS is 4 per free point, kept between 400 and 20 000.
mesh/uv.js
lines 140 to 180
ModelCraft 0.8.1
/**
* LSCM: wedge UVs from triangles (tri: wedge triples) with their own 2-D shapes (loc: 6 numbers per triangle). Two wedges
* pinned; `init` a starting guess (2 per wedge) the solver only improves.
*/
export function lscm(nW, tri, loc, pins, init) {
const pinned = new Map(pins.map(([w, u, v]) => [w, [u, v]]));
// unknowns: the free wedges' u and v
const col = new Int32Array(nW).fill(-1);
let n = 0;
for (let w = 0; w < nW; w++) if (!pinned.has(w)) col[w] = n++;
const rows = [];
for (let t = 0; t < tri.length / 3; t++) {
const [x0, y0, x1, y1, x2, y2] = loc.subarray(t * 6, t * 6 + 6);
const d = (x1 - x0) * (y2 - y0) - (x2 - x0) * (y1 - y0);
if (!(Math.abs(d) > 1e-30)) continue;
const s = 1 / Math.sqrt(Math.abs(d));
const W = [[(x2 - x1) * s, (y2 - y1) * s], [(x0 - x2) * s, (y0 - y2) * s], [(x1 - x0) * s, (y1 - y0) * s]];
// real part: Σ a u − b v; imaginary part: Σ b u + a v
for (const part of [0, 1]) {
const c = [], v = [];
let b = 0;
for (let j = 0; j < 3; j++) {
const w = tri[t * 3 + j], [a, bb] = W[j];
const cu = part ? bb : a, cv = part ? a : -bb;
const p = pinned.get(w);
if (p) { b -= cu * p[0] + cv * p[1]; continue; }
c.push(col[w] * 2, col[w] * 2 + 1); v.push(cu, cv);
}
if (c.length) rows.push({ c: Int32Array.from(c), v: Float64Array.from(v), b });
}
}
const x0 = new Float64Array(n * 2);
if (init) for (let w = 0; w < nW; w++) if (col[w] >= 0) { x0[col[w] * 2] = init[w * 2]; x0[col[w] * 2 + 1] = init[w * 2 + 1]; }
const x = cgls(rows, n * 2, x0, { iters: Math.max(400, Math.min(20000, n * 4)) });
const out = new Float64Array(nW * 2);
for (let w = 0; w < nW; w++) {
const p = pinned.get(w);
if (p) { out[w * 2] = p[0]; out[w * 2 + 1] = p[1]; } else { out[w * 2] = x[col[w] * 2]; out[w * 2 + 1] = x[col[w] * 2 + 1]; }
}
return out;
}
Angle based flattening
64-bit
Corner angles close to the surface's own that a flat triangulation can have.
minimize Σ (αᵢ − βᵢ)² ÷ βᵢ² subject to: each triangle sums to π, each inner wedge to 2π, Σ [ln sin α(next) − ln sin α(prev)] = 0 round each inner wedge
Target angles are clamped to 0.001…π − 0.001 and working angles to 0.0001…π − 0.0001 at every step, which keeps the sines and tangents away from zero. At most 30 Newton steps; it stops when the largest constraint residual is below 1e-9 or the largest step below 1e-12. The result is accepted only if the triangle and wedge sums hold within 1e-5; the wheel condition is not in that final test.
mesh/uv.js
lines 182 to 265
ModelCraft 0.8.1
/**
* ABF: triangle corner angles close to the surface's own (β, made to add up to 2π round inner wedges) that a flat
* triangulation can have: each triangle π, each inner wedge 2π, and each inner wedge's wheel closes (Π sin of the angles
* after = Π sin of the angles before). Newton steps; each solves the constraints' Schur complement by conjugate gradients.
* → Float64Array of angles (3 per triangle), or null when it did not settle.
*/
export function abfAngles(nW, tri, beta, inner) {
const T = tri.length / 3, n = T * 3;
// the corners of each inner wedge
const at = Array.from({ length: nW }, () => []);
for (let i = 0; i < n; i++) at[tri[i]].push(i);
const innerW = [];
for (let w = 0; w < nW; w++) if (inner[w] && at[w].length >= 3) innerW.push(w);
const b = Float64Array.from(beta);
for (const w of innerW) { let s = 0; for (const i of at[w]) s += b[i]; if (s > 0) for (const i of at[w]) b[i] *= TAU / s; }
const wt = new Float64Array(n);
for (let i = 0; i < n; i++) { b[i] = Math.min(Math.PI - 1e-3, Math.max(1e-3, b[i])); wt[i] = 1 / (b[i] * b[i]); }
const a = Float64Array.from(b);
const R = T + 2 * innerW.length;
let lam = new Float64Array(R);
const next = (i) => i - (i % 3) + ((i % 3) + 1) % 3, prev = (i) => i - (i % 3) + ((i % 3) + 2) % 3;
// the constraint rows as lists (rebuilt each step: the wheel rows' values change)
for (let iter = 0; iter < 30; iter++) {
const rows = [], C = new Float64Array(R);
for (let t = 0; t < T; t++) { rows.push({ c: [t * 3, t * 3 + 1, t * 3 + 2], v: [1, 1, 1] }); C[t] = a[t * 3] + a[t * 3 + 1] + a[t * 3 + 2] - Math.PI; }
innerW.forEach((w, k) => {
const cs = at[w];
rows.push({ c: cs.slice(), v: cs.map(() => 1) });
let s = 0; for (const i of cs) s += a[i];
C[T + k] = s - TAU;
});
innerW.forEach((w, k) => {
const c = [], v = [];
let s = 0;
for (const i of at[w]) {
const nI = next(i), pI = prev(i);
s += Math.log(Math.sin(a[nI])) - Math.log(Math.sin(a[pI]));
c.push(nI, pI); v.push(1 / Math.tan(a[nI]), -1 / Math.tan(a[pI]));
}
rows.push({ c, v });
C[T + innerW.length + k] = s;
});
let worst = 0; for (let r = 0; r < R; r++) worst = Math.max(worst, Math.abs(C[r]));
if (worst < 1e-9 && iter > 0) break;
const Dinv = new Float64Array(n);
for (let i = 0; i < n; i++) Dinv[i] = 1 / (2 * wt[i]);
// g + Jᵀλ
const gl = new Float64Array(n);
for (let i = 0; i < n; i++) gl[i] = 2 * wt[i] * (a[i] - b[i]);
rows.forEach((r, k) => { const l = lam[k]; if (l) for (let j = 0; j < r.c.length; j++) gl[r.c[j]] += r.v[j] * l; });
// rhs = C − J D gl
const rhs = new Float64Array(R);
rows.forEach((r, k) => { let s = 0; for (let j = 0; j < r.c.length; j++) s += r.v[j] * Dinv[r.c[j]] * gl[r.c[j]]; rhs[k] = C[k] - s; });
// S y = J D Jᵀ y, preconditioned by its diagonal
const tmp = new Float64Array(n);
const S = (y, out) => {
tmp.fill(0);
rows.forEach((r, k) => { const l = y[k]; if (l) for (let j = 0; j < r.c.length; j++) tmp[r.c[j]] += r.v[j] * l; });
for (let i = 0; i < n; i++) tmp[i] *= Dinv[i];
rows.forEach((r, k) => { let s = 0; for (let j = 0; j < r.c.length; j++) s += r.v[j] * tmp[r.c[j]]; out[k] = s; });
return out;
};
const diag = new Float64Array(R);
rows.forEach((r, k) => { let s = 0; for (let j = 0; j < r.c.length; j++) s += r.v[j] * r.v[j] * Dinv[r.c[j]]; diag[k] = s > 0 ? 1 / s : 0; });
const dl = pcg(S, rhs, diag, R);
// dα = −D (gl + Jᵀ dλ)
const jt = new Float64Array(n);
rows.forEach((r, k) => { const l = dl[k]; if (l) for (let j = 0; j < r.c.length; j++) jt[r.c[j]] += r.v[j] * l; });
let step = 0;
for (let i = 0; i < n; i++) {
const d = -Dinv[i] * (gl[i] + jt[i]);
a[i] = Math.min(Math.PI - 1e-4, Math.max(1e-4, a[i] + d));
step = Math.max(step, Math.abs(d));
}
for (let k = 0; k < R; k++) lam[k] += dl[k];
if (!a.every(Number.isFinite)) return null;
if (step < 1e-12) break;
}
// settled?
let worst = 0;
for (let t = 0; t < T; t++) worst = Math.max(worst, Math.abs(a[t * 3] + a[t * 3 + 1] + a[t * 3 + 2] - Math.PI));
for (const w of innerW) { let s = 0; for (const i of at[w]) s += a[i]; worst = Math.max(worst, Math.abs(s - TAU)); }
return worst < 1e-5 ? a : null;
}
Preconditioned conjugate gradients
64-bit
Solves the constraint system of each Newton step of the angle based method.
stop when r·r < 1e-14 × r₀·r₀
Diagonal preconditioning. The stop compares squared residuals, so it is reached when the residual is 1e-7 of its start; it also stops when p·q is not positive or after 2000 iterations. The starting r·r is floored at 1e-300.
mesh/uv.js
lines 266 to 285
ModelCraft 0.8.1
function pcg(A, b, Minv, n, { iters = 2000, tol = 1e-14 } = {}) {
const x = new Float64Array(n), r = Float64Array.from(b), z = new Float64Array(n), p = new Float64Array(n), q = new Float64Array(n);
for (let i = 0; i < n; i++) z[i] = r[i] * Minv[i];
p.set(z);
let rz = dot(r, z);
const r0 = Math.max(dot(r, r), 1e-300);
for (let it = 0; it < iters; it++) {
A(p, q);
const pq = dot(p, q);
if (!(pq > 0)) break;
const al = rz / pq;
for (let i = 0; i < n; i++) { x[i] += al * p[i]; r[i] -= al * q[i]; }
if (dot(r, r) < tol * r0) break;
for (let i = 0; i < n; i++) z[i] = r[i] * Minv[i];
const rz2 = dot(r, z), be = rz2 / rz;
rz = rz2;
for (let i = 0; i < n; i++) p[i] = z[i] + be * p[i];
}
return x;
}
Corner angles of a triangle
64-bit
The three corner angles of a triangle in 3D.
angle = acos(clamp(u·v ÷ (|u|·|v|), −1, 1))
The cosine is clamped to −1…1 before acos, which loses digits for angles near 0 and π. A corner with a side of zero length returns π ÷ 3.
mesh/uv.js
lines 293 to 297
ModelCraft 0.8.1
/** Corner angles of a 3-D triangle. */
function angles3(A, B, C) {
const ang = (p, q, r) => { const u = sub(q, p), v = sub(r, p), l = ln3(u) * ln3(v); return l > 0 ? Math.acos(Math.max(-1, Math.min(1, dt3(u, v) / l))) : Math.PI / 3; };
return [ang(A, B, C), ang(B, C, A), ang(C, A, B)];
}
A triangle's own 2D shape
64-bit
A triangle laid flat: A at the origin, B on the x axis, C where its side lengths put it.
B = (|AB|, 0); C = (AC·x, AC·y), x = AB ÷ |AB|, y = unit((AB × AC) × AB)
Two dot products with unit vectors. A zero length for AB or for the in-plane perpendicular is replaced by 1e-30 so the division does not fail.
mesh/uv.js
lines 298 to 303
ModelCraft 0.8.1
/** A triangle's own 2-D shape: A at 0, B on +x. */
function local2(A, B, C) {
const ab = sub(B, A), ac = sub(C, A), l = ln3(ab) || 1e-30, x = [ab[0] / l, ab[1] / l, ab[2] / l];
const n = crs(ab, ac), y = crs(n, ab), ly = ln3(y) || 1e-30;
return [0, 0, l, 0, dt3(ac, x), dt3(ac, [y[0] / ly, y[1] / ly, y[2] / ly])];
}
A triangle's shape from its angles
64-bit
A flat triangle with the given corner angles, sized by the law of sines.
A = (0, 0); B = (sin a₂, 0); C = (sin a₁·cos a₀, sin a₁·sin a₀)
Four trigonometric calls and two products. The triangle's size is set by the sines alone; lscm scales each triangle's rows by 1 ÷ √|d|, which removes the size.
mesh/uv.js
lines 304 to 308
ModelCraft 0.8.1
/** A triangle's shape from its angles (law of sines). */
function fromAngles(a0, a1, a2) {
const c = Math.sin(a2), b = Math.sin(a1);
return [0, 0, c, 0, b * Math.cos(a0), b * Math.sin(a0)];
}
Pinned points and first guess
64-bit
Chooses the two corners to pin and a starting layout for the solver.
pins: the two ends along the longest box axis at (0, 0) and (d, 0), d = their distance in 3D
The pin distance is the true 3D distance from Math.hypot (1 when it is zero), so an island unwrapped without packing is in mesh units. The first guess projects the corners on two unit directions; it only starts the iteration.
mesh/uv.js
lines 353 to 374
ModelCraft 0.8.1
function pinsAndGuess(nW, wp) {
// the longest axis of the box round the wedges, its two ends pinned at their true distance apart
let lo = [Infinity, Infinity, Infinity], hi = [-Infinity, -Infinity, -Infinity];
for (let w = 0; w < nW; w++) for (let k = 0; k < 3; k++) { lo[k] = Math.min(lo[k], wp[w * 3 + k]); hi[k] = Math.max(hi[k], wp[w * 3 + k]); }
const ax = [0, 1, 2].reduce((b, k) => (hi[k] - lo[k] > hi[b] - lo[b] ? k : b), 0);
let a = 0, b = 0;
for (let w = 0; w < nW; w++) { if (wp[w * 3 + ax] < wp[a * 3 + ax]) a = w; if (wp[w * 3 + ax] > wp[b * 3 + ax]) b = w; }
if (a === b) b = (a + 1) % Math.max(1, nW);
const d = Math.hypot(wp[b * 3] - wp[a * 3], wp[b * 3 + 1] - wp[a * 3 + 1], wp[b * 3 + 2] - wp[a * 3 + 2]) || 1;
// a first guess: the wedges projected on the plane that fits them best (Newell-like normal of the spread), placed like the pins
const c = [(lo[0] + hi[0]) / 2, (lo[1] + hi[1]) / 2, (lo[2] + hi[2]) / 2];
const e1 = norm3(sub([wp[b * 3], wp[b * 3 + 1], wp[b * 3 + 2]], [wp[a * 3], wp[a * 3 + 1], wp[a * 3 + 2]]));
// the direction of most spread square to e1
let best = null, bl = -1;
for (const k of [0, 1, 2]) { const t = [0, 0, 0]; t[k] = 1; const s = sub(t, e1.map((x) => x * dt3(t, e1))); const l = ln3(s); if (l > bl) { bl = l; best = s; } }
let e2 = norm3(best), spread = 0, spreadN = 0;
for (let w = 0; w < nW; w++) { const p = sub([wp[w * 3], wp[w * 3 + 1], wp[w * 3 + 2]], c); spread += Math.abs(dt3(p, e2)); spreadN += Math.abs(dt3(p, norm3(crs(e1, e2)))); }
if (spreadN > spread) e2 = norm3(crs(e1, e2));
const init = new Float64Array(nW * 2), A = [wp[a * 3], wp[a * 3 + 1], wp[a * 3 + 2]];
for (let w = 0; w < nW; w++) { const p = sub([wp[w * 3], wp[w * 3 + 1], wp[w * 3 + 2]], A); init[w * 2] = dt3(p, e1); init[w * 2 + 1] = dt3(p, e2); }
return { pins: [[a, 0, 0], [b, d, 0]], init };
}
In detail
Method
Least squares conformal maps (Lévy, Petitjean, Ray, Maillot, SIGGRAPH 2002) solved by CGLS with column scaling and two pinned points; angle based flattening (Sheffer and de Sturler 2001; ABF++, Sheffer, Lévy, Mogilnitsky, Bogomyakov 2005) by Newton steps on the Schur complement solved with preconditioned conjugate gradients. References as cited in the file's header.
CGLS: gradient norm below 1e-12 × max(1, start), at most 400…20 000 iterations by size. ABF: at most 30 Newton steps, stop at a constraint residual below 1e-9 or a step below 1e-12, accepted at 1e-5 on the triangle and wedge angle sums. PCG: squared residual below 1e-14 of its start, at most 2000 iterations. LSCM skips triangles with |2 × area| ≤ 1e-30.
Limits
If the iteration cap is reached the result is the last iterate, with no warning. When ABF does not settle the island is unwrapped from its own triangle shapes instead, and the note 'angle_fell' is returned.
Approach: 64-bit, stated limit
The arithmetic is 64-bit, but the method or the data sets the accuracy: a solver tolerance, a mesh, a step size, a datum shift, a table. The limit is stated where the result is shown.