Merge remote-tracking branch 'maxGimeno/PMP-compare_faces_from_meshes-maxGimeno' into gsoc2019-PMPHDist-martinskrodzki

This commit is contained in:
Dmitry Anisimov
2021-04-07 13:54:02 +02:00
357 changed files with 7048 additions and 10260 deletions
@@ -41,6 +41,7 @@
#include <unordered_map>
#include <utility>
#include <vector>
#include <bitset>
namespace CGAL{
namespace Polygon_mesh_processing {
@@ -48,24 +49,6 @@ namespace Polygon_mesh_processing {
namespace internal
{
template <class Geom_traits, class Plane_3, class Point_3>
int
inter_pt_index(int i, int j,
const Plane_3& plane,
std::vector<Point_3>& points,
std::map<std::pair<int,int>, int>& id_map)
{
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)
points.push_back(
typename Geom_traits::Construct_plane_line_intersection_point_3()
(plane, points[i], points[j]));
return res.first->second;
}
template <class Plane_3,
class TriangleMesh,
class NamedParameters>
@@ -106,19 +89,42 @@ clip_to_bbox(const Plane_3& plane,
}};
// description of faces of the bbox
std::array<int, 24> face_indices =
{{ 0, 1, 2, 3,
2, 1, 5, 6,
3, 2, 6, 7,
1, 0, 4, 5,
4, 0, 3, 7,
6, 5, 4, 7 }};
constexpr std::array<int, 24> face_indices
{ { 0, 1, 2, 3,
2, 1, 5, 6,
3, 2, 6, 7,
1, 0, 4, 5,
4, 0, 3, 7,
6, 5, 4, 7 } };
constexpr std::array<int, 24> edge_indices
{ { 0, 1, 2, 3,
1, 4, 5, 6,
2, 6, 7, 8,
0, 9, 10, 4,
9, 3, 8, 11,
5, 10, 11, 7 } };
std::array<int, 12> edge_ipt_id;
edge_ipt_id.fill(-1);
auto inter_pt_index =
[&plane, &corners, &edge_ipt_id](int i, int j, int edge_id)
{
if (edge_ipt_id[edge_id]==-1)
{
edge_ipt_id[edge_id] = static_cast<int> (corners.size());
corners.push_back(typename Geom_traits::Construct_plane_line_intersection_point_3()
(plane, corners[i], corners[j]));
}
return edge_ipt_id[edge_id];
};
std::map<std::pair<int,int>, int> id_map;
std::vector< std::vector<int> > output_faces(6);
bool all_in = true;
bool all_out = true;
std::set<int> in_point_ids; // to collect the set of points in the clipped bbox
std::bitset<14> in_point_bits; // 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)
@@ -127,6 +133,7 @@ clip_to_bbox(const Plane_3& plane,
{
int current_id = face_indices[4*i + k];
int next_id = face_indices[4*i + (k+1)%4];
int edge_id = edge_indices[4 * i + k];
switch(orientations[ current_id ])
{
@@ -135,13 +142,13 @@ 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());
in_point_bits.set(output_faces[i].back());
// check for intersection of the edge
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(current_id, next_id, edge_id));
in_point_bits.set(output_faces[i].back());
}
break;
}
@@ -152,15 +159,15 @@ clip_to_bbox(const Plane_3& plane,
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(current_id, next_id, edge_id));
in_point_bits.set(output_faces[i].back());
}
break;
}
case ON_ORIENTED_BOUNDARY:
{
output_faces[i].push_back(current_id);
in_point_ids.insert(output_faces[i].back());
in_point_bits.set(output_faces[i].back());
}
}
}
@@ -182,11 +189,14 @@ clip_to_bbox(const Plane_3& plane,
typedef typename graph_traits::face_descriptor face_descriptor;
std::map<int, vertex_descriptor> out_vertices;
for(int i : in_point_ids)
for(int i=0; i<14;++i)
{
vertex_descriptor v = add_vertex(tm_out);
out_vertices.insert(std::make_pair(i, v));
put(vpm_out, v, corners[i]);
if (in_point_bits.test(i))
{
vertex_descriptor v = add_vertex(tm_out);
out_vertices.insert(std::make_pair(i, v));
put(vpm_out, v, corners[i]);
}
}
std::map< std::pair<int,int>, halfedge_descriptor> hedge_map;
@@ -1371,7 +1371,17 @@ public:
CGAL::Bounded_side in_tm2 = is_tm2_inside_out
? ON_UNBOUNDED_SIDE : ON_BOUNDED_SIDE;
Side_of_triangle_mesh<TriangleMesh, Kernel, VertexPointMap2> inside_tm2(tm2, vpm2);
typedef typename Nodes_vector::Exact_kernel Exact_kernel;
typedef Side_of_helper<TriangleMesh,
Node_id_map,
VertexPointMap2,
Nodes_vector, Kernel> VPM_helper;
typedef typename VPM_helper::VPM SOTM_vpm2;
typedef typename VPM_helper::Tree_type Tree_type;
Tree_type tree;
VPM_helper::build_tree(tm2, tree, vertex_to_node_id2, fids2, vpm2, nodes);
Side_of_triangle_mesh<TriangleMesh, Exact_kernel, SOTM_vpm2, Tree_type> inside_tm2(tree);
for(face_descriptor f : faces(tm1))
{
@@ -1382,16 +1392,20 @@ public:
patch_status_not_set_tm1.reset( patch_id );
halfedge_descriptor h = halfedge(f, tm1);
Node_id index_p1 = get_node_id(target(h, tm1), vertex_to_node_id1);
std::array<Node_id, 3> fnids = { index_p1, index_p1, index_p1 };
if (index_p1 != NID)
{
h=next(h, tm1);
index_p1 = get_node_id(target(h, tm1), vertex_to_node_id1);
fnids[1]=index_p1;
if (index_p1 != NID)
{
h=next(h, tm1);
index_p1 = get_node_id(target(h, tm1), vertex_to_node_id1);
fnids[2]=index_p1;
}
}
if (index_p1 != NID)
{
if (coplanar_patches_of_tm1.test(patch_id))
@@ -1401,12 +1415,13 @@ public:
}
else
{
// triangle which is tangent at its 3 vertices
// \todo improve this part which is not robust with a kernel
// with inexact constructions.
Bounded_side position = inside_tm2(centroid(get(vpm1, source(h, tm1)),
get(vpm1, target(h, tm1)),
get(vpm1, target(next(h, tm1), tm1)) ));
typename Exact_kernel::Point_3 e_centroid =
centroid(nodes.exact_node(fnids[0]),
nodes.exact_node(fnids[1]),
nodes.exact_node(fnids[2]));
Bounded_side position = inside_tm2(e_centroid);
CGAL_assertion( position != ON_BOUNDARY);
if ( position == in_tm2 )
is_patch_inside_tm2.set(patch_id);
@@ -1414,9 +1429,7 @@ public:
}
else
{
// TODO: tm2 might have been modified and an inexact vpm will
// provide a non-robust result.
Bounded_side position = inside_tm2( get(vpm1, target(h, tm1)));
Bounded_side position = inside_tm2( nodes.to_exact(get(vpm1, target(h, tm1))));
CGAL_assertion( position != ON_BOUNDARY);
if ( position == in_tm2 )
is_patch_inside_tm2.set(patch_id);
@@ -1433,7 +1446,18 @@ public:
CGAL::Bounded_side in_tm1 = is_tm1_inside_out
? ON_UNBOUNDED_SIDE : ON_BOUNDED_SIDE;
Side_of_triangle_mesh<TriangleMesh, Kernel, VertexPointMap1> inside_tm1(tm1, vpm1);
typedef typename Nodes_vector::Exact_kernel Exact_kernel;
typedef Side_of_helper<TriangleMesh,
Node_id_map,
VertexPointMap1,
Nodes_vector, Kernel> VPM_helper;
typedef typename VPM_helper::VPM SOTM_vpm1;
typedef typename VPM_helper::Tree_type Tree_type;
Tree_type tree;
VPM_helper::build_tree(tm1, tree, vertex_to_node_id1, fids1, vpm1, nodes);
Side_of_triangle_mesh<TriangleMesh, Exact_kernel, SOTM_vpm1, Tree_type> inside_tm1(tree);
for(face_descriptor f : faces(tm2))
{
const std::size_t f_id = get(fids2, f);
@@ -1443,14 +1467,17 @@ public:
patch_status_not_set_tm2.reset( patch_id );
halfedge_descriptor h = halfedge(f, tm2);
Node_id index_p2 = get_node_id(target(h, tm2), vertex_to_node_id2);
std::array<Node_id, 3> fnids = { index_p2, index_p2, index_p2 };
if (index_p2 != NID)
{
h=next(h, tm2);
index_p2 = get_node_id(target(h, tm2), vertex_to_node_id2);
fnids[1]=index_p2;
if (index_p2 != NID)
{
h=next(h, tm2);
index_p2 = get_node_id(target(h, tm2), vertex_to_node_id2);
fnids[2]=index_p2;
}
}
if (index_p2 != NID)
@@ -1462,11 +1489,11 @@ public:
}
else
{
// triangle which is tangent at its 3 vertices
// \todo improve this part which is not robust with a kernel
// with inexact constructions.
Bounded_side position = inside_tm1(midpoint(get(vpm2, source(h, tm2)),
get(vpm2, target(h, tm2)) ));
typename Exact_kernel::Point_3 e_centroid =
centroid(nodes.exact_node(fnids[0]),
nodes.exact_node(fnids[1]),
nodes.exact_node(fnids[2]));
Bounded_side position = inside_tm1(e_centroid);
CGAL_assertion( position != ON_BOUNDARY);
if ( position == in_tm1 )
is_patch_inside_tm1.set(patch_id);
@@ -1474,9 +1501,7 @@ public:
}
else
{
// TODO: tm1 might have been modified and an inexact vpm will
// provide a non-robust result.
Bounded_side position = inside_tm1( get(vpm2, target(h, tm2)));
Bounded_side position = inside_tm1( nodes.to_exact(get(vpm2, target(h, tm2))));
CGAL_assertion( position != ON_BOUNDARY);
if ( position == in_tm1 )
is_patch_inside_tm1.set(patch_id);
@@ -146,6 +146,207 @@ void copy_edge_mark(G&,
No_mark<G>&)
{} // nothing to do
// For exact side_of_triangle_mesh
template <class Node_id_map,
class VertexPointMap,
class NodeVector>
struct Node_vector_exact_vertex_point_map
{
// map type definitions
typedef typename boost::property_traits<VertexPointMap>::key_type key_type;
typedef typename NodeVector::Exact_kernel Exact_kernel;
typedef typename Exact_kernel::Point_3 value_type;
typedef value_type reference;
typedef boost::readable_property_map_tag category;
// internal type definitions
typedef std::size_t Node_id;
Node_vector_exact_vertex_point_map(){}
Node_vector_exact_vertex_point_map(const Node_id_map& node_ids,
const VertexPointMap& vpm,
const NodeVector& node_vector)
: node_ids(&node_ids)
, vpm(&vpm)
, node_vector(&node_vector)
{}
friend value_type get(Node_vector_exact_vertex_point_map m, key_type k)
{
typename Node_id_map::const_iterator it = m.node_ids->find(k);
if (it == m.node_ids->end())
return m.node_vector->to_exact( get(*(m.vpm), k) );
return m.node_vector->exact_node(it->second);
}
const Node_id_map* node_ids;
const VertexPointMap* vpm;
const NodeVector* node_vector;
};
// For exact side_of_triangle_mesh
template <class TriangleMesh, class PPM, class TreeTraits>
struct Split_primitives
{
Split_primitives(TriangleMesh& tm, PPM ppm)
: tm(tm)
, ppm(ppm)
{}
template<typename PrimitiveIterator>
void operator()(PrimitiveIterator first,
PrimitiveIterator beyond,
const CGAL::Bbox_3& bbox) const
{
typedef typename std::iterator_traits<PrimitiveIterator>::value_type Prmtv;
PrimitiveIterator middle = first + (beyond - first)/2;
typedef typename std::iterator_traits<PrimitiveIterator>::value_type Prmtv;
switch(TreeTraits::longest_axis(bbox))
{
case TreeTraits::CGAL_AXIS_X: // sort along x
std::nth_element(first, middle, beyond, [this](const Prmtv& p1, const Prmtv& p2){ return get(ppm, p1.id()).x() < get(ppm,p2.id()).x(); });
break;
case TreeTraits::CGAL_AXIS_Y: // sort along y
std::nth_element(first, middle, beyond, [this](const Prmtv& p1, const Prmtv& p2){ return get(ppm, p1.id()).y() < get(ppm,p2.id()).y(); });
break;
case TreeTraits::CGAL_AXIS_Z: // sort along z
std::nth_element(first, middle, beyond, [this](const Prmtv& p1, const Prmtv& p2){ return get(ppm, p1.id()).z() < get(ppm,p2.id()).z(); });
break;
default:
CGAL_error();
}
}
TriangleMesh& tm;
PPM ppm;
};
// For exact side_of_triangle_mesh
template <class BPM>
struct Compute_bbox {
Compute_bbox(const BPM& bpm)
: bpm(bpm)
{}
template<typename ConstPrimitiveIterator>
CGAL::Bbox_3 operator()(ConstPrimitiveIterator first,
ConstPrimitiveIterator beyond) const
{
CGAL::Bbox_3 bbox = get(bpm, first->id());
for(++first; first != beyond; ++first)
{
bbox += get(bpm, first->id());
}
return bbox;
}
BPM bpm;
};
// For exact side_of_triangle_mesh
template <class TriangleMesh,
class Node_id_map,
class VertexPointMap,
class NodeVector,
class Input_Kernel>
struct Side_of_helper
{
typedef Node_vector_exact_vertex_point_map<Node_id_map, VertexPointMap, NodeVector> VPM;
typedef CGAL::AABB_face_graph_triangle_primitive<TriangleMesh, VPM> Primitive;
typedef CGAL::AABB_traits<typename NodeVector::Exact_kernel, Primitive> Traits;
typedef CGAL::AABB_tree<Traits> Tree_type;
static
VPM get_vpm(const Node_id_map& node_ids,
const VertexPointMap& vpm,
const NodeVector& node_vector)
{
return VPM(node_ids, vpm, node_vector);
}
template <class FaceIdMap>
static
void build_tree(TriangleMesh& tm,
Tree_type& tree,
const Node_id_map& node_ids,
FaceIdMap fid,
const VertexPointMap& vpm,
const NodeVector& node_vector)
{
typedef typename boost::graph_traits<TriangleMesh>::face_descriptor face_descriptor;
typedef typename boost::graph_traits<TriangleMesh>::halfedge_descriptor halfedge_descriptor;
typedef typename boost::graph_traits<TriangleMesh>::vertex_descriptor vertex_descriptor;
// add primitives
tree.insert(faces(tm).begin(), faces(tm).end(), tm, get_vpm(node_ids, vpm, node_vector));
// pre-build bboxes (using approximation)
std::vector<Bbox_3> face_bboxes(num_faces(tm));
auto get_v_box = [&node_ids, &node_vector, &vpm](vertex_descriptor v)
{
typename Node_id_map::const_iterator it = node_ids.find(v);
if (it == node_ids.end())
return get(vpm, v).bbox();
return approx(node_vector.exact_node(it->second)).bbox();
};
for (face_descriptor f : faces(tm))
{
halfedge_descriptor h = halfedge(f, tm);
face_bboxes[get(fid, f)] = get_v_box( source(h, tm) ) +
get_v_box( target(h, tm) ) +
get_v_box( target(next(h, tm), tm) );
}
typedef CGAL::Pointer_property_map<CGAL::Bbox_3>::type Id_to_box;
Id_to_box id_to_box = CGAL::make_property_map(face_bboxes);
typedef Property_map_binder<FaceIdMap, Id_to_box> BPM;
BPM bpm(fid, id_to_box);
Compute_bbox<BPM> compute_bbox(bpm);
typedef One_point_from_face_descriptor_map<TriangleMesh, VertexPointMap> PPM;
PPM ppm(&tm, vpm);
Split_primitives<TriangleMesh, PPM, Traits> split_primitives(tm, ppm);
tree.custom_build(compute_bbox, split_primitives);
}
};
template <class TriangleMesh,
class Node_id_map,
class VertexPointMap,
class NodeVector>
struct Side_of_helper<TriangleMesh, Node_id_map, VertexPointMap, NodeVector, typename NodeVector::Exact_kernel>
{
typedef VertexPointMap VPM;
static
VPM get_vpm(const Node_id_map&,
const VertexPointMap& vpm,
const NodeVector&)
{
return vpm;
}
typedef CGAL::AABB_face_graph_triangle_primitive<TriangleMesh, VPM> Primitive;
typedef CGAL::AABB_traits<typename NodeVector::Exact_kernel, Primitive> Traits;
typedef CGAL::AABB_tree<Traits> Tree_type;
template <class FaceIdMap>
static
void build_tree(TriangleMesh& tm,
Tree_type& tree,
const Node_id_map& /* node_ids */,
FaceIdMap /* fid */,
const VertexPointMap& vpm,
const NodeVector& /* node_vector */)
{
tree.insert(faces(tm).begin(), faces(tm).end(), tm, vpm);
tree.build();
}
};
// Parts to get default property maps for output meshes based on the value type
// of input vertex point maps.
template <typename Point_3, typename vertex_descriptor>
@@ -424,12 +625,9 @@ struct Patch_description{
// shared_edges will be filled by halfedges pointing in the patch
// that are inside `is_intersection_edge`, thus mesh boundary halfedges
// are not necessarily inside.
template <class PolygonMesh, class FaceIndexMap, class IsIntersectionEdge>
template <class PolygonMesh, class IsIntersectionEdge>
void extract_patch_simplices(
std::size_t patch_id,
PolygonMesh& pm,
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,
std::vector<typename boost::graph_traits<PolygonMesh>::halfedge_descriptor>& interior_edges,
@@ -441,22 +639,18 @@ void extract_patch_simplices(
typedef typename GT::vertex_descriptor vertex_descriptor;
typedef typename GT::face_descriptor face_descriptor;
for(face_descriptor f : faces(pm))
for(face_descriptor f : patch_faces)
{
if ( patch_ids[ get(fids, f) ]==patch_id )
for(halfedge_descriptor h :
halfedges_around_face(halfedge(f, pm),pm))
{
patch_faces.push_back( f );
for(halfedge_descriptor h :
halfedges_around_face(halfedge(f, pm),pm))
if ( !is_intersection_edge.count(edge(h, pm)) )
{
if ( !is_intersection_edge.count(edge(h, pm)) )
{
if ( h < opposite(h,pm) || is_border(opposite(h,pm),pm) )
interior_edges.push_back( h );
}
else
shared_edges.push_back(h);
if ( h < opposite(h,pm) || is_border(opposite(h,pm),pm) )
interior_edges.push_back( h );
}
else
shared_edges.push_back(h);
}
}
@@ -498,14 +692,19 @@ struct Patch_container{
, patch_ids(patch_ids)
, fids(fids)
, is_intersection_edge(is_intersection_edge)
{}
{
typedef boost::graph_traits<PolygonMesh> GT;
typedef typename GT::face_descriptor face_descriptor;
for(face_descriptor f : faces(pm))
patches[patch_ids[ get(fids, f) ]].faces.push_back( f );
}
Patch_description<PolygonMesh>& operator[](std::size_t i) {
if ( !patches[i].is_initialized )
{
extract_patch_simplices(
i, pm,
fids, patch_ids,
pm,
patches[i].faces, patches[i].interior_vertices,
patches[i].interior_edges, patches[i].shared_edges,
is_intersection_edge
@@ -411,7 +411,7 @@ class Intersection_of_triangle_meshes
}
Key key(ipt.type_1, ipt.type_2, h1, h2);
if (&tm1==&tm2 && h1>h2)
if (&tm1==&tm2 && h2<h1)
key=Key(ipt.type_2, ipt.type_1, h2, h1);
std::pair<typename std::map<Key,Node_id>::iterator,bool> res=
@@ -28,10 +28,13 @@
#include <CGAL/Lazy.h> // needed for CGAL::exact(FT)/CGAL::exact(Lazy_exact_nt<T>)
#include <boost/container/small_vector.hpp>
#include <boost/unordered_set.hpp>
#include <boost/graph/graph_traits.hpp>
#include <boost/dynamic_bitset.hpp>
#include <utility>
#include <algorithm>
#ifdef DOXYGEN_RUNNING
#define CGAL_PMP_NP_TEMPLATE_PARAMETERS NamedParameters
@@ -50,6 +53,14 @@ public:
namespace Polygon_mesh_processing {
namespace pmp_internal {
inline void rearrange_face_ids(boost::container::small_vector<std::size_t, 4>& ids)
{
auto min_elem = std::min_element(ids.begin(), ids.end());
std::rotate(ids.begin(), min_elem, ids.end());
}
}//end pmp_internal
/**
* \ingroup measure_grp
* computes the length of an edge of a given polygon mesh.
@@ -820,6 +831,176 @@ centroid(const TriangleMesh& tmesh)
return centroid(tmesh, CGAL::Polygon_mesh_processing::parameters::all_default());
}
/**
* \ingroup measure_grp
* identifies faces only present in `m1` and `m2` as well as the faces present
* in both polygon meshes. Two faces are identical if they have the same
* orientation and the same points.
*
* @tparam PolygonMesh a model of `HalfedgeListGraph` and `FaceListGraph`
* @tparam FaceOutputIterator model of `OutputIterator`
holding `boost::graph_traits<PolygonMesh>::%face_descriptor`.
* @tparam FacePairOutputIterator model of `OutputIterator`
holding `std::pair<boost::graph_traits<PolygonMesh>::%face_descriptor,
boost::graph_traits<PolygonMesh>::%face_descriptor`.
*
* @tparam NamedParameters1 a sequence of \ref bgl_namedparameters "Named Parameters"
* @tparam NamedParameters2 a sequence of \ref bgl_namedparameters "Named Parameters"
*
* @param m1 the first `PolygonMesh`
* @param m2 the second `PolygonMesh`
* @param common output iterator collecting the faces that are common to both meshes.
* @param m1_only output iterator collecting the faces that are only in `m1`
* @param m2_only output iterator collecting the faces that are only in `m2`
* @param np1 an optional sequence of \ref bgl_namedparameters "Named Parameters" among the ones listed below
* @param np2 an optional sequence of \ref bgl_namedparameters "Named Parameters" among the ones listed below
*
* \cgalNamedParamsBegin
* \cgalParamNBegin{vertex_point_map}
* \cgalParamDescription{a property map associating points to the vertices of `m1` (`m2`)}
* \cgalParamType{a class model of `ReadablePropertyMap` with `boost::graph_traits<PolygonMesh>::%vertex_descriptor`
* as key type and `%Point_3` as value type. `%Point_3` must be LessThanComparable.}
* \cgalParamDefault{`boost::get(CGAL::vertex_point, m1 (m2))`}
* \cgalParamNEnd
*
* \cgalParamNBegin{vertex_index_map}
* \cgalParamDescription{a property map associating to each vertex of `m1` (`m2`) a unique index between `0` and `num_vertices(m1 (m2)) - 1`}
* \cgalParamType{a class model of `ReadablePropertyMap` with `boost::graph_traits<Graph>::%vertex_descriptor`
* as key type and `std::size_t` as value type}
* \cgalParamDefault{an automatically indexed internal map}
* \cgalParamExtra{If this parameter is not passed, internal machinery will create and initialize
* a face index property map, either using the internal property map if it exists
* or using an external map. The latter might result in - slightly - worsened performance
* in case of non-constant complexity for index access.}
* \cgalParamNEnd
* \cgalNamedParamsEnd
*
*/
template<typename PolygonMesh, typename FaceOutputIterator, typename FacePairOutputIterator, typename NamedParameters1, typename NamedParameters2 >
void match_faces(const PolygonMesh& m1, const PolygonMesh& m2,
FacePairOutputIterator common, FaceOutputIterator m1_only, FaceOutputIterator m2_only,
const NamedParameters1& np1, const NamedParameters2& np2)
{
using parameters::choose_parameter;
using parameters::get_parameter;
typedef typename GetVertexPointMap < PolygonMesh, NamedParameters1>::const_type VPMap1;
typedef typename GetVertexPointMap < PolygonMesh, NamedParameters2>::const_type VPMap2;
typedef typename GetInitializedVertexIndexMap<PolygonMesh, NamedParameters1>::const_type VIMap1;
typedef typename GetInitializedVertexIndexMap<PolygonMesh, NamedParameters2>::const_type VIMap2;
VPMap1 vpm1 = choose_parameter(get_parameter(np1, internal_np::vertex_point),
get_const_property_map(vertex_point, m1));
VPMap2 vpm2 = choose_parameter(get_parameter(np2, internal_np::vertex_point),
get_const_property_map(vertex_point, m2));
VIMap1 vim1 = get_initialized_vertex_index_map(m1, np1);
VIMap2 vim2 = get_initialized_vertex_index_map(m2, np2);
typedef typename boost::property_traits<VPMap2>::value_type Point_3;
typedef typename boost::graph_traits<PolygonMesh>::face_descriptor face_descriptor;
std::map<Point_3, std::size_t> point_id_map;
std::vector<std::size_t> m1_vertex_id(vertices(m1).size(), -1);
std::vector<std::size_t> m2_vertex_id(vertices(m2).size(), -1);
boost::dynamic_bitset<> shared_vertices(m1_vertex_id.size() + m2_vertex_id.size());
//iterate both meshes to set ids of all points, and set vertex/point_id maps.
std::size_t id = 0;
for(auto v : vertices(m1))
{
const Point_3& p = get(vpm1, v);
auto res = point_id_map.insert(std::make_pair(p, id));
if(res.second)
++id;
m1_vertex_id[get(vim1, v)]=res.first->second;
}
for(auto v : vertices(m2))
{
const Point_3& p = get(vpm2, v);
auto res = point_id_map.insert(std::make_pair(p, id));
if(res.second)
++id;
else
shared_vertices.set(res.first->second);
m2_vertex_id[get(vim2, v)]=res.first->second;
}
//fill a set with the "faces point-ids" of m1 and then iterate faces of m2 to compare.
std::map<boost::container::small_vector<std::size_t, 4>, face_descriptor> m1_faces_map;
for(auto f : faces(m1))
{
bool all_shared = true;
boost::container::small_vector<std::size_t, 4> ids;
for(auto v : CGAL::vertices_around_face(halfedge(f, m1), m1))
{
std::size_t vid = m1_vertex_id[get(vim1, v)];
ids.push_back(vid);
if(!shared_vertices.test(vid))
{
all_shared = false;
break;
}
}
if(all_shared)
{
pmp_internal::rearrange_face_ids(ids);
m1_faces_map.insert({ids, f});
}
else
*m1_only++ = f;
}
for(auto f : faces(m2))
{
boost::container::small_vector<std::size_t, 4> ids;
bool all_shared = true;
for(auto v : CGAL::vertices_around_face(halfedge(f, m2), m2))
{
std::size_t vid = m2_vertex_id[get(vim2, v)];
ids.push_back(vid);
if(!shared_vertices.test(vid))
{
all_shared = false;
break;
}
}
if(all_shared)
{
pmp_internal::rearrange_face_ids(ids);
auto it = m1_faces_map.find(ids);
if(it != m1_faces_map.end())
{
*common++ = std::make_pair(it->second, f);
m1_faces_map.erase(it);
}
else
{
*m2_only++ = f;
}
}
else
*m2_only++ = f;
}
//all shared faces have been removed from the map, so all that remains must go in m1_only
for(const auto& it : m1_faces_map)
{
*m1_only++ = it.second;
}
}
template<typename PolygonMesh, typename FaceOutputIterator, typename FacePairOutputIterator, typename NamedParameters>
void match_faces(const PolygonMesh& m1, const PolygonMesh& m2,
FacePairOutputIterator common, FaceOutputIterator m1_only, FaceOutputIterator m2_only,
const NamedParameters& np)
{
match_faces(m1, m2, common, m1_only, m2_only, np, parameters::all_default());
}
template<typename PolygonMesh, typename FaceOutputIterator, typename FacePairOutputIterator>
void match_faces(const PolygonMesh& m1, const PolygonMesh& m2,
FacePairOutputIterator common, FaceOutputIterator m1_only, FaceOutputIterator m2_only)
{
match_faces(m1, m2, common, m1_only, m2_only, parameters::all_default(), parameters::all_default());
}
} // namespace Polygon_mesh_processing
} // namespace CGAL
@@ -53,8 +53,8 @@ namespace internal{
{}
typedef bool result_type;
template <class vertex_descriptor>
bool operator()(vertex_descriptor v1, vertex_descriptor v2) const
template <class vertex_descriptor1, class vertex_descriptor2>
bool operator()(vertex_descriptor1 v1, vertex_descriptor2 v2) const
{
return CGAL::SMALLER == compare_z(get(vpmap, v1), get(vpmap, v2));
}
@@ -50,11 +50,12 @@ namespace CGAL {
namespace Polygon_mesh_processing {
namespace internal {
template <typename TriangleMesh, typename VPM, typename ECM, typename Traits>
template <typename TriangleMesh, typename VPM, typename VCM, 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 VCM& vcm,
const ECM& ecm,
const Traits& gt,
const double cap_threshold, // angle over 160° ==> cap
@@ -70,58 +71,61 @@ is_badly_shaped(const typename boost::graph_traits<TriangleMesh>::face_descripto
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)))
if(res != null_h && (!get(vcm, source(res, tmesh)) || !get(vcm, target(res, tmesh))) )
{
// don't want to collapse edges that are too large
if(collapse_length_threshold == 0 ||
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);
}
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>
template <typename TriangleMesh, typename HalfedgeContainer,
typename VPM, typename VCM, typename ECM, typename Traits>
void collect_badly_shaped_triangles(const typename boost::graph_traits<TriangleMesh>::face_descriptor f,
TriangleMesh& tmesh,
const VPM& vpm,
const VCM& vcm,
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)
HalfedgeContainer& edges_to_collapse,
HalfedgeContainer& 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,
std::array<halfedge_descriptor, 2> res = is_badly_shaped(f, tmesh, vpm, vcm, ecm, gt, cap_threshold,
needle_threshold, collapse_length_threshold);
if(res[0] != boost::graph_traits<TriangleMesh>::null_halfedge())
{
#ifdef CGAL_PMP_DEBUG_REMOVE_DEGENERACIES
#ifdef CGAL_PMP_DEBUG_REMOVE_DEGENERACIES_EXTRA
std::cout << "add new needle: " << edge(res[0], tmesh) << std::endl;
#endif
edges_to_collapse.insert(edge(res[0], tmesh));
CGAL_assertion(!is_border(res[0], tmesh));
CGAL_assertion(!get(ecm, edge(res[0], tmesh)));
edges_to_collapse.insert(res[0]);
}
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
#ifdef CGAL_PMP_DEBUG_REMOVE_DEGENERACIES_EXTRA
std::cout << "add new cap: " << edge(res[1],tmesh) << std::endl;
#endif
edges_to_flip.insert(edge(res[1], tmesh));
CGAL_assertion(!is_border(res[1], tmesh));
CGAL_assertion(!get(ecm, edge(res[1], tmesh)));
edges_to_flip.insert(res[1]);
}
}
}
@@ -254,8 +258,6 @@ get_best_edge_orientation(typename boost::graph_traits<TriangleMesh>::edge_descr
halfedge_descriptor h = halfedge(e, tmesh), ho = opposite(h, tmesh);
CGAL_assertion(!get(vcm, source(h, tmesh)) || !get(vcm, target(h, tmesh)));
boost::optional<FT> dv1 = get_collapse_volume(h, tmesh, vpm, gt);
boost::optional<FT> dv2 = get_collapse_volume(ho, tmesh, vpm, gt);
@@ -381,6 +383,9 @@ bool remove_almost_degenerate_faces(const FaceRange& face_range,
typedef typename boost::property_map<TriangleMesh, Vertex_property_tag>::type DVCM;
DVCM vcm = get(Vertex_property_tag(), tmesh);
CGAL_precondition(is_valid_polygon_mesh(tmesh));
CGAL_precondition(is_triangle_mesh(tmesh));
for(face_descriptor f : face_range)
{
if(f == boost::graph_traits<TriangleMesh>::null_face())
@@ -401,15 +406,22 @@ bool remove_almost_degenerate_faces(const FaceRange& face_range,
}
// Start the process of removing bad elements
std::set<edge_descriptor> edges_to_collapse;
std::set<edge_descriptor> edges_to_flip;
std::set<halfedge_descriptor> edges_to_collapse;
std::set<halfedge_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,
{
internal::collect_badly_shaped_triangles(f, tmesh, vpm, vcm, ecm, gt,
cap_threshold, needle_threshold, collapse_length_threshold,
edges_to_collapse, edges_to_flip);
}
#ifdef CGAL_PMP_DEBUG_REMOVE_DEGENERACIES
std::cout << edges_to_collapse.size() << " to collapse" << std::endl;
std::cout << edges_to_flip.size() << " to flip" << std::endl;
#endif
#ifdef CGAL_PMP_DEBUG_REMOVE_DEGENERACIES
int iter = 0;
@@ -421,6 +433,7 @@ bool remove_almost_degenerate_faces(const FaceRange& face_range,
#ifdef CGAL_PMP_DEBUG_REMOVE_DEGENERACIES
std::cout << edges_to_collapse.size() << " needles and " << edges_to_flip.size() << " caps" << std::endl;
std::cout << "Iter: " << iter << std::endl;
std::ostringstream oss;
oss << "degen_cleaning_iter_" << iter++ << ".off";
CGAL::write_polygon_mesh(oss.str(), tmesh, CGAL::parameters::stream_precision(17));
@@ -430,102 +443,114 @@ bool remove_almost_degenerate_faces(const FaceRange& face_range,
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;
std::set<halfedge_descriptor> next_edges_to_collapse;
std::set<halfedge_descriptor> next_edges_to_flip;
// treat needles
// 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();
halfedge_descriptor h = *edges_to_collapse.begin();
edges_to_collapse.erase(edges_to_collapse.begin());
CGAL_assertion(!get(ecm, e));
CGAL_assertion(!is_border(h, tmesh));
if(get(vcm, source(e, tmesh)) && get(vcm, target(e, tmesh)))
const edge_descriptor e = edge(h, tmesh);
CGAL_assertion(!get(ecm, edge(h, tmesh)));
if(get(vcm, source(h, tmesh)) && get(vcm, target(h, 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;
#ifdef CGAL_PMP_DEBUG_REMOVE_DEGENERACIES_EXTRA
std::cout << " treat needle: " << e
<< " (" << source(e, tmesh) << " " << tmesh.point(source(h, tmesh))
<< " --- " << source(e, tmesh) << " " << tmesh.point(target(h, 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,
// Verify that the element is still badly shaped
const std::array<halfedge_descriptor, 2> nc =
internal::is_badly_shaped(face(h, tmesh), tmesh, vpm, vcm, 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;
#ifdef CGAL_PMP_DEBUG_REMOVE_DEGENERACIES_EXTRA
std::cout << "\t Needle criteria no longer verified" << 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;
continue;
}
// pick the orientation of edge to keep the vertex minimizing the volume variation
const halfedge_descriptor best_h = internal::get_best_edge_orientation(e, tmesh, vpm, vcm, gt);
if(best_h == boost::graph_traits<TriangleMesh>::null_halfedge())
{
#ifdef CGAL_PMP_DEBUG_REMOVE_DEGENERACIES_EXTRA
std::cout << "\t Geometrically invalid edge collapse!" << std::endl;
#endif
next_edges_to_collapse.insert(h);
continue;
}
// Proceeding with the collapse, purge the sets from halfedges being removed
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);
edges_to_flip.erase(h);
edges_to_collapse.erase(h);
next_edges_to_collapse.erase(h);
halfedge_descriptor rm_h = prev(h, tmesh);
if(get(ecm, edge(rm_h, tmesh)))
rm_h = next(h, tmesh);
edges_to_flip.erase(rm_h);
edges_to_collapse.erase(rm_h);
next_edges_to_collapse.erase(rm_h);
halfedge_descriptor opp_rm_h = opposite(rm_h, tmesh);
edges_to_flip.erase(opp_rm_h);
edges_to_collapse.erase(opp_rm_h);
next_edges_to_collapse.erase(opp_rm_h);
}
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;
std::cout << " " << kk << " -- Collapsing " << tmesh.point(source(best_h, tmesh)) << " "
<< tmesh.point(target(best_h, tmesh)) << std::endl;
#endif
CGAL_assertion(!get(vcm, source(best_h, tmesh)));
// The function get_best_edge_orientation() has ensured that get(vcm, source(h, tmesh))
// is not constrained, so get(ecm, e) is also not constrained for all e incident
// to source(h, tmesh).
//
// The function Euler::collapse_edge() removes edge(prev(h, tmesh)), but that edge
// might be constrained. In that case, next() must be removed instead.
vertex_descriptor v;
if(get(ecm, edge(prev(h, tmesh), tmesh)))
v = Euler::collapse_edge(edge(best_h, tmesh), tmesh, ecm);
else
v = Euler::collapse_edge(edge(best_h, tmesh), tmesh);
// 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)));
// tmesh.point(v) = mp;
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,
internal::collect_badly_shaped_triangles(face(hv, tmesh), tmesh, vpm, vcm, ecm, gt,
cap_threshold, needle_threshold, collapse_length_threshold,
edges_to_collapse, edges_to_flip);
}
@@ -541,92 +566,84 @@ bool remove_almost_degenerate_faces(const FaceRange& face_range,
#endif
something_was_done = true;
}
else
else // ! CGAL::Euler::does_satisfy_link_condition(e, tmesh)
{
#ifdef CGAL_PMP_DEBUG_REMOVE_DEGENERACIES
std::cerr << "Warning: uncollapsable edge! " << tmesh.point(source(e, tmesh)) << " "
<< tmesh.point(target(e, tmesh)) << std::endl;
#ifdef CGAL_PMP_DEBUG_REMOVE_DEGENERACIES_EXTRA
std::cout << "\t Uncollapsable edge!" << std::endl;
#endif
next_edges_to_collapse.insert(e);
next_edges_to_collapse.insert(h);
}
}
// treat caps
// Treat caps ==================================================================================
CGAL_assertion(next_edges_to_flip.empty());
#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();
halfedge_descriptor h = *edges_to_flip.begin();
edges_to_flip.erase(edges_to_flip.begin());
CGAL_assertion(!is_border(h, tmesh));
const edge_descriptor e = edge(h, tmesh);
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;
#ifdef CGAL_PMP_DEBUG_REMOVE_DEGENERACIES_EXTRA
std::cout << " treat cap: " << e
<< " (" << source(e, tmesh) << " " << tmesh.point(source(h, tmesh))
<< " --- " << target(e, tmesh) << " " << tmesh.point(target(h, 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,
std::array<halfedge_descriptor,2> nc = internal::is_badly_shaped(face(h, tmesh), tmesh, vpm, vcm, ecm, gt,
cap_threshold, needle_threshold,
collapse_length_threshold);
// First check the triangle is still a cap
// 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;
#ifdef CGAL_PMP_DEBUG_REMOVE_DEGENERACIES_EXTRA
std::cout << "\t Cap criteria no longer verified" << 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;
continue;
}
// special case on the border
// special case of `edge(h, tmesh)` being a border edge --> remove the face
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));
for(halfedge_descriptor hh : CGAL::halfedges_around_face(h, tmesh))
{
// Remove from even 'next_edges_to_flip' because it might have been re-added from a flip
edges_to_flip.erase(hh);
next_edges_to_flip.erase(hh);
next_edges_to_collapse.erase(hh);
}
Euler::remove_face(h, tmesh);
something_was_done = true;
continue;
}
CGAL_assertion(!is_border(e, tmesh));
// 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;
#ifdef CGAL_PMP_DEBUG_REMOVE_DEGENERACIES_EXTRA
std::cout << "\t Flipping prevented: not the best diagonal" << std::endl;
#endif
next_edges_to_flip.insert(e);
next_edges_to_flip.insert(h);
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;
std::cout << "\t step " << kk << " -- Flipping" << std::endl;
#endif
Euler::flip_edge(h, tmesh);
CGAL_assertion(edge(h, tmesh) == e);
@@ -636,33 +653,28 @@ bool remove_almost_degenerate_faces(const FaceRange& face_range,
{
CGAL_assertion(!is_border(h, tmesh));
std::array<halfedge_descriptor, 2> nc =
internal::is_badly_shaped(face(h, tmesh), tmesh, vpm, ecm, gt,
internal::is_badly_shaped(face(h, tmesh), tmesh, vpm, vcm, 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));
}
}
if(nc[1] != boost::graph_traits<TriangleMesh>::null_halfedge() && nc[1] != h)
next_edges_to_flip.insert(nc[1]);
else if(nc[0] != boost::graph_traits<TriangleMesh>::null_halfedge())
next_edges_to_collapse.insert(nc[0]);
h = opposite(h, tmesh);
}
something_was_done = true;
}
#ifdef CGAL_PMP_DEBUG_REMOVE_DEGENERACIES
else
else // flipped edge already exists in the mesh
{
std::cerr << "Warning: unflippable edge! " << tmesh.point(source(h, tmesh)) << " --- "
<< tmesh.point(target(h, tmesh)) << std::endl;
next_edges_to_flip.insert(e);
}
#ifdef CGAL_PMP_DEBUG_REMOVE_DEGENERACIES_EXTRA
std::cout << "\t Unflippable edge!" << std::endl;
#endif
CGAL_assertion(!is_border(h, tmesh));
next_edges_to_flip.insert(h);
}
#ifdef CGAL_PMP_DEBUG_REMOVE_DEGENERACIES_EXTRA
std::string nb = std::to_string(++kk);
if(kk<10) nb = std::string("0")+nb;
@@ -673,11 +685,11 @@ bool remove_almost_degenerate_faces(const FaceRange& face_range,
#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;
std::swap(edges_to_collapse, next_edges_to_collapse);
std::swap(edges_to_flip, next_edges_to_flip);
}
return false;
@@ -1495,7 +1507,7 @@ bool remove_degenerate_edges(const EdgeRange& edge_range,
face_set.insert(face(hd, tmesh));
}
CGAL_assertion(is_valid_polygon_mesh(tmesh));
CGAL_expensive_assertion(is_valid_polygon_mesh(tmesh));
}
}
}
@@ -1582,6 +1594,7 @@ bool remove_degenerate_faces(const FaceRange& face_range,
const NamedParameters& np)
{
CGAL_assertion(CGAL::is_triangle_mesh(tmesh));
CGAL_assertion(CGAL::is_valid_polygon_mesh(tmesh));
using parameters::get_parameter;
using parameters::choose_parameter;
@@ -47,6 +47,7 @@
#include <exception>
#include <sstream>
#include <type_traits>
#include <typeinfo>
#include <vector>
#ifdef DOXYGEN_RUNNING
@@ -503,11 +504,21 @@ bool does_self_intersect(const FaceRange& face_range,
CGAL::Emptyset_iterator unused_out;
internal::self_intersections_impl<ConcurrencyTag>(face_range, tmesh, unused_out, true /*throw*/, np);
}
catch(CGAL::internal::Throw_at_output_exception&)
catch (const CGAL::internal::Throw_at_output_exception&)
{
return true;
}
#if defined(CGAL_LINKED_WITH_TBB) && TBB_USE_CAPTURED_EXCEPTION
catch (const tbb::captured_exception& e)
{
const char* ti1 = e.name();
const char* ti2 = typeid(const CGAL::internal::Throw_at_output_exception&).name();
const std::string tn1(ti1);
const std::string tn2(ti2);
if (tn1 == tn2) return true;
else throw;
}
#endif
return false;
}
@@ -358,6 +358,8 @@ is_needle_triangle_face(typename boost::graph_traits<TriangleMesh>::face_descrip
const NamedParameters& np)
{
CGAL_precondition(threshold >= 1.);
CGAL_precondition(f != boost::graph_traits<TriangleMesh>::null_face());
CGAL_precondition(CGAL::is_triangle(halfedge(f, tm), tm));
using parameters::get_parameter;
using parameters::choose_parameter;
@@ -462,7 +464,8 @@ is_cap_triangle_face(typename boost::graph_traits<TriangleMesh>::face_descriptor
const double threshold,
const NamedParameters& np)
{
CGAL_precondition(CGAL::is_triangle_mesh(tm));
CGAL_precondition(f != boost::graph_traits<TriangleMesh>::null_face());
CGAL_precondition(CGAL::is_triangle(halfedge(f, tm), tm));
CGAL_precondition(threshold >= -1.);
CGAL_precondition(threshold <= 0.);
@@ -150,6 +150,7 @@ struct Dummy_cycle_rep_maintainer
}
// Dummies just to fit the API
void add_representative(const halfedge_descriptor) const { }
void remove_representative(const halfedge_descriptor) const { }
void clear_representatives() const { }
@@ -214,13 +215,13 @@ public:
{
typedef typename boost::property_traits<VPM>::reference Point_ref;
CGAL_assertion(!cycle_halfedges.empty());
#ifdef CGAL_PMP_STITCHING_DEBUG
std::cout << "update_representatives(" << cycle_halfedges.size() << ", "
<< filtered_stitchable_halfedges.size() << ")" << std::endl;
#endif
CGAL_assertion(!cycle_halfedges.empty());
for(const halfedge_descriptor h : cycle_halfedges)
put(m_candidate_halfedges, h, true);
@@ -332,17 +333,22 @@ void fill_pairs(const Halfedge& he,
bool insertion_ok;
std::tie(set_it, insertion_ok) = border_halfedge_map.emplace(he, std::make_pair(1,0));
if(!insertion_ok) // we found already a halfedge with the points
if(!insertion_ok) // there is already a halfedge with the points
{
++set_it->second.first; // increase the multiplicity
if(set_it->second.first == 2)
{
const Halfedge other_he = set_it->first;
set_it->second.second = halfedge_pairs.size(); // set the id of the pair in the vector
halfedge_pairs.emplace_back(set_it->first, he);
if(get(vpm, source(he,pmesh)) == get(vpm, target(set_it->first, pmesh)) &&
get(vpm, target(he,pmesh)) == get(vpm, source(set_it->first, pmesh)))
halfedge_pairs.emplace_back(other_he, he);
if(get(vpm, source(he,pmesh)) == get(vpm, target(other_he, pmesh)) &&
get(vpm, target(he,pmesh)) == get(vpm, source(other_he, pmesh)))
{
manifold_halfedge_pairs.push_back(true);
// Even if the halfedges are compatible, refuse to stitch if that would break the graph
if(face(opposite(he, pmesh), pmesh) == face(opposite(other_he, pmesh), pmesh))
manifold_halfedge_pairs.push_back(false);
else
manifold_halfedge_pairs.push_back(true);
}
else
{
@@ -783,12 +789,12 @@ filter_stitchable_pairs(PolygonMesh& pmesh,
}
template <typename HalfedgePair, typename CandidateHalfedgeRange, typename PolygonMesh,
typename MaintainerVisitor, typename VertexPointMap>
typename CycleRepMaintainer, typename VertexPointMap>
std::size_t stitch_halfedge_range(const std::vector<HalfedgePair>& to_stitch,
const CandidateHalfedgeRange& representative_candidates,
PolygonMesh& pmesh,
MaintainerVisitor& mv,
const VertexPointMap& vpm)
const VertexPointMap& vpm,
CycleRepMaintainer& cycle_reps_maintainer)
{
typedef typename boost::graph_traits<PolygonMesh>::vertex_descriptor vertex_descriptor;
@@ -818,7 +824,7 @@ std::size_t stitch_halfedge_range(const std::vector<HalfedgePair>& to_stitch,
const std::vector<HalfedgePair>& to_stitch_filtered =
filter_stitchable_pairs(pmesh, to_stitch, to_stitch_local, uf_vertices, uf_handles);
mv.update_representatives(representative_candidates, to_stitch_filtered, vpm);
cycle_reps_maintainer.update_representatives(representative_candidates, to_stitch_filtered, vpm);
// Actually stitching
run_stitch_borders(pmesh, to_stitch_filtered, vpm, uf_vertices, uf_handles);
@@ -831,11 +837,11 @@ std::size_t stitch_halfedge_range(const std::vector<HalfedgePair>& to_stitch,
PolygonMesh& pmesh,
const VertexPointMap& vpm)
{
Dummy_cycle_rep_maintainer<PolygonMesh> mv(pmesh);
return stitch_halfedge_range(to_stitch, halfedges(pmesh), pmesh, mv, vpm);
Dummy_cycle_rep_maintainer<PolygonMesh> cycle_reps_maintainer(pmesh);
return stitch_halfedge_range(to_stitch, halfedges(pmesh), pmesh, vpm, cycle_reps_maintainer);
}
//overload to avoid a useless copy
// overload to avoid a useless copy
template <typename HalfedgePair, typename PolygonMesh, typename VertexPointMap>
std::size_t stitch_halfedge_range_dispatcher(const std::vector<HalfedgePair>& to_stitch,
PolygonMesh& pmesh,
@@ -844,7 +850,7 @@ std::size_t stitch_halfedge_range_dispatcher(const std::vector<HalfedgePair>& to
return stitch_halfedge_range(to_stitch, pmesh, vpm);
}
//overload to doing the copy
// overload making a copy
template <typename HalfedgePairRange, typename PolygonMesh, typename VertexPointMap>
std::size_t stitch_halfedge_range_dispatcher(const HalfedgePairRange& to_stitch_const,
PolygonMesh& pmesh,
@@ -855,20 +861,182 @@ std::size_t stitch_halfedge_range_dispatcher(const HalfedgePairRange& to_stitch_
return stitch_halfedge_range(to_stitch, pmesh, vpm);
}
// collect_duplicated_stitchable_boundary_edges() cannot handle configurations with non-manifoldness.
// However, even if non-manifoldness exists within a loop, it is safe choice to stitch consecutive
// stitchable halfedges
template <typename HalfedgeRange,
typename PolygonMesh,
typename VPM,
typename HalfedgeKeeper>
std::size_t zip_boundary_cycle(typename boost::graph_traits<PolygonMesh>::halfedge_descriptor& bh,
const HalfedgeRange& cycle_halfedges,
PolygonMesh& pmesh,
const VPM vpm,
const HalfedgeKeeper& hd_kpr)
{
typedef typename boost::graph_traits<PolygonMesh>::halfedge_descriptor halfedge_descriptor;
std::size_t stitched_boundary_cycles_n = 0;
// Zipping cannot change the topology of the hole so the maintenance is trivial
internal::Dummy_cycle_rep_maintainer<PolygonMesh> dummy_maintainer(pmesh);
// A boundary cycle might need to be stitched starting from different extremities
//
// v11 ------ v10
// | |
// v0 --- v1(v13) === v2(v12) v5(v9) === v6(v8) --- v7
// | |
// v3 ------- v4
//
// As long as we find vertices on the boundary with both incident halfedges being compatible,
// we zip it up as much as possible.
// not everything is always stitchable
std::set<halfedge_descriptor> unstitchable_halfedges;
const halfedge_descriptor null_h = boost::graph_traits<PolygonMesh>::null_halfedge();
for(;;) // until there is nothing to stitch anymore
{
if(bh == null_h) // the complete boundary cycle is stitched
break;
#ifdef CGAL_PMP_STITCHING_DEBUG
std::cout << "Walking border from halfedge: " << edge(bh, pmesh) << std::endl;
#endif
CGAL_assertion(is_border(bh, pmesh));
halfedge_descriptor hn = next(bh, pmesh), start_h = null_h;
do
{
halfedge_descriptor hnn = next(hn, pmesh);
CGAL_assertion(get(vpm, target(hn, pmesh)) == get(vpm, source(hnn, pmesh)));
if(get(vpm, source(hn, pmesh)) == get(vpm, target(hnn, pmesh)) &&
!is_degenerate_edge(edge(hn, pmesh), pmesh, parameters::vertex_point_map(vpm)))
{
if(unstitchable_halfedges.count(hn) == 0)
{
start_h = hn;
break;
}
}
hn = hnn;
}
while(hn != bh);
if(start_h == null_h) // nothing to be stitched on this boundary cycle
break;
#ifdef CGAL_PMP_STITCHING_DEBUG_PP
std::cout << "Starting stitching from halfedge: "
<< get(vpm, source(edge(start_h, pmesh), pmesh)) << " "
<< get(vpm, target(edge(start_h, pmesh), pmesh)) << std::endl;
#endif
CGAL_assertion(is_border(start_h, pmesh));
// Associate as many consecutive halfedge pairs as possible ("zipping")
std::vector<std::pair<halfedge_descriptor, halfedge_descriptor> > hedges_to_stitch;
halfedge_descriptor curr_h = start_h;
halfedge_descriptor curr_hn = next(curr_h, pmesh);
for(;;) // while we can expand the zipping range
{
// Don't create an invalid polygon mesh, even if the geometry allows it
if(face(opposite(curr_h, pmesh), pmesh) == face(opposite(curr_hn, pmesh), pmesh))
{
unstitchable_halfedges.insert(curr_h);
bh = curr_hn;
break;
}
CGAL_assertion(is_border(curr_h, pmesh));
CGAL_assertion(is_border(curr_hn, pmesh));
if(hd_kpr(curr_h, curr_hn) == curr_h)
hedges_to_stitch.emplace_back(curr_h, curr_hn);
else
hedges_to_stitch.emplace_back(curr_hn, curr_h);
#ifdef CGAL_PMP_STITCHING_DEBUG_PP
std::cout << "expand zip with:\n"
<< edge(curr_h, pmesh) << "\n\t" << source(curr_h, pmesh) << "\t(" << get(vpm, source(curr_h, pmesh)) << ")"
<< "\n\t" << target(curr_h, pmesh) << "\t(" << get(vpm, target(curr_h, pmesh)) << ")\n"
<< edge(curr_hn, pmesh) << "\n\t" << source(curr_hn, pmesh) << "\t(" << get(vpm, source(curr_hn, pmesh)) << ")"
<< "\n\t" << target(curr_hn, pmesh) << "\t(" << get(vpm, target(curr_hn, pmesh)) << ")" << std::endl;
#endif
// check if we have reached the end of the boundary cycle
if(prev(curr_h, pmesh) == curr_hn || prev(curr_h, pmesh) == next(curr_hn, pmesh))
{
bh = null_h;
break;
}
curr_h = prev(curr_h, pmesh);
curr_hn = next(curr_hn, pmesh);
// check if the next two halfedges are not geometrically compatible
if(get(vpm, source(curr_h, pmesh)) != get(vpm, target(curr_hn, pmesh)) ||
is_degenerate_edge(edge(curr_hn, pmesh), pmesh, parameters::vertex_point_map(vpm)))
{
bh = curr_hn;
break;
}
}
// bh must be a boundary halfedge on the border that will not be impacted by any stitching
CGAL_assertion_code(if(bh != null_h) {)
CGAL_assertion_code( for(const auto& hp : hedges_to_stitch) {)
CGAL_assertion( bh != hp.first && bh != hp.second);
CGAL_assertion_code(}})
if(!hedges_to_stitch.empty())
{
#ifdef CGAL_PMP_STITCHING_DEBUG_PP
std::cout << hedges_to_stitch.size() " halfedge pairs to stitch on border containing:\n"
<< edge(h, pmesh) << "\n\t" << source(h, pmesh) << "\t(" << get(vpm, source(h, pmesh)) << ")"
<< "\n\t" << target(h, pmesh) << "\t(" << get(vpm, target(h, pmesh)) << ")" << std::endl;
#endif
std::size_t local_stitches = internal::stitch_halfedge_range(hedges_to_stitch, cycle_halfedges,
pmesh, vpm, dummy_maintainer);
stitched_boundary_cycles_n += local_stitches;
if(local_stitches == 0) // refused to stitch this halfedge pair range due to manifold issue
{
#ifdef CGAL_PMP_STITCHING_DEBUG_PP
std::cout << "Failed to stitch this range!" << std::endl;
#endif
for(const auto& hp : hedges_to_stitch)
{
unstitchable_halfedges.insert(hp.first);
unstitchable_halfedges.insert(hp.second);
}
}
}
}
return stitched_boundary_cycles_n;
}
/// High-level functions
template <typename PolygonMesh, typename MaintainerVisitor, typename CGAL_PMP_NP_TEMPLATE_PARAMETERS>
std::size_t stitch_boundary_cycle(const typename boost::graph_traits<PolygonMesh>::halfedge_descriptor bh,
template <typename PolygonMesh, typename CycleRepMaintainer, typename CGAL_PMP_NP_TEMPLATE_PARAMETERS>
std::size_t stitch_boundary_cycle(const typename boost::graph_traits<PolygonMesh>::halfedge_descriptor h,
PolygonMesh& pmesh,
MaintainerVisitor& mv,
CycleRepMaintainer& cycle_reps_maintainer,
const CGAL_PMP_NP_CLASS& np)
{
typedef typename boost::graph_traits<PolygonMesh>::halfedge_descriptor halfedge_descriptor;
typedef typename std::pair<halfedge_descriptor, halfedge_descriptor> halfedges_pair;
CGAL_precondition(bh != boost::graph_traits<PolygonMesh>::null_halfedge());
CGAL_precondition(is_border(bh, pmesh));
CGAL_precondition(h != boost::graph_traits<PolygonMesh>::null_halfedge());
CGAL_precondition(is_border(h, pmesh));
CGAL_precondition(is_valid(pmesh));
using parameters::choose_parameter;
@@ -884,18 +1052,38 @@ std::size_t stitch_boundary_cycle(const typename boost::graph_traits<PolygonMesh
const Halfedge_keeper hd_kpr = choose_parameter(get_parameter(np, internal_np::halfedges_keeper),
Default_halfedges_keeper<PolygonMesh>());
halfedge_descriptor bh = h, bh_mem = bh;
std::vector<halfedge_descriptor> cycle_halfedges;
for(halfedge_descriptor h : halfedges_around_face(bh, pmesh))
cycle_halfedges.push_back(h);
std::size_t res = internal::zip_boundary_cycle(bh, cycle_halfedges, pmesh, vpm, hd_kpr);
if(bh == boost::graph_traits<PolygonMesh>::null_halfedge()) // stitched everything
{
cycle_reps_maintainer.remove_representative(bh);
return res;
}
// Re-compute the range if something was stitched
if(res != 0)
{
cycle_reps_maintainer.remove_representative(bh_mem);
cycle_reps_maintainer.add_representative(bh);
cycle_halfedges.clear();
for(halfedge_descriptor h : halfedges_around_face(bh, pmesh))
cycle_halfedges.push_back(h);
}
std::vector<halfedges_pair> to_stitch;
internal::collect_duplicated_stitchable_boundary_edges(cycle_halfedges, pmesh,
hd_kpr, false /*per cc*/,
std::back_inserter(to_stitch), np);
mv.remove_representative(bh);
res += stitch_halfedge_range(to_stitch, cycle_halfedges, pmesh, vpm, cycle_reps_maintainer);
return stitch_halfedge_range(to_stitch, cycle_halfedges, pmesh, mv, vpm);
return res;
}
} //end of namespace internal
@@ -920,7 +1108,7 @@ std::size_t stitch_boundary_cycle(const typename boost::graph_traits<PolygonMesh
/// \cgalParamDescription{a property map associating points to the vertices of `pm`}
/// \cgalParamType{a class model of `ReadWritePropertyMap` with `boost::graph_traits<PolygonMesh>::%vertex_descriptor`
/// as key type and `%Point_3` as value type}
/// \cgalParamDefault{`boost::get(CGAL::vertex_point, pm)`}
/// \cgalParamDefault{`boost::get(CGAL::vertex_point, pmesh)`}
/// \cgalParamExtra{If this parameter is omitted, an internal property map for `CGAL::vertex_point_t`
/// must be available in `PolygonMesh`.}
/// \cgalParamNEnd
@@ -936,8 +1124,8 @@ std::size_t stitch_boundary_cycle(const typename boost::graph_traits<PolygonMesh
PolygonMesh& pmesh,
const CGAL_PMP_NP_CLASS& np)
{
internal::Dummy_cycle_rep_maintainer<PolygonMesh> mv(pmesh);
return internal::stitch_boundary_cycle(h, pmesh, mv, np);
internal::Dummy_cycle_rep_maintainer<PolygonMesh> dummy_maintainer(pmesh);
return internal::stitch_boundary_cycle(h, pmesh, dummy_maintainer, np);
}
template <typename PolygonMesh>
@@ -950,17 +1138,17 @@ std::size_t stitch_boundary_cycle(const typename boost::graph_traits<PolygonMesh
namespace internal {
template <typename BorderHalfedgeRange, typename PolygonMesh,
typename MaintainerVisitor, typename CGAL_PMP_NP_TEMPLATE_PARAMETERS>
typename CycleRepMaintainer, typename CGAL_PMP_NP_TEMPLATE_PARAMETERS>
std::size_t stitch_boundary_cycles(const BorderHalfedgeRange& boundary_cycle_representatives,
PolygonMesh& pmesh,
MaintainerVisitor& mv,
CycleRepMaintainer& cycle_reps_maintainer,
const CGAL_PMP_NP_CLASS& np)
{
typedef typename boost::graph_traits<PolygonMesh>::halfedge_descriptor halfedge_descriptor;
std::size_t stitched_boundary_cycles_n = 0;
for(const halfedge_descriptor h : boundary_cycle_representatives)
stitched_boundary_cycles_n += stitch_boundary_cycle(h, pmesh, mv, np);
stitched_boundary_cycles_n += stitch_boundary_cycle(h, pmesh, cycle_reps_maintainer, np);
return stitched_boundary_cycles_n;
}
@@ -988,7 +1176,7 @@ std::size_t stitch_boundary_cycles(const BorderHalfedgeRange& boundary_cycle_rep
/// \cgalParamDescription{a property map associating points to the vertices of `pm`}
/// \cgalParamType{a class model of `ReadWritePropertyMap` with `boost::graph_traits<PolygonMesh>::%vertex_descriptor`
/// as key type and `%Point_3` as value type}
/// \cgalParamDefault{`boost::get(CGAL::vertex_point, pm)`}
/// \cgalParamDefault{`boost::get(CGAL::vertex_point, pmesh)`}
/// \cgalParamExtra{If this parameter is omitted, an internal property map for `CGAL::vertex_point_t`
/// must be available in `PolygonMesh`.}
/// \cgalParamNEnd
@@ -1006,8 +1194,8 @@ std::size_t stitch_boundary_cycles(const BorderHalfedgeRange& boundary_cycle_rep
{
// If this API is called, we are not from stitch_borders() (otherwise there would be a maintainer)
// so there is only one pass and we don't carea bout maintaining the cycle subset
internal::Dummy_cycle_rep_maintainer<PolygonMesh> mv(pmesh);
return stitch_boundary_cycles(boundary_cycle_representatives, pmesh, mv, np);
internal::Dummy_cycle_rep_maintainer<PolygonMesh> dummy_maintainer(pmesh);
return stitch_boundary_cycles(boundary_cycle_representatives, pmesh, dummy_maintainer, np);
}
///\cond SKIP_IN_MANUAL
@@ -1102,11 +1290,11 @@ std::size_t stitch_borders(PolygonMesh& pmesh,
namespace internal {
template <typename BorderHalfedgeRange, typename PolygonMesh,
typename MaintainerVisitor,
typename CycleRepMaintainer,
typename CGAL_PMP_NP_TEMPLATE_PARAMETERS>
std::size_t stitch_borders(const BorderHalfedgeRange& boundary_cycle_representatives,
PolygonMesh& pmesh,
MaintainerVisitor& mv,
CycleRepMaintainer& cycle_maintainer,
const CGAL_PMP_NP_CLASS& np)
{
typedef typename boost::graph_traits<PolygonMesh>::halfedge_descriptor halfedge_descriptor;
@@ -1130,35 +1318,34 @@ std::size_t stitch_borders(const BorderHalfedgeRange& boundary_cycle_representat
bool per_cc = choose_parameter(get_parameter(np, internal_np::apply_per_connected_component), false);
#ifdef CGAL_PMP_STITCHING_DEBUG
std::cout << "------- Stitch cycles... (" << boundary_cycle_representatives.size() << " cycle(s))" << std::endl;
std::cout << "------- Stitch cycles (#1)... (" << boundary_cycle_representatives.size() << " cycle(s))" << std::endl;
#endif
std::size_t res = stitch_boundary_cycles(boundary_cycle_representatives, pmesh, mv, np);
std::size_t res = stitch_boundary_cycles(boundary_cycle_representatives, pmesh, cycle_maintainer, np);
#ifdef CGAL_PMP_STITCHING_DEBUG
std::cout << "------- Stitched " << res << " in boundary cycles" << std::endl;
std::cout << "------- Stitched " << res << " halfedge pairs in boundary cycles" << std::endl;
std::cout << "------- Stitch all..." << std::endl;
#endif
const auto& to_consider = mv.halfedges_to_consider();
mv.clear_representatives();
const auto& to_consider = cycle_maintainer.halfedges_to_consider();
cycle_maintainer.clear_representatives();
std::vector<std::pair<halfedge_descriptor, halfedge_descriptor> > to_stitch;
internal::collect_duplicated_stitchable_boundary_edges(to_consider, pmesh, hd_kpr, per_cc,
std::back_inserter(to_stitch), np);
res += stitch_halfedge_range(to_stitch, to_consider, pmesh, mv, vpm);
const auto& new_representatives = mv.cycle_representatives();
res += stitch_halfedge_range(to_stitch, to_consider, pmesh, vpm, cycle_maintainer);
#ifdef CGAL_PMP_STITCHING_DEBUG
std::cout << "------- Stitched " << res << " after cycles & general" << std::endl;
std::cout << "------- Stitch cycles (#2)... (" << new_representatives.size() << " cycles)" << std::endl;
std::cout << "------- Stitched " << res << " halfedge pairs after cycles & general" << std::endl;
std::cout << "------- Stitch cycles (#2)... (" << new_representatives.size() << " cycle(s))" << std::endl;
#endif
const auto& new_representatives = cycle_maintainer.cycle_representatives();
// Don't care about keeping track of the sub-cycles as this is the last pass
internal::Dummy_cycle_rep_maintainer<PolygonMesh> null_mv(pmesh);
res += stitch_boundary_cycles(new_representatives, pmesh, null_mv, np);
internal::Dummy_cycle_rep_maintainer<PolygonMesh> dummy_cycle_maintainer(pmesh);
res += stitch_boundary_cycles(new_representatives, pmesh, dummy_cycle_maintainer, np);
#ifdef CGAL_PMP_STITCHING_DEBUG
std::cout << "------- Stitched " << res << " (total)" << std::endl;
@@ -1228,8 +1415,8 @@ std::size_t stitch_borders(const BorderHalfedgeRange& boundary_cycle_representat
)
{
// Need to keep track of the cycles since we are working on a subset of all the boundary cycles
internal::Boundary_cycle_rep_maintainer<PolygonMesh> mv(pmesh);
return stitch_borders(boundary_cycle_representatives, pmesh, mv, np);
internal::Boundary_cycle_rep_maintainer<PolygonMesh> cycle_reps_maintainer(pmesh);
return stitch_borders(boundary_cycle_representatives, pmesh, cycle_reps_maintainer, np);
}
/// \cond SKIP_IN_MANUAL
@@ -1242,8 +1429,8 @@ std::size_t stitch_borders(const BorderHalfedgeRange& boundary_cycle_representat
>::type* = 0)
{
// Need to keep track of the cycles since we are working on a subset of all the boundary cycles
internal::Boundary_cycle_rep_maintainer<PolygonMesh> mv(pmesh);
return stitch_borders(boundary_cycle_representatives, pmesh, mv, parameters::all_default());
internal::Boundary_cycle_rep_maintainer<PolygonMesh> cycle_reps_maintainer(pmesh);
return stitch_borders(boundary_cycle_representatives, pmesh, cycle_reps_maintainer, parameters::all_default());
}
template <typename PolygonMesh, typename CGAL_PMP_NP_TEMPLATE_PARAMETERS>
@@ -1256,8 +1443,8 @@ std::size_t stitch_borders(PolygonMesh& pmesh,
extract_boundary_cycles(pmesh, std::back_inserter(boundary_cycle_representatives));
// We are working on all boundary cycles, so there is no need to keep track of any subset
internal::Dummy_cycle_rep_maintainer<PolygonMesh> mv(pmesh);
return stitch_borders(boundary_cycle_representatives, pmesh, mv, np);
internal::Dummy_cycle_rep_maintainer<PolygonMesh> dummy_maintainer(pmesh);
return stitch_borders(boundary_cycle_representatives, pmesh, dummy_maintainer, np);
}
template <typename PolygonMesh>
@@ -729,7 +729,7 @@ private:
for (unsigned int i = 0; i < cutp.size(); i++){
const Plane& plane_i = prism[cutp[i]];
boost::optional<ePoint_3> op = intersection_point(line, plane_i.eplane);
boost::optional<ePoint_3> op = intersection_point_for_polyhedral_envelope(line, plane_i.eplane);
if(! op){
std::cout << "there must be an intersection 2" << std::endl;
}
@@ -850,8 +850,8 @@ private:
}
for (unsigned int j = 0; j < cidl.size(); j++) {
boost::optional<ePoint_3> op = intersection_point(line,
halfspace[prismindex[queue[i]]][cidl[j]].eplane);
boost::optional<ePoint_3> op = intersection_point_for_polyhedral_envelope(line,
halfspace[prismindex[queue[i]]][cidl[j]].eplane);
const ePoint_3& ip = *op;
inter = Implicit_Seg_Facet_interpoint_Out_Prism_return_local_id
(ip, idlist, jump1, check_id);
@@ -1034,7 +1034,7 @@ private:
const Plane& plane_i = prism[cutp[i]];
const eLine_3& eline = *(seg[k]);
boost::optional<ePoint_3> op = intersection_point(eline, plane_i.eplane);
boost::optional<ePoint_3> op = intersection_point_for_polyhedral_envelope(eline, plane_i.eplane);
if(! op){
#ifdef CGAL_ENVELOPE_DEBUG
std::cout << "there must be an intersection 6" << std::endl;
@@ -1088,7 +1088,7 @@ private:
int inter = 0;
boost::optional<ePoint_3> ipp = intersection_point(tri_eplane, prism[cutp[i]].eplane, prism[cutp[j]].eplane);
boost::optional<ePoint_3> ipp = intersection_point_for_polyhedral_envelope(tri_eplane, prism[cutp[i]].eplane, prism[cutp[j]].eplane);
if(ipp){
inter = is_3_triangle_cut_float_fast(tri0, tri1, tri2,
n,
@@ -1528,8 +1528,8 @@ private:
if (!cut) continue;
for (unsigned int j = 0; j < cidl.size(); j++) {
boost::optional<ePoint_3> op = intersection_point(eline,
halfspace[prismindex[queue[i]]][cidl[j]].eplane);
boost::optional<ePoint_3> op = intersection_point_for_polyhedral_envelope(eline,
halfspace[prismindex[queue[i]]][cidl[j]].eplane);
const ePoint_3& ip = *op;
inter = Implicit_Seg_Facet_interpoint_Out_Prism_return_local_id(ip, idlist, jump1, check_id);
@@ -1612,8 +1612,8 @@ private:
}
// now we know that there exists an intesection point
boost::optional<ePoint_3> op = intersection_point(eline,
halfspace[filtered_intersection[queue[i]]][intersect_face[queue[i]][j]].eplane);
boost::optional<ePoint_3> op = intersection_point_for_polyhedral_envelope(eline,
halfspace[filtered_intersection[queue[i]]][intersect_face[queue[i]][j]].eplane);
const ePoint_3& ip = *op;
inter = Implicit_Seg_Facet_interpoint_Out_Prism_return_local_id_with_face_order(ip, idlist, idlistorder, jump1, check_id);
@@ -1696,9 +1696,9 @@ private:
// We moved the intersection here
// In case there is no intersection point we continue
boost::optional<ePoint_3>
op = intersection_point(etriangle_eplane,
halfspace[jump1][intersect_face[queue[i]][k]].eplane,
halfspace[jump2][intersect_face[queue[j]][h]].eplane);
op = intersection_point_for_polyhedral_envelope(etriangle_eplane,
halfspace[jump1][intersect_face[queue[i]][k]].eplane,
halfspace[jump2][intersect_face[queue[j]][h]].eplane);
if(! op){
continue;
}
@@ -176,6 +176,12 @@ public:
: m_free_id(0)
{}
//! move constructor
Rigid_triangle_mesh_collision_detection(Rigid_triangle_mesh_collision_detection&& other)
{
*this = std::move(other);
}
~Rigid_triangle_mesh_collision_detection()
{
for (std::size_t k=0; k<m_free_id; ++k)
@@ -185,6 +191,23 @@ public:
}
}
//! move assignment operator
Rigid_triangle_mesh_collision_detection& operator=(Rigid_triangle_mesh_collision_detection&& other)
{
m_own_aabb_trees = std::move(other.m_own_aabb_trees);
m_aabb_trees = std::move(other.m_aabb_trees);
m_is_closed = std::move(other.m_is_closed);
m_points_per_cc = std::move(other.m_points_per_cc);
m_traversal_traits = std::move(other.m_traversal_traits);
m_free_id = std::move(other.m_free_id);
m_id_pool = std::move(other.m_id_pool);
for(std::size_t i = 0; i< other.m_own_aabb_trees.size(); ++i)
other.m_own_aabb_trees[i]= false;
return *this;
}
/*!
* adds mesh `tm` to the set of meshes to be considered for intersection.
*