// STEP 2 (B) oracle: does exact-Hermite override + snap still sharpen a face when
// the density comes from a TRILINEAR BRICK (pre-sampled scalar grid), not the
// direct SDF? This is the dense-chunk path: brick supplies density/topology, the
// override supplies exact planes. If flat -> the 2-line change works.
//
// Compares, on a ROTATED cube:
// brick-only : trilinear brick density, NO override, NO snap (= today's soup)
// brick + override+snap : find crossing on brick, snap pos to exact plane, snap vertex
#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;};
static V3 rot(V3 v,float ax,float ay,float az){
float cx=cosf(ax),sx=sinf(ax),cy=cosf(ay),sy=sinf(ay),cz=cosf(az),sz=sinf(az);
float y1=cx*v.y-sx*v.z,z1=sx*v.y+cx*v.z,x1=v.x;
float x2=cy*x1+sy*z1,z2=-sy*x1+cy*z1,y2=y1;
return{cz*x2-sz*y2,sz*x2+cz*y2,z2};
}
static std::vector<Plane> cube(float h,bool r){
V3 ax[6]={{1,0,0},{-1,0,0},{0,1,0},{0,-1,0},{0,0,1},{0,0,-1}};std::vector<Plane>p;
for(auto a:ax){V3 n=r?rot(a,0.35f,0.6f,0):a;p.push_back({n,h});}return p;}
static float sdf(V3 p,const std::vector<Plane>&pl){float m=-1e30f;for(auto&q:pl)m=fmaxf(m,dot(q.n,p)-q.d);return m;}
static V3 qef(const std::vector<V3>&P,const std::vector<V3>&N,V3 cmin,V3 cmax,float eps){
float a[6]={0,0,0,0,0,0};V3 b{0,0,0},mass{0,0,0};int n=0;
for(size_t i=0;i<P.size();i++){V3 m=N[i];a[0]+=m.x*m.x;a[1]+=m.x*m.y;a[2]+=m.x*m.z;a[3]+=m.y*m.y;a[4]+=m.y*m.z;a[5]+=m.z*m.z;
float d=dot(m,P[i]);b.x+=d*m.x;b.y+=d*m.y;b.z+=d*m.z;mass=mass+P[i];n++;}
if(!n)return mass;mass=mass*(1.0f/n);
float m00=a[0]+eps,m01=a[1],m02=a[2],m11=a[3]+eps,m12=a[4],m22=a[5]+eps;
V3 amp{a[0]*mass.x+a[1]*mass.y+a[2]*mass.z,a[1]*mass.x+a[3]*mass.y+a[4]*mass.z,a[2]*mass.x+a[4]*mass.y+a[5]*mass.z};
V3 rhs=b-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=fmaxf(cmin.x,fminf(cmax.x,r.x));r.y=fmaxf(cmin.y,fminf(cmax.y,r.y));r.z=fmaxf(cmin.z,fminf(cmax.z,r.z));return r;}
static V3 snap(V3 v,const std::vector<V3>&P,const std::vector<V3>&N,float vs){
V3 cn[3];float cd[3];int nc=0;
for(size_t i=0;i<N.size();i++){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)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]),det=1-cc*cc;if(det>1e-4f){float r0=cd[0]-dot(cn[0],v),r1=cd[1]-dot(cn[1],v),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;return sqrtf(dx*dx+dy*dy+dz*dz)<vs*2.0f?s:v;}
static float fdev(V3 p,const std::vector<Plane>&pl){float b=1e30f;for(auto&q:pl)b=fminf(b,fabsf(dot(q.n,p)-q.d));return b;}
int main(){
float h=5,vs=0.224f;auto C=cube(h,true);
V3 org{-h-2,-h-2,-h-2};int G=(int)ceilf((2*h+4)/vs)+1;
// BAKE the brick: sample SDF at every grid point (this is brick_delta, here the
// full density since base=0 for an isolated cube in a carved void).
std::vector<float> grid((size_t)G*G*G);
auto GI=[&](int i,int j,int k){return((size_t)k*G+j)*G+i;};
for(int k=0;k<G;k++)for(int j=0;j<G;j++)for(int i=0;i<G;i++)grid[GI(i,j,k)]=sdf(org+V3{i*vs,j*vs,k*vs},C);
// trilinear brick sample (what the mesh reads on a dense chunk)
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-1||j>=G-1||k>=G-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 mode=0;mode<2;mode++){ // 0=brick-only, 1=brick+override+snap
std::vector<float> dev;
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=grid[GI(i+a,j+b,k+c)];if(d<0)neg=true;else pos=true;}
if(!(neg&&pos))continue;std::vector<V3>P,N;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=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;
// binary search on the TRILINEAR BRICK (what the mesh does)
float lo=0,hi=1,dlo=d0;for(int it=0;it<8;it++){float mid=(lo+hi)*.5f;float dm=BR(p0+(p1-p0)*mid);if((dlo<0)==(dm<0)){lo=mid;dlo=dm;}else hi=mid;}
float t=(lo+hi)*.5f;V3 ps=p0+(p1-p0)*t;
V3 nm;// normal via brick central diff
float eps=vs*.01f;nm={BR(ps+V3{eps,0,0})-BR(ps-V3{eps,0,0}),BR(ps+V3{0,eps,0})-BR(ps-V3{0,eps,0}),BR(ps+V3{0,0,eps})-BR(ps-V3{0,0,eps})};
float l=sqrtf(dot(nm,nm));if(l>1e-8f)nm=nm*(1/l);
if(mode==1){// override: argmax active face at ps, exact ray-plane if within 0.5vox
float best=-1e30f;Plane bp{{0,0,0},0};for(auto&q:C){float dp=dot(q.n,ps)-q.d;if(dp>best){best=dp;bp=q;}}
if(fabsf(best)<vs*0.5f){V3 ed=p1-p0;float den=dot(bp.n,ed);if(fabsf(den)>1e-8f){float te=fmaxf(0.f,fminf(1.f,(bp.d-dot(bp.n,p0))/den));ps=p0+ed*te;}
nm=dot(bp.n,nm)<0?bp.n*-1.f:bp.n;hard=true;}}
P.push_back(ps);N.push_back(nm);}
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(mode==1&&hard)v=snap(v,P,N,vs);
dev.push_back(fdev(v,C));}
std::sort(dev.begin(),dev.end());
printf("%-24s cells=%zu face-dev p50=%.4f p95=%.4f MAX=%.4f (voxel=%.3f)\n",
mode?"brick+override+snap":"brick-only (today)",dev.size(),
dev.empty()?0:dev[dev.size()/2],dev.empty()?0:dev[(size_t)(dev.size()*.95)],dev.empty()?0:dev.back(),vs);
}
return 0;
}