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

template<typename KernelT>
void Face_decimation<KernelT>::seek_and_destroy_vertices_of_valence2(Mesh &mesh, 
  Mesh_properties &properties)
{
  // Look for valence 2 vertices.
  typename Mesh::Vertex_range::iterator v_it = mesh.vertices().begin(),
    v_ite = mesh.vertices().end();

  for(; v_it != v_ite; ++v_it)
  {
    if(mesh.degree(*v_it) == 2)
    {
      // Check if the vertex collapse is possible.
      const bool current_vertex_collapse_possible = check_vertex_collapse(mesh,
        *v_it, properties);

      if(current_vertex_collapse_possible)
      {
        Face_index left_face, right_face;
        compute_vertex_collapse(mesh, *v_it, left_face, right_face);

        update_properties_after_vertex_collapse(mesh, *v_it, left_face, right_face,
          properties);
      }
      else // Check if vertex removal is possible.
      {
        const bool current_vertex_removal_possible = check_vertex_removal(mesh,
          *v_it, properties);

        if(current_vertex_removal_possible)
        {
          Face_index other_face, new_face;
          compute_vertex_removal(mesh, *v_it, other_face, new_face);

          update_properties_after_vertex_removal(mesh, *v_it, other_face, new_face, 
            properties);
        }
      }
    }
  }
}


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

  // Adjacent faces to vk.
  Face_index f0, f1;
  const Halfedge_index h_in = mesh.halfedge(vk);

  // Check if new edge after collapse is not an already existing edge.
  {
    const Halfedge_index h_next = mesh.next(h_in);
    
    // The two vertices neighbouring vk (which is of valence 2).
    const Vertex_index vj = mesh.source(h_in);
    const Vertex_index vl = mesh.target(h_next);

    // Get the faces.
    f0 = mesh.face(h_in);
    f1 = mesh.face(mesh.opposite(h_in));

    // Look for the vertices neighbouring vl and check
    // if vj is a neighbour.
    CGAL::Vertex_around_target_iterator<Mesh> vv_it, vv_ite;
    boost::tie(vv_it, vv_ite) = CGAL::vertices_around_target(h_next, mesh);

    for(; vv_it != vv_ite; ++vv_it)
    {
      if(*vv_it == vj)
      {
        return collapse_possible;
      }
    }

    existing_edge = false;
  }

  // Check the injectivity condition.
  {
    // Check if vertex vk can be collapsed with function 'join_vertex'.
    if( (mesh.degree(f0) < 4) || (mesh.degree(f1) < 4) )
    {
//std::cout << "check vertex collapse: can not use join_vertex\n";
      return collapse_possible;
    }

    // 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, inverse_faces_correspondance;
    Mesh partial_mesh = extract_mesh_faces(mesh, h_in, partial_h_in, partial_faces_correspondance);

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

    // Set the face correspondance following on which face h_in is incident to.
    inverse_faces_correspondance[partial_faces_correspondance[f0]] = f0;
    inverse_faces_correspondance[partial_faces_correspondance[f1]] = f1;

    // Collapse vertex vk.
    partial_h_in = CGAL::Euler::join_vertex(partial_mesh.opposite(partial_h_in), partial_mesh);
      // Returns: prev(opposite(h,partial_mesh))

    // Compute the two associated polygons to the two faces and check for injectivity.
    {
      typename Mesh::Face_range::iterator f_it = partial_mesh.faces().begin(),
        f_ite = partial_mesh.faces().end();

      for(; f_it != f_ite; ++f_it)
      {
        // Construct the projection plane.
        const Plane plane(properties.at(inverse_faces_correspondance[*f_it]).centroid, 
          properties.at(inverse_faces_correspondance[*f_it]).normal);
        std::vector<CGAL::Point_2<Kernel>> points_2d;

        // 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(partial_mesh.halfedge(*f_it),
          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 collapse_possible;
        }
      }

      injectivity = true;
    }

    collapse_possible = (injectivity && !existing_edge);
  }

  return collapse_possible;
}


template<typename KernelT>
typename Face_decimation<KernelT>::Mesh 
Face_decimation<KernelT>::extract_mesh_faces(const Mesh &mesh, 
  const Halfedge_index &h, Halfedge_index &new_h,
  std::map<Face_index, Face_index> &faces_correspondance) const
{
  Mesh out_mesh;
  std::map<Vertex_index, Vertex_index> vhandle;
  const Face_index &f0 = mesh.face(h);
  const Face_index &f1 = mesh.face(mesh.opposite(h));
  const Vertex_index &vk = mesh.target(h); // Vertex to collapse.

  // Add vertices of face 0.
  {
    CGAL::Vertex_around_face_iterator<Mesh> vf_it, vf_ite;
    boost::tie(vf_it, vf_ite) = CGAL::vertices_around_face(h, mesh);

    for(; vf_it != vf_ite; ++vf_it)
    {
      const Point &p = mesh.point(*vf_it);
      vhandle[*vf_it] = out_mesh.add_vertex(p);
    }
  }

  // Add vertices of face 1.
  {
    CGAL::Vertex_around_face_iterator<Mesh> vf_it, vf_ite;
    boost::tie(vf_it, vf_ite) = CGAL::vertices_around_face(mesh.opposite(h), mesh);

    for(; vf_it != vf_ite; ++vf_it)
    {
      if(vhandle.find(*vf_it) == vhandle.end()) // Avoid duplicating vertices.
      {
        const Point &p = mesh.point(*vf_it);
        vhandle[*vf_it] = out_mesh.add_vertex(p);
      }
    }
  }

  // Add face 0.
  Face_index out_f0, out_f1;

  {
    std::vector<Vertex_index> face_vhandles;
    CGAL::Vertex_around_face_iterator<Mesh> vf_it, vf_ite;
    boost::tie(vf_it, vf_ite) = CGAL::vertices_around_face(h, mesh);

    for(; vf_it != vf_ite; ++vf_it)
      face_vhandles.push_back(vhandle[*vf_it]);

    out_f0 = out_mesh.add_face(face_vhandles);
  }

  // Add face 1.
  {
    std::vector<Vertex_index> face_vhandles;
    CGAL::Vertex_around_face_iterator<Mesh> vf_it, vf_ite;
    boost::tie(vf_it, vf_ite) = CGAL::vertices_around_face(mesh.opposite(h), mesh);

    for(; vf_it != vf_ite; ++vf_it)
      face_vhandles.push_back(vhandle[*vf_it]);

    out_f1 = out_mesh.add_face(face_vhandles);
  }

  // Compute new_h.
  {
    new_h = Mesh::null_halfedge();

    // Loop over halfedges of face out_f0.
    CGAL::Halfedge_around_face_iterator<Mesh> hf_it, hf_ite;
    boost::tie(hf_it, hf_ite) = CGAL::halfedges_around_face(out_mesh.halfedge(out_f0), out_mesh);

    for(; hf_it != hf_ite; ++hf_it)
    {
      if(out_mesh.target(*hf_it) == vhandle[vk])
      {
        new_h = *hf_it;
        break;
      }
    }
  }

  if(new_h == Mesh::null_halfedge())
    throw std::runtime_error("Face_decimation::extract_mesh_faces: null halfedge.");

  // Set the correspondance between faces of the mesh and faces of the out mesh.
  faces_correspondance[f0] = out_f0;
  faces_correspondance[f1] = out_f1;

  // Strange topology. Do not deal with it. (Typically: merge of two faces leave one or several holes inside it.)
  if(out_mesh.degree(out_mesh.target(new_h)) != 2)
  {
    out_mesh.clear();
    faces_correspondance.clear();
  }

  return out_mesh;
}

template<typename KernelT>
void Face_decimation<KernelT>::compute_vertex_collapse(Mesh &mesh, const Vertex_index &vk,
  Face_index &left_face, Face_index &right_face)
{
  // Adjacent faces to vk.
  const Halfedge_index h_in = mesh.halfedge(vk);
  const Halfedge_index h_in_opp = mesh.opposite(h_in);

  right_face = mesh.face(h_in);
  left_face  = mesh.face(h_in_opp);

  // Only mark vertices and other elements as deleted.
  CGAL::Euler::join_vertex(h_in_opp, mesh);
}


template<typename KernelT>
void Face_decimation<KernelT>::update_properties_after_vertex_collapse(const Mesh &mesh, 
  const Vertex_index &vk, const Face_index &f0, const Face_index &f1, Mesh_properties &properties)
{
  constexpr FT eps = std::numeric_limits<FT>::epsilon();
  const Face_index f[2] = {f0, f1};

  //
  // Update properties of the two faces adjacent to vertex vk.
  //
  {
    for(int i = 0; i < 2; ++i)
    {
      // Centroid and normal are left unchanged after vertex collapse.

      // Construct the projection plane.
      const Plane plane(properties[f[i]].centroid, properties[f[i]].normal);
      std::vector<CGAL::Point_2<Kernel>> points_2d;
      std::vector<std::size_t> index;

      // 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(mesh.halfedge(f[i]), mesh);

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

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

        if(!injectivity_condition(polygon))
        {
          std::string error_message = "Face_decimation::update_properties_after_vertex_collapse: polygon associated to face " + std::to_string(f[i]) + " not simple where it should be.";
          throw std::runtime_error(error_message.c_str());
        }

        // Compute constrained Delaunay triangulation.
        CDT cdt;
        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];
          }
        }

        // Update properties.
        properties[f[i]].polygon = polygon;
        properties[f[i]].triangulation = cdt;
      }

      // Compute lifted triangulation and face area.
      {
        Mesh mesh_t;
        lift_triangulation(f[i], properties, mesh, mesh_t);

        const FT area = CGAL::Polygon_mesh_processing::area(mesh_t);

        // Update properties.
        properties[f[i]].area = area;
        properties[f[i]].lifted_triangulation = mesh_t;
      }
    } // For i = 0,1.
  }

  //
  // Update errors for all neighbouring faces of f0 and f1 and for f0 and f1.
  //
  {
    // Faces f0, f1 and their neighbours once.
    std::set<Face_index> faces_to_update;

    for(int i = 0; i < 2; ++i)
    {
      faces_to_update.insert(f[i]);

      CGAL::Face_around_face_iterator<Mesh> ff_it, ff_ite;
      boost::tie(ff_it, ff_ite) = CGAL::faces_around_face(mesh.halfedge(f[i]), 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);
    }
  }
}