Hyperdimensional Quadratic Julia Set Fractals
Jeff Kelly
CSE 168 Spring 2005
Final Rendering: 183 seconds at 1024x768 on A64 3200 
Introduction
A Julia set is defined roughly as the union of all points in an iterated complex function R(z) that don't diverge as the number of iterations of the function approaches infinity. This definition encompasses a huge number of functions, and analysis is generally restricted to a subset of complex functions. Here, we consider functions of the form:
z(n) = z(n1)^2 + c
where z(n) is the value of z after the nth iteration of the function R and c is a complex constant. The constant c defines the shape of the fractal, and highly different shapes can be observed by tweaking c. For each point in the complex plane, this function is iterated and tested for divergence. [1]
Extension to Higher Dimensions
Normal complex numbers of the form a+bi only allow for twodimensional visualization of a fractal. To render in three dimensions, we must use a different representation of points in complex space. Two common choices exist: quaternions and hypercomplex numbers. Hypercomplex numbers are a generalization of normal complex numbers, whereas quaternions are complex numbers designed for representing 3D rotations. Neither choice is perfect, as multiplication is noncommutative in quaternion algebra, whereas hypercomplex numbers may not have a multiplicative inverse. I chose to use quaternions as I am more familiar with them, and the operations necessary for evaluating a quadratic Julia set are very simple in quaternion algebra [2]:
z = (r, x*i, y*j, z*k)
z^2 = (r*rx*xy*yz*z, r*r*x, r*r*y, r*r*z)
z+c = (z.r+c.r, z.x+c.x, z.y+c.y, z.z+c.z)
Practical Considerations
Obviously, we can't iterate each point to infinity. Instead, we set a maximum number of iterations and see if a point diverges before the maximum number of iterations is reached. Additionally, a useful property of Julia sets is they are bounded by a sphere of radius 2 centered at the origin. Hence if an iterated point moves more than 2 units away from the origin, it is safe to reject it. In practice, fairly small numbers of iterations give nice results. Higher numbers of iterations result in a more accurate fractal surface, but the more accurate surface is generally less pleasing visually due to the smallscale turbulence.
10 Iterations  15 Iterations  20 Iterations 
Additionally, because we are moving from 4D complex space (represented by a quaternion) to 3D visual space, we must somehow sample the 4D space to get a 3D answer. Although this is difficult to visualize, we are essentially freezing a constantly evolving shape at a point in time. This is done by setting one of the dimensions of the quaternions to a constant value.
Animation showing time slices from 1.0 to 1.0 by 0.1 
Raytracing
A straightforward method of generating the points in the set is to discretize the reduced 4D space and test each point. Although it is simple, this method has the huge disadvantage of being incredibly slow for any reasonable granularity, as well as only being as accurate as the size of the grid cells.
Although this method was implemented, it was found to be too slow and too limited for interesting effects that can easily be achieved with raytracing. Instead, I and implemented a distance estimator for the Julia set [3,4,5]. The distance estimator gives a lowerbound on the distance to the fractal from any given point outside the fractal, and is represented as follows:
distance >= alpha * ( norm(z(n)) / norm(z'(n)) )
norm(z'(n)) = 2*norm(z(n1))*norm(z'(n1))
The alpha parameter is analogous to the grid size, and small values (< 0.01) are necessary for correct results. Higher values can result in rather stylistic fractals:


Alpha = 0.09  Alpha = 0.0015 
Once we have a lowerbound on the distance, we step along our ray until we hit the fractal or have stepped too many times without hitting the fractal. This is analogous to the idea of an "unbounding sphere."
Figure from [3] 
As a minor optimization, rays are first intersected with a sphere of radius 2 around the origin, then only tested against the fractal surface if they intersect this sphere. This helps speed up tracing of rays that wouldn't intersect the fractal anyway.
An additional challenge of raytracing these fractals is that lighting calculations require surface normals, yet fractals have infinitely perturbed surfaces (and hence, no normal). Hence, the best we can do is a rough approximation of the general direction of the surface. To do this, we sample the surface at surrounding points and take the crossproduct of the resulting vectors, which gives us a rough idea of how the surface is changing at a point. [4]
Some random examples 
Results
References
[1] Eric W. Weisstein. "Julia Set." From MathWorldA Wolfram Web Resource. http://mathworld.wolfram.com/JuliaSet.html
[2] Eric W. Weisstein. "Quaternion." From MathWorldA Wolfram Web Resource. http://mathworld.wolfram.com/Quaternion.html
[3] Hart, J., Sandin, D. & Kauffman, L. "Ray Tracing Deterministic 3D Fractals," SIGGRAPH, 1989, pp. 289296.
[4] Dang, Y., Kauffman, L.H. & Sandin, D. Hypercomplex Iterations, Distance Estimation and Higher Dimensional Fractals Chapters 3, 7, 9
[5] Peitgen, H.O., Saupe, D. The Science of Fractal Images. SpringerVerlag: New York, 1988.