// PHASE 0 repro: does the plane-set capture SHATTER on a DENSE overlapping-brush
// junction, meshed through the trilinear brick — and does cross-brush hidden-face
// removal fix it? The committed tests only cover ONE cube / ONE crease; this models
// what the red screenshots show: overlapping brushes at a junction, running the
// EXACT qef_solve Step-5b capture (T-radius, inside-same-brush-other-planes, dedup
// 0.999, 2-plane crease det>1e-4, 1-voxel snap-cap).
//
// Mechanism targeted: two brush faces ~5deg apart (just past the 2.56deg dedup
// threshold) both captured -> phantom near-parallel "crease" -> crease-line solve
// flings the vertex far -> capped (== the red cells). Hidden-face removal skips a
// face buried inside another brush -> only 1 plane -> no phantom -> no fly.
#include <cstdio>
#include <cmath>
#include <vector>
#include <algorithm>
struct V3{float x,y,z;};
static V3 operator-(V3 a,V3 b){return{a.x-b.x,a.y-b.y,a.z-b.z};}
static V3 operator+(V3 a,V3 b){return{a.x+b.x,a.y+b.y,a.z+b.z};}
static V3 operator*(V3 a,float s){return{a.x*s,a.y*s,a.z*s};}
static float dot(V3 a,V3 b){return a.x*b.x+a.y*b.y+a.z*b.z;}
static 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};}
struct Plane{V3 n;float d;}; // n.p <= d inside
struct Brush{std::vector<Plane> pl;};
static float bsdf(V3 p,const Brush&b){float m=-1e30f;for(auto&q:b.pl)m=fmaxf(m,dot(q.n,p)-q.d);return m;} // convex, neg inside
// union of solids = inside ANY brush = min of per-brush sdf
static float usdf(V3 p,const std::vector<Brush>&B){float m=1e30f;for(auto&b:B)m=fminf(m,bsdf(p,b));return m;}
static Brush mkbox(V3 lo,V3 hi){
Brush b; b.pl={{{1,0,0},hi.x},{{-1,0,0},-lo.x},{{0,1,0},hi.y},{{0,-1,0},-lo.y},{{0,0,1},hi.z},{{0,0,-1},-lo.z}}; return b;}
// a box whose TOP face is tilted by `ang` rad about Y (near-parallel-to-flat generator)
static Brush mktilt(V3 lo,V3 hi,float ang){
Brush b=mkbox(lo,hi);
V3 n={sinf(ang),0,cosf(ang)}; float d=n.x*((lo.x+hi.x)*0.5f)+n.z*hi.z; // through top centre
b.pl[4]={n,d}; return b;}
int main(){
float vs=0.224f;
// Dense junction: a flat slab + an overlapping slab whose top is tilted ~5deg.
// They overlap in x so a band of cells sees BOTH near-parallel top faces.
std::vector<Brush> B;
B.push_back(mkbox({-4,-3, -1},{ 1, 3, 0})); // A: flat floor slab, top z=0
B.push_back(mktilt({-1,-3, -1},{ 4, 3, 0}, 0.09f)); // B: overlaps A in x[-1,1], top tilted 5.2deg
// (0.09 rad = 5.16deg > 2.56deg dedup threshold -> NOT merged -> phantom crease)
V3 org{-5,-4,-2}; float ext[3]={10,8,4};
int G[3]; for(int a=0;a<3;a++)G[a]=(int)ceilf(ext[a]/vs)+1;
auto GI=[&](int i,int j,int k){return((size_t)k*G[1]+j)*G[0]+i;};
std::vector<float> grid((size_t)G[0]*G[1]*G[2]);
for(int k=0;k<G[2];k++)for(int j=0;j<G[1];j++)for(int i=0;i<G[0];i++)
grid[GI(i,j,k)]=usdf(org+V3{i*vs,j*vs,k*vs},B);
auto BR=[&](V3 p){V3 g=(p-org)*(1.0f/vs);int i=(int)floorf(g.x),j=(int)floorf(g.y),k=(int)floorf(g.z);
if(i<0||j<0||k<0||i>=G[0]-1||j>=G[1]-1||k>=G[2]-1)return 1e30f;float fx=g.x-i,fy=g.y-j,fz=g.z-k;
float c00=grid[GI(i,j,k)]*(1-fx)+grid[GI(i+1,j,k)]*fx,c10=grid[GI(i,j+1,k)]*(1-fx)+grid[GI(i+1,j+1,k)]*fx;
float c01=grid[GI(i,j,k+1)]*(1-fx)+grid[GI(i+1,j,k+1)]*fx,c11=grid[GI(i,j+1,k+1)]*(1-fx)+grid[GI(i+1,j+1,k+1)]*fx;
return (c00*(1-fy)+c10*fy)*(1-fz)+(c01*(1-fy)+c11*fy)*fz;};
int edir[3][3]={{1,0,0},{0,1,0},{0,0,1}};
int ce[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}};
for(int fix=0;fix<2;fix++){ // 0 = buggy (same-brush only), 1 = + hidden-face removal
int active=0, flew=0; std::vector<float> devs; std::vector<float> moved;
int red10=0, red17=0, outcell=0; std::vector<float> snapdev; // cap sweep + uncapped analysis
for(int k=0;k<G[2]-1;k++)for(int j=0;j<G[1]-1;j++)for(int i=0;i<G[0]-1;i++){
bool neg=false,pos=false;
for(int a=0;a<2;a++)for(int b2=0;b2<2;b2++)for(int c=0;c<2;c++){float d=grid[GI(i+a,j+b2,k+c)];if(d<0)neg=true;else pos=true;}
if(!(neg&&pos))continue;
// QEF start = the mass point of brick crossings (rounded), clamped to cell
V3 cmin=org+V3{i*vs,j*vs,k*vs}, cmax=cmin+V3{vs,vs,vs};
V3 mass{0,0,0}; int nc=0;
for(int e=0;e<12;e++){int ci=i+ce[e][0],cj=j+ce[e][1],ck=k+ce[e][2],dir=ce[e][3];
float d0=grid[GI(ci,cj,ck)],d1=grid[GI(ci+edir[dir][0],cj+edir[dir][1],ck+edir[dir][2])];
if((d0<0)==(d1<0))continue;
V3 p0=org+V3{ci*vs,cj*vs,ck*vs},p1=p0+V3{edir[dir][0]*vs,edir[dir][1]*vs,edir[dir][2]*vs};
float t=d0/(d0-d1); mass=mass+p0+(p1-p0)*t; nc++;}
if(!nc)continue; active++; V3 result=mass*(1.0f/nc);
V3 cell_ctr=cmin+V3{vs*0.5f,vs*0.5f,vs*0.5f};
// ---- qef_solve Step 5b plane-set capture (EXACT logic) ----
V3 cl_n[3]; float cl_d[3]; int cl_count=0; float T=vs*1.2f;
for(size_t bi=0;bi<B.size();bi++){ const Brush&br=B[bi];
for(size_t pi=0;pi<br.pl.size();pi++){ Plane pl=br.pl[pi];
float fd=dot(pl.n,cell_ctr)-pl.d; if(fabsf(fd)>=T)continue;
bool onbrush=true; // cell inside this brush's OTHER planes (~T)
for(size_t pj=0;pj<br.pl.size();pj++){ if(pj==pi)continue;
if(dot(br.pl[pj].n,cell_ctr)-br.pl[pj].d>T){onbrush=false;break;} }
if(!onbrush)continue;
if(fix){ // HIDDEN-FACE REMOVAL (csg4 ClipInside): F is a real union face
// only if a point just OUTSIDE F (across its normal) is outside every
// other brush. If another brush covers that point, F is buried -> skip.
V3 face_pt=cell_ctr-pl.n*fd; // project cell centre onto plane F
V3 q=face_pt+pl.n*(vs*0.5f); // step outward half a voxel
bool hidden=false;
for(size_t bk=0;bk<B.size();bk++){ if(bk==bi)continue;
if(bsdf(q,B[bk])<0.0f){hidden=true;break;} } // q inside Bk -> F covered
if(hidden)continue;
}
int matched=-1; for(int c=0;c<cl_count;c++) if(fabsf(dot(pl.n,cl_n[c]))>0.999f){matched=c;break;}
if(matched<0&&cl_count<3){cl_n[cl_count]=pl.n;cl_d[cl_count]=pl.d;cl_count++;}
}
}
V3 snapped=result;
if(cl_count==1) snapped=result-cl_n[0]*(dot(cl_n[0],result)-cl_d[0]);
else if(cl_count==2){ float cc=dot(cl_n[0],cl_n[1]),det=1-cc*cc;
if(det>1e-4f){float r0=cl_d[0]-dot(cl_n[0],result),r1=cl_d[1]-dot(cl_n[1],result);
float t0=(r0-cc*r1)/det,t1=(r1-cc*r0)/det; snapped=result+cl_n[0]*t0+cl_n[1]*t1;} }
else if(cl_count==3){ V3 a=cl_n[0],b3=cl_n[1],cn=cl_n[2]; float den=dot(a,cross3(b3,cn));
if(fabsf(den)>1e-5f) snapped=(cross3(b3,cn)*cl_d[0]+cross3(cn,a)*cl_d[1]+cross3(a,b3)*cl_d[2])*(1.0f/den); }
// ESCAPE FALLBACK: if a crease/corner snap leaves the cell, it's a far/degenerate
// intersection -> project onto the single NEAREST face instead (in-cell + on surface).
if(cl_count>=2){
V3 lo=cmin-V3{vs*0.5f,vs*0.5f,vs*0.5f}, hi=cmax+V3{vs*0.5f,vs*0.5f,vs*0.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){
int nb=0; float ba=1e30f;
for(int c=0;c<cl_count;c++){ float a2=fabsf(dot(cl_n[c],cell_ctr)-cl_d[c]); if(a2<ba){ba=a2;nb=c;} }
snapped=result-cl_n[nb]*(dot(cl_n[nb],result)-cl_d[nb]);
}
}
// clamp to cell (like the QEF) — final safety so nothing spikes
snapped.x=fmaxf(cmin.x,fminf(cmax.x,snapped.x));
snapped.y=fmaxf(cmin.y,fminf(cmax.y,snapped.y));
snapped.z=fmaxf(cmin.z,fminf(cmax.z,snapped.z));
float mv=0; if(cl_count>=1){ V3 dd=snapped-result; mv=sqrtf(dot(dd,dd)); }
V3 v=result;
if(cl_count>=1){ if(mv<vs*1.0f) v=snapped; else flew++; } // the 1-voxel snap-cap == red cell
if(cl_count>=1){
moved.push_back(mv);
if(mv>=vs*1.0f) red10++; // red at cap 1.0 (current)
if(mv>=vs*1.7f) red17++; // red at cap 1.7 (~cell diagonal)
// if we ALWAYS accepted the snap (uncapped), is it correct + inside the cell?
snapdev.push_back(fabsf(usdf(snapped,B))); // surf-dev of the UNCAPPED snap
V3 lo=cmin-V3{vs*0.5f,vs*0.5f,vs*0.5f}, hi=cmax+V3{vs*0.5f,vs*0.5f,vs*0.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)
outcell++; // snapped vertex escaped the cell (+half) = spike
}
devs.push_back(fabsf(usdf(v,B))); // dev from true union surface
}
std::sort(moved.begin(),moved.end()); std::sort(devs.begin(),devs.end()); std::sort(snapdev.begin(),snapdev.end());
printf("%-22s active=%d red@cap1.0=%d(%.1f%%) red@cap1.7=%d(%.1f%%) UNCAPPED: surf-dev p95=%.4f MAX=%.4f out-of-cell(spikes)=%d(%.1f%%) (voxel=%.3f)\n",
fix?"+ hidden-face fix":"buggy (same-brush)", active,
red10, active?100.0f*red10/active:0, red17, active?100.0f*red17/active:0,
snapdev.empty()?0:snapdev[(size_t)(snapdev.size()*0.95)], snapdev.empty()?0:snapdev.back(),
outcell, active?100.0f*outcell/active:0, vs);
}
return 0;
}