diff --git a/Convex_decomposition_3/doc/Convex_decomposition_3/Convex_decomposition_3.txt b/Convex_decomposition_3/doc/Convex_decomposition_3/Convex_decomposition_3.txt index af1dfe4f7cf..219f6cab9b6 100644 --- a/Convex_decomposition_3/doc/Convex_decomposition_3/Convex_decomposition_3.txt +++ b/Convex_decomposition_3/doc/Convex_decomposition_3/Convex_decomposition_3.txt @@ -130,12 +130,12 @@ The method employs a top-down splitting phase followed by a bottom-up merging to The volume calculation, convex hull computation and the concavity search is accelerated by a voxel grid. The grid is prepared before the splitting phase and is labelled into outside, inside or surface. The convex hulls are calculated from voxel corners. Thus, a mesh with a high resolution is less penalized by its number of vertices. The splitting phase typically results in a number of convex volumes larger than targeted. The second phase employs a bottom-up merging that reduces the number of convex volumes to the targeted number while aiming at maintaining a low volume difference between convex volumes and the input mesh. The greedy merging maintains a priority queue to incrementally merge the pair of convex volumes with the smallest increase of volume difference. -The splitting phase is not limited by the chosen `maximum_number_of_convex_hulls`, because a splitting into a larger number of more convex parts with a subsequent merging leads to better results. +The splitting phase is not limited by the chosen `maximum_number_of_convex_volumes`, because a splitting into a larger number of more convex parts with a subsequent merging leads to better results. \subsection Convex_decomposition_3ACD_Parameters Parameters Several parameters of the algorithm impact the quality of the result as well as the running time. -- `maximum_number_of_convex_hulls`: The maximum number of convex volumes output by the method. The actual number may be lower for mostly convex input meshes, e.g., a sphere. The impact on the running time is rather low. The default is 16. +- `maximum_number_of_convex_volumes`: The maximum number of convex volumes output by the method. The actual number may be lower for mostly convex input meshes, e.g., a sphere. The impact on the running time is rather low. The default is 16. - `maximum_depth`: The maximum depth for the hierarchical splitting phase. For complex meshes, a higher maximum depth is required to split small concavities into convex parts. The choice of `maximum_depth` has a larger impact on the running time. The default is 10. - `maximum_number_of_voxels`: This parameter controls the resolution of the voxel grid used for speed-up. Larger numbers result in a higher memory footprint and a higher running time. A small number also limits the `maximum_depth`. The voxel grid is isotropic and the longest axis of the bounding box will be split into a number of voxels equal to the cubic root of `maximum_number_of_voxels`. The default value is 1.000.000. - `volume_error`: The splitting of a convex volume into smaller parts is controlled by the `volume_error` which provides the tolerance for difference in volume. The difference is calculated by \f$ (|C_i| - |P_i|) / |P_i|\f$. The default value is 0.01. Thus, if a convex volume has 1 percent more volume that the part of the input mesh it approximates, it will be further divided. diff --git a/Lab/demo/Lab/Plugins/Convex_decomposition/Approximate_convex_decomposition_plugin.cpp b/Lab/demo/Lab/Plugins/Convex_decomposition/Approximate_convex_decomposition_plugin.cpp index 8e2671ff008..6d673295291 100644 --- a/Lab/demo/Lab/Plugins/Convex_decomposition/Approximate_convex_decomposition_plugin.cpp +++ b/Lab/demo/Lab/Plugins/Convex_decomposition/Approximate_convex_decomposition_plugin.cpp @@ -111,13 +111,13 @@ approximate_convex_decomposition() QApplication::setOverrideCursor(Qt::WaitCursor); - std::vector convex_hulls; - convex_hulls.reserve(9); + std::vector convex_volumes; + convex_volumes.reserve(9); - CGAL::Polygon_mesh_processing::approximate_convex_decomposition(*(sm_item->face_graph()), std::back_inserter(convex_hulls), + CGAL::Polygon_mesh_processing::approximate_convex_decomposition(*(sm_item->face_graph()), std::back_inserter(convex_volumes), CGAL::parameters::maximum_depth(maximumDepth) .volume_error(volumeError) - .maximum_number_of_convex_hulls(maximumConvexHulls) + .maximum_number_of_convex_volumes(maximumConvexHulls) .split_at_concavity(splitConcavity) .maximum_number_of_voxels(numVoxels) .concurrency_tag(CGAL::Parallel_if_available_tag())); @@ -126,10 +126,10 @@ approximate_convex_decomposition() std::vector distinct_colors; // the first color is either the background or the unique domain - compute_deterministic_color_map(QColor(80, 250, 80), convex_hulls.size(), std::back_inserter(distinct_colors)); + compute_deterministic_color_map(QColor(80, 250, 80), convex_volumes.size(), std::back_inserter(distinct_colors)); - for (std::size_t i = 0; i < convex_hulls.size(); i++) { - const Convex_hull& ch = convex_hulls[i]; + for (std::size_t i = 0; i < convex_volumes.size(); i++) { + const Convex_hull& ch = convex_volumes[i]; SMesh sm; CGAL::Polygon_mesh_processing::polygon_soup_to_polygon_mesh(ch.first, ch.second, sm); diff --git a/Polygon_mesh_processing/examples/Polygon_mesh_processing/CMakeLists.txt b/Polygon_mesh_processing/examples/Polygon_mesh_processing/CMakeLists.txt index 2fcaa4c6912..88430e902e0 100644 --- a/Polygon_mesh_processing/examples/Polygon_mesh_processing/CMakeLists.txt +++ b/Polygon_mesh_processing/examples/Polygon_mesh_processing/CMakeLists.txt @@ -59,7 +59,6 @@ create_single_source_cgal_program("isotropic_remeshing_with_allow_move.cpp") create_single_source_cgal_program("triangle_mesh_autorefinement.cpp") create_single_source_cgal_program("soup_autorefinement.cpp") create_single_source_cgal_program("snap_polygon_soup.cpp") -create_single_source_cgal_program("approximate_convex_decomposition.cpp") find_package(Eigen3 QUIET) include(CGAL_Eigen3_support) @@ -144,7 +143,6 @@ if(TARGET CGAL::TBB_support) target_link_libraries(hausdorff_bounded_error_distance_example PRIVATE CGAL::TBB_support) target_link_libraries(soup_autorefinement PRIVATE CGAL::TBB_support) target_link_libraries(snap_polygon_soup PRIVATE CGAL::TBB_support) - target_link_libraries(approximate_convex_decomposition PRIVATE CGAL::TBB_support) create_single_source_cgal_program("corefinement_parallel_union_meshes.cpp") target_link_libraries(corefinement_parallel_union_meshes PRIVATE CGAL::TBB_support) diff --git a/Polygon_mesh_processing/examples/Polygon_mesh_processing/approximate_convex_decomposition.cpp b/Polygon_mesh_processing/examples/Polygon_mesh_processing/approximate_convex_decomposition.cpp deleted file mode 100644 index 1b3e34005d4..00000000000 --- a/Polygon_mesh_processing/examples/Polygon_mesh_processing/approximate_convex_decomposition.cpp +++ /dev/null @@ -1,47 +0,0 @@ -#include -#include - -#include -#include - -#include -#include -#include - -using K = CGAL::Exact_predicates_inexact_constructions_kernel; - -using Point = K::Point_3; - -using Convex_hull = std::pair, std::vector > >; -using Mesh = CGAL::Surface_mesh; -namespace PMP = CGAL::Polygon_mesh_processing; - -int main(int argc, char* argv[]) -{ - const std::string filename = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/knot2.off"); - std::cout << filename << std::endl; - - Mesh mesh; - if(!PMP::IO::read_polygon_mesh(filename, mesh)) { - std::cerr << "Invalid input." << std::endl; - return 1; - } - - std::vector convex_hulls; - convex_hulls.reserve(9); - - PMP::approximate_convex_decomposition(mesh, std::back_inserter(convex_hulls), - CGAL::parameters::maximum_depth(10) - .volume_error(0.1) - .maximum_number_of_convex_hulls(9) - .split_at_concavity(true) - .maximum_number_of_voxels(1000000) - .concurrency_tag(CGAL::Parallel_if_available_tag())); - - for (std::size_t i = 0;i - -#include -#include - -#include - -#include -#include - -#include -#include -#include -#include -#include -#include -#include - -#include -#include -#include - -#ifdef CGAL_LINKED_WITH_TBB -#include -#include -#include -#include -#else -#include -#endif - -namespace CGAL { -namespace Polygon_mesh_processing { -namespace internal { - -using Vec3_uint = std::array; - -struct Bbox_uint { - Vec3_uint lower; - Vec3_uint upper; - Bbox_uint(const Vec3_uint &lower, const Vec3_uint &upper) : lower(lower), upper(upper) {} -}; - -enum Grid_cell : int8_t { - OUTSIDE = -1, - SURFACE = 0, - INSIDE = 1 -}; - -void export_grid(const std::string& filename, const Bbox_3& bb, std::vector& grid, const Vec3_uint& grid_size, double voxel_size) { - const auto vox = [&grid, &grid_size](unsigned int x, unsigned int y, unsigned int z) -> int8_t& { - return grid[z + (y * grid_size[2]) + (x * grid_size[1] * grid_size[2])]; - }; - std::ofstream stream(filename); - - stream << - "ply\n" << - "format ascii 1.0\n" << - "element vertex " << (grid_size[0] * grid_size[1] * grid_size[2]) << "\n" << - "property double x\n" << - "property double y\n" << - "property double z\n" << - "property uchar red\n" << - "property uchar green\n" << - "property uchar blue\n" << - "end_header\n"; - - for (unsigned int x = 0; x < grid_size[0]; x++) - for (unsigned int y = 0; y < grid_size[1]; y++) - for (unsigned int z = 0; z < grid_size[2]; z++) { - stream << (bb.xmin() + (x + 0.5) * voxel_size) << " " << (bb.ymin() + (y + 0.5) * voxel_size) << " " << (bb.zmin() + (z + 0.5) * voxel_size) << " "; - switch (vox(x, y, z)) { - case INSIDE: - stream << "175 175 100\n"; - break; - case OUTSIDE: - stream << "125 125 175\n"; - break; - case SURFACE: - stream << "200 100 100\n"; - break; - default: - stream << "0 0 0\n"; - break; - } - } - - stream.close(); -} - -template -void export_grid(const std::string& filename, const Bbox_3& bb, std::vector& grid, const Vec3_uint& grid_size, double voxel_size, Filter &filter) { - const auto vox = [&grid, &grid_size](unsigned int x, unsigned int y, unsigned int z) -> int8_t& { - return grid[z + (y * grid_size[2]) + (x * grid_size[1] * grid_size[2])]; - }; - std::ofstream stream(filename); - - std::size_t count = 0; - for (unsigned int x = 0; x < grid_size[0]; x++) - for (unsigned int y = 0; y < grid_size[1]; y++) - for (unsigned int z = 0; z < grid_size[2]; z++) - if (filter(vox(x, y, z))) count++; - - stream << - "ply\n" << - "format ascii 1.0\n" << - "element vertex " << count << "\n" << - "property double x\n" << - "property double y\n" << - "property double z\n" << - "property uchar red\n" << - "property uchar green\n" << - "property uchar blue\n" << - "end_header\n"; - - for (unsigned int x = 0; x < grid_size[0]; x++) - for (unsigned int y = 0; y < grid_size[1]; y++) - for (unsigned int z = 0; z < grid_size[2]; z++) { - if (!filter(vox(x, y, z))) continue; - stream << (bb.xmin() + (x + 0.5) * voxel_size) << " " << (bb.ymin() + (y + 0.5) * voxel_size) << " " << (bb.zmin() + (z + 0.5) * voxel_size) << " "; - switch (vox(x, y, z)) { - case INSIDE: - stream << "175 175 100\n"; - break; - case OUTSIDE: - stream << "125 125 175\n"; - break; - case SURFACE: - stream << "200 100 100\n"; - break; - default: - stream << "0 0 0\n"; - break; - } - } - - stream.close(); -} - -void export_voxels(const std::string& filename, const Bbox_3& bb, std::vector& voxels, double voxel_size) { - std::ofstream stream(filename); - - stream << - "ply\n" << - "format ascii 1.0\n" << - "element vertex " << voxels.size() << "\n" << - "property double x\n" << - "property double y\n" << - "property double z\n" << - "end_header\n"; - - for (const Vec3_uint& v : voxels) { - stream << (bb.xmin() + (v[0] + 0.5) * voxel_size) << " " << (bb.ymin() + (v[1] + 0.5) * voxel_size) << " " << (bb.zmin() + (v[2] + 0.5) * voxel_size) << "\n"; - } - - stream.close(); -} - -template -void export_points(const std::string& filename, const Bbox_3& bb, std::vector& points) { - std::ofstream stream(filename); - - stream << - "ply\n" << - "format ascii 1.0\n" << - "element vertex " << points.size() << "\n" << - "property double x\n" << - "property double y\n" << - "property double z\n" << - "end_header\n"; - - for (const Point_3& p : points) { - stream << (bb.xmin() + p.x()) << " " << (bb.ymin() + p.y()) << " " << (bb.zmin() + p.z()) << "\n"; - } - - stream.close(); -} - -template -void export_points(const std::string& filename, Range& points) { - std::ofstream stream(filename); - stream << std::setprecision(18); - - for (const auto& p : points) { - stream << p.x() << " " << p.y() << " " << p.z() << "\n"; - } - - stream.close(); -} - -template -Box box_union(const Box& a, const Box& b) { - using FT = decltype(a.xmin()); - return Box( - (std::min)(a.xmin(), b.xmin()), - (std::min)(a.ymin(), b.ymin()), - (std::min)(a.zmin(), b.zmin()), - (std::max)(a.xmax(), b.xmax()), - (std::max)(a.ymax(), b.ymax()), - (std::max)(a.zmax(), b.zmax())); -} - -template -std::tuple calculate_grid_size(Bbox_3& bb, const unsigned int number_of_voxels) { - std::size_t max_voxels_axis = std::cbrt(number_of_voxels); - assert(max_voxels_axis > 3); - // get longest axis - FT longest = 0; - - if (bb.x_span() >= bb.y_span() && bb.x_span() >= bb.z_span()) - longest = bb.x_span(); - else if (bb.y_span() >= bb.x_span() && bb.y_span() >= bb.z_span()) - longest = bb.y_span(); - else if (bb.z_span() >= bb.x_span() && bb.z_span() >= bb.y_span()) - longest = bb.z_span(); - - const FT voxel_size = longest * FT(1.0 / (max_voxels_axis - 3)); - - FT s = 1.5 * voxel_size; - - bb = Bbox_3(to_double(bb.xmin() - s), to_double(bb.ymin() - s), to_double(bb.zmin() - s), to_double(bb.xmax() + s), to_double(bb.ymax() + s), to_double(bb.zmax() + s)); - - return { Vec3_uint{static_cast(to_double(bb.x_span() / voxel_size + 0.5)), static_cast(to_double(bb.y_span() / voxel_size + 0.5)), static_cast(to_double(bb.z_span() / voxel_size + 0.5))}, voxel_size}; -} - -template -const typename GeomTraits::Point_3 &point(typename boost::graph_traits::face_descriptor fd, - const PolygonMesh& pmesh, - const NamedParameters& np = parameters::default_values()) -{ - using parameters::choose_parameter; - using parameters::get_parameter; - typename GetVertexPointMap::const_type - vpm = choose_parameter(get_parameter(np, internal_np::vertex_point), - get_const_property_map(CGAL::vertex_point, pmesh)); - - return get(vpm, target(halfedge(fd, pmesh), pmesh)); -} - -template -Bbox_uint grid_bbox_face(const FaceGraph& mesh, const typename boost::graph_traits::face_descriptor fd, const Bbox_3& bb, const FT& voxel_size) { - Bbox_3 face_bb = face_bbox(fd, mesh); - double vs = to_double(voxel_size); - return Bbox_uint({ - static_cast((face_bb.xmin() - bb.xmin()) / vs - 0.5), - static_cast((face_bb.ymin() - bb.ymin()) / vs - 0.5), - static_cast((face_bb.zmin() - bb.zmin()) / vs - 0.5) - }, { - static_cast((face_bb.xmax() - bb.xmin()) / vs + 0.5), - static_cast((face_bb.ymax() - bb.ymin()) / vs + 0.5), - static_cast((face_bb.zmax() - bb.zmin()) / vs + 0.5) - }); -} - -template -Iso_cuboid_3 bbox_voxel(const Vec3_uint& voxel, const Bbox_3& bb, const typename GeomTraits::FT& voxel_size) { - double vs = to_double(voxel_size); - return Bbox_3( - bb.xmin() + voxel[0] * vs, - bb.ymin() + voxel[1] * vs, - bb.zmin() + voxel[2] * vs, - bb.xmin() + (voxel[0] + 1) * vs, - bb.ymin() + (voxel[1] + 1) * vs, - bb.zmin() + (voxel[2] + 1) * vs - ); -} - -void scanline_floodfill(Grid_cell label, std::vector& grid, const Vec3_uint& grid_size, std::deque& todo) { - const auto vox = [&grid, &grid_size](unsigned int x, unsigned int y, unsigned int z) -> int8_t& { - return grid[z + (y * grid_size[2]) + (x * grid_size[1] * grid_size[2])]; - }; - - while (!todo.empty()) { - auto [x, y, z0] = todo.front(); - todo.pop_front(); - if (vox(x, y, z0) != INSIDE) - return; - - bool xneg = false, xpos = false; - bool yneg = false, ypos = false; - bool zneg = false, zpos = false; - - // positive direction - for (unsigned int z = z0; z < grid_size[2]; z++) { - if (vox(x, y, z) != INSIDE) - break; - - vox(x, y, z) = label; - if (x > 0) { - if (vox(x - 1, y, z) == INSIDE) { - if (!xneg) { - xneg = true; - todo.push_back({ x - 1, y, z }); - } - } - else xneg = false; - } - - if (x < grid_size[0] - 1) { - if (vox(x + 1, y, z) == INSIDE) { - if (!xpos) { - xpos = true; - todo.push_back({ x + 1, y, z }); - } - } - else xpos = false; - } - - if (y > 0) { - if (vox(x, y - 1, z) == INSIDE) { - if (!yneg) { - yneg = true; - todo.push_front({ x, y - 1, z }); - } - } - else yneg = false; - } - - if (y < grid_size[1] - 1) { - if (vox(x, y + 1, z) == INSIDE) { - if (!ypos) { - ypos = true; - todo.push_front({ x, y + 1, z }); - } - } - else ypos = false; - } - } - - xneg = xpos = yneg = ypos = zneg = zpos = false; - - if (z0 == 0) - continue; - - for (unsigned int z = z0 - 1; z > 0; z--) { - if (vox(x, y, z) != INSIDE) - break; - - vox(x, y, z) = label; - if (x > 0) { - if (vox(x - 1, y, z) == INSIDE) { - if (!xneg) { - xneg = true; - todo.push_back({ x - 1, y, z }); - } - } - else xneg = false; - } - - if (x < grid_size[0] - 1) { - if (vox(x + 1, y, z) == INSIDE) { - if (!xpos) { - xpos = true; - todo.push_back({ x + 1, y, z }); - } - } - else xpos = false; - } - - if (y > 0) { - if (vox(x, y - 1, z) == INSIDE) { - if (!yneg) { - yneg = true; - todo.push_front({ x, y - 1, z }); - } - } - else yneg = false; - } - - if (y < grid_size[1] - 1) { - if (vox(x, y + 1, z) == INSIDE) { - if (!ypos) { - ypos = true; - todo.push_front({ x, y + 1, z }); - } - } - else ypos = false; - } - } - } -} - -// Valid voxel grids separate OUTSIDE from INSIDE via SURFACE -bool is_valid(std::vector& grid, const Vec3_uint& grid_size) { - const auto vox = [&grid, &grid_size](unsigned int x, unsigned int y, unsigned int z) -> int8_t& { - return grid[z + (y * grid_size[2]) + (x * grid_size[1] * grid_size[2])]; - }; - - std::size_t in = 0, out = 0, surface = 0, other = 0, violations = 0; - - for (unsigned int x = 1; x < grid_size[0] - 1; x++) - for (unsigned int y = 1; y < grid_size[1] - 1; y++) - for (unsigned int z = 1; z < grid_size[2] - 1; z++) { - switch (vox(x, y, z)) { - case INSIDE: - if ( vox(x - 1, y, z) == OUTSIDE - || vox(x + 1, y, z) == OUTSIDE - || vox(x, y - 1, z) == OUTSIDE - || vox(x, y + 1, z) == OUTSIDE - || vox(x, y, z - 1) == OUTSIDE - || vox(x, y, z + 1) == OUTSIDE) { - std::cout << "touching I-O: " << x << " " << y << " " << z << std::endl; - violations++; - } - in++; - break; - case OUTSIDE: - if (vox(x - 1, y, z) == INSIDE - || vox(x + 1, y, z) == INSIDE - || vox(x, y - 1, z) == INSIDE - || vox(x, y + 1, z) == INSIDE - || vox(x, y, z - 1) == INSIDE - || vox(x, y, z + 1) == INSIDE) { - std::cout << "touching O-I: " << x << " " << y << " " << z << std::endl; - violations++; - } - out++; - break; - case SURFACE: - surface++; - break; - default: - std::cout << "other " << x << " " << y << " " << z << std::endl; - other++; - break; - } - } - std::cout << "i: " << in << " o: " << out << " s: " << surface << " " << other << std::endl; - std::cout << "violations: " << violations << std::endl; - - return violations == 0; -} - -// Only works for closed meshes -void label_floodfill(std::vector& grid, const Vec3_uint& grid_size) { - // Walk around boundary and start floodfill when voxel label is INSIDE - const auto vox = [&grid, &grid_size](unsigned int x, unsigned int y, unsigned int z) -> int8_t& { - return grid[z + (y * grid_size[2]) + (x * grid_size[1] * grid_size[2])]; - }; - - std::deque todo; - - // xmin/xmax - for (unsigned int y = 0; y < grid_size[1]; y++) - for (unsigned int z = 0; z < grid_size[2]; z++) { - if (vox(0, y, z) == INSIDE) { - todo.push_front({0, y, z}); - scanline_floodfill(OUTSIDE, grid, grid_size, todo); - } - - if (vox(grid_size[0] - 1, y, z) == INSIDE) { - todo.push_front({ grid_size[0] - 1, y, z }); - scanline_floodfill(OUTSIDE, grid, grid_size, todo); - } - } - - // ymin/ymax - for (unsigned int x = 0; x < grid_size[0]; x++) - for (unsigned int z = 0; z < grid_size[2]; z++) { - if (vox(x, 0, z) == INSIDE) { - todo.push_front({ x, 0, z }); - scanline_floodfill(OUTSIDE, grid, grid_size, todo); - } - - if (vox(x, grid_size[1] - 1, z) == INSIDE) { - todo.push_front({ x, grid_size[1] - 1, z }); - scanline_floodfill(OUTSIDE, grid, grid_size, todo); - } - } - - // ymin/ymax - for (unsigned int x = 0; x < grid_size[0]; x++) - for (unsigned int y = 0; y < grid_size[1]; y++) { - if (vox(x, y, 0) == INSIDE) { - todo.push_front({ x, y, 0 }); - scanline_floodfill(OUTSIDE, grid, grid_size, todo); - } - - if (vox(x, y, grid_size[2] - 1) == INSIDE) { - todo.push_front({ x, y, grid_size[2] - 1 }); - scanline_floodfill(OUTSIDE, grid, grid_size, todo); - } - } -} - -void naive_floodfill(std::vector& grid, const Vec3_uint& grid_size) { - const auto vox = [&grid, &grid_size](unsigned int x, unsigned int y, unsigned int z) -> int8_t& { - return grid[z + (y * grid_size[2]) + (x * grid_size[1] * grid_size[2])]; - }; - - std::queue queue; - queue.push({0, 0, 0}); - - while (!queue.empty()) { - auto [x, y, z] = queue.front(); - queue.pop(); - - if (vox(x, y, z) != INSIDE) - continue; - - vox(x, y, z) = OUTSIDE; - if (x > 0) - if (vox(x - 1, y, z) == INSIDE) { - queue.push({ x - 1, y, z }); - } - - if (x < grid_size[0] - 1) - if (vox(x + 1, y, z) == INSIDE) { - queue.push({ x + 1, y, z }); - } - - if (y > 0) - if (vox(x, y - 1, z) == INSIDE) { - queue.push({ x, y - 1, z }); - } - - if (y < grid_size[1] - 1) - if (vox(x, y + 1, z) == INSIDE) { - queue.push({ x, y + 1, z }); - } - - if (z > 0) - if (vox(x, y, z - 1) == INSIDE) { - queue.push({ x, y, z - 1 }); - } - - if (z < grid_size[2] - 1) - if (vox(x, y, z + 1) == INSIDE) { - queue.push({ x, y, z + 1 }); - } - } -} - -template -void rayshooting_fill(std::vector& grid, const Vec3_uint& grid_size, const Bbox_3& bb, const typename GeomTraits::FT& voxel_size, const FaceGraph& mesh, CGAL::Parallel_tag) { -#ifdef CGAL_LINKED_WITH_TBB - const auto vox = [&grid, &grid_size](unsigned int x, unsigned int y, unsigned int z) -> int8_t& { - return grid[z + (y * grid_size[2]) + (x * grid_size[1] * grid_size[2])]; - }; - using face_descriptor = typename boost::graph_traits::face_descriptor; - - using Point_3 = typename GeomTraits::Point_3; - using Vector_3 = typename GeomTraits::Vector_3; - using Ray_3 = typename GeomTraits::Ray_3; - - using Primitive = CGAL::AABB_face_graph_triangle_primitive; - using Traits = CGAL::AABB_traits_3; - using Tree = CGAL::AABB_tree; - using Ray_intersection = std::optional::Type>; - - Tree tree(faces(mesh).first, faces(mesh).second, mesh); - - std::array dirs = { Vector_3{1, 0, 0}, Vector_3{-1, 0, 0}, Vector_3{0, 1, 0}, Vector_3{0,-1, 0}, Vector_3{0, 0, 1}, Vector_3{0, 0,-1} }; - - tbb::parallel_for(std::size_t(0), std::size_t(grid_size[0]), [&](const std::size_t x) - { - for (std::size_t y = 0; y < grid_size[1]; y++) - for (std::size_t z = 0; z < grid_size[2]; z++) { - if (vox(x, y, z) == SURFACE) - continue; - - Point_3 c(bb.xmin() + (x + 0.5) * voxel_size, bb.ymin() + (y + 0.5) * voxel_size, bb.zmin() + (z + 0.5) * voxel_size); - - unsigned int inside = 0; - unsigned int outside = 0; - - for (std::size_t i = 0; i < 6; i++) { - Ray_intersection intersection = tree.first_intersection(Ray_3(c, dirs[i])); - if (intersection) { - // A segment intersection is not helpful as it means the triangle normal is orthogonal to the ray - if (std::get_if(&(intersection->first))) { - face_descriptor fd = intersection->second; - Vector_3 n = compute_face_normal(fd, mesh); - if (dirs[i] * n > 0) - inside++; - else - outside++; - } - } - } - - if (inside >= 3 && outside == 0) { - vox(x, y, z) = INSIDE; - } - else - vox(x, y, z) = OUTSIDE; - } - } - ); -#else - CGAL_USE(grid); - CGAL_USE(grid_size); - CGAL_USE(bb); - CGAL_USE(voxel_size); - CGAL_USE(mesh); -#endif -} - -template -void rayshooting_fill(std::vector& grid, const Vec3_uint& grid_size, const Bbox_3& bb, const typename GeomTraits::FT& voxel_size, const FaceGraph& mesh, CGAL::Sequential_tag) { - const auto vox = [&grid, &grid_size](unsigned int x, unsigned int y, unsigned int z) -> int8_t& { - return grid[z + (y * grid_size[2]) + (x * grid_size[1] * grid_size[2])]; - }; - using face_descriptor = typename boost::graph_traits::face_descriptor; - - using Point_3 = typename GeomTraits::Point_3; - using Vector_3 = typename GeomTraits::Vector_3; - using Ray_3 = typename GeomTraits::Ray_3; - - using Primitive = CGAL::AABB_face_graph_triangle_primitive; - using Traits = CGAL::AABB_traits_3; - using Tree = CGAL::AABB_tree; - using Ray_intersection = std::optional::Type>; - - Tree tree(faces(mesh).first, faces(mesh).second, mesh); - - std::array dirs = { Vector_3{1, 0, 0}, Vector_3{-1, 0, 0}, Vector_3{0, 1, 0}, Vector_3{0,-1, 0}, Vector_3{0, 0, 1}, Vector_3{0, 0,-1} }; - - for (std::size_t x = 0; x < grid_size[0]; x++) - { - for (std::size_t y = 0; y < grid_size[1]; y++) - for (std::size_t z = 0; z < grid_size[2]; z++) { - if (vox(x, y, z) == SURFACE) - continue; - - Point_3 c(bb.xmin() + (x + 0.5) * voxel_size, bb.ymin() + (y + 0.5) * voxel_size, bb.zmin() + (z + 0.5) * voxel_size); - - unsigned int inside = 0; - unsigned int outside = 0; - - for (std::size_t i = 0; i < 6; i++) { - Ray_intersection intersection = tree.first_intersection(Ray_3(c, dirs[i])); - if (intersection) { - // A segment intersection is not helpful as it means the triangle normal is orthogonal to the ray - if (std::get_if(&(intersection->first))) { - face_descriptor fd = intersection->second; - Vector_3 n = compute_face_normal(fd, mesh); - if (dirs[i] * n > 0) - inside++; - else - outside++; - } - } - } - - if (inside >= 3 && outside == 0) { - vox(x, y, z) = INSIDE; - } - else - vox(x, y, z) = OUTSIDE; - } - } -} - -template -struct Convex_hull_candidate { - using FT = typename GeomTraits::FT; - using Point_3 = typename GeomTraits::Point_3; - Iso_cuboid_3 bbox; - FT voxel_volume; - FT volume; - FT volume_error; - std::vector points; - std::vector> indices; - - Convex_hull_candidate() noexcept : voxel_volume(0), volume(0), volume_error(0) {} - Convex_hull_candidate(Convex_hull_candidate&& o) noexcept { - bbox = o.bbox; - voxel_volume = o.voxel_volume; - volume = o.volume; - volume_error = o.volume_error; - points = std::move(o.points); - indices = std::move(o.indices); - } - - Convex_hull_candidate& operator= (Convex_hull_candidate&& o) noexcept { - bbox = o.bbox; - voxel_volume = o.voxel_volume; - volume = o.volume; - volume_error = o.volume_error; - points = std::move(o.points); - indices = std::move(o.indices); - - return *this; - } -}; - -template -struct Candidate { - using FT = typename GeomTraits::FT; - using Point_3 = typename GeomTraits::Point_3; - std::size_t index; - std::vector surface; - std::vector new_surface; - std::vector inside; - std::size_t depth; - Bbox_uint bbox; - Convex_hull_candidate ch; - - Candidate() : depth(0), bbox({ 0, 0, 0 }, { 0, 0, 0 }) {index = cidx++;} - Candidate(std::size_t depth, Bbox_uint &bbox) : depth(depth), bbox(bbox) { index = cidx++; } -private: - inline static std::size_t cidx = 0; -}; - -template -typename GeomTraits::FT volume(const std::vector &pts, const std::vector> &indices) { - ::CGAL::internal::Evaluate evaluate; - - typename GeomTraits::FT vol = 0; - typename GeomTraits::Compute_volume_3 cv3; - - typename GeomTraits::Point_3 origin(0, 0, 0); - - for (const std::array &i : indices) { - vol += cv3(origin, pts[i[0]], pts[i[1]], pts[i[2]]); - evaluate(vol); - } - - return vol; -} - -void enlarge(Bbox_uint& bbox, const Vec3_uint& v) { - bbox.lower[0] = (std::min)(bbox.lower[0], v[0]); - bbox.lower[1] = (std::min)(bbox.lower[1], v[1]); - bbox.lower[2] = (std::min)(bbox.lower[2], v[2]); - bbox.upper[0] = (std::max)(bbox.upper[0], v[0]); - bbox.upper[1] = (std::max)(bbox.upper[1], v[1]); - bbox.upper[2] = (std::max)(bbox.upper[2], v[2]); -} - -template > -struct is_std_hashable : std::false_type {}; - -template -struct is_std_hashable>()(std::declval()))>> : std::true_type {}; - -template::type::FT>, typename std::is_same::type::Kernel_tag, CGAL::Cartesian_tag>::type>::type> -struct Unordered_set_if_available { -}; - -template -struct Unordered_set_if_available { - using type = std::unordered_set; -}; - -template -struct Unordered_set_if_available { - using type = std::set; -}; - -template -void compute_candidate(Candidate &c, const Bbox_3& bb, typename GeomTraits::FT voxel_size) { - using Point_3 = typename GeomTraits::Point_3; - using FT = typename GeomTraits::FT; - - c.bbox.lower = c.bbox.upper = c.surface[0]; - - typename Unordered_set_if_available::type voxel_points; - - for (const Vec3_uint& v : c.surface) { - enlarge(c.bbox, v); - FT xmin = bb.xmin() + FT(v[0]) * voxel_size; - FT ymin = bb.ymin() + FT(v[1]) * voxel_size; - FT zmin = bb.zmin() + FT(v[2]) * voxel_size; - FT xmax = bb.xmin() + FT(v[0] + 1) * voxel_size; - FT ymax = bb.ymin() + FT(v[1] + 1) * voxel_size; - FT zmax = bb.zmin() + FT(v[2] + 1) * voxel_size; - voxel_points.insert(Point_3(xmin, ymin, zmin)); - voxel_points.insert(Point_3(xmin, ymax, zmin)); - voxel_points.insert(Point_3(xmin, ymin, zmax)); - voxel_points.insert(Point_3(xmin, ymax, zmax)); - voxel_points.insert(Point_3(xmax, ymin, zmin)); - voxel_points.insert(Point_3(xmax, ymax, zmin)); - voxel_points.insert(Point_3(xmax, ymin, zmax)); - voxel_points.insert(Point_3(xmax, ymax, zmax)); - } - - for (const Vec3_uint& v : c.new_surface) { - enlarge(c.bbox, v); - FT xmin = bb.xmin() + FT(v[0]) * voxel_size; - FT ymin = bb.ymin() + FT(v[1]) * voxel_size; - FT zmin = bb.zmin() + FT(v[2]) * voxel_size; - FT xmax = bb.xmin() + FT(v[0] + 1) * voxel_size; - FT ymax = bb.ymin() + FT(v[1] + 1) * voxel_size; - FT zmax = bb.zmin() + FT(v[2] + 1) * voxel_size; - voxel_points.insert(Point_3(xmin, ymin, zmin)); - voxel_points.insert(Point_3(xmin, ymax, zmin)); - voxel_points.insert(Point_3(xmin, ymin, zmax)); - voxel_points.insert(Point_3(xmin, ymax, zmax)); - voxel_points.insert(Point_3(xmax, ymin, zmin)); - voxel_points.insert(Point_3(xmax, ymax, zmin)); - voxel_points.insert(Point_3(xmax, ymin, zmax)); - voxel_points.insert(Point_3(xmax, ymax, zmax)); - } - - convex_hull_3(voxel_points.begin(), voxel_points.end(), c.ch.points, c.ch.indices); - - c.ch.volume = volume(c.ch.points, c.ch.indices); - - CGAL_assertion(c.ch.volume > 0); - c.ch.bbox = bbox_3(c.ch.points.begin(), c.ch.points.end()); - - c.ch.voxel_volume = (voxel_size * voxel_size * voxel_size) * FT(double(c.inside.size() + c.surface.size() + c.new_surface.size())); - c.ch.volume_error = CGAL::abs(c.ch.volume - c.ch.voxel_volume) / c.ch.voxel_volume; -} - -template -void fill_grid(Candidate &c, std::vector &grid, const FaceGraph &mesh, const Bbox_3& bb, const Vec3_uint& grid_size, const typename GeomTraits::FT& voxel_size, Concurrency_tag tag) { - const auto vox = [&grid, &grid_size](unsigned int x, unsigned int y, unsigned int z) -> int8_t& { - return grid[z + (y * grid_size[2]) + (x * grid_size[1] * grid_size[2])]; - }; - - for (const typename boost::graph_traits::face_descriptor fd : faces(mesh)) { - Bbox_uint face_bb = grid_bbox_face(mesh, fd, bb, voxel_size); - CGAL_assertion(face_bb.lower[0] <= face_bb.upper[0]); - CGAL_assertion(face_bb.lower[1] <= face_bb.upper[1]); - CGAL_assertion(face_bb.lower[2] <= face_bb.upper[2]); - CGAL_assertion(face_bb.upper[0] < grid_size[0]); - CGAL_assertion(face_bb.upper[1] < grid_size[1]); - CGAL_assertion(face_bb.upper[2] < grid_size[2]); - for (unsigned int x = face_bb.lower[0]; x <= face_bb.upper[0]; x++) - for (unsigned int y = face_bb.lower[1]; y <= face_bb.upper[1]; y++) - for (unsigned int z = face_bb.lower[2]; z <= face_bb.upper[2]; z++) { - Iso_cuboid_3 box = bbox_voxel({ x, y, z }, bb, voxel_size); - const typename GeomTraits::Point_3 &p = point(fd, mesh); - if (do_intersect(triangle(fd, mesh), box) || box.has_on_bounded_side(p)) - vox(x, y, z) = Grid_cell::SURFACE; - } - } - - if (CGAL::is_closed(mesh)) - naive_floodfill(grid, grid_size); - else - rayshooting_fill(grid, grid_size, bb, voxel_size, mesh, tag); - - c.bbox.upper = {grid_size[0] - 1, grid_size[1] - 1, grid_size[2] - 1}; - - for (unsigned int x = 0; x < grid_size[0]; x++) - for (unsigned int y = 0; y < grid_size[1]; y++) - for (unsigned int z = 0; z < grid_size[2]; z++) - if (vox(x, y, z) == INSIDE) - c.inside.push_back({x, y, z}); - else if (vox(x, y, z) == SURFACE) - c.surface.push_back({x, y, z}); -} - -template -void init(Candidate &c, const FaceGraph& mesh, std::vector& grid, const Bbox_3& bb, const Vec3_uint& grid_size, const typename GeomTraits::FT& voxel_size, Concurrency_tag tag) { - internal::fill_grid(c, grid, mesh, bb, grid_size, voxel_size, tag); - compute_candidate(c, bb, voxel_size); -} - -template -void split(std::vector &candidates, Candidate_& c, unsigned int axis, unsigned int location) { - //Just split the voxel bbox along 'axis' after voxel index 'location' - Candidate_ upper(c.depth + 1, c.bbox); - Candidate_ lower(c.depth + 1, c.bbox); - - CGAL_assertion(c.bbox.lower[axis] < location); - CGAL_assertion(c.bbox.upper[axis] > location); - - upper.bbox.lower[axis] = location + 1; - lower.bbox.upper[axis] = location; - - for (const Vec3_uint& v : c.surface) { - CGAL_assertion(c.bbox.lower[0] <= v[0] && v[0] <= c.bbox.upper[0]); - CGAL_assertion(c.bbox.lower[1] <= v[1] && v[1] <= c.bbox.upper[1]); - CGAL_assertion(c.bbox.lower[2] <= v[2] && v[2] <= c.bbox.upper[2]); - if (location < v[axis]) - upper.surface.push_back(v); - else - lower.surface.push_back(v); - } - - for (const Vec3_uint& v : c.new_surface) { - CGAL_assertion(c.bbox.lower[0] <= v[0] && v[0] <= c.bbox.upper[0]); - CGAL_assertion(c.bbox.lower[1] <= v[1] && v[1] <= c.bbox.upper[1]); - CGAL_assertion(c.bbox.lower[2] <= v[2] && v[2] <= c.bbox.upper[2]); - if (location < v[axis]) - upper.new_surface.push_back(v); - else - lower.new_surface.push_back(v); - } - - for (const Vec3_uint& v : c.inside) { - CGAL_assertion(c.bbox.lower[0] <= v[0] && v[0] <= c.bbox.upper[0]); - CGAL_assertion(c.bbox.lower[1] <= v[1] && v[1] <= c.bbox.upper[1]); - CGAL_assertion(c.bbox.lower[2] <= v[2] && v[2] <= c.bbox.upper[2]); - if (location < v[axis]) { - if ((location + 1) == v[axis]) - upper.new_surface.push_back(v); - else - upper.inside.push_back(v); - } - else { - if (location == v[axis]) - lower.new_surface.push_back(v); - else - lower.inside.push_back(v); - } - } - - if (!upper.surface.empty()) - candidates.emplace_back(std::move(upper)); - - if (!lower.surface.empty()) - candidates.emplace_back(std::move(lower)); -} - -std::size_t concavity(const Vec3_uint& s, const Vec3_uint& e, int axis, std::vector& grid, const Vec3_uint& grid_size) { - const auto vox = [&grid, &grid_size](const Vec3_uint &v) -> int8_t& { - return grid[v[2] + (v[1] * grid_size[2]) + (v[0] * grid_size[1] * grid_size[2])]; - }; - std::size_t i; - for (i = s[axis];i s[axis]; j--) { - Vec3_uint v = s; - v[axis] = j; - if (vox(v) != OUTSIDE) - break; - } - - std::size_t res = (i - s[axis]) + (e[axis] - j); - if(res >= grid_size[axis]) - std::cout << "violation!" << std::endl; - return (i - s[axis]) + (e[axis] - j); -} - -void choose_splitting_location_by_concavity(unsigned int& axis, unsigned int& location, const Bbox_uint &bbox, std::vector& grid, const Vec3_uint& grid_size) { - std::size_t length = bbox.upper[axis] - bbox.lower[axis] + 1; - - CGAL_assertion(length >= 8); - CGAL_precondition(axis <= 2); - - std::size_t idx0 = 0, idx1 = 1, idx2 = 2; - - switch(axis) { - case 0: - idx0 = 1; - idx1 = 2; - idx2 = 0; - break; - case 1: - idx0 = 0; - idx1 = 2; - idx2 = 1; - break; - case 2: - idx0 = 0; - idx1 = 1; - idx2 = 2; - break; - } - - std::vector diam(length, 0), diam2(length, 0); - - for (std::size_t i = bbox.lower[idx2]; i <= bbox.upper[idx2]; i++) { - for (std::size_t j = bbox.lower[idx0]; j <= bbox.upper[idx0]; j++) { - Vec3_uint s, e; - s[idx2] = i; - s[idx1] = bbox.lower[idx1]; - s[idx0] = j; - e[idx2] = i; - e[idx1] = bbox.upper[idx1]; - e[idx0] = j; - diam[i - bbox.lower[idx2]] += concavity(s, e, idx1, grid, grid_size); - } - } - - for (std::size_t i = bbox.lower[idx2]; i <= bbox.upper[idx2]; i++) { - for (std::size_t j = bbox.lower[idx1]; j <= bbox.upper[idx1]; j++) { - Vec3_uint s, e; - s[idx0] = bbox.lower[idx0]; - s[idx2] = i; - s[idx1] = j; - e[idx0] = bbox.upper[idx0]; - e[idx2] = i; - e[idx1] = j; - diam2[i - bbox.lower[idx2]] += concavity(s, e, idx0, grid, grid_size); - } - } - - // Skip initial border - std::size_t border = (length / 10) + 0.5; - std::size_t pos1, end1 = length; - int grad = diam[0] - diam[1]; - for (pos1 = 2; pos1 < border; pos1++) { - int grad1 = diam[pos1 - 1] - diam[pos1]; - // Stop if the gradient flips or flattens significantly - if (!(grad * grad1 > 0 && grad1 > (grad>>1))) - break; - if (grad < grad1) - grad = grad1; - } - - grad = diam[length - 1] - diam[length - 2]; - for (end1 = length - 3; end1 > (length - border - 1); end1--) { - int grad1 = diam[end1 + 1] - diam[end1]; - // Stop if the gradient flips or flattens significantly - if (!(grad * grad1 > 0 && grad1 > (grad >> 1))) - break; - if (grad < grad1) - grad = grad1; - } - - std::size_t pos2, end2 = length; - grad = diam2[0] - diam2[1]; - for (pos2 = 2; pos2 < border; pos2++) { - int grad2 = diam[pos2 - 1] - diam[pos2]; - // Stop if the gradient flips or flattens significantly - if (!(grad * grad2 > 0 && grad2 > (grad >> 1))) - break; - if (grad < grad2) - grad = grad2; - } - - grad = diam2[length - 1] - diam2[length - 2]; - for (end2 = length - 3; end2 > (length - border - 1); end2--) { - int grad2 = diam[end2 + 1] - diam[end2]; - // Stop if the gradient flips or flattens significantly - if (!(grad * grad2 > 0 && grad2 > (grad >> 1))) - break; - if (grad < grad2) - grad = grad2; - } - - std::size_t conc1 = abs(diam[pos1 + 1] - diam[pos1]); - std::size_t conc2 = abs(diam2[pos2 + 1] - diam2[pos2]); - - for (std::size_t i = pos1;i conc1) { - pos1 = i - 1; - conc1 = abs(diam[i] - diam[i - 1]); - } - - for (std::size_t i = pos2; i < end2; i++) - if (unsigned(abs(diam2[i] - diam2[i - 1])) > conc2) { - pos2 = i - 1; - conc2 = abs(diam2[i] - diam2[i - 1]); - } - - - if (conc1 <= conc2) { - if (pos1 < 2 || (length - 3) < pos1) - location = (bbox.upper[axis] + bbox.lower[axis]) / 2; - else - location = pos1 + bbox.lower[axis]; - } - else { - if (pos2 < 2 || (length - 3) < pos2) - location = (bbox.upper[axis] + bbox.lower[axis]) / 2; - else - location = pos2 + bbox.lower[axis]; - } -} - -template -void choose_splitting_plane(Candidate& c, unsigned int &axis, unsigned int &location, std::vector& grid, const Vec3_uint& grid_size, const NamedParameters& np) { - const bool search_concavity = parameters::choose_parameter(parameters::get_parameter(np, internal_np::split_at_concavity), true); - const std::array span = {c.bbox.upper[0] - c.bbox.lower[0], c.bbox.upper[1] - c.bbox.lower[1], c.bbox.upper[2] - c.bbox.lower[2]}; - - // Split longest axis - axis = (span[0] >= span[1]) ? 0 : 1; - axis = (span[axis] >= span[2]) ? axis : 2; - - if (span[axis] >= 8 && search_concavity) - choose_splitting_location_by_concavity(axis, location, c.bbox, grid, grid_size); - else - location = (c.bbox.upper[axis] + c.bbox.lower[axis]) / 2; -} - -template -bool finished(Candidate &c, const NamedParameters& np) { - const typename GeomTraits::FT max_error = parameters::choose_parameter(parameters::get_parameter(np, internal_np::volume_error), 0.01); - - if (c.ch.volume_error <= max_error) - return true; - - std::size_t max_span = 0; - for (std::size_t i = 0;i<3;i++) { - const std::size_t span = c.bbox.upper[i] - c.bbox.lower[i]; - max_span = (std::max)(max_span, span); - } - - if (max_span <= 1) - return true; - - return false; -} - -template -void recurse(std::vector>& candidates, std::vector& grid, const Vec3_uint& grid_size, const Bbox_3& bbox, const typename GeomTraits::FT& voxel_size, const NamedParameters& np, CGAL::Parallel_tag) { -#ifdef CGAL_LINKED_WITH_TBB - const std::size_t max_depth = parameters::choose_parameter(parameters::get_parameter(np, internal_np::maximum_depth), 10); - - std::vector> final_candidates; - - std::size_t depth = 0; - - while (!candidates.empty() && depth < max_depth) { - depth++; - std::vector> former_candidates = std::move(candidates); - for (Candidate& c : former_candidates) { - - if (finished(c, np)) { - CGAL::Bbox_3 bbox = CGAL::bbox_3(c.ch.points.begin(), c.ch.points.end()); - bbox.scale(1.1); - c.ch.bbox = bbox; - final_candidates.push_back(std::move(c)); - continue; - } - unsigned int axis = 0, location = 0; - choose_splitting_plane(c, axis, location, grid, grid_size, np); - split(candidates, c, axis, location); - } - tbb::parallel_for_each(candidates, [&](Candidate& c) { - compute_candidate(c, bbox, voxel_size); - }); - } - - if (candidates.empty()) - std::swap(candidates, final_candidates); - else - std::move(final_candidates.begin(), final_candidates.end(), std::back_inserter(candidates)); - -#else - CGAL_USE(candidates); - CGAL_USE(grid); - CGAL_USE(grid_size); - CGAL_USE(bbox); - CGAL_USE(voxel_size); - CGAL_USE(np); -#endif -} - -template -void recurse(std::vector>& candidates, std::vector& grid, const Vec3_uint& grid_size, const Bbox_3& bbox, const typename GeomTraits::FT& voxel_size, const NamedParameters& np, CGAL::Sequential_tag) { - const std::size_t max_depth = parameters::choose_parameter(parameters::get_parameter(np, internal_np::maximum_depth), 10); - - std::vector> final_candidates; - - std::size_t depth = 0; - - while (!candidates.empty() && depth < max_depth) { - depth++; - std::vector> former_candidates = std::move(candidates); - for (Candidate& c : former_candidates) { - - if (finished(c, np)) { - CGAL::Bbox_3 bbox = CGAL::bbox_3(c.ch.points.begin(), c.ch.points.end()); - bbox.scale(1.1); - c.ch.bbox = bbox; - final_candidates.push_back(std::move(c)); - continue; - } - unsigned int axis = 0, location = 0; - choose_splitting_plane(c, axis, location, grid, grid_size, np); - split(candidates, c, axis, location); - } - - for (Candidate& c : candidates) - compute_candidate(c, bbox, voxel_size); - } - - if (candidates.empty()) - std::swap(candidates, final_candidates); - else - std::move(final_candidates.begin(), final_candidates.end(), std::back_inserter(candidates)); -} - -template -void merge(std::vector>& candidates, const typename GeomTraits::FT& hull_volume, const unsigned int max_convex_hulls, CGAL::Parallel_tag) { -#ifdef CGAL_LINKED_WITH_TBB - if (candidates.size() <= max_convex_hulls) - return; - - using FT = typename GeomTraits::FT; - using Point_3 = typename GeomTraits::Point_3; - - struct Merged_candidate { - std::size_t ch_a, ch_b; - int ch; - FT volume_error; - - bool operator < (const Merged_candidate& other) const { - return (volume_error > other.volume_error); - } - - Merged_candidate() : ch_a(-1), ch_b(-1) {} - Merged_candidate(std::size_t ch_a, std::size_t ch_b) : ch_a(ch_a), ch_b(ch_b) {} - }; - - tbb::concurrent_unordered_map> hulls; - std::atomic num_hulls = candidates.size(); - - std::unordered_set keep; - - for (std::size_t i = 0; i < candidates.size(); i++) { - hulls.emplace(i, std::move(candidates[i])); - keep.insert(i); - } - - candidates.clear(); - candidates.reserve(max_convex_hulls); - - std::vector todo; - tbb::concurrent_priority_queue queue; - - const auto do_merge = [hull_volume, &hulls, &num_hulls](Merged_candidate& m) { - Convex_hull_candidate& ci = hulls[m.ch_a]; - Convex_hull_candidate& cj = hulls[m.ch_b]; - m.ch = num_hulls.fetch_add(1); - Convex_hull_candidate& ch = hulls[m.ch]; - - ch.bbox = box_union(ci.bbox, cj.bbox); - std::vector pts(ci.points.begin(), ci.points.end()); - pts.reserve(pts.size() + cj.points.size()); - std::copy(cj.points.begin(), cj.points.end(), std::back_inserter(pts)); - convex_hull_3(pts.begin(), pts.end(), ch.points, ch.indices); - - ch.volume = volume(ch.points, ch.indices); - - ch.volume_error = m.volume_error = CGAL::abs(ci.volume + cj.volume - ch.volume) / hull_volume; - }; - - for (std::size_t i : keep) { - const Convex_hull_candidate& ci = hulls[i]; - for (std::size_t j : keep) { - if (j <= i) - continue; - const Convex_hull_candidate& cj = hulls[j]; - if (CGAL::do_intersect(ci.bbox, cj.bbox)) - todo.emplace_back(Merged_candidate(i, j)); - else { - Merged_candidate m(i, j); - Bbox_3 bbox = box_union(ci.bbox, cj.bbox).bbox(); - m.ch = -1; - m.volume_error = CGAL::abs(ci.volume + cj.volume - bbox.x_span() * bbox.y_span() * bbox.z_span()) / hull_volume; - queue.push(std::move(m)); - } - } - } - - // parallel for if available - tbb::parallel_for_each(todo, do_merge); - for (Merged_candidate& m : todo) - queue.push(std::move(m)); - todo.clear(); - - while (!queue.empty() && keep.size() > max_convex_hulls) { - Merged_candidate m; - while (!queue.try_pop(m) && !queue.empty()); - - auto ch_a = hulls.find(m.ch_a); - if (ch_a == hulls.end()) - continue; - - auto ch_b = hulls.find(m.ch_b); - if (ch_b == hulls.end()) - continue; - - if (m.ch == -1) - do_merge(m); - - keep.erase(m.ch_a); - keep.erase(m.ch_b); - - hulls.unsafe_erase(ch_a); - hulls.unsafe_erase(ch_b); - - const Convex_hull_candidate& cj = hulls[m.ch]; - - for (std::size_t id : keep) { - const Convex_hull_candidate& ci = hulls[id]; - if (CGAL::do_intersect(ci.bbox, cj.bbox)) - todo.emplace_back(Merged_candidate(id, m.ch)); - else { - Merged_candidate merged(id, m.ch); - Bbox_3 bbox = box_union(ci.bbox, cj.bbox).bbox(); - merged.ch = -1; - merged.volume_error = CGAL::abs(ci.volume + cj.volume - bbox.x_span() * bbox.y_span() * bbox.z_span()) / hull_volume; - queue.push(std::move(merged)); - } - } - - keep.insert(m.ch); - - tbb::parallel_for_each(todo, do_merge); - for (Merged_candidate& m : todo) - queue.push(std::move(m)); - todo.clear(); - } - - num_hulls = 0; - - for (std::size_t i : keep) - candidates.push_back(std::move(hulls[i])); -#else - CGAL_USE(candidates); - CGAL_USE(hull_volume); - assert(false); -#endif -} - -template -void merge(std::vector>& candidates, const typename GeomTraits::FT& hull_volume, const unsigned int max_convex_hulls, CGAL::Sequential_tag) { - if (candidates.size() <= max_convex_hulls) - return; - - using FT = typename GeomTraits::FT; - using Point_3 = typename GeomTraits::Point_3; - - struct Merged_candidate { - std::size_t ch_a, ch_b; - int ch; - FT volume_error; - - bool operator < (const Merged_candidate& other) const { - return (volume_error > other.volume_error); - } - - Merged_candidate() : ch_a(-1), ch_b(-1) {} - Merged_candidate(std::size_t ch_a, std::size_t ch_b) : ch_a(ch_a), ch_b(ch_b) {} - }; - - std::unordered_map> hulls; - std::size_t num_hulls = candidates.size(); - - std::unordered_set keep; - - for (std::size_t i = 0; i < candidates.size(); i++) { - hulls.emplace(i, std::move(candidates[i])); - keep.insert(i); - } - - candidates.clear(); - candidates.reserve(max_convex_hulls); - - std::priority_queue queue; - - const auto do_merge = [hull_volume, &hulls, &num_hulls](Merged_candidate& m) { - Convex_hull_candidate& ci = hulls[m.ch_a]; - Convex_hull_candidate& cj = hulls[m.ch_b]; - - m.ch = num_hulls++; - Convex_hull_candidate& ch = hulls[m.ch]; - - ch.bbox = box_union(ci.bbox, cj.bbox); - std::vector pts(ci.points.begin(), ci.points.end()); - pts.reserve(pts.size() + cj.points.size()); - std::copy(cj.points.begin(), cj.points.end(), std::back_inserter(pts)); - convex_hull_3(pts.begin(), pts.end(), ch.points, ch.indices); - - ch.volume = volume(ch.points, ch.indices); - - ch.volume_error = m.volume_error = CGAL::abs(ci.volume + cj.volume - ch.volume) / hull_volume; - }; - - for (std::size_t i : keep) { - const Convex_hull_candidate& ci = hulls[i]; - for (std::size_t j : keep) { - if (j <= i) - continue; - const Convex_hull_candidate& cj = hulls[j]; - if (CGAL::do_intersect(ci.bbox, cj.bbox)) { - Merged_candidate m(i, j); - - m.ch = num_hulls++; - Convex_hull_candidate& ch = hulls[m.ch]; - ch.bbox = box_union(ci.bbox, cj.bbox); - std::vector pts(ci.points.begin(), ci.points.end()); - pts.reserve(pts.size() + cj.points.size()); - std::copy(cj.points.begin(), cj.points.end(), std::back_inserter(pts)); - convex_hull_3(pts.begin(), pts.end(), ch.points, ch.indices); - - ch.volume = volume(ch.points, ch.indices); - - ch.volume_error = m.volume_error = CGAL::abs(ci.volume + cj.volume - ch.volume) / hull_volume; - queue.push(std::move(m)); - } - else { - Merged_candidate m(i, j); - Bbox_3 bbox = box_union(ci.bbox, cj.bbox).bbox(); - m.ch = -1; - m.volume_error = CGAL::abs(ci.volume + cj.volume - bbox.x_span() * bbox.y_span() * bbox.z_span()) / hull_volume; - queue.push(std::move(m)); - } - } - } - - while (!queue.empty() && keep.size() > max_convex_hulls) { - Merged_candidate m = queue.top(); - queue.pop(); - - auto ch_a = hulls.find(m.ch_a); - if (ch_a == hulls.end()) - continue; - - auto ch_b = hulls.find(m.ch_b); - if (ch_b == hulls.end()) - continue; - - if (m.ch == -1) - do_merge(m); - - keep.erase(m.ch_a); - keep.erase(m.ch_b); - - hulls.erase(ch_a); - hulls.erase(ch_b); - - const Convex_hull_candidate& cj = hulls[m.ch]; - - for (std::size_t id : keep) { - const Convex_hull_candidate& ci = hulls[id]; - if (CGAL::do_intersect(ci.bbox, cj.bbox)) { - Merged_candidate merged(id, m.ch); - - merged.ch = num_hulls++; - Convex_hull_candidate& ch = hulls[merged.ch]; - ch.bbox = box_union(ci.bbox, cj.bbox); - std::vector pts(ci.points.begin(), ci.points.end()); - pts.reserve(pts.size() + cj.points.size()); - std::copy(cj.points.begin(), cj.points.end(), std::back_inserter(pts)); - convex_hull_3(pts.begin(), pts.end(), ch.points, ch.indices); - - ch.volume = volume(ch.points, ch.indices); - - ch.volume_error = merged.volume_error = CGAL::abs(ci.volume + cj.volume - ch.volume) / hull_volume; - queue.push(std::move(merged)); - } - else { - Merged_candidate merged(id, m.ch); - Bbox_3 bbox = box_union(ci.bbox, cj.bbox).bbox(); - merged.ch = -1; - merged.volume_error = CGAL::abs(ci.volume + cj.volume - bbox.x_span() * bbox.y_span() * bbox.z_span()) / hull_volume; - queue.push(std::move(merged)); - } - } - - keep.insert(m.ch); - } - - num_hulls = 0; - - for (std::size_t i : keep) - candidates.push_back(std::move(hulls[i])); -} - -} // namespace internal - -/** - * \ingroup PMP_convex_decomposition_grp - * - * \brief approximates the input mesh by a number of convex hulls. The input mesh is voxelized and the voxels intersecting with the mesh are labeled as surface. - * The remaining voxels are labeled as outside or inside by flood fill, in case the input mesh is closed, or by axis-aligned ray shooting, i.e., along x, - * y and z-axis in positive and negative directions. A voxel is only labeled as inside if at least 3 faces facing away from the voxel have been hit and - * no face facing towards the voxel. In a next step, the convex hull of the mesh is hierarchically split until the `volume_error` threshold is satisfied. - * Afterwards, a greedy pair-wise merging combines smaller convex hulls until the given number of convex hulls is met. - * - * \tparam FaceGraph a model of `HalfedgeListGraph`, and `FaceListGraph` - * - * \tparam OutputIterator must be an output iterator accepting variables of type `std::pair, std::vector > >`. - * - * \tparam NamedParameters a sequence of \ref bgl_namedparameters "Named Parameters" - * - * \param tmesh the input triangle mesh to approximate by convex hulls - * \param out_hulls output iterator into which convex hulls are recorded - * \param np an optional sequence of \ref bgl_namedparameters "Named Parameters" among the ones listed below - * - * \cgalNamedParamsBegin - * - * \cgalParamNBegin{maximum_number_of_voxels} - * \cgalParamDescription{gives an upper bound of the number of voxels. The longest bounding box side will have a length of the cubic root of `maximum_number_of_voxels` rounded down. Cannot be smaller than 64. } - * \cgalParamType{unsigned int} - * \cgalParamDefault{1,000,000} - * \cgalParamNEnd - * - * \cgalParamNBegin{maximum_depth} - * \cgalParamDescription{maximum depth of hierarchical splits} - * \cgalParamType{unsigned int} - * \cgalParamDefault{10} - * \cgalParamNEnd - * - * \cgalParamNBegin{maximum_number_of_convex_hulls} - * \cgalParamDescription{maximum number of convex hulls produced by the method} - * \cgalParamType{unsigned int} - * \cgalParamDefault{16} - * \cgalParamNEnd - * - * \cgalParamNBegin{volume_error} - * \cgalParamDescription{maximum difference in fraction of volumes of the local convex hull with the sum of voxel volumes. If surpassed, the convex hull will be split if the `maximum_depth` has not been reached yet.} - * \cgalParamType{double} - * \cgalParamDefault{0.01} - * \cgalParamNEnd - * - * \cgalParamNBegin{split_at_concavity} - * \cgalParamDescription{If `true`, the local box of a convex hull is split at the concavity along the longest axis of the bounding box. Otherwise, it is split in the middle of the longest axis, which is faster, but less precise.} - * \cgalParamType{Boolean} - * \cgalParamDefault{true} - * \cgalParamNEnd - * - * \cgalParamNBegin{vertex_point_map} - * \cgalParamDescription{a property map associating points to the vertices of `mesh`} - * \cgalParamType{a class model of `ReadablePropertyMap` with `boost::graph_traits::%vertex_descriptor` - * as key type and `%Point_3` as value type} - * \cgalParamDefault{`boost::get(CGAL::vertex_point, mesh)`} - * \cgalParamExtra{If this parameter is omitted, an internal property map for `CGAL::vertex_point_t` - * must be available in `FaceGraph`.} - * \cgalParamNEnd - * - * \cgalParamNBegin{concurrency_tag} - * \cgalParamDescription{a tag indicating if the task should be done using one or several threads.} - * \cgalParamType{Either `CGAL::Sequential_tag`, or `CGAL::Parallel_tag`, or `CGAL::Parallel_if_available_tag`} - * \cgalParamDefault{`CGAL::Parallel_if_available_tag`} - * \cgalParamNEnd - * - * \cgalParamNBegin{geom_traits} - * \cgalParamDescription{an instance of a geometric traits class} - * \cgalParamType{a class model of `Kernel`} - * \cgalParamDefault{a \cgal Kernel deduced from the point type of `FaceGraph`, using `CGAL::Kernel_traits`} - * \cgalParamExtra{The geometric traits class must be compatible with the vertex point type.} - * \cgalParamNEnd - * - * \cgalNamedParamsEnd - * - * \return the number of convex hulls. Note that this value may be lower than the `maximum_number_of_convex_hulls`, for example if the specified `volume_error` is quickly met. - * - * \sa `CGAL::Polygon_mesh_processing::polygon_soup_to_polygon_mesh()` - */ -template -std::size_t approximate_convex_decomposition(const FaceGraph& tmesh, OutputIterator out_hulls, const NamedParameters& np = parameters::default_values()) { - using Geom_traits = typename GetGeomTraits::type; - - using FT = typename Geom_traits::FT; - const unsigned int num_voxels = parameters::choose_parameter(parameters::get_parameter(np, internal_np::maximum_number_of_voxels), 1000000); - using Concurrency_tag = typename internal_np::Lookup_named_param_def::type; - -#ifndef CGAL_LINKED_WITH_TBB - if constexpr (std::is_same_v) { - CGAL_error_msg("CGAL was not compiled with TBB support. Use Sequential_tag instead."); - return 0; - } -#endif - - const unsigned int max_convex_hulls = parameters::choose_parameter(parameters::get_parameter(np, internal_np::maximum_number_of_convex_hulls), 16); - assert(max_convex_hulls > 0); - - if (max_convex_hulls == 1) { - internal::Convex_hull_candidate ch; - - using parameters::choose_parameter; - using parameters::get_parameter; - using VPM = typename GetVertexPointMap::const_type; - typedef CGAL::Property_map_to_unary_function Vpmap_fct; - - VPM vpm = choose_parameter(get_parameter(np, internal_np::vertex_point), get_const_property_map(CGAL::vertex_point, tmesh)); - Vpmap_fct v2p(vpm); - - convex_hull_3(boost::make_transform_iterator(vertices(tmesh).begin(), v2p), boost::make_transform_iterator(vertices(tmesh).end(), v2p), ch.points, ch.indices); - - *out_hulls = std::make_pair(std::move(ch.points), std::move(ch.indices)); - return 1; - } - - Bbox_3 bb = bbox(tmesh); - const auto [grid_size, voxel_size] = internal::calculate_grid_size(bb, num_voxels); - - std::vector grid(grid_size[0] * grid_size[1] * grid_size[2], internal::Grid_cell::INSIDE); - - std::vector> candidates(1); - init(candidates[0], tmesh, grid, bb, grid_size, voxel_size, Concurrency_tag()); - - const FT hull_volume = candidates[0].ch.volume; - - recurse(candidates, grid, grid_size, bb, voxel_size, np, Concurrency_tag()); - - std::vector> hulls; - for (internal::Candidate &c : candidates) - hulls.emplace_back(std::move(c.ch)); - - candidates.clear(); - - // merge until target number is reached - merge(hulls, hull_volume, max_convex_hulls, Concurrency_tag()); - - for (std::size_t i = 0; i < hulls.size(); i++) - *out_hulls++ = std::make_pair(std::move(hulls[i].points), std::move(hulls[i].indices)); - - return hulls.size(); -} - -} // namespace Polygon_mesh_processing -} // namespace CGAL - -#endif // CGAL_POLYGON_MESH_PROCESSING_CONVEX_DECOMPOSITION_H diff --git a/STL_Extension/include/CGAL/STL_Extension/internal/parameters_interface.h b/STL_Extension/include/CGAL/STL_Extension/internal/parameters_interface.h index 6424c0e4448..f02f1fd4676 100644 --- a/STL_Extension/include/CGAL/STL_Extension/internal/parameters_interface.h +++ b/STL_Extension/include/CGAL/STL_Extension/internal/parameters_interface.h @@ -182,7 +182,7 @@ CGAL_add_named_parameter(snap_grid_size_t, snap_grid_size, snap_grid_size) CGAL_add_named_parameter(maximum_number_of_voxels_t, maximum_number_of_voxels, maximum_number_of_voxels) CGAL_add_named_parameter(maximum_depth_t, maximum_depth, maximum_depth) CGAL_add_named_parameter(volume_error_t, volume_error, volume_error) -CGAL_add_named_parameter(maximum_number_of_convex_hulls_t, maximum_number_of_convex_hulls, maximum_number_of_convex_hulls) +CGAL_add_named_parameter(maximum_number_of_convex_volumes_t, maximum_number_of_convex_volumes, maximum_number_of_convex_volumes) CGAL_add_named_parameter(split_at_concavity_t, split_at_concavity, split_at_concavity) diff --git a/Surface_mesh_decomposition/examples/Surface_mesh_decomposition/approximate_convex_decomposition.cpp b/Surface_mesh_decomposition/examples/Surface_mesh_decomposition/approximate_convex_decomposition.cpp index c7a002ad401..ea2518d1718 100644 --- a/Surface_mesh_decomposition/examples/Surface_mesh_decomposition/approximate_convex_decomposition.cpp +++ b/Surface_mesh_decomposition/examples/Surface_mesh_decomposition/approximate_convex_decomposition.cpp @@ -27,19 +27,19 @@ int main(int argc, char* argv[]) return 1; } - std::vector convex_hulls; - convex_hulls.reserve(9); + std::vector convex_volumes; + convex_volumes.reserve(9); - CGAL::approximate_convex_decomposition(mesh, std::back_inserter(convex_hulls), + CGAL::approximate_convex_decomposition(mesh, std::back_inserter(convex_volumes), CGAL::parameters::maximum_depth(10) .volume_error(0.1) - .maximum_number_of_convex_hulls(9) + .maximum_number_of_convex_volumes(9) .split_at_concavity(true) .maximum_number_of_voxels(1000000) .concurrency_tag(CGAL::Parallel_if_available_tag())); - for (std::size_t i = 0;i>& candidates, const typ * * \cgalNamedParamsEnd * - * \return the number of convex hulls. Note that this value may be lower than the `maximum_number_of_convex_hulls`, for example if the specified `volume_error` is quickly met. + * \return the number of convex hulls. Note that this value may be lower than the `maximum_number_of_convex_volumes`, for example if the specified `volume_error` is quickly met. * * \sa `CGAL::Polygon_mesh_processing::polygon_soup_to_polygon_mesh()` */ @@ -1626,7 +1626,7 @@ std::size_t approximate_convex_decomposition(const FaceGraph& tmesh, OutputItera } #endif - const unsigned int max_convex_hulls = parameters::choose_parameter(parameters::get_parameter(np, internal_np::maximum_number_of_convex_hulls), 16); + const unsigned int max_convex_hulls = parameters::choose_parameter(parameters::get_parameter(np, internal_np::maximum_number_of_convex_volumes), 16); assert(max_convex_hulls > 0); if (max_convex_hulls == 1) { diff --git a/Surface_mesh_decomposition/test/Surface_mesh_decomposition/test_acd.cpp b/Surface_mesh_decomposition/test/Surface_mesh_decomposition/test_acd.cpp index 93100cf5cc5..9bea6abac6e 100644 --- a/Surface_mesh_decomposition/test/Surface_mesh_decomposition/test_acd.cpp +++ b/Surface_mesh_decomposition/test/Surface_mesh_decomposition/test_acd.cpp @@ -23,27 +23,27 @@ int main(int argc, char* argv[]) Mesh mesh; - std::vector convex_hulls; + std::vector convex_volumes; // Try with empty mesh - CGAL::approximate_convex_decomposition(mesh, std::back_inserter(convex_hulls), + CGAL::approximate_convex_decomposition(mesh, std::back_inserter(convex_volumes), CGAL::parameters::maximum_depth(10) .volume_error(0.1) - .maximum_number_of_convex_hulls(9) + .maximum_number_of_convex_volumes(9) .split_at_concavity(true) .maximum_number_of_voxels(1000000) .concurrency_tag(CGAL::Parallel_if_available_tag())); - assert(convex_hulls.size() == 0); + assert(convex_volumes.size() == 0); if (!PMP::IO::read_polygon_mesh(filename, mesh)) { std::cerr << "Invalid input." << std::endl; return 1; } - convex_hulls.reserve(9); + convex_volumes.reserve(9); - CGAL::approximate_convex_decomposition(mesh, std::back_inserter(convex_hulls), + CGAL::approximate_convex_decomposition(mesh, std::back_inserter(convex_volumes), CGAL::parameters::maximum_depth(10) .volume_error(0.1) .maximum_number_of_convex_hulls(9) @@ -51,8 +51,8 @@ int main(int argc, char* argv[]) .maximum_number_of_voxels(1000000) .concurrency_tag(CGAL::Parallel_if_available_tag())); - assert(convex_hulls.size() > 0); - assert(convex_hulls.size() <= 9); + assert(convex_volumes.size() > 0); + assert(convex_volumes.size() <= 9); return 0; }