diff --git a/Isosurfacing_3/test/Isosurfacing_3/CMakeLists.txt b/Isosurfacing_3/test/Isosurfacing_3/CMakeLists.txt index 37e2542ccbf..622ce58d187 100644 --- a/Isosurfacing_3/test/Isosurfacing_3/CMakeLists.txt +++ b/Isosurfacing_3/test/Isosurfacing_3/CMakeLists.txt @@ -7,6 +7,7 @@ project( Isosurfacing_3_tests ) find_package(CGAL REQUIRED) create_single_source_cgal_program( "test.cpp" ) +create_single_source_cgal_program( "test_marching_cubes.cpp" ) find_package(OpenMP) @@ -15,4 +16,6 @@ include(CGAL_TBB_support) if(TARGET CGAL::TBB_support) target_link_libraries(test PUBLIC CGAL::TBB_support) target_link_libraries(test PRIVATE OpenMP::OpenMP_CXX) + + target_link_libraries(test_marching_cubes PUBLIC CGAL::TBB_support) endif() diff --git a/Isosurfacing_3/test/Isosurfacing_3/test.cpp b/Isosurfacing_3/test/Isosurfacing_3/test.cpp index 9adb2ba2821..9d045540806 100644 --- a/Isosurfacing_3/test/Isosurfacing_3/test.cpp +++ b/Isosurfacing_3/test/Isosurfacing_3/test.cpp @@ -11,7 +11,6 @@ #include #include #include - #include #include "Timer.h" @@ -27,21 +26,23 @@ typedef tbb::concurrent_vector Point_range; typedef tbb::concurrent_vector> Polygon_range; int main() { - const Vector spacing(0.002f, 0.002f, 0.02f); + const Vector spacing(0.002, 0.002, 0.02); const CGAL::Bbox_3 bbox = {-1, -1, -1, 1, 1, 1}; auto sphere_function = [](const Point& point) { return std::sqrt(point.x() * point.x() + point.y() * point.y() + point.z() * point.z()); }; - CGAL::Isosurfacing::Implicit_domain< - Kernel, decltype(sphere_function), - decltype(CGAL::Isosurfacing::Default_gradient(sphere_function))> - implicit_domain({-1, -1, -1, 1, 1, 1}, spacing, sphere_function, - CGAL::Isosurfacing::Default_gradient( - sphere_function)); // TODO: this is ugly + typedef CGAL::Isosurfacing::Finite_difference_gradient Gradient; - Grid grid(2.f / spacing.x(), 2.f / spacing.y(), 2.f / spacing.z(), bbox); + CGAL::Isosurfacing::Implicit_domain implicit_domain( + {-1, -1, -1, 1, 1, 1}, spacing, sphere_function, Gradient(sphere_function, 0.0001)); // TODO: this is ugly + + const std::size_t nx = static_cast(2.0 / spacing.x()); + const std::size_t ny = static_cast(2.0 / spacing.y()); + const std::size_t nz = static_cast(2.0 / spacing.z()); + + Grid grid(nx, ny, nz, bbox); for (std::size_t x = 0; x < grid.xdim(); x++) { for (std::size_t y = 0; y < grid.ydim(); y++) { @@ -71,15 +72,15 @@ int main() { { ScopeTimer timer; - CGAL::Isosurfacing::make_triangle_mesh_using_marching_cubes(implicit_domain, 0.8f, points, - polygons); + CGAL::Isosurfacing::make_quad_mesh_using_dual_contouring(implicit_domain, 0.8f, points, + polygons); } // TODO: compare results with mesh_3 - //Mesh mesh; - //CGAL::Polygon_mesh_processing::polygon_soup_to_polygon_mesh(points, polygons, mesh); + // Mesh mesh; + // CGAL::Polygon_mesh_processing::polygon_soup_to_polygon_mesh(points, polygons, mesh); - //CGAL::IO::write_OFF("result.off", mesh); + // CGAL::IO::write_OFF("result.off", mesh); CGAL::IO::write_OFF("result.off", points, polygons); } diff --git a/Isosurfacing_3/test/Isosurfacing_3/test_marching_cubes.cpp b/Isosurfacing_3/test/Isosurfacing_3/test_marching_cubes.cpp new file mode 100644 index 00000000000..7f6553510c0 --- /dev/null +++ b/Isosurfacing_3/test/Isosurfacing_3/test_marching_cubes.cpp @@ -0,0 +1,96 @@ + +#include + +#include "test_util.h" + +#define WRITE_OFF + +struct Sphere_function { + FT operator()(const Point& point) const { + return std::sqrt(point.x() * point.x() + point.y() * point.y() + point.z() * point.z()); + } +}; + +template +void run(const Domain_& domain, const FT iso_value, Point_range& points, Polygon_range& polygons) { + CGAL::Isosurfacing::make_triangle_mesh_using_marching_cubes(domain, iso_value, points, + polygons); +} + +void test_implicit_sphere() { + const std::string test_name = "test_implicit_sphere()"; + + const Vector spacing(0.2, 0.2, 0.2); + const CGAL::Bbox_3 bbox = {-1, -1, -1, 1, 1, 1}; + + CGAL::Isosurfacing::Implicit_domain domain(bbox, spacing, + Sphere_function()); // TODO: this is ugly + + Point_range points; + Polygon_range polygons; + run(domain, 0.8, points, polygons); + +#ifdef WRITE_OFF + CGAL::IO::write_OFF(test_name + ".off", points, polygons); +#endif + + assert(is_polygon_mesh(polygons)); + Mesh m = to_mesh(points, polygons); + + assert(is_manifold(m)); + assert(!has_degenerate_faces(m)); + + std::cout << "Test passed: " << test_name << std::endl; +} + +void test_grid_sphere(const std::size_t n) { + const std::string test_name = "test_grid_sphere(" + std::to_string(n) + ")"; + + const CGAL::Bbox_3 bbox = {-1, -1, -1, 1, 1, 1}; + const Vector spacing(2.0 / (n - 1), 2.0 / (n - 1), 2.0 / (n - 1)); + + Sphere_function sphere_function; + + Grid grid(n, n, n, bbox); + + for (std::size_t x = 0; x < grid.xdim(); x++) { + for (std::size_t y = 0; y < grid.ydim(); y++) { + for (std::size_t z = 0; z < grid.zdim(); z++) { + + const Point pos(x * spacing.x() + bbox.xmin(), y * spacing.y() + bbox.ymin(), + z * spacing.z() + bbox.zmin()); + + grid.value(x, y, z) = sphere_function(pos); + } + } + } + + CGAL::Isosurfacing::Cartesian_grid_domain domain(grid); + + Point_range points; + Polygon_range polygons; + run(domain, 0.777, points, polygons); + +#ifdef WRITE_OFF + CGAL::IO::write_OFF(test_name + ".off", points, polygons); +#endif + + assert(is_polygon_mesh(polygons)); + Mesh m = to_mesh(points, polygons); + + assert(is_manifold(m)); + assert(!has_degenerate_faces(m)); + + std::cout << "Test passed: " << test_name << std::endl; +} + +int main() { + test_implicit_sphere(); + test_grid_sphere(2); + test_grid_sphere(3); + test_grid_sphere(10); + test_grid_sphere(11); + test_grid_sphere(100); + + std::cout << "All tests passed" << std::endl; +} diff --git a/Isosurfacing_3/test/Isosurfacing_3/test_util.h b/Isosurfacing_3/test/Isosurfacing_3/test_util.h new file mode 100644 index 00000000000..388772af99d --- /dev/null +++ b/Isosurfacing_3/test/Isosurfacing_3/test_util.h @@ -0,0 +1,70 @@ +#ifndef CGAL_ISOSURFACING_TEST_UTIL_H +#define CGAL_ISOSURFACING_TEST_UTIL_H + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel; +typedef typename Kernel::FT FT; +typedef typename Kernel::Vector_3 Vector; +typedef typename Kernel::Point_3 Point; + +typedef CGAL::Surface_mesh Mesh; +typedef CGAL::Cartesian_grid_3 Grid; + +typedef std::vector Point_range; +typedef std::vector> Polygon_range; + +namespace PMP = CGAL::Polygon_mesh_processing; + +bool has_duplicate_points(Point_range points, Polygon_range polygons) { + return PMP::merge_duplicate_points_in_polygon_soup(points, polygons) != 0; +} + +bool has_duplicate_polygons(Point_range points, Polygon_range polygons) { + return PMP::merge_duplicate_polygons_in_polygon_soup(points, polygons) != 0; +} + +bool has_isolated_vertices(Point_range points, Polygon_range polygons) { + return PMP::remove_isolated_points_in_polygon_soup(points, polygons) != 0; +} + +bool is_polygon_mesh(const Polygon_range& polygons) { + return PMP::is_polygon_soup_a_polygon_mesh(polygons); +} + +Mesh to_mesh(const Point_range& points, const Polygon_range& polygons) { + Mesh m; + PMP::polygon_soup_to_polygon_mesh(points, polygons, m); + return m; +} + +bool is_manifold(Mesh& m) { + return PMP::duplicate_non_manifold_vertices(m, CGAL::parameters::dry_run(true)) == 0; +} + +bool has_degenerate_faces(Mesh& m) { + return PMP::remove_connected_components_of_negligible_size( + m, CGAL::parameters::dry_run(true).area_threshold(std::numeric_limits::epsilon())) != 0; +} + +bool check_mesh_distance(const Mesh& m0, const Mesh& m1) { + auto dist = PMP::approximate_Hausdorff_distance( + m0, m1, CGAL::parameters::number_of_points_per_area_unit(4000)); + std::cout << dist << std::endl; + return true; +} + +#endif // CGAL_ISOSURFACING_TEST_UTIL_H \ No newline at end of file