Merge pull request #4250 from maxGimeno/PMP-Side_of_tm-use_bbox-maxGimeno

PMP: Make Side_of_triangle_mesh faster for disjoint meshes

AABB_tree/include/CGAL/AABB_face_graph_triangle_primitive.h

@@ -38,7 +38,7 @@ namespace CGAL {
  *\tparam FaceGraph is a model of the face graph concept.
  *\tparam VertexPointPMap  is a property map with `boost::graph_traits<FaceGraph>::%vertex_descriptor`
  *   as key type and a \cgal Kernel `Point_3` as value type.
- *                         The default is `typename boost::property_map< FaceGraph,vertex_point_t>::%type`.
+ *                         The default is `typename boost::property_map< FaceGraph,vertex_point_t>::%const_type`.
  *\tparam OneFaceGraphPerTree is either `CGAL::Tag_true` or `CGAL::Tag_false`.
  * In the former case, we guarantee that all the primitives will be from a
  * common `FaceGraph` and some data will be factorized so that the size of

AABB_tree/include/CGAL/AABB_halfedge_graph_segment_primitive.h

@@ -49,7 +49,7 @@ namespace CGAL {
  * \tparam HalfedgeGraph is a model of the halfedge graph concept.
  *   as key type and a \cgal Kernel `Point_3` as value type. 
  * \tparam VertexPointPMap is a property map with `boost::graph_traits<HalfedgeGraph>::%vertex_descriptor`.
- *                         The default is `typename boost::property_map< HalfedgeGraph,vertex_point_t>::%type`.
+ *                         The default is `typename boost::property_map< HalfedgeGraph,vertex_point_t>::%const_type`.
  * \tparam OneHalfedgeGraphPerTree is either `CGAL::Tag_true` or `CGAL::Tag_false`.
  * In the former case, we guarantee that all the primitives will be from a
  * common `HalfedgeGraph` and some data will be factorized so that the size of
@@ -77,17 +77,17 @@ class AABB_halfedge_graph_segment_primitive
                               HalfedgeGraph,
                               typename Default::Get<VertexPointPMap,
                                                     typename boost::property_map< HalfedgeGraph,
-                                                                                  vertex_point_t>::type >::type >,
+                                                                                  vertex_point_t>::const_type >::type >,
                             Source_point_from_edge_descriptor_map<
                               HalfedgeGraph,
                               typename Default::Get<VertexPointPMap,
                                                     typename boost::property_map< HalfedgeGraph,
-                                                                                  vertex_point_t>::type >::type >,
+                                                                                  vertex_point_t>::const_type >::type >,
                             OneHalfedgeGraphPerTree,
                             CacheDatum >
 #endif
 {
-  typedef typename Default::Get<VertexPointPMap,typename boost::property_map< HalfedgeGraph,vertex_point_t>::type >::type  VertexPointPMap_;
+  typedef typename Default::Get<VertexPointPMap,typename boost::property_map< HalfedgeGraph,vertex_point_t>::const_type >::type  VertexPointPMap_;
   typedef typename boost::graph_traits<HalfedgeGraph>::edge_descriptor ED;
   typedef typename boost::mpl::if_<OneHalfedgeGraphPerTree, ED, std::pair<ED, const HalfedgeGraph*> >::type Id_;
 
@@ -118,7 +118,7 @@ public:
   /*!
   The point type.
   */
-  typedef boost::property_traits< boost::property_map< HalfedgeGraph, vertex_point_t>::type >::value_type Point;
+  typedef boost::property_traits< boost::property_map< HalfedgeGraph, vertex_point_t>::const_type >::value_type Point;
   /*!
   Geometric data type.
   */

Mesh_3/include/CGAL/Polyhedral_mesh_domain_3.h

@@ -206,7 +206,7 @@ public:
   struct Primitive_type {
       //setting OneFaceGraphPerTree to false transforms the id type into 
       //std::pair<FD, const FaceGraph*>.
-    typedef AABB_face_graph_triangle_primitive<P, typename boost::property_map<P,vertex_point_t>::type, CGAL::Tag_false> type;
+    typedef AABB_face_graph_triangle_primitive<P, typename boost::property_map<P,vertex_point_t>::const_type, CGAL::Tag_false> type;
 
     static
     typename IGT_::Triangle_3 datum(const typename type::Id primitive_id) {

Polygon_mesh_processing/include/CGAL/Side_of_triangle_mesh.h

@@ -19,6 +19,7 @@
 #include <CGAL/disable_warnings.h>
 
 #include <CGAL/Polygon_mesh_processing/internal/Side_of_triangle_mesh/Point_inside_vertical_ray_cast.h>
+#include <CGAL/Polygon_mesh_processing/bbox.h>
 #include <CGAL/AABB_face_graph_triangle_primitive.h>
 
 #include <CGAL/AABB_tree.h>
@@ -66,7 +67,7 @@ namespace CGAL {
  */
 template <class TriangleMesh,
           class GeomTraits,
-          class VertexPointMap = Default
+          class VertexPointMap_ = Default
 #ifndef DOXYGEN_RUNNING
           , class AABBTree = Default
 #endif
@@ -76,25 +77,35 @@ class Side_of_triangle_mesh
   // typedefs
   template <typename TriangleMesh_,
             typename GeomTraits_,
-            typename VertexPointMap_>
+            typename VertexPointMap__>
   struct AABB_tree_default {
     typedef CGAL::AABB_face_graph_triangle_primitive<TriangleMesh_,
-                                                     VertexPointMap_> Primitive;
+                                                     VertexPointMap__> Primitive;
     typedef CGAL::AABB_traits<GeomTraits_, Primitive> Traits;
     typedef CGAL::AABB_tree<Traits> type;
   };
   typedef typename Default::Lazy_get<AABBTree,
                                      AABB_tree_default<TriangleMesh,
                                                        GeomTraits,
-                                                       VertexPointMap>
+                                                       VertexPointMap_>
                                      >::type AABB_tree_;
+  typedef typename Default::Get<VertexPointMap_,
+                                typename boost::property_map<TriangleMesh,
+                                                             vertex_point_t>::const_type>::type
+                                  VertexPointMap;
   typedef typename GeomTraits::Point_3 Point;
 
   //members
   typename GeomTraits::Construct_ray_3     ray_functor;
   typename GeomTraits::Construct_vector_3  vector_functor;
-  const AABB_tree_* tree_ptr;
+  mutable const AABB_tree_* tree_ptr;
+  const TriangleMesh* tm_ptr;
+  boost::optional<VertexPointMap> opt_vpm;
   bool own_tree;
+  CGAL::Bbox_3 box;
+#ifdef CGAL_HAS_THREADS
+  mutable CGAL_MUTEX tree_mutex;
+#endif
 
 public:
 
@@ -118,14 +129,14 @@ public:
                         const GeomTraits& gt=GeomTraits())
   : ray_functor(gt.construct_ray_3_object())
   , vector_functor(gt.construct_vector_3_object())
+  , tree_ptr(nullptr)
+  , tm_ptr(&tmesh)
+  , opt_vpm(vpmap)
   , own_tree(true)
   {
     CGAL_assertion(CGAL::is_triangle_mesh(tmesh));
     CGAL_assertion(CGAL::is_closed(tmesh));
-
+    box = Polygon_mesh_processing::bbox(tmesh, parameters::vertex_point_map(vpmap));
-    tree_ptr = new AABB_tree(faces(tmesh).first,
-                             faces(tmesh).second,
-                             tmesh, vpmap);
   }
 
   /**
@@ -138,17 +149,8 @@ public:
   */
   Side_of_triangle_mesh(const TriangleMesh& tmesh,
                         const GeomTraits& gt=GeomTraits())
-  : ray_functor(gt.construct_ray_3_object())
+  : Side_of_triangle_mesh(tmesh, get(vertex_point, tmesh), gt)
-  , vector_functor(gt.construct_vector_3_object())
+  {}
-  , own_tree(true)
-  {
-    CGAL_assertion(CGAL::is_triangle_mesh(tmesh));
-    CGAL_assertion(CGAL::is_closed(tmesh));
-
-    tree_ptr = new AABB_tree(faces(tmesh).first,
-                             faces(tmesh).second,
-                             tmesh);
-  }
 
   /**
   * Constructor that takes a pre-built \cgal `AABB_tree`
@@ -166,11 +168,12 @@ public:
   , tree_ptr(&tree)
   , own_tree(false)
   {
+    box = tree.bbox();
   }
 
   ~Side_of_triangle_mesh()
   {
-    if (own_tree)
+    if (own_tree && tree_ptr!=nullptr)
       delete tree_ptr;
   }
 
@@ -186,9 +189,37 @@ public:
    */
   Bounded_side operator()(const Point& point) const
   {
+    if(point.x() < box.xmin()
+       || point.x() > box.xmax()
+       || point.y() < box.ymin()
+       || point.y() > box.ymax()
+       || point.z() < box.zmin()
+       || point.z() > box.zmax())
+    {
+      return CGAL::ON_UNBOUNDED_SIDE;
+    }
+    else
+    {
+      // Lazily build the tree only when needed
+      if (tree_ptr==nullptr)
+      {
+#ifdef CGAL_HAS_THREADS
+        CGAL_SCOPED_LOCK(tree_mutex);
+#endif
+        CGAL_assertion(tm_ptr != nullptr && opt_vpm!=boost::none);
+        if (tree_ptr==nullptr)
+        {
+          tree_ptr = new AABB_tree(faces(*tm_ptr).first,
+                                   faces(*tm_ptr).second,
+                                   *tm_ptr, *opt_vpm);
+          const_cast<AABB_tree_*>(tree_ptr)->build();
+        }
+      }
+
       return internal::Point_inside_vertical_ray_cast<GeomTraits, AABB_tree>()(
             point, *tree_ptr, ray_functor, vector_functor);
     }
+  }
 
 };