graphic design tool
git clone https://git.lucas.co/cce-designer.git
src/spatial.rs (13.2K)
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 /// Where a ray meets a triangle, as a distance along the ray.
63 ///
64 /// Möller–Trumbore. Lives here beside the other spatial queries because three
65 /// callers want it now: Collision's inside test, and the volume builder's sign
66 /// pass, which casts one ray per grid row.
67 pub fn ray_triangle(origin: Vec3, dir: Vec3, v0: Vec3, v1: Vec3, v2: Vec3) -> Option<f32> {
68 let edge1 = v1 - v0;
69 let edge2 = v2 - v0;
70 let h = dir.cross(edge2);
71 let a = edge1.dot(h);
72 if a.abs() < 1e-6 {
73 return None;
74 }
75 let f = 1.0 / a;
76 let s = origin - v0;
77 let u = f * s.dot(h);
78 if !(0.0..=1.0).contains(&u) {
79 return None;
80 }
81 let q = s.cross(edge1);
82 let v = f * dir.dot(q);
83 if v < 0.0 || u + v > 1.0 {
84 return None;
85 }
86 let t = f * edge2.dot(q);
87 (t > 1e-5).then_some(t)
88 }
89
90 /// Where things are, bucketed by cell.
91 ///
92 /// Shared by both grids: they differ only in what they store and how they
93 /// answer, not in how they divide space.
94 struct Grid {
95 min: Vec3,
96 cell: f32,
97 dims: [i32; 3],
98 buckets: Vec<Vec<u32>>,
99 }
100
101 impl Grid {
102 /// A grid over `bounds` with cells about `cell` across, capped so that a
103 /// pathological request cannot ask for a billion buckets.
104 fn new(min: Vec3, max: Vec3, cell: f32) -> Grid {
105 let span = (max - min).max(Vec3::splat(1e-6));
106 let cell = cell.max(span.max_element() / 128.0).max(1e-6);
107 let dims = [
108 ((span.x / cell).ceil() as i32 + 1).clamp(1, 256),
109 ((span.y / cell).ceil() as i32 + 1).clamp(1, 256),
110 ((span.z / cell).ceil() as i32 + 1).clamp(1, 256),
111 ];
112 let n = (dims[0] * dims[1] * dims[2]) as usize;
113 Grid { min, cell, dims, buckets: vec![Vec::new(); n] }
114 }
115
116 fn coord(&self, p: Vec3) -> [i32; 3] {
117 let rel = (p - self.min) / self.cell;
118 [
119 (rel.x.floor() as i32).clamp(0, self.dims[0] - 1),
120 (rel.y.floor() as i32).clamp(0, self.dims[1] - 1),
121 (rel.z.floor() as i32).clamp(0, self.dims[2] - 1),
122 ]
123 }
124
125 fn index(&self, c: [i32; 3]) -> usize {
126 ((c[2] * self.dims[1] + c[1]) * self.dims[0] + c[0]) as usize
127 }
128
129 fn insert(&mut self, p: Vec3, id: u32) {
130 let i = self.index(self.coord(p));
131 self.buckets[i].push(id);
132 }
133
134 /// Everything in the cells overlapping the box, deduplicated.
135 fn gather(&self, lo: Vec3, hi: Vec3, out: &mut Vec<u32>) {
136 out.clear();
137 let (a, b) = (self.coord(lo), self.coord(hi));
138 for z in a[2]..=b[2] {
139 for y in a[1]..=b[1] {
140 for x in a[0]..=b[0] {
141 out.extend_from_slice(&self.buckets[self.index([x, y, z])]);
142 }
143 }
144 }
145 out.sort_unstable();
146 out.dedup();
147 }
148 }
149
150 /// What a surface lookup found.
151 #[derive(Clone, Copy, Debug)]
152 pub struct Hit {
153 pub point: Vec3,
154 pub distance: f32,
155 /// The face normal of the triangle the hit is on, normalized. Lets a
156 /// caller tell inside from outside in constant time — the sign of
157 /// `(p - point) . normal` — instead of casting a ray through the whole
158 /// mesh. It reads the wrong way in a concave crease, where the nearest
159 /// face is not the one facing you, which is why the node that uses it
160 /// says so.
161 pub normal: Vec3,
162 }
163
164 /// A grid over a surface's triangles, for asking what the nearest surface
165 /// point is.
166 pub struct TriGrid {
167 tris: Vec<[Vec3; 3]>,
168 grid: Grid,
169 }
170
171 impl TriGrid {
172 pub fn build(d: &Detail) -> TriGrid {
173 let tris: Vec<[Vec3; 3]> = d
174 .triangulate(|pos, _| Vec3::from(pos))
175 .chunks_exact(3)
176 .map(|t| [t[0], t[1], t[2]])
177 .collect();
178 let (min, max) = d.bounds().unwrap_or((Vec3::ZERO, Vec3::ZERO));
179 // Cells about the size of a triangle: small enough that a cell holds
180 // few, large enough that one triangle spans few.
181 let mean = if tris.is_empty() {
182 1.0
183 } else {
184 tris.iter()
185 .map(|t| (t[1] - t[0]).length().max((t[2] - t[0]).length()))
186 .sum::<f32>()
187 / tris.len() as f32
188 };
189 let mut grid = Grid::new(min, max, mean.max(1e-5));
190 // A triangle goes in every cell its bounding box touches, so a lookup
191 // that finds a cell finds every triangle passing through it.
192 for (i, t) in tris.iter().enumerate() {
193 let lo = t[0].min(t[1]).min(t[2]);
194 let hi = t[0].max(t[1]).max(t[2]);
195 let (a, b) = (grid.coord(lo), grid.coord(hi));
196 for z in a[2]..=b[2] {
197 for y in a[1]..=b[1] {
198 for x in a[0]..=b[0] {
199 let idx = grid.index([x, y, z]);
200 grid.buckets[idx].push(i as u32);
201 }
202 }
203 }
204 }
205 TriGrid { tris, grid }
206 }
207
208 pub fn is_empty(&self) -> bool {
209 self.tris.is_empty()
210 }
211
212 /// The triangles behind this grid, for a caller that needs them directly —
213 /// the volume builder's scanline sign pass casts rays at all of them.
214 pub fn triangles(&self) -> &[[Vec3; 3]] {
215 &self.tris
216 }
217
218 /// The closest point on the surface, and its distance.
219 ///
220 /// Searches an expanding box until the best hit is closer than the box is
221 /// wide — at which point nothing outside can beat it, because anything out
222 /// there is at least that far away.
223 pub fn closest(&self, p: Vec3) -> Option<Hit> {
224 self.closest_within(p, f32::INFINITY)
225 }
226
227 /// [`TriGrid::closest`], giving up once the search passes `limit`.
228 ///
229 /// The unbounded form doubles its reach until it finds something, so a
230 /// query far from the surface ends up gathering every triangle in the mesh
231 /// and sorting them — which is fine for the handful of queries an operator
232 /// makes and ruinous for the hundred thousand a volume build makes, where
233 /// most samples are nowhere near the surface. A caller that only needs to
234 /// know "further than this" says so and pays for a few cells.
235 pub fn closest_within(&self, p: Vec3, limit: f32) -> Option<Hit> {
236 if self.tris.is_empty() {
237 return None;
238 }
239 if limit.is_finite() {
240 // ONE gather of exactly the box asked for, rather than doubling up
241 // to it: a bounded query knows how far it cares about, and growing
242 // into that size in stages means gathering and sorting the same
243 // cells over and over. This is the difference between a volume
244 // build taking thirty seconds and taking two.
245 let mut scratch = Vec::new();
246 self.grid
247 .gather(p - Vec3::splat(limit), p + Vec3::splat(limit), &mut scratch);
248 return scratch
249 .iter()
250 .map(|&i| {
251 let t = self.tris[i as usize];
252 let q = closest_point_on_triangle(p, t[0], t[1], t[2]);
253 Hit {
254 point: q,
255 distance: (q - p).length(),
256 normal: (t[1] - t[0]).cross(t[2] - t[0]).normalize_or_zero(),
257 }
258 })
259 .min_by(|a, b| {
260 a.distance.partial_cmp(&b.distance).unwrap_or(std::cmp::Ordering::Equal)
261 })
262 .filter(|h| h.distance <= limit);
263 }
264 let hit = |i: usize| {
265 let t = self.tris[i];
266 let q = closest_point_on_triangle(p, t[0], t[1], t[2]);
267 Hit {
268 point: q,
269 distance: (q - p).length(),
270 normal: (t[1] - t[0]).cross(t[2] - t[0]).normalize_or_zero(),
271 }
272 };
273 let nearer = |a: &Hit, b: &Hit| {
274 a.distance.partial_cmp(&b.distance).unwrap_or(std::cmp::Ordering::Equal)
275 };
276
277 let mut reach = self.grid.cell;
278 let mut scratch = Vec::new();
279 for _ in 0..12 {
280 self.grid
281 .gather(p - Vec3::splat(reach), p + Vec3::splat(reach), &mut scratch);
282 let best = scratch.iter().map(|&i| hit(i as usize)).min_by(nearer);
283 match best {
284 Some(h) if h.distance <= reach => return Some(h),
285 _ => reach *= 2.0,
286 }
287 }
288 // The grid is clamped to the geometry's bounds, so twelve doublings
289 // have gathered everything; whatever it found is the answer.
290 (0..self.tris.len()).map(hit).min_by(nearer)
291 }
292 }
293
294 /// A grid over a point set, for asking which points are near a place.
295 pub struct PointGrid {
296 points: Vec<Vec3>,
297 grid: Grid,
298 }
299
300 impl PointGrid {
301 pub fn build(points: &[Vec3], cell: f32) -> PointGrid {
302 let (min, max) = points.iter().fold(
303 (Vec3::splat(f32::MAX), Vec3::splat(f32::MIN)),
304 |(lo, hi), &p| (lo.min(p), hi.max(p)),
305 );
306 let (min, max) = if points.is_empty() { (Vec3::ZERO, Vec3::ZERO) } else { (min, max) };
307 let mut grid = Grid::new(min, max, cell);
308 for (i, &p) in points.iter().enumerate() {
309 grid.insert(p, i as u32);
310 }
311 PointGrid { points: points.to_vec(), grid }
312 }
313
314 pub fn is_empty(&self) -> bool {
315 self.points.is_empty()
316 }
317
318 /// The nearest point, and its distance.
319 ///
320 /// Expands the search box until the best hit is closer than the box is
321 /// wide, the same argument [`TriGrid::closest`] makes: anything outside a
322 /// box that wide is at least that far away, so nothing out there can beat
323 /// what is already in hand.
324 pub fn nearest(&self, p: Vec3) -> Option<(u32, f32)> {
325 if self.points.is_empty() {
326 return None;
327 }
328 let mut reach = self.grid.cell;
329 let mut scratch = Vec::new();
330 for _ in 0..12 {
331 self.grid
332 .gather(p - Vec3::splat(reach), p + Vec3::splat(reach), &mut scratch);
333 let best = scratch
334 .iter()
335 .map(|&i| (i, (self.points[i as usize] - p).length()))
336 .min_by(|a, b| a.1.partial_cmp(&b.1).unwrap_or(std::cmp::Ordering::Equal));
337 match best {
338 Some(hit) if hit.1 <= reach => return Some(hit),
339 _ => reach *= 2.0,
340 }
341 }
342 self.points
343 .iter()
344 .enumerate()
345 .map(|(i, q)| (i as u32, (*q - p).length()))
346 .min_by(|a, b| a.1.partial_cmp(&b.1).unwrap_or(std::cmp::Ordering::Equal))
347 }
348
349 /// Indices of points within `radius` of `p`, excluding nothing — the
350 /// caller decides what does not count as a neighbour, because "not
351 /// itself" and "not topologically adjacent" are different questions.
352 pub fn within(&self, p: Vec3, radius: f32, out: &mut Vec<u32>) {
353 let mut scratch = Vec::new();
354 self.grid
355 .gather(p - Vec3::splat(radius), p + Vec3::splat(radius), &mut scratch);
356 let r2 = radius * radius;
357 out.clear();
358 out.extend(
359 scratch
360 .into_iter()
361 .filter(|&i| (self.points[i as usize] - p).length_squared() <= r2),
362 );
363 }
364 }