diff --git a/Heat_method_3/doc/Heat_method_3/Concepts/HeatMethodTraits_3.h b/Heat_method_3/doc/Heat_method_3/Concepts/HeatMethodTraits_3.h index f33440e95cf..dcfe6d573e0 100644 --- a/Heat_method_3/doc/Heat_method_3/Concepts/HeatMethodTraits_3.h +++ b/Heat_method_3/doc/Heat_method_3/Concepts/HeatMethodTraits_3.h @@ -49,6 +49,8 @@ public: /// Functor with operator: `FT operator()(const Point_3& p, const Point_3& q) const` which computes the squared distance between `p` and `q`. typedef unspecified_type Compute_squared_distance_3; + /// Functor with operator: `FT operator()(const Vector_3& v) const` which computes the squared length of `v`. + typedef unspecified_type Compute_squared_length_3; /// @} diff --git a/Heat_method_3/include/CGAL/Heat_method_3/Surface_mesh_geodesic_distances_3.h b/Heat_method_3/include/CGAL/Heat_method_3/Surface_mesh_geodesic_distances_3.h index 96f7846baa6..d8592b1346f 100644 --- a/Heat_method_3/include/CGAL/Heat_method_3/Surface_mesh_geodesic_distances_3.h +++ b/Heat_method_3/include/CGAL/Heat_method_3/Surface_mesh_geodesic_distances_3.h @@ -29,6 +29,7 @@ #include #endif +#include #include #include @@ -361,7 +362,6 @@ private: void compute_divergence() { - typename Traits::Construct_cross_product_vector_3 cross_product = Traits().construct_cross_product_vector_3_object(); typename Traits::Compute_scalar_product_3 scalar_product = Traits().compute_scalar_product_3_object(); typename Traits::Construct_vector_3 construct_vector = Traits().construct_vector_3_object(); Matrix indexD(dimension,1); @@ -374,36 +374,30 @@ private: Index i = get(vertex_id_map, current); Index j = get(vertex_id_map, neighbor_one); Index k = get(vertex_id_map, neighbor_two); - VertexPointMap_reference p_i = get(vpm,current); + VertexPointMap_reference p_i = get(vpm, current); VertexPointMap_reference p_j = get(vpm, neighbor_one); VertexPointMap_reference p_k = get(vpm, neighbor_two); Index face_i = get(face_id_map, f); - Vector_3 v_ij = construct_vector(p_i,p_j); - Vector_3 v_ik = construct_vector(p_i,p_k); - Vector_3 cross = cross_product(v_ij, v_ik); - double norm_cross = CGAL::sqrt(to_double(scalar_product(cross,cross))); - double dot = to_double(scalar_product(v_ij, v_ik)); - double cotan_i = dot/norm_cross; + const Vector_3 v_ij = construct_vector(p_i, p_j); + const Vector_3 v_ik = construct_vector(p_i, p_k); + const Vector_3 v_ji = construct_vector(p_j, p_i); + const Vector_3 v_jk = construct_vector(p_j, p_k); + const Vector_3 v_ki = construct_vector(p_k, p_i); + const Vector_3 v_kj = construct_vector(p_k, p_j); - Vector_3 v_ji = construct_vector(p_j, p_i); - Vector_3 v_jk = construct_vector(p_j, p_k); + const Traits traits; + const FT cotan_i = CGAL::Weights::cotangent(p_k, p_i, p_j, traits); + const FT cotan_j = CGAL::Weights::cotangent(p_k, p_j, p_i, traits); + const FT cotan_k = CGAL::Weights::cotangent(p_j, p_k, p_i, traits); - cross = cross_product(v_ji, v_jk); - dot = to_double(scalar_product(v_ji, v_jk)); - double cotan_j = dot/norm_cross; - - Vector_3 v_ki = construct_vector(p_k,p_i); - Vector_3 v_kj = construct_vector(p_k,p_j); - - cross = cross_product(v_ki, v_kj); - dot = to_double(scalar_product(v_ki,v_kj)); - double cotan_k = dot/norm_cross; - - const Vector_3& a = m_X[face_i]; - double i_entry = (to_double(scalar_product(a,v_ij)) * cotan_k) + (to_double(scalar_product(a,v_ik)) * cotan_j); - double j_entry = (to_double(scalar_product(a,v_jk)) * cotan_i) + (to_double(scalar_product(a,v_ji)) * cotan_k); - double k_entry = (to_double(scalar_product(a,v_ki)) * cotan_j) + (to_double(scalar_product(a,v_kj)) * cotan_i); + const Vector_3& a = m_X[face_i]; + const double i_entry = (CGAL::to_double(scalar_product(a, v_ij) * cotan_k)) + + (CGAL::to_double(scalar_product(a, v_ik) * cotan_j)); + const double j_entry = (CGAL::to_double(scalar_product(a, v_jk) * cotan_i)) + + (CGAL::to_double(scalar_product(a, v_ji) * cotan_k)); + const double k_entry = (CGAL::to_double(scalar_product(a, v_ki) * cotan_j)) + + (CGAL::to_double(scalar_product(a, v_kj) * cotan_i)); indexD.add_coef(i, 0, (1./2)*i_entry); indexD.add_coef(j, 0, (1./2)*j_entry); @@ -549,47 +543,40 @@ private: Index k = get(vertex_id_map, neighbor_two); Point_3 pi, pj, pk; - VertexPointMap_reference p_i = get(vpm,current); + const Traits traits; + VertexPointMap_reference p_i = get(vpm, current); VertexPointMap_reference p_j = get(vpm, neighbor_one); VertexPointMap_reference p_k = get(vpm, neighbor_two); pi = p_i; pj = p_j; pk = p_k; - Vector_3 v_ij = construct_vector(p_i,p_j); - Vector_3 v_ik = construct_vector(p_i,p_k); + const double cotan_i = CGAL::to_double( + CGAL::Weights::cotangent(pk, pi, pj, traits)); + m_cotan_matrix.add_coef(j, k, -(1./2) * cotan_i); + m_cotan_matrix.add_coef(k, j, -(1./2) * cotan_i); + m_cotan_matrix.add_coef(j, j, (1./2) * cotan_i); + m_cotan_matrix.add_coef(k, k, (1./2) * cotan_i); - Vector_3 cross = cross_product(v_ij, v_ik); - double dot = to_double(scalar_product(v_ij,v_ik)); + const double cotan_j = CGAL::to_double( + CGAL::Weights::cotangent(pk, pj, pi, traits)); + m_cotan_matrix.add_coef(i, k, -(1./2) * cotan_j); + m_cotan_matrix.add_coef(k, i, -(1./2) * cotan_j); + m_cotan_matrix.add_coef(i, i, (1./2) * cotan_j); + m_cotan_matrix.add_coef(k, k, (1./2) * cotan_j); - double norm_cross = (CGAL::sqrt(to_double(scalar_product(cross,cross)))); + const double cotan_k = CGAL::to_double( + CGAL::Weights::cotangent(pj, pk, pi, traits)); + m_cotan_matrix.add_coef(i, j, -(1./2) * cotan_k); + m_cotan_matrix.add_coef(j, i, -(1./2) * cotan_k); + m_cotan_matrix.add_coef(i, i, (1./2) * cotan_k); + m_cotan_matrix.add_coef(j, j, (1./2) * cotan_k); - double cotan_i = dot/norm_cross; - m_cotan_matrix.add_coef(j,k ,-(1./2)*cotan_i); - m_cotan_matrix.add_coef(k,j,-(1./2)* cotan_i); - m_cotan_matrix.add_coef(j,j,(1./2)*cotan_i); - m_cotan_matrix.add_coef(k,k,(1./2)* cotan_i); - - Vector_3 v_ji = construct_vector(p_j,p_i); - Vector_3 v_jk = construct_vector(p_j,p_k); - - cross = cross_product(v_ji, v_jk); - dot = to_double(scalar_product(v_ji, v_jk)); - double cotan_j = dot/norm_cross; - m_cotan_matrix.add_coef(i,k ,-(1./2)*cotan_j); - m_cotan_matrix.add_coef(k,i,-(1./2)* cotan_j); - m_cotan_matrix.add_coef(i,i,(1./2)* cotan_j); - m_cotan_matrix.add_coef(k,k,(1./2)* cotan_j); - - Vector_3 v_ki = construct_vector(p_k,p_i); - Vector_3 v_kj = construct_vector(p_k,p_j); - cross = cross_product(v_ki, v_kj); - dot = to_double(scalar_product(v_ki,v_kj)); - double cotan_k = dot/norm_cross; - m_cotan_matrix.add_coef(i,j,-(1./2)*cotan_k); - m_cotan_matrix.add_coef(j,i,-(1./2)* cotan_k); - m_cotan_matrix.add_coef(i,i,(1./2)* cotan_k); - m_cotan_matrix.add_coef(j,j,(1./2)* cotan_k); + const Vector_3 v_ij = construct_vector(p_i, p_j); + const Vector_3 v_ik = construct_vector(p_i, p_k); + const Vector_3 cross = cross_product(v_ij, v_ik); + const double norm_cross = CGAL::sqrt( + CGAL::to_double(scalar_product(cross, cross))); //double area_face = CGAL::Polygon_mesh_processing::face_area(f,tm); //cross is 2*area diff --git a/Heat_method_3/package_info/Heat_method_3/dependencies b/Heat_method_3/package_info/Heat_method_3/dependencies index 90510db1ab2..50488d5544f 100644 --- a/Heat_method_3/package_info/Heat_method_3/dependencies +++ b/Heat_method_3/package_info/Heat_method_3/dependencies @@ -18,3 +18,4 @@ Random_numbers STL_Extension Solver_interface Stream_support +Weights diff --git a/Heat_method_3/test/Heat_method_3/Heat_method_traits_3.h b/Heat_method_3/test/Heat_method_3/Heat_method_traits_3.h index b64613c641b..8b8e28c76f4 100644 --- a/Heat_method_3/test/Heat_method_3/Heat_method_traits_3.h +++ b/Heat_method_3/test/Heat_method_3/Heat_method_traits_3.h @@ -59,6 +59,11 @@ struct Heat_method_traits_3 { return 0;} }; + struct Compute_squared_length_3 { + + double operator()(const Vector_3&) const + { return 0;} + }; Construct_vector_3 construct_vector_3_object() const { @@ -79,6 +84,11 @@ struct Heat_method_traits_3 return Compute_squared_distance_3(); } + Compute_squared_length_3 compute_squared_length_3_object() const + { + return Compute_squared_length_3(); + } + Compute_scalar_product_3 compute_scalar_product_3_object() const { return Compute_scalar_product_3(); diff --git a/Installation/CHANGES.md b/Installation/CHANGES.md index 27dc295dcf1..9ca67be1cc9 100644 --- a/Installation/CHANGES.md +++ b/Installation/CHANGES.md @@ -37,7 +37,7 @@ Release date: December 2021 ### [2D and 3D Linear Geometry Kernel](https://doc.cgal.org/5.4/Manual/packages.html#PkgKernel23) - Added `construct_centroid_2_object()` and `compute_determinant_2_object()` in `Projection_traits_xy_3`, `Projection_traits_xz_3`, - and`Projection_traits_yz_3` classes. + and `Projection_traits_yz_3` classes. - Added documentation for the class `Projection_traits_3`, which enables the use of 2D algorithms on the projections of 3D data onto an arbitrary plane. @@ -62,6 +62,10 @@ Release date: December 2021 - Added an option to [`CGAL::Polygon_mesh_processing::self_intersections()`](https://doc.cgal.org/5.4/Polygon_mesh_processing/group__PMP__intersection__grp.html#gaf19c80ec12cbff7ebe9e69453f1d40b8) to report only a limited number of intersections (`maximum_number()`) +### [The Heat Method](https://doc.cgal.org/5.4/Manual/packages.html#PkgHeatMethod) + +- **Breaking change**: Added the functor `Compute_squared_length_3` that has an operator `operator(const Vector_3& v)`, which computes the squared length of `v`, to the `HeatMethodTraits_3` concept. + ### [Shape Regularization](https://doc.cgal.org/5.4/Manual/packages.html#PkgShapeRegularization) (new package) - This package enables to regularize a set of segments and open or closed contours in 2D 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..b6e9d00b9d0 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,15 +163,14 @@ 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::Secure_cotangent_weight_with_voronoi_area Default_weight_calculator; VPMap vpmap_ = choose_parameter(get_parameter(np, internal_np::vertex_point), get_property_map(vertex_point, tmesh)); return internal::fair(tmesh, vertices, choose_parameter(get_parameter(np, internal_np::sparse_linear_solver)), - choose_parameter(get_parameter(np, internal_np::weight_calculator), Default_Weight_calculator(tmesh, vpmap_)), + choose_parameter(get_parameter(np, internal_np::weight_calculator), Default_weight_calculator(tmesh, vpmap_)), choose_parameter(get_parameter(np, internal_np::fairing_continuity), 1), vpmap_); } 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..fa600557af3 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_; + const CGAL::Weights::Edge_cotangent_weight 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..62e1ce5e5e8 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,7 @@ #include #include #include -#include +#include #ifdef CGAL_PMP_FAIR_DEBUG #include #endif @@ -65,7 +65,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 +95,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/package_info/Polygon_mesh_processing/dependencies b/Polygon_mesh_processing/package_info/Polygon_mesh_processing/dependencies index e0691224967..6d471658093 100644 --- a/Polygon_mesh_processing/package_info/Polygon_mesh_processing/dependencies +++ b/Polygon_mesh_processing/package_info/Polygon_mesh_processing/dependencies @@ -35,3 +35,4 @@ TDS_3 Triangulation_2 Triangulation_3 Union_find +Weights diff --git a/Polygon_mesh_processing/test/Polygon_mesh_processing/triangulate_hole_Polyhedron_3_no_delaunay_test.cpp b/Polygon_mesh_processing/test/Polygon_mesh_processing/triangulate_hole_Polyhedron_3_no_delaunay_test.cpp index 3dc4125306c..5e071d05055 100644 --- a/Polygon_mesh_processing/test/Polygon_mesh_processing/triangulate_hole_Polyhedron_3_no_delaunay_test.cpp +++ b/Polygon_mesh_processing/test/Polygon_mesh_processing/triangulate_hole_Polyhedron_3_no_delaunay_test.cpp @@ -16,6 +16,8 @@ #include +#include + #include #include #include @@ -336,16 +338,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(CGAL::Weights::Uniform_weight()). 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(CGAL::Weights::Uniform_weight())); // default solver and weight read_poly_with_borders("elephant_quad_hole.off", poly, border_reps); 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..a83038f56d3 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 @@ -11,6 +11,7 @@ #include #include #include +#include #include #include @@ -326,16 +327,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(CGAL::Weights::Uniform_weight()). + 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(CGAL::Weights::Uniform_weight())); // default solver and weight read_poly_with_borders("elephant_quad_hole.off", poly, border_reps); diff --git a/Polyhedron/demo/Polyhedron/Plugins/PMP/Fairing_plugin.cpp b/Polyhedron/demo/Polyhedron/Plugins/PMP/Fairing_plugin.cpp index 93c130d4854..d20d302e43f 100644 --- a/Polyhedron/demo/Polyhedron/Plugins/PMP/Fairing_plugin.cpp +++ b/Polyhedron/demo/Polyhedron/Plugins/PMP/Fairing_plugin.cpp @@ -13,6 +13,7 @@ typedef Scene_surface_mesh_item Scene_facegraph_item; #include #include #include +#include #include #include @@ -98,7 +99,8 @@ public Q_SLOTS: if(weight_index == 1) CGAL::Polygon_mesh_processing::fair(*selection_item->polyhedron(), selection_item->selected_vertices, - CGAL::Polygon_mesh_processing::parameters::weight_calculator(CGAL::internal::Uniform_weight_fairing(*selection_item->polyhedron())). + CGAL::Polygon_mesh_processing::parameters:: + weight_calculator(CGAL::Weights::Uniform_weight()). fairing_continuity(continuity)); if(weight_index == 0) CGAL::Polygon_mesh_processing::fair(*selection_item->polyhedron(), diff --git a/Polyhedron/demo/Polyhedron/Plugins/PMP/Hole_filling_plugin.cpp b/Polyhedron/demo/Polyhedron/Plugins/PMP/Hole_filling_plugin.cpp index d746af5fbad..c4a6207c968 100644 --- a/Polyhedron/demo/Polyhedron/Plugins/PMP/Hole_filling_plugin.cpp +++ b/Polyhedron/demo/Polyhedron/Plugins/PMP/Hole_filling_plugin.cpp @@ -19,6 +19,8 @@ #include #include #include +#include +#include #include #include @@ -716,16 +718,18 @@ bool Polyhedron_demo_hole_filling_plugin::fill if(weight_index == 0) { success = std::get<0>(CGAL::Polygon_mesh_processing::triangulate_refine_and_fair_hole(poly, it, std::back_inserter(patch), CGAL::Emptyset_iterator(), - CGAL::Polygon_mesh_processing::parameters::weight_calculator - (CGAL::internal::Uniform_weight_fairing(poly)). + CGAL::Polygon_mesh_processing::parameters:: + weight_calculator(CGAL::Weights::Uniform_weight()). density_control_factor(alpha). fairing_continuity(continuity). use_delaunay_triangulation(use_DT))); } else { + auto pmap = get_property_map(CGAL::vertex_point, poly); success = std::get<0>(CGAL::Polygon_mesh_processing::triangulate_refine_and_fair_hole(poly, it, std::back_inserter(patch), CGAL::Emptyset_iterator(), - CGAL::Polygon_mesh_processing::parameters::weight_calculator(CGAL::internal::Cotangent_weight_with_voronoi_area_fairing(poly)). + CGAL::Polygon_mesh_processing::parameters:: + weight_calculator(CGAL::Weights::Secure_cotangent_weight_with_voronoi_area(poly, pmap)). density_control_factor(alpha). fairing_continuity(continuity). use_delaunay_triangulation(use_DT))); diff --git a/Surface_mesh_deformation/include/CGAL/Surface_mesh_deformation.h b/Surface_mesh_deformation/include/CGAL/Surface_mesh_deformation.h index 06a8151e9b7..95956028eb4 100644 --- a/Surface_mesh_deformation/include/CGAL/Surface_mesh_deformation.h +++ b/Surface_mesh_deformation/include/CGAL/Surface_mesh_deformation.h @@ -21,8 +21,8 @@ #include #include #include -#include #include +#include #include #include @@ -54,12 +54,15 @@ enum Deformation_algorithm_tag /// @cond CGAL_DOCUMENT_INTERNAL namespace internal { + template struct Types_selectors; template struct Types_selectors { - typedef internal::Single_cotangent_weight_impl Weight_calculator; + + // Get weight from the weight interface. + typedef CGAL::Weights::Single_cotangent_weight Weight_calculator; struct ARAP_visitor{ template @@ -80,7 +83,9 @@ struct Types_selectors { template struct Types_selectors { - typedef internal::Cotangent_weight_impl Weight_calculator; + + // Get weight from the weight interface. + typedef CGAL::Weights::Cotangent_weight Weight_calculator; typedef typename Types_selectors ::ARAP_visitor ARAP_visitor; @@ -88,7 +93,9 @@ struct Types_selectors { template struct Types_selectors { - typedef internal::Cotangent_weight_impl Weight_calculator; + + // Get weight from the weight interface. + typedef CGAL::Weights::Cotangent_weight Weight_calculator; class ARAP_visitor{ double m_nb_edges_incident; diff --git a/Surface_mesh_deformation/package_info/Surface_mesh_deformation/dependencies b/Surface_mesh_deformation/package_info/Surface_mesh_deformation/dependencies index 974b4a8f5f5..e4abdd72dbd 100644 --- a/Surface_mesh_deformation/package_info/Surface_mesh_deformation/dependencies +++ b/Surface_mesh_deformation/package_info/Surface_mesh_deformation/dependencies @@ -11,7 +11,6 @@ Interval_support Kernel_23 Modular_arithmetic Number_types -Polygon_mesh_processing Profiling_tools Property_map Random_numbers @@ -19,3 +18,4 @@ STL_Extension Solver_interface Stream_support Surface_mesh_deformation +Weights diff --git a/Surface_mesh_parameterization/include/CGAL/Surface_mesh_parameterization/ARAP_parameterizer_3.h b/Surface_mesh_parameterization/include/CGAL/Surface_mesh_parameterization/ARAP_parameterizer_3.h index 3fb7edade43..33fef26c674 100644 --- a/Surface_mesh_parameterization/include/CGAL/Surface_mesh_parameterization/ARAP_parameterizer_3.h +++ b/Surface_mesh_parameterization/include/CGAL/Surface_mesh_parameterization/ARAP_parameterizer_3.h @@ -26,8 +26,8 @@ #include #include #include - #include +#include #include @@ -449,7 +449,7 @@ private: const Point_3& position_vj = get(ppmap, vj); const Point_3& position_vk = get(ppmap, vk); - NT cot = internal::cotangent(position_vi, position_vj, position_vk); + const NT cot = CGAL::Weights::cotangent(position_vi, position_vj, position_vk); put(ctmap, hd, cot); } @@ -489,13 +489,13 @@ private: // coefficient corresponding to the angle at vk if vk is the vertex before vj // while circulating around vi - NT c_k = get(ctmap, opposite(hd, mesh)); + const NT c_k = get(ctmap, opposite(hd, mesh)); // coefficient corresponding to the angle at vl if vl is the vertex after vj // while circulating around vi - NT c_l = get(ctmap, hd); + const NT c_l = get(ctmap, hd); - NT weight = c_k + c_l; + const NT weight = c_k + c_l; return weight; } diff --git a/Surface_mesh_parameterization/include/CGAL/Surface_mesh_parameterization/Barycentric_mapping_parameterizer_3.h b/Surface_mesh_parameterization/include/CGAL/Surface_mesh_parameterization/Barycentric_mapping_parameterizer_3.h index c53c4c56d20..f311cb9396a 100644 --- a/Surface_mesh_parameterization/include/CGAL/Surface_mesh_parameterization/Barycentric_mapping_parameterizer_3.h +++ b/Surface_mesh_parameterization/include/CGAL/Surface_mesh_parameterization/Barycentric_mapping_parameterizer_3.h @@ -19,6 +19,7 @@ #include #include #include +#include #include #include @@ -180,7 +181,7 @@ protected: Vertex_around_target_circulator /* neighbor_vertex_v_j */ ) const { /// In the Tutte Barycentric Mapping algorithm, we have `w_ij = 1`, for `j` neighbor vertex of `i`. - return 1.; + return NT(1); } }; diff --git a/Surface_mesh_parameterization/include/CGAL/Surface_mesh_parameterization/Discrete_authalic_parameterizer_3.h b/Surface_mesh_parameterization/include/CGAL/Surface_mesh_parameterization/Discrete_authalic_parameterizer_3.h index 61ca95aab95..aeb891cff92 100644 --- a/Surface_mesh_parameterization/include/CGAL/Surface_mesh_parameterization/Discrete_authalic_parameterizer_3.h +++ b/Surface_mesh_parameterization/include/CGAL/Surface_mesh_parameterization/Discrete_authalic_parameterizer_3.h @@ -16,11 +16,10 @@ #include -#include #include #include - #include +#include #include @@ -177,34 +176,18 @@ protected: Vertex_around_target_circulator neighbor_vertex_v_j) const { const PPM ppmap = get(vertex_point, mesh); - const Point_3& position_v_i = get(ppmap, main_vertex_v_i); const Point_3& position_v_j = get(ppmap, *neighbor_vertex_v_j); - // Compute the square norm of v_j -> v_i vector - Vector_3 edge = position_v_i - position_v_j; - double square_len = edge*edge; - - // Compute cotangent of (v_k,v_j,v_i) corner (i.e. cotan of v_j corner) - // if v_k is the vertex before v_j when circulating around v_i vertex_around_target_circulator previous_vertex_v_k = neighbor_vertex_v_j; - previous_vertex_v_k--; + --previous_vertex_v_k; const Point_3& position_v_k = get(ppmap, *previous_vertex_v_k); - NT cotg_psi_ij = internal::cotangent(position_v_k, position_v_j, position_v_i); - // Compute cotangent of (v_i,v_j,v_l) corner (i.e. cotan of v_j corner) - // if v_l is the vertex after v_j when circulating around v_i vertex_around_target_circulator next_vertex_v_l = neighbor_vertex_v_j; - next_vertex_v_l++; - const Point_3& position_v_l = get(ppmap,*next_vertex_v_l); - NT cotg_theta_ij = internal::cotangent(position_v_i, position_v_j, position_v_l); + ++next_vertex_v_l; + const Point_3& position_v_l = get(ppmap, *next_vertex_v_l); - NT weight = 0.0; - CGAL_assertion(square_len != 0.0); // two points are identical! - if(square_len != 0.0) - weight = (cotg_psi_ij + cotg_theta_ij) / square_len; - - return weight; + return CGAL::Weights::authalic_weight(position_v_k, position_v_j, position_v_l, position_v_i) / NT(2); } }; diff --git a/Surface_mesh_parameterization/include/CGAL/Surface_mesh_parameterization/Discrete_conformal_map_parameterizer_3.h b/Surface_mesh_parameterization/include/CGAL/Surface_mesh_parameterization/Discrete_conformal_map_parameterizer_3.h index 62bf4f29b3c..3d09fad861c 100644 --- a/Surface_mesh_parameterization/include/CGAL/Surface_mesh_parameterization/Discrete_conformal_map_parameterizer_3.h +++ b/Surface_mesh_parameterization/include/CGAL/Surface_mesh_parameterization/Discrete_conformal_map_parameterizer_3.h @@ -16,11 +16,10 @@ #include -#include #include #include - #include +#include #if defined(CGAL_EIGEN3_ENABLED) #include @@ -173,26 +172,18 @@ protected: Vertex_around_target_circulator neighbor_vertex_v_j) const // its target is main_vertex_v_i { const PPM ppmap = get(vertex_point, mesh); - const Point_3& position_v_i = get(ppmap, main_vertex_v_i); const Point_3& position_v_j = get(ppmap, *neighbor_vertex_v_j); - // Compute cotangent of (v_i,v_k,v_j) corner (i.e. cotan of v_k corner) - // if v_k is the vertex before v_j when circulating around v_i vertex_around_target_circulator previous_vertex_v_k = neighbor_vertex_v_j; - previous_vertex_v_k--; + --previous_vertex_v_k; const Point_3& position_v_k = get(ppmap, *previous_vertex_v_k); - NT cotg_beta_ij = internal::cotangent(position_v_i, position_v_k, position_v_j); - // Compute cotangent of (v_j,v_l,v_i) corner (i.e. cotan of v_l corner) - // if v_l is the vertex after v_j when circulating around v_i vertex_around_target_circulator next_vertex_v_l = neighbor_vertex_v_j; - next_vertex_v_l++; + ++next_vertex_v_l; const Point_3& position_v_l = get(ppmap, *next_vertex_v_l); - NT cotg_alpha_ij = internal::cotangent(position_v_j, position_v_l, position_v_i); - NT weight = cotg_beta_ij + cotg_alpha_ij; - return weight; + return CGAL::Weights::cotangent_weight(position_v_k, position_v_j, position_v_l, position_v_i) / NT(2); } }; diff --git a/Surface_mesh_parameterization/include/CGAL/Surface_mesh_parameterization/Fixed_border_parameterizer_3.h b/Surface_mesh_parameterization/include/CGAL/Surface_mesh_parameterization/Fixed_border_parameterizer_3.h index c3b9d8de956..4222d5f2244 100644 --- a/Surface_mesh_parameterization/include/CGAL/Surface_mesh_parameterization/Fixed_border_parameterizer_3.h +++ b/Surface_mesh_parameterization/include/CGAL/Surface_mesh_parameterization/Fixed_border_parameterizer_3.h @@ -16,10 +16,8 @@ #include -#include #include #include - #include #include diff --git a/Surface_mesh_parameterization/include/CGAL/Surface_mesh_parameterization/Iterative_authalic_parameterizer_3.h b/Surface_mesh_parameterization/include/CGAL/Surface_mesh_parameterization/Iterative_authalic_parameterizer_3.h index 38093049023..2e475b11db8 100644 --- a/Surface_mesh_parameterization/include/CGAL/Surface_mesh_parameterization/Iterative_authalic_parameterizer_3.h +++ b/Surface_mesh_parameterization/include/CGAL/Surface_mesh_parameterization/Iterative_authalic_parameterizer_3.h @@ -21,7 +21,6 @@ #include -#include #include #include #include @@ -34,7 +33,9 @@ #include #include #include -#include +// #include +#include +#include #include #if defined(CGAL_EIGEN3_ENABLED) @@ -629,7 +630,7 @@ private: /// computes `w_ij`, coefficient of matrix `A` for `j` neighbor vertex of `i`. /// - /// \param mesh a triangulated surface. + /// \param tmesh a triangulated surface. /// \param main_vertex_v_i the vertex of `mesh` with index `i` /// \param neighbor_vertex_v_j the vertex of `mesh` with index `j` NT compute_w_ij(const Triangle_mesh& tmesh, @@ -637,36 +638,28 @@ private: Vertex_around_target_circulator neighbor_vertex_v_j) const { const PPM ppmap = get(vertex_point, tmesh); - const PPM_ref position_v_i = get(ppmap, main_vertex_v_i); const PPM_ref position_v_j = get(ppmap, *neighbor_vertex_v_j); - // Compute the square norm of v_j -> v_i vector - Vector_3 edge = position_v_i - position_v_j; - NT square_len = edge*edge; + const Vector_3 edge = position_v_i - position_v_j; + const NT squared_length = edge * edge; - // Compute cotangent of (v_k,v_j,v_i) corner (i.e. cotan of v_j corner) - // if v_k is the vertex before v_j when circulating around v_i vertex_around_target_circulator previous_vertex_v_k = neighbor_vertex_v_j; --previous_vertex_v_k; const PPM_ref position_v_k = get(ppmap, *previous_vertex_v_k); -// NT cotg_psi_ij = internal::cotangent(position_v_k, position_v_j, position_v_i); - NT cotg_beta_ij = internal::cotangent(position_v_i, position_v_k, position_v_j); - // Compute cotangent of (v_i,v_j,v_l) corner (i.e. cotan of v_j corner) - // if v_l is the vertex after v_j when circulating around v_i vertex_around_target_circulator next_vertex_v_l = neighbor_vertex_v_j; ++next_vertex_v_l; + const PPM_ref position_v_l = get(ppmap, *next_vertex_v_l); - const Point_3 position_v_l = get(ppmap, *next_vertex_v_l); -// NT cotg_theta_ij = internal::cotangent(position_v_i, position_v_j, position_v_l); - NT cotg_alpha_ij = internal::cotangent(position_v_j, position_v_l, position_v_i); - - NT weight = 0; - CGAL_assertion(square_len > NT(0)); // two points are identical! - if(square_len != NT(0)) - weight = cotg_beta_ij + cotg_alpha_ij; - + NT weight = NT(0); + CGAL_assertion(squared_length > NT(0)); // two points are identical! + if(squared_length != NT(0)) { + // This version was commented out to be an alternative weight + // in the original code by authors. + // weight = CGAL::Weights::authalic_weight(position_v_k, position_v_j, position_v_l, position_v_i) / NT(2); + weight = CGAL::Weights::cotangent_weight(position_v_k, position_v_j, position_v_l, position_v_i) / NT(2); + } return weight; } @@ -730,7 +723,7 @@ private: VertexIndexMap& vimap) const { auto vpm = get_const_property_map(CGAL::vertex_point, tmesh); - CGAL::internal::Mean_value_weight compute_mvc(tmesh, vpm); + const CGAL::Weights::Edge_tangent_weight compute_mvc(tmesh, vpm); const int i = get(vimap, v); diff --git a/Surface_mesh_parameterization/include/CGAL/Surface_mesh_parameterization/MVC_post_processor_3.h b/Surface_mesh_parameterization/include/CGAL/Surface_mesh_parameterization/MVC_post_processor_3.h index 7803a1980a4..a8053c2a9bd 100644 --- a/Surface_mesh_parameterization/include/CGAL/Surface_mesh_parameterization/MVC_post_processor_3.h +++ b/Surface_mesh_parameterization/include/CGAL/Surface_mesh_parameterization/MVC_post_processor_3.h @@ -21,6 +21,7 @@ #include #include +#include #include #include #include @@ -354,22 +355,6 @@ private: return OK; } - // -> -> - // Return angle (in radians) of of (P,Q,R) corner (i.e. QP,QR angle). - double compute_angle_rad(const Point_2& P, - const Point_2& Q, - const Point_2& R) const - { - Vector_2 u = P - Q; - Vector_2 v = R - Q; - - double angle = std::atan2(v.y(), v.x()) - std::atan2(u.y(), u.x()); - if(angle < 0) - angle += 2 * CGAL_PI; - - return angle; - } - // Fix vertices that are on the convex hull. template @@ -385,23 +370,6 @@ private: } while (++vc != vend); } - NT compute_w_ij_mvc(const Point_2& pi, const Point_2& pj, const Point_2& pk) const - { - // -> -> - // Compute the angle (pj, pi, pk), the angle between the vectors ij and ik - NT angle = compute_angle_rad(pj, pi, pk); - - // For flipped triangles, the connectivity is inversed and thus the angle - // computed by the previous function is not the one we need. Instead, - // we need the explementary angle. - if(angle > CGAL_PI) { // flipped triangle - angle = 2 * CGAL_PI - angle; - } - NT weight = std::tan(0.5 * angle); - - return weight; - } - void fill_linear_system_matrix_mvc_from_points(const Point_2& pi, int i, const Point_2& pj, int j, const Point_2& pk, int k, @@ -418,27 +386,25 @@ private: // The other parts of A(i,j) and A(i,k) will be added when this function // is called from the neighboring faces of F_ijk that share the vertex i - // Compute: - tan(alpha / 2) - NT w_i_base = -1.0 * compute_w_ij_mvc(pi, pj, pk); - // @fixme unefficient: lengths are computed (and inversed!) twice per edge + // Set w_i_base: - tan(alpha / 2) + // Match order of the input points to the new weight implementation. + const Point_2& p = pk; + const Point_2& q = pi; + const Point_2& r = pj; + const CGAL::Weights::Tangent_weight tangent_weight(p, q, r); + // Set w_ij in matrix - Vector_2 edge_ij = pi - pj; - double len_ij = std::sqrt(edge_ij * edge_ij); - CGAL_assertion(len_ij != 0.0); // two points are identical! - NT w_ij = w_i_base / len_ij; + const NT w_ij = tangent_weight.get_w_r(); A.add_coef(i, j, w_ij); // Set w_ik in matrix - Vector_2 edge_ik = pi - pk; - double len_ik = std::sqrt(edge_ik * edge_ik); - CGAL_assertion(len_ik != 0.0); // two points are identical! - NT w_ik = w_i_base / len_ik; + const NT w_ik = tangent_weight.get_w_p(); A.add_coef(i, k, w_ik); // Add to w_ii (w_ii = - sum w_ij) - NT w_ii = - w_ij - w_ik; + const NT w_ii = - w_ij - w_ik; A.add_coef(i, i, w_ii); } diff --git a/Surface_mesh_parameterization/include/CGAL/Surface_mesh_parameterization/Mean_value_coordinates_parameterizer_3.h b/Surface_mesh_parameterization/include/CGAL/Surface_mesh_parameterization/Mean_value_coordinates_parameterizer_3.h index 796504c93f2..95593adf95e 100644 --- a/Surface_mesh_parameterization/include/CGAL/Surface_mesh_parameterization/Mean_value_coordinates_parameterizer_3.h +++ b/Surface_mesh_parameterization/include/CGAL/Surface_mesh_parameterization/Mean_value_coordinates_parameterizer_3.h @@ -17,9 +17,9 @@ #include #include - #include #include +#include #ifdef CGAL_EIGEN3_ENABLED #include @@ -195,35 +195,18 @@ protected: Vertex_around_target_circulator neighbor_vertex_v_j) const { const PPM ppmap = get(vertex_point, mesh); - const Point_3& position_v_i = get(ppmap, main_vertex_v_i); const Point_3& position_v_j = get(ppmap, *neighbor_vertex_v_j); - // Compute the norm of v_j -> v_i vector - Vector_3 edge = position_v_i - position_v_j; - NT len = std::sqrt(edge * edge); - - // Compute angle of (v_j,v_i,v_k) corner (i.e. angle of v_i corner) - // if v_k is the vertex before v_j when circulating around v_i vertex_around_target_circulator previous_vertex_v_k = neighbor_vertex_v_j; - previous_vertex_v_k--; + --previous_vertex_v_k; const Point_3& position_v_k = get(ppmap, *previous_vertex_v_k); - NT gamma_ij = internal::compute_angle_rad(position_v_j, position_v_i, position_v_k); - // Compute angle of (v_l,v_i,v_j) corner (i.e. angle of v_i corner) - // if v_l is the vertex after v_j when circulating around v_i vertex_around_target_circulator next_vertex_v_l = neighbor_vertex_v_j; - next_vertex_v_l++; + ++next_vertex_v_l; const Point_3& position_v_l = get(ppmap, *next_vertex_v_l); - NT delta_ij = internal::compute_angle_rad(position_v_l, position_v_i, position_v_j); - NT weight = 0.0; - CGAL_assertion(len != 0.0); // two points are identical! - if(len != 0.0) - weight = (std::tan(0.5*gamma_ij) + std::tan(0.5*delta_ij)) / len; - CGAL_assertion(weight > 0); - - return weight; + return CGAL::Weights::tangent_weight(position_v_k, position_v_j, position_v_l, position_v_i) / NT(2); } }; diff --git a/Surface_mesh_parameterization/include/CGAL/Surface_mesh_parameterization/Orbifold_Tutte_parameterizer_3.h b/Surface_mesh_parameterization/include/CGAL/Surface_mesh_parameterization/Orbifold_Tutte_parameterizer_3.h index 5d5cc9e6a54..3f8337d0660 100644 --- a/Surface_mesh_parameterization/include/CGAL/Surface_mesh_parameterization/Orbifold_Tutte_parameterizer_3.h +++ b/Surface_mesh_parameterization/include/CGAL/Surface_mesh_parameterization/Orbifold_Tutte_parameterizer_3.h @@ -16,17 +16,15 @@ #include -#include #include #include #include - #include #include #include - +#include +#include #include -#include #include #include @@ -702,17 +700,6 @@ private: CGAL_postcondition(current_line_id_in_M - initial_line_id == number_of_linear_constraints(mesh)); } - // MVC computations - NT compute_w_ij_mvc(const Point_3& pi, const Point_3& pj, const Point_3& pk) const - { - // -> -> - // Compute the angle (pj, pi, pk), the angle between the vectors ij and ik - const NT angle = internal::compute_angle_rad(pj, pi, pk); - const NT weight = std::tan(0.5 * angle); - - return weight; - } - // Computes the coefficients of the mean value Laplacian matrix for the edge. // `ij` in the face `ijk` void fill_mvc_matrix(const Point_3& pi, int i, @@ -730,29 +717,27 @@ private: // The other parts of M(i,j) and M(i,k) will be added when this function // is called from the neighboring faces of F_ijk that share the vertex i - // Compute: - tan(alpha / 2) - const NT w_i_base = 1.0 * compute_w_ij_mvc(pi, pj, pk); - // @fixme unefficient: lengths are computed (and inversed!) twice per edge + // Set w_i_base: - tan(alpha / 2) + const Point_3& p = pk; + const Point_3& q = pi; + const Point_3& r = pj; + const CGAL::Weights::Tangent_weight tangent_weight(p, q, r); + // Set w_ij in matrix - const Vector_3 edge_ij = pi - pj; - const NT len_ij = CGAL::sqrt(edge_ij * edge_ij); - CGAL_assertion(len_ij != 0.0); // two points are identical! - const NT w_ij = w_i_base / len_ij; + const NT w_ij = tangent_weight.get_w_r(); M.add_coef(2*i, 2*j, w_ij); - M.add_coef(2*i +1, 2*j + 1, w_ij); + M.add_coef(2*i + 1, 2*j + 1, w_ij); // Set w_ik in matrix - Vector_3 edge_ik = pi - pk; - const NT len_ik = CGAL::sqrt(edge_ik * edge_ik); - CGAL_assertion(len_ik != 0.0); // two points are identical! - const NT w_ik = w_i_base / len_ik; + const NT w_ik = tangent_weight.get_w_p(); M.add_coef(2*i, 2*k, w_ik); M.add_coef(2*i + 1, 2*k + 1, w_ik); // Add to w_ii (w_ii = - sum w_ij) const NT w_ii = - w_ij - w_ik; + M.add_coef(2*i, 2*i, w_ii); M.add_coef(2*i + 1, 2*i + 1, w_ii); } @@ -790,35 +775,25 @@ private: VertexIndexMap vimap, Matrix& M) const { - const PPM ppmap = get(vertex_point, mesh); + const PPM pmap = get(vertex_point, mesh); + for (const halfedge_descriptor hd : halfedges(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); - - for(halfedge_descriptor hd : halfedges(mesh)) { const vertex_descriptor vi = source(hd, mesh); const vertex_descriptor vj = target(hd, mesh); const int i = get(vimap, vi); const int j = get(vimap, vj); - if(i > j) - continue; - - // times 2 because Cotangent_weight returns 1/2 (cot alpha + cot beta)... - const NT w_ij = 2 * cotan_weight_calculator(hd); + if (i > j) continue; + const CGAL::Weights::Cotangent_weight cotangent_weight; + const NT w_ij = NT(2) * cotangent_weight(hd, mesh, pmap); // ij M.set_coef(2*i, 2*j, w_ij, true /* new coef */); - M.set_coef(2*i +1, 2*j + 1, w_ij, true /* new coef */); + M.set_coef(2*i + 1, 2*j + 1, w_ij, true /* new coef */); // ji M.set_coef(2*j, 2*i, w_ij, true /* new coef */); - M.set_coef(2*j +1, 2*i + 1, w_ij, true /* new coef */); + M.set_coef(2*j + 1, 2*i + 1, w_ij, true /* new coef */); // ii M.add_coef(2*i, 2*i, - w_ij); diff --git a/Surface_mesh_parameterization/include/CGAL/Surface_mesh_parameterization/internal/angles.h b/Surface_mesh_parameterization/include/CGAL/Surface_mesh_parameterization/internal/angles.h deleted file mode 100644 index abd987d33f9..00000000000 --- a/Surface_mesh_parameterization/include/CGAL/Surface_mesh_parameterization/internal/angles.h +++ /dev/null @@ -1,119 +0,0 @@ -// Copyright (c) 2016 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) : Mael Rouxel-Labbé - -#ifndef CGAL_SURFACE_MESH_PARAMETERIZATION_ANGLES_H -#define CGAL_SURFACE_MESH_PARAMETERIZATION_ANGLES_H - -#include - -#include -#include - -#include - -namespace CGAL { - -namespace Surface_mesh_parameterization { - -namespace internal { - -// -> -> -// Returns the cotangent of the corner (P,Q,R) (i.e. the cotan of the angle (QP, QR) ). -template -typename K::FT cotangent(const typename K::Point_3& P, - const typename K::Point_3& Q, - const typename K::Point_3& R) -{ - typedef typename K::FT NT; - typedef typename K::Vector_3 Vector_3; - - Vector_3 u = P - Q; - Vector_3 v = R - Q; - NT dot = (u * v); - Vector_3 cross_vector = CGAL::cross_product(u, v); - NT cross_norm = CGAL::sqrt(cross_vector * cross_vector); - if(cross_norm != NT(0)) - return (dot / cross_norm); - else - return 0; // undefined -} - -// -> -> -// Returns the tangent of the corner (P,Q,R) (i.e. the tangent of angle (QP, QR) ). -template -typename K::FT tangent(const typename K::Point_3& P, - const typename K::Point_3& Q, - const typename K::Point_3& R) -{ - typedef typename K::FT NT; - typedef typename K::Vector_3 Vector_3; - - Vector_3 u = P - Q; - Vector_3 v = R - Q; - NT dot = (u * v); - Vector_3 cross_vector = CGAL::cross_product(u, v); - NT cross_norm = CGAL::sqrt(cross_vector * cross_vector); - if(dot != NT(0)) - return (cross_norm / dot); - else - return 0; // undefined -} - -// Fixes the sine to be within [-1;1]. -template -typename K::FT fix_sine(typename K::FT sine) -{ - if(sine >= 1) - return 1; - else if(sine <= -1) - return -1; - else - return sine; -} - -// -> -> -// Returns the angle (in radians) of the corner (P,Q,R) (i.e. the angle (QP, QR) ). -template -typename K::FT compute_angle_rad(const typename K::Point_3& P, - const typename K::Point_3& Q, - const typename K::Point_3& R) -{ - typedef typename K::FT NT; - typedef typename K::Vector_3 Vector_3; - - Vector_3 u = P - Q; - Vector_3 v = R - Q; - - NT product = CGAL::sqrt(u * u) * CGAL::sqrt(v * v); - if(product == NT(0)) - return 0; - - // cosine - NT dot = (u * v); - NT cosine = dot / product; - - // sine - Vector_3 w = CGAL::cross_product(u, v); - NT abs_sine = CGAL::sqrt(w * w) / product; - - if(cosine >= NT(0)) - return std::asin(fix_sine(abs_sine)); - else - return CGAL_PI - std::asin(fix_sine(abs_sine)); -} - -} // namespace internal - -} // namespace Surface_mesh_parameterization - -} // namespace CGAL - -#endif // CGAL_SURFACE_MESH_PARAMETERIZATION_ANGLES_H diff --git a/Surface_mesh_parameterization/include/CGAL/Surface_mesh_parameterization/orbifold_shortest_path.h b/Surface_mesh_parameterization/include/CGAL/Surface_mesh_parameterization/orbifold_shortest_path.h index 4243374f9d3..4b1e4ab2cc0 100644 --- a/Surface_mesh_parameterization/include/CGAL/Surface_mesh_parameterization/orbifold_shortest_path.h +++ b/Surface_mesh_parameterization/include/CGAL/Surface_mesh_parameterization/orbifold_shortest_path.h @@ -190,8 +190,8 @@ void compute_shortest_paths_between_cones(const TriangleMesh& mesh, compute_shortest_paths_between_two_cones(mesh, *first, *next, std::back_inserter(seams)); } - std::ofstream out("shortest_path.selection.txt"); #ifdef CGAL_SMP_ORBIFOLD_DEBUG + std::ofstream out("shortest_path.selection.txt"); internal::output_shortest_paths_to_selection_file(mesh, seams, out); #endif } diff --git a/Surface_mesh_parameterization/package_info/Surface_mesh_parameterization/dependencies b/Surface_mesh_parameterization/package_info/Surface_mesh_parameterization/dependencies index 6dedee5435a..0d5931b478b 100644 --- a/Surface_mesh_parameterization/package_info/Surface_mesh_parameterization/dependencies +++ b/Surface_mesh_parameterization/package_info/Surface_mesh_parameterization/dependencies @@ -26,3 +26,4 @@ Stream_support Surface_mesh_parameterization TDS_2 Triangulation_2 +Weights diff --git a/Surface_mesh_skeletonization/include/CGAL/Mean_curvature_flow_skeletonization.h b/Surface_mesh_skeletonization/include/CGAL/Mean_curvature_flow_skeletonization.h index 81bacc0d8a8..3882180786a 100644 --- a/Surface_mesh_skeletonization/include/CGAL/Mean_curvature_flow_skeletonization.h +++ b/Surface_mesh_skeletonization/include/CGAL/Mean_curvature_flow_skeletonization.h @@ -35,7 +35,7 @@ #include // Compute cotangent Laplacian -#include +#include // Compute the vertex normal #include @@ -222,12 +222,8 @@ public: typedef typename boost::graph_traits::edge_descriptor edge_descriptor; typedef typename boost::graph_traits::edge_iterator edge_iterator; - // Cotangent weight calculator - typedef internal::Cotangent_weight::type, - internal::Cotangent_value_minimum_zero::type, - internal::Cotangent_value_Meyer_secure > > Weight_calculator; + // Get weight from the weight interface. + typedef CGAL::Weights::Cotangent_weight Weight_calculator; typedef internal::Curve_skeleton Iterative_authalic_parameterizer_3.h - template< - typename PolygonMesh, - typename VertexPointMap = typename boost::property_map::type> - class Mean_value_weight { - - using GeomTraits = typename CGAL::Kernel_traits< - typename boost::property_traits::value_type>::type; - using FT = typename GeomTraits::FT; - - const PolygonMesh& m_pmesh; - const VertexPointMap m_pmap; - const GeomTraits m_traits; - - public: - using vertex_descriptor = typename boost::graph_traits::vertex_descriptor; - using halfedge_descriptor = typename boost::graph_traits::halfedge_descriptor; - - Mean_value_weight(const PolygonMesh& pmesh, const VertexPointMap pmap) : - m_pmesh(pmesh), m_pmap(pmap), m_traits() { } - - // Returns the mean-value coordinate of the specified halfedge_descriptor. - // Returns different values for different edge orientations (which is normal - // behavior according to the formula). - FT operator()(const halfedge_descriptor he) const { - - const vertex_descriptor v0 = target(he, m_pmesh); - const vertex_descriptor v1 = source(he, m_pmesh); - const auto& p0 = get(m_pmap, v0); - const auto& p1 = get(m_pmap, v1); - const FT norm = internal::distance_3(m_traits, p0, p1); - - // Only one triangle for border edges. - if (is_border_edge(he, m_pmesh)) { - const halfedge_descriptor he_cw = opposite(next(he, m_pmesh), m_pmesh); - vertex_descriptor v2 = source(he_cw, m_pmesh); - if (is_border_edge(he_cw, m_pmesh)) { - const halfedge_descriptor he_ccw = prev(opposite(he, m_pmesh), m_pmesh); - v2 = source(he_ccw, m_pmesh); - } - return half_tan_value_2(v1, v0, v2) / norm; - } else { - const halfedge_descriptor he_cw = opposite(next(he, m_pmesh), m_pmesh); - const vertex_descriptor v2 = source(he_cw, m_pmesh); - const halfedge_descriptor he_ccw = prev(opposite(he, m_pmesh), m_pmesh); - const vertex_descriptor v3 = source(he_ccw, m_pmesh); - return ( - half_tan_value_2(v1, v0, v2) / norm + - half_tan_value_2(v1, v0, v3) / norm ); - } - } - - private: - // The authors deviation built on Meyer_02. - FT half_tan_value_2( - const vertex_descriptor v0, - const vertex_descriptor v1, - const vertex_descriptor v2) const { - - using Get_sqrt = internal::Get_sqrt; - const auto sqrt = Get_sqrt::sqrt_object(m_traits); - - const auto squared_length_3 = - m_traits.compute_squared_length_3_object(); - const auto construct_vector_3 = - m_traits.construct_vector_3_object(); - - const auto& p0 = get(m_pmap, v0); - const auto& p1 = get(m_pmap, v1); - const auto& p2 = get(m_pmap, v2); - - const auto a = construct_vector_3(p1, p0); - const auto b = construct_vector_3(p1, p2); - - const FT dot_ab = a[0] * b[0] + a[1] * b[1] + a[2] * b[2]; - const FT dot_aa = squared_length_3(a); - const FT dot_bb = squared_length_3(b); - const FT dot_aa_bb = dot_aa * dot_bb; - - const FT cos_rep = dot_ab; - const FT sin_rep = sqrt(dot_aa_bb - dot_ab * dot_ab); - const FT normalizer = sqrt(dot_aa_bb); // |a| * |b| - - // The formula from [Floater04] page 4: - // tan(Q / 2) = (1 - cos(Q)) / sin(Q). - return (normalizer - cos_rep) / sin_rep; - } - }; - /// \endcond /*! diff --git a/Weights/include/CGAL/Weights/tangent_weights.h b/Weights/include/CGAL/Weights/tangent_weights.h index a7e93109a02..47a7539ae8a 100644 --- a/Weights/include/CGAL/Weights/tangent_weights.h +++ b/Weights/include/CGAL/Weights/tangent_weights.h @@ -372,6 +372,67 @@ namespace Weights { return tangent_weight(t, r, p, q, traits); } + // Undocumented tangent weight class. + // Its constructor takes a polygon mesh and a vertex to point map + // and its operator() is defined based on the halfedge_descriptor only. + // This version is currently used in: + // Surface_mesh_parameterizer -> Iterative_authalic_parameterizer_3.h + template< + typename PolygonMesh, + typename VertexPointMap = typename boost::property_map::type> + class Edge_tangent_weight { + + using GeomTraits = typename CGAL::Kernel_traits< + typename boost::property_traits::value_type>::type; + using FT = typename GeomTraits::FT; + + const PolygonMesh& m_pmesh; + const VertexPointMap m_pmap; + const GeomTraits m_traits; + + public: + using vertex_descriptor = typename boost::graph_traits::vertex_descriptor; + using halfedge_descriptor = typename boost::graph_traits::halfedge_descriptor; + + Edge_tangent_weight(const PolygonMesh& pmesh, const VertexPointMap pmap) : + m_pmesh(pmesh), m_pmap(pmap), m_traits() { } + + FT operator()(const halfedge_descriptor he) const { + + FT weight = FT(0); + if (is_border_edge(he, m_pmesh)) { + const auto h1 = next(he, m_pmesh); + + const auto v0 = target(he, m_pmesh); + const auto v1 = source(he, m_pmesh); + const auto v2 = target(h1, m_pmesh); + + const auto& p0 = get(m_pmap, v0); + const auto& p1 = get(m_pmap, v1); + const auto& p2 = get(m_pmap, v2); + + weight = internal::tangent_3(m_traits, p0, p2, p1); + + } else { + const auto h1 = next(he, m_pmesh); + const auto h2 = prev(opposite(he, m_pmesh), m_pmesh); + + const auto v0 = target(he, m_pmesh); + const auto v1 = source(he, m_pmesh); + const auto v2 = target(h1, m_pmesh); + const auto v3 = source(h2, m_pmesh); + + const auto& p0 = get(m_pmap, v0); + const auto& p1 = get(m_pmap, v1); + const auto& p2 = get(m_pmap, v2); + const auto& p3 = get(m_pmap, v3); + + weight = tangent_weight(p2, p1, p3, p0) / FT(2); + } + return weight; + } + }; + // Undocumented tangent weight class. // Its constructor takes three points either in 2D or 3D. // This version is currently used in: diff --git a/Weights/include/CGAL/Weights/utils.h b/Weights/include/CGAL/Weights/utils.h index 919b2fc2bf2..4f68b34e1ef 100644 --- a/Weights/include/CGAL/Weights/utils.h +++ b/Weights/include/CGAL/Weights/utils.h @@ -18,7 +18,6 @@ // Internal includes. #include -#include namespace CGAL { namespace Weights {