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 98f1e2e2d6a..16e859e049d 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 @@ -6,19 +6,31 @@ namespace CGAL { \anchor Chapter_Convex_Decomposition_of_Polyhedra \anchor chapterConvexDecomposition3 \cgalAutoToc -\author Peter Hachenberger +\author Peter Hachenberger, Sven Oesau \section Convex_decomposition_3Introduction Introduction For many applications on non-convex polyhedra, there are efficient solutions that first decompose the polyhedron into convex pieces. As an example, the Minkowski sum of two polyhedra can be computed by -decomposing both polyhedra into convex pieces, compute pair-wise -Minkowski sums of the convex pieces, and unite the pair-wise sums. +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 their number. +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 + +The method 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. @@ -29,16 +41,16 @@ optimal \cgalCite{c-cpplb-84}. Vertical decomposition based on the insertion of vertical facets (viewed from the top). Left: Non-convex polyhedron. Middle: Non-vertical reflex edges have been resolved. Right: Vertical reflex edges have been resolved. The sub-volumes are convex. \cgalFigureEnd -Our decomposition runs in two steps. In the first step, each +The decomposition runs in two steps. In the first step, each non-vertical reflex edge \f$ e\f$ is resolved by insertion of vertical -facets through \f$ e\f$. In the second step, we do the same with the -vertical reflex edges. \cgalFigureRef{figverticalDecomposition} +facets through \f$ e\f$. In the second step, the +vertical reflex edges are handled in the same way. \cgalFigureRef{figverticalDecomposition} illustrates the two steps. -At the moment our implementation is restricted to the decomposition of -bounded polyhedra. An extension to unbounded polyhedra is planned. +At the moment the implementation is restricted to the decomposition of +bounded polyhedra. -\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 @@ -74,11 +86,122 @@ of accessing the convex pieces is to convert them into separate Nef polyhedra, as illustrated by the example code given below. Note that due to the restriction to bounded polyhedra, the use of -extended kernels is unnecessary and expensive. We therefore do not -support the use of extended kernels in the convex decomposition. +extended kernels is unnecessary and expensive. Therefore the use of +extended kernels in the convex decomposition is not supported. + +\subsection Convex_decomposition_3Example Example + +This example shows the usage of `CGAL::convex_decomposition_3(Nef_polyhedron_3&)` to decompose a Nef polyhedron into convex parts. \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: input mesh, 6, 9, and 12 convex volumes +\cgalFigureCaptionEnd + +The H-VACD method \cgalCite{cgal:mamou2016volumetric}, "Hierarchical volumetrix approximate convex decomposition", 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 it. +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 may slightly overlap and their union 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 voxel cells overlapping with triangles 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-axes 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. 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. + +To optionally improve the fit of the convex volumes, they can be refitted to the mesh before starting the second phase. 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_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_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. +- `refitting`: The convex hulls can be refitted after the splitting phase to more tightly enclose the input mesh. It increases the running time, but significantly reduces the overhead volume included by the computed convex volulmes. It is enabled by default. +- `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 and refitting to the convex hull of the input mesh on volume overhead: + + +| Data set | Split location | Refitting | 4 volumes | 6 volumes | 8 volumes | 10 volumes | 12 volumes | +| :---- | ----: | ----: | -: | -: | -: | -: | -: | +| Camel | Concavity | + | 0.6951 | 0.4316 | 0.3016 | 0.2173 | 0.1939 | +| Camel | Concavity | - | 0.9932 | 0.7482 | 0.6174 | 0.5507 | 0.5261 | +| Elephant | Concavity | + | 0.7897 | 0.6505 | 0.4973 | 0.3986 | 0.3299 | +| Elephant | Concavity | - | 1.2140 | 1.0028 | 0.8071 | 0.7290 | 0.6870 | +| Triceratops | Concavity | + | 0.5952 | 0.3978 | 0.3548 | 0.2385 | 0.2057 | +| Triceratops | Concavity | - | 1.0073 | 0.7490 | 0.7035 | 0.5966 | 0.5429 | + + +And by using the mid split: + + +| Data set | Split location | Refitting | 4 volumes | 6 volumes | 8 volumes | 10 volumes | 12 volumes | +| :---- | ----: | ----: | -: | -: | -: | -: | -: | +| Camel | Mid | + | 0.9970 | 0.5762 | 0.5332 | 0.4040 | 0.2390 | +| Camel | Mid | - | 1.0875 | 0.9280 | 0.8394 | 0.6801 | 0.5529 | +| Elephant | Mid | + | 0.7424 | 0.5672 | 0.4207 | 0.3619 | 0.3043 | +| Elephant | Mid | - | 1.2075 | 1.0410 | 0.8247 | 0.6950 | 0.6434 | +| Triceratops | Mid | + | 0.7671 | 0.5096 | 0.3965 | 0.3309 | 0.2494 | +| Triceratops | Mid | - | 1.0470 | 0.8990 | 0.6749 | 0.5803 | 0.5230 | + +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 + +The example shows the approximate convex decomposition of the knot2.off mesh into 9 convex volumes that are saved as off files. + +\cgalExample{Convex_decomposition_3/approximate_convex_decomposition.cpp } + + +\section Convex_decomposition_3_history Design and Implementation History + +This package was created by Peter Hachenberger with the `CGAL::convex_decomposition_3()` method. In 2025, it has been extended by Sven Oesau with the `CGAL::approximate_convex_decomposition()` method. + */ } /* 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..143c35c2ff6 100644 --- a/Convex_decomposition_3/doc/Convex_decomposition_3/Doxyfile.in +++ b/Convex_decomposition_3/doc/Convex_decomposition_3/Doxyfile.in @@ -3,3 +3,5 @@ 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 diff --git a/Convex_decomposition_3/doc/Convex_decomposition_3/PackageDescription.txt b/Convex_decomposition_3/doc/Convex_decomposition_3/PackageDescription.txt index fac4f8479c7..b3cd73b80fa 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 @@ -20,7 +21,8 @@ \cgalClassifedRefPages \cgalCRPSection{Functions} -- `CGAL::convex_decomposition_3(Nef_polyhedron_3& N)` +- `CGAL::approximate_convex_decomposition(FaceGraph)` +- `CGAL::convex_decomposition_3(NefPolyhedron_3)` */ diff --git a/Convex_decomposition_3/doc/Convex_decomposition_3/dependencies b/Convex_decomposition_3/doc/Convex_decomposition_3/dependencies index 94e7ccf5922..9f1e130e7fe 100644 --- a/Convex_decomposition_3/doc/Convex_decomposition_3/dependencies +++ b/Convex_decomposition_3/doc/Convex_decomposition_3/dependencies @@ -5,3 +5,5 @@ Algebraic_foundations Circulator Stream_support Nef_3 +BGL +Polygon_mesh_processing diff --git a/Convex_decomposition_3/doc/Convex_decomposition_3/examples.txt b/Convex_decomposition_3/doc/Convex_decomposition_3/examples.txt index b40169d2c81..3c4f462dfed 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/approximate_convex_decomposition.cpp \example Convex_decomposition_3/list_of_convex_parts.cpp */ diff --git a/Convex_decomposition_3/examples/Convex_decomposition_3/CMakeLists.txt b/Convex_decomposition_3/examples/Convex_decomposition_3/CMakeLists.txt index 3e4532f6ffc..04bc377f7cb 100644 --- a/Convex_decomposition_3/examples/Convex_decomposition_3/CMakeLists.txt +++ b/Convex_decomposition_3/examples/Convex_decomposition_3/CMakeLists.txt @@ -14,3 +14,12 @@ file( foreach(cppfile ${cppfiles}) create_single_source_cgal_program("${cppfile}") endforeach() + +find_package(TBB QUIET) +include(CGAL_TBB_support) + +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/Convex_decomposition_3/examples/Convex_decomposition_3/approximate_convex_decomposition.cpp b/Convex_decomposition_3/examples/Convex_decomposition_3/approximate_convex_decomposition.cpp new file mode 100644 index 00000000000..9cb86c54e3b --- /dev/null +++ b/Convex_decomposition_3/examples/Convex_decomposition_3/approximate_convex_decomposition.cpp @@ -0,0 +1,48 @@ +#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_volumes; + convex_volumes.reserve(9); + + CGAL::approximate_convex_decomposition(mesh, std::back_inserter(convex_volumes), + CGAL::parameters::maximum_depth(10) + .volume_error(0.1) + .maximum_number_of_convex_volumes(9) + .split_at_concavity(true) + .refitting(true) + .maximum_number_of_voxels(1000000) + .concurrency_tag(CGAL::Parallel_if_available_tag())); + + for (std::size_t i = 0;i` */ -template -void convex_decomposition_3(Nef_polyhedron& N) +template +void convex_decomposition_3(NefPolyhedron_3& N) { - typedef typename Nef_polyhedron::SNC_structure SNC_structure; + typedef typename NefPolyhedron_3::SNC_structure SNC_structure; typedef typename SNC_structure::Halfedge_handle Halfedge_handle; - typedef typename Nef_polyhedron::Point_3 Point_3; - typedef typename Nef_polyhedron::Vector_3 Vector_3; - typedef typename Nef_polyhedron::Sphere_point Sphere_point; - typedef typename Nef_polyhedron::FT FT; + typedef typename NefPolyhedron_3::Point_3 Point_3; + typedef typename NefPolyhedron_3::Vector_3 Vector_3; + typedef typename NefPolyhedron_3::Sphere_point Sphere_point; + typedef typename NefPolyhedron_3::FT FT; - typedef typename CGAL::Single_wall_creator Single_wall; - typedef typename CGAL::YVertical_wall_builder YVertical_wall_builder; - typedef typename CGAL::Reflex_vertex_searcher Reflex_vertex_searcher; - typedef typename CGAL::Ray_hit_generator2 Ray_hit2; - typedef typename CGAL::External_structure_builder External_structure_builder; - typedef typename CGAL::SFace_separator SFace_separator; + typedef typename CGAL::Single_wall_creator Single_wall; + typedef typename CGAL::YVertical_wall_builder YVertical_wall_builder; + typedef typename CGAL::Reflex_vertex_searcher Reflex_vertex_searcher; + typedef typename CGAL::Ray_hit_generator2 Ray_hit2; + typedef typename CGAL::External_structure_builder External_structure_builder; + typedef typename CGAL::SFace_separator SFace_separator; typedef Compare_halfedges_in_reflex_edge_sorter > Less_edges; @@ -86,11 +88,11 @@ void convex_decomposition_3(Nef_polyhedron& N) typedef typename Positively_sorted_set::const_iterator Positive_reflex_edge_iterator; typedef typename Negatively_sorted_set::const_iterator Negative_reflex_edge_iterator; - typedef typename CGAL::Reflex_edge_searcher + typedef typename CGAL::Reflex_edge_searcher Reflex_edge_searcher; - typedef typename CGAL::Edge_sorter, Negatively_sorted_set> Edge_sorter; - typedef typename CGAL::Edge_sorter, Positively_sorted_set> Edge_sorter2; + typedef typename CGAL::Edge_sorter, Negatively_sorted_set> Edge_sorter; + typedef typename CGAL::Edge_sorter, Positively_sorted_set> Edge_sorter2; External_structure_builder esb; SFace_separator sf_sep; @@ -176,10 +178,10 @@ void convex_decomposition_3(Nef_polyhedron& N) N.delegate(esb); - CGAL_assertion_code(typename Nef_polyhedron::SHalfedge_const_iterator cse); + CGAL_assertion_code(typename NefPolyhedron_3::SHalfedge_const_iterator cse); CGAL_assertion_code(CGAL_forall_shalfedges(cse, N) if(cse->incident_sface()->mark()) - CGAL_assertion(!CGAL::is_reflex_sedge_in_any_direction(cse))); + CGAL_assertion(!CGAL::is_reflex_sedge_in_any_direction(cse))); } } //namespace CGAL diff --git a/Documentation/doc/biblio/cgal_manual.bib b/Documentation/doc/biblio/cgal_manual.bib index a16cfe0015f..4813821987b 100644 --- a/Documentation/doc/biblio/cgal_manual.bib +++ b/Documentation/doc/biblio/cgal_manual.bib @@ -3516,6 +3516,16 @@ pages = "207--221" year = {2015} } +@article{cgal:mamou2016volumetric, + title={Volumetric hierarchical approximate convex decomposition}, + author={Mamou, Khaled and Lengyel, E and Peters, A}, + journal={Game engine gems}, + volume={3}, + pages={141--158}, + year={2016}, + publisher={AK Peters} +} + % ---------------------------------------------------------------------------- % END OF BIBFILE % ---------------------------------------------------------------------------- diff --git a/Installation/CHANGES.md b/Installation/CHANGES.md index 5d5986ab124..8cc9e7d8b5d 100644 --- a/Installation/CHANGES.md +++ b/Installation/CHANGES.md @@ -10,6 +10,10 @@ Release date: July 2026 - Made the `Arr_linear_traits_2` a model of the new concept. - Added overloads of `draw(Arrangement_on_surface_2& arr, Bbox& bbox, ...)` that enable the drawing of arrangements induced by unbounded curves. +### [Convex Decomposition of Polyhedra](https://doc.cgal.org/6.2/Manual/packages.html#PkgConvexDecomposition3) +- Added the function `CGAL::approximate_convex_decomposition()` that computes a set of convex volumes that cover an input mesh. + It performs a hierarchical splitting of the convex hull guided by an volumetric error and performs a subsequent merging of convex hulls to obtain the targeted count. + ### [Linear Cell Complex](https://doc.cgal.org/6.2/Manual/packages.html#PkgLinearCellComplex) - **API Changes**: The following import functions have been deprecated and renamed for better naming clarity and consistency: 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/Lab/demo/Lab/Plugins/Convex_decomposition/Approximate_convex_decomposition_dialog.ui b/Lab/demo/Lab/Plugins/Convex_decomposition/Approximate_convex_decomposition_dialog.ui new file mode 100644 index 00000000000..bb7668ebf5e --- /dev/null +++ b/Lab/demo/Lab/Plugins/Convex_decomposition/Approximate_convex_decomposition_dialog.ui @@ -0,0 +1,211 @@ + + + Approximate_convex_decomposition_dialog + + + + 0 + 0 + 448 + 225 + + + + Approximate Convex Decomposition + + + + + + + 15 + + + + Approximate Convex Decomposition + + + + + + + Maximum number of components + + + Qt::AlignmentFlag::AlignRight|Qt::AlignmentFlag::AlignTrailing|Qt::AlignmentFlag::AlignVCenter + + + + + + + Number of voxels + + + Qt::AlignmentFlag::AlignRight|Qt::AlignmentFlag::AlignTrailing|Qt::AlignmentFlag::AlignVCenter + + + + + + + Maximum decomposition depth + + + Qt::AlignmentFlag::AlignRight|Qt::AlignmentFlag::AlignTrailing|Qt::AlignmentFlag::AlignVCenter + + + + + + + 1 + + + 10000 + + + 16 + + + + + + + Volume error + + + Qt::AlignmentFlag::AlignRight|Qt::AlignmentFlag::AlignTrailing|Qt::AlignmentFlag::AlignVCenter + + + + + + + Qt::Orientation::Horizontal + + + QDialogButtonBox::StandardButton::Cancel|QDialogButtonBox::StandardButton::Ok + + + + + + + 2 + + + 50 + + + 10 + + + + + + + 0.01 + + + + + + + Qt::LayoutDirection::RightToLeft + + + Split at concavity + + + true + + + + + + + 512 + + + 1000000000 + + + 100000 + + + 1000000 + + + + + + + Qt::Orientation::Vertical + + + + 20 + 40 + + + + + + + + Qt::LayoutDirection::RightToLeft + + + Refit to mesh + + + true + + + + + + + + DoubleEdit + QLineEdit +
CGAL_double_edit.h
+ + + + + + buttonBox + accepted() + Approximate_convex_decomposition_dialog + accept() + + + 384 + 191 + + + 157 + 195 + + + + + buttonBox + rejected() + Approximate_convex_decomposition_dialog + reject() + + + 384 + 191 + + + 286 + 195 + + + + + 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 new file mode 100644 index 00000000000..8d31ef536dc --- /dev/null +++ b/Lab/demo/Lab/Plugins/Convex_decomposition/Approximate_convex_decomposition_plugin.cpp @@ -0,0 +1,157 @@ +#include "config.h" + +#include +#include +#include +#include "ui_Approximate_convex_decomposition_dialog.h" + +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include "Scene_surface_mesh_item.h" +#include "Color_map.h" +#include +#include +#include +#include + +#include +#include +#include +#include + +using namespace CGAL::Three; + +class CGAL_Lab_approximate_convex_decomposition_plugin + : public QObject, + protected CGAL_Lab_plugin_interface +{ + Q_OBJECT + Q_INTERFACES(CGAL::Three::CGAL_Lab_plugin_interface) + Q_PLUGIN_METADATA(IID "com.geometryfactory.CGALLab.PluginInterface/1.0") + + using Convex_hull = std::pair, std::vector > >; + +private: + QAction* actionApproximateConvexDecomposition; + + Scene_interface *scene; + QMainWindow *mw; + +public: + void init(QMainWindow* mainWindow, + Scene_interface* scene_interface, + Messages_interface*) + { + this->scene = scene_interface; + this->mw = mainWindow; + + actionApproximateConvexDecomposition = new QAction(tr("Approximate Convex Decomposition"), mw); + actionApproximateConvexDecomposition->setProperty("subMenuName", "Polygon Mesh Processing"); + connect(actionApproximateConvexDecomposition, SIGNAL(triggered()), + this, SLOT(approximate_convex_decomposition())); + } + + bool applicable(QAction* action) const + { + if(action == actionApproximateConvexDecomposition) + { + if(scene->selectionIndices().size() == 1) + { + const int index = scene->mainSelectionIndex(); + return (qobject_cast(scene->item(index))); + } + } + + return false; + } + + QList actions() const + { + return QList() << actionApproximateConvexDecomposition; + } + +public Q_SLOTS: + void approximate_convex_decomposition(); +}; // class CGAL_Lab_approximate_convex_decomposition_plugin + +void +CGAL_Lab_approximate_convex_decomposition_plugin:: +approximate_convex_decomposition() +{ + Scene_surface_mesh_item* sm_item = nullptr; + + if (scene->selectionIndices().size() != 1) + return; + + sm_item = qobject_cast(scene->item(scene->selectionIndices().front())); + if (sm_item == nullptr) + return; + + QDialog dialog(mw); + Ui::Approximate_convex_decomposition_dialog ui; + ui.setupUi(&dialog); + + int i = dialog.exec(); + if (i == QDialog::Rejected) + return; + + const unsigned int maximum_depth = static_cast(ui.maximumDepth->value()); + const unsigned int maximum_convex_hulls = static_cast(ui.maximumConvexHulls->value()); + const unsigned int num_voxels = static_cast(ui.numVoxels->value()); + const double volume_error = ui.volumeError->value(); + const bool split_concavity = ui.splitConcavity->isChecked(); + const bool refitting = ui.refitting->isChecked(); + + QApplication::setOverrideCursor(Qt::WaitCursor); + + std::vector convex_volumes; + convex_volumes.reserve(maximum_convex_hulls); + + CGAL::approximate_convex_decomposition(*(sm_item->face_graph()), std::back_inserter(convex_volumes), + CGAL::parameters::maximum_depth(maximum_depth) + .volume_error(volume_error) + .maximum_number_of_convex_volumes(maximum_convex_hulls) + .split_at_concavity(split_concavity) + .refitting(refitting) + .maximum_number_of_voxels(num_voxels) + .concurrency_tag(CGAL::Parallel_if_available_tag())); + + std::vector distinct_colors; + // the first color is either the background or the unique domain + + compute_deterministic_color_map(QColor(80, 250, 80), convex_volumes.size(), std::back_inserter(distinct_colors)); + + Scene_group_item* group = new Scene_group_item(tr("%1 %2 decomposition").arg(sm_item->name()).arg(maximum_convex_hulls)); + scene->addItem(group); + + for (std::size_t i = 0; i < convex_volumes.size(); i++) { + Convex_hull& ch = convex_volumes[i]; + SMesh sm; + CGAL::Polygon_mesh_processing::polygon_soup_to_polygon_mesh(ch.first, ch.second, sm); + + Scene_surface_mesh_item* component_item = new Scene_surface_mesh_item(std::move(sm)); + component_item->setName(tr("%1 %2").arg(sm_item->name()).arg(i)); + component_item->setColor(distinct_colors[i]); + Three::scene()->addItem(component_item); + Three::scene()->changeGroup(component_item, group); + } + + QApplication::restoreOverrideCursor(); + + QApplication::restoreOverrideCursor(); + + QApplication::setOverrideCursor(Qt::BusyCursor); + QApplication::restoreOverrideCursor(); +} + + +#include "Approximate_convex_decomposition_plugin.moc" diff --git a/Lab/demo/Lab/Plugins/Convex_decomposition/CMakeLists.txt b/Lab/demo/Lab/Plugins/Convex_decomposition/CMakeLists.txt index b17a2a80c14..d5e2be10a60 100644 --- a/Lab/demo/Lab/Plugins/Convex_decomposition/CMakeLists.txt +++ b/Lab/demo/Lab/Plugins/Convex_decomposition/CMakeLists.txt @@ -1,5 +1,17 @@ include(CGALlab_macros) -cgal_lab_plugin(nef_plugin Nef_plugin) -target_link_libraries(nef_plugin PRIVATE scene_nef_polyhedron_item scene_surface_mesh_item) +set(CMAKE_AUTOMOC ON) +cgal_lab_plugin(nef_plugin Nef_plugin) + +qt6_wrap_ui(acdUI_FILES Approximate_convex_decomposition_dialog.ui) +cgal_lab_plugin(acd_plugin Approximate_convex_decomposition_plugin ${acdUI_FILES}) + +target_link_libraries(nef_plugin PRIVATE scene_nef_polyhedron_item scene_surface_mesh_item) +target_link_libraries(acd_plugin PRIVATE scene_surface_mesh_item) + +find_package(TBB QUIET) +include(CGAL_TBB_support) +if(TARGET CGAL::TBB_support) + target_link_libraries(acd_plugin PRIVATE CGAL::TBB_support) +endif() diff --git a/Nef_3/doc/Nef_3/CGAL/Nef_nary_union_3.h b/Nef_3/doc/Nef_3/CGAL/Nef_nary_union_3.h index 3919d156f76..013846df886 100644 --- a/Nef_3/doc/Nef_3/CGAL/Nef_nary_union_3.h +++ b/Nef_3/doc/Nef_3/CGAL/Nef_nary_union_3.h @@ -20,7 +20,7 @@ polyhedra in a sorted way and add them one by one to `Nef_nary_union_3`. */ -template< typename Nef_polyhedron_3 > +template< typename NefPolyhedron_3 > class Nef_nary_union_3 { public: diff --git a/Nef_3/doc/Nef_3/CGAL/Nef_polyhedron_3.h b/Nef_3/doc/Nef_3/CGAL/Nef_polyhedron_3.h index 06f67260279..550f2c9c719 100644 --- a/Nef_3/doc/Nef_3/CGAL/Nef_polyhedron_3.h +++ b/Nef_3/doc/Nef_3/CGAL/Nef_polyhedron_3.h @@ -34,7 +34,7 @@ namespace CGAL { `Quotient` or any other number type modeling \f$\mathbb{Q}\f$. The second parameter and the third parameter are for future considerations. - Neither `Nef_polyhedronItems_3` nor `Nef_polyhedronMarks` is + Neither `NefPolyhedronItems_3` nor `Nef_polyhedronMarks` is specified, yet. Do not use any other than the default types for these two template parameters. @@ -49,8 +49,8 @@ namespace CGAL { \sa `CGAL::Polyhedron_3` */ -template< class Nef_polyhedronTraits_3, - class Nef_polyhedronItems_3 = CGAL::Default_items +template< class NefPolyhedronTraits_3, + class NefPolyhedronItems_3 = CGAL::Default_items class Nef_polyhedronMarks = bool > class Nef_polyhedron_3 { public: @@ -79,7 +79,7 @@ public: function `twin()` returns the opposite halfedge. Looking at the incidence structure on a sphere map, the member function - `out_sedge` returns the first outgoing shalfedge, and `incident_sface` + `out_sedge()` returns the first outgoing shalfedge, and `incident_sface()` returns the incident sface. \cgalHeading{Creation} @@ -946,7 +946,7 @@ public: }; /* end Volume */ /*! - traits class selected for `Nef_polyhedronTraits_3`. + traits class selected for `NefPolyhedronTraits_3`. */ typedef unspecified_type Traits; diff --git a/Nef_3/doc/Nef_3/CGAL/OFF_to_nef_3.h b/Nef_3/doc/Nef_3/CGAL/OFF_to_nef_3.h index ee4d429515b..5f1f5c2bf21 100644 --- a/Nef_3/doc/Nef_3/CGAL/OFF_to_nef_3.h +++ b/Nef_3/doc/Nef_3/CGAL/OFF_to_nef_3.h @@ -4,7 +4,7 @@ namespace CGAL { \ingroup PkgNef3IOFunctions This function creates a 3D Nef polyhedron from an OFF file which -is read from input stream `in`. The purpose of `OFF_to_nef_3` +is read from input stream `in`. The purpose of `OFF_to_nef_3()` is to create a Nef polyhedron from an OFF file that cannot be handled by the `Nef_polyhedron_3` constructors. It handles double coordinates while using a homogeneous kernel, non-coplanar facets, @@ -12,10 +12,12 @@ surfaces with boundaries, self-intersecting surfaces, and single facets. Every closed volume gets marked. The function returns the number of facets it could not handle. +@tparam NefPolyhedron_3 an object of type `Nef_polyhedron_3`. + \sa `CGAL::Nef_polyhedron_3` */ -template -std::size_t OFF_to_nef_3(std::istream& in, Nef_polyhedron_3& N); +template +std::size_t OFF_to_nef_3(std::istream& in, NefPolyhedron_3& N); } /* namespace CGAL */ diff --git a/Nef_3/doc/Nef_3/CGAL/boost/graph/convert_nef_polyhedron_to_polygon_mesh.h b/Nef_3/doc/Nef_3/CGAL/boost/graph/convert_nef_polyhedron_to_polygon_mesh.h index 558c6321280..712c287c0a4 100644 --- a/Nef_3/doc/Nef_3/CGAL/boost/graph/convert_nef_polyhedron_to_polygon_mesh.h +++ b/Nef_3/doc/Nef_3/CGAL/boost/graph/convert_nef_polyhedron_to_polygon_mesh.h @@ -1,12 +1,12 @@ namespace CGAL { /// \ingroup PkgNef3IOFunctions -/// Converts an object of type `Nef_polyhedron_3` into a polygon mesh model of `MutableFaceGraph`. +/// Converts an object of type `NefPolyhedron_3` into a polygon mesh model of `MutableFaceGraph`. /// Note that contrary to `Nef_polyhedron_3::convert_to_polyhedron()`, the output is not triangulated /// (but faces with more than one connected component of the boundary). -/// The polygon mesh can be triangulated by setting `triangulate_all_faces` to `true` or by calling the function `triangulate_faces()`. -/// @tparam Nef_polyhedron an object of type `Nef_polyhedron_3`. -/// @tparam Polygon_mesh a model of `MutableFaceGraph` with an internal property map for `CGAL::vertex_point_t`. +/// The polygon mesh can be triangulated by setting `triangulate_all_faces` to `true` or by calling the function `Polygon_mesh_processing::triangulate_faces()`. +/// @tparam NefPolyhedron_3 an object of type `Nef_polyhedron_3`. +/// @tparam PolygonMesh a model of `MutableFaceGraph` with an internal property map for `CGAL::vertex_point_t`. /// /// The points from `nef` to `pm` are converted using /// `CGAL::Cartesian_converter`. @@ -17,16 +17,16 @@ namespace CGAL { /// @param pm the output. /// @param triangulate_all_faces indicates whether all the faces must be triangulated. /// -/// \pre `Polygon_mesh` must have an internal point property map with value type being `Nef_polyhedron_3::Point_3`. +/// \pre `PolygonMesh` must have an internal point property map with value type being `NefPolyhedron_3::Point_3`. /// \pre `nef.simple()` - template - void convert_nef_polyhedron_to_polygon_mesh(const Nef_polyhedron& nef, Polygon_mesh& pm, bool triangulate_all_faces = false); + template + void convert_nef_polyhedron_to_polygon_mesh(const NefPolyhedron_3& nef, PolygonMesh& pm, bool triangulate_all_faces = false); /// \ingroup PkgNef3IOFunctions - /// Converts an object of type `Nef_polyhedron_3` into a polygon soup. + /// Converts an object of type `NefPolyhedron_3` into a polygon soup. /// The polygons can be triangulated by setting `triangulate_all_faces` to `true`. - /// @tparam Nef_polyhedron an object of type `Nef_polyhedron_3`. + /// @tparam NefPolyhedron_3 an object of type `Nef_polyhedron_3`. /// @tparam PointRange a model of the concept `BackInsertionSequence` /// whose `value_type` is the point type /// @tparam PolygonRange a model of the concept @@ -43,8 +43,8 @@ namespace CGAL { /// @param points the output points of the soup /// @param polygons the output polygons of the soup. /// @param triangulate_all_faces indicates whether all polygons must be triangulated. - template - void convert_nef_polyhedron_to_polygon_soup(const Nef_polyhedron& nef, + template + void convert_nef_polyhedron_to_polygon_soup(const NefPolyhedron_3& nef, PointRange& points, PolygonRange& polygons, bool triangulate_all_faces = false); diff --git a/Nef_3/doc/Nef_3/dependencies b/Nef_3/doc/Nef_3/dependencies index 467cacb769e..0ec6513d194 100644 --- a/Nef_3/doc/Nef_3/dependencies +++ b/Nef_3/doc/Nef_3/dependencies @@ -12,3 +12,4 @@ BGL Surface_mesh Polygon_mesh_processing GraphicsView +Basic_viewer diff --git a/Nef_S2/doc/Nef_S2/CGAL/Nef_polyhedron_S2.h b/Nef_S2/doc/Nef_S2/CGAL/Nef_polyhedron_S2.h index a8454f79d2e..0adcf175e01 100644 --- a/Nef_S2/doc/Nef_S2/CGAL/Nef_polyhedron_S2.h +++ b/Nef_S2/doc/Nef_S2/CGAL/Nef_polyhedron_S2.h @@ -16,8 +16,8 @@ under the topological operations `boundary`, `closure`, and \cgalHeading{Parameters} \code -template< class Nef_polyhedronTraits_S2, - class Nef_polyhedronItems_S2 = CGAL::SM_items, +template< class NefPolyhedronTraits_S2, + class NefPolyhedronItems_S2 = CGAL::SM_items, class Nef_polyhedronMarks = bool > class Nef_polyhedron_S2; \endcode @@ -30,7 +30,7 @@ modeling \f$\mathbb{Z}\f$, or `Cartesian`, `Simple_cartesian` parametrized with type modeling \f$\mathbb{Q}\f$. The second parameter and the third parameter are for future considerations. -Neither `Nef_polyhedronItems_S2` nor `Nef_polyhedronMarks` is +Neither `NefPolyhedronItems_S2` nor `Nef_polyhedronMarks` is specified, yet. Do not use other than the default types for these two template parameters. @@ -55,8 +55,8 @@ output operator is defined in the file Nef polyhedra are implemented on top of a halfedge data structure and use linear space in the number of vertices, edges and facets. Operations like `empty` take constant time. The operations -`clear`, `complement`, `interior`, `closure`, -`boundary`, `regularization`, input and output take linear +`clear()`, `complement()`, `interior()`, `closure()`, +`boundary()`, `regularization()`, input and output take linear time. All binary set operations and comparison operations take time \cgalBigO{n \log n} where \f$ n\f$ is the size of the output plus the size of the input. diff --git a/Polygon_mesh_processing/test/Polygon_mesh_processing/CMakeLists.txt b/Polygon_mesh_processing/test/Polygon_mesh_processing/CMakeLists.txt index 6fcadef5628..965500699de 100644 --- a/Polygon_mesh_processing/test/Polygon_mesh_processing/CMakeLists.txt +++ b/Polygon_mesh_processing/test/Polygon_mesh_processing/CMakeLists.txt @@ -73,6 +73,7 @@ create_single_source_cgal_program("test_isolevel_refinement.cpp") create_single_source_cgal_program("test_corefinement_nm_bo.cpp") create_single_source_cgal_program("test_corefinement_cavities.cpp") create_single_source_cgal_program("issue_8730.cpp") +create_single_source_cgal_program("test_approximate_convex_decomposition.cpp") create_single_source_cgal_program("issue_7164.cpp") # create_single_source_cgal_program("test_pmp_repair_self_intersections.cpp") @@ -120,6 +121,7 @@ if(TARGET CGAL::TBB_support) target_link_libraries(self_intersection_surface_mesh_test PRIVATE CGAL::TBB_support) target_link_libraries(test_autorefinement PRIVATE CGAL::TBB_support) target_link_libraries(test_snap_rounding PRIVATE CGAL::TBB_support) + target_link_libraries(test_approximate_convex_decomposition PRIVATE CGAL::TBB_support) else() message(STATUS "NOTICE: Intel TBB was not found. Tests will use sequential code.") endif() diff --git a/Polygon_mesh_processing/test/Polygon_mesh_processing/test_approximate_convex_decomposition.cpp b/Polygon_mesh_processing/test/Polygon_mesh_processing/test_approximate_convex_decomposition.cpp new file mode 100644 index 00000000000..71f5cec1fab --- /dev/null +++ b/Polygon_mesh_processing/test/Polygon_mesh_processing/test_approximate_convex_decomposition.cpp @@ -0,0 +1,66 @@ +#include +#include +#include +#include + +#include +#include +#include + +#include + +#include +#include +#include +#include + +using EPICK = CGAL::Exact_predicates_inexact_constructions_kernel; +using EPECK = CGAL::Exact_predicates_exact_constructions_kernel; + +namespace PMP = CGAL::Polygon_mesh_processing; + +template +void test(const std::string &filename, const std::vector n) { + using Point = typename K::Point_3; + + using Convex_hull = std::pair, std::vector > >; + using Mesh = CGAL::Surface_mesh; + + Mesh mesh; + if (!PMP::IO::read_polygon_mesh(filename, mesh)) { + std::cerr << "Invalid input." << std::endl; + assert(false); + } + + for (std::size_t i : n) { + std::vector convex_hulls; + PMP::approximate_convex_decomposition(mesh, std::back_inserter(convex_hulls), + CGAL::parameters::maximum_depth(10) + .volume_error(0.1) + .maximum_number_of_convex_hulls(i) + .split_at_concavity(true) + .maximum_number_of_voxels(1000000) + .concurrency_tag(CGAL::Parallel_if_available_tag())); + std::cout << convex_hulls.size() << std::endl; + assert(convex_hulls.size() == i); + for (std::size_t j = 0; j < convex_hulls.size(); j++) { + Mesh m; + PMP::polygon_soup_to_polygon_mesh(convex_hulls[j].first, convex_hulls[j].second, m); + + assert(CGAL::is_strongly_convex_3(m, typename CGAL::Convex_hull_3::internal::Default_traits_for_Chull_3::type())); + } + } +} + +int main(int argc, char* argv[]) +{ + const std::string filename = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/knot2.off"); + std::cout << "Testing approximate convex decomposition with EPICK Kernel of " << filename << std::endl; + test(filename, {3, 6, 9}); + std::cout << "Testing approximate convex decomposition with EPECK Kernel of " << filename << std::endl; + test(filename, {5}); + + std::cout << "done" << std::endl; + + return 0; +} 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 73c54858bf2..4695d38325e 100644 --- a/STL_Extension/include/CGAL/STL_Extension/internal/parameters_interface.h +++ b/STL_Extension/include/CGAL/STL_Extension/internal/parameters_interface.h @@ -179,6 +179,13 @@ CGAL_add_named_parameter(tiling_t, tiling, tiling) CGAL_add_named_parameter(dimension_t, dimension, dimension) CGAL_add_named_parameter(apply_iterative_snap_rounding_t, apply_iterative_snap_rounding, apply_iterative_snap_rounding) 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(refitting_t, refitting, refitting) +CGAL_add_named_parameter(volume_error_t, volume_error, volume_error) +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) + // List of named parameters that we use in the package 'Surface Mesh Simplification' CGAL_add_named_parameter(get_cost_policy_t, get_cost_policy, get_cost) diff --git a/Surface_mesh_decomposition/include/CGAL/approximate_convex_decomposition.h b/Surface_mesh_decomposition/include/CGAL/approximate_convex_decomposition.h new file mode 100644 index 00000000000..c63ca2f6acb --- /dev/null +++ b/Surface_mesh_decomposition/include/CGAL/approximate_convex_decomposition.h @@ -0,0 +1,1956 @@ +// Copyright (c) 2025 GeometryFactory (France). +// All rights reserved. +// +// This file is part of CGAL (www.cgal.org). +// +// $URL$ +// $Id$ +// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial +// +// Author(s) : Sven Oesau +// + +#ifndef CGAL_SURFACE_MESH_DECOMPOSITION_APPROXIMATE_CONVEX_DECOMPOSITION_H +#define CGAL_SURFACE_MESH_DECOMPOSITION_APPROXIMATE_CONVEX_DECOMPOSITION_H + +#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 + +#include + +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 +typename GeomTraits::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 + ); +} + + +template +typename GeomTraits::Iso_cuboid_3 bbox_voxel_bbox(const Bbox_uint& voxelbb, const Bbox_3& bb, const typename GeomTraits::FT& voxel_size) { + double vs = to_double(voxel_size); + return Bbox_3( + bb.xmin() + voxelbb.lower[0] * vs, + bb.ymin() + voxelbb.lower[1] * vs, + bb.zmin() + voxelbb.lower[2] * vs, + bb.xmin() + (voxelbb.upper[0] + 1) * vs, + bb.ymin() + (voxelbb.upper[1] + 1) * vs, + bb.zmin() + (voxelbb.upper[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 = 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; + } + } +} + +template +struct Convex_hull_candidate { + using FT = typename GeomTraits::FT; + using Point_3 = typename GeomTraits::Point_3; + typename GeomTraits::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 }) {} + Candidate(std::size_t depth, Bbox_uint &bbox) : depth(depth), bbox(bbox) {} +}; + +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++) { + typename GeomTraits::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 shrink_candidates(const FaceGraph& tmesh, std::vector>& candidates, const Bbox_3& bbox, const typename GeomTraits::FT& voxel_size, CGAL::Sequential_tag) { + using face_descriptor = typename boost::graph_traits::face_descriptor; + + using Point_3 = typename GeomTraits::Point_3; + using Segment_3 = typename GeomTraits::Segment_3; + using Triangle_3 = typename GeomTraits::Triangle_3; + using FT = typename GeomTraits::FT; + + using Primitive = CGAL::AABB_face_graph_triangle_primitive; + using Traits = CGAL::AABB_traits_3; + using Tree = CGAL::AABB_tree; + + Tree tree(faces(tmesh).first, faces(tmesh).second, tmesh); + + for (Candidate& c : candidates) { + std::vector pts; + pts.reserve(c.new_surface.size() * 8); + using Box = typename GeomTraits::Iso_cuboid_3; + using IP_id = typename Tree::template Intersection_and_primitive_id::Type; + std::vector intersections; + tree.all_intersections(c.ch.bbox, std::back_inserter(intersections)); + + std::vector corners(8); + std::vector taken(8, false); + + corners[0] = Point_3(c.ch.bbox.xmin(), c.ch.bbox.ymin(), c.ch.bbox.zmin()); + corners[1] = Point_3(c.ch.bbox.xmin(), c.ch.bbox.ymax(), c.ch.bbox.zmin()); + corners[2] = Point_3(c.ch.bbox.xmin(), c.ch.bbox.ymin(), c.ch.bbox.zmax()); + corners[3] = Point_3(c.ch.bbox.xmin(), c.ch.bbox.ymax(), c.ch.bbox.zmax()); + corners[4] = Point_3(c.ch.bbox.xmax(), c.ch.bbox.ymin(), c.ch.bbox.zmin()); + corners[5] = Point_3(c.ch.bbox.xmax(), c.ch.bbox.ymax(), c.ch.bbox.zmin()); + corners[6] = Point_3(c.ch.bbox.xmax(), c.ch.bbox.ymin(), c.ch.bbox.zmax()); + corners[7] = Point_3(c.ch.bbox.xmax(), c.ch.bbox.ymax(), c.ch.bbox.zmax()); + + for (auto& i : intersections) { + const Point_3* p; + const Segment_3* s; + const Triangle_3* t; + const std::vector* v; + if (p = std::get_if(&(i.first))) { + pts.push_back(*p); + } + else if (s = std::get_if(&(i.first))) { + pts.push_back(s->source()); + pts.push_back(s->target()); + } + else if (t = std::get_if(&(i.first))) { + pts.push_back((*t)[0]); + pts.push_back((*t)[1]); + pts.push_back((*t)[2]); + auto pl = t->supporting_plane(); + for (std::size_t c = 0; c < 8; c++) + if (!taken[c]) + if (pl.oriented_side(corners[c]) != CGAL::ON_POSITIVE_SIDE) + taken[c] = true; + } + else if (v = std::get_if>(&(i.first))) { + const Triangle_3 &t = CGAL::Polygon_mesh_processing::triangle(i.second, tmesh); + auto pl = t.supporting_plane(); + for (std::size_t c = 0; c < 8; c++) + if (!taken[c]) + if (pl.oriented_side(corners[c]) != CGAL::ON_POSITIVE_SIDE) + taken[c] = true; + std::copy(v->begin(), v->end(), std::back_inserter(pts)); + } + } + + for (std::size_t c = 0; c < 8; c++) + if (taken[c]) + pts.push_back(corners[c]); + + pts.reserve(pts.size() + c.new_surface.size() * 8); + + for (const Vec3_uint& v : c.new_surface) { + FT xmin = bbox.xmin() + FT(v[0]) * voxel_size; + FT ymin = bbox.ymin() + FT(v[1]) * voxel_size; + FT zmin = bbox.zmin() + FT(v[2]) * voxel_size; + FT xmax = bbox.xmin() + FT(v[0] + 1) * voxel_size; + FT ymax = bbox.ymin() + FT(v[1] + 1) * voxel_size; + FT zmax = bbox.zmin() + FT(v[2] + 1) * voxel_size; + pts.push_back(Point_3(xmin, ymin, zmin)); + pts.push_back(Point_3(xmin, ymax, zmin)); + pts.push_back(Point_3(xmin, ymin, zmax)); + pts.push_back(Point_3(xmin, ymax, zmax)); + pts.push_back(Point_3(xmax, ymin, zmin)); + pts.push_back(Point_3(xmax, ymax, zmin)); + pts.push_back(Point_3(xmax, ymin, zmax)); + pts.push_back(Point_3(xmax, ymax, zmax)); + } + + convex_hull_3(pts.begin(), pts.end(), c.ch.points, c.ch.indices); + c.ch.bbox = bbox_3(c.ch.points.begin(), c.ch.points.end()); + + if (c.ch.indices.size() <= 3 || (c.ch.volume = volume(c.ch.points, c.ch.indices)) == 0) + c.ch.volume_error = -1; + else + c.ch.volume_error = CGAL::abs(c.ch.volume - c.ch.voxel_volume) / c.ch.voxel_volume; + } +} + +template +void shrink_candidates(const FaceGraph& tmesh, std::vector>& candidates, const Bbox_3& bbox, const typename GeomTraits::FT& voxel_size, CGAL::Parallel_tag) { +#ifdef CGAL_LINKED_WITH_TBB + using face_descriptor = typename boost::graph_traits::face_descriptor; + + using Point_3 = typename GeomTraits::Point_3; + using Segment_3 = typename GeomTraits::Segment_3; + using Triangle_3 = typename GeomTraits::Triangle_3; + using FT = typename GeomTraits::FT; + + using Primitive = CGAL::AABB_face_graph_triangle_primitive; + using Traits = CGAL::AABB_traits_3; + using Tree = CGAL::AABB_tree; + + Tree tree(faces(tmesh).first, faces(tmesh).second, tmesh); + + const auto shrink = [&tree, voxel_size, &bbox, &tmesh](Candidate& c) { + std::vector pts; + pts.reserve(c.new_surface.size() * 8); + using Box = typename GeomTraits::Iso_cuboid_3; + using IP_id = typename Tree::template Intersection_and_primitive_id::Type; + std::vector intersections; + tree.all_intersections(c.ch.bbox, std::back_inserter(intersections)); + + std::vector corners(8); + std::vector taken(8, false); + + corners[0] = Point_3(c.ch.bbox.xmin(), c.ch.bbox.ymin(), c.ch.bbox.zmin()); + corners[1] = Point_3(c.ch.bbox.xmin(), c.ch.bbox.ymax(), c.ch.bbox.zmin()); + corners[2] = Point_3(c.ch.bbox.xmin(), c.ch.bbox.ymin(), c.ch.bbox.zmax()); + corners[3] = Point_3(c.ch.bbox.xmin(), c.ch.bbox.ymax(), c.ch.bbox.zmax()); + corners[4] = Point_3(c.ch.bbox.xmax(), c.ch.bbox.ymin(), c.ch.bbox.zmin()); + corners[5] = Point_3(c.ch.bbox.xmax(), c.ch.bbox.ymax(), c.ch.bbox.zmin()); + corners[6] = Point_3(c.ch.bbox.xmax(), c.ch.bbox.ymin(), c.ch.bbox.zmax()); + corners[7] = Point_3(c.ch.bbox.xmax(), c.ch.bbox.ymax(), c.ch.bbox.zmax()); + + for (auto& i : intersections) { + const Point_3* p; + const Segment_3* s; + const Triangle_3* t; + const std::vector* v; + if (p = std::get_if(&(i.first))) { + pts.push_back(*p); + } + else if (s = std::get_if(&(i.first))) { + pts.push_back(s->source()); + pts.push_back(s->target()); + } + else if (t = std::get_if(&(i.first))) { + pts.push_back((*t)[0]); + pts.push_back((*t)[1]); + pts.push_back((*t)[2]); + auto pl = t->supporting_plane(); + for (std::size_t c = 0; c < 8; c++) + if (!taken[c]) + if (pl.oriented_side(corners[c]) != CGAL::ON_POSITIVE_SIDE) + taken[c] = true; + } + else if (v = std::get_if>(&(i.first))) { + const Triangle_3& t = CGAL::Polygon_mesh_processing::triangle(i.second, tmesh); + auto pl = t.supporting_plane(); + for (std::size_t c = 0; c < 8; c++) + if (!taken[c]) + if (pl.oriented_side(corners[c]) != CGAL::ON_POSITIVE_SIDE) + taken[c] = true; + std::copy(v->begin(), v->end(), std::back_inserter(pts)); + } + } + + for (std::size_t c = 0;c<8;c++) + if (taken[c]) + pts.push_back(corners[c]); + + pts.reserve(pts.size() + c.new_surface.size() * 8); + + for (const Vec3_uint& v : c.new_surface) { + FT xmin = bbox.xmin() + FT(v[0]) * voxel_size; + FT ymin = bbox.ymin() + FT(v[1]) * voxel_size; + FT zmin = bbox.zmin() + FT(v[2]) * voxel_size; + FT xmax = bbox.xmin() + FT(v[0] + 1) * voxel_size; + FT ymax = bbox.ymin() + FT(v[1] + 1) * voxel_size; + FT zmax = bbox.zmin() + FT(v[2] + 1) * voxel_size; + pts.push_back(Point_3(xmin, ymin, zmin)); + pts.push_back(Point_3(xmin, ymax, zmin)); + pts.push_back(Point_3(xmin, ymin, zmax)); + pts.push_back(Point_3(xmin, ymax, zmax)); + pts.push_back(Point_3(xmax, ymin, zmin)); + pts.push_back(Point_3(xmax, ymax, zmin)); + pts.push_back(Point_3(xmax, ymin, zmax)); + pts.push_back(Point_3(xmax, ymax, zmax)); + } + + convex_hull_3(pts.begin(), pts.end(), c.ch.points, c.ch.indices); + c.ch.bbox = bbox_3(c.ch.points.begin(), c.ch.points.end()); + + if (c.ch.indices.size() <= 3 || (c.ch.volume = volume(c.ch.points, c.ch.indices)) == 0) + c.ch.volume_error = -1; + else + c.ch.volume_error = CGAL::abs(c.ch.volume - c.ch.voxel_volume) / c.ch.voxel_volume; + }; + + tbb::parallel_for_each(candidates, shrink); +#else + assert(false); + CGAL_USE(candidates); + CGAL_USE(bbox); + CGAL_USE(voxel_size); +#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::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 { + if (volume_error == other.volume_error) + return other.ch > ch; + else + return (volume_error > other.volume_error); + } + + Merged_candidate() : ch_a(-1), ch_b(-1), ch(-1) {} + Merged_candidate(std::size_t ch_a, std::size_t ch_b) : ch_a(ch_a), ch_b(ch_b), ch(-1) {} + }; + + 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; + std::priority_queue queue; + + const auto do_merge = [hull_volume, &hulls, &num_hulls](Merged_candidate& m) { + const Convex_hull_candidate& ci = hulls[m.ch_a]; + const Convex_hull_candidate& cj = hulls[m.ch_b]; + + 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); + + if (ch.indices.size() <= 3) { + m.volume_error = ch.volume_error = -1; + return; + } + + ch.volume = volume(ch.points, ch.indices); + if (ci.volume_error == -1 || cj.volume_error == -1) { + m.volume_error = -1; + ch.volume_error = 0; + } + else + 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)); + todo.back().ch = num_hulls++; + } + else { + Merged_candidate m(i, j); + Bbox_3 bbox = box_union(ci.bbox, cj.bbox).bbox(); + + if (ci.volume_error == -1 || cj.volume_error == -1) + m.volume_error = -1; + else + 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 = 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) { + m.ch = num_hulls++; + 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)); + todo.back().ch = num_hulls++; + } + else { + Merged_candidate merged(id, m.ch); + Bbox_3 bbox = box_union(ci.bbox, cj.bbox).bbox(); + if (ci.volume_error == -1 || cj.volume_error == -1) + merged.volume_error = -1; + else + 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; + + candidates.reserve(max_convex_hulls); + + 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 { + if (volume_error == other.volume_error) + return other.ch > ch; + else + return (volume_error > other.volume_error); + } + + Merged_candidate() : ch_a(-1), ch_b(-1), ch(-1) {} + Merged_candidate(std::size_t ch_a, std::size_t ch_b) : ch_a(ch_a), ch_b(ch_b), ch(-1) {} + }; + + 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]; + + 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); + + if (ch.indices.size() <= 3) { + m.volume_error = ch.volume_error = -1; + return; + } + + ch.volume = volume(ch.points, ch.indices); + if (ci.volume_error == -1 || cj.volume_error == -1) { + m.volume_error = -1; + ch.volume_error = 0; + } + else + 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(); + + if (ci.volume_error == -1 || cj.volume_error == -1) + m.volume_error = -1; + else + 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) { + m.ch = num_hulls++; + 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(); + + if (ci.volume_error == -1 || cj.volume_error == -1) + merged.volume_error = -1; + else + 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 volumes. + * 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. + * In a next step, the convex hull of the mesh is hierarchically split until the `volume_error` threshold is satisfied or the `maximum_depth` is reached. + * Finally, a greedy pair-wise merging combines smaller convex volumes until `maximum_number_of_convex_volumes` 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 volumes + * \param out_volumes output iterator into which convex volumes 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{refitting} + * \cgalParamDescription{refitting of convex volumes after split phase} + * \cgalParamType{Boolean} + * \cgalParamDefault{true} + * \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 `tmesh`} + * \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, tmesh)`} + * \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 value type of the vertex-point map, 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_volumes`, for example if the specified `volume_error` is quickly met. + * + * \pre `tmesh` is a triangle mesh. + * \pre `tmesh` is not self-intersecting. + * \sa `CGAL::Polygon_mesh_processing::polygon_soup_to_polygon_mesh()` + */ +template +std::size_t approximate_convex_decomposition(const FaceGraph& tmesh, OutputIterator out_volumes, 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); + const bool refitting = parameters::choose_parameter(parameters::get_parameter(np, internal_np::refitting), true); + 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_volumes = parameters::choose_parameter(parameters::get_parameter(np, internal_np::maximum_number_of_convex_volumes), 16); + assert(max_convex_volumes > 0); + + if (max_convex_volumes == 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_volumes = 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); + internal::init(candidates[0], tmesh, grid, bb, grid_size, voxel_size, Concurrency_tag()); + + const FT hull_volume = candidates[0].ch.volume; + + internal::recurse(candidates, grid, grid_size, bb, voxel_size, np, Concurrency_tag()); + + if (refitting) + internal::shrink_candidates(tmesh, candidates, bb, voxel_size, Concurrency_tag()); + + std::vector> volumes; + for (internal::Candidate &c : candidates) + volumes.emplace_back(std::move(c.ch)); + + candidates.clear(); + + // merge until target number is reached + internal::merge(volumes, hull_volume, max_convex_volumes, Concurrency_tag()); + + for (std::size_t i = 0; i < volumes.size(); i++) + *out_volumes++ = std::make_pair(std::move(volumes[i].points), std::move(volumes[i].indices)); + + return volumes.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..652e2ada534 --- /dev/null +++ b/Surface_mesh_decomposition/package_info/Surface_mesh_decomposition/dependencies @@ -0,0 +1,30 @@ +AABB_tree +Algebraic_foundations +Arithmetic_kernel +BGL +CGAL_Core +Cartesian_kernel +Circulator +Convex_hull_2 +Convex_hull_3 +Distance_2 +Distance_3 +Filtered_kernel +Hash_map +Homogeneous_kernel +Installation +Intersections_2 +Intersections_3 +Interval_support +Kernel_23 +Kernel_d +Modular_arithmetic +Number_types +Polygon_mesh_processing +Profiling_tools +Property_map +Spatial_searching +STL_Extension +Stream_support +TDS_2 +Random_numbers 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..8bb8efcb72b --- /dev/null +++ b/Surface_mesh_decomposition/package_info/Surface_mesh_decomposition/license.txt @@ -0,0 +1 @@ +GPL (v3 or later) 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..3aea2f17c57 --- /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_Tests) + +# 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..9bea6abac6e --- /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_volumes; + + // Try with empty mesh + CGAL::approximate_convex_decomposition(mesh, std::back_inserter(convex_volumes), + CGAL::parameters::maximum_depth(10) + .volume_error(0.1) + .maximum_number_of_convex_volumes(9) + .split_at_concavity(true) + .maximum_number_of_voxels(1000000) + .concurrency_tag(CGAL::Parallel_if_available_tag())); + + assert(convex_volumes.size() == 0); + + if (!PMP::IO::read_polygon_mesh(filename, mesh)) { + std::cerr << "Invalid input." << std::endl; + return 1; + } + + convex_volumes.reserve(9); + + 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) + .split_at_concavity(true) + .maximum_number_of_voxels(1000000) + .concurrency_tag(CGAL::Parallel_if_available_tag())); + + assert(convex_volumes.size() > 0); + assert(convex_volumes.size() <= 9); + + return 0; +}