Here is a zip containing files you may use to render the mandelbulb in 
POV-Ray.  It contains the following files, mbsdl_* for use with any 
version of POV-Ray, and mb_* for use with a patched POV-Ray.


mbsdl_standalone.pov and mb_standalone.pov are minimal files with no 
dependencies.

mbsdl.inc          This is the include file used by all of the mbsdl_* 
files below.

mbsdl_stage1.pov   mb_stage1.pov   Minimal files
mbsdl_basic.pov    mb_basic.pov    A mandelbulb isosurface colored by a 
'bulb function
mbsdl_powers.pov   mb_powers.pov   Six mandelbulbs, powers 2-7
mbsdl_multi_p.pov  mb_multi_p.pov  Six composite power bulbs
mbsdl_pattern.pov  mb_pattern.pov  Use f_mandelbulb() as a pigment pattern
mbsdl_media.pov    mb_media.pov    Use f_mandelbulb() as a media density

mbsdl_renderall.sh mb_renderall.sh Bash scripts to render them all. (Linux)


There is also a long README.TXT, patches and replacement files to add 
f_mandelbulb() internally, and a patch to extend the expiration of beta 34.

Post a reply to this message

Attachments:
Download 'f_mandelbulb.010.zip' (57 KB)

Here are small anti-aliased images of each of the scene files.

Post a reply to this message

Attachments:
Download 'mb_basic.png' (70 KB) Download 'mb_media.png' (57 KB) Download 'mb_multi_p.png' (50 KB) Download 'mb_pattern.png' (52 KB) Download 'mb_powers.png' (59 KB) Download 'mb_stage1.png' (11 KB) Download 'mb_standalone.png' (36 KB)

Preview of image 'mb_basic.png'

Preview of image 'mb_media.png'

Preview of image 'mb_multi_p.png'

Preview of image 'mb_pattern.png'

Preview of image 'mb_powers.png'

Preview of image 'mb_stage1.png'

Preview of image 'mb_standalone.png'

Here's the update.  [I see the .binaries, so I should be OK.]

I incremented straight to 0.020 since the new function has one more 
parameter than the old (making it incompatible with existing calls), a 
phase parameter between p_phi and i.

I also put in a clunky way to switch between using the SDL and internal 
function, if you compiled it into POV-Ray with the internal patch.  Just 
rename or remove the mbsdl.inc file and the sample scenes will switch to 
using the internal function.  I updated the README.TXT to try to explain 
these and other mostly minor changes.

Also, I added another sample showing the four major variations.  The 
mb_variations.pov scene file is set to show the four power 8 variations, 
but I also attached an image of the power two variants to highlight the 
differences.

Happy 'Bulbing!

~David

Post a reply to this message

Attachments:
Download 'f_mandelbulb.020.zip' (50 KB) Download 'mb_variations.png' (35 KB) Download 'mb2_variations.png' (34 KB)

Preview of image 'mb_variations.png'

Preview of image 'mb2_variations.png'

This looks an awful lot like a 3D version of the tricorn to me.

http://en.wikipedia.org/wiki/Tricorn_%28mathematics%29

It's made with the new function using a phase of 0, p_r=2, p_theta=-2, 
p_phi=2 at ten iterations.  (Note the negative p_theta.)

On another note, Paolo Bonzini recently published a paper describing how 
to express the mandelbulb function variations as a quaternion function 
mapping.  I generalized his expression to work with any scalar power, 
and besides being more flexible, it looks like it will be more 
computationally efficient than the non-integral exponent versions of the 
White, Nylander, et. al., formulae (the ones currently in use).

This quaternion version requires only 1 atan2, 1 cos, 1 sin, 2 sqrt, 1 
fractional pow(), and 3 scalar divisions per iteration (along with a 
bunch of fast multiplies and adds).  Unfortunately, the expression 
becomes extremely cumbersome when attempting to shoehorn it into a 
single SDL function, so no working code yet.

I also plan to flirt with insanity by taking a crack at reworking it for 
full quaternion power exponents.

~David

Post a reply to this message

Attachments:
Download 'mb_iso_tricorn.png' (130 KB)

Preview of image 'mb_iso_tricorn.png'

I previosly wrote:
> I also plan to flirt with insanity by taking a crack at reworking it for 
> full quaternion power exponents.

Well, something cracked, all right.  Since I had to recompile anyway, I 
went ahead and coded a full quaternion power version.  Once I worked out 
the math (or more accurately, it worked me out), the resulting algorithm 
isn't that much more computationally intensive than the old White and 
Naylander formulas, and may even be a touch faster.  (It would still be 
too much to code in SDL unless I can learn how to build functions from 
macros, or write an external script to make the SDL function.)

However, although it returns a reasonable straight power-two mandelbulb, 
the current code doesn't seem to be doing higher powers quite right. 
Regardless, what varying the imaginary parts of the quaternion power 
does is indescribable.

There are much weirder 'bulbs than the one attached, but this (2,0,1,0) 
    quaternion power beast seemed appropriate for the new year.

Wishing ya'll the best!

~David Wagner

Post a reply to this message

Attachments:
Download 'mbq2010_beast_512.png' (136 KB)

Preview of image 'mbq2010_beast_512.png'

Earlier, I blathered:
> However, although it returns a reasonable straight power-two mandelbulb, 
> the current code doesn't seem to be doing higher powers quite right. 
> Regardless, what varying the imaginary parts of the quaternion power 
> does is indescribable.
> 
> There are much weirder 'bulbs than the one attached, but this (2,0,1,0) 
>    quaternion power beast seemed appropriate for the new year.

I found and corrected one error that only applied to using a power with 
an imaginary part.  I attached the 2,0,1,0 mandelquat after rendering it 
with the updated code.

For pure real power mandelbulbs, I'm still only getting the expected 
results for power two, and still don't know why.

~David

Post a reply to this message

Attachments:
Download 'mbq2010_beast_640.png' (199 KB)

Preview of image 'mbq2010_beast_640.png'

waggy babbled on:
> For pure real power mandelbulbs, I'm still only getting the expected 
> results for power two, and still don't know why.

Bonzini updated his paper today, and my pure real power results are 
consistent with his.  The formulas needed to reproduce arbitrary real 
higher-power mandelbulbs using quaternions are quite different than 
originally described (except for power two, after trig reduction), and 
seem to involve even more trig functions than the classic formulation.

http://github.com/bonzini/mbulb/raw/master/mbulb.pdf

There still may or may not be errors in how my code handles the 
imaginary part of quaternion powers.  I think I need to take a step back 
and figure out what it's supposed to do rather than continue to muck 
about with variations.

If there's interest, I'll post what I have so far.  Otherwise, I'll hang 
onto it until I can describe what it does.

~David

Post a reply to this message

waggy <hon### [at] handbasketorg> wrote:
> I previosly wrote:
> > I also plan to flirt with insanity by taking a crack at reworking it for
> > full quaternion power exponents.
>
> Well, something cracked, all right.  Since I had to recompile anyway, I
> went ahead and coded a full quaternion power version.  Once I worked out
> the math (or more accurately, it worked me out), the resulting algorithm
> isn't that much more computationally intensive than the old White and
> Naylander formulas, and may even be a touch faster.  (It would still be
> too much to code in SDL unless I can learn how to build functions from
> macros, or write an external script to make the SDL function.)
....

David,

If you would like to see how hyou can build functions with macros,
then below is a macro that produces a function that is similar to
your recursive (non-quaternion) Mandelbrot function.


// ===== 1 ======= 2 ======= 3 ======= 4 ======= 5 ======= 6 ======= 7

#macro MandelBulbFunction(R_Pwr, R_BO, Ph_Pwr, Ph_Phase, Th_Pwr, Max_Iter)

  #local ln_R_BO = ln(R_BO);
  #local ln_R_Pwr = ln(R_Pwr);
  #local J = Max_Iter - 1;
  #local MB_Fns = array[Max_Iter][2]
  #local I = Max_Iter;
  #while (I > 0)
    #local I = I - 1;
    #local MB_Fns[I][1] =
      function(x, y, z, _r, _rp, _ph, _th) {
        #if (I = J)
          0
        #else
          select(
            _r - R_BO,
            MB_Fns[I+1][0](
              x,
              y,
              z,
              x + _rp*cos(_ph)*cos(_th),
              y + _rp*sin(_ph),
              z + _rp*cos(_ph)*sin(_th)
            ),
            1/(I - ln(ln(_r)/ln_R_BO)/ln_R_Pwr)
          )
        #end // if
      }
    #local MB_Fns[I][0] =
      function(x, y, z, _ix, _iy, _iz) {
        MB_Fns[I][1](
          x,
          y,
          z,
          f_r(_ix, _iy, _iz),
          pow(f_r(_ix, _iy, _iz), R_Pwr),
          Ph_Pwr*(f_ph(_ix, _iy, _iz) + Ph_Phase),
          Th_Pwr*f_th(_ix, _iy, _iz)
        )
      }
  #end // while

  function { MB_Fns[0][0](x, y, z, x, y, z) }

#end // macro MandelBulbFunction


// It can be used like this:
#declare MandelBulbFn = MandelBulbFunction(8, 3, 8, pi/2, 8, 5)

// ===== 1 ======= 2 ======= 3 ======= 4 ======= 5 ======= 6 ======= 7

--
Tor Olav
http://subcube.com

Post a reply to this message

Tor Olav Kristensen wrote:
> If you would like to see how hyou can build functions with macros,
> then below is a macro that produces a function that is similar to
> your recursive (non-quaternion) Mandelbrot function.

Thank you very much!

All things considered, unless there is a big performance hit for some 
reason, I'd prefer to use this method to build the user-defined function 
at parse time.  The recursive formulation always struck me as an odd way 
to iterate, and povray crashing when the function is defined at the 
"wrong" time is reason enough not to stick with it.

Besides, if the quaternion road leads somewhere fruitful, trying to 
combine all the intermediate expressions by hand would be a nightmare.

With Great Appreciation,
~David Wagner

Post a reply to this message

"Tor Olav Kristensen" <tor### [at] TOBEREMOVEDgmailcom> wrote:
> [SNIP]
> // It can be used like this:
> #declare MandelBulbFn = MandelBulbFunction(8, 3, 8, pi/2, 8, 5)
>

could you please make a pov file with that for all the newbies to learn from?

Post a reply to this message

"makc" <mak### [at] gmailspamcom> wrote:
> "Tor Olav Kristensen" <tor### [at] TOBEREMOVEDgmailcom> wrote:
> > [SNIP]
> > // It can be used like this:
> > #declare MandelBulbFn = MandelBulbFunction(8, 3, 8, pi/2, 8, 5)
> >
>
> could you please make a pov file with that for all the newbies to learn from?

Ok. See the attached pov-file.
But it may not be the easiest example for new users to learn from.

Note that if one want to see more details at higher powers,
then this smoothing expression should probably be changed:

  1/(I - ln(ln(_r)/ln_R_BO)/ln_R_Pwr)

--
Tor Olav
http://subcube.com

Post a reply to this message

Attachments:
Download 'mandelbulb.pov.txt' (3 KB)

Tor Olav Kristensen wrote:
 > Ok. See the attached pov-file.

That is very nice, indeed, and far easier to read than my formulation. 
It probably renders faster than my recursive SDL hack, too.

A few quick renders seems to indicate your formula may match one of the 
bizarre "negative power" variations related to the tricorn fractal 
rather than one of the standard mandelbulb functions.  (You should still 
be able to get the "classic" mandelbulbs by messing with the phase and 
the signs of the powers.)  I don't have time right now to examine it in 
detail, but I would take a look at the use of the built-in f_ph() and 
f_th() functions; they don't correspond directly to published mandelbulb 
functions[1], nor to the conventional spherical coordinate system[2].

> But it may not be the easiest example for new users to learn from.

Plain old isosurfaces are tough for new users to tackle at all, and more 
power to those that do!

> Note that if one want to see more details at higher powers,
> then this smoothing expression should probably be changed:
> 
>   1/(I - ln(ln(_r)/ln_R_BO)/ln_R_Pwr)

Can you explain what you mean here, or post an example?

The smoothing function should not smooth the surface itself, it should 
just decrease the actual max_gradient of the function and make it easier 
for the root solver to find the surface.  To get to higher iteration 
surfaces, lower the isosurface threshold to a value closer to zero (from 
1/2 as it is now to 1/3, 1/4, and so on to get more and more fractal 
detail.)

~David

[1]
http://www.skytopia.com/project/fractal/2mandelbulb.html#formula

[2] 
http://en.wikipedia.org/wiki/Spherical_coordinate_system#Coordinate_system_conversions

Post a reply to this message

I messed with the f_mandelbulb internal function, squeezed about 10-15% better
performance out of it, and posted what I have so far below.  (Note the handy
sincos() function in the standard library saving a couple of trig calls.)

Interestingly, I also found an additional 10-15% improvement in render time by
seriously overloading the processors with one to two hundred POV-Ray worker
threads on eight cores and using smaller render block size (16 for 640*360
images).  I'm guessing the extra switching cost is more than offset by keeping
the floating-point pipelines full.

~David Wagner

DBL f_mandelbulb(FPUContext *ctx, DBL *ptr, unsigned int) // 79
{  // Coded by David Wagner, and based on
   // Daniel White's original squaring formula,
   // Paul Nylander's generalization and phi phase shift formulation,
   // the Normalized Iteration Count Algorithm,
   // and Abram Hindle's POV-Ray internal function patch.
   DBL cx=PARAM_X, cy=PARAM_Y, cz=PARAM_Z;
   DBL halfp_r   = PARAM(0);// Radial power exponent.
   DBL p_theta   = PARAM(1);// Power of elevation theta.
   DBL p_phi     = PARAM(2);// Power of azimuth phi.
   DBL phase     = PARAM(3);// Phase shift of phi.
   int i_bailout = PARAM(4);// Maximum number of iterations.
   DBL r2_bailout = pow(PARAM(5),2);// Assumed divergence radius.
   DBL x = cx, y = cy, z = cz;
   DBL r2 = x*x + y*y + z*z;
   DBL rp, sinpphi, cospphi, sinptheta, cosptheta;
   int i=1;
   halfp_r*=0.5;
   while(i < i_bailout && r2 < r2_bailout) {
       ++i;
       rp = pow(r2,halfp_r);
       sincos(p_phi * atan2(z,sqrt(x*x + y*y)) + phase, &sinpphi, &cospphi);
       sincos(p_theta * atan2(y,x), &sinptheta, &cosptheta);
       z = cz + rp * sinpphi;
       x = cx + cosptheta * rp*cospphi;
       y = cy + sinptheta * rp*cospphi;
       r2 = x*x + y*y + z*z;
   }
   if (i < i_bailout) {
       return 1/( i + log(log(r2_bailout) / log(  r2))
          / (log(halfp_r)+0.693147180559945) );
   }//          ==log(2*halfp_r)
   return 1/( i + log(log(r2_bailout) / log(1+r2))
      / (log(halfp_r)+0.693147180559945) );
}

Post a reply to this message

On 06.02.10 07:40, waggy wrote:
> I messed with the f_mandelbulb internal function, squeezed about 10-15% better
> performance out of it, and posted what I have so far below.  (Note the handy
> sincos() function in the standard library saving a couple of trig calls.)

Except that it is non-standard...

	Thorsten

Post a reply to this message

On 06.02.10 07:40, waggy wrote:
> I'm guessing the extra switching cost is more than offset by keeping
> the floating-point pipelines full.

Your guess is wrong as long as you counted "cores" as the logical number of 
cores seen by the operating system - including simultaneous multi-threading 
capabilities your processor may have.

What you are most likely seeing is that your OS gives your application 
slightly higher priority as a whole if it is running more threads, either 
because it is designed to detect compute-heavy applications or because its 
internal scheduling is on a thread rather than process basis.

	Thorsten

Post a reply to this message

Thorsten Froehlich wrote:
> On 06.02.10 07:40, waggy wrote:
 >>[...]
>> sincos() function in the standard library saving a couple of trig calls.)
> 
> Except that it is non-standard...
> 
>     Thorsten

Ah, thanks.  It's a GNU C extension.  I'll make a note to guard its use 
with a macro, and use an alternative with other compilers.

Google does not make it easy to find anything about sincos() since it 
returns results for sin_cos (where _ is any or no character).

~David

Post a reply to this message

Thorsten Froehlich wrote:
> On 06.02.10 07:40, waggy wrote:
>> I messed with the f_mandelbulb internal function, squeezed about
>> 10-15% better
>> performance out of it, and posted what I have so far below.  (Note the
>> handy
>> sincos() function in the standard library saving a couple of trig calls.)
> 
> Except that it is non-standard...

  AFAIK it's part of the C99 standard.

Post a reply to this message

Thorsten Froehlich wrote:
> On 06.02.10 07:40, waggy wrote:
>> I'm guessing the extra switching cost is more than offset by keeping
>> the floating-point pipelines full.
> 
> Your guess is wrong as long as you counted "cores" as the logical number 
> of cores seen by the operating system - including simultaneous 
> multi-threading capabilities your processor may have.
> 
I'm counting the eight physical cores in the two quad-core AMD 2344 HE 
sockets F on the motherboard.

> What you are most likely seeing is that your OS gives your application 
> slightly higher priority as a whole if it is running more threads, 
> either because it is designed to detect compute-heavy applications or 
> because its internal scheduling is on a thread rather than process basis.
> 
Interesting.  This is a on recent install of Ubuntu 9.10 Server Edition 
with their "server" kernel.  I'll do some better timing tests after I 
set up another one of these boards headless.

~David

Post a reply to this message

Warp <war### [at] tagpovrayorg> wrote:
> Thorsten Froehlich wrote:
> > On 06.02.10 07:40, waggy wrote:
> >> I messed with the f_mandelbulb internal function, squeezed about
> >> 10-15% better
> >> performance out of it, and posted what I have so far below.  (Note the
> >> handy
> >> sincos() function in the standard library saving a couple of trig calls.)
> >
> > Except that it is non-standard...
>
>   AFAIK it's part of the C99 standard.

So far as I can find, Gnu[1], Intel, HP[2], and Sun[3] include sincos() in libm
for C99.  It may take a bit more digging to sort out what compiler flags (if
any) are needed to use it with each of these.  I haven't found an answer for
Apple.

~David

[1] http://linux.die.net/man/3/sincos
[2] http://docs.hp.com/en/5992-2854/ch09s12.html
[3] http://docs.sun.com/app/docs/doc/816-5172/sincos-3m?a=view

Post a reply to this message

waggy wrote:
> So far as I can find, Gnu[1], Intel, HP[2], and Sun[3] include sincos() in libm
> for C99.  It may take a bit more digging to sort out what compiler flags (if
> any) are needed to use it with each of these.  I haven't found an answer for
> Apple.

  Well, at least with gcc if your program is C, then the compiler option
would be "-std=c99" (I really can't understand why that's not the
default). With C++ you could try "-std=c++0x".

Post a reply to this message