// Included directly in another file. No safeguard includes.template<typename KernelT>

template<typename KernelT>
bool Face_decimation<KernelT>::check_vertex_removal(const Mesh &mesh, const Vertex_index &vk, 
  const Mesh_properties &properties) const
{
  constexpr FT eps = std::numeric_limits<FT>::epsilon();
  bool removal_possible = false;
  bool injectivity = false;
  bool can_use_cgal_function = false;

  // Check if the cgal function remove_center_vertex can be used.
  {
    // "There are at least two distinct faces incident to the faces that are incident to target(h,g)"
    // https://doc.cgal.org/latest/BGL group__PkgBGLEulerOperations.html#ga8ec295610396e837258997c435b7a75f

    const Halfedge_index h_in  = mesh.halfedge(vk);
    const Halfedge_index h_opp = mesh.opposite(h_in);

    const Face_index faces[2] = { mesh.face(h_in), mesh.face(h_opp) };

    // Check if the two faces are on the same boundaries.
    if(properties.at(faces[0]).boundary != properties.at(faces[1]).boundary)
    {
      return removal_possible = false;
    }

    // Get the faces incident to the faces incident to vertex vk.
    std::set<Face_index> incident_faces[2];

    for(int i = 0; i < 2; ++i)
    {
      std::set<Face_index> inc_faces;
      CGAL::Face_around_face_iterator<Mesh> ff_it, ff_ite;
      boost::tie(ff_it, ff_ite) = CGAL::faces_around_face(mesh.halfedge(faces[i]), mesh);

      for(; ff_it != ff_ite; ++ff_it)
        inc_faces.insert(*ff_it);

      incident_faces[i] = inc_faces;
    }

    std::vector<Face_index> differences[2];
    for(int i = 0; i < 2; ++i)
      std::set_difference(incident_faces[i].begin(), incident_faces[i].end(),
                          incident_faces[(i+1) % 2].begin(), incident_faces[(i+1) % 2].end(),
                          std::back_inserter(differences[i]));

    if((differences[0].size() < 2) || (differences[1].size() < 2))
    {
      return removal_possible;
    }

    can_use_cgal_function = true;
  }

  // Check injectivity condition.
  {
    // Extract the two faces f0 and f1 in an independent partial mesh.
    Halfedge_index partial_h_in;
    std::map<Face_index, Face_index> partial_faces_correspondance;
    Mesh partial_mesh = extract_mesh_faces(mesh, mesh.halfedge(vk), partial_h_in,
      partial_faces_correspondance);

    if(partial_mesh.is_empty())
    {
      return removal_possible;
    }

    const Halfedge_index h_in = mesh.halfedge(vk);
    const Face_index new_face = mesh.face(h_in);
    const Face_index other_face = mesh.face(mesh.opposite(h_in));

    const FT &area0 = properties.at(new_face).area;
    const FT &area1 = properties.at(other_face).area;

    const Vector &normal0 = properties.at(new_face).normal;
    const Vector &normal1 = properties.at(other_face).normal;

    const Point &centroid0 = properties.at(new_face).centroid;
    const Point &centroid1 = properties.at(other_face).centroid;

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

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

    if(normal_norm <= eps)
      throw std::runtime_error("Face_decimation::update_properties_after_vertex_removal: normal norm <= epsilon. Should not happen whith merging already approved.");

    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
    );

    // Construct projection plane.
    const Plane plane(centroid, normal);
    std::vector<CGAL::Point_2<Kernel>> points_2d;

    // Get projected 2d points into the plane.
    {
      CGAL::Vertex_around_face_iterator<Mesh> vf_it, vf_ite;
      boost::tie(vf_it, vf_ite) = CGAL::vertices_around_face(partial_h_in, partial_mesh);

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

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

    if(!injectivity_condition(polygon))
    {
      return removal_possible;
    }

    injectivity = true;
  }

  removal_possible = (can_use_cgal_function && injectivity);

  return removal_possible;
}


template<typename KernelT>
void Face_decimation<KernelT>::compute_vertex_removal(Mesh &mesh, const Vertex_index &vk, 
  Face_index &other_face, Face_index &new_face)
{
  Halfedge_index h_in = mesh.halfedge(vk);
  new_face = mesh.face(h_in);
  other_face = mesh.face(mesh.opposite(h_in));

  // Compute vertex removal.
  h_in = CGAL::Euler::remove_center_vertex(h_in, mesh);

  // Get the correct faces.
  if(new_face != mesh.face(h_in))
  {
    std::swap(new_face, other_face);
  }
}


template<typename KernelT>
void Face_decimation<KernelT>::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)
{
  // Recompute the properties of the new face after vertex removal.
  constexpr FT eps = std::numeric_limits<FT>::epsilon();  

  const FT &area0 = properties[new_face].area;
  const FT &area1 = properties[other_face].area;

  const Vector &normal0 = properties[new_face].normal;
  const Vector &normal1 = properties[other_face].normal;

  const Point &centroid0 = properties[new_face].centroid;
  const Point &centroid1 = properties[other_face].centroid;

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

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

  if(normal_norm <= eps)
    throw std::runtime_error("Face_decimation::update_properties_after_vertex_removal: normal norm <= epsilon. Should not happen whith merging already approved.");

  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
  );

  // Update properties.
  {
    properties[new_face].centroid = centroid;
    properties[new_face].normal = normal;
    properties[new_face].area = area;
  }

  // 2D merged polygon.
  {
    Polygon polygon;
    CDT cdt;

    compute_2d_polygon(new_face, mesh, properties, polygon, cdt);

    if(!injectivity_condition(polygon))
    {
      throw std::runtime_error("Face_decimation::update_properties_after_vertex_removal: 2D polygon not simple. Should not happen whith merging already approved.");
    }
    else
    {
      properties[new_face].polygon = polygon;
      properties[new_face].triangulation = cdt;
    }
  }

  // Compute lifted triangulation and face area.
  {
    Mesh mesh_t;
    lift_triangulation(new_face, properties, mesh, mesh_t);

    // Update properties.
    properties[new_face].lifted_triangulation = mesh_t;
  }

  //
  // Update errors for all neighbouring faces of f0 and f1 and for f0 and f1.
  //
  {
    // new_face and its neighbours.
    std::set<Face_index> faces_to_update;
    faces_to_update.insert(new_face);
    
    CGAL::Face_around_face_iterator<Mesh> ff_it, ff_ite;
    boost::tie(ff_it, ff_ite) = CGAL::faces_around_face(mesh.halfedge(new_face), mesh);

    for(; ff_it != ff_ite; ++ff_it)
      faces_to_update.insert(*ff_it);

    // Update errors for faces inf faces_to_update.
    for(const Face_index &f_it : faces_to_update)
    {
      compute_merging_errors(mesh, f_it, properties);
    }
  }
  
  // Remove the properties of the removed face.
  properties.erase(other_face);
}