Fun with Fractals
This week I decided to go over a short research essay I did for ‘Numerical Methods & Advanced Mathematical Modeling’ at Trinity College.
In it I go over fractals, their origin, use and so on, but I’ll spare you the boring details and go directly to the fun part, sticking it on the GPU (and on the 360), namely the Mandelbrot set.
Some history…
For those who don’t know, a fractal can be roughly described as “a fragmented geometric shape that can be split into parts, each of which is (at least approximately) a reduced-size copy of the whole”. This term was coined by Mandelbrot himself, who used them to study natural shapes, like in his paper “How Long is the Coast of Britain?” tackling the coastline paradox. In 1980, Mandelbrot defined the Mandelbrot set, a set of parameters for which a certain Julia set is connected. The Mandelbrot set is still one of the most popular fractals ever discovered.
The algorithm…
I’d like to start out by naming my sources. One of my first references for implementing the Mandelbrot on the GPU was Shawn Hargreaves (you can check his article here, great stuff). Another reference worth having a look was Linas Vepstas’ article on renormalizing the Mandelbrot escape (see here). Check the paper for further references.
So, how do we get the Mandelbrot on a shader in the first place? First let’s go over ways to approximate the Mandelbrot numerically. The most popular way of estimating and representing the Mandelbrot set is to use the escape time algorithm. The escape time algorithm calculates evaluates the connectivity of each set from the famous quadratic polynomial
on the plotted area by testing whether or not the orbit of its starting point tends towards infinity under iteration. To do this we test the values against an escape condition. If that condition is reached then the tested set is considered disconnected (its orbit tents towards infinity) and it’s disregarded from the set. If, on the other hand, the escape condition is never reached after a maximum number of iterations, it is assumed that its orbit does not tend towards infinity and it is thus included into the set.
There are different test conditions, but the most common condition (for its simplicity) bases itself on the fact that no complex number of ‘z = x + yi’ where ‘x > 2′ or ‘y > 2′ can be part of the set.
So now that we know how to estimate the set, let’s see how we can represent it using the escape time algorithm. One of the most common ways would be to assign a binary set of colours (usually black and white) depending on whether or not the point is considered part of the set. That’s not very exciting though, we want pretty pictures, with all the shiny colours, so instead we use the number of iterations (how fast a given set’s orbit tended towards infinity) as an index for our colouring from a colour map, this map can be defined by an artist to make it all ‘ooooh’.
Taking all the stuff we discussed above we can express an algorithm as such…
Where complex values z and c compose the above mentioned quadratic polynomial, e defines the escape radius and m defines the maximum number of iterations.
Complex numbers 101
That’s all very nice in paper but we can’t just write that down on our shader straight away, so for those without alot of background on complex numbers, let’s break them down to size.
A complex number is defined by two parts, real and imaginary. The usual notation for these numbers is ‘a + bi’ where a and b are real numbers and i represents the imaginary unit. We’re going to plot a certain area of the complex plane with the above algorithm. The complex plane is a two-dimensional Cartesian coordinate system where a complex number represents a point or vector on said plane.
Let’s look at some basic associative, commutative and distributive complex number operations to help us flesh out our final algorithm.
For a complex number z = a + bi and a second complex number
w = c + di we de fine the following operations…
- Addition:
- Multiplication:
from this we can conclude…
- Absolute value:
And this way we can get a shader friendly algorithm like the following…
Pretty colours…
This is great, but we can get even nicer results by getting rid of the step effects we got from using an integer value for mapping our colour values.
To solve this we use the renormalized escape time algorithm defined by Linas Vepstas where the colour index is defined by…
A great explanation of this is given by Earl L. Hinrichs (quoted from Linas Vepstas).
“Ignoring the +c in the usual formula, the orbit point grows by z = z^2. So we have the progression z, z^2, z^4, z^8, z^16, etc. Iteration counting amounts to assigning an integer to these values. Ignoring a multiplier, the log of this sequence is: 1, 2, 4, 8, 16. The log-log is 0, 1, 2, 3, 4 which matches the integer value from the iteration counting. So the appropriate generalization of the discrete iteration counting to a continuous function would be the double log.”.
The code…
Now its just a matter of converting the algorithm into code. One optimization that’s worth noting here is the use of a texture lookup as a colour map. This can make use of a wider range of colours than the ones available in the map itself through the graphics card’s filtering, which can be pretty much free. In my example I only make use of one dimension when applying the colour map, I left it at that for the sake of simplicity. Also, it would be extremely easy for an artist to edit the colour map applied to the fractal.
// Calculate initial point
float2 c = (texCoord - 0.5) * Zoom * float2(1, Aspect) - Pan;
float2 v = c;
// Check for connectivity
// Escape Time Method
// IF absolute of Z > Escape Radius
- Disconnected (Tending towards infinity)
// OR Iterations == Iteration Limit
- Connected (Orbit bounded, assumed to be not tending towards infinity)
int n = 0;
while( v.x*v.x + v.y*v.y <= EscapeRadiusSquared && n < Iterations )
{
// Evaluate function: Fc(x) = z^2 + c
v = float2(v.x * v.x - v.y * v.y, v.x * v.y * 2) + c;
n++;
}
// IF disconnected, calculate color from number of iterations
float4 color = 0;
if(n < Iterations)
{
// Renormalize interation count
// a couple of extra iterations helps
// decrease the size of the error term.
v = float2(v.x * v.x - v.y * v.y, v.x * v.y * 2) + c; n++;
v = float2(v.x * v.x - v.y * v.y, v.x * v.y * 2) + c; n++;
float absolute = sqrt(v.x*v.x + v.y*v.y);
float mu = n - log(log(absolute))/log(2); // double log (0,1,2,...)
// Fetch color from color map (texture)
// Bilinear filtering will increase number of possible colours for free!
color = tex2D(TextureSampler, float2(fmod(mu,NumberOfColors)/NumberOfColors,0));
}
// Final plot result (coloured)
return color;
And there we go, a nice and pretty Mandelbrot shader running on real time on the GPU. The video below gives you a nice look at how it turns out, a beauty.
Hope you enjoyed this little article and I will be back soon with more stuff to show.
See ya!
Loading ...