// Sphere topology: icosahedron construction, subdivision, and Goldberg merging. // All data structures live here. No rendering, no DOM. import { Vec3, normalize_to_sphere, midpoint } from './geometry.mjs'; // --------------------------------------------------------------------------- // Low-level data classes // --------------------------------------------------------------------------- export class Vertex { constructor(index, x, y, z) { this.index = index; this.x = x; this.y = y; this.z = z; } to_vec3() { return new Vec3(this.x, this.y, this.z); } } export class Edge { constructor(index, vertex_a, vertex_b) { // vertex_a < vertex_b always (canonical ordering) this.index = index; this.vertex_a = vertex_a; this.vertex_b = vertex_b; this.face_left = null; // face whose winding traverses a→b this.face_right = null; // face whose winding traverses b→a } } export class Face { constructor(index, vertices) { this.index = index; this.vertices = vertices; // ordered vertex indices this.edges = []; // edge index for vertices[i]→vertices[(i+1)%n] this.edge_neighbors = []; // face index across edges[i]; null = boundary (shouldn't happen on closed manifold) this.edge_directions = []; // +1 if this face is on face_left of edge[i], -1 if face_right } get size() { return this.vertices.length; } } // --------------------------------------------------------------------------- // Polyhedron — triangulated surface, used at icosahedron and subdivision levels // --------------------------------------------------------------------------- export class Polyhedron { constructor(vertices, faces) { // vertices: array of Vertex // faces: array of Face (edges/neighbors populated by build_topology) this.vertices = vertices; this.faces = faces; this.edges = []; this._edge_map = new Map(); // 'a,b' → Edge index } // Build edge list and face connectivity from vertex/face data. // Call this once after setting this.vertices and this.faces with vertex lists. build_topology() { this.edges = []; this._edge_map = new Map(); for (const face of this.faces) { const n = face.vertices.length; face.edges = []; face.edge_directions = []; for (let i = 0; i < n; i++) { const va = face.vertices[i]; const vb = face.vertices[(i + 1) % n]; const key = va < vb ? `${va},${vb}` : `${vb},${va}`; if (!this._edge_map.has(key)) { const edge = new Edge(this.edges.length, Math.min(va, vb), Math.max(va, vb)); this._edge_map.set(key, edge.index); this.edges.push(edge); } const edge_index = this._edge_map.get(key); face.edges.push(edge_index); // Determine if this face is left or right of the edge. // Convention: if the face traverses a→b (i.e. va < vb in canonical form // and that matches the face's traversal direction), it's face_left. const edge = this.edges[edge_index]; if (va === edge.vertex_a) { // Face traverses a→b → face_left face.edge_directions.push(1); edge.face_left = face.index; } else { // Face traverses b→a → face_right face.edge_directions.push(-1); edge.face_right = face.index; } } } // Populate edge_neighbors on each face for (const face of this.faces) { face.edge_neighbors = []; for (let i = 0; i < face.edges.length; i++) { const edge = this.edges[face.edges[i]]; if (face.edge_directions[i] === 1) { // This face is face_left; neighbor is face_right face.edge_neighbors.push(edge.face_right); } else { // This face is face_right; neighbor is face_left face.edge_neighbors.push(edge.face_left); } } } } get_edge(va, vb) { const key = va < vb ? `${va},${vb}` : `${vb},${va}`; const idx = this._edge_map.get(key); return idx !== undefined ? this.edges[idx] : null; } // Validate basic topological invariants. Returns array of error strings. validate() { const errors = []; const v = this.vertices.length; const e = this.edges.length; const f = this.faces.length; if (v - e + f !== 2) { errors.push(`Euler characteristic: V(${v}) - E(${e}) + F(${f}) = ${v - e + f}, expected 2`); } for (const edge of this.edges) { if (edge.face_left === null || edge.face_right === null) { errors.push(`Edge ${edge.index} missing a face (left=${edge.face_left}, right=${edge.face_right})`); } } for (const face of this.faces) { for (let i = 0; i < face.edge_neighbors.length; i++) { if (face.edge_neighbors[i] === null || face.edge_neighbors[i] === undefined) { errors.push(`Face ${face.index} edge_neighbor[${i}] is null`); } } } return errors; } // Subdivide every triangle into 4 triangles by inserting edge midpoints. // Midpoints are projected onto the unit sphere. // Returns a new Polyhedron. subdivide() { const new_vertices = this.vertices.map(v => new Vertex(v.index, v.x, v.y, v.z)); const midpoint_map = new Map(); // edge_index → new vertex index // Create midpoint vertices for every edge for (const edge of this.edges) { const va = this.vertices[edge.vertex_a].to_vec3(); const vb = this.vertices[edge.vertex_b].to_vec3(); const mid = normalize_to_sphere(midpoint(va, vb)); const new_index = new_vertices.length; new_vertices.push(new Vertex(new_index, mid.x, mid.y, mid.z)); midpoint_map.set(edge.index, new_index); } // Split each triangle into 4 const new_face_defs = []; for (const face of this.faces) { const [a, b, c] = face.vertices; const edge_ab = this.get_edge(a, b).index; const edge_bc = this.get_edge(b, c).index; const edge_ca = this.get_edge(c, a).index; const mab = midpoint_map.get(edge_ab); const mbc = midpoint_map.get(edge_bc); const mca = midpoint_map.get(edge_ca); // Preserve winding: original face is (a, b, c) CCW from outside. new_face_defs.push([a, mab, mca]); new_face_defs.push([mab, b, mbc]); new_face_defs.push([mca, mbc, c]); new_face_defs.push([mab, mbc, mca]); // center triangle } const new_faces = new_face_defs.map((verts, i) => new Face(i, verts)); const poly = new Polyhedron(new_vertices, new_faces); poly.build_topology(); return poly; } // --------------------------------------------------------------------------- // Icosahedron constructor — pentagonal antiprism method // --------------------------------------------------------------------------- static build_icosahedron() { // Edge length = 1. All geometry derived from the equal-edge constraint. // // Pentagon circumradius: R = 1 / (2 * sin(π/5)) // Ring separation: Δz = sqrt(1 - 2R²(1 - cos(π/5))) // Cap height above/below ring: h = sqrt(1 - R²) // // Vertex layout: // 0 : top cap // 1..5 : top ring, angle 2πk/5 // 6..10 : bottom ring, angle 2πk/5 + π/5 (rotated 36°) // 11 : bottom cap const R = 1 / (2 * Math.sin(Math.PI / 5)); const delta_z = Math.sqrt(1 - 2 * R * R * (1 - Math.cos(Math.PI / 5))); const h = Math.sqrt(1 - R * R); // Place rings symmetrically about z=0 const z_top = delta_z / 2; const z_bot = -delta_z / 2; const raw_vertices = []; // Top cap raw_vertices.push([0, 0, z_top + h]); // Top ring for (let k = 0; k < 5; k++) { const angle = (2 * Math.PI * k) / 5; raw_vertices.push([R * Math.cos(angle), R * Math.sin(angle), z_top]); } // Bottom ring (rotated 36°) for (let k = 0; k < 5; k++) { const angle = (2 * Math.PI * k) / 5 + Math.PI / 5; raw_vertices.push([R * Math.cos(angle), R * Math.sin(angle), z_bot]); } // Bottom cap raw_vertices.push([0, 0, z_bot - h]); // Normalize all vertices to unit sphere const vertices = raw_vertices.map(([x, y, z], i) => { const v = normalize_to_sphere(new Vec3(x, y, z)); return new Vertex(i, v.x, v.y, v.z); }); // Build faces from antiprism structure. // All faces wound CCW when viewed from outside (outward normal). // // Top cap: vertex 0, top ring k and k+1 // Antiprism upward triangles: top[k], bot[k+1], bot[k] (CCW from outside) // Antiprism downward triangles: top[k], top[k+1], bot[k+1] (CCW from outside) // Bottom cap: bot[k], bot[k+1], vertex 11 const top = k => 1 + (k % 5); const bot = k => 6 + (k % 5); const TOP_CAP = 0; const BOT_CAP = 11; const face_defs = []; for (let k = 0; k < 5; k++) { // Top cap triangle face_defs.push([TOP_CAP, top(k), top(k + 1)]); // Antiprism "up" triangle: top(k) flanked by bot(k-1) and bot(k), both 36° away. // bot(k-1) = bot(k+4) using mod-5 arithmetic. face_defs.push([top(k), bot(k + 4), bot(k)]); // Antiprism "down" triangle: top(k), top(k+1) share bot(k) (36° from each). // Winding reversed so outward normal points away from origin. face_defs.push([top(k + 1), top(k), bot(k)]); // Bottom cap triangle face_defs.push([BOT_CAP, bot(k + 1), bot(k)]); } const faces = face_defs.map((verts, i) => new Face(i, verts)); const poly = new Polyhedron(vertices, faces); poly.build_topology(); return poly; } } // --------------------------------------------------------------------------- // Goldberg polyhedron — hexagons and pentagons // --------------------------------------------------------------------------- export class Goldberg_Face { constructor(index, center_vertex, component_triangles, vertices_3d) { this.index = index; this.center_vertex = center_vertex; // vertex index in the subdivided Polyhedron this.component_triangles = component_triangles; // triangle face indices (ordered ring) this.vertices_3d = vertices_3d; // Vec3 centroids of component triangles, in ring order this.relaxed_vertices_3d = null; // set by Goldberg_Polyhedron.relax_sphere() this.edge_neighbors = []; // Goldberg_Face index across edge i this.is_pentagon = false; } get size() { return this.component_triangles.length; } } export class Goldberg_Polyhedron { constructor(faces) { this.faces = faces; // array of Goldberg_Face } // Build from a subdivided triangulated Polyhedron. // Each vertex of the triangulation becomes one Goldberg face. static from_subdivided(poly) { // For each vertex, collect the ring of surrounding triangular faces. const vertex_star = new Array(poly.vertices.length).fill(null).map(() => []); for (const face of poly.faces) { for (const vi of face.vertices) { vertex_star[vi].push(face.index); } } const goldberg_faces = []; // Map from triangle face index → which Goldberg face it belongs to const triangle_to_goldberg = new Map(); for (let vi = 0; vi < poly.vertices.length; vi++) { const star = vertex_star[vi]; if (star.length < 3) { continue; } // Order the star faces into a ring by following adjacency. let ordered = order_face_ring(star, vi, poly); if (ordered === null) { console.warn(`Could not order face ring for vertex ${vi}`); continue; } // Compute centroids of the component triangles as Goldberg face vertices const vertices_3d = ordered.map(fi => { const face = poly.faces[fi]; const va = poly.vertices[face.vertices[0]].to_vec3(); const vb = poly.vertices[face.vertices[1]].to_vec3(); const vc = poly.vertices[face.vertices[2]].to_vec3(); return normalize_to_sphere(va.add(vb).add(vc).scale(1 / 3)); }); // Enforce consistent CCW winding when viewed from outside the sphere. // The outward direction at this vertex is its own position (unit sphere). // Compute the ring's normal as the sum of cross products of consecutive // centroid vectors; if it points inward (dot < 0 with center vertex), reverse. { const center = poly.vertices[vi].to_vec3(); const n = vertices_3d.length; let normal_x = 0, normal_y = 0, normal_z = 0; for (let i = 0; i < n; i++) { const a = vertices_3d[i]; const b = vertices_3d[(i + 1) % n]; normal_x += a.y * b.z - a.z * b.y; normal_y += a.z * b.x - a.x * b.z; normal_z += a.x * b.y - a.y * b.x; } const dot = normal_x * center.x + normal_y * center.y + normal_z * center.z; if (dot < 0) { ordered.reverse(); vertices_3d.reverse(); } } const gf = new Goldberg_Face(goldberg_faces.length, vi, ordered, vertices_3d); goldberg_faces.push(gf); for (const fi of ordered) { triangle_to_goldberg.set(fi, gf.index); } } // Mark pentagons: vertices with a star of 5 faces. // In a subdivided icosahedron, the 12 original vertices remain degree-5. for (const gf of goldberg_faces) { gf.is_pentagon = (gf.size === 5); } // Build edge_neighbors using the dual graph relationship. // // Goldberg edge i runs from centroid(ring[i]) to centroid(ring[i+1]). // By duality this edge separates the Goldberg face centered on `center_vertex` // from the Goldberg face centered on the vertex shared by ring[i] and ring[i+1] // that is NOT `center_vertex`. So the neighbor across edge i is simply: // vertex_to_goldberg[ shared_outer_vertex(ring[i], ring[i+1], center_vertex) ] // // This avoids any triangle→Goldberg mapping ambiguity entirely. const vertex_to_goldberg = new Map(); for (const gf of goldberg_faces) { vertex_to_goldberg.set(gf.center_vertex, gf.index); } for (const gf of goldberg_faces) { gf.edge_neighbors = new Array(gf.size).fill(null); for (let i = 0; i < gf.size; i++) { const tri_a = poly.faces[gf.component_triangles[i]]; const tri_b = poly.faces[gf.component_triangles[(i + 1) % gf.size]]; const shared_vi = find_shared_outer_vertex(tri_a, tri_b, gf.center_vertex); if (shared_vi !== null && vertex_to_goldberg.has(shared_vi)) { gf.edge_neighbors[i] = vertex_to_goldberg.get(shared_vi); } } } return new Goldberg_Polyhedron(goldberg_faces); } validate() { const errors = []; let pentagon_count = 0; for (const gf of this.faces) { if (gf.is_pentagon) { pentagon_count++; } for (let i = 0; i < gf.size; i++) { const nbr_idx = gf.edge_neighbors[i]; if (nbr_idx === null || nbr_idx === undefined) { errors.push(`Face ${gf.index} edge_neighbor[${i}] is null`); continue; } const nbr = this.faces[nbr_idx]; // Neighbor symmetry: nbr should contain gf.index in its edge_neighbors if (!nbr.edge_neighbors.includes(gf.index)) { errors.push(`Face ${gf.index} → ${nbr_idx}, but ${nbr_idx} does not list ${gf.index} as neighbor`); } } } if (pentagon_count !== 12) { errors.push(`Expected 12 pentagons, found ${pentagon_count}`); } return errors; } // BFS from a face, up to `depth` steps. // Returns array of { face, depth, parent_index, parent_edge_index } in BFS order. // The root entry has parent_index = null, parent_edge_index = null. get_neighborhood(face_index, depth) { const result = []; const visited = new Set(); const queue = [{ face_index, depth: 0, parent_index: null, parent_edge_index: null }]; visited.add(face_index); while (queue.length > 0) { const entry = queue.shift(); result.push({ face: this.faces[entry.face_index], depth: entry.depth, parent_index: entry.parent_index, parent_edge_index: entry.parent_edge_index, }); if (entry.depth >= depth) { continue; } const face = this.faces[entry.face_index]; for (let i = 0; i < face.edge_neighbors.length; i++) { const nbr_idx = face.edge_neighbors[i]; if (nbr_idx !== null && nbr_idx !== undefined && !visited.has(nbr_idx)) { visited.add(nbr_idx); queue.push({ face_index: nbr_idx, depth: entry.depth + 1, parent_index: entry.face_index, parent_edge_index: i, }); } } } return result; } // Relax all Goldberg vertices onto the unit sphere. // // Two independent spring types, each with its own alpha (set to 0 to disable): // alpha_edge — push each shared edge toward a uniform target length. // Scales well; recommended default ≈ 0.1. // alpha_centroid — push each vertex toward its face's average circumradius. // Adds shape regularity but converges poorly at high // subdivision depths; keep small (≤ 0.02) if used. // // Both springs are applied each iteration; vertices are re-normalised to the // unit sphere afterwards. Stores results as relaxed_vertices_3d on each face. relax_sphere(iters = 150, alpha_edge = 0, alpha_centroid = 0.04) { const faces = this.faces; const nf = faces.length; // --- Global union-find to identify shared vertices across faces --- const face_base = new Array(nf); let uf_size = 0; for (let i = 0; i < nf; i++) { face_base[i] = uf_size; uf_size += faces[i].size; } const uf_parent = Array.from({ length: uf_size }, (_, i) => i); const uf_find = (i) => { while (uf_parent[i] !== i) { uf_parent[i] = uf_parent[uf_parent[i]]; i = uf_parent[i]; } return i; }; const uf_union = (i, j) => { uf_parent[uf_find(i)] = uf_find(j); }; const adj_seen = new Set(); for (let fi_a = 0; fi_a < nf; fi_a++) { const face_a = faces[fi_a]; const na = face_a.size; const base_a = face_base[fi_a]; for (let ea = 0; ea < na; ea++) { const fi_b = face_a.edge_neighbors[ea]; if (fi_b === null || fi_b === undefined) { continue; } const pair_key = fi_a < fi_b ? fi_a * 65536 + fi_b : fi_b * 65536 + fi_a; if (adj_seen.has(pair_key)) { continue; } adj_seen.add(pair_key); const face_b = faces[fi_b]; const nb = face_b.size; const base_b = face_base[fi_b]; let eb = -1; for (let j = 0; j < nb; j++) { if (face_b.edge_neighbors[j] === fi_a) { eb = j; break; } } if (eb === -1) { continue; } uf_union(base_a + ea, base_b + (eb + 1) % nb); uf_union(base_a + (ea + 1) % na, base_b + eb); } } // Compact position IDs. const canon_to_pid = new Map(); let pid_count = 0; const face_pids = new Array(nf); for (let fi = 0; fi < nf; fi++) { const n = faces[fi].size; const base = face_base[fi]; const pids = new Array(n); for (let i = 0; i < n; i++) { const canon = uf_find(base + i); if (!canon_to_pid.has(canon)) { canon_to_pid.set(canon, pid_count++); } pids[i] = canon_to_pid.get(canon); } face_pids[fi] = pids; } // Initialise positions from vertices_3d (already on unit sphere). const px = new Float64Array(pid_count); const py = new Float64Array(pid_count); const pz = new Float64Array(pid_count); const written = new Uint8Array(pid_count); for (let fi = 0; fi < nf; fi++) { const face = faces[fi]; const pids = face_pids[fi]; for (let i = 0; i < face.size; i++) { const pid = pids[i]; if (!written[pid]) { const v = face.vertices_3d[i]; px[pid] = v.x; py[pid] = v.y; pz[pid] = v.z; written[pid] = 1; } } } // Collect unique edges (needed for edge springs). const edge_seen = new Set(); const eu = [], ev = []; if (alpha_edge > 0) { for (let fi = 0; fi < nf; fi++) { const pids = face_pids[fi]; const n = faces[fi].size; for (let i = 0; i < n; i++) { const u = pids[i], v = pids[(i + 1) % n]; const ek = u < v ? u * 65536 + v : v * 65536 + u; if (!edge_seen.has(ek)) { edge_seen.add(ek); eu.push(u); ev.push(v); } } } } const ne = eu.length; // Target edge length: average of all initial edge lengths. let target_len = 0; for (let e = 0; e < ne; e++) { const u = eu[e], v = ev[e]; const dx = px[v] - px[u], dy = py[v] - py[u], dz = pz[v] - pz[u]; target_len += Math.sqrt(dx*dx + dy*dy + dz*dz); } if (ne > 0) { target_len /= ne; } // Target circumradius: average centroid-to-vertex distance (initial geometry). let target_r = 0, tr_count = 0; if (alpha_centroid > 0) { for (let fi = 0; fi < nf; fi++) { const pids = face_pids[fi]; const n = faces[fi].size; let fcx = 0, fcy = 0, fcz = 0; for (let i = 0; i < n; i++) { fcx += px[pids[i]]; fcy += py[pids[i]]; fcz += pz[pids[i]]; } fcx /= n; fcy /= n; fcz /= n; for (let i = 0; i < n; i++) { const p = pids[i]; const dx = px[p] - fcx, dy = py[p] - fcy, dz = pz[p] - fcz; target_r += Math.sqrt(dx*dx + dy*dy + dz*dz); tr_count++; } } target_r /= tr_count; } // Spring relaxation loop. const dx_buf = new Float64Array(pid_count); const dy_buf = new Float64Array(pid_count); const dz_buf = new Float64Array(pid_count); for (let iter = 0; iter < iters; iter++) { dx_buf.fill(0); dy_buf.fill(0); dz_buf.fill(0); // Edge-length springs. for (let e = 0; e < ne; e++) { const u = eu[e], v = ev[e]; const dx = px[v] - px[u], dy = py[v] - py[u], dz = pz[v] - pz[u]; const len = Math.sqrt(dx*dx + dy*dy + dz*dz); if (len < 1e-10) { continue; } const err = (len - target_len) / len * 0.5 * alpha_edge; dx_buf[u] += dx * err; dy_buf[u] += dy * err; dz_buf[u] += dz * err; dx_buf[v] -= dx * err; dy_buf[v] -= dy * err; dz_buf[v] -= dz * err; } // Centroid-to-vertex springs. if (alpha_centroid > 0) { for (let fi = 0; fi < nf; fi++) { const pids = face_pids[fi]; const n = faces[fi].size; let fcx = 0, fcy = 0, fcz = 0; for (let i = 0; i < n; i++) { fcx += px[pids[i]]; fcy += py[pids[i]]; fcz += pz[pids[i]]; } fcx /= n; fcy /= n; fcz /= n; for (let i = 0; i < n; i++) { const p = pids[i]; const dx = px[p] - fcx, dy = py[p] - fcy, dz = pz[p] - fcz; const dist = Math.sqrt(dx*dx + dy*dy + dz*dz); if (dist < 1e-10) { continue; } const err = (dist - target_r) / dist * alpha_centroid; dx_buf[p] -= dx * err; dy_buf[p] -= dy * err; dz_buf[p] -= dz * err; } } } // Apply deltas and project back to unit sphere. for (let i = 0; i < pid_count; i++) { const nx = px[i] + dx_buf[i]; const ny = py[i] + dy_buf[i]; const nz = pz[i] + dz_buf[i]; const r = Math.sqrt(nx*nx + ny*ny + nz*nz); if (r < 1e-10) { continue; } px[i] = nx / r; py[i] = ny / r; pz[i] = nz / r; } } // Write relaxed_vertices_3d back onto each face. for (let fi = 0; fi < nf; fi++) { const face = faces[fi]; const pids = face_pids[fi]; face.relaxed_vertices_3d = pids.map(pid => new Vec3(px[pid], py[pid], pz[pid])); } } } // --------------------------------------------------------------------------- // Internal helpers // --------------------------------------------------------------------------- // Order the faces in `star` (all containing vertex `vi`) into a ring. // Traces the ring purely through shared outer vertices — no edge_neighbors dependency. // // Each face in the ring is a triangle [vi, va, vb]. Two consecutive ring faces share // exactly one outer vertex (the one between them). We build a map from each outer // vertex to the (at most 2) star faces that contain it, then follow the chain. // Returns ordered array of face indices, or null on failure. function order_face_ring(star, vi, poly) { if (star.length === 0) { return []; } if (star.length === 1) { return [star[0]]; } // Map: outer_vertex → [face_indices in this star that contain it] const vertex_to_star_faces = new Map(); for (const fi of star) { for (const v of poly.faces[fi].vertices) { if (v === vi) { continue; } if (!vertex_to_star_faces.has(v)) { vertex_to_star_faces.set(v, []); } vertex_to_star_faces.get(v).push(fi); } } // Trace the ring: from current face [vi, v_prev, v_next], the forward outer vertex // is v_next. The next ring face is the other star face that contains vi and v_next. const ordered = [star[0]]; const visited = new Set([star[0]]); const first_face = poly.faces[star[0]]; const first_outers = first_face.vertices.filter(v => v !== vi); // Arbitrary starting direction; winding will be corrected after ring is built. let forward_vertex = first_outers[1]; while (ordered.length < star.length) { const candidates = vertex_to_star_faces.get(forward_vertex) || []; const next_fi = candidates.find(fi => !visited.has(fi)); if (next_fi === undefined) { return null; } ordered.push(next_fi); visited.add(next_fi); // The next forward vertex is the outer vertex of next_fi that isn't forward_vertex. const next_outers = poly.faces[next_fi].vertices.filter(v => v !== vi); forward_vertex = next_outers.find(v => v !== forward_vertex); if (forward_vertex === undefined) { return null; } } return ordered; } // Return the vertex that appears in both face_a and face_b but is not exclude_vi. // This is the "shared outer vertex" between two consecutive ring triangles. function find_shared_outer_vertex(face_a, face_b, exclude_vi) { for (const va of face_a.vertices) { if (va !== exclude_vi && face_b.vertices.includes(va)) { return va; } } return null; }