#ifndef DC_GENERATOR_MESH_SDF_H
#define DC_GENERATOR_MESH_SDF_H
#include "dc_generator.h"
#include <atomic>
#include <vector>
#include <string>
#include <cstdint>
#include <cmath>
#include <algorithm>
#include <map>
// DCGenerator that loads a triangle mesh (OBJ) and computes signed distance
// using the pseudonormal method (Baerentzen & Aanaes 2005).
//
// Convention: positive = solid (inside mesh), negative = air (outside mesh).
class DCGeneratorMeshSDF : public DCGenerator {
public:
bool load_obj(const std::string& path);
// Load from raw arrays: positions (x,y,z triples), indices (triangle vertex indices)
bool load_arrays(const float* positions, int num_verts,
const int* indices, int num_tris);
// Pack BVH data into flat GPU-friendly arrays.
// node_data: 32 bytes per node (aabb + left_or_start + count)
// tri_data: 128 bytes per triangle (verts + all pseudonormals + padding)
// Returns number of BVH nodes, or 0 on failure.
int flatten_for_gpu(std::vector<uint8_t>& node_data,
std::vector<uint8_t>& tri_data,
std::vector<uint32_t>& tri_indices_out) const;
float evaluate(float x, float y, float z, float lod_voxel_size) override;
void gradient(float x, float y, float z, float lod_voxel_size, float normal[3]) override;
bool get_bounds(float min[3], float max[3]) const override;
int triangle_count() const;
private:
struct Vec3 {
float x, y, z;
Vec3() : x(0), y(0), z(0) {}
Vec3(float x_, float y_, float z_) : x(x_), y(y_), z(z_) {}
Vec3 operator-(const Vec3& b) const { return {x - b.x, y - b.y, z - b.z}; }
Vec3 operator+(const Vec3& b) const { return {x + b.x, y + b.y, z + b.z}; }
Vec3 operator*(float s) const { return {x * s, y * s, z * s}; }
float dot(const Vec3& b) const { return x * b.x + y * b.y + z * b.z; }
Vec3 cross(const Vec3& b) const {
return {y * b.z - z * b.y, z * b.x - x * b.z, x * b.y - y * b.x};
}
float length_sq() const { return x * x + y * y + z * z; }
float length() const { return std::sqrt(length_sq()); }
Vec3 normalized() const {
float l = length();
return l > 1e-12f ? Vec3(x / l, y / l, z / l) : Vec3(0, 0, 0);
}
};
struct Triangle {
Vec3 v0, v1, v2;
Vec3 normal; // face normal
Vec3 edge_normals[3]; // pseudonormals at edges (e01, e12, e20)
Vec3 vert_normals[3]; // pseudonormals at vertices (v0, v1, v2)
};
struct BVHNode {
float aabb_min[3], aabb_max[3];
int left = -1, right = -1; // child indices (-1 = leaf)
int tri_start = 0, tri_count = 0;
};
std::vector<Triangle> triangles_;
std::vector<int> tri_indices_; // for BVH reordering
std::vector<BVHNode> bvh_nodes_;
// Nodes are handed out by an ATOMIC BUMP, not push_back: the build is
// parallel below the top few levels and push_back mutates the vector's
// size. bvh_nodes_ is resized to its worst case (2n-1 for n triangles)
// before the build so no allocation can move it under a running thread,
// then shrunk to what was used.
std::atomic<int> bvh_alloc_{0};
// Triangle centroids, precomputed once, 3 floats per triangle.
//
// nth_element's comparator used to recompute a centroid from the Triangle
// on EVERY comparison -- n log n random dereferences into a 5.16M-element
// array of 9-float structs, which is why threading the build only bought
// 1.42x: it is memory bound, not compute bound. Comparing floats in a flat
// array instead keeps the sort in cache.
std::vector<float> tri_centroids_;
// Morton code per triangle, and the LBVH build that consumes it.
std::vector<uint32_t> tri_morton_;
// Skip the Baerentzen-Aanaes pseudonormals and answer UNSIGNED distance.
//
// They exist only to give a query point a SIGN, and they cost 5554 ms of a
// 6663 ms load_arrays on a 5.16M triangle mesh -- five times the BVH build
// -- because they need per-vertex AND per-EDGE angle-weighted normals,
// ~15.5M edge entries here. A caller that already knows inside from outside
// does not need them. Mesh to SDF Grid does: it floods the outside from the
// grid boundary.
bool unsigned_only_ = false;
void build_lbvh(int node_idx, int start, int count, int depth);
public:
void set_unsigned_only(bool v) { unsigned_only_ = v; }
private:
float bounds_min_[3] = {};
float bounds_max_[3] = {};
bool loaded_ = false;
bool _build_from_verts_indices(const std::vector<Vec3>& vertices,
const std::vector<int>& face_indices);
void compute_pseudonormals(const std::vector<Vec3>& vertices,
const std::vector<int>& face_indices);
void build_bvh();
void build_bvh_recursive(int node_idx, int start, int count, int depth);
bool ray_parity_inside(float x, float y, float z) const;
float winding_number(float x, float y, float z) const;
float closest_distance_sq(float x, float y, float z, Vec3& closest_normal) const;
void query_bvh(int node_idx, float x, float y, float z,
float& best_dist_sq, Vec3& best_normal) const;
// Squared distance from point to AABB
static float aabb_dist_sq(float px, float py, float pz,
const float aabb_min[3], const float aabb_max[3]);
// Point-triangle closest point with feature identification for pseudonormal
static float point_triangle_dist_sq(float px, float py, float pz,
const Triangle& tri, Vec3& normal);
};
#endif // DC_GENERATOR_MESH_SDF_H