Originally I'd planned to just write a short guide to generating a fractal using complex numbers, but complex numbers are a huge field with many applications in computer graphics. For now I'll keep everything in this article, but at some point it'll probably be split into different articles.
Introduction to Complex Numbers
The best way to describe a complex number is to say that it's a 2 dimensional number. It has two components, a real value and an 'imaginary' value. It's not really important to know the differences between the two components, except to know that all our standard operations (like addition, multiplication...) might work a little differently when we're working with 2 dimensional numbers.
A complex number can be represented in two ways, the first way is called the Cartesian form and looks like a + bi, where a is the real part and b is the imaginary part (the i isn't really a variable, it's just there to say b is the imaginary component). This is very handy for us working with a 2D display like Milkdrop, because it means we can think of a point on the display, P=(x,y), as a complex number, x is the real part and y is the imaginary part. In this form point P becomes a two dimensional number, and we can use it in an equation as if it were a single number (as long as we remember that the operations work differently), like this:
F(P) = P^2 + c
which is the equation for the Julia fractal. c is a second complex number, which for us is just another float2 with arbitrary values. We'll often assign some oscillating values to c in the form of q variables that we define in the per-frame code, with equations similar to 0.5*sin(time), because this will cause the fractal to change over time. We'll come back to the Julia fractal a bit later.
So the important thing to remember is that for us a complex number is nothing more than a float2, with one variable of any equation we use being the current uv coordinates (uv.x,uv.y).
Operations with Complex Numbers
I'll give an example of each operation and how it might look in Milkdrop. The main thing to note with all these rules is that the form a + bi doesn't mean a plus bi, the + in there is just there to signify that this is a two dimensional number. So when we map this to a float2 we simply ignore the + sign, a becomes x and b becomes y.
For these examples we're going to use these two variables:
float2 uv1 = float2(uv.x, uv.y); float2 c = float2(q1, q2);
Note uv1 is just our normal uv coordinates with a different name, and c is just two arbitrary values that we're importing from the per-frame code. All of the following operations will create a third float2 that would normally get assigned to a variable.
Addition: (a + bi) + (c + di) = (a + c) + (b + d)i
uv1 + c => float2(uv1.x + c.x, uv1.y + c.y)
Subtraction: (a + bi) - (c + di) = (a - c) + (b - d)i
uv1 - c => float2(uv1.x - c.x, uv1.y - c.y)
So addition and subtraction is pretty straight forward, it works just like adding two float2's together because that's all we're doing!
Multiplication: (a + bi) * (c + di) = (ac - bd) + (bc + ad)i
uv1 * c => float2(uv1.x*c.x - uv1.y*c.y, uv1.y*c.x + uv1.x*c.y)
We'll see this one pop up quite a bit, because rotations on complex numbers are all done with multiplication.
Division: (a + bi) / (c + di) = ( ac + bd / c^2 + d^2 ) + ( bc - ad / c^2 + d^2 )i
uv1 / c => float2( (uv1.x*c.x + uv1.y*c.y)/(c.x*c.x + c.y*c.y), (uv1.y*c.x + uv1.x*c.y)/(c.x*c.x + c.y*c.y) )
Division is a little messy, but on the upside we won't use it very often.
Exponents are tricky with complex numbers in Cartesian form. To avoid a complete shit storm of math we're limited to integer components, and realistically we're limited to the exponents 2 or 3. We achieve this simply by treating it as a multiplication that gets repeated for each exponent, uv1^3 = uv1*uv1*uv1. Each time we multiply we need to use the multiplication rule from above and feed the answer into another iteration of the multiplication rule. This quickly gets out of control, and there's no easy way to calculate something like uv1^1.5. For that we'll need to express complex numbers in polar form, but before we get to that lets take a look at how the Julia fractal can be computed in Milkdrop.
The Julia Fractal
The equation for the Julia fractal is z' = z^2 + c, where z is our uv coordinates and c is some arbitrary numbers that control how the resulting fractal looks. Using the operator rules from above, it's pretty easy to work this out:
float1 zoom = 1.6; float2 uv1 = (uv-0.5)*zoom; float2 c = float2(0.2, 0.45); uv1 = float2(uv1.x*uv1.x - uv1.y*uv1.y, 2*uv1.x*uv1.y)+c; ret = tex2D(sampler_fc_main, uv1) + 0.01;
First off, there's a few things we need to say about our 3 variables,
- zoom: It's hard to describe what this controls, but basically it decides how 'open' or 'closed' our fractal will look. Smaller values will zoom in and tend to 'open' the fractal up, while larger values will zoom out and tend to 'close' the fractal in on itself. It can be very picky and you'll probably need to experiment with each new fractal equation to find a value that works. 1.6 happens to work well for the Julia fractal.
- uv1: (uv-0.5) is basically making the origin (point (0,0)) at the center of the screen as opposed to the top left corner. After this the screen can be thought of as a standard graph with x and y axis, both of which run from -0.5..0.5, just like the standard graphs from trigonometry class. Multiplying this by zoom adjusts the size of the graph in the x and y dimensions, since here zoom is 1.6 our graph now runs from -0.8..0.8. This isn't really an important concept, it's just important to know that you always have to adjust uv in this way for a fractal to work.
- c: c is another picky parameter, and the range of 'acceptable' values will be different for every fractal equation.
So what we're doing in the calculation is simple, just multiplying uv1 by itself to get z^2. The y parameter of the new float2 doesn't look like the multiplication rule above, but that's just because we've simplified it from (uv1.x*uv1.y + uv1.y*uv1.x) to (2*uv1.x*uv1.y). We can do that so long as we're multiplying a complex number by itself. Then at the end we just add c, as per the Julia equation.
All the ret line is doing is sampling sampler_main and adding 0.01. This is just the standard way of displaying the fractal if we want to use sampler_main as the source, since it makes each iteration of the fractal equation a little brighter than the last, which gives the fractal its layered look.
You can plug that simple bit of code into the warp shader and you'll get the Julia fractal. Like any fractal, we could loop the equation to be done more than once per pixel, which would give us a very different (often more complex looking) version of the same fractal. It's possible though that multiple iterations will throw off the c and zoom variables so that you'll need to tweak them to get the fractal to display nicely again.
There are literally hundreds of fractal equations out there and not all of them are as easy to write as the Julia equation. For instance, say we wanted to try z' = z^1.5 + c, how do we translate an exponent of 1.5 into our system of multiplying complex numbers by themselves? We can't, at least I don't know how. What we need is a different way of expressing complex numbers, and luckily one exists. It's called Polar Form, and while it looks complicated in the beginning it's actually quite a bit easier to work with. But there's a trade off, working with complex numbers in polar form usually means a lot more instruction slots, because we'll be working with sin/cos/atan functions. The 64 instruction slots in PS 2.0 just won't be enough, so polar form is pretty much limited to PS 2.X or above and will still take a toll on older video cards.
Since this is meant as an introduction to complex numbers we're not going to bother figuring out why polar form works, just how. The basic idea is simple, we turn z = (a + bi) into z = r(cos(θ) + sin(θ)i). In this form, r is called the modulus (or absolute) form of z, which is calculated with the equation: |z| = sqrt(x^2 + y^2) Since for us z is just the uv coordinates, in shader code the equation looks like this:
float1 modulus = sqrt(uv.x*uv.x + uv.y*uv.y);
You'll notice this is just the Pythagorean theorem of a right triangle. Our x and y values are the two right sides of the triangle, and the absolute form of z is just the hypotenuse of this triangle, which is the distance from the origin (the center of the screen) to the point z. The θ (theta) is called the argument (or angle) of z. To find it we need to think of our complex number z as being on a graph:
A complex number "A" on a graph. Trig 101 tells us that to find the angle θ we need to take the arc tangent of y/x, and the shader has a built in function for this which looks like:
float1 theta = atan2(uv.y, uv.x);
Now we have all the variables we need to put our complex number into polar form:
float2 uv1 = float2(modulus * cos(theta), modulus * sin(theta));
With this we can use something called De Moivre's Theorem to raise z to any power we want without using the multiplication rule. De Moivre's equation is almost the same as the normal polar form equation:
File:Dem5.gif To see it in action, lets do the equation we were having trouble with earlier, z' = z^1.5 + c:
uv1 = float2(pow(modulus, 1.5)*cos(theta*1.5), pow(modulus, 1.5)*sin(theta*1.5))+c;
And the Julia set equation would look like this:
uv1 = float2(pow(modulus, 2)*cos(theta*2), pow(modulus, 2)*sin(theta*2))+c;