## The All New FractalForums is now in Public Beta Testing! Visit FractalForums.org and check it out!

 Pages: [1] 2 3   Go Down
 Author Topic: Better DE estimate using orbit traps  (Read 8925 times) Description: Using orbit traps to address overstepping 0 Members and 1 Guest are viewing this topic.
knighty
Fractal Iambus

Posts: 819

 « on: January 23, 2010, 10:29:14 PM »

Hi!
First of all sorry for my terrible english
Ok! I'm not a mathematician. I found this using my very basic knowledge, trial and error and feeling .
It addresses overstepping of the distance estimator for the mandelbulb and Juliabulb.
There are two types of overstepping:
- For mandelbulb and juliabulb near the "stars"... I mean this:

The dust at the center of the above picture is due to overstepping. One can notice that these areas have the property that at some point of the Zn+1=Zn^p+C iteration, Zn comes close to the Z axis. Here is a picture that shows the areas where Zn comes close to x (red), z (green) and y (blue). note that unlike the above image, in this one "y" axis is up (and its order is 5).

The Idea is simple: just get the orbit trap of the z axis and use a function of that value to scale the estimated distance. Because we need that scaling Near the surface, It's important to scale the scale factor itself by a function of the estimated distance.

- For the juliabulb: (The overstepping is more embarrassing )

The areas where the estimated distance diverges are those where the orbit trap of (0,0,0) is very small. It happens that the radius of those regions is related to the distance of (cx,cy,cz) to the mandelbulb. ((cx,cy,cz) are the julia parameters.)

Here is how I do (hope you will understand this crap) :
Code:
// mi is the iterations number
// imax is the iterations number for the current point
// rxy is the z axis trap
// expo is the exponent (that is, "p" in Zn+1=Zn^p+C)
//First get the DE of (cx,cy,cz)
imax=0;rxy=0;//declarations
d2M=mandelbulb(cx,cy,cz,rxy,imax);
if(imax==mi) d2M=0;//(cx,cy,cz) is inside the mandelbulb
is=1-1/expo;//a magic number.
d2M=d2M^is;//a magic operation
//////////////////////////////////////////////////////////////////////////////////////////////////
mandelbulb(x,y,z,&kz,&imax){// returns: the DE, kz:correction to the "stars", imax:the number of iterations actually used
x0=x;y0=y;z0=z;
dr=1;
i=0;
rt=10000;// for (0,0,0) trap
rz=10000;// for z axis trap
while(i<mi) {
r2=x*x+y*y+z*z;
r =sqrt(r2);

rt=min(rt,r);
rz=min(rz,x*x+y*y);

if(r>2^10) break;

zo0 = asin( z/r );
zi0 = atan2( y,x );

zr = pow( r, expo-1 );
zo = zo0 * expo;
zi = zi0 * expo;

dr = zr*dr*expo + 1;
zr *= r;
x = x0 + zr*cos(zo)*cos(zi);
y = y0 + zr*cos(zo)*sin(zi);
z = z0 + zr*sin(zo);

i++;
}
imax=i;
kz=1-(1/(16*rz+1))*(1-1/((2*rt)^expo+1));//this will be used in the ray marcher
return 0.5*log(r)*r/dr;
}
//////////////////////////////////////////////////////////////////////////////////////////////////
juliabulb(x,y,z,&kz,&imax){
x0=cx;y0=cy;z0=cz;
dr=1;
i=0;
rt=10000;
rz=10000;
while(i<mi) {
r2=x*x+y*y+z*z;
r =sqrt(r2);
rt=min(rt,r);
rz=min(rz,x*x+y*y);

if(r>2^10) break;

zo0 = asin( z/r );
zi0 = atan2( y,x );

zr = pow( r, expo-1 );
zo = zo0 * expo;
zi = zi0 * expo;

dr = zr*dr*expo ;
zr *= r;
x = x0 + zr*cos(zo)*cos(zi);
y = y0 + zr*cos(zo)*sin(zi);
z = z0 + zr*sin(zo);

i++;
}
imax=i;
dist=0.5*log(r)*r/(dr+2^-53);//Estimated Distance
k=1-d2M/(rt^(expo-1)+d2M);//Correction to the estimated distance due to (0,0,0) trap
kz=1-(1/(16*rz+1))*(1-1/((2*rt)^expo+1));//Correction to the estimated distance due to the "z" axis trap. used in the ray marcher.
return dist*k;//This scaling can be done here
}
/////////////////////////////////////////////////////////////////////////////////////
// The ray marcher: returns the point on the ray ((ox,oy,oz),(vx,vy,vz)) where the intersection occures
nextd0(ox,oy,oz,vx,vy,vz,tmin,tmax)
{
ldist=min(0.05,max(errmax,param3*2^-8));//Don't care. ldist is "targeted" distance to the fractal.
rxy=0;imax=0;
t=tmin;
for(i=0;i<maxiter;i++){
d=juliabulb(ox+t*vx,oy+t*vy,oz+t*vz,rxy,imax);//or mandelbulb(x,y,z,rxy,imax)
if (imax==mi) d=0;
rxy-=1;rxy*=exp(-20*d^2);rxy+=1;// scale the "star" correction term by a function of d.
dscale=(0.95*rxy+0.05);//the final scale. here: 0.05 is the minimal step distance. 0.95+0.05=1.
//One can use ((1-min_step_distance)*rxy+min_step_distance).
if(d-t*ldist<errmax*t) {// Reached an intersection: Refine
for(j=0;j<5;j++){
t=0.5*(t0+t1);
d=juliabulb(ox+t*vx,oy+t*vy,oz+t*vz,rxy,imax);//or mandelbulb(x,y,z,rxy,imax)
if((d0-t0*ldist)*(d-t*ldist)<0){t1=t;d1=d;}
else{t0=t;d0=d;}
}
t=0.5*(t0+t1);//Why not? but t may be inside!
break;
}
if(t>tmax) {break;}
t+=d*dscale;//We step forward by the scaled distance
}
return t;
}

Unlike the "julia" correction scale, the "star" correction scale have a formula that seems to be ad-hoc. Well... it is . It's "very" empirical and there is room for improvement. The formula for the "julia" correction scale is cleaner. That doesn't mean that I really know the mathematical reasons behind it

The results:

 « Last Edit: January 23, 2010, 10:31:14 PM by knighty » Logged
David Makin
Global Moderator
Fractal Senior

Posts: 2286

 « Reply #1 on: January 23, 2010, 11:44:49 PM »

Awesome !

This (and something I mentioned to Jos Leys privately) gives me an idea.....will post it in this thread shortly
 Logged

The meaning and purpose of life is to give life purpose and meaning.

http://www.fractalgallery.co.uk/
"Makin' Magic Music" on Jango
David Makin
Global Moderator
Fractal Senior

Posts: 2286

 « Reply #2 on: January 24, 2010, 12:05:40 AM »

OK, here are two images, both solid based on DE alone.

Scaled smooth iteration colouring:

Orbit trap (to origin) colouring:

Spot the difference !!

So the suggestion of using orbit trap values to fix the DE overstepping along with the above amazing similarity between the two colouring methods prompted me to think that maybe we can get a distance estimate more simply by using the smooth iteration value and the orbit trap (to origin) value using something like:

Distance = a*f(smooth iteration) - g(orbit trap)

Obviously this requires more work but if suitable f() and g() can be found then it should be even faster than the normal analytical method since the only code in the iteration loop beyond the main iteration would be the orbit trap test.
 « Last Edit: January 24, 2010, 05:52:26 PM by David Makin » Logged

The meaning and purpose of life is to give life purpose and meaning.

http://www.fractalgallery.co.uk/
"Makin' Magic Music" on Jango
JosLeys
Strange Attractor

Posts: 258

 « Reply #3 on: January 24, 2010, 01:12:25 AM »

Code:
rxy-=1;rxy*=exp(-20*d^2);rxy+=1;// scale the "star" correction term by a function of d.

Could you please give the explicit expression for rxy?
I'm confused about what you do here.

I'm trying to implement this, and it looks like it works very well!
 Logged
knighty
Fractal Iambus

Posts: 819

 « Reply #4 on: January 24, 2010, 01:21:26 PM »

I forgot to say that even if I'm happy with the results, I'm not satisfyed by the formulas. I hope that someone with better math background will find the good formulas with their justification.

Code:
rxy-=1;rxy*=exp(-20*d^2);rxy+=1;// scale the "star" correction term by a function of d.

Could you please give the explicit expression for rxy?
I'm confused about what you do here.

I'm trying to implement this, and it looks like it works very well!

"rxy" is named "kz" in juliabulb(). its value is: 1-(1/(16*rz+1))*(1-1/((2*rt)^expo+1));
I could have done it more simple:
in juliabulb(): kz(aka rxy)=1/(16*rz+1)*(1-1/((2*rt)^expo+1));
in nextd0():rxy=1-rxy*exp(-20*d^2);

OK, here are two images, both solid based on DE alone.

Scaled smooth iteration:

<Quoted Image Removed>

Orbit trap (to origin):

<Quoted Image Removed>

Spot the difference !!

So the suggestion of using orbit trap values to fix the DE overstepping along with the above amazing similarity between the two colouring methods prompted me to think that maybe we can get a distance estimate more simply by using the smooth iteration value and the orbit trap (to origin) value using something like:

Distance = a*f(smooth iteration) - g(orbit trap)

Obviously this requires more work but if suitable f() and g() can be found then it should be even faster than the normal analytical method since the only code in the iteration loop beyond the main iteration would be the orbit trap test.

Interresting! What happens around the surface? maybe a 2D slice will show the relation between the two values.
 Logged
David Makin
Global Moderator
Fractal Senior

Posts: 2286

 « Reply #5 on: January 24, 2010, 02:49:57 PM »

One thing I discovered that I should also mention.

When changing the DE "solid threshold" the smooth iteration colouring normally "drifts" because obviously at smaller distances we get higher iteration counts - this is particularly a problem when using a dynamic distance threshold based on distance from viewpoint and image plane distance.
When investigating this I discovered that merely offsetting UF's palette as the solid threshold changed was fixing the problem and surmised that maybe simply adding a multiple of the distance threshold to the smooth iteration value would fix the problem.
It turns out I was correct and not only that but the scale required appears to be a universal constant of approx. 0.521 i.e. it works for the various different Mandelbulb powers and for all the different variations of the Mandelbulb formula that I have implimented (about 15).
Basically this means that you calculate your smooth iteration value for the solid point and then adjust as follows:

smiter = smiter + 0.521*threshold

To see this in action I just uploaded an animation to YouTube:

Edit: I should add that (surprisingly) the apparently universal nature of the constant 0.521 making the results consistent under changes to the solid threshold seems to hold for different bailout values (e.g. from testing the square of the magnitude against values ranging from 4 to 1e20) though changes in bailout do require a further offset to the smooth iteration value to make the results consistent with different bailouts.
 « Last Edit: January 24, 2010, 03:12:42 PM by David Makin » Logged

The meaning and purpose of life is to give life purpose and meaning.

http://www.fractalgallery.co.uk/
"Makin' Magic Music" on Jango
kram1032
Fractal Senior

Posts: 1863

 « Reply #6 on: January 24, 2010, 04:28:49 PM »

Interesting
Maybe, the exact value is represented by one of these?
http://www.wolframalpha.com/input/?i=0.521

$\pi\over6$ is just above that...
 Logged
David Makin
Global Moderator
Fractal Senior

Posts: 2286

 « Reply #7 on: January 24, 2010, 05:51:34 PM »

Interesting
Maybe, the exact value is represented by one of these?
http://www.wolframalpha.com/input/?i=0.521

<Quoted Image Removed> is just above that...

sqr(tan(5/8)) looks favourite
 Logged

The meaning and purpose of life is to give life purpose and meaning.

http://www.fractalgallery.co.uk/
"Makin' Magic Music" on Jango
kram1032
Fractal Senior

Posts: 1863

 « Reply #8 on: January 24, 2010, 09:37:28 PM »

if you could find the analythical value, it will quite likely be related to tangens, I guess

Note, that if you click "more", you'll get more possible values, one being even more exact, as in, accurate to also the 1, however, they also start to get more complex and I guess, it's rather unlikely that they are correct...
 Logged
knighty
Fractal Iambus

Posts: 819

 « Reply #9 on: January 25, 2010, 10:04:58 PM »

I've made some small improvements to the correction terms formulas.

1- For the juliabulb:
power: the fractal's exponent.
c: the Juliabulb seed.
d2M: the distance from c to the Mandelbulb.
rt: the origin trap. rt=min|zi|
ED: the estimated distance. ED = 0.5*r*log(r)/dr.
cd: the correction term. cd = 1-d2M^(1-1/power)/(rt^(power-1)+d2M^(1-1/power)) = 1-1/(rt^(power-1)/d2M^(1-1/power)+1).
The corrected distance estimation will be: ED*cd.
It's better (in order to avoid some floating point issues) to pack ED*cd.
After some basic algebra manipulation: ED*cd = 0.5*log(r)*r*rt/(dr*(rt^(expo-1)+d2M^(1-1/power))).

• Star overshoot:
rxy: the z axis trap. rxy = min(zi.x^2+zi.y^2).
cs: the correction term. cs =1- (1/(16*rz+1))*(1-d2M^(1-1/power)/(rt^(power-1)+d2M^(1-1/power)))*exp(-20*(ED-SOLID_THRESHOLD)^2).
The second term [(1-d2M^(1-1/power)/(rt^(power-1)] prevents using the correction "inside" the juliabulb (that is, the hollow spheres).
The third term [exp(-20*(ED-SOLID_THRESHOLD)^2)] prevents using the correction far from the rendered surface.

2- For the mandelbulb:
We need only the star-overshoot correction term, because there is no "inside". Thus, we also don't need the second term discussed above. Thus:
cs =1- (1/(16*rz+1))*exp(-20*(ED-SOLID_THRESHOLD)^2).

Note that the coefficients that are underlined may be changed.

I hope you'll find it useful.
 Logged
David Makin
Global Moderator
Fractal Senior

Posts: 2286

 « Reply #10 on: January 26, 2010, 02:04:30 AM »

2- For the mandelbulb:
We need only the star-overshoot correction term, because there is no "inside". Thus, we also don't need the second term discussed above.

Did you check for the Julia-style problem on perturbated Mandelbrots ?
 Logged

The meaning and purpose of life is to give life purpose and meaning.

http://www.fractalgallery.co.uk/
"Makin' Magic Music" on Jango
knighty
Fractal Iambus

Posts: 819

 « Reply #11 on: January 26, 2010, 05:37:28 PM »

Oh! I've forgotten to say that in the case of the mandelbulb It's better to skip the z axis trap in the first iteration.
Anyway the "star" correction term is not "perfect". It doesn't work well with power 2. I'll try to find better one.

2- For the mandelbulb:
We need only the star-overshoot correction term, because there is no "inside". Thus, we also don't need the second term discussed above.

Did you check for the Julia-style problem on perturbated Mandelbrots ?

By "perturbated Mandelbrots", do you mean the morphing between julia and mandelbrot bulbs described here?
http://www.fractalforums.com/mandelbulb-implementation/mandelbulb-ray-tracing-plugin-for-photoshop-aftereffects-and-quartzcomposer/msg10121/#msg10121

Well! in principle, one have to recompute d2M (distance to mandelbrot) for each point, but this is not sufficient.

It works well in general, when calculating d2M only for the julia seed, but there are some areas where the estimated distance is too large. I think this problem occures even if one doesn't use the "correction terms" because of the difference between the ways "dr" is calculated in mandelbrot and in julia (the +1).

The result seem to be perfect when recomputing d2M for each point AND adding this term:

cd*MJI

to dri. where:
• cd=1-d2M^(1-1/power)/(rt^(power-1)+d2M^(1-1/power)) . //the same as defined in my previous post.
• MJI is the morphing parameter. MJI=0 for pure julia and 1 for pure mandel.

Well! this is quite slow.

In the case where one compute d2M once for the julia seed and add cd*MJI to dri, the results seem to be good. For now, I only have done the tests on 2D slices. I'll see what they give in 3D and post the results.
 Logged
David Makin
Global Moderator
Fractal Senior

Posts: 2286

 « Reply #12 on: January 26, 2010, 06:21:53 PM »

By "perturbated Mandelbrots", do you mean the morphing between julia and mandelbrot bulbs described here?

I just meant Mandelbrots where the start value is not zero, which can have the same stepping problem that you get in Julia Sets.

Edit: More generally I should say I mean Mandelbrots where the starting value is not a critical value
 « Last Edit: January 26, 2010, 09:22:14 PM by David Makin » Logged

The meaning and purpose of life is to give life purpose and meaning.

http://www.fractalgallery.co.uk/
"Makin' Magic Music" on Jango
knighty
Fractal Iambus

Posts: 819

 « Reply #13 on: January 26, 2010, 10:11:05 PM »

sorry! I'm not familiar with fractals' terminology.
It seems that the origin trap is useless in this case. There is a a kind of shift (perhaps because the 0 period attractor is not 0?). Moreover I have no idea about the size of the regions where overshooting occures. I'll try to investigate further but it seem more difficult  . Any idea?  :

About my previous post. It's quite strange that "cd*MJI" works better than just adding "MJI" to "dri" because the former looks more logical.
Mistery!

What about the 0.521? I'm curious about its origin and closed form expression (if it exists).
 Logged
David Makin
Global Moderator
Fractal Senior

Posts: 2286

 « Reply #14 on: January 27, 2010, 12:08:11 AM »

What about the 0.521? I'm curious about its origin and closed form expression (if it exists).

The 0.521 is purely experimental/observational in origin - if you render a 3D+ fractal using solid based solely on distance estimate threshold and then colour it using (accurate) smooth iteration colouring for the surface points (at threshold distance) then it appears that adjusting the smooth iteration value by (approx) 0.521*threshold keeps the colouring consistent as threshold changes.
More accurately it keeps the "iteration densest" regions consistent, the less dense iteration regions (which appear as threshold is reduced) have higher iteration counts and hence higher colouring values - see the "Snow" animation - the dark areas only appear as solid threshold is reduced.
I find it astonishing that the 0.521 value appears consistent not only for the z^p+c Mandelbulbs/Juliabulbs (all "bulb-type" maths and seemingly for all values of p>=2) but also works on hypercomplex fractals (bicomplex, quaternion etc).
Actually that prompts a question - is it also consistent on 2D fractals ? Hmmmmmm......

Anyone fancy finding the analytical route to the 0.521 ?
 Logged

The meaning and purpose of life is to give life purpose and meaning.

http://www.fractalgallery.co.uk/
"Makin' Magic Music" on Jango
 Pages: [1] 2 3   Go Down