// ============================================================================
// START.MAP CREASE ORACLE — the multi-brush gate.
//
// tests/test_crease_segment_oracle.cpp measures a SINGLE isolated convex brush.
// That cannot see the things START.MAP actually does: 1164 overlapping brushes,
// cross-brush hidden-face removal (bake_planeset.glsl:106-129), capture pulling
// planes from DIFFERENT brushes, and — for the segment-clamp fix — the
// owner-brush lookup matching the WRONG brush when faces are coplanar.
//
// Input: /tmp/startmap_brushes.bin, written by
// hyper-tps/scenes/test/oracle_dump.gd (godot4 --headless --script ...)
// which records the REAL QuakeMapImporter output, so the brush set is identical
// to what the game stamps (same recenter, world_scale, per-face plane transform,
// material + projection ids).
//
// GROUND TRUTH — creases of the UNION boundary, computed exactly:
// For each brush, take its intra-brush crease segments (plane pair clipped by
// that brush's own half-spaces). A point on such a segment is on the union
// boundary only where it is NOT strictly inside any OTHER brush. "Inside brush
// Q" along a segment is the intersection of Q's half-space constraints = one
// interval in the segment parameter, so the buried part is subtracted exactly
// (no sampling).
//
// SCOPE LIMIT (stated, not hidden): this measures creases that are an EDGE OF A
// BRUSH and survive on the union boundary. It does NOT measure cross-brush
// creases — where a face of brush A meets a face of brush B to form a seam that
// is an edge of neither brush alone. Those exist in START.MAP and are a second
// increment.
//
// MODES: A = shipped pipeline. C = the segment clamp now in qef_solve.glsl,
// including the owner-brush lookup, ported line-for-line so a wrong-owner match
// shows up here rather than in the game.
//
// Build:
// g++ -O2 -std=c++17 -o /tmp/startmap_crease tests/test_startmap_crease.cpp
// /tmp/startmap_crease [voxel_size]
// ============================================================================
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <cmath>
#include <vector>
#include <algorithm>
#include <unordered_map>
// ---------------------------------------------------------------------------
struct V3 { float x, y, z; };
static inline V3 operator+(V3 a, V3 b) { return {a.x+b.x, a.y+b.y, a.z+b.z}; }
static inline V3 operator-(V3 a, V3 b) { return {a.x-b.x, a.y-b.y, a.z-b.z}; }
static inline V3 operator*(V3 a, float s) { return {a.x*s, a.y*s, a.z*s}; }
static inline float dot(V3 a, V3 b) { return a.x*b.x + a.y*b.y + a.z*b.z; }
static inline V3 cross3(V3 a, V3 b) {
return {a.y*b.z-a.z*b.y, a.z*b.x-a.x*b.z, a.x*b.y-a.y*b.x};
}
static inline float len3(V3 a) { return std::sqrt(dot(a,a)); }
static inline V3 norm3(V3 a) { float l = len3(a); return l > 1e-20f ? a*(1.0f/l) : a; }
static inline float clampf(float v, float lo, float hi) { return v<lo?lo:(v>hi?hi:v); }
static inline float compAt(V3 v, int i) { return i==0?v.x:(i==1?v.y:v.z); }
struct Plane { V3 n; float d; unsigned proj; };
struct Brush {
V3 center; float radius; float strength; unsigned material;
std::vector<Plane> planes;
V3 lo, hi; // AABB, for the spatial grid
};
// ---------------------------------------------------------------------------
// Dump loader
// ---------------------------------------------------------------------------
static bool load_dump(const char* path, std::vector<Brush>& out, float& world_scale) {
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 in %s\n", path); fclose(f); return false;
}
unsigned ver = 0; size_t rd = fread(&ver,4,1,f); (void)rd;
if (ver != 1) { fprintf(stderr,"unsupported version %u\n",ver); fclose(f); return false; }
rd = fread(&world_scale,4,1,f);
float carve[6]; rd = fread(carve,4,6,f); // carve half + center (unused here)
unsigned n = 0; rd = fread(&n,4,1,f);
out.clear(); out.reserve(n);
for (unsigned i = 0; i < n; i++) {
Brush b; float c[3];
rd = fread(c,4,3,f); b.center = {c[0],c[1],c[2]};
rd = fread(&b.radius,4,1,f);
rd = fread(&b.strength,4,1,f);
rd = fread(&b.material,4,1,f);
unsigned pc = 0; rd = fread(&pc,4,1,f);
b.planes.reserve(pc);
for (unsigned p = 0; p < pc; p++) {
float v[4]; unsigned pr;
rd = fread(v,4,4,f); rd = fread(&pr,4,1,f);
b.planes.push_back({ {v[0],v[1],v[2]}, v[3], pr });
}
// AABB from the bounding sphere (planes alone need an LP to bound tightly)
b.lo = b.center - V3{b.radius,b.radius,b.radius};
b.hi = b.center + V3{b.radius,b.radius,b.radius};
out.push_back(std::move(b));
}
fclose(f);
return true;
}
// convex SDF of one brush: max_i(n_i.p - d_i), negative inside
static inline float brush_f(V3 p, const Brush& b) {
float m = -1e30f;
for (const auto& q : b.planes) m = std::max(m, dot(q.n,p) - q.d);
return m;
}
// ---------------------------------------------------------------------------
// Spatial grid over brushes (the whole map at 0.224 voxels is far too big to
// scan; everything below is driven off the crease segments instead).
// ---------------------------------------------------------------------------
struct BrushGrid {
float cell = 4.0f;
V3 origin{0,0,0};
std::unordered_map<long long, std::vector<int>> buckets;
static long long key(int x,int y,int z) {
return ((long long)(x+4096)<<42) ^ ((long long)(y+4096)<<21) ^ (long long)(z+4096);
}
void build(const std::vector<Brush>& bs, float expand) {
for (int i = 0; i < (int)bs.size(); i++) {
V3 lo = bs[i].lo - V3{expand,expand,expand};
V3 hi = bs[i].hi + V3{expand,expand,expand};
int x0=(int)std::floor(lo.x/cell), x1=(int)std::floor(hi.x/cell);
int y0=(int)std::floor(lo.y/cell), y1=(int)std::floor(hi.y/cell);
int z0=(int)std::floor(lo.z/cell), z1=(int)std::floor(hi.z/cell);
for (int z=z0; z<=z1; z++)
for (int y=y0; y<=y1; y++)
for (int x=x0; x<=x1; x++)
buckets[key(x,y,z)].push_back(i);
}
}
const std::vector<int>* at(V3 p) const {
auto it = buckets.find(key((int)std::floor(p.x/cell),
(int)std::floor(p.y/cell),
(int)std::floor(p.z/cell)));
return it == buckets.end() ? nullptr : &it->second;
}
};
static BrushGrid g_grid;
static const std::vector<Brush>* g_brushes = nullptr;
// Union density as apply_gpu_edits builds it for ADD brushes:
// d = max_i(-f_i) = -min_i(f_i) [noise_sdf.glsl:177-193]
// Terrain noise and the carve box are omitted: inside the map volume the carve
// makes the base field air and the brushes add the solid back, so the union of
// the brushes IS the surface there. Stated assumption.
static float density_at(V3 p) {
const std::vector<int>* v = g_grid.at(p);
if (!v) return -1e30f;
float mn = 1e30f;
for (int i : *v) {
const Brush& b = (*g_brushes)[i];
V3 dc = p - b.center;
if (dot(dc,dc) > b.radius*b.radius) continue;
mn = std::min(mn, brush_f(p, b));
}
return (mn > 1e29f) ? -1e30f : -mn;
}
// ---------------------------------------------------------------------------
// Interval arithmetic for exact union-boundary clipping
// ---------------------------------------------------------------------------
struct Iv { float a, b; };
// t-interval where p0 + t*dir is strictly inside brush q (margin shrinks it)
static bool inside_interval(V3 p0, V3 dir, const Brush& q, float margin, Iv& out) {
float lo = -1e30f, hi = 1e30f;
for (const auto& pl : q.planes) {
float a = dot(pl.n, dir);
float b = (pl.d - margin) - dot(pl.n, p0);
if (std::fabs(a) < 1e-9f) { if (b < 0.0f) return false; }
else if (a > 0.0f) hi = std::min(hi, b/a);
else lo = std::max(lo, b/a);
}
if (lo >= hi) return false;
out = {lo, hi};
return true;
}
static void subtract(std::vector<Iv>& keep, Iv cut) {
std::vector<Iv> out;
out.reserve(keep.size() + 1);
for (const Iv& s : keep) {
if (cut.b <= s.a || cut.a >= s.b) { out.push_back(s); continue; }
if (cut.a > s.a) out.push_back({s.a, cut.a});
if (cut.b < s.b) out.push_back({cut.b, s.b});
}
keep.swap(out);
}
struct Segment { V3 a, b; V3 n0; float d0; V3 n1; float d1; int brush; };
// Intra-brush crease segments that survive on the UNION boundary.
static std::vector<Segment> union_crease_segments(const std::vector<Brush>& bs, float min_len) {
std::vector<Segment> segs;
for (int bi = 0; bi < (int)bs.size(); bi++) {
const Brush& B = bs[bi];
int n = (int)B.planes.size();
for (int i = 0; i < n; i++)
for (int j = i+1; j < n; j++) {
V3 ni = B.planes[i].n, nj = B.planes[j].n;
V3 dir = cross3(ni, nj);
if (len3(dir) < 1e-6f) continue;
dir = norm3(dir);
float cc = dot(ni,nj), det = 1.0f - cc*cc;
if (det < 1e-8f) continue;
V3 p0 = ni*((B.planes[i].d - cc*B.planes[j].d)/det)
+ nj*((B.planes[j].d - cc*B.planes[i].d)/det);
// clip to this brush's own half-spaces => the brush's real edge
float tmin = -1e30f, tmax = 1e30f; bool empty = false;
for (int k = 0; k < n && !empty; k++) {
if (k==i || k==j) continue;
float a = dot(B.planes[k].n, dir);
float b = B.planes[k].d - dot(B.planes[k].n, p0);
if (std::fabs(a) < 1e-9f) { if (b < -1e-6f) empty = true; }
else if (a > 0.0f) tmax = std::min(tmax, b/a);
else tmin = std::max(tmin, b/a);
}
if (empty || tmax - tmin <= min_len) continue;
// subtract the parts buried strictly inside any OTHER brush
std::vector<Iv> keep{{tmin, tmax}};
V3 mid = p0 + dir*(0.5f*(tmin+tmax));
const std::vector<int>* cand = g_grid.at(mid);
// walk buckets along the edge so long edges are not under-tested
std::vector<int> pool;
{
int steps = (int)std::ceil((tmax-tmin)/g_grid.cell) + 1;
for (int s = 0; s <= steps; s++) {
V3 p = p0 + dir*(tmin + (tmax-tmin)*(float)s/(float)steps);
const std::vector<int>* v = g_grid.at(p);
if (v) pool.insert(pool.end(), v->begin(), v->end());
}
std::sort(pool.begin(), pool.end());
pool.erase(std::unique(pool.begin(), pool.end()), pool.end());
}
(void)cand;
for (int oi : pool) {
if (oi == bi || keep.empty()) continue;
Iv iv;
if (inside_interval(p0, dir, bs[oi], 1e-4f, iv)) subtract(keep, iv);
}
for (const Iv& s : keep) {
if (s.b - s.a <= min_len) continue;
segs.push_back({ p0 + dir*s.a, p0 + dir*s.b,
ni, B.planes[i].d, nj, B.planes[j].d, bi });
}
}
}
return segs;
}
static float point_line_dist(V3 p, V3 a, V3 b) {
V3 dir = norm3(b - a), w = p - a;
return len3(w - dir*dot(w,dir));
}
static bool seg_hits_box(V3 A, V3 B, V3 lo, V3 hi) {
V3 d = B - A; float t0 = 0.0f, t1 = 1.0f;
for (int a = 0; a < 3; a++) {
float da = compAt(d,a), pa = compAt(A,a), l = compAt(lo,a), h = compAt(hi,a);
if (std::fabs(da) < 1e-12f) { if (pa < l || pa > h) return false; }
else {
float inv = 1.0f/da, ta = (l-pa)*inv, tb = (h-pa)*inv;
if (ta > tb) std::swap(ta, tb);
t0 = std::max(t0, ta); t1 = std::min(t1, tb);
if (t0 > t1) return false;
}
}
return true;
}
// ---------------------------------------------------------------------------
// Pipeline replay (multi-brush)
// ---------------------------------------------------------------------------
struct Hermite { V3 pos, n; bool hard; };
static Hermite edge_hermite(V3 p0, V3 p1, float d0, float vs) {
float lo = 0.0f, hi = 1.0f, d_lo = d0;
for (int i = 0; i < 8; i++) { // [edge_intersect.glsl:196-208]
float mid = (lo+hi)*0.5f;
float dm = density_at(p0 + (p1-p0)*mid);
if ((d_lo >= 0.0f) == (dm >= 0.0f)) { lo = mid; d_lo = dm; } else hi = mid;
}
float t = clampf((lo+hi)*0.5f, 0.001f, 0.999f);
V3 pos = p0 + (p1-p0)*t;
float eps = vs*0.01f;
V3 n = { density_at({pos.x+eps,pos.y,pos.z}) - density_at({pos.x-eps,pos.y,pos.z}),
density_at({pos.x,pos.y+eps,pos.z}) - density_at({pos.x,pos.y-eps,pos.z}),
density_at({pos.x,pos.y,pos.z+eps}) - density_at({pos.x,pos.y,pos.z-eps}) };
float L = len3(n);
n = (L > 1e-8f) ? n*(1.0f/L) : V3{0,1,0};
V3 cdn = n;
// exact analytic override, argmax face, nearest brush wins [:249-288]
bool hit = false; Plane hit_pl{{0,0,0},0,0};
float min_absf = vs*1.8f;
const std::vector<int>* v = g_grid.at(pos);
if (v) for (int bi : *v) {
const Brush& b = (*g_brushes)[bi];
V3 dc = pos - b.center;
if (dot(dc,dc) > b.radius*b.radius) continue;
float best = -1e30f; Plane bp{{0,0,0},0,0};
for (const auto& q : b.planes) { float dp = dot(q.n,pos)-q.d; if (dp > best) { best = dp; bp = q; } }
if (std::fabs(best) < min_absf) { min_absf = std::fabs(best); hit_pl = bp; hit = true; }
}
if (hit) {
V3 ed = p1 - p0;
float den = dot(hit_pl.n, ed);
if (std::fabs(den) > 1e-8f)
pos = p0 + ed*clampf((hit_pl.d - dot(hit_pl.n,p0))/den, 0.0f, 1.0f);
n = (dot(hit_pl.n, cdn) < 0.0f) ? hit_pl.n*-1.0f : hit_pl.n;
}
return { pos, n, hit };
}
struct Capture { V3 n[3]; float d[3]; int count; };
// bake_planeset.glsl main(), INCLUDING the cross-brush buried test [:106-129]
static Capture bake_capture(V3 cc, float vs) {
Capture c; c.count = 0;
float T = vs*1.2f;
const std::vector<int>* v = g_grid.at(cc);
if (!v) return c;
for (int ei : *v) {
const Brush& B = (*g_brushes)[ei];
V3 dcv = cc - B.center; float rr = B.radius + T;
if (dot(dcv,dcv) > rr*rr) continue;
for (int pi = 0; pi < (int)B.planes.size(); pi++) {
const Plane& pl = B.planes[pi];
float fd = dot(pl.n,cc) - pl.d;
if (std::fabs(fd) >= T) continue;
bool onbrush = true;
for (int pj = 0; pj < (int)B.planes.size(); pj++) {
if (pj == pi) continue;
if (dot(B.planes[pj].n,cc) - B.planes[pj].d > T) { onbrush = false; break; }
}
if (!onbrush) continue;
// csg4 cross-brush hidden-face removal
V3 fpt = cc - pl.n*fd;
V3 qp = fpt + pl.n*(vs*0.5f);
bool buried = false;
const std::vector<int>* v2 = g_grid.at(qp);
if (v2) for (int ej : *v2) {
if (ej == ei || buried) continue;
const Brush& B2 = (*g_brushes)[ej];
V3 dc2 = qp - B2.center;
if (dot(dc2,dc2) > B2.radius*B2.radius) continue;
float mx = -1e30f;
for (const auto& p2 : B2.planes) mx = std::max(mx, dot(p2.n,qp) - p2.d);
if (mx < -vs) buried = true; // [:127]
}
if (buried) continue;
int matched = -1;
for (int k = 0; k < c.count; k++)
if (dot(pl.n, c.n[k]) > 0.999f) { matched = k; break; }
if (matched < 0 && c.count < 3) {
c.n[c.count] = pl.n; c.d[c.count] = pl.d; c.count++;
}
}
}
return c;
}
enum Mode { MODE_A = 0, MODE_C = 1, MODE_F = 2, MODE_E = 3 };
struct Res { V3 vert; bool escaped=false; bool snap_ok=false; bool wrong_owner=false; };
static Res solve_cell(const std::vector<Hermite>& H, V3 cmin, float vs,
const Capture& cap, Mode mode, int true_brush) {
Res R;
V3 cmax = cmin + V3{vs,vs,vs};
float ata[6]={0,0,0,0,0,0}; V3 atb{0,0,0}, mass{0,0,0}; int cnt=0;
for (const auto& h : H) {
V3 nm = h.n;
ata[0]+=nm.x*nm.x; ata[1]+=nm.x*nm.y; ata[2]+=nm.x*nm.z;
ata[3]+=nm.y*nm.y; ata[4]+=nm.y*nm.z; ata[5]+=nm.z*nm.z;
atb = atb + nm*dot(nm,h.pos); mass = mass + h.pos; cnt++;
}
if (!cnt) { R.vert = cmin + V3{vs*.5f,vs*.5f,vs*.5f}; return R; }
mass = mass*(1.0f/cnt);
float eps = 0.01f;
float m00=ata[0]+eps,m01=ata[1],m02=ata[2],m11=ata[3]+eps,m12=ata[4],m22=ata[5]+eps;
V3 amp{ ata[0]*mass.x+ata[1]*mass.y+ata[2]*mass.z,
ata[1]*mass.x+ata[3]*mass.y+ata[4]*mass.z,
ata[2]*mass.x+ata[4]*mass.y+ata[5]*mass.z };
V3 rhs = atb - amp;
float c00=m11*m22-m12*m12,c01=m02*m12-m01*m22,c02=m01*m12-m02*m11;
float c11=m00*m22-m02*m02,c12=m01*m02-m00*m12,c22=m00*m11-m01*m01;
float det = m00*c00+m01*c01+m02*c02;
V3 result = mass;
if (std::fabs(det) > 1e-12f) {
float iv = 1.0f/det;
result = mass + V3{ (c00*rhs.x+c01*rhs.y+c02*rhs.z)*iv,
(c01*rhs.x+c11*rhs.y+c12*rhs.z)*iv,
(c02*rhs.x+c12*rhs.y+c22*rhs.z)*iv };
}
result = { clampf(result.x,cmin.x,cmax.x), clampf(result.y,cmin.y,cmax.y),
clampf(result.z,cmin.z,cmax.z) };
if (cap.count >= 1) {
V3 snapped = result;
bool snap_ok = false;
if (cap.count == 1) {
snapped = result - cap.n[0]*(dot(cap.n[0],result) - cap.d[0]);
} else if (mode == MODE_C || mode == MODE_F || mode == MODE_E) {
int ia = 0, ib = 1;
if (cap.count == 3) {
V3 a=cap.n[0],b=cap.n[1],cn=cap.n[2];
float den = dot(a, cross3(b,cn));
if (std::fabs(den) > 1e-5f) {
V3 corner = (cross3(b,cn)*cap.d[0]+cross3(cn,a)*cap.d[1]
+ cross3(a,b)*cap.d[2])*(1.0f/den);
if (corner.x>=cmin.x&&corner.x<=cmax.x&&corner.y>=cmin.y&&
corner.y<=cmax.y&&corner.z>=cmin.z&&corner.z<=cmax.z) {
snapped = corner; snap_ok = true;
}
}
if (!snap_ok) {
V3 ctr = cmin + V3{vs*.5f,vs*.5f,vs*.5f};
float best = 1e30f;
for (int k = 0; k < 3; k++) {
int p=k, q=(k+1)%3;
float cc2 = dot(cap.n[p],cap.n[q]), dd2 = 1.0f-cc2*cc2;
if (dd2 <= 1e-4f) continue;
float r0 = cap.d[p]-dot(cap.n[p],ctr), r1 = cap.d[q]-dot(cap.n[q],ctr);
V3 pt = ctr + cap.n[p]*((r0-cc2*r1)/dd2) + cap.n[q]*((r1-cc2*r0)/dd2);
float dist = len3(pt-ctr);
if (dist < best) { best = dist; ia=p; ib=q; }
}
}
}
if (!snap_ok) {
float cc2 = dot(cap.n[ia],cap.n[ib]), dd2 = 1.0f-cc2*cc2;
if (dd2 > 1e-4f) {
float r0 = cap.d[ia]-dot(cap.n[ia],result), r1 = cap.d[ib]-dot(cap.n[ib],result);
V3 pline = result + cap.n[ia]*((r0-cc2*r1)/dd2) + cap.n[ib]*((r1-cc2*r0)/dd2);
V3 ldir = cross3(cap.n[ia],cap.n[ib]);
if (len3(ldir) > 1e-9f) {
ldir = norm3(ldir);
float tmin=-1e30f, tmax=1e30f; bool ok = true;
for (int a2 = 0; a2 < 3; a2++) {
float pv = compAt(pline,a2), dv = compAt(ldir,a2);
float l = compAt(cmin,a2), h = compAt(cmax,a2);
if (std::fabs(dv) < 1e-9f) { if (pv < l || pv > h) ok = false; }
else {
float ta=(l-pv)/dv, tb=(h-pv)/dv;
tmin = std::max(tmin, std::min(ta,tb));
tmax = std::min(tmax, std::max(ta,tb));
}
}
// owner-brush lookup — the risky part on a real map
if (ok) {
const std::vector<int>* v = g_grid.at(pline);
if (v) for (int ei : *v) {
const Brush& B = (*g_brushes)[ei];
V3 dcv = pline - B.center;
if (dot(dcv,dcv) > B.radius*B.radius) continue;
// Match the PLANE (normal AND offset), not just the
// normal: START.MAP brushes are overwhelmingly axis
// aligned, so a normal-only test matches hundreds of
// unrelated brushes and the first one wins.
bool has_a=false, has_b=false;
for (const auto& q : B.planes) {
if (dot(q.n,cap.n[ia]) > 0.999f && std::fabs(q.d-cap.d[ia]) < 1e-3f) has_a = true;
if (dot(q.n,cap.n[ib]) > 0.999f && std::fabs(q.d-cap.d[ib]) < 1e-3f) has_b = true;
}
if (!has_a || !has_b) continue;
if (true_brush >= 0 && ei != true_brush) R.wrong_owner = true;
for (const auto& q : B.planes) {
if (dot(q.n,cap.n[ia]) > 0.999f && std::fabs(q.d-cap.d[ia]) < 1e-3f) continue;
if (dot(q.n,cap.n[ib]) > 0.999f && std::fabs(q.d-cap.d[ib]) < 1e-3f) continue;
float aa = dot(q.n,ldir), bb = q.d - dot(q.n,pline);
if (std::fabs(aa) < 1e-9f) { if (bb < -1e-6f) ok = false; }
else if (aa > 0.0f) tmax = std::min(tmax, bb/aa);
else tmin = std::max(tmin, bb/aa);
}
break;
}
}
if (ok && tmin <= tmax) {
snapped = pline + ldir*clampf(0.0f,tmin,tmax);
snap_ok = true;
} else if (mode == MODE_F) {
// Empty range PROVES no crease crosses this cell: the
// T=1.2*voxel centre test tagged a FALSE crease. Treat it
// as the FACE cell it actually is — a point on a face plane
// inside the cell is ON the surface and cannot protrude.
V3 ctr = cmin + V3{vs*.5f,vs*.5f,vs*.5f};
int nb=0; float ba=1e30f;
for (int c = 0; c < cap.count; c++) {
float aa = std::fabs(dot(cap.n[c],ctr)-cap.d[c]);
if (aa < ba) { ba = aa; nb = c; }
}
snapped = result - cap.n[nb]*(dot(cap.n[nb],result)-cap.d[nb]);
snap_ok = true;
}
}
}
}
} else { // MODE_A crease / corner snap
if (cap.count == 2) {
float cc2 = dot(cap.n[0],cap.n[1]), dd2 = 1.0f-cc2*cc2;
if (dd2 > 1e-4f) {
float r0 = cap.d[0]-dot(cap.n[0],result), r1 = cap.d[1]-dot(cap.n[1],result);
snapped = result + cap.n[0]*((r0-cc2*r1)/dd2) + cap.n[1]*((r1-cc2*r0)/dd2);
}
} else {
V3 a=cap.n[0],b=cap.n[1],cn=cap.n[2];
float den = dot(a,cross3(b,cn));
if (std::fabs(den) > 1e-5f)
snapped = (cross3(b,cn)*cap.d[0]+cross3(cn,a)*cap.d[1]+cross3(a,b)*cap.d[2])*(1.0f/den);
}
}
R.snap_ok = snap_ok;
if (cap.count >= 2 && !snap_ok) { // escape fallback
V3 lo = cmin - V3{vs*.5f,vs*.5f,vs*.5f}, hi = cmax + V3{vs*.5f,vs*.5f,vs*.5f};
if (snapped.x<lo.x||snapped.y<lo.y||snapped.z<lo.z||
snapped.x>hi.x||snapped.y>hi.y||snapped.z>hi.z) {
V3 ctr = cmin + V3{vs*.5f,vs*.5f,vs*.5f};
int nb=0; float ba=1e30f;
for (int c = 0; c < cap.count; c++) {
float aa = std::fabs(dot(cap.n[c],ctr)-cap.d[c]);
if (aa < ba) { ba = aa; nb = c; }
}
snapped = result - cap.n[nb]*(dot(cap.n[nb],result)-cap.d[nb]);
}
}
result = { clampf(snapped.x,cmin.x,cmax.x), clampf(snapped.y,cmin.y,cmax.y),
clampf(snapped.z,cmin.z,cmax.z) };
}
R.vert = result;
// ---- MODE E: CONTAINMENT. Nothing in the snap path keeps the vertex INSIDE
// the union. If it ends up outside, pull it back onto the nearest brush's
// active face, then re-clamp to the cell. (Surface Nets is immune to this by
// construction, which is what suggested the constraint.)
if (mode == MODE_E) {
if (density_at(R.vert) < 0.0f) { // negative = OUTSIDE the union
const std::vector<int>* v = g_grid.at(R.vert);
float bestf = 1e30f; V3 bn{0,0,0}; float bd = 0.0f; bool found = false;
if (v) for (int bi : *v) {
const Brush& B = (*g_brushes)[bi];
float f = brush_f(R.vert, B);
if (f < bestf) { // brush whose surface is nearest
bestf = f; found = true;
float best = -1e30f;
for (const auto& q : B.planes) {
float dp = dot(q.n, R.vert) - q.d;
if (dp > best) { best = dp; bn = q.n; bd = q.d; }
}
}
}
if (found) {
V3 proj = R.vert - bn * (dot(bn, R.vert) - bd);
R.vert = { clampf(proj.x,cmin.x,cmax.x), clampf(proj.y,cmin.y,cmax.y),
clampf(proj.z,cmin.z,cmax.z) };
}
}
}
const float E = 1e-4f;
if (R.vert.x<cmin.x-E||R.vert.x>cmax.x+E||R.vert.y<cmin.y-E||
R.vert.y>cmax.y+E||R.vert.z<cmin.z-E||R.vert.z>cmax.z+E) R.escaped = true;
return R;
}
struct Stats { float p50=0,p95=0,mx=0; int n=0; };
static Stats stat_of(std::vector<float> v) {
Stats s; s.n=(int)v.size();
if (v.empty()) return s;
std::sort(v.begin(), v.end());
s.p50 = v[v.size()/2];
size_t i = (size_t)(v.size()*0.95f); if (i >= v.size()) i = v.size()-1;
s.p95 = v[i]; s.mx = v.back();
return s;
}
int main(int argc, char** argv) {
float vs = 0.224f; // map_dc_terrain.gd:86
if (argc > 1) vs = (float)atof(argv[1]);
const char* dump = (argc > 2) ? argv[2] : "/tmp/startmap_brushes.bin";
std::vector<Brush> brushes; float world_scale = 0;
if (!load_dump(dump, brushes, world_scale)) return 1;
g_brushes = &brushes;
g_grid.build(brushes, vs*4.0f);
printf("START.MAP CREASE ORACLE — multi-brush gate\n");
printf("brushes=%zu world_scale=%.4f voxel=%.3f\n", brushes.size(), world_scale, vs);
auto segs = union_crease_segments(brushes, vs*0.5f);
double total_len = 0;
for (auto& s : segs) total_len += len3(s.b - s.a);
printf("union-boundary crease segments=%zu total length=%.1f m\n", segs.size(), total_len);
static const int EOFF[12][4] = {
{0,0,0,0},{0,1,0,0},{0,0,1,0},{0,1,1,0},
{0,0,0,1},{1,0,0,1},{0,0,1,1},{1,0,1,1},
{0,0,0,2},{1,0,0,2},{0,1,0,2},{1,1,0,2} };
static const int EDIR[3][3] = {{1,0,0},{0,1,0},{0,0,1}};
// Walk cells along each segment (plus a 1-cell shell for false-crease stats).
std::unordered_map<long long,int> visited; // cell -> segment index
auto ckey = [](int x,int y,int z){
return ((long long)(x+1000000)<<42) ^ ((long long)(y+1000000)<<21) ^ (long long)(z+1000000);
};
std::vector<float> devA, devC, devF;
std::vector<float> protA, protC, protF;
std::vector<float> protE, devE;
int escF=0, escE=0, pcA=0, pcC=0, pcF=0, pcE=0;
int crease_cells=0, escA=0, escC=0, wrong_owner=0, segfail=0;
int cl_hist[4]={0,0,0,0};
int false_cells=0; std::vector<float> falseA, falseC;
for (size_t si = 0; si < segs.size(); si++) {
const Segment& S = segs[si];
float L = len3(S.b - S.a);
int steps = (int)std::ceil(L/(vs*0.5f)) + 1;
for (int st = 0; st <= steps; st++) {
V3 p = S.a + (S.b - S.a)*((float)st/(float)steps);
int bx=(int)std::floor(p.x/vs), by=(int)std::floor(p.y/vs), bz=(int)std::floor(p.z/vs);
for (int dz=-1; dz<=1; dz++)
for (int dy=-1; dy<=1; dy++)
for (int dx=-1; dx<=1; dx++) {
int cx=bx+dx, cy=by+dy, cz=bz+dz;
long long k = ckey(cx,cy,cz);
if (visited.count(k)) continue;
visited[k] = (int)si;
V3 cmin{cx*vs, cy*vs, cz*vs};
V3 cmax = cmin + V3{vs,vs,vs};
V3 cctr = cmin + V3{vs*.5f,vs*.5f,vs*.5f};
std::vector<Hermite> H; H.reserve(12);
for (int e = 0; e < 12; e++) {
int ex=cx+EOFF[e][0], ey=cy+EOFF[e][1], ez=cz+EOFF[e][2], dir=EOFF[e][3];
V3 q0{ex*vs, ey*vs, ez*vs};
V3 q1 = q0 + V3{EDIR[dir][0]*vs, EDIR[dir][1]*vs, EDIR[dir][2]*vs};
float d0 = density_at(q0), d1 = density_at(q1);
if ((d0 >= 0.0f) == (d1 >= 0.0f)) continue;
H.push_back(edge_hermite(q0, q1, d0, vs));
}
if (H.empty()) continue;
Capture cap = bake_capture(cctr, vs);
// is a TRUE union crease in this cell?
int hit_seg = -1;
for (size_t sj = 0; sj < segs.size() && hit_seg < 0; sj++) {
// cheap reject: only test the segment that led us here + neighbours
if (sj != si) continue;
if (seg_hits_box(segs[sj].a, segs[sj].b, cmin, cmax)) hit_seg = (int)sj;
}
Res A = solve_cell(H, cmin, vs, cap, MODE_A, -1);
Res C = solve_cell(H, cmin, vs, cap, MODE_C, segs[si].brush);
Res F = solve_cell(H, cmin, vs, cap, MODE_F, segs[si].brush);
Res Emo = solve_cell(H, cmin, vs, cap, MODE_E, segs[si].brush);
if (Emo.escaped) escE++;
if (A.escaped) escA++;
if (C.escaped) escC++;
if (F.escaped) escF++;
if (C.wrong_owner) wrong_owner++;
// outward protrusion: density_at < 0 means OUTSIDE the union solid
{ float dA=density_at(A.vert), dC=density_at(C.vert), dF=density_at(F.vert);
if (dA < 0 && dA > -1e29f) { protA.push_back(-dA); if (-dA > vs*0.1f) pcA++; }
if (dC < 0 && dC > -1e29f) { protC.push_back(-dC); if (-dC > vs*0.1f) pcC++; }
if (dF < 0 && dF > -1e29f) { protF.push_back(-dF); if (-dF > vs*0.1f) pcF++; }
float dE=density_at(Emo.vert);
if (dE < 0 && dE > -1e29f) { protE.push_back(-dE); if (-dE > vs*0.1f) pcE++; } }
if (hit_seg >= 0) {
crease_cells++;
cl_hist[std::min(cap.count,3)]++;
devA.push_back(point_line_dist(A.vert, segs[hit_seg].a, segs[hit_seg].b));
devC.push_back(point_line_dist(C.vert, segs[hit_seg].a, segs[hit_seg].b));
devF.push_back(point_line_dist(F.vert, segs[hit_seg].a, segs[hit_seg].b));
devE.push_back(point_line_dist(Emo.vert, segs[hit_seg].a, segs[hit_seg].b));
if (cap.count >= 2 && !C.snap_ok) segfail++;
} else if (cap.count == 2) {
false_cells++;
// distance from the union surface (|density| ~ SDF near the surface)
falseA.push_back(std::fabs(density_at(A.vert)));
falseC.push_back(std::fabs(density_at(C.vert)));
}
}
}
}
Stats sA = stat_of(devA), sC = stat_of(devC);
Stats fA = stat_of(falseA), fC = stat_of(falseC);
printf("\ncells evaluated=%zu crease cells=%d\n", visited.size(), crease_cells);
printf(" captured cl_count on crease cells: 0=%d 1(MISSED)=%d 2=%d 3=%d (missed %.1f%%)\n",
cl_hist[0], cl_hist[1], cl_hist[2], cl_hist[3],
crease_cells ? 100.0f*(cl_hist[0]+cl_hist[1])/crease_cells : 0.0f);
printf("\nDEVIATION FROM THE TRUE UNION CREASE LINE (n=%d)\n", sA.n);
printf(" %-14s %-10s %-10s %-12s %-10s %s\n","mode","p50","p95","MAX","MAX/voxel","escapes");
printf(" %-14s %-10.5f %-10.5f %-12.5f %-10.2f %d\n","A CURRENT",sA.p50,sA.p95,sA.mx,sA.mx/vs,escA);
printf(" %-14s %-10.5f %-10.5f %-12.5f %-10.2f %d\n","C SEGMENT",sC.p50,sC.p95,sC.mx,sC.mx/vs,escC);
{ Stats sF = stat_of(devF);
printf(" %-14s %-10.5f %-10.5f %-12.5f %-10.2f %d\n","F FACEFALL",sF.p50,sF.p95,sF.mx,sF.mx/vs,escF);
Stats sE = stat_of(devE);
printf(" %-14s %-10.5f %-10.5f %-12.5f %-10.2f %d\n","E CONTAIN",sE.p50,sE.p95,sE.mx,sE.mx/vs,escE); }
{ Stats pA=stat_of(protA), pC=stat_of(protC), pF=stat_of(protF);
printf("\nOUTWARD PROTRUSION (vertex outside the union solid = visible notch)\n");
printf(" %-14s p50=%-9.5f p95=%-9.5f MAX=%-9.5f cells>0.1vox=%d\n","A CURRENT",pA.p50,pA.p95,pA.mx,pcA);
printf(" %-14s p50=%-9.5f p95=%-9.5f MAX=%-9.5f cells>0.1vox=%d\n","C SEGMENT",pC.p50,pC.p95,pC.mx,pcC);
printf(" %-14s p50=%-9.5f p95=%-9.5f MAX=%-9.5f cells>0.1vox=%d\n","F FACEFALL",pF.p50,pF.p95,pF.mx,pcF);
Stats pE=stat_of(protE);
printf(" %-14s p50=%-9.5f p95=%-9.5f MAX=%-9.5f cells>0.1vox=%d\n","E CONTAIN",pE.p50,pE.p95,pE.mx,pcE); }
printf("\nOWNER-BRUSH LOOKUP (the multi-brush risk)\n");
printf(" wrong-brush matches = %d segment-clamp could not resolve = %d\n",
wrong_owner, segfail);
printf("\nFALSE CREASE (cl==2, no true crease in cell): %d cells\n", false_cells);
printf(" |density| at the vertex (0 = on the union surface)\n");
printf(" %-14s p50=%-9.5f p95=%-9.5f MAX=%.5f\n","A CURRENT",fA.p50,fA.p95,fA.mx);
printf(" %-14s p50=%-9.5f p95=%-9.5f MAX=%.5f\n","C SEGMENT",fC.p50,fC.p95,fC.mx);
return 0;
}