#include "dc_contour.h"
#include "dc_edit.h"
#include <cmath>
#include <cstring>
#include <unordered_map>
#include <unordered_set>
#include <functional>
#include <algorithm>
// Define DC_DEBUG to enable diagnostic prints. Disabled by default.
// #define DC_DEBUG
#include <cstdio>
#ifdef DC_DEBUG
#define DC_PRINT(...) printf(__VA_ARGS__)
#else
#define DC_PRINT(...) ((void)0)
#endif
// ---------------------------------------------------------------------------
// Edge table: 12 edges of a unit cube (Gildea corner ordering).
// Bit layout: x = bit2, y = bit1, z = bit0
// ---------------------------------------------------------------------------
static const int edgeCorners[12][2] = {
// X-axis edges (differ in bit 2)
{0, 4}, {1, 5}, {2, 6}, {3, 7},
// Y-axis edges (differ in bit 1)
{0, 2}, {1, 3}, {4, 6}, {5, 7},
// Z-axis edges (differ in bit 0)
{0, 1}, {2, 3}, {4, 5}, {6, 7},
};
// ---------------------------------------------------------------------------
// Helpers
// ---------------------------------------------------------------------------
static inline float lerp(float a, float b, float t) {
return a + t * (b - a);
}
static inline float length3(const float v[3]) {
return std::sqrt(v[0]*v[0] + v[1]*v[1] + v[2]*v[2]);
}
static inline void normalize3(float v[3]) {
float len = length3(v);
if (len > 1e-10f) {
v[0] /= len; v[1] /= len; v[2] /= len;
}
}
// ---------------------------------------------------------------------------
// Hermite data collection
// ---------------------------------------------------------------------------
static float find_zero_crossing(
float p0[3], float p1[3],
float d0, float d1,
DensityFunc density_fn, void* userdata)
{
// Always bisect the full [0,1] range using endpoint signs.
// Caller guarantees d0 and d1 have opposite signs.
// 8 bisection steps give 1/256 precision (Gildea recommends high precision here).
float lo = 0.0f, hi = 1.0f;
float d_lo = d0;
for (int i = 0; i < 8; i++) {
float mid = 0.5f * (lo + hi);
float pm[3] = {
lerp(p0[0], p1[0], mid),
lerp(p0[1], p1[1], mid),
lerp(p0[2], p1[2], mid),
};
float dm = density_fn(pm[0], pm[1], pm[2], userdata);
if ((d_lo > 0.0f) == (dm > 0.0f)) {
lo = mid; d_lo = dm;
} else {
hi = mid;
}
}
return 0.5f * (lo + hi);
}
static void compute_normal(
float wx, float wy, float wz, float h,
DensityFunc density_fn, void* userdata,
float normal[3])
{
// Negate gradient: density is positive inside (solid), so gradient points
// inward. Negating gives outward-facing surface normals.
normal[0] = density_fn(wx - h, wy, wz, userdata) - density_fn(wx + h, wy, wz, userdata);
normal[1] = density_fn(wx, wy - h, wz, userdata) - density_fn(wx, wy + h, wz, userdata);
normal[2] = density_fn(wx, wy, wz - h, userdata) - density_fn(wx, wy, wz + h, userdata);
normalize3(normal);
}
static void collect_hermite_data(
OctreeNode* node,
float voxel_scale, float origin[3],
DensityFunc density_fn, void* userdata)
{
if (!node) return;
if (node->type == OctreeNode::INTERNAL) {
for (int i = 0; i < 8; i++)
collect_hermite_data(node->children[i], voxel_scale, origin,
density_fn, userdata);
return;
}
if (node->type != OctreeNode::LEAF) return;
node->qef.clear();
node->average_normal[0] = node->average_normal[1] = node->average_normal[2] = 0.0f;
float h = voxel_scale * 0.01f; // gradient step — 1% of cell size
int edge_count = 0;
for (int e = 0; e < 12; e++) {
int c0 = edgeCorners[e][0];
int c1 = edgeCorners[e][1];
bool s0 = (node->corners & (1 << c0)) != 0;
bool s1 = (node->corners & (1 << c1)) != 0;
if (s0 == s1) continue;
float wp0[3], wp1[3];
for (int a = 0; a < 3; a++) {
wp0[a] = origin[a] + (node->min[a] + CHILD_MIN_OFFSETS[c0][a] * node->size) * voxel_scale;
wp1[a] = origin[a] + (node->min[a] + CHILD_MIN_OFFSETS[c1][a] * node->size) * voxel_scale;
}
// Use cached corner densities from octree construction
float d0 = node->corner_densities[c0];
float d1 = node->corner_densities[c1];
float t = find_zero_crossing(wp0, wp1, d0, d1, density_fn, userdata);
float crossing[3] = {
lerp(wp0[0], wp1[0], t),
lerp(wp0[1], wp1[1], t),
lerp(wp0[2], wp1[2], t),
};
float normal[3];
compute_normal(crossing[0], crossing[1], crossing[2], h,
density_fn, userdata, normal);
node->qef.add(crossing, normal);
node->average_normal[0] += normal[0];
node->average_normal[1] += normal[1];
node->average_normal[2] += normal[2];
edge_count++;
}
if (edge_count > 0) {
normalize3(node->average_normal);
// Compute cell center for QEF center bias.
float cell_center[3], cell_min[3], cell_max[3];
for (int a = 0; a < 3; a++) {
cell_min[a] = origin[a] + node->min[a] * voxel_scale;
cell_max[a] = origin[a] + (node->min[a] + node->size) * voxel_scale;
cell_center[a] = 0.5f * (cell_min[a] + cell_max[a]);
}
// Biased QEF solve (Tikhonov regularization)
node->qef.solve_biased(node->position, cell_center, 0.02f);
// If QEF solution is outside cell, fall back to masspoint.
// Hard clamping places vertices on cell faces far from the isosurface,
// creating thin "fin" triangles. The masspoint (average of edge crossings)
// is always on/near the isosurface and inside the cell.
bool outside = false;
for (int a = 0; a < 3; a++) {
if (node->position[a] < cell_min[a] || node->position[a] > cell_max[a]) {
outside = true;
break;
}
}
if (outside) {
for (int a = 0; a < 3; a++)
node->position[a] = node->qef.masspoint[a];
// Clamp masspoint too (shouldn't be needed, but safety)
for (int a = 0; a < 3; a++) {
if (node->position[a] < cell_min[a]) node->position[a] = cell_min[a];
if (node->position[a] > cell_max[a]) node->position[a] = cell_max[a];
}
}
} else {
for (int a = 0; a < 3; a++)
node->position[a] = origin[a] + (node->min[a] + 0.5f * node->size) * voxel_scale;
}
}
// ---------------------------------------------------------------------------
// Grid-accelerated hermite data collection
// Uses pre-computed DensityGrid for bisection + normals (trilinear interp).
// ---------------------------------------------------------------------------
static float find_zero_crossing_grid(
float /*p0*/[3], float /*p1*/[3],
float d0, float d1,
const DensityGrid& /*grid*/)
{
// Pure linear interpolation — no bisection needed on a smooth
// trilinear grid. The grid is already a piecewise-linear field,
// so t = d0/(d0-d1) gives the exact zero crossing within the cell.
float t = d0 / (d0 - d1);
if (t < 0.0f) return 0.0f;
if (t > 1.0f) return 1.0f;
return t;
}
// Snap edge crossing to nearest box face if within tolerance.
// Returns true if snapped (out_normal and out_pos overridden with axis-aligned data).
// Snap edge crossing to the box face perpendicular to the edge axis.
// edge_axis: 0=X, 1=Y, 2=Z — only checks faces on this axis to prevent
// wrong snapping at box corners where multiple faces are equidistant.
static bool snap_to_box_face(const float crossing[3],
const std::vector<TerrainEdit>* edits,
float voxel_scale,
int edge_axis,
float out_normal[3], float out_pos[3]) {
if (!edits) return false;
float tol = voxel_scale * 0.6f;
for (const auto& e : *edits) {
if (e.shape != TerrainEdit::BOX) continue;
float center[3] = {e.x, e.y, e.z};
float local[3] = {
crossing[0] - center[0],
crossing[1] - center[1],
crossing[2] - center[2]
};
// Must be near the box
bool near = true;
for (int a = 0; a < 3; a++) {
if (std::fabs(local[a]) > e.half_ext[a] + tol) { near = false; break; }
}
if (!near) continue;
// Only check the face perpendicular to the edge axis
float d_pos = std::fabs(local[edge_axis] - e.half_ext[edge_axis]);
float d_neg = std::fabs(local[edge_axis] + e.half_ext[edge_axis]);
float best_dist = d_pos;
float best_sign = 1.0f;
if (d_neg < best_dist) { best_dist = d_neg; best_sign = -1.0f; }
if (best_dist < tol) {
out_pos[0] = crossing[0];
out_pos[1] = crossing[1];
out_pos[2] = crossing[2];
out_pos[edge_axis] = center[edge_axis] + best_sign * e.half_ext[edge_axis];
out_normal[0] = 0.0f; out_normal[1] = 0.0f; out_normal[2] = 0.0f;
out_normal[edge_axis] = best_sign;
return true;
}
}
return false;
}
// Mark corners that are inside box edits. This determines which edges
// get axis-aligned normals (box) vs grid-sampled normals (terrain).
static void mark_box_corners(OctreeNode* node, float voxel_scale, float origin[3],
const std::vector<TerrainEdit>* edits) {
if (!node || !edits) return;
if (node->type == OctreeNode::INTERNAL) {
for (int i = 0; i < 8; i++)
mark_box_corners(node->children[i], voxel_scale, origin, edits);
return;
}
if (node->type != OctreeNode::LEAF && node->type != OctreeNode::PSEUDO) return;
node->corner_materials = 0;
for (int c = 0; c < 8; c++) {
float wx = origin[0] + (node->min[0] + CHILD_MIN_OFFSETS[c][0] * node->size) * voxel_scale;
float wy = origin[1] + (node->min[1] + CHILD_MIN_OFFSETS[c][1] * node->size) * voxel_scale;
float wz = origin[2] + (node->min[2] + CHILD_MIN_OFFSETS[c][2] * node->size) * voxel_scale;
for (const auto& e : *edits) {
if (e.shape == TerrainEdit::MESH && e.mesh_sdf) {
// Mesh edit: corner is inside if mesh SDF evaluates positive
float mdx = wx - e.x, mdy = wy - e.y, mdz = wz - e.z;
if (mdx*mdx + mdy*mdy + mdz*mdz <= e.radius * e.radius) {
float lx = mdx / e.scale, ly = mdy / e.scale, lz = mdz / e.scale;
if (e.mesh_sdf->evaluate(lx, ly, lz, 0.0f) > 0.0f) {
node->corner_materials |= (1 << c);
break;
}
}
continue;
}
if (e.shape != TerrainEdit::BOX) continue;
if (std::fabs(wx - e.x) <= e.half_ext[0] &&
std::fabs(wy - e.y) <= e.half_ext[1] &&
std::fabs(wz - e.z) <= e.half_ext[2]) {
node->corner_materials |= (1 << c);
break;
}
}
}
}
static void collect_hermite_data_grid(
OctreeNode* node,
float voxel_scale, float origin[3],
const DensityGrid& grid,
const std::vector<TerrainEdit>* edits = nullptr)
{
if (!node) return;
if (node->type == OctreeNode::INTERNAL) {
for (int i = 0; i < 8; i++)
collect_hermite_data_grid(node->children[i], voxel_scale, origin, grid, edits);
return;
}
if (node->type != OctreeNode::LEAF) return;
node->qef.clear();
node->average_normal[0] = node->average_normal[1] = node->average_normal[2] = 0.0f;
int edge_count = 0;
for (int e = 0; e < 12; e++) {
int c0 = edgeCorners[e][0];
int c1 = edgeCorners[e][1];
bool s0 = (node->corners & (1 << c0)) != 0;
bool s1 = (node->corners & (1 << c1)) != 0;
if (s0 == s1) continue;
float wp0[3], wp1[3];
for (int a = 0; a < 3; a++) {
wp0[a] = origin[a] + (node->min[a] + CHILD_MIN_OFFSETS[c0][a] * node->size) * voxel_scale;
wp1[a] = origin[a] + (node->min[a] + CHILD_MIN_OFFSETS[c1][a] * node->size) * voxel_scale;
}
float d0 = node->corner_densities[c0];
float d1 = node->corner_densities[c1];
float t = find_zero_crossing_grid(wp0, wp1, d0, d1, grid);
float crossing[3] = {
lerp(wp0[0], wp1[0], t),
lerp(wp0[1], wp1[1], t),
lerp(wp0[2], wp1[2], t),
};
int edge_axis = (e < 4) ? 0 : (e < 8) ? 1 : 2;
float normal[3];
float snap_pos[3];
// Material-aware normals: if either corner is inside a box edit,
// try to snap to the box face. Uses material check (corner_materials)
// with OR to catch all box-adjacent edges, then proximity snap
// confirms and provides the exact face position + axis-aligned normal.
bool c0_is_box = (node->corner_materials & (1 << c0)) != 0;
bool c1_is_box = (node->corner_materials & (1 << c1)) != 0;
bool edge_near_box = c0_is_box || c1_is_box;
// Tolerance: material-confirmed edges use wider snap range
float snap_tol = edge_near_box ? voxel_scale * 2.0f : voxel_scale * 0.6f;
if (snap_to_box_face(crossing, edits, snap_tol, edge_axis, normal, snap_pos)) {
node->qef.add(snap_pos, normal);
} else if (edge_near_box && edits) {
// Mesh SDF edge: use exact crossing position + gradient.
// Grid linear interpolation staircases at voxel boundaries;
// bisecting the mesh SDF directly gives sub-voxel precision.
bool used_mesh = false;
for (const auto& edit : *edits) {
if (edit.shape != TerrainEdit::MESH || !edit.mesh_sdf) continue;
float mdx = crossing[0] - edit.x, mdy = crossing[1] - edit.y, mdz = crossing[2] - edit.z;
if (mdx*mdx + mdy*mdy + mdz*mdz > edit.radius * edit.radius) continue;
// Bisect along the edge using mesh SDF for exact crossing
float bp0[3] = {wp0[0], wp0[1], wp0[2]};
float bp1[3] = {wp1[0], wp1[1], wp1[2]};
for (int step = 0; step < 8; step++) {
float mid[3] = {(bp0[0]+bp1[0])*0.5f, (bp0[1]+bp1[1])*0.5f, (bp0[2]+bp1[2])*0.5f};
float lx = (mid[0]-edit.x)/edit.scale, ly = (mid[1]-edit.y)/edit.scale, lz = (mid[2]-edit.z)/edit.scale;
float d = edit.mesh_sdf->evaluate(lx, ly, lz, 0.0f);
// Combined with terrain: apply same CSG op
float td = grid.sample(mid[0], mid[1], mid[2]);
float combined = (edit.strength < 0) ? std::min(td, -d) : std::max(td, d);
if ((combined >= 0) == s0) { bp0[0]=mid[0]; bp0[1]=mid[1]; bp0[2]=mid[2]; }
else { bp1[0]=mid[0]; bp1[1]=mid[1]; bp1[2]=mid[2]; }
}
float exact[3] = {(bp0[0]+bp1[0])*0.5f, (bp0[1]+bp1[1])*0.5f, (bp0[2]+bp1[2])*0.5f};
// Exact gradient at the precise crossing point
float lx = (exact[0]-edit.x)/edit.scale, ly = (exact[1]-edit.y)/edit.scale, lz = (exact[2]-edit.z)/edit.scale;
edit.mesh_sdf->gradient(lx, ly, lz, 0.0f, normal);
node->qef.add(exact, normal);
used_mesh = true;
break;
}
if (!used_mesh) {
grid.sample_normal(crossing[0], crossing[1], crossing[2], normal);
node->qef.add(crossing, normal);
}
} else {
grid.sample_normal(crossing[0], crossing[1], crossing[2], normal);
node->qef.add(crossing, normal);
}
node->average_normal[0] += normal[0];
node->average_normal[1] += normal[1];
node->average_normal[2] += normal[2];
edge_count++;
}
if (edge_count > 0) {
normalize3(node->average_normal);
// Compute cell center for QEF center bias.
float cell_center[3];
for (int a = 0; a < 3; a++) {
cell_center[a] = origin[a] + (node->min[a] + 0.5f * node->size) * voxel_scale;
}
// Solve QEF with Tikhonov regularization (Boris the Brave center bias).
// The bias_weight pulls the solution toward cell center when normals are
// nearly colinear (shallow slopes, flat surfaces). This prevents the
// solution from flying outside the cell without needing hard clamping,
// which causes stair-stepping artifacts.
// Weight 0.1 is gentle enough to preserve sharp features while stabilizing
// degenerate cases. For comparison: edge_count normals contribute weight ~1
// each to ATA diagonal, so 0.1 is ~10% of a single normal's influence.
node->qef.solve_biased(node->position, cell_center, 0.02f);
// If QEF solution is outside cell, fall back to masspoint.
// Hard clamping places vertices on cell faces far from the isosurface,
// creating thin "fin" triangles. The masspoint (average of edge crossings)
// is always on/near the isosurface and inside the cell.
float cell_min[3], cell_max[3];
for (int a = 0; a < 3; a++) {
cell_min[a] = origin[a] + node->min[a] * voxel_scale;
cell_max[a] = origin[a] + (node->min[a] + node->size) * voxel_scale;
}
bool outside = false;
for (int a = 0; a < 3; a++) {
if (node->position[a] < cell_min[a] || node->position[a] > cell_max[a]) {
outside = true;
break;
}
}
if (outside) {
for (int a = 0; a < 3; a++)
node->position[a] = node->qef.masspoint[a];
for (int a = 0; a < 3; a++) {
if (node->position[a] < cell_min[a]) node->position[a] = cell_min[a];
if (node->position[a] > cell_max[a]) node->position[a] = cell_max[a];
}
}
} else {
for (int a = 0; a < 3; a++)
node->position[a] = origin[a] + (node->min[a] + 0.5f * node->size) * voxel_scale;
}
}
// ---------------------------------------------------------------------------
// Contouring: edge-per-leaf approach with spatial hash and edge dedup set.
//
// For each leaf, iterate its 12 edges. For edges with a sign change,
// look up the 4 cells sharing that edge. Emit a quad/triangles for the
// cells that have vertices. Use an edge ID set to avoid duplicate quads.
// ---------------------------------------------------------------------------
// ---------------------------------------------------------------------------
// Flat-array leaf map: O(1) lookup by voxel position, no hashing overhead.
// Grid is offset so min corner is at (0,0,0). Size known at gather time.
// ---------------------------------------------------------------------------
struct FlatLeafMap {
std::vector<OctreeNode*> data;
int grid_size = 0; // cells per axis
void allocate(int size) {
grid_size = size;
data.assign(size * size * size, nullptr);
}
inline void set(int x, int y, int z, OctreeNode* node) {
if (x >= 0 && x < grid_size && y >= 0 && y < grid_size && z >= 0 && z < grid_size)
data[z * grid_size * grid_size + y * grid_size + x] = node;
}
inline OctreeNode* get(int x, int y, int z) const {
if (x < 0 || x >= grid_size || y < 0 || y >= grid_size || z < 0 || z >= grid_size)
return nullptr;
return data[z * grid_size * grid_size + y * grid_size + x];
}
inline int size() const {
int count = 0;
for (auto* p : data) if (p) count++;
return count;
}
};
static void gather_leaves_flat(OctreeNode* node, FlatLeafMap& map) {
if (!node) return;
if (node->type == OctreeNode::LEAF || node->type == OctreeNode::PSEUDO) {
map.set(node->min[0], node->min[1], node->min[2], node);
return;
}
for (int i = 0; i < 8; i++)
gather_leaves_flat(node->children[i], map);
}
// Edge ID for dedup: pack (axis, ex, ey, ez) into uint64.
static inline uint64_t edge_id(int axis, int ex, int ey, int ez) {
uint64_t a = static_cast<uint64_t>(axis);
uint64_t x = static_cast<uint64_t>(ex + (1 << 19)) & 0xFFFFF;
uint64_t y = static_cast<uint64_t>(ey + (1 << 19)) & 0xFFFFF;
uint64_t z = static_cast<uint64_t>(ez + (1 << 19)) & 0xFFFFF;
return (a << 60) | (x << 40) | (y << 20) | z;
}
// Legacy hash-based leaf map (kept for dc_generate_mesh compatibility)
typedef std::unordered_map<uint64_t, OctreeNode*> LeafMap;
static inline uint64_t leaf_key(int x, int y, int z) {
return pack_voxel_key(x, y, z);
}
static void gather_leaves(OctreeNode* node, LeafMap& map) {
if (!node) return;
if (node->type == OctreeNode::LEAF || node->type == OctreeNode::PSEUDO) {
map[leaf_key(node->min[0], node->min[1], node->min[2])] = node;
return;
}
for (int i = 0; i < 8; i++)
gather_leaves(node->children[i], map);
}
static int get_vertex_index(OctreeNode* node, DCMeshData& mesh) {
if (node->vertex_index >= 0) return node->vertex_index;
int idx = (int)(mesh.positions.size() / 3);
mesh.positions.push_back(node->position[0]);
mesh.positions.push_back(node->position[1]);
mesh.positions.push_back(node->position[2]);
mesh.normals.push_back(node->average_normal[0]);
mesh.normals.push_back(node->average_normal[1]);
mesh.normals.push_back(node->average_normal[2]);
node->vertex_index = idx;
return idx;
}
// Offsets from the edge position to the 4 cells sharing that edge.
// The edge at (ex,ey,ez) along axis A is shared by 4 cells differing
// in the two transverse directions by 0 or -1 (cell size units).
//
static const int edgeCellOffsets[3][4][3] = {
// X-axis edge: transverse Y, Z (right-hand: Y×Z = +X)
{{ 0, 0, 0}, { 0,-1, 0}, { 0, 0,-1}, { 0,-1,-1}},
// Y-axis edge: transverse Z, X (right-hand: Z×X = +Y)
// q1=Z, q2=X — swapped from naive X,Z to fix winding
{{ 0, 0, 0}, { 0, 0,-1}, {-1, 0, 0}, {-1, 0,-1}},
// Z-axis edge: transverse X, Y (right-hand: X×Y = +Z)
{{ 0, 0, 0}, {-1, 0, 0}, { 0,-1, 0}, {-1,-1, 0}},
};
// Shared quad emission logic for contouring. Emits triangles for a sign-change
// edge shared by up to 4 cells. Called by both hash-map and flat-array paths.
static inline void emit_edge_quads(
int vert_indices[4], int nv, bool s0, DCMeshData& mesh)
{
if (nv < 3) return;
// Winding: when s0 is true (c0 solid, c1 air), the outward face
// should point from solid toward air. Empirically: !s0 gets the
// (0,1,3)+(0,3,2) winding (which produces +axis geometric normal).
if (nv == 4) {
// Shortest-diagonal triangulation: pick the shorter diagonal to
// reduce visual creases and avoid bowtie triangles on non-planar quads.
const float* p0 = &mesh.positions[vert_indices[0] * 3];
const float* p1 = &mesh.positions[vert_indices[1] * 3];
const float* p2 = &mesh.positions[vert_indices[2] * 3];
const float* p3 = &mesh.positions[vert_indices[3] * 3];
float d03_x = p0[0] - p3[0], d03_y = p0[1] - p3[1], d03_z = p0[2] - p3[2];
float d12_x = p1[0] - p2[0], d12_y = p1[1] - p2[1], d12_z = p1[2] - p2[2];
float dist03_sq = d03_x * d03_x + d03_y * d03_y + d03_z * d03_z;
float dist12_sq = d12_x * d12_x + d12_y * d12_y + d12_z * d12_z;
if (dist03_sq <= dist12_sq) {
// Split along 0-3 diagonal: triangles (0,1,3) and (0,3,2)
if (!s0) {
mesh.indices.push_back(vert_indices[0]);
mesh.indices.push_back(vert_indices[1]);
mesh.indices.push_back(vert_indices[3]);
mesh.indices.push_back(vert_indices[0]);
mesh.indices.push_back(vert_indices[3]);
mesh.indices.push_back(vert_indices[2]);
} else {
mesh.indices.push_back(vert_indices[0]);
mesh.indices.push_back(vert_indices[3]);
mesh.indices.push_back(vert_indices[1]);
mesh.indices.push_back(vert_indices[0]);
mesh.indices.push_back(vert_indices[2]);
mesh.indices.push_back(vert_indices[3]);
}
} else {
// Split along 1-2 diagonal: triangles (0,1,2) and (1,3,2)
if (!s0) {
mesh.indices.push_back(vert_indices[0]);
mesh.indices.push_back(vert_indices[1]);
mesh.indices.push_back(vert_indices[2]);
mesh.indices.push_back(vert_indices[1]);
mesh.indices.push_back(vert_indices[3]);
mesh.indices.push_back(vert_indices[2]);
} else {
mesh.indices.push_back(vert_indices[0]);
mesh.indices.push_back(vert_indices[2]);
mesh.indices.push_back(vert_indices[1]);
mesh.indices.push_back(vert_indices[1]);
mesh.indices.push_back(vert_indices[2]);
mesh.indices.push_back(vert_indices[3]);
}
}
} else {
// 3 of 4 cells present: emit 1 triangle with correct winding.
int tri[3]; int ti = 0;
for (int q = 0; q < 4 && ti < 3; q++) {
if (vert_indices[q] >= 0) tri[ti++] = vert_indices[q];
}
if (!s0) {
mesh.indices.push_back(tri[0]);
mesh.indices.push_back(tri[1]);
mesh.indices.push_back(tri[2]);
} else {
mesh.indices.push_back(tri[0]);
mesh.indices.push_back(tri[2]);
mesh.indices.push_back(tri[1]);
}
}
}
static void contour_leaves(LeafMap& leaves, DCMeshData& mesh) {
std::unordered_set<uint64_t> processed_edges;
for (auto& [key, node] : leaves) {
int sz = node->size;
for (int e = 0; e < 12; e++) {
int c0 = edgeCorners[e][0];
int c1 = edgeCorners[e][1];
bool s0 = (node->corners & (1 << c0)) != 0;
bool s1 = (node->corners & (1 << c1)) != 0;
if (s0 == s1) continue;
int axis = (e < 4) ? 0 : (e < 8) ? 1 : 2;
int ex = node->min[0] + CHILD_MIN_OFFSETS[c0][0] * sz;
int ey = node->min[1] + CHILD_MIN_OFFSETS[c0][1] * sz;
int ez = node->min[2] + CHILD_MIN_OFFSETS[c0][2] * sz;
uint64_t eid = edge_id(axis, ex, ey, ez);
if (!processed_edges.insert(eid).second) continue;
int vert_indices[4] = {-1, -1, -1, -1};
int nv = 0;
for (int q = 0; q < 4; q++) {
int cx = ex + edgeCellOffsets[axis][q][0] * sz;
int cy = ey + edgeCellOffsets[axis][q][1] * sz;
int cz = ez + edgeCellOffsets[axis][q][2] * sz;
auto it = leaves.find(leaf_key(cx, cy, cz));
if (it != leaves.end()) {
vert_indices[q] = get_vertex_index(it->second, mesh);
nv++;
}
}
emit_edge_quads(vert_indices, nv, s0, mesh);
}
}
DC_PRINT("[DCContour] leaves=%d verts=%d tris=%d\n",
(int)leaves.size(),
(int)(mesh.positions.size()/3), (int)(mesh.indices.size()/3));
}
// ---------------------------------------------------------------------------
// Reset vertex indices for re-generation
// ---------------------------------------------------------------------------
static void reset_vertex_indices(OctreeNode* node) {
if (!node) return;
node->vertex_index = -1;
for (int i = 0; i < 8; i++)
reset_vertex_indices(node->children[i]);
}
// ---------------------------------------------------------------------------
// Public API
// ---------------------------------------------------------------------------
void dc_generate_mesh(
OctreeNode* root,
float voxel_scale, float origin[3],
DensityFunc density_fn, void* userdata,
DCMeshData& out_mesh)
{
if (!root) return;
out_mesh.positions.clear();
out_mesh.normals.clear();
out_mesh.indices.clear();
reset_vertex_indices(root);
collect_hermite_data(root, voxel_scale, origin, density_fn, userdata);
LeafMap leaves;
gather_leaves(root, leaves);
contour_leaves(leaves, out_mesh);
}
// ---------------------------------------------------------------------------
// Flat-array contouring: same logic as contour_leaves but uses FlatLeafMap
// for O(1) neighbor lookups instead of hash map lookups.
// ---------------------------------------------------------------------------
static void contour_leaves_flat(FlatLeafMap& leaves, int grid_size, int overlap, DCMeshData& mesh) {
// Flat bit array for edge dedup: 3 axes * (grid_size+1)^3 edges.
// For grid_size=64 this is ~96 KB — much cheaper than a hash set.
int edge_stride = grid_size + 1;
int edge_total = 3 * edge_stride * edge_stride * edge_stride;
std::vector<uint8_t> edge_bits((edge_total + 7) / 8, 0);
auto edge_bit_idx = [edge_stride](int axis, int bx, int by, int bz) -> int {
return axis * edge_stride * edge_stride * edge_stride +
bz * edge_stride * edge_stride + by * edge_stride + bx;
};
for (int z = 0; z < grid_size; z++)
for (int y = 0; y < grid_size; y++)
for (int x = 0; x < grid_size; x++) {
OctreeNode* node = leaves.get(x, y, z);
if (!node) continue;
int sz = node->size;
for (int e = 0; e < 12; e++) {
int c0 = edgeCorners[e][0];
int c1 = edgeCorners[e][1];
bool s0 = (node->corners & (1 << c0)) != 0;
bool s1 = (node->corners & (1 << c1)) != 0;
if (s0 == s1) continue;
int axis = (e < 4) ? 0 : (e < 8) ? 1 : 2;
int ex = node->min[0] + CHILD_MIN_OFFSETS[c0][0] * sz;
int ey = node->min[1] + CHILD_MIN_OFFSETS[c0][1] * sz;
int ez = node->min[2] + CHILD_MIN_OFFSETS[c0][2] * sz;
int bidx = edge_bit_idx(axis, ex, ey, ez);
int byte_idx = bidx >> 3;
int bit_mask = 1 << (bidx & 7);
if (edge_bits[byte_idx] & bit_mask) continue;
edge_bits[byte_idx] |= bit_mask;
// Cell-ownership: nominal bounds with +1 upper extension.
// The extension keeps boundary faces that connect the last nominal
// cell to the first overlap cell. Adjacent chunks generate matching
// faces from their own overlap data. Without density smoothing,
// the vertices match exactly — watertight with minor z-fighting.
if (overlap > 0) {
int nom_min = overlap;
int nom_max = grid_size - overlap + 1;
if (ex < nom_min || ex >= nom_max ||
ey < nom_min || ey >= nom_max ||
ez < nom_min || ez >= nom_max)
continue;
}
int vert_indices[4] = {-1, -1, -1, -1};
int nv = 0;
for (int q = 0; q < 4; q++) {
int cx = ex + edgeCellOffsets[axis][q][0] * sz;
int cy = ey + edgeCellOffsets[axis][q][1] * sz;
int cz = ez + edgeCellOffsets[axis][q][2] * sz;
OctreeNode* neighbor = leaves.get(cx, cy, cz);
if (neighbor) {
vert_indices[q] = get_vertex_index(neighbor, mesh);
nv++;
}
}
emit_edge_quads(vert_indices, nv, s0, mesh);
}
}
DC_PRINT("[DCContour] leaves=%d verts=%d tris=%d\n",
leaves.size(),
(int)(mesh.positions.size()/3), (int)(mesh.indices.size()/3));
}
void dc_generate_mesh_grid(
OctreeNode* root,
float voxel_scale, float origin[3],
const DensityGrid& grid,
int overlap,
DCMeshData& out_mesh,
const std::vector<TerrainEdit>* edits)
{
if (!root) return;
out_mesh.positions.clear();
out_mesh.normals.clear();
out_mesh.indices.clear();
out_mesh.collision_indices.clear();
reset_vertex_indices(root);
if (edits) mark_box_corners(root, voxel_scale, origin, edits);
collect_hermite_data_grid(root, voxel_scale, origin, grid, edits);
FlatLeafMap leaves;
leaves.allocate(grid.grid_size);
gather_leaves_flat(root, leaves);
contour_leaves_flat(leaves, grid.grid_size, overlap, out_mesh);
}
// ---------------------------------------------------------------------------
// Auto-smooth: split vertices at sharp edges to reveal DC sharp features.
//
// Algorithm:
// 1. Compute face normal for each triangle (cross product of edges)
// 2. Build per-vertex face adjacency list
// 3. For each vertex, cluster its adjacent faces by normal similarity
// (flood fill: faces in same cluster have angle < threshold)
// 4. One cluster = smooth vertex (average normal)
// Multiple clusters = split vertex, each gets its cluster's avg normal
// ---------------------------------------------------------------------------
static void cross3(const float a[3], const float b[3], float out[3]) {
out[0] = a[1]*b[2] - a[2]*b[1];
out[1] = a[2]*b[0] - a[0]*b[2];
out[2] = a[0]*b[1] - a[1]*b[0];
}
static float dot3_as(const float a[3], const float b[3]) {
return a[0]*b[0] + a[1]*b[1] + a[2]*b[2];
}
void dc_auto_smooth(DCMeshData& mesh, float angle_threshold_degrees) {
int tri_count = (int)(mesh.indices.size() / 3);
int vert_count = (int)(mesh.positions.size() / 3);
if (tri_count == 0 || vert_count == 0) return;
float cos_threshold = std::cos(angle_threshold_degrees * 3.14159265f / 180.0f);
// Step 1: Compute face normals (area-weighted via unnormalized cross product)
std::vector<float> face_normals(tri_count * 3);
for (int f = 0; f < tri_count; f++) {
int i0 = mesh.indices[f*3+0], i1 = mesh.indices[f*3+1], i2 = mesh.indices[f*3+2];
const float* p0 = &mesh.positions[i0*3];
const float* p1 = &mesh.positions[i1*3];
const float* p2 = &mesh.positions[i2*3];
float e1[3] = { p1[0]-p0[0], p1[1]-p0[1], p1[2]-p0[2] };
float e2[3] = { p2[0]-p0[0], p2[1]-p0[1], p2[2]-p0[2] };
float* fn = &face_normals[f*3];
cross3(e2, e1, fn); // reversed for CW winding
float len = std::sqrt(fn[0]*fn[0] + fn[1]*fn[1] + fn[2]*fn[2]);
if (len < 1e-10f) { fn[0]=0; fn[1]=1; fn[2]=0; }
}
// Step 2: Build per-vertex face adjacency
std::vector<std::vector<int>> vert_faces(vert_count);
for (int f = 0; f < tri_count; f++) {
for (int j = 0; j < 3; j++) {
vert_faces[mesh.indices[f*3+j]].push_back(f);
}
}
// Step 3: Cluster faces per vertex using topology-aware flood fill
std::vector<float> new_positions;
std::vector<float> new_normals;
std::vector<float> new_materials; // route ② texturing: material must ride the split
std::vector<float> new_proj; // Stage-2 per-face projection candidate 0 rides the split
std::vector<float> new_proj2; // ...and candidate 1 (crease other face)
std::vector<float> new_curv; // curvature rides the split too
bool has_mat = ((int)mesh.materials.size() == vert_count);
bool has_proj = ((int)mesh.proj.size() == vert_count);
bool has_proj2 = ((int)mesh.proj2.size() == vert_count);
bool has_curv = ((int)mesh.curv.size() == vert_count);
std::vector<int> corner_remap(mesh.indices.size(), -1);
new_positions.reserve(mesh.positions.size());
new_normals.reserve(mesh.normals.size());
if (has_mat) new_materials.reserve(mesh.materials.size());
if (has_proj) new_proj.reserve(mesh.proj.size());
if (has_proj2) new_proj2.reserve(mesh.proj2.size());
if (has_curv) new_curv.reserve(mesh.curv.size());
std::vector<bool> visited;
for (int v = 0; v < vert_count; v++) {
const auto& faces = vert_faces[v];
int nf = (int)faces.size();
if (nf == 0) continue;
// Build local adjacency: two faces in faces[] are adjacent at v
// if they share another vertex (i.e., share an edge incident on v).
std::vector<int> other_a(nf), other_b(nf);
for (int i = 0; i < nf; i++) {
int fi = faces[i];
int va = mesh.indices[fi*3+0], vb = mesh.indices[fi*3+1], vc = mesh.indices[fi*3+2];
if (va == v) { other_a[i] = vb; other_b[i] = vc; }
else if (vb == v) { other_a[i] = va; other_b[i] = vc; }
else { other_a[i] = va; other_b[i] = vb; }
}
std::vector<std::vector<int>> adj(nf);
for (int i = 0; i < nf; i++) {
for (int j = i + 1; j < nf; j++) {
if (other_a[i] == other_a[j] || other_a[i] == other_b[j] ||
other_b[i] == other_a[j] || other_b[i] == other_b[j]) {
adj[i].push_back(j);
adj[j].push_back(i);
}
}
}
visited.assign(nf, false);
int cluster_count = 0;
struct Cluster {
int new_idx;
float normal[3];
};
std::vector<Cluster> clusters;
for (int start = 0; start < nf; start++) {
if (visited[start]) continue;
Cluster cl;
cl.new_idx = (int)(new_positions.size() / 3);
cl.normal[0] = cl.normal[1] = cl.normal[2] = 0.0f;
new_positions.push_back(mesh.positions[v*3+0]);
new_positions.push_back(mesh.positions[v*3+1]);
new_positions.push_back(mesh.positions[v*3+2]);
new_normals.push_back(0); new_normals.push_back(0); new_normals.push_back(0);
if (has_mat) new_materials.push_back(mesh.materials[v]); // split inherits material
if (has_proj) new_proj.push_back(mesh.proj[v]); // ...and both projection candidates
if (has_proj2) new_proj2.push_back(mesh.proj2[v]);
if (has_curv) new_curv.push_back(mesh.curv[v]);
std::vector<int> stack;
stack.push_back(start);
visited[start] = true;
while (!stack.empty()) {
int si = stack.back(); stack.pop_back();
int fi = faces[si];
const float* fn = &face_normals[fi*3];
// Weighted normal: face_area × corner_angle (Blender's "Face Area and Angle")
// fn is already area-weighted (unnormalized cross product, |cross| = 2*area).
// Compute corner angle at vertex v in this face for additional weighting.
int i0 = mesh.indices[fi*3+0], i1 = mesh.indices[fi*3+1], i2 = mesh.indices[fi*3+2];
const float* pv = &mesh.positions[v*3];
const float* pa;
const float* pb;
if (i0 == v) { pa = &mesh.positions[i1*3]; pb = &mesh.positions[i2*3]; }
else if (i1 == v) { pa = &mesh.positions[i0*3]; pb = &mesh.positions[i2*3]; }
else { pa = &mesh.positions[i0*3]; pb = &mesh.positions[i1*3]; }
float ea[3] = { pa[0]-pv[0], pa[1]-pv[1], pa[2]-pv[2] };
float eb[3] = { pb[0]-pv[0], pb[1]-pv[1], pb[2]-pv[2] };
float la = std::sqrt(ea[0]*ea[0]+ea[1]*ea[1]+ea[2]*ea[2]);
float lb = std::sqrt(eb[0]*eb[0]+eb[1]*eb[1]+eb[2]*eb[2]);
float corner_angle = 1.0f;
if (la > 1e-10f && lb > 1e-10f) {
float d = (ea[0]*eb[0]+ea[1]*eb[1]+ea[2]*eb[2]) / (la * lb);
d = std::max(-1.0f, std::min(1.0f, d));
corner_angle = std::acos(d);
}
cl.normal[0] += fn[0] * corner_angle;
cl.normal[1] += fn[1] * corner_angle;
cl.normal[2] += fn[2] * corner_angle;
for (int j = 0; j < 3; j++) {
if (mesh.indices[fi*3+j] == v) corner_remap[fi*3+j] = cl.new_idx;
}
float fn_norm[3] = {fn[0], fn[1], fn[2]};
float fn_len = std::sqrt(fn_norm[0]*fn_norm[0]+fn_norm[1]*fn_norm[1]+fn_norm[2]*fn_norm[2]);
if (fn_len > 1e-10f) { fn_norm[0]/=fn_len; fn_norm[1]/=fn_len; fn_norm[2]/=fn_len; }
// Only flood to topologically adjacent faces (shared edge at v)
for (int ni : adj[si]) {
if (visited[ni]) continue;
const float* fn2 = &face_normals[faces[ni]*3];
float fn2_norm[3] = {fn2[0], fn2[1], fn2[2]};
float fn2_len = std::sqrt(fn2_norm[0]*fn2_norm[0]+fn2_norm[1]*fn2_norm[1]+fn2_norm[2]*fn2_norm[2]);
if (fn2_len > 1e-10f) { fn2_norm[0]/=fn2_len; fn2_norm[1]/=fn2_len; fn2_norm[2]/=fn2_len; }
if (dot3_as(fn_norm, fn2_norm) >= cos_threshold) {
visited[ni] = true;
stack.push_back(ni);
}
}
}
clusters.push_back(std::move(cl));
cluster_count++;
}
// Always use weighted face normals (area × corner_angle) for ALL vertices.
// Previously: single-cluster used gradient normals, multi-cluster used face normals.
// This caused visible bands at transitions where adjacent vertices disagreed on
// which path to take. Weighted face normals are consistent everywhere — they
// converge to the smooth normal on flat surfaces and give correct per-face
// normals at sharp edges. Same approach as Blender's Weighted Normal modifier.
for (auto& cl : clusters) {
float len = std::sqrt(cl.normal[0]*cl.normal[0]+cl.normal[1]*cl.normal[1]+cl.normal[2]*cl.normal[2]);
if (len > 1e-10f) { cl.normal[0]/=len; cl.normal[1]/=len; cl.normal[2]/=len; }
new_normals[cl.new_idx*3+0] = cl.normal[0];
new_normals[cl.new_idx*3+1] = cl.normal[1];
new_normals[cl.new_idx*3+2] = cl.normal[2];
}
}
// Step 4: Rebuild mesh
std::vector<int> new_indices(mesh.indices.size());
for (size_t i = 0; i < mesh.indices.size(); i++) {
new_indices[i] = corner_remap[i];
}
int old_verts = vert_count;
int new_verts = (int)(new_positions.size() / 3);
mesh.positions = std::move(new_positions);
mesh.normals = std::move(new_normals);
mesh.indices = std::move(new_indices);
if (has_mat) mesh.materials = std::move(new_materials);
if (has_proj) mesh.proj = std::move(new_proj);
if (has_proj2) mesh.proj2 = std::move(new_proj2);
if (has_curv) mesh.curv = std::move(new_curv);
DC_PRINT("[DCAutoSmooth] %d -> %d verts (%d split, threshold=%.0f deg)\n",
old_verts, new_verts, new_verts - old_verts, angle_threshold_degrees);
}
// ---------------------------------------------------------------------------
// Laplacian smoothing: reduce noise/spikes by averaging vertex positions.
//
// For each vertex, compute the average position of all edge-connected
// neighbors, then move the vertex toward that average by `lambda`.
// Sharp edges (where adjacent face normals differ by > threshold) act
// as barriers — vertices on opposite sides of a sharp edge don't
// influence each other, preserving DC's feature detection.
//
// This is Cepero's recommended post-process for DC meshes.
// ---------------------------------------------------------------------------
void dc_laplacian_smooth(DCMeshData& mesh, int iterations, float lambda,
float angle_threshold_degrees) {
int vert_count = (int)(mesh.positions.size() / 3);
int tri_count = (int)(mesh.indices.size() / 3);
if (vert_count == 0 || tri_count == 0) return;
float cos_threshold = std::cos(angle_threshold_degrees * 3.14159265f / 180.0f);
// Build edge adjacency: for each vertex, which other vertices are connected
// by edges that are NOT sharp (i.e., both adjacent faces have similar normals).
// An edge is sharp if the two faces sharing it have normals with dot < threshold.
// Step 1: Build edge-to-face map
struct Edge {
int v0, v1;
bool operator==(const Edge& o) const { return v0 == o.v0 && v1 == o.v1; }
};
struct EdgeHash {
size_t operator()(const Edge& e) const {
return std::hash<int>()(e.v0) ^ (std::hash<int>()(e.v1) << 16);
}
};
// Compute face normals
std::vector<float> face_normals(tri_count * 3);
for (int f = 0; f < tri_count; f++) {
int i0 = mesh.indices[f*3], i1 = mesh.indices[f*3+1], i2 = mesh.indices[f*3+2];
const float* p0 = &mesh.positions[i0*3];
const float* p1 = &mesh.positions[i1*3];
const float* p2 = &mesh.positions[i2*3];
float e1[3] = { p1[0]-p0[0], p1[1]-p0[1], p1[2]-p0[2] };
float e2[3] = { p2[0]-p0[0], p2[1]-p0[1], p2[2]-p0[2] };
float* fn = &face_normals[f*3];
cross3(e2, e1, fn);
float len = std::sqrt(fn[0]*fn[0] + fn[1]*fn[1] + fn[2]*fn[2]);
if (len > 1e-10f) { fn[0]/=len; fn[1]/=len; fn[2]/=len; }
}
// Build edge → face list
std::unordered_map<uint64_t, std::vector<int>> edge_faces;
edge_faces.reserve(tri_count * 3);
for (int f = 0; f < tri_count; f++) {
for (int j = 0; j < 3; j++) {
int a = mesh.indices[f*3+j];
int b = mesh.indices[f*3+((j+1)%3)];
int lo = std::min(a, b), hi = std::max(a, b);
uint64_t key = ((uint64_t)lo << 32) | (uint64_t)hi;
edge_faces[key].push_back(f);
}
}
// Build smooth neighbor lists (only across non-sharp edges)
std::vector<std::vector<int>> neighbors(vert_count);
for (auto& [key, faces] : edge_faces) {
int v0 = (int)(key >> 32);
int v1 = (int)(key & 0xFFFFFFFF);
bool is_sharp = false;
if (faces.size() == 2) {
float d = dot3_as(&face_normals[faces[0]*3], &face_normals[faces[1]*3]);
is_sharp = (d < cos_threshold);
} else if (faces.size() != 2) {
// Boundary or non-manifold edge — treat as sharp
is_sharp = true;
}
if (!is_sharp) {
neighbors[v0].push_back(v1);
neighbors[v1].push_back(v0);
}
}
// Iterative smoothing
std::vector<float> new_pos(mesh.positions.size());
for (int iter = 0; iter < iterations; iter++) {
// Copy current positions
std::copy(mesh.positions.begin(), mesh.positions.end(), new_pos.begin());
for (int v = 0; v < vert_count; v++) {
if (neighbors[v].empty()) continue;
// Average neighbor positions
float avg[3] = {0, 0, 0};
for (int n : neighbors[v]) {
avg[0] += mesh.positions[n*3];
avg[1] += mesh.positions[n*3+1];
avg[2] += mesh.positions[n*3+2];
}
float inv = 1.0f / (float)neighbors[v].size();
avg[0] *= inv; avg[1] *= inv; avg[2] *= inv;
// Move toward average
new_pos[v*3] = mesh.positions[v*3] + lambda * (avg[0] - mesh.positions[v*3]);
new_pos[v*3+1] = mesh.positions[v*3+1] + lambda * (avg[1] - mesh.positions[v*3+1]);
new_pos[v*3+2] = mesh.positions[v*3+2] + lambda * (avg[2] - mesh.positions[v*3+2]);
}
mesh.positions = new_pos;
}
// Recompute normals from smoothed geometry
// (face normals changed, so vertex normals need updating)
// Simple approach: area-weighted face normals averaged per vertex
std::fill(mesh.normals.begin(), mesh.normals.end(), 0.0f);
for (int f = 0; f < tri_count; f++) {
int i0 = mesh.indices[f*3], i1 = mesh.indices[f*3+1], i2 = mesh.indices[f*3+2];
const float* p0 = &mesh.positions[i0*3];
const float* p1 = &mesh.positions[i1*3];
const float* p2 = &mesh.positions[i2*3];
float e1[3] = { p1[0]-p0[0], p1[1]-p0[1], p1[2]-p0[2] };
float e2[3] = { p2[0]-p0[0], p2[1]-p0[1], p2[2]-p0[2] };
float fn[3];
cross3(e2, e1, fn); // area-weighted (unnormalized)
for (int j = 0; j < 3; j++) {
int vi = mesh.indices[f*3+j];
mesh.normals[vi*3] += fn[0];
mesh.normals[vi*3+1] += fn[1];
mesh.normals[vi*3+2] += fn[2];
}
}
for (int v = 0; v < vert_count; v++) {
float* n = &mesh.normals[v*3];
float len = std::sqrt(n[0]*n[0] + n[1]*n[1] + n[2]*n[2]);
if (len > 1e-10f) { n[0]/=len; n[1]/=len; n[2]/=len; }
}
}
// ---------------------------------------------------------------------------
// ---------------------------------------------------------------------------
// dc_convex_decompose: flood-fill adjacent triangles by normal similarity,
// output one point cloud per group for convex hull collision shapes.
// Runs on worker thread — pure C++, no Godot types.
// ---------------------------------------------------------------------------
void dc_convex_decompose(const DCMeshData& mesh, ConvexHullSet& out,
float cos_threshold) {
out.hulls.clear();
int tri_count = static_cast<int>(mesh.indices.size()) / 3;
if (tri_count == 0) return;
const float* pos = mesh.positions.data();
// Compute face normals
std::vector<float> fnormals(tri_count * 3);
for (int i = 0; i < tri_count; i++) {
int i0 = mesh.indices[i*3+0], i1 = mesh.indices[i*3+1], i2 = mesh.indices[i*3+2];
float ax = pos[i1*3+0]-pos[i0*3+0], ay = pos[i1*3+1]-pos[i0*3+1], az = pos[i1*3+2]-pos[i0*3+2];
float bx = pos[i2*3+0]-pos[i0*3+0], by = pos[i2*3+1]-pos[i0*3+1], bz = pos[i2*3+2]-pos[i0*3+2];
float nx = ay*bz - az*by, ny = az*bx - ax*bz, nz = ax*by - ay*bx;
float len = std::sqrt(nx*nx + ny*ny + nz*nz);
if (len > 1e-10f) { nx /= len; ny /= len; nz /= len; }
fnormals[i*3+0] = nx; fnormals[i*3+1] = ny; fnormals[i*3+2] = nz;
}
// Build edge adjacency: edge key -> list of triangle indices
std::unordered_map<uint64_t, std::vector<int>> edge_tris;
edge_tris.reserve(tri_count * 3);
for (int i = 0; i < tri_count; i++) {
for (int e = 0; e < 3; e++) {
int a = mesh.indices[i*3 + e];
int b = mesh.indices[i*3 + ((e+1)%3)];
uint64_t key = (static_cast<uint64_t>(std::min(a,b)) << 32) |
static_cast<uint64_t>(std::max(a,b));
edge_tris[key].push_back(i);
}
}
// Flood-fill groups by normal similarity
std::vector<int> group(tri_count, -1);
int num_groups = 0;
for (int i = 0; i < tri_count; i++) {
if (group[i] >= 0) continue;
int gid = num_groups++;
group[i] = gid;
std::vector<int> stack = { i };
while (!stack.empty()) {
int cur = stack.back(); stack.pop_back();
float cnx = fnormals[cur*3+0], cny = fnormals[cur*3+1], cnz = fnormals[cur*3+2];
for (int e = 0; e < 3; e++) {
int a = mesh.indices[cur*3 + e];
int b = mesh.indices[cur*3 + ((e+1)%3)];
uint64_t key = (static_cast<uint64_t>(std::min(a,b)) << 32) |
static_cast<uint64_t>(std::max(a,b));
auto it = edge_tris.find(key);
if (it == edge_tris.end()) continue;
for (int adj : it->second) {
if (adj == cur || group[adj] >= 0) continue;
float dot = cnx*fnormals[adj*3+0] + cny*fnormals[adj*3+1] + cnz*fnormals[adj*3+2];
if (dot > cos_threshold) {
group[adj] = gid;
stack.push_back(adj);
}
}
}
}
}
// Compute average normal per group
std::vector<float> group_normals(num_groups * 3, 0.0f);
for (int i = 0; i < tri_count; i++) {
int g = group[i];
group_normals[g*3+0] += fnormals[i*3+0];
group_normals[g*3+1] += fnormals[i*3+1];
group_normals[g*3+2] += fnormals[i*3+2];
}
for (int g = 0; g < num_groups; g++) {
float len = std::sqrt(group_normals[g*3+0]*group_normals[g*3+0] +
group_normals[g*3+1]*group_normals[g*3+1] +
group_normals[g*3+2]*group_normals[g*3+2]);
if (len > 1e-10f) {
group_normals[g*3+0] /= len;
group_normals[g*3+1] /= len;
group_normals[g*3+2] /= len;
}
}
// For each group: compute AABB and check if coplanar (all vertices within
// a thin band along the normal). Coplanar groups get a box shape.
// Non-coplanar groups get a convex hull.
const float COPLANAR_THRESHOLD = 0.01f; // max spread along normal to be "flat"
const float BOX_DEPTH = 1.0f; // depth of box into solid for coplanar groups
out.hulls.resize(num_groups);
for (int g = 0; g < num_groups; g++) {
out.hulls[g].type = ConvexHull::HULL; // default
}
// Collect vertices per group and compute AABB + normal spread
struct GroupInfo {
float aabb_min[3] = { 1e30f, 1e30f, 1e30f };
float aabb_max[3] = { -1e30f, -1e30f, -1e30f };
float dot_min = 1e30f, dot_max = -1e30f;
};
std::vector<GroupInfo> ginfo(num_groups);
for (int i = 0; i < tri_count; i++) {
int g = group[i];
float gnx = group_normals[g*3+0];
float gny = group_normals[g*3+1];
float gnz = group_normals[g*3+2];
for (int j = 0; j < 3; j++) {
int vi = mesh.indices[i*3+j];
float px = pos[vi*3+0], py = pos[vi*3+1], pz = pos[vi*3+2];
ginfo[g].aabb_min[0] = std::min(ginfo[g].aabb_min[0], px);
ginfo[g].aabb_min[1] = std::min(ginfo[g].aabb_min[1], py);
ginfo[g].aabb_min[2] = std::min(ginfo[g].aabb_min[2], pz);
ginfo[g].aabb_max[0] = std::max(ginfo[g].aabb_max[0], px);
ginfo[g].aabb_max[1] = std::max(ginfo[g].aabb_max[1], py);
ginfo[g].aabb_max[2] = std::max(ginfo[g].aabb_max[2], pz);
float d = px*gnx + py*gny + pz*gnz;
ginfo[g].dot_min = std::min(ginfo[g].dot_min, d);
ginfo[g].dot_max = std::max(ginfo[g].dot_max, d);
}
}
// Classify groups and build shapes
for (int g = 0; g < num_groups; g++) {
float spread = ginfo[g].dot_max - ginfo[g].dot_min;
if (spread < COPLANAR_THRESHOLD) {
// Coplanar group: emit a box (AABB extended into solid along normal)
out.hulls[g].type = ConvexHull::BOX;
float gnx = group_normals[g*3+0];
float gny = group_normals[g*3+1];
float gnz = group_normals[g*3+2];
// Box surface side = AABB of vertices
for (int a = 0; a < 3; a++) {
out.hulls[g].box_min[a] = ginfo[g].aabb_min[a];
out.hulls[g].box_max[a] = ginfo[g].aabb_max[a];
}
// Extend into solid (ALONG normal — DC face normals point inward)
out.hulls[g].box_min[0] += std::min(0.0f, gnx * BOX_DEPTH);
out.hulls[g].box_min[1] += std::min(0.0f, gny * BOX_DEPTH);
out.hulls[g].box_min[2] += std::min(0.0f, gnz * BOX_DEPTH);
out.hulls[g].box_max[0] += std::max(0.0f, gnx * BOX_DEPTH);
out.hulls[g].box_max[1] += std::max(0.0f, gny * BOX_DEPTH);
out.hulls[g].box_max[2] += std::max(0.0f, gnz * BOX_DEPTH);
}
}
// Collect points for non-box (HULL) groups
for (int i = 0; i < tri_count; i++) {
int g = group[i];
if (out.hulls[g].type != ConvexHull::HULL) continue;
auto& pts = out.hulls[g].points;
for (int j = 0; j < 3; j++) {
int vi = mesh.indices[i*3+j];
pts.push_back(pos[vi*3+0]);
pts.push_back(pos[vi*3+1]);
pts.push_back(pos[vi*3+2]);
}
}
}
// dc_clip_to_aabb: discard faces where ALL 3 vertices are outside the AABB.
// Keeps faces that straddle the boundary (at least one vertex inside).
// Re-indexes vertices to compact the output.
// ---------------------------------------------------------------------------
void dc_clip_to_aabb(DCMeshData& mesh,
const float aabb_min[3], const float aabb_max[3]) {
int tri_count = static_cast<int>(mesh.indices.size()) / 3;
if (tri_count == 0) return;
// Small epsilon tolerance so vertices at exact boundary aren't rejected
// by floating-point rounding from QEF solver.
constexpr float eps = 1e-4f;
auto inside = [&](int vi) -> bool {
float x = mesh.positions[vi*3+0];
float y = mesh.positions[vi*3+1];
float z = mesh.positions[vi*3+2];
return x >= aabb_min[0] - eps && x <= aabb_max[0] + eps &&
y >= aabb_min[1] - eps && y <= aabb_max[1] + eps &&
z >= aabb_min[2] - eps && z <= aabb_max[2] + eps;
};
// Mark which triangles to keep (at least one vertex inside AABB)
DCMeshData out;
std::unordered_map<int, int> vert_remap;
auto get_or_add = [&](int old_vi) -> int {
auto it = vert_remap.find(old_vi);
if (it != vert_remap.end()) return it->second;
int new_vi = static_cast<int>(out.positions.size()) / 3;
vert_remap[old_vi] = new_vi;
out.positions.push_back(mesh.positions[old_vi*3+0]);
out.positions.push_back(mesh.positions[old_vi*3+1]);
out.positions.push_back(mesh.positions[old_vi*3+2]);
out.normals.push_back(mesh.normals[old_vi*3+0]);
out.normals.push_back(mesh.normals[old_vi*3+1]);
out.normals.push_back(mesh.normals[old_vi*3+2]);
return new_vi;
};
for (int t = 0; t < tri_count; t++) {
int i0 = mesh.indices[t*3+0];
int i1 = mesh.indices[t*3+1];
int i2 = mesh.indices[t*3+2];
// Keep face if at least one vertex is inside the nominal AABB
if (inside(i0) || inside(i1) || inside(i2)) {
out.indices.push_back(get_or_add(i0));
out.indices.push_back(get_or_add(i1));
out.indices.push_back(get_or_add(i2));
}
}
mesh = std::move(out);
}
// ---------------------------------------------------------------------------
// Degenerate triangle removal: remove triangles with extreme aspect ratio.
// These are the thin "fin" artifacts from poorly-constrained QEF cells.
// ---------------------------------------------------------------------------
void dc_remove_degenerate_tris(DCMeshData& mesh, float max_aspect_ratio) {
int tri_count = (int)mesh.indices.size() / 3;
if (tri_count == 0) return;
// Rebuild with only good triangles
DCMeshData out;
std::unordered_map<int, int> vert_remap;
auto get_or_add = [&](int old_vi) -> int {
auto it = vert_remap.find(old_vi);
if (it != vert_remap.end()) return it->second;
int new_vi = (int)out.positions.size() / 3;
vert_remap[old_vi] = new_vi;
out.positions.push_back(mesh.positions[old_vi*3+0]);
out.positions.push_back(mesh.positions[old_vi*3+1]);
out.positions.push_back(mesh.positions[old_vi*3+2]);
if (!mesh.normals.empty()) {
out.normals.push_back(mesh.normals[old_vi*3+0]);
out.normals.push_back(mesh.normals[old_vi*3+1]);
out.normals.push_back(mesh.normals[old_vi*3+2]);
}
return new_vi;
};
float max_ar_sq = max_aspect_ratio * max_aspect_ratio;
for (int t = 0; t < tri_count; t++) {
int i0 = mesh.indices[t*3+0];
int i1 = mesh.indices[t*3+1];
int i2 = mesh.indices[t*3+2];
const float* p0 = &mesh.positions[i0*3];
const float* p1 = &mesh.positions[i1*3];
const float* p2 = &mesh.positions[i2*3];
// Edge lengths squared
float e0x = p1[0]-p0[0], e0y = p1[1]-p0[1], e0z = p1[2]-p0[2];
float e1x = p2[0]-p1[0], e1y = p2[1]-p1[1], e1z = p2[2]-p1[2];
float e2x = p0[0]-p2[0], e2y = p0[1]-p2[1], e2z = p0[2]-p2[2];
float l0sq = e0x*e0x + e0y*e0y + e0z*e0z;
float l1sq = e1x*e1x + e1y*e1y + e1z*e1z;
float l2sq = e2x*e2x + e2y*e2y + e2z*e2z;
float longest_sq = std::max({l0sq, l1sq, l2sq});
// Cross product = 2 * area
float cx = e0y*(-e2z) - e0z*(-e2y);
float cy = e0z*(-e2x) - e0x*(-e2z);
float cz = e0x*(-e2y) - e0y*(-e2x);
float area_sq_4 = cx*cx + cy*cy + cz*cz;
// Aspect ratio approximation: longest_edge^2 / (4 * area^2 / longest_edge^2)
// = longest_edge^4 / (4 * area^2)
// We compare: longest^4 / area_sq_4 > max_ar_sq
if (area_sq_4 < 1e-20f) continue; // degenerate (zero area)
float ar_approx = (longest_sq * longest_sq) / area_sq_4;
if (ar_approx > max_ar_sq) continue; // skip degenerate
out.indices.push_back(get_or_add(i0));
out.indices.push_back(get_or_add(i1));
out.indices.push_back(get_or_add(i2));
}
mesh = std::move(out);
}
// ---------------------------------------------------------------------------
// Island removal: find connected components, remove small ones.
// Uses union-find on triangle adjacency (shared vertex = connected).
// ---------------------------------------------------------------------------
void dc_remove_islands(DCMeshData& mesh, int min_triangles) {
int tri_count = (int)mesh.indices.size() / 3;
if (tri_count < 2) return;
int vert_count = (int)mesh.positions.size() / 3;
// Union-Find
std::vector<int> parent(vert_count);
std::vector<int> rank(vert_count, 0);
for (int i = 0; i < vert_count; i++) parent[i] = i;
std::function<int(int)> uf_find = [&](int x) -> int {
while (parent[x] != x) { parent[x] = parent[parent[x]]; x = parent[x]; }
return x;
};
auto unite = [&](int a, int b) {
a = uf_find(a); b = uf_find(b);
if (a == b) return;
if (rank[a] < rank[b]) std::swap(a, b);
parent[b] = a;
if (rank[a] == rank[b]) rank[a]++;
};
// Union vertices that share a triangle
for (int t = 0; t < tri_count; t++) {
int i0 = mesh.indices[t*3+0];
int i1 = mesh.indices[t*3+1];
int i2 = mesh.indices[t*3+2];
unite(i0, i1);
unite(i1, i2);
}
// Count triangles per component
std::unordered_map<int, int> component_size;
for (int t = 0; t < tri_count; t++) {
int root = uf_find(mesh.indices[t*3]);
component_size[root]++;
}
// Rebuild mesh keeping only large components
DCMeshData out;
std::unordered_map<int, int> vert_remap;
auto get_or_add = [&](int old_vi) -> int {
auto it = vert_remap.find(old_vi);
if (it != vert_remap.end()) return it->second;
int new_vi = (int)out.positions.size() / 3;
vert_remap[old_vi] = new_vi;
out.positions.push_back(mesh.positions[old_vi*3+0]);
out.positions.push_back(mesh.positions[old_vi*3+1]);
out.positions.push_back(mesh.positions[old_vi*3+2]);
if (!mesh.normals.empty()) {
out.normals.push_back(mesh.normals[old_vi*3+0]);
out.normals.push_back(mesh.normals[old_vi*3+1]);
out.normals.push_back(mesh.normals[old_vi*3+2]);
}
return new_vi;
};
int removed = 0;
for (int t = 0; t < tri_count; t++) {
int root = uf_find(mesh.indices[t*3]);
if (component_size[root] >= min_triangles) {
out.indices.push_back(get_or_add(mesh.indices[t*3+0]));
out.indices.push_back(get_or_add(mesh.indices[t*3+1]));
out.indices.push_back(get_or_add(mesh.indices[t*3+2]));
} else {
removed++;
}
}
mesh = std::move(out);
}