diff --git a/Arrangement_on_surface_2/doc/Arrangement_on_surface_2/Arrangement_on_surface_2.txt b/Arrangement_on_surface_2/doc/Arrangement_on_surface_2/Arrangement_on_surface_2.txt index ad8d8120b55..dc4d957e931 100644 --- a/Arrangement_on_surface_2/doc/Arrangement_on_surface_2/Arrangement_on_surface_2.txt +++ b/Arrangement_on_surface_2/doc/Arrangement_on_surface_2/Arrangement_on_surface_2.txt @@ -2681,7 +2681,7 @@ manner whenever possible. Thus, it resorts to exact computations only when the approximate computation fails to produce an unambiguous result. Note that most arrangement vertices are therefore associated with approximated points. You cannot access the coordinates of such points and obtain them as -algebraic numbers, and only access to the approximate coordinates in possible. +algebraic numbers, and only access to the approximate coordinates is possible. See the Reference Manual for the exact interface of the `Point_2`, `Curve_2` and `X_monotone_curve_2` defined by the traits class. diff --git a/Arrangement_on_surface_2/include/CGAL/Arr_spherical_gaussian_map_3/Arr_polyhedral_sgm.h b/Arrangement_on_surface_2/include/CGAL/Arr_spherical_gaussian_map_3/Arr_polyhedral_sgm.h index b2be5687c10..8a57cd28808 100644 --- a/Arrangement_on_surface_2/include/CGAL/Arr_spherical_gaussian_map_3/Arr_polyhedral_sgm.h +++ b/Arrangement_on_surface_2/include/CGAL/Arr_spherical_gaussian_map_3/Arr_polyhedral_sgm.h @@ -242,6 +242,11 @@ private: /*! Set the marked-face index */ void set_marked_facet_index(size_type id) {m_marked_facet_index = id;} + /*! Add vertices to the current facet. */ + template + void add_vertices_to_facet(Iterator begin, Iterator end, Builder& B) + { for (Iterator it = begin; it != end; ++it) B.add_vertex_to_facet(*it); } + /*! builds the polyhedron */ void operator()(HDS& hds) { @@ -262,11 +267,11 @@ private: for (CoordIndexIter it = m_indices_begin; it != m_indices_end; ++it) { Polyhedron_facet_handle fh = B.begin_facet(); if (counter == m_marked_facet_index) fh->set_marked(true); - //! \todo EF: when upgrading to C++11 change the type of the following - // iterator to auto. Better yet, use for (auto blah : foo). - for (std::vector::const_iterator iit = it->begin(); - iit != it->end(); ++iit) - B.add_vertex_to_facet(*iit); + //! \todo EF: when upgrading to C++11 enable the following code and + // remove add_vertices_to_facet(). + // for (const auto& facet : *it) B.add_vertex_to_facet(facet); + // B.end_facet(); + add_vertices_to_facet(it->begin(), it->end(), B); B.end_facet(); ++counter; } diff --git a/BGL/include/CGAL/boost/graph/METIS/partition_dual_graph.h b/BGL/include/CGAL/boost/graph/METIS/partition_dual_graph.h index 8bfeb002ae1..7e2a15870a7 100644 --- a/BGL/include/CGAL/boost/graph/METIS/partition_dual_graph.h +++ b/BGL/include/CGAL/boost/graph/METIS/partition_dual_graph.h @@ -36,12 +36,15 @@ #include #include +#include + namespace CGAL { namespace METIS { template -void partition_dual_graph(const TriangleMesh& tm, int nparts, +void partition_dual_graph(const TriangleMesh& tm, + int nparts, METIS_options options, // options array const NamedParameters& np) { @@ -93,11 +96,11 @@ void partition_dual_graph(const TriangleMesh& tm, int nparts, idx_t objval; // partition info for the nodes - idx_t* npart = (idx_t*) calloc(nn, sizeof(idx_t)); + idx_t* npart = (idx_t*) calloc(num_vertices(tm), sizeof(idx_t)); CGAL_assertion(npart != NULL); // partition info for the elements - idx_t* epart = (idx_t*) calloc(ne, sizeof(idx_t)); + idx_t* epart = (idx_t*) calloc(num_faces(tm), sizeof(idx_t)); CGAL_assertion(epart != NULL); // do not support Fortran-style arrays @@ -118,6 +121,12 @@ void partition_dual_graph(const TriangleMesh& tm, int nparts, Output_face_partition_ids fo; vo(tm, indices, npart, get_param(np, internal_np::vertex_partition_id)); fo(tm, epart, get_param(np, internal_np::face_partition_id)); + + delete[] eptr; + delete[] eind; + + std::free(npart); + std::free(epart); } template diff --git a/BGL/include/CGAL/boost/graph/METIS/partition_graph.h b/BGL/include/CGAL/boost/graph/METIS/partition_graph.h index 65176d37131..d0f3e56b759 100644 --- a/BGL/include/CGAL/boost/graph/METIS/partition_graph.h +++ b/BGL/include/CGAL/boost/graph/METIS/partition_graph.h @@ -34,6 +34,8 @@ #include #include +#include + namespace CGAL { namespace METIS { @@ -76,7 +78,8 @@ struct Output_face_partition_ids }; template -void partition_graph(const TriangleMesh& tm, int nparts, +void partition_graph(const TriangleMesh& tm, + int nparts, METIS_options options, // pointer to the options array const NamedParameters& np) { @@ -125,11 +128,11 @@ void partition_graph(const TriangleMesh& tm, int nparts, idx_t objval; // partition info for the nodes - idx_t* npart = (idx_t*) calloc(nn, sizeof(idx_t)); + idx_t* npart = (idx_t*) calloc(num_vertices(tm), sizeof(idx_t)); CGAL_assertion(npart != NULL); // partition info for the elements - idx_t* epart = (idx_t*) calloc(ne, sizeof(idx_t)); + idx_t* epart = (idx_t*) calloc(num_faces(tm), sizeof(idx_t)); CGAL_assertion(epart != NULL); // do not support Fortran-style arrays @@ -150,6 +153,12 @@ void partition_graph(const TriangleMesh& tm, int nparts, Output_face_partition_ids fo; vo(tm, indices, npart, get_param(np, internal_np::vertex_partition_id)); fo(tm, epart, get_param(np, internal_np::face_partition_id)); + + delete[] eptr; + delete[] eind; + + std::free(npart); + std::free(epart); } template diff --git a/BGL/include/CGAL/boost/graph/named_params_helper.h b/BGL/include/CGAL/boost/graph/named_params_helper.h index cbf5a2fedce..e585269eba0 100644 --- a/BGL/include/CGAL/boost/graph/named_params_helper.h +++ b/BGL/include/CGAL/boost/graph/named_params_helper.h @@ -405,6 +405,29 @@ namespace CGAL { > ::type type; }; + template + class GetIsConstrainedMap + { + struct DummyConstrainedMap + { + typedef typename std::iterator_traits::value_type key_type; + typedef bool value_type; + typedef value_type reference; + typedef boost::readable_property_map_tag category; + + typedef DummyConstrainedMap Self; + friend reference get(const Self&, const key_type&) { return false; } + }; + + public: + typedef DummyConstrainedMap NoMap; + typedef typename boost::lookup_named_param_def < + internal_np::point_is_constrained_t, + NamedParameters, + DummyConstrainedMap //default + > ::type type; + }; + } // namespace Point_set_processing_3 template diff --git a/BGL/include/CGAL/boost/graph/parameters_interface.h b/BGL/include/CGAL/boost/graph/parameters_interface.h index d2dd2208420..9e38eec688d 100644 --- a/BGL/include/CGAL/boost/graph/parameters_interface.h +++ b/BGL/include/CGAL/boost/graph/parameters_interface.h @@ -117,3 +117,4 @@ CGAL_add_named_parameter(plane_t, plane_map, plane_map) CGAL_add_named_parameter(plane_index_t, plane_index_map, plane_index_map) CGAL_add_named_parameter(select_percentage_t, select_percentage, select_percentage) CGAL_add_named_parameter(require_uniform_sampling_t, require_uniform_sampling, require_uniform_sampling) +CGAL_add_named_parameter(point_is_constrained_t, point_is_constrained, point_is_constrained_map) diff --git a/Installation/CHANGES.md b/Installation/CHANGES.md index 86dd48044b6..9f2c4f38d50 100644 --- a/Installation/CHANGES.md +++ b/Installation/CHANGES.md @@ -32,6 +32,14 @@ Release date: March 2019 - Added the class `CGAL::Rigid_triangle_mesh_collision_detection` to detect intersections between meshes and volumes undergoing affine transformations. +### Point Set Processing + +- `CGAL::mst_orient_normals()` can now be called with a set of user-selected + seed points that are known to be already oriented. A new optional named + parameter `point_is_constrained_map` is added for this purpose. The + original behavior (using one unique and automatically selected seed) is + kept if this parameter is not used. + ### 3D Fast Intersection and Distance Computation - The primitives `AABB_face_graph_triangle_primitive` and diff --git a/Mesh_3/include/CGAL/Mesh_3/Mesher_level.h b/Mesh_3/include/CGAL/Mesh_3/Mesher_level.h index 932de589365..c5554f3e109 100644 --- a/Mesh_3/include/CGAL/Mesh_3/Mesher_level.h +++ b/Mesh_3/include/CGAL/Mesh_3/Mesher_level.h @@ -1048,6 +1048,11 @@ public: { // Lock the element area on the grid Element element = derivd.extract_element_from_container_value(ce); + + // This is safe to do with the concurrent compact container because even if the element `ce` + // gets removed from the TDS at this point, it is not actually deleted in the cells container and + // `ce->vertex(0-3)` still points to a vertex of the vertices container whose `.point()` + // can be safely accessed (even if that vertex has itself also been removed from the TDS). bool locked = derivd.try_lock_element(element, FIRST_GRID_LOCK_RADIUS); if( locked ) diff --git a/Nef_3/include/CGAL/boost/graph/convert_nef_polyhedron_to_polygon_mesh.h b/Nef_3/include/CGAL/boost/graph/convert_nef_polyhedron_to_polygon_mesh.h index 369db23902d..7393881de00 100644 --- a/Nef_3/include/CGAL/boost/graph/convert_nef_polyhedron_to_polygon_mesh.h +++ b/Nef_3/include/CGAL/boost/graph/convert_nef_polyhedron_to_polygon_mesh.h @@ -343,36 +343,45 @@ void collect_polygon_mesh_info( } //end of namespace nef_to_pm - -template -void convert_nef_polyhedron_to_polygon_mesh(const Nef_polyhedron& nef, Polygon_mesh& pm, bool triangulate_all_faces = false) +template +void convert_nef_polyhedron_to_polygon_soup(const Nef_polyhedron& nef, + std::vector& points, + std::vector< std::vector >& polygons, + bool triangulate_all_faces = false) { typedef typename Nef_polyhedron::Point_3 Point_3; - typedef typename boost::property_traits::type>::value_type PM_Point; - typedef typename Kernel_traits::Kernel PM_Kernel; typedef typename Kernel_traits::Kernel Nef_Kernel; - typedef Cartesian_converter Converter; - + typedef Cartesian_converter Converter; + typedef typename Output_kernel::Point_3 Out_point; typename Nef_polyhedron::Volume_const_iterator vol_it = nef.volumes_begin(), vol_end = nef.volumes_end(); if ( Nef_polyhedron::Infi_box::extended_kernel() ) ++vol_it; // skip Infi_box CGAL_assertion ( vol_it != vol_end ); ++vol_it; // skip unbounded volume + Converter to_output; + for (;vol_it!=vol_end;++vol_it) + nef_to_pm::collect_polygon_mesh_info(points, + polygons, + nef, + vol_it->shells_begin(), + to_output, + triangulate_all_faces); +} + +template +void convert_nef_polyhedron_to_polygon_mesh(const Nef_polyhedron& nef, Polygon_mesh& pm, bool triangulate_all_faces = false) +{ + typedef typename boost::property_traits::type>::value_type PM_Point; + typedef typename Kernel_traits::Kernel PM_Kernel; + std::vector points; std::vector< std::vector > polygons; - Converter to_inexact; - for (;vol_it!=vol_end;++vol_it) - nef_to_pm::collect_polygon_mesh_info(points, - polygons, - nef, - vol_it->shells_begin(), - to_inexact, - triangulate_all_faces); - + convert_nef_polyhedron_to_polygon_soup(nef, points, polygons, triangulate_all_faces); Polygon_mesh_processing::polygon_soup_to_polygon_mesh(points, polygons, pm); } + } //end of namespace CGAL diff --git a/Point_set_3/include/CGAL/Point_set_3/IO.h b/Point_set_3/include/CGAL/Point_set_3/IO.h index 6c534f4c84c..7ae5b46dafc 100644 --- a/Point_set_3/include/CGAL/Point_set_3/IO.h +++ b/Point_set_3/include/CGAL/Point_set_3/IO.h @@ -394,7 +394,10 @@ read_ply_point_set( internal::PLY::Point_set_3_filler filler(point_set); if (!(reader.init (stream))) + { + stream.setstate(std::ios::failbit); return false; + } if (comments != NULL) *comments = reader.comments(); @@ -415,8 +418,7 @@ read_ply_point_set( { internal::PLY::PLY_read_number* property = element.property(k); property->get (stream); - - if (stream.eof()) + if (stream.fail()) return false; } diff --git a/Point_set_processing_3/doc/Point_set_processing_3/NamedParameters.txt b/Point_set_processing_3/doc/Point_set_processing_3/NamedParameters.txt index dcb0f56f626..ee321dfd95e 100644 --- a/Point_set_processing_3/doc/Point_set_processing_3/NamedParameters.txt +++ b/Point_set_processing_3/doc/Point_set_processing_3/NamedParameters.txt @@ -31,17 +31,17 @@ typename CGAL::Kernel_traits::Kernel \cgalNPEnd \cgalNPBegin{point_map} \anchor PSP_point_map - is the property map containing the points associated to the iterators of the point range `points`.\n + is the property map containing the points associated to the elements of the point range `points`.\n \b Type: a class model of `ReadablePropertyMap` with -`PointRange::iterator` as key type and +`PointRange::iterator::value_type` as key type and `geom_traits::Point_3` as value type. \n Default value: \code CGAL::Identity_property_map\endcode \cgalNPEnd \cgalNPBegin{normal_map} \anchor PSP_normal_map - is the property map containing the normal vectors associated to the iterators of the point range `points`.\n + is the property map containing the normal vectors associated to the elements of the point range `points`.\n \b Type: a class model of `ReadablePropertyMap` with -`PointRange::iterator` as key type and +`PointRange::iterator::value_type` as key type and `geom_traits::Vector_3` as value type. \n No default value. \cgalNPEnd @@ -67,9 +67,9 @@ Note that when a callback is run on a parallelized algorithm with `CGAL::Paralle \cgalNPEnd \cgalNPBegin{query_point_map} \anchor PSP_query_point_map - is the property map containing the points associated to the iterators of the point range `queries`.\n + is the property map containing the points associated to the elements of the point range `queries`.\n \b Type: a class model of `ReadablePropertyMap` with -`PointRange::iterator` as key type and +`PointRange::iterator::value_type` as key type and `geom_traits::Point_3` as value type. \n Default value: \code CGAL::Identity_property_map\endcode \cgalNPEnd @@ -146,9 +146,9 @@ multiple of a tolerance `epsilon` used to connect simplices. \cgalNPEnd \cgalNPBegin{plane_map} \anchor PSP_plane_map - is the property map containing the planes associated to the iterators of the plane range `planes`.\n + is the property map containing the planes associated to the elements of the plane range `planes`.\n \b Type: a class model of `ReadablePropertyMap` with -`PlaneRange::iterator` as key type and +`PlaneRange::iterator::value_type` as key type and `geom_traits::Plane_3` as value type. \n Default value: \code CGAL::Identity_property_map\endcode \cgalNPEnd @@ -182,6 +182,15 @@ give better result if the distribution of the input points is highly non-uniform Default value: `false` \cgalNPEnd +\cgalNPBegin{point_is_constrained_map} \anchor PSP_point_is_constrained_map +is the property map containing information about points being constrained or not. +Constrained points are left unaltered and are used as seeds in `mst_orient_normals()`.\n +\b Type: a class model of `ReadablePropertyMap` with +`PointRange::iterator::value_type` as key type and +`bool` as value type. \n +Default value: a property map with only the highest point constrained. +\cgalNPEnd + \cgalNPTableEnd */ diff --git a/Point_set_processing_3/include/CGAL/IO/read_ply_points.h b/Point_set_processing_3/include/CGAL/IO/read_ply_points.h index 69202810990..1dd14548a94 100644 --- a/Point_set_processing_3/include/CGAL/IO/read_ply_points.h +++ b/Point_set_processing_3/include/CGAL/IO/read_ply_points.h @@ -765,7 +765,10 @@ bool read_ply_points_with_properties (std::istream& stream, internal::PLY::PLY_reader reader; if (!(reader.init (stream))) + { + stream.setstate(std::ios::failbit); return false; + } for (std::size_t i = 0; i < reader.number_of_elements(); ++ i) { @@ -778,7 +781,7 @@ bool read_ply_points_with_properties (std::istream& stream, internal::PLY::PLY_read_number* property = element.property(k); property->get (stream); - if (stream.eof()) + if (stream.fail()) return false; } diff --git a/Point_set_processing_3/include/CGAL/mst_orient_normals.h b/Point_set_processing_3/include/CGAL/mst_orient_normals.h index 4935cd7d208..72cd4e3b48f 100644 --- a/Point_set_processing_3/include/CGAL/mst_orient_normals.h +++ b/Point_set_processing_3/include/CGAL/mst_orient_normals.h @@ -121,6 +121,34 @@ public: const NormalMap m_normal_map; }; +template +class Default_constrained_map +{ +public: + + typedef boost::readable_property_map_tag category; + typedef typename ForwardIterator::value_type key_type; + typedef bool value_type; + typedef value_type reference; + +private: + + ForwardIterator m_source_point; + +public: + + Default_constrained_map () { } + Default_constrained_map (ForwardIterator source_point) + : m_source_point (source_point) { } + + /// Free function to access the map elements. + friend inline + reference get(const Default_constrained_map& map, key_type p) + { + return (p == *map.m_source_point); + } + +}; /// Helper class: Propagate_normal_orientation /// @@ -145,10 +173,12 @@ struct Propagate_normal_orientation : public boost::base_visitor< Propagate_normal_orientation > { typedef internal::MST_graph MST_graph; + typedef typename MST_graph::vertex_descriptor vertex_descriptor; typedef boost::on_examine_edge event_filter; - Propagate_normal_orientation(double angle_max = CGAL_PI/2.) ///< max angle to propagate the normal orientation (radians) - : m_angle_max(angle_max) + Propagate_normal_orientation(vertex_descriptor source, + double angle_max = CGAL_PI/2.) ///< max angle to propagate the normal orientation (radians) + : m_source(source), m_angle_max(angle_max) { // Precondition: 0 < angle_max <= PI/2 CGAL_point_set_processing_precondition(0 < angle_max && angle_max <= CGAL_PI/2.); @@ -158,16 +188,28 @@ struct Propagate_normal_orientation void operator()(Edge& edge, const MST_graph& mst_graph) { typedef typename boost::property_traits::reference Vector_ref; - typedef typename MST_graph::vertex_descriptor vertex_descriptor; + + // Gets source + vertex_descriptor source_vertex = source(edge, mst_graph); + + // Gets target + vertex_descriptor target_vertex = target(edge, mst_graph); + bool& target_normal_is_oriented = ((MST_graph&)mst_graph)[target_vertex].is_oriented; + + // special case if vertex is source vertex (and thus has no related point/normal) + if (source_vertex == m_source) + { + target_normal_is_oriented = true; + return; + } // Gets source normal - vertex_descriptor source_vertex = source(edge, mst_graph); Vector_ref source_normal = get(mst_graph.m_normal_map, *(mst_graph[source_vertex].input_point) ); const bool source_normal_is_oriented = mst_graph[source_vertex].is_oriented; - // Gets target normal - vertex_descriptor target_vertex = target(edge, mst_graph); + + // Gets target Vector_ref target_normal = get( mst_graph.m_normal_map, *(mst_graph[target_vertex].input_point) ); - bool& target_normal_is_oriented = ((MST_graph&)mst_graph)[target_vertex].is_oriented; + if ( ! target_normal_is_oriented ) { // -> -> @@ -188,6 +230,7 @@ struct Propagate_normal_orientation // Data // Implementation note: boost::breadth_first_search() makes copies of this object => data must be constant or shared. private: + vertex_descriptor m_source; const double m_angle_max; ///< max angle to propagate the normal orientation (radians). }; @@ -264,6 +307,7 @@ template Riemannian_graph @@ -273,6 +317,7 @@ create_riemannian_graph( PointMap point_map, ///< property map: value_type of ForwardIterator -> Point_3 NormalMap normal_map, ///< property map: value_type of ForwardIterator -> Vector_3 IndexMap index_map, ///< property map ForwardIterator -> index + ConstrainedMap constrained_map, ///< property map ForwardIterator -> bool unsigned int k, ///< number of neighbors const Kernel& /*kernel*/) ///< geometric traits. { @@ -337,6 +382,11 @@ create_riemannian_graph( CGAL_point_set_processing_assertion(v == get(index_map,it)); riemannian_graph[v].input_point = it; } + + // add source vertex (virtual, does not correspond to a point) + add_vertex(riemannian_graph); + std::size_t source_point_index = num_input_points; + // // add edges Riemannian_graph_weight_map riemannian_graph_weight_map = get(boost::edge_weight, riemannian_graph); @@ -384,6 +434,20 @@ create_riemannian_graph( search_iterator++; } + + // Check if point is source + if (get(constrained_map, *it)) + { + typename boost::graph_traits::edge_descriptor e; + bool inserted; + boost::tie(e, inserted) = add_edge(vertex(it_index, riemannian_graph), + vertex(source_point_index, riemannian_graph), + riemannian_graph); + CGAL_point_set_processing_assertion(inserted); + + riemannian_graph_weight_map[e] = 0.; + } + } return riemannian_graph; @@ -418,8 +482,7 @@ create_mst_graph( IndexMap index_map, ///< property map ForwardIterator -> index unsigned int k, ///< number of neighbors const Kernel& kernel, ///< geometric traits. - const Riemannian_graph& riemannian_graph, ///< graph connecting each vertex to its knn - ForwardIterator source_point) ///< source point (with an oriented normal) + const Riemannian_graph& riemannian_graph) ///< graph connecting each vertex to its knn { // prevents warnings CGAL_USE(point_map); @@ -440,16 +503,17 @@ create_mst_graph( CGAL_point_set_processing_precondition(first != beyond); // Number of input points - const std::size_t num_input_points = num_vertices(riemannian_graph); + const std::size_t num_input_points = num_vertices(riemannian_graph) - 1; std::size_t memory = CGAL::Memory_sizer().virtual_size(); CGAL_TRACE(" %ld Mb allocated\n", memory>>20); CGAL_TRACE(" Calls boost::prim_minimum_spanning_tree()\n"); // Computes Minimum Spanning Tree. - std::size_t source_point_index = get(index_map, source_point); + std::size_t source_point_index = num_input_points; + Riemannian_graph_weight_map riemannian_graph_weight_map = get(boost::edge_weight, riemannian_graph); typedef std::vector PredecessorMap; - PredecessorMap predecessor(num_input_points); + PredecessorMap predecessor(num_input_points + 1); boost::prim_minimum_spanning_tree(riemannian_graph, &predecessor[0], weight_map( riemannian_graph_weight_map ) .root_vertex( vertex(source_point_index, riemannian_graph) )); @@ -472,8 +536,13 @@ create_mst_graph( typename MST_graph::vertex_descriptor v = add_vertex(mst_graph); CGAL_point_set_processing_assertion(v == get(index_map,it)); mst_graph[v].input_point = it; - mst_graph[v].is_oriented = (it == source_point); + mst_graph[v].is_oriented = false; } + + typename MST_graph::vertex_descriptor v = add_vertex(mst_graph); + CGAL_point_set_processing_assertion(v == source_point_index); + mst_graph[v].is_oriented = true; + // add edges for (std::size_t i=0; i < predecessor.size(); i++) // add edges { @@ -525,6 +594,11 @@ create_mst_graph( If this parameter is omitted, `CGAL::Identity_property_map` is used.\cgalParamEnd \cgalParamBegin{normal_map} a model of `ReadWritePropertyMap` with value type `geom_traits::Vector_3`.\cgalParamEnd + \cgalParamBegin{point_is_constrained_map} a model of `ReadablePropertyMap` with value type + `bool`. Points with a `true` value will be used as seed points: their normal will be considered as already + oriented, it won't be altered and it will be propagated to its neighbors. If this parameter is omitted, + the highest point (highest Z coordinate) will be used as the unique seed with an upward oriented + normal\cgalParamEnd \cgalParamBegin{geom_traits} an instance of a geometric traits class, model of `Kernel`\cgalParamEnd \cgalNamedParamsEnd @@ -545,6 +619,7 @@ mst_orient_normals( typedef typename Point_set_processing_3::GetPointMap::type PointMap; typedef typename Point_set_processing_3::GetNormalMap::type NormalMap; typedef typename Point_set_processing_3::GetK::Kernel Kernel; + typedef typename Point_set_processing_3::GetIsConstrainedMap::type ConstrainedMap; CGAL_static_assertion_msg(!(boost::is_same::NoMap>::value), @@ -552,6 +627,7 @@ mst_orient_normals( PointMap point_map = choose_param(get_param(np, internal_np::point_map), PointMap()); NormalMap normal_map = choose_param(get_param(np, internal_np::normal_map), NormalMap()); + ConstrainedMap constrained_map = choose_param(get_param(np, internal_np::point_is_constrained), ConstrainedMap()); Kernel kernel; // Bring private stuff to scope @@ -584,37 +660,46 @@ mst_orient_normals( // and get() requires a lookup in the map. IndexMap index_map(points.begin(), points.end()); - // Orients the normal of the point with maximum Z towards +Z axis. - typename PointRange::iterator source_point - = mst_find_source(points.begin(), points.end(), - point_map, normal_map, - kernel); - // Iterates over input points and creates Riemannian Graph: // - vertices are numbered like the input points index. // - vertices are empty. // - we add the edge (i, j) if either vertex i is in the k-neighborhood of vertex j, // or vertex j is in the k-neighborhood of vertex i. - Riemannian_graph riemannian_graph - = create_riemannian_graph(points.begin(), points.end(), - point_map, normal_map, index_map, - k, - kernel); + Riemannian_graph riemannian_graph; + + if (boost::is_same::NoMap>::value) + riemannian_graph = create_riemannian_graph(points.begin(), points.end(), + point_map, normal_map, index_map, + Default_constrained_map + (mst_find_source(points.begin(), points.end(), + point_map, normal_map, + kernel)), + k, + kernel); + else + riemannian_graph = create_riemannian_graph(points.begin(), points.end(), + point_map, normal_map, index_map, + constrained_map, + k, + kernel); // Creates a Minimum Spanning Tree starting at source_point MST_graph mst_graph = create_mst_graph(points.begin(), points.end(), point_map, normal_map, index_map, k, kernel, - riemannian_graph, - source_point); + riemannian_graph); memory = CGAL::Memory_sizer().virtual_size(); CGAL_TRACE(" %ld Mb allocated\n", memory>>20); CGAL_TRACE(" Calls boost::breadth_first_search()\n"); + const std::size_t num_input_points = distance(points.begin(), points.end()); + std::size_t source_point_index = num_input_points; + // Traverse the point set along the MST to propagate source_point's orientation - Propagate_normal_orientation orienter; - std::size_t source_point_index = get(index_map, source_point); + Propagate_normal_orientation orienter(source_point_index); + boost::breadth_first_search(mst_graph, vertex(source_point_index, mst_graph), // source visitor(boost::make_bfs_visitor(orienter))); diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/bbox.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/bbox.h index c3207e9acbb..0419603d52c 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/bbox.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/bbox.h @@ -76,13 +76,12 @@ namespace CGAL { GT gt = choose_param(get_param(np, internal_np::geom_traits), GT()); typename GT::Construct_bbox_3 get_bbox = gt.construct_bbox_3_object(); - typedef typename boost::graph_traits::halfedge_descriptor halfedge_descriptor; + typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; - halfedge_descriptor h0 = *(halfedges(pmesh).first); - CGAL::Bbox_3 bb = get(vpm, target(h0, pmesh)).bbox(); - BOOST_FOREACH(halfedge_descriptor h, halfedges(pmesh)) + CGAL::Bbox_3 bb; + BOOST_FOREACH(vertex_descriptor v, vertices(pmesh)) { - bb += get_bbox( get(vpm, target(h, pmesh)) ); + bb += get_bbox( get(vpm, v) ); } return bb; } diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Isotropic_remeshing/remesh_impl.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Isotropic_remeshing/remesh_impl.h index 7dd14e43a7d..33d3b231cf0 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Isotropic_remeshing/remesh_impl.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Isotropic_remeshing/remesh_impl.h @@ -675,6 +675,7 @@ namespace internal { if (is_on_border(he) || is_on_mesh(he)) { he = opposite(he, mesh_); //he now is PATCH_BORDER + e = edge(he, mesh_); CGAL_assertion(is_on_patch_border(he)); } }//end if(not on PATCH) diff --git a/Polyhedron/demo/Polyhedron/Plugins/Point_set/Point_set_normal_estimation_plugin.cpp b/Polyhedron/demo/Polyhedron/Plugins/Point_set/Point_set_normal_estimation_plugin.cpp index ab53c6c887b..5193fec725f 100644 --- a/Polyhedron/demo/Polyhedron/Plugins/Point_set/Point_set_normal_estimation_plugin.cpp +++ b/Polyhedron/demo/Polyhedron/Plugins/Point_set/Point_set_normal_estimation_plugin.cpp @@ -18,6 +18,8 @@ #include #include +#include + #include "run_with_qprogressdialog.h" #include "ui_Point_set_normal_estimation_plugin.h" @@ -77,6 +79,25 @@ struct Jet_estimate_normals_functor } }; +struct Vector_to_pmap +{ + typedef boost::readable_property_map_tag category; + typedef Point_set::Index key_type; + typedef bool value_type; + typedef value_type reference; + + std::vector* vec; + + Vector_to_pmap (std::vector* vec = NULL) : vec (vec) { } + + friend inline + reference get(const Vector_to_pmap& map, key_type p) + { + return (*map.vec)[p]; + } + +}; + using namespace CGAL::Three; class Polyhedron_demo_point_set_normal_estimation_plugin : @@ -88,6 +109,7 @@ class Polyhedron_demo_point_set_normal_estimation_plugin : Q_PLUGIN_METADATA(IID "com.geometryfactory.PolyhedronDemo.PluginInterface/1.0") QAction* actionNormalEstimation; + QAction* actionNormalOrientation; QAction* actionNormalInversion; public: @@ -99,6 +121,10 @@ public: actionNormalEstimation->setObjectName("actionNormalEstimation"); actionNormalEstimation->setProperty("subMenuName","Point Set Processing"); + actionNormalOrientation = new QAction(tr("Normal Orientation"), mainWindow); + actionNormalOrientation->setObjectName("actionNormalOrientation"); + actionNormalOrientation->setProperty("subMenuName","Point Set Processing"); + actionNormalInversion = new QAction(tr("Inverse Normal Orientations"), mainWindow); actionNormalInversion->setObjectName("actionNormalInversion"); actionNormalInversion->setProperty("subMenuName","Point Set Processing"); @@ -106,7 +132,7 @@ public: } QList actions() const { - return QList() << actionNormalEstimation << actionNormalInversion; + return QList() << actionNormalEstimation << actionNormalOrientation << actionNormalInversion; } bool applicable(QAction* action) const { @@ -124,6 +150,7 @@ public: public Q_SLOTS: void on_actionNormalEstimation_triggered(); + void on_actionNormalOrientation_triggered(); void on_actionNormalInversion_triggered(); }; // end PS_demo_smoothing_plugin @@ -142,7 +169,6 @@ class Point_set_demo_normal_estimation_dialog : public QDialog, private Ui::Norm unsigned int convolution_neighbors() const { return m_convolution_neighbors->value(); } double convolution_radius() const { return m_convolution_radius->value(); } double offset_radius() const { return m_offset_radius->value(); } - int orient_neighbors() const { return m_orient_neighbors->value(); } unsigned int method () const { @@ -151,10 +177,6 @@ class Point_set_demo_normal_estimation_dialog : public QDialog, private Ui::Norm if (buttonVCM->isChecked ()) return 2; return -1; } - bool orient () const - { - return buttonOrient->isChecked(); - } bool use_convolution_radius () const { return buttonRadius->isChecked(); @@ -209,9 +231,6 @@ void Polyhedron_demo_point_set_normal_estimation_plugin::on_actionNormalEstimati QApplication::setOverrideCursor(Qt::BusyCursor); QApplication::processEvents(); - // First point to delete - Point_set::iterator first_unoriented_point = points->end(); - //*************************************** // normal estimation //*************************************** @@ -275,59 +294,112 @@ void Polyhedron_demo_point_set_normal_estimation_plugin::on_actionNormalEstimati << (memory>>20) << " Mb allocated" << std::endl; } + item->resetMenu(); item->setRenderingMode(ShadedPoints); + // Updates scene + item->invalidateOpenGLBuffers(); + scene->itemChanged(index); + + QApplication::restoreOverrideCursor(); + } +#endif // !CGAL_DISABLE_NORMAL_ESTIMATION_PLUGIN +} + +void Polyhedron_demo_point_set_normal_estimation_plugin::on_actionNormalOrientation_triggered() +{ + const CGAL::Three::Scene_interface::Item_id index = scene->mainSelectionIndex(); + + Scene_points_with_normal_item* item = + qobject_cast(scene->item(index)); + + if(item) + { + // Gets point set + Point_set* points = item->point_set(); + if(points == NULL) + return; + + // Gets options + QMultipleInputDialog dialog ("Normal Orientation", mw); + QSpinBox* neighborhood = dialog.add ("Neighborhood Size: "); + neighborhood->setRange (1, 10000000); + neighborhood->setValue (18); + + QRadioButton* use_seed_points = NULL; + QRadioButton* orient_selection = NULL; + + if (points->nb_selected_points() != 0) + { + use_seed_points = dialog.add ("Use selection as seed points and orient the unselected points"); + use_seed_points->setChecked(true); + orient_selection = dialog.add ("Orient selection"); + orient_selection->setChecked(false); + } + + if(!dialog.exec()) + return; + + QApplication::setOverrideCursor(Qt::BusyCursor); + QApplication::processEvents(); + + // First point to delete + Point_set::iterator first_unoriented_point = points->end(); + //*************************************** // normal orientation //*************************************** - if (dialog.orient ()) - { - CGAL::Timer task_timer; task_timer.start(); - std::cerr << "Orient normals with a Minimum Spanning Tree (k=" << dialog.orient_neighbors() << ")...\n"; + CGAL::Timer task_timer; task_timer.start(); + std::cerr << "Orient normals with a Minimum Spanning Tree (k=" << neighborhood->value() << ")...\n"; - // Tries to orient normals - first_unoriented_point = - CGAL::mst_orient_normals(points->all_or_selection_if_not_empty(), - dialog.orient_neighbors(), - points->parameters()); + // Tries to orient normals + if (points->nb_selected_points() != 0 && use_seed_points->isChecked()) + { + std::vector constrained_map (points->size(), false); - std::size_t nb_unoriented_normals = std::distance(first_unoriented_point, points->end()); - std::size_t memory = CGAL::Memory_sizer().virtual_size(); - std::cerr << "Orient normals: " << nb_unoriented_normals << " point(s) with an unoriented normal are selected (" - << task_timer.time() << " seconds, " - << (memory>>20) << " Mb allocated)" - << std::endl; - - // Selects points with an unoriented normal - points->set_first_selected (first_unoriented_point); - - // Updates scene - item->invalidateOpenGLBuffers(); - scene->itemChanged(index); - - QApplication::restoreOverrideCursor(); - - // Warns user - if (nb_unoriented_normals > 0) - { - QMessageBox::information(NULL, - tr("Points with an unoriented normal"), - tr("%1 point(s) with an unoriented normal are selected.\nPlease orient them or remove them before running Poisson reconstruction.") - .arg(nb_unoriented_normals)); - } - } + for (Point_set::iterator it = points->first_selected(); it != points->end(); ++ it) + constrained_map[*it] = true; + + first_unoriented_point = + CGAL::mst_orient_normals(*points, + std::size_t(neighborhood->value()), + points->parameters(). + point_is_constrained_map(Vector_to_pmap(&constrained_map))); + } else - { - // Updates scene - item->invalidateOpenGLBuffers(); - scene->itemChanged(index); + first_unoriented_point = + CGAL::mst_orient_normals(points->all_or_selection_if_not_empty(), + std::size_t(neighborhood->value()), + points->parameters()); - QApplication::restoreOverrideCursor(); - } + std::size_t nb_unoriented_normals = std::distance(first_unoriented_point, points->end()); + std::size_t memory = CGAL::Memory_sizer().virtual_size(); + std::cerr << "Orient normals: " << nb_unoriented_normals << " point(s) with an unoriented normal are selected (" + << task_timer.time() << " seconds, " + << (memory>>20) << " Mb allocated)" + << std::endl; + + // Selects points with an unoriented normal + points->set_first_selected (first_unoriented_point); + + // Updates scene + item->invalidateOpenGLBuffers(); + scene->itemChanged(index); + + QApplication::restoreOverrideCursor(); + + // Warns user + if (nb_unoriented_normals > 0) + { + QMessageBox::information(NULL, + tr("Points with an unoriented normal"), + tr("%1 point(s) with an unoriented normal are selected.\nPlease orient them or remove them before running Poisson reconstruction.") + .arg(nb_unoriented_normals)); + } } -#endif // !CGAL_DISABLE_NORMAL_ESTIMATION_PLUGIN } + #include "Point_set_normal_estimation_plugin.moc" diff --git a/Polyhedron/demo/Polyhedron/Plugins/Point_set/Point_set_normal_estimation_plugin.ui b/Polyhedron/demo/Polyhedron/Plugins/Point_set/Point_set_normal_estimation_plugin.ui index 905ba695a9d..8ceaeaea887 100644 --- a/Polyhedron/demo/Polyhedron/Plugins/Point_set/Point_set_normal_estimation_plugin.ui +++ b/Polyhedron/demo/Polyhedron/Plugins/Point_set/Point_set_normal_estimation_plugin.ui @@ -6,8 +6,8 @@ 0 0 - 400 - 473 + 301 + 372 @@ -201,74 +201,6 @@ - - - - Qt::Horizontal - - - - - - - Qt::Vertical - - - - 20 - 40 - - - - - - - - Orient Estimated Normals - - - true - - - - - - - QFrame::StyledPanel - - - QFrame::Raised - - - - QFormLayout::AllNonFixedFieldsGrow - - - - - Neighborhood Size - - - - - - - Nearest Neighbors - - - 6 - - - 100000000 - - - 18 - - - - - - @@ -279,19 +211,6 @@ - - - - Qt::Vertical - - - - 20 - 40 - - - - @@ -344,22 +263,6 @@ - - buttonOrient - toggled(bool) - frame_4 - setEnabled(bool) - - - 199 - 316 - - - 199 - 355 - - - buttonRadius toggled(bool) diff --git a/Polyhedron/demo/Polyhedron/include/QMultipleInputDialog.h b/Polyhedron/demo/Polyhedron/include/QMultipleInputDialog.h index 17147651f29..4825f59a6c0 100644 --- a/Polyhedron/demo/Polyhedron/include/QMultipleInputDialog.h +++ b/Polyhedron/demo/Polyhedron/include/QMultipleInputDialog.h @@ -9,6 +9,7 @@ #include #include #include +#include class QMultipleInputDialog { @@ -29,8 +30,19 @@ public: template QObjectType* add (const char* name) { - QObjectType* out = new QObjectType (dialog); - form->addRow (QString(name), out); + QObjectType* out = NULL; + + if (boost::is_same::value) + { + out = dynamic_cast(new QRadioButton (QString(name), dialog)); + form->addRow (out); + } + else + { + out = new QObjectType (dialog); + form->addRow (QString(name), out); + } + return out; } diff --git a/Polyhedron_IO/include/CGAL/IO/PLY_reader.h b/Polyhedron_IO/include/CGAL/IO/PLY_reader.h index 98f15c92460..bbebda87c8b 100644 --- a/Polyhedron_IO/include/CGAL/IO/PLY_reader.h +++ b/Polyhedron_IO/include/CGAL/IO/PLY_reader.h @@ -54,7 +54,7 @@ namespace CGAL{ internal::PLY::PLY_read_number* property = element.property(k); property->get (in); - if (in.eof()) + if (in.fail()) return false; } @@ -84,7 +84,7 @@ namespace CGAL{ polygons.back()[i] = std::size_t(get<0>(new_face)[i]); } - return !in.bad(); + return true; } } @@ -106,7 +106,10 @@ namespace CGAL{ internal::PLY::PLY_reader reader; if (!(reader.init (in))) + { + in.setstate(std::ios::failbit); return false; + } for (std::size_t i = 0; i < reader.number_of_elements(); ++ i) { @@ -121,7 +124,7 @@ namespace CGAL{ internal::PLY::PLY_read_number* property = element.property(k); property->get (in); - if (in.eof()) + if (in.fail()) return false; } @@ -160,7 +163,7 @@ namespace CGAL{ internal::PLY::PLY_read_number* property = element.property(k); property->get (in); - if (in.eof()) + if (in.fail()) return false; } } @@ -189,8 +192,10 @@ namespace CGAL{ internal::PLY::PLY_reader reader; if (!(reader.init (in))) + { + in.setstate(std::ios::failbit); return false; - + } for (std::size_t i = 0; i < reader.number_of_elements(); ++ i) { internal::PLY::PLY_element& element = reader.element(i); @@ -217,7 +222,7 @@ namespace CGAL{ internal::PLY::PLY_read_number* property = element.property(k); property->get (in); - if (in.eof()) + if (in.fail()) return false; } @@ -313,14 +318,12 @@ namespace CGAL{ { internal::PLY::PLY_read_number* property = element.property(k); property->get (in); - - if (in.eof()) + if (in.fail()) return false; } } } } - return !in.bad(); } diff --git a/QP_solver/doc/QP_solver/CGAL/QP_functions.h b/QP_solver/doc/QP_solver/CGAL/QP_functions.h index 1c8b614dd2b..f06b7f6f9a4 100644 --- a/QP_solver/doc/QP_solver/CGAL/QP_functions.h +++ b/QP_solver/doc/QP_solver/CGAL/QP_functions.h @@ -4,6 +4,8 @@ namespace CGAL { /// @{ /*! +\tparam LinearProgram a model of `LinearProgram`. + This function writes a linear program to an output stream (in `MPSFormat`). The time complexity is \f$ \Theta (mn)\f$, even if the program is very sparse. @@ -12,18 +14,9 @@ It writes the linear program `lp` to `out` in `MPSFormat`. The name of the program will be the one provided by `problem_name`. -Requirements --------------- - +\cgalHeading{Requirements} Output operators are defined for all entry types of `lp`. -Example --------------- - -\ref QP_solver/print_first_lp.cpp - -\sa The concept `LinearProgram` - */ template void print_linear_program @@ -32,6 +25,8 @@ const std::string& problem_name = std::string("MY_MPS")); /*! +\tparam NonnegativeLinearProgram a model of `NonnegativeLinearProgram` + This function writes a nonnegative linear program to an output stream (in `MPSFormat`). The time complexity is \f$ \Theta (mn)\f$, even if the program is very sparse. @@ -40,18 +35,9 @@ Writes the nonnegative linear program `lp` to `out` in `MPSFormat`. The name of the program will be the one provided by `problem_name`. -Requirements --------------- - +\cgalHeading{Requirements} Output operators are defined for all entry types of `lp`. -Example --------------- - -\ref QP_solver/print_first_nonnegative_lp.cpp - -\sa The concept `NonnegativeLinearProgram` - */ template void print_nonnegative_linear_program @@ -60,6 +46,8 @@ const std::string& problem_name = std::string("MY_MPS")); /*! +\tparam NonnegativeQuadraticProgram a model of `NonnegativeQuadraticProgram` + This function writes a nonnegative quadratic program to an output stream (in `MPSFormat`). The time complexity is \f$ \Theta (n^2 + mn)\f$, even if the program is very sparse. @@ -68,18 +56,9 @@ Writes the nonnegative quadratic program `qp` to `out` in `MPSFormat`. The name of the program will be the one provided by `problem_name`. -Requirements --------------- - +\cgalHeading{Requirements} Output operators are defined for all entry types of `qp`. -Example --------------- - -\ref QP_solver/print_first_nonnegative_qp.cpp - -\sa The concept `NonnegativeQuadraticProgram` - */ template void print_nonnegative_quadratic_program @@ -88,6 +67,9 @@ const std::string& problem_name = std::string("MY_MPS")); /*! +\tparam QuadraticProgram a model of `QuadraticProgram` + + This function writes a quadratic program to an output stream (in `MPSFormat`). The time complexity is \f$ \Theta (n^2 + mn)\f$, even if the program is very sparse. @@ -95,17 +77,9 @@ if the program is very sparse. Writes the quadratic program `qp` to `out` in `MPSFormat`. The name of the program will be the one provided by `problem_name`. -Requirements --------------- - +\cgalHeading{Requirements} Output operators are defined for all entry types of `qp`. -Example --------------- - -\ref QP_solver/print_first_qp.cpp - -\sa The concept `QuadraticProgram` */ template void print_quadratic_program @@ -118,10 +92,10 @@ This function solves a linear program, using some exact Integral Domain `ET` for its computations. Various options may be provided, see `Quadratic_program_options`. -Requirements --------------- +\tparam LinearProgram a model of `LinearProgram` such as `Quadratic_program`, +`Quadratic_program_from_mps`, or `Linear_program_from_iterators` -`ET` is a model of the concepts `IntegralDomain` and +\tparam ET a model of the concepts `IntegralDomain` and `RealEmbeddable`; it must be an exact type, and all entries of `qp` are convertible to `ET`. @@ -144,15 +118,6 @@ factor. For maximum efficiency, it is advisable to define the macros \returns the solution of the linear program `lp`, solved with exact number type `ET`. -Example --------------- - -\ref QP_solver/first_lp.cpp - -\sa `Quadratic_program` -\sa `Quadratic_program_from_mps` -\sa `Linear_program_from_iterators` - */ template Quadratic_program_solution solve_linear_program @@ -165,10 +130,10 @@ This function solves a nonnegative linear program, using some exact Integral Domain `ET` for its computations. Various options may be provided, see `Quadratic_program_options`. -Requirements --------------- +\tparam NonnegativeQuadraticProgram a model of `NonnegativeQuadraticProgram` such as +`Quadratic_program`, `Quadratic_program_from_mps`, or `Nonnegative_quadratic_program_from_iterators`. -`ET` is a model of the concepts `IntegralDomain` and +\tparam ET is a model of the concepts `IntegralDomain` and `RealEmbeddable`; it must be an exact type, and all entries of `qp` are convertible to `ET`. @@ -191,16 +156,6 @@ factor. For maximum efficiency, it is advisable to define the macros \returns the solution of the nonnegative linear program `lp`, solved with exact number type `ET`. -Example --------------- - -\ref QP_solver/first_nonnegative_lp.cpp - -The models of \ref NonnegativeLinearProgram\: -\sa `Quadratic_program` -\sa `Quadratic_program_from_mps` -\sa `Nonnegative_linear_program_from_iterators` - */ template Quadratic_program_solution solve_nonnegative_linear_program @@ -213,10 +168,10 @@ This function solves a nonnegative quadratic program, using some exact Integral Domain `ET` for its computations. Various options may be provided, see `Quadratic_program_options`. -Requirements --------------- +\tparam NonnegativeQuadraticProgram a model of `NonnegativeQuadraticProgram` such as +`Quadratic_program`, `Quadratic_program_from_mps`, or `Nonnegative_quadratic_program_from_iterators`. -`ET` is a model of the concepts `IntegralDomain` and +\tparam ET a model of the concepts `IntegralDomain` and `RealEmbeddable`; it must be an exact type, and all entries of `qp` are convertible to `ET`. @@ -239,15 +194,6 @@ factor. For maximum efficiency, it is advisable to define the macros \returns the solution of the nonnegative quadratic program `qp`, solved with exact number type `ET`. -Example --------------- - -\ref QP_solver/first_nonnegative_qp.cpp - -The models of \ref ::NonnegativeQuadraticProgram\: -\sa `Quadratic_program` -\sa `Quadratic_program_from_mps` -\sa `Nonnegative_quadratic_program_from_iterators` */ template Quadratic_program_solution solve_nonnegative_quadratic_program @@ -260,10 +206,11 @@ This function solves a quadratic program, using some exact Integral Domain `ET` for its computations. Various options may be provided, see `Quadratic_program_options`. -Requirements --------------- -`ET` is a model of the concepts `IntegralDomain` and +\tparam NonnegativeQuadraticProgram a model of `NonnegativeQuadraticProgram` such as +`Quadratic_program`, `Quadratic_program_from_mps`, or `Nonnegative_quadratic_program_from_iterators`. + +\tparam ET is a model of the concepts `IntegralDomain` and `RealEmbeddable`; it must be an exact type, and all entries of `qp` are convertible to `ET`. @@ -286,16 +233,6 @@ factor. For maximum efficiency, it is advisable to define the macros \returns the solution of the quadratic program `qp`, solved with exact number type `ET`. -Example --------------- - -\ref QP_solver/first_qp.cpp - -The models of \ref QuadraticProgram\: -\sa `Quadratic_program` -\sa `Quadratic_program_from_mps` -\sa `Quadratic_program_from_iterators` - */ template Quadratic_program_solution solve_quadratic_program diff --git a/QP_solver/doc/QP_solver/CGAL/QP_models.h b/QP_solver/doc/QP_solver/CGAL/QP_models.h index 1ede1029905..28748e7f103 100644 --- a/QP_solver/doc/QP_solver/CGAL/QP_models.h +++ b/QP_solver/doc/QP_solver/CGAL/QP_models.h @@ -57,8 +57,7 @@ namespace CGAL { \cgalModels `QuadraticProgram` \cgalModels `LinearProgram` - Example - -------------- + \cgalHeading{Example} \ref QP_solver/first_lp_from_iterators.cpp @@ -113,8 +112,7 @@ public: \returns an instance of `Linear_program_from_iterators`, constructed from the given iterators. - Example - -------------- + \cgalHeading{Example} The following example demonstrates the typical usage of makers with the simpler function `make_nonnegative_linear_program_from_iterators()`. @@ -160,8 +158,7 @@ make_linear_program_from_iterators ( \returns an instance of `Nonnegative_linear_program_from_iterators`, constructed from the given iterators. - Example - -------------- + \cgalHeading{Example} `QP_solver/solve_convex_hull_containment_lp2.h` @@ -198,8 +195,7 @@ make_nonnegative_linear_program_from_iterators ( `Nonnegative_quadratic_program_from_iterators`, constructed from the given iterators. - Example - -------------- + \cgalHeading{Example} The following example demonstrates the typical usage of makers with the simpler function `make_nonnegative_linear_program_from_iterators()`. @@ -240,8 +236,7 @@ make_nonnegative_quadratic_program_from_iterators ( \returns an instance of `Quadratic_program_from_iterators`, constructed from the given iterators. - Example - -------------- + \cgalHeading{Example} The following example demonstrates the typical usage of makers with the simpler function `make_nonnegative_linear_program_from_iterators()`. @@ -330,8 +325,7 @@ make_quadratic_program_from_iterators ( \cgalModels `NonnegativeQuadraticProgram` \cgalModels `NonnegativeLinearProgram` - Example - -------------- + \cgalHeading{Example} \ref QP_solver/first_nonnegative_lp_from_iterators.cpp @@ -422,8 +416,7 @@ public: \cgalModels `QuadraticProgram` \cgalModels `NonnegativeQuadraticProgram` - Example - -------------- + \cgalHeading{Example} \ref QP_solver/first_nonnegative_qp_from_iterators.cpp @@ -522,8 +515,7 @@ public: \cgalModels `QuadraticProgram` - Example - -------------- + \cgalHeading{Example} \ref QP_solver/first_qp_from_iterators.cpp @@ -630,8 +622,7 @@ public: \cgalModels `NonnegativeQuadraticProgram` \cgalModels `NonnegativeLinearProgram` - Example - -------------- + \cgalHeading{Example} \ref QP_solver/first_qp_from_mps.cpp @@ -854,8 +845,7 @@ namespace CGAL { \cgalModels `NonnegativeQuadraticProgram` \cgalModels `NonnegativeLinearProgram` - Example - -------------- + \cgalHeading{Example} \ref QP_solver/first_qp.cpp diff --git a/QP_solver/doc/QP_solver/CGAL/QP_options.h b/QP_solver/doc/QP_solver/CGAL/QP_options.h index 658363321e9..389f4dc81c2 100644 --- a/QP_solver/doc/QP_solver/CGAL/QP_options.h +++ b/QP_solver/doc/QP_solver/CGAL/QP_options.h @@ -15,8 +15,7 @@ options referring to The idea is that this list grows in the future. -Operations --------------- +\cgalHeading{Operations} Here we just have set/get pairs for any option type. diff --git a/QP_solver/doc/QP_solver/CGAL/QP_solution.h b/QP_solver/doc/QP_solver/CGAL/QP_solution.h index c323a4fd101..4335cfac32f 100644 --- a/QP_solver/doc/QP_solver/CGAL/QP_solution.h +++ b/QP_solver/doc/QP_solver/CGAL/QP_solution.h @@ -36,13 +36,11 @@ the four functions `solve_nonnegative_quadratic_program`, and `solve_nonnegative_linear_program`. -Example --------------- +\cgalHeading{Example} \ref QP_solver/first_qp.cpp -Terminology --------------- +\cgalHeading{Terminology} If there is no \f$ \qpx\f$ that satisfies all the (in)equalities, the program is called infeasible, otherwise, it is feasible, diff --git a/Surface_mesh/include/CGAL/Surface_mesh/Surface_mesh.h b/Surface_mesh/include/CGAL/Surface_mesh/Surface_mesh.h index 18f7575e5ab..26545fc394a 100644 --- a/Surface_mesh/include/CGAL/Surface_mesh/Surface_mesh.h +++ b/Surface_mesh/include/CGAL/Surface_mesh/Surface_mesh.h @@ -1112,28 +1112,33 @@ public: bool join(const Surface_mesh& other) { + // increase capacity const size_type nv = num_vertices(), nh = num_halfedges(), nf = num_faces(); resize(num_vertices()+ other.num_vertices(), num_edges()+ other.num_edges(), num_faces()+ other.num_faces()); + // append properties in the free space created by resize vprops_.transfer(other.vprops_); hprops_.transfer(other.hprops_); fprops_.transfer(other.fprops_); eprops_.transfer(other.eprops_); + // translate halfedge index in vertex -> halfedge for(size_type i = nv; i < nv+other.num_vertices(); i++){ Vertex_index vi(i); if(vconn_[vi].halfedge_ != null_halfedge()){ vconn_[vi].halfedge_ = Halfedge_index(size_type(vconn_[vi].halfedge_)+nh); } } + // translate halfedge index in face -> halfedge for(size_type i = nf; i < nf+other.num_faces(); i++){ Face_index fi(i); if(fconn_[fi].halfedge_ != null_halfedge()){ fconn_[fi].halfedge_ = Halfedge_index(size_type(fconn_[fi].halfedge_)+nh); } } + // translate indices in halfedge -> face, halfedge -> target, halfedge -> prev, and halfedge -> next for(size_type i = nh; i < nh+other.num_halfedges(); i++){ Halfedge_index hi(i); if(hconn_[hi].face_ != null_face()){ @@ -1150,43 +1155,50 @@ public: } } size_type inf_value = (std::numeric_limits::max)(); + + // merge vertex free list if(other.vertices_freelist_ != inf_value){ - if(vertices_freelist_ != inf_value){ - Vertex_index vi(nv+other.vertices_freelist_); - Halfedge_index inf((std::numeric_limits::max)()); - while(vconn_[vi].halfedge_ != inf){ - Vertex_index corrected_vi = Vertex_index(size_type(vconn_[vi].halfedge_)+nv-nh); - vconn_[vi].halfedge_ = Halfedge_index(corrected_vi); - vi = corrected_vi; - } - vconn_[vi].halfedge_ = Halfedge_index(vertices_freelist_); + Vertex_index vi(nv+other.vertices_freelist_); + Halfedge_index inf((std::numeric_limits::max)()); + // correct the indices in the linked list of free vertices copied (due to vconn_ translation) + while(vconn_[vi].halfedge_ != inf){ + Vertex_index corrected_vi = Vertex_index(size_type(vconn_[vi].halfedge_)+nv-nh); + vconn_[vi].halfedge_ = Halfedge_index(corrected_vi); + vi = corrected_vi; } + // append the vertex free linked list of `this` to the copy of `other` + vconn_[vi].halfedge_ = Halfedge_index(vertices_freelist_); + // update the begin of the vertex free linked list vertices_freelist_ = nv + other.vertices_freelist_; } + // merge face free list if(other.faces_freelist_ != inf_value){ - if(faces_freelist_ != inf_value){ - Face_index fi(nf+other.faces_freelist_); - Halfedge_index inf((std::numeric_limits::max)()); - while(fconn_[fi].halfedge_ != inf){ - Face_index corrected_fi = Face_index(size_type(fconn_[fi].halfedge_)+nf-nh); - fconn_[fi].halfedge_ = Halfedge_index(corrected_fi); - fi = corrected_fi; - } - fconn_[fi].halfedge_ = Halfedge_index(faces_freelist_); + Face_index fi(nf+other.faces_freelist_); + Halfedge_index inf((std::numeric_limits::max)()); + // correct the indices in the linked list of free faces copied (due to fconn_ translation) + while(fconn_[fi].halfedge_ != inf){ + Face_index corrected_fi = Face_index(size_type(fconn_[fi].halfedge_)+nf-nh); + fconn_[fi].halfedge_ = Halfedge_index(corrected_fi); + fi = corrected_fi; } + // append the face free linked list of `this` to the copy of `other` + fconn_[fi].halfedge_ = Halfedge_index(faces_freelist_); + // update the begin of the face free linked list faces_freelist_ = nf + other.faces_freelist_; } + // merge edge free list if(other.edges_freelist_ != inf_value){ - if(edges_freelist_ != inf_value){ - Halfedge_index hi(nh+other.edges_freelist_); - Halfedge_index inf((std::numeric_limits::max)()); - while(hconn_[hi].next_halfedge_ != inf){ - hi = hconn_[hi].next_halfedge_; - } - hconn_[hi].next_halfedge_ = Halfedge_index(edges_freelist_); + Halfedge_index hi(nh+other.edges_freelist_); + Halfedge_index inf((std::numeric_limits::max)()); + while(hconn_[hi].next_halfedge_ != inf){ + hi = hconn_[hi].next_halfedge_; } + // append the halfedge free linked list of `this` to the copy of `other` + hconn_[hi].next_halfedge_ = Halfedge_index(edges_freelist_); + // update the begin of the halfedge free linked list edges_freelist_ = nh + other.edges_freelist_; } + // update garbage infos garbage_ = garbage_ || other.garbage_; removed_vertices_ += other.removed_vertices_; removed_edges_ += other.removed_edges_; @@ -1427,6 +1439,33 @@ public: if(!valid && verbose){ std::cerr << "#faces: iterated: " << fcount << " vs number_of_faces(): " << number_of_faces()<< std::endl; } + + size_type inf = (std::numeric_limits::max)(); + size_type vfl = vertices_freelist_; + size_type rv = 0; + while(vfl != inf){ + vfl = (size_type)vconn_[Vertex_index(vfl)].halfedge_; + rv++; + } + valid = valid && ( rv == removed_vertices_ ); + + + size_type efl = edges_freelist_; + size_type re = 0; + while(efl != inf){ + efl = (size_type)hconn_[Halfedge_index(efl)].next_halfedge_; + re++; + } + valid = valid && ( re == removed_edges_ ); + + size_type ffl = faces_freelist_; + size_type rf = 0; + while(ffl != inf){ + ffl = (size_type)fconn_[Face_index(ffl)].halfedge_; + rf++; + } + valid = valid && ( rf == removed_faces_ ); + return valid; } diff --git a/Surface_mesh_parameterization/test/Surface_mesh_parameterization/CMakeLists.txt b/Surface_mesh_parameterization/test/Surface_mesh_parameterization/CMakeLists.txt index 9455b197844..e325bf4aafa 100644 --- a/Surface_mesh_parameterization/test/Surface_mesh_parameterization/CMakeLists.txt +++ b/Surface_mesh_parameterization/test/Surface_mesh_parameterization/CMakeLists.txt @@ -14,24 +14,7 @@ find_package(CGAL QUIET) if ( CGAL_FOUND ) - # VisualC++ optimization for applications dealing with large data - if (MSVC) - # Allow Windows applications to use up to 3GB of RAM - SET (CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} /LARGEADDRESSAWARE") - - # Turn off stupid VC++ warnings - SET (CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} /wd4267 /wd4311 /wd4800 /wd4503 /wd4244 /wd4345 /wd4996 /wd4396 /wd4018") - - # Print new compilation options - message( STATUS "USING DEBUG CXXFLAGS = '${CMAKE_CXX_FLAGS} ${CMAKE_CXX_FLAGS_DEBUG}'" ) - message( STATUS "USING DEBUG EXEFLAGS = '${CMAKE_EXE_LINKER_FLAGS} ${CMAKE_EXE_LINKER_FLAGS_DEBUG}'" ) - message( STATUS "USING RELEASE CXXFLAGS = '${CMAKE_CXX_FLAGS} ${CMAKE_CXX_FLAGS_RELEASE}'" ) - message( STATUS "USING RELEASE EXEFLAGS = '${CMAKE_EXE_LINKER_FLAGS} ${CMAKE_EXE_LINKER_FLAGS_RELEASE}'" ) - endif(MSVC) - - - #find Eigen find_package(Eigen3 3.1.0) #(requires 3.1.0 or greater) if(EIGEN3_FOUND) include( ${EIGEN3_USE_FILE} ) diff --git a/TDS_3/include/CGAL/Triangulation_data_structure_3.h b/TDS_3/include/CGAL/Triangulation_data_structure_3.h index f248dd783f3..ff1fce1d78e 100644 --- a/TDS_3/include/CGAL/Triangulation_data_structure_3.h +++ b/TDS_3/include/CGAL/Triangulation_data_structure_3.h @@ -78,9 +78,13 @@ class Triangulation_data_structure_3 typedef Triangulation_data_structure_3 Tds; public: - typedef Concurrency_tag_ Concurrency_tag; + // This tag is used in the parallel operations of RT_3 to access some functions + // of the TDS (tds.vertices().is_used(Vertex_handle)) that are much more efficient + // than what is exposed by the TDS concept (tds.is_vertex(Vertex_handle)). + typedef CGAL::Tag_true Is_CGAL_TDS_3; + // Tools to change the Vertex and Cell types of the TDS. template < typename Vb2 > struct Rebind_vertex { diff --git a/Triangulation_3/include/CGAL/Regular_triangulation_3.h b/Triangulation_3/include/CGAL/Regular_triangulation_3.h index 5170956468e..22095049387 100644 --- a/Triangulation_3/include/CGAL/Regular_triangulation_3.h +++ b/Triangulation_3/include/CGAL/Regular_triangulation_3.h @@ -776,6 +776,38 @@ public: return std::copy(vertices.begin(), vertices.end(), res); } + // In parallel operations, we need to be able to check the health of the 'hint' vertex handle, + // which might be invalided by other threads. One way to do that is the 'is_vertex()' function + // of the TDS, but it runs in O(sqrt(n)) complexity. When we are using our TDS, we can use + // a lower level function from the compact container, which runs in constant time. + BOOST_MPL_HAS_XXX_TRAIT_DEF(Is_CGAL_TDS_3) + + template ::value> + struct Is_CGAL_TDS_3 : public CGAL::Tag_false + { }; + + template + struct Is_CGAL_TDS_3 : public CGAL::Boolean_tag + { }; + + template ::value> + struct Vertex_validity_checker + { + bool operator()(const typename TDS_::Vertex_handle vh_, const TDS_& tds_) { + return tds_.is_vertex(vh_); + } + }; + + template + struct Vertex_validity_checker + { + bool operator()(const typename TDS_::Vertex_handle vh_, const TDS_& tds_) { + return tds_.vertices().is_used(vh_); + } + }; + void remove(Vertex_handle v); // Concurrency-safe // See Triangulation_3::remove for more information @@ -1359,14 +1391,40 @@ protected: #endif Vertex_handle& hint = m_tls_hint.local(); + Vertex_validity_checker vertex_validity_check; + for(size_t i_point = r.begin() ; i_point != r.end() ; ++i_point) { bool success = false; const Weighted_point& p = m_points[i_point]; while(!success) { - if(m_rt.try_lock_vertex(hint) && m_rt.try_lock_point(p)) + // The 'hint' is unsafe to use immediately because we are in a regular triangulation, + // and the insertion of a (weighted) point in another thread might have hidden (deleted) + // the hint. + if(!vertex_validity_check(hint, m_rt.tds())) { + hint = m_rt.finite_vertices_begin(); + continue; + } + + // We need to make sure that while are locking the position P1 := hint->point(), 'hint' + // does not get its position changed to P2 != P1. + const Weighted_point hint_point_mem = hint->point(); + + if(m_rt.try_lock_point(hint_point_mem) && m_rt.try_lock_point(p)) + { + // Make sure that the hint is still valid (so that we can safely take hint->cell()) and + // that its position hasn't changed to ensure that we will start the locate from where + // we have locked. + if(!vertex_validity_check(hint, m_rt.tds()) || + hint->point() != hint_point_mem) + { + hint = m_rt.finite_vertices_begin(); + m_rt.unlock_all_elements(); + continue; + } + bool could_lock_zone; Locate_type lt; int li, lj; @@ -1445,6 +1503,8 @@ protected: #endif Vertex_handle& hint = m_tls_hint.local(); + Vertex_validity_checker vertex_validity_check; + for(size_t i_idx = r.begin() ; i_idx != r.end() ; ++i_idx) { bool success = false; @@ -1452,14 +1512,37 @@ protected: const Weighted_point& p = m_points[i_point]; while(!success) { - if(m_rt.try_lock_vertex(hint) && m_rt.try_lock_point(p)) + // The 'hint' is unsafe to use immediately because we are in a regular triangulation, + // and the insertion of a (weighted) point in another thread might have hidden (deleted) + // the hint. + if(!vertex_validity_check(hint, m_rt.tds())) { + hint = m_rt.finite_vertices_begin(); + continue; + } + + // We need to make sure that while are locking the position P1 := hint->point(), 'hint' + // does not get its position changed to P2 != P1. + const Weighted_point hint_point_mem = hint->point(); + + if(m_rt.try_lock_point(hint_point_mem) && m_rt.try_lock_point(p)) + { + // Make sure that the hint is still valid (so that we can safely take hint->cell()) and + // that its position hasn't changed to ensure that we will start the locate from where + // we have locked. + if(!vertex_validity_check(hint, m_rt.tds()) || + hint->point() != hint_point_mem) + { + hint = m_rt.finite_vertices_begin(); + m_rt.unlock_all_elements(); + continue; + } + bool could_lock_zone; Locate_type lt; int li, lj; - Cell_handle c = m_rt.locate(p, lt, li, lj, hint->cell(), - &could_lock_zone); + Cell_handle c = m_rt.locate(p, lt, li, lj, hint->cell(), &could_lock_zone); Vertex_handle v; if(could_lock_zone) v = m_rt.insert(p, lt, c, li, lj, &could_lock_zone); @@ -2461,8 +2544,13 @@ remove(Vertex_handle v, bool *could_lock_zone) } else { - Vertex_handle hint = (v->cell()->vertex(0) == v ? - v->cell()->vertex(1) : v->cell()->vertex(0)); + Vertex_validity_checker vertex_validity_check; + + // Check that the vertex hasn't be deleted from the TDS while we were locking it + if(!vertex_validity_check(v, tds())) + return true; // vertex is already gone from the TDS, nothing to do + + Vertex_handle hint = v->cell()->vertex(0) == v ? v->cell()->vertex(1) : v->cell()->vertex(0); Self tmp; Vertex_remover remover(tmp); @@ -2470,21 +2558,62 @@ remove(Vertex_handle v, bool *could_lock_zone) if(*could_lock_zone && removed) { - // Re-insert the points that v was hiding. - for(typename Vertex_remover::Hidden_points_iterator - hi = remover.hidden_points_begin(); - hi != remover.hidden_points_end(); ++hi) - { - bool could_lock_zone = false; - Vertex_handle hv; - while(!could_lock_zone) - { - hv = insert(*hi, hint, &could_lock_zone); - } + // The vertex has been removed, re-insert the points that 'v' was hiding - if(hv != Vertex_handle()) - hint = hv; + // Start by unlocking the area of the removed vertex to avoid deadlocks + this->unlock_all_elements(); + + for(typename Vertex_remover::Hidden_points_iterator + hi = remover.hidden_points_begin(); + hi != remover.hidden_points_end(); ++hi) + { + const Weighted_point& wp = *hi; + + // try to lock the positions of the hint and the hidden point + bool success = false; + while(!success) + { + // The 'hint' is unsafe to use immediately because we are in a regular triangulation, + // and the insertion of a (weighted) point in another thread might have hidden (deleted) + // the hint. + if(!vertex_validity_check(hint, tds())) + { + hint = finite_vertices_begin(); + continue; + } + + // We need to make sure that while are locking the position P1 := hint->point(), 'hint' + // does not get its position changed to P2 != P1. + const Weighted_point hint_point_mem = hint->point(); + + if(this->try_lock_point(hint_point_mem) && this->try_lock_point(wp)) + { + // Make sure that the hint is still valid (so that we can safely take hint->cell()) and + // that its position hasn't changed to ensure that we will start the locate from where + // we have locked. + if(!vertex_validity_check(hint, tds()) || + hint->point() != hint_point_mem) + { + hint = finite_vertices_begin(); + this->unlock_all_elements(); + continue; + } + + Vertex_handle hv = insert(wp, hint, could_lock_zone); + + if(*could_lock_zone) + { + success = true; + if(hv != Vertex_handle()) + hint = hv; + } + } + + // This unlocks everything in all cases: partial lock failure, failed insertion, successful insertion + this->unlock_all_elements(); + } } + CGAL_triangulation_expensive_postcondition(is_valid()); } } diff --git a/Triangulation_3/test/Triangulation_3/include/CGAL/_test_cls_parallel_triangulation_3.h b/Triangulation_3/test/Triangulation_3/include/CGAL/_test_cls_parallel_triangulation_3.h index 59abb7b6da2..facbd70f2d7 100644 --- a/Triangulation_3/test/Triangulation_3/include/CGAL/_test_cls_parallel_triangulation_3.h +++ b/Triangulation_3/test/Triangulation_3/include/CGAL/_test_cls_parallel_triangulation_3.h @@ -27,42 +27,82 @@ #include +template +struct PTR_random_pts_generator +{ + typedef typename Parallel_triangulation::Point Point; + + void operator()(const int num, CGAL::Random& rnd, std::vector& points) const + { + CGAL::Random_points_in_cube_3 gen(1., rnd); + + points.reserve(num); + for(int i=0; i!=num; ++i) + points.push_back(*gen++); + } +}; + +template +struct PTR_random_pts_generator +{ + typedef typename Parallel_triangulation::Bare_point Bare_point; + typedef typename Parallel_triangulation::Weighted_point Weighted_point; + + void operator()(const int num, CGAL::Random& rnd, std::vector& points) const + { + CGAL::Random_points_in_cube_3 gen(1., rnd); + + points.reserve(num); + for(int i=0; i!=num; ++i) + points.push_back(Weighted_point(*gen++, rnd.get_double(-1., 1.))); + } +}; + template void _test_cls_parallel_triangulation_3(const Parallel_triangulation &) { - const int NUM_INSERTED_POINTS = 5000; - typedef Parallel_triangulation Cls; typedef typename Cls::Vertex_handle Vertex_handle; typedef typename Cls::Point Point; - CGAL::Random_points_in_cube_3 rnd(1.); + typedef std::vector Point_container; + + CGAL::Random rnd; + std::cout << "Seed: " << rnd.get_seed() << std::endl; + + const int num_insert = 5000; + Point_container points; + PTR_random_pts_generator points_gen; + points_gen(num_insert, rnd, points); - // Construction from a vector of points - std::vector points; - points.reserve(NUM_INSERTED_POINTS); - for (int i = 0; i != NUM_INSERTED_POINTS; ++i) - points.push_back(*rnd++); - // Construct the locking data-structure, using the bounding-box of the points - typename Cls::Lock_data_structure locking_ds( - CGAL::Bbox_3(-1., -1., -1., 1., 1., 1.), 50); + typename Cls::Lock_data_structure locking_ds(CGAL::Bbox_3(-1., -1., -1., 1., 1., 1.), 50); + // Contruct the triangulation in parallel std::cout << "Construction and parallel insertion" << std::endl; Cls tr(points.begin(), points.end(), &locking_ds); + std::cout << "Triangulation has " << tr.number_of_vertices() << " vertices" << std::endl; + assert(tr.is_valid()); std::cout << "Parallel removal" << std::endl; - // Remove the first 100,000 vertices std::vector vertices_to_remove; typename Cls::Finite_vertices_iterator vit = tr.finite_vertices_begin(); - for (int i = 0 ; i < NUM_INSERTED_POINTS/10 ; ++i) + + const std::size_t num_remove = tr.number_of_vertices() / 10; + std::cout << "Removing " << num_remove << " from " << tr.number_of_vertices() << " vertices" << std::endl; + + for(std::size_t i=0 ; i Tds_parallel; typedef CGAL::Regular_triangulation_3< traits, Tds_parallel, Lock_ds> RT_parallel; - // The following test won't do things in parallel since it doesn't provide - // a lock data structure + + // The following test won't do things in parallel since it doesn't provide a lock data structure test_RT(); + // This test performs parallel operations _test_cls_parallel_triangulation_3( RT_parallel() ); #endif