songsenand
/
cgal
mirror of https://github.com/CGAL/cgal
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Merge remote-tracking branch 'cgal/5.4.x-branch' into HEAD

This commit is contained in:
Sébastien Loriot 2022-07-12 18:21:16 +02:00
parent 48f5785198 de2cd5c83f
commit 12795ba29c
17 changed files with 763 additions and 287 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
AABB_tree/doc/AABB_tree/Concepts/AABBRayIntersectionTraits.h
View File

@ -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 <A
HREF="http://www.siggraph.org/education/materials/HyperGraph/raytrace/rtinter3.htm">the
HREF="https://education.siggraph.org/static/HyperGraph/raytrace/rtinter3.htm">the
slab method</A>.
*/
typedef unspecified_type Intersection_distance;

2
Arrangement_on_surface_2/demo/Arrangement_on_surface_2/ArrangementPainterOstream.cpp
View File

@ -289,7 +289,7 @@ ArrangementPainterOstream<CGAL::Arr_linear_traits_2<Kernel_>>::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; }

2
Documentation/doc/Documentation/windows.txt
View File

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

2
Kernel_23/include/CGAL/Kernel_23/internal/Projection_traits_3.h
View File

@ -979,7 +979,7 @@ public:
Equal_x_2 eqx;
Equal_y_2 eqy;
return eqx(p,q) & eqy(p,q);
return eqx(p,q) && eqy(p,q);
}
};

3
Mesh_3/include/CGAL/Mesh_3/search_for_connected_components_in_labeled_image.h
View File

@ -73,8 +73,7 @@ search_for_connected_components_in_labeled_image(const CGAL::Image_3& image,
for(uint i=0; i<nx; i++)
{
using CGAL::IMAGEIO::static_evaluate;
if(visited[voxel_index] | second_pass[voxel_index]) {
if(visited[voxel_index] || second_pass[voxel_index]) {
++voxel_index;
continue;
}

2
Number_types/include/CGAL/GMP/Gmpzf_type.h
View File

@ -157,7 +157,7 @@ public:
return;
}
const int p = std::numeric_limits<double>::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:

2
Number_types/include/CGAL/MP_Float.h
View File

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

136
Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/stitch_borders.h
View File

@ -58,29 +58,46 @@ namespace internal {
////// Helper structs
// Used to compare halfedges based on their geometry
template <typename PolygonMesh, typename VertexPointMap>
template <typename PolygonMesh, typename VertexPointMap, typename GeomTraits>
struct Less_for_halfedge
{
typedef typename boost::graph_traits<PolygonMesh>::vertex_descriptor vertex_descriptor;
typedef typename boost::graph_traits<PolygonMesh>::halfedge_descriptor halfedge_descriptor;
typedef typename boost::property_traits<VertexPointMap>::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 Halfedge,
typename Border_halfedge_map,
typename Halfedge_pair,
typename Manifold_halfedge_pair,
typename Mesh,
typename VPM,
typename Mesh>
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<PolygonMesh, CGAL_PMP_NP_CLASS>::type GT;
GT gt = choose_parameter<GT>(get_parameter(np, internal_np::geom_traits));
typedef CGAL::dynamic_face_property_t<int> Face_property_tag;
typedef typename boost::property_map<PolygonMesh, Face_property_tag>::type Face_cc_map;
@ -386,10 +410,10 @@ OutputIterator collect_duplicated_stitchable_boundary_edges(const HalfedgeRange&
std::size_t num_cc = 0;
std::vector<std::vector<halfedge_descriptor> > border_edges_per_cc;
typedef Less_for_halfedge<PolygonMesh, VPM> Less_hedge;
typedef Less_for_halfedge<PolygonMesh, VPM, GT> Less_hedge;
typedef std::map<halfedge_descriptor, std::pair<int, std::size_t>, 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<std::pair<halfedge_descriptor, halfedge_descriptor> > 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 HalfedgeRange,
typename HalfedgeKeeper,
typename PolygonMesh,
typename VPM,
typename HalfedgeKeeper>
typename GT>
std::size_t zip_boundary_cycle(typename boost::graph_traits<PolygonMesh>::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<PolygonMesh>::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<PolygonMesh>::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<PolygonMesh>::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<PolygonMesh
VPM vpm = choose_parameter(get_parameter(np, internal_np::vertex_point),
get_const_property_map(vertex_point, pmesh));
typedef typename GetGeomTraits<PolygonMesh, CGAL_PMP_NP_CLASS>::type GT;
GT gt = choose_parameter<GT>(get_parameter(np, internal_np::geom_traits));
typedef typename internal_np::Lookup_named_param_def<internal_np::halfedges_keeper_t,
CGAL_NP_CLASS,
Default_halfedges_keeper<PolygonMesh> >::type Halfedge_keeper;
@ -1056,7 +1089,7 @@ std::size_t stitch_boundary_cycle(const typename boost::graph_traits<PolygonMesh
for(halfedge_descriptor h : halfedges_around_face(bh, pmesh))
cycle_halfedges.push_back(h);
std::size_t res = internal::zip_boundary_cycle(bh, cycle_halfedges, pmesh, vpm, hd_kpr);
std::size_t res = internal::zip_boundary_cycle(bh, cycle_halfedges, hd_kpr, pmesh, vpm, gt);
if(bh == boost::graph_traits<PolygonMesh>::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<PolygonMesh
/// \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
///
/// \returns the number of pairs of halfedges that were stitched.
@ -1172,6 +1216,17 @@ std::size_t stitch_boundary_cycles(const BorderHalfedgeRange& boundary_cycle_rep
/// \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
///
/// \returns the number of pairs of halfedges that were stitched.
@ -1412,15 +1467,6 @@ std::size_t stitch_borders(PolygonMesh& pmesh,
/// \param np optional sequence of \ref bgl_namedparameters "Named Parameters" among the ones listed below
///
/// \cgalNamedParamsBegin
/// \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<PolygonMesh>::%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<PolygonMesh>::%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.

4
STL_Extension/doc/STL_Extension/Concepts/Hashable.h
View File

@ -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<Key,Mapped,Hash>`
\sa <A HREF="http://www.cplusplus.com/reference/unordered_set/unordered_set/">`std::unordered_set`</a>
\sa <A HREF="http://www.cplusplus.com/reference/unordered_set/unordered_map/">`std::unordered_map`</a>
\sa <A HREF="https://en.cppreference.com/w/cpp/container/unordered_set">`std::unordered_set`</a>
\sa <A HREF="https://en.cppreference.com/w/cpp/container/unordered_map">`std::unordered_map`</a>
\sa <A HREF="https://www.boost.org/libs/unordered/doc/html/unordered.html#unordered_set">`boost::unordered_set`</a>
\sa <A HREF="https://www.boost.org/libs/unordered/doc/html/unordered.html#unordered_map">`boost::unordered_map`</a>

9
STL_Extension/include/CGAL/Uncertain.h
View File

@ -289,6 +289,11 @@ Uncertain<bool> operator!(Uncertain<bool> a)
return Uncertain<bool>(!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<bool> operator|(Uncertain<bool> a, Uncertain<bool> b)
{
@ -324,7 +329,9 @@ Uncertain<bool> operator&(Uncertain<bool> a, bool b)
{
return Uncertain<bool>(a.inf() & b, a.sup() & b);
}
#ifdef __clang__
# pragma GCC diagnostic pop
#endif
// Equality operators

2
Straight_skeleton_2/include/CGAL/constructions/Straight_skeleton_cons_ftC2.h
View File

@ -45,7 +45,7 @@ bool are_parallel_edges_equally_oriented( Segment_2_with_ID<K> const& e0, Segmen
template<class K>
bool are_edges_orderly_collinear( Segment_2_with_ID<K> const& e0, Segment_2_with_ID<K> 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);
}

36
Surface_mesh_shortest_path/examples/Surface_mesh_shortest_path/shortest_path_sequence.cpp
View File

@ -13,9 +13,11 @@
typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
typedef CGAL::Surface_mesh<Kernel::Point_3> Triangle_mesh;
typedef CGAL::Surface_mesh_shortest_path_traits<Kernel, Triangle_mesh> Traits;
typedef CGAL::Surface_mesh_shortest_path<Traits> Surface_mesh_shortest_path;
typedef Traits::Barycentric_coordinates Barycentric_coordinates;
typedef boost::graph_traits<Triangle_mesh> 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<halfedge_descriptor,double>& 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<face_descriptor, Barycentric_coordinates>& 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<int>(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);

445
Surface_mesh_shortest_path/include/CGAL/Surface_mesh_shortest_path/Surface_mesh_shortest_path.h
View File

@ -14,36 +14,29 @@
#include <CGAL/license/Surface_mesh_shortest_path.h>
#include <CGAL/disable_warnings.h>
#include <iterator>
#include <vector>
#include <utility>
#include <queue>
#include <algorithm>
#include <cstddef>
#include <list>
#include <boost/array.hpp>
#include <boost/lexical_cast.hpp>
#include <boost/utility/enable_if.hpp>
#include <boost/type_traits/is_same.hpp>
#include <iterator>
#include <CGAL/assertions.h>
#include <CGAL/AABB_tree.h>
#include <CGAL/Default.h>
#include <CGAL/enum.h>
#include <CGAL/number_utils.h>
#include <CGAL/Surface_mesh_shortest_path/barycentric.h>
#include <CGAL/Surface_mesh_shortest_path/internal/Cone_tree.h>
#include <CGAL/Surface_mesh_shortest_path/internal/misc_functions.h>
#include <CGAL/assertions.h>
#include <CGAL/AABB_tree.h>
#include <CGAL/boost/graph/helpers.h>
#include <CGAL/boost/graph/iterator.h>
#include <CGAL/Default.h>
#include <CGAL/enum.h>
#include <CGAL/number_utils.h>
#include <boost/lexical_cast.hpp>
#include <boost/variant/get.hpp>
#include <algorithm>
#include <cstddef>
#include <iterator>
#include <list>
#include <queue>
#include <utility>
#include <vector>
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 <<<<<<<<<<<<<<<<<<<" <<std::endl;
}
CGAL_assertion(cone->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 <class Visitor>
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<Point_2>(&*cgalIntersection);
if (!result)
result = &currentSourceImage;
if (!result) result = &currentSourceImage;
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<Node_distance_pair, Barycentric_coordinates> result = nearest_to_location(f, location);
Cone_tree_node* current = result.first.first;
const std::pair<Node_distance_pair, Barycentric_coordinates>& 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 <CGAL/enable_warnings.h>
#endif // CGAL_SURFACE_MESH_SHORTEST_PATH_SURFACE_MESH_SHORTEST_PATH_H

2
Surface_mesh_shortest_path/include/CGAL/Surface_mesh_shortest_path/function_objects.h
View File

@ -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<FT>::epsilon());

44
Surface_mesh_shortest_path/include/CGAL/Surface_mesh_shortest_path/internal/Cone_tree.h
View File

@ -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())
{
}

273
Surface_mesh_shortest_path/test/Surface_mesh_shortest_path/Surface_mesh_shortest_path_test_6.cpp
View File

@ -1,36 +1,279 @@
#include <iostream>
#include <fstream>
#include <stdlib.h>
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Surface_mesh.h>
#include <CGAL/Surface_mesh_shortest_path.h>
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/AABB_face_graph_triangle_primitive.h>
#include <CGAL/AABB_traits.h>
#include <CGAL/AABB_tree.h>
#include <iostream>
#include <fstream>
#include <limits>
#include <stdlib.h>
typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
typedef CGAL::Surface_mesh<Kernel::Point_3> Triangle_mesh;
typedef CGAL::Surface_mesh_shortest_path_traits<Kernel, Triangle_mesh> Traits;
typedef CGAL::Surface_mesh_shortest_path<Traits> Surface_mesh_shortest_path;
typedef boost::graph_traits<Triangle_mesh> Graph_traits;
typedef Graph_traits::vertex_iterator vertex_iterator;
typedef Graph_traits::face_iterator face_iterator;
typedef boost::property_map<Triangle_mesh, CGAL::vertex_point_t>::const_type VPM;
typedef CGAL::AABB_face_graph_triangle_primitive<Triangle_mesh, VPM> AABB_face_graph_primitive;
typedef CGAL::AABB_traits<Kernel, AABB_face_graph_primitive> AABB_face_graph_traits;
int main()
void test_all_pairs()
{
CGAL::Surface_mesh<Kernel::Point_3> 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<Kernel::Point_3> 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<AABB_face_graph_traits>(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<AABB_face_graph_traits>(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<AABB_face_graph_traits>(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<AABB_face_graph_traits>(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<AABB_face_graph_traits>(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<AABB_face_graph_traits>(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<AABB_face_graph_traits>(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;
}

84
Surface_mesh_shortest_path/test/Surface_mesh_shortest_path/data/flat_donut.off Normal file
View File

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