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 3a36080cd66..bfd47c38e97 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 @@ -17,8 +17,20 @@ decomposing both polyhedra into convex pieces, compute pairwise Minkowski sums of the convex pieces, and unite the pairwise sums. While it is desirable to have a decomposition into a minimum number of -pieces, this problem is known to be NP-hard \cgalCite{c-cpplb-84}. Our -implementation decomposes a Nef polyhedron \f$ N\f$ into \cgalBigO{r^2} convex +pieces, this problem is known to be NP-hard \cgalCite{c-cpplb-84}. This +package offers two methods for decomposing polyhedra. The +\ref Convex_decomposition_3Nef "Convex Decomposition of Nef Polyhedra" +splits polyhedra into convex pieces with an upper bound on the number +of pieces. The \ref Convex_decomposition_3ACD_Intro +"Approximate convex decomposition" method offers a fast +approximate decomposition of the convex hull into convex volumes. While +any number of convex volumes can be generated, these convex volumes +are more compact than the convex hull, but still include additional +empty space than just the input polyhedron. + +\section Convex_decomposition_3Nef Convex Decomposition of Nef Polyhedra + +Our implementation decomposes a Nef polyhedron \f$ N\f$ into \cgalBigO{r^2} convex pieces, where \f$ r\f$ is the number of edges that have two adjacent facets that span an angle of more than 180 degrees with respect to the interior of the polyhedron. Those edges are also called reflex edges. @@ -38,7 +50,7 @@ illustrates the two steps. At the moment our implementation is restricted to the decomposition of bounded polyhedra. An extension to unbounded polyhedra is planned. -\section Convex_decomposition_3InterfaceandUsage Interface and Usage +\subsection Convex_decomposition_3InterfaceandUsage Interface and Usage An instance of `Nef_polyhedron_3` represents a subdivision of the three-dimensional space into vertices, edges, facets, and @@ -79,6 +91,271 @@ support the use of extended kernels in the convex decomposition. \cgalExample{Convex_decomposition_3/list_of_convex_parts.cpp} +\section Convex_decomposition_3ACD_Intro Approximate Convex Decomposition + +\cgalFigureAnchor{Acd_topfig} +
+ +
+\cgalFigureCaptionBegin{Acd_topfig} +Approximate convex decomposition of the elephant.off model.\n From left to right: 1. input mesh 2. 5 convex volumes 3. 8 convex volumes 4. 12 convex volumes +\cgalFigureCaptionEnd + +The H-VACD, introduced by computes a set of convex volumes that fit the polyhedron. Contrary to the decomposition of the polyhedron into convex parts, the convex volumes cover the polyhedron, but also include additional volume outside of the polyhedron. +A sufficiently tight enclosure of the polyhedron by several convex volumes allows fast intersection calculation and collision detection among polyhedra while offering a non-hierachical approach and may thus be easier to use in a parallel setting. +The resulting set of convex volumes minimizes the volume between their union and the polyhedron while fully including the input polyhedron. While the optimal solution with `n` convex hulls that cover the polydron with the smallest additional volume remains NP-hard, this method provides a fast error-driven approximation. + +\subsection Convex_decomposition_3ACD_Algorithm Algorithm +The algorithm computes a set of convex volumes \f$ C=\{C_i\f$ with \f$ i \in[0..n-1]\} \f$ that cover the input polyhedron while minimizing the additional covered volume: + +\f{equation}{ +\arg \min_C d(\bigcup_{C_i \in C} C_i, P) \\ +\f} +
with
+\f{equation}{ +d(A, B) = |A| - |B| +\f} + +Where \f$|A|\f$ is the volume of A, P is the input polyhedron and \f$C_i\f$ are convex volumes. The convex volumes \f$C_i\f$ are pairwise disjunct, i.e., \f$|C_i \cap C_j| = 0\f$ if \f$i \neq j\f$. And the union of convex volumes contain the input polyhedron \f$ P \subset \bigcup_{C_i \in C} \f$. + +\cgalFigureAnchor{Acd_pipelinefig} +
+ +
+\cgalFigureCaptionBegin{Acd_pipelinefig} +Approximate convex decomposition pipeline.\n From left to right: 1. input mesh 2. voxel grid 3. convex volumes after error-driven splitting 4. final convex volumes after merging +\cgalFigureCaptionEnd + +The method employs a top-down splitting phase followed by a bottom-up merging to achieve the target number of convex volumes. The splitting phase aims at decomposing the input mesh into smaller mostly convex parts. Each part of the input mesh is approximated with its convex hull. In a hierarchical manner, each part of the mesh is split into two parts when its convexity is low. The convexity is measured by the volume difference of the part and its convex hull. Splitting a part into two can be done by simply cutting the longest side of the bounding box in half. A better choice is often found by searching the longest side of the bounding box for a concave spot. However, it comes with a higher computational cost. The hierarchical splitting stops when either the convexity is sufficiently high or the maximum depth is reached. +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. + + +\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_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. +- `split_at_concavity`: The splitting can be either performed after searching a concavity on the longest side of the bounding box or simply by splitting the longest side of the bounding box in half. The default value is true, i.e., splitting at the concavity. + +\subsection Convex_decomposition_3ACD_Performance Performance + +Here will be images and more tables to show the impact of different parameters + +The method has been evaluated on several models: + + + + + + + + + +
+
+Data set + +Faces + +Volume + +Convex hull volume + +Overhead +
+
+Camel + +19.536 + +0.0468 + +0.15541 + +2.32388 +
+Elephant + +5.558 + +0.0462 + +0.12987 + +1.81087 +
+Triceratops + +5.660 + +136.732 + +336.925 + +1.46412 +
+
+ +If not mentioned otherwise, all tests used a volume error of 0.01, a maximum depth of 10, 1 million voxels and split at the concavity. + +Impact of varying the number of generated convex volumes with splitting at the concavity on volume overhead: + + + + + + + + + +
+
+Data set + +Split location + +5 volumes + +7 volumes + +9 volumes + +11 volumes + +12 volumes +
+
+Camel + +Concavity + +0.8006 + +0.6680 + +0.5871 + +0.5736 + +0.5463 + +
+Elephant + +Concavity + +1.1927 + +0.9731 + +0.84613 + +0.7506 + +0.6947 + +
+Triceratops + +Concavity + +0.9770 + +0.7676 + +0.6722 + +0.5971 + +0.5658 + +
+
+ +And by using the mid split: + + + + + + + + + +
+
+Data set + +Split location + +5 volumes + +7 volumes + +9 volumes + +11 volumes + +12 volumes +
+
+Camel + +Mid + +1.1158 + +1.0468 + +0.8660 + +0.6764 + +0.6057 + +
+Elephant + +Mid + +1.2121 + +1.00867 + +0.8444 + +0.7465 + +0.6931 + +
+Triceratops + +Mid + +1.2191 + +0.6870 + +0.8463 + +0.7375 + +0.6871 + +
+
+ +The running time for all cases in the above tables were between 1.4 and 3 seconds while being slightly lower when splitting at the concavity. Although searching the voxel grid for the concavity takes additional computational time, it is more than compensated by fewer splits. + +\subsection Convex_decomposition_3ACD_Example Example + +\cgalExample{Surface_mesh_decomposition/approximate_convex_decomposition.cpp } + */ } /* namespace CGAL */ diff --git a/Convex_decomposition_3/doc/Convex_decomposition_3/Doxyfile.in b/Convex_decomposition_3/doc/Convex_decomposition_3/Doxyfile.in index aa7cdb5e1b8..ffc85d262bd 100644 --- a/Convex_decomposition_3/doc/Convex_decomposition_3/Doxyfile.in +++ b/Convex_decomposition_3/doc/Convex_decomposition_3/Doxyfile.in @@ -3,3 +3,7 @@ PROJECT_NAME = "CGAL ${CGAL_DOC_VERSION} - Convex Decomposition of Polyhedra" EXTRACT_ALL = false HIDE_UNDOC_CLASSES = true + +INPUT += ${CGAL_Surface_mesh_decomposition_INCLUDE_DIR}/CGAL/approximate_convex_decomposition.h + +EXAMPLE_PATH += ${CGAL_Surface_mesh_decomposition_EXAMPLE_DIR} \ No newline at end of file diff --git a/Convex_decomposition_3/doc/Convex_decomposition_3/PackageDescription.txt b/Convex_decomposition_3/doc/Convex_decomposition_3/PackageDescription.txt index fac4f8479c7..4db19541a8e 100644 --- a/Convex_decomposition_3/doc/Convex_decomposition_3/PackageDescription.txt +++ b/Convex_decomposition_3/doc/Convex_decomposition_3/PackageDescription.txt @@ -5,8 +5,9 @@ \cgalPkgDescriptionBegin{Convex Decomposition of Polyhedra,PkgConvexDecomposition3} \cgalPkgPicture{Convex_decomposition_3/fig/Convex_decomposition_3-teaser.png} \cgalPkgSummaryBegin -\cgalPkgAuthor{Peter Hachenberger} -\cgalPkgDesc{This packages provides a function for decomposing a bounded polyhedron into convex sub-polyhedra. The decomposition yields \cgalBigO{r^2} convex pieces, where \f$ r\f$ is the number of edges, whose adjacent facets form an angle of more than 180 degrees with respect to the polyhedron's interior. This bound is worst-case optimal. } +\cgalPkgAuthor{Peter Hachenberger, Sven Oesau} +\cgalPkgDesc{This packages provides two functions for computing a set of convex volumes that cover a bounded polyhedron. The `CGAL::convex_decomposition_3()` splits provides a decomposition into \cgalBigO{r^2} convex pieces, where \f$ r\f$ is the number of edges, whose adjacent facets form an angle of more than 180 degrees with respect to the polyhedron's interior. This bound is worst-case optimal. +The `CGAL::approximate_convex_decomposition()` instead splits the convex hull of the polyhedron into a set of convex volumes. While these convex volumes cover additional space outside of the polyhedron, the computation is fast for any chosen number of convex volumes.} \cgalPkgManuals{Chapter_Convex_Decomposition_of_Polyhedra,PkgConvexDecomposition3Ref} \cgalPkgSummaryEnd \cgalPkgShortInfoBegin @@ -21,6 +22,7 @@ \cgalCRPSection{Functions} - `CGAL::convex_decomposition_3(Nef_polyhedron_3& N)` +- `CGAL::approximate_convex_decomposition()` */ diff --git a/Convex_decomposition_3/doc/Convex_decomposition_3/dependencies b/Convex_decomposition_3/doc/Convex_decomposition_3/dependencies index 94e7ccf5922..6d540982926 100644 --- a/Convex_decomposition_3/doc/Convex_decomposition_3/dependencies +++ b/Convex_decomposition_3/doc/Convex_decomposition_3/dependencies @@ -5,3 +5,4 @@ Algebraic_foundations Circulator Stream_support Nef_3 +BGL diff --git a/Convex_decomposition_3/doc/Convex_decomposition_3/examples.txt b/Convex_decomposition_3/doc/Convex_decomposition_3/examples.txt index b40169d2c81..ab8e7489906 100644 --- a/Convex_decomposition_3/doc/Convex_decomposition_3/examples.txt +++ b/Convex_decomposition_3/doc/Convex_decomposition_3/examples.txt @@ -1,3 +1,4 @@ /*! \example Convex_decomposition_3/list_of_convex_parts.cpp +\example Surface_mesh_decomposition/approximate_convex_decomposition.cpp */ diff --git a/Installation/include/CGAL/license/Surface_mesh_decomposition.h b/Installation/include/CGAL/license/Surface_mesh_decomposition.h new file mode 100644 index 00000000000..ee8697fe7df --- /dev/null +++ b/Installation/include/CGAL/license/Surface_mesh_decomposition.h @@ -0,0 +1,54 @@ +// Copyright (c) 2016 GeometryFactory SARL (France). +// All rights reserved. +// +// This file is part of CGAL (www.cgal.org) +// +// $URL$ +// $Id$ +// SPDX-License-Identifier: LGPL-3.0-or-later OR LicenseRef-Commercial +// +// Author(s) : Andreas Fabri +// +// Warning: this file is generated, see include/CGAL/license/README.md + +#ifndef CGAL_LICENSE_SURFACE_MESH_DECOMPOSITION_H +#define CGAL_LICENSE_SURFACE_MESH_DECOMPOSITION_H + +#include +#include + +#ifdef CGAL_SURFACE_MESH_DECOMPOSITION_COMMERCIAL_LICENSE + +# if CGAL_SURFACE_MESH_DECOMPOSITION_COMMERCIAL_LICENSE < CGAL_RELEASE_DATE + +# if defined(CGAL_LICENSE_WARNING) + + CGAL_pragma_warning("Your commercial license for CGAL does not cover " + "this release of the Triangulated Surface Mesh Decomposition package.") +# endif + +# ifdef CGAL_LICENSE_ERROR +# error "Your commercial license for CGAL does not cover this release \ + of the Triangulated Surface Mesh Decomposition package. \ + You get this error, as you defined CGAL_LICENSE_ERROR." +# endif // CGAL_LICENSE_ERROR + +# endif // CGAL_SURFACE_MESH_DECOMPOSITION_COMMERCIAL_LICENSE < CGAL_RELEASE_DATE + +#else // no CGAL_SURFACE_MESH_DECOMPOSITION_COMMERCIAL_LICENSE + +# if defined(CGAL_LICENSE_WARNING) + CGAL_pragma_warning("\nThe macro CGAL_SURFACE_MESH_DECOMPOSITION_COMMERCIAL_LICENSE is not defined." + "\nYou use the CGAL Triangulated Surface Mesh Decomposition package under " + "the terms of the GPLv3+.") +# endif // CGAL_LICENSE_WARNING + +# ifdef CGAL_LICENSE_ERROR +# error "The macro CGAL_SURFACE_MESH_DECOMPOSITION_COMMERCIAL_LICENSE is not defined.\ + You use the CGAL Triangulated Surface Mesh Decomposition package under the terms of \ + the GPLv3+. You get this error, as you defined CGAL_LICENSE_ERROR." +# endif // CGAL_LICENSE_ERROR + +#endif // no CGAL_SURFACE_MESH_DECOMPOSITION_COMMERCIAL_LICENSE + +#endif // CGAL_LICENSE_SURFACE_MESH_DECOMPOSITION_H diff --git a/Installation/include/CGAL/license/gpl_package_list.txt b/Installation/include/CGAL/license/gpl_package_list.txt index 791271f4ac9..a5e73be196b 100644 --- a/Installation/include/CGAL/license/gpl_package_list.txt +++ b/Installation/include/CGAL/license/gpl_package_list.txt @@ -93,6 +93,7 @@ Straight_skeleton_extrusion_2 2D Straight Skeleton Extrusion Stream_lines_2 2D Placement of Streamlines Surface_mesh_approximation Triangulated Surface Mesh Approximation Surface_mesh_deformation Triangulated Surface Mesh Deformation +Surface_mesh_decomposition Triangulated Surface Mesh Decomposition Surface_mesher 3D Surface Mesh Generation Surface_mesh Surface Mesh Surface_mesh_parameterization Triangulated Surface Mesh Parameterization diff --git a/Surface_mesh_decomposition/examples/Surface_mesh_decomposition/CMakeLists.txt b/Surface_mesh_decomposition/examples/Surface_mesh_decomposition/CMakeLists.txt new file mode 100644 index 00000000000..5e6d6cd5d2c --- /dev/null +++ b/Surface_mesh_decomposition/examples/Surface_mesh_decomposition/CMakeLists.txt @@ -0,0 +1,19 @@ +# Created by the script cgal_create_CMakeLists +# This is the CMake script for compiling a set of CGAL applications. + +cmake_minimum_required(VERSION 3.12...3.31) +project(Surface_mesh_decomposition_Examples) + +# CGAL and its components +find_package(CGAL REQUIRED) + +find_package(TBB QUIET) +include(CGAL_TBB_support) + +create_single_source_cgal_program("approximate_convex_decomposition.cpp") + +if(TARGET CGAL::TBB_support) + target_link_libraries(approximate_convex_decomposition PRIVATE CGAL::TBB_support) +else() + message(STATUS "NOTICE: Intel TBB was not found. Sequential code will be used.") +endif() 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 new file mode 100644 index 00000000000..c7a002ad401 --- /dev/null +++ b/Surface_mesh_decomposition/examples/Surface_mesh_decomposition/approximate_convex_decomposition.cpp @@ -0,0 +1,47 @@ +#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); + + CGAL::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 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(); +} + +template +void export_grid_voxels(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 * 8 << "\n" << + "property double x\n" << + "property double y\n" << + "property double z\n" << + "element face " << count * 6 << "\n" << + "property list uchar int vertex_indices\n" << + "end_header\n"; + + // export vertices + 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)*voxel_size) << " " << (bb.ymin() + (y)*voxel_size) << " " << (bb.zmin() + (z)*voxel_size) << "\n"; + stream << (bb.xmin() + (x + 1.0) * voxel_size) << " " << (bb.ymin() + (y)*voxel_size) << " " << (bb.zmin() + (z)*voxel_size) << "\n"; + stream << (bb.xmin() + (x)*voxel_size) << " " << (bb.ymin() + (y + 1.0) * voxel_size) << " " << (bb.zmin() + (z)*voxel_size) << "\n"; + stream << (bb.xmin() + (x + 1.0) * voxel_size) << " " << (bb.ymin() + (y + 1.0) * voxel_size) << " " << (bb.zmin() + (z)*voxel_size) << "\n"; + stream << (bb.xmin() + (x)*voxel_size) << " " << (bb.ymin() + (y)*voxel_size) << " " << (bb.zmin() + (z + 1.0)*voxel_size) << "\n"; + stream << (bb.xmin() + (x + 1.0) * voxel_size) << " " << (bb.ymin() + (y)*voxel_size) << " " << (bb.zmin() + (z + 1.0)*voxel_size) << "\n"; + stream << (bb.xmin() + (x)*voxel_size) << " " << (bb.ymin() + (y + 1.0) * voxel_size) << " " << (bb.zmin() + (z + 1.0)*voxel_size) << "\n"; + stream << (bb.xmin() + (x + 1.0) * voxel_size) << " " << (bb.ymin() + (y + 1.0) * voxel_size) << " " << (bb.zmin() + (z + 1.0)*voxel_size) << "\n"; + } + + // export faces + 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))) continue; + stream << "4 " << (count * 8 + 0) << " " << (count * 8 + 1) << " " << (count * 8 + 3) << " " << (count * 8 + 2) << "\n"; + stream << "4 " << (count * 8 + 4) << " " << (count * 8 + 5) << " " << (count * 8 + 7) << " " << (count * 8 + 6) << "\n"; + stream << "4 " << (count * 8 + 0) << " " << (count * 8 + 1) << " " << (count * 8 + 5) << " " << (count * 8 + 4) << "\n"; + stream << "4 " << (count * 8 + 2) << " " << (count * 8 + 3) << " " << (count * 8 + 7) << " " << (count * 8 + 6) << "\n"; + stream << "4 " << (count * 8 + 0) << " " << (count * 8 + 2) << " " << (count * 8 + 6) << " " << (count * 8 + 4) << "\n"; + stream << "4 " << (count * 8 + 1) << " " << (count * 8 + 3) << " " << (count * 8 + 7) << " " << (count * 8 + 5) << "\n"; + + count++; + } + + 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 = Polygon_mesh_processing::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 = Polygon_mesh_processing::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(Polygon_mesh_processing::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 PkgConvexDecomposition3Ref + * + * \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 512.} + * \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_volumes} + * \cgalParamDescription{maximum number of convex volumes 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 = Polygon_mesh_processing::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 CGAL + +#endif // CGAL_SURFACE_MESH_DECOMPOSITION_APPROXIMATE_CONVEX_DECOMPOSITION_H diff --git a/Surface_mesh_decomposition/package_info/Surface_mesh_decomposition/copyright b/Surface_mesh_decomposition/package_info/Surface_mesh_decomposition/copyright new file mode 100644 index 00000000000..82b520cf64c --- /dev/null +++ b/Surface_mesh_decomposition/package_info/Surface_mesh_decomposition/copyright @@ -0,0 +1 @@ +GeometryFactory \ No newline at end of file diff --git a/Surface_mesh_decomposition/package_info/Surface_mesh_decomposition/dependencies b/Surface_mesh_decomposition/package_info/Surface_mesh_decomposition/dependencies new file mode 100644 index 00000000000..88ebdd57cc6 --- /dev/null +++ b/Surface_mesh_decomposition/package_info/Surface_mesh_decomposition/dependencies @@ -0,0 +1,25 @@ +Algebraic_foundations +Arithmetic_kernel +BGL +CGAL_Core +Cartesian_kernel +Circulator +Convex_decomposition_3 +Convex_hull_3 +Distance_3 +HalfedgeDS +Hash_map +Homogeneous_kernel +Installation +Intersections_3 +Interval_support +Kernel_23 +Modifier +Number_types +Polyhedron +Polygon_mesh_processing +Profiling_tools +Property_map +STL_Extension +Spatial_sorting +Stream_support diff --git a/Surface_mesh_decomposition/package_info/Surface_mesh_decomposition/license.txt b/Surface_mesh_decomposition/package_info/Surface_mesh_decomposition/license.txt new file mode 100644 index 00000000000..ba1f6f2dc96 --- /dev/null +++ b/Surface_mesh_decomposition/package_info/Surface_mesh_decomposition/license.txt @@ -0,0 +1 @@ +GPL (v3 or later) \ No newline at end of file diff --git a/Surface_mesh_decomposition/package_info/Surface_mesh_decomposition/maintainer b/Surface_mesh_decomposition/package_info/Surface_mesh_decomposition/maintainer new file mode 100644 index 00000000000..82b520cf64c --- /dev/null +++ b/Surface_mesh_decomposition/package_info/Surface_mesh_decomposition/maintainer @@ -0,0 +1 @@ +GeometryFactory \ No newline at end of file diff --git a/Surface_mesh_decomposition/test/Surface_mesh_decomposition/CMakeLists.txt b/Surface_mesh_decomposition/test/Surface_mesh_decomposition/CMakeLists.txt new file mode 100644 index 00000000000..5be1aa914d4 --- /dev/null +++ b/Surface_mesh_decomposition/test/Surface_mesh_decomposition/CMakeLists.txt @@ -0,0 +1,10 @@ +# Created by the script cgal_create_CMakeLists +# This is the CMake script for compiling a set of CGAL applications. + +cmake_minimum_required(VERSION 3.12...3.31) +project(Surface_mesh_decomposition_Examples) + +# CGAL and its components +find_package(CGAL REQUIRED) +create_single_source_cgal_program("test_acd.cpp") + diff --git a/Surface_mesh_decomposition/test/Surface_mesh_decomposition/test_acd.cpp b/Surface_mesh_decomposition/test/Surface_mesh_decomposition/test_acd.cpp new file mode 100644 index 00000000000..93100cf5cc5 --- /dev/null +++ b/Surface_mesh_decomposition/test/Surface_mesh_decomposition/test_acd.cpp @@ -0,0 +1,58 @@ +#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; + + std::vector convex_hulls; + + // Try with empty mesh + CGAL::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())); + + assert(convex_hulls.size() == 0); + + if (!PMP::IO::read_polygon_mesh(filename, mesh)) { + std::cerr << "Invalid input." << std::endl; + return 1; + } + + convex_hulls.reserve(9); + + CGAL::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())); + + assert(convex_hulls.size() > 0); + assert(convex_hulls.size() <= 9); + + return 0; +}