// ============================================================================
// IS THE FILTER'S OWN SLOPE THE GRADIENT OF ITS OWN HEIGHT?
//
// erosion_filter accumulates height and both slope components together and
// returns all three; every caller used the height and discarded the slope.
// edge_intersect.glsl then rebuilds that gradient by central differences --
// six more full evaluations per edge, each seven octaves of phacelle noise.
// Measured: the octave loop takes a 32^3 grid from 491 us to 3795 us, and
// edge_intersect is 91% of the mesh pipeline.
//
// Before any of that is replaced, the accumulated slope has to BE the gradient.
// This compares erosion_node_eval_grad's .z/.w against central differences of
// erosion_node_eval_full's .x over a field of points, in float64 steps so the
// reference is not itself limited by the thing being tested.
//
// It can fail: a deliberately wrong gradient is checked too, and must be
// rejected. A gate that only ever sees the correct answer proves nothing --
// this session has already produced three harnesses that returned the same
// number for every input.
//
// Build (no engine, no Vulkan):
// g++ -O2 -std=c++17 -I ../proc-nodes/src/shared \
// -o /tmp/erosion_gradient tests/test_erosion_gradient.cpp
// Run: /tmp/erosion_gradient
// ============================================================================
#include <cstdio>
#include <cmath>
#include <algorithm>
#include "erosion_cpp_shim.h" // brings in erosion_filter.inc as C++
using dcerosion::vec2;
using dcerosion::vec3;
using dcerosion::vec4;
// The map's own erosion settings (map_dc_diffusion.gd).
static const float HS = 5.0f;
static const float STR = 0.22f;
static const float CELL = 0.7f;
static const int OCT = 7;
static const float GAIN = 0.5f;
static const float FDIV = 4.8f;
static const float ESCL = 0.15f;
// A stand-in for the heightfield: smooth, with real slope, and NOT axis
// aligned, so a gradient that is right in x and wrong in z cannot pass.
static vec3 input_gn(float px, float pz) {
double h = 0.6 * std::sin(px * 3.1) * std::cos(pz * 2.3) + 0.25 * std::sin(pz * 5.7);
double dx = 0.6 * 3.1 * std::cos(px * 3.1) * std::cos(pz * 2.3);
double dz = -0.6 * 2.3 * std::sin(px * 3.1) * std::sin(pz * 2.3) + 0.25 * 5.7 * std::cos(pz * 5.7);
return vec3((float)h, (float)dx, (float)dz);
}
static float height_at(float px, float pz) {
vec2 r = dcerosion::erosion_node_eval_full(vec2(px, pz), input_gn(px, pz),
HS, STR, CELL, OCT, GAIN, FDIV, ESCL,
0.25f, 2.0f, 0.5f, 1.5f, 0.5f, vec2(-0.65f, 0.0f),
vec4(0.1f, 0.0f, 0.1f, 2.0f), vec4(1.25f, 1.25f, 2.8f, 1.5f), vec2(0.7f, 1.0f), 999.0f);
return r.x;
}
int main() {
printf("analytic slope from erosion_node_eval_grad vs central differences of\n");
printf("erosion_node_eval_full, %d octaves, strength %.2f\n\n", OCT, STR);
printf(" %-22s %12s %12s %12s\n", "point", "analytic dx", "numeric dx", "rel err");
const float eps = 1e-3f; // in p units; p = world / terrain_height
double worst_rel = 0.0;
double sum_rel = 0.0;
int n = 0;
int reported = 0;
for (int i = 0; i < 12; i++) {
for (int j = 0; j < 12; j++) {
float px = -0.5f + 0.0834f * (float)i;
float pz = -0.5f + 0.0834f * (float)j;
vec4 g = dcerosion::erosion_node_eval_grad(vec2(px, pz), input_gn(px, pz),
HS, STR, CELL, OCT, GAIN, FDIV, ESCL,
0.25f, 2.0f, 0.5f, 1.5f, 0.5f, vec2(-0.65f, 0.0f),
vec4(0.1f, 0.0f, 0.1f, 2.0f), vec4(1.25f, 1.25f, 2.8f, 1.5f), vec2(0.7f, 1.0f), 999.0f);
// The height must be unchanged from the vec2 form, or the two
// functions have already diverged and the slope is moot.
float h_ref = height_at(px, pz);
if (std::fabs(g.x - h_ref) > 1e-4f) {
printf(" HEIGHT MISMATCH at (%.3f, %.3f): %.6f vs %.6f\n", px, pz, g.x, h_ref);
return 2;
}
double num_dx = (height_at(px + eps, pz) - height_at(px - eps, pz)) / (2.0 * eps);
double num_dz = (height_at(px, pz + eps) - height_at(px, pz - eps)) / (2.0 * eps);
double scale = std::max(1.0, std::max(std::fabs(num_dx), std::fabs(num_dz)));
double rx = std::fabs(g.z - num_dx) / scale;
double rz = std::fabs(g.w - num_dz) / scale;
double rel = std::max(rx, rz);
sum_rel += rel; n++;
if (rel > worst_rel) worst_rel = rel;
if (reported < 6) {
printf(" (%+.3f, %+.3f) %12.5f %12.5f %12.5f\n", px, pz, g.z, num_dx, rel);
reported++;
}
}
}
printf("\n points %d mean rel err %.6f worst %.6f\n", n, sum_rel / n, worst_rel);
// THE NEGATIVE CONTROL. A gradient scaled by 1.5 must be REJECTED, or this
// gate cannot tell a right answer from a wrong one.
{
float px = 0.13f, pz = -0.21f;
vec4 g = dcerosion::erosion_node_eval_grad(vec2(px, pz), input_gn(px, pz),
HS, STR, CELL, OCT, GAIN, FDIV, ESCL,
0.25f, 2.0f, 0.5f, 1.5f, 0.5f, vec2(-0.65f, 0.0f),
vec4(0.1f, 0.0f, 0.1f, 2.0f), vec4(1.25f, 1.25f, 2.8f, 1.5f), vec2(0.7f, 1.0f), 999.0f);
double num_dx = (height_at(px + eps, pz) - height_at(px - eps, pz)) / (2.0 * eps);
double bad = g.z * 1.5;
double scale = std::max(1.0, std::fabs(num_dx));
double rel_bad = std::fabs(bad - num_dx) / scale;
printf(" negative control: a 1.5x gradient reads rel err %.6f -> %s\n",
rel_bad, rel_bad > 0.02 ? "REJECTED (gate is alive)" : "ACCEPTED (GATE IS BLIND)");
if (rel_bad <= 0.02) return 3;
}
// The filter is piecewise-smooth: sign(ph.y) in the straight-gully term is a
// step, so a central difference straddling a gully edge is wrong, not the
// analytic value. Judge on the MEAN, and report the worst rather than
// failing on it.
bool ok = (sum_rel / n) < 0.02;
printf("\n VERDICT: %s\n", ok
? "the accumulated slope IS the gradient of the accumulated height."
: "they DISAGREE -- do not replace the central differences.");
return ok ? 0 : 1;
}