upgrade test to BGL API

Ridges_3/include/CGAL/Ridges.h

@@ -29,7 +29,6 @@
 #include <CGAL/barycenter.h>
 #include <CGAL/property_map.h>
 #include <CGAL/assertions.h>
-#include <CGAL/Vertex2Data_Property_Map_with_std_map.h>
 #include <boost/type_traits/is_same.hpp>
 
 namespace CGAL {

Ridges_3/test/Ridges_3/PolyhedralSurf.h

@@ -1,17 +1,19 @@
-#ifndef _POLYHEDRALSURF_H_
+#ifndef CGAL_POLYHEDRALSURF_H_
-#define _POLYHEDRALSURF_H_
+#define CGAL_POLYHEDRALSURF_H_
 
+#include <CGAL/Simple_cartesian.h>
+#include <CGAL/Polyhedron_3.h>
+#include <CGAL/IO/Polyhedron_iostream.h>
+#include <CGAL/boost/graph/graph_traits_Polyhedron_3.h>
+#include <CGAL/boost/graph/properties_Polyhedron_3.h>
+#include <CGAL/boost/graph/helpers.h>
 #include <cstdlib>
 #include <cstdio>
 #include <algorithm>
 #include <vector>
 #include <list>
 
-#include <CGAL/Cartesian.h>
+#include <boost/foreach.hpp>
-#include <CGAL/Polyhedron_3.h>
-#include <CGAL/IO/Polyhedron_iostream.h>
-
-#include "PolyhedralSurf_operations.h" 
 
 //----------------------------------------------------------------
 // A redefined items class for the Polyhedron_3 with
@@ -29,11 +31,11 @@ public:
  typedef typename FGeomTraits::Vector_3 Vector_3;
 
 protected:
-  Vector_3 normal;
+  //Vector_3 normal;
 public:
   My_facet() {}
-  const Vector_3 & getUnitNormal() const { return normal; }
+  //const Vector_3 & getUnitNormal() const { std::cerr << "coucou" << std::endl;return normal; }
-  void setNormal(Vector_3 n) { normal = n; }
+  //void setNormal(Vector_3 n) { normal = n; }
 };
 
 //------------------------------------------------
@@ -58,20 +60,39 @@ struct Wrappers_VFH:public CGAL::Polyhedron_items_3 {
 //PolyhedralSurf with facet normal operations
 //------------------------------------------------
 typedef double                FT;
-typedef CGAL::Cartesian<FT>  Kernel;
+typedef CGAL::Simple_cartesian<FT>  Kernel;
 typedef CGAL::Polyhedron_3 < Kernel, Wrappers_VFH > Polyhedron;
 typedef Kernel::Vector_3 Vector_3;
 
-class PolyhedralSurf:public Polyhedron {
+class PolyhedralSurf;
+
+namespace boost {
+  template <>
+  struct graph_traits<PolyhedralSurf> : public boost::graph_traits<Polyhedron>
+  {};  
+
+  template <>
+  struct graph_traits<PolyhedralSurf const> : public boost::graph_traits<Polyhedron>
+  {};
+
+  template <class Tag>
+  struct property_map<PolyhedralSurf,Tag> : public property_map<Polyhedron,Tag>
+  {};
+  template <class Tag>
+  struct property_map<const PolyhedralSurf,Tag> : public property_map<const Polyhedron,Tag>
+  {};
+}
+
+class PolyhedralSurf : public Polyhedron {
 public:
+  typedef boost::graph_traits<PolyhedralSurf>::vertex_descriptor vertex_descriptor;
+  typedef boost::graph_traits<PolyhedralSurf>::face_descriptor face_descriptor;
+  typedef boost::graph_traits<PolyhedralSurf>::halfedge_descriptor halfedge_descriptor;
+
   PolyhedralSurf() {}
-  void compute_facets_normals();
+
-  const Vector_3 computeFacetsAverageUnitNormal(const Vertex_const_handle v);
 };
 
+
+
 #endif
-
-
-
-
-

Ridges_3/test/Ridges_3/PolyhedralSurf_rings.h

@@ -1,28 +1,28 @@
-#ifndef _POLYSURF_RINGS_H_
+#ifndef CGAL_PSURF_RINGS_H_
-#define _POLYSURF_RINGS_H_
+#define CGAL_PSURF_RINGS_H_
+
+#include <cassert>
+#include <vector>
+#include <map>
 
-using namespace std;
 
 //---------------------------------------------------------------------------
 //T_PolyhedralSurf_rings
 //---------------------------------------------------------------------------
-template < class TPoly > class T_PolyhedralSurf_rings
+template < class TPoly >
+class T_PolyhedralSurf_rings
 {
+
+  const TPoly& P;
+
 protected:
   //Polyhedron
-  typedef typename TPoly::Vertex_const_handle                     Vertex_const_handle;
+  typedef typename boost::graph_traits<TPoly>::vertex_descriptor                     Vertex_const_handle;
-  typedef typename TPoly::Halfedge_const_handle                   Halfedge_const_handle;
+  typedef typename boost::graph_traits<TPoly>::halfedge_descriptor                   Halfedge_const_handle;
-  typedef typename TPoly::Facet_const_handle                      Facet_const_handle;
+  typedef typename boost::graph_traits<TPoly>::vertex_iterator                   Vertex_const_iterator;
-  typedef typename TPoly::Halfedge_around_vertex_const_circulator Halfedge_around_vertex_const_circulator;
+  typedef CGAL::Halfedge_around_target_circulator<TPoly> Halfedge_around_vertex_const_circulator;
-  typedef typename TPoly::Vertex_const_iterator                   Vertex_const_iterator;
 
-  //tag to visit vertices
+  typedef std::map<Vertex_const_handle, int> Vertex2int_map;
-  struct Vertex_cmp{//comparison is wrt vertex addresses
-    bool operator()(Vertex_const_handle a,  Vertex_const_handle b) const{
-      return &*a < &*b;
-    }
-  };
-  typedef std::map<Vertex_const_handle, int, Vertex_cmp> Vertex2int_map;
   Vertex2int_map ring_index_map;
 
   //vertex indices are initialised to -1
@@ -60,10 +60,12 @@ protected:
 template < class TPoly >
 T_PolyhedralSurf_rings <TPoly>::
 T_PolyhedralSurf_rings(const TPoly& P)
+  : P(P)
 {
   //init the ring_index_map
-  Vertex_const_iterator itb = P.vertices_begin(), ite = P.vertices_end();
+  Vertex_const_iterator itb, ite;
-  for(;itb!=ite;itb++) ring_index_map[itb] = -1; 
+  boost::tie(itb,ite) = vertices(P);
+  for(;itb!=ite;itb++) ring_index_map[*itb] = -1;
 }
 
 template < class TPoly >
@@ -74,11 +76,11 @@ push_neighbours_of(const Vertex_const_handle start, const int ith,
 {
   Vertex_const_handle v;
   Halfedge_around_vertex_const_circulator
-    hedgeb = start->vertex_begin(), hedgee = hedgeb;
+    hedgeb(halfedge(start,P),P), hedgee = hedgeb;
 
  CGAL_For_all(hedgeb, hedgee)
   {
-    v = hedgeb->opposite()->vertex();
+    v = target(opposite(*hedgeb,P),P);
     if (ring_index_map[v] != -1)  continue;//if visited: next
 
     ring_index_map[v] = ith;

Ridges_3/test/Ridges_3/compute_normals.h

@@ -0,0 +1,62 @@
+#ifndef COMPUTE_NORMALS_H
+#define COMPUTE_NORMALS_H
+
+#include <CGAL/boost/graph/helpers.h>
+
+#include <boost/foreach.hpp>
+
+template <typename TriangleMesh, typename FaceVectorMap, typename Kernel>
+const typename Kernel::Vector_3 
+computeFacetsAverageUnitNormal(const TriangleMesh& tm,
+                               typename boost::graph_traits<TriangleMesh>::vertex_descriptor v,
+                               FaceVectorMap fvm,
+                               const Kernel& )
+{
+  typename boost::graph_traits<TriangleMesh>::halfedge_descriptor h;
+  typename boost::graph_traits<TriangleMesh>::face_descriptor f;
+  typename Kernel::Vector_3 sum(0., 0., 0.), n;
+
+  CGAL::Halfedge_around_target_circulator<TriangleMesh> hedgeb(halfedge(v,tm),tm), hedgee = hedgeb;
+
+  do
+    {
+      h = *hedgeb;
+      if (is_border_edge(h,tm))
+	{
+	  hedgeb++;
+	  continue;
+	}
+
+      f =  face(h,tm);
+      n = get(fvm,f);
+      sum = (sum + n);
+      hedgeb++;
+    }
+  while (hedgeb != hedgee);
+  sum = sum / std::sqrt(sum * sum);
+  return sum;
+}
+
+template <typename TriangleMesh, typename FaceVectorMap,typename Kernel>
+void compute_facets_normals(const TriangleMesh& tm,
+                            FaceVectorMap fvm,
+                            const Kernel& k)
+{
+  typedef typename boost::property_traits<FaceVectorMap>::value_type Vector_3;
+
+  typedef typename boost::property_map<TriangleMesh,CGAL::vertex_point_t>::const_type VPM;
+  VPM vpm = get(CGAL::vertex_point,tm);
+  BOOST_FOREACH(typename boost::graph_traits<TriangleMesh>::face_descriptor f, faces(tm)){
+    typename boost::graph_traits<TriangleMesh>::halfedge_descriptor h = halfedge(f,tm);
+    Vector_3 normal = 
+      CGAL::cross_product(get(vpm, target(h,tm)) -
+			    get(vpm, target(opposite(h,tm),tm)),
+                          get(vpm, target(next(h,tm),tm)) -
+			    get(vpm, target(opposite(h,tm),tm)));
+      put(fvm, f, normal / CGAL::sqrt(normal * normal));
+  }
+}
+
+
+
+#endif

Ridges_3/test/Ridges_3/ridge_test.cpp

@@ -2,92 +2,64 @@
 #include <cassert>
 #include <fstream>
 #include <vector>
-#include <CGAL/Ridges.h> 
-#include <CGAL/Umbilics.h>
-#include <CGAL/Monge_via_jet_fitting.h> 
-
 //This Is an enriched Polyhedron with facets' normal
 #include "PolyhedralSurf.h"
 #include "PolyhedralSurf_rings.h"
+#include "compute_normals.h"
+#include <CGAL/Ridges.h> 
+#include <CGAL/Umbilics.h>
+#include <CGAL/Monge_via_jet_fitting.h> 
 
 // Functions declared in PolyhedralSurf.h
 // They were previously defined in a separate file PolyhedralSurf.cpp, 
 // but I prefere to avoid custom CMakeLists.txt files in the testsuite.
 // -- Laurent Rineau, 2008/11/10
-void PolyhedralSurf::compute_facets_normals()
-{
-  std::for_each(this->facets_begin(), this->facets_end(),
-		Facet_unit_normal()); 
-}
-
-const Vector_3 PolyhedralSurf::computeFacetsAverageUnitNormal(const Vertex_const_handle v)
-{
-  Halfedge_const_handle h;
-  Facet_const_handle f;
-  Vector_3 sum(0., 0., 0.), n;
-
-  Halfedge_around_vertex_const_circulator
-    hedgeb = v->vertex_begin(), hedgee = hedgeb;
-
-  do
-    {
-      h = hedgeb;
-      if (h->is_border_edge())
-	{
-	  hedgeb++;
-	  continue;
-	}
-
-      f =  h->facet();
-      n = f->getUnitNormal();
-      sum = (sum + n);
-      hedgeb++;
-    }
-  while (hedgeb != hedgee);
-  sum = sum / std::sqrt(sum * sum);
-  return sum;
-}
-
 
 typedef PolyhedralSurf::Traits          Kernel;
 typedef Kernel::FT                      FT;
 typedef Kernel::Point_3                 Point_3;
 typedef Kernel::Vector_3                Vector_3;
 
-typedef PolyhedralSurf::Vertex_const_handle   Vertex_const_handle;
+typedef boost::graph_traits<PolyhedralSurf>::vertex_descriptor  vertex_descriptor;
-typedef PolyhedralSurf::Vertex_const_iterator Vertex_const_iterator;
+typedef boost::graph_traits<PolyhedralSurf>::vertex_iterator vertex_iterator;
+typedef boost::graph_traits<PolyhedralSurf>::face_descriptor face_descriptor;
 
 typedef T_PolyhedralSurf_rings<PolyhedralSurf> Poly_rings;
 typedef CGAL::Monge_via_jet_fitting<Kernel>    Monge_via_jet_fitting;
 typedef Monge_via_jet_fitting::Monge_form      Monge_form;
 
-typedef CGAL::Vertex2Data_Property_Map_with_std_map<PolyhedralSurf> Vertex2Data_Property_Map_with_std_map;
-typedef Vertex2Data_Property_Map_with_std_map::Vertex2FT_map Vertex2FT_map;
-typedef Vertex2Data_Property_Map_with_std_map::Vertex2Vector_map Vertex2Vector_map;
-typedef Vertex2Data_Property_Map_with_std_map::Vertex2FT_property_map Vertex2FT_property_map;
-typedef Vertex2Data_Property_Map_with_std_map::Vertex2Vector_property_map Vertex2Vector_property_map;
 
+typedef std::map<vertex_descriptor, FT> VertexFT_map;
+typedef boost::associative_property_map< VertexFT_map > VertexFT_property_map;
+
+typedef std::map<vertex_descriptor, Vector_3> VertexVector_map;
+typedef boost::associative_property_map< VertexVector_map > VertexVector_property_map;
+
+typedef std::map<face_descriptor, Vector_3> Face2Vector_map;
+typedef boost::associative_property_map< Face2Vector_map > Face2Vector_property_map;
 //RIDGES
 typedef CGAL::Ridge_line<PolyhedralSurf> Ridge_line;
 typedef CGAL::Ridge_approximation < PolyhedralSurf,
-				    Vertex2FT_property_map,
+				    VertexFT_property_map,
-				    Vertex2Vector_property_map > Ridge_approximation;
+				    VertexVector_property_map > Ridge_approximation;
 //UMBILICS
 typedef CGAL::Umbilic<PolyhedralSurf> Umbilic;
 typedef CGAL::Umbilic_approximation < PolyhedralSurf,
-				      Vertex2FT_property_map, 
+				      VertexFT_property_map,
-				      Vertex2Vector_property_map > Umbilic_approximation;
+				      VertexVector_property_map > Umbilic_approximation;
 
 //create property maps, to be moved in main?
-Vertex2FT_map vertex2k1_map, vertex2k2_map, 
+VertexFT_map vertex_k1_map, vertex_k2_map,
-  vertex2b0_map, vertex2b3_map, 
+  vertex_b0_map, vertex_b3_map,
-  vertex2P1_map, vertex2P2_map;
+  vertex_P1_map, vertex_P2_map;
-Vertex2Vector_map vertex2d1_map, vertex2d2_map;
+VertexVector_map vertex_d1_map, vertex_d2_map;
+Face2Vector_map face2normal_map;
 
-Vertex2FT_property_map vertex2k1_pm(vertex2k1_map), vertex2k2_pm(vertex2k2_map), 
+VertexFT_property_map vertex_k1_pm(vertex_k1_map), vertex_k2_pm(vertex_k2_map),
-  vertex2b0_pm(vertex2b0_map), vertex2b3_pm(vertex2b3_map), 
+  vertex_b0_pm(vertex_b0_map), vertex_b3_pm(vertex_b3_map),
-  vertex2P1_pm(vertex2P1_map), vertex2P2_pm(vertex2P2_map);
+  vertex_P1_pm(vertex_P1_map), vertex_P2_pm(vertex_P2_map);
-Vertex2Vector_property_map vertex2d1_pm(vertex2d1_map), vertex2d2_pm(vertex2d2_map);
+VertexVector_property_map vertex_d1_pm(vertex_d1_map), vertex_d2_pm(vertex_d2_map);
+Face2Vector_property_map  face2normal_pm(face2normal_map);
 
   // default fct parameter values and global variables
   unsigned int d_fitting = 4;
@@ -99,18 +71,22 @@ Vertex2Vector_property_map vertex2d1_pm(vertex2d1_map), vertex2d2_pm(vertex2d2_m
   bool verbose = false;
   unsigned int min_nb_points = (d_fitting + 1) * (d_fitting + 2) / 2;
 
+
+
 /* gather points around the vertex v using rings on the
    polyhedralsurf. the collection of points resorts to 3 alternatives:
    1. the exact number of points to be used
    2. the exact number of rings to be used
    3. nothing is specified
 */
-void gather_fitting_points(Vertex_const_handle v, 
+template <typename VertexPointMap>
+void gather_fitting_points(vertex_descriptor v,
 			   std::vector<Point_3> &in_points,
-			   Poly_rings& poly_rings)
+			   Poly_rings& poly_rings,
+                           VertexPointMap vpm)
 {
   //container to collect vertices of v on the PolyhedralSurf
-  std::vector<Vertex_const_handle> gathered; 
+  std::vector<vertex_descriptor> gathered;
   //initialize
   in_points.clear();
 
@@ -118,7 +94,7 @@ void gather_fitting_points(Vertex_const_handle v,
   //enough rings and discard some points of the last collected ring to
   //get the exact "nb_points_to_use"
   if ( nb_points_to_use != 0 ) {
-    poly_rings.collect_enough_rings(v, nb_points_to_use, gathered);//, vpm);
+    poly_rings.collect_enough_rings(v, nb_points_to_use, gathered);
     if ( gathered.size() > nb_points_to_use ) gathered.resize(nb_points_to_use);
   }
   else { // nb_points_to_use=0, this is the default and the option -p is not considered;
@@ -126,14 +102,14 @@ void gather_fitting_points(Vertex_const_handle v,
     // enough rings to get the min_nb_points required for the fitting
     // else collect the nb_rings required
     if ( nb_rings == 0 )
-      poly_rings.collect_enough_rings(v, min_nb_points, gathered);//, vpm);
+      poly_rings.collect_enough_rings(v, min_nb_points, gathered);
-    else poly_rings.collect_i_rings(v, nb_rings, gathered);//, vpm);
+    else poly_rings.collect_i_rings(v, nb_rings, gathered);
   }
 
   //store the gathered points
-  std::vector<Vertex_const_handle>::const_iterator 
+  std::vector<vertex_descriptor>::const_iterator
     itb = gathered.begin(), ite = gathered.end();
-  CGAL_For_all(itb,ite) in_points.push_back((*itb)->point());
+  CGAL_For_all(itb,ite) in_points.push_back(get(vpm,*itb));
 }
 
 /* Use the jet_fitting package and the class Poly_rings to compute
@@ -144,17 +120,20 @@ void compute_differential_quantities(PolyhedralSurf& P, Poly_rings& poly_rings)
   //container for approximation points
   std::vector<Point_3> in_points;
 
+  typedef boost::property_map<PolyhedralSurf,CGAL::vertex_point_t>::type VPM;
+  VPM vpm = get(CGAL::vertex_point,P);
+
   //MAIN LOOP
-  Vertex_const_iterator vitb = P.vertices_begin(), vite = P.vertices_end();
+  vertex_iterator vitb = P.vertices_begin(), vite = P.vertices_end();
   for (; vitb != vite; vitb++) {
     //initialize
-    Vertex_const_handle v = vitb;
+    vertex_descriptor v = * vitb;
     in_points.clear();
     Monge_form monge_form;
     Monge_via_jet_fitting monge_fit;
 
     //gather points around the vertex using rings
-    gather_fitting_points(v, in_points, poly_rings);
+    gather_fitting_points(v, in_points, poly_rings, vpm);
 
     //exit if the nb of points is too small
     if ( in_points.size() < min_nb_points )
@@ -167,26 +146,26 @@ void compute_differential_quantities(PolyhedralSurf& P, Poly_rings& poly_rings)
 			   d_fitting, d_monge);
 
     //switch min-max ppal curv/dir wrt the mesh orientation
-    const Vector_3 normal_mesh = P.computeFacetsAverageUnitNormal(v);
+    const Vector_3 normal_mesh = computeFacetsAverageUnitNormal(P,v, face2normal_pm, Kernel());
     monge_form.comply_wrt_given_normal(normal_mesh);
 
     //Store monge data needed for ridge computations in property maps
-    vertex2d1_map[v] = monge_form.maximal_principal_direction();
+    vertex_d1_map[v] = monge_form.maximal_principal_direction();
-    vertex2d2_map[v] = monge_form.minimal_principal_direction();
+    vertex_d2_map[v] = monge_form.minimal_principal_direction();
-    vertex2k1_map[v] = monge_form.coefficients()[0];
+    vertex_k1_map[v] = monge_form.coefficients()[0];
-    vertex2k2_map[v] = monge_form.coefficients()[1];
+    vertex_k2_map[v] = monge_form.coefficients()[1];
-    vertex2b0_map[v] = monge_form.coefficients()[2];
+    vertex_b0_map[v] = monge_form.coefficients()[2];
-    vertex2b3_map[v] = monge_form.coefficients()[5];
+    vertex_b3_map[v] = monge_form.coefficients()[5];
     if ( d_monge >= 4) {
       //= 3*b1^2+(k1-k2)(c0-3k1^3)
-      vertex2P1_map[v] =
+      vertex_P1_map[v] =
 	3*monge_form.coefficients()[3]*monge_form.coefficients()[3]
 	+(monge_form.coefficients()[0]-monge_form.coefficients()[1])
 	*(monge_form.coefficients()[6]
 	  -3*monge_form.coefficients()[0]*monge_form.coefficients()[0]
 	  *monge_form.coefficients()[0]);
       //= 3*b2^2+(k2-k1)(c4-3k2^3)
-      vertex2P2_map[v] = 
+      vertex_P2_map[v] =
 	3*monge_form.coefficients()[4]*monge_form.coefficients()[4]
 	+(-monge_form.coefficients()[0]+monge_form.coefficients()[1])
 	*(monge_form.coefficients()[10]
@@ -210,7 +189,7 @@ int main()
     {std::cerr << "not enough points in the model" << std::endl;   return 1;}
 
   //initialize Polyhedral data : normal of facets
-  P.compute_facets_normals();
+  compute_facets_normals(P,face2normal_pm, Kernel());
 
   //create a Poly_rings object
   Poly_rings poly_rings(P);
@@ -223,28 +202,23 @@ int main()
   //---------------------------------------------------------------------------
   //Ridges
   //--------------------------------------------------------------------------
-  Ridge_approximation ridge_approximation_tag_3(P, 
+   Ridge_approximation ridge_approximation(P,
-						vertex2k1_pm, vertex2k2_pm,
+					  vertex_k1_pm, vertex_k2_pm,
-						vertex2b0_pm, vertex2b3_pm,
+					  vertex_b0_pm, vertex_b3_pm,
-						vertex2d1_pm, vertex2d2_pm,
+					  vertex_d1_pm, vertex_d2_pm,
-						Vertex2FT_property_map(), 
+					  vertex_P1_pm, vertex_P2_pm );
-						Vertex2FT_property_map());
   std::vector<Ridge_line*> ridge_lines;
-  back_insert_iterator<std::vector<Ridge_line*> > ii(ridge_lines);
+  std::back_insert_iterator<std::vector<Ridge_line*> > ii(ridge_lines);
   
   //Find MAX_RIDGE, RED_RIDGE, CREST or all ridges
-  ridge_approximation_tag_3.compute_max_ridges(ii, tag_order);  
+  ridge_approximation.compute_max_ridges(ii, tag_order);  
-  ridge_approximation_tag_3.compute_min_ridges(ii, tag_order);  
+  ridge_approximation.compute_min_ridges(ii, tag_order);  
-  ridge_approximation_tag_3.compute_crest_ridges(ii, tag_order);  
+  ridge_approximation.compute_crest_ridges(ii, tag_order);  
  
   std::cout << "Compute ridges with tag_4" << std::endl;
   tag_order =  CGAL::Ridge_order_4;
    //Find MAX_RIDGE, RED_RIDGE, CREST or all ridges
-  Ridge_approximation ridge_approximation(P, 
+ 
-					  vertex2k1_pm, vertex2k2_pm,
-					  vertex2b0_pm, vertex2b3_pm,
-					  vertex2d1_pm, vertex2d2_pm,
-					  vertex2P1_pm, vertex2P2_pm );
   ridge_approximation.compute_max_ridges(ii, tag_order);  
   ridge_approximation.compute_min_ridges(ii, tag_order);  
   ridge_approximation.compute_crest_ridges(ii, tag_order);  
@@ -257,10 +231,10 @@ int main()
   // UMBILICS
   //--------------------------------------------------------------------------
   Umbilic_approximation umbilic_approximation(P, 
-					      vertex2k1_pm, vertex2k2_pm,
+					      vertex_k1_pm, vertex_k2_pm,
-					      vertex2d1_pm, vertex2d2_pm);
+					      vertex_d1_pm, vertex_d2_pm);
   std::vector<Umbilic*> umbilics;
-  back_insert_iterator<std::vector<Umbilic*> > umb_it(umbilics);
+  std::back_insert_iterator<std::vector<Umbilic*> > umb_it(umbilics);
   std::cout << "compute umbilics u=1" << std::endl;
   umbilic_approximation.compute(umb_it, umb_size);
   umb_size=2;