| """ |
| Mesh geodetica sferica a tasselli triangolari. |
| |
| Parte da un icosaedro regolare e lo suddivide con frequenza `frequency` |
| (ogni faccia → frequency² triangoli). Espone: |
| - indici faccia → 3 vertici |
| - tabella di adiacenza: 3 vicini per faccia (ordine lati 0,1,2) |
| - packing a coppie di facce adiacenti per GEMM rettangolari |
| """ |
|
|
| from __future__ import annotations |
|
|
| from dataclasses import dataclass |
| from typing import Dict, List, Optional, Tuple |
|
|
| import numpy as np |
|
|
| _PHI = (1.0 + np.sqrt(5.0)) / 2.0 |
|
|
|
|
| def _normalize(v: np.ndarray) -> np.ndarray: |
| n = np.linalg.norm(v, axis=-1, keepdims=True) |
| return v / np.maximum(n, 1e-12) |
|
|
|
|
| def icosahedron_vertices() -> np.ndarray: |
| """12 vertici icosaedro (ordine classico three.js / Wikipedia).""" |
| verts = np.asarray( |
| [ |
| [-1, _PHI, 0], |
| [1, _PHI, 0], |
| [-1, -_PHI, 0], |
| [1, -_PHI, 0], |
| [0, -1, _PHI], |
| [0, 1, _PHI], |
| [0, -1, -_PHI], |
| [0, 1, -_PHI], |
| [_PHI, 0, -1], |
| [_PHI, 0, 1], |
| [-_PHI, 0, -1], |
| [-_PHI, 0, 1], |
| ], |
| dtype=np.float64, |
| ) |
| return _normalize(verts) |
|
|
|
|
| def icosahedron_faces() -> np.ndarray: |
| """20 facce triangolari.""" |
| return np.asarray( |
| [ |
| [0, 11, 5], |
| [0, 5, 1], |
| [0, 1, 7], |
| [0, 7, 10], |
| [0, 10, 11], |
| [1, 5, 9], |
| [5, 11, 4], |
| [11, 10, 2], |
| [10, 7, 6], |
| [7, 1, 8], |
| [3, 9, 4], |
| [3, 4, 2], |
| [3, 2, 6], |
| [3, 6, 8], |
| [3, 8, 9], |
| [4, 9, 5], |
| [2, 4, 11], |
| [6, 2, 10], |
| [8, 6, 7], |
| [9, 8, 1], |
| ], |
| dtype=np.int64, |
| ) |
|
|
|
|
| def _edge_key(a: int, b: int) -> Tuple[int, int]: |
| return (a, b) if a < b else (b, a) |
|
|
|
|
| @dataclass(frozen=True) |
| class GeodesicMesh: |
| frequency: int |
| vertices: np.ndarray |
| faces: np.ndarray |
| adjacency: np.ndarray |
| face_centroids: np.ndarray |
| pairs: np.ndarray |
| unpaired: np.ndarray |
| compass_spiral: bool = False |
| spiral_order: Optional[np.ndarray] = None |
| spiral_rank: Optional[np.ndarray] = None |
| spiral_band: Optional[np.ndarray] = None |
|
|
| @property |
| def num_faces(self) -> int: |
| return int(self.faces.shape[0]) |
|
|
| @property |
| def num_vertices(self) -> int: |
| return int(self.vertices.shape[0]) |
|
|
|
|
| def build_geodesic_mesh( |
| frequency: int = 2, |
| *, |
| compass_spiral: bool = False, |
| ) -> GeodesicMesh: |
| """ |
| Mesh geodetica a frequenza `frequency` (>= 1). |
| Faces = 20 * frequency². |
| Se `compass_spiral=True`, riordina i 3 vicini esistenti di ogni faccia |
| lungo la spirale da polo a polo: |
| porta 0 = avanti lungo la spirale, |
| porta 1 = indietro lungo la spirale, |
| porta 2 = cambio di fascia (nord/sud). |
| """ |
| if frequency < 1: |
| raise ValueError("frequency deve essere >= 1") |
|
|
| base_v = icosahedron_vertices() |
| base_f = icosahedron_faces() |
|
|
| vert_list: List[np.ndarray] = [base_v[i].copy() for i in range(len(base_v))] |
| |
| edge_cache: Dict[Tuple[int, int, int], int] = {} |
|
|
| def point_on_edge(i: int, j: int, steps_from_i: int) -> int: |
| """steps_from_i=0 → i, =frequency → j.""" |
| if steps_from_i == 0: |
| return i |
| if steps_from_i == frequency: |
| return j |
| lo, hi = _edge_key(i, j) |
| steps_from_lo = steps_from_i if i == lo else frequency - steps_from_i |
| key = (lo, hi, steps_from_lo) |
| if key in edge_cache: |
| return edge_cache[key] |
| t = steps_from_lo / float(frequency) |
| p = _normalize((1.0 - t) * base_v[lo] + t * base_v[hi]) |
| idx = len(vert_list) |
| vert_list.append(p) |
| edge_cache[key] = idx |
| return idx |
|
|
| faces: List[Tuple[int, int, int]] = [] |
|
|
| for v0, v1, v2 in base_f: |
| v0, v1, v2 = int(v0), int(v1), int(v2) |
| grid: Dict[Tuple[int, int], int] = {} |
|
|
| for a in range(frequency + 1): |
| for b in range(frequency + 1 - a): |
| c = frequency - a - b |
| |
| if a == 0 and b == 0: |
| grid[(a, b)] = v0 |
| elif b == 0: |
| grid[(a, b)] = point_on_edge(v0, v1, a) |
| elif a == 0: |
| grid[(a, b)] = point_on_edge(v0, v2, b) |
| elif c == 0: |
| grid[(a, b)] = point_on_edge(v1, v2, b) |
| else: |
| p = _normalize( |
| (c * base_v[v0] + a * base_v[v1] + b * base_v[v2]) / float(frequency) |
| ) |
| idx = len(vert_list) |
| vert_list.append(p) |
| grid[(a, b)] = idx |
|
|
| for a in range(frequency): |
| for b in range(frequency - a): |
| i00 = grid[(a, b)] |
| i10 = grid[(a + 1, b)] |
| i01 = grid[(a, b + 1)] |
| faces.append((i00, i10, i01)) |
| if a + b + 1 < frequency: |
| i11 = grid[(a + 1, b + 1)] |
| faces.append((i10, i11, i01)) |
|
|
| vertices = _normalize(np.stack(vert_list, axis=0)) |
| faces_arr = np.asarray(faces, dtype=np.int64) |
| expected = 20 * frequency * frequency |
| if faces_arr.shape[0] != expected: |
| raise RuntimeError(f"Attese {expected} facce, ottenute {faces_arr.shape[0]}") |
|
|
| adjacency = build_adjacency(faces_arr) |
| centroids = _normalize(vertices[faces_arr].mean(axis=1)) |
| spiral_adj, spiral_order, spiral_rank, spiral_band = reorder_adjacency_spiral_compass( |
| adjacency, centroids |
| ) |
| if compass_spiral: |
| adjacency = spiral_adj |
| pairs, unpaired = pack_face_pairs(adjacency) |
|
|
| return GeodesicMesh( |
| frequency=frequency, |
| vertices=vertices, |
| faces=faces_arr, |
| adjacency=adjacency, |
| face_centroids=centroids, |
| pairs=pairs, |
| unpaired=unpaired, |
| compass_spiral=bool(compass_spiral), |
| spiral_order=spiral_order, |
| spiral_rank=spiral_rank, |
| spiral_band=spiral_band, |
| ) |
|
|
|
|
| def reorder_adjacency_spiral_compass( |
| adjacency: np.ndarray, |
| face_centroids: np.ndarray, |
| ) -> Tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]: |
| """ |
| Ordina i centroidi lungo una spirale da polo a polo attorno all'asse polare |
| icosaedrico (vertice 0 → vertice 3, simmetria 5-fold) e permuta i 3 vicini |
| già presenti in `adjacency` per ogni faccia `f`: |
| - porta 0: vicino «avanti» lungo la spirale (+Δlon sulla fascia), |
| - porta 1: vicino «indietro» lungo la spirale (-Δlon sulla fascia), |
| - porta 2: cambio di fascia (il terzo vicino, attraversamento nord/sud). |
| |
| Non aggiunge né rimuove archi o facce: per ogni f, set(new_adj[f]) == set(adjacency[f]). |
| Restituisce (new_adj, spiral_order, spiral_rank, spiral_band). |
| """ |
| n_faces = int(face_centroids.shape[0]) |
| if adjacency.shape != (n_faces, 3): |
| raise ValueError(f"adjacency attesa ({n_faces}, 3), ricevuta {adjacency.shape}") |
|
|
| |
| pole = icosahedron_vertices()[0] |
| u3 = pole / np.linalg.norm(pole) |
| ref = np.array([0.0, 0.0, 1.0], dtype=np.float64) |
| u1 = ref - np.dot(ref, u3) * u3 |
| u1 /= np.linalg.norm(u1) |
| u2 = np.cross(u3, u1) |
|
|
| z = face_centroids @ u3 |
| lon = np.mod(np.arctan2(face_centroids @ u2, face_centroids @ u1), 2.0 * np.pi) |
|
|
| new_adj = np.empty_like(adjacency) |
| for f in range(n_faces): |
| nbs = adjacency[f] |
| |
| dlons = [ |
| float(((lon[int(nb)] - lon[f] + np.pi) % (2.0 * np.pi)) - np.pi) |
| for nb in nbs |
| ] |
| |
| |
| idx_band = int(np.argmin(np.abs(dlons))) |
| rem = [i for i in range(3) if i != idx_band] |
| if dlons[rem[0]] >= dlons[rem[1]]: |
| idx_fwd, idx_bwd = rem[0], rem[1] |
| else: |
| idx_fwd, idx_bwd = rem[1], rem[0] |
| new_adj[f, 0] = nbs[idx_fwd] |
| new_adj[f, 1] = nbs[idx_bwd] |
| new_adj[f, 2] = nbs[idx_band] |
|
|
| |
| visited = np.zeros(n_faces, dtype=bool) |
| bands: List[List[int]] = [] |
| for f in range(n_faces): |
| if visited[f]: |
| continue |
| cyc: List[int] = [] |
| cur = f |
| while not visited[cur]: |
| visited[cur] = True |
| cyc.append(cur) |
| cur = int(new_adj[cur, 0]) |
| |
| start_idx = int(np.argmin([lon[idx] for idx in cyc])) |
| cyc = cyc[start_idx:] + cyc[:start_idx] |
| bands.append(cyc) |
|
|
| |
| bands.sort(key=lambda cyc: -float(np.mean(z[cyc]))) |
|
|
| spiral_order_list: List[int] = [] |
| spiral_band = np.empty(n_faces, dtype=np.int64) |
| for b_idx, cyc in enumerate(bands): |
| for f_id in cyc: |
| spiral_band[f_id] = b_idx |
| spiral_order_list.append(f_id) |
|
|
| spiral_order = np.asarray(spiral_order_list, dtype=np.int64) |
| spiral_rank = np.empty(n_faces, dtype=np.int64) |
| spiral_rank[spiral_order] = np.arange(n_faces, dtype=np.int64) |
|
|
| return new_adj, spiral_order, spiral_rank, spiral_band |
|
|
|
|
| def build_adjacency(faces: np.ndarray) -> np.ndarray: |
| """adjacency[f, i] = faccia oltre lato i; lato0=(v0,v1), lato1=(v1,v2), lato2=(v2,v0).""" |
| fcount = faces.shape[0] |
| edge_to_face: Dict[Tuple[int, int], List[Tuple[int, int]]] = {} |
| for fi, (a, b, c) in enumerate(faces): |
| for li, (u, v) in enumerate(((int(a), int(b)), (int(b), int(c)), (int(c), int(a)))): |
| edge_to_face.setdefault(_edge_key(u, v), []).append((fi, li)) |
|
|
| adj = np.full((fcount, 3), -1, dtype=np.int64) |
| for occ in edge_to_face.values(): |
| if len(occ) != 2: |
| continue |
| (f0, l0), (f1, l1) = occ |
| adj[f0, l0] = f1 |
| adj[f1, l1] = f0 |
|
|
| if np.any(adj < 0): |
| bad = int(np.sum(adj < 0)) |
| raise RuntimeError(f"Adiacenza incompleta ({bad} lati aperti)") |
| return adj |
|
|
|
|
| def pack_face_pairs(adjacency: np.ndarray) -> Tuple[np.ndarray, np.ndarray]: |
| """Matching greedy di facce adiacenti → matrici 2×D contigue.""" |
| fcount = adjacency.shape[0] |
| used = np.zeros(fcount, dtype=bool) |
| pairs: List[Tuple[int, int]] = [] |
| for f in range(fcount): |
| if used[f]: |
| continue |
| partner = -1 |
| for n in adjacency[f]: |
| ni = int(n) |
| if not used[ni]: |
| partner = ni |
| break |
| if partner < 0: |
| continue |
| used[f] = True |
| used[partner] = True |
| pairs.append((f, partner) if f < partner else (partner, f)) |
| unpaired = np.asarray([i for i in range(fcount) if not used[i]], dtype=np.int64) |
| return np.asarray(pairs, dtype=np.int64), unpaired |
|
|
|
|
| def nearest_face(mesh: GeodesicMesh, directions: np.ndarray) -> np.ndarray: |
| """Argmax allineamento con centroidi faccia (prodotto scalare su sfera).""" |
| d = _normalize(np.asarray(directions, dtype=np.float64)) |
| scores = d @ mesh.face_centroids.T |
| return np.argmax(scores, axis=-1).astype(np.int64) |
|
|
|
|
| def describe_mesh(mesh: GeodesicMesh) -> str: |
| return ( |
| f"GeodesicMesh(freq={mesh.frequency}, V={mesh.num_vertices}, " |
| f"F={mesh.num_faces}, pairs={len(mesh.pairs)}, unpaired={len(mesh.unpaired)})" |
| ) |
|
|
