// CPU reference implementation + tests for edge_classify.glsl (Stage 2).
// Build: g++ -std=c++17 -O2 test/test_edge_classify.cpp -o test/test_edge_classify -lm && ./test/test_edge_classify
#include <cmath>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <cstdint>
#include <vector>
static int g_pass = 0;
static int g_fail = 0;
static void check(bool cond, const char* name, const char* detail = nullptr) {
if (cond) {
printf(" PASS %s\n", name);
g_pass++;
} else {
printf(" FAIL %s", name);
if (detail) printf(" (%s)", detail);
printf("\n");
g_fail++;
}
}
// -------------------------------------------------------------------
// CPU reference: mirrors edge_classify.glsl exactly
// -------------------------------------------------------------------
static constexpr int V = 32; // grid voxels per axis
static constexpr int S = V + 1; // density grid stride (33)
static constexpr int TOTAL = V * V * V; // 32768 voxels
// Density index: z * S * S + y * S + x
static inline int dens_idx(int x, int y, int z) {
return z * S * S + y * S + x;
}
// Run edge classification on density[], write to edge_flags[]
// density: S*S*S floats (33³ = 35937)
// edge_flags: TOTAL*3 uints (32768 * 3 = 98304)
static void classify_edges(const float* density, uint32_t* edge_flags) {
for (int gid = 0; gid < TOTAL; gid++) {
int z = gid / (V * V);
int y = (gid % (V * V)) / V;
int x = gid % V;
float d0 = density[dens_idx(x, y, z)];
// Edge 0: +X
float dx = density[dens_idx(x + 1, y, z)];
edge_flags[gid * 3 + 0] = ((d0 >= 0.0f) != (dx >= 0.0f)) ? 1u : 0u;
// Edge 1: +Y
float dy = density[dens_idx(x, y + 1, z)];
edge_flags[gid * 3 + 1] = ((d0 >= 0.0f) != (dy >= 0.0f)) ? 1u : 0u;
// Edge 2: +Z
float dz = density[dens_idx(x, y, z + 1)];
edge_flags[gid * 3 + 2] = ((d0 >= 0.0f) != (dz >= 0.0f)) ? 1u : 0u;
}
}
// Count total active edges
static int count_active(const uint32_t* edge_flags, int count) {
int total = 0;
for (int i = 0; i < count; i++)
total += edge_flags[i];
return total;
}
// Count active edges for a specific direction (0=X, 1=Y, 2=Z)
static int count_active_dir(const uint32_t* edge_flags, int dir) {
int total = 0;
for (int gid = 0; gid < TOTAL; gid++)
total += edge_flags[gid * 3 + dir];
return total;
}
// -------------------------------------------------------------------
// Test 1: All-positive density => no sign changes => no active edges
// -------------------------------------------------------------------
static void test_all_positive() {
printf("\n--- All-positive density ---\n");
std::vector<float> density(S * S * S, 1.0f);
std::vector<uint32_t> edge_flags(TOTAL * 3, 99u);
classify_edges(density.data(), edge_flags.data());
int active = count_active(edge_flags.data(), TOTAL * 3);
char buf[128];
snprintf(buf, sizeof(buf), "active edges = %d", active);
check(active == 0, "all-positive density => 0 active edges", buf);
// Verify all flags are exactly 0
bool all_zero = true;
for (int i = 0; i < TOTAL * 3; i++) {
if (edge_flags[i] != 0u) { all_zero = false; break; }
}
check(all_zero, "all flags are exactly 0");
}
// -------------------------------------------------------------------
// Test 2: All-negative density => no sign changes => no active edges
// -------------------------------------------------------------------
static void test_all_negative() {
printf("\n--- All-negative density ---\n");
std::vector<float> density(S * S * S, -1.0f);
std::vector<uint32_t> edge_flags(TOTAL * 3, 99u);
classify_edges(density.data(), edge_flags.data());
int active = count_active(edge_flags.data(), TOTAL * 3);
char buf[128];
snprintf(buf, sizeof(buf), "active edges = %d", active);
check(active == 0, "all-negative density => 0 active edges", buf);
}
// -------------------------------------------------------------------
// Test 3: Plane at z=16 => exactly 1024 Z-axis crossings
//
// density(x,y,z) = z - 16.0 (negative below z=16, positive at/above z=16)
// Sign change on +Z edge at (x,y,z) when sign(d[x,y,z]) != sign(d[x,y,z+1])
// This happens at z=15: d[15] = -1 (neg), d[16] = 0 (non-neg) => sign change
// For all 32*32 = 1024 (x,y) pairs, exactly one Z-axis crossing each.
// No X-axis or Y-axis crossings (density depends only on z).
// -------------------------------------------------------------------
static void test_plane_z16() {
printf("\n--- Plane at z=16 ---\n");
std::vector<float> density(S * S * S);
for (int z = 0; z < S; z++)
for (int y = 0; y < S; y++)
for (int x = 0; x < S; x++)
density[dens_idx(x, y, z)] = static_cast<float>(z) - 16.0f;
std::vector<uint32_t> edge_flags(TOTAL * 3, 0u);
classify_edges(density.data(), edge_flags.data());
int x_crossings = count_active_dir(edge_flags.data(), 0);
int y_crossings = count_active_dir(edge_flags.data(), 1);
int z_crossings = count_active_dir(edge_flags.data(), 2);
char buf[128];
snprintf(buf, sizeof(buf), "got %d", x_crossings);
check(x_crossings == 0, "plane z=16: 0 X-axis crossings", buf);
snprintf(buf, sizeof(buf), "got %d", y_crossings);
check(y_crossings == 0, "plane z=16: 0 Y-axis crossings", buf);
// z=15 is the only row where sign changes: d[15]=-1 (neg), d[16]=0 (non-neg)
// 32 * 32 = 1024 voxels at z=15
snprintf(buf, sizeof(buf), "got %d, expected 1024", z_crossings);
check(z_crossings == 1024, "plane z=16: 1024 Z-axis crossings", buf);
// Verify crossings are all at z=15
int z15_count = 0;
for (int iy = 0; iy < V; iy++) {
for (int ix = 0; ix < V; ix++) {
int gid = 15 * V * V + iy * V + ix;
if (edge_flags[gid * 3 + 2] == 1u) z15_count++;
}
}
snprintf(buf, sizeof(buf), "z=15 crossings = %d", z15_count);
check(z15_count == 1024, "all Z crossings at z=15", buf);
}
// -------------------------------------------------------------------
// Test 4: Plane at y=0 => Y-axis crossings, no X or Z
// -------------------------------------------------------------------
static void test_plane_y0() {
printf("\n--- Plane at y=0 ---\n");
// density = y - 0.5 => sign change between y=0 (d=-0.5) and y=1 (d=0.5)
// Actually: at y=0, d = -0.5 (neg); at y=1, d = 0.5 (pos) => change at y=0
std::vector<float> density(S * S * S);
for (int z = 0; z < S; z++)
for (int y = 0; y < S; y++)
for (int x = 0; x < S; x++)
density[dens_idx(x, y, z)] = static_cast<float>(y) - 0.5f;
std::vector<uint32_t> edge_flags(TOTAL * 3, 0u);
classify_edges(density.data(), edge_flags.data());
int x_crossings = count_active_dir(edge_flags.data(), 0);
int y_crossings = count_active_dir(edge_flags.data(), 1);
int z_crossings = count_active_dir(edge_flags.data(), 2);
char buf[128];
snprintf(buf, sizeof(buf), "got %d", x_crossings);
check(x_crossings == 0, "plane y=0: 0 X-axis crossings", buf);
// y=0: d=-0.5 (neg), y=1: d=0.5 (pos) => 32*32 = 1024 crossings at y=0
snprintf(buf, sizeof(buf), "got %d, expected 1024", y_crossings);
check(y_crossings == 1024, "plane y=0: 1024 Y-axis crossings", buf);
snprintf(buf, sizeof(buf), "got %d", z_crossings);
check(z_crossings == 0, "plane y=0: 0 Z-axis crossings", buf);
}
// -------------------------------------------------------------------
// Test 5: Sphere centered at (16,16,16) radius 8 => some edges, not all
// -------------------------------------------------------------------
static void test_sphere() {
printf("\n--- Sphere r=8 at (16,16,16) ---\n");
const float cx = 16.0f, cy = 16.0f, cz = 16.0f, r = 8.0f;
std::vector<float> density(S * S * S);
for (int z = 0; z < S; z++)
for (int y = 0; y < S; y++)
for (int x = 0; x < S; x++) {
float dx = static_cast<float>(x) - cx;
float dy = static_cast<float>(y) - cy;
float dz = static_cast<float>(z) - cz;
// Negative inside sphere (air), positive outside (solid)
density[dens_idx(x, y, z)] = dx * dx + dy * dy + dz * dz - r * r;
}
std::vector<uint32_t> edge_flags(TOTAL * 3, 0u);
classify_edges(density.data(), edge_flags.data());
int total_active = count_active(edge_flags.data(), TOTAL * 3);
int x_active = count_active_dir(edge_flags.data(), 0);
int y_active = count_active_dir(edge_flags.data(), 1);
int z_active = count_active_dir(edge_flags.data(), 2);
char buf[256];
// Sphere should have *some* active edges
snprintf(buf, sizeof(buf), "total = %d", total_active);
check(total_active > 0, "sphere has active edges", buf);
// But definitely not all edges
snprintf(buf, sizeof(buf), "total = %d of %d possible", total_active, TOTAL * 3);
check(total_active < TOTAL * 3 / 2, "sphere has fewer than half edges active", buf);
// All 3 directions should have crossings (sphere is symmetric)
snprintf(buf, sizeof(buf), "X=%d, Y=%d, Z=%d", x_active, y_active, z_active);
check(x_active > 0 && y_active > 0 && z_active > 0,
"sphere has crossings in all 3 directions", buf);
// Symmetry: X, Y, Z crossing counts should be equal (sphere centered on grid)
check(x_active == y_active && y_active == z_active,
"sphere crossings symmetric across axes", buf);
// Reasonable edge count for a sphere of radius 8 in a 32³ grid
// Surface area ~ 4*pi*r² ~ 804, each surface element crosses ~1 edge
// Expect roughly 800-2000 crossings per axis
snprintf(buf, sizeof(buf), "per-axis = %d", x_active);
check(x_active > 200 && x_active < 5000,
"sphere crossing count in reasonable range", buf);
}
// -------------------------------------------------------------------
// Test 6: Checkerboard pattern => maximum sign changes
// -------------------------------------------------------------------
static void test_checkerboard() {
printf("\n--- Checkerboard pattern ---\n");
std::vector<float> density(S * S * S);
for (int z = 0; z < S; z++)
for (int y = 0; y < S; y++)
for (int x = 0; x < S; x++) {
// Alternating positive/negative
bool even = ((x + y + z) % 2) == 0;
density[dens_idx(x, y, z)] = even ? 1.0f : -1.0f;
}
std::vector<uint32_t> edge_flags(TOTAL * 3, 0u);
classify_edges(density.data(), edge_flags.data());
int total_active = count_active(edge_flags.data(), TOTAL * 3);
// Every single edge should be a sign change in a checkerboard
char buf[128];
snprintf(buf, sizeof(buf), "active = %d of %d", total_active, TOTAL * 3);
check(total_active == TOTAL * 3, "checkerboard: all edges active", buf);
}
// -------------------------------------------------------------------
// Test 7: Zero density => no sign changes (>= 0 on both sides)
// -------------------------------------------------------------------
static void test_all_zero() {
printf("\n--- All-zero density ---\n");
std::vector<float> density(S * S * S, 0.0f);
std::vector<uint32_t> edge_flags(TOTAL * 3, 99u);
classify_edges(density.data(), edge_flags.data());
int active = count_active(edge_flags.data(), TOTAL * 3);
char buf[128];
snprintf(buf, sizeof(buf), "active edges = %d", active);
check(active == 0, "all-zero density => 0 active edges (0.0 >= 0.0 on both sides)", buf);
}
// -------------------------------------------------------------------
// main
// -------------------------------------------------------------------
int main() {
printf("=== Edge Classification Tests ===\n");
test_all_positive();
test_all_negative();
test_plane_z16();
test_plane_y0();
test_sphere();
test_checkerboard();
test_all_zero();
printf("\n=== Results: %d passed, %d failed ===\n", g_pass, g_fail);
return g_fail > 0 ? 1 : 0;
}