Welcome to Fractal Forums

//IFS distance estimation proof of concept.
//It's derived from Hart's well known algorithm.
//Instead of intersecting a ray with bounding volumes we compute distance to them
//The algorithm needs a priority queue.
//The one implemented here is cartainly the most stupid of all: just arrays. Insertions are O(n) (very bad) and deletions O(1) (optimal)
//A better priority queue would be heap based: O(log(n)) insertion and deletion.
//It's better also because we do much more insertions than deletions.
//An even priority queue would be the one where insertions are even less costly.
//This script shows that for our purpose a (very short) stack is sufficient.
//There should be a relation between number of transformations and 'optimal' queue size and scales of the transformations.
(x,y){
   static binit=1;
   static frm=-1,mx,my,mr2;
   if(frm<numframes){
      frm=numframes;
      mx=(mousx-0.5*xres)*0.01;
      my=(0.5*yres-mousy)*0.01;
      mr2=deifs(mx,my);mr2*=mr2;
   }
   if((x-mx)^2+(y-my)^2<mr2) return 0;
   //x=abs(x);
   if(binit){//initialization of the transformations. Called once.
      binit=0;
      settrns(0,-1+0.25,-0.04,-0.25*PI,sqrt(2));
      settrns(1, 0.25,-0.04,1.25*PI,sqrt(2));
      //settrns(2, 0+0.44,   1+0.05,0.5,1.7);
   }
   deifs(x,y)^0.25
}
enum {Ntrns=2};//number of transformations.
enum {tcx=0,tcy=1,trr=2,tir=3,tis=4};//names for transformations entries.
static trns[Ntrns][8];//transformations array
settrns(i,scx,scy,sa,ss){//transformation set up function
   if(i>=Ntrns) return;//should be an error
   trns[i][tcx]=scx;//x center
   trns[i][tcy]=scy;//y center
   trns[i][trr]=ss*cos(sa);//x similarity
   trns[i][tir]=ss*sin(sa);//y similarity
   trns[i][tis]=abs(1/ss);//inverse of scaling factor
}
enum {Npq=16};//Size of the priority queue

deifs(x,y){
   //priority queue
      auto pqx[Npq];//x
      auto pqy[Npq];//y
      auto pqs[Npq];//1/accumulated_scale
      auto pqd[Npq];//distance
   
   stptr=0;//pointer in the priority queue
   bvr=1.27;//bounding volume (a disc centered at 0,0) radius
   bvrm=(4*bvr)^2;//bailout: controls the smoothness of the DE field (lower values faster, higher values smoother)
   mm=0;//to compute mean movement number in the priority queu
   {//push first instance
      pqx[stptr]=x; pqy[stptr]=y;
      pqs[stptr]=1; pqd[stptr]=sqrt(x*x+y*y)-bvr;
      stptr++;
   }
   bNoExit=(x*x+y*y<bvrm);i=0;ii=0;
   while(bNoExit && i<30){
      //extract nearest element
      {
         stptr--;
         x=pqx[stptr]; y=pqy[stptr];s=pqs[stptr];ad=pqd[stptr];
      }
      //generate and insert its children
      for(j=0;j<Ntrns;j++){
         //generate child j
         pqx[stptr]=(x-trns[j][tcx])*trns[j][trr]-(y-trns[j][tcy])*trns[j][tir]+trns[j][tcx];
         pqy[stptr]=(x-trns[j][tcx])*trns[j][tir]+(y-trns[j][tcy])*trns[j][trr]+trns[j][tcy];
         pqs[stptr]=s*trns[j][tis];
         pqd[stptr]=(sqrt((pqx[stptr])^2+(pqy[stptr])^2)-bvr)*pqs[stptr];
         pqd[stptr]=max(ad,pqd[stptr]);//this removes discontinuities when the bounding volume doesn't cover it's images

         //insert it at its place
         ptr=stptr;
         while(ptr>0){
            if(pqd[ptr]>pqd[ptr-1]){//swap
               k=pqx[ptr]; pqx[ptr]=pqx[ptr-1]; pqx[ptr-1]=k;
               k=pqy[ptr]; pqy[ptr]=pqy[ptr-1]; pqy[ptr-1]=k;
               k=pqs[ptr]; pqs[ptr]=pqs[ptr-1]; pqs[ptr-1]=k;
               k=pqd[ptr]; pqd[ptr]=pqd[ptr-1]; pqd[ptr-1]=k;
               //k=pqh[ptr]; pqh[ptr]=pqh[ptr-1]; pqh[ptr-1]=k;
               ptr--;
            } else break;
         }
         stptr++;
         if(stptr>=Npq){//remove last elements... The dumb way ;)
            n2r=1;//number of elements to remove, usually 1 but here it's equal to number of transformations
            for(ii=0;ii<Npq-n2r;ii++){
               pqx[ii]=pqx[ii+n2r];
               pqy[ii]=pqy[ii+n2r];
               pqs[ii]=pqs[ii+n2r];
               pqd[ii]=pqd[ii+n2r];
            }
            stptr-=n2r;
         }
      }
      bNoExit=(x*x+y*y<bvrm);
      i++;
   };
   pqd[stptr-1]
}

(x,y,t){
   static binit=1;
   if(binit){
      //binit=0;
      settrns(0,-0.4,0.25,0.2*t*0.5*pi,2);
      settrns(1, 0.4,0.25,0.2*t*0.5*pi,2);
      settrns(2, 0+0.1,   1+0.0,0.5,1.7);
   }
   //rect(x,y,1,1)
   deifs(x,y)^0.25
}
static recta=0.62,rectb=0.15;
rect(x,y,a,b){
   x=abs(x)-a;y=abs(y)-b;
   if(min(x,y)>0 && max(x,y)>0) return sqrt(x*x+y*y);
   max(x,y)
}
enum {Ntrns=2};
enum {tcx=0,tcy=1,trr=2,tir=3,tis=4};
static trns[Ntrns][8];
settrns(i,scx,scy,sa,ss){
   if(i>=Ntrns) return;
   trns[i][tcx]=scx;
   trns[i][tcy]=scy;
   trns[i][trr]=ss*cos(sa);
   trns[i][tir]=ss*sin(sa);
   trns[i][tis]=abs(1/ss);
}
enum {Npq=256};
//enum {cx=0,cy=1,cs=2,cd=3};
deifs(x,y){
   maxdepth=5;
   //priority queu
   auto pqx[Npq];//x
   auto pqy[Npq];//y
   auto pqs[Npq];//1/accumulated_scale
   auto pqd[Npq];//distance
   auto pqh[Npq];//depth
   stptr=0;//pointer in the priority queu
   bvr=1.2;//bounding volume (a disc centered at 0,0 for now) radius
   bvrm=100;//bailout. it controls the smoothness of the DE field (lower values faster, higher values smoother)
   mm=0;//to compute mean movement number in the priority queu
   {//first instance
      pqx[stptr]=x; pqy[stptr]=y;
      pqs[stptr]=1; pqd[stptr]=sqrt(x*x+y*y)-bvr;
      pqh[stptr]=0;
      stptr++;
   }
   bNoExit=(x*x+y*y<bvrm*bvr*bvr);i=0;ii=0;
   while(bNoExit && i<100){
      //extract nearest element
      {
         depth=pqh[stptr-1];
         if(depth>maxdepth+1) break;
         stptr--;
         if(depth>maxdepth) if(stptr==0) break; else continue;
         x=pqx[stptr]; y=pqy[stptr];s=pqs[stptr];ad=pqd[stptr];
      }
      //generate and insert its children
      for(j=0;j<Ntrns;j++){
         //generate child j
         pqx[stptr]=(x-trns[j][tcx])*trns[j][trr]-(y-trns[j][tcy])*trns[j][tir]+trns[j][tcx];
         pqy[stptr]=(x-trns[j][tcx])*trns[j][tir]+(y-trns[j][tcy])*trns[j][trr]+trns[j][tcy];
         pqs[stptr]=s*trns[j][tis];
         pqd[stptr]=(sqrt((pqx[stptr])^2+(pqy[stptr])^2)-bvr)*pqs[stptr];
         pqd[stptr]=max(ad,pqd[stptr]);//this removes discontinuities when the bounding volume doesn't cover it's images
         pqh[stptr]=depth+1;
         
         //insert it at its place in the priority queue
         ptr=stptr;
         imm=0;
         while(ptr>0){
            if(pqd[ptr]>pqd[ptr-1]){//swap
               k=pqx[ptr]; pqx[ptr]=pqx[ptr-1]; pqx[ptr-1]=k;
               k=pqy[ptr]; pqy[ptr]=pqy[ptr-1]; pqy[ptr-1]=k;
               k=pqs[ptr]; pqs[ptr]=pqs[ptr-1]; pqs[ptr-1]=k;
               k=pqd[ptr]; pqd[ptr]=pqd[ptr-1]; pqd[ptr-1]=k;
               k=pqh[ptr]; pqh[ptr]=pqh[ptr-1]; pqh[ptr-1]=k;
               ptr--;imm++;
            } else break;
         }
         stptr++;
      }
         //generate internal chape
         pqd[stptr]=max(ad,s*rect(x,y,recta,rectb));
         pqh[stptr]=maxdepth+2;
         //insert it at its place in the priority queue
         ptr=stptr;
         imm=0;
         while(ptr>0){
            if(pqd[ptr]>pqd[ptr-1]){//swap
               k=pqx[ptr]; pqx[ptr]=pqx[ptr-1]; pqx[ptr-1]=k;
               k=pqy[ptr]; pqy[ptr]=pqy[ptr-1]; pqy[ptr-1]=k;
               k=pqs[ptr]; pqs[ptr]=pqs[ptr-1]; pqs[ptr-1]=k;
               k=pqd[ptr]; pqd[ptr]=pqd[ptr-1]; pqd[ptr-1]=k;
               k=pqh[ptr]; pqh[ptr]=pqh[ptr-1]; pqh[ptr-1]=k;
               ptr--;imm++;
            } else break;
         }
         mm=mm+imm;ii++;
         stptr++;

      //d=pqd[stptr-1];x=pqx[stptr-1];y=pqy[stptr-1];s=pqs[stptr-1];
      bNoExit=(x*x+y*y<bvrm*bvr*bvr);//(d<0.01*cs*bvr);//||((d>0)&&(d-ad>0.01*ad));//
      i++;
   };
   pqd[stptr-1]//minde//
}