graphic design tool
git clone https://git.lucas.co/cce-designer.git
src/spatial.rs (18.6K)
1 //! Uniform grids for the "what is near this?" questions.
2 //!
3 //! Three operators want one of these and all three arrived at once, which is
4 //! why it is built here rather than inside any of them: the remesher's
5 //! projection pass asks for the closest point on a surface, Detangle asks
6 //! which points are within a thickness, and Suture asks both.
7 //!
8 //! A uniform grid rather than a BVH because the geometry these run on is
9 //! already near-uniform — that is what remeshing is for — and a grid sized to
10 //! the mesh's own scale has no degenerate case on it. A surface with wildly
11 //! varying triangle sizes would want a tree, and that is the day to write one.
12
13 use crate::detail::Detail;
14 use glam::Vec3;
15
16 /// The closest point to `p` on triangle `(a, b, c)`.
17 ///
18 /// The Voronoi-region walk from Ericson's *Real-Time Collision Detection*
19 /// §5.1.5: test the three vertex regions, then the three edge regions, and
20 /// what is left is the face interior.
21 pub fn closest_point_on_triangle(p: Vec3, a: Vec3, b: Vec3, c: Vec3) -> Vec3 {
22 let (ab, ac, ap) = (b - a, c - a, p - a);
23 let (d1, d2) = (ab.dot(ap), ac.dot(ap));
24 if d1 <= 0.0 && d2 <= 0.0 {
25 return a;
26 }
27 let bp = p - b;
28 let (d3, d4) = (ab.dot(bp), ac.dot(bp));
29 if d3 >= 0.0 && d4 <= d3 {
30 return b;
31 }
32 let vc = d1 * d4 - d3 * d2;
33 if vc <= 0.0 && d1 >= 0.0 && d3 <= 0.0 {
34 let denom = d1 - d3;
35 let v = if denom.abs() < 1e-20 { 0.0 } else { d1 / denom };
36 return a + ab * v;
37 }
38 let cp = p - c;
39 let (d5, d6) = (ab.dot(cp), ac.dot(cp));
40 if d6 >= 0.0 && d5 <= d6 {
41 return c;
42 }
43 let vb = d5 * d2 - d1 * d6;
44 if vb <= 0.0 && d2 >= 0.0 && d6 <= 0.0 {
45 let denom = d2 - d6;
46 let w = if denom.abs() < 1e-20 { 0.0 } else { d2 / denom };
47 return a + ac * w;
48 }
49 let va = d3 * d6 - d5 * d4;
50 if va <= 0.0 && (d4 - d3) >= 0.0 && (d5 - d6) >= 0.0 {
51 let denom = (d4 - d3) + (d5 - d6);
52 let w = if denom.abs() < 1e-20 { 0.0 } else { (d4 - d3) / denom };
53 return b + (c - b) * w;
54 }
55 let denom = va + vb + vc;
56 if denom.abs() < 1e-20 {
57 return a;
58 }
59 a + ab * (vb / denom) + ac * (vc / denom)
60 }
61
62 /// [`closest_point_on_triangle`] as WEIGHTS: how much of the closest point
63 /// each corner is, summing to one. The point is `a * w[0] + b * w[1] +
64 /// c * w[2]`, and the weights are what lets a caller hand a push on that
65 /// point back to the corners that carry it.
66 ///
67 /// The same Voronoi-region walk, region for region.
68 pub fn closest_weights_on_triangle(p: Vec3, a: Vec3, b: Vec3, c: Vec3) -> [f32; 3] {
69 let (ab, ac, ap) = (b - a, c - a, p - a);
70 let (d1, d2) = (ab.dot(ap), ac.dot(ap));
71 if d1 <= 0.0 && d2 <= 0.0 {
72 return [1.0, 0.0, 0.0];
73 }
74 let bp = p - b;
75 let (d3, d4) = (ab.dot(bp), ac.dot(bp));
76 if d3 >= 0.0 && d4 <= d3 {
77 return [0.0, 1.0, 0.0];
78 }
79 let vc = d1 * d4 - d3 * d2;
80 if vc <= 0.0 && d1 >= 0.0 && d3 <= 0.0 {
81 let denom = d1 - d3;
82 let v = if denom.abs() < 1e-20 { 0.0 } else { d1 / denom };
83 return [1.0 - v, v, 0.0];
84 }
85 let cp = p - c;
86 let (d5, d6) = (ab.dot(cp), ac.dot(cp));
87 if d6 >= 0.0 && d5 <= d6 {
88 return [0.0, 0.0, 1.0];
89 }
90 let vb = d5 * d2 - d1 * d6;
91 if vb <= 0.0 && d2 >= 0.0 && d6 <= 0.0 {
92 let denom = d2 - d6;
93 let w = if denom.abs() < 1e-20 { 0.0 } else { d2 / denom };
94 return [1.0 - w, 0.0, w];
95 }
96 let va = d3 * d6 - d5 * d4;
97 if va <= 0.0 && (d4 - d3) >= 0.0 && (d5 - d6) >= 0.0 {
98 let denom = (d4 - d3) + (d5 - d6);
99 let w = if denom.abs() < 1e-20 { 0.0 } else { (d4 - d3) / denom };
100 return [0.0, 1.0 - w, w];
101 }
102 let denom = va + vb + vc;
103 if denom.abs() < 1e-20 {
104 return [1.0, 0.0, 0.0];
105 }
106 let (v, w) = (vb / denom, vc / denom);
107 [1.0 - v - w, v, w]
108 }
109
110 /// Whether the segment `a`-`b` passes through triangle `(v0, v1, v2)`,
111 /// strictly between its ends.
112 ///
113 /// Möller–Trumbore with the segment as the ray, and a tolerance RELATIVE to
114 /// the lengths involved, so the answer does not change with the model's
115 /// scale. A segment lying in the triangle's plane does not cross it.
116 pub fn segment_crosses_triangle(a: Vec3, b: Vec3, v0: Vec3, v1: Vec3, v2: Vec3) -> bool {
117 let (d, e1, e2) = (b - a, v1 - v0, v2 - v0);
118 let h = d.cross(e2);
119 let det = e1.dot(h);
120 if det.abs() <= 1e-7 * e1.length() * e2.length() * d.length() {
121 return false;
122 }
123 let f = 1.0 / det;
124 let s = a - v0;
125 let u = f * s.dot(h);
126 if !(0.0..=1.0).contains(&u) {
127 return false;
128 }
129 let q = s.cross(e1);
130 let v = f * d.dot(q);
131 if v < 0.0 || u + v > 1.0 {
132 return false;
133 }
134 let t = f * e2.dot(q);
135 t > 0.0 && t < 1.0
136 }
137
138 /// Where a ray meets a triangle, as a distance along the ray.
139 ///
140 /// Möller–Trumbore. Lives here beside the other spatial queries because three
141 /// callers want it now: Collision's inside test, and the volume builder's sign
142 /// pass, which casts one ray per grid row.
143 pub fn ray_triangle(origin: Vec3, dir: Vec3, v0: Vec3, v1: Vec3, v2: Vec3) -> Option<f32> {
144 let edge1 = v1 - v0;
145 let edge2 = v2 - v0;
146 let h = dir.cross(edge2);
147 let a = edge1.dot(h);
148 if a.abs() < 1e-6 {
149 return None;
150 }
151 let f = 1.0 / a;
152 let s = origin - v0;
153 let u = f * s.dot(h);
154 if !(0.0..=1.0).contains(&u) {
155 return None;
156 }
157 let q = s.cross(edge1);
158 let v = f * dir.dot(q);
159 if v < 0.0 || u + v > 1.0 {
160 return None;
161 }
162 let t = f * edge2.dot(q);
163 (t > 1e-5).then_some(t)
164 }
165
166 /// Where things are, bucketed by cell.
167 ///
168 /// Shared by both grids: they differ only in what they store and how they
169 /// answer, not in how they divide space.
170 struct Grid {
171 min: Vec3,
172 cell: f32,
173 dims: [i32; 3],
174 buckets: Vec<Vec<u32>>,
175 }
176
177 impl Grid {
178 /// A grid over `bounds` with cells about `cell` across, capped so that a
179 /// pathological request cannot ask for a billion buckets.
180 fn new(min: Vec3, max: Vec3, cell: f32) -> Grid {
181 let span = (max - min).max(Vec3::splat(1e-6));
182 let cell = cell.max(span.max_element() / 128.0).max(1e-6);
183 let dims = [
184 ((span.x / cell).ceil() as i32 + 1).clamp(1, 256),
185 ((span.y / cell).ceil() as i32 + 1).clamp(1, 256),
186 ((span.z / cell).ceil() as i32 + 1).clamp(1, 256),
187 ];
188 let n = (dims[0] * dims[1] * dims[2]) as usize;
189 Grid { min, cell, dims, buckets: vec![Vec::new(); n] }
190 }
191
192 fn coord(&self, p: Vec3) -> [i32; 3] {
193 let rel = (p - self.min) / self.cell;
194 [
195 (rel.x.floor() as i32).clamp(0, self.dims[0] - 1),
196 (rel.y.floor() as i32).clamp(0, self.dims[1] - 1),
197 (rel.z.floor() as i32).clamp(0, self.dims[2] - 1),
198 ]
199 }
200
201 fn index(&self, c: [i32; 3]) -> usize {
202 ((c[2] * self.dims[1] + c[1]) * self.dims[0] + c[0]) as usize
203 }
204
205 fn insert(&mut self, p: Vec3, id: u32) {
206 let i = self.index(self.coord(p));
207 self.buckets[i].push(id);
208 }
209
210 /// Everything in the cells overlapping the box, deduplicated.
211 fn gather(&self, lo: Vec3, hi: Vec3, out: &mut Vec<u32>) {
212 out.clear();
213 let (a, b) = (self.coord(lo), self.coord(hi));
214 for z in a[2]..=b[2] {
215 for y in a[1]..=b[1] {
216 for x in a[0]..=b[0] {
217 out.extend_from_slice(&self.buckets[self.index([x, y, z])]);
218 }
219 }
220 }
221 out.sort_unstable();
222 out.dedup();
223 }
224 }
225
226 /// What a surface lookup found.
227 #[derive(Clone, Copy, Debug)]
228 pub struct Hit {
229 pub point: Vec3,
230 pub distance: f32,
231 /// The face normal of the triangle the hit is on, normalized. Lets a
232 /// caller tell inside from outside in constant time — the sign of
233 /// `(p - point) . normal` — instead of casting a ray through the whole
234 /// mesh. It reads the wrong way in a concave crease, where the nearest
235 /// face is not the one facing you, which is why the node that uses it
236 /// says so.
237 pub normal: Vec3,
238 }
239
240 /// A grid over a surface's triangles, for asking what the nearest surface
241 /// point is.
242 pub struct TriGrid {
243 tris: Vec<[Vec3; 3]>,
244 grid: Grid,
245 }
246
247 impl TriGrid {
248 pub fn build(d: &Detail) -> TriGrid {
249 let tris: Vec<[Vec3; 3]> = d
250 .triangulate(|pos, _| Vec3::from(pos))
251 .chunks_exact(3)
252 .map(|t| [t[0], t[1], t[2]])
253 .collect();
254 let (min, max) = d.bounds().unwrap_or((Vec3::ZERO, Vec3::ZERO));
255 // Cells about the size of a triangle: small enough that a cell holds
256 // few, large enough that one triangle spans few.
257 let mean = if tris.is_empty() {
258 1.0
259 } else {
260 tris.iter()
261 .map(|t| (t[1] - t[0]).length().max((t[2] - t[0]).length()))
262 .sum::<f32>()
263 / tris.len() as f32
264 };
265 let mut grid = Grid::new(min, max, mean.max(1e-5));
266 // A triangle goes in every cell its bounding box touches, so a lookup
267 // that finds a cell finds every triangle passing through it.
268 for (i, t) in tris.iter().enumerate() {
269 let lo = t[0].min(t[1]).min(t[2]);
270 let hi = t[0].max(t[1]).max(t[2]);
271 let (a, b) = (grid.coord(lo), grid.coord(hi));
272 for z in a[2]..=b[2] {
273 for y in a[1]..=b[1] {
274 for x in a[0]..=b[0] {
275 let idx = grid.index([x, y, z]);
276 grid.buckets[idx].push(i as u32);
277 }
278 }
279 }
280 }
281 TriGrid { tris, grid }
282 }
283
284 pub fn is_empty(&self) -> bool {
285 self.tris.is_empty()
286 }
287
288 /// The triangles behind this grid, for a caller that needs them directly —
289 /// the volume builder's scanline sign pass casts rays at all of them.
290 pub fn triangles(&self) -> &[[Vec3; 3]] {
291 &self.tris
292 }
293
294 /// The closest point on the surface, and its distance.
295 ///
296 /// Searches an expanding box until the best hit is closer than the box is
297 /// wide — at which point nothing outside can beat it, because anything out
298 /// there is at least that far away.
299 pub fn closest(&self, p: Vec3) -> Option<Hit> {
300 self.closest_within(p, f32::INFINITY)
301 }
302
303 /// [`TriGrid::closest`], giving up once the search passes `limit`.
304 ///
305 /// The unbounded form doubles its reach until it finds something, so a
306 /// query far from the surface ends up gathering every triangle in the mesh
307 /// and sorting them — which is fine for the handful of queries an operator
308 /// makes and ruinous for the hundred thousand a volume build makes, where
309 /// most samples are nowhere near the surface. A caller that only needs to
310 /// know "further than this" says so and pays for a few cells.
311 pub fn closest_within(&self, p: Vec3, limit: f32) -> Option<Hit> {
312 if self.tris.is_empty() {
313 return None;
314 }
315 if limit.is_finite() {
316 // ONE gather of exactly the box asked for, rather than doubling up
317 // to it: a bounded query knows how far it cares about, and growing
318 // into that size in stages means gathering and sorting the same
319 // cells over and over. This is the difference between a volume
320 // build taking thirty seconds and taking two.
321 let mut scratch = Vec::new();
322 self.grid
323 .gather(p - Vec3::splat(limit), p + Vec3::splat(limit), &mut scratch);
324 return scratch
325 .iter()
326 .map(|&i| {
327 let t = self.tris[i as usize];
328 let q = closest_point_on_triangle(p, t[0], t[1], t[2]);
329 Hit {
330 point: q,
331 distance: (q - p).length(),
332 normal: (t[1] - t[0]).cross(t[2] - t[0]).normalize_or_zero(),
333 }
334 })
335 .min_by(|a, b| {
336 a.distance.partial_cmp(&b.distance).unwrap_or(std::cmp::Ordering::Equal)
337 })
338 .filter(|h| h.distance <= limit);
339 }
340 // The nearest of the triangles named, and of those equally near the
341 // one with the lowest index — which is what the first version's
342 // `min_by` over a sorted list chose. Only the winner's normal is
343 // worked out: it was a cross product and a square root for every
344 // candidate, and a candidate is a triangle for every cell it is
345 // filed in.
346 let nearest = |among: &mut dyn Iterator<Item = u32>| -> Option<(f32, u32, Vec3)> {
347 let mut best: Option<(f32, u32, Vec3)> = None;
348 for i in among {
349 if best.is_some_and(|b| b.1 == i) {
350 continue;
351 }
352 let t = self.tris[i as usize];
353 let q = closest_point_on_triangle(p, t[0], t[1], t[2]);
354 let d = (q - p).length();
355 let better = match best {
356 None => true,
357 Some((bd, bi, _)) => match d.partial_cmp(&bd) {
358 Some(std::cmp::Ordering::Less) => true,
359 Some(std::cmp::Ordering::Greater) => false,
360 _ => i < bi,
361 },
362 };
363 if better {
364 best = Some((d, i, q));
365 }
366 }
367 best
368 };
369 let hit = |(distance, i, point): (f32, u32, Vec3)| {
370 let t = self.tris[i as usize];
371 Hit { point, distance, normal: (t[1] - t[0]).cross(t[2] - t[0]).normalize_or_zero() }
372 };
373
374 // From a quarter of a cell, where the first version began at a
375 // whole one: a query ON the surface — a remesh projecting its
376 // points, which is most of what is asked — is answered from the
377 // cell it is in and not from the twenty-seven about it. Where the
378 // search begins does not change what it finds: it ends only with a
379 // hit nearer than the box searched is wide, which nothing outside
380 // the box can beat or tie.
381 let mut reach = self.grid.cell * 0.25;
382 for _ in 0..14 {
383 let (a, b) = (self.grid.coord(p - Vec3::splat(reach)), self.grid.coord(p + Vec3::splat(reach)));
384 let mut among = (a[2]..=b[2]).flat_map(|z| (a[1]..=b[1]).map(move |y| (y, z))).flat_map(|(y, z)| {
385 (a[0]..=b[0]).flat_map(move |x| self.grid.buckets[self.grid.index([x, y, z])].iter().copied())
386 });
387 match nearest(&mut among) {
388 Some(h) if h.0 <= reach => return Some(hit(h)),
389 _ => reach *= 2.0,
390 }
391 }
392 // The grid is clamped to the geometry's bounds, so the doublings
393 // have gathered everything; whatever it found is the answer.
394 nearest(&mut (0..self.tris.len() as u32)).map(hit)
395 }
396
397 /// [`TriGrid::closest`] as it was first written, kept to hold the
398 /// faster one to: the same hit, bit for bit.
399 #[cfg(test)]
400 pub fn closest_reference(&self, p: Vec3) -> Option<Hit> {
401 if self.tris.is_empty() {
402 return None;
403 }
404 let hit = |i: usize| {
405 let t = self.tris[i];
406 let q = closest_point_on_triangle(p, t[0], t[1], t[2]);
407 Hit {
408 point: q,
409 distance: (q - p).length(),
410 normal: (t[1] - t[0]).cross(t[2] - t[0]).normalize_or_zero(),
411 }
412 };
413 let nearer = |a: &Hit, b: &Hit| {
414 a.distance.partial_cmp(&b.distance).unwrap_or(std::cmp::Ordering::Equal)
415 };
416 let mut reach = self.grid.cell;
417 let mut scratch = Vec::new();
418 for _ in 0..12 {
419 self.grid
420 .gather(p - Vec3::splat(reach), p + Vec3::splat(reach), &mut scratch);
421 let best = scratch.iter().map(|&i| hit(i as usize)).min_by(nearer);
422 match best {
423 Some(h) if h.distance <= reach => return Some(h),
424 _ => reach *= 2.0,
425 }
426 }
427 (0..self.tris.len()).map(hit).min_by(nearer)
428 }
429 }
430
431 /// A grid over a point set, for asking which points are near a place.
432 pub struct PointGrid {
433 points: Vec<Vec3>,
434 grid: Grid,
435 }
436
437 impl PointGrid {
438 pub fn build(points: &[Vec3], cell: f32) -> PointGrid {
439 let (min, max) = points.iter().fold(
440 (Vec3::splat(f32::MAX), Vec3::splat(f32::MIN)),
441 |(lo, hi), &p| (lo.min(p), hi.max(p)),
442 );
443 let (min, max) = if points.is_empty() { (Vec3::ZERO, Vec3::ZERO) } else { (min, max) };
444 let mut grid = Grid::new(min, max, cell);
445 for (i, &p) in points.iter().enumerate() {
446 grid.insert(p, i as u32);
447 }
448 PointGrid { points: points.to_vec(), grid }
449 }
450
451 pub fn is_empty(&self) -> bool {
452 self.points.is_empty()
453 }
454
455 /// The nearest point, and its distance.
456 ///
457 /// Expands the search box until the best hit is closer than the box is
458 /// wide, the same argument [`TriGrid::closest`] makes: anything outside a
459 /// box that wide is at least that far away, so nothing out there can beat
460 /// what is already in hand.
461 pub fn nearest(&self, p: Vec3) -> Option<(u32, f32)> {
462 if self.points.is_empty() {
463 return None;
464 }
465 let mut reach = self.grid.cell;
466 let mut scratch = Vec::new();
467 for _ in 0..12 {
468 self.grid
469 .gather(p - Vec3::splat(reach), p + Vec3::splat(reach), &mut scratch);
470 let best = scratch
471 .iter()
472 .map(|&i| (i, (self.points[i as usize] - p).length()))
473 .min_by(|a, b| a.1.partial_cmp(&b.1).unwrap_or(std::cmp::Ordering::Equal));
474 match best {
475 Some(hit) if hit.1 <= reach => return Some(hit),
476 _ => reach *= 2.0,
477 }
478 }
479 self.points
480 .iter()
481 .enumerate()
482 .map(|(i, q)| (i as u32, (*q - p).length()))
483 .min_by(|a, b| a.1.partial_cmp(&b.1).unwrap_or(std::cmp::Ordering::Equal))
484 }
485
486 /// Indices of points within `radius` of `p`, excluding nothing — the
487 /// caller decides what does not count as a neighbour, because "not
488 /// itself" and "not topologically adjacent" are different questions.
489 pub fn within(&self, p: Vec3, radius: f32, out: &mut Vec<u32>) {
490 let mut scratch = Vec::new();
491 self.grid
492 .gather(p - Vec3::splat(radius), p + Vec3::splat(radius), &mut scratch);
493 let r2 = radius * radius;
494 out.clear();
495 out.extend(
496 scratch
497 .into_iter()
498 .filter(|&i| (self.points[i as usize] - p).length_squared() <= r2),
499 );
500 }
501 }