Note: The full source-code for this project is available here [source]. See Tiny Engine example 4.

To visualize the generated results for many of my projects, I usually write an OpenGL or SDL2 wrapper from scratch. I recently realized that most of these rendering frameworks have a very similar structure, so I decided to unify them.

To do this, I wrote a “Tiny Engine”, with the goal of having an intuitive wrapper for boiler-plate OpenGL3 while keeping it as small and general as possible. The details are on Github.

To test the engine’s GPU rendering, I wrote two small programs: Cellular-Automata in shaders and Multi-Julia Sets.

The programs are very small because all the boiler-plate OpenGL has been wrapped away, and they work quite well. I was particularly happy with the visuals of the Multi-Julia Set program, so I would like to present that here.

In the following, I will show how I produced some beautiful fractal animations in real-time. I do not claim to have an intuitive understanding of the mathematics behind these images (“complex analysis”), so I will link some resources and only briefly explain the math that I actually used.

Multi-Julia Sets

As I understand it, the Julia-Set of a function f is the set of all values z on the stability boundary of the iteration z = f(z).

A typical choice for this function f is a complex polynomial:

   

where i is the iteration index, n is an integer power (“order”), c is a constant “bias” and both c and z are complex numbers.

“Stability” in this context means that as we iterate (i tends towards infinity), the value of z does not “blow up” for a specific choice of initial value and bias.

Note: The “multi” part comes from the fact that we are using an order n larger than n = 2. This is just terminology I guess. See Wikipedia for more info.

Complex numbers can be described on a basic level as two-dimensional numbers with specific multiplication rules.

If we plot the stability behavior (stable or unstable) for every initial value of z (given by a pixel on screen), the boundary between the two regions that form is our Julia-Set!

Note: Numberphile on Youtube has a much better explanation than I can offer (this video inspired me to try this visualization).

The exact shape of the boundary is then dependent on the choice of bias c and the order n of our complex polynomial. This is the foundation of the animated visualizations.

Finite Precision Stability Computation

The easiest way to determine a point z‘s stability behavior is to iterate for a number of steps and determine whether the magnitude of z has crossed some threshold:

//Stability of a complex number z

bool stable(vec2 z){
  int iter = 0;
  while(length(p) < thresh && iter < maxiter){
      p = cpow(p, order)+bias; //Complex Exponentiation
      iter++;
  }
  if(iter == maxiter) return true;
  else return false;
}

Instead of determining stability as a binary choice, we can return the fraction of the maximum iterations after which the stability threshold was crossed. This gives us a “stability gradient” outside the stable region for a nicer visualization:

//Stability of a complex number z

float stable(vec2 z){
  int iter = 0;
  while(length(p) < thresh && iter < maxiter){
      p = cpow(p, order)+bias; //Complex Exponentiation
      iter++;
  }
  //"Stability Fraction"
  return float(iter)/float(maxiter);
}

A number of things can be noted about computing the stability boundary this way:

• The resolution of the stability boundary (i.e. the fractal) scales with the maximum number of iterations
• All values of z with a magnitude greater than 2 are unconditionally unstable (see Numberphile video)
• The order n determines the symmetry of the generated fractal (I can’t easily explain why, sorry)

The finite precision comes from both the limited number of iterations and the choice of the stability threshold. Overall, this is what introduces interesting visual artifacts.

Additional Implementation Notes

The Julia-Set is well suited to computation on the GPU, as stability is determined individually for every pixel. This program was implemented as a fragment shader, enabling real-time animation at high-resolution.

GLSL does not have complex multiplication or exponentiation, but we can implement this easily:

//Complex Product
vec2 cprod(vec2 a, vec2 b){
  return vec2(a.x*b.x - a.y*b.y, a.x*b.y + a.y*b.x);
}

//Complex Exponent
vec2 cpow(vec2 a, int n){
  vec2 b = a;
  if(n == 0) return vec2(1, 0);
  for(int i = 1; i < n; i++)
    b = cprod(b, a);
  return b;
}

The color in the fragment shader is chosen based on the stability fraction. The positions of the pixels are scaled and zoomed according to where one wants to view the fractal.

All of this was wrapped in a Tiny Engine program (see Tiny Engine example 4) to generate the visualizations. One can zoom by scrolling and traverse the Julia-Set using WASD.

Results and Animations

The effects introduced by “numerical artifacts” (i.e. finite precision) are the most interesting for me. Low threshold values (~1.0) give a nice “blob” artifact (don’t know why tho).

Choosing a low maximum iteration number gives a very beautiful continuous color gradient between the stable and unstable regions.

In combination, these two effects create beautiful abstract images, which are even better when animated. Here is a set of animations that I made quickly. I spent most of the time while making these picking a nice color palette.

I added a “random” button to the interface which generates a parameter set in a boundary which I found to work well. With a method for generating color palettes, the animations could be fully randomly generated. I have no experience with generated color palettes though, so it is up to you to set the colors if the animation looks good.

Generating these images is very fun. If you want to do this too, feel free to download the code and give it a try.