diff --git a/AABB_tree/doc/AABB_tree/Concepts/AABBGeomTraits.h b/AABB_tree/doc/AABB_tree/Concepts/AABBGeomTraits.h index b4b49a1f13b..a579453e9f2 100644 --- a/AABB_tree/doc/AABB_tree/Concepts/AABBGeomTraits.h +++ b/AABB_tree/doc/AABB_tree/Concepts/AABBGeomTraits.h @@ -65,7 +65,6 @@ A functor object to construct the sphere centered at one point and passing throu typedef unspecified_type Construct_sphere_3; /*! -\todo This is not correct! that is not used! A functor object to compute the point on a geometric primitive which is closest from a query. Provides the operator: `Point_3 operator()(const Type_2& type_2, const Point_3& p);` where `Type_2` is any type among `Segment_3` and `Triangle_3`. The operator returns the point on `type_2` which is closest to `p`. */ @@ -76,7 +75,7 @@ A functor object to compare the distance of two points wrt a third one. Provides the operator: `CGAL::Comparision_result operator()(const Point_3& p1, const Point_3& p2, const Point_3& p3)`. The operator compare the distance between `p1 and `p2`, and between `p2` and `p3`. */ -typedef unspecified_type Compare_distance_3 +typedef unspecified_type Compare_distance_3; /*! A functor object to detect if a point lies inside a sphere or not. diff --git a/Documentation/doc/biblio/geom.bib b/Documentation/doc/biblio/geom.bib index 63be4547305..0d135a1805f 100644 --- a/Documentation/doc/biblio/geom.bib +++ b/Documentation/doc/biblio/geom.bib @@ -151977,6 +151977,7 @@ pages = {179--189} year={2015} } + @incollection{bfhhm-epsph-15, title={Exact Minkowski sums of polygons with holes}, author={Baram, Alon and Fogel, Efi and Halperin, Dan and Hemmer, Michael and Morr, Sebastian}, @@ -151985,3 +151986,14 @@ pages = {179--189} year={2015}, publisher={Springer} } + +@inproceedings{cignoni1998metro, + title={Metro: measuring error on simplified surfaces}, + author={Cignoni, Paolo and Rocchini, Claudio and Scopigno, Roberto}, + booktitle={Computer Graphics Forum}, + volume={17}, + number={2}, + pages={167--174}, + year={1998}, + organization={Blackwell Publishers} +} diff --git a/Generator/doc/Generator/CGAL/point_generators_2.h b/Generator/doc/Generator/CGAL/point_generators_2.h index d73c1146cc9..9d5c1da1c07 100644 --- a/Generator/doc/Generator/CGAL/point_generators_2.h +++ b/Generator/doc/Generator/CGAL/point_generators_2.h @@ -377,7 +377,9 @@ typedef const Point_2& reference; \sa `std::random_shuffle` */ - template< typename Point_2, typename Triangulation > + template< typename Point_2, + typename Triangulation, + typename Creator = Creator_uniform_2::Kernel::RT,Point_2> > class Random_points_in_triangle_mesh_2 { public: @@ -446,7 +448,9 @@ get_default_random() ); \sa `std::random_shuffle` */ - template< typename Point_2 > + template< typename Point_2, + typename Triangle_2 = typename Kernel_traits::Kernel::Triangle_2, + typename Creator = Creator_uniform_2::Kernel::RT,Point_2> > class Random_points_in_triangles_2 { public: diff --git a/Generator/doc/Generator/CGAL/point_generators_3.h b/Generator/doc/Generator/CGAL/point_generators_3.h index 5691d6d983d..36dd1b5f0e2 100644 --- a/Generator/doc/Generator/CGAL/point_generators_3.h +++ b/Generator/doc/Generator/CGAL/point_generators_3.h @@ -248,7 +248,71 @@ get_default_random()); /// @} }; /* end Random_points_in_triangle_3 */ - +} /* end namespace CGAL */ + +namespace CGAL{ +/*! + +The class `Random_points_on_segment_3` is an input iterator creating points uniformly +distributed on a segment. The default `Creator` is +`Creator_uniform_3::Kernel::RT,Point_3>`. + +\cgalModels `InputIterator` +\cgalModels `PointGenerator` + +\sa `CGAL::cpp11::copy_n()` +\sa `CGAL::Counting_iterator` +\sa `std::random_shuffle` + +*/ +template< typename Point_3, typename Creator > +class Random_points_on_segment_3 { +public: + +/// \name Types +/// @{ + +/*! + +*/ +typedef std::input_iterator_tag iterator_category; + +/*! + +*/ +typedef Point_3 value_type; + +/*! + +*/ +typedef std::ptrdiff_t difference_type; + +/*! + +*/ +typedef const Point_3* pointer; + +/*! + +*/ +typedef const Point_3& reference; + + +/*! +creates an input iterator `g` generating points of type `Point_3` uniformly +distributed on the segment from \f$ p\f$ to \f$ q\f$ (excluding \f$ q\f$), +i.e.\ \f$ *g == (1-\lambda)\, p + \lambda q\f$ where \f$ 0 \le\lambda< 1\f$. +A single random number is needed from `rnd` for each point. +The expressions `to_double(p.x())`, `to_double(p.y())`, and `to_double(p.z())` must result +in the respective `double` representation of the coordinates of \f$ p\f$, and similarly for \f$ q\f$. +*/ +Random_points_on_segment_3( const Point_3& p, const Point_3& q, +Random& rnd = get_default_random()); + +/// @} + +}; /* end Random_points_on_segment_3 */ + } /* end namespace CGAL */ namespace CGAL { @@ -264,7 +328,6 @@ distributed inside a tetrahedron. The default `Creator` is \sa `CGAL::cpp11::copy_n()` \sa `CGAL::Counting_iterator` -\sa `CGAL::Random_points_in_disc_2` \sa `CGAL::Random_points_in_cube_3` \sa `CGAL::Random_points_in_triangle_3` \sa `CGAL::Random_points_on_sphere_3` @@ -344,13 +407,9 @@ The triangle range must be valid and unchanged while the iterator is used. \sa `CGAL::cpp11::copy_n()` \sa `CGAL::Counting_iterator` -\sa `CGAL::Points_on_segment_2` -\sa `CGAL::Random_points_in_disc_2` -\sa `CGAL::Random_points_on_segment_2` -\sa `CGAL::Random_points_on_square_2` \sa `CGAL::Random_points_in_cube_3` -\sa `CGAL::Random_points_in_triangle_3` -\sa `CGAL::Random_points_in_tetrahedron_3` +\sa `CGAL::Random_points_in_triangle_3` +\sa `CGAL::Random_points_in_tetrahedron_3` \sa `CGAL::Random_points_in_triangle_mesh_3` \sa `CGAL::Random_points_in_tetrahedral_mesh_boundary_3` \sa `CGAL::Random_points_in_tetrahedral_mesh_3` @@ -358,7 +417,10 @@ The triangle range must be valid and unchanged while the iterator is used. \sa `std::random_shuffle` */ -template< typename Point_3> +template< typename Point_3, + typename Triangle_3=typename Kernel_traits::Kernel::Triangle_3, + typename Creator = Creator_uniform_3< typename Kernel_traits< Point_3 >::Kernel::RT, + Point_3 > > class Random_points_in_triangles_3 { public: @@ -429,8 +491,12 @@ The triangle mesh must be valid and unchanged while the iterator is used. \sa `std::random_shuffle` */ -template < class TriangleMesh, class VertexPointMap = typename boost::property_map::type> > +template < class TriangleMesh, + class VertexPointMap = typename boost::property_map::type>, + class Creator = Creator_uniform_3< + typename Kernel_traits< typename boost::property_traits::value_type >::Kernel::RT, + typename boost::property_traits::value_type > > class Random_points_in_triangle_mesh_3 { public: @@ -476,7 +542,6 @@ Similar to the previous constructor using `get(vertex_point, mesh)` as vertex po */ Random_points_in_triangle_mesh_3(const TriangleMesh& mesh, Random& rnd = get_default_random() ); - /// @} }; /* end Random_points_in_triangle_mesh_3 */ @@ -509,7 +574,11 @@ C3T3 is a model of `Mesh_complex_3_in_triangulation_3` \sa `std::random_shuffle` */ -template +template ::Kernel::RT, + typename C3t3::Point > +> class Random_points_in_tetrahedral_mesh_boundary_3 { public: @@ -582,7 +651,10 @@ C3T3 is a model of `Mesh_complex_3_in_triangulation_3` \sa `std::random_shuffle` */ -template +template ::Kernel::RT, + typename C3t3::Point > > class Random_points_in_tetrahedral_mesh_3 { public: diff --git a/Generator/include/CGAL/internal/Generic_random_point_generator.h b/Generator/include/CGAL/internal/Generic_random_point_generator.h index 619c887c5db..0a5418cda68 100644 --- a/Generator/include/CGAL/internal/Generic_random_point_generator.h +++ b/Generator/include/CGAL/internal/Generic_random_point_generator.h @@ -115,7 +115,7 @@ template< class Functor > struct Apply_approx_sqrt: public Functor { template - double operator()(const T& t) const + typename cpp11::result_of::type operator()(const T& t) const { return approximate_sqrt( static_cast(*this)(t) ); } diff --git a/Generator/include/CGAL/point_generators_2.h b/Generator/include/CGAL/point_generators_2.h index a9ffb7b522f..25b94a2e0d7 100644 --- a/Generator/include/CGAL/point_generators_2.h +++ b/Generator/include/CGAL/point_generators_2.h @@ -505,7 +505,7 @@ class Random_points_in_triangle_2 : public Random_generator_base

{ void generate_point(); public: typedef P result_type; - typedef Random_points_in_triangle_2

This; + typedef Random_points_in_triangle_2 This; typedef typename Kernel_traits

::Kernel::Triangle_2 Triangle_2; Random_points_in_triangle_2() {} Random_points_in_triangle_2( const This& x,Random& rnd = get_default_random()) @@ -561,7 +561,11 @@ public: } }; }//end namespace internal -template +template ::Kernel::RT,P> +> class Random_points_in_triangle_mesh_2 : public Generic_random_point_generator< typename T::Face_handle , internal::Triangle_from_face_2, @@ -569,11 +573,11 @@ class Random_points_in_triangle_mesh_2 : public Generic_random_point_generator< public: typedef Generic_random_point_generator, - Random_points_in_triangle_2

, + Random_points_in_triangle_2, P> Base; typedef typename T::Face_handle Id; typedef P result_type; - typedef Random_points_in_triangle_mesh_2 This; + typedef Random_points_in_triangle_mesh_2 This; Random_points_in_triangle_mesh_2( const T& triangulation, Random& rnd = get_default_random()) @@ -621,7 +625,9 @@ struct Address_of { }//namesapce internal template ::Kernel::Triangle_2> + class Triangle_2=typename Kernel_traits::Kernel::Triangle_2, + class Creator = + Creator_uniform_2::Kernel::RT,Point_2> > struct Random_points_in_triangles_2 : public Generic_random_point_generator, @@ -630,11 +636,11 @@ struct Random_points_in_triangles_2 { typedef Generic_random_point_generator, - Random_points_in_triangle_2, + Random_points_in_triangle_2, Point_2> Base; typedef const Triangle_2* Id; typedef Point_2 result_type; - typedef Random_points_in_triangles_2 This; + typedef Random_points_in_triangles_2 This; template Random_points_in_triangles_2( const TriangleRange& triangles, Random& rnd = get_default_random()) diff --git a/Generator/include/CGAL/point_generators_3.h b/Generator/include/CGAL/point_generators_3.h index 395737b418c..31d6dee097b 100644 --- a/Generator/include/CGAL/point_generators_3.h +++ b/Generator/include/CGAL/point_generators_3.h @@ -205,7 +205,7 @@ class Random_points_in_triangle_3 : public Random_generator_base

{ void generate_point(); public: typedef P result_type; - typedef Random_points_in_triangle_3

This; + typedef Random_points_in_triangle_3 This; typedef typename Kernel_traits

::Kernel::Triangle_3 Triangle_3; Random_points_in_triangle_3() {} Random_points_in_triangle_3( const This& x,Random& rnd = get_default_random()) @@ -251,13 +251,78 @@ void Random_points_in_triangle_3::generate_point() { } template < class P, class Creator = + Creator_uniform_3::Kernel::RT,P> > +class Random_points_on_segment_3 : public Random_generator_base

{ + P _p; + P _q; + void generate_point(); +public: + typedef Random_points_on_segment_3 This; + Random_points_on_segment_3() {} + Random_points_on_segment_3( const P& p, + const P& q, + Random& rnd = CGAL::get_default_random()) + // g is an input iterator creating points of type `P' uniformly + // distributed on the segment from p to q except q, i.e. `*g' == + // \lambda p + (1-\lambda)\, q where 0 <= \lambda < 1 . A single + // random number is needed from `rnd' for each point. + : Random_generator_base

( (std::max)( (std::max)( (std::max)(to_double(p.x()), to_double(q.x())), + (std::max)(to_double(p.y()), to_double(q.y()))), + (std::max)(to_double(p.z()), to_double(q.z()))), + rnd) , _p(p), _q(q) + { + generate_point(); + } + + template + Random_points_on_segment_3( const Segment_3& s, + Random& rnd = CGAL::get_default_random()) + // g is an input iterator creating points of type `P' uniformly + // distributed on the segment from p to q except q, i.e. `*g' == + // \lambda p + (1-\lambda)\, q where 0 <= \lambda < 1 . A single + // random number is needed from `rnd' for each point. + : Random_generator_base

( (std::max)( (std::max)( (std::max)(to_double(s[0].x()), to_double(s[1].x())), + (std::max)(to_double(s[0].y()), to_double(s[1].y()))), + (std::max)(to_double(s[0].z()), to_double(s[1].z()))), + rnd) , _p(s[0]), _q(s[1]) + { + generate_point(); + } + + const P& source() const { return _p; } + const P& target() const { return _q; } + This& operator++() { + generate_point(); + return *this; + } + This operator++(int) { + This tmp = *this; + ++(*this); + return tmp; + } +}; + +template < class P, class Creator > +void +Random_points_on_segment_3:: +generate_point() { + typedef typename Creator::argument_type T; + double la = this->_rnd.get_double(); + double mu = 1.0 - la; + Creator creator; + this->d_item = creator(T(mu * to_double(_p.x()) + la * to_double(_q.x())), + T(mu * to_double(_p.y()) + la * to_double(_q.y())), + T(mu * to_double(_p.z()) + la * to_double(_q.z()))); +} + +template < class P, class Creator = Creator_uniform_3::Kernel::RT,P> > class Random_points_in_tetrahedron_3 : public Random_generator_base

{ P _p,_q,_r,_s; void generate_point(); public: typedef P result_type; - typedef Random_points_in_tetrahedron_3

This; + typedef Random_points_in_tetrahedron_3 This; typedef typename Kernel_traits

::Kernel::Tetrahedron_3 Tetrahedron_3; Random_points_in_tetrahedron_3() {} Random_points_in_tetrahedron_3( const This& x,Random& rnd = get_default_random()) @@ -306,42 +371,44 @@ void Random_points_in_tetrahedron_3::generate_point() { -template ::const_type> +template ::const_type, + class Creator = Creator_uniform_3< + typename Kernel_traits< typename boost::property_traits::value_type >::Kernel::RT, + typename boost::property_traits::value_type > +> struct Random_points_in_triangle_mesh_3 : public Generic_random_point_generator< typename boost::graph_traits ::face_descriptor , CGAL::Property_map_to_unary_function >, - Random_points_in_triangle_3::value_type>, + Random_points_in_triangle_3::value_type, Creator>, typename boost::property_traits::value_type> { typedef typename boost::property_traits::value_type P; typedef Generic_random_point_generator< typename boost::graph_traits ::face_descriptor , - CGAL::Property_map_to_unary_function >, - Random_points_in_triangle_3

, P> Base; + Random_points_in_triangle_3 , P> Base; typedef typename CGAL::Triangle_from_face_descriptor_map< - TriangleMesh,VertexPointMap> Pmap; - typedef typename CGAL::Triangle_from_face_descriptor_map< - TriangleMesh,VertexPointMap> Object_from_id_map; - typedef Random_points_in_triangle_3

Generator_on_object; + TriangleMesh,VertexPointMap> Object_from_id; typedef typename boost::graph_traits::face_descriptor Id; typedef P result_type; - typedef Random_points_in_triangle_mesh_3< TriangleMesh, VertexPointMap> This; + typedef Random_points_in_triangle_mesh_3< TriangleMesh, VertexPointMap, Creator> This; Random_points_in_triangle_mesh_3( const TriangleMesh& mesh,Random& rnd = get_default_random()) : Base( faces(mesh), - CGAL::Property_map_to_unary_function(Pmap(&mesh, get(vertex_point, mesh))), + CGAL::Property_map_to_unary_function(Object_from_id(&mesh, get(vertex_point, mesh))), internal::Apply_approx_sqrt::Kernel::Compute_squared_area_3>(), rnd ) { } Random_points_in_triangle_mesh_3( const TriangleMesh& mesh, VertexPointMap vpm, Random& rnd = get_default_random()) : Base( faces(mesh), - CGAL::Property_map_to_unary_function(Pmap(&mesh, vpm)), + CGAL::Property_map_to_unary_function(Object_from_id(&mesh, vpm)), internal::Apply_approx_sqrt::Kernel::Compute_squared_area_3>(), rnd ) { @@ -361,6 +428,62 @@ struct Random_points_in_triangle_mesh_3 } }; +template ::const_type, + class Creator = Creator_uniform_3< + typename Kernel_traits< typename boost::property_traits::value_type >::Kernel::RT, + typename boost::property_traits::value_type > +> +struct Random_points_on_edge_list_graph_3 + : public Generic_random_point_generator< + typename boost::graph_traits ::edge_descriptor, + CGAL::Property_map_to_unary_function >, + Random_points_on_segment_3::value_type, Creator>, + typename boost::property_traits::value_type> +{ + typedef typename boost::property_traits::value_type P; + typedef Generic_random_point_generator< + typename boost::graph_traits ::edge_descriptor, + CGAL::Property_map_to_unary_function >, + Random_points_on_segment_3 , P> Base; + typedef typename CGAL::Segment_from_edge_descriptor_map< + EdgeListGraph,VertexPointMap> Object_from_id; + typedef typename boost::graph_traits::edge_descriptor Id; + typedef P result_type; + typedef Random_points_on_edge_list_graph_3< EdgeListGraph, VertexPointMap, Creator> This; + + Random_points_on_edge_list_graph_3( const EdgeListGraph& mesh,Random& rnd = get_default_random()) + : Base( edges(mesh), + CGAL::Property_map_to_unary_function(Object_from_id(&mesh, get(vertex_point, mesh))), + internal::Apply_approx_sqrt::Kernel::Compute_squared_length_3>(), + rnd ) + { + } + Random_points_on_edge_list_graph_3( const EdgeListGraph& mesh, VertexPointMap vpm, Random& rnd = get_default_random()) + : Base( edges(mesh), + CGAL::Property_map_to_unary_function(Object_from_id(&mesh, vpm)), + internal::Apply_approx_sqrt::Kernel::Compute_squared_length_3>(), + rnd ) + { + } + This& operator++() { + Base::generate_point(); + return *this; + } + This operator++(int) { + This tmp = *this; + ++(*this); + return tmp; + } + double mesh_length() const + { + return this->sum_of_weights(); + } +}; + namespace internal { @@ -408,7 +531,11 @@ public: }; }//end namespace internal -template +template ::Kernel::RT, + typename C3t3::Point > +> struct Random_points_in_tetrahedral_mesh_boundary_3 : public Generic_random_point_generator< std::pair, @@ -419,11 +546,11 @@ struct Random_points_in_tetrahedral_mesh_boundary_3 typedef Generic_random_point_generator< std::pair, internal::Triangle_from_face_C3t3, - Random_points_in_triangle_3, + Random_points_in_triangle_3, typename C3t3::Point> Base; typedef std::pair Id; typedef typename C3t3::Point result_type; - typedef Random_points_in_tetrahedral_mesh_boundary_3 This; + typedef Random_points_in_tetrahedral_mesh_boundary_3 This; Random_points_in_tetrahedral_mesh_boundary_3( const C3t3& c3t3,Random& rnd = get_default_random()) @@ -445,7 +572,11 @@ struct Random_points_in_tetrahedral_mesh_boundary_3 } }; -template +template ::Kernel::RT, + typename C3t3::Point > +> struct Random_points_in_tetrahedral_mesh_3 : public Generic_random_point_generator< typename C3t3::Triangulation::Cell_handle, @@ -456,11 +587,11 @@ struct Random_points_in_tetrahedral_mesh_3 typedef Generic_random_point_generator< typename C3t3::Triangulation::Cell_handle, internal::Tetrahedron_from_cell_C3t3, - Random_points_in_tetrahedron_3, + Random_points_in_tetrahedron_3, typename C3t3::Point> Base; typedef typename C3t3::Triangulation::Cell_handle Id; typedef typename C3t3::Point result_type; - typedef Random_points_in_tetrahedral_mesh_3 This; + typedef Random_points_in_tetrahedral_mesh_3 This; Random_points_in_tetrahedral_mesh_3( const C3t3& c3t3,Random& rnd = get_default_random()) @@ -484,7 +615,11 @@ struct Random_points_in_tetrahedral_mesh_3 template ::Kernel::Triangle_3> + class Triangle_3=typename Kernel_traits::Kernel::Triangle_3, + class Creator = Creator_uniform_3< + typename Kernel_traits< Point_3 >::Kernel::RT, + Point_3 > + > struct Random_points_in_triangles_3 : public Generic_random_point_generator, @@ -493,11 +628,11 @@ struct Random_points_in_triangles_3 { typedef Generic_random_point_generator, - Random_points_in_triangle_3, + Random_points_in_triangle_3, Point_3> Base; typedef const Triangle_3* Id; typedef Point_3 result_type; - typedef Random_points_in_triangles_3 This; + typedef Random_points_in_triangles_3 This; template Random_points_in_triangles_3( const TriangleRange& triangles, Random& rnd = get_default_random()) diff --git a/Generator/test/Generator/test_generators.cpp b/Generator/test/Generator/test_generators.cpp index 9c4788ee9e3..e2c2598b5a8 100644 --- a/Generator/test/Generator/test_generators.cpp +++ b/Generator/test/Generator/test_generators.cpp @@ -115,13 +115,16 @@ void test_point_generators_3() { /* Create test point set. */ std::vector points; - points.reserve(500); + points.reserve(600); Random_points_in_sphere_3 g1( 100.0); CGAL::cpp11::copy_n( g1, 100, std::back_inserter(points)); Random_points_on_sphere_3 g2( 100.0); Random_points_in_cube_3 g3( 100.0); + Random_points_on_segment_3 g4( Point_3(-100,-100, -100), + Point_3( 100, 100, 100)); CGAL::cpp11::copy_n( g2, 100, std::back_inserter(points)); CGAL::cpp11::copy_n( g3, 100, std::back_inserter(points)); + CGAL::cpp11::copy_n( g4, 100, std::back_inserter(points)); points_on_cube_grid_3( 50.0, (std::size_t)1, std::back_inserter(points), Creator()); points_on_cube_grid_3( 50.0, (std::size_t)2, @@ -133,7 +136,7 @@ void test_point_generators_3() { random_selection( points.begin(), points.end(), 100, std::back_inserter(points)); - assert( points.size() == 500); + assert( points.size() == 600); for ( std::vector::iterator i = points.begin(); i != points.end(); i++){ assert( i->x() <= 100); diff --git a/Installation/changes.html b/Installation/changes.html index 45e144a1af8..da41ac19992 100644 --- a/Installation/changes.html +++ b/Installation/changes.html @@ -278,6 +278,16 @@ and src/ directories). point sets based on the assumption that they sample a curve in 2D or a surface in 3D. +

Polygon Mesh Processing

+
    +
  • + Add functions to compute approximated Hausdorff distances between two meshes, a mesh and a point set, or a point set and a mesh: + sample_triangle_mesh(), approximated_Hausdorff_distance(), + approximated_symmetric_Hausdorff_distance(), + max_distance_to_triangle_mesh(), + max_distance_to_point_set(). +
    • +
@@ -310,6 +320,7 @@ and src/ directories).
  • in a tetrahedral mesh (CGAL::Random_points_in_tetrahedral_mesh_3),
  • in a 2D triangle mesh (CGAL::Random_points_in_triangle_mesh_2),
  • in a range of 2D or 3D triangles (CGAL::Random_points_in_triangles_3 and CGAL::Random_points_in_triangles_2).
    • +
  • on a 3D segment (CGAL::Random_points_on_segment_3).
    • diff --git a/Poisson_surface_reconstruction_3/examples/Poisson_surface_reconstruction_3/poisson_reconstruction_example.cpp b/Poisson_surface_reconstruction_3/examples/Poisson_surface_reconstruction_3/poisson_reconstruction_example.cpp index 0453ddaf5bc..664f4de91c8 100644 --- a/Poisson_surface_reconstruction_3/examples/Poisson_surface_reconstruction_3/poisson_reconstruction_example.cpp +++ b/Poisson_surface_reconstruction_3/examples/Poisson_surface_reconstruction_3/poisson_reconstruction_example.cpp @@ -12,6 +12,8 @@ #include #include +#include + #include #include @@ -104,5 +106,15 @@ int main(void) CGAL::output_surface_facets_to_polyhedron(c2t3, output_mesh); out << output_mesh; + + /// [PMP_distance_snippet] + // computes the approximation error of the reconstruction + double max_dist = + CGAL::Polygon_mesh_processing::approximate_max_distance_to_point_set(output_mesh, + points, + 4000); + std::cout << "Max distance to point_set: " << max_dist << std::endl; + /// [PMP_distance_snippet] + return EXIT_SUCCESS; } diff --git a/Polygon_mesh_processing/doc/Polygon_mesh_processing/Concepts/PMPDistanceTraits.h b/Polygon_mesh_processing/doc/Polygon_mesh_processing/Concepts/PMPDistanceTraits.h new file mode 100644 index 00000000000..0009489dcec --- /dev/null +++ b/Polygon_mesh_processing/doc/Polygon_mesh_processing/Concepts/PMPDistanceTraits.h @@ -0,0 +1,63 @@ +/// \ingroup PkgPolygonMeshProcessingConcepts +/// \cgalConcept +/// +/// The concept `PMPDistanceTraits` is a refinement of the +/// concepts `AABBGeomTraits` and `SpatialSortingTraits_3`. In addition to the types required by +/// those concepts, it also requires types and functors needed by the functions `approximate_max_distance_to_point_set()`, +/// `sample_triangle_mesh()`, `approximate_Hausdorff_distance()` and `max_distance_to_triangle_mesh()` +/// +/// \cgalRefines `AABBGeomTraits` +/// \cgalRefines `SpatialSortingTraits_3` +/// \cgalHasModel Any 3D Kernel is a model of this concept. + +class PMPDistanceTraits{ +public: + /*! + * A number type model of `Field` and `RealEmbeddable` + */ + typedef unspecified_type FT; + /// + /*! 3D point type + * It must be default constructible, and can be constructed from 3 objects of type `FT`. + * `bool operator<(Point_3, Point_3)` to lexicographically compare two points must be available. + * Access to Cartesian coordinates must be possible using `Point_3::x()`, `Point_3::y(), Point_3::z()` and + * `FT operator[](int i)` with `0 <= i < 3`. + * + * There must be a specialization of `CGAL::Kernel_traits` such that + * `CGAL::Kernel_traits::%Kernel` is a model implementing this concept. + */ + typedef unspecified_type Point_3; + + /// 3D vector type + typedef unspecified_type Vector_3; + + /// @name Functors + /// @{ + + /// Functor for computing squared area of a triangle. + /// It provides `FT operator()(const Point_3&, const Point_3&, const Point_3&) const` + /// and `FT operator()(const Triangle_3&) const` and has `FT` as result_type. + typedef unspecified_type Compute_squared_area_3; + /// Functor for computing squared length of a segment. + /// and `FT operator()(const Segment_3&) const` and has `FT` as result_type. + typedef unspecified_type Compute_squared_length_3; + /// Functor for constructing translated points. + /// It provides `Point_3 operator()(const Point_3 &, const Vector_3 &)` + typedef unspecified_type Construct_translated_point_3; + /// Functor for constructing vectors. + /// It provides `Vector_3 operator()(const Point_3 &, const Point_3 &)` + typedef unspecified_type Construct_vector_3; + /// Functor for constructing scaled vectors. + /// It provides `Vector_3 operator()(const Vector_3 &, const FT &)` + typedef unspecified_type Construct_scaled_vector_3; + /// @} + + /// @name Functions + /// @{ + Compute_squared_area_3 compute_squared_area_3_object(); + Compute_squared_length_3 compute_squared_length_3_object(); + Construct_translated_point_3 construct_translated_point_3_object(); + Construct_vector_3 construct_vector_3_object(); + Construct_scaled_vector_3 construct_scaled_vector_3_object(); + /// @} +}; diff --git a/Polygon_mesh_processing/doc/Polygon_mesh_processing/Concepts/SelfIntersectionTraits.h b/Polygon_mesh_processing/doc/Polygon_mesh_processing/Concepts/PMPSelfIntersectionTraits.h similarity index 97% rename from Polygon_mesh_processing/doc/Polygon_mesh_processing/Concepts/SelfIntersectionTraits.h rename to Polygon_mesh_processing/doc/Polygon_mesh_processing/Concepts/PMPSelfIntersectionTraits.h index e50813019ad..f71561f91d9 100644 --- a/Polygon_mesh_processing/doc/Polygon_mesh_processing/Concepts/SelfIntersectionTraits.h +++ b/Polygon_mesh_processing/doc/Polygon_mesh_processing/Concepts/PMPSelfIntersectionTraits.h @@ -2,7 +2,7 @@ /// \cgalConcept /// /// Geometric traits concept for the functions `CGAL::self_intersections()` and `CGAL::does_self_intersect()`. -class SelfIntersectionTraits{ +class PMPSelfIntersectionTraits{ public: /// @name Geometric Types /// @{ diff --git a/Polygon_mesh_processing/doc/Polygon_mesh_processing/Doxyfile.in b/Polygon_mesh_processing/doc/Polygon_mesh_processing/Doxyfile.in index 71e185e06ef..b67c5ca856c 100644 --- a/Polygon_mesh_processing/doc/Polygon_mesh_processing/Doxyfile.in +++ b/Polygon_mesh_processing/doc/Polygon_mesh_processing/Doxyfile.in @@ -18,6 +18,8 @@ ALIASES += "cgalNPTableEnd= " ALIASES += "cgalNPBegin{1}=\1 " ALIASES += "cgalNPEnd=" +EXAMPLE_PATH += ${CGAL_Poisson_surface_reconstruction_3_EXAMPLE_DIR} + MACRO_EXPANSION = YES EXPAND_ONLY_PREDEF = YES EXPAND_AS_DEFINED = CGAL_PMP_NP_TEMPLATE_PARAMETERS \ diff --git a/Polygon_mesh_processing/doc/Polygon_mesh_processing/NamedParameters.txt b/Polygon_mesh_processing/doc/Polygon_mesh_processing/NamedParameters.txt index 2ce6592baef..d8aef11f859 100644 --- a/Polygon_mesh_processing/doc/Polygon_mesh_processing/NamedParameters.txt +++ b/Polygon_mesh_processing/doc/Polygon_mesh_processing/NamedParameters.txt @@ -210,7 +210,113 @@ set to `true`, this parameter is ignored. \n \b Type : `bool` \n \b Default value is `true` +\cgalNPEnd +\cgalNPBegin{use_random_uniform_sampling} \anchor PMP_use_random_uniform_sampling +Parameter used in `sample_triangle_mesh()` to indicate if points should be picked +in a random uniform way. +\n +\b Type : `bool` \n +\b Default value is `true` +\cgalNPEnd + +\cgalNPBegin{use_grid_sampling} \anchor PMP_use_grid_sampling +Parameter used in `sample_triangle_mesh()` to indicate if points should be picked +in on a grid in each face. +\n +\b Type : `bool` \n +\b Default value is `false` +\cgalNPEnd + +\cgalNPBegin{use_monte_carlo_sampling} \anchor PMP_use_monte_carlo_sampling +Parameter used in `sample_triangle_mesh()` to indicate if points should be picked +using a Monte-Carlo approach. +\n +\b Type : `bool` \n +\b Default value is `false` +\cgalNPEnd + +\cgalNPBegin{sample_edges} \anchor PMP_sample_edges +Parameter used in `sample_triangle_mesh()` to indicate if a dedicated sampling +of edges should be done. +\n +\b Type : `bool` \n +\b Default value is `true` +\cgalNPEnd + +\cgalNPBegin{sample_vertices} \anchor PMP_sample_vertices +Parameter used in `sample_triangle_mesh()` to indicate if triangle vertices should +be copied in the output iterator. +\n +\b Type : `bool` \n +\b Default value is `true` +\cgalNPEnd + +\cgalNPBegin{sample_faces} \anchor PMP_sample_faces +Parameter used in `sample_triangle_mesh()` to indicate if the interior of faces +should be considered for the sampling. +\n +\b Type : `bool` \n +\b Default value is `true` +\cgalNPEnd + +\cgalNPBegin{number_of_points_on_faces} \anchor PMP_number_of_points_on_faces +Parameter used in `sample_triangle_mesh()` to set the number of points picked +using the random uniform method on faces. +\n +\b Type : `std::size_t` \n +\b Default value is `0` +\cgalNPEnd + +\cgalNPBegin{number_of_points_on_edges} \anchor PMP_number_of_points_on_edges +Parameter used in `sample_triangle_mesh()` to set the number of points picked +using the random uniform method on edges. +\n +\b Type : `std::size_t` \n +\b Default value is `0` +\cgalNPEnd + +\cgalNPBegin{number_of_points_per_face} \anchor PMP_number_of_points_per_face +Parameter used in `sample_triangle_mesh()` to set the number of points picked +per face using the Monte-Carlo method. +\n +\b Type : `std::size_t` \n +\b Default value is `0` +\cgalNPEnd + +\cgalNPBegin{number_of_points_per_edge} \anchor PMP_number_of_points_per_edge +Parameter used in `sample_triangle_mesh()` to set the number of points picked +per edge using the Monte-Carlo method. +\n +\b Type : `std::size_t` \n +\b Default value is `0` +\cgalNPEnd + +\cgalNPBegin{grid_spacing} \anchor PMP_grid_spacing +Parameter used in `sample_triangle_mesh()` to set the grid spacing when using +the grid sampling method. +\n +\b Type : `double` \n +\b Default value is `0` +\cgalNPEnd + +\cgalNPBegin{number_of_points_per_area_unit} \anchor PMP_number_of_points_per_area_unit +Parameter used in `sample_triangle_mesh()` to set the number of points per +area unit to be picked up in faces for the random uniform sampling and +Monte-Carlo methods. +\n +\b Type : `double` \n +\b Default value is `0` +\cgalNPEnd + +\cgalNPBegin{number_of_points_per_distance_unit} \anchor PMP_number_of_points_per_distance_unit +Parameter used in `sample_triangle_mesh()` to set the number of points per +distance unit to be picked up on edges for the random uniform sampling and +Monte-Carlo methods. +\n +\b Type : `double` \n +\b Default value is `0` +\cgalNPEnd \cgalNPTableEnd diff --git a/Polygon_mesh_processing/doc/Polygon_mesh_processing/PackageDescription.txt b/Polygon_mesh_processing/doc/Polygon_mesh_processing/PackageDescription.txt index 2ab97117884..dcaec250401 100644 --- a/Polygon_mesh_processing/doc/Polygon_mesh_processing/PackageDescription.txt +++ b/Polygon_mesh_processing/doc/Polygon_mesh_processing/PackageDescription.txt @@ -40,6 +40,12 @@ /// Functions to orient polygon soups and to stitch geometrically identical boundaries. /// \ingroup PkgPolygonMeshProcessing +/// \defgroup PMP_distance_grp Distance Functions +/// Functions to compute the distance between meshes, between a mesh and a point set and between a point set and a mesh. +/// \ingroup PkgPolygonMeshProcessing + + + /*! \addtogroup PkgPolygonMeshProcessing @@ -126,6 +132,13 @@ and provides a list of the parameters that are used in this package. - \link measure_grp `CGAL::Polygon_mesh_processing::edge_length()` \endlink - \link measure_grp `CGAL::Polygon_mesh_processing::face_border_length()` \endlink +## Distance Functions ## +- `CGAL::Polygon_mesh_processing::approximate_Hausdorff_distance()` +- `CGAL::Polygon_mesh_processing::approximate_symmetric_Hausdorff_distance()` +- `CGAL::Polygon_mesh_processing::approximate_max_distance_to_point_set()` +- `CGAL::Polygon_mesh_processing::max_distance_to_triangle_mesh()` +- `CGAL::Polygon_mesh_processing::sample_triangle_mesh()` + ## Miscellaneous ## - `CGAL::Polygon_mesh_slicer` - `CGAL::Side_of_triangle_mesh` diff --git a/Polygon_mesh_processing/doc/Polygon_mesh_processing/Polygon_mesh_processing.txt b/Polygon_mesh_processing/doc/Polygon_mesh_processing/Polygon_mesh_processing.txt index 0ae1b22cb07..0ed5a16889b 100644 --- a/Polygon_mesh_processing/doc/Polygon_mesh_processing/Polygon_mesh_processing.txt +++ b/Polygon_mesh_processing/doc/Polygon_mesh_processing/Polygon_mesh_processing.txt @@ -481,11 +481,53 @@ the propagation of a connected component index to cross it. \cgalExample{Polygon_mesh_processing/connected_components_example.cpp} +\section PMPDistance Approximate Hausdorff Distance + +This package provides methods to compute (approximate) distances between meshes and point sets. + +The function \link CGAL::Polygon_mesh_processing::approximate_Hausdorff_distance() `approximate_Hausdorff_distance()`\endlink +computes an approximation of the Hausdorff distance from a mesh `tm1` to a mesh `tm2`. Given a +a sampling of `tm1`, it computes the distance to `tm2` of the farthest sample point to `tm2` \cgalCite{cignoni1998metro}. +The symmetric version (\link CGAL::Polygon_mesh_processing::approximate_symmetric_Hausdorff_distance() `approximate_symmetric_Hausdorff_distance()`\endlink) +is the maximum of the two non-symmetric distances. Internally, points are sampled using +\link CGAL::Polygon_mesh_processing::sample_triangle_mesh() `sample_triangle_mesh()`\endlink and the distance +to each sample point is computed using +\link CGAL::Polygon_mesh_processing::max_distance_to_triangle_mesh() `max_distance_to_triangle_mesh()`\endlink. +The quality of the approximation depends on the quality of the sampling and the runtime depends on the number of sample points. +Three sampling methods with different parameters are provided (see \cgalFigureRef{sampling_bunny}). + +\cgalFigureBegin{sampling_bunny, pmp_sampling_bunny.jpg} +Sampling of a triangle mesh using different sampling methods. From left to right: (a) Grid sampling, (b) Monte-Carlo sampling with fixed number of points per face and per edge, +(c) Monte-Carlo sampling with a number of points proportional to the area/length, and (d) Uniform random sampling. +The four pictures represent the sampling on the same portion of a mesh, parameters were adjusted so that the total number of points sampled in faces (blue points) and on +edges (red points) are roughly the same. +Note that when using the random uniform sampling some faces/edges may not contain any point, but this method is the only one that allows to exactly match a given number of points. +\cgalFigureEnd + +The function \link CGAL::Polygon_mesh_processing::approximate_max_distance_to_point_set() `approximate_max_distance_to_point_set()`\endlink +computes an approximation of the Hausdorff distance from a mesh to a point set. +For each triangle, a lower and upper bound of the Hausdorff distance to the point set are computed. +Triangles are refined until the difference between the bounds is lower than a user-defined precision threshold. + +\subsection AHDExample Approximate Hausdorff Distance Example +In the following example, a mesh is isotropically remeshed and the approximate distance between the input and the output is computed. + +\cgalExample{Polygon_mesh_processing/hausdorff_distance_remeshing_example.cpp} + +\subsection PoissonDistanceExample Max Distance Between Point Set and Surface Example +In \ref Poisson_surface_reconstruction_3/poisson_reconstruction_example.cpp, +a triangulated surface mesh is constructed from a point set using the +\link PkgPoissonSurfaceReconstructionSummary Poisson reconstruction algorithm \endlink, +and the distance between the point set and the reconstructed surface is computed +with the following code: + +\snippet Poisson_surface_reconstruction_3/poisson_reconstruction_example.cpp PMP_distance_snippet \section PMPHistory Implementation History A first version of this package was started by Ilker %O. Yaz and Sébastien Loriot. Jane Tournois worked on the finalization of the API, code, and documentation. + */ } /* namespace CGAL */ diff --git a/Polygon_mesh_processing/doc/Polygon_mesh_processing/dependencies b/Polygon_mesh_processing/doc/Polygon_mesh_processing/dependencies index f7787717148..52428e7009d 100644 --- a/Polygon_mesh_processing/doc/Polygon_mesh_processing/dependencies +++ b/Polygon_mesh_processing/doc/Polygon_mesh_processing/dependencies @@ -11,3 +11,4 @@ Surface_mesh Surface_mesh_deformation AABB_tree Triangulation_2 +Spatial_sorting diff --git a/Polygon_mesh_processing/doc/Polygon_mesh_processing/examples.txt b/Polygon_mesh_processing/doc/Polygon_mesh_processing/examples.txt index 60587796e37..2c1b138146e 100644 --- a/Polygon_mesh_processing/doc/Polygon_mesh_processing/examples.txt +++ b/Polygon_mesh_processing/doc/Polygon_mesh_processing/examples.txt @@ -13,4 +13,6 @@ \example Polygon_mesh_processing/mesh_slicer_example.cpp \example Polygon_mesh_processing/isotropic_remeshing_example.cpp \example Polygon_mesh_processing/compute_normals_example_Polyhedron.cpp +\example Polygon_mesh_processing/hausdorff_distance_remeshing_example.cpp +\example Poisson_surface_reconstruction_3/poisson_reconstruction_example.cpp */ diff --git a/Polygon_mesh_processing/doc/Polygon_mesh_processing/fig/pmp_sampling_bunny.jpg b/Polygon_mesh_processing/doc/Polygon_mesh_processing/fig/pmp_sampling_bunny.jpg new file mode 100644 index 00000000000..b81b3d286ad Binary files /dev/null and b/Polygon_mesh_processing/doc/Polygon_mesh_processing/fig/pmp_sampling_bunny.jpg differ diff --git a/Polygon_mesh_processing/examples/Polygon_mesh_processing/CMakeLists.txt b/Polygon_mesh_processing/examples/Polygon_mesh_processing/CMakeLists.txt index 0aea8b58750..0a838cc3091 100644 --- a/Polygon_mesh_processing/examples/Polygon_mesh_processing/CMakeLists.txt +++ b/Polygon_mesh_processing/examples/Polygon_mesh_processing/CMakeLists.txt @@ -58,11 +58,7 @@ endif() # include for local package include_directories( BEFORE ../../include ) -find_package(Eigen3 3.1.0) #(requires 3.1.0 or greater) - - if (EIGEN3_FOUND) - # Executables that require Eigen 3.1 - include( ${EIGEN3_USE_FILE} ) +find_package(Eigen3 3.2.0) #(requires 3.2.0 or greater) # Creating entries for all .cpp/.C files with "main" routine @@ -70,6 +66,23 @@ find_package(Eigen3 3.1.0) #(requires 3.1.0 or greater) include( CGAL_CreateSingleSourceCGALProgram ) +find_package( TBB ) +if( TBB_FOUND ) + include( ${TBB_USE_FILE} ) + list( APPEND CGAL_3RD_PARTY_LIBRARIES ${TBB_LIBRARIES} ) +else() + message( STATUS "NOTICE: Intel TBB was not found. hausdorff_distance_example will use sequential code." ) +endif() +create_single_source_cgal_program( "hausdorff_distance_remeshing_example.cpp") +if (EIGEN3_FOUND) + # Executables that require Eigen 3.2 + include( ${EIGEN3_USE_FILE} ) + + create_single_source_cgal_program( "hole_filling_example.cpp" ) + create_single_source_cgal_program( "hole_filling_example_SM.cpp" ) + create_single_source_cgal_program( "refine_fair_example.cpp") +endif(EIGEN3_FOUND) + create_single_source_cgal_program( "self_intersections_example.cpp" ) create_single_source_cgal_program( "stitch_borders_example.cpp" ) create_single_source_cgal_program( "compute_normals_example_Polyhedron.cpp" CXX_FEATURES cxx_range_for ) @@ -83,13 +96,6 @@ create_single_source_cgal_program( "mesh_slicer_example.cpp") #create_single_source_cgal_program( "remove_degeneracies_example.cpp") create_single_source_cgal_program( "isotropic_remeshing_example.cpp") -if(NOT (${EIGEN3_VERSION} VERSION_LESS 3.2.0)) -create_single_source_cgal_program( "hole_filling_example.cpp" ) -create_single_source_cgal_program( "hole_filling_example_SM.cpp" ) -create_single_source_cgal_program( "refine_fair_example.cpp") -else() - message(STATUS "NOTICE: Some examples require Eigen 3.2 (or greater) and will not be compiled.") -endif() if(OpenMesh_FOUND) create_single_source_cgal_program( "compute_normals_example_OM.cpp" ) @@ -111,7 +117,5 @@ create_single_source_cgal_program( "triangulate_faces_example_OM.cpp") target_link_libraries( triangulate_faces_example_OM ${OPENMESH_LIBRARIES} ) endif(OpenMesh_FOUND) -else(EIGEN3_FOUND) - message(STATUS "NOTICE: Some examples require Eigen 3.1 (or greater) and will not be compiled.") -endif(EIGEN3_FOUND) + diff --git a/Polygon_mesh_processing/examples/Polygon_mesh_processing/hausdorff_distance_remeshing_example.cpp b/Polygon_mesh_processing/examples/Polygon_mesh_processing/hausdorff_distance_remeshing_example.cpp new file mode 100644 index 00000000000..bbd7b92d3f0 --- /dev/null +++ b/Polygon_mesh_processing/examples/Polygon_mesh_processing/hausdorff_distance_remeshing_example.cpp @@ -0,0 +1,34 @@ +#include +#include +#include +#include +#include + +#if CGAL_LINKED_WITH_TBB +#define TAG CGAL::Parallel_tag +#else +#define TAG CGAL::Sequential_tag +#endif + +typedef CGAL::Exact_predicates_inexact_constructions_kernel K; +typedef K::Point_3 Point; +typedef CGAL::Surface_mesh Mesh; + +namespace PMP = CGAL::Polygon_mesh_processing; + +int main() +{ + Mesh tm1, tm2; + CGAL::make_tetrahedron(Point(.0,.0,.0), + Point(2,.0,.0), + Point(1,1,1), + Point(1,.0,2), + tm1); + tm2=tm1; + CGAL::Polygon_mesh_processing::isotropic_remeshing(tm2.faces(),.05, tm2); + + std::cout << "Approximated Hausdorff distance: " + << CGAL::Polygon_mesh_processing::approximate_Hausdorff_distance + (tm1, tm2, PMP::parameters::number_of_points_per_area_unit(4000)) + << std::endl; +} diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/distance.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/distance.h index ba90c8e4031..c3c0217f0dc 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/distance.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/distance.h @@ -16,7 +16,7 @@ // $Id$ // // -// Author(s) : Sebastien Loriot +// Author(s) : Maxime Gimeno and Sebastien Loriot #ifndef CGAL_POLYGON_MESH_PROCESSING_DISTANCE_H #define CGAL_POLYGON_MESH_PROCESSING_DISTANCE_H @@ -26,28 +26,42 @@ #include #include #include +#include #include -#include +#include +#include +#include -/// \cond SKIP_IN_MANUAL +#include +#include + +#include +#ifdef CGAL_LINKED_WITH_TBB +#include +#include +#include +#endif // CGAL_LINKED_WITH_TBB + +#include +#include namespace CGAL{ namespace Polygon_mesh_processing { namespace internal{ - template OutputIterator -triangle_grid_sampling(const typename Kernel::Point_3& p0, - const typename Kernel::Point_3& p1, - const typename Kernel::Point_3& p2, - double distance, - OutputIterator out) +triangle_grid_sampling( const typename Kernel::Point_3& p0, + const typename Kernel::Point_3& p1, + const typename Kernel::Point_3& p2, + double distance, + OutputIterator out) { - const double d_p0p1 = CGAL::sqrt( CGAL::squared_distance(p0, p1) ); - const double d_p0p2 = CGAL::sqrt( CGAL::squared_distance(p0, p2) ); + typename Kernel::Compute_squared_distance_3 squared_distance; + const double d_p0p1 = to_double(approximate_sqrt( squared_distance(p0, p1) )); + const double d_p0p2 = to_double(approximate_sqrt( squared_distance(p0, p2) )); const double n = (std::max)(std::ceil( d_p0p1 / distance ), - std::ceil( d_p0p2 / distance )); + std::ceil( d_p0p2 / distance )); for (double i=1; i -OutputIterator -triangle_grid_sampling(const typename Kernel::Triangle_3& t, double distance, OutputIterator out) -{ - return triangle_grid_sampling(t[0], t[1], t[2], distance, out); -} - } //end of namespace internal -template +template OutputIterator -sample_triangles(const TriangleRange& triangles, double distance, OutputIterator out) +sample_triangles(const FaceRange& triangles, + const TriangleMesh& tm, + VertexPointMap vpm, + double distance, + OutputIterator out, + bool sample_faces, + bool sample_edges, + bool add_vertices) { - typedef typename Kernel::Point_3 Point_3; + typedef typename boost::property_traits::reference Point_ref; typedef typename Kernel::Vector_3 Vector_3; - typedef typename Kernel::Triangle_3 Triangle_3; + typedef boost::graph_traits GT; + typedef typename GT::face_descriptor face_descriptor; + typedef typename GT::halfedge_descriptor halfedge_descriptor; - std::set< std::pair > sampled_edges; + boost::unordered_set sampled_edges; + boost::unordered_set endpoints; - // sample edges but skip endpoints - BOOST_FOREACH(const Triangle_3& t, triangles) + BOOST_FOREACH(face_descriptor fd, triangles) { + // sample edges but skip endpoints + halfedge_descriptor hd = halfedge(fd, tm); for (int i=0;i<3; ++i) { - const Point_3& p0=t[i]; - const Point_3& p1=t[(i+1)%3]; - if ( sampled_edges.insert(CGAL::make_sorted_pair(p0, p1)).second ) + if (sample_edges && sampled_edges.insert(edge(hd, tm)).second ) { - const double d_p0p1 = CGAL::sqrt( CGAL::squared_distance(p0, p1) ); + Point_ref p0 = get(vpm, source(hd, tm)); + Point_ref p1 = get(vpm, target(hd, tm)); + typename Kernel::Compute_squared_distance_3 squared_distance; + const double d_p0p1 = to_double(approximate_sqrt( squared_distance(p0, p1) )); const double nb_pts = std::ceil( d_p0p1 / distance ); - const Vector_3 step_vec = (p1 - p0) / nb_pts; + const Vector_3 step_vec = typename Kernel::Construct_scaled_vector_3()( + typename Kernel::Construct_vector_3()(p0, p1), + typename Kernel::FT(1)/typename Kernel::FT(nb_pts)); for (double i=1; i(p0, p1, p2, distance, out); + } + } + return out; +} + +#ifdef CGAL_LINKED_WITH_TBB +template +struct Distance_computation{ + const AABB_tree& tree; + const std::vector& sample_points; + Point_3 initial_hint; + cpp11::atomic* distance; + + Distance_computation( + const AABB_tree& tree, + const Point_3& p, + const std::vector& sample_points, + cpp11::atomic* d) + : tree(tree) + , sample_points(sample_points) + , initial_hint(p) + , distance(d) + {} + + void + operator()(const tbb::blocked_range& range) const + { + Point_3 hint = initial_hint; + double hdist = 0; + for( std::size_t i = range.begin(); i != range.end(); ++i) + { + hint = tree.closest_point(sample_points[i], hint); + typename Kernel_traits::Kernel::Compute_squared_distance_3 squared_distance; + double d = to_double(CGAL::approximate_sqrt( squared_distance(hint,sample_points[i]) )); + if (d>hdist) hdist=d; + } + + if (hdist > distance->load(cpp11::memory_order_acquire)) + distance->store(hdist, cpp11::memory_order_release); + } +}; +#endif + +/** \ingroup PMP_distance_grp + * generates points taken on `tm` and outputs them to `out`, the sampling method + * is selected using named parameters. + * @tparam TriangleMesh a model of the concept `FaceListGraph` + * @tparam OutputIterator a model of `OutputIterator` + * holding objects of the same point type as + * the value type of the internal vertex point map of `tm` + * + * @param tm the triangle mesh that will be sampled + * @param out output iterator to be filled with sampled points + * @param np an optional sequence of \ref namedparameters among the ones listed below + * + * \cgalNamedParamsBegin + * \cgalParamBegin{vertex_point_map} the property map with the points + * associated to the vertices of `tm`. If this parameter is omitted, + * an internal property map for `CGAL::vertex_point_t` + * should be available for `TriangleMesh`. + * \cgalParamEnd + * \cgalParamBegin{geom_traits} a model of `PMPDistanceTraits`. \cgalParamEnd + * \cgalParamBegin{use_random_uniform_sampling} + * if `true` is passed (the default), points are generated in a random + * and uniform way on the surface of `tm`, and/or on edges of `tm`. + * For faces, the number of sample points is the value passed to the named + * parameter `number_of_points_on_faces()`. If not set, + * the value passed to the named parameter `number_of_points_per_area_unit()` + * is multiplied by the area of `tm` to get the number of sample points. + * If none of these parameters is set, the number of points sampled is `num_vertices(tm)`. + * For edges, the number of the number of sample points is the value passed to the named + * parameter `number_of_points_on_edges()`. If not set, + * the value passed to the named parameter `number_of_points_per_distance_unit()` + * is multiplied by the sum of the length of edges of `tm` to get the number of sample points. + * If none of these parameters is set, the number of points sampled is `num_vertices(tm)`. + * \cgalParamEnd + * \cgalParamBegin{use_grid_sampling} + * if `true` is passed, points are generated on a grid in each triangle, + * with a minimum of one point per triangle. The distance between + * two consecutive points in the grid is that of the length of the + * smallest non-null edge of `tm` or the value passed to + * the named parameter `grid_spacing()`. Edges are also split using the + * same distance, if requested. + * \cgalParamEnd + * \cgalParamBegin{use_monte_carlo_sampling} + * if `true` is passed, points are generated randomly in each triangle and/or + * on each edge. + * For faces, the number of points per triangle is the value passed to the named + * parameter `number_of_points_per_face()`. If not set, the value passed + * to the named parameter `number_of_points_per_area_unit()` is + * used to pick a number of points per face proportional to the triangle + * area with a minimum of one point per face. If none of these parameters + * is set, 2 divided by the square of the length of the smallest non-null + * edge of `tm` is used as if it was passed to + * `number_of_points_per_area_unit()`. + * For edges, the number of points per edge is the value passed to the named + * parameter `number_of_points_per_edge()`. If not set, the value passed + * to the named parameter `number_of_points_per_distance_unit()` is + * used to pick a number of points per edge proportional to the length of + * the edge with a minimum of one point per face. If none of these parameters + * is set, 1 divided by the length of the smallest non-null edge of `tm` + * is used as if it was passed to `number_of_points_per_distance_unit()`. + * \cgalParamEnd + * \cgalParamBegin{sample_vertices} if `true` is passed (default value), + * vertices of `tm` are put into `out`.\cgalParamEnd + * \cgalParamBegin{sample_edges} if `true` is passed (default value), + * edges of `tm` are sampled.\cgalParamEnd + * \cgalParamBegin{sample_faces} if `true` is passed (default value), + * faces of `tm` are sampled.\cgalParamEnd + * \cgalParamBegin{grid_spacing} a double value used as the grid spacing + * for the grid sampling method. + * \cgalParamEnd + * \cgalParamBegin{number_of_points_on_edges} an unsigned integral value used + * for the random sampling method as the number of points to pick exclusively + * on edges. + * \cgalParamEnd + * \cgalParamBegin{number_of_points_on_faces} an unsigned integral value used + * for the random sampling method as the number of points to pick on the surface. + * \cgalParamEnd + * \cgalParamBegin{number_of_points_per_distance_unit} a double value + * used for the random sampling and the Monte Carlo sampling methods to + * repectively determine the total number of points on edges and the + * number of points per edge. + * \cgalParamEnd + * \cgalParamBegin{number_of_points_per_edge} an unsigned integral value + * used by the Monte-Carlo sampling method as the number of points per edge + * to pick. + * \cgalParamEnd + * \cgalParamBegin{number_of_points_per_face} an unsigned integral value + * used by the Monte-Carlo sampling method as the number of points per face + * to pick. + * \cgalParamEnd + * \cgalParamBegin{number_of_points_per_area_unit} a double value + * used for the random sampling and the Monte Carlo sampling methods to + * repectively determine the total number of points inside faces + * and the number of points per face. + * \cgalParamEnd + * \cgalNamedParamsEnd + */ +template +OutputIterator +sample_triangle_mesh(const TriangleMesh& tm, + OutputIterator out, + NamedParameters np) +{ + typedef typename GetGeomTraits::type Geom_traits; + + typedef typename GetVertexPointMap::const_type Vpm; + + typedef boost::graph_traits GT; + typedef typename GT::face_descriptor face_descriptor; + typedef typename GT::halfedge_descriptor halfedge_descriptor; + typedef typename GT::edge_descriptor edge_descriptor; + + using boost::choose_param; + using boost::get_param; + using boost::is_default_param; + + Vpm pmap = choose_param(get_param(np, vertex_point), + get_const_property_map(vertex_point, tm)); + typedef Creator_uniform_3 Creator; + + Geom_traits geomtraits = choose_param(get_param(np, geom_traits), Geom_traits()); + + + bool use_rs = choose_param(get_param(np, random_uniform_sampling), true); + bool use_gs = choose_param(get_param(np, grid_sampling), false); + bool use_ms = choose_param(get_param(np, monte_carlo_sampling), false); + + if (use_gs || use_ms) + if (is_default_param(get_param(np, random_uniform_sampling))) + use_rs=false; + + bool smpl_vrtcs = choose_param(get_param(np, do_sample_vertices), true); + bool smpl_dgs = choose_param(get_param(np, do_sample_edges), true); + bool smpl_fcs = choose_param(get_param(np, do_sample_faces), true); + + double nb_pts_a_u = choose_param(get_param(np, nb_points_per_area_unit), 0.); + double nb_pts_l_u = choose_param(get_param(np, nb_points_per_distance_unit), 0.); + + // sample vertices + if (smpl_vrtcs) + { + Property_map_to_unary_function unary(pmap); + out = std::copy( + boost::make_transform_iterator(boost::begin(vertices(tm)), unary), + boost::make_transform_iterator(boost::end(vertices(tm)), unary), + out); + } + + // grid sampling + if (use_gs) + { + double grid_spacing_ = choose_param(get_param(np, grid_spacing), 0.); + if (grid_spacing_==0.) + { + // set grid spacing to the shortest edge length + double grid_spacing_ = (std::numeric_limits::max)(); + typedef typename boost::graph_traits + ::edge_descriptor edge_descriptor; + BOOST_FOREACH(edge_descriptor ed, edges(tm)) + { + double el = std::sqrt( + to_double( typename Geom_traits::Compute_squared_distance_3()( + get(pmap, source(ed, tm)), get(pmap, target(ed, tm)) ))); + if (el > 0 && el < grid_spacing_) + grid_spacing_ = el; + } + } + out=sample_triangles( + faces(tm), tm, pmap, grid_spacing_, out,smpl_fcs, smpl_dgs, false); + } + + // monte carlo sampling + if (use_ms) + { + typename Geom_traits::Compute_squared_distance_3 squared_distance; + double min_edge_length = (std::numeric_limits::max)(); + + std::size_t nb_points_per_face = + choose_param(get_param(np, number_of_points_per_face), 0); + std::size_t nb_points_per_edge = + choose_param(get_param(np, number_of_points_per_edge), 0); + + if ((nb_points_per_face == 0 && nb_pts_a_u ==0.) || + (nb_points_per_edge == 0 && nb_pts_l_u ==0.) ) + { + typedef typename boost::graph_traits + ::edge_descriptor edge_descriptor; + BOOST_FOREACH(edge_descriptor ed, edges(tm)) + { + double el = std::sqrt( + to_double( squared_distance(get(pmap, source(ed, tm)), + get(pmap, target(ed, tm)) ))); + if (min_edge_length > 0 && el < min_edge_length) + min_edge_length = el; + } + } + + // sample faces + if (smpl_fcs) + { + // set default value + if (nb_points_per_face == 0 && nb_pts_a_u ==0.) + nb_pts_a_u = 2. / CGAL::square(min_edge_length); + + BOOST_FOREACH(face_descriptor f, faces(tm)) + { + std::size_t nb_points = nb_points_per_face; + if (nb_points == 0) + { + nb_points = (std::max)( + static_cast( + std::ceil(to_double( + face_area(f,tm,parameters::geom_traits(geomtraits)))*nb_pts_a_u)) + ,std::size_t(1)); + } + // extract triangle face points + typename Geom_traits::Point_3 points[3]; + halfedge_descriptor hd(halfedge(f,tm)); + for(int i=0; i<3; ++i) + { + points[i] = get(pmap, target(hd, tm)); + hd = next(hd, tm); + } + // sample the triangle face + Random_points_in_triangle_3 + g(points[0], points[1], points[2]); + out=CGAL::cpp11::copy_n(g, nb_points, out); + } + } + // sample edges + if (smpl_dgs) + { + if (nb_points_per_edge == 0 && nb_pts_l_u == 0) + nb_pts_l_u = 1. / min_edge_length; + BOOST_FOREACH(edge_descriptor ed, edges(tm)) + { + std::size_t nb_points = nb_points_per_edge; + if (nb_points == 0) + { + nb_points = (std::max)( + static_cast( std::ceil( std::sqrt( to_double( + squared_distance(get(pmap, source(ed, tm)), + get(pmap, target(ed, tm)) )) )*nb_pts_l_u ) ), + std::size_t(1)); + } + // now do the sampling of the edge + Random_points_on_segment_3 + g(get(pmap, source(ed,tm)), get(pmap, target(ed,tm))); + out=CGAL::cpp11::copy_n(g, nb_points, out); + } } } - // sample triangles - BOOST_FOREACH(const Triangle_3& t, triangles) - out=internal::triangle_grid_sampling(t, distance, out); - - //add endpoints - std::set< Point_3 > endpoints; - BOOST_FOREACH(const Triangle_3& t, triangles) - for(int i=0; i<3; ++i) + // random uniform sampling + if (use_rs) + { + // sample faces + if(smpl_fcs) { - if ( endpoints.insert(t[i]).second ) *out++=t[i]; + std::size_t nb_points = choose_param(get_param(np, number_of_points_on_faces), 0); + Random_points_in_triangle_mesh_3 g(tm, pmap); + if (nb_points == 0) + { + if (nb_pts_a_u == 0.) + nb_points = num_vertices(tm); + else + nb_points = static_cast( + std::ceil(g.mesh_area()*nb_pts_a_u) ); + } + out = CGAL::cpp11::copy_n(g, nb_points, out); } + // sample edges + if (smpl_dgs) + { + std::size_t nb_points = + choose_param(get_param(np, number_of_points_on_edges), 0); + Random_points_on_edge_list_graph_3 g(tm, pmap); + if (nb_points == 0) + { + if (nb_pts_l_u == 0) + nb_points = num_vertices(tm); + else + nb_points = static_cast( + std::ceil( g.mesh_length()*nb_pts_a_u) ); + } + out = CGAL::cpp11::copy_n(g, nb_points, out); + } + } + return out; } -template -double approximated_Hausdorff_distance( - const TriangleRange1& triangles_1, - const TriangleRange2& triangles_2, - double targeted_precision -) +template +OutputIterator +sample_triangle_mesh(const TriangleMesh& tm, + OutputIterator out) { - std::vector sample_points; - sample_triangles( triangles_1, - targeted_precision, - std::back_inserter(sample_points) ); - // std::ofstream out("/tmp/samples.xyz"); - // std::copy(sample_points.begin(), sample_points.end(), std::ostream_iterator(out,"\n")); - - typedef typename TriangleRange2::const_iterator Triangle_iterator; - typedef AABB_triangle_primitive Primitive; - typedef AABB_traits Traits; - typedef AABB_tree< Traits > Tree; - - Tree tree(triangles_2.begin(), triangles_2.end()); - tree.accelerate_distance_queries(); - - double hdist = 0; - BOOST_FOREACH(const typename Kernel::Point_3& pt, sample_points) - { - double d = CGAL::sqrt( tree.squared_distance(pt) ); - // if ( d > 1e-1 ) std::cout << pt << "\n"; - if (d>hdist) hdist=d; - } - - return hdist; + return sample_triangle_mesh(tm, out, parameters::all_default()); } -template -double approximated_symmetric_Hausdorff_distance( - const TriangleRange1& triangles_1, - const TriangleRange2& triangles_2, - double targeted_precision -) +template +double approximate_Hausdorff_distance( + const PointRange& original_sample_points, + const TriangleMesh& tm, + VertexPointMap vpm) +{ + CGAL_assertion_code( bool is_triangle = is_triangle_mesh(tm) ); + CGAL_assertion_msg (is_triangle, + "Mesh is not triangulated. Distance computing impossible."); + #ifdef CGAL_HAUSDORFF_DEBUG + std::cout << "Nb sample points " << sample_points.size() << "\n"; + #endif + typedef typename Kernel::Point_3 Point_3; + std::vector sample_points + (boost::begin(original_sample_points), boost::end(original_sample_points) ); + + spatial_sort(sample_points.begin(), sample_points.end()); + + typedef AABB_face_graph_triangle_primitive Primitive; + typedef AABB_tree< AABB_traits > Tree; + + Tree tree( faces(tm).first, faces(tm).second, tm); + tree.accelerate_distance_queries(); + tree.build(); + Point_3 hint = get(vpm, *vertices(tm).first); +#ifndef CGAL_LINKED_WITH_TBB + CGAL_static_assertion_msg (!(boost::is_convertible::value), + "Parallel_tag is enabled but TBB is unavailable."); +#else + if (boost::is_convertible::value) + { + cpp11::atomic distance(0); + Distance_computation f(tree, hint, sample_points, &distance); + tbb::parallel_for(tbb::blocked_range(0, sample_points.size()), f); + return distance; + } + else +#endif + { + double hdist = 0; + BOOST_FOREACH(const Point_3& pt, sample_points) + { + hint = tree.closest_point(pt, hint); + typename Kernel::Compute_squared_distance_3 squared_distance; + typename Kernel::FT dist = squared_distance(hint,pt); + double d = to_double(CGAL::approximate_sqrt(dist)); + if(d>hdist) + hdist=d; + } + return hdist; + } +} + +template +double approximate_Hausdorff_distance( + const TriangleMesh& tm1, + const TriangleMesh& tm2, + NamedParameters np, + VertexPointMap vpm_2) +{ + std::vector sample_points; + sample_triangle_mesh( + tm1, + std::back_inserter(sample_points), + np); + return approximate_Hausdorff_distance(sample_points, tm2, vpm_2); +} + +// documented functions + +/** + * \ingroup PMP_distance_grp + * computes the approximate Hausdorff distance from `tm1` to `tm2` by returning + * the distance of the farthest point from `tm2` amongst a sampling of `tm1` + * generated with the function `sample_triangle_mesh()` with + * `tm1` and `np1` as parameter. + * + * A parallel version is provided and requires the executable to be + * linked against the Intel TBB library. + * To control the number of threads used, the user may use the `tbb::task_scheduler_init` class. + * See the TBB documentation + * for more details. + * + * @tparam Concurrency_tag enables sequential versus parallel algorithm. + * Possible values are `Sequential_tag` + * and `Parallel_tag`. + * @tparam TriangleMesh a model of the concept `FaceListGraph` + * @tparam NamedParameters1 a sequence of \ref namedparameters for `tm1` + * @tparam NamedParameters2 a sequence of \ref namedparameters for `tm2` + * + * @param tm1 the triangle mesh that will be sampled + * @param tm2 the triangle mesh to compute the distance to + * @param np1 optional sequence of \ref namedparameters for `tm1` passed to `sample_triangle_mesh()`. + * + * @param np2 optional sequence of \ref namedparameters for `tm2` among the ones listed below + * + * \cgalNamedParamsBegin + * \cgalParamBegin{vertex_point_map} the property map with the points associated to the vertices of `tm2` + * If this parameter is omitted, an internal property map for CGAL::vertex_point_t should be available in `TriangleMesh` + * and in all places where vertex_point_map is used. + * \cgalParamEnd + * \cgalNamedParamsEnd + * The function `CGAL::Polygon_mesh_processing::params::all_default()` can be used to indicate to use the default values for + * `np1` and specify custom values for `np2` + */ +template< class Concurrency_tag, + class TriangleMesh, + class NamedParameters1, + class NamedParameters2> +double approximate_Hausdorff_distance( const TriangleMesh& tm1, + const TriangleMesh& tm2, + const NamedParameters1& np1, + const NamedParameters2& np2) +{ + typedef typename GetGeomTraits::type Geom_traits; + + return approximate_Hausdorff_distance( + tm1, tm2, np1, choose_param(get_param(np2, vertex_point), + get_const_property_map(vertex_point, tm2))); +} + +/** + * \ingroup PMP_distance_grp + * computes the approximate symmetric Hausdorff distance between `tm1` and `tm2`. + * It returns the maximum of `approximate_Hausdorff_distance(tm1, tm2, np1, np2)` + * and `approximate_Hausdorff_distance(tm2, tm1, np2, np1)`. + */ +template< class Concurrency_tag, + class TriangleMesh, + class NamedParameters1, + class NamedParameters2> +double approximate_symmetric_Hausdorff_distance( + const TriangleMesh& tm1, + const TriangleMesh& tm2, + const NamedParameters1& np1, + const NamedParameters2& np2) { return (std::max)( - approximated_Hausdorff_distance(triangles_1, triangles_2, targeted_precision), - approximated_Hausdorff_distance(triangles_2, triangles_1, targeted_precision) + approximate_Hausdorff_distance(tm1,tm2,np1,np2), + approximate_Hausdorff_distance(tm2,tm1,np2,np1) ); } +/** + * \ingroup PMP_distance_grp + * returns the distance to `tm` of the point from `points` + * that is the furthest from `tm`. + * @tparam PointRange a range of `Point_3`, model of `Range`. + * @tparam TriangleMesh a model of the concept `FaceListGraph` + * @tparam NamedParameters a sequence of \ref namedparameters + * @param points the range of points of interest + * @param tm the triangle mesh to compute the distance to + * @param np an optional sequence of \ref namedparameters among the ones listed below + * + * \cgalNamedParamsBegin + * \cgalParamBegin{vertex_point_map} + * the property map with the points associated to the vertices of `tm`. If this parameter is omitted, + * an internal property map for `CGAL::vertex_point_t` should be available for the + vertices of `tm` \cgalParamEnd + * \cgalParamBegin{geom_traits} an instance of a geometric traits class, model of `PMPDistanceTraits`\cgalParamEnd + * \cgalNamedParamsEnd + */ +template< class Concurrency_tag, + class TriangleMesh, + class PointRange, + class NamedParameters> +double max_distance_to_triangle_mesh(const PointRange& points, + const TriangleMesh& tm, + const NamedParameters& np) +{ + typedef typename GetGeomTraits::type Geom_traits; -///\todo add approximated_Hausdorff_distance(FaceGraph, FaceGraph) -///\todo add approximated_Hausdorff_distance(TriangleRange, FaceGraph) -///\todo add approximated_Hausdorff_distance(FaceGraph, TriangleRange) -///\todo add approximated_Hausdorff_distance(PointRange, FaceGraph) -///\todo add approximated_Hausdorff_distance(PointRange, TriangleRange) -///\todo add barycentric_triangle_sampling(Triangle, PointPerAreaUnit) -///\todo add random_triangle_sampling(Triangle, PointPerAreaUnit) -/// \todo maybe we should use a sampling using epec to avoid being too far from the sampled triangle - + return approximate_Hausdorff_distance + (points,tm,choose_param(get_param(np, vertex_point), + get_const_property_map(vertex_point, tm))); +} + +/*! + *\ingroup PMP_distance_grp + * returns an approximation of the distance between `points` and the point lying on `tm` that is the farthest from `points` + * @tparam PointRange a range of `Point_3`, model of `Range`. + * @tparam TriangleMesh a model of the concept `FaceListGraph` + * @tparam NamedParameters a sequence of \ref namedparameters + * @param tm a triangle mesh + * @param points a range of points + * @param precision for each triangle of `tm`, the distance of its farthest point from `points` is bounded. + * A triangle is subdivided into sub-triangles so that the difference of its distance bounds + * is smaller than `precision`. `precision` must be strictly positive to avoid infinite loops. + * @param np an optional sequence of \ref namedparameters among the ones listed below + * + * \cgalNamedParamsBegin + * \cgalParamBegin{vertex_point_map} + * the property map with the points associated to the vertices of `tm`. If this parameter is omitted, + * an internal property map for `CGAL::vertex_point_t` should be available for the + vertices of `tm` \cgalParamEnd + * \cgalParamBegin{geom_traits} an instance of a geometric traits class, model of `PMPDistanceTraits`. \cgalParamEnd + * \cgalNamedParamsEnd + */ +template< class TriangleMesh, + class PointRange, + class NamedParameters> +double approximate_max_distance_to_point_set(const TriangleMesh& tm, + const PointRange& points, + const double precision, + const NamedParameters& np) +{ + typedef typename GetGeomTraits::type Geom_traits; + typedef boost::graph_traits GT; + + typedef Orthogonal_k_neighbor_search > Knn; + typedef typename Knn::Tree Tree; + Tree tree(points.begin(), points.end()); + CRefiner ref; + BOOST_FOREACH(typename GT::face_descriptor f, faces(tm)) + { + typename Geom_traits::Point_3 points[3]; + typename GT::halfedge_descriptor hd(halfedge(f,tm)); + for(int i=0; i<3; ++i) + { + points[i] = get(choose_param(get_param(np, vertex_point), + get_const_property_map(vertex_point, tm)), + target(hd, tm)); + hd = next(hd, tm); + } + ref.add(points[0], points[1], points[2], tree); + } + return to_double(ref.refine(precision, tree)); +} + +// convenience functions with default parameters + +template< class Concurrency_tag, + class TriangleMesh, + class PointRange> +double max_distance_to_triangle_mesh(const PointRange& points, + const TriangleMesh& tm) +{ + return max_distance_to_triangle_mesh + (points, tm, parameters::all_default()); +} + +template< class TriangleMesh, + class PointRange> +double approximate_max_distance_to_point_set(const TriangleMesh& tm, + const PointRange& points, + const double precision) +{ + return approximate_max_distance_to_point_set(tm, points, precision, + parameters::all_default()); +} + +template< class Concurrency_tag, + class TriangleMesh, + class NamedParameters> +double approximate_Hausdorff_distance(const TriangleMesh& tm1, + const TriangleMesh& tm2, + const NamedParameters& np) +{ + return approximate_Hausdorff_distance( + tm1, tm2, np, parameters::all_default()); +} + +template< class Concurrency_tag, + class TriangleMesh> +double approximate_Hausdorff_distance(const TriangleMesh& tm1, + const TriangleMesh& tm2) +{ + return approximate_Hausdorff_distance( + tm1, tm2, parameters::all_default(), parameters::all_default()); +} +template< class Concurrency_tag, + class TriangleMesh, + class NamedParameters> +double approximate_symmetric_Hausdorff_distance(const TriangleMesh& tm1, + const TriangleMesh& tm2, + const NamedParameters& np) +{ + return approximate_symmetric_Hausdorff_distance( + tm1, tm2, np, parameters::all_default()); +} -} } // end of namespace CGAL::Polygon_mesh_processing +template< class Concurrency_tag, + class TriangleMesh> +double approximate_symmetric_Hausdorff_distance(const TriangleMesh& tm1, + const TriangleMesh& tm2) +{ + return approximate_symmetric_Hausdorff_distance( + tm1, tm2, parameters::all_default(), parameters::all_default()); +} + +} +} // end of namespace CGAL::Polygon_mesh_processing -/// \endcond #endif //CGAL_POLYGON_MESH_PROCESSING_DISTANCE_H diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/mesh_to_point_set_hausdorff_distance.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/mesh_to_point_set_hausdorff_distance.h new file mode 100644 index 00000000000..9ec3c3d553c --- /dev/null +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/mesh_to_point_set_hausdorff_distance.h @@ -0,0 +1,341 @@ +// Copyright (c) 2016 INRIA Sophia-Antipolis (France). +// All rights reserved. +// +// This file is part of CGAL (www.cgal.org). +// You can redistribute it and/or modify it under the terms of the GNU +// General Public License as published by the Free Software Foundation, +// either version 3 of the License, or (at your option) any later version. +// +// Licensees holding a valid commercial license may use this file in +// accordance with the commercial license agreement provided with the software. +// +// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE +// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. +// +// $URL$ +// $Id$ +// +// +// Author(s) : Simon Giraudot and Maxime Gimeno + +#ifndef CGAL_MESH_TO_POINT_SET_HAUSDORFF_DISTANCE_H +#define CGAL_MESH_TO_POINT_SET_HAUSDORFF_DISTANCE_H + +#include +#include +#include +#include +#include + +namespace CGAL{ +namespace internal{ +template +class CPointH +{ + +public: + typedef typename Kernel::Point_3 Point; + typedef typename Kernel::FT FT; + + CPointH () + { + m_point = Point (FT(0.), FT(0.), FT(0.)); + m_hausdorff = 0.; + } + CPointH (const Point& p, FT h = 0.) + { + m_point = p; + m_hausdorff = h; + } + CPointH (const CPointH& p) + { + m_point = p (); + m_hausdorff = p.hausdorff (); + } + + const Point& operator() () const { return m_point; } + Point& operator() () { return m_point; } + FT hausdorff () const { return m_hausdorff; } + + static Point mid_point (const CPointH& a, const CPointH& b) + { + return Point (FT(0.5) * (a().x() + b().x()), + FT(0.5) * (a().y() + b().y()), + FT(0.5) * (a().z() + b().z())); + } + + +private: + + Point m_point; + FT m_hausdorff; +}; + +template +class CRefTriangle +{ + +private: + + typedef typename Kernel::FT FT; + typedef CPointH PointH; + typedef typename PointH::Point Point; + PointH m_point[3]; + std::size_t m_edge; + FT m_upper_bound; + FT m_lower_bound; + Point m_bisector; + +public: + + CRefTriangle (const PointH& a, const PointH& b, const PointH& c) // Triangle + { + m_point[0] = a; + m_point[1] = b; + m_point[2] = c; + + m_edge = 0; + + typename Kernel::Compute_squared_distance_3 squared_distance; + FT length_max = squared_distance(m_point[1](), m_point[2]()); + FT length1 = squared_distance(m_point[2](), m_point[0]()); + if (length1 > length_max) + { + m_edge = 1; + length_max = length1; + } + FT length2 = squared_distance(m_point[0](), m_point[1]()); + if (length2 > length_max) + m_edge = 2; + + m_bisector = PointH::mid_point (m_point[(m_edge+1)%3], + m_point[(m_edge+2)%3]); + + m_lower_bound = 0.; + m_upper_bound = 0.; + for (unsigned int i = 0; i < 3; ++ i) + { + if (m_point[i].hausdorff () > m_lower_bound) + m_lower_bound = m_point[i].hausdorff (); + + FT up = m_point[i].hausdorff () + + CGAL::approximate_sqrt (squared_distance (m_point[i](), m_bisector)); + + if (up > m_upper_bound) + m_upper_bound = up; + } + } + + CRefTriangle (const CRefTriangle& t) + { + m_point[0] = t.points ()[0]; + m_point[1] = t.points ()[1]; + m_point[2] = t.points ()[2]; + m_edge = t.edge (); + m_lower_bound = t.lower_bound (); + m_upper_bound = t.upper_bound (); + m_bisector = t.bisector (); + } + + FT lower_bound () const + { + return m_lower_bound; + } + + FT upper_bound () const + { + return m_upper_bound; + } + + const Point& bisector () const + { + return m_bisector; + } + + friend bool operator< (const CRefTriangle& a, const CRefTriangle& b) + { + return a.upper_bound () < b.upper_bound (); + } + + const PointH* points () const { return m_point; } + std::size_t edge () const { return m_edge; } + + #ifdef CGAL_MTPS_HD_DEBUG + void print () const + { + std::cerr << "[Refinement triangle]" << std::endl + << " Bounds: " << m_lower_bound << " to " + << m_upper_bound << std::endl + << " * " << m_point[0]() << std::endl + << " * " << m_point[1]() << std::endl + << " * " << m_point[2]() << std::endl + << " -> " << m_bisector << std::endl; + } + #endif +}; +}//internal + +template +class CRefiner +{ +private: + + typedef typename Kernel::Point_3 Point; + typedef typename Kernel::Triangle_3 Triangle; + typedef typename Kernel::FT FT; + typedef internal::CPointH PointH; + typedef internal::CRefTriangle RefTriangle; + typedef CGAL::Orthogonal_k_neighbor_search > Knn; + typedef typename Knn::Tree Tree; + std::priority_queue m_queue; + FT m_lower_bound; + FT m_upper_bound; + Point m_point; + +public: + + CRefiner () + { + m_lower_bound = 0.; + m_upper_bound = (std::numeric_limits::max) (); + } + + bool empty () + { + return m_queue.empty (); + } + void reset (FT lower_bound = 0.) + { + m_queue = std::priority_queue (); + m_lower_bound = lower_bound; + m_upper_bound = (std::numeric_limits::max) (); + } + + inline FT uncertainty () const + { + return m_upper_bound - m_lower_bound; + } + + inline FT mid () const + { + return ((m_upper_bound + m_lower_bound) / 2.); + } + + inline FT lower_bound () const + { + return m_lower_bound; + } + + inline FT upper_bound () const + { + return m_upper_bound; + } + + inline FT relative_error () const + { + return 0.5 * uncertainty () / mid (); + } + + + const Point& hausdorff_point () const { return m_point; } + + bool add (const Point& a, const Point& b, const Point& c, const Tree& tree) + { + + RefTriangle r (PointH (a, CGAL::approximate_sqrt (Knn(tree, a, 1).begin()->second)), + PointH (b, CGAL::approximate_sqrt (Knn(tree, b, 1).begin()->second)), + PointH (c, CGAL::approximate_sqrt (Knn(tree, c, 1).begin()->second))); + + if (r.lower_bound () > m_lower_bound) + { + m_lower_bound = r.lower_bound (); + } + if (r.upper_bound () > m_lower_bound) + m_queue.push (r); + return true; + } + bool clean_up_queue () + { + std::size_t before = m_queue.size (); + std::vector to_keep; + while (!(m_queue.empty ())) + { + const RefTriangle& current = m_queue.top (); + if (current.upper_bound () > m_lower_bound) + to_keep.push_back (current); + m_queue.pop (); + } + + m_queue = std::priority_queue (); + BOOST_FOREACH(RefTriangle& r, to_keep) + m_queue.push (r); + + return (m_queue.size () < before); + } + + FT refine (FT limit, const Tree& tree, FT upper_bound = 1e30) + { + + unsigned int nb_clean = 0; + + while (uncertainty () > limit && !(m_queue.empty ())) + { + if (m_queue.size () > 100000) + { + #ifdef CGAL_MTPS_HD_DEBUG + m_queue.top ().print (); + #endif + ++ nb_clean; + if (nb_clean > 5) + return m_upper_bound; + if (!clean_up_queue ()) + return m_upper_bound; + } + + + const RefTriangle& current = m_queue.top (); + + m_upper_bound = current.upper_bound (); + + if (current.lower_bound () > m_lower_bound) + m_lower_bound = current.lower_bound (); + typename Kernel::Compute_squared_area_3 squared_area; + + if(squared_area(current.points()[0](), + current.points()[1](), + current.points()[2]() + ) < 1e-20) + { + m_queue.pop(); + continue; + } + const Point& bisector = current.bisector (); + m_point = bisector; + //squared distance between bisector and its closst point in the mesh + FT hausdorff = CGAL::approximate_sqrt (Knn(tree, bisector, 1).begin()->second); + if (hausdorff > m_lower_bound) + m_lower_bound = hausdorff; + + if (m_lower_bound > upper_bound) + return m_upper_bound; + + PointH new_point (bisector, hausdorff); + std::size_t i = current.edge (); + PointH p0 (current.points()[i]); + PointH p1 (current.points()[(i+1)%3]); + PointH p2 (current.points()[(i+2)%3]); + + m_queue.pop (); + + m_queue.push (RefTriangle (new_point, p0, p1)); + m_queue.push (RefTriangle (new_point, p0, p2)); + } + + + return m_upper_bound; + + } + +}; +}//CGAL +#endif // CGAL_MESH_TO_POINT_SET_HAUSDORFF_DISTANCE_H diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/named_function_params.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/named_function_params.h index 6695d9eb265..75c1f592143 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/named_function_params.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/named_function_params.h @@ -14,7 +14,7 @@ // // $URL$ // $Id$ -// +// // // Author(s) : Jane Tournois @@ -39,6 +39,19 @@ namespace CGAL{ enum relax_constraints_t { relax_constraints }; enum vertex_is_constrained_t { vertex_is_constrained }; enum face_patch_t { face_patch }; + enum random_uniform_sampling_t { random_uniform_sampling }; + enum grid_sampling_t { grid_sampling }; + enum monte_carlo_sampling_t { monte_carlo_sampling }; + enum do_sample_edges_t { do_sample_edges }; + enum do_sample_vertices_t { do_sample_vertices }; + enum do_sample_faces_t { do_sample_faces }; + enum number_of_points_on_faces_t { number_of_points_on_faces }; + enum number_of_points_per_face_t { number_of_points_per_face }; + enum grid_spacing_t { grid_spacing }; + enum nb_points_per_area_unit_t { nb_points_per_area_unit }; + enum number_of_points_per_edge_t { number_of_points_per_edge }; + enum number_of_points_on_edges_t { number_of_points_on_edges }; + enum nb_points_per_distance_unit_t{ nb_points_per_distance_unit }; //to be documented enum face_normal_t { face_normal }; @@ -216,6 +229,109 @@ namespace CGAL{ return Params(p, *this); } + template + pmp_bgl_named_params + use_random_uniform_sampling(const Boolean& p) const + { + typedef pmp_bgl_named_params Params; + return Params(p, *this); + } + + template + pmp_bgl_named_params + use_grid_sampling(const Boolean& p) const + { + typedef pmp_bgl_named_params Params; + return Params(p, *this); + } + + template + pmp_bgl_named_params + use_monte_carlo_sampling(const Boolean& p) const + { + typedef pmp_bgl_named_params Params; + return Params(p, *this); + } + + template + pmp_bgl_named_params + sample_edges(const Boolean& p) const + { + typedef pmp_bgl_named_params Params; + return Params(p, *this); + } + + template + pmp_bgl_named_params + sample_vertices(const Boolean& p) const + { + typedef pmp_bgl_named_params Params; + return Params(p, *this); + } + + template + pmp_bgl_named_params + sample_faces(const Boolean& p) const + { + typedef pmp_bgl_named_params Params; + return Params(p, *this); + } + + template + pmp_bgl_named_params + number_of_points_on_faces(const NT& n) const + { + typedef pmp_bgl_named_params Params; + return Params(n, *this); + } + + template + pmp_bgl_named_params + number_of_points_per_face(const NT& n) const + { + typedef pmp_bgl_named_params Params; + return Params(n, *this); + } + + template + pmp_bgl_named_params + grid_spacing(const NT& n) const + { + typedef pmp_bgl_named_params Params; + return Params(n, *this); + } + + template + pmp_bgl_named_params + number_of_points_per_area_unit(const NT& n) const + { + typedef pmp_bgl_named_params Params; + return Params(n, *this); + } + + template + pmp_bgl_named_params + number_of_points_per_edge(const NT& n) const + { + typedef pmp_bgl_named_params Params; + return Params(n, *this); + } + + template + pmp_bgl_named_params + number_of_points_on_edges(const NT& n) const + { + typedef pmp_bgl_named_params Params; + return Params(n, *this); + } + + template + pmp_bgl_named_params + number_of_points_per_distance_unit(const NT& n) const + { + typedef pmp_bgl_named_params Params; + return Params(n, *this); + } }; namespace Polygon_mesh_processing{ @@ -379,6 +495,123 @@ namespace parameters{ return Params(p); } + //overload + template + pmp_bgl_named_params + use_random_uniform_sampling(const Boolean& p) + { + typedef pmp_bgl_named_params Params; + return Params(p); + } + + //overload + template + pmp_bgl_named_params + use_grid_sampling(const Boolean& p) + { + typedef pmp_bgl_named_params Params; + return Params(p); + } + + //overload + template + pmp_bgl_named_params + use_monte_carlo_sampling(const Boolean& p) + { + typedef pmp_bgl_named_params Params; + return Params(p); + } + + //overload + template + pmp_bgl_named_params + sample_edges(const Boolean& p) + { + typedef pmp_bgl_named_params Params; + return Params(p); + } + + //overload + template + pmp_bgl_named_params + sample_vertices(const Boolean& p) + { + typedef pmp_bgl_named_params Params; + return Params(p); + } + + //overload + template + pmp_bgl_named_params + sample_faces(const Boolean& p) + { + typedef pmp_bgl_named_params Params; + return Params(p); + } + + //overload + template + pmp_bgl_named_params + number_of_points_on_faces(const NT& n) + { + typedef pmp_bgl_named_params Params; + return Params(n); + } + + //overload + template + pmp_bgl_named_params + number_of_points_per_face(const NT& n) + { + typedef pmp_bgl_named_params Params; + return Params(n); + } + + //overload + template + pmp_bgl_named_params + grid_spacing(const NT& n) + { + typedef pmp_bgl_named_params Params; + return Params(n); + } + + //overload + template + pmp_bgl_named_params + number_of_points_per_area_unit(const NT& n) + { + typedef pmp_bgl_named_params Params; + return Params(n); + } + + //overload + template + pmp_bgl_named_params + number_of_points_per_edge(const NT& n) + { + typedef pmp_bgl_named_params Params; + return Params(n); + } + + //overload + template + pmp_bgl_named_params + number_of_points_on_edges(const NT& n) + { + typedef pmp_bgl_named_params Params; + return Params(n); + } + + //overload + template + pmp_bgl_named_params + number_of_points_per_distance_unit(const NT& n) + { + typedef pmp_bgl_named_params Params; + return Params(n); + } + } //namespace parameters } //namespace Polygon_mesh_processing diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/measure.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/measure.h index 2ff93f14301..0fffe7af464 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/measure.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/measure.h @@ -14,7 +14,7 @@ // // $URL$ // $Id$ -// +// // // Author(s) : Andreas Fabri @@ -413,39 +413,43 @@ namespace Polygon_mesh_processing { #else typename GetGeomTraits::type::FT #endif - volume(const TriangleMesh& tmesh, const CGAL_PMP_NP_CLASS& np) +volume(const TriangleMesh& tmesh, const CGAL_PMP_NP_CLASS& np) +{ + CGAL_assertion(is_triangle_mesh(tmesh)); + CGAL_assertion(is_closed(tmesh)); + + using boost::choose_param; + using boost::get_param; + + typename GetVertexPointMap::const_type + vpm = choose_param(get_param(np, vertex_point), + get_const_property_map(CGAL::vertex_point, tmesh)); + typename GetGeomTraits::type::Point_3 + origin(0, 0, 0); + + typedef typename boost::graph_traits::face_descriptor face_descriptor; + + typename GetGeomTraits::type::FT volume = 0.; + typename CGAL::Kernel_traits::type>::Kernel::Compute_volume_3 cv3; + + BOOST_FOREACH(face_descriptor f, faces(tmesh)) { - CGAL_assertion(is_triangle_mesh(tmesh)); - CGAL_assertion(is_closed(tmesh)); - - using boost::choose_param; - using boost::get_param; - - typename GetVertexPointMap::const_type - vpm = choose_param(get_param(np, vertex_point), - get_const_property_map(CGAL::vertex_point, tmesh)); - typename GetGeomTraits::type::Point_3 - origin(0, 0, 0); - - typedef typename boost::graph_traits::face_descriptor face_descriptor; - - typename GetGeomTraits::type::FT volume = 0.; - BOOST_FOREACH(face_descriptor f, faces(tmesh)) - { - volume += CGAL::volume(origin, - get(vpm, target(halfedge(f, tmesh), tmesh)), - get(vpm, target(next(halfedge(f, tmesh), tmesh), tmesh)), - get(vpm, target(prev(halfedge(f, tmesh), tmesh), tmesh))); - exact(volume); - } - return volume; + volume += cv3(origin, + get(vpm, target(halfedge(f, tmesh), tmesh)), + get(vpm, target(next(halfedge(f, tmesh), tmesh), tmesh)), + get(vpm, target(prev(halfedge(f, tmesh), tmesh), tmesh))); + exact(volume); } + return volume; +} template typename CGAL::Kernel_traits::type>::Kernel::FT volume(const TriangleMesh& tmesh) { + return volume(tmesh, CGAL::Polygon_mesh_processing::parameters::all_default()); } diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/self_intersections.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/self_intersections.h index fe0e81f266c..8e4a68f970b 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/self_intersections.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/self_intersections.h @@ -241,7 +241,7 @@ self_intersections( const FaceRange& face_range, * \cgalParamBegin{vertex_point_map} the property map with the points associated to the vertices of `pmesh`. * If this parameter is omitted, an internal property map for * `CGAL::vertex_point_t` should be available in `TriangleMesh`\cgalParamEnd - * \cgalParamBegin{geom_traits} an instance of a geometric traits class, model of `SelfIntersectionTraits` \cgalParamEnd + * \cgalParamBegin{geom_traits} an instance of a geometric traits class, model of `PMPSelfIntersectionTraits` \cgalParamEnd * \cgalNamedParamsEnd * * @return `out` @@ -300,7 +300,7 @@ self_intersections(const TriangleMesh& tmesh, OutputIterator out) * \cgalParamBegin{vertex_point_map} the property map with the points associated to the vertices of `pmesh`. * If this parameter is omitted, an internal property map for * `CGAL::vertex_point_t` should be available in `TriangleMesh`\cgalParamEnd - * \cgalParamBegin{geom_traits} an instance of a geometric traits class, model of `SelfIntersectionTraits` \cgalParamEnd + * \cgalParamBegin{geom_traits} an instance of a geometric traits class, model of `PMPSelfIntersectionTraits` \cgalParamEnd * \cgalNamedParamsEnd */ @@ -386,7 +386,7 @@ OutputIterator self_intersections(const FaceRange& face_range, * \cgalParamBegin{vertex_point_map} the property map with the points associated to the vertices of `tmesh`. * If this parameter is omitted, an internal property map for * `CGAL::vertex_point_t` should be available in `TriangleMesh`\cgalParamEnd - * \cgalParamBegin{geom_traits} an instance of a geometric traits class, model of `SelfIntersectionTraits` \cgalParamEnd + * \cgalParamBegin{geom_traits} an instance of a geometric traits class, model of `PMPSelfIntersectionTraits` \cgalParamEnd * \cgalNamedParamsEnd * * @return true if `tmesh` self-intersects diff --git a/Polygon_mesh_processing/test/Polygon_mesh_processing/CMakeLists.txt b/Polygon_mesh_processing/test/Polygon_mesh_processing/CMakeLists.txt index 7eb4eaa9848..de04e69c3e3 100644 --- a/Polygon_mesh_processing/test/Polygon_mesh_processing/CMakeLists.txt +++ b/Polygon_mesh_processing/test/Polygon_mesh_processing/CMakeLists.txt @@ -49,16 +49,26 @@ endif() # include for local package include_directories( BEFORE ../../include ) -find_package(Eigen3 3.1.0) #(requires 3.1.0 or greater) +find_package(Eigen3 3.2.0) #(requires 3.2.0 or greater) +find_package( TBB ) +if( TBB_FOUND ) + include( ${TBB_USE_FILE} ) + list( APPEND CGAL_3RD_PARTY_LIBRARIES ${TBB_LIBRARIES} ) +else() + message( STATUS "NOTICE: Intel TBB was not found. test_pmp_distance will use sequential code." ) +endif() + +include( CGAL_CreateSingleSourceCGALProgram ) if (EIGEN3_FOUND) - # Executables that require Eigen 3.1 include( ${EIGEN3_USE_FILE} ) # Creating entries for all .cpp/.C files with "main" routine # ########################################################## - - include( CGAL_CreateSingleSourceCGALProgram ) + create_single_source_cgal_program("fairing_test.cpp") + create_single_source_cgal_program("triangulate_hole_Polyhedron_3_no_delaunay_test.cpp" ) + create_single_source_cgal_program("triangulate_hole_Polyhedron_3_test.cpp") +endif(EIGEN3_FOUND) create_single_source_cgal_program("connected_component_polyhedron.cpp") create_single_source_cgal_program("connected_component_surface_mesh.cpp") @@ -77,16 +87,5 @@ if (EIGEN3_FOUND) create_single_source_cgal_program("measures_test.cpp") create_single_source_cgal_program("triangulate_faces_test.cpp") create_single_source_cgal_program("test_pmp_remove_border_edge.cpp") - -if(NOT (${EIGEN3_VERSION} VERSION_LESS 3.2.0)) - create_single_source_cgal_program("fairing_test.cpp") - create_single_source_cgal_program("triangulate_hole_Polyhedron_3_no_delaunay_test.cpp" ) - create_single_source_cgal_program("triangulate_hole_Polyhedron_3_test.cpp") + create_single_source_cgal_program("test_pmp_distance.cpp") create_single_source_cgal_program("triangulate_hole_polyline_test.cpp") -else() - message(STATUS "NOTICE: Some tests require Eigen 3.2 (or greater) and will not be compiled.") -endif() - -else(EIGEN3_FOUND) - message(STATUS "NOTICE: Some examples require Eigen 3.1 (or greater) and will not be compiled.") -endif(EIGEN3_FOUND) diff --git a/Polygon_mesh_processing/test/Polygon_mesh_processing/test_pmp_distance.cmd b/Polygon_mesh_processing/test/Polygon_mesh_processing/test_pmp_distance.cmd new file mode 100644 index 00000000000..856bd25a5bf --- /dev/null +++ b/Polygon_mesh_processing/test/Polygon_mesh_processing/test_pmp_distance.cmd @@ -0,0 +1 @@ +data/elephant.off data/blobby_3cc.off diff --git a/Polygon_mesh_processing/test/Polygon_mesh_processing/test_pmp_distance.cpp b/Polygon_mesh_processing/test/Polygon_mesh_processing/test_pmp_distance.cpp new file mode 100644 index 00000000000..9064beb635a --- /dev/null +++ b/Polygon_mesh_processing/test/Polygon_mesh_processing/test_pmp_distance.cpp @@ -0,0 +1,309 @@ +#include +#include +#include +#include + + +#include + + +#include +#include + +struct Custom_traits_Hausdorff +{ +// New requirements { + struct FT + { + FT(){} + FT(double){} + FT operator/(FT)const{return FT();} + FT operator-(const FT)const{return FT();} + FT operator+(FT)const{return FT();} + FT operator*(FT)const{return FT();} + bool operator<=(FT)const{return false;} + bool operator>=(FT)const{return false;} + bool operator!=(FT)const{return false;} + bool operator<(FT)const{return false;} + bool operator>(FT)const{return false;} + FT& operator+=(FT){ return *this; } + }; + + struct Point_3 + { + Point_3(){} + Point_3(FT, FT, FT){} + FT operator[](int)const{return FT();} + bool operator<(Point_3)const{return false;} + double x()const{return 0;} + double y()const{return 0;} + double z()const{return 0;} + }; + + struct Vector_3{}; + + struct Triangle_3 + { + Triangle_3(const Point_3&, const Point_3&, const Point_3&){} + Point_3 operator[](int)const{return Point_3();} + CGAL::Bbox_3 bbox(){return CGAL::Bbox_3();} + }; + + struct Segment_3 + { + Segment_3(const Point_3&, const Point_3&){} + Point_3 operator[](int)const{return Point_3();} + }; + + struct Compute_squared_area_3 + { + typedef FT result_type; + FT operator()(Point_3,Point_3,Point_3)const{return FT();} + FT operator()(Triangle_3)const{return FT();} + }; + + struct Compute_squared_length_3 + { + typedef FT result_type; + FT operator()(Segment_3)const{return FT();} + }; + + struct Construct_translated_point_3 + { + Point_3 operator() (const Point_3 &, const Vector_3 &){return Point_3();} + }; + + struct Construct_vector_3{ + Vector_3 operator() (const Point_3 &, const Point_3 &){return Vector_3();} + }; + + struct Construct_scaled_vector_3 + { + Vector_3 operator() (const Vector_3 &, const FT &) + {return Vector_3();} + + }; + + Compute_squared_area_3 compute_squared_area_3_object(){return Compute_squared_area_3();} +// } end of new requirements + +// requirements from AABBGeomTraits { + struct Sphere_3 + {}; + + struct Equal_3 + { + bool operator()(Point_3,Point_3) const {return false;} + }; + + struct Do_intersect_3 + { + CGAL::Comparison_result operator()(const Sphere_3& , const CGAL::Bbox_3&){ return CGAL::ZERO;} + }; + + struct Intersect_3 + { + struct result_type{}; + }; + + struct Construct_sphere_3 + { + Sphere_3 operator()(const Point_3& , FT ){return Sphere_3();} + }; + + struct Construct_projected_point_3 + { + const Point_3 operator()(Triangle_3, Point_3)const{return Point_3();} + }; + + struct Compare_distance_3 + { + CGAL::Comparison_result operator()(Point_3, Point_3, Point_3) + {return CGAL::ZERO;} + }; + + struct Has_on_bounded_side_3 + { + //documented as Comparision_result + CGAL::Comparison_result operator()(const Point_3&, const Point_3&, const Point_3&) + {return CGAL::ZERO;} + bool operator()(const Sphere_3&, const Point_3&) + {return false;} + }; + + struct Compute_squared_radius_3 + { + FT operator()(const Sphere_3&){return FT();} + }; + + struct Compute_squared_distance_3 + { + FT operator()(const Point_3& , const Point_3& ){return FT();} + }; + + Compare_distance_3 compare_distance_3_object(){return Compare_distance_3();} + Construct_sphere_3 construct_sphere_3_object(){return Construct_sphere_3();} + Construct_projected_point_3 construct_projected_point_3_object(){return Construct_projected_point_3();} + Compute_squared_distance_3 compute_squared_distance_3_object(){return Compute_squared_distance_3();} + Do_intersect_3 do_intersect_3_object(){return Do_intersect_3();} + Equal_3 equal_3_object(){return Equal_3();} +// } end of requirments from AABBGeomTraits + + +// requirements from SearchGeomTraits_3 { + struct Iso_cuboid_3{}; + + struct Construct_iso_cuboid_3{}; + + struct Construct_min_vertex_3{}; + + struct Construct_max_vertex_3{}; + + struct Construct_center_3{}; + + + typedef const FT* Cartesian_const_iterator_3; + + + struct Construct_cartesian_const_iterator_3{ + Construct_cartesian_const_iterator_3(){} + Construct_cartesian_const_iterator_3(const Point_3&){} + const FT* operator()(const Point_3&) const + { return 0; } + + const FT* operator()(const Point_3&, int) const + { return 0; } + typedef const FT* result_type; + }; +// } end of requirements from SearchGeomTraits_3 + +// requirements from SpatialSortingTraits_3 { + struct Less_x_3{ + bool operator()(Point_3, Point_3){return false;} + }; + struct Less_y_3{ + bool operator()(Point_3, Point_3){return false;} + }; + struct Less_z_3{ + bool operator()(Point_3, Point_3){return false;} + }; + Less_x_3 less_x_3_object()const{return Less_x_3();} + Less_y_3 less_y_3_object()const{return Less_y_3();} + Less_z_3 less_z_3_object()const{return Less_z_3();} +// } end of requirements from SpatialSortingTraits_3 +}; + +namespace CGAL{ +template<>struct Kernel_traits +{ + typedef Custom_traits_Hausdorff Kernel; +}; + +template<> +struct Algebraic_structure_traits< Custom_traits_Hausdorff::FT > +: public Algebraic_structure_traits_base< Custom_traits_Hausdorff::FT, Field_tag > +{ +}; + + +template<> +struct Real_embeddable_traits< Custom_traits_Hausdorff::FT > + : public INTERN_RET::Real_embeddable_traits_base< Custom_traits_Hausdorff::FT , CGAL::Tag_true> +{ + class To_double + : public std::unary_function< Custom_traits_Hausdorff::FT, double > { + public: + double operator()( const Custom_traits_Hausdorff::FT& ) const { return 0; } + }; +}; + +}//end CGAL + + +void exact(Custom_traits_Hausdorff::FT){} + +typedef CGAL::Exact_predicates_inexact_constructions_kernel K; + +#include + +namespace PMP = CGAL::Polygon_mesh_processing; + + +typedef CGAL::Surface_mesh Mesh; + +template +void general_tests(const TriangleMesh& m1, + const TriangleMesh& m2 ) +{ + + if (m1.number_of_vertices()==0 || m2.number_of_vertices()==0) return; + + std::vector points; + points.push_back(typename GeomTraits::Point_3(0,0,0)); + points.push_back(typename GeomTraits::Point_3(0,1,0)); + points.push_back(typename GeomTraits::Point_3(1,0,0)); + points.push_back(typename GeomTraits::Point_3(0,0,1)); + points.push_back(typename GeomTraits::Point_3(0,2,0)); + + std::cout << "Symmetric distance between meshes (sequential) " + << PMP::approximate_symmetric_Hausdorff_distance( + m1,m2, + PMP::parameters::number_of_points_per_area_unit(4000), + PMP::parameters::number_of_points_per_area_unit(4000)) + << "\n"; + + std::cout << "Max distance to point set " + << PMP::approximate_max_distance_to_point_set(m1,points,4000) + << "\n"; + + std::cout << "Max distance to triangle mesh (sequential) " + << PMP::max_distance_to_triangle_mesh(points,m1) + << "\n"; +} + +void test_concept() +{ + typedef Custom_traits_Hausdorff CK; + CGAL::Surface_mesh m1, m2; + general_tests(m1,m2); +} + +int main(int, char** argv) +{ + Mesh m1,m2; + std::ifstream input(argv[1]); + input >> m1; + input.close(); + + input.open(argv[2]); + input >> m2; + + std::cout << "First mesh has " << num_faces(m1) << " faces\n"; + std::cout << "Second mesh has " << num_faces(m2) << " faces\n"; + + CGAL::Real_timer time; + #ifdef CGAL_LINKED_WITH_TBB + time.start(); + std::cout << "Distance between meshes (parallel) " + << PMP::approximate_Hausdorff_distance( + m1,m2,PMP::parameters::number_of_points_per_area_unit(4000)) + << "\n"; + time.stop(); + std::cout << "done in " << time.time() << "s.\n"; + #endif + + time.reset(); + time.start(); + std::cout << "Distance between meshes (sequential) " + << PMP::approximate_Hausdorff_distance( + m1,m2,PMP::parameters::number_of_points_per_area_unit(4000)) + << "\n"; + time.stop(); + std::cout << "done in " << time.time() << "s.\n"; + + general_tests(m1,m2); + + test_concept(); + + return 0; +} diff --git a/Polyhedron/demo/Polyhedron/Plugins/PCA/Edit_box_plugin.cpp b/Polyhedron/demo/Polyhedron/Plugins/PCA/Edit_box_plugin.cpp index 6252a0a2bfe..6fe2b67ef33 100644 --- a/Polyhedron/demo/Polyhedron/Plugins/PCA/Edit_box_plugin.cpp +++ b/Polyhedron/demo/Polyhedron/Plugins/PCA/Edit_box_plugin.cpp @@ -111,6 +111,7 @@ void Edit_box_plugin::exportToPoly() break; } } + Polyhedron::Point_3 points[8]; for(int i=0; i<8; ++i) { diff --git a/Polyhedron/demo/Polyhedron/Plugins/PMP/CMakeLists.txt b/Polyhedron/demo/Polyhedron/Plugins/PMP/CMakeLists.txt index 42af8d59154..fbb1fd7f12d 100644 --- a/Polyhedron/demo/Polyhedron/Plugins/PMP/CMakeLists.txt +++ b/Polyhedron/demo/Polyhedron/Plugins/PMP/CMakeLists.txt @@ -75,3 +75,6 @@ target_link_libraries(repair_polyhedron_plugin scene_polyhedron_item) qt5_wrap_ui( isotropicRemeshingUI_FILES Isotropic_remeshing_dialog.ui) polyhedron_demo_plugin(isotropic_remeshing_plugin Isotropic_remeshing_plugin ${isotropicRemeshingUI_FILES}) target_link_libraries(isotropic_remeshing_plugin scene_polyhedron_item scene_polyhedron_selection_item) + +polyhedron_demo_plugin(distance_plugin Distance_plugin) +target_link_libraries(distance_plugin scene_polyhedron_item scene_color_ramp) diff --git a/Polyhedron/demo/Polyhedron/Plugins/PMP/Distance_plugin.cpp b/Polyhedron/demo/Polyhedron/Plugins/PMP/Distance_plugin.cpp new file mode 100644 index 00000000000..41b29c8fada --- /dev/null +++ b/Polyhedron/demo/Polyhedron/Plugins/PMP/Distance_plugin.cpp @@ -0,0 +1,467 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include "Messages_interface.h" +#include "Scene_polyhedron_item.h" +#include "Color_ramp.h" +#include "triangulate_primitive.h" +#include +#include +#include +#include +#include +#include +#include + +using namespace CGAL::Three; +namespace PMP = CGAL::Polygon_mesh_processing; + +#ifdef CGAL_LINKED_WITH_TBB +template +struct Distance_computation{ + const AABB_tree& tree; + const std::vector& sample_points; + Point_3 initial_hint; + CGAL::cpp11::atomic* distance; + std::vector& output; + + Distance_computation(const AABB_tree& tree, + const Point_3 p, + const std::vector& sample_points, + CGAL::cpp11::atomic* d, + std::vector& out ) + : tree(tree) + , sample_points(sample_points) + , initial_hint(p) + , distance(d) + , output(out) + { + } + void + operator()(const tbb::blocked_range& range) const + { + Point_3 hint = initial_hint; + double hdist = 0; + for( std::size_t i = range.begin(); i != range.end(); ++i) + { + hint = tree.closest_point(sample_points[i], hint); + Kernel::FT dist = squared_distance(hint,sample_points[i]); + double d = CGAL::sqrt(dist); + output[i] = d; + if (d>hdist) hdist=d; + } + + if (hdist > distance->load(CGAL::cpp11::memory_order_acquire)) + distance->store(hdist, CGAL::cpp11::memory_order_release); + } +}; +#endif + +class Scene_distance_polyhedron_item: public Scene_item +{ + Q_OBJECT +public: + Scene_distance_polyhedron_item(Polyhedron* poly, Polyhedron* polyB, QString other_name, int sampling_pts) + :Scene_item(NbOfVbos,NbOfVaos), + poly(poly), + poly_B(polyB), + are_buffers_filled(false), + other_poly(other_name) + { + nb_pts_per_face = sampling_pts; + this->setRenderingMode(FlatPlusEdges); + thermal_ramp.build_thermal(); + } + bool supportsRenderingMode(RenderingMode m) const { + return (m == Flat || m == FlatPlusEdges); + } + Scene_item* clone() const {return 0;} + QString toolTip() const {return QString("Item %1 with color indicating distance with %2").arg(this->name()).arg(other_poly);} + void draw(Viewer_interface *viewer) const + { + if(!are_buffers_filled) + { + computeElements(); + initializeBuffers(viewer); + compute_bbox(); + } + vaos[Facets]->bind(); + attribBuffers(viewer, PROGRAM_WITH_LIGHT); + program = getShaderProgram(PROGRAM_WITH_LIGHT); + program->bind(); + program->setUniformValue("is_selected", false); + viewer->glDrawArrays(GL_TRIANGLES, 0, static_cast(nb_pos/3)); + program->release(); + vaos[Facets]->release(); + } + void drawEdges(Viewer_interface* viewer) const + { + vaos[Edges]->bind(); + + attribBuffers(viewer, PROGRAM_WITHOUT_LIGHT); + program = getShaderProgram(PROGRAM_WITHOUT_LIGHT); + program->bind(); + //draw the edges + program->setAttributeValue("colors", QColor(Qt::black)); + program->setUniformValue("is_selected", false); + viewer->glDrawArrays(GL_LINES, 0, static_cast(nb_edge_pos/3)); + vaos[Edges]->release(); + program->release(); + } + void compute_bbox() const { + const Kernel::Point_3& p = *(poly->points_begin()); + CGAL::Bbox_3 bbox(p.x(), p.y(), p.z(), p.x(), p.y(), p.z()); + for(Polyhedron::Point_iterator it = poly->points_begin(); + it != poly->points_end(); + ++it) { + bbox = bbox + it->bbox(); + } + _bbox = Bbox(bbox.xmin(),bbox.ymin(),bbox.zmin(), + bbox.xmax(),bbox.ymax(),bbox.zmax()); + } +private: + Polyhedron* poly; + Polyhedron* poly_B; + mutable bool are_buffers_filled; + QString other_poly; + mutable std::vector vertices; + mutable std::vector edge_vertices; + mutable std::vector normals; + mutable std::vector colors; + Color_ramp thermal_ramp; + int nb_pts_per_face; + + enum VAOs { + Facets=0, + Edges, + NbOfVaos}; + + enum VBOs { + Vertices=0, + Edge_vertices, + Normals, + Colors, + NbOfVbos}; + + mutable std::size_t nb_pos; + mutable std::size_t nb_edge_pos; + mutable QOpenGLShaderProgram *program; + + //fills 'out' and returns the hausdorff distance for calibration of the color_ramp. + + double compute_distances(const Polyhedron& m, const std::vector& sample_points, + std::vector& out)const + { + typedef CGAL::AABB_face_graph_triangle_primitive Primitive; + typedef CGAL::AABB_traits Traits; + typedef CGAL::AABB_tree< Traits > Tree; + + Tree tree( faces(m).first, faces(m).second, m); + tree.accelerate_distance_queries(); + tree.build(); + + typename Traits::Point_3 hint = m.vertices_begin()->point(); + +#ifndef CGAL_LINKED_WITH_TBB + double hdist = 0; + for(std::size_t i = 0; ihdist) hdist=d; + } + return hdist; +#else + CGAL::cpp11::atomic distance(0); + Distance_computation f(tree, hint, sample_points, &distance, out); + tbb::parallel_for(tbb::blocked_range(0, sample_points.size()), f); + return distance; +#endif + } + + void computeElements()const + { + QApplication::setOverrideCursor(Qt::WaitCursor); + vertices.resize(0); + edge_vertices.resize(0); + normals.resize(0); + colors.resize(0); + + typedef Polyhedron::Traits Kernel; + typedef Kernel::Vector_3 Vector; + typedef Polyhedron::Facet_iterator Facet_iterator; + typedef boost::graph_traits::face_descriptor face_descriptor; + typedef boost::graph_traits::vertex_descriptor vertex_descriptor; + + //facets + { + boost::container::flat_map face_normals_map; + boost::associative_property_map< boost::container::flat_map > + nf_pmap(face_normals_map); + boost::container::flat_map vertex_normals_map; + boost::associative_property_map< boost::container::flat_map > + nv_pmap(vertex_normals_map); + + PMP::compute_normals(*poly, nv_pmap, nf_pmap); + std::vector total_points(0); + Facet_iterator f = poly->facets_begin(); + for(f = poly->facets_begin(); + f != poly->facets_end(); + f++) + { + Vector nf = get(nf_pmap, f); + f->plane() = Kernel::Plane_3(f->halfedge()->vertex()->point(), nf); + typedef FacetTriangulator::vertex_descriptor> FT; + double diagonal; + if(this->diagonalBbox() != std::numeric_limits::infinity()) + diagonal = this->diagonalBbox(); + else + diagonal = 0.0; + + //compute distance with other polyhedron + //sample facet + std::vector sampled_points; + std::size_t nb_points = (std::max)((int)std::ceil(nb_pts_per_face * PMP::face_area(f,*poly,PMP::parameters::geom_traits(Kernel()))), + 1); + CGAL::Random_points_in_triangle_3 g(f->halfedge()->vertex()->point(), f->halfedge()->next()->vertex()->point(), + f->halfedge()->next()->next()->vertex()->point()); + CGAL::cpp11::copy_n(g, nb_points, std::back_inserter(sampled_points)); + sampled_points.push_back(f->halfedge()->vertex()->point()); + sampled_points.push_back(f->halfedge()->next()->vertex()->point()); + sampled_points.push_back(f->halfedge()->next()->next()->vertex()->point()); + + //triangle facets with sample points for color display + FT triangulation(f,sampled_points,nf,poly,diagonal); + + if(triangulation.cdt->dimension() != 2 ) + { + qDebug()<<"Error : cdt not right (dimension != 2). Facet not displayed"; + return; + } + + //iterates on the internal faces to add the vertices to the positions + //and the normals to the appropriate vectors + + for(FT::CDT::Finite_faces_iterator + ffit = triangulation.cdt->finite_faces_begin(), + end = triangulation.cdt->finite_faces_end(); + ffit != end; ++ffit) + { + if(ffit->info().is_external) + continue; + + for (int i = 0; i<3; ++i) + { + total_points.push_back(ffit->vertex(i)->point()); + vertices.push_back(ffit->vertex(i)->point().x()); + vertices.push_back(ffit->vertex(i)->point().y()); + vertices.push_back(ffit->vertex(i)->point().z()); + + normals.push_back(nf.x()); + normals.push_back(nf.y()); + normals.push_back(nf.z()); + } + } + } + //compute the distances + typedef CGAL::Spatial_sort_traits_adapter_3::type > Search_traits_3; + + std::vector distances(total_points.size()); + std::vector indices; + indices.reserve(total_points.size()); + std::copy(boost::counting_iterator(0), + boost::counting_iterator(total_points.size()), + std::back_inserter(indices)); + spatial_sort(indices.begin(), + indices.end(), + Search_traits_3(CGAL::make_property_map(total_points))); + std::vector sorted_points(total_points.size()); + for(std::size_t i = 0; i < sorted_points.size(); ++i) + { + sorted_points[i] = total_points[indices[i]]; + } + + double hausdorff = compute_distances(*poly_B, + sorted_points, + distances); + //compute the colors + colors.resize(sorted_points.size()*3); + for(std::size_t i=0; iedges_begin(); + he != poly->edges_end(); + he++) + { + const Point& a = he->vertex()->point(); + const Point& b = he->opposite()->vertex()->point(); + { + + edge_vertices.push_back(a.x()); + edge_vertices.push_back(a.y()); + edge_vertices.push_back(a.z()); + + edge_vertices.push_back(b.x()); + edge_vertices.push_back(b.y()); + edge_vertices.push_back(b.z()); + } + } + } + QApplication::restoreOverrideCursor(); + } + void initializeBuffers(Viewer_interface *viewer)const + { + + program = getShaderProgram(PROGRAM_WITH_LIGHT, viewer); + program->bind(); + vaos[Facets]->bind(); + buffers[Vertices].bind(); + buffers[Vertices].allocate(vertices.data(), + static_cast(vertices.size()*sizeof(float))); + program->enableAttributeArray("vertex"); + program->setAttributeBuffer("vertex",GL_FLOAT,0,3); + buffers[Vertices].release(); + buffers[Normals].bind(); + buffers[Normals].allocate(normals.data(), + static_cast(normals.size()*sizeof(float))); + program->enableAttributeArray("normals"); + program->setAttributeBuffer("normals",GL_FLOAT,0,3); + buffers[Normals].release(); + buffers[Colors].bind(); + buffers[Colors].allocate(colors.data(), + static_cast(colors.size()*sizeof(float))); + program->enableAttributeArray("colors"); + program->setAttributeBuffer("colors",GL_FLOAT,0,3); + buffers[Colors].release(); + vaos[Facets]->release(); + program->release(); + + program = getShaderProgram(PROGRAM_WITHOUT_LIGHT, viewer); + program->bind(); + vaos[Edges]->bind(); + buffers[Edge_vertices].bind(); + buffers[Edge_vertices].allocate(edge_vertices.data(), + static_cast(edge_vertices.size()*sizeof(float))); + program->enableAttributeArray("vertex"); + program->setAttributeBuffer("vertex",GL_FLOAT,0,3); + buffers[Edge_vertices].release(); + vaos[Facets]->release(); + program->release(); + + nb_pos = vertices.size(); + vertices.resize(0); + //"Swap trick" insures that the memory is indeed freed and not kept available + std::vector(vertices).swap(vertices); + nb_edge_pos = edge_vertices.size(); + edge_vertices.resize(0); + std::vector(edge_vertices).swap(edge_vertices); + normals.resize(0); + std::vector(normals).swap(normals); + colors.resize(0); + std::vector(colors).swap(colors); + are_buffers_filled = true; + } +}; +class DistancePlugin : + public QObject, + public Polyhedron_demo_plugin_interface +{ + Q_OBJECT + Q_INTERFACES(CGAL::Three::Polyhedron_demo_plugin_interface) + Q_PLUGIN_METADATA(IID "com.geometryfactory.PolyhedronDemo.PluginInterface/1.0") + + typedef Kernel::Point_3 Point_3; +public: + //decides if the plugin's actions will be displayed or not. + bool applicable(QAction*) const + { + + return scene->selectionIndices().size() == 2 && + qobject_cast(scene->item(scene->selectionIndices().first())) && + qobject_cast(scene->item(scene->selectionIndices().last())); + + } + //the list of the actions of the plugin. + QList actions() const + { + return _actions; + } + //this acts like a constructor for the plugin. It gets the references to the mainwindow and the scene, and connects the action. + void init(QMainWindow* mw, Scene_interface* sc, Messages_interface* mi) + { + //gets the reference to the message interface, to display text in the console widget + this->messageInterface = mi; + //get the references + this->scene = sc; + this->mw = mw; + //creates the action + QAction *actionComputeDistance= new QAction(QString("Compute Distance Between Polyhedra"), mw); + //specifies the subMenu + actionComputeDistance->setProperty("submenuName", "Polygon Mesh Processing"); + //links the action + if(actionComputeDistance) { + connect(actionComputeDistance, SIGNAL(triggered()), + this, SLOT(createDistanceItems())); + _actions << actionComputeDistance; + } + } +public Q_SLOTS: + void createDistanceItems() + { + bool ok = false; + nb_pts_per_face = QInputDialog::getInt(mw, tr("Sampling"), + tr("Number of points per face:"),400, 1,2147483647,1, &ok); + if (!ok) + return; + + //check the initial conditions + Scene_polyhedron_item* itemA = qobject_cast(scene->item(scene->selectionIndices().first())); + Scene_polyhedron_item* itemB = qobject_cast(scene->item(scene->selectionIndices().last())); + if(!itemA->polyhedron()->is_pure_triangle() || + !itemB->polyhedron()->is_pure_triangle() ){ + messageInterface->error(QString("Distance not computed. (Both polyhedra must be triangulated)")); + return; + } + QApplication::setOverrideCursor(Qt::WaitCursor); + Scene_distance_polyhedron_item* new_itemA = new Scene_distance_polyhedron_item(itemA->polyhedron(),itemB->polyhedron(), itemB->name(), nb_pts_per_face); + Scene_distance_polyhedron_item* new_itemB = new Scene_distance_polyhedron_item(itemB->polyhedron(),itemA->polyhedron(), itemA->name(), nb_pts_per_face); + itemA->setVisible(false); + itemB->setVisible(false); + new_itemA->setName(QString("%1 to %2").arg(itemA->name()).arg(itemB->name())); + new_itemB->setName(QString("%1 to %2").arg(itemB->name()).arg(itemA->name())); + scene->addItem(new_itemA); + scene->addItem(new_itemB); + QApplication::restoreOverrideCursor(); + } +private: + int nb_pts_per_face; + QList _actions; + Messages_interface* messageInterface; + //The reference to the scene + Scene_interface* scene; + //The reference to the main window + QMainWindow* mw; +}; +#include "Distance_plugin.moc" diff --git a/Polyhedron/demo/Polyhedron/Plugins/Point_set/Surface_reconstruction_plugin_impl.cpp b/Polyhedron/demo/Polyhedron/Plugins/Point_set/Surface_reconstruction_plugin_impl.cpp index e167db904cb..10c4fc22c8d 100644 --- a/Polyhedron/demo/Polyhedron/Plugins/Point_set/Surface_reconstruction_plugin_impl.cpp +++ b/Polyhedron/demo/Polyhedron/Plugins/Point_set/Surface_reconstruction_plugin_impl.cpp @@ -237,8 +237,6 @@ Polyhedron* poisson_reconstruct(Point_set& points, } } - - return output_mesh; } diff --git a/Polyhedron/demo/Polyhedron/triangulate_primitive.h b/Polyhedron/demo/Polyhedron/triangulate_primitive.h index 3b45a18b2b4..be9c0b6334c 100644 --- a/Polyhedron/demo/Polyhedron/triangulate_primitive.h +++ b/Polyhedron/demo/Polyhedron/triangulate_primitive.h @@ -7,7 +7,7 @@ #include #include #include - +#include #include @@ -67,6 +67,23 @@ public: } triangulate(idPoints, normal, item_diag); } + FacetTriangulator(typename boost::graph_traits::face_descriptor fd, + const std::vector& more_points, + const Vector& normal, + Mesh *poly, + const double item_diag) + { + std::vector idPoints; + BOOST_FOREACH(halfedge_descriptor he_circ, halfedges_around_face( halfedge(fd, *poly), *poly)) + { + PointAndId idPoint; + idPoint.point = get(boost::vertex_point,*poly,source(he_circ, *poly)); + idPoint.id = source(he_circ, *poly); + idPoints.push_back(idPoint); + + } + triangulate_with_points(idPoints,more_points, normal, item_diag); + } FacetTriangulator(std::vector &idPoints, const Vector& normal, @@ -74,6 +91,13 @@ public: { triangulate(idPoints, normal, item_diag); } + FacetTriangulator(std::vector &idPoints, + const std::vector& more_points, + const Vector& normal, + const double item_diag) + { + triangulate_with_points(idPoints, more_points, normal, item_diag); + } private: void triangulate( std::vector &idPoints, @@ -89,11 +113,16 @@ private: double min_sq_dist = CGAL::square(0.0001*item_diag); // Iterate the points of the facet and decide if they must be inserted in the CDT + typename Kernel::FT x(0), y(0), z(0); + BOOST_FOREACH(PointAndId idPoint, idPoints) { + x += idPoint.point.x(); + y += idPoint.point.y(); + z += idPoint.point.z(); typename CDT::Vertex_handle vh; //Always insert the first point, then only insert - // if the distance with the precedent is reasonable. + // if the distance with the previous is reasonable. if(first == 0 || CGAL::squared_distance(idPoint.point, previous->point()) > min_sq_dist) { vh = cdt->insert(idPoint.point); @@ -150,6 +179,86 @@ private: } } } + + void triangulate_with_points( std::vector &idPoints, + const std::vector& more_points, + const Vector& normal, + const double item_diag ) + { + P_traits cdt_traits(normal); + cdt = new CDT(cdt_traits); + typename CDT::Vertex_handle previous, first, last_inserted; + //Compute a reasonable precision level used to decide + //if two consecutive points in a facet can be estimated + //equal. + + double min_sq_dist = CGAL::square(0.0001*item_diag); + // Iterate the points of the facet and decide if they must be inserted in the CDT + BOOST_FOREACH(PointAndId idPoint, idPoints) + { + typename CDT::Vertex_handle vh; + //Always insert the first point, then only insert + // if the distance with the previous is reasonable. + if(first == 0 || CGAL::squared_distance(idPoint.point, previous->point()) > min_sq_dist) + { + vh = cdt->insert(idPoint.point); + v2v[vh] = idPoint.id; + if(first == 0) { + first = vh; + } + if(previous != 0 && previous != vh) { + double sq_dist = CGAL::squared_distance(previous->point(), vh->point()); + if(sq_dist > min_sq_dist) + { + cdt->insert_constraint(previous, vh); + sq_dist = CGAL::squared_distance(previous->point(), first->point()); + if(sq_dist > min_sq_dist) + { + last_inserted = previous; + } + } + } + previous = vh; + } + } + double sq_dist = CGAL::squared_distance(previous->point(), first->point()); + + if(sq_dist > min_sq_dist) + { + cdt->insert_constraint(previous, first); + } + else + { + cdt->insert_constraint(last_inserted, first); + } + BOOST_FOREACH(typename Kernel::Point_3 point, more_points) + { + cdt->insert(point); + } + // sets mark is_external + for(typename CDT::All_faces_iterator + fit2 = cdt->all_faces_begin(), + end = cdt->all_faces_end(); + fit2 != end; ++fit2) + { + fit2->info().is_external = false; + } + //check if the facet is external or internal + std::queue face_queue; + face_queue.push(cdt->infinite_vertex()->face()); + while(! face_queue.empty() ) { + typename CDT::Face_handle fh = face_queue.front(); + face_queue.pop(); + if(fh->info().is_external) continue; + fh->info().is_external = true; + for(int i = 0; i <3; ++i) { + if(!cdt->is_constrained(std::make_pair(fh, i))) + { + face_queue.push(fh->neighbor(i)); + } + } + } + } }; #endif // TRIANGULATE_PRIMITIVE diff --git a/Spatial_searching/include/CGAL/Orthogonal_k_neighbor_search.h b/Spatial_searching/include/CGAL/Orthogonal_k_neighbor_search.h index 308d2937247..4dc03058bf6 100644 --- a/Spatial_searching/include/CGAL/Orthogonal_k_neighbor_search.h +++ b/Spatial_searching/include/CGAL/Orthogonal_k_neighbor_search.h @@ -145,7 +145,7 @@ private: FT diff2 = val - node->upper_low_value(); if ( (diff1 + diff2 >= FT(0.0)) ) { - new_off = 2*val < node->upper_low_value()+node->upper_high_value() ? + new_off = node->upper_low_value()+node->upper_high_value() > val*2? val - node->upper_high_value(): val - node->upper_low_value(); bestChild = node->lower(); @@ -153,7 +153,7 @@ private: } else // compute new distance { - new_off = 2*val < node->lower_low_value()+node->lower_high_value() ? + new_off = node->lower_low_value()+node->lower_high_value() > val*2 ? val - node->lower_high_value(): val - node->lower_low_value(); bestChild = node->upper();