int w = 300, h = 300;
float[] phase, newPhase;
float speed = 0.3;
float[] sinLUT;
float LUTSize;
Thing[] things;

void setup() {
  size(w,h,P3D);
  phase = new float[w*h];
  newPhase = new float[w*h];
  
  int lutSize = 4000;
  LUTSize = (float) lutSize;
  sinLUT = new float[lutSize];
  for (int i=0; i<lutSize; i++) {
    sinLUT[i] = sin(TWO_PI / LUTSize * (float) i);
  }
  
  things = new Thing[35];
  for (int i=0; i<things.length; i++) {
    things[i] = new Thing(things);
  }
}


void draw() {
  if (mousePressed) {
    if (mouseX>=0 && mouseX<width && mouseY>=0 && mouseY<height) {
    final int pos = mouseY*w+mouseX;
    phase[pos] += 1000;
    }
  }
  
  for (int y=1; y<h-1; y++) {
    for (int x=1; x<w-1; x++) {
      final int pos = y*w+x;
      newPhase[pos] = ( phase[pos-1] + phase[pos+1] + phase[pos-w] + phase[pos+w] ) / (4.0 / 0.8) + ( phase[pos-1-w] + phase[pos+1-w] + phase[pos-1+w] + phase[pos+1+w] ) / (4.0 / 0.2);
      newPhase[pos] = ( newPhase[pos]+speed );
    }
  }

  // special cases
  for (int x=1; x<w-1; x++) {
      int pos = x;
      newPhase[pos] = ( phase[pos-1] + phase[pos+1] + phase[pos+w] ) / (3.0 / 0.8) + ( phase[pos-1+w] + phase[pos+1+w] ) / (2.0 / 0.2);
      newPhase[pos] = ( newPhase[pos]+speed );
      pos = (h-1)*w+x;
      newPhase[pos] = ( phase[pos-1] + phase[pos+1] + phase[pos-w] ) / (3.0 / 0.8) + ( phase[pos-1-w] + phase[pos+1-w] ) / (2.0 / 0.2);
      newPhase[pos] = ( newPhase[pos]+speed );
  }
  for (int y=1; y<w-1; y++) {
      int pos = y*w;
      newPhase[pos] = ( phase[pos-w] + phase[pos+1] + phase[pos+w] ) / (3.0 / 0.8) + ( phase[pos-w+1] + phase[pos+w+1] ) / (2.0 / 0.2);
      newPhase[pos] = ( newPhase[pos]+speed ) ;
      pos = y*w+w-1;
      newPhase[pos] = ( phase[pos-w] + phase[pos-1] + phase[pos+w] ) / (3.0 / 0.8) + ( phase[pos-w-1] + phase[pos+w-1] ) / (2.0 / 0.2);
      newPhase[pos] = ( newPhase[pos]+speed ) ;
  }
  int pos = 0;
  newPhase[pos] = ( phase[pos+1] + phase[pos+w] ) / (2.0 / 0.8) + phase[pos+w+1] * 0.2;
  newPhase[pos] = ( newPhase[pos]+speed ) ;
  pos = w-1;
  newPhase[pos] = ( phase[pos-1] + phase[pos+w] ) / (2.0 / 0.8) + phase[pos+w-1] * 0.2;
  newPhase[pos] = ( newPhase[pos]+speed ) ;
  pos = w*(h-1);
  newPhase[pos] = ( phase[pos+1] + phase[pos-w] ) / (2.0 / 0.8) + phase[pos-w+1] * 0.2;
  newPhase[pos] = ( newPhase[pos]+speed ) ;
  pos = w*h-1;
  newPhase[pos] = ( phase[pos-1] + phase[pos-w] ) / (2.0 / 0.8) + phase[pos-w-1] * 0.2;
  newPhase[pos] = ( newPhase[pos]+speed ) ;
  
    
  loadPixels();
  for (int i=0; i<pixels.length; i++) {
    int c = (int) (sinLUT[(int) (( ( ( (newPhase[i]%TWO_PI)+TWO_PI)%TWO_PI)/TWO_PI)*LUTSize)]*127.0) + 128;
    pixels[i] = 0xFF000000 | c<<16 | c<<8 | c;
  }
  updatePixels();
  for (int i=0; i<things.length; i++) {
    things[i].process();
  }
  
  // switch phase and newPhase
  float[] swtch = phase;
  phase = newPhase;
  newPhase = swtch;
}



class Thing {
  int posX, posY;
  Thing[] allThings;
  
  public Thing(Thing[] arr) {
    allThings = arr;
    randomPos();
    for (int i=0; i<allThings.length; i++) {
      if (allThings[i] == null) {
        allThings[i] = this;
        break;
      } else {
        if (allThings[i].posX == posX && allThings[i].posY == posY) {
          randomPos();
          i=-1;
        }
      }
    }
  }
  
  public void randomPos() {
    posX = (int) random(w);
    posY = (int) random(h);
  }
  
  public void process() {
    newPhase[posY*w+posX] = 0;
    float maxPhase = newPhase[posY*w+posX];
    int mPosX = posX;
    int mPosY = posY;
    float c;
    int minSqrDist = 400;
    if (posX>0 && posY>0) {
      c = newPhase[(posY-1)*w+posX-1];
      int nn = nearestNeighbor(posX-1,posY-1);
      if (c>maxPhase && (nn>minSqrDist ^ nn>nearestNeighbor(posX,posY)) ) {
        maxPhase = c;
        mPosX = posX-1;
        mPosY = posY-1;
      }
    }
    if (posY>0) {
      c = newPhase[(posY-1)*w+posX];
      int nn = nearestNeighbor(posX,posY-1);
      if (c>maxPhase && (nn>minSqrDist ^ nn>nearestNeighbor(posX,posY)) ) {
        maxPhase = c;
        mPosX = posX;
        mPosY = posY-1;
      }
    }
    if (posX<w-1 && posY>0) {
      c = newPhase[(posY-1)*w+posX+1];
      int nn = nearestNeighbor(posX+1,posY-1);
      if (c>maxPhase && (nn>minSqrDist ^ nn>nearestNeighbor(posX,posY)) ) {
        maxPhase = c;
        mPosX = posX+1;
        mPosY = posY-1;
      }
    }
    if (posX>0) {
      c = newPhase[(posY)*w+posX-1];
      int nn = nearestNeighbor(posX-1,posY);
      if (c>maxPhase && (nn>minSqrDist ^ nn>nearestNeighbor(posX,posY)) ) {
        maxPhase = c;
        mPosX = posX-1;
        mPosY = posY;
      }
    }
    if (posX<w-1) {
      c = newPhase[(posY)*w+posX+1];
      int nn = nearestNeighbor(posX+1,posY);
      if (c>maxPhase && (nn>minSqrDist ^ nn>nearestNeighbor(posX,posY)) ) {
        maxPhase = c;
        mPosX = posX+1;
        mPosY = posY;
      }
    }
    if (posX>0 && posY<h-1) {
      c = newPhase[(posY+1)*w+posX-1];
      int nn = nearestNeighbor(posX-1,posY+1);
      if (c>maxPhase && (nn>minSqrDist ^ nn>nearestNeighbor(posX,posY)) ) {
        maxPhase = c;
        mPosX = posX-1;
        mPosY = posY+1;
      }
    }
    if (posY<h-1) {
      c = newPhase[(posY+1)*w+posX];
      int nn = nearestNeighbor(posX,posY+1);
      if (c>maxPhase && (nn>minSqrDist ^ nn>nearestNeighbor(posX,posY)) ) {
        maxPhase = c;
        mPosX = posX;
        mPosY = posY+1;
      }
    }
    if (posX<w-1 && posY<h-1) {
      c = newPhase[(posY+1)*w+posX+1];
      int nn = nearestNeighbor(posX+1,posY+1);
      if (c>maxPhase && (nn>minSqrDist ^ nn>nearestNeighbor(posX,posY)) ) {
        maxPhase = c;
        mPosX = posX+1;
        mPosY = posY+1;
      }
    }
    
    posX = mPosX;
    posY = mPosY;
    
    fill(255,96,0);
    ellipseMode(CENTER);
    ellipse(posX,posY,6,6);
  }
  
  public int nearestNeighbor(int x, int y) {
    int m = 300000000;
    for (int i=0; i<things.length; i++) {
      int sql = (int) (pow(things[i].posX-x,2)+pow(things[i].posY-y,2));
      if (sql < m) m = sql;
    }
    
    return m;
  }
}