diff --git a/AABB_tree/doc/AABB_tree/Concepts/AABBRayIntersectionTraits.h b/AABB_tree/doc/AABB_tree/Concepts/AABBRayIntersectionTraits.h index 630a569ebfb..215cbd65924 100644 --- a/AABB_tree/doc/AABB_tree/Concepts/AABBRayIntersectionTraits.h +++ b/AABB_tree/doc/AABB_tree/Concepts/AABBRayIntersectionTraits.h @@ -37,7 +37,7 @@ public: operator()(const Ray_3& r, const Primitive& primitive)`. A common algorithm to compute the intersection between a bounding box and a ray is the + HREF="https://education.siggraph.org/static/HyperGraph/raytrace/rtinter3.htm">the slab method. */ typedef unspecified_type Intersection_distance; diff --git a/Arrangement_on_surface_2/demo/Arrangement_on_surface_2/ArrangementPainterOstream.cpp b/Arrangement_on_surface_2/demo/Arrangement_on_surface_2/ArrangementPainterOstream.cpp index 1fa59ad88de..bdc2c17f627 100644 --- a/Arrangement_on_surface_2/demo/Arrangement_on_surface_2/ArrangementPainterOstream.cpp +++ b/Arrangement_on_surface_2/demo/Arrangement_on_surface_2/ArrangementPainterOstream.cpp @@ -289,7 +289,7 @@ ArrangementPainterOstream>::operator<<( QRectF seg_bb = this->convert(seg.bbox()); if ( this->clippingRect.isValid() && - !this->clippingRect.intersects(seg_bb) & + !this->clippingRect.intersects(seg_bb) && (!seg.is_horizontal() && !seg.is_vertical())) { return *this; } diff --git a/Documentation/doc/Documentation/windows.txt b/Documentation/doc/Documentation/windows.txt index 3943da4b88f..02596bcd35d 100644 --- a/Documentation/doc/Documentation/windows.txt +++ b/Documentation/doc/Documentation/windows.txt @@ -3,7 +3,7 @@ \cgalAutoToc \cgal \cgalReleaseNumber is supported for the following \ms Visual `C++` compilers: -14.0, 15.9, 16.0 (\visualstudio 2015, 2017, and 2019). +14.0, 15.9, 16.0, 17.0 (\visualstudio 2015, 2017, 2019, and 2022). \cgal is a library that has mandatory dependencies that must be first installed: \ref thirdpartyBoost and \ref thirdpartyMPFR. diff --git a/Kernel_23/include/CGAL/Kernel_23/internal/Projection_traits_3.h b/Kernel_23/include/CGAL/Kernel_23/internal/Projection_traits_3.h index 7a8663127cc..120fa052a85 100644 --- a/Kernel_23/include/CGAL/Kernel_23/internal/Projection_traits_3.h +++ b/Kernel_23/include/CGAL/Kernel_23/internal/Projection_traits_3.h @@ -1005,7 +1005,7 @@ public: Equal_x_2 eqx; Equal_y_2 eqy; - return eqx(p,q) & eqy(p,q); + return eqx(p,q) && eqy(p,q); } }; diff --git a/Mesh_3/include/CGAL/Mesh_3/search_for_connected_components_in_labeled_image.h b/Mesh_3/include/CGAL/Mesh_3/search_for_connected_components_in_labeled_image.h index 7a666d7bed2..56c5d7bc6d5 100644 --- a/Mesh_3/include/CGAL/Mesh_3/search_for_connected_components_in_labeled_image.h +++ b/Mesh_3/include/CGAL/Mesh_3/search_for_connected_components_in_labeled_image.h @@ -73,8 +73,7 @@ search_for_connected_components_in_labeled_image(const CGAL::Image_3& image, for(uint i=0; i::digits; - CGAL_assertion(CGAL_NTS is_finite(d) & is_valid(d)); + CGAL_assertion(CGAL_NTS is_finite(d) && is_valid(d)); int exp; double x = std::frexp(d, &exp); // x in [1/2, 1], x*2^exp = d mpz_init_set_d (man(), // to the following integer: diff --git a/Number_types/include/CGAL/MP_Float.h b/Number_types/include/CGAL/MP_Float.h index 68cfdecde44..030367c4ee9 100644 --- a/Number_types/include/CGAL/MP_Float.h +++ b/Number_types/include/CGAL/MP_Float.h @@ -762,7 +762,7 @@ namespace INTERN_MP_FLOAT { while (true) { x = x % y; if (x == 0) { - CGAL_postcondition(internal::divides(y, a) & internal::divides(y, b)); + CGAL_postcondition(internal::divides(y, a) && internal::divides(y, b)); y.gcd_normalize(); return y; } diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/stitch_borders.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/stitch_borders.h index 49c3fa8b49a..d729cc2654a 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/stitch_borders.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/stitch_borders.h @@ -58,29 +58,46 @@ namespace internal { ////// Helper structs // Used to compare halfedges based on their geometry -template +template struct Less_for_halfedge { + typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; typedef typename boost::graph_traits::halfedge_descriptor halfedge_descriptor; typedef typename boost::property_traits::reference Point; - Less_for_halfedge(const PolygonMesh& pmesh_, const VertexPointMap& vpm_) - : pmesh(pmesh_), vpm(vpm_) + Less_for_halfedge(const PolygonMesh& pmesh_, const VertexPointMap vpm_, const GeomTraits& gt_) + : pmesh(pmesh_), vpm(vpm_), gt(gt_) {} bool operator()(const halfedge_descriptor h1, const halfedge_descriptor h2) const { - Point s1 = get(vpm,target(opposite(h1, pmesh), pmesh)); - Point t1 = get(vpm,target(h1, pmesh)); - Point s2 = get(vpm,target(opposite(h2, pmesh), pmesh)); - Point t2 = get(vpm,target(h2, pmesh)); + typename GeomTraits::Equal_3 equal = gt.equal_3_object(); + typename GeomTraits::Less_xyz_3 less = gt.less_xyz_3_object(); - return (s1 < t1 ? std::make_pair(s1,t1) : std::make_pair(t1, s1)) - < (s2 < t2 ? std::make_pair(s2,t2) : std::make_pair(t2, s2)); + vertex_descriptor vm1 = source(h1, pmesh); + vertex_descriptor vM1 = target(h1, pmesh); + vertex_descriptor vm2 = source(h2, pmesh); + vertex_descriptor vM2 = target(h2, pmesh); + + if(less(get(vpm, vM1), get(vpm, vm1))) + std::swap(vM1, vm1); + if(less(get(vpm, vM2), get(vpm, vm2))) + std::swap(vM2, vm2); + + Point pm1 = get(vpm, vm1); + Point pM1 = get(vpm, vM1); + Point pm2 = get(vpm, vm2); + Point pM2 = get(vpm, vM2); + + if(equal(pm1, pm2)) + return less(pM1, pM2); + + return less(pm1, pm2); } const PolygonMesh& pmesh; - const VertexPointMap& vpm; + const VertexPointMap vpm; + const GeomTraits& gt; }; // The following structs determine which of the two halfedges is kept when a pair is merged @@ -316,15 +333,19 @@ template + typename GT> void fill_pairs(const Halfedge& he, Border_halfedge_map& border_halfedge_map, Halfedge_pair& halfedge_pairs, Manifold_halfedge_pair& manifold_halfedge_pairs, + const Mesh& pmesh, VPM vpm, - const Mesh& pmesh) + const GT& gt) { + typename GT::Equal_3 equal = gt.equal_3_object(); + typename Border_halfedge_map::iterator set_it; bool insertion_ok; std::tie(set_it, insertion_ok) = border_halfedge_map.emplace(he, std::make_pair(1,0)); @@ -337,8 +358,8 @@ void fill_pairs(const Halfedge& he, 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(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))) + if(equal(get(vpm, source(he,pmesh)), get(vpm, target(other_he, pmesh))) && + equal(get(vpm, target(he,pmesh)), get(vpm, source(other_he, pmesh)))) { // 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)) @@ -379,6 +400,9 @@ OutputIterator collect_duplicated_stitchable_boundary_edges(const HalfedgeRange& VPM vpm = choose_parameter(get_parameter(np, internal_np::vertex_point), get_const_property_map(vertex_point, pmesh)); + typedef typename GetGeomTraits::type GT; + GT gt = choose_parameter(get_parameter(np, internal_np::geom_traits)); + typedef CGAL::dynamic_face_property_t Face_property_tag; typedef typename boost::property_map::type Face_cc_map; @@ -386,10 +410,10 @@ OutputIterator collect_duplicated_stitchable_boundary_edges(const HalfedgeRange& std::size_t num_cc = 0; std::vector > border_edges_per_cc; - typedef Less_for_halfedge Less_hedge; + typedef Less_for_halfedge Less_hedge; typedef std::map, Less_hedge> Border_halfedge_map; - Less_hedge less_hedge(pmesh, vpm); + Less_hedge less_hedge(pmesh, vpm, gt); Border_halfedge_map border_halfedge_map(less_hedge); std::vector > halfedge_pairs; @@ -420,7 +444,7 @@ OutputIterator collect_duplicated_stitchable_boundary_edges(const HalfedgeRange& if(per_cc) border_edges_per_cc[get(cc, face(opposite(he, pmesh), pmesh))].push_back(he); else - fill_pairs(he, border_halfedge_map, halfedge_pairs, manifold_halfedge_pairs, vpm, pmesh); + fill_pairs(he, border_halfedge_map, halfedge_pairs, manifold_halfedge_pairs, pmesh, vpm, gt); } if(per_cc) @@ -435,7 +459,7 @@ OutputIterator collect_duplicated_stitchable_boundary_edges(const HalfedgeRange& { halfedge_descriptor he = border_edges_per_cc[i][j]; fill_pairs(he, border_halfedge_map_in_cc, halfedge_pairs, - manifold_halfedge_pairs, vpm, pmesh); + manifold_halfedge_pairs, pmesh, vpm, gt); } // put in `out` only manifold edges from the set of edges to stitch. @@ -863,17 +887,21 @@ std::size_t stitch_halfedge_range_dispatcher(const HalfedgePairRange& to_stitch_ // However, even if non-manifoldness exists within a loop, it is safe choice to stitch consecutive // stitchable halfedges template + typename GT> std::size_t zip_boundary_cycle(typename boost::graph_traits::halfedge_descriptor& bh, const HalfedgeRange& cycle_halfedges, + const HalfedgeKeeper& hd_kpr, PolygonMesh& pmesh, const VPM vpm, - const HalfedgeKeeper& hd_kpr) + const GT& gt) { typedef typename boost::graph_traits::halfedge_descriptor halfedge_descriptor; + typename GT::Equal_3 equal = gt.equal_3_object(); + std::size_t stitched_boundary_cycles_n = 0; // Zipping cannot change the topology of the hole so the maintenance is trivial @@ -909,10 +937,11 @@ std::size_t zip_boundary_cycle(typename boost::graph_traits::halfed do { halfedge_descriptor hnn = next(hn, pmesh); - CGAL_assertion(get(vpm, target(hn, pmesh)) == get(vpm, source(hnn, pmesh))); + CGAL_assertion(equal(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(equal(get(vpm, source(hn, pmesh)), get(vpm, target(hnn, pmesh))) && + !is_degenerate_edge(edge(hn, pmesh), pmesh, + parameters::vertex_point_map(vpm).geom_traits(gt))) { if(unstitchable_halfedges.count(hn) == 0) { @@ -978,8 +1007,9 @@ std::size_t zip_boundary_cycle(typename boost::graph_traits::halfed 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))) + if(!equal(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).geom_traits(gt))) { bh = curr_hn; break; @@ -1044,6 +1074,9 @@ std::size_t stitch_boundary_cycle(const typename boost::graph_traits::type GT; + GT gt = choose_parameter(get_parameter(np, internal_np::geom_traits)); + typedef typename internal_np::Lookup_named_param_def >::type Halfedge_keeper; @@ -1056,7 +1089,7 @@ std::size_t stitch_boundary_cycle(const typename boost::graph_traits::null_halfedge()) // stitched everything { cycle_reps_maintainer.remove_representative(bh); @@ -1111,6 +1144,17 @@ std::size_t stitch_boundary_cycle(const typename boost::graph_traits::%vertex_descriptor` -/// as key type and `%Point_3` as value type} -/// \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 -/// /// \cgalParamNBegin{apply_per_connected_component} /// \cgalParamDescription{specifies if the borders should only be stitched only within their own connected component.} /// \cgalParamType{Boolean} @@ -1433,6 +1479,26 @@ std::size_t stitch_borders(PolygonMesh& pmesh, /// as key type and `std::size_t` as value type} /// \cgalParamDefault{an automatically indexed internal map} /// \cgalParamNEnd +/// +/// \cgalParamNBegin{vertex_point_map} +/// \cgalParamDescription{a property map associating points to the vertices of `pmesh`} +/// \cgalParamType{a class model of `ReadWritePropertyMap` with `boost::graph_traits::%vertex_descriptor` +/// as key type and `%Point_3` as value type} +/// \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 +/// +/// \cgalParamNBegin{geom_traits} +/// \cgalParamDescription{an instance of a geometric traits class} +/// \cgalParamType{The traits class must provide the nested type `Point_3`, +/// and the nested functors: +/// - `Less_xyz_3` to compare lexicographically two points +/// - `Equal_3` to check whether two points are identical. +/// For each functor `Foo`, a function `Foo foo_object()` must be provided.} +/// \cgalParamDefault{a \cgal Kernel deduced from the point type, using `CGAL::Kernel_traits`} +/// \cgalParamExtra{The geometric traits class must be compatible with the vertex point type.} +/// \cgalParamNEnd /// \cgalNamedParamsEnd /// /// \return the number of pairs of halfedges that were stitched. diff --git a/STL_Extension/doc/STL_Extension/Concepts/Hashable.h b/STL_Extension/doc/STL_Extension/Concepts/Hashable.h index 17fa2304230..3c7f0b5cf71 100644 --- a/STL_Extension/doc/STL_Extension/Concepts/Hashable.h +++ b/STL_Extension/doc/STL_Extension/Concepts/Hashable.h @@ -12,8 +12,8 @@ of the `graph_traits` header files provided by \cgal. They can be disables by defining the macro `CGAL_DISABLE_HASH_OPENMESH`. \sa `CGAL::Unique_hash_map` -\sa `std::unordered_set` -\sa `std::unordered_map` +\sa `std::unordered_set` +\sa `std::unordered_map` \sa `boost::unordered_set` \sa `boost::unordered_map` diff --git a/STL_Extension/include/CGAL/Uncertain.h b/STL_Extension/include/CGAL/Uncertain.h index 00a7eaaa82f..5e5003a5d5b 100644 --- a/STL_Extension/include/CGAL/Uncertain.h +++ b/STL_Extension/include/CGAL/Uncertain.h @@ -289,6 +289,11 @@ Uncertain operator!(Uncertain a) return Uncertain(!a.sup(), !a.inf()); } +#ifdef __clang__ +# pragma GCC diagnostic push +# pragma GCC diagnostic ignored "-Wunknown-warning-option" +# pragma GCC diagnostic ignored "-Wbitwise-instead-of-logical" +#endif inline Uncertain operator|(Uncertain a, Uncertain b) { @@ -324,7 +329,9 @@ Uncertain operator&(Uncertain a, bool b) { return Uncertain(a.inf() & b, a.sup() & b); } - +#ifdef __clang__ +# pragma GCC diagnostic pop +#endif // Equality operators diff --git a/Straight_skeleton_2/include/CGAL/constructions/Straight_skeleton_cons_ftC2.h b/Straight_skeleton_2/include/CGAL/constructions/Straight_skeleton_cons_ftC2.h index 4f872f6a779..a391e44029a 100644 --- a/Straight_skeleton_2/include/CGAL/constructions/Straight_skeleton_cons_ftC2.h +++ b/Straight_skeleton_2/include/CGAL/constructions/Straight_skeleton_cons_ftC2.h @@ -45,7 +45,7 @@ bool are_parallel_edges_equally_oriented( Segment_2_with_ID const& e0, Segmen template bool are_edges_orderly_collinear( Segment_2_with_ID const& e0, Segment_2_with_ID const& e1 ) { - return are_edges_collinear(e0,e1) & are_parallel_edges_equally_oriented(e0,e1); + return are_edges_collinear(e0,e1) && are_parallel_edges_equally_oriented(e0,e1); } diff --git a/Surface_mesh_shortest_path/examples/Surface_mesh_shortest_path/shortest_path_sequence.cpp b/Surface_mesh_shortest_path/examples/Surface_mesh_shortest_path/shortest_path_sequence.cpp index adb49aa7130..98340eb0000 100644 --- a/Surface_mesh_shortest_path/examples/Surface_mesh_shortest_path/shortest_path_sequence.cpp +++ b/Surface_mesh_shortest_path/examples/Surface_mesh_shortest_path/shortest_path_sequence.cpp @@ -13,9 +13,11 @@ typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel; typedef CGAL::Surface_mesh Triangle_mesh; + typedef CGAL::Surface_mesh_shortest_path_traits Traits; typedef CGAL::Surface_mesh_shortest_path Surface_mesh_shortest_path; typedef Traits::Barycentric_coordinates Barycentric_coordinates; + typedef boost::graph_traits Graph_traits; typedef Graph_traits::vertex_iterator vertex_iterator; typedef Graph_traits::face_iterator face_iterator; @@ -34,23 +36,24 @@ struct Sequence_collector void operator()(halfedge_descriptor he, double alpha) { - - sequence.push_back( std::make_pair(he, alpha) ); + sequence.push_back(std::make_pair(he, alpha)); } void operator()(vertex_descriptor v) { - sequence.push_back( v ); + sequence.push_back(v); } void operator()(face_descriptor f, Barycentric_coordinates alpha) { - sequence.push_back( std::make_pair(f, alpha) ); + sequence.push_back(std::make_pair(f, alpha)); } }; // A visitor to print what a variant contains using boost::apply_visitor -struct Print_visitor : public boost::static_visitor<> { +struct Print_visitor + : public boost::static_visitor<> +{ int i; Triangle_mesh& g; @@ -58,22 +61,25 @@ struct Print_visitor : public boost::static_visitor<> { void operator()(vertex_descriptor v) { - std::cout << "#" << ++i << " : Vertex : " << get(boost::vertex_index, g)[v] << "\n"; + std::cout << "#" << ++i << " Vertex: " << get(boost::vertex_index, g)[v]; + std::cout << " Position: " << Surface_mesh_shortest_path::point(v, g) << "\n"; } void operator()(const std::pair& h_a) { - std::cout << "#" << ++i << " : Edge : " << get(CGAL::halfedge_index, g)[h_a.first] << " , (" - << 1.0 - h_a.second << " , " - << h_a.second << ")\n"; + std::cout << "#" << ++i << " Edge: " << get(CGAL::halfedge_index, g)[h_a.first] << " , (" + << 1.0 - h_a.second << " , " + << h_a.second << ")"; + std::cout << " Position: " << Surface_mesh_shortest_path::point(h_a.first, h_a.second, g) << "\n"; } void operator()(const std::pair& f_bc) { - std::cout << "#" << ++i << " : Face : " << get(CGAL::face_index, g)[f_bc.first] << " , (" - << f_bc.second[0] << " , " - << f_bc.second[1] << " , " - << f_bc.second[2] << ")\n"; + std::cout << "#" << ++i << " Face: " << get(CGAL::face_index, g)[f_bc.first] << " , (" + << f_bc.second[0] << " , " + << f_bc.second[1] << " , " + << f_bc.second[2] << ")"; + std::cout << " Position: " << Surface_mesh_shortest_path::point(f_bc.first, f_bc.second, g) << "\n"; } }; @@ -100,12 +106,16 @@ int main(int argc, char** argv) // construct a shortest path query object and add a source point Surface_mesh_shortest_path shortest_paths(tmesh); + + std::cout << "Add source: " << Surface_mesh_shortest_path::point(*face_it, face_location, tmesh) << std::endl; shortest_paths.add_source_point(*face_it, face_location); // pick a random target point inside a face face_it = faces(tmesh).first; std::advance(face_it, rand.get_int(0, static_cast(num_faces(tmesh)))); + std::cout << "Target is: " << Surface_mesh_shortest_path::point(*face_it, face_location, tmesh) << std::endl; + // collect the sequence of simplicies crossed by the shortest path Sequence_collector sequence_collector; shortest_paths.shortest_path_sequence_to_source_points(*face_it, face_location, sequence_collector); diff --git a/Surface_mesh_shortest_path/include/CGAL/Surface_mesh_shortest_path/Surface_mesh_shortest_path.h b/Surface_mesh_shortest_path/include/CGAL/Surface_mesh_shortest_path/Surface_mesh_shortest_path.h index 270f003af03..f116f933c2e 100644 --- a/Surface_mesh_shortest_path/include/CGAL/Surface_mesh_shortest_path/Surface_mesh_shortest_path.h +++ b/Surface_mesh_shortest_path/include/CGAL/Surface_mesh_shortest_path/Surface_mesh_shortest_path.h @@ -14,36 +14,29 @@ #include -#include - -#include -#include -#include -#include -#include -#include -#include - -#include -#include -#include -#include -#include - -#include -#include -#include -#include -#include - #include #include #include +#include +#include #include #include +#include +#include +#include + +#include #include +#include +#include +#include +#include +#include +#include +#include + namespace CGAL { /*! @@ -276,8 +269,6 @@ public: /// @} private: - typedef typename Graph_traits::vertex_iterator vertex_iterator; - typedef typename Graph_traits::halfedge_iterator halfedge_iterator; typedef typename Graph_traits::face_iterator face_iterator; typedef typename Traits::Triangle_3 Triangle_3; @@ -354,13 +345,11 @@ private: Expansion_priqueue m_expansionPriqueue; #if !defined(NDEBUG) - std::size_t m_currentNodeCount; std::size_t m_peakNodeCount; std::size_t m_queueAtPeakNodes; std::size_t m_peakQueueSize; std::size_t m_nodesAtPeakQueue; - #endif #if !defined(NDEBUG) @@ -420,7 +409,6 @@ public: #endif - public: /// \cond @@ -533,7 +521,8 @@ private: } /* - Filtering algorithm described in Xin and Wang (2009) "Improving chen and han's algorithm on the discrete geodesic problem." + Filtering algorithm described in Xin and Wang (2009) + "Improving chen and han's algorithm on the discrete geodesic problem." https://dl.acm.org/citation.cfm?doid=1559755.1559761 */ bool window_distance_filter(Cone_tree_node* cone, @@ -659,7 +648,6 @@ private: } CGAL_assertion(cone->m_pendingLeftSubtree != nullptr); - cone->m_pendingLeftSubtree = nullptr; if (window_distance_filter(cone, windowSegment, false)) @@ -689,14 +677,13 @@ private: typename Traits::Construct_vertex_2 cv2(m_traits.construct_vertex_2_object()); typename Traits::Construct_triangle_3_along_segment_2_flattening ft3as2(m_traits.construct_triangle_3_along_segment_2_flattening_object()); - CGAL_assertion(cone->m_pendingRightSubtree != nullptr); - cone->m_pendingRightSubtree = nullptr; - if (m_debugOutput) { std::cout << std::endl << " >>>>>>>>>>>>>>>>>>> Expanding RIGHT CHILD <<<<<<<<<<<<<<<<<<<" <m_pendingRightSubtree != nullptr); + cone->m_pendingRightSubtree = nullptr; if (window_distance_filter(cone, windowSegment, true)) { @@ -729,7 +716,7 @@ private: std::size_t associatedEdge; CGAL::Surface_mesh_shortest_paths_3::Barycentric_coordinates_type type; - boost::tie(type, associatedEdge) = classify_barycentric_coordinates(location); + std::tie(type, associatedEdge) = classify_barycentric_coordinates(location); switch (type) { @@ -777,11 +764,13 @@ private: Cone_tree_node* faceRoot = new Cone_tree_node(m_traits, m_graph, m_rootNodes.size()); node_created(); - m_rootNodes.push_back(std::make_pair(faceRoot, sourcePointIt)); + m_rootNodes.emplace_back(faceRoot, sourcePointIt); if (m_debugOutput) { typename Traits::Construct_barycentric_coordinates_weight cbcw(m_traits.construct_barycentric_coordinates_weight_object()); + + std::cout << std::endl << " ~~~~~~~~~~~~~~~~~~~~~~~~~~~" << std::endl; std::cout << "\tFace Root Expansion: id = " << get(m_faceIndexMap, f) << " , Location = " << cbcw(faceLocation, 0) << " " << cbcw(faceLocation, 1) << " " << cbcw(faceLocation, 2) << " " << std::endl; @@ -794,8 +783,12 @@ private: const Barycentric_coordinates rotatedFaceLocation(shifted_coordinates(faceLocation, currentVertex)); const Point_2 sourcePoint(construct_barycenter_in_triangle_2(layoutFace, rotatedFaceLocation)); - Cone_tree_node* child = new Cone_tree_node(m_traits, m_graph, current, layoutFace, sourcePoint, - FT(0), cv2(layoutFace, 0), cv2(layoutFace, 2), + Cone_tree_node* child = new Cone_tree_node(m_traits, m_graph, + current /*entryEdge*/, + layoutFace, sourcePoint, + FT(0) /*pseudoSourceDistance*/, + cv2(layoutFace, 0) /*windowLeft*/, + cv2(layoutFace, 2) /*windowRight*/, Cone_tree_node::FACE_SOURCE); node_created(); faceRoot->push_middle_child(child); @@ -803,7 +796,8 @@ private: if (m_debugOutput) { std::cout << "\tExpanding face root #" << currentVertex << " : " << std::endl;; - std::cout << "\t\tFace = " << layoutFace << std::endl; + std::cout << "\t\t3D Face = " << face3d << std::endl; + std::cout << "\t\t2D Face = " << layoutFace << std::endl; std::cout << "\t\tLocation = " << sourcePoint << std::endl; } @@ -820,53 +814,127 @@ private: const FT t0, const FT t1, Source_point_iterator sourcePointIt) { - typename Traits::Construct_barycenter_2 cb2(m_traits.construct_barycenter_2_object()); + CGAL_precondition(!is_border(baseEdge, m_graph)); + CGAL_precondition(t0 + t1 == FT(1)); + typename Traits::Construct_vertex_2 cv2(m_traits.construct_vertex_2_object()); typename Traits::Construct_triangle_3_to_triangle_2_projection pt3t2(m_traits.construct_triangle_3_to_triangle_2_projection_object()); if (m_debugOutput) { + std::cout << std::endl << " ~~~~~~~~~~~~~~~~~~~~~~~~~~~" << std::endl; std::cout << "\tEdge Root Expansion: faceA = " << get(m_faceIndexMap, face(baseEdge, m_graph)) << " , faceB = " << get(m_faceIndexMap, face(opposite(baseEdge, m_graph), m_graph)) << " , t0 = " << t0 << " , t1 = " << t1 << std::endl; + std::cout << "\t\tBoundary: " << is_border_edge(baseEdge, m_graph) << std::endl; } + Cone_tree_node* edgeRoot = new Cone_tree_node(m_traits, m_graph, m_rootNodes.size()); + node_created(); + m_rootNodes.emplace_back(edgeRoot, sourcePointIt); + + /* If v0v1 is not a border edge: + * + * v2 + * / \ + * / \ + * / \ + * v0 - S - v1 + * \ / + * \ / + * \ / + * v3 + * The source S must reach all Vi, so for each side of the edge, there are two windwows being spawned: + * - v0v1 targetting v2 propagating only on the left (v0v2) + * - v2v0 targetting v1 propagating only on the left (v2v1) + * - v1v0 targetting v3 propagating only on the left (v1v3) + * - v3v1 targetting v0 propagating only on the left (v3v0) + * + * If v0v1 is a border edge, spawn 3 children in the face, and none on the other side + */ + + if(is_border_edge(baseEdge, m_graph)) + { + const Face_location edgeSourceLocation = face_location(baseEdge, t0); + return expand_face_root(face(baseEdge, m_graph), edgeSourceLocation.second, sourcePointIt); + } + + // From here on, it is not a border edge --> spawn 2 children on each side + halfedge_descriptor baseEdges[2]; baseEdges[0] = baseEdge; baseEdges[1] = opposite(baseEdge, m_graph); - Triangle_3 faces3d[2]; - Triangle_2 layoutFaces[2]; - - for (std::size_t i = 0; i < 2; ++i) - { - faces3d[i] = triangle_from_halfedge(baseEdges[i]); - layoutFaces[i] = pt3t2(faces3d[i]); - } - - Point_2 sourcePoints[2]; - sourcePoints[0] = cb2(cv2(layoutFaces[0], 0), t0, cv2(layoutFaces[0], 1), t1); - sourcePoints[1] = cb2(cv2(layoutFaces[1], 0), t1, cv2(layoutFaces[1], 1), t0); - - Cone_tree_node* edgeRoot = new Cone_tree_node(m_traits, m_graph, m_rootNodes.size()); - node_created(); - m_rootNodes.push_back(std::make_pair(edgeRoot, sourcePointIt)); + // shift is because the entry halfedge is not necessarily equal to halfedge(face(entry_h, g), g) + Barycentric_coordinates edgeSourceLocations[2]; + edgeSourceLocations[0] = shifted_coordinates(face_location(baseEdges[0], t0).second, + Surface_mesh_shortest_paths_3::internal::edge_index(baseEdges[0], m_graph)); + edgeSourceLocations[1] = shifted_coordinates(face_location(baseEdges[1], t1).second, + Surface_mesh_shortest_paths_3::internal::edge_index(baseEdges[1], m_graph)); for (std::size_t side = 0; side < 2; ++side) { + Triangle_3 face3d(triangle_from_halfedge(baseEdges[side])); + Triangle_2 layoutFace(pt3t2(face3d)); + Point_2 sourcePoint(construct_barycenter_in_triangle_2(layoutFace, edgeSourceLocations[side])); + + // v0v1 targetting v2 if (m_debugOutput) { - std::cout << "\tExpanding edge root #" << side << " : " << std::endl;; - std::cout << "\t\tFace = " << layoutFaces[side] << std::endl; - std::cout << "\t\tLocation = " << sourcePoints[side] << std::endl; + std::cout << std::endl << " ~~~~~~~~~~~~~~~~~~~~~~~~~~~" << std::endl; + std::cout << "\tExpanding edge root, side #" << side << ", targetting LOCAL 'v2'" << std::endl; + std::cout << "\t\t3D Face = " << face3d << std::endl; + std::cout << "\t\t2D Face = " << layoutFace << std::endl; + std::cout << "\t\tBarycentric coordinates: " << edgeSourceLocations[side][0] + << " " << edgeSourceLocations[side][1] + << " " << edgeSourceLocations[side][2] << std::endl; + std::cout << "\t\tLocation = " << sourcePoint << std::endl; } - Cone_tree_node* mainChild = new Cone_tree_node(m_traits, m_graph, baseEdges[side], layoutFaces[side], - sourcePoints[side], FT(0), cv2(layoutFaces[side], 0), - cv2(layoutFaces[side], 1), Cone_tree_node::EDGE_SOURCE); + Cone_tree_node* v2_Child = new Cone_tree_node(m_traits, m_graph, + baseEdges[side] /*entryEdge*/, + layoutFace, + sourcePoint /*sourceImage*/, + FT(0) /*pseudoSourceDistance*/, + cv2(layoutFace, 0) /*windowLeft*/, + cv2(layoutFace, 2) /*windowRight*/, + Cone_tree_node::EDGE_SOURCE); node_created(); - edgeRoot->push_middle_child(mainChild); - process_node(mainChild); + edgeRoot->push_middle_child(v2_Child); + process_node(v2_Child); + + // v2v0 targetting v1 + face3d = triangle_from_halfedge(prev(baseEdges[side], m_graph)); + layoutFace = pt3t2(face3d); + + // shift the barycentric coordinates to correspond to the new layout + std::swap(edgeSourceLocations[side][1], edgeSourceLocations[side][2]); + std::swap(edgeSourceLocations[side][0], edgeSourceLocations[side][1]); + sourcePoint = Point_2(construct_barycenter_in_triangle_2(layoutFace, edgeSourceLocations[side])); + + if (m_debugOutput) + { + std::cout << std::endl << " ~~~~~~~~~~~~~~~~~~~~~~~~~~~" << std::endl; + std::cout << "\tExpanding edge root, side #" << side << ", targetting LOCAL 'v1'" << std::endl; + std::cout << "\t\t3D Face = " << face3d << std::endl; + std::cout << "\t\t2D Face = " << layoutFace << std::endl; + std::cout << "\t\tBarycentric coordinates: " << edgeSourceLocations[side][0] + << " " << edgeSourceLocations[side][1] + << " " << edgeSourceLocations[side][2] << std::endl; + std::cout << "\t\tLocation = " << sourcePoint << std::endl; + } + + Cone_tree_node* v1_Child = new Cone_tree_node(m_traits, m_graph, + prev(baseEdges[side], m_graph) /*entryEdge*/, + layoutFace, + sourcePoint /*sourceImage*/, + FT(0) /*pseudoSourceDistance*/, + cv2(layoutFace, 0) /*windowLeft*/, + cv2(layoutFace, 2) /*windowRight*/, + Cone_tree_node::EDGE_SOURCE); + node_created(); + edgeRoot->push_middle_child(v1_Child); + process_node(v1_Child); } } @@ -878,6 +946,7 @@ private: { if (m_debugOutput) { + std::cout << std::endl << " ~~~~~~~~~~~~~~~~~~~~~~~~~~~" << std::endl; std::cout << "\tVertex Root Expansion: Vertex = " << get(m_vertexIndexMap, vertex) << std::endl; } @@ -885,7 +954,7 @@ private: prev(halfedge(vertex, m_graph), m_graph)); node_created(); - m_rootNodes.push_back(std::make_pair(vertexRoot, sourcePointIt)); + m_rootNodes.emplace_back(vertexRoot, sourcePointIt); m_closestToVertices[get(m_vertexIndexMap, vertex)] = Node_distance_pair(vertexRoot, FT(0)); @@ -894,6 +963,13 @@ private: /* Create child nodes for each face surrounding the vertex occupied by `parent`, and push them to the queue + + By convention, source windows are left & middle (no right) as to not create overlaps. + A child is also created for the border halfedge (if any) as to propagate distance + to the next vertex on the border (aka `target(next(border_h, g), g)`). + This creates a nonsensical triangle made of `source(border_h, g)`, `target(border_h, g)`, + and `target(next(border_h, g), g)` but propagation is only done to the vertex (vertex source: + no propagation on the right by convention, and left is a border halfedge so no propagation either). */ void expand_pseudo_source(Cone_tree_node* parent) { @@ -911,8 +987,8 @@ private: if (m_debugOutput) { - std::cout << "expansionVertex: " << get(m_vertexIndexMap, expansionVertex) << std::endl; - std::cout << "Distance from target to root: " << distanceFromTargetToRoot << std::endl; + std::cout << "Pseudo source: V" << get(m_vertexIndexMap, expansionVertex) << std::endl; + std::cout << "Distance from pseudo source to root: " << distanceFromTargetToRoot << std::endl; } // A potential optimization could be made by only expanding in the 'necessary' range (i.e. the range outside of geodesic visibility), but the @@ -931,7 +1007,7 @@ private: << get(m_vertexIndexMap, source(currentEdge, m_graph)) << " " << get(m_vertexIndexMap, target(currentEdge, m_graph)) << std::endl; std::cout << "face id = "; - if (face(currentEdge, m_graph) != Graph_traits::null_face()) + if (!is_border(currentEdge, m_graph)) { std::cout << get(m_faceIndexMap, face(currentEdge, m_graph)) << std::endl; @@ -945,9 +1021,11 @@ private: } } - Cone_tree_node* child = new Cone_tree_node(m_traits, m_graph, currentEdge, layoutFace, - cv2(layoutFace, 1), distanceFromTargetToRoot, - cv2(layoutFace, 0), cv2(layoutFace, 2), + Cone_tree_node* child = new Cone_tree_node(m_traits, m_graph, currentEdge /*entryEdge*/, + layoutFace, cv2(layoutFace, 1) /*sourceImage*/, + distanceFromTargetToRoot, + cv2(layoutFace, 0) /*windowLeft*/, + cv2(layoutFace, 2) /*windowRight*/, Cone_tree_node::VERTEX_SOURCE); node_created(); @@ -1131,18 +1209,18 @@ private: std::cout << "\tParent node: " << node->parent() << std::endl; std::cout << "\tParent node type: " << node->parent()->node_type() << std::endl; std::cout << "\tFace = " << node->layout_face() << std::endl; - std::cout << "\tVertices = "; + std::cout << "\tVertices ="; halfedge_descriptor current = node->entry_edge(); for (std::size_t i = 0; i<3; ++i) { - std::cout << get(m_vertexIndexMap, source(current, m_graph)) << " "; + std::cout << " " << get(m_vertexIndexMap, source(current, m_graph)); current = next(current, m_graph); } std::cout << std::endl; std::cout << "\tSource Image = " << node->source_image() << std::endl; - std::cout << "\tEntry Halfedge = (" << get(m_vertexIndexMap, source(node->entry_edge(), m_graph)) << " " + std::cout << "\tEntry Halfedge = (V" << get(m_vertexIndexMap, source(node->entry_edge(), m_graph)) << " V" << get(m_vertexIndexMap, target(node->entry_edge(), m_graph)) << ")" << std::endl; - std::cout << "\tTarget vertex = " << get(m_vertexIndexMap, node->target_vertex()) << std::endl; + std::cout << "\tTarget vertex = V" << get(m_vertexIndexMap, node->target_vertex()) << std::endl; std::cout << "\tWindow Left = " << node->window_left() << std::endl; std::cout << "\tWindow Right = " << node->window_right() << std::endl; @@ -1156,18 +1234,10 @@ private: leftSide = node->has_left_side(); rightSide = node->has_right_side(); } - else // node corresponds to a source + else // source nodes only have left sides { - if (node->node_type() == Cone_tree_node::EDGE_SOURCE) - { - leftSide = true; - rightSide = true; - } - else - { - leftSide = true; - rightSide = false; - } + leftSide = true; + rightSide = false; } if (m_debugOutput) @@ -1189,7 +1259,7 @@ private: std::size_t entryHalfEdgeIndex = get(m_halfedgeIndexMap, node->entry_edge()); - Node_distance_pair currentOccupier = m_vertexOccupiers[entryHalfEdgeIndex]; + const Node_distance_pair& currentOccupier = m_vertexOccupiers[entryHalfEdgeIndex]; FT currentNodeDistance = node->distance_from_target_to_root(); if (m_debugOutput) @@ -1231,11 +1301,11 @@ private: if (m_debugOutput) { - std::cout << "\t Current occupier, EH (" - << get(m_vertexIndexMap, source(currentOccupier.first->entry_edge(), m_graph)) << " " + std::cout << "\t Current occupier, EH (V" + << get(m_vertexIndexMap, source(currentOccupier.first->entry_edge(), m_graph)) << " V" << get(m_vertexIndexMap, target(currentOccupier.first->entry_edge(), m_graph)) << ")" << std::endl; std::cout << "\t Current occupier, Source = " << currentOccupier.first->source_image() << std::endl; - std::cout << "\t Current Occupier Distance = " << currentOccupier.second << std::endl; + std::cout << "\t Current occupier Distance = " << currentOccupier.second << std::endl; std::cout << "\t smaller (-1)/equal (0)/larger (1) comparison? " << c << std::endl; } } @@ -1271,7 +1341,7 @@ private: // This is a consequence of using the same basic node type for source and interval nodes // If this is a source node, it is only pointing to one of the two opposite edges (the left one by convention) - if (node->node_type() != Cone_tree_node::INTERVAL && node->node_type() != Cone_tree_node::EDGE_SOURCE) + if (node->is_source_node()) { propagateRight = false; @@ -1313,9 +1383,9 @@ private: } } - // Check if node is now the absolute closest node, and replace the current closest as appropriate + // Check if `node` is now the absolute closest node, and replace the current closest as appropriate std::size_t targetVertexIndex = get(m_vertexIndexMap, node->target_vertex()); - Node_distance_pair currentClosest = m_closestToVertices[targetVertexIndex]; + const Node_distance_pair& currentClosest = m_closestToVertices[targetVertexIndex]; if (m_debugOutput && currentClosest.first != nullptr) { @@ -1324,12 +1394,13 @@ private: // If equal times, give priority to vertex sources since it's cleaner and simpler to handle than interval windows if (currentClosest.first == nullptr || - currentClosest.second > currentNodeDistance || - (currentClosest.second == currentNodeDistance && node->node_type() == Cone_tree_node::VERTEX_SOURCE)) + currentClosest.second > currentNodeDistance || + (currentClosest.second == currentNodeDistance && + node->node_type() == Cone_tree_node::VERTEX_SOURCE)) { if (m_debugOutput) { - std::cout << "\t Current node is now the closest at target vertex " + std::cout << "\t Current node is now the closest at target vertex V" << get(m_vertexIndexMap, node->target_vertex()) << std::endl; } @@ -1338,7 +1409,7 @@ private: { if (m_debugOutput) { - std::cout << "\t Vertex " << targetVertexIndex << " is a pseudo-source" << std::endl; + std::cout << "\t Vertex V" << targetVertexIndex << " is a pseudo-source" << std::endl; } if (currentClosest.first != nullptr) @@ -1399,11 +1470,13 @@ private: { if (propagateLeft) { + CGAL_assertion(!node->is_null_face()); push_left_child(node); } - if (propagateRight && (!node->is_source_node() || node->node_type() == Cone_tree_node::EDGE_SOURCE)) + if (propagateRight) { + CGAL_assertion(!node->is_source_node()); push_right_child(node); } @@ -1421,9 +1494,17 @@ private: void push_left_child(Cone_tree_node* parent) { + if (m_debugOutput) + { + std::cout << "Tentative push of left child edge " + << " (V" << get(m_vertexIndexMap, source(parent->left_child_edge(), m_graph)) + << " V" << get(m_vertexIndexMap, target(parent->left_child_edge(), m_graph)) << ")" << std::endl; + std::cout << "Boundary? " << is_border(parent->left_child_edge(), m_graph) << std::endl; + } + typename Traits::Compute_squared_distance_2 csd2(m_traits.compute_squared_distance_2_object()); - if (face(parent->left_child_edge(), m_graph) != Graph_traits::null_face()) + if (!is_border(parent->left_child_edge(), m_graph)) { Segment_2 leftWindow; @@ -1467,12 +1548,15 @@ private: { if (m_debugOutput) { - std::cout << "Tentative push of right child..." << std::endl; + std::cout << "Tentative push of right child edge" + << " (V" << get(m_vertexIndexMap, source(parent->right_child_edge(), m_graph)) + << " V" << get(m_vertexIndexMap, target(parent->right_child_edge(), m_graph)) << ")" << std::endl; + std::cout << "Boundary? " << is_border(parent->right_child_edge(), m_graph) << std::endl; } typename Traits::Compute_squared_distance_2 csd2(m_traits.compute_squared_distance_2_object()); - if (face(parent->right_child_edge(), m_graph) != Graph_traits::null_face()) + if (!is_border(parent->right_child_edge(), m_graph)) { Segment_2 rightWindow; bool result = clip_to_bounds(parent->right_child_base_segment(), parent->left_boundary(), @@ -1605,14 +1689,11 @@ private: void set_vertex_types() { - vertex_iterator current, end; - - for (boost::tie(current, end) = vertices(m_graph); current != end; ++current) + for(vertex_descriptor v : vertices(m_graph)) { - if (halfedge(*current, m_graph)==Graph_traits::null_halfedge()) - continue; - std::size_t vertexIndex = get(m_vertexIndexMap, *current); - m_vertexIsPseudoSource[vertexIndex] = (is_saddle_vertex(*current) || is_boundary_vertex(*current)); + std::size_t vertexIndex = get(m_vertexIndexMap, v); + m_vertexIsPseudoSource[vertexIndex] = !internal::is_isolated(v, m_graph) && + (is_saddle_vertex(v) || is_boundary_vertex(v)); } } @@ -1623,21 +1704,7 @@ private: bool is_boundary_vertex(const vertex_descriptor v) const { - halfedge_descriptor h = halfedge(v, m_graph); - halfedge_descriptor first = h; - - do - { - if (face(h, m_graph) == Graph_traits::null_face() || face(opposite(h, m_graph), m_graph) == Graph_traits::null_face()) - { - return true; - } - - h = opposite(next(h, m_graph), m_graph); - } - while(h != first); - - return false; + return bool(is_border(v, m_graph)); } void delete_all_nodes() @@ -1650,11 +1717,8 @@ private: void reset_algorithm(const bool clearFaceLocations = true) { - Cone_tree_node* null_value=nullptr; - m_closestToVertices.resize(num_vertices(m_graph)); - std::fill(m_closestToVertices.begin(), m_closestToVertices.end(), Node_distance_pair(null_value, FT(0))); - m_vertexOccupiers.resize(num_halfedges(m_graph)); - std::fill(m_vertexOccupiers.begin(), m_vertexOccupiers.end(), Node_distance_pair(null_value, FT(0))); + m_closestToVertices.assign(num_vertices(m_graph), Node_distance_pair(nullptr, FT(-1))); + m_vertexOccupiers.assign(num_halfedges(m_graph), Node_distance_pair(nullptr, FT(-1))); while (!m_expansionPriqueue.empty()) { @@ -1671,7 +1735,7 @@ private: delete_all_nodes(); m_rootNodes.clear(); - m_vertexIsPseudoSource.resize(num_vertices(m_graph)); + m_vertexIsPseudoSource.assign(num_vertices(m_graph), false); #if !defined(NDEBUG) m_currentNodeCount = 0; @@ -1684,7 +1748,7 @@ private: } template - void visit_shortest_path(Cone_tree_node* startNode, + void visit_shortest_path(const Cone_tree_node* startNode, const Point_2& startLocation, Visitor& visitor) { @@ -1695,7 +1759,7 @@ private: typename Traits::Construct_target_2 construct_target_2(m_traits.construct_target_2_object()); typename Traits::Intersect_2 intersect_2(m_traits.intersect_2_object()); - Cone_tree_node* current = startNode; + const Cone_tree_node* current = startNode; Point_2 currentLocation(startLocation); while (!current->is_root_node()) @@ -1703,21 +1767,22 @@ private: switch (current->node_type()) { case Cone_tree_node::INTERVAL: - case Cone_tree_node::EDGE_SOURCE: { - Segment_2 entrySegment = current->entry_segment(); + const Segment_2& entrySegment = current->entry_segment(); const Point_2& currentSourceImage = current->source_image(); Ray_2 rayToLocation(construct_ray_2(currentSourceImage, currentLocation)); - const auto cgalIntersection = intersect_2(construct_line_2(entrySegment), construct_line_2(rayToLocation)); + const auto cgalIntersection = intersect_2(construct_line_2(entrySegment), + construct_line_2(rayToLocation)); CGAL_assertion(bool(cgalIntersection)); const Point_2* result = boost::get(&*cgalIntersection); + if (!result) + result = ¤tSourceImage; - if (!result) result = ¤tSourceImage; - - FT t0 = parametric_distance_along_segment_2(construct_source_2(entrySegment), construct_target_2(entrySegment), *result); + FT t0 = parametric_distance_along_segment_2(construct_source_2(entrySegment), + construct_target_2(entrySegment), *result); if (m_debugOutput) { @@ -1740,8 +1805,8 @@ private: std::cout << "Current Right Window: " << current->window_right() << " , " << m_traits.compute_parametric_distance_along_segment_2_object()(entrySegment.start(), entrySegment.end(), current->window_right()) << std::endl; std::cout << "Current Segment Intersection: " << *result << std::endl; - std::cout << "Edge: (" << get(m_vertexIndexMap, source(current->entry_edge(), m_graph)) - << "," << get(m_vertexIndexMap, target(current->entry_edge(), m_graph)) << ") : " << t0 << std::endl; + std::cout << "Edge: (V" << get(m_vertexIndexMap, source(current->entry_edge(), m_graph)) + << ", V" << get(m_vertexIndexMap, target(current->entry_edge(), m_graph)) << ") : " << t0 << std::endl; } visitor(current->entry_edge(), t0); @@ -1753,13 +1818,16 @@ private: } break; case Cone_tree_node::VERTEX_SOURCE: + // This might be a pseudo source visitor(target(current->entry_edge(), m_graph)); currentLocation = current->parent()->target_point(); current = current->parent(); break; + case Cone_tree_node::EDGE_SOURCE: case Cone_tree_node::FACE_SOURCE: // This is guaranteed to be the final node in any sequence - visitor(m_rootNodes[current->tree_id()].second->first, m_rootNodes[current->tree_id()].second->second); + visitor(m_rootNodes[current->tree_id()].second->first, + m_rootNodes[current->tree_id()].second->second); current = current->parent(); break; default: @@ -1866,7 +1934,7 @@ private: std::size_t associatedEdge; CGAL::Surface_mesh_shortest_paths_3::Barycentric_coordinates_type type; - boost::tie(type, associatedEdge) = classify_barycentric_coordinates(location); + std::tie(type, associatedEdge) = classify_barycentric_coordinates(location); switch (type) { @@ -1979,42 +2047,35 @@ private: reset_algorithm(false); set_vertex_types(); - m_vertexOccupiers.resize(num_halfedges(m_graph)); - m_closestToVertices.resize(num_vertices(m_graph)); + m_vertexOccupiers.assign(num_halfedges(m_graph), Node_distance_pair(nullptr, FT(-1))); + m_closestToVertices.assign(num_vertices(m_graph), Node_distance_pair(nullptr, FT(-1))); if (m_debugOutput) { - vertex_iterator current, end; - std::size_t numVertices = 0; - - for (boost::tie(current,end) = vertices(m_graph); current != end; ++current) + for (vertex_descriptor v : vertices(m_graph)) { std::cout << "Vertex#" << numVertices - << ": p = " << get(m_vertexPointMap,*current) - << " , Saddle Vertex: " << (is_saddle_vertex(*current) ? "yes" : "no") - << " , Boundary Vertex: " << (is_boundary_vertex(*current) ? "yes" : "no") << std::endl; + << ": p = " << get(m_vertexPointMap, v) + << " , Saddle Vertex: " << (is_saddle_vertex(v) ? "yes" : "no") + << " , Boundary Vertex: " << (is_boundary_vertex(v) ? "yes" : "no") << std::endl; ++numVertices; } } - face_iterator facesCurrent; - face_iterator facesEnd; - if (m_debugOutput) { std::size_t numFaces = 0; - - for (boost::tie(facesCurrent, facesEnd) = faces(m_graph); facesCurrent != facesEnd; ++facesCurrent) + for (face_descriptor f : faces(m_graph)) { std::cout << "Face#" << numFaces << ": Vertices = ("; ++numFaces; - halfedge_descriptor faceEdgesStart = halfedge(*facesCurrent, m_graph); + halfedge_descriptor faceEdgesStart = halfedge(f, m_graph); halfedge_descriptor faceEdgesCurrent = faceEdgesStart; do { - std::cout << get(m_vertexIndexMap, source(faceEdgesCurrent, m_graph)); + std::cout << "V" << get(m_vertexIndexMap, source(faceEdgesCurrent, m_graph)); faceEdgesCurrent = next(faceEdgesCurrent, m_graph); @@ -2031,15 +2092,15 @@ private: std::cout << std::endl; } - } for (typename Source_point_list::iterator it = m_faceLocations.begin(); it != m_faceLocations.end(); ++it) { if (m_debugOutput) { - std::cout << "Root: " << get(m_faceIndexMap, it->first) - << " , " << it->second[0] << " " << it->second[1] << " " << it->second[2] << " " << std::endl; + std::cout << "Root: F" << get(m_faceIndexMap, it->first) + << " , bar " << it->second[0] << " " << it->second[1] << " " << it->second[2] << " " + << ", pos " << point(it->first, it->second) << std::endl; } expand_root(it->first, it->second, Source_point_iterator(it)); @@ -2054,7 +2115,7 @@ private: } - while (m_expansionPriqueue.size() > 0) + while (!m_expansionPriqueue.empty()) { if (m_debugOutput) { @@ -2077,10 +2138,10 @@ private: if (!event->m_cancelled) { std::cout << " ------ Parent (" << event->m_parent << ") INFO: "; - std::cout << "EH = (" << get(m_vertexIndexMap, source(event->m_parent->entry_edge(), m_graph)) << " " - << get(m_vertexIndexMap, target(event->m_parent->entry_edge(), m_graph)) << ") "; - std::cout << "S = (" << event->m_parent->source_image() << ") "; - std::cout << "T = " << get(m_vertexIndexMap, target(next(event->m_parent->entry_edge(), m_graph), m_graph)); + std::cout << "EH = (V" << get(m_vertexIndexMap, source(event->m_parent->entry_edge(), m_graph)) << " V" + << get(m_vertexIndexMap, target(event->m_parent->entry_edge(), m_graph)) << ") "; + std::cout << "Src = (" << event->m_parent->source_image() << ") "; + std::cout << "Tar = V" << get(m_vertexIndexMap, target(next(event->m_parent->entry_edge(), m_graph), m_graph)); } std::cout << std::endl; @@ -2114,8 +2175,8 @@ private: if (m_debugOutput) { std::cout << "Left Expansion: Parent = " << parent - << " Edge = (" << get(m_vertexIndexMap, source(event->m_parent->left_child_edge(), m_graph)) - << "," << get(m_vertexIndexMap, target(event->m_parent->left_child_edge(), m_graph)) + << " Edge = (V" << get(m_vertexIndexMap, source(event->m_parent->left_child_edge(), m_graph)) + << ", V" << get(m_vertexIndexMap, target(event->m_parent->left_child_edge(), m_graph)) << ") , Distance = " << event->m_distanceEstimate << " , Level = " << event->m_parent->level() + 1 << std::endl; } @@ -2126,8 +2187,8 @@ private: if (m_debugOutput) { std::cout << "Right Expansion: Parent = " << parent - << " , Edge = (" << get(m_vertexIndexMap, source(event->m_parent->right_child_edge(), m_graph)) - << "," << get(m_vertexIndexMap, target(event->m_parent->right_child_edge(), m_graph)) + << " , Edge = (V" << get(m_vertexIndexMap, source(event->m_parent->right_child_edge(), m_graph)) + << ", V" << get(m_vertexIndexMap, target(event->m_parent->right_child_edge(), m_graph)) << ") , Distance = " << event->m_distanceEstimate << " , Level = " << event->m_parent->level() + 1 << std::endl; } @@ -2165,7 +2226,7 @@ private: for (std::size_t i = 0; i < m_closestToVertices.size(); ++i) { std::cout << "\tVertex = " << i << std::endl; - std::cout << "\tDistance = " << m_closestToVertices[i].second << std::endl; + std::cout << "\tDistance = " << m_closestToVertices[i].second << " to " << m_closestToVertices[i].first << std::endl; } std::cout << std::endl; @@ -2280,6 +2341,13 @@ public: Source_point_iterator add_source_point(vertex_descriptor v) { Face_location location = face_location(v); + + if (m_debugOutput) + { + std::cout << "Face location from V" << get(m_vertexIndexMap, v) << " is F" << get(m_faceIndexMap, location.first) << " " + << location.second[0] << " " << location.second[1] << " " << location.second[2] << std::endl; + } + return add_source_point(location); } @@ -2306,6 +2374,11 @@ public: */ Source_point_iterator add_source_point(const Face_location& location) { + if (m_debugOutput) + { + std::cout << "Add source point at position " << point(location.first, location.second) << std::endl; + } + Source_point_underlying_iterator added = m_faceLocations.insert(m_faceLocations.end(), location); if (m_firstNewSourcePoint == m_faceLocations.end()) @@ -2471,9 +2544,8 @@ public: { build_sequence_tree(); - Node_distance_pair result = m_closestToVertices[get(m_vertexIndexMap, v)]; - - Cone_tree_node* current = result.first; + const Node_distance_pair& result = m_closestToVertices[get(m_vertexIndexMap, v)]; + const Cone_tree_node* current = result.first; if (current) { @@ -2500,9 +2572,8 @@ public: { build_sequence_tree(); - std::pair result = nearest_to_location(f, location); - - Cone_tree_node* current = result.first.first; + const std::pair& result = nearest_to_location(f, location); + const Cone_tree_node* current = result.first.first; if (current) { @@ -2542,8 +2613,8 @@ public: { build_sequence_tree(); - Node_distance_pair result = m_closestToVertices[get(m_vertexIndexMap, v)]; - Cone_tree_node* current = result.first; + const Node_distance_pair& result = m_closestToVertices[get(m_vertexIndexMap, v)]; + const Cone_tree_node* current = result.first; if (current) { @@ -2737,13 +2808,23 @@ public: /*! \brief returns the 3-dimensional coordinates of the given vertex. - \param vertex A vertex of the input face graph + \param v A vertex of the input face graph */ - Point_3 point(const vertex_descriptor vertex) const + decltype(auto) point(const vertex_descriptor v) const { - return get(m_vertexPointMap, vertex); + return get(m_vertexPointMap, v); } + /// \cond + + static decltype(auto) point(const vertex_descriptor v, + const Triangle_mesh& tm) + { + return get(CGAL::vertex_point, tm, v); + } + + /// \endcond + /// @} /// \name Surface Face Location Constructions @@ -2770,7 +2851,7 @@ public: { typename Traits::Construct_barycentric_coordinates construct_barycentric_coordinates(traits.construct_barycentric_coordinates_object()); halfedge_descriptor hinit=halfedge(vertex, tm); - while (face(hinit, tm) == Graph_traits::null_face()) + while (is_border(hinit, tm)) hinit = opposite(next(hinit, tm), tm); halfedge_descriptor he = next(hinit, tm); @@ -3043,7 +3124,7 @@ public: Vertex_point_map vertexPointMap) { face_iterator facesStart, facesEnd; - boost::tie(facesStart, facesEnd) = faces(tm); + std::tie(facesStart, facesEnd) = faces(tm); outTree.rebuild(facesStart, facesEnd, tm, vertexPointMap); outTree.build(); } @@ -3054,6 +3135,4 @@ public: } // namespace CGAL -#include - #endif // CGAL_SURFACE_MESH_SHORTEST_PATH_SURFACE_MESH_SHORTEST_PATH_H diff --git a/Surface_mesh_shortest_path/include/CGAL/Surface_mesh_shortest_path/function_objects.h b/Surface_mesh_shortest_path/include/CGAL/Surface_mesh_shortest_path/function_objects.h index 6dd44eec600..19a2ed5105b 100644 --- a/Surface_mesh_shortest_path/include/CGAL/Surface_mesh_shortest_path/function_objects.h +++ b/Surface_mesh_shortest_path/include/CGAL/Surface_mesh_shortest_path/function_objects.h @@ -493,7 +493,7 @@ public: // In the case of multiple rays reaching the same target, we want to know their respective position // so that pruning of branches can be done according to the "one angle one split" idiom. // However, the orientation predicate is evaluated in the unfolded 2D plane, which is obtained - // via square roots; inconsisnties will exist. We don't want to prune in case it might be wrong, + // via square roots; inconsistencies will exist. We don't want to prune in case it might be wrong, // so we add a little bit of tolerance on the evaluation of the predicate. If it's almost collinear, // return 'collinear' (EQUAL). const FT eps = (FT(100) * std::numeric_limits::epsilon()); diff --git a/Surface_mesh_shortest_path/include/CGAL/Surface_mesh_shortest_path/internal/Cone_tree.h b/Surface_mesh_shortest_path/include/CGAL/Surface_mesh_shortest_path/internal/Cone_tree.h index b98f9d6c114..fe3161a5c24 100644 --- a/Surface_mesh_shortest_path/include/CGAL/Surface_mesh_shortest_path/internal/Cone_tree.h +++ b/Surface_mesh_shortest_path/include/CGAL/Surface_mesh_shortest_path/internal/Cone_tree.h @@ -95,15 +95,22 @@ private: public: Cone_tree_node(const Traits& traits, const Triangle_mesh& g, - const std::size_t treeId) + const halfedge_descriptor entryEdge, + const Triangle_2& layoutFace, + const Point_2& sourceImage, + const FT& pseudoSourceDistance, + const Point_2& windowLeft, + const Point_2& windowRight, + const Node_type nodeType = INTERVAL) : m_traits(traits) , m_graph(g) - , m_sourceImage(Point_2(CGAL::ORIGIN)) - , m_layoutFace(Point_2(CGAL::ORIGIN),Point_2(CGAL::ORIGIN),Point_2(CGAL::ORIGIN)) - , m_pseudoSourceDistance(0.0) - , m_level(0) - , m_treeId(treeId) - , m_nodeType(ROOT) + , m_entryEdge(entryEdge) + , m_sourceImage(sourceImage) + , m_layoutFace(layoutFace) + , m_pseudoSourceDistance(pseudoSourceDistance) + , m_windowLeft(windowLeft) + , m_windowRight(windowRight) + , m_nodeType(nodeType) , m_leftChild(nullptr) , m_rightChild(nullptr) , m_pendingLeftSubtree(nullptr) @@ -135,27 +142,8 @@ public: Cone_tree_node(const Traits& traits, const Triangle_mesh& g, - const halfedge_descriptor entryEdge, - const Triangle_2& layoutFace, - const Point_2& sourceImage, - const FT& pseudoSourceDistance, - const Point_2& windowLeft, - const Point_2& windowRight, - const Node_type nodeType = INTERVAL) - : m_traits(traits) - , m_graph(g) - , m_entryEdge(entryEdge) - , m_sourceImage(sourceImage) - , m_layoutFace(layoutFace) - , m_pseudoSourceDistance(pseudoSourceDistance) - , m_windowLeft(windowLeft) - , m_windowRight(windowRight) - , m_nodeType(nodeType) - , m_leftChild(nullptr) - , m_rightChild(nullptr) - , m_pendingLeftSubtree(nullptr) - , m_pendingRightSubtree(nullptr) - , m_pendingMiddleSubtree(nullptr) + const std::size_t treeId) + : Cone_tree_node(traits, g, treeId, Graph_traits::null_halfedge()) { } diff --git a/Surface_mesh_shortest_path/test/Surface_mesh_shortest_path/Surface_mesh_shortest_path_test_6.cpp b/Surface_mesh_shortest_path/test/Surface_mesh_shortest_path/Surface_mesh_shortest_path_test_6.cpp index 71d52570350..8ee21f598e9 100644 --- a/Surface_mesh_shortest_path/test/Surface_mesh_shortest_path/Surface_mesh_shortest_path_test_6.cpp +++ b/Surface_mesh_shortest_path/test/Surface_mesh_shortest_path/Surface_mesh_shortest_path_test_6.cpp @@ -1,36 +1,279 @@ -#include -#include -#include +#include #include #include -#include + +#include +#include +#include + +#include +#include +#include +#include typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel; + typedef CGAL::Surface_mesh Triangle_mesh; + typedef CGAL::Surface_mesh_shortest_path_traits Traits; typedef CGAL::Surface_mesh_shortest_path Surface_mesh_shortest_path; -typedef boost::graph_traits Graph_traits; -typedef Graph_traits::vertex_iterator vertex_iterator; -typedef Graph_traits::face_iterator face_iterator; +typedef boost::property_map::const_type VPM; +typedef CGAL::AABB_face_graph_triangle_primitive AABB_face_graph_primitive; +typedef CGAL::AABB_traits AABB_face_graph_traits; -int main() +void test_all_pairs() { CGAL::Surface_mesh mesh; std::ifstream input("data/test_mesh_6.off"); input >> mesh; input.close(); + std::cout << "Input mesh: " << num_vertices(mesh) << " nv" << std::endl; + for (Triangle_mesh::Vertex_index v1 : vertices(mesh)) + { + Surface_mesh_shortest_path shortest_paths(mesh); + shortest_paths.add_source_point(v1); + for (Triangle_mesh::Vertex_index v2 : vertices(mesh)) { - Surface_mesh_shortest_path shortest_paths(mesh); - shortest_paths.add_source_point(v1); - double dist = shortest_paths.shortest_distance_to_source_points(v2).first; - assert (dist==0 || v1!=v2); - CGAL_USE(dist); - } + const double sq_dist = CGAL::square(shortest_paths.shortest_distance_to_source_points(v2).first); + const double lower_bound = CGAL::squared_distance(mesh.point(v1), mesh.point(v2)); + std::cout << "sq_dist(" << v1 << ", " << v2 << ") = " << sq_dist << std::endl; + std::cout << "lower bound: " << lower_bound << std::endl; - return 0; + if(v1 == v2) + assert(sq_dist == 0.); + else + assert(sq_dist > (1. - 1e-7) * lower_bound); // numerical errors + + CGAL_USE(sq_dist); + } + } +} + +void test_flat_donut() +{ + CGAL::Surface_mesh mesh; + std::ifstream input("data/flat_donut.off"); + input >> mesh; + input.close(); + + std::cout << "Input mesh: " << num_vertices(mesh) << " nv" << std::endl; + + { + Triangle_mesh::Face_index f0(0), f16(16); + auto h0 = halfedge(f0, mesh), h16 = halfedge(f16, mesh); + + for(std::size_t i=0; i<3; ++i) + { + Triangle_mesh::Vertex_index v0(0); + Surface_mesh_shortest_path shortest_paths(mesh); + shortest_paths.add_source_point(v0); + + double dist = shortest_paths.shortest_distance_to_source_points(Triangle_mesh::Vertex_index(3)).first; + assert(dist == 3); + + dist = shortest_paths.shortest_distance_to_source_points(Triangle_mesh::Vertex_index(12)).first; + assert(dist == 3); + + dist = shortest_paths.shortest_distance_to_source_points(Triangle_mesh::Vertex_index(1)).first; + assert(dist == 1); + + dist = shortest_paths.shortest_distance_to_source_points(Triangle_mesh::Vertex_index(4)).first; + assert(dist == 1); + + dist = shortest_paths.shortest_distance_to_source_points(Triangle_mesh::Vertex_index(5)).first; + assert(CGAL::abs(dist - sqrt(2.)) < 1e-7); + + // change the canonical halfedges to test barycentric coordinates + h0 = next(h0, mesh); + set_halfedge(f0, h0, mesh); + h16 = next(h16, mesh); + set_halfedge(f16, h16, mesh); + } + } + + { + Triangle_mesh::Vertex_index v0(3); + + Surface_mesh_shortest_path shortest_paths(mesh); + shortest_paths.add_source_point(v0); + + double dist = shortest_paths.shortest_distance_to_source_points(Triangle_mesh::Vertex_index(0)).first; + assert(dist == 3); + + dist = shortest_paths.shortest_distance_to_source_points(Triangle_mesh::Vertex_index(15)).first; + assert(dist == 3); + + dist = shortest_paths.shortest_distance_to_source_points(Triangle_mesh::Vertex_index(2)).first; + assert(dist == 1); + + dist = shortest_paths.shortest_distance_to_source_points(Triangle_mesh::Vertex_index(7)).first; + assert(dist == 1); + + dist = shortest_paths.shortest_distance_to_source_points(Triangle_mesh::Vertex_index(6)).first; + assert(CGAL::abs(dist - sqrt(2.)) < 1e-7); + } + + { + Surface_mesh_shortest_path shortest_paths(mesh); + shortest_paths.add_source_point(Triangle_mesh::Vertex_index(0)); + shortest_paths.add_source_point(Triangle_mesh::Vertex_index(3)); + + double dist = shortest_paths.shortest_distance_to_source_points(Triangle_mesh::Vertex_index(0)).first; + assert(dist == 0); + + dist = shortest_paths.shortest_distance_to_source_points(Triangle_mesh::Vertex_index(15)).first; + assert(dist == 3); + + dist = shortest_paths.shortest_distance_to_source_points(Triangle_mesh::Vertex_index(2)).first; + assert(dist == 1); + + dist = shortest_paths.shortest_distance_to_source_points(Triangle_mesh::Vertex_index(7)).first; + assert(dist == 1); + + dist = shortest_paths.shortest_distance_to_source_points(Triangle_mesh::Vertex_index(6)).first; + assert(CGAL::abs(dist - sqrt(2.)) < 1e-7); + } + + { + Triangle_mesh::Face_index f16(16), f1(1), f29(29), f6(6); + auto h16 = halfedge(f16, mesh), h1 = halfedge(f1, mesh), h29 = halfedge(f29, mesh), h6 = halfedge(f6, mesh); + + for(std::size_t i=0; i<3; ++i) + { + Surface_mesh_shortest_path shortest_paths(mesh); + shortest_paths.add_source_point(Triangle_mesh::Vertex_index(1)); + + double dist = shortest_paths.shortest_distance_to_source_points(Triangle_mesh::Vertex_index(0)).first; + assert(dist == 1); + + dist = shortest_paths.shortest_distance_to_source_points(Triangle_mesh::Vertex_index(1)).first; + assert(dist == 0); + + auto loc = shortest_paths.locate(Kernel::Point_3(0.5, 0, 0)); + dist = shortest_paths.shortest_distance_to_source_points(loc.first, loc.second).first; + assert(dist == 0.5); + + loc = shortest_paths.locate(Kernel::Point_3(2.5, 0, 0)); + dist = shortest_paths.shortest_distance_to_source_points(loc.first, loc.second).first; + assert(dist == 1.5); + + dist = shortest_paths.shortest_distance_to_source_points(Triangle_mesh::Vertex_index(5)).first; + assert(dist == 1); + + dist = shortest_paths.shortest_distance_to_source_points(Triangle_mesh::Vertex_index(9)).first; + assert(dist == 2); + + dist = shortest_paths.shortest_distance_to_source_points(Triangle_mesh::Vertex_index(13)).first; + assert(dist == 3); + + dist = shortest_paths.shortest_distance_to_source_points(Triangle_mesh::Vertex_index(2)).first; + assert(dist == 1); + + h16 = next(h16, mesh); + set_halfedge(f16, h16, mesh); + h1 = next(h1, mesh); + set_halfedge(f1, h1, mesh); + h29 = next(h29, mesh); + set_halfedge(f29, h29, mesh); + h6 = next(h6, mesh); + set_halfedge(f6, h6, mesh); + } + } + + { + Triangle_mesh::Face_index f6(6); + auto h6 = halfedge(f6, mesh); + + for(std::size_t i=0; i<3; ++i) + { + Surface_mesh_shortest_path shortest_paths(mesh); + shortest_paths.add_source_point(shortest_paths.locate(Kernel::Point_3(1.5, 0, 0))); + + double dist = shortest_paths.shortest_distance_to_source_points(Triangle_mesh::Vertex_index(0)).first; + assert(dist == 1.5); + + dist = shortest_paths.shortest_distance_to_source_points(Triangle_mesh::Vertex_index(3)).first; + assert(dist == 1.5); + + dist = shortest_paths.shortest_distance_to_source_points(Triangle_mesh::Vertex_index(1)).first; + assert(dist == 0.5); + + dist = shortest_paths.shortest_distance_to_source_points(Triangle_mesh::Vertex_index(26)).first; + assert((dist - sqrt(2.)) < 1e-7); + + auto loc = shortest_paths.locate(Kernel::Point_3(1.5, 1, 0)); + dist = shortest_paths.shortest_distance_to_source_points(loc.first, loc.second).first; + assert(dist == 1); + + loc = shortest_paths.locate(Kernel::Point_3(0.5, 1, 0)); + dist = shortest_paths.shortest_distance_to_source_points(loc.first, loc.second).first; + assert(CGAL::abs(dist - sqrt(2.)) < 1e-7); + + h6 = next(h6, mesh); + set_halfedge(f6, h6, mesh); + } + } + + { + Triangle_mesh::Face_index f28(28); + auto h28 = halfedge(f28, mesh); + Triangle_mesh::Face_index f29(29); + auto h29 = halfedge(f29, mesh); + + auto h = halfedge(Triangle_mesh::Vertex_index(31), Triangle_mesh::Vertex_index(25), mesh).first; + + for(std::size_t side=0; side<2; ++side) + { + for(std::size_t i=0; i<3; ++i) + { + Surface_mesh_shortest_path shortest_paths(mesh); + shortest_paths.add_source_point(shortest_paths.face_location(h, 0.5)); + + double dist = shortest_paths.shortest_distance_to_source_points(Triangle_mesh::Vertex_index(16)).first; + assert(dist == 0.75); + + dist = shortest_paths.shortest_distance_to_source_points(Triangle_mesh::Vertex_index(31)).first; + assert(CGAL::abs(dist - 0.25) < 1e-7); // numerical issue (returns '0.25000000000000006') + + dist = shortest_paths.shortest_distance_to_source_points(Triangle_mesh::Vertex_index(25)).first; + assert(dist == 0.25); + + dist = shortest_paths.shortest_distance_to_source_points(Triangle_mesh::Vertex_index(23)).first; + assert(dist == 1.25); + + auto loc = shortest_paths.locate(Kernel::Point_3(1.25, 0, 0)); + dist = shortest_paths.shortest_distance_to_source_points(loc.first, loc.second).first; + assert(dist == 0.5); + + loc = shortest_paths.locate(Kernel::Point_3(1.25, 1, 0)); + dist = shortest_paths.shortest_distance_to_source_points(loc.first, loc.second).first; + assert(dist == 0.5); + + h28 = next(h28, mesh); + set_halfedge(f28, h28, mesh); + h29 = next(h29, mesh); + set_halfedge(f29, h29, mesh); + } + + h = opposite(h, mesh); + } + } +} + +int main(int, char**) +{ + std::cerr.precision(17); + std::cout.precision(17); + +// test_all_pairs(); + test_flat_donut(); + + std::cout << "Done!" << std::endl; + + return EXIT_SUCCESS; } diff --git a/Surface_mesh_shortest_path/test/Surface_mesh_shortest_path/data/flat_donut.off b/Surface_mesh_shortest_path/test/Surface_mesh_shortest_path/data/flat_donut.off new file mode 100644 index 00000000000..164a09d9a49 --- /dev/null +++ b/Surface_mesh_shortest_path/test/Surface_mesh_shortest_path/data/flat_donut.off @@ -0,0 +1,84 @@ +OFF +32 48 0 + +0 0 0 +1 0 0 +2 0 0 +3 0 0 +0 1 0 +1 1 0 +2 1 0 +3 1 0 +0 2 0 +1 2 0 +2 2 0 +3 2 0 +0 3 0 +1 3 0 +2 3 0 +3 3 0 +0.5 0.5 0 +0.5 1.5 0 +0.5 2.5 0 +1.5 2.5 0 +2.5 2.5 0 +2.5 1.5 0 +2.5 1 0 +2.5 0.5 0 +2 0.5 0 +1.5 0.5 0 +0.5 1 0 +0.5 2 0 +1 2.5 0 +2 2.5 0 +2.5 2 0 +1 0.5 0 +3 16 4 0 +3 31 16 1 +3 17 26 5 +3 9 17 5 +3 18 27 9 +3 28 18 9 +3 2 25 1 +3 24 25 2 +3 20 11 15 +3 20 30 11 +3 10 19 9 +3 29 19 10 +3 3 23 2 +3 7 23 3 +3 21 22 7 +3 11 21 7 +3 0 1 16 +3 31 26 16 +3 31 5 26 +3 26 4 16 +3 17 4 26 +3 17 8 4 +3 27 8 17 +3 9 27 17 +3 18 8 27 +3 18 12 8 +3 13 12 18 +3 28 13 18 +3 25 5 31 +3 1 25 31 +3 6 5 25 +3 24 6 25 +3 15 14 20 +3 20 29 30 +3 30 29 10 +3 20 14 29 +3 19 13 28 +3 9 19 28 +3 14 13 19 +3 29 14 19 +3 23 6 24 +3 2 23 24 +3 22 6 23 +3 7 22 23 +3 21 6 22 +3 21 10 6 +3 30 10 21 +3 11 30 21 +