PMP: Accelerate isotropic remeshing

BGL/include/CGAL/boost/graph/Euler_operations.h

@@ -737,7 +737,7 @@ add_face(const VertexRange& vr, Graph& g)
                       patch_start, patch_end;
   // cache for set_next and vertex' set_halfedge
   typedef std::pair<halfedge_descriptor, halfedge_descriptor> NextCacheEntry;
-  typedef std::vector<NextCacheEntry>    NextCache;
+  typedef boost::container::small_vector<NextCacheEntry,9>    NextCache;
   NextCache next_cache;
   next_cache.reserve(3 * n);
 

Intersections_3/include/CGAL/Intersections_3/internal/Bbox_3_Sphere_3_do_intersect.h

@@ -35,39 +35,59 @@ namespace internal {
         typedef typename K::Point_3 Point;
         FT d = FT(0);
         FT distance = FT(0);
+        FT sr = sphere.squared_radius();
+
                 Point center = sphere.center();
 
                 if(center.x() < (FT)bbox.xmin())
                 {
                         d = (FT)bbox.xmin() - center.x();
-                        distance += d * d;
+                        d = square(d);
+                        if(d > sr){
+                          return false;
+                        }
+                        distance = d;
                 }
                 else if(center.x() > (FT)bbox.xmax())
                 {
                         d = center.x() - (FT)bbox.xmax();
-                        distance += d * d;
+                        d = square(d);
+                        if(d > sr){
+                          return false;
+                        }
+                        distance = d;
                 }
 
                 if(center.y() < (FT)bbox.ymin())
                 {
                         d = (FT)bbox.ymin() - center.y();
-                        distance += d * d;
+                        d = square(d);
+                        if(d > sr){
+                          return false;
+                        }
+                        distance += d ;
                 }
                 else if(center.y() > (FT)bbox.ymax())
                 {
                         d = center.y() - (FT)bbox.ymax();
-                        distance += d * d;
+                        d = square(d);
+                        if(d > sr){
+                          return false;
+                        }
+                        distance += d;
                 }
 
                 if(center.z() < (FT)bbox.zmin())
                 {
                         d = (FT)bbox.zmin() - center.z();
-                        distance += d * d;
+                        d = square(d);
+                        distance += d;
                 }
                 else if(center.z() > (FT)bbox.zmax())
                 {
                         d = center.z() - (FT)bbox.zmax();
-                        distance += d * d;
+                        d = square(d);
+                        distance += d;
                 }
 
                 // For unknown reason this causes a syntax error on VC2005
@@ -87,7 +107,7 @@ namespace internal {
                 //        }
                 //}
 
-                return distance <= sphere.squared_radius();
+                return distance <= sr;
     }
 
     template <class K>

Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/compute_normal.h

@@ -30,6 +30,7 @@
 #include <limits>
 #include <utility>
 #include <vector>
+#include <unordered_map>
 
 #ifdef CGAL_PMP_COMPUTE_NORMAL_DEBUG_PP
 # ifndef CGAL_PMP_COMPUTE_NORMAL_DEBUG
@@ -506,6 +507,7 @@ compute_vertex_normal_most_visible_min_circle(typename boost::graph_traits<Polyg
   typedef typename GT::Vector_3                                            Vector_3;
 
   std::vector<face_descriptor> incident_faces;
+  incident_faces.reserve(8);
   for(face_descriptor f : CGAL::faces_around_target(halfedge(v, pmesh), pmesh))
   {
     if(f == boost::graph_traits<PolygonMesh>::null_face())
@@ -669,7 +671,7 @@ compute_vertex_normal(typename boost::graph_traits<PolygonMesh>::vertex_descript
   VPMap vpmap = choose_parameter(get_parameter(np, internal_np::vertex_point),
                                  get_const_property_map(vertex_point, pmesh));
 
-  typedef std::map<face_descriptor, Vector_3>                                 Face_vector_map;
+  typedef std::unordered_map<face_descriptor, Vector_3>                       Face_vector_map;
   typedef boost::associative_property_map<Face_vector_map>                    Default_map;
 
   typedef typename internal_np::Lookup_named_param_def<internal_np::face_normal_t,

Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Isotropic_remeshing/remesh_impl.h

@@ -43,6 +43,7 @@
 #include <boost/unordered_set.hpp>
 #include <boost/shared_ptr.hpp>
 #include <boost/container/flat_set.hpp>
+#include <boost/optional.hpp>
 
 #include <map>
 #include <list>
@@ -1393,6 +1394,7 @@ namespace internal {
 
     bool collapse_would_invert_face(const halfedge_descriptor& h) const
     {
+      vertex_descriptor tv = target(h, mesh_);
       typename boost::property_traits<VertexPointMap>::reference
         s = get(vpmap_, source(h, mesh_)); //s for source
       typename boost::property_traits<VertexPointMap>::reference
@@ -1409,16 +1411,21 @@ namespace internal {
         if (face(hd, mesh_) == boost::graph_traits<PM>::null_face())
           continue;
 
+        vertex_descriptor tnhd = target(next(hd, mesh_), mesh_);
+        vertex_descriptor tnnhd = target(next(next(hd, mesh_), mesh_), mesh_);
         typename boost::property_traits<VertexPointMap>::reference
-          p = get(vpmap_, target(next(hd, mesh_), mesh_));
+          p = get(vpmap_, tnhd);
         typename boost::property_traits<VertexPointMap>::reference
-          q = get(vpmap_, target(next(next(hd, mesh_), mesh_), mesh_));
+          q = get(vpmap_, tnnhd);
 
 #ifdef CGAL_PMP_REMESHING_DEBUG
         CGAL_assertion((Triangle_3(t, p, q).is_degenerate())
                      == GeomTraits().collinear_3_object()(t, p, q));
 #endif
 
+        if((tv == tnnhd) || (tv == tnhd))
+          continue;
+
         if ( GeomTraits().collinear_3_object()(s, p, q)
           || GeomTraits().collinear_3_object()(t, p, q))
           continue;
@@ -1879,8 +1886,8 @@ namespace internal {
     bool check_normals(const HalfedgeRange& hedges) const
     {
       std::size_t nb_patches = patch_id_to_index_map.size();
+      std::vector<boost::optional<Vector_3> > normal_per_patch(nb_patches,boost::none);
 
-      std::vector< std::vector<Vector_3> > normals_per_patch(nb_patches);
       for(halfedge_descriptor hd : hedges)
       {
         Halfedge_status s = status(hd);
@@ -1892,18 +1899,15 @@ namespace internal {
         if (n == CGAL::NULL_VECTOR) //for degenerate faces
           continue;
         Patch_id pid = get_patch_id(face(hd, mesh_));
-        normals_per_patch[patch_id_to_index_map.at(pid)].push_back(n);
-      }
-
-      //on each surface patch,
-      //check all normals have same orientation
-      for (std::size_t i=0; i < nb_patches; ++i)
-      {
-        const std::vector<Vector_3>& normals = normals_per_patch[i];
-        if (normals.empty()) continue;
-
-        if (!check_orientation(normals))
-          return false;
+        std::size_t index = patch_id_to_index_map.at(pid);
+        if(normal_per_patch[index]){
+          const Vector_3& vec = *normal_per_patch[index];
+          double dot = to_double(n * vec);
+          if (dot <= 0.){
+            return false;
+          }
+        }
+        normal_per_patch[index] = boost::make_optional(n);
       }
       return true;
     }
@@ -1916,20 +1920,7 @@ namespace internal {
       Vector_3 no = compute_normal(face(opposite(h, mesh_), mesh_));
       return n * no > 0.;
     }
-
-    bool check_orientation(const std::vector<Vector_3>& normals) const
-    {
-      if (normals.size() < 2)
-        return true;
-      for (std::size_t i = 1; i < normals.size(); ++i)/*start at 1 on purpose*/
-      {
-        double dot = to_double(normals[i - 1] * normals[i]);
-        if (dot <= 0.)
-          return false;
-      }
-      return true;
-    }
-
+
   public:
     const Triangle_list& input_triangles() const {
       return input_triangles_;

Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/polygon_soup_to_polygon_mesh.h

@@ -28,6 +28,7 @@
 #include <boost/range/size.hpp>
 #include <boost/range/value_type.hpp>
 #include <boost/range/reference.hpp>
+#include <boost/container/flat_set.hpp>
 
 #include <array>
 #include <set>
@@ -177,14 +178,16 @@ bool is_polygon_soup_a_polygon_mesh(const PolygonRange& polygons)
   //check there is no duplicated ordered edge, and
   //check there is no polygon with twice the same vertex
   std::set<std::pair<V_ID, V_ID> > edge_set;
+  boost::container::flat_set<V_ID> polygon_vertices;
   V_ID max_id = 0;
+
   for(const Polygon& polygon : polygons)
   {
     std::size_t nb_edges = boost::size(polygon);
     if(nb_edges < 3)
       return false;
 
-    std::set<V_ID> polygon_vertices;
+    polygon_vertices.clear();
     V_ID prev = *std::prev(boost::end(polygon));
     for(V_ID id : polygon)
     {