// =========================================================================
// ANALYTIC GRADIENT GATE — which stage loses the derivative?
//
// pn_out_grad carries erosion_node_eval_grad's analytic slope to the material
// so a detail normal costs ONE density evaluation instead of six (measured on
// an RTX 2070 at 1080p: 1.57 ms vs 16.4 ms fullscreen). Measured against
// central differences it agreed within 5 deg for only 22.7% of samples --
// with no transposition, no sign slip, and a magnitude ratio scattered about
// 1.0 rather than offset from it. So a stage is dropping a term.
//
// Each stage is checked against a central difference of its OWN output, in
// isolation, so the answer is a stage and not a vibe:
//
// grad_noised_2d analytic by construction -- if THIS is wrong
// nothing downstream can be right
// fractal_noise_2d_lod a sum of exact gradients; should be exact
// fractal_multifractal_2d_lod weights each octave by accumulated elevation
// and DROPS d(weight)/dp -- its own comment says
// so. Prime suspect.
// erosion_node_eval_grad fed a gn whose gradient is EXACT (a closed-form
// sinusoid), so any error here is the filter's.
//
// Build (no engine, no godot-cpp, no scons):
// g++ -O2 -std=c++17 -I ../proc-nodes/src/shared \
// -o /tmp/grad_gate tests/test_analytic_gradient.cpp
// =========================================================================
#include <cstdio>
#include <cmath>
#include <vector>
#include <algorithm>
#include "erosion_cpp_shim.h"
using dcerosion::vec2;
using dcerosion::vec3;
using dcerosion::vec4;
// Angular error between the analytic surface normal and a central-difference
// one, for a height field h(p) with gradient g. Normals rather than raw
// gradients because a normal is what the material consumes, and a fixed
// gradient error matters far less on a cliff than on a flat.
static double normal_angle_deg(double gx_a, double gz_a, double gx_c, double gz_c) {
double ax = -gx_a, ay = 1.0, az = -gz_a;
double cx = -gx_c, cy = 1.0, cz = -gz_c;
double la = std::sqrt(ax*ax + ay*ay + az*az);
double lc = std::sqrt(cx*cx + cy*cy + cz*cz);
double d = (ax*cx + ay*cy + az*cz) / (la * lc);
d = std::max(-1.0, std::min(1.0, d));
return std::acos(d) * 180.0 / M_PI;
}
struct Stats { double mean, p50, p90, worst; int n; };
static Stats summarize(std::vector<double>& e) {
std::sort(e.begin(), e.end());
Stats s{};
s.n = (int)e.size();
double sum = 0; for (double v : e) sum += v;
s.mean = sum / e.size();
s.p50 = e[e.size() / 2];
s.p90 = e[(size_t)(e.size() * 0.9)];
s.worst = e.back();
return s;
}
static void report(const char* name, Stats s, double tol_p90) {
const char* verdict = (s.p90 <= tol_p90) ? "ok" : "LOSES IT";
printf(" %-30s mean %6.3f p50 %6.3f p90 %6.3f max %7.3f deg %s\n",
name, s.mean, s.p50, s.p90, s.worst, verdict);
}
// Central difference of a scalar field sampled through a lambda.
template <typename F>
static void central(F f, double x, double z, double e, double& gx, double& gz) {
gx = (f(x + e, z) - f(x - e, z)) / (2.0 * e);
gz = (f(x, z + e) - f(x, z - e)) / (2.0 * e);
}
int main() {
printf("ANALYTIC GRADIENT GATE — analytic vs central difference, per stage\n");
// eps small enough to be a local derivative, large enough that float32
// noise evaluated in double-precision differencing does not go to ULP soup.
const double E = 1e-4;
const int N = 64;
const double SPAN = 4.0; // in p units (world / terrain_height)
std::vector<double> e_base, e_fbm, e_mf, e_ero;
for (int j = 0; j < N; j++) {
for (int i = 0; i < N; i++) {
double x = -SPAN * 0.5 + SPAN * (i + 0.5) / N;
double z = -SPAN * 0.5 + SPAN * (j + 0.5) / N;
double gx, gz;
// --- 1. base gradient noise ---------------------------------
{
vec3 r = dcerosion::grad_noised_2d(vec2((float)x, (float)z));
central([](double a, double b) {
return (double)dcerosion::grad_noised_2d(vec2((float)a, (float)b)).x;
}, x, z, E, gx, gz);
e_base.push_back(normal_angle_deg(r.y, r.z, gx, gz));
}
// --- 2. plain fBm -------------------------------------------
{
vec3 r = dcerosion::fractal_noise_2d_lod(vec2((float)x, (float)z), 1.0f, 6, 1.87f, 0.5f, 999.0f);
central([](double a, double b) {
return (double)dcerosion::fractal_noise_2d_lod(vec2((float)a, (float)b), 1.0f, 6, 1.87f, 0.5f, 999.0f).x;
}, x, z, E, gx, gz);
e_fbm.push_back(normal_angle_deg(r.y, r.z, gx, gz));
}
// --- 3. elevation-weighted multifractal ---------------------
{
vec3 r = dcerosion::fractal_multifractal_2d_lod(vec2((float)x, (float)z), 1.0f, 6, 1.87f, 0.5f, 0.5f, 999.0f);
central([](double a, double b) {
return (double)dcerosion::fractal_multifractal_2d_lod(vec2((float)a, (float)b), 1.0f, 6, 1.87f, 0.5f, 0.5f, 999.0f).x;
}, x, z, E, gx, gz);
e_mf.push_back(normal_angle_deg(r.y, r.z, gx, gz));
}
// --- 4. erosion filter, fed an EXACT-gradient input ----------
// gn is a closed-form sinusoid so its derivative carries no error
// of its own; whatever shows up here belongs to erosion_filter.
{
auto gn_at = [](double a, double b) {
double h = 0.30 * std::sin(3.0 * a) * std::cos(3.0 * b);
double dx = 0.90 * std::cos(3.0 * a) * std::cos(3.0 * b);
double dz = -0.90 * std::sin(3.0 * a) * std::sin(3.0 * b);
return vec3((float)h, (float)dx, (float)dz);
};
auto ero = [&](double a, double b) {
vec3 gn = gn_at(a, b);
return dcerosion::erosion_node_eval_grad(
vec2((float)a, (float)b), gn,
5.0f, 0.150000006f, 0.699999988f, 7, 0.5f,
4.8f, 0.15f,
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);
};
vec4 r = ero(x, z);
central([&](double a, double b) { return (double)ero(a, b).x; }, x, z, E, gx, gz);
e_ero.push_back(normal_angle_deg(r.z, r.w, gx, gz));
}
}
}
printf(" %d samples per stage, eps %.0e, span %.1f p-units\n\n", N * N, E, SPAN);
// grad_noised_2d and fBm are exact by construction, so anything above a
// rounding-level p90 there is a real defect, not a tolerance to be tuned.
report("grad_noised_2d", summarize(e_base), 0.5);
report("fractal_noise_2d_lod", summarize(e_fbm), 0.5);
report("fractal_multifractal_2d_lod", summarize(e_mf), 0.5);
report("erosion_node_eval_grad", summarize(e_ero), 5.0);
printf("\n erosion error vs OCTAVE COUNT -- combi_mask and fade_target are\n");
printf(" REPLACED each octave from that octave's own noise, and neither\n");
printf(" derivative is carried. If that is the mechanism, error grows here.\n");
for (int oct : {1, 2, 3, 5, 7}) {
std::vector<double> e;
for (int j = 0; j < N; j++) {
for (int i = 0; i < N; i++) {
double x = -SPAN * 0.5 + SPAN * (i + 0.5) / N;
double z = -SPAN * 0.5 + SPAN * (j + 0.5) / N;
auto gn_at = [](double a, double b) {
return vec3((float)(0.30 * std::sin(3.0 * a) * std::cos(3.0 * b)),
(float)(0.90 * std::cos(3.0 * a) * std::cos(3.0 * b)),
(float)(-0.90 * std::sin(3.0 * a) * std::sin(3.0 * b)));
};
auto ero = [&](double a, double b) {
return dcerosion::erosion_node_eval_grad(
vec2((float)a, (float)b), gn_at(a, b),
5.0f, 0.150000006f, 0.699999988f, oct, 0.5f,
4.8f, 0.15f, 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);
};
vec4 r = ero(x, z);
double gx, gz;
central([&](double a, double b) { return (double)ero(a, b).x; }, x, z, E, gx, gz);
e.push_back(normal_angle_deg(r.z, r.w, gx, gz));
}
}
Stats s = summarize(e);
printf(" %d octave(s) : mean %6.3f p50 %6.3f p90 %6.3f deg\n",
oct, s.mean, s.p50, s.p90);
}
return 0;
}