The Sunflower Algorithm (which is not a sunflower)

While testing the GPU particle system in Detramotor I needed something to put in the CPU emitter path. The GPU path simulates particles with a compute shader, but the CPU path just hands you a buffer and says "fill this with a million positions, good luck". Uploading a million random points looks like TV static, so I wanted something with structure.

I sat down for about twenty minutes and wrote this:

void generateIterate(SunflowerAlgorithmParameters parameters,
                     CGPUEmitter::Particle*& particleBuffer,
                     int& particleCount,
                     const CGPUEmitter::Particle& head,
                     int recursionDepth) {
  if (recursionDepth <= 0) return;
  if (particleCount <= 0) return;

  SunflowerAlgorithmParameters childParameters = parameters;
  childParameters.alpha = childParameters.alpha * parameters.beta;

  for (int i = 0; i < parameters.childCount && particleCount > 0; i++) {
    CGPUEmitter::Particle child = head;

    child.position = head.position + randomDirection() * parameters.alpha;
    child.color    = glm::packUnorm4x8(glm::unpackUnorm4x8(child.color) +
                     glm::vec4(randomDirection() * parameters.alpha, 0.0f));

    (*particleBuffer++) = child;
    particleCount--;

    generateIterate(childParameters, particleBuffer, particleCount, child, recursionDepth - 1);
  }
}

That is the entire thing. Take a particle, spawn childCount children around it at distance alpha, then recurse on each child with alpha multiplied by beta. Stop when you run out of particles or out of depth.

I called it the Sunflower generator because the first render looked vaguely like the seed head of a sunflower. That name is wrong (more on that in a second) but it stuck, and I am not renaming variables in a codebase for aesthetic honesty.

Here is the thing running in your browser. Drag to orbit, scroll to zoom, and mess with the sliders:

First, the naming apology

There is a classical "sunflower" algorithm in computer graphics, and it is not this. The classical one is Vogel's phyllotaxis model, which places seed n at

$$r = c\sqrt{n}, \qquad \theta = n \cdot 137.508°$$

That golden-angle spiral is what actual sunflowers do, and it is deterministic, 2D, and has nothing to do with recursion. I arrived at my thing from intuition without checking prior art, and picked the name because of what the render looked like. Classic.

So what did I actually write? Structurally it is a branching random walk with geometric self-similarity. If you want to place it on the family tree, it sits much closer to a stochastic Iterated Function System (think L-systems, or the Barnsley fern) than to anything phyllotactic.

The math it is actually doing

Each child is placed at

$$p_{child} = p_{parent} + \alpha_d \cdot u$$

where $u$ is a random unit direction and $\alpha_d$ is the step size at generation $d$:

$$\alpha_d = \alpha_0 \cdot \beta^d$$

Now unroll the recursion. A particle sitting at depth $d$ got there by taking $d+1$ steps, each one in an independent random direction, each one shorter than the last:

$$P_d = \sum_{k=0}^{d} \alpha_0 \beta^k u_k$$

That is a sum of independently-directed, geometrically-shrinking random steps. And this is where it gets nice, because the moment you write it that way you can bound it with the triangle inequality:

$$\|P_d\| \leq \sum_{k=0}^{d} \alpha_0 \beta^k = \alpha_0 \frac{1 - \beta^{d+1}}{1 - \beta} \xrightarrow[d \to \infty]{} \frac{\alpha_0}{1 - \beta}$$

If $\beta < 1$ the whole structure is bounded no matter how deep you recurse. You can ask for a million particles, ten million, whatever, and the cloud never grows past $\alpha_0 / (1-\beta)$. With the demo defaults ($\alpha_0 = 5$, $\beta = 0.72$) that ceiling is about 17.9 units. Ever. Drag the particle slider from 500 to 60000 and watch the "radius bound" readout barely move, because it is already sitting near its limit.

This is the opposite of what a plain random walk does. A plain random walk (constant step size, no decay) has a variance that grows without bound, so it diffuses outward forever and its radius creeps up like $\sqrt{d}$. Mine curls up on itself instead. Each generation adds detail at half the scale of the previous one and the thing converges rather than escapes.

Here is that difference plotted. The blue line is the worst case bound (every step somehow lands in the same direction), the orange line is the mean radius measured over 400 actual walks. Slide $\beta$ toward 1 and watch the ceiling run away:

Notice the measured mean sits way below the bound. That is because the directions are random, so most of the time consecutive steps partly cancel. The bound is a real guarantee, just a pessimistic one.

Fractal dimension

Every node spawns childCount copies of itself scaled by $\beta$. That is textbook self-similar IFS geometry, which means the similarity (Hausdorff) dimension is just

$$D = \frac{\ln(childCount)}{\ln(1/\beta)}$$

The demo above prints this live as you drag the sliders. Some values to get a feel for it:

Two knobs, one for how much new structure gets added per level and one for how fast the perturbation shrinks. The ratio of their logs is the dimension. That is the whole story of the shape.

The low-β cases are worth staring at for a second, because they look wrong the first time. At $\beta = 0.5$ each generation halves the step, so by generation 9 the step is $5 \cdot 2^{-9} \approx 0.01$ while the tree has $3^9$ nodes to place. Almost every particle you generated ends up crammed into a tiny ball at the tip of a branch. The picture is mostly empty space and a few dense specks, which is exactly what dimension 1.58 in a 3D box should look like.

When $D$ comes out above 3 the formula is telling you the set wants more room than it has. In 3D you just get a solid-looking cloud, since nothing can be more space-filling than space.

Picking the depth

int recursionDepth = std::log(particleCount) / std::log(parameters.childCount);

This solves $childCount^d = particleCount$ for $d$, which is the depth of a complete childCount-ary tree holding that many nodes. The point is to spend the whole particle budget evenly across the tree. Get this wrong in one direction and the recursion terminates while you still have half your buffer unwritten. Get it wrong in the other and the first branch eats the entire budget before the second branch ever starts, and you end up rendering one very detailed twig.

The particleCount > 0 check inside the loop is what saves you from the second case, but only in the sense that it stops you from writing past the end of the buffer. It does not make the picture good. The depth formula is what makes the picture good.

Color rides the same walk

This is my favourite part, and it was an accident:

child.color = glm::packUnorm4x8(glm::unpackUnorm4x8(child.color) +
              glm::vec4(randomDirection() * parameters.alpha, 0.0f));

The color gets the same treatment as the position. Same magnitude $\alpha_d$, same kind of random direction, applied to RGB instead of XYZ:

$$c_{child} = c_{parent} + \alpha_d \cdot u$$

So the color is a correlated random walk over the same fractal structure as the geometry. Particles that are spatially close (same branch, similar depth) inherited most of their color from the same ancestors, so they come out with similar colors. You get smooth gradients running along branches instead of per-particle noise, for free, without writing a single line of color-specific logic.

Set the "Color step" slider in the demo to 0 and everything goes flat grey. Push it up and you can watch the gradients follow the branch structure.

One caveat I did not think about at the time: in the original code the root particle starts at black and packUnorm4x8 clamps to $[0,1]$. With $\alpha_0 = 5$ the very first color step is enormous, so every channel immediately saturates to either 0 or 1 and you get the eight corner colors of the RGB cube and nothing in between. Very vivid, not very subtle. The demo starts from mid grey and scales the color step down, which is why it looks like a gradient instead of a rave.

The bug I found while writing this post

glm::vec3 randomDirection() {
  return glm::normalize(glm::vec3(randf() - 0.5, randf() - 0.5, randf() - 0.5));
}

Looks fine. It is not fine. This samples a point uniformly inside a cube and then normalizes it, which is not the same as sampling uniformly on a sphere. The corner of a cube is $\sqrt{3} \approx 1.73$ times further from the center than the face center is, so there is a lot more cube volume out along the diagonals than along the axes. Normalize all of it and the diagonals get more samples.

Result: branching is denser along the 8 diagonal octant directions than along the axes. There is a faint cubic bias baked into every cloud the generator has ever produced.

Here it is, measured. Both panels use the same set of directions, one plot in space and one as a density-per-angle histogram (red guides mark the diagonals):

The cube version shows clear lobes lining up with the red diagonal guides, and they stay put when you resample or crank the sample count up (which is how you know it is bias and not noise). The Gaussian version is flat and stays flat.

Two ways to fix it:

// Option 1: rejection sampling, keep only points inside the unit ball
glm::vec3 randomDirection() {
  glm::vec3 v;
  do { v = glm::vec3(randf(), randf(), randf()) * 2.0f - 1.0f; }
  while (glm::length2(v) > 1.0f || glm::length2(v) < 1e-8f);
  return glm::normalize(v);
}

// Option 2: Gaussian components, exactly uniform on the sphere, no loop
glm::vec3 randomDirection() {
  return glm::normalize(glm::vec3(gauss(), gauss(), gauss()));
}

Option 2 is the one I would use. A normalized vector of independent Gaussians is uniform on the sphere because the 3D Gaussian is rotationally symmetric, which is one of those facts that feels like a magic trick the first time you see it. Rejection sampling works too and its acceptance rate is $\pi/6 \approx 52\%$, so you throw away about half your samples, which is fine but not free.

Tick "Uniform sphere" in the big demo to compare. Honestly, at high particle counts the difference is subtle. The bias was never going to ruin a puff of smoke. But it is the kind of thing that quietly ruins your results if you ever reuse the sampler for something that actually cares, like hemisphere sampling in a path tracer, where a directional bias shows up as visible energy error.

Where this sits next to noise functions

While writing up the math I kept getting the feeling I had seen this bookkeeping before, and I had. It is octave noise. Fractal Brownian motion, the thing every terrain generator on earth is built on, looks like this:

$$F(x) = \sum_{k=0}^{N} A_0 g^k \cdot noise(x \cdot L_0 \lambda^k)$$

At each octave $k$ the amplitude shrinks by a gain $g$ and the frequency grows by a lacunarity $\lambda$ (usually 2). Compare that to $\alpha_d = \alpha_0 \beta^d$. Same idea: $g$ is exactly my $\beta$, and $\lambda$ plays the part of my childCount.

The difference is where the decay gets spent.

fBm's roughness is controlled by the Hurst exponent $H$, related to gain and lacunarity by $g = \lambda^{-H}$, and the resulting fractal dimension is $D = d + 1 - H$. Compare to my $D = \ln(childCount)/\ln(1/\beta)$. Both are the same shape of statement: one exponent controls fractal dimension via a log ratio of a growth rate to a decay rate.

Why is one additive and the other a tree? fBm has to be evaluatable anywhere and stay continuous, so a query point has to receive a contribution from every octave at once. It cannot look up a point on a random tree unless that exact branch happens to exist. The sunflower generator does the opposite: it commits to a single stochastic path per particle. It is a sample from a branching process, not an evaluation of a field. That is exactly why particles physically cluster along branches, while fBm-driven terrain heights vary smoothly with no branch structure at all.

There is also a nice middle case: midpoint displacement (Diamond-Square, the other terrain classic). Like fBm, its amplitude shrinks by a constant factor per level. Like the sunflower generator, it is recursive and spatial rather than purely additive, since each level subdivides existing points and perturbs new midpoints by a shrinking random offset. It is arguably the closest relative to what I wrote, just constrained to regular grid subdivision instead of free branching in continuous space.

The unifying idea, across all three:

Fractal detail is a geometric series of decaying random perturbations, one per scale level. The dimension is controlled by how much new structure you add per level against how much the perturbation shrinks. The only real difference is whether that decay is applied in physical space (branching walk) or in frequency space (additive noise field).

Which is a fairly grand statement to derive from twenty minutes of writing a test scene, but here we are.

Was it useful?

For its actual job, yes. It fills a particle buffer with something that has structure at multiple scales, it costs nothing (one pass, no neighbour queries, no sorting), it is trivially bounded so you can frame the camera analytically instead of computing a bounding box, and the two parameters give you a spectrum from spiky to fluffy without touching the code.

It also does the one thing you want from a test scene: when the renderer is broken, it looks broken. Static noise hides bugs because noise looks like noise no matter what you do to it. A branching structure has coherence, so if your buffer upload strides are wrong or your point sprites are misaligned, you see it immediately.

The code is in test_gfx_emitter.cpp in Detramotor, hooked up to a few ImGui sliders and a Generate button. About sixty lines including the parameter struct. Sometimes the throwaway code is the code worth writing about.