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

template<typename KernelT>
void Face_decimation<KernelT>::identify_faces_at_boundaries(const Mesh &mesh, 
  Mesh_properties &properties, std::array<FT, 3> &min_c, std::array<FT, 3> &max_c)
{
  // Look for the minimum and maximum coordinates.
  constexpr FT big = std::numeric_limits<FT>::max();
  constexpr FT small = std::numeric_limits<FT>::lowest();

  max_c = {small, small, small};
  min_c = {big, big, big};

  {
    typename Mesh::Vertex_range::iterator v_it = mesh.vertices().begin(),
      v_ite = mesh.vertices().end();

    for(; v_it != v_ite; ++v_it)
    {
      const Point &p = mesh.point(*v_it);

      for(int i = 0; i < 3; ++i)
      {
        if(max_c[i] < p[i]) max_c[i] = p[i];
        if(min_c[i] > p[i]) min_c[i] = p[i];
      }
    }
  }

  // Identify faces on boundaries.
  {
    const FT factor = 10;
    typename Mesh::Face_range::iterator f_it = mesh.faces().begin(),
      f_ite = mesh.faces().end();

    for(; f_it != f_ite; ++f_it)
    {
      // Test if all the vertices of the current surface are located on a boundary.
      Boundary_type state = Boundary_type::BOUNDARY_NONE;
      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);

      std::map<Vertex_index, std::set<Boundary_type>> v_boundaries;

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

        if(nearly_equal(p[0], min_c[0], factor)/* || (std::abs(p[0]-min_c[0]) <= 0.01)*/)
          v_boundaries[*vf_it].insert(Boundary_type::BOUNDARY_LEFT);

        if(nearly_equal(p[0], max_c[0], factor)/* || (std::abs(p[0]-max_c[0]) <= 0.01)*/)
          v_boundaries[*vf_it].insert(Boundary_type::BOUNDARY_RIGHT);

        if(nearly_equal(p[1], min_c[1], factor)/* || (std::abs(p[1]-min_c[1]) <= 0.01)*/)
          v_boundaries[*vf_it].insert(Boundary_type::BOUNDARY_FRONT);

        if(nearly_equal(p[1], max_c[1], factor)/* || (std::abs(p[1]-max_c[1]) <= 0.01)*/)
          v_boundaries[*vf_it].insert(Boundary_type::BOUNDARY_REAR);

        if(nearly_equal(p[2], min_c[2], factor)/* || (std::abs(p[2]-min_c[2]) <= 0.01)*/)
          v_boundaries[*vf_it].insert(Boundary_type::BOUNDARY_BOTTOM);

        if(nearly_equal(p[2], max_c[2], factor)/* || (std::abs(p[2]-max_c[2]) <= 0.01)*/)
          v_boundaries[*vf_it].insert(Boundary_type::BOUNDARY_TOP);

        if(v_boundaries.find(*vf_it) == v_boundaries.end())
          v_boundaries[*vf_it].insert(Boundary_type::BOUNDARY_NONE);
      } // Loop on vertices of current face.

      // For each vertex look for a common boundary (if any).
      {
        typename std::map<Vertex_index, std::set<Boundary_type>>::const_iterator
          it_map = v_boundaries.begin(), ite_map = v_boundaries.end();
        std::set<Boundary_type> v_intersection, v_start = {Boundary_type::BOUNDARY_LEFT, 
          Boundary_type::BOUNDARY_RIGHT, Boundary_type::BOUNDARY_FRONT, 
          Boundary_type::BOUNDARY_REAR, Boundary_type::BOUNDARY_BOTTOM,
          Boundary_type::BOUNDARY_TOP};

        while(it_map != ite_map)
        {
          std::set_intersection(it_map->second.begin(), it_map->second.end(),
                                v_start.begin(), v_start.end(),
                                std::inserter(v_intersection, v_intersection.end()));

          v_start = v_intersection;
          v_intersection.clear();
          ++it_map;
        }

        if(!v_boundaries.empty())
        {
          if(v_start.size() == 1)
            state = *v_start.begin();
          else if(!v_start.empty())
            throw std::runtime_error("Face_decimation::identify_faces_at_boundaries: face vertices belong to multiple boundaries");
        }
      }

      properties[*f_it].boundary = state;
    } // Loop on faces.
  }
}