songsenand
/
cgal
mirror of https://github.com/CGAL/cgal
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

updated weights in the pmp package

This commit is contained in:
Dmitry Anisimov 2021-06-03 11:33:18 +02:00
parent 277b804eab
commit 2a6aaddd8f
8 changed files with 61 additions and 991 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

917
Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/Weights.h
View File

@ -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 <CGAL/license/Polygon_mesh_processing/core.h>
/// @cond CGAL_DOCUMENT_INTERNAL
#include <CGAL/boost/graph/helpers.h>
#include <CGAL/Kernel/global_functions_3.h>
#include <CGAL/property_map.h>
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<class PolygonMesh>
struct Cotangent_value_Meyer_impl
{
typedef typename boost::graph_traits<PolygonMesh>::vertex_descriptor vertex_descriptor;
template <class VertexPointMap>
double operator()(vertex_descriptor v0, vertex_descriptor v1, vertex_descriptor v2, const VertexPointMap& ppmap)
{
typedef typename Kernel_traits<
typename boost::property_traits<VertexPointMap>::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<double>::max)() :
-(std::numeric_limits<double>::max)();
}
return dot_ab / divider;
}
};
// Same as above but with a different API
template<class PolygonMesh
, class VertexPointMap = typename boost::property_map<PolygonMesh, CGAL::vertex_point_t>::type>
class Cotangent_value_Meyer
{
protected:
typedef typename boost::graph_traits<PolygonMesh>::vertex_descriptor vertex_descriptor;
typedef VertexPointMap Point_property_map;
typedef typename boost::property_traits<Point_property_map>::value_type Point;
typedef typename Kernel_traits<Point>::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<PolygonMesh>()(v0, v1, v2, ppmap());
}
};
// imported from skeletonization
template<class PolygonMesh,
class VertexPointMap =
typename boost::property_map<PolygonMesh, CGAL::vertex_point_t>::type>
class Cotangent_value_Meyer_secure
{
typedef typename boost::graph_traits<PolygonMesh>::vertex_descriptor vertex_descriptor;
typedef VertexPointMap Point_property_map;
typedef typename boost::property_traits<Point_property_map>::value_type Point;
typedef typename Kernel_traits<Point>::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<class PolygonMesh
, class VertexPointMap = typename boost::property_map<PolygonMesh, CGAL::vertex_point_t>::type
, class CotangentValue = Cotangent_value_Meyer<PolygonMesh, VertexPointMap> >
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<PolygonMesh>::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<class PolygonMesh
, class VertexPointMap = typename boost::property_map<PolygonMesh, CGAL::vertex_point_t>::type
, class CotangentValue = Cotangent_value_Meyer<PolygonMesh, VertexPointMap> >
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<PolygonMesh>::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<class PolygonMesh
, class VertexPointMap = typename boost::property_map<PolygonMesh, CGAL::vertex_point_t>::type
, class CotangentValue = Cotangent_value_Meyer<PolygonMesh, VertexPointMap> >
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<PolygonMesh>::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<class PolygonMesh,
class CotangentValue = Cotangent_value_Meyer_impl<PolygonMesh> >
struct Cotangent_value_minimum_zero_impl : CotangentValue
{
typedef typename boost::graph_traits<PolygonMesh>::vertex_descriptor vertex_descriptor;
template <class VertexPointMap>
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<class PolygonMesh
, class VertexPointMap = typename boost::property_map<PolygonMesh, CGAL::vertex_point_t>::type
, class CotangentValue
= Cotangent_value_Meyer<PolygonMesh, VertexPointMap> >
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<PolygonMesh>::vertex_descriptor vertex_descriptor;
typedef typename boost::graph_traits<PolygonMesh>::in_edge_iterator in_edge_iterator;
typedef typename boost::graph_traits<PolygonMesh>::halfedge_descriptor halfedge_descriptor;
typedef VertexPointMap Point_property_map;
typedef typename boost::property_traits<Point_property_map>::value_type Point;
typedef typename Kernel_traits<Point>::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<class PolygonMesh
, class VertexPointMap = typename boost::property_map<PolygonMesh, CGAL::vertex_point_t>::type
, class CotangentValue = Cotangent_value_Meyer<PolygonMesh, VertexPointMap> >
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<PolygonMesh>::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<PolygonMesh> >
struct Cotangent_weight_impl : CotangentValue
{
typedef typename boost::graph_traits<PolygonMesh>::halfedge_descriptor halfedge_descriptor;
typedef typename boost::graph_traits<PolygonMesh>::vertex_descriptor vertex_descriptor;
// Returns the cotangent weight of specified halfedge_descriptor
// Edge orientation is trivial
template<class VertexPointMap>
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<class PolygonMesh
, class VertexPointMap = typename boost::property_map<PolygonMesh, vertex_point_t>::type
, class CotangentValue
= Cotangent_value_minimum_zero<PolygonMesh, VertexPointMap> >
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<PolygonMesh>::halfedge_descriptor halfedge_descriptor;
typedef typename boost::graph_traits<PolygonMesh>::vertex_descriptor vertex_descriptor;
typedef typename boost::property_map<PolygonMesh,vertex_point_t>::type Point_property_map;
typedef typename boost::property_traits<Point_property_map>::value_type Point;
typedef typename Kernel_traits<Point>::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<class PolygonMesh,
class CotangentValue = Cotangent_value_Meyer_impl<PolygonMesh> >
struct Single_cotangent_weight_impl : CotangentValue
{
typedef typename boost::graph_traits<PolygonMesh>::halfedge_descriptor halfedge_descriptor;
typedef typename boost::graph_traits<PolygonMesh>::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 <class VertexPointMap>
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<class PolygonMesh
, class VertexPointMap = typename boost::property_map<PolygonMesh, CGAL::vertex_point_t>::type
, class CotangentValue = Cotangent_value_Meyer<PolygonMesh, VertexPointMap> >
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<PolygonMesh>::halfedge_descriptor halfedge_descriptor;
typedef typename boost::graph_traits<PolygonMesh>::vertex_descriptor vertex_descriptor;
typedef typename boost::property_map<PolygonMesh,vertex_point_t>::type Point_property_map;
typedef typename boost::property_traits<Point_property_map>::value_type Point;
typedef typename Kernel_traits<Point>::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<class PolygonMesh
, class VertexPointMap = typename boost::property_map<PolygonMesh, vertex_point_t>::type
, class CotangentValue = Cotangent_value_Meyer<PolygonMesh, VertexPointMap>
>
class Cotangent_weight_with_triangle_area : CotangentValue
{
typedef PolygonMesh PM;
typedef VertexPointMap VPMap;
typedef typename boost::property_traits<VPMap>::value_type Point;
typedef typename boost::graph_traits<PM>::halfedge_descriptor halfedge_descriptor;
typedef typename boost::graph_traits<PM>::vertex_descriptor vertex_descriptor;
Cotangent_weight_with_triangle_area()
{}
public:
Cotangent_weight_with_triangle_area(PolygonMesh& pmesh_, VertexPointMap vpmap_)
: CotangentValue(pmesh_, vpmap_)
{}
PolygonMesh& pmesh()
{
return CotangentValue::pmesh();
}
VertexPointMap& ppmap()
{
return CotangentValue::ppmap();
}
double operator()(halfedge_descriptor he)
{
vertex_descriptor v0 = target(he, pmesh());
vertex_descriptor v1 = source(he, pmesh());
// Only one triangle for border edges
if (is_border_edge(he, pmesh()))
{
halfedge_descriptor he_cw = opposite( next(he, pmesh()) , pmesh() );
vertex_descriptor v2 = source(he_cw, pmesh());
if (is_border_edge(he_cw, pmesh()) )
{
halfedge_descriptor he_ccw = prev( opposite(he, pmesh()) , pmesh() );
v2 = source(he_ccw, pmesh());
}
const Point& v0_p = get(ppmap(), v0);
const Point& v1_p = get(ppmap(), v1);
const Point& v2_p = get(ppmap(), v2);
double area_t = to_double(CGAL::sqrt(squared_area(v0_p, v1_p, v2_p)));
return ( CotangentValue::operator()(v0, v2, v1) / area_t );
}
else
{
halfedge_descriptor he_cw = opposite( next(he, pmesh()) , pmesh() );
vertex_descriptor v2 = source(he_cw, pmesh());
halfedge_descriptor he_ccw = prev( opposite(he, pmesh()) , pmesh() );
vertex_descriptor v3 = source(he_ccw, pmesh());
const Point& v0_p = get(ppmap(), v0);
const Point& v1_p = get(ppmap(), v1);
const Point& v2_p = get(ppmap(), v2);
const Point& v3_p = get(ppmap(), v3);
double area_t1 = to_double(CGAL::sqrt(squared_area(v0_p, v1_p, v2_p)));
double area_t2 = to_double(CGAL::sqrt(squared_area(v0_p, v1_p, v3_p)));
return ( CotangentValue::operator()(v0, v2, v1) / area_t1
+ CotangentValue::operator()(v0, v3, v1) / area_t2);
}
return 0.;
}
};
// Mean value calculator described in -[Floater04] Mean Value Coordinates-
template<class PolygonMesh
, class VertexPointMap = typename boost::property_map<PolygonMesh, CGAL::vertex_point_t>::type>
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<PolygonMesh>::halfedge_descriptor halfedge_descriptor;
typedef typename boost::graph_traits<PolygonMesh>::vertex_descriptor vertex_descriptor;
typedef VertexPointMap Point_property_map;
typedef typename boost::property_traits<Point_property_map>::value_type Point;
typedef typename Kernel_traits<Point>::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<PolygonMesh>,
class SecondaryWeight = Mean_value_weight<PolygonMesh> >
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<PolygonMesh>::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 PolygonMesh>
class Uniform_weight
{
public:
typedef typename boost::graph_traits<PolygonMesh>::halfedge_descriptor halfedge_descriptor;
double operator()(halfedge_descriptor /*e*/)
{ return 1.0; }
};
////////////////////////////////////////////////////////////////////////////
// FAIRING //
template<class PolygonMesh>
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<PolygonMesh>::halfedge_descriptor halfedge_descriptor;
typedef typename boost::graph_traits<PolygonMesh>::vertex_descriptor vertex_descriptor;
typedef typename boost::property_map<PolygonMesh,vertex_point_t>::type Point_property_map;
typedef typename boost::property_traits<Point_property_map>::value_type Point;
typedef typename Kernel_traits<Point>::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<double>::max)();
}
return 1.0 / divider;
}
};
template<class PolygonMesh
, class VertexPointMap = typename boost::property_map<PolygonMesh, vertex_point_t>::type
>
class Cotangent_weight_with_voronoi_area_fairing
{
typedef PolygonMesh PM;
typedef VertexPointMap VPMap;
Voronoi_area<PM, VPMap> voronoi_functor;
Cotangent_weight<PM, VPMap, Cotangent_value_Meyer<PM, VPMap> > 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<PM>::halfedge_descriptor halfedge_descriptor;
typedef typename boost::graph_traits<PM>::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<PolygonMesh, vertex_point_t>::type>
class Cotangent_weight_with_voronoi_area_fairing_secure {
typedef PolygonMesh PM;
typedef VertexPointMap VPMap;
Voronoi_area<PM, VPMap> voronoi_functor;
Cotangent_weight<PM, VPMap, Cotangent_value_Meyer_secure<PM, VPMap> > 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<PM>::halfedge_descriptor halfedge_descriptor;
typedef typename boost::graph_traits<PM>::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 PolygonMesh>
class Uniform_weight_fairing
{
public:
typedef typename boost::graph_traits<PolygonMesh>::halfedge_descriptor halfedge_descriptor;
typedef typename boost::graph_traits<PolygonMesh>::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

5
Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/fair.h
View File

@ -20,6 +20,7 @@
#include <CGAL/Polygon_mesh_processing/internal/fair_impl.h>
#include <CGAL/Polygon_mesh_processing/internal/named_function_params.h>
#include <CGAL/Polygon_mesh_processing/internal/named_params_helper.h>
#include <CGAL/Weights/internal/tools.h>
#if defined(CGAL_EIGEN3_ENABLED)
#include <CGAL/Eigen_solver_traits.h> // 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<TriangleMesh, VPMap>
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));

45
Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Smoothing/curvature_flow_impl.h
View File

@ -17,10 +17,9 @@
#include <CGAL/license/Polygon_mesh_processing/meshing_hole_filling.h>
#include <CGAL/Polygon_mesh_processing/measure.h>
#include <CGAL/Polygon_mesh_processing/Weights.h>
#include <CGAL/Polygon_mesh_processing/internal/named_function_params.h>
#include <CGAL/Polygon_mesh_processing/internal/named_params_helper.h>
#include <CGAL/Weights/internal/tools.h>
#include <CGAL/Dynamic_property_map.h>
#include <CGAL/utility.h>
@ -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<typename TriangleMesh,
typename VertexPointMap,
typename CotangentValue = CGAL::internal::Cotangent_value_Meyer<TriangleMesh, VertexPointMap> >
struct Edge_cotangent_weight
: public CotangentValue
{
typedef typename boost::graph_traits<TriangleMesh>::halfedge_descriptor halfedge_descriptor;
typedef typename boost::graph_traits<TriangleMesh>::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<typename TriangleMesh,
typename VertexPointMap,
typename VertexConstraintMap,
@ -409,7 +369,8 @@ private:
std::vector<bool> constrained_flags_;
const GeomTraits& traits_;
Edge_cotangent_weight<TriangleMesh, VertexPointMap> weight_calculator_;
CGAL::Weights::internal::Edge_cotangent_weight_wrapper<
TriangleMesh, VertexPointMap> weight_calculator_;
};
} // internal

7
Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/fair_impl.h
View File

@ -19,7 +19,6 @@
#include <map>
#include <set>
#include <CGAL/assertions.h>
#include <CGAL/Polygon_mesh_processing/Weights.h>
#ifdef CGAL_PMP_FAIR_DEBUG
#include <CGAL/Timer.h>
#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);

18
Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/repair_self_intersections.h
View File

@ -111,7 +111,7 @@ FaceOutputIterator replace_faces_with_patch(const std::vector<typename boost::gr
std::map<Point, vertex_descriptor> 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<typename boost::graph_tra
}
std::vector<std::vector<Point> > 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<Point>{get(local_vpm, target(h, local_mesh)),
@ -702,7 +702,7 @@ bool construct_tentative_hole_patch(std::vector<typename boost::graph_traits<Tri
// Collect vertices and edges inside the current selection cc: first collect all vertices and
// edges incident to the faces to remove...
for(const face_descriptor f : cc_faces)
for(const face_descriptor& f : cc_faces)
{
for(halfedge_descriptor h : halfedges_around_face(halfedge(f, tmesh), tmesh))
{
@ -843,7 +843,7 @@ bool construct_tentative_hole_patch(std::vector<typename boost::graph_traits<Tri
hole_points.reserve(cc_border_hedges.size());
third_points.reserve(cc_border_hedges.size());
for(const halfedge_descriptor h : cc_border_hedges)
for(const halfedge_descriptor& h : cc_border_hedges)
{
const vertex_descriptor v = source(h, tmesh);
hole_points.push_back(get(vpm, v));
@ -887,7 +887,7 @@ bool construct_tentative_sub_hole_patch(std::vector<std::vector<typename boost::
// (pointing inside the domain to be remeshed)
std::set<halfedge_descriptor> internal_hedges;
std::vector<halfedge_descriptor> 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<std::vector<typename boost::
hole_points.reserve(cc_border_hedges.size());
third_points.reserve(cc_border_hedges.size());
for(const halfedge_descriptor h : cc_border_hedges)
for(const halfedge_descriptor& h : cc_border_hedges)
{
const vertex_descriptor v = source(h, tmesh);
hole_points.push_back(get(vpm, v));
@ -1091,7 +1091,7 @@ bool fill_hole(std::vector<typename boost::graph_traits<TriangleMesh>::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<typename boost::graph_traits<Triangl
CGAL_assertion(!does_self_intersect(new_faces, tmesh, parameters::vertex_point_map(vpm)));
// Update working range with the new faces
for(const face_descriptor f : cc_faces)
for(const face_descriptor& f : cc_faces)
working_face_range.erase(f);
working_face_range.insert(new_faces.begin(), new_faces.end());
@ -1582,7 +1582,7 @@ remove_self_intersections_one_step(std::set<typename boost::graph_traits<Triangl
#endif
// remove faces from the set to process
for(const face_descriptor f : cc_faces)
for(const face_descriptor& f : cc_faces)
faces_to_remove.erase(f);
#ifdef CGAL_PMP_REMOVE_SELF_INTERSECTION_OUTPUT

25
Polygon_mesh_processing/test/Polygon_mesh_processing/triangulate_hole_Polyhedron_3_no_delaunay_test.cpp
View File

@ -25,6 +25,19 @@
typedef CGAL::Exact_predicates_inexact_constructions_kernel Epic;
typedef CGAL::Exact_predicates_exact_constructions_kernel Epec;
template<
class GeomTraits,
class PolygonMesh>
class Uniform_weight_fairing {
public:
using FT = typename GeomTraits::FT;
using halfedge_descriptor = typename boost::graph_traits<PolygonMesh>::halfedge_descriptor;
using vertex_descriptor = typename boost::graph_traits<PolygonMesh>::vertex_descriptor;
FT w_ij(halfedge_descriptor) { return FT(1); }
FT w_i(vertex_descriptor) { return FT(1); }
};
template <class Polyhedron>
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 <class Polyhedron>
template <class Kernel, class Polyhedron>
void test_triangulate_refine_and_fair_hole_compile() {
typedef typename boost::graph_traits<Polyhedron>::halfedge_descriptor Halfedge_handle;
typedef typename boost::graph_traits<Polyhedron>::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<Polyhedron>(poly)).
CGAL::Polygon_mesh_processing::parameters::
weight_calculator(Uniform_weight_fairing<Kernel, Polyhedron>()).
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<Polyhedron>(poly)));
CGAL::Polygon_mesh_processing::parameters::
weight_calculator(Uniform_weight_fairing<Kernel, Polyhedron>()));
// default solver and weight
read_poly_with_borders("elephant_quad_hole.off", poly, border_reps);
@ -404,7 +417,7 @@ typedef CGAL::Surface_mesh<typename Kernel::Point_3> Polyhedron;
test_triangulate_hole_should_be_no_output<Polyhedron>("data/non_manifold_vertex.off");
test_triangulate_hole_should_be_no_output<Polyhedron>("data/two_tris_collinear.off");
test_triangulate_refine_and_fair_hole_compile<Polyhedron>();
test_triangulate_refine_and_fair_hole_compile<Kernel, Polyhedron>();
}

25
Polygon_mesh_processing/test/Polygon_mesh_processing/triangulate_hole_Polyhedron_3_test.cpp
View File

@ -30,6 +30,20 @@ typedef boost::graph_traits<Polyhedron>::halfedge_iterator Halfedge_iterator;
typedef CGAL::Halfedge_around_face_circulator<Polyhedron> Halfedge_around_facet_circulator;
typedef boost::property_map<Polyhedron,CGAL::vertex_point_t>::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<PolygonMesh>::halfedge_descriptor;
using vertex_descriptor = typename boost::graph_traits<PolygonMesh>::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<Polyhedron>(poly)).
sparse_linear_solver(Default_solver()).use_2d_constrained_delaunay_triangulation(false));
CGAL::Polygon_mesh_processing::parameters::
weight_calculator(Uniform_weight_fairing<Kernel, Polyhedron>()).
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<Polyhedron>(poly)));
CGAL::Polygon_mesh_processing::parameters::
weight_calculator(Uniform_weight_fairing<Kernel, Polyhedron>()));
// default solver and weight
read_poly_with_borders("elephant_quad_hole.off", poly, border_reps);

10
Weights/include/CGAL/Weights/internal/tools.h
View File

@ -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;