// ============================================================================
// OUTER SHELL ORACLE — level-mesh-implementation.md §5a.2/§5a.3, hyper-games #25.
//
// The volume to subtract from the terrain field is everything inside the level's
// SEALED BOUNDARY, and that boundary is also where §5a.3 wants the terrain's
// isosurface to terminate. One surface, both halves. This computes it and, more
// importantly, checks whether it actually closes.
//
// WHY THIS EXISTS SEPARATELY FROM test_interior_decomp.cpp. That harness did the
// same classification on RAW per-brush windings and measured a closure ratio of
// 0.6-0.7 where a closed surface must give 1.0 — roughly a third of the enclosing
// surface missing. The defect was point-vs-area: it tests ONE sample at the face
// centroid, so a half-exposed face (the common case wherever two brushes abut)
// has its centroid buried in the neighbour and is rejected whole, leaving a hole
// the size of the exposed part.
//
// The fix is the csg4 skin. The neighbour clip SPLITS a half-buried face into an
// exposed fragment and a buried one, so a per-fragment centroid test becomes
// sound — the fragment is exposed or it is not, with no half-measures. That is
// the step the previous harness tried to skip.
//
// METHOD
// 1. csg4 clip_all -> the exact visible skin (fragments, not whole faces)
// 2. rasterize brushes conservatively + flood fill from the border -> EXTERIOR
// 3. a fragment is SHELL if the cell just outside its centroid is EXTERIOR
// 4. closure: divergence theorem V = 1/3 sum (c.n)A over the shell, compared
// against the flood fill's own solid+interior volume
//
// (4) IS THE POINT. Face counts and areas cannot tell a closed surface from one
// with holes in it; two independent volume measurements can. A shell that does
// not close is a subtract that leaks, and the leak would show up in game as
// terrain inside rooms — intermittently, since what leaks depends on chunk and
// LOD.
//
// CONSERVATIVE rasterization is required, not a refinement: Quake walls are
// commonly 16 units = 0.448 m at world_scale 0.028, so sampling cell CENTRES at
// 1 m misses them, the flood pours through, and every room is marked exterior.
//
// SKY BRUSHES are required for the same reason from the other side: a sky brush
// is SOLID and only its rendering is special, so dropping the 19 of them leaves
// holes in the boundary and the fill leaks out through the ceiling.
//
// Input: /tmp/startmap_brushes_sky.bin
// (godot4 --headless --script res://scenes/test/oracle_dump.gd -- --sky)
// Build: g++ -O2 -std=c++17 -I libdc/include -o /tmp/outer_shell \
// tests/test_outer_shell.cpp libdc/src/dc_csg4.cpp
// ============================================================================
#include "dc_csg4.h"
#include "dc_generator_mesh_sdf.h"
#include <cstdio>
#include <cstdint>
#include <cstring>
#include <cstdlib>
#include <vector>
#include <algorithm>
#include <cmath>
#include <random>
#include <map>
#include <tuple>
struct DumpBrush {
float cx, cy, cz, radius, strength;
uint32_t material;
std::vector<ConvexPlane> planes;
};
// ---------------------------------------------------------------- dump reader
static bool read_dump(const char* path, std::vector<DumpBrush>& out, float& ws,
float ch[3], float cc[3]) {
FILE* f = fopen(path, "rb");
if (!f) { fprintf(stderr, "cannot open %s\n", path); return false; }
char magic[4];
if (fread(magic, 1, 4, f) != 4 || memcmp(magic, "DCOR", 4) != 0) {
fprintf(stderr, "bad magic\n"); fclose(f); return false;
}
uint32_t version = 0;
if (fread(&version, 4, 1, f) != 1) { fclose(f); return false; }
if (fread(&ws, 4, 1, f) != 1) { fclose(f); return false; }
if (fread(ch, 4, 3, f) != 3) { fclose(f); return false; }
if (fread(cc, 4, 3, f) != 3) { fclose(f); return false; }
uint32_t n = 0;
if (fread(&n, 4, 1, f) != 1) { fclose(f); return false; }
out.resize(n);
for (uint32_t i = 0; i < n; i++) {
DumpBrush& b = out[i];
if (fread(&b.cx, 4, 1, f) != 1) { fclose(f); return false; }
if (fread(&b.cy, 4, 1, f) != 1) { fclose(f); return false; }
if (fread(&b.cz, 4, 1, f) != 1) { fclose(f); return false; }
if (fread(&b.radius, 4, 1, f) != 1) { fclose(f); return false; }
if (fread(&b.strength, 4, 1, f) != 1) { fclose(f); return false; }
if (fread(&b.material, 4, 1, f) != 1) { fclose(f); return false; }
uint32_t pc = 0;
if (fread(&pc, 4, 1, f) != 1) { fclose(f); return false; }
b.planes.resize(pc);
for (uint32_t p = 0; p < pc; p++) {
if (fread(&b.planes[p].nx, 4, 1, f) != 1) { fclose(f); return false; }
if (fread(&b.planes[p].ny, 4, 1, f) != 1) { fclose(f); return false; }
if (fread(&b.planes[p].nz, 4, 1, f) != 1) { fclose(f); return false; }
if (fread(&b.planes[p].d, 4, 1, f) != 1) { fclose(f); return false; }
uint32_t proj = 0;
if (fread(&proj, 4, 1, f) != 1) { fclose(f); return false; }
}
}
fclose(f);
return true;
}
// ------------------------------------------------------------------ the grid
enum : uint8_t { AIR = 0, SOLID = 1, EXTERIOR = 2 };
struct Grid {
int nx = 0, ny = 0, nz = 0;
float lo[3] = {0,0,0};
float cell = 1.0f;
std::vector<uint8_t> g;
long long solid = 0, exterior = 0, interior = 0;
inline uint8_t at(int x, int y, int z) const {
if (x < 0 || y < 0 || z < 0 || x >= nx || y >= ny || z >= nz) return EXTERIOR;
return g[((size_t)z * ny + y) * nx + x];
}
};
// Overlap, not containment — see the header. The box corner furthest along -n is
// the support point: if even that is outside a plane, no part of the box is in.
static inline bool overlaps_cell(const DumpBrush& b, float x0, float y0, float z0, float c) {
const float x1 = x0 + c, y1 = y0 + c, z1 = z0 + c;
for (const ConvexPlane& p : b.planes) {
const float sx = (p.nx > 0.0f) ? x0 : x1;
const float sy = (p.ny > 0.0f) ? y0 : y1;
const float sz = (p.nz > 0.0f) ? z0 : z1;
if (p.nx * sx + p.ny * sy + p.nz * sz - p.d > 0.0f) return false;
}
return true;
}
static Grid build_grid(const std::vector<DumpBrush>& dump,
const std::vector<CSG4Brush>& cb, float cell) {
Grid r;
r.cell = cell;
float lo[3] = { 1e30f, 1e30f, 1e30f }, hi[3] = { -1e30f, -1e30f, -1e30f };
for (const CSG4Brush& b : cb) {
lo[0] = std::min(lo[0], b.lo.x); hi[0] = std::max(hi[0], b.hi.x);
lo[1] = std::min(lo[1], b.lo.y); hi[1] = std::max(hi[1], b.hi.y);
lo[2] = std::min(lo[2], b.lo.z); hi[2] = std::max(hi[2], b.hi.z);
}
for (int a = 0; a < 3; a++) { lo[a] -= cell * 2.0f; hi[a] += cell * 2.0f; }
r.nx = (int)std::ceil((hi[0] - lo[0]) / cell);
r.ny = (int)std::ceil((hi[1] - lo[1]) / cell);
r.nz = (int)std::ceil((hi[2] - lo[2]) / cell);
r.lo[0] = lo[0]; r.lo[1] = lo[1]; r.lo[2] = lo[2];
r.g.assign((size_t)r.nx * r.ny * r.nz, AIR);
auto idx = [&](int x, int y, int z) { return ((size_t)z * r.ny + y) * r.nx + x; };
for (size_t bi = 0; bi < dump.size(); bi++) {
const CSG4Brush& cbb = cb[bi];
int x0 = std::max(0, (int)std::floor((cbb.lo.x - lo[0]) / cell));
int y0 = std::max(0, (int)std::floor((cbb.lo.y - lo[1]) / cell));
int z0 = std::max(0, (int)std::floor((cbb.lo.z - lo[2]) / cell));
int x1 = std::min(r.nx - 1, (int)std::ceil((cbb.hi.x - lo[0]) / cell));
int y1 = std::min(r.ny - 1, (int)std::ceil((cbb.hi.y - lo[1]) / cell));
int z1 = std::min(r.nz - 1, (int)std::ceil((cbb.hi.z - lo[2]) / cell));
for (int z = z0; z <= z1; z++)
for (int y = y0; y <= y1; y++)
for (int x = x0; x <= x1; x++) {
size_t i = idx(x, y, z);
if (r.g[i] == SOLID) continue;
if (overlaps_cell(dump[bi], lo[0] + x*cell, lo[1] + y*cell, lo[2] + z*cell, cell))
r.g[i] = SOLID;
}
}
std::vector<int> st;
st.reserve(1 << 20);
auto push = [&](int x, int y, int z) {
if (x < 0 || y < 0 || z < 0 || x >= r.nx || y >= r.ny || z >= r.nz) return;
size_t i = idx(x, y, z);
if (r.g[i] != AIR) return;
r.g[i] = EXTERIOR;
st.push_back(x); st.push_back(y); st.push_back(z);
};
for (int z = 0; z < r.nz; z++)
for (int y = 0; y < r.ny; y++)
for (int x = 0; x < r.nx; x++)
if (x == 0 || y == 0 || z == 0 || x == r.nx-1 || y == r.ny-1 || z == r.nz-1)
push(x, y, z);
while (!st.empty()) {
int z = st.back(); st.pop_back();
int y = st.back(); st.pop_back();
int x = st.back(); st.pop_back();
push(x+1,y,z); push(x-1,y,z); push(x,y+1,z); push(x,y-1,z); push(x,y,z+1); push(x,y,z-1);
}
for (size_t i = 0; i < r.g.size(); i++) {
if (r.g[i] == SOLID) r.solid++;
else if (r.g[i] == EXTERIOR) r.exterior++;
else r.interior++;
}
return r;
}
// ------------------------------------------------------------------- the shell
int main(int argc, char** argv) {
const char* path = (argc > 1) ? argv[1] : "/tmp/startmap_brushes_sky.bin";
std::vector<DumpBrush> dump;
float ws = 0, ch[3] = {0,0,0}, cc[3] = {0,0,0};
if (!read_dump(path, dump, ws, ch, cc)) return 1;
std::vector<CSG4Brush> brushes(dump.size());
for (size_t i = 0; i < dump.size(); i++) {
brushes[i].planes = dump[i].planes.data();
brushes[i].plane_count = (uint32_t)dump[i].planes.size();
brushes[i].material_id = dump[i].material;
brushes[i].id = (uint32_t)i;
}
dc_csg4::compute_bounds(brushes);
std::vector<CSG4Face> skin;
CSG4Stats stats{};
dc_csg4::clip_all(brushes, skin, &stats);
printf("dump %s\n", path);
printf(" brushes %zu\n", dump.size());
printf(" csg4 skin: raw %u -> out %u faces (removed %u, fragmented %u, max frags %u)\n",
stats.raw_faces, stats.out_faces, stats.removed, stats.fragmented, stats.max_fragments);
printf(" skin area %.0f m^2 (raw %.0f)\n", stats.out_area, stats.raw_area);
// EXACT REFERENCE, no voxels involved. The csg4 skin is the closed boundary
// of the brush union, so its divergence volume IS the level's solid volume.
// The voxel denominator below is biased HIGH — conservative rasterization
// fattens SOLID by up to a cell — so this is the number to trust, and
// shell_enclosed - skin_solid = the interior air the shell wraps
// is a relation between two exact quantities.
double skin_vol = 0.0, skin_area_chk = 0.0;
for (const CSG4Face& f : skin) {
if (f.pts.size() < 3) continue;
const ConvexPlane& p = dump[f.brush_id].planes[f.face_index];
double cx = 0, cy = 0, cz = 0;
for (const CSG4Vec3& v : f.pts) { cx += v.x; cy += v.y; cz += v.z; }
cx /= f.pts.size(); cy /= f.pts.size(); cz /= f.pts.size();
double area = 0.0;
for (size_t i = 1; i + 1 < f.pts.size(); i++) {
double ux = f.pts[i].x - f.pts[0].x, uy = f.pts[i].y - f.pts[0].y, uz = f.pts[i].z - f.pts[0].z;
double vx = f.pts[i+1].x - f.pts[0].x, vy = f.pts[i+1].y - f.pts[0].y, vz = f.pts[i+1].z - f.pts[0].z;
double ax = uy*vz - uz*vy, ay = uz*vx - ux*vz, az = ux*vy - uy*vx;
area += 0.5 * std::sqrt(ax*ax + ay*ay + az*az);
}
skin_area_chk += area;
skin_vol += (cx*p.nx + cy*p.ny + cz*p.nz) * area / 3.0;
}
printf(" FULL SKIN (closed by construction): area %.0f m^2 volume %.0f m^3\n",
skin_area_chk, skin_vol);
printf(" -> the level's SOLID volume; shell_encl minus this = the air it wraps\n\n");
// OFFSET SWEEP. Conservative rasterization marks a cell SOLID if a brush
// overlaps it AT ALL, which inflates the solid by up to one cell in every
// direction. A sample 0.75 cells outside a face therefore often lands in a
// cell the face's OWN brush still overlaps -> reported SOLID -> the face is
// rejected, systematically, at every wall. If that is the defect the ratio
// must climb with the offset and then flatten; if it does not, the missing
// surface is something else and the offset is not worth tuning.
printf("%-8s %8s %10s %10s %8s %12s %14s %14s %8s\n",
"cell", "offset", "skin", "shell", "shell%", "shell m^2", "encl m^3", "voxel m^3", "ratio");
for (float cell : { 1.0f, 0.5f, 0.25f }) {
Grid grid = build_grid(dump, brushes, cell);
for (float mult : { 0.75f, 1.5f, 2.5f, 4.0f }) {
long long shell = 0;
double shell_area = 0.0, encl = 0.0;
for (const CSG4Face& f : skin) {
if (f.pts.size() < 3) continue;
const ConvexPlane& p = dump[f.brush_id].planes[f.face_index];
double cx = 0, cy = 0, cz = 0;
for (const CSG4Vec3& v : f.pts) { cx += v.x; cy += v.y; cz += v.z; }
cx /= f.pts.size(); cy /= f.pts.size(); cz /= f.pts.size();
double area = 0.0;
for (size_t i = 1; i + 1 < f.pts.size(); i++) {
double ux = f.pts[i].x - f.pts[0].x, uy = f.pts[i].y - f.pts[0].y, uz = f.pts[i].z - f.pts[0].z;
double vx = f.pts[i+1].x - f.pts[0].x, vy = f.pts[i+1].y - f.pts[0].y, vz = f.pts[i+1].z - f.pts[0].z;
double ax = uy*vz - uz*vy, ay = uz*vx - ux*vz, az = ux*vy - uy*vx;
area += 0.5 * std::sqrt(ax*ax + ay*ay + az*az);
}
if (area <= 0.0) continue;
double sx = cx + p.nx * cell * mult;
double sy = cy + p.ny * cell * mult;
double sz = cz + p.nz * cell * mult;
int gx = (int)std::floor((sx - grid.lo[0]) / cell);
int gy = (int)std::floor((sy - grid.lo[1]) / cell);
int gz = (int)std::floor((sz - grid.lo[2]) / cell);
if (grid.at(gx, gy, gz) != EXTERIOR) continue;
shell++;
shell_area += area;
encl += (cx*p.nx + cy*p.ny + cz*p.nz) * area / 3.0;
}
double voxel = (double)(grid.solid + grid.interior) * cell * cell * cell;
printf("%-8.3f %8.2f %10zu %10lld %7.1f%% %12.0f %14.0f %14.0f %8.3f\n",
cell, mult, skin.size(), shell,
100.0 * shell / (double)std::max<size_t>(1, skin.size()),
shell_area, encl, voxel, encl / std::max(1.0, voxel));
fflush(stdout);
}
printf("\n");
}
// ---- emit the shell as OBJ at the settled parameters ----
// cell 0.25 / offset 1.5: the offset sweep shows the ratio plateaus from 1.5
// onward, and a probe under one cell lies at every wall.
{
const float cell = 0.25f, mult = 1.5f;
Grid grid = build_grid(dump, brushes, cell);
std::vector<CSG4Vec3> verts;
std::vector<uint32_t> idx;
long long tris = 0;
for (const CSG4Face& f : skin) {
if (f.pts.size() < 3) continue;
const ConvexPlane& p = dump[f.brush_id].planes[f.face_index];
double cx = 0, cy = 0, cz = 0;
for (const CSG4Vec3& v : f.pts) { cx += v.x; cy += v.y; cz += v.z; }
cx /= f.pts.size(); cy /= f.pts.size(); cz /= f.pts.size();
double sx = cx + p.nx * cell * mult;
double sy = cy + p.ny * cell * mult;
double sz = cz + p.nz * cell * mult;
int gx = (int)std::floor((sx - grid.lo[0]) / cell);
int gy = (int)std::floor((sy - grid.lo[1]) / cell);
int gz = (int)std::floor((sz - grid.lo[2]) / cell);
if (grid.at(gx, gy, gz) != EXTERIOR) continue;
// Fan-triangulate the convex fragment, keeping csg4's winding so the
// mesh SDF's inside/outside test agrees with the plane normals.
uint32_t base = (uint32_t)verts.size();
for (const CSG4Vec3& v : f.pts) verts.push_back(v);
for (size_t i = 1; i + 1 < f.pts.size(); i++) {
idx.push_back(base);
idx.push_back(base + (uint32_t)i);
idx.push_back(base + (uint32_t)i + 1);
tris++;
}
}
const char* obj = "/tmp/outer_shell.obj";
FILE* o = fopen(obj, "w");
if (o) {
fprintf(o, "# START.MAP outer shell — cell %.3f offset %.2f\n", cell, mult);
for (const CSG4Vec3& v : verts) fprintf(o, "v %.6f %.6f %.6f\n", v.x, v.y, v.z);
for (size_t i = 0; i + 2 < idx.size(); i += 3)
fprintf(o, "f %u %u %u\n", idx[i]+1, idx[i+1]+1, idx[i+2]+1);
fclose(o);
}
printf("\n-- SHELL MESH (cell %.2f, offset %.2f) --\n", cell, mult);
printf(" vertices %zu triangles %lld -> %s\n", verts.size(), tris, obj);
printf(" BVH cost at 128 B/tri: %.0f KB (buffers resize, no cap)\n",
tris * 128.0 / 1024.0);
}
// ---- OPTION 2: the flood-fill boundary, watertight BY CONSTRUCTION ----
//
// The csg4 skin is polygon soup (T-junctions from independent per-face
// clipping), so a pseudonormal mesh SDF cannot use it. The boundary of a
// VOXEL SET has no such problem: emit a quad wherever a cell inside the
// sealed volume touches an EXTERIOR cell, key the corners on the INTEGER
// LATTICE so welding is exact rather than epsilon-based, and the result is
// closed and manifold with nothing to repair.
//
// The cost is that it is stair-stepped at the cell size. Snapping those
// vertices onto the brush planes is the next step and is deliberately NOT
// done here — watertightness is proven first, on its own, so that if snapping
// later breaks it the edge count says so immediately.
{
const float cell = 0.25f;
Grid grid = build_grid(dump, brushes, cell);
// Lattice corner -> vertex index. (nx+1)*(ny+1)*(nz+1) corners.
const int cx1 = grid.nx + 1, cy1 = grid.ny + 1;
std::vector<int> vid((size_t)(grid.nx+1) * (grid.ny+1) * (grid.nz+1), -1);
std::vector<CSG4Vec3> verts;
std::vector<uint32_t> idx;
auto corner = [&](int x, int y, int z) -> uint32_t {
size_t k = ((size_t)z * cy1 + y) * cx1 + x;
if (vid[k] < 0) {
vid[k] = (int)verts.size();
verts.push_back({ grid.lo[0] + x * cell,
grid.lo[1] + y * cell,
grid.lo[2] + z * cell });
}
return (uint32_t)vid[k];
};
// WINDING IS REVERSED from the obvious choice, and this was found by the
// sign check rather than by reasoning. Emitting a,b,c / a,c,d gives a
// divergence volume of +38805, which by the usual convention means
// outward-facing — yet DCGeneratorMeshSDF then reported INSIDE for the
// open void and OUTSIDE for the sealed rooms, i.e. exactly inverted
// (exterior 1.2% correct, interior 3.1%). The two disagree about
// handedness, and no amount of surface statistics would have shown it:
// area, closure and manifoldness are all invariant under a global flip.
// The same trap as hyper-games #9, where 0.768.0 shipped an inside-out
// level with a fully passing gate.
auto quad = [&](uint32_t a, uint32_t b, uint32_t c, uint32_t d) {
idx.push_back(a); idx.push_back(c); idx.push_back(b);
idx.push_back(a); idx.push_back(d); idx.push_back(c);
};
auto inside = [&](int x, int y, int z) { return grid.at(x,y,z) != EXTERIOR; };
long long quads = 0;
for (int z = 0; z < grid.nz; z++)
for (int y = 0; y < grid.ny; y++)
for (int x = 0; x < grid.nx; x++) {
if (!inside(x, y, z)) continue;
// Winding chosen so the normal points OUT of the sealed volume; the
// divergence volume below is the check that it does.
if (!inside(x+1,y,z)) { quad(corner(x+1,y,z), corner(x+1,y+1,z), corner(x+1,y+1,z+1), corner(x+1,y,z+1)); quads++; }
if (!inside(x-1,y,z)) { quad(corner(x,y,z), corner(x,y,z+1), corner(x,y+1,z+1), corner(x,y+1,z)); quads++; }
if (!inside(x,y+1,z)) { quad(corner(x,y+1,z), corner(x,y+1,z+1), corner(x+1,y+1,z+1), corner(x+1,y+1,z)); quads++; }
if (!inside(x,y-1,z)) { quad(corner(x,y,z), corner(x+1,y,z), corner(x+1,y,z+1), corner(x,y,z+1)); quads++; }
if (!inside(x,y,z+1)) { quad(corner(x,y,z+1), corner(x+1,y,z+1), corner(x+1,y+1,z+1), corner(x,y+1,z+1)); quads++; }
if (!inside(x,y,z-1)) { quad(corner(x,y,z), corner(x,y+1,z), corner(x+1,y+1,z), corner(x+1,y,z)); quads++; }
}
// Divergence volume over the emitted triangles.
double vol = 0.0, area = 0.0;
for (size_t i = 0; i + 2 < idx.size(); i += 3) {
const CSG4Vec3& A = verts[idx[i]];
const CSG4Vec3& B = verts[idx[i+1]];
const CSG4Vec3& C = verts[idx[i+2]];
double ux = B.x-A.x, uy = B.y-A.y, uz = B.z-A.z;
double vx = C.x-A.x, vy = C.y-A.y, vz = C.z-A.z;
double nx = uy*vz - uz*vy, ny = uz*vx - ux*vz, nz = ux*vy - uy*vx;
double a = 0.5 * std::sqrt(nx*nx + ny*ny + nz*nz);
if (a <= 0.0) continue;
area += a;
double cxx = (A.x+B.x+C.x)/3.0, cyy = (A.y+B.y+C.y)/3.0, czz = (A.z+B.z+C.z)/3.0;
vol += (cxx*nx + cyy*ny + czz*nz) / 6.0; // n is unnormalised = 2*area*unit
}
double voxel_vol = (double)(grid.solid + grid.interior) * cell * cell * cell;
printf("\n-- OPTION 2: flood-fill boundary (watertight by construction) --\n");
printf(" cell %.3f quads %lld verts %zu tris %zu\n",
cell, quads, verts.size(), idx.size()/3);
printf(" area %.0f m^2 divergence volume %.0f m^3 voxel volume %.0f m^3 ratio %.4f\n",
area, vol, voxel_vol, vol / std::max(1.0, voxel_vol));
// ---- SNAP: lattice corners onto the brush planes ----
//
// Connectivity is fixed by the lattice, so this CANNOT open the surface —
// watertightness survives by construction. What it can break is geometry:
// degenerate triangles, flipped normals, and the enclosed volume. All
// three are measured below rather than hoped for.
//
// The rule is #22 §2.5's: an element of dimension k is the meet of 3-k
// planes. Gather the planes this vertex is near, then
// 1 plane -> project onto it
// 2 planes -> project onto their line of intersection
// 3 planes -> the point where they meet
//
// A PLANE IS INFINITE AND A FACE IS NOT. Only planes of brushes whose
// expanded AABB contains the vertex are considered, or a vertex would
// snap to some wall's plane extended across the map.
const std::vector<CSG4Vec3> pre = verts;
const float tol = cell * 1.5f; // within 1.5 cells of a plane
printf(" %-7s %9s %9s %9s %10s %9s %9s %9s %12s\n",
"clamp","snapped","unmoved","flipped","dropped","tris","bound","nonman","volume m^3");
for (float clampx : { 1.50f, 1.25f, 1.00f, 0.75f, 0.50f }) {
verts = pre;
long long snapped_1 = 0, snapped_2 = 0, snapped_3 = 0, unmoved = 0;
double max_move = 0.0, sum_move = 0.0;
for (size_t vi2 = 0; vi2 < verts.size(); vi2++) {
const CSG4Vec3 v = verts[vi2];
// Candidate planes: near this vertex AND belonging to a brush whose
// expanded bounds contain it.
float pn[3][3]; float pd[3]; int np = 0;
for (size_t bi = 0; bi < brushes.size() && np < 3; bi++) {
const CSG4Brush& bb = brushes[bi];
if (v.x < bb.lo.x - tol || v.x > bb.hi.x + tol) continue;
if (v.y < bb.lo.y - tol || v.y > bb.hi.y + tol) continue;
if (v.z < bb.lo.z - tol || v.z > bb.hi.z + tol) continue;
for (const ConvexPlane& pl : dump[bi].planes) {
float dist = pl.nx*v.x + pl.ny*v.y + pl.nz*v.z - pl.d;
if (std::fabs(dist) > tol) continue;
// Skip a plane parallel to one already taken (same surface).
bool dup = false;
for (int k = 0; k < np; k++) {
float dot = pn[k][0]*pl.nx + pn[k][1]*pl.ny + pn[k][2]*pl.nz;
if (std::fabs(dot) > 0.999f) { dup = true; break; }
}
if (dup) continue;
if (np < 3) { pn[np][0]=pl.nx; pn[np][1]=pl.ny; pn[np][2]=pl.nz; pd[np]=pl.d; np++; }
if (np >= 3) break;
}
}
if (np == 0) { unmoved++; continue; }
CSG4Vec3 nv = v;
if (np == 1) {
float dist = pn[0][0]*v.x + pn[0][1]*v.y + pn[0][2]*v.z - pd[0];
nv = { v.x - pn[0][0]*dist, v.y - pn[0][1]*dist, v.z - pn[0][2]*dist };
snapped_1++;
} else if (np == 2) {
// Project onto the line where the two planes meet: subtract the
// component of the offset lying in the plane spanned by the normals.
float a11=1.0f, a12=pn[0][0]*pn[1][0]+pn[0][1]*pn[1][1]+pn[0][2]*pn[1][2];
float d0 = pn[0][0]*v.x + pn[0][1]*v.y + pn[0][2]*v.z - pd[0];
float d1 = pn[1][0]*v.x + pn[1][1]*v.y + pn[1][2]*v.z - pd[1];
float det = a11*a11 - a12*a12;
if (std::fabs(det) > 1e-6f) {
float u = (d0*a11 - d1*a12) / det;
float w = (d1*a11 - d0*a12) / det;
nv = { v.x - pn[0][0]*u - pn[1][0]*w,
v.y - pn[0][1]*u - pn[1][1]*w,
v.z - pn[0][2]*u - pn[1][2]*w };
snapped_2++;
} else { unmoved++; continue; }
} else {
// Three planes meet at a point: solve the 3x3 by Cramer.
float m[3][3] = {{pn[0][0],pn[0][1],pn[0][2]},
{pn[1][0],pn[1][1],pn[1][2]},
{pn[2][0],pn[2][1],pn[2][2]}};
float b[3] = { pd[0], pd[1], pd[2] };
float det = m[0][0]*(m[1][1]*m[2][2]-m[1][2]*m[2][1])
- m[0][1]*(m[1][0]*m[2][2]-m[1][2]*m[2][0])
+ m[0][2]*(m[1][0]*m[2][1]-m[1][1]*m[2][0]);
if (std::fabs(det) > 1e-6f) {
float ix = (b[0]*(m[1][1]*m[2][2]-m[1][2]*m[2][1])
- m[0][1]*(b[1]*m[2][2]-m[1][2]*b[2])
+ m[0][2]*(b[1]*m[2][1]-m[1][1]*b[2])) / det;
float iy = (m[0][0]*(b[1]*m[2][2]-m[1][2]*b[2])
- b[0]*(m[1][0]*m[2][2]-m[1][2]*m[2][0])
+ m[0][2]*(m[1][0]*b[2]-b[1]*m[2][0])) / det;
float iz = (m[0][0]*(m[1][1]*b[2]-b[1]*m[2][1])
- m[0][1]*(m[1][0]*b[2]-b[1]*m[2][0])
+ b[0]*(m[1][0]*m[2][1]-m[1][1]*m[2][0])) / det;
nv = { ix, iy, iz };
snapped_3++;
} else { unmoved++; continue; }
}
// CLAMP. A vertex that moves further than a cell has found a plane it
// does not belong to; keeping it would tear the surface across itself.
double mv = std::sqrt((nv.x-v.x)*(nv.x-v.x) + (nv.y-v.y)*(nv.y-v.y) + (nv.z-v.z)*(nv.z-v.z));
if (mv > cell * clampx) { unmoved++; continue; }
max_move = std::max(max_move, mv);
sum_move += mv;
verts[vi2] = nv;
}
// WELD then DROP. A collapsed stair riser merges two lattice corners onto
// the same point, so degenerate triangles can only be removed safely after
// re-welding by POSITION — dropping them off the lattice indices alone
// would strand edges and open the surface. Then re-audit the edges,
// because that is the property the whole approach exists to preserve.
std::map<std::tuple<long long,long long,long long>, uint32_t> wmap;
std::vector<uint32_t> wid(verts.size());
const double Q = 1e4;
for (size_t i = 0; i < verts.size(); i++) {
auto k = std::make_tuple((long long)llround(verts[i].x*Q),
(long long)llround(verts[i].y*Q),
(long long)llround(verts[i].z*Q));
auto it = wmap.find(k);
if (it == wmap.end()) { uint32_t n2 = (uint32_t)wmap.size(); wmap[k] = n2; wid[i] = n2; }
else wid[i] = it->second;
}
double vol2 = 0.0, area2 = 0.0; long long degen = 0, flipped = 0, kept = 0;
std::map<std::pair<uint32_t,uint32_t>, int> ecount;
for (size_t i = 0; i + 2 < idx.size(); i += 3) {
uint32_t ia = wid[idx[i]], ib = wid[idx[i+1]], ic = wid[idx[i+2]];
if (ia == ib || ib == ic || ia == ic) { degen++; continue; }
const CSG4Vec3& A = verts[idx[i]]; const CSG4Vec3& B = verts[idx[i+1]]; const CSG4Vec3& C = verts[idx[i+2]];
double ux=B.x-A.x, uy=B.y-A.y, uz=B.z-A.z, vx=C.x-A.x, vy=C.y-A.y, vz=C.z-A.z;
double nx=uy*vz-uz*vy, ny=uz*vx-ux*vz, nz=ux*vy-uy*vx;
double a = 0.5*std::sqrt(nx*nx+ny*ny+nz*nz);
if (a <= 1e-12) { degen++; continue; }
const CSG4Vec3& A0 = pre[idx[i]]; const CSG4Vec3& B0 = pre[idx[i+1]]; const CSG4Vec3& C0 = pre[idx[i+2]];
double u0x=B0.x-A0.x,u0y=B0.y-A0.y,u0z=B0.z-A0.z, v0x=C0.x-A0.x,v0y=C0.y-A0.y,v0z=C0.z-A0.z;
double m0x=u0y*v0z-u0z*v0y, m0y=u0z*v0x-u0x*v0z, m0z=u0x*v0y-u0y*v0x;
if (nx*m0x + ny*m0y + nz*m0z < 0.0) flipped++;
kept++;
area2 += a;
vol2 += ((A.x+B.x+C.x)/3.0*nx + (A.y+B.y+C.y)/3.0*ny + (A.z+B.z+C.z)/3.0*nz) / 6.0;
uint32_t t[3] = { ia, ib, ic };
for (int e = 0; e < 3; e++) {
uint32_t u = t[e], v2 = t[(e+1)%3];
ecount[{std::min(u,v2), std::max(u,v2)}]++;
}
}
long long bnd = 0, nonman = 0;
for (auto& kv : ecount) { if (kv.second == 1) bnd++; else if (kv.second > 2) nonman++; }
printf(" %-7.2f %9lld %9lld %9lld %10lld %9lld %9lld %9lld %12.0f\n",
clampx, snapped_1+snapped_2+snapped_3, unmoved, flipped, degen, kept, bnd, nonman, vol2);
fflush(stdout);
}
printf(" (exact solid 12600 + estimated air ~23000 => target ~35600; pre-snap %.0f)\n", vol);
FILE* o = fopen("/tmp/shell_voxel.obj", "w");
if (o) {
for (const CSG4Vec3& v : verts) fprintf(o, "v %.6f %.6f %.6f\n", v.x, v.y, v.z);
for (size_t i = 0; i + 2 < idx.size(); i += 3)
fprintf(o, "f %u %u %u\n", idx[i]+1, idx[i+1]+1, idx[i+2]+1);
fclose(o);
printf(" -> /tmp/shell_voxel.obj (snapped at clamp 1.5, degenerates KEPT)\n");
}
// ---- SIGN VALIDATION: the only test that checks what we need ----
//
// Every mesh statistic so far measures the SURFACE. What the subtract
// actually needs is that the SDF says INSIDE where the level is sealed and
// OUTSIDE where it is not. The flood fill already labelled every cell, so
// that is a direct comparison against an independent answer rather than
// against itself.
//
// It also converts the outstanding degenerate-normal risk from a theory
// into an observation: if a zero-area triangle is ever chosen as the
// closest primitive and hands back a zero normal, it shows up here as a
// wrong sign.
//
// Convention (dc_generator_mesh_sdf.h:15): positive = solid (inside mesh).
{
DCGeneratorMeshSDF sdf;
if (!sdf.load_obj("/tmp/shell_voxel.obj")) {
printf(" SIGN CHECK: could not load the mesh back\n");
} else {
std::mt19937 rng(4242);
long long ok_in = 0, bad_in = 0, ok_out = 0, bad_out = 0, ok_solid = 0, bad_solid = 0;
double bad_in_d = 0, bad_out_d = 0, bad_solid_d = 0;
long long bad_out_far = 0, bad_solid_far = 0;
std::uniform_int_distribution<int> rx(1, grid.nx-2), ry(1, grid.ny-2), rz(1, grid.nz-2);
int want = 20000, tries = 0;
long long got_in = 0, got_out = 0, got_solid = 0;
while (tries < 4000000 && (got_in < want || got_out < want || got_solid < want)) {
tries++;
int x = rx(rng), y = ry(rng), z = rz(rng);
uint8_t cls = grid.at(x, y, z);
float px = grid.lo[0] + (x + 0.5f) * cell;
float py = grid.lo[1] + (y + 0.5f) * cell;
float pz = grid.lo[2] + (z + 0.5f) * cell;
if (cls == AIR && got_in < want) { // sealed interior
got_in++;
float d2 = sdf.evaluate(px, py, pz, cell);
if (d2 > 0.0f) ok_in++; else { bad_in++; bad_in_d += std::fabs(d2); }
} else if (cls == EXTERIOR && got_out < want) {
got_out++;
float d2 = sdf.evaluate(px, py, pz, cell);
if (d2 < 0.0f) ok_out++; else { bad_out++; bad_out_d += std::fabs(d2);
if (std::fabs(d2) > cell * 2.0f) bad_out_far++; }
} else if (cls == SOLID && got_solid < want) {
got_solid++;
float d2 = sdf.evaluate(px, py, pz, cell);
if (d2 > 0.0f) ok_solid++; else { bad_solid++; bad_solid_d += std::fabs(d2);
if (std::fabs(d2) > cell * 2.0f) bad_solid_far++; }
}
}
printf("\n-- SIGN CHECK vs the flood fill's own labels --\n");
printf(" interior (sealed air) %7lld ok %7lld wrong %6.2f%%\n",
ok_in, bad_in, 100.0 * ok_in / std::max(1LL, ok_in + bad_in));
printf(" solid (inside brush) %7lld ok %7lld wrong %6.2f%%\n",
ok_solid, bad_solid, 100.0 * ok_solid / std::max(1LL, ok_solid + bad_solid));
printf(" exterior (open void) %7lld ok %7lld wrong %6.2f%%\n",
ok_out, bad_out, 100.0 * ok_out / std::max(1LL, ok_out + bad_out));
// A wrong sample sitting ~half a cell from the surface is the grid
// and the snapped mesh disagreeing about where the boundary is,
// which is expected — snapping MOVED it by up to 1.5 cells. A wrong
// sample far from the surface is a real defect. The two are
// indistinguishable in a percentage, so separate them.
printf(" mean |sdf| of wrong samples: interior %.3f m solid %.3f m exterior %.3f m (cell %.3f)\n",
bad_in_d / std::max(1LL, bad_in), bad_solid_d / std::max(1LL, bad_solid),
bad_out_d / std::max(1LL, bad_out), cell);
printf(" wrong by MORE than 2 cells: solid %lld/%lld exterior %lld/%lld <- the real defects\n",
bad_solid_far, bad_solid, bad_out_far, bad_out);
}
}
}
printf("\nclosure: a closed, outward-wound shell gives ratio 1.0.\n");
printf("below 1.0 = holes (missing faces); above = double-counted surface.\n");
return 0;
}