| … | … | @@ -0,0 +1,220 @@ |
| 1 | // ============================================================================ |
| 2 | // ERODE A HEIGHTMAP — run runevision's filter on ANY input height field, in |
| 3 | // particular a terrain-diffusion tile. |
| 4 | // |
| 5 | // WHY |
| 6 | // |
| 7 | // #46 §4 measured that the diffusion model produces landscape (coherent |
| 8 | // catchments, beta steepening 4.18 -> 5.10) where our fBm mountain layer |
| 9 | // produces texture (flat beta). The blog is explicit that the filter "can be |
| 10 | // applied on top of any height function". So: let the model supply the large |
| 11 | // forms it is demonstrably good at, and let the filter supply detail below its |
| 12 | // 30 m floor. This is #46 §4b option 1, made runnable. |
| 13 | // |
| 14 | // It needs no GPU descriptor work. A heightmap transform is ONE erosion |
| 15 | // evaluation per pixel with no surface search, so the CPU shim is enough -- |
| 16 | // unlike dump_erosion_heightfield, which pays 572 density evaluations per |
| 17 | // column bisecting for a surface. |
| 18 | // |
| 19 | // UNITS, which is the whole difficulty |
| 20 | // |
| 21 | // The filter is dimensionless. Its p is world position over a horizontal unit |
| 22 | // L, so the coarsest gully wavelength is L * erosion_scale * cell_scale metres. |
| 23 | // Its height is normalised by a vertical unit H, and its SLOPE must be in the |
| 24 | // same space: dh/dp = (dh/dx_world) * L / H. Get that factor wrong and the |
| 25 | // gullies come out oriented correctly but scaled absurdly, which is exactly the |
| 26 | // 60x wrapper bug #45 found. |
| 27 | // |
| 28 | // THE IDENTITY GATE |
| 29 | // |
| 30 | // With --strength 0 the filter contributes nothing, so the output must equal |
| 31 | // the input to float precision. --check runs that and refuses to proceed if it |
| 32 | // does not hold. A unit error that survives this cannot exist. |
| 33 | // |
| 34 | // Build: |
| 35 | // cd dc-terrain |
| 36 | // g++ -O2 -std=c++17 -I ../proc-nodes/src/shared -o /tmp/erode_heightmap \ |
| 37 | // ../terrain-lookdev/erode_heightmap.cpp |
| 38 | // |
| 39 | // Run: |
| 40 | // /tmp/erode_heightmap --check |
| 41 | // /tmp/erode_heightmap --in tile.f32 --n 512 --mpp 30 --out eroded.f32 \ |
| 42 | // --gully-metres 900 --strength 0.22 --octaves 5 |
| 43 | // ============================================================================ |
| 44 | #include <cstdio> |
| 45 | #include <cstdlib> |
| 46 | #include <cstring> |
| 47 | #include <cmath> |
| 48 | #include <vector> |
| 49 | #include <string> |
| 50 | |
| 51 | #include "erosion_cpp_shim.h" |
| 52 | |
| 53 | using dcerosion::vec2; |
| 54 | using dcerosion::vec3; |
| 55 | using dcerosion::vec4; |
| 56 | |
| 57 | struct Opts { |
| 58 | int n = 512; |
| 59 | float mpp = 30.0f; |
| 60 | // The coarsest gully, in METRES. Everything else about the horizontal |
| 61 | // layout is derived from it, because a raw "scale" number is unreadable -- |
| 62 | // the same lesson the lookdev panel relabelling just learned. |
| 63 | float gully_metres = 900.0f; |
| 64 | float strength = 0.22f; // runevision's reference |
| 65 | int octaves = 5; // runevision's reference |
| 66 | float gain = 0.5f; |
| 67 | float cell_scale = 0.7f; |
| 68 | float lacunarity = 2.0f; |
| 69 | float gully_weight = 0.5f; |
| 70 | float detail = 1.5f; |
| 71 | float normalization = 0.5f; |
| 72 | float rounding_input = 0.1f; |
| 73 | float height_offset_a = -0.65f; |
| 74 | // Vertical normalisation. The filter's strength is relative to a height of |
| 75 | // order 1, so H should be the input's own relief, not an arbitrary number. |
| 76 | // 0 means "use the tile's (max-min)/2". |
| 77 | float vscale = 0.0f; |
| 78 | }; |
| 79 | |
| 80 | // One evaluation. Kept separate from the loop so --check can call it directly. |
| 81 | static float erode_point(const Opts& o, float h_norm, float dhdx, float dhdz, |
| 82 | float px, float pz, float L, float H, float* ridge) { |
| 83 | // world slope -> filter slope: multiply by L/H, matching p = world/L and |
| 84 | // height = metres/H. |
| 85 | vec3 gn; |
| 86 | gn.x = h_norm; |
| 87 | gn.y = dhdx * (L / H); |
| 88 | gn.z = dhdz * (L / H); |
| 89 | |
| 90 | // height_amp 1 and height_scale 1 leave the filter's own [0,1] remap as the |
| 91 | // only vertical transform, so the inverse below is exact. |
| 92 | vec2 en = dcerosion::erosion_node_eval_full( |
| 93 | vec2(px / L, pz / L), gn, |
| 94 | /* height_scale */ 1.0f, |
| 95 | o.strength, |
| 96 | o.cell_scale, o.octaves, o.gain, |
| 97 | /* fade_div */ 1.0f, |
| 98 | /* erosion_scale */ o.gully_metres / (L * o.cell_scale), |
| 99 | /* height_amp */ 1.0f, |
| 100 | o.lacunarity, o.gully_weight, o.detail, o.normalization, |
| 101 | vec2(o.height_offset_a, 0.0f), |
| 102 | vec4(0.1f, 0.0f, o.rounding_input, 2.0f), |
| 103 | vec4(1.25f, 1.25f, 2.8f, 1.5f), |
| 104 | vec2(0.7f, 1.0f)); |
| 105 | if (ridge) *ridge = en.y; |
| 106 | // invert the remap: the filter returned n.x*0.5+0.5 plus its delta. |
| 107 | return (en.x - 0.5f) * 2.0f; |
| 108 | } |
| 109 | |
| 110 | int main(int argc, char** argv) { |
| 111 | Opts o; |
| 112 | const char* in_path = nullptr; |
| 113 | const char* out_path = "/tmp/eroded.f32"; |
| 114 | const char* ridge_path = nullptr; |
| 115 | bool check = false; |
| 116 | |
| 117 | for (int i = 1; i < argc; i++) { |
| 118 | const char* a = argv[i]; |
| 119 | auto next = [&]() { return argv[++i]; }; |
| 120 | if (!strcmp(a, "--in") && i+1 < argc) in_path = next(); |
| 121 | else if (!strcmp(a, "--out") && i+1 < argc) out_path = next(); |
| 122 | else if (!strcmp(a, "--ridge")&& i+1 < argc) ridge_path = next(); |
| 123 | else if (!strcmp(a, "--n") && i+1 < argc) o.n = atoi(next()); |
| 124 | else if (!strcmp(a, "--mpp") && i+1 < argc) o.mpp = (float)atof(next()); |
| 125 | else if (!strcmp(a, "--gully-metres") && i+1<argc) o.gully_metres = (float)atof(next()); |
| 126 | else if (!strcmp(a, "--strength") && i+1<argc) o.strength = (float)atof(next()); |
| 127 | else if (!strcmp(a, "--octaves") && i+1<argc) o.octaves = atoi(next()); |
| 128 | else if (!strcmp(a, "--gain") && i+1<argc) o.gain = (float)atof(next()); |
| 129 | else if (!strcmp(a, "--cell-scale")&& i+1<argc) o.cell_scale = (float)atof(next()); |
| 130 | else if (!strcmp(a, "--rounding-input") && i+1<argc) o.rounding_input = (float)atof(next()); |
| 131 | else if (!strcmp(a, "--vscale") && i+1<argc) o.vscale = (float)atof(next()); |
| 132 | else if (!strcmp(a, "--offset") && i+1<argc) o.height_offset_a = (float)atof(next()); |
| 133 | else if (!strcmp(a, "--check")) check = true; |
| 134 | else { fprintf(stderr, "unknown arg %s\n", a); return 2; } |
| 135 | } |
| 136 | |
| 137 | // ---- the identity gate ------------------------------------------------- |
| 138 | if (check) { |
| 139 | Opts z = o; z.strength = 0.0f; z.height_offset_a = 0.0f; |
| 140 | double worst = 0.0; |
| 141 | for (int k = 0; k < 200; k++) { |
| 142 | float hn = -1.0f + 2.0f * (float)k / 199.0f; |
| 143 | float g = erode_point(z, hn, 0.013f, -0.007f, |
| 144 | 100.0f * k, 70.0f * k, 6000.0f, 500.0f, nullptr); |
| 145 | double d = fabs((double)g - (double)hn); |
| 146 | if (d > worst) worst = d; |
| 147 | } |
| 148 | printf("IDENTITY GATE — strength 0 and no height offset must be a no-op\n"); |
| 149 | printf(" worst |out - in| over 200 heights : %.3e\n", worst); |
| 150 | bool ok = worst < 1e-5; |
| 151 | printf(" VERDICT: %s\n", ok ? "unit conversion is exact." |
| 152 | : "WRONG. The height mapping does not invert; " |
| 153 | "do not trust any eroded tile."); |
| 154 | if (!ok || !in_path) return ok ? 0 : 1; |
| 155 | } |
| 156 | |
| 157 | if (!in_path) { fprintf(stderr, "need --in\n"); return 2; } |
| 158 | |
| 159 | const int n = o.n; |
| 160 | std::vector<float> h((size_t)n * n); |
| 161 | FILE* f = fopen(in_path, "rb"); |
| 162 | if (!f) { perror("open --in"); return 1; } |
| 163 | if (fread(h.data(), sizeof(float), h.size(), f) != h.size()) { |
| 164 | fprintf(stderr, "short read: expected %zu floats\n", h.size()); return 1; |
| 165 | } |
| 166 | fclose(f); |
| 167 | |
| 168 | float lo = h[0], hi = h[0]; |
| 169 | for (float v : h) { if (v < lo) lo = v; if (v > hi) hi = v; } |
| 170 | const float mid = 0.5f * (lo + hi); |
| 171 | const float H = (o.vscale > 0.0f) ? o.vscale |
| 172 | : fmaxf(0.5f * (hi - lo), 1.0f); |
| 173 | // L is fixed by the gully band we asked for, so the number the user gives |
| 174 | // is in metres and everything else follows. |
| 175 | const float L = o.gully_metres / o.cell_scale; |
| 176 | |
| 177 | printf("erode heightmap: %dx%d @ %g m/px = %.2f km across\n", |
| 178 | n, n, o.mpp, n * o.mpp / 1000.0f); |
| 179 | printf(" input relief %.1f .. %.1f m (mid %.1f, H %.1f)\n", lo, hi, mid, H); |
| 180 | printf(" gullies %.0f m -> %.0f m (%d oct, lac %.2f)\n", |
| 181 | o.gully_metres, o.gully_metres / powf(o.lacunarity, (float)(o.octaves - 1)), |
| 182 | o.octaves, o.lacunarity); |
| 183 | printf(" strength %.3g rounding_input %.3g offset %.3g\n", |
| 184 | o.strength, o.rounding_input, o.height_offset_a); |
| 185 | |
| 186 | std::vector<float> out(h.size()), ridge(h.size()); |
| 187 | for (int j = 0; j < n; j++) { |
| 188 | for (int i = 0; i < n; i++) { |
| 189 | // central differences, clamped at the border |
| 190 | int im = (i > 0) ? i - 1 : i, ip = (i < n - 1) ? i + 1 : i; |
| 191 | int jm = (j > 0) ? j - 1 : j, jp = (j < n - 1) ? j + 1 : j; |
| 192 | float dhdx = (h[(size_t)j * n + ip] - h[(size_t)j * n + im]) / |
| 193 | (o.mpp * (float)(ip - im)); |
| 194 | float dhdz = (h[(size_t)jp * n + i] - h[(size_t)jm * n + i]) / |
| 195 | (o.mpp * (float)(jp - jm)); |
| 196 | float hn = (h[(size_t)j * n + i] - mid) / H; |
| 197 | float r = 0.0f; |
| 198 | float g = erode_point(o, hn, dhdx, dhdz, |
| 199 | (float)i * o.mpp, (float)j * o.mpp, L, H, &r); |
| 200 | out[(size_t)j * n + i] = mid + g * H; |
| 201 | ridge[(size_t)j * n + i] = r; |
| 202 | } |
| 203 | } |
| 204 | |
| 205 | float olo = out[0], ohi = out[0]; |
| 206 | for (float v : out) { if (v < olo) olo = v; if (v > ohi) ohi = v; } |
| 207 | printf(" output relief %.1f .. %.1f m\n", olo, ohi); |
| 208 | |
| 209 | f = fopen(out_path, "wb"); |
| 210 | if (!f) { perror("open --out"); return 1; } |
| 211 | fwrite(out.data(), sizeof(float), out.size(), f); |
| 212 | fclose(f); |
| 213 | printf(" wrote %s\n", out_path); |
| 214 | if (ridge_path) { |
| 215 | f = fopen(ridge_path, "wb"); |
| 216 | if (f) { fwrite(ridge.data(), sizeof(float), ridge.size(), f); fclose(f); |
| 217 | printf(" wrote %s (ridge: -1 creases, +1 ridges)\n", ridge_path); } |
| 218 | } |
| 219 | return 0; |
| 220 | } |