// Faithful mini-pipeline repro: why do ROTATED (oblique) brush faces come out
// bumpy while AXIS-aligned faces are razor-sharp, under exact-Hermite + low-eps QEF?
//
// Mirrors the shader path: sample a convex cube's SDF on a voxel grid, classify
// edge sign changes, compute EXACT ray-vs-plane Hermite per active edge (argmax
// active face at the crossing), run the qef_solve QEF per cell, then measure how
// far each face-cell vertex deviates from the TRUE cube face (flatness).
#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;}
struct Plane{V3 n;float d;};
// Rotate a vector by Euler (mimics Basis.from_euler in the GD test).
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);
// Z*Y*X
V3 r;
float x=v.x,y=v.y,z=v.z;
// Rx
float y1=cx*y - sx*z, z1=sx*y + cx*z; float x1=x;
// Ry
float x2=cy*x1 + sy*z1, z2=-sy*x1 + cy*z1; float y2=y1;
// Rz
r.x=cz*x2 - sz*y2; r.y=sz*x2 + cz*y2; r.z=z2;
return r;
}
static std::vector<Plane> make_cube(V3 c,float h,bool rotate){
V3 axes[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:axes){
V3 n = rotate? rot(a,0.35f,0.6f,0.0f) : a;
p.push_back({n, dot(n,c)+h}); // inside: n.x<=d
}
return p;
}
static float sdf(V3 p,const std::vector<Plane>&pl){ // convex: max(n.p-d), neg inside
float m=-1e30f; for(auto&q:pl) m=std::max(m,dot(q.n,p)-q.d); return m;
}
// qef_solve.glsl 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;
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 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;
}
// PLANE-SNAP (qef_solve Step 5b), now fed EXACT normals: cluster the cell's
// normals into <=3 distinct planes and project the (clamped) QEF vertex onto the
// face / crease / corner. Safety net: reject snaps that move it > 2 voxels.
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};}
static V3 plane_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; }}
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;
}
// distance from a point to the nearest cube face plane it's ~on (flatness metric)
static float face_dev(V3 p,const std::vector<Plane>&pl){
float best=1e30f; for(auto&q:pl){float e=fabsf(dot(q.n,p)-q.d); best=std::min(best,e);} return best;
}
static void run(bool rotate,float eps,float noise,bool diag=false,bool snap=false){
V3 c{0,0,0}; float h=5.0f; auto cube=make_cube(c,h,rotate);
float vs=0.224f; // voxel size (deployed)
// grid spanning the cube +margin
int G=(int)ceilf((2*h+4)/vs);
V3 org{-h-2,-h-2,-h-2};
auto D=[&](int i,int j,int k){V3 p=org+V3{i*vs,j*vs,k*vs};
float s=sdf(p,cube);
if(noise>0){ // deterministic pseudo-noise ~ detail noise on the surface
float nn=sinf(p.x*12.3f)*sinf(p.y*11.1f)*sinf(p.z*13.7f);
s+=nn*noise; }
return s; };
int edir[3][3]={{1,0,0},{0,1,0},{0,0,1}};
// per cell, gather exact-Hermite from its 12 edges, QEF, measure flatness
std::vector<float> devs, devs_unsorted; std::vector<int> planecount; int hardcells=0; float maxdev=0;
for(int k=0;k<G-1;k++)for(int j=0;j<G-1;j++)for(int i=0;i<G-1;i++){
// 8 corners
float dv[2][2][2];
for(int a=0;a<2;a++)for(int b=0;b<2;b++)for(int cc2=0;cc2<2;cc2++)
dv[a][b][cc2]=D(i+a,j+b,k+cc2);
// is this an active cell? any corner sign differs
bool anyneg=false,anypos=false;
for(int a=0;a<2;a++)for(int b=0;b<2;b++)for(int cc2=0;cc2<2;cc2++){
if(dv[a][b][cc2]<0)anyneg=true; else anypos=true;}
if(!(anyneg&&anypos))continue;
// 12 edges: enumerate as (corner, dir)
std::vector<V3> P,N; bool hard=false;
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 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};
V3 p1=p0+V3{edir[dir][0]*vs,edir[dir][1]*vs,edir[dir][2]*vs};
float d0=D(ci,cj,ck);
float d1=D(ci+edir[dir][0],cj+edir[dir][1],ck+edir[dir][2]);
if((d0<0)==(d1<0))continue; // no sign change
// binary search crossing on sampled field (8 iters like shader)
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(pm,cube)+(noise>0? sinf(pm.x*12.3f)*sinf(pm.y*11.1f)*sinf(pm.z*13.7f)*noise:0);
if((dlo<0)==(dm<0)){lo=mid;dlo=dm;}else hi=mid;}
float t=(lo+hi)*0.5f; V3 pos=p0+(p1-p0)*t;
// argmax active face at pos (within half a voxel => hard)
float best=-1e30f; Plane bpl{{0,0,0},0};
for(auto&q:cube){float dp=dot(q.n,pos)-q.d; if(dp>best){best=dp;bpl=q;}}
if(fabsf(best) < vs*0.5f){
// EXACT ray-plane
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); hard=true;
} else {
// fallback: sampled pos + central-diff normal
float eps2=vs*0.01f;
V3 nm{sdf(pos+V3{eps2,0,0},cube)-sdf(pos-V3{eps2,0,0},cube),
sdf(pos+V3{0,eps2,0},cube)-sdf(pos-V3{0,eps2,0},cube),
sdf(pos+V3{0,0,eps2},cube)-sdf(pos-V3{0,0,eps2},cube)};
float l=sqrtf(dot(nm,nm)); if(l>1e-8f)nm=nm*(1.0f/l);
P.push_back(pos); N.push_back(nm);
}
}
if(P.empty())continue;
float use_eps = hard? eps : 0.01f;
V3 cmin=org+V3{i*vs,j*vs,k*vs}, cmax=cmin+V3{vs,vs,vs};
V3 v=qef(P,N,cmin,cmax,use_eps);
if(snap&&hard) v=plane_snap(v,P,N,vs);
if(hard){hardcells++; float dv2=face_dev(v,cube); devs.push_back(dv2); devs_unsorted.push_back(dv2); maxdev=std::max(maxdev,dv2);
int np=0; V3 seen[8];
for(auto&nm:N){bool m=false; for(int q=0;q<np;q++) if(fabsf(dot(nm,seen[q]))>0.99f){m=true;break;} if(!m&&np<8)seen[np++]=nm;}
planecount.push_back(np);
}
}
// breakdown: for each distinct-plane bucket, count + max dev (1=face,2=edge,3=corner)
if(diag){
int cnt[5]={0,0,0,0,0}; float mx[5]={0,0,0,0,0}; int single_bad=0;
for(size_t z=0;z<planecount.size();z++){int b=planecount[z]; if(b>4)b=4; cnt[b]++; mx[b]=std::max(mx[b],devs_unsorted[z]);
if(planecount[z]==1&&devs_unsorted[z]>0.05f)single_bad++;}
printf(" breakdown by #distinct planes: face(1) n=%d maxdev=%.4f | edge(2) n=%d maxdev=%.4f | corner(3) n=%d maxdev=%.4f | 4+ n=%d maxdev=%.4f\n",
cnt[1],mx[1],cnt[2],mx[2],cnt[3],mx[3],cnt[4],mx[4]);
printf(" *** SINGLE-plane FACE cells with dev>0.05 (these SHOULD be flat): %d ***\n", single_bad);
}
std::sort(devs.begin(),devs.end());
float p50=devs.empty()?0:devs[devs.size()/2];
float p95=devs.empty()?0:devs[(size_t)(devs.size()*0.95)];
printf("%-8s eps=%-7g noise=%-5g hardcells=%-5d face-dev p50=%.4f p95=%.4f MAX=%.4f (voxel=%.3f)\n",
snap?(rotate?"ROT+SNAP":"AX+SNAP"):(rotate?"ROTATED":"AXIS"), eps, noise, hardcells, p50, p95, maxdev, vs);
}
int main(){
printf("Flatness = distance of each hard-cell vertex from the true cube face (0 = perfect).\n\n");
printf("-- no noise --\n");
run(false,1e-4f,0.0f); run(true,1e-4f,0.0f,true);
run(false,0.01f,0.0f); run(true,0.01f,0.0f);
printf("\n-- with surface noise (0.05, ~mimics detail noise) --\n");
run(false,1e-4f,0.05f); run(true,1e-4f,0.05f);
run(false,0.01f,0.05f); run(true,0.01f,0.05f);
printf("\n-- EXACT HERMITE + SNAP (project onto exact plane/crease/corner) --\n");
run(true,1e-4f,0.0f,false,true); run(true,1e-4f,0.05f,false,true);
run(true,0.01f,0.0f,false,true); run(true,0.01f,0.05f,false,true);
run(false,1e-4f,0.05f,false,true);
return 0;
}