| … | … | @@ -0,0 +1,550 @@ |
| 1 | """ |
| 2 | Proves the bake is an IDENTITY map in world space, end to end, against a fractal |
| 3 | whose position is known analytically. |
| 4 | |
| 5 | `test_volume_axes.py` proves `write_vdb`'s conventions in ISOLATION: it hands the |
| 6 | writer a synthetic `[z][y][x]` array and checks the marker comes out at the right |
| 7 | world position. That is necessary and not sufficient, because it ASSUMES the array |
| 8 | `accumulate` returns is indexed `[z][y][x]` in the first place. That assumption was |
| 9 | measured once by hand and never asserted, and everything downstream — the transpose |
| 10 | in `write_vdb`, the axis reorder in `_mass_bounds` — is built on it. If the readback |
| 11 | order were different, both would be consistently wrong and both tests would pass. |
| 12 | |
| 13 | So this starts at the GPU and finishes in world coordinates: |
| 14 | |
| 15 | document -> chaos game on the GPU -> voxel grid -> readback -> .vdb -> world |
| 16 | |
| 17 | and compares the result to where the fractal provably is. |
| 18 | |
| 19 | ## The ground truth |
| 20 | |
| 21 | The document is four Affine iterators, `p := 0.5*p + 0.5*F_i`, so the attractor is |
| 22 | the Sierpinski tetrahedron with vertices `F_i` — and TWO facts about it are exact, |
| 23 | not estimated: |
| 24 | |
| 25 | * its **convex hull** is the tetrahedron, so the per-axis extent is just the range |
| 26 | of the `F_i` components; |
| 27 | * its **centroid** is the mean of the `F_i`. The invariant measure satisfies |
| 28 | `mu = E[0.5*p + 0.5*F_i] = 0.5*mu + 0.5*mean(F)`, hence `mu = mean(F)`. |
| 29 | |
| 30 | The vertices are chosen so every axis is distinguishable from every other — three |
| 31 | different spans, three different centroid components, none symmetric about zero. |
| 32 | A transposed, mirrored or displaced bake therefore cannot land on the right answer |
| 33 | by coincidence, which is exactly the failure mode a symmetric test document hides. |
| 34 | `test_the_document_can_detect_a_swap` guards that property. |
| 35 | |
| 36 | ## Why the still render is here too |
| 37 | |
| 38 | The CPU chaos game in this file is the reference for the bake, so it has to be the |
| 39 | RIGHT reference. Projecting those same CPU points through the engine's camera and |
| 40 | comparing them to what the GPU actually plots is what earns it: if the two agree, |
| 41 | the CPU model reproduces the kernel, and its verdict on the bake means something. |
| 42 | The camera is rotated 30 degrees about Y on purpose, so screen x mixes world x and |
| 43 | z and the check is sensitive to all three axes rather than two. |
| 44 | |
| 45 | Needs the GPU and openvdb, so it runs inside Blender: |
| 46 | blender -b --factory-startup -noaudio --python blender/test_bake_orientation.py |
| 47 | """ |
| 48 | |
| 49 | import math |
| 50 | import os |
| 51 | import sys |
| 52 | import tempfile |
| 53 | import unittest |
| 54 | |
| 55 | sys.path.insert(0, os.path.dirname(os.path.abspath(__file__))) |
| 56 | |
| 57 | # `ifs_volume` / `ifs_render` pull in `gpu`, which only exists inside Blender, so |
| 58 | # they are imported per-test rather than at module scope — see test_volume_axes. |
| 59 | |
| 60 | S = 0.5 |
| 61 | # Deliberately asymmetric in all three axes. Extents: x 1.2, y 0.8, z 0.4; |
| 62 | # centroid (0.15, 0.10, 0.05). No two axes share a span or a centre. |
| 63 | FIXED_POINTS = ( |
| 64 | (0.9, 0.0, 0.0), |
| 65 | (0.0, 0.6, 0.0), |
| 66 | (0.0, 0.0, 0.3), |
| 67 | (-0.3, -0.2, -0.1), |
| 68 | ) |
| 69 | # Non-cubic AND off-centre, so a per-axis mix-up distorts as well as rotates. |
| 70 | BOUNDS = ((-0.6, -0.5, -0.4), (1.2, 0.9, 0.6)) |
| 71 | |
| 72 | VOXELS = 64 |
| 73 | BAKE_DISPATCHES = 20 |
| 74 | BAKE_ITERS = 100 # keeps the fixed-point histogram far from uint32 wrap |
| 75 | |
| 76 | STILL_SIZE = 192 |
| 77 | STILL_DISPATCHES = 40 |
| 78 | STILL_ITERS = 200 |
| 79 | |
| 80 | CPU_POINTS = 120_000 |
| 81 | CPU_STEPS = 60 # 0.5**60 is far below float precision: fully converged |
| 82 | |
| 83 | # 30 degrees about Y, so the camera is oblique to every world axis. |
| 84 | CAM_YAW = 30.0 |
| 85 | CAM_DISTANCE = 3.0 |
| 86 | |
| 87 | |
| 88 | def _mean(points, axis): |
| 89 | return sum(p[axis] for p in points) / len(points) |
| 90 | |
| 91 | |
| 92 | def _xyz(v): |
| 93 | return {"X": float(v[0]), "Y": float(v[1]), "Z": float(v[2])} |
| 94 | |
| 95 | |
| 96 | def _camera(): |
| 97 | """Position and orientation, in the engine's convention. |
| 98 | |
| 99 | `view_projection` builds the view matrix looking along `position - forward`, |
| 100 | which inverts the usual sense: the points a camera SEES are the ones along |
| 101 | `+forward` from it (that is the documented negative-w quirk in ifs_render). |
| 102 | So framing the centroid at distance d means stepping BACK along forward. |
| 103 | """ |
| 104 | half = math.radians(CAM_YAW) * 0.5 |
| 105 | quat = (0.0, math.sin(half), 0.0, math.cos(half)) # (x, y, z, w) |
| 106 | forward = (math.sin(math.radians(CAM_YAW)), 0.0, math.cos(math.radians(CAM_YAW))) |
| 107 | centroid = tuple(_mean(FIXED_POINTS, a) for a in range(3)) |
| 108 | position = tuple(centroid[a] - CAM_DISTANCE * forward[a] for a in range(3)) |
| 109 | return position, quat |
| 110 | |
| 111 | |
| 112 | def _document_dict(): |
| 113 | position, quat = _camera() |
| 114 | iterators = [] |
| 115 | for index, fixed in enumerate(FIXED_POINTS): |
| 116 | iterators.append({ |
| 117 | "Transform": {"Item1": "Affine"}, |
| 118 | "Vec3Params": { |
| 119 | # Fixed point F of `p := T + S*p` is T/(1-S), so T = F*(1-S). |
| 120 | "Translate": _xyz(tuple(c * (1.0 - S) for c in fixed)), |
| 121 | "Scale": _xyz((S, S, S)), |
| 122 | "Rotate": _xyz((0.0, 0.0, 0.0)), |
| 123 | }, |
| 124 | "Opacity": 1.0, |
| 125 | "ColorIndex": index / (len(FIXED_POINTS) - 1.0), |
| 126 | "ColorSpeed": 0.5, |
| 127 | "BaseWeight": 1.0, |
| 128 | "StartWeight": 1.0, |
| 129 | }) |
| 130 | return { |
| 131 | "Title": "orientation-probe", |
| 132 | "Iterators": iterators, |
| 133 | "Camera": { |
| 134 | "Position": _xyz(position), |
| 135 | "Orientation": {"X": quat[0], "Y": quat[1], "Z": quat[2], "W": quat[3]}, |
| 136 | "FieldOfView": 60.0, |
| 137 | "FocusDistance": CAM_DISTANCE, |
| 138 | }, |
| 139 | # White, so the still render's brightness tracks density rather than hue. |
| 140 | "Palette": {"Colors": [{"X": 1.0, "Y": 1.0, "Z": 1.0, "W": 1.0}]}, |
| 141 | "Brightness": 4.0, |
| 142 | "Warmup": 20, |
| 143 | "Entropy": 0.01, |
| 144 | "ImageResolution": f"{STILL_SIZE}, {STILL_SIZE}", |
| 145 | } |
| 146 | |
| 147 | |
| 148 | def _cpu_attractor(seed=20260817): |
| 149 | """ |
| 150 | The same chaos game, on the CPU, in world coordinates. |
| 151 | |
| 152 | Independent of the GPU, the readback order, openvdb and the camera — which is |
| 153 | what makes it usable as ground truth for all of them. The RNG does not need to |
| 154 | match the kernel's: an IFS attractor is the invariant SET, so any sequence of |
| 155 | map choices converges to the same distribution. |
| 156 | """ |
| 157 | import numpy as np |
| 158 | rng = np.random.default_rng(seed) |
| 159 | fixed = np.array(FIXED_POINTS, dtype=np.float64) |
| 160 | points = rng.normal(size=(CPU_POINTS, 3)) |
| 161 | for _ in range(CPU_STEPS): |
| 162 | choice = rng.integers(0, len(fixed), size=CPU_POINTS) |
| 163 | points = S * points + (1.0 - S) * fixed[choice] |
| 164 | return points |
| 165 | |
| 166 | |
| 167 | def _read_vdb(path): |
| 168 | """Every lit voxel of a written `.vdb`, as (world position, density). |
| 169 | |
| 170 | Positions come through the grid's OWN transform, which is what Cycles reads — |
| 171 | not a re-derivation from BOUNDS, which would only re-state the assumption under |
| 172 | test. |
| 173 | |
| 174 | `indexToWorld(i)` is the voxel's CENTRE, not the corner of its cell. That was |
| 175 | measured, not assumed: a grid with one lit voxel at index (0,0,0) and voxel |
| 176 | size 1 imports into Blender with a bounding box starting at -0.5. So the world |
| 177 | position of a voxel's mass IS its indexToWorld, and no correction belongs here |
| 178 | — `accumulate_voxel` bins a sample into the cell it falls inside, so the grid |
| 179 | transform has to place that cell's CENTRE, which is what `write_vdb` does. |
| 180 | """ |
| 181 | import numpy as np |
| 182 | import openvdb |
| 183 | grid = openvdb.read(path, "density") |
| 184 | array = np.zeros((VOXELS, VOXELS, VOXELS), dtype=np.float32) |
| 185 | grid.copyToArray(array) |
| 186 | lit = np.argwhere(array > 0.0) |
| 187 | weights = array[lit[:, 0], lit[:, 1], lit[:, 2]].astype(np.float64) |
| 188 | world = np.array([grid.transform.indexToWorld(tuple(int(c) for c in idx)) |
| 189 | for idx in lit]) |
| 190 | return world, weights |
| 191 | |
| 192 | |
| 193 | def _weighted_stats(positions, weights): |
| 194 | """Centroid and per-axis mass percentiles, the two statistics compared.""" |
| 195 | import numpy as np |
| 196 | total = float(weights.sum()) |
| 197 | centroid = (positions * weights[:, None]).sum(axis=0) / total |
| 198 | percentiles = [] |
| 199 | for axis in range(3): |
| 200 | order = np.argsort(positions[:, axis]) |
| 201 | values = positions[order, axis] |
| 202 | cumulative = np.cumsum(weights[order]) / total |
| 203 | percentiles.append([float(values[np.searchsorted(cumulative, q)]) |
| 204 | for q in (0.05, 0.5, 0.95)]) |
| 205 | return centroid, percentiles |
| 206 | |
| 207 | |
| 208 | class Orientation(unittest.TestCase): |
| 209 | doc = None |
| 210 | grid = None |
| 211 | stats = None |
| 212 | |
| 213 | @classmethod |
| 214 | def setUpClass(cls): |
| 215 | try: |
| 216 | import numpy # noqa: F401 |
| 217 | import openvdb # noqa: F401 |
| 218 | import ifs_volume # noqa: F401 |
| 219 | except ImportError as ex: |
| 220 | raise unittest.SkipTest(f"needs Blender's gpu/openvdb modules ({ex})") |
| 221 | |
| 222 | import json |
| 223 | import ifs_document |
| 224 | import ifs_render |
| 225 | import ifs_volume |
| 226 | |
| 227 | cls.tmp = tempfile.mkdtemp() |
| 228 | cls.doc_path = os.path.join(cls.tmp, "orientation.ifsjson") |
| 229 | with open(cls.doc_path, "w", encoding="utf-8") as handle: |
| 230 | json.dump(_document_dict(), handle) |
| 231 | |
| 232 | library = ifs_render.load_library( |
| 233 | os.path.join(ifs_render.REPO, "Transforms"), |
| 234 | os.path.join(ifs_render.REPO, "WpfDisplay", "Library", "Includes")) |
| 235 | cls.doc = ifs_document.load(cls.doc_path, library) |
| 236 | cls.points = _cpu_attractor() |
| 237 | |
| 238 | # Bake with EXPLICIT bounds: this test is about orientation, so the box is |
| 239 | # held fixed rather than fitted. `test_auto_bounds_brackets_the_fractal` |
| 240 | # covers the fitted path separately. |
| 241 | cls.vdb_path = os.path.join(cls.tmp, "orientation.vdb") |
| 242 | _path, cls.stats = ifs_volume.bake( |
| 243 | cls.doc_path, cls.vdb_path, size=VOXELS, dispatches=BAKE_DISPATCHES, |
| 244 | bounds=BOUNDS, invocation_iters=BAKE_ITERS) |
| 245 | |
| 246 | # --------------------------------------------------------------- the guards |
| 247 | |
| 248 | def test_the_document_can_detect_a_swap(self): |
| 249 | """Guards the guard: a symmetric document would make a transpose invisible.""" |
| 250 | lo = [min(f[a] for f in FIXED_POINTS) for a in range(3)] |
| 251 | hi = [max(f[a] for f in FIXED_POINTS) for a in range(3)] |
| 252 | spans = [hi[a] - lo[a] for a in range(3)] |
| 253 | self.assertEqual(3, len(set(spans)), f"spans {spans} are not all distinct") |
| 254 | centres = [_mean(FIXED_POINTS, a) for a in range(3)] |
| 255 | self.assertEqual(3, len(set(centres)), |
| 256 | f"centroid {centres} cannot distinguish an axis swap") |
| 257 | for axis in range(3): |
| 258 | self.assertNotAlmostEqual( |
| 259 | lo[axis], -hi[axis], places=6, |
| 260 | msg=f"axis {axis} is symmetric about zero, so a mirror is invisible") |
| 261 | |
| 262 | def test_the_cpu_reference_reproduces_the_analytic_attractor(self): |
| 263 | """The CPU chaos game must land on the tetrahedron, or it is no reference.""" |
| 264 | import numpy as np |
| 265 | expected = np.array([_mean(FIXED_POINTS, a) for a in range(3)]) |
| 266 | measured = self.points.mean(axis=0) |
| 267 | for axis in range(3): |
| 268 | self.assertAlmostEqual( |
| 269 | expected[axis], measured[axis], places=2, |
| 270 | msg=f"axis {axis}: CPU centroid {measured[axis]:.4f} is not the " |
| 271 | f"analytic {expected[axis]:.4f}") |
| 272 | # The hull is the tetrahedron, so nothing may sit outside it by more than |
| 273 | # float slop. |
| 274 | for axis in range(3): |
| 275 | self.assertGreaterEqual( |
| 276 | self.points[:, axis].min(), |
| 277 | min(f[axis] for f in FIXED_POINTS) - 1e-6) |
| 278 | self.assertLessEqual( |
| 279 | self.points[:, axis].max(), |
| 280 | max(f[axis] for f in FIXED_POINTS) + 1e-6) |
| 281 | |
| 282 | # ------------------------------------------------- leg 1: is the bake right |
| 283 | |
| 284 | def test_the_bake_lands_where_the_fractal_is(self): |
| 285 | """ |
| 286 | The centroid of the baked density, in world coordinates, against the |
| 287 | analytic centroid. An x/z transpose moves it by 0.10 world units — more |
| 288 | than three voxels — so this cannot pass on a transposed bake. |
| 289 | |
| 290 | The tolerance is a QUARTER of a voxel, which is placement accuracy rather |
| 291 | than merely "not transposed". It can afford to be: the residual measured on |
| 292 | a correct bake is 0.01-0.03 voxels, so this still has an order of magnitude |
| 293 | of headroom, and a half-voxel convention slip — the grid transform placing |
| 294 | cell CORNERS where centres belong — shows up as a clean -0.5 on all three |
| 295 | axes rather than passing quietly inside a loose bound. |
| 296 | """ |
| 297 | import numpy as np |
| 298 | positions, weights = _read_vdb(self.vdb_path) |
| 299 | centroid, _pct = _weighted_stats(positions, weights) |
| 300 | expected = np.array([_mean(FIXED_POINTS, a) for a in range(3)]) |
| 301 | voxel = np.array([(BOUNDS[1][a] - BOUNDS[0][a]) / VOXELS for a in range(3)]) |
| 302 | offset = (centroid - expected) / voxel |
| 303 | print(f"[bake] centroid {np.round(centroid, 4)} expected " |
| 304 | f"{np.round(expected, 4)} offset {np.round(offset, 3)} voxels") |
| 305 | for axis in range(3): |
| 306 | self.assertLess( |
| 307 | abs(offset[axis]), 0.25, |
| 308 | f"axis {axis}: baked centroid {centroid[axis]:.4f} is " |
| 309 | f"{offset[axis]:.2f} voxels from the analytic {expected[axis]:.4f} " |
| 310 | f"— the bake does not sit where the fractal is") |
| 311 | |
| 312 | def test_the_bake_has_the_right_shape_on_every_axis(self): |
| 313 | """ |
| 314 | A centroid alone cannot see a swap of two axes with equal means, so compare |
| 315 | the mass distribution per axis: the 5th, 50th and 95th percentiles of the |
| 316 | baked density against the same percentiles of the CPU attractor. |
| 317 | """ |
| 318 | import numpy as np |
| 319 | positions, weights = _read_vdb(self.vdb_path) |
| 320 | _centroid, baked = _weighted_stats(positions, weights) |
| 321 | reference = _weighted_stats( |
| 322 | self.points, np.ones(len(self.points)))[1] |
| 323 | voxel = [(BOUNDS[1][a] - BOUNDS[0][a]) / VOXELS for a in range(3)] |
| 324 | for axis in range(3): |
| 325 | for label, got, want in zip(("p05", "p50", "p95"), |
| 326 | baked[axis], reference[axis]): |
| 327 | drift = (got - want) / voxel[axis] |
| 328 | print(f"[bake] axis {axis} {label}: {got:+.4f} vs {want:+.4f} " |
| 329 | f"({drift:+.2f} voxels)") |
| 330 | self.assertLess( |
| 331 | abs(drift), 2.0, |
| 332 | f"axis {axis} {label}: baked {got:.4f} against the CPU " |
| 333 | f"attractor's {want:.4f}, {drift:.2f} voxels apart") |
| 334 | |
| 335 | def test_auto_bounds_brackets_the_fractal(self): |
| 336 | """ |
| 337 | The fitted box is a second consumer of the readback order (`_mass_bounds` |
| 338 | reorders array axes into world xyz), so a wrong order fits a transposed box |
| 339 | while the bake inside it still looks like a fractal. |
| 340 | """ |
| 341 | import numpy as np |
| 342 | import ifs_volume |
| 343 | shader, _library = ifs_volume.build_shader() |
| 344 | ubos = ifs_volume._document_buffers(self.doc) |
| 345 | alias = ifs_volume.ifs_host.build_alias_tables(self.doc.iterators) |
| 346 | n = len(self.doc.iterators) |
| 347 | import gpu |
| 348 | alias_tex = gpu.types.GPUTexture( |
| 349 | (n, n), format='RGBA32F', |
| 350 | data=gpu.types.Buffer('FLOAT', len(alias), alias)) |
| 351 | lo, hi, _fit = ifs_volume.auto_bounds(shader, self.doc, ubos, alias_tex) |
| 352 | |
| 353 | body_lo = np.percentile(self.points, 1.0, axis=0) |
| 354 | body_hi = np.percentile(self.points, 99.0, axis=0) |
| 355 | span = body_hi - body_lo |
| 356 | print(f"[auto_bounds] lo {np.round(lo, 4)} hi {np.round(hi, 4)} " |
| 357 | f"vs body {np.round(body_lo, 4)}..{np.round(body_hi, 4)}") |
| 358 | for axis in range(3): |
| 359 | self.assertLess( |
| 360 | lo[axis], body_lo[axis], |
| 361 | f"axis {axis}: fitted lo {lo[axis]:.4f} clips the fractal, whose " |
| 362 | f"1st percentile is {body_lo[axis]:.4f}") |
| 363 | self.assertGreater( |
| 364 | hi[axis], body_hi[axis], |
| 365 | f"axis {axis}: fitted hi {hi[axis]:.4f} clips the fractal, whose " |
| 366 | f"99th percentile is {body_hi[axis]:.4f}") |
| 367 | # Generous, but a transposed fit on this document is off by 2-3x. |
| 368 | self.assertLess( |
| 369 | (hi[axis] - lo[axis]), 3.0 * span[axis], |
| 370 | f"axis {axis}: fitted span {(hi[axis] - lo[axis]):.4f} is far wider " |
| 371 | f"than the fractal's {span[axis]:.4f} on that axis") |
| 372 | |
| 373 | # ------------------------------- leg 2: is the CPU reference the right one |
| 374 | |
| 375 | def _still_density(self): |
| 376 | """ |
| 377 | The still path's dispatch loop, for the raw density it tonemaps away. |
| 378 | |
| 379 | `render()` returns an image datablock, and the tonemap plus the 8-bit |
| 380 | buffer destroy exactly the faint tail this comparison needs. The SHADER is |
| 381 | the same one `render()` builds, so nothing about the projection is |
| 382 | reproduced here — only the loop around it. |
| 383 | """ |
| 384 | import gpu |
| 385 | import numpy as np |
| 386 | import ifs_host |
| 387 | import ifs_render |
| 388 | |
| 389 | camera = self.doc.camera |
| 390 | camera.view_proj = ifs_render.view_projection( |
| 391 | camera.position[:3], camera.forward[:3], camera.up, |
| 392 | fov_degrees=camera.field_of_view) |
| 393 | |
| 394 | kernel, drain, _library = ifs_render._build_shaders( |
| 395 | STILL_SIZE, STILL_SIZE, ifs_render.DEFAULT_HIST_SCALE) |
| 396 | ubos = { |
| 397 | 'settings': gpu.types.GPUUniformBuf( |
| 398 | ifs_host.pack_settings(camera, self.doc.settings)), |
| 399 | 'iterators_buf': gpu.types.GPUUniformBuf( |
| 400 | ifs_host.pack_iterators(self.doc.iterators)), |
| 401 | 'palette_buf': gpu.types.GPUUniformBuf( |
| 402 | ifs_host.pack_palette( |
| 403 | self.doc.palette(self.doc.settings.gradient_resolution))), |
| 404 | 'real_params_buf': gpu.types.GPUUniformBuf( |
| 405 | ifs_host.pack_real_params( |
| 406 | ifs_host.flatten_real_params(self.doc.iterators))), |
| 407 | 'vec3_params_buf': gpu.types.GPUUniformBuf( |
| 408 | ifs_host.pack_vec3_params( |
| 409 | ifs_host.flatten_vec3_params(self.doc.iterators))), |
| 410 | } |
| 411 | alias = ifs_host.build_alias_tables(self.doc.iterators) |
| 412 | n = len(self.doc.iterators) |
| 413 | alias_tex = gpu.types.GPUTexture( |
| 414 | (n, n), format='RGBA32F', |
| 415 | data=gpu.types.Buffer('FLOAT', len(alias), alias)) |
| 416 | |
| 417 | histogram = gpu.types.GPUTexture((STILL_SIZE * 4, STILL_SIZE), format='R32UI') |
| 418 | histogram.clear(format='UINT', value=(0, 0, 0, 0)) |
| 419 | accum = gpu.types.GPUTexture((STILL_SIZE, STILL_SIZE), format='RGBA32F') |
| 420 | accum.clear(format='FLOAT', value=(0.0, 0.0, 0.0, 0.0)) |
| 421 | rows = (ifs_render.INVOCATIONS + ifs_render.STATE_W - 1) // ifs_render.STATE_W |
| 422 | state = gpu.types.GPUTexture( |
| 423 | (ifs_render.STATE_W * 2, rows), format='RGBA32F') |
| 424 | state.clear(format='FLOAT', value=(0.0, 0.0, 0.0, 0.0)) |
| 425 | |
| 426 | for index in range(STILL_DISPATCHES): |
| 427 | frame = gpu.types.GPUUniformBuf(ifs_host.pack_frame( |
| 428 | seed=0x9E3779B9, width=STILL_SIZE, height=STILL_SIZE, |
| 429 | dispatch_cnt=index, reset_points_state=1 if index == 0 else 0, |
| 430 | invocation_iters=STILL_ITERS)) |
| 431 | kernel.bind() |
| 432 | kernel.image('histogram_img', histogram) |
| 433 | kernel.image('state_img', state) |
| 434 | kernel.uniform_sampler('alias_tables_tex', alias_tex) |
| 435 | kernel.uniform_block('frame', frame) |
| 436 | for name, ubo in ubos.items(): |
| 437 | kernel.uniform_block(name, ubo) |
| 438 | gpu.compute.dispatch(kernel, ifs_render.WORKGROUPS, 1, 1) |
| 439 | |
| 440 | drain.bind() |
| 441 | drain.image('histogram_img', histogram) |
| 442 | drain.image('accum_img', accum) |
| 443 | drain.uniform_int('out_size', (STILL_SIZE, STILL_SIZE)) |
| 444 | gpu.compute.dispatch(drain, (STILL_SIZE + 7) // 8, (STILL_SIZE + 7) // 8, 1) |
| 445 | |
| 446 | pixels = np.array(accum.read(), dtype=np.float32).reshape( |
| 447 | STILL_SIZE, STILL_SIZE, 4) |
| 448 | return pixels[:, :, 3] |
| 449 | |
| 450 | def _projected_density(self): |
| 451 | """The CPU attractor through the engine's own projection, same binning.""" |
| 452 | import numpy as np |
| 453 | import ifs_render |
| 454 | view_proj = self.doc.camera.view_proj |
| 455 | out = np.zeros((STILL_SIZE, STILL_SIZE), dtype=np.float64) |
| 456 | for point in self.points: |
| 457 | result = ifs_render.project_point( |
| 458 | view_proj, (float(point[0]), float(point[1]), float(point[2])), |
| 459 | STILL_SIZE, STILL_SIZE) |
| 460 | if result is None: |
| 461 | continue |
| 462 | x, y, _w = result |
| 463 | out[int(round(y)), int(round(x))] += 1.0 |
| 464 | return out |
| 465 | |
| 466 | def test_the_still_render_matches_the_cpu_reference(self): |
| 467 | """ |
| 468 | Earns the reference. If the CPU chaos game and the engine's own projection |
| 469 | put the fractal in the same pixels, the CPU model reproduces the kernel — |
| 470 | which is what makes its verdict on the bake, above, mean anything. |
| 471 | """ |
| 472 | import numpy as np |
| 473 | rendered = self._still_density() |
| 474 | predicted = self._projected_density() |
| 475 | self.assertGreater(float(rendered.sum()), 0.0, "the still render is empty") |
| 476 | self.assertGreater(float(predicted.sum()), 0.0, |
| 477 | "the CPU points all projected off screen") |
| 478 | |
| 479 | def centroid(image): |
| 480 | ys, xs = np.nonzero(image) |
| 481 | weights = image[ys, xs] |
| 482 | total = weights.sum() |
| 483 | return np.array([(xs * weights).sum() / total, |
| 484 | (ys * weights).sum() / total]) |
| 485 | |
| 486 | got, want = centroid(rendered), centroid(predicted) |
| 487 | print(f"[still] centroid {np.round(got, 2)} vs projected CPU " |
| 488 | f"{np.round(want, 2)} px") |
| 489 | for axis, name in enumerate(("x", "y")): |
| 490 | self.assertLess( |
| 491 | abs(got[axis] - want[axis]), 4.0, |
| 492 | f"screen {name}: the render's centroid is {got[axis]:.1f} px and " |
| 493 | f"the projected CPU attractor's is {want[axis]:.1f} px") |
| 494 | |
| 495 | # Shape, not just position: a mirror about the centroid leaves it put. |
| 496 | overlap = np.logical_and(rendered > 0, predicted > 0).sum() |
| 497 | union = np.logical_or(rendered > 0, predicted > 0).sum() |
| 498 | iou = overlap / float(union) |
| 499 | print(f"[still] lit-pixel IoU {iou:.3f}") |
| 500 | self.assertGreater( |
| 501 | iou, 0.5, |
| 502 | f"only {iou:.1%} of lit pixels agree — the render and the CPU " |
| 503 | "attractor are not the same shape in the same place") |
| 504 | |
| 505 | # ------------------------------------------------------------- the control |
| 506 | |
| 507 | def test_a_transposed_bake_would_be_caught(self): |
| 508 | """ |
| 509 | The control. Feeds `write_vdb` an x/z-swapped grid and asserts the centroid |
| 510 | check REJECTS it — otherwise a green suite proves only that the assertions |
| 511 | are loose enough to pass anything. |
| 512 | """ |
| 513 | import numpy as np |
| 514 | import openvdb |
| 515 | import ifs_volume |
| 516 | |
| 517 | # Reconstruct an accumulator from the good bake, then swap x and z. |
| 518 | grid = openvdb.read(self.vdb_path, "density") |
| 519 | good = np.zeros((VOXELS, VOXELS, VOXELS), dtype=np.float32) |
| 520 | grid.copyToArray(good) |
| 521 | # `write_vdb` expects [z][y][x] and transposes; handing it the array that |
| 522 | # is ALREADY [x][y][z] reproduces exactly the bug this file guards. |
| 523 | swapped = np.zeros((VOXELS, VOXELS, VOXELS, 4), dtype=np.float32) |
| 524 | swapped[:, :, :, 3] = good |
| 525 | swapped[:, :, :, :3] = good[:, :, :, None] |
| 526 | path = os.path.join(self.tmp, "transposed.vdb") |
| 527 | ifs_volume.write_vdb(swapped, BOUNDS, path) |
| 528 | |
| 529 | bad = openvdb.read(path, "density") |
| 530 | array = np.zeros((VOXELS, VOXELS, VOXELS), dtype=np.float32) |
| 531 | bad.copyToArray(array) |
| 532 | lit = np.argwhere(array > 0.0) |
| 533 | weights = array[lit[:, 0], lit[:, 1], lit[:, 2]].astype(np.float64) |
| 534 | world = np.array([bad.transform.indexToWorld(tuple(int(c) for c in idx)) |
| 535 | for idx in lit]) |
| 536 | centroid, _pct = _weighted_stats(world, weights) |
| 537 | expected = np.array([_mean(FIXED_POINTS, a) for a in range(3)]) |
| 538 | voxel = np.array([(BOUNDS[1][a] - BOUNDS[0][a]) / VOXELS for a in range(3)]) |
| 539 | offset = np.abs(centroid - expected) / voxel |
| 540 | print(f"[control] transposed centroid {np.round(centroid, 4)} " |
| 541 | f"offset {np.round(offset, 2)} voxels") |
| 542 | self.assertGreater( |
| 543 | offset.max(), 1.5, |
| 544 | "a transposed bake landed within tolerance of the correct centroid — " |
| 545 | "the orientation assertions above cannot detect a transpose") |
| 546 | |
| 547 | |
| 548 | if __name__ == "__main__": |
| 549 | result = unittest.main(argv=sys.argv[:1], exit=False, verbosity=2).result |
| 550 | sys.exit(0 if result.wasSuccessful() else 1) |