Iterative Fractals on the GPU

CS 493 Lecture, Dr. Lawlor

It's often the case that we want to use a small simple shader to simulate unlimited detail. One way to get this is via "fractals", self-similar geometric objects. For example, this PixAnvil program renders the Mandelbrot Set:
GPU Fractal Rendering demo

The basic rationale is to pick a constant c, and see how long this iteration lasts before exploding to infinity:
z = z² + c
This turns out to be really boring for real numbers z and c, so the tradition is to pick them from the complex plane. This is easy in 2D, where you just make the x axis the real axis, and y axis the imaginary axis. (You can also pick z and c to be 4D quaternions, but this gives you Mandelbrots in the axis aligned cross sections through the origin, and not too much in between.)

Here's the iterative code, run in a pixel shader:

/* We start with two complex numbers:
    c is a constant, giving your location onscreen.
    z starts at c, and we'll iterate it.
*/
float cr=-0.123, ci=0.745;
float zr=texcoords.x, zi=texcoords.y;
float mag=0.0;

/* Mandelbrot Set fractal iteration: */
for (int a=0;a<50;a++) // <- higher constant gives more detail, but is slower.
{
    // complex magnitude
    mag=zr*zr+zi*zi; 
    if (mag>4.0) break;

    // complex multiply: n=z*z
    float nr = zr*zr - zi*zi; 
    float ni = 2.0*zr*zi;

    // z=n+c
    zr = nr + cr; 
    zi = ni + ci;
}

// Colorize according to the final z 
return vec3(
   /* red: */   fract(0.2*time+zr*0.5),
   /* green: */ fract(0.13*time+zi*0.5),
   /* blue:  */ mag/16.0
);

That final colorization is my little addition. The time-dependent fract is what causes the color animation. It's more traditional to use the final value of "a" to decide the color, but this is difficult in OpenGL, and to my eye doesn't look as interesting.

Modifying the initial condition setup to fix c and vary z onscreen results in the closely related "Julia Set", with this particular c value giving Douady's Rabbit (also known as the Dragon).

float cr=-0.123, ci=0.745;
float zr=texcoords.x, zi=texcoords.y;

Actually, modifying nearly *anything* results in a new fractal, many of them quite interesting looking! For example, taking the absolute value of zr and zi before squaring results in the burning ship fractal. Modifying the constant 2.0 as a function of time also creates interesting behavior.

Note that the basic model here is a combination of (1) a nonlinear operation, and (2) repetition. These two almost always combine to give interesting behavior. Here are some examples:

Fract + rotate

float angle=time*0.08;
float g=1.2;
float s=g*sin(angle), c=g*cos(angle);

/* fractal iteration: */
for (int a=0;a<5;a++)
{
    zr=fract(zr);
    zi=fract(zi);

    // coordinate system rotation
    float nr = c*zr + s*zi; 
    float ni =-s*zr + c*zi;

    // z=n+c
    zr = nr + cr; 
    zi = ni + ci;
}

$Fract + rotate fractal image$

Multi-octave version of above

vec3 ret=vec3(0.0);

/* Multi-octave fractal iteration: */
for (int a=0;a<10;a++)
{
    zr=fract(zr);
    zi=fract(zi);

    // coordinate system rotation
    float nr = c*zr + s*zi; 
    float ni =-s*zr + c*zi;

    // z=n+c
    zr = nr + cr; 
    zi = ni + ci;


// Colorize according to each z after the first few
   if (a>3)
   ret+=(1.0/7.0) * vec3(
   /* red: */   fract(0.2*time+zr*0.5),
   /* green: */ fract(0.13*time+zi*0.5),
   /* blue:  */ mag/16.0
   );
}

return ret;

The nonlinear operation here is a multiply:

for (int a=0;a<10;a++)
{
    // Limit range of z
    zr=fract(zr); zi=fract(zi);

    // Nonlinear smooshing
    float nr = (0.5+fract(time*0.1))*zi*zr;
    float ni = zr;

    // z=n+c
    zr = nr + cr; 
    zi = ni + ci;
}

I've modified the usual keyboard user interface so the WASD pans you around in the plane, and Z and Q zoom in and out. Eventually, if you keep zooming in, you'll run out of numerical precision. This happens pretty quickly with the usual single-precision floats of the graphics card.