diff --git a/Mesh_3/include/CGAL/Polyhedral_mesh_domain_3.h b/Mesh_3/include/CGAL/Polyhedral_mesh_domain_3.h index 933c1dc21ce..8d445b29f81 100644 --- a/Mesh_3/include/CGAL/Polyhedral_mesh_domain_3.h +++ b/Mesh_3/include/CGAL/Polyhedral_mesh_domain_3.h @@ -27,6 +27,8 @@ #ifndef CGAL_POLYHEDRAL_MESH_DOMAIN_3_H #define CGAL_POLYHEDRAL_MESH_DOMAIN_3_H +#include + #include #include #include @@ -718,47 +720,12 @@ Polyhedral_mesh_domain_3:: Is_in_domain::operator()(const Point_3& p) const { if(r_domain_.bounding_tree_ == 0) return Subdomain(); - const Bounding_box& bbox = r_domain_.bounding_tree_->bbox(); - if( p.x() < bbox.xmin() || p.x() > bbox.xmax() - || p.y() < bbox.ymin() || p.y() > bbox.ymax() - || p.z() < bbox.zmin() || p.z() > bbox.zmax() ) - { - return Subdomain(); - } - -#ifdef CGAL_POLYHEDRAL_MESH_DOMAIN_USE_GRID - Vector_3 v = p - r_domain_.grid_base; - int i = boost::math::round(v.x() / r_domain_.grid_dx); - int j = boost::math::round(v.y() / r_domain_.grid_dy); - int k = boost::math::round(v.z() / r_domain_.grid_dz); - if(i>19)i=19; - if(j>19)j=19; - if(k>19)k=19; - int index = i*400 + j*20 + k; - - const std::pair& close_point = r_domain_.grid[index]; - typename IGT::Construct_segment_3 segment = IGT().construct_segment_3_object(); - const Segment_3 query = segment(p, close_point.first); - typename AABB_tree::size_type M = (close_point.second)? 0 : 1; - -#else - - typename IGT::Construct_ray_3 ray = IGT().construct_ray_3_object(); - typename IGT::Construct_vector_3 vector = IGT().construct_vector_3_object(); - - Random_points_on_sphere_3 random_point(1.); - - const Ray_3 query = ray(p, vector(CGAL::ORIGIN,*random_point)); - typename AABB_tree::size_type M = 1; - -#endif - - if ( (r_domain_.bounding_tree_->number_of_intersected_primitives(query)&1) == M ) - return Subdomain(Subdomain_index(1)); - else - return Subdomain(); + internal::Point_inside_vertical_ray_cast inside_functor; + Bounded_side side = inside_functor(p, *(r_domain_.bounding_tree_)); + if(side == CGAL::ON_UNBOUNDED_SIDE) { return Subdomain(); } + else { return Subdomain(Subdomain_index(1)); } // case ON_BOUNDARY && ON_BOUNDED_SIDE } diff --git a/Operations_on_polyhedra/include/CGAL/Point_inside_polyhedron_3.h b/Operations_on_polyhedra/include/CGAL/Point_inside_polyhedron_3.h index 5b9da2df4d8..a2b03e0b35a 100644 --- a/Operations_on_polyhedra/include/CGAL/Point_inside_polyhedron_3.h +++ b/Operations_on_polyhedra/include/CGAL/Point_inside_polyhedron_3.h @@ -22,16 +22,13 @@ #ifndef CGAL_POINT_INSIDE_POLYHEDRON_H #define CGAL_POINT_INSIDE_POLYHEDRON_H -#include +#include #include #include #include -#include #include -#include - namespace CGAL { /** @@ -57,17 +54,14 @@ class Point_inside_polyhedron_3{ // typedefs typedef CGAL::internal::AABB_triangle_accessor_3_primitive Primitive; typedef CGAL::AABB_traits Traits; - typedef typename Traits::Bounding_box Bounding_box; typedef CGAL::AABB_tree Tree; typedef typename Kernel::Point_3 Point; - typedef typename Kernel::Ray_3 Ray; + //members typename Kernel::Construct_ray_3 ray_functor; typename Kernel::Construct_vector_3 vector_functor; Tree tree; - const static unsigned int seed = 1340818006; - public: /** * Default constructor. The domain is considered to be empty. @@ -137,53 +131,8 @@ public: */ Bounded_side operator()(const Point& point) const { - const Bounding_box& bbox = tree.bbox(); - - if( point.x() < bbox.xmin() || point.x() > bbox.xmax() - || point.y() < bbox.ymin() || point.y() > bbox.ymax() - || point.z() < bbox.zmin() || point.z() > bbox.zmax() ) - { - return ON_UNBOUNDED_SIDE; - } - - //the direction of the vertical ray depends on the position of the point in the bbox - //in order to limit the expected number of nodes visited. - Ray query = ray_functor(point, vector_functor(0,0,(2*point.z() < tree.bbox().zmax()+tree.bbox().zmin()?-1:1))); - boost::optional res = is_inside_ray_tree_traversal(query); - - if(!res) { - CGAL::Random rg(seed); // seed some value for make it easy to debug - Random_points_on_sphere_3 random_point(1.,rg); - - do { //retry with a random ray - query = ray_functor(point, vector_functor(CGAL::ORIGIN,*random_point++)); - res = is_inside_ray_tree_traversal(query); - } while (!res); - } - return *res; + return internal::Point_inside_vertical_ray_cast()(point, tree, ray_functor, vector_functor); } - -private: - template - boost::optional - is_inside_ray_tree_traversal(const Query& query) const - { - std::pair status( boost::logic::tribool(boost::logic::indeterminate), 0); - - internal::Ray_3_Triangle_3_traversal_traits > traversal_traits(status); - tree.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; - } - return boost::optional(); // indeterminate - } - }; } // namespace CGAL diff --git a/Operations_on_polyhedra/include/CGAL/internal/Operations_on_polyhedra/Point_inside_vertical_ray_cast.h b/Operations_on_polyhedra/include/CGAL/internal/Operations_on_polyhedra/Point_inside_vertical_ray_cast.h new file mode 100644 index 00000000000..12803897282 --- /dev/null +++ b/Operations_on_polyhedra/include/CGAL/internal/Operations_on_polyhedra/Point_inside_vertical_ray_cast.h @@ -0,0 +1,80 @@ +#ifndef CGAL_POINT_INSIDE_POLYHEDRON_POINT_INSIDE_VERTICAL_RAY_CAST_H +#define CGAL_POINT_INSIDE_POLYHEDRON_POINT_INSIDE_VERTICAL_RAY_CAST_H + +#include +#include + +#include + +namespace CGAL { +namespace internal { + +// internal class for point inside test, using existing AABB tree +template +class Point_inside_vertical_ray_cast +{ + typedef typename Kernel::Point_3 Point; + typedef typename Kernel::Ray_3 Ray; + typedef typename AABBTree::AABB_traits Traits; + + const static unsigned int seed = 1340818006; + +public: + Bounded_side operator()( + const Point& point, + const AABBTree& tree, + typename Kernel::Construct_ray_3 ray_functor = Kernel().construct_ray_3_object(), + typename Kernel::Construct_vector_3 vector_functor = Kernel().construct_vector_3_object() ) const + { + const typename Traits::Bounding_box& bbox = tree.bbox(); + + if( point.x() < bbox.xmin() || point.x() > bbox.xmax() + || point.y() < bbox.ymin() || point.y() > bbox.ymax() + || point.z() < bbox.zmin() || point.z() > bbox.zmax() ) + { + return ON_UNBOUNDED_SIDE; + } + + //the direction of the vertical ray depends on the position of the point in the bbox + //in order to limit the expected number of nodes visited. + Ray query = ray_functor(point, vector_functor(0,0,(2*point.z() < tree.bbox().zmax()+tree.bbox().zmin()?-1:1))); + boost::optional res = is_inside_ray_tree_traversal(query, tree); + + if(!res) { + CGAL::Random rg(seed); // seed some value for make it easy to debug + Random_points_on_sphere_3 random_point(1.,rg); + + do { //retry with a random ray + query = ray_functor(point, vector_functor(CGAL::ORIGIN,*random_point++)); + res = is_inside_ray_tree_traversal(query, tree); + } while (!res); + } + return *res; + } + +private: + template + boost::optional + is_inside_ray_tree_traversal(const Ray& ray, const AABBTree& tree) const + { + std::pair status( boost::logic::tribool(boost::logic::indeterminate), 0); + + Ray_3_Triangle_3_traversal_traits > traversal_traits(status); + tree.traversal(ray, 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; + } + return boost::optional(); // indeterminate + } +}; + +}// namespace internal +}// namespace CGAL + +#endif