diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Side_of_triangle_mesh/Ray_3_Triangle_3_traversal_traits.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Side_of_triangle_mesh/Ray_3_Triangle_3_traversal_traits.h index 285017791a5..130c2732ea1 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Side_of_triangle_mesh/Ray_3_Triangle_3_traversal_traits.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Side_of_triangle_mesh/Ray_3_Triangle_3_traversal_traits.h @@ -216,6 +216,60 @@ public: } }; +//special case when ray query is from another Kernel K1 is the kernel compatible with the AABB-tree +template +class K2_Ray_3_K1_Triangle_3_traversal_traits +{ + //the status indicates whether the query point is strictly inside the polyhedron, and the number of intersected triangles if yes + std::pair& m_status; + bool m_stop; + const AABBTraits& m_aabb_traits; + typedef typename AABBTraits::Primitive Primitive; + typedef CGAL::AABB_node Node; + Helper m_helper; + CGAL::Cartesian_converter to_K2; + +public: + K2_Ray_3_K1_Triangle_3_traversal_traits(std::pair& status, + const AABBTraits& aabb_traits, + const Helper& h) + :m_status(status), m_stop(false), m_aabb_traits(aabb_traits), m_helper(h) + {m_status.first=true;} + + bool go_further() const { return !m_stop; } + + template + void intersection(const Query& query, const Primitive& primitive) + { + Intersections::internal::r3t3_do_intersect_endpoint_position_visitor visitor; + std::pair res= + Intersections::internal::do_intersect(to_K2(m_helper.get_primitive_datum(primitive, m_aabb_traits)), + query, K2(), visitor); + + if (res.first){ + switch (res.second){ + case Intersections::internal::R3T3_intersection::CROSS_FACET: + ++m_status.second; + break; + case Intersections::internal::R3T3_intersection::ENDPOINT_IN_TRIANGLE: + m_status.first=false; + m_stop=true; + break; + default: + m_status.first=boost::logic::indeterminate; + m_stop=true; + } + } + } + + template + bool do_intersect(const Query& query, const Node& node) const + { + return CGAL::do_intersect(query, m_helper.get_node_bbox(node)); + } +}; + + }// namespace internal }// namespace CGAL diff --git a/Polygon_mesh_processing/include/CGAL/Side_of_triangle_mesh.h b/Polygon_mesh_processing/include/CGAL/Side_of_triangle_mesh.h index 20ec968d50a..f0c23942f60 100644 --- a/Polygon_mesh_processing/include/CGAL/Side_of_triangle_mesh.h +++ b/Polygon_mesh_processing/include/CGAL/Side_of_triangle_mesh.h @@ -242,6 +242,81 @@ public: } } +#ifndef DOXYGEN_RUNNING + template + Bounded_side operator()(const typename K2::Point_3& point, const K2& k2) const + { + if(point.x() < box.xmin() + || point.x() > box.xmax() + || point.y() < box.ymin() + || point.y() > box.ymax() + || point.z() < box.zmin() + || point.z() > box.zmax()) + { + return CGAL::ON_UNBOUNDED_SIDE; + } + +#ifdef CGAL_HAS_THREADS + AABB_tree_* tree_ptr = + const_cast(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(atomic_tree_ptr.load(std::memory_order_relaxed)); +#endif + CGAL_assertion(tm_ptr != nullptr && opt_vpm!=boost::none); + if (tree_ptr==nullptr) + { + tree_ptr = new AABB_tree(faces(*tm_ptr).first, + faces(*tm_ptr).second, + *tm_ptr, *opt_vpm); + const_cast(tree_ptr)->build(); +#ifdef CGAL_HAS_THREADS + atomic_tree_ptr.store(tree_ptr, std::memory_order_release); +#endif + } + } + + typedef typename Kernel_traits::Kernel K1; + typedef typename AABB_tree::AABB_traits AABB_traits; + typedef internal::Default_tree_helper Helper; + Helper helper; + typename AABB_traits::Bounding_box bbox = helper.get_tree_bbox(*tree_ptr); + + static const unsigned int seed = 1340818006; + CGAL::Random rg(seed); // seed some value for make it easy to debug + Random_points_on_sphere_3 random_point(1.,rg); + + typename K2::Construct_ray_3 ray = k2.construct_ray_3_object(); + typename K2::Construct_vector_3 vector = k2.construct_vector_3_object(); + + do { //retry with a random ray + typename K2::Ray_3 query = ray(point, vector(CGAL::ORIGIN,*random_point++)); + + std::pair + status( boost::logic::tribool(boost::logic::indeterminate), 0); + + internal::K2_Ray_3_K1_Triangle_3_traversal_traits + traversal_traits(status, tree_ptr->traits(), helper); + + tree_ptr->traversal(query, traversal_traits); + + if ( !boost::logic::indeterminate(status.first) ) + { + if (status.first) + return (status.second&1) == 1 ? ON_BOUNDED_SIDE : ON_UNBOUNDED_SIDE; + //otherwise the point is on the facet + return ON_BOUNDARY; + } + } while (true); + return ON_BOUNDARY; // should never be reached + } + +#endif + }; } // namespace CGAL