// STEP 2a oracle: dense abutting brushes. Two boxes forming an L-beam share a
// concave edge (the notch). Question: does CROSS-BRUSH clustering (drop the
// per-brush confinement) snap that shared edge sharp, where PER-BRUSH confinement
// (current shader) rounds it? And does it REGRESS two separate close boxes
// (must NOT merge — the reason confinement was added)?
//
// Faithful mini-pipeline: grid-sample union SDF -> exact Hermite per crossing
// (tagged with the brush whose face it lies on, argmin|f| across brushes) -> QEF
// -> snap in one of two modes. Measures concave-edge deviation from the true
// crease line, and separate-box "bridge" leakage.
#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;};
struct Box{std::vector<Plane> pl;};
static Box make_box(V3 c,V3 h){
Box b; V3 ax[6]={{1,0,0},{-1,0,0},{0,1,0},{0,-1,0},{0,0,1},{0,0,-1}};
float hs[6]={h.x,h.x,h.y,h.y,h.z,h.z};
for(int i=0;i<6;i++) b.pl.push_back({ax[i], dot(ax[i],c)+hs[i]});
return b;
}
static float sdf_box(V3 p,const Box&b){float m=-1e30f;for(auto&q:b.pl)m=std::max(m,dot(q.n,p)-q.d);return m;}
static float sdf_union(V3 p,const std::vector<Box>&bs){float m=1e30f;for(auto&b:bs)m=std::min(m,sdf_box(p,b));return m;} // neg inside
// QEF (qef_solve Step 3): regularise toward mass, Cramer, clamp to cell.
static V3 qef(const std::vector<V3>&P,const std::vector<V3>&N,V3 cmin,V3 cmax,float eps){
float ata[6]={0,0,0,0,0,0};V3 atb{0,0,0},mass{0,0,0};int n=0;
for(size_t i=0;i<P.size();i++){V3 nm=N[i];
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;
float d=dot(nm,P[i]);atb.x+=d*nm.x;atb.y+=d*nm.y;atb.z+=d*nm.z;mass=mass+P[i];n++;}
if(!n)return mass;mass=mass*(1.0f/n);
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,c11=m00*m22-m02*m02,c12=m01*m02-m00*m12,c22=m00*m11-m01*m01;
float det=m00*c00+m01*c01+m02*c02;V3 r;
if(fabsf(det)>1e-12f){float iv=1.0f/det;r=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};}else r=mass;
r.x=std::max(cmin.x,std::min(cmax.x,r.x));r.y=std::max(cmin.y,std::min(cmax.y,r.y));r.z=std::max(cmin.z,std::min(cmax.z,r.z));
return r;
}
// snap: mode 0 = per-brush confinement (only first brush's planes), mode 1 = cross-brush (all)
static V3 snap(V3 v,const std::vector<V3>&P,const std::vector<V3>&N,const std::vector<int>&B,float vs,int mode){
V3 cn[3];float cd[3];int nc=0;int confine=-1;
for(size_t i=0;i<N.size();i++){
if(mode==0){ if(confine<0)confine=B[i]; if(B[i]!=confine)continue; }
int m=-1;for(int q=0;q<nc;q++)if(fabsf(dot(N[i],cn[q]))>0.999f){m=q;break;}
if(m<0){if(nc<3){cn[nc]=N[i];cd[nc]=dot(N[i],P[i]);nc++;}else return v;}
}
if(nc==0)return v;
V3 s=v;
if(nc==1)s=v-cn[0]*(dot(cn[0],v)-cd[0]);
else if(nc==2){float cc=dot(cn[0],cn[1]);float det=1-cc*cc;if(det>1e-4f){float r0=cd[0]-dot(cn[0],v),r1=cd[1]-dot(cn[1],v);float t0=(r0-cc*r1)/det,t1=(r1-cc*r0)/det;s=v+cn[0]*t0+cn[1]*t1;}}
else{V3 a=cn[0],b=cn[1],c=cn[2];float den=dot(a,cross3(b,c));if(fabsf(den)>1e-5f)s=(cross3(b,c)*cd[0]+cross3(c,a)*cd[1]+cross3(a,b)*cd[2])*(1.0f/den);}
float dx=s.x-v.x,dy=s.y-v.y,dz=s.z-v.z;
if(sqrtf(dx*dx+dy*dy+dz*dz)<vs*2.0f)return s;
return v;
}
// process the grid; returns nothing but fills callback via measuring closure
static void run(const char*name,const std::vector<Box>&bs,int mode,
float creaseX,float creaseY,float zlo,float zhi, // concave-edge line x=creaseX,y=creaseY
float* out_edge_p95,float* out_edge_max,int* out_edge_n){
float vs=0.224f;
V3 org{-3,-3,-1}; int G=(int)ceilf((16.0f)/vs);
auto D=[&](int i,int j,int k){return sdf_union(org+V3{i*vs,j*vs,k*vs},bs);};
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}};
std::vector<float> edev;
for(int k=0;k<G-1;k++)for(int j=0;j<G-1;j++)for(int i=0;i<G-1;i++){
bool neg=false,pos=false;
for(int a=0;a<2;a++)for(int b=0;b<2;b++)for(int c=0;c<2;c++){float d=D(i+a,j+b,k+c);if(d<0)neg=true;else pos=true;}
if(!(neg&&pos))continue;
std::vector<V3>P,N;std::vector<int>Bid;bool hard=false;
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];
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 d0=D(ci,cj,ck),d1=D(ci+edir[dir][0],cj+edir[dir][1],ck+edir[dir][2]);
if((d0<0)==(d1<0))continue;
float lo=0,hi=1,dlo=d0;
for(int it=0;it<8;it++){float mid=(lo+hi)*0.5f;V3 pm=p0+(p1-p0)*mid;float dm=sdf_union(pm,bs);if((dlo<0)==(dm<0)){lo=mid;dlo=dm;}else hi=mid;}
float t=(lo+hi)*0.5f;V3 pos=p0+(p1-p0)*t;
// cross-brush: pick the brush whose face is closest (argmin |maxplane|), and its plane
float bestf=1e30f;int bb=-1;Plane bpl{{0,0,0},0};
for(int bi=0;bi<(int)bs.size();bi++){float best=-1e30f;Plane bp{{0,0,0},0};
for(auto&q:bs[bi].pl){float dp=dot(q.n,pos)-q.d;if(dp>best){best=dp;bp=q;}}
if(fabsf(best)<bestf){bestf=fabsf(best);bb=bi;bpl=bp;}}
if(bestf<vs*0.5f){V3 ed=p1-p0;float den=dot(bpl.n,ed);if(fabsf(den)>1e-8f){float te=(bpl.d-dot(bpl.n,p0))/den;te=std::max(0.f,std::min(1.f,te));pos=p0+ed*te;}
P.push_back(pos);N.push_back(bpl.n);Bid.push_back(bb);hard=true;}
}
if(P.empty())continue;
V3 cmin=org+V3{i*vs,j*vs,k*vs},cmax=cmin+V3{vs,vs,vs};
V3 v=qef(P,N,cmin,cmax,0.01f);
if(hard)v=snap(v,P,N,Bid,vs,mode);
// concave-edge metric: cells whose true nearest surface is the crease line
// (both creaseX and creaseY within the cell, z in range)
if(cmin.x<=creaseX&&creaseX<=cmax.x && cmin.y<=creaseY&&creaseY<=cmax.y && v.z>zlo&&v.z<zhi){
float dx=v.x-creaseX,dy=v.y-creaseY;edev.push_back(sqrtf(dx*dx+dy*dy));
}
}
std::sort(edev.begin(),edev.end());
*out_edge_n=(int)edev.size();
*out_edge_p95=edev.empty()?0:edev[(size_t)(edev.size()*0.95)];
*out_edge_max=edev.empty()?0:edev.back();
printf("%-26s mode=%s concave-edge cells=%d dev-from-crease p95=%.4f MAX=%.4f (voxel=%.3f)\n",
name, mode?"CROSS-BRUSH":"per-brush ", *out_edge_n, *out_edge_p95, *out_edge_max, vs);
}
int main(){
printf("L-beam junction: two abutting boxes share the concave edge at x=2,y=2.\n");
printf("Good = vertex ON the crease (dev ~0). per-brush confinement should ROUND it.\n\n");
// L-beam: horizontal arm A (x0..10,y0..2) + vertical arm B (x0..2,y0..10), z0..6
std::vector<Box> L = { make_box({5,1,3},{5,1,3}), make_box({1,5,3},{1,5,3}) };
float p95,mx;int n;
run("L-beam junction", L, 0, 2.0f,2.0f, 0.5f,5.5f, &p95,&mx,&n);
run("L-beam junction", L, 1, 2.0f,2.0f, 0.5f,5.5f, &p95,&mx,&n);
printf("\nRegression: two SEPARATE boxes with a %.2f-voxel gap must NOT merge.\n", 0.6f);
printf("(measure the 'crease' between their FACING faces at x=2.3 mid-gap; a real\n merge would place vertices there — dev-from-that-plane near 0 = BAD leak)\n\n");
// Box C: x[-2,2], Box D: x[2.6,6] (gap 0.6 in x), both y[0,8] z[0,6]
std::vector<Box> S = { make_box({0,4,3},{2,4,3}), make_box({4.3f,4,3},{1.7f,4,3}) };
// there is no true crease between them; measure cells in the gap region x~2.3
// by reusing the metric against a fake crease line at (2.3, 4) — vertices there = leak
run("separate boxes (gap)", S, 1, 2.3f,4.0f, 0.5f,5.5f, &p95,&mx,&n);
printf(" ^ CROSS-BRUSH on separate boxes: cells near mid-gap should be FEW/none\n (no surface there); if many with low dev, the safety net failed to protect.\n");
return 0;
}