dmc/
dmc.rs

1//! Contains the Dual Marching Cubes implementation, along with a function
2//! ([`mesh_from_octree`]) to interface with it.
3//!
4//! # Explanation
5//! The goal here is to adapt the Marching Cubes algorithm to hashed octrees using Morton keys.
6//! The [original Marching Cubes algorithm](http://paulbourke.net/geometry/polygonise/) used an
7//! uniform grid to generate the final mesh. However, with this approach, we use the dual grid
8//! generated from the octree given. The dual grid is composed of vertices (One at the center of
9//! each leaf node) and edges, which connect vertices of adjacent leaves.
10//!
11//! # References
12//! Refer to the comments at the start of [`src/lib.rs`](https://github.com/alexdevteam/dmc/blob/main/src/lib.rs).
13
14use crate::duals::DualGrid;
15use crate::octree::*;
16use crate::tables;
17use crate::util::*;
18
19use cgmath::{Point3, Vector3};
20use rayon::prelude::*;
21use std::mem::MaybeUninit;
22
23/// Creates a mesh from an octree full of values sampled from a Signed Distance Function (SDF).
24/// The triangles of the resulting mesh will be in **counter-clockwise** order.
25/// # Panics
26/// Should never panic. If it does, it's an error in the crate; please report it.
27pub fn mesh_from_octree<T: McNode>(octree: &HashedOctree<T>, scale: f32) -> Mesh {
28    let dual = DualGrid::from_octree(octree);
29
30    // Let's process the dual volumes one by one, concurrently.
31    let vertices = dual
32        .volumes
33        .into_par_iter()
34        .map(|volume| {
35            let VolumeData { nodes, cube_index } = fetch_volume_data(volume, &octree);
36            (volume, nodes, cube_index)
37        })
38        .filter(|(_, _, cube_index)| is_cell_visible(*cube_index))
39        .map(|(volume, nodes, cube_index)| {
40            let cell_data = calculate_mc_vertex_data(volume, nodes, cube_index, scale);
41            (cell_data, cube_index)
42        })
43        .flat_map_iter(|(cell_data, cube_index)| obtain_mc_vertices(cell_data, cube_index))
44        .collect::<Vec<_>>();
45
46    // Since the algorithm doesn't merge vertices yet, the indices to use are just the items in
47    // the `vertices` vector incrementally.
48    let indices = (0..vertices.len() as u32).collect();
49
50    Mesh { indices, vertices }
51}
52
53/// A simple mesh type that holds a vector of vertices and another one of indices.
54/// This type is meant to be converted to your own mesh type via the [`std::convert::From`] trait.
55#[derive(Clone, Debug, PartialEq)]
56pub struct Mesh {
57    /// The vertices of the mesh.
58    pub vertices: Vec<Vertex>,
59    /// The vertex indices of the mesh.
60    pub indices: Vec<u32>,
61}
62
63/// A 3D vertex that holds a position and a normal.
64#[derive(Clone, Copy, Debug, PartialEq)]
65pub struct Vertex {
66    /// The position of the vertex.
67    pub position: cgmath::Point3<f32>,
68    /// The normalized normal of the vertex.
69    pub normal: cgmath::Vector3<f32>,
70}
71
72/// Provides the functions that a node must define to create a mesh out of a group of them.
73pub trait McNode: Copy + Send + Sync {
74    /// This function should return the distance from this node's center to its contour. This can be
75    /// precalculated using a Signed Distance Function (SDF).
76    ///
77    /// # Resources
78    /// [Here](https://iquilezles.org/www/articles/distfunctions/distfunctions.htm) is a link with
79    /// basic SDFs and operations over them.
80    fn density(&self) -> f32;
81
82    /// This function should return the gradient of the SDF sampled, that is, the normal of the node
83    /// if it were on the surface of the SDF.
84    ///
85    /// # How To Calculate The Gradient
86    /// The gradient of a SDF is its derivative. However, if you don't want to calculate it (Which
87    /// you won't, because it will get complex really quickly), you can approximate it with the
88    /// following formula:
89    /// ```rust
90    /// use cgmath::*;
91    ///
92    /// fn sdf(pos: Point3<f32>) -> f32 {
93    ///     // Calculate SDF here...
94    ///     pos.distance(point3(0.,0.,0.)) - 0.4
95    /// }
96    ///
97    /// pub fn sdf_gradient(pos: Point3<f32>) -> Vector3<f32> {
98    ///     vec3(
99    ///        sdf(point3(pos.x + f32::EPSILON, pos.y, pos.z))
100    ///            - sdf(point3(pos.x - f32::EPSILON, pos.y, pos.z)),
101    ///        sdf(point3(pos.x, pos.y + f32::EPSILON, pos.z))
102    ///             - sdf(point3(pos.x, pos.y - f32::EPSILON, pos.z)),
103    ///        sdf(point3(pos.x, pos.y, pos.z + f32::EPSILON))
104    ///             - sdf(point3(pos.x, pos.y, pos.z - f32::EPSILON)),
105    ///     )
106    ///     .normalize()
107    /// }
108    /// ```
109    fn normal(&self) -> Vector3<f32>;
110}
111
112struct CellVertexData {
113    pub positions: [Point3<f32>; 12],
114    pub normals: [Vector3<f32>; 12],
115}
116
117fn obtain_mc_vertices(
118    cell_data: CellVertexData,
119    cube_index: u8,
120) -> impl std::iter::Iterator<Item = Vertex> {
121    tables::TRI_TABLE[cube_index as usize]
122        .chunks(3)
123        .take_while(|chunk| chunk[0] != -1)
124        .flat_map(move |chunk| {
125            let (vtx1, vtx2, vtx3) = (
126                cell_data.positions[chunk[0] as usize],
127                cell_data.positions[chunk[1] as usize],
128                cell_data.positions[chunk[2] as usize],
129            );
130            let (n1, n2, n3) = (
131                cell_data.normals[chunk[0] as usize],
132                cell_data.normals[chunk[1] as usize],
133                cell_data.normals[chunk[2] as usize],
134            );
135
136            std::array::IntoIter::new([(vtx1, n1), (vtx2, n2), (vtx3, n3)]).map(
137                move |(vtx, norm)| Vertex {
138                    position: vtx.into(),
139                    normal: norm.into(),
140                },
141            )
142        })
143}
144
145struct VolumeData {
146    pub nodes: [NodeData; 8],
147    pub cube_index: u8,
148}
149
150struct NodeData {
151    pub density: f32,
152    pub normal: Vector3<f32>,
153}
154
155fn calculate_mc_vertex_data(
156    volume: [MortonKey; 8],
157    nodes: [NodeData; 8],
158    cube_index: u8,
159    vertex_position_scale: f32,
160) -> CellVertexData {
161    /// Bindings from volume vertex (index) to Marching Cubes corner vertex (item).
162    const MC_CUBE_EDGES: [(usize, usize); 12] = [
163        (4, 5),
164        (1, 5),
165        (0, 1),
166        (0, 4),
167        (6, 7),
168        (3, 7),
169        (2, 3),
170        (2, 6),
171        (4, 6),
172        (5, 7),
173        (1, 3),
174        (0, 2),
175    ];
176
177    let mut positions: [Point3<f32>; 12] = unsafe { std::mem::MaybeUninit::uninit().assume_init() };
178    let mut normals: [Vector3<f32>; 12] = unsafe { std::mem::MaybeUninit::uninit().assume_init() };
179    let edges = tables::EDGE_TABLE[cube_index as usize];
180
181    (0..12)
182        .into_iter()
183        .filter(|i| (edges & (1 << i)) > 0)
184        .for_each(|i| {
185            let (v1_i, v2_i) = MC_CUBE_EDGES[i];
186            let factor = interpolation_factor(nodes[v1_i].density, nodes[v2_i].density);
187            positions[i] = (volume[v1_i]
188                .position()
189                .interpolate(volume[v2_i].position(), factor)
190                * vertex_position_scale)
191                .into();
192            normals[i] = nodes[v1_i]
193                .normal
194                .interpolate(nodes[v2_i].normal, factor)
195                .into();
196        });
197
198    CellVertexData { positions, normals }
199}
200
201fn is_cell_visible(cube_index: u8) -> bool {
202    cube_index != 0 && cube_index != 0xFF
203}
204
205fn fetch_volume_data<T: McNode>(volume: [MortonKey; 8], octree: &HashedOctree<T>) -> VolumeData {
206    const VOLUME_TO_MC_VERTEX: [usize; 8] = [3, 2, 7, 6, 0, 1, 4, 5];
207
208    let mut nodes: [MaybeUninit<NodeData>; 8] = unsafe { MaybeUninit::uninit().assume_init() };
209
210    let mut cube_index = 0;
211    volume.iter().enumerate().for_each(|(volume_i, &key)| {
212        let mc_i = VOLUME_TO_MC_VERTEX[volume_i];
213
214        if key == MortonKey::none() {
215            // If any of the vertices of the volume are null, that means that the volume is
216            // in the edge of the octree, and as such should not be processed since there
217            // is not enough data to form a marching cubes mesh out of it.
218            return;
219        }
220
221        match octree.value(key) {
222            Some(node) => {
223                let (density, normal) = (node.density(), node.normal());
224                if density > 0. {
225                    cube_index |= 1 << mc_i;
226                }
227                nodes[volume_i] = MaybeUninit::new(NodeData { density, normal });
228            }
229
230            // This node should ALWAYS exist, since the duals algorithm always returns
231            // either 0 (MortonKey::none()) or valid nodes as volume vertices.
232            None => unreachable!(),
233        }
234    });
235
236    let nodes: [NodeData; 8] = unsafe { std::mem::transmute(nodes) };
237
238    VolumeData { nodes, cube_index }
239}
240
241/// Given two values p1 and p2, estimates a relative value R such that `p1 + (p2 - p1) * R == 0`.
242/// Used to estimate where to place marching cube vertices to match up the isosurface.
243/// [Visual example](https://editor.p5js.org/alexinfdev/sketches/ldqHrNcr8)
244fn interpolation_factor(p1: f32, p2: f32) -> f32 {
245    if p1.abs() < 0.00001 {
246        p1
247    } else if p2.abs() < 0.00001 {
248        p2
249    } else if (p1 - p2).abs() < 0.00001 {
250        p1
251    } else {
252        (0. - p1) / (p2 - p1)
253    }
254}
dmc/dmc.rs

dmc/
dmc.rs
