/**==========================================================
 *
 * \file OpenMeshVoronoicellInterface.hpp
 * \brief Creates a Voro++ Voronoi cell from OpenMesh, cut it
 * and return the cutted cell to OpenMesh.
 *
 * \author Leblanc Christophe.
 * \date 02/06/2014
 * \email cleblancad@gmail.com
 *
 * DISCLAIMER:
 *
 * THERE IS NO WARRANTY FOR THE PROGRAM, TO THE EXTENT PERMITTED BY APPLICABLE LAW. 
 * EXCEPT WHEN OTHERWISE STATED IN WRITING THE COPYRIGHT HOLDERS AND/OR OTHER PARTIES 
 * PROVIDE THE PROGRAM “AS IS” WITHOUT WARRANTY OF ANY KIND, EITHER EXPRESSED OR IMPLIED, 
 * INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND 
 * FITNESS FOR A PARTICULAR PURPOSE. THE ENTIRE RISK AS TO THE QUALITY AND PERFORMANCE 
 * OF THE PROGRAM IS WITH YOU. SHOULD THE PROGRAM PROVE DEFECTIVE, YOU ASSUME 
 * THE COST OF ALL NECESSARY SERVICING, REPAIR OR CORRECTION.
 *
 *==========================================================*/
 	
#ifndef INTERSECTION_OPENMESH_VORONOICELL_INTERFACE_HPP
#define INTERSECTION_OPENMESH_VORONOICELL_INTERFACE_HPP
 	
#include <vector>
#include <map>
 	
#include "openMesh/openMesh.hh"
#include "voropp/voro++.hh"
  	
using namespace voro;
 	
namespace intersection {
 	
/** \class OpenMeshVoronoicellInterface.
  * \brief Creates a Voro++ Voronoi cell from OpenMesh.
  *
  * Uses voro++ and OpenMesh.
  *
  * \ingroup voronoicell.
  */
template<class MeshTypeT>
class OpenMeshVoronoicellInterface: public voronoicell
{
public:
  typedef voronoicell Superclass; //!< Parent class.
  typedef MeshTypeT   MeshType;   //!< Type of mesh.
  	
  /// \brief Constructor.
  /// \param[in] mesh reference on mesh.
  OpenMeshVoronoicellInterface(const MeshType &mesh);
 	
  /// \brief Copy constructor.
  /// \param[in] rhs constant reference to OpenMeshVoronoicellInterface object.
  OpenMeshVoronoicellInterface(const OpenMeshVoronoicellInterface &rhs);
  	
  /// \brief Destructor.
  virtual ~OpenMeshVoronoicellInterface();
  	
  /// \brief Copy operator.
  /// \param[in] rhs constant reference to OpenMeshVoronoicellInterface object.
  /// \return Reference to OpenMeshVoronoicellInterface object.
  OpenMeshVoronoicellInterface& operator=(const OpenMeshVoronoicellInterface &rhs);
 	
  /// \brief Create a voronoicell from a mesh.
  void Create();
 	
  /// \brief Update the mesh accordingly to the voronoicell.
  /// \warning Properties of cutted meshes are lost!
  void Update();
  	
    /// \brief Return a voronoicell.
  /// \return Reference on voronoicell.
  voronoicell& GetVoronoiCell()
  { return (*this); }
 	
  /// \brief Return a voronoicell.
  /// \return Constant reference on voronoicell.
  const voronoicell& GetVoronoiCell() const
  { return (*this); }
 	
  /// \brief Return the mesh.
  /// \return Constant reference on mesh.
  const MeshType& GetMesh() const
  { return *m_Mesh; }
 	
  /// \brief Return the mesh.
  /// \return Mesh by value.
  MeshType GetMesh()
  { return *m_Mesh; }
 	
protected:
  MeshType *m_Mesh;      //!< The mesh.
  bool      m_Allocated; //!< True if memory allocated for voronoicell.
 	
private:
  /// \brief Check if there is enough memory for the current cell. Increase if if necessary.
  void CheckMemory();
 	
  /// \brief Get maximum vertex order of mesh.
  /// \return Maximum vertex order of the mesh.
  inline int GetMaximumMeshVertexOrder() const;
  	
  /// \brief Get the number of vertex of order p.
  /// \param[in] p vertex order.
  /// \return Number of vertex of order p.
  inline int GetNumberOfVertexOfOrderP(int p) const;
 	
  /// \brief Doubles the maximum allowed vertex order.
  void AddMemoryVorder();
  	
  /// \brief Increase the memory storage for a particular vertex order.
  void AddMemory(int i, int *stackp2);
 	
  /// \brief Doubles the maximum number of vertices allowed.
  void AddMemoryVertices();
  	
  /// \brief Reset visited edges to positive values.
  inline void ResetEdges();
};
  	
//
// Public functions.
//
 	
template<class MeshTypeT>
OpenMeshVoronoicellInterface<MeshTypeT>::OpenMeshVoronoicellInterface(const MeshType &mesh):
  voronoicell(), m_Allocated(false)
{ m_Mesh = new MeshType(mesh); }
 	
template<class MeshTypeT>
OpenMeshVoronoicellInterface<MeshTypeT>::OpenMeshVoronoicellInterface(
  const OpenMeshVoronoicellInterface &rhs): voronoicell()
{
  m_Allocated = rhs.m_Allocated;
  m_Mesh = new MeshType(*rhs.m_Mesh);
 	
  // voronoicell does not have proper copy constructor.
  const voronoicell *vconst = static_cast<const voronoicell*>(&rhs);
  voronoicell *v = const_cast<voronoicell*>(vconst);
  voronoicell::operator=(*v);
}
 	
template<class MeshTypeT>
OpenMeshVoronoicellInterface<MeshTypeT>::~OpenMeshVoronoicellInterface()
{ delete m_Mesh; }
 	
template<class MeshTypeT>
OpenMeshVoronoicellInterface<MeshTypeT>& OpenMeshVoronoicellInterface<MeshTypeT>::
operator=(const OpenMeshVoronoicellInterface &rhs)
{
  if(this != &rhs) {
    m_Allocated = rhs.m_Allocated;
 	
    delete m_Mesh; m_Mesh = NULL;
    if(rhs.m_Mesh != NULL)
      m_Mesh = new MeshType(*rhs.m_Mesh);
  	
    // voronoicell does not have proper copy operator.
    const voronoicell *vconst = static_cast<const voronoicell*>(&rhs);
    voronoicell *v = const_cast<voronoicell*>(vconst);
    voronoicell::operator=(*v);
  }
  	
  return (*this);
}
 	
template<class MeshTypeT>
void OpenMeshVoronoicellInterface<MeshTypeT>::Create()
{
  // First, check if there is enough memory for the current cell.
  // Increase if if necessary.
  this->CheckMemory();
  int i, j;
  	
  // Set the number of vertex in cell.
  this->p = m_Mesh->n_vertices();
 	
  // Starting vertex index.
  this->up = 0;
 	
  // Set the number of vertex of order i.
  for(i = 0; i < this->current_vertex_order; i++)
    this->mec[i] = this->GetNumberOfVertexOfOrderP(i);
 	
  // Initialize the order of each vertex to zero.
  for(i = 0; i < this->current_vertices; i++)
    this->nu[i] = 0;
  	
  std::map<int, typename MeshType::VertexHandle> vertexHandleMap;
 	
  // Set order of vertex i and vertex coordinates.
  {
    typename MeshType::ConstVertexIter cv_it = m_Mesh->vertices_begin(),
      cv_ite = m_Mesh->vertices_end();
 	
    for(; cv_it != cv_ite; ++cv_it, ++i) {
      i = cv_it.handle().idx();
      vertexHandleMap[i] = cv_it.handle();
      this->nu[i] = m_Mesh->valence(cv_it.handle());
 	
      typename MeshType::Point point = m_Mesh->point(cv_it.handle());
      this->pts[3*i]   = point[0] * 2.0;  // Voro++ divides points' coordinates by 2.
      this->pts[3*i+1] = point[1] * 2.0;
      this->pts[3*i+2] = point[2] * 2.0;      
    }
  }
 	
  // Set the vertex connectivity.
  std::vector<int> vectorCount;
  vectorCount.assign(this->current_vertex_order, 0);
  j = 0;
  	
  for(i = 0; i < this->current_vertices; i++) {
    if(nu[i] > 0) {
      int *q = mep[nu[i]];
 	
      int k = 0;
      for(typename MeshType::ConstVertexVertexIter itVertexVertex = m_Mesh->cvv_iter(
	    vertexHandleMap[i]); itVertexVertex; ++itVertexVertex) {
	q[vectorCount[nu[i]]+k] = itVertexVertex.handle().idx();
	k++;
      }
 	
      ed[j] = q + vectorCount[nu[i]];
      vectorCount[nu[i]] += 2*nu[i]+1;
      j++;
    }
  }
  	
  for(i = 0; i < p; i++) {
    for(j = 0; j < nu[i]; j++) {
      // Determine ed[i][j+nu[i]] such that ed[ ed[i][j] ] [ ed[i][j+nu[i]] ] = i.
      int k1 = ed[i][j];
      int k2 = 0;
 	
      while(ed[k1][k2] != i) {
	k2++;
      }
  	
      ed[i][nu[i]+j] = k2;
    }
  	
    // Set the back pointer.
    ed[i][2*nu[i]] = i;
  }
 	
  m_Allocated = true;
}
 	
template<class MeshTypeT>
void OpenMeshVoronoicellInterface<MeshTypeT>::Update()
{
  if(!m_Allocated)
    voro_fatal_error("voronoicell not allocated in vorop::Update()", VOROPP_MEMORY_ERROR);
  	
  // Temporary copy of mesh for copying its properties.
  MeshType tempMesh = *m_Mesh;
 	
  /// Delete mesh.
  delete m_Mesh;
  m_Mesh = new MeshType();
 	
  // Add new vertex.
  typedef typename MeshType::VertexHandle VertexHandleType;
  int i, j, k, l, m;
  double *ptsp = (pts);
  std::vector< VertexHandleType > vhandles;
  vhandles.resize(p);
 	
  // Code inspired (copied) from voro++-0.4.5/cell.cc line 1542.
  for(i = 0; i < p; i++, ptsp += 3)
    vhandles[i] = m_Mesh->add_vertex(typename MeshType::Point(
      static_cast<float>(ptsp[0]*0.5), static_cast<float>(ptsp[1]*0.5), 
	  static_cast<float>(ptsp[2]*0.5)));
 	
  // Add new faces.
  // Code inspired (copied) from voro++-0.4.5/cell.cc line 1839.
  std::vector< VertexHandleType > face_vhandles;
 	
  for(i = 1; i < p; i++) {
    for(j = 0; j < nu[i]; j++) {
      k = ed[i][j];
 	
      if(k >= 0) {
	face_vhandles.push_back(vhandles[i]);
	ed[i][j] = -1-k;
	l = cycle_up(ed[i][nu[i]+j], k);
  	
	do {
	  face_vhandles.push_back(vhandles[k]);
	  m = ed[k][l];
	  ed[k][l] = -1-m;
	  l = cycle_up(ed[k][nu[k]+l], m);
	  k = m;
	} while(k != i);
 	
	m_Mesh->add_face(face_vhandles);
	face_vhandles.clear();
      }
    }
  }
  	
  this->ResetEdges();
}
 	
//
// Private functions.
//
  	
template<class MeshTypeT>
void OpenMeshVoronoicellInterface<MeshTypeT>::CheckMemory()
{
  if(m_Mesh == NULL)
    voro_fatal_error("NULL mesh in vorop::CheckMemory()", VOROPP_MEMORY_ERROR);
  	
  const int maxVertexOrder = this->GetMaximumMeshVertexOrder();
  while(this->current_vertex_order < maxVertexOrder)
    this->AddMemoryVorder();
 	
  for(int i = 0; i < (this->current_vertex_order); i++) {
    while(mem[i] < (this->GetNumberOfVertexOfOrderP(i)))
      this->AddMemory(i, ds2);
  }
 	
  while(this->current_vertices < static_cast<int>(m_Mesh->n_vertices()))
    this->AddMemoryVertices();
}
 	
/* Increases the memory storage for a particular vertex order, by increasing
 * the size of the of the corresponding mep array. If the arrays already exist,
 * their size is doubled; if they don't exist, then new ones of size
 * init_n_vertices are allocated. The routine also ensures that the pointers in
 * the ed array are updated, by making use of the back pointers. For the cases
 * where the back pointer has been temporarily overwritten in the marginal
 * vertex code, the auxiliary delete stack is scanned to find out how to update
 * the ed value.
 */
template<class MeshTypeT>
void OpenMeshVoronoicellInterface<MeshTypeT>::AddMemory(int i, int *stackp2)
{
  // Code copied from voro++0.4.5 cell.cpp line 129.
  int s = (i << 1) + 1;
 	
  if(mem[i] == 0) {
    this->n_allocate(i, init_n_vertices);
    mep[i] = new int[init_n_vertices*s];
    mem[i] = init_n_vertices;
  }
  else {
    int k, *l, j = 0;
    mem[i] <<= 1;
 	
    if(mem[i] > max_n_vertices)
      voro_fatal_error("Point memory allocation exceeded absolute maximum", VOROPP_MEMORY_ERROR);
 	
    l = new int[s*mem[i]];
    int m = 0;
    this->n_allocate_aux1(i);
  	
    while(j < s*mec[i]) {
      k = mep[i][j + (i << 1)];
  	
      if(k >= 0) {
	ed[k] = l+j;
	this->n_set_to_aux1_offset(k, m);
      }
      else {
	int *dsp;
	for(dsp = ds2; dsp < stackp2; dsp++) {
	  if(ed[*dsp] == mep[i] + j) {
	    ed[*dsp] = l+j;
	    this->n_set_to_aux1_offset(*dsp, m);
	    break;
	  }
	}
 	
	if(dsp == stackp2)
	  voro_fatal_error("Couldn't relocate dangling pointer", VOROPP_INTERNAL_ERROR);
      }
  	
      for(k = 0; k < s; k++, j++) l[j] = mep[i][j];
      for(k = 0; k < i; k++, m++) this->n_copy_to_aux1(i, m);
    }
 	
    delete[] mep[i];
    mep[i] = l;
    this->n_switch_to_aux1(i);
  }
}
 	
template<class MeshTypeT>
inline int OpenMeshVoronoicellInterface<MeshTypeT>::
GetMaximumMeshVertexOrder() const
{
  int count(0), current_count;
  typename MeshType::ConstVertexIter v_it = m_Mesh->vertices_begin(),
    v_ite = m_Mesh->vertices_end();
 	
  for(; v_it != v_ite; ++v_it) {
    current_count = m_Mesh->valence(v_it.handle());
  	
    if(count < current_count)
      count = current_count;
  }
 	
  return count;
}
 	
template<class MeshTypeT>
inline int OpenMeshVoronoicellInterface<MeshTypeT>::
GetNumberOfVertexOfOrderP(int p) const
{
  int count(0);
  typename MeshType::ConstVertexIter v_it = m_Mesh->vertices_begin(),
    v_ite = m_Mesh->vertices_end();
 	
  for(; v_it != v_ite; ++v_it) {
    if(m_Mesh->valence(v_it.handle()) == static_cast<unsigned int>(p))
      count++;
  }
 	
  return count;
}
  	
template<class MeshTypeT>
void OpenMeshVoronoicellInterface<MeshTypeT>::AddMemoryVorder()
{
  // Code copied from voro++-0.4.5 line 204.
  int i = (current_vertex_order << 1), j, *p1, **p2;
 	
  p1 = new int[i];
  for(j = 0; j < current_vertex_order; j++) p1[j] = mem[j]; while(j < i) p1[j++] = 0;
  delete[] mem ;mem = p1;
 	
  p2 = new int*[i];
  for(j = 0; j < current_vertex_order; j++) p2[j] = mep[j];
  delete[] mep; mep = p2;
  	
  p1 = new int[i];
  for(j = 0; j < current_vertex_order; j++) p1[j] = mec[j]; while(j < i) p1[j++] = 0;
  delete[] mec; mec = p1;
  this->n_add_memory_vorder(i);
  current_vertex_order = i;
}
 	
template<class MeshTypeT>
void OpenMeshVoronoicellInterface<MeshTypeT>::AddMemoryVertices()
{
  // Code copied from voro++-0.4.5 line 178.
  int i = (current_vertices << 1), j, **pp, *pnu;
  	
  double *ppts;
  pp = new int*[i];
  for(j = 0; j < current_vertices; j++) pp[j] = ed[j];
  delete[] ed; ed = pp;
  this->n_add_memory_vertices(i);
  	
  pnu = new int[i];
  for(j = 0; j < current_vertices; j++) pnu[j] = nu[j];
  delete[] nu; nu = pnu;
  ppts = new double[3*i];
  for(j = 0; j < 3*current_vertices; j++) ppts[j] = pts[j];
  delete[] pts;pts = ppts;
  current_vertices = i;
}
 	
/* Several routines in the class that gather cell-based statistics internally
 * track their progress by flipping edges to negative so that they know what
 * parts of the cell have already been tested. This function resets them back
 * to positive. When it is called, it assumes that every edge in the routine
 * should have already been flipped to negative, and it bails out with an
 * internal error if it encounters a positive edge. 
 *
 * Code copied from voro++-0.4.5/cell.cc line 1572.
 */
template<class MeshTypeT>
inline void OpenMeshVoronoicellInterface<MeshTypeT>::
ResetEdges()
{
  int i, j;
  for(i = 0; i < p; i++) {
    for(j = 0; j < nu[i]; j++) {
      ed[i][j] = -1-ed[i][j];
    }
  }
}
 	
} // namespace intersection.
  	
#endif // INTERSECTION_OPENMESH_VORONOICELL_INTERFACE_HPP