diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/Weights.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/Weights.h deleted file mode 100644 index 22c308757e1..00000000000 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/Weights.h +++ /dev/null @@ -1,917 +0,0 @@ -// Copyright (c) 2015 GeometryFactory (France). -// All rights reserved. -// -// This file is part of CGAL (www.cgal.org). -// -// $URL$ -// $Id$ -// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial -// -// -// Author(s) : Yin Xu, Andreas Fabri and Ilker O. Yaz - -#ifndef CGAL_PMP_WEIGHTS_H -#define CGAL_PMP_WEIGHTS_H - -#include - -/// @cond CGAL_DOCUMENT_INTERNAL - -#include -#include -#include - -namespace CGAL { -namespace internal { - -// 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 -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 = to_double(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 = to_double(CGAL::approximate_sqrt(cross_ab*cross_ab)); - - if(divider == 0 /*|| divider != divider*/) - { - CGAL::collinear(get(ppmap, v0), get(ppmap, v1), get(ppmap, v2)) ? - CGAL_warning_msg(false, "Infinite Cotangent value with degenerate triangle!") : - CGAL_warning_msg(false, "Infinite Cotangent value due to floating point arithmetic!"); - - - return dot_ab > 0 ? (std::numeric_limits::max)() : - -(std::numeric_limits::max)(); - } - - return dot_ab / divider; - } -}; - -// Same as above but with a different API -template::type> -class Cotangent_value_Meyer -{ -protected: - typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; - typedef VertexPointMap Point_property_map; - typedef typename boost::property_traits::value_type Point; - typedef typename Kernel_traits::Kernel::Vector_3 Vector; - - PolygonMesh& pmesh_; - Point_property_map ppmap_; - -public: - - Cotangent_value_Meyer(PolygonMesh& pmesh_, VertexPointMap vpmap_) - : pmesh_(pmesh_) - , ppmap_(vpmap_) - {} - - PolygonMesh& pmesh() - { - 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()); - } -}; - - -// imported from skeletonization -template::type> -class Cotangent_value_Meyer_secure -{ - typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; - typedef VertexPointMap Point_property_map; - typedef typename boost::property_traits::value_type Point; - typedef typename Kernel_traits::Kernel::Vector_3 Vector; - - PolygonMesh& pmesh_; - Point_property_map ppmap_; - -public: - - Cotangent_value_Meyer_secure(PolygonMesh& pmesh_, VertexPointMap vpmap_) - : pmesh_(pmesh_) - , ppmap_(vpmap_) - {} - - PolygonMesh& pmesh() - { - return pmesh_; - } - - Point_property_map& ppmap() - { - return ppmap_; - } - - 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 = CGAL::to_double(a * b); - double dot_aa = CGAL::to_double(a.squared_length()); - double dot_bb = CGAL::to_double(b.squared_length()); - double lb = -0.999, ub = 0.999; - double cosine = dot_ab / CGAL::sqrt(dot_aa) / CGAL::sqrt(dot_bb); - cosine = (cosine < lb) ? lb : cosine; - cosine = (cosine > ub) ? ub : cosine; - double sine = std::sqrt(1.0 - cosine * cosine); - return cosine / sine; - } -}; - -// Returns the cotangent value of half angle v0 v1 v2 by clamping between [1, 89] degrees -// as suggested by -[Friedel] Unconstrained Spherical Parameterization- -template::type - , class CotangentValue = Cotangent_value_Meyer > -class Cotangent_value_clamped : CotangentValue -{ - Cotangent_value_clamped() - {} -public: - - Cotangent_value_clamped(PolygonMesh& pmesh_, VertexPointMap vpmap_) - : CotangentValue(pmesh_, vpmap_) - {} - - PolygonMesh& pmesh() - { - 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) - { - const double cot_1 = 57.289962; - const double cot_89 = 0.017455; - double value = CotangentValue::operator()(v0, v1, v2); - return (std::max)(cot_89, (std::min)(value, cot_1)); - } -}; - -template::type - , class CotangentValue = Cotangent_value_Meyer > -class Cotangent_value_clamped_2 : CotangentValue -{ - Cotangent_value_clamped_2() - {} - -public: - - Cotangent_value_clamped_2(PolygonMesh& pmesh_, VertexPointMap vpmap_) - : CotangentValue(pmesh_, vpmap_) - {} - - PolygonMesh& pmesh() - { - 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) - { - const double cot_5 = 5.671282; - const double cot_175 = -cot_5; - double value = CotangentValue::operator()(v0, v1, v2); - return (std::max)(cot_175, (std::min)(value, cot_5)); - } -}; - -template::type - , class CotangentValue = Cotangent_value_Meyer > -class Cotangent_value_minimum_zero : CotangentValue -{ - Cotangent_value_minimum_zero() - {} -public: - Cotangent_value_minimum_zero(PolygonMesh& pmesh_, VertexPointMap vpmap_) - : CotangentValue(pmesh_, vpmap_) - {} - - PolygonMesh& pmesh() - { - 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) - { - double value = CotangentValue::operator()(v0, v1, v2); - return (std::max)(0.0, value); - } -}; - -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 - = Cotangent_value_Meyer > -class Voronoi_area : CotangentValue -{ -public: - Voronoi_area(PolygonMesh& pmesh_, VertexPointMap vpmap_) - : CotangentValue(pmesh_, vpmap_) - {} - - PolygonMesh& pmesh() - { - 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; - - typedef VertexPointMap Point_property_map; - typedef typename boost::property_traits::value_type Point; - typedef typename Kernel_traits::Kernel::Vector_3 Vector; - - double operator()(vertex_descriptor v0) { - - //return 1.0; - double voronoi_area = 0.0; - for(halfedge_descriptor he : - halfedges_around_target( halfedge(v0,pmesh()), pmesh()) ) - { - if( is_border(he,pmesh()) ) { continue; } - - CGAL_expensive_assertion(CGAL::is_valid_polygon_mesh(pmesh())); - CGAL_expensive_assertion(CGAL::is_triangle_mesh(pmesh())); - CGAL_assertion( v0 == target(he, pmesh()) ); - vertex_descriptor v1 = source(he, pmesh()); - vertex_descriptor v_op = target(next(he, pmesh()), pmesh()); - - const Point& v0_p = get(ppmap(), v0); - const Point& v1_p = get(ppmap(), v1); - const Point& v_op_p = get(ppmap(), v_op); - - // (?) check if there is a better way to predicate triangle is obtuse or not - CGAL::Angle angle0 = CGAL::angle(v1_p, v0_p, v_op_p); - CGAL::Angle angle1 = CGAL::angle(v_op_p, v1_p, v0_p); - CGAL::Angle angle_op = CGAL::angle(v0_p, v_op_p, v1_p); - - bool obtuse = (angle0 == CGAL::OBTUSE) || (angle1 == CGAL::OBTUSE) || (angle_op == CGAL::OBTUSE); - - if(!obtuse) { - double cot_v1 = CotangentValue::operator()(v_op, v1, v0); - double cot_v_op = CotangentValue::operator()(v0, v_op, v1); - - double term1 = cot_v1 * to_double((v_op_p - v0_p).squared_length()); - double term2 = cot_v_op * to_double((v1_p - v0_p).squared_length()); - voronoi_area += (1.0 / 8.0) * (term1 + term2); - } - else { - double area_t = to_double(CGAL::approximate_sqrt(squared_area(v0_p, v1_p, v_op_p))); - if(angle0 == CGAL::OBTUSE) { - voronoi_area += area_t / 2.0; - } - else { - voronoi_area += area_t / 4.0; - } - } - } - CGAL_warning_msg(voronoi_area != 0, "Zero voronoi area!"); - return voronoi_area; - } -}; -// Returns the cotangent value of half angle v0 v1 v2 by dividing the triangle area -// as suggested by -[Mullen08] Spectral Conformal Parameterization- -template::type - , class CotangentValue = Cotangent_value_Meyer > -class Cotangent_value_area_weighted : CotangentValue -{ - Cotangent_value_area_weighted() - {} - -public: - - Cotangent_value_area_weighted(PolygonMesh& pmesh_, VertexPointMap vpmap_) - : CotangentValue(pmesh_, vpmap_) - {} - - PolygonMesh& pmesh() - { - 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) - { - return CotangentValue::operator()(v0, v1, v2) / - CGAL::sqrt(squared_area(get(this->ppmap(), v0), - get(this->ppmap(), v1), - get(this->ppmap(), v2))); - } -}; -///////////////////////////////////////////////////////////////////////////////////////// - -///////////////////////////// Edge 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 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 - = Cotangent_value_minimum_zero > -class Cotangent_weight : CotangentValue -{ - Cotangent_weight() - {} - -public: - Cotangent_weight(PolygonMesh& pmesh_, VertexPointMap vpmap_) - : CotangentValue(pmesh_, vpmap_) - {} - - Cotangent_weight(PolygonMesh& pmesh_) - : CotangentValue(pmesh_, get(CGAL::vertex_point, pmesh_)) - {} - - PolygonMesh& pmesh() - { - 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; - - typedef typename boost::property_map::type Point_property_map; - typedef typename boost::property_traits::value_type Point; - typedef typename Kernel_traits::Kernel::Vector_3 Vector; - - // Returns the cotangent weight of specified halfedge_descriptor - // Edge orientation is trivial - 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()); - } - return ( CotangentValue::operator()(v0, v2, v1)/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)/2.0 + CotangentValue::operator()(v0, v3, v1)/2.0 ); - } - } -}; - -// 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 > -class Single_cotangent_weight : CotangentValue -{ - Single_cotangent_weight() - {} -public: - Single_cotangent_weight(PolygonMesh& pmesh_, VertexPointMap vpmap_) - : CotangentValue(pmesh_, vpmap_) - {} - - PolygonMesh& pmesh() - { - 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; - - typedef typename boost::property_map::type Point_property_map; - typedef typename boost::property_traits::value_type Point; - typedef typename Kernel_traits::Kernel::Vector_3 Vector; - - // Returns the cotangent of the opposite angle of the edge - // 0 for border edges (which does not have an opposite angle) - double operator()(halfedge_descriptor he) - { - 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); - } -}; - - -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> -class Mean_value_weight -{ - //Mean_value_weight() - //{} - - PolygonMesh& pmesh_; - VertexPointMap vpmap_; - -public: - Mean_value_weight(PolygonMesh& pmesh_, VertexPointMap vpmap) - : pmesh_(pmesh_), vpmap_(vpmap) - {} - - PolygonMesh& pmesh() - { - return pmesh_; - } - - typedef typename boost::graph_traits::halfedge_descriptor halfedge_descriptor; - typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; - - typedef VertexPointMap Point_property_map; - typedef typename boost::property_traits::value_type Point; - typedef typename Kernel_traits::Kernel::Vector_3 Vector; - - // Returns the mean-value coordinate of specified halfedge_descriptor - // Returns different value for different edge orientation (which is a normal behaivour according to formula) - double operator()(halfedge_descriptor he) - { - vertex_descriptor v0 = target(he, pmesh()); - vertex_descriptor v1 = source(he, pmesh()); - Vector vec = get(vpmap_, v0) - get(vpmap_, v1); - double norm = CGAL::sqrt( vec.squared_length() ); - - // 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 ( half_tan_value_2(v1, v0, v2)/norm); - } - 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 ( 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) - { - Vector vec0 = get(vpmap_, v1) - get(vpmap_, v2); - Vector vec1 = get(vpmap_, v2) - get(vpmap_, v0); - Vector vec2 = get(vpmap_, v0) - get(vpmap_, v1); - double e0_square = vec0.squared_length(); - double e1_square = vec1.squared_length(); - double e2_square = vec2.squared_length(); - double e0 = CGAL::sqrt(e0_square); - double e2 = CGAL::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) - { - Vector a = get(vpmap_, v0) - get(vpmap_, v1); - Vector b = get(vpmap_, v2) - get(vpmap_, v1); - 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 = CGAL::sqrt(dot_aa_bb - dot_ab * dot_ab); - double normalizer = CGAL::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 PolygonMesh, - class PrimaryWeight = Cotangent_weight, - class SecondaryWeight = Mean_value_weight > -class Hybrid_weight : public PrimaryWeight, SecondaryWeight -{ - PrimaryWeight primary; - SecondaryWeight secondary; - - Hybrid_weight() - {} - -public: - Hybrid_weight(PolygonMesh& pmesh_) - : primary(pmesh_), secondary(pmesh_) - {} - - PolygonMesh& pmesh() - { - return primary.pmesh(); - } - - typedef typename boost::graph_traits::halfedge_descriptor halfedge_descriptor; - - double operator()(halfedge_descriptor he) - { - double weight = primary(he); - //if(weight < 0) { std::cout << "Negative weight" << std::endl; } - return (weight >= 0) ? weight : secondary(he); - } -}; - -// Trivial uniform weights (created for test purposes) -template -class Uniform_weight -{ -public: - typedef typename boost::graph_traits::halfedge_descriptor halfedge_descriptor; - - double operator()(halfedge_descriptor /*e*/) - { return 1.0; } -}; - -//////////////////////////////////////////////////////////////////////////// -// FAIRING // -template -class Scale_dependent_weight_fairing -{ - PolygonMesh& pmesh_; -public: - Scale_dependent_weight_fairing(PolygonMesh& pmesh_) - : pmesh_(pmesh_) - {} - - PolygonMesh& pmesh() - { - return pmesh_; - } - - typedef typename boost::graph_traits::halfedge_descriptor halfedge_descriptor; - typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; - - typedef typename boost::property_map::type Point_property_map; - typedef typename boost::property_traits::value_type Point; - typedef typename Kernel_traits::Kernel::Vector_3 Vector; - - double w_i(vertex_descriptor /*v_i*/) { return 1.0; } - - double w_ij(halfedge_descriptor he) - { - Vector v = target(he, pmesh())->point() - source(he, pmesh())->point(); - double divider = CGAL::sqrt(v.squared_length()); - if(divider == 0.0) { - CGAL_warning_msg(false, "Scale dependent weight - zero length edge."); - return (std::numeric_limits::max)(); - } - return 1.0 / divider; - } -}; - -template::type -> -class Cotangent_weight_with_voronoi_area_fairing -{ - typedef PolygonMesh PM; - typedef VertexPointMap VPMap; - Voronoi_area voronoi_functor; - Cotangent_weight > cotangent_functor; - -public: - Cotangent_weight_with_voronoi_area_fairing(PM& pmesh_) - : voronoi_functor(pmesh_, get(CGAL::vertex_point, pmesh_)) - , cotangent_functor(pmesh_, get(CGAL::vertex_point, pmesh_)) - {} - - Cotangent_weight_with_voronoi_area_fairing(PM& pmesh_, VPMap vpmap_) - : voronoi_functor(pmesh_, vpmap_) - , cotangent_functor(pmesh_, vpmap_) - {} - - PM& pmesh() - { - return voronoi_functor.pmesh(); - } - - typedef typename boost::graph_traits::halfedge_descriptor halfedge_descriptor; - typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; - - double w_i(vertex_descriptor v_i) - { - return 0.5 / voronoi_functor(v_i); - } - - double w_ij(halfedge_descriptor he) { - - return cotangent_functor(he) * 2.0; - } -}; - -// Cotangent_value_Meyer has been changed to the version: -// Cotangent_value_Meyer_secure to avoid imprecisions from -// the issue #4706 - https://github.com/CGAL/cgal/issues/4706. -template< -class PolygonMesh, class VertexPointMap = typename boost::property_map::type> -class Cotangent_weight_with_voronoi_area_fairing_secure { - - typedef PolygonMesh PM; - typedef VertexPointMap VPMap; - Voronoi_area voronoi_functor; - Cotangent_weight > cotangent_functor; - -public: - Cotangent_weight_with_voronoi_area_fairing_secure(PM& pmesh_) : - voronoi_functor(pmesh_, get(CGAL::vertex_point, pmesh_)), - cotangent_functor(pmesh_, get(CGAL::vertex_point, pmesh_)) - { } - - Cotangent_weight_with_voronoi_area_fairing_secure(PM& pmesh_, VPMap vpmap_) : - voronoi_functor(pmesh_, vpmap_), - cotangent_functor(pmesh_, vpmap_) - { } - - PM& pmesh() { - return voronoi_functor.pmesh(); - } - - typedef typename boost::graph_traits::halfedge_descriptor halfedge_descriptor; - typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; - - double w_i(vertex_descriptor v_i) { - return 0.5 / voronoi_functor(v_i); - } - - double w_ij(halfedge_descriptor he) { - return cotangent_functor(he) * 2.0; - } -}; - -template -class Uniform_weight_fairing -{ -public: - typedef typename boost::graph_traits::halfedge_descriptor halfedge_descriptor; - typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; - - Uniform_weight_fairing(PolygonMesh&) - {} - - double w_ij(halfedge_descriptor /*e*/) { return 1.0; } - - double w_i(vertex_descriptor /*v_i*/) { return 1.0; } -}; -//////////////////////////////////////////////////////////////////////////// - -}//namespace internal - - -}//namespace CGAL -/// @endcond -#endif //CGAL_PMP_WEIGHTS_H diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/fair.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/fair.h index 727e007995a..1aa166f8b6f 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/fair.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/fair.h @@ -20,6 +20,7 @@ #include #include #include +#include #if defined(CGAL_EIGEN3_ENABLED) #include // for sparse linear system solver @@ -162,8 +163,8 @@ namespace internal { // Cotangent_weight_with_voronoi_area_fairing has been changed to the version: // Cotangent_weight_with_voronoi_area_fairing_secure to avoid imprecisions from // the issue #4706 - https://github.com/CGAL/cgal/issues/4706. - typedef CGAL::internal::Cotangent_weight_with_voronoi_area_fairing_secure - Default_Weight_calculator; + typedef CGAL::Weights::internal::Cotangent_weight_wrapper_with_voronoi_secure< + TriangleMesh, VPMap> Default_Weight_calculator; VPMap vpmap_ = choose_parameter(get_parameter(np, internal_np::vertex_point), get_property_map(vertex_point, tmesh)); diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Smoothing/curvature_flow_impl.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Smoothing/curvature_flow_impl.h index 0c9e2d9fa0e..2d75758ecae 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Smoothing/curvature_flow_impl.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Smoothing/curvature_flow_impl.h @@ -17,10 +17,9 @@ #include #include -#include #include #include - +#include #include #include @@ -41,45 +40,6 @@ namespace CGAL { namespace Polygon_mesh_processing { namespace internal { -// Empirically, _Meyer seems to produce the best results from the various weights available in Weights.h -template > -struct Edge_cotangent_weight - : public CotangentValue -{ - typedef typename boost::graph_traits::halfedge_descriptor halfedge_descriptor; - typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; - - Edge_cotangent_weight(TriangleMesh& pmesh_, VertexPointMap vpmap_) : CotangentValue(pmesh_, vpmap_) {} - - TriangleMesh& pmesh() { return CotangentValue::pmesh(); } - - double operator()(halfedge_descriptor he) - { - if(is_border_edge(he, pmesh())) - { - halfedge_descriptor h1 = next(he, pmesh()); - vertex_descriptor vs = source(he, pmesh()); - vertex_descriptor vt = target(he, pmesh()); - vertex_descriptor v1 = target(h1, pmesh()); - - return CotangentValue::operator()(vs, v1, vt); - } - else - { - halfedge_descriptor h1 = next(he, pmesh()); - halfedge_descriptor h2 = prev(opposite(he, pmesh()), pmesh()); - vertex_descriptor vs = source(he, pmesh()); - vertex_descriptor vt = target(he, pmesh()); - vertex_descriptor v1 = target(h1, pmesh()); - vertex_descriptor v2 = source(h2, pmesh()); - - return CotangentValue::operator()(vs, v1, vt) + CotangentValue::operator()(vs, v2, vt); - } - } -}; - template constrained_flags_; const GeomTraits& traits_; - Edge_cotangent_weight weight_calculator_; + CGAL::Weights::internal::Edge_cotangent_weight_wrapper< + TriangleMesh, VertexPointMap> weight_calculator_; }; } // internal diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/fair_impl.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/fair_impl.h index a45ff4b55e9..71b325f5d88 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/fair_impl.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/fair_impl.h @@ -19,7 +19,6 @@ #include #include #include -#include #ifdef CGAL_PMP_FAIR_DEBUG #include #endif @@ -65,7 +64,7 @@ private: double weight = 0; Halfedge_around_vertex_circulator circ(halfedge(v,pmesh),pmesh), done(circ); do { - weight += weight_calculator.w_ij(*circ); + weight += CGAL::to_double(weight_calculator.w_ij(*circ)); } while(++circ != done); return weight; } @@ -95,11 +94,11 @@ private: } } else { - double w_i = weight_calculator.w_i(v); + double w_i = CGAL::to_double(weight_calculator.w_i(v)); Halfedge_around_vertex_circulator circ(halfedge(v,pmesh),pmesh), done(circ); do { - double w_i_w_ij = w_i * weight_calculator.w_ij(*circ) ; + double w_i_w_ij = w_i * CGAL::to_double(weight_calculator.w_ij(*circ)) ; vertex_descriptor nv = target(opposite(*circ,pmesh),pmesh); compute_row(nv, row_id, matrix, x, y, z, -w_i_w_ij*multiplier, vertex_id_map, depth-1); diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/repair_self_intersections.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/repair_self_intersections.h index f5fdc39b066..5bc9190b83d 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/repair_self_intersections.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/repair_self_intersections.h @@ -111,7 +111,7 @@ FaceOutputIterator replace_faces_with_patch(const std::vector point_to_vs; // first, add those for which the vertex will not change - for(const vertex_descriptor v : border_vertices) + for(const vertex_descriptor& v : border_vertices) point_to_vs[get(vpm, v)] = v; // now build a correspondence map and the faces with vertices @@ -505,7 +505,7 @@ bool remove_self_intersections_with_smoothing(std::set > patch; - for(const face_descriptor f : faces(local_mesh)) + for(const face_descriptor& f : faces(local_mesh)) { halfedge_descriptor h = halfedge(f, local_mesh); patch.emplace_back(std::initializer_list{get(local_vpm, target(h, local_mesh)), @@ -702,7 +702,7 @@ bool construct_tentative_hole_patch(std::vector internal_hedges; std::vector cc_border_hedges; - for(const face_descriptor fd : sub_cc_faces) + for(const face_descriptor& fd : sub_cc_faces) { halfedge_descriptor h = halfedge(fd, tmesh); for(int i=0; i<3;++i) @@ -933,7 +933,7 @@ bool construct_tentative_sub_hole_patch(std::vector::halfedge_ // Could renew the range directly within the patch replacement function // to avoid erasing and re-adding the same face - for(const face_descriptor f : cc_faces) + for(const face_descriptor& f : cc_faces) working_face_range.erase(f); // Plug the new triangles in the mesh, reusing previous edges and faces @@ -1256,7 +1256,7 @@ bool fill_hole_with_constraints(std::vector +class Uniform_weight_fairing { + +public: + using FT = typename GeomTraits::FT; + using halfedge_descriptor = typename boost::graph_traits::halfedge_descriptor; + using vertex_descriptor = typename boost::graph_traits::vertex_descriptor; + + FT w_ij(halfedge_descriptor) { return FT(1); } + FT w_i(vertex_descriptor) { return FT(1); } +}; template void read_poly(const char* file_name, Polyhedron& poly) { @@ -315,7 +328,7 @@ void test_ouput_iterators_triangulate_and_refine_hole(const char* file_name) { std::cout << " Done!" << std::endl; } -template +template void test_triangulate_refine_and_fair_hole_compile() { typedef typename boost::graph_traits::halfedge_descriptor Halfedge_handle; typedef typename boost::graph_traits::face_descriptor Facet_handle; @@ -336,16 +349,16 @@ void test_triangulate_refine_and_fair_hole_compile() { read_poly_with_borders("elephant_quad_hole.off", poly, border_reps); CGAL::Polygon_mesh_processing::triangulate_refine_and_fair_hole (poly, border_reps[0], back_inserter(patch_facets), back_inserter(patch_vertices), - CGAL::Polygon_mesh_processing::parameters::weight_calculator( - CGAL::internal::Uniform_weight_fairing(poly)). + CGAL::Polygon_mesh_processing::parameters:: + weight_calculator(Uniform_weight_fairing()). sparse_linear_solver(Default_solver())); // default solver read_poly_with_borders("elephant_quad_hole.off", poly, border_reps); CGAL::Polygon_mesh_processing::triangulate_refine_and_fair_hole (poly, border_reps[0], back_inserter(patch_facets), back_inserter(patch_vertices), - CGAL::Polygon_mesh_processing::parameters::weight_calculator( - CGAL::internal::Uniform_weight_fairing(poly))); + CGAL::Polygon_mesh_processing::parameters:: + weight_calculator(Uniform_weight_fairing())); // default solver and weight read_poly_with_borders("elephant_quad_hole.off", poly, border_reps); @@ -404,7 +417,7 @@ typedef CGAL::Surface_mesh Polyhedron; test_triangulate_hole_should_be_no_output("data/non_manifold_vertex.off"); test_triangulate_hole_should_be_no_output("data/two_tris_collinear.off"); - test_triangulate_refine_and_fair_hole_compile(); + test_triangulate_refine_and_fair_hole_compile(); } diff --git a/Polygon_mesh_processing/test/Polygon_mesh_processing/triangulate_hole_Polyhedron_3_test.cpp b/Polygon_mesh_processing/test/Polygon_mesh_processing/triangulate_hole_Polyhedron_3_test.cpp index fd8e8448a7b..6d4ee44ef03 100644 --- a/Polygon_mesh_processing/test/Polygon_mesh_processing/triangulate_hole_Polyhedron_3_test.cpp +++ b/Polygon_mesh_processing/test/Polygon_mesh_processing/triangulate_hole_Polyhedron_3_test.cpp @@ -30,6 +30,20 @@ typedef boost::graph_traits::halfedge_iterator Halfedge_iterator; typedef CGAL::Halfedge_around_face_circulator Halfedge_around_facet_circulator; typedef boost::property_map::type Point_property_map; +template< +class GeomTraits, +class PolygonMesh> +class Uniform_weight_fairing { + +public: + using FT = typename GeomTraits::FT; + using halfedge_descriptor = typename boost::graph_traits::halfedge_descriptor; + using vertex_descriptor = typename boost::graph_traits::vertex_descriptor; + + FT w_ij(halfedge_descriptor) { return FT(1); } + FT w_i(vertex_descriptor) { return FT(1); } +}; + void read_poly(const char* file_name, Polyhedron& poly) { poly.clear(); @@ -326,16 +340,17 @@ void test_triangulate_refine_and_fair_hole_compile() { read_poly_with_borders("elephant_quad_hole.off", poly, border_reps); CGAL::Polygon_mesh_processing::triangulate_refine_and_fair_hole (poly, border_reps[0], back_inserter(patch_facets), back_inserter(patch_vertices), - CGAL::Polygon_mesh_processing::parameters::weight_calculator( - CGAL::internal::Uniform_weight_fairing(poly)). - sparse_linear_solver(Default_solver()).use_2d_constrained_delaunay_triangulation(false)); + CGAL::Polygon_mesh_processing::parameters:: + weight_calculator(Uniform_weight_fairing()). + sparse_linear_solver(Default_solver()). + use_2d_constrained_delaunay_triangulation(false)); // default solver read_poly_with_borders("elephant_quad_hole.off", poly, border_reps); CGAL::Polygon_mesh_processing::triangulate_refine_and_fair_hole (poly, border_reps[0], back_inserter(patch_facets), back_inserter(patch_vertices), - CGAL::Polygon_mesh_processing::parameters::weight_calculator( - CGAL::internal::Uniform_weight_fairing(poly))); + CGAL::Polygon_mesh_processing::parameters:: + weight_calculator(Uniform_weight_fairing())); // default solver and weight read_poly_with_borders("elephant_quad_hole.off", poly, border_reps); diff --git a/Weights/include/CGAL/Weights/internal/tools.h b/Weights/include/CGAL/Weights/internal/tools.h index a7c87769079..430fbf5f773 100644 --- a/Weights/include/CGAL/Weights/internal/tools.h +++ b/Weights/include/CGAL/Weights/internal/tools.h @@ -39,8 +39,6 @@ decltype(auto) cotangent_3_secure( const auto dot_product_3 = traits.compute_scalar_product_3_object(); - const auto cross_product_3 = - traits.construct_cross_product_vector_3_object(); const auto construct_vector_3 = traits.construct_vector_3_object(); @@ -78,11 +76,11 @@ public: CGAL_assertion(m_d_r != FT(0)); // two points are identical! m_d_p = distance(q, p); CGAL_assertion(m_d_p != FT(0)); // two points are identical! - const auto area = area(p, q, r); - CGAL_assertion(area != FT(0)); // three points are identical! + const auto A = area(p, q, r); + CGAL_assertion(A != FT(0)); // three points are identical! const auto scalar = scalar_product(p, q, r); - m_w_base = -tangent_half_angle(m_d_r, m_d_p, area, scalar); + m_w_base = -tangent_half_angle(m_d_r, m_d_p, A, scalar); } FT get_w_r() const { @@ -256,7 +254,7 @@ private: FT voronoi_area = FT(0); CGAL_assertion(CGAL::is_triangle_mesh(m_pmesh)); - for (const auto he : halfedges_around_target(halfedge(v0, m_pmesh), m_pmesh)) { + for (const auto& he : halfedges_around_target(halfedge(v0, m_pmesh), m_pmesh)) { CGAL_assertion(v0 == target(he, m_pmesh)); if (is_border(he, m_pmesh)) { continue;