diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/Weights.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/Weights.h index 0732931dd5d..fec0c593f9d 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/Weights.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/Weights.h @@ -87,13 +87,13 @@ protected: typedef typename Kernel_traits::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()(v0,v1,v2,ppmap); + return Cotangent_value_Meyer_impl()(v0, v1, v2, ppmap()); } }; @@ -120,13 +125,13 @@ class Cotangent_value_Meyer_secure typedef typename Kernel_traits::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::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::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::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::vertex_descriptor vertex_descriptor; typedef typename boost::graph_traits::in_edge_iterator in_edge_iterator; typedef typename boost::graph_traits::halfedge_descriptor halfedge_descriptor; @@ -351,6 +376,11 @@ public: return CotangentValue::pmesh(); } + VertexPointMap& ppmap() + { + return CotangentValue::ppmap(); + } + typedef typename boost::graph_traits::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::halfedge_descriptor halfedge_descriptor; typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; @@ -506,6 +541,11 @@ public: return CotangentValue::pmesh(); } + VertexPointMap& ppmap() + { + return CotangentValue::ppmap(); + } + typedef typename boost::graph_traits::halfedge_descriptor halfedge_descriptor; typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; @@ -527,6 +567,82 @@ public: } }; + +template::type + , class CotangentValue = Cotangent_value_Meyer +> +class Cotangent_weight_with_triangle_area : CotangentValue +{ + typedef PolygonMesh PM; + typedef VertexPointMap VPMap; + typedef typename boost::property_traits::value_type Point; + + typedef typename boost::graph_traits::halfedge_descriptor halfedge_descriptor; + typedef typename boost::graph_traits::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::type> diff --git a/Surface_mesh_parameterization/include/CGAL/Surface_mesh_parameterization/Orbital_Tutte_parameterizer_3.h b/Surface_mesh_parameterization/include/CGAL/Surface_mesh_parameterization/Orbital_Tutte_parameterizer_3.h index fe8f77174b0..17c3eab166f 100644 --- a/Surface_mesh_parameterization/include/CGAL/Surface_mesh_parameterization/Orbital_Tutte_parameterizer_3.h +++ b/Surface_mesh_parameterization/include/CGAL/Surface_mesh_parameterization/Orbital_Tutte_parameterizer_3.h @@ -28,6 +28,8 @@ #include +#include + #include #include #include @@ -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 + 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 Cotan_weights; +// typedef CGAL::internal::Cotangent_weight_with_triangle_area 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 void assign_solution(const TriangleMesh& mesh, @@ -734,7 +782,10 @@ public: // %%%%%%%%%%%%%%%%%%%% Matrix L(2 * nbVertices, 2 * nbVertices); - mean_value_laplacian(mesh, vimap, L); + 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 std::ofstream outL("matrices/L.txt"); @@ -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