#ifndef EVALUATIONTREE
#define EVALUATIONTREE

#include "BasisFunction.h"
#include "ImplicitFunction.h"
#include <stdio.h>

class EvaluationTree : public ImplicitFunction{
  
private:
  class EvaluationCell{
  private:
    //start <= i < middle: evalate at this cell
    //middle <= i < end: evaluate at child cells
    int start, middle, end;
    BasisFunction** bfs;
    float cx, cy, cz;
    float size;
    float call;
    
    EvaluationCell** child;
    
  public:
    EvaluationCell(BasisFunction** bfs, int start, int end,
                   float cx, float cy, float cz, float size){
      this->bfs = bfs;
      this->start = middle = start;
      this->end = end;
      
      this->cx = cx;
      this->cy = cy;
      this->cz = cz;
      this->size = size;
      
      child = NULL;
      
      assignBF();
    }
  
    ~EvaluationCell(){
      if(child != NULL){
        for(int i=0; i<8; i++){
          if(child[i] != NULL)
            delete child[i];
        }
        delete child;
      }
    }
    
    double value(float x, float y, float z){
      double f = 0;
      if(child != NULL){
        float minB[3] = {x-call, y-call, z-call};
        float maxB[3] = {x+call, y+call, z+call};
        for(int i=0; i<8; i++){
          EvaluationCell* c = child[i];
          if(c != NULL && c->touch(minB, maxB))
            f += c->value(x, y, z);
        }
      }
      
      for(int i=start; i<middle; i++)
        f += bfs[i]->value(x, y, z);
      
      return f;
    }
    
    void gradient(float g[3], float x, float y, float z){
    }
    
    void computeCallSize(){
      call = 0.25f*size;
      
      if(child != NULL){
        for(int i=0; i<8; i++){
          if(child[i] == NULL)
            continue;
          
          child[i]->computeCallSize();
          for(int j=child[i]->start; j<child[i]->middle; j++){
            float r = bfs[j]->support;
            if(call < r)
              call = r;
          } 
        }
      }
    }
    
  private:
    bool inline touch(float minB[3], float maxB[3]){
      if(maxB[0] < cx-size || maxB[1] < cy-size || maxB[2] < cz-size ||
         minB[0] > cx+size || minB[1] > cy+size || minB[2] > cz+size)
        return false;
      else
        return true;
    }
    
    void assignBF(){
      for(int i=start; i<end; i++){
        BasisFunction* bfi = bfs[i];
        float r = bfi->support;
        if(2*r > size){
          bfs[i] = bfs[middle];
          bfs[middle++] = bfi;
        }
      }
      
      //all functions are assigned?
      if(middle == end)
        return;
      
      int sZ = splitZ(middle, end, cz);
      
      int sY1 = splitY(middle, sZ, cy);
      int sY2 = splitY(sZ, end, cy);
      
      int sX11 = splitX(middle, sY1, cx);
      int sX12 = splitX(sY1, sZ, cx);
      int sX21 = splitX(sZ, sY2, cx);
      int sX22 = splitX(sY2, end, cx);
      
      child = new EvaluationCell*[8];
      
      float s = 0.5f*size;
      
      if(middle != sX11)
        child[0] = new EvaluationCell(bfs, middle, sX11, cx-s, cy-s, cz-s, s);
      else
        child[0] = NULL;
      
      if(sX11 != sY1)
        child[1] = new EvaluationCell(bfs, sX11, sY1, cx+s, cy-s, cz-s, s);
      else
        child[1] = NULL;
      
      if(sY1 != sX12)
        child[2] = new EvaluationCell(bfs, sY1, sX12, cx-s, cy+s, cz-s, s);
      else
        child[2] = NULL;
      
      if(sX12 != sZ)
        child[3] = new EvaluationCell(bfs, sX12, sZ, cx+s, cy+s, cz-s, s);
      else
        child[3] = NULL;
      
      if(sZ != sX21)
        child[4] = new EvaluationCell(bfs, sZ, sX21, cx-s, cy-s, cz+s, s);
      else
        child[4] = NULL;
      
      if(sX21 != sY2)
        child[5] = new EvaluationCell(bfs, sX21, sY2, cx+s, cy-s, cz+s, s);
      else
        child[5] = NULL;
      
      if(sY2 != sX22)
        child[6] = new EvaluationCell(bfs, sY2, sX22, cx-s, cy+s, cz+s, s);
      else
        child[6] = NULL;
      
      if(sX22 != end)
        child[7] = new EvaluationCell(bfs, sX22, end, cx+s, cy+s, cz+s, s);
      else
        child[7] = NULL;
    }
    
    int splitX(int s, int e, float m){
      int i = s;
      int j = e-1;
      while(i <= j){
        BasisFunction* bfi = bfs[i];
        if(bfi->centerX < m)
          i++;
        else{
          bfs[i] = bfs[j];
          bfs[j--] = bfi;
        }
      }
      return i;
    }
    
    int splitY(int s, int e, float m){
      int i = s;
      int j = e-1;
      while(i <= j){
        BasisFunction* bfi = bfs[i];
        if(bfi->centerY < m)
          i++;
        else{
          bfs[i] = bfs[j];
          bfs[j--] = bfi;
        }
      }
      return i;
    }
    
    int splitZ(int s, int e, float m){
      int i = s;
      int j = e-1;
      while(i <= j){
        BasisFunction* bfi = bfs[i];
        if(bfi->centerZ < m)
          i++;
        else{
          bfs[i] = bfs[j];
          bfs[j--] = bfi;
        }
      }
      return i;
    }
  };
  
private:
  EvaluationCell* root;
  BasisFunction** bfs;
  int bfN;
  
public:
  float max[3], min[3];
  float out;
  
  EvaluationTree(BasisFunction** bfs, int bfN, float min[3], float max[3]){
    this->bfs = bfs;
    this->bfN = bfN;
    
    this->min[0] = min[0];
    this->min[1] = min[1];
    this->min[2] = min[2];
    
    this->max[0] = max[0];
    this->max[1] = max[1];
    this->max[2] = max[2];
    
    float mx = 0.5f*(min[0] + max[0]);
    float my = 0.5f*(min[1] + max[1]);
    float mz = 0.5f*(min[2] + max[2]);
    float s = max[0] - min[0];
    if(s < max[1] - min[1])
      s = max[1] - min[1];
    if(s < max[2] - min[2])
      s = max[2] - min[2];
    s *= 0.5f;
    
    printf("Octree construction is started.\n");
    
    root = new EvaluationCell(bfs, 0, bfN, mx, my, mz, s);
    root->computeCallSize();
    
    printf("Octree construction is Fnished.\n\n");
    
    out = -10;
  }
  
  ~EvaluationTree(){
    for(int i=0; i<bfN; i++)
      delete bfs[i];
    delete bfs;
  }
  
  float value(float x, float y, float z){
    if(isIn(x,y,z))
      return (float)root->value(x, y, z)+out;
    else
      return out;
  }
  
  void gradient(float g[3], float x, float y, float z){
    if(isIn(x,y,z))
      root->gradient(g, x, y, z);
    else
      g[0] = g[1] = g[2] = out;
  }
  
private:
  bool inline isIn(float x, float y, float z){
    if(x < min[0] || y < min[1] || z < min[2] ||
       x > max[0] || y > max[1] || z > max[2])
      return false;
    else
      return true;
  }
};

#endif