/**==========================================================
 *
 * \file Intersection.hpp
 * \brief Given two convex cells, returns the intersection (if any).
 *
 * \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_INTERSECTION_HPP
#define INTERSECTION_INTERSECTION_HPP
 	
#include <cmath>
#include <stdexcept>
#include "OpenMeshVoronoicellInterface.hpp"
#include "MeshToVTKWriter.hpp"
 	
namespace intersection {
  	
/** \class Intersection.
  * \brief Given two convex cells, returns the intersection (if any).
  * \warning REQUIRED: Normals have to point inside the polytope.
  */
template<typename MeshTypeT>
class Intersection
{
public:
  typedef MeshTypeT   MeshType;   //!< Type of mesh.
 	
  /// \brief Constructor.
  /// \param[in] mesh1: first cell. (Cell to cut).
  /// \param[in] mesh2: second cell. (Cutting cell).
  Intersection(MeshType &mesh1, MeshType &mesh2);
 	
  /// \brief Destructor.
  ~Intersection();
  	
  /// \brief Copy constructor.
  /// \param[in] rhs Intersection object.
  Intersection(const Intersection &rhs);
 	
  /// \brief Copy operator.
  /// \param[in] rhs Intersection object.
  /// \return Copy of Intersection object.
  Intersection& operator=(const Intersection &rhs);
  	
  /// \brief Tell if intersection is empty or not.
  /// \return True if intersection empty, false otherwise.
  inline bool EmptyIntersection() const
  { return !m_IntersectionNotEmpty; }
  	
  /// \brief Return the intersection mesh.
  /// \return Constant reference on mesh.
  inline const MeshType& GetIntersectionMesh() const
  { 
    if(EmptyIntersection()) {
      throw std::runtime_error("Intersection:: error: no intersection to return.");
    }
 	
    return *m_Intersection; 
  }

  /// \brief Return the intersection mesh.
  /// \return Mesh by value.
  inline MeshType GetMesh()
  { 
    if(EmptyIntersection()) {
      throw std::runtime_error("Intersection:: error: no intersection to return.");
    }
  	
    return *m_Intersection; 
  }
 	
  /// \brief Return the intersection mesh.
  /// \return Mesh by value.
  inline MeshType GetIntersectionMesh()
  { 
    if(EmptyIntersection()) {
      throw std::runtime_error("Intersection:: error: no intersection to return.");
    }
  	
    return *m_Intersection; 
  }
 	
  /// \brief Compute the intersection.
  /// \return True if intersection not empty, false otherwise.
  bool Compute();
 	
protected:
  typedef OpenMeshVoronoicellInterface<MeshType> OpenMeshVoronoicellInterfaceType; //!< Interface type betwenn OpenMesh and voronoicell.
 	
  MeshType *m_Mesh1, *m_Mesh2; //!< The two input mesh.
  MeshType *m_Intersection;    //!< The intersection mesh.
  bool m_IntersectionNotEmpty; //!< True if intersection is not empty, false otherwise.
 	
private:
 	
  /// \brief Compute the cutting plane parameters for voro++
  ///  given a mesh and face handle of it.
  /// \param[in] mesh considered mesh.
  /// \param[in] fh face handle.
  /// \param[in] c centroid of the mesh/cell.
  /// \param[out] (x, y, z) plane parameters for voro++.
  void ComputePlaneParameters(const MeshType &mesh, 
    const typename MeshType::FaceHandle fh, const typename OpenMesh::VectorT<double, 3> &c, 
    double &x, double &y, double &z) const;
};
 	
//
// Public functions.
//
  	
template<class MeshTypeT>
Intersection<MeshTypeT>::Intersection(MeshType &mesh1, MeshType &mesh2):
  m_Mesh1(&mesh1), m_Mesh2(&mesh2), m_Intersection(NULL),
  m_IntersectionNotEmpty(false)
{}
 	
template<class MeshTypeT>
Intersection<MeshTypeT>::Intersection(const Intersection &rhs):
  m_Mesh1(rhs.m_Mesh1), m_Mesh2(rhs.m_Mesh2), m_Intersection(NULL),
  m_IntersectionNotEmpty(rhs.m_IntersectionNotEmpty)
{
  if(rhs.m_Intersection != NULL)
    m_Intersection = new MeshType(*rhs.m_Intersection);
}
  	
template<class MeshTypeT>
Intersection<MeshTypeT>& Intersection<MeshTypeT>::
  operator=(const Intersection &rhs)
{
  if(this != &rhs) {
    m_Mesh1 = rhs.m_Mesh1;
    m_Mesh2 = rhs.m_Mesh2;
    m_IntersectionNotEmpty = rhs.m_IntersectionNotEmpty;
 	
    delete m_Intersection; m_Intersection = NULL;
    if(rhs.m_Intersection != NULL)
      m_Intersection = new MeshType(*rhs.m_Intersection);
  }
  	
  return (*this);
}
  	
template<class MeshTypeT>
Intersection<MeshTypeT>::~Intersection()
{ delete m_Intersection; }
 	
 	
template<class MeshTypeT>
bool Intersection<MeshTypeT>::Compute()
{
  // Check if we have face normals for meshes. Otherwise we request them.
  if(!(m_Mesh1->has_face_normals())) {
    m_Mesh1->request_face_normals();
  }
 	
  if(!(m_Mesh2->has_face_normals())) {
    m_Mesh2->request_face_normals();
  }
 	
  // Update face normals.
  m_Mesh1->update_face_normals();
  m_Mesh2->update_face_normals();

  // Compute the voronoi cell 2 from mesh 2.
  OpenMeshVoronoicellInterface<MeshType> cell2(*m_Mesh2);
  cell2.Create();
 	
  // Compute the centroid of cell 2.
  typename OpenMesh::VectorT<double, 3> c2;
  cell2.centroid(c2[0], c2[1], c2[2]);
 	
  // Shift a copy of mesh 1 by c2.
  m_Intersection = new MeshType(*m_Mesh1);
  typename MeshType::VertexIter v_it  = m_Intersection->vertices_begin();
  typename MeshType::VertexIter v_ite = m_Intersection->vertices_end();
  	
  for(; v_it != v_ite; ++v_it) {
    typename MeshType::Point &p = m_Intersection->point(v_it.handle());
    p[0] -= static_cast<float>(c2[0]); 
	p[1] -= static_cast<float>(c2[1]); 
	p[2] -= static_cast<float>(c2[2]);
  }
  	
 	
  // Compute the voronoicell 1 from the mesh copy.
  OpenMeshVoronoicellInterface<MeshType> cell1(*m_Intersection);
  cell1.Create();
  	
  // Cut cell 1 by the faces of cell 2.
  // Loop on faces of cell 2 and cut.
  m_IntersectionNotEmpty = true;
  typename MeshType::ConstFaceIter f_it  = m_Mesh2->faces_begin();
  typename MeshType::ConstFaceIter f_ite = m_Mesh2->faces_end();
 	
  for(; m_IntersectionNotEmpty && (f_it != f_ite); ++f_it) {
    // Compute plane parameters for voro++ for the current face.
    double x, y, z;
    ComputePlaneParameters(*m_Mesh2, f_it.handle(), c2, x, y, z);
 	
    // Perform cut.
    m_IntersectionNotEmpty &= cell1.nplane(x, y, z, f_it.handle().idx());
  }
 	
  delete m_Intersection;
  m_Intersection = NULL;
  	
  if(m_IntersectionNotEmpty) {
    // Get cutted cell, which is the intersection.
    cell1.Update();
    m_Intersection = new MeshType(cell1.GetMesh());
 	
    // Shift mesh backward.
    v_it  = m_Intersection->vertices_begin();
    v_ite = m_Intersection->vertices_end();
 	
    for(; v_it != v_ite; ++v_it) {
      typename MeshType::Point &p = m_Intersection->point(v_it.handle());
 	
      // Rounding.
      for(int i = 0; i < 3; i++)
        if(fabs(p[i]) <= voro::tolerance) p[i] = 0.0;
 	
      // Shift.
      p[0] += static_cast<float>(c2[0]); 
	  p[1] += static_cast<float>(c2[1]); 
	  p[2] += static_cast<float>(c2[2]);
    }
  }
  	
  // Done.
  return m_IntersectionNotEmpty;
}
 	
 	
//
// Private functions.
//
  	
template<class MeshTypeT>
void Intersection<MeshTypeT>::ComputePlaneParameters(const MeshType &mesh, 
  const typename MeshType::FaceHandle fh, const typename OpenMesh::VectorT<double, 3> &c, 
  double &x, double &y, double &z) const
{
  // Get normal to face.
  const typename MeshType::Normal n = (mesh.normal(fh));
 	
  // Compute plane bisectrix.
  double meanScalarProduct(0.0);
  unsigned int nbVerticesInFace(0);
  typename MeshType::ConstFaceVertexIter fv_it = mesh.cfv_iter(fh);
  	
  for(; fv_it; ++fv_it, ++nbVerticesInFace) {
    typename MeshType::Point p = mesh.point(fv_it.handle());
    p[0] -= static_cast<float>(c[0]); 
	p[1] -= static_cast<float>(c[1]); 
	p[2] -= static_cast<float>(c[2]);
    meanScalarProduct += (p | n);
  }
  	
  meanScalarProduct /= static_cast<double>(nbVerticesInFace);
  x = 2.0*meanScalarProduct*n[0];
  y = 2.0*meanScalarProduct*n[1];
  z = 2.0*meanScalarProduct*n[2];
}
 	
} // namespace intersection.
 	
#endif // INTERSECTION_INTERSECTION_HPP