User Manual

Authors

Sven Oesau, Yannick Verdie, Clément Jamin, Pierre Alliez

Introduction

This CGAL component implements the efficient RANSAC method for shape detection, contributed by Schnabel et al. [1]. From an unstructured point set with unoriented normals, the algorithm detects a set of shapes (see Figure 69.1). Five types of primitive shapes are provided by this package: plane, sphere, cylinder, cone and torus. Detecting other types of shapes requires implementing a class derived from a base shape.

Figure 69.1 (a) Input point set. (b) Point set depicted with one color per detected shape.

Method

The method takes as input a point set with unoriented normals and provides as output a set of detected shapes with associated input points. The shapes are detected via a RANSAC-type approach, i.e., a random sample consensus. The basic RANSAC approach repeats the following steps: (1) Randomly select samples from the input points. (2) Fit a shape to the selected samples. (3) Count the number of inliers to the shape, inliers being within a user-specified error tolerance to the shape. Steps 1-3 are repeated for a prescribed number of iterations and the shape with highest number of inliers, referred to as largest shape, is kept.

In our context, the error between a point and a shape is defined by its distance and normal deviation to the shape. A random subset corresponds to the minimal number of points (with normals) required to uniquely define a primitive.

For very large point sets the basic RANSAC method is not practical when testing all possible shape candidates against the input data in order to find the largest shape. The main idea behind the efficient RANSAC method consists in testing shape candidates against subsets of the input data. Shape candidates are constructed until the probability to miss the largest candidate is lower than a user-specified parameter. The largest shape is repeatedly extracted until no more shapes, restricted to cover a minimum number of points, can be extracted. An additional gain in efficiency is achieved through exploiting the normal attributes during initial shape construction and enumeration of inliers.

The support of a shape refers to the footprint of the points covered by the primitive. To avoid generating shapes with fragmented support we enforce a connectivity constraint by considering only one connected component, referred to as cluster, selected as the one covering the largest number of inliers. See Section Parameters for more details.

The output of the shape detection algorithm is a set of detected shapes with assigned points, and all remaining points not covered by shapes. Each input point can be assigned to at most one detected shape.

Parameters

The algorithm has five parameters:

• epsilon and normal_threshold: The error between a point-with-normal \(p\) and a shape \(S\) is defined by its Euclidean distance and normal deviation to \(S\). The normal deviation is computed between the normal at \(p\) and the normal of \(S\) at the closest projection of \(p\) onto \(S\). Parameter epsilon defines the absolute maximum tolerance Euclidean distance between a point and a shape. A high value for epsilon leads to the detection of fewer large shapes and hence a less detailed detection. A low value for epsilon yields a more detailed detection, but may lead to either lower coverage or over-segmentation. Over-segmentation translates into detection of fragmented shapes when epsilon is within or below the noise level. When the input point set is made of free-form parts, a higher tolerance epsilon allows for detecting more primitive shapes that approximate some of the free-form surfaces. The impact of this parameter is depicted by Figure Figure 69.2. Its impact on performance is evaluated in Section Performance.
  
  Figure 69.2 Impact of epsilon over level of details of the detection. (a) Input point set. (b) Detection of planar shapes with epsilon set to 2.0 (one color per detected shape). Most details such as chimneys on the roof are not distinguished. (c) Detection with epsilon set to 0.5. The facades are correctly detected and some details of the roof are detected. (d) Setting epsilon to 0.25 yields a more detailed but slightly over-segmented detection.
  
  
• cluster_epsilon: Clustering of the points into connected components covered by a detected shape is controlled via parameter cluster_epsilon. For developable shapes that admit a trivial planar parameterization (plane, cylinder, cone) the points covered by a shape are mapped to a 2D parameter space chosen to minimize distortion and best preserve arc-length distances. This 2D parameter space is discretized using a regular grid, and a connected component search is performed to identify the largest cluster. Parameter cluster_epsilon defines the spacing between two cells of the regular grid, so that two points separated by a distance of at most \(2\sqrt{2}\), cluster_epsilon are considered adjacent. For non-developable shapes the connected components are identified by computing a neighboring graph in 3D and walking in the graph. The impact of parameter cluster_epsilon is depicted by Figure Figure 69.3.
  
  Figure 69.3 Parameter cluster_epsilon controls the connectivity of the points covered by a detected shape. The input point set is sampled on four coplanar squares. (a) A large value for cluster_epsilon leads to detecting a single planar shape. (b) A moderate value for cluster_epsilon yields the detection of four squares. Notice that few points within the squares are not detected as not connected. (c) A small value for cluster_epsilon leads to over-segmentation.
  
  
• probability: This parameter defines the probability to miss the largest candidate shape. A lower probability provides a higher reliability and determinism at the cost of longer running times due to higher search endurance.
• min_points: This minimum number of points controls the termination of the algorithm. The shape search is iterated until no further shapes can be found with a higher support. Note that this parameter is not strict: depending on the chosen probability, shapes may be extracted with a number of points lower than the specified parameter.

Examples

The main class Shape_detection_3::Efficient_RANSAC takes a template parameter Shape_detection_3::Efficient_RANSAC_traits that defines the geometric types and input format. Property maps provide a means to interface with user-specific data structures. The first parameter of the Shape_detection_3::Efficient_RANSAC_traits class is the common Kernel. In order to match the constraints of property maps, an iterator type and two maps that map an iterator to a point and a normal are specified in the Shape_detection_3::Efficient_RANSAC_traits class. The concept behind property maps is detailed in a dedicated chapter.

Typical usage consists in five steps:

1. Provide input data via a range iterator
2. Register shape factories
3. Choose parameters
4. Detect
5. Retrieve detected shapes

Basic Planar Shape Detection

The following minimal example reads a point set from a file and detects only planar shapes. The default parameters are used for detection.

File Point_set_shape_detection_3/efficient_RANSAC_basic.cpp

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/IO/read_xyz_points.h>
#include <CGAL/Point_with_normal_3.h>
#include <CGAL/property_map.h>

#include <CGAL/Shape_detection_3.h>

#include <iostream>
#include <fstream>

// Type declarations
typedef CGAL::Exact_predicates_inexact_constructions_kernel  Kernel;
typedef std::pair<Kernel::Point_3, Kernel::Vector_3>         Point_with_normal;
typedef std::vector<Point_with_normal>                       Pwn_vector;
typedef CGAL::First_of_pair_property_map<Point_with_normal>  Point_map;
typedef CGAL::Second_of_pair_property_map<Point_with_normal> Normal_map;

// In Efficient_RANSAC_traits the basic types, i.e., Point and Vector types
// as well as iterator type and property maps, are defined.
typedef CGAL::Shape_detection_3::Efficient_RANSAC_traits
  <Kernel, Pwn_vector, Point_map, Normal_map>                Traits;
typedef CGAL::Shape_detection_3::Efficient_RANSAC<Traits>    Efficient_ransac;
typedef CGAL::Shape_detection_3::Plane<Traits>               Plane;

int main() 
{
  // Points with normals.
  Pwn_vector points;

  // Loads point set from a file. 
  // read_xyz_points_and_normals takes an OutputIterator for storing the points
  // and a property map to store the normal vector with each point.
  std::ifstream stream("data/cube.pwn");

  if (!stream || 
    !CGAL::read_xyz_points_and_normals(stream,
      std::back_inserter(points),
      Point_map(),
      Normal_map()))
  {
      std::cerr << "Error: cannot read file cube.pwn" << std::endl;
      return EXIT_FAILURE;
  }

  // Instantiates shape detection engine.
  Efficient_ransac ransac;

  // Provides the input data.
  ransac.set_input(points);

  // Registers planar shapes via template method.
  ransac.add_shape_factory<Plane>();

  // Detects registered shapes with default parameters.
  ransac.detect();

  // Prints number of detected shapes.
  std::cout << ransac.shapes().end() - ransac.shapes().begin() << " shapes detected." << std::endl;

  return EXIT_SUCCESS;
}

Setting Parameters and Using Different Types of Shape

This example illustrates the user selection of parameters using the Shape_detection_3::Efficient_RANSAC::Parameters struct. Shape detection is performed on five types of shapes (plane, cylinder, sphere, cone, torus). Basic information of the detected shapes is written to the standard output. The input point set is sampled on a surface mostly composed of piece-wise planar and cylindrical parts, in addition to free-form parts.

File Point_set_shape_detection_3/efficient_RANSAC_parameters.cpp

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/IO/read_xyz_points.h>
#include <CGAL/Point_with_normal_3.h>
#include <CGAL/property_map.h>

#include <CGAL/Shape_detection_3.h>

#include <iostream>
#include <fstream>


// Type declarations
typedef CGAL::Exact_predicates_inexact_constructions_kernel  Kernel;
typedef Kernel::FT                                           FT;
typedef std::pair<Kernel::Point_3, Kernel::Vector_3>         Point_with_normal;
typedef std::vector<Point_with_normal>                       Pwn_vector;
typedef CGAL::First_of_pair_property_map<Point_with_normal>  Point_map;
typedef CGAL::Second_of_pair_property_map<Point_with_normal> Normal_map;

// In Efficient_RANSAC_traits the basic types, i.e., Point and Vector types
// as well as iterator type and property maps, are defined.
typedef CGAL::Shape_detection_3::Efficient_RANSAC_traits<Kernel,
  Pwn_vector, Point_map, Normal_map>            Traits;
typedef CGAL::Shape_detection_3::Efficient_RANSAC<Traits> Efficient_ransac;
typedef CGAL::Shape_detection_3::Cone<Traits>             Cone;
typedef CGAL::Shape_detection_3::Cylinder<Traits>         Cylinder;
typedef CGAL::Shape_detection_3::Plane<Traits>            Plane;
typedef CGAL::Shape_detection_3::Sphere<Traits>           Sphere;
typedef CGAL::Shape_detection_3::Torus<Traits>            Torus;


int main() 
{
  // Points with normals.
  Pwn_vector points;

  // Loads point set from a file. 
  // read_xyz_points_and_normals takes an OutputIterator for storing the points
  // and a property map to store the normal vector with each point.
  std::ifstream stream("data/cube.pwn");

  if (!stream ||
    !CGAL::read_xyz_points_and_normals(stream,
      std::back_inserter(points),
      Point_map(),
      Normal_map()))
  {
    std::cerr << "Error: cannot read file cube.pwn" << std::endl;
    return EXIT_FAILURE;
  }

  std::cout << points.size() << " points" << std::endl;

  // Instantiates shape detection engine.
  Efficient_ransac ransac;

  // Provides the input data.
  ransac.set_input(points);
    
  // Register shapes for detection
  ransac.add_shape_factory<Plane>();
  ransac.add_shape_factory<Sphere>();
  ransac.add_shape_factory<Cylinder>();
  ransac.add_shape_factory<Cone>();
  ransac.add_shape_factory<Torus>();

  // Sets parameters for shape detection.
  Efficient_ransac::Parameters parameters;

  // Sets probability to miss the largest primitive at each iteration.
  parameters.probability = 0.05;
 
  // Detect shapes with at least 500 points.
  parameters.min_points = 200;

  // Sets maximum Euclidean distance between a point and a shape.
  parameters.epsilon = 0.002;
 
  // Sets maximum Euclidean distance between points to be clustered.
  parameters.cluster_epsilon = 0.01;
 
  // Sets maximum normal deviation.
  // 0.9 < dot(surface_normal, point_normal); 
  parameters.normal_threshold = 0.9;   
  
  // Detects shapes
  ransac.detect(parameters);

  // Prints number of detected shapes and unassigned points.
   std::cout << ransac.shapes().end() - ransac.shapes().begin() << " detected shapes, "
     << ransac.number_of_unassigned_points()
     << " unassigned points." << std::endl;
  
  // Efficient_ransac::shapes() provides
  // an iterator range to the detected shapes.
  Efficient_ransac::Shape_range shapes = ransac.shapes();
  Efficient_ransac::Shape_range::iterator it = shapes.begin();

  while (it != shapes.end()) {
    // Prints the parameters of the detected shape.
    std::cout << (*it)->info() << std::endl;

    // Proceeds with next detected shape.
    it++;
  }

  return EXIT_SUCCESS;
}

Retrieving Points Assigned to Shapes

This example illustrates how to access the points assigned to each shape and compute the mean error. A timer measures the running performance.

File Point_set_shape_detection_3/efficient_RANSAC_point_access.cpp

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/IO/read_xyz_points.h>
#include <CGAL/Point_with_normal_3.h>
#include <CGAL/property_map.h>
#include <CGAL/Timer.h>
#include <CGAL/number_utils.h>

#include <CGAL/Shape_detection_3.h>

#include <iostream>
#include <fstream>

// Type declarations
typedef CGAL::Exact_predicates_inexact_constructions_kernel  Kernel;
typedef Kernel::FT                                           FT;
typedef std::pair<Kernel::Point_3, Kernel::Vector_3>         Point_with_normal;
typedef std::vector<Point_with_normal>                       Pwn_vector;
typedef CGAL::First_of_pair_property_map<Point_with_normal>  Point_map;
typedef CGAL::Second_of_pair_property_map<Point_with_normal> Normal_map;

// In Efficient_RANSAC_traits the basic types, i.e., Point and Vector types
// as well as iterator type and property maps, are defined.
typedef CGAL::Shape_detection_3::Efficient_RANSAC_traits<Kernel,
  Pwn_vector, Point_map, Normal_map>            Traits;
typedef CGAL::Shape_detection_3::Efficient_RANSAC<Traits>   Efficient_ransac;
typedef CGAL::Shape_detection_3::Plane<Traits>              Plane;


int main()
{
  // Points with normals.
  Pwn_vector points;

  // Loads point set from a file. 
  // read_xyz_points_and_normals takes an OutputIterator for storing the points
  // and a property map to store the normal vector with each point.
  std::ifstream stream("data/cube.pwn");

  if (!stream ||
    !CGAL::read_xyz_points_and_normals(stream,
      std::back_inserter(points),
      Point_map(),
      Normal_map()))
  {
      std::cerr << "Error: cannot read file cube.pwn" << std::endl;
      return EXIT_FAILURE;
  }

  // Instantiates shape detection engine.
  Efficient_ransac ransac;

  // Provides the input data.
  ransac.set_input(points);

  // Registers detection of planes
  ransac.add_shape_factory<Plane>();

  // Measures time before setting up the shape detection.
  CGAL::Timer time;
  time.start();

  // Build internal data structures.
  ransac.preprocess();

  // Measures time after preprocessing.
  time.stop();

  std::cout << "preprocessing took: " << time.time() * 1000 << "ms" << std::endl;

  // Perform detection several times and choose result with highest coverage.
  Efficient_ransac::Shape_range shapes = ransac.shapes();
  FT best_coverage = 0;

  for (size_t i = 0;i<3;i++) {
    // Reset timer.
    time.reset();
    time.start();

    // Detects shapes.
    ransac.detect();

    // Measures time after detection.
    time.stop();

    // Compute coverage, i.e. ratio of the points assigned to a shape.
    FT coverage = FT(points.size() - ransac.number_of_unassigned_points())
                  / FT(points.size());

    // Prints number of assigned shapes and unsassigned points.
    std::cout << "time: " << time.time() * 1000 << "ms" << std::endl;
    std::cout << ransac.shapes().end() - ransac.shapes().begin() << " primitives, "
      << coverage << " coverage" << std::endl;
    
    // Choose result with highest coverage.
    if (coverage > best_coverage) {
      best_coverage = coverage;
      // Efficient_ransac::shapes() provides
      // an iterator range to the detected shapes. 
      shapes = ransac.shapes();
    }
  }

  Efficient_ransac::Shape_range::iterator it = shapes.begin();

  while (it != shapes.end()) {
    boost::shared_ptr<Efficient_ransac::Shape> shape = *it;
    // Using Shape_base::info() for printing 
    // the parameters of the detected shape.
    std::cout << (*it)->info();

    // Sums distances of points to detected shapes.
    FT sum_distances = 0;

    // Iterates through point indices assigned to each detected shape.
    std::vector<std::size_t>::const_iterator
      index_it = (*it)->indices_of_assigned_points().begin();

    while (index_it != (*it)->indices_of_assigned_points().end()) {
      
      // Retrieves point
      const Point_with_normal &p = *(points.begin() + (*index_it));

      // Adds Euclidean distance between point and shape.
      sum_distances += CGAL::sqrt((*it)->squared_distance(p.first));

      // Proceeds with next point.
      index_it++;
    }

    // Computes and prints average distance.
    FT average_distance = sum_distances / shape->indices_of_assigned_points().size();
    std::cout << " average distance: " << average_distance << std::endl;

    // Proceeds with next detected shape.
    it++;
  }

  return EXIT_SUCCESS;
}

Adding More Shapes

Other types of shapes can be detected by implementing a shape class derived from the class Shape_base and registering it to the shape detection factory. This class must provide the following functions: construct a shape from a small set of given points, compute the squared distance from a query point and the shape and compute the normal deviation between a query point with normal and the normal to the shape at the closest point from the query. The used shape parameters are added as members to the derived class.

Note that the RANSAC approach is efficient for shapes that are uniquely defined by a small number of points, denoted by the number of required samples. The algorithm aims at detected the largest shape via many random samples, and the combinatorial complexity of possible samples increases rapidly with the number of required samples.

More specifically, the functions to be implemented are defined in the base class Shape_detection_3::Shape_base:

• Shape_detection_3::Shape_base::minimum_sample_size() const: Returns the minimal number of required samples.
• Shape_detection_3::Shape_base::create_shape(const std::vector<size_t> &indices): The randomly generated samples are provided via a vector of indices. Shape_detection_3::Shape_base::point(size_t i) and Shape_detection_3::Shape_base::normal(size_t i) are used to retrieve the actual points and normals (see example below). The provided number of samples might actually be larger than the set minimal number of required samples, depending on the other types of shape types. If the provided samples are not sufficient for define a unique shape, e.g., in a degenerated case the shape is considered invalid.
• Shape_detection_3::Shape_base::squared_distance(const Point &p) const: This function computes the squared distance from a query point to the shape. It is used for traversing the hierarchical spatial data structure.
• Shape_detection_3::Shape_base::squared_distance(std::vector<FT> &dists, const std::vector<size_t> &indices)
• Shape_detection_3::Shape_base::cos_to_normal(const std::vector<size_t> &indice, sstd::vector<FT> &angles) const: These functions are used to determine the number of inlier points to the shape. They compute respectively the squared distance from a set of points to the shape, and the dot product between the point's normals and the normals at the shape for the closest points on the shape. The access to the actual point and normal data is carried out via Shape_detection_3::Shape_base::point(size_t i) and Shape_detection_3::Shape_base::normal(size_t i) (see example below). The resulting squared distance/dot product is stored in the vector provided as the first argument.

By default the connected component is detected via the neighbor graph as mentioned above. However, for shapes that admit a faster approach to detect a connected component, the user can provide his/her own implementation to extract the connected component via:

• Shape_detection_3::Shape_base::connected_component(std::vector<std::size_t>& indices, FT cluster_epsilon): The indices of all supporting points are stored in the vector indices. All points not belonging to the largest cluster of points are removed from the vector indices.

Another optional method can be implemented to provide a helper function providing the shape parameters written to a string:

• Shape_detection_3::Shape_base::info(): This function returns a string suitable for printing the shape parameters into a log/console. The default solution provides an empty string.

The property maps are used to map the indices to the corresponding points and normals. The following example shows an implementation of a planar shape primitive, which is used by the example Point_set_shape_detection_3/efficient_RANSAC_custom_shape.cpp.

File Point_set_shape_detection_3/efficient_RANSAC_custom_shape.h

#ifndef MY_PLANE_SHAPE_H
#define MY_PLANE_SHAPE_H

#include <CGAL/Shape_detection_3.h>
#include <CGAL/number_utils.h>

/*
My_Plane derives from Shape_base. The plane is represented by
its normal vector and distance to the origin.
*/
template <class Traits>
class My_Plane : public CGAL::Shape_detection_3::Shape_base<Traits> {
public:
  typedef typename Traits::FT FT;
  typedef typename Traits::Point_3 Point;
  typedef typename Traits::Vector_3 Vector;

public:
  My_Plane()
    : CGAL::Shape_detection_3::Shape_base<Traits>()
  {}

  //  Computes squared Euclidean distance from query point to the shape.
  virtual FT squared_distance(const Point &p) const {
    const FT sd = (this->constr_vec(m_point_on_primitive, p)) * m_normal;
    return sd * sd;
  }

  Vector plane_normal() const {
    return m_normal;
  }

  FT d() const {
    return m_d;
  }

  // Returns a string with shape parameters.
  virtual std::string info() const {
    std::stringstream sstr;
    sstr << "Type: plane (" << this->get_x(m_normal) << ", " 
      << this->get_y(m_normal) << ", " << this->get_z(m_normal) << ")x - " <<
      m_d << " = 0" << " #Pts: " << this->m_indices.size();

    return sstr.str();
  }

protected:
  // Constructs shape based on minimal set of samples from the input data.    
  virtual void create_shape(const std::vector<std::size_t> &indices) {
    const Point p1 = this->point(indices[0]);
    const Point p2 = this->point(indices[1]);
    const Point p3 = this->point(indices[2]);

    m_normal = this->cross_pdct(p1 - p2, p1 - p3);

    m_normal = m_normal * (1.0 / sqrt(this->sqlen(m_normal)));
    m_d = -(p1[0] * m_normal[0] + p1[1] * m_normal[1] + p1[2] * m_normal[2]);

    m_point_on_primitive = p1;

    this->m_is_valid = true;
  }

  // Computes squared Euclidean distance from a set of points.
  virtual void squared_distance(const std::vector<std::size_t> &indices,
                                std::vector<FT> &dists) const {
      for (std::size_t i = 0; i < indices.size(); i++) {
        const FT sd = (this->point(indices[i])
          - m_point_on_primitive) * m_normal;
        dists[i] = sd * sd;
      }
  }

  /*
  Computes the normal deviation between shape and
  a set of points with normals.
  */
  virtual void cos_to_normal(const std::vector<std::size_t> &indices,
                             std::vector<FT> &angles) const {
      for (std::size_t i = 0; i < indices.size(); i++)
        angles[i] = CGAL::abs(this->normal(indices[i]) * m_normal);
  }

  // Returns the number of required samples for construction.
  virtual std::size_t minimum_sample_size() const {
    return 3;
  }

private:
  Point m_point_on_primitive;
  Vector m_normal;
  FT m_d;
};
#endif

Performance

The running time and detection performance depend on the chosen parameters. A selective error tolerance parameter leads to higher running times and smaller shapes, as many shape candidates are generated to find the largest shape. We plot the detection performance against the epsilon error tolerance parameter for detecting planes in a complex scene with 5M points, see Figure 69.4. The probability parameter controls the endurance when searching for the largest candidate at each iteration. It barely impacts the number of detected shapes, has a moderate impact on the size of the detected shapes and increases the running times. We plot the performance against the probability parameter, see Figure 69.5.

Figure 69.4 The graph depicts the number of detected shapes (purple) and the coverage (green), i.e., the ratio assignedPoints / totalPoints, against the epsilon tolerance parameter. A higher value for epsilon, i.e., a more tolerant error, leads to fewer but larger shapes and shorter running times.

Figure 69.5 The graph depicts the time, coverage and number of detected primitives against the search endurance parameter, i.e., the probability to miss the largest shape at each iteration. The number of shapes is stable and the coverage increases when the probability is lowered. The running times increase significantly as many more candidates are generated during each iteration of the algorithm.