#include "dc_level_shell.h"
#include <algorithm>
#include <cmath>
#include <unordered_map>
namespace dc_level_shell {
namespace {
// CLEAR is a fourth state and not a shade of SOLID, for one reason: sealed() is
// interior_volume > solid_volume * 0.5, and clearance is AIR THE LEVEL OWNS.
// Counted as solid it inflates the wall term, so a level that is perfectly
// sealed reports otherwise and drops to the moat fallback — and it would do that
// at the exact moment someone built a generous enough passage, which is the
// worst possible time to find it.
//
// Everything downstream then falls out with no further change: the flood
// advances only through AIR so CLEAR blocks it, the cell census counts anything
// that is neither SOLID nor EXTERIOR as interior, and `inside` is != EXTERIOR so
// the boundary is emitted around it.
enum : uint8_t { AIR = 0, SOLID = 1, EXTERIOR = 2, CLEAR = 3 };
// Does the brush overlap the CELL, rather than contain its centre? Centre
// sampling misses Quake's thin walls and the exterior flood pours through them —
// 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 inside the brush.
inline bool overlaps_cell(const CSG4Brush& b, float x0, float y0, float z0, float c) {
const float x1 = x0 + c, y1 = y0 + c, z1 = z0 + c;
for (uint32_t i = 0; i < b.plane_count; i++) {
const ConvexPlane& p = b.planes[i];
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;
}
} // namespace
bool build(const std::vector<CSG4Brush>& brushes, ShellMesh& out,
float cell, float clamp_cells) {
out = ShellMesh{};
if (brushes.empty() || cell <= 0.0f) return false;
// ---- bounds, padded so the fill always has an outside to start from ----
float lo[3] = { 1e30f, 1e30f, 1e30f }, hi[3] = { -1e30f, -1e30f, -1e30f };
for (const CSG4Brush& b : brushes) {
if (!b.bounds_valid) continue;
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);
}
if (lo[0] > hi[0]) return false; // no brush had valid bounds
for (int a = 0; a < 3; a++) { lo[a] -= cell * 2.0f; hi[a] += cell * 2.0f; }
const int nx = (int)std::ceil((hi[0] - lo[0]) / cell);
const int ny = (int)std::ceil((hi[1] - lo[1]) / cell);
const int nz = (int)std::ceil((hi[2] - lo[2]) / cell);
if (nx < 3 || ny < 3 || nz < 3) return false;
out.grid[0] = nx; out.grid[1] = ny; out.grid[2] = nz;
std::vector<uint8_t> g((size_t)nx * ny * nz, AIR);
auto idx3 = [&](int x, int y, int z) -> size_t {
return ((size_t)z * ny + y) * nx + x;
};
auto at = [&](int x, int y, int z) -> uint8_t {
if (x < 0 || y < 0 || z < 0 || x >= nx || y >= ny || z >= nz) return EXTERIOR;
return g[idx3(x, y, z)];
};
// ---- 1. rasterize, each brush over its own bounds only ----
//
// TWO PASSES, SOLID FIRST, so a cell claimed by both real geometry and a
// clearance volume is wall. They overlap by construction — an archway's
// clearance encloses the ring it belongs to — and one pass would decide it by
// whichever brush the store happened to list first, which is import order:
// stable, arbitrary, and different once anything is edited.
auto rasterize = [&](const CSG4Brush& b, uint8_t mark) {
if (!b.bounds_valid) return;
int x0 = std::max(0, (int)std::floor((b.lo.x - lo[0]) / cell));
int y0 = std::max(0, (int)std::floor((b.lo.y - lo[1]) / cell));
int z0 = std::max(0, (int)std::floor((b.lo.z - lo[2]) / cell));
int x1 = std::min(nx - 1, (int)std::ceil((b.hi.x - lo[0]) / cell));
int y1 = std::min(ny - 1, (int)std::ceil((b.hi.y - lo[1]) / cell));
int z1 = std::min(nz - 1, (int)std::ceil((b.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 = idx3(x, y, z);
if (g[i] != AIR) continue;
if (overlaps_cell(b, lo[0] + x * cell, lo[1] + y * cell, lo[2] + z * cell, cell))
g[i] = mark;
}
};
for (const CSG4Brush& b : brushes) if (!b.clearance) rasterize(b, SOLID);
for (const CSG4Brush& b : brushes) if (b.clearance) rasterize(b, CLEAR);
// ---- 2. flood fill EXTERIOR from the border (qbsp's own leak test) ----
std::vector<int> st;
st.reserve(1 << 16);
auto push = [&](int x, int y, int z) {
if (x < 0 || y < 0 || z < 0 || x >= nx || y >= ny || z >= nz) return;
size_t i = idx3(x, y, z);
if (g[i] != AIR) return;
g[i] = EXTERIOR;
st.push_back(x); st.push_back(y); st.push_back(z);
};
for (int z = 0; z < nz; z++)
for (int y = 0; y < ny; y++)
for (int x = 0; x < nx; x++)
if (x == 0 || y == 0 || z == 0 || x == nx - 1 || y == ny - 1 || z == 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 < g.size(); i++) {
if (g[i] == SOLID) out.cells_solid++;
else if (g[i] == EXTERIOR) out.cells_exterior++;
else out.cells_interior++;
}
// ---- 3. emit the boundary, corners keyed on the integer lattice ----
const int cx1 = nx + 1, cy1 = ny + 1;
std::vector<int> vid((size_t)(nx + 1) * (ny + 1) * (nz + 1), -1);
auto corner = [&](int x, int y, int z) -> int {
size_t k = ((size_t)z * cy1 + y) * cx1 + x;
if (vid[k] < 0) {
vid[k] = (int)(out.vertices.size() / 3);
out.vertices.push_back(lo[0] + x * cell);
out.vertices.push_back(lo[1] + y * cell);
out.vertices.push_back(lo[2] + z * cell);
}
return vid[k];
};
// Reversed from the obvious order — see the header. This winding is what
// makes the SDF report inside for the sealed volume.
auto quad = [&](int a, int b, int c, int d) {
out.indices.push_back(a); out.indices.push_back(c); out.indices.push_back(b);
out.indices.push_back(a); out.indices.push_back(d); out.indices.push_back(c);
};
auto inside = [&](int x, int y, int z) { return at(x, y, z) != EXTERIOR; };
for (int z = 0; z < nz; z++)
for (int y = 0; y < ny; y++)
for (int x = 0; x < nx; x++) {
if (!inside(x, y, z)) continue;
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));
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));
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));
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));
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));
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));
}
// ---- 4. snap corners onto the brush planes ----
// #22 §2.5: an element of dimension k is the meet of 3-k planes. A plane is
// infinite and a face is not, so only planes of brushes whose expanded bounds
// contain the corner are considered — otherwise a corner snaps to some far
// wall's plane extended across the map.
const float tol = cell * 1.5f;
const float maxmove = cell * clamp_cells;
const int nverts = (int)(out.vertices.size() / 3);
for (int vi = 0; vi < nverts; vi++) {
const float vx = out.vertices[vi*3+0];
const float vy = out.vertices[vi*3+1];
const float vz = out.vertices[vi*3+2];
float pn[3][3]; float pd[3]; int np = 0;
for (size_t bi = 0; bi < brushes.size() && np < 3; bi++) {
const CSG4Brush& b = brushes[bi];
if (!b.bounds_valid) continue;
if (vx < b.lo.x - tol || vx > b.hi.x + tol) continue;
if (vy < b.lo.y - tol || vy > b.hi.y + tol) continue;
if (vz < b.lo.z - tol || vz > b.hi.z + tol) continue;
for (uint32_t pi = 0; pi < b.plane_count && np < 3; pi++) {
const ConvexPlane& pl = b.planes[pi];
if (std::fabs(pl.nx*vx + pl.ny*vy + pl.nz*vz - pl.d) > tol) continue;
bool dup = false;
for (int k = 0; k < np; k++)
if (std::fabs(pn[k][0]*pl.nx + pn[k][1]*pl.ny + pn[k][2]*pl.nz) > 0.999f) { dup = true; break; }
if (dup) continue;
pn[np][0] = pl.nx; pn[np][1] = pl.ny; pn[np][2] = pl.nz; pd[np] = pl.d; np++;
}
}
if (np == 0) { out.unmoved++; continue; }
float ox = vx, oy = vy, oz = vz;
if (np == 1) {
float dist = pn[0][0]*vx + pn[0][1]*vy + pn[0][2]*vz - pd[0];
ox = vx - pn[0][0]*dist; oy = vy - pn[0][1]*dist; oz = vz - pn[0][2]*dist;
} else if (np == 2) {
float a12 = pn[0][0]*pn[1][0] + pn[0][1]*pn[1][1] + pn[0][2]*pn[1][2];
float d0 = pn[0][0]*vx + pn[0][1]*vy + pn[0][2]*vz - pd[0];
float d1 = pn[1][0]*vx + pn[1][1]*vy + pn[1][2]*vz - pd[1];
float det = 1.0f - a12 * a12;
if (std::fabs(det) < 1e-6f) { out.unmoved++; continue; }
float u = (d0 - d1 * a12) / det;
float w = (d1 - d0 * a12) / det;
ox = vx - pn[0][0]*u - pn[1][0]*w;
oy = vy - pn[0][1]*u - pn[1][1]*w;
oz = vz - pn[0][2]*u - pn[1][2]*w;
} else {
float m00=pn[0][0],m01=pn[0][1],m02=pn[0][2];
float m10=pn[1][0],m11=pn[1][1],m12=pn[1][2];
float m20=pn[2][0],m21=pn[2][1],m22=pn[2][2];
float det = m00*(m11*m22-m12*m21) - m01*(m10*m22-m12*m20) + m02*(m10*m21-m11*m20);
if (std::fabs(det) < 1e-6f) { out.unmoved++; continue; }
float b0=pd[0], b1=pd[1], b2=pd[2];
ox = (b0*(m11*m22-m12*m21) - m01*(b1*m22-m12*b2) + m02*(b1*m21-m11*b2)) / det;
oy = (m00*(b1*m22-m12*b2) - b0*(m10*m22-m12*m20) + m02*(m10*b2-b1*m20)) / det;
oz = (m00*(m11*b2-b1*m21) - m01*(m10*b2-b1*m20) + b0*(m10*m21-m11*m20)) / det;
}
float dx = ox - vx, dy = oy - vy, dz = oz - vz;
if (std::sqrt(dx*dx + dy*dy + dz*dz) > maxmove) { out.unmoved++; continue; }
out.vertices[vi*3+0] = ox;
out.vertices[vi*3+1] = oy;
out.vertices[vi*3+2] = oz;
out.snapped++;
}
// ---- watertightness, measured on what we are actually about to hand out ----
// Keyed on the WELDED lattice index, which is exact: corners come from the
// integer grid, so two triangles meeting at a corner reference the same
// vertex id with no epsilon involved. Snapping moves positions and never
// touches connectivity, so this is valid after it.
{
std::unordered_map<uint64_t, uint32_t> ecount;
ecount.reserve(out.indices.size());
auto key = [](int a, int b) -> uint64_t {
uint32_t lo2 = (uint32_t)(a < b ? a : b), hi2 = (uint32_t)(a < b ? b : a);
return ((uint64_t)hi2 << 32) | lo2;
};
for (size_t i = 0; i + 2 < out.indices.size(); i += 3) {
const int a = out.indices[i], b = out.indices[i+1], c = out.indices[i+2];
if (a == b || b == c || a == c) continue; // degenerate: no edges
ecount[key(a,b)]++; ecount[key(b,c)]++; ecount[key(c,a)]++;
}
for (const auto& kv : ecount) {
if (kv.second == 1) out.boundary_edges++;
else if (kv.second > 2) out.nonmanifold_edges++;
}
}
out.solid_volume = (double)out.cells_solid * cell * cell * cell;
out.interior_volume = (double)out.cells_interior * cell * cell * cell;
// ---- diagnostics: area and the divergence volume of what we emitted ----
for (size_t i = 0; i + 2 < out.indices.size(); i += 3) {
const int ia = out.indices[i]*3, ib = out.indices[i+1]*3, ic = out.indices[i+2]*3;
const double ax = out.vertices[ia], ay = out.vertices[ia+1], az = out.vertices[ia+2];
double ux = out.vertices[ib] - ax, uy = out.vertices[ib+1] - ay, uz = out.vertices[ib+2] - az;
double vx2 = out.vertices[ic] - ax, vy2 = out.vertices[ic+1] - ay, vz2 = out.vertices[ic+2] - az;
double cxn = uy*vz2 - uz*vy2, cyn = uz*vx2 - ux*vz2, czn = ux*vy2 - uy*vx2;
double a = 0.5 * std::sqrt(cxn*cxn + cyn*cyn + czn*czn);
if (a <= 0.0) continue;
out.area += a;
double gx = (ax + out.vertices[ib] + out.vertices[ic]) / 3.0;
double gy = (ay + out.vertices[ib+1] + out.vertices[ic+1]) / 3.0;
double gz = (az + out.vertices[ib+2] + out.vertices[ic+2]) / 3.0;
out.volume += (gx*cxn + gy*cyn + gz*czn) / 6.0;
}
// Negative because the winding is reversed (see the header) — report the
// magnitude so callers compare it against a volume rather than a signed one.
out.volume = std::fabs(out.volume);
return true;
}
// ---------------------------------------------------------------------------
// Islands
// ---------------------------------------------------------------------------
namespace {
// Do two AABBs come within `pad` of each other on every axis?
inline bool near(const CSG4Brush& a, const CSG4Brush& b, float pad) {
return a.lo.x - pad <= b.hi.x && b.lo.x - pad <= a.hi.x &&
a.lo.y - pad <= b.hi.y && b.lo.y - pad <= a.hi.y &&
a.lo.z - pad <= b.hi.z && b.lo.z - pad <= a.hi.z;
}
int find_root(std::vector<int>& parent, int i) {
while (parent[i] != i) { parent[i] = parent[parent[i]]; i = parent[i]; }
return i;
}
} // namespace
bool build_islands(const std::vector<CSG4Brush>& brushes, ShellMesh& out,
float cell, float clamp_cells, float merge,
int* island_count) {
out = ShellMesh{};
if (brushes.empty() || cell <= 0.0f) return false;
// ---- 1. connected components over the brush AABBs ----
// O(n^2) and deliberately so: START.MAP is 1183 brushes, which is 700 K
// cheap comparisons — far below the cost of a single grid pass, and a
// spatial index here would be an optimization of the part that is not slow.
const int n = (int)brushes.size();
std::vector<int> parent(n);
for (int i = 0; i < n; i++) parent[i] = i;
for (int i = 0; i < n; i++) {
if (!brushes[i].bounds_valid) continue;
for (int j = i + 1; j < n; j++) {
if (!brushes[j].bounds_valid) continue;
if (!near(brushes[i], brushes[j], merge)) continue;
const int ri = find_root(parent, i), rj = find_root(parent, j);
if (ri != rj) parent[ri] = rj;
}
}
std::unordered_map<int, std::vector<CSG4Brush>> groups;
for (int i = 0; i < n; i++) {
if (!brushes[i].bounds_valid) continue;
groups[find_root(parent, i)].push_back(brushes[i]);
}
if (groups.empty()) return false;
if (island_count) *island_count = (int)groups.size();
// ---- 2. one grid per island, over its own extent ----
bool any = false;
for (auto& kv : groups) {
ShellMesh part;
if (!build(kv.second, part, cell, clamp_cells)) continue;
// CONCATENATED WITH AN INDEX OFFSET. A mesh SDF is happy with disjoint
// closed components — that is already what a detached clearance volume
// produces — so the pieces need no stitching, only renumbering.
const int base = (int)(out.vertices.size() / 3);
out.vertices.insert(out.vertices.end(), part.vertices.begin(), part.vertices.end());
for (int idx : part.indices) out.indices.push_back(idx + base);
// Diagnostics SUM, so `sealed()` still means what it meant: closed if
// every island is closed, and more air than wall over the whole level.
// The grid dims become the LARGEST island's, because there is no longer
// one grid to report and the biggest is the one that cost the most.
out.cells_solid += part.cells_solid;
out.cells_exterior += part.cells_exterior;
out.cells_interior += part.cells_interior;
out.snapped += part.snapped;
out.unmoved += part.unmoved;
out.volume += part.volume;
out.area += part.area;
out.boundary_edges += part.boundary_edges;
out.nonmanifold_edges += part.nonmanifold_edges;
out.solid_volume += part.solid_volume;
out.interior_volume += part.interior_volume;
if ((long long)part.grid[0] * part.grid[1] * part.grid[2] >
(long long)out.grid[0] * out.grid[1] * out.grid[2]) {
out.grid[0] = part.grid[0];
out.grid[1] = part.grid[1];
out.grid[2] = part.grid[2];
}
any = true;
}
return any;
}
} // namespace dc_level_shell