qsm_core/bet/
icosphere.rs

1//! Icosphere mesh generation
2
3use std::collections::HashMap;
4
5/// Create a tessellated icosphere mesh
6///
7/// # Arguments
8/// * `subdivisions` - Number of subdivision levels (4 gives 2562 vertices)
9///
10/// # Returns
11/// Tuple of (vertices, faces) where vertices is Vec<[f64; 3]> and faces is Vec<[usize; 3]>
12pub fn create_icosphere(subdivisions: usize) -> (Vec<[f64; 3]>, Vec<[usize; 3]>) {
13    let phi = (1.0 + 5.0_f64.sqrt()) / 2.0;
14
15    // Initial icosahedron vertices
16    let mut vertices: Vec<[f64; 3]> = vec![
17        [-1.0,  phi, 0.0], [ 1.0,  phi, 0.0], [-1.0, -phi, 0.0], [ 1.0, -phi, 0.0],
18        [ 0.0, -1.0,  phi], [ 0.0,  1.0,  phi], [ 0.0, -1.0, -phi], [ 0.0,  1.0, -phi],
19        [ phi, 0.0, -1.0], [ phi, 0.0,  1.0], [-phi, 0.0, -1.0], [-phi, 0.0,  1.0],
20    ];
21
22    // Normalize to unit sphere
23    let norm = (vertices[0][0].powi(2) + vertices[0][1].powi(2) + vertices[0][2].powi(2)).sqrt();
24    for v in vertices.iter_mut() {
25        v[0] /= norm;
26        v[1] /= norm;
27        v[2] /= norm;
28    }
29
30    // Initial icosahedron faces
31    let mut faces: Vec<[usize; 3]> = vec![
32        [0, 11, 5], [0, 5, 1], [0, 1, 7], [0, 7, 10], [0, 10, 11],
33        [1, 5, 9], [5, 11, 4], [11, 10, 2], [10, 7, 6], [7, 1, 8],
34        [3, 9, 4], [3, 4, 2], [3, 2, 6], [3, 6, 8], [3, 8, 9],
35        [4, 9, 5], [2, 4, 11], [6, 2, 10], [8, 6, 7], [9, 8, 1],
36    ];
37
38    // Subdivide
39    for _ in 0..subdivisions {
40        let (new_vertices, new_faces) = subdivide_icosphere(&vertices, &faces);
41        vertices = new_vertices;
42        faces = new_faces;
43    }
44
45    (vertices, faces)
46}
47
48/// Subdivide each triangle into 4 triangles
49fn subdivide_icosphere(vertices: &[[f64; 3]], faces: &[[usize; 3]]) -> (Vec<[f64; 3]>, Vec<[usize; 3]>) {
50    let mut new_vertices: Vec<[f64; 3]> = vertices.to_vec();
51    let mut edge_midpoints: HashMap<(usize, usize), usize> = HashMap::new();
52    let mut new_faces: Vec<[usize; 3]> = Vec::with_capacity(faces.len() * 4);
53
54    let mut get_midpoint = |i1: usize, i2: usize, verts: &mut Vec<[f64; 3]>| -> usize {
55        let key = if i1 < i2 { (i1, i2) } else { (i2, i1) };
56        if let Some(&idx) = edge_midpoints.get(&key) {
57            return idx;
58        }
59
60        let v1 = verts[i1];
61        let v2 = verts[i2];
62        let mut mid = [
63            (v1[0] + v2[0]) / 2.0,
64            (v1[1] + v2[1]) / 2.0,
65            (v1[2] + v2[2]) / 2.0,
66        ];
67
68        // Normalize to unit sphere
69        let norm = (mid[0].powi(2) + mid[1].powi(2) + mid[2].powi(2)).sqrt();
70        mid[0] /= norm;
71        mid[1] /= norm;
72        mid[2] /= norm;
73
74        let idx = verts.len();
75        verts.push(mid);
76        edge_midpoints.insert(key, idx);
77        idx
78    };
79
80    for &[v0, v1, v2] in faces {
81        let m01 = get_midpoint(v0, v1, &mut new_vertices);
82        let m12 = get_midpoint(v1, v2, &mut new_vertices);
83        let m20 = get_midpoint(v2, v0, &mut new_vertices);
84
85        new_faces.push([v0, m01, m20]);
86        new_faces.push([v1, m12, m01]);
87        new_faces.push([v2, m20, m12]);
88        new_faces.push([m01, m12, m20]);
89    }
90
91    (new_vertices, new_faces)
92}
93
94#[cfg(test)]
95mod tests {
96    use super::*;
97
98    #[test]
99    fn test_icosphere_subdivisions() {
100        // subdivisions=0: 12 vertices, 20 faces
101        let (v0, f0) = create_icosphere(0);
102        assert_eq!(v0.len(), 12);
103        assert_eq!(f0.len(), 20);
104
105        // subdivisions=1: 42 vertices, 80 faces
106        let (v1, f1) = create_icosphere(1);
107        assert_eq!(v1.len(), 42);
108        assert_eq!(f1.len(), 80);
109
110        // subdivisions=4: 2562 vertices
111        let (v4, f4) = create_icosphere(4);
112        assert_eq!(v4.len(), 2562);
113        assert_eq!(f4.len(), 5120);
114    }
115
116    #[test]
117    fn test_vertices_on_unit_sphere() {
118        let (vertices, _) = create_icosphere(2);
119        for v in vertices {
120            let norm = (v[0].powi(2) + v[1].powi(2) + v[2].powi(2)).sqrt();
121            assert!((norm - 1.0).abs() < 1e-10, "Vertex not on unit sphere: norm = {}", norm);
122        }
123    }
124}
qsm_core/bet/icosphere.rs

qsm_core/bet/
icosphere.rs
