Logo by Kaliyuga - Contribute your own Logo!

END OF AN ERA, FRACTALFORUMS.COM IS CONTINUED ON FRACTALFORUMS.ORG

it was a great time but no longer maintainable by c.Kleinhuis contact him for any data retrieval,
thanks and see you perhaps in 10 years again

this forum will stay online for reference
News: Did you know ? you can use LaTex inside Postings on fractalforums.com!
 
*
Welcome, Guest. Please login or register. May 18, 2022, 10:15:38 PM


Login with username, password and session length


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


Pages: [1]   Go Down
  Print  
Share this topic on DiggShare this topic on FacebookShare this topic on GoogleShare this topic on RedditShare this topic on StumbleUponShare this topic on Twitter
Author Topic: Potential Optimizations for DE based Stepping  (Read 2460 times)
0 Members and 1 Guest are viewing this topic.
hobold
Fractal Bachius
*
Posts: 573


« on: February 09, 2010, 04:11:12 PM »

As I am working on my own rendering code, I keep reinventing various wheels. However, some of them don't seem to be part of the standard literature yet. So I decided to present them here, either to be used, or to be shot down, or to be recognized as a derivative of prior art. They are untested as of yet, because in the context of my own rendering paradigm, DE is just a fast lookup. Extra sophistication is likely to cost more time than it can save. But when every DE sample is being computed on demand, these things could well be worth it.


Idea number one concerns the standard practice of using distance estimates to step along a ray through a series of unbounding volumes. Typically, every next sample point is chosen to be at border of the current unbounding sphere (fudge factors notwithstanding). It seems to me that this algorithm is simple and elegant, but does not reap the largest possible benefit from each DE sample. That's because every DE sample, by its nature as an isotropic value, not only looks forward along the ray, but backward, too.

In other words, if we carefully overstep a bit, the next DE sample could be centered outside the current unbounding sphere, as long as the next unbounding sphere is large enough to touch or intersect the current one. This would not always succeed, but if the next unbounding sphere is too small, we can reliably detect that. In such a case, I would remember the two recently computed DEs for later, and compute a third one centered at the midpoint of the ray interval that was just skipped. If that third DE can close the gap: fine! If not, the third DE takes on the role of "next" DE and we recurse.

So how to choose the intentional overstepping in a meaningful way? I think a first order approximation model might do. Basically, I'd pretend that the ray is headed for an infinite plane that it eventually intersects (because we tend to aim the camera where there is something to see, and because we render smoothed surfaces, that appear flatter and flatter the closer we approach). So the DE values linearly grow or shrink as we move along the ray. The last two DEs would lead us to an "estimate of the estimate", i.e. a guess how far we can overstep for the next DE sample.

That was idea number one: a method that can potentially halve the number of steps along a ray, hence halve the number of DE computations.


Idea number two arises from the isotropic quality of distance estimates as well. Whenever we compute a distance estimate, we don't just learn something about the distance along the direction of the one ray we are computing. Instead, we learn something about a volume of space around the center point of the distance estimate.

So we could do the following: rather than compute each ray in isolation, we could regard bundles of nearby rays, centered around a "main ray" that we are currently stepping along as usual. Whenever we compute a DE for the main ray, we are really carving an unbounding sphere away from the unknown space around the fractal shape. The union of unbounding spheres along the main ray will effectively be a solid tunnel.

We could intersect the secondary rays (i.e. the immediate neighbours of the main ray) with that tunnel. Probably by intersecting them with one unbounding sphere after the other. Sooner or later the unbounding spheres become too small to intersect the secondary rays. But by then, the secondary rays will have skipped a (hopefully!) considerable portion of the empty space around the fractal object, and we saved ourselves the effort to compute quite a few DE values.

That was idea number two: a method that can advance several rays towards the surface in parallel, without computing more DE values.


I repeat: these ideas are untested. Their usefulness greatly depends on the actual cost of computing a DE value, and the overhead that arises from the extra sophistication/complexity. Is anybody aware of prior experiments with the above algorithmic tricks?
Logged
JosLeys
Strange Attractor
***
Posts: 258


WWW
« Reply #1 on: February 09, 2010, 08:30:47 PM »

I use the secondary rays to calculate the normal. What I do is to start along the secondary rays just a small distance back as measured by the main ray.
In other words, if the main ray gives a distance d from the viewpoint, I  start the process for the secondary ray at d minus something small. This avoids running the secondary rays as 'fresh' ones.
The other thing is that I just use the DE found at that  point, although this does not work for all cases. In the cases where it does, I only need one calculation per secondary ray..
Logged
knighty
Fractal Iambus
***
Posts: 819


« Reply #2 on: February 09, 2010, 09:36:20 PM »

Hi,
Cool ideas but,
Nowadays, it's hard to find new ones. wink

Something similar to your first idea is described in this article:
http://graphics.cs.uiuc.edu/~jch/papers/hyper.pdf
See title 3.4: Overshooting.

Same thing for the second:
http://graphics.cs.uiuc.edu/~jch/papers/zeno.pdf
See title 2.5.1: Image coherence. (2nd paragraph)

These articles are from the inventor of sphere tracing himself so... he devised other improvements  smiley The coolest is perhaps antialiasing with unbounding spheres!

Sometime ago I had an idea similar to your second one and have tested it using evaldraw. Thre are two stages:
1- (edit) Carve Dig (/edit) the Z-buffer using the unbounding spheres until their radius is too small.
2- improve the result using raytracing.

The script follows:
Code:
() //press F1 for help
{
   //static rparam[7][3];//for future use
   static colr,colg,colb;
   static rmin,xyinc;
   static juliadist=2.1,bradius=1.5;
   static juliac[4]={-0.25,0.9,0.0,0};
   static juliao[4]={0,0,0,0};
   static juliax[4]={1,0,0,0};
   static juliay[4]={0,1,0,0};
   static juliaz[4]={0,0,1,0};
   static p[4],tmp1[4],tmp2[4],tmp3[4];
   if (numframes==0)
   {
      cls(0);//cls(0x808080);
      clz(1e16);
      setcol(80,80,100);
      drawsph(0,0,juliadist,bradius);
      rmin=3;
      xyinc=256;
   }
   //srand(klock());
   //setcol(255,255,255);
 if (rmin>1e-3){  
  rmin*=0.5;
  xyinc=ceil(xyinc*0.5);
  if(xyinc>1){
   for (vy=0; vy<yres;vy+=xyinc)
   for (vx=0; vx<xres;vx+=xyinc)
   {
      do
      {
         z=getpix(vx,vy,colr,colg,colb);
         if (z<juliadist+bradius){
            x=(2*vx-xres)*z/xres;
            y=(2*vy-yres)*z/xres;
            jx=x;jy=y;jz=z-juliadist;
            //qrmul(jx,juliax,tmp1);qadd(juliao,tmp1,tmp2);
            //qrmul(jy,juliay,tmp1);qadd(tmp1,tmp2,tmp3);
            //qrmul(jz,juliaz,tmp1);qadd(tmp1,tmp3,p);
            qequate(juliao,p);
            p[0]+=jx;p[1]+=jy;p[2]+=jz;
            D=getJuliaDist(p,juliac,100,16);
            c=sin(60*z)*16+(z-juliadist+bradius)*128;
            setcol(c,c,c);
            drawhsph(x,y,z,D);
         } else {
            D=-100; setcol(0);
            setpix(vx,vy,1e16);
         };
      }while(D>rmin);
   }
  }else{
   for (vy=0; vy<yres;vy+=xyinc)
   for (vx=0; vx<xres;vx+=xyinc)
   {
      i=100;
      do
      {
         i--;
         z=getpix(vx,vy,colr,colg,colb);
         if (z<juliadist+bradius){
            x=(2*vx-xres)*z/xres;
            y=(2*vy-yres)*z/xres;
            jx=x;jy=y;jz=z-juliadist;
            //qrmul(jx,juliax,tmp1);qadd(juliao,tmp1,tmp2);
            //qrmul(jy,juliay,tmp1);qadd(tmp1,tmp2,tmp3);
            //qrmul(jz,juliaz,tmp1);qadd(tmp1,tmp3,p);
            qequate(juliao,p);
            p[0]+=jx;p[1]+=jy;p[2]+=jz;
            D=getJuliaDist(p,juliac,100,16);
            c=sin(60*z)*16+(z-juliadist+bradius)*128;
            setcol(c,c,c);
            setpix(vx,vy,z+D);
            //drawhsph(x,y,z,D);
         } else {
            D=-100; setcol(0);
            setpix(vx,vy,1e16);
         };
      }while((D>rmin) && i>0);
   }
  }
 }
 shadeit();  
}
shadeit()
{
   static lum[3]={0,0,-1};
   static colr,colb,colg;
   lum[0]=(mousx-0.5*xres)/xres;
   lum[1]=(mousy-0.5*yres)/yres;
   d=1/sqrt(lum[0]*lum[0]+lum[1]*lum[1]+lum[2]*lum[2]);
   lum[0]*=d;lum[1]*=d;lum[2]*=d;
   for (vy=0; vy<yres;vy++)
   for (vx=0; vx<xres;vx++)
   {
     z=getpix(vx,vy,colr,colg,colb);
     if (z<10)
     {
      x=(2*vx-xres)*z/xres;
      y=(2*vy-yres)*z/xres;
      
      zx=getpix(vx+1,vy,colr,colg,colb);
      xx=(2*(vx+1)-xres)*zx/xres;
      yx=(2*vy-yres)*zx/xres;
      
      zy=getpix(vx,vy+1,colr,colg,colb);
      xy=(2*vx-xres)*zy/xres;
      yy=(2*(vy+1)-yres)*zy/xres;
      
      zx-=z;xx-=x;yx-=y;
      zy-=z;xy-=x;yy-=y;
      
      if (0){
         tmp=zx;zx=-xx;xx=tmp;
         tmp=zy;zy=-yy;yy=tmp;
         x=xx+xy;y=yx+yy;z=zx+zy;
      } else
         {x=-yx*zy+zx*yy; y=-zx*xy+xx*zy; z=-xx*yy+yx*xy;}
      
      d=1/sqrt(x*x+y*y+z*z);
      c=(x*lum[0]+y*lum[1]+z*lum[2])*d;
      c=0.5+0.5*c+0.5*abs(c)^20;
      c*=255;
      setcol(c,c,c);
      setpix(vx,vy);
      
     }
     else {
      setcol(50,50,50);
      setpix(vx,vy);
     }
   }
}
qadd(in1[4],in2[4],out[4])
{
   out[0]=in1[0]+in2[0];out[1]=in1[1]+in2[1];
   out[2]=in1[2]+in2[2];out[3]=in1[3]+in2[3];
}
qsub(in1[4],in2[4],out[4])
{
   out[0]=in1[0]-in2[0];out[1]=in1[1]-in2[1];
   out[2]=in1[2]-in2[2];out[3]=in1[3]-in2[3];
}
qmul(in1[4],in2[4],out[4])
{
   out[0]=in1[0]*in2[0]-in1[1]*in2[1]-in1[2]*in2[2]-in1[3]*in2[3];
   out[1]=in1[0]*in2[1]+in1[1]*in2[0]+in1[2]*in2[3]-in1[3]*in2[2];
   out[2]=in1[0]*in2[2]+in1[2]*in2[0]+in1[3]*in2[1]-in1[1]*in2[3];
   out[3]=in1[0]*in2[3]+in1[3]*in2[0]+in1[1]*in2[2]-in1[2]*in2[1];
}
qrmul(r,in[4],out[4])
{
   out[0]=r*in[0];out[1]=r*in[1];out[2]=r*in[2];out[3]=r*in[3];
}
qnormsquared(in[4])
{
   return in[0]*in[0]+in[1]*in[1]+in[2]*in[2]+in[3]*in[3];
}
qnorm(in[4])
{
   return sqrt(in[0]*in[0]+in[1]*in[1]+in[2]*in[2]+in[3]*in[3]);
}
qconjugate(in[4],out[4])
{
   out[0]=in[0];out[1]=-in[1];out[2]=-in[2];out[3]=-in[3];
}
qequate(in[4],out[4])
{
   out[0]=in[0];out[1]=in[1];out[2]=in[2];out[3]=in[3];
}
getjuliadist(p[4],c[4],esc_r,max_iter)
{
   static tmp[4],g[4];
   g[0]=1; g[1]=0; g[2]=0; g[3]=0;
   i=0;
   while(i<max_iter){
      if ((qnormsquared(p)>esc_r)){
         z=qnorm(p);
         return 0.5*z*log(z)/qnorm(g);
      }    
      qmul(p,g,tmp);qrmul(2,tmp,g);
      qmul(p,p,tmp);qadd(tmp,c,p);
      i++;
   }
   return 0;
}
//slow and dirty sphere hollow
drawhsph(xc,yc,zc,r)
{
   r=abs(r);
   if (zc<r) return;//pour l'instant
   D=xres*0.5;D1=yres*0.5;
  
   //Bounding rectangle
   xmax=-1e30; ymax=-1e30;
   A=r*r-zc*zc;
   B=zc*xc;
   if (abs(A)<1e-10){
      xmin=D*(r*r-xc*xc)*0.5/B;
   } else {
      Delta=B*B-A*(r*r-xc*xc);
      A=1/A; Delta=sqrt(Delta);
      xmin=ceil(-D*(B-Delta)*A);
      xmax=ceil(-D*(B+Delta)*A);
      if (xmin<-D) xmin=-D;
      if (xmax>D) xmax=D;
   }
   A=r*r-zc*zc; B=zc*yc;
   if (abs(A)<1e-10){
      ymin=D*(r*r-yc*yc)*0.5/B;
   } else {
      Delta=B*B-A*(r*r-yc*yc);
      Delta=sqrt(Delta);
      A=1/A;
      ymin=ceil(-D*(B-Delta)*A);
      ymax=ceil(-D*(B+Delta)*A);
      if (ymin<-D1) ymin=-D1;
      if (ymax>D1) ymax=D1;
   }
  
   //A little Hack to obtain current color
   static rcol,gcol,bcol,rc,gc,bc;
   getpix(xmin+D,ymin+D1,rc,gc,bc);
   setpix(xmin+D,ymin+D1);
   getpix(xmin+D,ymin+D1,rcol,gcol,bcol);
   setcol(rc,gc,bc);
   setpix(xmin+D,ymin+D1);
   setcol(rcol,gcol,bcol);
   //let's rock (slowly!)
   C=xc*xc+yc*yc+zc*zc-r*r;r1=1/sqrt(3)/r;
   for (yp=ymin;yp<ymax;yp++){
      A1=yp*yp+D*D;B1=yp*yc+D*zc;
      for (xp=xmin;xp<xmax;xp++){
         A=A1+xp*xp; B=B1+xp*xc;
         Delta=B*B-A*C;
         if (Delta>0){
            Delta=sqrt(Delta);A=1/A;
            tmin=(B-Delta)*A;zmin=D*tmin;
            tmax=(B+Delta)*A;zmax=D*tmax;
            z=getpix(xp+D,yp+D1,rc,gc,bc);
            if (zmin<z)
               if (zmax>z){
                  //col=(r1*((xp+yp)*tmax-(xc+yc+zc)+zmax));
                  //if (col<0) col=0;
                  //col=1+col*0.5;//+col^20;
                  //setcol(col*rcol,col*gcol,col*bcol);
                  setpix(xp+D,yp+D1,zmax);
               }
         }
      }
   }
   setcol(rcol,gcol,bcol);
}
(just copy-past into evaldraw's editor then press F2. If it hangs or displays nothing -there is a bug;)- do [Esc] -options>select compiler>Ken's Eval then [Ctrl]-[Enter])

It's also possible to do coarse to fine rendering. Say... begin with 1/16 resolution with min unbounding sphere apparent radius >= 16x16 pixels then 1/4 resolution begining at previously reached  distance... etc.
« Last Edit: February 10, 2010, 09:19:17 PM by knighty » Logged
hobold
Fractal Bachius
*
Posts: 573


« Reply #3 on: February 10, 2010, 09:09:13 AM »

Nowadays, it's hard to find new ones. wink
Thank you for the links! I didn't expect my ideas to be really new, when the basic algorithm is over ten years old and many significant minds have had the opportunity to spend brain cycles on it. It's just the usual problem of not knowing the right buzzwords to google for, until you reinvented most of the details. And by then, it's really too late for googling ...
Logged
knighty
Fractal Iambus
***
Posts: 819


« Reply #4 on: February 10, 2010, 09:26:21 PM »

It's just the usual problem of not knowing the right buzzwords to google for, until you reinvented most of the details. And by then, it's really too late for googling ...

This has already happend to me too... More than once grin
Logged
Pages: [1]   Go Down
  Print  
 
Jump to:  

Related Topics
Subject Started by Replies Views Last post
Optimizations for rendering the M-set on GPU? Programming laser blaster 5 1131 Last post April 18, 2013, 10:58:07 PM
by cKleinhuis
Optimizations for bulbs other than main cardioid and period-1 bulb? Mandelbrot & Julia Set Levi 2 1783 Last post May 24, 2013, 06:37:57 PM
by Levi
Mandelbulb3DV189: Stepping the Parameters not working? Mandelbulb 3d HeartNotes3 1 1030 Last post October 19, 2013, 01:50:16 PM
by Sockratease
Fractal Flotsam - Stepping Out Mandelbulber Gallery mclarekin 0 267 Last post September 14, 2014, 05:52:50 AM
by mclarekin
Stepping Stone Mandelbulb3D Gallery 1Bryan1 0 381 Last post May 25, 2016, 08:15:22 AM
by 1Bryan1

Powered by MySQL Powered by PHP Powered by SMF 1.1.21 | SMF © 2015, Simple Machines

Valid XHTML 1.0! Valid CSS! Dilber MC Theme by HarzeM
Page created in 0.161 seconds with 26 queries. (Pretty URLs adds 0.007s, 2q)