conspire/geometry/ntree/from/tessellation/
mod.rs

1use crate::{
2    geometry::{
3        Coordinate, CoordinateList,
4        bbox::BoundingBox,
5        mesh::Tessellation,
6        ntree::{
7            Octree,
8            balance::Balancing,
9            node::{Kind, Node, split::Split},
10            pair::Pairing,
11            rescale::Rescaling,
12        },
13    },
14    math::{Scalar, Tensor, TensorVec},
15};
16use std::{array::from_fn, f64::consts::FRAC_PI_3, ops::Add};
17
18const D: usize = 3;
19const M: usize = 6;
20
21impl<T, U> Octree<T, U>
22where
23    T: Add<Output = T> + Copy + From<u16> + Into<Scalar> + Into<usize> + PartialOrd + Split,
24    U: Copy + From<usize> + Into<usize>,
25{
26    pub fn from_sdf(tessellation: &Tessellation, scale: Scalar) -> Self {
27        let sdf = tessellation.shape_diameter_function(FRAC_PI_3, 3, 8);
28        let coordinates = tessellation.mesh().coordinates();
29        if coordinates.is_empty() {
30            return Self {
31                balanced: Balancing::None,
32                nodes: vec![Node {
33                    corner: from_fn(|_| T::from(0)),
34                    length: T::from(1),
35                    facets: [None; M],
36                    kind: Kind::Leaf,
37                    value: None,
38                }],
39                paired: Pairing::None,
40                rescale: Rescaling {
41                    center: [0.0; D],
42                    cell: 1.0,
43                    half: 0.0,
44                },
45            };
46        }
47        let mut min_coord: [f64; D] = from_fn(|_| f64::INFINITY);
48        let mut max_coord: [f64; D] = from_fn(|_| f64::NEG_INFINITY);
49        for point in coordinates {
50            for ax in 0..D {
51                min_coord[ax] = min_coord[ax].min(point[ax]);
52                max_coord[ax] = max_coord[ax].max(point[ax]);
53            }
54        }
55        let max_extent = (0..D)
56            .map(|ax| max_coord[ax] - min_coord[ax])
57            .fold(0.0f64, f64::max);
58        let min_sdf = sdf
59            .iter()
60            .copied()
61            .filter(|value| *value > 0.0)
62            .fold(f64::INFINITY, f64::min);
63        let min_length = if min_sdf.is_finite() {
64            min_sdf / scale
65        } else {
66            max_extent
67        };
68        let levels = if max_extent <= 0.0 || min_length <= 0.0 {
69            0u32
70        } else {
71            (max_extent / min_length).log2().ceil().max(0.0) as u32
72        };
73        let root_length: u16 = 1u16.checked_shl(levels).unwrap_or(u16::MAX);
74        let center: [f64; D] = from_fn(|ax| (min_coord[ax] + max_coord[ax]) / 2.0);
75        let mut tree = Self {
76            balanced: Balancing::None,
77            rescale: Rescaling {
78                center,
79                cell: min_length,
80                half: root_length as Scalar / 2.0,
81            },
82            nodes: vec![Node {
83                corner: from_fn(|_| T::from(0)),
84                length: T::from(root_length),
85                facets: [None; M],
86                kind: Kind::Leaf,
87                value: None,
88            }],
89            paired: Pairing::None,
90        };
91        let elements: Vec<&[usize]> = tessellation
92            .mesh()
93            .connectivities()
94            .iter()
95            .flatten()
96            .collect();
97        let targets: Vec<Scalar> = elements
98            .iter()
99            .map(|element| sdf[element[0]].min(sdf[element[1]]).min(sdf[element[2]]))
100            .collect();
101        let half = root_length as Scalar / 2.0;
102        let overlaps = |bbox: &BoundingBox<3>, triangle: usize| {
103            let element = elements[triangle];
104            bbox.overlaps_triangle(
105                &coordinates[element[0]],
106                &coordinates[element[1]],
107                &coordinates[element[2]],
108            )
109        };
110        let mut stack: Vec<(usize, Vec<usize>)> = vec![(0, (0..elements.len()).collect())];
111        while let Some((index, overlapping)) = stack.pop() {
112            let cells: usize = tree.nodes[index].length.into();
113            let extent: Scalar = tree.nodes[index].length.into();
114            let target = overlapping
115                .iter()
116                .map(|&triangle| targets[triangle])
117                .fold(f64::INFINITY, f64::min);
118            if cells <= 1 || (extent * min_length) * scale <= target {
119                continue;
120            }
121            if tree.subdivide(U::from(index)).is_err() {
122                continue;
123            }
124            let children: Vec<usize> = tree.nodes[index]
125                .orthants()
126                .unwrap()
127                .iter()
128                .map(|&child| child.into())
129                .collect();
130            for child in children {
131                let corner = tree.nodes[child].corner;
132                let child_extent: Scalar = tree.nodes[child].length.into();
133                let minimum = Coordinate::const_from(from_fn(|ax| {
134                    center[ax] + (Into::<Scalar>::into(corner[ax]) - half) * min_length
135                }));
136                let maximum =
137                    Coordinate::const_from(from_fn(|ax| minimum[ax] + child_extent * min_length));
138                let bbox = BoundingBox::from(CoordinateList::const_from([minimum, maximum]));
139                let inside: Vec<usize> = overlapping
140                    .iter()
141                    .copied()
142                    .filter(|&triangle| overlaps(&bbox, triangle))
143                    .collect();
144                if !inside.is_empty() {
145                    stack.push((child, inside));
146                }
147            }
148        }
149        tree
150    }
151}
conspire/geometry/ntree/from/tessellation/mod.rs

conspire/geometry/ntree/from/tessellation/
mod.rs
