// ============================================================================
// CREASE-SEGMENT ORACLE — grade DC vertex placement against the TRUE crease
// segment, and A/B the candidate fixes offline.
//
// WHY THIS EXISTS
// ---------------
// test_planeset.cpp measures "does the vertex land on the planes this cell
// captured" — ~0 by construction, because it grades the cell against its own
// decision. test_rotated_face.cpp measures flatness against the nearest FACE.
// Neither computes the true crease SEGMENT, so neither can see the defect the
// user reports: subtle stair-stepping along the straight edges of a ROTATED
// cube (the N-key diagnostic, hyper-tps map_dc_terrain.gd:326).
//
// The pipeline has no edge object. bake_planeset.glsl classifies each cell
// independently by distance from the CELL CENTRE to each brush plane
// (T = voxel_size * 1.2); qef_solve.glsl then snaps to the INFINITE line where
// two captured planes meet, and finally clamps the 3D point to the cell's
// axis-aligned box. Nothing knows where a crease starts or ends.
//
// THREE MODES, measured on identical cells:
// A CURRENT — the shipped pipeline, verbatim.
// B SURFACE NETS (naive) — vertex = centroid of the cell's edge crossings.
// This is qef_solve.glsl's own `mass_point` (:230/:249). No
// normals, no QEF, no snap. Always in-cell by construction,
// so the clamp can never fire — but no sharp features either.
// C SEGMENT-CLAMP — snap to the crease LINE, then constrain the coordinate
// ALONG the line to (line n cell) n (line n brush half-spaces)
// instead of clamping the 3D point to the box. Derived only
// from data qef_solve already has bound (gpu_edits binding 9 +
// convex_planes binding 10) — NOT from this harness's ground
// truth, so the number is what the shader could actually get.
//
// WHAT IS REPLAYED (byte-faithful ports, file:line of the original)
// density -f(p), f = max_i(n_i.p - d_i) [noise_sdf.glsl:177-193]
// edge active (d0>=0) != (d1>=0) [edge_classify.glsl:52]
// edge hermite 8-iter bisect + central diff [edge_intersect.glsl:196-226]
// exact hermite argmax face, |best| < vs*1.8 [edge_intersect.glsl:239-288]
// plane-set capture T = vs*1.2 from CELL CENTRE [bake_planeset.glsl:84-154]
// QEF eps=0.01, Cramer, clamp to cell [qef_solve.glsl:190-383]
// plane snap face / crease-LINE / corner [qef_solve.glsl:515-538]
// escape fallback nearest face if out of cell+0.5 [qef_solve.glsl:546-559]
// final clamp clamp(snapped, cmin, cmax) [qef_solve.glsl:560]
// The polish (qef_solve.glsl:390-474) is SKIPPED for hard cells in the real
// shader, so it is skipped here too.
//
// CONTROL: the same cube AXIS-ALIGNED. The harness is only trustworthy if it
// reproduces the asymmetry the user sees — rotated edges jagged, axis edges
// razor-sharp. If the control is not clean, the harness is wrong.
//
// Build: g++ -O2 -std=c++17 -o /tmp/crease_oracle tests/test_crease_segment_oracle.cpp && /tmp/crease_oracle
// ============================================================================
#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <vector>
#include <algorithm>
// ---------------------------------------------------------------------------
// Minimal vector math
// ---------------------------------------------------------------------------
struct V3 { float x, y, z; };
static inline V3 operator+(V3 a, V3 b) { return {a.x+b.x, a.y+b.y, a.z+b.z}; }
static inline V3 operator-(V3 a, V3 b) { return {a.x-b.x, a.y-b.y, a.z-b.z}; }
static inline V3 operator*(V3 a, float s) { return {a.x*s, a.y*s, a.z*s}; }
static inline float dot(V3 a, V3 b) { return a.x*b.x + a.y*b.y + a.z*b.z; }
static inline 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 inline float len3(V3 a) { return std::sqrt(dot(a, a)); }
static inline V3 norm3(V3 a) { float l = len3(a); return l > 1e-20f ? a * (1.0f/l) : a; }
static inline float clampf(float v, float lo, float hi) {
return v < lo ? lo : (v > hi ? hi : v);
}
struct Plane { V3 n; float d; }; // inside: n.p <= d
// ---------------------------------------------------------------------------
// Godot 4 Basis.from_euler(v) — default order EULER_ORDER_YXZ (ymat*xmat*zmat).
// map_dc_terrain.gd:340 uses Basis.from_euler(Vector3(0.35, 0.6, 0.0)); with
// z = 0 the zmat is identity, so B = Ry(0.6) * Rx(0.35).
// ---------------------------------------------------------------------------
static V3 basis_yxz(V3 v, float ex, float ey) {
float cx = std::cos(ex), sx = std::sin(ex);
float cy = std::cos(ey), sy = std::sin(ey);
return { cy*v.x + (sy*sx)*v.y + (sy*cx)*v.z,
cx*v.y + (-sx)*v.z,
-sy*v.x + (cy*sx)*v.y + (cy*cx)*v.z };
}
// The N-key brush (map_dc_terrain.gd:338-353)
static std::vector<Plane> make_cube(V3 center, float half, bool rotate) {
const V3 axes[6] = {{1,0,0},{-1,0,0},{0,1,0},{0,-1,0},{0,0,1},{0,0,-1}};
std::vector<Plane> out;
for (int i = 0; i < 6; i++) {
V3 n = norm3(rotate ? basis_yxz(axes[i], 0.35f, 0.6f) : axes[i]);
out.push_back({ n, dot(n, center) + half });
}
return out;
}
static float convex_f(V3 p, const std::vector<Plane>& pl) {
float m = -1e30f;
for (const auto& q : pl) m = std::max(m, dot(q.n, p) - q.d);
return m;
}
// ADD brush into deep air: d = max(d_terrain, -f) => -f. Only sign CHANGES are
// ever tested downstream, so the polarity convention affects nothing here.
static float density_at(V3 p, const std::vector<Plane>& pl) { return -convex_f(p, pl); }
// ---------------------------------------------------------------------------
// GROUND TRUTH: exact crease segments of the convex polytope.
// For each plane pair (i,j): line L = P_i ^ P_j, clipped by every other
// half-space. A cube yields 12.
// ---------------------------------------------------------------------------
struct Segment { V3 a, b; int pi, pj; };
static std::vector<Segment> true_crease_segments(const std::vector<Plane>& pl) {
std::vector<Segment> segs;
int n = (int)pl.size();
for (int i = 0; i < n; i++)
for (int j = i + 1; j < n; j++) {
V3 ni = pl[i].n, nj = pl[j].n;
V3 dir = cross3(ni, nj);
if (len3(dir) < 1e-6f) continue;
dir = norm3(dir);
float cc = dot(ni, nj), det = 1.0f - cc*cc;
if (det < 1e-8f) continue;
V3 p0 = ni * ((pl[i].d - cc*pl[j].d)/det) + nj * ((pl[j].d - cc*pl[i].d)/det);
float tmin = -1e30f, tmax = 1e30f; bool empty = false;
for (int k = 0; k < n && !empty; k++) {
if (k == i || k == j) continue;
float a = dot(pl[k].n, dir);
float b = pl[k].d - dot(pl[k].n, p0);
if (std::fabs(a) < 1e-9f) { if (b < -1e-6f) empty = true; }
else if (a > 0.0f) tmax = std::min(tmax, b/a);
else tmin = std::max(tmin, b/a);
}
if (empty || tmax - tmin <= 1e-5f) continue;
segs.push_back({ p0 + dir*tmin, p0 + dir*tmax, i, j });
}
return segs;
}
static float point_line_dist(V3 p, const Segment& s) {
V3 dir = norm3(s.b - s.a), w = p - s.a;
return len3(w - dir * dot(w, dir));
}
static float point_seg_dist(V3 p, const Segment& s) {
V3 ab = s.b - s.a; float L2 = dot(ab, ab);
if (L2 < 1e-20f) return len3(p - s.a);
return len3(p - (s.a + ab * clampf(dot(p - s.a, ab)/L2, 0.0f, 1.0f)));
}
static bool seg_hits_box(const Segment& s, V3 lo, V3 hi) {
V3 d = s.b - s.a;
float t0 = 0.0f, t1 = 1.0f;
const float* pa=&s.a.x; const float* pd=&d.x; const float* pl=&lo.x; const float* ph=&hi.x;
for (int a = 0; a < 3; a++) {
if (std::fabs(pd[a]) < 1e-12f) { if (pa[a] < pl[a] || pa[a] > ph[a]) return false; }
else {
float inv = 1.0f/pd[a], ta = (pl[a]-pa[a])*inv, tb = (ph[a]-pa[a])*inv;
if (ta > tb) std::swap(ta, tb);
t0 = std::max(t0, ta); t1 = std::min(t1, tb);
if (t0 > t1) return false;
}
}
return true;
}
// ---------------------------------------------------------------------------
// PIPELINE REPLAY
// ---------------------------------------------------------------------------
struct Hermite { V3 pos, n; bool hard; };
static Hermite edge_hermite(V3 p0, V3 p1, float d0, const std::vector<Plane>& pl,
V3 brush_c, float brush_r, float vs) {
float lo = 0.0f, hi = 1.0f, d_lo = d0;
for (int i = 0; i < 8; i++) { // [edge_intersect.glsl:196-208]
float mid = (lo + hi) * 0.5f;
float dm = density_at(p0 + (p1 - p0)*mid, pl);
if ((d_lo >= 0.0f) == (dm >= 0.0f)) { lo = mid; d_lo = dm; } else hi = mid;
}
float t = clampf((lo + hi)*0.5f, 0.001f, 0.999f); // [:211]
V3 pos = p0 + (p1 - p0)*t;
float eps = vs * 0.01f; // [:215-226]
V3 n = { density_at({pos.x+eps,pos.y,pos.z},pl) - density_at({pos.x-eps,pos.y,pos.z},pl),
density_at({pos.x,pos.y+eps,pos.z},pl) - density_at({pos.x,pos.y-eps,pos.z},pl),
density_at({pos.x,pos.y,pos.z+eps},pl) - density_at({pos.x,pos.y,pos.z-eps},pl) };
float L = len3(n);
n = (L > 1e-8f) ? n * (1.0f/L) : V3{0,1,0};
V3 cdn = n;
bool hit = false; Plane hit_pl{{0,0,0},0}; // [:239-288]
V3 dcv = pos - brush_c;
if (dot(dcv, dcv) <= brush_r*brush_r) {
float best = -1e30f; Plane best_pl{{0,0,0},0};
for (const auto& q : pl) { float dp = dot(q.n,pos)-q.d; if (dp > best) { best = dp; best_pl = q; } }
if (std::fabs(best) < vs*1.8f) { hit = true; hit_pl = best_pl; }
}
if (hit) {
V3 edir = p1 - p0;
float denom = dot(hit_pl.n, edir);
if (std::fabs(denom) > 1e-8f)
pos = p0 + edir * clampf((hit_pl.d - dot(hit_pl.n,p0))/denom, 0.0f, 1.0f);
n = (dot(hit_pl.n, cdn) < 0.0f) ? hit_pl.n * -1.0f : hit_pl.n;
}
return { pos, n, hit };
}
struct Capture { V3 n[3]; float d[3]; int count; int in_range; };
static Capture bake_capture(V3 cc, const std::vector<Plane>& pl,
V3 brush_c, float brush_r, float vs) {
Capture c; c.count = 0; c.in_range = 0;
float T = vs * 1.2f; // [bake_planeset.glsl:84]
V3 dcv = cc - brush_c; float rr = brush_r + T;
if (dot(dcv,dcv) > rr*rr) return c; // [:90-91]
V3 seen[8]; int ns = 0;
for (int pi = 0; pi < (int)pl.size(); pi++) {
if (std::fabs(dot(pl[pi].n,cc) - pl[pi].d) >= T) continue; // [:98]
bool onbrush = true; // [:99-104]
for (int pj = 0; pj < (int)pl.size(); pj++) {
if (pj == pi) continue;
if (dot(pl[pj].n,cc) - pl[pj].d > T) { onbrush = false; break; }
}
if (!onbrush) continue;
bool dup = false;
for (int q = 0; q < ns; q++) if (dot(pl[pi].n, seen[q]) > 0.999f) { dup = true; break; }
if (!dup && ns < 8) seen[ns++] = pl[pi].n;
int matched = -1; // [:131-135]
for (int k = 0; k < c.count; k++)
if (dot(pl[pi].n, c.n[k]) > 0.999f) { matched = k; break; }
if (matched < 0 && c.count < 3) { // [:136-140]
c.n[c.count] = pl[pi].n; c.d[c.count] = pl[pi].d; c.count++;
}
}
c.in_range = ns;
return c;
}
enum Mode { MODE_CURRENT = 0, MODE_SURFACE_NETS = 1, MODE_SEGMENT = 2,
MODE_SEG_CELLONLY = 3, MODE_SEG_CONTAIN = 4, MODE_SEG_FACEFALL = 5 };
#define NMODES 6
struct CellResult {
V3 vert;
V3 pre_clamp;
float clamp_dist = 0.0f;
int cl_count = 0;
bool hard = false;
bool fallback_fired= false;
bool clamp_moved = false;
bool det_degenerate= false;
bool overcap = false;
bool escaped_cell = false; // final vertex outside its own cell (safety)
bool seg_empty = false; // MODE_SEGMENT: constrained range was empty
};
// Constrain a point to the crease line, with the ALONG-line coordinate limited
// to (line n cell) n (line n the brush's own half-spaces).
// Everything here is available to qef_solve: the two captured planes come from
// the plane-set brick, the remaining half-spaces from gpu_edits/convex_planes.
static bool segment_clamp(V3 result, V3 n0, float d0, V3 n1, float d1,
const std::vector<Plane>& pl, V3 cmin, V3 cmax,
V3& out, bool use_brush_planes) {
float cc = dot(n0, n1), dd = 1.0f - cc*cc;
if (dd <= 1e-4f) return false; // near-parallel
float r0 = d0 - dot(n0, result), r1 = d1 - dot(n1, result);
V3 pline = result + n0*((r0 - cc*r1)/dd) + n1*((r1 - cc*r0)/dd); // closest pt on line
V3 dir = cross3(n0, n1);
if (len3(dir) < 1e-9f) return false;
dir = norm3(dir);
float tmin = -1e30f, tmax = 1e30f;
// clip to the cell box (slabs), t measured from pline
const float* pp=&pline.x; const float* pdv=&dir.x; const float* lo=&cmin.x; const float* hi=&cmax.x;
for (int a = 0; a < 3; a++) {
if (std::fabs(pdv[a]) < 1e-9f) { if (pp[a] < lo[a] || pp[a] > hi[a]) return false; }
else {
float ta = (lo[a]-pp[a])/pdv[a], tb = (hi[a]-pp[a])/pdv[a];
if (ta > tb) std::swap(ta, tb);
tmin = std::max(tmin, ta); tmax = std::min(tmax, tb);
}
}
// clip to the brush's other half-spaces => stay on the REAL finite edge
if (use_brush_planes) for (const auto& q : pl) {
if (dot(q.n, n0) > 0.999f || dot(q.n, n1) > 0.999f) continue;
float a = dot(q.n, dir), b = q.d - dot(q.n, pline);
if (std::fabs(a) < 1e-9f) { if (b < -1e-6f) return false; }
else if (a > 0.0f) tmax = std::min(tmax, b/a);
else tmin = std::max(tmin, b/a);
}
if (tmin > tmax) return false;
out = pline + dir * clampf(0.0f, tmin, tmax);
return true;
}
static CellResult solve_cell(const std::vector<Hermite>& H, V3 cmin, float vs,
const Capture& cap, const std::vector<Plane>& pl,
float eps, Mode mode) {
CellResult R;
V3 cmax = cmin + V3{vs, vs, vs};
float ata[6] = {0,0,0,0,0,0};
V3 atb{0,0,0}, mass{0,0,0};
int cnt = 0;
for (const auto& h : H) { // [qef_solve.glsl:213-241]
V3 nm = h.n;
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;
atb = atb + nm * dot(nm, h.pos);
mass = mass + h.pos;
cnt++;
if (h.hard) R.hard = true;
}
if (cnt == 0) { R.vert = R.pre_clamp = cmin + V3{vs*.5f,vs*.5f,vs*.5f}; return R; }
mass = mass * (1.0f / cnt);
// ---- MODE B: naive Surface Nets = the centroid, i.e. qef_solve's own
// mass_point (:230/:249). No normals, no QEF, no snap, no clamp needed.
if (mode == MODE_SURFACE_NETS) {
R.cl_count = cap.count;
R.vert = R.pre_clamp = mass;
return R;
}
float m00=ata[0]+eps, m01=ata[1], m02=ata[2]; // [:316-378]
float 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 result = mass;
if (std::fabs(det) > 1e-12f) {
float iv = 1.0f/det;
result = 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 };
}
result = { clampf(result.x,cmin.x,cmax.x), // [:383]
clampf(result.y,cmin.y,cmax.y),
clampf(result.z,cmin.z,cmax.z) };
R.cl_count = cap.count;
if (cap.count >= 1) {
// ---- MODE C: constrain ALONG the crease instead of clamping the point.
if ((mode == MODE_SEGMENT || mode == MODE_SEG_CELLONLY ||
mode == MODE_SEG_CONTAIN || mode == MODE_SEG_FACEFALL) && cap.count >= 2) {
int a = 0, b = 1;
if (cap.count == 3) {
// A cell near, but not at, a corner is really on an edge: pick the
// crease pair whose line passes closest to the cell centre.
V3 ctr = cmin + V3{vs*.5f, vs*.5f, vs*.5f};
float best = 1e30f;
const int PR[3][2] = {{0,1},{0,2},{1,2}};
for (int k = 0; k < 3; k++) {
V3 na = cap.n[PR[k][0]], nb = cap.n[PR[k][1]];
float cc2 = dot(na,nb), dd2 = 1.0f - cc2*cc2;
if (dd2 <= 1e-4f) continue;
float ra = cap.d[PR[k][0]] - dot(na,ctr), rb = cap.d[PR[k][1]] - dot(nb,ctr);
V3 p = ctr + na*((ra-cc2*rb)/dd2) + nb*((rb-cc2*ra)/dd2);
float dist = len3(p - ctr);
if (dist < best) { best = dist; a = PR[k][0]; b = PR[k][1]; }
}
// If the true 3-plane corner is inside this cell, keep the corner.
V3 A=cap.n[0], B=cap.n[1], C=cap.n[2];
float den = dot(A, cross3(B, C));
if (std::fabs(den) > 1e-5f) {
V3 corner = (cross3(B,C)*cap.d[0] + cross3(C,A)*cap.d[1]
+ cross3(A,B)*cap.d[2]) * (1.0f/den);
if (corner.x >= cmin.x && corner.x <= cmax.x &&
corner.y >= cmin.y && corner.y <= cmax.y &&
corner.z >= cmin.z && corner.z <= cmax.z) {
R.vert = R.pre_clamp = corner;
return R;
}
}
}
V3 out;
if (segment_clamp(result, cap.n[a], cap.d[a], cap.n[b], cap.d[b],
pl, cmin, cmax, out, mode == MODE_SEGMENT)) {
R.vert = R.pre_clamp = out;
return R;
}
R.seg_empty = true;
// ---- MODE F: an empty range PROVES the crease line does not cross this
// cell, i.e. the T=1.2*voxel centre test tagged a FALSE crease. Falling
// back to the crease snap (mode C) is what pushes the vertex outside the
// solid: measured, 139/139 of the remaining protruding cells are exactly
// this population. Treat the cell as what it actually is — a FACE cell —
// and snap to the nearest captured plane. A point on a face plane inside
// the cell is ON the surface, so it cannot protrude.
if (mode == MODE_SEG_FACEFALL) {
V3 ctr = cmin + V3{vs*.5f, vs*.5f, vs*.5f};
int nb = 0; float ba = 1e30f;
for (int c = 0; c < cap.count; c++) {
float aa = std::fabs(dot(cap.n[c], ctr) - cap.d[c]);
if (aa < ba) { ba = aa; nb = c; }
}
V3 fs = result - cap.n[nb] * (dot(cap.n[nb], result) - cap.d[nb]);
R.vert = R.pre_clamp = { clampf(fs.x, cmin.x, cmax.x),
clampf(fs.y, cmin.y, cmax.y),
clampf(fs.z, cmin.z, cmax.z) };
return R;
}
// otherwise fall through to CURRENT behaviour
}
// ---- MODE A (and MODE C fallback): the shipped snap. [:515-561]
V3 snapped = result;
if (cap.count == 1) {
snapped = result - cap.n[0]*(dot(cap.n[0],result) - cap.d[0]);
} else if (cap.count == 2) {
float cc2 = dot(cap.n[0], cap.n[1]), dd2 = 1.0f - cc2*cc2;
if (dd2 > 1e-4f) {
float r0 = cap.d[0]-dot(cap.n[0],result), r1 = cap.d[1]-dot(cap.n[1],result);
snapped = result + cap.n[0]*((r0-cc2*r1)/dd2) + cap.n[1]*((r1-cc2*r0)/dd2);
} else R.det_degenerate = true;
} else {
V3 A=cap.n[0], B=cap.n[1], C=cap.n[2];
float den = dot(A, cross3(B,C));
if (std::fabs(den) > 1e-5f)
snapped = (cross3(B,C)*cap.d[0]+cross3(C,A)*cap.d[1]+cross3(A,B)*cap.d[2])*(1.0f/den);
}
if (cap.count >= 2) { // [:546-559]
V3 lo = cmin - V3{vs*.5f,vs*.5f,vs*.5f}, hi = cmax + V3{vs*.5f,vs*.5f,vs*.5f};
if (snapped.x<lo.x||snapped.y<lo.y||snapped.z<lo.z||
snapped.x>hi.x||snapped.y>hi.y||snapped.z>hi.z) {
R.fallback_fired = true;
V3 ctr = cmin + V3{vs*.5f,vs*.5f,vs*.5f};
int nb = 0; float ba = 1e30f;
for (int c = 0; c < cap.count; c++) {
float aa = std::fabs(dot(cap.n[c],ctr) - cap.d[c]);
if (aa < ba) { ba = aa; nb = c; }
}
snapped = result - cap.n[nb]*(dot(cap.n[nb],result) - cap.d[nb]);
}
}
R.pre_clamp = snapped;
V3 clamped = { clampf(snapped.x,cmin.x,cmax.x), // [:560]
clampf(snapped.y,cmin.y,cmax.y),
clampf(snapped.z,cmin.z,cmax.z) };
R.clamp_dist = len3(clamped - snapped);
if (R.clamp_dist > 1e-6f) R.clamp_moved = true;
result = clamped;
} else {
R.pre_clamp = result;
}
R.vert = result;
// ---- MODE E: CONTAINMENT. Nothing anywhere in the snap path constrains the
// vertex to stay INSIDE the union. Surface Nets is immune (a centroid of
// on-cell-edge crossings is necessarily inside the local solid), which is why
// B measures 0 protrusion at every voxel size while A/C/D do not. Enforce that
// property explicitly: if the snapped vertex is outside, project it onto the
// violated (argmax) face, then re-clamp to the cell.
if (mode == MODE_SEG_CONTAIN) {
float f = convex_f(R.vert, pl);
if (f > 1e-6f) {
float best = -1e30f; V3 bn{0,0,0}; float bd = 0.0f;
for (const auto& q : pl) {
float dp = dot(q.n, R.vert) - q.d;
if (dp > best) { best = dp; bn = q.n; bd = q.d; }
}
V3 proj = R.vert - bn * (dot(bn, R.vert) - bd);
R.vert = { clampf(proj.x, cmin.x, cmax.x),
clampf(proj.y, cmin.y, cmax.y),
clampf(proj.z, cmin.z, cmax.z) };
}
}
return R;
}
// ---------------------------------------------------------------------------
struct Stats { float p50=0, p95=0, mx=0; int n=0; };
static Stats stat_of(std::vector<float> v) {
Stats s; s.n = (int)v.size();
if (v.empty()) return s;
std::sort(v.begin(), v.end());
s.p50 = v[v.size()/2];
size_t i95 = (size_t)(v.size()*0.95f); if (i95 >= v.size()) i95 = v.size()-1;
s.p95 = v[i95]; s.mx = v.back();
return s;
}
struct ModeReport {
Stats dev; // deviation from the true crease LINE (crease cells)
Stats dev_corner; // ...restricted to cl==3 cells
Stats surf_false; // false-crease cells: distance from the true surface
int escapes = 0; // vertices outside their own cell (containment safety)
int seg_empty = 0;
};
static void run(const char* label, bool rotate, float vs, float half) {
V3 center{0,0,0};
auto pl = make_cube(center, half, rotate);
float brush_r = half * 1.9f; // map_dc_terrain.gd:353
auto segs = true_crease_segments(pl);
float margin = 2.0f;
V3 org{ -half-margin, -half-margin, -half-margin };
int G = (int)std::ceil((2*half + 2*margin)/vs) + 2;
auto DP = [&](int i,int j,int k){ return density_at(org + V3{i*vs,j*vs,k*vs}, pl); };
static const int EOFF[12][4] = { // [qef_solve.glsl:172-188]
{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} };
static const int EDIR[3][3] = {{1,0,0},{0,1,0},{0,0,1}};
ModeReport rep[NMODES];
std::vector<float> dev[NMODES], devc[NMODES], surf[NMODES];
// SIGNED OUTWARD PROTRUSION — the metric that corresponds to a visible
// silhouette notch. A vertex can sit exactly on the crease LINE and still
// bulge out of the solid; deviation-from-line cannot see that. Measured on
// EVERY active cell, not just crease cells, because a face cell next to the
// crease protrudes into the same silhouette.
std::vector<float> prot[NMODES]; // max(0, f(vert)) over ALL active cells
std::vector<float> prot_crease[NMODES]; // ...restricted to crease cells
int prot_cells[NMODES] = {0,0,0,0,0,0}; // count with protrusion > 0.1 voxel
int active=0, on_seg=0, cl_hist[4]={0,0,0,0}, off_seg=0, false_cl2=0;
int lbl_cl[4]={0,0,0,0}; int lbl_crease=0,lbl_false=0,lbl_other=0,lbl_fallback=0,lbl_nosnap=0;
int fallback_on=0, clamp_on=0, degen_on=0, overcap_on=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++) {
std::vector<Hermite> H; H.reserve(12);
for (int e = 0; e < 12; e++) {
int ex=i+EOFF[e][0], ey=j+EOFF[e][1], ez=k+EOFF[e][2], dir=EOFF[e][3];
if (ex<0||ey<0||ez<0||ex>=G-1||ey>=G-1||ez>=G-1) continue;
int nx=ex+EDIR[dir][0], ny=ey+EDIR[dir][1], nz=ez+EDIR[dir][2];
float d0=DP(ex,ey,ez), d1=DP(nx,ny,nz);
if ((d0>=0.0f) == (d1>=0.0f)) continue; // [edge_classify.glsl:52]
H.push_back(edge_hermite(org+V3{ex*vs,ey*vs,ez*vs},
org+V3{nx*vs,ny*vs,nz*vs}, d0, pl, center, brush_r, vs));
}
if (H.empty()) continue;
active++;
V3 cmin = org + V3{i*vs,j*vs,k*vs};
V3 cmax = cmin + V3{vs,vs,vs};
V3 cctr = cmin + V3{vs*.5f,vs*.5f,vs*.5f};
Capture cap = bake_capture(cctr, pl, center, brush_r, vs);
int best_seg = -1; float best_d = 1e30f;
for (int s = 0; s < (int)segs.size(); s++) {
if (!seg_hits_box(segs[s], cmin, cmax)) continue;
float dd = point_seg_dist(cctr, segs[s]);
if (dd < best_d) { best_d = dd; best_seg = s; }
}
if (best_seg >= 0) { on_seg++; cl_hist[std::min(cap.count,3)]++; }
else off_seg++;
for (int m = 0; m < NMODES; m++) {
CellResult R = solve_cell(H, cmin, vs, cap, pl, 0.01f, (Mode)m);
const float E = 1e-4f;
if (R.vert.x < cmin.x-E || R.vert.x > cmax.x+E ||
R.vert.y < cmin.y-E || R.vert.y > cmax.y+E ||
R.vert.z < cmin.z-E || R.vert.z > cmax.z+E) rep[m].escapes++;
if (R.seg_empty) rep[m].seg_empty++;
// outward protrusion: convex_f > 0 means the vertex is OUTSIDE the solid
float f_out = convex_f(R.vert, pl);
float pr = f_out > 0.0f ? f_out : 0.0f;
prot[m].push_back(pr);
if (pr > vs*0.1f) {
prot_cells[m]++;
if (m == MODE_SEGMENT) { // WHICH cells still protrude under the shipped fix?
lbl_cl[std::min(cap.count,3)]++;
if (best_seg >= 0) lbl_crease++;
else if (cap.count == 2) lbl_false++;
else lbl_other++;
if (R.fallback_fired) lbl_fallback++;
if (R.seg_empty) lbl_nosnap++; // segment clamp gave up -> fell back to A
}
}
if (best_seg >= 0) prot_crease[m].push_back(pr);
if (best_seg >= 0) {
float d = point_line_dist(R.vert, segs[best_seg]);
dev[m].push_back(d);
if (cap.count == 3) devc[m].push_back(d);
if (m == 0) {
if (R.fallback_fired) fallback_on++;
if (R.clamp_moved) clamp_on++;
if (R.det_degenerate) degen_on++;
if (cap.in_range > 3) overcap_on++;
}
} else if (cap.count == 2) {
if (m == 0) false_cl2++;
float fd = 1e30f;
for (const auto& q : pl) fd = std::min(fd, std::fabs(dot(q.n,R.vert)-q.d));
surf[m].push_back(fd);
}
}
}
for (int m = 0; m < NMODES; m++) {
rep[m].dev = stat_of(dev[m]);
rep[m].dev_corner = stat_of(devc[m]);
rep[m].surf_false = stat_of(surf[m]);
}
printf("\n================ %s ================\n", label);
printf("voxel=%.3f half=%.1f true segments=%d active cells=%d\n", vs, half, (int)segs.size(), active);
printf(" crease cells=%d captured cl: 0=%d 1(MISSED)=%d 2=%d 3=%d (missed rate %.1f%%)\n",
on_seg, cl_hist[0], cl_hist[1], cl_hist[2], cl_hist[3],
on_seg ? 100.0f*(cl_hist[0]+cl_hist[1])/on_seg : 0.0f);
printf(" [mode A] fallback=%d clamp-displaced=%d det-degen=%d >3planes=%d\n",
fallback_on, clamp_on, degen_on, overcap_on);
const char* names[NMODES] = { "A CURRENT ", "B SURFACE NETS ", "C SEG+BRUSHCLIP ", "D SEG CELL-ONLY ", "E SEG+CONTAIN ", "F SEG+FACEFALL " };
printf("\n DEVIATION FROM THE TRUE CREASE LINE (n=%d crease cells)\n", on_seg);
printf(" %-16s %-10s %-10s %-12s %-9s %s\n", "mode", "p50", "p95", "MAX", "MAX/voxel", "escapes cell");
for (int m = 0; m < NMODES; m++)
printf(" %-16s %-10.5f %-10.5f %-12.5f %-9.2f %d\n",
names[m], rep[m].dev.p50, rep[m].dev.p95, rep[m].dev.mx,
rep[m].dev.mx/vs, rep[m].escapes);
printf("\n corner cells only (cl==3, n=%d):\n", rep[0].dev_corner.n);
for (int m = 0; m < NMODES; m++)
printf(" %-16s p95=%-10.5f MAX=%.5f\n", names[m], rep[m].dev_corner.p95, rep[m].dev_corner.mx);
printf("\n OUTWARD PROTRUSION — vertex OUTSIDE the true solid (visible silhouette notch)\n");
printf(" all active cells n=%d:\n", (int)prot[0].size());
for (int m = 0; m < NMODES; m++) {
Stats ps = stat_of(prot[m]);
printf(" %-16s p50=%-9.5f p95=%-9.5f MAX=%-9.5f (%.2f vox MAX) cells>0.1vox=%d\n",
names[m], ps.p50, ps.p95, ps.mx, ps.mx/vs, prot_cells[m]);
}
printf(" WHICH cells still protrude under mode C (>0.1 voxel):\n");
printf(" by cl_count: 0=%d 1(face)=%d 2(crease)=%d 3(corner)=%d\n",
lbl_cl[0], lbl_cl[1], lbl_cl[2], lbl_cl[3]);
printf(" true-crease=%d false-crease=%d other=%d | escape-fallback=%d snap-unresolved=%d\n",
lbl_crease, lbl_false, lbl_other, lbl_fallback, lbl_nosnap);
printf(" crease cells only:\n");
for (int m = 0; m < NMODES; m++) {
Stats ps = stat_of(prot_crease[m]);
printf(" %-16s p50=%-9.5f p95=%-9.5f MAX=%.5f\n", names[m], ps.p50, ps.p95, ps.mx);
}
printf("\n FALSE CREASE (cl==2, no true crease in cell): %d cells of %d off-crease\n", false_cl2, off_seg);
printf(" distance from the true surface:\n");
for (int m = 0; m < NMODES; m++)
printf(" %-16s p50=%-9.5f p95=%-9.5f MAX=%.5f\n",
names[m], rep[m].surf_false.p50, rep[m].surf_false.p95, rep[m].surf_false.mx);
printf("\n range-empty -> fell back to A: C=%d cells D=%d cells\n", rep[2].seg_empty, rep[3].seg_empty);
}
int main(int argc, char** argv) {
float vs = 0.224f; // map_dc_terrain.gd:86 — the deployed value
if (argc > 1) vs = (float)atof(argv[1]);
printf("CREASE-SEGMENT ORACLE — DC vertex placement vs the TRUE crease segment.\n");
printf("A=shipped B=naive Surface Nets (centroid) C=segment-clamp candidate.\n");
printf("Perfect: deviation 0, escapes 0. The AXIS control must be clean or the\n");
printf("harness is wrong, not the pipeline.\n");
run("CONTROL — AXIS cube, half=5", false, vs, 5.0f);
run("SUBJECT — ROTATED cube, half=5", true, vs, 5.0f);
run("CONTROL — AXIS cube, half=2", false, vs, 2.0f);
run("SUBJECT — ROTATED cube, half=2", true, vs, 2.0f);
return 0;
}