GPU-Accelerated Rendering of Unbounded Nonlinear Iterated Function System Fixed Points

Algorithm 3

GLSL source code generated to render Figure 1 IFS using the per-pixel inversion method. An IFS with more or different map functions will generate a different listing. Affine parameters are uniforms, to allow animation.

varying  vec2  destcoords;  //  destination  location

uniform  sampler2D  src;  //  previous  attractor

/*––  convert plane-to-texture  –––-*/

uniform  float  texscale,texscalei;

float  T(float x)  {  //  plane  to  texture

 x*=texscale;

 return  x/sqrt(1.0+x*x)*0.5+0.5;

}

vec2  T(vec2  p)  {  return  vec2(T(p·x),T(p·y));  }

float  P(float  s)  {  //  texture  to  plane

 float  u=2.0*s-1.0;

 return  texscalei*u/sqrt(1.0-u*u);

}

vec2  P(vec2 s)  { return  vec2(P(s·x),P(s·y));  }

/*–––––––– w0 ––––––––-*/

uniform  vec4  color0;

uniform  vec2  w0x,w0y,w0o;  //  w0’s  parameters

uniform  float  w0j,w0v;

float  jacobian0(vec2 t)  {return 1;}

vec4  f0(vec2 inv)  {  //  runs  at  each  inverse

 float  area=1e-2+abs(w0j*jacobian0(inv));

 vec2  t=T(w0x*inv·x+w0y*inv·y+w0o);

 return  (exp2(texture2D(src,t)*20)-1)/area;

}

vec4  nonlinear_inverse0(vec2 p)  {

 p=p*w0v;  //  Draves’  variation  coefficient

 return  f0(p);

}

/*–––––––– w1 ––––––––-*/

uniform  vec4  color1;

uniform  vec2  w1x,w1y,w1o;  //  w1’s parameters

uniform  float  w1j,w1v;

float  jacobian1(vec2 t)  {

 float  r2=dot(t,t);

 return  (1-2*t·y*t·y/r2)/r2;

}

vec4  f1(vec2 inv)  {  //  runs at each inverse

 float  area=1e-2+abs(w1j*jacobian1(inv));

 vec2  t=T(w1x*inv·x+w1y*inv.y+w1o);

 return  (exp2(texture2D(src,t)*20)-1)/area;

}

vec4 nonlinear_inverse1(vec2 p) {

 p=p*w1v; //Draves’ variation coefficient

 float ix = 0.5/p·x;

 float det=1 -  4*p·x*p·x*p·y*p·y;

 if (det>=0) {

  float sq = sqrt(det);

  return f1(vec2(ix*(1 - sq),p·y))

     +f1(vec2(ix*(1 + sq),p·y));

 } else { return vec4(0); }

}

/* Combined inverse-sampling function */

vec4 sum_inverses(vec2 p) {

 vec4 sum=vec4(0);

 sum+=color0*nonlinear_inverse0(p);

 sum+=color1*nonlinear_inverse1(p);

 return log2(sum+1)*(1.0/20);

}

void main(void) {

 gl_FragColor = sum_inverses(P(destcoords));

}