Merge pull request #5509 from sloriot/PMP-coref_robust_side_of

Robustify internal side-of-triangle-mesh calls
This commit is contained in:
Laurent Rineau 2021-03-03 17:18:11 +01:00
commit fd353e84e9
3 changed files with 250 additions and 26 deletions

View File

@ -428,6 +428,11 @@ public:
Closest_point closest_point_object() const {return Closest_point(*this);}
Compare_distance compare_distance_object() const {return Compare_distance();}
typedef enum { CGAL_AXIS_X = 0,
CGAL_AXIS_Y = 1,
CGAL_AXIS_Z = 2} Axis;
static Axis longest_axis(const Bounding_box& bbox);
private:
/**
@ -446,13 +451,6 @@ private:
return internal::Primitive_helper<AT>::get_datum(pr,*this).bbox();
}
typedef enum { CGAL_AXIS_X = 0,
CGAL_AXIS_Y = 1,
CGAL_AXIS_Z = 2} Axis;
static Axis longest_axis(const Bounding_box& bbox);
/// Comparison functions
static bool less_x(const Primitive& pr1, const Primitive& pr2,const AABB_traits<GeomTraits,AABBPrimitive, BboxMap>& traits)
{

View File

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

View File

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