Mandelbrot fractal

Question

This is a code that outputs a .ppm file with the mandelbrot fractal.

How can I optimize this?

#include<bits/stdc++.h>
using namespace std;

int findMandelbrot(double cr, double ci, int max_iterations)
{
    int i = 0;
    double zr = 0.0, zi = 0.0;
    while (i < max_iterations && zr * zr + zi * zi < 4.0)
    {
        double temp = zr * zr - zi * zi + cr;
        zi = 2.0 * zr * zi + ci;
        zr = temp;
        ++i;
    }
    return i;
}

double mapToReal(int x, int imageWidth, double minR, double maxR)
{
    double range = maxR - minR;
    return x * (range / imageWidth) + minR;
}

double mapToImaginary(int y, int imageHeight, double minI, double maxI)
{
    double range = maxI - minI;
    return y * (range / imageHeight) + minI;
}

int main()
{
    ifstream f("input.txt");
    int imageWidth, imageHeight, maxN;
    double minR, maxR, minI, maxI;

    if (!f)
    {
        cout << "Could not open file!" << endl;
        return 1;
    }

    f >> imageWidth >> imageHeight >> maxN;
    f >> minR >> maxR >> minI >> maxI;

    ofstream g("output_image.ppm");
    g << "P3" << endl;
    g << imageWidth << " " << imageHeight << endl;
    g << "255" << endl;


    double start = clock();

    for (int i = 0; i < imageHeight; i++)
    {
        for (int j = 0; j < imageWidth; j++)
        {
            double cr = mapToReal(j, imageWidth, minR, maxR);
            double ci = mapToImaginary(i, imageHeight, minI, maxI);

            int n = findMandelbrot(cr, ci, maxN);

            int r = ((int)sqrt(n) % 256);
            int gr = (2*n % 256);
            int b = (n % 256);

            g << r << " " << gr << " " << b << " ";
        }
        g << endl;

        if(i == imageHeight / 2) break;
    }

    cout << "Finished!" << endl;

    double stop = clock();

    cout << (stop-start)/CLOCKS_PER_SEC;
    return 0;
}

I go until imageHeight / 2 because in photoshop I can just copy the other half.

I was thinking about logarthmic power but tried something and only works with integers... At the end I am showing the time it takes to run some data. I tried some things from the internet but nothing really worked out or maybe I didn't how to put it right..

this is the output with the following input: 512 512 512 -1.5 0.7 -1.0 1.0

and it took 0.315s to produce this.

Answer 1

Your code produces a "plain" PPM file where all sample values are represented as ASCII decimal numbers (separated by white space). As a consequence, a substantial amount of time is spent with formatting and writing the data.

An alternative is to write a "binary" PPM file. It has the magic "P6" instead of "P3", and all sample values are represented as a single byte. (See PPM Format Specification for the details.)

So you would create the file and write the header with (using a more descriptive variable name for the file stream):

ofstream ppmFile("output_image.ppm", ios::out | ios::binary);
ppmFile << "P6" << endl;
ppmFile << imageWidth << " " << imageHeight << endl;
ppmFile << "255" << endl;

(endl flushes the output file in addition to writing a newline character, which is not necessary here. On the other hand, it does no harm because this part is not performance critical.)

In the inner loop, an RGB triplet is appended with

int n = findMandelbrot(cr, ci, maxN);

uint8_t rgbTriplet[3];
rgbTriplet[0] = (int)sqrt(n);
rgbTriplet[1] = 2 * n;
rgbTriplet[2] = n;
ppmFile.write((char *)rgbTriplet, 3);

and no newlines are written after each row. Truncating of int to the unsigned type uint8_t is well-defined so that the explicit remainder operation % 256 is not needed.

In my test on a MacBook this reduces the running time for the input data "512 512 512 -1.5 0.7 -1.0 1.0" from ≈0.25 seconds to ≈0.14 seconds.

Answer 2

I agree with most of what @user140417 said (though I agree with Deduplicator on the call to pow). So sticking strictly to performance, I have the following recommendations.

Precompute Everything You Can

You are calculating the real and imaginary parts on every iteration of your inner loop in main(). This is costly. Every row will have the same real values and every column will have the same imaginary values. So pre-calculate 1 row of reals and 1 row of imaginaries before starting the main loop.

It's the same situation with your coloring. Since you know maxN, you can have 3 arrays - one for the red channel, one for the green channel, and one for the blue channel. Just pre-calculate all maxN possible values. Then in the inner loop, it's just a table look-up rather than a calculation.

Use SIMD Techniques

Much of your math is computing the same loop for different sets of floating point numbers. This is exactly the kind of thing that SIMD instruction sets like SSE, AVX, etc. were invented for. You can compute 2 or 4 values per instruction, allowing you to significantly improve the speed of your implementation. You can take it even further with…

Multiple Threads

If you break the input area up into either tiles or strips, you can have multiple threads, each working on a different tile or strip at the same time. On multi-core machines this will speed things up in a nearly linear fashion (in my experience) - up to a point. Combine that with SIMD instructions and you can really get things moving!

Use the GPU

If you really want to make things fast, you can write a GPU implementation using OpenGL compute shaders, or OpenCL kernels, or whatever other library you like to use. GPUs generally have hundreds to thousands of cores of floating point processors making this a really good fit.

Answer 3

Readability

Don't include bits/stdc++.h

This is a non-portable header file that pulls in way more than needed, increasing compilation time and binary size.

Don't use using namespace std;

This combined with the above greatly pollutes the global namespace. If you ever decide to use a third party library (for precision math, image drawing, or what have you) there is a high probability you'll run into conflicts.

Use more descriptive variable names

All of the zr, zi, et al. is not only hard to read but easy to make a fatal typo that's hard to debug. You may as well use std::complex or at the least a std::pair. If you later decide to add the ability to zoom or pan your code readability will rapidly diminish if you continue to use those types of variable names.

Rename mapToReal and mapToImaginary

Maybe you are using the mathematically correct definition, but from what I can tell a better name for the function would be "normalize" or "scale", as "map" is too abstract/vague.

Optimization

Factor out the arithmetic inside findMandelbrot

You are doing a total of (I count) 6 multiplications which can be reduced to three if you refactor the code slightly:

std::complex z{0, 0};
std::complex square_temp{0, 0};
while (square_tmp.real + square_tmp.imag <= 4.0) {
    z.imag = std::pow(z.real + z.imag, 2.0) - square_temp.real - square_temp.imag;
    z.imag += c.imag;
    z.real = square_temp.real - square_temp.imag + c.real;
    square_temp.real = std::pow(z.real, 2);
    square_temp.imag = std::pow(z.imag, 2);
}
// Note code untested

Consider using threads/workers

It's possible to gain a speed increase if you pass chunks of the image to render to multiple threads. How you decide to implement this up to you.