Today we will construct the Mandelbrot set. We will also color it and its neighborhood a little bit – with the degrees of gray. I will try to explain it in details.
In my previous post I was discussing the transformation
T(z) = z2+c
where z and c are complex numbers. We were starting with z=0 and then applying T many times, while c was fixed. We were interested at which iteration our point escaped the enchanted circle of radius 2. That is for which n T(T(T(....(T(z))))) (n times) the complex number obtained this way will have square of its modulus greater than 4.
But were investigating only the constant c of the form:
c=-0.75+0.1j
c=-0.75+0.01j
...
c=-075+0.0000001j
While making these tests we observed an interesting phenomenon – appearance of the number Pi.
Today we will increase the range of the tested values of c. The real part of c will be tested between -2,5 and 2, while the imaginary part of c will be tested between -2 to 2.
Why such a range? Simply because, as you will soon see, it is interesting. In any case only the interior of the enchanted circle is of interest. Therefore I could as well restrict myself to the range (-2,2) for the real part of c. Why -2.5, which is outside the circle? The reason is simple: I want to have a square 4 by 4, and I want to have the fractal nicely centered within this square. Simply by trial error one comes to just the selected range. You will be able to change these numbers, when you will get the method. Yet I must start with the explanation of how we will make the graphics. And I will restrict myself to the graphics of FreeBasic. Every programming language, every dialect, has its own set of commands for drawing pictures on the screen. So, let me describe how it is with FreeBasic, Version 0.20.0 Beta. A different version may need somewhat different syntax.
The graphics window consists of pixels. Certain number of pixels horizontal, and a certain number in vertical. Suppose our window has the width of pw pixels and the height of ph pixels.
In the program below I will set pw = 500 and ph = 500. It should fit within almost every screen. But you will be able to change these numbers in the Editor. Just pay attention not to exceed the resolution of your screen, as it may cause a disaster – the blue screen of death and the lost of your unsaved data. You also need to take into account the fact, the with the increase of the resolution, the time taken by the program to produce the image will increase. There will be more pixels to work upon!
First the listing of the complete program, then the explanations.
''''''''
Code:
#define stopat 512
#define pw 500
#define ph 500
screenres pw,ph,32
Declare Sub drawfractal
drawfractal
Sub drawfractal
Dim As Integer i, i_out, gray
Dim As Integer m,n
Dim As Double cr,ci,x,y,dum,a
Dim As Double delta_h=4/pw,c0_x=-2.5+2/pw, delta_v=-4/ph, c0_y=2-2/pw
Dim As Integer Ptr sp=screenptr
screenlock
For n=0 To ph-1
For m=0 To pw-1
cr=m*delta_h+c0_x
ci=n*delta_v+c0_y
x=0
y=0
i=0
i_out=stopat
Do While i< stopat
i=i+1
dum = x*x-y*y+cr
y=2*x*y+ci
x=dum
a=x*x+y*y
if a>4 then
i_out=i
Exit Do
end if
Loop
If i_out<stopat Then
gray=Abs((i_out And 63)-31)*7+31
*(sp+n*pw+m)=rgb(gray,gray,gray)
Else
*(sp+n*pw+m)=0
End If
Next
Next
screenunlock
Sleep
End Sub
'''''''''''''''''''''''''''
In the program there are two integer variables m,n. The variable m takes pw values from 0 to pw-1, while n takes ph values from 0 to ph-1.
m=0, n=0 is the pixel in the left upper corner of the graphics window
m=pw-1, n=ph-1 is the pixel in the right lower corner
We will column by column, row by row, from the left to the right, like it is when we are reading text.
n=0, m=0,1,2,...,pw-1 is the the first row at the top
n=1, m=0,1,2,..., pw-1 is the second row
and so on. Finally
n=ph-1, m=0,1,2,..., pw-1 is the last row of the graphics window. The last pixel.
Now we have to translate each pixel into the complex number c represented by this pixel. How to do it? A simple calculation is needed We have pw pixels in a row. And we have chosen 4 to be the with of our window in real numbers. Therofore each pixel is 4/pw wide on the x-axis. Since we have chosen the square 4x4, each pixel will have the height 4/ph on the z-axis.
Now I will decsribe the program drawing the fractal. We start with the declaration of the constants – parameters that will not change during the run. In the Pi program we were declaring then via „const”. But we can declare them also this way:
#define stopat 512
#define pw 500
#define ph 500
The constant sopat declares the maximal number of iterations of the transformation T. It tells us how many times we are going to repeat the transformation until we give up if the point is still inside the circle. If the point does not escape after 512 iterations we will treat it as „never escapes” and we will paint black the appropriate pixel, representing the value of c leading to such a behavior. These points will form the fractal – the Mandelbrot set. You may say it is somewhat „risky”. Well, it is risky, but it works. You can try the value stopat equal 5, then 50, then 500. You can even go beyond 500, the the execution time will be longer. But not only that, since the precision of the calculations will start to play a role. To relay on calculations where a given number is taken to a power of 1000 would be unreasonable with the finite precision of our hardware and software.
In the next line we declare the screen resolution and the color depth.
screenres pw,ph,32
We already know what are pw and ph. The number 32 means that we will be working with the 32 bit TrueColor. That means, in parctice, that each of the colors Red, Green, Blue will be taking values between 0 and 253. Nowadays almost every graphics card has such a depth. The coloring (by the degrees of gray alone) will take place at the end of the pixel loop.
One comment: since our program will operate on the computer's memory, there is always a danger that something unexpected will happen. For instance your computer may have a smaller color depth than the one assumed. It is hard to predict what will happen in such a case. Therefore, before running the program, please save all of your unsaved data, those that you would not want to loose in the case of a crash, when you will have to restart your PC. By the way this warning applies to ANY program that you are running the first time. Even when it is the program originating from a renowned company.
The calculations and the coloring will be contained in a subroutine. This is not really necessary now, but it will be handy in the future, when we will learn how to ZOOM at the particular region of the fractal boundary with a click of the mouse. Therefore, having in mind the future applicationsm we declare the subroutine:
Declare Sub drawfractal
This subroutine does not need any parameters. But in the future it will have the parameters of the mouse click point.
Then we call the subroutine, we tell it to run.
drawfractal
Now I will describe what this subroutine
Sub drawfractal
.....
End Sub
is doing. This subroutine is doing all the work. It starts, agin, with declaration of variables that are being used by the subroutine.
Dim As Integer i, i_out, gray
Dim As Integer m,n
Dim As Double cr,ci,x,y,dum,a
Dim As Double delta_h=4/pw,c0_x=-2.5+2/pw, delta_v=-4/ph, c0_y=2-2/pw
Dim As Integer Ptr sp=screenptr
The variable i will take integer values and will number the iterations.
The variable i_out, also integer, will note the actual value of i when the escape happens.
The variable gray will assign the gray level to each pixel. It will depend on the value of i_out.
Roughly, the sooner the point escapes, the lighter will be the color. But you will see at the end, that we will, in fact, circulate colors. For a good reason that will be explained.
The variables m,n – integers will number the columns and the rowas of the window – as I have described above.
The variable cr,ci will be real numbers in double precision (16 significant digits) and will serve us for storing the real and the imaginary part of the constant c. (c will be a constant for each pixel, but when we will move to the next pixel, the value of c will change).
The variables x,y will be used for storing the real and the imaginary part of the transformed complex number z.
The varaibale dum will be like in the Pi program.
Next we have the lines:
Dim As Double delta_h = 4/pw, delta_v = -4/ph
Dim As Double cr = -2.5+2/pw, ci = 2-2/pw
Here we declare four more variables, initializing them at the same time, that is giving them the starting values.
delta_h=4/pw and delta_v=-4/ph have already been explained above. But why the minus sign in front of delta_v? The point is that rows of the graphics window are being numbered from the top to the bottom. But the y-coordinate, normally, increases from the the bottom to the top. Therefore with the increase of the row number n, the value of y will decrease. Thus the minus.
What remains to be explained are the mysterious c0_x=-2.5+2/pw and c0_y=2-2/pw. These are the Cartesian coorinates assigned to the first top left pixel. The top left corner of the window has cartesian coordinates (-2,5,2). The bottom right of this pixel has the coordinates (-2.5+delta_x, 2-delta_y). Now, you may just calculate the coordinates of the center of this pixel using the values delta_w=4/pw and delta_h=4/ph – what will be the coordinates of this center? This is a very small correction, but it should improve the precision of the shape of our fractal at the pixel level. One needs to keep it in mind that some points are very sensitive to the value of the constant c. One such example we met in the previous program, in the neighborhood of c=-075+0.1j. We have seen that the change of only 1/100 in the imaginary part changes the number of interations needed for the escape 100 times. It is the points near the fractal boundary that are so sensitive. And the boundary is not smnooth at all – it is a fractal!
Finally we need to explain the line:
Dim As Integer Ptr sp=screenptr
I do not want to go into the details here, but this line semply reserves a region in the RAM of the machine where the image will be stored before all the calculations will be done and the whole picture shown at once. We could draw on the screen pixel after pixel, but it would slow down the program. We will see below how this „drawing in the memory” happens.
But first we lock the graphics window with the command
screenlock
We will unlock it when all calculatiuons are done. Then the whole picture almost instantly will appear on our screen.
Next we start working on the pixels, one by one. For this a double loop is being used:
For n=0 To ph-1
For m=0 To pw-1
........
Next
Next
After we are done with all our pixels, we will unlock the screen with the command
screenunlock
And we will put the image to sleep, so that it would disappear right after showing us the shape.
Sleep
This way of producing graphics is particular for the FreeBasic dialect. Each programming language has its own ways of dealing with the graphics.
Now let us move inside the loop and see what is going on there. With m,n given we are working on the pixel in the m-th column and n-th row of the graphics window. What are the Cartesian coordinates assigned to the center of this pixel? By making a simple sketch you will convince yourself that these coordinates are
cr=m*delta_h+c0_x
ci=n*delta_v+c0_y
And this is the value of c assigned to this pixel. Next we go like in the previous Pi-program, with slight changes. First we assign the initial values to the real and imaginary part of the starting point of our iteration, z. We also initialize the iteration variable i:
x=0
y=0
i=0
At the beginning we have 0 iterations.
The we have one more initializing line:
i_out = stopat.
We assign to the variable i_out the maximal allowed value. That means that a priori we consider each pixel black – as belonging to the fractal. Then we start the iterartion loop and if the point will escape before, we will change the value of i_out and the color of the pixel. The sooner it will escape the lighter will be the color (yet not exactly, as I will explain later).
Then we start the iteration loop, the same way as in the previous program, but with slight changes:
Do While i< stopat
i=i+1
dum = x*x-y*y+cr
y=2*x*y+ci
x=dum
a=x*x+y*y
if a>4 then
i_out=i
Exit Do
end if
Loop
First of all the loop starts with
Do While i< stopat
i=i+1
insetad, as it was before, with
Do
i=i+1
Why? The point is that the previous Do was not safe. We were asking the prgram to increase the value of i without telling it when to stop. And what if it will not stop? In such a case the value of i can become greater than the maximal integer allowed by our software and hardware and the program would crash. Or, if the calcualations are slow, after one hour still nothing will appear and we will have to force the program to shut down or to reset the machine. This time we are taking the safety measures. And anyway, we are going to consider „stopat” iterations as „infinity”. Nothing happend till i=stopat? We paint the pixel black and move to the next pixel.
Then it goes like in the previous listing, with some changes
Dalej leci jak w poprzednim programie, z mała tylko zmianą:
Do While i< stopat
i=i+1
dum = x*x-y*y+cr
y=2*x*y+ci
x=dum
a=x*x+y*y
if a>4 then
i_out=i
Exit Do
end if
Loop
When the point jumps outside the circle, we note the value of i, store it in i_out, and stop the iterations. What remains is given a color to the pixel (m,n). This is done in the piece:
If i_out<>stopat Then
gray=Abs((i_out And 63)-31)*7+31
*(sp+n*pw+m)=rgb(gray,gray,gray)
Else
*(sp+n*pw+m)=0
End If
This piece is not so important for the beginners. In fact I took this part from another program and I was trying to understand what it is doing. I think I understand it now, I could change it, but I decided to leave it as it is. If you know what are bitwise operations, and if you understand what you are doing – you can play with it. Otherwise leave it unchanged.
Let me say that it is not quite true that the more iterations it takes for the point to escape, the darker a given pixel is being painted. That would be unreasonable. We want to see sharply the boundary of the fractal. But close to the boundary the value of i_out changes very fast. On the other hand the human eye will not distinguish between black and „almost black”. So the boundary of the fractal would look fuzzy. Therefore, instead of making pixels darker and darker with approaching the boundary from the outside, we are circulating 33 gray levels. These levels are given by the numbers:
31,38,45,52,59,66,73,80,87,94,101,108,115,122,129,136,143,150,157,164,171,178,185,192,199,206,213,220,227,234,241,248,255
As you can see we change them by seven. 255 is white, 0 is black. Black is reserved for the points that will „never” escape. When increasing i_out by one, we move to the next gray level from the list. When we come to the end, we start again from the beginning. This is what the formula
gray=Abs((i_out And 63)-31)*7+31
does. I do not want to enter into the „bitwise And”. Enough to note that the number 63 in the binary form is 111111 – all ones. That ends the part describing how gray levels are being assigned to the values of i_out.
Finally the last line:
*(sp+n*pw+m)=rgb(gray,gray,gray)
Here sp is address of the memory part reserved for the our image. Not going the details, enough to say that sp+n*pw+m is the address of the pixel we are working on. We assign to this pixel the calculated gray level. Gray means the same amount of Red, Green, and Blue – rgb.
And that is all for now. You can play changing the values of pw, ph and stopat. You make sure never exceed the actual resolution of your screen.
And o not forget that after copy/paste operation of the listing into your Editor, you will have to save it on the disk and give it a name. I called mine "mymandel.bas”. Have fun.
And if have made some mistake, please correct.
Here is my result:
, but it's fun anyway!
