More DC tests

Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/implicit_shapes_helper.h

@@ -61,26 +61,6 @@ box(const typename K::Point_3& b,
   return inside ? - d : d;
 }
 
-// template <typename K>
-// typename K::FT
-// disk(const typename K::Point_3& c,
-//      const typename K::Vector_3& n,
-//      const typename K::FT r,
-//      const typename K::Point_3& q)
-// {
-//   typename K::Plane_3 pl(c, n);
-//   typename K::Point_3 pq = pl.projection(q);
-
-//   typename K::FT sq_dpl = CGAL::squared_distance(q, pl);
-
-//   if(CGAL::squared_distance(pq, c) < CGAL::square(r))
-//     return CGAL::approximate_sqrt(sq_dpl);
-//   else
-//     return CGAL::approximate_sqrt(CGAL::square(CGAL::approximate_sqrt(CGAL::squared_distance(pq, c)) - r) + sq_dpl);
-// }
-
-
-
 // p is the center of the base disk
 // q is the center of the top disk
 template<typename K>
@@ -95,7 +75,6 @@ infinite_cylinder(const typename K::Point_3& b,
   return CGAL::approximate_sqrt(CGAL::squared_distance(pq, b)) - r;
 }
 
-
 // c is the center of the torus
 // n is the normal of the plane containing all centers of the tube
 // r is the small radius

Isosurfacing_3/test/Isosurfacing_3/test_dual_contouring.cpp

@@ -17,15 +17,17 @@
 #include <CGAL/Polygon_mesh_processing/triangulate_faces.h>
 #include <CGAL/Image_3.h>
 #include <CGAL/boost/graph/IO/OFF.h>
+#include <CGAL/Random.h>
 
 #include <array>
 #include <deque>
 #include <iostream>
+#include <functional>
 
 namespace IS = CGAL::Isosurfacing;
 
 template <typename K>
-void test_cube()
+void test_implicit_shape(std::function<typename K::FT(const typename K::Point_3&)> implicit_function)
 {
   using FT = typename K::FT;
   using Point = typename K::Point_3;
@@ -40,40 +42,6 @@ void test_cube()
   using Point_range = std::vector<Point>;
   using Polygon_range = std::vector<std::array<std::size_t, 4> >;
 
-  using Mesh = CGAL::Surface_mesh<Point>;
-
-auto implicit_function = [](const Point& q) -> FT
-{
-  return IS::Shapes::box<K>(Point(-1,-1,-1), Point(1,1,1), q);
-
-  // ---
-  auto cyl = [](const Point& q) { return IS::Shapes::infinite_cylinder<K>(Point(0,0,0), Vector(0,0,1), 0.5, q); };
-  auto cube = [](const Point& q) { return IS::Shapes::box<K>(Point(-0.5,-0.5,-0.5), Point(0.5,0.5,0.5), q); };
-  auto cyl_and_cube = [&](const Point& q) { return IS::Shapes::shape_union<K>(cyl, cube, q); };
-
-  auto sphere = [](const Point& q) { return IS::Shapes::sphere<K>(Point(0,0,0.5), 1, q); };
-  return IS::Shapes::shape_difference<K>(cyl_and_cube, sphere, q);
-
-  // ---
-  auto box_1 = [](const Point& q) { return IS::Shapes::box<K>(Point(0,0,0), Point(1,1,1), q); };
-  auto box_2 = [](const Point& q) { return IS::Shapes::box<K>(Point(0.5,0.5,0.5), Point(1.5,1.5,1.5), q); };
-  return IS::Shapes::shape_union<K>(box_1, box_2, q);
-
-  // ---
-  return IS::Shapes::box<K>(Point(0,0,0), Point(1,2,3), q);
-
-  // ---
-  return IS::Shapes::torus<K>(Point(0,0,0), Vector(0,0,1), 0.25, 1, q);
-
-  // ---
-  return IS::Shapes::infinite_cylinder<K>(Point(0,0,0), Vector(1,1,1), 1, q);
-
-  // ---
-  auto sphere_1 = [](const Point& q) { return IS::Shapes::sphere<K>(Point(0,0,0), 1, q); };
-  auto sphere_2 = [](const Point& q) { return IS::Shapes::sphere<K>(Point(0,0,0.5), 1, q); };
-  return IS::Shapes::shape_union<K>(sphere_1, sphere_2, q);
-};
-
   const CGAL::Bbox_3 bbox = {-2., -2., -2., 2., 2., 2.};
   const Vector spacing { 0.05, 0.05, 0.05 };
   const FT isovalue = 0;
@@ -83,13 +51,13 @@ auto implicit_function = [](const Point& q) -> FT
   std::cout << "Kernel: " << typeid(K).name() << std::endl;
 
   Grid grid { bbox, spacing };
-  Values values {implicit_function , grid };
+  Values values { implicit_function, grid };
   Gradients gradients { values, 0.01 * spacing[0] };
   Domain domain { grid, values, gradients };
 
   Point_range points;
   Polygon_range polygons;
-  IS::dual_contouring<CGAL::Parallel_if_available_tag>(domain, isovalue, points, polygons, CGAL::parameters::do_not_triangulate_faces(true)); // if you change that, change the array
+  IS::dual_contouring<CGAL::Parallel_if_available_tag>(domain, isovalue, points, polygons, CGAL::parameters::do_not_triangulate_faces(true));
 
   std::cout << "Output #vertices: " << points.size() << std::endl;
   std::cout << "Output #polygons: " << polygons.size() << std::endl;
@@ -121,8 +89,142 @@ auto implicit_function = [](const Point& q) -> FT
   std::cout << "Check for degenerate faces..." << std::endl;
   assert(!has_degenerate_faces(mesh));
 
-  std::cout << "Volume: " << CGAL::Polygon_mesh_processing::volume(mesh) << std::endl;
+  const FT computed_volume = CGAL::Polygon_mesh_processing::volume(mesh);
-  assert(CGAL::abs(CGAL::Polygon_mesh_processing::volume(mesh) - FT(8)) < 1e-2);
+  std::cout << "Computed Volume: " << computed_volume << std::endl;
+
+  // Brute force estimation of the volume
+  std::cout << "Estimating volume using Monte Carlo..." << std::endl;
+  const std::size_t num_samples = 1000000;
+  CGAL::Random rng(0);
+
+  std::size_t inside_count = 0;
+  for (std::size_t i = 0; i < num_samples; ++i) {
+    Point sample(rng.get_double(bbox.xmin(), bbox.xmax()),
+                 rng.get_double(bbox.ymin(), bbox.ymax()),
+                 rng.get_double(bbox.zmin(), bbox.zmax()));
+    if (implicit_function(sample) <= isovalue) {
+      ++inside_count;
+    }
+  }
+
+  const FT bbox_volume = (bbox.xmax() - bbox.xmin()) * (bbox.ymax() - bbox.ymin()) * (bbox.zmax() - bbox.zmin());
+  const FT monte_carlo_volume = bbox_volume * static_cast<FT>(inside_count) / static_cast<FT>(num_samples);
+
+  std::cout << "Monte Carlo Estimated Volume: " << monte_carlo_volume << std::endl;
+
+  // Compare computed volume with Monte Carlo estimation
+  const FT tolerance = 1e-2 * bbox_volume;
+  assert(CGAL::abs(computed_volume - monte_carlo_volume) < tolerance);
+  std::cout << "Volume check passed!" << std::endl;
+}
+
+template <typename K>
+void test_box()
+{
+  auto box_function = [](const typename K::Point_3& q) {
+    return IS::Shapes::box<K>(typename K::Point_3(-1, -1, -1), typename K::Point_3(1, 1, 1), q);
+  };
+  test_implicit_shape<K>(box_function);
+}
+
+template <typename K>
+void test_union_of_boxes()
+{
+  auto union_of_boxes_function = [](const typename K::Point_3& q) {
+    auto box_1 = [](const typename K::Point_3& q) {
+      return IS::Shapes::box<K>(typename K::Point_3(0, 0, 0), typename K::Point_3(1, 1, 1), q);
+    };
+    auto box_2 = [](const typename K::Point_3& q) {
+      return IS::Shapes::box<K>(typename K::Point_3(0.5, 0.5, 0.5), typename K::Point_3(1.5, 1.5, 1.5), q);
+    };
+    return IS::Shapes::shape_union<K>(box_1, box_2, q);
+  };
+  test_implicit_shape<K>(union_of_boxes_function);
+}
+
+template <typename K>
+void test_union_of_boxes_minus_sphere()
+{
+  auto b1_b2_ms1_function = [](const typename K::Point_3& q) {
+    auto box_1 = [](const typename K::Point_3& q) {
+      return IS::Shapes::box<K>(typename K::Point_3(0, 0, 0), typename K::Point_3(1, 1, 1), q);
+    };
+    auto box_2 = [](const typename K::Point_3& q) {
+      return IS::Shapes::box<K>(typename K::Point_3(0.5, 0.5, 0.5), typename K::Point_3(1.5, 1.5, 1.5), q);
+    };
+    auto b1_and_b2 = [&](const typename K::Point_3& q) {
+      return IS::Shapes::shape_union<K>(box_1, box_2, q);
+    };
+    auto s1 = [](const typename K::Point_3& q) {
+      return IS::Shapes::sphere<K>(typename K::Point_3(0, 0, 0.5), 1, q);
+    };
+    return IS::Shapes::shape_difference<K>(b1_and_b2, s1, q);
+  };
+  test_implicit_shape<K>(b1_b2_ms1_function);
+}
+
+template <typename K>
+void test_union_of_spheres()
+{
+  auto union_of_spheres_function = [](const typename K::Point_3& q) {
+    auto sphere_1 = [](const typename K::Point_3& q) {
+      return IS::Shapes::sphere<K>(typename K::Point_3(0, 0, 0), 1, q);
+    };
+    auto sphere_2 = [](const typename K::Point_3& q) {
+      return IS::Shapes::sphere<K>(typename K::Point_3(0, 0, 0.5), 1, q);
+    };
+    return IS::Shapes::shape_union<K>(sphere_1, sphere_2, q);
+  };
+  test_implicit_shape<K>(union_of_spheres_function);
+}
+
+template <typename K>
+void test_torus()
+{
+  auto torus_function = [](const typename K::Point_3& q) {
+    return IS::Shapes::torus<K>(typename K::Point_3(0, 0, 0), typename K::Vector_3(0, 0, 1), 0.25, 1, q);
+  };
+  test_implicit_shape<K>(torus_function);
+}
+
+// template <typename K>
+// void test_infinite_cylinder()
+// {
+//   auto cylinder_function = [](const typename K::Point_3& q) {
+//     return IS::Shapes::infinite_cylinder<K>(typename K::Point_3(0, 0, 0), typename K::Vector_3(1, 1, 1), 1, q);
+//   };
+//   test_implicit_shape<K>(cylinder_function);
+// }
+
+// template <typename K>
+// void test_cylinder_and_cube_minus_sphere()
+// {
+//   auto cyl_cube_function = [](const typename K::Point_3& q) {
+//     auto cyl = [](const typename K::Point_3& q) {
+//       return IS::Shapes::infinite_cylinder<K>(typename K::Point_3(0, 0, 0), typename K::Vector_3(0, 0, 1), 0.5, q);
+//     };
+//     auto cube = [](const typename K::Point_3& q) {
+//       return IS::Shapes::box<K>(typename K::Point_3(-0.5, -0.5, -0.5), typename K::Point_3(0.5, 0.5, 0.5), q);
+//     };
+//     auto cyl_and_cube = [&](const typename K::Point_3& q) {
+//       return IS::Shapes::shape_union<K>(cyl, cube, q);
+//     };
+//     auto sphere = [](const typename K::Point_3& q) {
+//       return IS::Shapes::sphere<K>(typename K::Point_3(0, 0, 0.5), 1, q);
+//     };
+//     return IS::Shapes::shape_difference<K>(cyl_and_cube, sphere, q);
+//   };
+//   test_implicit_shape<K>(cyl_cube_function);
+// }
+
+template <typename K>
+void test_implicit_shapes()
+{
+  test_box<K>();
+  test_union_of_boxes<K>();
+  test_union_of_boxes_minus_sphere<K>();
+  test_union_of_spheres<K>();
+  test_torus<K>();
 }
 
 template <typename K>
@@ -190,10 +292,10 @@ void test_image()
 
 int main(int, char**)
 {
-  test_cube<CGAL::Simple_cartesian<double> >();
+  test_implicit_shapes<CGAL::Simple_cartesian<double>>();
-  test_image<CGAL::Simple_cartesian<double> >();
+  test_image<CGAL::Simple_cartesian<double>>();
 
-  test_cube<CGAL::Exact_predicates_inexact_constructions_kernel>();
+  test_implicit_shapes<CGAL::Exact_predicates_inexact_constructions_kernel>();
   test_image<CGAL::Exact_predicates_inexact_constructions_kernel>();
 
   std::cout << "Done" << std::endl;