Region growing segmentation

In this tutorial we will learn how to use the region growing algorithm implemented in the pcl::RegionGrowing class. The purpose of the said algorithm is to merge the points that are close enough in terms of the smoothness constraint. Thereby, the output of this algorithm is the set of clusters, where each cluster is a set of points that are considered to be a part of the same smooth surface. The work of this algorithm is based on the comparison of the angles between the points normals.

Theoretical Primer

Let’s take a look at how the algorithm works.

First of all it sorts the points by their curvature value. It needs to be done because the region begins its growth from the point that has the minimum curvature value. The reason for this is that the point with the minimum curvature is located in the flat area (growth from the flattest area allows to reduce the total number of segments).

So we have the sorted cloud. Until there are no unlabeled points in the cloud, the algorithm picks up the point with minimum curvature value and starts the growth of the region. This process occurs as follows:

The picked point is added to the set called seeds.

For every seed point, the algorithm finds its neighbouring points.

Every neighbour is tested for the angle between its normal and normal of the current seed point. If the angle is less than the threshold value then current point is added to the current region.

After that every neighbour is tested for the curvature value. If the curvature is less than the threshold value then this point is added to the seeds.

Current seed is removed from the seeds.

If the seeds set becomes empty this means that the algorithm has grown the region and the process is repeated from the beginning. You can find the pseudocode for the said algorithm below.

Inputs:

Point cloud = $\{P\}$

Point normals = $\{N\}$

Points curvatures = $\{c\}$

Neighbour finding function $\Omega(.)$

Curvature threshold $c_{th}$

Angle threshold $\theta_{th}$

Initialize:

Region list ${R}\leftarrow{\O}$

Available points list $\{A\}\leftarrow\{1,...,|P|\}$

Algorithm:

While $\{A\}$ is not empty do

Current region $\{R_c\}\leftarrow{\O}$

Current seeds $\{S_c\}\leftarrow{\O}$

Point with minimum curvature in $\{A\}\rightarrow P_{min}$

$\{S_c\}\leftarrow\{S_c\}\cup P_{min}$

$\{R_c\}\leftarrow\{R_c\}\cup P_{min}$

$\{A\}\leftarrow\{A\}\setminus P_{min}$

for $i=0$ to size ( $\{S_c\}$ ) do

Find nearest neighbours of current seed point $\{B_c\}\leftarrow\Omega(S_c\{i\})$

for $j=0$ to size ( $\{B_c\}$ ) do

Current neighbour point $P_j\leftarrow B_c\{j\}$

If $\{A\}$ contains $P_j$ and $cos^{-1}(|(N\{S_c\{i\}\},N\{S_c\{j\}\})|)<\theta_{th}$ then

$\{R_c\}\leftarrow\{R_c\}\cup P_j$

$\{A\}\leftarrow\{A\}\setminus P_j$

If $c\{P_j\}<c_{th}$ then

$\{S_c\}\leftarrow\{S_c\}\cup P_j$

end if

end if

end for

end for

Add current region to global segment list $\{R\}\leftarrow\{R\}\cup\{R_c\}$

end while

Return $\{R\}$

The code

First of all you will need the point cloud for this tutorial. This is a good one for the purposes of the algorithm. Next what you need to do is to create a file region_growing_segmentation.cpp in any editor you prefer and copy the following code inside of it:

#include <iostream>
#include <vector>
#include <pcl/point_types.h>
#include <pcl/io/pcd_io.h>
#include <pcl/search/search.h>
#include <pcl/search/kdtree.h>
#include <pcl/features/normal_3d.h>
#include <pcl/visualization/cloud_viewer.h>
#include <pcl/filters/filter_indices.h> // for pcl::removeNaNFromPointCloud
#include <pcl/segmentation/region_growing.h>

int
main ()
{
  pcl::PointCloud<pcl::PointXYZ>::Ptr cloud (new pcl::PointCloud<pcl::PointXYZ>);
  if ( pcl::io::loadPCDFile <pcl::PointXYZ> ("region_growing_tutorial.pcd", *cloud) == -1)
  {
    std::cout << "Cloud reading failed." << std::endl;
    return (-1);
  }

  pcl::search::Search<pcl::PointXYZ>::Ptr tree (new pcl::search::KdTree<pcl::PointXYZ>);
  pcl::PointCloud <pcl::Normal>::Ptr normals (new pcl::PointCloud <pcl::Normal>);
  pcl::NormalEstimation<pcl::PointXYZ, pcl::Normal> normal_estimator;
  normal_estimator.setSearchMethod (tree);
  normal_estimator.setInputCloud (cloud);
  normal_estimator.setKSearch (50);
  normal_estimator.compute (*normals);

  pcl::IndicesPtr indices (new std::vector <int>);
  pcl::removeNaNFromPointCloud(*cloud, *indices);

  pcl::RegionGrowing<pcl::PointXYZ, pcl::Normal> reg;
  reg.setMinClusterSize (50);
  reg.setMaxClusterSize (1000000);
  reg.setSearchMethod (tree);
  reg.setNumberOfNeighbours (30);
  reg.setInputCloud (cloud);
  reg.setIndices (indices);
  reg.setInputNormals (normals);
  reg.setSmoothnessThreshold (3.0 / 180.0 * M_PI);
  reg.setCurvatureThreshold (1.0);

  std::vector <pcl::PointIndices> clusters;
  reg.extract (clusters);

  std::cout << "Number of clusters is equal to " << clusters.size () << std::endl;
  std::cout << "First cluster has " << clusters[0].indices.size () << " points." << std::endl;
  std::cout << "These are the indices of the points of the initial" <<
    std::endl << "cloud that belong to the first cluster:" << std::endl;
  std::size_t counter = 0;
  while (counter < clusters[0].indices.size ())
  {
    std::cout << clusters[0].indices[counter] << ", ";
    counter++;
    if (counter % 10 == 0)
      std::cout << std::endl;
  }
  std::cout << std::endl;

  pcl::PointCloud <pcl::PointXYZRGB>::Ptr colored_cloud = reg.getColoredCloud ();
  pcl::visualization::CloudViewer viewer ("Cluster viewer");
  viewer.showCloud(colored_cloud);
  while (!viewer.wasStopped ())
  {
  }

  return (0);
}

The explanation

Now let’s study out what is the purpose of this code. First few lines will be omitted, because they are obvious.

First lines that are of interest are these:

  pcl::PointCloud<pcl::PointXYZ>::Ptr cloud (new pcl::PointCloud<pcl::PointXYZ>);
  if ( pcl::io::loadPCDFile <pcl::PointXYZ> ("region_growing_tutorial.pcd", *cloud) == -1)
  {
    std::cout << "Cloud reading failed." << std::endl;
    return (-1);
  }

They are simply loading the cloud from the .pcd file. No doubt that you saw how it is done hundreds of times, so let’s move on.

  pcl::search::Search<pcl::PointXYZ>::Ptr tree (new pcl::search::KdTree<pcl::PointXYZ>);
  pcl::PointCloud <pcl::Normal>::Ptr normals (new pcl::PointCloud <pcl::Normal>);
  pcl::NormalEstimation<pcl::PointXYZ, pcl::Normal> normal_estimator;
  normal_estimator.setSearchMethod (tree);
  normal_estimator.setInputCloud (cloud);
  normal_estimator.setKSearch (50);
  normal_estimator.compute (*normals);

As mentioned before, the algorithm requires normals. Here the pcl::NormalEstimation class is used to compute them. To learn more about how it is done you should take a look at the Estimating Surface Normals in a PointCloud tutorial in the Features section.

  pcl::IndicesPtr indices (new std::vector <int>);
  pcl::removeNaNFromPointCloud(*cloud, *indices);

Insofar as pcl::RegionGrowing is derived from pcl::PCLBase, it can work with indices. It means you can instruct that you only segment those points that are listed in the indices array instead of the whole point cloud. Here, only the valid points are chosen for segmentation.

  pcl::RegionGrowing<pcl::PointXYZ, pcl::Normal> reg;
  reg.setMinClusterSize (50);
  reg.setMaxClusterSize (1000000);

You have finally reached the part where pcl::RegionGrowing is instantiated. It is a template class that has two parameters:

PointT - type of points to use(in the given example it is pcl::PointXYZ)
NormalT - type of normals to use(in the given example it is pcl::Normal)

After that minimum and maximum cluster sizes are set. It means that after the segmentation is done all clusters that have less points than minimum(or have more than maximum) will be discarded. The default values for minimum and maximum are 1 and ‘as much as possible’ respectively.

  reg.setSearchMethod (tree);
  reg.setNumberOfNeighbours (30);
  reg.setInputCloud (cloud);
  reg.setIndices (indices);
  reg.setInputNormals (normals);

The algorithm needs K nearest search in its internal structure, so here is the place where a search method is provided and number of neighbours is set. After that it receives the cloud that must be segmented, point indices and normals.

  reg.setSmoothnessThreshold (3.0 / 180.0 * M_PI);
  reg.setCurvatureThreshold (1.0);

These two lines are the most important part in the algorithm initialization, because they are responsible for the mentioned smoothness constraint. First method sets the angle in radians that will be used as the allowable range for the normals deviation. If the deviation between points normals is less than the smoothness threshold then they are suggested to be in the same cluster (new point - the tested one - will be added to the cluster). The second one is responsible for the curvature threshold. If two points have a small normals deviation then the disparity between their curvatures is tested. And if this value is less than the curvature threshold then the algorithm will continue the growth of the cluster using the newly added point.

  std::vector <pcl::PointIndices> clusters;
  reg.extract (clusters);

This method simply launches the segmentation algorithm. After its work it will return clusters array.

  std::cout << "Number of clusters is equal to " << clusters.size () << std::endl;
  std::cout << "First cluster has " << clusters[0].indices.size () << " points." << std::endl;
  std::cout << "These are the indices of the points of the initial" <<
    std::endl << "cloud that belong to the first cluster:" << std::endl;
  std::size_t counter = 0;
  while (counter < clusters[0].indices.size ())
  {
    std::cout << clusters[0].indices[counter] << ", ";
    counter++;
    if (counter % 10 == 0)
      std::cout << std::endl;
  }
  std::cout << std::endl;

These lines are simple enough, so they won’t be commented. They are intended for those who are not familiar with how to work with pcl::PointIndices and how to access its elements.

  pcl::PointCloud <pcl::PointXYZRGB>::Ptr colored_cloud = reg.getColoredCloud ();
  pcl::visualization::CloudViewer viewer ("Cluster viewer");
  viewer.showCloud(colored_cloud);
  while (!viewer.wasStopped ())
  {
  }

  return (0);
}

The pcl::RegionGrowing class provides a method that returns the colored cloud where each cluster has its own color. So in this part of code the pcl::visualization::CloudViewer is instantiated for viewing the result of the segmentation - the same colored cloud. You can learn more about cloud visualization in the The CloudViewer tutorial.

Compiling and running the program

Add the following lines to your CMakeLists.txt file:

cmake_minimum_required(VERSION 3.5 FATAL_ERROR)

project(region_growing_segmentation)

find_package(PCL 1.5 REQUIRED)

include_directories(${PCL_INCLUDE_DIRS})
link_directories(${PCL_LIBRARY_DIRS})
add_definitions(${PCL_DEFINITIONS})

add_executable (region_growing_segmentation region_growing_segmentation.cpp)
target_link_libraries (region_growing_segmentation ${PCL_LIBRARIES})

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

$ ./region_growing_segmentation

After the segmentation the cloud viewer window will be opened and you will see something similar to those images:

_images/region_growing_segmentation_1.jpg

_images/region_growing_segmentation_2.jpg

In the last image you can see that the colored cloud has many red points. This means that these points belong to the clusters that were rejected, because they had too much/little points.