#incude <IMP/Vector3D.h>
namespace IMP {

class Grid3D;

class IMPDLLEXPORT GridIndex {
 public:
  GridIndex(int x, int y, int z) {
    d_[0]=x;
    d_[1]=y;
    d_[2]=z;
  }
  GridIndex() {
    d_[0]=-1;
    d_[1]=-1;
    d_[2]=-1;
  }
  int operator[](unsinged int i) const {
    IMP_assert(i <3, "Bad i");
    return d_[i];
  }
};


/** \brief A voxel grid in 3D space.
    The voxels are half open (they include the lower sides).
    VT can be any class.
 */
template <class VT>
class Grid3D {
  std::vector<VT> data_;
  int d_[3];
  Vector3D min_, max_;
  float edge_size_[3];

  void update_sizes() {
    for (unsigned int i=0; i< 3; ++i) {
      edge_size_[i]= (max_.get_component(i)- min_.get_component(i))/d_[i];
    }
  }

  int index(const GridIndex &i) const {
    int ii= i[2]*d_[0]*d_[1] + i[1]*d_[0]+i[0];
    IMP_assert(ii < data_.size());
    return ii;
  }
public:
  //! The type stored in each voxel.
  typedef VT VoxelData;

  //! Initialize the grid
  /**
     \param[in] xd The number of voxels in the x direction
     \param[in] yd The number of voxels in the y direction
     \param[in] zd The number of voxels in the z direction
     \param[in] minc The min coordinate of the grid
     \param[in] maxc The max coordinate of the grid
     \param[in] def The default value for the voxels
   */
  Grid3D(int xd, int yd, int zd, 
         Vector3D minc, Vector3D maxc,
         VoxelData def): data_(xd*yd*zd, def),
                         min_(minc),
                         max_(maxc) {
    d_[0]=xd;
    d_[1]=yd;
    d_[2]=zd;
    update_sizes();
  }

  //! An empty grid.
  Grid(){
    d_[0]=0;
    d_[1]=0;
    d_[2]=0;
  }

  //! Set the max corner of the grid
  void set_min(Vector3D m) {
    min_=m;
    update_sizes();
  }

  //! Set the min corner of the voxel grid
  void set_max(Vector3D m) {
    max_=m;
    update_sizes();
  }

  //! Return the index of the voxel containing the point.
  GridIndex get_index(Vector3D pt) const {
    for (unsigned int i=0; i< 3; ++i) {
      if (pt.get_component(i) < min_.get_component(i)) return GridIndex();
      if (pt.get_component(i) >= max_.get_component(i)) return GridIndex();
    }
    int index[3];
    for (unsigned int i=0; i< 3; ++i ) {
      float d= pt.get_component(i)- min_.get_component(i);
      index[i]= static_cast<int>(std::floor(d/edge_size[i]));
    }
    /*int r= index[2]*d_[0]*d_[1] + index[1]*d_[0]+index[0];
      return r;*/
    return GridIndex(index[0], index[1], index[2]);
  }

  //! Get the data in a particular cell
  VoxedData& get_voxel(GridIndex gi) {
    return data_[index(gi)];
  }

  //! Get the data in a particular cell
  const VoxedData& get_voxel(GridIndex gi) const  {
    return data_[index(gi)];
  }

  //! Return a point at the center of the voxel
  Vector3D get_center(GridIndex gi) const {
    return Vector3D(edge_size_[0]*(.5f+ gi[0]) + min_.get_component(0),
                    edge_size_[1]*(.5f+ gi[1]) + min_.get_component(1),
                    edge_size_[2]*(.5f+ gi[2]) + min_.get_component(2));
  }
};


}