Added cotangent weights for Orbifold embeddings

Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/Weights.h

@@ -87,13 +87,13 @@ protected:
   typedef typename Kernel_traits<Point>::Kernel::Vector_3  Vector;
 
   PolygonMesh& pmesh_;
-  Point_property_map ppmap;
+  Point_property_map ppmap_;
 
 public:
 
   Cotangent_value_Meyer(PolygonMesh& pmesh_, VertexPointMap vpmap_)
     : pmesh_(pmesh_)
-    , ppmap(vpmap_)
+    , ppmap_(vpmap_)
   {}
 
   PolygonMesh& pmesh()
@@ -101,9 +101,14 @@ public:
     return pmesh_;
   }
 
+  Point_property_map& ppmap()
+  {
+    return ppmap_;
+  }
+
   double operator()(vertex_descriptor v0, vertex_descriptor v1, vertex_descriptor v2)
   {
-    return Cotangent_value_Meyer_impl<PolygonMesh>()(v0,v1,v2,ppmap);
+    return Cotangent_value_Meyer_impl<PolygonMesh>()(v0, v1, v2, ppmap());
   }
 };
 
@@ -120,13 +125,13 @@ class Cotangent_value_Meyer_secure
   typedef typename Kernel_traits<Point>::Kernel::Vector_3  Vector;
 
   PolygonMesh& pmesh_;
-  Point_property_map ppmap;
+  Point_property_map ppmap_;
 
 public:
 
   Cotangent_value_Meyer_secure(PolygonMesh& pmesh_, VertexPointMap vpmap_)
     : pmesh_(pmesh_)
-    , ppmap(vpmap_)
+    , ppmap_(vpmap_)
   {}
 
   PolygonMesh& pmesh()
@@ -171,6 +176,11 @@ public:
     return CotangentValue::pmesh();
   }
 
+  VertexPointMap& ppmap()
+  {
+    return CotangentValue::ppmap();
+  }
+
   typedef typename boost::graph_traits<PolygonMesh>::vertex_descriptor vertex_descriptor;
 
   double operator()(vertex_descriptor v0, vertex_descriptor v1, vertex_descriptor v2)
@@ -201,6 +211,11 @@ public:
     return CotangentValue::pmesh();
   }
 
+  VertexPointMap& ppmap()
+  {
+    return CotangentValue::ppmap();
+  }
+
   typedef typename boost::graph_traits<PolygonMesh>::vertex_descriptor vertex_descriptor;
 
   double operator()(vertex_descriptor v0, vertex_descriptor v1, vertex_descriptor v2)
@@ -229,6 +244,11 @@ public:
     return CotangentValue::pmesh();
   }
 
+  VertexPointMap& ppmap()
+  {
+    return CotangentValue::ppmap();
+  }
+
   typedef typename boost::graph_traits<PolygonMesh>::vertex_descriptor vertex_descriptor;
 
   double operator()(vertex_descriptor v0, vertex_descriptor v1, vertex_descriptor v2)
@@ -275,6 +295,11 @@ public:
     return CotangentValue::pmesh();
   }
 
+  VertexPointMap& ppmap()
+  {
+    return CotangentValue::ppmap();
+  }
+
   typedef typename boost::graph_traits<PolygonMesh>::vertex_descriptor vertex_descriptor;
   typedef typename boost::graph_traits<PolygonMesh>::in_edge_iterator in_edge_iterator;
   typedef typename boost::graph_traits<PolygonMesh>::halfedge_descriptor halfedge_descriptor;
@@ -351,6 +376,11 @@ public:
     return CotangentValue::pmesh();
   }
 
+  VertexPointMap& ppmap()
+  {
+    return CotangentValue::ppmap();
+  }
+
   typedef typename boost::graph_traits<PolygonMesh>::vertex_descriptor vertex_descriptor;
 
   double operator()(vertex_descriptor v0, vertex_descriptor v1, vertex_descriptor v2)
@@ -428,6 +458,11 @@ public:
     return CotangentValue::pmesh();
   }
 
+  VertexPointMap& ppmap()
+  {
+    return CotangentValue::ppmap();
+  }
+
   typedef typename boost::graph_traits<PolygonMesh>::halfedge_descriptor   halfedge_descriptor;
   typedef typename boost::graph_traits<PolygonMesh>::vertex_descriptor vertex_descriptor;
 
@@ -506,6 +541,11 @@ public:
     return CotangentValue::pmesh();
   }
 
+  VertexPointMap& ppmap()
+  {
+    return CotangentValue::ppmap();
+  }
+
   typedef typename boost::graph_traits<PolygonMesh>::halfedge_descriptor   halfedge_descriptor;
   typedef typename boost::graph_traits<PolygonMesh>::vertex_descriptor vertex_descriptor;
 
@@ -527,6 +567,82 @@ public:
   }
 };
 
+
+template<class PolygonMesh
+  , class VertexPointMap = typename boost::property_map<PolygonMesh, vertex_point_t>::type
+  , class CotangentValue = Cotangent_value_Meyer<PolygonMesh, VertexPointMap>
+>
+class Cotangent_weight_with_triangle_area : CotangentValue
+{
+  typedef PolygonMesh                                         PM;
+  typedef VertexPointMap                                      VPMap;
+  typedef typename boost::property_traits<VPMap>::value_type  Point;
+
+  typedef typename boost::graph_traits<PM>::halfedge_descriptor   halfedge_descriptor;
+  typedef typename boost::graph_traits<PM>::vertex_descriptor     vertex_descriptor;
+
+  Cotangent_weight_with_triangle_area()
+  {}
+public:
+  Cotangent_weight_with_triangle_area(PolygonMesh& pmesh_, VertexPointMap vpmap_)
+    : CotangentValue(pmesh_, vpmap_)
+  {}
+
+  PolygonMesh& pmesh()
+  {
+    return CotangentValue::pmesh();
+  }
+
+  VertexPointMap& ppmap()
+  {
+    return CotangentValue::ppmap();
+  }
+
+  double operator()(halfedge_descriptor he)
+  {
+    vertex_descriptor v0 = target(he, pmesh());
+    vertex_descriptor v1 = source(he, pmesh());
+
+     // Only one triangle for border edges
+    if (is_border_edge(he, pmesh()))
+    {
+      halfedge_descriptor he_cw = opposite( next(he, pmesh()) , pmesh() );
+      vertex_descriptor v2 = source(he_cw, pmesh());
+      if (is_border_edge(he_cw, pmesh()) )
+      {
+        halfedge_descriptor he_ccw = prev( opposite(he, pmesh()) , pmesh() );
+        v2 = source(he_ccw, pmesh());
+      }
+
+      const Point& v0_p = get(ppmap(), v0);
+      const Point& v1_p = get(ppmap(), v1);
+      const Point& v2_p = get(ppmap(), v2);
+      double area_t = to_double(CGAL::sqrt(squared_area(v0_p, v1_p, v2_p)));
+
+      return ( CotangentValue::operator()(v0, v2, v1) / area_t );
+    }
+    else
+    {
+      halfedge_descriptor he_cw = opposite( next(he, pmesh()) , pmesh() );
+      vertex_descriptor v2 = source(he_cw, pmesh());
+      halfedge_descriptor he_ccw = prev( opposite(he, pmesh()) , pmesh() );
+      vertex_descriptor v3 = source(he_ccw, pmesh());
+
+      const Point& v0_p = get(ppmap(), v0);
+      const Point& v1_p = get(ppmap(), v1);
+      const Point& v2_p = get(ppmap(), v2);
+      const Point& v3_p = get(ppmap(), v3);
+      double area_t1 = to_double(CGAL::sqrt(squared_area(v0_p, v1_p, v2_p)));
+      double area_t2 = to_double(CGAL::sqrt(squared_area(v0_p, v1_p, v3_p)));
+
+      return ( CotangentValue::operator()(v0, v2, v1) / area_t1
+             + CotangentValue::operator()(v0, v3, v1) / area_t2);
+     }
+
+    return 0.;
+  }
+};
+
 // Mean value calculator described in -[Floater04] Mean Value Coordinates-
 template<class PolygonMesh
        , class VertexPointMap = typename boost::property_map<PolygonMesh, CGAL::vertex_point_t>::type>

Surface_mesh_parameterization/include/CGAL/Surface_mesh_parameterization/Orbital_Tutte_parameterizer_3.h

@@ -28,6 +28,8 @@
 
 #include <CGAL/Surface_mesh_parameterization/Error_code.h>
 
+#include <CGAL/Polygon_mesh_processing/Weights.h>
+
 #include <CGAL/assertions.h>
 #include <CGAL/circulator.h>
 #include <CGAL/Eigen_solver_traits.h>
@@ -74,6 +76,12 @@ enum Orbifold_type
   Parallelogram
 };
 
+enum Weight_type
+{
+  Cotangent = 0,
+  Mean_value
+};
+
 const char* get_orbifold_type(int orb_type)
 {
   // Messages corresponding to Error_code list above. Must be kept in sync!
@@ -90,12 +98,6 @@ const char* get_orbifold_type(int orb_type)
     return type[orb_type];
 }
 
-enum Weight_type
-{
-  Cotan,
-  Mvc
-};
-
 template
 <
   typename TriangleMesh,
@@ -131,6 +133,7 @@ private:
   typedef typename Kernel::Point_3                                  Point_3;
 
   Orbifold_type orb_type;
+  Weight_type weight_type;
 
 private:
   // check input
@@ -534,6 +537,51 @@ private:
     }
   }
 
+  template<typename VertexIndexMap>
+  void cotangent_laplacien(TriangleMesh& mesh,
+                           VertexIndexMap vimap,
+                           Matrix& L) const
+  {
+    const PPM ppmap = get(vertex_point, mesh);
+
+    // not exactly sure which cotan weights should be used:
+    // 0.5 (cot a + cot b) ? 1/T1 cot a + 1/T2 cot b ? 1/Vor(i) (cot a + cot b?)
+    // Comparing to the matlab code, the basic Cotangent_weight gives the same results.
+    typedef CGAL::internal::Cotangent_weight<TriangleMesh>                      Cotan_weights;
+//    typedef CGAL::internal::Cotangent_weight_with_triangle_area<TriangleMesh>   Cotan_weights;
+
+    Cotan_weights cotan_weight_calculator(mesh, ppmap);
+
+    BOOST_FOREACH(halfedge_descriptor hd, halfedges(mesh)) {
+      vertex_descriptor vi = source(hd, mesh);
+      vertex_descriptor vj = target(hd, mesh);
+      int i = get(vimap, vi);
+      int j = get(vimap, vj);
+
+      if(i > j)
+        continue;
+
+      // times 2 because Cotangent_weight returns 1/2 (cot alpha + cot beta)...
+      double w_ij = 2 * cotan_weight_calculator(hd);
+
+      // ij
+      L.set_coef(2*i, 2*j, w_ij, true /* new coef */);
+      L.set_coef(2*i +1, 2*j + 1, w_ij, true /* new coef */);
+
+      // ji
+      L.set_coef(2*j, 2*i, w_ij, true /* new coef */);
+      L.set_coef(2*j +1, 2*i + 1, w_ij, true /* new coef */);
+
+      // ii
+      L.add_coef(2*i, 2*i, - w_ij);
+      L.add_coef(2*i + 1, 2*i + 1, - w_ij);
+
+      // jj
+      L.add_coef(2*j, 2*j, - w_ij);
+      L.add_coef(2*j + 1, 2*j + 1, - w_ij);
+    }
+  }
+
   /// Copy the solution into the UV property map.
   template <typename VertexIndexMap, typename VertexUVMap>
   void assign_solution(const TriangleMesh& mesh,
@@ -734,6 +782,9 @@ public:
     // %%%%%%%%%%%%%%%%%%%%
 
     Matrix L(2 * nbVertices, 2 * nbVertices);
+    if(weight_type == Cotangent)
+      cotangent_laplacien(mesh, vimap, L);
+    else // weight_type == Mean_value
       mean_value_laplacian(mesh, vimap, L);
 
 #ifdef CGAL_SMP_OUTPUT_ORBITAL_MATRICES
@@ -769,7 +820,12 @@ public:
 
 /// Constructor
 public:
-  Orbital_Tutte_parameterizer_3(Orbifold_type orb_type) : orb_type(orb_type) { }
+  Orbital_Tutte_parameterizer_3(Orbifold_type orb_type = Square,
+                                Weight_type weight_type = Cotangent)
+    :
+      orb_type(orb_type),
+      weight_type(weight_type)
+  { }
 };
 
 } // namespace Surface_mesh_parameterization