From aa32e718769a5dd8814ed19083cdf00bd23778b4 Mon Sep 17 00:00:00 2001 From: Pierre Alliez Date: Sun, 24 Dec 2023 18:40:18 +0100 Subject: [PATCH] multiple offsets (for the teaser) --- .../examples/Isosurfacing_3/CMakeLists.txt | 3 + .../marching_cubes_multiple_mesh_offsets.cpp | 117 ++++++++++++++++++ 2 files changed, 120 insertions(+) create mode 100644 Isosurfacing_3/examples/Isosurfacing_3/marching_cubes_multiple_mesh_offsets.cpp diff --git a/Isosurfacing_3/examples/Isosurfacing_3/CMakeLists.txt b/Isosurfacing_3/examples/Isosurfacing_3/CMakeLists.txt index f87d4485e63..d32bd8a1962 100644 --- a/Isosurfacing_3/examples/Isosurfacing_3/CMakeLists.txt +++ b/Isosurfacing_3/examples/Isosurfacing_3/CMakeLists.txt @@ -9,8 +9,10 @@ find_package(CGAL REQUIRED) create_single_source_cgal_program( "marching_cubes_implicit_sphere.cpp" ) create_single_source_cgal_program( "marching_cubes_Cartesian_grid_sphere.cpp" ) create_single_source_cgal_program( "marching_cubes_signed_mesh_offset.cpp" ) +create_single_source_cgal_program( "marching_cubes_multiple_mesh_offsets.cpp" ) create_single_source_cgal_program( "marching_cubes_inrimage.cpp" ) + find_package(Eigen3 3.1.0 QUIET) #(3.1.0 or greater) include(CGAL_Eigen3_support) if(TARGET CGAL::Eigen3_support) @@ -28,6 +30,7 @@ if(TARGET CGAL::Eigen3_support) create_single_source_cgal_program( "dual_contouring_implicit_iwp.cpp" ) target_link_libraries(dual_contouring_implicit_iwp PRIVATE CGAL::Eigen3_support) + else() message(STATUS "NOTICE: Some examples use Eigen, and will not be compiled.") endif() diff --git a/Isosurfacing_3/examples/Isosurfacing_3/marching_cubes_multiple_mesh_offsets.cpp b/Isosurfacing_3/examples/Isosurfacing_3/marching_cubes_multiple_mesh_offsets.cpp new file mode 100644 index 00000000000..aeb76749520 --- /dev/null +++ b/Isosurfacing_3/examples/Isosurfacing_3/marching_cubes_multiple_mesh_offsets.cpp @@ -0,0 +1,117 @@ +#include +#include + +#include +#include +#include + +#include +#include +#include +#include + +#include + +#include +#include + +using Kernel = CGAL::Simple_cartesian; +using FT = typename Kernel::FT; +using Point = typename Kernel::Point_3; +using Vector = typename Kernel::Vector_3; + +using Grid = CGAL::Isosurfacing::Cartesian_grid_3; + +using Mesh = CGAL::Surface_mesh; + +using Primitive = CGAL::AABB_face_graph_triangle_primitive; +using Traits = CGAL::AABB_traits; +using Tree = CGAL::AABB_tree; + +using Point_range = std::vector; +using Polygon_range = std::vector >; + +// computes the Euclidean distance from query point p to the mesh +// via the AABB tree data structure +inline Kernel::FT distance_to_mesh(const Tree& tree, + const Point& p) +{ + const Point x = tree.closest_point(p); + return sqrt((p - x).squared_length()); +} + +int main(int argc, char **argv) +{ + const std::string input_name(argv[1]); + + std::cout << "Input file: " << input_name << std::endl; + + // load input mesh + Mesh mesh_input; + if(!CGAL::IO::read_OFF(input_name, mesh_input)) + { + std::cerr << "Could not read input mesh" << std::endl; + return EXIT_FAILURE; + } + + // construct loose bounding box from input mesh + CGAL::Bbox_3 aabb_grid = CGAL::Polygon_mesh_processing::bbox(mesh_input); + const FT loose_offset = 3.0; + Vector aabb_increase_vec = Vector(loose_offset, loose_offset, loose_offset); + aabb_grid += (Point(aabb_grid.xmax(), aabb_grid.ymax(), aabb_grid.zmax()) + aabb_increase_vec).bbox(); + aabb_grid += (Point(aabb_grid.xmin(), aabb_grid.ymin(), aabb_grid.zmin()) - aabb_increase_vec).bbox(); + + // construct AABB tree and functor to address inside/outside point queries + Tree tree(mesh_input.faces_begin(), mesh_input.faces_end(), mesh_input); + CGAL::Side_of_triangle_mesh::type> sotm(mesh_input); + + // create grid + const int n_voxels = 250; + Grid grid { n_voxels, n_voxels, n_voxels, aabb_grid }; + + for(std::size_t z = 0; z < grid.zdim(); ++z) { + for(std::size_t y = 0; y < grid.ydim(); ++y) { + for(std::size_t x = 0; x < grid.xdim(); ++x) + { + const FT pos_x = x * grid.spacing()[0] + grid.bbox().xmin(); + const FT pos_y = y * grid.spacing()[1] + grid.bbox().ymin(); + const FT pos_z = z * grid.spacing()[2] + grid.bbox().zmin(); + const Point p(pos_x, pos_y, pos_z); + + // compute unsigned distance to input mesh + grid.value(x, y, z) = distance_to_mesh(tree, p); + + // sign distance so that it is negative inside the mesh + const bool is_inside = (sotm(p) == CGAL::ON_BOUNDED_SIDE); + if(is_inside) + grid.value(x, y, z) *= -1.0; + } + } + } + + // create domain from the grid + auto domain = CGAL::Isosurfacing::create_explicit_Cartesian_grid_domain(grid); + + int index = 0; + for(FT offset = 0.0; offset < 0.3; offset += 0.01, index++) + { + // containers for the triangle soup output + Point_range points; + Polygon_range polygons; + + // execute marching cubes with an isovalue equating offset + std::cout << "Marching cubes with offset " << offset << "..."; + CGAL::Isosurfacing::marching_cubes(domain, offset, points, polygons); + std::cout << "done" << std::endl; + + // save the output + std::string filename("output-"); + filename.append(std::to_string(index)); + filename.append(std::string(".off")); + std::cout << "Save to file " << filename << "..."; + CGAL::IO::write_polygon_soup(filename, points, polygons); + std::cout << "done" << std::endl; + } + + return EXIT_SUCCESS; +}