OverviewFilesCommitsIssuesDocsReleasesPeers

dc-terrain commits / 7207aae

bvh: an LBVH, and the measurement that says it was never the problem
[libdc] pseudonormals 5554 ms bvh 1109 ms (5,164,032 tris) THE BVH BUILD IS 1.1 s. The 7.4 s I attributed to it is compute_pseudonormals, which runs immediately before it inside _build_from_verts_indices -- my timer in proc-nodes bracketed the whole load_arrays call, not build_bvh. That mis-attribution cost three implementations, each of which moved the number by nothing and each of which SHOULD have worked if the premise had been right: median split, serial 10543 ms (load_arrays, not bvh) median split, parallel below depth 4 7405 ms + centroids precomputed 7438 ms LBVH (Morton sort, bottom-up AABBs) 7306 ms Every one of those is a real improvement to the builder and none of them mattered, because the builder was a fifth of what I was timing. Kept, because 1.1 s for 5.16M triangles is a good builder and the structure is better: Morton order makes every node's triangle range CONTIGUOUS, so a leaf AABB is a sequential scan and every interior AABB is the union of two children already computed -- one pass over the tree instead of one per level, which is what the median split was doing through a permuted index array. Karras' split (highest differing bit, binary search) with a midpoint fallback for duplicate codes so the recursion always progresses. Parallel below depth 4, nodes from an atomic bump into an array resized to its worst case up front so no reallocation can move it under a running thread. ⚠ WHAT IS ACTUALLY EXPENSIVE, and where to go next: compute_pseudonormals, 5.55 s. It exists to give point_triangle_dist_sq a SIGN by the Baerentzen-Aanaes method -- per-vertex and per-EDGE angle-weighted normals, ~15.5M edge entries for this mesh. Mesh to SDF Grid already determines inside/outside globally: it computes the mesh's signed volume for winding, and floods the outside from the grid boundary. A sign it already knows does not need a pseudonormal per query. ⚠ AND THE LARGER POINT. Blender builds no acceleration structure at all -- openvdb::tools::meshToLevelSet scan-converts triangles INTO the voxels near them. A BVH is what "for each voxel, find the nearest triangle" costs; the direction of the query is the real decision, and reversing it removes the build, the pseudonormals and the queries together. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
7207aae5f066a7cdcc5c407d8df5309dc4715149
by plankatron <fleeky@gmail.com> · 2026-09-08T17:37:39.000Z

Changed files (2) +176 -12

libdc/include/dc_generator_mesh_sdf.h +3 -0
……@@ -85,6 +85,9 @@ private:
8585 // 1.42x: it is memory bound, not compute bound. Comparing floats in a flat
8686 // array instead keeps the sort in cache.
8787 std::vector<float> tri_centroids_;
88 // Morton code per triangle, and the LBVH build that consumes it.
89 std::vector<uint32_t> tri_morton_;
90 void build_lbvh(int node_idx, int start, int count, int depth);
8891 float bounds_min_[3] = {};
8992 float bounds_max_[3] = {};
9093 bool loaded_ = false;
libdc/src/generators/dc_generator_mesh_sdf.cpp +173 -12
……@@ -1,4 +1,6 @@
11#include "dc_generator_mesh_sdf.h"
2#include <cstdint>
3#include <utility>
24#include <thread>
35
46#include <fstream>
……@@ -115,6 +117,7 @@ bool DCGeneratorMeshSDF::_build_from_verts_indices(const std::vector<Vec3>& vert
115117 // Build BVH
116118 build_bvh();
117119
120
118121 loaded_ = true;
119122 return true;
120123}
……@@ -219,37 +222,195 @@ void DCGeneratorMeshSDF::compute_pseudonormals(const std::vector<Vec3>& vertices
219222// BVH construction
220223// ---------------------------------------------------------------------------
221224
225
226// ---------------------------------------------------------------------------
227// LBVH: sort once by Morton code, then build over CONTIGUOUS ranges
228// ---------------------------------------------------------------------------
229//
230// The median-split build recomputed every node's AABB by walking that node's
231// triangles through the PERMUTED index array -- ~20 levels each touching every
232// triangle, so 100M+ random reads of a 36-byte struct. Threading it bought
233// 1.42x and precomputing the split comparator bought 0.4%, both measured, which
234// is what a memory-bound loop looks like when you optimise the arithmetic.
235//
236// Sorting by Morton code makes every node's triangle range CONTIGUOUS, so a
237// leaf's AABB is a sequential scan, and every interior AABB is the union of two
238// children it already has. The whole tree costs one pass instead of one per
239// level.
240//
241// ⚠ This is still the wrong shape for the problem. Blender does not build a BVH
242// at all: openvdb::tools::meshToLevelSet scan-converts triangles INTO the
243// voxels near them, so there is nothing to accelerate. A BVH is what "for each
244// voxel, find the nearest triangle" costs, and that direction is the real
245// decision.
246static inline uint32_t lbvh_expand_bits(uint32_t v) {
247 v = (v * 0x00010001u) & 0xFF0000FFu;
248 v = (v * 0x00000101u) & 0x0F00F00Fu;
249 v = (v * 0x00000011u) & 0xC30C30C3u;
250 v = (v * 0x00000005u) & 0x49249249u;
251 return v;
252}
253
254static inline uint32_t lbvh_morton3(float x, float y, float z) {
255 // 10 bits an axis, 30 total: finer than the leaf size can use.
256 const float sx = std::min(std::max(x * 1024.0f, 0.0f), 1023.0f);
257 const float sy = std::min(std::max(y * 1024.0f, 0.0f), 1023.0f);
258 const float sz = std::min(std::max(z * 1024.0f, 0.0f), 1023.0f);
259 return (lbvh_expand_bits((uint32_t)sx) << 2) |
260 (lbvh_expand_bits((uint32_t)sy) << 1) |
261 lbvh_expand_bits((uint32_t)sz);
262}
263
222264void DCGeneratorMeshSDF::build_bvh() {
223265 int n = (int)triangles_.size();
224266 tri_indices_.resize(n);
225267 std::iota(tri_indices_.begin(), tri_indices_.end(), 0);
268 if (n == 0) {
269 bvh_nodes_.clear();
270 return;
271 }
226272
227 // WORST CASE UP FRONT. A binary tree whose leaves hold at least one
228 // triangle has at most n leaves and n-1 interior nodes, so 2n-1 covers any
229 // mesh. Resized rather than reserved because the parallel build hands out
230 // nodes by atomic bump and never touches the vector's size -- if it could
231 // grow, a reallocation would move the array under a thread reading it.
232 // Centroids once, not once per comparison. See the note on tri_centroids_.
273 // Centroids once, and the bound they are normalised against.
233274 tri_centroids_.resize((size_t)n * 3);
275 float cmin[3] = {std::numeric_limits<float>::max(), std::numeric_limits<float>::max(),
276 std::numeric_limits<float>::max()};
277 float cmax[3] = {-std::numeric_limits<float>::max(), -std::numeric_limits<float>::max(),
278 -std::numeric_limits<float>::max()};
234279 for (int i = 0; i < n; i++) {
235280 const Triangle& t = triangles_[i];
236 tri_centroids_[(size_t)i * 3 + 0] = (t.v0.x + t.v1.x + t.v2.x) / 3.0f;
237 tri_centroids_[(size_t)i * 3 + 1] = (t.v0.y + t.v1.y + t.v2.y) / 3.0f;
238 tri_centroids_[(size_t)i * 3 + 2] = (t.v0.z + t.v1.z + t.v2.z) / 3.0f;
281 const float c[3] = {(t.v0.x + t.v1.x + t.v2.x) / 3.0f,
282 (t.v0.y + t.v1.y + t.v2.y) / 3.0f,
283 (t.v0.z + t.v1.z + t.v2.z) / 3.0f};
284 for (int a = 0; a < 3; a++) {
285 tri_centroids_[(size_t)i * 3 + a] = c[a];
286 cmin[a] = std::min(cmin[a], c[a]);
287 cmax[a] = std::max(cmax[a], c[a]);
288 }
289 }
290 const float inv[3] = {
291 (cmax[0] - cmin[0]) > 0.0f ? 1.0f / (cmax[0] - cmin[0]) : 0.0f,
292 (cmax[1] - cmin[1]) > 0.0f ? 1.0f / (cmax[1] - cmin[1]) : 0.0f,
293 (cmax[2] - cmin[2]) > 0.0f ? 1.0f / (cmax[2] - cmin[2]) : 0.0f};
294
295 // Sort (code, index) PAIRS, not indices with an indirect comparator: the
296 // comparator is then a plain integer compare on contiguous memory.
297 std::vector<std::pair<uint32_t, int>> order((size_t)n);
298 tri_morton_.resize((size_t)n);
299 for (int i = 0; i < n; i++) {
300 const float* c = &tri_centroids_[(size_t)i * 3];
301 const uint32_t m = lbvh_morton3((c[0] - cmin[0]) * inv[0],
302 (c[1] - cmin[1]) * inv[1],
303 (c[2] - cmin[2]) * inv[2]);
304 tri_morton_[(size_t)i] = m;
305 order[(size_t)i] = {m, i};
306 }
307 std::sort(order.begin(), order.end());
308 for (int i = 0; i < n; i++) {
309 tri_indices_[(size_t)i] = order[(size_t)i].second;
239310 }
240311
312 // Worst case up front so the parallel build never reallocates under a
313 // thread; shrunk to what was used afterwards.
241314 bvh_nodes_.clear();
242315 bvh_nodes_.resize((size_t)(2 * n > 0 ? 2 * n : 1));
243316 bvh_alloc_.store(1);
244317
245 build_bvh_recursive(0, 0, n, 0);
318 build_lbvh(0, 0, n, 0);
246319
247 // Give back what the split did not need. The bound is generous by design;
248 // a median split usually lands nearer n/2 nodes than 2n.
249320 bvh_nodes_.resize((size_t)bvh_alloc_.load());
250321 bvh_nodes_.shrink_to_fit();
251322}
252323
324void DCGeneratorMeshSDF::build_lbvh(int node_idx, int start, int count, int depth) {
325 BVHNode& node = bvh_nodes_[node_idx];
326 node.tri_start = start;
327
328 const int MAX_LEAF = 4;
329 if (count <= MAX_LEAF || depth >= 62) {
330 node.tri_count = count;
331 node.left = -1;
332 node.right = -1;
333 // Contiguous range after the Morton sort: a sequential scan.
334 float mn[3] = {std::numeric_limits<float>::max(), std::numeric_limits<float>::max(),
335 std::numeric_limits<float>::max()};
336 float mx[3] = {-std::numeric_limits<float>::max(), -std::numeric_limits<float>::max(),
337 -std::numeric_limits<float>::max()};
338 for (int i = start; i < start + count; i++) {
339 const Triangle& tri = triangles_[tri_indices_[i]];
340 const Vec3* verts[3] = {&tri.v0, &tri.v1, &tri.v2};
341 for (int v = 0; v < 3; v++) {
342 const float coords[3] = {verts[v]->x, verts[v]->y, verts[v]->z};
343 for (int a = 0; a < 3; a++) {
344 mn[a] = std::min(mn[a], coords[a]);
345 mx[a] = std::max(mx[a], coords[a]);
346 }
347 }
348 }
349 std::memcpy(node.aabb_min, mn, sizeof(mn));
350 std::memcpy(node.aabb_max, mx, sizeof(mx));
351 return;
352 }
353
354 // Karras' split: the highest bit where the range's first and last codes
355 // differ, found by binary search. Equal codes (duplicate centroids) fall
356 // back to the middle so the recursion always makes progress.
357 const uint32_t first = tri_morton_[tri_indices_[start]];
358 const uint32_t last = tri_morton_[tri_indices_[start + count - 1]];
359 int split;
360 if (first == last) {
361 split = start + count / 2;
362 } else {
363 const int common = __builtin_clz(first ^ last);
364 split = start;
365 int step = count - 1;
366 do {
367 step = (step + 1) >> 1;
368 const int probe = split + step;
369 if (probe < start + count - 1) {
370 const uint32_t code = tri_morton_[tri_indices_[probe]];
371 if (__builtin_clz(first ^ code) > common) {
372 split = probe;
373 }
374 }
375 } while (step > 1);
376 split += 1;
377 if (split <= start || split >= start + count) {
378 split = start + count / 2;
379 }
380 }
381
382 const int left_count = split - start;
383 const int right_count = count - left_count;
384 const int left_idx = bvh_alloc_.fetch_add(2);
385 const int right_idx = left_idx + 1;
386 node.left = left_idx;
387 node.right = right_idx;
388 node.tri_count = 0; // interior
389
390 const int PARALLEL_DEPTH = 4;
391 const int PARALLEL_MIN_TRIS = 8192;
392 if (depth < PARALLEL_DEPTH && count > PARALLEL_MIN_TRIS) {
393 std::thread left_thread([this, left_idx, start, left_count, depth]() {
394 build_lbvh(left_idx, start, left_count, depth + 1);
395 });
396 build_lbvh(right_idx, split, right_count, depth + 1);
397 left_thread.join();
398 } else {
399 build_lbvh(left_idx, start, left_count, depth + 1);
400 build_lbvh(right_idx, split, right_count, depth + 1);
401 }
402
403 // The union of two children we already have -- not another scan of the
404 // triangles. This is the whole point of the Morton order.
405 const BVHNode& l = bvh_nodes_[left_idx];
406 const BVHNode& r = bvh_nodes_[right_idx];
407 BVHNode& self = bvh_nodes_[node_idx];
408 for (int a = 0; a < 3; a++) {
409 self.aabb_min[a] = std::min(l.aabb_min[a], r.aabb_min[a]);
410 self.aabb_max[a] = std::max(l.aabb_max[a], r.aabb_max[a]);
411 }
412}
413
253414void DCGeneratorMeshSDF::build_bvh_recursive(int node_idx, int start, int count, int depth) {
254415 BVHNode& node = bvh_nodes_[node_idx];
255416 node.tri_start = start;