RoPs (Rotational Projection Statistics) feature

In this tutorial we will learn how to use the pcl::ROPSEstimation class in order to extract points features. The feature extraction method implemented in this class was proposed by Yulan Guo, Ferdous Sohel, Mohammed Bennamoun, Min Lu and Jianwei Wanalso in their article “Rotational Projection Statistics for 3D Local Surface Description and Object Recognition”

Theoretical Primer

The idea of the feature extraction method is as follows. Having a mesh and a set of points for which feature must be computed we perform some simple steps. First of all for a given point of interest the local surface is cropped. Local surface consists of the points and triangles that are within the given support radius. For the given local surface LRF (Local Reference Frame) is computed. LRF is simply a triplet of vectors, the comprehensive information about how these vectors are computed you can find in the article. What is really important is that using these vectors we can provide the invariance to the rotation of the cloud. To do that, we simply translate points of the local surface in such way that point of interest became the origin, after that we rotate local surface so that the LRF vectors were aligned with the Ox, Oy and Oz axes. Having this done, we then start the feature extraction. For every axis Ox, Oy and Oz the following steps are performed, we will refer to these axes as current axis:

  • local surface is rotated around the current axis by a given angle;

  • points of the rotated local surface are projected onto three planes XY, XZ and YZ;

  • for each projection distribution matrix is built, this matrix simply shows how much points fall onto each bin. Number of bins represents the matrix dimension and is the parameter of the algorithm, as well as the support radius;

  • for each distribution matrix central moments are calculated: M11, M12, M21, M22, E. Here E is the Shannon entropy;

  • calculated values are then concatenated to form the sub-feature.

We iterate through these steps several times. Number of iterations depends on the given number of rotations. Sub-features for different axes are concatenated to form the final RoPS descriptor.

The code

For this tutorial we will use the model from the Queen’s Dataset. You can choose any other point cloud, but in order to make the code work you will need to use the triangulation algorithm in order to obtain polygons. You can find the proposed model here:

  • points - contains the point cloud

  • indices - contains indices of the points for which RoPs must be computed

  • triangles - contains the polygons

Next what you need to do is to create a file rops_feature.cpp in any editor you prefer and copy the following code inside of it:

 1#include <pcl/features/rops_estimation.h>
 2#include <pcl/io/pcd_io.h>
 3
 4int main (int argc, char** argv)
 5{
 6  if (argc != 4)
 7    return (-1);
 8
 9  pcl::PointCloud<pcl::PointXYZ>::Ptr cloud (new pcl::PointCloud<pcl::PointXYZ> ());
10  if (pcl::io::loadPCDFile (argv[1], *cloud) == -1)
11    return (-1);
12
13  pcl::PointIndicesPtr indices (new pcl::PointIndices);
14  std::ifstream indices_file;
15  indices_file.open (argv[2], std::ifstream::in);
16  for (std::string line; std::getline (indices_file, line);)
17  {
18    std::istringstream in (line);
19    unsigned int index = 0;
20    in >> index;
21    indices->indices.push_back (index - 1);
22  }
23  indices_file.close ();
24
25  std::vector <pcl::Vertices> triangles;
26  std::ifstream triangles_file;
27  triangles_file.open (argv[3], std::ifstream::in);
28  for (std::string line; std::getline (triangles_file, line);)
29  {
30    pcl::Vertices triangle;
31    std::istringstream in (line);
32    unsigned int vertex = 0;
33    in >> vertex;
34    triangle.vertices.push_back (vertex - 1);
35    in >> vertex;
36    triangle.vertices.push_back (vertex - 1);
37    in >> vertex;
38    triangle.vertices.push_back (vertex - 1);
39    triangles.push_back (triangle);
40  }
41
42  float support_radius = 0.0285f;
43  unsigned int number_of_partition_bins = 5;
44  unsigned int number_of_rotations = 3;
45
46  pcl::search::KdTree<pcl::PointXYZ>::Ptr search_method (new pcl::search::KdTree<pcl::PointXYZ>);
47  search_method->setInputCloud (cloud);
48
49  pcl::ROPSEstimation <pcl::PointXYZ, pcl::Histogram <135> > feature_estimator;
50  feature_estimator.setSearchMethod (search_method);
51  feature_estimator.setSearchSurface (cloud);
52  feature_estimator.setInputCloud (cloud);
53  feature_estimator.setIndices (indices);
54  feature_estimator.setTriangles (triangles);
55  feature_estimator.setRadiusSearch (support_radius);
56  feature_estimator.setNumberOfPartitionBins (number_of_partition_bins);
57  feature_estimator.setNumberOfRotations (number_of_rotations);
58  feature_estimator.setSupportRadius (support_radius);
59
60  pcl::PointCloud<pcl::Histogram <135> >::Ptr histograms (new pcl::PointCloud <pcl::Histogram <135> > ());
61  feature_estimator.compute (*histograms);
62
63  return (0);
64}

The explanation

Now let’s study out what is the purpose of this code.

  pcl::PointCloud<pcl::PointXYZ>::Ptr cloud (new pcl::PointCloud<pcl::PointXYZ> ());
  if (pcl::io::loadPCDFile (argv[1], *cloud) == -1)
    return (-1);

These lines are simply loading the cloud from the .pcd file.

  pcl::PointIndicesPtr indices (new pcl::PointIndices);
  std::ifstream indices_file;
  indices_file.open (argv[2], std::ifstream::in);
  for (std::string line; std::getline (indices_file, line);)
  {
    std::istringstream in (line);
    unsigned int index = 0;
    in >> index;
    indices->indices.push_back (index - 1);
  }
  indices_file.close ();

Here the indices of points for which RoPS feature must be computed are loaded. You can comment it and compute features for every single point in the cloud. if you want.

  std::vector <pcl::Vertices> triangles;
  std::ifstream triangles_file;
  triangles_file.open (argv[3], std::ifstream::in);
  for (std::string line; std::getline (triangles_file, line);)
  {
    pcl::Vertices triangle;
    std::istringstream in (line);
    unsigned int vertex = 0;
    in >> vertex;
    triangle.vertices.push_back (vertex - 1);
    in >> vertex;
    triangle.vertices.push_back (vertex - 1);
    in >> vertex;
    triangle.vertices.push_back (vertex - 1);
    triangles.push_back (triangle);
  }

These lines are loading the information about the polygons. You can replace them with the code for the triangulation if you have only the point cloud instead of the mesh.

  float support_radius = 0.0285f;
  unsigned int number_of_partition_bins = 5;
  unsigned int number_of_rotations = 3;

These code defines important algorithm parameters: support radius for local surface cropping, number of partition bins used to form the distribution matrix and the number of rotations. The last parameter affects the length of the descriptor.

  pcl::search::KdTree<pcl::PointXYZ>::Ptr search_method (new pcl::search::KdTree<pcl::PointXYZ>);
  search_method->setInputCloud (cloud);

These lines set up the search method that will be used by the algorithm.

  pcl::ROPSEstimation <pcl::PointXYZ, pcl::Histogram <135> > feature_estimator;
  feature_estimator.setSearchMethod (search_method);
  feature_estimator.setSearchSurface (cloud);
  feature_estimator.setInputCloud (cloud);
  feature_estimator.setIndices (indices);
  feature_estimator.setTriangles (triangles);
  feature_estimator.setRadiusSearch (support_radius);
  feature_estimator.setNumberOfPartitionBins (number_of_partition_bins);
  feature_estimator.setNumberOfRotations (number_of_rotations);
  feature_estimator.setSupportRadius (support_radius);

Here is the place where the instantiation of the pcl::ROPSEstimation class takes place. It has two parameters:

  • PointInT - type of the input points;

  • PointOutT - type of the output points.

Immediately after that we set the input all the necessary data needed for the feature computation.

  pcl::PointCloud<pcl::Histogram <135> >::Ptr histograms (new pcl::PointCloud <pcl::Histogram <135> > ());
  feature_estimator.compute (*histograms);

Here is the place where the computational process is launched.

Compiling and running the program

Add the following lines to your CMakeLists.txt file:

 1cmake_minimum_required(VERSION 3.5 FATAL_ERROR)
 2
 3project(rops_feature)
 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 (rops_feature rops_feature.cpp)
12target_link_libraries (rops_feature ${PCL_LIBRARIES})

After you have made the executable, you can run it. Simply do:

$ ./rops_feature points.pcd indices.txt triangles.txt