From 67a9aa8da01384a699b67c3d4cc47456fe48d2c4 Mon Sep 17 00:00:00 2001 From: Andreas Fabri Date: Fri, 12 Apr 2024 18:10:25 +0100 Subject: [PATCH] Fixed code and moved into a header --- .../boolean_constrained_plus.cpp | 196 +-------- .../include/CGAL/Triangulation_2/Boolean.h | 378 ++++++++++++++++++ 2 files changed, 397 insertions(+), 177 deletions(-) create mode 100644 Triangulation_2/include/CGAL/Triangulation_2/Boolean.h diff --git a/Triangulation_2/examples/Triangulation_2/boolean_constrained_plus.cpp b/Triangulation_2/examples/Triangulation_2/boolean_constrained_plus.cpp index d79fde5b05c6..c14ff562b177 100644 --- a/Triangulation_2/examples/Triangulation_2/boolean_constrained_plus.cpp +++ b/Triangulation_2/examples/Triangulation_2/boolean_constrained_plus.cpp @@ -3,214 +3,56 @@ #include #include -#include -#include -#include -#include -#include +#include +#include +#include #include -#include -#include + #include #include #include - -struct Boolean_cdt_2 { - - -struct FaceInfo { - - FaceInfo() - {} - - int nesting_level[2]; - - bool in_domain(int i) const - { - return nesting_level[i] % 2 == 1; - } - - template - bool in_domain(const Fct& fct) const - { - return fct(in_domain(0), in_domain(1)); - } - -}; - -typedef CGAL::Exact_predicates_exact_constructions_kernel K; -typedef K::Point_2 Point_2; -using Polygon_with_holes_2 = CGAL::Polygon_with_holes_2; -using Multipolygon_with_holes_2 = CGAL::Multipolygon_with_holes_2; - - -typedef CGAL::Exact_intersections_tag Itag; -typedef CGAL::Triangulation_vertex_base_2 Vb; -typedef CGAL::Triangulation_face_base_with_info_2 Fbb; -typedef CGAL::Constrained_triangulation_face_base_2 Fb; -typedef CGAL::Triangulation_data_structure_2 Tds; -typedef CGAL::Constrained_Delaunay_triangulation_2 CDT; -typedef CGAL::Constrained_triangulation_plus_2 CDTplus; -typedef CDTplus::Face_handle Face_handle; -typedef CDTplus::Constraint_id Constraint_id; -typedef CDTplus::Edge Edge; -typedef CDTplus::Context Context; - - - CDTplus cdt; - std::set idA, idB; - - Boolean_cdt_2() = default; - - void operator()(const Multipolygon_with_holes_2& pA, const Multipolygon_with_holes_2& pB) - { - for(const auto& pwh : pA.polygons_with_holes()){ - Constraint_id cidA = cdt.insert_constraint(pwh.outer_boundary().vertices_begin(), pwh.outer_boundary().vertices_end(), true); - idA.insert(cidA); - for(auto const& hole : pwh.holes()){ - cidA = cdt.insert_constraint(hole.vertices_begin(), hole.vertices_end(), true); - idA.insert(cidA); - } - } - - for(const auto& pwh : pB.polygons_with_holes()){ - Constraint_id cidB = cdt.insert_constraint(pwh.outer_boundary().vertices_begin(), pwh.outer_boundary().vertices_end(), true); - idB.insert(cidB); - for(auto const& hole : pwh.holes()){ - cidB = cdt.insert_constraint(hole.vertices_begin(), hole.vertices_end(), true); - idB.insert(cidB); - } - } - - mark_domains(idA, 0); - mark_domains(idB, 1); - } - -void -mark_domains(Face_handle start, - int index, - std::list& border, - const std::set& cids, - int aorb) -{ - if(start->info().nesting_level[aorb] != -1){ - return; - } - std::list queue; - queue.push_back(start); - - while(! queue.empty()){ - Face_handle fh = queue.front(); - queue.pop_front(); - if(fh->info().nesting_level[aorb] == -1){ - fh->info().nesting_level[aorb] = index; - for(int i = 0; i < 3; i++){ - Edge e(fh,i); - Face_handle n = fh->neighbor(i); - if(n->info().nesting_level[aorb] == -1){ - if(cdt.is_constrained(e)){ - for(Context c : cdt.contexts(e.first->vertex(cdt.cw(e.second)), - e.first->vertex(cdt.ccw(e.second)))){ - if(cids.find(c.id()) != cids.end()){ - border.push_back(e); - } - } - }else{ - queue.push_back(n); - } - } - } - } - } -} - - - - -void -mark_domains(const std::set& cids, int aorb) -{ - for(Face_handle f : cdt.all_face_handles()){ - f->info().nesting_level[aorb] = -1; - } - - std::list border; - mark_domains(cdt.infinite_face(), 0, border, cids, aorb); - while(! border.empty()){ - Edge e = border.front(); - border.pop_front(); - Face_handle n = e.first->neighbor(e.second); - if(n->info().nesting_level[aorb] == -1){ - mark_domains(n, e.first->info().nesting_level[aorb]+1, border, cids, aorb); - } - } -} - -template -void to_stl(const Fct& fct) -{ - std::ofstream out("cdt.stl"); - out.precision(17); - out << "solid outside\n"; - for(auto fit = cdt.finite_faces_begin(); fit != cdt.finite_faces_end(); ++fit){ - if(fit->info().in_domain(fct)){ - out << "facet normal 0 0 0\n" - << "outer loop\n"; - for(int i=0; i < 3; ++i){ - out << "vertex " << fit->vertex(i)->point() << " 0\n"; - } - out << "endloop\n" - << "endfacet\n"; - } - } - out << "endsolid" << std::endl; - out.close(); -} - -}; - using K = CGAL::Exact_predicates_exact_constructions_kernel; using Polygon_with_holes_2 = CGAL::Polygon_with_holes_2; using Multipolygon_with_holes_2 = CGAL::Multipolygon_with_holes_2; - - int main( ) { - - Boolean_cdt_2 bcdt; - + CGAL::Triangulations::Boolean bops; Multipolygon_with_holes_2 pA, pB; - { std::istringstream is("MULTIPOLYGON( ((0 0, 20 0, 20 30, 0 30, 0 0), (1 1, 1 2, 2 2, 2 1, 1 1 ) ) , (( 50 0, 60 0, 60 60, 50 60)) )"); + //std::istringstream is("MULTIPOLYGON( ((0 0, 2 0, 2 3, 0 3) ) )"); // (0.1 0.1, 0.1 0.4, 0.4 0.1) CGAL::IO::read_multi_polygon_WKT(is, pA); } { std::istringstream is("MULTIPOLYGON( ((10 1, 30 1, 30 2, 20 2, 20 4, 10 4)) )"); + //std::istringstream is("MULTIPOLYGON( ((2 1, 3 1, 3 2, 2 2)) "); CGAL::IO::read_multi_polygon_WKT(is, pB); } - bcdt(pA,pB); - - std::map map; - for(auto fh : bcdt.cdt.finite_face_handles()){ - map[fh] = fh->info().in_domain(0) || fh->info().in_domain(1); - } - - bcdt.to_stl([](bool a, bool b){ return a || b;} ); + bops.insert(pA,pB); + Multipolygon_with_holes_2 mpwh = bops([](bool a, bool b){ return a || b;}); + CGAL::IO::write_multi_polygon_WKT(std::cout, mpwh); + CGAL::draw(mpwh); - CGAL::draw(bcdt.cdt, boost::make_assoc_property_map(map)); + /* + std::map map; + for(auto fh : bops.cdt.finite_face_handles()){ + map[fh] = fh->info().in_domain(0) || fh->info().in_domain(1); + assert(map[fh] == (fh->info().label != 0)); + } + CGAL::draw(bops.cdt, boost::make_assoc_property_map(map)); + */ return 0; } diff --git a/Triangulation_2/include/CGAL/Triangulation_2/Boolean.h b/Triangulation_2/include/CGAL/Triangulation_2/Boolean.h new file mode 100644 index 000000000000..4886eace0714 --- /dev/null +++ b/Triangulation_2/include/CGAL/Triangulation_2/Boolean.h @@ -0,0 +1,378 @@ +// Copyright (c) 2024 GeometryFactory SARL (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) : Andreas Fabri + +#ifndef CGAL_TRIANGULATION_2_BOOLEAN_H +#define CGAL_TRIANGULATION_2_BOOLEAN_H + +#include + +#include + +#include +#include + +#include +#include +#include +#include + +#include +#include + +#include + +namespace CGAL { +namespace Triangulations { + +template +class Boolean { + +private: + struct FaceInfo { + + FaceInfo() + {} + + int label; + int nesting_level[2]; + bool processed; + + bool + in_domain(int i) const + { + return nesting_level[i] % 2 == 1; + } + + template + bool + in_domain(const Fct& fct) const + { + return fct(in_domain(0), in_domain(1)); + } + + }; + + using K = Kernel; + using Point_2 = typename K::Point_2; + using Polygon_2 = CGAL::Polygon_2; + using Polygon_with_holes_2 = CGAL::Polygon_with_holes_2; + using Multipolygon_with_holes_2 = CGAL::Multipolygon_with_holes_2; + + + using Itag = CGAL::Exact_intersections_tag; + using Vb = CGAL::Triangulation_vertex_base_2; + using Fbb = CGAL::Triangulation_face_base_with_info_2; + using Fb = CGAL::Constrained_triangulation_face_base_2; + using Tds = CGAL::Triangulation_data_structure_2; + using CDT = CGAL::Constrained_Delaunay_triangulation_2; + using CDTplus = CGAL::Constrained_triangulation_plus_2; + using Face_handle = typename CDTplus::Face_handle; + using Face_circulator = typename CDTplus::Face_circulator; + using Vertex_handle = typename CDTplus::Vertex_handle; + using Constraint_id = typename CDTplus::Constraint_id; + using Edge = typename CDTplus::Edge; + using Context = typename CDTplus::Context; + + + // @todo taken from Polygon_repair should be factorized + struct Polygon_less { + + bool + operator()(const Polygon_2& pa, const Polygon_2& pb) const + { + typename Polygon_2::Vertex_iterator va = pa.vertices_begin(); + typename Polygon_2::Vertex_iterator vb = pb.vertices_begin(); + while (va != pa.vertices_end() && vb != pb.vertices_end()) { + if (*va != *vb) return *va < *vb; + ++va; + ++vb; + } + if (vb == pb.vertices_end()) return false; + return true; + } + + }; + + + // @todo taken from Polygon_repair should be factorized + struct Polygon_with_holes_less { + Polygon_less pl; + + bool + operator()(const Polygon_with_holes_2& pa, const Polygon_with_holes_2& pb) const + { + if (pl(pa.outer_boundary(), pb.outer_boundary())) return true; + if (pl(pb.outer_boundary(), pa.outer_boundary())) return false; + typename Polygon_with_holes_2::Hole_const_iterator ha = pa.holes_begin(); + typename Polygon_with_holes_2::Hole_const_iterator hb = pb.holes_begin(); + while (ha != pa.holes_end() && hb != pb.holes_end()) { + if (pl(*ha, *hb)) return true; + if (pl(*hb, *ha)) return false; + } + if (hb == pb.holes_end()) return false; + return true; + } + + }; + + + CDTplus cdt; + std::set idA, idB; + +public: + + Boolean() = default; + + + void + insert(const Multipolygon_with_holes_2& pA, const Multipolygon_with_holes_2& pB) + { + for(const auto& pwh : pA.polygons_with_holes()){ + Constraint_id cidA = cdt.insert_constraint(pwh.outer_boundary().vertices_begin(), pwh.outer_boundary().vertices_end(), true); + idA.insert(cidA); + for(auto const& hole : pwh.holes()){ + cidA = cdt.insert_constraint(hole.vertices_begin(), hole.vertices_end(), true); + idA.insert(cidA); + } + } + + for(const auto& pwh : pB.polygons_with_holes()){ + Constraint_id cidB = cdt.insert_constraint(pwh.outer_boundary().vertices_begin(), pwh.outer_boundary().vertices_end(), true); + idB.insert(cidB); + for(auto const& hole : pwh.holes()){ + cidB = cdt.insert_constraint(hole.vertices_begin(), hole.vertices_end(), true); + idB.insert(cidB); + } + } + + mark_domains(idA, 0); + mark_domains(idB, 1); + } + +private: + + void + mark_domains(Face_handle start, + int index, + std::list& border, + const std::set& cids, + int aorb) + { + if(start->info().nesting_level[aorb] != -1){ + return; + } + std::list queue; + queue.push_back(start); + + while(! queue.empty()){ + Face_handle fh = queue.front(); + queue.pop_front(); + if(fh->info().nesting_level[aorb] == -1){ + fh->info().nesting_level[aorb] = index; + for(int i = 0; i < 3; i++){ + Edge e(fh,i); + Face_handle n = fh->neighbor(i); + if(n->info().nesting_level[aorb] == -1){ + if(cdt.is_constrained(e)){ + bool found = false; + for(Context c : cdt.contexts(e.first->vertex(cdt.cw(e.second)), + e.first->vertex(cdt.ccw(e.second)))){ + if(cids.find(c.id()) != cids.end()){ + found = true; + break; + } + } + if (found) { + border.push_back(e); + } else { + queue.push_back(n); + } + }else{ + queue.push_back(n); + } + } + } + } + } + } + + + // this marks the domains for either the first or the second multipolygon + void + mark_domains(const std::set& cids, int aorb) + { + for(Face_handle f : cdt.all_face_handles()){ + f->info().nesting_level[aorb] = -1; + } + + std::list border; + mark_domains(cdt.infinite_face(), 0, border, cids, aorb); + + while(! border.empty()){ + Edge e = border.front(); + border.pop_front(); + Face_handle n = e.first->neighbor(e.second); + if(n->info().nesting_level[aorb] == -1){ + mark_domains(n, e.first->info().nesting_level[aorb]+1, border, cids, aorb); + } + } + } + + + template + void + label_domains(Face_handle start, int label, const Fct& fct) + { + std::list queue; + start->info().label = label; + queue.push_back(start); + + while(! queue.empty()){ + Face_handle fh = queue.front(); + queue.pop_front(); + + for(int i = 0; i < 3; i++){ + Face_handle n = fh->neighbor(i); + if(n->info().in_domain(fct)){ + if(n->info().label == 0){ + n->info().label = label; + queue.push_back(n); + } + } + } + } + } + + // this marks the domain for the Boolean operation applied on the two multipolygons + template + int + label_domains(const Fct& fct) + { + for (auto const face: cdt.all_face_handles()) { + face->info().processed = false; + face->info().label = 0; + } + int label = 1; + for (auto const face: cdt.all_face_handles()) { + if(face->info().in_domain(fct) && face->info().label == 0){ + label_domains(face, label, fct); + ++label; + } + } + return label; + } + + + + // @todo taken from Polygon_repair and adapted; might be factorized + // Reconstruct ring boundary starting from an edge (face + opposite vertex) that is part of it + template + void + reconstruct_ring(std::list& ring, + Face_handle face_adjacent_to_boundary, + int opposite_vertex, + const Fct& fct) + { + // Create ring + Face_handle current_face = face_adjacent_to_boundary; + int current_opposite_vertex = opposite_vertex; + do { + CGAL_assertion(current_face->info().in_domain(fct) == face_adjacent_to_boundary->info().in_domain(fct)); + current_face->info().processed = true; + Vertex_handle pivot_vertex = current_face->vertex(current_face->cw(current_opposite_vertex)); + // std::cout << "\tAdding point " << pivot_vertex->point() << std::endl; + ring.push_back(pivot_vertex->point()); + Face_circulator fc = cdt.incident_faces(pivot_vertex, current_face); + do { + ++fc; + } while (fc->info().in_domain(fct) != current_face->info().in_domain(fct)); + current_face = fc; + current_opposite_vertex = fc->cw(fc->index(pivot_vertex)); + } while (current_face != face_adjacent_to_boundary || + current_opposite_vertex != opposite_vertex); + + // Start at lexicographically smallest vertex + typename std::list::iterator smallest_vertex = ring.begin(); + for (typename std::list::iterator current_vertex = ring.begin(); + current_vertex != ring.end(); ++current_vertex) { + if (*current_vertex < *smallest_vertex) smallest_vertex = current_vertex; + } + if (ring.front() != *smallest_vertex) { + ring.splice(ring.begin(), ring, smallest_vertex, ring.end()); + } + } + + +public: + + // @todo taken from Polygon_repair and adapted; might be factorized + template + Multipolygon_with_holes_2 + operator()(const Fct& fct) + { + int number_of_polygons = label_domains(fct) - 1; + + Multipolygon_with_holes_2 mp; + std::vector polygons; // outer boundaries + std::vector> holes; // holes are ordered (per polygon) + polygons.resize(number_of_polygons); + holes.resize(number_of_polygons); + + for (auto const face: cdt.all_face_handles()) { + face->info().processed = false; + } + for (auto const &face: cdt.finite_face_handles()) { + if (! face->info().in_domain(fct)) continue; // exterior triangle + if (face->info().processed) continue; // already reconstructed + for (int opposite_vertex = 0; opposite_vertex < 3; ++opposite_vertex) { + if (face->info().in_domain(fct) == face->neighbor(opposite_vertex)->info().in_domain(fct)) continue; // not adjacent to boundary + + // Reconstruct ring + // @todo less copying + std::list ring; + reconstruct_ring(ring, face, opposite_vertex, fct); + + // Put ring in polygons + Polygon_2 polygon(ring.begin(), ring.end()); + if (polygon.orientation() == CGAL::COUNTERCLOCKWISE) { + polygons[face->info().label-1] = polygon; + } else { + holes[face->info().label-1].insert(polygon); + } break; + } + } + + // Create polygons with holes and put in multipolygon + std::set ordered_polygons; + for (std::size_t i = 0; i < polygons.size(); ++i) { + ordered_polygons.insert(Polygon_with_holes_2(polygons[i], holes[i].begin(), holes[i].end())); + } + for (auto const& polygon: ordered_polygons) { + mp.add_polygon_with_holes(polygon); + } + return mp; + } + + + const CDTplus& + triangulation() const + { + return cdt; + } + +}; + + +} // namespace Triangulations +} //namespace CGAL + +#endif CGAL_TRIANGULATION_2_BOOLEAN_H