From d7b0fcc4f60b004c92f35313bafb45605bacc828 Mon Sep 17 00:00:00 2001 From: Simon Giraudot Date: Tue, 23 Jul 2019 12:49:59 +0200 Subject: [PATCH] WIP: simplify/reorganize reconstruction plugin --- .../Plugins/Point_set/CMakeLists.txt | 2 +- ...ce_reconstruction_advancing_front_impl.cpp | 235 +++ .../Surface_reconstruction_plugin.cpp | 1269 ++--------------- .../Surface_reconstruction_plugin.ui | 627 ++++---- .../Surface_reconstruction_plugin_impl.cpp | 262 ---- .../Surface_reconstruction_poisson_impl.cpp | 238 ++++ ...urface_reconstruction_scale_space_impl.cpp | 228 +++ 7 files changed, 1049 insertions(+), 1812 deletions(-) create mode 100644 Polyhedron/demo/Polyhedron/Plugins/Point_set/Surface_reconstruction_advancing_front_impl.cpp delete mode 100644 Polyhedron/demo/Polyhedron/Plugins/Point_set/Surface_reconstruction_plugin_impl.cpp create mode 100644 Polyhedron/demo/Polyhedron/Plugins/Point_set/Surface_reconstruction_poisson_impl.cpp create mode 100644 Polyhedron/demo/Polyhedron/Plugins/Point_set/Surface_reconstruction_scale_space_impl.cpp diff --git a/Polyhedron/demo/Polyhedron/Plugins/Point_set/CMakeLists.txt b/Polyhedron/demo/Polyhedron/Plugins/Point_set/CMakeLists.txt index df97b278ef6..818ecbd1082 100644 --- a/Polyhedron/demo/Polyhedron/Plugins/Point_set/CMakeLists.txt +++ b/Polyhedron/demo/Polyhedron/Plugins/Point_set/CMakeLists.txt @@ -1,7 +1,7 @@ include( polyhedron_demo_macros ) if(EIGEN3_FOUND) qt5_wrap_ui( surface_reconstructionUI_FILES Surface_reconstruction_plugin.ui) - polyhedron_demo_plugin(surface_reconstruction_plugin Surface_reconstruction_plugin Surface_reconstruction_plugin_impl ${surface_reconstructionUI_FILES} KEYWORDS PointSetProcessing) + polyhedron_demo_plugin(surface_reconstruction_plugin Surface_reconstruction_plugin Surface_reconstruction_poisson_impl Surface_reconstruction_advancing_front_impl Surface_reconstruction_scale_space_impl ${surface_reconstructionUI_FILES} KEYWORDS PointSetProcessing) target_link_libraries(surface_reconstruction_plugin PUBLIC scene_polygon_soup_item scene_surface_mesh_item scene_points_with_normal_item) qt5_wrap_ui( point_set_normal_estimationUI_FILES Point_set_normal_estimation_plugin.ui) diff --git a/Polyhedron/demo/Polyhedron/Plugins/Point_set/Surface_reconstruction_advancing_front_impl.cpp b/Polyhedron/demo/Polyhedron/Plugins/Point_set/Surface_reconstruction_advancing_front_impl.cpp new file mode 100644 index 00000000000..58bfbeff846 --- /dev/null +++ b/Polyhedron/demo/Polyhedron/Plugins/Point_set/Surface_reconstruction_advancing_front_impl.cpp @@ -0,0 +1,235 @@ +#include +#include +#include + +#include "Kernel_type.h" +#include "SMesh_type.h" +#include "Scene_points_with_normal_item.h" + +typedef std::array Facet; +typedef CGAL::Point_set_with_structure Structuring; +typedef Kernel::Plane_3 Plane_3; + +struct Construct{ + typedef std::array Facet; + typedef typename boost::property_map::type VPmap; + SMesh& mesh; + VPmap vpmap; + std::vector::vertex_descriptor> vertices; + + template + Construct(SMesh& mesh, const PointRange& points) + : mesh(mesh) + { + vpmap = get(boost::vertex_point, mesh); + for (const auto& p : points) + { + typename boost::graph_traits::vertex_descriptor v; + v = add_vertex(mesh); + vertices.push_back(v); + put(vpmap, v, p); + } + } + + Construct& operator=(const Facet f) + { + typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; + std::vector facet; + facet.resize(3); + facet[0]=vertices[f[0]]; + facet[1]=vertices[f[1]]; + facet[2]=vertices[f[2]]; + CGAL::Euler::add_face(facet, mesh); + return *this; + } + Construct& + operator*() { return *this; } + Construct& + operator++() { return *this; } + Construct + operator++(int) { return *this; } +}; + +struct Radius { + + double bound; + + Radius(double bound) + : bound(bound) + {} + + template + double operator() (const AdvancingFront& adv, Cell_handle& c, + const int& index) const + { + // bound == 0 is better than bound < infinity + // as it avoids the distance computations + if(bound == 0){ + return adv.smallest_radius_delaunay_sphere (c, index); + } + + // If radius > bound, return infinity so that facet is not used + double d = 0; + d = sqrt(squared_distance(c->vertex((index+1)%4)->point(), + c->vertex((index+2)%4)->point())); + if(d>bound) return adv.infinity(); + d = sqrt(squared_distance(c->vertex((index+2)%4)->point(), + c->vertex((index+3)%4)->point())); + if(d>bound) return adv.infinity(); + d = sqrt(squared_distance(c->vertex((index+1)%4)->point(), + c->vertex((index+3)%4)->point())); + if(d>bound) return adv.infinity(); + + // Otherwise, return usual priority value: smallest radius of + // delaunay sphere + return adv.smallest_radius_delaunay_sphere (c, index); + } +}; + + +template +struct Priority_with_structure_coherence { + + Structuring& structuring; + double bound; + + Priority_with_structure_coherence(Structuring& structuring, + double bound) + : structuring (structuring), bound (bound) + {} + + template + double operator() (AdvancingFront& adv, Cell_handle& c, + const int& index) const + { + // If perimeter > bound, return infinity so that facet is not used + if (bound != 0) + { + double d = 0; + d = sqrt(squared_distance(c->vertex((index+1)%4)->point(), + c->vertex((index+2)%4)->point())); + if(d>bound) return adv.infinity(); + d += sqrt(squared_distance(c->vertex((index+2)%4)->point(), + c->vertex((index+3)%4)->point())); + if(d>bound) return adv.infinity(); + d += sqrt(squared_distance(c->vertex((index+1)%4)->point(), + c->vertex((index+3)%4)->point())); + if(d>bound) return adv.infinity(); + } + + Facet f = {{ c->vertex ((index + 1) % 4)->info (), + c->vertex ((index + 2) % 4)->info (), + c->vertex ((index + 3) % 4)->info () }}; + + double weight = 100. * (5 - structuring.facet_coherence (f)); + + return weight * adv.smallest_radius_delaunay_sphere (c, index); + } + +}; + +void get_planes_from_shape_map (const Point_set& points, + Point_set::Property_map shape_map, + std::vector& planes) +{ + std::vector sorted_indices; + sorted_indices.reserve (points.size()); + std::copy (points.begin(), points.end(), std::back_inserter (sorted_indices)); + + std::sort (sorted_indices.begin(), sorted_indices.end(), + [&](const Point_set::Index& a, const Point_set::Index& b) -> bool + { + return shape_map[a] < shape_map[b]; + }); + + std::size_t nb_planes = shape_map[sorted_indices.back()] + 1; + planes.reserve (nb_planes); + + std::vector::iterator begin = sorted_indices.end(); + int plane_idx = shape_map[sorted_indices.front()]; + if (plane_idx != -1) + begin = sorted_indices.begin(); + + for (std::vector::iterator it = sorted_indices.begin(); + it != sorted_indices.end(); ++ it) + { + if (shape_map[*it] != plane_idx) + { + if (begin != sorted_indices.end()) + { + Plane_3 plane; + CGAL::linear_least_squares_fitting_3 + (boost::make_transform_iterator + (begin, CGAL::Property_map_to_unary_function(points.point_map())), + boost::make_transform_iterator + (it, CGAL::Property_map_to_unary_function(points.point_map())), + plane, CGAL::Dimension_tag<0>()); + planes.push_back (plane); + } + begin = it; + plane_idx = shape_map[*it]; + } + } + Plane_3 plane; + CGAL::linear_least_squares_fitting_3 + (boost::make_transform_iterator + (begin, CGAL::Property_map_to_unary_function(points.point_map())), + boost::make_transform_iterator + (sorted_indices.end(), CGAL::Property_map_to_unary_function(points.point_map())), + plane, CGAL::Dimension_tag<0>()); + planes.push_back (plane); +} + +SMesh* advancing_front (const Point_set& points, + double longest_edge, + double radius_ratio_bound, + double beta, + bool structuring) +{ + SMesh* mesh = new SMesh; + + if (structuring) // todo + { + Point_set::Property_map shape_map + = points.property_map("shape").first; + + typedef CGAL::Point_set_with_structure Structuring; + std::vector planes; + + get_planes_from_shape_map (points, shape_map, planes); + Structuring structuring + (points, planes, + longest_edge, + points.parameters(). + plane_map(CGAL::Identity_property_map()). + plane_index_map(shape_map)); + + std::vector structured; + structured.reserve (structuring.size()); + for (std::size_t i = 0; i < structuring.size(); ++ i) + structured.push_back (structuring.point(i)); + + Priority_with_structure_coherence priority + (structuring, longest_edge); + Construct construct(*mesh, structured); + CGAL::advancing_front_surface_reconstruction(points.points().begin(), + points.points().end(), + construct, + priority, + radius_ratio_bound, + beta); + } + else + { + Radius filter(longest_edge); + Construct construct (*mesh, points.points()); + CGAL::advancing_front_surface_reconstruction(points.points().begin(), + points.points().end(), + construct, + filter, + radius_ratio_bound, + beta); + } + + return mesh; +} diff --git a/Polyhedron/demo/Polyhedron/Plugins/Point_set/Surface_reconstruction_plugin.cpp b/Polyhedron/demo/Polyhedron/Plugins/Point_set/Surface_reconstruction_plugin.cpp index f98cd56a243..783dc5a4e82 100644 --- a/Polyhedron/demo/Polyhedron/Plugins/Point_set/Surface_reconstruction_plugin.cpp +++ b/Polyhedron/demo/Polyhedron/Plugins/Point_set/Surface_reconstruction_plugin.cpp @@ -14,674 +14,40 @@ #include #include -#include -#include -#include #include #include -#include #include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#include -#include -#include #include "ui_Surface_reconstruction_plugin.h" #include "CGAL/Kernel_traits.h" -// Concurrency -#ifdef CGAL_LINKED_WITH_TBB -typedef CGAL::Parallel_tag Concurrency_tag; -#else -typedef CGAL::Sequential_tag Concurrency_tag; -#endif - typedef Kernel::Point_3 Point; typedef Kernel::Vector_3 Vector; -// types for K nearest neighbors search -typedef CGAL::Search_traits_3 Tree_traits; -typedef CGAL::Orthogonal_k_neighbor_search Neighbor_search; -typedef Neighbor_search::Tree Tree; -typedef Neighbor_search::iterator Search_iterator; - -typedef CGAL::Scale_space_surface_reconstruction_3 ScaleSpace; -typedef CGAL::Scale_space_reconstruction_3::Advancing_front_mesher ScaleSpaceAFM; -typedef CGAL::Scale_space_reconstruction_3::Alpha_shape_mesher ScaleSpaceASM; -typedef CGAL::Scale_space_reconstruction_3::Jet_smoother ScaleSpaceJS; -typedef CGAL::Scale_space_reconstruction_3::Weighted_PCA_smoother ScaleSpaceWPS; - -typedef std::array Facet; -template -struct Construct{ - typedef std::array Facet; - typedef typename Traits::Point_3 Point_3; - typedef typename boost::property_map::type VPmap; - Mesh& mesh; - VPmap vpmap; - std::vector::vertex_descriptor> vertices; - template < typename PointIterator> - Construct(Mesh& mesh,PointIterator b, PointIterator e) - : mesh(mesh) - { - vpmap = get(boost::vertex_point, mesh); - for(; b!=e; ++b){ - typename boost::graph_traits::vertex_descriptor v; - v = add_vertex(mesh); - vertices.push_back(v); - put(vpmap, v, *b); - } - } - - Construct& operator=(const Facet f) - { - typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; - std::vector facet; - facet.resize(3); - facet[0]=vertices[f[0]]; - facet[1]=vertices[f[1]]; - facet[2]=vertices[f[2]]; - CGAL::Euler::add_face(facet, mesh); - return *this; - } - Construct& - operator*() { return *this; } - Construct& - operator++() { return *this; } - Construct - operator++(int) { return *this; } -}; - - -// Poisson reconstruction method: -// Reconstructs a surface mesh from a point set and returns it. -SMesh* poisson_reconstruct_sm(Point_set& points, - Kernel::FT sm_angle, // Min triangle angle (degrees). - Kernel::FT sm_radius, // Max triangle size w.r.t. point set average spacing. - Kernel::FT sm_distance, // Approximation error w.r.t. point set average spacing. - const QString& solver_name, // solver name - bool use_two_passes, - bool do_not_fill_holes); - - -struct Radius { - - double bound; - - Radius(double bound) - : bound(bound) - {} - - template - double operator() (const AdvancingFront& adv, Cell_handle& c, - const int& index) const - { - // bound == 0 is better than bound < infinity - // as it avoids the distance computations - if(bound == 0){ - return adv.smallest_radius_delaunay_sphere (c, index); - } - - // If radius > bound, return infinity so that facet is not used - double d = 0; - d = sqrt(squared_distance(c->vertex((index+1)%4)->point(), - c->vertex((index+2)%4)->point())); - if(d>bound) return adv.infinity(); - d = sqrt(squared_distance(c->vertex((index+2)%4)->point(), - c->vertex((index+3)%4)->point())); - if(d>bound) return adv.infinity(); - d = sqrt(squared_distance(c->vertex((index+1)%4)->point(), - c->vertex((index+3)%4)->point())); - if(d>bound) return adv.infinity(); - - // Otherwise, return usual priority value: smallest radius of - // delaunay sphere - return adv.smallest_radius_delaunay_sphere (c, index); - } -}; - - -template -struct Priority_with_structure_coherence { - - Structuring& structuring; - double bound; - - Priority_with_structure_coherence(Structuring& structuring, - double bound) - : structuring (structuring), bound (bound) - {} - - template - double operator() (AdvancingFront& adv, Cell_handle& c, - const int& index) const - { - // If perimeter > bound, return infinity so that facet is not used - if (bound != 0) - { - double d = 0; - d = sqrt(squared_distance(c->vertex((index+1)%4)->point(), - c->vertex((index+2)%4)->point())); - if(d>bound) return adv.infinity(); - d += sqrt(squared_distance(c->vertex((index+2)%4)->point(), - c->vertex((index+3)%4)->point())); - if(d>bound) return adv.infinity(); - d += sqrt(squared_distance(c->vertex((index+1)%4)->point(), - c->vertex((index+3)%4)->point())); - if(d>bound) return adv.infinity(); - } - - Facet f = {{ c->vertex ((index + 1) % 4)->info (), - c->vertex ((index + 2) % 4)->info (), - c->vertex ((index + 3) % 4)->info () }}; - - double weight = 100. * (5 - structuring.facet_coherence (f)); - - return weight * adv.smallest_radius_delaunay_sphere (c, index); - } - -}; - - -struct On_the_fly_pair{ - const Point_set& points; - typedef std::pair result_type; - - On_the_fly_pair(const Point_set& points) : points(points) {} - - result_type - operator()(const Point_set::Index& i) const - { - return result_type(points.point(i), (std::size_t)i); - } - -}; - - -namespace SurfaceReconstruction -{ - typedef Kernel::Point_3 Point; - typedef Kernel::Vector_3 Vector; - // types for K nearest neighbors search - typedef CGAL::Search_traits_3 SearchTraits_3; - typedef CGAL::Search_traits_adapter Search_traits; - typedef CGAL::Orthogonal_k_neighbor_search Neighbor_search; - typedef Neighbor_search::Tree Tree; - typedef Neighbor_search::Distance Distance; - typedef Neighbor_search::iterator Search_iterator; - - template - void generate_scales (OutputIterator out, - const unsigned int scale_min = 6, - const double factor = 1.15, - unsigned int nb_scales = 30) - { - unsigned int prev = -1; - - for (unsigned int i = 0; i < nb_scales; ++ i) - { - unsigned int current = static_cast(scale_min * std::pow (factor, i)); - if (current != prev) - { - *(out ++) = current; - prev = current; - } - else - ++ nb_scales; - } - } - - unsigned int scale_of_anisotropy (const Point_set& points, double& size) - { - Tree tree(points.begin_or_selection_begin(), points.end(), - Tree::Splitter(), Search_traits (points.point_map())); - - double ratio_kept = (points.size() < 1000) - ? 1. : 1000. / (points.size()); - - std::vector subset; - for (Point_set::const_iterator it = points.begin(); it != points.end(); ++ it) - if (CGAL::get_default_random().get_double() < ratio_kept) - subset.push_back (points.point(*it)); - - std::vector scales; - generate_scales (std::back_inserter (scales)); - - std::vector chosen; - Distance tr_dist (points.point_map()); - - for (std::size_t i = 0; i < subset.size (); ++ i) - { - Neighbor_search search(tree, subset[i],scales.back(), 0, true, tr_dist); - double current = 0.; - unsigned int nb = 0; - std::size_t index = 0; - double maximum = 0.; - unsigned int c = 0; - - for (Search_iterator search_iterator = search.begin(); - search_iterator != search.end (); ++ search_iterator, ++ nb) - { - current += search_iterator->second; - - if (nb + 1 == scales[index]) - { - double score = std::sqrt (current / scales[index]) - / std::pow (scales[index], 0.75); // NB ^ (3/4) - - if (score > maximum) - { - maximum = score; - c = scales[index]; - } - - ++ index; - if (index == scales.size ()) - break; - } - } - chosen.push_back (c); - } - - double mean = 0.; - for (std::size_t i = 0; i < chosen.size(); ++ i) - mean += chosen[i]; - mean /= chosen.size(); - unsigned int aniso_scale = static_cast(mean); - - size = 0.; - for (std::size_t i = 0; i < subset.size (); ++ i) - { - Neighbor_search search(tree, subset[i], aniso_scale, 0, true, tr_dist); - size += std::sqrt ((-- search.end())->second); - } - size /= subset.size(); - - return aniso_scale; - } - - - unsigned int scale_of_noise (const Point_set& points, double& size) - { - Tree tree(points.begin_or_selection_begin(), points.end(), - Tree::Splitter(), Search_traits (points.point_map())); - Distance tr_dist (points.point_map()); - - double ratio_kept = (points.size() < 1000) - ? 1. : 1000. / (points.size()); - - std::vector subset; - for (Point_set::const_iterator it = points.begin(); it != points.end(); ++ it) - if (CGAL::get_default_random().get_double() < ratio_kept) - subset.push_back (points.point(*it)); - - std::vector scales; - generate_scales (std::back_inserter (scales)); - - std::vector chosen; - - for (std::size_t i = 0; i < subset.size (); ++ i) - { - Neighbor_search search(tree, subset[i],scales.back(), 0, true, tr_dist); - double current = 0.; - unsigned int nb = 0; - std::size_t index = 0; - double minimum = (std::numeric_limits::max)(); - unsigned int c = 0; - - for (Search_iterator search_iterator = search.begin(); - search_iterator != search.end (); ++ search_iterator, ++ nb) - { - current += search_iterator->second; - - if (nb + 1 == scales[index]) - { - double score = std::sqrt (current / scales[index]) - / std::pow (scales[index], 0.375); // NB ^ (5/12) - - if (score < minimum) - { - minimum = score; - c = scales[index]; - } - - ++ index; - if (index == scales.size ()) - break; - } - } - chosen.push_back (c); - } - - std::sort (chosen.begin (), chosen.end()); - - unsigned int noise_scale = chosen[chosen.size() / 2]; - size = 0.; - for (std::size_t i = 0; i < subset.size (); ++ i) - { - Neighbor_search search(tree, subset[i], noise_scale, 0, true, tr_dist); - size += std::sqrt ((-- search.end())->second); - } - size /= subset.size(); - - - return noise_scale; - } - - void simplify_point_set (Point_set& points, double size) - { - points.set_first_selected (CGAL::grid_simplify_point_set (points, size)); - points.delete_selection(); - } - - void smooth_point_set (Point_set& points, unsigned int scale) - { - CGAL::jet_smooth_point_set(points, scale); - } - - template - void scale_space (const Point_set& points, ItemsInserter items, - bool jet_smoother, - unsigned int iterations, - unsigned int neighbors, unsigned int fitting, unsigned int monge, - unsigned int neighborhood_size, unsigned int samples, - bool advancing_front_mesher, - bool generate_smooth, - double longest_edge, double radius_ratio_bound, double beta_angle, - bool separate_shells, bool force_manifold) - { - ScaleSpace reconstruct (points.points().begin(), points.points().end()); - - double squared_radius = 0.; - if (jet_smoother) - { - ScaleSpaceJS smoother(neighbors, fitting, monge); - reconstruct.increase_scale(iterations, smoother); - if (!advancing_front_mesher) - squared_radius = CGAL::compute_average_spacing (points, neighbors); - } - else - { - ScaleSpaceWPS smoother(neighborhood_size, samples); - reconstruct.increase_scale(iterations, smoother); - squared_radius = smoother.squared_radius(); - - } - - if (advancing_front_mesher) - { - ScaleSpaceAFM mesher (longest_edge, radius_ratio_bound, beta_angle); - reconstruct.reconstruct_surface (mesher); - - Scene_polygon_soup_item* new_item - = new Scene_polygon_soup_item (); - new_item->setColor(Qt::lightGray); - new_item->setRenderingMode(FlatPlusEdges); - new_item->init_polygon_soup(points.size(), reconstruct.number_of_facets ()); - - Scene_polygon_soup_item* smooth_item = NULL; - if (generate_smooth) - { - smooth_item = new Scene_polygon_soup_item (); - smooth_item->setColor(Qt::lightGray); - smooth_item->setRenderingMode(FlatPlusEdges); - smooth_item->init_polygon_soup(points.size(), reconstruct.number_of_facets ()); - } - - std::map map_i2i; - std::size_t current_index = 0; - - for (ScaleSpace::Facet_iterator it = reconstruct.facets_begin(); - it != reconstruct.facets_end(); ++ it) - { - for (unsigned int ind = 0; ind < 3; ++ ind) - { - if (map_i2i.find ((*it)[ind]) == map_i2i.end ()) - { - map_i2i.insert (std::make_pair ((*it)[ind], current_index ++)); - Point p = points.point(*(points.begin_or_selection_begin() + (*it)[ind])); - new_item->new_vertex (p.x (), p.y (), p.z ()); - - if (generate_smooth) - { - p = *(reconstruct.points_begin() + (*it)[ind]); - smooth_item->new_vertex (p.x (), p.y (), p.z ()); - } - } - } - new_item->new_triangle( map_i2i[(*it)[0]], - map_i2i[(*it)[1]], - map_i2i[(*it)[2]] ); - if (generate_smooth) - smooth_item->new_triangle( map_i2i[(*it)[0]], - map_i2i[(*it)[1]], - map_i2i[(*it)[2]] ); - - } - - *(items ++) = new_item; - if (generate_smooth) - *(items ++) = smooth_item; - } - else - { - ScaleSpaceASM mesher (squared_radius, separate_shells, force_manifold); - reconstruct.reconstruct_surface (mesher); - - for( unsigned int sh = 0; sh < mesher.number_of_shells(); ++sh ) - { - Scene_polygon_soup_item* new_item - = new Scene_polygon_soup_item (); - new_item->setColor(Qt::lightGray); - new_item->setRenderingMode(FlatPlusEdges); - new_item->init_polygon_soup(points.size(), mesher.number_of_triangles ()); - - Scene_polygon_soup_item* smooth_item = NULL; - if (generate_smooth) - { - smooth_item = new Scene_polygon_soup_item (); - smooth_item->setColor(Qt::lightGray); - smooth_item->setRenderingMode(FlatPlusEdges); - smooth_item->init_polygon_soup(points.size(), mesher.number_of_triangles ()); - } - - std::map map_i2i; - unsigned int current_index = 0; - - for (ScaleSpaceASM::Facet_iterator it = mesher.shell_begin (sh); - it != mesher.shell_end (sh); ++ it) - { - for (unsigned int ind = 0; ind < 3; ++ ind) - { - if (map_i2i.find ((*it)[ind]) == map_i2i.end ()) - { - map_i2i.insert (std::make_pair ((*it)[ind], current_index ++)); - Point p = points.point(*(points.begin_or_selection_begin() + (*it)[ind])); - new_item->new_vertex (p.x (), p.y (), p.z ()); - - if (generate_smooth) - { - p = *(reconstruct.points_begin() + (*it)[ind]); - smooth_item->new_vertex (p.x (), p.y (), p.z ()); - } - } - } - new_item->new_triangle( map_i2i[(*it)[0]], - map_i2i[(*it)[1]], - map_i2i[(*it)[2]] ); - if (generate_smooth) - smooth_item->new_triangle( map_i2i[(*it)[0]], - map_i2i[(*it)[1]], - map_i2i[(*it)[2]] ); - - } - - *(items ++) = new_item; - if (generate_smooth) - *(items ++) = smooth_item; - } - - if (force_manifold) - { - std::ptrdiff_t num = std::distance( mesher.garbage_begin( ), - mesher.garbage_end( ) ); - - Scene_polygon_soup_item* new_item - = new Scene_polygon_soup_item (); - new_item->setColor(Qt::blue); - new_item->setRenderingMode(FlatPlusEdges); - new_item->init_polygon_soup(points.size(), num); - - Scene_polygon_soup_item* smooth_item = NULL; - if (generate_smooth) - { - smooth_item = new Scene_polygon_soup_item (); - smooth_item->setColor(Qt::blue); - smooth_item->setRenderingMode(FlatPlusEdges); - smooth_item->init_polygon_soup(points.size(), num); - } - - std::map map_i2i; - - std::size_t current_index = 0; - for (ScaleSpaceASM::Facet_iterator it=mesher.garbage_begin(), - end=mesher.garbage_end();it!=end;++it) - { - for (unsigned int ind = 0; ind < 3; ++ ind) - { - if (map_i2i.find ((*it)[ind]) == map_i2i.end ()) - { - map_i2i.insert (std::make_pair ((*it)[ind], current_index ++)); - Point p = points.point(*(points.begin_or_selection_begin() + (*it)[ind])); - new_item->new_vertex (p.x (), p.y (), p.z ()); - - if (generate_smooth) - { - p = *(reconstruct.points_begin() + (*it)[ind]); - smooth_item->new_vertex (p.x (), p.y (), p.z ()); - } - } - - } - new_item->new_triangle( map_i2i[(*it)[0]], - map_i2i[(*it)[1]], - map_i2i[(*it)[2]] ); - if (generate_smooth) - smooth_item->new_triangle( map_i2i[(*it)[0]], - map_i2i[(*it)[1]], - map_i2i[(*it)[2]] ); - } - - *(items ++) = new_item; - if (generate_smooth) - *(items ++) = smooth_item; - - } - } - } - - struct Point_set_make_pair_point_index - : public CGAL::cpp98::unary_function > - { - const Point_set& point_set; - Point_set_make_pair_point_index (const Point_set& point_set) : point_set (point_set) { } - std::pair operator() (const Point_set::Index& i) const - { - return std::make_pair (point_set.point (i), i); - } - }; - - - struct Radius { - - double bound; - - Radius(double bound) - : bound(bound) - {} - - template - double operator() (const AdvancingFront& adv, Cell_handle& c, - const int& index) const - { - // bound == 0 is better than bound < infinity - // as it avoids the distance computations - if(bound == 0){ - return adv.smallest_radius_delaunay_sphere (c, index); - } - - // If radius > bound, return infinity so that facet is not used - double d = 0; - d = sqrt(squared_distance(c->vertex((index+1)%4)->point(), - c->vertex((index+2)%4)->point())); - if(d>bound) return adv.infinity(); - d = sqrt(squared_distance(c->vertex((index+2)%4)->point(), - c->vertex((index+3)%4)->point())); - if(d>bound) return adv.infinity(); - d = sqrt(squared_distance(c->vertex((index+1)%4)->point(), - c->vertex((index+3)%4)->point())); - if(d>bound) return adv.infinity(); - - // Otherwise, return usual priority value: smallest radius of - // delaunay sphere - return adv.smallest_radius_delaunay_sphere (c, index); - } - }; - - template - void advancing_front (const Point_set& points, FaceGraphItem* new_item, double size, - double radius_ratio_bound = 5., double beta = 0.52) - { - - // TODO: build DT with indices - typedef typename FaceGraphItem::Face_graph FaceGraph; - typedef typename CGAL::Kernel_traits::type>::value_type>::Kernel Traits; - - FaceGraph& P = * const_cast(new_item->face_graph()); - Radius filter(size); - Construct construct(P,points.points().begin(),points.points().end()); - CGAL::advancing_front_surface_reconstruction(points.points().begin(), - points.points().end(), - construct, - filter, - radius_ratio_bound, - beta); - - } - - void compute_normals (Point_set& points, unsigned int neighbors) - { - CGAL::jet_estimate_normals(points, 2 * neighbors); - - points.set_first_selected (CGAL::mst_orient_normals (points, 2 * neighbors)); - points.delete_selection(); - } - - struct build_from_pair - { - Point_set& m_pts; - - build_from_pair (Point_set& pts) : m_pts (pts) { } - - void operator() (const std::pair& pair) - { - m_pts.insert (pair.first, pair.second); - } - - - }; -} - +SMesh* advancing_front (const Point_set& points, + double longest_edge, + double radius_ratio_bound, + double beta, + bool structuring); + +SMesh* poisson_reconstruct (Point_set& points, + Kernel::FT sm_angle, // Min triangle angle (degrees). + Kernel::FT sm_radius, // Max triangle size w.r.t. point set average spacing. + Kernel::FT sm_distance, // Approximation error w.r.t. point set average spacing. + bool conjugate_gradient, + bool use_two_passes, + bool do_not_fill_holes); + +void scale_space (const Point_set& points, + std::vector& items, + bool jet_smoother, + unsigned int iterations, + unsigned int neighbors, unsigned int fitting, unsigned int monge, + unsigned int neighborhood_size, unsigned int samples, + bool advancing_front_mesher, + bool generate_smooth, + double longest_edge, double radius_ratio_bound, double beta_angle, + bool separate_shells, bool force_manifold); class Polyhedron_demo_surface_reconstruction_plugin_dialog : public QDialog, private Ui::SurfaceReconstructionDialog { @@ -690,26 +56,30 @@ public: Polyhedron_demo_surface_reconstruction_plugin_dialog(QWidget* /*parent*/ = 0) { setupUi(this); - -#ifdef CGAL_EIGEN3_ENABLED - m_inputSolver->addItem("Eigen - built-in CG"); - m_inputSolver->addItem("Eigen - built-in simplicial LDLt"); -#endif } unsigned int method () const { return tabWidget->currentIndex(); } - - // Auto - bool boundaries () const { return m_boundaries->isChecked (); } - bool interpolate () const { return m_interpolate->isChecked (); } + void disable_poisson() + { + tabWidget->setTabEnabled(1, false); + tabWidget->setTabText(1, QString("Poisson (requires oriented normals)")); + } + + void enable_structuring() + { + m_use_structuring->setEnabled(true); + m_use_structuring->setText(QString("Use Point Set Structuring")); + } + // Advancing front double longest_edge () const { return m_longestEdge->value (); } double radius_ratio_bound () const { return m_radiusRatioBound->value (); } double beta_angle () const { return m_betaAngle->value (); } + bool structuring() const { return m_use_structuring->isChecked(); } // Scale Space bool scalespace_js() const { return m_scalespace_jet->isChecked(); } @@ -727,27 +97,16 @@ public: bool separate_shells () const { return m_genShells->isChecked (); } bool force_manifold () const { return m_forceManifold->isChecked (); } - // Poisson double angle () const { return m_inputAngle->value (); } double radius () const { return m_inputRadius->value (); } double distance () const { return m_inputDistance->value (); } + bool conjugate_gradient() const { return m_conjugate_gradient->isChecked(); } bool two_passes () const { return m_inputTwoPasses->isChecked (); } bool do_not_fill_holes () const { return m_doNotFillHoles->isChecked (); } - - // RANSAC - bool region_growing () const { return m_regionGrowing->isChecked (); } - double connectivity_tolerance () const { return m_connectivityTolerance->value (); } - double noise_tolerance () const { return m_noiseTolerance->value (); } - double normal_tolerance () const { return m_normalTolerance->value (); } - unsigned int min_size_subset () const { return m_minSizeSubset->value (); } - bool generate_structured () const { return m_generateStructured->isChecked (); } - QString solver () const { return m_inputSolver->currentText (); } }; -#include - class Polyhedron_demo_surface_reconstruction_plugin : public QObject, public CGAL::Three::Polyhedron_demo_plugin_helper @@ -768,11 +127,9 @@ public: } - void automatic_reconstruction (const Polyhedron_demo_surface_reconstruction_plugin_dialog& dialog); void advancing_front_reconstruction (const Polyhedron_demo_surface_reconstruction_plugin_dialog& dialog); void scale_space_reconstruction (const Polyhedron_demo_surface_reconstruction_plugin_dialog& dialog); void poisson_reconstruction (const Polyhedron_demo_surface_reconstruction_plugin_dialog& dialog); - void ransac_reconstruction (const Polyhedron_demo_surface_reconstruction_plugin_dialog& dialog); //! Applicate for Point_sets with normals. bool applicable(QAction*) const { @@ -785,265 +142,6 @@ public: private: - void region_growing_reconstruction_impl (const Polyhedron_demo_surface_reconstruction_plugin_dialog& dialog) { - - using Point_map = typename Point_set::Point_map; - using Normal_map = typename Point_set::Vector_map; - - using Neighbor_query = - CGAL::Shape_detection::Point_set::Sphere_neighbor_query; - using Region_type = - CGAL::Shape_detection::Point_set::Least_squares_plane_fit_region; - using Region_growing = - CGAL::Shape_detection::Region_growing; - - const CGAL::Three::Scene_interface::Item_id index = scene->mainSelectionIndex(); - Scene_points_with_normal_item* point_set_item = - qobject_cast(scene->item(index)); - - if (point_set_item) { - - // Get a point set. - Point_set* points = point_set_item->point_set(); - if(!points) return; - QApplication::setOverrideCursor(Qt::WaitCursor); - - CGAL::Timer global_timer; - global_timer.start(); - CGAL::Timer local_timer; - - // Estimate normals. - if (!(point_set_item->has_normals())) { - local_timer.start(); - - std::cerr << "Estimation of normal vectors... "; - points->add_normal_map(); - CGAL::jet_estimate_normals(*points, 12); - local_timer.stop(); - point_set_item->setRenderingMode(PointsPlusNormals); - - std::cerr << "done in " << local_timer.time() << " second(s)" << std::endl; - local_timer.reset(); - } - - local_timer.start(); - - // Set parameters. - const double search_sphere_radius = - dialog.connectivity_tolerance(); - const double max_distance_to_plane = - dialog.noise_tolerance(); - const double max_accepted_angle = - (std::acos(0.7) * 180.0) / CGAL_PI; - const std::size_t min_region_size = - dialog.min_size_subset(); - - // Region growing. - Neighbor_query neighbor_query( - *points, - search_sphere_radius, - points->point_map()); - - Region_type region_type( - *points, - max_distance_to_plane, max_accepted_angle, min_region_size, - points->point_map(), points->normal_map()); - - Region_growing region_growing( - *points, neighbor_query, region_type); - - std::vector< std::vector > regions; - region_growing.detect(std::back_inserter(regions)); - - local_timer.stop(); - std::cout << regions.size() << " plane(s) found in " - << local_timer.time() << " second(s)" << std::endl; - local_timer.reset(); - - std::cout << "Structuring point set... " << std::endl; - using Plane_3 = typename Kernel::Plane_3; - typedef CGAL::Point_set_with_structure Structuring; - std::vector planes; - CGAL::Shape_detection::internal::create_planes_from_points( - *points, points->point_map(), regions, planes); - - local_timer.start(); - Structuring structuring( - *points, - planes, - search_sphere_radius, - points->parameters(). - plane_map(CGAL::Identity_property_map()). - plane_index_map(CGAL::Shape_detection::RG::Point_to_shape_index_map(*points, regions))); - - Scene_points_with_normal_item *structured = new Scene_points_with_normal_item; - structured->point_set()->add_normal_map(); - for (std::size_t i = 0; i < structuring.size(); ++i) - structured->point_set()->insert(structuring.point(i), structuring.normal(i)); - - local_timer.stop(); - std::cerr << structured->point_set()->size() << " point(s) generated in " - << local_timer.time() << std::endl; - local_timer.reset(); - - std::cerr << "Reconstructing... "; - local_timer.start(); - - Priority_with_structure_coherence priority( - structuring, 10.0 * search_sphere_radius); - - Scene_surface_mesh_item* reco_item = new Scene_surface_mesh_item(SMesh()); - SMesh& P = * const_cast(reco_item->polyhedron()); - Construct construct(P, - structured->point_set()->points().begin(), - structured->point_set()->points().end()); - CGAL::advancing_front_surface_reconstruction( - structured->point_set()->points().begin(), - structured->point_set()->points().end(), - construct, - priority, - 5.0, - 0.52); - - local_timer.stop(); - std::cerr << "done in " << local_timer.time() << " second(s)" << std::endl; - reco_item->setName(tr("%1 (Region-Growing-based reconstruction)").arg(scene->item(index)->name())); - reco_item->setColor(Qt::magenta); - reco_item->setRenderingMode(FlatPlusEdges); - reco_item->invalidateOpenGLBuffers(); - scene->addItem(reco_item); - - if (dialog.generate_structured()) { - structured->setName(tr("%1 (structured)").arg(point_set_item->name())); - structured->setRenderingMode(PointsPlusNormals); - structured->setColor(Qt::blue); - scene->addItem (structured); - - } else delete structured; - - std::cerr << "All done in " << global_timer.time() << " seconds." << std::endl; - QApplication::restoreOverrideCursor(); - } - } - - void ransac_reconstruction_impl (const Polyhedron_demo_surface_reconstruction_plugin_dialog& dialog) - { - - typedef Point_set::Point_map PointMap; - typedef Point_set::Vector_map NormalMap; - typedef CGAL::Shape_detection::Efficient_RANSAC_traits Traits; - typedef CGAL::Shape_detection::Efficient_RANSAC Efficient_RANSAC; - - const CGAL::Three::Scene_interface::Item_id index = scene->mainSelectionIndex(); - - Scene_points_with_normal_item* point_set_item = - qobject_cast(scene->item(index)); - - if(point_set_item) - { - // Gets point set - Point_set* points = point_set_item->point_set(); - if(!points) return; - - QApplication::setOverrideCursor(Qt::WaitCursor); - - CGAL::Timer global_timer; - global_timer.start(); - - CGAL::Timer local_timer; - - if (!(point_set_item->has_normals())) - { - local_timer.start(); - - std::cerr << "Estimation of normal vectors... "; - points->add_normal_map(); - CGAL::jet_estimate_normals(*points, 12); - local_timer.stop(); - point_set_item->setRenderingMode(PointsPlusNormals); - - std::cerr << "done in " << local_timer.time() << " second(s)" << std::endl; - local_timer.reset(); - } - - local_timer.start(); - Efficient_RANSAC ransac; - ransac.set_input(*points, points->point_map(), points->normal_map()); - - ransac.template add_shape_factory >(); - - typename Efficient_RANSAC::Parameters op; - op.min_points = dialog.min_size_subset(); - op.epsilon = dialog.noise_tolerance(); - op.cluster_epsilon = dialog.connectivity_tolerance(); - op.normal_threshold = 0.7; - - ransac.detect(op); - local_timer.stop(); - std::cout << ransac.shapes().size() << " plane(s) found in " - << local_timer.time() << " second(s)" << std::endl; - local_timer.reset(); - - std::cout << "Structuring point set... " << std::endl; - typedef CGAL::Point_set_with_structure Structuring; - typename Efficient_RANSAC::Plane_range planes = ransac.planes(); - - local_timer.start(); - Structuring structuring (*points, - planes, - op.cluster_epsilon, - points->parameters(). - plane_map(CGAL::Shape_detection::Plane_map()). - plane_index_map(CGAL::Shape_detection::Point_to_shape_index_map(*points, planes))); - - - Scene_points_with_normal_item *structured = new Scene_points_with_normal_item; - structured->point_set()->add_normal_map(); - for (std::size_t i = 0; i < structuring.size(); ++ i) - structured->point_set()->insert (structuring.point(i), structuring.normal(i)); - - local_timer.stop (); - std::cerr << structured->point_set()->size() << " point(s) generated in " - << local_timer.time() << std::endl; - local_timer.reset(); - - std::cerr << "Reconstructing... "; - local_timer.start(); - - Priority_with_structure_coherence priority (structuring, 10. * op.cluster_epsilon); - Scene_surface_mesh_item* reco_item = new Scene_surface_mesh_item(SMesh()); - SMesh& P = * const_cast(reco_item->polyhedron()); - Construct construct(P,structured->point_set()->points().begin(),structured->point_set()->points().end()); - CGAL::advancing_front_surface_reconstruction(structured->point_set()->points().begin(), - structured->point_set()->points().end(), - construct, - priority, - 5., - 0.52); - local_timer.stop(); - std::cerr << "done in " << local_timer.time() << " second(s)" << std::endl; - reco_item->setName(tr("%1 (RANSAC-based reconstruction)").arg(scene->item(index)->name())); - reco_item->setColor(Qt::magenta); - reco_item->setRenderingMode(FlatPlusEdges); - reco_item->invalidateOpenGLBuffers(); - scene->addItem(reco_item); - - if (dialog.generate_structured ()) - { - structured->setName(tr("%1 (structured)").arg(point_set_item->name())); - structured->setRenderingMode(PointsPlusNormals); - structured->setColor(Qt::blue); - scene->addItem (structured); - } - else - delete structured; - - std::cerr << "All done in " << global_timer.time() << " seconds." << std::endl; - - QApplication::restoreOverrideCursor(); - } - } - public Q_SLOTS: void on_actionSurfaceReconstruction_triggered(); }; // end class Polyhedron_surface_reconstruction_plugin @@ -1061,8 +159,14 @@ void Polyhedron_demo_surface_reconstruction_plugin::on_actionSurfaceReconstructi //generate the dialog box to set the options Polyhedron_demo_surface_reconstruction_plugin_dialog dialog; dialog.setWindowFlags(Qt::Dialog|Qt::CustomizeWindowHint|Qt::WindowCloseButtonHint); + + if (!pts_item->point_set()->has_normal_map()) + dialog.disable_poisson(); + if (pts_item->point_set()->has_property_map ("shape")) + dialog.enable_structuring(); + if(!dialog.exec()) - return; + return; unsigned int method = dialog.method (); switch (method) @@ -1076,12 +180,6 @@ void Polyhedron_demo_surface_reconstruction_plugin::on_actionSurfaceReconstructi case 2: scale_space_reconstruction (dialog); break; - case 3: - ransac_reconstruction (dialog); - break; - case 4: - automatic_reconstruction (dialog); - break; default: std::cerr << "Error: unkown method." << std::endl; return; @@ -1091,191 +189,6 @@ void Polyhedron_demo_surface_reconstruction_plugin::on_actionSurfaceReconstructi } } -void Polyhedron_demo_surface_reconstruction_plugin::automatic_reconstruction -(const Polyhedron_demo_surface_reconstruction_plugin_dialog& dialog) -{ - const CGAL::Three::Scene_interface::Item_id index = scene->mainSelectionIndex(); - - Scene_points_with_normal_item* pts_item = - qobject_cast(scene->item(index)); - - if(pts_item) - { - // Gets point set - Point_set* points = pts_item->point_set(); - // wait cursor - QApplication::setOverrideCursor(Qt::WaitCursor); - - QTime time; - std::cout << "Meta surface reconstruction with the following requirements:" << std::endl - << (dialog.boundaries() ? " * Output shape has boundaries" : " * Output shape is closed") << std::endl - << (dialog.interpolate() ? " * Output shape passes through input points" - : " * Output shape approximates input points") << std::endl; - - Scene_points_with_normal_item* new_item = NULL; - if (!(dialog.interpolate())) - { - new_item = new Scene_points_with_normal_item(*pts_item); - new_item->setName(QString("%1 (preprocessed)").arg(pts_item->name())); - new_item->resetSelection(); - new_item->invalidateOpenGLBuffers(); - - points = new_item->point_set(); - } - - std::cerr << "Analysing isotropy of point set... "; - time.start(); - double aniso_size; - unsigned int aniso_scale = SurfaceReconstruction::scale_of_anisotropy (*points, aniso_size); - std::cerr << "ok (" << time.elapsed() << " ms)" << std::endl; - - std::cerr << " -> Scale / size = " << aniso_scale << " / " << aniso_size << std::endl; - - bool isotropic = (aniso_scale == 6); - std::cerr << (isotropic ? " -> Point set is isotropic" : " -> Point set is anisotropic") << std::endl; - if (!(dialog.interpolate()) && !isotropic) - { - std::cerr << "Correcting anisotropy of point set... "; - time.restart(); - std::size_t prev_size = points->size (); - SurfaceReconstruction::simplify_point_set (*points, aniso_size); - std::cerr << "ok (" << time.elapsed() << " ms)" << std::endl; - std::cerr << " -> " << prev_size - points->size() << " point(s) removed (" - << 100. * (prev_size - points->size()) / (double)(prev_size) - << "%)" << std::endl; - } - - std::cerr << "Analysing noise of point set... "; - time.restart(); - double noise_size; - unsigned int noise_scale = SurfaceReconstruction::scale_of_noise (*points, noise_size); - std::cerr << "ok (" << time.elapsed() << " ms)" << std::endl; - - std::cerr << " -> Scale / size = " << noise_scale << " / " << noise_size << std::endl; - - bool noisy = (noise_scale > 6); - std::cerr << (noisy ? " -> Point set is noisy" : " -> Point set is noise-free") << std::endl; - - if (!(dialog.interpolate()) && noisy) - { - std::cerr << "Denoising point set... "; - time.restart(); - SurfaceReconstruction::smooth_point_set (*points, noise_scale); - new_item->point_set()->remove_normal_map(); - std::cerr << "ok (" << time.elapsed() << " ms)" << std::endl; - } - - if (dialog.interpolate()) - { - if (noisy) - { - std::cerr << "Scale space reconstruction... "; - time.restart(); - - std::vector reco_items; - - SurfaceReconstruction::scale_space (*points, std::back_inserter (reco_items), - true, - 4, - (std::max)(noise_scale, aniso_scale), 2, 2, - 0, 0, - true, - false, - 10. * (std::max)(noise_size, aniso_size), 5., 0.52, - false, false); - - for (std::size_t i = 0; i < reco_items.size (); ++ i) - { - reco_items[i]->setName(tr("%1 (scale space)").arg(scene->item(index)->name())); - scene->addItem (reco_items[i]); - } - - std::cerr << "ok (" << time.elapsed() << " ms)" << std::endl; - } - else - { - std::cerr << "Advancing front reconstruction... "; - time.restart(); - - Scene_surface_mesh_item* reco_item = new Scene_surface_mesh_item(SMesh()); - SurfaceReconstruction::advancing_front (*points, reco_item, 10. * (std::max)(noise_size, aniso_size)); - - reco_item->setName(tr("%1 (advancing front)").arg(scene->item(index)->name())); - reco_item->setColor(Qt::lightGray); - reco_item->setRenderingMode(FlatPlusEdges); - reco_item->invalidateOpenGLBuffers(); - scene->addItem(reco_item); - std::cerr << "ok (" << time.elapsed() << " ms)" << std::endl; - } - } - else - { - if (dialog.boundaries()) - { - std::cerr << "Advancing front reconstruction... "; - time.restart(); - Scene_surface_mesh_item* reco_item = new Scene_surface_mesh_item(SMesh()); - SurfaceReconstruction::advancing_front (*points, reco_item, 10. * (std::max)(noise_size, aniso_size)); - - reco_item->setName(tr("%1 (advancing front)").arg(scene->item(index)->name())); - reco_item->setColor(Qt::lightGray); - reco_item->setRenderingMode(FlatPlusEdges); - reco_item->invalidateOpenGLBuffers(); - scene->addItem(reco_item); - std::cerr << "ok (" << time.elapsed() << " ms)" << std::endl; - } - else - { - if (!(new_item->has_normals())) - { - std::cerr << "Estimation of normal vectors... "; - time.restart(); - - SurfaceReconstruction::compute_normals (*points, noise_scale); - - new_item->point_set()->add_normal_map(); - new_item->setRenderingMode(PointsPlusNormals); - - std::cerr << "ok (" << time.elapsed() << " ms)" << std::endl; - } - - std::cerr << "Poisson reconstruction... "; - time.restart(); - SMesh* smRemesh = - poisson_reconstruct_sm(*points, - 20, - 100 * (std::max)(noise_size, aniso_size), - (std::max)(noise_size, aniso_size), - QString ("Eigen - built-in CG"), false, false); - if(smRemesh) - { - // Add polyhedron to scene - Scene_surface_mesh_item* reco_item = new Scene_surface_mesh_item(smRemesh); - reco_item->setName(tr("%1 (poisson)").arg(pts_item->name())); - reco_item->setColor(Qt::lightGray); - reco_item->invalidateOpenGLBuffers(); - scene->addItem(reco_item); - } - - - std::cerr << "ok (" << time.elapsed() << " ms)" << std::endl; - } - } - - if (!(dialog.interpolate())) - { - if (noisy || !isotropic) - scene->addItem(new_item); - else - delete new_item; - } - - // default cursor - QApplication::restoreOverrideCursor(); - } -} - - void Polyhedron_demo_surface_reconstruction_plugin::advancing_front_reconstruction (const Polyhedron_demo_surface_reconstruction_plugin_dialog& dialog) { @@ -1293,17 +206,27 @@ void Polyhedron_demo_surface_reconstruction_plugin::advancing_front_reconstructi QApplication::setOverrideCursor(Qt::WaitCursor); std::cerr << "Advancing front reconstruction... "; - Scene_surface_mesh_item* reco_item = new Scene_surface_mesh_item(SMesh()); - SurfaceReconstruction::advancing_front (*points, reco_item, - dialog.longest_edge (), - dialog.radius_ratio_bound (), - CGAL_PI * dialog.beta_angle () / 180.); + + // Reconstruct point set as a polyhedron + SMesh* mesh = advancing_front (*points, + dialog.longest_edge(), + dialog.radius_ratio_bound(), + CGAL_PI * dialog.beta_angle() / 180., + dialog.structuring()); + if (mesh) + { + // Add polyhedron to scene + Scene_surface_mesh_item* new_item = new Scene_surface_mesh_item(mesh); + new_item->setName(tr("%1 (advancing front)").arg(pts_item->name())); + new_item->setColor(Qt::lightGray); + new_item->invalidateOpenGLBuffers(); + scene->addItem(new_item); + + // Hide point set + pts_item->setVisible(false); + scene->itemChanged(index); + } - reco_item->setName(tr("%1 (advancing front)").arg(scene->item(index)->name())); - reco_item->setColor(Qt::lightGray); - reco_item->setRenderingMode(FlatPlusEdges); - reco_item->invalidateOpenGLBuffers(); - scene->addItem(reco_item); QApplication::restoreOverrideCursor(); } } @@ -1329,16 +252,16 @@ void Polyhedron_demo_surface_reconstruction_plugin::scale_space_reconstruction std::vector reco_items; - SurfaceReconstruction::scale_space (*points, std::back_inserter (reco_items), - dialog.scalespace_js(), - dialog.iterations(), - dialog.neighbors(), dialog.fitting(), dialog.monge(), - dialog.neighborhood_size (), dialog.samples(), - dialog.scalespace_af(), - dialog.generate_smoothed (), - dialog.longest_edge_2(), dialog.radius_ratio_bound_2(), - CGAL_PI * dialog.beta_angle_2 () / 180., - dialog.separate_shells (), dialog.force_manifold ()); + scale_space (*points, reco_items, + dialog.scalespace_js(), + dialog.iterations(), + dialog.neighbors(), dialog.fitting(), dialog.monge(), + dialog.neighborhood_size (), dialog.samples(), + dialog.scalespace_af(), + dialog.generate_smoothed (), + dialog.longest_edge_2(), dialog.radius_ratio_bound_2(), + CGAL_PI * dialog.beta_angle_2 () / 180., + dialog.separate_shells (), dialog.force_manifold ()); for (std::size_t i = 0; i < reco_items.size (); ++ i) { @@ -1401,41 +324,24 @@ void Polyhedron_demo_surface_reconstruction_plugin::poisson_reconstruction const double sm_angle = dialog.angle (); const double sm_radius = dialog.radius (); const double sm_distance = dialog.distance (); - const QString sm_solver = dialog.solver(); + bool conjugate_gradient = dialog.conjugate_gradient(); bool use_two_passes = dialog.two_passes(); bool do_not_fill_holes = dialog.do_not_fill_holes(); QApplication::setOverrideCursor(Qt::WaitCursor); - if (!(point_set_item->has_normals())) - { - std::cerr << "Estimation of normal vectors... "; - points->add_normal_map(); - SurfaceReconstruction::compute_normals (*points, 12); - - - point_set_item->setRenderingMode(PointsPlusNormals); - - } - - // Reconstruct point set as a polyhedron - SMesh* smRemesh= NULL; - smRemesh = poisson_reconstruct_sm(*points, sm_angle, sm_radius, sm_distance, sm_solver, use_two_passes, - do_not_fill_holes); - if(smRemesh) + SMesh* mesh = poisson_reconstruct (*points, sm_angle, sm_radius, sm_distance, + conjugate_gradient, use_two_passes, + do_not_fill_holes); + if (mesh) { // Add polyhedron to scene - Scene_surface_mesh_item* new_item = new Scene_surface_mesh_item(smRemesh); - new_item->setName(tr("%1 Poisson (%2 %3 %4)") - .arg(point_set_item->name()) - .arg(sm_angle) - .arg(sm_radius) - .arg(sm_distance)); + Scene_surface_mesh_item* new_item = new Scene_surface_mesh_item(mesh); + new_item->setName(tr("%1 Poisson").arg(point_set_item->name())); new_item->setColor(Qt::lightGray); scene->addItem(new_item); - // Hide point set point_set_item->setVisible(false); scene->itemChanged(index); @@ -1445,15 +351,4 @@ void Polyhedron_demo_surface_reconstruction_plugin::poisson_reconstruction } } -void Polyhedron_demo_surface_reconstruction_plugin::ransac_reconstruction -(const Polyhedron_demo_surface_reconstruction_plugin_dialog& dialog) -{ - if (dialog.region_growing()) - region_growing_reconstruction_impl(dialog); - else - ransac_reconstruction_impl(dialog); -} - - - #include "Surface_reconstruction_plugin.moc" diff --git a/Polyhedron/demo/Polyhedron/Plugins/Point_set/Surface_reconstruction_plugin.ui b/Polyhedron/demo/Polyhedron/Plugins/Point_set/Surface_reconstruction_plugin.ui index 3cdd1a75902..b7b49c2dc31 100644 --- a/Polyhedron/demo/Polyhedron/Plugins/Point_set/Surface_reconstruction_plugin.ui +++ b/Polyhedron/demo/Polyhedron/Plugins/Point_set/Surface_reconstruction_plugin.ui @@ -6,8 +6,8 @@ 0 0 - 725 - 429 + 702 + 463 @@ -23,10 +23,20 @@ Advancing Front - - - QFormLayout::AllNonFixedFieldsGrow - + + + + + Qt::Vertical + + + + 20 + 40 + + + + @@ -41,14 +51,14 @@ - + Radius ratio bound - + 0.010000000000000 @@ -61,14 +71,14 @@ - + Beta angle - + ° @@ -81,18 +91,15 @@ - - - - Qt::Vertical + + + + false - - - 20 - 40 - + + Use Point Set Structuring (run Shape Detection first) - + @@ -113,117 +120,8 @@ Poisson - - - - - Min triangle angle: - - - - - - - ° - - - 1.000000000000000 - - - 30.000000000000000 - - - 20.000000000000000 - - - - - - - Max triangle size: - - - - - - - * average spacing - - - 0 - - - 1.000000000000000 - - - 1000.000000000000000 - - - 1.000000000000000 - - - 100.000000000000000 - - - - - - - Approximation error: - - - - - - - * average spacing - - - 6 - - - 0.010000000000000 - - - 100.000000000000000 - - - 0.010000000000000 - - - 0.250000000000000 - - - - - - - Solver: - - - - - - - - - - - - - - Perform two passes - - - - - - - Do not fill holes - - - - + + Qt::Vertical @@ -236,7 +134,198 @@ - + + + + + + + + Mesher: + + + + + + + + + Marching Tetrahedra + + + true + + + + + + + Surface Mesher + + + + + + + + + + + Conjugate Gradient + + + true + + + + + + + Simplicial LDLT + + + + + + + + + Solver: + + + + + + + + + Qt::Horizontal + + + + 40 + 20 + + + + + + + + false + + + + + + Surface Mesher Parameters + + + + + + + + Approximation error: + + + + + + + * average spacing + + + 6 + + + 0.010000000000000 + + + 100.000000000000000 + + + 0.010000000000000 + + + 0.250000000000000 + + + + + + + Min triangle angle: + + + + + + + ° + + + 1.000000000000000 + + + 30.000000000000000 + + + 20.000000000000000 + + + + + + + * average spacing + + + 0 + + + 1.000000000000000 + + + 1000.000000000000000 + + + 1.000000000000000 + + + 100.000000000000000 + + + + + + + Max triangle size: + + + + + + + + + + + + + + Perform two passes + + + + + + + Do not fill holes + + + + Qt::Vertical @@ -662,208 +751,6 @@ - - - Feature Preserving - - - - QFormLayout::AllNonFixedFieldsGrow - - - - - Qt::Vertical - - - - 20 - 40 - - - - - - - - Region growing - - - true - - - - - - - RANSAC - - - - - - - Connectivity tolerance - - - - - - - 3 - - - 0.001000000000000 - - - 0.010000000000000 - - - 0.010000000000000 - - - - - - - Noise tolerance - - - - - - - 3 - - - 0.001000000000000 - - - 0.001000000000000 - - - - - - - Normal tolerance - - - - - - - 0.010000000000000 - - - 1.000000000000000 - - - 0.010000000000000 - - - 0.900000000000000 - - - - - - - Minimum size of subset - - - - - - - 1 - - - 1000000000 - - - 10 - - - 100 - - - - - - - Generate structured point set - - - - - - - Qt::Vertical - - - - 20 - 40 - - - - - - - - - Automatic - - - - - - Qt::Vertical - - - - 20 - 40 - - - - - - - - Output surface has boundaries - - - true - - - - - - - Output surface must pass exactly through input points - - - true - - - - - - - Qt::Vertical - - - - 20 - 40 - - - - - - @@ -910,34 +797,18 @@ - m_scalespace_jet + m_scalespace_as toggled(bool) - frame_jet + frame_as setEnabled(bool) - 215 - 87 + 565 + 246 - 215 - 158 - - - - - m_scalespace_pca - toggled(bool) - frame_pca - setEnabled(bool) - - - 215 - 236 - - - 215 - 307 + 565 + 308 @@ -958,18 +829,50 @@ - m_scalespace_as + m_scalespace_pca toggled(bool) - frame_as + frame_pca setEnabled(bool) - 565 - 246 + 215 + 236 - 565 - 308 + 215 + 307 + + + + + m_scalespace_jet + toggled(bool) + frame_jet + setEnabled(bool) + + + 215 + 87 + + + 215 + 158 + + + + + m_surface_mesher + toggled(bool) + m_surface_mesh_parameters + setEnabled(bool) + + + 149 + 251 + + + 532 + 200 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 deleted file mode 100644 index 16e090f5b00..00000000000 --- a/Polyhedron/demo/Polyhedron/Plugins/Point_set/Surface_reconstruction_plugin_impl.cpp +++ /dev/null @@ -1,262 +0,0 @@ -//---------------------------------------------------------- -// Poisson reconstruction method: -// Reconstructs a surface mesh from a point set and returns it as a polyhedron. -//---------------------------------------------------------- - - -// CGAL -#include // must be included before kernel -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#ifdef CGAL_EIGEN3_ENABLED -#include -#endif - -#include - -#include "Kernel_type.h" -#include "SMesh_type.h" -#include "Scene_points_with_normal_item.h" - - - -// Concurrency -#ifdef CGAL_LINKED_WITH_TBB -typedef CGAL::Parallel_tag Concurrency_tag; -#else -typedef CGAL::Sequential_tag Concurrency_tag; -#endif - - -// Poisson reconstruction method: -// Reconstructs a surface mesh from a point set and returns it as a polyhedron. -template -bool poisson_reconstruct(FaceGraph* graph, - Point_set& points, - typename Traits::FT sm_angle, // Min triangle angle (degrees). - typename Traits::FT sm_radius, // Max triangle size w.r.t. point set average spacing. - typename Traits::FT sm_distance, // Approximation error w.r.t. point set average spacing. - const QString& solver_name, // solver name - bool use_two_passes, - bool do_not_fill_holes) -{ - // Poisson implicit function - typedef CGAL::Poisson_reconstruction_function Poisson_reconstruction_function; - - // Surface mesher - typedef CGAL::Surface_mesh_default_triangulation_3 STr; - typedef CGAL::Surface_mesh_complex_2_in_triangulation_3 C2t3; - typedef CGAL::Implicit_surface_3 Surface_3; - - // AABB tree - typedef CGAL::AABB_face_graph_triangle_primitive Primitive; - typedef CGAL::AABB_traits AABB_traits; - typedef CGAL::AABB_tree AABB_tree; - CGAL::Timer task_timer; task_timer.start(); - - //*************************************** - // Checks requirements - //*************************************** - - if (points.size() == 0) - { - std::cerr << "Error: empty point set" << std::endl; - return false; - } - - bool points_have_normals = points.has_normal_map(); - if ( ! points_have_normals ) - { - std::cerr << "Input point set not supported: this reconstruction method requires oriented normals" << std::endl; - return false; - } - - CGAL::Timer reconstruction_timer; reconstruction_timer.start(); - - //*************************************** - // Computes implicit function - //*************************************** - - - std::cerr << "Computes Poisson implicit function " - << "using " << solver_name.toLatin1().data() << " solver...\n"; - - - // Creates implicit function from the point set. - // Note: this method requires an iterator over points - // + property maps to access each point's position and normal. - Poisson_reconstruction_function function(points.begin_or_selection_begin(), points.end(), - points.point_map(), points.normal_map()); - - bool ok = false; - #ifdef CGAL_EIGEN3_ENABLED - if(solver_name=="Eigen - built-in simplicial LDLt") - { - CGAL::Eigen_solver_traits::EigenType> > solver; - ok = function.compute_implicit_function(solver, use_two_passes); - } - if(solver_name=="Eigen - built-in CG") - { - CGAL::Eigen_solver_traits::EigenType> > solver; - solver.solver().setTolerance(1e-6); - solver.solver().setMaxIterations(1000); - ok = function.compute_implicit_function(solver, use_two_passes); - } - #endif - - // Computes the Poisson indicator function f() - // at each vertex of the triangulation. - if ( ! ok ) - { - std::cerr << "Error: cannot compute implicit function" << std::endl; - return false; - } - - // Prints status - std::cerr << "Total implicit function (triangulation+refinement+solver): " << task_timer.time() << " seconds\n"; - task_timer.reset(); - - //*************************************** - // Surface mesh generation - //*************************************** - - std::cerr << "Surface meshing...\n"; - - // Computes average spacing - Kernel::FT average_spacing = CGAL::compute_average_spacing(points.all_or_selection_if_not_empty(), - 6 /* knn = 1 ring */, - points.parameters()); - - // Gets one point inside the implicit surface - Kernel::Point_3 inner_point = function.get_inner_point(); - Kernel::FT inner_point_value = function(inner_point); - if(inner_point_value >= 0.0) - { - std::cerr << "Error: unable to seed (" << inner_point_value << " at inner_point)" << std::endl; - return false; - } - - // Gets implicit function's radius - Kernel::Sphere_3 bsphere = function.bounding_sphere(); - Kernel::FT radius = std::sqrt(bsphere.squared_radius()); - - // Defines the implicit surface: requires defining a - // conservative bounding sphere centered at inner point. - Kernel::FT sm_sphere_radius = 5.0 * radius; - Kernel::FT sm_dichotomy_error = sm_distance*average_spacing/1000.0; // Dichotomy error must be << sm_distance - Surface_3 surface(function, - Kernel::Sphere_3(inner_point,sm_sphere_radius*sm_sphere_radius), - sm_dichotomy_error/sm_sphere_radius); - - // Defines surface mesh generation criteria - CGAL::Surface_mesh_default_criteria_3 criteria(sm_angle, // Min triangle angle (degrees) - sm_radius*average_spacing, // Max triangle size - sm_distance*average_spacing); // Approximation error - - CGAL_TRACE_STREAM << " make_surface_mesh(sphere center=("<::face_descriptor> faces_to_keep; - - for (Point_set::const_iterator p=points.begin_or_selection_begin(); p!=points.end(); p++) - { - typename AABB_traits::Point_and_primitive_id pap = tree.closest_point_and_primitive (points.point (*p)); - double distance = std::sqrt(CGAL::squared_distance (pap.first, points.point(*p))); - - max_distance = (std::max)(max_distance, distance); - avg_distance += distance; - - typename boost::graph_traits::face_descriptor f = pap.second; - faces_to_keep.insert (f); - } - avg_distance /= double(points.size()); - - std::cerr << "Reconstruction error:\n" - << " max = " << max_distance << " = " << max_distance/average_spacing << " * average spacing\n" - << " avg = " << avg_distance << " = " << avg_distance/average_spacing << " * average spacing\n"; - - if (do_not_fill_holes) - { - typename boost::graph_traits::face_iterator it = faces(*graph).begin (); - while (it != faces(*graph).end ()) - { - typename boost::graph_traits::face_iterator current = it ++; - - if (faces_to_keep.find (*current) == faces_to_keep.end ()) - { - CGAL::Euler::remove_face(halfedge (*current, *graph), *graph); - } - - } - - } - return true; -} - -SMesh* poisson_reconstruct_sm(Point_set& points, - Kernel::FT sm_angle, // Min triangle angle (degrees). - Kernel::FT sm_radius, // Max triangle size w.r.t. point set average spacing. - Kernel::FT sm_distance, // Approximation error w.r.t. point set average spacing. - const QString& solver_name, // solver name - bool use_two_passes, - bool do_not_fill_holes) -{ - SMesh* res = new SMesh; - if(poisson_reconstruct(res, points, sm_angle, sm_radius, sm_distance, - solver_name, use_two_passes, do_not_fill_holes)) - return res; - else - { - delete res; - return NULL; - } -} diff --git a/Polyhedron/demo/Polyhedron/Plugins/Point_set/Surface_reconstruction_poisson_impl.cpp b/Polyhedron/demo/Polyhedron/Plugins/Point_set/Surface_reconstruction_poisson_impl.cpp new file mode 100644 index 00000000000..3ad588c776f --- /dev/null +++ b/Polyhedron/demo/Polyhedron/Plugins/Point_set/Surface_reconstruction_poisson_impl.cpp @@ -0,0 +1,238 @@ +//---------------------------------------------------------- +// Poisson reconstruction method: +// Reconstructs a surface mesh from a point set and returns it as a polyhedron. +//---------------------------------------------------------- + +// CGAL +#include // must be included before kernel +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#include + +#include "Kernel_type.h" +#include "SMesh_type.h" +#include "Scene_points_with_normal_item.h" + + + +// Concurrency +#ifdef CGAL_LINKED_WITH_TBB +typedef CGAL::Parallel_tag Concurrency_tag; +#else +typedef CGAL::Sequential_tag Concurrency_tag; +#endif + + +// Poisson reconstruction method: +// Reconstructs a surface mesh from a point set and returns it as a polyhedron. +SMesh* poisson_reconstruct_sm(Point_set& points, + Kernel::FT sm_angle, // Min triangle angle (degrees). + Kernel::FT sm_radius, // Max triangle size w.r.t. point set average spacing. + Kernel::FT sm_distance, // Approximation error w.r.t. point set average spacing. + bool conjugate_gradient, + bool use_two_passes, + bool do_not_fill_holes) +{ + // Poisson implicit function + typedef CGAL::Poisson_reconstruction_function Poisson_reconstruction_function; + + // Surface mesher + typedef CGAL::Surface_mesh_default_triangulation_3 STr; + typedef CGAL::Surface_mesh_complex_2_in_triangulation_3 C2t3; + typedef CGAL::Implicit_surface_3 Surface_3; + + // AABB tree + typedef CGAL::AABB_face_graph_triangle_primitive Primitive; + typedef CGAL::AABB_traits AABB_traits; + typedef CGAL::AABB_tree AABB_tree; + CGAL::Timer task_timer; task_timer.start(); + + //*************************************** + // Checks requirements + //*************************************** + + if (points.size() == 0) + { + std::cerr << "Error: empty point set" << std::endl; + return nullptr; + } + + bool points_have_normals = points.has_normal_map(); + if ( ! points_have_normals ) + { + std::cerr << "Input point set not supported: this reconstruction method requires oriented normals" << std::endl; + return nullptr; + } + + CGAL::Timer reconstruction_timer; reconstruction_timer.start(); + + //*************************************** + // Computes implicit function + //*************************************** + + + std::cerr << "Computes Poisson implicit function " + << "using " << (conjugate_gradient ? "Conjugate Gradient" : "Simplicial LDLT") << "..." << std::endl; + + + // Creates implicit function from the point set. + // Note: this method requires an iterator over points + // + property maps to access each point's position and normal. + Poisson_reconstruction_function function(points.begin_or_selection_begin(), points.end(), + points.point_map(), points.normal_map()); + + bool ok = false; + if(conjugate_gradient) + { + CGAL::Eigen_solver_traits::EigenType> > solver; + ok = function.compute_implicit_function(solver, use_two_passes); + } + else + { + CGAL::Eigen_solver_traits::EigenType> > solver; + solver.solver().setTolerance(1e-6); + solver.solver().setMaxIterations(1000); + ok = function.compute_implicit_function(solver, use_two_passes); + } + + // Computes the Poisson indicator function f() + // at each vertex of the triangulation. + if ( ! ok ) + { + std::cerr << "Error: cannot compute implicit function" << std::endl; + return nullptr; + } + + // Prints status + std::cerr << "Total implicit function (triangulation+refinement+solver): " << task_timer.time() << " seconds\n"; + task_timer.reset(); + + //*************************************** + // Surface mesh generation + //*************************************** + + std::cerr << "Surface meshing...\n"; + + // Computes average spacing + Kernel::FT average_spacing = CGAL::compute_average_spacing(points.all_or_selection_if_not_empty(), + 6 /* knn = 1 ring */, + points.parameters()); + + // Gets one point inside the implicit surface + Kernel::Point_3 inner_point = function.get_inner_point(); + Kernel::FT inner_point_value = function(inner_point); + if(inner_point_value >= 0.0) + { + std::cerr << "Error: unable to seed (" << inner_point_value << " at inner_point)" << std::endl; + return nullptr; + } + + // Gets implicit function's radius + Kernel::Sphere_3 bsphere = function.bounding_sphere(); + Kernel::FT radius = std::sqrt(bsphere.squared_radius()); + + // Defines the implicit surface: requires defining a + // conservative bounding sphere centered at inner point. + Kernel::FT sm_sphere_radius = 5.0 * radius; + Kernel::FT sm_dichotomy_error = sm_distance*average_spacing/1000.0; // Dichotomy error must be << sm_distance + Surface_3 surface(function, + Kernel::Sphere_3(inner_point,sm_sphere_radius*sm_sphere_radius), + sm_dichotomy_error/sm_sphere_radius); + + // Defines surface mesh generation criteria + CGAL::Surface_mesh_default_criteria_3 criteria(sm_angle, // Min triangle angle (degrees) + sm_radius*average_spacing, // Max triangle size + sm_distance*average_spacing); // Approximation error + + CGAL_TRACE_STREAM << " make_surface_mesh(sphere center=("<::face_descriptor> faces_to_keep; + + for (Point_set::const_iterator p=points.begin_or_selection_begin(); p!=points.end(); p++) + { + typename AABB_traits::Point_and_primitive_id pap = tree.closest_point_and_primitive (points.point (*p)); + double distance = std::sqrt(CGAL::squared_distance (pap.first, points.point(*p))); + + max_distance = (std::max)(max_distance, distance); + avg_distance += distance; + + typename boost::graph_traits::face_descriptor f = pap.second; + faces_to_keep.insert (f); + } + avg_distance /= double(points.size()); + + std::cerr << "Reconstruction error:\n" + << " max = " << max_distance << " = " << max_distance/average_spacing << " * average spacing\n" + << " avg = " << avg_distance << " = " << avg_distance/average_spacing << " * average spacing\n"; + + if (do_not_fill_holes) + { + typename boost::graph_traits::face_iterator it = faces(*mesh).begin (); + while (it != faces(*mesh).end ()) + { + typename boost::graph_traits::face_iterator current = it ++; + + if (faces_to_keep.find (*current) == faces_to_keep.end ()) + { + CGAL::Euler::remove_face(halfedge (*current, *mesh), *mesh); + } + + } + + } + return mesh; +} + diff --git a/Polyhedron/demo/Polyhedron/Plugins/Point_set/Surface_reconstruction_scale_space_impl.cpp b/Polyhedron/demo/Polyhedron/Plugins/Point_set/Surface_reconstruction_scale_space_impl.cpp new file mode 100644 index 00000000000..582bad37820 --- /dev/null +++ b/Polyhedron/demo/Polyhedron/Plugins/Point_set/Surface_reconstruction_scale_space_impl.cpp @@ -0,0 +1,228 @@ +#include "Kernel_type.h" +#include "SMesh_type.h" +#include "Scene_points_with_normal_item.h" +#include "Scene_polygon_soup_item.h" +#include "Point_set_3.h" + +#include +#include +#include +#include +#include +#include + + +// Concurrency +#ifdef CGAL_LINKED_WITH_TBB +typedef CGAL::Parallel_tag Concurrency_tag; +#else +typedef CGAL::Sequential_tag Concurrency_tag; +#endif + +typedef CGAL::Scale_space_surface_reconstruction_3 ScaleSpace; +typedef CGAL::Scale_space_reconstruction_3::Advancing_front_mesher ScaleSpaceAFM; +typedef CGAL::Scale_space_reconstruction_3::Alpha_shape_mesher ScaleSpaceASM; +typedef CGAL::Scale_space_reconstruction_3::Jet_smoother ScaleSpaceJS; +typedef CGAL::Scale_space_reconstruction_3::Weighted_PCA_smoother ScaleSpaceWPS; + +void scale_space (const Point_set& points, + std::vector& items, + bool jet_smoother, + unsigned int iterations, + unsigned int neighbors, unsigned int fitting, unsigned int monge, + unsigned int neighborhood_size, unsigned int samples, + bool advancing_front_mesher, + bool generate_smooth, + double longest_edge, double radius_ratio_bound, double beta_angle, + bool separate_shells, bool force_manifold) +{ + ScaleSpace reconstruct (points.points().begin(), points.points().end()); + + double squared_radius = 0.; + if (jet_smoother) + { + ScaleSpaceJS smoother(neighbors, fitting, monge); + reconstruct.increase_scale(iterations, smoother); + if (!advancing_front_mesher) + squared_radius = CGAL::compute_average_spacing (points, neighbors); + } + else + { + ScaleSpaceWPS smoother(neighborhood_size, samples); + reconstruct.increase_scale(iterations, smoother); + squared_radius = smoother.squared_radius(); + + } + + if (advancing_front_mesher) + { + ScaleSpaceAFM mesher (longest_edge, radius_ratio_bound, beta_angle); + reconstruct.reconstruct_surface (mesher); + + Scene_polygon_soup_item* new_item + = new Scene_polygon_soup_item (); + new_item->setColor(Qt::lightGray); + new_item->setRenderingMode(FlatPlusEdges); + new_item->init_polygon_soup(points.size(), reconstruct.number_of_facets ()); + + Scene_polygon_soup_item* smooth_item = NULL; + if (generate_smooth) + { + smooth_item = new Scene_polygon_soup_item (); + smooth_item->setColor(Qt::lightGray); + smooth_item->setRenderingMode(FlatPlusEdges); + smooth_item->init_polygon_soup(points.size(), reconstruct.number_of_facets ()); + } + + std::map map_i2i; + std::size_t current_index = 0; + + for (ScaleSpace::Facet_iterator it = reconstruct.facets_begin(); + it != reconstruct.facets_end(); ++ it) + { + for (unsigned int ind = 0; ind < 3; ++ ind) + { + if (map_i2i.find ((*it)[ind]) == map_i2i.end ()) + { + map_i2i.insert (std::make_pair ((*it)[ind], current_index ++)); + Point_3 p = points.point(*(points.begin_or_selection_begin() + (*it)[ind])); + new_item->new_vertex (p.x (), p.y (), p.z ()); + + if (generate_smooth) + { + p = *(reconstruct.points_begin() + (*it)[ind]); + smooth_item->new_vertex (p.x (), p.y (), p.z ()); + } + } + } + new_item->new_triangle( map_i2i[(*it)[0]], + map_i2i[(*it)[1]], + map_i2i[(*it)[2]] ); + if (generate_smooth) + smooth_item->new_triangle( map_i2i[(*it)[0]], + map_i2i[(*it)[1]], + map_i2i[(*it)[2]] ); + + } + + items.push_back(new_item); + if (generate_smooth) + items.push_back (smooth_item); + } + else + { + ScaleSpaceASM mesher (squared_radius, separate_shells, force_manifold); + reconstruct.reconstruct_surface (mesher); + + for( unsigned int sh = 0; sh < mesher.number_of_shells(); ++sh ) + { + Scene_polygon_soup_item* new_item + = new Scene_polygon_soup_item (); + new_item->setColor(Qt::lightGray); + new_item->setRenderingMode(FlatPlusEdges); + new_item->init_polygon_soup(points.size(), mesher.number_of_triangles ()); + + Scene_polygon_soup_item* smooth_item = NULL; + if (generate_smooth) + { + smooth_item = new Scene_polygon_soup_item (); + smooth_item->setColor(Qt::lightGray); + smooth_item->setRenderingMode(FlatPlusEdges); + smooth_item->init_polygon_soup(points.size(), mesher.number_of_triangles ()); + } + + std::map map_i2i; + unsigned int current_index = 0; + + for (ScaleSpaceASM::Facet_iterator it = mesher.shell_begin (sh); + it != mesher.shell_end (sh); ++ it) + { + for (unsigned int ind = 0; ind < 3; ++ ind) + { + if (map_i2i.find ((*it)[ind]) == map_i2i.end ()) + { + map_i2i.insert (std::make_pair ((*it)[ind], current_index ++)); + Point_3 p = points.point(*(points.begin_or_selection_begin() + (*it)[ind])); + new_item->new_vertex (p.x (), p.y (), p.z ()); + + if (generate_smooth) + { + p = *(reconstruct.points_begin() + (*it)[ind]); + smooth_item->new_vertex (p.x (), p.y (), p.z ()); + } + } + } + new_item->new_triangle( map_i2i[(*it)[0]], + map_i2i[(*it)[1]], + map_i2i[(*it)[2]] ); + if (generate_smooth) + smooth_item->new_triangle( map_i2i[(*it)[0]], + map_i2i[(*it)[1]], + map_i2i[(*it)[2]] ); + + } + + items.push_back (new_item); + if (generate_smooth) + items.push_back (smooth_item); + } + + if (force_manifold) + { + std::ptrdiff_t num = std::distance( mesher.garbage_begin( ), + mesher.garbage_end( ) ); + + Scene_polygon_soup_item* new_item + = new Scene_polygon_soup_item (); + new_item->setColor(Qt::blue); + new_item->setRenderingMode(FlatPlusEdges); + new_item->init_polygon_soup(points.size(), num); + + Scene_polygon_soup_item* smooth_item = NULL; + if (generate_smooth) + { + smooth_item = new Scene_polygon_soup_item (); + smooth_item->setColor(Qt::blue); + smooth_item->setRenderingMode(FlatPlusEdges); + smooth_item->init_polygon_soup(points.size(), num); + } + + std::map map_i2i; + + std::size_t current_index = 0; + for (ScaleSpaceASM::Facet_iterator it=mesher.garbage_begin(), + end=mesher.garbage_end();it!=end;++it) + { + for (unsigned int ind = 0; ind < 3; ++ ind) + { + if (map_i2i.find ((*it)[ind]) == map_i2i.end ()) + { + map_i2i.insert (std::make_pair ((*it)[ind], current_index ++)); + Point_3 p = points.point(*(points.begin_or_selection_begin() + (*it)[ind])); + new_item->new_vertex (p.x (), p.y (), p.z ()); + + if (generate_smooth) + { + p = *(reconstruct.points_begin() + (*it)[ind]); + smooth_item->new_vertex (p.x (), p.y (), p.z ()); + } + } + + } + new_item->new_triangle( map_i2i[(*it)[0]], + map_i2i[(*it)[1]], + map_i2i[(*it)[2]] ); + if (generate_smooth) + smooth_item->new_triangle( map_i2i[(*it)[0]], + map_i2i[(*it)[1]], + map_i2i[(*it)[2]] ); + } + + items.push_back (new_item); + if (generate_smooth) + items.push_back (smooth_item); + + } + } +} +