diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Hole_filling/Weights.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Hole_filling/Weights.h index e0579dfc1fc..69aaa7837d3 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Hole_filling/Weights.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Hole_filling/Weights.h @@ -63,6 +63,44 @@ public: // The potential problem with previous one (Cotangent_value) is that it does not produce symmetric results // (i.e. for v0, v1, v2 and v2, v1, v0 returned cot weights can be slightly different) // This one provides stable results. +template +struct Cotangent_value_Meyer_impl +{ + typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; + + template + double operator()(vertex_descriptor v0, vertex_descriptor v1, vertex_descriptor v2, const VertexPointMap& ppmap) + { + typedef typename Kernel_traits< + typename boost::property_traits::value_type >::Kernel::Vector_3 Vector; + + Vector a = get(ppmap, v0) - get(ppmap, v1); + Vector b = get(ppmap, v2) - get(ppmap, v1); + + double dot_ab = a*b; + // rewritten for safer fp arithmetic + //double dot_aa = a.squared_length(); + //double dot_bb = b.squared_length(); + //double divider = CGAL::sqrt( dot_aa * dot_bb - dot_ab * dot_ab ); + + Vector cross_ab = CGAL::cross_product(a, b); + double divider = CGAL::sqrt(cross_ab*cross_ab); + + if(divider == 0 /*|| divider != divider*/) + { + CGAL::collinear(get(ppmap, v0), get(ppmap, v1), get(ppmap, v2)) ? + CGAL_warning(!"Infinite Cotangent value with degenerate triangle!") : + CGAL_warning(!"Infinite Cotangent value due to floating point arithmetic!"); + + + return dot_ab > 0 ? (std::numeric_limits::max)() : + -(std::numeric_limits::max)(); + } + + return dot_ab / divider; + } +}; + template::type> class Cotangent_value_Meyer @@ -90,30 +128,7 @@ public: double operator()(vertex_descriptor v0, vertex_descriptor v1, vertex_descriptor v2) { - Vector a = get(ppmap, v0) - get(ppmap, v1); - Vector b = get(ppmap, v2) - get(ppmap, v1); - - double dot_ab = a*b; - // rewritten for safer fp arithmetic - //double dot_aa = a.squared_length(); - //double dot_bb = b.squared_length(); - //double divider = CGAL::sqrt( dot_aa * dot_bb - dot_ab * dot_ab ); - - Vector cross_ab = CGAL::cross_product(a, b); - double divider = CGAL::sqrt(cross_ab*cross_ab); - - if(divider == 0 /*|| divider != divider*/) - { - CGAL::collinear(get(ppmap, v0), get(ppmap, v1), get(ppmap, v2)) ? - CGAL_warning(!"Infinite Cotangent value with degenerate triangle!") : - CGAL_warning(!"Infinite Cotangent value due to floating point arithmetic!"); - - - return dot_ab > 0 ? (std::numeric_limits::max)() : - -(std::numeric_limits::max)(); - } - - return dot_ab / divider; + return Cotangent_value_Meyer_impl()(v0,v1,v2,ppmap); } }; @@ -248,6 +263,20 @@ public: } }; +template > +struct Cotangent_value_minimum_zero_impl : CotangentValue +{ + typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; + + template + double operator()(vertex_descriptor v0, vertex_descriptor v1, vertex_descriptor v2, const VertexPointMap& ppmap) + { + double value = CotangentValue::operator()(v0, v1, v2,ppmap); + return (std::max)(0.0, value); + } +}; + template::type , class CotangentValue @@ -362,6 +391,46 @@ public: // Cotangent weight calculator // Cotangent_value: as suggested by -[Sorkine07] ARAP Surface Modeling- // Cotangent_value_area_weighted: as suggested by -[Mullen08] Spectral Conformal Parameterization- +template< class PolygonMesh, + class CotangentValue = Cotangent_value_minimum_zero_impl > +struct Cotangent_weight_impl : CotangentValue +{ + typedef typename boost::graph_traits::halfedge_descriptor halfedge_descriptor; + typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; + + // Returns the cotangent weight of specified halfedge_descriptor + // Edge orientation is trivial + template + double operator()(halfedge_descriptor he, PolygonMesh& pmesh, const VertexPointMap& ppmap) + { + 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); + } + return ( CotangentValue::operator()(v0, v2, v1, ppmap)/2.0 ); + } + 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); + + return ( CotangentValue::operator()(v0, v2, v1, ppmap)/2.0 + + CotangentValue::operator()(v0, v3, v1, ppmap)/2.0 ); + } + } +}; + template::type , class CotangentValue @@ -424,6 +493,28 @@ public: }; // Single cotangent from -[Chao10] Simple Geometric Model for Elastic Deformation +template > +struct Single_cotangent_weight_impl : CotangentValue +{ + typedef typename boost::graph_traits::halfedge_descriptor halfedge_descriptor; + typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; + + // Returns the cotangent of the opposite angle of the edge + // 0 for border edges (which does not have an opposite angle) + template + double operator()(halfedge_descriptor he, PolygonMesh& pmesh, const VertexPointMap& ppmap) + { + if(is_border(he, pmesh)) { return 0.0;} + + vertex_descriptor v0 = target(he, pmesh); + vertex_descriptor v1 = source(he, pmesh); + + vertex_descriptor v_op = target(next(he, pmesh), pmesh); + return CotangentValue::operator()(v0, v_op, v1, ppmap); + } +}; + template::type , class CotangentValue = Cotangent_value_Meyer > diff --git a/Surface_modeling/examples/Surface_modeling/all_roi_assign_example_custom_polyhedron.cpp b/Surface_modeling/examples/Surface_modeling/all_roi_assign_example_custom_polyhedron.cpp index 562b996ff2f..2160083ab9f 100644 --- a/Surface_modeling/examples/Surface_modeling/all_roi_assign_example_custom_polyhedron.cpp +++ b/Surface_modeling/examples/Surface_modeling/all_roi_assign_example_custom_polyhedron.cpp @@ -2,6 +2,7 @@ #include #include #include +#include struct Custom_point_3{ // Required by File_scanner_OFF diff --git a/Surface_modeling/include/CGAL/Surface_mesh_deformation.h b/Surface_modeling/include/CGAL/Surface_mesh_deformation.h index 64c505d56d8..c83f6b552b3 100644 --- a/Surface_modeling/include/CGAL/Surface_mesh_deformation.h +++ b/Surface_modeling/include/CGAL/Surface_mesh_deformation.h @@ -21,10 +21,12 @@ #define CGAL_SURFACE_MESH_DEFORMATION_H #include -#include #include #include +#include +#include + #include #include #include @@ -56,19 +58,41 @@ enum Deformation_algorithm_tag /// @cond CGAL_DOCUMENT_INTERNAL namespace internal { template -struct Weight_calculator_selector { - typedef Uniform_weight weight_calculator; -}; +struct Weight_calculator_selector; template struct Weight_calculator_selector { - typedef Single_cotangent_weight weight_calculator; + typedef Single_cotangent_weight_impl weight_calculator; }; template struct Weight_calculator_selector { - typedef Cotangent_weight weight_calculator; + typedef Cotangent_weight_impl weight_calculator; }; + +// property map that create a Simple_cartesian::Point_3 +// on the fly in order the deformation class to be used +// with points with minimal requirements +template +struct SC_on_the_fly_pmap: public Vertex_point_map{ + typedef boost::readable_property_map_tag category; + typedef CGAL::Simple_cartesian::Point_3 value_type; + typedef value_type reference; + typedef typename boost::property_traits::key_type key_type; + + SC_on_the_fly_pmap(Vertex_point_map base): + Vertex_point_map(base) {} + + friend value_type + get(const SC_on_the_fly_pmap map, key_type k) + { + typename boost::property_traits::reference base= + get(static_cast(map), k); + return value_type(base[0], base[1], base[2]); + } +}; + + }//namespace internal /// @endcond @@ -367,13 +391,14 @@ public: private: void init() { + typedef internal::SC_on_the_fly_pmap Wrapper; // compute halfedge weights halfedge_iterator eb, ee; hedge_weight.reserve(2*num_edges(m_halfedge_graph)); for(cpp11::tie(eb, ee) = halfedges(m_halfedge_graph); eb != ee; ++eb) { hedge_weight.push_back( - this->weight_calculator(*eb, m_halfedge_graph, vertex_point_map)); + this->weight_calculator(*eb, m_halfedge_graph, Wrapper(vertex_point_map))); } #ifdef CGAL_DEFORM_MESH_USE_EXPERIMENTAL_SR_ARAP m_sr_arap_alpha=2; @@ -761,6 +786,7 @@ public: */ void overwrite_initial_geometry() { + typedef internal::SC_on_the_fly_pmap Wrapper; if(roi.empty()) { return; } // no ROI to overwrite region_of_solution(); // the roi should be preprocessed since we are using original_position vec @@ -781,13 +807,13 @@ public: std::size_t id_e = id(he); if(is_weight_computed[id_e]) { continue; } - hedge_weight[id_e] = weight_calculator(he, m_halfedge_graph, vertex_point_map); + hedge_weight[id_e] = weight_calculator(he, m_halfedge_graph, Wrapper(vertex_point_map)); is_weight_computed[id_e] = true; halfedge_descriptor e_opp = opposite(he, m_halfedge_graph); std::size_t id_e_opp = id(e_opp); - hedge_weight[id_e_opp] = weight_calculator(e_opp, m_halfedge_graph, vertex_point_map); + hedge_weight[id_e_opp] = weight_calculator(e_opp, m_halfedge_graph, Wrapper(vertex_point_map)); is_weight_computed[id_e_opp] = true; } } diff --git a/Surface_modeling/include/CGAL/internal/Surface_modeling/Weights.h b/Surface_modeling/include/CGAL/internal/Surface_modeling/Weights.h deleted file mode 100644 index 406a099ec67..00000000000 --- a/Surface_modeling/include/CGAL/internal/Surface_modeling/Weights.h +++ /dev/null @@ -1,313 +0,0 @@ -// Copyright (c) 2014 GeometryFactory -// All rights reserved. -// -// This file is part of CGAL (www.cgal.org). -// You can redistribute it and/or modify it under the terms of the GNU -// General Public License as published by the Free Software Foundation, -// either version 3 of the License, or (at your option) any later version. -// -// Licensees holding a valid commercial license may use this file in -// accordance with the commercial license agreement provided with the software. -// -// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE -// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. -// -// $URL$ -// $Id$ -// -// Author(s) : Ilker O. Yaz - -#ifndef CGAL_SURFACE_MODELING_WEIGHTS_H -#define CGAL_SURFACE_MODELING_WEIGHTS_H -/// @cond CGAL_DOCUMENT_INTERNAL - -#include -#include - - -namespace CGAL { -namespace internal { - -namespace Surface_modeling{ -typedef CGAL::Simple_cartesian::Vector_3 Vector; -} //end of namespace Surface_modeling - -template -Surface_modeling::Vector to_vector(const Point& b, const Point& a) { - return Surface_modeling::Vector(a[0] - b[0], a[1] - b[1], a[2] - b[2]); -} - -// Returns the cotangent value of half angle v0 v1 v2 -// using formula in -[Meyer02] Discrete Differential-Geometry Operators for- page 19 -// The potential problem with previous one (Cotangent_value) is that it does not produce symmetric results -// (i.e. for v0, v1, v2 and v2, v1, v0 returned cot weights can be slightly different) -// This one provides stable results. -template -class Cotangent_value_Meyer -{ -public: - typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; - - typedef HalfedgeGraph Halfedge_graph; - - template - double operator()(vertex_descriptor v0, vertex_descriptor v1, vertex_descriptor v2, VertexPointMap vpm) - { - const Surface_modeling::Vector& a = to_vector(get(vpm, v1), get(vpm, v0)); - const Surface_modeling::Vector& b = to_vector(get(vpm, v1), get(vpm, v2)); - - double dot_ab = a*b; - Surface_modeling::Vector cross_ab = CGAL::cross_product(a, b); - double divider = std::sqrt(cross_ab*cross_ab); - - if(divider == 0 /*|| divider != divider*/) - { - this->collinear(get(vpm, v0), get(vpm, v1), get(vpm, v2)) ? - CGAL_warning(!"Infinite Cotangent value with degenerate triangle!") : - CGAL_warning(!"Infinite Cotangent value due to floating point arithmetic!"); - - return dot_ab > 0 ? (std::numeric_limits::max)() : - -(std::numeric_limits::max)(); - } - - return dot_ab / divider; - } - - /////////////////////////////////////////////////////////////////////////////////////// - // WARNING: this two functions are just used when cotangent weight turns out to be +-inf, - // just for raising a proper warning message (i.e nothing functional) - template - bool collinear(const Point&, const Point&, const Point&) { - return true; - } - template - bool collinear(const CGAL::Point_3& a, const CGAL::Point_3& b, const CGAL::Point_3& c) { - return CGAL::collinear(a, b, c); - } - /////////////////////////////////////////////////////////////////////////////////////// - -}; - -// Returns the cotangent value of half angle v0 v1 v2 by clamping between [1, 89] degrees -// as suggested by -[Friedel] Unconstrained Spherical Parameterization- -template > -class Cotangent_value_clamped : CotangentValue -{ -public: - typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; - - typedef HalfedgeGraph Halfedge_graph; - - template - double operator()(vertex_descriptor v0, vertex_descriptor v1, vertex_descriptor v2, VertexPointMap vpm) - { - const double cot_1 = 57.289962; - const double cot_89 = 0.017455; - double value = CotangentValue::operator()(v0, v1, v2, vpm); - return (std::max)(cot_89, (std::min)(value, cot_1)); - } -}; - -template > -class Cotangent_value_minimum_zero : CotangentValue -{ -public: - typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; - - template - double operator()(vertex_descriptor v0, vertex_descriptor v1, vertex_descriptor v2, VertexPointMap vpm) - { - double value = CotangentValue::operator()(v0, v1, v2, vpm); - return (std::max)(0.0, value); - } -}; - - -///////////////////////////// Halfedge Weight Calculators /////////////////////////////////// -// Cotangent weight calculator -// Cotangent_value: as suggested by -[Sorkine07] ARAP Surface Modeling- -// Cotangent_value_area_weighted: as suggested by -[Mullen08] Spectral Conformal Parameterization- -template > -class Cotangent_weight : CotangentValue -{ -public: - typedef typename boost::graph_traits::halfedge_descriptor halfedge_descriptor; - typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; - typedef HalfedgeGraph Halfedge_graph; - // Returns the cotangent weight of specified halfedge_descriptor - // Edge orientation is trivial - template - double operator()(halfedge_descriptor he, HalfedgeGraph& halfedge_graph, VertexPointMap vpm) - { - vertex_descriptor v0 = target(he, halfedge_graph); - vertex_descriptor v1 = source(he, halfedge_graph); - // Only one triangle for border edges - if ( is_border(he, halfedge_graph) || - is_border(opposite(he, halfedge_graph), halfedge_graph) ) - { - halfedge_descriptor he_cw = opposite( next(he, halfedge_graph), halfedge_graph ); - vertex_descriptor v2 = source(he_cw, halfedge_graph); - if ( is_border(he_cw, halfedge_graph) || - is_border(opposite(he_cw, halfedge_graph), halfedge_graph) ) - { - halfedge_descriptor he_ccw = prev( opposite(he, halfedge_graph), halfedge_graph ); - v2 = source(he_ccw, halfedge_graph); - } - return ( CotangentValue::operator()(v0, v2, v1, vpm)/2.0 ); - } - else - { - halfedge_descriptor he_cw = opposite( next(he, halfedge_graph), halfedge_graph ); - vertex_descriptor v2 = source(he_cw, halfedge_graph); - halfedge_descriptor he_ccw = prev( opposite(he, halfedge_graph), halfedge_graph ); - vertex_descriptor v3 = source(he_ccw, halfedge_graph); - - return ( CotangentValue::operator()(v0, v2, v1, vpm)/2.0 + CotangentValue::operator()(v0, v3, v1, vpm)/2.0 ); - } - } -}; - -// Single cotangent from -[Chao10] Simple Geometric Model for Elastic Deformation -template > -class Single_cotangent_weight : CotangentValue -{ -public: - typedef typename boost::graph_traits::halfedge_descriptor halfedge_descriptor; - typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; - typedef HalfedgeGraph Halfedge_graph; - - // Returns the cotangent of the opposite angle of the halfedge - // 0 for border edges (which does not have an opposite angle) - template - double operator()(halfedge_descriptor he, HalfedgeGraph& halfedge_graph, VertexPointMap vpm) - { - if(is_border(he, halfedge_graph)) { return 0.0;} - - vertex_descriptor v0 = target(he, halfedge_graph); - vertex_descriptor v1 = source(he, halfedge_graph); - - vertex_descriptor v_op = target(next(he, halfedge_graph), halfedge_graph); - return CotangentValue::operator()(v0, v_op, v1, vpm); - } -}; - -// Mean value calculator described in -[Floater04] Mean Value Coordinates- -// WARNING: Need to be updated to use point pmap -template -class Mean_value_weight -{ -public: - typedef typename boost::graph_traits::halfedge_descriptor halfedge_descriptor; - typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; - typedef HalfedgeGraph Halfedge_graph; - - // Returns the mean-value coordinate of specified halfedge_descriptor - // Returns different value for different halfedge orientation (which is a normal behaviour according to formula) - double operator()(halfedge_descriptor he, HalfedgeGraph& halfedge_graph) - { - vertex_descriptor v0 = target(he, halfedge_graph); - vertex_descriptor v1 = source(he, halfedge_graph); - Surface_modeling::Vector vec(v1->point(), v0->point()); - double norm = std::sqrt( vec.squared_length() ); - - // Only one triangle for border edges - if ( is_border(he, halfedge_graph) || - is_border( opposite(he, halfedge_graph), halfedge_graph) ) - { - halfedge_descriptor he_cw = opposite( next(he, halfedge_graph), halfedge_graph ); - vertex_descriptor v2 = source(he_cw, halfedge_graph); - if ( is_border(he_cw, halfedge_graph) || - is_border(opposite(he_cw, halfedge_graph), halfedge_graph) ) - { - halfedge_descriptor he_ccw = prev( opposite(he, halfedge_graph), halfedge_graph ); - v2 = source(he_ccw, halfedge_graph); - } - - return ( half_tan_value_2(v1, v0, v2)/norm); - } - else - { - halfedge_descriptor he_cw = opposite( next(he, halfedge_graph), halfedge_graph ); - vertex_descriptor v2 = source(he_cw, halfedge_graph); - halfedge_descriptor he_ccw = prev( opposite(he, halfedge_graph), halfedge_graph ); - vertex_descriptor v3 = source(he_ccw, halfedge_graph); - - return ( half_tan_value_2(v1, v0, v2)/norm + half_tan_value_2(v1, v0, v3)/norm); - } - } - -private: - // Returns the tangent value of half angle v0_v1_v2/2 - double half_tan_value(vertex_descriptor v0, vertex_descriptor v1, vertex_descriptor v2) - { - Surface_modeling::Vector vec0(v2->point(), v1->point()); - Surface_modeling::Vector vec1(v0->point(), v2->point()); - Surface_modeling::Vector vec2(v1->point(), v0->point()); - double e0_square = vec0.squared_length(); - double e1_square = vec1.squared_length(); - double e2_square = vec2.squared_length(); - double e0 = std::sqrt(e0_square); - double e2 = std::sqrt(e2_square); - double cos_angle = ( e0_square + e2_square - e1_square ) / 2.0 / e0 / e2; - cos_angle = (std::max)(-1.0, (std::min)(1.0, cos_angle)); // clamp into [-1, 1] - double angle = acos(cos_angle); - - return ( tan(angle/2.0) ); - } - - // My deviation built on Meyer_02 - double half_tan_value_2(vertex_descriptor v0, vertex_descriptor v1, vertex_descriptor v2) - { - Surface_modeling::Vector a(v1->point(), v0->point()); - Surface_modeling::Vector b(v1->point(), v2->point()); - double dot_ab = a[0]*b[0] + a[1]*b[1] + a[2]*b[2]; - double dot_aa = a.squared_length(); - double dot_bb = b.squared_length(); - double dot_aa_bb = dot_aa * dot_bb; - - double cos_rep = dot_ab; - double sin_rep = std::sqrt(dot_aa_bb - dot_ab * dot_ab); - double normalizer = std::sqrt(dot_aa_bb); // |a|*|b| - - return (normalizer - cos_rep) / sin_rep; // formula from [Floater04] page 4 - // tan(Q/2) = (1 - cos(Q)) / sin(Q) - } -}; - -template< class HalfedgeGraph, - class PrimaryWeight = Cotangent_weight, - class SecondaryWeight = Mean_value_weight > -class Hybrid_weight : public PrimaryWeight, SecondaryWeight -{ -public: - typedef typename boost::graph_traits::halfedge_descriptor halfedge_descriptor; - typedef HalfedgeGraph Halfedge_graph; - - double operator()(halfedge_descriptor he, HalfedgeGraph& halfedge_graph) - { - double weight = PrimaryWeight::operator()(he, halfedge_graph); - //if(weight < 0) { std::cout << "Negative weight" << std::endl; } - return (weight >= 0) ? weight : SecondaryWeight::operator()(he, halfedge_graph); - } -}; - -// Trivial uniform weights (created for test purposes) -template -class Uniform_weight -{ -public: - typedef HalfedgeGraph Halfedge_graph; - typedef typename boost::graph_traits::halfedge_descriptor halfedge_descriptor; - - double operator()(halfedge_descriptor /*he*/, HalfedgeGraph& /*halfedge_graph*/) - { return 1.0; } -}; - - - -}//namespace internal -/// @endcond -}//namespace CGAL -#endif //CGAL_SURFACE_MODELING_WEIGHTS_H