From 76dbf05e36f8d0fa641821ad32f95920dea29f2d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Mon, 24 Mar 2025 15:34:21 +0100 Subject: [PATCH] More DC tests --- .../internal/implicit_shapes_helper.h | 21 -- .../Isosurfacing_3/test_dual_contouring.cpp | 186 ++++++++++++++---- 2 files changed, 144 insertions(+), 63 deletions(-) diff --git a/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/implicit_shapes_helper.h b/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/implicit_shapes_helper.h index dd613d93fc8..83d2276b0f7 100644 --- a/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/implicit_shapes_helper.h +++ b/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::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 @@ -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 diff --git a/Isosurfacing_3/test/Isosurfacing_3/test_dual_contouring.cpp b/Isosurfacing_3/test/Isosurfacing_3/test_dual_contouring.cpp index c37cda9130e..5c717d122be 100644 --- a/Isosurfacing_3/test/Isosurfacing_3/test_dual_contouring.cpp +++ b/Isosurfacing_3/test/Isosurfacing_3/test_dual_contouring.cpp @@ -17,15 +17,17 @@ #include #include #include +#include #include #include #include +#include namespace IS = CGAL::Isosurfacing; template -void test_cube() +void test_implicit_shape(std::function 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; using Polygon_range = std::vector >; - using Mesh = CGAL::Surface_mesh; - -auto implicit_function = [](const Point& q) -> FT -{ - return IS::Shapes::box(Point(-1,-1,-1), Point(1,1,1), q); - - // --- - auto cyl = [](const Point& q) { return IS::Shapes::infinite_cylinder(Point(0,0,0), Vector(0,0,1), 0.5, q); }; - auto cube = [](const Point& q) { return IS::Shapes::box(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(cyl, cube, q); }; - - auto sphere = [](const Point& q) { return IS::Shapes::sphere(Point(0,0,0.5), 1, q); }; - return IS::Shapes::shape_difference(cyl_and_cube, sphere, q); - - // --- - auto box_1 = [](const Point& q) { return IS::Shapes::box(Point(0,0,0), Point(1,1,1), q); }; - auto box_2 = [](const Point& q) { return IS::Shapes::box(Point(0.5,0.5,0.5), Point(1.5,1.5,1.5), q); }; - return IS::Shapes::shape_union(box_1, box_2, q); - - // --- - return IS::Shapes::box(Point(0,0,0), Point(1,2,3), q); - - // --- - return IS::Shapes::torus(Point(0,0,0), Vector(0,0,1), 0.25, 1, q); - - // --- - return IS::Shapes::infinite_cylinder(Point(0,0,0), Vector(1,1,1), 1, q); - - // --- - auto sphere_1 = [](const Point& q) { return IS::Shapes::sphere(Point(0,0,0), 1, q); }; - auto sphere_2 = [](const Point& q) { return IS::Shapes::sphere(Point(0,0,0.5), 1, q); }; - return IS::Shapes::shape_union(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(domain, isovalue, points, polygons, CGAL::parameters::do_not_triangulate_faces(true)); // if you change that, change the array + IS::dual_contouring(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; - assert(CGAL::abs(CGAL::Polygon_mesh_processing::volume(mesh) - FT(8)) < 1e-2); + const FT computed_volume = CGAL::Polygon_mesh_processing::volume(mesh); + 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(inside_count) / static_cast(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 +void test_box() +{ + auto box_function = [](const typename K::Point_3& q) { + return IS::Shapes::box(typename K::Point_3(-1, -1, -1), typename K::Point_3(1, 1, 1), q); + }; + test_implicit_shape(box_function); +} + +template +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(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(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(box_1, box_2, q); + }; + test_implicit_shape(union_of_boxes_function); +} + +template +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(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(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(box_1, box_2, q); + }; + auto s1 = [](const typename K::Point_3& q) { + return IS::Shapes::sphere(typename K::Point_3(0, 0, 0.5), 1, q); + }; + return IS::Shapes::shape_difference(b1_and_b2, s1, q); + }; + test_implicit_shape(b1_b2_ms1_function); +} + +template +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(typename K::Point_3(0, 0, 0), 1, q); + }; + auto sphere_2 = [](const typename K::Point_3& q) { + return IS::Shapes::sphere(typename K::Point_3(0, 0, 0.5), 1, q); + }; + return IS::Shapes::shape_union(sphere_1, sphere_2, q); + }; + test_implicit_shape(union_of_spheres_function); +} + +template +void test_torus() +{ + auto torus_function = [](const typename K::Point_3& q) { + return IS::Shapes::torus(typename K::Point_3(0, 0, 0), typename K::Vector_3(0, 0, 1), 0.25, 1, q); + }; + test_implicit_shape(torus_function); +} + +// template +// void test_infinite_cylinder() +// { +// auto cylinder_function = [](const typename K::Point_3& q) { +// return IS::Shapes::infinite_cylinder(typename K::Point_3(0, 0, 0), typename K::Vector_3(1, 1, 1), 1, q); +// }; +// test_implicit_shape(cylinder_function); +// } + +// template +// 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(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(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(cyl, cube, q); +// }; +// auto sphere = [](const typename K::Point_3& q) { +// return IS::Shapes::sphere(typename K::Point_3(0, 0, 0.5), 1, q); +// }; +// return IS::Shapes::shape_difference(cyl_and_cube, sphere, q); +// }; +// test_implicit_shape(cyl_cube_function); +// } + +template +void test_implicit_shapes() +{ + test_box(); + test_union_of_boxes(); + test_union_of_boxes_minus_sphere(); + test_union_of_spheres(); + test_torus(); } template @@ -190,10 +292,10 @@ void test_image() int main(int, char**) { - test_cube >(); - test_image >(); + test_implicit_shapes>(); + test_image>(); - test_cube(); + test_implicit_shapes(); test_image(); std::cout << "Done" << std::endl;