Add support for reentrant qhull

cmake/Modules/FindQhull.cmake

@@ -9,14 +9,14 @@
 # If QHULL_USE_STATIC is specified then look for static libraries ONLY else
 # look for shared ones
 
-set(QHULL_MAJOR_VERSION 6)
+set(QHULL_MAJOR_VERSION 8)
 
 if(QHULL_USE_STATIC)
   set(QHULL_RELEASE_NAME qhullstatic)
   set(QHULL_DEBUG_NAME qhullstatic_d)
 else()
-  set(QHULL_RELEASE_NAME qhull_p qhull${QHULL_MAJOR_VERSION} qhull)
-  set(QHULL_DEBUG_NAME qhull_p_d qhull${QHULL_MAJOR_VERSION}_d qhull_d${QHULL_MAJOR_VERSION} qhull_d)
+  set(QHULL_RELEASE_NAME qhull_r qhull${QHULL_MAJOR_VERSION} qhull)
+  set(QHULL_DEBUG_NAME qhull_r_d qhull${QHULL_MAJOR_VERSION}_d qhull_d${QHULL_MAJOR_VERSION} qhull_d)
 endif()
 
 find_file(QHULL_HEADER

surface/include/pcl/surface/convex_hull.h

@@ -88,11 +88,11 @@ namespace pcl
       using PointCloudConstPtr = typename PointCloud::ConstPtr;
 
       /** \brief Empty constructor. */
-      ConvexHull () : compute_area_ (false), total_area_ (0), total_volume_ (0), dimension_ (0), 
+      ConvexHull () : compute_area_ (false), total_area_ (0), total_volume_ (0), dimension_ (0),
                       projection_angle_thresh_ (std::cos (0.174532925) ), qhull_flags ("qhull "),
                       x_axis_ (1.0, 0.0, 0.0), y_axis_ (0.0, 1.0, 0.0), z_axis_ (0.0, 0.0, 1.0)
       {
-      };
+      }
 
       /** \brief Empty destructor */
       ~ConvexHull () {}

surface/include/pcl/surface/impl/concave_hull.hpp

@@ -194,8 +194,13 @@ pcl::ConcaveHull<PointInT>::performReconstruction (PointCloud &alpha_shape, std:
       points[i * dim_ + 2] = static_cast<coordT> (cloud_transformed[i].z);
   }
 
+  qhT qh_qh;
+  qhT* qh = &qh_qh;
+  QHULL_LIB_CHECK
+  qh_zero(qh, errfile);
+
   // Compute concave hull
-  exitcode = qh_new_qhull (dim_, static_cast<int> (cloud_transformed.size ()), points, ismalloc, flags, outfile, errfile);
+  exitcode = qh_new_qhull (qh, dim_, static_cast<int> (cloud_transformed.size ()), points, ismalloc, flags, outfile, errfile);
 
   if (exitcode != 0)
   {
@@ -227,16 +232,16 @@ pcl::ConcaveHull<PointInT>::performReconstruction (PointCloud &alpha_shape, std:
     alpha_shape.width = alpha_shape.height = 0;
     polygons.resize (0);
 
-    qh_freeqhull (!qh_ALL);
+    qh_freeqhull (qh, !qh_ALL);
     int curlong, totlong;
-    qh_memfreeshort (&curlong, &totlong);
+    qh_memfreeshort (qh, &curlong, &totlong);
 
     return;
   }
 
-  qh_setvoronoi_all ();
+  qh_setvoronoi_all (qh);
 
-  int num_vertices = qh num_vertices;
+  int num_vertices = qh->num_vertices;
   alpha_shape.resize (num_vertices);
 
   vertexT *vertex;
@@ -253,11 +258,11 @@ pcl::ConcaveHull<PointInT>::performReconstruction (PointCloud &alpha_shape, std:
   ++max_vertex_id;
   std::vector<int> qhid_to_pcidx (max_vertex_id);
 
-  int num_facets = qh num_facets;
+  int num_facets = qh->num_facets;
 
   if (dim_ == 3)
   {
-    setT *triangles_set = qh_settemp (4 * num_facets);
+    setT *triangles_set = qh_settemp (qh, 4 * num_facets);
     if (voronoi_centers_)
       voronoi_centers_->points.resize (num_facets);
 
@@ -283,29 +288,29 @@ pcl::ConcaveHull<PointInT>::performReconstruction (PointCloud &alpha_shape, std:
         if (r <= alpha_)
         {
           // all triangles in tetrahedron are good, add them all to the alpha shape (triangles_set)
-          qh_makeridges (facet);
+          qh_makeridges (qh, facet);
           facet->good = true;
-          facet->visitid = qh visit_id;
+          facet->visitid = qh->visit_id;
           ridgeT *ridge, **ridgep;
           FOREACHridge_ (facet->ridges)
           {
             facetT *neighb = otherfacet_ (ridge, facet);
-            if ((neighb->visitid != qh visit_id))
-              qh_setappend (&triangles_set, ridge);
+            if ((neighb->visitid != qh->visit_id))
+              qh_setappend (qh, &triangles_set, ridge);
           }
         }
         else
         {
           // consider individual triangles from the tetrahedron...
           facet->good = false;
-          facet->visitid = qh visit_id;
-          qh_makeridges (facet);
+          facet->visitid = qh->visit_id;
+          qh_makeridges (qh, facet);
           ridgeT *ridge, **ridgep;
           FOREACHridge_ (facet->ridges)
           {
             facetT *neighb;
             neighb = otherfacet_ (ridge, facet);
-            if ((neighb->visitid != qh visit_id))
+            if ((neighb->visitid != qh->visit_id))
             {
               // check if individual triangle is good and add it to triangles_set
 
@@ -322,7 +327,7 @@ pcl::ConcaveHull<PointInT>::performReconstruction (PointCloud &alpha_shape, std:
 
               double r = pcl::getCircumcircleRadius (a, b, c);
               if (r <= alpha_)
-                qh_setappend (&triangles_set, ridge);
+                qh_setappend (qh, &triangles_set, ridge);
             }
           }
         }
@@ -354,7 +359,7 @@ pcl::ConcaveHull<PointInT>::performReconstruction (PointCloud &alpha_shape, std:
       {
         polygons[triangles].vertices.resize (3);
         int vertex_n, vertex_i;
-        FOREACHvertex_i_ ((*ridge).vertices)  //3 vertices per ridge!
+        FOREACHvertex_i_ (qh, (*ridge).vertices)  //3 vertices per ridge!
         {
           if (!added_vertices[vertex->id])
           {
@@ -383,7 +388,7 @@ pcl::ConcaveHull<PointInT>::performReconstruction (PointCloud &alpha_shape, std:
   {
     // Compute the alpha complex for the set of points
     // Filters the delaunay triangles
-    setT *edges_set = qh_settemp (3 * num_facets);
+    setT *edges_set = qh_settemp (qh, 3 * num_facets);
     if (voronoi_centers_)
       voronoi_centers_->points.resize (num_facets);
 
@@ -403,12 +408,12 @@ pcl::ConcaveHull<PointInT>::performReconstruction (PointCloud &alpha_shape, std:
         if (r <= alpha_)
         {
           pcl::Vertices facet_vertices;   //TODO: is not used!!
-          qh_makeridges (facet);
+          qh_makeridges (qh, facet);
           facet->good = true;
 
           ridgeT *ridge, **ridgep;
           FOREACHridge_ (facet->ridges)
-          qh_setappend (&edges_set, ridge);
+          qh_setappend (qh, &edges_set, ridge);
 
           if (voronoi_centers_)
           {
@@ -438,7 +443,7 @@ pcl::ConcaveHull<PointInT>::performReconstruction (PointCloud &alpha_shape, std:
         std::vector<int> pcd_indices;
         pcd_indices.resize (2);
 
-        FOREACHvertex_i_ ((*ridge).vertices)  //in 2-dim, 2 vertices per ridge!
+        FOREACHvertex_i_ (qh, (*ridge).vertices)  //in 2-dim, 2 vertices per ridge!
         {
           if (!added_vertices[vertex->id])
           {
@@ -540,9 +545,9 @@ pcl::ConcaveHull<PointInT>::performReconstruction (PointCloud &alpha_shape, std:
       voronoi_centers_->points.resize (dd);
   }
 
-  qh_freeqhull (!qh_ALL);
+  qh_freeqhull (qh, !qh_ALL);
   int curlong, totlong;
-  qh_memfreeshort (&curlong, &totlong);
+  qh_memfreeshort (qh, &curlong, &totlong);
 
   Eigen::Affine3d transInverse = transform1.inverse ();
   pcl::transformPointCloud (alpha_shape, alpha_shape, transInverse);

surface/include/pcl/surface/impl/convex_hull.hpp

@@ -180,40 +180,45 @@ pcl::ConvexHull<PointInT>::performReconstruction2D (PointCloud &hull, std::vecto
     // This should only happen if we had invalid input
     PCL_ERROR ("[pcl::%s::performReconstruction2D] Invalid input!\n", getClassName ().c_str ());
   }
+
+  qhT qh_qh;
+  qhT* qh = &qh_qh;
+  QHULL_LIB_CHECK
+  qh_zero(qh, errfile);
 
   // Compute convex hull
-  int exitcode = qh_new_qhull (dimension, static_cast<int> (indices_->size ()), points, ismalloc, const_cast<char*> (flags), outfile, errfile);
+  int exitcode = qh_new_qhull (qh, dimension, static_cast<int> (indices_->size ()), points, ismalloc, const_cast<char*> (flags), outfile, errfile);
 #ifdef HAVE_QHULL_2011
   if (compute_area_)
   {
-    qh_prepare_output();
+    qh_prepare_output(qh);
   }
 #endif
 
   // 0 if no error from qhull or it doesn't find any vertices
-  if (exitcode != 0 || qh num_vertices == 0)
+  if (exitcode != 0 || qh->num_vertices == 0)
   {
     PCL_ERROR ("[pcl::%s::performReconstrution2D] ERROR: qhull was unable to compute a convex hull for the given point cloud (%lu)!\n", getClassName ().c_str (), indices_->size ());
 
     hull.resize (0);
     hull.width = hull.height = 0;
     polygons.resize (0);
 
-    qh_freeqhull (!qh_ALL);
+    qh_freeqhull (qh, !qh_ALL);
     int curlong, totlong;
-    qh_memfreeshort (&curlong, &totlong);
+    qh_memfreeshort (qh, &curlong, &totlong);
 
     return;
   }
 
   // Qhull returns the area in volume for 2D
   if (compute_area_)
   {
-    total_area_ = qh totvol;
+    total_area_ = qh->totvol;
     total_volume_ = 0.0;
   }
 
-  int num_vertices = qh num_vertices;
+  int num_vertices = qh->num_vertices;
   hull.resize (num_vertices);
   memset (&hull.points[0], hull.size (), sizeof (PointInT));
 
@@ -226,8 +231,8 @@ pcl::ConvexHull<PointInT>::performReconstruction2D (PointCloud &hull, std::vecto
 
   FORALLvertices
   {
-    hull[i] = (*input_)[(*indices_)[qh_pointid (vertex->point)]];
-    idx_points[i].first = qh_pointid (vertex->point);
+    hull[i] = (*input_)[(*indices_)[qh_pointid (qh, vertex->point)]];
+    idx_points[i].first = qh_pointid (qh, vertex->point);
     ++i;
   }
 
@@ -274,9 +279,9 @@ pcl::ConvexHull<PointInT>::performReconstruction2D (PointCloud &hull, std::vecto
     polygons[0].vertices[j] = static_cast<unsigned int> (j);
   }
 
-  qh_freeqhull (!qh_ALL);
+  qh_freeqhull (qh, !qh_ALL);
   int curlong, totlong;
-  qh_memfreeshort (&curlong, &totlong);
+  qh_memfreeshort (qh, &curlong, &totlong);
 
   hull.width = hull.size ();
   hull.height = 1;
@@ -320,12 +325,19 @@ pcl::ConvexHull<PointInT>::performReconstruction3D (
     points[j + 2] = static_cast<coordT> ((*input_)[(*indices_)[i]].z);
   }
 
+  static FILE* null = fopen("/dev/null", "w");
+
+  qhT qh_qh;
+  qhT* qh = &qh_qh;
+  QHULL_LIB_CHECK
+  qh_zero(qh, null);
+
   // Compute convex hull
-  int exitcode = qh_new_qhull (dimension, static_cast<int> (indices_->size ()), points, ismalloc, const_cast<char*> (flags), outfile, errfile);
+  int exitcode = qh_new_qhull (qh, dimension, static_cast<int> (indices_->size ()), points, ismalloc, const_cast<char*> (flags), outfile, errfile);
 #ifdef HAVE_QHULL_2011
   if (compute_area_)
   {
-    qh_prepare_output();
+    qh_prepare_output(qh);
   }
 #endif
 
@@ -341,18 +353,18 @@ pcl::ConvexHull<PointInT>::performReconstruction3D (
     hull.width = hull.height = 0;
     polygons.resize (0);
 
-    qh_freeqhull (!qh_ALL);
+    qh_freeqhull (qh, !qh_ALL);
     int curlong, totlong;
-    qh_memfreeshort (&curlong, &totlong);
+    qh_memfreeshort (qh, &curlong, &totlong);
 
     return;
   }
 
-  qh_triangulate ();
+  qh_triangulate (qh);
 
-  int num_facets = qh num_facets;
+  int num_facets = qh->num_facets;
 
-  int num_vertices = qh num_vertices;
+  int num_vertices = qh->num_vertices;
   hull.resize (num_vertices);
 
   vertexT * vertex;
@@ -375,7 +387,7 @@ pcl::ConvexHull<PointInT>::performReconstruction3D (
   FORALLvertices
   {
     // Add vertices to hull point_cloud and store index
-    hull_indices_.indices.push_back ((*indices_)[qh_pointid (vertex->point)]);
+    hull_indices_.indices.push_back ((*indices_)[qh_pointid (qh, vertex->point)]);
     hull[i] = (*input_)[hull_indices_.indices.back ()];
 
     qhid_to_pcidx[vertex->id] = i; // map the vertex id of qhull to the point cloud index
@@ -384,8 +396,8 @@ pcl::ConvexHull<PointInT>::performReconstruction3D (
 
   if (compute_area_)
   {
-    total_area_  = qh totarea;
-    total_volume_ = qh totvol;
+    total_area_  = qh->totarea;
+    total_volume_ = qh->totvol;
   }
 
   if (fill_polygon_data)
@@ -400,16 +412,16 @@ pcl::ConvexHull<PointInT>::performReconstruction3D (
 
       // Needed by FOREACHvertex_i_
       int vertex_n, vertex_i;
-      FOREACHvertex_i_ ((*facet).vertices)
+      FOREACHvertex_i_ (qh, (*facet).vertices)
       //facet_vertices.vertices.push_back (qhid_to_pcidx[vertex->id]);
       polygons[dd].vertices[vertex_i] = qhid_to_pcidx[vertex->id];
       ++dd;
     }
   }
   // Deallocates memory (also the points)
-  qh_freeqhull (!qh_ALL);
+  qh_freeqhull (qh, !qh_ALL);
   int curlong, totlong;
-  qh_memfreeshort (&curlong, &totlong);
+  qh_memfreeshort (qh, &curlong, &totlong);
 
   hull.width = hull.size ();
   hull.height = 1;

surface/include/pcl/surface/qhull.h

@@ -49,14 +49,14 @@
 extern "C"
 {
 #ifdef HAVE_QHULL_2011
-#  include "libqhull/libqhull.h"
-#  include "libqhull/mem.h"
-#  include "libqhull/qset.h"
-#  include "libqhull/geom.h"
-#  include "libqhull/merge.h"
-#  include "libqhull/poly.h"
-#  include "libqhull/io.h"
-#  include "libqhull/stat.h"
+#  include "libqhull_r/libqhull_r.h"
+#  include "libqhull_r/mem_r.h"
+#  include "libqhull_r/qset_r.h"
+#  include "libqhull_r/geom_r.h"
+#  include "libqhull_r/merge_r.h"
+#  include "libqhull_r/poly_r.h"
+#  include "libqhull_r/io_r.h"
+#  include "libqhull_r/stat_r.h"
 #else
 #  include "qhull/qhull.h"
 #  include "qhull/mem.h"