From f0e2bb33c91a85b317af2d7a0add9104d166bb5f Mon Sep 17 00:00:00 2001 From: Julian Stahl Date: Thu, 19 Dec 2024 19:01:06 +0100 Subject: [PATCH] Add topology verifier --- .../test/Isosurfacing_3/CMakeLists.txt | 2 + .../test/Isosurfacing_3/test_topology.cpp | 50 ---- .../test/Isosurfacing_3/test_util.h | 66 +++++ .../test/Isosurfacing_3/verifier.cpp | 238 ++++++++++++++++++ 4 files changed, 306 insertions(+), 50 deletions(-) create mode 100644 Isosurfacing_3/test/Isosurfacing_3/verifier.cpp diff --git a/Isosurfacing_3/test/Isosurfacing_3/CMakeLists.txt b/Isosurfacing_3/test/Isosurfacing_3/CMakeLists.txt index 0138f1e7fc1..67116810931 100644 --- a/Isosurfacing_3/test/Isosurfacing_3/CMakeLists.txt +++ b/Isosurfacing_3/test/Isosurfacing_3/CMakeLists.txt @@ -12,6 +12,7 @@ include(CGAL_TBB_support) create_single_source_cgal_program("test_marching_cubes.cpp") create_single_source_cgal_program("test_topology.cpp") +create_single_source_cgal_program("verifier.cpp") create_single_source_cgal_program("test_tmc_csg.cpp") if(TARGET CGAL::Eigen3_support) @@ -45,4 +46,5 @@ endif() if(TARGET CGAL::TBB_support) target_link_libraries(test_marching_cubes PRIVATE CGAL::TBB_support) target_link_libraries(test_topology PRIVATE CGAL::TBB_support) + target_link_libraries(verifier PRIVATE CGAL::TBB_support) endif() diff --git a/Isosurfacing_3/test/Isosurfacing_3/test_topology.cpp b/Isosurfacing_3/test/Isosurfacing_3/test_topology.cpp index 5ca139605b7..eea9463ff8b 100644 --- a/Isosurfacing_3/test/Isosurfacing_3/test_topology.cpp +++ b/Isosurfacing_3/test/Isosurfacing_3/test_topology.cpp @@ -27,33 +27,6 @@ namespace IS = CGAL::Isosurfacing; -template -std::size_t connected_components(PolygonMesh& mesh) -{ - using face_descriptor = typename boost::graph_traits::face_descriptor; - - auto fccmap = mesh.template add_property_map("f:CC").first; - - return CGAL::Polygon_mesh_processing::connected_components(mesh, fccmap); -} - -template -std::size_t euler_characteristic(const PolygonMesh& mesh) -{ - return mesh.number_of_vertices() - mesh.number_of_edges() + mesh.number_of_faces(); -} - -template -std::size_t boundary_components(const PolygonMesh& mesh) -{ - using halfedge_descriptor = typename boost::graph_traits::halfedge_descriptor; - std::vector border_cycles; - - CGAL::Polygon_mesh_processing::extract_boundary_cycles(mesh, std::back_inserter(border_cycles)); - - return border_cycles.size(); -} - template bool check_closed_not_empty(const Grid& grid, const IS::Interpolated_discrete_values_3& values, const typename Grid::Geom_traits::FT iso = 0) { @@ -441,29 +414,6 @@ FT generate_predefined_plane(const PlaneCase topology_case, std::array& v return iso; } -template -void write_debug_grid(const Domain& domain, const std::string& filename) { - using Point = typename Domain::Geom_traits::Point_3; - using Mesh = CGAL::Surface_mesh; - - Mesh debug_grid; - auto debug_grid_creator = [&](const typename Domain::cell_descriptor& c) - { - std::vector cell_vertices; - for (const auto& v : domain.cell_vertices(c)) { - cell_vertices.push_back(debug_grid.add_vertex(domain.point(v))); - } - debug_grid.add_face(cell_vertices[6], cell_vertices[2], cell_vertices[0], cell_vertices[4]); - debug_grid.add_face(cell_vertices[1], cell_vertices[3], cell_vertices[7], cell_vertices[5]); - debug_grid.add_face(cell_vertices[0], cell_vertices[1], cell_vertices[5], cell_vertices[4]); - debug_grid.add_face(cell_vertices[6], cell_vertices[7], cell_vertices[3], cell_vertices[2]); - debug_grid.add_face(cell_vertices[2], cell_vertices[3], cell_vertices[1], cell_vertices[0]); - debug_grid.add_face(cell_vertices[4], cell_vertices[5], cell_vertices[7], cell_vertices[6]); - }; - domain.template for_each_cell(debug_grid_creator); - CGAL::IO::write_OFF(filename, debug_grid); -} - template void compare_tmc_mc_trilinear(const std::array& case_values, typename KERNEL::FT iso) { diff --git a/Isosurfacing_3/test/Isosurfacing_3/test_util.h b/Isosurfacing_3/test/Isosurfacing_3/test_util.h index 0e819b37150..e311a3f2f0c 100644 --- a/Isosurfacing_3/test/Isosurfacing_3/test_util.h +++ b/Isosurfacing_3/test/Isosurfacing_3/test_util.h @@ -7,6 +7,7 @@ #include #include #include +#include #include #include @@ -59,6 +60,48 @@ bool has_degenerate_faces(PolygonMesh& pmesh) return !dfs.empty(); } +template +int connected_components(PolygonMesh& mesh) +{ + using face_descriptor = typename boost::graph_traits::face_descriptor; + + auto fccmap = mesh.template add_property_map("f:CC").first; + + return CGAL::Polygon_mesh_processing::connected_components(mesh, fccmap); +} + +template +int euler_characteristic(const PolygonMesh& mesh) +{ + return mesh.number_of_vertices() - mesh.number_of_edges() + mesh.number_of_faces(); +} + +template +int boundary_components(const PolygonMesh& mesh) +{ + using halfedge_descriptor = typename boost::graph_traits::halfedge_descriptor; + std::vector border_cycles; + + CGAL::Polygon_mesh_processing::extract_boundary_cycles(mesh, std::back_inserter(border_cycles)); + + return border_cycles.size(); +} + +template +int betti_0(PolygonMesh& mesh) { + return connected_components(mesh); +} + +template +int betti_2(PolygonMesh& mesh) { + return connected_components(mesh); // only for closed surfaces +} + +template +int betti_1(PolygonMesh& mesh) { + return betti_0(mesh) + betti_2(mesh) - euler_characteristic(mesh); +} + // template // bool check_mesh_distance(const PolygonMesh& m0, const PolygonMesh& m1) // { @@ -68,4 +111,27 @@ bool has_degenerate_faces(PolygonMesh& pmesh) // return true; // } +template +void write_debug_grid(const Domain& domain, const std::string& filename) { + using Point = typename Domain::Geom_traits::Point_3; + using Mesh = CGAL::Surface_mesh; + + Mesh debug_grid; + auto debug_grid_creator = [&](const typename Domain::cell_descriptor& c) + { + std::vector cell_vertices; + for (const auto& v : domain.cell_vertices(c)) { + cell_vertices.push_back(debug_grid.add_vertex(domain.point(v))); + } + debug_grid.add_face(cell_vertices[6], cell_vertices[2], cell_vertices[0], cell_vertices[4]); + debug_grid.add_face(cell_vertices[1], cell_vertices[3], cell_vertices[7], cell_vertices[5]); + debug_grid.add_face(cell_vertices[0], cell_vertices[1], cell_vertices[5], cell_vertices[4]); + debug_grid.add_face(cell_vertices[6], cell_vertices[7], cell_vertices[3], cell_vertices[2]); + debug_grid.add_face(cell_vertices[2], cell_vertices[3], cell_vertices[1], cell_vertices[0]); + debug_grid.add_face(cell_vertices[4], cell_vertices[5], cell_vertices[7], cell_vertices[6]); + }; + domain.template for_each_cell(debug_grid_creator); + CGAL::IO::write_OFF(filename, debug_grid); +} + #endif // CGAL_ISOSURFACING_TEST_UTIL_H diff --git a/Isosurfacing_3/test/Isosurfacing_3/verifier.cpp b/Isosurfacing_3/test/Isosurfacing_3/verifier.cpp new file mode 100644 index 00000000000..410dcfcd44c --- /dev/null +++ b/Isosurfacing_3/test/Isosurfacing_3/verifier.cpp @@ -0,0 +1,238 @@ + +#include "test_util.h" + +#include +#include + +#include +#include +#include +#include +#include +#include + +/* +Verifier for topological correctness +Based on: Topology Verification for Isosurface Extraction, Tiago Etiene +Data sets: [Marching Cubes cases] [Randomly generated grids] from http://liscustodio.github.io/C_MC33/ +*/ + +namespace IS = CGAL::Isosurfacing; + +template +void read_iso_volume(const std::string& filename, Grid& grid, Values& values) { + using Geom_traits = typename Grid::Geom_traits; + using FT = typename Geom_traits::FT; + using Iso_cuboid_3 = typename Geom_traits::Iso_cuboid_3; + typename Geom_traits::Construct_point_3 point = grid.geom_traits().construct_point_3_object(); + typename Geom_traits::Construct_iso_cuboid_3 iso_cuboid = grid.geom_traits().construct_iso_cuboid_3_object(); + + std::ifstream file(filename, std::ios::binary); + if (!file.is_open()) { + throw std::runtime_error("Cannot open file: " + filename); + } + + // Read the three dimensions + int nx, ny, nz; + file.read(reinterpret_cast(&nx), sizeof(int)); + file.read(reinterpret_cast(&ny), sizeof(int)); + file.read(reinterpret_cast(&nz), sizeof(int)); + + // Read the bounding box + // (xmin, xmax, ymin, ymax, zmin, zmax) + std::array bbox; + for (int i = 0; i < 6; ++i) { + file.read(reinterpret_cast(&bbox[i]), sizeof(float)); + } + + Iso_cuboid_3 span = iso_cuboid(point(bbox[0], bbox[2], bbox[4]), + point(bbox[1], bbox[3], bbox[5])); + grid = Grid { span, CGAL::make_array(nx, ny, nz) }; + + // Calculate total number of data points + const std::size_t total_points = nx * ny * nz; + + // Read the volume data (32-bit floats) + std::vector volume_data(total_points); + for (std::size_t i = 0; i < total_points; ++i) { + file.read(reinterpret_cast(&volume_data[i]), sizeof(float)); + std::cout << volume_data[i] << std::endl; + } + + for(int x=0; x> euler; + return euler; +} + +std::array read_betti(const std::string& filename) { + std::ifstream file(filename); + if (!file.is_open()) { + throw std::runtime_error("Cannot open file: " + filename); + } + int b0, b1; + file >> b0; + file >> b1; + return {b0, b1}; +} + +template +void verify_euler() { + using FT = typename K::FT; + using Point = typename K::Point_3; + using Vector = typename K::Vector_3; + + using Mesh = CGAL::Surface_mesh; + using Grid = IS::Cartesian_grid_3; + using Values = IS::Interpolated_discrete_values_3; + using Domain = IS::Marching_cubes_domain_3; + + using Point_range = std::vector; + using Polygon_range = std::vector >; + + using Mesh = CGAL::Surface_mesh; + + const std::size_t num_tests = 10000; + + for (std::size_t i = 0; i < num_tests; i++) { + + Grid grid; + Values values { grid }; + Domain domain { grid, values }; + + read_iso_volume("data/MarchingCubes_cases/Grids/" + std::to_string(i) + "-scalar_field.iso", grid, values); + + Point_range points; + Polygon_range triangles; + IS::marching_cubes(domain, 0, points, triangles, CGAL::parameters::use_topologically_correct_marching_cubes(true)); + + assert(points.size() && triangles.size()); + assert(!has_duplicate_points(points, triangles)); + assert(!has_duplicate_polygons(points, triangles)); + assert(!has_isolated_vertices(points, triangles)); + + assert(CGAL::Polygon_mesh_processing::is_polygon_soup_a_polygon_mesh(triangles)); + Mesh m; + CGAL::Polygon_mesh_processing::polygon_soup_to_polygon_mesh(points, triangles, m); + + const int euler = euler_characteristic(m); + const int solution = read_euler("data/MarchingCubes_cases/Grid_invariants/" + std::to_string(i) + "-euler_number.txt"); + + if (euler != solution) + std::cout << "error in test " << i << ": euler " << euler << " != " << solution << std::endl; + } +} + +template +void verify_betti() { + using FT = typename K::FT; + using Point = typename K::Point_3; + using Vector = typename K::Vector_3; + + using Mesh = CGAL::Surface_mesh; + using Grid = IS::Cartesian_grid_3; + using Values = IS::Interpolated_discrete_values_3; + using Domain = IS::Marching_cubes_domain_3; + + using Point_range = std::vector; + using Polygon_range = std::vector >; + + using Mesh = CGAL::Surface_mesh; + + const std::size_t num_tests = 10000; + + for (std::size_t i = 0; i < num_tests; i++) { + + Grid grid; + Values values { grid }; + Domain domain { grid, values }; + + read_iso_volume("data/Closed_Surfaces/Grids/" + std::to_string(i) + "-scalar_field.iso", grid, values); + + Point_range points; + Polygon_range triangles; + IS::marching_cubes(domain, 0, points, triangles, CGAL::parameters::use_topologically_correct_marching_cubes(true)); + + assert(points.size() && triangles.size()); + assert(!has_duplicate_points(points, triangles)); + assert(!has_duplicate_polygons(points, triangles)); + assert(!has_isolated_vertices(points, triangles)); + + assert(CGAL::Polygon_mesh_processing::is_polygon_soup_a_polygon_mesh(triangles)); + Mesh m; + CGAL::Polygon_mesh_processing::polygon_soup_to_polygon_mesh(points, triangles, m); + + const int b0 = betti_0(m); + const int b1 = betti_1(m); + const auto solution = read_betti("data/Closed_Surfaces/InvariantsGrid/" + std::to_string(i) + "-invariant_grid.txt"); + + if (b0 != solution[0]) + std::cout << "error in test " << i << ": b0 " << b0 << " != " << solution[0] << std::endl; + if (b1 != solution[1]) + std::cout << "error in test " << i << ": b1 " << b1 << " != " << solution[1] << std::endl; + } +} + +template +void compare_to_reference(const std::string& filename) { + using FT = typename K::FT; + using Point = typename K::Point_3; + using Vector = typename K::Vector_3; + + using Mesh = CGAL::Surface_mesh; + using Grid = IS::Cartesian_grid_3; + using Values = IS::Interpolated_discrete_values_3; + using Domain = IS::Marching_cubes_domain_3; + + using Point_range = std::vector; + using Polygon_range = std::vector >; + + using Mesh = CGAL::Surface_mesh; + + Grid grid; + Values values { grid }; + Domain domain { grid, values }; + + read_iso_volume(filename, grid, values); + + Point_range points; + Polygon_range triangles; + IS::marching_cubes(domain, 0, points, triangles, CGAL::parameters::use_topologically_correct_marching_cubes(true)); + + CGAL::IO::write_polygon_soup("verify_tmc.off", points, triangles, CGAL::parameters::stream_precision(17)); + + Grid grid_high_res { grid.span().min(), grid.span().max(), std::array{151, 151, 151} }; + IS::Value_function_3 values_high_res { values, grid_high_res }; + IS::Marching_cubes_domain_3> domain_high_res { grid_high_res, values_high_res }; + + Point_range points_high_res; + Polygon_range triangles_high_res; + IS::marching_cubes(domain_high_res, 0, points_high_res, triangles_high_res, CGAL::parameters::use_topologically_correct_marching_cubes(false)); + + CGAL::IO::write_polygon_soup("verify_reference.off", points_high_res, triangles_high_res); + + write_debug_grid(domain, "verify_cell.off"); +} + +int main(int, char**) +{ + using K = CGAL::Simple_cartesian; + using FT = typename K::FT; + + // verify_euler(); + // verify_betti(); + compare_to_reference("data/MarchingCubes_cases/Grids/" + std::to_string(100) + "-scalar_field.iso"); + // compare_to_reference("data/Closed_Surfaces/Grids/" + std::to_string(0) + "-scalar_field.iso"); +} \ No newline at end of file