1use std::collections::{BTreeMap, HashMap};
2
3use crate::math::{add, dot, length, normalize, scale, sub};
4use crate::mesh::Mesh;
5
6const BLOCK: i32 = 8;
7const SIDE: usize = BLOCK as usize + 1;
8const REACH: f32 = 1.5;
9
10#[derive(Clone, Copy, Debug, PartialEq)]
11pub struct Surface {
12 pub voxel: f32,
13 pub relax: u32,
14}
15
16impl Surface {
17 pub fn new(voxel: f32, relax: u32) -> Self {
18 Self { voxel, relax }
19 }
20
21 pub fn mesh<F: Fn([f32; 3]) -> f32>(&self, field: &F, lo: [f32; 3], hi: [f32; 3]) -> Mesh {
22 mesh(field, lo, hi, self.voxel, self.relax)
23 }
24}
25
26pub fn smooth_min(a: f32, b: f32, blend: f32) -> f32 {
27 if blend <= 0.0 {
28 return a.min(b);
29 }
30 let h = (blend - (a - b).abs()).max(0.0) / blend;
31 a.min(b) - h * h * blend * 0.25
32}
33
34pub fn gradient<F: Fn([f32; 3]) -> f32>(field: &F, p: [f32; 3], eps: f32) -> [f32; 3] {
35 let axis = |i: usize| {
36 let mut plus = p;
37 let mut minus = p;
38 plus[i] += eps;
39 minus[i] -= eps;
40 (field(plus) - field(minus)) / (2.0 * eps)
41 };
42 [axis(0), axis(1), axis(2)]
43}
44
45struct Grid {
46 origin: [f32; 3],
47 voxel: f32,
48 blocks: BTreeMap<[i32; 3], Vec<f32>>,
49}
50
51impl Grid {
52 fn point(&self, cell: [i32; 3]) -> [f32; 3] {
53 [
54 self.origin[0] + cell[0] as f32 * self.voxel,
55 self.origin[1] + cell[1] as f32 * self.voxel,
56 self.origin[2] + cell[2] as f32 * self.voxel,
57 ]
58 }
59
60 fn value(&self, block: &[f32], local: [i32; 3]) -> f32 {
61 block[(local[2] as usize * SIDE + local[1] as usize) * SIDE + local[0] as usize]
62 }
63}
64
65pub fn mesh<F: Fn([f32; 3]) -> f32>(
66 field: &F,
67 lo: [f32; 3],
68 hi: [f32; 3],
69 voxel: f32,
70 relax: u32,
71) -> Mesh {
72 let voxel = voxel.max(1e-6);
73 let pad = 2.0 * voxel;
74 let origin = sub(lo, [pad; 3]);
75 let counts = [0, 1, 2].map(|i| (((hi[i] - lo[i] + 2.0 * pad) / voxel).ceil() as i32).max(1));
76 let block_counts = counts.map(|c| (c + BLOCK - 1) / BLOCK);
77 let reach = REACH * 3f32.sqrt() * (BLOCK as f32 * 0.5 + 1.0) * voxel;
78 let mut grid = Grid {
79 origin,
80 voxel,
81 blocks: BTreeMap::new(),
82 };
83 for bz in 0..block_counts[2] {
84 for by in 0..block_counts[1] {
85 for bx in 0..block_counts[0] {
86 let base = [bx * BLOCK, by * BLOCK, bz * BLOCK];
87 let center = grid.point(base.map(|v| v + BLOCK / 2));
88 if field(center).abs() > reach {
89 continue;
90 }
91 let mut values = Vec::with_capacity(SIDE * SIDE * SIDE);
92 for z in 0..SIDE as i32 {
93 for y in 0..SIDE as i32 {
94 for x in 0..SIDE as i32 {
95 values.push(field(grid.point([base[0] + x, base[1] + y, base[2] + z])));
96 }
97 }
98 }
99 grid.blocks.insert([bx, by, bz], values);
100 }
101 }
102 }
103 let mut positions: Vec<[f32; 3]> = Vec::new();
104 let mut cells: Vec<([i32; 3], [i32; 3])> = Vec::new();
105 let mut index: HashMap<[i32; 3], u32> = HashMap::new();
106 for (key, values) in &grid.blocks {
107 for z in 0..BLOCK {
108 for y in 0..BLOCK {
109 for x in 0..BLOCK {
110 let local = [x, y, z];
111 let mut corners = [0.0f32; 8];
112 let mut inside = 0;
113 for (corner, slot) in corners.iter_mut().enumerate() {
114 let offset = [
115 (corner & 1) as i32,
116 ((corner >> 1) & 1) as i32,
117 (corner >> 2) as i32,
118 ];
119 *slot = grid.value(values, [x + offset[0], y + offset[1], z + offset[2]]);
120 if *slot < 0.0 {
121 inside += 1;
122 }
123 }
124 if inside == 0 || inside == 8 {
125 continue;
126 }
127 let cell = [key[0] * BLOCK + x, key[1] * BLOCK + y, key[2] * BLOCK + z];
128 let corner_at = |c: usize| -> [f32; 3] {
129 [(c & 1) as f32, ((c >> 1) & 1) as f32, (c >> 2) as f32]
130 };
131 let mut sum = [0.0f32; 3];
132 let mut count = 0.0f32;
133 for (a, b) in EDGES {
134 let (va, vb) = (corners[a], corners[b]);
135 if (va < 0.0) == (vb < 0.0) {
136 continue;
137 }
138 let t = va / (va - vb);
139 let pa = corner_at(a);
140 let pb = corner_at(b);
141 sum = add(sum, add(pa, scale(sub(pb, pa), t)));
142 count += 1.0;
143 }
144 let unit = scale(sum, 1.0 / count);
145 let position = add(grid.point(cell), scale(unit, voxel));
146 index.insert(cell, positions.len() as u32);
147 positions.push(position);
148 cells.push((*key, local));
149 }
150 }
151 }
152 }
153 let mut quads: Vec<[u32; 4]> = Vec::new();
154 for (key, local) in &cells {
155 let values = &grid.blocks[key];
156 let cell = [
157 key[0] * BLOCK + local[0],
158 key[1] * BLOCK + local[1],
159 key[2] * BLOCK + local[2],
160 ];
161 let start = grid.value(values, *local);
162 for axis in 0..3 {
163 let (b, c) = ((axis + 1) % 3, (axis + 2) % 3);
164 let mut end_local = *local;
165 end_local[axis] += 1;
166 let end = grid.value(values, end_local);
167 if (start < 0.0) == (end < 0.0) {
168 continue;
169 }
170 let step = |cell: [i32; 3], along: usize| {
171 let mut out = cell;
172 out[along] -= 1;
173 out
174 };
175 let around = [cell, step(cell, b), step(step(cell, b), c), step(cell, c)];
176 let found = around.map(|cell| index.get(&cell).copied());
177 if let [Some(q0), Some(q1), Some(q2), Some(q3)] = found {
178 if start < 0.0 {
179 quads.push([q0, q1, q2, q3]);
180 } else {
181 quads.push([q0, q3, q2, q1]);
182 }
183 }
184 }
185 }
186 let eps = voxel * 0.1;
187 let mut anchors = positions.clone();
188 for (p, anchor) in positions.iter_mut().zip(&anchors) {
189 *p = project(field, *p, *anchor, voxel, eps, 3);
190 }
191 if relax > 0 {
192 let mut neighbours: Vec<Vec<u32>> = vec![Vec::new(); positions.len()];
193 for quad in &quads {
194 for i in 0..4 {
195 let (a, b) = (quad[i], quad[(i + 1) % 4]);
196 neighbours[a as usize].push(b);
197 neighbours[b as usize].push(a);
198 }
199 }
200 for list in &mut neighbours {
201 list.sort_unstable();
202 list.dedup();
203 }
204 anchors.clone_from(&positions);
205 for _ in 0..relax {
206 let previous = positions.clone();
207 for (i, list) in neighbours.iter().enumerate() {
208 if list.is_empty() {
209 continue;
210 }
211 let mut sum = [0.0f32; 3];
212 for &n in list {
213 sum = add(sum, previous[n as usize]);
214 }
215 let average = scale(sum, 1.0 / list.len() as f32);
216 let moved = add(previous[i], scale(sub(average, previous[i]), 0.5));
217 positions[i] = project(field, moved, anchors[i], voxel, eps, 1);
218 }
219 }
220 }
221 let normals = positions
222 .iter()
223 .map(|p| normalize(gradient(field, *p, eps)))
224 .collect::<Vec<_>>();
225 let mut indices = Vec::with_capacity(quads.len() * 6);
226 for q in &quads {
227 let p = q.map(|i| positions[i as usize]);
228 if length(sub(p[0], p[2])) <= length(sub(p[1], p[3])) {
229 indices.extend([q[0], q[1], q[2], q[0], q[2], q[3]]);
230 } else {
231 indices.extend([q[0], q[1], q[3], q[1], q[2], q[3]]);
232 }
233 }
234 let uvs = planar_uvs(&positions);
235 Mesh::new(positions, normals, Vec::new(), uvs, indices)
236}
237
238const EDGES: [(usize, usize); 12] = [
239 (0, 1),
240 (2, 3),
241 (4, 5),
242 (6, 7),
243 (0, 2),
244 (1, 3),
245 (4, 6),
246 (5, 7),
247 (0, 4),
248 (1, 5),
249 (2, 6),
250 (3, 7),
251];
252
253fn project<F: Fn([f32; 3]) -> f32>(
254 field: &F,
255 mut p: [f32; 3],
256 anchor: [f32; 3],
257 voxel: f32,
258 eps: f32,
259 steps: u32,
260) -> [f32; 3] {
261 for _ in 0..steps {
262 let d = field(p);
263 let g = gradient(field, p, eps);
264 let g2 = dot(g, g);
265 if g2 <= 1e-12 {
266 break;
267 }
268 p = sub(p, scale(g, d / g2));
269 let offset = sub(p, anchor);
270 let far = length(offset);
271 if far > voxel {
272 p = add(anchor, scale(offset, voxel / far));
273 }
274 }
275 p
276}
277
278pub(crate) fn planar_uvs(positions: &[[f32; 3]]) -> Vec<[f32; 2]> {
279 let mut lo = [f32::INFINITY; 2];
280 let mut hi = [f32::NEG_INFINITY; 2];
281 for p in positions {
282 for i in 0..2 {
283 lo[i] = lo[i].min(p[i]);
284 hi[i] = hi[i].max(p[i]);
285 }
286 }
287 let span = (hi[0] - lo[0]).max(hi[1] - lo[1]).max(1e-9);
288 positions
289 .iter()
290 .map(|p| [(p[0] - lo[0]) / span, (hi[1] - p[1]) / span])
291 .collect()
292}