Skip to main content

cranpose_ui_graphics/space/
sdf.rs

1//! Meshes from signed distance fields: the surface where a field is zero,
2//! found on a grid by surface nets.
3
4use std::f32::consts::{PI, TAU};
5
6use super::{Mesh, MeshError, Vec3};
7
8/// The corners of a grid cell, as offsets along x, y and z.
9const CORNERS: [[usize; 3]; 8] = [
10    [0, 0, 0],
11    [1, 0, 0],
12    [0, 1, 0],
13    [1, 1, 0],
14    [0, 0, 1],
15    [1, 0, 1],
16    [0, 1, 1],
17    [1, 1, 1],
18];
19
20/// The edges of a grid cell, as pairs of its corners.
21const EDGES: [[usize; 2]; 12] = [
22    [0, 1],
23    [2, 3],
24    [4, 5],
25    [6, 7],
26    [0, 2],
27    [1, 3],
28    [4, 6],
29    [5, 7],
30    [0, 4],
31    [1, 5],
32    [2, 6],
33    [3, 7],
34];
35
36/// For each axis, the step along a grid edge, and the four cells around the
37/// edge, as offsets back from its start, in turn around the edge.
38const AROUND_EDGES: [([usize; 3], [[usize; 3]; 4]); 3] = [
39    ([1, 0, 0], [[0, 1, 1], [0, 0, 1], [0, 0, 0], [0, 1, 0]]),
40    ([0, 1, 0], [[1, 0, 1], [1, 0, 0], [0, 0, 0], [0, 0, 1]]),
41    ([0, 0, 1], [[1, 1, 0], [0, 1, 0], [0, 0, 0], [1, 0, 0]]),
42];
43
44/// A field sampled on the points of a grid of `cells` cells along each side
45/// of the box from `low` with cells of `step`.
46struct Samples {
47    values: Vec<f32>,
48    cells: usize,
49    low: Vec3,
50    step: Vec3,
51}
52
53impl Samples {
54    fn of(field: &impl Fn(Vec3) -> f32, (low, high): (Vec3, Vec3), cells: usize) -> Self {
55        let side = cells as f32;
56        let step = Vec3::new(
57            (high.x - low.x) / side,
58            (high.y - low.y) / side,
59            (high.z - low.z) / side,
60        );
61        let points = cells + 1;
62        let mut values = Vec::with_capacity(points * points * points);
63        for k in 0..points {
64            for j in 0..points {
65                for i in 0..points {
66                    values.push(field(Self::place(low, step, [i, j, k])));
67                }
68            }
69        }
70        Self {
71            values,
72            cells,
73            low,
74            step,
75        }
76    }
77
78    fn place(low: Vec3, step: Vec3, [i, j, k]: [usize; 3]) -> Vec3 {
79        Vec3::new(
80            low.x + step.x * i as f32,
81            low.y + step.y * j as f32,
82            low.z + step.z * k as f32,
83        )
84    }
85
86    fn point(&self, at: [usize; 3]) -> Vec3 {
87        Self::place(self.low, self.step, at)
88    }
89
90    fn value(&self, [i, j, k]: [usize; 3]) -> f32 {
91        let points = self.cells + 1;
92        self.values[(k * points + j) * points + i]
93    }
94
95    /// Where the surface crosses the cell at `cell`: the mean of the points
96    /// where it crosses the cell's edges, or `None` when it misses the cell.
97    fn cell_point(&self, [i, j, k]: [usize; 3]) -> Option<Vec3> {
98        let corner = |index: usize| {
99            let [di, dj, dk] = CORNERS[index];
100            [i + di, j + dj, k + dk]
101        };
102        let (mut sum, mut count) = (Vec3::ZERO, 0.0);
103        for [a, b] in EDGES {
104            let (va, vb) = (self.value(corner(a)), self.value(corner(b)));
105            if (va < 0.0) == (vb < 0.0) {
106                continue;
107            }
108            let (pa, pb) = (self.point(corner(a)), self.point(corner(b)));
109            sum = sum + pa + (pb - pa) * (va / (va - vb));
110            count += 1.0;
111        }
112        (count > 0.0).then(|| sum * (1.0 / count))
113    }
114}
115
116impl Mesh {
117    /// The surface where `field` is zero within the box from `bounds.0` to
118    /// `bounds.1`: a signed distance field, negative inside the shape,
119    /// sampled on a grid of `resolution` cells along each side of the box.
120    /// The normals follow the field's gradient, and the texture wraps the
121    /// shape as it wraps [`Mesh::sphere`], seen from the middle of the box:
122    /// its middle faces the viewer, its left and right edges meet at the
123    /// back and its top and bottom edges are the poles. A field that changes
124    /// gives a shape that changes: build the mesh again for each new shape.
125    pub fn from_sdf(
126        field: impl Fn(Vec3) -> f32,
127        bounds: (Vec3, Vec3),
128        resolution: u32,
129    ) -> Result<Self, MeshError> {
130        let cells = resolution.max(2) as usize;
131        let samples = Samples::of(&field, bounds, cells);
132        let middle = (bounds.0 + bounds.1) * 0.5;
133        let reach = samples.step.x.min(samples.step.y).min(samples.step.z) * 0.5;
134        let mut net = Net::default();
135        let mut cell_vertices = vec![u32::MAX; cells * cells * cells];
136        for k in 0..cells {
137            for j in 0..cells {
138                for i in 0..cells {
139                    if let Some(point) = samples.cell_point([i, j, k]) {
140                        cell_vertices[(k * cells + j) * cells + i] =
141                            net.vertex(point, gradient(&field, point, reach), middle);
142                    }
143                }
144            }
145        }
146        let cell = |[i, j, k]: [usize; 3]| cell_vertices[(k * cells + j) * cells + i];
147        // Each grid edge the surface crosses joins the four cells around it.
148        for (step, around) in AROUND_EDGES {
149            for k in 0..cells {
150                for j in 0..cells {
151                    for i in 0..cells {
152                        let start = [i, j, k];
153                        let back = |axis: usize| around.iter().any(|offset| offset[axis] > 0);
154                        if (0..3).any(|axis| back(axis) && start[axis] == 0) {
155                            continue;
156                        }
157                        let end = [i + step[0], j + step[1], k + step[2]];
158                        let cells_around =
159                            around.map(|[di, dj, dk]| cell([i - di, j - dj, k - dk]));
160                        net.join(&samples, (start, end), cells_around);
161                    }
162                }
163            }
164        }
165        Self::new(net.positions, net.uvs, Some(net.normals), net.indices)
166    }
167}
168
169/// The vertices and triangles of a surface net as it grows.
170#[derive(Default)]
171struct Net {
172    positions: Vec<[f32; 3]>,
173    normals: Vec<[f32; 3]>,
174    uvs: Vec<[f32; 2]>,
175    indices: Vec<u32>,
176    /// For each vertex, the copy of it whose texture coordinate lies one
177    /// turn on, for the triangles that cross the seam at the back.
178    turned: Vec<u32>,
179}
180
181impl Net {
182    fn vertex(&mut self, point: Vec3, normal: Vec3, middle: Vec3) -> u32 {
183        let index = self.positions.len() as u32;
184        self.positions.push([point.x, point.y, point.z]);
185        self.normals.push([normal.x, normal.y, normal.z]);
186        self.uvs.push(wrapped_uv(point - middle));
187        self.turned.push(u32::MAX);
188        index
189    }
190
191    /// Two triangles over the four cells around the grid edge from `edge.0`
192    /// to `edge.1`, when the surface crosses it.
193    fn join(&mut self, samples: &Samples, edge: ([usize; 3], [usize; 3]), around: [u32; 4]) {
194        let (from, to) = (samples.value(edge.0), samples.value(edge.1));
195        if (from < 0.0) == (to < 0.0) || around.contains(&u32::MAX) {
196            return;
197        }
198        // The quad turns to face the outside, where the field grows.
199        let [a, b, c, d] = if from < 0.0 {
200            around
201        } else {
202            [around[0], around[3], around[2], around[1]]
203        };
204        self.triangle([a, b, c]);
205        self.triangle([a, c, d]);
206    }
207
208    /// A triangle; one that crosses the texture's seam takes the copies of
209    /// its corners on the far side of the seam.
210    fn triangle(&mut self, corners: [u32; 3]) {
211        let us = corners.map(|corner| self.uvs[corner as usize][0]);
212        let crosses = us.iter().copied().fold(f32::NEG_INFINITY, f32::max)
213            - us.iter().copied().fold(f32::INFINITY, f32::min)
214            > 0.5;
215        for (corner, u) in corners.into_iter().zip(us) {
216            let index = if crosses && u < 0.5 {
217                self.turned_copy(corner)
218            } else {
219                corner
220            };
221            self.indices.push(index);
222        }
223    }
224
225    fn turned_copy(&mut self, corner: u32) -> u32 {
226        let at = corner as usize;
227        if self.turned[at] == u32::MAX {
228            self.turned[at] = self.positions.len() as u32;
229            self.positions.push(self.positions[at]);
230            self.normals.push(self.normals[at]);
231            let [u, v] = self.uvs[at];
232            self.uvs.push([u + 1.0, v]);
233            self.turned.push(u32::MAX);
234        }
235        self.turned[at]
236    }
237}
238
239/// The texture coordinate of the point in `direction` from a shape's
240/// middle, as [`Mesh::sphere`] lays its texture.
241fn wrapped_uv(direction: Vec3) -> [f32; 2] {
242    let Some(direction) = direction.normalized() else {
243        return [0.5, 0.5];
244    };
245    [
246        direction.x.atan2(-direction.z) / TAU + 0.5,
247        (-direction.y).clamp(-1.0, 1.0).acos() / PI,
248    ]
249}
250
251/// The direction the field grows fastest at `point`, from steps of `reach`.
252fn gradient(field: &impl Fn(Vec3) -> f32, point: Vec3, reach: f32) -> Vec3 {
253    let along = |axis: Vec3| field(point + axis * reach) - field(point - axis * reach);
254    Vec3::new(
255        along(Vec3::new(1.0, 0.0, 0.0)),
256        along(Vec3::new(0.0, 1.0, 0.0)),
257        along(Vec3::new(0.0, 0.0, 1.0)),
258    )
259    .normalized()
260    .unwrap_or(Vec3::new(0.0, 0.0, -1.0))
261}