More efficient generation of Mandelbulbs in voxel environments

[Mator]

Hi there!

You may have seen my introduction topic, if not you can check that out here.

I currently have a program written in javascript which generates Menger Sponges and Mandelbulbs in Minecraft, a voxel-based game.  This script currently generates these objects while the game/server is running in real-time.  I'm eventually going to work my way up to generating other fractals as well.

Yesterday I optimized the Menger Sponge generation code by using a recursive function with some cleverly placed continue;'s, so as to make generation more efficient.  After this optimization the script could generate an 81x81x81 (level 4) Menger Sponge in under 3 seconds.  I felt pretty good about this, as I had come up with the generation code and math all by myself (and it was a great improvement over giant mess of an if statement I was using previously).

So now I'm trying to optimize my Mandelbulb generation code.  I'd like to be able to tell if a block isn't going to reach the cutoff value of 1024 before going through all 21 iterations of applying the vector formula and all that jazz, but I'm not sure how to do this.  I'm pretty certain I should be able to tell by using spherical coordinates to exclude all locations in a certain range, but I'm not entirely certain how to do this safely without accidentally clipping the Mandelbulb in some way.

Perhaps I could just shift the origin to the center of the Mandelbulb and close out of the loop if r <= some value?  r being the radius in spherical coordiantes: r = sqrt(x^2 + y^2 + z^2)?

Edit: Implemented this, I estimate a 16% reduction in execution time.  I've now managed to generate 150x150x150 Mandelbulbs.  Pretty huge!  (295259 blocks)

I'd still like to optimize things some more if possible, suggestions are welcome!   (code updated)

Edit 2: I've been thinking a bit more about optimization of Mandelbulb generation.  My previous optimization only skipped 16% of the empty blocks.  I think, by defining a more complex shape than a sphere using spherical coordinates, I can skip a larger percentage of the empty blocks and optimize things further.  I expect that I can skip more like 60% of the air blocks.  I also realize that by using these min/max radius values I'm causing the generation to lose some generality; viz., changing the formulas could lead to clipping.  I might be able to find a way to use the formulas to generate the min/max radius values if I can identify the theta/phi location of the minimum/maximum radius blocks generally.  Until then I'll probably just have a secondary mandelbulb generator which won't use these values.

Code:

// set global variables
var tfMult = 2.0;
var tfSub = 1.0;
var cutoff = 1024;
var itMin = 5;
var itMax = 21;
var amplitude = 8;
var thetaAdj = 0.000001;
var phiAdj = 0.0000001;
var radiusMin = 0.32 * d;
var radiusMax = 0.56 * d;

// making a mandelbulb
function bulb() {
    var volume = 0;
    for (var x = 0; x < d; x++) {
        var xc = tf(x);
        for (var y = 0; y < d; y++) {
            var yc = tf(y);
            for (var z = 0; z < d; z++) {
                /*
                * Continue if not in range radiusMin -> radiusMax.
                * This saves a lot of time by not performing 21 iterations
                * on locations that clearly won't have blocks. */
                var r = rfc(x,y,z);
                if (r < radiusMin)
                    continue;
                if (r > radiusMax)
                    continue;
               
                var zc = tf(z);
               
                var iterations = -1;
                var C = new Vector(xc, yc, zc);
                var Z = new Vector(0, 0, 0);
               
                // iterate over vectors
                while ((mag(Z) <= cutoff) && (iterations < itMax)) {
                    Z = add(formula(Z, amplitude), C);
                    iterations++;
                }
               
                // place block if cutoff reached in itMin -> itMax range
                if ((iterations >= itMin) && (iterations < itMax)) {
                    if (r < radiusMin)
                        radiusMin = r;
                    if (r > radiusMax)
                        radiusMax = r;
                    volume++;
                    if (volume % 10000 == 0)
                        context.print("    "+volume+" blocks placed... \n");
                    var vec = new Vector(
                        region.getMinimumPoint().getX() + x,
                        region.getMaximumPoint().getY() + y,
                        region.getMinimumPoint().getZ() + z);
                    if (bt.equals("default"))
                        blocktype = context.getBlock(defaultBlockList[iterations - itMin])
                    else if (bt.equals("wool"))
                        blocktype = context.getBlock(woolBlockList[iterations - itMin])
                    else if (bt.equals("plant"))
                        blocktype = context.getBlock(plantBlockList[iterations - itMin])
                    else if (bt.equals("nether"))
                        blocktype = context.getBlock(netherBlockList[iterations - itMin])
                    else if (bt.equals("jungle"))
                        blocktype = context.getBlock(jungleBlockList[iterations - itMin])
                    else if (bt.equals("mine"))
                        blocktype = context.getBlock(mineBlockList[iterations - itMin])
                    else if (bt.equals("desert"))
                        blocktype = context.getBlock(desertBlockList[iterations - itMin])
                    else if (bt.equals("city"))
                        blocktype = context.getBlock(cityBlockList[iterations - itMin])
                    else if (bt.equals("cabin"))
                        blocktype = context.getBlock(cabinBlockList[iterations - itMin])
                    else if (bt.equals("test")) {
                        blocktype = context.getBlock(testBlockList[iterations - itMin]);
                        //test[iterations - itMin]++;
                    }
                    blocks.setBlock(vec, blocktype);
                }
            }
        }
    }
    context.print("Placed "+volume+" blocks!");
}

// the coordinate transform function
// returns a number between -1 and 1
function tf(v) {
    return (tfMult * v)/d - tfSub;
}

// vector magnitude function
// returns the length of the vector
function mag(vec) {
    var x = vec.getX(), y = vec.getY(), z = vec.getZ();
    return Math.sqrt(x * x + y * y + z * z);
}

// vector addition function
// adds two vectors together
function add(v1, v2) {
    var x1 = v1.getX(), y1 = v1.getY(), z1 = v1.getZ();
    var x2 = v2.getX(), y2 = v2.getY(), z2 = v2.getZ();
    var output = new Vector(
        x1+x2,
        y1+y2,
        z1+z2);
   
    return output;
}

// vector growth formula
// if current vector magnitude > 1 and t is not too close to a multiple of (pi/2)
// this will give a vector of greater magnitude than the vector put in
// n, which is param, is the rate at which vector magnitude grows.
function formula(vec, n) {
    var t = theta(vec);
    var p = phi(vec);
    var k = Math.pow(mag(vec), n);
    var output = new Vector(
        k*Math.sin(n*t)*Math.cos(n*p),
        k*Math.sin(n*t)*Math.sin(n*p),
        k*Math.cos(n*t));
   
    return output;
}

// theta vector value (arccos)
function theta(vec) {
    var x = vec.getX(), y = vec.getY(), z = vec.getZ();
    return (Math.acos(z/(mag(vec) + thetaAdj)));
}

// phi vector value (arctan)
function phi(vec) {
    var x = vec.getX(), y = vec.getY(), z = vec.getZ();
    return (Math.atan(y/(x + phiAdj)));
}

// distance from center
// moves vector origin to center and calculates its magnitude
function rfc(x, y, z) {
    var vec = new Vector(x - (d/2), y - (d/2), z - (d/2));
    return mag(vec);
}

Here's a few pictures of some of the fractals I've generated so far. 

[Mator]

Sorry for the double post, but I feel this deserves a separate post.

In looking into Mandelbulb generation efficiency I've continued to learn more about how they operate.  Essentially, the closer you get to the center of the Mandelbulb the more "chaotic" things become.  I've found that the first 2 iterations can't yield a value greater than my cutoff of 1024, and that the third iteration accounts for 33.176% of the volume of the cube.  successive iterations account for smaller and smaller amounts of volume.

Iterations proceed from the border of the cube inward.  e.g.


gray = iteration 3


black = iteration 4


red = iteration 5


orange = iteration 6


yellow = iteration 7

Like peeling back layers of an onion!

(This is with a 100x100x100 Mandelbulb in Minecraft with my specified settings, so doesn't apply generally.)

Code:

the third iteration yields
331176/1000000 = 33.1176% of the volume

fourth iteration yields
212260/1000000 = 21.2260% of the volume

fifth iteration yields
59761/1000000 = 5.9761% of the volume

sixth iteration yields
27566/1000000 = 2.7566% of the volume

seventh iteration yields
14745/1000000 = 1.4745% of the volume

eight iteration yields
10010/1000000 = 1.0010% of the volume

ninth iteration yields
6521/1000000 = 0.6521% of the volume

tenth iteration yields
5447/1000000 = 0.5447% of the volume

eleventh iteration yields
4015/1000000 = 0.4015% of the volume

twelfth iteration yields
3294/1000000 = 0.3294% of the volume

thirteenth iteration yields
2587/1000000 = 0.2587% of the volume

fourteenth iteration yields
2540/1000000 = 0.2540% of the volume

fifteenth iteration yields
2058/1000000 = 0.2058% of the volume

the tenth-fifteenth iterations account for only 1.99% of the volume (higher iterations would account for an even lower % of the total volume.

using these data I've learned that the majority of the execution time is being spent on the interior of the Mandelbulb where all 20 iterations are being used up without exceeding the cutoff value.  Using a 100x100x100 Mandelbulb I roughly calculated the minimum radius at 100 different theta-phi ranges in the Mandelbulb to range from 0.32*d to 0.37*d (d being the length of a side of the cube containing the Mandelbulb).  Storing these values in an array before execution will allow me to cut down execution time even more than I previously had.  I'll perform some tests to see how much time this actually saves, but this should be nice for basic Mandelbulb generation in Minecraft.