1 // Copyright (c) 2014 INRIA Sophia-Antipolis (France), INRIA Lorraine LORIA. 2 // All rights reserved. 3 // 4 // This file is part of CGAL (www.cgal.org). 5 // 6 // $URL: https://github.com/CGAL/cgal/blob/v5.3/Optimal_transportation_reconstruction_2/include/CGAL/Optimal_transportation_reconstruction_2.h $ 7 // $Id: Optimal_transportation_reconstruction_2.h 0779373 2020-03-26T13:31:46+01:00 Sébastien Loriot 8 // SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial 9 // 10 // Author(s) : Fernando de Goes, Pierre Alliez, Ivo Vigan, Clément Jamin 11 12 #ifndef CGAL_OPTIMAL_TRANSPORTATION_RECONSTRUCTION_2_H_ 13 #define CGAL_OPTIMAL_TRANSPORTATION_RECONSTRUCTION_2_H_ 14 15 #include <CGAL/license/Optimal_transportation_reconstruction_2.h> 16 17 18 #include <CGAL/OTR_2/Reconstruction_triangulation_2.h> 19 #include <CGAL/OTR_2/Reconstruction_edge_2.h> 20 21 #include <CGAL/property_map.h> 22 #include <CGAL/Real_timer.h> 23 #include <CGAL/Random.h> 24 25 #include <iterator> 26 #include <iostream> 27 #include <vector> 28 #include <list> 29 #include <algorithm> 30 #include <utility> // std::pair 31 32 #include <boost/multi_index_container.hpp> 33 #include <boost/multi_index/mem_fun.hpp> 34 #include <boost/multi_index/ordered_index.hpp> 35 #include <boost/multi_index/identity.hpp> 36 #include <CGAL/boost/iterator/transform_iterator.hpp> 37 #include <boost/type_traits/is_float.hpp> 38 39 namespace CGAL { 40 41 42 /*! 43 \ingroup PkgOptimalTransportationReconstruction2Classes 44 45 This class provides a means to reconstruct a 1-dimensional shape from a set of 2D points with masses. 46 The algorithm computes an initial 2D Delaunay triangulation from the input points, 47 and performs a simplification of the triangulation by performing half edge collapses, edge flips and vertex relocations. 48 49 The edges are either processed in the order imposed by an priority queue, or 50 in an order based on random selection of edge collapse operators. 51 As the exhaustive priority queue guarantees a higher quality it is the default. 52 The user can switch to the other method, for example for an initial 53 simplification round, by calling `set_random_sample_size()`. 54 55 By default edge flip operators are applied to ensure that every edge of the 56 triangulation are candidate to be collapsed, while preserving a valid embedding 57 of the triangulation. This option can be disabled by calling 58 \link set_use_flip() `set_use_flip(false)`\endlink to reduce the running times. 59 60 By default the vertices are not relocated after each half edge collapse. 61 This option can be changed by setting the number of vertex relocation steps 62 performed between two edge collapse operators. 63 64 The simplification is performed by calling either 65 \link run_until() `run_until(n)`\endlink or \link run() `run(steps)`\endlink. 66 The former simplifies the triangulation until n points remain, while the latter 67 stops after `steps` edge collapse operators have been performed. 68 Furthermore, we can relocate the vertices by calling `relocate_all_points()`. 69 70 \tparam Traits a model of the concept `OptimalTransportationReconstructionTraits_2`. 71 72 \tparam PointPMap a model of `ReadablePropertyMap` with value type `Traits::Point_2`. 73 Defaults to <a href="https://www.boost.org/doc/libs/release/libs/property_map/doc/identity_property_map.html">`boost::typed_identity_property_map<Traits::Point_2>`</a> 74 (for the case the input is points without mass). 75 76 \tparam MassPMap a model of `ReadablePropertyMap` with value type `Traits::FT` 77 Defaults to <a href="https://www.boost.org/doc/libs/release/libs/property_map/doc/static_property_map.html">`boost::static_property_map<Traits::FT>`</a> 78 (for the case the input is points without mass). 79 80 */ 81 template< 82 class Traits, 83 class PointPMap = boost::typed_identity_property_map <typename Traits::Point_2>, 84 class MassPMap = boost::static_property_map <typename Traits::FT> > 85 class Optimal_transportation_reconstruction_2 86 { 87 public: 88 89 /// \name Types 90 /// @{ 91 /*! 92 Number type. 93 */ 94 typedef typename Traits::FT FT; 95 96 /*! 97 Point type. 98 */ 99 typedef typename Traits::Point_2 Point; 100 101 /*! 102 Segment type. 103 */ 104 typedef typename Traits::Segment_2 Segment; 105 106 /// \cond SKIP_IN_MANUAL 107 /*! 108 Vector type. 109 */ 110 typedef typename Traits::Vector_2 Vector; 111 112 typedef typename std::pair<Point, FT> PointMassPair; 113 typedef typename std::vector<PointMassPair> PointMassList; 114 115 116 /*! 117 The Output simplex. 118 */ 119 typedef OTR_2::Reconstruction_triangulation_2<Traits> Triangulation; 120 121 typedef typename Triangulation::Vertex Vertex; 122 typedef typename Triangulation::Vertex_handle Vertex_handle; 123 typedef typename Triangulation::Vertex_iterator Vertex_iterator; 124 typedef typename Triangulation::Vertex_circulator Vertex_circulator; 125 typedef typename Triangulation::Finite_vertices_iterator 126 Finite_vertices_iterator; 127 128 typedef typename Triangulation::Edge Edge; 129 typedef typename Triangulation::Edge_circulator Edge_circulator; 130 typedef typename Triangulation::Finite_edges_iterator Finite_edges_iterator; 131 132 typedef typename Triangulation::Face_handle Face_handle; 133 typedef typename Triangulation::Face_circulator Face_circulator; 134 typedef typename Triangulation::Finite_faces_iterator Finite_faces_iterator; 135 136 typedef typename Triangulation::Vertex_handle_map Vertex_handle_map; 137 typedef typename Triangulation::Face_handle_map Face_handle_map; 138 139 typedef typename Triangulation::Vertex_handle_set Vertex_handle_set; 140 typedef typename Triangulation::Edge_set Edge_set; 141 142 typedef typename Triangulation::Edge_vector Edge_vector; 143 typedef std::list<Edge> Edge_list; 144 145 typedef typename Triangulation::Cost_ Cost_; 146 typedef typename Triangulation::Sample_ Sample_; 147 typedef typename Triangulation::Sample_vector Sample_vector; 148 typedef typename Triangulation::Sample_vector_const_iterator 149 Sample_vector_const_iterator; 150 151 typedef typename Triangulation::PSample PSample; 152 typedef typename Triangulation::SQueue SQueue; 153 154 typedef typename Triangulation::Rec_edge_2 Rec_edge_2; 155 156 typedef typename Triangulation::MultiIndex MultiIndex; 157 158 /// @} 159 160 protected: 161 Triangulation m_dt; 162 Traits const& m_traits; 163 MultiIndex m_mindex; 164 int m_ignore; 165 int m_verbose; 166 std::size_t m_mchoice; // # Edges 167 bool m_use_flip; 168 FT m_alpha; // [0, 1] 169 FT m_ghost; // ghost vs solid 170 unsigned int m_relocation; // # relocations 171 FT m_tolerance; 172 173 PointPMap point_pmap; 174 MassPMap mass_pmap; 175 176 /// \endcond 177 178 public: 179 180 /// \name Initialization 181 /// @{ 182 183 /*! 184 Constructor of the optimal transportation reconstruction class. 185 It builds an initial simplicial complex 186 for a given range of point-mass pairs. 187 188 \tparam InputRange is a model of `Range` with forward iterators, 189 providing input points and point masses through the 190 `PointPMap` and `MassPMap` property maps. 191 192 \param input_range Range of input data. 193 \param point_map A `ReadablePropertyMap` used to access the input points. 194 \param mass_map A `ReadablePropertyMap` used to access the input 195 points' masses. 196 \param sample_size If `sample_size != 0`, the size of the random sample 197 which replaces the exhaustive priority queue. 198 \param use_flip If `true` the edge flipping procedure is used to ensure 199 that every edge can be made collapsible. 200 \param relocation The number of point relocations that are performed 201 between two edge collapses. 202 \param verbose Controls how much console output is produced by 203 the algorithm. The values are 0, 1, or > 1. 204 \param traits The traits class. 205 */ 206 template <class InputRange> 207 Optimal_transportation_reconstruction_2( 208 const InputRange& input_range, 209 PointPMap point_map = PointPMap(), 210 MassPMap mass_map = MassPMap(1), 211 std::size_t sample_size = 0, 212 bool use_flip = true, 213 unsigned int relocation = 2, 214 int verbose = 0, 215 Traits traits = Traits()) m_dt(traits)216 : m_dt(traits), 217 m_traits(m_dt.geom_traits()), 218 m_ignore(0), 219 m_verbose(verbose), 220 m_mchoice(sample_size), 221 m_use_flip(use_flip), 222 m_alpha(0.5), 223 m_ghost(1.0), 224 m_relocation(relocation), 225 m_tolerance (FT(-1.)), 226 point_pmap(point_map), 227 mass_pmap(mass_map) 228 { 229 initialize(input_range.begin(), input_range.end()); 230 } 231 232 /// @} 233 234 /// \name Settting Parameters 235 /// @{ 236 /*! 237 If `sample_size == 0`, the simplification is performed using an exhaustive priority queue. 238 If `sample_size` is stricly positive the simplification is performed using a 239 multiple choice approach, ie, a best-choice selection in a random sample of 240 edge collapse operators, of size `sample_size`. A typical value for the sample 241 size is 15, but this value must be enlarged when targeting a very coarse simplification. 242 \param sample_size If `sample_size != 0`, the size of the random sample replaces the priority queue. 243 */ set_random_sample_size(std::size_t sample_size)244 void set_random_sample_size(std::size_t sample_size) { 245 m_mchoice = sample_size; 246 } 247 248 /*! 249 Determines how much console output the algorithm generates. 250 If set to a value larger than 0 251 details about the reconstruction process are written to `std::cerr`. 252 253 \param verbose The verbosity level. 254 */ set_verbose(int verbose)255 void set_verbose(int verbose) { 256 m_verbose = verbose; 257 } 258 259 260 261 262 /*! 263 The use_flip parameter determines whether the edge flipping procedure 264 is used for the half-edge collapse. 265 */ set_use_flip(const bool use_flip)266 void set_use_flip(const bool use_flip) { 267 m_use_flip = use_flip; 268 } 269 270 271 /*! 272 Sets the number of vertex relocations 273 that are performed between two edge collapses. 274 */ set_relocation(unsigned int relocation)275 void set_relocation(unsigned int relocation) { 276 m_relocation = relocation; 277 } 278 279 /// \cond SKIP_IN_MANUAL relocation()280 unsigned int relocation() const { 281 return m_relocation; 282 } 283 /// \endcond 284 285 286 /*! 287 \param relevance The relevance threshold used for filtering the edges. 288 An edge is relevant from the approximation point of view 289 if it is long, covers a large mass (or equivalently the 290 number of points when all masses are equal), and has a 291 small transport cost. This notion is defined as 292 \f$ m(e) * |e|^2 / cost(e) \f$, where \f$ m(e) \f$ 293 denotes the mass of the points approximated by the edge, 294 \f$ |e| \f$ denotes the edge length and \f$ cost(e) \f$ 295 its approximation error. 296 As the cost is defined by mass time squared distance the 297 relevance is unitless. 298 299 The default value is 1, so that all edges receiving some mass 300 are considered relevant. 301 Setting a large relevance value is used to get robustness to a 302 large amount of outliers. 303 */ set_relevance(const FT relevance)304 void set_relevance(const FT relevance) { 305 m_ghost = relevance; 306 m_dt.ghost_factor() = m_ghost; 307 } 308 309 310 /// \cond SKIP_IN_MANUAL ghost()311 FT ghost() { 312 return m_ghost; 313 } 314 tolerance()315 FT tolerance() const { return m_tolerance; } 316 317 /// @} 318 319 /// \cond SKIP_IN_MANUAL 320 Optimal_transportation_reconstruction_2()321 Optimal_transportation_reconstruction_2() 322 : m_traits(m_dt.geom_traits()) 323 { 324 initialize_parameters(); 325 } 326 327 ~Optimal_transportation_reconstruction_2()328 ~Optimal_transportation_reconstruction_2() { 329 clear(); 330 } 331 initialize_parameters()332 void initialize_parameters() { 333 m_verbose = 0; 334 m_mchoice = 0; 335 m_use_flip = true; 336 m_alpha = FT(0.5); 337 m_ghost = FT(1); 338 m_relocation = 0; 339 340 m_ignore = 0; 341 } 342 343 //Function if one wants to create a Optimal_transportation_reconstruction_2 344 //without yet specifying the input in the constructor. 345 template <class InputIterator> initialize(InputIterator start_itr,InputIterator beyond_itr,PointPMap point_map,MassPMap mass_map)346 void initialize( 347 InputIterator start_itr, 348 InputIterator beyond_itr, 349 PointPMap point_map, 350 MassPMap mass_map) 351 { 352 point_pmap = point_map; 353 mass_pmap = mass_map; 354 355 initialize(start_itr, beyond_itr); 356 } 357 358 359 template <class InputIterator> initialize(InputIterator start,InputIterator beyond)360 void initialize(InputIterator start, InputIterator beyond) { 361 362 clear(); 363 Property_map_to_unary_function<PointPMap> get_point(point_pmap); 364 365 Bbox_2 bbox = bbox_2( 366 boost::make_transform_iterator(start,get_point), 367 boost::make_transform_iterator(beyond,get_point)); 368 369 insert_loose_bbox(bbox); 370 init(start, beyond); 371 372 std::vector<Sample_*> m_samples; 373 for (InputIterator it = start; it != beyond; it++) { 374 Point point = get(point_pmap, *it); 375 FT mass = get( mass_pmap, *it); 376 Sample_* s = new Sample_(point, mass); 377 m_samples.push_back(s); 378 } 379 assign_samples(m_samples.begin(), m_samples.end()); 380 } 381 382 template <class InputIterator> initialize_with_custom_vertices(InputIterator samples_start,InputIterator samples_beyond,InputIterator vertices_start,InputIterator vertices_beyond,PointPMap point_map,MassPMap mass_map)383 void initialize_with_custom_vertices(InputIterator samples_start, 384 InputIterator samples_beyond, 385 InputIterator vertices_start, 386 InputIterator vertices_beyond, 387 PointPMap point_map, 388 MassPMap mass_map) { 389 point_pmap = point_map; 390 mass_pmap = mass_map; 391 clear(); 392 Property_map_to_unary_function<PointPMap> get_point(point_pmap); 393 394 Bbox_2 bbox = bbox_2( 395 boost::make_transform_iterator(samples_start,get_point), 396 boost::make_transform_iterator(samples_beyond,get_point)); 397 398 insert_loose_bbox(bbox); 399 init(vertices_start, vertices_beyond); 400 401 std::vector<Sample_*> m_samples; 402 for (InputIterator it = samples_start; it != samples_beyond; it++) { 403 #ifdef CGAL_USE_PROPERTY_MAPS_API_V1 404 Point point = get(point_pmap, it); 405 FT mass = get( mass_pmap, it); 406 #else 407 Point point = get(point_pmap, *it); 408 FT mass = get( mass_pmap, *it); 409 #endif 410 Sample_* s = new Sample_(point, mass); 411 m_samples.push_back(s); 412 } 413 assign_samples(m_samples.begin(), m_samples.end()); 414 } 415 416 417 template <class Vector> random_vec(const FT scale)418 Vector random_vec(const FT scale) const 419 { 420 FT dx = -scale + get_default_random().get_double() * 2* scale; 421 FT dy = -scale + get_default_random().get_double() * 2* scale; 422 return m_traits.construct_vector_2_object()(dx, dy); 423 } 424 clear()425 void clear() { 426 Sample_vector samples; 427 m_dt.collect_all_samples(samples); 428 // Deallocate samples 429 for (Sample_vector_const_iterator s_it = samples.begin(); 430 s_it != samples.end(); ++s_it) 431 { 432 delete *s_it; 433 } 434 } 435 436 437 // INIT // insert_loose_bbox(const Bbox_2 & bbox)438 void insert_loose_bbox(const Bbox_2& bbox) { 439 CGAL::Real_timer timer; 440 if (m_verbose > 0) 441 std::cerr << "insert loose bbox..."; 442 443 double dl = (std::max)((bbox.xmax()-bbox.xmin()) / 2., 444 (bbox.ymax()-bbox.ymin()) / 2.); 445 446 timer.start(); 447 int nb = static_cast<int>(m_dt.number_of_vertices()); 448 typename Traits::Construct_point_2 point_2 449 = m_traits.construct_point_2_object(); 450 insert_point(point_2(bbox.xmin()-dl, bbox.ymin()-dl), true, nb++); 451 insert_point(point_2(bbox.xmin()-dl, bbox.ymax()+dl), true, nb++); 452 insert_point(point_2(bbox.xmax()+dl, bbox.ymax()+dl), true, nb++); 453 insert_point(point_2(bbox.xmax()+dl, bbox.ymin()-dl), true, nb++); 454 455 if (m_verbose > 0) 456 std::cerr << "done (" << nb << " vertices, " 457 << timer.time() << " s)" << std::endl; 458 } 459 460 template<class Iterator> // value_type = Point* init(Iterator begin,Iterator beyond)461 void init(Iterator begin, Iterator beyond) { 462 CGAL::Real_timer timer; 463 if (m_verbose > 0) 464 std::cerr << "init..."; 465 466 timer.start(); 467 int nb = static_cast<int>(m_dt.number_of_vertices()); 468 m_dt.infinite_vertex()->pinned() = true; 469 for (Iterator it = begin; it != beyond; it++) { 470 Point point = get(point_pmap, *it); 471 insert_point(point, false, nb++); 472 } 473 474 if (m_verbose > 0) 475 std::cerr << "done (" << nb << " vertices, " 476 << timer.time() << " s)" 477 << std::endl; 478 } 479 480 private: insert_point(const Point & point,const bool pinned,const int id)481 Vertex_handle insert_point( 482 const Point& point, const bool pinned, const int id) 483 { 484 Vertex_handle v = m_dt.insert(point); 485 v->pinned() = pinned; 486 v->id() = id; 487 return v; 488 } 489 public: 490 491 // ASSIGNMENT // 492 cleanup_assignments()493 void cleanup_assignments() { 494 m_dt.cleanup_assignments(); 495 } 496 497 template<class Iterator> // value_type = Sample_* assign_samples(Iterator begin,Iterator end)498 void assign_samples(Iterator begin, Iterator end) { 499 CGAL::Real_timer timer; 500 if (m_verbose > 0) 501 std::cerr << "assign samples..."; 502 503 timer.start(); 504 m_dt.assign_samples(begin, end); 505 m_dt.reset_all_costs(); 506 507 if (m_verbose > 0) 508 std::cerr << "done (" << timer.time() << " s)" << std::endl; 509 } 510 reassign_samples()511 void reassign_samples() { 512 Sample_vector samples; 513 m_dt.collect_all_samples(samples); 514 m_dt.cleanup_assignments(); 515 m_dt.assign_samples(samples.begin(), samples.end()); 516 m_dt.reset_all_costs(); 517 } 518 reassign_samples_around_vertex(Vertex_handle vertex)519 void reassign_samples_around_vertex(Vertex_handle vertex) { 520 Sample_vector samples; 521 m_dt.collect_samples_from_vertex(vertex, samples, true); 522 m_dt.assign_samples(samples.begin(), samples.end()); 523 524 Edge_vector hull; 525 m_dt.get_edges_from_star_minus_link(vertex, hull, true); 526 update_cost(hull.begin(), hull.end()); 527 } 528 529 do_collapse(Edge edge)530 bool do_collapse(Edge edge) { 531 Vertex_handle s = m_dt.source_vertex(edge); 532 Vertex_handle t = m_dt.target_vertex(edge); 533 534 if (m_verbose > 1) { 535 std::cerr << std::endl << "do collapse (" 536 << s->id() << "->" << t->id() << ") ... " << std::endl; 537 } 538 539 Sample_vector samples; 540 m_dt.collect_samples_from_vertex(s, samples, true); 541 542 Edge_vector hull; 543 m_dt.get_edges_from_star_minus_link(s, hull, true); 544 545 if (m_mchoice == 0) 546 remove_stencil_from_pqueue(hull.begin(), hull.end()); 547 548 if (m_use_flip) 549 m_dt.make_collapsible(edge, hull.begin(), hull.end(), m_verbose); 550 551 // debug test 552 bool ok = m_dt.check_kernel_test(edge); 553 if (!ok) { 554 if (m_verbose > 1) 555 std::cerr << "do_collapse: kernel test failed: " << std::endl; 556 return false; 557 } 558 // 559 560 m_dt.collapse(edge, m_verbose); 561 562 m_dt.assign_samples(samples.begin(), samples.end()); 563 564 update_cost(hull.begin(), hull.end()); 565 566 if (m_mchoice == 0) 567 push_stencil_to_pqueue(hull.begin(), hull.end()); 568 569 for (unsigned int i = 0; i < m_relocation; ++i) { 570 relocate_one_ring(hull.begin(), hull.end()); 571 } 572 573 if (m_verbose > 1) { 574 std::cerr << "done" << std::endl; 575 } 576 577 return true; 578 } 579 simulate_collapse(const Edge & edge,Cost_ & cost)580 bool simulate_collapse(const Edge& edge, Cost_& cost) { 581 bool ok; 582 Vertex_handle s = m_dt.source_vertex(edge); 583 Vertex_handle t = m_dt.target_vertex(edge); 584 585 if (m_verbose > 1) { 586 std::cerr << "simulate collapse (" 587 << s->id() << "->" << t->id() << ") ... " << std::endl; 588 } 589 590 Triangulation copy; 591 Edge copy_edge = copy_star(edge, copy); 592 Vertex_handle copy_source = copy.source_vertex(copy_edge); 593 594 if (m_use_flip) { 595 Edge_vector copy_hull; 596 copy.get_edges_from_star_minus_link(copy_source, copy_hull, true); 597 ok = copy.make_collapsible(copy_edge, copy_hull.begin(), 598 copy_hull.end(), m_verbose); 599 if (!ok) { 600 if (m_verbose > 1) 601 std::cerr << "simulation: failed (make collapsible)" << std::endl; 602 return false; 603 } 604 ok = copy.check_validity_test(); 605 if (!ok) { 606 if (m_verbose > 1) 607 std::cerr << "simulation: failed (validity test)" << std::endl; 608 return false; 609 } 610 } 611 612 ok = copy.check_kernel_test(copy_edge); 613 if (!ok) { 614 if (m_verbose > 1) 615 std::cerr << "simulation: failed (kernel test)" << std::endl; 616 return false; 617 } 618 619 copy.collapse(copy_edge, m_verbose); 620 621 ok = copy.check_validity_test(); 622 if (!ok) { 623 if (m_verbose > 1) 624 std::cerr << "simulation: failed (validity test)" << std::endl; 625 return false; 626 } 627 628 Sample_vector samples; 629 m_dt.collect_samples_from_vertex(s, samples, false); 630 631 backup_samples(samples.begin(), samples.end()); 632 copy.assign_samples_brute_force(samples.begin(), samples.end()); 633 copy.reset_all_costs(); 634 cost = copy.compute_total_cost(); 635 cost.set_total_weight (samples); 636 restore_samples(samples.begin(), samples.end()); 637 638 if (m_verbose > 1) { 639 std::cerr << "done" << std::endl; 640 } 641 642 return true; 643 } 644 645 template<class Iterator> // value_type = Sample_* backup_samples(Iterator begin,Iterator end)646 void backup_samples(Iterator begin, Iterator end) const { 647 for (Iterator it = begin; it != end; ++it) { 648 Sample_* sample = *it; 649 sample->backup(); 650 } 651 } 652 653 template<class Iterator> // value_type = Sample_* restore_samples(Iterator begin,Iterator end)654 void restore_samples(Iterator begin, Iterator end) const { 655 for (Iterator it = begin; it != end; ++it) { 656 Sample_* sample = *it; 657 sample->restore(); 658 } 659 } 660 661 // PEDGE // 662 decimate()663 bool decimate() { 664 bool ok; 665 Rec_edge_2 pedge; 666 ok = pick_edge(m_mchoice, pedge); 667 if (!ok) 668 return false; 669 670 ok = do_collapse(pedge.edge()); 671 if (!ok) 672 return false; 673 return true; 674 } 675 is_above_tolerance(const Rec_edge_2 & pedge)676 bool is_above_tolerance (const Rec_edge_2& pedge) 677 { 678 if (m_tolerance == (FT)(-1.)) 679 return false; 680 FT cost = CGAL::approximate_sqrt (FT(pedge.after() / pedge.total_weight())); 681 return cost > m_tolerance; 682 } 683 create_pedge(const Edge & edge,Rec_edge_2 & pedge)684 bool create_pedge(const Edge& edge, Rec_edge_2& pedge) { 685 Cost_ after_cost; 686 bool ok = simulate_collapse(edge, after_cost); 687 if (!ok) 688 return false; 689 690 Vertex_handle source = m_dt.source_vertex(edge); 691 Cost_ before_cost = m_dt.compute_cost_around_vertex(source); 692 693 FT before = before_cost.finalize(m_alpha); 694 FT after = after_cost.finalize(m_alpha); 695 pedge = Rec_edge_2(edge, before, after, after_cost.total_weight()); 696 697 if (is_above_tolerance (pedge)) 698 return false; 699 700 return true; 701 } 702 703 704 // COST // 705 init_cost()706 void init_cost() { 707 m_dt.reset_all_costs(); 708 } 709 710 template<class Iterator> // value_type = Edge update_cost(Iterator begin,Iterator end)711 void update_cost(Iterator begin, Iterator end) { 712 Edge_vector edges; 713 collect_cost_stencil(m_dt, begin, end, edges); 714 715 typename Edge_vector::iterator ei; 716 for (ei = edges.begin(); ei != edges.end(); ++ei) { 717 Edge edge = *ei; 718 m_dt.update_cost(edge); 719 } 720 } 721 722 template<class Iterator> // value_type = Edge collect_cost_stencil(const Triangulation & mesh,Iterator begin,Iterator end,Edge_vector & edges)723 void collect_cost_stencil( 724 const Triangulation& mesh, Iterator begin, Iterator end, 725 Edge_vector& edges) const 726 { 727 Edge_set done; 728 Edge_list fifo; 729 for (Iterator it = begin; it != end; ++it) { 730 Edge edge = *it; 731 fifo.push_back(edge); 732 done.insert(edge); 733 } 734 735 while (!fifo.empty()) { 736 Edge edge = fifo.front(); 737 fifo.pop_front(); 738 739 edge = mesh.twin_edge(edge); 740 edges.push_back(edge); 741 742 Edge next = mesh.next_edge(edge); 743 if (done.insert(next).second) 744 fifo.push_back(next); 745 746 Edge prev = mesh.prev_edge(edge); 747 if (done.insert(prev).second) 748 fifo.push_back(prev); 749 } 750 } 751 752 // PQUEUE (MCHOICE or EXHAUSTIVE) // 753 pick_edge(std::size_t nb,Rec_edge_2 & best_pedge)754 bool pick_edge(std::size_t nb, Rec_edge_2& best_pedge) { 755 if (m_dt.number_of_faces() < 2) 756 return false; 757 758 std::size_t ne = 2 * m_dt.tds().number_of_edges(); 759 if (nb > ne) 760 nb = ne; 761 762 bool ok; 763 if (nb == 0) { 764 ok = pick_edge_from_pqueue(best_pedge); 765 return ok; 766 } 767 m_mindex.clear(); 768 769 if (nb == ne) { 770 ok = pick_edge_brute_force(best_pedge); 771 return ok; 772 } 773 774 ok = pick_edge_randomly(nb, best_pedge); 775 return ok; 776 } 777 pick_edge_from_pqueue(Rec_edge_2 & best_pedge)778 bool pick_edge_from_pqueue(Rec_edge_2& best_pedge) { 779 if (m_mindex.empty()) 780 populate_pqueue(); 781 if (m_mindex.empty()) 782 return false; 783 best_pedge = *(m_mindex.template get<1>()).begin(); 784 (m_mindex.template get<0>()).erase(best_pedge); 785 return true; 786 } 787 pick_edge_brute_force(Rec_edge_2 & best_pedge)788 bool pick_edge_brute_force(Rec_edge_2& best_pedge) { 789 MultiIndex mindex; 790 Finite_edges_iterator ei; 791 for (ei = m_dt.finite_edges_begin(); ei != m_dt.finite_edges_end(); 792 ++ei) { 793 Edge edge = *ei; 794 push_to_mindex(edge, mindex); 795 796 edge = m_dt.twin_edge(edge); 797 push_to_mindex(edge, mindex); 798 } 799 if (mindex.empty()) 800 return false; 801 best_pedge = *(mindex.template get<1>()).begin(); 802 return true; 803 } 804 pick_edge_randomly(std::size_t nb,Rec_edge_2 & best_pedge)805 bool pick_edge_randomly(std::size_t nb, Rec_edge_2& best_pedge) { 806 MultiIndex mindex; 807 for (std::size_t i = 0; i < nb; ++i) { 808 Rec_edge_2 pedge; 809 if (random_pedge(pedge)) 810 mindex.insert(pedge); 811 } 812 if (mindex.empty()) 813 return false; 814 best_pedge = *(mindex.template get<1>()).begin(); 815 return true; 816 } 817 populate_pqueue()818 void populate_pqueue() { 819 Finite_edges_iterator ei; 820 for (ei = m_dt.finite_edges_begin(); ei != m_dt.finite_edges_end(); 821 ++ei) { 822 Edge edge = *ei; 823 push_to_mindex(edge, m_mindex); 824 825 edge = m_dt.twin_edge(edge); 826 push_to_mindex(edge, m_mindex); 827 } 828 } 829 830 push_to_mindex(const Edge & edge,MultiIndex & mindex)831 bool push_to_mindex(const Edge& edge, MultiIndex& mindex) { 832 if (m_dt.is_pinned(edge)) 833 return false; 834 if (m_dt.is_target_cyclic(edge)) 835 return false; 836 837 Rec_edge_2 pedge; 838 bool ok = create_pedge(edge, pedge); 839 if (!ok) 840 return false; 841 mindex.insert(pedge); 842 return true; 843 } 844 845 846 random_pedge(Rec_edge_2 & pedge)847 bool random_pedge(Rec_edge_2& pedge) { 848 for (unsigned int i = 0; i < 10; ++i) { 849 Edge edge = m_dt.random_finite_edge(); 850 if (m_dt.is_pinned(edge)) 851 continue; 852 if (m_dt.is_target_cyclic(edge)) 853 continue; 854 bool ok = create_pedge(edge, pedge); 855 if (ok) 856 return true; 857 } 858 return false; 859 } 860 861 template<class Iterator> // value_type = Edge remove_stencil_from_pqueue(Iterator begin,Iterator end)862 void remove_stencil_from_pqueue(Iterator begin, Iterator end) 863 { 864 if (m_mindex.empty()) 865 return; 866 867 Edge_vector edges; 868 collect_pqueue_stencil(m_dt, begin, end, edges); 869 870 typename Edge_vector::const_iterator ei; 871 for (ei = edges.begin(); ei != edges.end(); ++ei) { 872 Edge edge = *ei; 873 (m_mindex.template get<0>()).erase(Rec_edge_2(edge)); 874 } 875 } 876 877 template<class Iterator> // value_type = Edge push_stencil_to_pqueue(Iterator begin,Iterator end)878 void push_stencil_to_pqueue(Iterator begin, Iterator end) { 879 Edge_vector edges; 880 collect_pqueue_stencil(m_dt, begin, end, edges); 881 882 typename Edge_vector::const_iterator ei; 883 for (ei = edges.begin(); ei != edges.end(); ++ei) { 884 Edge edge = *ei; 885 push_to_mindex(edge, m_mindex); 886 } 887 } 888 889 template<class Iterator> // value_type = Edge collect_pqueue_stencil(const Triangulation & mesh,Iterator begin,Iterator end,Edge_vector & edges)890 void collect_pqueue_stencil( 891 const Triangulation& mesh, Iterator begin, Iterator end, 892 Edge_vector& edges) const 893 { 894 Vertex_handle_set vertex_set; 895 for (Iterator it = begin; it != end; ++it) { 896 Edge edge = *it; 897 Edge twin = mesh.twin_edge(edge); 898 899 Vertex_handle s = mesh.source_vertex(edge); 900 if (!s->pinned()) 901 vertex_set.insert(s); 902 903 Vertex_handle t = mesh.target_vertex(edge); 904 if (!t->pinned()) 905 vertex_set.insert(t); 906 907 Vertex_handle f = mesh.opposite_vertex(edge); 908 if (!f->pinned()) 909 vertex_set.insert(f); 910 911 Vertex_handle b = mesh.opposite_vertex(twin); 912 if (!b->pinned()) 913 vertex_set.insert(b); 914 } 915 916 typename Vertex_handle_set::const_iterator vi; 917 for (vi = vertex_set.begin(); vi != vertex_set.end(); ++vi) { 918 Vertex_handle v = *vi; 919 Edge_circulator ecirc = mesh.incident_edges(v); 920 Edge_circulator eend = ecirc; 921 CGAL_For_all(ecirc, eend) 922 { 923 Edge edge = *ecirc; 924 if (mesh.source_vertex(edge) != v) 925 edge = mesh.twin_edge(edge); 926 edges.push_back(edge); 927 } 928 } 929 } 930 931 // COPY STAR // 932 933 // edge must not be pinned or have cyclic target copy_star(const Edge & edge,Triangulation & copy)934 Edge copy_star(const Edge& edge, Triangulation& copy) { 935 copy.tds().clear(); 936 Vertex_handle vinf = copy.tds().create_vertex(); 937 copy.set_infinite_vertex (vinf); 938 copy.tds().set_dimension(2); 939 copy.infinite_vertex()->pinned() = true; 940 941 // copy vertices 942 Vertex_handle_map cvmap; 943 944 Vertex_handle s = m_dt.source_vertex(edge); 945 CGAL_assertion(s != m_dt.infinite_vertex() ); 946 947 Vertex_handle cs = copy.tds().create_vertex(); 948 cvmap[s] = copy_vertex(s, cs); 949 950 Vertex_circulator vcirc = m_dt.incident_vertices(s); 951 Vertex_circulator vend = vcirc; 952 CGAL_For_all(vcirc, vend) 953 { 954 Vertex_handle v = vcirc; 955 CGAL_assertion(v!=m_dt.infinite_vertex()); 956 CGAL_assertion (cvmap.find(v) == cvmap.end()); 957 Vertex_handle cv = copy.tds().create_vertex(); 958 cvmap[v] = copy_vertex(v, cv); 959 } 960 961 // copy faces 962 Face_handle_map cfmap; 963 Face_circulator fcirc = m_dt.incident_faces(s); 964 Face_circulator fend = fcirc; 965 CGAL_For_all(fcirc, fend) 966 { 967 Face_handle f = fcirc; 968 Face_handle cf = copy.tds().create_face(); 969 cfmap[f] = copy_face(f, cf, cvmap); 970 } 971 972 // set neighbors 973 fcirc = m_dt.incident_faces(s); 974 fend = fcirc; 975 CGAL_For_all(fcirc, fend) 976 { 977 Face_handle f = fcirc; 978 copy_neighbors(f, s, cfmap); 979 } 980 981 // make copy homeomorphic to S^2 982 close_copy_mesh(cs, copy); 983 984 // copy samples surrounding star 985 copy_samples(s, cs, cfmap, copy); 986 987 // get copy of edge 988 Edge copy_edge = get_copy_edge(edge, cvmap, cfmap); 989 return copy_edge; 990 } 991 copy_vertex(Vertex_handle v0,Vertex_handle v1)992 Vertex_handle copy_vertex(Vertex_handle v0, Vertex_handle v1) const { 993 v1->id() = v0->id(); 994 v1->set_point(v0->point()); 995 v1->pinned() = v0->pinned(); 996 v1->set_sample(v0->sample()); 997 return v1; 998 } 999 copy_face(Face_handle f0,Face_handle f1,Vertex_handle_map & vmap)1000 Face_handle copy_face( 1001 Face_handle f0, Face_handle f1, Vertex_handle_map& vmap) const 1002 { 1003 for (unsigned int i = 0; i < 3; ++i) { 1004 Vertex_handle v0i = f0->vertex(i); 1005 CGAL_assertion (vmap.find(v0i) != vmap.end()); 1006 Vertex_handle v1i = vmap[v0i]; 1007 CGAL_assertion (v1i != Vertex_handle()); 1008 f1->set_vertex(i, v1i); 1009 v1i->set_face(f1); 1010 } 1011 return f1; 1012 } 1013 copy_neighbors(Face_handle f,Vertex_handle v,Face_handle_map & fmap)1014 void copy_neighbors( 1015 Face_handle f, Vertex_handle v, 1016 Face_handle_map& fmap) const 1017 { 1018 int i = f->index(v); 1019 Face_handle cf = fmap[f]; 1020 1021 if (fmap.find(f->neighbor(i)) != fmap.end()) { 1022 Face_handle fi = f->neighbor(i); 1023 Face_handle cfi = fmap[fi]; 1024 cf->set_neighbor(i, cfi); 1025 } 1026 1027 for (unsigned int j = 0; j < 2; ++j) { 1028 i = (i + 1) % 3; 1029 Face_handle fi = f->neighbor(i); 1030 Face_handle cfi = fmap[fi]; 1031 cf->set_neighbor(i, cfi); 1032 } 1033 } 1034 close_copy_mesh(Vertex_handle vertex,Triangulation & copy)1035 void close_copy_mesh(Vertex_handle vertex, Triangulation& copy) const { 1036 std::vector<Face_handle> outer_faces; 1037 1038 Face_circulator fcirc = copy.incident_faces(vertex); 1039 Face_circulator fend = fcirc; 1040 CGAL_For_all(fcirc, fend) 1041 { 1042 Face_handle face = fcirc; 1043 int i = face->index(vertex); 1044 1045 if (face->neighbor(i) != Face_handle()) 1046 continue; 1047 1048 Vertex_handle v1 = face->vertex((i + 1) % 3); 1049 Vertex_handle v2 = face->vertex((i + 2) % 3); 1050 1051 Face_handle outer = copy.tds().create_face(); 1052 outer->set_vertex(0, copy.infinite_vertex()); 1053 outer->set_vertex(1, v2); 1054 outer->set_vertex(2, v1); 1055 1056 face->set_neighbor(i, outer); 1057 outer->set_neighbor(0, face); 1058 1059 outer_faces.push_back(outer); 1060 } 1061 1062 for (unsigned int i = 0; i < outer_faces.size(); ++i) { 1063 unsigned int j = (i + 1) % outer_faces.size(); 1064 outer_faces[i]->set_neighbor(2, outer_faces[j]); 1065 outer_faces[j]->set_neighbor(1, outer_faces[i]); 1066 } 1067 1068 if (!outer_faces.empty()) 1069 copy.infinite_vertex()->set_face(outer_faces[0]); 1070 } 1071 copy_samples(Vertex_handle vertex,Vertex_handle copy_vertex,Face_handle_map & fmap,Triangulation & copy)1072 void copy_samples( 1073 Vertex_handle vertex, Vertex_handle copy_vertex, 1074 Face_handle_map& fmap, Triangulation& copy) const 1075 { 1076 Face_circulator fcirc = m_dt.incident_faces(vertex); 1077 Face_circulator fend = fcirc; 1078 CGAL_For_all(fcirc, fend) 1079 { 1080 Face_handle face = fcirc; 1081 int index = face->index(vertex); 1082 Edge twin = m_dt.twin_edge(Edge(face, index)); 1083 1084 Face_handle copy_face = fmap[face]; 1085 index = copy_face->index(copy_vertex); 1086 Edge copy_twin = copy.twin_edge(Edge(copy_face, index)); 1087 1088 Sample_vector samples; 1089 m_dt.collect_samples_from_edge(twin, samples); 1090 copy_twin.first->samples(copy_twin.second) = samples; 1091 } 1092 copy_vertex->set_sample(nullptr); 1093 } 1094 get_copy_edge(const Edge & edge,Vertex_handle_map & vmap,Face_handle_map & fmap)1095 Edge get_copy_edge( 1096 const Edge& edge, Vertex_handle_map& vmap, Face_handle_map& fmap) const 1097 { 1098 Face_handle f = edge.first; 1099 Vertex_handle v = f->vertex(edge.second); 1100 1101 Face_handle cf = fmap[f]; 1102 Vertex_handle cv = vmap[v]; 1103 1104 return Edge(cf, cf->index(cv)); 1105 } 1106 1107 // RELOCATION // 1108 relocate_one_vertex(Vertex_handle vertex)1109 void relocate_one_vertex(Vertex_handle vertex) { 1110 std::swap(vertex->point(), vertex->relocated()); 1111 reassign_samples_around_vertex(vertex); 1112 } 1113 1114 template<class Iterator> // value_type = Edge relocate_one_ring(Iterator begin,Iterator end)1115 void relocate_one_ring(Iterator begin, Iterator end) { 1116 Vertex_handle_set vertices; 1117 for (Iterator it = begin; it != end; ++it) { 1118 Edge edge = *it; 1119 vertices.insert(m_dt.source_vertex(edge)); 1120 vertices.insert(m_dt.target_vertex(edge)); 1121 } 1122 1123 typename Vertex_handle_set::const_iterator vi; 1124 for (vi = vertices.begin(); vi != vertices.end(); ++vi) { 1125 Vertex_handle v = *vi; 1126 if (v->pinned()) 1127 continue; 1128 v->relocated() = compute_relocation(v); 1129 } 1130 1131 for (vi = vertices.begin(); vi != vertices.end(); ++vi) { 1132 Vertex_handle v = *vi; 1133 if (v->pinned()) 1134 continue; 1135 if (v->point() == v->relocated()) 1136 continue; 1137 1138 Edge_vector hull; 1139 m_dt.get_edges_from_star_minus_link(v, hull, false); 1140 bool ok = m_dt.is_in_kernel(v->relocated(), hull.begin(), 1141 hull.end()); 1142 1143 if (ok) { 1144 // do relocation 1145 FT norm_bef = m_dt.compute_cost_around_vertex(v).norm(); 1146 relocate_one_vertex(v); 1147 FT norm_aft = m_dt.compute_cost_around_vertex(v).norm(); 1148 1149 if (norm_bef < norm_aft) { 1150 // undo relocation 1151 relocate_one_vertex(v); 1152 } else if (m_mchoice == 0) { 1153 // update queue 1154 hull.clear(); 1155 m_dt.get_edges_from_star_minus_link(v, hull, true); 1156 remove_stencil_from_pqueue(hull.begin(), hull.end()); 1157 push_stencil_to_pqueue(hull.begin(), hull.end()); 1158 } 1159 } 1160 } 1161 } 1162 1163 /// \endcond 1164 1165 1166 /// \cond SKIP_IN_MANUAL compute_gradient(Vertex_handle vertex)1167 Vector compute_gradient(Vertex_handle vertex) const { 1168 Vector grad = m_traits.construct_vector_2_object()(FT(0), FT(0)); 1169 Edge_circulator ecirc = m_dt.incident_edges(vertex); 1170 Edge_circulator eend = ecirc; 1171 CGAL_For_all(ecirc, eend) 1172 { 1173 Edge edge = *ecirc; 1174 if (m_dt.source_vertex(edge) != vertex) 1175 edge = m_dt.twin_edge(edge); 1176 1177 if (m_dt.get_plan(edge) == 0) 1178 grad = m_traits.construct_sum_of_vectors_2_object()( 1179 grad, compute_gradient_for_plan0(edge)); 1180 else 1181 grad = m_traits.construct_sum_of_vectors_2_object()( 1182 grad, compute_gradient_for_plan1(edge)); 1183 } 1184 return grad; 1185 } 1186 1187 // If the underlying number type used is not a floating point base 1188 // number type (like a multiprecision), the coordinates of the points 1189 // will increase a lot due to the relocation step. These functions 1190 // simply turn a relocated point to a rounded to double version. relocate_on_the_double_grid(Point &,boost::true_type)1191 void relocate_on_the_double_grid(Point&, boost::true_type) const 1192 {} relocate_on_the_double_grid(Point & p,boost::false_type)1193 void relocate_on_the_double_grid(Point& p, boost::false_type) const 1194 { 1195 double x=to_double(m_traits.compute_x_2_object()(p)); 1196 double y=to_double(m_traits.compute_y_2_object()(p)); 1197 p=m_traits.construct_point_2_object()(FT(x),FT(y)); 1198 } relocate_on_the_double_grid(Point & p)1199 void relocate_on_the_double_grid(Point& p) const 1200 { 1201 relocate_on_the_double_grid(p, 1202 typename boost::is_float<typename Traits::FT>::type()); 1203 } 1204 compute_relocation(Vertex_handle vertex)1205 Point compute_relocation(Vertex_handle vertex) const { 1206 FT coef = FT(0); 1207 Vector rhs = m_traits.construct_vector_2_object()(FT(0), FT(0)); 1208 1209 Edge_circulator ecirc = m_dt.incident_edges(vertex); 1210 Edge_circulator eend = ecirc; 1211 CGAL_For_all(ecirc, eend) 1212 { 1213 Edge edge = *ecirc; 1214 if (m_dt.source_vertex(edge) != vertex) 1215 edge = m_dt.twin_edge(edge); 1216 1217 if (m_dt.get_plan(edge) == 0) 1218 compute_relocation_for_plan0(edge, coef, rhs); 1219 else 1220 compute_relocation_for_plan1(edge, coef, rhs); 1221 } 1222 compute_relocation_for_vertex(vertex, coef, rhs); 1223 1224 if (coef == FT(0)) 1225 return vertex->point(); 1226 1227 Point res = m_traits.construct_translated_point_2_object()( 1228 CGAL::ORIGIN, 1229 m_traits.construct_scaled_vector_2_object()(rhs, FT(1) / coef)); 1230 relocate_on_the_double_grid(res); 1231 return res; 1232 } 1233 compute_relocation_for_vertex(Vertex_handle vertex,FT & coef,Vector & rhs)1234 void compute_relocation_for_vertex( 1235 Vertex_handle vertex, FT& coef, Vector& rhs) const 1236 { 1237 Sample_* sample = vertex->sample(); 1238 if (sample) { 1239 const FT m = sample->mass(); 1240 const Point& ps = sample->point(); 1241 rhs = m_traits.construct_sum_of_vectors_2_object()(rhs, 1242 m_traits.construct_scaled_vector_2_object()( 1243 m_traits.construct_vector_2_object()(CGAL::ORIGIN, ps), m)); 1244 coef += m; 1245 } 1246 } 1247 compute_gradient_for_plan0(const Edge & edge)1248 Vector compute_gradient_for_plan0(const Edge& edge) const { 1249 Edge twin = m_dt.twin_edge(edge); 1250 const Point& pa = m_dt.source_vertex(edge)->point(); 1251 const Point& pb = m_dt.target_vertex(edge)->point(); 1252 1253 Sample_vector samples; 1254 m_dt.collect_samples_from_edge(edge, samples); 1255 m_dt.collect_samples_from_edge(twin, samples); 1256 1257 Vector grad = m_traits.construct_vector_2_object()(FT(0), FT(0)); 1258 Sample_vector_const_iterator it; 1259 for (it = samples.begin(); it != samples.end(); ++it) { 1260 Sample_* sample = *it; 1261 const FT m = sample->mass(); 1262 const Point& ps = sample->point(); 1263 1264 FT Da = m_traits.compute_squared_distance_2_object()(ps, pa); 1265 FT Db = m_traits.compute_squared_distance_2_object()(ps, pb); 1266 if (Da < Db) 1267 grad = m_traits.construct_sum_of_vectors_2_object()( 1268 grad, 1269 m_traits.construct_scaled_vector_2_object()( 1270 m_traits.construct_vector_2_object()(ps, pa), m)); 1271 } 1272 return grad; 1273 } 1274 compute_relocation_for_plan0(const Edge & edge,FT & coef,Vector & rhs)1275 void compute_relocation_for_plan0( 1276 const Edge& edge, FT& coef, Vector& rhs) const 1277 { 1278 Edge twin = m_dt.twin_edge(edge); 1279 const Point& pa = m_dt.source_vertex(edge)->point(); 1280 const Point& pb = m_dt.target_vertex(edge)->point(); 1281 1282 Sample_vector samples; 1283 m_dt.collect_samples_from_edge(edge, samples); 1284 m_dt.collect_samples_from_edge(twin, samples); 1285 1286 Sample_vector_const_iterator it; 1287 for (it = samples.begin(); it != samples.end(); ++it) { 1288 Sample_* sample = *it; 1289 const FT m = sample->mass(); 1290 const Point& ps = sample->point(); 1291 1292 FT Da = m_traits.compute_squared_distance_2_object()(ps, pa); 1293 FT Db = m_traits.compute_squared_distance_2_object()(ps, pb); 1294 1295 if (Da < Db) { 1296 rhs = m_traits.construct_sum_of_vectors_2_object()(rhs, 1297 m_traits.construct_scaled_vector_2_object()( 1298 m_traits.construct_vector_2_object()(CGAL::ORIGIN, ps), m)); 1299 coef += m; 1300 } 1301 } 1302 } 1303 compute_gradient_for_plan1(const Edge & edge)1304 Vector compute_gradient_for_plan1(const Edge& edge) const { 1305 //FT M = m_dt.get_mass(edge); 1306 const Point& pa = m_dt.source_vertex(edge)->point(); 1307 const Point& pb = m_dt.target_vertex(edge)->point(); 1308 1309 SQueue queue; 1310 m_dt.sort_samples_from_edge(edge, queue); 1311 1312 //FT start = FT(0); 1313 Vector grad = m_traits.construct_vector_2_object()(FT(0), FT(0)); 1314 while (!queue.empty()) { 1315 PSample psample = queue.top(); 1316 queue.pop(); 1317 1318 const FT m = psample.sample()->mass(); 1319 const Point& ps = psample.sample()->point(); 1320 1321 // normal + tangnetial 1322 const FT coord = psample.priority(); 1323 Point pf = m_traits.construct_translated_point_2_object()( 1324 CGAL::ORIGIN, 1325 m_traits.construct_sum_of_vectors_2_object()( 1326 m_traits.construct_scaled_vector_2_object()( 1327 m_traits.construct_vector_2_object()(CGAL::ORIGIN, pa), 1328 1.0 - coord), 1329 m_traits.construct_scaled_vector_2_object()( 1330 m_traits.construct_vector_2_object()(CGAL::ORIGIN, pb), 1331 coord))); 1332 grad = m_traits.construct_sum_of_vectors_2_object()( 1333 grad, 1334 m_traits.construct_scaled_vector_2_object()( 1335 m_traits.construct_vector_2_object()(ps, pf), m * (1.0 - coord))); 1336 1337 /* 1338 // only normal 1339 FT bin = m/M; 1340 FT center = start + 0.5*bin; 1341 Point pc = CGAL::ORIGIN + (1.0-center)*(pa - CGAL::ORIGIN) + center*(pb - CGAL::ORIGIN); 1342 start += bin; 1343 grad = grad + m*(bin*bin/12.0)*(pa - pb); 1344 grad = grad + m*(1.0-center)*(pc - pf); 1345 */ 1346 } 1347 return grad; 1348 } 1349 compute_relocation_for_plan1(const Edge & edge,FT & coef,Vector & rhs)1350 void compute_relocation_for_plan1( 1351 const Edge& edge, FT& coef, Vector& rhs) const 1352 { 1353 //FT M = m_dt.get_mass(edge); 1354 const Point& pb = m_dt.target_vertex(edge)->point(); 1355 1356 SQueue queue; 1357 m_dt.sort_samples_from_edge(edge, queue); 1358 1359 //FT start = FT(0); 1360 while (!queue.empty()) { 1361 PSample psample = queue.top(); 1362 queue.pop(); 1363 1364 const FT m = psample.sample()->mass(); 1365 const Point& ps = psample.sample()->point(); 1366 1367 const FT coord = psample.priority(); 1368 const FT one_minus_coord = 1.0 - coord; 1369 1370 // normal + tangential 1371 coef += m * one_minus_coord * one_minus_coord; 1372 rhs = m_traits.construct_sum_of_vectors_2_object()( 1373 rhs, 1374 m_traits.construct_scaled_vector_2_object()( 1375 m_traits.construct_sum_of_vectors_2_object()( 1376 m_traits.construct_vector_2_object()(CGAL::ORIGIN, ps), 1377 m_traits.construct_scaled_vector_2_object()( 1378 m_traits.construct_vector_2_object()(CGAL::ORIGIN, pb), -coord)), 1379 m * one_minus_coord)); 1380 1381 /* 1382 // only normal 1383 FT bin = m/M; 1384 FT center = start + 0.5*bin; 1385 Point pc = CGAL::ORIGIN + (1.0-center)*(pa - CGAL::ORIGIN) + center*(pb - CGAL::ORIGIN); 1386 start += bin; 1387 grad = grad + m*(bin*bin/12.0)*(pa - pb); 1388 grad = grad + m*(1.0-center)*(pc - pf); 1389 */ 1390 1391 /* 1392 // only normal 1393 FT bin = m/M; 1394 FT center = start + 0.5*bin; 1395 FT one_minus_center = 1.0 - center; 1396 start += bin; 1397 1398 coef += m*bin*bin/12.0; 1399 rhs = rhs + m*(bin*bin/12.0)*(pb - CGAL::ORIGIN); 1400 1401 coef += m*one_minus_center*(coord - center); 1402 rhs = rhs + m*one_minus_center*(coord - center)*(pb - CGAL::ORIGIN); 1403 */ 1404 } 1405 } 1406 print_stats_debug()1407 void print_stats_debug() const 1408 { 1409 int nb_solid = 0; 1410 int nb_ghost = 0; 1411 for (Finite_edges_iterator ei = m_dt.finite_edges_begin(); 1412 ei != m_dt.finite_edges_end(); ++ei) 1413 { 1414 Edge edge = *ei; 1415 if (m_dt.is_ghost(edge)) 1416 nb_ghost++; 1417 else 1418 nb_solid++; 1419 } 1420 1421 std::cerr << "STATS" << std::endl; 1422 std::cerr << "# vertices : " << m_dt.number_of_vertices()-4 << std::endl; 1423 std::cerr << "# isolated vertices : " << number_of_isolated_vertices() << std::endl; 1424 std::cerr << "# triangles: " << m_dt.number_of_faces() << std::endl; 1425 std::cerr << "# edges: " << m_dt.tds().number_of_edges() << std::endl; 1426 std::cerr << "# solid: " << nb_solid << std::endl; 1427 std::cerr << "# ghost: " << nb_ghost << std::endl; 1428 } 1429 1430 1431 /*! 1432 Returns the number of vertices present in the reconstructed triangulation. 1433 */ number_of_vertices()1434 std::size_t number_of_vertices() const { 1435 return m_dt.number_of_vertices() - 4; 1436 1437 } 1438 1439 /*! 1440 Returns the number of isolated vertices present in the reconstructed triangulation. 1441 */ number_of_isolated_vertices()1442 int number_of_isolated_vertices () const 1443 { 1444 int nb_isolated = 0; 1445 for (Vertex_iterator vi = m_dt.vertices_begin(); 1446 vi != m_dt.vertices_end(); ++vi) 1447 { 1448 if (!((*vi).has_sample_assigned())) 1449 continue; 1450 1451 typename Triangulation::Edge_circulator start = m_dt.incident_edges(vi); 1452 typename Triangulation::Edge_circulator cur = start; 1453 1454 do { 1455 if (!m_dt.is_ghost(*cur)) { 1456 ++nb_isolated; 1457 break; 1458 } 1459 ++cur; 1460 } while (cur != start); 1461 } 1462 return nb_isolated; 1463 } 1464 1465 /*! 1466 Returns the number of (solid) edges present in the reconstructed triangulation. 1467 */ number_of_edges()1468 int number_of_edges() const { 1469 int nb_solid = 0; 1470 for (Finite_edges_iterator ei = m_dt.finite_edges_begin(); 1471 ei != m_dt.finite_edges_end(); ++ei) 1472 { 1473 Edge edge = *ei; 1474 if (m_dt.is_ghost(edge)) 1475 continue; 1476 nb_solid++; 1477 } 1478 return nb_solid; 1479 } 1480 1481 1482 /*! 1483 Returns the cost of the (solid) edges present in the 1484 reconstructed triangulation. 1485 */ total_edge_cost()1486 FT total_edge_cost() const { 1487 FT total_cost = 0; 1488 for (Finite_edges_iterator ei = m_dt.finite_edges_begin(); 1489 ei != m_dt.finite_edges_end(); ++ei) 1490 { 1491 Edge edge = *ei; 1492 if (m_dt.is_ghost(edge)) 1493 continue; 1494 1495 total_cost += m_dt.get_cost(edge).finalize(); 1496 } 1497 return total_cost; 1498 } 1499 1500 /// \endcond 1501 1502 1503 /// \name Simplification 1504 /// You can freely mix calls of the following functions. 1505 /// @{ 1506 /*! 1507 Computes a shape consisting of `np` points, reconstructing the input 1508 points. 1509 \param np The number of points which will be present in the output. 1510 \return `true` if the number of points `np` was reached, `false` 1511 if the algorithm was prematurely ended because no more edge 1512 collapse was possible. 1513 */ run_until(std::size_t np)1514 bool run_until(std::size_t np) { 1515 m_tolerance = (FT)(-1.); 1516 CGAL::Real_timer timer; 1517 if (m_verbose > 0) 1518 std::cerr << "reconstruct until " << np << " V"; 1519 1520 timer.start(); 1521 std::size_t N = np + 4; 1522 std::size_t performed = 0; 1523 while (m_dt.number_of_vertices() > N) { 1524 bool ok = decimate(); 1525 if (!ok) 1526 break; 1527 performed++; 1528 } 1529 1530 if (m_verbose) 1531 std::cerr << " done" << " (" << performed 1532 << " iters, " << m_dt.number_of_vertices() - 4 << " V " 1533 << timer.time() << " s)" 1534 << std::endl; 1535 1536 return (m_dt.number_of_vertices() <= N); 1537 } 1538 1539 /*! 1540 Computes a shape, reconstructing the input, by performing `steps` 1541 edge collapse operators on the output simplex. 1542 \param steps The number of edge collapse operators to be performed. 1543 \return `true` if the required number of steps was performed, 1544 `false` if the algorithm was prematurely ended because no more 1545 edge collapse was possible. 1546 */ run(const unsigned int steps)1547 bool run(const unsigned int steps) { 1548 m_tolerance = (FT)(-1.); 1549 CGAL::Real_timer timer; 1550 if (m_verbose > 0) 1551 std::cerr << "reconstruct " << steps; 1552 1553 timer.start(); 1554 unsigned int performed = 0; 1555 for (unsigned int i = 0; i < steps; ++i) { 1556 bool ok = decimate(); 1557 if (!ok) 1558 break; 1559 performed++; 1560 } 1561 1562 if (m_verbose > 0) 1563 std::cerr << " done" << " (" << performed << "/" 1564 << steps << " iters, " << m_dt.number_of_vertices() - 4 1565 << " V, " << timer.time() << " s)" 1566 << std::endl; 1567 return (performed == steps); 1568 } 1569 1570 1571 /*! 1572 Computes a shape, reconstructing the input, by performing edge 1573 collapse operators on the output simplex until the user-defined 1574 tolerance is reached. 1575 1576 \note The tolerance is given in the sense of the Wasserstein 1577 distance. It is _not_ a Hausdorff tolerance: it does not mean that 1578 the distance between the input samples and the output polyline is 1579 guaranteed to be less than `tolerance`. It means that the square 1580 root of transport cost per mass (homogeneous to a distance) is at 1581 most `tolerance`. 1582 1583 \param tolerance Tolerance on the Wasserstein distance. 1584 */ run_under_wasserstein_tolerance(const FT tolerance)1585 void run_under_wasserstein_tolerance (const FT tolerance) { 1586 m_tolerance = tolerance; 1587 CGAL::Real_timer timer; 1588 if (m_verbose > 0) 1589 std::cerr << "reconstruct under tolerance " << tolerance; 1590 1591 timer.start(); 1592 unsigned int performed = 0; 1593 while (decimate ()) 1594 performed++; 1595 1596 if (m_verbose > 0) 1597 std::cerr << " done" << " (" << performed 1598 << " iters, " << m_dt.number_of_vertices() - 4 1599 << " V, " << timer.time() << " s)" 1600 << std::endl; 1601 } 1602 1603 1604 /*! 1605 Since noise and missing data may prevent the reconstructed shape to have sharp corners well located, the algorithm offers the possibility to automatically relocate points after each edge collapse. The new location of the points is chosen such that the fitting of the output segments to the input points is improved. 1606 */ relocate_all_points()1607 void relocate_all_points() { 1608 CGAL::Real_timer timer; 1609 if (m_verbose > 0) 1610 std::cerr << "relocate all points" << "..."; 1611 1612 timer.start(); 1613 m_mindex.clear(); // pqueue must be recomputed 1614 1615 for (Finite_vertices_iterator v = m_dt.finite_vertices_begin(); 1616 v != m_dt.finite_vertices_end(); ++v) { 1617 if (v->pinned()) 1618 continue; 1619 v->relocated() = compute_relocation(v); 1620 } 1621 1622 for (Finite_vertices_iterator v = m_dt.finite_vertices_begin(); 1623 v != m_dt.finite_vertices_end(); ++v) { 1624 if (v->pinned()) 1625 continue; 1626 if (v->point() == v->relocated()) 1627 continue; 1628 1629 Edge_vector hull; 1630 m_dt.get_edges_from_star_minus_link(v, hull, false); 1631 bool ok = m_dt.is_in_kernel(v->relocated(), hull.begin(), 1632 hull.end()); 1633 1634 if (ok) { 1635 // do relocation 1636 FT norm_bef = m_dt.compute_cost_around_vertex(v).norm(); 1637 relocate_one_vertex(v); 1638 FT norm_aft = m_dt.compute_cost_around_vertex(v).norm(); 1639 1640 // undo relocation 1641 if (norm_bef < norm_aft) 1642 relocate_one_vertex(v); 1643 } 1644 } 1645 1646 if (m_verbose > 0) 1647 std::cerr << "done" << " (" << timer.time() << " s)" << std::endl; 1648 } 1649 1650 /// @} 1651 1652 /// \name Output 1653 /// @{ 1654 1655 /*! 1656 Writes the points and segments of the output simplex in an indexed format into output iterators. 1657 \tparam PointOutputIterator An output iterator with value type 1658 \link Optimal_transportation_reconstruction_2::Point Point \endlink. 1659 \tparam IndexOutputIterator An output iterator with value type 1660 `std::size_t`. 1661 \tparam IndexPairOutputIterator An output iterator with value type 1662 `std::pair<std::size_t, std::size_t>`. 1663 1664 \param points The output iterator for all points. 1665 \param isolated_points The output iterator for the indices of isolated points. 1666 \param segments The output iterator for the pairs of segment indices. 1667 */ 1668 template < 1669 typename PointOutputIterator, 1670 typename IndexOutputIterator, 1671 typename IndexPairOutputIterator> 1672 std::tuple< 1673 PointOutputIterator, 1674 IndexOutputIterator, 1675 IndexPairOutputIterator> indexed_output(PointOutputIterator points,IndexOutputIterator isolated_points,IndexPairOutputIterator segments)1676 indexed_output( 1677 PointOutputIterator points, 1678 IndexOutputIterator isolated_points, 1679 IndexPairOutputIterator segments) const 1680 { 1681 std::vector<Point> isolated_points_; 1682 std::vector<Segment> edges; 1683 1684 list_output ( 1685 std::back_inserter(isolated_points_), std::back_inserter(edges)); 1686 1687 // vertices_of_edges 1688 std::set<Point> edge_vertices; 1689 for (typename std::vector<Segment>::iterator it = edges.begin(); 1690 it != edges.end(); it++) { 1691 1692 Point a = (*it).source(); 1693 Point b = (*it).target(); 1694 1695 edge_vertices.insert(a); 1696 edge_vertices.insert(b); 1697 } 1698 1699 std::size_t count_points = 0; 1700 for (typename std::set<Point>::iterator it = edge_vertices.begin(); 1701 it != edge_vertices.end(); it++) { 1702 1703 *points++ = *it; 1704 ++count_points; 1705 } 1706 1707 for (typename std::vector<Point>::iterator it = isolated_points_.begin(); 1708 it != isolated_points_.end(); it++) { 1709 1710 *isolated_points++ = count_points; 1711 *points++ = *it; 1712 ++count_points; 1713 } 1714 1715 for (typename std::vector<Segment>::iterator it = edges.begin(); 1716 it != edges.end(); it++) { 1717 1718 Point const& a = it->source(); 1719 Point const& b = it->target(); 1720 1721 typename std::set<Point>::iterator it_a = edge_vertices.find(a); 1722 typename std::set<Point>::iterator it_b = edge_vertices.find(b); 1723 1724 std::size_t pos_a = std::distance(edge_vertices.begin(), it_a); 1725 std::size_t pos_b = std::distance(edge_vertices.begin(), it_b); 1726 1727 *segments++ = std::make_pair(pos_a, pos_b); 1728 } 1729 1730 return std::make_tuple(points, isolated_points, segments); 1731 } 1732 1733 /*! 1734 Returns the solid edges and vertices present after the reconstruction 1735 process finished. 1736 1737 \details It takes two output iterators, one for storing the 1738 isolated points and one for storing the edges of the reconstructed shape. 1739 1740 \tparam PointOutputIterator An output iterator with value type 1741 \link Optimal_transportation_reconstruction_2::Point Point \endlink. 1742 \tparam SegmentOutputIterator An output iterator with value type 1743 \link Optimal_transportation_reconstruction_2::Segment Segment \endlink. 1744 */ 1745 template<class PointOutputIterator, class SegmentOutputIterator> list_output(PointOutputIterator v_it,SegmentOutputIterator e_it)1746 void list_output (PointOutputIterator v_it, SegmentOutputIterator e_it) const 1747 { 1748 for (Vertex_iterator vi = m_dt.vertices_begin(); 1749 vi != m_dt.vertices_end(); ++vi) 1750 { 1751 bool incident_edges_have_sample = false; 1752 typename Triangulation::Edge_circulator start = m_dt.incident_edges(vi); 1753 typename Triangulation::Edge_circulator cur = start; 1754 1755 do { 1756 if (!m_dt.is_ghost(*cur)) { 1757 incident_edges_have_sample = true; 1758 break; 1759 } 1760 ++cur; 1761 } while (cur != start); 1762 1763 if (!incident_edges_have_sample) { 1764 if ((*vi).has_sample_assigned()) { 1765 Point p = (*vi).point(); 1766 *v_it = p; 1767 v_it++; 1768 } 1769 } 1770 } 1771 1772 for (Finite_edges_iterator ei = m_dt.finite_edges_begin(); ei != m_dt.finite_edges_end(); ++ei) 1773 { 1774 Edge edge = *ei; 1775 if (m_dt.is_ghost(edge)) 1776 continue; 1777 1778 int index = edge.second; 1779 Vertex_handle source = edge.first->vertex( (index+1)%3 ); 1780 Vertex_handle target = edge.first->vertex( (index+2)%3 ); 1781 1782 Segment s = m_traits.construct_segment_2_object()( 1783 source->point(), target->point()); 1784 *e_it = s; 1785 e_it++; 1786 } 1787 } 1788 /// \endcond 1789 1790 1791 /// \cond SKIP_IN_MANUAL tds()1792 const Triangulation& tds() const { return m_dt; } 1793 extract_tds_output(Triangulation & rt2)1794 void extract_tds_output(Triangulation& rt2) const { 1795 rt2 = m_dt; 1796 //mark vertices 1797 for (Vertex_iterator vi = rt2.vertices_begin(); 1798 vi != rt2.vertices_end(); ++vi) 1799 { 1800 bool incident_edges_have_sample = false; 1801 typename Triangulation::Edge_circulator start = rt2.incident_edges(vi); 1802 typename Triangulation::Edge_circulator cur = start; 1803 1804 do { 1805 if (!rt2.is_ghost(*cur)) { 1806 incident_edges_have_sample = true; 1807 break; 1808 } 1809 ++cur; 1810 } while (cur != start); 1811 1812 if (!incident_edges_have_sample) { 1813 if ((*vi).has_sample_assigned()) 1814 (*vi).set_relevance(1); 1815 } 1816 } 1817 1818 // mark edges 1819 for (Finite_edges_iterator ei = rt2.finite_edges_begin(); ei != rt2.finite_edges_end(); ++ei) 1820 { 1821 Edge edge = *ei; 1822 FT relevance = 0; 1823 if (!rt2.is_ghost(edge)) { 1824 relevance = rt2.get_edge_relevance(edge); // >= 0 1825 } 1826 edge.first->relevance(edge.second) = relevance; 1827 } 1828 } 1829 1830 1831 /// \endcond 1832 /// @} 1833 1834 }; 1835 } // namespace 1836 1837 #endif // CGAL_OPTIMAL_TRANSPORTATION_RECONSTRUCTION_2_H_ 1838