| … | … | @@ -358,13 +358,30 @@ void accumulate_voxel(vec3 pos, vec4 color) |
| 358 | 358 | \tvec3 t = (pos - volume_min) / span; |
| 359 | 359 | \tif (any(lessThan(t, vec3(0.0))) || any(greaterThanEqual(t, vec3(1.0)))) |
| 360 | 360 | \t\treturn; // outside the bake bounds; the chaos game is unbounded |
| 361 | \timageAtomicAdd(volume_img, ivec3(t * vec3(volume_size)), | |
| 361 | \t// Stored TRANSPOSED, deliberately: voxel (x,y,z) goes to texel (z,y,x). | |
| 362 | \t// | |
| 363 | \t// GPUTexture.read() hands back a 3D texture as [z][y][x], while openvdb's | |
| 364 | \t// copyFromArray reads index (i,j,k) as world (x,y,z) — so the host used to fix | |
| 365 | \t// this up with `ascontiguousarray(density.transpose(2,1,0))`. That is a full | |
| 366 | \t// strided copy of the whole grid with cache-hostile access: measured 17.6s at | |
| 367 | \t// N=512, about 70% of the entire write, against 0.19s for copyFromArray itself. | |
| 368 | \t// Permuting the index HERE is free — it is three registers — and the readback | |
| 369 | \t// then arrives in exactly the order openvdb wants. | |
| 370 | \t// | |
| 371 | \t// (Passing the transposed VIEW to copyFromArray instead does NOT work: it reads | |
| 372 | \t// the underlying buffer in memory order and silently ignores the strides, which | |
| 373 | \t// writes the fractal transposed. test_volume_axes and test_bake_orientation both | |
| 374 | \t// catch it; activeVoxelCount does not, being invariant under permutation.) | |
| 375 | \tivec3 v = ivec3(t * vec3(volume_size)); | |
| 376 | \timageAtomicAdd(volume_img, ivec3(v.z, v.y, v.x), | |
| 362 | 377 | \t uint(max(color.w, 0.0) * HIST_SCALE)); |
| 363 | 378 | |
| 364 | 379 | \t// Three contiguous texels per colour voxel, the same coalescing argument the |
| 365 | 380 | \t// old four-channel layout made. |
| 381 | \t// Same permutation as density, and the channel triple stays on the FIRST axis | |
| 382 | \t// so a colour voxel is still three contiguous texels. | |
| 366 | 383 | \tivec3 c = ivec3(t * vec3(volume_size >> 1)); |
| 367 | \tivec3 cbase = ivec3(c.x * 3, c.y, c.z); | |
| 384 | \tivec3 cbase = ivec3(c.z * 3, c.y, c.x); | |
| 368 | 385 | \timageAtomicAdd(colour_img, cbase + ivec3(0, 0, 0), uint(max(color.r, 0.0) * HIST_SCALE)); |
| 369 | 386 | \timageAtomicAdd(colour_img, cbase + ivec3(1, 0, 0), uint(max(color.g, 0.0) * HIST_SCALE)); |
| 370 | 387 | \timageAtomicAdd(colour_img, cbase + ivec3(2, 0, 0), uint(max(color.b, 0.0) * HIST_SCALE)); |
| … | … | @@ -215,8 +215,10 @@ def auto_bounds(shader, doc, ubos, alias_tex, probe_size=64, dispatches=8, |
| 215 | 215 | # toward a zero-size voxel, so this stops when an axis is thinner than one voxel |
| 216 | 216 | # of the LARGEST axis: past that it is planar in any sense the bake can represent. |
| 217 | 217 | for _ in range(refine_passes): |
| 218 | # `axis` is a world axis and the array is [x][y][z], so they coincide; this | |
| 219 | # used to be `2 - axis` to undo the old [z][y][x] readback. | |
| 218 | 220 | occupied = [int((density.sum(axis=tuple(a for a in range(3) |
| 219 | if a != 2 - axis)) > 0).sum()) for axis in range(3)] | |
| 221 | if a != axis)) > 0).sum()) for axis in range(3)] | |
| 220 | 222 | span = hi_w - lo_w |
| 221 | 223 | floor = float(span.max()) / probe_size |
| 222 | 224 | collapsed = [a for a in range(3) |
| … | … | @@ -261,7 +263,9 @@ def _mass_bounds(density, probe_size, percentile): |
| 261 | 263 | where the mass is instead, which is what the bake should spend its resolution on. |
| 262 | 264 | """ |
| 263 | 265 | lo_idx, hi_idx = [], [] |
| 264 | for array_axis in (2, 1, 0): # array is [z][y][x]; iterate world x, y, z | |
| 266 | # The array is [x][y][z] now (the kernel stores transposed), so world axis and | |
| 267 | # array axis coincide — this used to iterate (2, 1, 0) to undo [z][y][x]. | |
| 268 | for array_axis in (0, 1, 2): | |
| 265 | 269 | others = tuple(a for a in range(3) if a != array_axis) |
| 266 | 270 | marginal = density.sum(axis=others) |
| 267 | 271 | total = float(marginal.sum()) |
| … | … | @@ -371,11 +375,16 @@ def write_vdb(density_grid: np.ndarray, colour_grid: np.ndarray, bounds, |
| 371 | 375 | peak = float(density.max()) |
| 372 | 376 | if peak <= 0: |
| 373 | 377 | raise RuntimeError("nothing was accumulated - the bounds may exclude the fractal") |
| 374 | # Normalise IN PLACE on a single transposed copy rather than building one | |
| 375 | # array per step. The old chain (divide -> power -> transpose -> contiguous) | |
| 376 | # held up to four full-size temporaries at once; at N=384 each is 226 MB. | |
| 377 | density = np.ascontiguousarray(density.transpose(2, 1, 0), dtype=np.float32) | |
| 378 | density /= peak | |
| 378 | # NO transpose. The kernel now stores voxel (x,y,z) at texel (z,y,x), so the | |
| 379 | # readback already arrives as [x][y][z] — exactly what copyFromArray wants. | |
| 380 | # Doing it here instead cost 17.6s at N=512, ~70% of the whole write, because | |
| 381 | # a (N,N,N) transpose is a full strided copy with cache-hostile access. | |
| 382 | # | |
| 383 | # Still a COPY, not in place: `density_grid` is read again below for the colour | |
| 384 | # blocks and again by bake_document for its occupancy and exposure statistics, | |
| 385 | # so normalising the caller's array would silently double-normalise. A | |
| 386 | # contiguous divide is ~0.2s; the strided copy was the expensive part. | |
| 387 | density = density_grid / peak | |
| 379 | 388 | if gamma and gamma != 1.0: |
| 380 | 389 | np.power(density, 1.0 / gamma, out=density) |
| 381 | 390 | |
| … | … | @@ -397,11 +406,10 @@ def write_vdb(density_grid: np.ndarray, colour_grid: np.ndarray, bounds, |
| 397 | 406 | half, 2, half, 2, half, 2).sum(axis=(1, 3, 5)) |
| 398 | 407 | safe = np.maximum(block[:, :, :, None], 1e-9) |
| 399 | 408 | |
| 400 | # Same transpose as density - the spatial axes only, leaving the component axis | |
| 401 | # last. A colour grid transposed differently from density would tint the wrong | |
| 402 | # parts of the fractal, which is harder to spot than a whole-volume rotation. | |
| 403 | colour = np.ascontiguousarray( | |
| 404 | np.clip(colour_grid / safe, 0.0, 1.0).transpose(2, 1, 0, 3)) | |
| 409 | # No transpose here either, for the same reason — and the two grids must agree, | |
| 410 | # because a colour grid oriented differently from density tints the wrong parts | |
| 411 | # of the fractal, which is far harder to spot than a whole-volume rotation. | |
| 412 | colour = np.clip(colour_grid / safe, 0.0, 1.0) | |
| 405 | 413 | |
| 406 | 414 | try: |
| 407 | 415 | cgrid = openvdb.Vec3SGrid() |
| … | … | @@ -437,10 +437,16 @@ class Orientation(unittest.TestCase): |
| 437 | 437 | grid.copyToArray(good) |
| 438 | 438 | # `write_vdb` expects [z][y][x] and transposes; handing it the array that |
| 439 | 439 | # is ALREADY [x][y][z] reproduces exactly the bug this file guards. |
| 440 | # Reproduce the bug by handing write_vdb data in the OLD [z][y][x] order. | |
| 441 | # It no longer transposes — the kernel stores voxel (x,y,z) at texel (z,y,x) | |
| 442 | # so the readback already arrives [x][y][z] — which means feeding it a | |
| 443 | # transposed array is now what produces a transposed bake. Contiguous, not a | |
| 444 | # view: copyFromArray reads memory order and ignores strides. | |
| 440 | 445 | half = VOXELS // 2 |
| 441 | 446 | swapped_colour = np.zeros((half, half, half, 3), dtype=np.float32) |
| 442 | 447 | path = os.path.join(self.tmp, "transposed.vdb") |
| 443 | ifs_volume.write_vdb(good, swapped_colour, BOUNDS, path) | |
| 448 | ifs_volume.write_vdb(np.ascontiguousarray(good.transpose(2, 1, 0)), | |
| 449 | swapped_colour, BOUNDS, path) | |
| 444 | 450 | |
| 445 | 451 | bad = openvdb.read(path, "density") |
| 446 | 452 | array = np.zeros((VOXELS, VOXELS, VOXELS), dtype=np.float32) |
| … | … | @@ -44,7 +44,12 @@ NONCUBIC = ((-1.0, -2.0, -3.0), (1.0, 2.0, 3.0)) |
| 44 | 44 | |
| 45 | 45 | def _grid_with_marker(): |
| 46 | 46 | """ |
| 47 | The pair `accumulate` returns: full-res density and HALF-res colour. | |
| 47 | The pair `accumulate` returns: full-res density and HALF-res colour, indexed | |
| 48 | [x][y][z]. | |
| 49 | ||
| 50 | The kernel stores voxel (x,y,z) at texel (z,y,x) precisely so the readback | |
| 51 | arrives in this order and needs no transpose — doing it host-side cost 17.6s | |
| 52 | at N=512. So the fixture is built [x][y][z] too. | |
| 48 | 53 | |
| 49 | 54 | Colour is premultiplied by density, as the kernel accumulates it, so the |
| 50 | 55 | marker's colour voxel carries density x hue and write_vdb divides it back out |
| … | … | @@ -52,9 +57,9 @@ def _grid_with_marker(): |
| 52 | 57 | """ |
| 53 | 58 | import numpy as np |
| 54 | 59 | density = np.zeros((S, S, S), dtype=np.float32) |
| 55 | density[MARK_Z, MARK_Y, MARK_X] = 1.0 | |
| 60 | density[MARK_X, MARK_Y, MARK_Z] = 1.0 | |
| 56 | 61 | colour = np.zeros((S // 2, S // 2, S // 2, 3), dtype=np.float32) |
| 57 | colour[MARK_Z // 2, MARK_Y // 2, MARK_X // 2] = (1.0, 0.0, 0.0) | |
| 62 | colour[MARK_X // 2, MARK_Y // 2, MARK_Z // 2] = (1.0, 0.0, 0.0) | |
| 58 | 63 | return density, colour |
| 59 | 64 | |
| 60 | 65 | |
| … | … | @@ -91,9 +96,10 @@ class BakeAxes(unittest.TestCase): |
| 91 | 96 | |
| 92 | 97 | self.assertGreater( |
| 93 | 98 | correct, 0.5, |
| 94 | "the lit voxel is not at (x, y, z) = (%d, %d, %d); the array is indexed " | |
| 95 | "[z][y][x] and copyFromArray reads (i, j, k) as (x, y, z), so it needs " | |
| 96 | "transposing before the write" % (MARK_X, MARK_Y, MARK_Z)) | |
| 99 | "the lit voxel is not at (x, y, z) = (%d, %d, %d); the kernel stores " | |
| 100 | "voxel (x,y,z) at texel (z,y,x) so the readback is [x][y][z], which is " | |
| 101 | "what copyFromArray reads as (x, y, z) — no transpose should be needed" | |
| 102 | % (MARK_X, MARK_Y, MARK_Z)) | |
| 97 | 103 | self.assertLess( |
| 98 | 104 | swapped, 0.5, |
| 99 | 105 | "the lit voxel is at the x/z-SWAPPED position — the bake is transposed") |