From 5cdfd8217c2c6d6ba4bfbb19b7193eb9d3877df8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Thu, 25 Apr 2019 16:07:21 +0200 Subject: [PATCH] Fixes to PMP::locate() --- .../CGAL/Polygon_mesh_processing/locate.h | 511 +++++++++++------- 1 file changed, 320 insertions(+), 191 deletions(-) diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/locate.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/locate.h index 158195af63f..d69a6a32d88 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/locate.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/locate.h @@ -27,7 +27,6 @@ #include #include #include -#include #include #include #include @@ -38,13 +37,16 @@ #include #include -#include #include +#include #include #include #include +#include +#include #include +#include // Everywhere in this file: // If `tm` is the input triangulated surface mesh and given the pair (`f`, `bc`) @@ -72,20 +74,20 @@ struct Locate_types halfedge_descriptor, face_descriptor> descriptor_variant; - typedef typename boost::mpl::if_< - boost::is_same< - NamedParameters, Default>, + typedef typename std::conditional< + std::is_same< + NamedParameters, Default>::value, typename boost::property_map::const_type, + boost::vertex_point_t>::const_type, typename GetVertexPointMap::const_type - >::type VertexPointMap; - typedef typename boost::property_traits::value_type Point; + >::type VPM; + typedef typename boost::property_traits::value_type Point; typedef typename CGAL::Kernel_traits::type Kernel; typedef typename Kernel::FT FT; - typedef CGAL::cpp11::array Barycentric_coordinates; + typedef std::array Barycentric_coordinates; typedef std::pair Face_location; }; @@ -156,10 +158,11 @@ snap_coordinates_to_border(typename Locate_types::Barycentric_coor { typedef typename internal::Locate_types::FT FT; - // @tmp clean that or protect it with a macro/variable -// std::cout << "Pre-snapping: " << coords[0] << " " << coords[1] << " " << coords[2] << std::endl; -// std::cout << "Sum: " << coords[0] + coords[1] + coords[2] << std::endl; -// std::cout << "tolerance: " << tolerance << std::endl; +#ifdef CGAL_PMP_LOCATE_DEBUG + std::cout << "Pre-snapping: " << coords[0] << " " << coords[1] << " " << coords[2] << std::endl; + std::cout << "Sum: " << coords[0] + coords[1] + coords[2] << std::endl; + std::cout << "tolerance: " << tolerance << std::endl; +#endif // To still keep a sum roughly equals to 1, keep in memory the small changes FT residue = 0.; @@ -191,11 +194,12 @@ snap_coordinates_to_border(typename Locate_types::Barycentric_coor } } - // @tmp clean that or protect it with a macro/variable -// std::cout << "Post-snapping: " << coords[0] << " " -// << coords[1] << " " -// << coords[2] << std::endl; -// std::cout << "Sum: " << coords[0] + coords[1] + coords[2] << std::endl; +#ifdef CGAL_PMP_LOCATE_DEBUG + std::cout << "Post-snapping: " << coords[0] << " " + << coords[1] << " " + << coords[2] << std::endl; + std::cout << "Sum: " << coords[0] + coords[1] + coords[2] << std::endl; +#endif return snapped; } @@ -230,6 +234,114 @@ random_entity_in_range(const CGAL::Iterator_range& range, return random_entity_in_range(range.begin(), range.end(), rnd); } +template +struct Barycentric_coordinate_calculator // 2D version +{ + CGAL::array + operator()(const P& ip, const P& iq, const P& ir, const P& iquery, const K& k) const + { + typedef typename K::FT FT; + typedef typename K::Vector_2 Vector_2; + + typename K::Construct_point_2 cp2 = k.construct_point_2_object(); + typename K::Construct_vector_2 cv2 = k.construct_vector_2_object(); + typename K::Compute_scalar_product_2 csp2 = k.compute_scalar_product_2_object(); + + const typename K::Point_2& p = cp2(ip); + const typename K::Point_2& q = cp2(iq); + const typename K::Point_2& r = cp2(ir); + const typename K::Point_2& query = cp2(iquery); + + Vector_2 v0 = cv2(p, q); + Vector_2 v1 = cv2(p, r); + Vector_2 v2 = cv2(p, query); + + FT d00 = csp2(v0, v0); + FT d01 = csp2(v0, v1); + FT d11 = csp2(v1, v1); + FT d20 = csp2(v2, v0); + FT d21 = csp2(v2, v1); + + FT denom = d00 * d11 - d01 * d01; + + FT v = (d11 * d20 - d01 * d21) / denom; + FT w = (d00 * d21 - d01 * d20) / denom; + + return CGAL::make_array(FT(FT(1) - v - w), v, w); + } +}; + +template +struct Barycentric_coordinate_calculator +{ + CGAL::array + operator()(const P& ip, const P& iq, const P& ir, const P& iquery, const K& k) const + { + typedef typename K::FT FT; + typedef typename K::Vector_3 Vector_3; + + typename K::Construct_point_3 cp3 = k.construct_point_3_object(); + typename K::Construct_vector_3 cv3 = k.construct_vector_3_object(); + typename K::Compute_scalar_product_3 csp3 = k.compute_scalar_product_3_object(); + + const typename K::Point_3& p = cp3(ip); + const typename K::Point_3& q = cp3(iq); + const typename K::Point_3& r = cp3(ir); + const typename K::Point_3& query = cp3(iquery); + + Vector_3 v0 = cv3(p, q); + Vector_3 v1 = cv3(p, r); + Vector_3 v2 = cv3(p, query); + + FT d00 = csp3(v0, v0); + FT d01 = csp3(v0, v1); + FT d11 = csp3(v1, v1); + FT d20 = csp3(v2, v0); + FT d21 = csp3(v2, v1); + + CGAL_assertion((d00 * d11 - d01 * d01) != FT(0)); // denom != 0. + FT denom_inv = 1. / (d00 * d11 - d01 * d01); + + FT v = (d11 * d20 - d01 * d21) * denom_inv; + FT w = (d00 * d21 - d01 * d20) * denom_inv; + + return CGAL::make_array(FT(FT(1) - v - w), v, w); + } +}; + +template +struct Barycentric_point_constructor // 2D version +{ + typedef typename K::FT FT; + + P operator()(const P& p, const FT wp, const P& q, const FT wq, const P& r, const FT wr) const + { + FT sum = wp + wq + wr; + CGAL_assertion(sum != 0); + FT x = (wp * p.x() + wq * q.x() + wr * r.x()) / sum; + FT y = (wp * p.y() + wq * q.y() + wr * r.y()) / sum; + + return P(x, y); + } +}; + +template +struct Barycentric_point_constructor // 3D version +{ + typedef typename K::FT FT; + + P operator()(const P& p, const FT wp, const P& q, const FT wq, const P& r, const FT wr) const + { + FT sum = wp + wq + wr; + CGAL_assertion(sum != 0); + FT x = (wp * p.x() + wq * q.x() + wr * r.x()) / sum; + FT y = (wp * p.y() + wq * q.y() + wr * r.y()) / sum; + FT z = (wp * p.z() + wq * q.z() + wr * r.z()) / sum; + + return P(x, y, z); + } +}; + } // namespace internal /// \brief returns a random non-null vertex incident to the face `fd` of the polygon mesh `pm`. @@ -323,8 +435,11 @@ int vertex_index_in_face(const typename boost::graph_traits::vertex halfedge_descriptor current = start; int counter = 0; + std::cout << "vd to find: " << &*vd << " pt: " << vd->point() << std::endl; + do { + std::cout << "current source: " << &*source(current, pm) << " pt: " << source(current, pm)->point() << std::endl; if(source(current, pm) == vd) break; @@ -370,81 +485,6 @@ int halfedge_index_in_face(typename boost::graph_traits::halfedge_d return count; } -template -struct Barycentric_coordinate_calculator // 2D version -{ - CGAL::array - operator()(const P& ip, const P& iq, const P& ir, const P& iquery, const K& k) const - { - typedef typename K::FT FT; - typedef typename K::Vector_2 Vector_2; - - typename K::Construct_point_2 cp2 = k.construct_point_2_object(); - typename K::Construct_vector_2 cv2 = k.construct_vector_2_object(); - typename K::Compute_scalar_product_2 csp2 = k.compute_scalar_product_2_object(); - - const typename K::Point_2& p = cp2(ip); - const typename K::Point_2& q = cp2(iq); - const typename K::Point_2& r = cp2(ir); - const typename K::Point_2& query = cp2(iquery); - - Vector_2 v0 = cv2(p, q); - Vector_2 v1 = cv2(p, r); - Vector_2 v2 = cv2(p, query); - - FT d00 = csp2(v0, v0); - FT d01 = csp2(v0, v1); - FT d11 = csp2(v1, v1); - FT d20 = csp2(v2, v0); - FT d21 = csp2(v2, v1); - - FT denom = d00 * d11 - d01 * d01; - - FT v = (d11 * d20 - d01 * d21) / denom; - FT w = (d00 * d21 - d01 * d20) / denom; - - return CGAL::make_array(FT(1) - v - w, v, w); - } -}; - -template -struct Barycentric_coordinate_calculator -{ - CGAL::array - operator()(const P& ip, const P& iq, const P& ir, const P& iquery, const K& k) const - { - typedef typename K::FT FT; - typedef typename K::Vector_3 Vector_3; - - typename K::Construct_point_3 cp3 = k.construct_point_3_object(); - typename K::Construct_vector_3 cv3 = k.construct_vector_3_object(); - typename K::Compute_scalar_product_3 csp3 = k.compute_scalar_product_3_object(); - - const typename K::Point_3& p = cp3(ip); - const typename K::Point_3& q = cp3(iq); - const typename K::Point_3& r = cp3(ir); - const typename K::Point_3& query = cp3(iquery); - - Vector_3 v0 = cv3(p, q); - Vector_3 v1 = cv3(p, r); - Vector_3 v2 = cv3(p, query); - - FT d00 = csp3(v0, v0); - FT d01 = csp3(v0, v1); - FT d11 = csp3(v1, v1); - FT d20 = csp3(v2, v0); - FT d21 = csp3(v2, v1); - - CGAL_assertion((d00 * d11 - d01 * d01) != FT(0)); // denom != 0. - FT denom_inv = 1. / (d00 * d11 - d01 * d01); - - FT v = (d11 * d20 - d01 * d21) * denom_inv; - FT w = (d00 * d21 - d01 * d20) * denom_inv; - - return CGAL::make_array(FT(1) - v - w, v, w); - } -}; - /// \brief Given a set of three points and a query point, computes the barycentric /// coordinates of the query point with respect to the first three points. /// @@ -453,15 +493,21 @@ struct Barycentric_coordinate_calculator /// (this is the case for all standard %CGAL point types). /// \pre `query` lies on the plane defined by `p`, `q`, and `r`. /// +template +CGAL::array +barycentric_coordinates(const P& p, const P& q, const P& r, const P& query, const K& k) +{ + internal::Barycentric_coordinate_calculator calculator; + return calculator(p, q, r, query, k); +} + template CGAL::array::type::FT, 3> barycentric_coordinates(const P& p, const P& q, const P& r, const P& query) { typedef typename CGAL::Kernel_traits

::type Kernel; - Barycentric_coordinate_calculator calculator; - - return calculator(p, q, r, query, Kernel()); + return barycentric_coordinates(p, q, r, query, Kernel()); } // Random locations @@ -487,7 +533,7 @@ random_location_on_halfedge(typename boost::graph_traits::halfedge CGAL_precondition(CGAL::is_triangle_mesh(tm)); - FT t = rnd.uniform_real(FT(0), FT(1)); + FT t(rnd.uniform_real(0., 1.)); return locate_in_face(hd, t, tm); } @@ -513,8 +559,10 @@ random_location_on_face(typename boost::graph_traits::face_descrip CGAL_precondition(CGAL::is_triangle_mesh(tm)); CGAL_precondition(fd != boost::graph_traits::null_face()); - FT u = rnd.uniform_real(FT(0), FT(1)); - FT v = rnd.uniform_real(FT(0), FT(FT(1) - u)); + // calling 'rnd.uniform real' with double in case FT comes from an EPECK kernel (which doesn't seem to work too well) + FT u(rnd.uniform_real(0., 1.)); + FT v(rnd.uniform_real(0., double(FT(1) - u))); + return std::make_pair(fd, CGAL::make_array(u, v, FT(FT(1) - u - v))); } @@ -613,7 +661,7 @@ get_descriptor_from_location(const typename internal::Locate_types /// \cgalParamBegin{vertex_point_map} /// the property map with the points associated to the vertices of `tm`. /// If this parameter is omitted, an internal property map for -/// `CGAL::vertex_point_t` should be available in `TriangleMesh`. +/// `boost::vertex_point_t` should be available in `TriangleMesh`. /// \cgalParamEnd /// \cgalNamedParamsEnd /// @@ -630,24 +678,24 @@ location_to_point(const typename internal::Locate_types::Face_loca typedef typename GetVertexPointMap::const_type VertexPointMap; typedef typename boost::property_traits::value_type Point; - using boost::choose_param; - using boost::get_param; + typedef typename internal::Locate_types::Kernel Kernel; CGAL_precondition(CGAL::is_triangle_mesh(tm)); - VertexPointMap vpm = choose_param(get_param(np, internal_np::vertex_point), - get_const_property_map(boost::vertex_point, tm)); + VertexPointMap vpm = boost::choose_param(boost::get_param(np, internal_np::vertex_point), + get_const_property_map(boost::vertex_point, tm)); halfedge_descriptor hd = halfedge(loc.first, tm); - const Point p0 = get(vpm, source(hd, tm)); - const Point p1 = get(vpm, target(hd, tm)); - const Point p2 = get(vpm, target(next(hd, tm), tm)); + const Point& p0 = get(vpm, source(hd, tm)); + const Point& p1 = get(vpm, target(hd, tm)); + const Point& p2 = get(vpm, target(next(hd, tm), tm)); - return CGAL::barycenter(p0, loc.second[0], p1, loc.second[1], p2, loc.second[2]); + internal::Barycentric_point_constructor bp_constructor; + return bp_constructor(p0, loc.second[0], p1, loc.second[1], p2, loc.second[2]); } template -typename property_map_value::type +typename property_map_value::type location_to_point(const typename internal::Locate_types::Face_location& loc, const TriangleMesh& tm) { @@ -1069,12 +1117,12 @@ locate_in_face(const typename boost::graph_traits::halfedge_descri /// \cgalParamBegin{vertex_point_map} /// the property map with the points associated to the vertices of `tm`. /// If this parameter is omitted, an internal property map for -/// `CGAL::vertex_point_t` should be available in `TriangleMesh`. +/// `boost::vertex_point_t` should be available in `TriangleMesh`. /// \cgalParamEnd /// \cgalNamedParamsEnd /// template -typename internal::Locate_types::Face_location +typename internal::Locate_types::Face_location locate_in_face(const typename internal::Locate_types::Point& query, const typename boost::graph_traits::face_descriptor fd, const TriangleMesh& tm, @@ -1089,11 +1137,8 @@ locate_in_face(const typename internal::Locate_types::const_type VertexPointMap; typedef typename boost::property_traits::value_type Point; - using boost::choose_param; - using boost::get_param; - - VertexPointMap vpm = choose_param(get_param(np, internal_np::vertex_point), - get_const_property_map(boost::vertex_point, tm)); + VertexPointMap vpm = boost::choose_param(boost::get_param(np, internal_np::vertex_point), + get_const_property_map(boost::vertex_point, tm)); vertex_descriptor vd0 = source(halfedge(fd, tm), tm); vertex_descriptor vd1 = target(halfedge(fd, tm), tm); @@ -1103,7 +1148,7 @@ locate_in_face(const typename internal::Locate_types coords = barycentric_coordinates(p0, p1, p2, query); + std::array coords = barycentric_coordinates(p0, p1, p2, query); if(!is_in_face(coords, tm)) { @@ -1120,7 +1165,7 @@ locate_in_face(const typename internal::Locate_types typename internal::Locate_types::Face_location -locate_in_face(const typename property_map_value::type& query, +locate_in_face(const typename property_map_value::type& query, const typename boost::graph_traits::face_descriptor f, const TriangleMesh& tm) { @@ -1247,7 +1292,7 @@ locate_in_common_face(typename internal::Locate_types::Face_locati if(fd == boost::graph_traits::null_face()) continue; - // check if query can be found in that face + // check if 'query' can be found in that face query_location = locate_in_face(query, fd, tm); internal::snap_location_to_border(query_location, tolerance); // @tmp keep or not ? @@ -1367,32 +1412,55 @@ namespace internal { template::value> -struct Point_to_Point_3 +struct Point_to_Point_3 // 2D case { - typedef typename CGAL::Kernel_traits< - typename property_map_value::type>::Kernel K; - typedef typename K::Point_3 Point_3; + typedef typename internal::Locate_types::Kernel K; + typedef typename K::Point_3 Point_3; - Point_3 operator()(const Point& p) const { return Point_3(p.x(), p.y(), 0.); } + Point_3 operator()(const Point& p) const { return Point_3(p.x(), p.y(), 0); } +}; + +template +struct Point_to_Point_3::Kernel::Point_3, + 3> // 3D case with nothing to do +{ + typedef typename internal::Locate_types::Kernel K; + typedef typename K::Point_3 Point_3; + + const Point_3& operator()(const Point_3& p) const { return p; } }; template -struct Point_to_Point_3 +struct Point_to_Point_3 // Generic 3D case { - typedef typename CGAL::Kernel_traits< - typename property_map_value::type>::Kernel K; - typedef typename K::Point_3 Point_3; + typedef typename internal::Locate_types::Kernel K; + typedef typename K::Point_3 Point_3; - const Point_3& operator()(const Point_3& p) const { return p; } - Point_3 operator()(const Point& p) { return Point_3(p.x(), p.y(), p.z()); } + Point_3 operator()(const Point& p) const { return Point_3(p.x(), p.y(), p.z()); } +}; + +template +struct Ray_to_Ray_3 // 2D case +{ + typedef typename internal::Locate_types::Kernel K; + typedef typename K::Ray_2 Ray_2; + typedef typename K::Ray_3 Ray_3; + typedef Point_to_Point_3 P2_to_P3; + + Ray_3 operator()(const Ray_2& r) const + { + P2_to_P3 to_p3; + return Ray_3(to_p3(r.source()), to_p3(r.second_point())); + } + + const Ray_3& operator()(const Ray_3& r) const { return r; } }; /// Readable property map that converts the output of a given vertex point map to a 3D point template::const_type> + boost::vertex_point_t>::const_type> struct Point_to_Point_3_VPM { private: @@ -1408,14 +1476,11 @@ public: typedef typename K::Point_3 Point_3; // required typedefs - typedef vertex_descriptor key_type; + typedef typename boost::property_traits::key_type key_type; typedef Point_3 value_type; typedef value_type reference; typedef boost::readable_property_map_tag category; - CGAL_static_assertion((boost::is_same::key_type, - vertex_descriptor>::value)); - // Constructors Point_to_Point_3_VPM() : conv_(), vpm_() { } // required for compilation by AABBtraits Point_to_Point_3_VPM(const VertexPointMap vpm) : conv_(), vpm_(vpm) { } @@ -1438,6 +1503,38 @@ private: VertexPointMap vpm_; }; +// Two different functions, because the AABB's traits' VPM must match the passed VPM (that is, +// the original VPM wrapped with P_to_P3_VPM if the VPM's value_type was not Point_3) +template +void build_AABB_tree(const TriangleMesh& tm, + AABB_tree& outTree, + const VPM& wrapped_vpm, + typename std::enable_if< + std::is_same< + typename AABBTraits::Point_3, typename boost::property_traits::value_type + >::value>::type* = 0) +{ + typename boost::graph_traits::face_iterator ffirst, fbeyond; + boost::tie(ffirst, fbeyond) = faces(tm); + outTree.rebuild(ffirst, fbeyond, tm, wrapped_vpm); + outTree.build(); +} + +template +void build_AABB_tree(const TriangleMesh& tm, + AABB_tree& outTree, + const VPM& vpm, + typename std::enable_if< + !std::is_same< + typename AABBTraits::Point_3, typename boost::property_traits::value_type + >::value>::type* = 0) +{ + typedef internal::Point_to_Point_3_VPM Wrapped_VPM; + const Wrapped_VPM wrapped_vpm(vpm); + + return internal::build_AABB_tree(tm, outTree, wrapped_vpm); +} + } // namespace internal /// \name Nearest Face Location Queries @@ -1460,7 +1557,7 @@ private: /// \cgalParamBegin{vertex_point_map} /// the property map with the points associated to the vertices of `tm`. /// If this parameter is omitted, an internal property map for -/// `CGAL::vertex_point_t` should be available in `TriangleMesh`. +/// `boost::vertex_point_t` should be available in `TriangleMesh`. /// \cgalParamEnd /// \cgalNamedParamsEnd /// @@ -1469,18 +1566,12 @@ void build_AABB_tree(const TriangleMesh& tm, AABB_tree& outTree, const NamedParameters& np) { - typedef typename GetVertexPointMap::const_type VertexPointMap; + typedef typename GetVertexPointMap::const_type VertexPointMap; - using boost::choose_param; - using boost::get_param; + const VertexPointMap vpm = boost::choose_param(boost::get_param(np, internal_np::vertex_point), + get_const_property_map(boost::vertex_point, tm)); - VertexPointMap vpm = choose_param(get_param(np, internal_np::vertex_point), - get_const_property_map(boost::vertex_point, tm)); - - typename boost::graph_traits::face_iterator facesStart, facesEnd; - boost::tie(facesStart, facesEnd) = faces(tm); - outTree.rebuild(facesStart, facesEnd, tm, vpm); - outTree.build(); + return internal::build_AABB_tree(tm, outTree, vpm); } template @@ -1505,25 +1596,46 @@ void build_AABB_tree(const TriangleMesh& tm, /// \cgalParamBegin{vertex_point_map} /// the property map with the points associated to the vertices of `tm`. /// If this parameter is omitted, an internal property map for -/// `CGAL::vertex_point_t` should be available in `TriangleMesh`. +/// `boost::vertex_point_t` should be available in `TriangleMesh`. /// \cgalParamEnd /// \cgalNamedParamsEnd /// template -typename internal::Locate_types::Face_location -locate_with_AABB_tree(const typename AABBTraits::Point_3& p, +typename internal::Locate_types::Face_location +locate_with_AABB_tree(const typename internal::Locate_types::Point& p, const AABB_tree& tree, const TriangleMesh& tm, const NamedParameters& np) { - typename AABB_tree::Point_and_primitive_id result = tree.closest_point_and_primitive(p); + typedef typename internal::Locate_types::Point Point; + typedef internal::Point_to_Point_3 P_to_P3; + typedef typename AABBTraits::Point_3 Point_3; + CGAL_static_assertion((std::is_same::value)); - return locate_in_face(result.first, result.second, tm, np); + typedef typename GetVertexPointMap::const_type VertexPointMap; + typedef internal::Point_to_Point_3_VPM WrappedVPM; + + const Point_3& p3 = P_to_P3()(p); + typename AABB_tree::Point_and_primitive_id result = tree.closest_point_and_primitive(p3); + + // The VPM might return a point of any dimension, but the AABB tree necessarily returns + // a Point_3. So, wrap the VPM (again) to give a Point_3. Even if it's already wrapped, we're just + // forwarding a const& anyway. + const VertexPointMap vpm = boost::choose_param(boost::get_param(np, internal_np::vertex_point), + get_const_property_map(boost::vertex_point, tm)); + const WrappedVPM wrapped_vpm(vpm); + + std::cout << "Result: " << result.first << std::endl; + std::cout << "in face: " << result.second->vertex(0)->point() << std::endl; + std::cout << "in face: " << result.second->vertex(1)->point() << std::endl; + std::cout << "in face: " << result.second->vertex(2)->point() << std::endl; + + return locate_in_face(result.first, result.second, tm, CGAL::parameters::vertex_point_map(wrapped_vpm)); } template typename internal::Locate_types::Face_location -locate_with_AABB_tree(const typename AABBTraits::Point_3& p, +locate_with_AABB_tree(const typename internal::Locate_types::Point& p, const AABB_tree& tree, const TriangleMesh& tm) { @@ -1549,42 +1661,47 @@ locate_with_AABB_tree(const typename AABBTraits::Point_3& p, /// \cgalParamBegin{vertex_point_map} /// the property map with the points associated to the vertices of `tm`. /// If this parameter is omitted, an internal property map for -/// `CGAL::vertex_point_t` should be available in `TriangleMesh`. +/// `boost::vertex_point_t` should be available in `TriangleMesh`. /// \cgalParamEnd /// \cgalNamedParamsEnd /// template -typename internal::Locate_types::Face_location -locate(const typename internal::Locate_types::Point& p, +typename internal::Locate_types::Face_location +locate(const typename internal::Locate_types::Point& p, const TriangleMesh& tm, const NamedParameters& np) { - typedef typename GetVertexPointMap::const_type VertexPointMap; - // Wrap the input VPM with a one converting to 3D (costs nothing if the input VPM // already has value type Kernel::Point_3) - typedef internal::Point_to_Point_3_VPM VPM; + typedef typename GetVertexPointMap::const_type VertexPointMap; + typedef internal::Point_to_Point_3_VPM WrappedVPM; + typedef typename internal::Locate_types::Point Intrinsic_point; - typedef AABB_face_graph_triangle_primitive AABB_face_graph_primitive; - typedef CGAL::AABB_traits::Kernel, - AABB_face_graph_primitive> AABB_face_graph_traits; + typedef AABB_face_graph_triangle_primitive AABB_face_graph_primitive; + typedef CGAL::AABB_traits< + typename internal::Locate_types::Kernel, + AABB_face_graph_primitive> AABB_face_graph_traits; - using boost::get_param; - using boost::choose_param; + typedef internal::Point_to_Point_3 P_to_P3; + typedef typename AABB_face_graph_traits::Point_3 Point_3; - const VertexPointMap vpm = choose_param(get_param(np, internal_np::vertex_point), - get_const_property_map(boost::vertex_point, tm)); - const VPM wrapped_vpm(vpm); + CGAL_static_assertion((std::is_same::value)); + + const VertexPointMap vpm = boost::choose_param(boost::get_param(np, internal_np::vertex_point), + get_const_property_map(boost::vertex_point, tm)); + const WrappedVPM wrapped_vpm(vpm); AABB_tree tree; build_AABB_tree(tm, tree, parameters::vertex_point_map(wrapped_vpm)); - return locate_with_AABB_tree(p, tree, tm, parameters::vertex_point_map(wrapped_vpm)); + const Point_3& p3 = P_to_P3()(p); + + return locate_with_AABB_tree(p3, tree, tm, parameters::vertex_point_map(wrapped_vpm)); } template typename internal::Locate_types::Face_location -locate(const typename property_map_value::type& p, +locate(const typename property_map_value::type& p, const TriangleMesh& tm) { return locate(p, tm, parameters::all_default()); @@ -1592,7 +1709,7 @@ locate(const typename property_map_value::ty namespace internal { -// The Ray must have the same ambient dimension as the property map's value type (point type) +// The Ray must have the same ambient dimension as the property map's value type (aka, the point type) template ::value> @@ -1626,32 +1743,39 @@ struct Ray_type_selector /// \cgalParamBegin{vertex_point_map} /// the property map with the points associated to the vertices of `tm`. /// If this parameter is omitted, an internal property map for -/// `CGAL::vertex_point_t` should be available in `TriangleMesh`. +/// `boost::vertex_point_t` should be available in `TriangleMesh`. /// \cgalParamEnd /// \cgalNamedParamsEnd /// template typename internal::Locate_types::Face_location -locate_with_AABB_tree(const typename CGAL::Kernel_traits::type::Ray_3& ray, +locate_with_AABB_tree(const typename internal::Ray_type_selector< + typename internal::Locate_types::Point>::type& ray, const AABB_tree& tree, const TriangleMesh& tm, const NamedParameters& np) { - typedef typename CGAL::Kernel_traits::type Kernel; + typedef typename internal::Locate_types::Kernel Kernel; - typedef typename Kernel::FT FT; - typedef typename Kernel::Point_3 Point_3; - typedef typename Kernel::Ray_3 Ray_3; + typedef typename Kernel::FT FT; + typedef typename Kernel::Point_3 Point_3; + typedef typename Kernel::Ray_3 Ray_3; - typedef typename boost::graph_traits::face_descriptor face_descriptor; + typedef typename boost::graph_traits::face_descriptor face_descriptor; - typedef AABB_tree AABB_face_graph_tree; + typedef typename GetVertexPointMap::const_type VertexPointMap; + typedef internal::Point_to_Point_3_VPM WrappedVPM; + typedef internal::Ray_to_Ray_3 R_to_R3; + typedef AABB_tree AABB_face_graph_tree; typedef typename AABB_face_graph_tree::template Intersection_and_primitive_id::Type Intersection_type; - typedef boost::optional Ray_intersection; + typedef boost::optional Ray_intersection; + + // First, transform the ray into a 3D ray if needed + Ray_3 ray_3 = R_to_R3()(ray); std::vector intersections; - tree.all_intersections(ray, std::back_inserter(intersections)); + tree.all_intersections(ray_3, std::back_inserter(intersections)); bool found = false; FT nearest_distance = 0; @@ -1666,7 +1790,7 @@ locate_with_AABB_tree(const typename CGAL::Kernel_traits::null_face(), CGAL::make_array(FT(0), FT(0), FT(0))); @@ -1688,7 +1819,8 @@ locate_with_AABB_tree(const typename CGAL::Kernel_traits typename internal::Locate_types::Face_location -locate_with_AABB_tree(const typename CGAL::Kernel_traits::type::Ray_3& ray, +locate_with_AABB_tree(const typename internal::Ray_type_selector< + typename internal::Locate_types::Point>::type& ray, const AABB_tree& tree, const TriangleMesh& tm) { @@ -1714,12 +1846,12 @@ locate_with_AABB_tree(const typename CGAL::Kernel_traits -typename internal::Locate_types::Face_location +typename internal::Locate_types::Face_location locate(const typename internal::Ray_type_selector< typename internal::Locate_types::Point>::type& ray, const TriangleMesh& tm, @@ -1735,11 +1867,8 @@ locate(const typename internal::Ray_type_selector< typedef CGAL::AABB_traits::Kernel, AABB_face_graph_primitive> AABB_face_graph_traits; - using boost::get_param; - using boost::choose_param; - - const VertexPointMap vpm = choose_param(get_param(np, internal_np::vertex_point), - get_const_property_map(boost::vertex_point, tm)); + const VertexPointMap vpm = boost::choose_param(boost::get_param(np, internal_np::vertex_point), + get_const_property_map(boost::vertex_point, tm)); const VPM wrapped_vpm(vpm); AABB_tree tree;