fast_surface_nets/lib.rs
1//! A fast, chunk-friendly implementation of Naive Surface Nets on regular grids.
2//!
3//! 
5//!
6//! Surface Nets is an algorithm for extracting an isosurface mesh from a [signed distance
7//! field](https://en.wikipedia.org/wiki/Signed_distance_function) sampled on a regular grid. It is nearly the same as Dual
8//! Contouring, but instead of using hermite (derivative) data to estimate surface points, Surface Nets will do a simpler form
9//! of interpolation (average) between points where the isosurface crosses voxel cube edges.
10//!
11//! Benchmarks show that [`surface_nets`] generates about 20 million triangles per second on a single core
12//! of a 2.5 GHz Intel Core i7. This implementation achieves high performance by using small lookup tables and SIMD acceleration
13//! provided by `glam` when doing 3D floating point vector math. (Users are not required to use `glam` types in any API
14//! signatures.) To run the benchmarks yourself, `cd bench/ && cargo bench`.
15//!
16//! High-quality surface normals are estimated by:
17//!
18//! 1. calculating SDF derivatives using central differencing
19//! 2. using bilinear interpolation of SDF derivatives along voxel cube edges
20//!
21//! When working with sparse data sets, [`surface_nets`] can generate meshes for array chunks that fit
22//! together seamlessly. This works because faces are not generated on the positive boundaries of a chunk. One must only apply a
23//! translation of the mesh into proper world coordinates for the given chunk.
24//!
25//! # Example Code
26//!
27//! ```
28//! use fast_surface_nets::ndshape::{ConstShape, ConstShape3u32};
29//! use fast_surface_nets::{surface_nets, SurfaceNetsBuffer};
30//!
31//! // A 16^3 chunk with 1-voxel boundary padding.
32//! type ChunkShape = ConstShape3u32<18, 18, 18>;
33//!
34//! // This chunk will cover just a single octant of a sphere SDF (radius 15).
35//! let mut sdf = [1.0; ChunkShape::USIZE];
36//! for i in 0u32..ChunkShape::SIZE {
37//! let [x, y, z] = ChunkShape::delinearize(i);
38//! sdf[i as usize] = ((x * x + y * y + z * z) as f32).sqrt() - 15.0;
39//! }
40//!
41//! let mut buffer = SurfaceNetsBuffer::default();
42//! surface_nets(&sdf, &ChunkShape {}, [0; 3], [17; 3], &mut buffer);
43//!
44//! // Some triangles were generated.
45//! assert!(!buffer.indices.is_empty());
46//! ```
47
48pub use glam;
49pub use ndshape;
50
51use glam::{Vec3A, Vec3Swizzles};
52use ndshape::Shape;
53
54pub trait SignedDistance: Into<f32> + Copy {
55 fn is_negative(self) -> bool;
56}
57
58impl SignedDistance for f32 {
59 fn is_negative(self) -> bool {
60 self < 0.0
61 }
62}
63
64/// The output buffers used by [`surface_nets`]. These buffers can be reused to avoid reallocating memory.
65#[derive(Default, Clone)]
66pub struct SurfaceNetsBuffer {
67 /// The triangle mesh positions.
68 ///
69 /// These are in array-local coordinates, i.e. at array position `(x, y, z)`, the vertex position would be `(x, y, z) +
70 /// centroid` if the isosurface intersects that voxel.
71 pub positions: Vec<[f32; 3]>,
72 /// The triangle mesh normals.
73 ///
74 /// The normals are **not** normalized, since that is done most efficiently on the GPU.
75 pub normals: Vec<[f32; 3]>,
76 /// The triangle mesh indices.
77 pub indices: Vec<u32>,
78
79 /// Local 3D array coordinates of every voxel that intersects the isosurface.
80 pub surface_points: Vec<[u32; 3]>,
81 /// Stride of every voxel that intersects the isosurface. Can be used for efficient post-processing.
82 pub surface_strides: Vec<u32>,
83 /// Used to map back from voxel stride to vertex index.
84 pub stride_to_index: Vec<u32>,
85}
86
87impl SurfaceNetsBuffer {
88 /// Clears all of the buffers, but keeps the memory allocated for reuse.
89 fn reset(&mut self, array_size: usize) {
90 self.positions.clear();
91 self.normals.clear();
92 self.indices.clear();
93 self.surface_points.clear();
94 self.surface_strides.clear();
95
96 // Just make sure this buffer is big enough, whether or not we've used it before.
97 self.stride_to_index.resize(array_size, NULL_VERTEX);
98 }
99}
100
101/// This stride of the SDF array did not produce a vertex.
102pub const NULL_VERTEX: u32 = u32::MAX;
103
104/// The Naive Surface Nets smooth voxel meshing algorithm.
105///
106/// Extracts an isosurface mesh from the [signed distance field](https://en.wikipedia.org/wiki/Signed_distance_function) `sdf`.
107/// Each value in the field determines how close that point is to the isosurface. Negative values are considered "interior" of
108/// the surface volume, and positive values are considered "exterior." These lattice points will be considered corners of unit
109/// cubes. For each unit cube, at most one isosurface vertex will be estimated, as below, where `p` is a positive corner value,
110/// `n` is a negative corner value, `s` is an isosurface vertex, and `|` or `-` are mesh polygons connecting the vertices.
111///
112/// ```text
113/// p p p p
114/// s---s
115/// p | n | p p
116/// s s---s
117/// p | n n | p
118/// s---s---s
119/// p p p p
120/// ```
121///
122/// The set of corners sampled is exactly the set of points in `[min, max]`. `sdf` must contain all of those points.
123///
124/// Note that the scheme illustrated above implies that chunks must be padded with a 1-voxel border copied from neighboring
125/// voxels in order to connect seamlessly.
126pub fn surface_nets<T, S>(
127 sdf: &[T],
128 shape: &S,
129 min: [u32; 3],
130 max: [u32; 3],
131 output: &mut SurfaceNetsBuffer,
132) where
133 T: SignedDistance,
134 S: Shape<3, Coord = u32>,
135{
136 // SAFETY
137 // Make sure the slice matches the shape before we start using get_unchecked.
138 assert!(shape.linearize(min) <= shape.linearize(max));
139 assert!((shape.linearize(max) as usize) < sdf.len());
140
141 output.reset(sdf.len());
142
143 estimate_surface(sdf, shape, min, max, output);
144 make_all_quads(sdf, shape, min, max, output);
145}
146
147// Find all vertex positions and normals. Also generate a map from grid position to vertex index to be used to look up vertices
148// when generating quads.
149fn estimate_surface<T, S>(
150 sdf: &[T],
151 shape: &S,
152 [minx, miny, minz]: [u32; 3],
153 [maxx, maxy, maxz]: [u32; 3],
154 output: &mut SurfaceNetsBuffer,
155) where
156 T: SignedDistance,
157 S: Shape<3, Coord = u32>,
158{
159 for z in minz..maxz {
160 for y in miny..maxy {
161 for x in minx..maxx {
162 let stride = shape.linearize([x, y, z]);
163 let p = Vec3A::from([x as f32, y as f32, z as f32]);
164 if estimate_surface_in_cube(sdf, shape, p, stride, output) {
165 output.stride_to_index[stride as usize] = output.positions.len() as u32 - 1;
166 output.surface_points.push([x, y, z]);
167 output.surface_strides.push(stride);
168 } else {
169 output.stride_to_index[stride as usize] = NULL_VERTEX;
170 }
171 }
172 }
173 }
174}
175
176// Consider the grid-aligned cube where `p` is the minimal corner. Find a point inside this cube that is approximately on the
177// isosurface.
178//
179// This is done by estimating, for each cube edge, where the isosurface crosses the edge (if it does at all). Then the estimated
180// surface point is the average of these edge crossings.
181fn estimate_surface_in_cube<T, S>(
182 sdf: &[T],
183 shape: &S,
184 p: Vec3A,
185 min_corner_stride: u32,
186 output: &mut SurfaceNetsBuffer,
187) -> bool
188where
189 T: SignedDistance,
190 S: Shape<3, Coord = u32>,
191{
192 // Get the signed distance values at each corner of this cube.
193 let mut corner_dists = [0f32; 8];
194 let mut num_negative = 0;
195 for (i, dist) in corner_dists.iter_mut().enumerate() {
196 let corner_stride = min_corner_stride + shape.linearize(CUBE_CORNERS[i]);
197 let d = *unsafe { sdf.get_unchecked(corner_stride as usize) };
198 *dist = d.into();
199 if d.is_negative() {
200 num_negative += 1;
201 }
202 }
203
204 if num_negative == 0 || num_negative == 8 {
205 // No crossings.
206 return false;
207 }
208
209 let c = centroid_of_edge_intersections(&corner_dists);
210
211 output.positions.push((p + c).into());
212 output.normals.push(sdf_gradient(&corner_dists, c).into());
213
214 true
215}
216
217fn centroid_of_edge_intersections(dists: &[f32; 8]) -> Vec3A {
218 let mut count = 0;
219 let mut sum = Vec3A::ZERO;
220 for &[corner1, corner2] in CUBE_EDGES.iter() {
221 let d1 = dists[corner1 as usize];
222 let d2 = dists[corner2 as usize];
223 if (d1 < 0.0) != (d2 < 0.0) {
224 count += 1;
225 sum += estimate_surface_edge_intersection(corner1, corner2, d1, d2);
226 }
227 }
228
229 sum / count as f32
230}
231
232// Given two cube corners, find the point between them where the SDF is zero. (This might not exist).
233fn estimate_surface_edge_intersection(
234 corner1: u32,
235 corner2: u32,
236 value1: f32,
237 value2: f32,
238) -> Vec3A {
239 let interp1 = value1 / (value1 - value2);
240 let interp2 = 1.0 - interp1;
241
242 interp2 * CUBE_CORNER_VECTORS[corner1 as usize]
243 + interp1 * CUBE_CORNER_VECTORS[corner2 as usize]
244}
245
246/// Calculate the normal as the gradient of the distance field. Don't bother making it a unit vector, since we'll do that on the
247/// GPU.
248///
249/// For each dimension, there are 4 cube edges along that axis. This will do bilinear interpolation between the differences
250/// along those edges based on the position of the surface (s).
251fn sdf_gradient(dists: &[f32; 8], s: Vec3A) -> Vec3A {
252 let p00 = Vec3A::from([dists[0b001], dists[0b010], dists[0b100]]);
253 let n00 = Vec3A::from([dists[0b000], dists[0b000], dists[0b000]]);
254
255 let p10 = Vec3A::from([dists[0b101], dists[0b011], dists[0b110]]);
256 let n10 = Vec3A::from([dists[0b100], dists[0b001], dists[0b010]]);
257
258 let p01 = Vec3A::from([dists[0b011], dists[0b110], dists[0b101]]);
259 let n01 = Vec3A::from([dists[0b010], dists[0b100], dists[0b001]]);
260
261 let p11 = Vec3A::from([dists[0b111], dists[0b111], dists[0b111]]);
262 let n11 = Vec3A::from([dists[0b110], dists[0b101], dists[0b011]]);
263
264 // Each dimension encodes an edge delta, giving 12 in total.
265 let d00 = p00 - n00; // Edges (0b00x, 0b0y0, 0bz00)
266 let d10 = p10 - n10; // Edges (0b10x, 0b0y1, 0bz10)
267 let d01 = p01 - n01; // Edges (0b01x, 0b1y0, 0bz01)
268 let d11 = p11 - n11; // Edges (0b11x, 0b1y1, 0bz11)
269
270 let neg = Vec3A::ONE - s;
271
272 // Do bilinear interpolation between 4 edges in each dimension.
273 neg.yzx() * neg.zxy() * d00
274 + neg.yzx() * s.zxy() * d10
275 + s.yzx() * neg.zxy() * d01
276 + s.yzx() * s.zxy() * d11
277}
278
279// For every edge that crosses the isosurface, make a quad between the "centers" of the four cubes touching that surface. The
280// "centers" are actually the vertex positions found earlier. Also make sure the triangles are facing the right way. See the
281// comments on `maybe_make_quad` to help with understanding the indexing.
282fn make_all_quads<T, S>(
283 sdf: &[T],
284 shape: &S,
285 [minx, miny, minz]: [u32; 3],
286 [maxx, maxy, maxz]: [u32; 3],
287 output: &mut SurfaceNetsBuffer,
288) where
289 T: SignedDistance,
290 S: Shape<3, Coord = u32>,
291{
292 let xyz_strides = [
293 shape.linearize([1, 0, 0]) as usize,
294 shape.linearize([0, 1, 0]) as usize,
295 shape.linearize([0, 0, 1]) as usize,
296 ];
297
298 for (&[x, y, z], &p_stride) in output
299 .surface_points
300 .iter()
301 .zip(output.surface_strides.iter())
302 {
303 let p_stride = p_stride as usize;
304 let eval_max_plane = cfg!(feature = "eval-max-plane");
305
306 // Do edges parallel with the X axis
307 if y != miny && z != minz && (eval_max_plane || x != maxx - 1) {
308 maybe_make_quad(
309 sdf,
310 &output.stride_to_index,
311 &output.positions,
312 p_stride,
313 p_stride + xyz_strides[0],
314 xyz_strides[1],
315 xyz_strides[2],
316 &mut output.indices,
317 );
318 }
319 // Do edges parallel with the Y axis
320 if x != minx && z != minz && (eval_max_plane || y != maxy - 1) {
321 maybe_make_quad(
322 sdf,
323 &output.stride_to_index,
324 &output.positions,
325 p_stride,
326 p_stride + xyz_strides[1],
327 xyz_strides[2],
328 xyz_strides[0],
329 &mut output.indices,
330 );
331 }
332 // Do edges parallel with the Z axis
333 if x != minx && y != miny && (eval_max_plane || z != maxz - 1) {
334 maybe_make_quad(
335 sdf,
336 &output.stride_to_index,
337 &output.positions,
338 p_stride,
339 p_stride + xyz_strides[2],
340 xyz_strides[0],
341 xyz_strides[1],
342 &mut output.indices,
343 );
344 }
345 }
346}
347
348// Construct a quad in the dual graph of the SDF lattice.
349//
350// The surface point s was found somewhere inside of the cube with minimal corner p1.
351//
352// x ---- x
353// / /|
354// x ---- x |
355// | s | x
356// | |/
357// p1 --- p2
358//
359// And now we want to find the quad between p1 and p2 where s is a corner of the quad.
360//
361// s
362// /|
363// / |
364// | |
365// p1 | | p2
366// | /
367// |/
368//
369// If A is (of the three grid axes) the axis between p1 and p2,
370//
371// A
372// p1 ---> p2
373//
374// then we must find the other 3 quad corners by moving along the other two axes (those orthogonal to A) in the negative
375// directions; these are axis B and axis C.
376#[allow(clippy::too_many_arguments)]
377fn maybe_make_quad<T>(
378 sdf: &[T],
379 stride_to_index: &[u32],
380 positions: &[[f32; 3]],
381 p1: usize,
382 p2: usize,
383 axis_b_stride: usize,
384 axis_c_stride: usize,
385 indices: &mut Vec<u32>,
386) where
387 T: SignedDistance,
388{
389 let d1 = unsafe { sdf.get_unchecked(p1) };
390 let d2 = unsafe { sdf.get_unchecked(p2) };
391 let negative_face = match (d1.is_negative(), d2.is_negative()) {
392 (true, false) => false,
393 (false, true) => true,
394 _ => return, // No face.
395 };
396
397 // The triangle points, viewed face-front, look like this:
398 // v1 v3
399 // v2 v4
400 let v1 = stride_to_index[p1];
401 let v2 = stride_to_index[p1 - axis_b_stride];
402 let v3 = stride_to_index[p1 - axis_c_stride];
403 let v4 = stride_to_index[p1 - axis_b_stride - axis_c_stride];
404 let (pos1, pos2, pos3, pos4) = (
405 Vec3A::from(positions[v1 as usize]),
406 Vec3A::from(positions[v2 as usize]),
407 Vec3A::from(positions[v3 as usize]),
408 Vec3A::from(positions[v4 as usize]),
409 );
410 // Split the quad along the shorter axis, rather than the longer one.
411 let quad = if pos1.distance_squared(pos4) < pos2.distance_squared(pos3) {
412 if negative_face {
413 [v1, v4, v2, v1, v3, v4]
414 } else {
415 [v1, v2, v4, v1, v4, v3]
416 }
417 } else if negative_face {
418 [v2, v3, v4, v2, v1, v3]
419 } else {
420 [v2, v4, v3, v2, v3, v1]
421 };
422 indices.extend_from_slice(&quad);
423}
424
425const CUBE_CORNERS: [[u32; 3]; 8] = [
426 [0, 0, 0],
427 [1, 0, 0],
428 [0, 1, 0],
429 [1, 1, 0],
430 [0, 0, 1],
431 [1, 0, 1],
432 [0, 1, 1],
433 [1, 1, 1],
434];
435const CUBE_CORNER_VECTORS: [Vec3A; 8] = [
436 Vec3A::from_array([0.0, 0.0, 0.0]),
437 Vec3A::from_array([1.0, 0.0, 0.0]),
438 Vec3A::from_array([0.0, 1.0, 0.0]),
439 Vec3A::from_array([1.0, 1.0, 0.0]),
440 Vec3A::from_array([0.0, 0.0, 1.0]),
441 Vec3A::from_array([1.0, 0.0, 1.0]),
442 Vec3A::from_array([0.0, 1.0, 1.0]),
443 Vec3A::from_array([1.0, 1.0, 1.0]),
444];
445const CUBE_EDGES: [[u32; 2]; 12] = [
446 [0b000, 0b001],
447 [0b000, 0b010],
448 [0b000, 0b100],
449 [0b001, 0b011],
450 [0b001, 0b101],
451 [0b010, 0b011],
452 [0b010, 0b110],
453 [0b011, 0b111],
454 [0b100, 0b101],
455 [0b100, 0b110],
456 [0b101, 0b111],
457 [0b110, 0b111],
458];