Mandelbrot image generator

Question

I have written a C program to generate a PPM (Portable Pixmap) image of the Mandelbrot set. The program implements many things in C I am fairly unfamiliar with (Structures, error handling, use of new libraries (complex.h), and files (I'm completely unfamiliar with image files)) and the purpose of writing this program is primarily to test my knowledge.

If possible, could someone tell me things I've done well here, and things I've done badly, as well as a brief summary of how to do things better in the future? I'd also quite like to know any ways I could reduce execution time. All feedback is greatly appreciated.

I've added comments to give a brief summary of each function and their purposes, as well as make some possibly unclear areas clearer.

//mandelbrot.c - generates a .PPM (Portable Pixmap format) file of the Mandelbrot set with shading
//Still to add: Implement a better file format, optimise to reduce time, 
#include "stdio.h"
#include "complex.h"
#include "math.h"
#define MAX_TESTS 650

int error(const char *message);
struct colour mandelbrot_test(double complex c);
struct colour rgb_gen(int iterations, double complex num);
struct dimensions dim_gen(int height);

struct dimensions
{
    double hinc;
    double winc;
    unsigned int height;
    unsigned int width;
};

struct colour
{
    unsigned int red;
    unsigned int green;
    unsigned int blue;
};

int main(int argc, char *argv[])
{
    if (argc < 3)
    {
        return error("Too few args!\nCorrect usage: program_name image_name image_height");
    }

    char *file_name = argv[1];
    FILE *file;
    double i, j;
    double complex num;
    struct colour rgb;
    struct dimensions dim;

    dim.height = atoi(argv[2]);
    file = fopen(file_name, "w");

    if (!file)
    {
        return error("Unable to access file!\n");
    }
    else if (dim.height < 1000) //values under ~1000px cause scaling issues with the file
    {
        return error("Image cannot be less than 1000 px in height");
    }

    dim = dim_gen(dim.height);
    fprintf(file, "P3\n%d %d\n255\n", dim.width + 1, dim.height); //Magic number, image dimensions and max RGB value

    for (j = -1.0; j <= 1.0; j += dim.hinc)
    {
        for (i = -2.0; i <= 0.5; i += dim.winc)
        {
            num = i + j * I; //Generate the complex number
            rgb = mandelbrot_test(num); //Generate an RGB value for that number
            fprintf(file, "%d %d %d ", rgb.red, rgb.green, rgb.blue); //Print it to the file
        }
        fprintf(file, "\n");
    }

    fclose(file);

    return 0;
}

struct colour mandelbrot_test(double complex c) //Function to test for divergence of a given number
{
    double complex x = 0;
    int i;
    double r_squared = cimag(c) * cimag(c) + creal(c) * creal(c);

    if (r_squared * (8.0 * r_squared - 3.0) < 3.0/32.0 - creal(c)) //Quick test to see if we can bail out early
    {
        return rgb_gen(MAX_TESTS, x);
    } 

    for (i = 1; i < MAX_TESTS; i++) //Actual testing loop
    {
        x *= x;
        x += c;

        if (cimag(x) * cimag(x) + creal(x) * creal(x) >= 4)
        {
            return rgb_gen(i, x);
        }
    }

    return rgb_gen(MAX_TESTS, x);
}

struct colour rgb_gen(int iterations, double complex num) //Function to generate the RGB values for a given complex number
{
    struct colour rgb;  

    if (iterations == MAX_TESTS)
    {
        rgb.red = 0;
        rgb.green = 0;
        rgb.blue = 0;
    }
    else
    {
        double z = sqrt(creal(num) * creal(num) + cimag(num) * cimag(num));
        int brightness = 256.0 * log2(2.75 + iterations - log2(log2(z))) / log2((double)MAX_TESTS); //this generates shading based on how fast the number diverges
        rgb.red = brightness;
        rgb.green = brightness;
        rgb.blue = 255;
    }
    return rgb;
}

struct dimensions dim_gen(int height) //Function to generate 5:2 scale dimensions and width/height increments based on the user given height
{
    struct dimensions dim;
    dim.height = height;
    dim.width = 5.0 * (float)height / 2.0;
    dim.hinc = 1.0 / (0.5 * (float)dim.height);
    dim.winc = dim.hinc / 2.0;
    return dim;
}

int error(const char *message) //Rudimentary error handling
{
    printf("%s", message);
    return 1;
}

Answer 1

Several things:

1. I compiled your code on my Linux host (Centos 7) using gcc version 4.8.3 via the command line: gcc -std=c99 -pedantic -Wall mandlebrot.c -o mandlebrot -lm and I got the following warning:

   gcc -std=c99 -pedantic -Wall mandlebrot.c -o mandlebrot -lm
   mandlebrot.c: In function ‘main’:
   mandlebrot.c:42:5: warning: implicit declaration of function ‘atoi’ [-Wimplicit-function-declaration]
   dim.height = atoi(argv[2]);
   

You need to add #include <stdlib.h> to your program.

2. System include files should have angle brackets instead of quotes, as others have pointed out.

3. My personal taste is to initialize all variable prior to use, this I would write `FILE* file=NULL;

4. I dislike positional command line arguments (i.e. file name first and then height); you might want to look at the getopt library if you are using Linux\Unix (even if you are using Windows it is still worth a look).

5. You can "leak resources", you never close the file you open in 43 if the height is less than 1000. Probably not a big issue in this code sample, but you should get in the habit of releasing any resources that you acquire. For this program you can do it by rearranging your tests:

   if(dim.height < 1000)
   {
         return error("image cannot be less than 1000px in height");
   }
   
   if(NULL != (file = fopen(file_name, "w")))
   {
        // ... generate mandlebrote image here....
   
        fclose(file);
   }
   else
   {
       // ... handle failure to open file
   }
   

6. As others have mentioned, typedef'ing your structures can make your code more readable, and save you on typing.

7. Assuming you are using a c99 compliant compiler (and seeing as how you are using the // single line comment it appears you are), you can write your for-loops like:

   for (double j = -1.0; j <= 1.0; j += dim.hinc)
   {
       for (double i = -2.0; i <= 0.5; i += dim.winc)
       {
           num = i + j * I; 
           rgb = mandelbrot_test(num); 
           fprintf(file, "%d %d %d ", rgb.red, rgb.green, rgb.blue); 
        }
        fprintf(file, "\n");
   }
   

This limits the scope of your index variables to just the loops.

Answer 2

• A P3 magic number indicates a plain PPM, which has a serious limitation:

  No line should be longer than 70 characters.

  Your code obviously violates it. I recommend to make it into raw PPM (magic number P6, and binary representation of colors).

• complex implements cabs(). Use it instead of spelling out cimag(c) * cimag(c) + creal(c) * creal(c).

• mandelbrot_test calculates number of iterations and calculates the resulting color. These activities are unrelated and shall be done in separate functions.

• The code is overcommented to my taste. The only comment I think has a right to exist is about a quick convergence test, and I believe it should be expanded. It took me a while to understand why this test works.

Answer 3

I'm not familiar with the algorithm you are implementing, so I can only comment on the C code.

---

#include "stdio.h"
#include "complex.h"
#include "math.h"

Those are standard header files, so the correct way of referencing them is between <>, not quotes. E.g.:

#include <stdio.h>

When you reference a file between <> you tell the compiler to look first for that file in the standard include paths. Using "" tells the compiler to look for files in the project directory, that's why you use quotes for local project files.

---

It is usually better to avoid function prototypes if you can and declare your helper functions before main(). The issue with function prototypes is that they produce code duplication and thus maintenance overhead.

---

You could typedef your structs to avoid having to supply the struct tag whenever you declare a colour or dimension. I tend to prefer that as I don't see a point in repeating yourself every time, but YMMV.

typedef struct
{
    double hinc;
    double winc;
    unsigned int height;
    unsigned int width;
} dimensions;

typedef struct
{
    unsigned int red;
    unsigned int green;
    unsigned int blue;
} colour;

// 'colour' and 'dimensions' are now first class names and 
// can be used directly to declare instances.

---

When if clauses return at the end, you should avoid chaining with else. It is cleaner to just return. E.g.:

if (!file)
{
    return error("Unable to access file!\n");
}
else if (dim.height < 1000) //values under ~1000px cause scaling issues with the file
{
    return error("Image cannot be less than 1000 px in height");
}

You don't need that else since both return. It could be broken into two independent ifs. In this case, this is a minor nitpicking, but if you had a longer chain with every entry returning, then it would be a nice improvement to simplify it.

---

I personally find it more readable to align assignment blocks like the following by the = sign:

struct dimensions dim;
dim.height = height;
dim.width  = 5.0 * (float)height / 2.0;
dim.hinc   = 1.0 / (0.5 * (float)dim.height);
dim.winc   = dim.hinc / 2.0;
//         ^
//    aligned here

BTW, those casts above should be doubles, since you are not using floats anywhere else. However, since one side of the expressions is already a double, the casts are superfluous.

---

Your error() function should probably print to stderr, which is the standard error output. stdout (printf) is assumed to print normal program output.

---

If memory and storage ever becomes a concern to you, you can store unsigned chars (bytes) in your colour struct, since the RGB values are between 0 and 255.

Answer 4

Brightness calculation

Regarding the following code:

double z = sqrt(creal(num) * creal(num) + cimag(num) * cimag(num));
int brightness = 256.0 * log2(2.75 + iterations - log2(log2(z))) / log2((double)MAX_TESTS);

You don't need to calculate creal(num) * creal(num) + cimag(num) * cimag(num), since whenever you call this function, it's within a line or two of already having calculated it (I'm going to call this value r). Just pass the already calculated value in.

You also don't need to use sqrt here, since log(sqrt(x)) == 0.5*log(x). Also, since log(a*b) == log(a) * log(b) and log2(1/2) == -1, you can simplify it further:

int brightness = 256.0 * log2(3.75 + iterations - log2(log2(r))) / log2((double)MAX_TESTS);

With that, you don't have to calculate r or its square root, but I still don't like the log2(log2(r)) part. I'm not sure why it's there (there is a comment on that line, but it doesn't explain where the formula came from). Taking a log of a log makes a number really small really fast -- log2(log2(11,224)) is 3.75. We know r >= 4, and since we're just squaring and adding a small constant, it won't get much higher than 16.

Unless I've misunderstood or missed something (in which case, there should be a comment explaining it), this line will be approximately equivalent:

int brightness = 256.0 * log2(2.0 + iterations) / log2((double)MAX_TESTS);

but I think this next one would be slightly better, since we can tell that it changes the multiple of 256 from 0.0 to 1.0 exactly:

int brightness = 256.0 * log2(iterations) / log2((double)MAX_TESTS-1);

Symmetry

Complex numbers have a symmetry due to the impossibility of distinguishing between i and -i (There is no way to show that we put the minus sign on the 'correct' one. They would behave in the exact same way if we renamed i to -i and -i to i.)

So whenever you calculate some pixel's color, you've also calculated the color of another pixel on the opposite side of the real number line. If that pixel also happens to be in the image, you don't need to recalculate it later.

Bailout test

if (r_squared * (8.0 * r_squared - 3.0) < 3.0/32.0 - creal(c))

This is hard to read, and the comment didn't help. If the bailout test could be improved to cover a larger portion of the points in the Mandelbrot set, the code would run faster, since all of those points would otherwise run through the loop the maximum number of times.

Since I don't know how it works, I can't make any suggestions on improving it, but if you want to see how much it covers right now, change the code to set the bailout points to some color that doesn't appear elsewhere in the image. Seeing the parts it doesn't cover may give you an idea of how to improve it, or you may end up satisfied that it is good enough as is.

Answer 5

1. Use perror instead of defining your own error function. It will automatically report system error messages in addition to whatever string you supply, and it prints to stderr instead of stdout.

2. Prefer static global variables to locals in the outer scope of main. You can still implement encapsulation and pass the globals to functions. The dimensions of the image might as well be a global. But rgb really should be local to the inner loop.

3. Likewise, the actual complex math should be in a separate function from main. i, j perhaps should be a complex number in the first place, instead of regenerating num with every inner iteration.

4. The formula c * conj(c) is simpler than writing out a sum of squares, without sacrificing performance.

5. I think x = ( x * x ) + c; would be more readable than x *= x; x += c;. It's more like a formula and less like an algorithm.

6. Prefer strtol over atoi because it validates the input. It's nice to print a usage reminder instead of complaining about the height if the first argument is not a number.

7. Bit of a nitpick, but you can break multi-line strings like so:

   perror("Too few args!\n"
          "Correct usage: program_name image_name image_height");
   

   and

   fprintf(file,
       "P3\n" //Magic number
       "%d %d\n" // image dimensions
       "255\n", // max RGB value
       dim.width + 1, dim.height);
   

   The latter perhaps should be several statements instead.