An intro to Julia sets, and a description of how to render the 4D version of them, follows!
What is a Julia set?
Julia sets consist of a set of points in the complex plane, $z_0$, where the recurrence $z_n+1 = z_n^2 + c$ has a finite limit. That is – you put values into that formula, and some will go to infinity, while some will not. The ones that don’t are part of a Julia set, namely the one for that particular value of c.
The colorful representations you see online are formed by taking $z_0$ as pixel on the screen (where x is the real component and y the imaginary), then repeatedly applying that formula. The color represents the number of iterations used to guess if $z_0$ is in the set. Black is when we’ve reached the max number of iterations.
Quaternion Julia Sets
Quaternions are generalizations of the complex numbers, and so the Julia set notion can also be applied to them. We can set $z% to $a + bi + c j + dk\$ for some a, b, c, and d, just like we chose x and y in the 2D case. To render this, we can take a 3D slice by fixing one of the four dimensions.
To brute force the calculation of this set like we did in the 3D case is more difficult: we now have to consider $n^4$ vs $n^2$ points. And to ray trace the scene, you still have to advance one step at a time, evaluating the recurrence, until you hit a surface that’s part of the set. There’s no obvious way to get the surface normal for shading, either, as you’ve only learned whether the point is in the set or not.
There is a formula that will give you the closest point in the Julia set, though. And you can use this process to safely take larger steps while raytracing – if the closest volume is x distance away, you can advance x along the ray without hitting anything. Most of the time you will not hit the volume by doing that, but you avoid having to evaluate any steps in between.
This formula can also yield an estimate of the normal, using finite distance approximation. We know how to sample the field of distances from the Julia set, and getting the gradient of that gives us the normal. It points in the direction the distance field is increasing the most, i.e. away from the Julia set.
Raytracing Julia sets is simplier than raytracing normal geometries because you don’t need to worry about storing the geometries in memory. Instead, for each pixel, you just go through the unbounding volumes process described above. There is the matter of certain pixels taking longer (because of more unbounding volumes steps), and also of branching.
While I’m not familiar with GPU architectures at the time this code was written (2005), the branching penalty could potentially be very severe.