05b8b3a · terrain-lookdev

Changed files (8) +756 -19

.gitignore +1 -0

…	…	@@ -0,0 +1 @@
	1	__pycache__/

compare_engine_vs_lookdev.py +86 -0

…	…	@@ -0,0 +1,86 @@
	1	#!/usr/bin/env python3
	2	"""
	3	Put the ENGINE's terrain next to the LOOKDEV's, same ground, same colours.
	4
	5	Every check in this session has produced a number, and a number has not once
	6	told anybody what they were looking at. This writes one PNG: lookdev on the
	7	left, engine on the right, difference on the far right, hypsometric tint shared
	8	so the two are directly comparable by eye.
	9
	10	The engine half comes from map_dc_diffusion writing /tmp/engine_gpu_surface.f32
	11	on load -- a 64^3 debug_gpu_density dispatch, so it is the GPU's own field, not
	12	the CPU proxy that misled this session repeatedly.
	13
	14	python3 compare_engine_vs_lookdev.py
	15	"""
	16	import os, subprocess, sys
	17	import numpy as np
	18
	19	DC = os.path.abspath(os.path.join(os.path.dirname(__file__), "..", "dc-terrain"))
	20	HBLK = os.path.abspath(os.path.join(os.path.dirname(__file__), "..", "hyper-tps",
	21	"assets", "terrain", "diffusion_gully_30km.hblk"))
	22	ENG = "/tmp/engine_gpu_surface.f32"
	23	SPAN, N = 7680.0, 64
	24	CX, CZ = 53760.0, 38400.0
	25	OX, OZ = CX - SPAN / 2.0, CZ - SPAN / 2.0
	26	STEP = SPAN / float(N - 1)
	27
	28
	29	def tint(h, lo, hi):
	30	t = np.clip((h - lo) / max(hi - lo, 1e-6), 0, 1)
	31	stops = np.array([[0.05,0.12,0.28],[0.20,0.35,0.45],[0.30,0.42,0.22],
	32	[0.55,0.50,0.30],[0.60,0.55,0.52],[0.92,0.92,0.94]])
	33	idx = t * (len(stops) - 1)
	34	i0 = np.floor(idx).astype(int); i1 = np.minimum(i0 + 1, len(stops) - 1)
	35	f = (idx - i0)[..., None]
	36	rgb = stops[i0] * (1 - f) + stops[i1] * f
	37	# hillshade so shape reads, not just altitude
	38	gy, gx = np.gradient(h, STEP)
	39	n = 1.0 / np.sqrt(gxgx + gygy + 1.0)
	40	shade = np.clip(0.45 + 0.75 * ((-gx - gy) * n * 0.5 + 0.5), 0, 1.35)
	41	return np.clip(rgb * shade[..., None], 0, 1)
	42
	43
	44	def main():
	45	if not os.path.exists(ENG):
	46	print(f"missing {ENG} — load map_dc_diffusion once; it writes this on load")
	47	return 2
	48	eng = np.fromfile(ENG, dtype="<f4", count=N*N).astype(np.float64).reshape(N, N)
	49	miss = eng < -99998.0
	50	print(f"engine: {N}x{N}, {(~miss).sum()} columns with a surface, "
	51	f"{eng[~miss].min():.0f}..{eng[~miss].max():.0f} m")
	52
	53	out = "/tmp/lookdev_surface.f32"
	54	r = subprocess.run(["/tmp/erosion_render4", "--heightmap", HBLK,
	55	"--ox", str(OX), "--oz", str(OZ), "--n", str(N),
	56	"--mpp", str(STEP), "--out", out], cwd=DC,
	57	capture_output=True, text=True)
	58	if r.returncode != 0:
	59	print(r.stdout, r.stderr); return 1
	60	look = np.fromfile(out, dtype="<f4", count=N*N).astype(np.float64).reshape(N, N)
	61	print(f"lookdev: {look.min():.0f}..{look.max():.0f} m")
	62
	63	valid = ~miss
	64	d = np.abs(eng - look)
	65	print(f"\n max \|engine - lookdev\| = {d[valid].max():8.1f} m")
	66	print(f" rms = {np.sqrt((d[valid]**2).mean()):8.1f} m")
	67
	68	lo = min(eng[valid].min(), look.min()); hi = max(eng[valid].max(), look.max())
	69	a = tint(look, lo, hi); b = tint(np.where(miss, lo, eng), lo, hi)
	70	dm = np.clip(d / max(d[valid].max(), 1e-6), 0, 1)
	71	c = np.stack([dm, dm*0.3, 1.0-dm], -1)
	72	gap = np.ones((N, 4, 3))
	73	img = (np.concatenate([a, gap, b, gap, c], axis=1) * 255).astype(np.uint8)
	74	try:
	75	from PIL import Image
	76	Image.fromarray(img).resize((img.shape[1]6, img.shape[0]6), Image.NEAREST)\
	77	.save("/tmp/engine_vs_lookdev.png")
	78	print("\n wrote /tmp/engine_vs_lookdev.png (lookdev \| engine \| difference)")
	79	except ImportError:
	80	img.tofile("/tmp/engine_vs_lookdev.raw")
	81	print("\n no PIL; wrote raw RGB to /tmp/engine_vs_lookdev.raw")
	82	return 0
	83
	84
	85	if __name__ == "__main__":
	86	sys.exit(main())

erosion_render.cpp +121 -19

…	…	@@ -59,6 +59,66 @@
59	59
60	60	static const float NO_SURFACE = -3.0e38f;
61	61
	62	// ---------------------------------------------------------------------------
	63	// EXTERNAL HEIGHTFIELD (.hblk), written by terrain-lookdev/stitch_block.py.
	64	// Row = world Z, column = world X -- the convention erode_heightmap.cpp uses on
	65	// its raw input, so the CPU reference and this renderer address the field
	66	// identically rather than by coincidence.
	67	struct HeightBlock {
	68	int n = 0;
	69	float mpp = 30.0f;
	70	float ox = 0.0f, oz = 0.0f;
	71	std::vector<float> h;
	72	};
	73
	74	static HeightBlock read_hblk(const char* path) {
	75	FILE* f = fopen(path, "rb");
	76	if (!f) { fprintf(stderr, "FATAL: cannot open %s\n", path); exit(1); }
	77	char magic[4]; uint32_t ver, n, r0, r1; float mpp, ox, oz;
	78	if (fread(magic,1,4,f)!=4 \|\| memcmp(magic,"HBLK",4)!=0) {
	79	fprintf(stderr, "FATAL: %s is not an .hblk\n", path); exit(1); }
	80	if (fread(&ver,4,1,f)!=1 \|\| fread(&n,4,1,f)!=1 \|\| fread(&mpp,4,1,f)!=1 \|\|
	81	fread(&ox,4,1,f)!=1 \|\| fread(&oz,4,1,f)!=1 \|\|
	82	fread(&r0,4,1,f)!=1 \|\| fread(&r1,4,1,f)!=1) {
	83	fprintf(stderr, "FATAL: %s: short header\n", path); exit(1); }
	84	if (ver != 1) { fprintf(stderr, "FATAL: %s: version %u\n", path, ver); exit(1); }
	85	HeightBlock b; b.n=(int)n; b.mpp=mpp; b.ox=ox; b.oz=oz;
	86	std::vector<int16_t> raw((size_t)n*n);
	87	if (fread(raw.data(), sizeof(int16_t), raw.size(), f) != raw.size()) {
	88	fprintf(stderr, "FATAL: %s: short data\n", path); exit(1); }
	89	fclose(f);
	90	b.h.resize(raw.size());
	91	for (size_t i = 0; i < raw.size(); i++) b.h[i] = (float)raw[i];
	92	return b;
	93	}
	94
	95	// The lookup the shared wrapper's heightmap_gn requires. Bilinear, clamped at
	96	// the edge -- the apron ring means the clamp is only reached OUTSIDE the block's
	97	// footprint, so the central difference inside it always has real neighbours.
	98	static std::string hblk_sampler_glsl(const HeightBlock& b) {
	99	char buf[2048];
	100	snprintf(buf, sizeof(buf),
	101	"layout(set=0, binding=1) readonly buffer HBlk { float custom_data[]; };\n"
	102	"const int HBLK_N = %d;\n"
	103	"const float HBLK_OX = %.9g;\n"
	104	"const float HBLK_OZ = %.9g;\n"
	105	"const float HBLK_MPP = %.9g;\n"
	106	"float hblk_height_m(float x, float z) {\n"
	107	" float u = (x - HBLK_OX) / HBLK_MPP;\n"
	108	" float v = (z - HBLK_OZ) / HBLK_MPP;\n"
	109	" u = clamp(u, 0.0, float(HBLK_N - 1));\n"
	110	" v = clamp(v, 0.0, float(HBLK_N - 1));\n"
	111	" int i0 = int(floor(u)); int j0 = int(floor(v));\n"
	112	" int i1 = min(i0 + 1, HBLK_N - 1); int j1 = min(j0 + 1, HBLK_N - 1);\n"
	113	" float fu = u - float(i0); float fv = v - float(j0);\n"
	114	" float h00 = hblk_data[j0HBLK_N + i0]; float h10 = hblk_data[j0HBLK_N + i1];\n"
	115	" float h01 = hblk_data[j1HBLK_N + i0]; float h11 = hblk_data[j1HBLK_N + i1];\n"
	116	" return mix(mix(h00, h10, fu), mix(h01, h11, fu), fv);\n"
	117	"}\n",
	118	b.n, (double)b.ox, (double)b.oz, (double)b.mpp);
	119	return std::string(buf);
	120	}
	121
62	122	struct PushC {
63	123	uint32_t w, h;
64	124	float ox, oz, mpp;
…	…	@@ -71,7 +131,8 @@ struct PushC {
71	131	// --verify is comparing the FIELD and not two different search policies.
72	132	static std::string build_render_shader(const std::string& inc,
73	133	const ErosionWrapperParams& wp,
74		int steps, int bisect) {
	134	int steps, int bisect,
	135	const HeightBlock* hb = nullptr) {
75	136	std::string s =
76	137	"#version 450\n"
77	138	"layout(local_size_x = 8, local_size_y = 8) in;\n"
…	…	@@ -82,6 +143,8 @@ static std::string build_render_shader(const std::string& inc,
82	143	"\n";
83	144	s += inc;
84	145	s += "\n";
	146	// Must precede the wrapper: heightmap_gn calls hblk_height_m.
	147	if (hb) s += hblk_sampler_glsl(*hb);
85	148	s += erosion_wrapper_glsl(wp);
86	149
87	150	char tail[2048];
…	…	@@ -130,6 +193,7 @@ struct Renderer {
130	193	VkPipeline pipe = VK_NULL_HANDLE;
131	194	VkDescriptorPool dpool = VK_NULL_HANDLE;
132	195	VkDescriptorSet dset = VK_NULL_HANDLE;
	196	std::vector<float> hblk; // the heightfield, empty when not in use
133	197
134	198	void build(const std::string& src) {
135	199	init_vulkan(c);
…	…	@@ -139,13 +203,15 @@ struct Renderer {
139	203	smi.codeSize = spv.size() * 4; smi.pCode = spv.data();
140	204	VK(vkCreateShaderModule(c.dev, &smi, nullptr, &mod));
141	205
142		VkDescriptorSetLayoutBinding bind{};
143		bind.binding = 0;
144		bind.descriptorType = VK_DESCRIPTOR_TYPE_STORAGE_BUFFER;
145		bind.descriptorCount = 1;
146		bind.stageFlags = VK_SHADER_STAGE_COMPUTE_BIT;
	206	VkDescriptorSetLayoutBinding bind[2]{};
	207	for (int k = 0; k < 2; k++) {
	208	bind[k].binding = (uint32_t)k;
	209	bind[k].descriptorType = VK_DESCRIPTOR_TYPE_STORAGE_BUFFER;
	210	bind[k].descriptorCount = 1;
	211	bind[k].stageFlags = VK_SHADER_STAGE_COMPUTE_BIT;
	212	}
147	213	VkDescriptorSetLayoutCreateInfo dli{VK_STRUCTURE_TYPE_DESCRIPTOR_SET_LAYOUT_CREATE_INFO};
148		dli.bindingCount = 1; dli.pBindings = &bind;
	214	dli.bindingCount = 2; dli.pBindings = bind;
149	215	VK(vkCreateDescriptorSetLayout(c.dev, &dli, nullptr, &dsl));
150	216
151	217	VkPushConstantRange pcr{VK_SHADER_STAGE_COMPUTE_BIT, 0, sizeof(PushC)};
…	…	@@ -161,7 +227,7 @@ struct Renderer {
161	227	cpi.layout = playout;
162	228	VK(vkCreateComputePipelines(c.dev, VK_NULL_HANDLE, 1, &cpi, nullptr, &pipe));
163	229
164		VkDescriptorPoolSize ps{VK_DESCRIPTOR_TYPE_STORAGE_BUFFER, 1};
	230	VkDescriptorPoolSize ps{VK_DESCRIPTOR_TYPE_STORAGE_BUFFER, 2};
165	231	VkDescriptorPoolCreateInfo dpi{VK_STRUCTURE_TYPE_DESCRIPTOR_POOL_CREATE_INFO};
166	232	dpi.maxSets = 1; dpi.poolSizeCount = 1; dpi.pPoolSizes = &ps;
167	233	VK(vkCreateDescriptorPool(c.dev, &dpi, nullptr, &dpool));
…	…	@@ -179,12 +245,25 @@ struct Renderer {
179	245	out.assign((size_t)w * h, NO_SURFACE);
180	246	Buf buf = make_host_buf(c, (size_t)w * band * sizeof(float));
181	247
182		VkDescriptorBufferInfo bi{buf.buf, 0, VK_WHOLE_SIZE};
183		VkWriteDescriptorSet wr{VK_STRUCTURE_TYPE_WRITE_DESCRIPTOR_SET};
184		wr.dstSet = dset; wr.dstBinding = 0; wr.descriptorCount = 1;
185		wr.descriptorType = VK_DESCRIPTOR_TYPE_STORAGE_BUFFER;
186		wr.pBufferInfo = &bi;
187		vkUpdateDescriptorSets(c.dev, 1, &wr, 0, nullptr);
	248	// Binding 1 always exists in the layout; with no heightfield it is a
	249	// one-float placeholder the shader never declares. Leaving it UNWRITTEN
	250	// would be undefined behaviour rather than merely unused.
	251	size_t hb_bytes = hblk.empty() ? sizeof(float) : hblk.size()*sizeof(float);
	252	Buf hbuf = make_host_buf(c, hb_bytes);
	253	if (!hblk.empty()) memcpy(hbuf.mapped, hblk.data(), hb_bytes);
	254	else ((float*)hbuf.mapped)[0] = 0.0f;
	255
	256	VkDescriptorBufferInfo bi[2] = {{buf.buf, 0, VK_WHOLE_SIZE},
	257	{hbuf.buf, 0, VK_WHOLE_SIZE}};
	258	VkWriteDescriptorSet wr[2]{};
	259	for (int k = 0; k < 2; k++) {
	260	wr[k].sType = VK_STRUCTURE_TYPE_WRITE_DESCRIPTOR_SET;
	261	wr[k].dstSet = dset; wr[k].dstBinding = (uint32_t)k;
	262	wr[k].descriptorCount = 1;
	263	wr[k].descriptorType = VK_DESCRIPTOR_TYPE_STORAGE_BUFFER;
	264	wr[k].pBufferInfo = &bi[k];
	265	}
	266	vkUpdateDescriptorSets(c.dev, 2, wr, 0, nullptr);
188	267
189	268	for (uint32_t row0 = 0; row0 < h; row0 += band) {
190	269	uint32_t rows = (row0 + band <= h) ? band : (h - row0);
…	…	@@ -273,6 +352,8 @@ int main(int argc, char** argv) {
273	352	int n = 512, steps = 512, bisect = 60, band = 32;
274	353	float mpp = 30.0f, ox = 0.0f, oz = 0.0f;
275	354	bool verify = false;
	355	const char* hblk_path = nullptr;
	356	bool have_geom = false;
276	357
277	358	// NOT ErosionWrapperParams{} -- the struct's own defaults are a convenience
278	359	// and are not authoritative. The map's values come from the oracle's single
…	…	@@ -285,10 +366,11 @@ int main(int argc, char** argv) {
285	366	const char* a = argv[i];
286	367	if (!strcmp(a, "--inc") && i+1 < argc) inc_path = argv[++i];
287	368	else if (!strcmp(a, "--out") && i+1 < argc) out_path = argv[++i];
288		else if (!strcmp(a, "--n") && i+1 < argc) n = atoi(argv[++i]);
289		else if (!strcmp(a, "--mpp") && i+1 < argc) mpp = (float)atof(argv[++i]);
290		else if (!strcmp(a, "--ox") && i+1 < argc) ox = (float)atof(argv[++i]);
291		else if (!strcmp(a, "--oz") && i+1 < argc) oz = (float)atof(argv[++i]);
	369	else if (!strcmp(a, "--n") && i+1 < argc) { n = atoi(argv[++i]); have_geom = true; }
	370	else if (!strcmp(a, "--mpp") && i+1 < argc) { mpp = (float)atof(argv[++i]); have_geom = true; }
	371	else if (!strcmp(a, "--ox") && i+1 < argc) { ox = (float)atof(argv[++i]); have_geom = true; }
	372	else if (!strcmp(a, "--oz") && i+1 < argc) { oz = (float)atof(argv[++i]); have_geom = true; }
	373	else if (!strcmp(a, "--heightmap") && i+1 < argc) hblk_path = argv[++i];
292	374	else if (!strcmp(a, "--steps") && i+1 < argc) steps = atoi(argv[++i]);
293	375	else if (!strcmp(a, "--bisect") && i+1 < argc) bisect = atoi(argv[++i]);
294	376	else if (!strcmp(a, "--band") && i+1 < argc) band = atoi(argv[++i]);
…	…	@@ -298,9 +380,29 @@ int main(int argc, char** argv) {
298	380	}
299	381
300	382	std::string inc = read_file(inc_path);
301		std::string src = build_render_shader(inc, wp, steps, bisect);
	383
	384	HeightBlock hb;
	385	bool use_hblk = (hblk_path != nullptr);
	386	if (use_hblk) {
	387	hb = read_hblk(hblk_path);
	388	wp.heightmap_source = true;
	389	wp.heightmap_step = hb.mpp;
	390	if (!have_geom) {
	391	n = hb.n - 2; // render the interior; the ring is apron
	392	mpp = hb.mpp;
	393	ox = hb.ox + hb.mpp;
	394	oz = hb.oz + hb.mpp;
	395	}
	396	printf("heightfield: %s\n", hblk_path);
	397	printf(" %d x %d px at %g m/px, world x %.0f z %.0f\n",
	398	hb.n, hb.n, hb.mpp, hb.ox, hb.oz);
	399	}
	400
	401	std::string src = build_render_shader(inc, wp, steps, bisect,
	402	use_hblk ? &hb : nullptr);
302	403
303	404	Renderer r;
	405	if (use_hblk) r.hblk = hb.h;
304	406	r.build(src);
305	407
306	408	if (verify) {

gate_dc_vs_lookdev.py +82 -0

…	…	@@ -0,0 +1,82 @@
	1	#!/usr/bin/env python3
	2	"""
	3	Compare dc-terrain's surface (from Godot) against the lookdev harness.
	4
	5	The Godot half is scenes/test/gate_dc_vs_lookdev.gd, which walks
	6	DCTerrain.sample_density and writes /tmp/dc_surface.txt. This renders the same
	7	columns with erosion_render and compares. Two independent implementations --
	8	GDScript-driven C++ inside the engine, and a standalone Vulkan compute shader --
	9	so agreement is evidence rather than a tautology.
	10
	11	THE CONTROL renders the same columns with erosion turned OFF and requires the
	12	comparison to REJECT it. Without that, a comparison that passed everything
	13	would read as success.
	14
	15	python3 gate_dc_vs_lookdev.py [--bin /tmp/erosion_render4]
	16	"""
	17	import argparse, os, subprocess, sys
	18	import numpy as np
	19
	20	DC = os.path.abspath(os.path.join(os.path.dirname(__file__), "..", "dc-terrain"))
	21	HBLK = os.path.abspath(os.path.join(os.path.dirname(__file__), "..",
	22	"hyper-tps", "assets", "terrain",
	23	"diffusion_gully_30km.hblk"))
	24
	25
	26	def render(binpath, ox, oz, n, mpp, out, extra=()):
	27	cmd = [binpath, "--heightmap", HBLK, "--ox", str(ox), "--oz", str(oz),
	28	"--n", str(n), "--mpp", str(mpp), "--out", out] + list(extra)
	29	r = subprocess.run(cmd, cwd=DC, capture_output=True, text=True)
	30	if r.returncode != 0:
	31	print(r.stdout, r.stderr); raise SystemExit("renderer failed")
	32	return np.fromfile(out, dtype="<f4", count=n * n).astype(np.float64).reshape(n, n)
	33
	34
	35	def main():
	36	ap = argparse.ArgumentParser()
	37	ap.add_argument("--bin", default="/tmp/erosion_render4")
	38	ap.add_argument("--dump", default="/tmp/dc_surface.txt")
	39	ap.add_argument("--tol", type=float, default=2.0, help="metres")
	40	a = ap.parse_args()
	41
	42	rows = []
	43	for line in open(a.dump):
	44	if line.startswith("#"):
	45	continue
	46	x, z, y = line.split()
	47	rows.append((float(x), float(z), float(y)))
	48	xs = sorted(set(r[0] for r in rows)); zs = sorted(set(r[1] for r in rows))
	49	n = len(xs)
	50	mpp = xs[1] - xs[0]
	51	godot = np.zeros((len(zs), n))
	52	for x, z, y in rows:
	53	godot[zs.index(z)][xs.index(x)] = y
	54	print(f"Godot surface: {n}x{len(zs)} columns from ({xs[0]:.0f}, {zs[0]:.0f}) "
	55	f"at {mpp:.0f} m spacing")
	56
	57	look = render(a.bin, xs[0], zs[0], n, mpp, "/tmp/lk_surface.f32")
	58	d = np.abs(godot - look)
	59	ok = d.max() <= a.tol
	60	print(f"\n MATCH max \|godot - lookdev\| = {d.max():8.3f} m "
	61	f"rms {np.sqrt((d**2).mean()):7.3f} m {'PASS' if ok else 'FAIL'}")
	62
	63	zero = render(a.bin, xs[0], zs[0], n, mpp, "/tmp/lk_zero.f32",
	64	["--set", "en_strength", "0"])
	65	dz = np.abs(godot - zero)
	66	ctl = dz.max() > a.tol
	67	print(f" CONTROL erosion off: max diff {dz.max():8.1f} m "
	68	f"{'REJECTED (gate discriminates)' if ctl else 'ACCEPTED — GATE IS BLIND'}")
	69
	70	print()
	71	if not ctl:
	72	print("VERDICT: the comparison accepts a different field, so MATCH means nothing.")
	73	return 2
	74	if not ok:
	75	print("VERDICT: dc-terrain and the lookdev DISAGREE about the terrain.")
	76	return 1
	77	print("VERDICT: dc-terrain's field is the lookdev's field.")
	78	return 0
	79
	80
	81	if __name__ == "__main__":
	82	sys.exit(main())

gate_heightmap.py +127 -0

…	…	@@ -0,0 +1,127 @@
	1	#!/usr/bin/env python3
	2	"""
	3	GATE: the erosion filter running on a diffusion heightfield, in GLSL.
	4
	5	Three rows, and the second and third exist so the first cannot pass by accident.
	6
	7	IDENTITY strength 0, no height offset, no detail layer. The surface the
	8	renderer finds MUST equal the input heightfield. This is the check
	9	erode_heightmap.cpp's --check does for the CPU shim, moved to the
	10	GPU path: a unit error in the metres -> filter-space conversion
	11	cannot survive it, and that conversion is exactly where #45 lost a
	12	session to a 60x factor.
	13
	14	CONTROL the same run at the REAL strength must NOT match. An identity gate
	15	that passes whatever you feed it is measuring nothing; this shows
	16	it discriminates.
	17
	18	SEAM the eroded field's step across an internal tile join must be no
	19	worse than the steps just inside the neighbouring tiles. This is
	20	the one that would catch the per-tile normalisation bug -- and it
	21	needs a block of at least 2x2 tiles to be able to fail at all,
	22	which is why stitch_block defaults to 4.
	23
	24	python gate_heightmap.py /tmp/block.hblk --bin /tmp/erosion_render2
	25	"""
	26	import argparse, os, subprocess, sys
	27	import numpy as np
	28	from stitch_block import read_hblk
	29
	30	DC = os.path.abspath(os.path.join(os.path.dirname(__file__), "..", "dc-terrain"))
	31	TILE_PX = 256
	32
	33
	34	def run(binpath, hblk, out, extra):
	35	cmd = [binpath, "--heightmap", hblk, "--out", out] + extra
	36	r = subprocess.run(cmd, cwd=DC, capture_output=True, text=True)
	37	if r.returncode != 0:
	38	print(r.stdout, r.stderr); raise SystemExit(f"renderer failed: {cmd}")
	39	return r.stdout
	40
	41
	42	def load(out, n):
	43	a = np.fromfile(out, dtype="<f4", count=n * n).reshape(n, n)
	44	return a.astype(np.float64)
	45
	46
	47	def main():
	48	ap = argparse.ArgumentParser()
	49	ap.add_argument("hblk"); ap.add_argument("--bin", default="/tmp/erosion_render2")
	50	ap.add_argument("--strength", default="0.15")
	51	ap.add_argument("--tol", type=float, default=1.0, help="identity tolerance, metres")
	52	a = ap.parse_args()
	53
	54	arr, meta = read_hblk(a.hblk)
	55	n = meta["side"] - 2
	56	ref = arr[1:-1, 1:-1].astype(np.float64) # interior = what is rendered
	57	print(f"block {meta['side']}x{meta['side']} px, rendering interior {n}x{n}"
	58	f" ({n*meta['mpp']/1000:.2f} km)\n")
	59
	60	# ---- IDENTITY ---------------------------------------------------------
	61	run(a.bin, a.hblk, "/tmp/gate_id.f32",
	62	["--set", "en_strength", "0", "--set", "en_height_offset_a", "0",
	63	"--set", "en_height_offset_b", "0", "--set", "detail_scale", "0"])
	64	idf = load("/tmp/gate_id.f32", n)
	65	d = np.abs(idf - ref)
	66	id_ok = d.max() <= a.tol
	67	print(f" IDENTITY maxdiff {d.max():8.3f} m rmse {np.sqrt((d**2).mean()):8.4f} m"
	68	f" {'PASS' if id_ok else 'FAIL'}")
	69
	70	# ---- CONTROL ----------------------------------------------------------
	71	run(a.bin, a.hblk, "/tmp/gate_er.f32",
	72	["--set", "en_strength", a.strength, "--set", "detail_scale", "0"])
	73	er = load("/tmp/gate_er.f32", n)
	74	de = np.abs(er - ref)
	75	ctl_ok = de.max() > a.tol
	76	print(f" CONTROL strength {a.strength}: maxdiff {de.max():8.1f} m "
	77	f"rmse {np.sqrt((de**2).mean()):7.1f} m "
	78	f"{'DIFFERS (gate discriminates)' if ctl_ok else 'IDENTICAL -- GATE IS BLIND'}")
	79
	80	# ---- SEAM -------------------------------------------------------------
	81	# MEAN SIGNED step across the join, not mean absolute.
	82	#
	83	# The first version of this check used mean \|step\| and was BLIND: a 60 m
	84	# per-tile offset injected on purpose moved it from 1.007 to 1.088 against
	85	# a 1.5 threshold, because the terrain's own pixel-to-pixel roughness is
	86	# far larger than the seam and drowns it. A normalisation error is a
	87	# CONSISTENT SIGNED jump; natural roughness averages to about zero. So the
	88	# signed mean isolates exactly the thing being looked for, and is compared
	89	# against the spread of the same statistic at nearby interior columns.
	90	def col_signed_step(f, b):
	91	return float(np.mean(f[:, b + 1] - f[:, b]))
	92
	93	bounds = [k * TILE_PX - 1 for k in range(1, n // TILE_PX)]
	94	if not bounds:
	95	print(" SEAM SKIPPED -- block is 1 tile wide, this check cannot fail")
	96	seam_ok = False
	97	else:
	98	worst_z, rows = 0.0, []
	99	for b in bounds:
	100	at = col_signed_step(er, b)
	101	near = [col_signed_step(er, c) for c in range(b - 24, b + 25)
	102	if abs(c - b) > 2 and 0 <= c < n - 1]
	103	sd = float(np.std(near)) or 1e-9
	104	z = abs(at) / sd
	105	worst_z = max(worst_z, z)
	106	rows.append((b, at, sd, z))
	107	for b, at, sd, z in rows:
	108	print(f" SEAM col {b:>5}: signed step {at:+8.3f} m "
	109	f"neighbours sd {sd:6.3f} m z {z:6.2f}")
	110	seam_ok = bool(worst_z < 5.0)
	111	print(f" worst z {worst_z:.2f} "
	112	f"{'PASS (<5)' if seam_ok else 'FAIL -- the join is a systematic step'}")
	113
	114	print()
	115	if not ctl_ok:
	116	print("VERDICT: the control did not move, so IDENTITY proves nothing."); return 2
	117	if not id_ok:
	118	print("VERDICT: units are WRONG. Do not trust any eroded field."); return 1
	119	if not seam_ok:
	120	print("VERDICT: identity holds but the tile joins are visible."); return 1
	121	print("VERDICT: units exact, erosion active, joins invisible. The shared "
	122	"wrapper's heightmap path is correct on the GPU.")
	123	return 0
	124
	125
	126	if __name__ == "__main__":
	127	sys.exit(main())

gate_heightmap_lod.py +106 -0

…	…	@@ -0,0 +1,106 @@
	1	#!/usr/bin/env python3
	2	"""
	3	GATE: is the baked heightfield BAND-LIMITED per LOD?
	4
	5	The property the octree clipmap needs is not "the coarse surface is close to the
	6	fine one" -- it is "the coarse surface is a SMOOTHED version of the fine one".
	7	Every other layer in this generator gets that from nyquist_octaves(..., voxel_scale,
	8	..), which drops octaves as the voxel grows. A sampled heightfield has no such
	9	property unless somebody gives it one.
	10
	11	Without it, LOD k samples a 30 m field at 2^k * 30 m spacing. At LOD 7 that is
	12	points 3840 m apart in a field whose sd is 1167 m, so neighbouring voxels differ by
	13	about 1.6 km and the coarse mesh is noise rather than terrain. geometric_error then
	14	never falls (dc_octree_clipmap.cpp:328-351 uses it for both the content-aware split
	15	boost and the screen-space test), so the clipmap splits forever and never converges.
	16
	17	THE CHECK: render the surface at voxel V on a V-spaced grid, and compare against
	18	the FINE surface box-averaged over V. A band-limited field matches; a point-sampled
	19	one does not.
	20
	21	ALIGNMENT, which this gate got wrong first time round and reported as a code bug.
	22	Level-k texel j is the average of level-0 texels [j2^k, (j+1)2^k), so its CENTRE
	23	sits at ox + (j2^k + (2^k-1)/2)mpp -- half a coarse texel to the right of
	24	ox + jV. Sampling the renderer at ox + jV and comparing against box j compares
	25	two different places, and reads as 57 m of error at LOD 4 where there is 13. The
	26	coarse render is therefore offset onto the box centres.
	27
	28	TOLERANCE IS RELATIVE TO THE VOXEL. A difference far below the voxel size cannot
	29	be represented by that LOD's mesh and is not a defect; what matters is whether the
	30	coarse surface is the fine one smoothed, to within what the LOD can express.
	31
	32	THE CONTROL is V = 30, the field's own pitch, where the two are the same thing by
	33	construction. If the control ever fails the comparison itself is wrong.
	34
	35	python gate_heightmap_lod.py assets/terrain/diffusion_gully_30km.hblk
	36	"""
	37	import argparse, os, subprocess, sys
	38	import numpy as np
	39
	40	DC = os.path.abspath(os.path.join(os.path.dirname(__file__), "..", "dc-terrain"))
	41	MPP = 30.0
	42
	43
	44	def render(binpath, hblk, ox, oz, n, mpp, voxel, out):
	45	cmd = [binpath, "--heightmap", hblk, "--ox", str(ox), "--oz", str(oz),
	46	"--n", str(n), "--mpp", str(mpp), "--out", out,
	47	"--set", "voxel_size", str(voxel)]
	48	r = subprocess.run(cmd, cwd=DC, capture_output=True, text=True)
	49	if r.returncode != 0:
	50	print(r.stdout, r.stderr); raise SystemExit("renderer failed")
	51	return np.fromfile(out, dtype="<f4", count=n * n).astype(np.float64).reshape(n, n)
	52
	53
	54	def main():
	55	ap = argparse.ArgumentParser()
	56	ap.add_argument("hblk"); ap.add_argument("--bin", default="/tmp/erosion_render2")
	57	ap.add_argument("--ox", type=float, default=45000.0)
	58	ap.add_argument("--oz", type=float, default=30000.0)
	59	ap.add_argument("--fine-n", type=int, default=256)
	60	ap.add_argument("--tol-frac", type=float, default=0.10,
	61	help="rms tolerance as a fraction of the voxel size")
	62	a = ap.parse_args()
	63
	64	hb = os.path.abspath(a.hblk)
	65	span = a.fine_n * MPP
	66	print(f"window {span/1000:.2f} km at ({a.ox:.0f}, {a.oz:.0f}), "
	67	f"fine {a.fine_n}x{a.fine_n} @ {MPP:g} m\n")
	68
	69	fine = render(a.bin, hb, a.ox, a.oz, a.fine_n, MPP, 1.0, "/tmp/lod_fine.f32")
	70
	71	print(f"{'voxel m':>8} {'LOD':>4} {'rms vs box-avg':>15} {'tol':>8} verdict")
	72	fails = 0
	73	for V in [30.0, 120.0, 480.0, 1920.0]:
	74	s = int(round(V / MPP))
	75	n = a.fine_n // s
	76	if n < 4:
	77	continue
	78	# onto the box centres — see ALIGNMENT above
	79	shift = (s - 1) * MPP * 0.5
	80	coarse = render(a.bin, hb, a.ox + shift, a.oz + shift, n, V, V,
	81	f"/tmp/lod_{int(V)}.f32")
	82	m = n * s
	83	box = fine[:m, :m].reshape(n, s, n, s).mean(axis=(1, 3))
	84	rms = float(np.sqrt(((coarse - box) ** 2).mean()))
	85	tol = max(a.tol_frac * V, 5.0)
	86	ok = rms <= tol
	87	tag = "PASS" if ok else "FAIL"
	88	if V == MPP:
	89	tag = ("CONTROL ok" if ok else "CONTROL BROKEN — the comparison is wrong")
	90	elif not ok:
	91	fails += 1
	92	print(f"{V:>8.0f} {int(np.log2(V/MPP)):>4} {rms:>15.1f} {tol:>8.1f} {tag}")
	93
	94	print()
	95	if fails == 0:
	96	print("VERDICT: the field is band-limited — every LOD is a smoothed version")
	97	print(" of the one below it, which is what the clipmap needs to settle.")
	98	return 0
	99	print("VERDICT: the field is NOT band-limited. Coarse LODs point-sample a 30 m")
	100	print(" field, so they are noise rather than smoothed terrain, and the")
	101	print(" clipmap has nothing to converge to.")
	102	return 1
	103
	104
	105	if __name__ == "__main__":
	106	sys.exit(main())

gate_stitch.py +83 -0

…	…	@@ -0,0 +1,83 @@
	1	#!/usr/bin/env python3
	2	"""
	3	GATE: does the stitched .hblk describe the same world region the service does?
	4
	5	This is the check that would have caught #50 section 2.4 -- the transposed tile
	6	layout that cost a session, where the instrument and the page shared the same
	7	axis bug and therefore agreed with each other.
	8
	9	It re-fetches the block's exact world footprint from the live service in ONE
	10	request and compares to the stitched array. Tolerance is 1 m: tiles are stored
	11	int16 and tiles.py documents that fetching a region as one request or as
	12	quadrants differs by at most that.
	13
	14	AND IT PROVES IT CAN FAIL. A gate that cannot fail is not evidence, so the
	15	CONTROL row feeds the same comparison a transposed copy of the same data and
	16	requires it to be REJECTED. If the control passes, the gate is broken and the
	17	main row means nothing.
	18
	19	python gate_stitch.py /tmp/block.hblk [--url http://127.0.0.1:8817]
	20	"""
	21	import argparse, sys, time, urllib.request
	22	import numpy as np
	23	from stitch_block import read_hblk, MPP
	24
	25
	26	def fetch(url, i1, j1, n):
	27	u = f"{url}/terrain?i1={i1}&j1={j1}&i2={i1+n}&j2={j1+n}&scale=1"
	28	t0 = time.time()
	29	with urllib.request.urlopen(u, timeout=7200) as r:
	30	h = int(r.headers["X-Height"]); w = int(r.headers["X-Width"]); b = r.read()
	31	return (np.frombuffer(b[:hw2], dtype="<i2").reshape(h, w).astype(np.int64),
	32	time.time() - t0)
	33
	34
	35	def compare(a, b, tol=1):
	36	d = np.abs(a - b)
	37	return dict(maxdiff=int(d.max()), rmse=float(np.sqrt((d**2).mean())),
	38	frac_over=float((d > tol).mean()), ok=bool(d.max() <= tol))
	39
	40
	41	def main():
	42	ap = argparse.ArgumentParser()
	43	ap.add_argument("hblk"); ap.add_argument("--url", default="http://127.0.0.1:8817")
	44	ap.add_argument("--tol", type=int, default=1)
	45	a = ap.parse_args()
	46
	47	arr, meta = read_hblk(a.hblk)
	48	side = meta["side"]
	49	# The header stores world metres; the service is indexed in pixels, and
	50	# row = z, column = x (see stitch_block.py's axis note).
	51	i1 = int(round(meta["origin_z"] / meta["mpp"])) # first axis
	52	j1 = int(round(meta["origin_x"] / meta["mpp"])) # second axis
	53	print(f"block {side}x{side} px at world x {meta['origin_x']:,.0f} "
	54	f"z {meta['origin_z']:,.0f} -> service i1={i1} j1={j1} n={side}")
	55
	56	live, dt = fetch(a.url, i1, j1, side)
	57	print(f"fetched in {dt:.1f}s")
	58	disk = arr.astype(np.int64)
	59
	60	r = compare(disk, live, a.tol)
	61	print(f"\n STITCH maxdiff {r['maxdiff']:>5} m rmse {r['rmse']:>7.3f} m "
	62	f"over-tol {100*r['frac_over']:.4f}% {'PASS' if r['ok'] else 'FAIL'}")
	63
	64	# CONTROL: the same comparison must REJECT a transposed block.
	65	c = compare(disk.T, live, a.tol)
	66	ctrl_ok = not c["ok"]
	67	print(f" CONTROL transposed copy -> maxdiff {c['maxdiff']:>5} m "
	68	f"{'REJECTED (gate can fail)' if ctrl_ok else 'ACCEPTED -- GATE IS BROKEN'}")
	69
	70	if not ctrl_ok:
	71	print("\nVERDICT: the gate accepted a transposed block, so it is not "
	72	"measuring layout at all. The PASS above means nothing.")
	73	return 2
	74	if not r["ok"]:
	75	print("\nVERDICT: stitched block does NOT match the service.")
	76	return 1
	77	print("\nVERDICT: stitched block matches the service, and the gate rejects "
	78	"a layout error. The .hblk is the world the lookdev showed.")
	79	return 0
	80
	81
	82	if __name__ == "__main__":
	83	sys.exit(main())

stitch_block.py +150 -0

…	…	@@ -0,0 +1,150 @@
	1	#!/usr/bin/env python3
	2	"""
	3	STITCH a solid block of lookdev cache tiles into one .hblk heightfield.
	4
	5	WHY THIS EXISTS
	6
	7	The lookdev cache is 756 int16 tiles of 256 px + a 1 px apron, written by
	8	tiles.py while someone flew around /explore. They are a verified-single-world
	9	heightfield (8 samples re-fetched from the live fp32 seed-1234 service matched
	10	to <= 1 m, which is the int16 quantisation tiles.py already documents), and
	11	they are the cheapest possible input for proving the engine path: no service,
	12	no PyTorch, no network.
	13
	14	AXES, WHICH ARE THE EASY THING TO GET WRONG
	15
	16	#50 section 2.4 cost a session to this. tiles.py fetches a tile as
	17	fetch(i1 = tx256, j1 = ty256) and the API returns [i, j] -- so:
	18
	19	tx indexes the FIRST array axis -> WORLD Z
	20	ty indexes the SECOND array axis -> WORLD X
	21
	22	erode_heightmap.cpp uses the same convention on its raw input: it indexes
	23	h[jn + i] with px = impp and pz = j*mpp, i.e. column = x, row = z. The .hblk
	24	below is row-major with row = z and column = x, so both agree by construction
	25	rather than by luck.
	26
	27	THE APRON IS KEPT, AND IS THE POINT
	28
	29	Each tile carries 1 px of its neighbour on every side. Stitching interiors
	30	only would leave the block's outer ring without the neighbours the central
	31	difference needs, and the slope error there would look exactly like a seam.
	32	The block therefore spans N*256 + 2 px and the consumer treats the outer ring
	33	as apron: valid heights, not valid slopes.
	34
	35	Usage:
	36	python stitch_block.py --n 4 --out /tmp/block.hblk # auto-pick
	37	python stitch_block.py --tx -3 --ty -1 --n 4 --out ... # explicit
	38	python stitch_block.py --list # what's solid
	39	"""
	40	import argparse, os, re, struct, sys
	41	import numpy as np
	42
	43	CACHE = os.environ.get("LOOKDEV_TILE_CACHE",
	44	os.path.expanduser("~/.cache/terrain-lookdev/tiles"))
	45	TILE_PX = 256
	46	APRON = 1
	47	MPP = 30.0
	48	MAGIC = b"HBLK"
	49	VERSION = 1
	50	HEADER = 32 # magic4 + version4 + n4 + mpp4 + origin_x4 + origin_z4 + 8 reserved
	51
	52
	53	def cached_tiles(cache=CACHE, tile_px=TILE_PX):
	54	pat = re.compile(rf"^t2_{tile_px}_(-?\d+)_(-?\d+)\.npy$")
	55	out = set()
	56	for f in os.listdir(cache):
	57	m = pat.match(f)
	58	if m:
	59	out.add((int(m.group(1)), int(m.group(2))))
	60	return out
	61
	62
	63	def solid_origins(tiles, n):
	64	"""Tile origins (tx, ty) where the whole n x n block is present."""
	65	return sorted(o for o in tiles
	66	if all((o[0] + i, o[1] + j) in tiles
	67	for i in range(n) for j in range(n)))
	68
	69
	70	def stitch(tx0, ty0, n, cache=CACHE, tile_px=TILE_PX):
	71	"""Return (array, origin_x_m, origin_z_m). Array is (n*tile_px+2)^2 int16,
	72	row = z, column = x, with a 1 px apron ring."""
	73	side = n * tile_px + 2 * APRON
	74	out = np.zeros((side, side), dtype=np.int16)
	75	for i in range(n): # tx -> z
	76	for j in range(n): # ty -> x
	77	p = os.path.join(cache, f"t2_{tile_px}_{tx0+i}_{ty0+j}.npy")
	78	t = np.load(p) # (tile_px+2, tile_px+2)
	79	if t.shape != (tile_px + 2, tile_px + 2):
	80	raise SystemExit(f"{p}: expected {(tile_px+2,)*2}, got {t.shape}")
	81	# Write the tile INCLUDING its apron. Interiors never overlap;
	82	# aprons overwrite each other with identical values (same world
	83	# pixel, same world), so order does not matter -- and the block's
	84	# outer ring ends up filled by the edge tiles' aprons, which is
	85	# the only way to get it at all.
	86	r0, c0 = i * tile_px, j * tile_px
	87	out[r0:r0 + tile_px + 2, c0:c0 + tile_px + 2] = t
	88	origin_z = (tx0 * tile_px - APRON) * MPP
	89	origin_x = (ty0 * tile_px - APRON) * MPP
	90	return out, origin_x, origin_z
	91
	92
	93	def write_hblk(path, arr, origin_x, origin_z, mpp=MPP):
	94	side = arr.shape[0]
	95	hdr = struct.pack("<4sIIfffII", MAGIC, VERSION, side, mpp,
	96	origin_x, origin_z, 0, 0)
	97	assert len(hdr) == HEADER, len(hdr)
	98	with open(path, "wb") as f:
	99	f.write(hdr)
	100	f.write(np.ascontiguousarray(arr, dtype="<i2").tobytes())
	101
	102
	103	def read_hblk(path):
	104	with open(path, "rb") as f:
	105	hdr = f.read(HEADER)
	106	magic, ver, side, mpp, ox, oz, _, _ = struct.unpack("<4sIIfffII", hdr)
	107	if magic != MAGIC:
	108	raise SystemExit(f"{path}: not an hblk (magic {magic!r})")
	109	a = np.frombuffer(f.read(side * side * 2), dtype="<i2").reshape(side, side)
	110	return a, dict(version=ver, side=side, mpp=mpp, origin_x=ox, origin_z=oz)
	111
	112
	113	def main():
	114	ap = argparse.ArgumentParser()
	115	ap.add_argument("--n", type=int, default=4, help="block side in TILES")
	116	ap.add_argument("--tx", type=int, default=None)
	117	ap.add_argument("--ty", type=int, default=None)
	118	ap.add_argument("--out", default="/tmp/block.hblk")
	119	ap.add_argument("--list", action="store_true")
	120	a = ap.parse_args()
	121
	122	tiles = cached_tiles()
	123	if a.list:
	124	for n in (1, 2, 4, 8, 16):
	125	s = solid_origins(tiles, n)
	126	print(f"solid {n:>2}x{n:<2} blocks: {len(s):>4}"
	127	f" span {nTILE_PXMPP/1000:>6.1f} km"
	128	+ (f" e.g. --tx {s[len(s)//2][0]} --ty {s[len(s)//2][1]}" if s else ""))
	129	return
	130
	131	if a.tx is None or a.ty is None:
	132	s = solid_origins(tiles, a.n)
	133	if not s:
	134	raise SystemExit(f"no solid {a.n}x{a.n} block in {CACHE}")
	135	a.tx, a.ty = s[len(s) // 2]
	136	print(f"auto-picked block origin tile ({a.tx}, {a.ty})")
	137
	138	arr, ox, oz = stitch(a.tx, a.ty, a.n)
	139	write_hblk(a.out, arr, ox, oz)
	140	side = arr.shape[0]
	141	print(f"wrote {a.out}")
	142	print(f" {side} x {side} px ({a.n}x{a.n} tiles + 1 px apron)"
	143	f" = {(side-2)*MPP/1000:.2f} km span")
	144	print(f" world origin x {ox:,.0f} m z {oz:,.0f} m at {MPP:g} m/px")
	145	print(f" elevation {arr.min()} .. {arr.max()} m")
	146	print(f" {os.path.getsize(a.out)/1024:.0f} KB")
	147
	148
	149	if __name__ == "__main__":
	150	main()

terrain-lookdev commits / 05b8b3a

Changed files (8) +756 -19