Clustering of Pointclouds into Supervoxels - Theoretical primer

In this tutorial, we show how to divide a pointcloud into a number of supervoxel clusters using pcl::SupervoxelClustering, and then how to use and visualize the adjacency information and supervoxels themselves.

_images/supervoxel_clustering_example.jpg

An example of supervoxels and adjacency graph generated for a cloud

Segmentation algorithms aim to group pixels in images into perceptually meaningful regions which conform to object boundaries. Graph-based approaches, such as Markov Random Field (MRF) and Conditional Random Field (CRF), have become popular, as they merge relational low-level context within the image with object level class knowledge. The cost of solving pixel-level graphs led to the development of mid-level inference schemes which do not use pixels directly, but rather use groupings of pixels, known as superpixels, as the base level for nodes. Superpixels are formed by over-segmenting the image into small regions based on local low-level features, reducing the number of nodes which must be considered for inference.

Due to their strong impact on the quality of the eventual segmentation, it is important that superpixels have certain characteristics. Of these, avoiding violating object boundaries is the most vital, as failing to do so will decrease the accuracy of classifiers used later - since they will be forced to consider pixels which belong to more than one class. Additionally, even if the classifier does manage a correct output, the final pixel level segmentation will necessarily contain errors. Another useful quality is regular distribution over the area being segmented, as this will produce a simpler graph for later steps.

Voxel Cloud Connectivity Segmentation (VCCS) is a recent “superpixel” method which generates volumetric over-segmentations of 3D point cloud data, known as supervoxels. Supervoxels adhere to object boundaries better than state-of-the-art 2D methods, while remaining efficient enough to use in online applications. VCCS uses a region growing variant of k-means clustering for generating its labeling of points directly within a voxel octree structure. Supervoxels have two important properties; they are evenly distributed across the 3D space, and they cannot cross boundaries unless the underlying voxels are spatial connected. The former is accomplished by seeding supervoxels directly in the cloud, rather than the projected plane, while the latter uses an octree structure which maintains adjacency information of leaves. Supervoxels maintain adjacency relations in voxelized 3D space; specifically, 26-adjacency- that is neighboring voxels are those that share a face, edge, or vertex, as seen below.

_images/supervoxel_clustering_adjacency.jpg

From right to left, 6 (faces), 18 (faces,edges), and 26 (faces, edges, vertices) adjacency

The adjacency graph of supervoxels (and the underlying voxels) is maintained efficiently within the octree by specifying that neighbors are voxels within R_voxel of one another, where R_voxel specifies the octree leaf resolution. This adjacency graph is used extensively for both the region growing used to generate the supervoxels, as well as determining adjacency of the resulting supervoxels themselves.

VCCS is a region growing method which incrementally expand supervoxels from a set of seed points distributed evenly in space on a grid with resolution R_seed. To maintain efficiency, VCCS does not search globally, but rather only considers points within R_seed of the seed center. Additionally, seeds which are isolated are filtered out by establishing a small search radius R_search around each seed and removing seeds which do not have sufficient neighbor voxels connected to them.

_images/supervoxel_clustering_parameters.jpg

The various sizing parameters which affect supervoxel clustering. R_seed and R_voxel must both be set by the user.

Expansion from the seed points is governed by a distance measure calculated in a feature space consisting of spatial extent, color, and normals. The spatial distance D_s is normalized by the seeding resolution, color distance D_c is the euclidean distance in normalized RGB space, and normal distance D_n measures the angle between surface normal vectors.

_images/supervoxel_clustering_distance_eqn.png

Weighting equation used in supervoxel clustering. w_c, w_s, and w_n, the color, spatial, and normal weights, respectively, are user controlled parameters.

Supervoxels are grown iteratively, using a local k-means clustering which considers connectivity and flow. The general process is as follows. Beginning at the voxel nearest the cluster center, we flow outward to adjacent voxels and compute the distance from each of these to the supervoxel center using the distance equation above. If the distance is the smallest this voxel has seen, its label is set, and using the adjacency graph, we add its neighbors which are further from the center to our search queue for this label. We then proceed to the next supervoxel, so that each level outwards from the center is considered at the same time for all supervoxels (a 2d version of this is seen in the figure below). We proceed iteratively outwards until we have reached the edge of the search volume for each supervoxel (or have no more neighbors to check).

_images/supervoxel_clustering_search_order.jpg

Search order in the adjacency octree for supervoxel cluster expansion. Dotted edges in the adjacency graph are not searched, since they have already been considered earlier in the queue.

Alright, let’s get to the code… but if you want further details on how supervoxels work (and if you use them in an academic work) please reference the following publication:

@InProceedings{Papon13CVPR,
  author={Jeremie Papon and Alexey Abramov and Markus Schoeler and Florentin W\"{o}rg\"{o}tter},
  title={Voxel Cloud Connectivity Segmentation - Supervoxels for Point Clouds},
  booktitle={Computer Vision and Pattern Recognition (CVPR), 2013 IEEE Conference on},
  month     = {June 22-27},
  year      = {2013},
  address   = {Portland, Oregon},
}

Oh, and for a more complicated example which uses Supervoxels, see pcl/examples/segmentation/supervoxel_clustering.cpp.

The code

First, grab a pcd file made from a kinect or similar device - here we shall use milk_cartoon_all_small_clorox.pcd which is available in the pcl git here). Next, copy and paste the following code into your editor and save it as supervoxel_clustering.cpp (or download the source file here).

  1#include <pcl/console/parse.h>
  2#include <pcl/point_cloud.h>
  3#include <pcl/point_types.h>
  4#include <pcl/io/pcd_io.h>
  5#include <pcl/visualization/pcl_visualizer.h>
  6#include <pcl/segmentation/supervoxel_clustering.h>
  7
  8//VTK include needed for drawing graph lines
  9#include <vtkPolyLine.h>
 10
 11// Types
 12typedef pcl::PointXYZRGBA PointT;
 13typedef pcl::PointCloud<PointT> PointCloudT;
 14typedef pcl::PointNormal PointNT;
 15typedef pcl::PointCloud<PointNT> PointNCloudT;
 16typedef pcl::PointXYZL PointLT;
 17typedef pcl::PointCloud<PointLT> PointLCloudT;
 18
 19void addSupervoxelConnectionsToViewer (PointT &supervoxel_center,
 20                                       PointCloudT &adjacent_supervoxel_centers,
 21                                       std::string supervoxel_name,
 22                                       pcl::visualization::PCLVisualizer::Ptr & viewer);
 23
 24
 25int
 26main (int argc, char ** argv)
 27{
 28  if (argc < 2)
 29  {
 30    pcl::console::print_error ("Syntax is: %s <pcd-file> \n "
 31                                "--NT Dsables the single cloud transform \n"
 32                                "-v <voxel resolution>\n-s <seed resolution>\n"
 33                                "-c <color weight> \n-z <spatial weight> \n"
 34                                "-n <normal_weight>\n", argv[0]);
 35    return (1);
 36  }
 37
 38
 39  PointCloudT::Ptr cloud (new PointCloudT);
 40  pcl::console::print_highlight ("Loading point cloud...\n");
 41  if (pcl::io::loadPCDFile<PointT> (argv[1], *cloud))
 42  {
 43    pcl::console::print_error ("Error loading cloud file!\n");
 44    return (1);
 45  }
 46
 47
 48  bool disable_transform = pcl::console::find_switch (argc, argv, "--NT");
 49
 50  float voxel_resolution = 0.008f;
 51  bool voxel_res_specified = pcl::console::find_switch (argc, argv, "-v");
 52  if (voxel_res_specified)
 53    pcl::console::parse (argc, argv, "-v", voxel_resolution);
 54
 55  float seed_resolution = 0.1f;
 56  bool seed_res_specified = pcl::console::find_switch (argc, argv, "-s");
 57  if (seed_res_specified)
 58    pcl::console::parse (argc, argv, "-s", seed_resolution);
 59
 60  float color_importance = 0.2f;
 61  if (pcl::console::find_switch (argc, argv, "-c"))
 62    pcl::console::parse (argc, argv, "-c", color_importance);
 63
 64  float spatial_importance = 0.4f;
 65  if (pcl::console::find_switch (argc, argv, "-z"))
 66    pcl::console::parse (argc, argv, "-z", spatial_importance);
 67
 68  float normal_importance = 1.0f;
 69  if (pcl::console::find_switch (argc, argv, "-n"))
 70    pcl::console::parse (argc, argv, "-n", normal_importance);
 71
 72  //////////////////////////////  //////////////////////////////
 73  ////// This is how to use supervoxels
 74  //////////////////////////////  //////////////////////////////
 75
 76  pcl::SupervoxelClustering<PointT> super (voxel_resolution, seed_resolution);
 77  if (disable_transform)
 78    super.setUseSingleCameraTransform (false);
 79  super.setInputCloud (cloud);
 80  super.setColorImportance (color_importance);
 81  super.setSpatialImportance (spatial_importance);
 82  super.setNormalImportance (normal_importance);
 83
 84  std::map <std::uint32_t, pcl::Supervoxel<PointT>::Ptr > supervoxel_clusters;
 85
 86  pcl::console::print_highlight ("Extracting supervoxels!\n");
 87  super.extract (supervoxel_clusters);
 88  pcl::console::print_info ("Found %d supervoxels\n", supervoxel_clusters.size ());
 89
 90  pcl::visualization::PCLVisualizer::Ptr viewer (new pcl::visualization::PCLVisualizer ("3D Viewer"));
 91  viewer->setBackgroundColor (0, 0, 0);
 92
 93  PointCloudT::Ptr voxel_centroid_cloud = super.getVoxelCentroidCloud ();
 94  viewer->addPointCloud (voxel_centroid_cloud, "voxel centroids");
 95  viewer->setPointCloudRenderingProperties (pcl::visualization::PCL_VISUALIZER_POINT_SIZE,2.0, "voxel centroids");
 96  viewer->setPointCloudRenderingProperties (pcl::visualization::PCL_VISUALIZER_OPACITY,0.95, "voxel centroids");
 97
 98  PointLCloudT::Ptr labeled_voxel_cloud = super.getLabeledVoxelCloud ();
 99  viewer->addPointCloud (labeled_voxel_cloud, "labeled voxels");
100  viewer->setPointCloudRenderingProperties (pcl::visualization::PCL_VISUALIZER_OPACITY,0.8, "labeled voxels");
101
102  PointNCloudT::Ptr sv_normal_cloud = super.makeSupervoxelNormalCloud (supervoxel_clusters);
103  //We have this disabled so graph is easy to see, uncomment to see supervoxel normals
104  //viewer->addPointCloudNormals<PointNormal> (sv_normal_cloud,1,0.05f, "supervoxel_normals");
105
106  pcl::console::print_highlight ("Getting supervoxel adjacency\n");
107  std::multimap<std::uint32_t, std::uint32_t> supervoxel_adjacency;
108  super.getSupervoxelAdjacency (supervoxel_adjacency);
109  //To make a graph of the supervoxel adjacency, we need to iterate through the supervoxel adjacency multimap
110  for (auto label_itr = supervoxel_adjacency.cbegin (); label_itr != supervoxel_adjacency.cend (); )
111  {
112    //First get the label
113    std::uint32_t supervoxel_label = label_itr->first;
114    //Now get the supervoxel corresponding to the label
115    pcl::Supervoxel<PointT>::Ptr supervoxel = supervoxel_clusters.at (supervoxel_label);
116
117    //Now we need to iterate through the adjacent supervoxels and make a point cloud of them
118    PointCloudT adjacent_supervoxel_centers;
119    for (auto adjacent_itr = supervoxel_adjacency.equal_range (supervoxel_label).first; adjacent_itr!=supervoxel_adjacency.equal_range (supervoxel_label).second; ++adjacent_itr)
120    {
121      pcl::Supervoxel<PointT>::Ptr neighbor_supervoxel = supervoxel_clusters.at (adjacent_itr->second);
122      adjacent_supervoxel_centers.push_back (neighbor_supervoxel->centroid_);
123    }
124    //Now we make a name for this polygon
125    std::stringstream ss;
126    ss << "supervoxel_" << supervoxel_label;
127    //This function is shown below, but is beyond the scope of this tutorial - basically it just generates a "star" polygon mesh from the points given
128    addSupervoxelConnectionsToViewer (supervoxel->centroid_, adjacent_supervoxel_centers, ss.str (), viewer);
129    //Move iterator forward to next label
130    label_itr = supervoxel_adjacency.upper_bound (supervoxel_label);
131  }
132
133  while (!viewer->wasStopped ())
134  {
135    viewer->spinOnce (100);
136  }
137  return (0);
138}
139
140void
141addSupervoxelConnectionsToViewer (PointT &supervoxel_center,
142                                  PointCloudT &adjacent_supervoxel_centers,
143                                  std::string supervoxel_name,
144                                  pcl::visualization::PCLVisualizer::Ptr & viewer)
145{
146  vtkSmartPointer<vtkPoints> points = vtkSmartPointer<vtkPoints>::New ();
147  vtkSmartPointer<vtkCellArray> cells = vtkSmartPointer<vtkCellArray>::New ();
148  vtkSmartPointer<vtkPolyLine> polyLine = vtkSmartPointer<vtkPolyLine>::New ();
149
150  //Iterate through all adjacent points, and add a center point to adjacent point pair
151  for (auto adjacent_itr = adjacent_supervoxel_centers.begin (); adjacent_itr != adjacent_supervoxel_centers.end (); ++adjacent_itr)
152  {
153    points->InsertNextPoint (supervoxel_center.data);
154    points->InsertNextPoint (adjacent_itr->data);
155  }
156  // Create a polydata to store everything in
157  vtkSmartPointer<vtkPolyData> polyData = vtkSmartPointer<vtkPolyData>::New ();
158  // Add the points to the dataset
159  polyData->SetPoints (points);
160  polyLine->GetPointIds  ()->SetNumberOfIds(points->GetNumberOfPoints ());
161  for(unsigned int i = 0; i < points->GetNumberOfPoints (); i++)
162    polyLine->GetPointIds ()->SetId (i,i);
163  cells->InsertNextCell (polyLine);
164  // Add the lines to the dataset
165  polyData->SetLines (cells);
166  viewer->addModelFromPolyData (polyData,supervoxel_name);
167}

The explanation

We start by defining convenience types in order not to clutter the code.

typedef pcl::PointXYZRGBA PointT;
typedef pcl::PointCloud<PointT> PointCloudT;
typedef pcl::PointNormal PointNT;
typedef pcl::PointCloud<PointNT> PointNCloudT;
typedef pcl::PointXYZL PointLT;
typedef pcl::PointCloud<PointLT> PointLCloudT;

Then we load the input cloud based on the input argument

  PointCloudT::Ptr cloud (new PointCloudT);
  pcl::console::print_highlight ("Loading point cloud...\n");
  if (pcl::io::loadPCDFile<PointT> (argv[1], *cloud))
  {
    pcl::console::print_error ("Error loading cloud file!\n");
    return (1);
  }

Next we check the input arguments and set default values. You can play with the various parameters to see how they affect the supervoxels, but briefly:

  • --NT Disables the single-view transform (this is the default for unorganized clouds, only affects organized clouds)

  • -v Sets the voxel size, which determines the leaf size of the underlying octree structure (in meters)

  • -s Sets the seeding size, which determines how big the supervoxels will be (in meters)

  • -c Sets the weight for color - how much color will influence the shape of the supervoxels

  • -z Sets the weight for spatial term - higher values will result in supervoxels with very regular shapes (lower will result in supervoxels which follow normals and/or colors, but are not very regular)

  • -n Sets the weight for normal - how much surface normals will influence the shape of the supervoxels

  bool disable_transform = pcl::console::find_switch (argc, argv, "--NT");

  float voxel_resolution = 0.008f;
  bool voxel_res_specified = pcl::console::find_switch (argc, argv, "-v");
  if (voxel_res_specified)
    pcl::console::parse (argc, argv, "-v", voxel_resolution);

  float seed_resolution = 0.1f;
  bool seed_res_specified = pcl::console::find_switch (argc, argv, "-s");
  if (seed_res_specified)
    pcl::console::parse (argc, argv, "-s", seed_resolution);

  float color_importance = 0.2f;
  if (pcl::console::find_switch (argc, argv, "-c"))
    pcl::console::parse (argc, argv, "-c", color_importance);

  float spatial_importance = 0.4f;
  if (pcl::console::find_switch (argc, argv, "-z"))
    pcl::console::parse (argc, argv, "-z", spatial_importance);

  float normal_importance = 1.0f;
  if (pcl::console::find_switch (argc, argv, "-n"))
    pcl::console::parse (argc, argv, "-n", normal_importance);

We are now ready to setup the supervoxel clustering. We use the class SupervoxelClustering, which implements the clustering process and give it the parameters.

Important

By default, the algorithm will use a special transform compressing the depth in Z if your input cloud is organized (eg, from an RGBD sensor like the Kinect). You MUST set use_transform to false if you are using an organized cloud which doesn’t have the camera at (0,0,0) and depth in positive Z. The transform is specifically designed to help improve Kinect data by increasing voxel bin size as distance from the camera increases. If your cloud is unorganized, this transform will not be used by default, but can be enabled by using setUseSingleCameraTransform(true).

  pcl::SupervoxelClustering<PointT> super (voxel_resolution, seed_resolution);
  if (disable_transform)
    super.setUseSingleCameraTransform (false);
  super.setInputCloud (cloud);
  super.setColorImportance (color_importance);
  super.setSpatialImportance (spatial_importance);
  super.setNormalImportance (normal_importance);

Then we initialize the data structure which will be used to extract the supervoxels, and run the algorithm. The data structure is a map from labels to shared pointers of Supervoxel templated on the input point type. Supervoxels have the following fields:

  • normal_ The normal calculated for the voxels contained in the supervoxel

  • centroid_ The centroid of the supervoxel - average voxel

  • voxels_ A Pointcloud of the voxels in the supervoxel

  • normals_ A Pointcloud of the normals for the points in the supervoxel

  std::map <std::uint32_t, pcl::Supervoxel<PointT>::Ptr > supervoxel_clusters;

  pcl::console::print_highlight ("Extracting supervoxels!\n");
  super.extract (supervoxel_clusters);
  pcl::console::print_info ("Found %d supervoxels\n", supervoxel_clusters.size ());

We then load a viewer and use some of the getter functions of SupervoxelClustering to pull out clouds to display. voxel_centroid_cloud contains the voxel centroids coming out of the octree (basically the downsampled original cloud), and colored_voxel_cloud are the voxels colored according to their supervoxel labels (random colors). sv_normal_cloud contains a cloud of the supervoxel normals, but we don’t display it here so that the graph is visible.

  pcl::visualization::PCLVisualizer::Ptr viewer (new pcl::visualization::PCLVisualizer ("3D Viewer"));
  viewer->setBackgroundColor (0, 0, 0);

  PointCloudT::Ptr voxel_centroid_cloud = super.getVoxelCentroidCloud ();
  viewer->addPointCloud (voxel_centroid_cloud, "voxel centroids");
  viewer->setPointCloudRenderingProperties (pcl::visualization::PCL_VISUALIZER_POINT_SIZE,2.0, "voxel centroids");
  viewer->setPointCloudRenderingProperties (pcl::visualization::PCL_VISUALIZER_OPACITY,0.95, "voxel centroids");

  PointLCloudT::Ptr labeled_voxel_cloud = super.getLabeledVoxelCloud ();
  viewer->addPointCloud (labeled_voxel_cloud, "labeled voxels");
  viewer->setPointCloudRenderingProperties (pcl::visualization::PCL_VISUALIZER_OPACITY,0.8, "labeled voxels");

  PointNCloudT::Ptr sv_normal_cloud = super.makeSupervoxelNormalCloud (supervoxel_clusters);
  //We have this disabled so graph is easy to see, uncomment to see supervoxel normals
  //viewer->addPointCloudNormals<PointNormal> (sv_normal_cloud,1,0.05f, "supervoxel_normals");

Finally, we extract the supervoxel adjacency list (in the form of a multimap of label adjacencies).

  pcl::console::print_highlight ("Getting supervoxel adjacency\n");
  std::multimap<std::uint32_t, std::uint32_t> supervoxel_adjacency;
  super.getSupervoxelAdjacency (supervoxel_adjacency);

Then we iterate through the multimap, creating a point cloud of the centroids of each supervoxel’s neighbors.

  for (auto label_itr = supervoxel_adjacency.cbegin (); label_itr != supervoxel_adjacency.cend (); )
  {
    //First get the label
    std::uint32_t supervoxel_label = label_itr->first;
    //Now get the supervoxel corresponding to the label
    pcl::Supervoxel<PointT>::Ptr supervoxel = supervoxel_clusters.at (supervoxel_label);

    //Now we need to iterate through the adjacent supervoxels and make a point cloud of them
    PointCloudT adjacent_supervoxel_centers;
    for (auto adjacent_itr = supervoxel_adjacency.equal_range (supervoxel_label).first; adjacent_itr!=supervoxel_adjacency.equal_range (supervoxel_label).second; ++adjacent_itr)
    {
      pcl::Supervoxel<PointT>::Ptr neighbor_supervoxel = supervoxel_clusters.at (adjacent_itr->second);
      adjacent_supervoxel_centers.push_back (neighbor_supervoxel->centroid_);
    }

Then we create a string label for the supervoxel graph we will draw and call addSupervoxelConnectionsToViewer, a drawing helper function implemented later in the tutorial code. The details of addSupervoxelConnectionsToViewer are beyond the scope of this tutorial, but all it does is draw a star polygon mesh of the supervoxel centroid to all of its neighbors centroids. We need to do this like this because adding individual lines using the addLine functionality of pcl_visualizer is too slow for large numbers of lines.

    //Now we make a name for this polygon
    std::stringstream ss;
    ss << "supervoxel_" << supervoxel_label;
    //This function is shown below, but is beyond the scope of this tutorial - basically it just generates a "star" polygon mesh from the points given
    addSupervoxelConnectionsToViewer (supervoxel->centroid_, adjacent_supervoxel_centers, ss.str (), viewer);
    //Move iterator forward to next label
    label_itr = supervoxel_adjacency.upper_bound (supervoxel_label);

This results in a supervoxel graph that looks like this for seed size of 0.1m (top) and 0.05m (middle). The bottom is the original cloud, given for reference.:

_images/supervoxel_clustering_results.jpg

Compiling and running the program

Create a CMakeLists.txt file with the following content (or download it here):

 1cmake_minimum_required(VERSION 3.5 FATAL_ERROR)
 2
 3project(supervoxel_clustering)
 4
 5find_package(PCL 1.8 REQUIRED)
 6
 7include_directories(${PCL_INCLUDE_DIRS})
 8link_directories(${PCL_LIBRARY_DIRS})
 9add_definitions(${PCL_DEFINITIONS})
10
11add_executable (supervoxel_clustering supervoxel_clustering.cpp)
12target_link_libraries (supervoxel_clustering ${PCL_LIBRARIES})

After you have made the executable, you can run it like so, assuming the pcd file is in the same folder as the executable:

$ ./supervoxel_clustering milk_cartoon_all_small_clorox.pcd --NT

Don’t be afraid to play around with the parameters (especially the seed size, -s) to see what happens. The pcd file name should always be the first parameter!