"""
Proves the bake is an IDENTITY map in world space, end to end, against a fractal
whose position is known analytically.
`test_volume_axes.py` proves `write_vdb`'s conventions in ISOLATION: it hands the
writer a synthetic `[z][y][x]` array and checks the marker comes out at the right
world position. That is necessary and not sufficient, because it ASSUMES the array
`accumulate` returns is indexed `[z][y][x]` in the first place. That assumption was
measured once by hand and never asserted, and everything downstream — the transpose
in `write_vdb`, the axis reorder in `_mass_bounds` — is built on it. If the readback
order were different, both would be consistently wrong and both tests would pass.
So this starts at the GPU and finishes in world coordinates:
document -> chaos game on the GPU -> voxel grid -> readback -> .vdb -> world
and compares the result to where the fractal provably is.
## The ground truth
The document is four Affine iterators, `p := 0.5*p + 0.5*F_i`, so the attractor is
the Sierpinski tetrahedron with vertices `F_i` — and TWO facts about it are exact,
not estimated:
* its **convex hull** is the tetrahedron, so the per-axis extent is just the range
of the `F_i` components;
* its **centroid** is the mean of the `F_i`. The invariant measure satisfies
`mu = E[0.5*p + 0.5*F_i] = 0.5*mu + 0.5*mean(F)`, hence `mu = mean(F)`.
The vertices are chosen so every axis is distinguishable from every other — three
different spans, three different centroid components, none symmetric about zero.
A transposed, mirrored or displaced bake therefore cannot land on the right answer
by coincidence, which is exactly the failure mode a symmetric test document hides.
`test_the_document_can_detect_a_swap` guards that property.
## Why the still render is here too
The CPU chaos game in this file is the reference for the bake, so it has to be the
RIGHT reference. Projecting those same CPU points through the engine's camera and
comparing them to what the GPU actually plots is what earns it: if the two agree,
the CPU model reproduces the kernel, and its verdict on the bake means something.
The camera is rotated 30 degrees about Y on purpose, so screen x mixes world x and
z and the check is sensitive to all three axes rather than two.
Needs the GPU and openvdb, so it runs inside Blender:
blender -b --factory-startup -noaudio --python blender/test_bake_orientation.py
"""
import math
import os
import sys
import tempfile
import unittest
sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
# `ifs_volume` / `ifs_render` pull in `gpu`, which only exists inside Blender, so
# they are imported per-test rather than at module scope — see test_volume_axes.
S = 0.5
# Deliberately asymmetric in all three axes. Extents: x 1.2, y 0.8, z 0.4;
# centroid (0.15, 0.10, 0.05). No two axes share a span or a centre.
FIXED_POINTS = (
(0.9, 0.0, 0.0),
(0.0, 0.6, 0.0),
(0.0, 0.0, 0.3),
(-0.3, -0.2, -0.1),
)
# Non-cubic AND off-centre, so a per-axis mix-up distorts as well as rotates.
BOUNDS = ((-0.6, -0.5, -0.4), (1.2, 0.9, 0.6))
VOXELS = 64
BAKE_DISPATCHES = 20
BAKE_ITERS = 100 # keeps the fixed-point histogram far from uint32 wrap
STILL_SIZE = 192
STILL_DISPATCHES = 40
STILL_ITERS = 200
CPU_POINTS = 120_000
CPU_STEPS = 60 # 0.5**60 is far below float precision: fully converged
# 30 degrees about Y, so the camera is oblique to every world axis.
CAM_YAW = 30.0
CAM_DISTANCE = 3.0
def _mean(points, axis):
return sum(p[axis] for p in points) / len(points)
def _xyz(v):
return {"X": float(v[0]), "Y": float(v[1]), "Z": float(v[2])}
def _camera():
"""Position and orientation, in the engine's convention.
`view_projection` builds the view matrix looking along `position - forward`,
which inverts the usual sense: the points a camera SEES are the ones along
`+forward` from it (that is the documented negative-w quirk in ifs_render).
So framing the centroid at distance d means stepping BACK along forward.
"""
half = math.radians(CAM_YAW) * 0.5
quat = (0.0, math.sin(half), 0.0, math.cos(half)) # (x, y, z, w)
forward = (math.sin(math.radians(CAM_YAW)), 0.0, math.cos(math.radians(CAM_YAW)))
centroid = tuple(_mean(FIXED_POINTS, a) for a in range(3))
position = tuple(centroid[a] - CAM_DISTANCE * forward[a] for a in range(3))
return position, quat
def _document_dict():
position, quat = _camera()
iterators = []
for index, fixed in enumerate(FIXED_POINTS):
iterators.append({
"Transform": {"Item1": "Affine"},
"Vec3Params": {
# Fixed point F of `p := T + S*p` is T/(1-S), so T = F*(1-S).
"Translate": _xyz(tuple(c * (1.0 - S) for c in fixed)),
"Scale": _xyz((S, S, S)),
"Rotate": _xyz((0.0, 0.0, 0.0)),
},
"Opacity": 1.0,
"ColorIndex": index / (len(FIXED_POINTS) - 1.0),
"ColorSpeed": 0.5,
"BaseWeight": 1.0,
"StartWeight": 1.0,
})
return {
"Title": "orientation-probe",
"Iterators": iterators,
"Camera": {
"Position": _xyz(position),
"Orientation": {"X": quat[0], "Y": quat[1], "Z": quat[2], "W": quat[3]},
"FieldOfView": 60.0,
"FocusDistance": CAM_DISTANCE,
},
# White, so the still render's brightness tracks density rather than hue.
"Palette": {"Colors": [{"X": 1.0, "Y": 1.0, "Z": 1.0, "W": 1.0}]},
"Brightness": 4.0,
"Warmup": 20,
"Entropy": 0.01,
"ImageResolution": f"{STILL_SIZE}, {STILL_SIZE}",
}
def _cpu_attractor(seed=20260817):
"""
The same chaos game, on the CPU, in world coordinates.
Independent of the GPU, the readback order, openvdb and the camera — which is
what makes it usable as ground truth for all of them. The RNG does not need to
match the kernel's: an IFS attractor is the invariant SET, so any sequence of
map choices converges to the same distribution.
"""
import numpy as np
rng = np.random.default_rng(seed)
fixed = np.array(FIXED_POINTS, dtype=np.float64)
points = rng.normal(size=(CPU_POINTS, 3))
for _ in range(CPU_STEPS):
choice = rng.integers(0, len(fixed), size=CPU_POINTS)
points = S * points + (1.0 - S) * fixed[choice]
return points
def _read_vdb(path):
"""Every lit voxel of a written `.vdb`, as (world position, density).
Positions come through the grid's OWN transform, which is what Cycles reads —
not a re-derivation from BOUNDS, which would only re-state the assumption under
test.
`indexToWorld(i)` is the voxel's CENTRE, not the corner of its cell. That was
measured, not assumed: a grid with one lit voxel at index (0,0,0) and voxel
size 1 imports into Blender with a bounding box starting at -0.5. So the world
position of a voxel's mass IS its indexToWorld, and no correction belongs here
— `accumulate_voxel` bins a sample into the cell it falls inside, so the grid
transform has to place that cell's CENTRE, which is what `write_vdb` does.
"""
import numpy as np
import openvdb
grid = openvdb.read(path, "density")
array = np.zeros((VOXELS, VOXELS, VOXELS), dtype=np.float32)
grid.copyToArray(array)
lit = np.argwhere(array > 0.0)
weights = array[lit[:, 0], lit[:, 1], lit[:, 2]].astype(np.float64)
world = np.array([grid.transform.indexToWorld(tuple(int(c) for c in idx))
for idx in lit])
return world, weights
def _weighted_stats(positions, weights):
"""Centroid and per-axis mass percentiles, the two statistics compared."""
import numpy as np
total = float(weights.sum())
centroid = (positions * weights[:, None]).sum(axis=0) / total
percentiles = []
for axis in range(3):
order = np.argsort(positions[:, axis])
values = positions[order, axis]
cumulative = np.cumsum(weights[order]) / total
percentiles.append([float(values[np.searchsorted(cumulative, q)])
for q in (0.05, 0.5, 0.95)])
return centroid, percentiles
class Orientation(unittest.TestCase):
doc = None
grid = None
stats = None
@classmethod
def setUpClass(cls):
try:
import numpy # noqa: F401
import openvdb # noqa: F401
import ifs_volume # noqa: F401
except ImportError as ex:
raise unittest.SkipTest(f"needs Blender's gpu/openvdb modules ({ex})")
import json
import ifs_document
import ifs_render
import ifs_volume
cls.tmp = tempfile.mkdtemp()
cls.doc_path = os.path.join(cls.tmp, "orientation.ifsjson")
with open(cls.doc_path, "w", encoding="utf-8") as handle:
json.dump(_document_dict(), handle)
library = ifs_render.load_library(
os.path.join(ifs_render.REPO, "Transforms"),
os.path.join(ifs_render.REPO, "WpfDisplay", "Library", "Includes"))
cls.doc = ifs_document.load(cls.doc_path, library)
cls.points = _cpu_attractor()
# Bake with EXPLICIT bounds: this test is about orientation, so the box is
# held fixed rather than fitted. `test_auto_bounds_brackets_the_fractal`
# covers the fitted path separately.
cls.vdb_path = os.path.join(cls.tmp, "orientation.vdb")
_path, cls.stats = ifs_volume.bake(
cls.doc_path, cls.vdb_path, size=VOXELS, dispatches=BAKE_DISPATCHES,
bounds=BOUNDS, invocation_iters=BAKE_ITERS)
# --------------------------------------------------------------- the guards
def test_the_document_can_detect_a_swap(self):
"""Guards the guard: a symmetric document would make a transpose invisible."""
lo = [min(f[a] for f in FIXED_POINTS) for a in range(3)]
hi = [max(f[a] for f in FIXED_POINTS) for a in range(3)]
spans = [hi[a] - lo[a] for a in range(3)]
self.assertEqual(3, len(set(spans)), f"spans {spans} are not all distinct")
centres = [_mean(FIXED_POINTS, a) for a in range(3)]
self.assertEqual(3, len(set(centres)),
f"centroid {centres} cannot distinguish an axis swap")
for axis in range(3):
self.assertNotAlmostEqual(
lo[axis], -hi[axis], places=6,
msg=f"axis {axis} is symmetric about zero, so a mirror is invisible")
def test_the_cpu_reference_reproduces_the_analytic_attractor(self):
"""The CPU chaos game must land on the tetrahedron, or it is no reference."""
import numpy as np
expected = np.array([_mean(FIXED_POINTS, a) for a in range(3)])
measured = self.points.mean(axis=0)
for axis in range(3):
self.assertAlmostEqual(
expected[axis], measured[axis], places=2,
msg=f"axis {axis}: CPU centroid {measured[axis]:.4f} is not the "
f"analytic {expected[axis]:.4f}")
# The hull is the tetrahedron, so nothing may sit outside it by more than
# float slop.
for axis in range(3):
self.assertGreaterEqual(
self.points[:, axis].min(),
min(f[axis] for f in FIXED_POINTS) - 1e-6)
self.assertLessEqual(
self.points[:, axis].max(),
max(f[axis] for f in FIXED_POINTS) + 1e-6)
# ------------------------------------------------- leg 1: is the bake right
def test_the_bake_lands_where_the_fractal_is(self):
"""
The centroid of the baked density, in world coordinates, against the
analytic centroid. An x/z transpose moves it by 0.10 world units — more
than three voxels — so this cannot pass on a transposed bake.
The tolerance is a QUARTER of a voxel, which is placement accuracy rather
than merely "not transposed". It can afford to be: the residual measured on
a correct bake is 0.01-0.03 voxels, so this still has an order of magnitude
of headroom, and a half-voxel convention slip — the grid transform placing
cell CORNERS where centres belong — shows up as a clean -0.5 on all three
axes rather than passing quietly inside a loose bound.
"""
import numpy as np
positions, weights = _read_vdb(self.vdb_path)
centroid, _pct = _weighted_stats(positions, weights)
expected = np.array([_mean(FIXED_POINTS, a) for a in range(3)])
voxel = np.array([(BOUNDS[1][a] - BOUNDS[0][a]) / VOXELS for a in range(3)])
offset = (centroid - expected) / voxel
print(f"[bake] centroid {np.round(centroid, 4)} expected "
f"{np.round(expected, 4)} offset {np.round(offset, 3)} voxels")
for axis in range(3):
self.assertLess(
abs(offset[axis]), 0.25,
f"axis {axis}: baked centroid {centroid[axis]:.4f} is "
f"{offset[axis]:.2f} voxels from the analytic {expected[axis]:.4f} "
f"— the bake does not sit where the fractal is")
def test_the_bake_has_the_right_shape_on_every_axis(self):
"""
A centroid alone cannot see a swap of two axes with equal means, so compare
the mass distribution per axis: the 5th, 50th and 95th percentiles of the
baked density against the same percentiles of the CPU attractor.
"""
import numpy as np
positions, weights = _read_vdb(self.vdb_path)
_centroid, baked = _weighted_stats(positions, weights)
reference = _weighted_stats(
self.points, np.ones(len(self.points)))[1]
voxel = [(BOUNDS[1][a] - BOUNDS[0][a]) / VOXELS for a in range(3)]
for axis in range(3):
for label, got, want in zip(("p05", "p50", "p95"),
baked[axis], reference[axis]):
drift = (got - want) / voxel[axis]
print(f"[bake] axis {axis} {label}: {got:+.4f} vs {want:+.4f} "
f"({drift:+.2f} voxels)")
self.assertLess(
abs(drift), 2.0,
f"axis {axis} {label}: baked {got:.4f} against the CPU "
f"attractor's {want:.4f}, {drift:.2f} voxels apart")
def test_auto_bounds_brackets_the_fractal(self):
"""
The fitted box is a second consumer of the readback order (`_mass_bounds`
reorders array axes into world xyz), so a wrong order fits a transposed box
while the bake inside it still looks like a fractal.
"""
import numpy as np
import ifs_volume
shader, _library = ifs_volume.build_shader()
ubos = ifs_volume._document_buffers(self.doc)
alias = ifs_volume.ifs_host.build_alias_tables(self.doc.iterators)
n = len(self.doc.iterators)
import gpu
alias_tex = gpu.types.GPUTexture(
(n, n), format='RGBA32F',
data=gpu.types.Buffer('FLOAT', len(alias), alias))
lo, hi, _fit = ifs_volume.auto_bounds(shader, self.doc, ubos, alias_tex)
body_lo = np.percentile(self.points, 1.0, axis=0)
body_hi = np.percentile(self.points, 99.0, axis=0)
span = body_hi - body_lo
print(f"[auto_bounds] lo {np.round(lo, 4)} hi {np.round(hi, 4)} "
f"vs body {np.round(body_lo, 4)}..{np.round(body_hi, 4)}")
for axis in range(3):
self.assertLess(
lo[axis], body_lo[axis],
f"axis {axis}: fitted lo {lo[axis]:.4f} clips the fractal, whose "
f"1st percentile is {body_lo[axis]:.4f}")
self.assertGreater(
hi[axis], body_hi[axis],
f"axis {axis}: fitted hi {hi[axis]:.4f} clips the fractal, whose "
f"99th percentile is {body_hi[axis]:.4f}")
# Generous, but a transposed fit on this document is off by 2-3x.
self.assertLess(
(hi[axis] - lo[axis]), 3.0 * span[axis],
f"axis {axis}: fitted span {(hi[axis] - lo[axis]):.4f} is far wider "
f"than the fractal's {span[axis]:.4f} on that axis")
# ------------------------------- leg 2: is the CPU reference the right one
def _still_density(self):
return still_density(self.doc)
def _projected_density(self):
"""The CPU attractor through the engine's own projection, same binning."""
return projected_density(self.doc.camera.view_proj, self.points)
def test_the_still_render_matches_the_cpu_reference(self):
"""
Earns the reference. If the CPU chaos game and the engine's own projection
put the fractal in the same pixels, the CPU model reproduces the kernel —
which is what makes its verdict on the bake, above, mean anything.
"""
import numpy as np
rendered = self._still_density()
predicted = self._projected_density()
self.assertGreater(float(rendered.sum()), 0.0, "the still render is empty")
self.assertGreater(float(predicted.sum()), 0.0,
"the CPU points all projected off screen")
def centroid(image):
ys, xs = np.nonzero(image)
weights = image[ys, xs]
total = weights.sum()
return np.array([(xs * weights).sum() / total,
(ys * weights).sum() / total])
got, want = centroid(rendered), centroid(predicted)
print(f"[still] centroid {np.round(got, 2)} vs projected CPU "
f"{np.round(want, 2)} px")
for axis, name in enumerate(("x", "y")):
self.assertLess(
abs(got[axis] - want[axis]), 4.0,
f"screen {name}: the render's centroid is {got[axis]:.1f} px and "
f"the projected CPU attractor's is {want[axis]:.1f} px")
# Shape, not just position: a mirror about the centroid leaves it put.
overlap = np.logical_and(rendered > 0, predicted > 0).sum()
union = np.logical_or(rendered > 0, predicted > 0).sum()
iou = overlap / float(union)
print(f"[still] lit-pixel IoU {iou:.3f}")
self.assertGreater(
iou, 0.5,
f"only {iou:.1%} of lit pixels agree — the render and the CPU "
"attractor are not the same shape in the same place")
# ------------------------------------------------------------- the control
def test_a_transposed_bake_would_be_caught(self):
"""
The control. Feeds `write_vdb` an x/z-swapped grid and asserts the centroid
check REJECTS it — otherwise a green suite proves only that the assertions
are loose enough to pass anything.
"""
import numpy as np
import openvdb
import ifs_volume
# Reconstruct an accumulator from the good bake, then swap x and z.
grid = openvdb.read(self.vdb_path, "density")
good = np.zeros((VOXELS, VOXELS, VOXELS), dtype=np.float32)
grid.copyToArray(good)
# `write_vdb` expects [z][y][x] and transposes; handing it the array that
# is ALREADY [x][y][z] reproduces exactly the bug this file guards.
# Reproduce the bug by handing write_vdb data in the OLD [z][y][x] order.
# It no longer transposes — the kernel stores voxel (x,y,z) at texel (z,y,x)
# so the readback already arrives [x][y][z] — which means feeding it a
# transposed array is now what produces a transposed bake. Contiguous, not a
# view: copyFromArray reads memory order and ignores strides.
half = VOXELS // 2
swapped_colour = np.zeros((half, half, half, 3), dtype=np.float32)
path = os.path.join(self.tmp, "transposed.vdb")
ifs_volume.write_vdb(np.ascontiguousarray(good.transpose(2, 1, 0)),
swapped_colour, BOUNDS, path)
bad = openvdb.read(path, "density")
array = np.zeros((VOXELS, VOXELS, VOXELS), dtype=np.float32)
bad.copyToArray(array)
lit = np.argwhere(array > 0.0)
weights = array[lit[:, 0], lit[:, 1], lit[:, 2]].astype(np.float64)
world = np.array([bad.transform.indexToWorld(tuple(int(c) for c in idx))
for idx in lit])
centroid, _pct = _weighted_stats(world, weights)
expected = np.array([_mean(FIXED_POINTS, a) for a in range(3)])
voxel = np.array([(BOUNDS[1][a] - BOUNDS[0][a]) / VOXELS for a in range(3)])
offset = np.abs(centroid - expected) / voxel
print(f"[control] transposed centroid {np.round(centroid, 4)} "
f"offset {np.round(offset, 2)} voxels")
self.assertGreater(
offset.max(), 1.5,
"a transposed bake landed within tolerance of the correct centroid — "
"the orientation assertions above cannot detect a transpose")
def still_density(doc, size=STILL_SIZE, dispatches=STILL_DISPATCHES,
iters=STILL_ITERS):
"""
The still path's dispatch loop, for the raw density it tonemaps away.
`render()` returns an image datablock, and the tonemap plus the 8-bit buffer
destroy exactly the faint tail this comparison needs. The SHADER is the same one
`render()` builds, so nothing about the projection is reproduced here — only the
loop around it.
Module level rather than a method so `test_points_emit` can measure the emitted
cloud against the same still render instead of re-deriving one.
"""
import gpu
import numpy as np
import ifs_host
import ifs_render
camera = doc.camera
camera.view_proj = ifs_render.view_projection(
camera.position[:3], camera.forward[:3], camera.up,
fov_degrees=camera.field_of_view)
kernel, drain, _library = ifs_render._build_shaders(
size, size, ifs_render.DEFAULT_HIST_SCALE)
ubos = {
'settings': gpu.types.GPUUniformBuf(
ifs_host.pack_settings(camera, doc.settings)),
'iterators_buf': gpu.types.GPUUniformBuf(
ifs_host.pack_iterators(doc.iterators)),
'palette_buf': gpu.types.GPUUniformBuf(
ifs_host.pack_palette(doc.palette(doc.settings.gradient_resolution))),
'real_params_buf': gpu.types.GPUUniformBuf(
ifs_host.pack_real_params(ifs_host.flatten_real_params(doc.iterators))),
'vec3_params_buf': gpu.types.GPUUniformBuf(
ifs_host.pack_vec3_params(ifs_host.flatten_vec3_params(doc.iterators))),
}
alias = ifs_host.build_alias_tables(doc.iterators)
n = len(doc.iterators)
alias_tex = gpu.types.GPUTexture(
(n, n), format='RGBA32F',
data=gpu.types.Buffer('FLOAT', len(alias), alias))
histogram = gpu.types.GPUTexture((size * 4, size), format='R32UI')
histogram.clear(format='UINT', value=(0, 0, 0, 0))
accum = gpu.types.GPUTexture((size, size), format='RGBA32F')
accum.clear(format='FLOAT', value=(0.0, 0.0, 0.0, 0.0))
state_rows = (ifs_render.INVOCATIONS + ifs_render.STATE_W - 1) // ifs_render.STATE_W
state = gpu.types.GPUTexture(
(ifs_render.STATE_W * 2, state_rows), format='RGBA32F')
state.clear(format='FLOAT', value=(0.0, 0.0, 0.0, 0.0))
for index in range(dispatches):
frame = gpu.types.GPUUniformBuf(ifs_host.pack_frame(
seed=0x9E3779B9, width=size, height=size,
dispatch_cnt=index, reset_points_state=1 if index == 0 else 0,
invocation_iters=iters))
kernel.bind()
kernel.image('histogram_img', histogram)
kernel.image('state_img', state)
kernel.uniform_sampler('alias_tables_tex', alias_tex)
kernel.uniform_block('frame', frame)
for name, ubo in ubos.items():
kernel.uniform_block(name, ubo)
gpu.compute.dispatch(kernel, ifs_render.WORKGROUPS, 1, 1)
drain.bind()
drain.image('histogram_img', histogram)
drain.image('accum_img', accum)
drain.uniform_int('out_size', (size, size))
gpu.compute.dispatch(drain, (size + 7) // 8, (size + 7) // 8, 1)
pixels = np.array(accum.read(), dtype=np.float32).reshape(size, size, 4)
return pixels[:, :, 3]
def projected_density(view_proj, points, size=STILL_SIZE):
"""A point set through the engine's own projection, binned like the render."""
import numpy as np
import ifs_render
out = np.zeros((size, size), dtype=np.float64)
for point in points:
result = ifs_render.project_point(
view_proj, (float(point[0]), float(point[1]), float(point[2])),
size, size)
if result is None:
continue
x, y, _w = result
out[int(round(y)), int(round(x))] += 1.0
return out
if __name__ == "__main__":
result = unittest.main(argv=sys.argv[:1], exit=False, verbosity=2).result
sys.exit(0 if result.wasSuccessful() else 1)