diff --git a/Surface_reconstruction_points_3/include/CGAL/Poisson_mesh_cell_criteria_3.h b/Surface_reconstruction_points_3/include/CGAL/Poisson_mesh_cell_criteria_3.h index f2e24e31948..3df0fb7b257 100644 --- a/Surface_reconstruction_points_3/include/CGAL/Poisson_mesh_cell_criteria_3.h +++ b/Surface_reconstruction_points_3/include/CGAL/Poisson_mesh_cell_criteria_3.h @@ -24,11 +24,39 @@ namespace CGAL { +namespace internal { +namespace Poisson { + template +class Constant_sizing_field { + double sq_radius_bound; +public: + double cell_radius_bound() const { return CGAL::sqrt(sq_radius_bound); } + + Constant_sizing_field(double sq_radius_bound = 0.) + : sq_radius_bound(sq_radius_bound) {} + + template + double operator()(const Point&) const { return sq_radius_bound; } +}; // end class Constant_sizing_field + +} // end namespace Poisson +} // end namespace internal + +template < + class Tr, + typename Sizing_field = internal::Poisson::Constant_sizing_field, + typename Second_sizing_field = internal::Poisson::Constant_sizing_field + > class Poisson_mesh_cell_criteria_3 { - double squared_radius_bound_; + Sizing_field sizing_field; + Second_sizing_field second_sizing_field; double radius_edge_bound_; + + typedef Poisson_mesh_cell_criteria_3 Self; public: struct Cell_quality : public std::pair { @@ -56,25 +84,34 @@ class Poisson_mesh_cell_criteria_3 } }; - inline - double squared_radius_bound() const - { - return squared_radius_bound_; - } + // inline + // double squared_radius_bound() const + // { + // return squared_radius_bound_; + // } typedef typename Tr::Cell_handle Cell_handle; Poisson_mesh_cell_criteria_3(const double radius_edge_bound = 2, //< radius edge ratio bound (ignored if zero) const double radius_bound = 0) //< cell radius bound (ignored if zero) - : squared_radius_bound_(radius_bound*radius_bound), + : sizing_field(radius_bound*radius_bound), + second_sizing_field(), radius_edge_bound_(radius_edge_bound) {} - inline - void set_squared_radius_bound(const double squared_radius_bound) - { - squared_radius_bound_ = squared_radius_bound; - } + Poisson_mesh_cell_criteria_3(const double radius_edge_bound, //< radius edge ratio bound (ignored if zero) + const Sizing_field& sizing_field, + const Second_sizing_field second_sizing_field = Second_sizing_field()) + : sizing_field(sizing_field), + second_sizing_field(second_sizing_field), + radius_edge_bound_(radius_edge_bound) + {} + + // inline + // void set_squared_radius_bound(const double squared_radius_bound) + // { + // squared_radius_bound_ = squared_radius_bound; + // } inline double radius_edge_bound() const @@ -91,15 +128,12 @@ class Poisson_mesh_cell_criteria_3 class Is_bad { protected: - const double radius_edge_bound_; - const double squared_radius_bound_; + const Self* cell_criteria; public: typedef typename Tr::Point Point_3; - Is_bad(const double radius_edge_bound, - const double squared_radius_bound) - : radius_edge_bound_(radius_edge_bound), - squared_radius_bound_(squared_radius_bound) {} + Is_bad(const Self* cell_criteria) + : cell_criteria(cell_criteria) {} bool operator()(const Cell_handle& c, Cell_quality& qual) const @@ -114,23 +148,30 @@ class Poisson_mesh_cell_criteria_3 typedef typename Geom_traits::Compute_squared_distance_3 Distance; typedef typename Geom_traits::FT FT; - Radius radius = Geom_traits().compute_squared_radius_3_object(); + Radius sq_radius = Geom_traits().compute_squared_radius_3_object(); Distance distance = Geom_traits().compute_squared_distance_3_object(); - double size = to_double(radius(p, q, r, s)); - - if( squared_radius_bound_ != 0 ) - { - qual.second = size / squared_radius_bound_; - // normalized by size bound to deal - // with size field + double size = to_double(sq_radius(p, q, r, s)); + + // If the tetrahedron is small enough according to the second sizing + // field, then let's say the tetrahedron is OK according to the + // sizing fields. + if( size > cell_criteria->second_sizing_field(p) ) + { + // first sizing field + const double size_bound = cell_criteria->sizing_field(p); + if(size_bound > 0.) { + qual.second = size / size_bound; + // normalized by size bound to deal + // with size field if( qual.sq_size() > 1 ) - { - qual.first = 1; // (do not compute aspect) - return true; - } + { + qual.first = 1; // (do not compute aspect) + return true; + } } - if( radius_edge_bound_ == 0 ) + } + if( cell_criteria->radius_edge_bound_ == 0 ) { qual = Cell_quality(0,1); return false; @@ -145,14 +186,14 @@ class Poisson_mesh_cell_criteria_3 qual.first = size / min_sq_length; - return (qual.first > radius_edge_bound_); + return (qual.first > cell_criteria->radius_edge_bound_); } }; // end Is_bad Is_bad is_bad_object() const - { return Is_bad(radius_edge_bound_, squared_radius_bound_); } + { return Is_bad(this); } }; // end Poisson_mesh_cell_criteria_3 diff --git a/Surface_reconstruction_points_3/include/CGAL/Poisson_reconstruction_function.h b/Surface_reconstruction_points_3/include/CGAL/Poisson_reconstruction_function.h index a1bcb7c5541..0d8ec5a2819 100644 --- a/Surface_reconstruction_points_3/include/CGAL/Poisson_reconstruction_function.h +++ b/Surface_reconstruction_points_3/include/CGAL/Poisson_reconstruction_function.h @@ -19,6 +19,12 @@ #ifndef CGAL_POISSON_RECONSTRUCTION_FUNCTION_H #define CGAL_POISSON_RECONSTRUCTION_FUNCTION_H +#ifndef CGAL_DIV_NORMALIZED +# ifndef CGAL_DIV_NON_NORMALIZED +# define CGAL_DIV_NON_NORMALIZED 1 +# endif +#endif + #include #include #include @@ -27,18 +33,60 @@ #include #include #include -#include #include #include #include -#include #include #include +#include #include +#include namespace CGAL { + namespace internal { +template +bool +invert( + const RT& a0, const RT& a1, const RT& a2, + const RT& a3, const RT& a4, const RT& a5, + const RT& a6, const RT& a7, const RT& a8, + RT& i0, RT& i1, RT& i2, + RT& i3, RT& i4, RT& i5, + RT& i6, RT& i7, RT& i8) +{ + // Compute the adjoint. + i0 = a4*a8 - a5*a7; + i1 = a2*a7 - a1*a8; + i2 = a1*a5 - a2*a4; + i3 = a5*a6 - a3*a8; + i4 = a0*a8 - a2*a6; + i5 = a2*a3 - a0*a5; + i6 = a3*a7 - a4*a6; + i7 = a1*a6 - a0*a7; + i8 = a0*a4 - a1*a3; + + RT det = a0*i0 + a1*i3 + a2*i6; + + if(det != 0) { + RT idet = (RT(1.0))/det; + i0 *= idet; + i1 *= idet; + i2 *= idet; + i3 *= idet; + i4 *= idet; + i5 *= idet; + i6 *= idet; + i7 *= idet; + i8 *= idet; + return true; + } + + return false; +} + + } /// Given a set of 3D points with oriented normals sampled on the boundary of a 3D solid, /// the Poisson Surface Reconstruction method [Kazhdan06] solves for an approximate indicator function @@ -55,13 +103,51 @@ namespace CGAL { /// @heading Parameters: /// @param Gt Geometric traits class. -template + +struct Poisson_visitor { + void before_insertion() const + {} +}; + +struct Poisson_skip_vertices { + double ratio; + Poisson_skip_vertices(const double ratio) : ratio(ratio) {} + + template + bool operator()(Iterator) const { + return CGAL::default_random.get_double() < ratio; + } +}; + +template +struct Special_wrapper_of_two_functions_keep_pointers { + F1 *f1; + F2 *f2; + Special_wrapper_of_two_functions_keep_pointers(F1* f1, F2* f2) + : f1(f1), f2(f2) {} + + template + double operator()(const X& x) const { + return (std::max)((*f1)(x), CGAL::square((*f2)(x))); + } + + template + double operator()(const X& x) { + return (std::max)((*f1)(x), CGAL::square((*f2)(x))); + } +}; // end struct Special_wrapper_of_two_functions_keep_pointers + + template class Poisson_reconstruction_function { + // Public types public: typedef Gt Geom_traits; ///< Geometric traits class + typedef Reconstruction_triangulation_3 > + Triangulation; + typedef typename Triangulation::Cell_handle Cell_handle; // Geometric types typedef typename Geom_traits::FT FT; ///< typedef to Geom_traits::FT @@ -74,8 +160,6 @@ class Poisson_reconstruction_function /// Internal 3D triangulation, of type Reconstruction_triangulation_3. // Note: poisson_refine_triangulation() requires a robust circumcenter computation. - typedef Reconstruction_triangulation_3 > - Triangulation; // Repeat Triangulation types typedef typename Triangulation::Triangulation_data_structure Triangulation_data_structure; @@ -84,7 +168,6 @@ class Poisson_reconstruction_function typedef typename Geom_traits::Segment_3 Segment; typedef typename Geom_traits::Triangle_3 Triangle; typedef typename Geom_traits::Tetrahedron_3 Tetrahedron; - typedef typename Triangulation::Cell_handle Cell_handle; typedef typename Triangulation::Vertex_handle Vertex_handle; typedef typename Triangulation::Cell Cell; typedef typename Triangulation::Vertex Vertex; @@ -113,10 +196,17 @@ class Poisson_reconstruction_function // the Poisson equation Laplacian(f) = divergent(normals field). boost::shared_ptr m_tr; + mutable boost::shared_ptr > > m_Bary; + mutable std::vector Dual; + mutable std::vector Normal; + // contouring and meshing Point m_sink; // Point with the minimum value of operator() mutable Cell_handle m_hint; // last cell found = hint for next search + FT average_spacing; + + // Public methods public: @@ -131,14 +221,17 @@ class Poisson_reconstruction_function // This variant requires all parameters. template Poisson_reconstruction_function( InputIterator first, ///< iterator over the first input point. InputIterator beyond, ///< past-the-end iterator over the input points. PointPMap point_pmap, ///< property map to access the position of an input point. - NormalPMap normal_pmap) ///< property map to access the *oriented* normal of an input point. - : m_tr(new Triangulation) + NormalPMap normal_pmap, ///< property map to access the *oriented* normal of an input point. + Visitor visitor) + : m_tr(new Triangulation), m_Bary(new std::vector > ) + , average_spacing(CGAL::compute_average_spacing(first, beyond, 6)) { CGAL::Timer task_timer; task_timer.start(); CGAL_TRACE_STREAM << "Creates Poisson triangulation...\n"; @@ -147,24 +240,28 @@ class Poisson_reconstruction_function m_tr->insert( first,beyond, point_pmap, - normal_pmap); + normal_pmap, + Triangulation::INPUT, + visitor); // Prints status CGAL_TRACE_STREAM << "Creates Poisson triangulation: " << task_timer.time() << " seconds, " - << (CGAL::Memory_sizer().virtual_size()>>20) << " Mb allocated" << std::endl; } /// @cond SKIP_IN_MANUAL // This variant creates a default point property map = Dereference_property_map. template Poisson_reconstruction_function( InputIterator first, ///< iterator over the first input point. InputIterator beyond, ///< past-the-end iterator over the input points. - NormalPMap normal_pmap) ///< property map to access the *oriented* normal of an input point. - : m_tr(new Triangulation) + NormalPMap normal_pmap, ///< property map to access the *oriented* normal of an input point. + Visitor visitor) + : m_tr(new Triangulation), m_Bary(new std::vector > ) + , average_spacing(CGAL::compute_average_spacing(first, beyond, 6)) { CGAL::Timer task_timer; task_timer.start(); CGAL_TRACE_STREAM << "Creates Poisson triangulation...\n"; @@ -172,19 +269,27 @@ class Poisson_reconstruction_function // Inserts points in triangulation m_tr->insert( first,beyond, - normal_pmap); + normal_pmap, + Triangulation::INPUT, + visitor); // Prints status CGAL_TRACE_STREAM << "Creates Poisson triangulation: " << task_timer.time() << " seconds, " - << (CGAL::Memory_sizer().virtual_size()>>20) << " Mb allocated" << std::endl; } /// @endcond + + ~Poisson_reconstruction_function() + {} /// Returns a sphere bounding the inferred surface. Sphere bounding_sphere() const { - return m_tr->input_points_bounding_sphere(); + return m_tr->bounding_sphere(); + } + + const Triangulation& tr() const { + return *m_tr; } /// The function compute_implicit_function() must be called @@ -201,9 +306,13 @@ class Poisson_reconstruction_function /// @return false if the linear solver fails. // This variant requires all parameters. - template + template bool compute_implicit_function( - SparseLinearAlgebraTraits_d solver = SparseLinearAlgebraTraits_d()) ///< sparse linear solver + SparseLinearAlgebraTraits_d solver,// = SparseLinearAlgebraTraits_d(), + Visitor visitor, + double approximation_ratio = 0, + double average_spacing_ratio = 5) { CGAL::Timer task_timer; task_timer.start(); CGAL_TRACE_STREAM << "Delaunay refinement...\n"; @@ -214,16 +323,80 @@ class Poisson_reconstruction_function const FT enlarge_ratio = 1.5; const FT radius = sqrt(bounding_sphere().squared_radius()); // get triangulation's radius const FT cell_radius_bound = radius/5.; // large - unsigned int nb_vertices_added = delaunay_refinement(radius_edge_ratio_bound,cell_radius_bound,max_vertices,enlarge_ratio); + internal::Poisson::Constant_sizing_field sizing_field(CGAL::square(cell_radius_bound)); + + std::vector NB; + + NB.push_back( delaunay_refinement(radius_edge_ratio_bound,sizing_field,max_vertices,enlarge_ratio)); + + while(m_tr->insert_fraction(visitor)){ + + NB.push_back( delaunay_refinement(radius_edge_ratio_bound,sizing_field,max_vertices,enlarge_ratio)); + } + + if(approximation_ratio > 0. && + approximation_ratio * std::distance(m_tr->input_points_begin(), + m_tr->input_points_end()) > 20) { + + typedef Filter_iterator Some_points_iterator; + Poisson_skip_vertices skip(1.-approximation_ratio); + + CGAL_TRACE_STREAM << "SPECIAL PASS that uses an approximation of the result (approximation ratio: " + << approximation_ratio << ")" << std::endl; + CGAL::Timer approximation_timer; approximation_timer.start(); + + CGAL::Timer sizing_field_timer; sizing_field_timer.start(); + Poisson_reconstruction_function + coarse_poisson_function(Some_points_iterator(m_tr->input_points_end(), + skip, + m_tr->input_points_begin()), + Some_points_iterator(m_tr->input_points_end(), + skip), + Normal_of_point_with_normal_pmap(), + Poisson_visitor()); + coarse_poisson_function.compute_implicit_function(solver, Poisson_visitor(), + 0.); + internal::Poisson::Constant_sizing_field + min_sizing_field(CGAL::square(average_spacing)); + internal::Poisson::Constant_sizing_field + sizing_field_ok(CGAL::square(average_spacing*average_spacing_ratio)); + + Special_wrapper_of_two_functions_keep_pointers< + internal::Poisson::Constant_sizing_field, + Poisson_reconstruction_function > sizing_field2(&min_sizing_field, + &coarse_poisson_function); + + sizing_field_timer.stop(); + std::cerr << "Construction time of the sizing field: " << sizing_field_timer.time() + << " seconds" << std::endl; + + NB.push_back( delaunay_refinement(radius_edge_ratio_bound, + sizing_field2, + max_vertices, + enlarge_ratio, + sizing_field_ok) ); + approximation_timer.stop(); + CGAL_TRACE_STREAM << "SPECIAL PASS END (" << approximation_timer.time() << " seconds)" << std::endl; + } + + // Prints status - CGAL_TRACE_STREAM << "Delaunay refinement: " << "added " << nb_vertices_added << " Steiner points, " - << task_timer.time() << " seconds, " - << (CGAL::Memory_sizer().virtual_size()>>20) << " Mb allocated" - << std::endl; + CGAL_TRACE_STREAM << "Delaunay refinement: " << "added "; + for(int i = 0; i < NB.size()-1; i++){ + CGAL_TRACE_STREAM << NB[i] << " + "; + } + CGAL_TRACE_STREAM << NB.back() << " Steiner points, " + << task_timer.time() << " seconds, " + << std::endl; task_timer.reset(); +#ifdef CGAL_DIV_NON_NORMALIZED + CGAL_TRACE_STREAM << "Solve Poisson equation with non-normalized divergence...\n"; +#else CGAL_TRACE_STREAM << "Solve Poisson equation with normalized divergence...\n"; +#endif // Computes the Poisson indicator function operator() // at each vertex of the triangulation. @@ -241,25 +414,37 @@ class Poisson_reconstruction_function // Prints status CGAL_TRACE_STREAM << "Solve Poisson equation: " << task_timer.time() << " seconds, " - << (CGAL::Memory_sizer().virtual_size()>>20) << " Mb allocated" << std::endl; task_timer.reset(); return true; } - /// @cond SKIP_IN_MANUAL - // This variant provides the default sparse linear traits class = Taucs_symmetric_solver_traits. - bool compute_implicit_function() + + boost::tuple special_func(const Point& p) const { - return compute_implicit_function< Taucs_symmetric_solver_traits >(); + m_hint = m_tr->locate(p ,m_hint ); // no hint when we use hierarchy + + if(m_tr->is_infinite(m_hint)) { + int i = m_hint->index(m_tr->infinite_vertex()); + return boost::make_tuple(m_hint->vertex((i+1)&3)->f(), + m_hint, true); + } + + FT a,b,c,d; + barycentric_coordinates(p,m_hint,a,b,c,d); + return boost::make_tuple(a * m_hint->vertex(0)->f() + + b * m_hint->vertex(1)->f() + + c * m_hint->vertex(2)->f() + + d * m_hint->vertex(3)->f(), + m_hint, false); } - /// @endcond /// 'ImplicitFunction' interface: evaluates the implicit function at a given 3D query point. FT operator()(const Point& p) const { - m_hint = m_tr->locate(p,m_hint); + static int I=0; + m_hint = m_tr->locate(p ,m_hint); if(m_tr->is_infinite(m_hint)) { int i = m_hint->index(m_tr->infinite_vertex()); @@ -274,6 +459,84 @@ class Poisson_reconstruction_function d * m_hint->vertex(3)->f(); } + void initialize_cell_indices() + { + int i=0; + for(Finite_cells_iterator fcit = m_tr->finite_cells_begin(); + fcit != m_tr->finite_cells_end(); + ++fcit){ + fcit->info()= i++; + } + } + + void initialize_barycenters() const + { + m_Bary->resize(m_tr->number_of_cells()); + + for(int i=0; i< m_Bary->size();i++){ + (*m_Bary)[i][0]=-1; + } + } + + + + + void initialize_cell_normals() const + { + Normal.resize(m_tr->number_of_cells()); + int i = 0; + int N = 0; + for(Finite_cells_iterator fcit = m_tr->finite_cells_begin(); + fcit != m_tr->finite_cells_end(); + ++fcit){ + Normal[i] = cell_normal(fcit); + if(Normal[i] == NULL_VECTOR){ + N++; + } + ++i; + } + std::cerr << N << " out of " << i << " cells have NULL_VECTOR as normal" << std::endl; + } + + void initialize_duals() const + { + Dual.resize(m_tr->number_of_cells()); + int i = 0; + for(Finite_cells_iterator fcit = m_tr->finite_cells_begin(); + fcit != m_tr->finite_cells_end(); + ++fcit){ + Dual[i++] = m_tr->dual(fcit); + } + } + + void clear_duals() const + { + Dual.clear(); + } + + void clear_normals() const + { + Normal.clear(); + } + + void initialize_matrix_entry(Cell_handle ch) const + { + boost::array & entry = (*m_Bary)[ch->info()]; + const Point& pa = ch->vertex(0)->point(); + const Point& pb = ch->vertex(1)->point(); + const Point& pc = ch->vertex(2)->point(); + const Point& pd = ch->vertex(3)->point(); + + Vector va = pa - pd; + Vector vb = pb - pd; + Vector vc = pc - pd; + + internal::invert(va.x(), va.y(), va.z(), + vb.x(), vb.y(), vb.z(), + vc.x(), vc.y(), vc.z(), + entry[0],entry[1],entry[2],entry[3],entry[4],entry[5],entry[6],entry[7],entry[8]); + } + /// Returns a point located inside the inferred surface. Point get_inner_point() const { @@ -286,20 +549,32 @@ class Poisson_reconstruction_function /// Delaunay refinement (break bad tetrahedra, where /// bad means badly shaped or too big). The normal of - /// Steiner points is set to zero. - /// Returns the number of vertices inserted. + /// Steiner points is set to zero. + /// Returns the number of vertices inserted. + + template unsigned int delaunay_refinement(FT radius_edge_ratio_bound, ///< radius edge ratio bound (ignored if zero) - FT cell_radius_bound, ///< cell radius bound (ignored if zero) + Sizing_field sizing_field, ///< cell radius bound (ignored if zero) unsigned int max_vertices, ///< number of vertices bound FT enlarge_ratio) ///< bounding box enlarge ratio { - CGAL_TRACE("Calls delaunay_refinement(radius_edge_ratio_bound=%lf, cell_radius_bound=%lf, max_vertices=%u, enlarge_ratio=%lf)\n", - radius_edge_ratio_bound, cell_radius_bound, max_vertices, enlarge_ratio); + return delaunay_refinement(radius_edge_ratio_bound, + sizing_field, + max_vertices, + enlarge_ratio, + internal::Poisson::Constant_sizing_field()); + } + template + unsigned int delaunay_refinement(FT radius_edge_ratio_bound, ///< radius edge ratio bound (ignored if zero) + Sizing_field sizing_field, ///< cell radius bound (ignored if zero) + unsigned int max_vertices, ///< number of vertices bound + FT enlarge_ratio, ///< bounding box enlarge ratio + Second_sizing_field second_sizing_field) + { Sphere elarged_bsphere = enlarged_bounding_sphere(enlarge_ratio); - unsigned int nb_vertices_added = poisson_refine_triangulation(*m_tr,radius_edge_ratio_bound,cell_radius_bound,max_vertices,elarged_bsphere); - - CGAL_TRACE("End of delaunay_refinement()\n"); + unsigned int nb_vertices_added = poisson_refine_triangulation(*m_tr,radius_edge_ratio_bound,sizing_field,second_sizing_field,max_vertices,elarged_bsphere); return nb_vertices_added; } @@ -321,39 +596,46 @@ class Poisson_reconstruction_function double duration_assembly = 0.0; double duration_solve = 0.0; - CGAL_TRACE(" Creates matrix...\n"); + + initialize_cell_indices(); + initialize_barycenters(); // get #variables - unsigned int nb_variables = m_tr->index_unconstrained_vertices(); + constrain_one_vertex_on_convex_hull(); + m_tr->index_unconstrained_vertices(); + unsigned int nb_variables = m_tr->number_of_vertices()-1; - // at least one vertex must be constrained - if(nb_variables == m_tr->number_of_vertices()) - { - constrain_one_vertex_on_convex_hull(); - nb_variables = m_tr->index_unconstrained_vertices(); - } + CGAL_TRACE(" Number of variables: %ld\n", (long)(nb_variables)); // Assemble linear system A*X=B typename SparseLinearAlgebraTraits_d::Matrix A(nb_variables); // matrix is symmetric definite positive typename SparseLinearAlgebraTraits_d::Vector X(nb_variables), B(nb_variables); + initialize_duals(); +#ifndef CGAL_DIV_NON_NORMALIZED + initialize_cell_normals(); +#endif Finite_vertices_iterator v, e; for(v = m_tr->finite_vertices_begin(), e = m_tr->finite_vertices_end(); v != e; ++v) { - if(!v->constrained()) - { + if(!m_tr->is_constrained(v)) { +#ifdef CGAL_DIV_NON_NORMALIZED B[v->index()] = div(v); // rhs -> divergent +#else // not defined(CGAL_DIV_NORMALIZED) + B[v->index()] = div_normalized(v); // rhs -> divergent +#endif // not defined(CGAL_DIV_NORMALIZED) assemble_poisson_row(A,v,B,lambda); } } - + + clear_duals(); + clear_normals(); duration_assembly = (clock() - time_init)/CLOCKS_PER_SEC; CGAL_TRACE(" Creates matrix: done (%.2lf s)\n", duration_assembly); - CGAL_TRACE(" %ld Mb allocated\n", long(CGAL::Memory_sizer().virtual_size()>>20)); CGAL_TRACE(" Solve sparse linear system...\n"); // Solve "A*X = B". On success, solution is (1/D) * X. @@ -369,10 +651,9 @@ class Poisson_reconstruction_function // copy function's values to vertices unsigned int index = 0; for (v = m_tr->finite_vertices_begin(), e = m_tr->finite_vertices_end(); v!= e; ++v) - if(!v->constrained()) + if(!m_tr->is_constrained(v)) v->f() = X[index++]; - CGAL_TRACE(" %ld Mb allocated\n", long(CGAL::Memory_sizer().virtual_size()>>20)); CGAL_TRACE("End of solve_poisson()\n"); return true; @@ -431,16 +712,43 @@ class Poisson_reconstruction_function FT& c, FT& d) const { - const Point& pa = cell->vertex(0)->point(); - const Point& pb = cell->vertex(1)->point(); - const Point& pc = cell->vertex(2)->point(); - const Point& pd = cell->vertex(3)->point(); + // const Point& pa = cell->vertex(0)->point(); + // const Point& pb = cell->vertex(1)->point(); + // const Point& pc = cell->vertex(2)->point(); + const Point& pd = cell->vertex(3)->point(); +#if 1 + //Vector va = pa - pd; + //Vector vb = pb - pd; + //Vector vc = pc - pd; + Vector vp = p - pd; + + //FT i00, i01, i02, i10, i11, i12, i20, i21, i22; + //internal::invert(va.x(), va.y(), va.z(), + // vb.x(), vb.y(), vb.z(), + // vc.x(), vc.y(), vc.z(), + // i00, i01, i02, i10, i11, i12, i20, i21, i22); + const boost::array & i = (*m_Bary)[cell->info()]; + if(i[0]==-1){ + initialize_matrix_entry(cell); + } + // UsedBary[cell->info()] = true; + a = i[0] * vp.x() + i[3] * vp.y() + i[6] * vp.z(); + b = i[1] * vp.x() + i[4] * vp.y() + i[7] * vp.z(); + c = i[2] * vp.x() + i[5] * vp.y() + i[8] * vp.z(); + d = 1 - ( a + b + c); +#else FT v = volume(pa,pb,pc,pd); - a = CGAL::abs(volume(pb,pc,pd,p) / v); - b = CGAL::abs(volume(pa,pc,pd,p) / v); - c = CGAL::abs(volume(pb,pa,pd,p) / v); - d = CGAL::abs(volume(pb,pc,pa,p) / v); + a = std::fabs(volume(pb,pc,pd,p) / v); + b = std::fabs(volume(pa,pc,pd,p) / v); + c = std::fabs(volume(pb,pa,pd,p) / v); + d = std::fabs(volume(pb,pc,pa,p) / v); + + std::cerr << "_________________________________\n"; + std::cerr << aa << " " << bb << " " << cc << " " << dd << std::endl; + std::cerr << a << " " << b << " " << c << " " << d << std::endl; + +#endif } FT find_sink() @@ -484,31 +792,29 @@ class Poisson_reconstruction_function Vertex_handle any_vertex_on_convex_hull() { - // TODO: return NULL if none and assert - std::vector vertices; - vertices.reserve(32); - m_tr->adjacent_vertices(m_tr->infinite_vertex(),std::back_inserter(vertices)); - typename std::vector::iterator it = vertices.begin(); - return *it; + Cell_handle ch = m_tr->infinite_vertex()->cell(); + return ch->vertex( (ch->index( m_tr->infinite_vertex())+1)%4); } + void constrain_one_vertex_on_convex_hull(const FT value = 0.0) { Vertex_handle v = any_vertex_on_convex_hull(); - v->constrained() = true; + m_tr->constrain(v); v->f() = value; } - // divergent - FT div(Vertex_handle v) + // TODO: Some entities are computed too often + // - nn and area should not be computed for the face and its opposite face + // + // divergent + FT div_normalized(Vertex_handle v) { std::vector cells; cells.reserve(32); m_tr->incident_cells(v,std::back_inserter(cells)); - if(cells.size() == 0) - return 0.0; - - FT div = 0.0; + + FT div = 0; typename std::vector::iterator it; for(it = cells.begin(); it != cells.end(); it++) { @@ -517,35 +823,84 @@ class Poisson_reconstruction_function continue; // compute average normal per cell - Vector n = cell_normal(cell); + Vector n = get_cell_normal(cell); // zero normal - no need to compute anything else if(n == CGAL::NULL_VECTOR) continue; + // compute n' int index = cell->index(v); + const Point& x = cell->vertex(index)->point(); const Point& a = cell->vertex((index+1)%4)->point(); const Point& b = cell->vertex((index+2)%4)->point(); const Point& c = cell->vertex((index+3)%4)->point(); Vector nn = (index%2==0) ? CGAL::cross_product(b-a,c-a) : CGAL::cross_product(c-a,b-a); - div += n * nn; + nn = nn / std::sqrt(nn*nn); // normalize + Vector p = a - x; + Vector q = b - x; + Vector r = c - x; + FT p_n = std::sqrt(p*p); + FT q_n = std::sqrt(q*q); + FT r_n = std::sqrt(r*r); + FT solid_angle = p*(CGAL::cross_product(q,r)); + solid_angle = std::abs(solid_angle / (p_n*q_n*r_n + (p*q)*r_n + (q*r)*p_n + (r*p)*q_n)); + + FT area = std::sqrt(squared_area(a,b,c)); + FT length = p_n + q_n + r_n; + div += n * nn * area / length ; + } + return div * FT(3.0); + } + + FT div(Vertex_handle v) + { + std::vector cells; + cells.reserve(32); + m_tr->incident_cells(v,std::back_inserter(cells)); + + FT div = 0.0; + typename std::vector::iterator it; + for(it = cells.begin(); it != cells.end(); it++) + { + Cell_handle cell = *it; + if(m_tr->is_infinite(cell)) + continue; + + const int index = cell->index(v); + const Point& a = cell->vertex(m_tr->vertex_triple_index(index, 0))->point(); + const Point& b = cell->vertex(m_tr->vertex_triple_index(index, 1))->point(); + const Point& c = cell->vertex(m_tr->vertex_triple_index(index, 2))->point(); + const Vector nn = CGAL::cross_product(b-a,c-a); + + div+= nn * (//v->normal() + + cell->vertex((index+1)%4)->normal() + + cell->vertex((index+2)%4)->normal() + + cell->vertex((index+3)%4)->normal()); } return div; } - Vector cell_normal(Cell_handle cell) + Vector get_cell_normal(Cell_handle cell) + { + return Normal[cell->info()]; + } + + Vector cell_normal(Cell_handle cell) const { const Vector& n0 = cell->vertex(0)->normal(); const Vector& n1 = cell->vertex(1)->normal(); const Vector& n2 = cell->vertex(2)->normal(); const Vector& n3 = cell->vertex(3)->normal(); Vector n = n0 + n1 + n2 + n3; - FT sq_norm = n*n; - if(sq_norm != 0.0) - return n / std::sqrt(sq_norm); // normalize - else - return CGAL::NULL_VECTOR; + if(n != NULL_VECTOR){ + FT sq_norm = n*n; + if(sq_norm != 0.0){ + return n / std::sqrt(sq_norm); // normalize + } + } + return NULL_VECTOR; } // cotan formula as area(voronoi face) / len(primal edge) @@ -575,7 +930,7 @@ class Poisson_reconstruction_function { Cell_handle cell = circ; if(!m_tr->is_infinite(cell)) - voronoi_points.push_back(m_tr->dual(cell)); + voronoi_points.push_back(Dual[cell->info()]); else // one infinite tet, switch to another calculation return area_voronoi_face_boundary(edge); circ++; @@ -596,8 +951,7 @@ class Poisson_reconstruction_function { const Point& b = voronoi_points[i]; const Point& c = voronoi_points[i+1]; - Triangle triangle(a,b,c); - area += std::sqrt(triangle.squared_area()); + area += std::sqrt(squared_area(a,b,c)); } return area; } @@ -622,13 +976,14 @@ class Poisson_reconstruction_function if(!m_tr->is_infinite(cell)) { // circumcenter of cell - Point c = m_tr->dual(cell); + Point c = Dual[cell->info()]; Tetrahedron tet = m_tr->tetrahedron(cell); int i = cell->index(vi); int j = cell->index(vj); - int k = -1, l = -1; - other_two_indices(i,j, &k,&l); + int k = Triangulation_utils_3::next_around_edge(i,j); + int l = Triangulation_utils_3::next_around_edge(j,i); + Vertex_handle vk = cell->vertex(k); Vertex_handle vl = cell->vertex(l); @@ -646,11 +1001,8 @@ class Poisson_reconstruction_function Point ck = CGAL::circumcenter(pi,pj,pk); Point cl = CGAL::circumcenter(pi,pj,pl); - Triangle mcck(m,c,ck); - Triangle mccl(m,c,cl); - - area += std::sqrt(mcck.squared_area()); - area += std::sqrt(mccl.squared_area()); + area += std::sqrt(squared_area(m,c,ck)); + area += std::sqrt(squared_area(m,c,cl)); } circ++; } @@ -658,32 +1010,6 @@ class Poisson_reconstruction_function return area; } - // Gets indices different from i and j - void other_two_indices(int i, int j, int* k, int* l) - { - CGAL_surface_reconstruction_points_assertion(i != j); - bool k_done = false; - bool l_done = false; - for(int index=0;index<4;index++) - { - if(index != i && index != j) - { - if(!k_done) - { - *k = index; - k_done = true; - } - else - { - *l = index; - l_done = true; - } - } - } - CGAL_surface_reconstruction_points_assertion(k_done); - CGAL_surface_reconstruction_points_assertion(l_done); - } - /// Assemble vi's row of the linear system A*X=B /// /// @commentheading Template parameters: @@ -700,39 +1026,51 @@ class Poisson_reconstruction_function double diagonal = 0.0; - for(typename std::vector::iterator it = edges.begin(); + for(typename std::vector::iterator it = edges.begin(); it != edges.end(); it++) - { - Vertex_handle vj = it->first->vertex(it->third); - if(vj == vi){ - vj = it->first->vertex(it->second); - } - if(m_tr->is_infinite(vj)) - continue; + { + Vertex_handle vj = it->first->vertex(it->third); + if(vj == vi){ + vj = it->first->vertex(it->second); + } + if(m_tr->is_infinite(vj)) + continue; - // get corresponding edge - Edge edge( it->first, it->first->index(vi), it->first->index(vj)); - if(vi->index() < vj->index()){ - std::swap(edge.second, edge.third); - } + // get corresponding edge + Edge edge( it->first, it->first->index(vi), it->first->index(vj)); + if(vi->index() < vj->index()){ + std::swap(edge.second, edge.third); + } - double cij = cotan_geometric(edge); - if(vj->constrained()) - B[vi->index()] -= cij * vj->f(); // change rhs - else - A.set_coef(vi->index(),vj->index(), -cij, true /*new*/); // off-diagonal coefficient + double cij = cotan_geometric(edge); + + if(m_tr->is_constrained(vj)){ + if(! is_valid(vj->f())){ + std::cerr << "vj->f() = " << vj->f() << " is not valid" << std::endl; + } + B[vi->index()] -= cij * vj->f(); // change rhs + if(! is_valid( B[vi->index()])){ + std::cerr << " B[vi->index()] = " << B[vi->index()] << " is not valid" << std::endl; + } + + } else { + if(! is_valid(cij)){ + std::cerr << "cij = " << cij << " is not valid" << std::endl; + } + A.set_coef(vi->index(),vj->index(), -cij, true /*new*/); // off-diagonal coefficient + } - diagonal += cij; - } + diagonal += cij; + } // diagonal coefficient - if (vi->type() == Triangulation::INPUT) + if (vi->type() == Triangulation::INPUT){ A.set_coef(vi->index(),vi->index(), diagonal + lambda, true /*new*/) ; - else + } else{ A.set_coef(vi->index(),vi->index(), diagonal, true /*new*/); - + } } - + /// Computes enlarged geometric bounding sphere of the embedded triangulation. Sphere enlarged_bounding_sphere(FT ratio) const diff --git a/Surface_reconstruction_points_3/include/CGAL/Reconstruction_triangulation_3.h b/Surface_reconstruction_points_3/include/CGAL/Reconstruction_triangulation_3.h index 672c0b23c9d..650c742a8d9 100644 --- a/Surface_reconstruction_points_3/include/CGAL/Reconstruction_triangulation_3.h +++ b/Surface_reconstruction_points_3/include/CGAL/Reconstruction_triangulation_3.h @@ -27,10 +27,14 @@ #include #include +#include #include -#include +#include #include +#include +#include + #include namespace CGAL { @@ -89,26 +93,21 @@ class Reconstruction_vertex_base_3 : public Vb public: Reconstruction_vertex_base_3() - : Vb(), m_f(FT(0.0)), m_constrained(false), m_type(0), m_index(0) + : Vb(), m_f(FT(0.0)), m_type(0), m_index(0) {} Reconstruction_vertex_base_3(const Point_with_normal& p) - : Vb(p), m_f(FT(0.0)), m_constrained(false), m_type(0), m_index(0) + : Vb(p), m_f(FT(0.0)), m_type(0), m_index(0) {} Reconstruction_vertex_base_3(const Point_with_normal& p, Cell_handle c) - : Vb(p,c), m_f(FT(0.0)), m_constrained(false), m_type(0), m_index(0) + : Vb(p,c), m_f(FT(0.0)), m_type(0), m_index(0) {} Reconstruction_vertex_base_3(Cell_handle c) - : Vb(c), m_f(FT(0.0)), m_constrained(false), m_type(0), m_index(0) + : Vb(c), m_f(FT(0.0)), m_type(0), m_index(0) {} - /// Is vertex constrained, i.e. - /// does it contribute to the right or left member of the linear system? - /// Default value is false. - bool constrained() const { return m_constrained; } - bool& constrained() { return m_constrained; } /// Gets/sets the value of the implicit function. /// Default value is 0.0. @@ -169,7 +168,7 @@ struct Reconstruction_triangulation_default_geom_traits_3 : public BaseGt template , - class Tds_ = Triangulation_data_structure_3 > > + class Tds_ = Triangulation_data_structure_3, Triangulation_cell_base_with_info_3 > > class Reconstruction_triangulation_3 : public Delaunay_triangulation_3 { // Private types @@ -236,8 +235,8 @@ class Reconstruction_triangulation_3 : public Delaunay_triangulation_3 /// Point type enum Point_type { - INPUT, ///< Input point. - STEINER ///< Steiner point created by Delaunay refinement. + INPUT=0, ///< Input point. + STEINER=1 ///< Steiner point created by Delaunay refinement. }; /// Iterator over input vertices. @@ -248,13 +247,21 @@ class Reconstruction_triangulation_3 : public Delaunay_triangulation_3 typedef Iterator_project > Input_point_iterator; -// Public methods + mutable Sphere sphere; + std::vector points; + std::size_t fraction; + std::list fractions; + Vertex_handle constrained_vertex; + + public: /// Default constructor. Reconstruction_triangulation_3() - { - } + {} + + ~Reconstruction_triangulation_3() + {} // Default copy constructor and operator =() are fine. @@ -294,62 +301,55 @@ class Reconstruction_triangulation_3 : public Delaunay_triangulation_3 return Input_point_iterator(input_vertices_end()); } - /// Gets the bounding sphere of all points. - Sphere bounding_sphere() const - { - typedef Min_sphere_of_spheres_d_traits_3 Traits; - typedef Min_sphere_of_spheres_d Min_sphere; - typedef typename Traits::Sphere Traits_sphere; - - // Represents *all* points by a set of spheres with 0 radius - std::vector spheres; - spheres.reserve(number_of_vertices()); + /// Gets the bounding sphere of input points. - for (Point_iterator it=points_begin(), eit=points_end(); - it != eit; ++it) - spheres.push_back(Traits_sphere(*it,0)); - // Computes min sphere - Min_sphere ms(spheres.begin(),spheres.end()); - typename Min_sphere::Cartesian_const_iterator coord = ms.center_cartesian_begin(); - FT cx = *coord++; - FT cy = *coord++; - FT cz = *coord++; - return Sphere(Point(cx,cy,cz), ms.radius()*ms.radius()); + Sphere bounding_sphere() const + { + return sphere; } - /// Gets the bounding sphere of input points. - Sphere input_points_bounding_sphere() const + void initialize_bounding_sphere() const { - typedef Min_sphere_of_spheres_d_traits_3 Traits; + typedef Min_sphere_of_points_d_traits_3 Traits; typedef Min_sphere_of_spheres_d Min_sphere; - typedef typename Traits::Sphere Traits_sphere; - - // Represents *input* points by a set of spheres with 0 radius - std::vector spheres; - for (Input_point_iterator it=input_points_begin(), eit=input_points_end(); - it != eit; - ++it) - spheres.push_back(Traits_sphere(*it,0)); - + // Computes min sphere - Min_sphere ms(spheres.begin(),spheres.end()); + Min_sphere ms(points.begin(),points.end()); + typename Min_sphere::Cartesian_const_iterator coord = ms.center_cartesian_begin(); FT cx = *coord++; FT cy = *coord++; - FT cz = *coord++; - return Sphere(Point(cx,cy,cz), ms.radius()*ms.radius()); + FT cz = *coord; + + sphere = Sphere(Point(cx,cy,cz), ms.radius()*ms.radius()); } /// Insert point in the triangulation. /// Default type is INPUT. + template Vertex_handle insert(const Point_with_normal& p, - Point_type type = INPUT, - Cell_handle start = Cell_handle()) + Point_type type,// = INPUT, + Cell_handle start,// = Cell_handle(), + Visitor visitor) { - Vertex_handle v = Base::insert(p, start); + + if(type == INPUT){ + visitor.before_insertion(); + } + if(this->dimension() < 3){ + Vertex_handle v = Base::insert(p, start); + v->type() = type; + return v; + } + typename Base::Locate_type lt; + int li, lj; + Cell_handle ch = Base::locate(p, lt, li, lj, start); + + Vertex_handle v = Base::insert(p, lt, ch, li, lj); v->type() = type; return v; + } /// Insert the [first, beyond) range of points in the triangulation using a spatial sort. @@ -366,57 +366,94 @@ class Reconstruction_triangulation_3 : public Delaunay_triangulation_3 // This variant requires all parameters. template int insert( InputIterator first, ///< iterator over the first input point. InputIterator beyond, ///< past-the-end iterator over the input points. PointPMap point_pmap, ///< property map to access the position of an input point. NormalPMap normal_pmap, ///< property map to access the *oriented* normal of an input point. - Point_type type = INPUT) + Point_type type, // = INPUT, + Visitor visitor) { - int n = number_of_vertices(); - + if(! points.empty()){ + std::cerr << "WARNING: not all points inserted yet" << std::endl; + } // Convert input points to Point_with_normal_3 - std::vector points; + //std::vector points; for (InputIterator it = first; it != beyond; ++it) { Point_with_normal pwn(get(point_pmap,it), get(normal_pmap,it)); points.push_back(pwn); } + int n = points.size(); + + initialize_bounding_sphere(); + + boost::rand48 random; + boost::random_number_generator rng(random); + std::random_shuffle (points.begin(), points.end(), rng); + fraction = 0; + + fractions.clear(); + fractions.push_back(1.0); + + double m = n; + + while(m > 500){ + m /= 2; + fractions.push_front(m/n); + } + + insert_fraction(visitor); + return 0; + } - // Spatial sorting - std::random_shuffle (points.begin(), points.end()); - spatial_sort (points.begin(), points.end(), geom_traits()); - - // Insert in triangulation + template + bool insert_fraction(Visitor visitor) + { + if(fractions.empty()){ + points.clear(); + return false; + } + double frac = fractions.front(); + fractions.pop_front(); + std::size_t more = (std::size_t)(points.size() * frac) - fraction; + if((fraction+more) > points.size()){ + more = points.size() - fraction; + } Cell_handle hint; - for (typename std::vector::const_iterator p = points.begin(); - p != points.end(); ++p) + spatial_sort (points.begin()+fraction, points.begin()+fraction+more, geom_traits()); + for (typename std::vector::const_iterator p = points.begin()+fraction; + p != points.begin()+fraction+more; ++p) { - Vertex_handle v = insert(*p, type, hint); + Vertex_handle v = insert(*p, INPUT, hint, visitor); hint = v->cell(); } - - return number_of_vertices() - n; + fraction += more; + return true; } /// @cond SKIP_IN_MANUAL // This variant creates a default point property map = Dereference_property_map. template int insert( InputIterator first, ///< iterator over the first input point. InputIterator beyond, ///< past-the-end iterator over the input points. NormalPMap normal_pmap, ///< property map to access the *oriented* normal of an input point. - Point_type type = INPUT) + Point_type type,// = INPUT, + Visitor visitor) { return insert( first,beyond, make_dereference_property_map(first), normal_pmap, - type); + type, + visitor); } /// Delaunay refinement callback: @@ -442,12 +479,26 @@ class Reconstruction_triangulation_3 : public Delaunay_triangulation_3 v!= e; ++v) { - if(!v->constrained()) + if(! is_constrained(v)) v->index() = index++; } return index; } + /// Is vertex constrained, i.e. + /// does it contribute to the right or left member of the linear system? + + bool is_constrained(Vertex_handle v) const + { + return v == constrained_vertex; + } + + void constrain(Vertex_handle v) + { + constrained_vertex = v; + } + + }; // end of Reconstruction_triangulation_3 diff --git a/Surface_reconstruction_points_3/include/CGAL/poisson_refine_triangulation.h b/Surface_reconstruction_points_3/include/CGAL/poisson_refine_triangulation.h index 8fd2f65c62a..6e487d4e964 100644 --- a/Surface_reconstruction_points_3/include/CGAL/poisson_refine_triangulation.h +++ b/Surface_reconstruction_points_3/include/CGAL/poisson_refine_triangulation.h @@ -67,7 +67,6 @@ class Poisson_mesher_level_impl_base : Poisson_mesher_level_impl_base(Tr& t, Criteria crit, unsigned int max_vert, Surface& surface, Oracle& oracle) : Base(t, crit, surface, oracle), max_vertices(max_vert) ///< number of vertices bound (ignored if zero) - {} protected: @@ -192,15 +191,17 @@ class Poisson_mesher_level : /// @param Tr 3D Delaunay triangulation. /// @param Surface Sphere_3 or Iso_cuboid_3. template + typename Surface, + typename Sizing_field, + typename Second_sizing_field> unsigned int poisson_refine_triangulation( Tr& tr, double radius_edge_ratio_bound, ///< radius edge ratio bound (>= 1.0) - double cell_radius_bound, ///< cell radius bound (ignored if zero) + const Sizing_field& sizing_field, ///< cell radius bound (ignored if zero) + const Second_sizing_field& second_sizing_field, ///< cell radius bound (ignored if zero) unsigned int max_vertices, ///< number of vertices bound (ignored if zero) Surface& enlarged_bbox) ///< new bounding sphere or box { - CGAL_TRACE("Calls poisson_refine_triangulation()\n"); // Convergence is guaranteed if radius_edge_ratio_bound >= 1.0 CGAL_surface_reconstruction_points_precondition(radius_edge_ratio_bound >= 1.0); @@ -211,16 +212,17 @@ unsigned int poisson_refine_triangulation( typedef typename Gt::Point_3 Point; // Mesher_level types - typedef Poisson_mesh_cell_criteria_3 Tets_criteria; + typedef Poisson_mesh_cell_criteria_3 Tets_criteria; typedef typename CGAL::Surface_mesh_traits_generator_3::type Oracle; typedef Poisson_mesher_level Refiner; - CGAL_TRACE(" Creates queue\n"); int nb_vertices = tr.number_of_vertices(); // get former #vertices // Delaunay refinement - Tets_criteria tets_criteria(radius_edge_ratio_bound*radius_edge_ratio_bound, cell_radius_bound); + Tets_criteria tets_criteria(radius_edge_ratio_bound*radius_edge_ratio_bound, + sizing_field, + second_sizing_field); Oracle oracle; Null_mesher_level null_mesher_level; Refiner refiner(tr, tets_criteria, max_vertices, enlarged_bbox, oracle, null_mesher_level); @@ -229,7 +231,6 @@ unsigned int poisson_refine_triangulation( int nb_vertices_added = tr.number_of_vertices() - nb_vertices; - CGAL_TRACE("End of poisson_refine_triangulation()\n"); return (unsigned int) nb_vertices_added; }