| … | … | @@ -123,7 +123,15 @@ def accumulate(shader, doc, ubos, alias_tex, size, bounds, dispatches, |
| 123 | 123 | |
| 124 | 124 | raw = np.array(volume.read(), dtype=np.uint32) |
| 125 | 125 | # (z, y, 4x) as read back; split the interleaved channels out. |
| 126 | | grid = raw.reshape(size, size, size, 4).astype(np.float32) / HIST_SCALE |
| 126 | # |
| 127 | # ONE allocation, not two. `raw.astype(np.float32) / HIST_SCALE` materialises |
| 128 | # the float copy AND the quotient — at N=384 that is two 906 MB arrays where |
| 129 | # one will do, and measured peak RSS was 2.4x the grid size because of chains |
| 130 | # like this. np.divide with an explicit dtype converts and scales in a single |
| 131 | # pass. |
| 132 | grid = np.divide(raw.reshape(size, size, size, 4), HIST_SCALE, |
| 133 | dtype=np.float32) |
| 134 | del raw |
| 127 | 135 | return grid |
| 128 | 136 | |
| 129 | 137 | |
| … | … | @@ -353,11 +361,13 @@ def write_vdb(grid: np.ndarray, bounds, path: str, name: str = "density", |
| 353 | 361 | peak = float(density.max()) |
| 354 | 362 | if peak <= 0: |
| 355 | 363 | raise RuntimeError("nothing was accumulated - the bounds may exclude the fractal") |
| 356 | | normalised = density / peak |
| 364 | # Normalise IN PLACE on a single transposed copy rather than building one |
| 365 | # array per step. The old chain (divide -> power -> transpose -> contiguous) |
| 366 | # held up to four full-size temporaries at once; at N=384 each is 226 MB. |
| 367 | density = np.ascontiguousarray(density.transpose(2, 1, 0), dtype=np.float32) |
| 368 | density /= peak |
| 357 | 369 | if gamma and gamma != 1.0: |
| 358 | | normalised = np.power(normalised, 1.0 / gamma) |
| 359 | | density = np.ascontiguousarray( |
| 360 | | normalised.transpose(2, 1, 0), dtype=np.float32) |
| 370 | np.power(density, 1.0 / gamma, out=density) |
| 361 | 371 | |
| 362 | 372 | grids = [] |
| 363 | 373 | dgrid = openvdb.FloatGrid() |
| … | … | @@ -367,13 +377,21 @@ def write_vdb(grid: np.ndarray, bounds, path: str, name: str = "density", |
| 367 | 377 | grids.append(dgrid) |
| 368 | 378 | |
| 369 | 379 | # Colour is accumulated weighted by density, so divide it back out. |
| 370 | | safe = np.maximum(grid[:, :, :, 3:4], 1e-9) |
| 371 | 380 | # Same transpose as density - the spatial axes only, leaving the component axis |
| 372 | 381 | # last. A colour grid transposed differently from density would tint the wrong |
| 373 | 382 | # parts of the fractal, which is harder to spot than a whole-volume rotation. |
| 374 | | colour = np.ascontiguousarray( |
| 375 | | np.clip(grid[:, :, :, :3] / safe, 0.0, 1.0) |
| 376 | | .transpose(2, 1, 0, 3).astype(np.float32)) |
| 383 | # |
| 384 | # Built per SLAB. The whole-array form (divide -> clip -> transpose -> astype) |
| 385 | # peaked at four simultaneous copies of a 3-channel volume; slabs bound that to |
| 386 | # one plane at a time while producing byte-identical output. |
| 387 | colour = np.empty((size, size, size, 3), dtype=np.float32) |
| 388 | for z in range(size): |
| 389 | plane = grid[z] # (y, x, 4) |
| 390 | safe = np.maximum(plane[:, :, 3:4], 1e-9) |
| 391 | tinted = np.clip(plane[:, :, :3] / safe, 0.0, 1.0) |
| 392 | # grid is [z][y][x]; the target is [x][y][z], so this plane becomes a |
| 393 | # column across the first axis. |
| 394 | colour[:, :, z, :] = tinted.transpose(1, 0, 2) |
| 377 | 395 | try: |
| 378 | 396 | cgrid = openvdb.Vec3SGrid() |
| 379 | 397 | cgrid.copyFromArray(colour) |