websperiments/standalone/goldberg-sphere/topology.mjs

// 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;
}

// ---------------------------------------------------------------------------
// Face adjacency (shared-edge graph)
// ---------------------------------------------------------------------------

// Build adjacency list for Goldberg faces via shared edges.
// Uses unrelaxed vertex positions as canonical keys — topology is stable across stages.
export function build_goldberg_adjacency(goldberg) {
	const edge_to_entries = new Map();
	for (let fi = 0; fi < goldberg.faces.length; fi++) {
		const verts = goldberg.faces[fi].vertices_3d;
		const n = verts.length;
		for (let i = 0; i < n; i++) {
			const a = verts[i], b = verts[(i + 1) % n];
			const ka = `${a.x.toFixed(5)},${a.y.toFixed(5)},${a.z.toFixed(5)}`;
			const kb = `${b.x.toFixed(5)},${b.y.toFixed(5)},${b.z.toFixed(5)}`;
			const key = ka < kb ? `${ka}|${kb}` : `${kb}|${ka}`;
			if (!edge_to_entries.has(key)) { edge_to_entries.set(key, []); }
			edge_to_entries.get(key).push({ fi, edge_idx: i });
		}
	}
	const adj = Array.from({ length: goldberg.faces.length }, () => []);
	for (const entries of edge_to_entries.values()) {
		if (entries.length === 2) {
			const [e0, e1] = entries;
			adj[e0.fi].push({ fj: e1.fi, fi_edge_idx: e0.edge_idx });
			adj[e1.fi].push({ fj: e0.fi, fi_edge_idx: e1.edge_idx });
		}
	}
	return adj;
}