diff --git a/Generator/include/CGAL/point_generators_3.h b/Generator/include/CGAL/point_generators_3.h index 3f5d27c475f..79bb2f54525 100644 --- a/Generator/include/CGAL/point_generators_3.h +++ b/Generator/include/CGAL/point_generators_3.h @@ -552,20 +552,20 @@ public: template class Triangle_from_soup { - typedef typename PointRange::value_type Point_3; - typedef typename Kernel_traits::Kernel Kernel; + typedef typename boost::range_value::type Point_3; + typedef typename Kernel_traits::Kernel Kernel; +public: + typedef typename Kernel::Triangle_3 result_type; +private: const PointRange& points; public: Triangle_from_soup(const PointRange& pts) :points(pts) {} - typedef typename Kernel::Triangle_3 result_type; result_type operator()(const Triangle& t) const { - return result_type(points[t[0]], - points[t[1]], - points[t[2]]); + return result_type(points[t[0]], points[t[1]], points[t[2]]); } }; @@ -720,7 +720,7 @@ struct Random_points_in_triangles_3 template , + class Triangle = std::vector, class Creator = Creator_uniform_3< typename Kernel_traits< typename PointRange::value_type >::Kernel::RT, typename PointRange::value_type> @@ -734,26 +734,29 @@ struct Random_points_in_triangle_soup typedef Generic_random_point_generator, Random_points_in_triangle_3, - typename PointRange::value_type> Base; - - typedef typename PointRange::value_type Point_3; - typedef typename Kernel_traits::Kernel Kernel; - typedef Triangle Id; - typedef Point_3 result_type; - typedef Random_points_in_triangle_soup This; + typename PointRange::value_type> Base; + typedef typename PointRange::value_type Point_3; + typedef typename Kernel_traits::Kernel Kernel; + typedef Triangle Id; + typedef Point_3 result_type; + typedef Random_points_in_triangle_soup This; template - Random_points_in_triangle_soup( const TriangleRange& triangles, const PointRange& points, Random& rnd = get_default_random()) + Random_points_in_triangle_soup(const TriangleRange& triangles, + const PointRange& points, + Random& rnd = get_default_random()) : Base(triangles, internal::Triangle_from_soup(points), internal::Apply_approx_sqrt::Kernel::Compute_squared_area_3>() ,rnd ) { } + This& operator++() { Base::generate_point(); return *this; } + This operator++(int) { This tmp = *this; ++(*this); 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 e99e6b14782..6c6149037e2 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/distance.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/distance.h @@ -160,56 +160,85 @@ double approximate_Hausdorff_distance_impl( } template -struct Triangle_structure_sampler_base{ +struct Triangle_structure_sampler_base +{ NamedParameters np; - Geom_traits geomtraits; + GeomTraits geomtraits; OutputIterator& out; + Triangle_structure_sampler_base(OutputIterator& out, NamedParameters np) :np(np), out(out) {} + void sample_points(); + double get_minimum_edge_length(); + template double get_tr_area(const Tr&); + template - void get_tr_points(const Tr& tr, typename Geom_traits::Point_3 points[]); - void ms_edges_sample(std::size_t nb_points_per_edge, - std::size_t nb_pts_l_u); + void get_tr_points(const Tr& tr, typename GeomTraits::Point_3 points[]); + + void ms_edges_sample(const std::size_t& nb_points_per_edge, + const std::size_t& nb_pts_l_u); + void ru_edges_sample(); + Randomizer get_randomizer(); + void internal_sample_triangles(double, bool, bool); + std::pair get_range(); + std::size_t get_points_size(); + void procede() { using parameters::choose_parameter; using parameters::get_parameter; using parameters::is_default_parameter; - geomtraits = choose_parameter(get_parameter(np, internal_np::geom_traits), Geom_traits()); + geomtraits = choose_parameter(get_parameter(np, internal_np::geom_traits), GeomTraits()); - bool use_rs = choose_parameter(get_parameter(np, internal_np::random_uniform_sampling), true); - bool use_gs = choose_parameter(get_parameter(np, internal_np::grid_sampling), false); + bool use_rs = choose_parameter + (get_parameter(np, internal_np::random_uniform_sampling), true); + + bool use_gs = choose_parameter( + get_parameter(np, internal_np::grid_sampling), false); + bool use_ms = choose_parameter(get_parameter(np, internal_np::monte_carlo_sampling), false); if (use_gs || use_ms) + { if (is_default_parameter(get_parameter(np, internal_np::random_uniform_sampling))) + { use_rs=false; + } + } - bool smpl_vrtcs = choose_parameter(get_parameter(np, internal_np::do_sample_vertices), true); - bool smpl_dgs = choose_parameter(get_parameter(np, internal_np::do_sample_edges), true); - bool smpl_fcs = choose_parameter(get_parameter(np, internal_np::do_sample_faces), true); + bool smpl_vrtcs + = choose_parameter(get_parameter(np, internal_np::do_sample_vertices), true); - double nb_pts_a_u = choose_parameter(get_parameter(np, internal_np::nb_points_per_area_unit), 0.); - double nb_pts_l_u = choose_parameter(get_parameter(np, internal_np::nb_points_per_distance_unit), 0.); + bool smpl_dgs + = choose_parameter(get_parameter(np, internal_np::do_sample_edges), true); + + bool smpl_fcs + = choose_parameter(get_parameter(np, internal_np::do_sample_faces), true); + + double nb_pts_a_u + = choose_parameter(get_parameter(np, internal_np::nb_points_per_area_unit), 0.); + + double nb_pts_l_u + = choose_parameter(get_parameter(np, internal_np::nb_points_per_distance_unit), 0.); // sample vertices if (smpl_vrtcs) @@ -220,12 +249,15 @@ struct Triangle_structure_sampler_base{ // grid sampling if (use_gs) { - double grid_spacing_ = choose_parameter(get_parameter(np, internal_np::grid_spacing), 0.); + double grid_spacing_ + = choose_parameter(get_parameter(np, internal_np::grid_spacing), 0.); + if (grid_spacing_==0.) { // set grid spacing to the shortest edge length grid_spacing_ = static_cast(this)->get_minimum_edge_length(); } + static_cast(this)->internal_sample_triangles( grid_spacing_, smpl_fcs, smpl_dgs); } @@ -237,6 +269,7 @@ struct Triangle_structure_sampler_base{ std::size_t nb_points_per_face = choose_parameter(get_parameter(np, internal_np::number_of_points_per_face), 0); + std::size_t nb_points_per_edge = choose_parameter(get_parameter(np, internal_np::number_of_points_per_edge), 0); @@ -263,17 +296,18 @@ struct Triangle_structure_sampler_base{ nb_points = (std::max)( static_cast( std::ceil(static_cast(this)->get_tr_area(tr)) - *nb_pts_a_u) + *nb_pts_a_u) ,std::size_t(1)); } // extract triangle face points - typename Geom_traits::Point_3 points[3]; + typename GeomTraits::Point_3 points[3]; static_cast(this)->get_tr_points(tr, points); - Random_points_in_triangle_3 + 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) { @@ -287,18 +321,25 @@ struct Triangle_structure_sampler_base{ // sample faces if(smpl_fcs) { - std::size_t nb_points = choose_parameter(get_parameter(np, internal_np::number_of_points_on_faces), 0); + std::size_t nb_points + = choose_parameter(get_parameter(np, internal_np::number_of_points_on_faces), 0); + typename Derived::Randomizer g = static_cast(this)->get_randomizer(); if (nb_points == 0) { if (nb_pts_a_u == 0.) + { nb_points = static_cast(this)->get_points_size(); + } else + { nb_points = static_cast( std::ceil(g.sum_of_weights()*nb_pts_a_u) ); + } } out = CGAL::cpp11::copy_n(g, nb_points, out); } + // sample edges if (smpl_dgs) { @@ -383,33 +424,33 @@ namespace internal{ template struct Triangle_structure_sampler_for_triangle_mesh : Triangle_structure_sampler_base< OutputIterator, - Geom_traits, + GeomTraits, NamedParameters, typename boost::graph_traits::face_iterator, Random_points_in_triangle_mesh_3, Creator, Triangle_structure_sampler_for_triangle_mesh > { typedef Triangle_structure_sampler_for_triangle_mesh This; typedef Triangle_structure_sampler_base< OutputIterator, - Geom_traits, + GeomTraits, NamedParameters, typename boost::graph_traits::face_iterator, Random_points_in_triangle_mesh_3, @@ -463,7 +504,7 @@ struct Triangle_structure_sampler_for_triangle_mesh BOOST_FOREACH(edge_descriptor ed, edges(tm)) { double el = std::sqrt( - to_double( typename Geom_traits::Compute_squared_distance_3()( + to_double( typename GeomTraits::Compute_squared_distance_3()( get(pmap, source(ed, tm)), get(pmap, target(ed, tm)) ))); if (el > 0 && el < min_edge_length_) min_edge_length_ = el; @@ -477,7 +518,7 @@ struct Triangle_structure_sampler_for_triangle_mesh } template//tr = face_descriptor here - void get_tr_points(const Tr& tr, typename Geom_traits::Point_3 points[]) + void get_tr_points(const Tr& tr, typename GeomTraits::Point_3 points[]) { halfedge_descriptor hd(halfedge(tr,tm)); for(int i=0; i<3; ++i) @@ -489,7 +530,7 @@ struct Triangle_structure_sampler_for_triangle_mesh void ms_edges_sample(std::size_t nb_points_per_edge, double nb_pts_l_u) { - typename Geom_traits::Compute_squared_distance_3 squared_distance = this->geomtraits.compute_squared_distance_3_object(); + typename GeomTraits::Compute_squared_distance_3 squared_distance = this->geomtraits.compute_squared_distance_3_object(); 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)) @@ -504,7 +545,7 @@ struct Triangle_structure_sampler_for_triangle_mesh std::size_t(1)); } // now do the sampling of the edge - Random_points_on_segment_3 + Random_points_on_segment_3 g(get(pmap, source(ed,tm)), get(pmap, target(ed,tm))); this->out=CGAL::cpp11::copy_n(g, nb_points, this->out); } @@ -533,7 +574,7 @@ struct Triangle_structure_sampler_for_triangle_mesh } void internal_sample_triangles(double grid_spacing_, bool smpl_fcs, bool smpl_dgs) { - this->out = sample_triangles( + this->out = sample_triangles( faces(tm), tm, pmap, grid_spacing_, this->out,smpl_fcs, smpl_dgs, false); } @@ -548,56 +589,55 @@ struct Triangle_structure_sampler_for_triangle_mesh template struct Triangle_structure_sampler_for_triangle_soup : Triangle_structure_sampler_base< OutputIterator, - Geom_traits, + GeomTraits, NamedParameters, typename TriangleRange::const_iterator, Random_points_in_triangle_soup, Creator, Triangle_structure_sampler_for_triangle_soup + TriangleRange, + OutputIterator, + GeomTraits, + Creator, + NamedParameters> > { typedef typename TriangleRange::value_type TriangleType; typedef Triangle_structure_sampler_for_triangle_soup This; typedef Triangle_structure_sampler_base< OutputIterator, - Geom_traits, + GeomTraits, NamedParameters, typename TriangleRange::const_iterator, Random_points_in_triangle_soup, Creator, - This> Base; - + This> Base; typedef Random_points_in_triangle_soup Randomizer; typedef typename TriangleRange::const_iterator TriangleIterator; double min_edge_length_; - const TriangleRange& triangles; const PointRange& points; + const TriangleRange& triangles; + Triangle_structure_sampler_for_triangle_soup(const PointRange& pts, const TriangleRange& trs, OutputIterator& out, NamedParameters np) - :Base(out, np), - points(pts), triangles(trs) + :Base(out, np), points(pts), triangles(trs) { using parameters::choose_parameter; using parameters::get_parameter; @@ -609,6 +649,7 @@ struct Triangle_structure_sampler_for_triangle_soup { return std::make_pair(triangles.begin(), triangles.end()); } + void sample_points() { this->out = std::copy( @@ -616,19 +657,20 @@ struct Triangle_structure_sampler_for_triangle_soup points.end(), this->out); } + double get_minimum_edge_length() { if(min_edge_length_ != (std::numeric_limits::max)()) return min_edge_length_; - for(auto tr : triangles) + for(const auto& tr : triangles) { for(std::size_t i = 0; i< 3; ++i) { - typename Geom_traits::Point_3 a(points[tr[i]]), + typename GeomTraits::Point_3 a(points[tr[i]]), b(points[tr[(i+1)%3]]); double el = std::sqrt( - to_double( typename Geom_traits::Compute_squared_distance_3()( + to_double( typename GeomTraits::Compute_squared_distance_3()( a, b))); if (el > 0 && el < min_edge_length_) min_edge_length_ = el; @@ -636,6 +678,7 @@ struct Triangle_structure_sampler_for_triangle_soup } return min_edge_length_; } + template double get_tr_area(const Tr& tr) { @@ -645,35 +688,39 @@ struct Triangle_structure_sampler_for_triangle_soup } template - void get_tr_points(const Tr& tr, typename Geom_traits::Point_3 points[]) + void get_tr_points(const Tr& tr, typename GeomTraits::Point_3 points[]) { for(int i=0; i<3; ++i) { points[i] = this->points[tr[i] ]; } } + void ms_edges_sample(std::size_t , double ) { //don't sample edges in soup. } + void ru_edges_sample(double, double) { //don't sample edges in soup. } + Randomizer get_randomizer() { return Randomizer(triangles, points); } + void internal_sample_triangles(double distance, bool, bool) { - for(auto tr : triangles) + for(const auto& tr : triangles) { - typename Geom_traits::Point_3 p0 = points[tr[0] ]; - typename Geom_traits::Point_3 p1 = points[tr[1] ]; - typename Geom_traits::Point_3 p2 = points[tr[2] ]; - this->out=internal::triangle_grid_sampling(p0, p1, p2, distance, + typename GeomTraits::Point_3 p0 = points[tr[0] ]; + typename GeomTraits::Point_3 p1 = points[tr[1] ]; + typename GeomTraits::Point_3 p2 = points[tr[2] ]; + this->out=internal::triangle_grid_sampling(p0, p1, p2, distance, this->out); } } @@ -682,6 +729,7 @@ struct Triangle_structure_sampler_for_triangle_soup { return points.size(); } + }; }//end internal @@ -790,14 +838,14 @@ sample_triangle_mesh(const TriangleMesh& tm, NamedParameters np) { typedef typename GetGeomTraits::type Geom_traits; + NamedParameters>::type GeomTraits; typedef typename GetVertexPointMap::const_type Vpm; internal::Triangle_structure_sampler_for_triangle_mesh, + GeomTraits, + Creator_uniform_3, Vpm, NamedParameters> performer(tm, out, np); @@ -810,7 +858,7 @@ sample_triangle_mesh(const TriangleMesh& tm, * generates points taken on `triangles` and outputs them to `out`, the sampling method * is selected using named parameters. * @tparam PointRange a model of the concept `RandomAccessContainer` whose value type is the point type. - * @tparam TriangleMesh a model of the concept `RandomAccessContainer` + * @tparam TriangleRange a model of the concept `RandomAccessContainer` * whose value_type is itself a model of the concept `RandomAccessContainer` * whose value_type is `std::size_t`. * @tparam OutputIterator a model of `OutputIterator` @@ -885,20 +933,21 @@ sample_triangle_soup(const PointRange& points, NamedParameters np) { typedef typename PointRange::value_type Point_3; - typedef typename Kernel_traits::Kernel Geom_traits; + typedef typename Kernel_traits::Kernel GeomTraits; internal::Triangle_structure_sampler_for_triangle_soup, + GeomTraits, + Creator_uniform_3, NamedParameters> performer(points, triangles, out, np); performer.procede(); return performer.out; } + template OutputIterator sample_triangle_mesh(const TriangleMesh& tm, @@ -907,6 +956,18 @@ sample_triangle_mesh(const TriangleMesh& tm, return sample_triangle_mesh(tm, out, parameters::all_default()); } +template +OutputIterator +sample_triangle_soup(const PointRange& points, + const TriangleRange& triangles, + OutputIterator out) +{ + return sample_triangle_soup(points, triangles, out, parameters::all_default()); +} + + template ::type Geom_traits; + NamedParameters1>::type GeomTraits; - return approximate_Hausdorff_distance( + return approximate_Hausdorff_distance( tm1, tm2, np1, parameters::choose_parameter(parameters::get_parameter(np2, internal_np::vertex_point), get_const_property_map(vertex_point, tm2))); } @@ -1062,9 +1123,9 @@ double max_distance_to_triangle_mesh(const PointRange& points, const NamedParameters& np) { typedef typename GetGeomTraits::type Geom_traits; + NamedParameters>::type GeomTraits; - return approximate_Hausdorff_distance + return approximate_Hausdorff_distance (points,tm,parameters::choose_parameter(parameters::get_parameter(np, internal_np::vertex_point), get_const_property_map(vertex_point, tm))); } @@ -1099,16 +1160,16 @@ double approximate_max_distance_to_point_set(const TriangleMesh& tm, const NamedParameters& np) { typedef typename GetGeomTraits::type Geom_traits; + NamedParameters>::type GeomTraits; typedef boost::graph_traits GT; - typedef Orthogonal_k_neighbor_search > Knn; + typedef Orthogonal_k_neighbor_search > Knn; typedef typename Knn::Tree Tree; Tree tree(points.begin(), points.end()); - CRefiner ref; + CRefiner ref; for(typename GT::face_descriptor f : faces(tm)) { - typename Geom_traits::Point_3 points[3]; + typename GeomTraits::Point_3 points[3]; typename GT::halfedge_descriptor hd(halfedge(f,tm)); for(int i=0; i<3; ++i) { 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 index 8a36f88cfe8..2a121e6771e 100644 --- a/Polygon_mesh_processing/test/Polygon_mesh_processing/test_pmp_distance.cpp +++ b/Polygon_mesh_processing/test/Polygon_mesh_processing/test_pmp_distance.cpp @@ -1,6 +1,7 @@ #include #include #include +#include #include @@ -258,6 +259,11 @@ void general_tests(const TriangleMesh& m1, std::cout << "Max distance to triangle mesh (sequential) " << PMP::max_distance_to_triangle_mesh(points,m1) << "\n"; + + std::vector samples; + PMP::sample_triangle_mesh(m1, std::back_inserter(samples)); + std::cout << samples.size()<<" points sampled on mesh."<> m1; input.close(); - input.open(argv[2]); input >> m2; - + input.close(); std::cout << "First mesh has " << num_faces(m1) << " faces\n"; std::cout << "Second mesh has " << num_faces(m2) << " faces\n"; @@ -305,5 +310,16 @@ int main(int, char** argv) test_concept(); + std::vector > faces; + std::vector points; + input.open(argv[1]); + CGAL::read_OFF(input, points, faces); + input.close(); + + std::vector samples; + PMP::sample_triangle_soup(points, faces, std::back_inserter(samples)); + std::cout<