The mathematical background on julia
The Julia-set(s) is generated from the points (in the complex plane)
generated by the equation:
zn = zn-12+c
where c is any complex number.
The points that dont escape (ie they do not grow out of bounds)
are those that belong to the Julia-set.
The normal coloring-scheme is to set the color according to how
many iterations it takes before the we reach a point we now will
escape.
So you want to generate the julia-set.
First try:
Some investigation shows that the rect ((-2,-2),(2,2)) is a
good one for the julia.
Choose a complex number c0 (this is the value that specify this julia-set)
For each point in the rect get its logical cordinates
then run this point through iterate and see how many
loops we get (the number of loops is the color-value).
So what does iterate do?
Here we must introduce a value: MAX-ITERATIONS,
if we loop through iterate MAX-ITERATIONS times we say
that the point is in the julia-set. A value of 80+
seems to be good enought for a reasonably fast generation. A
Higher value gives a more precise Julia-set, but takes more time
to generate.
function
iterate (z,iterations) {
if (iterations > MAX-ITERATIONS) // to many iterations ?
return (iterations - 1) // we add one when we get back..
z = z**2 + C0 // square z and add c0
if (abs(z) > SOME_VALUE) // growing out of bounds?
return 1 // yes it does.
return (1 + iterate (z, iterations+1)) // nop, try again
}
note abs(z) = squareroot(z.re2+z.im2)
SOME_VALUE can be set to 8 for starters.
Conclusions
This works but is painfully slow.
What can we do better? After some investigation we see that
what iterate really does is this: It jumps from one point in the
plane to another. How would it bee if we already have a value for
that point? If we had a precalculated value we could then return
that value + 1 to the point we just visited and save this value.
"But I work with float for precision"
I hear you say, and its true that If we choose to work with
floats it is not possible to tabulate all values.
Go with Integers if you want speed I say.
2:nd try
to work with integers we need to do some work, we cant simply use
the integers from our rect. If we did we would get a 5 * 5 pixel julia
(That would indeed be fast to generate but not so big as we want).
If we want julia of size 256 * 256 pixels we see that 1 in rect
is 64 pixels on the screen (256/(2-(-2)) = 64).
When we take the square of z we now get a value that is to big, no
problem we divide (or shift actually) it to the correct size again.
Here we will get a roundoff error but the speed gain
we will get compromises this very well (I think).
Using only integers we now tabulate all values we reach (within rect).
Iterate now looks something like this:
function
iterate (z, iterations) {
if (tabell(z.re, z.im) != 0) { // do we have a previous value
if (iterations < MAX-ITERATIONS) { // to many iterations
z = z**2 + c0
if (abs(z) > SOME_VALUE) // growing out of bounds
it = 1 // yes it does.
else
it = (1+iterate (z, iterations+1)) // nop, try again
}
tabell(z.re, z.im) = it
}
if (it == MAX-ITERATIONS)
it = it - 1 // previous iterate add one to this
return it;
}
conclusions
This is much faster. A shame its a littel jagged in the edge (try setting
colors as dark,light,dark,light... to see, or do as I do set them
near each other so one doesnt notice.
Further into wizardry
Have you noticed that the julia-set is mirrored around the axis?
so when you calculate one value you are actually calculating two.
This means that we only need to calculate half of the values (quite
a great improvement).
Precalculate the outher bound of the julia-set. This is nothing strange
as you can see we take the absolute value of a imaginary number and compare it
to another. This is actually to check if it is within a circle (in the
complex plane) with radius SOME_VALUE.
By doing this we save ourself some work all those points outside this circle
is spinning off to infinity fast so we dont need to calculate their values.
One other improvement we can use if we have this circle calculated is that
we dont have to calculate the whole of z2. Start by calculating
the new imaginary part and then only if zn+1.im is within range
do we calculate the real part (I do it this way because I store the circle
as a colum of startvalues, indexed by the imaginary part).
in c iterate now looks like this:
int
iterate(int x, int y, int iter){
int it; // return value
int nx, ny; // z(n+1)
int v = x - (y * 256); // for fast index.
if ((it = pix[v]) == 0){ // prev value
if (iter < iterations) { // nop :(
ny = (x * y + im) >> 5; // imaginary part
//(2*re*im/64)=(re*im)/32
it = 1;
if (ny < 128 && ny > -128) // whitin range?
{
nx = (sqr[x] - sqr[y] + re) >> 6; // tabulated squares.
if (nx <= hww[ny] && nx >= -hww[ny]) // whitin circle
it += iterate (nx, ny, iter+1) ; // yes, so calc new value
}
} else {
it = iterations;
}
pix[v] = pix[-v] = it;
}
return it == iterations ? it - 1 : it; // iterate add one.
}
Arent you going to do anything? okey here is how you build
the hww-table(the circle-bounds table) and use it.
int sqrs[256];
int *sqr = sqrs + 127; // from the middle
int hw[256];
int *hww = hw + 127; // from the middle
for (i = 0; i < 128; i++) {
hw[i] = hw[255 - i] = // trigonemetrics..
(int) (128 * sin(acos ((128-i)/(float)128)));
sqrs[127 - i] = sqrs [128 + i] = (i * i);
}
re=2048*cos(0.11); // c0.re
im=1024*sin(0.4711); // c0.im
memset (f->buffer, 0, 256*256); // clear
for (y = 0; y < 128; y++){ // mirrored, remeber?
for (x = -hww[y]; x < hww[y] + 1; x++){ // from circle-edge,to
iterate (x, y, 0); // cicle-edge, iterate
}
}
Show(pix); // yupp...
Now go on and play with this and have fun. For more Source see
the full source of my FastJulia here