#include "dc_generator_mesh_sdf.h"
#include <cstdint>
#include <utility>
#include <thread>
#include <fstream>
#include <sstream>
#include <cstring>
#include <limits>
#include <numeric>
// ---------------------------------------------------------------------------
// OBJ loading
// ---------------------------------------------------------------------------
bool DCGeneratorMeshSDF::load_obj(const std::string& path) {
fprintf(stderr, "[MESH-SDF] load_obj: opening %s\n", path.c_str());
std::ifstream file(path);
fprintf(stderr, "[MESH-SDF] file opened: %d\n", file.is_open());
if (!file.is_open()) return false;
std::vector<Vec3> vertices;
std::vector<int> face_indices; // flat triplets of vertex indices
std::string line;
while (std::getline(file, line)) {
if (line.empty() || line[0] == '#') continue;
std::istringstream iss(line);
std::string token;
iss >> token;
if (token == "v") {
Vec3 v;
iss >> v.x >> v.y >> v.z;
vertices.push_back(v);
} else if (token == "f") {
// Support: "f 1 2 3", "f 1//1 2//2 3//3", "f 1/1 2/2 3/3", "f 1/1/1 ..."
std::vector<int> face_verts;
std::string vert_token;
while (iss >> vert_token) {
int vi = 0;
size_t slash = vert_token.find('/');
if (slash != std::string::npos) {
vi = std::stoi(vert_token.substr(0, slash));
} else {
vi = std::stoi(vert_token);
}
face_verts.push_back(vi - 1); // OBJ is 1-based
}
// Triangulate polygon fans
for (size_t i = 2; i < face_verts.size(); i++) {
face_indices.push_back(face_verts[0]);
face_indices.push_back(face_verts[i - 1]);
face_indices.push_back(face_verts[i]);
}
}
}
if (vertices.empty() || face_indices.empty()) return false;
return _build_from_verts_indices(vertices, face_indices);
}
bool DCGeneratorMeshSDF::load_arrays(const float* positions, int num_verts,
const int* indices, int num_tris) {
if (!positions || !indices || num_verts < 3 || num_tris < 1) return false;
std::vector<Vec3> vertices(num_verts);
for (int i = 0; i < num_verts; i++) {
vertices[i] = {positions[i * 3], positions[i * 3 + 1], positions[i * 3 + 2]};
}
std::vector<int> face_indices(num_tris * 3);
for (int i = 0; i < num_tris * 3; i++) {
face_indices[i] = indices[i];
}
return _build_from_verts_indices(vertices, face_indices);
}
bool DCGeneratorMeshSDF::_build_from_verts_indices(const std::vector<Vec3>& vertices,
const std::vector<int>& face_indices) {
// Build triangles with face normals
int num_tris = (int)face_indices.size() / 3;
triangles_.resize(num_tris);
for (int i = 0; i < num_tris; i++) {
int i0 = face_indices[i * 3 + 0];
int i1 = face_indices[i * 3 + 1];
int i2 = face_indices[i * 3 + 2];
triangles_[i].v0 = vertices[i0];
triangles_[i].v1 = vertices[i1];
triangles_[i].v2 = vertices[i2];
Vec3 e1 = triangles_[i].v1 - triangles_[i].v0;
Vec3 e2 = triangles_[i].v2 - triangles_[i].v0;
triangles_[i].normal = e1.cross(e2).normalized();
}
// Compute proper pseudonormals for robust sign determination
if (!unsigned_only_) {
compute_pseudonormals(vertices, face_indices);
}
// Compute bounds
bounds_min_[0] = bounds_min_[1] = bounds_min_[2] = std::numeric_limits<float>::max();
bounds_max_[0] = bounds_max_[1] = bounds_max_[2] = -std::numeric_limits<float>::max();
for (const auto& v : vertices) {
bounds_min_[0] = std::min(bounds_min_[0], v.x);
bounds_min_[1] = std::min(bounds_min_[1], v.y);
bounds_min_[2] = std::min(bounds_min_[2], v.z);
bounds_max_[0] = std::max(bounds_max_[0], v.x);
bounds_max_[1] = std::max(bounds_max_[1], v.y);
bounds_max_[2] = std::max(bounds_max_[2], v.z);
}
// Build BVH
build_bvh();
loaded_ = true;
return true;
}
// ---------------------------------------------------------------------------
// Pseudonormal computation (Baerentzen & Aanaes 2005)
// ---------------------------------------------------------------------------
void DCGeneratorMeshSDF::compute_pseudonormals(const std::vector<Vec3>& vertices,
const std::vector<int>& face_indices) {
int num_tris = (int)triangles_.size();
int num_verts = (int)vertices.size();
// Vertex pseudonormals: accumulate face normals weighted by angle at vertex
std::vector<Vec3> vert_normals(num_verts);
for (int t = 0; t < num_tris; t++) {
int i0 = face_indices[t * 3 + 0];
int i1 = face_indices[t * 3 + 1];
int i2 = face_indices[t * 3 + 2];
const Vec3& fn = triangles_[t].normal;
Vec3 e01 = vertices[i1] - vertices[i0];
Vec3 e02 = vertices[i2] - vertices[i0];
Vec3 e10 = vertices[i0] - vertices[i1];
Vec3 e12 = vertices[i2] - vertices[i1];
Vec3 e20 = vertices[i0] - vertices[i2];
Vec3 e21 = vertices[i1] - vertices[i2];
float l01 = e01.length(), l02 = e02.length();
float l10 = e10.length(), l12 = e12.length();
float l20 = e20.length(), l21 = e21.length();
auto safe_acos = [](float x) -> float {
return std::acos(std::max(-1.0f, std::min(1.0f, x)));
};
float angle0 = (l01 > 0 && l02 > 0) ? safe_acos(e01.dot(e02) / (l01 * l02)) : 0.0f;
float angle1 = (l10 > 0 && l12 > 0) ? safe_acos(e10.dot(e12) / (l10 * l12)) : 0.0f;
float angle2 = (l20 > 0 && l21 > 0) ? safe_acos(e20.dot(e21) / (l20 * l21)) : 0.0f;
vert_normals[i0] = vert_normals[i0] + fn * angle0;
vert_normals[i1] = vert_normals[i1] + fn * angle1;
vert_normals[i2] = vert_normals[i2] + fn * angle2;
}
for (int i = 0; i < num_verts; i++) {
vert_normals[i] = vert_normals[i].normalized();
}
// Edge pseudonormals: average face normals of the two adjacent triangles
struct EdgeKey {
int v0, v1;
bool operator<(const EdgeKey& o) const {
return v0 < o.v0 || (v0 == o.v0 && v1 < o.v1);
}
};
std::map<EdgeKey, Vec3> edge_normal_map;
for (int t = 0; t < num_tris; t++) {
int idx[3] = {face_indices[t * 3 + 0], face_indices[t * 3 + 1], face_indices[t * 3 + 2]};
const Vec3& fn = triangles_[t].normal;
for (int e = 0; e < 3; e++) {
int ea = idx[e], eb = idx[(e + 1) % 3];
EdgeKey key = {std::min(ea, eb), std::max(ea, eb)};
auto it = edge_normal_map.find(key);
if (it == edge_normal_map.end()) {
edge_normal_map[key] = fn;
} else {
it->second = it->second + fn;
}
}
}
for (auto& pair : edge_normal_map) {
pair.second = pair.second.normalized();
}
// Assign pseudonormals to each triangle
for (int t = 0; t < num_tris; t++) {
int i0 = face_indices[t * 3 + 0];
int i1 = face_indices[t * 3 + 1];
int i2 = face_indices[t * 3 + 2];
EdgeKey k01 = {std::min(i0, i1), std::max(i0, i1)};
EdgeKey k12 = {std::min(i1, i2), std::max(i1, i2)};
EdgeKey k20 = {std::min(i2, i0), std::max(i2, i0)};
triangles_[t].edge_normals[0] = edge_normal_map[k01];
triangles_[t].edge_normals[1] = edge_normal_map[k12];
triangles_[t].edge_normals[2] = edge_normal_map[k20];
triangles_[t].vert_normals[0] = vert_normals[i0];
triangles_[t].vert_normals[1] = vert_normals[i1];
triangles_[t].vert_normals[2] = vert_normals[i2];
}
}
// ---------------------------------------------------------------------------
// BVH construction
// ---------------------------------------------------------------------------
// ---------------------------------------------------------------------------
// LBVH: sort once by Morton code, then build over CONTIGUOUS ranges
// ---------------------------------------------------------------------------
//
// The median-split build recomputed every node's AABB by walking that node's
// triangles through the PERMUTED index array -- ~20 levels each touching every
// triangle, so 100M+ random reads of a 36-byte struct. Threading it bought
// 1.42x and precomputing the split comparator bought 0.4%, both measured, which
// is what a memory-bound loop looks like when you optimise the arithmetic.
//
// Sorting by Morton code makes every node's triangle range CONTIGUOUS, so a
// leaf's AABB is a sequential scan, and every interior AABB is the union of two
// children it already has. The whole tree costs one pass instead of one per
// level.
//
// ⚠ This is still the wrong shape for the problem. Blender does not build a BVH
// at all: openvdb::tools::meshToLevelSet scan-converts triangles INTO the
// voxels near them, so there is nothing to accelerate. A BVH is what "for each
// voxel, find the nearest triangle" costs, and that direction is the real
// decision.
static inline uint32_t lbvh_expand_bits(uint32_t v) {
v = (v * 0x00010001u) & 0xFF0000FFu;
v = (v * 0x00000101u) & 0x0F00F00Fu;
v = (v * 0x00000011u) & 0xC30C30C3u;
v = (v * 0x00000005u) & 0x49249249u;
return v;
}
static inline uint32_t lbvh_morton3(float x, float y, float z) {
// 10 bits an axis, 30 total: finer than the leaf size can use.
const float sx = std::min(std::max(x * 1024.0f, 0.0f), 1023.0f);
const float sy = std::min(std::max(y * 1024.0f, 0.0f), 1023.0f);
const float sz = std::min(std::max(z * 1024.0f, 0.0f), 1023.0f);
return (lbvh_expand_bits((uint32_t)sx) << 2) |
(lbvh_expand_bits((uint32_t)sy) << 1) |
lbvh_expand_bits((uint32_t)sz);
}
void DCGeneratorMeshSDF::build_bvh() {
int n = (int)triangles_.size();
tri_indices_.resize(n);
std::iota(tri_indices_.begin(), tri_indices_.end(), 0);
if (n == 0) {
bvh_nodes_.clear();
return;
}
// Centroids once, and the bound they are normalised against.
tri_centroids_.resize((size_t)n * 3);
float cmin[3] = {std::numeric_limits<float>::max(), std::numeric_limits<float>::max(),
std::numeric_limits<float>::max()};
float cmax[3] = {-std::numeric_limits<float>::max(), -std::numeric_limits<float>::max(),
-std::numeric_limits<float>::max()};
for (int i = 0; i < n; i++) {
const Triangle& t = triangles_[i];
const float c[3] = {(t.v0.x + t.v1.x + t.v2.x) / 3.0f,
(t.v0.y + t.v1.y + t.v2.y) / 3.0f,
(t.v0.z + t.v1.z + t.v2.z) / 3.0f};
for (int a = 0; a < 3; a++) {
tri_centroids_[(size_t)i * 3 + a] = c[a];
cmin[a] = std::min(cmin[a], c[a]);
cmax[a] = std::max(cmax[a], c[a]);
}
}
const float inv[3] = {
(cmax[0] - cmin[0]) > 0.0f ? 1.0f / (cmax[0] - cmin[0]) : 0.0f,
(cmax[1] - cmin[1]) > 0.0f ? 1.0f / (cmax[1] - cmin[1]) : 0.0f,
(cmax[2] - cmin[2]) > 0.0f ? 1.0f / (cmax[2] - cmin[2]) : 0.0f};
// Sort (code, index) PAIRS, not indices with an indirect comparator: the
// comparator is then a plain integer compare on contiguous memory.
std::vector<std::pair<uint32_t, int>> order((size_t)n);
tri_morton_.resize((size_t)n);
for (int i = 0; i < n; i++) {
const float* c = &tri_centroids_[(size_t)i * 3];
const uint32_t m = lbvh_morton3((c[0] - cmin[0]) * inv[0],
(c[1] - cmin[1]) * inv[1],
(c[2] - cmin[2]) * inv[2]);
tri_morton_[(size_t)i] = m;
order[(size_t)i] = {m, i};
}
std::sort(order.begin(), order.end());
for (int i = 0; i < n; i++) {
tri_indices_[(size_t)i] = order[(size_t)i].second;
}
// Worst case up front so the parallel build never reallocates under a
// thread; shrunk to what was used afterwards.
bvh_nodes_.clear();
bvh_nodes_.resize((size_t)(2 * n > 0 ? 2 * n : 1));
bvh_alloc_.store(1);
build_lbvh(0, 0, n, 0);
bvh_nodes_.resize((size_t)bvh_alloc_.load());
bvh_nodes_.shrink_to_fit();
}
void DCGeneratorMeshSDF::build_lbvh(int node_idx, int start, int count, int depth) {
BVHNode& node = bvh_nodes_[node_idx];
node.tri_start = start;
const int MAX_LEAF = 4;
if (count <= MAX_LEAF || depth >= 62) {
node.tri_count = count;
node.left = -1;
node.right = -1;
// Contiguous range after the Morton sort: a sequential scan.
float mn[3] = {std::numeric_limits<float>::max(), std::numeric_limits<float>::max(),
std::numeric_limits<float>::max()};
float mx[3] = {-std::numeric_limits<float>::max(), -std::numeric_limits<float>::max(),
-std::numeric_limits<float>::max()};
for (int i = start; i < start + count; i++) {
const Triangle& tri = triangles_[tri_indices_[i]];
const Vec3* verts[3] = {&tri.v0, &tri.v1, &tri.v2};
for (int v = 0; v < 3; v++) {
const float coords[3] = {verts[v]->x, verts[v]->y, verts[v]->z};
for (int a = 0; a < 3; a++) {
mn[a] = std::min(mn[a], coords[a]);
mx[a] = std::max(mx[a], coords[a]);
}
}
}
std::memcpy(node.aabb_min, mn, sizeof(mn));
std::memcpy(node.aabb_max, mx, sizeof(mx));
return;
}
// Karras' split: the highest bit where the range's first and last codes
// differ, found by binary search. Equal codes (duplicate centroids) fall
// back to the middle so the recursion always makes progress.
const uint32_t first = tri_morton_[tri_indices_[start]];
const uint32_t last = tri_morton_[tri_indices_[start + count - 1]];
int split;
if (first == last) {
split = start + count / 2;
} else {
const int common = __builtin_clz(first ^ last);
split = start;
int step = count - 1;
do {
step = (step + 1) >> 1;
const int probe = split + step;
if (probe < start + count - 1) {
const uint32_t code = tri_morton_[tri_indices_[probe]];
if (__builtin_clz(first ^ code) > common) {
split = probe;
}
}
} while (step > 1);
split += 1;
if (split <= start || split >= start + count) {
split = start + count / 2;
}
}
const int left_count = split - start;
const int right_count = count - left_count;
const int left_idx = bvh_alloc_.fetch_add(2);
const int right_idx = left_idx + 1;
node.left = left_idx;
node.right = right_idx;
node.tri_count = 0; // interior
const int PARALLEL_DEPTH = 4;
const int PARALLEL_MIN_TRIS = 8192;
if (depth < PARALLEL_DEPTH && count > PARALLEL_MIN_TRIS) {
std::thread left_thread([this, left_idx, start, left_count, depth]() {
build_lbvh(left_idx, start, left_count, depth + 1);
});
build_lbvh(right_idx, split, right_count, depth + 1);
left_thread.join();
} else {
build_lbvh(left_idx, start, left_count, depth + 1);
build_lbvh(right_idx, split, right_count, depth + 1);
}
// The union of two children we already have -- not another scan of the
// triangles. This is the whole point of the Morton order.
const BVHNode& l = bvh_nodes_[left_idx];
const BVHNode& r = bvh_nodes_[right_idx];
BVHNode& self = bvh_nodes_[node_idx];
for (int a = 0; a < 3; a++) {
self.aabb_min[a] = std::min(l.aabb_min[a], r.aabb_min[a]);
self.aabb_max[a] = std::max(l.aabb_max[a], r.aabb_max[a]);
}
}
void DCGeneratorMeshSDF::build_bvh_recursive(int node_idx, int start, int count, int depth) {
BVHNode& node = bvh_nodes_[node_idx];
node.tri_start = start;
node.tri_count = count;
// Compute AABB
float mn[3] = {std::numeric_limits<float>::max(),
std::numeric_limits<float>::max(),
std::numeric_limits<float>::max()};
float mx[3] = {-std::numeric_limits<float>::max(),
-std::numeric_limits<float>::max(),
-std::numeric_limits<float>::max()};
for (int i = start; i < start + count; i++) {
const Triangle& tri = triangles_[tri_indices_[i]];
const Vec3* verts[3] = {&tri.v0, &tri.v1, &tri.v2};
for (int v = 0; v < 3; v++) {
float coords[3] = {verts[v]->x, verts[v]->y, verts[v]->z};
for (int a = 0; a < 3; a++) {
mn[a] = std::min(mn[a], coords[a]);
mx[a] = std::max(mx[a], coords[a]);
}
}
}
std::memcpy(node.aabb_min, mn, sizeof(mn));
std::memcpy(node.aabb_max, mx, sizeof(mx));
const int MAX_LEAF = 4;
const int MAX_DEPTH = 64;
if (count <= MAX_LEAF || depth >= MAX_DEPTH) {
return;
}
// Split along axis with largest extent at centroid median
float extent[3] = {mx[0] - mn[0], mx[1] - mn[1], mx[2] - mn[2]};
int axis = 0;
if (extent[1] > extent[axis]) axis = 1;
if (extent[2] > extent[axis]) axis = 2;
const float* cent = tri_centroids_.data();
int mid = start + count / 2;
std::nth_element(tri_indices_.begin() + start,
tri_indices_.begin() + mid,
tri_indices_.begin() + start + count,
[cent, axis](int a, int b) {
return cent[(size_t)a * 3 + axis] < cent[(size_t)b * 3 + axis];
});
int left_count = mid - start;
int right_count = count - left_count;
if (left_count == 0 || right_count == 0) {
return; // Can't split, keep as leaf
}
const int left_idx = bvh_alloc_.fetch_add(2);
const int right_idx = left_idx + 1;
bvh_nodes_[node_idx].left = left_idx;
bvh_nodes_[node_idx].right = right_idx;
bvh_nodes_[node_idx].tri_count = 0; // interior node
// ── PARALLEL BELOW THE TOP LEVELS ───────────────────────────────────────
//
// Measured 2026-09-08 on the bridged voronoi mountain: this build is 10.5 s
// of a 12.7 s Mesh to SDF Grid, over 5.16M triangles, on one core. The
// narrow band had already taken the QUERIES to 1.2 s, which is what made
// the build visible as the remaining cost.
//
// Safe because the two subtrees own DISJOINT ranges of tri_indices_ --
// nth_element has already partitioned [start, start+count) about mid, and
// neither side touches the other's half again -- triangles_ is immutable
// after load_arrays, and nodes come from the atomic bump above into an
// array that cannot move.
//
// Only the first few levels spawn: depth 4 is up to 16 subtrees, which
// matches a typical core count, and below that a thread costs more than
// the subtree does. The count guard keeps small meshes entirely serial.
const int PARALLEL_DEPTH = 4;
const int PARALLEL_MIN_TRIS = 8192;
if (depth < PARALLEL_DEPTH && count > PARALLEL_MIN_TRIS) {
std::thread left_thread([this, left_idx, start, left_count, depth]() {
build_bvh_recursive(left_idx, start, left_count, depth + 1);
});
build_bvh_recursive(right_idx, mid, right_count, depth + 1);
left_thread.join();
} else {
build_bvh_recursive(left_idx, start, left_count, depth + 1);
build_bvh_recursive(right_idx, mid, right_count, depth + 1);
}
}
// ---------------------------------------------------------------------------
// SDF evaluation
// ---------------------------------------------------------------------------
// Ray-parity inside/outside test: cast ray along +X, count intersections
// Odd = inside, even = outside. Robust, no pseudonormals needed.
bool DCGeneratorMeshSDF::ray_parity_inside(float x, float y, float z) const {
int crossings = 0;
for (int ti = 0; ti < (int)triangles_.size(); ti++) {
const auto& t = triangles_[ti];
// Ray from (x,y,z) in +X direction. Test intersection with triangle.
// Shift triangle to ray origin
float v0y = t.v0.y - y, v0z = t.v0.z - z;
float v1y = t.v1.y - y, v1z = t.v1.z - z;
float v2y = t.v2.y - y, v2z = t.v2.z - z;
// Check if ray (in YZ plane) passes through triangle projection
// Using sign-of-cross-product edge tests
bool b0 = (v0y * v1z - v0z * v1y) > 0;
bool b1 = (v1y * v2z - v1z * v2y) > 0;
bool b2 = (v2y * v0z - v2z * v0y) > 0;
if (b0 == b1 && b1 == b2) {
// Ray intersects triangle in YZ — check X coordinate
float v0x = t.v0.x - x;
float v1x = t.v1.x - x;
float v2x = t.v2.x - x;
// Barycentric interpolation for intersection X
float denom = (v1y - v2y) * (v0z - v2z) + (v2z - v1z) * (v0y - v2y);
if (std::fabs(denom) < 1e-12f) continue;
float inv = 1.0f / denom;
float u = ((v1y - v2y) * (-v2z) + (v2z - v1z) * (-v2y)) * inv;
float v = ((v2y - v0y) * (-v2z) + (v0z - v2z) * (-v2y)) * inv;
float w = 1.0f - u - v;
float ix = u * v0x + v * v1x + w * v2x;
if (ix > 0) crossings++;
}
}
return (crossings & 1) != 0;
}
// Generalized Winding Number: sum of solid angles subtended by each triangle.
// Result > 0.5 = inside mesh. Robust for non-watertight meshes.
// Reference: Jacobson et al. 2013, "Robust Inside-Outside Segmentation using
// Generalized Winding Numbers"
float DCGeneratorMeshSDF::winding_number(float x, float y, float z) const {
double winding = 0.0;
for (const auto& tri : triangles_) {
// Vectors from query point to triangle vertices
double ax = tri.v0.x - x, ay = tri.v0.y - y, az = tri.v0.z - z;
double bx = tri.v1.x - x, by = tri.v1.y - y, bz = tri.v1.z - z;
double cx = tri.v2.x - x, cy = tri.v2.y - y, cz = tri.v2.z - z;
double la = std::sqrt(ax*ax + ay*ay + az*az);
double lb = std::sqrt(bx*bx + by*by + bz*bz);
double lc = std::sqrt(cx*cx + cy*cy + cz*cz);
if (la < 1e-10 || lb < 1e-10 || lc < 1e-10) {
winding += 0.5; // on vertex = ambiguous, treat as inside
continue;
}
// Solid angle via Van Oosterom and Strackee formula
double numerator = ax*(by*cz - bz*cy) - ay*(bx*cz - bz*cx) + az*(bx*cy - by*cx);
double denominator = la*lb*lc
+ (ax*bx + ay*by + az*bz) * lc
+ (ax*cx + ay*cy + az*cz) * lb
+ (bx*cx + by*cy + bz*cz) * la;
winding += std::atan2(numerator, denominator);
}
return (float)(winding / (2.0 * 3.14159265358979323846));
}
float DCGeneratorMeshSDF::evaluate(float x, float y, float z, float /*lod_voxel_size*/) {
if (!loaded_) return -1.0f;
Vec3 closest_normal;
float signed_dist_sq = closest_distance_sq(x, y, z, closest_normal);
float dist = std::sqrt(std::fabs(signed_dist_sq));
if (unsigned_only_) {
// No pseudonormals were built, so the sign carried in signed_dist_sq is
// meaningless. The caller supplies it.
return dist;
}
float sign = (signed_dist_sq >= 0.0f) ? 1.0f : -1.0f;
// Convention: positive = solid (inside), negative = air (outside)
return dist * sign;
}
void DCGeneratorMeshSDF::gradient(float x, float y, float z, float /*lod_voxel_size*/, float normal[3]) {
if (!loaded_) { normal[0] = 0; normal[1] = 1; normal[2] = 0; return; }
// Use exact pseudonormal from BVH closest-point query.
// At sharp creases, this returns the face normal (not averaged),
// giving the QEF solver the exact constraint for sharp feature preservation.
Vec3 closest_normal;
closest_distance_sq(x, y, z, closest_normal);
normal[0] = closest_normal.x;
normal[1] = closest_normal.y;
normal[2] = closest_normal.z;
}
bool DCGeneratorMeshSDF::get_bounds(float min[3], float max[3]) const {
if (!loaded_) return false;
std::memcpy(min, bounds_min_, sizeof(bounds_min_));
std::memcpy(max, bounds_max_, sizeof(bounds_max_));
return true;
}
// ---------------------------------------------------------------------------
// Closest distance query (returns signed squared distance)
// ---------------------------------------------------------------------------
float DCGeneratorMeshSDF::closest_distance_sq(float x, float y, float z, Vec3& closest_normal) const {
float best_dist_sq = std::numeric_limits<float>::max();
closest_normal = Vec3(0, 1, 0);
if (bvh_nodes_.empty()) return best_dist_sq;
query_bvh(0, x, y, z, best_dist_sq, closest_normal);
return best_dist_sq;
}
// ---------------------------------------------------------------------------
// BVH traversal
// ---------------------------------------------------------------------------
float DCGeneratorMeshSDF::aabb_dist_sq(float px, float py, float pz,
const float aabb_min[3], const float aabb_max[3]) {
float dist_sq = 0.0f;
float coords[3] = {px, py, pz};
for (int i = 0; i < 3; i++) {
if (coords[i] < aabb_min[i]) {
float d = aabb_min[i] - coords[i];
dist_sq += d * d;
} else if (coords[i] > aabb_max[i]) {
float d = coords[i] - aabb_max[i];
dist_sq += d * d;
}
}
return dist_sq;
}
void DCGeneratorMeshSDF::query_bvh(int node_idx, float x, float y, float z,
float& best_dist_sq, Vec3& best_normal) const {
const BVHNode& node = bvh_nodes_[node_idx];
// Prune: AABB distance must be less than current best unsigned distance
float aabb_dsq = aabb_dist_sq(x, y, z, node.aabb_min, node.aabb_max);
if (aabb_dsq > std::fabs(best_dist_sq)) return;
if (node.tri_count > 0) {
// Leaf: test all triangles
for (int i = node.tri_start; i < node.tri_start + node.tri_count; i++) {
Vec3 normal;
float dsq = point_triangle_dist_sq(x, y, z, triangles_[tri_indices_[i]], normal);
if (std::fabs(dsq) < std::fabs(best_dist_sq)) {
best_dist_sq = dsq;
best_normal = normal;
}
}
return;
}
// Interior: traverse closer child first
if (node.left >= 0 && node.right >= 0) {
float d_left = aabb_dist_sq(x, y, z,
bvh_nodes_[node.left].aabb_min,
bvh_nodes_[node.left].aabb_max);
float d_right = aabb_dist_sq(x, y, z,
bvh_nodes_[node.right].aabb_min,
bvh_nodes_[node.right].aabb_max);
int first = node.left, second = node.right;
if (d_right < d_left) std::swap(first, second);
query_bvh(first, x, y, z, best_dist_sq, best_normal);
query_bvh(second, x, y, z, best_dist_sq, best_normal);
} else if (node.left >= 0) {
query_bvh(node.left, x, y, z, best_dist_sq, best_normal);
} else if (node.right >= 0) {
query_bvh(node.right, x, y, z, best_dist_sq, best_normal);
}
}
// ---------------------------------------------------------------------------
// Point-triangle closest point with pseudonormal sign determination
//
// Returns SIGNED squared distance: positive = outside, negative = inside.
// Sign determined by dot(P - closest_point, pseudonormal_at_feature).
// Based on Ericson's "Real-Time Collision Detection" for the closest point,
// with Baerentzen & Aanaes pseudonormal method for sign.
// ---------------------------------------------------------------------------
float DCGeneratorMeshSDF::point_triangle_dist_sq(float px, float py, float pz,
const Triangle& tri, Vec3& normal) {
Vec3 P(px, py, pz);
Vec3 A = tri.v0, B = tri.v1, C = tri.v2;
Vec3 AB = B - A;
Vec3 AC = C - A;
Vec3 AP = P - A;
float d1 = AB.dot(AP);
float d2 = AC.dot(AP);
// Vertex A region
if (d1 <= 0.0f && d2 <= 0.0f) {
normal = tri.vert_normals[0];
Vec3 diff = P - A;
float sign_dot = diff.dot(normal);
return sign_dot >= 0.0f ? diff.length_sq() : -diff.length_sq();
}
Vec3 BP = P - B;
float d3 = AB.dot(BP);
float d4 = AC.dot(BP);
// Vertex B region
if (d3 >= 0.0f && d4 <= d3) {
normal = tri.vert_normals[1];
Vec3 diff = P - B;
float sign_dot = diff.dot(normal);
return sign_dot >= 0.0f ? diff.length_sq() : -diff.length_sq();
}
// Edge AB region
float vc = d1 * d4 - d3 * d2;
if (vc <= 0.0f && d1 >= 0.0f && d3 <= 0.0f) {
float v = d1 / (d1 - d3);
Vec3 closest = A + AB * v;
normal = tri.edge_normals[0];
Vec3 diff = P - closest;
float sign_dot = diff.dot(normal);
return sign_dot >= 0.0f ? diff.length_sq() : -diff.length_sq();
}
Vec3 CP = P - C;
float d5 = AB.dot(CP);
float d6 = AC.dot(CP);
// Vertex C region
if (d6 >= 0.0f && d5 <= d6) {
normal = tri.vert_normals[2];
Vec3 diff = P - C;
float sign_dot = diff.dot(normal);
return sign_dot >= 0.0f ? diff.length_sq() : -diff.length_sq();
}
// Edge AC region
float vb = d5 * d2 - d1 * d6;
if (vb <= 0.0f && d2 >= 0.0f && d6 <= 0.0f) {
float w = d2 / (d2 - d6);
Vec3 closest = A + AC * w;
normal = tri.edge_normals[2]; // edge v2-v0 (AC)
Vec3 diff = P - closest;
float sign_dot = diff.dot(normal);
return sign_dot >= 0.0f ? diff.length_sq() : -diff.length_sq();
}
// Edge BC region
float va = d3 * d6 - d5 * d4;
if (va <= 0.0f && (d4 - d3) >= 0.0f && (d5 - d6) >= 0.0f) {
float w = (d4 - d3) / ((d4 - d3) + (d5 - d6));
Vec3 closest = B + (C - B) * w;
normal = tri.edge_normals[1]; // edge v1-v2 (BC)
Vec3 diff = P - closest;
float sign_dot = diff.dot(normal);
return sign_dot >= 0.0f ? diff.length_sq() : -diff.length_sq();
}
// Interior face region
float denom = 1.0f / (va + vb + vc);
float v = vb * denom;
float w = vc * denom;
Vec3 closest = A + AB * v + AC * w;
normal = tri.normal;
Vec3 diff = P - closest;
float sign_dot = diff.dot(normal);
return sign_dot >= 0.0f ? diff.length_sq() : -diff.length_sq();
}
int DCGeneratorMeshSDF::triangle_count() const {
return static_cast<int>(triangles_.size());
}
// ---------------------------------------------------------------------------
// GPU BVH flattening
// ---------------------------------------------------------------------------
int DCGeneratorMeshSDF::flatten_for_gpu(std::vector<uint8_t>& node_data,
std::vector<uint8_t>& tri_data,
std::vector<uint32_t>& tri_indices_out) const {
if (!loaded_ || bvh_nodes_.empty()) return 0;
// GPU BVH node: 32 bytes
// [aabb_min(3f), left_or_start(u), aabb_max(3f), count(u)]
// Interior: count=0, left_or_start=left child, right=left+1
// Leaf: count>0, left_or_start=tri_start
struct GpuBVHNodePacked {
float aabb_min[3];
uint32_t left_or_start;
float aabb_max[3];
uint32_t count;
};
static_assert(sizeof(GpuBVHNodePacked) == 32, "GpuBVHNode must be 32 bytes");
int num_nodes = static_cast<int>(bvh_nodes_.size());
node_data.resize(num_nodes * 32);
for (int i = 0; i < num_nodes; i++) {
const BVHNode& n = bvh_nodes_[i];
GpuBVHNodePacked packed;
std::memcpy(packed.aabb_min, n.aabb_min, sizeof(float) * 3);
std::memcpy(packed.aabb_max, n.aabb_max, sizeof(float) * 3);
if (n.tri_count > 0) {
packed.left_or_start = static_cast<uint32_t>(n.tri_start);
packed.count = static_cast<uint32_t>(n.tri_count);
} else {
packed.left_or_start = static_cast<uint32_t>(n.left);
packed.count = 0;
}
std::memcpy(node_data.data() + i * 32, &packed, 32);
}
// GPU triangle: 128 bytes (30 floats + 2 padding)
// [v0(3), v1(3), v2(3), face_normal(3),
// edge_normal0(3), edge_normal1(3), edge_normal2(3),
// vert_normal0(3), vert_normal1(3), vert_normal2(3), pad(2)]
int num_tris = static_cast<int>(triangles_.size());
tri_data.resize(num_tris * 128);
std::memset(tri_data.data(), 0, tri_data.size());
for (int i = 0; i < num_tris; i++) {
const Triangle& t = triangles_[i];
float* dst = reinterpret_cast<float*>(tri_data.data() + i * 128);
// v0
dst[0] = t.v0.x; dst[1] = t.v0.y; dst[2] = t.v0.z;
// v1
dst[3] = t.v1.x; dst[4] = t.v1.y; dst[5] = t.v1.z;
// v2
dst[6] = t.v2.x; dst[7] = t.v2.y; dst[8] = t.v2.z;
// face normal
dst[9] = t.normal.x; dst[10] = t.normal.y; dst[11] = t.normal.z;
// edge normals
for (int e = 0; e < 3; e++) {
dst[12 + e*3 + 0] = t.edge_normals[e].x;
dst[12 + e*3 + 1] = t.edge_normals[e].y;
dst[12 + e*3 + 2] = t.edge_normals[e].z;
}
// vertex normals
for (int v = 0; v < 3; v++) {
dst[21 + v*3 + 0] = t.vert_normals[v].x;
dst[21 + v*3 + 1] = t.vert_normals[v].y;
dst[21 + v*3 + 2] = t.vert_normals[v].z;
}
// dst[30], dst[31] = padding (already zeroed)
}
// Triangle indices
tri_indices_out.resize(tri_indices_.size());
for (size_t i = 0; i < tri_indices_.size(); i++) {
tri_indices_out[i] = static_cast<uint32_t>(tri_indices_[i]);
}
return num_nodes;
}