Plotting the Mandelbrot set efficiently

Question

This is my Mandelbrot set program - it includes Smooth Coloring, perodicity checks, and my approach to biomorphs.

I know of Edge Detection, but I couldn't find a simple explanation that I can understand. Any help appreciated!

How do I optimize the Mandelbrot set plotting algorithm implemented in escapesmooth() and display() further?

#include <GL/glut.h>
#include <iostream>
#include <cmath>
float lerp(float a,float b,float t){return a+t*(b-a);}
const int WIDTH =800;
const int HEIGHT=600;
float r1,r2,g1,g2,b1,b2;
bool biomorph=false;
double xPos=-3.2;
double yPos=-2.0;
double zoom=150;
double startX,startY,start_X,start_Y;
int maxIter =50;
int maxIter2=50;
double zoomFactor=1.0;
float period=0;
bool render=true;
float escapesmooth(double real, double imag){ // Optimized Escape Time Algorithm + smooth colors + perodicity checks
    double r0=real;
    double i0=imag;
    float iter=0.0;
    double tempreal,io,ro;
    while(iter<maxIter&&real*real+imag*imag<=16){
        tempreal=real*real-imag*imag+r0;
        imag=2*real*imag+i0;
        real=tempreal;
        iter++;
        if(real==ro&&imag==io){
            iter=maxIter;
            break;
        }
        //interesting effect
        if(biomorph){
            if(abs(real)>0.75&&abs(imag)>0.75){
                iter=real+imag;
                break;
            }
        }
        period++;
        if(period>20){
            period=0;
            ro=real;
            io=imag;
        }
    }
    if(iter<maxIter){
        float log_zn=log((real*real)+(imag*imag))/2;
        float nu=log(log_zn/log(2))/log(2);
        iter=iter+1.0-nu;
    }
    return iter;
}

void display(){
    if(render){
        for(int x=0;x<WIDTH;x++){
            glBegin(GL_POINTS);
            for(int y=0;y<HEIGHT;y++){
                float iterations=escapesmooth(x/zoom+xPos,y/zoom+yPos);
                glVertex2i(x,y);
                if(iterations==maxIter){
                    glColor3f(0,0,0);
                }
                else{
                    float r1=0.5+sin(floor(iterations)/10)/2;
                    float g1=floor(iterations)/maxIter;
                    float b1=1-floor(iterations)/maxIter;

                    float r2=0.5+sin((floor(iterations)+1)/10)/2;
                    float g2=  (floor(iterations)+1)/maxIter;
                    float b2=1-(floor(iterations)+1)/maxIter;
                    /*Sunburst:*//*
                    float r1=floor(iterations)/32;
                    float r2=floor(iterations+1)/32;
                    float g1=floor(iterations)/64;
                    float g2=floor(iterations+1)/64;
                    float b1=floor(iterations)/128;
                    float b2=floor(iterations+1)/128;
                    /*test*//*
                    float r1=floor(iterations)/32;
                    float r2=floor(iterations+1)/32;
                    float g1=floor(iterations)/64;
                    float g2=floor(iterations+1)/64;
                    float b1=floor(iterations)/128;
                    float b2=floor(iterations+1)/128;*/

                    float t=fmod(iterations,1);
                    glColor3f(lerp(r1,r2,t),lerp(g1,g2,t),lerp(b1,b2,t));
                }
            }
            glEnd();
            glFlush();}
        }
    //maxIter*=1.2;
}
void reshape(int w,int h){glutReshapeWindow(WIDTH,HEIGHT);}
bool isDragging = false;
int lastX,lastY;

void mouse(int button, int state, int x, int y){
    double mouseX=static_cast<double>(x)/WIDTH;
    double mouseY=static_cast<double>(y)/HEIGHT;
    if(button==GLUT_LEFT_BUTTON&&state==GLUT_DOWN){
        isDragging=true;
        lastX=x;
        lastY=y;
    } else if(button==GLUT_LEFT_BUTTON&&state==GLUT_UP){
        isDragging=false;
    }
}
void motion(int x, int y) {
    if (isDragging) {
        double deltaX=(x-lastX)/zoom;
        double deltaY=(y-lastY)/zoom;
        xPos-=deltaX;
        yPos+=deltaY;
        lastX=x;
        lastY=y;
        maxIter=maxIter2;
        glutPostRedisplay();
    }
}
void key(unsigned char key, int x, int y) {
    switch (key){
    case 'q': {
        zoom *= 1.1;
        maxIter = maxIter2;
        render = true;
        glutPostRedisplay();
        break;
    }

    case 'e': {
        zoom *= 0.9;
        maxIter = maxIter2;
        render = true;
        glutPostRedisplay();
        break;
    }
        case '=': {maxIter2 += 1;maxIter = maxIter2;render = true;glutPostRedisplay();break;}
        case '+': {maxIter2 += 1;maxIter = maxIter2;render = true;glutPostRedisplay();break;}
        case '-': {maxIter2 -= 1;maxIter = maxIter2;render = true;glutPostRedisplay();break;}
        case 'b': {if(biomorph){biomorph=false;}else{biomorph=true;}glutPostRedisplay();break;}
        case 'r': {double xPos=-3.2;double yPos=-2.0;double zoom=150;maxIter = 50;render = true;glutPostRedisplay();break;}
    }
}

int main(int argc, char** argv){
    glutInit(&argc,argv);
    glutInitDisplayMode(GLUT_SINGLE | GLUT_RGB);
    glutInitWindowSize(WIDTH,HEIGHT);
    glutCreateWindow("mandel badandel");

    glMatrixMode(GL_PROJECTION);
    glLoadIdentity();
    gluOrtho2D(0, WIDTH, 0, HEIGHT);
    glMatrixMode(GL_MODELVIEW);
    glLoadIdentity();
    glClearColor(0.0,0.0,0.0,1.0);

    glutDisplayFunc(display);
    glutReshapeFunc(reshape);
    glutMouseFunc(mouse);
    glutMotionFunc(motion);
    glutKeyboardFunc(key);

    glutMainLoop();
    return 0;
}

Answer 1

Unnecessary use of OpenGL

Using OpenGL to draw individual pixels is the completely wrong thing to do. Apart from requiring OpenGL support which might not be present on all systems, this has a huge overhead.

Instead of using the rather low-level GLUT library, I recommend that you use a library like SDL to draw to a 2D image (in SDL, that is done using a Renderer) and to handle user input.

float or double?

I see both float and double being used in the same expressions. Usually you want to pick one type and stick with it, otherwise unnecessary conversions between them will happen.

Make better use of C++

As Toby Speight already mentioned, use std::complex<> types to store complex numbers. Most math functions will also directly work on those types, greatly simplifying your equations.

Since C++20 there is std::lerp(), so you don't have to write your own.

Use a code formatting tool

Toby Speight already mentioned the high code density. It would be much more readable if you had consistent whitespace around operators, empty lines between functions, some more empty lines to separate blocks of code and to make the structure more clear.

Instead of doing that all manually, consider using a code formatting tool to automate this for you. Some IDEs and code editors might also have a formatting function built-in.

Answer 2

First impression is that the code is very dense. Reading would be a bit easier if you give the code room to breathe!

Variable names in particular are short and uninformative (why do we have both i0 and io in the same scope, for example). Some of the longest names are real and imag, which we wouldn't even need if we included <complex> and used a std::complex type.

Though we include the C++ <cmath> header, we use the C names log, abs, etc. to call its functions. That's not portable - write std::log, std::abs, etc. in full.

Can you do anything to reduce the number of global variables? Shared state always makes it harder to reason about the behaviour of functions in a program.

Some code is unnecessarily convoluted - e.g. if(biomorph){biomorph=false;}else{biomorph=true;} where simple biomorph = !biomorph; would be clearer.

Answer 3

General Observations

I would not want to maintain this code or delegate this code to someone else to maintain.

It is difficult to read and even more difficult to modify.

Before trying to optimize the code the code should be easy to read, write and maintain.

Commented out code such as the Sunburst code in the `display indicates that the code is not ready for review, nor is it ready to be maintained.

Avoid Global Variables

It is very difficult to read, write, debug and maintain programs that use global variables. Global variables can be modified by any function within the program and therefore require each function to be examined before making changes in the code. In C and C++ global variables impact the namespace and they can cause linking errors if they are defined in multiple files. The answers in this stackoverflow question provide a fuller explanation.

It is important to note that not all of the global variables are being used. All of the variables on the following lines are not used in the code.

float r1, r2, g1, g2, b1, b2;
double startX, startY, start_X, start_Y;

It is possible that some of the global variables are of the wrong type, for instance

double xPos = -3.2;

should probably be declared as

GLdouble xPos = -3.2;

for portability.

Declare the Variables as Needed

In the original version of C back in the 1970s and 1980s variables had to be declared at the top of the function. That is no longer the case, and a recommended programming practice to declare the variable as needed. In C the language doesn't provide a default initialization of the variable so variables should be initialized as part of the declaration. For readability and maintainability each variable should be declared and initialized on its own line.

Vertical Spacing

While there are some blank lines in the code, the first section of code is almost unreadable due to the lack of vertical spacing.

Rather than cram multiple statements onto a single line the code would be more maintainable if each statement was on a separate line.

    case '=': {maxIter2 += 1; maxIter = maxIter2; render = true; glutPostRedisplay(); break; }

would be much easier to read and maintain if it looked like this:

    case '=':
        {
            maxIter2 += 1;
            maxIter = maxIter2;
            render = true;
            glutPostRedisplay();
            break;
        }

DRY Code

There is a programming principle called the Don't Repeat Yourself Principle sometimes referred to as DRY code. If you find yourself repeating the same code multiple times it is better to encapsulate it in a function. If it is possible to loop through the code that can reduce repetition as well.

In the switch statement in the key() function there is duplication of code:

    case '=': {maxIter2 += 1; maxIter = maxIter2; render = true; glutPostRedisplay(); break; }
    case '+': {maxIter2 += 1; maxIter = maxIter2; render = true; glutPostRedisplay(); break; }

In this case the code duplication is not necessary, because you can put both of the cases together:

    case '=':
    case '+':
        {
            maxIter2 += 1;
            maxIter = maxIter2;
            render = true;
            glutPostRedisplay();
            break;
        }

The could be more reduction of code duplication if the maxIter2 assignment was removed from the block.

    case '=':
    case '+':
        maxIter2 += 1;
    case '-':
        maxIter2 -= 1;
        {
            maxIter = maxIter2;
            render = true;
            glutPostRedisplay();
            break;
        }

It is possible that there could be even more reduction of duplicate code in this function.

Answer 4

Since you want to do this efficiently - measure how long each OpenGL call takes. I wouldn't be surprised at all if that takes most of the time. Find out how to store various pixel values int a buffer and draw them with one call.

Measure how much time you actually save by checking that the starting value repeats. I wouldn't be surprised if the amount is negative.

Figure out where the main black blob is and find a formula to check whether a value is inside it. You can also check for the next black circle. If you think that takes too long do the check only if the previous pointer hit maxIter. So if you check points horizontally, you do this check after the first black point until you leave the black area again.

Don't check for escape on every iteration. Do say 8 iterations, then check. If you have escaped, redo the last 8 iterations to find exactly when you escaped. That saves lots of checks and lots of looping.

Use fma() whenever possible. It is faster and more precise (but not 100% identical, so the compiler doesn't do it automatically). fma(a, b, c) = a * b + c with a single rounding, and usually takes the same time as the multiplication on its own.

Check which is faster - fabs(x) > 0.75 or x * x > 0.75 * 0.75.

Use vector registers to calculate multiple values at the same time. And then try doing two or for sets of points at the same time.

float vs double: If you use vector registers then float is faster BUT float isn't every precise. In areas with 200 or 500 iterations, your images will suffer badly.