"""
Emits the fractal as a PointCloud, so Cycles and EEVEE both render it as scene data.
This is the third way of getting a flame out of the chaos game, and the first with
no resolution ceiling.
* The **live viewport** is a screen-space overlay: no depth of its own, invisible
to F12. It is a preview.
* The **volume bake** is real scene data, but a voxel grid is 4N x N x N uint32 on
the GPU - 268 MB at N=256 and 2.1 GB at N=512, which does not fit. Detail is
bought cubically, so a flame's thin filaments end up as mush.
* **Points** store what the chaos game actually produces: samples, at float
precision. Nothing is quantised, so there is no grid to be finer than. The cost
is linear in the number of samples rather than cubic in the resolution, and the
same representation renders in Cycles properly and in EEVEE as an in-scene
approximation - one representation, both engines.
## How the points get off the GPU
There is nowhere to append them to: Blender's GPU module has no SSBOs. The
substitute is an RGBA32F image used as a flat array, one texel per point, with an
atomic counter in a 1x1 R32UI handing out slots - see `ifs_layout`'s points mode.
The host allocates `points_width x ceil(capacity/points_width)`, dispatches until
the counter says the array is full, and reads the two images back.
Because slots are claimed in atomic order, a full buffer holds the FIRST `capacity`
samples of the run - and WHICH samples those are is a sampling-quality question, not
just a performance one. See `iters_for`: the dispatch is sized so that one pass over
all 16384 chains fills the buffer exactly once, which is what makes those samples
`capacity / 16384` iterations from EVERY chain rather than a long correlated run
from whichever few thousand the GPU scheduled first.
## Radius
Cycles draws each point as a sphere, so the cloud only reads as a continuous object
if the radius is near the local spacing between samples. That spacing is MEASURED
rather than assumed - see `estimate_spacing`, which also falls out as an estimate of
the attractor's fractal dimension. It is a starting point, not a truth: how big a
point should look also depends on how close the camera is, so the operator exposes a
multiplier.
"""
from __future__ import annotations
import os
import gpu
import numpy as np
import ifs_codegen
import ifs_document
import ifs_host
import ifs_layout
import ifs_render
# One texel per point. 4096 is comfortably inside every GL_MAX_TEXTURE_SIZE worth
# worrying about (16384 is the common floor) and keeps the row count small enough
# that a 16M-point buffer is 4096 rows.
POINTS_WIDTH = 4096
# NOT "color". That name is poisoned for points and meshes in Blender - a per-point
# attribute called `color` renders black - while being the EXPECTED name for a
# volume grid. Both halves of that cost a wrong conclusion during the volume work.
COLOR_ATTRIBUTE = "ifs_color"
# The degree-0 spherical-harmonics basis constant, 1/(2*sqrt(pi)).
#
# Verified against THREE primary sources that agree: graphdeco-inria's
# `gaussian_model.py`, antimatter15/splat's `convert.py`, and Blender's own EEVEE
# splat shader on the `gsplat` branch, which evaluates
# `srgb_to_linear(0.282094792 * color.rgb + ... + 0.5)`.
SH_C0 = 0.28209479177387814
def encode_sh_dc(colour):
"""
Colour -> the RAW degree-0 SH coefficient that splat renderers store.
Degree-0 SH is not colour, and this is the mapping that catches people out:
the renderer computes `0.5 + C0 * coefficient`, so going the other way needs
BOTH the offset and the constant. An encoder that drops the 0.5 still round
trips against its own decoder while rendering everything washed toward grey,
which is why the test checks a hand-computed value as well as a round trip.
"""
return (np.asarray(colour, dtype=np.float32) - 0.5) / SH_C0
def decode_sh_dc(coefficient):
"""The inverse of `encode_sh_dc`, i.e. what the shader does."""
return 0.5 + SH_C0 * np.asarray(coefficient, dtype=np.float32)
def build_shader():
"""Compiles the points variant of the kernel."""
library = ifs_codegen.load_library(
os.path.join(ifs_render.REPO, "Transforms"),
os.path.join(ifs_render.REPO, "WpfDisplay", "Library", "Includes"))
template_path = os.path.join(
ifs_render.REPO, "IFSEngine", "Rendering", "Shaders", "ifs_kernel.comp.shader")
with open(template_path, "r", encoding="utf-8-sig") as handle:
kernel = ifs_codegen.generate_kernel(handle.read(), library)
kernel, _ = ifs_codegen.apply_portability_fixes(kernel)
body, typedefs, report = ifs_layout.rewrite(
kernel, ifs_layout.DEFAULT_HIST_SCALE, ifs_render.STATE_W, points=True)
if not report.ok:
raise RuntimeError("points rewrite incomplete: " + "; ".join(report.missing))
info = gpu.types.GPUShaderCreateInfo()
info.typedef_source(typedefs)
info.local_group_size(ifs_render.LOCAL_SIZE)
rw = {'NO_RESTRICT', 'READ', 'WRITE'}
info.image(0, 'RGBA32F', 'FLOAT_2D', 'points_img', qualifiers=rw)
info.image(1, 'RGBA32F', 'FLOAT_2D', 'colors_img', qualifiers=rw)
info.image(2, 'R32UI', 'UINT_2D', 'counter_img', qualifiers=rw)
info.image(3, 'RGBA32F', 'FLOAT_2D', 'state_img', qualifiers=rw)
info.sampler(0, 'FLOAT_2D', 'alias_tables_tex')
info.push_constant('INT', 'point_capacity')
info.push_constant('INT', 'points_width')
info.uniform_buf(0, 'FrameParams', 'frame')
info.uniform_buf(1, 'SettingsParams', 'settings')
info.uniform_buf(2, 'IteratorsParams', 'iterators_buf')
info.uniform_buf(3, 'PaletteParams', 'palette_buf')
info.uniform_buf(4, 'RealParams', 'real_params_buf')
info.uniform_buf(5, 'Vec3Params', 'vec3_params_buf')
info.compute_source(body)
return gpu.shader.create_from_info(info), library
def _document_buffers(doc):
"""The five per-document UBOs.
Deliberately a copy of `ifs_volume._document_buffers` rather than an import:
it is mechanical std140 packing, and coupling the two bake modes through a
private name would mean a change for one silently reshapes the other.
"""
return {
'settings': gpu.types.GPUUniformBuf(
ifs_host.pack_settings(doc.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))),
}
def iters_for(capacity, warmup):
"""
How many iterations each invocation should run to fill `capacity` ONCE.
This is a sampling-quality decision, not a performance one. A dispatch at the
engine's default 500 iterations emits ~8M samples from 16384 chains; against a
1M buffer, seven eighths of them hit a full array and are dropped. That would be
merely wasteful if the survivors were a fair sample, but they are not: slots are
claimed in atomic order, so the buffer fills from whichever invocations the GPU
happens to schedule FIRST. Only a few thousand chains are resident at once, so a
500-iteration dispatch fills the array from a fraction of the chains, each
contributing a long correlated run.
Sizing the dispatch to the buffer instead means every one of the 16384 chains
contributes its share - `capacity / 16384` samples each - which is the widest
front the kernel can offer. The warmup is added on top because those iterations
are discarded before anything is plotted.
"""
per_chain = (capacity + ifs_render.INVOCATIONS - 1) // ifs_render.INVOCATIONS
return int(max(warmup + 8, warmup + per_chain))
def emit(shader, doc, ubos, alias_tex, capacity, max_dispatches=400,
invocation_iters=None, poll_every=1):
"""
Runs the chaos game until `capacity` points are stored or the budget runs out.
Returns `(positions, colors, stats)` with positions (N, 3) and colors (N, 4),
both float32 and both in the order the GPU claimed slots.
The counter is polled rather than trusted at the end, because once the array is
full every further dispatch is pure waste and the counter is a uint32 that a
long enough run would wrap. With `invocation_iters` derived from the capacity a
single dispatch very nearly fills the buffer, so polling every dispatch costs
one readback rather than serialising anything.
"""
if invocation_iters is None:
invocation_iters = iters_for(capacity, doc.settings.warmup)
rows = max(1, (capacity + POINTS_WIDTH - 1) // POINTS_WIDTH)
allocated = rows * POINTS_WIDTH
points_tex = gpu.types.GPUTexture((POINTS_WIDTH, rows), format='RGBA32F')
points_tex.clear(format='FLOAT', value=(0.0, 0.0, 0.0, 0.0))
colors_tex = gpu.types.GPUTexture((POINTS_WIDTH, rows), format='RGBA32F')
colors_tex.clear(format='FLOAT', value=(0.0, 0.0, 0.0, 0.0))
counter = gpu.types.GPUTexture((1, 1), format='R32UI')
counter.clear(format='UINT', value=(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))
def claimed():
return int(np.array(counter.read(), dtype=np.uint32).ravel()[0])
dispatched = 0
requested = 0
for index in range(max_dispatches):
frame = gpu.types.GPUUniformBuf(ifs_host.pack_frame(
seed=0x9E3779B9, width=POINTS_WIDTH, height=rows, dispatch_cnt=index,
reset_points_state=1 if index == 0 else 0,
invocation_iters=invocation_iters))
shader.bind()
shader.image('points_img', points_tex)
shader.image('colors_img', colors_tex)
shader.image('counter_img', counter)
shader.image('state_img', state)
shader.uniform_sampler('alias_tables_tex', alias_tex)
shader.uniform_int('point_capacity', (allocated,))
shader.uniform_int('points_width', (POINTS_WIDTH,))
shader.uniform_block('frame', frame)
for name, ubo in ubos.items():
shader.uniform_block(name, ubo)
gpu.compute.dispatch(shader, ifs_render.WORKGROUPS, 1, 1)
dispatched = index + 1
if dispatched % poll_every == 0:
requested = claimed()
if requested >= allocated:
break
requested = max(requested, claimed())
stored = min(requested, allocated)
if stored == 0:
raise RuntimeError(
"no points were emitted - every iterator may have zero opacity, or the "
"document may produce nothing")
positions = np.array(points_tex.read(), dtype=np.float32).reshape(-1, 4)
colors = np.array(colors_tex.read(), dtype=np.float32).reshape(-1, 4)
return positions[:stored, :3].copy(), colors[:stored].copy(), {
"points": stored,
"requested": requested,
"capacity": allocated,
"dispatches": dispatched,
"filled": bool(requested >= allocated),
}
def estimate_spacing(points, sample=2048, seed=0):
"""
The typical distance between neighbouring samples, MEASURED, plus the fractal
dimension it implies.
A radius has to come from somewhere, and "somewhere" is usually a number typed
once and tuned by eye - which is exactly what went wrong with the volume path's
emission scale, where the right value differed by 200x between documents.
The measurement: mean nearest-neighbour distance over a random subsample. That
distance shrinks with sample count as `d ~ n**(-1/D)` for a set of dimension D,
so measuring it at TWO sizes gives D without assuming one, and extrapolating to
the full count gives the spacing of the cloud actually being built. A flame is
typically somewhere between a surface and a solid, so D lands near 2 - and it is
clamped to [1, 3] because a degenerate document (every sample in one spot) would
otherwise produce a nonsense exponent and a nonsense radius.
"""
rng = np.random.default_rng(seed)
total = len(points)
big = min(sample, total)
small = max(4, big // 4)
if total < 16:
span = float(np.ptp(points, axis=0).max()) if total else 1.0
return max(span, 1e-6) / max(total, 1), 2.0
def nearest(count):
index = rng.choice(total, size=count, replace=False)
chosen = points[index].astype(np.float64)
square = (chosen ** 2).sum(axis=1)
# |a-b|^2 via the dot-product identity: an (n, n) matrix instead of the
# (n, n, 3) a broadcast subtraction would build.
distance = square[:, None] + square[None, :] - 2.0 * (chosen @ chosen.T)
np.fill_diagonal(distance, np.inf)
return float(np.sqrt(np.maximum(distance.min(axis=1), 0.0)).mean())
d_big = nearest(big)
d_small = nearest(small)
ratio = d_small / d_big if d_big > 0 else 1.0
if ratio > 1.0 and big > small:
dimension = float(np.log(big / small) / np.log(ratio))
else:
dimension = 2.0
dimension = float(min(3.0, max(1.0, dimension)))
spacing = d_big * (total / float(big)) ** (-1.0 / dimension)
return max(spacing, 1e-9), dimension
def local_radius(points, dimension=None, cells=64, k=0.5, clamp=(0.02, 50.0)):
"""
A radius PER POINT, from the local sample density.
A single global radius cannot work on a fractal, and this is measurable rather
than aesthetic: within one document the nearest-neighbour distance varies by
16-18x between the p10 and p90 point (light_rods 0.90 -> 14.19,
nebula 0.0175 -> 0.322). Any single value therefore leaves the sparse regions
as separate dots while over-covering the dense ones — which is exactly what a
flame looks like when it renders as fuzzy beads instead of filaments.
Exact nearest-neighbour distances are O(n^2) and out of the question at a
million points, so density is read off a uniform grid: bin the points, and
inside a cell of side `s` holding `c` samples of an attractor with dimension
`D`, the typical spacing is `s / c**(1/D)`. That is the same reasoning
`estimate_spacing` uses globally, applied locally and for free — the whole
thing is a bincount.
The grid is fitted to the MASS, not the extent. light_rods spans 3236 units
while its 99.5% body spans ~1167 and its dense sheet far less; gridding the
full extent would put the entire fractal in one cell and measure nothing.
"""
points = np.asarray(points, dtype=np.float32)
count = len(points)
if count == 0:
return np.zeros(0, dtype=np.float32)
if dimension is None:
_spacing, dimension = estimate_spacing(points)
dimension = float(min(3.0, max(1.0, dimension)))
# Percentile bounds, so a handful of far-flung samples cannot define the grid.
lo = np.percentile(points, 0.5, axis=0)
hi = np.percentile(points, 99.5, axis=0)
span = np.maximum(hi - lo, 1e-9)
cell = span / float(cells)
index = np.floor((points - lo) / cell).astype(np.int64)
np.clip(index, 0, cells - 1, out=index)
flat = (index[:, 0] * cells + index[:, 1]) * cells + index[:, 2]
occupancy = np.bincount(flat, minlength=cells ** 3)
# Cell side as a single length: the geometric mean, so a flat fractal whose
# cells are slabs is not judged by its thinnest axis.
side = float(np.exp(np.log(cell).mean()))
per_point = occupancy[flat].astype(np.float32)
spacing = side / np.maximum(per_point, 1.0) ** (1.0 / dimension)
# CALIBRATE the magnitude against measured nearest-neighbour distances.
#
# The grid gets the SHAPE of the density field right and the scale wrong, and
# both halves were measured. Uncalibrated, the median radius came out 0.15x
# (light_rods) and 0.49x (nebula) of half the true nearest-neighbour spacing —
# i.e. still too small to touch, still beads — while correlating 0.64/0.82 with
# the real thing. The error is systematic, not random: cells are slabs on a
# sheet-like fractal, so a cell's geometric-mean side does not describe the
# spacing inside it.
#
# So the grid supplies the variation and an exact measurement supplies the
# scale. `k` is then what it claims to be — 0.5 means neighbouring splats just
# touch — rather than a number to tune by eye.
if count >= 16:
sample = min(2000, count)
pick = np.random.default_rng(0).choice(count, size=sample, replace=False)
chosen = points[pick].astype(np.float64)
square = (chosen ** 2).sum(axis=1)
distance = square[:, None] + square[None, :] - 2.0 * (chosen @ chosen.T)
np.fill_diagonal(distance, np.inf)
measured = np.sqrt(np.maximum(distance.min(axis=1), 0.0))
target = float(np.median(measured))
here = float(np.median(spacing[pick]))
if here > 0 and np.isfinite(target) and target > 0:
spacing = spacing * (target / here)
radius = (k * spacing).astype(np.float32)
typical = float(np.median(radius)) if count else 1.0
return np.clip(radius, typical * clamp[0], typical * clamp[1]).astype(np.float32)
def emit_document(doc, capacity=1_000_000, max_dispatches=400,
invocation_iters=None, shader=None):
"""
Emits a document that is ALREADY loaded, which is what makes the cloud follow
what is on screen.
Same contract as `ifs_volume.bake_document`, and for the same reason: the node
graph and the panel sliders write into the live session's document, so anything
that reloaded from disk would silently discard every edit made inside Blender.
The document's camera is unused - points mode replaces the projection block
outright, so nothing here depends on where the camera is looking.
"""
gpu.init()
if shader is None:
shader, _library = build_shader()
ubos = _document_buffers(doc)
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))
positions, colors, stats = emit(
shader, doc, ubos, alias_tex, capacity, max_dispatches, invocation_iters)
spacing, dimension = estimate_spacing(positions)
lo = positions.min(axis=0)
hi = positions.max(axis=0)
stats.update({
"spacing": spacing,
# Spheres of half the spacing just touch, which is the smallest radius that
# still reads as a continuous surface rather than dust.
"radius": 0.5 * spacing,
"dimension": dimension,
"bounds": (tuple(float(v) for v in lo), tuple(float(v) for v in hi)),
})
return positions, colors, stats
def emit_from_file(document_path, capacity=1_000_000, max_dispatches=400,
invocation_iters=None):
"""Loads a document from disk and emits it. Returns (positions, colors, stats)."""
gpu.init()
shader, library = build_shader()
doc = ifs_document.load(document_path, library)
return emit_document(doc, capacity, max_dispatches, invocation_iters,
shader=shader)
def to_pointcloud(positions, colors, radius, name="IFS Points",
splat_contract=True):
"""
Builds a Blender PointCloud from the emitted arrays.
Everything goes through the ATTRIBUTE api rather than the `points` collection:
measured at 1M points, `attributes['position'].data.foreach_set('vector', ...)`
takes 1 ms where `points.foreach_set('co', ...)` takes 46 and radius through the
collection takes 90. Same data, same result, forty times cheaper - the
collection path wraps every element in RNA on the way through.
"""
import bpy
count = len(positions)
cloud = bpy.data.pointclouds.new(name)
cloud.resize(count) # creates the `position` attribute
cloud.attributes['position'].data.foreach_set(
'vector', np.ascontiguousarray(positions, dtype=np.float32).ravel())
def attribute(name, type_):
existing = cloud.attributes.get(name)
if existing is not None:
return existing
return cloud.attributes.new(name=name, type=type_, domain='POINT')
colour = attribute(COLOR_ATTRIBUTE, 'FLOAT_COLOR')
rgba = np.ones((count, 4), dtype=np.float32)
rgba[:, :colors.shape[1]] = colors[:, :4]
# The colour's alpha stays 1. Opacity rides in the w channel on the GPU side,
# but it belongs in its own attribute (below) rather than in the colour: the
# material multiplies nothing by this alpha, so letting w leak in here would
# dim the emission for reasons no slider explains.
opacity = (colors[:, 3].astype(np.float32).copy() if colors.shape[1] > 3
else np.ones(count, dtype=np.float32))
rgba[:, 3] = 1.0
colour.data.foreach_set('color', rgba.ravel())
# Scalar or per-point. `local_radius` returns one value per sample, because a
# single radius cannot span the 16-18x density range inside one fractal.
radius_per_point = np.broadcast_to(
np.asarray(radius, dtype=np.float32), (count,)).astype(np.float32)
radii = attribute('radius', 'FLOAT')
radii.data.foreach_set('value', radius_per_point)
# ── the gaussian-splat attributes ───────────────────────────────────────
#
# Blender's native splat work (design task #159470) represents a splat as a
# POINT CLOUD carrying `position`, `scale`, `rotation` (quaternion) and
# `opacity` — the container this function already builds. Writing them costs
# one array each and makes the cloud readable as a splat by that path when it
# lands, and by the addons that exist now, while changing nothing about how
# it renders as points today.
#
# Two of the three are not new measurements, they are data we were DISCARDING:
# `opacity` is the w the kernel already emits (0.5..1.0 across
# light_rods_sample, so genuinely per-point), and `scale` is the radius that
# was already derived from the cloud's own measured nearest-neighbour spacing.
#
# `rotation` is the identity, and that is the correct value rather than a
# placeholder: these splats are isotropic. An oriented splat would need a
# covariance we do not compute — the chaos game gives us sample positions,
# not local anisotropy.
#
# NOT written here: `sh_0_0`. Degree-0 spherical harmonics are not raw colour
# — the 3DGS convention stores a coefficient offset and scaled by the SH basis
# constant — and writing that mapping from memory is how you get a file that
# loads and looks wrong. Colour stays in `ifs_color` until that encoding is
# confirmed against a primary source.
attribute('opacity', 'FLOAT').data.foreach_set('value', opacity)
scale = np.empty((count, 3), dtype=np.float32)
scale[:] = radius_per_point[:, None]
attribute('scale', 'FLOAT_VECTOR').data.foreach_set('vector', scale.ravel())
# (w, x, y, z) — Blender's quaternion order is w-first, and so is the 3DGS
# PLY convention, so the identity is 1 followed by three zeros in both.
rotation = np.zeros((count, 4), dtype=np.float32)
rotation[:, 0] = 1.0
attribute('rotation', 'QUATERNION').data.foreach_set('value', rotation.ravel())
# ── Blender's implemented splat contract, opt-in ────────────────────────
#
# Read from the `gsplat` branch (2026-08-18) rather than from design task
# #159470's prose, which they DISAGREE on. What both its EEVEE and Cycles
# paths actually read is:
#
# color float4 .rgb = the RAW degree-0 SH coefficient (f_dc),
# .a = opacity, sigmoid ALREADY applied
# scale float3 activated (exp already applied) <- ours already is
# rotation quat w-first <- ours already is
# sh_0..14 float3 degrees 1-3, flat running index
#
# There is NO `opacity` attribute and no `sh_0_0` in the implementation.
#
# ON BY DEFAULT, and the "poisoned name" worry was MEASURED rather than
# inherited. `color` renders black when a MATERIAL READS it; merely carrying
# the attribute changes nothing. A/B on a cloud known to render, same material
# reading `ifs_color`, with and without this block:
#
# without color: lit=69 max=0.8157 mean=0.00103
# with color: lit=69 max=0.8157 mean=0.00103 (byte-identical)
#
# So the two contracts do not actually collide, and writing both means the
# cloud is renderable by the points path today and by the splat path when it
# ships. The invariant that DOES matter is that our material keeps reading
# `ifs_color` — test_points_operator pins that.
#
# `opacity` is kept alongside regardless: Blender does not read it, but it is
# OUR per-point data surfaced under an obvious name, and it costs one array.
if splat_contract:
splat_colour = np.empty((count, 4), dtype=np.float32)
splat_colour[:, :3] = encode_sh_dc(colors[:, :3])
splat_colour[:, 3] = opacity
attribute('color', 'FLOAT_COLOR').data.foreach_set(
'color', splat_colour.ravel())
return cloud