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

Merge branch 'releases/CGAL-4.11-branch'

This commit is contained in:
Laurent Rineau 2017-12-19 16:29:10 +01:00
parent 8f3e229802 5181cd7474
commit 1c2e9df8cf
28 changed files with 442 additions and 509 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
BGL/doc/BGL/Concepts/HalfedgeGraph.h
View File

@ -49,7 +49,7 @@ class HalfedgeGraph {
};
/*! \relates HalfedgeGraph
returns the edge corresponding to halfedges `h` and `opposite(h,g)`.
returns the edge corresponding to halfedges `h` and `opposite(h,g)`, with the following invariant `halfedge(edge(h,g),g)==h`.
*/
template <typename HalfedgeGraph>
boost::graph_traits<HalfedgeGraph>::edge_descriptor

9
Convex_hull_3/include/CGAL/Convex_hull_3/dual/halfspace_intersection_3.h
View File

@ -67,8 +67,10 @@ namespace CGAL
typedef typename Polyhedron_dual::Vertex_const_iterator
Vertex_const_iterator;
// typedefs for primal
// typedef and type for primal
typedef typename boost::graph_traits<Polyhedron>::vertex_descriptor vertex_descriptor;
typename boost::property_map<Polyhedron, vertex_point_t>::type vpm =
get(CGAL::vertex_point, primal);
// Typedefs for intersection
typedef typename Kernel::Plane_3 Plane_3;
@ -115,8 +117,9 @@ namespace CGAL
origin.y() + pp->y(),
origin.z() + pp->z());
primal_vertices[it] = add_vertex(ppp, primal);
vertex_descriptor vd = add_vertex(primal);
primal_vertices[it] = vd;
put(vpm, vd, ppp);
}
// Then, add facets to the primal polyhedron

13
Convex_hull_3/include/CGAL/Convex_hull_3/dual/halfspace_intersection_with_constructions_3.h
View File

@ -56,22 +56,25 @@ namespace CGAL
typedef typename boost::graph_traits<Polyhedron>::face_descriptor face_descriptor;
typedef typename boost::graph_traits<Polyhedron>::vertex_descriptor vertex_descriptor;
typename boost::property_map<Polyhedron, vertex_point_t>::const_type vpmap = get(CGAL::vertex_point, primal);
typename boost::property_map<Polyhedron, vertex_point_t>::const_type vpm_primal = get(CGAL::vertex_point, primal);
typename boost::property_map<Polyhedron, vertex_point_t>::type vpm_dual = get(CGAL::vertex_point, dual);
// compute coordinates of extreme vertices in the dual polyhedron
// from primal faces
boost::unordered_map<face_descriptor, vertex_descriptor> extreme_points;
BOOST_FOREACH (face_descriptor fd , faces( primal)){
halfedge_descriptor h = halfedge(fd,primal);
Plane_3 p (get(vpmap, target(h, primal)),
get(vpmap, target(next(h, primal), primal)),
get(vpmap, target(next(next(h, primal), primal), primal)));
Plane_3 p (get(vpm_primal, target(h, primal)),
get(vpm_primal, target(next(h, primal), primal)),
get(vpm_primal, target(next(next(h, primal), primal), primal)));
// translate extreme vertex
Point_3 extreme_p = CGAL::ORIGIN + p.orthogonal_vector () / (-p.d());
Point_3 translated_extreme_p(extreme_p.x() + origin.x(),
extreme_p.y() + origin.y(),
extreme_p.z() + origin.z());
extreme_points[fd] = add_vertex(translated_extreme_p,dual);
vertex_descriptor vd = add_vertex(dual);
extreme_points[fd] = vd;
put(vpm_dual, vd, translated_extreme_p);
}
// build faces

17
Convex_hull_3/include/CGAL/convex_hull_3.h
View File

@ -316,13 +316,26 @@ void coplanar_3_hull(InputIterator first, InputIterator beyond,
typename boost::property_map<Polyhedron_3, CGAL::vertex_point_t>::type vpm
= get(CGAL::vertex_point, P);
std::vector<typename boost::graph_traits<Polyhedron_3>::vertex_descriptor> vertices;
typedef boost::graph_traits<Polyhedron_3> Graph_traits;
typedef typename Graph_traits::vertex_descriptor vertex_descriptor;
typedef typename Graph_traits::halfedge_descriptor halfedge_descriptor;
typedef typename Graph_traits::face_descriptor face_descriptor;
std::vector<vertex_descriptor> vertices;
vertices.reserve(CH_2.size());
BOOST_FOREACH(const Point_3& p, CH_2){
vertices.push_back(add_vertex(P));
put(vpm, vertices.back(),p);
}
Euler::add_face(vertices, P);
face_descriptor f = Euler::add_face(vertices, P);
// Then triangulate that face
const halfedge_descriptor he = halfedge(f, P);
halfedge_descriptor other_he = next(next(he, P), P);
for(std::size_t i = 3, end = vertices.size(); i < end; ++i) {
const halfedge_descriptor next_he = next(other_he, P);
Euler::split_face(other_he, he, P);
other_he = next_he;
}
}

18
Convex_hull_3/test/Convex_hull_3/quickhull_test_3.cpp
View File

@ -106,6 +106,22 @@ void test_small_hull()
polyhedron2.size_of_facets() == 6 );
}
void test_coplanar_hull()
{
std::vector<Point_3> points;
points.push_back(Point_3(0,0,0));
points.push_back(Point_3(1,0,0));
points.push_back(Point_3(0,1,0));
points.push_back(Point_3(1,1,0));
points.push_back(Point_3(1.5,0.5,0));
points.push_back(Point_3(0.5,1.1,0));
Polyhedron_3 polyhedron;
CGAL::convex_hull_3(points.begin(), points.end(), polyhedron, Traits());
assert( polyhedron.is_valid() );
assert( polyhedron.size_of_facets() == 4 );
assert( polyhedron.is_pure_triangle() );
}
int main()
{
@ -115,6 +131,8 @@ int main()
test_tetrahedron_convexity();
std::cerr << "Testing small hull" << std::endl;
test_small_hull();
std::cerr << "Testing coplanar hull" << std::endl;
test_coplanar_hull();
std::cerr << "Testing 500 random points" << std::endl;
std::vector<Point_3> points;

4
Mesh_2/doc/Mesh_2/CGAL/lloyd_optimize_mesh_2.h
View File

@ -101,13 +101,13 @@ lloyd_optimize_mesh_2(cdt,
\endcode
\sa `Mesh_optimization_return_code`
\sa `CGAL::Mesh_optimization_return_code`
\sa `CGAL::refine_Delaunay_mesh_2()`
*/
template<typename CDT, typename PointIterator>
Mesh_optimization_return_code
CGAL::Mesh_optimization_return_code
lloyd_optimize_mesh_2(CDT& cdt,
double parameters::time_limit=0,
std::size_t parameters::max_iteration_number=0,

2
Mesh_2/include/CGAL/Mesh_2/Mesh_global_optimizer_2.h
View File

@ -308,7 +308,7 @@ private:
// Find the minimum value
do
{
if(face->is_in_domain())
if (!cdt_.is_infinite(face))
min_sqr = (std::min)(min_sqr, sq_circumradius(face));
face++;
}

4
Point_set_3/include/CGAL/Point_set_3.h
View File

@ -413,7 +413,9 @@ public:
else
{
-- m_nb_removed;
return m_indices.end() - m_nb_removed - 1;
iterator out = m_indices.end() - m_nb_removed - 1;
m_base.reset(*out);
return out;
}
}

12
Point_set_processing_3/include/CGAL/estimate_scale.h
View File

@ -71,7 +71,8 @@ class Quick_multiscale_approximate_knn_distance<Kernel, typename Kernel::Point_3
{
PointPMap point_pmap;
Pmap_unary_function (PointPMap point_pmap) : point_pmap (point_pmap) { }
const typename Kernel::Point_3& operator() (const ValueType& v) const { return get(point_pmap, v); }
typename boost::property_traits<PointPMap>::reference
operator() (const ValueType& v) const { return get(point_pmap, v); }
};
std::size_t m_cluster_size;
@ -227,7 +228,8 @@ class Quick_multiscale_approximate_knn_distance<Kernel, typename Kernel::Point_2
{
PointPMap point_pmap;
Pmap_unary_function (PointPMap point_pmap) : point_pmap (point_pmap) { }
const typename Kernel::Point_2& operator() (const ValueType& v) const { return get(point_pmap, v); }
typename boost::property_traits<PointPMap>::reference
operator() (const ValueType& v) const { return get(point_pmap, v); }
};
template <typename PointPMap>
@ -243,7 +245,8 @@ class Quick_multiscale_approximate_knn_distance<Kernel, typename Kernel::Point_2
: point_pmap (point_pmap) { }
friend inline value_type get (const Pmap_to_3d& ppmap, key_type i)
{
typename Kernel::Point_2 p2 = get(ppmap.point_pmap, i);
typename boost::property_traits<PointPMap>::reference
p2 = get(ppmap.point_pmap, i);
return value_type (p2.x(), p2.y(), 0.);
}
@ -368,7 +371,8 @@ public:
FT nb = 0.;
std::size_t index = 0;
const typename Kernel::Point_2& pquery = get(point_pmap, *query);
typename boost::property_traits<PointPMap>::reference
pquery = get(point_pmap, *query);
for (std::size_t t = 0; t < m_point_sets.size(); ++ t)
{
std::size_t size = ((t == m_point_sets.size() - 1)

6
Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/parameters_interface.h
View File

@ -32,9 +32,9 @@ CGAL_add_named_parameter(protect_constraints_t, protect_constraints, protect_con
CGAL_add_named_parameter(relax_constraints_t, relax_constraints, relax_constraints)
CGAL_add_named_parameter(vertex_is_constrained_t, vertex_is_constrained, vertex_is_constrained_map)
CGAL_add_named_parameter(face_patch_t, face_patch, face_patch_map)
CGAL_add_named_parameter(random_uniform_sampling_t, random_uniform_sampling, random_uniform_sampling)
CGAL_add_named_parameter(grid_sampling_t, grid_sampling, grid_sampling)
CGAL_add_named_parameter(monte_carlo_sampling_t, monte_carlo_sampling, monte_carlo_sampling)
CGAL_add_named_parameter(random_uniform_sampling_t, random_uniform_sampling, use_random_uniform_sampling)
CGAL_add_named_parameter(grid_sampling_t, grid_sampling, use_grid_sampling)
CGAL_add_named_parameter(monte_carlo_sampling_t, monte_carlo_sampling, use_monte_carlo_sampling)
CGAL_add_named_parameter(do_sample_edges_t, do_sample_edges, do_sample_edges)
CGAL_add_named_parameter(do_sample_vertices_t, do_sample_vertices, do_sample_vertices)
CGAL_add_named_parameter(do_sample_faces_t, do_sample_faces, do_sample_faces)

2
Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/measure.h
View File

@ -36,6 +36,8 @@
#include <CGAL/squared_distance_3.h>
#include <CGAL/Kernel/global_functions_3.h>
#include <CGAL/Lazy.h> // needed for CGAL::exact(FT)/CGAL::exact(Lazy_exact_nt<T>)
#include <utility>
#ifdef DOXYGEN_RUNNING

2
Polygon_mesh_processing/test/Polygon_mesh_processing/measures_test.cpp
View File

@ -1,10 +1,10 @@
#include <CGAL/Polygon_mesh_processing/measure.h>
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
#include <CGAL/Polyhedron_3.h>
#include <CGAL/Surface_mesh.h>
#include <CGAL/Polygon_mesh_processing/measure.h>
#include <CGAL/Polygon_mesh_processing/bbox.h>
#include <CGAL/Bbox_3.h>

10
Ridges_3/doc/Ridges_3/CGAL/Ridges.h
View File

@ -170,21 +170,21 @@ Ridge_approximation(const TriangleMesh &tm,
Outputs ridges of types `MAX_ELLIPTIC_RIDGE` and `MAX_HYPERBOLIC_RIDGE`.
\tparam OutputIterator an output iterator wìth value type `Ridge_line*`.
*/
template <class OutputIterator> OutputIterator compute_max_ridges(OutputIterator it, Rdge_order ord = Ridge_order_3);
template <class OutputIterator> OutputIterator compute_max_ridges(OutputIterator it, Ridge_order ord = Ridge_order_3);
/*!
Outputs ridges of types `MIN_ELLIPTIC_RIDGE` and `MIN_HYPERBOLIC_RIDGE`.
\tparam OutputIterator an output iterator with
value type `Ridge_line*`.
*/
template <class OutputIterator> OutputIterator compute_min_ridges(OutputIterator it, Rdge_order ord = Ridge_order_3);
template <class OutputIterator> OutputIterator compute_min_ridges(OutputIterator it, Ridge_order ord = Ridge_order_3);
/*!
Outputs ridges of types `MAX_CREST_RIDGE` and `MIN_CREST_RIDGE`.
\tparam OutputIterator is an output iterator with
value type `Ridge_line*`.
*/
template <class OutputIterator> OutputIterator compute_crest_ridges(OutputIterator it, Rdge_order ord = Ridge_order_3);
template <class OutputIterator> OutputIterator compute_crest_ridges(OutputIterator it, Ridge_order ord = Ridge_order_3);
/// @}
@ -226,7 +226,7 @@ typedef typename TriangleMesh::Traits::FT FT;
A halfedge crossed by a ridge is paired with the barycentric
coordinate of the crossing point.
*/
typedef std::pair< halfedge_descriptor, FT> Ridge_halfhedge;
typedef std::pair< halfedge_descriptor, FT> Ridge_halfedge;
/// @}
@ -261,7 +261,7 @@ FT sharpness() const;
/*!
*/
const std::list<Ridge_halfhedge>* line() const;
const std::list<Ridge_halfedge>* line() const;
/// @}

106
Ridges_3/include/CGAL/Ridges.h
View File

@ -69,7 +69,8 @@ public:
typedef typename Kernel::FT FT;
typedef typename boost::graph_traits<TriangleMesh>::halfedge_descriptor halfedge_descriptor;
typedef std::pair< halfedge_descriptor, FT> ridge_halfhedge;
typedef std::pair< halfedge_descriptor, FT> Ridge_halfedge;
typedef Ridge_halfedge ridge_halfhedge; //kept for backward compatibility
Ridge_type line_type() const {return m_line_type;}
Ridge_type& line_type() {return m_line_type;}
@ -80,8 +81,8 @@ public:
const FT sharpness() const {return m_sharpness;}
FT& sharpness() {return m_sharpness;}
const std::list<ridge_halfhedge>* line() const { return &m_line;}
std::list<ridge_halfhedge>* line() { return &m_line;}
const std::list<Ridge_halfedge>* line() const { return &m_line;}
std::list<Ridge_halfedge>* line() { return &m_line;}
//constructor
Ridge_line(const TriangleMesh& P);
@ -100,7 +101,7 @@ protected:
//one of MAX_ELLIPTIC_RIDGE, MAX_HYPERBOLIC_RIDGE, MAX_CREST_RIDGE,
//MIN_ELLIPTIC_RIDGE, MIN_HYPERBOLIC_RIDGE or MIN_CREST_RIDGE
Ridge_type m_line_type;
std::list<ridge_halfhedge> m_line;
std::list<Ridge_halfedge> m_line;
FT m_strength;// = integral of ppal curvature along the line
FT m_sharpness;// = (integral of second derivative of curvature
// along the line) multiplied by the squared of
@ -131,7 +132,7 @@ dump_4ogl(std::ostream& out_stream,
<< strength() << " "
<< sharpness() << " ";
typedef typename boost::property_traits<VertexPointMap>::value_type Point_3;
typename std::list<ridge_halfhedge >::const_iterator
typename std::list<Ridge_halfedge >::const_iterator
iter = line()->begin(),
ite = line()->end();
for (;iter!=ite;iter++){
@ -156,7 +157,7 @@ dump_verbose(std::ostream& out_stream, VertexPointMap vpm) const
<< "Sharpness is : " << sharpness() << std::endl
<< "Polyline point coordinates are : " << std::endl;
typename std::list<ridge_halfhedge>::const_iterator
typename std::list<Ridge_halfedge>::const_iterator
iter = line()->begin(),
ite = line()->end();
for (;iter!=ite;iter++){
@ -204,7 +205,8 @@ class Ridge_approximation
CGAL_static_assertion((boost::is_same<FT, typename VertexFTMap::value_type>::value));
CGAL_static_assertion((boost::is_same<Vector_3, typename VertexVectorMap::value_type>::value));
typedef std::pair< halfedge_descriptor, FT> Ridge_halfhedge;
typedef std::pair< halfedge_descriptor, FT> Ridge_halfedge;
typedef Ridge_halfedge Ridge_halfhedge; // kept for backward compatibility
typedef CGAL::Ridge_line<TriangleMesh> Ridge_line;
Ridge_approximation(const TriangleMesh &P,
@ -502,7 +504,7 @@ facet_ridge_type(const face_descriptor f, halfedge_descriptor& he1, halfedge_des
//check for regular facet
//i.e. if there is a coherent orientation of ppal dir at the facet vertices
if ( d1[v1]*d1[v2] * d1[v1]*d1[v3] * d1[v2]*d1[v3] < 0 )
if ( get(d1,v1)*get(d1,v2) * get(d1,v1)*get(d1,v3) * get(d1,v2)*get(d1,v3) < 0 )
return NO_RIDGE;
//determine potential crest color
@ -511,11 +513,11 @@ facet_ridge_type(const face_descriptor f, halfedge_descriptor& he1, halfedge_des
Ridge_type crest_color = NO_RIDGE;
if (r_type == CREST_RIDGE)
{
if ( CGAL::abs(k1[v1]+k1[v2]+k1[v3]) > CGAL::abs(k2[v1]+k2[v2]+k2[v3]) )
if ( CGAL::abs(get(k1,v1)+get(k1,v2)+get(k1,v3)) > CGAL::abs(get(k2, v1)+get(k2,v2)+get(k2,v3)) )
crest_color = MAX_CREST_RIDGE;
if ( CGAL::abs(k1[v1]+k1[v2]+k1[v3]) < CGAL::abs(k2[v1]+k2[v2]+k2[v3]) )
if ( CGAL::abs(get(k1,v1)+get(k1,v2)+get(k1,v3)) < CGAL::abs(get(k2,v1)+get(k2,v2)+get(k2,v3)) )
crest_color = MIN_CREST_RIDGE;
if ( CGAL::abs(k1[v1]+k1[v2]+k1[v3]) == CGAL::abs(k2[v1]+k2[v2]+k2[v3]) )
if ( CGAL::abs(get(k1,v1)+get(k1,v2)+get(k1,v3)) == CGAL::abs(get(k2,v1)+get(k2,v2)+get(k2,v3)) )
return NO_RIDGE;
}
@ -587,15 +589,15 @@ xing_on_edge(const halfedge_descriptor he, bool& is_crossed, Ridge_interrogation
is_crossed = false;
FT sign = 0;
FT b_p, b_q; // extremalities at p and q for he: p->q
Vector_3 d_p = d1[target(opposite(he,P),P)],
d_q = d1[target(he,P)]; //ppal dir
Vector_3 d_p = get(d1,target(opposite(he,P),P)),
d_q = get(d1,target(he,P)); //ppal dir
if ( color == MAX_RIDGE ) {
b_p = b0[target(opposite(he,P),P)];
b_q = b0[target(he,P)];
b_p = get(b0,target(opposite(he,P),P));
b_q = get(b0,target(he,P));
}
else {
b_p = b3[target(opposite(he,P),P)];
b_q = b3[target(he,P)];
b_p = get(b3,target(opposite(he,P),P));
b_q = get(b3,target(he,P));
}
if ( b_p == 0 && b_q == 0 ) return;
if ( b_p == 0 && b_q !=0 ) sign = d_p*d_q * b_q;
@ -620,13 +622,13 @@ tag_as_elliptic_hyperbolic(const Ridge_interrogation_type color,
FT coord1, coord2;
if (color == MAX_RIDGE)
{
coord1 = CGAL::abs(b0[v_q1]) / ( CGAL::abs(b0[v_p1]) + CGAL::abs(b0[v_q1]) );
coord2 = CGAL::abs(b0[v_q2]) / ( CGAL::abs(b0[v_p2]) + CGAL::abs(b0[v_q2]) );
coord1 = CGAL::abs(get(b0,v_q1)) / ( CGAL::abs(get(b0,v_p1)) + CGAL::abs(get(b0,v_q1)) );
coord2 = CGAL::abs(get(b0,v_q2)) / ( CGAL::abs(get(b0,v_p2)) + CGAL::abs(get(b0,v_q2)) );
}
else
{
coord1 = CGAL::abs(b3[v_q1]) / ( CGAL::abs(b3[v_p1]) + CGAL::abs(b3[v_q1]) );
coord2 = CGAL::abs(b3[v_q2]) / ( CGAL::abs(b3[v_p2]) + CGAL::abs(b3[v_q2]) );
coord1 = CGAL::abs(get(b3,v_q1)) / ( CGAL::abs(get(b3,v_p1)) + CGAL::abs(get(b3,v_q1)) );
coord2 = CGAL::abs(get(b3,v_q2)) / ( CGAL::abs(get(b3,v_p2)) + CGAL::abs(get(b3,v_q2)) );
}
if ( tag_order == Ridge_order_3 ) {
@ -649,10 +651,10 @@ tag_as_elliptic_hyperbolic(const Ridge_interrogation_type color,
// of Pi at the two crossing points
FT sign_P;
if (color == MAX_RIDGE)
sign_P = P1[v_p1]*coord1 + P1[v_q1]*(1-coord1)
+ P1[v_p2]*coord2 + P1[v_q2]*(1-coord2);
else sign_P = P2[v_p1]*coord1 + P2[v_q1]*(1-coord1)
+ P2[v_p2]*coord2 + P2[v_q2]*(1-coord2);
sign_P = get(P1,v_p1)*coord1 + get(P1,v_q1)*(1-coord1)
+ get(P1,v_p2)*coord2 + get(P1,v_q2)*(1-coord2);
else sign_P = get(P2,v_p1)*coord1 + get(P2,v_q1)*(1-coord1)
+ get(P2,v_p2)*coord2 + get(P2,v_q2)*(1-coord2);
if ( sign_P < 0 ) return true; else return false;
}
@ -672,20 +674,20 @@ b_sign_pointing_to_ridge(const vertex_descriptor v1,
Vector_3 r = r2 - r1, dv1, dv2, dv3;
FT bv1, bv2, bv3;
if ( color == MAX_RIDGE ) {
bv1 = b0[v1];
bv2 = b0[v2];
bv3 = b0[v3];
dv1 = d1[v1];
dv2 = d1[v2];
dv3 = d1[v3];
bv1 = get(b0,v1);
bv2 = get(b0,v2);
bv3 = get(b0,v3);
dv1 = get(d1,v1);
dv2 = get(d1,v2);
dv3 = get(d1,v3);
}
else {
bv1 = b3[v1];
bv2 = b3[v2];
bv3 = b3[v3];
dv1 = d2[v1];
dv2 = d2[v2];
dv3 = d2[v3];
bv1 = get(b3,v1);
bv2 = get(b3,v2);
bv3 = get(b3,v3);
dv1 = get(d2,v1);
dv2 = get(d2,v2);
dv3 = get(d2,v3);
}
if ( r != CGAL::NULL_VECTOR ) r = r/CGAL::sqrt(r*r);
FT sign1, sign2, sign3;
@ -712,7 +714,7 @@ init_ridge_line(Ridge_line* ridge_line,
const Ridge_type r_type)
{
ridge_line->line_type() = r_type;
ridge_line->line()->push_back(Ridge_halfhedge(h1, bary_coord(h1,r_type)));
ridge_line->line()->push_back(Ridge_halfedge(h1, bary_coord(h1,r_type)));
addback(ridge_line, h2, r_type);
}
@ -735,8 +737,8 @@ addback(Ridge_line* ridge_line, const halfedge_descriptor he,
FT k1x, k2x; //abs value of the ppal curvatures at the Xing point on he.
FT k_second = 0; // abs value of the second derivative of the curvature
// along the line of curvature
k1x = CGAL::abs(k1[v_p]) * coord + CGAL::abs(k1[v_q]) * (1-coord) ;
k2x = CGAL::abs(k2[v_p]) * coord + CGAL::abs(k2[v_q]) * (1-coord) ;
k1x = CGAL::abs(get(k1,v_p)) * coord + CGAL::abs(get(k1,v_q)) * (1-coord) ;
k2x = CGAL::abs(get(k2,v_p)) * coord + CGAL::abs(get(k2,v_q)) * (1-coord) ;
if ( (ridge_line->line_type() == MAX_ELLIPTIC_RIDGE)
|| (ridge_line->line_type() == MAX_HYPERBOLIC_RIDGE)
@ -744,7 +746,7 @@ addback(Ridge_line* ridge_line, const halfedge_descriptor he,
ridge_line->strength() += k1x * CGAL::sqrt(segment * segment);
if (tag_order == Ridge_order_4) {
if (k1x != k2x)
k_second =CGAL::abs(( CGAL::abs(P1[v_p]) * coord + CGAL::abs(P1[v_q]) * (1-coord) )/(k1x-k2x));
k_second =CGAL::abs(( CGAL::abs(get(P1,v_p)) * coord + CGAL::abs(get(P1,v_q)) * (1-coord) )/(k1x-k2x));
ridge_line->sharpness() += k_second * CGAL::sqrt(segment * segment) * squared_model_size; }
}
if ( (ridge_line->line_type() == MIN_ELLIPTIC_RIDGE)
@ -753,10 +755,10 @@ addback(Ridge_line* ridge_line, const halfedge_descriptor he,
ridge_line->strength() += k2x * CGAL::sqrt(segment * segment);
if (tag_order == Ridge_order_4) {
if (k1x != k2x)
k_second =CGAL::abs(( CGAL::abs(P2[v_p]) * coord + CGAL::abs(P2[v_q]) * (1-coord) )/(k1x-k2x));
k_second =CGAL::abs(( CGAL::abs(get(P2,v_p)) * coord + CGAL::abs(get(P2,v_q)) * (1-coord) )/(k1x-k2x));
ridge_line->sharpness() += k_second * CGAL::sqrt(segment * segment) * squared_model_size; }
}
ridge_line->line()->push_back( Ridge_halfhedge(he, coord));
ridge_line->line()->push_back( Ridge_halfedge(he, coord));
}
@ -779,8 +781,8 @@ addfront(Ridge_line* ridge_line,
FT k1x, k2x; //abs value of the ppal curvatures at the Xing point on he.
FT k_second = 0.; // abs value of the second derivative of the curvature
// along the line of curvature
k1x = CGAL::abs(k1[v_p]) * coord + CGAL::abs(k1[v_q]) * (1-coord) ;
k2x = CGAL::abs(k2[v_p]) * coord + CGAL::abs(k2[v_q]) * (1-coord) ;
k1x = CGAL::abs(get(k1,v_p)) * coord + CGAL::abs(get(k1,v_q)) * (1-coord) ;
k2x = CGAL::abs(get(k2,v_p)) * coord + CGAL::abs(get(k2,v_q)) * (1-coord) ;
if ( (ridge_line->line_type() == MAX_ELLIPTIC_RIDGE)
|| (ridge_line->line_type() == MAX_HYPERBOLIC_RIDGE)
@ -788,7 +790,7 @@ addfront(Ridge_line* ridge_line,
ridge_line->strength() += k1x * CGAL::sqrt(segment * segment);
if (tag_order == Ridge_order_4) {
if (k1x != k2x)
k_second =CGAL::abs(( CGAL::abs(P1[v_p]) * coord + CGAL::abs(P1[v_q]) * (1-coord) )/(k1x-k2x));
k_second =CGAL::abs(( CGAL::abs(get(P1,v_p)) * coord + CGAL::abs(get(P1,v_q)) * (1-coord) )/(k1x-k2x));
ridge_line->sharpness() += k_second * CGAL::sqrt(segment * segment) * squared_model_size; }
}
if ( (ridge_line->line_type() == MIN_ELLIPTIC_RIDGE)
@ -797,10 +799,10 @@ addfront(Ridge_line* ridge_line,
ridge_line->strength() += k2x * CGAL::sqrt(segment * segment);
if (tag_order == Ridge_order_4) {
if (k1x != k2x)
k_second =CGAL::abs(( CGAL::abs(P2[v_p]) * coord + CGAL::abs(P2[v_q]) * (1-coord) )/(k1x-k2x));
k_second =CGAL::abs(( CGAL::abs(get(P2,v_p)) * coord + CGAL::abs(get(P2,v_q)) * (1-coord) )/(k1x-k2x));
ridge_line->sharpness() += k_second * CGAL::sqrt(segment * segment) * squared_model_size; }
}
ridge_line->line()->push_front( Ridge_halfhedge(he, coord));
ridge_line->line()->push_front( Ridge_halfedge(he, coord));
}
@ -815,14 +817,14 @@ bary_coord(const halfedge_descriptor he, const Ridge_type r_type)
if ( (r_type == MAX_ELLIPTIC_RIDGE)
|| (r_type == MAX_HYPERBOLIC_RIDGE)
|| (r_type == MAX_CREST_RIDGE) ) {
b_p = b0[target(opposite(he,P),P)];
b_q = b0[target(he,P)];
b_p = get(b0,target(opposite(he,P),P));
b_q = get(b0,target(he,P));
}
if ( (r_type == MIN_ELLIPTIC_RIDGE)
|| (r_type == MIN_HYPERBOLIC_RIDGE)
|| (r_type == MIN_CREST_RIDGE) ) {
b_p = b3[target(opposite(he,P),P)];
b_q = b3[target(he,P)];
b_p = get(b3,target(opposite(he,P),P));
b_q = get(b3,target(he,P));
}
//denominator cannot be 0 since there is no crossing when both extremalities are 0
return CGAL::abs(b_q) / ( CGAL::abs(b_q) + CGAL::abs(b_p) );

14
Ridges_3/include/CGAL/Umbilics.h
View File

@ -194,7 +194,7 @@ compute(OutputIterator umbilics_it, FT size)
boost::tie(itb,ite) = vertices(P);
for (;itb != ite; itb++) {
vertex_descriptor vh = *itb;
umbilicEstimatorVertex = cgal_abs(k1[vh]-k2[vh]);
umbilicEstimatorVertex = cgal_abs(get(k1,vh)-get(k2,vh));
//reset vector, list and bool
vces.clear();
contour.clear();
@ -218,7 +218,7 @@ compute(OutputIterator umbilics_it, FT size)
itev = vces.end();
itbv++;
for (; itbv != itev; itbv++)
{ umbilicEstimatorNeigh = cgal_abs( k1[*itbv] - k2[*itbv] );
{ umbilicEstimatorNeigh = cgal_abs( get(k1,*itbv) - get(k2,*itbv) );
if ( umbilicEstimatorNeigh < umbilicEstimatorVertex )
{is_umbilic = false; break;}
}
@ -245,14 +245,14 @@ compute_type(Umbilic& umb)
itlast = --umb.contour_list().end();
v = target(*itb, P);
dir = d1[v];
normal = CGAL::cross_product(d1[v], d2[v]);
dir = get(d1,v);
normal = CGAL::cross_product(get(d1,v), get(d2,v));
//sum angles along the contour
do{
itb++;
v = target(*itb, P);
dirnext = d1[v];
dirnext = get(d1,v);
cosinus = To_double(dir*dirnext);
if (cosinus < 0) {dirnext = dirnext*(-1); cosinus *= -1;}
if (cosinus>1) cosinus = 1;
@ -261,13 +261,13 @@ compute_type(Umbilic& umb)
else angle = -acos(cosinus);
angleSum += angle;
dir = dirnext;
normal = CGAL::cross_product(d1[v], d2[v]);
normal = CGAL::cross_product(get(d1,v), get(d2,v));
}
while (itb != (itlast));
//angle (v_last, v_0)
v = target(*umb.contour_list().begin(), P);
dirnext = d1[v];
dirnext = get(d1,v);
cosinus = To_double(dir*dirnext);
if (cosinus < 0) {dirnext = dirnext*(-1); cosinus *= -1;}
if (cosinus>1) cosinus = 1;

77
Subdivision_method_3/doc/Subdivision_method_3/Subdivision_method_3.txt
View File

@ -127,12 +127,12 @@ There is only one line deserving a detailed explanation:
\code{.cpp}
Subdivision_method_3::CatmullClark_subdivision(pmesh, d);
Subdivision_method_3::CatmullClark_subdivision(pmesh, params::number_of_iterations(d));
\endcode
`Subdivision_method_3` specifies the namespace of the
subdivision functions. `CatmullClark_subdivision(P, d)` computes the
subdivision functions. `CatmullClark_subdivision(P, params::number_of_iterations(d))` computes the
Catmull-Clark subdivision surface of the polygon mesh `pmesh` after
`d` iterations of the refinements. The polygon mesh `pmesh` is
passed by reference, and is modified (i.e.\ subdivided) by the
@ -150,12 +150,6 @@ in the internal container are sequentially ordered (e.g.
This implies that the iterators traverse the primitives in
the order of their creations/insertions.
\attention The four subdivision techniques have been made to work with
`Surface_mesh` despite `Surface_mesh` not necessarily inserting new elements at
the end of its containers (see Section \ref sectionSurfaceMesh_properties of the
`Surface_mesh`). However, it is required that the `Surface_mesh` passed as input
is garbage-free, which can be achieved by first calling `Surface_mesh::collect_garbage()`.
Section \ref secRefHost gives detailed explanations on this
two restrictions.
@ -187,8 +181,8 @@ three components of the Catmull-Clark subdivision method.
\code{.cpp}
template <class PolygonMesh, class Mask>
void PQQ(PolygonMesh& p, Mask mask, int depth)
template <class PolygonMesh, class Mask, class NamedParameters>
void PQQ(PolygonMesh& p, Mask mask, NamedParameters np)
\endcode
@ -203,8 +197,7 @@ new points by cooperating with the `mask`.
To implement Catmull-Clark subdivision,
`Mask`, the <I>geometry policy</I>, has to realize the
geometry masks of Catmull-Clark subdivision.
The parameter `depth` specifies the iterations of the refinement
on the control mesh.
The number of iterations as well as the vertex point map can be specified using the named parameter `np`.
To implement the geometry masks, we need to know how
a refinement host communicates with its geometry masks.
@ -258,8 +251,8 @@ class CatmullClark_mask_3 {
};
void edge_node(halfedge_descriptor edge, Point_3& pt) {
Point_3& p1 = get(vpm,target(edge, pmesh));
Point_3& p2 = get(vpm,source(edge, pmesh));
Point_3 p1 = get(vpm,target(edge, pmesh));
Point_3 p2 = get(vpm,source(edge, pmesh));
Point_3 f1, f2;
face_node(face(edge,pmesh), f1);
face_node(face(opposite(edge,pmesh),pmesh), f2);
@ -273,11 +266,11 @@ class CatmullClark_mask_3 {
typename boost::graph_traits<PolygonMesh>::degree_size_type n = degree(vertex,pmesh);
FT Q[] = {0.0, 0.0, 0.0}, R[] = {0.0, 0.0, 0.0};
Point_3& S = get(vpm,vertex);
Point_3 S = get(vpm,vertex);
Point_3 q;
for (int i = 0; i < n; i++, ++vcir) {
Point_3& p2 = get(vpm, target(opposite(*vcir,pmesh),pmesh));
Point_3 p2 = get(vpm, target(opposite(*vcir,pmesh),pmesh));
R[0] += (S[0]+p2[0])/2;
R[1] += (S[1]+p2[1])/2;
R[2] += (S[2]+p2[2])/2;
@ -305,7 +298,7 @@ call `PQQ()` with the Catmull-Clark masks that we have just defined.
\code{.cpp}
PQQ(pmesh, CatmullClark_mask_3(pmesh), depth);
PQQ(pmesh, CatmullClark_mask_3(pmesh), params::number_of_iterations(depth));
\endcode
@ -334,17 +327,17 @@ and \f$ \sqrt{3}\f$ subdivision.
\code{.cpp}
namespace Subdivision_method_3 {
template <class PolygonMesh, class Mask>
void PQQ(PolygonMesh& pmesh, Mask mask, int step);
template <class PolygonMesh, class Mask, class NamedParameters>
void PQQ(PolygonMesh& pmesh, Mask mask, NamedParameters np);
template <class PolygonMesh, class Mask>
void PTQ(PolygonMesh& pmesh, Mask mask, int step);
template <class PolygonMesh, class Mask, class NamedParameters>
void PTQ(PolygonMesh& pmesh, Mask mask, NamedParameters np);
template <class PolygonMesh, class Mask>
void DQQ(PolygonMesh& pmesh, Mask mask, int step)
template <class PolygonMesh, class Mask, class NamedParameters>
void DQQ(PolygonMesh& pmesh, Mask mask, NamedParameters np)
template <class PolygonMesh, class Mask>
void Sqrt3(PolygonMesh& pmesh, Mask mask, int step)
template <class PolygonMesh, class Mask, class NamedParameters>
void Sqrt3(PolygonMesh& pmesh, Mask mask, NamedParameters np)
}
\endcode
@ -461,15 +454,15 @@ is listed below, where `ept` returns the new point splitting
\code{.cpp}
void border_node(halfedge_descriptor edge, Point_3& ept, Point_3& vpt) {
Point_3& ep1 = get(vpm, target(edge, pmesh));
Point_3& ep2 = get(vpm, target(opposite(edge, pmesh), pmesh));
Point_3 ep1 = get(vpm, target(edge, pmesh));
Point_3 ep2 = get(vpm, target(opposite(edge, pmesh), pmesh));
ept = Point_3((ep1[0]+ep2[0])/2, (ep1[1]+ep2[1])/2, (ep1[2]+ep2[2])/2);
Halfedge_around_target_circulator<Poly> vcir(edge, pmesh);
Point_3& vp1 = get(vpm,target(opposite(*vcir, pmesh ), pmesh));
Point_3& vp0 = get(vpm, target(*vcir, pmesh));
Point_3 vp1 = get(vpm,target(opposite(*vcir, pmesh ), pmesh));
Point_3 vp0 = get(vpm, target(*vcir, pmesh));
--vcir;
Point_3& vp_1 = get(vpm, target(opposite(*vcir, pmesh), pmesh));
Point_3 vp_1 = get(vpm, target(opposite(*vcir, pmesh), pmesh));
vpt = Point_3((vp_1[0] + 6*vp0[0] + vp1[0])/8,
(vp_1[1] + 6*vp0[1] + vp1[1])/8,
(vp_1[2] + 6*vp0[2] + vp1[2])/8 );
@ -533,24 +526,24 @@ based on that kernel.
\code{.cpp}
namespace Subdivision_method_3 {
template <class PolygonMesh>
void CatmullClark_subdivision(PolygonMesh& pmesh, int step = 1) {
PQQ(pmesh, CatmullClark_mask_3<PolygonMesh>(pmesh), step);
template <class PolygonMesh, class NamedParameters>
void CatmullClark_subdivision(PolygonMesh& pmesh, NamedParameters np) {
PQQ(pmesh, CatmullClark_mask_3<PolygonMesh>(pmesh), np);
}
template <class PolygonMesh>
void Loop_subdivision(PolygonMesh& pmesh, int step = 1) {
PTQ(pmesh, Loop_mask_3<PolygonMesh>(pmesh) , step);
template <class PolygonMesh, class NamedParameters>
void Loop_subdivision(PolygonMesh& pmesh, NamedParameters np) {
PTQ(pmesh, Loop_mask_3<PolygonMesh>(pmesh) , np);
}
template <class PolygonMesh>
void DooSabin_subdivision(PolygonMesh& pmesh, int step = 1) {
DQQ(pmesh, DooSabin_mask_3<PolygonMesh>(pmesh), step);
template <class PolygonMesh, class NamedParameters>
void DooSabin_subdivision(PolygonMesh& pmesh, NamedParameters np) {
DQQ(pmesh, DooSabin_mask_3<PolygonMesh>(pmesh), np);
}
template <class PolygonMesh>
void Sqrt3_subdivision(PolygonMesh& pmesh, int step = 1) {
Sqrt3(pmesh, Sqrt3_mask_3<PolygonMesh>(pmesh), step);
template <class PolygonMesh, class NamedParameters>
void Sqrt3_subdivision(PolygonMesh& pmesh, NamedParameters np) {
Sqrt3(pmesh, Sqrt3_mask_3<PolygonMesh>(pmesh), np);
}
}

3
Subdivision_method_3/examples/Subdivision_method_3/CatmullClark_subdivision.cpp
View File

@ -16,6 +16,7 @@ typedef CGAL::Surface_mesh<Kernel::Point_3> PolygonMesh;
using namespace std;
using namespace CGAL;
namespace params = CGAL::Polygon_mesh_processing::parameters;
int main(int argc, char** argv) {
if (argc > 4) {
@ -40,7 +41,7 @@ int main(int argc, char** argv) {
Timer t;
t.start();
Subdivision_method_3::CatmullClark_subdivision(pmesh, d);
Subdivision_method_3::CatmullClark_subdivision(pmesh, params::number_of_iterations(d));
std::cerr << "Done (" << t.time() << " s)" << std::endl;
std::ofstream out(out_file);

26
Subdivision_method_3/examples/Subdivision_method_3/Customized_subdivision.cpp
View File

@ -16,6 +16,7 @@ typedef CGAL::Simple_cartesian<double> Kernel;
typedef CGAL::Surface_mesh<Kernel::Point_3> PolygonMesh;
using namespace std;
using namespace CGAL;
namespace params = CGAL::Polygon_mesh_processing::parameters;
// ======================================================================
template <class Poly>
@ -27,6 +28,7 @@ class WLoop_mask_3 {
typedef typename boost::property_map<PolygonMesh, vertex_point_t>::type Vertex_pmap;
typedef typename boost::property_traits<Vertex_pmap>::value_type Point;
typedef typename boost::property_traits<Vertex_pmap>::reference Point_ref;
PolygonMesh& pmesh;
Vertex_pmap vpm;
@ -37,10 +39,10 @@ public:
{}
void edge_node(halfedge_descriptor hd, Point& pt) {
Point& p1 = get(vpm, target(hd,pmesh));
Point& p2 = get(vpm, target(opposite(hd,pmesh),pmesh));
Point& f1 = get(vpm, target(next(hd,pmesh),pmesh));
Point& f2 = get(vpm, target(next(opposite(hd,pmesh),pmesh),pmesh));
Point_ref p1 = get(vpm, target(hd,pmesh));
Point_ref p2 = get(vpm, target(opposite(hd,pmesh),pmesh));
Point_ref f1 = get(vpm, target(next(hd,pmesh),pmesh));
Point_ref f2 = get(vpm, target(next(opposite(hd,pmesh),pmesh),pmesh));
pt = Point((3*(p1[0]+p2[0])+f1[0]+f2[0])/8,
(3*(p1[1]+p2[1])+f1[1]+f2[1])/8,
@ -48,12 +50,12 @@ public:
}
void vertex_node(vertex_descriptor vd, Point& pt) {
double R[] = {0.0, 0.0, 0.0};
Point& S = get(vpm,vd);
Point_ref S = get(vpm,vd);
std::size_t n = 0;
BOOST_FOREACH(halfedge_descriptor hd, halfedges_around_target(vd, pmesh)){
++n;
Point& p = get(vpm, target(opposite(hd,pmesh),pmesh));
Point_ref p = get(vpm, target(opposite(hd,pmesh),pmesh));
R[0] += p[0]; R[1] += p[1]; R[2] += p[2];
}
@ -71,15 +73,15 @@ public:
}
void border_node(halfedge_descriptor hd, Point& ept, Point& vpt) {
Point& ep1 = get(vpm, target(hd,pmesh));
Point& ep2 = get(vpm, target(opposite(hd,pmesh),pmesh));
Point_ref ep1 = get(vpm, target(hd,pmesh));
Point_ref ep2 = get(vpm, target(opposite(hd,pmesh),pmesh));
ept = Point((ep1[0]+ep2[0])/2, (ep1[1]+ep2[1])/2, (ep1[2]+ep2[2])/2);
Halfedge_around_target_circulator<Poly> vcir(hd,pmesh);
Point& vp1 = get(vpm, target(opposite(*vcir,pmesh),pmesh));
Point& vp0 = get(vpm, target(*vcir,pmesh));
Point_ref vp1 = get(vpm, target(opposite(*vcir,pmesh),pmesh));
Point_ref vp0 = get(vpm, target(*vcir,pmesh));
--vcir;
Point& vp_1 = get(vpm,target(opposite(*vcir,pmesh),pmesh));
Point_ref vp_1 = get(vpm,target(opposite(*vcir,pmesh),pmesh));
vpt = Point((vp_1[0] + 6*vp0[0] + vp1[0])/8,
(vp_1[1] + 6*vp0[1] + vp1[1])/8,
(vp_1[2] + 6*vp0[2] + vp1[2])/8 );
@ -109,7 +111,7 @@ int main(int argc, char **argv) {
Timer t;
t.start();
Subdivision_method_3::PTQ(pmesh, WLoop_mask_3<PolygonMesh>(pmesh), d);
Subdivision_method_3::PTQ(pmesh, WLoop_mask_3<PolygonMesh>(pmesh), params::number_of_iterations(d));
std::cerr << "Done (" << t.time() << " s)" << std::endl;
std::ofstream out(out_file);

3
Subdivision_method_3/examples/Subdivision_method_3/DooSabin_subdivision.cpp
View File

@ -18,6 +18,7 @@ typedef CGAL::Polyhedron_3<Kernel> PolygonMesh;
using namespace std;
using namespace CGAL;
namespace params = CGAL::Polygon_mesh_processing::parameters;
int main(int argc, char **argv) {
if (argc > 4) {
@ -42,7 +43,7 @@ int main(int argc, char **argv) {
Timer t;
t.start();
Subdivision_method_3::DooSabin_subdivision(pmesh, d);
Subdivision_method_3::DooSabin_subdivision(pmesh, params::number_of_iterations(d));
std::cerr << "Done (" << t.time() << " s)" << std::endl;
std::ofstream out(out_file);

3
Subdivision_method_3/examples/Subdivision_method_3/Loop_subdivision.cpp
View File

@ -16,6 +16,7 @@ typedef CGAL::Surface_mesh<Kernel::Point_3> PolygonMesh;
using namespace std;
using namespace CGAL;
namespace params = CGAL::Polygon_mesh_processing::parameters;
int main(int argc, char **argv) {
if (argc > 4) {
@ -40,7 +41,7 @@ int main(int argc, char **argv) {
Timer t;
t.start();
Subdivision_method_3::Loop_subdivision(pmesh,d);
Subdivision_method_3::Loop_subdivision(pmesh, params::number_of_iterations(d));
std::cerr << "Done (" << t.time() << " s)" << std::endl;
std::ofstream out(out_file);

3
Subdivision_method_3/examples/Subdivision_method_3/Sqrt3_subdivision.cpp
View File

@ -16,6 +16,7 @@ typedef CGAL::Surface_mesh<Kernel::Point_3> PolygonMesh;
using namespace std;
using namespace CGAL;
namespace params = CGAL::Polygon_mesh_processing::parameters;
int main(int argc, char **argv) {
if (argc > 4) {
@ -40,7 +41,7 @@ int main(int argc, char **argv) {
Timer t;
t.start();
Subdivision_method_3::Sqrt3_subdivision(pmesh,d);
Subdivision_method_3::Sqrt3_subdivision(pmesh, params::number_of_iterations(d));
std::cerr << "Done (" << t.time() << " s)" << std::endl;
std::ofstream out(out_file);

355
Subdivision_method_3/include/CGAL/Subdivision_method_3/internal/subdivision_hosts_impl_3.h
View File

@ -33,11 +33,8 @@
#include <CGAL/circulator.h>
#include <CGAL/tags.h>
#include <boost/array.hpp>
#include <boost/foreach.hpp>
#include <boost/mpl/if.hpp>
#include <boost/unordered_map.hpp>
#include <boost/unordered_set.hpp>
#include <iterator>
#include <list>
@ -49,18 +46,34 @@ namespace Subdivision_method_3 {
namespace internal {
template <class PolygonMesh>
void call_reserve(PolygonMesh& pm, std::size_t v, std::size_t e, std::size_t f)
{
typedef boost::graph_traits<PolygonMesh> GT;
reserve(pm,
static_cast<typename GT::vertices_size_type>(v),
static_cast<typename GT::edges_size_type>(e),
static_cast<typename GT::faces_size_type>(f));
}
template <class Poly, class VertexPointMap, class Mask>
void PQQ_1step(Poly& p, VertexPointMap vpm, Mask mask) {
typedef typename boost::graph_traits<Poly>::vertex_descriptor vertex_descriptor;
typedef typename boost::graph_traits<Poly>::halfedge_descriptor halfedge_descriptor;
typedef typename boost::graph_traits<Poly>::edge_descriptor edge_descriptor;
typedef typename boost::graph_traits<Poly>::vertex_iterator vertex_iterator;
typedef typename boost::graph_traits<Poly>::edge_iterator edge_iterator;
typedef typename boost::graph_traits<Poly>::face_iterator face_iterator;
typedef typename boost::graph_traits<Poly>::face_descriptor face_descriptor;
typedef Halfedge_around_face_circulator<Poly> Halfedge_around_facet_circulator;
//First back up initial vertices/faces/edges
std::vector<vertex_descriptor> p_vertices(vertices(p).first, vertices(p).second);
std::vector<face_descriptor> p_faces(faces(p).first, faces(p).second);
std::vector<edge_descriptor> p_edges(edges(p).first, edges(p).second);
std::size_t num_v = p_vertices.size();
std::size_t num_e = p_edges.size();
std::size_t num_f = p_faces.size();
// Build a new vertices buffer has the following structure
//
// 0 1 ... e_begin ... f_begin ... (end_of_buffer)
@ -69,34 +82,29 @@ void PQQ_1step(Poly& p, VertexPointMap vpm, Mask mask) {
// f_begin ... (end) : store the positions of the face-vertices
// The index of the vertices buffer should 1-1 map to the distance
// of the corresponding iterator to the begin of the iterator.
typename boost::graph_traits<Poly>::vertices_size_type num_vertex = num_vertices(p);
typename boost::graph_traits<Poly>::halfedges_size_type num_edge = num_halfedges(p)/2;
typename boost::graph_traits<Poly>::faces_size_type num_facet = num_faces(p);
// We need to reserve the memory to prevent reallocation.
reserve(p,num_vertex+num_edge+num_facet, 4*2*num_edge, 4*num_edge/2);
call_reserve(p,num_v+num_e+num_f, 4*2*num_e, 4*num_e/2);
typedef typename boost::property_traits<VertexPointMap>::value_type Point;
Point* vertex_point_buffer = new Point[num_vertex + num_edge + num_facet];
Point* edge_point_buffer = vertex_point_buffer + num_vertex;
Point* face_point_buffer = edge_point_buffer + num_edge;
Point* vertex_point_buffer = new Point[num_v + num_e + num_f];
Point* edge_point_buffer = vertex_point_buffer + num_v;
Point* face_point_buffer = edge_point_buffer + num_e;
int i=0;
boost::unordered_map<vertex_descriptor,int> v_index;
BOOST_FOREACH(vertex_descriptor vh, vertices(p)){
BOOST_FOREACH(vertex_descriptor vh, p_vertices){
v_index[vh]= i++;
}
std::vector<bool> v_onborder(num_vertex);
face_iterator fitr = faces(p).first;
for (size_t i = 0; i < num_facet; i++, ++fitr)
std::vector<bool> v_onborder(num_v);
typename std::vector<face_descriptor>::iterator fitr=p_faces.begin();
for (std::size_t i = 0; i < num_f; i++, ++fitr)
mask.face_node(*fitr, face_point_buffer[i]);
{
std::size_t i = 0;
BOOST_FOREACH(edge_descriptor ed, edges(p)){
BOOST_FOREACH(edge_descriptor ed, p_edges){
if(is_border(ed,p)){
halfedge_descriptor h=halfedge(ed,p);
if (is_border(h,p)) h=opposite(h,p);
@ -111,8 +119,8 @@ void PQQ_1step(Poly& p, VertexPointMap vpm, Mask mask) {
}
}
vertex_iterator vitr = vertices(p).first;
for (size_t i = 0; i < num_vertex; i++, ++vitr)
typename std::vector<vertex_descriptor>::iterator vitr = p_vertices.begin();
for (std::size_t i = 0; i < num_v; i++, ++vitr)
if (!v_onborder[v_index[*vitr]]) mask.vertex_node(*vitr, vertex_point_buffer[i]);
// Build the connectivity using insert_vertex() and insert_edge()
@ -123,16 +131,16 @@ void PQQ_1step(Poly& p, VertexPointMap vpm, Mask mask) {
// 4. insert_edge() between all other new inserted vertices of step 1 and
// the new inserted vertex of step 3
// Step 1.
edge_iterator eitr = edges(p).first;
for (size_t i = 0; i < num_edge; i++, ++eitr) {
typename std::vector<edge_descriptor>::iterator eitr = p_edges.begin();
for (std::size_t i = 0; i < num_e; i++, ++eitr) {
vertex_descriptor vh = insert_vertex(p, halfedge(*eitr,p));
put(vpm, vh, edge_point_buffer[i]);
}
fitr = faces(p).first;
// TODO: the topoloy modification can be done by a template function
// and that gives the user a chance to create new topological masks.
for (size_t i = 0; i < num_facet; i++, ++fitr) {
fitr = p_faces.begin();
for (std::size_t i = 0; i < num_f; i++, ++fitr) {
// Step 2.
Halfedge_around_facet_circulator hcir_begin(halfedge(*fitr,p),p);
Halfedge_around_facet_circulator hcir = hcir_begin;
@ -161,8 +169,8 @@ void PQQ_1step(Poly& p, VertexPointMap vpm, Mask mask) {
// Update the geometry data of the newly inserted vertices by the
// vertices buffer
vitr = vertices(p).first;
for (size_t i = 0; i < num_vertex; i++, ++vitr)
vitr = p_vertices.begin();
for (std::size_t i = 0; i < num_v; i++, ++vitr)
put(vpm, *vitr, vertex_point_buffer[i]);
// CGAL_postcondition(p.is_valid());
@ -175,15 +183,21 @@ void PTQ_1step(Poly& p, VertexPointMap vpm, Mask mask) {
typedef typename boost::graph_traits<Poly>::vertex_descriptor vertex_descriptor;
typedef typename boost::graph_traits<Poly>::halfedge_descriptor halfedge_descriptor;
typedef typename boost::graph_traits<Poly>::edge_descriptor edge_descriptor;
typedef typename boost::graph_traits<Poly>::vertex_iterator vertex_iterator;
typedef typename boost::graph_traits<Poly>::edge_iterator edge_iterator;
typedef typename boost::graph_traits<Poly>::face_iterator face_iterator;
typedef typename boost::graph_traits<Poly>::face_descriptor face_descriptor;
typedef Halfedge_around_face_circulator<Poly> Halfedge_around_face_circulator;
typedef typename boost::property_traits<VertexPointMap>::value_type Point;
//First back up initial vertices/faces/edges
std::vector<vertex_descriptor> p_vertices(vertices(p).first, vertices(p).second);
std::vector<face_descriptor> p_faces(faces(p).first, faces(p).second);
std::vector<edge_descriptor> p_edges(edges(p).first, edges(p).second);
std::size_t num_v = p_vertices.size();
std::size_t num_e = p_edges.size();
std::size_t num_f = p_faces.size();
// Build a new vertices buffer has the following structure
//
// 0 1 ... e_begin ... f_begin ... (end_of_buffer)
@ -191,26 +205,23 @@ void PTQ_1step(Poly& p, VertexPointMap vpm, Mask mask) {
// e_begin ... (end) : store the positions of the edge-vertices
// The index of the vertices buffer should 1-1 map to the distance
// of the corresponding iterator to the begin of the iterator.
typename boost::graph_traits<Poly>::vertices_size_type num_vertex = num_vertices(p);
typename boost::graph_traits<Poly>::halfedges_size_type num_edge = num_halfedges(p)/2;
typename boost::graph_traits<Poly>::faces_size_type num_facet = num_faces(p);
// We need to reserve the memory to prevent reallocation.
reserve(p,num_vertex+num_edge, 2*2*num_edge, 4*num_edge/2);
call_reserve(p,num_v + num_e, 2*2*num_e, 4*num_e/2);
Point* vertex_point_buffer = new Point[num_vertex + num_edge];
Point* edge_point_buffer = vertex_point_buffer + num_vertex;
Point* vertex_point_buffer = new Point[num_v + num_e];
Point* edge_point_buffer = vertex_point_buffer + num_v;
int i=0;
boost::unordered_map<vertex_descriptor,int> v_index;
BOOST_FOREACH(vertex_descriptor vh, vertices(p)){
BOOST_FOREACH(vertex_descriptor vh, p_vertices){
v_index[vh]= i++;
}
std::vector<bool> v_onborder(num_vertex);
std::vector<bool> v_onborder(num_v);
{
std::size_t i = 0;
BOOST_FOREACH(edge_descriptor ed, edges(p)){
BOOST_FOREACH(edge_descriptor ed, p_edges){
if(! is_border(ed,p)){
mask.edge_node(halfedge(ed,p), edge_point_buffer[i]);
} else{
@ -223,8 +234,8 @@ void PTQ_1step(Poly& p, VertexPointMap vpm, Mask mask) {
++i;
}
}
vertex_iterator vitr = vertices(p).first;
for (size_t i = 0; i < num_vertex; i++, ++vitr)
typename std::vector<vertex_descriptor>::iterator vitr = p_vertices.begin();
for (std::size_t i = 0; i < num_v; i++, ++vitr)
if (!v_onborder[i]) mask.vertex_node(*vitr, vertex_point_buffer[i]);
// Build the connectivity using insert_vertex() and insert_edge()
@ -235,13 +246,14 @@ void PTQ_1step(Poly& p, VertexPointMap vpm, Mask mask) {
// 4. insert_edge() between all other new inserted vertices of step 1 and
// the new inserted vertex of step 3
// Step 1.
edge_iterator eitr = edges(p).first;
for (size_t i = 0; i < num_edge; i++, ++eitr) {
vertex_descriptor vh = insert_vertex(p, halfedge(*eitr,p));
typename std::vector<edge_descriptor>::iterator eit = p_edges.begin();
for (std::size_t i = 0; i < num_e; i++, ++eit) {
vertex_descriptor vh = insert_vertex(p, halfedge(*eit,p));
put(vpm,vh, edge_point_buffer[i]);
}
face_iterator fitr = faces(p).first;
for (size_t i = 0; i < num_facet; i++, ++fitr) {
typename std::vector<face_descriptor>::iterator fitr = p_faces.begin();
for (std::size_t i = 0; i < num_f; i++, ++fitr) {
// Step 2.
Halfedge_around_face_circulator hcir_begin(halfedge(*fitr,p),p);
Halfedge_around_face_circulator hcir = hcir_begin;
@ -261,8 +273,8 @@ void PTQ_1step(Poly& p, VertexPointMap vpm, Mask mask) {
// Update the geometry data of the newly inserted vertices by the
// vertices buffer
vitr = vertices(p).first;
for (size_t i = 0; i < num_vertex; i++, ++vitr)
vitr = p_vertices.begin();
for (std::size_t i = 0; i < num_v; i++, ++vitr)
put(vpm, *vitr, vertex_point_buffer[i]);
// CGAL_postcondition(p.is_valid());
@ -272,54 +284,60 @@ void PTQ_1step(Poly& p, VertexPointMap vpm, Mask mask) {
// ======================================================================
template <class Poly, class VertexPointMap, class Mask>
void DQQ_1step_impl(Poly& p, VertexPointMap vpm, Mask mask, CGAL::Tag_false) {
void DQQ_1step(Poly& p, VertexPointMap vpm, Mask mask) {
typedef typename boost::graph_traits<Poly>::vertex_descriptor vertex_descriptor;
typedef typename boost::graph_traits<Poly>::halfedge_descriptor halfedge_descriptor;
typedef typename boost::graph_traits<Poly>::edge_descriptor edge_descriptor;
typedef typename boost::graph_traits<Poly>::face_descriptor face_descriptor;
typedef typename boost::graph_traits<Poly>::vertex_iterator vertex_iterator;
typedef typename boost::graph_traits<Poly>::edge_iterator edge_iterator;
typedef typename boost::property_traits<VertexPointMap>::value_type Point;
typename boost::graph_traits<Poly>::vertices_size_type num_v = num_vertices(p);
typename boost::graph_traits<Poly>::halfedges_size_type num_e = num_halfedges(p)/2;
typename boost::graph_traits<Poly>::faces_size_type num_f = num_faces(p);
//First back up initial vertices/faces/edges
std::vector<vertex_descriptor> p_vertices(vertices(p).first, vertices(p).second);
std::vector<face_descriptor> p_faces(faces(p).first, faces(p).second);
std::vector<edge_descriptor> p_edges(edges(p).first, edges(p).second);
std::size_t num_v = p_vertices.size();
std::size_t num_e = p_edges.size();
std::size_t num_f = p_faces.size();
std::vector<halfedge_descriptor> border_halfedges;
size_t num_be = 0 ;
BOOST_FOREACH(edge_descriptor ed, edges(p)){
BOOST_FOREACH(edge_descriptor ed, p_edges){
if(is_border(ed,p)){
++num_be;
border_halfedges.push_back(halfedge(ed,p));
}
}
Point* point_buffer = new Point[num_e*2];
// Build the point_buffer
vertex_iterator vitr, vitr_end;
boost::tie(vitr,vitr_end) = vertices(p);
int pi = 0;
BOOST_FOREACH(vertex_descriptor vd, vertices(p)){
int vi=0;
std::vector<bool> is_border_vertex(num_v, false);
BOOST_FOREACH(vertex_descriptor vd, p_vertices){
BOOST_FOREACH(halfedge_descriptor hd, halfedges_around_target(vd,p)){
if (! is_border(hd,p)){
mask.corner_node(hd, point_buffer[pi++]);
}
else
is_border_vertex[vi]=true;
}
++vi;
}
// Reserve to avoid rellocations during insertions
reserve(p,num_v+num_e+num_f, 2*num_e, (2+4+2)*num_e);
call_reserve(p,num_v+num_e+num_f, 2*num_e, (2+4+2)*num_e);
// Build the connectivity using insert_vertex() and insert_edge()
pi = 0;
for (size_t i = 0; i < num_v; ++i) {
typename std::vector<vertex_descriptor>::iterator vitr=p_vertices.begin();
for (std::size_t i = 0; i < num_v; ++i) {
vertex_descriptor vh = *vitr;
++vitr;
Halfedge_around_target_circulator<Poly> vcir(vh,p);
size_t vn = degree(vh,p);
for (size_t j = 0; j < vn; ++j) {
typename boost::graph_traits<Poly>::degree_size_type vn = degree(vh,p);
for (typename boost::graph_traits<Poly>::degree_size_type j = 0; j < vn; ++j) {
halfedge_descriptor e = *vcir;
++vcir;
if (! is_border(e,p)) {
@ -329,7 +347,7 @@ void DQQ_1step_impl(Poly& p, VertexPointMap vpm, Mask mask, CGAL::Tag_false) {
}
vcir = Halfedge_around_target_circulator<Poly>(vh,p);
for (size_t j = 0; j < vn; ++j) {
for (typename boost::graph_traits<Poly>::vertices_size_type j = 0; j < vn; ++j) {
if (! is_border(*vcir,p)) {
halfedge_descriptor e1 = prev(*vcir, p);
++vcir;
@ -343,8 +361,8 @@ void DQQ_1step_impl(Poly& p, VertexPointMap vpm, Mask mask, CGAL::Tag_false) {
}
}
edge_iterator eitr = edges(p).first;
for (size_t i = 0; i < num_e; ++i) {
typename std::vector<edge_descriptor>::iterator eitr = p_edges.begin();
for (typename boost::graph_traits<Poly>::edges_size_type i = 0; i < num_e; ++i) {
halfedge_descriptor eh = halfedge(*eitr,p);
++eitr;
if (! is_border(edge(eh,p),p)) {
@ -361,7 +379,7 @@ void DQQ_1step_impl(Poly& p, VertexPointMap vpm, Mask mask, CGAL::Tag_false) {
}
}
// After this point, the original border edges are in front!
BOOST_FOREACH(halfedge_descriptor eeh, border_halfedges){
halfedge_descriptor eh = eeh;
if (is_border(eh,p)){
@ -378,169 +396,17 @@ void DQQ_1step_impl(Poly& p, VertexPointMap vpm, Mask mask, CGAL::Tag_false) {
Euler::remove_face(ehe,p);
}
vitr = vertices(p).first;
for (size_t i = 0; i < num_v-num_be; ++i) {
vitr = p_vertices.begin();
for (typename boost::graph_traits<Poly>::vertices_size_type i = 0; i < num_v; ++i) {
vertex_descriptor vh = *vitr;
++vitr;
Euler::remove_center_vertex(halfedge(vh,p),p);
if (!is_border_vertex[i])
Euler::remove_center_vertex(halfedge(vh,p),p);
}
delete []point_buffer;
}
template <class Poly, class VertexPointMap, class Mask>
void DQQ_1step_impl(Poly& p, VertexPointMap vpm, Mask mask, CGAL::Tag_true) {
typedef typename boost::graph_traits<Poly>::vertex_descriptor vertex_descriptor;
typedef typename boost::graph_traits<Poly>::halfedge_descriptor halfedge_descriptor;
typedef typename boost::graph_traits<Poly>::face_descriptor face_descriptor;
typedef typename boost::property_traits<VertexPointMap>::value_type Point;
// Note that for types like 'Surface_mesh', num_vertices() and other similar functions
// return the TOTAL number of elements, which may include removed vertices.
typename boost::graph_traits<Poly>::vertices_size_type num_v = num_vertices(p);
typename boost::graph_traits<Poly>::edges_size_type num_e = num_edges(p);
typename boost::graph_traits<Poly>::faces_size_type num_f = num_faces(p);
// Move `p` into `moved_p`, and build the subdivided mesh from scratch in `p`.
// This is done to make the algorithm work with CGAL::Surface_mesh,
// even though CGAL::Surface_mesh does not insert elements at the end (due to removed elements).
// The DooSabin subdivision of a mesh is a completely different mesh so there
// is no additional cost to rebuild from scratch (but there is a bit from
// using `copy_face_graph`).
Poly moved_p;
reserve(moved_p,num_v, num_e, num_f);
// We must use copy_face_graph rather than an assignement operator because
// we need the correspondence between vertex_descriptors
boost::unordered_map<vertex_descriptor, vertex_descriptor> v2v(num_v);
CGAL::copy_face_graph(p, moved_p, std::inserter(v2v, v2v.end()));
VertexPointMap moved_vpm = get(vertex_point, moved_p);
// Move the position information to the internal property map of moved_p
typename boost::unordered_map<vertex_descriptor, vertex_descriptor>::iterator it = v2v.begin(),
end = v2v.end();
for(; it!=end; ++it) {
put(moved_vpm, it->second, get(vpm, it->first));
}
// Temporarily change the members of the mask to `moved_p`
mask.pmesh = &moved_p;
mask.vpm = moved_vpm;
clear(p);
reserve(p,num_v+num_e+num_f, 2*num_e, (2+4+2)*num_e);
// Correspondence between halfedges of the original mesh and some of the
// halfedges in the subdivided mesh. Since we have halfedge_index_t,
// we can simply use a vector!
// Note: need to make sure the halfedge index map is initialized
// @fixme this overwrites previous initilizations...
helpers::init_halfedge_indices(const_cast<Poly&>(moved_p),
get(boost::halfedge_index, moved_p));
std::vector<halfedge_descriptor> old_to_new(2 * num_e);
Property_map_binder<typename boost::property_map<Poly, boost::halfedge_index_t>::type,
typename Pointer_property_map<halfedge_descriptor>::type>
hmap = bind_property_maps(get(boost::halfedge_index, moved_p),
make_property_map(old_to_new));
// Build new n-faces
BOOST_FOREACH(face_descriptor fd, faces(moved_p)) {
halfedge_descriptor hd = halfedge(fd, moved_p);
std::list<vertex_descriptor> vertices_of_new_face;
// Keep the first outside; it will be used to build the correspondence
// between old and new halfedges
Point first_pt;
mask.corner_node(hd, first_pt);
vertex_descriptor first_v = add_vertex(p);
put(vpm, first_v, first_pt);
vertices_of_new_face.push_back(first_v);
// Loop normally and add the rest of the vertices
halfedge_descriptor done = hd;
hd = next(hd, moved_p);
while(hd != done) {
Point pt;
mask.corner_node(hd, pt);
vertex_descriptor v = add_vertex(p);
put(vpm, v, pt);
vertices_of_new_face.push_back(v);
hd = next(hd, moved_p);
}
face_descriptor new_face = Euler::add_face(vertices_of_new_face, p);
// Find the starting halfedge in the new face that corresponds to halfedge(fd, p)
halfedge_descriptor nf_hd = halfedge(new_face, p);
while(target(nf_hd, p) != first_v) {
nf_hd = next(nf_hd, p);
}
// Build the correspondence between old and new halfedges
hd = halfedge(fd, moved_p);
done = nf_hd;
do {
put(hmap, hd, nf_hd);
hd = next(hd, moved_p);
nf_hd = next(nf_hd, p);
} while (nf_hd != done);
}
// Build new edge-faces
BOOST_FOREACH(halfedge_descriptor hd, halfedges(moved_p)) {
if(is_border(hd, moved_p))
continue;
halfedge_descriptor hd_opp = opposite(hd, moved_p);
if(is_border(hd_opp, moved_p))
continue;
if(hd > hd_opp)
continue;
halfedge_descriptor new_hd = opposite(get(hmap, hd), p);
halfedge_descriptor new_hd_opp = opposite(get(hmap, hd_opp), p);
boost::array<vertex_descriptor, 4> v = {{source(new_hd, p),
target(new_hd, p),
source(new_hd_opp, p),
target(new_hd_opp, p)}};
Euler::add_face(v, p);
}
// Build new vertex-faces
BOOST_FOREACH(vertex_descriptor vd, vertices(moved_p)) {
if(is_border(vd, moved_p))
continue;
halfedge_descriptor hd = halfedge(vd, moved_p);
halfedge_descriptor new_hd = opposite(get(hmap, hd), p);
halfedge_descriptor new_face_hd = opposite(prev(new_hd, p), p), done = new_face_hd;
std::list<vertex_descriptor> vertices_of_new_faces;
do {
vertices_of_new_faces.push_back(source(new_face_hd, p));
new_face_hd = next(new_face_hd, p);
} while(new_face_hd != done);
Euler::add_face(vertices_of_new_faces, p);
}
// Reset the members of the mask
mask.pmesh = &p;
mask.vpm = vpm;
}
template <class Poly, class VertexPointMap, class Mask>
void DQQ_1step(Poly& p, VertexPointMap vpm, Mask mask) {
// Check if halfedges are index-based, which allows to use vectors instead of maps
DQQ_1step_impl(p, vpm, mask,
boost::graph_has_property<Poly, boost::halfedge_index_t>());
// CGAL_postcondition(p.is_valid());
}
// ======================================================================
template <class Poly, class VertexPointMap, class Mask>
void Sqrt3_1step(Poly& p, VertexPointMap vpm, Mask mask,
@ -556,19 +422,21 @@ void Sqrt3_1step(Poly& p, VertexPointMap vpm, Mask mask,
typedef typename boost::graph_traits<Poly>::edge_descriptor edge_descriptor;
typedef typename boost::graph_traits<Poly>::face_descriptor face_descriptor;
typedef typename boost::graph_traits<Poly>::edge_iterator edge_iterator;
typedef typename boost::graph_traits<Poly>::face_iterator face_iterator;
typedef typename boost::property_traits<VertexPointMap>::value_type Point;
typename boost::graph_traits<Poly>::vertices_size_type num_v = num_vertices(p);
typename boost::graph_traits<Poly>::halfedges_size_type num_e = num_halfedges(p)/2;
typename boost::graph_traits<Poly>::faces_size_type num_f = num_faces(p);
//First back up initial vertices/faces/edges
std::vector<vertex_descriptor> p_vertices(vertices(p).first, vertices(p).second);
std::vector<face_descriptor> p_faces(faces(p).first, faces(p).second);
std::vector<edge_descriptor> p_edges(edges(p).first, edges(p).second);
std::size_t num_v = p_vertices.size();
std::size_t num_e = p_edges.size();
std::size_t num_f = p_faces.size();
// reserve enough size for the new points
typename boost::graph_traits<Poly>::faces_size_type new_pts_size = num_f;
std::size_t new_pts_size = num_f;
if(refine_border) {
BOOST_FOREACH(edge_descriptor ed, edges(p)){
BOOST_FOREACH(edge_descriptor ed, p_edges){
if(is_border(ed, p))
++new_pts_size;
}
@ -576,7 +444,7 @@ void Sqrt3_1step(Poly& p, VertexPointMap vpm, Mask mask,
Point* cpt = new Point[new_pts_size];
// size of the subdivided mesh
reserve(p,num_v + new_pts_size, (num_e + 2*num_f + new_pts_size)*2, 3*num_f);
call_reserve(p,num_v + new_pts_size, (num_e + 2*num_f + new_pts_size)*2, 3*num_f);
// keep in memory whether a face is incident to the border and, if so, which
// halfedge corresponds to THE (there can only be one) border edge.
@ -587,7 +455,7 @@ void Sqrt3_1step(Poly& p, VertexPointMap vpm, Mask mask,
// compute the positions of new points
std::size_t i = 0;
std::size_t face_id = 0;
BOOST_FOREACH (face_descriptor fd, faces(p)) {
BOOST_FOREACH (face_descriptor fd, p_faces) {
//ASSERTION_MSG(circulator_size(fitr->facet_begin())==3, "(ERROR) Non-triangle facet!");
if(refine_border) {
BOOST_FOREACH(halfedge_descriptor hd, halfedges_around_face(halfedge(fd, p), p)) {
@ -614,7 +482,7 @@ void Sqrt3_1step(Poly& p, VertexPointMap vpm, Mask mask,
}
// smooth the position of existing vertices
BOOST_FOREACH(vertex_descriptor vd, vertices(p)){
BOOST_FOREACH(vertex_descriptor vd, p_vertices){
Point pt;
if(!is_border(vd, p)) {
mask.vertex_node(vd, pt);
@ -623,10 +491,9 @@ void Sqrt3_1step(Poly& p, VertexPointMap vpm, Mask mask,
}
// insert the new subdividing points
face_iterator b,e;
boost::tie(b,e) = faces(p);
for(std::size_t i=0, cpt_id=0; i < num_f; ++i, ++b){
face_descriptor fd = *b;
typename std::vector<face_descriptor>::iterator fit=p_faces.begin();
for(std::size_t i=0, cpt_id=0; i < num_f; ++i, ++fit){
face_descriptor fd = *fit;
halfedge_descriptor hd = face_halfedge_border[i];
if(refine_border && hd != boost::graph_traits<Poly>::null_halfedge()) {
halfedge_descriptor hd_next = next(hd, p);
@ -650,8 +517,8 @@ void Sqrt3_1step(Poly& p, VertexPointMap vpm, Mask mask,
}
// flip the old edges (except the border edges)
edge_iterator eitr = edges(p).first;
for (size_t i = 0; i < num_e; ++i) {
typename std::vector<edge_descriptor>::iterator eitr = p_edges.begin();
for (std::size_t i = 0; i < num_e; ++i) {
halfedge_descriptor e = halfedge(*eitr,p);
++eitr; // move to next edge before flip since flip destroys current edge
if (! is_border(edge(e,p),p)) {

56
Subdivision_method_3/include/CGAL/Subdivision_method_3/subdivision_masks_3.h
View File

@ -93,6 +93,8 @@ public:
typedef typename Base::FT FT;
typedef typename Base::Point Point;
typedef typename Base::Vector Vector;
typedef typename boost::property_traits<VertexPointMap>::reference Point_ref;
#endif
public:
@ -117,8 +119,8 @@ public:
}
void edge_node(halfedge_descriptor edge, Point& pt) {
const Point& p1 = get(this->vpmap, target(edge, *(this->pmesh)));
const Point& p2 = get(this->vpmap, source(edge, *(this->pmesh)));
Point_ref p1 = get(this->vpmap, target(edge, *(this->pmesh)));
Point_ref p2 = get(this->vpmap, source(edge, *(this->pmesh)));
pt = Point((p1[0]+p2[0])/2, (p1[1]+p2[1])/2, (p1[2]+p2[2])/2);
}
@ -173,6 +175,8 @@ public:
typedef typename Base::FT FT;
typedef typename Base::Point Point;
typedef typename Base::Vector Vector;
typedef typename boost::property_traits<VertexPointMap>::reference Point_ref;
#endif
public:
@ -197,8 +201,8 @@ public:
/// computes the Catmull-Clark edge-point `pt` of the edge `edge`.
void edge_node(halfedge_descriptor edge, Point& pt) {
const Point& p1 = get(this->vpmap,target(edge, *(this->pmesh)));
const Point& p2 = get(this->vpmap,source(edge, *(this->pmesh)));
Point_ref p1 = get(this->vpmap,target(edge, *(this->pmesh)));
Point_ref p2 = get(this->vpmap,source(edge, *(this->pmesh)));
Point f1, f2;
this->face_node(face(edge, *(this->pmesh)), f1);
this->face_node(face(opposite(edge, *(this->pmesh)), *(this->pmesh)), f2);
@ -213,11 +217,11 @@ public:
typename boost::graph_traits<Mesh>::degree_size_type n = degree(vertex, *(this->pmesh));
FT Q[] = {0.0, 0.0, 0.0}, R[] = {0.0, 0.0, 0.0};
Point& S = get(this->vpmap,vertex);
Point_ref S = get(this->vpmap,vertex);
Point q;
for (unsigned int i = 0; i < n; i++, ++vcir) {
const Point& p2 = get(this->vpmap, target(opposite(*vcir, *(this->pmesh)), *(this->pmesh)));
for (typename boost::graph_traits<Mesh>::degree_size_type i = 0; i < n; i++, ++vcir) {
Point_ref p2 = get(this->vpmap, target(opposite(*vcir, *(this->pmesh)), *(this->pmesh)));
R[0] += (S[0] + p2[0]) / 2;
R[1] += (S[1] + p2[1]) / 2;
R[2] += (S[2] + p2[2]) / 2;
@ -237,15 +241,15 @@ public:
/// computes the Catmull-Clark edge-point `ept` and the Catmull-Clark
/// vertex-point `vpt` of the border edge `edge`.
void border_node(halfedge_descriptor edge, Point& ept, Point& vpt) {
const Point& ep1 = get(this->vpmap,target(edge, *(this->pmesh)));
const Point& ep2 = get(this->vpmap,target(opposite(edge, *(this->pmesh)), *(this->pmesh)));
Point_ref ep1 = get(this->vpmap,target(edge, *(this->pmesh)));
Point_ref ep2 = get(this->vpmap,target(opposite(edge, *(this->pmesh)), *(this->pmesh)));
ept = Point((ep1[0]+ep2[0])/2, (ep1[1]+ep2[1])/2, (ep1[2]+ep2[2])/2);
Halfedge_around_target_circulator<Mesh> vcir(edge, *(this->pmesh));
const Point& vp1 = get(this->vpmap,target(opposite(*vcir, *(this->pmesh)), *(this->pmesh)));
const Point& vp0 = get(this->vpmap, target(*vcir, *(this->pmesh)));
Point_ref vp1 = get(this->vpmap,target(opposite(*vcir, *(this->pmesh)), *(this->pmesh)));
Point_ref vp0 = get(this->vpmap, target(*vcir, *(this->pmesh)));
--vcir;
const Point& vp_1 = get(this->vpmap, target(opposite(*vcir, *(this->pmesh)), *(this->pmesh)));
Point_ref vp_1 = get(this->vpmap, target(opposite(*vcir, *(this->pmesh)), *(this->pmesh)));
vpt = Point((vp_1[0] + 6*vp0[0] + vp1[0])/8,
(vp_1[1] + 6*vp0[1] + vp1[1])/8,
(vp_1[2] + 6*vp0[2] + vp1[2])/8 );
@ -296,6 +300,7 @@ public:
typedef typename Base::FT FT;
typedef typename Base::Point Point;
typedef typename Base::Vector Vector;
typedef typename boost::property_traits<VertexPointMap>::reference Point_ref;
#endif
typedef Halfedge_around_face_circulator<Mesh> Halfedge_around_facet_circulator;
@ -323,10 +328,10 @@ public:
/// computes the Loop edge-point `pt` of the edge `edge`.
void edge_node(halfedge_descriptor edge, Point& pt) {
const Point& p1 = get(this->vpmap,target(edge, *(this->pmesh)));
const Point& p2 = get(this->vpmap,target(opposite(edge, *(this->pmesh)), *(this->pmesh)));
const Point& f1 = get(this->vpmap,target(next(edge, *(this->pmesh)), *(this->pmesh)));
const Point& f2 = get(this->vpmap,target(next(opposite(edge, *(this->pmesh)), *(this->pmesh)), *(this->pmesh)));
Point_ref p1 = get(this->vpmap,target(edge, *(this->pmesh)));
Point_ref p2 = get(this->vpmap,target(opposite(edge, *(this->pmesh)), *(this->pmesh)));
Point_ref f1 = get(this->vpmap,target(next(edge, *(this->pmesh)), *(this->pmesh)));
Point_ref f2 = get(this->vpmap,target(next(opposite(edge, *(this->pmesh)), *(this->pmesh)), *(this->pmesh)));
pt = Point((3*(p1[0]+p2[0])+f1[0]+f2[0])/8,
(3*(p1[1]+p2[1])+f1[1]+f2[1])/8,
@ -339,10 +344,10 @@ public:
size_t n = circulator_size(vcir);
FT R[] = {0.0, 0.0, 0.0};
const Point& S = get(this->vpmap,vertex);
Point_ref S = get(this->vpmap,vertex);
for (size_t i = 0; i < n; i++, ++vcir) {
const Point& p = get(this->vpmap,target(opposite(*vcir, *(this->pmesh)), *(this->pmesh)));
Point_ref p = get(this->vpmap,target(opposite(*vcir, *(this->pmesh)), *(this->pmesh)));
R[0] += p[0]; R[1] += p[1]; R[2] += p[2];
}
if (n == 6) {
@ -361,15 +366,15 @@ public:
/// computes the Loop edge-point `ept` and the Loop vertex-point `vpt` of the border edge `edge`.
void border_node(halfedge_descriptor edge, Point& ept, Point& vpt) {
const Point& ep1 = get(this->vpmap,target(edge, *(this->pmesh)));
const Point& ep2 = get(this->vpmap,target(opposite(edge, *(this->pmesh)), *(this->pmesh)));
Point_ref ep1 = get(this->vpmap,target(edge, *(this->pmesh)));
Point_ref ep2 = get(this->vpmap,target(opposite(edge, *(this->pmesh)), *(this->pmesh)));
ept = Point((ep1[0]+ep2[0])/2, (ep1[1]+ep2[1])/2, (ep1[2]+ep2[2])/2);
Halfedge_around_vertex_circulator vcir(edge, *(this->pmesh));
const Point& vp1 = get(this->vpmap,target(opposite(*vcir, *(this->pmesh) ), *(this->pmesh)));
const Point& vp0 = get(this->vpmap,target(*vcir, *(this->pmesh)));
Point_ref vp1 = get(this->vpmap,target(opposite(*vcir, *(this->pmesh) ), *(this->pmesh)));
Point_ref vp0 = get(this->vpmap,target(*vcir, *(this->pmesh)));
--vcir;
const Point& vp_1 = get(this->vpmap,target(opposite(*vcir, *(this->pmesh)), *(this->pmesh)));
Point_ref vp_1 = get(this->vpmap,target(opposite(*vcir, *(this->pmesh)), *(this->pmesh)));
vpt = Point((vp_1[0] + 6*vp0[0] + vp1[0])/8,
(vp_1[1] + 6*vp0[1] + vp1[1])/8,
(vp_1[2] + 6*vp0[2] + vp1[2])/8 );
@ -457,6 +462,7 @@ public:
typedef typename Base::FT FT;
typedef typename Base::Point Point;
typedef typename Base::Vector Vector;
typedef typename boost::property_traits<VertexPointMap>::reference Point_ref;
#endif
public:
@ -574,12 +580,12 @@ public:
/// computes the \f$ \sqrt{3}\f$ vertex-point `pt` of the vertex `vd`.
void vertex_node(vertex_descriptor vertex, Point& pt) {
Halfedge_around_target_circulator<Mesh> vcir(vertex, *(this->pmesh));
const size_t n = degree(vertex, *(this->pmesh));
const typename boost::graph_traits<Mesh>::degree_size_type n = degree(vertex, *(this->pmesh));
const FT a = (FT) ((4.0-2.0*std::cos(2.0*CGAL_PI/(double)n))/9.0);
Vector cv = ((FT)(1.0-a)) * (get(this->vpmap, vertex) - CGAL::ORIGIN);
for (size_t i = 1; i <= n; ++i, --vcir) {
for (typename boost::graph_traits<Mesh>::degree_size_type i = 1; i <= n; ++i, --vcir) {
cv = cv + (a/FT(n))*(get(this->vpmap, target(opposite(*vcir, *(this->pmesh)), *(this->pmesh)))-CGAL::ORIGIN);
}

2
Subdivision_method_3/test/Subdivision_method_3/test_Subdivision_method_3.cpp
View File

@ -338,7 +338,7 @@ void test_Subdivision_surface_3_SM_NP() {
// some arbitrary new coordinates (different from the internal vpm)
BOOST_FOREACH(vertex_descriptor vd, vertices(P)) {
const Point& pt = get(vpm, vd);
boost::property_traits<VPM>::reference pt = get(vpm, vd);
Vector v = pt - Point(0., 0., -3.);
put(apm, vd, pt + 0.5*v);
}

2
Surface_mesh_shortest_path/include/CGAL/Surface_mesh_shortest_path/Surface_mesh_shortest_path.h
View File

@ -1269,7 +1269,7 @@ private:
propagateRight = rightSide;
}
if (node->level() <= num_faces(m_graph))
if (node->level() <= static_cast<std::size_t>(num_faces(m_graph)))
{
if (propagateLeft)
{

4
Triangulation_2/include/CGAL/Constrained_Delaunay_triangulation_2.h
View File

@ -750,9 +750,9 @@ propagating_flip(Face_handle f,int i, int depth)
Face_handle ni = f->neighbor(i);
flip(f, i); // flip for constrained triangulations
propagating_flip(f,i);
propagating_flip(f,i, depth+1);
i = ni->index(f->vertex(i));
propagating_flip(ni,i);
propagating_flip(ni,i, depth+1);
#endif
}
#else

72
Triangulation_2/include/CGAL/Constrained_triangulation_2.h
View File

@ -691,45 +691,53 @@ insert_constraint(Vertex_handle vaa, Vertex_handle vbb)
// if a vertex vc of t lies on segment ab
// or if ab intersect some constrained edges
{
CGAL_triangulation_precondition( vaa != vbb);
Vertex_handle vi;
std::stack<std::pair<Vertex_handle, Vertex_handle> > stack;
stack.push(std::make_pair(vaa,vbb));
Face_handle fr;
int i;
if(includes_edge(vaa,vbb,vi,fr,i)) {
mark_constraint(fr,i);
if (vi != vbb) {
insert_constraint(vi,vbb);
while(! stack.empty()){
boost::tie(vaa,vbb) = stack.top();
stack.pop();
CGAL_triangulation_precondition( vaa != vbb);
Vertex_handle vi;
Face_handle fr;
int i;
if(includes_edge(vaa,vbb,vi,fr,i)) {
mark_constraint(fr,i);
if (vi != vbb) {
stack.push(std::make_pair(vi,vbb));
}
continue;
}
return;
}
List_faces intersected_faces;
List_edges conflict_boundary_ab, conflict_boundary_ba;
List_faces intersected_faces;
List_edges conflict_boundary_ab, conflict_boundary_ba;
bool intersection = find_intersected_faces( vaa, vbb,
intersected_faces,
conflict_boundary_ab,
conflict_boundary_ba,
vi);
if ( intersection) {
if (vi != vaa && vi != vbb) {
insert_constraint(vaa,vi);
insert_constraint(vi,vbb);
}
else insert_constraint(vaa,vbb);
return;
}
bool intersection = find_intersected_faces( vaa, vbb,
intersected_faces,
conflict_boundary_ab,
conflict_boundary_ba,
vi);
if ( intersection) {
if (vi != vaa && vi != vbb) {
stack.push(std::make_pair(vaa,vi));
stack.push(std::make_pair(vi,vbb));
}
else{
stack.push(std::make_pair(vaa,vbb));
}
continue;
}
//no intersection
triangulate_hole(intersected_faces,
conflict_boundary_ab,
conflict_boundary_ba);
//no intersection
triangulate_hole(intersected_faces,
conflict_boundary_ab,
conflict_boundary_ba);
if (vi != vbb) {
insert_constraint(vi,vbb);
if (vi != vbb) {
stack.push(std::make_pair(vi,vbb));
}
}
return;
}

124
Triangulation_2/include/CGAL/Constrained_triangulation_plus_2.h
View File

@ -781,82 +781,88 @@ insert_subconstraint(Vertex_handle vaa,
// insert the subconstraint [vaa vbb]
// it will eventually be splitted into several subconstraints
{
CGAL_triangulation_precondition( vaa != vbb);
Vertex_handle vi;
std::stack<std::pair<Vertex_handle, Vertex_handle> > stack;
stack.push(std::make_pair(vaa,vbb));
Face_handle fr;
int i;
if(this->includes_edge(vaa,vbb,vi,fr,i)) {
this->mark_constraint(fr,i);
if (vi != vbb) {
hierarchy.split_constraint(vaa,vbb,vi);
insert_subconstraint(vi,vbb, out);
while(! stack.empty()){
boost::tie(vaa,vbb) = stack.top();
stack.pop();
CGAL_triangulation_precondition( vaa != vbb);
Vertex_handle vi;
Face_handle fr;
int i;
if(this->includes_edge(vaa,vbb,vi,fr,i)) {
this->mark_constraint(fr,i);
if (vi != vbb) {
hierarchy.split_constraint(vaa,vbb,vi);
stack.push(std::make_pair(vi,vbb));
}
continue;
}
return;
}
List_faces intersected_faces;
List_edges conflict_boundary_ab, conflict_boundary_ba;
List_faces intersected_faces;
List_edges conflict_boundary_ab, conflict_boundary_ba;
bool intersection = this->find_intersected_faces(
vaa, vbb,
intersected_faces,
conflict_boundary_ab,
conflict_boundary_ba,
vi);
bool intersection = this->find_intersected_faces(
vaa, vbb,
intersected_faces,
conflict_boundary_ab,
conflict_boundary_ba,
vi);
if ( intersection) {
if (vi != vaa && vi != vbb) {
hierarchy.split_constraint(vaa,vbb,vi);
insert_subconstraint(vaa,vi, out);
insert_subconstraint(vi,vbb, out);
}
else insert_subconstraint(vaa,vbb,out);
if ( intersection) {
if (vi != vaa && vi != vbb) {
hierarchy.split_constraint(vaa,vbb,vi);
stack.push(std::make_pair(vaa,vi));
stack.push(std::make_pair(vi,vbb));
}
else stack.push(std::make_pair(vaa,vbb));
return;
}
continue;
}
//no intersection
//no intersection
List_edges edges(conflict_boundary_ab);
std::copy(conflict_boundary_ba.begin(), conflict_boundary_ba.end(), std::back_inserter(edges));
List_edges edges(conflict_boundary_ab);
std::copy(conflict_boundary_ba.begin(), conflict_boundary_ba.end(), std::back_inserter(edges));
// edges may contain mirror edges. They no longer exist after triangulate_hole
// so we have to remove them before calling get_bounded_faces
if(! edges.empty()){
// edges may contain mirror edges. They no longer exist after triangulate_hole
// so we have to remove them before calling get_bounded_faces
if(! edges.empty()){
#if defined(BOOST_MSVC) && (BOOST_VERSION == 105500)
std::set<Face_handle> faces(intersected_faces.begin(), intersected_faces.end());
std::set<Face_handle> faces(intersected_faces.begin(), intersected_faces.end());
#else
boost::container::flat_set<Face_handle> faces(intersected_faces.begin(), intersected_faces.end());
boost::container::flat_set<Face_handle> faces(intersected_faces.begin(), intersected_faces.end());
#endif
typename List_edges::iterator it2;
for(typename List_edges::iterator it = edges.begin(); it!= edges.end();){
if(faces.find(it->first) != faces.end()){
typename List_edges::iterator it2 = it;
++it;
edges.erase(it2);
}else {
++it;
typename List_edges::iterator it2;
for(typename List_edges::iterator it = edges.begin(); it!= edges.end();){
if(faces.find(it->first) != faces.end()){
typename List_edges::iterator it2 = it;
++it;
edges.erase(it2);
}else {
++it;
}
}
}
this->triangulate_hole(intersected_faces,
conflict_boundary_ab,
conflict_boundary_ba);
this->get_bounded_faces(edges.begin(),
edges.end(),
out);
if (vi != vbb) {
hierarchy.split_constraint(vaa,vbb,vi);
stack.push(std::make_pair(vi,vbb));
}
}
this->triangulate_hole(intersected_faces,
conflict_boundary_ab,
conflict_boundary_ba);
this->get_bounded_faces(edges.begin(),
edges.end(),
out);
if (vi != vbb) {
hierarchy.split_constraint(vaa,vbb,vi);
insert_subconstraint(vi,vbb, out);
}
return;
}