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

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

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