#ifndef _GRID_H
#define _GRID_H

#include <math.h>
#include <vector>
#include "numerics.h"
#include "./Surface.h"
#include <iostream>
#include <fstream>

using std::ofstream;
using std::cerr;


/*
CLASS
  MoleculeGrid

  This class provides discrete representation of the molecule.
  The molecule is represented by a grid.
  Each voxel of the grid is marked as follows:
  <UL>
  <LI> Negative value for the inside of the molecule
  <LI> Zero for the surface
  <LI> Positive for outside voxels
  </UL>
  In addition the distance function can be calculated for each voxel.
  The distances are zero at the surface voxels and change as the distance
  from the surface increases/decreases. The distance function is negative
  for inside molecule voxels and positive for outside voxels.
  
KEYWORDS
  Grid, Molecule, Surface

AUTHORS
  Yuval Inbar <inbaryuv@post.tau.ac.il>
  Dina Duhovny <duhovka@post.tau.ac.il>,

  Copyright: SAMBA group, Tel-Aviv Univ. Israel, 2001.

CHANGES LOG
<UL>
<LI>
12/01/04 - Dina
added method for volume computation
</UL>

GOALS
  The class is designed to represent a molecule as accurate as possible
  using small amount of space and fast algorithms for grid calculations.

  Distance function calculation:
    initialization:
      The grid voxels are initialized to MAX_VALUE.
      Then the surface voxels are marked with zero.
    The algorithm used for distance function computation is grass-fire algorithm.
    The analogy is that a fire spreads out from some voxels (surface voxels in this case),
    igniting each voxel in the fire's path.
    The distance function is computed for each layer starting from
    surface voxels. The radius for distance function computation is specified
    by the user. There is also a possibility to compute exact distance function:
    each voxel will have the exact distance to the closest surface point.
    Since the algorithm is iterative, the time complexity is linear to the
    number of voxels.

  Marking the inside of the molecule:
    The same algorithm of grass-fire is used here.
    The signs of the voxel distance function is changed layer by layer, starting
    from the centers of molecule atoms. The procedure stops when it reaches
    surface voxels (distance function is zero or very small).
    
  Space Complexity of the Grid: 
    The density of the grid: delta, is specified by the user.
    Assume that the width of the molecule along each of the axis is:
    x_width, y_width, z_width. The space used by the grid will be:
    x_width*y_width*z_width/(delta^3).

  Volume function and Volume normal:
    The class also supports volume function and normal calculation for a point
    on molecule surface.
    The volume function is the ratio of voxels outside the
    molecule in the sphere of radius R (in most articles it is suggested to use
    radius of 5A~), and the total number of voxels in the sphere.
    The volume normal is the average of the points outside the sphere of radius R.
    This normal is more accurate than the one given by MS or Shuo program.
    
USAGE

  The MoleculeGrid class uses Surface representation of the molecule
  generated by Connolly MS program. The density of the points needed
  for dense grid is 5A or more. The density of the grid (delta) should
  be  0.5-1.0. As the density of Connolly points increases, you can decrease
  the delta value of the grid.
  
  EXAMPLE

  Surface msSurface;
  Molecule<Atom> mol;
  // ... read msSurface
  // ... read mol
  // initialization of grid
  MoleculeGrid* grid = new MoleculeGrid(msSurface, delta, maxDist) ;
  // distance function calculation
  grid->computeDistFromSurface(msSurface);
  // marking the molecule
  grid->markTheInside(mol);

  Vector3 point(1.2, 2.3, 3.4);
  float vf = grid->calcVolumeFunc(point,vfRadius) ;
  Vector3 normal = grid->calcVolumeNormal(point,vfRadius) ;
   
  END
*/

class MoleculeGrid
{
 public:
  // GROUP: Constructors
  //// The MoleculeGrid constructor receives Connoly MS Surface, delta - grid
  // density and maxRadius - maximum radius for distance function calculation.
  MoleculeGrid(const Surface& surface, const float inDelta, const float maxRadius);

  // GROUP: Functions
  //// Returns distance from surface for a given point.
  float getDist(const Vector3& point) const {
    int index = getIndexForPoint(point);
    if (!isValidIndex(index))
      return MAX_FLOAT;
    return grid[index];
  }

  //// Calculates volume normal for Vector3 center in a given radius.
  Vector3 calcVolumeNormal(const Vector3& center, const float radius = 5.0) const;

  //// Calculates volume function for Vector3 center in a given radius.
  float calcVolumeFunc(const Vector3& center, const float radius = 5.0) const;

  //// Calculates volume function and normal at one pass.
  // Usefull when both of them need to be calculated.
  void calcVolumeFuncAndNormal(const Vector3& center, const float radius,
			       float& volFunc, Vector3& volNormal) const;
  
  void calcVolumeFuncAndNormal(const Vector3& center,
			       float& volFunc, Vector3& volNormal) const {
    calcVolumeFuncAndNormal(center, 5.0, volFunc, volNormal);
  }
  //// Computes distance function for the molecule.
  // The voxels corresponding to surface points are given zero distance
  // function value and serve as a first layer.
  // If precise value is true, than the exact dist. function is calculated.
  void computeDistFromSurface(const Surface& surface, bool precise = false);

  //// Marks the inside of the molecule by changing the sign of
  // the distance function for inner voxels to negative.
  // the molecule M shouls contain the coordinates of the centers of all atoms
  // of the molecule.
  template< class MoleculeT>
  void markTheInside(const MoleculeT& M);

  //// Computes the volume of the molecule based on grid
  float volume() const;
  
  
  //// Prints each surface voxel as a coordinate in PDB format.
  // Usefull for grid visualization.
  void printGrid(ofstream& file) const;

 protected:
  //// computes index in grid for a point.
  int getIndexForPoint( const Vector3& point) const {
    if (xMin > point[0] || yMin > point[1] || zMin > point[2] || xMax < point[0] || yMax < point[1] || zMax < point[2])
      return -1;
    return ((int)( (point[0] - xMin ) / delta + 0.5)) +
      xGridNum * ((int) ( (point[1] - yMin ) / delta + 0.5)) +
      xyGridNum * ((int) ( (point[2] - zMin ) / delta + 0.5)) ;    
  }

  //// computes a point corresponding to an index.
  Vector3 getPointForIndex(int index) const {
    int x = index % xGridNum;
    index /= xGridNum;
    int y = index % yGridNum;
    int z = index / yGridNum;
    return Vector3(x * delta + xMin, y * delta + yMin, z * delta + zMin);
  }

  //// checks if the index is valid for the grid
  bool isValidIndex(int index) const { return (index >= 0 && index < maxEntry); }

  //// convert radius to grid valid radius (number of voxels)
  int getIntGridRadius( float radius ) const { return (int) (radius / delta + 0.5) ; }

 private:
  // computes xMin,yMin,zMin,xMax,yMax,zMax
  void computeBoundingValues(const Surface& surface);

  // sets distance func. falue for index to be the distance from p.
  bool setDIST(int index, Vector3& p) {
    float d = getPointForIndex(index).dist(p);
    if (d >= grid[index]) return false;
    grid[index] = d;
    return true;
  }
  
 protected:
  //// grid resolution
  float delta ;
  
  int maxGridRadius ;
  int maxEntry ;

  int xGridNum ;
  int yGridNum ;
  int xyGridNum ;
  int zGridNum ;
  float xMin ;
  float yMin ;
  float zMin ;
  float xMax ;
  float yMax ;
  float zMax ;

  //// cube neighbors
  // there are 26 neighbors for each voxel
  vector < int > neigbor;

  //// cube neighbors' distances 
  vector < float > neigborDist;
  
  //// distances grid
  vector < float > grid;
} ;

template< class MoleculeT>
void MoleculeGrid::markTheInside(const MoleculeT& M)
{
  vector< int> layer1, layer2, *curr = &layer1, *next = &layer2, *tmp;
  const float SURFACE_THRESHOLD = delta * sqrt(3.0);
  layer1.reserve(100);
  layer2.reserve(100);
  for (typename MoleculeT::const_iterator it = M.begin(); it!= M.end(); ++it) {
    int index = getIndexForPoint(it->position());
    if (isValidIndex(index) && grid[index] > 0) {
      layer1.push_back(index);
      grid[index] = -grid[index];
    } else {
      cerr << "isValidIndex(index)" << isValidIndex(index) << " index " << index << " grid[index] " << grid[index] << endl;
    }
  }
  int layer = 0;
  //  cerr << "molecule atoms count: " << curr->size() << endl;
  while(!curr->empty()) {
    layer++;
    cerr << "   markTheInside: in Layer " << layer << " out of " << maxGridRadius << " curr->size() " << curr->size() << endl;
    for (vector<int>::iterator it = curr->begin(); it!= curr->end(); ++it) {
      for (int j = 0; j < 26; j++) {
        int nIndex = *it + neigbor[j];
        if (isValidIndex(nIndex) && grid[nIndex] > 0) {
          if (grid[nIndex] > SURFACE_THRESHOLD) next->push_back(nIndex);
          grid[nIndex] = -grid[nIndex];
        }
      }
    }
    curr->clear();
    tmp = curr; curr = next; next = tmp;
  }
}

#endif