// ============================================================================
// EROSION RENDERER — the map's surface as a heightfield, on the GPU, from the
// same erosion_filter.inc and the same wrapper text the parity oracle guards.
//
// WHY IT EXISTS
//
// The CPU dumper takes 13.5 minutes for a 512^2 tile. That is fine for a
// measurement and useless for looking at something while turning a knob. This
// does the same job on the GPU.
//
// WHAT IT SHARES, AND WHAT IT DOES NOT
//
// SHARED, and this is the point:
// - erosion_filter.inc, read off disk at run time and compiled by glslc.
// - the map WRAPPER, from proc-nodes/src/shared/erosion_wrapper_glsl.h --
// the same function test_erosion_parity.cpp calls, which is why adding
// this second GPU consumer did not create a third copy of the wrapper.
// - the Vulkan setup, by including the oracle's translation unit.
//
// NOT SHARED:
// - the root finder. find_surface() is C++ in the oracle; the scan-and-
// bisect below is GLSL, because doing it per column inside the shader is
// the whole speed win. Two implementations of a generic root finder is a
// real duplication, but a much smaller one than the wrapper, and it
// fails LOUDLY: --verify checks this renderer against the oracle's own
// published surface heights, and those came from the C++ finder.
//
// VERIFY BEFORE BELIEVING A PICTURE (#45 section 5 -- instruments are layers)
//
// erosion_render --verify
//
// evaluates the nine columns the oracle prints and compares. The oracle's
// values at MAP defaults are, from a green run:
// (-10000,-10000) 5430.60 (0,0) -1059.81 (10000,10000) -5038.45
// A picture from a renderer that has not passed --verify is decoration.
//
// Build (no engine, no godot-cpp, no scons, nothing into bin/):
// cd dc-terrain
// g++ -O2 -std=c++17 -I ../proc-nodes/src/shared -I . \
// -o /tmp/erosion_render ../terrain-lookdev/erosion_render.cpp -lvulkan
//
// Run (from the dc-terrain root, so the default --inc path resolves):
// /tmp/erosion_render --n 512 --mpp 30 --out /tmp/h.f32
// /tmp/erosion_render --verify
//
// Output: N*N float32 LE surface heights in metres, row-major (z outer, x
// inner). Columns with no surface are -3.0e38.
// ============================================================================
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <cmath>
#include <string>
#include <vector>
#define main parity_main_unused
#include "tests/test_erosion_parity.cpp"
#undef main
static const float NO_SURFACE = -3.0e38f;
// ---------------------------------------------------------------------------
// EXTERNAL HEIGHTFIELD (.hblk), written by terrain-lookdev/stitch_block.py.
// Row = world Z, column = world X -- the convention erode_heightmap.cpp uses on
// its raw input, so the CPU reference and this renderer address the field
// identically rather than by coincidence.
struct HeightBlock {
int n = 0;
float mpp = 30.0f;
float ox = 0.0f, oz = 0.0f;
std::vector<float> h;
};
static HeightBlock read_hblk(const char* path) {
FILE* f = fopen(path, "rb");
if (!f) { fprintf(stderr, "FATAL: cannot open %s\n", path); exit(1); }
char magic[4]; uint32_t ver, n, r0, r1; float mpp, ox, oz;
if (fread(magic,1,4,f)!=4 || memcmp(magic,"HBLK",4)!=0) {
fprintf(stderr, "FATAL: %s is not an .hblk\n", path); exit(1); }
if (fread(&ver,4,1,f)!=1 || fread(&n,4,1,f)!=1 || fread(&mpp,4,1,f)!=1 ||
fread(&ox,4,1,f)!=1 || fread(&oz,4,1,f)!=1 ||
fread(&r0,4,1,f)!=1 || fread(&r1,4,1,f)!=1) {
fprintf(stderr, "FATAL: %s: short header\n", path); exit(1); }
if (ver != 1) { fprintf(stderr, "FATAL: %s: version %u\n", path, ver); exit(1); }
HeightBlock b; b.n=(int)n; b.mpp=mpp; b.ox=ox; b.oz=oz;
std::vector<int16_t> raw((size_t)n*n);
if (fread(raw.data(), sizeof(int16_t), raw.size(), f) != raw.size()) {
fprintf(stderr, "FATAL: %s: short data\n", path); exit(1); }
fclose(f);
b.h.resize(raw.size());
for (size_t i = 0; i < raw.size(); i++) b.h[i] = (float)raw[i];
return b;
}
// THE STORAGE the shared sampler reads -- and nothing else.
//
// The split of duties moved. The consumer declares STORAGE; the wrapper
// declares the SAMPLER: erosion_heightmap_sampler_glsl() emits hblk_texel_l()
// and hblk_sample(), and heightmap_gn calls those, not any function of ours.
// So the old flat bilinear hblk_height_m is gone rather than fixed -- nothing
// called it. (erosion_wrapper_glsl.h:115 still documents the pre-mip contract
// "the consumer must define float hblk_height_m(...)" and is stale.) Its body
// also still indexed hblk_data[] after the buffer was renamed to custom_data,
// which is what produced the undeclared-identifier errors that found this.
//
// THE WHOLE PYRAMID, not level 0. custom_data holds every mip level end to
// end and HBLK_OFF/HBLK_SIDE index into it. Uploading only level 0 is exactly
// #55 section 2.1: HBLK_OFF[1] then points one float past the end of the
// buffer, every sample coarser than the block's pitch reads out of bounds and
// comes back zero, and the terrain is flat wherever the LOD is coarse.
//
// voxel_scale is the sampler's LOD selector and has no declaration of its own
// anywhere: the wrapper emits M_VOXEL, a different name, and emits it AFTER
// this text, so it cannot be aliased here. It comes from the same
// wp.voxel_size that feeds M_VOXEL, so the two cannot disagree.
static std::string hblk_sampler_glsl(const HeightBlock& b,
const HeightBlockMips& m,
float voxel_size) {
char buf[1024];
snprintf(buf, sizeof(buf),
"layout(set=0, binding=1) readonly buffer HBlk { float custom_data[]; };\n"
"const int HBLK_N = %d;\n"
"const float HBLK_OX = %.9g;\n"
"const float HBLK_OZ = %.9g;\n"
"const float HBLK_MPP = %.9g;\n"
"const float voxel_scale = %.9g;\n",
b.n, (double)b.ox, (double)b.oz, (double)b.mpp, (double)voxel_size);
return std::string(buf) + erosion_heightmap_levels_glsl(m);
}
struct PushC {
uint32_t w, h;
float ox, oz, mpp;
uint32_t row0;
};
// The root finder, in GLSL. Deliberately the same shape as find_surface() in
// test_erosion_parity.cpp: a top-down scan for a sign change, then bisection
// inside the bracket that found it. Same STEPS and same bisection count, so
// --verify is comparing the FIELD and not two different search policies.
static std::string build_render_shader(const std::string& inc,
const ErosionWrapperParams& wp,
int steps, int bisect,
const HeightBlock* hb = nullptr,
const HeightBlockMips* hbm = nullptr) {
std::string s =
"#version 450\n"
"layout(local_size_x = 8, local_size_y = 8) in;\n"
"layout(set=0, binding=0) writeonly buffer OutH { float heights[]; };\n"
"layout(push_constant) uniform PC {\n"
" uint w; uint h; float ox; float oz; float mpp; uint row0;\n"
"} pc;\n"
"\n";
s += inc;
s += "\n";
// Must precede the wrapper: the sampler the wrapper emits reads these.
if (hb) s += hblk_sampler_glsl(*hb, *hbm, wp.voxel_size);
s += erosion_wrapper_glsl(wp);
char tail[2048];
snprintf(tail, sizeof(tail),
"const int SCAN_STEPS = %d;\n"
"const int BISECT = %d;\n"
"void main() {\n"
" uint gx = gl_GlobalInvocationID.x;\n"
" uint gy = gl_GlobalInvocationID.y;\n"
" if (gx >= pc.w || gy >= pc.h) return;\n"
" float x = pc.ox + float(gx) * pc.mpp;\n"
" float z = pc.oz + float(pc.row0 + gy) * pc.mpp;\n"
" float lo = M_FLOOR_Y + 1.0;\n"
" float hi = M_CEIL_Y - 1.0;\n"
" float prev_y = hi;\n"
" float prev_d = gpu_wrapper(x, hi, z);\n"
" float found = %.9g;\n"
" for (int i = 1; i <= SCAN_STEPS; i++) {\n"
" float y = hi + (lo - hi) * float(i) / float(SCAN_STEPS);\n"
" float d = gpu_wrapper(x, y, z);\n"
" if ((prev_d < 0.0) != (d < 0.0)) {\n"
" float a = y, b = prev_y;\n"
" for (int k = 0; k < BISECT; k++) {\n"
" float m = 0.5 * (a + b);\n"
" if ((gpu_wrapper(x, m, z) < 0.0) == (prev_d < 0.0)) b = m;\n"
" else a = m;\n"
" }\n"
" found = 0.5 * (a + b);\n"
" break;\n"
" }\n"
" prev_y = y; prev_d = d;\n"
" }\n"
" heights[gy * pc.w + gx] = found;\n"
"}\n",
steps, bisect, (double)NO_SURFACE);
s += tail;
return s;
}
// ---------------------------------------------------------------------------
struct Renderer {
Ctx c{};
VkShaderModule mod = VK_NULL_HANDLE;
VkDescriptorSetLayout dsl = VK_NULL_HANDLE;
VkPipelineLayout playout = VK_NULL_HANDLE;
VkPipeline pipe = VK_NULL_HANDLE;
VkDescriptorPool dpool = VK_NULL_HANDLE;
VkDescriptorSet dset = VK_NULL_HANDLE;
std::vector<float> hblk; // the heightfield, empty when not in use
void build(const std::string& src) {
init_vulkan(c);
std::vector<uint32_t> spv = compile_glsl(src, "erosion render");
VkShaderModuleCreateInfo smi{VK_STRUCTURE_TYPE_SHADER_MODULE_CREATE_INFO};
smi.codeSize = spv.size() * 4; smi.pCode = spv.data();
VK(vkCreateShaderModule(c.dev, &smi, nullptr, &mod));
VkDescriptorSetLayoutBinding bind[2]{};
for (int k = 0; k < 2; k++) {
bind[k].binding = (uint32_t)k;
bind[k].descriptorType = VK_DESCRIPTOR_TYPE_STORAGE_BUFFER;
bind[k].descriptorCount = 1;
bind[k].stageFlags = VK_SHADER_STAGE_COMPUTE_BIT;
}
VkDescriptorSetLayoutCreateInfo dli{VK_STRUCTURE_TYPE_DESCRIPTOR_SET_LAYOUT_CREATE_INFO};
dli.bindingCount = 2; dli.pBindings = bind;
VK(vkCreateDescriptorSetLayout(c.dev, &dli, nullptr, &dsl));
VkPushConstantRange pcr{VK_SHADER_STAGE_COMPUTE_BIT, 0, sizeof(PushC)};
VkPipelineLayoutCreateInfo pli{VK_STRUCTURE_TYPE_PIPELINE_LAYOUT_CREATE_INFO};
pli.setLayoutCount = 1; pli.pSetLayouts = &dsl;
pli.pushConstantRangeCount = 1; pli.pPushConstantRanges = &pcr;
VK(vkCreatePipelineLayout(c.dev, &pli, nullptr, &playout));
VkComputePipelineCreateInfo cpi{VK_STRUCTURE_TYPE_COMPUTE_PIPELINE_CREATE_INFO};
cpi.stage.sType = VK_STRUCTURE_TYPE_PIPELINE_SHADER_STAGE_CREATE_INFO;
cpi.stage.stage = VK_SHADER_STAGE_COMPUTE_BIT;
cpi.stage.module = mod; cpi.stage.pName = "main";
cpi.layout = playout;
VK(vkCreateComputePipelines(c.dev, VK_NULL_HANDLE, 1, &cpi, nullptr, &pipe));
VkDescriptorPoolSize ps{VK_DESCRIPTOR_TYPE_STORAGE_BUFFER, 2};
VkDescriptorPoolCreateInfo dpi{VK_STRUCTURE_TYPE_DESCRIPTOR_POOL_CREATE_INFO};
dpi.maxSets = 1; dpi.poolSizeCount = 1; dpi.pPoolSizes = &ps;
VK(vkCreateDescriptorPool(c.dev, &dpi, nullptr, &dpool));
VkDescriptorSetAllocateInfo dai{VK_STRUCTURE_TYPE_DESCRIPTOR_SET_ALLOCATE_INFO};
dai.descriptorPool = dpool; dai.descriptorSetCount = 1; dai.pSetLayouts = &dsl;
VK(vkAllocateDescriptorSets(c.dev, &dai, &dset));
}
// One band of rows per submit. A single dispatch over a whole tile is
// hundreds of millions of erosion evaluations, and a compute submit that
// long on the DISPLAY gpu trips the driver's watchdog and takes the
// desktop with it. Banding keeps each submit short.
void render(float ox, float oz, uint32_t w, uint32_t h, float mpp,
uint32_t band, std::vector<float>& out) {
out.assign((size_t)w * h, NO_SURFACE);
Buf buf = make_host_buf(c, (size_t)w * band * sizeof(float));
// Binding 1 always exists in the layout; with no heightfield it is a
// one-float placeholder the shader never declares. Leaving it UNWRITTEN
// would be undefined behaviour rather than merely unused.
size_t hb_bytes = hblk.empty() ? sizeof(float) : hblk.size()*sizeof(float);
Buf hbuf = make_host_buf(c, hb_bytes);
if (!hblk.empty()) memcpy(hbuf.mapped, hblk.data(), hb_bytes);
else ((float*)hbuf.mapped)[0] = 0.0f;
VkDescriptorBufferInfo bi[2] = {{buf.buf, 0, VK_WHOLE_SIZE},
{hbuf.buf, 0, VK_WHOLE_SIZE}};
VkWriteDescriptorSet wr[2]{};
for (int k = 0; k < 2; k++) {
wr[k].sType = VK_STRUCTURE_TYPE_WRITE_DESCRIPTOR_SET;
wr[k].dstSet = dset; wr[k].dstBinding = (uint32_t)k;
wr[k].descriptorCount = 1;
wr[k].descriptorType = VK_DESCRIPTOR_TYPE_STORAGE_BUFFER;
wr[k].pBufferInfo = &bi[k];
}
vkUpdateDescriptorSets(c.dev, 2, wr, 0, nullptr);
for (uint32_t row0 = 0; row0 < h; row0 += band) {
uint32_t rows = (row0 + band <= h) ? band : (h - row0);
PushC pc{w, rows, ox, oz, mpp, row0};
VkCommandBufferAllocateInfo cbi{VK_STRUCTURE_TYPE_COMMAND_BUFFER_ALLOCATE_INFO};
cbi.commandPool = c.cmdpool; cbi.level = VK_COMMAND_BUFFER_LEVEL_PRIMARY;
cbi.commandBufferCount = 1;
VkCommandBuffer cmd;
VK(vkAllocateCommandBuffers(c.dev, &cbi, &cmd));
VkCommandBufferBeginInfo cbb{VK_STRUCTURE_TYPE_COMMAND_BUFFER_BEGIN_INFO};
cbb.flags = VK_COMMAND_BUFFER_USAGE_ONE_TIME_SUBMIT_BIT;
VK(vkBeginCommandBuffer(cmd, &cbb));
vkCmdBindPipeline(cmd, VK_PIPELINE_BIND_POINT_COMPUTE, pipe);
vkCmdBindDescriptorSets(cmd, VK_PIPELINE_BIND_POINT_COMPUTE, playout,
0, 1, &dset, 0, nullptr);
vkCmdPushConstants(cmd, playout, VK_SHADER_STAGE_COMPUTE_BIT, 0,
sizeof(PushC), &pc);
vkCmdDispatch(cmd, (w + 7) / 8, (rows + 7) / 8, 1);
VK(vkEndCommandBuffer(cmd));
VkSubmitInfo si{VK_STRUCTURE_TYPE_SUBMIT_INFO};
si.commandBufferCount = 1; si.pCommandBuffers = &cmd;
VK(vkQueueSubmit(c.queue, 1, &si, VK_NULL_HANDLE));
VK(vkQueueWaitIdle(c.queue));
vkFreeCommandBuffers(c.dev, c.cmdpool, 1, &cmd);
memcpy(out.data() + (size_t)row0 * w, buf.mapped,
(size_t)rows * w * sizeof(float));
}
}
};
// ---------------------------------------------------------------------------
static void set_param(ErosionWrapperParams& p, const char* k, const char* v) {
float f = (float)atof(v);
int i = atoi(v);
if (!strcmp(k, "floor_y")) p.floor_y = f;
else if (!strcmp(k, "ceiling_y")) p.ceiling_y = f;
else if (!strcmp(k, "terrain_height")) p.terrain_height = f;
else if (!strcmp(k, "mountain_scale")) p.mountain_scale = f;
else if (!strcmp(k, "mountain_height")) p.mountain_height = f;
else if (!strcmp(k, "mountain_gain")) p.mountain_gain = f;
else if (!strcmp(k, "mountain_variance")) p.mountain_variance = f;
else if (!strcmp(k, "mountain_octaves")) p.mountain_octaves = i;
else if (!strcmp(k, "relief_scale")) p.relief_scale = f;
else if (!strcmp(k, "relief_floor")) p.relief_floor = f;
else if (!strcmp(k, "detail_scale")) p.detail_scale = f;
else if (!strcmp(k, "detail_height")) p.detail_height = f;
else if (!strcmp(k, "detail_gain")) p.detail_gain = f;
else if (!strcmp(k, "detail_octaves")) p.detail_octaves = i;
else if (!strcmp(k, "noise_lacunarity")) p.noise_lacunarity = f;
else if (!strcmp(k, "voxel_size")) p.voxel_size = f;
else if (!strcmp(k, "gn_scale")) p.gn_scale = f;
else if (!strcmp(k, "gn_gain")) p.gn_gain = f;
else if (!strcmp(k, "gn_detail")) p.gn_detail = i;
else if (!strcmp(k, "en_height_scale")) p.en_height_scale = f;
else if (!strcmp(k, "en_strength")) p.en_strength = f;
else if (!strcmp(k, "en_cell_scale")) p.en_cell_scale = f;
else if (!strcmp(k, "en_gain")) p.en_gain = f;
else if (!strcmp(k, "en_octaves")) p.en_octaves = i;
else if (!strcmp(k, "en_erosion_scale")) p.en_erosion_scale = f;
else if (!strcmp(k, "en_fade_div")) p.en_fade_div = f;
else if (!strcmp(k, "en_height_amp")) p.en_height_amp = f;
else if (!strcmp(k, "en_erosion_lac")) p.en_erosion_lac = f;
else if (!strcmp(k, "en_gully_weight")) p.en_gully_weight = f;
else if (!strcmp(k, "en_detail")) p.en_detail = f;
else if (!strcmp(k, "en_normalization")) p.en_normalization = f;
else if (!strcmp(k, "en_height_offset_a")) p.en_height_offset_a = f;
else if (!strcmp(k, "en_height_offset_b")) p.en_height_offset_b = f;
else if (!strcmp(k, "rounding_ridges")) p.rounding_ridges = f;
else if (!strcmp(k, "rounding_creases")) p.rounding_creases = f;
else if (!strcmp(k, "rounding_input")) p.rounding_input = f;
else if (!strcmp(k, "rounding_per_oct")) p.rounding_per_oct = f;
else if (!strcmp(k, "onset_input")) p.onset_input = f;
else if (!strcmp(k, "onset_per_oct")) p.onset_per_oct = f;
else if (!strcmp(k, "onset_ridge_input")) p.onset_ridge_input = f;
else if (!strcmp(k, "onset_ridge_oct")) p.onset_ridge_oct = f;
else if (!strcmp(k, "assumed_slope_mag")) p.assumed_slope_mag = f;
else if (!strcmp(k, "assumed_slope_mix")) p.assumed_slope_mix = f;
else if (!strcmp(k, "field_half_range")) p.field_half_range = f;
else { fprintf(stderr, "unknown parameter '%s'\n", k); exit(2); }
}
int main(int argc, char** argv) {
const char* inc_path = "../proc-nodes/src/shared/erosion_filter.inc";
const char* out_path = "/tmp/h.f32";
int n = 512, steps = 512, bisect = 60, band = 32;
float mpp = 30.0f, ox = 0.0f, oz = 0.0f;
bool verify = false;
const char* hblk_path = nullptr;
bool have_geom = false;
// NOT ErosionWrapperParams{} -- the struct's own defaults are a convenience
// and are not authoritative. The map's values come from the oracle's single
// mapping, so a constant that moves there moves here in the same commit.
// Hand-copying them instead is what put this renderer 7.3 km out on its
// first run, and the --verify gate is the only reason that was noticed.
ErosionWrapperParams wp = erosion_wrapper_params_from_map();
for (int i = 1; i < argc; i++) {
const char* a = argv[i];
if (!strcmp(a, "--inc") && i+1 < argc) inc_path = argv[++i];
else if (!strcmp(a, "--out") && i+1 < argc) out_path = argv[++i];
else if (!strcmp(a, "--n") && i+1 < argc) { n = atoi(argv[++i]); have_geom = true; }
else if (!strcmp(a, "--mpp") && i+1 < argc) { mpp = (float)atof(argv[++i]); have_geom = true; }
else if (!strcmp(a, "--ox") && i+1 < argc) { ox = (float)atof(argv[++i]); have_geom = true; }
else if (!strcmp(a, "--oz") && i+1 < argc) { oz = (float)atof(argv[++i]); have_geom = true; }
else if (!strcmp(a, "--heightmap") && i+1 < argc) hblk_path = argv[++i];
else if (!strcmp(a, "--steps") && i+1 < argc) steps = atoi(argv[++i]);
else if (!strcmp(a, "--bisect") && i+1 < argc) bisect = atoi(argv[++i]);
else if (!strcmp(a, "--band") && i+1 < argc) band = atoi(argv[++i]);
else if (!strcmp(a, "--verify")) verify = true;
else if (!strcmp(a, "--set") && i+2 < argc) { set_param(wp, argv[i+1], argv[i+2]); i += 2; }
else { fprintf(stderr, "unknown arg %s\n", a); return 2; }
}
std::string inc = read_file(inc_path);
HeightBlock hb;
HeightBlockMips hbm;
bool use_hblk = (hblk_path != nullptr);
if (use_hblk) {
hb = read_hblk(hblk_path);
wp.heightmap_source = true;
wp.heightmap_step = hb.mpp;
if (!have_geom) {
n = hb.n - 2; // render the interior; the ring is apron
mpp = hb.mpp;
ox = hb.ox + hb.mpp;
oz = hb.oz + hb.mpp;
}
hbm = erosion_heightmap_build_mips(hb.h.data(), hb.n);
printf("heightfield: %s\n", hblk_path);
printf(" %d x %d px at %g m/px, world x %.0f z %.0f\n",
hb.n, hb.n, hb.mpp, hb.ox, hb.oz);
}
std::string src = build_render_shader(inc, wp, steps, bisect,
use_hblk ? &hb : nullptr,
use_hblk ? &hbm : nullptr);
Renderer r;
if (use_hblk) r.hblk = hbm.data;
r.build(src);
if (verify) {
// The oracle's own columns, at MAP defaults, from a green run.
struct { float x, z, expect; } cols[] = {
{-10000.0f, -10000.0f, 5430.60f},
{-10000.0f, 0.0f, 4179.92f},
{-10000.0f, 10000.0f, -353.49f},
{ 0.0f, -10000.0f, 4299.36f},
{ 0.0f, 0.0f, -1059.81f},
{ 0.0f, 10000.0f, -5454.78f},
{ 10000.0f, -10000.0f, 215.53f},
{ 10000.0f, 0.0f, -4136.03f},
{ 10000.0f, 10000.0f, -5038.45f},
};
printf("VERIFY against test_erosion_parity's published columns\n");
printf(" scan steps %d, bisect %d (the oracle uses 512 / 60)\n\n", steps, bisect);
printf(" %11s %11s %13s %13s %12s\n", "x", "z", "oracle", "renderer", "delta");
double worst = 0.0;
for (auto& col : cols) {
std::vector<float> one;
r.render(col.x, col.z, 1, 1, 1.0f, 1, one);
double d = (double)one[0] - (double)col.expect;
if (fabs(d) > worst) worst = fabs(d);
printf(" %11.1f %11.1f %13.2f %13.2f %12.4f\n",
col.x, col.z, col.expect, one[0], d);
}
printf("\n worst |delta| : %.4f m\n", worst);
// The oracle prints 2 decimals, so agreement can only be shown to
// about half a centimetre. Anything above that is a real difference.
bool ok = worst < 0.01;
printf(" VERDICT: %s\n", ok ? "renderer MATCHES the oracle."
: "renderer DISAGREES with the oracle. Do not "
"trust any picture it makes.");
return ok ? 0 : 1;
}
printf("erosion render: %dx%d px @ %g m/px = %.2f km across\n",
n, n, mpp, n * mpp / 1000.0f);
printf(" origin (%g, %g) m voxel %g scan %d/%d band %d rows\n",
ox, oz, wp.voxel_size, steps, bisect, band);
std::vector<float> h;
r.render(ox, oz, (uint32_t)n, (uint32_t)n, mpp, (uint32_t)band, h);
long hits = 0; double lo = 1e30, hi = -1e30, sum = 0;
for (float v : h) {
if (v == NO_SURFACE) continue;
hits++; sum += v;
if (v < lo) lo = v;
if (v > hi) hi = v;
}
printf(" surface found in %ld of %d columns\n", hits, n * n);
if (hits) printf(" height min %.1f m max %.1f m mean %.1f m\n", lo, hi, sum / hits);
FILE* f = fopen(out_path, "wb");
if (!f) { perror("fopen"); return 1; }
fwrite(h.data(), sizeof(float), h.size(), f);
fclose(f);
printf(" wrote %zu floats to %s\n", h.size(), out_path);
return hits > 0 ? 0 : 1;
}