//W:Mute 2004 - Abstract
//www.wmute.org

//This is based on Diffusion Limited Aggregation (DLA), a rough aproximation of slow additive growth processes
//like icedrawings on a window. One particle at a time walks around randomly on a lattice until it hits a growth
//core. Here DLA is extended to 3D, making it both interesting and tremendously slow. All the rest is just graphical fooling around.

int siz=80; //dimension of the cube matrix of sites
int hsiz=siz/2;
int cores=10; //number of growth cores
int index=0;
int sites[] = new int[siz*siz*siz]; //3D lattice
int coreSites[][] = new int[cores][3];//growth core positions
int curveCPoints[][]= new int[cores][3];//stores handles for core curve
boolean curveExists[][]= new boolean[cores][cores];//connection between cores exists?
boolean curveDirection[][]= new boolean[cores][cores];//direction of connection between cores
float curveOffset[][]=new float[cores][cores]; // offset of motion along curves
float curveSpeed[][]=new float[cores][cores];// speed of motion along curves
float ax[]=new float[siz*siz*siz]; //stores rotation angles of the elements
float ay[]=new float[siz*siz*siz]; //
int indices[][][]= new int[siz][siz][siz];
int trialCounter = 0; // termination condition
int trialLimit=100;
int cubeCounter = 0; // termination condition - your ninja machine is no match for D4LA :-)
int cubeLimit=10000;
float a=0.0f; //view rotation
float ac=0.0f; //position along curve
int w=400;
int hw=w/2;
float sscale = w/siz;
int output[]=new int[640*360];

void setup()
{
  size(640,360);
  colorMode(HSB);
  background(0,0,255);
  
  for (int i=0; i<siz; i++) //Empty all sites
  {
    for (int j=0;j<siz;j++)
    {
      for (int k=0;k<siz;k++)
      {
        sites[index]=0;
        ax[index]=random(TWO_PI);
        ay[index]=random(PI);
        indices[i][j][k]=index;
        index++;
      }
    }
  }
  index=0;
  for(int i=0; i<height; i++) { 
    for(int j=0; j<width; j++) { 
      output[index] = color(0,0,255);    
      index++;
    } 
  } 
  for (int i=0; i<cores; i++) //Create growth cores
  {
    int rx=int(random(siz));
    int ry=int(random(siz));
    int rz=int(random(siz));
    sites[indices[rx][ry][rz]]=int(random(21));
    coreSites[i][0]=rx;
    coreSites[i][1]=ry;
    coreSites[i][2]=rz;
    //The following stores parameters for the Bezier curves I use to spice things up
    curveCPoints[i][0]=rx+int(random(siz))-hsiz;
    curveCPoints[i][1]=ry+int(random(siz))-hsiz;
    curveCPoints[i][2]=rz+int(random(siz))-hsiz;
    curveExists[i][i]=false;
    curveDirection[i][i]=false;
    curveOffset[i][i]=0.0f;
    curveSpeed[i][i]=1.0f;
    for (int j=i; j<cores;j++)
    {
      curveExists[i][j]=curveExists[j][i]=(random(100)<80);
      curveDirection[i][j]=curveDirection[j][i]=(random(100)<50);
      curveOffset[i][j]=curveOffset[j][i]=random(1.0);
      curveSpeed[i][j]=curveSpeed[j][i]=0.6f+random(0.8);
    }
  }

}

void loop()
{
  background(0,0,255);
  //Reset on mousepressed or termination condition met
  if ((mousePressed == true)||(trialCounter>trialLimit)||(cubeCounter>cubeLimit))
  {
    trialCounter=0;
    index=0;
    for (int i=0; i<siz; i++) //Empty all
    {
      for (int j=0;j<siz;j++)
      {
        for (int k=0;k<siz;k++)
        {
          sites[index]=0;
          ax[index]=random(TWO_PI);
          ay[index]=random(PI);
          index++;
        }
      }
    }
    for (int i=0; i<cores; i++) //Growth cores
    {
      int rx=int(random(siz));
      int ry=int(random(siz));
      int rz=int(random(siz));
      sites[indices[rx][ry][rz]]=int(random(21));
      coreSites[i][0]=rx;
      coreSites[i][1]=ry;
      coreSites[i][2]=rz;
          //The following stores parameters for the Bezier curves I use to spice things up
      curveCPoints[i][0]=rx+int(random(siz))-hsiz;
      curveCPoints[i][1]=ry+int(random(siz))-hsiz;
      curveCPoints[i][2]=rz+int(random(siz))-hsiz;
      curveExists[i][i]=false;
      curveDirection[i][i]=false;
      curveOffset[i][i]=0.0f;
      curveSpeed[i][i]=1.0f;
      for (int j=i; j<cores;j++)
      {
        curveExists[i][j]=curveExists[j][i]=(random(100)<80);
        curveDirection[i][j]=curveDirection[j][i]=(random(100)<50);
        curveOffset[i][j]=curveOffset[j][i]=random(1.0);
        curveSpeed[i][j]=curveSpeed[j][i]=0.6f+random(0.8);
      }
    }
  }
  boolean mobile=true; // New mobile particle
  int cx=int(random(siz)); //Random startin position
  int cy=int(random(siz));
  int cz=int(random(siz));
  if (sites[indices[cx][cy][cz]]>0) //Starting position already occupied, try again
  {
    mobile=false;
    trialCounter++;
  }
  //Move particle randomly until collision occurs
  //Periodic boundary conditions enforced
  while (mobile) 
  {
    int movement = int(random(6)); //Pick random direction
    if (movement==0)
    {
      cx++;
      if (cx==siz) {cx=0;};
    }
    if (movement==1)
    {
      cx--;
      if (cx==-1) {cx=siz-1;};
    }
    if (movement==2)
    {
      cy++;
      if (cy==siz) {cy=0;};
    }
    if (movement==3)
    {
      cy--;
      if (cy==-1) {cy=siz-1;};
    }
    if (movement==4)
    {
      cz++;
      if (cz==siz) {cz=0;};
    }
    if (movement==5)
    {
      cz--;
      if (cz==-1) {cz=siz-1;};
    }
    //Collision check, examine neighboring sites
    //Periodic boundary conditions
    int ux = cx+1;
    if (ux==siz) {ux=0;};
    int uy = cy+1;
    if (uy==siz) {uy=0;};
    int uz = cz+1;
    if (uz==siz) {uz=0;};
    int lx = cx -1;
    if (lx==-1){lx=siz-1;};
    int ly = cy-1;
    if (ly==-1) {ly=siz-1;};
    int lz = cz-1;
    if (lz==-1) {lz=siz-1;};

    if ((sites[indices[ux][cy][cz]]>0)||(sites[indices[lx][cy][cz]]>0)||(sites[indices[cx][uy][cz]]>0)||(sites[indices[cx][ly][cz]]>0)||(sites[indices[cx][cy][lz]]>0)||(sites[indices[cx][cy][uz]]>0))
    {
      mobile=false;
      cubeCounter++;
      sites[indices[cx][cy][cz]]=sites[indices[ux][cy][cz]]+sites[indices[lx][cy][cz]]+sites[indices[cx][uy][cz]]+sites[indices[cx][ly][cz]]+sites[indices[cx][cy][lz]]+sites[indices[cx][cy][uz]];
    if (sites[indices[ux][cy][cz]]>0) {sites[indices[ux][cy][cz]]++;};
    if (sites[indices[lx][cy][cz]]>0) {sites[indices[lx][cy][cz]]++;};
    if (sites[indices[cx][uy][cz]]>0) {sites[indices[cx][uy][cz]]++;};
    if (sites[indices[cx][ly][cz]]>0) {sites[indices[cx][ly][cz]]++;};
    if (sites[indices[cx][cy][uz]]>0) {sites[indices[cx][cy][uz]]++;};
    if (sites[indices[cx][cy][lz]]>0) {sites[indices[cx][cy][lz]]++;};
    }
    index++;
  }
  //Drawing stuff, pretty straightforward
  push();
  translate(300,hw,0);
  push();
  rotateY(a+=0.002);
  rotateZ(a);
  if ((ac+=0.01)>1.0){ac-=1.0;};
  stroke(141,119,250);
 for (int i=0; i<cores;i++)
  {
    for (int j=i+1; j<cores;j++)
    {
      if (curveExists[i][j])
      {
        float lac = (ac+0.01*(curveSpeed[i][j]-1.0f))+ curveOffset[i][j];
        if (lac>1.0f){lac-=1.0f;};
        if (curveDirection[i][j]){lac=1.0f-lac;};
        float llac = lac-curveSpeed[i][j]*0.02f;
        float ulac = lac+curveSpeed[i][j]*0.02f;
        if (llac<0.0){llac=0.0f;};
        if (ulac>1.0){ulac=1.0f;};            
        //beginShape(LINE_STRIP);
        //bezierVertex(sscale*coreSites[i][0]-hw, sscale*coreSites[i][1]-hw,sscale*coreSites[i][2]-hw);  
        //bezierVertex(sscale*curveCPoints[i][0]-hw,sscale*curveCPoints[i][1]-hw,sscale*curveCPoints[i][2]-hw);
        //bezierVertex(sscale*curveCPoints[j][0]-hw,sscale*curveCPoints[j][1]-hw,sscale*curveCPoints[j][2]-hw);
        //bezierVertex(sscale*coreSites[j][0]-hw, sscale*coreSites[j][1]-hw,sscale*coreSites[j][2]-hw);    
        //endShape();              
        
          
              
        for (int k = 0; k <= 11; k++)
        {  
          float dk=k/5.0f;
          if (k>5) {dk=2.0f-dk;};
          float lt =0;
           if (curveDirection[i][j]){
          lt= llac + k * k *(ulac-llac)/ 121.0f;}
          else{
          lt= llac + k * k *(llac-ulac)/ 121.0f;}
          float lx=bezierPoint(sscale*coreSites[i][0]-hw,sscale*curveCPoints[i][0]-hw,sscale*curveCPoints[j][0]-hw,sscale*coreSites[j][0]-hw,lt);
          float ly=bezierPoint(sscale*coreSites[i][1]-hw,sscale*curveCPoints[i][1]-hw,sscale*curveCPoints[j][1]-hw,sscale*coreSites[j][1]-hw,lt);
          float lz=bezierPoint(sscale*coreSites[i][2]-hw,sscale*curveCPoints[i][2]-hw,sscale*curveCPoints[j][2]-hw,sscale*coreSites[j][2]-hw,lt);
          beginShape(LINE_STRIP);
          vertex(lx-dk, ly-dk,lz);  
          vertex(lx-dk, ly+dk,lz);
          vertex(lx+dk, ly+dk,lz);
          vertex(lx+dk, ly-dk,lz);
          vertex(lx-dk, ly-dk,lz);  
          endShape();              
          }  

       }
    }
  }
 
  index=0;
  noStroke();
  for (int i=0; i<siz; i++)
  {
    for (int j=0; j<siz; j++)
    {
      for(int k=0; k<siz;k++)
      {
        if (sites[index]>0)
        {
          fill(141,119,sites[index]);
          push();
          translate(sscale*i-hw, sscale*j-hw, sscale*k-hw);     
          push();
          rotateX(ax[index]+=0.005);
          rotateY(ay[index]+=0.005);
          box(0.7*sscale);
          pop();
          pop();
        }
        index++;
      }
    }
  }
  pop();
  pop();

  index=0;
  for(int i=0; i<height; i++) { 
    for(int j=0; j<width; j++) { 
      pixels[index] = blend(pixels[index], output[index],185) ; 
      output[index] = pixels[index];    
      index++;
    } 
  } 
}

int mix(int a,int b, int f) 
  { 
    return a+(((b-a)*f)>>8); 
  } 
int blend(int a, int b, int f) 
  { 
    return mix(a>>16,b>>16,f)<<16 | mix(a>>8&0xff,b>>8&0xff,f)<<8 | mix(a&0xff,b&0xff,f); 
  }