do not modify the list of halfedges while iterating over it

Mean_curvature_skeleton/include/CGAL/Mean_curvature_flow_skeletonization.h

@@ -1159,50 +1159,59 @@ private:
   /// Split triangles with an angle greater than `alpha_TH`.
   std::size_t split_flat_triangles()
   {
-    int ne = 2 * static_cast<int>(num_edges(m_tmesh));
     compute_incident_angle();
 
-    std::size_t cnt = 0;
-    /// \todo this is unsafe, we loop over a sequence that we modify!!!
-    BOOST_FOREACH(halfedge_descriptor hd, halfedges(m_tmesh))
+    // collect edges that should be split because
+    // both opposite angle are larger than the
+    // threshold
+    std::vector<edge_descriptor> edges_to_split;
+    BOOST_FOREACH(edge_descriptor ed, edges(m_tmesh))
     {
-      halfedge_descriptor ei = hd;
+      halfedge_descriptor ei = halfedge(ed, m_tmesh);
+      halfedge_descriptor ej = opposite(ei, m_tmesh);
+      int ei_id = static_cast<int>(get(m_hedge_id_pmap, ei));
+      int ej_id = static_cast<int>(get(m_hedge_id_pmap, ej));
+
+      double angle_i = m_halfedge_angle[ei_id];
+      double angle_j = m_halfedge_angle[ej_id];
+      if (angle_i >= m_alpha_TH && angle_j >= m_alpha_TH)
+        edges_to_split.push_back( ed );
+    }
+
+    // now split the edge
+    std::size_t cnt = 0;
+    BOOST_FOREACH(edge_descriptor ed, edges_to_split)
+    {
+      halfedge_descriptor ei = halfedge(ed, m_tmesh);
       halfedge_descriptor ej = opposite(ei, m_tmesh);
       int ei_id = static_cast<int>(get(m_hedge_id_pmap, ei));
       int ej_id = static_cast<int>(get(m_hedge_id_pmap, ej));
-      if (ei_id < 0 || ei_id >= ne
-       || ej_id < 0 || ej_id >= ne)
-      {
-        continue;
-      }
 
       vertex_descriptor vs = source(ei, m_tmesh);
       vertex_descriptor vt = target(ei, m_tmesh);
 
       double angle_i = m_halfedge_angle[ei_id];
       double angle_j = m_halfedge_angle[ej_id];
-      if (angle_i < m_alpha_TH || angle_j < m_alpha_TH)
+
+      halfedge_descriptor ek = next(angle_i > angle_j ? ei : ej, m_tmesh);
+      vertex_descriptor vk = target(ek, m_tmesh);
+
+      // split the edge
+      halfedge_descriptor en = m_tmesh.split_edge(ei);
+      // split the incident faces
+      Euler::split_face(en, next(ei,m_tmesh), m_tmesh);
+      if (! is_border(ej,m_tmesh))
       {
-        continue;
+        Euler::split_face(ej, next(next(ej,m_tmesh),m_tmesh), m_tmesh);
       }
 
-      halfedge_descriptor ek;
-      if (angle_i > angle_j)
-      {
-        ek = next(ei, m_tmesh);
-      }
-      else
-      {
-        ek = next(ej, m_tmesh);
-      }
-      vertex_descriptor vk = target(ek, m_tmesh);
-      halfedge_descriptor en = internal::mesh_split(m_tmesh, ei);
       // set id for new vertex
       put(m_vertex_id_pmap, target(en,m_tmesh), m_vertex_id_count++);
       Point pn = project_vertex(vs, vt, vk, target(en, m_tmesh));
       // set point of new vertex
       put(m_tmesh_point_pmap, target(en,m_tmesh), pn);
-      cnt++;
+      target(en,m_tmesh)->vertices.clear(); // do no copy the info
+      ++cnt;
     }
     return cnt;
   }

Mean_curvature_skeleton/include/CGAL/extract_mean_curvature_flow_skeleton.h

@@ -51,7 +51,6 @@ namespace CGAL{
 ///        and the set of input vertices that contracted to `vd` can be retrieved
 ///        using `skeleton[vd].point` and `skeleton[vd].vertices` respectively.
 ///
-/// \todo I need to tweak SkeletonVertexVerticesMap to match the documentation
 template <class TriangleMesh>
 void extract_mean_curvature_flow_skeleton(const TriangleMesh& tmesh,
                                           typename Mean_curvature_flow_skeletonization<TriangleMesh>::Skeleton& skeleton)

Mean_curvature_skeleton/include/CGAL/internal/Mean_curvature_skeleton/Utility.h

@@ -35,46 +35,6 @@
 namespace CGAL {
 namespace internal {
 
-/**
-* Split the edge
-* @param hg the mesh containing the given edge
-* @param ei the edge to be split
-*/
-template<class TriangleMesh>
-typename boost::graph_traits<TriangleMesh>::halfedge_descriptor
-mesh_split(TriangleMesh& hg,
-           typename boost::graph_traits<TriangleMesh>::halfedge_descriptor ei)
-{
-  typedef typename boost::graph_traits<TriangleMesh>::halfedge_descriptor            halfedge_descriptor;
-
-
-  // halfedge_descriptor en = Euler::split_edge(ei, hg); // there is an issue in this function for now use the polyhedron version in the meantime
-  halfedge_descriptor en = hg.split_edge(ei);
-  en->vertex()->vertices.clear();
-  Euler::split_face(en, next(ei,hg), hg);
-
-  en->id() = -1;
-  opposite(en,hg)->id() = -1;
-  ei->id() = -1;
-  opposite(ei,hg)->id() = -1;
-  next(en,hg)->id() = -1;
-  opposite(next(en,hg),hg)->id() = -1;
-  next(next(en,hg),hg)->id() = -1;
-  next(en,hg)->id() = -1; // AF: the same as 3 lines above? error or duplicate?
-  halfedge_descriptor ej = opposite(en, hg);
-  if (! is_border(ej,hg))
-  {
-    Euler::split_face(opposite(ei,hg), next(ej,hg), hg);
-    next(ej,hg)->id() = -1;
-    halfedge_descriptor ei_op_next = next(opposite(ei,hg),hg);
-    ei_op_next->id() = -1;
-    opposite(ei_op_next,hg)->id() = -1;
-    next(ei_op_next,hg)->id() = -1;
-  }
-
-  return en;
-}
-
 template<class TriangleMesh, class TriangleMeshPointPMap, class Traits>
 double get_surface_area(TriangleMesh& hg, TriangleMeshPointPMap& hg_point_pmap, const Traits& traits)
 {