| … | … | @@ -366,6 +366,90 @@ VOLUME_LAYOUT: tuple[Resource, ...] = ( |
| 366 | 366 | ) |
| 367 | 367 | |
| 368 | 368 | |
| 369 | # ------------------------------------------------------------------- points mode | |
| 370 | # | |
| 371 | # Emit the chaos game's samples as POINTS instead of binning them, and hand them to | |
| 372 | # Cycles as a PointCloud. Where volume mode trades resolution against memory - the | |
| 373 | # grid is 4N x N x N, so 512 is 2.1 GB and does not fit - points sit at float | |
| 374 | # precision and have no resolution ceiling at all. The cost moves from "how finely | |
| 375 | # can the grid resolve it" to "how many samples can be stored", which is a straight | |
| 376 | # linear trade instead of a cubic one. | |
| 377 | # | |
| 378 | # The obstacle is that there is no buffer to append to: Blender's GPU module has no | |
| 379 | # SSBOs (module docstring, reason 1). An image is the substitute - one texel per | |
| 380 | # point in an RGBA32F image - and the append is an atomic counter in a 1x1 R32UI, | |
| 381 | # which is the viewport's auto-exposure max pass turned from a reduction into an | |
| 382 | # allocator. | |
| 383 | # | |
| 384 | # The plot block is replaced wholesale for the same reason volume mode replaces it: | |
| 385 | # project() can `continue` on out-of-frame, so anything that merely ignored the | |
| 386 | # projection would still drop every sample outside the current view and emit a | |
| 387 | # view-dependent cloud. | |
| 388 | ||
| 389 | # The kernel draws the chaos game's map choice from `gl_WorkGroupID.x`, so all 64 | |
| 390 | # invocations in a workgroup follow the SAME sequence of maps. That is deliberate and | |
| 391 | # it is a good trade for a histogram: every thread in the workgroup takes the same | |
| 392 | # branch through the transform switch, which is a large SIMD win, and the duplicate | |
| 393 | # trajectories it produces merely add density to the same bins. Nothing is lost when | |
| 394 | # the output is a histogram accumulated over hundreds of dispatches. | |
| 395 | # | |
| 396 | # For a point cloud it is quietly ruinous, because duplicates become actual duplicate | |
| 397 | # points and one dispatch is the whole sample. MEASURED on the orientation document at | |
| 398 | # 200k points: 27% of emitted points were exact duplicates, the four maps came out | |
| 399 | # [50216, 46516, 50951, 53021] against 50176 expected - a 7.3% deficit on one of them, | |
| 400 | # which showed up as the attractor's centroid sitting 0.014 low on the axis that map | |
| 401 | # is responsible for. The draw pool is only `workgroups x iterations` values wide, | |
| 402 | # so a single dispatch never averages it out. | |
| 403 | # | |
| 404 | # Points mode therefore draws per INVOCATION. It costs branch divergence inside each | |
| 405 | # workgroup - the throughput trade the engine deliberately made in the other | |
| 406 | # direction - and buys 16384 independent chains instead of 256. | |
| 407 | _WORKGROUP_DRAW = re.compile(r"f_hash4\(seed, gl_WorkGroupID\.x,") | |
| 408 | ||
| 409 | _PLOT_POINTS = """ vec4 color = vec4(getPaletteColor(p.color_index), selected_iterator.opacity); | |
| 410 | if (color.w == 0.0) | |
| 411 | continue; | |
| 412 | emit_point(p.pos.xyz, color);""" | |
| 413 | ||
| 414 | _EMIT_HELPER = """ | |
| 415 | // An append-only list built out of an image, because there are no SSBOs to append | |
| 416 | // to. `imageAtomicAdd` returns the PREVIOUS value, so every invocation claims a | |
| 417 | // unique slot with no lock and no ordering assumption. | |
| 418 | // | |
| 419 | // The counter deliberately keeps counting PAST the end of the image and the store | |
| 420 | // is guarded instead. That costs one comparison and buys the host the difference | |
| 421 | // between "this cloud is the whole run" and "this cloud is the first million of | |
| 422 | // twelve" - a truncated point cloud is still a plausible fractal, so silence there | |
| 423 | // would be the same class of quiet wrongness as a transposed bake. | |
| 424 | void emit_point(vec3 pos, vec4 color) | |
| 425 | { | |
| 426 | uint slot = imageAtomicAdd(counter_img, ivec2(0, 0), 1u); | |
| 427 | if (slot >= uint(point_capacity)) | |
| 428 | return; | |
| 429 | ivec2 texel = ivec2(int(slot % uint(points_width)), int(slot / uint(points_width))); | |
| 430 | imageStore(points_img, texel, vec4(pos, 1.0)); | |
| 431 | imageStore(colors_img, texel, color); | |
| 432 | } | |
| 433 | """ | |
| 434 | ||
| 435 | POINTS_LAYOUT: tuple[Resource, ...] = ( | |
| 436 | Resource("image", 0, "points_img", "FLOAT_2D", | |
| 437 | "RGBA32F, points_width x ceil(capacity/points_width). xyz + 1 per point."), | |
| 438 | Resource("image", 1, "colors_img", "FLOAT_2D", | |
| 439 | "RGBA32F, same shape. Palette rgb + the iterator's opacity."), | |
| 440 | Resource("image", 2, "counter_img", "UINT_2D", | |
| 441 | "R32UI 1x1. imageAtomicAdd allocator; counts past capacity on purpose."), | |
| 442 | Resource("image", 3, "state_img", "FLOAT_2D", "as the 2D path"), | |
| 443 | Resource("sampler", 0, "alias_tables_tex", "FLOAT_2D", ""), | |
| 444 | Resource("uniform_buf", 0, "frame_ubo", "FrameParams", "", instance="frame"), | |
| 445 | Resource("uniform_buf", 1, "settings_ubo", "SettingsParams", "", instance="settings"), | |
| 446 | Resource("uniform_buf", 2, "iterators_ubo", "IteratorsParams", "", instance="iterators_buf"), | |
| 447 | Resource("uniform_buf", 3, "palette_ubo", "PaletteParams", "", instance="palette_buf"), | |
| 448 | Resource("uniform_buf", 4, "real_params_ubo", "RealParams", "", instance="real_params_buf"), | |
| 449 | Resource("uniform_buf", 5, "vec3_params_ubo", "Vec3Params", "", instance="vec3_params_buf"), | |
| 450 | ) | |
| 451 | ||
| 452 | ||
| 369 | 453 | @dataclass |
| 370 | 454 | class RewriteReport: |
| 371 | 455 | """What the rewrite did, so a silent no-op is impossible to miss.""" |
| … | … | @@ -381,7 +465,8 @@ class RewriteReport: |
| 381 | 465 | def rewrite(source: str, hist_scale: float = DEFAULT_HIST_SCALE, |
| 382 | 466 | state_width: int = 256, |
| 383 | 467 | volume: bool = False, |
| 384 | depth: bool = False) -> tuple[str, str, RewriteReport]: | |
| 468 | depth: bool = False, | |
| 469 | points: bool = False) -> tuple[str, str, RewriteReport]: | |
| 385 | 470 | """ |
| 386 | 471 | Rewrites an engine kernel onto the Blender resource model. |
| 387 | 472 | |
| … | … | @@ -392,10 +477,22 @@ def rewrite(source: str, hist_scale: float = DEFAULT_HIST_SCALE, |
| 392 | 477 | `depth` adds the nearest-sample depth channel the viewport overlay needs to be |
| 393 | 478 | occluded by scene geometry. It is meaningless with `volume`, which replaces the |
| 394 | 479 | 2D projection wholesale and so has no screen pixel to be in front of. |
| 480 | ||
| 481 | `volume` and `points` are the two ways of turning the chaos game into real scene | |
| 482 | data - a density grid or a list of samples. They are alternatives, not layers: | |
| 483 | each replaces the plot block with its own write, and each binds image slot 0 to | |
| 484 | something the other would read as a different type. | |
| 395 | 485 | """ |
| 396 | 486 | if volume and depth: |
| 397 | 487 | raise ValueError("depth is a 2D-projection concept; it has no meaning in " |
| 398 | 488 | "volume mode, which never projects to a screen pixel") |
| 489 | chosen = [name for name, on in | |
| 490 | (("volume", volume), ("points", points), ("depth", depth)) if on] | |
| 491 | if len(chosen) > 1: | |
| 492 | raise ValueError( | |
| 493 | f"{' and '.join(chosen)} cannot be combined: each replaces the plot " | |
| 494 | "block with a different write and binds the image slots under " | |
| 495 | "different types") | |
| 399 | 496 | applied: list[str] = [] |
| 400 | 497 | missing: list[str] = [] |
| 401 | 498 | |
| … | … | @@ -427,6 +524,14 @@ def rewrite(source: str, hist_scale: float = DEFAULT_HIST_SCALE, |
| 427 | 524 | "replace the projection block with a voxel write", out) |
| 428 | 525 | out = sub(_ACCUMULATE, lambda _m: _VOXEL_HELPER, |
| 429 | 526 | "swap accumulate_hit for accumulate_voxel", out) |
| 527 | elif points: | |
| 528 | out = sub(_PLOT_BLOCK, lambda _m: _PLOT_POINTS, | |
| 529 | "replace the projection block with a point emit", out) | |
| 530 | out = sub(_ACCUMULATE, lambda _m: _EMIT_HELPER, | |
| 531 | "swap accumulate_hit for emit_point", out) | |
| 532 | out = sub(_WORKGROUP_DRAW, "f_hash4(seed, gl_GlobalInvocationID.x,", | |
| 533 | "draw the chaos game per invocation, not per workgroup", | |
| 534 | out, expected=3) | |
| 430 | 535 | elif depth: |
| 431 | 536 | out = sub(_ACCUMULATE, lambda _m: _ACCUMULATE_DEPTH, |
| 432 | 537 | "rewrite accumulate_hit with depth", out) |
| … | … | @@ -0,0 +1,355 @@ |
| 1 | """ | |
| 2 | Emits the fractal as a PointCloud, so Cycles and EEVEE both render it as scene data. | |
| 3 | ||
| 4 | This is the third way of getting a flame out of the chaos game, and the first with | |
| 5 | no resolution ceiling. | |
| 6 | ||
| 7 | * The **live viewport** is a screen-space overlay: no depth of its own, invisible | |
| 8 | to F12. It is a preview. | |
| 9 | * The **volume bake** is real scene data, but a voxel grid is 4N x N x N uint32 on | |
| 10 | the GPU - 268 MB at N=256 and 2.1 GB at N=512, which does not fit. Detail is | |
| 11 | bought cubically, so a flame's thin filaments end up as mush. | |
| 12 | * **Points** store what the chaos game actually produces: samples, at float | |
| 13 | precision. Nothing is quantised, so there is no grid to be finer than. The cost | |
| 14 | is linear in the number of samples rather than cubic in the resolution, and the | |
| 15 | same representation renders in Cycles properly and in EEVEE as an in-scene | |
| 16 | approximation - one representation, both engines. | |
| 17 | ||
| 18 | ## How the points get off the GPU | |
| 19 | ||
| 20 | There is nowhere to append them to: Blender's GPU module has no SSBOs. The | |
| 21 | substitute is an RGBA32F image used as a flat array, one texel per point, with an | |
| 22 | atomic counter in a 1x1 R32UI handing out slots - see `ifs_layout`'s points mode. | |
| 23 | The host allocates `points_width x ceil(capacity/points_width)`, dispatches until | |
| 24 | the counter says the array is full, and reads the two images back. | |
| 25 | ||
| 26 | Because slots are claimed in atomic order, a full buffer holds the FIRST `capacity` | |
| 27 | samples of the run - and WHICH samples those are is a sampling-quality question, not | |
| 28 | just a performance one. See `iters_for`: the dispatch is sized so that one pass over | |
| 29 | all 16384 chains fills the buffer exactly once, which is what makes those samples | |
| 30 | `capacity / 16384` iterations from EVERY chain rather than a long correlated run | |
| 31 | from whichever few thousand the GPU scheduled first. | |
| 32 | ||
| 33 | ## Radius | |
| 34 | ||
| 35 | Cycles draws each point as a sphere, so the cloud only reads as a continuous object | |
| 36 | if the radius is near the local spacing between samples. That spacing is MEASURED | |
| 37 | rather than assumed - see `estimate_spacing`, which also falls out as an estimate of | |
| 38 | the attractor's fractal dimension. It is a starting point, not a truth: how big a | |
| 39 | point should look also depends on how close the camera is, so the operator exposes a | |
| 40 | multiplier. | |
| 41 | """ | |
| 42 | ||
| 43 | from __future__ import annotations | |
| 44 | ||
| 45 | import os | |
| 46 | ||
| 47 | import gpu | |
| 48 | import numpy as np | |
| 49 | ||
| 50 | import ifs_codegen | |
| 51 | import ifs_document | |
| 52 | import ifs_host | |
| 53 | import ifs_layout | |
| 54 | import ifs_render | |
| 55 | ||
| 56 | # One texel per point. 4096 is comfortably inside every GL_MAX_TEXTURE_SIZE worth | |
| 57 | # worrying about (16384 is the common floor) and keeps the row count small enough | |
| 58 | # that a 16M-point buffer is 4096 rows. | |
| 59 | POINTS_WIDTH = 4096 | |
| 60 | ||
| 61 | # NOT "color". That name is poisoned for points and meshes in Blender - a per-point | |
| 62 | # attribute called `color` renders black - while being the EXPECTED name for a | |
| 63 | # volume grid. Both halves of that cost a wrong conclusion during the volume work. | |
| 64 | COLOR_ATTRIBUTE = "ifs_color" | |
| 65 | ||
| 66 | ||
| 67 | def build_shader(): | |
| 68 | """Compiles the points variant of the kernel.""" | |
| 69 | library = ifs_codegen.load_library( | |
| 70 | os.path.join(ifs_render.REPO, "Transforms"), | |
| 71 | os.path.join(ifs_render.REPO, "WpfDisplay", "Library", "Includes")) | |
| 72 | template_path = os.path.join( | |
| 73 | ifs_render.REPO, "IFSEngine", "Rendering", "Shaders", "ifs_kernel.comp.shader") | |
| 74 | with open(template_path, "r", encoding="utf-8-sig") as handle: | |
| 75 | kernel = ifs_codegen.generate_kernel(handle.read(), library) | |
| 76 | kernel, _ = ifs_codegen.apply_portability_fixes(kernel) | |
| 77 | body, typedefs, report = ifs_layout.rewrite( | |
| 78 | kernel, ifs_layout.DEFAULT_HIST_SCALE, ifs_render.STATE_W, points=True) | |
| 79 | if not report.ok: | |
| 80 | raise RuntimeError("points rewrite incomplete: " + "; ".join(report.missing)) | |
| 81 | ||
| 82 | info = gpu.types.GPUShaderCreateInfo() | |
| 83 | info.typedef_source(typedefs) | |
| 84 | info.local_group_size(ifs_render.LOCAL_SIZE) | |
| 85 | rw = {'NO_RESTRICT', 'READ', 'WRITE'} | |
| 86 | info.image(0, 'RGBA32F', 'FLOAT_2D', 'points_img', qualifiers=rw) | |
| 87 | info.image(1, 'RGBA32F', 'FLOAT_2D', 'colors_img', qualifiers=rw) | |
| 88 | info.image(2, 'R32UI', 'UINT_2D', 'counter_img', qualifiers=rw) | |
| 89 | info.image(3, 'RGBA32F', 'FLOAT_2D', 'state_img', qualifiers=rw) | |
| 90 | info.sampler(0, 'FLOAT_2D', 'alias_tables_tex') | |
| 91 | info.push_constant('INT', 'point_capacity') | |
| 92 | info.push_constant('INT', 'points_width') | |
| 93 | info.uniform_buf(0, 'FrameParams', 'frame') | |
| 94 | info.uniform_buf(1, 'SettingsParams', 'settings') | |
| 95 | info.uniform_buf(2, 'IteratorsParams', 'iterators_buf') | |
| 96 | info.uniform_buf(3, 'PaletteParams', 'palette_buf') | |
| 97 | info.uniform_buf(4, 'RealParams', 'real_params_buf') | |
| 98 | info.uniform_buf(5, 'Vec3Params', 'vec3_params_buf') | |
| 99 | info.compute_source(body) | |
| 100 | return gpu.shader.create_from_info(info), library | |
| 101 | ||
| 102 | ||
| 103 | def _document_buffers(doc): | |
| 104 | """The five per-document UBOs. | |
| 105 | ||
| 106 | Deliberately a copy of `ifs_volume._document_buffers` rather than an import: | |
| 107 | it is mechanical std140 packing, and coupling the two bake modes through a | |
| 108 | private name would mean a change for one silently reshapes the other. | |
| 109 | """ | |
| 110 | return { | |
| 111 | 'settings': gpu.types.GPUUniformBuf( | |
| 112 | ifs_host.pack_settings(doc.camera, doc.settings)), | |
| 113 | 'iterators_buf': gpu.types.GPUUniformBuf( | |
| 114 | ifs_host.pack_iterators(doc.iterators)), | |
| 115 | 'palette_buf': gpu.types.GPUUniformBuf( | |
| 116 | ifs_host.pack_palette(doc.palette(doc.settings.gradient_resolution))), | |
| 117 | 'real_params_buf': gpu.types.GPUUniformBuf( | |
| 118 | ifs_host.pack_real_params(ifs_host.flatten_real_params(doc.iterators))), | |
| 119 | 'vec3_params_buf': gpu.types.GPUUniformBuf( | |
| 120 | ifs_host.pack_vec3_params(ifs_host.flatten_vec3_params(doc.iterators))), | |
| 121 | } | |
| 122 | ||
| 123 | ||
| 124 | def iters_for(capacity, warmup): | |
| 125 | """ | |
| 126 | How many iterations each invocation should run to fill `capacity` ONCE. | |
| 127 | ||
| 128 | This is a sampling-quality decision, not a performance one. A dispatch at the | |
| 129 | engine's default 500 iterations emits ~8M samples from 16384 chains; against a | |
| 130 | 1M buffer, seven eighths of them hit a full array and are dropped. That would be | |
| 131 | merely wasteful if the survivors were a fair sample, but they are not: slots are | |
| 132 | claimed in atomic order, so the buffer fills from whichever invocations the GPU | |
| 133 | happens to schedule FIRST. Only a few thousand chains are resident at once, so a | |
| 134 | 500-iteration dispatch fills the array from a fraction of the chains, each | |
| 135 | contributing a long correlated run. | |
| 136 | ||
| 137 | Sizing the dispatch to the buffer instead means every one of the 16384 chains | |
| 138 | contributes its share - `capacity / 16384` samples each - which is the widest | |
| 139 | front the kernel can offer. The warmup is added on top because those iterations | |
| 140 | are discarded before anything is plotted. | |
| 141 | """ | |
| 142 | per_chain = (capacity + ifs_render.INVOCATIONS - 1) // ifs_render.INVOCATIONS | |
| 143 | return int(max(warmup + 8, warmup + per_chain)) | |
| 144 | ||
| 145 | ||
| 146 | def emit(shader, doc, ubos, alias_tex, capacity, max_dispatches=400, | |
| 147 | invocation_iters=None, poll_every=1): | |
| 148 | """ | |
| 149 | Runs the chaos game until `capacity` points are stored or the budget runs out. | |
| 150 | ||
| 151 | Returns `(positions, colors, stats)` with positions (N, 3) and colors (N, 4), | |
| 152 | both float32 and both in the order the GPU claimed slots. | |
| 153 | ||
| 154 | The counter is polled rather than trusted at the end, because once the array is | |
| 155 | full every further dispatch is pure waste and the counter is a uint32 that a | |
| 156 | long enough run would wrap. With `invocation_iters` derived from the capacity a | |
| 157 | single dispatch very nearly fills the buffer, so polling every dispatch costs | |
| 158 | one readback rather than serialising anything. | |
| 159 | """ | |
| 160 | if invocation_iters is None: | |
| 161 | invocation_iters = iters_for(capacity, doc.settings.warmup) | |
| 162 | rows = max(1, (capacity + POINTS_WIDTH - 1) // POINTS_WIDTH) | |
| 163 | allocated = rows * POINTS_WIDTH | |
| 164 | ||
| 165 | points_tex = gpu.types.GPUTexture((POINTS_WIDTH, rows), format='RGBA32F') | |
| 166 | points_tex.clear(format='FLOAT', value=(0.0, 0.0, 0.0, 0.0)) | |
| 167 | colors_tex = gpu.types.GPUTexture((POINTS_WIDTH, rows), format='RGBA32F') | |
| 168 | colors_tex.clear(format='FLOAT', value=(0.0, 0.0, 0.0, 0.0)) | |
| 169 | counter = gpu.types.GPUTexture((1, 1), format='R32UI') | |
| 170 | counter.clear(format='UINT', value=(0, 0, 0, 0)) | |
| 171 | ||
| 172 | state_rows = (ifs_render.INVOCATIONS + ifs_render.STATE_W - 1) // ifs_render.STATE_W | |
| 173 | state = gpu.types.GPUTexture( | |
| 174 | (ifs_render.STATE_W * 2, state_rows), format='RGBA32F') | |
| 175 | state.clear(format='FLOAT', value=(0.0, 0.0, 0.0, 0.0)) | |
| 176 | ||
| 177 | def claimed(): | |
| 178 | return int(np.array(counter.read(), dtype=np.uint32).ravel()[0]) | |
| 179 | ||
| 180 | dispatched = 0 | |
| 181 | requested = 0 | |
| 182 | for index in range(max_dispatches): | |
| 183 | frame = gpu.types.GPUUniformBuf(ifs_host.pack_frame( | |
| 184 | seed=0x9E3779B9, width=POINTS_WIDTH, height=rows, dispatch_cnt=index, | |
| 185 | reset_points_state=1 if index == 0 else 0, | |
| 186 | invocation_iters=invocation_iters)) | |
| 187 | shader.bind() | |
| 188 | shader.image('points_img', points_tex) | |
| 189 | shader.image('colors_img', colors_tex) | |
| 190 | shader.image('counter_img', counter) | |
| 191 | shader.image('state_img', state) | |
| 192 | shader.uniform_sampler('alias_tables_tex', alias_tex) | |
| 193 | shader.uniform_int('point_capacity', (allocated,)) | |
| 194 | shader.uniform_int('points_width', (POINTS_WIDTH,)) | |
| 195 | shader.uniform_block('frame', frame) | |
| 196 | for name, ubo in ubos.items(): | |
| 197 | shader.uniform_block(name, ubo) | |
| 198 | gpu.compute.dispatch(shader, ifs_render.WORKGROUPS, 1, 1) | |
| 199 | dispatched = index + 1 | |
| 200 | ||
| 201 | if dispatched % poll_every == 0: | |
| 202 | requested = claimed() | |
| 203 | if requested >= allocated: | |
| 204 | break | |
| 205 | ||
| 206 | requested = max(requested, claimed()) | |
| 207 | stored = min(requested, allocated) | |
| 208 | if stored == 0: | |
| 209 | raise RuntimeError( | |
| 210 | "no points were emitted - every iterator may have zero opacity, or the " | |
| 211 | "document may produce nothing") | |
| 212 | ||
| 213 | positions = np.array(points_tex.read(), dtype=np.float32).reshape(-1, 4) | |
| 214 | colors = np.array(colors_tex.read(), dtype=np.float32).reshape(-1, 4) | |
| 215 | return positions[:stored, :3].copy(), colors[:stored].copy(), { | |
| 216 | "points": stored, | |
| 217 | "requested": requested, | |
| 218 | "capacity": allocated, | |
| 219 | "dispatches": dispatched, | |
| 220 | "filled": bool(requested >= allocated), | |
| 221 | } | |
| 222 | ||
| 223 | ||
| 224 | def estimate_spacing(points, sample=2048, seed=0): | |
| 225 | """ | |
| 226 | The typical distance between neighbouring samples, MEASURED, plus the fractal | |
| 227 | dimension it implies. | |
| 228 | ||
| 229 | A radius has to come from somewhere, and "somewhere" is usually a number typed | |
| 230 | once and tuned by eye - which is exactly what went wrong with the volume path's | |
| 231 | emission scale, where the right value differed by 200x between documents. | |
| 232 | ||
| 233 | The measurement: mean nearest-neighbour distance over a random subsample. That | |
| 234 | distance shrinks with sample count as `d ~ n**(-1/D)` for a set of dimension D, | |
| 235 | so measuring it at TWO sizes gives D without assuming one, and extrapolating to | |
| 236 | the full count gives the spacing of the cloud actually being built. A flame is | |
| 237 | typically somewhere between a surface and a solid, so D lands near 2 - and it is | |
| 238 | clamped to [1, 3] because a degenerate document (every sample in one spot) would | |
| 239 | otherwise produce a nonsense exponent and a nonsense radius. | |
| 240 | """ | |
| 241 | rng = np.random.default_rng(seed) | |
| 242 | total = len(points) | |
| 243 | big = min(sample, total) | |
| 244 | small = max(4, big // 4) | |
| 245 | if total < 16: | |
| 246 | span = float(np.ptp(points, axis=0).max()) if total else 1.0 | |
| 247 | return max(span, 1e-6) / max(total, 1), 2.0 | |
| 248 | ||
| 249 | def nearest(count): | |
| 250 | index = rng.choice(total, size=count, replace=False) | |
| 251 | chosen = points[index].astype(np.float64) | |
| 252 | square = (chosen ** 2).sum(axis=1) | |
| 253 | # |a-b|^2 via the dot-product identity: an (n, n) matrix instead of the | |
| 254 | # (n, n, 3) a broadcast subtraction would build. | |
| 255 | distance = square[:, None] + square[None, :] - 2.0 * (chosen @ chosen.T) | |
| 256 | np.fill_diagonal(distance, np.inf) | |
| 257 | return float(np.sqrt(np.maximum(distance.min(axis=1), 0.0)).mean()) | |
| 258 | ||
| 259 | d_big = nearest(big) | |
| 260 | d_small = nearest(small) | |
| 261 | ratio = d_small / d_big if d_big > 0 else 1.0 | |
| 262 | if ratio > 1.0 and big > small: | |
| 263 | dimension = float(np.log(big / small) / np.log(ratio)) | |
| 264 | else: | |
| 265 | dimension = 2.0 | |
| 266 | dimension = float(min(3.0, max(1.0, dimension))) | |
| 267 | spacing = d_big * (total / float(big)) ** (-1.0 / dimension) | |
| 268 | return max(spacing, 1e-9), dimension | |
| 269 | ||
| 270 | ||
| 271 | def emit_document(doc, capacity=1_000_000, max_dispatches=400, | |
| 272 | invocation_iters=None, shader=None): | |
| 273 | """ | |
| 274 | Emits a document that is ALREADY loaded, which is what makes the cloud follow | |
| 275 | what is on screen. | |
| 276 | ||
| 277 | Same contract as `ifs_volume.bake_document`, and for the same reason: the node | |
| 278 | graph and the panel sliders write into the live session's document, so anything | |
| 279 | that reloaded from disk would silently discard every edit made inside Blender. | |
| 280 | ||
| 281 | The document's camera is unused - points mode replaces the projection block | |
| 282 | outright, so nothing here depends on where the camera is looking. | |
| 283 | """ | |
| 284 | gpu.init() | |
| 285 | if shader is None: | |
| 286 | shader, _library = build_shader() | |
| 287 | ubos = _document_buffers(doc) | |
| 288 | ||
| 289 | alias = ifs_host.build_alias_tables(doc.iterators) | |
| 290 | n = len(doc.iterators) | |
| 291 | alias_tex = gpu.types.GPUTexture( | |
| 292 | (n, n), format='RGBA32F', | |
| 293 | data=gpu.types.Buffer('FLOAT', len(alias), alias)) | |
| 294 | ||
| 295 | positions, colors, stats = emit( | |
| 296 | shader, doc, ubos, alias_tex, capacity, max_dispatches, invocation_iters) | |
| 297 | ||
| 298 | spacing, dimension = estimate_spacing(positions) | |
| 299 | lo = positions.min(axis=0) | |
| 300 | hi = positions.max(axis=0) | |
| 301 | stats.update({ | |
| 302 | "spacing": spacing, | |
| 303 | # Spheres of half the spacing just touch, which is the smallest radius that | |
| 304 | # still reads as a continuous surface rather than dust. | |
| 305 | "radius": 0.5 * spacing, | |
| 306 | "dimension": dimension, | |
| 307 | "bounds": (tuple(float(v) for v in lo), tuple(float(v) for v in hi)), | |
| 308 | }) | |
| 309 | return positions, colors, stats | |
| 310 | ||
| 311 | ||
| 312 | def emit_from_file(document_path, capacity=1_000_000, max_dispatches=400, | |
| 313 | invocation_iters=None): | |
| 314 | """Loads a document from disk and emits it. Returns (positions, colors, stats).""" | |
| 315 | gpu.init() | |
| 316 | shader, library = build_shader() | |
| 317 | doc = ifs_document.load(document_path, library) | |
| 318 | return emit_document(doc, capacity, max_dispatches, invocation_iters, | |
| 319 | shader=shader) | |
| 320 | ||
| 321 | ||
| 322 | def to_pointcloud(positions, colors, radius, name="IFS Points"): | |
| 323 | """ | |
| 324 | Builds a Blender PointCloud from the emitted arrays. | |
| 325 | ||
| 326 | Everything goes through the ATTRIBUTE api rather than the `points` collection: | |
| 327 | measured at 1M points, `attributes['position'].data.foreach_set('vector', ...)` | |
| 328 | takes 1 ms where `points.foreach_set('co', ...)` takes 46 and radius through the | |
| 329 | collection takes 90. Same data, same result, forty times cheaper - the | |
| 330 | collection path wraps every element in RNA on the way through. | |
| 331 | """ | |
| 332 | import bpy | |
| 333 | ||
| 334 | count = len(positions) | |
| 335 | cloud = bpy.data.pointclouds.new(name) | |
| 336 | cloud.resize(count) # creates the `position` attribute | |
| 337 | ||
| 338 | cloud.attributes['position'].data.foreach_set( | |
| 339 | 'vector', np.ascontiguousarray(positions, dtype=np.float32).ravel()) | |
| 340 | ||
| 341 | colour = cloud.attributes.get(COLOR_ATTRIBUTE) | |
| 342 | if colour is None: | |
| 343 | colour = cloud.attributes.new( | |
| 344 | name=COLOR_ATTRIBUTE, type='FLOAT_COLOR', domain='POINT') | |
| 345 | rgba = np.ones((count, 4), dtype=np.float32) | |
| 346 | rgba[:, :colors.shape[1]] = colors[:, :4] | |
| 347 | rgba[:, 3] = 1.0 # opacity rides in the colour's w on the GPU side | |
| 348 | colour.data.foreach_set('color', rgba.ravel()) | |
| 349 | ||
| 350 | radii = cloud.attributes.get('radius') | |
| 351 | if radii is None: | |
| 352 | radii = cloud.attributes.new(name='radius', type='FLOAT', domain='POINT') | |
| 353 | radii.data.foreach_set('value', np.full(count, radius, dtype=np.float32)) | |
| 354 | ||
| 355 | return cloud |
| … | … | @@ -373,95 +373,12 @@ class Orientation(unittest.TestCase): |
| 373 | 373 | # ------------------------------- leg 2: is the CPU reference the right one |
| 374 | 374 | |
| 375 | 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] | |
| 376 | return still_density(self.doc) | |
| 449 | 377 | |
| 450 | 378 | def _projected_density(self): |
| 451 | 379 | """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 | |
| 380 | return projected_density(self.doc.camera.view_proj, self.points) | |
| 381 | ||
| 465 | 382 | |
| 466 | 383 | def test_the_still_render_matches_the_cpu_reference(self): |
| 467 | 384 | """ |
| … | … | @@ -545,6 +462,98 @@ class Orientation(unittest.TestCase): |
| 545 | 462 | "the orientation assertions above cannot detect a transpose") |
| 546 | 463 | |
| 547 | 464 | |
| 465 | def still_density(doc, size=STILL_SIZE, dispatches=STILL_DISPATCHES, | |
| 466 | iters=STILL_ITERS): | |
| 467 | """ | |
| 468 | The still path's dispatch loop, for the raw density it tonemaps away. | |
| 469 | ||
| 470 | `render()` returns an image datablock, and the tonemap plus the 8-bit buffer | |
| 471 | destroy exactly the faint tail this comparison needs. The SHADER is the same one | |
| 472 | `render()` builds, so nothing about the projection is reproduced here — only the | |
| 473 | loop around it. | |
| 474 | ||
| 475 | Module level rather than a method so `test_points_emit` can measure the emitted | |
| 476 | cloud against the same still render instead of re-deriving one. | |
| 477 | """ | |
| 478 | import gpu | |
| 479 | import numpy as np | |
| 480 | import ifs_host | |
| 481 | import ifs_render | |
| 482 | ||
| 483 | camera = doc.camera | |
| 484 | camera.view_proj = ifs_render.view_projection( | |
| 485 | camera.position[:3], camera.forward[:3], camera.up, | |
| 486 | fov_degrees=camera.field_of_view) | |
| 487 | ||
| 488 | kernel, drain, _library = ifs_render._build_shaders( | |
| 489 | size, size, ifs_render.DEFAULT_HIST_SCALE) | |
| 490 | ubos = { | |
| 491 | 'settings': gpu.types.GPUUniformBuf( | |
| 492 | ifs_host.pack_settings(camera, doc.settings)), | |
| 493 | 'iterators_buf': gpu.types.GPUUniformBuf( | |
| 494 | ifs_host.pack_iterators(doc.iterators)), | |
| 495 | 'palette_buf': gpu.types.GPUUniformBuf( | |
| 496 | ifs_host.pack_palette(doc.palette(doc.settings.gradient_resolution))), | |
| 497 | 'real_params_buf': gpu.types.GPUUniformBuf( | |
| 498 | ifs_host.pack_real_params(ifs_host.flatten_real_params(doc.iterators))), | |
| 499 | 'vec3_params_buf': gpu.types.GPUUniformBuf( | |
| 500 | ifs_host.pack_vec3_params(ifs_host.flatten_vec3_params(doc.iterators))), | |
| 501 | } | |
| 502 | alias = ifs_host.build_alias_tables(doc.iterators) | |
| 503 | n = len(doc.iterators) | |
| 504 | alias_tex = gpu.types.GPUTexture( | |
| 505 | (n, n), format='RGBA32F', | |
| 506 | data=gpu.types.Buffer('FLOAT', len(alias), alias)) | |
| 507 | ||
| 508 | histogram = gpu.types.GPUTexture((size * 4, size), format='R32UI') | |
| 509 | histogram.clear(format='UINT', value=(0, 0, 0, 0)) | |
| 510 | accum = gpu.types.GPUTexture((size, size), format='RGBA32F') | |
| 511 | accum.clear(format='FLOAT', value=(0.0, 0.0, 0.0, 0.0)) | |
| 512 | state_rows = (ifs_render.INVOCATIONS + ifs_render.STATE_W - 1) // ifs_render.STATE_W | |
| 513 | state = gpu.types.GPUTexture( | |
| 514 | (ifs_render.STATE_W * 2, state_rows), format='RGBA32F') | |
| 515 | state.clear(format='FLOAT', value=(0.0, 0.0, 0.0, 0.0)) | |
| 516 | ||
| 517 | for index in range(dispatches): | |
| 518 | frame = gpu.types.GPUUniformBuf(ifs_host.pack_frame( | |
| 519 | seed=0x9E3779B9, width=size, height=size, | |
| 520 | dispatch_cnt=index, reset_points_state=1 if index == 0 else 0, | |
| 521 | invocation_iters=iters)) | |
| 522 | kernel.bind() | |
| 523 | kernel.image('histogram_img', histogram) | |
| 524 | kernel.image('state_img', state) | |
| 525 | kernel.uniform_sampler('alias_tables_tex', alias_tex) | |
| 526 | kernel.uniform_block('frame', frame) | |
| 527 | for name, ubo in ubos.items(): | |
| 528 | kernel.uniform_block(name, ubo) | |
| 529 | gpu.compute.dispatch(kernel, ifs_render.WORKGROUPS, 1, 1) | |
| 530 | ||
| 531 | drain.bind() | |
| 532 | drain.image('histogram_img', histogram) | |
| 533 | drain.image('accum_img', accum) | |
| 534 | drain.uniform_int('out_size', (size, size)) | |
| 535 | gpu.compute.dispatch(drain, (size + 7) // 8, (size + 7) // 8, 1) | |
| 536 | ||
| 537 | pixels = np.array(accum.read(), dtype=np.float32).reshape(size, size, 4) | |
| 538 | return pixels[:, :, 3] | |
| 539 | ||
| 540 | ||
| 541 | def projected_density(view_proj, points, size=STILL_SIZE): | |
| 542 | """A point set through the engine's own projection, binned like the render.""" | |
| 543 | import numpy as np | |
| 544 | import ifs_render | |
| 545 | out = np.zeros((size, size), dtype=np.float64) | |
| 546 | for point in points: | |
| 547 | result = ifs_render.project_point( | |
| 548 | view_proj, (float(point[0]), float(point[1]), float(point[2])), | |
| 549 | size, size) | |
| 550 | if result is None: | |
| 551 | continue | |
| 552 | x, y, _w = result | |
| 553 | out[int(round(y)), int(round(x))] += 1.0 | |
| 554 | return out | |
| 555 | ||
| 556 | ||
| 548 | 557 | if __name__ == "__main__": |
| 549 | 558 | result = unittest.main(argv=sys.argv[:1], exit=False, verbosity=2).result |
| 550 | 559 | sys.exit(0 if result.wasSuccessful() else 1) |
| … | … | @@ -0,0 +1,293 @@ |
| 1 | """ | |
| 2 | Proves the emitted point cloud IS the fractal — same ground truth as the bake. | |
| 3 | ||
| 4 | `test_bake_orientation` established a document whose attractor is known | |
| 5 | analytically rather than measured: four Affine iterators `p := 0.5p + 0.5F_i`, so | |
| 6 | the attractor is the Sierpinski tetrahedron with vertices `F_i`, its convex hull is | |
| 7 | exactly that tetrahedron and its centroid is exactly `mean(F)`. This reuses it, so | |
| 8 | the volume path and the points path are judged against one reference instead of | |
| 9 | each against its own. | |
| 10 | ||
| 11 | The tetrahedron pays a second dividend here. Four maps at ratio 0.5 give it a | |
| 12 | Hausdorff dimension of exactly `ln 4 / ln 2 = 2`, which is a real prediction for | |
| 13 | `estimate_spacing` — the radius derivation measures the attractor's dimension on | |
| 14 | the way to a spacing, so there is a known number for it to hit rather than a | |
| 15 | plausible-looking one to accept. | |
| 16 | ||
| 17 | Three failure modes this is built to catch, none of which look like a failure: | |
| 18 | ||
| 19 | * **A truncated cloud.** The array fills and the counter keeps going; a cloud of the | |
| 20 | first million of twelve million samples is still a perfectly convincing fractal. | |
| 21 | * **Uninitialised texels.** Reading past the counter yields (0, 0, 0), and the origin | |
| 22 | sits INSIDE this attractor's hull, so a hull test alone would not notice. The | |
| 23 | kernel treats `dot(pos, pos) == 0` as an invalid point and never emits it, so an | |
| 24 | exact zero in the output can only be a slot that was never written. | |
| 25 | * **A cloud that is not what the renderer draws.** The emitted points are projected | |
| 26 | through the engine's own camera and compared to the actual still render. | |
| 27 | ||
| 28 | Needs the GPU, so it runs inside Blender: | |
| 29 | blender -b --factory-startup -noaudio --python blender/test_points_emit.py | |
| 30 | """ | |
| 31 | ||
| 32 | import json | |
| 33 | import os | |
| 34 | import sys | |
| 35 | import tempfile | |
| 36 | import unittest | |
| 37 | ||
| 38 | sys.path.insert(0, os.path.dirname(os.path.abspath(__file__))) | |
| 39 | ||
| 40 | from test_bake_orientation import ( | |
| 41 | FIXED_POINTS, | |
| 42 | STILL_SIZE, | |
| 43 | _cpu_attractor, | |
| 44 | _document_dict, | |
| 45 | _mean, | |
| 46 | ) | |
| 47 | ||
| 48 | CAPACITY = 200_000 | |
| 49 | SMALL_CAPACITY = 8192 # exactly two rows of the points image | |
| 50 | ||
| 51 | # The single white key in the orientation document would make every point the same | |
| 52 | # colour, which cannot tell "colour is carried" from "colour is a constant". | |
| 53 | PALETTE = [ | |
| 54 | {"X": 1.0, "Y": 0.2, "Z": 0.1, "W": 1.0}, | |
| 55 | {"X": 0.2, "Y": 1.0, "Z": 0.3, "W": 1.0}, | |
| 56 | {"X": 0.1, "Y": 0.3, "Z": 1.0, "W": 1.0}, | |
| 57 | ] | |
| 58 | ||
| 59 | ||
| 60 | def _outside_hull(points, vertices, tolerance): | |
| 61 | """How many points fall outside the tetrahedron, and by how far at worst.""" | |
| 62 | import numpy as np | |
| 63 | verts = np.array(vertices, dtype=np.float64) | |
| 64 | outside, worst = 0, 0.0 | |
| 65 | for corner in range(4): | |
| 66 | face = np.array([verts[j] for j in range(4) if j != corner]) | |
| 67 | normal = np.cross(face[1] - face[0], face[2] - face[0]) | |
| 68 | normal = normal / np.linalg.norm(normal) | |
| 69 | offset = float(normal @ face[0]) | |
| 70 | # Orient the plane so "inside" is the side the fourth vertex is on. | |
| 71 | inside = 1.0 if (normal @ verts[corner] - offset) > 0 else -1.0 | |
| 72 | signed = inside * ((points @ normal) - offset) | |
| 73 | outside += int((signed < -tolerance).sum()) | |
| 74 | worst = min(worst, float(signed.min())) | |
| 75 | return outside, worst | |
| 76 | ||
| 77 | ||
| 78 | class Emit(unittest.TestCase): | |
| 79 | @classmethod | |
| 80 | def setUpClass(cls): | |
| 81 | try: | |
| 82 | import numpy # noqa: F401 | |
| 83 | import ifs_points # noqa: F401 | |
| 84 | except ImportError as ex: | |
| 85 | raise unittest.SkipTest(f"needs Blender's gpu module ({ex})") | |
| 86 | ||
| 87 | import ifs_document | |
| 88 | import ifs_points | |
| 89 | import ifs_render | |
| 90 | ||
| 91 | document = _document_dict() | |
| 92 | document["Palette"] = {"Colors": PALETTE} | |
| 93 | cls.tmp = tempfile.mkdtemp() | |
| 94 | cls.doc_path = os.path.join(cls.tmp, "points.ifsjson") | |
| 95 | with open(cls.doc_path, "w", encoding="utf-8") as handle: | |
| 96 | json.dump(document, handle) | |
| 97 | ||
| 98 | library = ifs_render.load_library( | |
| 99 | os.path.join(ifs_render.REPO, "Transforms"), | |
| 100 | os.path.join(ifs_render.REPO, "WpfDisplay", "Library", "Includes")) | |
| 101 | cls.doc = ifs_document.load(cls.doc_path, library) | |
| 102 | cls.positions, cls.colors, cls.stats = ifs_points.emit_document( | |
| 103 | cls.doc, capacity=CAPACITY) | |
| 104 | print(f"[points] {cls.stats}") | |
| 105 | ||
| 106 | # ------------------------------------------------- is it the right fractal | |
| 107 | ||
| 108 | def test_the_points_land_on_the_attractor(self): | |
| 109 | import numpy as np | |
| 110 | expected = np.array([_mean(FIXED_POINTS, a) for a in range(3)]) | |
| 111 | measured = self.positions.mean(axis=0).astype(np.float64) | |
| 112 | print(f"[points] centroid {np.round(measured, 5)} expected " | |
| 113 | f"{np.round(expected, 5)}") | |
| 114 | for axis in range(3): | |
| 115 | self.assertAlmostEqual( | |
| 116 | expected[axis], measured[axis], delta=0.01, | |
| 117 | msg=f"axis {axis}: emitted centroid {measured[axis]:.5f} against " | |
| 118 | f"the analytic {expected[axis]:.5f}") | |
| 119 | ||
| 120 | def test_no_point_is_outside_the_convex_hull(self): | |
| 121 | """ | |
| 122 | The attractor is contained in the tetrahedron by construction, so anything | |
| 123 | outside it is not a sample of this fractal. The tolerance is float slop, not | |
| 124 | latitude — the maps are exact contractions toward the vertices. | |
| 125 | """ | |
| 126 | import numpy as np | |
| 127 | outside, worst = _outside_hull( | |
| 128 | self.positions.astype(np.float64), FIXED_POINTS, 1e-4) | |
| 129 | print(f"[points] outside hull: {outside} of {len(self.positions)}, " | |
| 130 | f"worst {worst:.2e}") | |
| 131 | self.assertEqual(0, outside, | |
| 132 | f"{outside} points lie outside the attractor's hull, " | |
| 133 | f"the worst by {abs(worst):.4f}") | |
| 134 | ||
| 135 | def test_no_point_sits_at_the_origin(self): | |
| 136 | """ | |
| 137 | `ivec3` of an unwritten texel is (0, 0, 0), and the origin is INSIDE this | |
| 138 | hull, so a slot never written would pass every other check here. The kernel | |
| 139 | rejects `dot(pos, pos) == 0` as an invalid point, so a zero can only mean the | |
| 140 | readback ran past the counter. | |
| 141 | """ | |
| 142 | import numpy as np | |
| 143 | zeros = int((np.abs(self.positions).sum(axis=1) == 0.0).sum()) | |
| 144 | self.assertEqual(0, zeros, | |
| 145 | f"{zeros} points are exactly at the origin — the readback " | |
| 146 | "is including slots the kernel never wrote") | |
| 147 | ||
| 148 | def test_the_cloud_reproduces_the_still_render(self): | |
| 149 | """ | |
| 150 | The checkpoint the plan asked for: does the emitted cloud have the same shape, | |
| 151 | in the same place, as what the renderer draws? Projected through the engine's | |
| 152 | own camera and compared to an actual GPU still render. | |
| 153 | """ | |
| 154 | import numpy as np | |
| 155 | from test_bake_orientation import projected_density, still_density | |
| 156 | ||
| 157 | rendered = still_density(self.doc) | |
| 158 | # A subsample keeps the pure-Python projection loop honest about runtime; | |
| 159 | # 60k points over a 192x192 frame still covers every lit pixel many times. | |
| 160 | index = np.random.default_rng(7).choice( | |
| 161 | len(self.positions), size=min(60_000, len(self.positions)), | |
| 162 | replace=False) | |
| 163 | predicted = projected_density( | |
| 164 | self.doc.camera.view_proj, self.positions[index]) | |
| 165 | ||
| 166 | def centroid(image): | |
| 167 | ys, xs = np.nonzero(image) | |
| 168 | weights = image[ys, xs] | |
| 169 | total = weights.sum() | |
| 170 | return np.array([(xs * weights).sum() / total, | |
| 171 | (ys * weights).sum() / total]) | |
| 172 | ||
| 173 | got, want = centroid(predicted), centroid(rendered) | |
| 174 | overlap = np.logical_and(rendered > 0, predicted > 0).sum() | |
| 175 | union = np.logical_or(rendered > 0, predicted > 0).sum() | |
| 176 | iou = overlap / float(union) | |
| 177 | print(f"[points] projected centroid {np.round(got, 2)} vs still render " | |
| 178 | f"{np.round(want, 2)} px, IoU {iou:.3f}") | |
| 179 | ||
| 180 | for axis, name in enumerate(("x", "y")): | |
| 181 | self.assertLess( | |
| 182 | abs(got[axis] - want[axis]), 4.0, | |
| 183 | f"screen {name}: the cloud projects to {got[axis]:.1f} px and the " | |
| 184 | f"still render sits at {want[axis]:.1f} px") | |
| 185 | self.assertGreater( | |
| 186 | iou, 0.5, | |
| 187 | f"only {iou:.1%} of lit pixels agree — the cloud is not the shape the " | |
| 188 | "renderer draws") | |
| 189 | ||
| 190 | # --------------------------------------------------- is it carried correctly | |
| 191 | ||
| 192 | def test_colour_is_carried_per_point(self): | |
| 193 | import numpy as np | |
| 194 | rgb = self.colors[:, :3] | |
| 195 | self.assertTrue(np.isfinite(self.colors).all(), "colour has NaN or inf") | |
| 196 | self.assertGreaterEqual(float(rgb.min()), 0.0) | |
| 197 | self.assertLessEqual(float(rgb.max()), 1.0) | |
| 198 | spread = float(rgb.std(axis=0).max()) | |
| 199 | print(f"[points] colour spread {spread:.4f}, mean {np.round(rgb.mean(0), 3)}") | |
| 200 | self.assertGreater( | |
| 201 | spread, 0.01, | |
| 202 | "every point has the same colour — the palette is not reaching the emit") | |
| 203 | # The w channel carries the iterator's opacity, which is 1.0 throughout here. | |
| 204 | self.assertAlmostEqual(1.0, float(self.colors[:, 3].min()), places=5) | |
| 205 | ||
| 206 | def test_the_dimension_estimate_recovers_the_known_value(self): | |
| 207 | """ | |
| 208 | Four maps at ratio 0.5 put the Sierpinski tetrahedron at ln4/ln2 = 2 exactly. | |
| 209 | The radius derivation measures dimension on its way to a spacing, so there is | |
| 210 | a right answer for it to land on — which is the difference between a derived | |
| 211 | number and a plausible one. | |
| 212 | """ | |
| 213 | import ifs_points | |
| 214 | spacing, dimension = ifs_points.estimate_spacing(self.positions) | |
| 215 | print(f"[points] spacing {spacing:.6f}, dimension {dimension:.3f} " | |
| 216 | f"(exact 2.0)") | |
| 217 | self.assertAlmostEqual( | |
| 218 | 2.0, dimension, delta=0.3, | |
| 219 | msg=f"measured dimension {dimension:.3f} against the tetrahedron's " | |
| 220 | "exact 2.0") | |
| 221 | self.assertGreater(spacing, 0.0) | |
| 222 | # A spacing wider than the fractal itself would mean one blob per point. | |
| 223 | self.assertLess(spacing, 0.1) | |
| 224 | ||
| 225 | # ------------------------------------------------------ capacity and honesty | |
| 226 | ||
| 227 | def test_capacity_is_respected_and_overflow_reported(self): | |
| 228 | """ | |
| 229 | A truncated cloud looks exactly like a complete one, so the host has to say | |
| 230 | so. The counter runs past the end on purpose to make that possible. | |
| 231 | """ | |
| 232 | import ifs_points | |
| 233 | _positions, _colors, stats = ifs_points.emit_document( | |
| 234 | self.doc, capacity=SMALL_CAPACITY) | |
| 235 | print(f"[points] small emit {stats}") | |
| 236 | self.assertEqual(SMALL_CAPACITY, stats["capacity"]) | |
| 237 | self.assertEqual(SMALL_CAPACITY, stats["points"]) | |
| 238 | self.assertEqual(SMALL_CAPACITY, len(_positions)) | |
| 239 | self.assertTrue(stats["filled"]) | |
| 240 | self.assertGreaterEqual(stats["requested"], stats["capacity"]) | |
| 241 | ||
| 242 | def test_the_arrays_agree_with_the_reported_count(self): | |
| 243 | self.assertEqual(self.stats["points"], len(self.positions)) | |
| 244 | self.assertEqual(self.stats["points"], len(self.colors)) | |
| 245 | self.assertEqual(3, self.positions.shape[1]) | |
| 246 | self.assertEqual(4, self.colors.shape[1]) | |
| 247 | ||
| 248 | def test_the_emitted_cloud_matches_a_cpu_chaos_game(self): | |
| 249 | """Per-axis mass distribution against an independent CPU run of the same IFS.""" | |
| 250 | import numpy as np | |
| 251 | reference = _cpu_attractor() | |
| 252 | for axis in range(3): | |
| 253 | for q in (5, 50, 95): | |
| 254 | got = float(np.percentile(self.positions[:, axis], q)) | |
| 255 | want = float(np.percentile(reference[:, axis], q)) | |
| 256 | self.assertAlmostEqual( | |
| 257 | want, got, delta=0.02, | |
| 258 | msg=f"axis {axis} p{q:02d}: emitted {got:+.4f} against the CPU " | |
| 259 | f"attractor's {want:+.4f}") | |
| 260 | ||
| 261 | # ------------------------------------------------------------------ ingest | |
| 262 | ||
| 263 | def test_the_pointcloud_ingests_into_blender(self): | |
| 264 | import numpy as np | |
| 265 | import ifs_points | |
| 266 | take = 5000 | |
| 267 | positions = self.positions[:take] | |
| 268 | colors = self.colors[:take] | |
| 269 | cloud = ifs_points.to_pointcloud(positions, colors, radius=0.001, | |
| 270 | name="IFS Points Test") | |
| 271 | self.assertEqual(take, len(cloud.points)) | |
| 272 | ||
| 273 | out = np.zeros(take * 3, dtype=np.float32) | |
| 274 | cloud.attributes['position'].data.foreach_get('vector', out) | |
| 275 | self.assertLess(float(np.abs(out.reshape(take, 3) - positions).max()), 1e-6) | |
| 276 | ||
| 277 | self.assertIn(ifs_points.COLOR_ATTRIBUTE, cloud.attributes) | |
| 278 | self.assertNotIn( | |
| 279 | "color", [a.name for a in cloud.attributes], | |
| 280 | "`color` is a poisoned attribute name for points — it renders black") | |
| 281 | got = np.zeros(take * 4, dtype=np.float32) | |
| 282 | cloud.attributes[ifs_points.COLOR_ATTRIBUTE].data.foreach_get('color', got) | |
| 283 | self.assertLess( | |
| 284 | float(np.abs(got.reshape(take, 4)[:, :3] - colors[:, :3]).max()), 1e-6) | |
| 285 | ||
| 286 | radii = np.zeros(take, dtype=np.float32) | |
| 287 | cloud.attributes['radius'].data.foreach_get('value', radii) | |
| 288 | self.assertAlmostEqual(0.001, float(radii.min()), places=6) | |
| 289 | ||
| 290 | ||
| 291 | if __name__ == "__main__": | |
| 292 | result = unittest.main(argv=sys.argv[:1], exit=False, verbosity=2).result | |
| 293 | sys.exit(0 if result.wasSuccessful() else 1) |
| … | … | @@ -0,0 +1,105 @@ |
| 1 | """ | |
| 2 | Tests for the POINTS variant of the kernel rewrite. | |
| 3 | ||
| 4 | Volume mode swaps the plot block for a voxel write; points mode swaps it for an | |
| 5 | append to a flat list. The list is the interesting part, because Blender's GPU | |
| 6 | module has no SSBOs (`ifs_layout`'s module docstring, reason 1) — so there is no | |
| 7 | buffer to append to. The substitute is an RGBA32F image used as an array, one | |
| 8 | texel per point, indexed by an atomic counter in a 1x1 R32UI image. That is the | |
| 9 | same trick the viewport's auto-exposure max pass uses, turned from a reduction | |
| 10 | into an allocator. | |
| 11 | ||
| 12 | What can go wrong here is not that the arithmetic is subtly off; it is that a | |
| 13 | substitution silently stops matching and leaves the 2D projection in place, so | |
| 14 | the kernel still compiles and still runs and quietly plots to a histogram nobody | |
| 15 | reads. So these assert what the body must and must NOT contain, and that the | |
| 16 | three modes stay mutually exclusive — a kernel that tried to be two of them at | |
| 17 | once would bind images under two different meanings. | |
| 18 | ||
| 19 | Run: python3 -m unittest discover -s blender | |
| 20 | """ | |
| 21 | ||
| 22 | import os | |
| 23 | import unittest | |
| 24 | ||
| 25 | from ifs_codegen import apply_portability_fixes, generate_kernel, load_library | |
| 26 | from ifs_layout import POINTS_LAYOUT, rewrite | |
| 27 | ||
| 28 | REPO = os.path.dirname(os.path.dirname(os.path.abspath(__file__))) | |
| 29 | TRANSFORMS = os.path.join(REPO, "Transforms") | |
| 30 | INCLUDES = os.path.join(REPO, "WpfDisplay", "Library", "Includes") | |
| 31 | TEMPLATE = os.path.join(REPO, "IFSEngine", "Rendering", "Shaders", "ifs_kernel.comp.shader") | |
| 32 | ||
| 33 | ||
| 34 | class PointsRewrite(unittest.TestCase): | |
| 35 | @classmethod | |
| 36 | def setUpClass(cls): | |
| 37 | if not os.path.isfile(TEMPLATE): | |
| 38 | raise unittest.SkipTest("kernel template not found") | |
| 39 | library = load_library(TRANSFORMS, INCLUDES) | |
| 40 | with open(TEMPLATE, "r", encoding="utf-8-sig") as handle: | |
| 41 | kernel = generate_kernel(handle.read(), library) | |
| 42 | kernel, _ = apply_portability_fixes(kernel) | |
| 43 | cls.body, cls.typedefs, cls.report = rewrite(kernel, points=True) | |
| 44 | ||
| 45 | def test_every_substitution_applied(self): | |
| 46 | self.assertTrue(self.report.ok, | |
| 47 | "rewrite incomplete: " + "; ".join(self.report.missing)) | |
| 48 | ||
| 49 | def test_the_projection_is_gone_entirely(self): | |
| 50 | """ | |
| 51 | Not bypassed — REPLACED. `project()` can `continue` on out-of-frame, so a | |
| 52 | points mode that merely ignored the result would still drop every sample | |
| 53 | outside the current view and emit a view-dependent cloud. That is the same | |
| 54 | reasoning volume mode is built on. | |
| 55 | """ | |
| 56 | self.assertNotIn("projf = project(", self.body) | |
| 57 | self.assertNotIn("accumulate_hit", self.body) | |
| 58 | ||
| 59 | def test_points_are_appended_through_the_atomic_counter(self): | |
| 60 | self.assertIn("emit_point", self.body) | |
| 61 | self.assertIn("imageAtomicAdd(counter_img", self.body) | |
| 62 | self.assertIn("imageStore(points_img", self.body) | |
| 63 | self.assertIn("imageStore(colors_img", self.body) | |
| 64 | ||
| 65 | def test_the_emit_is_bounded_by_capacity(self): | |
| 66 | """ | |
| 67 | The counter keeps counting past the end of the image on purpose, so the host | |
| 68 | can report how many points were WANTED. Without a capacity guard that would | |
| 69 | be an out-of-bounds imageStore on every sample past the end. | |
| 70 | """ | |
| 71 | self.assertIn("point_capacity", self.body) | |
| 72 | self.assertIn("points_width", self.body) | |
| 73 | ||
| 74 | def test_colour_is_not_premultiplied(self): | |
| 75 | """ | |
| 76 | The volume path multiplies colour by opacity because it accumulates a | |
| 77 | density-weighted sum and divides it back out. A point is not accumulated — | |
| 78 | it carries its own colour — so premultiplying would darken every point by | |
| 79 | its iterator's opacity and there is nothing to divide it out again. | |
| 80 | """ | |
| 81 | self.assertNotIn("color.xyz *= color.w", self.body) | |
| 82 | ||
| 83 | def test_the_modes_are_mutually_exclusive(self): | |
| 84 | with self.assertRaises(ValueError): | |
| 85 | rewrite("", points=True, volume=True) | |
| 86 | with self.assertRaises(ValueError): | |
| 87 | rewrite("", points=True, depth=True) | |
| 88 | ||
| 89 | def test_the_layout_fits_blenders_image_budget(self): | |
| 90 | """Blender caps images at 8 and the budget is STATIC — see ifs_layout.""" | |
| 91 | images = [r for r in POINTS_LAYOUT if r.kind == "image"] | |
| 92 | self.assertLessEqual(len(images), 8) | |
| 93 | self.assertEqual(len(images), len({r.slot for r in images}), | |
| 94 | "two images share a binding slot") | |
| 95 | ||
| 96 | def test_the_layout_names_match_the_glsl(self): | |
| 97 | """A declared name that the body never uses binds nothing and wastes a slot.""" | |
| 98 | for resource in POINTS_LAYOUT: | |
| 99 | if resource.kind == "image": | |
| 100 | self.assertIn(resource.name, self.body, | |
| 101 | f"{resource.name} is declared but never used") | |
| 102 | ||
| 103 | ||
| 104 | if __name__ == "__main__": | |
| 105 | unittest.main() |