subsphere/
hex.rs

1//! Contains types related to [`HexSphere`].
2use crate::prelude::*;
3use crate::proj::BaseTriProjector;
4use crate::tri::{self, BaseRegion, BaseRegionType, TriSphere};
5use std::num::NonZero;
6
7/// A tessellation of the unit sphere into *mostly* hexagonal [`Face`]s, constructed by grouping
8/// triangles of a [`TriSphere`].
9///
10/// Equivalently, this is a tessellation formed by projecting a
11/// [Goldberg polyhedron](https://en.wikipedia.org/wiki/Goldberg_polyhedron) onto
12/// the sphere.
13#[derive(Clone, Copy, Debug, PartialEq, Eq, Hash)]
14pub struct HexSphere<Proj> {
15    kis: TriSphere<Proj>,
16}
17
18impl<Proj: Clone + BaseTriProjector> TriSphere<Proj> {
19    /// [Truncates](https://en.wikipedia.org/wiki/Truncation_(geometry)) the vertices of this
20    /// [`TriSphere`] to create a [`HexSphere`].
21    ///
22    /// This will convert all existing faces into hexagons, and create a new face at each vertex.
23    pub fn truncate(self) -> HexSphere<Proj> {
24        HexSphere {
25            kis: self.subdivide_edge(NonZero::new(3).unwrap()),
26        }
27    }
28}
29
30impl<Proj> HexSphere<Proj> {
31    /// Attempts to construct a [`HexSphere`] whose [`kis`](HexSphere::kis) is the given
32    /// [`TriSphere`].
33    ///
34    /// For this to succeed, the [`TriSphere`] parameters must satisfy `b % 3 == c % 3`.
35    /// This will return [`None`] otherwise.
36    pub fn new(kis: TriSphere<Proj>) -> Option<Self> {
37        if (kis.b() + 2 * kis.c()) % 3 == 0 {
38            Some(Self { kis })
39        } else {
40            None
41        }
42    }
43
44    /// Performs the
45    /// [kis](https://en.wikipedia.org/wiki/Conway_polyhedron_notation#Original_operations)
46    /// operation on this [`HexSphere`].
47    ///
48    /// This creates a [`Face`](crate::Face) for each [`HalfEdge`] in the [`TriSphere`].
49    pub fn kis(self) -> TriSphere<Proj> {
50        self.kis
51    }
52
53    /// [`TriSphere::num_divisions`] for the dual of this [`HexSphere`].
54    const fn dual_num_divisions(&self) -> usize {
55        let b = self.kis.b() as usize;
56        let c = self.kis.c() as usize;
57        let dual_c = (b - c) / 3;
58        let dual_b = c + dual_c;
59        dual_b * (dual_b + dual_c) + dual_c * dual_c
60    }
61
62    /// The number of faces per edge [`BaseRegion`].
63    fn num_faces_per_edge_region(&self) -> usize {
64        let b = self.kis.b() as usize;
65        let c = self.kis.c() as usize;
66        num_faces_on_edge(b, c + 1)
67    }
68
69    /// The number of faces per interior [`BaseRegion`].
70    fn num_faces_per_interior_region(&self) -> usize {
71        let b = self.kis.b() as usize;
72        let c = self.kis.c() as usize;
73        let n = b - c;
74        if n == 0 {
75            (c % 3 == 0) as usize
76        } else {
77            num_faces_on_interior(c, n, n)
78        }
79    }
80
81    /// The number of vertices per edge [`BaseRegion`], assuming that the origin vertex is not
82    /// counted.
83    fn num_vertices_per_edge_region(&self) -> usize {
84        self.kis.num_vertices_per_edge_region() - self.num_faces_per_edge_region()
85    }
86
87    /// The number of vertices per interior [`BaseRegion`].
88    fn num_vertices_per_interior_region(&self) -> usize {
89        self.kis.num_vertices_per_interior_region() - self.num_faces_per_interior_region()
90    }
91}
92
93/// Counts the number of points `(u, v)` such that `0 < u < b`, `v < d` and `u % 3 == v % 3`.
94fn num_faces_on_edge(b: usize, d: usize) -> usize {
95    ((b - 1) * d + (b % 3) / 2) / 3
96}
97
98#[test]
99fn test_num_faces_on_edge() {
100    for b in 1..10 {
101        for d in 0..10 {
102            let mut count = 0;
103            for u in 1..b {
104                for v in 0..d {
105                    if u % 3 == v % 3 {
106                        count += 1;
107                    }
108                }
109            }
110            assert_eq!(
111                num_faces_on_edge(b, d),
112                count,
113                "failed for b = {}, d = {}",
114                b,
115                d
116            );
117        }
118    }
119}
120
121/// Counts the number of points `(u, v)` such that `0 < u + v < n`, `v < d` and
122/// `u % 3 == (v + k) % 3`. Assumes that `n` is a multiple of 3.
123fn num_faces_on_interior(k: usize, n: usize, d: usize) -> usize {
124    (d - 1) * (2 * n - d - 2) / 6 + ((k % 3 == 0) & (d % 3 != 1)) as usize
125}
126
127#[test]
128fn test_num_faces_on_interior() {
129    for k in 0..2 {
130        for n_3 in 1..10 {
131            let n = n_3 * 3;
132            for d in 1..n {
133                let mut count = 0;
134                for u in 1..n {
135                    for v in 1..d {
136                        if u + v < n && u % 3 == (v + k) % 3 {
137                            count += 1;
138                        }
139                    }
140                }
141                assert_eq!(
142                    num_faces_on_interior(k, n, d),
143                    count,
144                    "failed for k = {}, n = {}, d = {}",
145                    k,
146                    n,
147                    d
148                );
149            }
150        }
151    }
152}
153
154impl<Proj: Eq + Clone + BaseTriProjector> crate::Sphere for HexSphere<Proj> {
155    type Face = Face<Proj>;
156    type Vertex = Vertex<Proj>;
157    type HalfEdge = HalfEdge<Proj>;
158
159    fn num_faces(&self) -> usize {
160        self.num_vertices() / 2 + 2
161    }
162
163    fn face(&self, index: usize) -> Face<Proj> {
164        // TODO: "Real" implementation with constant-time performance
165        self.faces().nth(index).expect("face index out of bounds")
166    }
167
168    fn faces(&self) -> impl Iterator<Item = Face<Proj>> {
169        FaceIter {
170            sphere: self.clone(),
171            region: self.kis.base().first_region(),
172            u: 0,
173            u_end: self.kis.b(),
174            v: 0,
175        }
176    }
177
178    fn face_at(&self, point: [f64; 3]) -> Face<Proj> {
179        unsafe { Face::from_kis(self.kis.face_at(point)).unwrap_unchecked() }
180    }
181
182    fn num_vertices(&self) -> usize {
183        self.dual_num_divisions() * self.kis.base().num_faces()
184    }
185
186    fn vertex(&self, index: usize) -> Vertex<Proj> {
187        todo!()
188    }
189
190    fn vertices(&self) -> impl Iterator<Item = Vertex<Proj>> {
191        VertexIter {
192            sphere: self.clone(),
193            region: self.kis.base().first_region(),
194            u: 1,
195            u_end: self.kis.b(),
196            v: 0,
197            skip: false,
198        }
199    }
200}
201
202/// Represents a face on a [`HexSphere`].
203#[derive(Clone, Copy, Debug, PartialEq, Eq, Hash)]
204pub struct Face<Proj> {
205    center: tri::Vertex<Proj>,
206}
207
208impl<Proj: Eq + Clone + BaseTriProjector> Face<Proj> {
209    /// Gets the hexsphere face whose central vertex, in the [kised](HexSphere::kis) [`TriSphere`],
210    /// is the given [`Vertex`](tri::Vertex).
211    ///
212    /// This will return [`None`] if the given vertex does not correspond to a face on any
213    /// [`HexSphere`].
214    ///
215    /// # Examples
216    /// ```
217    /// # use subsphere::prelude::*;
218    /// let sphere = subsphere::icosphere().truncate();
219    /// let face = sphere.face(0);
220    /// let center = face.center();
221    /// assert_eq!(subsphere::hex::Face::from_center(center), Some(face));
222    /// ```
223    pub fn from_center(center: tri::Vertex<Proj>) -> Option<Self> {
224        if center.sphere.b() % 3 != center.sphere.c() % 3 {
225            return None;
226        }
227        if center.u % 3 != center.adjusted_v() % 3 {
228            return None;
229        }
230        Some(Self { center })
231    }
232
233    /// Gets the hexsphere face which, when [kised](HexSphere::kis), produces the given
234    /// triangular [`Face`](tri::Face).
235    ///
236    /// This will return [`None`] if there is no way to produce the given face by performing a
237    /// kis operation.
238    ///
239    /// # Examples
240    /// ```
241    /// # use subsphere::prelude::*;
242    /// let sphere = subsphere::icosphere().truncate();
243    /// let face = sphere.face(0);
244    /// let kis_face = face.side(0).kis().inside();
245    /// assert_eq!(subsphere::hex::Face::from_kis(kis_face), Some(face));
246    /// ```
247    pub fn from_kis(kis: tri::Face<Proj>) -> Option<Self> {
248        if kis.sphere.b() % 3 != kis.sphere.c() % 3 {
249            return None;
250        }
251        let adj_v = if !kis.region.ty().is_edge() {
252            kis.sphere.c()
253        } else {
254            0
255        };
256        let [u, v] = match ((kis.u_0 + 2 * (kis.v_0 + adj_v)) % 3, kis.boundary_along_v) {
257            (0, _) => [kis.u_0, kis.v_0],
258            (1, _) => [kis.u_0, kis.v_0 + 1],
259            (_, false) => [kis.u_0 + 1, kis.v_0],
260            (_, true) => [kis.u_0 - 1, kis.v_0 + 1],
261        };
262        Some(unsafe {
263            Self::from_center(tri::Vertex::new(kis.sphere, kis.region, u, v)).unwrap_unchecked()
264        })
265    }
266
267    /// The [`HexSphere`] that this [`Face`] belongs to.
268    pub fn sphere(&self) -> HexSphere<Proj> {
269        HexSphere {
270            kis: self.center.sphere.clone(),
271        }
272    }
273
274    /// The central vertex of this face on the [kised](HexSphere::kis) [`TriSphere`].
275    pub fn center(self) -> tri::Vertex<Proj> {
276        self.center
277    }
278
279    /// Indicates whether this is a hexagonal face.
280    pub fn is_hex(&self) -> bool {
281        self.center.as_base().is_none()
282    }
283}
284
285impl<Proj: Eq + Clone + BaseTriProjector> crate::Face for Face<Proj> {
286    type Vertex = Vertex<Proj>;
287    type HalfEdge = HalfEdge<Proj>;
288
289    fn index(&self) -> usize {
290        let b = self.center.sphere.b() as usize;
291        let c = self.center.sphere.c() as usize;
292        if self.center.region.ty().is_edge() {
293            self.sphere().base_face_index(self.center.region)
294                + num_faces_on_edge(b, self.center.v as usize)
295                + (self.center.u as usize).div_ceil(3)
296        } else if c < b {
297            self.sphere().base_face_index(self.center.region)
298                + num_faces_on_interior(c, b - c, self.center.v as usize)
299                + (self.center.u as usize).div_ceil(3)
300        } else {
301            self.sphere().base_face_index(self.center.region) + 1
302        }
303    }
304
305    fn area(&self) -> f64 {
306        crate::util::poly_area(self.vertices().map(|v| v.pos()))
307    }
308
309    fn num_sides(&self) -> usize {
310        if let Some(base) = self.center.as_base() {
311            base.degree()
312        } else {
313            6
314        }
315    }
316
317    fn side(&self, index: usize) -> HalfEdge<Proj> {
318        HalfEdge {
319            kis: self.center.outgoing(index).next(),
320        }
321    }
322}
323
324/// Represents a vertex on a [`HexSphere`].
325#[derive(Clone, Copy, Debug, PartialEq, Eq, Hash)]
326pub struct Vertex<Proj> {
327    kis: tri::Vertex<Proj>,
328}
329
330impl<Proj> Vertex<Proj> {
331    /// Constructs a [`Vertex`] from the given [`tri::Vertex`].
332    fn new(kis: tri::Vertex<Proj>) -> Self {
333        debug_assert_ne!(kis.u % 3, kis.adjusted_v() % 3);
334        Self { kis }
335    }
336}
337
338impl<Proj> tri::Vertex<Proj> {
339    /// The `v` coordinate of this vertex, measured from the base vertex of its region, rather
340    /// than the local origin.
341    ///
342    /// The combination of `u` and `adjusted_v` can be used to tell whether a [`tri::Vertex`]
343    /// corresponds to a face or a vertex on a [`HexSphere`].
344    fn adjusted_v(&self) -> u32 {
345        let mut v = self.v;
346        if !self.region.ty().is_edge() {
347            v += self.sphere.c();
348        }
349        v
350    }
351}
352
353impl<Proj: Clone> Vertex<Proj> {
354    /// The [`HexSphere`] that this [`Vertex`] belongs to.
355    pub fn sphere(&self) -> HexSphere<Proj> {
356        HexSphere {
357            kis: self.kis.sphere.clone(),
358        }
359    }
360
361    /// The corresponding vertex on the [kised](HexSphere::kis) [`TriSphere`].
362    pub fn kis(self) -> tri::Vertex<Proj> {
363        self.kis
364    }
365}
366
367impl<Proj: Eq + Clone + BaseTriProjector> crate::Vertex for Vertex<Proj> {
368    type Face = Face<Proj>;
369    type HalfEdge = HalfEdge<Proj>;
370
371    fn index(&self) -> usize {
372        let b = self.kis.sphere.b() as usize;
373        let c = self.kis.sphere.c() as usize;
374        if self.kis.region.ty().is_edge() {
375            self.sphere().base_vertex_index(self.kis.region) + self.kis.v as usize * (b - 1)
376                - num_faces_on_edge(b, self.kis.v as usize)
377                + (self.kis.u as usize * 2 - 1 - (self.kis.v % 3 == 1) as usize) / 3
378        } else if c < b {
379            let n = b - c;
380            self.sphere().base_vertex_index(self.kis.region)
381                + (self.kis.v as usize - 1) * (2 * n - self.kis.v as usize - 2) / 2
382                - num_faces_on_interior(c, n, self.kis.v as usize)
383                + (self.kis.u as usize * 2 - 1 - ((self.kis.v as usize + c) % 3 == 1) as usize) / 3
384        } else {
385            self.sphere().base_vertex_index(self.kis.region)
386        }
387    }
388
389    fn pos(&self) -> [f64; 3] {
390        self.kis.pos()
391    }
392
393    fn degree(&self) -> usize {
394        3
395    }
396
397    fn outgoing(&self, index: usize) -> HalfEdge<Proj> {
398        assert!(index < 3, "index out of bounds");
399        let phase = (self.kis.u + 2 * self.kis.adjusted_v()) % 3;
400        debug_assert_ne!(phase, 0);
401        HalfEdge {
402            kis: tri::HalfEdge::new(
403                self.kis.sphere.clone(),
404                self.kis.region,
405                self.kis.u,
406                self.kis.v,
407                tri::HalfEdgeDir::from_index(2 * index + phase as usize - 1),
408            ),
409        }
410    }
411}
412
413/// Represents one "side" of an edge on a [`HexSphere`].
414#[derive(Clone, Copy, Debug, PartialEq, Eq, Hash)]
415pub struct HalfEdge<Proj> {
416    kis: tri::HalfEdge<Proj>,
417}
418
419impl<Proj> HalfEdge<Proj> {
420    /// The corresponding half-edge on the [kised](HexSphere::kis) [`TriSphere`].
421    pub fn kis(self) -> tri::HalfEdge<Proj> {
422        self.kis
423    }
424}
425
426impl<Proj: Eq + Clone + BaseTriProjector> crate::HalfEdge for HalfEdge<Proj> {
427    type Face = Face<Proj>;
428    type Vertex = Vertex<Proj>;
429
430    fn side_index(&self) -> usize {
431        self.kis.prev().outgoing_index()
432    }
433
434    fn length(&self) -> f64 {
435        self.kis.length()
436    }
437
438    fn angle(&self) -> f64 {
439        let v_a = self.prev().start().pos();
440        let v_b = self.start().pos();
441        let v_c = self.next().start().pos();
442        crate::util::angle(v_a, v_b, v_c)
443    }
444
445    fn inside(&self) -> Face<Proj> {
446        Face {
447            center: self.kis.prev().start(),
448        }
449    }
450
451    fn start(&self) -> Vertex<Proj> {
452        Vertex {
453            kis: self.kis.start(),
454        }
455    }
456
457    fn complement(&self) -> Self {
458        Self {
459            kis: self.kis.complement(),
460        }
461    }
462
463    fn prev(&self) -> Self {
464        Self {
465            kis: self.kis.prev().complement().prev(),
466        }
467    }
468
469    fn next(&self) -> Self {
470        Self {
471            kis: self.kis.next().complement().next(),
472        }
473    }
474}
475
476impl<Proj> HexSphere<Proj> {
477    /// Gets the index of the first face which is owned by the given base region,
478    /// not counting the one corresponding to the first vertex of the base shape.
479    fn base_face_index(&self, region: BaseRegion) -> usize {
480        let face = region.owner();
481        let num_owned_vertices_before = face.num_owned_vertices_before()
482            + (face.owns_vertex_1() && region.ty() > BaseRegionType::Edge0) as usize;
483        let num_edge_regions_before =
484            face.num_owned_edges_before() + (region.ty() > BaseRegionType::Edge0) as usize;
485        let num_interior_regions_before =
486            face.index() + (region.ty() > BaseRegionType::Interior) as usize;
487        num_owned_vertices_before
488            + num_edge_regions_before * self.num_faces_per_edge_region()
489            + num_interior_regions_before * self.num_faces_per_interior_region()
490    }
491
492    /// Gets the index of the first vertex which is owned by the given base region.
493    fn base_vertex_index(&self, region: BaseRegion) -> usize {
494        let face = region.owner();
495        let num_edge_regions_before =
496            face.num_owned_edges_before() + (region.ty() > BaseRegionType::Edge0) as usize;
497        let num_interior_regions_before =
498            face.index() + (region.ty() > BaseRegionType::Interior) as usize;
499        num_edge_regions_before * self.num_vertices_per_edge_region()
500            + num_interior_regions_before * self.num_vertices_per_interior_region()
501    }
502}
503
504impl<Proj> HexSphere<Proj> {
505    /// Iterates over the faces of this [`HexSphere`], starting with the given region.
506    fn faces_from(self, region: BaseRegion) -> FaceIter<Proj> {
507        let b = self.kis.b();
508        let c = self.kis.c();
509        if region.ty().is_edge() {
510            FaceIter {
511                sphere: self,
512                region,
513                u: 3,
514                u_end: b,
515                v: 0,
516            }
517        } else if c < b {
518            let n = b - c;
519            FaceIter {
520                sphere: self,
521                region,
522                u: 1 + c % 3,
523                u_end: n - 1,
524                v: 1,
525            }
526        } else {
527            // This is a special case. When `b == c` and `c % 3 == 0`, there is exactly one
528            // interior face, at (0, 0).
529            FaceIter {
530                sphere: self,
531                region,
532                u: c % 3,
533                u_end: 1,
534                v: 0,
535            }
536        }
537    }
538}
539
540/// An iterator over the faces of an [`TriSphere`].
541#[derive(Clone, Debug)]
542pub struct FaceIter<Proj> {
543    sphere: HexSphere<Proj>,
544    region: BaseRegion,
545    u: u32,
546    u_end: u32,
547    v: u32,
548}
549
550impl<Proj: Eq + Clone + BaseTriProjector> Iterator for FaceIter<Proj> {
551    type Item = Face<Proj>;
552    fn next(&mut self) -> Option<Face<Proj>> {
553        loop {
554            if self.u < self.u_end {
555                let res = unsafe {
556                    Face::from_center(tri::Vertex {
557                        sphere: self.sphere.kis.clone(),
558                        region: self.region,
559                        u: self.u,
560                        v: self.v,
561                    })
562                    .unwrap_unchecked()
563                };
564                self.u += 3;
565                return Some(res);
566            } else if self.region.ty().is_edge() {
567                if self.v < self.sphere.kis.c() {
568                    self.v += 1;
569                    self.u = 1 + (self.v + 2) % 3;
570                    continue;
571                } else if self.u <= self.u_end
572                    && self.region.ty() == BaseRegionType::Edge0
573                    && self.region.owner().owns_vertex_1()
574                {
575                    let res = unsafe {
576                        Face::from_center(tri::Vertex {
577                            sphere: self.sphere.kis.clone(),
578                            region: self.region,
579                            u: self.u,
580                            v: self.v,
581                        })
582                        .unwrap_unchecked()
583                    };
584                    self.u += 3;
585                    return Some(res);
586                }
587            } else {
588                let n = self.sphere.kis.b() - self.sphere.kis.c();
589                if self.v + 1 < n {
590                    self.v += 1;
591                    self.u = 1 + (self.v + self.sphere.kis.c() + 2) % 3;
592                    self.u_end = n - self.v;
593                    continue;
594                }
595            }
596            if let Some(region) = self.sphere.kis.base().next_region(self.region) {
597                *self = self.sphere.clone().faces_from(region);
598                continue;
599            } else {
600                return None;
601            }
602        }
603    }
604}
605
606impl<Proj: Eq + Clone + BaseTriProjector> HexSphere<Proj> {
607    /// Iterates over the vertices of this [`HexSphere`], starting with the given region.
608    fn vertices_from(self, region: BaseRegion) -> VertexIter<Proj> {
609        let b = self.kis.b();
610        let c = self.kis.c();
611        if region.ty().is_edge() {
612            VertexIter {
613                sphere: self,
614                region,
615                u: 1,
616                u_end: b,
617                v: 0,
618                skip: false,
619            }
620        } else if c < b {
621            let n = b - c;
622            VertexIter {
623                sphere: self,
624                region,
625                u: 1 + (c % 3 == 0) as u32,
626                u_end: n.saturating_sub(1),
627                v: 1,
628                skip: c % 3 == 1,
629            }
630        } else {
631            // This is a special case. When `b == c` and `c % 3 > 0`, there is exactly one
632            // interior vertex, at (0, 0).
633            VertexIter {
634                sphere: self,
635                region,
636                u: (c % 3 == 0) as u32,
637                u_end: 1,
638                v: 0,
639                skip: c % 3 == 2,
640            }
641        }
642    }
643}
644
645/// An iterator over the vertices of an [`HexSphere`].
646#[derive(Clone, Debug)]
647pub struct VertexIter<Proj> {
648    sphere: HexSphere<Proj>,
649    region: BaseRegion,
650    u: u32,
651    u_end: u32,
652    v: u32,
653    skip: bool,
654}
655
656impl<Proj: Eq + Clone + BaseTriProjector> Iterator for VertexIter<Proj> {
657    type Item = Vertex<Proj>;
658    fn next(&mut self) -> Option<Vertex<Proj>> {
659        loop {
660            if self.u < self.u_end {
661                let res = Vertex::new(tri::Vertex {
662                    sphere: self.sphere.kis.clone(),
663                    region: self.region,
664                    u: self.u,
665                    v: self.v,
666                });
667                self.u += 1;
668                self.u += self.skip as u32;
669                self.skip = !self.skip;
670                return Some(res);
671            } else if self.region.ty().is_edge() {
672                if self.v < self.sphere.kis.c() {
673                    self.v += 1;
674                    (self.u, self.skip) = match self.v % 3 {
675                        0 => (1, false),
676                        1 => (2, false),
677                        _ => (1, true),
678                    };
679                    continue;
680                }
681            } else {
682                let n = self.sphere.kis.b() - self.sphere.kis.c();
683                if self.v + 1 < n {
684                    self.v += 1;
685                    self.u_end = n - self.v;
686                    (self.u, self.skip) = match (self.v + self.sphere.kis.c()) % 3 {
687                        0 => (1, false),
688                        1 => (2, false),
689                        _ => (1, true),
690                    };
691                    continue;
692                }
693            }
694            if let Some(region) = self.sphere.kis.base().next_region(self.region) {
695                *self = self.sphere.clone().vertices_from(region);
696                continue;
697            } else {
698                return None;
699            }
700        }
701    }
702}
subsphere/hex.rs

subsphere/
hex.rs
