Merge remote-tracking branch 'cgal/master' into PMP-Halfedge_comparer_for_stitch_borders-GF

This commit is contained in:
Mael Rouxel-Labbé
2020-05-12 19:19:41 +02:00
6402 changed files with 271147 additions and 198533 deletions
@@ -6,7 +6,7 @@
// $URL$
// $Id$
// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
//
//
//
// Author(s) : Yin Xu, Andreas Fabri and Ilker O. Yaz
@@ -39,10 +39,10 @@ struct Cotangent_value_Meyer_impl
{
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();
@@ -52,17 +52,17 @@ struct Cotangent_value_Meyer_impl
Vector cross_ab = CGAL::cross_product(a, b);
double divider = to_double(CGAL::approximate_sqrt(cross_ab*cross_ab));
if(divider == 0 /*|| divider != divider*/)
if(divider == 0 /*|| divider != divider*/)
{
CGAL::collinear(get(ppmap, v0), get(ppmap, v1), get(ppmap, v2)) ?
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;
}
};
@@ -384,7 +384,7 @@ public:
/////////////////////////////////////////////////////////////////////////////////////////
///////////////////////////// Edge Weight Calculators ///////////////////////////////////
// Cotangent weight calculator
// 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,
@@ -404,7 +404,7 @@ struct Cotangent_weight_impl : CotangentValue
// 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) )
@@ -417,7 +417,7 @@ struct Cotangent_weight_impl : CotangentValue
else
{
halfedge_descriptor he_cw = opposite( next(he, pmesh) , pmesh );
vertex_descriptor v2 = source(he_cw, pmesh);
vertex_descriptor v2 = source(he_cw, pmesh);
halfedge_descriptor he_ccw = prev( opposite(he, pmesh) , pmesh );
vertex_descriptor v3 = source(he_ccw, pmesh);
@@ -471,7 +471,7 @@ public:
// 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()) )
@@ -484,7 +484,7 @@ public:
else
{
halfedge_descriptor he_cw = opposite( next(he, pmesh()) , pmesh() );
vertex_descriptor v2 = source(he_cw, pmesh());
vertex_descriptor v2 = source(he_cw, pmesh());
halfedge_descriptor he_ccw = prev( opposite(he, pmesh()) , pmesh() );
vertex_descriptor v3 = source(he_ccw, pmesh());
@@ -507,7 +507,7 @@ struct Single_cotangent_weight_impl : CotangentValue
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);
@@ -550,7 +550,7 @@ public:
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());
@@ -688,7 +688,7 @@ public:
else
{
halfedge_descriptor he_cw = opposite( next(he, pmesh()) , pmesh() );
vertex_descriptor v2 = source(he_cw, pmesh());
vertex_descriptor v2 = source(he_cw, pmesh());
halfedge_descriptor he_ccw = prev( opposite(he, pmesh()) , pmesh() );
vertex_descriptor v3 = source(he_ccw, pmesh());
@@ -711,7 +711,7 @@ private:
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) );
}
@@ -734,7 +734,7 @@ private:
}
};
template< class PolygonMesh,
template< class PolygonMesh,
class PrimaryWeight = Cotangent_weight<PolygonMesh>,
class SecondaryWeight = Mean_value_weight<PolygonMesh> >
class Hybrid_weight : public PrimaryWeight, SecondaryWeight
@@ -68,12 +68,13 @@ namespace CGAL {
{
using parameters::choose_parameter;
using parameters::get_parameter;
typename GetVertexPointMap<PolygonMesh, CGAL_PMP_NP_CLASS>::const_type
vpm = choose_parameter(get_parameter(np, internal_np::vertex_point),
get_const_property_map(CGAL::vertex_point, pmesh));
get_const_property_map(CGAL::vertex_point, pmesh));
typedef typename GetGeomTraits<PolygonMesh, CGAL_PMP_NP_CLASS>::type GT;
GT gt = choose_parameter(get_parameter(np, internal_np::geom_traits), GT());
GT gt = choose_parameter<GT>(get_parameter(np, internal_np::geom_traits));
typename GT::Construct_bbox_3 get_bbox = gt.construct_bbox_3_object();
typedef typename boost::graph_traits<PolygonMesh>::vertex_descriptor vertex_descriptor;
@@ -120,10 +121,10 @@ namespace CGAL {
using parameters::get_parameter;
typename GetVertexPointMap<PolygonMesh, NamedParameters>::const_type
vpm = choose_parameter(get_parameter(np, internal_np::vertex_point),
get_const_property_map(CGAL::vertex_point, pmesh));
get_const_property_map(CGAL::vertex_point, pmesh));
typedef typename GetGeomTraits<PolygonMesh, NamedParameters>::type GT;
GT gt = choose_parameter(get_parameter(np, internal_np::geom_traits), GT());
GT gt = choose_parameter<GT>(get_parameter(np, internal_np::geom_traits));
typename GT::Construct_bbox_3 get_bbox = gt.construct_bbox_3_object();
return get_bbox( get(vpm, vd) );
@@ -161,12 +162,13 @@ namespace CGAL {
{
using parameters::choose_parameter;
using parameters::get_parameter;
typename GetVertexPointMap<PolygonMesh, NamedParameters>::const_type
vpm = choose_parameter(get_parameter(np, internal_np::vertex_point),
get_const_property_map(CGAL::vertex_point, pmesh));
get_const_property_map(CGAL::vertex_point, pmesh));
typedef typename GetGeomTraits<PolygonMesh, NamedParameters>::type GT;
GT gt = choose_parameter(get_parameter(np, internal_np::geom_traits), GT());
GT gt = choose_parameter<GT>(get_parameter(np, internal_np::geom_traits));
typename GT::Construct_bbox_3 get_bbox = gt.construct_bbox_3_object();
return get_bbox( get(vpm, source(ed, pmesh)) ) +
@@ -205,12 +207,13 @@ namespace CGAL {
{
using parameters::choose_parameter;
using parameters::get_parameter;
typename GetVertexPointMap<PolygonMesh, NamedParameters>::const_type
vpm = choose_parameter(get_parameter(np, internal_np::vertex_point),
get_const_property_map(CGAL::vertex_point, pmesh));
get_const_property_map(CGAL::vertex_point, pmesh));
typedef typename GetGeomTraits<PolygonMesh, NamedParameters>::type GT;
GT gt = choose_parameter(get_parameter(np, internal_np::geom_traits), GT());
GT gt = choose_parameter<GT>(get_parameter(np, internal_np::geom_traits));
typename GT::Construct_bbox_3 get_bbox = gt.construct_bbox_3_object();
typedef typename boost::graph_traits<PolygonMesh>::halfedge_descriptor halfedge_descriptor;
@@ -53,7 +53,7 @@ std::size_t border_size(typename boost::graph_traits<PolygonMesh>::halfedge_desc
template<typename PM
, typename FaceRange
, typename HalfedgeOutputIterator>
HalfedgeOutputIterator border_halfedges_impl(const FaceRange& faces
HalfedgeOutputIterator border_halfedges_impl(const FaceRange& face_range
, HalfedgeOutputIterator out
, const PM& pmesh)
{
@@ -64,7 +64,7 @@ std::size_t border_size(typename boost::graph_traits<PolygonMesh>::halfedge_desc
// the bool is true if the halfedge stored is the one of the face,
// false if it is its opposite
std::map<halfedge_descriptor, bool> border;
for(face_descriptor f : faces)
for(face_descriptor f : face_range)
{
for(halfedge_descriptor h :
halfedges_around_face(halfedge(f, pmesh), pmesh))
@@ -94,13 +94,13 @@ std::size_t border_size(typename boost::graph_traits<PolygonMesh>::halfedge_desc
, typename FaceRange
, typename HalfedgeOutputIterator
, typename NamedParameters>
HalfedgeOutputIterator border_halfedges_impl(const FaceRange& faces
HalfedgeOutputIterator border_halfedges_impl(const FaceRange& face_range
, typename boost::cgal_no_property::type
, HalfedgeOutputIterator out
, const PM& pmesh
, const NamedParameters& /* np */)
{
return border_halfedges_impl(faces, out, pmesh);
return border_halfedges_impl(face_range, out, pmesh);
}
template<typename PM
@@ -108,7 +108,7 @@ std::size_t border_size(typename boost::graph_traits<PolygonMesh>::halfedge_desc
, typename FaceIndexMap
, typename HalfedgeOutputIterator
, typename NamedParameters>
HalfedgeOutputIterator border_halfedges_impl(const FaceRange& faces
HalfedgeOutputIterator border_halfedges_impl(const FaceRange& face_range
, const FaceIndexMap& fmap
, HalfedgeOutputIterator out
, const PM& pmesh
@@ -117,25 +117,13 @@ std::size_t border_size(typename boost::graph_traits<PolygonMesh>::halfedge_desc
typedef typename boost::graph_traits<PM>::halfedge_descriptor halfedge_descriptor;
typedef typename boost::graph_traits<PM>::face_descriptor face_descriptor;
//make a minimal check that it's properly initialized :
//if the 2 first faces have the same id, we know the property map is not initialized
if (boost::is_same<typename GetFaceIndexMap<PM, NamedParameters>::Is_internal_map,
boost::true_type>::value)
{
typename boost::range_iterator<const FaceRange>::type it = boost::const_begin(faces);
if (get(fmap, *it) == get(fmap, *std::next(it)))
{
std::cerr << "WARNING : the internal property map for CGAL::face_index_t" << std::endl
<< " is not properly initialized." << std::endl
<< " Initialize it before calling border_halfedges()" << std::endl;
}
}
CGAL_assertion(BGL::internal::is_index_map_valid(fmap, num_faces(pmesh), faces(pmesh)));
std::vector<bool> present(num_faces(pmesh), false);
for(face_descriptor fd : faces)
for(face_descriptor fd : face_range)
present[get(fmap, fd)] = true;
for(face_descriptor fd : faces)
for(face_descriptor fd : face_range)
for(halfedge_descriptor hd :
halfedges_around_face(halfedge(fd, pmesh), pmesh))
{
@@ -161,21 +149,16 @@ std::size_t border_size(typename boost::graph_traits<PolygonMesh>::halfedge_desc
* For each returned halfedge `h`, `opposite(h, pmesh)` belongs to a face of the patch,
* but `face(h, pmesh)` does not belong to the patch.
*
* @tparam PolygonMesh model of `HalfedgeGraph`. If `PolygonMesh`
* has an internal property map
* for `CGAL::face_index_t` and no `face_index_map` is given
* as a named parameter, then the internal one must be initialized
* @tparam FaceRange range of
`boost::graph_traits<PolygonMesh>::%face_descriptor`, model of `Range`.
Its iterator type is `InputIterator`.
* @tparam PolygonMesh model of `HalfedgeGraph`
* @tparam FaceRange a model of `Range` with value type `boost::graph_traits<PolygonMesh>::%face_descriptor`.
* @tparam HalfedgeOutputIterator model of `OutputIterator`
holding `boost::graph_traits<PolygonMesh>::%halfedge_descriptor`
for patch border
* @tparam NamedParameters a sequence of \ref pmp_namedparameters "Named Parameters"
*
* @param pmesh the polygon mesh to which `faces` belong
* @param faces the range of faces defining the patch whose border halfedges
* are collected
* @param pmesh the polygon mesh to which the faces in `face_range` belong
* @param face_range the range of faces defining the patch whose border halfedges
* are collected
* @param out the output iterator that collects the border halfedges of the patch,
* seen from outside.
* @param np optional sequence of \ref pmp_namedparameters "Named Parameters" among the ones listed below
@@ -190,29 +173,18 @@ std::size_t border_size(typename boost::graph_traits<PolygonMesh>::halfedge_desc
, typename FaceRange
, typename HalfedgeOutputIterator
, typename NamedParameters>
HalfedgeOutputIterator border_halfedges(const FaceRange& faces
HalfedgeOutputIterator border_halfedges(const FaceRange& face_range
, const PolygonMesh& pmesh
, HalfedgeOutputIterator out
, const NamedParameters& np)
{
if (faces.empty()) return out;
if (face_range.empty())
return out;
typedef PolygonMesh PM;
typedef typename GetFaceIndexMap<PM, NamedParameters>::const_type FIMap;
typedef typename boost::property_map<typename internal::Dummy_PM,
CGAL::face_index_t>::type Unset_FIMap;
typedef typename CGAL::GetInitializedFaceIndexMap<PolygonMesh, NamedParameters>::const_type FIMap;
FIMap fim = CGAL::get_initialized_face_index_map(pmesh, np);
if (boost::is_same<FIMap, Unset_FIMap>::value || faces.size() == 1)
{
//face index map is not given in named parameters, nor as an internal property map
return internal::border_halfedges_impl(faces, out, pmesh);
}
//face index map given as a named parameter, or as an internal property map
FIMap fim = parameters::choose_parameter(parameters::get_parameter(np, internal_np::face_index),
get_const_property_map(CGAL::face_index, pmesh));
return internal::border_halfedges_impl(faces, fim, out, pmesh, np);
return internal::border_halfedges_impl(face_range, fim, out, pmesh, np);
}
template<typename PolygonMesh
@@ -231,11 +203,11 @@ std::size_t border_size(typename boost::graph_traits<PolygonMesh>::halfedge_desc
template<typename PolygonMesh
, typename FaceRange
, typename HalfedgeOutputIterator>
HalfedgeOutputIterator border_halfedges(const FaceRange& faces
HalfedgeOutputIterator border_halfedges(const FaceRange& face_range
, const PolygonMesh& pmesh
, HalfedgeOutputIterator out)
{
return border_halfedges(faces, pmesh, out,
return border_halfedges(face_range, pmesh, out,
CGAL::Polygon_mesh_processing::parameters::all_default());
}
@@ -244,7 +216,7 @@ std::size_t border_size(typename boost::graph_traits<PolygonMesh>::halfedge_desc
//
// @tparam PolygonMesh model of `HalfedgeGraph`.
//
// @param pmesh the polygon mesh to which `faces` belong
// @param pmesh the polygon mesh to which `face_range` belong
//
template<typename PolygonMesh>
unsigned int number_of_borders(const PolygonMesh& pmesh)
@@ -15,15 +15,25 @@
#include <CGAL/license/Polygon_mesh_processing/corefinement.h>
#include <CGAL/Polygon_mesh_processing/corefinement.h>
#include <CGAL/Polygon_mesh_processing/connected_components.h>
#include <CGAL/Polygon_mesh_processing/bbox.h>
#include <CGAL/Polygon_mesh_processing/triangulate_faces.h>
#include <CGAL/Polygon_mesh_processing/orientation.h>
#include <CGAL/Polygon_mesh_processing/triangulate_hole.h>
#include <CGAL/Polygon_mesh_processing/border.h>
#include <CGAL/Polygon_mesh_processing/repair.h>
#include <CGAL/iterator.h>
#include <CGAL/AABB_triangle_primitive.h>
#include <CGAL/boost/graph/properties.h>
#include <CGAL/boost/graph/Face_filtered_graph.h>
#include <boost/property_map/property_map.hpp>
#include <unordered_map>
namespace CGAL{
namespace Polygon_mesh_processing {
@@ -40,10 +50,10 @@ inter_pt_index(int i, int j,
std::pair<std::map<std::pair<int,int>, int>::iterator, bool> res =
id_map.insert(std::make_pair(make_sorted_pair(i,j),
static_cast<int> (points.size())));
if (res.second)
if(res.second)
points.push_back(
typename Geom_traits::Construct_plane_line_intersection_point_3()
(plane, points[i], points[j]) );
(plane, points[i], points[j]));
return res.first->second;
}
@@ -55,7 +65,7 @@ Oriented_side
clip_to_bbox(const Plane_3& plane,
const Bbox_3& bbox,
TriangleMesh& tm_out,
const NamedParameters& np )
const NamedParameters& np)
{
typedef typename GetGeomTraits<TriangleMesh, NamedParameters>::type Geom_traits;
typedef typename Geom_traits::Point_3 Point_3;
@@ -63,7 +73,7 @@ clip_to_bbox(const Plane_3& plane,
NamedParameters>::type Vpm;
Vpm vpm_out = parameters::choose_parameter(parameters::get_parameter(np, internal_np::vertex_point),
get_property_map(boost::vertex_point, tm_out));
get_property_map(boost::vertex_point, tm_out));
std::vector<Point_3> corners(8);
@@ -103,9 +113,9 @@ clip_to_bbox(const Plane_3& plane,
std::set<int> in_point_ids; // to collect the set of points in the clipped bbox
// for each face of the bbox, we look for intersection of the plane with its edges
for (int i=0; i<6; ++i)
for(int i=0; i<6; ++i)
{
for (int k=0; k< 4; ++k)
for(int k=0; k< 4; ++k)
{
int current_id = face_indices[4*i + k];
int next_id = face_indices[4*i + (k+1)%4];
@@ -116,14 +126,14 @@ clip_to_bbox(const Plane_3& plane,
{
all_out=false;
// point on or on the negative side
output_faces[i].push_back( current_id );
in_point_ids.insert( output_faces[i].back() );
output_faces[i].push_back(current_id);
in_point_ids.insert(output_faces[i].back());
// check for intersection of the edge
if (orientations[ next_id ] == ON_POSITIVE_SIDE)
if(orientations[ next_id ] == ON_POSITIVE_SIDE)
{
output_faces[i].push_back(
inter_pt_index<Geom_traits>(current_id, next_id, plane, corners, id_map) );
in_point_ids.insert( output_faces[i].back() );
inter_pt_index<Geom_traits>(current_id, next_id, plane, corners, id_map));
in_point_ids.insert(output_faces[i].back());
}
break;
}
@@ -131,31 +141,31 @@ clip_to_bbox(const Plane_3& plane,
{
all_in = false;
// check for intersection of the edge
if ( orientations[ next_id ] == ON_NEGATIVE_SIDE )
if(orientations[ next_id ] == ON_NEGATIVE_SIDE)
{
output_faces[i].push_back(
inter_pt_index<Geom_traits>(current_id, next_id, plane, corners, id_map) );
in_point_ids.insert( output_faces[i].back() );
inter_pt_index<Geom_traits>(current_id, next_id, plane, corners, id_map));
in_point_ids.insert(output_faces[i].back());
}
break;
}
case ON_ORIENTED_BOUNDARY:
{
output_faces[i].push_back( current_id );
in_point_ids.insert( output_faces[i].back() );
output_faces[i].push_back(current_id);
in_point_ids.insert(output_faces[i].back());
}
}
}
if (output_faces[i].size() < 3){
if(output_faces[i].size() < 3){
CGAL_assertion(output_faces[i].empty() ||
(output_faces[i].front()<8 && output_faces[i].back()<8) );
(output_faces[i].front()<8 && output_faces[i].back()<8));
output_faces[i].clear(); // edge of the bbox included in the plane
}
}
// the intersection is the full bbox
if (all_in) return ON_NEGATIVE_SIDE;
if (all_out) return ON_POSITIVE_SIDE;
if(all_in) return ON_NEGATIVE_SIDE;
if(all_out) return ON_POSITIVE_SIDE;
// build the clipped bbox
typedef boost::graph_traits<TriangleMesh> graph_traits;
@@ -167,7 +177,7 @@ clip_to_bbox(const Plane_3& plane,
for(int i : in_point_ids)
{
vertex_descriptor v = add_vertex(tm_out);
out_vertices.insert( std::make_pair(i, v ) );
out_vertices.insert(std::make_pair(i, v));
put(vpm_out, v, corners[i]);
}
@@ -176,7 +186,7 @@ clip_to_bbox(const Plane_3& plane,
const face_descriptor null_fd = graph_traits::null_face();
for(const std::vector<int>& findices : output_faces)
{
if (findices.empty()) continue;
if(findices.empty()) continue;
const face_descriptor fd=add_face(tm_out);
int prev_id = findices.back();
@@ -189,12 +199,12 @@ clip_to_bbox(const Plane_3& plane,
std::pair<typename std::map< std::pair<int,int>,
halfedge_descriptor>::iterator, bool> res =
hedge_map.insert( std::make_pair(std::make_pair(prev_id, current_id), null_hedge) );
if (res.second)
hedge_map.insert(std::make_pair(std::make_pair(prev_id, current_id), null_hedge));
if(res.second)
{
res.first->second = halfedge( add_edge(tm_out), tm_out);
hedge_map.insert( std::make_pair(std::make_pair(current_id, prev_id),
opposite(res.first->second, tm_out) ) );
res.first->second = halfedge(add_edge(tm_out), tm_out);
hedge_map.insert(std::make_pair(std::make_pair(current_id, prev_id),
opposite(res.first->second, tm_out)));
set_face(opposite(res.first->second, tm_out), null_fd, tm_out);
}
@@ -208,7 +218,7 @@ clip_to_bbox(const Plane_3& plane,
set_halfedge(src, opposite(hedges.back(), tm_out), tm_out);
set_halfedge(tgt, hedges.back(), tm_out);
if (current_id==findices.front())
if(current_id==findices.front())
set_halfedge(fd, hedges.back(), tm_out);
prev_id = current_id;
@@ -230,7 +240,7 @@ clip_to_bbox(const Plane_3& plane,
// until we reach another border halfedge
for(halfedge_descriptor h : halfedges(tm_out))
{
if (face(h, tm_out) == null_fd)
if(face(h, tm_out) == null_fd)
{
face_descriptor fd = add_face(tm_out);
set_halfedge(fd, h, tm_out);
@@ -244,7 +254,7 @@ clip_to_bbox(const Plane_3& plane,
} while(face(h_curr, tm_out) != null_fd && h_curr!=h);
set_face(h_prev, fd, tm_out);
set_next(h_prev, h_curr, tm_out);
if (h_curr==h)
if(h_curr==h)
break;
h_prev=h_curr;
} while(true);
@@ -259,6 +269,137 @@ clip_to_bbox(const Plane_3& plane,
return ON_ORIENTED_BOUNDARY;
}
template <class TriangleMesh, class Ecm, class VPM>
void split_along_edges(TriangleMesh& tm,
Ecm ecm,
VPM vpm)
{
typedef boost::graph_traits<TriangleMesh> GT;
typedef typename GT::face_descriptor face_descriptor;
typedef typename GT::edge_descriptor edge_descriptor;
typedef typename GT::vertex_descriptor vertex_descriptor;
typedef typename GT::halfedge_descriptor halfedge_descriptor;
std::vector<edge_descriptor> shared_edges;
for(edge_descriptor e : edges(tm))
if(get(ecm, e))
shared_edges.push_back(e);
std::size_t nb_shared_edges = shared_edges.size();
std::vector<halfedge_descriptor> hedges_to_update;
typedef CGAL::dynamic_halfedge_property_t<bool> H_tag;
typename boost::property_map<TriangleMesh, H_tag>::type
no_target_update = get(H_tag(), tm);
std::vector< std::pair<halfedge_descriptor, vertex_descriptor> > vertices_to_duplicate;
//collect border halfedges having as target one of the edge endpoints
std::set<halfedge_descriptor> extra_border_hedges;
for(std::size_t k=0; k<nb_shared_edges; ++k)
{
for(halfedge_descriptor h : halfedges_around_target(target(shared_edges[k], tm), tm))
if(is_border(h, tm))
extra_border_hedges.insert(h);
for(halfedge_descriptor h : halfedges_around_target(source(shared_edges[k], tm), tm))
if(is_border(h, tm))
extra_border_hedges.insert(h);
}
for(halfedge_descriptor h : extra_border_hedges)
{
put(no_target_update, h, true);
set_halfedge(target(h, tm), h, tm);
hedges_to_update.push_back(h);
}
// now duplicate the edge and set its pointers
for(std::size_t k=0; k<nb_shared_edges; ++k)
{
halfedge_descriptor h = halfedge(shared_edges[k], tm);
face_descriptor fh = face(h, tm);
//add edge
halfedge_descriptor new_hedge = halfedge(add_edge(tm), tm),
new_opp = opposite(new_hedge,tm);
vertex_descriptor vt = target(h, tm);
vertex_descriptor vs = source(h, tm);
//replace h with new_hedge
set_next(new_hedge, next(h, tm), tm);
set_next(prev(h, tm), new_hedge, tm);
set_face(new_hedge, fh, tm);
set_halfedge(fh, new_hedge, tm);
set_target(new_hedge, vt, tm);
set_target(new_opp, vs, tm);
set_face(new_opp, GT::null_face(), tm);
set_face(h, GT::null_face(), tm);
// handle vertices to duplicate
halfedge_descriptor h_vt = halfedge(vt, tm);
if(get(no_target_update, h_vt))
vertices_to_duplicate.push_back(std::make_pair(h, vt));
else
set_halfedge(vt, h, tm);
halfedge_descriptor h_vs = halfedge(vs, tm);
if(get(no_target_update, h_vs))
vertices_to_duplicate.push_back(std::make_pair(new_opp, vs));
else
set_halfedge(vs, new_opp, tm);
hedges_to_update.push_back(h);
put(no_target_update, h, true);
hedges_to_update.push_back(new_opp);
put(no_target_update, new_opp, true);
CGAL_assertion(next(prev(new_hedge, tm), tm) == new_hedge);
CGAL_assertion(prev(next(new_hedge, tm), tm) == new_hedge);
}
// update next/prev relationship
for(halfedge_descriptor h : hedges_to_update)
{
CGAL_assertion(is_border(h, tm));
halfedge_descriptor h_opp = opposite(h, tm);
// set next pointer of h, visiting faces inside the patch we consider
halfedge_descriptor candidate = opposite(prev(h_opp, tm), tm);
while (!is_border(candidate, tm))
candidate = opposite(prev(candidate, tm), tm);
set_next(h, candidate, tm);
CGAL_assertion(prev(next(h_opp, tm), tm)==h_opp);
CGAL_assertion(prev(next(h, tm), tm) == h);
CGAL_assertion(is_border(next(h, tm), tm));
}
for(const std::pair<halfedge_descriptor, vertex_descriptor>& p : vertices_to_duplicate)
{
vertex_descriptor nv = add_vertex(tm);
put(vpm, nv, get(vpm, p.second));
for(halfedge_descriptor h : halfedges_around_target(p.first, tm))
set_target(h, nv, tm);
set_halfedge(nv, p.first, tm);
}
CGAL_assertion_code(for(halfedge_descriptor h : hedges_to_update))
{
CGAL_assertion(next(prev(h, tm), tm) == h);
CGAL_assertion(prev(next(h, tm), tm) == h);
}
for(halfedge_descriptor h : hedges_to_update)
{
for(halfedge_descriptor hh : halfedges_around_target(h, tm))
if(h!=hh)
set_target(hh, target(h, tm), tm);
}
CGAL_assertion(is_valid_polygon_mesh(tm));
}
} // end of internal namespace
/**
@@ -269,14 +410,11 @@ clip_to_bbox(const Plane_3& plane,
* See Subsection \ref coref_clip for more details.
* \attention With the current implementation, `clipper` will be modified (refined with the intersection with `tm`).
*
* \pre \link CGAL::Polygon_mesh_processing::does_self_intersect() `!CGAL::Polygon_mesh_processing::does_self_intersect(tm1)` \endlink
* \pre \link CGAL::Polygon_mesh_processing::does_self_intersect() `!CGAL::Polygon_mesh_processing::does_self_intersect(tm)` \endlink
* \pre \link CGAL::Polygon_mesh_processing::does_self_intersect() `!CGAL::Polygon_mesh_processing::does_self_intersect(clipper)` \endlink
* \pre \link CGAL::Polygon_mesh_processing::does_bound_a_volume() `CGAL::Polygon_mesh_processing::does_bound_a_volume(clipper)` \endlink
*
* @tparam TriangleMesh a model of `MutableFaceGraph`, `HalfedgeListGraph` and `FaceListGraph`.
* If `TriangleMesh` has an internal property map for `CGAL::face_index_t`,
* as a named parameter, then it must be initialized.
*
* @tparam NamedParameters1 a sequence of \ref pmp_namedparameters "Named Parameters"
* @tparam NamedParameters2 a sequence of \ref pmp_namedparameters "Named Parameters"
*
@@ -291,9 +429,8 @@ clip_to_bbox(const Plane_3& plane,
* If this parameter is omitted, an internal property map for
* `CGAL::vertex_point_t` must be available in `TriangleMesh`
* \cgalParamEnd
* \cgalParamBegin{face_index_map} a property map containing the index of each face of `tm` (`clipper`).
* Note that if the property map is writable, the indices of the faces
* of `tm` and `clipper` will be set after refining `tm` with the intersection with `clipper`.
* \cgalParamBegin{face_index_map}
* a property map containing a unique index for each face of `tm` (`clipper`).
* \cgalParamEnd
* \cgalParamBegin{visitor} a class model of `PMPCorefinementVisitor`
* that is used to track the creation of new faces.
@@ -301,14 +438,14 @@ clip_to_bbox(const Plane_3& plane,
* \cgalParamBegin{throw_on_self_intersection} if `true`,
* the set of triangles closed to the intersection of `tm` and `clipper` will be
* checked for self-intersections and `CGAL::Polygon_mesh_processing::Corefinement::Self_intersection_exception`
* will be thrown if at least one is found.
* will be thrown if at least one is found. Default value is `false`.
* \cgalParamEnd
* \cgalParamBegin{clip_volume} if `true` and `tm` is closed, the clipping will be done on
* the volume \link coref_def_subsec bounded \endlink by `tm` rather than on its surface
* (i.e., `tm` will be kept closed).
* (i.e., `tm` will be kept closed). Default value is `false`.
* \cgalParamEnd
* \cgalParamBegin{use_compact_clipper} if `false`, the parts of `tm` coplanar with `clipper` will not be part of the output.
* This option has an effect only if a surface and not a volume is clipped (i.e. if `clip_volume` is `false` or if `tm` is open).
* \cgalParamBegin{use_compact_clipper} if `false`, the parts of `tm` coplanar with `clipper`
* will not be part of the output. Default value is `true`.
* \cgalParamEnd
* \cgalNamedParamsEnd
*
@@ -319,41 +456,21 @@ template <class TriangleMesh,
class NamedParameters1,
class NamedParameters2>
bool
clip( TriangleMesh& tm,
TriangleMesh& clipper,
clip(TriangleMesh& tm,
TriangleMesh& clipper,
const NamedParameters1& np_tm,
const NamedParameters2& np_c)
{
const bool clip_volume =
parameters::choose_parameter(parameters::get_parameter(np_tm, internal_np::clip_volume), false);
if (clip_volume && is_closed(tm))
if(clip_volume && is_closed(tm))
return corefine_and_compute_intersection(tm, clipper, tm, np_tm, np_c);
return corefine_and_compute_intersection(tm, clipper, tm,
np_tm.use_bool_op_to_clip_surface(true),
np_c);
}
namespace internal{
template <class TriangleMesh, class NamedParameters>
bool dispatch_clip_call(TriangleMesh& tm, TriangleMesh& clipper,
const NamedParameters& np, Tag_false)
{
return clip(tm, clipper,
np.face_index_map(get(CGAL::dynamic_face_property_t<std::size_t>(), tm)),
parameters::face_index_map(get(CGAL::dynamic_face_property_t<std::size_t>(), clipper)));
}
template <class TriangleMesh, class NamedParameters>
bool dispatch_clip_call(TriangleMesh& tm, TriangleMesh& clipper,
const NamedParameters& np, Tag_true)
{
return clip(tm, clipper,
np.face_index_map(get(face_index, tm)),
parameters::face_index_map(get(face_index, clipper)));
}
}
/**
* \ingroup PMP_corefinement_grp
* clips `tm` by keeping the part that is on the negative side of `plane` (side opposite to its normal vector).
@@ -365,8 +482,6 @@ bool dispatch_clip_call(TriangleMesh& tm, TriangleMesh& clipper,
* \pre \link CGAL::Polygon_mesh_processing::does_self_intersect() `!CGAL::Polygon_mesh_processing::does_self_intersect(tm)` \endlink
*
* @tparam TriangleMesh a model of `MutableFaceGraph`, `HalfedgeListGraph` and `FaceListGraph`.
* If `TriangleMesh` has an internal property map for `CGAL::face_index_t`,
* as a named parameter, then it must be initialized.
* An internal property map for `CGAL::vertex_point_t` must be available.
*
* @tparam NamedParameters a sequence of \ref pmp_namedparameters "Named Parameters"
@@ -383,14 +498,14 @@ bool dispatch_clip_call(TriangleMesh& tm, TriangleMesh& clipper,
* \cgalParamBegin{throw_on_self_intersection} if `true`,
* the set of triangles closed to the intersection of `tm` and `plane` will be
* checked for self-intersections and `CGAL::Polygon_mesh_processing::Corefinement::Self_intersection_exception`
* will be thrown if at least one is found.
* will be thrown if at least one is found. Default value is `false`.
* \cgalParamEnd
* \cgalParamBegin{clip_volume} if `true` and `tm` is closed, the clipping will be done on
* the volume \link coref_def_subsec bounded \endlink by `tm` rather than on its surface
* (i.e., `tm` will be kept closed).
* (i.e., `tm` will be kept closed). Default value is `false`.
* \cgalParamEnd
* \cgalParamBegin{use_compact_clipper} if `false` and `clip_volume` is `false` and `tm` is open, the parts of `tm` coplanar with `plane`
* will not be part of the output.
* \cgalParamBegin{use_compact_clipper} if `false` the parts of `tm` coplanar with `plane`
* will not be part of the output. Default value is `true`.
* \cgalNamedParamsEnd
*
* @return `true` if the output surface mesh is manifold.
@@ -398,15 +513,15 @@ bool dispatch_clip_call(TriangleMesh& tm, TriangleMesh& clipper,
*/
template <class TriangleMesh,
class NamedParameters>
bool clip( TriangleMesh& tm,
#ifdef DOXYGEN_RUNNING
bool clip(TriangleMesh& tm,
#ifdef DOXYGEN_RUNNING
const Plane_3& plane,
#else
#else
const typename GetGeomTraits<TriangleMesh, NamedParameters>::type::Plane_3& plane,
#endif
#endif
const NamedParameters& np)
{
if( boost::begin(faces(tm))==boost::end(faces(tm)) ) return true;
if(boost::begin(faces(tm))==boost::end(faces(tm))) return true;
CGAL::Bbox_3 bbox = ::CGAL::Polygon_mesh_processing::bbox(tm);
@@ -428,9 +543,8 @@ bool clip( TriangleMesh& tm,
default:
break;
}
// dispatch is needed because face index map for tm and clipper have to be of the same time
return internal::dispatch_clip_call(tm, clipper,
np, CGAL::graph_has_property<TriangleMesh, CGAL::face_index_t>());
return clip(tm, clipper, np, parameters::all_default());
}
/**
@@ -444,8 +558,6 @@ bool clip( TriangleMesh& tm,
* \pre \link CGAL::Polygon_mesh_processing::does_self_intersect() `!CGAL::Polygon_mesh_processing::does_self_intersect(tm)` \endlink
*
* @tparam TriangleMesh a model of `MutableFaceGraph`, `HalfedgeListGraph` and `FaceListGraph`.
* If `TriangleMesh` has an internal property map for `CGAL::face_index_t`,
* as a named parameter, then it must be initialized.
* An internal property map for `CGAL::vertex_point_t` must be available.
*
* @tparam NamedParameters a sequence of \ref pmp_namedparameters "Named Parameters"
@@ -476,15 +588,15 @@ bool clip( TriangleMesh& tm,
*/
template <class TriangleMesh,
class NamedParameters>
bool clip( TriangleMesh& tm,
#ifdef DOXYGEN_RUNNING
bool clip(TriangleMesh& tm,
#ifdef DOXYGEN_RUNNING
const Iso_cuboid_3& iso_cuboid,
#else
#else
const typename GetGeomTraits<TriangleMesh, NamedParameters>::type::Iso_cuboid_3& iso_cuboid,
#endif
#endif
const NamedParameters& np)
{
if( boost::begin(faces(tm))==boost::end(faces(tm)) ) return true;
if(boost::begin(faces(tm))==boost::end(faces(tm))) return true;
TriangleMesh clipper;
make_hexahedron(iso_cuboid[0], iso_cuboid[1], iso_cuboid[2], iso_cuboid[3],
@@ -492,16 +604,192 @@ bool clip( TriangleMesh& tm,
clipper);
triangulate_faces(clipper);
// dispatch is needed because face index map for tm and clipper have to be of the same time
return internal::dispatch_clip_call(tm, clipper,
np, CGAL::graph_has_property<TriangleMesh, CGAL::face_index_t>());
return clip(tm, clipper, np, parameters::all_default());
}
/*!
* \ingroup PMP_corefinement_grp
* corefines `tm` and `splitter` and duplicates edges in `tm` that are on the intersection with `splitter`.
*
* \pre \link CGAL::Polygon_mesh_processing::does_self_intersect() `!CGAL::Polygon_mesh_processing::does_self_intersect(tm)` \endlink
* \pre \link CGAL::Polygon_mesh_processing::does_self_intersect() `!CGAL::Polygon_mesh_processing::does_self_intersect(splitter)` \endlink
*
* @tparam TriangleMesh a model of `MutableFaceGraph`, `HalfedgeListGraph` and `FaceListGraph`.
*
* @tparam NamedParameters1 a sequence of \ref pmp_namedparameters "Named Parameters"
* @tparam NamedParameters2 a sequence of \ref pmp_namedparameters "Named Parameters"
*
* @param tm input triangulated surface mesh
* @param splitter triangulated surface mesh used to split `tm`
* @param np_tm optional sequence of \ref pmp_namedparameters "Named Parameters" among the ones listed below
* @param np_s optional sequence of \ref pmp_namedparameters "Named Parameters" among the ones listed below
*
* \cgalNamedParamsBegin
* \cgalParamBegin{vertex_point_map}
* the property map with the points associated to the vertices of `tm` (`splitter`).
* If this parameter is omitted, an internal property map for
* `CGAL::vertex_point_t` must be available in `TriangleMesh`
* \cgalParamEnd
* \cgalParamBegin{visitor} a class model of `PMPCorefinementVisitor`
* that is used to track the creation of new faces.
* \cgalParamEnd
* \cgalParamBegin{throw_on_self_intersection} if `true`,
* the set of triangles closed to the intersection of `tm` and `splitter` will be
* checked for self-intersections and `CGAL::Polygon_mesh_processing::Corefinement::Self_intersection_exception`
* will be thrown if at least one is found.
* \cgalParamEnd
* \cgalNamedParamsEnd
*/
template <class TriangleMesh,
class NamedParameters1,
class NamedParameters2>
void split(TriangleMesh& tm,
TriangleMesh& splitter,
const NamedParameters1& np_tm,
const NamedParameters2& np_s)
{
namespace PMP = CGAL::Polygon_mesh_processing;
using parameters::get_parameter;
using parameters::choose_parameter;
typedef typename GetVertexPointMap<TriangleMesh, NamedParameters1>::type VPM1;
typedef typename GetVertexPointMap<TriangleMesh, NamedParameters2>::type VPM2;
typedef typename boost::template property_map<TriangleMesh, CGAL::dynamic_edge_property_t<bool> >::type Ecm;
VPM1 vpm_tm = choose_parameter(get_parameter(np_tm, internal_np::vertex_point),
get_property_map(vertex_point, tm));
VPM2 vpm_s = choose_parameter(get_parameter(np_s, internal_np::vertex_point),
get_property_map(vertex_point, splitter));
Ecm ecm = get(CGAL::dynamic_edge_property_t<bool>(), tm);
// create a constrained edge map and corefine input mesh with the splitter,
// and mark edges
PMP::corefine(tm, splitter,
CGAL::parameters::vertex_point_map(vpm_tm).edge_is_constrained_map(ecm),
CGAL::parameters::vertex_point_map(vpm_s));
//split mesh along marked edges
internal::split_along_edges(tm, ecm, vpm_tm);
}
/**
* \ingroup PMP_corefinement_grp
* adds intersection edges of `plane` and `tm` in `tm` and duplicates those edges.
*
* \note In the current implementation it is not possible to set the vertex point map and the default will be used.
* \pre \link CGAL::Polygon_mesh_processing::does_self_intersect() `!CGAL::Polygon_mesh_processing::does_self_intersect(tm)` \endlink
*
* @tparam TriangleMesh a model of `MutableFaceGraph`, `HalfedgeListGraph` and `FaceListGraph`
* An internal property map for `CGAL::vertex_point_t` must be available.
*
* @tparam NamedParameters a sequence of \ref pmp_namedparameters "Named Parameters"
*
* @param tm input triangulated surface mesh
* @param plane the plane that will be used to split `tm`.
* `Plane_3` is the plane type for the same CGAL kernel as the point of the vertex point map of `tm`.
* @param np optional sequence of \ref pmp_namedparameters "Named Parameters" among the ones listed below
*
* \cgalNamedParamsBegin
* \cgalParamBegin{visitor} a class model of `PMPCorefinementVisitor`
* that is used to track the creation of new faces.
* \cgalParamEnd
* \cgalParamBegin{throw_on_self_intersection} if `true`,
* the set of triangles closed to the intersection of `tm` and `plane` will be
* checked for self-intersections and `CGAL::Polygon_mesh_processing::Corefinement::Self_intersection_exception`
* will be thrown if at least one is found.
* \cgalParamEnd
* \cgalNamedParamsEnd
*/
template <class TriangleMesh,
class NamedParameters>
void split(TriangleMesh& tm,
#ifdef DOXYGEN_RUNNING
const Plane_3& plane,
#else
const typename GetGeomTraits<TriangleMesh, NamedParameters>::type::Plane_3& plane,
#endif
const NamedParameters& np)
{
using parameters::get_parameter;
using parameters::choose_parameter;
namespace PMP = CGAL::Polygon_mesh_processing;
namespace params = PMP::parameters;
// create a splitter mesh for the splitting plane using an internal CGAL function
CGAL::Bbox_3 bbox = ::CGAL::Polygon_mesh_processing::bbox(tm, np);
double xd = (std::max)(1., 0.01 * (bbox.xmax() - bbox.xmin()));
double yd = (std::max)(1., 0.01 * (bbox.ymax() - bbox.ymin()));
double zd = (std::max)(1., 0.01 * (bbox.zmax() - bbox.zmin()));
bbox = CGAL::Bbox_3(bbox.xmin()-xd, bbox.ymin()-yd, bbox.zmin()-zd,
bbox.xmax()+xd, bbox.ymax()+yd, bbox.zmax()+zd);
TriangleMesh splitter;
CGAL::Oriented_side os = PMP::internal::clip_to_bbox(plane, bbox, splitter, PMP::parameters::all_default());
if(os == CGAL::ON_ORIENTED_BOUNDARY)
return split(tm, splitter, np, params::all_default());
//else nothing to do, no intersection.
}
/**
* \ingroup PMP_corefinement_grp
* adds intersection edges of `iso_cuboid` and `tm` in `tm` and duplicates those edges.
*
* \note In the current implementation it is not possible to set the vertex point map and the default will be used.
* \note `Iso_cuboid_3` must be from the same %Kernel as the point of the vertex point map.
* \pre \link CGAL::Polygon_mesh_processing::does_self_intersect() `!CGAL::Polygon_mesh_processing::does_self_intersect(tm)` \endlink
*
* @tparam TriangleMesh a model of `MutableFaceGraph`, `HalfedgeListGraph` and `FaceListGraph`
* An internal property map for `CGAL::vertex_point_t` must be available.
*
* @tparam NamedParameters a sequence of \ref pmp_namedparameters "Named Parameters"
*
* @param tm input triangulated surface mesh
* @param iso_cuboid iso-cuboid used to split `tm`.
* @param np optional sequence of \ref pmp_namedparameters "Named Parameters" among the ones listed below
*
* \cgalNamedParamsBegin
* \cgalParamBegin{visitor} a class model of `PMPCorefinementVisitor`
* that is used to track the creation of new faces.
* \cgalParamEnd
* \cgalParamBegin{throw_on_self_intersection} if `true`,
* the set of triangles closed to the intersection of `tm` and `iso_cuboid` will be
* checked for self-intersections and `CGAL::Polygon_mesh_processing::Corefinement::Self_intersection_exception`
* will be thrown if at least one is found.
* \cgalParamEnd
* \cgalNamedParamsEnd
*/
template <class TriangleMesh,
class NamedParameters>
void split(TriangleMesh& tm,
#ifdef DOXYGEN_RUNNING
const Iso_cuboid_3& iso_cuboid,
#else
const typename GetGeomTraits<TriangleMesh, NamedParameters>::type::Iso_cuboid_3& iso_cuboid,
#endif
const NamedParameters& np)
{
TriangleMesh splitter;
make_hexahedron(iso_cuboid[0], iso_cuboid[1], iso_cuboid[2], iso_cuboid[3],
iso_cuboid[4], iso_cuboid[5], iso_cuboid[6], iso_cuboid[7],
splitter);
triangulate_faces(splitter);
split(tm, splitter, np, parameters::all_default());
}
/// \cond SKIP_IN_MANUAL
// convenience overloads
template <class TriangleMesh>
bool clip( TriangleMesh& tm,
bool clip(TriangleMesh& tm,
const typename GetGeomTraits<TriangleMesh>::type::Plane_3& plane)
{
return clip(tm, plane, parameters::all_default());
@@ -509,7 +797,7 @@ bool clip( TriangleMesh& tm,
// convenience overloads
template <class TriangleMesh>
bool clip( TriangleMesh& tm,
bool clip(TriangleMesh& tm,
const typename GetGeomTraits<TriangleMesh>::type::Iso_cuboid_3& iso_cuboid)
{
return clip(tm, iso_cuboid, parameters::all_default());
@@ -519,8 +807,8 @@ bool clip( TriangleMesh& tm,
template <class TriangleMesh,
class NamedParameters1>
bool
clip( TriangleMesh& tm,
TriangleMesh& clipper,
clip(TriangleMesh& tm,
TriangleMesh& clipper,
const NamedParameters1& np_tm)
{
return clip(tm, clipper, np_tm, parameters::all_default());
@@ -529,11 +817,47 @@ clip( TriangleMesh& tm,
// convenience overload
template <class TriangleMesh>
bool
clip( TriangleMesh& tm,
TriangleMesh& clipper)
clip(TriangleMesh& tm,
TriangleMesh& clipper)
{
return clip(tm, clipper, parameters::all_default());
}
// convenience overload
template <class TriangleMesh,
class NamedParameters1>
void
split(TriangleMesh& tm,
TriangleMesh& splitter,
const NamedParameters1& np_tm)
{
split(tm, splitter, np_tm, parameters::all_default());
}
// convenience overload
template <class TriangleMesh>
void
split(TriangleMesh& tm,
TriangleMesh& splitter)
{
split(tm, splitter, parameters::all_default());
}
template <class TriangleMesh>
void split(TriangleMesh& tm,
const typename GetGeomTraits<TriangleMesh>::type::Plane_3& plane)
{
split(tm, plane, parameters::all_default());
}
template <class TriangleMesh>
void split(TriangleMesh& tm,
const typename GetGeomTraits<TriangleMesh>::type::Iso_cuboid_3& iso_cuboid)
{
split(tm, iso_cuboid, parameters::all_default());
}
/// \endcond
} } //end of namespace CGAL::Polygon_mesh_processing
@@ -98,6 +98,12 @@ void sum_normals(const PM& pmesh,
const Point_ref pvnn = get(vpmap, target(next(he, pmesh), pmesh));
const Vector n = internal::triangle_normal(pv, pvn, pvnn, traits);
#ifdef CGAL_PMP_COMPUTE_NORMAL_DEBUG_PP
std::cout << "Normal of " << f << " pts: " << pv << " ; " << pvn << " ; " << pvnn << std::endl;
std::cout << " --> " << n << std::endl;
#endif
sum = traits.construct_sum_of_vectors_3_object()(sum, n);
he = next(he, pmesh);
@@ -146,11 +152,11 @@ compute_face_normal(typename boost::graph_traits<PolygonMesh>::face_descriptor f
using parameters::get_parameter;
typedef typename GetGeomTraits<PolygonMesh, NamedParameters>::type GT;
GT traits = choose_parameter(get_parameter(np, internal_np::geom_traits), GT());
GT traits = choose_parameter<GT>(get_parameter(np, internal_np::geom_traits));
typedef typename GetVertexPointMap<PolygonMesh, NamedParameters>::const_type VPMap;
VPMap vpmap = choose_parameter(get_parameter(np, internal_np::vertex_point),
get_const_property_map(vertex_point, pmesh));
get_const_property_map(vertex_point, pmesh));
typedef typename GT::Point_3 Point;
typedef typename GT::Vector_3 Vector_3;
@@ -206,7 +212,7 @@ void compute_face_normals(const PolygonMesh& pmesh,
{
typename Kernel::Vector_3 vec = compute_face_normal(f, pmesh, np);
put(face_normals, f, vec);
#ifdef CGAL_PMP_COMPUTE_NORMAL_DEBUG
#ifdef CGAL_PMP_COMPUTE_NORMAL_DEBUG_PP
std::cout << "normal at face " << f << " is " << get(face_normals, f) << std::endl;
#endif
}
@@ -522,6 +528,10 @@ compute_vertex_normal_as_sum_of_weighted_normals(typename boost::graph_traits<Po
typename GT::Construct_vector_3 cv_3 = traits.construct_vector_3_object();
typename GT::Compute_squared_length_3 csl_3 = traits.compute_squared_length_3_object();
#ifdef CGAL_PMP_COMPUTE_NORMAL_DEBUG_PP
std::cout << "Compute normal as weighted sums; type: " << vn_type << std::endl;
#endif
Vector_3 normal = cv_3(CGAL::NULL_VECTOR);
halfedge_descriptor h = halfedge(v, pmesh);
@@ -548,7 +558,13 @@ compute_vertex_normal_as_sum_of_weighted_normals(typename boost::graph_traits<Po
const FT den = CGAL::approximate_sqrt(csl_3(v1) * csl_3(v2));
if(den == FT(0))
{
#ifdef CGAL_PMP_COMPUTE_NORMAL_DEBUG_PP
std::cout << "Null denominator, switching to no weights" << std::endl;
#endif
return compute_vertex_normal_as_sum_of_weighted_normals(v, NO_WEIGHT, face_normals, vpmap, pmesh, traits);
}
n = traits.construct_scaled_vector_3_object()(n, FT(1) / den);
normal = traits.construct_sum_of_vectors_3_object()(normal, n);
@@ -614,11 +630,11 @@ compute_vertex_normal(typename boost::graph_traits<PolygonMesh>::vertex_descript
typedef typename GetGeomTraits<PolygonMesh, NamedParameters>::type GT;
typedef typename GT::Vector_3 Vector_3;
GT traits = choose_parameter(get_parameter(np, internal_np::geom_traits), GT());
GT traits = choose_parameter<GT>(get_parameter(np, internal_np::geom_traits));
typedef typename GetVertexPointMap<PolygonMesh, NamedParameters>::const_type VPMap;
VPMap vpmap = choose_parameter(get_parameter(np, internal_np::vertex_point),
get_const_property_map(vertex_point, pmesh));
get_const_property_map(vertex_point, pmesh));
typedef std::map<face_descriptor, Vector_3> Face_vector_map;
typedef boost::associative_property_map<Face_vector_map> Default_map;
@@ -632,9 +648,7 @@ compute_vertex_normal(typename boost::graph_traits<PolygonMesh>::vertex_descript
const bool must_compute_face_normals = is_default_parameter(get_parameter(np, internal_np::face_normal));
#ifdef CGAL_PMP_COMPUTE_NORMAL_DEBUG_PP
std::cout << std::endl << std::endl;
std::cout << "----------------------------------------------------------------------" << std::endl;
std::cout << "compute vertex at " << get(vpmap, v)
std::cout << "<----- compute vertex normal at " << get(vpmap, v)
<< ", must compute face normals? " << must_compute_face_normals << std::endl;
#endif
@@ -654,7 +668,7 @@ compute_vertex_normal(typename boost::graph_traits<PolygonMesh>::vertex_descript
}
}
#ifdef CGAL_PMP_COMPUTE_NORMAL_DEBUG
#ifdef CGAL_PMP_COMPUTE_NORMAL_DEBUG_PP
std::cout << "Incident face normals:" << std::endl;
for(halfedge_descriptor h : CGAL::halfedges_around_target(v, pmesh))
{
@@ -729,7 +743,7 @@ void compute_vertex_normals(const PolygonMesh& pmesh,
typedef typename boost::property_map<PolygonMesh, Face_normal_tag>::const_type Face_normal_dmap;
#ifdef CGAL_PMP_COMPUTE_NORMAL_DEBUG_PP
GT traits = choose_parameter(get_parameter(np, internal_np::geom_traits), GT());
GT traits = choose_parameter<GT>(get_parameter(np, internal_np::geom_traits));
typedef typename GetVertexPointMap<PolygonMesh, NamedParameters>::const_type VPMap;
VPMap vpmap = choose_parameter(get_parameter(np, internal_np::vertex_point),
@@ -28,6 +28,12 @@
#include <CGAL/assertions.h>
#include <CGAL/boost/graph/iterator.h>
#include <CGAL/boost/graph/helpers.h>
#include <CGAL/boost/graph/Face_filtered_graph.h>
#include <CGAL/boost/graph/copy_face_graph.h>
#include <CGAL/assertions.h>
#include <CGAL/tuple.h>
#include <CGAL/boost/graph/Dual.h>
#include <CGAL/Default.h>
#include <CGAL/Dynamic_property_map.h>
@@ -55,7 +61,7 @@ namespace internal {
}
};
// A property map
// A property map
template <typename G>
struct No_constraint {
friend bool get(No_constraint<G>, typename boost::graph_traits<G>::edge_descriptor)
@@ -130,8 +136,7 @@ connected_component(typename boost::graph_traits<PolygonMesh>::face_descriptor s
internal::No_constraint<PolygonMesh>//default
> ::type EdgeConstraintMap;
EdgeConstraintMap ecmap
= choose_parameter(get_parameter(np, internal_np::edge_is_constrained),
internal::No_constraint<PolygonMesh>());
= choose_parameter<EdgeConstraintMap>(get_parameter(np, internal_np::edge_is_constrained));
typedef typename boost::graph_traits<PolygonMesh>::face_descriptor face_descriptor;
typedef typename boost::graph_traits<PolygonMesh>::halfedge_descriptor halfedge_descriptor;
@@ -145,7 +150,7 @@ connected_component(typename boost::graph_traits<PolygonMesh>::face_descriptor s
if (!already_processed.insert(seed_face).second) continue;
*out++=seed_face;
for(halfedge_descriptor hd :
CGAL::halfedges_around_face(halfedge(seed_face, pmesh), pmesh) )
halfedges_around_face(halfedge(seed_face, pmesh), pmesh) )
{
if(! get(ecmap, edge(hd, pmesh))){
face_descriptor neighbor = face( opposite(hd, pmesh), pmesh );
@@ -171,9 +176,6 @@ connected_component(typename boost::graph_traits<PolygonMesh>::face_descriptor s
* \ingroup keep_connected_components_grp
* computes for each face the index of the corresponding connected component.
*
* A property map for `CGAL::face_index_t` must be either available as an internal property map
* to `pmesh` or provided as one of the \ref pmp_namedparameters "Named Parameters".
*
* \tparam PolygonMesh a model of `FaceListGraph`
* \tparam FaceComponentMap a model of `WritablePropertyMap` with
`boost::graph_traits<PolygonMesh>::%face_descriptor` as key type and
@@ -213,13 +215,12 @@ connected_components(const PolygonMesh& pmesh,
NamedParameters,
internal::No_constraint<PolygonMesh>//default
> ::type EdgeConstraintMap;
EdgeConstraintMap ecmap
= choose_parameter(get_parameter(np, internal_np::edge_is_constrained),
internal::No_constraint<PolygonMesh>());
typedef typename GetFaceIndexMap<PolygonMesh, NamedParameters>::const_type FaceIndexMap;
FaceIndexMap fimap = choose_parameter(get_parameter(np, internal_np::face_index),
get_const_property_map(boost::face_index, pmesh));
EdgeConstraintMap ecmap
= choose_parameter<EdgeConstraintMap>(get_parameter(np, internal_np::edge_is_constrained));
typedef typename GetInitializedFaceIndexMap<PolygonMesh, NamedParameters>::const_type FaceIndexMap;
FaceIndexMap fimap = get_initialized_face_index_map(pmesh, np);
typename boost::property_traits<FaceComponentMap>::value_type i=0;
std::vector<bool> handled(num_faces(pmesh), false);
@@ -298,7 +299,8 @@ template <typename PolygonMesh,
std::size_t number_of_connected_components(const PolygonMesh& pmesh,
const CGAL_PMP_NP_CLASS& np)
{
typedef CGAL::dynamic_face_property_t<std::size_t> Face_property_tag;
typedef typename boost::graph_traits<PolygonMesh>::faces_size_type faces_size_type;
typedef CGAL::dynamic_face_property_t<faces_size_type> Face_property_tag;
typedef typename boost::property_map<PolygonMesh, Face_property_tag >::const_type Patch_ids_map;
Patch_ids_map patch_ids_map = get(Face_property_tag(), pmesh);
@@ -323,10 +325,6 @@ std::size_t number_of_connected_components(const PolygonMesh& pmesh)
* By default, the size of a face is `1` (and thus the size of a connected component is the number
* of faces it contains), but it is also possible to pass custom sizes, such as the area of the face.
*
* Property maps for `CGAL::face_index_t` and `CGAL::vertex_index_t`
* must be either available as internal property maps
* to `pmesh` or provided as \ref pmp_namedparameters "Named Parameters".
*
* \tparam PolygonMesh a model of `FaceListGraph` and `MutableFaceGraph`
* \tparam NamedParameters a sequence of \ref pmp_namedparameters "Named Parameters"
*
@@ -368,10 +366,8 @@ std::size_t keep_largest_connected_components(PolygonMesh& pmesh,
using parameters::choose_parameter;
using parameters::get_parameter;
// FaceIndexMap
typedef typename GetFaceIndexMap<PM, NamedParameters>::type FaceIndexMap;
FaceIndexMap fimap = choose_parameter(get_parameter(np, internal_np::face_index),
get_property_map(boost::face_index, pmesh));
typedef typename CGAL::GetInitializedFaceIndexMap<PolygonMesh, NamedParameters>::type FaceIndexMap;
FaceIndexMap fimap = CGAL::get_initialized_face_index_map(pmesh, np);
// FaceSizeMap
typedef typename internal_np::Lookup_named_param_def<internal_np::face_size_map_t,
@@ -388,7 +384,7 @@ std::size_t keep_largest_connected_components(PolygonMesh& pmesh,
typedef typename internal_np::Lookup_named_param_def<internal_np::output_iterator_t,
NamedParameters,
Emptyset_iterator>::type Output_iterator;
Output_iterator out = choose_parameter(get_parameter(np, internal_np::output_iterator), Emptyset_iterator());
Output_iterator out = choose_parameter<Output_iterator>(get_parameter(np, internal_np::output_iterator));
// vector_property_map
boost::vector_property_map<std::size_t, FaceIndexMap> face_cc(fimap);
@@ -456,10 +452,6 @@ std::size_t keep_largest_connected_components(PolygonMesh& pmesh,
* the size of a connected component is the number of faces it contains), but it is also possible
* to pass custom sizes, such as the area of the face.
*
* Property maps for `CGAL::face_index_t` and `CGAL::vertex_index_t`
* must be either available as internal property maps
* to `pmesh` or provided as \ref pmp_namedparameters "Named Parameters".
*
* \tparam PolygonMesh a model of `FaceListGraph` and `MutableFaceGraph`
* \tparam ThresholdValueType the type of the threshold value
* \tparam NamedParameters a sequence of \ref pmp_namedparameters "Named Parameters"
@@ -505,15 +497,13 @@ std::size_t keep_large_connected_components(PolygonMesh& pmesh,
using parameters::choose_parameter;
using parameters::get_parameter;
// FaceIndexMap
typedef typename GetFaceIndexMap<PM, NamedParameters>::type FaceIndexMap;
FaceIndexMap fim = choose_parameter(get_parameter(np, internal_np::face_index),
get_property_map(boost::face_index, pmesh));
typedef typename CGAL::GetInitializedFaceIndexMap<PolygonMesh, NamedParameters>::type FaceIndexMap;
FaceIndexMap fim = CGAL::get_initialized_face_index_map(pmesh, np);
typedef typename internal_np::Lookup_named_param_def<internal_np::face_size_map_t,
NamedParameters,
Constant_property_map<face_descriptor, std::size_t> // default
>::type FaceSizeMap;
NamedParameters,
Constant_property_map<face_descriptor, std::size_t> // default
>::type FaceSizeMap;
typedef typename boost::property_traits<FaceSizeMap>::value_type Face_size;
CGAL_static_assertion((std::is_convertible<ThresholdValueType, Face_size>::value));
@@ -523,9 +513,9 @@ std::size_t keep_large_connected_components(PolygonMesh& pmesh,
Emptyset_iterator>::type Output_iterator;
FaceSizeMap face_size_pmap = choose_parameter(get_parameter(np, internal_np::face_size_map),
Constant_property_map<face_descriptor, std::size_t>(1));
Constant_property_map<face_descriptor, std::size_t>(1));
const bool dry_run = choose_parameter(get_parameter(np, internal_np::dry_run), false);
Output_iterator out = choose_parameter(get_parameter(np, internal_np::output_iterator), Emptyset_iterator());
Output_iterator out = choose_parameter<Output_iterator>(get_parameter(np, internal_np::output_iterator));
// vector_property_map
boost::vector_property_map<std::size_t, FaceIndexMap> face_cc(fim);
@@ -587,7 +577,6 @@ void keep_or_remove_connected_components(PolygonMesh& pmesh
, bool keep
, const NamedParameters& np)
{
typedef PolygonMesh PM;
using parameters::choose_parameter;
using parameters::get_parameter;
@@ -599,10 +588,8 @@ void keep_or_remove_connected_components(PolygonMesh& pmesh
typedef typename boost::graph_traits<PolygonMesh>::edge_descriptor edge_descriptor;
typedef typename boost::graph_traits<PolygonMesh>::edge_iterator edge_iterator;
//VertexIndexMap
typedef typename GetVertexIndexMap<PM, NamedParameters>::type VertexIndexMap;
VertexIndexMap vim = choose_parameter(get_parameter(np, internal_np::vertex_index),
get_const_property_map(boost::vertex_index, pmesh));
typedef typename GetInitializedVertexIndexMap<PolygonMesh, NamedParameters>::type VertexIndexMap;
VertexIndexMap vim = get_initialized_vertex_index_map(pmesh, np);
std::set<std::size_t> cc_to_keep;
for(std::size_t i : components_to_keep)
@@ -721,10 +708,6 @@ void keep_or_remove_connected_components(PolygonMesh& pmesh
* \note If the removal of the connected components makes `pmesh` a non-manifold surface,
* then the behavior of this function is undefined.
*
* Property maps for `CGAL::vertex_index_t`
* must be either available as internal property map
* to `pmesh` or provided as \ref pmp_namedparameters "Named Parameters".
*
* \tparam PolygonMesh a model of `FaceListGraph` and `MutableFaceGraph`
* \tparam NamedParameters a sequence of \ref pmp_namedparameters "Named Parameters"
* \tparam ComponentRange a range of ids convertible to `std::size`
@@ -764,11 +747,6 @@ void keep_connected_components(PolygonMesh& pmesh
* \note If the removal of the connected components makes `pmesh` a non-manifold surface,
* then the behavior of this function is undefined.
*
* Property maps for `CGAL::vertex_index_t`
* must be either available as internal property map
* to `pmesh` or provided as \ref pmp_namedparameters "Named Parameters".
*
*
* \tparam PolygonMesh a model of `FaceListGraph` and `MutableFaceGraph`
* \tparam NamedParameters a sequence of \ref pmp_namedparameters "Named Parameters"
* \tparam ComponentRange a range of ids convertible to `std::size`
@@ -805,10 +783,6 @@ void remove_connected_components(PolygonMesh& pmesh
* keeps the connected components not designated by the faces in `components_to_remove`,
* and removes the other connected components and all isolated vertices.
*
* Property maps for `CGAL::face_index_t` and `CGAL::vertex_index_t`
* must be either available as internal property maps
* to `pmesh` or provided as \ref pmp_namedparameters "Named Parameters".
*
* \note If the removal of the connected components makes `pmesh` a non-manifold surface,
* then the behavior of this function is undefined.
*
@@ -835,20 +809,19 @@ void remove_connected_components(PolygonMesh& pmesh
, const FaceRange& components_to_remove
, const CGAL_PMP_NP_CLASS& np)
{
if (components_to_remove.empty()) return;
typedef PolygonMesh PM;
typedef typename boost::graph_traits<PM>::face_descriptor face_descriptor;
using parameters::choose_parameter;
using parameters::get_parameter;
//FaceIndexMap
typedef typename GetFaceIndexMap<PM, CGAL_PMP_NP_CLASS>::type FaceIndexMap;
FaceIndexMap fim = choose_parameter(get_parameter(np, internal_np::face_index),
get_property_map(boost::face_index, pmesh));
if (components_to_remove.empty())
return;
typedef PolygonMesh PM;
typedef typename boost::graph_traits<PM>::face_descriptor face_descriptor;
typedef typename CGAL::GetInitializedFaceIndexMap<PolygonMesh, CGAL_PMP_NP_CLASS>::type FaceIndexMap;
FaceIndexMap fim = CGAL::get_initialized_face_index_map(pmesh, np);
//vector_property_map
boost::vector_property_map<std::size_t, FaceIndexMap> face_cc(fim);
connected_components(pmesh, face_cc, np);
std::vector<std::size_t> cc_to_remove;
@@ -863,10 +836,6 @@ void remove_connected_components(PolygonMesh& pmesh
* keeps the connected components designated by the faces in `components_to_keep`,
* and removes the other connected components and all isolated vertices.
*
* Property maps for `CGAL::face_index_t` and `CGAL::vertex_index_t`
* must be either available as internal property maps
* to `pmesh` or provided as \ref pmp_namedparameters "Named Parameters".
*
* \note If the removal of the connected components makes `pmesh` a non-manifold surface,
* then the behavior of this function is undefined.
*
@@ -899,14 +868,10 @@ void keep_connected_components(PolygonMesh& pmesh
using parameters::choose_parameter;
using parameters::get_parameter;
//FaceIndexMap
typedef typename GetFaceIndexMap<PM, CGAL_PMP_NP_CLASS>::type FaceIndexMap;
FaceIndexMap fim = choose_parameter(get_parameter(np, internal_np::face_index),
get_property_map(boost::face_index, pmesh));
typedef typename CGAL::GetInitializedFaceIndexMap<PolygonMesh, CGAL_PMP_NP_CLASS>::type FaceIndexMap;
FaceIndexMap fim = CGAL::get_initialized_face_index_map(pmesh, np);
//vector_property_map
boost::vector_property_map<std::size_t, FaceIndexMap> face_cc(fim);
connected_components(pmesh, face_cc, np);
std::vector<std::size_t> cc_to_keep;
@@ -956,8 +921,137 @@ void keep_connected_components(PolygonMesh& pmesh
CGAL::Polygon_mesh_processing::parameters::all_default());
}
} // namespace Polygon_mesh_processing
namespace internal {
template <typename G>
struct No_mark
{
friend bool get(No_mark<G>, typename boost::graph_traits<G>::edge_descriptor) { return false; }
friend void put(No_mark<G>, typename boost::graph_traits<G>::edge_descriptor, bool) { }
};
template < class PolygonMesh, class PolygonMeshRange,
class FIMap, class VIMap,
class HIMap, class Ecm, class NamedParameters >
void split_connected_components_impl(FIMap fim,
HIMap him,
VIMap vim,
Ecm ecm,
PolygonMeshRange& range,
const PolygonMesh& tm,
const NamedParameters& np)
{
typedef typename boost::graph_traits<PolygonMesh>::faces_size_type faces_size_type;
typedef typename internal_np::Lookup_named_param_def <
internal_np::face_patch_t,
NamedParameters,
typename boost::template property_map<
PolygonMesh, CGAL::dynamic_face_property_t<faces_size_type > >::const_type> ::type
Fpm;
using parameters::choose_parameter;
using parameters::get_parameter;
using parameters::is_default_parameter;
Fpm pidmap = choose_parameter(get_parameter(np, internal_np::face_patch),
get(CGAL::dynamic_face_property_t<faces_size_type>(), tm));
faces_size_type nb_patches = 0;
if(is_default_parameter(get_parameter(np, internal_np::face_patch)))
{
nb_patches = CGAL::Polygon_mesh_processing::connected_components(
tm, pidmap, CGAL::parameters::face_index_map(fim)
.edge_is_constrained_map(ecm));
}
else
{
for(const auto& f : faces(tm))
{
faces_size_type patch_id = get(pidmap, f);
if(patch_id > nb_patches)
nb_patches = patch_id;
}
nb_patches+=1;
}
for(faces_size_type i=0; i<nb_patches; ++i)
{
CGAL::Face_filtered_graph<PolygonMesh, FIMap, VIMap, HIMap>
filter_graph(tm, i, pidmap, CGAL::parameters::face_index_map(fim)
.halfedge_index_map(him)
.vertex_index_map(vim));
range.push_back(PolygonMesh());
PolygonMesh& new_graph = range.back();
CGAL::copy_face_graph(filter_graph, new_graph);
}
}
}//internal
/*!
* \ingroup keep_connected_components_grp
* identifies the connected components of `pm` and pushes back a new `PolygonMesh` for each connected component in `cc_meshes`.
*
*
* \tparam PolygonMesh a model of `FaceListGraph`
* \tparam PolygonMeshRange a model of `SequenceContainer` with `PolygonMesh` as value type.
*
* \tparam NamedParameters a sequence of Named Parameters
*
* \param pm the polygon mesh
* \param cc_meshes container that is filled with the extracted connected components.
* \param np an optional sequence of Named Parameters among the ones listed below
*
* \cgalNamedParamsBegin
* \cgalParamBegin{edge_is_constrained_map} a property map containing the constrained-or-not status of each edge of `pm` \cgalParamEnd
* \cgalParamBegin{face_index_map}
* a property map containing a unique index for each face initialized from 0 to `num_faces(pm)`
* \cgalParamEnd
* \cgalParamBegin{vertex_index_map}
* a property map containing a unique index for each vertex initialized 0 to `num_vertices(pm)`
* \cgalParamEnd
* \cgalNPBegin{halfedge_index_map}
* a property map containing a unique index for each halfedge initialized 0 to `num_halfedges(pm)`
* \cgalNPEnd
* \cgalParamBegin{face_patch_map} a property map with the patch id's associated to the
faces of `pm`. Instance of a class model of `ReadPropertyMap`.
If not provided, an internal map will be filled with a call to
`connected_components()` with `edge_is_constrained_map()` (if provided).
* \cgalParamEnd
* \cgalNamedParamsEnd
*
*/
template <class PolygonMesh, class PolygonMeshRange, class NamedParameters>
void split_connected_components(const PolygonMesh& pm,
PolygonMeshRange& cc_meshes,
const NamedParameters& np)
{
typedef typename internal_np::Lookup_named_param_def <
internal_np::edge_is_constrained_t,
NamedParameters,
internal::No_mark<PolygonMesh>//default
> ::type Ecm;
using parameters::choose_parameter;
using parameters::get_parameter;
Ecm ecm = choose_parameter(get_parameter(np, internal_np::edge_is_constrained),
internal::No_mark<PolygonMesh>());
internal::split_connected_components_impl(CGAL::get_initialized_face_index_map(pm, np),
CGAL::get_initialized_halfedge_index_map(pm, np),
CGAL::get_initialized_vertex_index_map(pm, np),
ecm, cc_meshes, pm, np);
}
template <class PolygonMesh, class PolygonMeshRange>
void split_connected_components(const PolygonMesh& pm,
PolygonMeshRange& cc_meshes)
{
split_connected_components(pm, cc_meshes, parameters::all_default());
}
} // namespace Polygon_mesh_processing
} // namespace CGAL
#include <CGAL/enable_warnings.h>
@@ -17,11 +17,12 @@
#include <CGAL/disable_warnings.h>
#include <CGAL/boost/graph/copy_face_graph.h>
#include <CGAL/boost/graph/named_params_helper.h>
#include <CGAL/Polygon_mesh_processing/intersection.h>
#include <CGAL/Polygon_mesh_processing/internal/Corefinement/Visitor.h>
#include <CGAL/Polygon_mesh_processing/internal/Corefinement/Face_graph_output_builder.h>
#include <CGAL/Polygon_mesh_processing/internal/Corefinement/Output_builder_for_autorefinement.h>
#include <CGAL/boost/graph/copy_face_graph.h>
#include <CGAL/iterator.h>
namespace CGAL {
@@ -34,95 +35,6 @@ using Polygon_mesh_processing::Corefinement::Self_intersection_exception;
namespace Polygon_mesh_processing {
namespace internal {
template <class Kernel, class TriangleMesh, class VD, class Fid_map, class Vpm>
bool recursive_does_bound_a_volume(const TriangleMesh& tm,
Vpm& vpm,
Fid_map& fid_map,
const std::vector<VD>& xtrm_vertices,
boost::dynamic_bitset<>& cc_handled,
const std::vector<std::size_t>& face_cc,
std::size_t xtrm_cc_id,
bool is_parent_outward_oriented)
{
typedef boost::graph_traits<TriangleMesh> GT;
typedef typename GT::face_descriptor face_descriptor;
typedef Side_of_triangle_mesh<TriangleMesh, Kernel, Vpm> Side_of_tm;
// first check that the orientation of the current cc is consistant with its
// parent cc containing it
bool new_is_parent_outward_oriented = internal::is_outward_oriented(
xtrm_vertices[xtrm_cc_id], tm, parameters::vertex_point_map(vpm));
if (new_is_parent_outward_oriented==is_parent_outward_oriented)
return false;
cc_handled.set(xtrm_cc_id);
std::size_t nb_cc = cc_handled.size();
// get all cc that are inside xtrm_cc_id
std::vector<face_descriptor> cc_faces;
for(face_descriptor fd : faces(tm))
{
if(face_cc[get(fid_map, fd)]==xtrm_cc_id)
cc_faces.push_back(fd);
}
typename Side_of_tm::AABB_tree aabb_tree(cc_faces.begin(), cc_faces.end(),
tm, vpm);
Side_of_tm side_of_cc(aabb_tree);
std::vector<std::size_t> cc_inside;
for(std::size_t id=0; id<nb_cc; ++id)
{
if (cc_handled.test(id)) continue;
if (side_of_cc(get(vpm,xtrm_vertices[id]))==ON_BOUNDED_SIDE)
cc_inside.push_back(id);
}
// check whether we need another recursion for cc inside xtrm_cc_id
if (!cc_inside.empty())
{
std::size_t new_xtrm_cc_id = cc_inside.front();
boost::dynamic_bitset<> new_cc_handled(nb_cc,0);
new_cc_handled.set();
new_cc_handled.reset(new_xtrm_cc_id);
cc_handled.set(new_xtrm_cc_id);
std::size_t nb_candidates = cc_inside.size();
for (std::size_t i=1;i<nb_candidates;++i)
{
std::size_t candidate = cc_inside[i];
if(get(vpm,xtrm_vertices[candidate]).z() >
get(vpm,xtrm_vertices[new_xtrm_cc_id]).z()) new_xtrm_cc_id=candidate;
new_cc_handled.reset(candidate);
cc_handled.set(candidate);
}
if ( !internal::recursive_does_bound_a_volume<Kernel>(
tm, vpm, fid_map, xtrm_vertices, new_cc_handled, face_cc,
new_xtrm_cc_id, new_is_parent_outward_oriented) ) return false;
}
// now explore remaining cc included in the same cc as xtrm_cc_id
boost::dynamic_bitset<> cc_not_handled = ~cc_handled;
std::size_t new_xtrm_cc_id = cc_not_handled.find_first();
if (new_xtrm_cc_id == cc_not_handled.npos) return true;
for (std::size_t candidate = cc_not_handled.find_next(new_xtrm_cc_id);
candidate < cc_not_handled.npos;
candidate = cc_not_handled.find_next(candidate))
{
if(get(vpm,xtrm_vertices[candidate]).z() > get(vpm,xtrm_vertices[new_xtrm_cc_id]).z())
new_xtrm_cc_id = candidate;
}
return internal::recursive_does_bound_a_volume<Kernel>(
tm, vpm, fid_map, xtrm_vertices, cc_handled, face_cc,
new_xtrm_cc_id, is_parent_outward_oriented);
}
} //end of namespace internal
namespace Corefinement
{
/** \ingroup PMP_corefinement_grp
@@ -143,101 +55,6 @@ enum Boolean_operation_type {UNION = 0, INTERSECTION=1,
#endif
}
/** \ingroup PMP_corefinement_grp
*
* indicates if `tm` bounds a volume.
* See \ref coref_def_subsec for details.
*
* @tparam TriangleMesh a model of `MutableFaceGraph`, `HalfedgeListGraph` and `FaceListGraph`.
* @tparam NamedParameters a sequence of \ref pmp_namedparameters "Named Parameters"
*
* @param tm a closed triangulated surface mesh
* @param np optional sequence of \ref pmp_namedparameters "Named Parameters" among the ones listed below
*
* @pre `CGAL::is_closed(tm)`
*
* \cgalNamedParamsBegin
* \cgalParamBegin{vertex_point_map}
* the property map with the points associated to the vertices of `tm`.
* If this parameter is omitted, an internal property map for
* `CGAL::vertex_point_t` must be available in `TriangleMesh`
* \cgalParamEnd
* \cgalParamBegin{face_index_map}
* a property map containing the index of each face of `tm`.
* \cgalParamEnd
* \cgalNamedParamsEnd
*
* \see `CGAL::Polygon_mesh_processing::orient_to_bound_a_volume()`
*/
template <class TriangleMesh, class NamedParameters>
bool does_bound_a_volume(const TriangleMesh& tm, const NamedParameters& np)
{
typedef boost::graph_traits<TriangleMesh> GT;
typedef typename GT::vertex_descriptor vertex_descriptor;
typedef typename GetVertexPointMap<TriangleMesh,
NamedParameters>::const_type Vpm;
typedef typename GetFaceIndexMap<TriangleMesh,
NamedParameters>::const_type Fid_map;
typedef typename Kernel_traits<
typename boost::property_traits<Vpm>::value_type >::Kernel Kernel;
if (!is_closed(tm)) return false;
if (!is_triangle_mesh(tm)) return false;
Vpm vpm = parameters::choose_parameter(parameters::get_parameter(np, internal_np::vertex_point),
get_const_property_map(boost::vertex_point, tm));
Fid_map fid_map = parameters::choose_parameter(parameters::get_parameter(np, internal_np::face_index),
get_const_property_map(boost::face_index, tm));
std::vector<std::size_t> face_cc(num_faces(tm), std::size_t(-1));
// set the connected component id of each face
std::size_t nb_cc = connected_components(tm,
bind_property_maps(fid_map,make_property_map(face_cc)),
parameters::face_index_map(fid_map));
if (nb_cc == 1)
return true;
boost::dynamic_bitset<> cc_handled(nb_cc, 0);
// extract a vertex with max z coordinate for each connected component
std::vector<vertex_descriptor> xtrm_vertices(nb_cc, GT::null_vertex());
for(vertex_descriptor vd : vertices(tm))
{
std::size_t cc_id = face_cc[get(fid_map, face(halfedge(vd, tm), tm))];
if (xtrm_vertices[cc_id]==GT::null_vertex())
xtrm_vertices[cc_id]=vd;
else
if (get(vpm, vd).z()>get(vpm,xtrm_vertices[cc_id]).z())
xtrm_vertices[cc_id]=vd;
}
//extract a vertex with max z amongst all components
std::size_t xtrm_cc_id=0;
for(std::size_t id=1; id<nb_cc; ++id)
if (get(vpm, xtrm_vertices[id]).z()>get(vpm,xtrm_vertices[xtrm_cc_id]).z())
xtrm_cc_id=id;
bool is_parent_outward_oriented =
!internal::is_outward_oriented(xtrm_vertices[xtrm_cc_id], tm, np);
return internal::recursive_does_bound_a_volume<Kernel>(tm, vpm, fid_map,
xtrm_vertices,
cc_handled,
face_cc,
xtrm_cc_id,
is_parent_outward_oriented);
}
/// \cond SKIP_IN_MANUAL
template <class TriangleMesh>
bool does_bound_a_volume(const TriangleMesh& tm)
{
return does_bound_a_volume(tm, parameters::all_default());
}
/// \endcond
#define CGAL_COREF_SET_OUTPUT_EDGE_MARK_MAP(I) \
typedef typename internal_np::Lookup_named_param_def < \
@@ -246,9 +63,7 @@ bool does_bound_a_volume(const TriangleMesh& tm)
Corefinement::No_mark<TriangleMesh> \
> ::type Ecm_out_##I; \
Ecm_out_##I ecm_out_##I = \
parameters::choose_parameter( parameters::get_parameter(std::get<I>(nps_out), internal_np::edge_is_constrained), \
Corefinement::No_mark<TriangleMesh>() );
parameters::choose_parameter<Ecm_out_##I>(parameters::get_parameter(std::get<I>(nps_out), internal_np::edge_is_constrained));
/**
* \ingroup PMP_corefinement_grp
@@ -351,8 +166,11 @@ corefine_and_compute_boolean_operations(
NamedParametersOut2,
NamedParametersOut3>& nps_out)
{
using parameters::choose_parameter;
using parameters::get_parameter;
const bool throw_on_self_intersection =
parameters::choose_parameter(parameters::get_parameter(np1, internal_np::throw_on_self_intersection), false);
choose_parameter(get_parameter(np1, internal_np::throw_on_self_intersection), false);
// Vertex point maps
//for input meshes
@@ -365,11 +183,11 @@ corefine_and_compute_boolean_operations(
static const bool same_vpm = (boost::is_same<Vpm,Vpm2>::value); )
CGAL_static_assertion(same_vpm);
Vpm vpm1 = parameters::choose_parameter(parameters::get_parameter(np1, internal_np::vertex_point),
get_property_map(boost::vertex_point, tm1));
Vpm vpm1 = choose_parameter(get_parameter(np1, internal_np::vertex_point),
get_property_map(boost::vertex_point, tm1));
Vpm vpm2 = parameters::choose_parameter(parameters::get_parameter(np2, internal_np::vertex_point),
get_property_map(boost::vertex_point, tm2));
Vpm vpm2 = choose_parameter(get_parameter(np2, internal_np::vertex_point),
get_property_map(boost::vertex_point, tm2));
typedef typename boost::property_traits<Vpm>::value_type Point_3;
@@ -497,10 +315,8 @@ corefine_and_compute_boolean_operations(
Corefinement::No_mark<TriangleMesh>//default
> ::type Ecm2;
Ecm1 ecm1 = parameters::choose_parameter( parameters::get_parameter(np1, internal_np::edge_is_constrained),
Corefinement::No_mark<TriangleMesh>() );
Ecm2 ecm2 = parameters::choose_parameter( parameters::get_parameter(np2, internal_np::edge_is_constrained),
Corefinement::No_mark<TriangleMesh>() );
Ecm1 ecm1 = choose_parameter<Ecm1>(get_parameter(np1, internal_np::edge_is_constrained));
Ecm2 ecm2 = choose_parameter<Ecm2>(get_parameter(np2, internal_np::edge_is_constrained));
typedef Corefinement::Ecm_bind<TriangleMesh, Ecm1, Ecm2> Ecm_in;
@@ -514,34 +330,28 @@ corefine_and_compute_boolean_operations(
typedef std::tuple<Ecm_out_0, Ecm_out_1, Ecm_out_2, Ecm_out_3>
Edge_mark_map_tuple;
// Face index point maps
typedef typename GetFaceIndexMap<TriangleMesh,
NamedParameters1>::type Fid_map;
typedef typename GetFaceIndexMap<TriangleMesh,
NamedParameters2>::type Fid_map2;
CGAL_USE_TYPE(Fid_map2);
CGAL_assertion_code(
static const bool same_fidmap = (boost::is_same<Fid_map,Fid_map2>::value);)
CGAL_static_assertion(same_fidmap);
// Face index point maps
typedef typename CGAL::GetInitializedFaceIndexMap<TriangleMesh, NamedParameters1>::type FaceIndexMap1;
typedef typename CGAL::GetInitializedFaceIndexMap<TriangleMesh, NamedParameters2>::type FaceIndexMap2;
Fid_map fid_map1 = parameters::choose_parameter(parameters::get_parameter(np1, internal_np::face_index),
get_property_map(boost::face_index, tm1));
Fid_map fid_map2 = parameters::choose_parameter(parameters::get_parameter(np2, internal_np::face_index),
get_property_map(boost::face_index, tm2));
// User visitor
FaceIndexMap1 fid_map1 = get_initialized_face_index_map(tm1, np1);
FaceIndexMap2 fid_map2 = get_initialized_face_index_map(tm2, np2);
// User visitor
typedef typename internal_np::Lookup_named_param_def <
internal_np::graph_visitor_t,
NamedParameters1,
Corefinement::Default_visitor<TriangleMesh>//default
> ::type User_visitor;
User_visitor uv( parameters::choose_parameter( parameters::get_parameter(np1, internal_np::graph_visitor),
Corefinement::Default_visitor<TriangleMesh>() ) );
User_visitor uv(choose_parameter<User_visitor>(get_parameter(np1, internal_np::graph_visitor)));
// surface intersection algorithm call
typedef Corefinement::Face_graph_output_builder<TriangleMesh,
Vpm,
Vpm_out_tuple,
Fid_map,
FaceIndexMap1,
FaceIndexMap2,
Default,
Ecm_in,
Edge_mark_map_tuple,
@@ -551,20 +361,17 @@ corefine_and_compute_boolean_operations(
TriangleMesh, Vpm, Ob, Ecm_in, User_visitor> Algo_visitor;
Ecm_in ecm_in(tm1,tm2,ecm1,ecm2);
Edge_mark_map_tuple ecms_out(ecm_out_0, ecm_out_1, ecm_out_2, ecm_out_3);
Ob ob(tm1, tm2, vpm1, vpm2, fid_map1, fid_map2, ecm_in,
vpm_out_tuple, ecms_out, uv, output);
Ob ob(tm1, tm2, vpm1, vpm2, fid_map1, fid_map2, ecm_in, vpm_out_tuple, ecms_out, uv, output);
// special case used for clipping open meshes
if ( parameters::choose_parameter( parameters::get_parameter(np1, internal_np::use_bool_op_to_clip_surface),
false) )
if (choose_parameter(get_parameter(np1, internal_np::use_bool_op_to_clip_surface), false))
{
CGAL_assertion(output[Corefinement::INTERSECTION] != boost::none);
CGAL_assertion(output[Corefinement::UNION] == boost::none);
CGAL_assertion(output[Corefinement::TM1_MINUS_TM2] == boost::none);
CGAL_assertion(output[Corefinement::TM2_MINUS_TM1] == boost::none);
const bool use_compact_clipper =
parameters::choose_parameter( parameters::get_parameter(np1, internal_np::use_compact_clipper),
true);
choose_parameter(get_parameter(np1, internal_np::use_compact_clipper), true);
ob.setup_for_clipping_a_surface(use_compact_clipper);
}
@@ -822,17 +629,20 @@ corefine_and_compute_difference( TriangleMesh& tm1,
* \cgalNamedParamsEnd
*
*/
template <class TriangleMesh,
class NamedParameters1,
class NamedParameters2>
void
corefine( TriangleMesh& tm1,
TriangleMesh& tm2,
const NamedParameters1& np1,
const NamedParameters2& np2)
template <class TriangleMesh,
class NamedParameters1,
class NamedParameters2>
void
corefine( TriangleMesh& tm1,
TriangleMesh& tm2,
const NamedParameters1& np1,
const NamedParameters2& np2)
{
using parameters::choose_parameter;
using parameters::get_parameter;
const bool throw_on_self_intersection =
parameters::choose_parameter(parameters::get_parameter(np1, internal_np::throw_on_self_intersection), false);
choose_parameter(get_parameter(np1, internal_np::throw_on_self_intersection), false);
// Vertex point maps
typedef typename GetVertexPointMap<TriangleMesh,
@@ -844,11 +654,11 @@ corefine_and_compute_difference( TriangleMesh& tm1,
static const bool same_vpm = (boost::is_same<Vpm,Vpm2>::value);)
CGAL_static_assertion(same_vpm);
Vpm vpm1 = parameters::choose_parameter(parameters::get_parameter(np1, internal_np::vertex_point),
get_property_map(boost::vertex_point, tm1));
Vpm vpm1 = choose_parameter(get_parameter(np1, internal_np::vertex_point),
get_property_map(boost::vertex_point, tm1));
Vpm vpm2 = parameters::choose_parameter(parameters::get_parameter(np2, internal_np::vertex_point),
get_property_map(boost::vertex_point, tm2));
Vpm vpm2 = choose_parameter(get_parameter(np2, internal_np::vertex_point),
get_property_map(boost::vertex_point, tm2));
// Edge is-constrained maps
typedef typename internal_np::Lookup_named_param_def <
@@ -863,10 +673,8 @@ corefine_and_compute_difference( TriangleMesh& tm1,
Corefinement::No_mark<TriangleMesh>//default
> ::type Ecm2;
Ecm1 ecm1 = parameters::choose_parameter( parameters::get_parameter(np1, internal_np::edge_is_constrained),
Corefinement::No_mark<TriangleMesh>() );
Ecm2 ecm2 = parameters::choose_parameter( parameters::get_parameter(np2, internal_np::edge_is_constrained),
Corefinement::No_mark<TriangleMesh>() );
Ecm1 ecm1 = choose_parameter<Ecm1>(get_parameter(np1, internal_np::edge_is_constrained));
Ecm2 ecm2 = choose_parameter<Ecm2>(get_parameter(np2, internal_np::edge_is_constrained));
typedef Corefinement::Ecm_bind<TriangleMesh, Ecm1, Ecm2> Ecm;
@@ -883,8 +691,7 @@ corefine_and_compute_difference( TriangleMesh& tm1,
NamedParameters1,
Corefinement::Default_visitor<TriangleMesh>//default
> ::type User_visitor;
User_visitor uv( parameters::choose_parameter( parameters::get_parameter(np1, internal_np::graph_visitor),
Corefinement::Default_visitor<TriangleMesh>() ) );
User_visitor uv(choose_parameter<User_visitor>(get_parameter(np1, internal_np::graph_visitor)));
// surface intersection algorithm call
typedef Corefinement::No_extra_output_from_corefinement<TriangleMesh> Ob;
@@ -926,18 +733,20 @@ namespace experimental {
* \cgalNamedParamsEnd
*
*/
template <class TriangleMesh,
class NamedParameters>
void
autorefine( TriangleMesh& tm,
const NamedParameters& np)
template <class TriangleMesh,
class NamedParameters>
void
autorefine( TriangleMesh& tm,
const NamedParameters& np)
{
// Vertex point maps
typedef typename GetVertexPointMap<TriangleMesh,
NamedParameters>::type Vpm;
using parameters::choose_parameter;
using parameters::get_parameter;
Vpm vpm = parameters::choose_parameter(parameters::get_parameter(np, internal_np::vertex_point),
get_property_map(boost::vertex_point, tm));
// Vertex point maps
typedef typename GetVertexPointMap<TriangleMesh, NamedParameters>::type Vpm;
Vpm vpm = choose_parameter(get_parameter(np, internal_np::vertex_point),
get_property_map(boost::vertex_point, tm));
// Edge is-constrained maps
typedef typename internal_np::Lookup_named_param_def <
@@ -945,10 +754,7 @@ namespace experimental {
NamedParameters,
Corefinement::No_mark<TriangleMesh>//default
> ::type Ecm;
Ecm ecm = parameters::choose_parameter( parameters::get_parameter(np, internal_np::edge_is_constrained),
Corefinement::No_mark<TriangleMesh>() );
Ecm ecm = choose_parameter<Ecm>(get_parameter(np, internal_np::edge_is_constrained));
// User visitor
typedef typename internal_np::Lookup_named_param_def <
@@ -956,8 +762,7 @@ namespace experimental {
NamedParameters,
Corefinement::Default_visitor<TriangleMesh>//default
> ::type User_visitor;
User_visitor uv( parameters::choose_parameter( parameters::get_parameter(np, internal_np::graph_visitor),
Corefinement::Default_visitor<TriangleMesh>() ) );
User_visitor uv(choose_parameter<User_visitor>(get_parameter(np, internal_np::graph_visitor)));
// surface intersection algorithm call
@@ -1002,38 +807,40 @@ namespace experimental {
* \cgalNamedParamsEnd
*
*/
template <class TriangleMesh,
class NamedParameters>
bool
autorefine_and_remove_self_intersections( TriangleMesh& tm,
const NamedParameters& np)
template <class TriangleMesh,
class NamedParameters>
bool
autorefine_and_remove_self_intersections( TriangleMesh& tm,
const NamedParameters& np)
{
using parameters::choose_parameter;
using parameters::get_parameter;
// Vertex point maps
typedef typename GetVertexPointMap<TriangleMesh,
NamedParameters>::type Vpm;
Vpm vpm = parameters::choose_parameter(parameters::get_parameter(np, internal_np::vertex_point),
get_property_map(boost::vertex_point, tm));
typedef typename GetVertexPointMap<TriangleMesh, NamedParameters>::type Vpm;
Vpm vpm = choose_parameter(get_parameter(np, internal_np::vertex_point),
get_property_map(boost::vertex_point, tm));
// Face index map
typedef typename GetFaceIndexMap<TriangleMesh,
NamedParameters>::type Fid_map;
Fid_map fid_map = parameters::choose_parameter(parameters::get_parameter(np, internal_np::face_index),
get_property_map(boost::face_index, tm));
typedef typename GetInitializedFaceIndexMap<TriangleMesh, NamedParameters>::type Fid_map;
Fid_map fid_map = get_initialized_face_index_map(tm, np);
// Edge is-constrained maps
typedef typename internal_np::Lookup_named_param_def <
internal_np::edge_is_constrained_t,
NamedParameters,
Corefinement::No_mark<TriangleMesh>//default
> ::type Ecm;
Ecm ecm = parameters::choose_parameter( parameters::get_parameter(np, internal_np::edge_is_constrained),
Corefinement::No_mark<TriangleMesh>() );
Ecm ecm = choose_parameter<Ecm>(get_parameter(np, internal_np::edge_is_constrained));
// User visitor
typedef typename internal_np::Lookup_named_param_def <
internal_np::graph_visitor_t,
NamedParameters,
Corefinement::Default_visitor<TriangleMesh>//default
> ::type User_visitor;
User_visitor uv( parameters::choose_parameter( parameters::get_parameter(np, internal_np::graph_visitor),
Corefinement::Default_visitor<TriangleMesh>() ) );
User_visitor uv(choose_parameter<User_visitor>(get_parameter(np, internal_np::graph_visitor)));
// surface intersection algorithm call
typedef Corefinement::Output_builder_for_autorefinement<TriangleMesh,
@@ -138,23 +138,18 @@ detect_surface_patches(PolygonMesh& p,
EdgeIsFeatureMap eif,
const NamedParameters& np)
{
//extract types from NPs
typename GetFaceIndexMap<PolygonMesh, NamedParameters>::const_type
fimap = parameters::choose_parameter(parameters::get_parameter(np, internal_np::face_index),
get_const_property_map(boost::face_index, p));
int offset = static_cast<int>(
parameters::choose_parameter(parameters::get_parameter(np, internal_np::first_index),
1));
parameters::choose_parameter(parameters::get_parameter(np, internal_np::first_index), 1));
internal::PatchIdMapWrapper<PatchIdMap,
typename boost::property_traits<PatchIdMap>::value_type>
wrapmap(patch_id_map, offset);
return connected_components(p, wrapmap,
parameters::edge_is_constrained_map(eif)
.face_index_map(fimap));
.face_index_map(CGAL::get_initialized_face_index_map(p, np)));
}
template <typename PolygonMesh, typename EdgeIsFeatureMap, typename PatchIdMap>
typename boost::graph_traits<PolygonMesh>::faces_size_type
detect_surface_patches(PolygonMesh& p,
@@ -373,8 +368,6 @@ namespace internal
* computing a
* surface patch id for each face.
*
* A property map for `CGAL::face_index_t` must be either available
* as an internal property map to `pmesh` or provided as one of the Named Parameters.
*
* \tparam PolygonMesh a model of `FaceGraph`
* \tparam FT a number type. It is
File diff suppressed because it is too large Load Diff
@@ -178,15 +178,16 @@ void extrude_mesh(const InputMesh& input,
using parameters::choose_parameter;
VPMap output_vpm = choose_parameter(get_parameter(np_out, internal_np::vertex_point),
get_property_map(vertex_point, output));
get_property_map(vertex_point, output));
IVPMap input_vpm = choose_parameter(get_parameter(np_in, internal_np::vertex_point),
get_const_property_map(vertex_point, input));
get_const_property_map(vertex_point, input));
std::vector<std::pair<input_vertex_descriptor, output_vertex_descriptor> > bottom_v2v;
std::vector<std::pair<input_halfedge_descriptor, output_halfedge_descriptor> > bottom_h2h;
copy_face_graph(input, output, std::back_inserter(bottom_v2v),
std::back_inserter(bottom_h2h), Emptyset_iterator(),
input_vpm, output_vpm);
copy_face_graph(input, output, parameters::vertex_point_map(input_vpm)
.vertex_to_vertex_output_iterator(std::back_inserter(bottom_v2v))
.halfedge_to_halfedge_output_iterator(std::back_inserter(bottom_h2h)),
parameters::vertex_point_map(output_vpm));
// create the offset for the other side
for(std::size_t i = 0; i< bottom_v2v.size(); ++i)
@@ -198,9 +199,11 @@ void extrude_mesh(const InputMesh& input,
// collect border halfedges for the creation of the triangle strip
std::vector<std::pair<input_vertex_descriptor, output_vertex_descriptor> > top_v2v;
std::vector<std::pair<input_halfedge_descriptor, output_halfedge_descriptor> > top_h2h;
copy_face_graph(input, output, std::inserter(top_v2v, top_v2v.end()),
std::inserter(top_h2h, top_h2h.end()), Emptyset_iterator(),
input_vpm, output_vpm);
copy_face_graph(input, output, parameters::vertex_point_map(input_vpm)
.vertex_to_vertex_output_iterator(std::inserter(top_v2v, top_v2v.end()))
.halfedge_to_halfedge_output_iterator(std::inserter(top_h2h, top_h2h.end())),
parameters::vertex_point_map(output_vpm));
for(std::size_t i = 0; i< top_v2v.size(); ++i)
{
top(top_v2v[i].first, top_v2v[i].second);
@@ -6,7 +6,7 @@
// $URL$
// $Id$
// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
//
//
//
// Author(s) : Ilker O. Yaz
@@ -91,7 +91,7 @@ namespace internal {
\cgalParamBegin{fairing_continuity} tangential continuity of the output surface patch. The larger `fairing_continuity` gets, the more fixed vertices are required \cgalParamEnd
\cgalParamBegin{sparse_linear_solver} an instance of the sparse linear solver used for fairing \cgalParamEnd
\cgalNamedParamsEnd
@return `true` if fairing is successful, otherwise no vertices are relocated
@pre `is_triangle_mesh(tmesh)`
@@ -141,14 +141,13 @@ namespace internal {
Default_Weight_calculator;
VPMap vpmap_ = choose_parameter(get_parameter(np, internal_np::vertex_point),
get_property_map(vertex_point, tmesh));
get_property_map(vertex_point, tmesh));
return internal::fair(tmesh, vertices,
choose_parameter(get_parameter(np, internal_np::sparse_linear_solver), Default_solver()),
choose_parameter<Default_solver>(get_parameter(np, internal_np::sparse_linear_solver)),
choose_parameter(get_parameter(np, internal_np::weight_calculator), Default_Weight_calculator(tmesh, vpmap_)),
choose_parameter(get_parameter(np, internal_np::fairing_continuity), 1),
vpmap_
);
vpmap_);
}
template<typename TriangleMesh, typename VertexRange>
@@ -51,9 +51,9 @@ class Traversal_traits_with_transformation_helper
AK::Aff_transformation_3 a_at = c2f(at);
AK::FT xtrm[6] = { c2f(bbox.min(0)), c2f(bbox.max(0)),
c2f(bbox.min(1)), c2f(bbox.max(1)),
c2f(bbox.min(2)), c2f(bbox.max(2)) };
AK::FT xtrm[6] = { c2f((bbox.min)(0)), c2f((bbox.max)(0)),
c2f((bbox.min)(1)), c2f((bbox.max)(1)),
c2f((bbox.min)(2)), c2f((bbox.max)(2)) };
typename AK::Point_3 ps[8];
ps[0] = a_at( AK::Point_3(xtrm[0], xtrm[2], xtrm[4]) );
@@ -83,9 +83,9 @@ class Traversal_traits_with_transformation_helper
AK::Aff_transformation_3 a_at = c2f(at);
AK::FT xtrm[6] = { c2f(bbox.min(0)), c2f(bbox.max(0)),
c2f(bbox.min(1)), c2f(bbox.max(1)),
c2f(bbox.min(2)), c2f(bbox.max(2)) };
AK::FT xtrm[6] = { c2f((bbox.min)(0)), c2f((bbox.max)(0)),
c2f((bbox.min)(1)), c2f((bbox.max)(1)),
c2f((bbox.min)(2)), c2f((bbox.max)(2)) };
typename AK::Point_3 ps[2];
ps[0] = a_at( AK::Point_3(xtrm[0], xtrm[2], xtrm[4]) );
@@ -60,7 +60,8 @@ namespace params=PMP::parameters;
template <class TriangleMesh,
class VertexPointMap,
class VpmOutTuple,
class FaceIdMap,
class FaceIdMap1,
class FaceIdMap2,
class Kernel_=Default,
class EdgeMarkMapBind_ = Default,
class EdgeMarkMapTuple_ = Default,
@@ -110,8 +111,10 @@ class Face_graph_output_builder
//Data members
TriangleMesh &tm1, &tm2;
// property maps of input meshes
const VertexPointMap &vpm1, &vpm2;
const FaceIdMap &fids1, &fids2;
const VertexPointMap vpm1;
const VertexPointMap vpm2;
FaceIdMap1 fids1;
FaceIdMap2 fids2;
EdgeMarkMapBind& marks_on_input_edges;
// property maps of output meshes
const VpmOutTuple& output_vpms;
@@ -215,7 +218,7 @@ class Face_graph_output_builder
// detect if a polyline is incident to two patches that won't be imported
// for the current operation (polylines skipt are always incident to a
// coplanar patch)
template <class TM, class FIM>
template <class TM, class FIM1, class FIM2>
static
void fill_polylines_to_skip(
Intersection_polylines& polylines,
@@ -223,8 +226,8 @@ class Face_graph_output_builder
const std::vector<std::size_t>& tm2_patch_ids,
const boost::dynamic_bitset<>& patches_of_tm1_used,
const boost::dynamic_bitset<>& patches_of_tm2_used,
const FIM& fids1,
const FIM& fids2,
const FIM1 fids1,
const FIM2 fids2,
const TM& tm1,
const TM& tm2)
{
@@ -341,18 +344,17 @@ class Face_graph_output_builder
public:
Face_graph_output_builder( TriangleMesh& tm1,
TriangleMesh& tm2,
const VertexPointMap &vpm1,
const VertexPointMap &vpm2,
const FaceIdMap& fids1,
const FaceIdMap& fids2,
EdgeMarkMapBind& marks_on_input_edges,
Face_graph_output_builder(TriangleMesh& tm1,
TriangleMesh& tm2,
const VertexPointMap vpm1,
const VertexPointMap vpm2,
FaceIdMap1 fids1,
FaceIdMap2 fids2,
EdgeMarkMapBind& marks_on_input_edges,
const VpmOutTuple& output_vpms,
EdgeMarkMapTuple& out_edge_mark_maps,
UserVisitor& user_visitor,
const std::array<
boost::optional<TriangleMesh*>, 4 >& requested_output)
EdgeMarkMapTuple& out_edge_mark_maps,
UserVisitor& user_visitor,
const std::array<boost::optional<TriangleMesh*>, 4 >& requested_output)
: tm1(tm1), tm2(tm2)
, vpm1(vpm1), vpm2(vpm2)
, fids1(fids1), fids2(fids2)
@@ -460,9 +462,20 @@ public:
Intersection_edge_map& intersection_edges1 = mesh_to_intersection_edges[&tm1];
Intersection_edge_map& intersection_edges2 = mesh_to_intersection_edges[&tm2];
// this will initialize face indices if the face index map is writable.
helpers::init_face_indices(tm1, fids1);
helpers::init_face_indices(tm2, fids2);
// The property map must be either writable or well-initialized
if( CGAL::internal::Is_writable_property_map<FaceIdMap1>::value &&
!BGL::internal::is_index_map_valid(fids1, num_faces(tm1), faces(tm1)) )
{
BGL::internal::initialize_face_index_map(fids1, tm1);
}
CGAL_assertion(BGL::internal::is_index_map_valid(fids1, num_faces(tm1), faces(tm1)));
if( CGAL::internal::Is_writable_property_map<FaceIdMap2>::value &&
!BGL::internal::is_index_map_valid(fids2, num_faces(tm2), faces(tm2)) )
{
BGL::internal::initialize_face_index_map(fids2, tm2);
}
CGAL_assertion(BGL::internal::is_index_map_valid(fids2, num_faces(tm2), faces(tm2)));
// bitset to identify coplanar faces
boost::dynamic_bitset<> tm1_coplanar_faces(num_faces(tm1), 0);
@@ -620,9 +633,8 @@ public:
std::size_t nb_patches_tm1 =
PMP::connected_components(tm1,
bind_property_maps(fids1,make_property_map(&tm1_patch_ids[0])),
params::edge_is_constrained_map(
is_marked_1)
.face_index_map(fids1));
params::edge_is_constrained_map(is_marked_1)
.face_index_map(fids1));
std::vector <std::size_t> tm1_patch_sizes(nb_patches_tm1, 0);
for(std::size_t i : tm1_patch_ids)
@@ -634,9 +646,8 @@ public:
std::size_t nb_patches_tm2 =
PMP::connected_components(tm2,
bind_property_maps(fids2,make_property_map(&tm2_patch_ids[0])),
params::edge_is_constrained_map(
is_marked_2)
.face_index_map(fids2));
params::edge_is_constrained_map(is_marked_2)
.face_index_map(fids2));
std::vector <std::size_t> tm2_patch_sizes(nb_patches_tm2, 0);
for(Node_id i : tm2_patch_ids)
@@ -1245,9 +1256,8 @@ public:
polyline_lengths.push_back(polyline_info.second+1);
}
typedef Patch_container<TriangleMesh,
FaceIdMap,
Intersection_edge_map> Patches;
typedef Patch_container<TriangleMesh, FaceIdMap1, Intersection_edge_map> Patches1;
typedef Patch_container<TriangleMesh, FaceIdMap2, Intersection_edge_map> Patches2;
boost::unordered_set<vertex_descriptor> border_nm_vertices; // only used if used_to_clip_a_surface == true
if (used_to_clip_a_surface)
@@ -1279,8 +1289,8 @@ public:
}
//store the patch description in a container to avoid recomputing it several times
Patches patches_of_tm1( tm1, tm1_patch_ids, fids1, intersection_edges1, nb_patches_tm1),
patches_of_tm2( tm2, tm2_patch_ids, fids2, intersection_edges2, nb_patches_tm2);
Patches1 patches_of_tm1(tm1, tm1_patch_ids, fids1, intersection_edges1, nb_patches_tm1);
Patches2 patches_of_tm2(tm2, tm2_patch_ids, fids2, intersection_edges2, nb_patches_tm2);
// for each boolean operation, define two bitsets of patches contributing
// to the result
@@ -1460,11 +1470,12 @@ public:
// operation in tm1 with removal (and optionally inside-out) delayed
// First backup the border edges of patches to be used
Patches tmp_patches_of_tm1(tm1,
patches_of_tm1.patch_ids,
patches_of_tm1.fids,
patches_of_tm1.is_intersection_edge,
patches_of_tm1.patches.size());
Patches1 tmp_patches_of_tm1(tm1,
patches_of_tm1.patch_ids,
patches_of_tm1.fids,
patches_of_tm1.is_intersection_edge,
patches_of_tm1.patches.size());
boost::dynamic_bitset<> patches_of_tm1_removed =
~patches_of_tm1_used[inplace_operation_tm1];
for (std::size_t i = patches_of_tm1_removed.find_first();
@@ -214,9 +214,6 @@ public:
const boost::dynamic_bitset<>& is_node_of_degree_one,
const Mesh_to_map_node&)
{
// this will initialize face indices if the face index map is writable.
helpers::init_face_indices(tm, fids);
// first build an unordered_map mapping a vertex to its node id + a set
// of all intersection edges
typedef boost::unordered_set<edge_descriptor> Intersection_edge_map;
@@ -239,8 +236,13 @@ public:
intersection_edges.insert(edge(p.second.h2, tm));
}
// this will initialize face indices if the face index map is writable.
helpers::init_face_indices(tm, fids);
// The property map must be either writable or well-initialized
if( CGAL::internal::Is_writable_property_map<FaceIdMap>::value &&
!BGL::internal::is_index_map_valid(fids, num_faces(tm), faces(tm)) )
{
BGL::internal::initialize_face_index_map(fids, tm);
}
CGAL_assertion(BGL::internal::is_index_map_valid(fids, num_faces(tm), faces(tm)));
// bitset to identify coplanar faces
boost::dynamic_bitset<> tm_coplanar_faces(num_faces(tm), 0);
@@ -424,7 +424,7 @@ template <class PolygonMesh, class FaceIndexMap, class IsIntersectionEdge>
void extract_patch_simplices(
std::size_t patch_id,
PolygonMesh& pm,
const FaceIndexMap& fids,
const FaceIndexMap fids,
const std::vector<std::size_t>& patch_ids,
std::vector<typename boost::graph_traits<PolygonMesh>::face_descriptor>& patch_faces,
std::set<typename boost::graph_traits<PolygonMesh>::vertex_descriptor>& interior_vertices,
@@ -480,13 +480,13 @@ struct Patch_container{
// external data members
PolygonMesh& pm;
const std::vector<std::size_t>& patch_ids;
const FaceIndexMap& fids;
const FaceIndexMap fids;
const IsIntersectionEdge& is_intersection_edge;
// constructor
Patch_container(
PolygonMesh& pm,
const std::vector<std::size_t>& patch_ids,
const FaceIndexMap& fids,
const FaceIndexMap fids,
const IsIntersectionEdge& is_intersection_edge,
std::size_t nb_patches
) : patches(nb_patches)
@@ -1010,14 +1010,15 @@ template < class TriangleMesh,
class EdgeMarkMap2,
class EdgeMarkMapOut,
class IntersectionPolylines,
class PatchContainer,
class PatchContainer1,
class PatchContainer2,
class UserVisitor>
void fill_new_triangle_mesh(
TriangleMesh& output,
const boost::dynamic_bitset<>& patches_of_tm1_to_import,
const boost::dynamic_bitset<>& patches_of_tm2_to_import,
PatchContainer& patches_of_tm1,
PatchContainer& patches_of_tm2,
PatchContainer1& patches_of_tm1,
PatchContainer2& patches_of_tm2,
bool reverse_orientation_of_patches_from_tm1,
bool reverse_orientation_of_patches_from_tm2,
const IntersectionPolylines& polylines,
@@ -1256,7 +1257,8 @@ void disconnect_patches(
}
template <class TriangleMesh,
class PatchContainer,
class PatchContainer1,
class PatchContainer2,
class IntersectionPolylines,
class EdgeMap,
class VertexPointMap,
@@ -1269,8 +1271,8 @@ void compute_inplace_operation_delay_removal_and_insideout(
TriangleMesh& tm2,
const boost::dynamic_bitset<>& patches_of_tm1_to_keep,
const boost::dynamic_bitset<>& patches_of_tm2_to_import,
PatchContainer& patches_of_tm1,
PatchContainer& patches_of_tm2,
PatchContainer1& patches_of_tm1,
PatchContainer2& patches_of_tm2,
bool reverse_patch_orientation_tm2,
const IntersectionPolylines& polylines,
const VertexPointMap& vpm1,
@@ -1416,7 +1418,8 @@ remove_patches(TriangleMesh& tm,
}
template <class TriangleMesh,
class PatchContainer,
class PatchContainer1,
class PatchContainer2,
class VertexPointMap,
class EdgeMarkMapIn1,
class EdgeMarkMapIn2,
@@ -1427,8 +1430,8 @@ void compute_inplace_operation(
const TriangleMesh& /*tm2*/,
const boost::dynamic_bitset<>& patches_of_tm1_to_keep,
const boost::dynamic_bitset<>& patches_of_tm2_to_import,
PatchContainer& patches_of_tm1,
PatchContainer& patches_of_tm2,
PatchContainer1& patches_of_tm1,
PatchContainer2& patches_of_tm2,
bool reverse_patch_orientation_tm1,
bool reverse_patch_orientation_tm2,
const VertexPointMap& vpm1,
@@ -1485,14 +1488,15 @@ void compute_inplace_operation(
template <class TriangleMesh,
class IntersectionPolylines,
class PatchContainer,
class PatchContainer1,
class PatchContainer2,
class EdgeMap>
void compute_border_edge_map(
const TriangleMesh& tm1,
const TriangleMesh& tm2,
const IntersectionPolylines& polylines,
PatchContainer& patches_of_tm1,
PatchContainer& patches_of_tm2,
PatchContainer1& patches_of_tm1,
PatchContainer2& patches_of_tm2,
EdgeMap& tm2_edge_to_tm1_edge)
{
typedef boost::graph_traits<TriangleMesh> GT;
@@ -1520,7 +1524,8 @@ void compute_border_edge_map(
template <class TriangleMesh,
class PatchContainer,
class PatchContainer1,
class PatchContainer2,
class IntersectionPolylines,
class VertexPointMap,
class EdgeMarkMapIn1,
@@ -1532,8 +1537,8 @@ void compute_inplace_operation(
const TriangleMesh& tm2,
const boost::dynamic_bitset<>& patches_of_tm1_to_keep,
const boost::dynamic_bitset<>& patches_of_tm2_to_import,
PatchContainer& patches_of_tm1,
PatchContainer& patches_of_tm2,
PatchContainer1& patches_of_tm1,
PatchContainer2& patches_of_tm2,
bool reverse_patch_orientation_tm1,
bool reverse_patch_orientation_tm2,
const VertexPointMap& vpm1,
@@ -304,7 +304,7 @@ public:
}
if ( abcq==COPLANAR &&
is_edge_target_incident_to_face(opposite(eh, tm), fh) )
is_edge_target_incident_to_face(opposite(eh, tm), fh) )
{
return; // no intersection (incident edge)
}
@@ -1268,6 +1268,8 @@ public:
: nodes(tm1, tm2, vpm1, vpm2)
, visitor(v)
{
CGAL_precondition(is_triangle_mesh(tm1));
CGAL_precondition(is_triangle_mesh(tm2));
CGAL_assertion_code( doing_autorefinement=false; )
}
@@ -1278,6 +1280,7 @@ public:
: nodes(tm, tm, vpm, vpm)
, visitor(v)
{
CGAL_precondition(is_triangle_mesh(tm));
CGAL_assertion_code( doing_autorefinement=true; )
}
@@ -6,7 +6,7 @@
// $URL$
// $Id$
// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
//
//
//
// Author(s) : Ilker O. Yaz
@@ -30,7 +30,7 @@ namespace internal {
template<class PolygonMesh, class OutputIterator>
struct Tracer_polyhedron
struct Tracer_polyhedron
{
typedef typename boost::graph_traits<PolygonMesh>::halfedge_descriptor halfedge_descriptor;
@@ -41,8 +41,8 @@ struct Tracer_polyhedron
{ }
template <class LookupTable>
halfedge_descriptor
operator()(const LookupTable& lambda,
halfedge_descriptor
operator()(const LookupTable& lambda,
int i, int k,
bool last = true)
{
@@ -54,14 +54,14 @@ struct Tracer_polyhedron
{
h = P[i+1];
Euler::fill_hole(h,pmesh); }
else
else
{ h = Euler::add_face_to_border(prev(P[i+1],pmesh), P[i+2/*k*/], pmesh); }
CGAL_assertion(face(h,pmesh) != boost::graph_traits<PolygonMesh>::null_face());
*out++ = face(h,pmesh);
return opposite(h,pmesh);
}
else
}
else
{
int la = lambda.get(i, k);
h = operator()(lambda, i, la, false);
@@ -72,7 +72,7 @@ struct Tracer_polyhedron
h = g;
Euler::fill_hole(g,pmesh);
}
else
else
{ h = Euler::add_face_to_border(prev(h,pmesh), g, pmesh); }
CGAL_assertion(face(h,pmesh) != boost::graph_traits<PolygonMesh>::null_face());
@@ -88,8 +88,8 @@ struct Tracer_polyhedron
// This function is used in test cases (since it returns not just OutputIterator but also Weight)
template<class PolygonMesh, class OutputIterator, class VertexPointMap, class Kernel>
std::pair<OutputIterator, CGAL::internal::Weight_min_max_dihedral_and_area>
triangulate_hole_polygon_mesh(PolygonMesh& pmesh,
std::pair<OutputIterator, CGAL::internal::Weight_min_max_dihedral_and_area>
triangulate_hole_polygon_mesh(PolygonMesh& pmesh,
typename boost::graph_traits<PolygonMesh>::halfedge_descriptor border_halfedge,
OutputIterator out,
VertexPointMap vpmap,
@@ -100,7 +100,7 @@ triangulate_hole_polygon_mesh(PolygonMesh& pmesh,
typedef typename boost::graph_traits<PolygonMesh>::vertex_descriptor vertex_descriptor;
typedef typename boost::graph_traits<PolygonMesh>::halfedge_descriptor halfedge_descriptor;
typedef typename Kernel::Point_3 Point_3;
typedef std::map<vertex_descriptor, int> Vertex_map;
typedef typename Vertex_map::iterator Vertex_map_it;
@@ -128,7 +128,7 @@ triangulate_hole_polygon_mesh(PolygonMesh& pmesh,
CGAL::internal::Weight_min_max_dihedral_and_area::NOT_VALID());
}
} while (++circ != done);
// existing_edges contains neighborhood information between boundary vertices
// more precisely if v_i is neighbor to any other vertex than v_(i-1) and v_(i+1),
// this edge is put into existing_edges
@@ -173,7 +173,7 @@ triangulate_hole_polygon_mesh(PolygonMesh& pmesh,
// fill hole using polyline function, with custom tracer for PolygonMesh
Tracer_polyhedron<PolygonMesh, OutputIterator>
tracer(out, pmesh, P_edges);
CGAL::internal::Weight_min_max_dihedral_and_area weight =
CGAL::internal::Weight_min_max_dihedral_and_area weight =
triangulate_hole_polyline(P, Q, tracer, WC(is_valid),
use_delaunay_triangulation, k)
#ifdef CGAL_USE_WEIGHT_INCOMPLETE
@@ -6,7 +6,7 @@
// $URL$
// $Id$
// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
//
//
//
// Author(s) : Ilker O. Yaz
@@ -62,11 +62,11 @@ private:
bool bound_check(int i, int j) const {
CGAL_assertion(i >= 0 && i < n);
CGAL_assertion(j >= 0 && j < n);
CGAL_assertion(i < j);
CGAL_assertion(i < j);
CGAL_USE(i);
CGAL_USE(j);
// previous implementation was based on directly vector and i supposed to be always smaller than j.
// this check actually can be removed and i =min(i,j) j = max(i,j) can be used for reflexive access
// this check actually can be removed and i =min(i,j) j = max(i,j) can be used for reflexive access
return true;
}
std::vector<T> table;
@@ -81,12 +81,12 @@ public:
void put(int i, int j, const T& t) {
CGAL_assertion(bound_check(i,j));
if(t == default_) {
table.erase(std::make_pair(i,j));
return;
}
std::pair<typename Map::iterator, bool> inserted = table.insert(std::make_pair(std::make_pair(i,j), t));
if(!inserted.second) { inserted.first->second = t;}
}
@@ -107,14 +107,14 @@ public:
// 5-6
typename Map::iterator it;
if(b == 0) { it = table.begin(); }
else {
else {
it = table.upper_bound(std::make_pair(b-1, n)); // to find first entry where entry.first == b
}
while(it != table.end() && it->first.first != e) {
if(it->first.second <= e)
if(it->first.second <= e)
{ table.erase(it++); }
else
else
{ ++it; }
}
}
@@ -140,8 +140,8 @@ private:
struct Is_valid_existing_edges
{
typedef std::vector<std::pair<int, int> > Edge_container;
Is_valid_existing_edges(Edge_container& existing_edges)
Is_valid_existing_edges(Edge_container& existing_edges)
: existing_edges(existing_edges)
{
std::sort(existing_edges.begin(), existing_edges.end());
@@ -155,7 +155,7 @@ struct Is_valid_existing_edges
template<class Point_3>
bool operator()(const std::vector<Point_3>&,
int v0, int v1, int v2) const
int v0, int v1, int v2) const
{
CGAL_assertion(v0 < v1 && v1 < v2);
@@ -166,7 +166,7 @@ struct Is_valid_existing_edges
if(v1 + 1 != v2 &&
std::binary_search(existing_edges.begin(), existing_edges.end(), std::make_pair(v1,v2)) )
{ return false; }
if(std::binary_search(existing_edges.begin(), existing_edges.end(), std::make_pair(v2,v0)) )
{ return false; }
@@ -180,7 +180,7 @@ struct Is_not_degenerate_triangle
{
template<class Point_3>
bool operator()(const std::vector<Point_3>& P,
int v0, int v1, int v2) const
int v0, int v1, int v2) const
{
return !CGAL::collinear(P[v0], P[v1], P[v2]);
}
@@ -194,7 +194,7 @@ struct Is_valid_existing_edges_and_degenerate_triangle
template<class Point_3>
bool operator()(const std::vector<Point_3>& P,
int v0, int v1, int v2) const
int v0, int v1, int v2) const
{
return Is_not_degenerate_triangle()(P,v0,v1,v2)
&& is_valid_edges(P,v0,v1,v2);
@@ -208,7 +208,7 @@ struct Is_valid_existing_edges_and_degenerate_triangle
************************************************************************/
// First minimizes the worst dihedral angle between patch triangles, then the total surface area as a tiebreaker.
class Weight_min_max_dihedral_and_area
class Weight_min_max_dihedral_and_area
{
template<class Weight_, class IsValid>
friend struct Weight_calculator;
@@ -224,22 +224,22 @@ public:
// below required by Weight concept
private:
template<class Point_3, class LookupTable>
Weight_min_max_dihedral_and_area(const std::vector<Point_3>& P,
const std::vector<Point_3>& Q,
int i, int j, int k,
Weight_min_max_dihedral_and_area(const std::vector<Point_3>& P,
const std::vector<Point_3>& Q,
int i, int j, int k,
const LookupTable& lambda)
{
CGAL_assertion(i < j);
CGAL_assertion(j < k);
int n = static_cast<int>(P.size()) -1; // because the first and last point are equal
// The CGAL::dihedral angle is measured between the oriented triangles, that is it goes from [-pi, pi]
// What we need is the angle between the normals of the triangles between [0, pi]
double ang_max = 0;
// Test each edge
int vertices[] = {i, j, k};
for(int e = 0; e < 3; ++e)
for(int e = 0; e < 3; ++e)
{
int v0 = vertices[e];
int v1 = vertices[(e+1)%3];
@@ -247,25 +247,25 @@ private:
double angle = 0;
// check whether the edge is border
if( (v0 + 1 == v1 || (v0 == n-1 && v1 == 0) ) && !Q.empty() ) {
angle = 180 - CGAL::abs(
angle = 180 - CGAL::abs(
to_double(CGAL::approximate_dihedral_angle(P[v0],P[v1],P[v_other],Q[v0])) );
}
else {
if(e == 2) { continue; }
if(lambda.get(v0, v1) != -1){
const Point_3& p01 = P[lambda.get(v0, v1)];
angle = 180 - CGAL::abs(
angle = 180 - CGAL::abs(
to_double(CGAL::approximate_dihedral_angle(P[v0],P[v1],P[v_other],p01)) );
}
}
ang_max = (std::max)(ang_max, angle);
}
w = std::make_pair(ang_max, to_double(CGAL::approximate_sqrt(CGAL::squared_area(P[i],P[j],P[k]))));
}
public:
Weight_min_max_dihedral_and_area operator+(const Weight_min_max_dihedral_and_area& w2) const
Weight_min_max_dihedral_and_area operator+(const Weight_min_max_dihedral_and_area& w2) const
{
CGAL_assertion((*this) != NOT_VALID());
CGAL_assertion(w2 != NOT_VALID());
@@ -283,15 +283,15 @@ public:
return w.first < w2.w.first;
}
bool operator==(const Weight_min_max_dihedral_and_area& w2) const
bool operator==(const Weight_min_max_dihedral_and_area& w2) const
{ return w.first == w2.w.first && w.second == w2.w.second; }
bool operator!=(const Weight_min_max_dihedral_and_area& w2) const
bool operator!=(const Weight_min_max_dihedral_and_area& w2) const
{ return !(*this == w2); }
static const Weight_min_max_dihedral_and_area DEFAULT() // rule: x + DEFAULT() == x
{ return Weight_min_max_dihedral_and_area(0,0); }
static const Weight_min_max_dihedral_and_area NOT_VALID()
static const Weight_min_max_dihedral_and_area NOT_VALID()
{ return Weight_min_max_dihedral_and_area(-1,-1); }
friend std::ostream& operator<<(std::ostream& out, const Weight_min_max_dihedral_and_area& w) {
@@ -311,9 +311,9 @@ private:
Weight_total_edge(double total_length = 0) : total_length(total_length) { }
template<class Point_3, class LookupTable>
Weight_total_edge(const std::vector<Point_3>& P,
const std::vector<Point_3>&,
int i, int j, int k,
Weight_total_edge(const std::vector<Point_3>& P,
const std::vector<Point_3>&,
int i, int j, int k,
const LookupTable&)
: total_length(0)
{
@@ -323,11 +323,11 @@ private:
// Test each edge
int vertices[] = {i, j, k};
for(int e = 0; e < 3; ++e)
for(int e = 0; e < 3; ++e)
{
int v0 = vertices[e];
int v1 = vertices[(e+1)%3];
// check whether the edge is border
bool border = (v0 + 1 == v1) || (v0 == n-1 && v1 == 0);
if(!border) {
@@ -337,7 +337,7 @@ private:
}
public:
Weight_total_edge operator+(const Weight_total_edge& w2) const
Weight_total_edge operator+(const Weight_total_edge& w2) const
{
CGAL_assertion((*this) != NOT_VALID());
CGAL_assertion(w2 != NOT_VALID());
@@ -345,15 +345,15 @@ public:
}
bool operator<(const Weight_total_edge& w2) const
{
{
CGAL_assertion((*this) != NOT_VALID());
CGAL_assertion(w2 != NOT_VALID());
return total_length < w2.total_length;
return total_length < w2.total_length;
}
bool operator==(const Weight_total_edge& w2) const
bool operator==(const Weight_total_edge& w2) const
{ return total_length == w2.total_length; }
bool operator!=(const Weight_total_edge& w2) const
bool operator!=(const Weight_total_edge& w2) const
{ return !(*this == w2); }
static const Weight_total_edge DEFAULT() { return Weight_total_edge(0); } // rule: x + DEFAULT() == x
@@ -373,19 +373,19 @@ class Weight_incomplete
private:
template<class Point_3, class LookupTable>
Weight_incomplete(const std::vector<Point_3>& P,
const std::vector<Point_3>& Q,
int i, int j, int k,
Weight_incomplete(const std::vector<Point_3>& P,
const std::vector<Point_3>& Q,
int i, int j, int k,
const LookupTable& lambda)
: weight(P,Q,i,j,k,lambda), patch_size(1)
{ }
Weight_incomplete(const ActualWeight& weight, int patch_size)
Weight_incomplete(const ActualWeight& weight, int patch_size)
: weight(weight), patch_size(patch_size)
{ }
public:
Weight_incomplete operator+(const Weight_incomplete& w2) const
Weight_incomplete operator+(const Weight_incomplete& w2) const
{
CGAL_assertion((*this) != NOT_VALID());
CGAL_assertion(w2 != NOT_VALID());
@@ -404,10 +404,10 @@ public:
return patch_size > w2.patch_size; // if patch size is larger, then the weight is smaller
}
bool operator==(const Weight_incomplete& w2) const
bool operator==(const Weight_incomplete& w2) const
{ return weight == w2.weight && patch_size == w2.patch_size; }
bool operator!=(const Weight_incomplete& w2) const
bool operator!=(const Weight_incomplete& w2) const
{ return !(*this == w2); }
static const Weight_incomplete DEFAULT() // rule: x + DEFAULT() == x
@@ -426,18 +426,18 @@ public:
// Weight calculator class is both responsible from calculating weights, and checking validity of triangle
template<class Weight_, class IsValid>
struct Weight_calculator
struct Weight_calculator
{
typedef Weight_ Weight;
Weight_calculator(const IsValid& is_valid = IsValid()) : is_valid(is_valid) { }
template<class Point_3, class LookupTable>
Weight operator()(const std::vector<Point_3>& P,
const std::vector<Point_3>& Q,
int i, int j, int k,
const LookupTable& lambda) const
Weight operator()(const std::vector<Point_3>& P,
const std::vector<Point_3>& Q,
int i, int j, int k,
const LookupTable& lambda) const
{
if( !is_valid(P,i,j,k) )
if( !is_valid(P,i,j,k) )
{ return Weight::NOT_VALID(); }
return Weight(P, Q, i,j,k, lambda);
}
@@ -450,7 +450,7 @@ struct Weight_calculator
// It can produce a patch from both complete and incomplete lambda
template<class OutputIteratorValueType, class OutputIteratorPatch, class OutputIteratorHole>
struct Tracer_polyline_incomplete {
Tracer_polyline_incomplete(OutputIteratorPatch out, OutputIteratorHole out_hole)
Tracer_polyline_incomplete(OutputIteratorPatch out, OutputIteratorHole out_hole)
: out(out), out_hole(out_hole)
{ }
@@ -463,7 +463,7 @@ struct Tracer_polyline_incomplete {
ranges.push(std::make_pair(v0, v1));
while(!ranges.empty()) {
std::pair<int, int> r = ranges.top();
std::pair<int, int> r = ranges.top();
ranges.pop();
CGAL_assertion(r.first >= 0 && r.first < n);
CGAL_assertion(r.second >= 0 && r.second < n);
@@ -498,7 +498,7 @@ template<unsigned int Dimension, class Triangulator>
struct Incident_facet_circulator;
template<class Triangulator>
struct Incident_facet_circulator_base
struct Incident_facet_circulator_base
{
typedef typename Triangulator::Facet Facet;
typedef typename Triangulator::Edge Edge;
@@ -511,7 +511,7 @@ struct Incident_facet_circulator_base
{ return (std::min)(e.first->vertex(e.second)->info(), e.first->vertex(e.third)->info()); }
int vertex_second()
{ return (std::max)(e.first->vertex(e.second)->info(), e.first->vertex(e.third)->info()); }
Edge e;
};
@@ -525,7 +525,7 @@ struct Incident_facet_circulator_base
if(i == f.second) { continue; } // skip the vertex which is not on `f`
int f3 = f.first->vertex(i)->info();
if(f3 != v0_info && f3 != v1_info) {
return f3;
return f3;
}
}
CGAL_assertion(false);
@@ -535,8 +535,8 @@ struct Incident_facet_circulator_base
Edge_wrapper edge_first(Facet f, Edge e) {
int v0_info = (std::min)(e.first->vertex(e.second)->info(),
e.first->vertex(e.third)->info());
return Edge(f.first,
get_vertex_index(f.first, v0_info) ,
return Edge(f.first,
get_vertex_index(f.first, v0_info) ,
get_vertex_index(f.first, get_third(f,e)));
}
@@ -587,9 +587,9 @@ struct Incident_facet_circulator<2, Triangulator>
int get_third()
{ return Incident_facet_circulator_base<Triangulator>::get_third(it, e); }
Edge_wrapper edge_first()
Edge_wrapper edge_first()
{ return Incident_facet_circulator_base<Triangulator>::edge_first(it, e); }
Edge_wrapper edge_second()
Edge_wrapper edge_second()
{ return Incident_facet_circulator_base<Triangulator>::edge_second(it, e); }
Facet f1, f2, it;
@@ -615,12 +615,12 @@ struct Incident_facet_circulator<3, Triangulator>
return *this;
}
operator bool() const { return it != end; }
int get_third()
int get_third()
{ return Incident_facet_circulator_base<Triangulator>::get_third(*it, e); }
Edge_wrapper edge_first()
Edge_wrapper edge_first()
{ return Incident_facet_circulator_base<Triangulator>::edge_first(*it, e); }
Edge_wrapper edge_second()
Edge_wrapper edge_second()
{ return Incident_facet_circulator_base<Triangulator>::edge_second(*it, e); }
Facet_circulator it;
@@ -630,7 +630,7 @@ struct Incident_facet_circulator<3, Triangulator>
// Another DS for search space, which can be used in triangulate_DT
// It is useful for extending the search space of 3D Triangulation by appending new triangles
struct Edge_graph
struct Edge_graph
{
struct Edge_comp {
bool operator()(std::pair<int, int> p0, std::pair<int, int> p1) const {
@@ -683,7 +683,7 @@ struct Edge_graph
typedef typename Triangulation::Finite_edges_iterator Finite_edges_iterator;
n = static_cast<int>(edge_exist.size());
for(Finite_edges_iterator eb = tr.finite_edges_begin(); eb != tr.finite_edges_end(); ++eb)
for(Finite_edges_iterator eb = tr.finite_edges_begin(); eb != tr.finite_edges_end(); ++eb)
{
int v0 = eb->first->vertex(eb->second)->info();
int v1 = eb->first->vertex(eb->third )->info();
@@ -729,14 +729,14 @@ class Triangulate_hole_polyline;
#ifndef CGAL_HOLE_FILLING_DO_NOT_USE_DT3
// By default Lookup_table_map is used, since Lookup_table requires n*n mem.
// Performance decrease is nearly 2x (for n = 10,000, for larger n Lookup_table just goes out of mem)
// Performance decrease is nearly 2x (for n = 10,000, for larger n Lookup_table just goes out of mem)
template<
class Kernel,
class Tracer,
class WeightCalculator,
template <class> class LookupTable = Lookup_table_map
>
class Triangulate_hole_polyline_DT
class Triangulate_hole_polyline_DT
{
struct Auto_count {
typedef std::pair<typename Kernel::Point_3, int> result_type;
@@ -767,7 +767,7 @@ public:
typedef Incident_facet_circulator<2, Triangulate_hole_polyline_DT> IFC_2;
typedef Incident_facet_circulator<3, Triangulate_hole_polyline_DT> IFC_3;
Weight operator()(const Polyline_3& P,
Weight operator()(const Polyline_3& P,
const Polyline_3& Q,
Tracer& tracer,
const WeightCalculator& WC) const
@@ -796,7 +796,7 @@ public:
LookupTable<int> lambda(n,-1);
typename Incident_facet_circulator_base<Triangulate_hole_polyline_DT>::Edge_wrapper
e_start(*boost::get<0>(res));
e_start(*boost::get<0>(res));
if(tr.dimension() == 3) {
triangulate_DT<IFC_3>(P, Q, W, lambda, e_start, tr, WC, false);
}
@@ -817,12 +817,12 @@ public:
tracer(lambda, 0, n-1);
return W.get(0,n-1);
}
// How to handle missing border edges
#if 1
return fill_by_extra_triangles(tr, edge_exist, P, Q, tracer, WC);
#else
// This approach produce better patches when used with Weight_incomplete
// This approach produce better patches when used with Weight_incomplete
// (which should be arranged in internal::triangulate_hole_Polyhedron, triangulate_polyline)
return fill_by_incomplete_patches(tr, res.get<0>(), edge_exist, P, Q, tracer, WC);
#endif
@@ -842,10 +842,10 @@ private:
* - edge_first() and edge_second() neighbor edges to get_third() vertex
************************************************************************/
template<class IncidentFacetCirculator, class Edge_DT, class Triangulation_DT>
void triangulate_DT(const Polyline_3& P,
const Polyline_3& Q,
LookupTable<Weight>& W,
LookupTable<int>& lambda,
void triangulate_DT(const Polyline_3& P,
const Polyline_3& Q,
LookupTable<Weight>& W,
LookupTable<int>& lambda,
Edge_DT e,
const Triangulation_DT& tr,
const WeightCalculator& WC,
@@ -877,7 +877,7 @@ private:
if(WC(P,Q, v0,v2,v1, lambda) == Weight::NOT_VALID())
{ continue; } // computed weight in here is not correct weight
// since max dih angle requires neighbor ranges to be already computed. It is just for checking validity.
Weight w = Weight::DEFAULT();
Edge_DT e0 = fb.edge_first(); // edge v0-v2
@@ -885,7 +885,7 @@ private:
triangulate_DT<IncidentFacetCirculator>(P, Q, W, lambda, e0, tr, WC, produce_incomplete); // region v0-v2
const Weight& we0 = W.get(v0, v2);
if(!produce_incomplete && we0 == Weight::NOT_VALID())
if(!produce_incomplete && we0 == Weight::NOT_VALID())
{ continue; } // not producing incomplete patches and failed to fill sub-range v0-v2, so no reason to proceed
if(we0 != Weight::NOT_VALID()) // to not consider we0 if it is NOT_VALID (it is valid when produce_incomplete = true)
{ w = w + we0; }
@@ -895,7 +895,7 @@ private:
triangulate_DT<IncidentFacetCirculator>(P, Q, W, lambda, e1, tr, WC, produce_incomplete); // region v2-v1
const Weight& we1 = W.get(v2, v1);
if(!produce_incomplete && we1 == Weight::NOT_VALID())
if(!produce_incomplete && we1 == Weight::NOT_VALID())
{ continue; }
if(we1 != Weight::NOT_VALID()) // to not consider we1 if it is NOT_VALID (it is valid when produce_incomplete = true)
{ w = w + we1; }
@@ -917,7 +917,7 @@ private:
construct_3D_triangulation(const Polyline_3& P,
std::pair<int,int> h,
Triangulation& tr,
std::vector<bool>& edge_exist) const
std::vector<bool>& edge_exist) const
{
// construct 3D tr with P[h.first], P[h.second] also assign ids from h.first to h.second
boost::optional<Edge> e;
@@ -934,7 +934,7 @@ private:
Finite_edges_iterator v_first_v_second_edge; // range.first - range.second edge
for(Finite_edges_iterator eb = tr.finite_edges_begin();
eb != tr.finite_edges_end();
++eb)
++eb)
{
int v0_id = eb->first->vertex(eb->second)->info();
int v1_id = eb->first->vertex(eb->third )->info();
@@ -964,13 +964,13 @@ private:
* Try to construct hole part by part.
*
* What need to be improved:
* + if 3D triangulation does not contain the start-edge (edge between v0 vn-1) we directly switch to all space.
* + if 3D triangulation does not contain the start-edge (edge between v0 vn-1) we directly switch to all space.
* + when switched to all-space, we use map based lookup tables.
************************************************************************/
Weight fill_by_incomplete_patches(Triangulation& tr,
boost::optional<Edge> start_edge,
std::vector<bool>& edge_exist,
const Polyline_3& P,
const Polyline_3& P,
const Polyline_3& Q,
Tracer& tracer,
const WeightCalculator& WC) const
@@ -978,11 +978,11 @@ private:
typedef std::pair<int, int> Range;
typedef std::back_insert_iterator<std::vector<Range> > Output_hole_iterator;
typedef Tracer_polyline_incomplete<boost::tuple<int, int, int>, Emptyset_iterator, Output_hole_iterator> Remaining_holes_tracer;
std::vector<Range> remaining_holes;
int n_all = P.size()-1;// because the first point and last point are equal
remaining_holes.push_back(Range(0, n_all-1)); // corresponds to start_edge
remaining_holes.push_back(Range(0, n_all-1)); // corresponds to start_edge
LookupTable<Weight> W(n_all, Weight::DEFAULT()); // do not forget that these default values are not changed for [i, i+1]
LookupTable<int> lambda(n_all,-1);
@@ -990,15 +990,15 @@ private:
while(true) {
Range h = remaining_holes.back();
remaining_holes.pop_back();
if(start_edge) {
typename IFC_3::Edge_wrapper e(*start_edge);
CGAL_assertion(h.first == e.vertex_first() &&
h.second == e.vertex_second());
CGAL_assertion(h.first == e.vertex_first() &&
h.second == e.vertex_second());
}
if(!start_edge) {
// switch to brute force
// switch to brute force
Triangulate_hole_polyline<Kernel, Tracer, WeightCalculator, LookupTable> all_space;
all_space.triangulate_all(P, Q, WC, std::make_pair(h.first, h.second), W, lambda);
if(W.get(h.first, h.second) == Weight::NOT_VALID()) {
@@ -1009,7 +1009,7 @@ private:
else {
// run the algorithm
typename IFC_3::Edge_wrapper e(*start_edge);
if(tr.dimension() == 3)
if(tr.dimension() == 3)
{
triangulate_DT<IFC_3>(P, Q, W, lambda, e, tr, WC, true);
}
@@ -1019,7 +1019,7 @@ private:
}
// check whether there is any improvement (at least we should construct one triangle)
if(W.get(h.first, h.second) == Weight::NOT_VALID()) {
// switch to brute force
// switch to brute force
Triangulate_hole_polyline<Kernel, Tracer, WeightCalculator, LookupTable> all_space;
all_space.triangulate_all(P, Q, WC, std::make_pair(h.first, h.second), W, lambda);
if(W.get(h.first, h.second) == Weight::NOT_VALID()) {
@@ -1048,13 +1048,13 @@ private:
}
tracer(lambda, 0, n_all-1);
// W.get(0, n_all -1) is not correct weight (since we do not update weights while we are filling remaining holes),
// W.get(0, n_all -1) is not correct weight (since we do not update weights while we are filling remaining holes),
// we need to recalculate it
std::stack<std::pair<int, int> > ranges;
ranges.push(std::make_pair(0, n_all-1));
Weight total_weight = Weight::DEFAULT();
while(!ranges.empty()) {
std::pair<int, int> r = ranges.top();
std::pair<int, int> r = ranges.top();
ranges.pop();
if(r.first + 1 == r.second) { continue; }
int la = lambda.get(r.first, r.second);
@@ -1064,20 +1064,20 @@ private:
}
return total_weight;
}
/************************************************************************
* This approach extends the search space by adding extra triangles.
*
* + Initial search space is 3D Triangulation.
* + For each border edge which is not inside 3DT, we add all possible triangles containing that edge to the search space.
* Example: say border edge [4-5] is not found in 3DT, then triangles = { [0,4,5] [1,4,5] [2,4,5] ... [ n-1,4,5] }
* Example: say border edge [4-5] is not found in 3DT, then triangles = { [0,4,5] [1,4,5] [2,4,5] ... [ n-1,4,5] }
* are added to the search space.
*
* I guess this approach does not make the search space complete, since there are some cases that it still returns no patch.
************************************************************************/
Weight fill_by_extra_triangles(const Triangulation& tr,
Weight fill_by_extra_triangles(const Triangulation& tr,
const std::vector<bool>& edge_exist,
const Polyline_3& P,
const Polyline_3& P,
const Polyline_3& Q,
Tracer& tracer,
const WeightCalculator& WC) const
@@ -1088,15 +1088,15 @@ private:
Edge_graph edge_graph;
if(tr.dimension() == 3)
if(tr.dimension() == 3)
{ edge_graph.init<IFC_3>(tr, edge_exist); }
else
else
{ edge_graph.init<IFC_2>(tr, edge_exist); }
Edge_graph::Edge_wrapper e_start(std::make_pair(0, n-1));
triangulate_DT<Edge_graph::Incident_facet_circulator>
(P, Q, W, lambda, e_start, edge_graph, WC, false);
if(W.get(0, n-1) == Weight::NOT_VALID()) {
#ifndef CGAL_TEST_SUITE
CGAL_warning_msg(false, "Returning no output using Delaunay triangulation.\n Falling back to the general Triangulation framework.");
@@ -1136,11 +1136,11 @@ public:
CGAL_assertion(P.front() == P.back());
CGAL_assertion(Q.empty() || (Q.front() == Q.back()));
CGAL_assertion(Q.empty() || (P.size() == Q.size()));
int n = static_cast<int>(P.size()) - 1; // because the first and last point are equal
LookupTable<Weight> W(n,Weight::DEFAULT()); // do not forget that these default values are not changed for [i, i+1]
LookupTable<int> lambda(n,-1);
triangulate_all(P, Q, WC, std::make_pair(0,n-1), W, lambda);
if(W.get(0,n-1) == Weight::NOT_VALID() || n <= 2) {
@@ -1160,23 +1160,23 @@ public:
const Polyline_3& Q,
const WeightCalculator& WC,
std::pair<int, int> range,
LookupTable<Weight>& W,
LookupTable<Weight>& W,
LookupTable<int>& lambda) const
{
for(int j = 2; j<= range.second; ++j) { // determines range (2 - 3 - 4 )
for(int i=range.first; i<= range.second-j; ++i) { // iterates over ranges and find min triangulation in those ranges
for(int i=range.first; i<= range.second-j; ++i) { // iterates over ranges and find min triangulation in those ranges
int k = i+j; // like [0-2, 1-3, 2-4, ...], [0-3, 1-4, 2-5, ...]
int m_min = -1;
Weight w_min = Weight::NOT_VALID();
// i is the range start (e.g. 1) k is the range end (e.g. 5) -> [1-5]. Now subdivide the region [1-5] with m -> 2,3,4
for(int m = i+1; m<k; ++m) {
for(int m = i+1; m<k; ++m) {
// now the regions i-m and m-k might be valid(constructed) patches,
if( W.get(i,m) == Weight::NOT_VALID() || W.get(m,k) == Weight::NOT_VALID() )
if( W.get(i,m) == Weight::NOT_VALID() || W.get(m,k) == Weight::NOT_VALID() )
{ continue; }
const Weight& w_imk = WC(P,Q,i,m,k, lambda);
if(w_imk == Weight::NOT_VALID())
if(w_imk == Weight::NOT_VALID())
{ continue; }
const Weight& w = W.get(i,m) + W.get(m,k) + w_imk;
@@ -6,7 +6,7 @@
// $URL$
// $Id$
// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
//
//
//
// Author(s) : Jane Tournois
@@ -6,7 +6,7 @@
// $URL$
// $Id$
// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
//
//
//
// Author(s) : Ilker O. Yaz
@@ -17,10 +17,10 @@
// It can produce a patch from both complete and incomplete lambda
// WARNING: Not working good for all cases
// For holes, this code first close them then erase them.
// For holes, this code first close them then erase them.
// However the algorithm might produce holes which are invalid to close (closing them breaks edge manifoldness, so erasing doesn't work)
template<class Polyhedron, class OutputIteratorPatch, class OutputIteratorHole>
struct Tracer_polyhedron_incomplete
struct Tracer_polyhedron_incomplete
{
typedef typename Polyhedron::Halfedge_handle Halfedge_handle;
typedef typename Polyhedron::Facet_handle Facet_handle;
@@ -33,12 +33,12 @@ struct Tracer_polyhedron_incomplete
{ }
template <class LookupTable>
void
void
operator()(const LookupTable& lambda, int i, int k)
{
std::vector<Facet_handle> facets_to_delete;
(*this)(lambda, i, k, facets_to_delete, true);
for(typename std::vector<Facet_handle>::iterator it = facets_to_delete.begin();
for(typename std::vector<Facet_handle>::iterator it = facets_to_delete.begin();
it != facets_to_delete.end(); ++it)
{
*out_hole++=(*it)->halfedge(); // each deleted facet corresponds to a new hole
@@ -48,8 +48,8 @@ struct Tracer_polyhedron_incomplete
private:
template <class LookupTable>
Halfedge_handle
operator()(const LookupTable& lambda,
Halfedge_handle
operator()(const LookupTable& lambda,
int i, int k,
std::vector<Facet_handle>& facets_to_delete,
bool last)
@@ -61,7 +61,7 @@ private:
if(last)
{ h = polyhedron.fill_hole(P[i+1]); }
else
else
{ h = polyhedron.add_facet_to_border(P[i+1]->prev(), P[i+2/*k*/]); }
CGAL_assertion(h->facet() != Facet_handle());
@@ -74,14 +74,14 @@ private:
*out++ = h->facet();
}
return h->opposite();
}
else
}
else
{
int la = lambda.get(i, k);
if(la == -1) {
if(last)
{ h = polyhedron.fill_hole(P[i+1]); }
else
else
{ h = polyhedron.add_facet_to_border(P[i+1]->prev(), P[i+2/*k*/]); }
facets_to_delete.push_back(h->facet());
return h->opposite();
@@ -92,7 +92,7 @@ private:
if(last)
{ h = polyhedron.fill_hole(g); }
else
else
{ h = polyhedron.add_facet_to_border(h->prev(), g); }
CGAL_assertion(h->facet() != Facet_handle());
@@ -112,11 +112,11 @@ public:
// Try closing holes by gathering incomplete patches together (an external approach)
template <typename OutputIteratorValueType, typename InputIterator, typename OutputIterator>
OutputIterator
triangulate_hole_polyline_incomplete(InputIterator pbegin, InputIterator pend,
InputIterator qbegin, InputIterator qend,
triangulate_hole_polyline_incomplete(InputIterator pbegin, InputIterator pend,
InputIterator qbegin, InputIterator qend,
OutputIterator out)
{
typedef typename std::iterator_traits<InputIterator>::value_type Point_3;
typedef Weight_incomplete<Weight_min_max_dihedral_and_area> Weight;
typedef Weight_calculator<Weight, Is_valid_degenerate_triangle> WC;
@@ -125,7 +125,7 @@ triangulate_hole_polyline_incomplete(InputIterator pbegin, InputIterator pend,
typedef std::back_insert_iterator<Facet_vector> OutIt;
typedef Tracer_polyline_incomplete<Facet_vector::value_type, OutIt> Tracer;
typedef std::pair<int, int> Range;
std::vector<Point_3> P(pbegin, pend);
std::vector<Point_3> Q(qbegin, qend);
@@ -135,7 +135,7 @@ triangulate_hole_polyline_incomplete(InputIterator pbegin, InputIterator pend,
Q.push_back(Q.front());
}
}
std::vector<OutputIteratorValueType> patch_facets;
std::stack<Range> remaining_holes;
@@ -184,8 +184,8 @@ triangulate_hole_polyline_incomplete(InputIterator pbegin, InputIterator pend,
// (for Polyhedron_3) Try closing holes by gathering incomplete patches together (an external approach)
template<class Polyhedron, class OutputIterator>
std::pair<OutputIterator, Weight_min_max_dihedral_and_area>
triangulate_hole_Polyhedron_incomplete(Polyhedron& polyhedron,
typename Polyhedron::Halfedge_handle border_halfedge,
triangulate_hole_Polyhedron_incomplete(Polyhedron& polyhedron,
typename Polyhedron::Halfedge_handle border_halfedge,
OutputIterator out)
{
typedef typename Polyhedron::Halfedge_handle Halfedge_handle;
@@ -349,7 +349,7 @@ namespace internal {
Patch_id pid = get_patch_id(f);
input_triangles_.push_back(triangle(f));
input_patch_ids_.push_back(pid);
std::pair<typename Patch_id_to_index_map::iterator, bool>
std::pair<typename Patch_id_to_index_map::iterator, bool>
res = patch_id_to_index_map.insert(std::make_pair(pid,0));
if(res.second){
res.first->second = patch_id_to_index_map.size()-1;
@@ -1361,7 +1361,7 @@ private:
do
{
if (is_on_patch_border(nxt))
{
{
CGAL_assertion(get_patch_id(face(nxt, mesh_)) == pid);
return nxt;
}
@@ -1513,7 +1513,9 @@ private:
// tag patch border halfedges
for(halfedge_descriptor h : halfedges(mesh_))
{
if (status(h)==PATCH && status(opposite(h, mesh_))!=PATCH)
if (status(h) == PATCH
&& ( status(opposite(h, mesh_)) != PATCH
|| get_patch_id(face(h, mesh_)) != get_patch_id(face(opposite(h, mesh_), mesh_))))
{
set_status(h, PATCH_BORDER);
has_border_ = true;
@@ -170,7 +170,7 @@ public:
// calls compute once to factorize with the preconditioner
if(!solver.factor(A, D))
{
#ifdef CGAL_PMP_SMOOTHING_VERBOSE
#ifdef CGAL_PMP_SMOOTHING_DEBUG
std::cerr << "Could not factorize linear system with preconditioner." << std::endl;
#endif
return false;
@@ -180,7 +180,7 @@ public:
!solver.linear_solver(by, Xy) ||
!solver.linear_solver(bz, Xz))
{
#ifdef CGAL_PMP_SMOOTHING_VERBOSE
#ifdef CGAL_PMP_SMOOTHING_DEBUG
std::cerr << "Could not solve linear system." << std::endl;
#endif
return false;
@@ -21,7 +21,7 @@
#endif
#include <CGAL/Polygon_mesh_processing/compute_normal.h>
#include <CGAL/Polygon_mesh_processing/repair.h>
#include <CGAL/Polygon_mesh_processing/shape_predicates.h>
#include <CGAL/AABB_tree.h>
#include <CGAL/AABB_traits.h>
@@ -50,9 +50,11 @@ namespace Polygon_mesh_processing {
namespace internal {
template <typename V, typename GT>
double get_radian_angle(const V& v1, const V& v2, const GT& gt)
typename GT::FT get_radian_angle(const V& v1, const V& v2, const GT& gt)
{
return gt.compute_approximate_angle_3_object()(v1, v2) * CGAL_PI / 180.;
typedef typename GT::FT FT;
return gt.compute_approximate_angle_3_object()(v1, v2) * CGAL_PI / FT(180);
}
// super naive for now. Not sure it even makes sense to do something like that for surfaces
@@ -68,6 +70,7 @@ class Delaunay_edge_flipper
typedef typename boost::graph_traits<TriangleMesh>::face_descriptor face_descriptor;
typedef typename boost::property_traits<VertexPointMap>::reference Point_ref;
typedef typename GeomTraits::FT FT;
typedef typename GeomTraits::Vector_3 Vector;
public:
@@ -86,20 +89,13 @@ public:
const halfedge_descriptor h = halfedge(e, mesh_);
const halfedge_descriptor opp_h = opposite(h, mesh_);
vertex_descriptor v0 = source(h, mesh_);
vertex_descriptor v1 = target(h, mesh_);
vertex_descriptor v2 = target(next(h, mesh_), mesh_);
vertex_descriptor v3 = target(next(opp_h, mesh_), mesh_);
const Point_ref p0 = get(vpmap_, v0);
const Point_ref p1 = get(vpmap_, v1);
const Point_ref p2 = get(vpmap_, v2);
const Point_ref p3 = get(vpmap_, v3);
const vertex_descriptor v0 = source(h, mesh_);
const vertex_descriptor v1 = target(h, mesh_);
const vertex_descriptor v2 = target(next(h, mesh_), mesh_);
const vertex_descriptor v3 = target(next(opp_h, mesh_), mesh_);
double alpha = get_radian_angle(Vector(p0 - p2), Vector(p1 - p2), traits_);
double beta = get_radian_angle(Vector(p1 - p3), Vector(p0 - p3), traits_);
// not local Delaunay if the sum of the angles is greater than pi
if(alpha + beta <= CGAL_PI)
std::set<vertex_descriptor> unique_vs { v0, v1, v2, v3 };
if(unique_vs.size() != 4)
return false;
// Don't want to flip if the other diagonal already exists
@@ -108,7 +104,16 @@ public:
if(other_hd_already_exists.second)
return false;
return true;
// not local Delaunay := sum of the opposite angles is greater than pi
const Point_ref p0 = get(vpmap_, v0);
const Point_ref p1 = get(vpmap_, v1);
const Point_ref p2 = get(vpmap_, v2);
const Point_ref p3 = get(vpmap_, v3);
FT alpha = get_radian_angle(Vector(p0 - p2), Vector(p1 - p2), traits_);
FT beta = get_radian_angle(Vector(p1 - p3), Vector(p0 - p3), traits_);
return (alpha + beta > CGAL_PI);
}
template <typename Marked_edges_map, typename EdgeRange>
@@ -127,6 +132,10 @@ public:
template <typename FaceRange>
void operator()(const FaceRange& face_range)
{
#ifdef CGAL_PMP_SMOOTHING_DEBUG
std::cout << "Flipping edges" << std::endl;
#endif
// edges to consider
std::vector<edge_descriptor> edge_range;
edge_range.reserve(3 * face_range.size());
@@ -166,6 +175,10 @@ public:
++flipped_n;
halfedge_descriptor h = halfedge(e, mesh_);
#ifdef CGAL_PMP_SMOOTHING_DEBUG_PP
std::cout << "Flipping " << edge(h, mesh_) << std::endl;
#endif
Euler::flip_edge(h, mesh_);
add_to_stack_if_unmarked(edge(next(h, mesh_), mesh_), marks, edge_range);
@@ -174,6 +187,10 @@ public:
add_to_stack_if_unmarked(edge(prev(opposite(h, mesh_), mesh_), mesh_), marks, edge_range);
}
}
#ifdef CGAL_PMP_SMOOTHING_DEBUG
std::cout << flipped_n << " flips" << std::endl;
#endif
}
private:
@@ -190,6 +207,7 @@ class Angle_smoother
typedef typename boost::graph_traits<TriangleMesh>::halfedge_descriptor halfedge_descriptor;
typedef typename boost::property_traits<VertexPointMap>::reference Point_ref;
typedef typename GeomTraits::FT FT;
typedef typename GeomTraits::Vector_3 Vector;
typedef std::pair<halfedge_descriptor, halfedge_descriptor> He_pair;
@@ -231,7 +249,7 @@ public:
Vector operator()(const vertex_descriptor v) const
{
Vector move = CGAL::NULL_VECTOR;
double weights_sum = 0.;
FT weights_sum = FT(0);
for(halfedge_descriptor main_he : halfedges_around_source(v, mesh_))
{
@@ -249,26 +267,30 @@ public:
Vector right_v(pt, right_pt);
// rotate
double angle = get_radian_angle(right_v, left_v, traits_);
CGAL_warning(angle != 0.); // no degenerate faces is a precondition
if(angle == 0.)
continue;
Vector bisector = rotate_edge(main_he, incident_pair);
double scaling_factor = CGAL::approximate_sqrt(
traits_.compute_squared_distance_3_object()(get(vpmap_, source(main_he, mesh_)),
get(vpmap_, target(main_he, mesh_))));
FT scaling_factor = CGAL::approximate_sqrt(
traits_.compute_squared_distance_3_object()(get(vpmap_, source(main_he, mesh_)),
get(vpmap_, target(main_he, mesh_))));
bisector = traits_.construct_scaled_vector_3_object()(bisector, scaling_factor);
Vector ps_psi(ps, traits_.construct_translated_point_3_object()(pt, bisector));
FT angle = get_radian_angle(right_v, left_v, traits_);
if(angle == FT(0))
{
// no degenerate faces is a precondition, angle can be 0 but it should be a numerical error
CGAL_warning(!is_degenerate_triangle_face(face(main_he, mesh_), mesh_));
return ps_psi; // since a small angle gives more weight, a null angle give priority (?)
}
// small angles carry more weight
double weight = 1. / (angle*angle);
FT weight = 1. / CGAL::square(angle);
weights_sum += weight;
move += weight * ps_psi;
}
if(weights_sum != 0.)
if(weights_sum != FT(0))
move /= weights_sum;
return move;
@@ -288,6 +310,7 @@ class Area_smoother
typedef typename boost::property_traits<VertexPointMap>::value_type Point;
typedef typename boost::property_traits<VertexPointMap>::reference Point_ref;
typedef typename GeomTraits::FT FT;
typedef typename GeomTraits::Vector_3 Vector;
public:
@@ -298,27 +321,27 @@ public:
{ }
private:
double element_area(const vertex_descriptor v1,
const vertex_descriptor v2,
const vertex_descriptor v3) const
FT element_area(const vertex_descriptor v1,
const vertex_descriptor v2,
const vertex_descriptor v3) const
{
return CGAL::to_double(CGAL::approximate_sqrt(traits_.compute_squared_area_3_object()(get(vpmap_, v1),
get(vpmap_, v2),
get(vpmap_, v3))));
return CGAL::approximate_sqrt(traits_.compute_squared_area_3_object()(get(vpmap_, v1),
get(vpmap_, v2),
get(vpmap_, v3)));
}
double element_area(const Point& P,
const vertex_descriptor v2,
const vertex_descriptor v3) const
FT element_area(const Point& P,
const vertex_descriptor v2,
const vertex_descriptor v3) const
{
return CGAL::to_double(CGAL::approximate_sqrt(traits_.compute_squared_area_3_object()(P,
get(vpmap_, v2),
get(vpmap_, v3))));
return CGAL::approximate_sqrt(traits_.compute_squared_area_3_object()(P,
get(vpmap_, v2),
get(vpmap_, v3)));
}
double compute_average_area_around(const vertex_descriptor v) const
FT compute_average_area_around(const vertex_descriptor v) const
{
double sum_areas = 0.;
FT sum_areas = 0;
unsigned int number_of_edges = 0;
for(halfedge_descriptor h : halfedges_around_source(v, mesh_))
@@ -327,7 +350,7 @@ private:
vertex_descriptor vi = source(next(h, mesh_), mesh_);
vertex_descriptor vj = target(next(h, mesh_), mesh_);
double S = element_area(v, vi, vj);
FT S = element_area(v, vi, vj);
sum_areas += S;
++number_of_edges;
}
@@ -337,7 +360,7 @@ private:
struct Face_energy
{
Face_energy(const Point& pi, const Point& pj, const double s_av)
Face_energy(const Point& pi, const Point& pj, const FT s_av)
:
qx(pi.x()), qy(pi.y()), qz(pi.z()),
rx(pj.x()), ry(pj.y()), rz(pj.z()),
@@ -346,7 +369,7 @@ private:
// next two functions are just for convencience, the only thing ceres cares about is the operator()
template <typename T>
double area(const T x, const T y, const T z) const
FT area(const T x, const T y, const T z) const
{
return CGAL::approximate_sqrt(CGAL::squared_area(Point(x, y, z),
Point(qx, qy, qz),
@@ -354,7 +377,7 @@ private:
}
template <typename T>
double evaluate(const T x, const T y, const T z) const { return area(x, y, z) - s_av; }
FT evaluate(const T x, const T y, const T z) const { return area(x, y, z) - s_av; }
template <typename T>
bool operator()(const T* const x, const T* const y, const T* const z,
@@ -390,9 +413,9 @@ private:
}
private:
const double qx, qy, qz;
const double rx, ry, rz;
const double s_av;
const FT qx, qy, qz;
const FT rx, ry, rz;
const FT s_av;
};
public:
@@ -401,12 +424,12 @@ public:
#ifdef CGAL_PMP_USE_CERES_SOLVER
const Point_ref vp = get(vpmap_, v);
const double S_av = compute_average_area_around(v);
const FT S_av = compute_average_area_around(v);
const double initial_x = vp.x();
const double initial_y = vp.y();
const double initial_z = vp.z();
double x = initial_x, y = initial_y, z = initial_z;
const FT initial_x = vp.x();
const FT initial_y = vp.y();
const FT initial_z = vp.z();
FT x = initial_x, y = initial_y, z = initial_z;
ceres::Problem problem;
@@ -523,7 +546,8 @@ public:
Optimizer compute_move(mesh_, vpmap_, traits_);
#ifdef CGAL_PMP_SMOOTHING_DEBUG
double total_displacement = 0;
FT total_displacement = 0;
std::cout << "apply_moves_in_single_batch: " << apply_moves_in_single_batch << std::endl;
#endif
std::size_t moved_points = 0;
@@ -532,6 +556,10 @@ public:
if(is_border(v, mesh_) || is_constrained(v))
continue;
#ifdef CGAL_PMP_SMOOTHING_DEBUG_PP
std::cout << "Considering " << v << " pos: " << get(vpmap_, v) << std::endl;
#endif
// compute normal to v
Vector vn = compute_vertex_normal(v, mesh_, CGAL::parameters::vertex_point_map(vpmap_)
.geom_traits(traits_));
@@ -542,17 +570,17 @@ public:
// Gram Schmidt so that the new location is on the tangent plane of v (i.e. do mv -= (mv*n)*n)
const FT sp = traits_.compute_scalar_product_3_object()(vn, move);
move = traits_.construct_sum_of_vectors_3_object()(move,
traits_.construct_scaled_vector_3_object()(vn, - sp));
move = traits_.construct_sum_of_vectors_3_object()(
move, traits_.construct_scaled_vector_3_object()(vn, - sp));
Point new_pos = pos + move;
if((!use_sanity_checks || !does_move_create_bad_faces(v, new_pos)) &&
const Point new_pos = pos + move;
if(move != CGAL::NULL_VECTOR &&
!does_move_create_degenerate_faces(v, new_pos) &&
(!use_sanity_checks || !does_move_create_bad_faces(v, new_pos)) &&
(!enforce_no_min_angle_regression || does_improve_min_angle_in_star(v, new_pos)))
{
#ifdef CGAL_PMP_SMOOTHING_DEBUG
#ifdef CGAL_PMP_SMOOTHING_DEBUG_PP
std::cout << "moving " << get(vpmap_, v) << " to " << new_pos << std::endl;
total_displacement += CGAL::approximate_sqrt(traits_.compute_squared_length_3_object()(move));
#endif
if(apply_moves_in_single_batch)
@@ -560,11 +588,15 @@ public:
else
put(vpmap_, v, new_pos);
#ifdef CGAL_PMP_SMOOTHING_DEBUG
total_displacement += CGAL::approximate_sqrt(traits_.compute_squared_length_3_object()(move));
#endif
++moved_points;
}
else // some sanity check failed
{
#ifdef CGAL_PMP_SMOOTHING_DEBUG
#ifdef CGAL_PMP_SMOOTHING_DEBUG_PP
std::cout << "move is rejected!" << std::endl;
#endif
if(apply_moves_in_single_batch)
@@ -607,6 +639,10 @@ public:
Point_ref p_query = get(vpmap_, v);
const Point projected = tree.closest_point(p_query);
#ifdef CGAL_PMP_SMOOTHING_DEBUG_PP
std::cout << p_query << " to " << projected << std::endl;
#endif
put(vpmap_, v, projected);
}
}
@@ -617,11 +653,10 @@ private:
return get(vcmap_, v);
}
// check for degenerate or inversed faces
bool does_move_create_bad_faces(const vertex_descriptor v,
const Point& new_pos) const
// Null faces are bad because they make normal computation difficult
bool does_move_create_degenerate_faces(const vertex_descriptor v,
const Point& new_pos) const
{
// check for null faces and face inversions
for(halfedge_descriptor main_he : halfedges_around_source(v, mesh_))
{
const halfedge_descriptor prev_he = prev(main_he, mesh_);
@@ -630,6 +665,23 @@ private:
if(traits_.collinear_3_object()(lpt, rpt, new_pos))
return true;
}
return false;
}
// check for degenerate or inversed faces
bool does_move_create_bad_faces(const vertex_descriptor v,
const Point& new_pos) const
{
// check for face inversions
for(halfedge_descriptor main_he : halfedges_around_source(v, mesh_))
{
const halfedge_descriptor prev_he = prev(main_he, mesh_);
const Point_ref lpt = get(vpmap_, target(main_he, mesh_));
const Point_ref rpt = get(vpmap_, source(prev_he, mesh_));
CGAL_assertion(!traits_.collinear_3_object()(lpt, rpt, new_pos)); // checked above
const Point_ref old_pos = get(vpmap_, v);
Vector ov_1 = traits_.construct_vector_3_object()(old_pos, lpt);
@@ -641,7 +693,7 @@ private:
if(!is_positive(traits_.compute_scalar_product_3_object()(old_n, new_n)))
{
#ifdef CGAL_PMP_SMOOTHING_DEBUG
#ifdef CGAL_PMP_SMOOTHING_DEBUG_PP
std::cout << "Moving vertex would result in the inversion of a face normal!" << std::endl;
#endif
return true;
@@ -655,7 +707,7 @@ private:
const Point& new_pos) const
{
// check if the minimum angle of the star has not deteriorated
double old_min_angle = CGAL_PI;
FT old_min_angle = CGAL_PI;
for(halfedge_descriptor main_he : halfedges_around_source(v, mesh_))
{
const Point_ref old_pos = get(vpmap_, v);
@@ -683,7 +735,7 @@ private:
get_radian_angle(Vector(lpt, rpt), Vector(lpt, new_pos), traits_) < old_min_angle ||
get_radian_angle(Vector(rpt, new_pos), Vector(rpt, lpt), traits_) < old_min_angle)
{
#ifdef CGAL_PMP_SMOOTHING_DEBUG
#ifdef CGAL_PMP_SMOOTHING_DEBUG_PP
const Point_ref old_pos = get(vpmap_, v);
std::cout << "deterioration of min angle in the star!" << std::endl;
@@ -61,7 +61,7 @@ public:
angles_.push_back(traits_.compute_approximate_angle_3_object()(a, b, c));
}
#ifdef CGAL_PMP_SMOOTHING_VERBOSE
#ifdef CGAL_PMP_SMOOTHING_DEBUG
std::cout << "angles_ size = " << angles_.size() << std::endl;
#endif
}
@@ -79,7 +79,7 @@ public:
for(face_descriptor f : faces(mesh_))
areas_.push_back(face_area(f, mesh_));
#ifdef CGAL_PMP_SMOOTHING_VERBOSE
#ifdef CGAL_PMP_SMOOTHING_DEBUG
std::cout << "areas_ size = " << areas_.size() << std::endl;
#endif
}
@@ -97,7 +97,7 @@ public:
for(face_descriptor f : faces(mesh_))
aspect_ratios_.push_back(CGAL::Polygon_mesh_processing::face_aspect_ratio(f, mesh_));
#ifdef CGAL_PMP_SMOOTHING_VERBOSE
#ifdef CGAL_PMP_SMOOTHING_DEBUG
std::cout << "aspect_ratios_ size = " << aspect_ratios_.size() << std::endl;
#endif
}
@@ -58,7 +58,7 @@ void assign_tolerance_with_local_edge_length_bound(const HalfedgeRange& halfedge
using parameters::get_parameter;
using parameters::choose_parameter;
GT gt = choose_parameter(get_parameter(np, internal_np::geom_traits), GT());
GT gt = choose_parameter<GT>(get_parameter(np, internal_np::geom_traits));
VPM vpm = choose_parameter(get_parameter(np, internal_np::vertex_point),
get_property_map(vertex_point, mesh));
@@ -97,7 +97,7 @@ void simplify_range(HalfedgeRange& halfedge_range,
using parameters::get_parameter;
using parameters::choose_parameter;
const GT gt = choose_parameter(get_parameter(np, internal_np::geom_traits), GT());
const GT gt = choose_parameter<GT>(get_parameter(np, internal_np::geom_traits));
VPM vpm = choose_parameter(get_parameter(np, internal_np::vertex_point), get_property_map(vertex_point, tm));
typedef CGAL::dynamic_halfedge_property_t<bool> Halfedge_bool_tag;
@@ -794,7 +794,7 @@ std::size_t snap_non_conformal_one_way(const HalfedgeRange& halfedge_range_S,
VPMS vpm_S = choose_parameter(get_parameter(snp, internal_np::vertex_point), get_property_map(vertex_point, tm_S));
VPMT vpm_T = choose_parameter(get_parameter(tnp, internal_np::vertex_point), get_property_map(vertex_point, tm_T));
const GT gt = choose_parameter(get_parameter(snp, internal_np::geom_traits), GT());
const GT gt = choose_parameter<GT>(get_parameter(snp, internal_np::geom_traits));
#ifdef CGAL_PMP_SNAP_DEBUG
std::cout << "Gather unique points in source range..." << std::endl;
@@ -1185,13 +1185,13 @@ std::size_t snap_borders(TriangleMesh& tm_A,
{
typedef typename boost::graph_traits<TriangleMesh>::halfedge_descriptor halfedge_descriptor;
std::vector<halfedge_descriptor> border_vertices_1;
border_halfedges(tm_A, std::back_inserter(border_vertices_1));
std::vector<halfedge_descriptor> border_vertices_2;
border_halfedges(tm_B, std::back_inserter(border_vertices_2));
std::vector<halfedge_descriptor> border_vertices_A;
border_halfedges(tm_A, std::back_inserter(border_vertices_A));
std::vector<halfedge_descriptor> border_vertices_B;
border_halfedges(tm_B, std::back_inserter(border_vertices_B));
return internal::snap_non_conformal<ConcurrencyTag>(border_vertices_1, tm_A, tolerance_map_A,
border_vertices_2, tm_B, tolerance_map_B,
return internal::snap_non_conformal<ConcurrencyTag>(border_vertices_A, tm_A, tolerance_map_A,
border_vertices_B, tm_B, tolerance_map_B,
false /*not self snapping*/, np_A, np_B);
}
@@ -1211,19 +1211,19 @@ std::size_t snap_borders(TriangleMesh& tm_A,
typedef CGAL::dynamic_vertex_property_t<FT> Vertex_property_tag;
typedef typename boost::property_map<TriangleMesh, Vertex_property_tag>::type Tolerance_map;
std::vector<halfedge_descriptor> border_vertices_1;
std::vector<halfedge_descriptor> border_vertices_2;
border_halfedges(tm_A, std::back_inserter(border_vertices_1));
border_halfedges(tm_B, std::back_inserter(border_vertices_2));
std::vector<halfedge_descriptor> border_vertices_A;
std::vector<halfedge_descriptor> border_vertices_B;
border_halfedges(tm_A, std::back_inserter(border_vertices_A));
border_halfedges(tm_B, std::back_inserter(border_vertices_B));
const FT tol_mx(std::numeric_limits<double>::max());
const FT tol_mx((std::numeric_limits<double>::max)());
Tolerance_map tolerance_map_A = get(Vertex_property_tag(), tm_A);
internal::assign_tolerance_with_local_edge_length_bound(border_vertices_1, tolerance_map_A, tol_mx, tm_A, np_A);
Tolerance_map tolerance_map_B = get(Vertex_property_tag(), tm_A);
internal::assign_tolerance_with_local_edge_length_bound(border_vertices_2, tolerance_map_B, tol_mx, tm_B, np_B);
internal::assign_tolerance_with_local_edge_length_bound(border_vertices_A, tolerance_map_A, tol_mx, tm_A, np_A);
Tolerance_map tolerance_map_B = get(Vertex_property_tag(), tm_B);
internal::assign_tolerance_with_local_edge_length_bound(border_vertices_B, tolerance_map_B, tol_mx, tm_B, np_B);
return internal::snap_non_conformal<ConcurrencyTag>(border_vertices_1, tm_A, tolerance_map_A,
border_vertices_2, tm_B, tolerance_map_B,
return internal::snap_non_conformal<ConcurrencyTag>(border_vertices_A, tm_A, tolerance_map_A,
border_vertices_B, tm_B, tolerance_map_B,
false /*no self snapping*/, np_A, np_B);
}
@@ -1296,7 +1296,7 @@ std::size_t snap_borders(TriangleMesh& tm,
std::vector<halfedge_descriptor> border_vertices;
border_halfedges(tm, std::back_inserter(border_vertices));
const FT tol_mx(std::numeric_limits<double>::max());
const FT tol_mx((std::numeric_limits<double>::max)());
Tolerance_map tolerance_map = get(Vertex_property_tag(), tm);
internal::assign_tolerance_with_local_edge_length_bound(border_vertices, tolerance_map, tol_mx, tm, np);
@@ -299,7 +299,7 @@ std::size_t snap_vertices_two_way(const HalfedgeRange_A& halfedge_range_A,
CGAL_static_assertion((std::is_same<Point, typename GT::Point_3>::value));
GT gt = choose_parameter(get_parameter(np_A, internal_np::geom_traits), GT());
GT gt = choose_parameter<GT>(get_parameter(np_A, internal_np::geom_traits));
VPM_A vpm_A = choose_parameter(get_parameter(np_A, internal_np::vertex_point),
get_property_map(vertex_point, tm_A));
VPM_B vpm_B = choose_parameter(get_parameter(np_B, internal_np::vertex_point),
@@ -868,7 +868,7 @@ std::size_t snap_vertices(const HalfedgeRange_A& halfedge_range_A,
typedef CGAL::dynamic_vertex_property_t<FT> Vertex_property_tag;
typedef typename boost::property_map<PolygonMesh, Vertex_property_tag>::type Tolerance_map;
const FT max_tol(std::numeric_limits<double>::max());
const FT max_tol((std::numeric_limits<double>::max)());
Tolerance_map tolerance_map_A = get(Vertex_property_tag(), tm_A);
internal::assign_tolerance_with_local_edge_length_bound(halfedge_range_A, tolerance_map_A, max_tol, tm_A, np_A);
@@ -6,7 +6,7 @@
// $URL$
// $Id$
// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
//
//
//
// Author(s) : Sebastien Loriot
@@ -6,7 +6,7 @@
// $URL$
// $Id$
// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
//
//
//
// Author(s) : Jane Tournois
@@ -59,7 +59,7 @@ public:
, weight_calculator(weight_calculator)
, ppmap(vpmap)
{ }
private:
double sum_weight(vertex_descriptor v) {
double weight = 0;
@@ -75,7 +75,7 @@ private:
void compute_row(
vertex_descriptor v,
int row_id, // which row to insert in [ frees stay left-hand side ]
Solver_matrix& matrix,
Solver_matrix& matrix,
double& x, double& y, double& z, // constants transfered to right-hand side
double multiplier,
const std::map<vertex_descriptor, std::size_t>& vertex_id_map,
@@ -87,7 +87,7 @@ private:
int col_id = static_cast<int>(vertex_id_it->second);
matrix.add_coef(row_id, col_id, multiplier);
}
else {
else {
typename boost::property_traits<VertexPointMap>::reference p = get(ppmap, v);
x += multiplier * - to_double(p.x());
y += multiplier * - to_double(p.y());
@@ -119,7 +119,7 @@ public:
int depth = static_cast<int>(fc) + 1;
if(depth < 0 || depth > 3) {
CGAL_warning_msg(false, "Continuity should be between 0 and 2 inclusively!");
return false;
return false;
}
std::set<vertex_descriptor> interior_vertices(boost::begin(vertices),
@@ -170,25 +170,25 @@ public:
// solve
bool is_all_solved = solver.linear_solver(Bx, X) && solver.linear_solver(By, Y) && solver.linear_solver(Bz, Z);
if(!is_all_solved) {
CGAL_warning_msg(false, "linear_solver failed!");
return false;
CGAL_warning_msg(false, "linear_solver failed!");
return false;
}
#ifdef CGAL_PMP_FAIR_DEBUG
std::cerr << "**Timer** System solver: " << timer.time() << std::endl; timer.reset();
#endif
/* This relative error is to large for cases that the results are not good */
/* This relative error is to large for cases that the results are not good */
/*
double rel_err_x = (A.eigen_object()*X - Bx).norm() / Bx.norm();
double rel_err_y = (A.eigen_object()*Y - By).norm() / By.norm();
double rel_err_z = (A.eigen_object()*Z - Bz).norm() / Bz.norm();
CGAL_TRACE_STREAM << "rel error: " << rel_err_x
CGAL_TRACE_STREAM << "rel error: " << rel_err_x
<< " " << rel_err_y
<< " " << rel_err_z << std::endl;
*/
// update
// update
id = 0;
for(vertex_descriptor vd : interior_vertices)
{
@@ -86,7 +86,7 @@ public:
m_point[2] = c;
m_edge = 0;
typename Kernel::Compute_squared_distance_3 squared_distance;
FT length_max = squared_distance(m_point[1](), m_point[2]());
FT length1 = squared_distance(m_point[2](), m_point[0]());
@@ -6,7 +6,7 @@
// $URL$
// $Id$
// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
//
//
//
// Author(s) : Jane Tournois
@@ -52,7 +52,7 @@ class Refine_Polyhedron_3 {
private:
PolygonMesh& pmesh;
VertexPointMap vpmap;
bool flippable(halfedge_descriptor h) {
// this check is added so that edge flip does not break manifoldness
// it might happen when there is an edge where flip_edge(h) will be placed (i.e. two edges collide after flip)
@@ -62,7 +62,7 @@ private:
do {
if(target(opposite(*v_cir, pmesh),pmesh) == v_tip_1) { return false; }
} while(++v_cir != v_end);
// also eliminate collinear triangle generation
if( CGAL::collinear(get(vpmap, v_tip_0), get(vpmap, v_tip_1), get(vpmap, target(h, pmesh))) ||
CGAL::collinear(get(vpmap, v_tip_0), get(vpmap, v_tip_1), get(vpmap, target(opposite(h, pmesh),pmesh))) ) {
@@ -95,7 +95,7 @@ private:
class FaceRange>
bool subdivide(const FaceRange& faces,
const std::set<halfedge_descriptor>& border_edges,
std::map<vertex_descriptor, double>& scale_attribute,
std::map<vertex_descriptor, double>& scale_attribute,
VertexOutputIterator& vertex_out,
FaceOutputIterator& facet_out,
std::vector<face_descriptor>& new_faces,
@@ -125,12 +125,12 @@ private:
scale_attribute[target(h,pmesh)] = sac;
*vertex_out++ = target(h,pmesh);
// collect 2 new facets for next round
// collect 2 new facets for next round
face_descriptor h1 = face(opposite(next(h,pmesh),pmesh),pmesh);
face_descriptor h2 = face(opposite(h,pmesh),pmesh);
new_faces.push_back(h1); new_faces.push_back(h2);
*facet_out++ = h1; *facet_out++ = h2;
// relax edges of the patching mesh
// relax edges of the patching mesh
halfedge_descriptor e_ij = prev(h,pmesh);
halfedge_descriptor e_ik = next(opposite(h,pmesh),pmesh);
halfedge_descriptor e_jk = prev(opposite(next(h,pmesh),pmesh),pmesh);
@@ -156,7 +156,7 @@ private:
{
int flips = 0;
std::list<halfedge_descriptor> interior_edges;
std::set<halfedge_descriptor> included_map;
std::set<halfedge_descriptor> included_map;
collect_interior_edges(faces, border_edges, interior_edges, included_map);
collect_interior_edges(new_faces, border_edges, interior_edges, included_map);
@@ -190,7 +190,7 @@ private:
do {
halfedge_descriptor h = *circ;
if (border_edges.find(h) == border_edges.end()){
// do not remove included_map and use if(&*h < &*oh) { interior_edges.push_back(h) }
// do not remove included_map and use if(&*h < &*oh) { interior_edges.push_back(h) }
// which will change the order of edges from run to run
halfedge_descriptor oh = opposite(h, pmesh);
halfedge_descriptor h_rep = (h < oh) ? h : oh; // AF: was &*h < &*oh
@@ -203,7 +203,7 @@ private:
}
double average_length(vertex_descriptor vh,
const std::set<face_descriptor>& interior_map,
const std::set<face_descriptor>& interior_map,
bool accept_internal_facets)
{
const Point_3& vp = get(vpmap, vh);
@@ -231,14 +231,14 @@ private:
void calculate_scale_attribute(const FaceRange& faces,
const std::set<face_descriptor>& interior_map,
std::map<vertex_descriptor, double>& scale_attribute,
bool accept_internal_facets)
bool accept_internal_facets)
{
for(face_descriptor fd : faces)
{
Halfedge_around_face_circulator<PolygonMesh> circ(halfedge(fd,pmesh),pmesh), done(circ);
do {
vertex_descriptor v = target(*circ,pmesh);
std::pair<typename std::map<vertex_descriptor, double>::iterator, bool> v_insert
std::pair<typename std::map<vertex_descriptor, double>::iterator, bool> v_insert
= scale_attribute.insert(std::make_pair(v, 0));
if(!v_insert.second) { continue; } // already calculated
v_insert.first->second = average_length(v, interior_map, accept_internal_facets);
@@ -262,7 +262,7 @@ private:
if(interior_map.find(f) == interior_map.end() || interior_map.find(f_op) == interior_map.end()) {
internal_v = false;
break;
}
}
} while(++circ_v != done_v);
if(internal_v) { return true; }
@@ -287,7 +287,7 @@ public:
// do not use just std::set, the order effects the output (for the same input we want to get same output)
std::set<face_descriptor> interior_map(boost::begin(faces), boost::end(faces));
// store boundary edges - to be used in relax
// store boundary edges - to be used in relax
std::set<halfedge_descriptor> border_edges;
for(face_descriptor f : faces)
{
@@ -1,709 +0,0 @@
// Copyright (c) 2019 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) : Sebastien Loriot,
// Mael Rouxel-Labbé
#ifndef CGAL_POLYGON_MESH_PROCESSING_REMOVE_DEGENERACIES_H
#define CGAL_POLYGON_MESH_PROCESSING_REMOVE_DEGENERACIES_H
#include <CGAL/license/Polygon_mesh_processing/repair.h>
#include <CGAL/Polygon_mesh_processing/shape_predicates.h>
#include <CGAL/Polygon_mesh_processing/measure.h>
#include <CGAL/Dynamic_property_map.h>
#include <CGAL/boost/graph/Euler_operations.h>
#include <CGAL/boost/graph/Named_function_parameters.h>
#include <CGAL/property_map.h>
#include <boost/graph/graph_traits.hpp>
#include <iostream>
#include <set>
#ifdef CGAL_PMP_DEBUG_REMOVE_DEGENERACIES
#include <sstream>
#include <fstream>
#endif
namespace CGAL {
namespace Polygon_mesh_processing {
namespace internal {
template <typename TriangleMesh, typename VPM, typename ECM, typename Traits>
std::array<typename boost::graph_traits<TriangleMesh>::halfedge_descriptor, 2>
is_badly_shaped(const typename boost::graph_traits<TriangleMesh>::face_descriptor f,
TriangleMesh& tmesh,
const VPM& vpm,
const ECM& ecm,
const Traits& gt,
const double cap_threshold, // angle over 160° ==> cap
const double needle_threshold, // longest edge / shortest edge over this ratio ==> needle
const double collapse_length_threshold) // max length of edges allowed to be collapsed
{
namespace PMP = CGAL::Polygon_mesh_processing;
typedef typename boost::graph_traits<TriangleMesh>::halfedge_descriptor halfedge_descriptor;
const halfedge_descriptor null_h = boost::graph_traits<TriangleMesh>::null_halfedge();
halfedge_descriptor res = PMP::is_needle_triangle_face(f, tmesh, needle_threshold,
parameters::vertex_point_map(vpm)
.geom_traits(gt));
if(res != null_h && !get(ecm, edge(res, tmesh)))
{
// don't want to collapse edges that are too large
if(collapse_length_threshold == 0 ||
edge_length(res, tmesh, parameters::vertex_point_map(vpm).geom_traits(gt)) <= collapse_length_threshold)
{
return make_array(res, null_h);
}
}
else // let's not make it possible to have a face be both a cap and a needle (for now)
{
res = PMP::is_cap_triangle_face(f, tmesh, cap_threshold, parameters::vertex_point_map(vpm).geom_traits(gt));
if(res != null_h && !get(ecm, edge(res, tmesh)))
return make_array(null_h, res);
}
return make_array(null_h, null_h);
}
template <typename TriangleMesh, typename EdgeContainer,
typename VPM, typename ECM, typename Traits>
void collect_badly_shaped_triangles(const typename boost::graph_traits<TriangleMesh>::face_descriptor f,
TriangleMesh& tmesh,
const VPM& vpm,
const ECM& ecm,
const Traits& gt,
const double cap_threshold, // angle over this threshold (as a cosine) ==> cap
const double needle_threshold, // longest edge / shortest edge over this ratio ==> needle
const double collapse_length_threshold, // max length of edges allowed to be collapsed
EdgeContainer& edges_to_collapse,
EdgeContainer& edges_to_flip)
{
typedef typename boost::graph_traits<TriangleMesh>::halfedge_descriptor halfedge_descriptor;
std::array<halfedge_descriptor, 2> res = is_badly_shaped(f, tmesh, vpm, ecm, gt, cap_threshold,
needle_threshold, collapse_length_threshold);
if(res[0] != boost::graph_traits<TriangleMesh>::null_halfedge())
{
#ifdef CGAL_PMP_DEBUG_REMOVE_DEGENERACIES
std::cout << "add new needle: " << edge(res[0], tmesh) << std::endl;
#endif
edges_to_collapse.insert(edge(res[0], tmesh));
}
else // let's not make it possible to have a face be both a cap and a needle (for now)
{
if(res[1] != boost::graph_traits<TriangleMesh>::null_halfedge())
{
#ifdef CGAL_PMP_DEBUG_REMOVE_DEGENERACIES
std::cout << "add new cap: " << edge(res[1],tmesh) << std::endl;
#endif
edges_to_flip.insert(edge(res[1], tmesh));
}
}
}
/*
// Following Ronfard et al. 96 we look at variation of the normal after the collapse
// the collapse must be topologically valid
template <class TriangleMesh, class NamedParameters>
bool is_collapse_geometrically_valid(typename boost::graph_traits<TriangleMesh>::halfedge_descriptor h,
const TriangleMesh& tmesh,
const NamedParameters& np)
{
using CGAL::parameters::choose_parameter;
using CGAL::parameters::get_parameter;
typedef typename boost::graph_traits<TriangleMesh>::halfedge_descriptor halfedge_descriptor;
typedef typename GetVertexPointMap<TriangleMesh, NamedParameters>::const_type VPM;
typedef typename boost::property_traits<VPM>::reference Point_ref;
typedef typename GetGeomTraits<TriangleMesh, NamedParameters>::type Traits;
VPM vpm = choose_parameter(get_parameter(np, internal_np::vertex_point),
get_const_property_map(vertex_point, tmesh));
Traits gt = choose_parameter(get_parameter(np, internal_np::geom_traits), Traits());
/// @todo handle boundary edges
h = opposite(h, tmesh); // Euler::collapse edge keeps the target and removes the source
// source is kept, target is removed
CGAL_assertion(target(h, tmesh) == vertex_removed);
Point_ref kept = get(vpm, source(h, tmesh));
Point_ref removed= get(vpm, target(h, tmesh));
// consider triangles incident to the vertex removed
halfedge_descriptor stop = prev(opposite(h, tmesh), tmesh);
halfedge_descriptor hi = opposite(next(h, tmesh), tmesh);
std::vector<halfedge_descriptor> triangles;
while(hi != stop)
{
if(!is_border(hi, tmesh))
{
Point_ref a = get(vpm, target(next(hi, tmesh), tmesh));
Point_ref b = get(vpm, source(hi, tmesh));
//ack a-b-point_remove and a-b-point_kept has a compatible orientation
/// @todo use a predicate
typename Traits::Vector_3 n1 = gt.construct_cross_product_vector_3_object()(removed-a, b-a);
typename Traits::Vector_3 n2 = gt.construct_cross_product_vector_3_object()(kept-a, b-a);
if(gt.compute_scalar_product_3_object()(n1, n2) <= 0)
return false;
}
hi = opposite(next(hi, tmesh), tmesh);
}
return true;
}
*/
template <class TriangleMesh, typename VPM, typename Traits>
boost::optional<double> get_collapse_volume(typename boost::graph_traits<TriangleMesh>::halfedge_descriptor h,
const TriangleMesh& tmesh,
const VPM& vpm,
const Traits& gt)
{
typedef typename boost::graph_traits<TriangleMesh>::halfedge_descriptor halfedge_descriptor;
typedef typename boost::property_traits<VPM>::reference Point_ref;
typedef typename Traits::Vector_3 Vector_3;
const typename Traits::Point_3 origin(ORIGIN);
/// @todo handle boundary edges
h = opposite(h, tmesh); // Euler::collapse edge keeps the target and removes the source
// source is kept, target is removed
Point_ref kept = get(vpm, source(h, tmesh));
Point_ref removed= get(vpm, target(h, tmesh));
// init volume with incident triangles (reversed orientation
double delta_vol = volume(removed, kept, get(vpm, target(next(h, tmesh), tmesh)), origin) +
volume(kept, removed, get(vpm, target(next(opposite(h, tmesh), tmesh), tmesh)), origin);
// consider triangles incident to the vertex removed
halfedge_descriptor stop = prev(opposite(h, tmesh), tmesh);
halfedge_descriptor hi = opposite(next(h, tmesh), tmesh);
std::vector<halfedge_descriptor> triangles;
while(hi != stop)
{
if(!is_border(hi, tmesh))
{
Point_ref a = get(vpm, target(next(hi, tmesh), tmesh));
Point_ref b = get(vpm, source(hi, tmesh));
//ack a-b-point_remove and a-b-point_kept has a compatible orientation
/// @todo use a predicate
Vector_3 v_ab = gt.construct_vector_3_object()(a, b);
Vector_3 v_ar = gt.construct_vector_3_object()(a, removed);
Vector_3 v_ak = gt.construct_vector_3_object()(a, kept);
Vector_3 n1 = gt.construct_cross_product_vector_3_object()(v_ar, v_ab);
Vector_3 n2 = gt.construct_cross_product_vector_3_object()(v_ak, v_ab);
if(gt.compute_scalar_product_3_object()(n1, n2) <= 0)
return boost::none;
delta_vol += volume(b, a, removed, origin) + volume(a, b, kept, origin); // opposite orientation
}
hi = opposite(next(hi, tmesh), tmesh);
}
return CGAL::abs(delta_vol);
}
template <typename TriangleMesh, typename VPM, typename VCM, typename Traits>
typename boost::graph_traits<TriangleMesh>::halfedge_descriptor
get_best_edge_orientation(typename boost::graph_traits<TriangleMesh>::edge_descriptor e,
const TriangleMesh& tmesh,
const VPM& vpm,
const VCM& vcm,
const Traits& gt)
{
typedef typename boost::graph_traits<TriangleMesh>::halfedge_descriptor halfedge_descriptor;
halfedge_descriptor h = halfedge(e, tmesh), ho = opposite(h, tmesh);
CGAL_assertion(!get(vcm, source(h, tmesh)) || !get(vcm, target(h, tmesh)));
boost::optional<double> dv1 = get_collapse_volume(h, tmesh, vpm, gt);
boost::optional<double> dv2 = get_collapse_volume(ho, tmesh, vpm, gt);
// the resulting point of the collapse of a halfedge is the target of the halfedge before collapse
if(get(vcm, source(h, tmesh)))
return dv2 != boost::none ? ho
: boost::graph_traits<TriangleMesh>::null_halfedge();
if(get(vcm, target(h, tmesh)))
return dv1 != boost::none ? h
: boost::graph_traits<TriangleMesh>::null_halfedge();
if(dv1 != boost::none)
{
if(dv2 != boost::none)
return (*dv1 < *dv2) ? h : ho;
return h;
}
if(dv2 != boost::none)
return ho;
return boost::graph_traits<TriangleMesh>::null_halfedge();
}
// adapted from triangulate_faces
template <typename TriangleMesh, typename VPM, typename Traits>
bool should_flip(typename boost::graph_traits<TriangleMesh>::edge_descriptor e,
const TriangleMesh& tmesh,
const VPM& vpm,
const Traits& gt)
{
typedef typename boost::graph_traits<TriangleMesh>::halfedge_descriptor halfedge_descriptor;
typedef typename boost::property_traits<VPM>::reference Point_ref;
typedef typename Traits::Vector_3 Vector_3;
CGAL_precondition(!is_border(e, tmesh));
halfedge_descriptor h = halfedge(e, tmesh);
Point_ref p0 = get(vpm, target(h, tmesh));
Point_ref p1 = get(vpm, target(next(h, tmesh), tmesh));
Point_ref p2 = get(vpm, source(h, tmesh));
Point_ref p3 = get(vpm, target(next(opposite(h, tmesh), tmesh), tmesh));
/* Chooses the diagonal that will split the quad in two triangles that maximize
* the scalar product of of the un-normalized normals of the two triangles.
* The lengths of the un-normalized normals (computed using cross-products of two vectors)
* are proportional to the area of the triangles.
* Maximize the scalar product of the two normals will avoid skinny triangles,
* and will also taken into account the cosine of the angle between the two normals.
* In particular, if the two triangles are oriented in different directions,
* the scalar product will be negative.
*/
// CGAL::cross_product(p2-p1, p3-p2) * CGAL::cross_product(p0-p3, p1-p0);
// CGAL::cross_product(p1-p0, p1-p2) * CGAL::cross_product(p3-p2, p3-p0);
const Vector_3 v01 = gt.construct_vector_3_object()(p0, p1);
const Vector_3 v12 = gt.construct_vector_3_object()(p1, p2);
const Vector_3 v23 = gt.construct_vector_3_object()(p2, p3);
const Vector_3 v30 = gt.construct_vector_3_object()(p3, p0);
const double p1p3 = gt.compute_scalar_product_3_object()(
gt.construct_cross_product_vector_3_object()(v12, v23),
gt.construct_cross_product_vector_3_object()(v30, v01));
const Vector_3 v21 = gt.construct_opposite_vector_3_object()(v12);
const Vector_3 v03 = gt.construct_opposite_vector_3_object()(v30);
const double p0p2 = gt.compute_scalar_product_3_object()(
gt.construct_cross_product_vector_3_object()(v01, v21),
gt.construct_cross_product_vector_3_object()(v23, v03));
return p0p2 <= p1p3;
}
} // namespace internal
namespace experimental {
// @todo check what to use as priority queue with removable elements, set might not be optimal
template <typename FaceRange, typename TriangleMesh, typename NamedParameters>
bool remove_almost_degenerate_faces(const FaceRange& face_range,
TriangleMesh& tmesh,
const double cap_threshold,
const double needle_threshold,
const double collapse_length_threshold,
const NamedParameters& np)
{
using CGAL::parameters::choose_parameter;
using CGAL::parameters::get_parameter;
typedef typename boost::graph_traits<TriangleMesh>::vertex_descriptor vertex_descriptor;
typedef typename boost::graph_traits<TriangleMesh>::halfedge_descriptor halfedge_descriptor;
typedef typename boost::graph_traits<TriangleMesh>::edge_descriptor edge_descriptor;
typedef typename boost::graph_traits<TriangleMesh>::face_descriptor face_descriptor;
typedef Constant_property_map<vertex_descriptor, bool> Default_VCM;
typedef typename internal_np::Lookup_named_param_def<internal_np::vertex_is_constrained_t,
NamedParameters,
Default_VCM>::type VCM;
VCM vcm_np = choose_parameter(get_parameter(np, internal_np::vertex_is_constrained), Default_VCM(false));
typedef Constant_property_map<edge_descriptor, bool> Default_ECM;
typedef typename internal_np::Lookup_named_param_def<internal_np::edge_is_constrained_t,
NamedParameters,
Default_ECM>::type ECM;
ECM ecm = choose_parameter(get_parameter(np, internal_np::edge_is_constrained), Default_ECM(false));
typedef typename GetVertexPointMap<TriangleMesh, NamedParameters>::const_type VPM;
VPM vpm = choose_parameter(get_parameter(np, internal_np::vertex_point),
get_const_property_map(vertex_point, tmesh));
typedef typename GetGeomTraits<TriangleMesh, NamedParameters>::type Traits;
Traits gt = choose_parameter(get_parameter(np, internal_np::geom_traits), Traits());
// Vertex property map that combines the VCM and the fact that extremities of a constrained edge should be constrained
typedef CGAL::dynamic_vertex_property_t<bool> Vertex_property_tag;
typedef typename boost::property_map<TriangleMesh, Vertex_property_tag>::type DVCM;
DVCM vcm = get(Vertex_property_tag(), tmesh);
for(face_descriptor f : face_range)
{
if(f == boost::graph_traits<TriangleMesh>::null_face())
continue;
for(halfedge_descriptor h : CGAL::halfedges_around_face(halfedge(f, tmesh), tmesh))
{
if(get(ecm, edge(h, tmesh)))
{
put(vcm, source(h, tmesh), true);
put(vcm, target(h, tmesh), true);
}
else if(get(vcm_np, target(h, tmesh)))
{
put(vcm, target(h, tmesh), true);
}
}
}
// Start the process of removing bad elements
std::set<edge_descriptor> edges_to_collapse;
std::set<edge_descriptor> edges_to_flip;
// @todo could probably do something a bit better by looping edges, consider the incident faces
// f1 / f2 and look at f1 if f1<f2, and the edge is smaller than the two other edges...
for(face_descriptor f : face_range)
internal::collect_badly_shaped_triangles(f, tmesh, vpm, ecm, gt,
cap_threshold, needle_threshold, collapse_length_threshold,
edges_to_collapse, edges_to_flip);
#ifdef CGAL_PMP_DEBUG_REMOVE_DEGENERACIES
int iter = 0;
#endif
for(;;)
{
bool something_was_done = false;
#ifdef CGAL_PMP_DEBUG_REMOVE_DEGENERACIES
std::cout << edges_to_collapse.size() << " needles and " << edges_to_flip.size() << " caps" << std::endl;
std::ostringstream oss;
oss << "degen_cleaning_iter_" << iter++ << ".off";
std::ofstream out(oss.str().c_str());
out << std::setprecision(17);
out << tmesh;
out.close();
#endif
if(edges_to_collapse.empty() && edges_to_flip.empty())
return true;
/// @todo maybe using a priority queue handling the more almost degenerate elements should be used
std::set<edge_descriptor> next_edges_to_collapse;
std::set<edge_descriptor> next_edges_to_flip;
// treat needles
#ifdef CGAL_PMP_DEBUG_REMOVE_DEGENERACIES_EXTRA
int kk=0;
std::ofstream(std::string("tmp/n-00000.off")) << tmesh;
#endif
while(!edges_to_collapse.empty())
{
edge_descriptor e = *edges_to_collapse.begin();
edges_to_collapse.erase(edges_to_collapse.begin());
CGAL_assertion(!get(ecm, e));
if(get(vcm, source(e, tmesh)) && get(vcm, target(e, tmesh)))
continue;
#ifdef CGAL_PMP_DEBUG_REMOVE_DEGENERACIES
std::cout << " treat needle: " << e << " (" << tmesh.point(source (e, tmesh))
<< " --- " << tmesh.point(target(e, tmesh)) << ")" << std::endl;
#endif
if(CGAL::Euler::does_satisfy_link_condition(e, tmesh))
{
// the following edges are removed by the collapse
halfedge_descriptor h = halfedge(e, tmesh);
CGAL_assertion(!is_border(h, tmesh)); // because extracted from a face
std::array<halfedge_descriptor, 2> nc =
internal::is_badly_shaped(face(h, tmesh), tmesh, vpm, ecm, gt,
cap_threshold, needle_threshold, collapse_length_threshold);
if(nc[0] != h)
{
#ifdef CGAL_PMP_DEBUG_REMOVE_DEGENERACIES
std::cerr << "Warning: Needle criteria no longer verified " << tmesh.point(source(e, tmesh)) << " "
<< tmesh.point(target(e, tmesh)) << std::endl;
#endif
// the opposite edge might also have been inserted in the set and might still be a needle
h = opposite(h, tmesh);
if(is_border(h, tmesh))
continue;
nc = internal::is_badly_shaped(face(h, tmesh), tmesh, vpm, ecm, gt,
cap_threshold, needle_threshold,
collapse_length_threshold);
if(nc[0] != h)
continue;
}
for(int i=0; i<2; ++i)
{
if(!is_border(h, tmesh))
{
edge_descriptor pe = edge(prev(h, tmesh), tmesh);
edges_to_flip.erase(pe);
next_edges_to_collapse.erase(pe);
edges_to_collapse.erase(pe);
}
h = opposite(h, tmesh);
}
// pick the orientation of edge to keep the vertex minimizing the volume variation
h = internal::get_best_edge_orientation(e, tmesh, vpm, vcm, gt);
if(h == boost::graph_traits<TriangleMesh>::null_halfedge())
{
#ifdef CGAL_PMP_DEBUG_REMOVE_DEGENERACIES
std::cerr << "Warning: geometrically invalid edge collapse! "
<< tmesh.point(source(e, tmesh)) << " "
<< tmesh.point(target(e, tmesh)) << std::endl;
#endif
next_edges_to_collapse.insert(e);
continue;
}
edges_to_flip.erase(e);
next_edges_to_collapse.erase(e); // for edges added in faces incident to a vertex kept after a collapse
#ifdef CGAL_PMP_DEBUG_REMOVE_DEGENERACIES_EXTRA
std::cerr << " " << kk << " -- Collapsing " << tmesh.point(source(h, tmesh)) << " "
<< tmesh.point(target(h, tmesh)) << std::endl;
#endif
// moving to the midpoint is not a good idea. On a circle for example you might endpoint with
// a bad geometry because you iteratively move one point
// auto mp = midpoint(tmesh.point(source(h, tmesh)), tmesh.point(target(h, tmesh)));
vertex_descriptor v = Euler::collapse_edge(edge(h, tmesh), tmesh);
//tmesh.point(v) = mp;
// examine all faces incident to the vertex kept
for(halfedge_descriptor hv : halfedges_around_target(v, tmesh))
{
if(!is_border(hv, tmesh))
{
internal::collect_badly_shaped_triangles(face(hv, tmesh), tmesh, vpm, ecm, gt,
cap_threshold, needle_threshold, collapse_length_threshold,
edges_to_collapse, edges_to_flip);
}
}
#ifdef CGAL_PMP_DEBUG_REMOVE_DEGENERACIES_EXTRA
std::string nb = std::to_string(++kk);
if(kk<10) nb = std::string("0")+nb;
if(kk<100) nb = std::string("0")+nb;
if(kk<1000) nb = std::string("0")+nb;
if(kk<10000) nb = std::string("0")+nb;
std::ofstream(std::string("tmp/n-")+nb+std::string(".off")) << tmesh;
#endif
something_was_done = true;
}
else
{
#ifdef CGAL_PMP_DEBUG_REMOVE_DEGENERACIES
std::cerr << "Warning: uncollapsable edge! " << tmesh.point(source(e, tmesh)) << " "
<< tmesh.point(target(e, tmesh)) << std::endl;
#endif
next_edges_to_collapse.insert(e);
}
}
// treat caps
#ifdef CGAL_PMP_DEBUG_REMOVE_DEGENERACIES_EXTRA
kk=0;
std::ofstream(std::string("tmp/c-000.off")) << tmesh;
#endif
while(!edges_to_flip.empty())
{
edge_descriptor e = *edges_to_flip.begin();
edges_to_flip.erase(edges_to_flip.begin());
CGAL_assertion(!get(ecm, e));
if(get(vcm, source(e, tmesh)) && get(vcm, target(e, tmesh)))
continue;
#ifdef CGAL_PMP_DEBUG_REMOVE_DEGENERACIES
std::cout << "treat cap: " << e << " (" << tmesh.point(source(e, tmesh))
<< " --- " << tmesh.point(target(e, tmesh)) << ")" << std::endl;
#endif
halfedge_descriptor h = halfedge(e, tmesh);
std::array<halfedge_descriptor,2> nc = internal::is_badly_shaped(face(h, tmesh), tmesh, vpm, ecm, gt,
cap_threshold, needle_threshold,
collapse_length_threshold);
// First check the triangle is still a cap
if(nc[1] != h)
{
#ifdef CGAL_PMP_DEBUG_REMOVE_DEGENERACIES
std::cerr << "Warning: Cap criteria no longer verified " << tmesh.point(source(e, tmesh)) << " --- "
<< tmesh.point(target(e, tmesh)) << std::endl;
#endif
// the opposite edge might also have been inserted in the set and might still be a cap
h = opposite(h, tmesh);
if(is_border(h, tmesh))
continue;
nc = internal::is_badly_shaped(face(h, tmesh), tmesh, vpm, ecm, gt,
cap_threshold, needle_threshold, collapse_length_threshold);
if(nc[1] != h)
continue;
}
// special case on the border
if(is_border(opposite(h, tmesh), tmesh))
{
// remove the triangle
edges_to_flip.erase(edge(prev(h, tmesh), tmesh));
edges_to_flip.erase(edge(next(h, tmesh), tmesh));
next_edges_to_collapse.erase(edge(prev(h, tmesh), tmesh));
next_edges_to_collapse.erase(edge(next(h, tmesh), tmesh));
Euler::remove_face(h, tmesh);
something_was_done = true;
continue;
}
// condition for the flip to be valid (the edge to be created does not already exist)
if(!halfedge(target(next(h, tmesh), tmesh),
target(next(opposite(h, tmesh), tmesh), tmesh), tmesh).second)
{
if(!internal::should_flip(e, tmesh, vpm, gt))
{
#ifdef CGAL_PMP_DEBUG_REMOVE_DEGENERACIES
std::cout << "Flipping prevented: not the best diagonal" << std::endl;
#endif
next_edges_to_flip.insert(e);
continue;
}
#ifdef CGAL_PMP_DEBUG_REMOVE_DEGENERACIES
std::cout << "Flipping" << std::endl;
#endif
#ifdef CGAL_PMP_DEBUG_REMOVE_DEGENERACIES_EXTRA
std::cerr << "step " << kk << "\n";
std::cerr << " Flipping " << tmesh.point(source(h, tmesh)) << " "
<< tmesh.point(target(h, tmesh)) << std::endl;
#endif
Euler::flip_edge(h, tmesh);
CGAL_assertion(edge(h, tmesh) == e);
// handle face updates
for(int i=0; i<2; ++i)
{
CGAL_assertion(!is_border(h, tmesh));
std::array<halfedge_descriptor, 2> nc =
internal::is_badly_shaped(face(h, tmesh), tmesh, vpm, ecm, gt,
cap_threshold, needle_threshold, collapse_length_threshold);
if(nc[1] != boost::graph_traits<TriangleMesh>::null_halfedge())
{
if(edge(nc[1], tmesh) != e)
next_edges_to_flip.insert(edge(nc[1], tmesh));
}
else
{
if(nc[0] != boost::graph_traits<TriangleMesh>::null_halfedge())
{
next_edges_to_collapse.insert(edge(nc[0], tmesh));
}
}
h = opposite(h, tmesh);
}
something_was_done = true;
}
#ifdef CGAL_PMP_DEBUG_REMOVE_DEGENERACIES
else
{
std::cerr << "Warning: unflippable edge! " << tmesh.point(source(h, tmesh)) << " --- "
<< tmesh.point(target(h, tmesh)) << std::endl;
next_edges_to_flip.insert(e);
}
#endif
#ifdef CGAL_PMP_DEBUG_REMOVE_DEGENERACIES_EXTRA
std::string nb = std::to_string(++kk);
if(kk<10) nb = std::string("0")+nb;
if(kk<100) nb = std::string("0")+nb;
if(kk<1000) nb = std::string("0")+nb;
if(kk<10000) nb = std::string("0")+nb;
std::ofstream(std::string("tmp/c-")+nb+std::string(".off")) << tmesh;
#endif
}
std::swap(edges_to_collapse, next_edges_to_collapse);
std::swap(edges_to_flip, next_edges_to_flip);
if(!something_was_done)
return false;
}
return false;
}
template <typename FaceRange, typename TriangleMesh>
bool remove_almost_degenerate_faces(const FaceRange& face_range,
TriangleMesh& tmesh,
const double cap_threshold,
const double needle_threshold,
const double collapse_length_threshold)
{
return remove_almost_degenerate_faces(face_range, tmesh,
cap_threshold, needle_threshold, collapse_length_threshold,
CGAL::parameters::all_default());
}
template <typename TriangleMesh, typename CGAL_PMP_NP_TEMPLATE_PARAMETERS>
bool remove_almost_degenerate_faces(TriangleMesh& tmesh,
const double cap_threshold,
const double needle_threshold,
const double collapse_length_threshold,
const CGAL_PMP_NP_CLASS& np)
{
return remove_almost_degenerate_faces(faces(tmesh), tmesh, cap_threshold, needle_threshold,
collapse_length_threshold, np);
}
template<class TriangleMesh>
bool remove_almost_degenerate_faces(TriangleMesh& tmesh,
const double cap_threshold,
const double needle_threshold,
const double collapse_length_threshold)
{
return remove_almost_degenerate_faces(faces(tmesh), tmesh,
cap_threshold, needle_threshold, collapse_length_threshold,
CGAL::parameters::all_default());
}
} // namespace experimental
} // namespace Polygon_mesh_processing
} // namespace CGAL
#endif // CGAL_POLYGON_MESH_PROCESSING_REMOVE_DEGENERACIES_H
@@ -414,7 +414,7 @@ compute_face_face_intersection(const FaceRange& face_range1,
CGAL::Bbox_3 b1 = CGAL::Polygon_mesh_processing::bbox(tm1, np1),
b2 = CGAL::Polygon_mesh_processing::bbox(tm2, np2);
if(!CGAL::do_overlap(b1, b2))
{
return out;
@@ -457,7 +457,7 @@ compute_face_face_intersection(const FaceRange& face_range1,
// compute intersections filtered out by boxes
typedef typename GetGeomTraits<TM, NamedParameters1>::type GeomTraits;
GeomTraits gt = choose_parameter(get_parameter(np1, internal_np::geom_traits), GeomTraits());
GeomTraits gt = choose_parameter<GeomTraits>(get_parameter(np1, internal_np::geom_traits));
internal::Intersect_faces<TM,
GeomTraits,
@@ -575,7 +575,7 @@ compute_face_polyline_intersection( const FaceRange& face_range,
// compute intersections filtered out by boxes
typedef typename GetGeomTraits<TM, NamedParameters>::type GeomTraits;
GeomTraits gt = choose_parameter(get_parameter(np, internal_np::geom_traits), GeomTraits());
GeomTraits gt = choose_parameter<GeomTraits>(get_parameter(np, internal_np::geom_traits));
internal::Intersect_face_polyline<TM,
GeomTraits,
@@ -643,7 +643,7 @@ compute_face_polylines_intersection(const FaceRange& face_range,
using parameters::get_parameter;
CGAL_precondition(CGAL::is_triangle_mesh(tm));
CGAL::Bbox_3 b1,b2;
b1 = CGAL::Polygon_mesh_processing::bbox(tm, np);
for(std::size_t i =0; i< polyline_range.size(); ++i)
@@ -651,10 +651,10 @@ compute_face_polylines_intersection(const FaceRange& face_range,
b2 += CGAL::bbox_3(polyline_range[i].begin(),
polyline_range[i].end());
}
if(!CGAL::do_overlap(b1,b2))
return out;
typedef TriangleMesh TM;
typedef typename boost::graph_traits<TM>::face_descriptor face_descriptor;
typedef typename GetVertexPointMap<TM, NamedParameters>::const_type VertexPointMap;
@@ -711,7 +711,7 @@ compute_face_polylines_intersection(const FaceRange& face_range,
// compute intersections filtered out by boxes
typedef typename GetGeomTraits<TM, NamedParameters>::type GeomTraits;
GeomTraits gt = choose_parameter(get_parameter(np, internal_np::geom_traits), GeomTraits());
GeomTraits gt = choose_parameter<GeomTraits>(get_parameter(np, internal_np::geom_traits));
internal::Intersect_face_polylines<TM,
GeomTraits,
@@ -1207,6 +1207,7 @@ bool do_intersect(const TriangleMesh& tm1,
{
using parameters::choose_parameter;
using parameters::get_parameter;
bool test_overlap = choose_parameter(get_parameter(np1, internal_np::overlap_test),false) ||
choose_parameter(get_parameter(np2, internal_np::overlap_test),false);
@@ -1232,7 +1233,7 @@ bool do_intersect(const TriangleMesh& tm1,
VertexPointMap2 vpm2 = choose_parameter(get_parameter(np2, internal_np::vertex_point),
get_const_property_map(boost::vertex_point, tm2));
typedef typename GetGeomTraits<TriangleMesh, NamedParameters1>::type GeomTraits;
GeomTraits gt = choose_parameter(get_parameter(np1, internal_np::geom_traits), GeomTraits());
GeomTraits gt = choose_parameter<GeomTraits>(get_parameter(np1, internal_np::geom_traits));
return internal::is_mesh2_in_mesh1(tm1, tm2, vpm1, vpm2, gt) ||
internal::is_mesh2_in_mesh1(tm2, tm1, vpm2, vpm1, gt);
@@ -1590,7 +1591,7 @@ OutputIterator intersecting_meshes(const TriangleMeshRange& range,
typedef typename boost::range_value<NamedParametersRange>::type NP_rng;
typedef typename boost::range_value<TriangleMeshRange>::type TriangleMesh;
typedef typename GetGeomTraits<TriangleMesh, NamedParameters, NP_rng>::type GT;
GT gt = choose_parameter(get_parameter(np, internal_np::geom_traits), GT());
GT gt = choose_parameter<GT>(get_parameter(np, internal_np::geom_traits));
//get all the pairs of meshes intersecting (no strict inclusion test)
std::ptrdiff_t cutoff = 2000;
@@ -611,7 +611,7 @@ construct_point(const std::pair<typename boost::graph_traits<TriangleMesh>::face
VertexPointMap vpm = parameters::choose_parameter(parameters::get_parameter(np, internal_np::vertex_point),
get_const_property_map(boost::vertex_point, tm));
Geom_traits gt = choose_parameter(get_parameter(np, internal_np::geom_traits), Geom_traits());
Geom_traits gt = choose_parameter<Geom_traits>(get_parameter(np, internal_np::geom_traits));
halfedge_descriptor hd = halfedge(loc.first, tm);
const Point_reference p0 = get(vpm, source(hd, tm));
@@ -1126,7 +1126,7 @@ locate_in_face(const typename internal::Location_traits<TriangleMesh, NamedParam
VertexPointMap vpm = parameters::choose_parameter(parameters::get_parameter(np, internal_np::vertex_point),
get_const_property_map(boost::vertex_point, tm));
Geom_traits gt = choose_parameter(get_parameter(np, internal_np::geom_traits), Geom_traits());
Geom_traits gt = choose_parameter<Geom_traits>(get_parameter(np, internal_np::geom_traits));
FT snap_tolerance = choose_parameter(get_parameter(np, internal_np::snapping_tolerance), 0);
@@ -0,0 +1,461 @@
// Copyright (c) 2015-2019 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) : Sebastien Loriot,
// Mael Rouxel-Labbé
//
#ifndef CGAL_POLYGON_MESH_PROCESSING_MANIFOLDNESS_H
#define CGAL_POLYGON_MESH_PROCESSING_MANIFOLDNESS_H
#include <CGAL/license/Polygon_mesh_processing/repair.h>
#include <CGAL/Polygon_mesh_processing/internal/named_function_params.h>
#include <CGAL/Polygon_mesh_processing/internal/named_params_helper.h>
#include <CGAL/algorithm.h>
#include <CGAL/assertions.h>
#include <CGAL/boost/graph/helpers.h>
#include <CGAL/Dynamic_property_map.h>
#include <CGAL/iterator.h>
#include <CGAL/property_map.h>
#include <iterator>
#include <map>
#include <utility>
#include <vector>
namespace CGAL {
namespace Polygon_mesh_processing {
/// \ingroup PMP_repairing_grp
/// checks whether a vertex of a polygon mesh is non-manifold.
///
/// @tparam PolygonMesh a model of `HalfedgeListGraph`
///
/// @param v a vertex of `pm`
/// @param pm a triangle mesh containing `v`
///
/// \warning This function has linear runtime with respect to the size of the mesh.
///
/// \sa `duplicate_non_manifold_vertices()`
///
/// \return `true` if the vertex is non-manifold, `false` otherwise.
template <typename PolygonMesh>
bool is_non_manifold_vertex(typename boost::graph_traits<PolygonMesh>::vertex_descriptor v,
const PolygonMesh& pm)
{
typedef typename boost::graph_traits<PolygonMesh>::halfedge_descriptor halfedge_descriptor;
typedef CGAL::dynamic_halfedge_property_t<bool> Halfedge_property_tag;
typedef typename boost::property_map<PolygonMesh, Halfedge_property_tag>::const_type Visited_halfedge_map;
// Dynamic pmaps do not have default initialization values (yet)
Visited_halfedge_map visited_halfedges = get(Halfedge_property_tag(), pm);
for(halfedge_descriptor h : halfedges(pm))
put(visited_halfedges, h, false);
std::size_t incident_null_faces_counter = 0;
for(halfedge_descriptor h : halfedges_around_target(v, pm))
{
put(visited_halfedges, h, true);
if(CGAL::is_border(h, pm))
++incident_null_faces_counter;
}
if(incident_null_faces_counter > 1)
{
// The vertex is the sole connection between two connected components --> non-manifold
return true;
}
for(halfedge_descriptor h : halfedges(pm))
{
if(v == target(h, pm))
{
// Haven't seen that halfedge yet ==> more than one umbrella incident to 'v' ==> non-manifold
if(!get(visited_halfedges, h))
return true;
}
}
return false;
}
namespace internal {
template <typename G>
struct Vertex_collector
{
typedef typename boost::graph_traits<G>::vertex_descriptor vertex_descriptor;
bool has_old_vertex(const vertex_descriptor v) const { return collections.count(v) != 0; }
void tag_old_vertex(const vertex_descriptor v)
{
CGAL_precondition(!has_old_vertex(v));
collections[v];
}
void collect_vertices(vertex_descriptor v1, vertex_descriptor v2)
{
std::vector<vertex_descriptor>& verts = collections[v1];
if(verts.empty())
verts.push_back(v1);
verts.push_back(v2);
}
template<typename OutputIterator>
void dump(OutputIterator out)
{
typedef std::pair<const vertex_descriptor, std::vector<vertex_descriptor> > Pair_type;
for(const Pair_type& p : collections)
*out++ = p.second;
}
void dump(Emptyset_iterator) { }
std::map<vertex_descriptor, std::vector<vertex_descriptor> > collections;
};
template <typename PolygonMesh, typename VPM, typename ConstraintMap>
typename boost::graph_traits<PolygonMesh>::vertex_descriptor
create_new_vertex_for_sector(typename boost::graph_traits<PolygonMesh>::halfedge_descriptor sector_begin_h,
typename boost::graph_traits<PolygonMesh>::halfedge_descriptor sector_last_h,
PolygonMesh& pm,
const VPM& vpm,
const ConstraintMap& cmap)
{
typedef typename boost::graph_traits<PolygonMesh>::vertex_descriptor vertex_descriptor;
typedef typename boost::graph_traits<PolygonMesh>::halfedge_descriptor halfedge_descriptor;
vertex_descriptor old_vd = target(sector_begin_h, pm);
vertex_descriptor new_vd = add_vertex(pm);
put(vpm, new_vd, get(vpm, old_vd));
put(cmap, new_vd, true);
set_halfedge(new_vd, sector_begin_h, pm);
halfedge_descriptor h = sector_begin_h;
do
{
set_target(h, new_vd, pm);
if(h == sector_last_h)
break;
else
h = prev(opposite(h, pm), pm);
}
while(h != sector_begin_h); // for safety
CGAL_assertion(h != sector_begin_h);
return new_vd;
}
template <typename PolygonMesh, typename NamedParameters>
std::size_t make_umbrella_manifold(typename boost::graph_traits<PolygonMesh>::halfedge_descriptor h,
PolygonMesh& pm,
internal::Vertex_collector<PolygonMesh>& dmap,
const NamedParameters& np)
{
typedef typename boost::graph_traits<PolygonMesh>::vertex_descriptor vertex_descriptor;
typedef typename boost::graph_traits<PolygonMesh>::halfedge_descriptor halfedge_descriptor;
using parameters::get_parameter;
using parameters::choose_parameter;
typedef typename GetVertexPointMap<PolygonMesh, NamedParameters>::type VertexPointMap;
VertexPointMap vpm = choose_parameter(get_parameter(np, internal_np::vertex_point),
get_property_map(vertex_point, pm));
typedef typename internal_np::Lookup_named_param_def<internal_np::vertex_is_constrained_t,
NamedParameters,
Constant_property_map<vertex_descriptor, bool> // default (no constraint pmap)
>::type VerticesMap;
VerticesMap cmap = choose_parameter(get_parameter(np, internal_np::vertex_is_constrained),
Constant_property_map<vertex_descriptor, bool>(false));
std::size_t nb_new_vertices = 0;
vertex_descriptor old_v = target(h, pm);
put(cmap, old_v, true); // store the duplicates
// count the number of borders
int border_counter = 0;
halfedge_descriptor ih = h, done = ih, border_h = h;
do
{
if(is_border(ih, pm))
{
border_h = ih;
++border_counter;
}
ih = prev(opposite(ih, pm), pm);
}
while(ih != done);
bool is_non_manifold_within_umbrella = (border_counter > 1);
if(!is_non_manifold_within_umbrella)
{
const bool first_time_meeting_v = !dmap.has_old_vertex(old_v);
if(first_time_meeting_v)
{
// The star is manifold, so if it is the first time we have met that vertex,
// there is nothing to do, we just keep the same vertex.
set_halfedge(old_v, h, pm); // to ensure halfedge(old_v, pm) stays valid
dmap.tag_old_vertex(old_v); // so that we know we have met old_v already, next time, we'll have to duplicate
}
else
{
// This is not the canonical star associated to 'v'.
// Create a new vertex, and move the whole star to that new vertex
halfedge_descriptor last_h = opposite(next(h, pm), pm);
vertex_descriptor new_v = create_new_vertex_for_sector(h, last_h, pm, vpm, cmap);
dmap.collect_vertices(old_v, new_v);
nb_new_vertices = 1;
}
}
// if there is more than one sector, look at each sector and split them away from the main one
else
{
// the first manifold sector, described by two halfedges
halfedge_descriptor sector_start_h = border_h;
CGAL_assertion(is_border(border_h, pm));
bool should_stop = false;
bool is_main_sector = true;
do
{
CGAL_assertion(is_border(sector_start_h, pm));
// collect the sector and split it away if it must be
halfedge_descriptor sector_last_h = sector_start_h;
do
{
halfedge_descriptor next_h = prev(opposite(sector_last_h, pm), pm);
if(is_border(next_h, pm))
break;
sector_last_h = next_h;
}
while(sector_last_h != sector_start_h);
CGAL_assertion(!is_border(sector_last_h, pm));
CGAL_assertion(sector_last_h != sector_start_h);
halfedge_descriptor next_start_h = prev(opposite(sector_last_h, pm), pm);
// there are multiple CCs incident to this particular vertex, and we should create a new vertex
// if it's not the first umbrella around 'old_v' or not the first sector, but only not if it's
// both the first umbrella and first sector.
bool must_create_new_vertex = (!is_main_sector || dmap.has_old_vertex(old_v));
// In any case, we must set up the next pointer correctly
set_next(sector_start_h, opposite(sector_last_h, pm), pm);
if(must_create_new_vertex)
{
vertex_descriptor new_v = create_new_vertex_for_sector(sector_start_h, sector_last_h, pm, vpm, cmap);
dmap.collect_vertices(old_v, new_v);
++nb_new_vertices;
}
else
{
// We are in the first sector and first star, ensure that halfedge(old_v, pm) stays valid
set_halfedge(old_v, sector_start_h, pm);
}
is_main_sector = false;
sector_start_h = next_start_h;
should_stop = (sector_start_h == border_h);
}
while(!should_stop);
}
return nb_new_vertices;
}
} // end namespace internal
/// \ingroup PMP_repairing_grp
/// collects the non-manifold vertices (if any) present in the mesh. A non-manifold vertex `v` is returned
/// via one incident halfedge `h` such that `target(h, pm) = v` for all the umbrellas that `v` apppears in
/// (an <i>umbrella</i> being the set of faces incident to all the halfedges reachable by walking around `v`
/// using `hnext = prev(opposite(h, pm), pm)`, starting from `h`).
///
/// @tparam PolygonMesh a model of `HalfedgeListGraph`
/// @tparam OutputIterator a model of `OutputIterator` holding objects of type
/// `boost::graph_traits<PolygonMesh>::%halfedge_descriptor`
///
/// @param pm a triangle mesh
/// @param out the output iterator that collects halfedges incident to `v`
///
/// \sa `is_non_manifold_vertex()`
/// \sa `duplicate_non_manifold_vertices()`
///
/// \return the output iterator.
template <typename PolygonMesh, typename OutputIterator>
OutputIterator non_manifold_vertices(const PolygonMesh& pm,
OutputIterator out)
{
// Non-manifoldness can appear either:
// - if 'pm' is pinched at a vertex. While traversing the incoming halfedges at this vertex,
// we will meet strictly more than one border halfedge.
// - if there are multiple umbrellas around a vertex. In that case, we will find a non-visited
// halfedge that has for target a vertex that is already visited.
typedef typename boost::graph_traits<PolygonMesh>::vertex_descriptor vertex_descriptor;
typedef typename boost::graph_traits<PolygonMesh>::halfedge_descriptor halfedge_descriptor;
typedef CGAL::dynamic_vertex_property_t<bool> Vertex_bool_tag;
typedef typename boost::property_map<PolygonMesh, Vertex_bool_tag>::const_type Known_manifold_vertex_map;
typedef CGAL::dynamic_vertex_property_t<halfedge_descriptor> Vertex_halfedge_tag;
typedef typename boost::property_map<PolygonMesh, Vertex_halfedge_tag>::const_type Visited_vertex_map;
typedef CGAL::dynamic_halfedge_property_t<bool> Halfedge_property_tag;
typedef typename boost::property_map<PolygonMesh, Halfedge_property_tag>::const_type Visited_halfedge_map;
Known_manifold_vertex_map known_nm_vertices = get(Vertex_bool_tag(), pm);
Visited_vertex_map visited_vertices = get(Vertex_halfedge_tag(), pm);
Visited_halfedge_map visited_halfedges = get(Halfedge_property_tag(), pm);
halfedge_descriptor null_h = boost::graph_traits<PolygonMesh>::null_halfedge();
// Dynamic pmaps do not have default initialization values (yet)
for(vertex_descriptor v : vertices(pm))
{
put(known_nm_vertices, v, false);
put(visited_vertices, v, null_h);
}
for(halfedge_descriptor h : halfedges(pm))
put(visited_halfedges, h, false);
for(halfedge_descriptor h : halfedges(pm))
{
// If 'h' is not visited yet, we walk around the target of 'h' and mark these
// halfedges as visited. Thus, if we are here and the target is already marked as visited,
// it means that the vertex is non manifold.
if(!get(visited_halfedges, h))
{
put(visited_halfedges, h, true);
bool is_non_manifold = false;
vertex_descriptor v = target(h, pm);
if(get(visited_vertices, v) != null_h) // already seen this vertex, but not from this star
{
is_non_manifold = true;
// if this is the second time we visit that vertex and the first star was manifold, we have
// never reported the first star, but we must now
if(!get(known_nm_vertices, v))
*out++ = get(visited_vertices, v); // that's a halfedge of the first star we've seen 'v' in
}
else
{
// first time we meet this vertex, just mark it so, and keep the halfedge we found the vertex with in memory
put(visited_vertices, v, h);
}
// While walking the star of this halfedge, if we meet a border halfedge more than once,
// it means the mesh is pinched and we are also in the case of a non-manifold situation
halfedge_descriptor ih = h, done = ih;
int border_counter = 0;
do
{
put(visited_halfedges, ih, true);
if(is_border(ih, pm))
++border_counter;
ih = prev(opposite(ih, pm), pm);
}
while(ih != done);
if(border_counter > 1)
is_non_manifold = true;
if(is_non_manifold)
{
*out++ = h;
put(known_nm_vertices, v, true);
}
}
}
return out;
}
/// \ingroup PMP_repairing_grp
/// duplicates all the non-manifold vertices of the input mesh.
///
/// @tparam PolygonMesh a model of `HalfedgeListGraph` and `MutableHalfedgeGraph`
/// @tparam NamedParameters a sequence of \ref pmp_namedparameters "Named Parameters"
///
/// @param pm the surface mesh to be repaired
/// @param np optional \ref pmp_namedparameters "Named Parameters" described below
///
/// \cgalNamedParamsBegin
/// \cgalParamBegin{vertex_point_map} the property map with the points associated to the vertices of `pmesh`.
/// The type of this map is model of `ReadWritePropertyMap`.
/// If this parameter is omitted, an internal property map for
/// `CGAL::vertex_point_t` should be available in `PolygonMesh`
/// \cgalParamEnd
/// \cgalParamBegin{vertex_is_constrained_map} a writable property map with `vertex_descriptor`
/// as key and `bool` as `value_type`. `put(pmap, v, true)` will be called for each duplicated
/// vertices, as well as the original non-manifold vertex in the input mesh.
/// \cgalParamEnd
/// \cgalParamBegin{output_iterator} a model of `OutputIterator` with value type
/// `std::vector<vertex_descriptor>`. The first vertex of each vector is a non-manifold vertex
/// of the input mesh, followed by the new vertices that were created to fix this precise
/// non-manifold configuration.
/// \cgalParamEnd
/// \cgalNamedParamsEnd
///
/// \return the number of vertices created.
template <typename PolygonMesh, typename NamedParameters>
std::size_t duplicate_non_manifold_vertices(PolygonMesh& pm,
const NamedParameters& np)
{
using parameters::get_parameter;
using parameters::choose_parameter;
typedef boost::graph_traits<PolygonMesh> GT;
typedef typename GT::halfedge_descriptor halfedge_descriptor;
typedef typename internal_np::Lookup_named_param_def<internal_np::output_iterator_t,
NamedParameters,
Emptyset_iterator>::type Output_iterator;
Output_iterator out = choose_parameter(get_parameter(np, internal_np::output_iterator),
Emptyset_iterator());
std::vector<halfedge_descriptor> non_manifold_cones;
non_manifold_vertices(pm, std::back_inserter(non_manifold_cones));
internal::Vertex_collector<PolygonMesh> dmap;
std::size_t nb_new_vertices = 0;
if(!non_manifold_cones.empty())
{
for(halfedge_descriptor h : non_manifold_cones)
nb_new_vertices += internal::make_umbrella_manifold(h, pm, dmap, np);
dmap.dump(out);
}
return nb_new_vertices;
}
template <class PolygonMesh>
std::size_t duplicate_non_manifold_vertices(PolygonMesh& pm)
{
return duplicate_non_manifold_vertices(pm, parameters::all_default());
}
} // namespace Polygon_mesh_processing
} // namespace CGAL
#endif // CGAL_POLYGON_MESH_PROCESSING_MANIFOLDNESS_H
@@ -327,7 +327,7 @@ face_area(typename boost::graph_traits<TriangleMesh>::face_descriptor f,
halfedge_descriptor nhd = next(hd, tmesh);
typedef typename GetGeomTraits<TriangleMesh, CGAL_PMP_NP_CLASS>::type GT;
GT traits = choose_parameter(get_parameter(np, internal_np::geom_traits), GT());
GT traits = choose_parameter<GT>(get_parameter(np, internal_np::geom_traits));
return approximate_sqrt(traits.compute_squared_area_3_object()(get(vpm, source(hd, tmesh)),
get(vpm, target(hd, tmesh)),
@@ -501,7 +501,7 @@ volume(const TriangleMesh& tmesh,
typename GetVertexPointMap<TriangleMesh, CGAL_PMP_NP_CLASS>::const_type
vpm = choose_parameter(get_parameter(np, internal_np::vertex_point),
get_const_property_map(CGAL::vertex_point, tmesh));
get_const_property_map(CGAL::vertex_point, tmesh));
typename GetGeomTraits<TriangleMesh, CGAL_PMP_NP_CLASS>::type::Point_3 origin(0, 0, 0);
typedef typename boost::graph_traits<TriangleMesh>::face_descriptor face_descriptor;
@@ -578,11 +578,11 @@ face_aspect_ratio(typename boost::graph_traits<TriangleMesh>::face_descriptor f,
typename GetVertexPointMap<TriangleMesh, CGAL_PMP_NP_CLASS>::const_type
vpm = choose_parameter(get_parameter(np, internal_np::vertex_point),
get_const_property_map(CGAL::vertex_point, tmesh));
get_const_property_map(CGAL::vertex_point, tmesh));
halfedge_descriptor h = halfedge(f, tmesh);
Geom_traits gt = choose_parameter(get_parameter(np, internal_np::geom_traits), Geom_traits());
Geom_traits gt = choose_parameter<Geom_traits>(get_parameter(np, internal_np::geom_traits));
#if 0
const FT sq_triangle_area = gt.compute_squared_area_3_object()(get(vpm, source(h, tmesh)),
@@ -680,7 +680,7 @@ centroid(const TriangleMesh& tmesh, const CGAL_PMP_NP_CLASS& np)
typedef typename GetVertexPointMap<TriangleMesh, CGAL_PMP_NP_CLASS>::const_type Vpm;
Vpm vpm = choose_parameter(get_parameter(np, internal_np::vertex_point),
get_const_property_map(CGAL::vertex_point, tmesh));
get_const_property_map(CGAL::vertex_point, tmesh));
typedef typename GetGeomTraits<TriangleMesh, CGAL_PMP_NP_CLASS>::type Kernel;
typedef typename Kernel::Point_3 Point_3;
@@ -22,8 +22,11 @@
#include <CGAL/Polygon_mesh_processing/internal/named_params_helper.h>
#include <CGAL/Polygon_mesh_processing/stitch_borders.h>
#include <boost/unordered_set.hpp>
#include <boost/bind.hpp>
#include <algorithm>
#include <iostream>
#include <iterator>
#include <vector>
#include <utility>
namespace CGAL {
@@ -58,6 +61,87 @@ struct Less_on_point_of_target
const VertexPointMap& vpm;
};
// Given a container of vectors of halfedges whose target are geometrically indentical,
// check that the intervals described by these pairs are either disjoint or nested.
// This is done to ensure valid combinatorics when we merge the vertices.
// If incompatible (overlapping) intervals are found, the pair representating the longest
// interval (arbitrary choice) is removed from the candidate list.
template <typename VPM, typename PolygonMesh>
void sanitize_candidates(const std::vector<std::pair<typename boost::graph_traits<PolygonMesh>::halfedge_descriptor, std::size_t> >& cycle_hedges,
std::vector<std::vector<std::size_t> >& candidate_hedges_with_id,
const VPM vpm,
const PolygonMesh& pm)
{
if(candidate_hedges_with_id.empty())
return;
std::size_t nm_vertices_n = candidate_hedges_with_id.size();
for(std::size_t fr_id=0, fr_end=nm_vertices_n-1; fr_id<fr_end; ++fr_id)
{
std::vector<std::size_t>& first_candidates = candidate_hedges_with_id[fr_id];
CGAL_assertion(first_candidates.size() >= 2);
for(std::size_t i=0, ie=first_candidates.size()-1; i<ie; ++i)
{
const std::size_t first_left = cycle_hedges[first_candidates[i]].second;
const std::size_t first_right = cycle_hedges[first_candidates[i+1]].second;
CGAL_assertion(first_left < first_right);
for(std::size_t sr_id=i+1, sr_end=nm_vertices_n; sr_id<sr_end; ++sr_id)
{
std::vector<std::size_t>& second_candidates = candidate_hedges_with_id[sr_id];
CGAL_assertion(second_candidates.size() >= 2);
for(std::size_t j=0, je=second_candidates.size()-1; j<je; ++j)
{
const std::size_t second_left = cycle_hedges[second_candidates[j]].second;
const std::size_t second_right = cycle_hedges[second_candidates[j+1]].second;
CGAL_assertion(second_left < second_right);
// The pair of intervals should be either disjoint or nested
// so reject:
// sl -- fl -- sr -- fr and fl -- sl -- fr -- sr
if((second_left < first_left && first_left < second_right && second_right < first_right) ||
(first_left < second_left && second_left < first_right && first_right < second_right))
{
// Remove the candidate with largest range
const std::size_t first_candidates_range =
cycle_hedges[first_candidates.back()].second - cycle_hedges[first_candidates.front()].second;
const std::size_t second_candidates_range =
cycle_hedges[second_candidates.back()].second - cycle_hedges[second_candidates.front()].second;
CGAL_assertion(first_candidates_range <= cycle_hedges.size());
CGAL_assertion(second_candidates_range <= cycle_hedges.size());
#ifdef CGAL_PMP_MERGE_BORDER_VERTICES_DEBUG
std::cout << "Incompatible ranges:\n";
std::cout << "first range: " << first_left << " to " << first_right << std::endl;
std::cout << "second range: " << second_left << " to " << second_right << std::endl;
std::cout << "Full ranges:" << std::endl;
std::cout << cycle_hedges[first_candidates.front()].second << " to " << cycle_hedges[first_candidates.back()].second;
std::cout << " (" << first_candidates.size() << " halfedges)";
std::cout << " at " << get(vpm, target(cycle_hedges[first_candidates.front()].first, pm)) << std::endl;
std::cout << cycle_hedges[second_candidates.front()].second << " to " << cycle_hedges[second_candidates.back()].second;
std::cout << " (" << second_candidates.size() << " halfedges)";
std::cout << " at " << get(vpm, target(cycle_hedges[second_candidates.front()].first, pm)) << std::endl;
#endif
std::vector<std::vector<std::size_t> >::iterator to_remove_iter = candidate_hedges_with_id.begin();
if(first_candidates_range > second_candidates_range)
std::advance(to_remove_iter, fr_id);
else
std::advance(to_remove_iter, sr_id);
candidate_hedges_with_id.erase(to_remove_iter);
// restart the whole thing
return sanitize_candidates(cycle_hedges, candidate_hedges_with_id, vpm, pm);
}
} // entries of the second range
} // second range
} // entries of the first range
} // first range
}
// warning: cycle_hedges will be altered (sorted)
template <class PolygonMesh, class Vpm, class halfedge_descriptor>
@@ -75,26 +159,23 @@ void detect_identical_mergeable_vertices(
std::size_t nbv=cycle_hedges.size();
std::size_t i=1;
std::set< std::pair<std::size_t, std::size_t> > intervals;
// IDs of cycle_hedges
std::vector<std::vector<std::size_t> > candidate_hedges_with_id;
while(i!=nbv)
while(i != nbv)
{
if ( get(vpm, target(cycle_hedges[i].first, pm)) ==
get(vpm, target(cycle_hedges[i-1].first, pm)) )
if(get(vpm, target(cycle_hedges[i].first, pm)) ==
get(vpm, target(cycle_hedges[i-1].first, pm)) )
{
hedges_with_identical_point_target.push_back( std::vector<halfedge_descriptor>() );
hedges_with_identical_point_target.back().push_back(cycle_hedges[i-1].first);
hedges_with_identical_point_target.back().push_back(cycle_hedges[i].first);
intervals.insert( std::make_pair(cycle_hedges[i-1].second, cycle_hedges[i].second) );
std::size_t previous = cycle_hedges[i].second;
while(++i!=nbv)
candidate_hedges_with_id.resize(candidate_hedges_with_id.size() + 1);
candidate_hedges_with_id.back().push_back(i-1);
candidate_hedges_with_id.back().push_back(i);
while(++i != nbv)
{
if ( get(vpm, target(cycle_hedges[i].first, pm)) ==
get(vpm, target(cycle_hedges[i-1].first, pm)) )
if(get(vpm, target(cycle_hedges[i].first, pm)) ==
get(vpm, target(cycle_hedges[i-1].first, pm)))
{
hedges_with_identical_point_target.back().push_back(cycle_hedges[i].first);
intervals.insert( std::make_pair(previous, cycle_hedges[i].second) );
previous = cycle_hedges[i].second;
candidate_hedges_with_id.back().push_back(i);
}
else
{
@@ -104,29 +185,22 @@ void detect_identical_mergeable_vertices(
}
}
else
{
++i;
}
}
// check that intervals are disjoint or strictly nested
// if there is only one issue we drop the whole cycle.
// @todo shall we try to be more conservative?
if (hedges_with_identical_point_target.empty()) return;
std::set< std::pair<std::size_t, std::size_t> >::iterator it1 = intervals.begin(),
end2 = intervals.end(),
end1 = std::prev(end2),
it2;
for (; it1!=end1; ++it1)
for(it2=std::next(it1); it2!= end2; ++it2 )
// Check that intervals are disjoint or strictly nested
sanitize_candidates(cycle_hedges, candidate_hedges_with_id, vpm, pm);
for(const std::vector<std::size_t>& candidates : candidate_hedges_with_id)
{
hedges_with_identical_point_target.resize(hedges_with_identical_point_target.size() + 1);
for(const std::size_t hid : candidates)
{
CGAL_assertion(it1->first<it2->first);
CGAL_assertion(it1->first < it1->second && it2->first < it2->second);
if (it1->second > it2->first && it2->second > it1->second)
{
std::cerr << "Merging is skipt to avoid bad cycle connections\n";
hedges_with_identical_point_target.clear();
return;
}
hedges_with_identical_point_target.back().push_back(cycle_hedges[hid].first);
}
}
}
// \ingroup PMP_repairing_grp
@@ -447,7 +447,7 @@ struct Polygon_soup_orienter
* amounts to duplicate the polygon to which it belongs.
*
* These points are either an endpoint of an edge incident to more
* than two polygons, an endpoint of an edge between
* than two polygons, an endpoint of an edge between
* two polygons with incompatible orientations (during the re-orientation process),
* or more generally a point \a p at which the intersection
* of an infinitesimally small ball centered at \a p
@@ -0,0 +1,196 @@
// Copyright (c) 2019 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) : Sebastien Loriot and Maxime Gimeno
#ifndef CGAL_ORIENT_POLYGON_SOUP_EXTENSION_H
#define CGAL_ORIENT_POLYGON_SOUP_EXTENSION_H
#include <CGAL/license/Polygon_mesh_processing/repair.h>
#include <CGAL/AABB_tree.h>
#include <CGAL/AABB_face_graph_triangle_primitive.h>
#include <CGAL/AABB_traits.h>
#include <CGAL/Polygon_mesh_processing/shape_predicates.h>
#include <CGAL/Polygon_mesh_processing/compute_normal.h>
#include <CGAL/Polygon_mesh_processing/orient_polygon_soup.h>
#include <boost/iterator/counting_iterator.hpp>
#include <boost/iterator/filter_iterator.hpp>
#ifdef CGAL_LINKED_WITH_TBB
#include <tbb/parallel_for.h>
#endif // CGAL_LINKED_WITH_TBB
namespace CGAL {
namespace Polygon_mesh_processing {
/*!
* \ingroup PMP_orientation_grp
* duplicates each point \a p at which the intersection
* of an infinitesimally small ball centered at \a p
* with the polygons incident to it is not a topological disk.
*
* @tparam PointRange a model of the concepts `RandomAccessContainer`
* and `BackInsertionSequence` whose `value_type` is the point type.
* @tparam PolygonRange a model of the concept `RandomAccessContainer`
* whose `value_type` is a model of the concept `RandomAccessContainer`
* whose `value_type` is `std::size_t`, and is also a model of `BackInsertionSequence`.
*
* @param points points of the soup of polygons. Some additional points might be pushed back to resolve
* non-manifoldness or non-orientability issues.
* @param polygons each element in the vector describes a polygon using the indices of the points in `points`.
* If needed the order of the indices of a polygon might be reversed.
* @return `false` if some points were duplicated, thus producing a self-intersecting surface mesh.
* @return `true` otherwise.
* @sa `orient_polygon_soup()`
*/
template <class PointRange, class PolygonRange>
bool
duplicate_non_manifold_edges_in_polygon_soup(PointRange& points,
PolygonRange& polygons)
{
std::size_t inital_nb_pts = points.size();
typedef CGAL::Polygon_mesh_processing::internal::
Polygon_soup_orienter<PointRange, PolygonRange> Orienter;
Orienter orienter(points, polygons);
orienter.fill_edge_map();
// make edges to duplicate
for(std::size_t i1=0;i1<points.size();++i1)
for(const typename Orienter::Internal_map_type::value_type& i2_and_pids : orienter.edges[i1])
if (i2_and_pids.second.size() > 1)
orienter.set_edge_marked(i1,i2_and_pids.first,orienter.marked_edges);
orienter.duplicate_singular_vertices();
return inital_nb_pts==points.size();
}
/*!
* \ingroup PMP_orientation_grp
* orients each triangle of a triangle soup using the orientation of its
* closest non degenerate triangle in `tm_ref`.
* \tparam Concurrency_tag enables sequential versus parallel orientation.
Possible values are `Sequential_tag` (the default),
`Parallel_if_available_tag`, and `Parallel_tag`.
* \tparam PointRange a model of the concepts `RandomAccessContainer`
* and `BackInsertionSequence` whose value type is the point type.
* @tparam TriangleRange a model of the concept `RandomAccessContainer`
* whose `value_type` is a model of the concept `RandomAccessContainer`
* whose `value_type` is `std::size_t`and of size 3.
* @tparam TriangleMesh a model of `FaceListGraph` and `MutableFaceGraph` .
*
* \param tm_ref the reference triangle_mesh.
* \param points the points of the soup.
* \param triangles the triangles of the soup.
* @param np optional sequence of \ref pmp_namedparameters among the ones listed below
*
* \cgalNamedParamsBegin
* \cgalParamBegin{vertex_point_map}
* the property map with the points associated to the vertices of `tm_ref`.
* If this parameter is omitted, an internal property map for
* `CGAL::vertex_point_t` must be available in `TriangleMesh`
* \cgalParamEnd
* \cgalParamBegin{geom_traits} a geometric traits class instance.
* The traits class must provide the nested functor `Collinear_3`
* to check whether three points are collinear.
* \cgalParamEnd
* \cgalNamedParamsEnd
*
* \attention The types of points in `PointRange`, `geom_traits` and `vertex_point_map` must be the same.
*/
template <class Concurrency_tag = Sequential_tag, class PointRange, class TriangleRange,
class TriangleMesh, class NamedParameters>
void
orient_triangle_soup_with_reference_triangle_mesh(
const TriangleMesh& tm_ref,
PointRange& points,
TriangleRange& triangles,
const NamedParameters& np)
{
namespace PMP = CGAL::Polygon_mesh_processing;
typedef boost::graph_traits<TriangleMesh> GrT;
typedef typename GrT::face_descriptor face_descriptor;
typedef typename PointRange::value_type Point_3;
typedef typename GetGeomTraits<TriangleMesh, NamedParameters>::type K;
typedef typename
GetVertexPointMap<TriangleMesh, NamedParameters>::const_type Vpm;
Vpm vpm = parameters::choose_parameter(parameters::get_parameter(np, internal_np::vertex_point),
get_const_property_map(CGAL::vertex_point, tm_ref));
typedef std::function<bool(face_descriptor)> Face_predicate;
Face_predicate is_not_deg =
[&tm_ref, np](face_descriptor f)
{
return !PMP::is_degenerate_triangle_face(f, tm_ref, np);
};
// build a tree filtering degenerate faces
typedef CGAL::AABB_face_graph_triangle_primitive<TriangleMesh, Vpm> Primitive;
typedef CGAL::AABB_traits<K, Primitive> Tree_traits;
boost::filter_iterator<Face_predicate, typename GrT::face_iterator>
begin(is_not_deg, faces(tm_ref).begin(), faces(tm_ref).end()),
end(is_not_deg, faces(tm_ref).end(), faces(tm_ref).end());
CGAL::AABB_tree<Tree_traits> tree(begin, end, tm_ref, vpm);
// now orient the faces
tree.build();
tree.accelerate_distance_queries();
auto process_facet =
[&points, &tree, &tm_ref, &triangles](std::size_t fid) {
const Point_3& p0 = points[triangles[fid][0]];
const Point_3& p1 = points[triangles[fid][1]];
const Point_3& p2 = points[triangles[fid][2]];
const Point_3 mid = CGAL::centroid(p0, p1, p2);
std::pair<Point_3, face_descriptor> pt_and_f =
tree.closest_point_and_primitive(mid);
auto face_ref_normal = PMP::compute_face_normal(pt_and_f.second, tm_ref);
if(face_ref_normal * cross_product(p1-p0, p2-p0) < 0) {
std::swap(triangles[fid][1], triangles[fid][2]);
}
};
#if !defined(CGAL_LINKED_WITH_TBB)
CGAL_static_assertion_msg (!(boost::is_convertible<Concurrency_tag, CGAL::Parallel_tag>::value),
"Parallel_tag is enabled but TBB is unavailable.");
#else
if (boost::is_convertible<Concurrency_tag,CGAL::Parallel_tag>::value)
tbb::parallel_for(std::size_t(0), triangles.size(), std::size_t(1), process_facet);
else
#endif
std::for_each(
boost::counting_iterator<std::size_t> (0),
boost::counting_iterator<std::size_t> (triangles.size()),
process_facet);
}
template <class Concurrency_tag = Sequential_tag, class PointRange, class TriangleRange,
class TriangleMesh>
void
orient_triangle_soup_with_reference_triangle_mesh(
const TriangleMesh& tm_ref,
PointRange& points,
TriangleRange& triangles)
{
orient_triangle_soup_with_reference_triangle_mesh<Concurrency_tag>(tm_ref, points, triangles, CGAL::parameters::all_default());
}
}}//end namespace CGAL::Polygon_mesh_processing
#endif // CGAL_ORIENT_POLYGON_SOUP_EXTENSION_H
File diff suppressed because it is too large Load Diff
@@ -0,0 +1,129 @@
// Copyright (c) 2020 GeometryFactory (France).
// All rights reserved.
//
// This file is part of CGAL (www.cgal.org).
//
// $URL$
// $Id$
// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
//
//
// Author(s) : Mael Rouxel-Labbé
#ifndef CGAL_POLYGON_MESH_PROCESSING_POLYGON_MESH_TO_POLYGON_SOUP_H
#define CGAL_POLYGON_MESH_PROCESSING_POLYGON_MESH_TO_POLYGON_SOUP_H
#include <CGAL/license/Polygon_mesh_processing/repair.h>
#include <CGAL/algorithm.h>
#include <CGAL/boost/graph/iterator.h>
#include <CGAL/boost/graph/named_params_helper.h>
#include <CGAL/Container_helper.h>
#include <CGAL/Dynamic_property_map.h>
#include <CGAL/property_map.h>
#include <boost/range/value_type.hpp>
#include <boost/range/reference.hpp>
namespace CGAL {
namespace Polygon_mesh_processing {
/// \ingroup PMP_repairing_grp
///
/// Adds the vertices and faces of a mesh into a (possibly non-empty) polygon soup.
///
/// \tparam PolygonMesh a model of `FaceListGraph`
/// \tparam PointRange a model of the concepts `RandomAccessContainer` and
/// `BackInsertionSequence` whose value type can be constructed from
/// the point type of the polygon mesh
/// \tparam PolygonRange a model of the concepts `RandomAccessContainer` and `BackInsertionSequence` whose
/// value type is itself a model of the concepts `RandomAccessContainer` and
/// `BackInsertionSequence` whose value type is `std::size_t`
/// \tparam NamedParameters a sequence of \ref pmp_namedparameters "Named Parameters"
///
/// \param mesh the mesh whose faces are being put in the polygon soup
/// \param points points making the polygons of the soup
/// \param polygons each element in the vector describes a polygon using the indices of the points in `points`
/// \param np optional sequence of \ref pmp_namedparameters "Named Parameters" among the ones listed below
///
/// \cgalNamedParamsBegin
/// \cgalParamBegin{vertex_point_map}
/// a model of `ReadablePropertyMap`, the property map with the points associated to the vertices of `mesh`.
/// If this parameter is omitted, an internal property map for
/// `CGAL::vertex_point_t` must be available in `PolygonMesh`.
/// \cgalParamEnd
/// \cgalNamedParamsEnd
///
/// \cgalAdvancedBegin
/// `PolygonRange` can also be a model of the concepts `RandomAccessContainer` and `BackInsertionSequence`
/// whose value type is an array, but it is the user's responsability to ensure that
/// all faces have the same number of vertices, and that this number is equal to the size of the array.
/// \cgalAdvancedEnd
///
/// \sa `CGAL::Polygon_mesh_processing::orient_polygon_soup()`
/// \sa `CGAL::Polygon_mesh_processing::is_polygon_soup_a_polygon_mesh()`
/// \sa `CGAL::Polygon_mesh_processing::polygon_soup_to_polygon_mesh()`
///
template<typename PolygonMesh,
typename PointRange, typename PolygonRange,
typename NamedParameters>
void polygon_mesh_to_polygon_soup(const PolygonMesh& mesh,
PointRange& points,
PolygonRange& polygons,
const NamedParameters& np)
{
typedef typename boost::graph_traits<PolygonMesh>::vertex_descriptor vertex_descriptor;
typedef typename boost::graph_traits<PolygonMesh>::halfedge_descriptor halfedge_descriptor;
typedef typename boost::graph_traits<PolygonMesh>::face_descriptor face_descriptor;
using parameters::choose_parameter;
using parameters::get_parameter;
typedef typename GetVertexPointMap<PolygonMesh, NamedParameters>::const_type VPM;
VPM vpm = choose_parameter(get_parameter(np, internal_np::vertex_point),
get_const_property_map(vertex_point, mesh));
typedef CGAL::dynamic_vertex_property_t<std::size_t> Vertex_index;
typedef typename boost::property_map<PolygonMesh, Vertex_index>::const_type VIM;
VIM vim = get(Vertex_index(), mesh);
typedef typename boost::range_value<PolygonRange>::type Polygon;
std::size_t index = points.size(); // so that multiple meshes can be put into the same soup
CGAL::internal::reserve(points, points.size() + vertices(mesh).size());
CGAL::internal::reserve(polygons, polygons.size() + faces(mesh).size());
for(const vertex_descriptor v : vertices(mesh))
{
points.emplace_back(get(vpm, v));
put(vim, v, index++);
}
for(const face_descriptor f : faces(mesh))
{
CGAL::Iterator_range<CGAL::Halfedge_around_face_iterator<PolygonMesh> > incident_halfedges =
CGAL::halfedges_around_face(halfedge(f, mesh), mesh);
Polygon polygon;
CGAL::internal::resize(polygon, incident_halfedges.size());
std::size_t pos = 0;
for(halfedge_descriptor h : incident_halfedges)
polygon[pos++] = get(vim, target(h, mesh));
polygons.push_back(polygon);
}
}
template<typename PolygonMesh, typename PointRange, typename PolygonRange>
void polygon_mesh_to_polygon_soup(const PolygonMesh& mesh,
PointRange& points,
PolygonRange& polygons)
{
return polygon_mesh_to_polygon_soup(mesh, points, polygons, CGAL::parameters::all_default());
}
} // namespace Polygon_mesh_processing
} // namespace CGAL
#endif // CGAL_POLYGON_MESH_PROCESSING_POLYGON_MESH_TO_POLYGON_SOUP_H
@@ -10,46 +10,40 @@
//
// Author(s) : Laurent Rineau and Ilker O. Yaz
#ifndef CGAL_POLYGON_MESH_PROCESSING_POLYGON_SOUP_TO_POLYGON_MESH
#define CGAL_POLYGON_MESH_PROCESSING_POLYGON_SOUP_TO_POLYGON_MESH
#ifndef CGAL_POLYGON_MESH_PROCESSING_POLYGON_SOUP_TO_POLYGON_MESH_H
#define CGAL_POLYGON_MESH_PROCESSING_POLYGON_SOUP_TO_POLYGON_MESH_H
#include <CGAL/license/Polygon_mesh_processing/repair.h>
#include <CGAL/disable_warnings.h>
#include <CGAL/boost/graph/Euler_operations.h>
#include <CGAL/property_map.h>
#include <CGAL/Polygon_mesh_processing/orient_polygon_soup.h>
#include <CGAL/algorithm.h>
#include <set>
#include <boost/dynamic_bitset.hpp>
#include <CGAL/algorithm.h>
#include <CGAL/boost/graph/Euler_operations.h>
#include <CGAL/boost/graph/iterator.h>
#include <CGAL/boost/graph/named_params_helper.h>
#include <CGAL/Dynamic_property_map.h>
#include <CGAL/property_map.h>
#include <boost/dynamic_bitset.hpp>
#include <boost/range/size.hpp>
#include <boost/range/value_type.hpp>
#include <boost/range/reference.hpp>
namespace CGAL
{
namespace Polygon_mesh_processing
{
namespace internal
{
template <typename PM
, typename PointRange
, typename PolygonRange>
class Polygon_soup_to_polygon_mesh
{
const PointRange& _points;
const PolygonRange& _polygons;
#include <set>
#include <type_traits>
#include <vector>
typedef typename boost::property_map<PM, CGAL::vertex_point_t>::type Vpmap;
typedef typename boost::property_traits<Vpmap>::value_type Point_3;
namespace CGAL {
namespace Polygon_mesh_processing {
namespace internal {
typedef typename boost::graph_traits<PM>::vertex_descriptor vertex_descriptor;
typedef typename boost::graph_traits<PM>::halfedge_descriptor halfedge_descriptor;
typedef typename boost::graph_traits<PM>::face_descriptor face_descriptor;
typedef typename PolygonRange::value_type Polygon;
typedef typename PointRange::value_type Point;
template <typename PointRange,
typename PolygonRange,
typename PointMap = typename CGAL::GetPointMap<PointRange>::const_type>
class PS_to_PM_converter
{
typedef typename boost::range_value<PolygonRange>::type Polygon;
typedef typename boost::property_traits<PointMap>::value_type Point;
public:
/**
@@ -57,178 +51,237 @@ public:
* @param points points of the soup of polygons.
* @param polygons each element in the range describes a polygon using the index of the points in the range.
*/
Polygon_soup_to_polygon_mesh(const PointRange& points,
const PolygonRange& polygons)
: _points(points),
_polygons(polygons)
PS_to_PM_converter(const PointRange& points,
const PolygonRange& polygons,
const PointMap pm = PointMap())
: m_points(points),
m_polygons(polygons),
m_pm(pm)
{ }
void operator()(PM& pmesh, const bool insert_isolated_vertices = true)
template <typename PolygonMesh, typename VertexPointMap>
void operator()(PolygonMesh& pmesh,
VertexPointMap vpm,
const bool insert_isolated_vertices = true)
{
reserve(pmesh, static_cast<typename boost::graph_traits<PM>::vertices_size_type>(_points.size()),
static_cast<typename boost::graph_traits<PM>::edges_size_type>(2*_polygons.size()),
static_cast<typename boost::graph_traits<PM>::faces_size_type>(_polygons.size()) );
typedef typename boost::graph_traits<PolygonMesh>::vertex_descriptor vertex_descriptor;
Vpmap vpmap = get(CGAL::vertex_point, pmesh);
typedef typename boost::property_traits<VertexPointMap>::value_type PM_Point;
CGAL_static_assertion((std::is_convertible<Point, PM_Point>::value));
reserve(pmesh, static_cast<typename boost::graph_traits<PolygonMesh>::vertices_size_type>(m_points.size()),
static_cast<typename boost::graph_traits<PolygonMesh>::edges_size_type>(2*m_polygons.size()),
static_cast<typename boost::graph_traits<PolygonMesh>::faces_size_type>(m_polygons.size()));
boost::dynamic_bitset<> not_isolated;
if (!insert_isolated_vertices)
if(!insert_isolated_vertices)
{
not_isolated.resize(_points.size());
for (std::size_t i = 0, end = _polygons.size(); i < end; ++i)
not_isolated.resize(m_points.size());
for(std::size_t i = 0, end = m_polygons.size(); i < end; ++i)
{
const Polygon& polygon = _polygons[i];
const Polygon& polygon = m_polygons[i];
const std::size_t size = polygon.size();
for (std::size_t j = 0; j < size; ++j)
for(std::size_t j = 0; j < size; ++j)
not_isolated.set(polygon[j], true);
}
}
std::vector<vertex_descriptor> vertices(_points.size());
for (std::size_t i = 0, end = _points.size(); i < end; ++i)
std::vector<vertex_descriptor> vertices(m_points.size());
for(std::size_t i = 0, end = m_points.size(); i < end; ++i)
{
if (!insert_isolated_vertices && !not_isolated.test(i))
if(!insert_isolated_vertices && !not_isolated.test(i))
continue;
Point_3 pi(_points[i][0], _points[i][1], _points[i][2]);
vertices[i] = add_vertex(pmesh);
put(vpmap, vertices[i], pi);
PM_Point pi(get(m_pm, m_points[i]));
put(vpm, vertices[i], pi);
}
for (std::size_t i = 0, end = _polygons.size(); i < end; ++i)
for(std::size_t i = 0, end = m_polygons.size(); i < end; ++i)
{
const Polygon& polygon = _polygons[i];
const Polygon& polygon = m_polygons[i];
const std::size_t size = polygon.size();
std::vector<vertex_descriptor> vr(size); //vertex range
vr.resize(size);
for (std::size_t j = 0; j < size; ++j)
for(std::size_t j = 0; j < size; ++j)
vr[j] = vertices[polygon[j] ];
CGAL_assertion_code(face_descriptor fd = )
CGAL::Euler::add_face(vr, pmesh);
CGAL_assertion(fd != boost::graph_traits<PM>::null_face());
CGAL_assertion_code(typename boost::graph_traits<PolygonMesh>::face_descriptor fd =)
CGAL::Euler::add_face(vr, pmesh);
CGAL_assertion(fd != boost::graph_traits<PolygonMesh>::null_face());
}
}
template <typename PolygonMesh>
void operator()(PolygonMesh& pmesh,
const bool insert_isolated_vertices = true)
{
return operator()(pmesh, get(CGAL::vertex_point, pmesh), insert_isolated_vertices);
}
private:
const PointRange& m_points;
const PolygonRange& m_polygons;
const PointMap m_pm;
};
}//end namespace internal
} // namespace internal
/**
* \ingroup PMP_repairing_grp
*
* returns `true` if the soup of polygons defines a valid polygon
* mesh that can be handled by
* `CGAL::Polygon_mesh_processing::polygon_soup_to_polygon_mesh()`.
* It checks that each edge has at most two incident faces and such an edge
* is visited in opposite direction along the two face boundaries,
* no polygon has twice the same vertex,
* and the polygon soup describes a manifold surface.
* This function does not require a range of points as an argument
* since the check is purely topological. To each vertex of the mesh
* is associated an index that is used in the description of the
* boundaries of the polygons provided in `polygons`.
*
* @tparam PolygonRange a model of the concept `RandomAccessContainer`
* whose value_type is a model of the concept `RandomAccessContainer`
* whose value_type is `std::size_t`.
*
* @param polygons each element in the range describes a polygon
* using the indices of the vertices.
*
* @sa `orient_polygon_soup()`
*/
template<class PolygonRange>
bool is_polygon_soup_a_polygon_mesh(const PolygonRange& polygons)
/**
* \ingroup PMP_repairing_grp
*
* returns `true` if the soup of polygons defines a valid polygon
* mesh that can be handled by
* `CGAL::Polygon_mesh_processing::polygon_soup_to_polygon_mesh()`.
* It checks that each edge has at most two incident faces and such an edge
* is visited in opposite direction along the two face boundaries,
* no polygon has twice the same vertex,
* and the polygon soup describes a manifold surface.
* This function does not require a range of points as an argument
* since the check is purely topological. To each vertex of the mesh
* is associated an index that is used in the description of the
* boundaries of the polygons provided in `polygons`.
*
* @tparam PolygonRange a model of the concept `RandomAccessContainer`
* whose value_type is a model of the concept `RandomAccessContainer`
* whose value_type is `std::size_t`.
*
* @param polygons each element in the range describes a polygon
* using the indices of the vertices.
*
* @sa `orient_polygon_soup()`
*/
template<typename PolygonRange>
bool is_polygon_soup_a_polygon_mesh(const PolygonRange& polygons)
{
typedef typename boost::range_value<PolygonRange>::type Polygon;
typedef typename boost::range_value<Polygon>::type V_ID;
if(boost::begin(polygons) == boost::end(polygons))
return true;
//check there is no duplicated ordered edge, and
//check there is no polygon with twice the same vertex
std::set<std::pair<V_ID, V_ID> > edge_set;
V_ID max_id = 0;
for(const Polygon& polygon : polygons)
{
typedef typename boost::range_value<
typename boost::range_value<
PolygonRange>::type >::type V_ID;
typedef typename boost::range_value<
PolygonRange>::type Polygon;
if(boost::begin(polygons) == boost::end(polygons)){
return true;
}
//check there is no duplicated ordered edge, and
//check there is no polygon with twice the same vertex
std::set< std::pair<V_ID, V_ID> > edge_set;
V_ID max_id=0;
for(const Polygon& polygon : polygons)
{
std::size_t nb_edges = boost::size(polygon);
if (nb_edges<3) return false;
std::set<V_ID> polygon_vertices;
V_ID prev= *std::prev(boost::end(polygon));
for(V_ID id : polygon)
{
if (max_id<id) max_id=id;
if (! edge_set.insert(std::pair<V_ID, V_ID>(prev,id)).second )
return false;
else
prev=id;
if (!polygon_vertices.insert(id).second)
return false;//vertex met twice in the same polygon
}
}
//check manifoldness
typedef std::vector<V_ID> PointRange;
typedef internal::Polygon_soup_orienter<PointRange, PolygonRange> Orienter;
typename Orienter::Edge_map edges(max_id+1);
typename Orienter::Marked_edges marked_edges;
Orienter::fill_edge_map(edges, marked_edges, polygons);
//returns false if duplication is necessary
if (!marked_edges.empty())
std::size_t nb_edges = boost::size(polygon);
if(nb_edges < 3)
return false;
return Orienter::has_singular_vertices(static_cast<std::size_t>(max_id+1),polygons,edges,marked_edges);
std::set<V_ID> polygon_vertices;
V_ID prev = *std::prev(boost::end(polygon));
for(V_ID id : polygon)
{
if(max_id<id)
max_id = id;
if(! edge_set.insert(std::pair<V_ID, V_ID>(prev,id)).second)
return false;
else
prev = id;
if(!polygon_vertices.insert(id).second)
return false; // vertex met twice in the same polygon
}
}
/**
* \ingroup PMP_repairing_grp
* builds a polygon mesh from a soup of polygons.
* @pre the input polygon soup describes a consistently oriented
* polygon mesh.
*
* @tparam PolygonMesh a model of `MutableFaceGraph` with an internal point
* property map
* @tparam PointRange a model of the concepts `RandomAccessContainer` and
* `BackInsertionSequence` whose value type is the point type
* @tparam PolygonRange a model of the concept `RandomAccessContainer` whose
* `value_type` is a model of the concept `RandomAccessContainer` whose `value_type` is `std::size_t`.
//check manifoldness
typedef std::vector<V_ID> PointRange;
typedef internal::Polygon_soup_orienter<PointRange, PolygonRange> Orienter;
*
* @param points points of the soup of polygons
* @param polygons each element in the vector describes a polygon using the index of the points in `points`
* @param out the polygon mesh to be built
*
* @pre \link CGAL::Polygon_mesh_processing::is_polygon_soup_a_polygon_mesh()
* CGAL::Polygon_mesh_processing::is_polygon_soup_a_polygon_mesh(polygons) \endlink
*
* \sa `CGAL::Polygon_mesh_processing::orient_polygon_soup()`
* \sa `CGAL::Polygon_mesh_processing::is_polygon_soup_a_polygon_mesh()`
*
*/
template<class PolygonMesh, class PointRange, class PolygonRange>
void polygon_soup_to_polygon_mesh(
const PointRange& points,
const PolygonRange& polygons,
PolygonMesh& out)
{
CGAL_precondition_msg(is_polygon_soup_a_polygon_mesh(polygons),
"Input soup needs to be a polygon mesh!");
typename Orienter::Edge_map edges(max_id+1);
typename Orienter::Marked_edges marked_edges;
Orienter::fill_edge_map(edges, marked_edges, polygons);
internal::Polygon_soup_to_polygon_mesh<PolygonMesh, PointRange, PolygonRange>
converter(points, polygons);
converter(out);
}
//returns false if duplication is necessary
if(!marked_edges.empty())
return false;
}//end namespace Polygon_mesh_processing
return Orienter::has_singular_vertices(static_cast<std::size_t>(max_id+1),polygons,edges,marked_edges);
}
}// end namespace CGAL
/**
* \ingroup PMP_repairing_grp
* builds a polygon mesh from a soup of polygons.
*
* @pre the input polygon soup describes a consistently oriented
* polygon mesh. This can be checked using the function
* \link CGAL::Polygon_mesh_processing::is_polygon_soup_a_polygon_mesh()
* `CGAL::Polygon_mesh_processing::is_polygon_soup_a_polygon_mesh(polygons)` \endlink.
*
* @tparam PolygonMesh a model of `MutableFaceGraph`
* @tparam PointRange a model of the concept `RandomAccessContainer`
* whose value type is the point type
* @tparam PolygonRange a model of the concept `RandomAccessContainer` whose
* value type is a model of the concept `RandomAccessContainer` whose value type is `std::size_t`
* @tparam NamedParameters_PS a sequence of \ref pmp_namedparameters "Named Parameters"
* @tparam NamedParameters_PM a sequence of \ref pmp_namedparameters "Named Parameters"
*
* @param points points of the soup of polygons
* @param polygons each element in the vector describes a polygon using the indices of the points in `points`
* @param out the polygon mesh to be built
* @param np_ps optional sequence of \ref pmp_namedparameters "Named Parameters" among the ones listed below
*
* \cgalNamedParamsBegin
* \cgalParamBegin{point_map}
* a model of `ReadablePropertyMap` whose value type is a point type convertible to the point type
* of the vertex point map associated to the polygon mesh. If this parameter is omitted, `CGAL::Identity_property_map` is used.
* \cgalParamEnd
* \cgalNamedParamsEnd
*
* @param np_pm optional sequence of \ref pmp_namedparameters "Named Parameters" among the ones listed below
*
* \cgalNamedParamsBegin
* \cgalParamBegin{vertex_point_map}
* a model of `WritablePropertyMap`, the property map with the points associated to the vertices of `out`.
* If this parameter is omitted, an internal property map for
* `CGAL::vertex_point_t` must be available in `PolygonMesh`.
* \cgalParamEnd
* \cgalNamedParamsEnd
*
* \sa `CGAL::Polygon_mesh_processing::orient_polygon_soup()`
* \sa `CGAL::Polygon_mesh_processing::is_polygon_soup_a_polygon_mesh()`
* \sa `CGAL::Polygon_mesh_processing::polygon_mesh_to_polygon_soup()`
*/
template<typename PolygonMesh,
typename PointRange, typename PolygonRange,
typename NamedParameters_PS, typename NamedParameters_PM>
void polygon_soup_to_polygon_mesh(const PointRange& points,
const PolygonRange& polygons,
PolygonMesh& out,
const NamedParameters_PS& np_ps,
const NamedParameters_PM& np_pm)
{
CGAL_precondition_msg(is_polygon_soup_a_polygon_mesh(polygons),
"Input soup needs to define a valid polygon mesh! See is_polygon_soup_a_polygon_mesh() for further information.");
#include <CGAL/enable_warnings.h>
using parameters::choose_parameter;
using parameters::get_parameter;
#endif // CGAL_POLYGON_MESH_PROCESSING_POLYGON_SOUP_TO_POLYGON_MESH
typedef typename CGAL::GetPointMap<PointRange, NamedParameters_PS>::const_type Point_map;
Point_map pm = choose_parameter<Point_map>(get_parameter(np_ps, internal_np::point_map));
typedef typename CGAL::GetVertexPointMap<PolygonMesh, NamedParameters_PM>::type Vertex_point_map;
Vertex_point_map vpm = choose_parameter(get_parameter(np_pm, internal_np::vertex_point),
get_property_map(CGAL::vertex_point, out));
internal::PS_to_PM_converter<PointRange, PolygonRange, Point_map> converter(points, polygons, pm);
converter(out, vpm);
}
template<typename PolygonMesh, typename PointRange, typename PolygonRange>
void polygon_soup_to_polygon_mesh(const PointRange& points,
const PolygonRange& polygons,
PolygonMesh& out)
{
return polygon_soup_to_polygon_mesh(points, polygons, out, parameters::all_default(), parameters::all_default());
}
} // namespace Polygon_mesh_processing
} // namespace CGAL
#endif // CGAL_POLYGON_MESH_PROCESSING_POLYGON_SOUP_TO_POLYGON_MESH_H
@@ -104,7 +104,7 @@ namespace internal {
*
* @tparam VertexRange model of `Range`, holding
* vertices of type `boost::graph_traits<TriangleMesh>::%vertex_descriptor`.
* Its iterator type is `ForwardIterator`.
* Its iterator type is `ForwardIterator`.
* @tparam TriangleMesh model of `MutableFaceGraph`.
* @tparam NamedParameters a sequence of \ref pmp_namedparameters "Named Parameters"
*
@@ -154,11 +154,11 @@ void random_perturbation(VertexRange vertices
#endif
typedef typename GetGeomTraits<PM, NamedParameters>::type GT;
GT gt = choose_parameter(get_parameter(np, internal_np::geom_traits), GT());
GT gt = choose_parameter<GT>(get_parameter(np, internal_np::geom_traits));
typedef typename GetVertexPointMap<PM, NamedParameters>::type VPMap;
VPMap vpmap = choose_parameter(get_parameter(np, internal_np::vertex_point),
get_property_map(vertex_point, tmesh));
get_property_map(vertex_point, tmesh));
typedef typename internal_np::Lookup_named_param_def <
internal_np::vertex_is_constrained_t,
@@ -6,7 +6,7 @@
// $URL$
// $Id$
// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
//
//
//
// Author(s) : Ilker O. Yaz
@@ -41,9 +41,7 @@ namespace Polygon_mesh_processing {
* The descriptor types `boost::graph_traits<PolygonMesh>::%face_descriptor`
* and `boost::graph_traits<PolygonMesh>::%halfedge_descriptor` must be
* models of `Hashable`.
* If `PolygonMesh` has an internal property map for `CGAL::face_index_t`,
* and no `face_index_map` is given
* as a named parameter, then the internal one must be initialized
*
* @tparam FaceRange range of `boost::graph_traits<PolygonMesh>::%face_descriptor`,
model of `Range`. Its iterator type is `ForwardIterator`.
* @tparam NamedParameters a sequence of \ref pmp_namedparameters "Named Parameters"
@@ -139,6 +137,7 @@ void isotropic_remeshing(const FaceRange& faces
typedef PolygonMesh PM;
typedef typename boost::graph_traits<PM>::vertex_descriptor vertex_descriptor;
typedef typename boost::graph_traits<PM>::edge_descriptor edge_descriptor;
using parameters::get_parameter;
using parameters::choose_parameter;
@@ -154,15 +153,14 @@ void isotropic_remeshing(const FaceRange& faces
parameters::is_default_parameter(get_parameter(np, internal_np::projection_functor));
typedef typename GetGeomTraits<PM, NamedParameters>::type GT;
GT gt = choose_parameter(get_parameter(np, internal_np::geom_traits), GT());
GT gt = choose_parameter<GT>(get_parameter(np, internal_np::geom_traits));
typedef typename GetVertexPointMap<PM, NamedParameters>::type VPMap;
VPMap vpmap = choose_parameter(get_parameter(np, internal_np::vertex_point),
get_property_map(vertex_point, pmesh));
get_property_map(vertex_point, pmesh));
typedef typename GetFaceIndexMap<PM, NamedParameters>::type FIMap;
FIMap fimap = choose_parameter(get_parameter(np, internal_np::face_index),
get_property_map(face_index, pmesh));
typedef typename GetInitializedFaceIndexMap<PolygonMesh, NamedParameters>::type FIMap;
FIMap fimap = CGAL::get_initialized_face_index_map(pmesh, np);
typedef typename internal_np::Lookup_named_param_def <
internal_np::edge_is_constrained_t,
@@ -332,15 +330,14 @@ void split_long_edges(const EdgeRange& edges
using parameters::get_parameter;
typedef typename GetGeomTraits<PM, NamedParameters>::type GT;
GT gt = choose_parameter(get_parameter(np, internal_np::geom_traits), GT());
GT gt = choose_parameter<GT>(get_parameter(np, internal_np::geom_traits));
typedef typename GetVertexPointMap<PM, NamedParameters>::type VPMap;
VPMap vpmap = choose_parameter(get_parameter(np, internal_np::vertex_point),
get_property_map(vertex_point, pmesh));
get_property_map(vertex_point, pmesh));
typedef typename GetFaceIndexMap<PM, NamedParameters>::type FIMap;
FIMap fimap = choose_parameter(get_parameter(np, internal_np::face_index),
get_property_map(face_index, pmesh));
typedef typename GetInitializedFaceIndexMap<PolygonMesh, NamedParameters>::type FIMap;
FIMap fimap = CGAL::get_initialized_face_index_map(pmesh, np);
typedef typename internal_np::Lookup_named_param_def <
internal_np::edge_is_constrained_t,
@@ -349,7 +346,7 @@ void split_long_edges(const EdgeRange& edges
> ::type ECMap;
ECMap ecmap = choose_parameter(get_parameter(np, internal_np::edge_is_constrained),
Constant_property_map<edge_descriptor, bool>(false));
typename internal::Incremental_remesher<PM, VPMap, GT, ECMap,
Constant_property_map<vertex_descriptor, bool>, // no constraint pmap
internal::Connected_components_pmap<PM, FIMap>,
File diff suppressed because it is too large Load Diff
@@ -513,7 +513,7 @@ std::size_t merge_duplicate_points_in_polygon_soup(PointRange& points,
using parameters::choose_parameter;
typedef typename internal::GetPolygonGeomTraits<PointRange, PolygonRange, NamedParameters>::type Traits;
Traits traits = choose_parameter(get_parameter(np, internal_np::geom_traits), Traits());
Traits traits = choose_parameter<Traits>(get_parameter(np, internal_np::geom_traits));
typedef typename Traits::Less_xyz_3 Less_xyz_3;
@@ -661,7 +661,7 @@ Polygon construct_canonical_polygon(const PointRange& points,
reversed = false;
return polygon;
}
#ifdef CGAL_PMP_REPAIR_POLYGON_SOUP_VERBOSE_PP
std::cout << "Input polygon:";
@@ -915,7 +915,7 @@ std::size_t merge_duplicate_polygons_in_polygon_soup(const PointRange& points,
#endif
typedef typename internal::GetPolygonGeomTraits<PointRange, PolygonRange, NamedParameters>::type Traits;
Traits traits = choose_parameter(get_parameter(np, internal_np::geom_traits), Traits());
Traits traits = choose_parameter<Traits>(get_parameter(np, internal_np::geom_traits));
std::vector<std::vector<P_ID> > all_duplicate_polygons;
internal::collect_duplicate_polygons(points, polygons, std::back_inserter(all_duplicate_polygons), traits, same_orientation);
@@ -1062,7 +1062,7 @@ void repair_polygon_soup(PointRange& points,
using parameters::choose_parameter;
typedef typename internal::GetPolygonGeomTraits<PointRange, PolygonRange, NamedParameters>::type Traits;
Traits traits = choose_parameter(get_parameter(np, internal_np::geom_traits), Traits());
Traits traits = choose_parameter<Traits>(get_parameter(np, internal_np::geom_traits));
#ifdef CGAL_PMP_REPAIR_POLYGON_SOUP_VERBOSE
std::cout << "Repairing soup with " << points.size() << " points and " << polygons.size() << " polygons" << std::endl;
@@ -220,7 +220,7 @@ self_intersections_impl(const FaceRange& face_range,
typedef CGAL::Box_intersection_d::Box_with_info_d<double, 3, face_descriptor, Box_policy> Box;
typedef typename GetGeomTraits<TM, NamedParameters>::type GT;
GT gt = choose_parameter(get_parameter(np, internal_np::geom_traits), GT());
GT gt = choose_parameter<GT>(get_parameter(np, internal_np::geom_traits));
typedef typename GetVertexPointMap<TM, NamedParameters>::const_type VPM;
VPM vpmap = choose_parameter(get_parameter(np, internal_np::vertex_point),
@@ -23,6 +23,7 @@
#include <CGAL/boost/graph/iterator.h>
#include <CGAL/boost/graph/helpers.h>
#include <boost/range/has_range_iterator.hpp>
#include <boost/graph/graph_traits.hpp>
#include <limits>
@@ -70,10 +71,10 @@ bool is_degenerate_edge(typename boost::graph_traits<PolygonMesh>::edge_descript
typedef typename GetVertexPointMap<PolygonMesh, NamedParameters>::const_type VertexPointMap;
VertexPointMap vpmap = choose_parameter(get_parameter(np, internal_np::vertex_point),
get_const_property_map(vertex_point, pm));
get_const_property_map(vertex_point, pm));
typedef typename GetGeomTraits<PolygonMesh, NamedParameters>::type Traits;
Traits traits = choose_parameter(get_parameter(np, internal_np::geom_traits), Traits());
Traits traits = choose_parameter<Traits>(get_parameter(np, internal_np::geom_traits));
return traits.equal_3_object()(get(vpmap, source(e, pm)), get(vpmap, target(e, pm)));
}
@@ -85,6 +86,81 @@ bool is_degenerate_edge(typename boost::graph_traits<PolygonMesh>::edge_descript
return is_degenerate_edge(e, pm, parameters::all_default());
}
/// \ingroup PMP_repairing_grp
/// collects the degenerate edges within a given range of edges.
///
/// @tparam EdgeRange a model of `Range` with value type `boost::graph_traits<TriangleMesh>::%edge_descriptor`
/// @tparam TriangleMesh a model of `HalfedgeGraph`
/// @tparam NamedParameters a sequence of \ref pmp_namedparameters "Named Parameters"
///
/// @param edges a subset of edges of `tm`
/// @param tm a triangle mesh
/// @param out an output iterator in which the degenerate edges are written
/// @param np optional \ref pmp_namedparameters "Named Parameters" described below
///
/// \cgalNamedParamsBegin
/// \cgalParamBegin{vertex_point_map} the property map with the points associated to the vertices of `tm`.
/// The type of this map is model of `ReadWritePropertyMap`.
/// If this parameter is omitted, an internal property map for
/// `CGAL::vertex_point_t` should be available in `TriangleMesh`
/// \cgalParamEnd
/// \cgalParamBegin{geom_traits} a geometric traits class instance.
/// The traits class must provide the nested type `Point_3`,
/// and the nested functor `Equal_3` to check whether two points are identical.
/// \cgalParamEnd
/// \cgalNamedParamsEnd
template <typename EdgeRange, typename TriangleMesh, typename OutputIterator, typename NamedParameters>
OutputIterator degenerate_edges(const EdgeRange& edges,
const TriangleMesh& tm,
OutputIterator out,
const NamedParameters& np)
{
typedef typename boost::graph_traits<TriangleMesh>::edge_descriptor edge_descriptor;
for(edge_descriptor ed : edges)
if(is_degenerate_edge(ed, tm, np))
*out++ = ed;
return out;
}
template <typename EdgeRange, typename TriangleMesh, typename OutputIterator>
OutputIterator degenerate_edges(const EdgeRange& edges,
const TriangleMesh& tm,
OutputIterator out,
typename boost::enable_if<
typename boost::has_range_iterator<EdgeRange>
>::type* = 0)
{
return degenerate_edges(edges, tm, out, CGAL::parameters::all_default());
}
/// \ingroup PMP_repairing_grp
/// calls the function `degenerate_edges()` with the range: `edges(tm)`.
///
/// See above for the comprehensive description of the parameters.
///
template <typename TriangleMesh, typename OutputIterator, typename NamedParameters>
OutputIterator degenerate_edges(const TriangleMesh& tm,
OutputIterator out,
const NamedParameters& np
#ifndef DOXYGEN_RUNNING
,
typename boost::disable_if<
boost::has_range_iterator<TriangleMesh>
>::type* = 0
#endif
)
{
return degenerate_edges(edges(tm), tm, out, np);
}
template <typename TriangleMesh, typename OutputIterator>
OutputIterator degenerate_edges(const TriangleMesh& tm, OutputIterator out)
{
return degenerate_edges(edges(tm), tm, out, CGAL::parameters::all_default());
}
/// \ingroup PMP_repairing_grp
/// checks whether a triangle face is degenerate.
/// A triangle face is considered degenerate if the geometric positions of its vertices are collinear.
@@ -123,10 +199,10 @@ bool is_degenerate_triangle_face(typename boost::graph_traits<TriangleMesh>::fac
typedef typename GetVertexPointMap<TriangleMesh, NamedParameters>::const_type VertexPointMap;
VertexPointMap vpmap = choose_parameter(get_parameter(np, internal_np::vertex_point),
get_const_property_map(vertex_point, tm));
get_const_property_map(vertex_point, tm));
typedef typename GetGeomTraits<TriangleMesh, NamedParameters>::type Traits;
Traits traits = choose_parameter(get_parameter(np, internal_np::geom_traits), Traits());
Traits traits = choose_parameter<Traits>(get_parameter(np, internal_np::geom_traits));
typename boost::graph_traits<TriangleMesh>::halfedge_descriptor h = halfedge(f, tm);
@@ -142,6 +218,83 @@ bool is_degenerate_triangle_face(typename boost::graph_traits<TriangleMesh>::fac
return CGAL::Polygon_mesh_processing::is_degenerate_triangle_face(f, tm, parameters::all_default());
}
/// \ingroup PMP_repairing_grp
/// collects the degenerate faces within a given range of faces.
///
/// @tparam FaceRange a model of `Range` with value type `boost::graph_traits<TriangleMesh>::%face_descriptor`
/// @tparam TriangleMesh a model of `FaceGraph`
/// @tparam NamedParameters a sequence of \ref pmp_namedparameters "Named Parameters"
///
/// @param faces a subset of faces of `tm`
/// @param tm a triangle mesh
/// @param out an output iterator in which the degenerate faces are put
/// @param np optional \ref pmp_namedparameters "Named Parameters" described below
///
/// \cgalNamedParamsBegin
/// \cgalParamBegin{vertex_point_map} the property map with the points associated to the vertices of `tm`.
/// The type of this map is model of `ReadWritePropertyMap`.
/// If this parameter is omitted, an internal property map for
/// `CGAL::vertex_point_t` should be available in `TriangleMesh`
/// \cgalParamEnd
/// \cgalParamBegin{geom_traits} a geometric traits class instance.
/// The traits class must provide the nested functor `Collinear_3`
/// to check whether three points are collinear.
/// \cgalParamEnd
/// \cgalNamedParamsEnd
///
template <typename FaceRange, typename TriangleMesh, typename OutputIterator, typename NamedParameters>
OutputIterator degenerate_faces(const FaceRange& faces,
const TriangleMesh& tm,
OutputIterator out,
const NamedParameters& np)
{
typedef typename boost::graph_traits<TriangleMesh>::face_descriptor face_descriptor;
for(face_descriptor fd : faces)
{
if(is_degenerate_triangle_face(fd, tm, np))
*out++ = fd;
}
return out;
}
template <typename FaceRange, typename TriangleMesh, typename OutputIterator>
OutputIterator degenerate_faces(const FaceRange& faces,
const TriangleMesh& tm,
OutputIterator out,
typename boost::enable_if<
boost::has_range_iterator<FaceRange>
>::type* = 0)
{
return degenerate_faces(faces, tm, out, CGAL::parameters::all_default());
}
/// \ingroup PMP_repairing_grp
/// calls the function `degenerate_faces()` with the range: `faces(tm)`.
///
/// See above for the comprehensive description of the parameters.
///
template <typename TriangleMesh, typename OutputIterator, typename NamedParameters>
OutputIterator degenerate_faces(const TriangleMesh& tm,
OutputIterator out,
const NamedParameters& np
#ifndef DOXYGEN_RUNNING
,
typename boost::disable_if<
boost::has_range_iterator<TriangleMesh>
>::type* = 0
#endif
)
{
return degenerate_faces(faces(tm), tm, out, np);
}
template <typename TriangleMesh, typename OutputIterator>
OutputIterator degenerate_faces(const TriangleMesh& tm, OutputIterator out)
{
return degenerate_faces(faces(tm), tm, out, CGAL::parameters::all_default());
}
/// \ingroup PMP_repairing_grp
/// checks whether a triangle face is needle.
/// A triangle is said to be a <i>needle</i> if its longest edge is much longer than its shortest edge.
@@ -184,10 +337,10 @@ is_needle_triangle_face(typename boost::graph_traits<TriangleMesh>::face_descrip
typedef typename GetVertexPointMap<TriangleMesh, NamedParameters>::const_type VertexPointMap;
VertexPointMap vpmap = choose_parameter(get_parameter(np, internal_np::vertex_point),
get_const_property_map(vertex_point, tm));
get_const_property_map(vertex_point, tm));
typedef typename GetGeomTraits<TriangleMesh, NamedParameters>::type Traits;
Traits traits = choose_parameter(get_parameter(np, internal_np::geom_traits), Traits());
Traits traits = choose_parameter<Traits>(get_parameter(np, internal_np::geom_traits));
typedef typename Traits::FT FT;
@@ -287,10 +440,10 @@ is_cap_triangle_face(typename boost::graph_traits<TriangleMesh>::face_descriptor
typedef typename GetVertexPointMap<TriangleMesh, NamedParameters>::const_type VertexPointMap;
VertexPointMap vpmap = choose_parameter(get_parameter(np, internal_np::vertex_point),
get_const_property_map(vertex_point, tm));
get_const_property_map(vertex_point, tm));
typedef typename GetGeomTraits<TriangleMesh, NamedParameters>::type Traits;
Traits traits = choose_parameter(get_parameter(np, internal_np::geom_traits), Traits());
Traits traits = choose_parameter<Traits>(get_parameter(np, internal_np::geom_traits));
typedef typename Traits::FT FT;
typedef typename Traits::Vector_3 Vector_3;
@@ -150,10 +150,9 @@ void smooth_mesh(const FaceRange& faces,
using parameters::get_parameter;
// named parameters
GeomTraits gt = choose_parameter(get_parameter(np, internal_np::geom_traits),
GeomTraits());
GeomTraits gt = choose_parameter<GeomTraits>(get_parameter(np, internal_np::geom_traits));
VertexPointMap vpmap = choose_parameter(get_parameter(np, internal_np::vertex_point),
get_property_map(CGAL::vertex_point, tmesh));
get_property_map(CGAL::vertex_point, tmesh));
const bool use_angle_smoothing = choose_parameter(get_parameter(np, internal_np::use_angle_smoothing), true);
bool use_area_smoothing = choose_parameter(get_parameter(np, internal_np::use_area_smoothing), true);
@@ -173,7 +172,7 @@ void smooth_mesh(const FaceRange& faces,
const bool use_Delaunay_flips = choose_parameter(get_parameter(np, internal_np::use_Delaunay_flips), true);
VCMap vcmap = choose_parameter(get_parameter(np, internal_np::vertex_is_constrained),
get(Vertex_property_tag(), tmesh));
get(Vertex_property_tag(), tmesh));
// If it's the default vcmap, manually set everything to false because the dynamic pmap has no default initialization
if((std::is_same<VCMap, Default_VCMap>::value))
@@ -183,7 +182,7 @@ void smooth_mesh(const FaceRange& faces,
}
ECMap ecmap = choose_parameter(get_parameter(np, internal_np::edge_is_constrained),
Constant_property_map<edge_descriptor, bool>(false));
Constant_property_map<edge_descriptor, bool>(false));
// a constrained edge has constrained extremities
for(face_descriptor f : faces)
@@ -235,8 +234,16 @@ void smooth_mesh(const FaceRange& faces,
for(unsigned int i=0; i<nb_iterations; ++i)
{
#ifdef CGAL_PMP_SMOOTHING_DEBUG
std::cout << "Iteration #" << i << std::endl;
#endif
if(use_area_smoothing)
{
#ifdef CGAL_PMP_SMOOTHING_DEBUG
std::cout << "Smooth areas..." << std::endl;
#endif
// First apply area smoothing...
area_smoother.optimize(use_safety_constraints /*check for bad faces*/,
false /*apply moves as soon as they're calculated*/,
@@ -245,7 +252,7 @@ void smooth_mesh(const FaceRange& faces,
{
if(use_safety_constraints && does_self_intersect(tmesh))
{
#ifdef CGAL_PMP_SMOOTHING_VERBOSE
#ifdef CGAL_PMP_SMOOTHING_DEBUG
std::cerr << "Cannot re-project as there are self-intersections in the mesh!\n";
#endif
break;
@@ -261,6 +268,10 @@ void smooth_mesh(const FaceRange& faces,
// ... then angle smoothing
if(use_angle_smoothing)
{
#ifdef CGAL_PMP_SMOOTHING_DEBUG
std::cout << "Smooth angles..." << std::endl;
#endif
angle_smoother.optimize(use_safety_constraints /*check for bad faces*/,
true /*apply all moves at once*/,
use_safety_constraints /*check if the min angle is improved*/);
@@ -269,7 +280,7 @@ void smooth_mesh(const FaceRange& faces,
{
if(use_safety_constraints && does_self_intersect(tmesh))
{
#ifdef CGAL_PMP_SMOOTHING_VERBOSE
#ifdef CGAL_PMP_SMOOTHING_DEBUG
std::cerr << "Can't do re-projection, there are self-intersections in the mesh!\n";
#endif
break;
@@ -97,11 +97,11 @@ void smooth_shape(const FaceRange& faces,
using parameters::choose_parameter;
using parameters::get_parameter;
GeomTraits gt = choose_parameter(get_parameter(np, internal_np::geom_traits), GeomTraits());
GeomTraits gt = choose_parameter<GeomTraits>(get_parameter(np, internal_np::geom_traits));
VertexPointMap vpmap = choose_parameter(get_parameter(np, internal_np::vertex_point),
get_property_map(CGAL::vertex_point, tmesh));
get_property_map(CGAL::vertex_point, tmesh));
VCMap vcmap = choose_parameter(get_parameter(np, internal_np::vertex_is_constrained),
Constant_property_map<vertex_descriptor, bool>(false));
Constant_property_map<vertex_descriptor, bool>(false));
const unsigned int nb_iterations = choose_parameter(get_parameter(np, internal_np::number_of_iterations), 1);
#if defined(CGAL_EIGEN3_ENABLED)
@@ -132,7 +132,7 @@ void smooth_shape(const FaceRange& faces,
typedef typename Sparse_solver::Matrix Eigen_matrix;
typedef typename Sparse_solver::Vector Eigen_vector;
Sparse_solver solver = choose_parameter(get_parameter(np, internal_np::sparse_linear_solver), Default_solver());
Sparse_solver solver = choose_parameter<Default_solver>(get_parameter(np, internal_np::sparse_linear_solver));
std::size_t n = vertices(tmesh).size();
Eigen_matrix A(n, n);
@@ -149,7 +149,7 @@ void smooth_shape(const FaceRange& faces,
for(unsigned int iter=0; iter<nb_iterations; ++iter)
{
#ifdef CGAL_PMP_SMOOTHING_VERBOSE
#ifdef CGAL_PMP_SMOOTHING_DEBUG
std::cout << "iteration #" << iter << std::endl;
#endif
@@ -161,7 +161,7 @@ void smooth_shape(const FaceRange& faces,
}
else
{
#ifdef CGAL_PMP_SMOOTHING_VERBOSE
#ifdef CGAL_PMP_SMOOTHING_DEBUG
std::cerr << "Failed to solve system!" << std::endl;
#endif
break;
@@ -234,10 +234,8 @@ collect_duplicated_stitchable_boundary_edges(PolygonMesh& pmesh,
if(per_cc)
{
cc = get(Face_property_tag(), pmesh);
typedef typename GetFaceIndexMap<PolygonMesh, CGAL_PMP_NP_CLASS>::const_type FIMap;
FIMap fim = choose_parameter(get_parameter(np, internal_np::face_index),
get_const_property_map(face_index, pmesh));
num_component = num_component_wrapper(pmesh, cc, fim);
num_component = num_component_wrapper(pmesh, cc, CGAL::get_initialized_face_index_map(pmesh, np));
border_edges_per_cc.resize(num_component);
}
@@ -927,7 +925,7 @@ std::size_t stitch_borders(PolygonMesh& pmesh,
template <typename PolygonMesh,
typename HalfedgePairsRange>
std::size_t stitch_borders(PolygonMesh& pmesh,
const HalfedgePairsRange& hedge_pairs_to_stitch)
const HalfedgePairsRange& hedge_pairs_to_stitch)
{
return stitch_borders(pmesh, hedge_pairs_to_stitch, CGAL::parameters::all_default());
}
@@ -946,16 +944,18 @@ std::size_t stitch_borders(PolygonMesh& pmesh,
/// @param np optional sequence of \ref pmp_namedparameters "Named Parameters" among the ones listed below
///
/// \cgalNamedParamsBegin
/// \cgalParamBegin{vertex_point_map} the property map with the points associated to the vertices of `pmesh`.
/// If this parameter is omitted, an internal property map for
/// `CGAL::vertex_point_t` must be available in `PolygonMesh`.\cgalParamEnd
/// \cgalParamBegin{apply_per_connected_component}
/// specifies if the borders should only be stitched inside their own connected component.
/// In that case, a property map for `CGAL::face_index_t` should be either available as an internal property map
/// to `pmesh` or provided as the \ref pmp_namedparameters "Named Parameter" `face_index_map`. If this is not the case,
/// a default map will be created on the fly.
/// Default value is `false`.\cgalParamEnd
/// \cgalParamBegin{face_index_map} a property map containing the index of each face of `pmesh` \cgalParamEnd
/// \cgalParamBegin{vertex_point_map}
/// the property map with the points associated to the vertices of `pmesh`.
/// If this parameter is omitted, an internal property map for
/// `CGAL::vertex_point_t` must be available in `PolygonMesh`.
/// \cgalParamEnd
/// \cgalParamBegin{apply_per_connected_component}
/// specifies if the borders should only be stitched inside their own connected component.
/// Default value is `false`.
/// \cgalParamEnd
/// \cgalParamBegin{face_index_map}
/// a property map containing for each face of `pmesh` a unique index between `0` and `num_faces(pmesh)-1`
/// \cgalParamEnd
/// \cgalNamedParamsEnd
///
/// @return the number of pairs of halfedges that were stitched.
@@ -21,33 +21,33 @@ namespace Polygon_mesh_processing{
/**
* \ingroup PkgPolygonMeshProcessingRef
* applies a transformation to every vertex of a `PolygonMesh`.
*
*
* @tparam Transformation a functor that has an `operator()(Point_3)`, with `Point_3`
* the `value_type` of `vertex_point_map` (see below). Such a functor can be
* `CGAL::Aff_transformation_3` for example.
* @tparam PolygonMesh a model of `VertexListGraph`
* @tparam NamedParameters a sequence of \ref pmp_namedparameters "Named Parameters"
*
*
* @param transformation the transformation functor to apply to the points of `mesh`.
* @param mesh the `PolygonMesh` to transform.
* @param np optional sequence of \ref pmp_namedparameters for `mesh`, among the ones listed below
*
*
* * \cgalNamedParamsBegin
* \cgalParamBegin{vertex_point_map} the property map with the points associated to the vertices of `mesh`.
* If this parameter is omitted, an internal property map for
* `CGAL::vertex_point_t` must be available in `PolygonMesh`\cgalParamEnd
* \cgalNamedParamsEnd
*
*
*/
template<class Transformation, class PolygonMesh,class NamedParameters>
void transform(const Transformation& transformation,
void transform(const Transformation& transformation,
PolygonMesh& mesh,
const NamedParameters& np)
{
typedef typename GetVertexPointMap<PolygonMesh, NamedParameters>::type VPMap;
VPMap vpm = parameters::choose_parameter(parameters::get_parameter(np, internal_np::vertex_point),
get_property_map(vertex_point, mesh));
for(typename boost::graph_traits<PolygonMesh>::vertex_descriptor vd : vertices(mesh))
{
put(vpm, vd, transformation(get(vpm, vd)));
@@ -431,7 +431,7 @@ bool triangulate_face(typename boost::graph_traits<PolygonMesh>::face_descriptor
//Kernel
typedef typename GetGeomTraits<PolygonMesh, NamedParameters>::type Kernel;
Kernel traits = choose_parameter(get_parameter(np, internal_np::geom_traits), Kernel());
Kernel traits = choose_parameter<Kernel>(get_parameter(np, internal_np::geom_traits));
//Option
bool use_cdt = choose_parameter(get_parameter(np, internal_np::use_delaunay_triangulation), true);
@@ -486,7 +486,7 @@ bool triangulate_faces(FaceRange face_range,
//Kernel
typedef typename GetGeomTraits<PolygonMesh, NamedParameters>::type Kernel;
Kernel traits = choose_parameter(get_parameter(np, internal_np::geom_traits), Kernel());
Kernel traits = choose_parameter<Kernel>(get_parameter(np, internal_np::geom_traits));
//Option
bool use_cdt = choose_parameter(get_parameter(np, internal_np::use_delaunay_triangulation), true);
@@ -6,7 +6,7 @@
// $URL$
// $Id$
// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
//
//
// Author(s) : Ilker O. Yaz
#ifndef CGAL_POLYGON_MESH_PROCESSING_TRIANGULATE_HOLE_H
@@ -76,7 +76,7 @@ namespace Polygon_mesh_processing {
@todo Then, insert the holes vertices in the set of possibilities
for connecting vertices together
@todo handle the case where an island is reduced to a point
*/
*/
template<typename PolygonMesh,
typename OutputIterator,
typename NamedParameters>
@@ -89,6 +89,8 @@ namespace Polygon_mesh_processing {
using parameters::choose_parameter;
using parameters::get_parameter;
typedef typename GetGeomTraits<PolygonMesh,NamedParameters>::type GeomTraits;
bool use_dt3 =
#ifdef CGAL_HOLE_FILLING_DO_NOT_USE_DT3
false;
@@ -103,8 +105,7 @@ namespace Polygon_mesh_processing {
out,
choose_parameter(get_parameter(np, internal_np::vertex_point), get_property_map(vertex_point, pmesh)),
use_dt3,
choose_parameter(get_parameter(np, internal_np::geom_traits), typename GetGeomTraits<PolygonMesh,NamedParameters>::type()))
.first;
choose_parameter<GeomTraits>(get_parameter(np, internal_np::geom_traits))).first;
}
template<typename PolygonMesh, typename OutputIterator>
@@ -360,14 +361,14 @@ namespace Polygon_mesh_processing {
typedef typename value_type_traits<OutputIterator>::type OutputIteratorValueType;
CGAL::internal::Tracer_polyline_incomplete<OutputIteratorValueType, OutputIterator, Holes_out>
tracer(out, Holes_out(holes));
typedef typename PointRange1::iterator InIterator;
typedef typename std::iterator_traits<InIterator>::value_type Point;
typedef typename CGAL::Kernel_traits<Point>::Kernel Kernel;
triangulate_hole_polyline(points, third_points, tracer, WC(),
use_dt3,
choose_parameter(get_parameter(np, internal_np::geom_traits),
typename CGAL::Kernel_traits<Point>::Kernel()));
use_dt3,
choose_parameter<Kernel>(get_parameter(np, internal_np::geom_traits)));
CGAL_assertion(holes.empty());
return tracer.out;
@@ -388,7 +389,7 @@ namespace Polygon_mesh_processing {
/*!
\ingroup hole_filling_grp
same as above but the range of third points is omitted. They are not
taken into account in the cost computation that leads the hole filling.
taken into account in the cost computation that leads the hole filling.
*/
template <typename PointRange,
typename OutputIterator,
@@ -225,7 +225,7 @@ public:
Vpm vpm =
parameters::choose_parameter(parameters::get_parameter(np, internal_np::vertex_point),
get_const_property_map(boost::vertex_point, tm) );
get_const_property_map(boost::vertex_point, tm) );
// now add the mesh
std::size_t id = get_id_for_new_mesh();
CGAL_assertion( m_aabb_trees[id] == nullptr );
@@ -529,7 +529,7 @@ public:
Vpm vpm =
parameters::choose_parameter(parameters::get_parameter(np, internal_np::vertex_point),
get_const_property_map(boost::vertex_point, tm) );
get_const_property_map(boost::vertex_point, tm) );
if (maybe_several_cc)
{
@@ -537,11 +537,8 @@ public:
std::vector<std::size_t> cc_ids(num_faces(tm));
// face index map
typedef typename CGAL::GetFaceIndexMap<TriangleMesh, NamedParameters>::type Fid_map;
Fid_map fid_map =
parameters::choose_parameter(parameters::get_parameter(np, internal_np::face_index),
get_const_property_map(boost::face_index, tm));
typedef typename GetInitializedFaceIndexMap<TriangleMesh, NamedParameters>::const_type FaceIndexMap;
FaceIndexMap fid_map = CGAL::get_initialized_face_index_map(tm, np);
std::size_t nb_cc =
Polygon_mesh_processing::connected_components(
@@ -6,7 +6,7 @@
// $URL$
// $Id$
// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
//
//
//
// Author(s) : Sebastien Loriot and Ilker O. Yaz
@@ -28,7 +28,7 @@
namespace CGAL {
/**
/**
* \ingroup PkgPolygonMeshProcessingRef
* This class provides an efficient point location functionality with respect to a domain bounded
* by one or several disjoint closed triangle meshes.
@@ -98,13 +98,15 @@ class Side_of_triangle_mesh
//members
typename GeomTraits::Construct_ray_3 ray_functor;
typename GeomTraits::Construct_vector_3 vector_functor;
mutable const AABB_tree_* tree_ptr;
const TriangleMesh* tm_ptr;
boost::optional<VertexPointMap> opt_vpm;
bool own_tree;
CGAL::Bbox_3 box;
#ifdef CGAL_HAS_THREADS
mutable CGAL_MUTEX tree_mutex;
mutable std::atomic<const AABB_tree_*> atomic_tree_ptr;
#else
mutable const AABB_tree_* tree_ptr;
#endif
public:
@@ -129,10 +131,14 @@ public:
const GeomTraits& gt=GeomTraits())
: ray_functor(gt.construct_ray_3_object())
, vector_functor(gt.construct_vector_3_object())
, tree_ptr(nullptr)
, tm_ptr(&tmesh)
, opt_vpm(vpmap)
, own_tree(true)
#ifdef CGAL_HAS_THREADS
, atomic_tree_ptr(nullptr)
#else
, tree_ptr(nullptr)
#endif
{
CGAL_assertion(CGAL::is_triangle_mesh(tmesh));
CGAL_assertion(CGAL::is_closed(tmesh));
@@ -165,16 +171,24 @@ public:
const GeomTraits& gt = GeomTraits())
: ray_functor(gt.construct_ray_3_object())
, vector_functor(gt.construct_vector_3_object())
, tree_ptr(&tree)
, own_tree(false)
#ifdef CGAL_HAS_THREADS
, atomic_tree_ptr(&tree)
#else
, tree_ptr(&tree)
#endif
{
box = tree.bbox();
}
~Side_of_triangle_mesh()
{
if (own_tree && tree_ptr!=nullptr)
if (own_tree)
#ifdef CGAL_HAS_THREADS
delete atomic_tree_ptr.load();
#else
delete tree_ptr;
#endif
}
public:
@@ -182,7 +196,7 @@ public:
* returns the location of a query point
* @param point the query point to be located with respect to the input
polyhedral surface
* @return
* @return
* - `CGAL::ON_BOUNDED_SIDE` if the point is inside the volume bounded by the input triangle mesh
* - `CGAL::ON_BOUNDARY` if the point is on triangle mesh
* - `CGAL::ON_UNBOUNDED_SIDE` if the point is outside triangle mesh
@@ -200,11 +214,16 @@ public:
}
else
{
#ifdef CGAL_HAS_THREADS
AABB_tree_* tree_ptr =
const_cast<AABB_tree_*>(atomic_tree_ptr.load(std::memory_order_acquire));
#endif
// Lazily build the tree only when needed
if (tree_ptr==nullptr)
{
#ifdef CGAL_HAS_THREADS
CGAL_SCOPED_LOCK(tree_mutex);
tree_ptr = const_cast<AABB_tree_*>(atomic_tree_ptr.load(std::memory_order_relaxed));
#endif
CGAL_assertion(tm_ptr != nullptr && opt_vpm!=boost::none);
if (tree_ptr==nullptr)
@@ -213,9 +232,11 @@ public:
faces(*tm_ptr).second,
*tm_ptr, *opt_vpm);
const_cast<AABB_tree_*>(tree_ptr)->build();
#ifdef CGAL_HAS_THREADS
atomic_tree_ptr.store(tree_ptr, std::memory_order_release);
#endif
}
}
return internal::Point_inside_vertical_ray_cast<GeomTraits, AABB_tree>()(
point, *tree_ptr, ray_functor, vector_functor);
}
@@ -6,7 +6,7 @@
// $URL$
// $Id$
// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
//
//
//
// Author(s) : Jane Tournois
@@ -20,6 +20,7 @@
#include <CGAL/Polygon_mesh_processing/orient_polygon_soup.h>
#include <CGAL/Polygon_mesh_processing/self_intersections.h>
#include <CGAL/Polygon_mesh_processing/polygon_soup_to_polygon_mesh.h>
#include <CGAL/Polygon_mesh_processing/polygon_mesh_to_polygon_soup.h>
#include <CGAL/Polygon_mesh_processing/triangulate_faces.h>
#include <CGAL/Polygon_mesh_processing/triangulate_hole.h>
#include <CGAL/Polygon_mesh_processing/compute_normal.h>
@@ -32,6 +33,8 @@
#include <CGAL/Polygon_mesh_processing/bbox.h>
#include <CGAL/Polygon_mesh_processing/border.h>
#include <CGAL/Polygon_mesh_processing/repair.h>
#include <CGAL/Polygon_mesh_processing/repair_degeneracies.h>
#include <CGAL/Polygon_mesh_processing/repair_self_intersections.h>
#include <CGAL/Polygon_mesh_processing/remesh.h>
#include <CGAL/Polygon_mesh_processing/corefinement.h>
#include <CGAL/Polygon_mesh_processing/detect_features.h>
@@ -45,7 +48,7 @@
#include <CGAL/Polygon_mesh_processing/merge_border_vertices.h>
#include <CGAL/Polygon_mesh_processing/smooth_mesh.h>
#include <CGAL/Polygon_mesh_processing/smooth_shape.h>
#include <CGAL/Polygon_mesh_processing/internal/remove_degeneracies.h>
#include <CGAL/Polygon_mesh_processing/manifoldness.h>
// the named parameter header being not documented the doc is put here for now
#ifdef DOXYGEN_RUNNING