//*********************************
//
// Depth Of Field Simulation v1.0
//
// Author: Mario Klingemann
// mario@quasimondo.com
// http://incubator.quasimondo.com
//
//
// This is the first attempt and
// not yet documented. 
//
//********************************

Convolver c;
BImage[] layers;
int layerCount;
int maxLayer;
int maxRadius;
float fLayerCount;
boolean initDone;
float zmin;
float zmax;
float ax;
float ay;

void setup()
{
  size(200, 200);
  background(255);
  stroke(200);
  smooth();
  c=new Convolver(1);
  layerCount=3;
  maxLayer=layerCount-1;
  fLayerCount=float(layerCount);
  layers=new BImage[layerCount];
  maxRadius=layerCount*2;
  zmin=.94;
  zmax=.05;
  ax=HALF_PI;
  ay=0;
}

void loop()
{

  background(255);

  translate(100,100,-250);
  rotateY(ay);
  rotateX(ax);
  ax+=.1;
  ay+=.03;

  fill(240,230,220);
  box(100);
  fill(230,230,230);
  push();
  translate(0,-150,0);
  box(50);
  pop();
  push();
  translate(0,150,0);
  box(50);
  pop();
  push();
  translate(150,0,0);
  box(50);
  pop();
  push();
  translate(-150,0,0);
  box(50);
  pop();
  push();
  translate(0,0,150);
  box(50);
  pop();
  translate(0,0,-150);
  box(50);
  zBlur();
}

void zBlur(){

  int x,y,i,l,m1,m2,rv,gv,bv;
  float d;
  layers[0]=new BImage(width, height);
  layers[0].pixels=pixels;

  for (i=1;i<layerCount;i++){
    layers[i]=new BImage(width, height);
    System.arraycopy( layers[i-1].pixels,0,layers[i].pixels,0,layers[i-1].pixels.length);
    c.setRadius(i*2);
    c.blur(layers[i]);
  }

  float[] zb = g.zbuffer;

  i=0;
  for (y=0;y<height;y++){
    for (x=0;x<width;x++){
      if (zb[i]==MAX_FLOAT){
        d=((getMaxZ(x,y,maxRadius)-zmin)/zmax)*fLayerCount;
      } else {
        d=((zb[i]-zmin)/zmax)*fLayerCount;
      }
      l=max(0,min(int(d),maxLayer));
      if (l==maxLayer){
        pixels[i]=layers[l].pixels[i];
      } else {
        m1=int (256*(d-floor(d)));
        m2=256-m1;
        int p1=layers[l].pixels[i];
        int p2=layers[l+1].pixels[i];
        rv=((((p1 & 0xff0000)>>16)*m2 + ((p2 & 0xff0000)>>16)*m1) << 8) & 0xff0000;
        gv=((((p1 & 0x00ff00)>>8 )*m2 + ((p2 & 0x00ff00)>>8 )*m1)) & 0x00ff00;
        bv=((((p1 & 0x0000ff)     *m2 +  (p2 & 0x0000ff)     *m1)) >> 8) & 0x0000ff;
        pixels[i]=rv | gv | bv;
      }
      i++;
    }
  }

}

float getMaxZ(int x,int y,int r){
  x=max(0,x-r);
  y=max(0,y-r);
  r*=2;
  int xu=min(x+r,width-1);
  int yu=min(y+r,height-1);

  float  maxZ=0;

  float[] gz=g.zbuffer;
  float gzi;

  for (int ym=y;ym<yu;ym++){
    int i=ym*width+x;
    for (int xm=x;xm<xu;xm++){
      gzi=gz[i];
      maxZ=max(maxZ,gzi==MAX_FLOAT ? maxZ : gzi);
      i++;
    }
  }
  return maxZ;
}

class Convolver{
  int radius;
  int kernelSize;
  int[] kernel;
  int[][] mult;

  Convolver(int sz){
    this.setRadius(sz);
  }

  void setRadius(int sz){

    int i,j;
    sz=min(max(1,sz),248);
    if (radius==sz) return;
    kernelSize=1+sz*2;
    radius=sz;
    kernel= new int[1+sz*2];
    mult=new int[1+sz*2][256];

    int sum=0;
    for (i=1;i<sz;i++){
      int szi=sz-i;
      kernel[sz+i]=kernel[szi]=szi*szi;
      sum+=kernel[szi]+kernel[szi];
      for (j=0;j<256;j++){
        mult[sz+i][j]=mult[szi][j]=kernel[szi]*j;
      }
    }
    kernel[sz]=sz*sz;
    sum+=kernel[sz];
    for (j=0;j<256;j++){
      mult[sz][j]=kernel[sz]*j;
    }
  }

  void blur(BImage img){

    int sum,cr,cg,cb,k;
    int pixel,read,i,ri,xl,yl,yi,ym,riw;
    int[] pix=img.pixels;
    int w=img.width;
    int h=img.height;
    int wh=w*h;

    int r[]=new int[wh];
    int g[]=new int[wh];
    int b[]=new int[wh];

    for (i=0;i<wh;i++){
      ri=pix[i];
      r[i]=(ri&0xff0000)>>16;
      g[i]=(ri&0x00ff00)>>8;
      b[i]=(ri&0x0000ff);
    }

    int r2[]=new int[wh];
    int g2[]=new int[wh];
    int b2[]=new int[wh];

    yi=0;
   for (yl=0;yl<h;yl++){
      for (xl=0;xl<w;xl++){
        cb=cg=cr=sum=0;
        ri=xl-radius;
        for (i=0;i<kernelSize;i++){
          read=ri+i;
          if (read>=0 && read<w){
            read+=yi;
            cr+=mult[i][r[read]];
            cg+=mult[i][g[read]];
            cb+=mult[i][b[read]];
            sum+=kernel[i];
          }
        }
        ri=yi+xl;
        r2[ri]=cr/sum;
        g2[ri]=cg/sum;
        b2[ri]=cb/sum;
      }
      yi+=w;
    }
    yi=0;

    for (yl=0;yl<h;yl++){
      ym=yl-radius;
      riw=ym*w;
      for (xl=0;xl<w;xl++){
        cb=cg=cr=sum=0;
        ri=ym;
        read=xl+riw;
        for (i=0;i<kernelSize;i++){
          if (ri<h && ri>=0){
            cr+=mult[i][r2[read]];
            cg+=mult[i][g2[read]];
            cb+=mult[i][b2[read]];
            sum+=kernel[i];
          }
          ri++;
          read+=w;
        }
        pix[xl+yi]=0xff000000 | (cr/sum)<<16 | (cg/sum)<<8 | (cb/sum);
      }
      yi+=w;
    }
  }
}