Simple "Mandelbulb" C source code.

Compile with

gcc -o Mandelbulb Mandelbulb.c -lm -O3


/*
 * Mandelbulb
 *
 * Written by PM 2Ring November 2009.
 * Converted from Python to C, 2009.11.21
 *
 */

/*

From http://www.skytopia.com/project/fractal/mandelbulb.html

Similar to the original 2D Mandelbrot , the 3D formula is defined by:

z -> z^n + c

....but where 'z' and 'c' are hypercomplex ('triplex') numbers, representing
Cartesian x, y, and z coordinates. The exponentiation term is defined by:

{x,y,z}^n = r^n {sin(theta*n) * cos(phi*n), sin(theta*n) * sin(phi*n),
cos(theta*n)}
where
r = sqrt(x^2 + y^2 + z^2)
theta = atan2(sqrt(x^2+y^2), z)
phi = atan2(y, x)
n is the order of the 3D Mandelbulb. Use n=8 to find the exact object in this
article.

The addition term in z -> z^n + c is similar to standard complex addition,
and is simply defined by:

{x,y,z}+{a,b,c} = {x+a, y+b, z+c}

The rest of the algorithm is similar to the 2D Mandelbrot!

*/

#include <math.h>
#include <stdio.h>
#include <time.h>

//Compute sin(8t) & cos(8t) given r*sin(t), r*cos(t) & r
void trig8(double *ps, double *pc, double r)
{
    //int i;
    double c = *pc, s = *ps, t;

    if (r > 0.0)
    {
        s /= r;
        c /= r;
        //Iterate formulae for sin(2t) & cos(2t) 3 times
        //to get sin(8t) & cos(8t)
        //for (i=0; i<3; i++)
        //{
            t = 2.0*s*c; c = c*c - s*s; s = t;
            t = 2.0*s*c; c = c*c - s*s; s = t;
            t = 2.0*s*c; c = c*c - s*s; s = t;
        //}

        *ps = s;
        *pc = c;
    }
}

//Determine if voxel (cx, cy, cz) is in the Mandelbulb
int Mandel3D(double cx, double cy, double cz, int iters)
{
    double x = cx, y = cy, z = cz, w, w2, r, r2, r8,
    sth8, cth8, sph8, cph8;
    int i;

    for (i=0; i<iters; i++)
    {
        //Convert (x, y, z) coordinates to polar (r, th, ph)
        w2 = x*x + y*y;
        w = sqrt(w2);
        r2 = w2 + z*z;
        r = sqrt(r2);
        if (r >= 2.0)
            return 0;

        //Raise (r, th, ph) to the 8th power
        sth8 = w; cth8 = z; trig8(&sth8, &cth8, r);
        sph8 = y; cph8 = x; trig8(&sph8, &cph8, w);

        //(r, th, ph) <- c + (r, th, ph)^8
        r8 = r2 * r2 * r2 * r2;
        x = cx + r8 * sth8 * cph8;
        y = cy + r8 * sth8 * sph8;
        z = cz + r8 * cth8;
    }
    return 1;
}

int Mandelbulb(double vx, double vy, double vz, double scale, int side, int
iters, char *fname)
{
    FILE *f;
    unsigned char *row;
    int i, j, k;
    time_t itime;
    double x, y, z;

    if(!(f = fopen(fname, "wb")))
        return 1;

    if(!(row = (unsigned char *) malloc(side)))
        return 1;

    //Write header: voxel grid dimensions in little endian
    {
        unsigned char lo = side & 0xff, hi = side >> 8;

        fprintf(f, "%c%c%c%c%c%c", hi, lo, hi, lo, hi, lo);
    }

    itime = time(NULL);
    for (k=0; k<side; k++)
    {
        //Planes
        z = vz + scale * k;
        for (j=0; j<side; j++)
        {
            //Rows
            y = vy + scale * j;
            for (i=0; i<side; i++)
                //Voxels
                row[i] = Mandel3D(vx + scale * i, y, z, iters) ? 255 : 0;
            fwrite(row, 1, side, f);
        }

        if (k%10 == 0)
        {
            int n = (int) difftime(time(NULL), itime);

            //fputc('.', stderr);
            fprintf(stderr, "\r %2d:%02d:%02d Plane %d of %d ",
                n/3600, (n/60)%60, n%60, k, side);
        }
    }

    free(row);
    fclose(f);

    return 0;
}

int main(int argc, char *argv[])
{
    int side, iters, rc;
    char *fname;
    double r, scale;

    side = argc > 1 ? atoi(argv[1]) : 100;
    fname = argc > 2 ? argv[2] : "Mandelbulb.df3";
    iters = argc > 3 ? atoi(argv[3]) : 7;

    r = 1.1;    //Mandelbulb Set "radius"
    //r = 1.05 * (1.0 + 2.0 / side);
    scale = 2.0*r / (side - 1);

    rc = Mandelbulb(-r, -r, -r, scale, side, iters, fname);

    fputc('\n', stderr);
    return rc;
}

Post a reply to this message