The Mandelbulb fractal has been the fractal of the year, yet the hype around it lasted for no more than a month. I guess taht gives an idea of how out of fashion fractals
are today. Nevertheless it's an interesting fractal that might encourages the most brave fractal resarchers left to further explore this branch of maths and computer
graphics. I don't consider myself a researcher in the field of fractals as I do recognize my lack of proper academic backgroud, depth of understanding of the field and
lack of dedication unlike others who know little and do bad yet they pretend being experts. Nevertheless, just as any other fractal hobbist that deserves such a geeky
title, I did my experiments on rendering Mandelbulb (and other further math investigations) when the innertia of the hype moment told me to do so. I must say it was a fun
time and I truly thank those who started the wave in the first place, that was really fun.
|A zoom in the Mandelbulb fractal|
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
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.
Because 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:
float wr = sqrt(dot(w,w));
float wo = acos(w.y/wr);
float wi = atan(w.x,w.z);
wr = pow( wr, 8.0 );"
wo = wo * 8.0;"
wi = wi * 8.0;"
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.
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 no 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;
Rendering 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.
|Orbit trap for coloring. Note that thare are no lighting computations of any sort going on|
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.