#version 450
// ============================================================================
// QEF solve compute shader for Dual Contouring terrain (Stage 8).
//
// For each active voxel, gathers hermite data (intersection positions and
// normals) from its 12 cube edges, solves the Quadratic Error Function to
// find the optimal vertex position, and writes the result to the compacted
// vertex buffer.
//
// The QEF minimises E(x) = sum_i (n_i . (x - p_i))^2
// by accumulating the normal equations A^T A x = A^T b and solving via
// Cramer's rule on the regularised system (A^T A + eps*I).
// Falls back to the mass point when the determinant is near zero.
//
// Dispatch: ceil(32768 / 256), 1, 1
// ============================================================================
#define THREADS 256
layout(local_size_x = THREADS, local_size_y = 1, local_size_z = 1) in;
// ---- Buffers ----
layout(set = 0, binding = 0, std430) readonly buffer DensityBuf { float density[]; };
layout(set = 0, binding = 1, std430) readonly buffer EdgeFlagBuf { uint edge_flags[]; };
layout(set = 0, binding = 2, std430) readonly buffer EdgeScanBuf { uint edge_scan[]; };
layout(set = 0, binding = 3, std430) readonly buffer EdgeDataBuf { float edge_data[]; };
layout(set = 0, binding = 4, std430) readonly buffer VoxelFlagBuf { uint voxel_flags[]; };
layout(set = 0, binding = 5, std430) readonly buffer VoxelScanBuf { uint voxel_scan[]; };
layout(set = 0, binding = 6, std430) writeonly buffer VertexBuf { float vertices[]; };
// Phase E: Compacted active voxel index list + counts
layout(set = 0, binding = 7, std430) readonly buffer ActiveVoxelsBuf { uint active_voxel_indices[]; };
layout(set = 0, binding = 8, std430) readonly buffer CountBuf { uint counts[]; };
// Convex-brush edits (binding 9) + planes (binding 10) — the SAME RIDs bound to
// edge_intersect. Used for the feature-aware per-cell PLANE-SET snap: capture the
// exact brush faces near a hard cell (face/crease/corner) even when the grid
// crossings only sampled one face — fixes the angled-edge sawtooth.
struct GpuEdit {
uint edit_type;
float pos_x, pos_y, pos_z;
float half_x, half_y, half_z;
float radius;
float strength;
uint _pad0, _pad1, _pad2; // _pad0=material_id, _pad1=plane offset, _pad2=count
};
layout(set = 0, binding = 9, std430) readonly buffer EditBuf {
uint edit_count;
uint _edit_pad0, _edit_pad1, _edit_pad2;
GpuEdit gpu_edits[];
};
layout(set = 0, binding = 10, std430) readonly buffer ConvexPlaneBuf { vec4 convex_planes[]; };
// Plane-set BRICK (baked per cell by bake_planeset.glsl this pass). qef READS it
// instead of re-deriving the capture — O(voxels) not O(edits), and the mutability
// substrate: settled architecture is baked, remeshes read it.
struct CellPlaneSet {
uint cl_count;
uint material;
uint proj0; // two candidate face projection indices (Stage-2); the fragment picks
uint proj1; // between them by normal so sharp creases don't stretch one across.
vec4 planes[3];
};
layout(set = 0, binding = 11, std430) readonly buffer PlaneSetBuf { CellPlaneSet cell_planeset[]; };
// ---- Push constants ----
layout(push_constant) uniform PushConstants {
vec3 cell_origin; // 0: world-space origin of this chunk
float voxel_size; // 12: size of one voxel in world units
ivec3 grid_dims; // 16: voxel grid dimensions (32, 32, 32)
uint total_elements; // 28: grid_dims.x * grid_dims.y * grid_dims.z (32768)
uint mesh_quality_flags;
uint _pc_pad0;
uint _pc_pad1;
uint _pc_pad2; // 32: bit 1 = adaptive_qef
};
// ---- Constants ----
// Tikhonov regularization (eps*I added to ATA diagonal).
// NOT the same as Gildea's pseudoinverse threshold (0.1) — Tikhonov shifts
// all eigenvalues, while pseudoinverse zeros small ones. 0.01 is the right
// weight for Cramer's rule: it stabilizes degenerate directions without
// over-smoothing sharp CSG features.
#define QEF_REGULARISATION 0.01
#define DET_THRESHOLD 1e-12
// ============================================================================
// Helpers
// ============================================================================
// Linear index into the V^3 voxel grid.
int voxel_index(int x, int y, int z, int V) {
return z * V * V + y * V + x;
}
// Linear index into the (V+1)^3 density grid.
int density_index(int x, int y, int z, int S) {
return z * S * S + y * S + x;
}
// Get edge flag for edge at (ex, ey, ez) in direction dir (0=X, 1=Y, 2=Z).
// Returns 0 if out of bounds.
uint get_edge_flag(int ex, int ey, int ez, int dir, int V) {
if (ex < 0 || ey < 0 || ez < 0 || ex >= V || ey >= V || ez >= V)
return 0u;
return edge_flags[uint(voxel_index(ex, ey, ez, V)) * 3u + uint(dir)];
}
// Get edge scan (compacted index) for edge at (ex, ey, ez) in direction dir.
// Caller must ensure coordinates are in bounds.
uint get_edge_scan(int ex, int ey, int ez, int dir, int V) {
return edge_scan[uint(voxel_index(ex, ey, ez, V)) * 3u + uint(dir)];
}
// ============================================================================
// Main
// ============================================================================
void main() {
uint tid = gl_GlobalInvocationID.x;
if (tid >= counts[2]) // active voxel count from prepare_dispatch (counts[1] = quad_emit index_count)
return;
uint gid = active_voxel_indices[tid]; // original voxel grid index
int V = grid_dims.x; // 32
int S = V + 1; // 33 (density grid stride)
// Decode linear voxel index to (x, y, z)
int z = int(gid) / (V * V);
int y = (int(gid) % (V * V)) / V;
int x = int(gid) % V;
// ---- Step 1: Compute 8 corner signs ----
// Corners are at (x..x+1, y..y+1, z..z+1) in the density grid.
uint corner_signs = 0u;
for (int cz = 0; cz <= 1; cz++) {
for (int cy = 0; cy <= 1; cy++) {
for (int cx = 0; cx <= 1; cx++) {
float d = density[density_index(x + cx, y + cy, z + cz, S)];
uint bit = uint(cz * 4 + cy * 2 + cx);
if (d >= 0.0)
corner_signs |= (1u << bit);
}
}
}
// ---- Step 2: Gather hermite data from 12 edges ----
// Accumulate A^T A (symmetric 3x3, 6 values) and A^T b (3 values).
//
// Symmetric matrix layout (matching CPU dc_qef.h):
// ata[0]=xx ata[1]=xy ata[2]=xz ata[3]=yy ata[4]=yz ata[5]=zz
float ata[6] = float[6](0.0, 0.0, 0.0, 0.0, 0.0, 0.0);
float atb[3] = float[3](0.0, 0.0, 0.0);
vec3 mass_point = vec3(0.0);
int point_count = 0;
bool hard_cell = false; // any edge crossing on a convex brush face = exact Hermite
uint cell_material = 0u; // first brush material_id in this cell (route ② texturing)
uint cell_proj0 = 0u; // two candidate face projection indices (Stage-2 per-face texturing)
uint cell_proj1 = 0u;
// 12 edges of the cube at (x, y, z):
//
// 4 X-edges (dir=0):
// (x, y, z ) (x, y+1, z ) (x, y, z+1) (x, y+1, z+1)
// 4 Y-edges (dir=1):
// (x, y, z ) (x+1,y, z ) (x, y, z+1) (x+1,y, z+1)
// 4 Z-edges (dir=2):
// (x, y, z ) (x+1,y, z ) (x, y+1, z ) (x+1,y+1,z )
// Edge table: [12][4] = { ex_offset, ey_offset, ez_offset, direction }
const int edge_offsets[12 * 4] = int[12 * 4](
// 4 X-edges (dir=0)
0, 0, 0, 0,
0, 1, 0, 0,
0, 0, 1, 0,
0, 1, 1, 0,
// 4 Y-edges (dir=1)
0, 0, 0, 1,
1, 0, 0, 1,
0, 0, 1, 1,
1, 0, 1, 1,
// 4 Z-edges (dir=2)
0, 0, 0, 2,
1, 0, 0, 2,
0, 1, 0, 2,
1, 1, 0, 2
);
for (int e = 0; e < 12; e++) {
int ex = x + edge_offsets[e * 4 + 0];
int ey = y + edge_offsets[e * 4 + 1];
int ez = z + edge_offsets[e * 4 + 2];
int dir = edge_offsets[e * 4 + 3];
// Check if this edge has a sign change (with bounds check)
uint flag = get_edge_flag(ex, ey, ez, dir, V);
if (flag == 0u)
continue;
// Read intersection position and normal from compacted edge data.
// edge_data layout: 8 floats per edge [pos.x, pos.y, pos.z, norm.x, norm.y, norm.z, edge_idx, pad]
uint compact_idx = get_edge_scan(ex, ey, ez, dir, V);
uint base = compact_idx * 8u;
vec3 pos = vec3(edge_data[base + 0u],
edge_data[base + 1u],
edge_data[base + 2u]);
vec3 norm = vec3(edge_data[base + 3u],
edge_data[base + 4u],
edge_data[base + 5u]);
// Accumulate A^T A += outer(normal, normal)
ata[0] += norm.x * norm.x;
ata[1] += norm.x * norm.y;
ata[2] += norm.x * norm.z;
ata[3] += norm.y * norm.y;
ata[4] += norm.y * norm.z;
ata[5] += norm.z * norm.z;
// d = dot(normal, point)
float d = dot(norm, pos);
// A^T b += d * normal
atb[0] += d * norm.x;
atb[1] += d * norm.y;
atb[2] += d * norm.z;
// Accumulate mass point (running sum, divide later)
mass_point += pos;
point_count++;
// STEP 1: edge_data[base+7] packs (brush id << 12 | material_id). Nonzero =>
// this edge carries exact analytic Hermite (mark the cell hard); capture the
// first material for the vertex (route ② texturing).
uint packed7 = floatBitsToUint(edge_data[base + 7u]);
if (packed7 != 0u) {
hard_cell = true;
if (cell_material == 0u) cell_material = packed7 & 0xFFFu;
}
}
// No intersections gathered -- should not happen for an active voxel,
// but guard against it.
if (point_count == 0)
return;
// Compute mass point average
mass_point /= float(point_count);
// ---- Step 3: Solve QEF ----
// Regularise toward mass point: (A^T A + eps*I) x = A^T b + eps * mass_point
//
// When adaptive QEF is enabled (mesh_quality_flags bit 1), compute eigenvalues
// of ATA and apply strong regularization only in underconstrained directions.
// This fixes the rank-deficiency on flat/heightfield surfaces where XZ position
// is unconstrained (all normals point up).
float eps = float(QEF_REGULARISATION);
// Eigenvalue-adaptive regularization (Cardano's formula for 3x3 symmetric)
if ((mesh_quality_flags & 2u) != 0u) {
// Compute eigenvalues of ATA using the characteristic polynomial.
// For symmetric 3x3: lambda^3 - tr*lambda^2 + k*lambda - det = 0
// where tr = trace, k = sum of 2x2 minors on diagonal, det = determinant.
float tr = ata[0] + ata[3] + ata[5];
float k = ata[0]*ata[3] - ata[1]*ata[1]
+ ata[0]*ata[5] - ata[2]*ata[2]
+ ata[3]*ata[5] - ata[4]*ata[4];
float det_ata = ata[0]*(ata[3]*ata[5] - ata[4]*ata[4])
- ata[1]*(ata[1]*ata[5] - ata[4]*ata[2])
+ ata[2]*(ata[1]*ata[4] - ata[3]*ata[2]);
// Shift to depressed cubic: t^3 + p*t + q = 0 where t = lambda - tr/3
float tr3 = tr / 3.0;
float p = (3.0*k - tr*tr) / 3.0;
float q = (2.0*tr*tr*tr - 9.0*tr*k + 27.0*det_ata) / 27.0;
// Solve via trigonometric method (real roots guaranteed for symmetric matrix)
float disc = q*q/4.0 + p*p*p/27.0;
float lambda_max, lambda_min;
if (disc < 0.0 && p < -1e-10) {
float m = sqrt(-p / 3.0);
float theta = acos(clamp(-q / (2.0 * m*m*m), -1.0, 1.0)) / 3.0;
float e0 = tr3 + 2.0*m*cos(theta);
float e1 = tr3 + 2.0*m*cos(theta - 2.094395); // 2*pi/3
float e2 = tr3 + 2.0*m*cos(theta + 2.094395);
lambda_max = max(e0, max(e1, e2));
lambda_min = min(e0, min(e1, e2));
} else {
// Degenerate (all eigenvalues equal or nearly so)
lambda_max = tr3;
lambda_min = tr3;
}
// Condition number = lambda_max / lambda_min.
// When poorly conditioned (ratio > 10), boost regularization.
// This strongly biases toward mass point in underconstrained directions
// while preserving accuracy in well-constrained directions.
if (lambda_max > 1e-6) {
float condition = lambda_min / lambda_max;
// Scale eps from 0.01 (well-conditioned) to 1.0 (rank-deficient)
eps = mix(1.0, 0.01, clamp(condition * 10.0, 0.0, 1.0));
}
}
// NOTE: hard cells keep the STANDARD eps (NOT a tiny 1e-4). A near-zero eps is
// near-singular in the QEF's underconstrained (tangential) direction and
// AMPLIFIES float-rounding on OBLIQUE face normals into a large offset that the
// axis-aligned cell clamp then pushes OFF the plane — the "rotated-cube
// sandpaper" (verified: tests/test_rotated_face.cpp, single-plane face cells
// dev up to 0.16 at eps=1e-4). Sharpness comes from the plane-snap (Step 5b)
// projecting onto the EXACT plane/crease/corner, not from shrinking eps.
// Regularised matrix: M = A^T A + eps*I
float m00 = ata[0] + eps;
float m01 = ata[1];
float m02 = ata[2];
float m11 = ata[3] + eps;
float m12 = ata[4];
float m22 = ata[5] + eps;
// Right-hand side: r = A^T b - A^T A * mass_point
// A^T A * mass_point (using original unregularised ata):
vec3 ata_mp = vec3(
ata[0] * mass_point.x + ata[1] * mass_point.y + ata[2] * mass_point.z,
ata[1] * mass_point.x + ata[3] * mass_point.y + ata[4] * mass_point.z,
ata[2] * mass_point.x + ata[4] * mass_point.y + ata[5] * mass_point.z
);
vec3 rhs = vec3(atb[0], atb[1], atb[2]) - ata_mp;
// ---- Cramer's rule for 3x3 symmetric system ----
// M = [m00 m01 m02]
// [m01 m11 m12]
// [m02 m12 m22]
//
// det(M) using cofactor expansion along first row:
// det = m00*(m11*m22 - m12*m12) - m01*(m01*m22 - m12*m02) + m02*(m01*m12 - m11*m02)
//
// Cofactor matrix (for symmetric M):
// C00 = m11*m22 - m12*m12
// C01 = m02*m12 - m01*m22 (negated minor)
// C02 = m01*m12 - m02*m11
// C10 = C01 (symmetric cofactor for symmetric matrix)
// C11 = m00*m22 - m02*m02
// C12 = m01*m02 - m00*m12 (negated minor)
// C20 = C02
// C21 = C12
// C22 = m00*m11 - m01*m01
float c00 = m11 * m22 - m12 * m12;
float c01 = m02 * m12 - m01 * m22;
float c02 = m01 * m12 - m02 * m11;
float c11 = m00 * m22 - m02 * m02;
float c12 = m01 * m02 - m00 * m12;
float c22 = m00 * m11 - m01 * m01;
float det = m00 * c00 + m01 * c01 + m02 * c02;
vec3 result;
if (abs(det) > DET_THRESHOLD) {
// M^{-1} = (1/det) * cofactor^T (= (1/det) * cofactor for symmetric)
float inv_det = 1.0 / det;
// x_offset = M^{-1} * rhs
vec3 offset = vec3(
(c00 * rhs.x + c01 * rhs.y + c02 * rhs.z) * inv_det,
(c01 * rhs.x + c11 * rhs.y + c12 * rhs.z) * inv_det,
(c02 * rhs.x + c12 * rhs.y + c22 * rhs.z) * inv_det
);
result = mass_point + offset;
} else {
// Degenerate -- fall back to mass point
result = mass_point;
}
// ---- Step 4: Clamp vertex to voxel cell bounds ----
vec3 cell_min = cell_origin + vec3(float(x), float(y), float(z)) * voxel_size;
vec3 cell_max = cell_min + vec3(voxel_size);
result = clamp(result, cell_min, cell_max);
// ---- Step 4b: Vertex polish (fused from Stage 9) ----
// Clamped gradient descent toward the isosurface using the density grid.
// 3 iterations of Newton steps, clamped to 0.5*voxel_size for stability.
// STEP 1: SKIP for hard cells — polishing toward the SAMPLED density field
// would re-round the exact-Hermite vertex the QEF just placed exactly.
if (!hard_cell) {
int dim = V + 1; // density grid stride (33)
float max_step = 0.5 * voxel_size;
vec3 local = result - cell_origin;
for (int iter = 0; iter < 3; iter++) {
// Trilinear density sample
vec3 gc = clamp(local / voxel_size, vec3(0.0), vec3(float(dim - 2)));
ivec3 i0 = ivec3(floor(gc));
vec3 f = gc - vec3(i0);
float c000 = density[i0.z*dim*dim + i0.y*dim + i0.x];
float c100 = density[i0.z*dim*dim + i0.y*dim + (i0.x+1)];
float c010 = density[i0.z*dim*dim + (i0.y+1)*dim + i0.x];
float c110 = density[i0.z*dim*dim + (i0.y+1)*dim + (i0.x+1)];
float c001 = density[(i0.z+1)*dim*dim + i0.y*dim + i0.x];
float c101 = density[(i0.z+1)*dim*dim + i0.y*dim + (i0.x+1)];
float c011 = density[(i0.z+1)*dim*dim + (i0.y+1)*dim + i0.x];
float c111 = density[(i0.z+1)*dim*dim + (i0.y+1)*dim + (i0.x+1)];
float c00 = mix(c000, c100, f.x);
float c10 = mix(c010, c110, f.x);
float c01 = mix(c001, c101, f.x);
float c11 = mix(c011, c111, f.x);
float c0_v = mix(c00, c10, f.y);
float c1_v = mix(c01, c11, f.y);
float d = mix(c0_v, c1_v, f.z);
// Central differences gradient
float geps = voxel_size * 0.5;
vec3 grad;
// X gradient
vec3 gxp = clamp((local + vec3(geps,0,0)) / voxel_size, vec3(0.0), vec3(float(dim-2)));
vec3 gxn = clamp((local - vec3(geps,0,0)) / voxel_size, vec3(0.0), vec3(float(dim-2)));
ivec3 ixp = ivec3(floor(gxp)); vec3 fxp = gxp - vec3(ixp);
ivec3 ixn = ivec3(floor(gxn)); vec3 fxn = gxn - vec3(ixn);
float dxp = mix(mix(mix(density[ixp.z*dim*dim+ixp.y*dim+ixp.x], density[ixp.z*dim*dim+ixp.y*dim+ixp.x+1], fxp.x),
mix(density[ixp.z*dim*dim+(ixp.y+1)*dim+ixp.x], density[ixp.z*dim*dim+(ixp.y+1)*dim+ixp.x+1], fxp.x), fxp.y),
mix(mix(density[(ixp.z+1)*dim*dim+ixp.y*dim+ixp.x], density[(ixp.z+1)*dim*dim+ixp.y*dim+ixp.x+1], fxp.x),
mix(density[(ixp.z+1)*dim*dim+(ixp.y+1)*dim+ixp.x], density[(ixp.z+1)*dim*dim+(ixp.y+1)*dim+ixp.x+1], fxp.x), fxp.y), fxp.z);
float dxn = mix(mix(mix(density[ixn.z*dim*dim+ixn.y*dim+ixn.x], density[ixn.z*dim*dim+ixn.y*dim+ixn.x+1], fxn.x),
mix(density[ixn.z*dim*dim+(ixn.y+1)*dim+ixn.x], density[ixn.z*dim*dim+(ixn.y+1)*dim+ixn.x+1], fxn.x), fxn.y),
mix(mix(density[(ixn.z+1)*dim*dim+ixn.y*dim+ixn.x], density[(ixn.z+1)*dim*dim+ixn.y*dim+ixn.x+1], fxn.x),
mix(density[(ixn.z+1)*dim*dim+(ixn.y+1)*dim+ixn.x], density[(ixn.z+1)*dim*dim+(ixn.y+1)*dim+ixn.x+1], fxn.x), fxn.y), fxn.z);
grad.x = dxp - dxn;
// Y gradient
vec3 gyp = clamp((local + vec3(0,geps,0)) / voxel_size, vec3(0.0), vec3(float(dim-2)));
vec3 gyn = clamp((local - vec3(0,geps,0)) / voxel_size, vec3(0.0), vec3(float(dim-2)));
ivec3 iyp = ivec3(floor(gyp)); vec3 fyp = gyp - vec3(iyp);
ivec3 iyn = ivec3(floor(gyn)); vec3 fyn = gyn - vec3(iyn);
float dyp = mix(mix(mix(density[iyp.z*dim*dim+iyp.y*dim+iyp.x], density[iyp.z*dim*dim+iyp.y*dim+iyp.x+1], fyp.x),
mix(density[iyp.z*dim*dim+(iyp.y+1)*dim+iyp.x], density[iyp.z*dim*dim+(iyp.y+1)*dim+iyp.x+1], fyp.x), fyp.y),
mix(mix(density[(iyp.z+1)*dim*dim+iyp.y*dim+iyp.x], density[(iyp.z+1)*dim*dim+iyp.y*dim+iyp.x+1], fyp.x),
mix(density[(iyp.z+1)*dim*dim+(iyp.y+1)*dim+iyp.x], density[(iyp.z+1)*dim*dim+(iyp.y+1)*dim+iyp.x+1], fyp.x), fyp.y), fyp.z);
float dyn = mix(mix(mix(density[iyn.z*dim*dim+iyn.y*dim+iyn.x], density[iyn.z*dim*dim+iyn.y*dim+iyn.x+1], fyn.x),
mix(density[iyn.z*dim*dim+(iyn.y+1)*dim+iyn.x], density[iyn.z*dim*dim+(iyn.y+1)*dim+iyn.x+1], fyn.x), fyn.y),
mix(mix(density[(iyn.z+1)*dim*dim+iyn.y*dim+iyn.x], density[(iyn.z+1)*dim*dim+iyn.y*dim+iyn.x+1], fyn.x),
mix(density[(iyn.z+1)*dim*dim+(iyn.y+1)*dim+iyn.x], density[(iyn.z+1)*dim*dim+(iyn.y+1)*dim+iyn.x+1], fyn.x), fyn.y), fyn.z);
grad.y = dyp - dyn;
// Z gradient
vec3 gzp = clamp((local + vec3(0,0,geps)) / voxel_size, vec3(0.0), vec3(float(dim-2)));
vec3 gzn = clamp((local - vec3(0,0,geps)) / voxel_size, vec3(0.0), vec3(float(dim-2)));
ivec3 izp = ivec3(floor(gzp)); vec3 fzp = gzp - vec3(izp);
ivec3 izn = ivec3(floor(gzn)); vec3 fzn = gzn - vec3(izn);
float dzp = mix(mix(mix(density[izp.z*dim*dim+izp.y*dim+izp.x], density[izp.z*dim*dim+izp.y*dim+izp.x+1], fzp.x),
mix(density[izp.z*dim*dim+(izp.y+1)*dim+izp.x], density[izp.z*dim*dim+(izp.y+1)*dim+izp.x+1], fzp.x), fzp.y),
mix(mix(density[(izp.z+1)*dim*dim+izp.y*dim+izp.x], density[(izp.z+1)*dim*dim+izp.y*dim+izp.x+1], fzp.x),
mix(density[(izp.z+1)*dim*dim+(izp.y+1)*dim+izp.x], density[(izp.z+1)*dim*dim+(izp.y+1)*dim+izp.x+1], fzp.x), fzp.y), fzp.z);
float dzn = mix(mix(mix(density[izn.z*dim*dim+izn.y*dim+izn.x], density[izn.z*dim*dim+izn.y*dim+izn.x+1], fzn.x),
mix(density[izn.z*dim*dim+(izn.y+1)*dim+izn.x], density[izn.z*dim*dim+(izn.y+1)*dim+izn.x+1], fzn.x), fzn.y),
mix(mix(density[(izn.z+1)*dim*dim+izn.y*dim+izn.x], density[(izn.z+1)*dim*dim+izn.y*dim+izn.x+1], fzn.x),
mix(density[(izn.z+1)*dim*dim+(izn.y+1)*dim+izn.x], density[(izn.z+1)*dim*dim+(izn.y+1)*dim+izn.x+1], fzn.x), fzn.y), fzn.z);
grad.z = dzp - dzn;
float grad_sq = dot(grad, grad);
if (grad_sq < 1e-12) break;
float step = clamp(d / grad_sq, -max_step, max_step);
local -= step * grad;
}
// Re-clamp to cell bounds after polish
local = clamp(local, cell_min - cell_origin, cell_max - cell_origin);
result = local + cell_origin;
}
// ---- Step 5: shading normal (all edges) ----
vec3 avg_normal = vec3(0.0);
for (int e = 0; e < 12; e++) {
int ex = x + edge_offsets[e * 4 + 0];
int ey = y + edge_offsets[e * 4 + 1];
int ez = z + edge_offsets[e * 4 + 2];
int dir = edge_offsets[e * 4 + 3];
if (get_edge_flag(ex, ey, ez, dir, V) == 0u) continue;
uint base = get_edge_scan(ex, ey, ez, dir, V) * 8u;
avg_normal += vec3(edge_data[base + 3u], edge_data[base + 4u], edge_data[base + 5u]);
}
avg_normal = normalize(avg_normal);
// ---- Step 5b: FEATURE-AWARE PLANE-SET SNAP ----
// The old cluster (grid crossings) missed the crease on cells the edge only
// grazed on one face → angled-edge SAWTOOTH. Instead, capture the EXACT local
// plane-set straight from the brush geometry: brush faces within T of the cell
// CENTRE where the cell is inside all OTHER planes (within T). That yields the
// face (1) / crease (2) / corner (3) even for cells the grid crossings under-
// sampled. Snap the QEF vertex to it, and take the material from the nearest
// face (snap brush → no bleed). Verified: tests/test_planeset.cpp (0 err, edge
// cells get the crease). Only convex-brush (hard) cells; organic untouched.
vec3 cl_n[3];
float cl_d[3];
int cl_count = 0;
if (hard_cell) {
// READ the baked plane-set brick (bake_planeset.glsl captured it this pass,
// with csg4 hidden-face removal). No per-cell re-derivation here anymore.
int Vg = grid_dims.x;
uint ci = uint(z) * uint(Vg * Vg) + uint(y) * uint(Vg) + uint(x);
CellPlaneSet ps = cell_planeset[ci];
cl_count = int(ps.cl_count);
cell_material = ps.material; // dominant-brush material (bake picks smallest union-SDF)
cell_proj0 = ps.proj0; // two candidate face projections (Stage-2 per-face texturing)
cell_proj1 = ps.proj1;
if (cl_count > 0) { cl_n[0] = ps.planes[0].xyz; cl_d[0] = ps.planes[0].w; }
if (cl_count > 1) { cl_n[1] = ps.planes[1].xyz; cl_d[1] = ps.planes[1].w; }
if (cl_count > 2) { cl_n[2] = ps.planes[2].xyz; cl_d[2] = ps.planes[2].w; }
}
if (cl_count >= 1) {
vec3 snapped = result;
bool snap_ok = false; // a SEGMENT-CLAMPed (or in-cell corner) vertex needs
// no escape fallback — it is already in-cell AND
// exactly on the feature.
if (cl_count == 1) {
// Face: project onto the plane.
snapped = result - cl_n[0] * (dot(cl_n[0], result) - cl_d[0]);
} else {
// ---- SEGMENT CLAMP (replaces "snap to the INFINITE crease line, then
// clamp the 3D point to the axis-aligned cell box").
//
// The old path was exact right up to the clamp and then thrown off the
// surface BY it: on an OBLIQUE crease the correct point leaves its own
// cell in ~22% of crease cells (by up to 0.52 voxels), and clamping a
// 3D point to the box pins it to a voxel-grid-aligned face — which IS
// the stair-step, quantised to the voxel grid. Axis-aligned creases run
// parallel to the grid and never leave their cell, which is exactly why
// the artifact only ever showed on rotated geometry.
//
// Fix: stay ON the crease line and constrain only the coordinate ALONG
// it, to (line n cell) n (line n the owning brush's other half-spaces).
// In-cell by construction AND exactly on the feature.
//
// Oracle: tests/test_crease_segment_oracle.cpp (rotated cube, half=5,
// voxel 0.224, 604 crease cells): p95 0.01011 -> 0.00000,
// MAX 0.08990 -> 0.00000, containment violations 0. Axis control stays
// 0.00000. Corner cells (the worst case) 0.08990 -> 0.00000.
int ia = 0, ib = 1;
if (cl_count == 3) {
// If the true 3-plane corner lies in THIS cell, it IS the feature.
vec3 a = cl_n[0], b = cl_n[1], cn = cl_n[2];
float denom = dot(a, cross(b, cn));
if (abs(denom) > 1e-5) {
vec3 corner = (cl_d[0] * cross(b, cn) + cl_d[1] * cross(cn, a)
+ cl_d[2] * cross(a, b)) / denom;
if (all(greaterThanEqual(corner, cell_min)) &&
all(lessThanEqual(corner, cell_max))) {
snapped = corner;
snap_ok = true;
}
}
if (!snap_ok) {
// Near a corner but not AT it => this cell is really on an EDGE.
// Snapping all three planes to the single corner point was the
// worst case measured (cl==3 p95 0.08990 vs 0.00857 for cl==2).
// Take the crease pair whose line passes nearest the centre.
vec3 ctr = cell_min + vec3(voxel_size * 0.5);
float best = 1e30;
for (int k = 0; k < 3; k++) {
int p = k, q = (k + 1) % 3;
float cc = dot(cl_n[p], cl_n[q]);
float dd = 1.0 - cc * cc;
if (dd <= 1e-4) continue;
float r0 = cl_d[p] - dot(cl_n[p], ctr);
float r1 = cl_d[q] - dot(cl_n[q], ctr);
vec3 pt = ctr + cl_n[p] * ((r0 - cc * r1) / dd)
+ cl_n[q] * ((r1 - cc * r0) / dd);
float dist = length(pt - ctr);
if (dist < best) { best = dist; ia = p; ib = q; }
}
}
}
if (!snap_ok) {
float cc = dot(cl_n[ia], cl_n[ib]);
float dd = 1.0 - cc * cc;
if (dd > 1e-4) {
float r0 = cl_d[ia] - dot(cl_n[ia], result);
float r1 = cl_d[ib] - dot(cl_n[ib], result);
vec3 pline = result + cl_n[ia] * ((r0 - cc * r1) / dd)
+ cl_n[ib] * ((r1 - cc * r0) / dd);
vec3 ldir = cross(cl_n[ia], cl_n[ib]);
float ll = length(ldir);
if (ll > 1e-9) {
ldir /= ll;
float tmin = -1e30, tmax = 1e30;
bool ok = true;
// clip the line to the cell box (t measured from pline)
for (int a2 = 0; a2 < 3; a2++) {
float pv = pline[a2], dv = ldir[a2];
if (abs(dv) < 1e-9) {
if (pv < cell_min[a2] || pv > cell_max[a2]) ok = false;
} else {
float ta = (cell_min[a2] - pv) / dv;
float tb = (cell_max[a2] - pv) / dv;
tmin = max(tmin, min(ta, tb));
tmax = min(tmax, max(ta, tb));
}
}
// Clip to the OWNING brush's remaining half-spaces so the
// vertex cannot slide past the end of a finite edge. Uses
// the SAME bindings edge_intersect already reads (9 =
// gpu_edits, 10 = convex_planes) — no new bake data, no
// buffer-layout change. Measured worth it on tight brushes:
// corner MAX 0.084 -> 0.045 at half=2 (identical at half=5).
if (ok) {
for (uint ei = 0u; ei < edit_count; ei++) {
GpuEdit ce = gpu_edits[ei];
if (ce.edit_type != 8u && ce.edit_type != 9u) continue;
vec3 dcv = pline - vec3(ce.pos_x, ce.pos_y, ce.pos_z);
if (dot(dcv, dcv) > ce.radius * ce.radius) continue;
uint poff = ce._pad1, pcnt = ce._pad2;
// owner = the brush carrying BOTH captured planes.
// Match the PLANE (normal AND offset), not just the
// normal: Quake brushes are overwhelmingly axis
// aligned, so a normal-only test matches hundreds of
// unrelated brushes and the first one in the list
// wins. Measured on START.MAP
// (tests/test_startmap_crease.cpp): normal-only gave
// 123214 wrong-brush matches and left the segment
// clamp unresolved on 64538 crease cells (p95 0.125,
// i.e. no better than the old path); matching the
// plane gives 14872 / 3755 and p95 0.00000.
bool has_a = false, has_b = false;
for (uint pi = 0u; pi < pcnt; pi++) {
vec4 q = convex_planes[poff + pi];
if (dot(q.xyz, cl_n[ia]) > 0.999 && abs(q.w - cl_d[ia]) < 1e-3) has_a = true;
if (dot(q.xyz, cl_n[ib]) > 0.999 && abs(q.w - cl_d[ib]) < 1e-3) has_b = true;
}
if (!has_a || !has_b) continue;
for (uint pi = 0u; pi < pcnt; pi++) {
vec4 q = convex_planes[poff + pi];
if (dot(q.xyz, cl_n[ia]) > 0.999 && abs(q.w - cl_d[ia]) < 1e-3) continue;
if (dot(q.xyz, cl_n[ib]) > 0.999 && abs(q.w - cl_d[ib]) < 1e-3) continue;
float aa = dot(q.xyz, ldir);
float bb = q.w - dot(q.xyz, pline);
if (abs(aa) < 1e-9) { if (bb < -1e-6) ok = false; }
else if (aa > 0.0) { tmax = min(tmax, bb / aa); }
else { tmin = max(tmin, bb / aa); }
}
break;
}
}
if (ok && tmin <= tmax) {
snapped = pline + ldir * clamp(0.0, tmin, tmax);
snap_ok = true;
}
}
}
}
}
// ESCAPE FALLBACK + CLAMP (Phase 1 fix — oracle tests/test_brick_dense:
// red 15.8% -> 0.04%, out-of-cell spikes 14.6% -> 0, surf-dev ~0). A crease/
// corner snap whose target leaves the cell is a far/degenerate intersection
// (the OLD 1-voxel cap turned these into jagged QEF verts = sawtooth, and an
// uncapped snap would SPIKE outside the cell). Instead project onto the single
// NEAREST face — stays IN-cell AND on the true surface — then clamp so nothing
// can spike. Always accept: fallback+clamp guarantees a valid in-cell vertex.
//
// Now reached ONLY when the segment clamp above could not produce a vertex
// (near-parallel planes, or the crease line does not actually cross this
// cell — the latter is the FALSE-CREASE case: cl_count==2 captured by the
// T=1.2*voxel centre test on a cell that has no crease in it). A successful
// segment clamp is already in-cell and exactly on the feature, so running
// this would only move it back off.
if (cl_count >= 2 && !snap_ok) {
vec3 lo = cell_min - vec3(voxel_size * 0.5);
vec3 hi = cell_max + vec3(voxel_size * 0.5);
if (any(lessThan(snapped, lo)) || any(greaterThan(snapped, hi))) {
vec3 cc2 = cell_min + vec3(voxel_size * 0.5);
int nb = 0;
float ba = 1e30;
for (int c = 0; c < cl_count; c++) {
float aa = abs(dot(cl_n[c], cc2) - cl_d[c]);
if (aa < ba) { ba = aa; nb = c; }
}
snapped = result - cl_n[nb] * (dot(cl_n[nb], result) - cl_d[nb]);
}
}
result = clamp(snapped, cell_min, cell_max);
}
// ---- Step 5c: CONTAINMENT ----
// Nothing above keeps the vertex INSIDE the union solid. A vertex that pokes
// out is exactly what shows as a silhouette notch, and deviation-from-the-
// crease-line is structurally blind to it — which is why the crease fix
// measured perfect while the artifact was still visible in game.
//
// Naive Surface Nets is immune to this (a centroid of on-cell-edge crossings
// is necessarily inside the local solid, measured 0 protrusion at every voxel
// size) but has no sharp features. This enforces the same property explicitly
// while keeping them: if the vertex is outside the union, pull it back onto
// the nearest brush's active face and re-clamp to the cell.
//
// Oracle (tests/test_startmap_crease.cpp, real 1164-brush START.MAP), cells
// protruding >0.1 voxel: 1070 before the crease fix, 275 after it, 27 with
// this. p95 protrusion 0.00090 -> 0.00056 -> 0.00013. Crease-line deviation
// is UNCHANGED (p50/p95 0.00000), so this costs no feature accuracy.
// NOTE the MAX (0.224) does not move — the count collapses, the tail does not.
if (hard_cell && cl_count >= 1) {
float f_union = 1e30;
int best_brush = -1;
for (uint ei = 0u; ei < edit_count; ei++) {
GpuEdit ce = gpu_edits[ei];
if (ce.edit_type != 8u && ce.edit_type != 9u) continue;
uint poff = ce._pad1, pcnt = ce._pad2;
float f = -1e30; // convex SDF: max_j(n.p - d)
for (uint pi = 0u; pi < pcnt; pi++) {
vec4 pl = convex_planes[poff + pi];
f = max(f, dot(pl.xyz, result) - pl.w);
}
if (f < f_union) { f_union = f; best_brush = int(ei); }
}
// f_union > 0 => outside every brush => outside the union.
if (best_brush >= 0 && f_union > 1e-6) {
GpuEdit cb = gpu_edits[uint(best_brush)];
uint poff = cb._pad1, pcnt = cb._pad2;
float best = -1e30;
vec4 bp = vec4(0.0);
for (uint pi = 0u; pi < pcnt; pi++) {
vec4 pl = convex_planes[poff + pi];
float dp = dot(pl.xyz, result) - pl.w;
if (dp > best) { best = dp; bp = pl; } // active face
}
if (best > -1e29) {
vec3 proj = result - bp.xyz * (dot(bp.xyz, result) - bp.w);
result = clamp(proj, cell_min, cell_max);
}
}
}
// ---- Step 6: Seam flag ----
// Boundary voxel if any coordinate is 0 or V-2 (VOXELS_PER_CELL - 2)
uint seam_flag = 0u;
if (x == 0 || y == 0 || z == 0 ||
x == V - 2 || y == V - 2 || z == V - 2) {
seam_flag = 1u;
}
// ---- Step 7: Write output ----
// 12 floats per vertex: pos(3) + normal(3) + corner_signs + seam_flag +
// material_id + voxel_idx + pad(2)
uint out_idx = voxel_scan[gid] * 12u;
vertices[out_idx + 0u] = result.x;
vertices[out_idx + 1u] = result.y;
vertices[out_idx + 2u] = result.z;
vertices[out_idx + 3u] = avg_normal.x;
vertices[out_idx + 4u] = avg_normal.y;
vertices[out_idx + 5u] = avg_normal.z;
vertices[out_idx + 6u] = uintBitsToFloat(corner_signs);
vertices[out_idx + 7u] = uintBitsToFloat(seam_flag);
vertices[out_idx + 8u] = uintBitsToFloat(cell_material); // material_id (route ②)
vertices[out_idx + 9u] = uintBitsToFloat(gid); // voxel_idx
vertices[out_idx + 10u] = uintBitsToFloat(cell_proj0); // per-face projection candidate 0 (Stage-2)
vertices[out_idx + 11u] = uintBitsToFloat(cell_proj1); // per-face projection candidate 1 (crease other face)
}