/**============================================================================
 *
 * \file   halfedge.hpp
 * \brief  Declarations for the  Face_decimation class.
 *
 * \author Leblanc Christophe.
 * \date   05/09/2019
 * \email cleblancad@gmail.com
 *
 * \source: "Automatic Generation of Structure Preserving Multiresolution Models",
 *  Marinov, M., Kobbelt, L., Computer Graphics Froum, vol 24, issue 3, pp. 479-486 (2005)
 *  doi: 10.1111/j.1467-8659.2005.00873.x
 *
 *===========================================================================*/

#ifndef MARINOV_FACE_DECIMATION_HPP
#define MARINOV_FACE_DECIMATION_HPP

#include <string>
#include <vector>
#include <map>
#include <set>
#include <list>
#include <stdexcept>
#include <algorithm>
#include <cmath>
#include <utility>
#include <array>

#include <CGAL/Surface_mesh.h>
#include <CGAL/Polygon_mesh_processing/measure.h>
#include <CGAL/Polygon_mesh_processing/compute_normal.h>
#include <CGAL/linear_least_squares_fitting_3.h>
#include <CGAL/Polygon_2.h>
#include <CGAL/Plane_3.h>
#include <CGAL/squared_distance_2.h>

#include <CGAL/Triangulation_vertex_base_with_info_2.h>
#include <CGAL/Constrained_Delaunay_triangulation_2.h>
#include <CGAL/Triangulation_conformer_2.h>

#include <marinov/progress_bar.hpp>


namespace marinov {

template<class KernelT>
class Face_decimation
{
private:
  // Forward declaration.
  struct Face_properties;

public:
  typedef KernelT Kernel;
  typedef typename Kernel::Point_3 Point;
  typedef CGAL::Surface_mesh<Point> Mesh;

  typedef typename Mesh::Vertex_index Vertex_index;
  typedef typename Mesh::Halfedge_index Halfedge_index;
  typedef typename Mesh::Edge_index Edge_index;
  typedef typename Mesh::Face_index Face_index;

  typedef std::map<Face_index, Face_properties> Mesh_properties;


  // Enums.
  enum Boundary_type
  {
    BOUNDARY_UNKNOWN = 0, ///< Do not known if belongs to a boundary or not.
    BOUNDARY_LEFT = 10,
    BOUNDARY_RIGHT = 11,
    BOUNDARY_FRONT = 20,
    BOUNDARY_REAR = 21,
    BOUNDARY_BOTTOM = 30,
    BOUNDARY_TOP = 31,
    BOUNDARY_NONE = 100 ///< Do not belong to a boundary.
  };


private:
  typedef typename Kernel::FT FT;
  typedef typename Kernel::Vector_3 Vector;
  typedef CGAL::Plane_3<Kernel> Plane;

  typedef CGAL::Triangulation_vertex_base_with_info_2<std::size_t, Kernel> Vertex_info;
  typedef CGAL::Triangulation_face_base_2<Kernel> Face_info;
  typedef CGAL::Constrained_triangulation_face_base_2<Kernel, Face_info> Face_base;
  typedef CGAL::Triangulation_data_structure_2<Vertex_info, Face_base> TDS;

  typedef CGAL::Constrained_Delaunay_triangulation_2<Kernel, TDS, CGAL::No_intersection_tag>
    CDT;

  typedef CGAL::Polygon_2<Kernel> Polygon;
  typedef typename Mesh::template Property_map<Face_index, Polygon> Polygon_map;

public:
  Face_decimation(const Mesh &mesh, std::ostream &os = std::cerr);

  void decimate();

  const Mesh& get_mesh() const
  { return m_Mesh; }

  const Mesh& get_compound_mesh() const
  { return m_CompoundMesh; }

  const Mesh_properties& get_mesh_properties() const
  { return m_MeshProperties; }

  const Mesh_properties& get_compound_mesh_properties() const
  { return m_CompoundMeshProperties; }

  void set_target_number_of_faces(const std::size_t n)
  { m_TargetNumberOfFaces = n; }

  void set_triangulate_mesh(const bool tri)
  { m_TriangulateMesh = tri; }

  const std::size_t get_target_number_of_faces() const
  { return m_TargetNumberOfFaces; }

  const std::ostream* get_error_stream() const
  { return m_Os; }

  const bool get_triangulate_mesh() const
  { return m_TriangulateMesh; }

private:
  // Typedefs.
  struct Face_properties
  {
    bool selection;
    bool mergeable;
    Point centroid;
    Vector normal;
    FT area;
    Polygon polygon;
    CDT triangulation;
    Mesh lifted_triangulation;
    FT error;
    Face_index face_to_merge;
    Boundary_type boundary;

    Face_properties()
    { clear(); }

    void clear()
    {
      constexpr FT nan = std::numeric_limits<FT>::quiet_NaN();

      selection = true;
      mergeable = false;
      centroid = Point(nan, nan, nan);
      normal = Vector(nan, nan, nan);
      area = nan;
      polygon.clear();
      triangulation.clear();
      lifted_triangulation.clear();
      error = nan;
      face_to_merge = Mesh::null_face();
      boundary = Boundary_type::BOUNDARY_UNKNOWN;
    }
  };

  // Variables.
  Mesh m_Mesh;
  Mesh m_CompoundMesh;
  Mesh_properties m_MeshProperties;
  Mesh_properties m_CompoundMeshProperties;
  std::size_t m_TargetNumberOfFaces;
  std::ostream* m_Os;
  bool m_TriangulateMesh;

  //
  // Marinov algorithm.
  //

  // Initialization.
  void compute_initial_properties(const Mesh &mesh, Mesh_properties &properties);
  void compute_face_centroid(const Face_index &f, const Mesh &mesh, Mesh_properties &properties);
  void compute_2d_polygon(const Face_index &f, const Mesh &mesh,
    const Mesh_properties &properties, Polygon &polygon, CDT &cdt);
  void lift_triangulation(const Face_index &f, const Mesh_properties &properties,
    const Mesh &mesh, Mesh &mesh_t) const;
  void compute_merging_errors(const Mesh &mesh, Mesh_properties &properties);
  void compute_merging_errors(const Mesh &mesh, const Face_index &f,
    Mesh_properties &properties);
  FT compute_merging_error(const Mesh &mesh, const Mesh_properties &properties,
    const Face_index &f0, const Face_index &f1, bool &mergeable) const;
  FT compute_l2_error(const Mesh &mesh, const Plane &plane) const;
  std::vector<Vertex_index> merge_face(const Mesh &mesh,
    const Face_index &f0, const Face_index &f1, const Vector &normal, std::vector<FT> &alpha) const;
  bool injectivity_condition(const Polygon &polygon) const;

  // Computation.
  bool compute_next_mesh(Mesh &mesh, Mesh_properties &properties);
  bool compute_next_merged_face(Mesh &mesh, const Mesh_properties &properties,
    Face_index &face_to_merge, Face_index &new_face);
  void update_properties_after_face_merge(const Mesh &mesh, Mesh_properties &properties,
    const Face_index &face_to_merge, const Face_index &new_face);

  // Vertex collapse.
  void seek_and_destroy_vertices_of_valence2(Mesh &mesh, Mesh_properties &properties);
  bool check_vertex_collapse(const Mesh &mesh, const Vertex_index &v0,
    const Mesh_properties &properties);
  Mesh extract_mesh_faces(const Mesh &mesh, const Halfedge_index &h, Halfedge_index &new_h,
    std::map<Face_index, Face_index> &faces_correspondance) const;
  void compute_vertex_collapse(Mesh &mesh, const Vertex_index &vk,
    Face_index &left_face, Face_index &right_face);
  void update_properties_after_vertex_collapse(const Mesh &mesh, const Vertex_index &vk,
    const Face_index &f0, const Face_index &f1, Mesh_properties &properties);

  // Vertex removal.
  bool check_vertex_removal(const Mesh &mesh, const Vertex_index &vk,
    const Mesh_properties &properties) const;
  void compute_vertex_removal(Mesh &mesh, const Vertex_index &vk, Face_index &other_face,
    Face_index &new_face);
  void update_properties_after_vertex_removal(const Mesh &mesh, const Vertex_index &vk,
    const Face_index &other_face, const Face_index &new_face, Mesh_properties &properties);

  // Boundary identification.
  void identify_faces_at_boundaries(const Mesh &mesh, Mesh_properties &properties,
    std::array<FT, 3> &min_c, std::array<FT, 3> &max_c);

  // Compound mesh.
  void compute_compound_mesh(const Mesh &mesh, const Mesh_properties &properties,
    Mesh &triangulated_mesh);

  // The machine epsilon has to be scaled to the magnitude of the values used
  // and multiplied by the desired precision in ULPs (units in the last place)
  // Source: https://en.cppreference.com/w/cpp/types/numeric_limits/epsilon
  /// \brief Return true if x and y are "nearly" equal, false otherwise.
  template<class T>
  inline bool nearly_equal(const T &x, const T &y, const T factor = 1.0, const int ulp = 2) const
  {
    return (std::abs(x-y) <= std::numeric_limits<T>::epsilon() * factor * std::abs(x+y) * ulp)
      // unless the result is subnormal
           || (std::abs(x-y) < std::numeric_limits<T>::min());
  }
};


template<typename KernelT>
Face_decimation<KernelT>::Face_decimation(const Mesh &mesh, std::ostream &os):
  m_Os(&os), m_TriangulateMesh(true)
{
  m_Mesh = mesh;
}


template<typename KernelT>
void Face_decimation<KernelT>::decimate()
{
  std::array<FT, 3> min_c, max_c;
  identify_faces_at_boundaries(m_Mesh, m_MeshProperties, min_c, max_c);

  (*m_Os) << "\nFace decimation with " << m_TargetNumberOfFaces << " target number of faces, "
          << "starting from: " << m_Mesh.number_of_faces() << " faces.\n";

  (*m_Os) << "\nCompute initial face properties...\n";
  compute_initial_properties(m_Mesh, m_MeshProperties);

  std::size_t number_of_faces = m_Mesh.number_of_faces();
  const double a = m_TargetNumberOfFaces /
                   (static_cast<double>(m_TargetNumberOfFaces) - number_of_faces);
  const double b = -a*number_of_faces;

  ProgressBar progress_bar(*m_Os, 50);
  (*m_Os) << "\nFace decimation...\n";
  progress_bar.Init();

  while(number_of_faces > m_TargetNumberOfFaces)
  {
    const bool merge_done = compute_next_mesh(m_Mesh, m_MeshProperties);

    number_of_faces = m_Mesh.number_of_faces();

    progress_bar.Update(static_cast<int>(a*number_of_faces+b), m_TargetNumberOfFaces);

    if(!merge_done)
      break;
  }

  progress_bar.Complete();

  // Extract simplified compound mesh.
  // Must be called before any garbage collection.
  compute_compound_mesh(m_Mesh, m_MeshProperties, m_CompoundMesh);

  // Identify faces at boundaries by type.
  identify_faces_at_boundaries(m_CompoundMesh, m_CompoundMeshProperties,
    min_c, max_c);

  // Collect garbage.
  m_CompoundMesh.collect_garbage();
  m_Mesh.collect_garbage();
}


//
// Included functions belonging to the class.
//

#include "marinov_init.hpp"
#include "marinov1.hpp"
#include "vertex_collapse.hpp"
#include "vertex_removal.hpp"
#include "boundaries.hpp"
#include "compound_mesh.hpp"

} // namespace marinov

#endif // MARINOV_FACE_DECIMATION_HPP