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

template<typename KernelT>
bool Face_decimation<KernelT>::compute_next_merged_face(Mesh &mesh, 
  const Mesh_properties &properties,
  Face_index &face_to_merge, Face_index &new_face)
{
  std::set<Face_index> excluded_faces;

SEARCH_FOR_FACE_TO_MERGE:

  // Look for the first legal face with the smallest error.
  bool face_found = false;
  bool merge_done = false;
  face_to_merge = new_face = Mesh::null_face();
  FT error = std::numeric_limits<FT>::max();

  for(const auto &it : properties)
  {
    // Face can not be merged.
    if(!it.second.mergeable) continue;

    if( (error > it.second.error) &&
        (excluded_faces.find(it.first) == excluded_faces.end()) )
    {
      error = it.second.error;
      face_to_merge = it.first;
      face_found = true;
    }
  }

  if(!face_found)
    return merge_done;

  // Look for an edge shared between the two faces to merge and merge them, if possible.
  Halfedge_index new_h;
  {
    const Face_index &other_face = properties.at(face_to_merge).face_to_merge;
    CGAL::Halfedge_around_face_iterator<Mesh> hf_it, hf_ite;
    boost::tie(hf_it, hf_ite) = CGAL::halfedges_around_face(mesh.halfedge(face_to_merge), mesh);

    for(; hf_it != hf_ite; ++hf_it)
    {
      const Halfedge_index h_opp = mesh.opposite(*hf_it);

      // Shared edge found. Join the faces, if possible.
      if( (mesh.face(h_opp) == other_face) &&
          (mesh.degree(mesh.source(*hf_it)) >= 3) &&
          (mesh.degree(mesh.target(*hf_it)) >= 3) )
      {
        new_h = CGAL::Euler::join_face(*hf_it, mesh);
        merge_done = true;
        break;
      }
    }
  }

  // Merging failed. Try to do merging with the next smallest error.
  if(!merge_done)
  {
    excluded_faces.insert(face_to_merge);
    goto SEARCH_FOR_FACE_TO_MERGE;
  }

  // Note 1: do not collect garbage here, as the face index 'new_face' will
  // probabbly be pointing to a wrong element and cause undefined behavior.

  // Note 2: iterators iterate over non deleted elements.

  // Return the new face index.
  if(merge_done) 
  {
    new_face = mesh.face(new_h);
    (*m_Os) << "Current number of faces: " << mesh.number_of_faces() 
            << " (target: " << m_TargetNumberOfFaces << ", error: " << (error-1) << ")\n";
  }

  return merge_done;
}

template<typename KernelT>
void Face_decimation<KernelT>::update_properties_after_face_merge(const Mesh &mesh,
  Mesh_properties &properties, const Face_index &face_to_merge, const Face_index &new_face)
{
  if(face_to_merge != new_face)
    throw std::runtime_error("Face_decimation::update_properties_after_face_merge: merged face mismatch.");

  // Update the properties of the current merged face.
  const Face_index other_face = properties[new_face].face_to_merge; // Do not put a reference here...
                                                                    // Because this properties is updated and we will need the non-updated value after the update.

  {
    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_face_merge: 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].area = area;
      properties[new_face].normal = normal;
      properties[new_face].centroid = centroid;
    }

    // 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_face_merge: 2D polygon not simple. Should not happen whith merging already approved.");
      }
      else
      {
        properties[new_face].polygon = polygon;
        properties[new_face].triangulation = cdt;
      }
    }

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

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

  // Update the errors of the current merged face and its neighbouring faces.
  {
    compute_merging_errors(mesh, new_face, properties);

    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)
    {
      compute_merging_errors(mesh, *ff_it, properties);
    }
  }

  // Remove the properties of the removed face.
  properties.erase(other_face);
}