diff --git a/Apollonius_graph_2/include/CGAL/Parabola_segment_2.h b/Apollonius_graph_2/include/CGAL/Parabola_segment_2.h
index 20883169ba7..63faa2d9c92 100644
--- a/Apollonius_graph_2/include/CGAL/Parabola_segment_2.h
+++ b/Apollonius_graph_2/include/CGAL/Parabola_segment_2.h
@@ -63,7 +63,7 @@ struct Parabola_segment_2 : public Parabola_2< Gt >
   }
 
   int compute_k(const FT tt, const FT STEP) const {
-    return int(CGAL::to_double(CGAL::sqrt(tt / STEP)));
+    return int(CGAL::to_double(CGAL::approximate_sqrt(tt / STEP)));
   }
 
   // s0 and s1 define a desired drawing "range"
diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/distance.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/distance.h
index 2bded6d6373..29d4d5810bb 100644
--- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/distance.h
+++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/distance.h
@@ -1689,8 +1689,9 @@ bounded_error_squared_Hausdorff_distance_impl(const TriangleMesh1& tm1,
       // Thus, subdivision can only decrease the min, and the upper bound.
       Local_bounds<Kernel, Face_handle_1, Face_handle_2> bounds(triangle_bounds.upper);
 
-      // Ensure 'uface' is initialized in case the upper bound is not changed by the subdivision
+      // Ensure 'lface' and 'uface' are initialized in case the bounds are not changed by the subdivision
       bounds.tm2_uface = triangle_bounds.tm2_uface;
+      bounds.tm2_lface = triangle_bounds.tm2_lface;
 
       TM2_hd_traits traversal_traits_tm2(sub_t1_bbox, tm2, vpm2, bounds, global_bounds, infinity_value);
       tm2_tree.traversal_with_priority(sub_triangles[i], traversal_traits_tm2);
diff --git a/Polygon_mesh_processing/test/Polygon_mesh_processing/CMakeLists.txt b/Polygon_mesh_processing/test/Polygon_mesh_processing/CMakeLists.txt
index 51e2e1491a2..55ebce6c368 100644
--- a/Polygon_mesh_processing/test/Polygon_mesh_processing/CMakeLists.txt
+++ b/Polygon_mesh_processing/test/Polygon_mesh_processing/CMakeLists.txt
@@ -72,6 +72,7 @@ create_single_source_cgal_program("test_isolevel_refinement.cpp")
 create_single_source_cgal_program("test_corefinement_nm_bo.cpp")
 create_single_source_cgal_program("test_corefinement_cavities.cpp")
 create_single_source_cgal_program("issue_8730.cpp")
+create_single_source_cgal_program("issue_7164.cpp")
 #  create_single_source_cgal_program("test_pmp_repair_self_intersections.cpp")
 
 find_package(Eigen3 3.2.0 QUIET) #(requires 3.2.0 or greater)
diff --git a/Polygon_mesh_processing/test/Polygon_mesh_processing/issue_7164.cpp b/Polygon_mesh_processing/test/Polygon_mesh_processing/issue_7164.cpp
new file mode 100644
index 00000000000..de4e16a74e6
--- /dev/null
+++ b/Polygon_mesh_processing/test/Polygon_mesh_processing/issue_7164.cpp
@@ -0,0 +1,61 @@
+#include <CGAL/Surface_mesh.h>
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+#include <CGAL/Polygon_mesh_processing/polygon_soup_to_polygon_mesh.h>
+#include <CGAL/Polygon_mesh_processing/distance.h>
+#include <CGAL/Random.h>
+using Kernel = CGAL::Exact_predicates_inexact_constructions_kernel;
+using Point_3 = Kernel::Point_3;
+using Mesh = CGAL::Surface_mesh<Point_3>;
+namespace PMP = CGAL::Polygon_mesh_processing;
+
+int main(/*int argc, char** argv*/)
+{
+  // A simple triangle
+  std::vector<Point_3> pts_A;
+  std::vector<std::vector<size_t>> trs_A;
+  pts_A.emplace_back( 0.26641936088212415, 0.2664193608821242, 0.73358063911787585);
+  pts_A.emplace_back(-0.14011519816541251, 0.6017979969632727, 1.1810107045967466);
+  pts_A.emplace_back(-0.14011519816541279,-0.1810107045967464, 0.39820200303672726);
+  trs_A.emplace_back(std::vector<size_t>{0,1,2});
+  Mesh A;
+  PMP::polygon_soup_to_polygon_mesh(pts_A, trs_A, A);
+
+  // An open tetrahedron
+  std::vector<Point_3> pts_B;
+  std::vector<std::vector<size_t>> trs_B;
+  pts_B.emplace_back(0,0,0);
+  pts_B.emplace_back(1,1,0);
+  pts_B.emplace_back(1,0,1);
+  pts_B.emplace_back(0,1,1);
+  trs_B.emplace_back(std::vector<size_t>{0,1,2});
+  trs_B.emplace_back(std::vector<size_t>{3,1,0});
+  trs_B.emplace_back(std::vector<size_t>{3,2,1});
+  Mesh B;
+  PMP::polygon_soup_to_polygon_mesh(pts_B, trs_B, B);
+
+  double bound = 0.01 * 0.42149467833714593;
+  PMP::bounded_error_Hausdorff_distance<CGAL::Sequential_tag>(A, B, bound);
+  PMP::bounded_error_Hausdorff_distance<CGAL::Sequential_tag>(B, A, bound);
+
+  // The bug was possible with closed models
+  std::vector<Point_3> pts_C;
+  std::vector<std::vector<size_t>> trs_C;
+  pts_C.emplace_back(0,0,0);
+  pts_C.emplace_back(1,1,0);
+  pts_C.emplace_back(1,0,1);
+  pts_C.emplace_back(0,1,1);
+  pts_C.emplace_back(0.75,0.75,0);
+  trs_C.emplace_back(std::vector<size_t>{0,1,2});
+  trs_C.emplace_back(std::vector<size_t>{3,1,0});
+  trs_C.emplace_back(std::vector<size_t>{3,2,1});
+  trs_C.emplace_back(std::vector<size_t>{0,2,4});
+  trs_C.emplace_back(std::vector<size_t>{3,0,4});
+  trs_C.emplace_back(std::vector<size_t>{3,4,2});
+  Mesh C;
+  PMP::polygon_soup_to_polygon_mesh(pts_C, trs_C, C);
+
+  PMP::bounded_error_Hausdorff_distance<CGAL::Sequential_tag>(A, C, bound);
+  PMP::bounded_error_Hausdorff_distance<CGAL::Sequential_tag>(C, A, bound);
+
+  return EXIT_SUCCESS;
+}