Introduction

Quaternion Julia fractals are created by the same principle as the more traditional Julia set except that it uses 4 dimensional complex numbers instead of 2 dimensional complex numbers. A 2D complex number is written as z = r + a i where i² = -1. A quaternion has two more complex components and might be written as q = r + a i + b j + c k where r, a, b, and c are real numbers. There are slightly more complicated relationships between i, j, and k

To generate a quaternion fractal a function is iterated z_n+1 = f(z_n) and if it tends to infinity then it is outside the Julia set, if it is bounded then it is inside the set. Non linear functions are by far the more interesting, in what follows we use the simplest non linear function z_n+1 = z_n² + c where c and z are both quaternions. z₀ is the point in quaternion space being considered and c is a constant that identifies the particular quaternion. Another function that is often used is a cubic, namely z_n+1 = z_n³ + c.

Just as with 2D traditional Julia fractals which are either connected or not depending on the constant c, the same applies to quaternion Julia sets. Some examples of the range of Julia set shapes is given in the column on the right. The disconnected sets aren't represented because they are hard to visualise in 2D, they really require a stereoscopic 3D display and/or projection system. The images created here were rendered using custom software distributed on N separate machines. The goal was to create interactive quaternion sets, at the time of writing 16 machines could render a typical set at about 1 frame every 2 seconds.

The fundamental issue is whether the series "escapes" or not. This cannot always be determined without generating the series to an infinite length, something undesirable in practice. There are two criteria that if they are met allow one to make a decision in a finite time. The first is to assume that if the modulus (length) of a term in the series is larger than some value then the series will escape. This is typically taken to be 4. Another criteria is to limit the number of iterations and if the series hasn't escaped by then, then the point is considered to be inside the set. This maximum number of iterations is often as low as 50, the higher it is the longer the fractal takes to create but the more accurate it is.

Another approximation is the precision at which the 4D quaternion space is sampled. One approach is to relate the pixel resolution of the image to quaternion space and sample one or two times per pixel. Too coarse a sampling risks missing parts of the fractal, a very fine sampling results in a very long compute time. Since these are fractal objects there are good reasons to sample them very finely, this raises interesting questions when one tries to antialias rendering of quaternions by supersampling. Because these are fractal in nature, instead of supersampling generating less noisy surface it can generate more noisy surfaces since there is variation at infinitely small scales.

Since it is rather hard to draw 4 dimensional objects one needs a way of rendering 4D quaternion fractals on a 2D screen. The approach used here is to intersect the 4D solid with a plane, in essence this makes one of the quaternion components dependent on the other three. To get a feel for the true nature of the quaternion fractal one needs create a whole series of slices along an axis perpendicular to the slice plane. This is the same as what one does when drawing contour lines to visualise a landscape, each contour represents a slice of the landscape by a plane perpendicular to the vertical axis. By "stacking" the contours together we gain an appreciation of the surface. Unfortunately in the case of a 4D object we need to "stack" 3D solid objects along a 4th axis which is a little more difficult for our limited 3D visual system.
The following shows 6 slices moving the cutting plane along one of the quaternion axes from the origin in steps of 0.1, c = (-0.08,0.0,-0.8,-0.03).

Smoothly changing the slice plane results in very attractive animation sequences.

If you wish to experiment with quaternion Julia fractals without writing your own code then there is perhaps no better way than by using PovRay, a raytracing rendering program that contains quaternion Julia fractals as a primitive. As a starting point this PovRay scene file may be helpful.

c = (-1,0.2,0,0)	c = (-0.291,-0.399,0.339,0.437)
c = (-0.2,0.4,-0.4,-0.4)	c = (-0.213,-0.0410,-0.563,-0.560)
c = (-0.2,0.6,0.2,0.2)	c = (-0.162,0.163,0.560,-0.599)
c = (-0.2,0.8,0,0)	c = (-0.445,0.339,-0.0889,-0.562)
c = (0.185,0.478,0.125,-0.392)	c = (-0.450,-0.447,0.181,0.306)
c = (-0.218,-0.113,-0.181,-0.496)	c = (-0.137,-0.630,-0.475,-0.046)
c = (-0.125,-0.256,0.847,0.0895)

Most students of algebra are familiar with complex numbers, that is, numbers that have two components called real and imaginary that can often be interpreted as 2 dimensions. A complex number is normally written as a + ib where i² = -1 and a and b are two real values quantities. This idea can be extended to higher dimensions but it turns out that 4 components have useful properties. These are called quaternions and are attributed to Sir William Rowan Hamilton who published a major analysis in 1844 called "On a Species of Imaginary Quantities Connected with a Theory of Quaternions" in the Proceedings of the Royal Irish Academ (2, pp 424-434)

In this discussion we will write a quaternion "Q" as

When performing operations on complex numbers whenever one encounters i² then one knows that is equal to the simpler -1. There are similar but slightly more complicated relationships between i,j,k in quaternion space. They are as follows:

Note that the order in which i,j,k appears in an expression is important. Also note that there is no linear relationship between i,j,k.

Addition (or subtraction) of two quaternions Q₁ = r₁ + a₁ i + b₁ j + c₁ k and Q₂ = r₂ + a₂ i + b₂ j + c₂ k is performed as follows.

The conjugate of Q = Q^* = r - a i - b j - c k.

Multiplication of two quaternions is somewhat involved but follows directly from the relationships above.

Note that quaternion multiplication is not commutative, that is, Q₁ Q₂ is NOT the same as Q₂ Q₁

The length (magnitude) of a quaternion is the familiar coordinate length in 4 dimensional space.

where Q* is the conjugate (see later) which expands to

and

Inverse

The inverse of a quaternion Q^-1 such that Q Q^-1 = 1 is given by

The inverse of a normalised quaternion is simply the conjugate, otherwise the magnitude of the inverse is 1/|Q|. So the above expression normalises the quaternion and then scales by 1/|Q|.

Division of Q₁ by Q₂ is as follows

If m = sqrt(a² + b² + c²) and v is the unit vector (a,b,c) / m then the exponential of the quaternion Q is

The equivalent to polar coordinates in quaternion space are

r = |Q| cos(theta1)
a = |Q| sin(theta1) cos(theta2)
b = |Q| sin(theta1) sin(theta2) cos(theta3)
c = |Q| sin(theta1) sin(theta2) sin(theta3)

theta1 is known as the amplitude of the quaternion, theta2 and theta3 are the latitude (or co-latitude) and longitude respectively. The representative point of a quaternion is the normalised vector (a,b,c), that is, where (a,b,c) intersects the unit sphere centered at the origin.

To rotate a 3D vector "p" by angle theta about a (unit) axis "r" one forms the quaternion

and the rotation quaternion

The rotated vector is the last three components of the quaternion

It is easy to see that rotation in the opposite direction (-theta) can be achieved by reversing the order of the multiplication.

Note also that the quaternion Q₂ is of unit magnitude, and needs to be in order to be a valid rotation.

Converting a quaternion to a matrix

Given a quaternion rotation the corresponding 3x3 rotation matrix M is given by

M =

1 - 2 b2 - 2 c2 2 a b - 2 r c 2 a c + 2 r b 2 a b + 2 r c 1 - 2 a2 - 2 c2 2 b c - 2 r c 2 a c - 2 r b 2 b c + 2 r c 1 - 2 a2 - 2 b2