// Included directly in another file. No safeguard includes.

template<typename KernelT>
void Face_decimation<KernelT>::compute_initial_properties(const Mesh &mesh, 
  Mesh_properties &properties)
{
  typename Mesh::Face_range::iterator f_it = mesh.faces().begin(),
    f_ite = mesh.faces().end();

  for(; f_it != f_ite; ++f_it)
  {
    // Compute face area.
    properties[*f_it].area = CGAL::Polygon_mesh_processing::face_area(*f_it, mesh);

    // Compute face normal.
    properties[*f_it].normal = CGAL::Polygon_mesh_processing::compute_face_normal(*f_it, mesh);

    // Compute face centroid.
    compute_face_centroid(*f_it, mesh, properties);

    // Compute the 2D polygon projected into the plane defined by the centroid
    // and the normal.
    Polygon polygon;
    CDT cdt;

    compute_2d_polygon(*f_it, mesh, properties, polygon, cdt);
    properties[*f_it].polygon = polygon;
    properties[*f_it].triangulation = cdt;

    // Compute the lifted triangulation.
    Mesh mesh_t;
    lift_triangulation(*f_it, properties, mesh, mesh_t);
    properties[*f_it].lifted_triangulation = mesh_t;
  }

  // Compute the minimal merging errors for each face.
  compute_merging_errors(mesh, properties);
}


template<typename KernelT>
void Face_decimation<KernelT>::compute_face_centroid(const Face_index &f, const Mesh &mesh, 
  Mesh_properties &properties)
{
  Vector centroid(0.0, 0.0, 0.0);
  std::size_t valence = 0;

  CGAL::Vertex_around_face_iterator<Mesh> vf_it, vf_ite;
  boost::tie(vf_it, vf_ite) = CGAL::vertices_around_face(mesh.halfedge(f), mesh);
      
  for(; vf_it != vf_ite; ++vf_it, ++valence)
  {
    const Point &p = mesh.point(*vf_it);
    centroid += Vector(p[0], p[1], p[2]);
  }

  centroid /= valence;
  properties[f].centroid = Point(centroid[0], centroid[1], centroid[2]);
}


template<typename KernelT>
void Face_decimation<KernelT>::compute_2d_polygon(const Face_index &f, const Mesh &mesh, 
  const Mesh_properties &properties, Polygon &polygon, CDT &cdt)
{
  // Construct the projection plane.
  const Plane plane(properties.at(f).centroid, properties.at(f).normal);
  std::vector<CGAL::Point_2<Kernel>> points_2d;
  std::vector<std::size_t> index;

  // Get projected 2d point into the plane.
  CGAL::Vertex_around_face_iterator<Mesh> vf_it, vf_ite;
  boost::tie(vf_it, vf_ite) = CGAL::vertices_around_face(mesh.halfedge(f), mesh);

  for(; vf_it != vf_ite; ++vf_it)
  {
    points_2d.push_back(plane.to_2d(mesh.point(*vf_it)));
    index.push_back(*vf_it);
  }

  // Construct and return the polygon.
  polygon = Polygon(points_2d.begin(), points_2d.end());

  cdt.insert_constraint(polygon.vertices_begin(),
                        polygon.vertices_end());

  // Associate vertices of the 2D triangulation to the vertices of the mesh.
  std::size_t i = 0;
  for(typename CDT::Finite_vertices_iterator cdt_vit = cdt.finite_vertices_begin();
      cdt_vit != cdt.finite_vertices_end(); ++cdt_vit, ++i)
  {
    cdt_vit->info() = index[i];
  }
}


template<typename KernelT>
void Face_decimation<KernelT>::lift_triangulation(const Face_index &f,
  const Mesh_properties &properties, const Mesh &mesh, Mesh &mesh_t) const
{
  mesh_t.clear();
  const CDT &cdt = properties.at(f).triangulation;
  std::map<std::size_t, Vertex_index> vhandle;

  // Vertices.
  for(typename CDT::Finite_vertices_iterator cdt_vit = cdt.finite_vertices_begin();
      cdt_vit != cdt.finite_vertices_end(); ++cdt_vit)
  {
    const Point &p = mesh.point(Vertex_index(cdt_vit->info()));
    vhandle[cdt_vit->info()] = mesh_t.add_vertex(p);
  }

  // Faces.
  for(typename CDT::Finite_faces_iterator cdt_fit = cdt.finite_faces_begin();
      cdt_fit != cdt.finite_faces_end(); ++cdt_fit)
  {
    // Check if current Delaunay triangle is inside the polygon.
    const typename Kernel::Triangle_2 tri = cdt.triangle(cdt_fit);

    // Check if the centroid of the current triangle is inside the polygon.
    if(properties.at(f).polygon.has_on_bounded_side(centroid(tri)))
    {
      std::vector<Vertex_index> face_vhandles(3);

      for(int i = 0; i < 3; ++i)
        face_vhandles[i] = vhandle[cdt_fit->vertex(i)->info()];

      mesh_t.add_face(face_vhandles);
    }
  }
}


template<typename KernelT>
void Face_decimation<KernelT>::compute_merging_errors(const Mesh &mesh, Mesh_properties &properties)
{
  typename Mesh::Face_range::iterator f_it = mesh.faces().begin(),
    f_ite = mesh.faces().end();

  for(; f_it != f_ite; ++f_it)
  {
    compute_merging_errors(mesh, *f_it, properties);
  }
}


template<typename KernelT>
void Face_decimation<KernelT>::compute_merging_errors(const Mesh &mesh, const Face_index &f,
  Mesh_properties &properties)
{
    // Set initial error to the max possible value.
    properties[f].error = std::numeric_limits<FT>::max();
    properties[f].mergeable = false;
    properties[f].face_to_merge = Mesh::null_face();

    // Loop around the neighbouring faces of the current face.
    CGAL::Face_around_face_iterator<Mesh> ff_it, ff_ite;
    boost::tie(ff_it, ff_ite) = CGAL::faces_around_face(mesh.halfedge(f), mesh);

    for(; ff_it != ff_ite; ++ff_it)
    {
      // Compute the merging error between f_it and ff_it.
      bool mergeable;
      const FT error = compute_merging_error(mesh, properties, f, *ff_it, mergeable);

      // The resulting merged face does not respect the injectivity condition.
      // Reject merge.
      if(!mergeable) continue;

      if(properties[f].error > error)
      {
        properties[f].error = error;
        properties[f].mergeable = true;
        properties[f].face_to_merge = *ff_it;
      }
    }
}


template<typename KernelT>
typename Face_decimation<KernelT>::FT 
Face_decimation<KernelT>::compute_merging_error(const Mesh &mesh, const Mesh_properties &properties,
  const Face_index &f0, const Face_index &f1, bool &mergeable) const
{
  if(mesh.is_removed(f0) || mesh.is_removed(f1))
    throw std::runtime_error("Face_decimation::compute_merging_error: can not merge at least one face has been removed.");

  constexpr FT eps = std::numeric_limits<FT>::epsilon();
  FT error = std::numeric_limits<FT>::max();

  const FT &area0 = properties.at(f0).area;
  const FT &area1 = properties.at(f1).area;

  const Vector &normal0 = properties.at(f0).normal;
  const Vector &normal1 = properties.at(f1).normal;

  const Point &centroid0 = properties.at(f0).centroid;
  const Point &centroid1 = properties.at(f1).centroid;

  const Mesh &mesh_t0 = properties.at(f0).lifted_triangulation;
  const Mesh &mesh_t1 = properties.at(f1).lifted_triangulation;

  // Do not merge faces belonging to different boundaries.
  if(properties.at(f0).boundary != properties.at(f1).boundary)
  {
    mergeable = false;
    return error;
  }

  // Mergeable.
  mergeable = true;

  // Face area.
  const FT area = area0 + area1;

  // Face normal.
  Vector normal = (area0*normal0 + area1*normal1);
  const FT normal_norm = sqrt(normal.squared_length());

  // If the normal norm is small, this indicates a big L2,1 error.
  // The face will not be selected for merging. So return here.
  if(normal_norm <= eps)
  {
    mergeable = false;
    return error;
  }

  normal /= normal_norm;

  // Face centroid.
  Point centroid(
    (area0*centroid0[0] + area1*centroid1[0]) / area,
    (area0*centroid0[1] + area1*centroid1[1]) / area,
    (area0*centroid0[2] + area1*centroid1[2]) / area
  );
  
  // Merge faces f0 and f1.
  // (Does not actually merge the faces in the mesh).
  std::vector<FT> inner_angles;
  std::vector<Vertex_index> merged_face;
  merged_face = merge_face(mesh, f0, f1, normal, inner_angles);

  if(merged_face.empty()) // Face merge not possible (typically hole).
  {
    mergeable = false;
    return error;
  }

  // Projection plane.
  const Plane plane(centroid, normal);

  // 2D merged polygon.
  {
    std::vector<typename Kernel::Point_2> points_2d;
    points_2d.reserve(merged_face.size());

    for(const Vertex_index &v : merged_face)
    {
      points_2d.push_back(plane.to_2d(mesh.point(v)));
    }

    const Polygon polygon(points_2d.begin(), points_2d.end());

    // Test for injectivity condition.
    if(!injectivity_condition(polygon))
    {
      mergeable = false;
      return error;
    }
  }

  // Compute the L2,1 error.
  const FT l21_error = area0 * (normal0 - normal).squared_length()
                     + area1 * (normal1 - normal).squared_length();
  
  // Compute the L2 error with respect to the triangulations of f0 and f1.
  FT l2_error;
  {
    const FT l2_error0 = compute_l2_error(mesh_t0, plane);
    const FT l2_error1 = compute_l2_error(mesh_t1, plane);
    
    l2_error = l2_error0 + l2_error1;
  }

  // Compute the total error.
  {
    error = (1.0 + l2_error) * (1.0 + l21_error);

    // Modification of the total error computation in order to favour
    // decimation of low valence faces  and penalize non-regular angles (Marinov et al.)
    const std::size_t m = merged_face.size();
    error *= (1 + 0.01*((m-4)*(m-4)));

    FT inner_angles_error(0);
    for(std::size_t i = 0; i < inner_angles.size(); ++i)
      inner_angles_error += (inner_angles[i] - M_PI/2.0) * (inner_angles[i] - M_PI/2.0);

    error *= (1 + 0.01*inner_angles_error);
  }

  return error;
}


template<typename KernelT>
typename Face_decimation<KernelT>::FT 
Face_decimation<KernelT>::compute_l2_error(const Mesh &mesh, const Plane &plane) const
{
  FT error = 0.0;

  typename Mesh::Face_range::iterator f_it = mesh.faces().begin(),
    f_ite = mesh.faces().end();
  for(; f_it != f_ite; ++f_it)
  {
    const FT face_area = CGAL::Polygon_mesh_processing::face_area(*f_it, mesh); // OK.
    FT d[3]; // Distances of the vertices of the current triangle from the projection plane.

    CGAL::Vertex_around_face_iterator<Mesh> vf_it, vf_ite;
    boost::tie(vf_it, vf_ite) = CGAL::vertices_around_face(mesh.halfedge(*f_it), mesh);

    for(int i = 0; vf_it != vf_ite; ++vf_it, ++i)
    {
      // Project mesh point onto the plane.
      const Point &p = mesh.point(*vf_it);
      const Point p_proj = plane.projection(p);
      d[i] = sqrt((p - p_proj).squared_length());
    }

    error += face_area/6.0 * (d[0]*d[0] + d[1]*d[1] + d[2]*d[2] +
                              d[0]*d[1] + d[0]*d[2] + d[1]*d[2]);
  }

  return error;
}

// Note: this merging algorithm ignores faces with holes.
// Faces must be touching by at least one edge and non-overlap.
// Faces may be non-convex.
// Does not actually merge the faces in the mesh.
template<typename KernelT>
std::vector<typename Face_decimation<KernelT>::Vertex_index>
Face_decimation<KernelT>::merge_face(const Mesh &mesh,
  const Face_index &f0, const Face_index &f1,
  const Vector &normal, std::vector<FT> &alpha) const
{
  if(mesh.is_removed(f0) || mesh.is_removed(f1))
    throw std::runtime_error("Face_decimation::merge_face: can not merge at least one face has been removed.");

  alpha.clear();
  std::vector<Vertex_index> merged_face;
  bool have_hole = false;

  // Determine common edges between the faces f0 and f1.
  std::vector<Edge_index> edges0, edges1;
  std::set<Edge_index> intersection_edges, difference_edges; // Inner and boundary edges.
  {
    CGAL::Edge_around_face_iterator<Mesh> ef_it, ef_ite;
    boost::tie(ef_it, ef_ite) = CGAL::edges_around_face(mesh.halfedge(f0), mesh);
    for(; ef_it != ef_ite; ++ef_it)
      edges0.push_back(*ef_it);
  }

  {
    CGAL::Edge_around_face_iterator<Mesh> ef_it, ef_ite;
    boost::tie(ef_it, ef_ite) = CGAL::edges_around_face(mesh.halfedge(f1), mesh);
    for(; ef_it != ef_ite; ++ef_it)
      edges1.push_back(*ef_it);
  }

  std::sort(edges0.begin(), edges0.end());
  std::sort(edges1.begin(), edges1.end());
  std::set_intersection(edges0.begin(), edges0.end(),
                        edges1.begin(), edges1.end(),
                        std::inserter(intersection_edges, intersection_edges.end()));
  std::set_symmetric_difference(edges0.begin(), edges0.end(),
                                edges1.begin(), edges1.end(),
                                std::inserter(difference_edges, difference_edges.end()));
  
  // Loop over the difference edges and check if there only one loop
  // or if there is a hole.
  {
    const Vertex_index start_vertex = mesh.vertex(*difference_edges.begin(), 0);
    Vertex_index current_vertex = start_vertex;
    std::set<Edge_index> visited_edges;

    do
    {
      // Look for a vertex adjacent to the current vertex
      // whose shared edge is not a common edge between f0 and f1
      // and belongs to f0 or f1.
      CGAL::Halfedge_around_source_iterator<Mesh> hs_it, hs_ite;
      boost::tie(hs_it, hs_ite) = CGAL::halfedges_around_source(current_vertex, mesh);

      for(; hs_it != hs_ite; ++hs_it)
      {
        const Edge_index e = mesh.edge(*hs_it);

        if( (difference_edges.find(e) != difference_edges.end()) &&
            (visited_edges.find(e) == visited_edges.end()) ) // Avoid to follow an already followed edge.
        {
          current_vertex = mesh.target(*hs_it);
          merged_face.push_back(current_vertex);
          visited_edges.insert(e);
          break;
        }
      }
    } while(current_vertex != start_vertex);

    // We do not traverse all the boundary edges, meaning that there is
    // a least one hole.
    have_hole = (visited_edges.size() != difference_edges.size());
  }

  // Face merging generates a hole. Stop here.
  if(have_hole)
  {
    merged_face.clear();
    return merged_face;
  }

  // Compute the inner angles at the vertices incident to the removed edges.
  {
    // Look for vertices coincident to a boundary edge and an inner edge.
    std::set<Vertex_index> inner_vertices;
    const Face_index faces[2] = {f0, f1};

    for(int i = 0; i < 2; ++i)
    {
      CGAL::Vertex_around_face_iterator<Mesh> vf_it, vf_ite;
      boost::tie(vf_it, vf_ite) = CGAL::vertices_around_face(mesh.halfedge(faces[i]), mesh);

      for(; vf_it != vf_ite; ++vf_it)
      {
        // Check if current vertex is incident to a boundary and an inner edge.
        bool incident_inner_edge_found = false;
        bool incident_boundary_edge_found = false;

        CGAL::Halfedge_around_source_iterator<Mesh> hs_it, hs_ite;
        boost::tie(hs_it, hs_ite) = CGAL::halfedges_around_source(*vf_it, mesh);

        for(; hs_it != hs_ite; ++hs_it)
        {
          const Edge_index e = mesh.edge(*hs_it);

          if(intersection_edges.find(e) != intersection_edges.end())
          {
            incident_inner_edge_found = true;
          }
          else if(difference_edges.find(e) != difference_edges.end())
          {
            incident_boundary_edge_found = true;
          }
        } // Loop on outgoing halfedges of the current vertex.

        if(incident_inner_edge_found && incident_boundary_edge_found)
           inner_vertices.insert(*vf_it);
      } // Loop on vertices.
    } // For i.

    // For each inner vertex, get the two outgoing halfedge belonging to f0 or f1.
    // Then measure the angle between these halfedges.

    // Look the outgoing halfedges around the inner vertices.
    for(const Vertex_index &v : inner_vertices)
    {
      Halfedge_index boundary_halfedges[2];
      int i = 0;

      CGAL::Halfedge_around_source_iterator<Mesh> hs_it, hs_ite;
      boost::tie(hs_it, hs_ite) = CGAL::halfedges_around_source(v, mesh);

      for(; hs_it != hs_ite; ++hs_it)
      {
        const Edge_index e = mesh.edge(*hs_it);

        if( (difference_edges.find(e) != difference_edges.end()) &&
            (intersection_edges.find(e) == intersection_edges.end()) )
        {
          boundary_halfedges[i] = *hs_it;
          ++i;
        }
      }

      if(i != 2)
        throw std::runtime_error("Face_decimation::merge_face: inner vertex adjacent to less than two boundary edges, unable to compute error on inner angles.");

      Vector vec0(mesh.point(mesh.source(boundary_halfedges[0])),
                  mesh.point(mesh.target(boundary_halfedges[0])));
      Vector vec1(mesh.point(mesh.source(boundary_halfedges[1])),
                  mesh.point(mesh.target(boundary_halfedges[1])));

      // Orient vectors.
      if((mesh.face(boundary_halfedges[0])) != f0 && (mesh.face(boundary_halfedges[0]) != f1))
        std::swap(vec0, vec1);

      // Normalize.
      vec0 /= sqrt(vec0.squared_length());
      vec1 /= sqrt(vec1.squared_length());

      // Compute angle.
      const FT x = (vec0 * vec1);
      const FT y = cross_product(vec0, vec1) * normal;

      alpha.push_back(atan2(y, x));
    }
  }

  if(merged_face.size() < 3)
    throw std::runtime_error("Face_decimation::merge_face: vector < 3.");

  return merged_face;
}


template<typename KernelT>
bool Face_decimation<KernelT>::injectivity_condition(const Polygon &polygon) const
{
  const bool polygon_simple = polygon.is_simple();
  bool polygon_collinear = true;

  if(polygon_simple)
  {
    const std::size_t s = polygon.size();

    for(std::size_t i = 0; i < s; ++i)
    {
      bool have_collinearity = false;
      const typename Kernel::Point_2 &p0 = polygon[i % s];
      const typename Kernel::Point_2 &p1 = polygon[(i+1) % s];
      const typename Kernel::Point_2 &p2 = polygon[(i+2) % s];

      if(collinear(p0, p1, p2))
      {
        have_collinearity = true;
      }

      polygon_collinear &= have_collinearity;
    }
  }

  const bool ok = polygon_simple && (!polygon_collinear);
  return ok;
}