qsm_core/bet/
icosphere.rs

1//! Icosphere mesh generation
2
3use std::collections::HashMap;
4
5#[cfg(test)]
6/// Create a tessellated icosphere mesh
7///
8/// # Arguments
9/// * `subdivisions` - Number of subdivision levels (4 gives 2562 vertices)
10///
11/// # Returns
12/// Tuple of (vertices, faces) where vertices is Vec<[f64; 3]> and faces is Vec<[usize; 3]>
13pub fn create_icosphere(subdivisions: usize) -> (Vec<[f64; 3]>, Vec<[usize; 3]>) {
14    let phi = (1.0 + 5.0_f64.sqrt()) / 2.0;
15
16    // Initial icosahedron vertices
17    let mut vertices: Vec<[f64; 3]> = vec![
18        [-1.0,  phi, 0.0], [ 1.0,  phi, 0.0], [-1.0, -phi, 0.0], [ 1.0, -phi, 0.0],
19        [ 0.0, -1.0,  phi], [ 0.0,  1.0,  phi], [ 0.0, -1.0, -phi], [ 0.0,  1.0, -phi],
20        [ phi, 0.0, -1.0], [ phi, 0.0,  1.0], [-phi, 0.0, -1.0], [-phi, 0.0,  1.0],
21    ];
22
23    // Normalize to unit sphere
24    let norm = (vertices[0][0].powi(2) + vertices[0][1].powi(2) + vertices[0][2].powi(2)).sqrt();
25    for v in vertices.iter_mut() {
26        v[0] /= norm;
27        v[1] /= norm;
28        v[2] /= norm;
29    }
30
31    // Initial icosahedron faces
32    let mut faces: Vec<[usize; 3]> = vec![
33        [0, 11, 5], [0, 5, 1], [0, 1, 7], [0, 7, 10], [0, 10, 11],
34        [1, 5, 9], [5, 11, 4], [11, 10, 2], [10, 7, 6], [7, 1, 8],
35        [3, 9, 4], [3, 4, 2], [3, 2, 6], [3, 6, 8], [3, 8, 9],
36        [4, 9, 5], [2, 4, 11], [6, 2, 10], [8, 6, 7], [9, 8, 1],
37    ];
38
39    // Subdivide
40    for _ in 0..subdivisions {
41        let (new_vertices, new_faces) = subdivide_icosphere(&vertices, &faces);
42        vertices = new_vertices;
43        faces = new_faces;
44    }
45
46    (vertices, faces)
47}
48
49#[cfg(test)]
50/// Subdivide each triangle into 4 triangles
51fn subdivide_icosphere(vertices: &[[f64; 3]], faces: &[[usize; 3]]) -> (Vec<[f64; 3]>, Vec<[usize; 3]>) {
52    let mut new_vertices: Vec<[f64; 3]> = vertices.to_vec();
53    let mut edge_midpoints: HashMap<(usize, usize), usize> = HashMap::new();
54    let mut new_faces: Vec<[usize; 3]> = Vec::with_capacity(faces.len() * 4);
55
56    let mut get_midpoint = |i1: usize, i2: usize, verts: &mut Vec<[f64; 3]>| -> usize {
57        let key = if i1 < i2 { (i1, i2) } else { (i2, i1) };
58        if let Some(&idx) = edge_midpoints.get(&key) {
59            return idx;
60        }
61
62        let v1 = verts[i1];
63        let v2 = verts[i2];
64        let mut mid = [
65            (v1[0] + v2[0]) / 2.0,
66            (v1[1] + v2[1]) / 2.0,
67            (v1[2] + v2[2]) / 2.0,
68        ];
69
70        // Normalize to unit sphere
71        let norm = (mid[0].powi(2) + mid[1].powi(2) + mid[2].powi(2)).sqrt();
72        mid[0] /= norm;
73        mid[1] /= norm;
74        mid[2] /= norm;
75
76        let idx = verts.len();
77        verts.push(mid);
78        edge_midpoints.insert(key, idx);
79        idx
80    };
81
82    for &[v0, v1, v2] in faces {
83        let m01 = get_midpoint(v0, v1, &mut new_vertices);
84        let m12 = get_midpoint(v1, v2, &mut new_vertices);
85        let m20 = get_midpoint(v2, v0, &mut new_vertices);
86
87        new_faces.push([v0, m01, m20]);
88        new_faces.push([v1, m12, m01]);
89        new_faces.push([v2, m20, m12]);
90        new_faces.push([m01, m12, m20]);
91    }
92
93    (new_vertices, new_faces)
94}
95
96/// FSL-BET2's `make_mesh_from_icosa(n, m)` (Bostrix/FSL-BET2 d5acd7f, `meshclass/mesh.cpp`):
97/// its icosahedron (vertex coordinates `tau = 0.8506508084`, `one = 0.5257311121`, in its
98/// orientation and order), every triangle's last two vertices swapped, then `n − 1` passes of
99/// `Mesh::retessellate` (each triangle `(v0, v1, v2)` split at its edge midpoints `p0 = (v1+v2)/2`,
100/// `p1 = (v0+v2)/2`, `p2 = (v0+v1)/2` into `(p2, p0, p1)`, `(p1, v0, p2)`, `(p0, v2, p1)`,
101/// `(p2, v1, p0)`, a midpoint shared with an earlier triangle reused) with every point scaled
102/// back onto the unit sphere by `1/norm` after each pass. `n = 5` (bet2's) gives 2562 vertices.
103pub fn bet2_icosa(n: usize) -> (Vec<[f64; 3]>, Vec<[usize; 3]>) {
104    const TAU: f64 = 0.8506508084;
105    const ONE: f64 = 0.5257311121;
106    // ZA ZB ZC ZD YA YB YC YD XA XB XC XD
107    let mut pts: Vec<[f64; 3]> = vec![
108        [TAU, ONE, 0.0], [-TAU, ONE, 0.0], [-TAU, -ONE, 0.0], [TAU, -ONE, 0.0],
109        [ONE, 0.0, TAU], [ONE, 0.0, -TAU], [-ONE, 0.0, -TAU], [-ONE, 0.0, TAU],
110        [0.0, TAU, ONE], [0.0, -TAU, ONE], [0.0, -TAU, -ONE], [0.0, TAU, -ONE],
111    ];
112    let (za, zb, zc, zd, ya, yb, yc, yd, xa, xb, xc, xd) = (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11);
113    let mut tris: Vec<[usize; 3]> = vec![
114        [ya, xa, yd], [ya, yd, xb], [yb, yc, xd], [yb, xc, yc], [za, ya, zd],
115        [za, zd, yb], [zc, yd, zb], [zc, zb, yc], [xa, za, xd], [xa, xd, zb],
116        [xb, xc, zd], [xb, zc, xc], [xa, ya, za], [xd, za, yb], [ya, xb, zd],
117        [yb, zd, xc], [yd, xa, zb], [yc, zb, xd], [yd, zc, xb], [yc, xc, zc],
118    ];
119    for t in tris.iter_mut() {
120        t.swap(1, 2); // Triangle::swap
121    }
122    for _ in 1..n {
123        let mut added: HashMap<[u64; 3], usize> = HashMap::new();
124        let mut new_tris = Vec::with_capacity(tris.len() * 4);
125        let mid = |a: [f64; 3], b: [f64; 3]| [(a[0] + b[0]) / 2.0, (a[1] + b[1]) / 2.0, (a[2] + b[2]) / 2.0];
126        for &[v0, v1, v2] in &tris {
127            let pt0 = mid(pts[v1], pts[v2]);
128            let pt2 = mid(pts[v0], pts[v1]);
129            let pt1 = mid(pts[v0], pts[v2]);
130            let mut get = |pt: [f64; 3], pts: &mut Vec<[f64; 3]>| -> usize {
131                // bet2 matches an earlier midpoint to 1e-8; the same edge gives the same bits.
132                let key = [pt[0].to_bits(), pt[1].to_bits(), pt[2].to_bits()];
133                *added.entry(key).or_insert_with(|| {
134                    pts.push(pt);
135                    pts.len() - 1
136                })
137            };
138            let p0 = get(pt0, &mut pts);
139            let p1 = get(pt1, &mut pts);
140            let p2 = get(pt2, &mut pts);
141            new_tris.push([p2, p0, p1]);
142            new_tris.push([p1, v0, p2]);
143            new_tris.push([p0, v2, p1]);
144            new_tris.push([p2, v1, p0]);
145        }
146        tris = new_tris;
147        for p in pts.iter_mut() {
148            let norm = (p[0] * p[0] + p[1] * p[1] + p[2] * p[2]).sqrt();
149            let inv = 1.0 / norm;
150            *p = [p[0] * inv, p[1] * inv, p[2] * inv];
151        }
152    }
153    (pts, tris)
154}
155
156#[cfg(test)]
157mod tests {
158    use super::*;
159
160    #[test]
161    fn test_icosphere_subdivisions() {
162        // subdivisions=0: 12 vertices, 20 faces
163        let (v0, f0) = create_icosphere(0);
164        assert_eq!(v0.len(), 12);
165        assert_eq!(f0.len(), 20);
166
167        // subdivisions=1: 42 vertices, 80 faces
168        let (v1, f1) = create_icosphere(1);
169        assert_eq!(v1.len(), 42);
170        assert_eq!(f1.len(), 80);
171
172        // subdivisions=4: 2562 vertices
173        let (v4, f4) = create_icosphere(4);
174        assert_eq!(v4.len(), 2562);
175        assert_eq!(f4.len(), 5120);
176    }
177
178    #[test]
179    fn test_vertices_on_unit_sphere() {
180        let (vertices, _) = create_icosphere(2);
181        for v in vertices {
182            let norm = (v[0].powi(2) + v[1].powi(2) + v[2].powi(2)).sqrt();
183            assert!((norm - 1.0).abs() < 1e-10, "Vertex not on unit sphere: norm = {}", norm);
184        }
185    }
186}
qsm_core/bet/icosphere.rs

qsm_core/bet/
icosphere.rs
