diff --git a/cmake/Modules/FindQhull.cmake b/cmake/Modules/FindQhull.cmake
index c9c72f09020..015cee9ba46 100644
--- a/cmake/Modules/FindQhull.cmake
+++ b/cmake/Modules/FindQhull.cmake
@@ -6,44 +6,38 @@
 # QHULL_INCLUDE_DIRS - Directories containing the QHULL include files.
 # QHULL_LIBRARIES - Libraries needed to use QHULL.
 # QHULL_DEFINITIONS - Compiler flags for QHULL.
-# If QHULL_USE_STATIC is specified then look for static libraries ONLY else 
-# look for shared ones
+# If QHULL_USE_SHARED is specified then look for shared libraries ONLY else
+# look for static ones
 
-set(QHULL_MAJOR_VERSION 6)
+set(QHULL_MAJOR_VERSION 7)
 
-if(QHULL_USE_STATIC)
-  set(QHULL_RELEASE_NAME qhullstatic)
-  set(QHULL_DEBUG_NAME qhullstatic_d)
-else(QHULL_USE_STATIC)
-  set(QHULL_RELEASE_NAME qhull_p qhull${QHULL_MAJOR_VERSION} qhull)
-  set(QHULL_DEBUG_NAME qhull_pd qhull${QHULL_MAJOR_VERSION}_d qhull_d${QHULL_MAJOR_VERSION} qhull_d)
-endif(QHULL_USE_STATIC)
+if(QHULL_USE_SHARED)
+  set(QHULL_RELEASE_NAME qhull_r)
+  set(QHULL_DEBUG_NAME qhull_rd)
+else(QHULL_USE_SHARED)
+  set(QHULL_RELEASE_NAME qhullstatic_r)
+  set(QHULL_DEBUG_NAME qhullstatic_rd)
+endif(QHULL_USE_SHARED)
 
 find_file(QHULL_HEADER
-          NAMES libqhull/libqhull.h qhull.h
+          NAMES libqhull_r/libqhull_r.h
           HINTS "${QHULL_ROOT}" "$ENV{QHULL_ROOT}" "${QHULL_INCLUDE_DIR}"
           PATHS "$ENV{PROGRAMFILES}/QHull" "$ENV{PROGRAMW6432}/QHull" 
-          PATH_SUFFIXES qhull src/libqhull libqhull include)
+          PATH_SUFFIXES src/libqhull_r libqhull_r qhull src/libqhull libqhull include)
 
 set(QHULL_HEADER "${QHULL_HEADER}" CACHE INTERNAL "QHull header" FORCE )
 
 if(QHULL_HEADER)
   get_filename_component(qhull_header ${QHULL_HEADER} NAME_WE)
-  if("${qhull_header}" STREQUAL "qhull")
-    set(HAVE_QHULL_2011 OFF)
-    get_filename_component(QHULL_INCLUDE_DIR ${QHULL_HEADER} PATH)
-  elseif("${qhull_header}" STREQUAL "libqhull")
-    set(HAVE_QHULL_2011 ON)
-    get_filename_component(QHULL_INCLUDE_DIR ${QHULL_HEADER} PATH)
-    get_filename_component(QHULL_INCLUDE_DIR ${QHULL_INCLUDE_DIR} PATH)
-  endif()
+  get_filename_component(QHULL_INCLUDE_DIR ${QHULL_HEADER} PATH)
+  get_filename_component(QHULL_INCLUDE_DIR ${QHULL_INCLUDE_DIR} PATH)
 else(QHULL_HEADER)
   set(QHULL_INCLUDE_DIR "QHULL_INCLUDE_DIR-NOTFOUND")
 endif(QHULL_HEADER)
 
 set(QHULL_INCLUDE_DIR "${QHULL_INCLUDE_DIR}" CACHE PATH "QHull include dir." FORCE)
 
-find_library(QHULL_LIBRARY 
+find_library(QHULL_LIBRARY
              NAMES ${QHULL_RELEASE_NAME}
              HINTS "${QHULL_ROOT}" "$ENV{QHULL_ROOT}"
              PATHS "$ENV{PROGRAMFILES}/QHull" "$ENV{PROGRAMW6432}/QHull" 
@@ -51,20 +45,33 @@ find_library(QHULL_LIBRARY
 
 get_filename_component(QHULL_LIBRARY_NAME "${QHULL_LIBRARY}" NAME)
 
-find_library(QHULL_LIBRARY_DEBUG 
+find_library(QHULL_LIBRARY_DEBUG
              NAMES ${QHULL_DEBUG_NAME} ${QHULL_RELEASE_NAME}
              HINTS "${QHULL_ROOT}" "$ENV{QHULL_ROOT}"
              PATHS "$ENV{PROGRAMFILES}/QHull" "$ENV{PROGRAMW6432}/QHull" 
              PATH_SUFFIXES project build bin lib debug/lib)
 
+find_library(QHULL_CPP_LIBRARY
+               NAMES qhullcpp
+               HINTS "${QHULL_ROOT}" "$ENV{QHULL_ROOT}"
+               PATHS "$ENV{PROGRAMFILES}/QHull" "$ENV{PROGRAMW6432}/QHull"
+               PATH_SUFFIXES project build bin lib)
+
+find_library(QHULL_CPP_LIBRARY_DEBUG
+        NAMES qhullcpp_d qhullcpp
+        HINTS "${QHULL_ROOT}" "$ENV{QHULL_ROOT}"
+        PATHS "$ENV{PROGRAMFILES}/QHull" "$ENV{PROGRAMW6432}/QHull"
+        PATH_SUFFIXES project build bin lib)
+
 if(NOT QHULL_LIBRARY_DEBUG)
+  set(QHULL_CPP_LIBRARY_DEBUG ${QHULL_CPP_LIBRARY})
   set(QHULL_LIBRARY_DEBUG ${QHULL_LIBRARY})
 endif(NOT QHULL_LIBRARY_DEBUG)
 
 get_filename_component(QHULL_LIBRARY_DEBUG_NAME "${QHULL_LIBRARY_DEBUG}" NAME)
 
 set(QHULL_INCLUDE_DIRS ${QHULL_INCLUDE_DIR})
-set(QHULL_LIBRARIES optimized ${QHULL_LIBRARY} debug ${QHULL_LIBRARY_DEBUG})
+set(QHULL_LIBRARIES optimized ${QHULL_CPP_LIBRARY} ${QHULL_LIBRARY} debug ${QHULL_CPP_LIBRARY_DEBUG} ${QHULL_LIBRARY_DEBUG})
 
 include(FindPackageHandleStandardArgs)
 find_package_handle_standard_args(Qhull DEFAULT_MSG QHULL_LIBRARY QHULL_INCLUDE_DIR)
@@ -73,11 +80,5 @@ mark_as_advanced(QHULL_LIBRARY QHULL_LIBRARY_DEBUG QHULL_INCLUDE_DIR)
 
 if(QHULL_FOUND)
   set(HAVE_QHULL ON)
-  if(NOT QHULL_USE_STATIC)
-    add_definitions("-Dqh_QHpointer")
-    if(MSVC)
-      add_definitions("-Dqh_QHpointer_dllimport")
-    endif(MSVC)
-  endif(NOT QHULL_USE_STATIC)
   message(STATUS "QHULL found (include: ${QHULL_INCLUDE_DIRS}, lib: ${QHULL_LIBRARIES})")
 endif(QHULL_FOUND)
diff --git a/pcl_config.h.in b/pcl_config.h.in
index a654294f7be..cc354047c1f 100644
--- a/pcl_config.h.in
+++ b/pcl_config.h.in
@@ -21,8 +21,6 @@
 
 #cmakedefine HAVE_QHULL 1
 
-#cmakedefine HAVE_QHULL_2011 1
-
 #cmakedefine HAVE_CUDA 1
 
 #cmakedefine HAVE_FZAPI 1
diff --git a/surface/include/pcl/surface/impl/concave_hull.hpp b/surface/include/pcl/surface/impl/concave_hull.hpp
index faaa96fa4ea..bba80de569a 100644
--- a/surface/include/pcl/surface/impl/concave_hull.hpp
+++ b/surface/include/pcl/surface/impl/concave_hull.hpp
@@ -172,416 +172,412 @@ pcl::ConcaveHull<PointInT>::performReconstruction (PointCloud &alpha_shape, std:
   pcl::demeanPointCloud (*input_, *indices_, xyz_centroid, cloud_transformed);
   pcl::transformPointCloud (cloud_transformed, cloud_transformed, transform1);
 
-  // True if qhull should free points in qh_freeqhull() or reallocation
-  boolT ismalloc = True;
   // option flags for qhull, see qh_opt.htm
-  char flags[] = "qhull d QJ";
-  // output from qh_produce_output(), use NULL to skip qh_produce_output()
-  FILE *outfile = NULL;
-  // error messages from qhull code
-  FILE *errfile = stderr;
-  // 0 if no error from qhull
-  int exitcode;
+  std::string flags = "d QJ";
 
   // Array of coordinates for each point
-  coordT *points = reinterpret_cast<coordT*> (calloc (cloud_transformed.points.size () * dim_, sizeof(coordT)));
+  std::vector<double> points (cloud_transformed.points.size () * dim_);
 
   for (size_t i = 0; i < cloud_transformed.points.size (); ++i)
   {
-    points[i * dim_ + 0] = static_cast<coordT> (cloud_transformed.points[i].x);
-    points[i * dim_ + 1] = static_cast<coordT> (cloud_transformed.points[i].y);
+    points[i * dim_ + 0] = cloud_transformed.points[i].x;
+    points[i * dim_ + 1] = cloud_transformed.points[i].y;
 
     if (dim_ > 2)
-      points[i * dim_ + 2] = static_cast<coordT> (cloud_transformed.points[i].z);
+      points[i * dim_ + 2] = cloud_transformed.points[i].z;
   }
 
   // Compute concave hull
-  exitcode = qh_new_qhull (dim_, static_cast<int> (cloud_transformed.points.size ()), points, ismalloc, flags, outfile, errfile);
-
-  if (exitcode != 0)
-  {
-    PCL_ERROR ("[pcl::%s::performReconstrution] ERROR: qhull was unable to compute a concave hull for the given point cloud (%lu)!\n", getClassName ().c_str (), cloud_transformed.points.size ());
-
-    //check if it fails because of NaN values...
-    if (!cloud_transformed.is_dense)
+  try {
+    orgQhull::Qhull qhull;
+    qhull.runQhull("", dim_, cloud_transformed.size (), points.data(), flags.c_str());
+  
+    if (qhull.qhullStatus() != 0)
     {
-      bool NaNvalues = false;
-      for (size_t i = 0; i < cloud_transformed.size (); ++i)
-      {
-        if (!pcl_isfinite (cloud_transformed.points[i].x) ||
-            !pcl_isfinite (cloud_transformed.points[i].y) ||
-            !pcl_isfinite (cloud_transformed.points[i].z))
-        {
-          NaNvalues = true;
-          break;
-        }
-      }
-
-      if (NaNvalues)
-        PCL_ERROR ("[pcl::%s::performReconstruction] ERROR: point cloud contains NaN values, consider running pcl::PassThrough filter first to remove NaNs!\n", getClassName ().c_str ());
+      throw std::runtime_error(qhull.qhullMessage());
     }
-
-    alpha_shape.points.resize (0);
-    alpha_shape.width = alpha_shape.height = 0;
-    polygons.resize (0);
-
-    qh_freeqhull (!qh_ALL);
-    int curlong, totlong;
-    qh_memfreeshort (&curlong, &totlong);
-
-    return;
-  }
-
-  qh_setvoronoi_all ();
-
-  int num_vertices = qh num_vertices;
-  alpha_shape.points.resize (num_vertices);
-
-  vertexT *vertex;
-  // Max vertex id
-  int max_vertex_id = 0;
-  FORALLvertices
-  {
-    if (vertex->id + 1 > max_vertex_id)
-      max_vertex_id = vertex->id + 1;
-  }
-
-  facetT *facet;    // set by FORALLfacets
-
-  ++max_vertex_id;
-  std::vector<int> qhid_to_pcidx (max_vertex_id);
-
-  int num_facets = qh num_facets;
-  int dd = 0;
-
-  if (dim_ == 3)
-  {
-    setT *triangles_set = qh_settemp (4 * num_facets);
-    if (voronoi_centers_)
-      voronoi_centers_->points.resize (num_facets);
-
-    int non_upper = 0;
-    FORALLfacets
+  
+    qh_setvoronoi_all(qhull.qh());
+  
+    int num_vertices = qhull.vertexCount();
+    alpha_shape.points.resize (num_vertices);
+  
+    // Max vertex id
+    int max_vertex_id = 0;
+    orgQhull::QhullVertexList vertexs = qhull.vertexList();
+    for(orgQhull::QhullVertexList::iterator it = vertexs.begin();it!=vertexs.end();it++)
+    {
+      if (it->id() + 1 > max_vertex_id)
+        max_vertex_id = it->id() + 1;
+    }
+  
+    ++max_vertex_id;
+    std::vector<int> qhid_to_pcidx (max_vertex_id);
+  
+    int num_facets = qhull.facetCount();
+    int dd = 0;
+  
+    orgQhull::QhullFacetList facets = qhull.facetList();
+  
+    if (dim_ == 3)
     {
-      // Facets are tetrahedrons (3d)
-      if (!facet->upperdelaunay)
+      std::vector<orgQhull::QhullRidge> triangles_set;
+      if (voronoi_centers_)
+        voronoi_centers_->points.resize (num_facets);
+  
+      int non_upper = 0;
+      for(orgQhull::QhullFacetList::iterator it = facets.begin();it!=facets.end();it++)
       {
-        vertexT *anyVertex = static_cast<vertexT*> (facet->vertices->e[0].p);
-        double *center = facet->center;
-        double r = qh_pointdist (anyVertex->point,center,dim_);
-        facetT *neighb;
-
-        if (voronoi_centers_)
-        {
-          voronoi_centers_->points[non_upper].x = static_cast<float> (facet->center[0]);
-          voronoi_centers_->points[non_upper].y = static_cast<float> (facet->center[1]);
-          voronoi_centers_->points[non_upper].z = static_cast<float> (facet->center[2]);
-        }
-
-        non_upper++;
-
-        if (r <= alpha_)
+        // Facets are tetrahedrons (3d)
+        if (!it->isUpperDelaunay())
         {
-          // all triangles in tetrahedron are good, add them all to the alpha shape (triangles_set)
-          qh_makeridges (facet);
-          facet->good = true;
-          facet->visitid = qh visit_id;
-          ridgeT *ridge, **ridgep;
-          FOREACHridge_ (facet->ridges)
+          orgQhull::QhullPoint anyVertex = it->vertices().first().point();
+          orgQhull::QhullPoint center = it->getCenter();
+          anyVertex.setDimension(center.dimension());
+          double r = center.distance(anyVertex);
+  
+          if (voronoi_centers_)
           {
-            neighb = otherfacet_ (ridge, facet);
-            if ((neighb->visitid != qh visit_id))
-              qh_setappend (&triangles_set, ridge);
+            voronoi_centers_->points[non_upper].x = static_cast<float> (center[0]);
+            voronoi_centers_->points[non_upper].y = static_cast<float> (center[1]);
+            voronoi_centers_->points[non_upper].z = static_cast<float> (center[2]);
           }
-        }
-        else
-        {
-          // consider individual triangles from the tetrahedron...
-          facet->good = false;
-          facet->visitid = qh visit_id;
-          qh_makeridges (facet);
-          ridgeT *ridge, **ridgep;
-          FOREACHridge_ (facet->ridges)
+  
+          non_upper++;
+  
+          qh_makeridges(qhull.qh(), it->getFacetT());
+          orgQhull::QhullRidgeSet ridges = it->ridges();
+  
+          if (r <= alpha_)
           {
-            facetT *neighb;
-            neighb = otherfacet_ (ridge, facet);
-            if ((neighb->visitid != qh visit_id))
+            // all triangles in tetrahedron are good, add them all to the alpha shape (triangles_set)
+            it->getFacetT()->good = true;
+            it->getFacetT()->visitid = qhull.qh()->visit_id;
+            for(orgQhull::QhullRidgeSet::iterator rt = ridges.begin(); rt != ridges.end(); rt++)
             {
-              // check if individual triangle is good and add it to triangles_set
-
-              PointInT a, b, c;
-              a.x = static_cast<float> ((static_cast<vertexT*>(ridge->vertices->e[0].p))->point[0]);
-              a.y = static_cast<float> ((static_cast<vertexT*>(ridge->vertices->e[0].p))->point[1]);
-              a.z = static_cast<float> ((static_cast<vertexT*>(ridge->vertices->e[0].p))->point[2]);
-              b.x = static_cast<float> ((static_cast<vertexT*>(ridge->vertices->e[1].p))->point[0]);
-              b.y = static_cast<float> ((static_cast<vertexT*>(ridge->vertices->e[1].p))->point[1]);
-              b.z = static_cast<float> ((static_cast<vertexT*>(ridge->vertices->e[1].p))->point[2]);
-              c.x = static_cast<float> ((static_cast<vertexT*>(ridge->vertices->e[2].p))->point[0]);
-              c.y = static_cast<float> ((static_cast<vertexT*>(ridge->vertices->e[2].p))->point[1]);
-              c.z = static_cast<float> ((static_cast<vertexT*>(ridge->vertices->e[2].p))->point[2]);
-
-              double r = pcl::getCircumcircleRadius (a, b, c);
-              if (r <= alpha_)
-                qh_setappend (&triangles_set, ridge);
+              orgQhull::QhullRidge ridge = *rt;
+              orgQhull::QhullFacet neighb = ridge.otherFacet(*it);
+              if (neighb.getFacetT()->visitid != qhull.qh()->visit_id)
+                triangles_set.push_back(ridge);
+            }
+          }
+          else
+          {
+            // consider individual triangles from the tetrahedron...
+            it->getFacetT()->good = false;
+            it->getFacetT()->visitid = qhull.qh()->visit_id;
+            for(orgQhull::QhullRidgeSet::iterator rt = ridges.begin(); rt != ridges.end(); rt++)
+            {
+              orgQhull::QhullRidge ridge = *rt;
+              orgQhull::QhullFacet neighb = ridge.otherFacet(*it);
+              if (neighb.getFacetT()->visitid != qhull.qh()->visit_id)
+              {
+                // check if individual triangle is good and add it to triangles_set
+                orgQhull::QhullVertexSet ridgeVertices = ridge.vertices();
+  
+                PointInT a, b, c;
+                a.x = static_cast<float> (ridgeVertices[0].point()[0]);
+                a.y = static_cast<float> (ridgeVertices[0].point()[1]);
+                a.z = static_cast<float> (ridgeVertices[0].point()[2]);
+                b.x = static_cast<float> (ridgeVertices[1].point()[0]);
+                b.y = static_cast<float> (ridgeVertices[1].point()[1]);
+                b.z = static_cast<float> (ridgeVertices[1].point()[2]);
+                c.x = static_cast<float> (ridgeVertices[2].point()[0]);
+                c.y = static_cast<float> (ridgeVertices[2].point()[1]);
+                c.z = static_cast<float> (ridgeVertices[2].point()[2]);
+  
+                double r = pcl::getCircumcircleRadius (a, b, c);
+                if (r <= alpha_)
+                  triangles_set.push_back(ridge);
+              }
             }
           }
         }
       }
-    }
-
-    if (voronoi_centers_)
-      voronoi_centers_->points.resize (non_upper);
-
-    // filter, add points to alpha_shape and create polygon structure
-
-    int num_good_triangles = 0;
-    ridgeT *ridge, **ridgep;
-    FOREACHridge_ (triangles_set)
-    {
-      if (ridge->bottom->upperdelaunay || ridge->top->upperdelaunay || !ridge->top->good || !ridge->bottom->good)
-        num_good_triangles++;
-    }
-
-    polygons.resize (num_good_triangles);
-
-    int vertices = 0;
-    std::vector<bool> added_vertices (max_vertex_id, false);
-
-    int triangles = 0;
-    FOREACHridge_ (triangles_set)
-    {
-      if (ridge->bottom->upperdelaunay || ridge->top->upperdelaunay || !ridge->top->good || !ridge->bottom->good)
+  
+      if (voronoi_centers_)
+        voronoi_centers_->points.resize (non_upper);
+  
+      // filter, add points to alpha_shape and create polygon structure
+  
+      int num_good_triangles = 0;
+      for(std::vector<orgQhull::QhullRidge>::iterator it = triangles_set.begin();it!=triangles_set.end();it++)
       {
-        polygons[triangles].vertices.resize (3);
-        int vertex_n, vertex_i;
-        FOREACHvertex_i_ ((*ridge).vertices)  //3 vertices per ridge!
+        orgQhull::QhullRidge& ridge = *it;
+        if (ridge.bottomFacet().isUpperDelaunay() || ridge.topFacet().isUpperDelaunay() || !ridge.topFacet().isGood() || !ridge.bottomFacet().isGood())
+          num_good_triangles++;
+      }
+  
+      polygons.resize (num_good_triangles);
+  
+      int vertices = 0;
+      std::vector<bool> added_vertices (max_vertex_id, false);
+  
+      int triangles = 0;
+      for(std::vector<orgQhull::QhullRidge>::iterator it = triangles_set.begin();it!=triangles_set.end();it++)
+      {
+        orgQhull::QhullRidge& ridge = *it;
+        if (ridge.bottomFacet().isUpperDelaunay() || ridge.topFacet().isUpperDelaunay() || !ridge.topFacet().isGood() || !ridge.bottomFacet().isGood())
         {
-          if (!added_vertices[vertex->id])
+          polygons[triangles].vertices.resize (3);
+          int vertex_n, vertex_i;
+          orgQhull::QhullVertexSet ridgeVertices = ridge.vertices();
+          for(orgQhull::QhullVertexSet::iterator jt = ridgeVertices.begin();jt!=ridgeVertices.end();jt++)  //3 vertices per ridge!
           {
-            alpha_shape.points[vertices].x = static_cast<float> (vertex->point[0]);
-            alpha_shape.points[vertices].y = static_cast<float> (vertex->point[1]);
-            alpha_shape.points[vertices].z = static_cast<float> (vertex->point[2]);
-
-            qhid_to_pcidx[vertex->id] = vertices;   //map the vertex id of qhull to the point cloud index
-            added_vertices[vertex->id] = true;
-            vertices++;
+            orgQhull::QhullVertex vertex = *jt;
+            if (!added_vertices[vertex.id()])
+            {
+              orgQhull::QhullPoint point = vertex.point();
+              alpha_shape.points[vertices].x = static_cast<float> (point[0]);
+              alpha_shape.points[vertices].y = static_cast<float> (point[1]);
+              alpha_shape.points[vertices].z = static_cast<float> (point[2]);
+  
+              qhid_to_pcidx[vertex.id()] = vertices;   //map the vertex id of qhull to the point cloud index
+              added_vertices[vertex.id()] = true;
+              vertices++;
+            }
+  
+            polygons[triangles].vertices[vertex_i] = qhid_to_pcidx[vertex.id()];
+  
           }
-
-          polygons[triangles].vertices[vertex_i] = qhid_to_pcidx[vertex->id];
-
+  
+          triangles++;
         }
-
-        triangles++;
       }
+  
+      alpha_shape.points.resize (vertices);
+      alpha_shape.width = static_cast<uint32_t> (alpha_shape.points.size ());
+      alpha_shape.height = 1;
     }
-
-    alpha_shape.points.resize (vertices);
-    alpha_shape.width = static_cast<uint32_t> (alpha_shape.points.size ());
-    alpha_shape.height = 1;
-  }
-  else
-  {
-    // Compute the alpha complex for the set of points
-    // Filters the delaunay triangles
-    setT *edges_set = qh_settemp (3 * num_facets);
-    if (voronoi_centers_)
-      voronoi_centers_->points.resize (num_facets);
-
-    FORALLfacets
+    else
     {
-      // Facets are the delaunay triangles (2d)
-      if (!facet->upperdelaunay)
+      // Compute the alpha complex for the set of points
+      // Filters the delaunay triangles
+      std::vector<orgQhull::QhullRidge> edges_set;
+      if (voronoi_centers_)
+        voronoi_centers_->points.resize (num_facets);
+  
+      for(orgQhull::QhullFacetList::iterator it = facets.begin();it!=facets.end();it++)
       {
-        // Check if the distance from any vertex to the facet->center
-        // (center of the voronoi cell) is smaller than alpha
-        vertexT *anyVertex = static_cast<vertexT*>(facet->vertices->e[0].p);
-        double r = (sqrt ((anyVertex->point[0] - facet->center[0]) *
-                          (anyVertex->point[0] - facet->center[0]) +
-                          (anyVertex->point[1] - facet->center[1]) *
-                          (anyVertex->point[1] - facet->center[1])));
-        if (r <= alpha_)
+        // Facets are the delaunay triangles (2d)
+        if (!it->isUpperDelaunay())
         {
-          pcl::Vertices facet_vertices;   //TODO: is not used!!
-          qh_makeridges (facet);
-          facet->good = true;
-
-          ridgeT *ridge, **ridgep;
-          FOREACHridge_ (facet->ridges)
-          qh_setappend (&edges_set, ridge);
-
-          if (voronoi_centers_)
+          orgQhull::QhullVertexSet vertices = it->vertices();
+  
+          // Check if the distance from any vertex to the facet->center
+          // (center of the voronoi cell) is smaller than alpha
+          orgQhull::QhullPoint anyVertex = it->vertices().first().point();
+          orgQhull::QhullPoint center = it->getCenter();
+          anyVertex.setDimension(center.dimension());
+          double r = center.distance(anyVertex);
+  
+          if (r <= alpha_)
           {
-            voronoi_centers_->points[dd].x = static_cast<float> (facet->center[0]);
-            voronoi_centers_->points[dd].y = static_cast<float> (facet->center[1]);
-            voronoi_centers_->points[dd].z = 0.0f;
+            qh_makeridges(qhull.qh(), it->getFacetT());
+            orgQhull::QhullRidgeSet ridges = it->ridges();
+            for(orgQhull::QhullRidgeSet::iterator rt = ridges.begin(); rt != ridges.end(); rt++) {
+              edges_set.push_back(*rt);
+            }
+  
+            if (voronoi_centers_)
+            {
+              voronoi_centers_->points[dd].x = static_cast<float> (center[0]);
+              voronoi_centers_->points[dd].y = static_cast<float> (center[1]);
+              voronoi_centers_->points[dd].z = 0.0f;
+            }
+  
+            ++dd;
           }
-
-          ++dd;
+          else
+            it->getFacetT()->good = false;
         }
-        else
-          facet->good = false;
       }
-    }
-
-    int vertices = 0;
-    std::vector<bool> added_vertices (max_vertex_id, false);
-    std::map<int, std::vector<int> > edges;
-
-    ridgeT *ridge, **ridgep;
-    FOREACHridge_ (edges_set)
-    {
-      if (ridge->bottom->upperdelaunay || ridge->top->upperdelaunay || !ridge->top->good || !ridge->bottom->good)
+  
+      int vertices = 0;
+      std::vector<bool> added_vertices (max_vertex_id, false);
+      std::map<int, std::vector<int> > edges;
+  
+      for(std::vector<orgQhull::QhullRidge>::iterator rt = edges_set.begin(); rt != edges_set.end(); rt++)
       {
-        int vertex_n, vertex_i;
-        int vertices_in_ridge=0;
-        std::vector<int> pcd_indices;
-        pcd_indices.resize (2);
-
-        FOREACHvertex_i_ ((*ridge).vertices)  //in 2-dim, 2 vertices per ridge!
+        orgQhull::QhullRidge& ridge = *rt;
+        if (ridge.bottomFacet().isUpperDelaunay() || ridge.topFacet().isUpperDelaunay() || !ridge.topFacet().isGood() || !ridge.bottomFacet().isGood())
         {
-          if (!added_vertices[vertex->id])
+          orgQhull::QhullVertexSet ridgeVertices = ridge.vertices();
+          int vertices_in_ridge=0;
+          std::vector<int> pcd_indices;
+          pcd_indices.resize (2);
+  
+          for(orgQhull::QhullVertexSet::iterator jt = ridgeVertices.begin();jt!=ridgeVertices.end();jt++)  //in 2-dim, 2 vertices per ridge!
           {
-            alpha_shape.points[vertices].x = static_cast<float> (vertex->point[0]);
-            alpha_shape.points[vertices].y = static_cast<float> (vertex->point[1]);
-
-            if (dim_ > 2)
-              alpha_shape.points[vertices].z = static_cast<float> (vertex->point[2]);
+            orgQhull::QhullVertex vertex = *jt;
+            if (!added_vertices[vertex.id()])
+            {
+              orgQhull::QhullPoint point = vertex.point();
+              alpha_shape.points[vertices].x = static_cast<float> (point[0]);
+              alpha_shape.points[vertices].y = static_cast<float> (point[1]);
+  
+              if (dim_ > 2)
+                alpha_shape.points[vertices].z = static_cast<float> (point[2]);
+              else
+                alpha_shape.points[vertices].z = 0;
+  
+              qhid_to_pcidx[vertex.id()] = vertices;   //map the vertex id of qhull to the point cloud index
+              added_vertices[vertex.id()] = true;
+              pcd_indices[vertices_in_ridge] = vertices;
+              vertices++;
+            }
             else
-              alpha_shape.points[vertices].z = 0;
-
-            qhid_to_pcidx[vertex->id] = vertices;   //map the vertex id of qhull to the point cloud index
-            added_vertices[vertex->id] = true;
-            pcd_indices[vertices_in_ridge] = vertices;
-            vertices++;
-          }
-          else
-          {
-            pcd_indices[vertices_in_ridge] = qhid_to_pcidx[vertex->id];
+            {
+              pcd_indices[vertices_in_ridge] = qhid_to_pcidx[vertex.id()];
+            }
+  
+            vertices_in_ridge++;
           }
-
-          vertices_in_ridge++;
+  
+          // make edges bidirectional and pointing to alpha_shape pointcloud...
+          edges[pcd_indices[0]].push_back (pcd_indices[1]);
+          edges[pcd_indices[1]].push_back (pcd_indices[0]);
         }
-
-        // make edges bidirectional and pointing to alpha_shape pointcloud...
-        edges[pcd_indices[0]].push_back (pcd_indices[1]);
-        edges[pcd_indices[1]].push_back (pcd_indices[0]);
       }
-    }
-
-    alpha_shape.points.resize (vertices);
-
-    std::vector<std::vector<int> > connected;
-    PointCloud alpha_shape_sorted;
-    alpha_shape_sorted.points.resize (vertices);
-
-    // iterate over edges until they are empty!
-    std::map<int, std::vector<int> >::iterator curr = edges.begin ();
-    int next = - 1;
-    std::vector<bool> used (vertices, false);   // used to decide which direction should we take!
-    std::vector<int> pcd_idx_start_polygons;
-    pcd_idx_start_polygons.push_back (0);
-
-    // start following edges and removing elements
-    int sorted_idx = 0;
-    while (!edges.empty ())
-    {
-      alpha_shape_sorted.points[sorted_idx] = alpha_shape.points[(*curr).first];
-      // check where we can go from (*curr).first
-      for (size_t i = 0; i < (*curr).second.size (); i++)
+  
+      alpha_shape.points.resize (vertices);
+  
+      std::vector<std::vector<int> > connected;
+      PointCloud alpha_shape_sorted;
+      alpha_shape_sorted.points.resize (vertices);
+  
+      // iterate over edges until they are empty!
+      std::map<int, std::vector<int> >::iterator curr = edges.begin ();
+      int next = - 1;
+      std::vector<bool> used (vertices, false);   // used to decide which direction should we take!
+      std::vector<int> pcd_idx_start_polygons;
+      pcd_idx_start_polygons.push_back (0);
+  
+      // start following edges and removing elements
+      int sorted_idx = 0;
+      while (!edges.empty ())
       {
-        if (!used[((*curr).second)[i]])
+        alpha_shape_sorted.points[sorted_idx] = alpha_shape.points[(*curr).first];
+        // check where we can go from (*curr).first
+        for (size_t i = 0; i < (*curr).second.size (); i++)
         {
-          // we can go there
-          next = ((*curr).second)[i];
+          if (!used[((*curr).second)[i]])
+          {
+            // we can go there
+            next = ((*curr).second)[i];
+            break;
+          }
+        }
+  
+        used[(*curr).first] = true;
+        edges.erase (curr);   // remove edges starting from curr
+  
+        sorted_idx++;
+  
+        if (edges.empty ())
           break;
+  
+        // reassign current
+        curr = edges.find (next);   // if next is not found, then we have unconnected polygons.
+        if (curr == edges.end ())
+        {
+          // set current to any of the remaining in edge!
+          curr = edges.begin ();
+          pcd_idx_start_polygons.push_back (sorted_idx);
         }
       }
-
-      used[(*curr).first] = true;
-      edges.erase (curr);   // remove edges starting from curr
-
-      sorted_idx++;
-
-      if (edges.empty ())
-        break;
-
-      // reassign current
-      curr = edges.find (next);   // if next is not found, then we have unconnected polygons.
-      if (curr == edges.end ())
+  
+      pcd_idx_start_polygons.push_back (sorted_idx);
+  
+      alpha_shape.points = alpha_shape_sorted.points;
+  
+      polygons.reserve (pcd_idx_start_polygons.size () - 1);
+  
+      for (size_t poly_id = 0; poly_id < pcd_idx_start_polygons.size () - 1; poly_id++)
       {
-        // set current to any of the remaining in edge!
-        curr = edges.begin ();
-        pcd_idx_start_polygons.push_back (sorted_idx);
+        // Check if we actually have a polygon, and not some degenerated output from QHull
+        if (pcd_idx_start_polygons[poly_id + 1] - pcd_idx_start_polygons[poly_id] >= 3)
+        {
+          pcl::Vertices vertices;
+          vertices.vertices.resize (pcd_idx_start_polygons[poly_id + 1] - pcd_idx_start_polygons[poly_id]);
+          // populate points in the corresponding polygon
+          for (int j = pcd_idx_start_polygons[poly_id]; j < pcd_idx_start_polygons[poly_id + 1]; ++j)
+            vertices.vertices[j - pcd_idx_start_polygons[poly_id]] = static_cast<uint32_t> (j);
+  
+          polygons.push_back (vertices);
+        }
       }
+  
+      if (voronoi_centers_)
+        voronoi_centers_->points.resize (dd);
     }
-
-    pcd_idx_start_polygons.push_back (sorted_idx);
-
-    alpha_shape.points = alpha_shape_sorted.points;
-
-    polygons.reserve (pcd_idx_start_polygons.size () - 1);
-
-    for (size_t poly_id = 0; poly_id < pcd_idx_start_polygons.size () - 1; poly_id++)
+  
+    Eigen::Affine3d transInverse = transform1.inverse ();
+    pcl::transformPointCloud (alpha_shape, alpha_shape, transInverse);
+    xyz_centroid[0] = - xyz_centroid[0];
+    xyz_centroid[1] = - xyz_centroid[1];
+    xyz_centroid[2] = - xyz_centroid[2];
+    pcl::demeanPointCloud (alpha_shape, xyz_centroid, alpha_shape);
+  
+    // also transform voronoi_centers_...
+    if (voronoi_centers_)
+    {
+      pcl::transformPointCloud (*voronoi_centers_, *voronoi_centers_, transInverse);
+      pcl::demeanPointCloud (*voronoi_centers_, xyz_centroid, *voronoi_centers_);
+    }
+  
+    if (keep_information_)
     {
-      // Check if we actually have a polygon, and not some degenerated output from QHull
-      if (pcd_idx_start_polygons[poly_id + 1] - pcd_idx_start_polygons[poly_id] >= 3)
+      // build a tree with the original points
+      pcl::KdTreeFLANN<PointInT> tree (true);
+      tree.setInputCloud (input_, indices_);
+  
+      std::vector<int> neighbor;
+      std::vector<float> distances;
+      neighbor.resize (1);
+      distances.resize (1);
+  
+      // for each point in the concave hull, search for the nearest neighbor in the original point cloud
+      hull_indices_.header = input_->header;
+      hull_indices_.indices.clear ();
+      hull_indices_.indices.reserve (alpha_shape.points.size ());
+  
+      for (size_t i = 0; i < alpha_shape.points.size (); i++)
       {
-        pcl::Vertices vertices;
-        vertices.vertices.resize (pcd_idx_start_polygons[poly_id + 1] - pcd_idx_start_polygons[poly_id]);
-        // populate points in the corresponding polygon
-        for (int j = pcd_idx_start_polygons[poly_id]; j < pcd_idx_start_polygons[poly_id + 1]; ++j)
-          vertices.vertices[j - pcd_idx_start_polygons[poly_id]] = static_cast<uint32_t> (j);
-
-        polygons.push_back (vertices);
+        tree.nearestKSearch (alpha_shape.points[i], 1, neighbor, distances);
+        hull_indices_.indices.push_back (neighbor[0]);
       }
+  
+      // replace point with the closest neighbor in the original point cloud
+      pcl::copyPointCloud (*input_, hull_indices_.indices, alpha_shape);
     }
-
-    if (voronoi_centers_)
-      voronoi_centers_->points.resize (dd);
-  }
-
-  qh_freeqhull (!qh_ALL);
-  int curlong, totlong;
-  qh_memfreeshort (&curlong, &totlong);
-
-  Eigen::Affine3d transInverse = transform1.inverse ();
-  pcl::transformPointCloud (alpha_shape, alpha_shape, transInverse);
-  xyz_centroid[0] = - xyz_centroid[0];
-  xyz_centroid[1] = - xyz_centroid[1];
-  xyz_centroid[2] = - xyz_centroid[2];
-  pcl::demeanPointCloud (alpha_shape, xyz_centroid, alpha_shape);
-
-  // also transform voronoi_centers_...
-  if (voronoi_centers_)
-  {
-    pcl::transformPointCloud (*voronoi_centers_, *voronoi_centers_, transInverse);
-    pcl::demeanPointCloud (*voronoi_centers_, xyz_centroid, *voronoi_centers_);
   }
-
-  if (keep_information_)
+  catch (const std::exception& e)
   {
-    // build a tree with the original points
-    pcl::KdTreeFLANN<PointInT> tree (true);
-    tree.setInputCloud (input_, indices_);
-
-    std::vector<int> neighbor;
-    std::vector<float> distances;
-    neighbor.resize (1);
-    distances.resize (1);
-
-    // for each point in the concave hull, search for the nearest neighbor in the original point cloud
-    hull_indices_.header = input_->header;
-    hull_indices_.indices.clear ();
-    hull_indices_.indices.reserve (alpha_shape.points.size ());
-
-    for (size_t i = 0; i < alpha_shape.points.size (); i++)
+    PCL_ERROR ("[pcl::%s::performReconstrution] ERROR: qhull was unable to compute a concave hull for the given point cloud (%lu)!\nError: %s\n", getClassName ().c_str (), cloud_transformed.points.size (), e.what());
+  
+    //check if it fails because of NaN values...
+    if (!cloud_transformed.is_dense)
     {
-      tree.nearestKSearch (alpha_shape.points[i], 1, neighbor, distances);
-      hull_indices_.indices.push_back (neighbor[0]);
+      bool NaNvalues = false;
+      for (size_t i = 0; i < cloud_transformed.size (); ++i)
+      {
+        if (!pcl_isfinite (cloud_transformed.points[i].x) ||
+            !pcl_isfinite (cloud_transformed.points[i].y) ||
+            !pcl_isfinite (cloud_transformed.points[i].z))
+        {
+          NaNvalues = true;
+          break;
+        }
+      }
+    
+      if (NaNvalues)
+        PCL_ERROR ("[pcl::%s::performReconstruction] ERROR: point cloud contains NaN values, consider running pcl::PassThrough filter first to remove NaNs!\n", getClassName ().c_str ());
     }
-
-    // replace point with the closest neighbor in the original point cloud
-    pcl::copyPointCloud (*input_, hull_indices_.indices, alpha_shape);
+    
+    alpha_shape.points.resize (0);
+    alpha_shape.width = alpha_shape.height = 0;
+    polygons.resize (0);
   }
 }
 #ifdef __GNUC__
diff --git a/surface/include/pcl/surface/impl/convex_hull.hpp b/surface/include/pcl/surface/impl/convex_hull.hpp
index fb0af36f32a..b4d6c5b593b 100644
--- a/surface/include/pcl/surface/impl/convex_hull.hpp
+++ b/surface/include/pcl/surface/impl/convex_hull.hpp
@@ -131,135 +131,77 @@ pcl::ConvexHull<PointInT>::performReconstruction2D (PointCloud &hull, std::vecto
     yz_proj_safe = false;
   }
 
-  // True if qhull should free points in qh_freeqhull() or reallocation
-  boolT ismalloc = True;
-  // output from qh_produce_output(), use NULL to skip qh_produce_output()
-  FILE *outfile = NULL;
-
-#ifndef HAVE_QHULL_2011
-  if (compute_area_)
-    outfile = stderr;
-#endif
-
-  // option flags for qhull, see qh_opt.htm
-  const char* flags = qhull_flags.c_str ();
-  // error messages from qhull code
-  FILE *errfile = stderr;
-
   // Array of coordinates for each point
-  coordT *points = reinterpret_cast<coordT*> (calloc (indices_->size () * dimension, sizeof (coordT)));
+  std::vector<double> points (2 * indices_->size ());
 
-  // Build input data, using appropriate projection
-  int j = 0;
-  if (xy_proj_safe)
-  {
-    for (size_t i = 0; i < indices_->size (); ++i, j+=dimension)
-    {
-      points[j + 0] = static_cast<coordT> (input_->points[(*indices_)[i]].x);
-      points[j + 1] = static_cast<coordT> (input_->points[(*indices_)[i]].y);
-    }
-  } 
-  else if (yz_proj_safe)
-  {
-    for (size_t i = 0; i < indices_->size (); ++i, j+=dimension)
-    {
-      points[j + 0] = static_cast<coordT> (input_->points[(*indices_)[i]].y);
-      points[j + 1] = static_cast<coordT> (input_->points[(*indices_)[i]].z);
-    }
-  }
-  else if (xz_proj_safe)
-  {
-    for (size_t i = 0; i < indices_->size (); ++i, j+=dimension)
-    {
-      points[j + 0] = static_cast<coordT> (input_->points[(*indices_)[i]].x);
-      points[j + 1] = static_cast<coordT> (input_->points[(*indices_)[i]].z);
-    }
-  }
-  else
+  if (!xy_proj_safe && !yz_proj_safe && !xz_proj_safe)
   {
     // This should only happen if we had invalid input
     PCL_ERROR ("[pcl::%s::performReconstruction2D] Invalid input!\n", getClassName ().c_str ());
+    return;
+  }
+
+  // works because all _proj_safe are mutually exclusive
+  const uint8_t c0 = yz_proj_safe;
+  const uint8_t c1 = !xy_proj_safe + 1;
+
+  for (size_t i = 0, j = 0; i < indices_->size ();++i)
+  {
+    points[j++].input_->points[(*indices_)[i]].data[c0];
+    points[j++].input_->points[(*indices_)[i]].data[c1];
   }
-   
+
+  orgQhull::Qhull qhull;
   // Compute convex hull
-  int exitcode = qh_new_qhull (dimension, static_cast<int> (indices_->size ()), points, ismalloc, const_cast<char*> (flags), outfile, errfile);
-#ifdef HAVE_QHULL_2011
-  if (compute_area_)
+  try
   {
-    qh_prepare_output();
+    qhull.runQhull("", dimension, indices_->size (), points.data(), "");
+
+    // 0 if no error from qhull or it doesn't find any vertices
+    if (qhull.qhullStatus() != 0 || qhull.vertexCount() == 0)
+    {
+      throw std::runtime_error(qhull.qhullMessage());
+    }
   }
-#endif
-    
-  // 0 if no error from qhull or it doesn't find any vertices
-  if (exitcode != 0 || qh num_vertices == 0)
+  catch (const std::exception& e)
   {
-    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 ());
+    PCL_ERROR ("[pcl::%s::performReconstrution2D] ERROR: qhull was unable to compute a convex hull for the given point cloud (%lu)!\nError: %s\n", getClassName ().c_str (), indices_->size (), e.what());
 
     hull.points.resize (0);
     hull.width = hull.height = 0;
     polygons.resize (0);
 
-    qh_freeqhull (!qh_ALL);
-    int curlong, totlong;
-    qh_memfreeshort (&curlong, &totlong);
-
     return;
   }
 
   // Qhull returns the area in volume for 2D
   if (compute_area_)
   {
-    total_area_ = qh totvol;
+    total_area_ = qhull.volume();
     total_volume_ = 0.0;
   }
 
-  int num_vertices = qh num_vertices;
-  hull.points.resize (num_vertices);
-  memset (&hull.points[0], static_cast<int> (hull.points.size ()), sizeof (PointInT));
-
-  vertexT * vertex;
-  int i = 0;
-
-  std::vector<std::pair<int, Eigen::Vector4f>, Eigen::aligned_allocator<std::pair<int, Eigen::Vector4f> > > idx_points (num_vertices);
-  idx_points.resize (hull.points.size ());
-  memset (&idx_points[0], static_cast<int> (hull.points.size ()), sizeof (std::pair<int, Eigen::Vector4f>));
+  orgQhull::QhullVertexList vertexs = qhull.vertexList();
 
-  FORALLvertices
+  std::vector<std::pair<int, Eigen::Vector4f>, Eigen::aligned_allocator<std::pair<int, Eigen::Vector4f> > > idx_points (vertexs.size ());
+  for(orgQhull::QhullVertexList::iterator it = vertexs.begin (); it != vertexs.end (); ++it)
   {
-    hull.points[i] = input_->points[(*indices_)[qh_pointid (vertex->point)]];
-    idx_points[i].first = qh_pointid (vertex->point);
-    ++i;
+    hull.push_back(input_->points[(*indices_)[it->point().id ()]]);
+    idx_points.push_back(std::make_pair (it->point ().id (), Eigen::Vector4f ()));
   }
 
   // Sort
   Eigen::Vector4f centroid;
   pcl::compute3DCentroid (hull, centroid);
-  if (xy_proj_safe)
+
+  for (size_t j = 0; j < hull.points.size (); ++j)
   {
-    for (size_t j = 0; j < hull.points.size (); j++)
-    {
-      idx_points[j].second[0] = hull.points[j].x - centroid[0];
-      idx_points[j].second[1] = hull.points[j].y - centroid[1];
-    }
-  }
-  else if (yz_proj_safe)
-  {
-    for (size_t j = 0; j < hull.points.size (); j++)
-    {
-      idx_points[j].second[0] = hull.points[j].y - centroid[1];
-      idx_points[j].second[1] = hull.points[j].z - centroid[2];
-    }
-  }
-  else if (xz_proj_safe)
-  {
-    for (size_t j = 0; j < hull.points.size (); j++)
-    {
-      idx_points[j].second[0] = hull.points[j].x - centroid[0];
-      idx_points[j].second[1] = hull.points[j].z - centroid[2];
-    }
+    idx_points[j].second[0] = hull.points[j][c0] - centroid[c0];
+    idx_points[j].second[1] = hull.points[j][c1] - centroid[c1];
   }
+
   std::sort (idx_points.begin (), idx_points.end (), comparePoints2D);
-    
+
   polygons.resize (1);
   polygons[0].vertices.resize (hull.points.size ());
 
@@ -267,21 +209,16 @@ pcl::ConvexHull<PointInT>::performReconstruction2D (PointCloud &hull, std::vecto
   hull_indices_.indices.clear ();
   hull_indices_.indices.reserve (hull.points.size ());
 
-  for (int j = 0; j < static_cast<int> (hull.points.size ()); j++)
+  for (unsigned int j = 0; j < hull.points.size (); ++j)
   {
     hull_indices_.indices.push_back ((*indices_)[idx_points[j].first]);
     hull.points[j] = input_->points[(*indices_)[idx_points[j].first]];
-    polygons[0].vertices[j] = static_cast<unsigned int> (j);
+    polygons[0].vertices[j] = j;
   }
-    
-  qh_freeqhull (!qh_ALL);
-  int curlong, totlong;
-  qh_memfreeshort (&curlong, &totlong);
 
   hull.width = static_cast<uint32_t> (hull.points.size ());
   hull.height = 1;
   hull.is_dense = false;
-  return;
 }
 
 #ifdef __GNUC__
@@ -294,119 +231,90 @@ pcl::ConvexHull<PointInT>::performReconstruction3D (
 {
   int dimension = 3;
 
-  // True if qhull should free points in qh_freeqhull() or reallocation
-  boolT ismalloc = True;
-  // output from qh_produce_output(), use NULL to skip qh_produce_output()
-  FILE *outfile = NULL;
-
-#ifndef HAVE_QHULL_2011
-  if (compute_area_)
-    outfile = stderr;
-#endif
-
-  // option flags for qhull, see qh_opt.htm
-  const char *flags = qhull_flags.c_str ();
-  // error messages from qhull code
-  FILE *errfile = stderr;
-
   // Array of coordinates for each point
-  coordT *points = reinterpret_cast<coordT*> (calloc (indices_->size () * dimension, sizeof (coordT)));
+  std::vector<double> points (3 * indices_->size ());
 
-  int j = 0;
-  for (size_t i = 0; i < indices_->size (); ++i, j+=dimension)
+  for (size_t i = 0; i < indices_->size (); ++i)
   {
-    points[j + 0] = static_cast<coordT> (input_->points[(*indices_)[i]].x);
-    points[j + 1] = static_cast<coordT> (input_->points[(*indices_)[i]].y);
-    points[j + 2] = static_cast<coordT> (input_->points[(*indices_)[i]].z);
+    points[3*i].push_back(input_->points[(*indices_)[i]].x);
+    points[3*i + 1].push_back(input_->points[(*indices_)[i]].y);
+    points[3*i + 2].push_back(input_->points[(*indices_)[i]].z);
   }
 
+  orgQhull::Qhull qhull;
   // Compute convex hull
-  int exitcode = qh_new_qhull (dimension, static_cast<int> (indices_->size ()), points, ismalloc, const_cast<char*> (flags), outfile, errfile);
-#ifdef HAVE_QHULL_2011
-  if (compute_area_)
-  {
-    qh_prepare_output();
-  }
-#endif
+  try {
+    qhull.runQhull("", dimension, indices_->size (), points.data(), "");
+
+    // 0 if no error from qhull or it doesn't find any vertices
+    if (qhull.qhullStatus() != 0 || qhull.vertexCount() == 0)
+    {
+      throw std::runtime_error(qhull.qhullMessage());
+    }
 
-  // 0 if no error from qhull
-  if (exitcode != 0)
+    qh_triangulate (qhull.qh());
+  }
+  catch (const std::exception& e)
   {
-    PCL_ERROR ("[pcl::%s::performReconstrution3D] ERROR: qhull was unable to compute a convex hull for the given point cloud (%lu)!\n", getClassName ().c_str (), input_->points.size ());
+    PCL_ERROR ("[pcl::%s::performReconstrution3D] ERROR: qhull was unable to compute a convex hull for the given point cloud (%lu)!\nError: %s\n", getClassName ().c_str (), input_->points.size (), e.what());
 
     hull.points.resize (0);
     hull.width = hull.height = 0;
     polygons.resize (0);
-
-    qh_freeqhull (!qh_ALL);
-    int curlong, totlong;
-    qh_memfreeshort (&curlong, &totlong);
-
-    return;
   }
 
-  qh_triangulate ();
-
-  int num_facets = qh num_facets;
-
-  int num_vertices = qh num_vertices;
+  int num_facets = qhull.facetCount();
+  int num_vertices = qhull.vertexCount();
   hull.points.resize (num_vertices);
 
-  vertexT * vertex;
-  int i = 0;
   // Max vertex id
   unsigned int max_vertex_id = 0;
-  FORALLvertices
+  orgQhull::QhullVertexList vertexs = qhull.vertexList();
+  for(orgQhull::QhullVertexList::iterator it = vertexs.begin (); it != vertexs.end (); ++it)
   {
-    if (vertex->id + 1 > max_vertex_id)
-      max_vertex_id = vertex->id + 1;
+    if (it->id() + 1 > max_vertex_id)
+      max_vertex_id = it->id() + 1;
   }
 
-  ++max_vertex_id;
-  std::vector<int> qhid_to_pcidx (max_vertex_id);
+  std::vector<uint32_t> qhid_to_pcidx (++max_vertex_id);
 
   hull_indices_.header = input_->header;
   hull_indices_.indices.clear ();
   hull_indices_.indices.reserve (num_vertices);
 
-  FORALLvertices
+  uint32_t i = 0;
+  for(orgQhull::QhullVertexList::iterator it = vertexs.begin (); it != vertexs.end (); ++it, ++i)
   {
     // Add vertices to hull point_cloud and store index
-    hull_indices_.indices.push_back ((*indices_)[qh_pointid (vertex->point)]);
-    hull.points[i] = input_->points[(*indices_)[hull_indices_.indices.back ()]];
+    hull_indices_.indices.push_back ((*indices_)[it->point().id()]);
+    hull.points[i] = input_->points[hull_indices_.indices.back ()];
 
-    qhid_to_pcidx[vertex->id] = i; // map the vertex id of qhull to the point cloud index
-    ++i;
+    qhid_to_pcidx[it->id()] = i; // map the vertex id of qhull to the point cloud index
   }
 
   if (compute_area_)
   {
-    total_area_  = qh totarea;
-    total_volume_ = qh totvol;
+    total_area_  = qhull.area();
+    total_volume_ = qhull.volume();
   }
 
   if (fill_polygon_data)
   {
     polygons.resize (num_facets);
-    int dd = 0;
+    size_t dd = 0;
 
-    facetT * facet;
-    FORALLfacets
+    orgQhull::QhullFacetList facets = qhull.facetList();
+    for(orgQhull::QhullFacetList::iterator it = facets.begin (); it != facets.end (); ++it, ++dd)
     {
       polygons[dd].vertices.resize (3);
 
-      // Needed by FOREACHvertex_i_
-      int vertex_n, vertex_i;
-      FOREACHvertex_i_ ((*facet).vertices)
-      //facet_vertices.vertices.push_back (qhid_to_pcidx[vertex->id]);
-      polygons[dd].vertices[vertex_i] = qhid_to_pcidx[vertex->id];
-      ++dd;
+      orgQhull::QhullVertexSet vertices = it->vertices();
+      size_t vertex_i = 0;
+      for(orgQhull::QhullVertexSet::iterator jt = vertices.begin (); jt != vertices.end (); ++jt, ++vertex_i) {
+        polygons[dd].vertices[vertex_i] = qhid_to_pcidx[(*jt).id()];
+      }
     }
   }
-  // Deallocates memory (also the points)
-  qh_freeqhull (!qh_ALL);
-  int curlong, totlong;
-  qh_memfreeshort (&curlong, &totlong);
 
   hull.width = static_cast<uint32_t> (hull.points.size ());
   hull.height = 1;
diff --git a/surface/include/pcl/surface/qhull.h b/surface/include/pcl/surface/qhull.h
index ac738991e71..f5293f8e561 100644
--- a/surface/include/pcl/surface/qhull.h
+++ b/surface/include/pcl/surface/qhull.h
@@ -43,32 +43,11 @@
 #ifndef PCL_QHULL_H
 #define PCL_QHULL_H
 
-#if defined __GNUC__
-#  pragma GCC system_header 
-#endif
-
-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"
-#else
-#  include "qhull/qhull.h"
-#  include "qhull/mem.h"
-#  include "qhull/qset.h"
-#  include "qhull/geom.h"
-#  include "qhull/merge.h"
-#  include "qhull/poly.h"
-#  include "qhull/io.h"
-#  include "qhull/stat.h"
-#endif
-}
+#include <libqhullcpp/Qhull.h>
+#include <libqhullcpp/QhullFacetList.h>
+#include <libqhullcpp/QhullVertexSet.h>
+#include <libqhullcpp/QhullRidge.h>
+#include <libqhullcpp/QhullError.h>
 
 #endif      // PCL_QHULL_H
 #endif