Ray Tracing Quaternion Julia Sets on the GPU Keenan Crane (http://graphics.cs.uiuc.edu/\~kcrane)
This project takes advantage of the floating point power of recent GPUs to quickly visualize quaternion Julia sets. Pictures above were ray traced in less than a second on a pair of GeForce 7800 GTX graphics cards running in SLI (top picture originally rendered at 1280 x 1024 and downsampled 50% to reduce aliasing).
First, a few "frequently asked questions" about the project:
Q: What in the name of Gaston Maurice Julia are quaternion Julia sets?!
A: A Julia set is a kind of fractal, or object that looks somewhat like itself at many different scales. The quaternion Julia set is a four-dimensional version of this fractal (more detail below). Playing with a small number of parameters for a set results in a huge number of complex, beautiful shapes such as the one shown above.
Q: Why would you want to render these things on the GPU?
A: In general there are two problems with quaternion Julia sets:
- They take forever to visualize
- They're not useful for anything
We can solve problem number one using ray tracing and the GPU: the GPU is really good at doing lots of things in parallel as long as none of those things depend on each other. Ray tracing is a great application of the GPU because each ray can work independently. Problem number two, however, is an open research problem.
Q: But can't you use the quaternion Julia set to render taffy??
A: Yes! In fact, here's a link to a movie comparing a taffy machine on the Santa Cruz boardwalk to a rendering of several Julia sets (taffy footage courtesy of Yisong Yue):
The quaternion Julia set does a screen test (movie)
Alas, Pixar is not expected to cast taffy in a lead role anytime in the near future.
Julia Sets in 2D
Before explaining Julia sets in the quaternions, let's take a look at the traditional two-dimensional Julia sets. These eponymous fractals were the invention of French mathematician Gaston Julia. The fractal exists in the complex plane, a coordinate system where the "x" component of a point's location corresponds to a real number, and the "y" component corresponds to an imaginary number (i.e., a single number x such that x2 is less than zero). Each point in the complex plane is a complex number of the form z = a + bi. A Julia set (technically a filled-in Julia set) is the set of all points z0 in the complex plane for which the sequence
zn+1 = zn2 + c
has a finite limit (i.e., does not get arbitrarily large as n approaches infinity). Geometrically this recurrence pushes points around on the complex plane, since each iteration will produces a new complex number / point on the plane. Points included in a Julia set will hover around the origin, while points not in the set will shoot off to infinity. Different constants c specify a particular Julia set, and can be thought of as placing differently shaped boundaries on the plane which prevent points inside the boundary from escaping.
In practice, to determine whether a given point diverges we start to compute points in the sequence above and see how quickly their magnitude (geometrically: their distance from the origin) gets large. In fact, if the magnitude of any point in the sequence exceeds a value of 2 (geometrically: the point leaves a circular region of radius 2 centered at the origin), the sequence diverges and z0 is not in the set. By applying the convergence test to all the points corresponding to pixels in an image, we get an image like the one below. This particular set is known as the "Rabbit Ears" Julia set (c = -0.12 + 0.75i). Grey values in the image correspond to the number of iterations used to guess if a point is in the set.
"Rabbit Ears" Julia set
Julia Sets in 4D
The quaternion Julia sets are almost exactly like the original Julia sets, except that points in a set are from a four-dimensional complex space called the quaternions. This variation on the fractal was first explored by Alan Norton at the IBM T.J.Watson Research Center . The quaternions are similar to the complex numbers, except that a quaternion has three imaginary components: z = a + bi + cj + dk. Plugging one of these points into the same recurrence will give you similar behavior: some of the points will "escape" (their magnitude will approach infinity), and some will move around inside a region close to the origin. However, visualizing this fractal becomes tougher in 4D for a couple reasons.
First, we can't visualize the set of 4D points directly as we can with the 2D fractal. Instead, we have to come up with a way to "project" the set into a lower dimension. One way to look at a high dimensional object in a lower dimension is to take slices of the object by finding its intersection with several lower-dimensional planes. For example, if we intersect a three-dimensional sphere with a two-dimensional plane, we get a two-dimensional circle whose radius depends on the location of the plane relative to the sphere. By looking at a number of these cross sections we could infer the true shape of the sphere in 3D. Putting an egg through an egg slicer gives a similar result: each blade defines a slicing plane through the egg, and we can infer the shape of the yolk from the shape of yellow regions in each slice. So if we take 3D slices of a 4D Julia set we can get some idea of what it really looks like.
2D slices of a 3D egg
In the case of the Julia set, we specify a 3D slice by picking three basis quaternions. The idea of a basis in a 2D plane is simple: pick two basis vectors which point in different directions in the plane. We can then get to any point in the plane by starting at the origin and moving some distance in the direction of the first basis vector and some other distance in the direction second basis vector. (Think about where we would be able to go if we had only picked one basis vector). If we instead pick two (and only two) 3D vectors, then we still have a basis for a 2D plane, but this is a 2D plane sitting in a 3D space. We can now get to any point in a slice of 3D space by moving different distances along our two 3D basis vectors. If we need to know where a point (x, y) from our 2D plane sits in the 3D space, we simply start at the origin of the 3D space, move a distance x along the first basis vector, and a distance y along the second basis vector. Similarly, if at any time we need to know where a 3D point (x, y, z) is in the 4D quaternion space, we start at the origin and move a distance x along our first basis quaternion, a distance y along our second basis quaternion, and a distance z along our third basis quaternion.
Left: two 2D vectors which define a basis for the 2D plane. Middle: any point in the 2D plane can be thought of as some combination of these two vectors. Right: two 3D vectors which define a basis for a 2D subset of points in 3D space.
The second problem with the 4D version of the Julia set is that our scheme for generating an image of the set is much less efficient in 4D. While we could test the convergence of every point on a 4D grid, the cost would be the same as rendering n stacks of n 2D images, or n2 times more expensive, and would consume n2 times as much memory. Even if we used a grid of only three dimensions and tested the convergence of points in the corresponding slice there would be wasted effort: in the end we only want to look at the surface of the 3D projection, and we don't care about interior points.
Since we are only interested in looking at the surface of a 3D object, ray tracing seems like it's worth investigating: ray tracing is often useful for detecting the thing that's closest to the viewer. However, there is no known way to analytically find the nearest point along a ray which intersects a Julia set (this is the usual method for ray tracing simple objects like spheres or cones). Instead, we could take small steps along the ray, testing the convergence at each point, and stop when we reach the first point which does not diverge. This might work, but we would have to evaluate a huge number of points for every pixel rendered. Additionally, there is no obvious way to shade the surface once we find an intersection, because the fractal surface does not yield a surface normal.
Fortunately, there is a distance estimator which will tell us, given any point z in the quaternion space, the distance to the closest point in the Julia set. This distance estimator can be used to accelerate the ray tracing process using unbounding volumes, a method presented by professor John Hart et al in Ray Tracing Deterministic 3D Fractals . Note that the formula given for the distance estimator in this paper is a slight misprint - it should be:
where z' is just another sequence of points:
z'n+1 = 2 zn zn'
which can be computed alongside zn. The first point in this sequence should always be z0' = 1 + 0i + 0j + 0k. Just as the magnitude or norm of a complex number is its distance from the origin of the complex plane, the magnitude of a quaternion z = a + bi + cj + dk is simply its distance from the origin of quaternion space sqrt(a2 + b2 + c2 + d2).
The idea of unbounding volumes is simple: imagine we're tracing a ray through quaternion space, looking for the first point in the set along this ray. At any time we can query the distance estimator, and it will tell us the distance from our ray's origin to the closest point in the Julia set (in any direction - not necessarily in the direction of the ray). We can then safely take a step of this same distance along the ray direction, because we are guaranteed not to skip over any points in the set (for if we passed a point, that point would have been the closest one, and the distance estimate would have been smaller).
The animation below illustrates several steps of this process. If we pay attention to the first ray tracing step in the animation, we see the benefit of this method: a ray originating far from the Julia set will move much more quickly than if we had taken many small, uniform steps. This is because the initial few unbounding steps account for a large part of the distance between the ray origin and the surface.
Several steps of ray tracing using unbounding volumes
Because we are always taking steps which only cover part of the remaining distance, we can never actually reach the surface (except in the unusual case that the closest point in the set is directly along our ray). This phenomenon is known as Zeno's Paradox of Motion. Zeno was an ancient Greek philosopher whose Paradox of Motion stated something like the following: In order to get from one point to another you must first travel half the total distance. However, to travel the remaining distance, you must first cover half the distance between the midpoint and the destination. Since there will always be another half to cover (he claims), one can never reach the final destination.
In reality we will eventually reach a stopping point because of the limited precision of floating point values. However, getting to that point would probably take a huge number of steps, and we would like to be able to stop before then for the sake of computational efficiency. In practice we specify some small value epsilon and say that if the minimum distance to the set (given by our distance estimator) is less than or equal to epsilon, then we have reached the surface.
In this case we are actually rendering an isosurface of the distance estimator rather than the Julia set itself. An isosurface of a function is simply all the points where that function has the same value. An example of this would be rings on a tree stump: each ring can be thought of as an isosurface of a function mapping cells of the tree to their age in years, since all cells in a ring were grown in the same year.
A tree trunk displays age-isosurfaces of its cells
The isosurface we're rendering will not be as detailed as the Julia set itself -- the distance function looks more like a smoothed version of the set. Detail can be increased by making epsilon smaller, but this also increases rendering time due to the greater number of steps required to reach the isosurface. There is one redeeming property of rendering an isosurface, however: since the isosurface is a continuous function we can generate normals for shading the surface. Surface normals are impossible to define on the Julia set itself because there is detail at every level. More details on generating normals from the isosurface can be found in .
Ray Tracing on the GPU
The GPU was never meant to do strange stuff like rendering fractal sets or tracing rays in four-dimensional space, but in the last few years the GPU has become much more powerful than the CPU, and many people are becoming interested in its use for applications outside of typical video game graphics. This topic is known as GPGPU or General-Purpose computation on Graphics Processing Units, and encompasses anything from protein folding to options pricing.
Most GPGPU applications ignore the majority of graphics card features (vertex processing, fog, lighting, etc.) and instead run a complex fragment program on a single rectangle which covers the entire screen. Each fragment in this rectangle can be thought of as a separate process which works independently of all other processes, and only differs from other processes as a result of its input data.
To ray trace a Julia set, we run a fragment program which steps a single ray through quaternion space using unbounding volumes as described above. In a high-level shader language like Cg, the fragment code looks similar to CPU ray tracing code, and can be examined here: QJuliaFragment.cg. The only information that needs to be sent to the program is the origin (starting point) and direction of the rays which make up the image. Since the rays vary linearly across the screen, we need only to specify this information at vertices and the correct ray values for each fragment will be interpolated across the rectangle. We can use, e.g., the three color components (red, green, blue), to send three spatial coordinates (x, y, z), resulting in input rectangles like the ones below.
Left: Input ray origins. Right: Input ray directions
Although the GPU is good at quickly doing many floating point operations relative to the CPU, it still suffers from coarse control flow across processes (fragments). Because the GPU is built to render simple graphics, it expects that nearby pixels will usually follow the same branches in a program. All pixels within a relatively large block must wait for each other to finish before moving on to the next batch, but this wait time is typically small (or zero) for the simple shaders used in video games. Unfortunately, many GPGPU applications do not behave this way. In the ray tracer, for instance, two pixels in a block may differ greatly in their distance from the viewer, meaning closer rays effectively take the same amount of time to render as distant ones. The figures below illustrate the general difficulties of branching on the GPU: if the black tree on the left represents all possible branches a process can take, and the colored lines indicate the paths actually taken by three different processes from the same block, then the three lines on the right represent the time needed to actually run each process, highlighting the region where useful work is done. Despite all this, a GPU quaternion Julia ray tracer still far outperforms a CPU implementation.
Left: Potential control flow decisions in a fragment program, and paths taken by different processes in a block.
Right: Effective running times of processes from the same block
Q: The pictures are neat, but this quaternion-fractal-GPU stuff is truly esoteric and convoluted. Is there some code I can just play with?
A:Yes! Well-documented source code and executables for a couple GPU implementations of the ray tracer are available for download from here.
Q: 2D, 4D, ...what about 3D? 5D? n-D?
A: A three-dimensional Julia set cannot be defined in the same way as the 2D and 4D versions because multiplication (which is needed to generate a sequence of points) in a three-dimensional complex space is ill-defined. However, the Cayley-Dickson construction gives us a method of generating any 2n-dimensional complex space, and these hypercomplex spaces have a sum, product, and norm defined, which are all we need to evaluate the distance estimator. The first few hypercomplex spaces are the quaternions (4D), the octonions (8D), and the sedenions (16D). Unfortunately, as we continue to increase the dimension of a Julia set, a conservative distance estimator becomes less and less useful: the closest point in the entire 2n-dimensional set probably won't be a good estimate for the closest point in the slice we are rendering.
- Norton, A. 1989. Julia Sets in the Quaternions. In Computers and Graphics 13:2, 1989, 267-278.
- Hart, J. C., Sandin, D. J., and Kauffman, L. H. 1989. Ray tracing deterministic 3D fractals. In Proceedings of the 16th Annual Conference on Computer Graphics and interactive Techniques SIGGRAPH '89. ACM Press, New York, NY, 289-296.
- great introduction to complex numbers and quaternions
Great Images, Great Work, and, definitely, really GREAT documentation!
but the BEST thing is, just the same evening, just about 1 hour before i saw your picture and read your text, i had a talk about something with some friends. and guess what (actually the discussion was about racism, politics, and the generic "how the world sucks today ") happened to get part of the discussion? exactly: the Zeno's Paradox of Motion. We didn't knew it's name. Now we know, hehe.
Sometimes life is strange. Over years you haven't heard about it anymore, and then, twice in about an hour at completely different places..
Oh, and, the way to trace towards the surface is great. Similar implementations exist to render heightmaps (on hardware to do the better bumpmapping/parallax mapping. I think ATI has a paper about this way to accurately trace the bumpy displaced surface, i'll check that out). It's definitely a great technique. if you can accurately measure the closest distance (or estimate a value below or equal lowest distance), you can trace very fast towards the surface, not missing any detail. and it works even for "Teapot in a Stadium" situation.
Perfect for the quaternion julia sets, wich i love. espencially when animated