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

template<typename KernelT>
bool Face_decimation<KernelT>::compute_next_mesh(Mesh &mesh, Mesh_properties &properties)
{
  bool merge_done = false;
  Face_index face_to_merge, new_face;

  // Merge faces in mesh with the least error.
  merge_done = compute_next_merged_face(mesh, properties, face_to_merge, new_face);

  // Update properties for the merged face and its neighbouring faces.
  if(merge_done)
    update_properties_after_face_merge(mesh, properties, face_to_merge, new_face);

  seek_and_destroy_vertices_of_valence2(mesh, properties);

  return merge_done;
}


// Note: do not triangulate plane faces.
template<typename KernelT>
void Face_decimation<KernelT>::compute_compound_mesh(const Mesh &mesh, 
  const Mesh_properties &properties, Mesh &compound_mesh)
{
  constexpr FT angle_threshold = 0.0; //1e-1; // Degree.
  constexpr FT quality_threshold = 1.0; //0.99;
  compound_mesh.clear();

  //
  // Construct compound mesh.
  //

  // Vertices.
  std::map<std::size_t, Vertex_index> vhandle;

  for(const std::pair<Face_index, Face_properties> &prop : properties)
  {
    const CDT &cdt = prop.second.triangulation;
    for(typename CDT::Finite_vertices_iterator cdt_vit = cdt.finite_vertices_begin();
        cdt_vit != cdt.finite_vertices_end(); ++cdt_vit)
    {
      const Vertex_index v_index(cdt_vit->info());

      // Dot not insert already inserted vertices.
      if(vhandle.find(v_index) == vhandle.end())
      {
        const Point &p = mesh.point(v_index);
        vhandle[cdt_vit->info()] = compound_mesh.add_vertex(p);
      }
    }
  }

  // Faces.
  std::size_t id = 0;
  for(const std::pair<Face_index, Face_properties> &prop : properties)
  {
    // Do not triangulate any face.
    if(!m_TriangulateMesh)
    {
      CGAL::Vertex_around_face_iterator<Mesh> vf_it, vf_ite;
      boost::tie(vf_it, vf_ite) = CGAL::vertices_around_face(mesh.halfedge(prop.first), mesh);
      std::vector<Vertex_index> face_vhandles;

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

      const Face_index face = compound_mesh.add_face(face_vhandles);
      continue;
    }

    // Test whether the face is plane or not.
    bool is_face_plane = false;
    
    // Faces with 3 points are automatically plane.
    is_face_plane = (mesh.degree(prop.first) == 3);

    // Face is plane. Do not triangulate and mark it as plane.
    if(is_face_plane)
    {
      CGAL::Vertex_around_face_iterator<Mesh> vf_it, vf_ite;
      boost::tie(vf_it, vf_ite) = CGAL::vertices_around_face(mesh.halfedge(prop.first), mesh);
      std::vector<Vertex_index> face_vhandles;

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

      const Face_index face = compound_mesh.add_face(face_vhandles);
      continue;
    }

    // Compare the face normals of the lifted triangulation.
    {
      const CDT &cdt = prop.second.triangulation;
      std::vector<Vector> normals;

      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(prop.second.polygon.has_on_bounded_side(centroid(tri)))
        {
          std::vector<Point> points(3);

          for(int i = 0; i < 3; ++i)
            points[i] = mesh.point(Vertex_index(cdt_fit->vertex(i)->info()));

          const Plane plane(points[0], points[1], points[2]);
          Vector normal = plane.orthogonal_vector();
          normal /= sqrt(normal.squared_length());
          normals.push_back(normal);
        }
      }

      // Compute the maximum angle between normals.
      FT max_angle = 0.0;
      for(std::size_t i = 0; i < normals.size(); ++i)
      {
        for(std::size_t j = i+1; j < normals.size(); ++j)
        {
          const FT angle = acos(normals[i]*normals[j]);

          if(max_angle < angle)
            max_angle = angle;
        }
      }
      
      is_face_plane = (max_angle*180.0/M_PI <= angle_threshold);
    }

    // Face is (almost) plane. Do not triangulate and mark it as plane.
    if(is_face_plane)
    {
      CGAL::Vertex_around_face_iterator<Mesh> vf_it, vf_ite;
      boost::tie(vf_it, vf_ite) = CGAL::vertices_around_face(mesh.halfedge(prop.first), mesh);
      std::vector<Vertex_index> face_vhandles;

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

      const Face_index face = compound_mesh.add_face(face_vhandles);
      continue;
    }

    // Try to fit a plane to the vertices of the face.
    {
      std::vector<Point> points;
      CGAL::Vertex_around_face_iterator<Mesh> vf_it, vf_ite;
      boost::tie(vf_it, vf_ite) = CGAL::vertices_around_face(mesh.halfedge(prop.first), mesh);

      for(; vf_it != vf_ite; ++vf_it)
      {
        const Point &p = mesh.point(*vf_it);
        points.push_back(p);
      }

      Plane plane;
      const FT quality = linear_least_squares_fitting_3(
        points.begin(), points.end(), plane, CGAL::Dimension_tag<0>());

      is_face_plane = (quality >= quality_threshold);
    }

    // Face is (almost) plane. Do not triangulate and mark it as plane.
    if(is_face_plane)
    {
      CGAL::Vertex_around_face_iterator<Mesh> vf_it, vf_ite;
      boost::tie(vf_it, vf_ite) = CGAL::vertices_around_face(mesh.halfedge(prop.first), mesh);
      std::vector<Vertex_index> face_vhandles;

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

      const Face_index face = compound_mesh.add_face(face_vhandles);
      continue;
    }

    // Face not plane. Construct triangulated compound face.
    {
      const CDT &cdt = prop.second.triangulation;

      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(prop.second.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()];

          const Face_index face = compound_mesh.add_face(face_vhandles);
        }
      }
    }

    ++id;
  }
}