// ISOLATE the persisting sawtooth. Hypothesis: along a crease (two brush faces
// meeting at an edge), DC cells ALTERNATE between "captures both planes" (vertex
// snaps to the crease, correct) and "captures only one plane" (vertex on the face,
// off the crease → sticks out → sawtooth). If so, it's a feature-edge TOPOLOGY
// problem — NOT fixable by min_absf or the plane-set brick (both fix positions).
//
// Model a convex wedge = two planes meeting at a 90° crease (like ceiling∧wall).
// Grid-sample the exact union SDF, mesh (per active cell: gather edge-crossing
// normals, cluster into distinct planes, snap), and for the cells ALONG the crease
// report the distribution of distinct-plane-count and each vertex's deviation from
// the true crease LINE.
#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 rot(V3 v,float a){ // rotate about Z by a (tilt the crease off-grid)
float c=cosf(a),s=sinf(a); return {c*v.x-s*v.y, s*v.x+c*v.y, v.z};
}
struct Plane{V3 n;float d;};
int main(){
float vs=0.224f, tilt=0.4f; // tilt so the crease is NOT grid-aligned (worst case)
// Wedge: solid where BOTH n1·p<=d1 AND n2·p<=d2 (a convex corner). Crease at
// the line where both planes meet. Two faces: top (n1) and side (n2), tilted.
Plane p1{ rot({0,1,0}, tilt), 0.0f }; // "ceiling" (normal ~+y), through origin
Plane p2{ rot({1,0,0}, tilt), 0.0f }; // "wall" (normal ~+x), through origin
auto sdf=[&](V3 p){ return std::max(dot(p1.n,p)-p1.d, dot(p2.n,p)-p2.d); }; // convex, neg inside
V3 org{-3,-3,-1}; int G=(int)ceilf(6.0f/vs)+1;
auto D=[&](int i,int j,int k){return sdf(org+V3{i*vs,j*vs,k*vs});};
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}};
// crease line: intersection of p1 and p2 planes (both through origin) = the z-axis
// rotated... both planes pass through origin, their intersection is a line. Since
// p1.n and p2.n both lie in the xy-plane (z-component 0), the crease line is the
// z-axis direction through origin: points (0,0,t). Deviation from crease = sqrt(x^2+y^2).
int plane_hist[5]={0,0,0,0,0};
std::vector<float> crease_dev;
int crease_cells=0, off_crease=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++){
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;
// gather crossing normals (exact plane the crossing lies on = argmax active face)
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],dir=ce[e][3];
V3 a=org+V3{ci*vs,cj*vs,ck*vs},b=a+V3{edir[dir][0]*vs,edir[dir][1]*vs,edir[dir][2]*vs};
float da=D(ci,cj,ck),db=D(ci+edir[dir][0],cj+edir[dir][1],ck+edir[dir][2]);
if((da<0)==(db<0))continue;
float t=da/(da-db); V3 ps=a+(b-a)*t;
// exact active plane at ps (the wedge's argmax face)
V3 nrm = (dot(p1.n,ps)-p1.d > dot(p2.n,ps)-p2.d) ? p1.n : p2.n;
P.push_back(ps); N.push_back(nrm);
}
if(P.empty())continue;
// is this a crease cell? (cell straddles the crease line: x and y both span 0)
V3 cmin=org+V3{i*vs,j*vs,k*vs}, cmax=cmin+V3{vs,vs,vs};
// project crease line (x=0,y=0 in rotated frame): crease passes through the cell
// if the origin's (x,y) is within the cell's xy footprint (approx; crease is
// near x=y=0 line). Use: cell contains a point with small |rotated-xy|.
// Simplest: cell is a crease cell if it captures BOTH planes OR its center is
// within ~1 voxel of the crease line (x^2+y^2 small at center).
V3 cc=(cmin+cmax)*0.5f;
bool near_crease = (cc.x*cc.x+cc.y*cc.y) < (1.5f*vs)*(1.5f*vs);
if(!near_crease)continue;
crease_cells++;
// count distinct plane normals
int np=0; V3 seen[4];
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<4)seen[np++]=nm;}
plane_hist[np>4?4:np]++;
// vertex: if 2 planes → snap to crease; if 1 plane → stays on that face (mass point)
V3 mass{0,0,0}; for(auto&p:P)mass=mass+p; mass=mass*(1.0f/P.size());
V3 v;
if(np>=2){ // snap to crease line (x=0,y=0 rotated): project mass onto the two planes
// both planes through origin → crease is the z-axis in rotated frame.
// deviation after ideal crease-snap = 0. Represent as on-crease.
v = {0,0,mass.z};
} else { // 1 plane → vertex stays on the face at the mass point (OFF the crease)
v = mass;
}
float dev=sqrtf(v.x*v.x+v.y*v.y);
crease_dev.push_back(dev);
if(dev>0.05f)off_crease++;
}
// ===== BRICK APPROACH: at each near-crease cell, capture the brush faces within
// ~1 voxel of the CELL (from exact brush geometry, NOT grid crossings). A cell
// near the brush edge sees BOTH planes → snaps to the crease even if its grid
// crossings only hit one face. This is what a plane-set brick would bake. =====
std::vector<float> brick_dev; int brick_two=0, brick_one=0, brick_off=0, brick_cells=0;
float T = 1.2f * vs; // bake capture radius (cell within this of a face → store it)
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}, cc=cmin+V3{vs*0.5f,vs*0.5f,vs*0.5f};
bool near_crease = (cc.x*cc.x+cc.y*cc.y) < (1.5f*vs)*(1.5f*vs);
if(!near_crease)continue;
brick_cells++;
// capture brush faces within T of the cell centre (exact geometry)
bool haveA = fabsf(dot(p1.n,cc)-p1.d) < T;
bool haveB = fabsf(dot(p2.n,cc)-p2.d) < T;
V3 v;
if(haveA && haveB){ brick_two++; v={0,0,cc.z}; } // crease → on the crease line
else if(haveA){ brick_one++; v=cc-p1.n*(dot(p1.n,cc)-p1.d); } // face A
else { brick_one++; v=cc-p2.n*(dot(p2.n,cc)-p2.d); } // face B
float dev=sqrtf(v.x*v.x+v.y*v.y); brick_dev.push_back(dev);
if(dev>0.05f)brick_off++;
}
std::sort(brick_dev.begin(),brick_dev.end());
printf("=== BRICK (exact per-cell faces, T=%.2f) ===\n", T);
printf(" near-crease cells=%d: 2-plane(crease)=%d 1-plane(face)=%d OFF-crease=%d (%.0f%%)\n",
brick_cells, brick_two, brick_one, brick_off, brick_cells?100.0f*brick_off/brick_cells:0);
printf(" vertex dev-from-crease p95=%.4f MAX=%.4f\n\n",
brick_dev.empty()?0:brick_dev[(size_t)(brick_dev.size()*0.95)], brick_dev.empty()?0:brick_dev.back());
std::sort(crease_dev.begin(),crease_dev.end());
printf("=== GRID-CROSSING (current) ===\n");
printf("Crease (tilted %.2f rad, off-grid) — cells along the edge:\n", tilt);
printf(" distinct-plane count: 1-plane(face)=%d 2-plane(crease)=%d 3+=%d\n",
plane_hist[1], plane_hist[2], plane_hist[3]+plane_hist[4]);
printf(" crease cells=%d, OFF-crease(dev>0.05)=%d (%.0f%%)\n",
crease_cells, off_crease, crease_cells?100.0f*off_crease/crease_cells:0);
printf(" vertex dev-from-crease p50=%.4f p95=%.4f MAX=%.4f (voxel=%.3f)\n",
crease_dev.empty()?0:crease_dev[crease_dev.size()/2],
crease_dev.empty()?0:crease_dev[(size_t)(crease_dev.size()*0.95)],
crease_dev.empty()?0:crease_dev.back(), vs);
printf("\nIf many 1-plane cells sit ON the crease line → they place their vertex on\n");
printf("the FACE (off-crease) → sawtooth. That's feature-edge TOPOLOGY, which the\n");
printf("plane-set brick does NOT fix (needs finer res / manifold DC / crease detect).\n");
return 0;
}