// Reproduce the THIN-SLAB "one face smooth, one face stair-stepped" asymmetry — CPU oracle.
//
// A tilted slab has two big PARALLEL faces (opposite normals). When the slab gets thin
// (thickness within ~2*T of the plane-set window), a cell near one face also captures the
// OTHER face. The bake dedups by |dot(n, cl_n)| > 0.999 — which also matches ANTI-parallel
// normals — so the two faces collapse to ONE plane (the first-captured). Cells then snap to
// the wrong face → one face is fine, the other aliases (stair-steps). This is distinct from
// the angled-crease sawtooth (that one is pure resolution).
//
// Measures, per face (front/back, classified by the face the cell centre is nearest), the
// deviation of the final DC vertex from THAT face's plane. Low = on its face (smooth);
// ~thickness = snapped to the other face (collapsed/stairstep).
#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 rotY(V3 v, float a) { float c = cosf(a), s = sinf(a); return {c * v.x + s * v.z, v.y, -s * v.x + c * v.z}; }
struct Plane { V3 n; float d; };
struct Brush { std::vector<Plane> pl; int mat; };
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; }
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 slab(V3 lo, V3 hi, float a, int mat) { // box rotated about Y (through origin)
Brush b; b.mat = mat;
b.pl = { {rotY({1,0,0}, a), hi.x}, {rotY({-1,0,0}, a), -lo.x},
{{0,1,0}, hi.y}, {{0,-1,0}, -lo.y},
{rotY({0,0,1}, a), hi.z}, {rotY({0,0,-1}, a), -lo.z} };
return b;
}
static void activeFace(V3 p, const std::vector<Brush>& B, V3& n) {
float best = 1e30f; int bi = 0; for (size_t i = 0; i < B.size(); i++) { float s = bsdf(p, B[i]); if (s < best) { best = s; bi = (int)i; } }
float bp = -1e30f; n = {0,0,1}; for (auto& q : B[bi].pl) { float dp = dot(q.n, p) - q.d; if (dp > bp) { bp = dp; n = q.n; } }
}
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 = b + m * d; 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 int edir[3][3] = {{1,0,0},{0,1,0},{0,0,1}};
static 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}};
// Deployed method = plane-set snap. `dedup_antiparallel`: if false, DON'T merge anti-parallel
// faces (the candidate fix) so a thin slab keeps both faces distinct.
static void run(const char* label, float vs, float thick, float ang, bool dedup_antiparallel) {
Brush s = slab({-2, -2, -thick * 0.5f}, {2, 2, thick * 0.5f}, ang, 0);
std::vector<Brush> B = {s};
// The two big faces (±z rotated). Classify a cell to the nearer one; measure dev to it.
V3 nf = rotY({0,0,1}, ang); float df = thick * 0.5f; // front face n·p = df
V3 nb = rotY({0,0,-1}, ang); float db = thick * 0.5f; // back face n·p = db
V3 org{-3,-3,-3}; int G = (int)ceilf(6.0f / vs) + 1;
auto D = [&](int i, int j, int k) { return usdf(org + V3{i*vs, j*vs, k*vs}, B); };
std::vector<float> devF, devB;
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;
V3 cmin=org+V3{i*vs,j*vs,k*vs}, cmax=cmin+V3{vs,vs,vs}, cc=cmin+V3{vs*.5f,vs*.5f,vs*.5f};
// Measure only cells firmly on a BIG face — skip cells near the slab's thin ±x/±y
// faces (edge creases, which are resolution-bound and not what this oracle isolates).
float xr = dot(rotY({1,0,0}, ang), cc);
if (fabsf(xr) > 2.0f - 2.0f * vs) continue;
if (fabsf(cc.y) > 2.0f - 2.0f * vs) continue;
std::vector<V3> P, N;
for (int e=0;e<12;e++){int ci=i+ce[e][0],cj=j+ce[e][1],ck=k+ce[e][2],dr=ce[e][3];
float d0=D(ci,cj,ck),d1=D(ci+edir[dr][0],cj+edir[dr][1],ck+edir[dr][2]); if((d0<0)==(d1<0))continue;
V3 p0=org+V3{ci*vs,cj*vs,ck*vs},p1=p0+V3{edir[dr][0]*vs,edir[dr][1]*vs,edir[dr][2]*vs};
float t=d0/(d0-d1); V3 ps=p0+(p1-p0)*t; V3 nm; activeFace(ps,B,nm); P.push_back(ps); N.push_back(nm);}
if (P.empty()) continue;
V3 v = qef(P, N, cmin, cmax, 0.01f);
// plane-set snap (mode A), with the dedup rule under test
V3 cl_n[3]; float cl_d[3]; int cl=0; float T=vs*1.2f;
for (size_t bi=0; bi<B.size(); bi++) for (size_t pi=0; pi<B[bi].pl.size(); pi++) { Plane pl=B[bi].pl[pi];
float fd=dot(pl.n,cc)-pl.d; if(fabsf(fd)>=T)continue; bool ob=true;
for(size_t pj=0;pj<B[bi].pl.size();pj++){if(pj==pi)continue;if(dot(B[bi].pl[pj].n,cc)-B[bi].pl[pj].d>T){ob=false;break;}}
if(!ob)continue;
int mm=-1; for(int q=0;q<cl;q++){ float dp=dot(pl.n,cl_n[q]);
bool par = dedup_antiparallel ? (fabsf(dp)>0.999f) : (dp>0.999f); // FIX: only merge SAME-facing
if(par){mm=q;break;} }
if(mm<0&&cl<3){cl_n[cl]=pl.n;cl_d[cl]=pl.d;cl++;} }
V3 sn=v;
if(cl==1)sn=v-cl_n[0]*(dot(cl_n[0],v)-cl_d[0]);
else if(cl==2){float cco=dot(cl_n[0],cl_n[1]),det=1-cco*cco;if(det>1e-4f){float r0=cl_d[0]-dot(cl_n[0],v),r1=cl_d[1]-dot(cl_n[1],v);sn=v+cl_n[0]*((r0-cco*r1)/det)+cl_n[1]*((r1-cco*r0)/det);}}
if(cl>=2){V3 lo=cmin-V3{vs*.5f,vs*.5f,vs*.5f},hi=cmax+V3{vs*.5f,vs*.5f,vs*.5f};
if(sn.x<lo.x||sn.y<lo.y||sn.z<lo.z||sn.x>hi.x||sn.y>hi.y||sn.z>hi.z){int nb2=0;float ba=1e30f;for(int q=0;q<cl;q++){float aa=fabsf(dot(cl_n[q],cc)-cl_d[q]);if(aa<ba){ba=aa;nb2=q;}}sn=v-cl_n[nb2]*(dot(cl_n[nb2],v)-cl_d[nb2]);}}
sn.x=fmaxf(cmin.x,fminf(cmax.x,sn.x));sn.y=fmaxf(cmin.y,fminf(cmax.y,sn.y));sn.z=fmaxf(cmin.z,fminf(cmax.z,sn.z));
if(cl>=1)v=sn;
// classify to the nearer face and record deviation to it
float ff = fabsf(dot(nf, cc) - df), fb = fabsf(dot(nb, cc) - db);
if (ff <= fb) devF.push_back(fabsf(dot(nf, v) - df));
else devB.push_back(fabsf(dot(nb, v) - db));
}
auto p95 = [](std::vector<float>& d) { if (d.empty()) return 0.0f; std::sort(d.begin(), d.end()); return d[(size_t)(d.size() * 0.95f)]; };
printf(" %-30s vox=%.3f thick=%.2f(%.1fvox) frontP95=%.4f backP95=%.4f\n",
label, vs, thick, thick / vs, p95(devF), p95(devB));
}
int main() {
float a = 0.4f; // tilt so faces are off-grid (the hard case)
printf("Thin-slab per-face deviation (tilt=%.1f rad). Low = on its own face (smooth).\n", a);
printf("-- current dedup (merges anti-parallel too) --\n");
run("thick slab", 0.224f, 1.20f, a, true);
run("medium slab", 0.224f, 0.50f, a, true);
run("thin slab", 0.224f, 0.25f, a, true);
printf("-- FIX: don't merge anti-parallel faces --\n");
run("thick slab", 0.224f, 1.20f, a, false);
run("medium slab", 0.224f, 0.50f, a, false);
run("thin slab", 0.224f, 0.25f, a, false);
printf("-- resolution (thin slab, current dedup) --\n");
run("thin @2x", 0.112f, 0.25f, a, true);
return 0;
}