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) :
// 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
d0=ad;d1=d;
t0=t-ad;t1=t;
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:
There is practically very small overhead when using these "correction scales". Here are some evaldraw (
http://advsys.net/ken/download.htm) scripts to try:
http://implicit.googlecode.com/files/tests.zip