1use std::f32::consts::{PI, TAU};
5
6use super::{Mesh, MeshError, Vec3};
7
8const 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
20const 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
36const 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
44struct 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 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 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 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#[derive(Default)]
171struct Net {
172 positions: Vec<[f32; 3]>,
173 normals: Vec<[f32; 3]>,
174 uvs: Vec<[f32; 2]>,
175 indices: Vec<u32>,
176 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 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 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 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
239fn 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
251fn 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}