Project of C/C++ 2011 Spring, CS, HIT
Interactively Exploring the Mandelbrot Set
The exercise in a nutshell
You will write a program to draw the Mandelbrot Set; you will be able to interactively select a region of the Mandelbrot Set to zoom in on, to view it in greater detail.
This lab only requires the use of very simple C. You will be writing some mathematical expressions and some control structures, manipulating some variables, and calling functions that have been provided for you, in a library called OpenCV. You will need the skills you picked up in the lab classes to compile, link and execute your program.
The instructions for this exercise come in three parts: first comes a general introduction and background to the lab, and you are expected to read this thoroughly before you attend your scheduled lab session; the second part is a ‘manual’ for the library functions you have been provided with, and you should make sure you are familiar with this as well before turning up; the third and final part is a detailed description of the individual tasks that comprise the lab. Again it is important that you read all of this at least once before you turn up for your lab, so that there are no surprises when you sit down at your machine. (Some of the tasks—such as planning the mappings between co-ordinate systems—are much easier to figure out on paper than they are to do when you are sitting in front of a computer.)
Part 1: Introduction
Before your scheduled lab session, you should spend some time reading, in conjunction with your lecture notes, the following description of the Mandelbrot Set. It is very important that you do this: the time allocated for this lab assumes you are ready to start on the programming part of the exercise as soon as you turn up.
1
Project of C/C++ 2011 Spring, CS, HIT
About the Mandelbrot set
The Mandelbrot Set – named after the mathematician Benoit Mandelbrot, who discovered it on 1st March 1980 – is a graph of a mathematical formula. It is also a fractal, the name given to shapes that have a property known as ‘self similarity’. Examining fractals at different levels of detail reveals similar characteristics. A good example of this is a piece of yummy broccoli: if you divide a head of broccoli into individual florets, you’ll see that they look pretty much the same as the whole head. Many natural phenomena such as mountains or coastlines have self-similar properties. Self-similarity in the natural world, however, has its limits – you can only split a broccoli floret up a certain number of times before the little bits of green vegetable you are left with bear no relationship to their original shape. The Mandelbrot Set however, being a purely mathematical concept, has no such limits, and you can examine it at ever-increasing levels of detail and it will always look intricate. (In purely practical terms, the only limit to the level of detail at which you can explore the Mandelbrot Set is the precision level of the real numbers of your computer.)
In order to complete this exercise, you will need to know a little about the mathematics of the Mandelbrot Set itself, since you will be writing some code to calculate it. First, you need to know a little about complex numbers.
Complex numbers
A complex number comprises a pair of values: a real part and an imaginary part. The real part is an ordinary number, for example, -3. The imaginary part is a real number multiplied by a special number called i and written, for example, as 2i. An example of a complex number would be -3 + 2i. The number i was invented to make it possible to manipulate expressions that involve the square root of negative numbers; i is defined as being the square root of –1, so i2 is equal to –1, and the square of, (for example) 3i is –9 (3 * 3 * i * i).
By the way, don’t worry about trying to understand what the square root of a negative number actually means.
Real numbers can quite easily be represented on a one-dimensional line (called, not surprisingly, the real number line). Negative numbers are plotted on the left of the zero origin, and positive numbers are plotted to the right. Complex numbers have two components (real and imaginary), so need to be plotted on a two-dimensional plane instead, with the real values along the x axis, and the imaginary values on the y axis, as we can see in Figure 1. This is called the complex plane.
2
Figure 1: The Complex Plane
The mathematics of the Mandelbrot Set
Figure 2 shows a visualization of the Mandelbrot Set.
Project of C/C++ 2011 Spring, CS, HIT
Figure 2: The Mandelbrot Set
The Mandelbrot Set is a set of complex numbers. To determine whether a given complex number C belongs to the set or not, we need to perform a simple test on it. This test involves plugging C into the following equation:
Zn+1 = Zn2 + C
Here, Zn is another complex number, which starts off as Z0 = 0 + 0i (this is the complex number equivalent to the real number 0). The above equation is an iterative equation; in other words, after computing the equation we get a new value of Z, which we then plug back into the equation, and repeat the process.
To illustrate, let’s look at how this equation behaves, for some value of C (we’ll choose C = 1.3 + 0.4i, as a random example). Remember we always initialise Z0 to be
3
Project of C/C++ 2011 Spring, CS, HIT
0 + 0i:
Z0 = 0 + 0i Z1 = Z02 + C
= Z02 + (1.3 + 0.4i) = 0 + 1.3 + 0.4i
= 1.3 + 0.4i
Z2 = Z12 + (1.3 + 0.4i)
= (1.3 + 0.4i)2 + (1.3 + 0.4i) = 1.53 + 1.04i + (1.3 + 0.4i) = 2.83 + 1.44i
Z3 = Z22 + (1.3 + 0.4i) = 7.2353 + 8.5504i
As you can see, the value of Z changes each time we iterate the function. This doesn’t seem to be very interesting, but in fact this simple equation behaves in a rather unexpected way. We can see this, if we examine what happens to the modulus of Z, after each iteration. The modulus (sometimes also called magnitude, or absolute value) of a complex number Z = a + ib, is written as |Z|, and is defined as the distance of Z from the origin (0, 0) of the complex plane. It’s easy to calculate: |Z| = sqrt (a2 + b2) – it’s just the Pythagoras rule.
As we iterate the function, a number of things could potentially happen to the value of |Z|:
1) it could jump around chaotically
2) it could get bigger and bigger without limit
3) it could tend to zero or stabilize at some other constant value
It turns out that it is the particular value of C which is the deciding factor for the
behavior of |Z| (remember, Z always starts off as Z = 0 + 0i).
And this leads to the definition of the Mandelbrot Set: Values of C that cause |Z|
to tend to some fixed value, or to jump around chaotically, are members of the Mandelbrot Set; values of C that cause |Z| to get bigger and bigger and tend towards infinity are not members of the Mandelbrot Set.
In fact, when we iterate our equation, it turns out that one of only two specific things happens to the value of |Z|: it either stays equal to or below 2 forever, or it eventually becomes greater than 2, and will then continue to tend to infinity. This might seem a little odd, and indeed it is, but we don’t need to worry about why for the purposes of this exercise.
In case you’re wondering, in the example above, the successive values of |Z| were (to 2 decimal places): 0, 1.36, 2.32, and 4.21. So in this case, the complex number C = 1.3 + 0.4i was not a member of the Mandelbrot Set. This may give you a clue as to why the Mandelbrot Set didn’t become an object of interest until the computer age – the calculations required are long and tedious to do by hand!
So, now we can turn our attention to visualizing the Mandelbrot Set. All we need to do is to test enough complex numbers in the complex plane, and for each we shall note whether it is in the Mandelbrot Set or not. If it is, we shall plot a pixel for it; if it isn’t, we won’t.
4
(0, 0)
Pixel (i, j)
Project of C/C++ 2011 Spring, CS, HIT
Part 2: Drawing the Mandelbrot Set
Coordinate systems and pixels
The first thing to note is that the complex plane is a two-dimensional space, so this will map very conveniently onto a two-dimensional window on our desktop display.
Further, we can easily find a mapping between every point in our screen window (we’ll call these points ‘pixels’ for now), and the corresponding complex number in the complex plane, as Figure 3 shows.
(0, 0)
Complex number (re, im)
The pixel window, where the The pixel window, where the origin The complex plane, where origin is on the top left corner. is the middle of the window. complex numbers live
Figure 3: mapping from complex numbers to pixels
In order to display the Mandelbrot Set, we therefore have three main tasks:
1. find some way of drawing to our screen
2. choose a range of complex numbers and to calculate whether they belong to
the Mandelbrot Set or not; the range of values we will choose is the square
region of the complex plane defined by (–2, –2i) to (2, 2i)
3. map between our set of values and the screen co-ordinates
The first of these tasks—that of getting the entire infrastructure in place to create
a window on the screen—turns out to be a fairly complex job, so we have provided you with a library of functions that simplifies this, called the OpenCV library.
One of your tasks will be to understand how to link your code that you have written in C with the library we have provided. The main body of this exercise revolves around writing your own code to evaluate the Mandelbrot Set and to map your data structures onto the library we have provided in order to make the set appear on the screen.
You are provided with functions to:
create the window –cvNamedWindow create image- cvCreateImage
show image in the window-cvShowImage
(0, 0)
Pixel (i, j)
5
Project of C/C++ 2011 Spring, CS, HIT
set the color of the pixel-cvSet2D save image-cvSaveImage
copy image-cvCloneImage
Using the OpenCV library
Here, we’ll describe each of the data structures, constants and functions OpenCV have provided you with. You don’t need to understand how these functions actually work internally – you just need to understand what they do, and how to use them.
First the handy data structures:
typedef struct CvPoint {
int x;
int y; }CvPoint;
typedef struct CvScalar {
double val[4]; }CvScalar;
typedef struct _IplImage {
int nSize;
int ID;
int nChannels;
int alphaChannel;
int depth;
char colorModel[4]; char channelSeq[4]; int dataOrder;
int origin;
int align;
int width; int height;
/* sizeof(IplImage) */
/* version (=0)*/
/* Most of OpenCV functions support 1,2,3 or 4 channels */
/* ignored by OpenCV */
/* pixel depth in bits: IPL_DEPTH_8U, IPL_DEPTH_8S, IPL_DEPTH_16S, IPL_DEPTH_32S, IPL_DEPTH_32F and IPL_DEPTH_64F are supported */
/* ignored by OpenCV */
/* ditto */
/* 0 – interleaved color channels, 1 – separate color channels. cvCreateImage can only create interleaved images */
/* 0 – top-left origin,1 – bottom-left origin (Windows bitmaps style) */
/* Alignment of image rows (4 or 8).OpenCV ignores it and uses widthStep instead */
/* image width in pixels */ /* image height in pixels */
6
struct _IplROI *roi;
struct_IplImage*maskROI; void *imageId;
struct _IplTileInfo *tileInfo; int imageSize;
char *imageData;
int widthStep;
int BorderMode[4]; int BorderConst[4]; char *imageDataOrigin;
/* image ROI. if NULL, the whole image is selected */
/*mustbeNULL*/
/* ditto */
/* ditto */
/* image data size in bytes (==image->height*image->widthStep
in case of interleaved data)*/
/* pointer to aligned image data */
/* size of aligned image row in bytes */
/* ignored by OpenCV */
/* ditto */
/*pointer to very origin of image data (not necessarily aligned) -needed for correct deallocation */
Project of C/C++ 2011 Spring, CS, HIT
}/*here, the items in grey font color will not be used in this work*/
Now the functions:
cvCreateImage
The function cvCreateImage creates the header and allocates data.
IplImage* cvCreateImage (CvSize size, int depth, int channels)
size
Image width and height
depth
Bit depth of image elements. Can be one of:
IPL_DEPTH_8U – unsigned 8-bit integers
IPL_DEPTH_8S – signed 8-bit integers
IPL_DEPTH_16U – unsigned 16-bit integers IPL_DEPTH_16S – signed 16-bit integers IPL_DEPTH_32S – signed 32-bit integers IPL_DEPTH_32F – single precision floating-point numbers IPL_DEPTH_64F – double precision floating-point numbers
channels
Number of channels per element (pixel).Can be 1, 2, 3 or 4.The channels are interleaved, for example the usual data layout of a color image is:b0 g0 r0 b1 g1 r1…
Although in general IPL image format can store non-interleaved images as well and some of OpenCV can process it, this function can create interleaved images only
cvRectangle
The function cvRectangle draws a rectangle with two opposite corners pt1 and pt2.
void cvRectangle (CvArr*img, CvPoint pt1, CvPoint pt2, CvScalar color,int thickness=1, int line_type=8, int shift=0 )
img
7
One of the rectangle vertices
Opposite rectangle vertex
Project of C/C++ 2011 Spring, CS, HIT
Image,
pt1
pt2
color
Line color (RGB) or brightness (grayscale image)
Thickness
Thickness of lines that make up the rectangle. Negative values, e.g.
CV_FILLED, make the function to draw a filled rectangle
line_type
Type of the line, see cvLine description
shift
Number of fractional bits in the point coordinates
cvCloneImage
The function cvCloneImage makes a full copy of the image including header, ROI and data
IplImage* cvCloneImage( const IplImage* image );
image
Original image.
cvNamedWindow
The function cvNamedWindow creates a window which can be used as a placeholder for images and trackbars. Created windows are referred by their names. If the window with such a name already exists, the function does nothing.
int cvNamedWindow( const char* name, int flags=CV_WINDOW_AUTOSIZE );
name
Name of the window which is used as window identifier and appears in the
window caption
flags
Flags of the window. Currently the only supported flag is CV_WINDOW_AUTOSIZE. If it is set, window size is automatically adjusted to fit the displayed image, while user can not change the window size manually.
cvShowImage
The function cvShowImage shows the image in the specified window. If the window was created with CV_WINDOW_AUTOSIZE flag then the image is shown with its original size, otherwise the image is scaled to fit the window.
void cvShowImage( const char* name, const CvArr* image );
name
Name of the window
image
Image to be shown
cvSaveImage
The function cvSaveImage saves the image to the specified file. The image format is
8
Project of C/C++ 2011 Spring, CS, HIT
chosen depending on the filename extension. Only 8-bit single-channel or 3-channel (with ‘BGR’ channel order) images can be saved using this function. If the format, depth or channel order is different, use cvCvtScale and cvCvtColor to convert it before saving, or use universal cvSave to save the image to XML or YAML format.
int cvSaveImage( const char* filename, const CvArr* image );
filename
Name of the file
image
Image to be saved
cvDestroyWindow
The function cvDestroyWindow destroys the window with a given name.
void cvDestroyWindow( const char* name );
name
Name of the window to be destroyed
cvReleaseImage
The function cvReleaseImage releases the header and the image data
void cvReleaseImage( IplImage** image );
image
Double pointer to the header of the deallocated image
Part 3: The lab exercise
We’re finally at the point where we can detail your tasks in the laboratory exercise.
Task 1
Create a file called ‘mandelbrot.c’ or ‘mandelbrot.cpp’ using whichever editor you prefer. Write and compile yourself a ‘main’ function that prints out the name of your favorite root vegetable just to be sure that everything is working ok.
Task 2
Add the following function:
int mSetTest(doublec_re, doublec_im) {
/*your stuff goes here*/ }
This function will take the complex number C, whose real part is c_re, and whose imaginary part is c_im, and will determine whether or not C is a member of the Mandelbrot Set. As you will recall from the discussion above, this will involve
9
Project of C/C++ 2011 Spring, CS, HIT
iterating the following equation, where Z0 is defined to be (0 + 0i): Zn+1 =(Zn)2+C
Let’s recall the definition of the Mandelbrot Set given above: if, during all the function iterations, |Z| stays less than 2.0, then C is in the Mandelbrot Set; but as soon as |Z| exceeds 2.0, we immediately know that C is not in the Mandelbrot Set.
This immediately leads to a practical programming problem. Suppose we are testing a particular C that just happens to be a member of the Mandelbrot Set. Then, no matter how many iterations we do, |Z| will always be less than 2.0 (by definition). So, should we continue iterating the function forever…? Obviously that would be silly. What we will do is to make a pragmatic decision, and say that if, after 200 iterations, |Z| is still less than 2.0, then it’s a pretty good bet that C is a member of the Mandelbrot set. In such a case, your function should return a value of 0.
But what if a particular C is not in the Mandelbrot Set? In this case, |Z| will soon exceed 2.0, after which there is no point doing any further iterations. In this case, your function should return a value which corresponds to the number of iterations performed before |Z| exceeded 2.0. Of course, this value will always lie between 1 and 200.
Once you have written this function, you must test it on the following values of C:
C
(2.0, 2.0) 1 (0.0, 0.0) 0 (-1.0, -1.0) 3 (1.0, 1.0) 2
Task 3
Now you’re ready to draw the Mandelbrot Set, so you’ll need to link your program with the OpenCV library we’ve given you:
http://cms.hit.edu.cn >>I Can C You>>”Share Your Code” Event–2>> OpenCV
Result
A. Create an image and show it in a window
In your main function, add the following lines.
IplImage *pImg =
cvCreateImage (/*size*/cvSize (IMAGE_SIZE, IMAGE_SIZE), /*depth*/ 8, /*n Channels*/3);
cvNamedWindow (/*name of the window*/”mandelbrot”, 1);
cvShowImage (“mandbrot”, pImg);
cvWaitKey (/*delay*/0);
cvDestroyWindow (“mandlbrot”);
You will also need to add the #include directive to the top of the file to get
cvRelease(&pImg);
the 10
Project of C/C++ 2011 Spring, CS, HIT
prototype for this function from the header file (cv.h, highgui.h), and you’ll need to link your compiled executable with the OpenCV library. Compile this program, and run it. A grey window should appear on the screen. You’ll have to cause the program to exit by pressing any key (act by cvWaitKey).
B. Mapping from pixels to complex numbers
You’re now going to actually color in the Mandelbrot Set. You will need two nested ‘for’ loops, one that iterates over ‘x’ dimension pixels and another that iterates over ‘y’ dimension pixels. For the time being, vary both these values from 0 to IMAGE_SIZE in both dimensions. Within the loops, you will have to map from these pixel co-ordinates to the coordinates of the Mandelbrot Set. Remember your pixel values are being varied between 0 and IMAGE_SIZE, and your Mandelbrot Set values vary between –2 and +2 in both the real and imaginary dimensions.
C. Color the Mandelbrot Set
Once you’ve written an expression for the mapping, you can plug the real and imaginary co-ordinates into your mSetTest function and get a value representing whether the co-ordinate is part of The Set. The library function, cvSet2D (cxcore.h), is usually used to set the color of a pixel in an image. You can use this function to color the whole set.
// an example of black/white color scheme
CvScalar sca;
if the coordinate of pixel( i,j) in complex plane is in the Mandelbrot Set {
sca.val[0]=0; sca.val[1]=0; sca.val[2]=0; cvSet2D(pImg,j,i,sca);
} else {
sca.val[0]=255; sca.val[1]=255; sca.val[2]=255; cvSet2D(pImg,j,i,sca);
}
//blue //green //red
Task 4
You’re now going to start making the program interactive, and you’ll start by allowing the user to get the coordinate by repeatedly clicking the mouse on a point
11
Project of C/C++ 2011 Spring, CS, HIT
that is to become the new origin. Remember, we are currently drawing the set centered on the imaginary number (0 + 0i) … in this task you are allowing the user to choose a new point to centre our drawing around.
At first, you have to assign the callback function on_mouse which will be used to deal with the mouse message from our system. In you program, you don’t have to call on_mouse function, and it will be called by the system.
The on_mouse function will wait for the user to click on a point in the drawing window, and will get the position of this point via x and y in pixel window co-ordinates (i.e. NOT in Mandelbrot Set co-ordinates.)
You have to work out what these values mean in Mandelbrot co-ordinates (this is the inverse of the calculation you did previously to turn screen co-ordinates into Set co-ordinates of course!) and use these new values as the new centre of your plot.
cvSetMouseCallback(“mandelbrot”,on_mouse,0); // put where you want to deal with the //mouse event
void on_mouse( int event, int x, int y, int flags,void* param) {
// you have to implement this function switch (event)
{
case CV_EVENT_RBUTTONDBLCLK:
} }
break;
case CV_EVENT_LBUTTONDOWN:
break;
case CV_EVENT_MOUSEMOVE:
if((flags & CV_EVENT_FLAG_LBUTTON)) {}
break;
Task 5
The final bit of interactivity is to allow the user to zoom in on their favorite part of the Set to see that it really is a fractal that goes on forever and ever (or, at least, until our floating point calculations run out of precision.). It allows the user to draw a square box on the Mandelbrot Set, and returns the top left and bottom right co-ordinates of that box.
Your program works by allowing the user to click a point on the screen that is to be one corner of the square and then—whilst holding down the mouse button—to drag out a box. The function returns when the user lets go of the mouse button.
Remember, all that you are doing is selecting a new origin around which to centre your calls to mSetTest, and now also changing the distance between the values you are passing to mSetTest.
12
Project of C/C++ 2011 Spring, CS, HIT
*Task 6(optional)
You’ll notice that some parts of the Set aren’t actually very interesting… there are big smooth areas, that if you zoom into them turn into… well, big smooth areas. The real beauty of the set lies in the lovely intricate crinkly bits around the edges.
Automate the process of finding these more chaotic (crinkly!) bits, and ‘steer’ your zoom function so that it continuously selects an interesting area and zooms slowly into it. Very hypnotic. Also, very hard.
Marking Scheme
The project is an optional one for you. Once you decide to do it, please do it entirely (Task 6 may not included). And we may have your code or demo being showed in due course and proper manner.
This lab is designed to get tougher as it goes on. There is very little room for ‘subjectivity’ in the marking scheme… things either work (in which case you get a mark), or they don’t (in which case you don’t get a mark!)
The one subjective mark is given for ‘good coding style’. Any decent, consistent, clear coding style is fine, with any good choices of variable names – we are not looking to penalized you unnecessarily here (though you wont get this mark if your code is a mess, if any of the variables have dreadful names, or if you have made something work but in a completely mad way.)
The following mark rules are for your information.
100 marks for must-do + no marks for optional
10 marks for getting the mSetTest function to work [Task 2]
10 marks for linking with the library and getting a window to appear [Task 3-A] 10 marks for drawing the Set in white and black [Task 3-B]
15 marks for coloring the Mandelbrot Set [Task 3-C]
15 marks for response the mouse event [Task 4-A]
marks for redrawing the Set with the new centre designated by mouse [Task 4-B] 15 marks for getting the zoom to work [Task 5]
10 marks for decent coding style, good variable names etc.
optional: the automated interesting-area zoom function [Task 6]
Acknowledgment:
The original script of the project is kindly provided by Dr. Aphrodite Galata, Manchester, U.K.
13