From 58abe5a32a02c7ae42bd4db2044b519b23f07b48 Mon Sep 17 00:00:00 2001 From: Laurent Rineau Date: Mon, 30 Jan 2017 17:09:11 +0100 Subject: [PATCH 01/10] Do not use the nonlinear strategy from inside refine_balls() --- .../CGAL/Mesh_3/Protect_edges_sizing_field.h | 27 +++++++++++++++---- 1 file changed, 22 insertions(+), 5 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 988f08bb1ee..f594e1c481f 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 @@ -47,6 +47,19 @@ namespace internal { const double weight_modifier = .81; //0.9025;//0.81; const double distance_divisor = 2.1; const int max_nb_vertices_to_reevaluate_size = 10; + + const int refine_balls_max_nb_of_loops = 29; +// for the origins of `refine_balls_max_nb_of_loops`, that dates from the +// very beginning of this file: +// +// commit e9b3ff3e5730dab319a8cd581e3eb191559c98db +// Author: Stéphane Tayeb +// Date: Tue Apr 20 14:53:11 2010 +0000 +// +// Add draft of _with_features classes. +// +// That constant has had different values in the Git history: 9, 99, and now 29. + } // end namespace internal } // end namespace Mesh_3 } // end namespace CGAL @@ -345,6 +358,7 @@ private: Weight minimal_weight_; std::set treated_edges_; std::set unchecked_vertices_; + int refine_balls_iteration_nb; bool nonlinear_growth_of_balls; }; @@ -358,6 +372,7 @@ Protect_edges_sizing_field(C3T3& c3t3, const MD& domain, , size_(size) , minimal_size_(minimal_size) , minimal_weight_(CGAL::square(minimal_size)) + , refine_balls_iteration_nb(0) , nonlinear_growth_of_balls(false) { #ifndef CGAL_MESH_3_NO_PROTECTION_NON_LINEAR @@ -967,7 +982,7 @@ insert_balls(const Vertex_handle& vp, int n = static_cast(std::floor(FT(2)*(d-sq) / (sp+sq))+.5); // if( minimal_weight_ != 0 && n == 0 ) return; - if(nonlinear_growth_of_balls) { + if(nonlinear_growth_of_balls && refine_balls_iteration_nb == 0) { // This block tries not to apply the general rule that the size of // protecting balls is a linear interpolation of the size of protecting // balls at corner. When the curve segment is long enough, pick a point @@ -1113,14 +1128,16 @@ refine_balls() // Loop bool restart = true; - int nb=0; - while ( (!unchecked_vertices_.empty() || restart) && nb<29) + using CGAL::Mesh_3::internal::refine_balls_max_nb_of_loops; + this->refine_balls_iteration_nb = 0; + while ( (!unchecked_vertices_.empty() || restart) && + this->refine_balls_iteration_nb < refine_balls_max_nb_of_loops) { #if CGAL_MESH_3_PROTECTION_DEBUG & 1 - std::cerr << "RESTART REFINE LOOP (" << nb << ")\n" + std::cerr << "RESTART REFINE LOOP (" << refine_balls_iteration_nb << ")\n" << "\t unchecked_vertices size: " << unchecked_vertices_.size() <<"\n"; #endif - ++nb; + ++refine_balls_iteration_nb; restart = false; std::map new_sizes; From 2b932b5b71b5aa7e6dd54ef2d950a98adbe23e62 Mon Sep 17 00:00:00 2001 From: Jane Tournois Date: Thu, 2 Feb 2017 16:51:49 +0100 Subject: [PATCH 02/10] add dump of c3t3 at every step of protecting balls placement --- .../include/CGAL/Mesh_3/Protect_edges_sizing_field.h | 11 +++++++++++ 1 file changed, 11 insertions(+) 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 f594e1c481f..0d6101281ce 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 @@ -40,6 +40,10 @@ #if ! defined(CGAL_NO_PRECONDITIONS) # include #endif +#ifdef CGAL_MESH_3_DUMP_FEATURES_PROTECTION_ITERATIONS +#include +#endif + namespace CGAL { namespace Mesh_3 { namespace internal { @@ -1133,6 +1137,13 @@ refine_balls() while ( (!unchecked_vertices_.empty() || restart) && this->refine_balls_iteration_nb < refine_balls_max_nb_of_loops) { +#ifdef CGAL_MESH_3_DUMP_FEATURES_PROTECTION_ITERATIONS + std::ostringstream oss; + oss << "dump_protecting_balls_" << refine_balls_iteration_nb << ".cgal"; + std::ofstream outfile(oss.str(), std::ios_base::binary | std::ios_base::out); + CGAL::Mesh_3::save_binary_file(outfile, c3t3_, true); +#endif //CGAL_MESH_3_DUMP_FEATURES_PROTECTION_ITERATIONS + #if CGAL_MESH_3_PROTECTION_DEBUG & 1 std::cerr << "RESTART REFINE LOOP (" << refine_balls_iteration_nb << ")\n" << "\t unchecked_vertices size: " << unchecked_vertices_.size() <<"\n"; From 056309de5e929e5cf24de994921909e357cb669a Mon Sep 17 00:00:00 2001 From: Jane Tournois Date: Thu, 2 Feb 2017 17:00:38 +0100 Subject: [PATCH 03/10] allow to use the sizing field of curves more than once the protecting balls placement algorithm is now allowed to use the sizing field during 3 iterations, instead of 1. This allows it to fix most of the forbidden intersections of spheres before switching to the constant size case this commit relaxes the condition set in commit 051c55b08f2bfe5047eef1be1f70f0192b1c5c60 --- Mesh_3/include/CGAL/Mesh_3/Protect_edges_sizing_field.h | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) 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 0d6101281ce..e1eacc75d98 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 @@ -986,7 +986,8 @@ insert_balls(const Vertex_handle& vp, int n = static_cast(std::floor(FT(2)*(d-sq) / (sp+sq))+.5); // if( minimal_weight_ != 0 && n == 0 ) return; - if(nonlinear_growth_of_balls && refine_balls_iteration_nb == 0) { + if(nonlinear_growth_of_balls && refine_balls_iteration_nb < 3) + { // This block tries not to apply the general rule that the size of // protecting balls is a linear interpolation of the size of protecting // balls at corner. When the curve segment is long enough, pick a point From 83f74c04484d2d3f4b002782da6a649fea65dfba Mon Sep 17 00:00:00 2001 From: Laurent Rineau Date: Thu, 2 Feb 2017 16:22:24 +0100 Subject: [PATCH 04/10] Fix the sampling of protecting balls The previous code never verified that the curve is inside the union of balls. Now it does. --- .../CGAL/Mesh_3/Protect_edges_sizing_field.h | 32 ++++++++++++++++++- 1 file changed, 31 insertions(+), 1 deletion(-) 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 e1eacc75d98..8aa2c7060d4 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 @@ -1487,8 +1487,38 @@ is_sampling_dense_enough(const Vertex_handle& v1, const Vertex_handle& v2) const // Get sizes FT size_v1 = get_size(v1); FT size_v2 = get_size(v2); - FT distance_v1v2 = compute_distance(v1,v2); + + Curve_segment_index curve_index = Curve_segment_index(); + if(get_dimension(v1) == 1) { + curve_index = domain_.curve_segment_index(v1->index()); + } else if(get_dimension(v2) == 1) { + curve_index = domain_.curve_segment_index(v2->index()); + } + if(curve_index != Curve_segment_index()) { + FT geodesic_distance = domain_.geodesic_distance(v1->point().point(), + v2->point().point(), + curve_index); + if(domain_.is_cycle(v1->point().point(), curve_index)) { + geodesic_distance = + (std::min)(geodesic_distance, + domain_.geodesic_distance(v2->point().point(), + v1->point().point(), + curve_index)); + } else { + geodesic_distance = CGAL::abs(geodesic_distance); + } + // Sufficient condition so that the curve portion between v1 and v2 is + // inside the union of the two balls. + if(geodesic_distance > (size_v1 + size_v2)) { + std::cerr << "Note: on curve #" << curve_index << ", between (" + << v1->point() << ") and (" << v2->point() << "), the " + << "geodesic distance is " << geodesic_distance << " and the" + << " sum of radii is " << size_v1 + size_v2 << std::endl; + return false; + } + } + const FT distance_v1v2 = compute_distance(v1,v2); // Ensure size_v1 > size_v2 if ( size_v1 < size_v2 ) { std::swap(size_v1, size_v2); } From e023fc9157a5302c91f43f59bc6b64fa8565654f Mon Sep 17 00:00:00 2001 From: Jane Tournois Date: Thu, 2 Feb 2017 17:26:05 +0100 Subject: [PATCH 05/10] protect verbose code with macro CGAL_MESH_3_PROTECTION_DEBUG --- Mesh_3/include/CGAL/Mesh_3/Protect_edges_sizing_field.h | 2 ++ 1 file changed, 2 insertions(+) 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 8aa2c7060d4..d2462ee20f1 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 @@ -1510,10 +1510,12 @@ is_sampling_dense_enough(const Vertex_handle& v1, const Vertex_handle& v2) const // Sufficient condition so that the curve portion between v1 and v2 is // inside the union of the two balls. if(geodesic_distance > (size_v1 + size_v2)) { +#if CGAL_MESH_3_PROTECTION_DEBUG & 1 std::cerr << "Note: on curve #" << curve_index << ", between (" << v1->point() << ") and (" << v2->point() << "), the " << "geodesic distance is " << geodesic_distance << " and the" << " sum of radii is " << size_v1 + size_v2 << std::endl; +#endif return false; } } From e28b298c99c04d2891ed017e8c000c17385d9f97 Mon Sep 17 00:00:00 2001 From: Laurent Rineau Date: Fri, 17 Feb 2017 12:17:18 +0100 Subject: [PATCH 06/10] Refactoring of Sizing_field_with_aabb_tree (API breakage too) The refactoring allows to compute a better size of corners in a cycle. ... But the bug is still there on the nasty data set! --- .../Mesh_domain_with_polyline_features_3.h | 45 ++++++++++++++++--- 1 file changed, 39 insertions(+), 6 deletions(-) diff --git a/Mesh_3/include/CGAL/Mesh_domain_with_polyline_features_3.h b/Mesh_3/include/CGAL/Mesh_domain_with_polyline_features_3.h index fa6a8308b06..313f462f1e2 100644 --- a/Mesh_3/include/CGAL/Mesh_domain_with_polyline_features_3.h +++ b/Mesh_3/include/CGAL/Mesh_domain_with_polyline_features_3.h @@ -99,6 +99,16 @@ public: return start_point() == end_point(); } + /// Returns the angle at the first point. + /// \pre The polyline must be a cycle. + Angle angle_at_first_point() const { + CGAL_precondition(is_cycle()); + const Point_3& first = points_.front(); + const Point_3& next = points_[1]; + const Point_3& prev = points_[points_.size() - 2]; + return angle(prev, first, next); + } + /// Returns the length of the polyline FT length() const { @@ -576,7 +586,11 @@ public: template IndicesOutputIterator - get_corner_incidences(Curve_segment_index id, IndicesOutputIterator out) const; + get_corner_incidences(Corner_index id, IndicesOutputIterator out) const; + + template + IndicesOutputIterator + get_corner_incident_curves(Corner_index id, IndicesOutputIterator out) const; typedef std::set Surface_patch_index_set; @@ -862,6 +876,19 @@ get_corner_incidences(Corner_index id, return std::copy(incidences.begin(), incidences.end(), indices_out); } +template +template +IndicesOutputIterator +Mesh_domain_with_polyline_features_3:: +get_corner_incident_curves(Corner_index id, + IndicesOutputIterator indices_out) const +{ + typename Corners_tmp_incidences::const_iterator it = + corners_tmp_incidences_.find(id); + const std::set& incidences = it->second; + return std::copy(incidences.begin(), incidences.end(), indices_out); +} + namespace internal { namespace Mesh_3 { template @@ -960,14 +987,20 @@ compute_corners_incidences() corner_tmp_incidences = corners_tmp_incidences_[id]; // If the corner is incident to only one curve, and that curve is a - // cycle, then remove the corner from the set. + // cycle, then remove the corner from the set, only if the angle is not + // acute. If the angle is acute, the corner must remain as a corner, + // to deal correctly with the angle. if(corner_tmp_incidences.size() == 1 && is_cycle(Point_3(), *corner_tmp_incidences.begin())) { - typename Corners::iterator to_erase = cit; - ++cit; - corners_.erase(to_erase); - continue; + const Curve_segment_index curve_id = *corner_tmp_incidences.begin(); + const Polyline& polyline = edges_[curve_id]; + if(polyline.angle_at_first_point() == OBTUSE) { + typename Corners::iterator to_erase = cit; + ++cit; + corners_.erase(to_erase); + continue; + } } Surface_patch_index_set& incidences = corners_incidences_[id]; From 9668b1b6dd8f42eb28077f84660e26e61b75fb4c Mon Sep 17 00:00:00 2001 From: Jane Tournois Date: Fri, 17 Feb 2017 17:22:17 +0100 Subject: [PATCH 07/10] add missing abs --- Mesh_3/include/CGAL/Mesh_3/Protect_edges_sizing_field.h | 6 +++--- 1 file changed, 3 insertions(+), 3 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 d2462ee20f1..f203db51c8f 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 @@ -1500,10 +1500,10 @@ is_sampling_dense_enough(const Vertex_handle& v1, const Vertex_handle& v2) const curve_index); if(domain_.is_cycle(v1->point().point(), curve_index)) { geodesic_distance = - (std::min)(geodesic_distance, - domain_.geodesic_distance(v2->point().point(), + (std::min)(CGAL::abs(geodesic_distance), + CGAL::abs(domain_.geodesic_distance(v2->point().point(), v1->point().point(), - curve_index)); + curve_index))); } else { geodesic_distance = CGAL::abs(geodesic_distance); } From 9c88d16e7b37cc8d458236e716764dc6e3f2b1fb Mon Sep 17 00:00:00 2001 From: Laurent Rineau Date: Tue, 21 Feb 2017 11:53:46 +0100 Subject: [PATCH 08/10] Bug fix! --- Mesh_3/include/CGAL/Mesh_3/Protect_edges_sizing_field.h | 6 ++++++ 1 file changed, 6 insertions(+) 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 f203db51c8f..a916d0916ad 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 @@ -1773,6 +1773,12 @@ repopulate_edges_around_corner(const Vertex_handle& v, ErasedVeOutIt out) const Vertex_handle& next = vit->first; const Curve_segment_index& curve_index = vit->second; + // if `v` is incident to a cycle, it might be that the full cycle, + // including the edge `[next, v]`, has already been processed by + // `analyze_and_repopulate()` walking in the other direction. + if(domain_.is_cycle(v->point(), curve_index) && + !c3t3_.is_in_complex(v, next)) continue; + // Walk along each incident edge of the corner Vertex_vector to_repopulate; to_repopulate.push_back(v); From 58481453f0e1508a4aad745e1eb2f1546c4d47eb Mon Sep 17 00:00:00 2001 From: Laurent Rineau Date: Tue, 21 Feb 2017 11:54:15 +0100 Subject: [PATCH 09/10] More debug possibilities --- Mesh_3/include/CGAL/Mesh_3/Dump_c3t3.h | 16 ++++++++++++++++ .../CGAL/Mesh_3/Protect_edges_sizing_field.h | 16 ++++++++++++++++ 2 files changed, 32 insertions(+) diff --git a/Mesh_3/include/CGAL/Mesh_3/Dump_c3t3.h b/Mesh_3/include/CGAL/Mesh_3/Dump_c3t3.h index 2da23f6f4b0..6081a22a9bb 100644 --- a/Mesh_3/include/CGAL/Mesh_3/Dump_c3t3.h +++ b/Mesh_3/include/CGAL/Mesh_3/Dump_c3t3.h @@ -103,6 +103,22 @@ struct Dump_c3t3 { } }; // end struct template specialization Dump_c3t3 +template +void dump_c3t3_edges(const C3t3& c3t3, std::string prefix) { + std::ofstream file((prefix+".polylines.txt").c_str()); + file.precision(17); + for(typename C3t3::Edges_in_complex_iterator + edge_it = c3t3.edges_in_complex_begin(), + end = c3t3.edges_in_complex_end(); + edge_it != end; ++edge_it) + { + const typename C3t3::Triangulation::Cell_handle c = edge_it->first; + const int i = edge_it->second; + const int j = edge_it->third; + file << "2 " << c->vertex(i)->point().point() + << " " << c->vertex(j)->point().point() << "\n"; + } +} template void dump_c3t3(const C3t3& c3t3, std::string prefix) { if(!prefix.empty()) { 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 a916d0916ad..dc8a29898f7 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 @@ -1129,6 +1129,10 @@ void Protect_edges_sizing_field:: refine_balls() { +#if CGAL_MESH_3_PROTECTION_DEBUG & 4 + dump_c3t3(c3t3_, "dump-before-refine_balls"); + dump_c3t3_edges(c3t3_, "dump-before-refine_balls"); +#endif Triangulation& tr = c3t3_.triangulation(); // Loop @@ -1217,6 +1221,10 @@ refine_balls() } } +#if CGAL_MESH_3_PROTECTION_DEBUG & 4 + dump_c3t3(c3t3_, "dump-before-check_and_repopulate_edges"); + dump_c3t3_edges(c3t3_, "dump-before-check_and_repopulate_edges"); +#endif // Check edges check_and_repopulate_edges(); } @@ -1537,6 +1545,14 @@ walk_along_edge(const Vertex_handle& start, const Vertex_handle& next, bool /*test_sampling*/, ErasedVeOutIt out) const { +#if CGAL_MESH_3_PROTECTION_DEBUG & 4 + if(!c3t3_.is_in_complex(start, next)) { + std::cerr << "ERROR: the edge ( " << start->point() << " , " + << next->point() << " ) is not in complex!\n"; + dump_c3t3(c3t3_, "dump-bug"); + dump_c3t3_edges(c3t3_, "dump-bug-c3t3"); + } +#endif CGAL_precondition( c3t3_.is_in_complex(start, next) ); Vertex_handle previous = start; From e8a7391dd5ad88d974e5c224fb771af40e095238 Mon Sep 17 00:00:00 2001 From: Laurent Rineau Date: Wed, 1 Mar 2017 19:05:30 +0100 Subject: [PATCH 10/10] Fix testsuite That is a followup of e28b298c99c04d2891ed017e8c000c17385d9f97. --- Mesh_3/test/Mesh_3/test_domain_with_polyline_features.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Mesh_3/test/Mesh_3/test_domain_with_polyline_features.cpp b/Mesh_3/test/Mesh_3/test_domain_with_polyline_features.cpp index fb520ba3f40..5077d715ac0 100644 --- a/Mesh_3/test/Mesh_3/test_domain_with_polyline_features.cpp +++ b/Mesh_3/test/Mesh_3/test_domain_with_polyline_features.cpp @@ -105,7 +105,7 @@ public: Corners_vector corners; domain_.get_corners(std::back_inserter(corners)); - assert(corners.empty()); + assert(corners.size() == 1); } void test_curve_segments() const