From 4bdafe6908c67866c62c1fa0d12bb2c9fe5e447f Mon Sep 17 00:00:00 2001 From: Laurent Rineau Date: Wed, 26 Nov 2014 15:56:37 +0100 Subject: [PATCH 1/2] Bug-fix in Mesh_3 with time stamps That is a follow-up to that commit: | commit eefe2012c01a5d1c364b20ce5746f7e9285a37b8 | Author: Laurent Rineau | Date: Thu Oct 2 16:59:19 2014 +0200 | | Bug-fix in Mesh_3 with time stamps | | The implementation of the class template Protect_edges_sizing_field was | not correct when time stamps are used for Vertex_handle. The issue is | that the time stamp of a vertex change when one do: | tr.remove(v); | v = tr.insert(point); | The vertex pointed by 'v' is re-used by the TDS, but the time stamp of | the vertex does change. | The fix was not complete. --- .../CGAL/Mesh_3/Protect_edges_sizing_field.h | 246 +++++++++++------- 1 file changed, 153 insertions(+), 93 deletions(-) diff --git a/Mesh_3/include/CGAL/Mesh_3/Protect_edges_sizing_field.h b/Mesh_3/include/CGAL/Mesh_3/Protect_edges_sizing_field.h index fac4a50bdaf..3b6c4b2ac25 100644 --- a/Mesh_3/include/CGAL/Mesh_3/Protect_edges_sizing_field.h +++ b/Mesh_3/include/CGAL/Mesh_3/Protect_edges_sizing_field.h @@ -148,31 +148,38 @@ private: * It also ensures that no point of the triangulation will be inside its * sphere, by decreasing w. */ - Vertex_handle smart_insert_point(const Bare_point& p, - Weight w, - int dim, - const Index& index); + template + std::pair + smart_insert_point(const Bare_point& p, + Weight w, + int dim, + const Index& index, + ErasedVeOutIt out); /// Insert balls between points which are pointed by handles \c vp and \c vq /// on curve identified by \c curve_index - void insert_balls(const Vertex_handle& vp, - const Vertex_handle& vq, - const Curve_segment_index& curve_index); + template + ErasedVeOutIt insert_balls(const Vertex_handle& vp, + const Vertex_handle& vq, + const Curve_segment_index& curve_index, + ErasedVeOutIt out); /** * Insert balls * Preconditions: * - size_p < size_q * - pq_geodesic > 0 - */ - void insert_balls(const Vertex_handle& vp, - const Vertex_handle& vq, - const FT size_p, - const FT size_q, - const FT pq_geodesic, - const CGAL::Sign distance_sign, - const Curve_segment_index& curve_index); + */ + template + ErasedVeOutIt insert_balls(const Vertex_handle& vp, + const Vertex_handle& vq, + const FT size_p, + const FT size_q, + const FT pq_geodesic, + const CGAL::Sign distance_sign, + const Curve_segment_index& curve_index, + ErasedVeOutIt out); /// Returns true if balls of \c va and \c vb intersect, and (va,vb) is not /// an edge of the complex @@ -199,19 +206,19 @@ private: /// Checks if vertex \c v is well sampled, and if its not the case, fix it. /// Fills out with deleted vertices during this process. out value type /// is Vertex_handle. - template - OutputIterator - check_and_fix_vertex_along_edge(const Vertex_handle& v, OutputIterator out); + template + ErasedVeOutIt + check_and_fix_vertex_along_edge(const Vertex_handle& v, ErasedVeOutIt out); /// Walk along edge from \c start, following the direction \c start to /// \c next, and fills \c out with the vertices which do not fullfill /// the sampling conditions - template - OutputIterator + template + ErasedVeOutIt walk_along_edge(const Vertex_handle& start, const Vertex_handle& next, const bool test_sampling, - OutputIterator out) const; + ErasedVeOutIt out) const; /// Returns next vertex along edge, i.e vertex after \c start, following /// the direction from \c previous to \c start @@ -222,16 +229,18 @@ private: /// Replace vertices between ]begin,last[ by new vertices, along curve /// identified by \c curve_index /// The value type of InputIterator is Vertex_handle. - template - void repopulate(InputIterator begin, - InputIterator last, - const Curve_segment_index& index); + template + ErasedVeOutIt repopulate(InputIterator begin, + InputIterator last, + const Curve_segment_index& index, + ErasedVeOutIt out); - template - void + template + ErasedVeOutIt analyze_and_repopulate(InputIterator begin, InputIterator last, - const Curve_segment_index& index); + const Curve_segment_index& index, + ErasedVeOutIt out); /// Checks if \c v2 size is compatible (i.e. greater) with the linear /// interpolation of sizes of \c v1 and \c v3 @@ -241,9 +250,9 @@ private: /// Repopulate all incident curve around corner \c v /// \pre \c v is a corner of c3t3 - template - OutputIterator - repopulate_edges_around_corner(const Vertex_handle& v, OutputIterator out); + template + ErasedVeOutIt + repopulate_edges_around_corner(const Vertex_handle& v, ErasedVeOutIt out); /// Returns true if edge with index \c curve_index is already treated bool is_treated(const Curve_segment_index& curve_index) const @@ -368,6 +377,7 @@ operator()(const bool refine) std::cerr << "refine_balls() done. Nb of points in triangulation: " << c3t3_.triangulation().number_of_vertices() << std::endl; #endif + debug_dump_c3t3("dump-mesh-after-protect_edges.binary.cgal", c3t3_); CGAL_assertion(minimal_size_ > 0 || c3t3_.is_valid()); } @@ -442,7 +452,8 @@ insert_corners() } // Insert corner with ball (dim is zero because p is a corner) - Vertex_handle v = smart_insert_point(p, w, 0, p_index); + Vertex_handle v = smart_insert_point(p, w, 0, p_index, + CGAL::Emptyset_iterator()).first; CGAL_assertion(v != Vertex_handle()); // As C3t3::add_to_complex modifies the 'in_dimension' of the vertex, @@ -453,7 +464,7 @@ insert_corners() set_special(v); } } -} +} //end insert_corners() template @@ -516,9 +527,12 @@ insert_point(const Bare_point& p, const Weight& w, int dim, const Index& index, template -typename Protect_edges_sizing_field::Vertex_handle +template +std::pair::Vertex_handle, + ErasedVeOutIt> Protect_edges_sizing_field:: -smart_insert_point(const Bare_point& p, Weight w, int dim, const Index& index) +smart_insert_point(const Bare_point& p, Weight w, int dim, const Index& index, + ErasedVeOutIt out) { #ifdef CGAL_MESH_3_PROTECTION_DEBUG std::cerr << "smart_insert_point( (" << p @@ -558,6 +572,7 @@ smart_insert_point(const Bare_point& p, Weight w, int dim, const Index& index) } // Adapt size + *out++ = nearest_vh; Vertex_handle new_vh = change_ball_size(nearest_vh, CGAL::sqrt(sq_d), special_ball); ch = tr.locate(p, lt, li, lj, new_vh); @@ -575,8 +590,9 @@ smart_insert_point(const Bare_point& p, Weight w, int dim, const Index& index) // Change w in order to be sure that no existing point will be included // in (p,w) - std::set vertices_in_conflict_zone; + std::vector vertices_in_conflict_zone; { // fill vertices_in_conflict_zone + std::set vertices_in_conflict_zone_set; std::vector cells_in_conflicts; tr.find_conflicts(Weighted_point(p, w), ch, CGAL::Emptyset_iterator(), @@ -590,16 +606,19 @@ smart_insert_point(const Bare_point& p, Weight w, int dim, const Index& index) for(int i = 0, d = tr.dimension(); i <= d; ++i) { const Vertex_handle v = (*it)->vertex(i); if( ! c3t3_.triangulation().is_infinite(v) ) { - vertices_in_conflict_zone.insert(v); + vertices_in_conflict_zone_set.insert(v); } } } + vertices_in_conflict_zone.insert(vertices_in_conflict_zone.end(), + vertices_in_conflict_zone_set.begin(), + vertices_in_conflict_zone_set.end()); } FT min_sq_d = w; #ifdef CGAL_MESH_3_PROTECTION_DEBUG typename Tr::Point nearest_point; #endif - for(typename std::set::const_iterator + for(typename std::vector::const_iterator it = vertices_in_conflict_zone.begin(), end = vertices_in_conflict_zone.end(); it != end ; ++it ) { @@ -611,6 +630,7 @@ smart_insert_point(const Bare_point& p, Weight w, int dim, const Index& index) #endif min_sq_d = minimal_weight_; if( ! is_special(*it) ) { + *out++ = *it; ch = change_ball_size(*it, minimal_size_, true)->cell(); // special ball } } else { @@ -666,6 +686,7 @@ smart_insert_point(const Bare_point& p, Weight w, int dim, const Index& index) insert_a_special_ball = true; } if( ! is_special(it) ) { + *out++ = it; change_ball_size(it, CGAL::sqrt(sq_d), special_ball); restart = true; } @@ -714,9 +735,10 @@ smart_insert_point(const Bare_point& p, Weight w, int dim, const Index& index) insert_a_special_ball = true; } Vertex_handle v = insert_point(p,w,dim,index, insert_a_special_ball); + /// @TODO `insert_point` does insert in unchecked_vertices anyway! if ( add_handle_to_unchecked ) { unchecked_vertices_.insert(v); } - return v; + return std::pair(v, out); } @@ -779,7 +801,8 @@ insert_balls_on_edges() vp = smart_insert_point(p, CGAL::square(p_size), 1, - p_index); + p_index, + CGAL::Emptyset_iterator()).first; } // No 'else' because in that case 'is_vertex(..)' already filled // the variable 'vp'. @@ -792,7 +815,7 @@ insert_balls_on_edges() // } // else { - insert_balls(vp, vq, curve_index); + insert_balls(vp, vq, curve_index, Emptyset_iterator()); } set_treated(curve_index); } @@ -800,7 +823,7 @@ insert_balls_on_edges() // s << "dump-mesh-curve-" << curve_index << ".binary.cgal"; // debug_dump_c3t3(s.str(), c3t3_); } -} +} //end insert_balls_on_edges() template @@ -820,11 +843,13 @@ get_vertex_corner_from_point(const Bare_point& p, const Index&) const template -void +template +ErasedVeOutIt Protect_edges_sizing_field:: insert_balls(const Vertex_handle& vp, const Vertex_handle& vq, - const Curve_segment_index& curve_index) + const Curve_segment_index& curve_index, + ErasedVeOutIt out) { // Get size of p & q const Bare_point& p = vp->point().point(); @@ -839,14 +864,17 @@ insert_balls(const Vertex_handle& vp, const FT pq_geo = CGAL::abs(pq_geo_signed); // Insert balls - return (sp <= sq) ? insert_balls(vp, vq, sp, sq, pq_geo, d_sign, curve_index) - : insert_balls(vq, vp, sq, sp, pq_geo, -d_sign, curve_index); + return + (sp <= sq) ? + insert_balls(vp, vq, sp, sq, pq_geo, d_sign, curve_index, out) : + insert_balls(vq, vp, sq, sp, pq_geo, -d_sign, curve_index, out); } template -void +template +ErasedVeOutIt Protect_edges_sizing_field:: insert_balls(const Vertex_handle& vp, const Vertex_handle& vq, @@ -854,7 +882,8 @@ insert_balls(const Vertex_handle& vp, const FT sq, const FT d, const CGAL::Sign d_sign, - const Curve_segment_index& curve_index) + const Curve_segment_index& curve_index, + ErasedVeOutIt out) { #ifdef CGAL_MESH_3_PROTECTION_DEBUG std::cerr << "insert_balls(vp=" << (void*)(&*vp) << " (" << vp->point() << "),\n" @@ -935,22 +964,26 @@ insert_balls(const Vertex_handle& vp, std::cerr << " middle point: " << new_point << std::endl; std::cerr << " new weight: " << point_weight << std::endl; #endif - const Vertex_handle new_vertex = smart_insert_point(new_point, - point_weight, - dim, - index); + std::pair pair = + smart_insert_point(new_point, + point_weight, + dim, + index, + out); + const Vertex_handle new_vertex = pair.first; + out = pair.second; const FT sn = get_size(new_vertex); if(sp <= sn) { - insert_balls(vp, new_vertex, sp, sn, d/2, d_sign, curve_index); + out=insert_balls(vp, new_vertex, sp, sn, d/2, d_sign, curve_index, out); } else { - insert_balls(new_vertex, vp, sn, sp, d/2, -d_sign, curve_index); + out=insert_balls(new_vertex, vp, sn, sp, d/2, -d_sign, curve_index, out); } if(sn <= sq) { - insert_balls(new_vertex, vq, sn, sq, d/2, d_sign, curve_index); + out=insert_balls(new_vertex, vq, sn, sq, d/2, d_sign, curve_index, out); } else { - insert_balls(vq, new_vertex, sq, sn, d/2, -d_sign, curve_index); + out=insert_balls(vq, new_vertex, sq, sn, d/2, -d_sign, curve_index, out); } - return; + return out; } #endif // not CGAL_MESH_3_NO_PROTECTION_NON_LINEAR @@ -1009,7 +1042,10 @@ insert_balls(const Vertex_handle& vp, int dim = 1; // new_point is on edge // Insert point into c3t3 - Vertex_handle new_vertex = smart_insert_point(new_point, point_weight, dim, index); + std::pair pair = + smart_insert_point(new_point, point_weight, dim, index, out); + Vertex_handle new_vertex = pair.first; + out = pair.second; // Add edge to c3t3 if(!c3t3_.is_in_complex(prev, new_vertex)) { @@ -1034,6 +1070,7 @@ insert_balls(const Vertex_handle& vp, c3t3_.add_to_complex(prev, vq, curve_index); } } + return out; } @@ -1070,6 +1107,7 @@ refine_balls() // Compute correct size of balls const FT ab = compute_distance(va,vb); + /// @TOTO pb: get_size(va) is called several times FT sa_new = (std::min)(ab/distance_divisor, get_size(va)); FT sb_new = (std::min)(ab/distance_divisor, get_size(vb)); @@ -1088,21 +1126,34 @@ refine_balls() { new_sizes[vb] = sb_new; } } } + CGAL_assertion(std::is_sorted(new_sizes.begin(), + new_sizes.end())); + + // The std::map with Vertex_handle as the key is not robust, because + // the time stamp of vertices can change during the following loop. The + // solution is to copy it in a vector. + std::vector > + new_sizes_copy(new_sizes.begin(), new_sizes.end()); + new_sizes.clear(); // Update size of balls - for ( typename std::map::iterator it = new_sizes.begin(), - end = new_sizes.end() ; it != end ; ++it) + for ( typename std::vector >::iterator + it = new_sizes_copy.begin(), + end = new_sizes_copy.end(); + it != end ; ++it ) { + const Vertex_handle v = it->first; + const FT new_size = it->second; // Set size of the ball to new value - if(minimal_size_ != FT() && it->second < minimal_size_) { - if(!is_special(it->first)) { - change_ball_size(it->first,minimal_size_,true); // special ball + if(minimal_size_ != FT() && new_size < minimal_size_) { + if(!is_special(v)) { + change_ball_size(v,minimal_size_,true); // special ball // Loop will have to be run again restart = true; } } else { - change_ball_size(it->first,it->second); + change_ball_size(v,new_size); // Loop will have to be run again restart = true; @@ -1112,7 +1163,7 @@ refine_balls() // Check edges check_and_repopulate_edges(); } -} +} // end refine_balls() template @@ -1188,7 +1239,8 @@ change_ball_size(const Vertex_handle& v, const FT size, const bool special_ball) Bare_point p = v->point().point(); // Remove v from corners - boost::optional corner_index; + boost::optional corner_index = + boost::make_optional(false, Corner_index()); if ( c3t3_.is_in_complex(v) ) { corner_index = c3t3_.corner_index(v); @@ -1282,15 +1334,17 @@ check_and_repopulate_edges() check_and_fix_vertex_along_edge(v, boost::make_function_output_iterator(erase_from_vertices)); + CGAL_assertion(std::is_sorted(vertices.begin(), + vertices.end())); } } template -template -OutputIterator +template +ErasedVeOutIt Protect_edges_sizing_field:: -check_and_fix_vertex_along_edge(const Vertex_handle& v, OutputIterator out) +check_and_fix_vertex_along_edge(const Vertex_handle& v, ErasedVeOutIt out) { #ifdef CGAL_MESH_3_PROTECTION_DEBUG std::cerr << "check_and_fix_vertex_along_edge(" @@ -1353,14 +1407,15 @@ check_and_fix_vertex_along_edge(const Vertex_handle& v, OutputIterator out) } // Store erased vertices - std::copy(to_repopulate.begin(), to_repopulate.end(), out); + // out = std::copy(to_repopulate.begin(), to_repopulate.end(), out); // Repopulate edge const Curve_segment_index& curve_index = incident_vertices.front().second; - analyze_and_repopulate(to_repopulate.begin(), - --to_repopulate.end(), - curve_index); + out = analyze_and_repopulate(to_repopulate.begin(), + --to_repopulate.end(), + curve_index, + out); return out; } @@ -1388,12 +1443,12 @@ is_sampling_dense_enough(const Vertex_handle& v1, const Vertex_handle& v2) const template -template -OutputIterator +template +ErasedVeOutIt Protect_edges_sizing_field:: walk_along_edge(const Vertex_handle& start, const Vertex_handle& next, bool /*test_sampling*/, - OutputIterator out) const + ErasedVeOutIt out) const { CGAL_precondition( c3t3_.is_in_complex(start, next) ); @@ -1448,10 +1503,11 @@ next_vertex_along_edge(const Vertex_handle& start, template -template -void +template +ErasedVeOutIt Protect_edges_sizing_field:: -repopulate(InputIterator begin, InputIterator last, const Curve_segment_index& index) +repopulate(InputIterator begin, InputIterator last, + const Curve_segment_index& index, ErasedVeOutIt out) { #ifdef CGAL_MESH_3_PROTECTION_DEBUG std::cerr << "repopulate(begin=" << (void*)(&**begin) @@ -1465,7 +1521,7 @@ repopulate(InputIterator begin, InputIterator last, const Curve_segment_index& i CGAL_assertion( std::distance(begin,last) >= 0 ); // May happen - if ( begin == last ) { return; } + if ( begin == last ) { return out; } // Valid because begin < last InputIterator current = begin; @@ -1502,6 +1558,7 @@ repopulate(InputIterator begin, InputIterator last, const Curve_segment_index& i } std::cerr << c3t3_.index(*current) << std::endl; #endif // CGAL_MESH_3_PROTECTION_DEBUG + *out++ = *current; c3t3_.triangulation().remove(*current); } @@ -1517,15 +1574,16 @@ repopulate(InputIterator begin, InputIterator last, const Curve_segment_index& i } // Repopulate edge - insert_balls(*begin, *last, index); + return insert_balls(*begin, *last, index, out); } template -template -void +template +ErasedVeOutIt Protect_edges_sizing_field:: -analyze_and_repopulate(InputIterator begin, InputIterator last, const Curve_segment_index& index) +analyze_and_repopulate(InputIterator begin, InputIterator last, + const Curve_segment_index& index, ErasedVeOutIt out) { #ifdef CGAL_MESH_3_PROTECTION_DEBUG std::cerr << "analyze_and_repopulate(begin=" << (void*)(&**begin) @@ -1539,11 +1597,11 @@ analyze_and_repopulate(InputIterator begin, InputIterator last, const Curve_segm CGAL_assertion( std::distance(begin,last) >= 0 ); // May happen - if ( begin == last ) { return; } + if ( begin == last ) { return out; } if ( std::distance(begin,last) == 1 ) { - repopulate(begin, last, index); - return; + out = repopulate(begin, last, index, out); + return out; } // Here std::distance(begin,last) > 1 @@ -1588,9 +1646,10 @@ analyze_and_repopulate(InputIterator begin, InputIterator last, const Curve_segm InputIterator next = ch_stack.top(); ch_stack.pop(); // Iterators are on the reverse order in the stack, thus use [next,current] - repopulate(next, current, index); + out = repopulate(next, current, index, out); current = next; - } + } + return out; } template @@ -1611,10 +1670,10 @@ is_sizing_field_correct(const Vertex_handle& v1, template -template -OutputIterator +template +ErasedVeOutIt Protect_edges_sizing_field:: -repopulate_edges_around_corner(const Vertex_handle& v, OutputIterator out) +repopulate_edges_around_corner(const Vertex_handle& v, ErasedVeOutIt out) { CGAL_precondition(c3t3_.is_in_complex(v)); @@ -1636,7 +1695,8 @@ repopulate_edges_around_corner(const Vertex_handle& v, OutputIterator out) std::copy(to_repopulate.begin(), to_repopulate.end(), out); // Repopulate - analyze_and_repopulate(to_repopulate.begin(), --to_repopulate.end(), curve_index); + out = analyze_and_repopulate(to_repopulate.begin(), --to_repopulate.end(), + curve_index, out); } return out; From 0c6153ce8c521fe31693e1283e06d03e43d3ad18 Mon Sep 17 00:00:00 2001 From: Laurent Rineau Date: Thu, 27 Nov 2014 10:49:17 +0100 Subject: [PATCH 2/2] cleanup I forgot to remove extra expensive assertions that I have used using the debugging! I clean up also a few spaces at end of lines. --- .../CGAL/Mesh_3/Protect_edges_sizing_field.h | 13 ++++--------- 1 file changed, 4 insertions(+), 9 deletions(-) diff --git a/Mesh_3/include/CGAL/Mesh_3/Protect_edges_sizing_field.h b/Mesh_3/include/CGAL/Mesh_3/Protect_edges_sizing_field.h index 3b6c4b2ac25..500153c9dfc 100644 --- a/Mesh_3/include/CGAL/Mesh_3/Protect_edges_sizing_field.h +++ b/Mesh_3/include/CGAL/Mesh_3/Protect_edges_sizing_field.h @@ -234,7 +234,7 @@ private: InputIterator last, const Curve_segment_index& index, ErasedVeOutIt out); - + template ErasedVeOutIt analyze_and_repopulate(InputIterator begin, @@ -377,7 +377,6 @@ operator()(const bool refine) std::cerr << "refine_balls() done. Nb of points in triangulation: " << c3t3_.triangulation().number_of_vertices() << std::endl; #endif - debug_dump_c3t3("dump-mesh-after-protect_edges.binary.cgal", c3t3_); CGAL_assertion(minimal_size_ > 0 || c3t3_.is_valid()); } @@ -618,7 +617,7 @@ smart_insert_point(const Bare_point& p, Weight w, int dim, const Index& index, #ifdef CGAL_MESH_3_PROTECTION_DEBUG typename Tr::Point nearest_point; #endif - for(typename std::vector::const_iterator + for(typename std::vector::const_iterator it = vertices_in_conflict_zone.begin(), end = vertices_in_conflict_zone.end(); it != end ; ++it ) { @@ -964,8 +963,8 @@ insert_balls(const Vertex_handle& vp, std::cerr << " middle point: " << new_point << std::endl; std::cerr << " new weight: " << point_weight << std::endl; #endif - std::pair pair = - smart_insert_point(new_point, + std::pair pair = + smart_insert_point(new_point, point_weight, dim, index, @@ -1126,8 +1125,6 @@ refine_balls() { new_sizes[vb] = sb_new; } } } - CGAL_assertion(std::is_sorted(new_sizes.begin(), - new_sizes.end())); // The std::map with Vertex_handle as the key is not robust, because // the time stamp of vertices can change during the following loop. The @@ -1334,8 +1331,6 @@ check_and_repopulate_edges() check_and_fix_vertex_along_edge(v, boost::make_function_output_iterator(erase_from_vertices)); - CGAL_assertion(std::is_sorted(vertices.begin(), - vertices.end())); } }