diff --git a/Installation/CHANGES.md b/Installation/CHANGES.md index b0ac14e4654c..19c2303a3a2e 100644 --- a/Installation/CHANGES.md +++ b/Installation/CHANGES.md @@ -26,6 +26,9 @@ Release date: December 2026 - The corefinement based operations (including Boolean operations) has been optimized to better handle cases when some identical faces are shared between the input meshes. This leads to a significant speed up in those cases. +- Added new multithreading compatible functions operating directly on triangle soups, enabling the handling of non-manifold input and output volumes: + `CGAL::Polygon_mesh_processing::compute_union()`, `CGAL::Polygon_mesh_processing::compute_intersection()`, and + `CGAL::Polygon_mesh_processing::compute_difference()`. ## [Release 6.2](https://github.com/CGAL/cgal/releases/tag/v6.2) diff --git a/PMP_Boolean_operations/doc/PMP_Boolean_operations/Doxyfile.in b/PMP_Boolean_operations/doc/PMP_Boolean_operations/Doxyfile.in index ef1b78846943..eaf44efac471 100644 --- a/PMP_Boolean_operations/doc/PMP_Boolean_operations/Doxyfile.in +++ b/PMP_Boolean_operations/doc/PMP_Boolean_operations/Doxyfile.in @@ -8,5 +8,3 @@ HIDE_UNDOC_CLASSES = true WARN_IF_UNDOCUMENTED = false EXCLUDE_SYMBOLS += experimental - -HTML_EXTRA_FILES = ${CGAL_PACKAGE_DOC_DIR}/fig/kernel.png diff --git a/PMP_Boolean_operations/doc/PMP_Boolean_operations/PMP_Boolean_operations.txt b/PMP_Boolean_operations/doc/PMP_Boolean_operations/PMP_Boolean_operations.txt index 8da42ef9df0e..164f3575ef37 100644 --- a/PMP_Boolean_operations/doc/PMP_Boolean_operations/PMP_Boolean_operations.txt +++ b/PMP_Boolean_operations/doc/PMP_Boolean_operations/PMP_Boolean_operations.txt @@ -18,7 +18,7 @@ of the geometry processing toolkit, with numerous applications in CAD and computer graphics. As an example, it constitutes the essential operation of Constructive Solid Geometry (CSG). This \cgal package offers functions to compute the union, intersection, and difference -of volumes bounded by triangle meshes. The core of the method is the corefinement technique, +of volumes bounded by triangle meshes and triangle soups. The core of the method is the corefinement technique, which refines two meshes such that they become conformal. This technique is re-used to also propose functions that enable trimming a mesh ("clipping") and bisecting a mesh ("splitting") with a plane, a box, or another mesh. @@ -145,7 +145,7 @@ possible, the input mesh will still be corefined. If the input meshes contain triangles that are coplanar, it is guaranteed that the triangulations of those triangles in the two corefined meshes will be identical. Also all edges of those refined triangles will be reported as intersection -edges (both in the function `CGAL::Polygon_mesh_processing::surface_intersection()` and in the named parameter +edges (both in the function `CGAL::Polygon_mesh_processing::intersection_polylines()` and in the named parameter `edge_is_constrained_map()` of corefine prefixed functions). \subsubsection coref_bolop_input_autorefsub Autorefinement @@ -154,6 +154,22 @@ There is no restriction on the input for the function `CGAL::Polygon_mesh_proces Degenerate triangles will be eliminated, duplicated points will be merged. If several triangles are coplanar, the output will contain as many copies of the triangulation of those refined triangles. +\subsection tsoup_bolop_subsec Computing Boolean Operations on Triangle Soups +The Boolean operation functions based on corefine operate on triangle meshes. This data structure is not meant to +handle non-manifold edges. As it is sometime desirable to handle such cases, we also provide operations that are +taking as input triangle soups, representing possibly non-manifold volumes, and providing the result of a Boolean +operation as a triangle soup. In the following example, we use the function `CGAL::Polygon_mesh_processing::compute_difference()` +to iteratively subtract pieces of a puzzle. + +\cgalFigureAnchor{nm_no_fig} +
+ +
+\cgalFigureCaptionBegin{nm_no_fig} +Iterative subtraction of the pieces of a puzzle from David's torso. The result of the second operation is non-manifold and +so are all subsequent operations, until the last one. +\cgalFigureCaptionEnd + \subsection coref_valid_subsec Kernel Choice and Output Validity @@ -323,5 +339,7 @@ The kernel functions and the specialization of clipping and plane refinement to inputs were implemented in 2025 by Léo Valque, inspired by the work of Sorgente et al. \cgalCite{sorgente2022polyhedron}. +The illustrations of the package pay homage to French sculptor and painter Sacha Sosno. + */ } /* namespace CGAL */ diff --git a/PMP_Boolean_operations/doc/PMP_Boolean_operations/PackageDescription.txt b/PMP_Boolean_operations/doc/PMP_Boolean_operations/PackageDescription.txt index 0e9d4cd3af18..2f2252fce1a5 100644 --- a/PMP_Boolean_operations/doc/PMP_Boolean_operations/PackageDescription.txt +++ b/PMP_Boolean_operations/doc/PMP_Boolean_operations/PackageDescription.txt @@ -6,8 +6,12 @@ /// Functions to corefine triangulated surface meshes and triangle soups. /// \ingroup PkgPMPBooleanOperationsRef -/// \defgroup PMP_boolop_grp Boolean Operations -/// Functions to compute triangulated the union, difference, and intersection of meshes. +/// \defgroup PMP_boolop_grp Boolean Operations on Triangle Meshes +/// Functions to compute the union, difference, and intersection of triangle meshes. +/// \ingroup PkgPMPBooleanOperationsRef + +/// \defgroup PMP_boolop_soup_grp Boolean Operations on Triangle Soups +/// Functions to compute triangulated the union, difference, and intersection of triangle soups (possibly non-manifold). /// \ingroup PkgPMPBooleanOperationsRef /// \defgroup PMP_clip_grp Clipping and Splitting @@ -66,6 +70,11 @@ The page \ref bgl_namedparameters "Named Parameters" describes their usage. - `CGAL::Polygon_mesh_processing::surface_intersection()` (deprecated) - `CGAL::Polygon_mesh_processing::intersection_polylines()` +\cgalCRPSection{General Boolean Operations} +- `CGAL::Polygon_mesh_processing::compute_union()` +- `CGAL::Polygon_mesh_processing::compute_difference()` +- `CGAL::Polygon_mesh_processing::compute_intersection()` + \cgalCRPSection{Clipping and Splitting} - `CGAL::Polygon_mesh_processing::clip()` - `CGAL::Polygon_mesh_processing::split()` diff --git a/PMP_Boolean_operations/doc/PMP_Boolean_operations/fig/i_nm_bo.png b/PMP_Boolean_operations/doc/PMP_Boolean_operations/fig/i_nm_bo.png new file mode 100644 index 000000000000..6d650f92c82f Binary files /dev/null and b/PMP_Boolean_operations/doc/PMP_Boolean_operations/fig/i_nm_bo.png differ diff --git a/PMP_Boolean_operations/examples/PMP_Boolean_operations/CMakeLists.txt b/PMP_Boolean_operations/examples/PMP_Boolean_operations/CMakeLists.txt index ae035b0858de..1f6e746c1f03 100644 --- a/PMP_Boolean_operations/examples/PMP_Boolean_operations/CMakeLists.txt +++ b/PMP_Boolean_operations/examples/PMP_Boolean_operations/CMakeLists.txt @@ -18,6 +18,7 @@ create_single_source_cgal_program("corefinement_mesh_union_with_attributes.cpp") create_single_source_cgal_program("corefinement_polyhedron_union.cpp") create_single_source_cgal_program("corefinement_LCC.cpp") create_single_source_cgal_program("soup_autorefinement.cpp") +create_single_source_cgal_program("soup_union.cpp") create_single_source_cgal_program("surface_mesh_intersection.cpp") create_single_source_cgal_program("triangle_mesh_autorefinement.cpp") create_single_source_cgal_program("snap_polygon_soup.cpp") @@ -49,6 +50,7 @@ find_package(TBB QUIET) include(CGAL_TBB_support) if(TARGET CGAL::TBB_support) target_link_libraries(soup_autorefinement PRIVATE CGAL::TBB_support) + target_link_libraries(soup_union PRIVATE CGAL::TBB_support) target_link_libraries(snap_polygon_soup PRIVATE CGAL::TBB_support) create_single_source_cgal_program("corefinement_parallel_union_meshes.cpp") diff --git a/PMP_Boolean_operations/examples/PMP_Boolean_operations/soup_union.cpp b/PMP_Boolean_operations/examples/PMP_Boolean_operations/soup_union.cpp new file mode 100644 index 000000000000..b05f0899bbf8 --- /dev/null +++ b/PMP_Boolean_operations/examples/PMP_Boolean_operations/soup_union.cpp @@ -0,0 +1,52 @@ +#include + +#include +#include + +#include +#include +#include + +typedef CGAL::Exact_predicates_exact_constructions_kernel K; +typedef K::Point_3 Point; + +namespace PMP = CGAL::Polygon_mesh_processing; + +int main(int argc, char* argv[]) +{ + const std::string filename1 = + (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/blobby.off"); + const std::string filename2 = + (argc > 2) ? argv[2] : CGAL::data_file_path("meshes/eight.off"); + + std::vector points1, points2; + std::vector> triangles1, triangles2; + + if(!CGAL::IO::read_polygon_soup(filename1, points1, triangles1)) + { + std::cerr << "Failed to read file: " << filename1 << std::endl; + return 1; + } + + if(!CGAL::IO::read_polygon_soup(filename2, points2, triangles2)) + { + std::cerr << "Failed to read file: " << filename2 << std::endl; + return 1; + } + + std::vector points_output; + std::vector> triangles_output; + + PMP::compute_union(points1, triangles1, + points2, triangles2, + points_output, triangles_output, + CGAL::parameters::concurrency_tag(CGAL::Parallel_if_available_tag())); + + std::cout << "Union was successfully computed\n"; + + CGAL::IO::write_polygon_soup( + "union.off", points_output, triangles_output, + CGAL::parameters::stream_precision(17)); + + return 0; +} \ No newline at end of file diff --git a/PMP_Boolean_operations/include/CGAL/Polygon_mesh_processing/autorefinement.h b/PMP_Boolean_operations/include/CGAL/Polygon_mesh_processing/autorefinement.h index 1340dc5aeff8..6801f7265468 100644 --- a/PMP_Boolean_operations/include/CGAL/Polygon_mesh_processing/autorefinement.h +++ b/PMP_Boolean_operations/include/CGAL/Polygon_mesh_processing/autorefinement.h @@ -65,10 +65,6 @@ //#define CGAL_AUTOREFINE_DEBUG_COUNTERS //#define CGAL_AUTOREF_DEBUG_DEPTH -#ifdef CGAL_AUTOREF_USE_DEBUG_PARALLEL_TIMERS -#include -#endif - #if defined(CGAL_AUTOREFINE_DEBUG_COUNTERS) || defined(CGAL_AUTOREF_USE_DEBUG_PARALLEL_TIMERS) #include #endif @@ -977,6 +973,75 @@ void generate_subtriangles(std::size_t ti, #endif namespace autorefine_impl{ + +// +template +void compute_intersecting_pairs(const PointRange& points, const PolygonRange& triangles, + std::size_t split_id, std::vector& si_pairs, + const PointMap& pm) +{ + typedef CGAL::Box_intersection_d::ID_FROM_BOX_ADDRESS Box_policy; + typedef CGAL::Box_intersection_d::Box_with_info_d Box; + + std::vector boxes_1, boxes_2; + std::vector boxes_1_ptr, boxes_2_ptr; + + boxes_1.reserve(split_id); + boxes_1_ptr.reserve(boxes_1.size()); + for (std::size_t i=0;i; + using TBB_vector = tbb::concurrent_vector; + std::conditional_t > collecter; +#else + std::vector& collecter=si_pairs; +#endif + auto callback = [&points, &triangles, &collecter, &pm] + (const Box* b1, const Box* b2) + { + std::size_t i1=b1->info(), i2=b2->info(); + + if (do_intersect(typename Kernel::Triangle_3(get(pm, points[triangles[i1][0]]),get(pm, points[triangles[i1][1]]),get(pm, points[triangles[i1][2]])), + typename Kernel::Triangle_3(get(pm, points[triangles[i2][0]]),get(pm, points[triangles[i2][1]]),get(pm, points[triangles[i2][2]])))) + { + collecter.emplace_back(i1,i2); + } + }; + + //using pointers in box_intersection_d is about 10% faster + CGAL::box_intersection_d(boxes_1_ptr.begin(), boxes_1_ptr.end(), + boxes_2_ptr.begin(), boxes_2_ptr.end(), + callback, cutoff ); + +#ifdef CGAL_LINKED_WITH_TBB + if constexpr (need_copy) + { + si_pairs.reserve(collecter.size()); + std::copy(collecter.begin(), collecter.end(), std::back_inserter(si_pairs)); + } + else + si_pairs.swap(collecter); +#endif +} + + // Forward declaration struct Wrap_snap_visitor; @@ -992,6 +1057,7 @@ bool autorefine_triangle_soup(PointRange& soup_points, { using parameters::choose_parameter; using parameters::get_parameter; + using parameters::is_default_parameter; typedef typename GetPolygonSoupGeomTraits::type GT; typedef typename GetPointMap::const_type Point_map; @@ -1027,17 +1093,27 @@ bool autorefine_triangle_soup(PointRange& soup_points, // collect intersecting pairs of triangles CGAL_PMP_AUTOREFINE_VERBOSE("collect intersecting pairs"); -#ifdef CGAL_AUTOREF_USE_DEBUG_PARALLEL_TIMERS - Real_timer t; - t.start(); -#endif - triangle_soup_self_intersections(soup_points, soup_triangles, std::back_inserter(si_pairs), np); -#ifdef CGAL_AUTOREF_USE_DEBUG_PARALLEL_TIMERS + +#if defined(CGAL_AUTOREFINE_DEBUG_COUNTERS) || defined(CGAL_AUTOREF_USE_DEBUG_PARALLEL_TIMERS) std::string mode = parallel_execution ? "parallel" : "sequential"; - t.stop(); - std::cout << t.time() << " sec. for #0 (" << mode << ")" << std::endl; - t.reset(); + CGAL::Real_timer timer; + timer.start(); +#endif + + if constexpr(!parameters::is_default_parameter::value) + { + std::size_t split_id=get_parameter(np, internal_np::split_triangle_range_at); + compute_intersecting_pairs(soup_points, soup_triangles, split_id, si_pairs, pm); + } + else + triangle_soup_self_intersections(soup_points, soup_triangles, std::back_inserter(si_pairs), np); + +#if defined(CGAL_AUTOREFINE_DEBUG_COUNTERS) || defined(CGAL_AUTOREF_USE_DEBUG_PARALLEL_TIMERS) + timer.stop(); + std::cout << "collecting intersecting pairs took (" << mode << ") " << timer.time() << "\n"; + timer.reset(); #endif + if (si_pairs.empty()) { if constexpr (!std::is_same_v) @@ -1103,7 +1179,7 @@ bool autorefine_triangle_soup(PointRange& soup_points, CGAL_PMP_AUTOREFINE_VERBOSE("compute intersections ("+std::to_string(si_pairs.size())+")"); #ifdef CGAL_AUTOREF_USE_DEBUG_PARALLEL_TIMERS - t.start(); + timer.start(); #endif std::set > intersecting_triangles; std::set > coplanar_triangles; @@ -1192,13 +1268,13 @@ bool autorefine_triangle_soup(PointRange& soup_points, } } #ifdef CGAL_AUTOREF_USE_DEBUG_PARALLEL_TIMERS - t.stop(); - std::cout << t.time() << " sec. for #2" << std::endl; - t.reset(); + timer.stop(); + std::cout << timer.time() << " sec. for #2" << std::endl; + timer.reset(); #endif #ifdef CGAL_AUTOREF_USE_DEBUG_PARALLEL_TIMERS - t.start(); + timer.start(); #endif // deduplicate inserted segments @@ -1304,9 +1380,9 @@ bool autorefine_triangle_soup(PointRange& soup_points, } #ifdef CGAL_AUTOREF_USE_DEBUG_PARALLEL_TIMERS - t.stop(); - std::cout << t.time() << " sec. for #3" << std::endl; - t.reset(); + timer.stop(); + std::cout << timer.time() << " sec. for #3" << std::endl; + timer.reset(); #endif CGAL_PMP_AUTOREFINE_VERBOSE("triangulate faces ("+std::to_string(triangles.size())+")"); @@ -1355,7 +1431,7 @@ bool autorefine_triangle_soup(PointRange& soup_points, }; #ifdef CGAL_AUTOREF_USE_DEBUG_PARALLEL_TIMERS - t.start(); + timer.start(); #endif #ifdef CGAL_LINKED_WITH_TBB if (parallel_execution) @@ -1374,9 +1450,10 @@ bool autorefine_triangle_soup(PointRange& soup_points, } #ifdef CGAL_AUTOREF_USE_DEBUG_PARALLEL_TIMERS - t.stop(); - std::cout << t.time() << " sec. for #1" << std::endl; - t.reset(); + timer.stop(); + std::cout << timer.time() << " sec. for #1" << std::endl; + timer.reset(); + timer.start(); #endif // brute force output: create a soup, orient and to-mesh @@ -1394,23 +1471,55 @@ bool autorefine_triangle_soup(PointRange& soup_points, Point_id_map point_id_map; #if ! defined(CGAL_NDEBUG) || defined(CGAL_DEBUG_PMP_AUTOREFINE) - std::vector exact_soup_points; +#ifdef CGAL_LINKED_WITH_TBB + typedef std::conditional_t, + std::vector> Exact_point_range; +#else + typedef std::vector Exact_point_range; +#endif + Exact_point_range exact_soup_points; +#endif + + +#if ! defined(CGAL_NDEBUG) || defined(CGAL_DEBUG_PMP_AUTOREFINE) + exact_soup_points.resize(soup_points.size()); #endif - // TODO: parallel_for? // for input points, we on purpose keep duplicated points and isolated points - for (std::size_t pid = 0; pidfirst); + exact_soup_points[pid]=insert_res.first->first; +#endif + }; +#ifdef CGAL_LINKED_WITH_TBB + if constexpr (parallel_execution) + { + tbb::parallel_for(tbb::blocked_range(0, soup_points.size()), + [&](const tbb::blocked_range& r) { + for (std::size_t pid = r.begin(); pid != r.end(); ++pid) + insert_point_id(pid); + }); + } + else #endif + for (std::size_t pid = 0; pid(exact_soup_points, soup_triangles_out) ); + // If split_triangle_range is defined, intersection can still be present inside one of the range + if constexpr(parameters::is_default_parameter::value){ + CGAL_assertion( !does_triangle_soup_self_intersect(exact_soup_points, soup_triangles_out) ); + } #else #ifdef CGAL_DEBUG_PMP_AUTOREFINE CGAL_PMP_AUTOREFINE_VERBOSE("check soup"); @@ -1757,10 +1875,12 @@ bool autorefine_triangle_soup(PointRange& soup_points, /** * \ingroup PMP_corefinement_grp - * refines a triangle mesh so that no triangles intersects in their interior. + * refines a triangle mesh to ensure no pair of triangles intersects, except at a shared edge or vertex. + * If the input does self-intersect, then the intersection will be resolved in the output but the + * mesh will still be self-intersecting. * * Note that this function is only provided as a shortcut for calling `autorefine_triangle_soup()` - * with a mesh. For any advance usage the aforementioned function should be called directly. + * with a mesh. For any advanced usage the aforementioned function should be called directly. * * @tparam TriangleMesh a model of `HalfedgeListGraph`, `FaceListGraph`, and `MutableFaceGraph` * @tparam NamedParameters a sequence of \ref bgl_namedparameters "Named Parameters" diff --git a/PMP_Boolean_operations/include/CGAL/Polygon_mesh_processing/internal/Corefinement/Visitor.h b/PMP_Boolean_operations/include/CGAL/Polygon_mesh_processing/internal/Corefinement/Visitor.h index 3794336dad27..19bdcc5a7d6f 100644 --- a/PMP_Boolean_operations/include/CGAL/Polygon_mesh_processing/internal/Corefinement/Visitor.h +++ b/PMP_Boolean_operations/include/CGAL/Polygon_mesh_processing/internal/Corefinement/Visitor.h @@ -1344,7 +1344,18 @@ class Surface_intersection_visitor_for_corefinement{ q = nodes.to_exact(get(vpm,f_vertices[1])), r = nodes.to_exact(get(vpm,f_vertices[2])); ///TODO use a positive normal and remove all workaround to guarantee that triangulation of coplanar patches are compatible - CDT_traits traits(typename EK::Construct_normal_3()(p,q,r)); + typename EK::Vector_3 n = typename EK::Construct_normal_3()(p,q,r); + typename EK::Point_3 o(CGAL::ORIGIN); + + bool orientation_flipped = false; + if ( typename EK::Less_xyz_3()(o+n,o) ) + { + n=-n; + orientation_flipped = true; + } + + + CDT_traits traits(n); CDT cdt(traits); // insert triangle points @@ -1352,7 +1363,7 @@ class Surface_intersection_visitor_for_corefinement{ //we can do this to_exact because these are supposed to be input points. triangle_vertices[0]=cdt.insert_outside_affine_hull(p); triangle_vertices[1]=cdt.insert_outside_affine_hull(q); - triangle_vertices[2]=cdt.tds().insert_dim_up(cdt.infinite_vertex(), false); + triangle_vertices[2]=cdt.tds().insert_dim_up(cdt.infinite_vertex(), orientation_flipped); triangle_vertices[2]->set_point(r); triangle_vertices[0]->info()=f_indices[0]; @@ -1378,6 +1389,7 @@ class Surface_intersection_visitor_for_corefinement{ // is already tight and the call in Intersection_nodes::finalize() will not fix anything } } + //insert points on edges if (it_fb!=face_boundaries.end()) //if f not a triangle? { @@ -1386,9 +1398,18 @@ class Surface_intersection_visitor_for_corefinement{ for (int i=0;i<3;++i) { int oi=-1; - CGAL_assertion_code(bool is_edge = ) - cdt.is_edge(triangle_vertices[i], triangle_vertices[(i+1)%3], infinite_faces[i], oi); - CGAL_assertion(is_edge); + if (orientation_flipped) + { + CGAL_assertion_code(bool is_edge = ) + cdt.is_edge(triangle_vertices[(i+1)%3], triangle_vertices[i], infinite_faces[i], oi); + CGAL_assertion(is_edge); + } + else + { + CGAL_assertion_code(bool is_edge = ) + cdt.is_edge(triangle_vertices[i], triangle_vertices[(i+1)%3], infinite_faces[i], oi); + CGAL_assertion(is_edge); + } CGAL_assertion( cdt.is_infinite( infinite_faces[i]->vertex(oi) ) ); } @@ -1400,22 +1421,30 @@ class Surface_intersection_visitor_for_corefinement{ //handle case of halfedge starting at triangle_vertices[i] // and ending at triangle_vertices[(i+1)%3] - const Node_ids& ids_on_edge=f_boundary.node_ids_array[i]; + Node_ids& ids_on_edge=f_boundary.node_ids_array[i]; CDT_Vertex_handle previous=triangle_vertices[i]; - Node_id prev_index=f_indices[i];// node-id of the mesh vertex + Node_id prev_index=orientation_flipped?f_indices[(i+1)%3]:f_indices[i];// node-id of the mesh vertex halfedge_descriptor hedge = next(f_boundary.halfedges[(i+2)%3],tm); - CGAL_assertion( source(hedge,tm)==f_boundary.vertices[i] ); + if (orientation_flipped) + for (std::size_t l=0;lvertex(cw(infinite_vertex)) to fh->vertex(ccw(infinite_vertex)) + if (orientation_flipped) + std::reverse(ids_on_edge.begin(), ids_on_edge.end()); for(Node_id id : ids_on_edge) { CDT_Vertex_handle vh=insert_point_on_ch_edge(cdt,infinite_faces[i],nodes.exact_node(id)); vh->info()=id; id_to_CDT_vh.insert(std::make_pair(id,vh)); - edge_to_hedge[std::make_pair(prev_index,id)]=hedge; + auto key = orientation_flipped ? std::make_pair(id, prev_index) + : std::make_pair(prev_index,id); + edge_to_hedge[key]=hedge; previous=vh; - hedge=next(hedge,tm); + hedge=orientation_flipped?prev(hedge,tm):next(hedge,tm); prev_index=id; } } @@ -1424,9 +1453,11 @@ class Surface_intersection_visitor_for_corefinement{ CGAL_assertion( target(hd,tm) == f_boundary.vertices[(i+1)%3] ); CGAL_assertion( source(hd,tm) == f_boundary.vertices[ i ] ); } - CGAL_assertion(hedge==f_boundary.halfedges[i]); - edge_to_hedge[std::make_pair(prev_index,f_indices[(i+1)%3])] = - it_fb->second.halfedges[i]; + CGAL_assertion(orientation_flipped || hedge==f_boundary.halfedges[i]); + CGAL_assertion(!orientation_flipped || hedge==next(f_boundary.halfedges[(i+2)%3],tm)); + auto key = orientation_flipped ? std::make_pair(f_indices[i], prev_index) + : std::make_pair(prev_index,f_indices[(i+1)%3]); + edge_to_hedge[key] = hedge/* it_fb->second.halfedges[i] */; } } @@ -1507,7 +1538,7 @@ class Surface_intersection_visitor_for_corefinement{ // import the triangle in `cdt` in the face `f` of `tm` triangulate_a_face(f, tm, nodes, node_ids, node_id_to_vertex, - edge_to_hedge, cdt, vpm, output_builder, user_visitor); + edge_to_hedge, cdt, vpm, orientation_flipped, output_builder, user_visitor); // TODO Here we do the update only for internal edges. // Update for border halfedges could be done during the split diff --git a/PMP_Boolean_operations/include/CGAL/Polygon_mesh_processing/internal/Corefinement/face_graph_utils.h b/PMP_Boolean_operations/include/CGAL/Polygon_mesh_processing/internal/Corefinement/face_graph_utils.h index 4276641ce7c7..e65b13ad5b11 100644 --- a/PMP_Boolean_operations/include/CGAL/Polygon_mesh_processing/internal/Corefinement/face_graph_utils.h +++ b/PMP_Boolean_operations/include/CGAL/Polygon_mesh_processing/internal/Corefinement/face_graph_utils.h @@ -526,6 +526,7 @@ triangulate_a_face( ::halfedge_descriptor>& edge_to_hedge, const CDT& cdt, const VertexPointMap& vpm, + bool orientation_flipped, OutputBuilder& output_builder, UserVisitor& user_visitor) { @@ -592,6 +593,9 @@ triangulate_a_face( typename CDT::Vertex_handle cdt_v1=it->vertex(1); typename CDT::Vertex_handle cdt_v2=it->vertex(2); + if (orientation_flipped) + std::swap(cdt_v0, cdt_v2); + Node_id i0=cdt_v0->info(), i1=cdt_v1->info(), i2=cdt_v2->info(); CGAL_assertion(edge_to_hedge.count(std::make_pair(i0,i1))!= 0); diff --git a/PMP_Boolean_operations/include/CGAL/Polygon_mesh_processing/triangle_soup_boolean_operations_3.h b/PMP_Boolean_operations/include/CGAL/Polygon_mesh_processing/triangle_soup_boolean_operations_3.h new file mode 100644 index 000000000000..7311f36e4354 --- /dev/null +++ b/PMP_Boolean_operations/include/CGAL/Polygon_mesh_processing/triangle_soup_boolean_operations_3.h @@ -0,0 +1,1359 @@ +// Copyright (c) 2026 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) : Sébastien Loriot + +#ifndef CGAL_POLYGON_MESH_PROCESSING_TRIANGLE_SOUP_BOOLEAN_OPERATIONS_3_H +#define CGAL_POLYGON_MESH_PROCESSING_TRIANGLE_SOUP_BOOLEAN_OPERATIONS_3_H + +#include + +#include +#include +#include +#include + +#include + +#include + +#ifdef CGAL_LINKED_WITH_TBB +#include +#endif + +#if defined(CGAL_BO3_TIMERS) || defined(CGAL_BO3_SELECT_TIMERS) +#include +#endif + +namespace CGAL { +namespace boolops_3 { + +//TODO: update nm_edges for identical coplanar status (less test below then) +template +std::vector +connected_components(std::size_t begin, + std::size_t end, + int offset, + int increment, + const EdgeConnectivity& adjacent_edges, + const FaceConnectivity& adjacent_faces, + std::vector& ccids, + NM_edges& nm_edges) +{ + auto is_handled = [&](std::size_t tid){ return ccids[tid]!=std::size_t(-1); }; + std::vector one_triangle_per_cc; + + std::size_t i=offset; + + std::vector queue; + for (std::size_t tid=begin; tid 0); + switch( tids.size() ) + { + case 1: + // boundary edge + break; + case 2: + { + std::size_t n_tid = tids[0]!=q_tid ? tids[0] : tids[1]; + if (!is_handled(n_tid)){ + ccids[n_tid]=i; + queue.push_back(n_tid); + } + } + break; + default: + // non-manifold edge + nm_edges.emplace_back(ei); + break; + } + } + } + one_triangle_per_cc.push_back(tid); + i += increment; + } + return one_triangle_per_cc; +} + +template +std::vector +connected_components(const EdgeConnectivity& adjacent_edges, + const FaceConnectivity& adjacent_faces, + std::vector& ccids, + NM_edges& nm_edges) +{ + std::size_t nb_poly = adjacent_edges.size(); + ccids.resize(nb_poly, std::size_t(-1)); + return connected_components(0, nb_poly, 0, 1, adjacent_edges, adjacent_faces, ccids, nm_edges); +} + +#ifdef CGAL_LINKED_WITH_TBB +template +std::vector +parallel_connected_components(const size_t split_range_at, + const size_t new_triangles_at, + const EdgeConnectivity& adjacent_edges, + const FaceConnectivity& adjacent_faces, + std::vector& ccids, + NM_edges& nm_edges) +{ + std::size_t nb_poly = adjacent_edges.size(); + ccids.resize(nb_poly, std::size_t(-1)); + std::vector trs_per_cc_1, trs_per_cc_2, trs_per_cc_end, merged_trs_per_cc; + + tbb::task_group task_group; + task_group.run([&](){ + trs_per_cc_1 = connected_components(0, split_range_at, 0, 2, adjacent_edges, adjacent_faces, ccids, nm_edges); + }); + NM_edges nm_edges_cc_2; + trs_per_cc_2 = connected_components(split_range_at, new_triangles_at, 1, 2, adjacent_edges, adjacent_faces, ccids, nm_edges_cc_2); + task_group.wait(); + nm_edges.insert(nm_edges.end(), nm_edges_cc_2.begin(), nm_edges_cc_2.end()); + std::size_t max_size = (std::max)(trs_per_cc_1.size(), trs_per_cc_2.size()); + trs_per_cc_end = connected_components(new_triangles_at, nb_poly, max_size, 1, adjacent_edges, adjacent_faces, ccids, nm_edges); + + merged_trs_per_cc.resize(2 * max_size + trs_per_cc_end.size(), std::size_t(-1)); + for(std::size_t i=0; i ccids; + std::vector coplanar_patches; + std::vector coplanar_patches_for_union_and_intersection; + std::vector< std::vector > cc_is_inside_per_mesh; // TODO: use an optional + std::vector input_ids; +}; + +template +Selection_data +select(const std::vector& points, + const std::vector& triangles, + std::vector& input_ids, +#ifdef CGAL_LINKED_WITH_TBB + std::size_t split_at, + std::size_t new_triangles_at, +#else + std::size_t /* split_at */, + std::size_t /* new_triangles_at */, +#endif + const unsigned int nb_meshes) +{ +#ifdef CGAL_BO3_SELECT_TIMERS + Real_timer t; + t.start(); +#endif + + namespace PMP = CGAL::Polygon_mesh_processing; + namespace pred = PMP::Corefinement; + Selection_data data_out; + + // Create connectivity information + using PID = std::size_t; + using EID = std::size_t; + using TID = std::size_t; + + using Edge_record = std::tuple; + std::vector partial_edges; + partial_edges.resize(3 * triangles.size()); +#ifdef CGAL_LINKED_WITH_TBB + if (std::is_same_v) + { + tbb::parallel_for(std::size_t(0), triangles.size(), [&](std::size_t tid) + { + auto add_edge = [&](std::size_t idx, PID a, PID b){ + if (b < a) std::swap(a, b); + partial_edges[idx] = Edge_record{a, b, tid}; + }; + add_edge(3 * tid + 0, triangles[tid][0], triangles[tid][1]); + add_edge(3 * tid + 1, triangles[tid][1], triangles[tid][2]); + add_edge(3 * tid + 2, triangles[tid][2], triangles[tid][0]); + }); + } + else +#endif + { + for (TID tid = 0; tid < triangles.size(); ++tid) + { + auto add_edge = [&](std::size_t idx, PID a, PID b){ + if (b < a) std::swap(a, b); + partial_edges[idx] = Edge_record{a, b, tid}; + }; + add_edge(3 * tid + 0, triangles[tid][0], triangles[tid][1]); + add_edge(3 * tid + 1, triangles[tid][1], triangles[tid][2]); + add_edge(3 * tid + 2, triangles[tid][2], triangles[tid][0]); + } + } + +#ifdef CGAL_BO3_SELECT_TIMERS + std::cout << t.time() << " sec. to collect edge records\n"; +#endif + +#ifdef CGAL_LINKED_WITH_TBB + if (std::is_same_v) + { + tbb::parallel_sort(partial_edges.begin(), partial_edges.end()); + } + else +#endif + { + std::sort(partial_edges.begin(), partial_edges.end()); + } + +#ifdef CGAL_BO3_SELECT_TIMERS + std::cout << t.time() << " sec. to sort edge records\n"; +#endif + + // Build adjacencies between edges and faces + std::vector> adjacent_edges(triangles.size()); + std::vector> adjacent_faces; + +#ifdef CGAL_LINKED_WITH_TBB + if constexpr(std::is_same_v) + { + // Cut partial_edges in blocks and compute the offset on edges for each block + const std::size_t block_size = 4096; + const std::size_t nb_blocks = (partial_edges.size() + block_size - 1) / block_size; + + std::vector counts(nb_blocks, 0); + + // count edge groups whose first record belongs to each block. + tbb::parallel_for(tbb::blocked_range(0, nb_blocks, 1), [&](const tbb::blocked_range& r) + { + for (std::size_t bi = r.begin(); bi != r.end(); ++bi){ + const std::size_t begin = bi * block_size; + const std::size_t end = (std::min)(begin + block_size, partial_edges.size()); + + for (std::size_t i = begin; i < end; ++i) + if (i == 0 || std::get<0>(partial_edges[i]) != std::get<0>(partial_edges[i-1]) || + std::get<1>(partial_edges[i]) != std::get<1>(partial_edges[i-1])) + ++counts[bi]; + } + }); + + // Compute offset for each block + std::vector offsets(nb_blocks+1, 0); + for (std::size_t i=0; i < nb_blocks; ++i) + offsets[i+1] = offsets[i] + counts[i]; + + adjacent_faces.resize(offsets.back()); + + tbb::parallel_for(tbb::blocked_range(0, nb_blocks, 1), [&](const tbb::blocked_range& r) + { + for (std::size_t bi = r.begin(); bi != r.end(); ++bi) + { + std::size_t begin = bi * block_size; + const std::size_t end = (std::min)(begin + block_size, partial_edges.size()); + std::size_t idx = offsets[bi]; + + // Skip a group that started in the previous block. + if(begin != 0) + while(std::get<0>(partial_edges[begin-1]) == std::get<0>(partial_edges[begin]) && + std::get<1>(partial_edges[begin-1]) == std::get<1>(partial_edges[begin])) + ++begin; + + for (std::size_t i = begin; i < end; ) + { + const PID a = std::get<0>(partial_edges[i]); + const PID b = std::get<1>(partial_edges[i]); + + auto& faces = adjacent_faces[idx]; + CGAL_assertion(faces.empty()); + + // may extend beyond the block's end, but is processed only once. + do + { + const TID tid = std::get<2>(partial_edges[i]); + const auto& t = triangles[tid]; + faces.push_back(tid); + + const std::size_t k = (t[0] == a || t[0] == b) ? ((t[1] == a || t[1] == b) ? 2 : 1) : 0; + CGAL_assertion((t[(k+1)%3] == a && t[(k+2)%3] == b) || (t[(k+1)%3] == b && t[(k+2)%3] == a)); + adjacent_edges[tid][k] = idx; + ++i; + } while (i < partial_edges.size() && std::get<0>(partial_edges[i]) == a && + std::get<1>(partial_edges[i]) == b); + ++idx; + } + CGAL_assertion(idx == offsets[bi+1]); + } + }); + } + else +#endif + { + adjacent_faces.reserve(partial_edges.size()); + for (std::size_t i = 0; i < partial_edges.size(); ) + { + const PID a = std::get<0>(partial_edges[i]); + const PID b = std::get<1>(partial_edges[i]); + const std::size_t idx = adjacent_faces.size(); + adjacent_faces.emplace_back(); + CGAL_assertion(adjacent_faces.back().empty()); + do + { + const TID tid = std::get<2>(partial_edges[i]); + const auto& t = triangles[tid]; + adjacent_faces.back().push_back(tid); + ++i; + + std::size_t k = (t[0] == a || t[0] == b) ? ((t[1] == a || t[1] == b) ? 2 : 1) : 0; + CGAL_assertion((t[(k+1)%3] == a && t[(k+2)%3] == b) || (t[(k+1)%3] == b && t[(k+2)%3] == a)); + adjacent_edges[tid][k] = idx; + } while (i < partial_edges.size() && std::get<0>(partial_edges[i]) == a && std::get<1>(partial_edges[i]) == b); + CGAL_assertion(!adjacent_faces.back().empty()); + } + } + +#ifdef CGAL_BO3_SELECT_TIMERS + std::cout << t.time() << " sec. to build adjacencies\n"; +#endif + + // extract connected_components bounded by non-manifold edges + std::vector& ccids = data_out.ccids; + std::vector< EID > nm_edges; + std::vector one_triangle_per_cc; +#ifdef CGAL_LINKED_WITH_TBB + if constexpr(std::is_same_v) + { + one_triangle_per_cc = boolops_3::parallel_connected_components(split_at, new_triangles_at, adjacent_edges, adjacent_faces, ccids, nm_edges); + } + else +#endif + { + one_triangle_per_cc = boolops_3::connected_components(adjacent_edges, adjacent_faces, ccids, nm_edges); + } + std::size_t nb_cc = one_triangle_per_cc.size(); + +#ifdef CGAL_BO3_SELECT_TIMERS + std::cout << t.time() << " sec. to compute connected components\n"; +#endif + + // Erase duplicates + std::sort(nm_edges.begin(), nm_edges.end()); + nm_edges.erase(std::unique(nm_edges.begin(), nm_edges.end()), nm_edges.end()); + +#ifdef CGAL_BO_AUTOREF_DEBUG + std::cout << "nb_cc " << nb_cc << "\n"; + std::vector< std::vector< boost::container::small_vector > > triangles_per_cc(nb_cc); + for (std::size_t tid=0; tid cc_points = points; + PMP::remove_isolated_points_in_polygon_soup(cc_points, triangles_per_cc[i]); + CGAL::IO::write_OFF("cc_"+std::to_string(i)+".off", cc_points, triangles_per_cc[i], CGAL::parameters::stream_precision(17)); + } +#endif + + std::size_t nb_handled = 0; + std::vector cc_handled(nb_cc, false); + std::vector& coplanar_patches=data_out.coplanar_patches; + std::vector& coplanar_patches_for_union_and_intersection=data_out.coplanar_patches_for_union_and_intersection; + std::vector< std::vector >& cc_is_inside_per_mesh=data_out.cc_is_inside_per_mesh; + + for(std::size_t cc=0; cc > cc_is_inside_per_mesh_init(nb_meshes, std::vector(nb_cc, false)); + int case_id=-1; +#endif + cc_is_inside_per_mesh.assign(nb_meshes, std::vector(nb_cc, false)); // TODO: use an optional + + for (auto e : nm_edges) + { + if (nb_handled==nb_cc) break; + + bool all_handled=true; + auto &tids = adjacent_faces[e]; + for (TID tid : tids) + { + if (!cc_handled[ccids[tid]]) + { + all_handled=false; + break; + } + } + + if (all_handled) continue; + + bool all_from_same_mesh = true; + + for (TID tid : tids) + { + if (input_ids[tid]!=input_ids[tids.front()]) + { + all_from_same_mesh=false; + break; + } + } + if (all_from_same_mesh) continue; +#ifdef CGAL_BO_AUTOREF_DEBUG + std::cout << "case #" << ++case_id << " looking at edge " << pid << " " << e << "\n"; + std::cout << " incident CCs:"; + for (TID tid : tids) std::cout << " " << ccids[tid] << "("<< input_ids[tid] <<")"; + std::cout << "\n"; +#endif + for (TID tid : tids) + { + if (!cc_handled[ccids[tid]]) ++nb_handled; + cc_handled[ccids[tid]] = true; + } + + std::vector> to_sort; + to_sort.reserve(tids.size()); + + for (TID tid : tids){ + for (int i=0; i<3; ++i) + if( adjacent_edges[tid][i] == e ) + to_sort.emplace_back(i, tid); + } + + // we arbitrarily take the last triangle + std::pair tref=to_sort.front(); + std::size_t i0=triangles[tref.second][(tref.first+1)%3], + i1=triangles[tref.second][(tref.first+2)%3]; + + const Exact_point& pref=points[triangles[tref.second][tref.first]]; + const Exact_point& p0=points[i0]; + const Exact_point& p1=points[i1]; + + //Considering the plane with normal vector [o_prime,o] and containing o. + //We define the counterclockwise order around o when looking from + //the side of the plane containing o_prime. + //We consider the portion of the plane defined by rotating a ray starting at o + //from the planar projection of p1 to the planar projection of p2 in + //counterclockwise order. + //The predicates indicates whether the planar projection of point q lies in this + //portion of the plane. + //Preconditions: + // o_prime,o,p1 are not collinear + // o_prime,o,p2 are not collinear + // o_prime,o,q are not collinear + // o_prime,o,p1,q are not coplanar or coplanar_orientation(o,o_prime,p1,q)==NEGATIVE + // o_prime,o,p2,q are not coplanar or coplanar_orientation(o,o_prime,p2,q)==NEGATIVE + // template + // bool sorted_around_edge( + // const typename Kernel::Point_3& o_prime, const typename Kernel::Point_3& o, + // const typename Kernel::Point_3& p1, const typename Kernel::Point_3& p2, + // const typename Kernel::Point_3& q) + auto less = [&tref, &pref, &points, &triangles, &p0, &p1](const std::pair& e1, + const std::pair& e2) + { + using EPECK = CGAL::Exact_predicates_exact_constructions_kernel; + // handling coplanar + //TODO: we will make several time the comparision with the same point, not sure if it is coslty. Maybe group them? + if (triangles[e1.second][e1.first]==triangles[e2.second][e2.first]) + return e1.second < e2.second; + if (triangles[e1.second][e1.first]==triangles[tref.second][tref.first]) + return true; + if (triangles[e2.second][e2.first]==triangles[tref.second][tref.first]) + return false; + + return pred::sorted_around_edge(p0, p1, pref, + points[triangles[e2.second][e2.first]], + points[triangles[e1.second][e1.first]]); + }; + +#ifdef CGAL_BO_AUTOREF_DEBUG + std::set third_point; + for (const std::pair& p : to_sort) + third_point.insert(triangles[p.second][p.first]); + + if (third_point.size()!=to_sort.size()+1) + std::cerr << " presence of coplanar triangles\n"; +#endif + + std::sort(std::next(to_sort.begin()), to_sort.end(), less); + boost::dynamic_bitset<> is_inside_mesh(nb_meshes, 0); //indicate the status while turning around the edge + + std::vector init(nb_meshes, false); + init[input_ids[tref.second]]=true; + is_inside_mesh[input_ids[tref.second]]=true; + for (std::size_t pos=1; pos& p = to_sort[pos]; + if (init[input_ids[p.second]]) continue; + init[input_ids[p.second]]=true; + CGAL_assertion(pos!=to_sort.size()-1); // not possible since there must be at least 2 triangles per mesh + + //check if collinear with ref + if (triangles[p.second][p.first]==triangles[tref.second][tref.first]) + { + // incompatible orientation -> flip the bit + if (triangles[p.second][(p.first+1)%3]==triangles[tref.second][(tref.first+1)%3]) + is_inside_mesh.set(input_ids[p.second]); + } + else + // incompatible orientation -> flip the bit + if (triangles[p.second][(p.first+1)%3]==triangles[tref.second][(tref.first+1)%3]) + is_inside_mesh.set(input_ids[p.second]); + } +#ifdef CGAL_BO_AUTOREF_DEBUG + auto print_t_info = [&](std::size_t tid) + { + std::cout << " " << tid << " from mesh " << input_ids[tid] << " part of CC " << ccids[tid] << " " + << points[triangles[tid][0]] << " " << points[triangles[tid][1]] << " " << points[triangles[tid][2]] << "\n"; + }; + + int nb=-1; + for (const std::pair& p : to_sort) + { + std::ofstream debug("case_"+ std::to_string(case_id)+"_"+std::to_string(++nb)+".polylines.txt"); + debug << 4 << " " << points[triangles[p.second][0]] << " "<< " " << points[triangles[p.second][1]] << " " << points[triangles[p.second][2]] << " " << points[triangles[p.second][0]] << "\n"; + debug.close(); + print_t_info(p.second); + } + std::cout << " is_inside_mesh = "; + for (unsigned int i=0; i& p = to_sort[pos]; + ++pos; + std::vector mesh_ids; + mesh_ids.push_back(input_ids[p.second]); + + while (pos1) + { + coplanar_patches[ccids[p.second]]=true; + std::sort(mesh_ids.begin(), mesh_ids.end()); + + //check orientation -> mark only if all the orientations are identical + bool same_orientation=true; + for (std::size_t k=start_pos+1;k > triangles_of_intersect_per_cc(nb_cc); + std::vector< std::pair > points_and_cc; + points_and_cc.reserve(nb_cc-nb_handled); + for (std::size_t cc=0; cc bb_per_cc(nb_cc); + for (std::size_t tid=0; tid pt.x() || t_bbox.xmax() < pt.x() || + t_bbox.ymin() > pt.y() || t_bbox.ymax() < pt.y() || + t_bbox.zmax() < pt.z() ) continue; + triangles_of_intersect_per_cc[cc].push_back(tid); + } + } + + //bbox of 2 CCs that are not nested --> cc are disjoint (because points are shared) + //TODO: this part works as is only for 2 meshes + for (std::size_t cc=0; cc& subtri=triangles_of_intersect_per_cc[cc]; + const Exact_point& pt = points[triangles[one_triangle_per_cc[cc]][0]]; + Exact_predicates_exact_constructions_kernel::Ray_3 ray(pt, Exact_predicates_exact_constructions_kernel::Vector_3(0.,0.,1.)); + using T3 = Exact_predicates_exact_constructions_kernel::Triangle_3; + std::vector> to_sort; + for (std::size_t tid : subtri) + { + T3 tr(points[triangles[tid][0]],points[triangles[tid][1]],points[triangles[tid][2]]); + if (do_intersect(ray, tr)) + to_sort.emplace_back(tid, tr); + } + + std::sort(to_sort.begin(), to_sort.end(), + [&pt](const std::pair& p1, const std::pair& p2) + { + return compare_distance(pt, p1.second, p2.second) == CGAL::SMALLER; + } + ); + + bool is_outside = true; + for (const std::pair& p : to_sort) + { + CGAL::Orientation ori = orientation(p.second[0], p.second[1], p.second[2], pt); + if (ori == CGAL::COPLANAR) continue; + is_outside = (ori == POSITIVE); + + break; + } + + cc_handled[cc]=true; + if (!is_outside) + // TODO: clearly working only for 2 meshes + cc_is_inside_per_mesh[ (input_ids[one_triangle_per_cc[cc]]+1)%2 ][cc]=true; + } + } + } + + + data_out.input_ids.swap(input_ids); +#ifdef CGAL_BO3_SELECT_TIMERS + std::cout << "select done in " << t.time() << "\n"; +#endif + return data_out; +} + +struct Autoref_visitor +{ + Autoref_visitor(std::vector& input_ids, + std::size_t nbt1, std::size_t nbt2) + : input_ids(input_ids) + , nbt1(nbt1) + , nbt2(nbt2) + {} + + inline void number_of_output_triangles(std::size_t nbt) + { + input_ids.resize(nbt, -1); + } + inline void verbatim_triangle_copy(std::size_t tgt_id, std::size_t src_id) + { + input_ids[tgt_id] = src_id < nbt1? 0 : 1; + } + inline void new_subtriangle(std::size_t tgt_id, std::size_t src_id) + { + input_ids[tgt_id] = src_id < nbt1? 0 : 1; + } + inline void delete_triangle(std::size_t /*src_id*/) {} + + std::vector& input_ids; + std::size_t nbt1; + std::size_t nbt2; +}; + +template +Selection_data +prepare_data(const PointRange& points_1, const TriangleRange& triangles_1, + const PointRange& points_2, const TriangleRange& triangles_2, + PointRange& out_points, TriangleRange& out_triangles) +{ + namespace PMP = CGAL::Polygon_mesh_processing; + +#ifdef CGAL_BO3_TIMERS + Real_timer timer; + timer.start(); +#endif + out_points.reserve(points_1.size()+points_2.size()); + out_points.insert(out_points.end(), points_1.begin(), points_1.end()); + out_points.insert(out_points.end(), points_2.begin(), points_2.end()); + out_triangles.reserve(triangles_1.size()+triangles_2.size()); + out_triangles.insert(out_triangles.end(), triangles_1.begin(), triangles_1.end()); + std::size_t nbt_1=triangles_1.size(), nbt_2=triangles_2.size(); + + for (auto t : triangles_2) + { + for (int i=0;i<3;++i) + t[i]+=points_1.size(); + out_triangles.push_back(t); + } + + std::vector input_ids(nbt_1+nbt_2, 0); + for (std::size_t i=nbt_1; i(out_points, out_triangles, input_ids, nbt_1, nbt_1+nbt_2, 2); +#ifdef CGAL_BO3_TIMERS + timer.stop(); + std::cout << "select done in " << timer.time() << "\n"; +#endif + return res; +} + +template +void +update_output(PointRange& out_points, TriangleRange& out_triangles, std::size_t nbt, const std::vector& to_keep, const std::vector& to_flip=std::vector()) +{ + if (nbt==out_triangles.size() && to_flip.empty()) return; + + TriangleRange triangles; + triangles.reserve(nbt); + + PointRange points; + points.reserve(nbt<3?3:((nbt-2)/2+1)); + std::vector point_map(out_points.size(), out_points.size()); + + auto point_id = [&](std::size_t i) + { + if (point_map[i]==out_points.size()) + { + point_map[i]=points.size(); + points.push_back(out_points[i]); + } + return point_map[i]; + }; + + for (std::size_t ti=0; ti +void compute_self_union(PointRange& points, TriangleRange& triangles, + const NamedParameters& np = parameters::default_values()) +{ + CGAL_USE(np); + + typedef typename internal_np::Lookup_named_param_def < + internal_np::concurrency_tag_t, + NamedParameters, + Sequential_tag + > ::type Concurrency_tag; + + namespace PMP = CGAL::Polygon_mesh_processing; + using namespace boolops_3; + +#ifdef CGAL_BO3_TIMERS + Real_timer timer; + timer.start(); +#endif + // TODO: convert points to EPECK if not already the case + using Exact_point = typename std::iterator_traits::value_type; + PMP::autorefine_triangle_soup(points, triangles, CGAL::parameters::concurrency_tag(Concurrency_tag())); +#ifdef CGAL_BO3_TIMERS + timer.stop(); + std::cout << "autorefine done in " << timer.time() << "\n"; + timer.reset(); +#endif +#ifdef CGAL_BO_AUTOREF_DEBUG + CGAL::IO::write_OFF("autorefined.off", out_points, out_triangles, CGAL::parameters::stream_precision(17)); +#endif + + //TODO: use Concurrency_tagf + namespace PMP = CGAL::Polygon_mesh_processing; + namespace pred = PMP::Corefinement; + //~ Selection_data data_out; + + // fill edge map (TODO: move to a function) + using PID = std::size_t; + using EID = std::size_t; + using TID = std::size_t; + using Edge_record = std::tuple; + + std::vector partial_edges(3 * triangles.size()); + for (TID tid = 0; tid < triangles.size(); ++tid) + { + auto add_edge = [&](std::size_t idx, PID a, PID b) + { + if (b < a) std::swap(a, b); + partial_edges[idx] = Edge_record{a, b, tid}; + }; + add_edge(3 * tid, triangles[tid][0], triangles[tid][1]); + add_edge(3 * tid + 1, triangles[tid][1], triangles[tid][2]); + add_edge(3 * tid + 2, triangles[tid][2], triangles[tid][0]); + } + std::sort(partial_edges.begin(), partial_edges.end()); + + std::vector> adjacent_edges(triangles.size()); + std::vector> adjacent_faces; + adjacent_faces.reserve(partial_edges.size()); + for (std::size_t i = 0; i < partial_edges.size(); ) + { + const PID a = std::get<0>(partial_edges[i]); + const PID b = std::get<1>(partial_edges[i]); + const EID eid = adjacent_faces.size(); + adjacent_faces.emplace_back(); + do + { + const TID tid = std::get<2>(partial_edges[i++]); + adjacent_faces.back().push_back(tid); + for (std::size_t j = 0; j < 3; ++j) + if ((triangles[tid][j] == a && triangles[tid][(j+1)%3] == b) || + (triangles[tid][j] == b && triangles[tid][(j+1)%3] == a)) + { + adjacent_edges[tid][(j+2)%3] = eid; + break; + } + } while (i < partial_edges.size() && std::get<0>(partial_edges[i]) == a && std::get<1>(partial_edges[i]) == b); + } + + // extract connected_components bounded by non-manifold edges + std::vector ccids; + std::vector nm_edges; + std::vector one_triangle_per_cc = + boolops_3::connected_components(adjacent_edges, adjacent_faces, ccids, nm_edges); + std::size_t nb_cc = one_triangle_per_cc.size(); + + std::sort(nm_edges.begin(), nm_edges.end()); + nm_edges.erase(std::unique(nm_edges.begin(), nm_edges.end()), nm_edges.end()); + + // Get the vertices associate to an edge + auto edge_vertices = [&](EID eid) { + std::size_t tid = adjacent_faces[eid].front(); + std::size_t k = adjacent_edges[tid][0] == eid ? 0 :adjacent_edges[tid][1] == eid ? 1 : 2; + CGAL_assertion(adjacent_edges[tid][k] == eid); + return std::make_pair(triangles[tid][(k+1)%3], triangles[tid][(k+2)%3]); + }; + + + using UF_CC = CGAL::Union_find; + using UF_handle = UF_CC::handle; + + UF_CC uf_cc; + std::vector uf_handles(nb_cc); + for (std::size_t i=0; i edges + // TODO we need split graph into polylines I guess + std::vector > cc_to_edges(nb_cc); + std::vector pts_on_edge(points.size(), false); // TODO: fill it in autoref visitor? + + for (EID eid : nm_edges) + { + std::set ccs; + const auto [pid, other_pid] = edge_vertices(eid); + pts_on_edge[pid]=true; + pts_on_edge[other_pid]=true; + for (TID tid : adjacent_faces[eid]) + ccs.insert(ccids[tid]); + for (std::size_t ccid : ccs) + cc_to_edges[ccid].push_back(eid); + } + + std::size_t min_pt_id=0; + for (std::size_t i=1; ipoints[i].z()) min_pt_id=i; + } + else + if (points[min_pt_id].y()>points[i].y()) min_pt_id=i; + } + else + if (points[min_pt_id].x()>points[i].x()) min_pt_id=i; + } + + //TODO: temporary condition assuming extreme point is manifold, need another selection method otherwise + CGAL_precondition(!pts_on_edge[min_pt_id]); + + std::size_t selected_tid=triangles.size(); + for (TID tid=0; tid ccs_handled(nb_cc, false); + std::vector cc_stack; + cc_stack.push_back(ccids[selected_tid]); + + // main loop to classify cc + while(!cc_stack.empty()) + { + std::size_t ccid = cc_stack.back(); + cc_stack.pop_back(); + + if (ccs_handled[ccid]) continue; + ccs_handled[ccid]=true; + + // iterate on all non-manifold edges of the current cc + for (EID eid : cc_to_edges[ccid]) + { + const auto [pid, other_pid] = edge_vertices(eid); +#ifdef CGAL_BO_AUTOREF_DEBUG + std::cout << "case #" << ++case_id << " looking at edge " << pid << " " << other_pid << "\n"; + std::cout << " incident CCs:"; + for (TID tid : adjacent_faces[eid]) std::cout << " " << ccids[tid] << "("<< input_ids[tid] <<")"; + std::cout << "\n"; +#endif + + std::vector> to_sort; + to_sort.reserve(adjacent_faces[eid].size()); + + std::size_t i0=pid, i1=other_pid; + for (TID tid : adjacent_faces[eid]) + { + for (int i=0; i<3; ++i) + { + if (triangles[tid][i]!=i0 && triangles[tid][i]!=i1) + to_sort.emplace_back(i, tid); + } + if (ccids[tid]==ccid && to_sort.size()>1) + std::swap(to_sort.front(), to_sort.back()); + } + + std::pair tref=to_sort.front(); // tid is from ccid + if (triangles[tref.second][(tref.first+1)%3]!=i0) + std::swap(i0,i1); + + const Exact_point& pref=points[triangles[tref.second][tref.first]]; + const Exact_point& p0=points[i0]; + const Exact_point& p1=points[i1]; + + //Considering the plane with normal vector [o_prime,o] and containing o. + //We define the counterclockwise order around o when looking from + //the side of the plane containing o_prime. + //We consider the portion of the plane defined by rotating a ray starting at o + //from the planar projection of p1 to the planar projection of p2 in + //counterclockwise order. + //The predicates indicates whether the planar projection of point q lies in this + //portion of the plane. + //Preconditions: + // o_prime,o,p1 are not collinear + // o_prime,o,p2 are not collinear + // o_prime,o,q are not collinear + // o_prime,o,p1,q are not coplanar or coplanar_orientation(o,o_prime,p1,q)==NEGATIVE + // o_prime,o,p2,q are not coplanar or coplanar_orientation(o,o_prime,p2,q)==NEGATIVE + // template + // bool sorted_around_edge( + // const typename Kernel::Point_3& o_prime, const typename Kernel::Point_3& o, + // const typename Kernel::Point_3& p1, const typename Kernel::Point_3& p2, + // const typename Kernel::Point_3& q) + auto less = [&tref, &pref, &points, &triangles, &p0, &p1](const std::pair& e1, + const std::pair& e2) + { + using EPECK = CGAL::Exact_predicates_exact_constructions_kernel; + // handling coplanar + //TODO: we will make several time the comparision with the same point, not sure if it is coslty. Maybe group them? + if (triangles[e1.second][e1.first]==triangles[e2.second][e2.first]) + return e1.second < e2.second; + if (triangles[e1.second][e1.first]==triangles[tref.second][tref.first]) + return true; + if (triangles[e2.second][e2.first]==triangles[tref.second][tref.first]) + return false; + + return pred::sorted_around_edge(p0, p1, pref, + points[triangles[e2.second][e2.first]], + points[triangles[e1.second][e1.first]]); + }; + + #ifdef CGAL_BO_AUTOREF_DEBUG + std::set third_point; + for (const std::pair& p : to_sort) + third_point.insert(triangles[p.second][p.first]); + + if (third_point.size()!=to_sort.size()+1) + std::cerr << " presence of coplanar triangles\n"; + #endif + + std::sort(std::next(to_sort.begin()), to_sort.end(), less); + + // TODO: we need to make sure orentations are compatible + + std::size_t n_ccid = ccids[to_sort[1].second]; + uf_handles[n_ccid] = uf_cc.find(uf_handles[n_ccid]); + uf_handles[ccid] = uf_cc.find(uf_handles[ccid]); + uf_cc.unify_sets(uf_handles[n_ccid], uf_handles[ccid]); + if (!ccs_handled[n_ccid]) + cc_stack.push_back(n_ccid); + } + } + + // mark ccids to keep + std::size_t ref_cc = ccids[selected_tid]; + UF_handle ref_handle = uf_cc.find(uf_handles[ref_cc]); + std::vector selected_ccs(nb_cc); + for (std::size_t ccid=0; ccid to_keep(triangles.size(), false); + std::size_t nb_out_triangles=0; + for (std::size_t tid=0; tid(points, triangles, nb_out_triangles, to_keep); +} +#endif + +/** + * \ingroup PMP_boolop_soup_grp + * + * computes the union of two volumes, each bounded by a triangle soup. + * + * \pre \link CGAL::Polygon_mesh_processing::does_triangle_soup_self_intersect() `!CGAL::Polygon_mesh_processing::does_triangle_soup_self_intersect(points_1, triangles_1)` \endlink + * \pre \link CGAL::Polygon_mesh_processing::does_triangle_soup_self_intersect() `!CGAL::Polygon_mesh_processing::does_triangle_soup_self_intersect(points_2, triangles_2)` \endlink + * + * @tparam PointRange a model of the concepts `RandomAccessContainer` whose value type is `GeomTraits::Point_3` as value type, `GeomTraits` being the type of the named parameter `geom_traits`. + * @tparam TriangleRange a model of the concept `RandomAccessContainer` whose `value_type` is a model of the concept `RandomAccessContainer` whose `value_type` is `std::size_t`. + * + * @param points_1 points of the first triangle soup. + * @param triangles_1 triangles of the first triangle soup, each element in the vector describes a triangle using the indices of the points in `points_1`. + * + * @param points_2 points of the second triangle soup. + * @param triangles_2 triangles of the second triangle soup, each element in the vector describes a triangle using the indices of the points in `points_2`. + * + * @param out_points points of the output triangle soup bounding the volume result of the Boolean operation. + * @param out_triangles triangles of the ouput triangle soup, each element in the vector describes a triangle using the indices of the points in `out_points`. + * + * @param np optional sequence of named parameters among the ones listed below + * + * \cgalNamedParamsBegin + * \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::Sequential_tag`} + * \cgalParamNEnd + * \cgalParamNBegin{geom_traits} + * \cgalParamDescription{an instance of a geometric traits class} + * \cgalParamType{a class model of `Kernel`} + * \cgalParamDefault{a \cgal Kernel deduced from the point type, using `CGAL::Kernel_traits`} + * \cgalParamExtra{The geometric traits class must be compatible with the vertex point type.} + * \cgalParamNEnd + * \cgalNamedParamsEnd + */ +template +void compute_union(const PointRange& points_1, const TriangleRange& triangles_1, + const PointRange& points_2, const TriangleRange& triangles_2, + PointRange& out_points, TriangleRange& out_triangles, + const NamedParameters& np = parameters::default_values()) +{ + CGAL_USE(np); + + typedef typename internal_np::Lookup_named_param_def < + internal_np::concurrency_tag_t, + NamedParameters, + Sequential_tag + > ::type Concurrency_tag; + +#ifdef CGAL_BO3_TIMERS + Real_timer t; + t.start(); +#endif + + using namespace boolops_3; + Selection_data data_out=prepare_data(points_1, triangles_1, points_2, triangles_2, out_points, out_triangles); + const unsigned int nb_meshes=2; + +#ifdef CGAL_BO3_TIMERS + std::cout << "Prepare_data + autoref: " << t.time() << "\n"; +#endif + std::vector to_keep(out_triangles.size(), false); + std::size_t nb_out_triangles=0; + + for (std::size_t ti=0; ti(out_points, out_triangles, nb_out_triangles, to_keep); +#ifdef CGAL_BO3_TIMERS + std::cout << "BO3 Building output: " << t.time() << "\n"; +#endif +} + +/** + * \ingroup PMP_boolop_soup_grp + * computes the intersection of two volumes, each bounded by a triangle soup. + * \copydetails CGAL::Polygon_mesh_processing::compute_union() + */ +template +void compute_intersection(const PointRange& points_1, const TriangleRange& triangles_1, + const PointRange& points_2, const TriangleRange& triangles_2, + PointRange& out_points, TriangleRange& out_triangles, + const NamedParameters& np = parameters::default_values()) +{ + CGAL_USE(np); + + typedef typename internal_np::Lookup_named_param_def < + internal_np::concurrency_tag_t, + NamedParameters, + Sequential_tag + > ::type Concurrency_tag; + + using namespace boolops_3; + Selection_data data_out=prepare_data(points_1, triangles_1, points_2, triangles_2, out_points, out_triangles); + const unsigned int nb_meshes=2; + + std::vector to_keep(out_triangles.size(), false); + std::size_t nb_out_triangles=0; + + for (std::size_t ti=0; ti(out_points, out_triangles, nb_out_triangles, to_keep); +} + +/** + * \ingroup PMP_boolop_soup_grp + * computes the difference of two volumes, each bounded by a triangle soup. + * \copydetails CGAL::Polygon_mesh_processing::compute_union() + */ +template +void compute_difference(const PointRange& points_1, const TriangleRange& triangles_1, + const PointRange& points_2, const TriangleRange& triangles_2, + PointRange& out_points, TriangleRange& out_triangles, + const NamedParameters& np = parameters::default_values()) +{ + CGAL_USE(np); + + typedef typename internal_np::Lookup_named_param_def < + internal_np::concurrency_tag_t, + NamedParameters, + Sequential_tag + > ::type Concurrency_tag; + + using namespace boolops_3; + Selection_data data_out=prepare_data(points_1, triangles_1, points_2, triangles_2, out_points, out_triangles); + // const unsigned int nb_meshes=2; + + std::vector to_keep(out_triangles.size(), false); + std::vector to_flip(out_triangles.size(), false); + std::size_t nb_out_triangles=0; + + for (std::size_t ti=0; ti(out_points, out_triangles, nb_out_triangles, to_keep, to_flip); +} + +} } // CGAL::Polygon_mesh_processing namespace + +#endif // CGAL_POLYGON_MESH_PROCESSING_TRIANGLE_SOUP_BOOLEAN_OPERATIONS_3_H diff --git a/PMP_Boolean_operations/test/PMP_Boolean_operations/CMakeLists.txt b/PMP_Boolean_operations/test/PMP_Boolean_operations/CMakeLists.txt index 82dd242e9a8c..a52f5b7347ae 100644 --- a/PMP_Boolean_operations/test/PMP_Boolean_operations/CMakeLists.txt +++ b/PMP_Boolean_operations/test/PMP_Boolean_operations/CMakeLists.txt @@ -8,6 +8,7 @@ project(PMP_Boolean_operations_Tests) find_package(CGAL REQUIRED) create_single_source_cgal_program("autorefinement_sm.cpp") +create_single_source_cgal_program("self_union.cpp") create_single_source_cgal_program("corefine_non_manifold.cpp" ) create_single_source_cgal_program("polygon_mesh_slicer_test.cpp") create_single_source_cgal_program("surface_intersection_sm_poly.cpp") @@ -16,6 +17,7 @@ create_single_source_cgal_program("test_corefine.cpp") create_single_source_cgal_program("test_coref_epic_points_identity.cpp") create_single_source_cgal_program("test_corefinement_and_constraints.cpp") create_single_source_cgal_program("test_corefinement_and_constraints_nearest.cpp") +create_single_source_cgal_program("test_corefinement_vs_autorefine.cpp") create_single_source_cgal_program("test_corefinement_bool_op.cpp") create_single_source_cgal_program("test_corefinement_cavities.cpp") create_single_source_cgal_program("test_corefinement_nm_bo.cpp") @@ -24,12 +26,15 @@ create_single_source_cgal_program("test_isolevel_refinement.cpp") create_single_source_cgal_program("test_pmp_clip.cpp") create_single_source_cgal_program("test_snap_rounding.cpp") create_single_source_cgal_program("test_mesh_kernel.cpp") +create_single_source_cgal_program("test_boolean_on_triangle_soup.cpp") find_package(TBB QUIET) include(CGAL_TBB_support) if(TARGET CGAL::TBB_support) target_link_libraries(test_autorefinement PRIVATE CGAL::TBB_support) + target_link_libraries(self_union PRIVATE CGAL::TBB_support) target_link_libraries(test_snap_rounding PRIVATE CGAL::TBB_support) + target_link_libraries(test_boolean_on_triangle_soup PRIVATE CGAL::TBB_support) else() message(STATUS "NOTICE: Intel TBB was not found. Tests will use sequential code.") endif() diff --git a/PMP_Boolean_operations/test/PMP_Boolean_operations/data-nm_boolop/cube3_minus_prism.off b/PMP_Boolean_operations/test/PMP_Boolean_operations/data-nm_boolop/cube3_minus_prism.off new file mode 100644 index 000000000000..30228694bd49 --- /dev/null +++ b/PMP_Boolean_operations/test/PMP_Boolean_operations/data-nm_boolop/cube3_minus_prism.off @@ -0,0 +1,48 @@ +OFF +16 28 0 + +-1 0 -1 +1 0 -1 +0 2 -1 +1 0 1 +0 2 1 +-1 0 1 +-2 2 -1 +2 2 -1 +-2 -2 -1 +-2 -2 1 +-2 2 1 +2 -2 1 +2 -2 -1 +2 2 1 +0 2 -1 +0 2 1 +3 9 3 5 +3 8 6 0 +3 8 9 10 +3 12 7 13 +3 8 12 11 +3 0 6 14 +3 2 15 13 +3 4 14 6 +3 11 3 9 +3 4 5 0 +3 2 1 3 +3 3 1 0 +3 12 0 1 +3 8 0 12 +3 12 1 7 +3 1 2 7 +3 0 5 3 +3 10 6 8 +3 5 4 10 +3 3 15 2 +3 13 11 12 +3 0 14 4 +3 11 9 8 +3 3 13 15 +3 11 13 3 +3 13 7 2 +3 6 10 4 +3 9 5 10 + diff --git a/PMP_Boolean_operations/test/PMP_Boolean_operations/data-nm_boolop/cube4.off b/PMP_Boolean_operations/test/PMP_Boolean_operations/data-nm_boolop/cube4.off new file mode 100644 index 000000000000..023b160982b3 --- /dev/null +++ b/PMP_Boolean_operations/test/PMP_Boolean_operations/data-nm_boolop/cube4.off @@ -0,0 +1,22 @@ +OFF +8 12 0 +0 1 -1 +0 2 -1 +1 2 -1 +1 1 -1 +0 1 1 +0 2 1 +1 2 1 +1 1 1 +3 0 1 3 +3 3 1 2 +3 0 4 1 +3 1 4 5 +3 3 2 7 +3 7 2 6 +3 4 0 3 +3 7 4 3 +3 6 4 7 +3 6 5 4 +3 1 5 6 +3 2 1 6 diff --git a/PMP_Boolean_operations/test/PMP_Boolean_operations/data-nm_boolop/cube_0.off b/PMP_Boolean_operations/test/PMP_Boolean_operations/data-nm_boolop/cube_0.off new file mode 100644 index 000000000000..3d7bf7101998 --- /dev/null +++ b/PMP_Boolean_operations/test/PMP_Boolean_operations/data-nm_boolop/cube_0.off @@ -0,0 +1,24 @@ +OFF +8 12 0 + +-0.5 -0.5 0.5 +0.5 -0.5 0.5 +0.5 -0.5 -0.5 +-0.5 -0.5 -0.5 +-0.5 0.5 -0.5 +-0.5 0.5 0.5 +0.5 0.5 0.5 +0.5 0.5 -0.5 +3 4 5 6 +3 1 0 3 +3 6 1 2 +3 0 1 6 +3 3 0 5 +3 7 2 3 +3 6 7 4 +3 3 2 1 +3 2 7 6 +3 6 5 0 +3 5 4 3 +3 3 4 7 + diff --git a/PMP_Boolean_operations/test/PMP_Boolean_operations/data-nm_boolop/cube_1.off b/PMP_Boolean_operations/test/PMP_Boolean_operations/data-nm_boolop/cube_1.off new file mode 100644 index 000000000000..534270d2bc86 --- /dev/null +++ b/PMP_Boolean_operations/test/PMP_Boolean_operations/data-nm_boolop/cube_1.off @@ -0,0 +1,24 @@ +OFF +8 12 0 + +0.054291486740112305 0.26874995231628418 -0.34933948516845703 +1.0542914867401123 0.26874995231628418 -0.34933948516845703 +1.0542914867401123 0.26874995231628418 -1.349339485168457 +0.054291486740112305 0.26874995231628418 -1.349339485168457 +0.054291486740112305 1.2687499523162842 -1.349339485168457 +0.054291486740112305 1.2687499523162842 -0.34933948516845703 +1.0542914867401123 1.2687499523162842 -0.34933948516845703 +1.0542914867401123 1.2687499523162842 -1.349339485168457 +3 4 5 6 +3 1 0 3 +3 6 1 2 +3 0 1 6 +3 3 0 5 +3 7 2 3 +3 6 7 4 +3 3 2 1 +3 2 7 6 +3 6 5 0 +3 5 4 3 +3 3 4 7 + diff --git a/PMP_Boolean_operations/test/PMP_Boolean_operations/data-nm_boolop/cube_2.off b/PMP_Boolean_operations/test/PMP_Boolean_operations/data-nm_boolop/cube_2.off new file mode 100644 index 000000000000..bc4ce81514fc --- /dev/null +++ b/PMP_Boolean_operations/test/PMP_Boolean_operations/data-nm_boolop/cube_2.off @@ -0,0 +1,24 @@ +OFF +8 12 0 + +-0.109773188829422 -0.065066307783126831 1.1743011474609375 +0.89022684097290039 -0.065066307783126831 1.1743011474609375 +0.89022684097290039 -0.065066307783126831 0.1743011474609375 +-0.109773188829422 -0.065066307783126831 0.1743011474609375 +-0.109773188829422 0.93493366241455078 0.1743011474609375 +-0.109773188829422 0.93493366241455078 1.1743011474609375 +0.89022684097290039 0.93493366241455078 1.1743011474609375 +0.89022684097290039 0.93493366241455078 0.1743011474609375 +3 4 5 6 +3 1 0 3 +3 6 1 2 +3 0 1 6 +3 3 0 5 +3 7 2 3 +3 6 7 4 +3 3 2 1 +3 2 7 6 +3 6 5 0 +3 5 4 3 +3 3 4 7 + diff --git a/PMP_Boolean_operations/test/PMP_Boolean_operations/data-nm_boolop/nested_cubes_1.off b/PMP_Boolean_operations/test/PMP_Boolean_operations/data-nm_boolop/nested_cubes_1.off new file mode 100644 index 000000000000..824fad12f3b2 --- /dev/null +++ b/PMP_Boolean_operations/test/PMP_Boolean_operations/data-nm_boolop/nested_cubes_1.off @@ -0,0 +1,24 @@ +OFF +8 12 0 + +-2 -2 2 +2 -2 2 +2 -2 -2 +-2 -2 -2 +-2 2 -2 +-2 2 2 +2 2 2 +2 2 -2 +3 4 5 6 +3 1 0 3 +3 6 1 2 +3 0 1 6 +3 3 0 5 +3 7 2 3 +3 6 7 4 +3 3 2 1 +3 2 7 6 +3 6 5 0 +3 5 4 3 +3 3 4 7 + diff --git a/PMP_Boolean_operations/test/PMP_Boolean_operations/data-nm_boolop/nested_cubes_2.off b/PMP_Boolean_operations/test/PMP_Boolean_operations/data-nm_boolop/nested_cubes_2.off new file mode 100644 index 000000000000..cb0bd77ae734 --- /dev/null +++ b/PMP_Boolean_operations/test/PMP_Boolean_operations/data-nm_boolop/nested_cubes_2.off @@ -0,0 +1,24 @@ +OFF +8 12 0 + +-1 -1 1 +1 -1 1 +1 -1 -1 +-1 -1 -1 +-1 1 -1 +-1 1 1 +1 1 1 +1 1 -1 +3 4 5 6 +3 1 0 3 +3 6 1 2 +3 0 1 6 +3 3 0 5 +3 7 2 3 +3 6 7 4 +3 3 2 1 +3 2 7 6 +3 6 5 0 +3 5 4 3 +3 3 4 7 + diff --git a/PMP_Boolean_operations/test/PMP_Boolean_operations/data-nm_boolop/pair_of_nm_cubes.off b/PMP_Boolean_operations/test/PMP_Boolean_operations/data-nm_boolop/pair_of_nm_cubes.off new file mode 100644 index 000000000000..6bf14d1d0ed4 --- /dev/null +++ b/PMP_Boolean_operations/test/PMP_Boolean_operations/data-nm_boolop/pair_of_nm_cubes.off @@ -0,0 +1,42 @@ +OFF +14 24 0 + +-0.5 -0.5 0.5 +0.5 -0.5 0.5 +0.5 -0.5 -0.5 +-0.5 -0.5 -0.5 +-0.5 0.5 -0.5 +-0.5 0.5 0.5 +0.5 0.5 0.5 +0.5 0.5 -0.5 +0.5 1.5 0.5 +1.5 1.5 0.5 +1.5 1.5 -0.5 +0.5 1.5 -0.5 +1.5 0.5 -0.5 +1.5 0.5 0.5 +3 4 5 6 +3 1 0 3 +3 6 1 2 +3 0 1 6 +3 3 0 5 +3 7 2 3 +3 9 10 11 +3 13 6 7 +3 9 13 12 +3 9 8 6 +3 8 11 7 +3 10 12 7 +3 6 7 4 +3 3 2 1 +3 2 7 6 +3 6 5 0 +3 5 4 3 +3 3 4 7 +3 11 8 9 +3 7 12 13 +3 12 10 9 +3 6 13 9 +3 7 6 8 +3 7 11 10 + diff --git a/PMP_Boolean_operations/test/PMP_Boolean_operations/data-nm_boolop/tangent_cube1.off b/PMP_Boolean_operations/test/PMP_Boolean_operations/data-nm_boolop/tangent_cube1.off new file mode 100644 index 000000000000..3d7bf7101998 --- /dev/null +++ b/PMP_Boolean_operations/test/PMP_Boolean_operations/data-nm_boolop/tangent_cube1.off @@ -0,0 +1,24 @@ +OFF +8 12 0 + +-0.5 -0.5 0.5 +0.5 -0.5 0.5 +0.5 -0.5 -0.5 +-0.5 -0.5 -0.5 +-0.5 0.5 -0.5 +-0.5 0.5 0.5 +0.5 0.5 0.5 +0.5 0.5 -0.5 +3 4 5 6 +3 1 0 3 +3 6 1 2 +3 0 1 6 +3 3 0 5 +3 7 2 3 +3 6 7 4 +3 3 2 1 +3 2 7 6 +3 6 5 0 +3 5 4 3 +3 3 4 7 + diff --git a/PMP_Boolean_operations/test/PMP_Boolean_operations/data-nm_boolop/tangent_cube2.off b/PMP_Boolean_operations/test/PMP_Boolean_operations/data-nm_boolop/tangent_cube2.off new file mode 100644 index 000000000000..5870a6f385f8 --- /dev/null +++ b/PMP_Boolean_operations/test/PMP_Boolean_operations/data-nm_boolop/tangent_cube2.off @@ -0,0 +1,24 @@ +OFF +8 12 0 + +-0.25 -0.25 1.5 +0.25 -0.25 1.5 +0.25 -0.25 0.5 +-0.25 -0.25 0.5 +-0.25 0.25 0.5 +-0.25 0.25 1.5 +0.25 0.25 1.5 +0.25 0.25 0.5 +3 4 5 6 +3 1 0 3 +3 6 1 2 +3 0 1 6 +3 3 0 5 +3 7 2 3 +3 6 7 4 +3 3 2 1 +3 2 7 6 +3 6 5 0 +3 5 4 3 +3 3 4 7 + diff --git a/PMP_Boolean_operations/test/PMP_Boolean_operations/self_union.cpp b/PMP_Boolean_operations/test/PMP_Boolean_operations/self_union.cpp new file mode 100644 index 000000000000..2c941b2d95b4 --- /dev/null +++ b/PMP_Boolean_operations/test/PMP_Boolean_operations/self_union.cpp @@ -0,0 +1,47 @@ +#include +#include +#include +#include + +#include + +#include + +#include + +#ifdef CGAL_LINKED_WITH_TBB +using Concurrency_tag = CGAL::Parallel_tag; +#else +using Concurrency_tag = CGAL::Sequential_tag; +#endif + +typedef CGAL::Exact_predicates_exact_constructions_kernel EPECK; +typedef EPECK::Point_3 Exact_point; +typedef CGAL::Surface_mesh Exact_mesh; +using Triangle = boost::container::small_vector; +namespace PMP = CGAL::Polygon_mesh_processing; +namespace params = CGAL::parameters; + +int main(int argc, char** argv) +{ + std::vector points; + std::vector triangles; + + auto filename = argc == 1 ? CGAL::data_file_path("meshes/elephant.off") : argv[1]; + + if (!CGAL::IO::read_polygon_soup(filename, points, triangles)) + { + std::cerr << "Error reading " << filename << "\n"; + return 1; + } + + CGAL::Real_timer timer; + timer.start(); + PMP::compute_self_union(points, triangles, CGAL::parameters::concurrency_tag(Concurrency_tag())); + timer.stop(); + std::cout << "Done in " << timer.time() << "\n"; + + CGAL::IO::write_OFF("self_union.off", points, triangles, params::stream_precision(17)); + + return 0; +} diff --git a/PMP_Boolean_operations/test/PMP_Boolean_operations/test_boolean_on_triangle_soup.cmd b/PMP_Boolean_operations/test/PMP_Boolean_operations/test_boolean_on_triangle_soup.cmd new file mode 100644 index 000000000000..ee48564603cd --- /dev/null +++ b/PMP_Boolean_operations/test/PMP_Boolean_operations/test_boolean_on_triangle_soup.cmd @@ -0,0 +1,6 @@ +data-nm_boolop/cube_0.off data-nm_boolop/cube_1.off +data-nm_boolop/cube3_minus_prism.off data-nm_boolop/cube4.off +data-nm_boolop/pair_of_nm_cubes.off data-nm_boolop/cube_2.off +data-nm_boolop/nested_cube1.off data-nm_boolop/nested_cube2.off +data-nm_boolop/nested_cube2.off data-nm_boolop/nested_cube1.off +data-nm_boolop/tangent_cube1.off data-nm_boolop/tangent_cube2.off diff --git a/PMP_Boolean_operations/test/PMP_Boolean_operations/test_boolean_on_triangle_soup.cpp b/PMP_Boolean_operations/test/PMP_Boolean_operations/test_boolean_on_triangle_soup.cpp new file mode 100644 index 000000000000..32dc568a3ef8 --- /dev/null +++ b/PMP_Boolean_operations/test/PMP_Boolean_operations/test_boolean_on_triangle_soup.cpp @@ -0,0 +1,77 @@ +#include +#include +#include +#include + +#include + +#include + +#include + +#ifdef CGAL_LINKED_WITH_TBB +using Concurrency_tag = CGAL::Parallel_tag; +#else +using Concurrency_tag = CGAL::Sequential_tag; +#endif + +typedef CGAL::Exact_predicates_exact_constructions_kernel EPECK; +typedef EPECK::Point_3 Exact_point; +typedef CGAL::Surface_mesh Exact_mesh; +namespace PMP = CGAL::Polygon_mesh_processing; +namespace params = CGAL::parameters; + + +void run(std::string f1, std::string f2) +{ + std::cout << "Running with " << f1 << " and " << f2 << "\n"; + using Triangle = boost::container::small_vector; + + std::vector points_1, points_2; + std::vector triangles_1, triangles_2; + + CGAL::IO::read_polygon_soup(f1, points_1, triangles_1); + CGAL::IO::read_polygon_soup(f2, points_2, triangles_2); + + auto bo_np = CGAL::parameters::concurrency_tag(Concurrency_tag()); + +{ + std::vector points_res; + std::vector triangles_res; + CGAL::Real_timer timer; + timer.start(); + PMP::compute_union(points_1, triangles_1, points_2, triangles_2, points_res, triangles_res, bo_np); + timer.stop(); + std::cout << "Done in " << timer.time() << "\n"; + + CGAL::IO::write_OFF("union.off", points_res, triangles_res, params::stream_precision(17)); +} +{ + std::vector points_res; + std::vector triangles_res; + CGAL::Real_timer timer; + timer.start(); + PMP::compute_intersection(points_1, triangles_1, points_2, triangles_2, points_res, triangles_res, bo_np); + timer.stop(); + std::cout << "Done in " << timer.time() << "\n"; + + CGAL::IO::write_OFF("intersection.off", points_res, triangles_res, params::stream_precision(17)); +} +{ + std::vector points_res; + std::vector triangles_res; + CGAL::Real_timer timer; + timer.start(); + PMP::compute_difference(points_1, triangles_1, points_2, triangles_2, points_res, triangles_res, bo_np); + timer.stop(); + std::cout << "Done in " << timer.time() << "\n"; + + CGAL::IO::write_OFF("difference.off", points_res, triangles_res, params::stream_precision(17)); +} +} + +int main(int argc, char** argv) +{ + for (int i=0; i< (argc-1)/2; ++i) + run(argv[2*i+1], argv[2*i+2]); +} diff --git a/PMP_Boolean_operations/test/PMP_Boolean_operations/test_corefinement_vs_autorefine.cmd b/PMP_Boolean_operations/test/PMP_Boolean_operations/test_corefinement_vs_autorefine.cmd new file mode 100644 index 000000000000..ff852e0e0235 --- /dev/null +++ b/PMP_Boolean_operations/test/PMP_Boolean_operations/test_corefinement_vs_autorefine.cmd @@ -0,0 +1,21 @@ +${CGAL_DATA_DIR}/meshes/elephant.off ${CGAL_DATA_DIR}/meshes/sphere.off +${CGAL_DATA_DIR}/meshes/open_cube.off data-coref/incompatible_with_open_cube.off +data-coref/floating_squares.off data-coref/hexa.off +data-coref/star_tgt1-0.off data-coref/star_tgt2-0.off +data-coref/star_tgt1-0.off data-coref/star_tgt2-1.off +data-coref/star_tgt1-0.off data-coref/star_tgt2-2.off +data-coref/star_tgt1-0.off data-coref/star_tgt2-3.off +data-coref/star_tgt1-0.off data-coref/star_tgt2-4.off +data-coref/star_tgt1-0.off data-coref/star_tgt2-5.off +data-coref/star_tgt1-1.off data-coref/star_tgt2-0.off +data-coref/star_tgt1-1.off data-coref/star_tgt2-1.off +data-coref/star_tgt1-1.off data-coref/star_tgt2-2.off +data-coref/star_tgt1-1.off data-coref/star_tgt2-3.off +data-coref/star_tgt1-1.off data-coref/star_tgt2-4.off +data-coref/star_tgt1-1.off data-coref/star_tgt2-5.off +data-coref/star_tgt1-2.off data-coref/star_tgt2-0.off +data-coref/star_tgt1-2.off data-coref/star_tgt2-1.off +data-coref/star_tgt1-2.off data-coref/star_tgt2-2.off +data-coref/star_tgt1-2.off data-coref/star_tgt2-3.off +data-coref/star_tgt1-2.off data-coref/star_tgt2-4.off +data-coref/star_tgt1-2.off data-coref/star_tgt2-5.off diff --git a/PMP_Boolean_operations/test/PMP_Boolean_operations/test_corefinement_vs_autorefine.cpp b/PMP_Boolean_operations/test/PMP_Boolean_operations/test_corefinement_vs_autorefine.cpp new file mode 100644 index 000000000000..4cb019011630 --- /dev/null +++ b/PMP_Boolean_operations/test/PMP_Boolean_operations/test_corefinement_vs_autorefine.cpp @@ -0,0 +1,132 @@ +#include +#include + +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +namespace PMP = CGAL::Polygon_mesh_processing; +namespace CFR = PMP::Corefinement; + +using Kernel = CGAL::Exact_predicates_exact_constructions_kernel; +using Point_3 = Kernel::Point_3; +using Mesh = CGAL::Surface_mesh; +using Polygon_3 = std::vector; +using Polygons = std::vector; +using Points = std::vector; +using face_descriptor = boost::graph_traits::face_descriptor; + +using Triangle = std::array; + +void test(const std::string& p_filename, const std::string& q_filename) +{ + Mesh P, Q; + + const std::filesystem::path p_path(p_filename); + const std::filesystem::path q_path(q_filename); + std::cout << "Testing " << p_path.filename().string() << " and " << q_path.filename().string() << '\n'; + + if (!CGAL::IO::read_polygon_mesh(p_filename, P)) + { + std::cerr << "Failed to read polygon mesh: " << p_filename << '\n'; + std::exit(EXIT_FAILURE); + } + if (!CGAL::IO::read_polygon_mesh(q_filename, Q)) + { + std::cerr << "Failed to read polygon mesh: " << q_filename << '\n'; + std::exit(EXIT_FAILURE); + } + + // Convert the inputs to triangle soups. + Points p_points, q_points; + Polygons p_polygons, q_polygons; + PMP::polygon_mesh_to_polygon_soup(P, p_points, p_polygons); + PMP::polygon_mesh_to_polygon_soup(Q, q_points, q_polygons); + + // Results computed directly on triangle soups. + const std::array names = {{ + "union", "intersection", "P - Q", "Q - P" + }}; + + std::array soup_points; + std::array soup_polygons; + + PMP::compute_union(p_points, p_polygons, q_points, q_polygons, soup_points[CFR::UNION], soup_polygons[CFR::UNION]); + PMP::compute_intersection(p_points, p_polygons, q_points, q_polygons, soup_points[CFR::INTERSECTION], soup_polygons[CFR::INTERSECTION]); + PMP::compute_difference(p_points, p_polygons, q_points, q_polygons, soup_points[CFR::TM1_MINUS_TM2], soup_polygons[CFR::TM1_MINUS_TM2]); + PMP::compute_difference(q_points, q_polygons, p_points, p_polygons, soup_points[CFR::TM2_MINUS_TM1], soup_polygons[CFR::TM2_MINUS_TM1]); + + // Results computed by the existing mesh API. + Mesh union_mesh, intersection_mesh, p_minus_q_mesh, q_minus_p_mesh; + std::array, 4> corefine_outputs{&union_mesh, &intersection_mesh, &p_minus_q_mesh, &q_minus_p_mesh}; + const std::array corefine_success = PMP::corefine_and_compute_boolean_operations(P, Q, corefine_outputs); + + for (std::size_t i = 0; i < 4; ++i) + { + if (!corefine_success[i] ){ + std::cout << names[i] << " is not manifold\n"; + continue; + } + std::cout << "Checking " << names[i] << "... "; + + // Convert the mesh result to a soup for comparison. + Mesh M; + PMP::duplicate_non_manifold_edges_in_polygon_soup(soup_points[i], soup_polygons[i]); + PMP::polygon_soup_to_polygon_mesh(soup_points[i], soup_polygons[i], M); + + std::vector m1_only; + std::vector m2_only; + PMP::match_faces(*(corefine_outputs[i].value()), M, boost::make_function_output_iterator([](const auto&) {}), + std::back_inserter(m1_only), + std::back_inserter(m2_only)); + const bool equal = m1_only.empty() && m2_only.empty(); + + if (!equal) + { + std::cerr << "FAILED: results are different\n"; + assert(false); + } + else + std::cout << " OK\n"; + } +} + +int main(int argc, char** argv) +{ + const std::filesystem::path program_path(argv[0]); + + if (argc < 3) + { + std::cerr << "Usage: " << program_path.filename().string() << " P.off Q.off\n"; + return EXIT_FAILURE; + } + + if(argc == 3) + { + test(argv[1], argv[2]); + return EXIT_SUCCESS; + } + + if ((argc-1)%2==0) + { + // cmd mode + for (int i=1; i