// CPU oracle for STEP 1: exact analytic Hermite -> QEF -> exact vertex,
// WITHOUT vertex-polish and WITHOUT plane-snap.
//
// Hypothesis (Cepero / Ju 2002): if the edge crossings + normals fed to the QEF
// are EXACT (analytic ray-vs-plane, not sampled off a blended density field),
// then Dual Contouring's QEF ALONE places the cell vertex exactly on the
// face / crease / corner — no downstream snapping needed. This mirrors the
// EXACT math in qef_solve.glsl Step 3 (ATA/atb, eps=0.01 regularise toward mass
// point, Cramer's rule) but skips polish + snap. If this passes, the shader
// change is: (a) make edge_intersect positions analytic, (b) skip polish/snap
// for hard cells -> razor cubes through the clean path.
#include <cstdio>
#include <cmath>
#include <vector>
struct V3 { float x, y, z; };
static V3 sub(V3 a, V3 b){ return {a.x-b.x, a.y-b.y, a.z-b.z}; }
static V3 add(V3 a, V3 b){ return {a.x+b.x, a.y+b.y, a.z+b.z}; }
static V3 mul(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; }; // n·x = d (n unit)
struct Edge { V3 p0, p1; }; // a cell edge (two grid endpoints)
// ---- EXACT analytic crossing: ray p0->p1 vs plane n·x=d (edge_intersect Step1)
// t = (d - n·p0) / (n·(p1-p0)); pos = p0 + t*(p1-p0). Returns false if parallel.
static bool exact_cross(Edge e, Plane pl, V3& pos, V3& nrm) {
V3 dir = sub(e.p1, e.p0);
float denom = dot(pl.n, dir);
if (std::fabs(denom) < 1e-8f) return false; // parallel: no clean crossing
float t = (pl.d - dot(pl.n, e.p0)) / denom;
if (t < -0.001f || t > 1.001f) return false; // crossing not on this edge
pos = add(e.p0, mul(dir, t));
nrm = pl.n;
return true;
}
// ---- QEF exactly as qef_solve.glsl Step 3 (NO polish, NO snap) --------------
struct Hermite { V3 pos, nrm; };
static V3 qef_solve(const std::vector<Hermite>& hs, V3 cell_min, V3 cell_max) {
float ata[6] = {0,0,0,0,0,0};
V3 atb = {0,0,0};
V3 mass = {0,0,0};
int n = 0;
for (auto& h : hs) {
V3 nm = h.nrm;
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, h.pos);
atb.x+=d*nm.x; atb.y+=d*nm.y; atb.z+=d*nm.z;
mass = add(mass, h.pos); n++;
}
if (n == 0) return mass;
mass = mul(mass, 1.0f/float(n));
float eps = 0.01f;
float m00=ata[0]+eps, m01=ata[1], m02=ata[2];
float m11=ata[3]+eps, m12=ata[4], m22=ata[5]+eps;
// rhs = atb - ATA*mass (unregularised ATA)
V3 atamp = {
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 = sub(atb, atamp);
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 result;
if (std::fabs(det) > 1e-12f) {
float inv = 1.0f/det;
V3 off = {
(c00*rhs.x + c01*rhs.y + c02*rhs.z)*inv,
(c01*rhs.x + c11*rhs.y + c12*rhs.z)*inv,
(c02*rhs.x + c12*rhs.y + c22*rhs.z)*inv
};
result = add(mass, off);
} else result = mass;
// clamp to cell
result.x = std::fmax(cell_min.x, std::fmin(cell_max.x, result.x));
result.y = std::fmax(cell_min.y, std::fmin(cell_max.y, result.y));
result.z = std::fmax(cell_min.z, std::fmin(cell_max.z, result.z));
return result;
}
// Build a cell's Hermite by intersecting its 12 edges with the cube's planes,
// selecting per edge the ACTIVE plane (max n·p-d among the planes the edge
// actually crosses) — mirrors edge_intersect's argmax face pick.
static std::vector<Hermite> cell_hermite(V3 cmin, V3 cmax, const std::vector<Plane>& planes) {
// 8 corners
V3 c[8]; int k=0;
for(int zi=0;zi<2;zi++)for(int yi=0;yi<2;yi++)for(int xi=0;xi<2;xi++)
c[k++] = { xi?cmax.x:cmin.x, yi?cmax.y:cmin.y, zi?cmax.z:cmin.z };
// 12 edges as corner-index pairs
int E[12][2] = {{0,1},{2,3},{4,5},{6,7},{0,2},{1,3},{4,6},{5,7},{0,4},{1,5},{2,6},{3,7}};
// A corner is "solid" if inside ALL planes (n·x <= d for every plane) => convex brush.
auto inside=[&](V3 p){ for(auto&pl:planes) if(dot(pl.n,p)-pl.d > 1e-6f) return false; return true; };
std::vector<Hermite> out;
for (int e=0;e<12;e++){
V3 a=c[E[e][0]], b=c[E[e][1]];
bool ia=inside(a), ib=inside(b);
if (ia==ib) continue; // no sign change on this edge
Edge ed{a,b};
// pick the plane this edge crosses with the crossing on-segment; among
// valid, the ACTIVE face is the one the solid endpoint is ON/closest to.
float best_t_key = 1e30f; V3 bpos, bnrm; bool found=false;
for (auto& pl : planes) {
V3 pos,nrm;
if (!exact_cross(ed, pl, pos, nrm)) continue;
// must be a real boundary: the crossing point lies on the brush surface
// (inside all OTHER planes). Reject crossings outside the convex body.
bool onsurf=true;
for (auto& q : planes){ if(&q==&pl) continue; if(dot(q.n,pos)-q.d > 1e-4f){onsurf=false;break;} }
if(!onsurf) continue;
// among on-surface crossings keep the first (there should be exactly one)
float key = 0.0f;
if (key < best_t_key){ best_t_key=key; bpos=pos; bnrm=nrm; found=true; }
}
if(found) out.push_back({bpos,bnrm});
}
return out;
}
static float plane_err(V3 x, Plane p){ return std::fabs(dot(p.n,x)-p.d); }
int main(){
int fails=0;
// Axis-aligned cube, faces at ±1. Planes: n·x <= d inside.
std::vector<Plane> cube = {
{{ 1,0,0}, 1},{{-1,0,0}, 1},
{{0, 1,0}, 1},{{0,-1,0}, 1},
{{0,0, 1}, 1},{{0,0,-1}, 1},
};
Plane px{{1,0,0},1}, py{{0,1,0},1}, pz{{0,0,1},1};
// voxel_size 0.4; cube much bigger than a voxel (isolated-cube regime).
// FACE cell straddling x=1 (well away from edges/corners):
{
V3 cmin={0.8f,0.0f,0.0f}, cmax={1.2f,0.4f,0.4f};
auto h = cell_hermite(cmin,cmax,cube);
V3 r = qef_solve(h,cmin,cmax);
float er = plane_err(r,px);
printf("FACE n=%zu vert=(%.4f,%.4f,%.4f) x-plane_err=%.5f %s\n",
h.size(), r.x,r.y,r.z, er, er<1e-3f?"OK":"FAIL");
if(er>=1e-3f) fails++;
}
// EDGE cell straddling x=1,y=1 crease:
{
V3 cmin={0.8f,0.8f,0.0f}, cmax={1.2f,1.2f,0.4f};
auto h = cell_hermite(cmin,cmax,cube);
V3 r = qef_solve(h,cmin,cmax);
float er = plane_err(r,px)+plane_err(r,py);
printf("EDGE n=%zu vert=(%.4f,%.4f,%.4f) crease_err=%.5f %s\n",
h.size(), r.x,r.y,r.z, er, er<1e-3f?"OK":"FAIL");
if(er>=1e-3f) fails++;
}
// CORNER cell at x=1,y=1,z=1:
{
V3 cmin={0.8f,0.8f,0.8f}, cmax={1.2f,1.2f,1.2f};
auto h = cell_hermite(cmin,cmax,cube);
V3 r = qef_solve(h,cmin,cmax);
float er = plane_err(r,px)+plane_err(r,py)+plane_err(r,pz);
printf("CORNER n=%zu vert=(%.4f,%.4f,%.4f) expect(1,1,1) err=%.5f %s\n",
h.size(), r.x,r.y,r.z, er, er<1e-3f?"OK":"FAIL");
if(er>=1e-3f) fails++;
}
printf("\n%s (%d failures)\n", fails==0?"ALL PASS":"FAILURES", fails);
return fails?1:0;
}