diff --git a/Isosurfacing_3/examples/Isosurfacing_3/CMakeLists.txt b/Isosurfacing_3/examples/Isosurfacing_3/CMakeLists.txt index a3633c99adc..a07b86852e0 100644 --- a/Isosurfacing_3/examples/Isosurfacing_3/CMakeLists.txt +++ b/Isosurfacing_3/examples/Isosurfacing_3/CMakeLists.txt @@ -6,51 +6,51 @@ project( Isosurfacing_3_Examples ) find_package(CGAL REQUIRED) -create_single_source_cgal_program( "marching_cubes_implicit_sphere.cpp" ) -create_single_source_cgal_program( "marching_cubes_seq_vs_parallel.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) + +find_package(TBB QUIET) +include(CGAL_TBB_support) + +create_single_source_cgal_program("marching_cubes.cpp") + if(TARGET CGAL::Eigen3_support) - create_single_source_cgal_program( "dual_contouring_Cartesian_grid.cpp" ) - target_link_libraries(dual_contouring_Cartesian_grid PRIVATE CGAL::Eigen3_support) + create_single_source_cgal_program("dual_contouring.cpp") + create_single_source_cgal_program("contouring_discrete_data.cpp") + create_single_source_cgal_program("contouring_image.cpp") + create_single_source_cgal_program("contouring_implicit_data.cpp") + create_single_source_cgal_program("contouring_mesh_offset.cpp") + create_single_source_cgal_program("contouring_octree.cpp") - create_single_source_cgal_program( "dual_contouring_mesh_offset.cpp" ) - target_link_libraries(dual_contouring_mesh_offset PRIVATE CGAL::Eigen3_support) + # undocumented + create_single_source_cgal_program("dual_contouring_strategies.cpp") + create_single_source_cgal_program("dual_contouring_intersection_oracles.cpp") - create_single_source_cgal_program( "dual_contouring_octree.cpp" ) - target_link_libraries(dual_contouring_octree PRIVATE CGAL::Eigen3_support) + target_link_libraries(dual_contouring PRIVATE CGAL::Eigen3_support) + target_link_libraries(contouring_discrete_data PRIVATE CGAL::Eigen3_support) + target_link_libraries(contouring_image PRIVATE CGAL::Eigen3_support) + target_link_libraries(contouring_implicit_data PRIVATE CGAL::Eigen3_support) + target_link_libraries(contouring_mesh_offset PRIVATE CGAL::Eigen3_support) + target_link_libraries(contouring_octree PRIVATE CGAL::Eigen3_support) - create_single_source_cgal_program( "all_Cartesian_cube.cpp" ) - target_link_libraries(all_Cartesian_cube PRIVATE CGAL::Eigen3_support) + target_link_libraries(dual_contouring_strategies PRIVATE CGAL::Eigen3_support) + target_link_libraries(dual_contouring_intersection_oracles PRIVATE CGAL::Eigen3_support) - create_single_source_cgal_program( "dual_contouring_implicit_iwp.cpp" ) - target_link_libraries(dual_contouring_implicit_iwp PRIVATE CGAL::Eigen3_support) + if(TARGET CGAL::TBB_support) + target_link_libraries(dual_contouring PRIVATE CGAL::TBB_support) + target_link_libraries(contouring_discrete_data PRIVATE CGAL::TBB_support) + target_link_libraries(contouring_image PRIVATE CGAL::TBB_support) + target_link_libraries(contouring_implicit_data PRIVATE CGAL::TBB_support) + target_link_libraries(contouring_mesh_offset PRIVATE CGAL::TBB_support) + target_link_libraries(contouring_octree PRIVATE CGAL::TBB_support) + target_link_libraries(dual_contouring_strategies PRIVATE CGAL::TBB_support) + target_link_libraries(dual_contouring_intersection_oracles PRIVATE CGAL::TBB_support) + endif() else() message(STATUS "NOTICE: Some examples use Eigen, and will not be compiled.") endif() -find_package(TBB QUIET) -include(CGAL_TBB_support) if(TARGET CGAL::TBB_support) - target_link_libraries(marching_cubes_seq_vs_parallel PRIVATE CGAL::TBB_support) - target_link_libraries(marching_cubes_implicit_sphere PRIVATE CGAL::TBB_support) - target_link_libraries(marching_cubes_Cartesian_grid_sphere PRIVATE CGAL::TBB_support) - target_link_libraries(marching_cubes_signed_mesh_offset PRIVATE CGAL::TBB_support) - target_link_libraries(marching_cubes_multiple_mesh_offsets PRIVATE CGAL::TBB_support) - target_link_libraries(marching_cubes_inrimage PRIVATE CGAL::TBB_support) - - if(TARGET CGAL::Eigen3_support) - target_link_libraries(dual_contouring_Cartesian_grid PRIVATE CGAL::TBB_support) - target_link_libraries(dual_contouring_mesh_offset PRIVATE CGAL::TBB_support) - target_link_libraries(dual_contouring_octree PRIVATE CGAL::TBB_support) - target_link_libraries(all_Cartesian_cube PRIVATE CGAL::TBB_support) - target_link_libraries(dual_contouring_implicit_iwp PRIVATE CGAL::TBB_support) - endif() + target_link_libraries(marching_cubes PRIVATE CGAL::TBB_support) endif() diff --git a/Isosurfacing_3/examples/Isosurfacing_3/all_Cartesian_cube.cpp b/Isosurfacing_3/examples/Isosurfacing_3/all_Cartesian_cube.cpp deleted file mode 100644 index 721480f4a8b..00000000000 --- a/Isosurfacing_3/examples/Isosurfacing_3/all_Cartesian_cube.cpp +++ /dev/null @@ -1,87 +0,0 @@ -#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 Point_range = std::vector; -using Polygon_range = std::vector >; - -// return 1.0 if value has positive sign, and -1.0 otherwise -FT sign(FT value) -{ - return (value > 0.0) - (value < 0.0); -} - -int main(int, char**) -{ - // create a Cartesian grid with 7^3 grid points and the bounding box [-1, 1]^3 - const CGAL::Bbox_3 bbox{-1., -1., -1., 1., 1., 1.}; - Grid grid { 7, 7, 7, bbox }; - - // calculate the value at all grid points - for(std::size_t i=0; i 0.00001) - g /= CGAL::approximate_sqrt(length_sq); - - return g; - }; - - // create domain from given grid and gradient - auto domain = CGAL::Isosurfacing::create_explicit_Cartesian_grid_domain(grid, cube_gradient); - - // containers for output indexed surface meshes - Point_range points_mc, points_dc; - Polygon_range polygons_mc, polygons_dc; - - // run topologically correct Marching Cubes and Dual Contouring with given isovalue - const FT isovalue = 0.88; - CGAL::Isosurfacing::marching_cubes(domain, isovalue, points_mc, polygons_mc); - CGAL::Isosurfacing::dual_contouring(domain, isovalue, points_dc, polygons_dc); - - // save output indexed meshes to files, in the OFF format @fixme these are not meshes... - CGAL::IO::write_OFF("output_mc.off", points_mc, polygons_mc); - CGAL::IO::write_OFF("output_dc.off", points_dc, polygons_dc); - - return EXIT_SUCCESS; -} diff --git a/Isosurfacing_3/examples/Isosurfacing_3/contouring_discrete_data.cpp b/Isosurfacing_3/examples/Isosurfacing_3/contouring_discrete_data.cpp new file mode 100644 index 00000000000..12d282288f8 --- /dev/null +++ b/Isosurfacing_3/examples/Isosurfacing_3/contouring_discrete_data.cpp @@ -0,0 +1,124 @@ +#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 Values = CGAL::Isosurfacing::Interpolated_discrete_values_3; +using Gradients = CGAL::Isosurfacing::Interpolated_discrete_gradients_3; + +using Point_range = std::vector; +using Polygon_range = std::vector >; + +namespace IS = CGAL::Isosurfacing; + +void run_marching_cubes(const Grid& grid, + const FT isovalue) +{ + using Domain = IS::Marching_cubes_domain_3; + + std::cout << "\n ---- " << std::endl; + std::cout << "Running Marching Cubes with isovalue = " << isovalue << std::endl; + + // fill up values + Values values { grid }; + + for(std::size_t i=0; i; + + std::cout << "\n ---- " << std::endl; + std::cout << "Running Dual Contouring with isovalue = " << isovalue << std::endl; + + // fill up values and gradients + Values values { grid }; + Gradients gradients { grid }; + + for(std::size_t i=0; i 1) ? std::stod(argv[1]) : 0.8; + + // create bounding box and grid + const CGAL::Bbox_3 bbox { -1., -1., -1., 1., 1., 1. }; + Grid grid { bbox, 30, 30, 30 }; + + std::cout << "Bbox: " << grid.bbox() << std::endl; + std::cout << "Cell dimensions: " << grid.spacing()[0] << " " << grid.spacing()[1] << " " << grid.spacing()[2] << std::endl; + std::cout << "Cell #: " << grid.xdim() << ", " << grid.ydim() << ", " << grid.zdim() << std::endl; + + run_marching_cubes(grid, isovalue); + + run_dual_contouring(grid, isovalue); + + std::cout << "Done" << std::endl; + + return EXIT_SUCCESS; +} diff --git a/Isosurfacing_3/examples/Isosurfacing_3/contouring_image.cpp b/Isosurfacing_3/examples/Isosurfacing_3/contouring_image.cpp new file mode 100644 index 00000000000..18d02a792b7 --- /dev/null +++ b/Isosurfacing_3/examples/Isosurfacing_3/contouring_image.cpp @@ -0,0 +1,114 @@ +#include + +#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 Grid = CGAL::Isosurfacing::Cartesian_grid_3; +using Values = CGAL::Isosurfacing::Interpolated_discrete_values_3; +using Gradients = CGAL::Isosurfacing::Finite_difference_gradient_3; + +using Point_range = std::vector; +using Polygon_range = std::vector >; + +namespace IS = CGAL::Isosurfacing; + +void run_marching_cubes(const Grid& grid, + const FT isovalue, + const Values& values) +{ + using Domain = IS::Marching_cubes_domain_3; + + std::cout << "\n ---- " << std::endl; + std::cout << "Running Marching Cubes with isovalue = " << isovalue << std::endl; + + // fill up values + + // create a domain from the grid + Domain domain { grid, values }; + + // prepare collections for the output indexed soup + Point_range points; + Polygon_range triangles; + + // execute marching cubes + IS::marching_cubes(domain, isovalue, points, triangles); + + std::cout << "Output #vertices: " << points.size() << std::endl; + std::cout << "Output #triangles: " << triangles.size() << std::endl; + + // save output indexed mesh to a file, in the OFF format + CGAL::IO::write_polygon_soup("marching_cubes_image.off", points, triangles); +} + +void run_dual_contouring(const Grid& grid, + const FT isovalue, + const Values& values) +{ + using Domain = IS::Dual_contouring_domain_3; + + std::cout << "\n ---- " << std::endl; + std::cout << "Running Dual Contouring with isovalue = " << isovalue << std::endl; + + // fill up values and gradients + Gradients gradients { values }; + Domain domain { grid, values, gradients }; + + Point_range points; + Polygon_range triangles; + + // run dual contouring isosurfacing + IS::dual_contouring(domain, isovalue, points, triangles); + + std::cout << "Output #vertices: " << points.size() << std::endl; + std::cout << "Output #triangles: " << triangles.size() << std::endl; + CGAL::IO::write_polygon_soup("dual_contouring_image.off", points, triangles); +} + +int main(int argc, char* argv[]) +{ + const std::string fname = (argc > 1) ? argv[1] : CGAL::data_file_path("images/skull_2.9.inr"); + const FT isovalue = (argc > 2) ? std::stod(argv[2]) : - 2.9; + + // load volumetric image from a file + CGAL::Image_3 image; + if(!image.read(fname)) + { + std::cerr << "Error: Cannot read image file " << fname << std::endl; + return EXIT_FAILURE; + } + + // convert image to a Cartesian grid + auto [grid, values] = IS::IO::read_Image_3(image); + + for (std::size_t i=0; i + +#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 Values = CGAL::Isosurfacing::Value_function_3; +using Gradients = CGAL::Isosurfacing::Gradient_function_3; + +using Point_range = std::vector; +using Polygon_range = std::vector >; + +// --- +const FT alpha = 5.01; + +auto iwp_value = [](const Point& point) +{ + const FT x = alpha * (point.x() + FT(1.0)) * CGAL_PI; + const FT y = alpha * (point.y() + FT(1.0)) * CGAL_PI; + const FT z = alpha * (point.z() + FT(1.0)) * CGAL_PI; + return cos(x)*cos(y) + cos(y)*cos(z) + cos(z)*cos(x) - cos(x)*cos(y)*cos(z); // isovalue = 0 +}; +auto iwp_gradient = [](const Point& point) +{ + const FT x = alpha * (point.x() + FT(1.0)) * CGAL_PI; + const FT y = alpha * (point.y() + FT(1.0)) * CGAL_PI; + const FT z = alpha * (point.z() + FT(1.0)) * CGAL_PI; + + const FT gx = CGAL_PI * alpha * sin(x) * (cos(y) * (cos(z) - FT(1.0)) - cos(z)); + const FT gy = CGAL_PI * alpha * sin(y) * (cos(x) * (cos(z) - FT(1.0)) - cos(z)); + const FT gz = CGAL_PI * alpha * sin(z) * (cos(x) * (cos(y) - FT(1.0)) - cos(y)); + return Vector(gx, gy, gz); +}; + +void run_marching_cubes(const Grid& grid, + const FT isovalue) +{ + using Domain = CGAL::Isosurfacing::Marching_cubes_domain_3; + + std::cout << "\n ---- " << std::endl; + std::cout << "Running Marching Cubes with isovalue = " << isovalue << std::endl; + + // fill up values + Values values { iwp_value, grid }; + Domain domain { grid, values }; + + // output containers + Point_range points; + Polygon_range triangles; + + // run Marching Cubes + CGAL::Isosurfacing::marching_cubes(domain, isovalue, points, triangles); + + std::cout << "Output #vertices (MC): " << points.size() << std::endl; + std::cout << "Output #triangles (MC): " << triangles.size() << std::endl; + CGAL::IO::write_polygon_soup("marching_cubes_implicit.off", points, triangles); +} + +void run_dual_contouring(const Grid& grid, + const FT isovalue) +{ + using Domain = CGAL::Isosurfacing::Dual_contouring_domain_3; + + std::cout << "\n ---- " << std::endl; + std::cout << "Running Dual Contouring with isovalue = " << isovalue << std::endl; + + // fill up values and gradients + Values values { iwp_value, grid }; + Gradients gradients { iwp_gradient, grid }; + Domain domain { grid, values, gradients }; + + // output containers + Point_range points; + Polygon_range triangles; + + // run Dual Contouring + CGAL::Isosurfacing::dual_contouring(domain, isovalue, points, triangles); + + std::cout << "Output #vertices (DC): " << points.size() << std::endl; + std::cout << "Output #triangles (DC): " << triangles.size() << std::endl; + CGAL::IO::write_polygon_soup("dual_contouring_implicit.off", points, triangles); +} + +int main(int argc, char** argv) +{ + const FT isovalue = (argc > 1) ? std::stod(argv[1]) : 0.; + + const CGAL::Bbox_3 bbox{-1, -1, -1, 1, 1, 1}; + const FT step = 0.0078125; // 0.02 + const std::array spacing { step, step, step }; + const Grid grid { bbox, spacing }; + + std::cout << "Bbox: " << grid.bbox() << std::endl; + std::cout << "Cell dimensions: " << grid.spacing()[0] << " " << grid.spacing()[1] << " " << grid.spacing()[2] << std::endl; + std::cout << "Cell #: " << grid.xdim() << ", " << grid.ydim() << ", " << grid.zdim() << std::endl; + + run_marching_cubes(grid, isovalue); + + run_dual_contouring(grid, isovalue); + + std::cout << "Done" << std::endl; + + return EXIT_SUCCESS; +} diff --git a/Isosurfacing_3/examples/Isosurfacing_3/contouring_mesh_offset.cpp b/Isosurfacing_3/examples/Isosurfacing_3/contouring_mesh_offset.cpp new file mode 100644 index 00000000000..4c6aed8bbb8 --- /dev/null +++ b/Isosurfacing_3/examples/Isosurfacing_3/contouring_mesh_offset.cpp @@ -0,0 +1,171 @@ +#include +#include + +#include +#include +#include +#include +#include +#include +#include + +#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 Mesh = CGAL::Surface_mesh; + +using Grid = CGAL::Isosurfacing::Cartesian_grid_3; +using Values = CGAL::Isosurfacing::Value_function_3; +using Gradients = CGAL::Isosurfacing::Finite_difference_gradient_3; + +using Point_range = std::vector; +using Polygon_range = std::vector >; + +struct Offset_oracle +{ + using Primitive = CGAL::AABB_face_graph_triangle_primitive; + using Traits = CGAL::AABB_traits; + using Tree = CGAL::AABB_tree; + +private: + const bool is_closed; + const Tree tree; + CGAL::Side_of_triangle_mesh sotm; + +public: + Offset_oracle(const Mesh& mesh) + : is_closed(CGAL::is_closed(mesh)), tree(mesh.faces_begin(), mesh.faces_end(), mesh), sotm(mesh) + { } + + FT distance(const Point& p) const + { + const Point cp = tree.closest_point(p); + FT d = sqrt((p - cp).squared_length()); + + if(is_closed && sotm(p) == (CGAL::ON_BOUNDED_SIDE)) + d *= -1; + + return d; + } +}; + +void run_marching_cubes(const Grid& grid, + const FT offset_value, + const Offset_oracle& offset_oracle) +{ + using Domain = CGAL::Isosurfacing::Marching_cubes_domain_3; + + std::cout << "\n ---- " << std::endl; + std::cout << "Running Marching Cubes with offset value = " << offset_value << std::endl; + + // fill up values + auto mesh_distance = [&offset_oracle](const Point& p) { return offset_oracle.distance(p); }; + Values values { mesh_distance, grid }; + Domain domain { grid, values }; + + Point_range points; + Polygon_range triangles; + + std::cout << "Output #vertices (MC): " << points.size() << std::endl; + std::cout << "Output #triangles (MC): " << triangles.size() << std::endl; + CGAL::IO::write_polygon_soup("marching_cubes_offsets.off", points, triangles); + +} + +void run_dual_contouring(const Grid& grid, + const FT offset_value, + const Offset_oracle& offset_oracle) +{ + using Domain = CGAL::Isosurfacing::Dual_contouring_domain_3; + + std::cout << "\n ---- " << std::endl; + std::cout << "Running Dual Contouring with offset value = " << offset_value << std::endl; + + // fill up values and gradients + auto mesh_distance = [&offset_oracle](const Point& p) { return offset_oracle.distance(p); }; + + Values values { mesh_distance, grid }; + Gradients gradients { values }; + + Domain domain { grid, values, gradients }; + + // output containers + Point_range points; + Polygon_range triangles; + + // run dual contouring + std::cout << "Running Dual Contouring with isovalue = " << offset_value << std::endl; + CGAL::Isosurfacing::dual_contouring(domain, offset_value, points, triangles); + + std::cout << "Output #vertices (DC): " << points.size() << std::endl; + std::cout << "Output #triangles (DC): " << triangles.size() << std::endl; + CGAL::IO::write_polygon_soup("dual_contouring_mesh_offset.off", points, triangles); +} + +int main(int argc, char** argv) +{ + const std::string filename = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/cross.off"); + const FT offset_value = (argc > 2) ? std::stod(argv[2]) : 0.2; + + if(offset_value <= 0) + { + std::cerr << "Offset value must be positive" << std::endl; + return EXIT_FAILURE; + } + + Mesh mesh; + if(!CGAL::IO::read_polygon_mesh(filename, mesh)) + { + std::cerr << "Could not read input mesh" << std::endl; + return EXIT_FAILURE; + } + + if(CGAL::is_closed(mesh)) + std::cout << "Input mesh is closed - using signed distance offset" << std::endl; + else + std::cout << "Input mesh is not closed - using unsigned distance offset" << std::endl; + + // construct loose bounding box from input mesh + CGAL::Bbox_3 bbox = CGAL::Polygon_mesh_processing::bbox(mesh); + + const FT diag_length = sqrt(CGAL::square(bbox.xmax() - bbox.xmin()) + + CGAL::square(bbox.ymax() - bbox.ymin()) + + CGAL::square(bbox.zmax() - bbox.zmin())); + const FT loose_offset = offset_value + 0.1 * diag_length; + + Vector aabb_increase_vec = Vector(loose_offset, loose_offset, loose_offset); + bbox += (Point(bbox.xmax(), bbox.ymax(), bbox.zmax()) + aabb_increase_vec).bbox(); + bbox += (Point(bbox.xmin(), bbox.ymin(), bbox.zmin()) - aabb_increase_vec).bbox(); + + const int n_voxels = 250; + Grid grid { bbox, n_voxels, n_voxels, n_voxels }; + + std::cout << "Bbox: " << grid.bbox() << std::endl; + std::cout << "Cell dimensions: " << grid.spacing()[0] << " " << grid.spacing()[1] << " " << grid.spacing()[2] << std::endl; + std::cout << "Cell #: " << grid.xdim() << ", " << grid.ydim() << ", " << grid.zdim() << std::endl; + + Offset_oracle offset_oracle(mesh); + + run_marching_cubes(grid, offset_value, offset_oracle); + + run_dual_contouring(grid, offset_value, offset_oracle); + + return EXIT_SUCCESS; +} diff --git a/Isosurfacing_3/examples/Isosurfacing_3/dual_contouring_octree.cpp b/Isosurfacing_3/examples/Isosurfacing_3/contouring_octree.cpp similarity index 81% rename from Isosurfacing_3/examples/Isosurfacing_3/dual_contouring_octree.cpp rename to Isosurfacing_3/examples/Isosurfacing_3/contouring_octree.cpp index 5bd54657d86..fd8b4774129 100644 --- a/Isosurfacing_3/examples/Isosurfacing_3/dual_contouring_octree.cpp +++ b/Isosurfacing_3/examples/Isosurfacing_3/contouring_octree.cpp @@ -1,10 +1,11 @@ #include #include -#include +#include +#include #include -#include +#include #include #include @@ -19,6 +20,9 @@ using Point_range = std::vector; using Polygon_range = std::vector >; using Octree_wrapper = CGAL::Isosurfacing::internal::Octree_wrapper; +using Values = CGAL::Isosurfacing::Value_function_3; +using Gradients = CGAL::Isosurfacing::Gradient_function_3; +using Domain = CGAL::Isosurfacing::Isosurfacing_domain_3; struct Refine_one_eighth { @@ -27,6 +31,14 @@ struct Refine_one_eighth std::size_t octree_dim_; + Refine_one_eighth(std::size_t min_depth, + std::size_t max_depth) + : min_depth_(min_depth), + max_depth_(max_depth) + { + octree_dim_ = std::size_t(1) << max_depth_; + } + Octree_wrapper::Uniform_coords uniform_coordinates(const Octree_wrapper::Octree::Node& node) const { auto coords = node.global_coordinates(); @@ -37,14 +49,6 @@ struct Refine_one_eighth return coords; } - Refine_one_eighth(std::size_t min_depth, - std::size_t max_depth) - : min_depth_(min_depth), - max_depth_(max_depth) - { - octree_dim_ = std::size_t(1) << max_depth_; - } - bool operator()(const Octree_wrapper::Octree::Node& n) const { if(n.depth() < min_depth_) @@ -87,16 +91,15 @@ int main(int, char**) return g / std::sqrt(g.squared_length()); }; - auto domain = CGAL::Isosurfacing::create_implicit_octree_domain(octree_wrap, - sphere_function, - sphere_gradient); + Data data(sphere_function, sphere_gradient); + Domain domain(octree_wrap, data); Point_range points; Polygon_range polygons; CGAL::Isosurfacing::dual_contouring(domain, 0.8, points, polygons); - CGAL::IO::write_OFF("dual_contouring_octree.off", points, polygons); + CGAL::IO::write_polygon_soup("dual_contouring_octree.off", points, polygons); return EXIT_SUCCESS; } diff --git a/Isosurfacing_3/examples/Isosurfacing_3/dual_contouring.cpp b/Isosurfacing_3/examples/Isosurfacing_3/dual_contouring.cpp new file mode 100644 index 00000000000..3823714069e --- /dev/null +++ b/Isosurfacing_3/examples/Isosurfacing_3/dual_contouring.cpp @@ -0,0 +1,73 @@ +#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 Values = CGAL::Isosurfacing::Value_function_3; +using Gradients = CGAL::Isosurfacing::Gradient_function_3; +using Domain = CGAL::Isosurfacing::Dual_contouring_domain_3; + +using Point_range = std::vector; +using Polygon_range = std::vector >; + +int main(int argc, char** argv) +{ + const FT isovalue = (argc > 1) ? std::stod(argv[1]) : 0.8; + + // create bounding box and grid + const CGAL::Bbox_3 bbox { -1., -1., -1., 1., 1., 1. }; + Grid grid { bbox, 30, 30, 30 }; + + std::cout << "Bbox: " << grid.bbox() << std::endl; + std::cout << "Cell dimensions: " << grid.spacing()[0] << " " << grid.spacing()[1] << " " << grid.spacing()[2] << std::endl; + std::cout << "Cell #: " << grid.xdim() << ", " << grid.ydim() << ", " << grid.zdim() << std::endl; + + // fill up values and gradients + auto sphere_value_fn = [](const Point& p) -> FT + { + return sqrt(p.x()*p.x() + p.y()*p.y() + p.z()*p.z()); + }; + + auto sphere_gradient_fn = [](const Point& p) -> Vector + { + const FT d = sqrt(p.x()*p.x() + p.y()*p.y() + p.z()*p.z()); + return Vector(CGAL::ORIGIN, p) / d; + }; + + Values values { sphere_value_fn, grid }; + Gradients gradients { sphere_gradient_fn, grid }; + + Domain domain { grid, values, gradients }; + // the domain could also be created with: + // auto domain = CGAL::Isosurfacing::create_dual_contouring_domain_3(grid, values, gradients); + + Point_range points; + Polygon_range triangles; + + // run dual contouring isosurfacing + std::cout << "Running Dual Contouring with isovalue = " << isovalue << std::endl; + CGAL::Isosurfacing::dual_contouring(domain, isovalue, points, triangles, + CGAL::parameters::do_not_triangulate_faces(true) + .constrain_to_cell(true)); + + std::cout << "Output #vertices: " << points.size() << std::endl; + std::cout << "Output #triangles: " << triangles.size() << std::endl; + CGAL::IO::write_polygon_soup("dual_contouring.off", points, triangles); + + std::cout << "Done" << std::endl; + + return EXIT_SUCCESS; +} diff --git a/Isosurfacing_3/examples/Isosurfacing_3/dual_contouring_Cartesian_grid.cpp b/Isosurfacing_3/examples/Isosurfacing_3/dual_contouring_Cartesian_grid.cpp deleted file mode 100644 index 22830782e67..00000000000 --- a/Isosurfacing_3/examples/Isosurfacing_3/dual_contouring_Cartesian_grid.cpp +++ /dev/null @@ -1,64 +0,0 @@ -#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 Point_range = std::vector; -using Polygon_range = std::vector >; - -int main(int, char**) -{ - // create bounding box and grid - const CGAL::Bbox_3 bbox{-1., -1., -1., 1., 1., 1.}; - Grid grid { 30, 30, 30, bbox }; - - // compute field values and gradients - 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 FT pos_x = x * grid.spacing()[0] + bbox.xmin(); - const FT pos_y = y * grid.spacing()[1] + bbox.ymin(); - const FT pos_z = z * grid.spacing()[2] + bbox.zmin(); - - const Vector direction(pos_x, pos_y, pos_z); - const FT distance = CGAL::approximate_sqrt(direction.squared_length()); - - grid.value(x, y, z) = distance; - - if(distance != 0) - grid.gradient(x, y, z) = direction / distance; - else - grid.gradient(x, y, z) = CGAL::NULL_VECTOR; - } - } - } - - // gradient field - CGAL::Isosurfacing::Explicit_Cartesian_grid_gradient_3 gradient(grid); - - // create domain from scalar and gradient fields - auto domain = CGAL::Isosurfacing::create_explicit_Cartesian_grid_domain(grid, gradient); - - Point_range points; - Polygon_range polygons; - - // run dual contouring isosurfacing - CGAL::Isosurfacing::dual_contouring(domain, 0.8, points, polygons); - - // write output indexed surface mesh to file, in OFF format - CGAL::IO::write_OFF("dual_contouring_Cartesian_grid.off", points, polygons); - - return EXIT_SUCCESS; -} diff --git a/Isosurfacing_3/examples/Isosurfacing_3/dual_contouring_implicit_iwp.cpp b/Isosurfacing_3/examples/Isosurfacing_3/dual_contouring_implicit_iwp.cpp deleted file mode 100644 index 8f58a1baeb3..00000000000 --- a/Isosurfacing_3/examples/Isosurfacing_3/dual_contouring_implicit_iwp.cpp +++ /dev/null @@ -1,59 +0,0 @@ -#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 Point_range = std::vector; -using Polygon_range = std::vector >; - -int main(int, char**) -{ - const FT alpha = 5.01; - - auto iwp_value = [alpha](const Point& point) - { - const FT x = alpha * (point.x() + FT(1.0)) * CGAL_PI; - const FT y = alpha * (point.y() + FT(1.0)) * CGAL_PI; - const FT z = alpha * (point.z() + FT(1.0)) * CGAL_PI; - return cos(x) * cos(y) + cos(y) * cos(z) + cos(z) * cos(x) - cos(x) * cos(y) * cos(z); // isovalue = 0 - }; - - auto iwp_gradient = [alpha](const Point& point) - { - const FT x = alpha * (point.x() + FT(1.0)) * CGAL_PI; - const FT y = alpha * (point.y() + FT(1.0)) * CGAL_PI; - const FT z = alpha * (point.z() + FT(1.0)) * CGAL_PI; - - const FT gx = CGAL_PI * alpha * sin(x) * (cos(y) * (cos(z) - FT(1.0)) - cos(z)); - const FT gy = CGAL_PI * alpha * sin(y) * (cos(x) * (cos(z) - FT(1.0)) - cos(z)); - const FT gz = CGAL_PI * alpha * sin(z) * (cos(x) * (cos(y) - FT(1.0)) - cos(y)); - return Vector(gx, gy, gz); - }; - - const CGAL::Bbox_3 bbox{-1.0, -1.0, -1.0, 1.0, 1.0, 1.0}; - const FT spacing = 0.02; - const Vector vec_spacing(spacing, spacing, spacing); - - // create a domain with given bounding box and grid spacing - auto domain = CGAL::Isosurfacing::create_implicit_Cartesian_grid_domain(bbox, vec_spacing, iwp_value, iwp_gradient); - - // prepare collections for the result - Point_range points; - Polygon_range polygons; - - // run dual contouring with isovalue set to 0 - CGAL::Isosurfacing::dual_contouring(domain, 0., points, polygons); - - // save output to the OFF format - CGAL::IO::write_OFF("dual_contouring_implicit_iwp.off", points, polygons); - - return EXIT_SUCCESS; -} diff --git a/Isosurfacing_3/examples/Isosurfacing_3/dual_contouring_intersection_oracles.cpp b/Isosurfacing_3/examples/Isosurfacing_3/dual_contouring_intersection_oracles.cpp new file mode 100644 index 00000000000..68e56894907 --- /dev/null +++ b/Isosurfacing_3/examples/Isosurfacing_3/dual_contouring_intersection_oracles.cpp @@ -0,0 +1,111 @@ +#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 Values = CGAL::Isosurfacing::Value_function_3; +using Gradients = CGAL::Isosurfacing::Finite_difference_gradient_3; + +using Point_range = std::vector; +using Polygon_range = std::vector >; + +namespace IS = CGAL::Isosurfacing; + +auto implicit_function = [](const Point& q) -> FT +{ + 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), 0.75, q); }; + return IS::Shapes::shape_difference(cyl_and_cube, sphere, q); +}; + +// The example shows that this is a bad choice to use a linear interpolant +// because the implicit function gives much more information than just using the values at grid vertices. +// @todo implement the SDF interpolant and show that it's as good as the dichotomy, but faster +int main(int argc, char** argv) +{ + const FT isovalue = (argc > 1) ? std::stod(argv[1]) : 0.; + std::cout << "Isovalue: " << isovalue << std::endl; + + // create bounding box and grid + const CGAL::Bbox_3 bbox = {-2., -2., -2., 2., 2., 2.}; + Grid grid { bbox, 50, 50, 50 }; + + std::cout << "Bbox: " << grid.bbox() << std::endl; + std::cout << "Cell dimensions: " << grid.spacing()[0] << " " << grid.spacing()[1] << " " << grid.spacing()[2] << std::endl; + std::cout << "Cell #: " << grid.xdim() << ", " << grid.ydim() << ", " << grid.zdim() << std::endl; + + // fill up values and gradients + Values values { implicit_function, grid }; + Gradients gradients { values }; + + const bool triangulate_faces = false; + + // Compute edge intersections with dichotomy + { + using Domain = IS::Dual_contouring_domain_3; + Domain domain { grid, values, gradients }; + + CGAL::Real_timer timer; + timer.start(); + + Point_range points; + Polygon_range triangles; + + std::cout << "--- Dual Contouring (Dichotomy)" << std::endl; + IS::dual_contouring(domain, isovalue, points, triangles, + CGAL::parameters::do_not_triangulate_faces(!triangulate_faces)); + + timer.stop(); + + std::cout << "Output #vertices: " << points.size() << std::endl; + std::cout << "Output #triangles: " << triangles.size() << std::endl; + std::cout << "Elapsed time: " << timer.time() << " seconds" << std::endl; + CGAL::IO::write_polygon_soup("dual_contouring-dichotomy_oracle.off", points, triangles); + } + + // Compute edge intersections with linear interpolation + { + using Domain = IS::Dual_contouring_domain_3; + Domain domain { grid, values, gradients }; + + CGAL::Real_timer timer; + timer.start(); + + Point_range points; + Polygon_range triangles; + + std::cout << "--- Dual Contouring (Linear Interpolation)" << std::endl; + IS::dual_contouring(domain, isovalue, points, triangles, + CGAL::parameters::do_not_triangulate_faces(!triangulate_faces)); + + timer.stop(); + + std::cout << "Output #vertices: " << points.size() << std::endl; + std::cout << "Output #triangles: " << triangles.size() << std::endl; + std::cout << "Elapsed time: " << timer.time() << " seconds" << std::endl; + CGAL::IO::write_polygon_soup("dual_contouring-linear_oracle.off", points, triangles); + } + + std::cout << "Done" << std::endl; + + return EXIT_SUCCESS; +} diff --git a/Isosurfacing_3/examples/Isosurfacing_3/dual_contouring_mesh_offset.cpp b/Isosurfacing_3/examples/Isosurfacing_3/dual_contouring_mesh_offset.cpp deleted file mode 100644 index 7f2fd02d70a..00000000000 --- a/Isosurfacing_3/examples/Isosurfacing_3/dual_contouring_mesh_offset.cpp +++ /dev/null @@ -1,87 +0,0 @@ -#include -#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 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 >; - -int main(int argc, char* argv[]) -{ - - const std::string input_name = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/cross.off"); - const Vector grid_spacing(0.1, 0.1, 0.1); - const FT offset_value = 0.2; - CGAL_assertion(offset_value > 0); - - Mesh input_mesh; - if(!CGAL::IO::read_OFF(input_name, input_mesh)) - { - std::cerr << "Could not read input mesh" << std::endl; - return EXIT_FAILURE; - } - - // compute bounding box - CGAL::Bbox_3 bbox = CGAL::Polygon_mesh_processing::bbox(input_mesh); - Vector aabb_increase_vec = Vector(offset_value + 0.01, offset_value + 0.01, offset_value + 0.01); - bbox += (Point(bbox.xmax(), bbox.ymax(), bbox.zmax()) + aabb_increase_vec).bbox(); - bbox += (Point(bbox.xmin(), bbox.ymin(), bbox.zmin()) - aabb_increase_vec).bbox(); - - // construct AABB tree - Tree tree(input_mesh.faces_begin(), input_mesh.faces_end(), input_mesh); - - CGAL::Side_of_triangle_mesh::type> sotm(input_mesh); - - // functors for addressing distance and normal queries - auto mesh_distance = [&tree](const Point& p) - { - const Point x = tree.closest_point(p); - return sqrt((p - x).squared_length()); - }; - - auto mesh_normal = [&tree](const Point& p) - { - const Point x = tree.closest_point(p); - const Vector n = p - x; - return n / sqrt(n.squared_length()); // normalize output vector - }; - - // create a domain with given bounding box and grid spacing - auto domain = CGAL::Isosurfacing::create_implicit_Cartesian_grid_domain(bbox, grid_spacing, - mesh_distance, mesh_normal); - // containers for output indexed surface mesh - Point_range points; - Polygon_range polygons; - - // run dual contouring - CGAL::Isosurfacing::dual_contouring(domain, offset_value, points, polygons); - - // save output indexed mesh to a file, in the OFF format - CGAL::IO::write_OFF("dual_contouring_mesh_offset.off", points, polygons); - - return EXIT_SUCCESS; -} diff --git a/Isosurfacing_3/examples/Isosurfacing_3/dual_contouring_strategies.cpp b/Isosurfacing_3/examples/Isosurfacing_3/dual_contouring_strategies.cpp new file mode 100644 index 00000000000..4c8a355ea20 --- /dev/null +++ b/Isosurfacing_3/examples/Isosurfacing_3/dual_contouring_strategies.cpp @@ -0,0 +1,154 @@ +#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 Values = CGAL::Isosurfacing::Value_function_3; +using Gradients = CGAL::Isosurfacing::Finite_difference_gradient_3; +using Domain = CGAL::Isosurfacing::Dual_contouring_domain_3; + +using Point_range = std::vector; +using Polygon_range = std::vector >; + +namespace IS = CGAL::Isosurfacing; + +auto implicit_function = [](const Point& q) -> FT +{ + 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), 0.75, q); }; + return IS::Shapes::shape_difference(cyl_and_cube, sphere, q); +}; + +int main(int argc, char** argv) +{ + const FT isovalue = (argc > 1) ? std::stod(argv[1]) : 0.; + std::cout << "Isovalue: " << isovalue << std::endl; + + // create bounding box and grid + const CGAL::Bbox_3 bbox = {-2., -2., -2., 2., 2., 2.}; + Grid grid { bbox, 50, 50, 50 }; + + std::cout << "Bbox: " << grid.bbox() << std::endl; + std::cout << "Cell dimensions: " << grid.spacing()[0] << " " << grid.spacing()[1] << " " << grid.spacing()[2] << std::endl; + std::cout << "Cell #: " << grid.xdim() << ", " << grid.ydim() << ", " << grid.zdim() << std::endl; + + // fill up values and gradients + Values values { implicit_function, grid }; + Gradients gradients { values }; + Domain domain { grid, values, gradients }; + + const bool triangulate_faces = false; + + // unconstrained QEM Placement strategy (default) + { + CGAL::Real_timer timer; + timer.start(); + + Point_range points; + Polygon_range triangles; + + std::cout << "--- Dual Contouring (QEM - unconstrained)" << std::endl; + IS::internal::Dual_contourer contourer; + contourer(domain, isovalue, points, triangles, + CGAL::parameters::do_not_triangulate_faces(!triangulate_faces)); + + timer.stop(); + + std::cout << "Output #vertices: " << points.size() << std::endl; + std::cout << "Output #triangles: " << triangles.size() << std::endl; + std::cout << "Elapsed time: " << timer.time() << " seconds" << std::endl; + CGAL::IO::write_polygon_soup("dual_contouring_QEM-unconstrained.off", points, triangles); + } + + // constrained QEM Placement strategy + { + CGAL::Real_timer timer; + timer.start(); + + Point_range points; + Polygon_range triangles; + + std::cout << "--- Dual Contouring (QEM - constrained)" << std::endl; + IS::internal::Dual_contourer contourer; + contourer(domain, isovalue, points, triangles, + CGAL::parameters::do_not_triangulate_faces(!triangulate_faces) + .constrain_to_cell(true)); + + timer.stop(); + + std::cout << "Output #vertices: " << points.size() << std::endl; + std::cout << "Output #triangles: " << triangles.size() << std::endl; + std::cout << "Elapsed time: " << timer.time() << " seconds" << std::endl; + CGAL::IO::write_polygon_soup("dual_contouring_QEM-constrained.off", points, triangles); + } + + + // Centroid of Edge Intersections strategy + { + CGAL::Real_timer timer; + timer.start(); + + Point_range points; + Polygon_range triangles; + + std::cout << "--- Dual Contouring (Centroid of Edge Intersections)" << std::endl; + IS::internal::Dual_contourer contourer; + contourer(domain, isovalue, points, triangles, + CGAL::parameters::do_not_triangulate_faces(!triangulate_faces)); + + timer.stop(); + + std::cout << "Output #vertices: " << points.size() << std::endl; + std::cout << "Output #triangles: " << triangles.size() << std::endl; + std::cout << "Elapsed time: " << timer.time() << " seconds" << std::endl; + CGAL::IO::write_polygon_soup("dual_contouring_CEI.off", points, triangles); + } + + // Cell Center strategy + { + CGAL::Real_timer timer; + timer.start(); + + Point_range points; + Polygon_range triangles; + + std::cout << "--- Dual Contouring (Cell Center)" << std::endl; + IS::internal::Dual_contourer contourer; + contourer(domain, isovalue, points, triangles, + CGAL::parameters::do_not_triangulate_faces(!triangulate_faces)); + + timer.stop(); + + std::cout << "Output #vertices: " << points.size() << std::endl; + std::cout << "Output #triangles: " << triangles.size() << std::endl; + std::cout << "Elapsed time: " << timer.time() << " seconds" << std::endl; + CGAL::IO::write_polygon_soup("dual_contouring_CC.off", points, triangles); + } + + std::cout << "Done" << std::endl; + + return EXIT_SUCCESS; +} diff --git a/Isosurfacing_3/examples/Isosurfacing_3/marching_cubes.cpp b/Isosurfacing_3/examples/Isosurfacing_3/marching_cubes.cpp new file mode 100644 index 00000000000..a5e3fb38bb9 --- /dev/null +++ b/Isosurfacing_3/examples/Isosurfacing_3/marching_cubes.cpp @@ -0,0 +1,62 @@ +#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 Values = CGAL::Isosurfacing::Value_function_3; +using Domain = CGAL::Isosurfacing::Marching_cubes_domain_3; + +using Point_range = std::vector; +using Polygon_range = std::vector >; + +int main(int argc, char** argv) +{ + const FT isovalue = (argc > 1) ? std::stod(argv[1]) : 0.8; + + // create bounding box and grid + const CGAL::Bbox_3 bbox { -1., -1., -1., 1., 1., 1. }; + Grid grid { bbox, 30, 30, 30 }; + + std::cout << "Bbox: " << grid.bbox() << std::endl; + std::cout << "Cell dimensions: " << grid.spacing()[0] << " " << grid.spacing()[1] << " " << grid.spacing()[2] << std::endl; + std::cout << "Cell #: " << grid.xdim() << ", " << grid.ydim() << ", " << grid.zdim() << std::endl; + + // fill up values + auto sphere_value_fn = [](const Point& p) -> FT + { + return sqrt(p.x()*p.x() + p.y()*p.y() + p.z()*p.z()); + }; + + Values values { sphere_value_fn, grid }; + Domain domain { grid, values }; + // the domain could also be created with: + // auto domain = CGAL::Isosurfacing::create_marching_cubes_domain_3(grid, values); + + Point_range points; + Polygon_range triangles; + + // run marching cubes isosurfacing + std::cout << "Running Marching Cubes with isovalue = " << isovalue << std::endl; + CGAL::Isosurfacing::marching_cubes(domain, isovalue, points, triangles, + CGAL::parameters::use_topologically_correct_marching_cubes(true)); + + std::cout << "Output #vertices: " << points.size() << std::endl; + std::cout << "Output #triangles: " << triangles.size() << std::endl; + CGAL::IO::write_polygon_soup("marching_cubes.off", points, triangles); + + std::cout << "Done" << std::endl; + + return EXIT_SUCCESS; +} diff --git a/Isosurfacing_3/examples/Isosurfacing_3/marching_cubes_Cartesian_grid_sphere.cpp b/Isosurfacing_3/examples/Isosurfacing_3/marching_cubes_Cartesian_grid_sphere.cpp deleted file mode 100644 index 51ac920815e..00000000000 --- a/Isosurfacing_3/examples/Isosurfacing_3/marching_cubes_Cartesian_grid_sphere.cpp +++ /dev/null @@ -1,54 +0,0 @@ -#include - -#include -#include -#include - -#include - -#include - -using Kernel = CGAL::Simple_cartesian; -using FT = typename Kernel::FT; -using Point = typename Kernel::Point_3; - -using Grid = CGAL::Isosurfacing::Cartesian_grid_3; -using Point_range = std::vector; -using Polygon_range = std::vector >; - -int main(int, char**) -{ - // 3D bounding box [-1, 1]^3 and Cartesian grid - const CGAL::Bbox_3 bbox{-1., -1., -1., 1., 1., 1.}; - Grid grid { 50, 50, 50, bbox }; - - // compute and store function values at all grid points - 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 FT pos_x = x * grid.spacing()[0] + bbox.xmin(); - const FT pos_y = y * grid.spacing()[1] + bbox.ymin(); - const FT pos_z = z * grid.spacing()[2] + bbox.zmin(); - - // Euclidean distance to the origin - grid.value(x, y, z) = sqrt(pos_x * pos_x + pos_y * pos_y + pos_z * pos_z); - } - } - } - - // create a domain from the grid - auto domain = CGAL::Isosurfacing::create_explicit_Cartesian_grid_domain(grid); - - // prepare collections for the result - Point_range points; - Polygon_range triangles; - - // run marching cubes with an isovalue of 0.8 - CGAL::Isosurfacing::marching_cubes(domain, 0.8, points, triangles); - - // save output indexed surface mesh to file, in the OFF format - CGAL::IO::write_OFF("marching_cubes_Cartesian_grid_sphere.off", points, triangles); - - return EXIT_SUCCESS; -} diff --git a/Isosurfacing_3/examples/Isosurfacing_3/marching_cubes_implicit_sphere.cpp b/Isosurfacing_3/examples/Isosurfacing_3/marching_cubes_implicit_sphere.cpp deleted file mode 100644 index 2c651ad808c..00000000000 --- a/Isosurfacing_3/examples/Isosurfacing_3/marching_cubes_implicit_sphere.cpp +++ /dev/null @@ -1,50 +0,0 @@ -#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 Point_range = std::vector; -using Triangle_range = std::vector >; - -// Sphere function = Euclidean distance function to the origin -auto sphere_function = [](const Point& p) -> FT -{ - return std::sqrt(p.x() * p.x() + p.y() * p.y() + p.z() * p.z()); -}; - -int main(int, char**) -{ - // box domain and spacing vector - const CGAL::Bbox_3 bbox{ -1.0, -1.0, -1.0, 1.0, 1.0, 1.0 }; - const FT spacing = 0.05; - const Vector vec_spacing(spacing, spacing, spacing); - - // create domain with sphere function - auto domain = CGAL::Isosurfacing::create_implicit_Cartesian_grid_domain - (bbox, vec_spacing, sphere_function); - - // points and triangles for the output indexed soup - Point_range points; - Triangle_range triangles; - - // run marching cubes with given isovalue - std::cout << "marching cubes..."; - const FT isovalue = 0.8; - CGAL::Timer timer; - timer.start(); - CGAL::Isosurfacing::marching_cubes(domain, isovalue, points, triangles); - timer.stop(); - std::cout << "done (" << timer.time() << "s, " << triangles.size() << " triangles)" << std::endl; - - // save ouput indexed mesh to a file, in the OFF format - CGAL::IO::write_OFF("marching_cubes_implicit_sphere.off", points, triangles); - - return EXIT_SUCCESS; -} diff --git a/Isosurfacing_3/examples/Isosurfacing_3/marching_cubes_inrimage.cpp b/Isosurfacing_3/examples/Isosurfacing_3/marching_cubes_inrimage.cpp deleted file mode 100644 index e76beb480b1..00000000000 --- a/Isosurfacing_3/examples/Isosurfacing_3/marching_cubes_inrimage.cpp +++ /dev/null @@ -1,52 +0,0 @@ -#include - -#include -#include -#include - -#include - -#include - -using Kernel = CGAL::Simple_cartesian; -using Point = typename Kernel::Point_3; -using Grid = CGAL::Isosurfacing::Cartesian_grid_3; - -using Point_range = std::vector; -using Polygon_range = std::vector >; - -int main(int argc, char* argv[]) -{ - const std::string fname = (argc > 1) ? argv[1] : CGAL::data_file_path("images/skull_2.9.inr"); - - // load volumetric image from a file - CGAL::Image_3 image; - if(!image.read(fname)) - { - std::cerr << "Error: Cannot read image file " << fname << std::endl; - return EXIT_FAILURE; - } - - // convert image to a Cartesian grid - Grid grid{image}; - - for (std::size_t i=0; i -#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 std::sqrt((p - x).squared_length()); -} - -// Usage : marching_cubes_multiple_mesh_offsets input.off -int main(int argc, char *argv[]) -{ - const std::string input_name = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/cross.off"); - - 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); - - // run Marching cubes with a range of offsets, - // and save all output meshes to files "output-index.off" - 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 triangles; - - // execute marching cubes with an isovalue equating offset - std::cout << "Marching cubes with offset " << offset << "..."; - CGAL::Isosurfacing::marching_cubes(domain, offset, points, triangles); - 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 << "Saving to file " << filename << "..."; - CGAL::IO::write_polygon_soup(filename, points, triangles); - std::cout << "done" << std::endl; - } - - return EXIT_SUCCESS; -} diff --git a/Isosurfacing_3/examples/Isosurfacing_3/marching_cubes_seq_vs_parallel.cpp b/Isosurfacing_3/examples/Isosurfacing_3/marching_cubes_seq_vs_parallel.cpp deleted file mode 100644 index 553ccc92be9..00000000000 --- a/Isosurfacing_3/examples/Isosurfacing_3/marching_cubes_seq_vs_parallel.cpp +++ /dev/null @@ -1,62 +0,0 @@ -#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 Point_range = std::vector; -using Triangle_range = std::vector >; - -// Sphere function = Euclidean distance function to the origin -auto sphere_function = [](const Point& p) -> FT -{ - return std::sqrt(p.x() * p.x() + p.y() * p.y() + p.z() * p.z()); -}; - -int main(int, char**) -{ - // box domain and spacing vector - const CGAL::Bbox_3 bbox{ -1.0, -1.0, -1.0, 1.0, 1.0, 1.0 }; - const FT spacing = 0.002; - const Vector vec_spacing(spacing, spacing, spacing); - - // create domain with sphere function - auto domain = CGAL::Isosurfacing::create_implicit_Cartesian_grid_domain - (bbox, vec_spacing, sphere_function); - - // points and triangles for the output indexed soup - Point_range points; - Triangle_range triangles; - - // run marching cubes sequential - std::cout << "marching cubes sequential..."; - const FT isovalue = 0.8; - CGAL::Real_timer timer; - timer.start(); - CGAL::Isosurfacing::marching_cubes(domain, isovalue, points, triangles); - timer.stop(); - std::cout << "done (" << timer.time() << "s, " << triangles.size() << " triangles)" << std::endl; - CGAL::IO::write_OFF("output-sequential.off", points, triangles); - - // clear points and triangles - points.clear(); - triangles.clear(); - - // run marching cubes parallel - std::cout << "marching cubes parallel..."; - timer.reset(); - timer.start(); - CGAL::Isosurfacing::marching_cubes(domain, isovalue, points, triangles); - timer.stop(); - std::cout << "done (" << timer.time() << "s, " << triangles.size() << " triangles)" << std::endl; - CGAL::IO::write_OFF("output-parallel.off", points, triangles); - - return EXIT_SUCCESS; -} diff --git a/Isosurfacing_3/examples/Isosurfacing_3/marching_cubes_signed_mesh_offset.cpp b/Isosurfacing_3/examples/Isosurfacing_3/marching_cubes_signed_mesh_offset.cpp deleted file mode 100644 index 432e0995e78..00000000000 --- a/Isosurfacing_3/examples/Isosurfacing_3/marching_cubes_signed_mesh_offset.cpp +++ /dev/null @@ -1,105 +0,0 @@ -#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 cp = tree.closest_point(p); - return sqrt((p - cp).squared_length()); -} - -int main(int argc, char* argv[]) -{ - const std::string input_name = (argc > 1) ? argv[1] :CGAL::data_file_path("meshes/cross.off"); - const int n_voxels = 20; - const FT offset_value = 0.2; - - // 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; - } - - // compute loose bounding box of the mesh - CGAL::Bbox_3 aabb_grid = CGAL::Polygon_mesh_processing::bbox(mesh_input); - const FT loose_offset = offset_value + 0.01; - 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 - Grid grid { n_voxels, n_voxels, n_voxels, aabb_grid }; - - for(std::size_t k=0; k