// An accumulation bin counts the number of hits in a certain area. This area
// is determined by the range of the accumulation grid and its resolution.

// The accumulation allows fractional counts, for example a hit can result in a 0.25/0.75 distribution over two bins.
// This generalisation helps in reducing sharp pixel artefacts in very high density areas.

// The accumulation grid serves two purposes. It collects an array of bins and
// assigns each a certain physical area.
// The physical limits (e.g. -1.0f to 1.0f) are determined by the attractor system.
// The resolution (e.g. 100 by 100 bins) is defined here.

class AccumulationGrid{

  float[][] count; // bins
  float[][] fastCount; // simple bins without fractiobal hits
  float maxCount;
  float minCount;
  float invLogCount;
  
  GridParameters parameters;

  float idx, idy;// inverse dimensions of bin area

  AccumulationGrid(GridParameters GP){ 
    parameters = GP.get();
    idx=parameters.resX/(parameters.maxVisX-parameters.minVisX);// inverse width of a bin
    idy=parameters.resY/(parameters.maxVisY-parameters.minVisY);// inverse height of a bin
    count = new float[parameters.resX][parameters.resY];
    fastCount = new float[parameters.resX][parameters.resY];
    maxCount=minCount=invLogCount=0f;
  }

  void reset(boolean resized){    
    idx=parameters.resX/(parameters.maxVisX-parameters.minVisX);
    idy=parameters.resY/(parameters.maxVisY-parameters.minVisY);
    if(resized){
         count = new float[parameters.resX][parameters.resY];
         fastCount = new float[parameters.resX][parameters.resY];
         maxCount=minCount=invLogCount=0f;
    }
    else{
      for(int i=0;i<parameters.resX;i++){
        for(int j=0;j<parameters.resY;j++){
          count[i][j]=0f;
          fastCount[i][j]=0f;
        }
      }
    }
    maxCount=minCount=invLogCount=0f;
  }

  // The first important part of the accumulation grid is the add method.
  // It converts physical coordinates (for example (-0.5f,0.3f)) to a range of bins.
  // Depending on the rendermode the added value can be spread over adjacent bins.
  // renderMode 0, "fast" : full value in closest bin.
  // renderMode 1, "bilinear" : bilinear interpolation over 4 bins.
  // renderMode 2: "Gaussian" : Gaussian distribution over several bins. The width is determined by
  //   the local normalized logDensity and the focus parameter.
   

  void add(float x, float y, float c, int quality,  float focus, GaussianConvolutionKernel gck){

    if((x>=parameters.minVisX)&&(x<parameters.maxVisX)&&(y>=parameters.minVisY)&&(y<parameters.maxVisY)){
      float gridfx =(x-parameters.minVisX)*idx;
      float gridfy =(y-parameters.minVisY)*idy;
      int gridX =constrain((gridfx<0f)?1+(int)gridfx:(int)gridfx,0,parameters.resX-1);
      int gridY = constrain((gridfy<0f)?1+(int)gridfy:(int)gridfy,0,parameters.resY-1);

      // Density dependent calculation of the spread
      // For low density areas logCount is very small: large spread.
      // High density areas (count>10^rc) have no spread, only bilinear interpolation
            
     fastCount[gridX][gridY]++;
      
      int range =0;
      if(quality==2){
        if(fastCount[gridX][gridY]<10){// The first few hits in Gaussian mode are used to estimate the local density
          range=-1;   
        }
        else{
          range=max(0, gck.maxRange-(int)(2f*gck.maxRange*FastFunctions.log2(fastCount[gridX][gridY]+1f)*invLogCount));
        }
      }


      if( quality==0){
        count[gridX][gridY]+=c; 
        updateStatBins(gridX,gridY);
      }
      else if(range==0){
        float residuX=gridfx-gridX;
        float residuY=gridfy-gridY;

        count[gridX][gridY]+=(1f-residuX)*(1f-residuY)*c; 
        updateStatBins(gridX,gridY);

        if(gridX+1<parameters.resX){ 
          count[gridX+1][gridY]+=residuX*(1f-residuY)*c; 
          updateStatBins(gridX+1,gridY);
        }
        if(gridY+1<parameters.resY){ 
          count[gridX][gridY+1]+=(1f-residuX)*residuY*c;  
          updateStatBins(gridX,gridY+1);
        }
        if((gridX+1<parameters.resX)&&(gridY+1<parameters.resY)){ 
          count[gridX+1][gridY+1]+=residuX*residuY*c; 
          updateStatBins(gridX+1,gridY+1);
        }
      }
      else if(range==-1){
        //do nothing, except add to fastGrid
      }
      else{
        convAdd(gridX,gridY,c,range,gck.limits[range],gck);
      }
    }
  }

  void convAdd(int X, int Y,float c, int range, int lim,GaussianConvolutionKernel gck){
    for(int i=-lim;i<lim+1;i++){
      for(int j=-lim;j<lim+1;j++){
        if((X+i>=0)&&(X+i<parameters.resX)&&(Y+j>=0)&&(Y+j<parameters.resY)){
          //Normalized 2D Gaussian function, the total cumulative value of a hit stays 1 independent of range.
          //This is necessary to ensure that after long runs there's no overexposure of the low density regions
         count[X+i][Y+j]+=c*gck.v[gck.maxSize+i][gck.maxSize+j][range];
          updateStatBins(X+i,Y+j);
        }
      }
    }
  }

  void updateStatBins(int i, int j){
    boolean recalcLogCount=false;
  
    
    if(count[i][j]>maxCount){
      maxCount=count[i][j];
      recalcLogCount=true;
    }
    
    if(count[i][j]<minCount){
      minCount=count[i][j];
      recalcLogCount=true;
    }
   
    
    if(recalcLogCount)invLogCount=((maxCount-minCount)>0f)?1f/FastFunctions.log2(1f+maxCount-minCount):0f;

  }

  // This is the second important part, used for displaying the grid. It resamples the accumulation grid
  //  to arbitrary resolution, e.g 200x200 bins to 500x500 pixels or 2000x2000 bins to 1000x1000 pixels.
  //  Both downscaling and upscaling are supported. Non-integer ratios,
  //  e.g. 1.5x1.5 bins per pixel are handled by weighing values for partial overlap, half a bin would add half its value. 
  //  Furthermore the range of the returned values should be independent of the requested resolution. For this the returned values
  //  are averaged over the sampling area.
  
  ResultBin getValues(int i, int j, int w, int h){
    float binsPerPixelx = (float)parameters.resX/(float)w;
    float binsPerPixely = (float)parameters.resY/(float)h;  
    float startx=max(0f,i*binsPerPixelx);
    int sx = (int)startx;
    float residusx = 1f-startx+sx;
    float starty=max(0f,j*binsPerPixely);
    int sy = (int)starty;
    float residusy = 1f-starty+sy;
    float endx=startx+binsPerPixelx-0.0001f;
    int ex = (int)endx;
    float residuex = endx-ex;
    float endy=starty+binsPerPixely-0.0001f;
    int ey = (int)endy;
    float residuey = endy-ey; 

    float c=0f; // Running total of count
    float area=0f; // Running total of the number of bins
    float factorx=0f;
    float factory=0f;
    float phi;

    for(int x=sx; x<ex+1;x++){
      factorx=1f;
      if(x==sx) factorx=residusx; //Take into account partial overlap of first column of bin
      if(x==ex) factorx=residuex;//Take into account partial overlap of last column of bins
      for (int y=sy; y<ey+1;y++){
        factory=1f;
        if(y==sy) factory=residusy; //Take into account partial overlap of first row of bins
        if(y==ey) factory=residuey; //Take into account partial overlap of last column of bins
        if((x<parameters.resX)&&(y<parameters.resY)){
          phi=factorx*factory*factorx*factory;
          c+=phi*count[x][y]; 
          area+=phi;                              
        }
      }
    }
    //Weighted average value over involved bins, this preserves the range of the values. 
    c/=area; 

    ResultBin result = new ResultBin();    
    result.rawCount=c;
    result.normCount=(c-minCount)/(maxCount-minCount);
    result.logNormCount=FastFunctions.log2(c-minCount+1f)*invLogCount;
    return result;
  }

}


// Auxiliary structure that contains all values for a pixel
class ResultBin{
  float rawCount;
  float normCount;
  float logNormCount;
 
  ResultBin(){
  }
}