conspire/geometry/mesh/tessellation/sdf/
mod.rs

1#[cfg(test)]
2#[cfg(feature = "netcdf")]
3mod test;
4
5use std::{
6    f64::consts::TAU,
7    thread::{available_parallelism, scope},
8};
9
10use crate::{
11    geometry::{Coordinate, mesh::tessellation::Tessellation},
12    math::{Scalar, Tensor, Vector},
13};
14
15impl Tessellation {
16    pub fn shape_diameter_function(
17        &self,
18        half_angle: Scalar,
19        rings: usize,
20        azimuthal: usize,
21    ) -> Vector {
22        let mesh = self.mesh();
23        let bvh = self.bvh();
24        let elements: Vec<&[usize]> = mesh.connectivities().iter().flatten().collect();
25        let coordinates = mesh.coordinates();
26        let centroids = mesh.centroids();
27        let normals: Vec<&Coordinate<3>> = self.normals.iter().flatten().collect();
28        let number_of_faces = normals.len();
29        let mut face_diameters = vec![0.0; number_of_faces];
30        let threads = available_parallelism().map_or(1, |threads| threads.get());
31        let chunk_size = number_of_faces.div_ceil(threads).max(1);
32        scope(|scope| {
33            let (bvh, elements, centroids, normals) = (bvh, &elements, &centroids, &normals);
34            face_diameters
35                .chunks_mut(chunk_size)
36                .enumerate()
37                .for_each(|(chunk, diameters)| {
38                    scope.spawn(move || {
39                        let offset = chunk * chunk_size;
40                        diameters
41                            .iter_mut()
42                            .enumerate()
43                            .for_each(|(local, diameter)| {
44                                let face = offset + local;
45                                let samples =
46                                    cone_directions(&-normals[face], half_angle, rings, azimuthal)
47                                        .into_iter()
48                                        .filter_map(|(direction, weight)| {
49                                            let ray = (centroids[face].clone(), direction).into();
50                                            bvh.intersect(&ray, coordinates, elements)
51                                                .filter(|hit| hit.index() != face)
52                                                .map(|hit| (hit.distance(), weight))
53                                        })
54                                        .collect();
55                                *diameter = weighted_diameter(samples);
56                            });
57                    });
58                });
59        });
60        interpolate_to_nodes(face_diameters.into(), elements, coordinates.len())
61    }
62}
63
64fn interpolate_to_nodes(
65    face_diameters: Vector,
66    elements: Vec<&[usize]>,
67    number_of_nodes: usize,
68) -> Vector {
69    let mut nodal = Vector::zero(number_of_nodes);
70    let mut counts = vec![0; number_of_nodes];
71    elements
72        .into_iter()
73        .zip(face_diameters)
74        .for_each(|(element, diameter)| {
75            element.iter().for_each(|&node| {
76                nodal[node] += diameter;
77                counts[node] += 1;
78            })
79        });
80    nodal.iter_mut().zip(counts).for_each(|(value, count)| {
81        if count > 0 {
82            *value /= count as Scalar
83        }
84    });
85    nodal
86}
87
88fn cone_directions(
89    axis: &Coordinate<3>,
90    half_angle: Scalar,
91    rings: usize,
92    azimuthal: usize,
93) -> Vec<(Coordinate<3>, Scalar)> {
94    let basis = axis.orthonormal_basis();
95    let (axis, tangent_1, tangent_2) = (&basis[0], &basis[1], &basis[2]);
96    let mut directions = Vec::with_capacity(1 + rings * azimuthal);
97    directions.push((axis.clone(), 1.0));
98    for ring in 1..=rings {
99        let polar = half_angle * ring as Scalar / rings as Scalar;
100        let (sin_polar, cos_polar) = polar.sin_cos();
101        for sample in 0..azimuthal {
102            let (sin_azimuth, cos_azimuth) =
103                (TAU * sample as Scalar / azimuthal as Scalar).sin_cos();
104            let direction = axis * cos_polar
105                + tangent_1 * (sin_polar * cos_azimuth)
106                + tangent_2 * (sin_polar * sin_azimuth);
107            directions.push((direction, cos_polar));
108        }
109    }
110    directions
111}
112
113fn weighted_diameter(samples: Vec<(Scalar, Scalar)>) -> Scalar {
114    if samples.is_empty() {
115        return 0.0;
116    }
117    let mut distances: Vec<Scalar> = samples.iter().map(|&(distance, _)| distance).collect();
118    distances.sort_by(|a, b| a.partial_cmp(b).unwrap());
119    let median = distances[distances.len() / 2];
120    let mean = distances.iter().sum::<Scalar>() / distances.len() as Scalar;
121    let standard_deviation = (distances
122        .iter()
123        .map(|distance| (distance - mean).powi(2))
124        .sum::<Scalar>()
125        / distances.len() as Scalar)
126        .sqrt();
127    let (numerator, denominator) = samples
128        .into_iter()
129        .filter(|&(distance, _)| (distance - median).abs() <= standard_deviation)
130        .fold(
131            (0.0, 0.0),
132            |(numerator, denominator), (distance, weight)| {
133                (numerator + weight * distance, denominator + weight)
134            },
135        );
136    if denominator > 0.0 {
137        numerator / denominator
138    } else {
139        median
140    }
141}
conspire/geometry/mesh/tessellation/sdf/mod.rs

conspire/geometry/mesh/tessellation/sdf/
mod.rs
