// IFS System

int kIterationsPerPass = 1000;
int kMinimumPasses = 20;
int passes = 0;

int     kNbrFunctions = 3;

float fParams[][] = new float[3][6];
float weights[] = new float[5];

float   bx,by;
int     map[][];
int     warpJ = 0;
float    rad,rad2;
float   gx,gy;
boolean initedIFS=false;
float   hueMult,huePhase;
float   satMult, satPhase;

void setup()
{
  size(500,500);
  rad = width/2;
  rad2 = rad*rad;
  
  // noLoop();
  map = new int[width][height];
  colorMode(HSB, 1);
  background(0);
  initIFS();
}

void initIFS()
{
      for (int y = 0; y < height; ++y)
    {
      for (int x = 0; x < width; ++x)
      {
        map[x][y] = 0;
      }
    }
    for (int i = 0; i < kNbrFunctions; ++i)
    {
      for (int j = 0; j < 6; ++j) {
        fParams[i][j] = -1 + random(2);
      }
    }
    weights[0] = random(1);
    weights[1] = random(1-weights[0]);
    weights[2] = 1 - (weights[0]+weights[1]);
    passes = 0;
    bx = -1 + random(2);
    by = -1 + random(2);
    warpJ = int(random(8));
    background(0);
    hueMult = random(8);
    huePhase = random(2*PI);
    satMult = random(8);
    satPhase = random(2*PI);
    initedIFS = true;
}

void draw()
{
  if (mousePressed) {
    if (!initedIFS)
      initIFS();
  }
  else
    initedIFS = false;
  for (int k = 0; k < kIterationsPerPass; ++k)
  {
     passes++;
     float s = 0;

     // choose transform params r
     int r = 0;
     float rr = random(1);
     for (int i = 0; i < kNbrFunctions; ++i) 
     {
       s += weights[i];
       if (rr < s) {
          r = i;
          break;       
       }
     }
     
     float tp[];
     float tx = fParams[r][0]*bx + fParams[r][1]*by + fParams[r][2];
     float ty = fParams[r][3]*bx + fParams[r][4]*by + fParams[r][5];
     switch (warpJ) 
     {
     case 0:  warp0(tx,ty);  break;
     case 1:  warp1(tx,ty);  break;
     case 2:  warp2(tx,ty);  break;
     case 3:  warp3(tx,ty);  break;
     case 4:  warp3(tx,ty);  break;
     case 5:  warp3(tx,ty);  break;
     case 6:  warp3(tx,ty);  break;
     case 7:  warp3(tx,ty);  break;
     }
     bx = gx;
     by = gy;
     if (passes > kMinimumPasses) 
     {
       int px = (int) (rad+rad*bx);
       int py = (int)(rad+rad*by);
       if (px >= 0 && px < width && py >= 0 && py < height)
       {
         set(px,py,colorFunc(++map[px][py]));
         // symetry
         px = width-(px+1);
         set(px,py,colorFunc(++map[px][py]));
//         py = height-(py+1);
//         set(px,py,colorFunc(++map[px][py]));
//         px = width-(px+1);
//         set(px,py,colorFunc(++map[px][py]));
       }
     }
  }
}

// identity
void warp0(float x, float y)
{
  gx = x;
  gy = y;
}

// sinusoidal
void warp1(float x, float y)
{
  gx = sin(x);
  gy = sin(y);
}

// spherical
void warp2(float x, float y)
{
  float r = sqrt(x*x+y*y);
  gx = x/(r*r);
  gy = y/(r*r);
}

// swirl
void warp3(float x, float y)
{
  float r = sqrt(x*x+y*y);
  float theta = atan2(y,x);
  gx = r*cos(theta*r);
  gy = r*sin(theta*r);
}

// horsehoe
void warp4(float x, float y)
{
  float r = sqrt(x*x+y*y);
  float theta = atan2(y,x);
  gx = r*cos(theta*2);
  gy = r*sin(theta*2);
}

// polar
void warp5(float x, float y)
{
  float r = sqrt(x*x+y*y);
  float theta = atan2(y,x);
  gx = theta/PI;
  gy = r-1;
}

// hankerchief
void warp6(float x, float y)
{
  float r = sqrt(x*x+y*y);
  float theta = atan2(y,x);
  gx = r*sin(theta+r);
  gy = r*cos(theta-r);
}

// heart
void warp7(float x, float y)
{
  float r = sqrt(x*x+y*y);
  float theta = atan2(y,x);
  gx = r*sin(theta*r);
  gy = -r*cos(theta*r);
}


static float cFreq = 0.25;
color colorFunc(int x)
{
  float fx = log(x)*cFreq;
  return color(sin(hueMult*fx+huePhase)/2+.5, sin(satMult*fx+satPhase)/2+.5, fx);
}