The Mandelbrot set is a very well known and interesting mathematical topic that is fun to explore. It's also a very computationally expensive one as well. To produce a satisfying image of the Mandelbrot set requires a number of pixels in the order of millions, each requiring thousands of calculations. Even then, the produced image is only an approximation. The image can be enhanced, zoomed, and panned to relay more of the infinite information to the user who is exploring it. There is even more information when one considers the Julia Set, which is related to the Mandelbrot Set. Displaying that information adds more computations but adds even more data for the user to explore. My goal was to make this process as pain free and fast as possible.
The program I made to achieve this goal was written in C using the Simple DirectMedia Library (libsdl.org) library. Both C and SDL have a good reputation of being fast and portable. The program should run on Windows, Mac OS X, and Linux machines for the foreseeable future. In addition, the program can optionally make use of OpenCL (khronos.org) (wikipedia.org) to take advantage of the parallel processing nature of today's processors (CPU) and especially graphics cards (GPU). At the time of writing, OpenCL is rather hard to get working on any computer and requires a scrumptious new graphics card, so it is turned off by default, but the code can still use multiple CPU cores without OpenCL. Without OpenCL, the program can not take advantage of the graphics card's (often considerable) computing power.
Here's what it looks like rendered using OpenCL.
Keep in mind that the Julia set (the one on the right) is never rendered outside a single CPU core using standard C code because that process could not be parallelized by me at this time.
Math behind the Mandelbrot Set and Julia is very easy to come by, as are algorithms for producing the images in the program. Therefore, I will keep it short because others can likely explain it better than I.
A complex number c is represented by the following equation:
where:
The magnitude of a complex number is:
The square of a complex number is:
The Mandelbrot set is defined as the set of complex numbers, where any complex number in the set, let's say c, keeps z bounded, where z is recursively defined as:
Bounded in this case means that it never goes off to infinity. For the set builder notation lovers:
At each pixel we render on the Mandelbrot, we have a complex constant that represents that pixel's position in the complex plain. Let's call that c. First we define that:
We compute the next value of z using the following equation.
Then we check that the magnitude of this iteration of z is less than 2. This is a cheap way of performing a bounded check.
If so, we compute the next value z and repeat the steps following that computation, including this one, until n reaches adjustable, yet constant maximum iteration count. If the max is reached, c is assumed to be in the Mandelbrot set and the corresponding pixel is colored white. If the magnitude check fails, the n that it divided by the max number of iterations. The result of that is multiplied by an arbitrary color and that product is set as the color of the corresponding pixel. The image produced is a color based display of how close a map of complex numbers are to the Mandelbrot Set. Brighter pixels are closer to the Mandelbrot Set. White pixels are so close to being in the Mandelbrot set that the computer can not tell without further iterations.
The Julia Set is harder to explain and I'm afraid I don't fully understand it. I will not poison your mind with my guesses. The method of computing I used is finicky and took trial and error to get it. For the morbidly curious I present it here.
Let's define the square root of a complex number as:
First choose any complex number c. Let z be defined as the following.
The negative or positive square root is chosen at random. Each iteration of z is plotted on a pixel corresponding to the point on the complex plane. The maximum number of iterations should be high, but certain points that should be there will never (as far as I can tell) be touched by z.
For your convenience, gentle reader, I have included ways of configuring the program in the following table.
| Property | Description | How to Configure |
|---|---|---|
| Thread Count | The number of threads to use when rendering the Mandelbrot set when OpenCL is not used. It has no effect otherwise. It will improve performance to have more threads than actual cores. My machine had Intel Core i7 machine with 4 hyper-threaded cores worked best at around 16 threads. Tune this value to your own machine for best performance. | Run the program in the command line with the number of threads you desire as the first argument. An example is: ./mj 16 |
| Iteration Count | The number of maximum number of iterations to use when rendering the Mandelbrot set. A higher number produces a darker image and takes a longer amount of time, but it is more accurate and it is the only way to see more detail upon very deep zooming. If this value is unset, it starts at 50 and increments by one after every render to produce a more detailed image over time. A value of 2000 will run slowly even on a graphics card, but produces a very detailed image on max zoom. A value of 50 looks nice for unzoomed viewing. I encourage you to experiment with this value because it's really run. Just keep the value positive and below 5000 to avoid turning your computer into a slide show. | Run the program in the command line with the number of threads (even when using OpenCL; it's a bug) you desire as the first argument and then enter the iteration count as the second argument. An example is: ./mj 16 50 |
| Width * | The width of the entire rendered window in pixels. The Mandelbrot will be half this width. | Change line 13 of mandyjule.c to #define WIDTH INSERT_WIDTH_HERE |
| Height * | The height of the entire rendered window in pixels. The Mandelbrot will be this height. | Change line 17 of mandyjule.c to #define HEIGHT INSERT_HEIGHT_HERE |
| Fullscreen * | Render the entire window in fullscreen. The given Width and Height will be ignored and replaced with the screen resolution. | By default it is disabled (except when using the makefile like in Windows). To enable it add -DFULLSCREENto the command line that starts with gcc. (If you are using the makefile, remove -DFULLSCREEN from line 4 of makefile to disable fullscreen.) |
| OpenCL * | The choice to use OpenCL or not. OpenCL must be installed on the system. This has only been tested on Windows. | To enable it use make useclin windows instead of the normal make |
* indicates that the project must be rebuilt for the configuration to take effect.
Note: Omit the ./ at the beginning of the command lines in the examples in the table if you are using Windows, unless you are using the msys shell as your current command line
The program can be built from source on Windows, Mac OS X, and Linux as well as some other Unix-like systems supported by SDL. As far as I know, my code is standard C99 and will compile with minimal modification on whatever compile you happen to be using on your platform. I have personally compiled the code on a Windows 7 64 bit machine using the MinGW (mingw.org) compiler and using gcc version 4.0 on a Mac Mini running Mac OS X 10.5.8. For those of you reading this decades from now, those Mac OS X machines were what we had to use in the labs, and even then, they were old. The code will very likely compile elsewhere, but don't be afraid to wrench around the innards if it doesn't compile the first try. I tried to keep the C code simpilish for your benefit.
The program requires SDL version 1.2.x to be installed. The latest source code can be downloaded from: http://www.libsdl.org/release/SDL-1.2.14.tar.gz. Run the following to download, unpack, configure, compile, and install SDL.
$ wget http://www.libsdl.org/release/SDL-1.2.14.tar.gz $ tar -xvf SDL-1.2.14.tar.gz $ cd SDL-1.2.14 $ ./configure $ make $ make installTo confirm that the install worked, run:
$ sdl-config --static-libs --cflags
That should spit out useful compiler flags. That confirms SDL is installed on your system.
Once that is done. You should be able to compile the mandelbrot.c file with:
$ gcc mandyjule.c -std=gnu99 `sdl-config --static-libs --cflags` -o mj $ ./mj
That last line will run the program.
Now some bad news, gentle reader. I had no access to a Mac or Linux machine with OpenCL on it. I know Mac OS X 10.6 and higher includes OpenCL though. There is no reason you can not get it working given you know you're way around gcc, or at least Google. If you get it installed on Linux or have a compatible Mac the following compile command should get you started:
$ gcc mandyjule.c -std=gnu99 -DUSE_OPENCL `sdl-config --static-libs --cflags` -o mj
I doubt it will work the first try because you'll have to link the OpenCL library, but be sure to include the
-DUSE_OPENCLfor that part enables the OpenCL code within the C source file.
In Windows I made sure it would build using the compiler provided with aszguard. If you don't have it you can get it from here, but it links to ~gfrancis's site, so don't expect it to remain a live link or for it to remain up to date. If you prefer to not use aszguard, a normal MinGW (mingw.org) will work, but this option is more complicated and I trust that if you use that, you know what you're doing. From here I will assume you are using aszguard.
First thing we must do is determine if you want to enable OpenCL support. If your card supports it, it will really speed up the application, so it is encouraged. If you have no graphics card or it does not support OpenCL, the code can take advantage of the multiple cores in your system without OpenCL, so there is no reason to go to the trouble of enabling it in that case.
I have only personally had it working on an AMD Radeon 6950 inside my desktop. Aside from the drivers, I also needed the AMD Accelerated Parallel Processing SDK for headers and libraries. The name may have changed by the time you read this. Google is your friend in such a case. In the case of an NVIDIA card try the Cuda Zone for relevant headers and libraries. If you have something else, Google it. Remember, you will need a cl.h and an libOpenCL.a or OpenCL.a from whichever SDK you acquire. I know for a fact that the AMD SDK I mentioned installs these.
No matter what path you take, you will need SDL. Grab the SDL 1.2.x development libraries from the SDL website. I used SDL-devel-1.2.14-mingw32.tar.gz. Unpack that into your aszguard directory using your favorite unzipping program and fire up msys_shell.bat from the top level of the aszguard directory. Then type the following the command windows that opens.
$ cd SDL-1.2.14 $ make install
To confirm that SDL installed run:
$ sdl-config --static-libs --cflags
Hopefully this output compiler flags. If so, the install worked.
Now type the following to compile and run the program.
$ cd /path/to/project/folder # Replace that with the path to this project's folder $ make usecl # Omit the usecl if you are choosing to not use OpenCL $ ./mj
If the build was successful the program should now be running. If it fails because it can not find a library or an CL/cl.h, edit the makefile in the project directory and uncomment (remove the #) on lines 6 and 7. Replace the paths with the correct paths from your SDK. Be sure to keep the -L and the -I parts in the makefile.
Once the program is up and running, using it is very straightforward. It can best be summed up with the following table
| Action | Input |
|---|---|
| To Exit | Press ESC |
| To Zoom In | Click to where you wish to zoom to |
| To Zoom Out | Right Click anywhere (be patient, it takes 10 second to zoom out no matter how zoomed in you are) |
| To Reset Iterations Count to 0 | Press SPACE |
| To Configure | See Configuration Above |
For those curious about GPU programming with OpenCL, the following code listing is the kernel that runs on the GPU that computes the number of iterations for each point in the complex plane that is displayed.
#pragma OPENCL EXTENSION cl_amd_fp64 : enable kernel void mandelbrot( const double2 min, const double2 max, const int2 steps, const int max_iterations, global int* buffer ) { double2 global_pos = (double2)(get_global_id(0)/(double)steps.x, get_global_id(1)/(double)steps.y); double2 complex_pos = mix(min, max, global_pos); double2 current_pos = (double2)(0,0); int iterations = 2; while (dot(current_pos, current_pos) < 4 && iterations < max_iterations) { double newReal = current_pos.x * current_pos.x - current_pos.y * current_pos.y + complex_pos.x; current_pos.y = 2 * current_pos.x * current_pos.y + complex_pos.y; current_pos.x = newReal; iterations++; } if (dot(current_pos, current_pos) < 4) iterations = -1; buffer[get_global_id(0) + get_global_id(1) * steps.x] = iterations; };