1 // Copyright (c) 2018 Liangliang Nan. All rights reserved. 2 // 3 // This file is part of CGAL (www.cgal.org) 4 // 5 // $URL: https://github.com/CGAL/cgal/blob/v5.3/Polygonal_surface_reconstruction/include/CGAL/internal/hypothesis.h $ 6 // $Id: hypothesis.h 6b87fe3 2020-12-05T11:11:33+01:00 Mael Rouxel-Labbé 7 // SPDX-License-Identifier: LGPL-3.0-or-later OR LicenseRef-Commercial 8 // 9 // Author(s) : Liangliang Nan 10 11 #ifndef CGAL_POLYGONAL_SURFACE_RECONSTRUCTION_HYPOTHESIS_H 12 #define CGAL_POLYGONAL_SURFACE_RECONSTRUCTION_HYPOTHESIS_H 13 14 #include <CGAL/license/Polygonal_surface_reconstruction.h> 15 16 #include <CGAL/Surface_mesh.h> 17 #include <CGAL/Projection_traits_xy_3.h> 18 #include <CGAL/convex_hull_2.h> 19 #include <CGAL/bounding_box.h> 20 #include <CGAL/intersections.h> 21 #include <CGAL/assertions.h> 22 #include <CGAL/internal/parameters.h> 23 #include <CGAL/internal/point_set_with_planes.h> 24 25 #include <set> 26 #include <unordered_map> 27 28 29 namespace CGAL { 30 31 namespace internal { 32 33 /** 34 * Generates candidate faces by pairwise intersecting of the supporting planes of 35 * the planar segments. 36 */ 37 38 template <typename Kernel> 39 class Hypothesis 40 { 41 private: 42 typedef typename Kernel::FT FT; 43 typedef typename Kernel::Point_3 Point; 44 typedef typename Kernel::Point_2 Point2; 45 typedef typename Kernel::Vector_3 Vector; 46 typedef typename Kernel::Line_3 Line; 47 typedef typename Kernel::Segment_3 Segment; 48 typedef typename Kernel::Plane_3 Plane; 49 typedef internal::Planar_segment<Kernel> Planar_segment; 50 typedef internal::Point_set_with_planes<Kernel> Point_set_with_planes; 51 52 typedef CGAL::Surface_mesh<Point> Polygon_mesh; 53 typedef typename Polygon_mesh::Face_index Face_descriptor; 54 typedef typename Polygon_mesh::Edge_index Edge_descriptor; 55 typedef typename Polygon_mesh::Vertex_index Vertex_descriptor; 56 57 typedef typename Polygon_mesh::Halfedge_index Halfedge_descriptor; 58 59 public: 60 Hypothesis(); 61 ~Hypothesis(); 62 63 void generate(Point_set_with_planes& point_set, Polygon_mesh& candidate_faces); 64 65 /// 'Intersection' represents a set of faces intersecting at a common edge. 66 /// \note The faces are represented by their halfedges. 67 struct Intersection : public std::vector<Halfedge_descriptor> { 68 const Point* s; 69 const Point* t; 70 }; 71 typedef typename std::vector<Intersection> Adjacency; 72 /// Extracts the adjacency of the pairwise intersection. 73 /// The extracted adjacency will be used to formulate the hard constraints 74 /// in the face selection stage. 75 Adjacency extract_adjacency(const Polygon_mesh& candidate_faces); 76 77 private: 78 // Merges near co-planar segments 79 void refine_planes(); 80 81 // Constructs a mesh representing the bounding box of the point set 82 void construct_bbox_mesh(Polygon_mesh& bbox_mesh); 83 84 // Construct a mesh from the segments bounded by the bounding box mesh 85 void construct_proxy_mesh(const Polygon_mesh& bbox_mesh, Polygon_mesh& candidate_faces); 86 87 // Pairwise intersection 88 void pairwise_intersection(Polygon_mesh& candidate_faces); 89 90 // Counts the number of points that are with the dist_threshold to its supporting plane 91 std::size_t number_of_points_on_plane(const Planar_segment* s, const Plane* plane, FT dist_threshold); 92 93 // Merges two planar segments; 94 void merge(Planar_segment* s1, Planar_segment* s2); 95 96 // Pre-computes all potential intersections of plane triplets 97 void compute_triplet_intersections(); 98 99 // Queries the intersecting point for a plane triplet 100 const Point* query_intersection(const Plane* plane1, const Plane* plane2, const Plane* plane3); 101 102 bool halfedge_exists(Vertex_descriptor v1, Vertex_descriptor v2, const Polygon_mesh& mesh); 103 104 // Tests if 'face' insects 'plane' 105 bool do_intersect(const Polygon_mesh& mesh, Face_descriptor face, const Plane* plane); 106 107 // Cuts face using the cutting_plane and returns the new faces 108 std::vector<Face_descriptor> cut(Face_descriptor face, const Plane* cutting_plane, Polygon_mesh& mesh); 109 110 // Collects all faces in 'mesh' that intersect 'face'. Store in std::set() so easier to erase 111 std::set<Face_descriptor> collect_intersecting_faces(Face_descriptor face, const Polygon_mesh& mesh); 112 113 // Represents an intersecting point at an edge 114 struct EdgePos { EdgePosEdgePos115 EdgePos(Edge_descriptor e, const Point* p) : edge(e), pos(p) {} 116 Edge_descriptor edge; 117 const Point* pos; 118 }; 119 120 // Computes the intersecting points of face and cutting_plane. The intersecting points are returned 121 // by 'existing_vts' (if the plane intersects the face at its vertices) and 'new_vts' (if the plane 122 // intersects the face at its edges). 123 void compute_intersections(const Polygon_mesh& mesh, 124 Face_descriptor face, const Plane* cutting_plane, 125 std::vector<Vertex_descriptor>& existing_vts, 126 std::vector<EdgePos>& new_vts 127 ); 128 129 // This function will 130 // - split an edge denoted by 'ep' 131 // - assign the new edges the supporting faces 132 // - return the halfedge pointing to the new vertex 133 // Internally it uses Euler split_edge(). 134 Halfedge_descriptor split_edge(Polygon_mesh& mesh, const EdgePos& ep, const Plane* cutting_plane); 135 136 // Clears cached intermediate results 137 void clear(); 138 139 private: 140 // The input point cloud with planes 141 Point_set_with_planes * point_set_; 142 143 // The intersection of the planes can be unreliable when the planes are near parallel. 144 // Here are the tricks we use in our implementation: 145 // - We first test if an intersection exists for every pair of planes. We then collect 146 // plane triplets such that every pair in the plane triplet intersect. This is achieved 147 // by testing each plane against the known intersecting pairs. 148 // - The 3D vertices of the final faces are obtained by computing the intersections of 149 // the plane triplets. To cope with limited floating point precision, each vertex is 150 // identified by the pointers of (in an increasing order) of the three planes from 151 // which it is computed. By doing so, two vertices with almost identical positions can 152 // be distinguished. This turned out to be quite robust in handling very close and near 153 // parallel planes. 154 155 // The supporting planes of all planar segments and the bounding box faces 156 std::vector<const Plane*> supporting_planes_; 157 158 // Precomputed intersecting points of all plane triplets 159 std::vector<const Point*> intersecting_points_; 160 161 typedef typename std::unordered_map<const Plane*, const Point*> Plane_to_point_map; 162 typedef typename std::unordered_map<const Plane*, Plane_to_point_map> Two_planes_to_point_map; 163 typedef typename std::unordered_map<const Plane*, Two_planes_to_point_map> Planes_to_point_map; 164 Planes_to_point_map triplet_intersections_; 165 }; 166 167 168 ////////////////////////////////////////////////////////////////////////// 169 170 // implementation 171 172 173 template <typename Kernel> Hypothesis()174 Hypothesis<Kernel>::Hypothesis() 175 { 176 } 177 178 179 template <typename Kernel> ~Hypothesis()180 Hypothesis<Kernel>::~Hypothesis() 181 { 182 clear(); 183 } 184 185 186 template <typename Kernel> clear()187 void Hypothesis<Kernel>::clear() { 188 for (std::size_t i = 0; i < supporting_planes_.size(); ++i) 189 delete supporting_planes_[i]; 190 supporting_planes_.clear(); 191 192 for (std::size_t i = 0; i < intersecting_points_.size(); ++i) 193 delete intersecting_points_[i]; 194 intersecting_points_.clear(); 195 196 triplet_intersections_.clear(); 197 } 198 199 200 template <typename Kernel> generate(Point_set_with_planes & point_set,Polygon_mesh & candidate_faces)201 void Hypothesis<Kernel>::generate(Point_set_with_planes& point_set, Polygon_mesh& candidate_faces) { 202 point_set_ = &point_set; 203 204 refine_planes(); 205 206 Polygon_mesh bbox_mesh; 207 construct_bbox_mesh(bbox_mesh); 208 209 construct_proxy_mesh(bbox_mesh, candidate_faces); 210 211 pairwise_intersection(candidate_faces); 212 } 213 214 215 /// \cond SKIP_IN_MANUAL 216 217 template <typename Planar_segment> 218 class SegmentSizeIncreasing 219 { 220 public: SegmentSizeIncreasing()221 SegmentSizeIncreasing() {} operator()222 bool operator()(const Planar_segment* s0, const Planar_segment* s1) const { 223 return s0->size() < s1->size(); 224 } 225 }; 226 227 template <typename Planar_segment> 228 class SegmentSizeDecreasing 229 { 230 public: SegmentSizeDecreasing()231 SegmentSizeDecreasing() {} operator()232 bool operator()(const Planar_segment* s0, const Planar_segment* s1) const { 233 return s0->size() > s1->size(); 234 } 235 }; 236 237 template <typename FT, typename Vector> normalize(Vector & v)238 void normalize(Vector& v) { 239 FT s = std::sqrt(v.squared_length()); 240 if (s > 1e-30) 241 s = FT(1) / s; 242 v *= s; 243 } 244 245 template <typename BBox> bbox_radius(const BBox & box)246 typename BBox::FT bbox_radius(const BBox& box) { 247 typedef typename BBox::FT FT; 248 FT dx = box.xmax() - box.xmin(); 249 FT dy = box.ymax() - box.ymin(); 250 FT dz = box.zmax() - box.zmin(); 251 return FT(0.5) * std::sqrt(dx * dx + dy * dy + dz * dz); 252 } 253 254 // Computes the intersection of a plane triplet 255 // Returns true if the intersection exists (p returns the point) 256 template <typename Plane, typename Point> intersect_plane_triplet(const Plane * plane1,const Plane * plane2,const Plane * plane3,Point & p)257 bool intersect_plane_triplet(const Plane* plane1, const Plane* plane2, const Plane* plane3, Point& p) { 258 if (plane1 == plane2 || plane1 == plane3 || plane2 == plane3) 259 return false; 260 261 CGAL::Object obj = CGAL::intersection(*plane1, *plane2, *plane3); 262 263 // pt is the intersection point of the 3 planes 264 if (const Point* pt = CGAL::object_cast<Point>(&obj)) { 265 p = *pt; 266 return true; 267 } 268 else { 269 // If reached here, the reason might be: 270 // (1) two or more are parallel; 271 // (2) they intersect at the same line 272 // We can simply ignore these cases 273 return false; 274 } 275 } 276 277 278 template <typename VT> sort_increasing(VT & v1,VT & v2,VT & v3)279 void sort_increasing(VT& v1, VT& v2, VT& v3) { 280 VT vmin = 0; 281 if (v1 < v2 && v1 < v3) 282 vmin = v1; 283 else if (v2 < v1 && v2 < v3) 284 vmin = v2; 285 else 286 vmin = v3; 287 288 VT vmid = 0; 289 if ((v1 > v2 && v1 < v3) || (v1 < v2 && v1 > v3)) 290 vmid = v1; 291 else if ((v2 > v1 && v2 < v3) || (v2 < v1 && v2 > v3)) 292 vmid = v2; 293 else 294 vmid = v3; 295 296 VT vmax = 0; 297 if (v1 > v2 && v1 > v3) 298 vmax = v1; 299 else if (v2 > v1 && v2 > v3) 300 vmax = v2; 301 else 302 vmax = v3; 303 304 v1 = vmin; 305 v2 = vmid; 306 v3 = vmax; 307 } 308 /// \endcond 309 310 311 template <typename Kernel> number_of_points_on_plane(const Planar_segment * s,const Plane * plane,FT dist_threshold)312 std::size_t Hypothesis<Kernel>::number_of_points_on_plane(const Planar_segment* s, const Plane* plane, FT dist_threshold) { 313 CGAL_assertion(const_cast<Planar_segment*>(s)->point_set() == point_set_); 314 315 std::size_t count = 0; 316 const typename Point_set_with_planes::Point_map& points = point_set_->point_map(); 317 for (std::size_t i = 0; i < s->size(); ++i) { 318 std::size_t idx = s->at(i); 319 const Point& p = points[idx]; 320 321 FT sdist = CGAL::squared_distance(*plane, p); 322 FT dist = std::sqrt(sdist); 323 if (dist < dist_threshold) 324 ++count; 325 } 326 return count; 327 } 328 329 template <typename Kernel> merge(Planar_segment * s1,Planar_segment * s2)330 void Hypothesis<Kernel>::merge(Planar_segment* s1, Planar_segment* s2) { 331 CGAL_assertion(const_cast<Planar_segment*>(s1)->point_set() == point_set_); 332 CGAL_assertion(const_cast<Planar_segment*>(s2)->point_set() == point_set_); 333 std::vector< Planar_segment* >& segments = point_set_->planar_segments(); 334 335 std::vector<std::size_t> points_indices; 336 points_indices.insert(points_indices.end(), s1->begin(), s1->end()); 337 points_indices.insert(points_indices.end(), s2->begin(), s2->end()); 338 339 Planar_segment* s = new Planar_segment(point_set_); 340 s->insert(s->end(), points_indices.begin(), points_indices.end()); 341 s->fit_supporting_plane(); 342 segments.push_back(s); 343 344 typename std::vector< Planar_segment* >::iterator pos = std::find(segments.begin(), segments.end(), s1); 345 if (pos != segments.end()) { 346 Planar_segment* tmp = *pos; 347 const Plane* plane = tmp->supporting_plane(); 348 segments.erase(pos); 349 delete tmp; 350 delete plane; 351 } 352 else 353 std::cerr << "Fatal error: should not reach here" << std::endl; 354 355 pos = std::find(segments.begin(), segments.end(), s2); 356 if (pos != segments.end()) { 357 Planar_segment* tmp = *pos; 358 const Plane* plane = tmp->supporting_plane(); 359 segments.erase(pos); 360 delete tmp; 361 delete plane; 362 } 363 else 364 std::cerr << "Fatal error: should not reach here" << std::endl; 365 } 366 367 template <typename Kernel> refine_planes()368 void Hypothesis<Kernel>::refine_planes() { 369 std::vector< Planar_segment* >& segments = point_set_->planar_segments(); 370 const typename Point_set_with_planes::Point_map& points = point_set_->point_map(); 371 372 FT avg_max_dist = 0; 373 for (std::size_t i = 0; i < segments.size(); ++i) { 374 Planar_segment* s = segments[i]; 375 const Plane* plane = s->fit_supporting_plane(); // user may provide invalid plane fitting (we always fit) 376 377 FT max_dist = -(std::numeric_limits<FT>::max)(); 378 for (std::size_t j = 0; j < s->size(); ++j) { 379 std::size_t idx = s->at(j); 380 const Point& p = points[idx]; 381 FT sdist = CGAL::squared_distance(*plane, p); 382 max_dist = (std::max)(max_dist, std::sqrt(sdist)); 383 } 384 385 avg_max_dist += max_dist; 386 } 387 avg_max_dist /= segments.size(); 388 avg_max_dist /= FT(2.0); 389 390 FT theta = static_cast<FT>(CGAL_PI * 10.0 / FT(180.0)); // in radian 391 bool merged = false; 392 do { 393 merged = false; 394 // Segments with less points have less confidences and thus should be merged first. 395 // So we sort the segments according to their sizes. 396 std::sort(segments.begin(), segments.end(), internal::SegmentSizeIncreasing<Planar_segment>()); 397 398 for (std::size_t i = 0; i < segments.size(); ++i) { 399 Planar_segment* s1 = segments[i]; 400 const Plane* plane1 = s1->supporting_plane(); 401 Vector n1 = plane1->orthogonal_vector(); 402 internal::normalize<FT, Vector>(n1); 403 404 FT num_threshold = s1->size() / FT(5.0); 405 for (std::size_t j = i + 1; j < segments.size(); ++j) { 406 Planar_segment* s2 = segments[j]; 407 const Plane* plane2 = s2->supporting_plane(); 408 Vector n2 = plane2->orthogonal_vector(); 409 internal::normalize<FT, Vector>(n2); 410 411 if (std::abs(n1 * n2) > std::cos(theta)) { 412 std::size_t set1on2 = number_of_points_on_plane(s1, plane2, avg_max_dist); 413 std::size_t set2on1 = number_of_points_on_plane(s2, plane1, avg_max_dist); 414 if (set1on2 > num_threshold || set2on1 > num_threshold) { 415 merge(s1, s2); 416 merged = true; 417 break; 418 } 419 } 420 } 421 if (merged) 422 break; 423 } 424 } while (merged); 425 426 std::sort(segments.begin(), segments.end(), internal::SegmentSizeDecreasing<Planar_segment>()); 427 428 // Stores all the supporting planes 429 for (std::size_t i = 0; i < segments.size(); ++i) { 430 Planar_segment* s = segments[i]; 431 const Plane* plane = s->supporting_plane(); 432 supporting_planes_.push_back(plane); 433 } 434 } 435 436 437 template <typename Kernel> construct_bbox_mesh(Polygon_mesh & mesh)438 void Hypothesis<Kernel>::construct_bbox_mesh(Polygon_mesh& mesh) { 439 const typename Point_set_with_planes::Point_map& points = point_set_->point_map(); 440 441 typedef typename Kernel::Iso_cuboid_3 BBox; 442 const BBox& box = CGAL::bounding_box(points.begin(), points.end()); 443 444 FT dx = box.xmax() - box.xmin(); 445 FT dy = box.ymax() - box.ymin(); 446 FT dz = box.zmax() - box.zmin(); 447 FT radius = FT(0.5) * std::sqrt(dx * dx + dy * dy + dz * dz); 448 FT offset = radius * FT(0.05); 449 450 // make the box larger to ensure all points are enclosed. 451 FT xmin = box.xmin() - offset, xmax = box.xmax() + offset; 452 FT ymin = box.ymin() - offset, ymax = box.ymax() + offset; 453 FT zmin = box.zmin() - offset, zmax = box.zmax() + offset; 454 455 mesh.clear(); 456 457 Vertex_descriptor v0 = mesh.add_vertex(Point(xmin, ymin, zmin)); // 0 458 Vertex_descriptor v1 = mesh.add_vertex(Point(xmax, ymin, zmin)); // 1 459 Vertex_descriptor v2 = mesh.add_vertex(Point(xmax, ymin, zmax)); // 2 460 Vertex_descriptor v3 = mesh.add_vertex(Point(xmin, ymin, zmax)); // 3 461 Vertex_descriptor v4 = mesh.add_vertex(Point(xmax, ymax, zmax)); // 4 462 Vertex_descriptor v5 = mesh.add_vertex(Point(xmax, ymax, zmin)); // 5 463 Vertex_descriptor v6 = mesh.add_vertex(Point(xmin, ymax, zmin)); // 6 464 Vertex_descriptor v7 = mesh.add_vertex(Point(xmin, ymax, zmax)); // 7 465 466 mesh.add_face(v0, v1, v2, v3); 467 mesh.add_face(v1, v5, v4, v2); 468 mesh.add_face(v1, v0, v6, v5); 469 mesh.add_face(v4, v5, v6, v7); 470 mesh.add_face(v0, v3, v7, v6); 471 mesh.add_face(v2, v4, v7, v3); 472 473 // The supporting plane of each face 474 typename Polygon_mesh::template Property_map<Face_descriptor, const Plane*> face_supporting_planes = 475 mesh.template add_property_map<Face_descriptor, const Plane*>("f:supp_plane").first; 476 477 // The supporting planes of each edge 478 typename Polygon_mesh::template Property_map<Edge_descriptor, std::set<const Plane*> > edge_supporting_planes 479 = mesh.template add_property_map<Edge_descriptor, std::set<const Plane*> >("e:supp_plane").first; 480 481 // The supporting planes of each vertex 482 typename Polygon_mesh::template Property_map<Vertex_descriptor, std::set<const Plane*> > vertex_supporting_planes 483 = mesh.template add_property_map<Vertex_descriptor, std::set<const Plane*> >("v:supp_plane").first; 484 485 // Assigns the original plane for each face 486 const typename Polygon_mesh::template Property_map<Vertex_descriptor, Point>& coords = mesh.points(); 487 for(auto fd : mesh.faces()) { 488 Halfedge_descriptor h = mesh.halfedge(fd); 489 Vertex_descriptor va = mesh.target(h); const Point& pa = coords[va]; h = mesh.next(h); 490 Vertex_descriptor vb = mesh.target(h); const Point& pb = coords[vb]; h = mesh.next(h); 491 Vertex_descriptor vc = mesh.target(h); const Point& pc = coords[vc]; 492 const Plane* plane = new Plane(pa, pb, pc); 493 supporting_planes_.push_back(plane); 494 face_supporting_planes[fd] = plane; 495 } 496 497 // Assigns the original planes for each edge 498 for( auto ed : mesh.edges()) { 499 Halfedge_descriptor h1 = mesh.halfedge(ed); 500 Halfedge_descriptor h2 = mesh.opposite(h1); 501 502 Face_descriptor f1 = mesh.face(h1); 503 Face_descriptor f2 = mesh.face(h2); 504 CGAL_assertion(f1 != Polygon_mesh::null_face()); // the bbox mesh is closed 505 CGAL_assertion(f2 != Polygon_mesh::null_face()); // the bbox mesh is closed 506 507 const Plane* plane1 = face_supporting_planes[f1]; 508 const Plane* plane2 = face_supporting_planes[f2]; 509 CGAL_assertion(plane1 && plane2 && plane1 != plane2); 510 511 edge_supporting_planes[ed].insert(plane1); 512 edge_supporting_planes[ed].insert(plane2); 513 CGAL_assertion(edge_supporting_planes[ed].size() == 2); 514 } 515 516 // Assigns the original planes for each vertex 517 for(auto vd : mesh.vertices()) { 518 CGAL_assertion(vertex_supporting_planes[vd].size() == 0); 519 CGAL::Halfedge_around_target_circulator<Polygon_mesh> hbegin(vd, mesh), done(hbegin); 520 do { 521 Halfedge_descriptor h = *hbegin; 522 Face_descriptor f = mesh.face(h); 523 const Plane* plane = face_supporting_planes[f]; 524 vertex_supporting_planes[vd].insert(plane); 525 ++hbegin; 526 } while (hbegin != done); 527 CGAL_assertion(vertex_supporting_planes[vd].size() == 3); 528 } 529 530 std::sort(supporting_planes_.begin(), supporting_planes_.end()); 531 532 CGAL_assertion(mesh.is_valid()); 533 } 534 535 536 template <typename Kernel> construct_proxy_mesh(const Polygon_mesh & bbox_mesh,Polygon_mesh & candidate_faces)537 void Hypothesis<Kernel>::construct_proxy_mesh(const Polygon_mesh& bbox_mesh, Polygon_mesh& candidate_faces) { 538 539 // Properties of the bbox_mesh 540 541 typename Polygon_mesh::template Property_map<Edge_descriptor, std::set<const Plane*> > bbox_edge_supporting_planes 542 = bbox_mesh.template property_map<Edge_descriptor, std::set<const Plane*> >("e:supp_plane").first; 543 typename Polygon_mesh::template Property_map<Vertex_descriptor, std::set<const Plane*> > bbox_vertex_supporting_planes 544 = bbox_mesh.template property_map<Vertex_descriptor, std::set<const Plane*> >("v:supp_plane").first; 545 546 // The properties of the proxy mesh 547 candidate_faces.clear(); 548 549 // The supporting plane of each face 550 typename Polygon_mesh::template Property_map<Face_descriptor, const Plane*> face_supporting_planes 551 = candidate_faces.template add_property_map<Face_descriptor, const Plane*>("f:supp_plane").first; 552 553 // The supporting planar segment of each face 554 typename Polygon_mesh::template Property_map<Face_descriptor, Planar_segment*> face_supporting_segments 555 = candidate_faces.template add_property_map<Face_descriptor, Planar_segment*>("f:supp_segment").first; 556 557 // The supporting planes of each edge 558 typename Polygon_mesh::template Property_map<Edge_descriptor, std::set<const Plane*> > edge_supporting_planes 559 = candidate_faces.template add_property_map<Edge_descriptor, std::set<const Plane*> >("e:supp_plane").first; 560 561 // The supporting planes of each vertex 562 typename Polygon_mesh::template Property_map<Vertex_descriptor, std::set<const Plane*> > vertex_supporting_planes 563 = candidate_faces.template add_property_map<Vertex_descriptor, std::set<const Plane*> >("v:supp_plane").first; 564 565 const std::vector<Planar_segment*>& segments = point_set_->planar_segments(); 566 const typename Polygon_mesh::template Property_map<Vertex_descriptor, Point>& coords = bbox_mesh.points(); 567 568 for (std::size_t i = 0; i < segments.size(); ++i) { 569 Planar_segment* g = segments[i]; 570 const Plane* cutting_plane = g->supporting_plane(); 571 572 std::vector<Point> intersecting_points; 573 std::vector< std::set<const Plane*> > intersecting_points_source_planes; 574 575 for(auto ed : bbox_mesh.edges()) { 576 Vertex_descriptor sd = bbox_mesh.vertex(ed, 0); 577 Vertex_descriptor td = bbox_mesh.vertex(ed, 1); 578 const Point& s = coords[sd]; 579 const Point& t = coords[td]; 580 581 CGAL::Oriented_side ss = cutting_plane->oriented_side(s); 582 CGAL::Oriented_side st = cutting_plane->oriented_side(t); 583 584 if ((ss == CGAL::ON_POSITIVE_SIDE && st == CGAL::ON_NEGATIVE_SIDE) || (ss == ON_NEGATIVE_SIDE && st == CGAL::ON_POSITIVE_SIDE)) { 585 CGAL::Object obj = CGAL::intersection(*cutting_plane, Line(s, t)); 586 if (const Point* p = CGAL::object_cast<Point>(&obj)) { 587 intersecting_points.push_back(*p); 588 std::set<const Plane*> planes = bbox_edge_supporting_planes[ed]; 589 590 planes.insert(cutting_plane); 591 CGAL_assertion(planes.size() == 3); 592 intersecting_points_source_planes.push_back(planes); 593 } 594 else 595 std::cerr << "Fatal error: should have intersection" << std::endl; 596 } 597 else { 598 if (ss == CGAL::ON_ORIENTED_BOUNDARY && st != CGAL::ON_ORIENTED_BOUNDARY) { 599 intersecting_points.push_back(s); 600 const std::set<const Plane*>& planes = bbox_vertex_supporting_planes[sd]; 601 CGAL_assertion(planes.size() == 3); 602 intersecting_points_source_planes.push_back(planes); 603 } 604 else if (st == CGAL::ON_ORIENTED_BOUNDARY && ss != CGAL::ON_ORIENTED_BOUNDARY) { 605 intersecting_points.push_back(t); 606 const std::set<const Plane*>& planes = bbox_vertex_supporting_planes[td]; 607 CGAL_assertion(planes.size() == 3); 608 intersecting_points_source_planes.push_back(planes); 609 } 610 else { 611 // The intersection is the current edge, nothing to do 612 } 613 } 614 } 615 616 // Decides the orientation of the points 617 if (intersecting_points.size() >= 3) { 618 std::list<Point> pts; 619 for (std::size_t i = 0; i < intersecting_points.size(); ++i) { 620 const Point& p = intersecting_points[i]; 621 const Point2& q = cutting_plane->to_2d(p); 622 pts.push_back(Point(q.x(), q.y(), FT(i))); // the z component stores the point index 623 } 624 625 typedef CGAL::Projection_traits_xy_3<Kernel> Projection; 626 std::list<Point> hull; 627 CGAL::convex_hull_2(pts.begin(), pts.end(), std::back_inserter(hull), Projection()); 628 629 std::vector<Point> ch; 630 std::vector< std::set<const Plane*> > ch_source_planes; 631 for (typename std::list<Point>::iterator it = hull.begin(); it != hull.end(); ++it) { 632 std::size_t idx = std::size_t(it->z()); 633 ch.push_back(intersecting_points[idx]); 634 ch_source_planes.push_back(intersecting_points_source_planes[idx]); 635 } 636 637 if (ch.size() >= 3) { 638 std::vector<Vertex_descriptor> descriptors; 639 for (std::size_t j = 0; j < ch.size(); ++j) { 640 Vertex_descriptor vd = candidate_faces.add_vertex(ch[j]); 641 descriptors.push_back(vd); 642 vertex_supporting_planes[vd] = ch_source_planes[j]; 643 CGAL_assertion(vertex_supporting_planes[vd].size() == 3); 644 } 645 646 Face_descriptor fd = candidate_faces.add_face(descriptors); 647 face_supporting_segments[fd] = g; 648 face_supporting_planes[fd] = cutting_plane; 649 650 // Assigns each edge the supporting planes 651 CGAL::Halfedge_around_face_circulator<Polygon_mesh> hbegin(candidate_faces.halfedge(fd), candidate_faces), done(hbegin); 652 do { 653 Halfedge_descriptor hd = *hbegin; 654 Edge_descriptor ed = candidate_faces.edge(hd); 655 656 Vertex_descriptor s_vd = candidate_faces.source(hd); 657 Vertex_descriptor t_vd = candidate_faces.target(hd); 658 const std::set<const Plane*>& s_planes = vertex_supporting_planes[s_vd]; 659 const std::set<const Plane*>& t_planes = vertex_supporting_planes[t_vd]; 660 std::set<const Plane*> common_planes; 661 std::set_intersection(s_planes.begin(), s_planes.end(), t_planes.begin(), t_planes.end(), std::inserter(common_planes, common_planes.begin())); 662 if (common_planes.size() == 2) { 663 CGAL_assertion(edge_supporting_planes[ed].size() == 0); 664 edge_supporting_planes[ed] = common_planes; 665 CGAL_assertion(edge_supporting_planes[ed].size() == 2); 666 } 667 else // If reached here, there must be topological errors. 668 std::cerr << "topological error" << std::endl; 669 670 ++hbegin; 671 } while (hbegin != done); 672 } 673 } 674 } 675 676 CGAL_assertion(candidate_faces.is_valid()); 677 } 678 679 680 template <typename Kernel> compute_triplet_intersections()681 void Hypothesis<Kernel>::compute_triplet_intersections() { 682 triplet_intersections_.clear(); 683 if (supporting_planes_.size() < 4) // no closed surface will be constructed from fewer than 4 planes 684 return; 685 686 for (std::size_t i = 0; i < supporting_planes_.size(); ++i) { 687 const Plane* plane1 = supporting_planes_[i]; 688 for (std::size_t j = i + 1; j < supporting_planes_.size(); ++j) { 689 const Plane* plane2 = supporting_planes_[j]; 690 for (std::size_t k = j + 1; k < supporting_planes_.size(); ++k) { 691 const Plane* plane3 = supporting_planes_[k]; 692 CGAL_assertion(plane1 < plane2 && plane2 < plane3); 693 Point p; 694 if (internal::intersect_plane_triplet<Plane, Point>(plane1, plane2, plane3, p)) { 695 // Stores the intersection for future query 696 Point* new_point = new Point(p); 697 triplet_intersections_[plane1][plane2][plane3] = new_point; 698 intersecting_points_.push_back(new_point); 699 } 700 } 701 } 702 } 703 } 704 705 706 template <typename Kernel> 707 const typename Hypothesis<Kernel>::Point* query_intersection(const Plane * min_plane,const Plane * mid_plane,const Plane * max_plane)708 Hypothesis<Kernel>::query_intersection(const Plane* min_plane, const Plane* mid_plane, const Plane* max_plane) { 709 CGAL_assertion(min_plane < mid_plane); 710 CGAL_assertion(mid_plane < max_plane); 711 712 if (triplet_intersections_.find(min_plane) == triplet_intersections_.end()) 713 return nullptr; 714 715 Two_planes_to_point_map& map2 = triplet_intersections_[min_plane]; 716 if (map2.find(mid_plane) == map2.end()) 717 return nullptr; 718 719 Plane_to_point_map& map1 = map2[mid_plane]; 720 if (map1.find(max_plane) == map1.end()) 721 return nullptr; 722 723 return map1[max_plane]; 724 } 725 726 727 template <typename Kernel> 728 typename Hypothesis<Kernel>::Halfedge_descriptor split_edge(Polygon_mesh & mesh,const EdgePos & ep,const Plane * cutting_plane)729 Hypothesis<Kernel>::split_edge(Polygon_mesh& mesh, const EdgePos& ep, const Plane* cutting_plane) { 730 // The supporting planes of each edge 731 typename Polygon_mesh::template Property_map<Edge_descriptor, std::set<const Plane*> > edge_supporting_planes = 732 mesh.template property_map<Edge_descriptor, std::set<const Plane*> >("e:supp_plane").first; 733 734 // The supporting planes of each vertex 735 typename Polygon_mesh::template Property_map<Vertex_descriptor, std::set<const Plane*> > vertex_supporting_planes 736 = mesh.template property_map<Vertex_descriptor, std::set<const Plane*> >("v:supp_plane").first; 737 738 // We cannot use const reference, because it will become invalid after splitting 739 std::set<const Plane*> sfs = edge_supporting_planes[ep.edge]; 740 CGAL_assertion(sfs.size() == 2); 741 742 Halfedge_descriptor h = Euler::split_edge(mesh.halfedge(ep.edge), mesh); 743 if (h == Polygon_mesh::null_halfedge()) // failed splitting edge 744 return h; 745 746 Vertex_descriptor v = mesh.target(h); 747 if (v == Polygon_mesh::null_vertex()) // failed splitting edge 748 return Polygon_mesh::null_halfedge(); 749 750 typename Polygon_mesh::template Property_map<Vertex_descriptor, Point>& coords = mesh.points(); 751 coords[v] = *ep.pos; 752 753 Edge_descriptor e1 = mesh.edge(h); 754 edge_supporting_planes[e1] = sfs; 755 Edge_descriptor e2 = mesh.edge(mesh.next(h)); 756 edge_supporting_planes[e2] = sfs; 757 758 vertex_supporting_planes[v] = sfs; 759 vertex_supporting_planes[v].insert(cutting_plane); 760 CGAL_assertion(vertex_supporting_planes[v].size() == 3); 761 762 return h; 763 } 764 765 766 // Cuts f using the cutter and returns the new faces 767 template <typename Kernel> 768 std::vector<typename Hypothesis<Kernel>::Face_descriptor> cut(Face_descriptor face,const Plane * cutting_plane,Polygon_mesh & mesh)769 Hypothesis<Kernel>::cut(Face_descriptor face, const Plane* cutting_plane, Polygon_mesh& mesh) { 770 std::vector<Face_descriptor> new_faces; 771 772 // The supporting plane of each face 773 typename Polygon_mesh::template Property_map<Face_descriptor, const Plane*> face_supporting_planes = 774 mesh.template property_map<Face_descriptor, const Plane*>("f:supp_plane").first; 775 const Plane* supporting_plane = face_supporting_planes[face]; 776 777 if (supporting_plane == cutting_plane) 778 return new_faces; 779 780 // The supporting planar segment of each face 781 typename Polygon_mesh::template Property_map<Face_descriptor, Planar_segment*> face_supporting_segments = 782 mesh.template property_map<Face_descriptor, Planar_segment*>("f:supp_segment").first; 783 784 // The supporting planes of each edge 785 typename Polygon_mesh::template Property_map<Edge_descriptor, std::set<const Plane*> > edge_supporting_planes = 786 mesh.template property_map<Edge_descriptor, std::set<const Plane*> >("e:supp_plane").first; 787 788 Planar_segment* supporting_segment = face_supporting_segments[face]; 789 790 std::vector<Vertex_descriptor> existing_vts; 791 std::vector<EdgePos> new_vts; 792 compute_intersections(mesh, face, cutting_plane, existing_vts, new_vts); 793 794 // We need to check here because new faces are emerging 795 if (existing_vts.size() + new_vts.size() != 2) 796 return new_faces; 797 798 else if (existing_vts.size() == 2) { 799 // Tests if the two intersecting points are both very close to an existing vertex. 800 // Since we allow snapping, we test if the two intersecting points are the same. 801 if (existing_vts[0] == existing_vts[1]) 802 return new_faces; 803 804 // Tests if an edge already exists, i.e., the plane cuts at this edge 805 if (halfedge_exists(existing_vts[0], existing_vts[1], mesh)) 806 return new_faces; 807 } 808 809 Halfedge_descriptor h0 = Polygon_mesh::null_halfedge(); 810 Halfedge_descriptor h1 = Polygon_mesh::null_halfedge(); 811 812 if (existing_vts.size() == 2) { // cutting_plane cuts the face at two existing vertices (not an edge) 813 h0 = mesh.halfedge(existing_vts[0]); 814 h1 = mesh.halfedge(existing_vts[1]); 815 } 816 else if (existing_vts.size() == 1) { 817 h0 = mesh.halfedge(existing_vts[0]); 818 h1 = split_edge(mesh, new_vts[0], cutting_plane); 819 } 820 else if (new_vts.size() == 2) { 821 h0 = split_edge(mesh, new_vts[0], cutting_plane); 822 h1 = split_edge(mesh, new_vts[1], cutting_plane); 823 } 824 CGAL_assertion(h0 != Polygon_mesh::null_halfedge()); 825 CGAL_assertion(h1 != Polygon_mesh::null_halfedge()); 826 827 // To split the face, `h0` and `h1` must be incident to the same face 828 if (mesh.face(h0) != face) { 829 Halfedge_descriptor end = h0; 830 do { 831 h0 = mesh.opposite(mesh.next(h0)); // h0 = h0->next()->opposite(); 832 if (mesh.face(h0) == face) 833 break; 834 } while (h0 != end); 835 } 836 CGAL_assertion(mesh.face(h0) == face); 837 838 if (mesh.face(h1) != face) { 839 Halfedge_descriptor end = h1; 840 do { 841 h1 = mesh.opposite(mesh.next(h1)); // h1 = h1->next()->opposite(); 842 if (mesh.face(h1) == face) 843 break; 844 } while (h1 != end); 845 } 846 CGAL_assertion(mesh.face(h1) == face); 847 848 Halfedge_descriptor h = Euler::split_face(h0, h1, mesh); 849 if (h == Polygon_mesh::null_halfedge() || mesh.face(h) == Polygon_mesh::null_face()) { 850 std::cerr << "Fatal error. could not split face" << std::endl; 851 return new_faces; 852 } 853 854 Edge_descriptor e = mesh.edge(h); 855 edge_supporting_planes[e].insert(supporting_plane); 856 edge_supporting_planes[e].insert(cutting_plane); 857 CGAL_assertion(edge_supporting_planes[e].size() == 2); 858 859 // Now the two faces 860 Face_descriptor f1 = mesh.face(h); 861 face_supporting_segments[f1] = supporting_segment; 862 face_supporting_planes[f1] = supporting_plane; 863 new_faces.push_back(f1); 864 865 Face_descriptor f2 = mesh.face(mesh.opposite(h)); 866 face_supporting_segments[f2] = supporting_segment; 867 face_supporting_planes[f2] = supporting_plane; 868 new_faces.push_back(f2); 869 870 return new_faces; 871 } 872 873 874 template <typename Kernel> compute_intersections(const Polygon_mesh & mesh,Face_descriptor face,const Plane * cutting_plane,std::vector<Vertex_descriptor> & existing_vts,std::vector<EdgePos> & new_vts)875 void Hypothesis<Kernel>::compute_intersections(const Polygon_mesh& mesh, 876 Face_descriptor face, const Plane* cutting_plane, 877 std::vector<Vertex_descriptor>& existing_vts, 878 std::vector<EdgePos>& new_vts) 879 { 880 existing_vts.clear(); 881 new_vts.clear(); 882 883 // The supporting plane of each face 884 typename Polygon_mesh::template Property_map<Face_descriptor, const Plane*> face_supporting_planes = 885 mesh.template property_map<Face_descriptor, const Plane*>("f:supp_plane").first; 886 const Plane* supporting_plane = face_supporting_planes[face]; 887 if (supporting_plane == cutting_plane) 888 return; 889 890 typename Polygon_mesh::template Property_map<Edge_descriptor, std::set<const Plane*> > edge_supporting_planes 891 = mesh.template property_map<Edge_descriptor, std::set<const Plane*> >("e:supp_plane").first; 892 893 const typename Polygon_mesh::template Property_map<Vertex_descriptor, Point>& coords = mesh.points(); 894 895 Halfedge_descriptor cur = mesh.halfedge(face); 896 Halfedge_descriptor end = cur; 897 do { 898 Edge_descriptor ed = mesh.edge(cur); 899 const std::set<const Plane*>& supporting_planes = edge_supporting_planes[ed]; 900 if (supporting_planes.find(cutting_plane) != supporting_planes.end()) // the edge lies on the cutting plane 901 return; 902 903 Vertex_descriptor s_vd = mesh.source(cur); 904 Vertex_descriptor t_vd = mesh.target(cur); 905 const Point& s = coords[s_vd]; 906 const Point& t = coords[t_vd]; 907 908 CGAL::Oriented_side s_side = cutting_plane->oriented_side(s); 909 CGAL::Oriented_side t_side = cutting_plane->oriented_side(t); 910 911 if (t_side == CGAL::ON_ORIENTED_BOUNDARY) { 912 if (s_side == CGAL::ON_ORIENTED_BOUNDARY) // the edge lies on the cutting plane 913 return; 914 else 915 existing_vts.push_back(t_vd); 916 } 917 918 else if ( 919 (s_side == CGAL::ON_POSITIVE_SIDE && t_side == CGAL::ON_NEGATIVE_SIDE) || 920 (s_side == CGAL::ON_NEGATIVE_SIDE && t_side == CGAL::ON_POSITIVE_SIDE)) // intersects at the interior of the edge 921 { 922 FT s_sdist = CGAL::squared_distance(*cutting_plane, s); 923 FT t_sdist = CGAL::squared_distance(*cutting_plane, t); 924 925 if (s_sdist <= CGAL::snap_squared_distance_threshold<FT>()) // plane cuts at vertex 's' 926 existing_vts.push_back(s_vd); 927 else if (t_sdist <= CGAL::snap_squared_distance_threshold<FT>()) // plane cuts at vertex 't' 928 existing_vts.push_back(t_vd); 929 else { 930 const Plane* plane1 = *(supporting_planes.begin()); 931 const Plane* plane2 = *(supporting_planes.rbegin()); 932 const Plane* plane3 = const_cast<const Plane*>(cutting_plane); 933 934 if (plane3 != plane1 && plane3 != plane2) { 935 internal::sort_increasing(plane1, plane2, plane3); 936 const Point* p = query_intersection(plane1, plane2, plane3); 937 if (p) { 938 if (CGAL::squared_distance(*p, s) <= CGAL::snap_squared_distance_threshold<FT>()) // snap to 's' 939 existing_vts.push_back(s_vd); 940 else if (CGAL::squared_distance(*p, t) <= CGAL::snap_squared_distance_threshold<FT>()) // snap to 't' 941 existing_vts.push_back(t_vd); 942 else 943 new_vts.push_back(EdgePos(ed, p)); 944 } 945 else 946 std::cerr << "Fatal error: should have intersection" << std::endl; 947 } 948 else 949 std::cerr << "Fatal error: should not have duplicated planes." << std::endl; 950 } 951 } 952 953 else { 954 // Nothing needs to do here, we will test the next edge 955 } 956 957 cur = mesh.next(cur); 958 } while (cur != end); 959 } 960 961 962 template <typename Kernel> halfedge_exists(Vertex_descriptor v1,Vertex_descriptor v2,const Polygon_mesh & mesh)963 bool Hypothesis<Kernel>::halfedge_exists(Vertex_descriptor v1, Vertex_descriptor v2, const Polygon_mesh& mesh) { 964 Halfedge_descriptor h = mesh.halfedge(v1); 965 Halfedge_descriptor end = h; 966 do { 967 Halfedge_descriptor opp = mesh.opposite(h); 968 if (mesh.target(opp) == v2) 969 return true; 970 h = mesh.prev(opp); 971 } while (h != end); 972 return false; 973 } 974 975 976 // Tests if face 'f' insects plane 'plane' 977 template <typename Kernel> do_intersect(const Polygon_mesh & mesh,Face_descriptor f,const Plane * plane)978 bool Hypothesis<Kernel>::do_intersect(const Polygon_mesh& mesh, Face_descriptor f, const Plane* plane) { 979 std::vector<Vertex_descriptor> existing_vts; 980 std::vector<EdgePos> new_vts; 981 compute_intersections(mesh, f, plane, existing_vts, new_vts); 982 983 if (existing_vts.size() == 2) { 984 if (!halfedge_exists(existing_vts[0], existing_vts[1], mesh)) 985 return true; 986 } 987 else if (existing_vts.size() + new_vts.size() == 2) 988 return true; 989 990 return false; 991 } 992 993 994 // Collects all faces in 'mesh' that intersect 'face', stored in std::set() so easier to erase 995 template <typename Kernel> collect_intersecting_faces(Face_descriptor face,const Polygon_mesh & mesh)996 std::set<typename Hypothesis<Kernel>::Face_descriptor> Hypothesis<Kernel>::collect_intersecting_faces(Face_descriptor face, const Polygon_mesh& mesh) 997 { 998 // The supporting plane of each face 999 typename Polygon_mesh::template Property_map<Face_descriptor, const Plane*> face_supporting_planes = 1000 mesh.template property_map<Face_descriptor, const Plane*>("f:supp_plane").first; 1001 1002 // The supporting planar segment of each face 1003 typename Polygon_mesh::template Property_map<Face_descriptor, Planar_segment*> face_supporting_segments = 1004 mesh.template property_map<Face_descriptor, Planar_segment*>("f:supp_segment").first; 1005 1006 std::set<Face_descriptor> intersecting_faces; 1007 for(auto f : mesh.faces()) { 1008 if (f == face || 1009 face_supporting_segments[f] == face_supporting_segments[face] || 1010 face_supporting_planes[f] == face_supporting_planes[face]) 1011 { 1012 continue; 1013 } 1014 1015 const Plane* plane = face_supporting_planes[face]; 1016 CGAL_assertion(plane != nullptr); 1017 if (do_intersect(mesh, f, plane)) 1018 intersecting_faces.insert(f); 1019 } 1020 return intersecting_faces; 1021 } 1022 1023 1024 template <typename Kernel> pairwise_intersection(Polygon_mesh & candidate_faces)1025 void Hypothesis<Kernel>::pairwise_intersection(Polygon_mesh& candidate_faces) 1026 { 1027 // Pre-computes all potential intersection of plane triplets 1028 compute_triplet_intersections(); 1029 1030 // Since we are going to split faces, we can not use the reference. 1031 // const Polygon_mesh::Face_range& all_faces = mesh.faces(); 1032 // So we make a local copy 1033 std::vector<Face_descriptor> all_faces(candidate_faces.faces().begin(), candidate_faces.faces().end()); 1034 1035 // The supporting plane of each face 1036 typename Polygon_mesh::template Property_map<Face_descriptor, const Plane*> face_supporting_planes 1037 = candidate_faces.template property_map<Face_descriptor, const Plane*>("f:supp_plane").first; 1038 1039 for (std::size_t i = 0; i < all_faces.size(); ++i) { 1040 Face_descriptor face = all_faces[i]; 1041 const Plane* face_plane = face_supporting_planes[face]; 1042 1043 std::set<Face_descriptor> intersecting_faces = collect_intersecting_faces(face, candidate_faces); 1044 if (intersecting_faces.empty()) 1045 continue; 1046 1047 std::vector<Face_descriptor> cutting_faces(intersecting_faces.begin(), intersecting_faces.end()); 1048 1049 // 1. 'face' will be cut by all the intersecting faces 1050 // \note After each cut, the original face doesn't exist any more and it is replaced by multiple pieces. 1051 // then each piece will be cut by another face. 1052 std::vector<Face_descriptor> faces_to_be_cut; 1053 faces_to_be_cut.push_back(face); 1054 while (!intersecting_faces.empty()) { 1055 Face_descriptor cutting_face = *(intersecting_faces.begin()); 1056 const Plane* cutting_plane = face_supporting_planes[cutting_face]; 1057 1058 std::set<Face_descriptor> new_faces; // stores the new faces 1059 std::set<Face_descriptor> remained_faces; // faces that will be cut later 1060 for (std::size_t j = 0; j < faces_to_be_cut.size(); ++j) { 1061 Face_descriptor current_face = faces_to_be_cut[j]; 1062 std::vector<Face_descriptor> tmp = cut(current_face, cutting_plane, candidate_faces); 1063 new_faces.insert(tmp.begin(), tmp.end()); 1064 if (tmp.empty()) { // no actual cut occurred. The face will be cut later 1065 remained_faces.insert(current_face); 1066 } 1067 } 1068 1069 // The new faces might be cut by other faces 1070 faces_to_be_cut = std::vector<Face_descriptor>(new_faces.begin(), new_faces.end()); 1071 1072 // Don't forget the remained faces 1073 faces_to_be_cut.insert(faces_to_be_cut.end(), remained_faces.begin(), remained_faces.end()); 1074 1075 // The job of cutting_face is done, remove it 1076 intersecting_faces.erase(cutting_face); 1077 } 1078 1079 // 2. All the cutting_faces will be cut by f. 1080 for (std::size_t j = 0; j < cutting_faces.size(); ++j) { 1081 cut(cutting_faces[j], face_plane, candidate_faces); 1082 } 1083 } 1084 1085 CGAL_assertion(candidate_faces.is_valid()); 1086 } 1087 1088 1089 template <class Kernel> extract_adjacency(const Polygon_mesh & candidate_faces)1090 typename Hypothesis<Kernel>::Adjacency Hypothesis<Kernel>::extract_adjacency(const Polygon_mesh& candidate_faces) 1091 { 1092 typename Polygon_mesh::template Property_map<Vertex_descriptor, std::set<const Plane*> > vertex_supporting_planes 1093 = candidate_faces.template property_map<Vertex_descriptor, std::set<const Plane*> >("v:supp_plane").first; 1094 1095 // An edge is denoted by its two end points 1096 typedef typename std::unordered_map<const Point*, std::set<Halfedge_descriptor> > Edge_map; 1097 typedef typename std::unordered_map<const Point*, Edge_map > Face_pool; 1098 Face_pool face_pool; 1099 1100 for(auto h : candidate_faces.halfedges()) { 1101 Face_descriptor f = candidate_faces.face(h); 1102 if (f == Polygon_mesh::null_face()) 1103 continue; 1104 1105 Vertex_descriptor sd = candidate_faces.source(h); 1106 Vertex_descriptor td = candidate_faces.target(h); 1107 const std::set<const Plane*>& set_s = vertex_supporting_planes[sd]; 1108 const std::set<const Plane*>& set_t = vertex_supporting_planes[td]; 1109 CGAL_assertion(set_s.size() == 3); 1110 CGAL_assertion(set_t.size() == 3); 1111 1112 std::vector<const Plane*> s_planes(set_s.begin(), set_s.end()); 1113 CGAL_assertion(s_planes[0] < s_planes[1]); 1114 CGAL_assertion(s_planes[1] < s_planes[2]); 1115 const Point* s = triplet_intersections_[s_planes[0]][s_planes[1]][s_planes[2]]; 1116 1117 std::vector<const Plane*> t_planes(set_t.begin(), set_t.end()); 1118 CGAL_assertion(t_planes[0] < t_planes[1]); 1119 CGAL_assertion(t_planes[1] < t_planes[2]); 1120 const Point* t = triplet_intersections_[t_planes[0]][t_planes[1]][t_planes[2]]; 1121 1122 if (s > t) 1123 std::swap(s, t); 1124 face_pool[s][t].insert(candidate_faces.halfedge(f)); 1125 } 1126 1127 Adjacency fans; 1128 typename Face_pool::const_iterator it = face_pool.begin(); 1129 for (; it != face_pool.end(); ++it) { 1130 const Point* s = it->first; 1131 const Edge_map& tmp = it->second; 1132 typename Edge_map::const_iterator cur = tmp.begin(); 1133 for (; cur != tmp.end(); ++cur) { 1134 const Point* t = cur->first; 1135 const std::set<Halfedge_descriptor>& faces = cur->second; 1136 Intersection fan; 1137 fan.s = s; 1138 fan.t = t; 1139 fan.insert(fan.end(), faces.begin(), faces.end()); 1140 fans.push_back(fan); 1141 } 1142 } 1143 1144 return fans; 1145 } 1146 1147 } //namespace internal 1148 1149 } //namespace CGAL 1150 1151 1152 #endif // CGAL_POLYGONAL_SURFACE_RECONSTRUCTION_HYPOTHESIS_H 1153