mirror of https://github.com/CGAL/cgal
Add topology verifier
This commit is contained in:
Julian Stahl
parent
a534ef374a
commit
f0e2bb33c9
|
|
@ -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()
|
||||
|
|
|
|||
|
|
@ -27,33 +27,6 @@
|
|||
|
||||
namespace IS = CGAL::Isosurfacing;
|
||||
|
||||
template <typename PolygonMesh>
|
||||
std::size_t connected_components(PolygonMesh& mesh)
|
||||
{
|
||||
using face_descriptor = typename boost::graph_traits<PolygonMesh>::face_descriptor;
|
||||
|
||||
auto fccmap = mesh.template add_property_map<face_descriptor, std::size_t>("f:CC").first;
|
||||
|
||||
return CGAL::Polygon_mesh_processing::connected_components(mesh, fccmap);
|
||||
}
|
||||
|
||||
template <typename PolygonMesh>
|
||||
std::size_t euler_characteristic(const PolygonMesh& mesh)
|
||||
{
|
||||
return mesh.number_of_vertices() - mesh.number_of_edges() + mesh.number_of_faces();
|
||||
}
|
||||
|
||||
template <typename PolygonMesh>
|
||||
std::size_t boundary_components(const PolygonMesh& mesh)
|
||||
{
|
||||
using halfedge_descriptor = typename boost::graph_traits<PolygonMesh>::halfedge_descriptor;
|
||||
std::vector<halfedge_descriptor> border_cycles;
|
||||
|
||||
CGAL::Polygon_mesh_processing::extract_boundary_cycles(mesh, std::back_inserter(border_cycles));
|
||||
|
||||
return border_cycles.size();
|
||||
}
|
||||
|
||||
template <typename Grid>
|
||||
bool check_closed_not_empty(const Grid& grid, const IS::Interpolated_discrete_values_3<Grid>& values, const typename Grid::Geom_traits::FT iso = 0)
|
||||
{
|
||||
|
|
@ -441,29 +414,6 @@ FT generate_predefined_plane(const PlaneCase topology_case, std::array<FT, 8>& v
|
|||
return iso;
|
||||
}
|
||||
|
||||
template<typename Domain>
|
||||
void write_debug_grid(const Domain& domain, const std::string& filename) {
|
||||
using Point = typename Domain::Geom_traits::Point_3;
|
||||
using Mesh = CGAL::Surface_mesh<Point>;
|
||||
|
||||
Mesh debug_grid;
|
||||
auto debug_grid_creator = [&](const typename Domain::cell_descriptor& c)
|
||||
{
|
||||
std::vector<typename Mesh::Vertex_index> 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<CGAL::Sequential_tag>(debug_grid_creator);
|
||||
CGAL::IO::write_OFF(filename, debug_grid);
|
||||
}
|
||||
|
||||
template<typename KERNEL>
|
||||
void compare_tmc_mc_trilinear(const std::array<typename KERNEL::FT, 8>& case_values, typename KERNEL::FT iso)
|
||||
{
|
||||
|
|
|
|||
|
|
@ -7,6 +7,7 @@
|
|||
#include <CGAL/Polygon_mesh_processing/repair.h>
|
||||
#include <CGAL/Polygon_mesh_processing/repair_degeneracies.h>
|
||||
#include <CGAL/Polygon_mesh_processing/repair_polygon_soup.h>
|
||||
#include <CGAL/Surface_mesh.h>
|
||||
|
||||
#include <CGAL/Named_function_parameters.h>
|
||||
#include <CGAL/boost/graph/named_params_helper.h>
|
||||
|
|
@ -59,6 +60,48 @@ bool has_degenerate_faces(PolygonMesh& pmesh)
|
|||
return !dfs.empty();
|
||||
}
|
||||
|
||||
template <typename PolygonMesh>
|
||||
int connected_components(PolygonMesh& mesh)
|
||||
{
|
||||
using face_descriptor = typename boost::graph_traits<PolygonMesh>::face_descriptor;
|
||||
|
||||
auto fccmap = mesh.template add_property_map<face_descriptor, std::size_t>("f:CC").first;
|
||||
|
||||
return CGAL::Polygon_mesh_processing::connected_components(mesh, fccmap);
|
||||
}
|
||||
|
||||
template <typename PolygonMesh>
|
||||
int euler_characteristic(const PolygonMesh& mesh)
|
||||
{
|
||||
return mesh.number_of_vertices() - mesh.number_of_edges() + mesh.number_of_faces();
|
||||
}
|
||||
|
||||
template <typename PolygonMesh>
|
||||
int boundary_components(const PolygonMesh& mesh)
|
||||
{
|
||||
using halfedge_descriptor = typename boost::graph_traits<PolygonMesh>::halfedge_descriptor;
|
||||
std::vector<halfedge_descriptor> border_cycles;
|
||||
|
||||
CGAL::Polygon_mesh_processing::extract_boundary_cycles(mesh, std::back_inserter(border_cycles));
|
||||
|
||||
return border_cycles.size();
|
||||
}
|
||||
|
||||
template <typename PolygonMesh>
|
||||
int betti_0(PolygonMesh& mesh) {
|
||||
return connected_components(mesh);
|
||||
}
|
||||
|
||||
template <typename PolygonMesh>
|
||||
int betti_2(PolygonMesh& mesh) {
|
||||
return connected_components(mesh); // only for closed surfaces
|
||||
}
|
||||
|
||||
template <typename PolygonMesh>
|
||||
int betti_1(PolygonMesh& mesh) {
|
||||
return betti_0(mesh) + betti_2(mesh) - euler_characteristic(mesh);
|
||||
}
|
||||
|
||||
// template <typename PolygonMesh>
|
||||
// bool check_mesh_distance(const PolygonMesh& m0, const PolygonMesh& m1)
|
||||
// {
|
||||
|
|
@ -68,4 +111,27 @@ bool has_degenerate_faces(PolygonMesh& pmesh)
|
|||
// return true;
|
||||
// }
|
||||
|
||||
template<typename Domain>
|
||||
void write_debug_grid(const Domain& domain, const std::string& filename) {
|
||||
using Point = typename Domain::Geom_traits::Point_3;
|
||||
using Mesh = CGAL::Surface_mesh<Point>;
|
||||
|
||||
Mesh debug_grid;
|
||||
auto debug_grid_creator = [&](const typename Domain::cell_descriptor& c)
|
||||
{
|
||||
std::vector<typename Mesh::Vertex_index> 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<CGAL::Sequential_tag>(debug_grid_creator);
|
||||
CGAL::IO::write_OFF(filename, debug_grid);
|
||||
}
|
||||
|
||||
#endif // CGAL_ISOSURFACING_TEST_UTIL_H
|
||||
|
|
|
|||
|
|
@ -0,0 +1,238 @@
|
|||
|
||||
#include "test_util.h"
|
||||
|
||||
#include <CGAL/Simple_cartesian.h>
|
||||
#include <CGAL/Surface_mesh.h>
|
||||
|
||||
#include <CGAL/Isosurfacing_3/Cartesian_grid_3.h>
|
||||
#include <CGAL/Isosurfacing_3/marching_cubes_3.h>
|
||||
#include <CGAL/Isosurfacing_3/Marching_cubes_domain_3.h>
|
||||
#include <CGAL/Isosurfacing_3/Interpolated_discrete_values_3.h>
|
||||
#include <CGAL/Isosurfacing_3/Value_function_3.h>
|
||||
#include <CGAL/IO/polygon_soup_io.h>
|
||||
|
||||
/*
|
||||
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 <typename Grid, typename Values>
|
||||
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<char*>(&nx), sizeof(int));
|
||||
file.read(reinterpret_cast<char*>(&ny), sizeof(int));
|
||||
file.read(reinterpret_cast<char*>(&nz), sizeof(int));
|
||||
|
||||
// Read the bounding box
|
||||
// (xmin, xmax, ymin, ymax, zmin, zmax)
|
||||
std::array<float, 6> bbox;
|
||||
for (int i = 0; i < 6; ++i) {
|
||||
file.read(reinterpret_cast<char*>(&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<std::size_t>(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<float> volume_data(total_points);
|
||||
for (std::size_t i = 0; i < total_points; ++i) {
|
||||
file.read(reinterpret_cast<char*>(&volume_data[i]), sizeof(float));
|
||||
std::cout << volume_data[i] << std::endl;
|
||||
}
|
||||
|
||||
for(int x=0; x<nx; ++x)
|
||||
for(int y=0; y<ny; ++y)
|
||||
for(int z=0; z<nz; ++z)
|
||||
values(x, y, z) = volume_data[x + y * nx + z * nx * ny];
|
||||
|
||||
file.close();
|
||||
}
|
||||
|
||||
int read_euler(const std::string& filename) {
|
||||
std::ifstream file(filename);
|
||||
if (!file.is_open()) {
|
||||
throw std::runtime_error("Cannot open file: " + filename);
|
||||
}
|
||||
int euler;
|
||||
file >> euler;
|
||||
return euler;
|
||||
}
|
||||
|
||||
std::array<int, 2> 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 <typename K>
|
||||
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<Point>;
|
||||
using Grid = IS::Cartesian_grid_3<K>;
|
||||
using Values = IS::Interpolated_discrete_values_3<Grid>;
|
||||
using Domain = IS::Marching_cubes_domain_3<Grid, Values>;
|
||||
|
||||
using Point_range = std::vector<Point>;
|
||||
using Polygon_range = std::vector<std::vector<std::size_t> >;
|
||||
|
||||
using Mesh = CGAL::Surface_mesh<Point>;
|
||||
|
||||
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<CGAL::Sequential_tag>(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 <typename K>
|
||||
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<Point>;
|
||||
using Grid = IS::Cartesian_grid_3<K>;
|
||||
using Values = IS::Interpolated_discrete_values_3<Grid>;
|
||||
using Domain = IS::Marching_cubes_domain_3<Grid, Values>;
|
||||
|
||||
using Point_range = std::vector<Point>;
|
||||
using Polygon_range = std::vector<std::vector<std::size_t> >;
|
||||
|
||||
using Mesh = CGAL::Surface_mesh<Point>;
|
||||
|
||||
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<CGAL::Sequential_tag>(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 <typename K>
|
||||
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<Point>;
|
||||
using Grid = IS::Cartesian_grid_3<K>;
|
||||
using Values = IS::Interpolated_discrete_values_3<Grid>;
|
||||
using Domain = IS::Marching_cubes_domain_3<Grid, Values>;
|
||||
|
||||
using Point_range = std::vector<Point>;
|
||||
using Polygon_range = std::vector<std::vector<std::size_t> >;
|
||||
|
||||
using Mesh = CGAL::Surface_mesh<Point>;
|
||||
|
||||
Grid grid;
|
||||
Values values { grid };
|
||||
Domain domain { grid, values };
|
||||
|
||||
read_iso_volume(filename, grid, values);
|
||||
|
||||
Point_range points;
|
||||
Polygon_range triangles;
|
||||
IS::marching_cubes<CGAL::Sequential_tag>(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<std::size_t, 3>{151, 151, 151} };
|
||||
IS::Value_function_3<Grid> values_high_res { values, grid_high_res };
|
||||
IS::Marching_cubes_domain_3<Grid, IS::Value_function_3<Grid>> domain_high_res { grid_high_res, values_high_res };
|
||||
|
||||
Point_range points_high_res;
|
||||
Polygon_range triangles_high_res;
|
||||
IS::marching_cubes<CGAL::Parallel_if_available_tag>(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<double>;
|
||||
using FT = typename K::FT;
|
||||
|
||||
// verify_euler<K>();
|
||||
// verify_betti<K>();
|
||||
compare_to_reference<K>("data/MarchingCubes_cases/Grids/" + std::to_string(100) + "-scalar_field.iso");
|
||||
// compare_to_reference<K>("data/Closed_Surfaces/Grids/" + std::to_string(0) + "-scalar_field.iso");
|
||||
}
|
||||
Loading…
Reference in New Issue