 |
| A zoom in the Mandelbulb fractal |
The idea
The idea that originated the wave of Mandelbulb fractal was that of trying, once more, to extend the definition of the classic Mandelbrot set to three dimensions. Quaterinons, hypercomplex numbers and all other sort of (often inconsistent) algebras had been used as an attempt to generate the all-time-dream of a true 3D Mandelbrot set. Up to this point the closest thing to a 3D Mandelbrot set were very symetric extrapolations of the 2D set, like the once I did myself back in 2001.
However the construction of the Mandelbulb has a different approach. Instead of thinking in algebra and in the usual iteration formula w->w^2+c, this time we think on the classic 2D set as the result of a geometric process of iterating points by squaring their distance to the origin, rotating them by an angle equal to the current angle with the positive x axis, and then translating by c. Now it's possible to try to extend this geometric process to three dimensions regardless of any algebraic correctness. And in fact the Mandelbulb formula is completelly wrong and incoherent in terms of a dynamic system, but it produces beautifull images, which is in itslef a sign that, perhaps, algebraic correctness is not what we always want. In fact, most of the proto-Mandelbulb images shown in the forums were produced with buggy code that implemented the geometric transformations completelly wrong, with incorrect derivatives calculation for the distance estimation, etc, yet the images were looking just right. The lesson learnt with this is that there is something about duplicating lengths and angles that for some reason makes sense, not only regardless of algebraic correctness of meaning, but even regardless of any geometric interpretation (remember wrong geometric formulas led to nice images anyway). This idea is reinforced by the fact that when breaking the rule of equal length exponentiation to equal angle multiplication, then the images do indeed (finally!) look wrong.
The algorithmBecause in the end any formula that follows the equal-length-exponentiation-to-equal-angle-multiplication will do it, I will choose the one that is geometrically correct. If you are unsure the code you read in other forums or websites is correct, follow the one here, I ensure the one here it is correct. The basic formulation of the Mandelbulb is derived from extracting the polar coordinates of a 3d point and doubling its angles and squaring it's length. The idea can be generalized to other numbers like three, or the more popular case of eight. We will in fact choose the arbitrary value of eight, because for higher powers the asymptotic behaviour of the formulas tends to produce more symmetric shapes (which is another sign that indeed the starting point of multiplying these arbitrary angles is wrong). So, lets call w to our 3D point, then choose eigth as out Mandelbrot power, and so multiply the polar angles of our 3D point by eigth and expand it's modulo by a power of eigth: // extract polar coordinates float wr = sqrt(dot(w,w)); float wo = acos(w.y/wr); float wi = atan(w.x,w.z); // scale and rotate the point wr = pow( wr, 8.0 );" wo = wo * 8.0;" wi = wi * 8.0;" // convert back to cartesian coordinates w.x = sin(wo)*sin(wi); w.y = cos(wo); w.z = sin(wo)*cos(wi); We only have to add c to w now and iterate it in the regular way. You can pretty much then take the code and insert it in your favourite raymarching engine (like this very old but simple one), unless you care about rendering speed. |
 |
Optimizing
The above code is mathematically correct, and even fast enough to run in realtime in modern hardware (as shown thru this website, or in the Shader Toy section). However, one can still gain a 5x speed factor when rendering in some cases (like when rendering in the CPU in pure C or C++).
The first thing to do is, as usual, to get rid of all those trascendental functions (the trigonometric ones) which not only are slow but also introducee unnecesary errors in the computations. By using the basic trigonometric identities of the cosinus and sinus of a doubled angle, one can replace the 8-times angle computations by repeatedly applying the identities (three times). The result is a polynomial which has to trigonometric functions and runs much faster:
float x = w.x; float x2 = x*x; float x4 = x2*x2; float y = w.y; float y2 = y*y; float y4 = y2*y2; float z = w.z; float z2 = z*z; float z4 = z2*z2; float k3 = x2 + z2; float k2 = inversesqrt( k3*k3*k3*k3*k3*k3*k3 ); float k1 = x4 + y4 + z4 - 6.0*y2*z2 - 6.0*x2*y2 + 2.0*z2*x2; float k4 = x2 - y2 + z2; w.x = 64.0*x*y*z*(x2-z2)*k4*(x4-6.0*x2*z2+z4)*k1*k2; w.y = -16.0*y2*k3*k4*k4 + k1*k1; w.z = -8.0*y*k4*(x4*x4 - 28.0*x4*x2*z2 + 70.0*x4*z4 - 28.0*x2*z2*z4 + z4*z4)*k1*k2;
RenderingRendering happens as usual, like for any other 3D fractal. My implementation of the Mandelbulb runs realtime for moderate screen resolutions (say, 800x600) with shadows, but no antialiasing when run on a GPU, and takes a few seconds to render in the CPU. To get this speed one has to use distance based raymarching. In the implementation I used to compute the images and videos in this page there was a subtility, tho. For raymarching with distance fields one needs to be able to compute (an estimation of) the distance from any point to the surface of the set. The well-kown distance formulation for polynomial Julia (and Mandelbrot) sets involves the derivatives of the function being iterated. At the time of implementation I didn't want to extract the Jacobian of the formulas above which would then give me the derivatives, as the formulas where complex and I didn't have any copy of mathematica handy. Others (like in the forum linked above) have proposed formulas for the derivatives which are all wrong, and so far I have not seen any implementation of the Mandelbulb out there using the right code (but many of them believe are doing the right thing cause they simply copied the code from old Quaternion based Julia set renderers). So, I preferred to take another approach while waiting to have a copy of mathematica or some free time to extract the Jacobians manually. I based my solution in the Hubbard-Douady potential theory, which is in turn the principple used to extract the official distance estimation formula. But instead of using the distance estimation formula I applied the raw definition of distance, which is  where G(c) is the H-D potential as is computed as for the regular 2D Mandelbrot set or any other polynomial:  with p being the power of out set (eight in out case, and two for the classic Mandelbrit set). G'(c) is of course the gradient of the potential, and can be approximated by central differences as usual, and that canbe used as shading normal once the intersection of the ray with the set is found. |
|
Experiments
The first video is a close zoom to the surface of the fractal object, and the second one is a morphing of the associated Julia sets of the Mandelbulb set.
 Click image to enlarge |