Dual numbers are explained here: http://en.wikipedia.org/wiki/Automatic_differentiation#Automatic_differentiation_using_dual_numbers

If your code is written in C++, you can get derivatives for "free" by just using operator overloading and a templated distance estimation formula.


I've looked at this and a few other pages - but none gave a clear explanation of how the formulas are derived. Perhaps I'm missing something. I found the same as you - 0.5 works in most cases (presumably because it is a lower bound). I was hoping that there was a more precise lower bound available that could let me take larger steps.

hi fpsunflower,


The log is the natural logarithm (base e¡2.71...), yes.

The formula comes, more or less, from this sequence of steps: (assume we have the regular quadratic case)

1. We first construct a uniformization function ph(c) (also known as Bottcher coordinates) with will map the exterior od the M set to the unit circle. It's like a deformation function, a smooth deformation in fact (with no wholes or spikes, or as mathematicians call it, a holomorphica map). Since z->z^2 + c and as z grows c has less importance, in the limit (far from the origin) z->z^2. So the uniformization map is:

phi(c) = Zn ^ (1/2^n)

where n is the number of iterations (ideally infinite). Far from the origin the map is almost identity. You can see a picture of phi(c) in the first image here: http://iquilezles.org/trastero/fieldlines/

2. The Hubbard-Douady potential is the next idea. It's defined as G(c) = log(phi(c)) with natural logarithm. This function is useful for demonstrating the connecteness of the M set (and is related to the theory of external rays).

G(c) = log |phi(c)| = log( Zn ^ (1/2^n) ).

You can think of G as a magneitc filed emanating from the M set, but increasing as you move away from it. It is zero inside M. As we said pchi(c) = c when far from M, so G = log|c| as you move away from M. You can see how it looks here: http://iquilezles.org/trastero/potential.

G = log|Zn^(1/2^n))| = log |Zn| / 2^n

As with any scalar field, you can compute the gradient of it, or in other words, the maximun slope of the field, if it was interpreted as a terrain. The gradient is related to the derivatives, so you must know how to take derivatives of logarithms. Basically,

grad(G) = |dZn| / (|Zn| * 2^n)

Now, you can approximate a smooth function with a linear version of it (a first order Taylor series) with the gradient (first ortder derivative), like you can say that G(C+eps) = G(C) + eps·|grad(G)|. If we take C inside the M set and we want to compute the distance to it in a close point (C+eps), we can isolate eps int he previous formula and get

eps = distance to M = (G(C+eps)-G(C) )/grad(G(C)). But G(C) is zero because it is inside M, so  distance = G/drad(G).  Replacing now G and grad G with what we got before, and canceling out the 2^n terms, we get:

dist = |Z|·log |Z| / |dZ|

Now, I'm missing the 2 factor, and I'm not sure why. In fact, apparently the real formula is dist = 2·|Z|·log |Z| / |dZ| and then by something called the "Koebe 1/4 theorem" we know the distance must be not more than a quarted of the previous quantity, so:

dist >= 0.5 · |Z|·log |Z| / |dZ|

Apparently the 1/4 factor holds for 3D too. The 2 that I'm missing in the numerator factor Im not sure where it comes from, probably from one of those 2^n thingies, so I would say it is related to the power of the M set indeed. In that case, the distance estimator would be

DE = (power/4) · |Z|·log |Z| / |dZ|

but Im not sure. Perhapse somebody here can confirm?

Wow, I'm away for a couple of days, and then I miss all this!

Welcome, Daniel, Iñigo, Aexion and fpsunflower. Some cracking images from you guys, I particularly enjoyed the metallic brot render by fpsunflower, and the two images by Aexion. Iñigo, 15 minutes is fast for that, gotta make an animation!

Aexion, you've got to one of my fave fractal artists ever - I love your recent "The Hexahedral Puzzle" - it has a lovely Aztec feel about it. Still think your Sunset Castle (and Return thereof) is awesome.

Paul, wow you did powers 4-8 - that must've taken a while. Look forward to doing some speed comparisons with that along with the other techniques here...


Hey, I saw that too a while back - very cool!

Bib, I'm not sure what to make of that latest blue image - maybe render a different angle and zoom in. It could be very cool...

@daniel_p i think there should be at least one scientific description of what we they are doing here, it would be a pleasure
for me to collect all the informations we have here, but it would take at least 2 or 3 monthes ( hello bachelor thesis, actually i am trying to find someone who would find that interesting )

cheers 

I try to. But navigating using the parameters variations instead of "natural" 2D zooming in UF is quite tough.

I know there is nothing revolutionnary in the images I post in this thread compared to all the ground-breaking innovations we can see here, but it's just to show that zooming into 3D fractals that look like whipped cream at first sight can reveal interesting spirals and so on...

Just a note to anyone reading the thread who may be unfamilier with |z| as used in standard maths notation and/or as used in Fractint and Ultra Fractal.

If z is complex as (x+iy) then in standard maths notation |z| = sqrt(x^2 + y^2)

BUT in Ultra Fractal and Fractint then |z| = x^2 + y^2 and cabs(z) is used as sqrt(x^2 + y^2) (cabs for complex absolute).

I was prompted to mention this after IQ's query regarding the correct scaling of the DE value, I also couldn't work out how "0.5" was arrived at. In UF and Fractint:

0.5*log(|z|) is equivalent to log(|z|) in standard math notation

And at one point I considered that maybe the 0.5 came from mistranslations between fractal program and standard math notations.

But I don't think the 0.5 comes from that, I think it really comes from the fact that the value calculated is only an *estimate* and the 0.5 is simply used to avoid cases where the estimate is larger than the true distance - but that's just my opinion
In fact the DE value I use for the Normals is assumed to be correct without the "0.5", I do however multiply the DE value by 0.5 before using the value as a step distance during ray stepping.
Also I should mention that (based on how I had to scale the DE values used to compute the normals to get correct lighting) I found (experimentally) that the visibly correct (Ultra Fractal) final DE calculation was:

            dist = 0.5*sqrt(@mpwr-1.0)*log(magn)*sqrt(magn/|dzri|)

for a Julibrot where @mpwr is the power in z^power+c, magn is |final z| (i.e. x^2+y^2) and dzri is the final derivative value.
However looking at IQ's derivation I suspect that:

            dist = 0.25*@mpwr*log(magn)*sqrt(magn/|dzri|)

May be more accurate Edit: Just tried it and it doesn't seem to be better, in fact it definitely overestimates at higher powers - I not only got gradually darker and darker shading but started getting gross overstepping as well when the method using sqrt(@mpwr-1) still worked fine.
Just to add I actually checked it on quaternions and the White/Nylander.

I'm not sure what went wrong with the z^3+c when I tried it (or the other higher powers where you use abs and sign) maybe I'm missing something obvious.
Anyway the formula I quoted here http://www.fractalforums.com/3d-fractal-generation/true-3d-mandlebrot-type-fractal/msg8475/#msg8475 can be expanded to use reals for a given power and then optimised - that formula definitely exactly matches the trig version (at least for integer powers from 2 to 60, I haven't checked fractional or negative powers).

Here is a very simple general procedure to obtain sin(p.a) and cos(p.a) if you know the value for sin(a) and cos(a) :

Let S=sin(a) and C=cos(a) and make a complex number z=S+iC, and put z1=z

Do z1=z1.z,   now real(z1)=cos(2a) and imag(z1)=sin(2a)  (z1=C^2-S^2 +i.2.S.C)

Do again z1=z1.z now real(z1)=cos(3a) and imag(z1)=sin(3a)

(because (sin(2a)+i.cos(2a))(sin(a)+i.cos(a) gives as real part sin(2a)cos(a)+cos(2a)sin(a)=sin(3a) and as imaginary part cos(2a)cos(a)-sin(2a)sin(a)=cos(3a) )

...
..(the (p-1)-th time)  real(z1)=cos(pa) and imag(z1)=sin(pa)

So no need for very long complicated formulas if you want do to a non-trig p-th power...
The drawback is that you need sin(a) and cos(a) which involves a square root.

For even powers this can be adapted so that you only need (sin(a))^2 and (cos(a))^2, which does away with the square root. However, when using a distance estimate based on the derivative, one always needs also the value for sin((p-1)a) and cos((p-1)a), so getting a value for sin(a) and cos(a) cannot be avoided then.

Also, if you want to do powers<1, it is no good of course.