CGAL 6.1.1 - Shape Detection

Loading...
Searching...
No Matches
User Manual

Authors
  Sven Oesau, Yannick Verdie, Clément Jamin, Pierre Alliez, Florent Lafarge, Simon Giraudot, Thien Hoang, and Dmitry Anisimov

# Introduction

This CGAL component implements two algorithms for shape detection:

- the Efficient RANSAC (RANdom SAmple Consensus) method, contributed by Schnabel et al. [2];
- the Region Growing method, contributed by Lafarge and Mallet [1].

# Efficient RANSAC

From an unstructured point set with unoriented normals, this algorithm detects a set of shapes (see Figure 88.1). Five types of primitive shapes are provided by this package: plane, sphere, cylinder, cone, and torus. Other primitive shapes can be easily added by the user (see Section Custom Shapes).

Figure 88.1 Input and output of the Efficient RANSAC method. (a) Input point set. (b) Point set depicted with one color per detected shape.

This method takes as input a point set with unoriented normals and provides as output a set of detected shapes with associated input points. The output of the algorithm is a set of detected shapes with assigned points and all remaining points not covered by these shapes. Each input point can be assigned to at most one detected shape.

The shapes are detected via a RANSAC-type approach, that is 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 the highest number of inliers, referred to as the 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 minimum 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 is 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 threshold. 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 the 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).

## 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\). The parameter `epsilon` defines the absolute maximum tolerance Euclidean distance between a point and a shape. A high value of `epsilon` leads to the detection of fewer large shapes and hence a less detailed detection. A low value of `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` enables to detect more primitive shapes that approximate some of the free-form surfaces. The impact of this parameter is depicted by Figure 88.2. Its impact on performance is evaluated in Section Performance.

Figure 88.2 Impact of the epsilon parameter over the levels of detail 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`: The Efficient RANSAC uses this parameter to cluster the points into connected components covered by a detected shape. 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. The 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 the parameter `cluster_epsilon` is depicted in Figure 88.3.

Figure 88.3 The 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 of `cluster_epsilon` leads to detecting a single planar shape. (b) A moderate value of `cluster_epsilon` yields the detection of four squares. Notice that a few points within the squares are not detected as not connected. (c) A small value of `cluster_epsilon` leads to over-segmentation.

- `min_points`: The 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.
- `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 time due to a higher search endurance.

## Examples

The main class `Shape_detection::Efficient_RANSAC` takes a template parameter `Shape_detection::Efficient_RANSAC_traits` that defines the geometric types and input format. Property maps provide a means to interface with the user-specific data structures. The first parameter of the `Shape_detection::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::Efficient_RANSAC_traits` class. The concept behind property maps is detailed in Manual CGAL and Property Maps.

Typical usage consists of five steps:

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

### Basic Plane Detection

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

File Shape_detection/efficient_RANSAC_basic.cpp

#if defined (_MSC_VER) && !defined (_WIN64)
#pragma warning(disable:4244) // boost::number_distance::distance()
// converts 64 to 32 bits integers
#endif
#include <fstream>
#include <iostream>
#include <CGAL/property_map.h>
#include <CGAL/IO/read_points.h>
#include <CGAL/Point_with_normal_3.h>
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Shape_detection/Efficient_RANSAC.h>
// 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;
typedef CGAL::Shape_detection::Efficient_RANSAC_traits
<Kernel, Pwn_vector, Point_map, Normal_map> Traits;
typedef CGAL::Shape_detection::Efficient_RANSAC<Traits> Efficient_ransac;
typedef CGAL::Shape_detection::Plane<Traits> Plane;
int main (int argc, char** argv) {
std::cout << "Efficient RANSAC" << std::endl;
const std::string filename = (argc > 1) ? argv[1] : CGAL::data_file_path("points_3/cube.pwn");
// Points with normals.
Pwn_vector points;
// Load point set from a file.
if (!CGAL::IO::read_points(
filename,
std::back_inserter(points),
CGAL::parameters::point_map(Point_map()).
normal_map(Normal_map()))) {
std::cerr << "Error: cannot read file cube.pwn!" << std::endl;
return EXIT_FAILURE;
}
// Instantiate shape detection engine.
Efficient_ransac ransac;
// Provide input data.
ransac.set_input(points);
// Register planar shapes via template method.
ransac.add_shape_factory<Plane>();
// Detect registered shapes with default parameters.
ransac.detect();
// Print number of detected shapes.
std::cout << ransac.shapes().end() - ransac.shapes().begin()
<< " shapes detected." << std::endl;
return EXIT_SUCCESS;
}
Efficient_RANSAC.h
A convenience header that includes all classes related to the efficient RANSAC algorithm.
CGAL::Exact_predicates_inexact_constructions_kernel
CGAL::Shape_detection::Efficient_RANSAC
Shape detection algorithm based on the RANSAC method.
Definition: Efficient_RANSAC.h:60
CGAL::Shape_detection::Plane
Plane implements Shape_base.
Definition: Plane.h:33
CGAL::IO::read_points
bool read_points(const std::string &fname, PointOutputIterator output, const NamedParameters &np=parameters::default_values())
Kernel
CGAL::First_of_pair_property_map
CGAL::Second_of_pair_property_map
CGAL::Shape_detection::Efficient_RANSAC_traits
Default traits class for the CGAL::Shape_detection::Efficient_RANSAC.
Definition: Efficient_RANSAC_traits.h:41

### Plane Detection With Callback

The Efficient RANSAC class provides a callback mechanism that enables the user to track the progress of the algorithm. It can be used, for example, to terminate the algorithm based on a timeout. In the following example, the algorithm stops if it takes more than half a second and prints out the progress made.

File Shape_detection/efficient_RANSAC_with_callback.cpp

#if defined (_MSC_VER) && !defined (_WIN64)
#pragma warning(disable:4244) // boost::number_distance::distance()
// converts 64 to 32 bits integers
#endif
#include <fstream>
#include <iostream>
#include <CGAL/Timer.h>
#include <CGAL/property_map.h>
#include <CGAL/IO/read_points.h>
#include <CGAL/Point_with_normal_3.h>
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Shape_detection/Efficient_RANSAC.h>
// 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;
typedef CGAL::Shape_detection::Efficient_RANSAC_traits
<Kernel, Pwn_vector, Point_map, Normal_map> Traits;
typedef CGAL::Shape_detection::Efficient_RANSAC<Traits> Efficient_ransac;
typedef CGAL::Shape_detection::Plane<Traits> Plane;
struct Timeout_callback {
mutable int nb;
mutable CGAL::Timer timer;
const double limit;
Timeout_callback(double limit) :
nb(0), limit(limit) {
timer.start();
}
bool operator()(double advancement) const {
// Avoid calling time() at every single iteration, which could
// impact performances very badly.
++nb;
if (nb % 1000 != 0)
return true;
// If the limit is reached, interrupt the algorithm.
if (timer.time() > limit) {
std::cerr << "Algorithm takes too long, exiting ("
<< 100.0 * advancement << "% done)" << std::endl;
return false;
}
return true;
}
};
int main (int argc, char** argv) {
std::cout << "Efficient RANSAC" << std::endl;
const std::string filename = (argc > 1) ? argv[1] : CGAL::data_file_path("points_3/cube.pwn");
Pwn_vector points;
if (!CGAL::IO::read_points(
filename,
std::back_inserter(points),
CGAL::parameters::point_map(Point_map()).
normal_map(Normal_map()))) {
std::cerr << "Error: cannot read file cube.pwn!" << std::endl;
return EXIT_FAILURE;
}
Efficient_ransac ransac;
ransac.set_input(points);
ransac.add_shape_factory<Plane>();
// Create callback that interrupts the algorithm
// if it takes more than half a second.
Timeout_callback timeout_callback(0.5);
// Detect registered shapes with the default parameters.
ransac.detect(Efficient_ransac::Parameters(), timeout_callback);
return EXIT_SUCCESS;
}

### Setting Parameters And Using Different Shape Types

This example illustrates the user selection of parameters using the `Shape_detection::Efficient_RANSAC::Parameters` class. Shape detection is performed on five shape types (plane, cylinder, sphere, cone, and torus). The input point set is sampled on a surface mostly composed of piecewise planar and cylindrical parts, in addition to free-form parts.

Basic information of the detected shapes is written to the standard output: if the shape is either a plane or a cylinder, specific parameters are recovered, otherwise the general method `info()` is used to get the shape parameters in a string object. Note that specific parameters can be recovered for any of the provided shapes.

File Shape_detection/efficient_RANSAC_with_parameters.cpp

#include <fstream>
#include <iostream>
#include <CGAL/property_map.h>
#include <CGAL/IO/read_points.h>
#include <CGAL/Point_with_normal_3.h>
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Shape_detection/Efficient_RANSAC.h>
// 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;
typedef CGAL::Shape_detection::Efficient_RANSAC_traits
<Kernel, Pwn_vector, Point_map, Normal_map> Traits;
typedef CGAL::Shape_detection::Efficient_RANSAC<Traits> Efficient_ransac;
typedef CGAL::Shape_detection::Cone<Traits> Cone;
typedef CGAL::Shape_detection::Cylinder<Traits> Cylinder;
typedef CGAL::Shape_detection::Plane<Traits> Plane;
typedef CGAL::Shape_detection::Sphere<Traits> Sphere;
typedef CGAL::Shape_detection::Torus<Traits> Torus;
int main(int argc, char** argv) {
// Points with normals.
Pwn_vector points;
// Load point set from a file.
if (!CGAL::IO::read_points(
((argc > 1) ? argv[1] : CGAL::data_file_path("points_3/cube.pwn")),
std::back_inserter(points),
CGAL::parameters::point_map(Point_map()).
normal_map(Normal_map()))) {
std::cerr << "Error: cannot read file cube.pwn!" << std::endl;
return EXIT_FAILURE;
}
std::cout << points.size() << " points" << std::endl;
// Instantiate shape detection engine.
Efficient_ransac ransac;
// Provide 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>();
// Set parameters for shape detection.
Efficient_ransac::Parameters parameters;
// Set probability to miss the largest primitive at each iteration.
parameters.probability = 0.05;
// Detect shapes with at least 200 points.
parameters.min_points = 200;
// Set maximum Euclidean distance between a point and a shape.
parameters.epsilon = 0.002;
// Set maximum Euclidean distance between points to be clustered.
parameters.cluster_epsilon = 0.01;
// Set maximum normal deviation.
// 0.9 < dot(surface_normal, point_normal);
parameters.normal_threshold = 0.9;
// Detect shapes.
ransac.detect(parameters);
// Print 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()) {
// Get specific parameters depending on the detected shape.
if (Plane* plane = dynamic_cast<Plane*>(it->get())) {
Kernel::Vector_3 normal = plane->plane_normal();
std::cout << "Plane with normal " << normal << std::endl;
// Plane shape can also be converted to the Kernel::Plane_3.
std::cout << "Kernel::Plane_3: " <<
static_cast<Kernel::Plane_3>(*plane) << std::endl;
} else if (Cylinder* cyl = dynamic_cast<Cylinder*>(it->get())) {
Kernel::Line_3 axis = cyl->axis();
FT radius = cyl->radius();
std::cout << "Cylinder with axis "
<< axis << " and radius " << radius << std::endl;
} else {
// Print the parameters of the detected shape.
// This function is available for any type of shape.
std::cout << (*it)->info() << std::endl;
}
// Proceed with the next detected shape.
it++;
}
return EXIT_SUCCESS;
}
CGAL::Shape_detection::Cone
Cone implements Shape_base.
Definition: Cone.h:36
CGAL::Shape_detection::Cylinder
Cylinder implements Shape_base.
Definition: Cylinder.h:37
CGAL::Shape_detection::Sphere
Sphere implements Shape_base.
Definition: Sphere.h:32
CGAL::Shape_detection::Torus
Torus implements Shape_base.
Definition: Torus.h:35
Kernel::Line_3
Kernel::Plane_3
Kernel::Vector_3
CGAL::normal
CGAL::Vector_3< Kernel > normal(const CGAL::Point_3< Kernel > &p, const CGAL::Point_3< Kernel > &q, const CGAL::Point_3< Kernel > &r)

### 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 Shape_detection/efficient_RANSAC_with_point_access.cpp

#include <fstream>
#include <iostream>
#include <CGAL/Timer.h>
#include <CGAL/number_utils.h>
#include <CGAL/property_map.h>
#include <CGAL/IO/read_points.h>
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Shape_detection/Efficient_RANSAC.h>
// 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;
typedef CGAL::Shape_detection::Efficient_RANSAC_traits
<Kernel, Pwn_vector, Point_map, Normal_map> Traits;
typedef CGAL::Shape_detection::Efficient_RANSAC<Traits> Efficient_ransac;
typedef CGAL::Shape_detection::Plane<Traits> Plane;
int main(int argc, char** argv) {
// Points with normals.
Pwn_vector points;
// Load point set from a file.
if (!CGAL::IO::read_points(
((argc > 1) ? argv[1] : CGAL::data_file_path("points_3/cube.pwn")),
std::back_inserter(points),
CGAL::parameters::point_map(Point_map()).
normal_map(Normal_map()))) {
std::cerr << "Error: cannot read file cube.pwn!" << std::endl;
return EXIT_FAILURE;
}
// Instantiate shape detection engine.
Efficient_ransac ransac;
// Provide input data.
ransac.set_input(points);
// Register detection of planes.
ransac.add_shape_factory<Plane>();
// Measure time before setting up the shape detection.
CGAL::Timer time;
time.start();
// Build internal data structures.
ransac.preprocess();
// Measure time after preprocessing.
time.stop();
std::cout << "preprocessing took: " << time.time() * 1000 << "ms" << std::endl;
// Perform detection several times and choose result with the highest coverage.
Efficient_ransac::Shape_range shapes = ransac.shapes();
FT best_coverage = 0;
for (std::size_t i = 0; i < 3; ++i) {
// Reset timer.
time.reset();
time.start();
// Detect shapes.
ransac.detect();
// Measure 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());
// Print number of assigned shapes and unassigned 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 the 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()) {
std::shared_ptr<Efficient_ransac::Shape> shape = *it;
// Use Shape_base::info() to print the parameters of the detected shape.
std::cout << (*it)->info();
// Sums distances of points to the detected shapes.
FT sum_distances = 0;
// Iterate 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()) {
// Retrieve 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));
// Proceed with the next point.
index_it++;
}
// Compute and print the average distance.
FT average_distance = sum_distances / shape->indices_of_assigned_points().size();
std::cout << " average distance: " << average_distance << std::endl;
// Proceed with the next detected shape.
it++;
}
return EXIT_SUCCESS;
}
CGAL::sqrt
NT sqrt(const NT &x)

### Custom Shapes

Other shape types can be detected by implementing a shape class derived from the class `Shape_detection::Shape_base` and registering it to the shape detection factory of the Efficient RANSAC object. 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 to the shape, and compute the normal deviation between a query point with the 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 detecting 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::Shape_base`:

- `Shape_detection::Shape_base::minimum_sample_size()` const: Returns the minimum number of required samples.
- `Shape_detection::Shape_base::create_shape(const std::vector<size_t>& indices)`: The randomly generated samples are provided via a vector of indices. `Shape_detection::Shape_base::point(std::size_t index)` and `Shape_detection::Shape_base::normal(std::size_t index)` are used to retrieve the actual points and normals (see the example below). The provided number of samples might actually be larger than the above minimum number of required samples, depending on the other shape types. If the provided samples are not sufficient to define a unique shape, for example in a degenerated case, the shape is considered invalid.
- `Shape_detection::Shape_base::squared_distance(const Point& point)` 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::Shape_base::squared_distance(std::vector<FT>& distances, const std::vector<size_t>& indices)` and
- `Shape_detection::Shape_base::cos_to_normal(const std::vector<size_t>& indices, std::vector<FT>& angles)` const.

The last two 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 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::Shape_base::point(std::size_t index)` and `Shape_detection::Shape_base::normal(std::size_t index)` (see the 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::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 that do not belong 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::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 header shows an implementation of a planar shape primitive, which is used by the example Shape_detection/efficient_RANSAC_with_custom_shape.cpp.

File Shape_detection/include/efficient_RANSAC_with_custom_shape.h

#ifndef MY_PLANE_SHAPE_H
#define MY_PLANE_SHAPE_H
#include <CGAL/number_utils.h>
#include <CGAL/Shape_detection/Efficient_RANSAC.h>
// My_Plane is derived 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::Shape_base<Traits> {
public:
typedef typename Traits::FT FT; // number type
typedef typename Traits::Point_3 Point; // point type
typedef typename Traits::Vector_3 Vector; // vector type
My_Plane() :
CGAL::Shape_detection::Shape_base<Traits>()
{ }
// Compute 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;
}
// Return 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:
// Construct shape base on a 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;
}
// Compute 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;
}
}
// Compute the normal deviation between a 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);
}
// Return 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 // MY_PLANE_SHAPE_H
CGAL::Shape_detection::Shape_base
Base class for shape types that defines an interface to construct a shape from a set of points and to...
Definition: Shape_base.h:49
CGAL::abs
NT abs(const NT &x)
CGAL::squared_distance
Kernel::FT squared_distance(Type1< Kernel > obj1, Type2< Kernel > obj2)
CGAL

## Performance

The running time and detection performance of the Efficient RANSAC depend on the chosen parameters. A selective error tolerance parameter leads to higher running time and fewer 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 88.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 time. We plot the performance against the `probability` parameter (see Figure 88.5).

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

Figure 88.5 The graph depicts the runtime, coverage, and the number of detected primitives against the search endurance parameter, that is `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 time increases significantly as many more candidates are generated during each iteration of the algorithm.

# Region Growing

This shape detection component is based on the region growing algorithm applied to a set of user-specified items. Shapes are detected by growing regions from seed items, where each region is created as follows:

1. Pick the next available item;
2. Find its neighbors in the data set;
3. Include those neighbors, which satisfy the region requirements;
4. Repeat the procedure for all included neighbors;
5. If no further neighbor satisfies the requirements, start a new region.

Together with the generic algorithm's implementation `CGAL::Shape_detection::Region_growing`, five particular instances of this algorithm are provided:

- Line and circle detection in a 2D point set;
- Line detection in a 2D/3D segment set;
- Plane, sphere, and cylinder detection in a 3D point set;
- Plane detection on a polygon mesh.

Other instances can be easily added by the user, as explained below.

## Framework

The class `CGAL::Shape_detection::Region_growing` provides an implementation of the region growing algorithm. This algorithm can detect geometric primitives among a set of input geometric objects. Each such geometric object is identified using an item. The item is, in general but not necessarily, a lightweight class that uniquely identify a geometric object like an iterator over a container of geometric object, an integer pointing to the cell of a vector, a `face_descriptor` representing a face in a model of `FaceGraph`,... The geometric object is retrieved using a property map with the item as key type and the geometric object as value type. The class `CGAL::Shape_detection::Region_growing` is constructed by providing an input range with each element being converted to an item thanks to another property map passed to the constructor. For most cases the default is fine and conveniences alias and functions are provided for the `CGAL::Point_set_3` class (see Convenience Aliases and Functions for Point_set_3).

The algorithm provided by `CGAL::Shape_detection::Region_growing` works using two classes that are respectively models of the following concepts:

- `NeighborQuery`: responsible for providing the neighbors of an item;
- `RegionType`: responsible for the definition of a region from a list of items as well as for validating the addition of an item to an already defined region;

Additionally, a range of items can be provided as input seeds when constructing the `Region_growing` class. It defines the seeding order of items that is which items are used first to grow regions from. Such items are referred to as seed items. When not provided, the order used is that of the input range. Also note that the seed range may not contain all items of the `input_range`. In such case, items not provided and not reached by the region growing algorithm will have not region assigned.

Using this generic framework, users can grow any type of regions on a set of arbitrary items with their own propagation and seeding conditions (see an example).

### Examples

This toy example shows how to define a custom models of`NeighborQuery` and `RegionType` concept, which are used to parameterize the `CGAL::Shape_detection::Region_growing`.

We choose a simple custom object item. We define three such objects, where for each object, we manually provide the iterators to its neighbors. The operator `NeighborQuery::operator()()` does nothing but accessing those neighbors. The `Region_type` class defines the necessary types and functions:

- `Region_type::Primitive` This type represents the parameters of the primitive, e.g., in case of a sphere it could be a struct containing a `Point_3` for the center and a floating point number for the radius.
- `Region_type::Item` - a `const_iterator` of `std::vector<Object>`
- `Region_type::Region_index_map` - an unordered map from `Region_type::Item` to `std::size_t` encapsuled in a `boost::associative_property_map`.
- `Region_type::Primitive` - a `std::size_t`
- `Region_type::is_part_of_region()` - `true` if the first and second objects are neighbors,
- `Region_type::is_valid_region()` - always `true` after the first call to the function `update()`,
- `Region_type::update()` - updates the internal flag from the default `false` to `true`.

The result of using these classes with the region growing main class is that the first two objects form the first region, the third object forms the second region.

File Shape_detection/region_growing_with_custom_classes.cpp

// STL includes.
#include <unordered_map>
#include <list>
#include <vector>
#include <string>
#include <iostream>
#include <iterator>
#include <cassert>
// CGAL includes.
#include "include/utils.h"
#include <CGAL/Simple_cartesian.h>
#include <CGAL/Shape_detection/Region_growing/Region_growing.h>
// Custom Neighbor_query and Region_type classes for region growing.
namespace Custom {
// An object that stores indices of all its neighbors.
struct Object {
std::vector<std::vector<Object>::const_iterator> neighbors;
bool operator==(const Object& obj) const {
return neighbors == obj.neighbors;
}
};
// A range of objects.
using Objects = std::vector<Object>; // std::list<Object> works as well
// The Neighbor_query functor that accesses neighbors stored in
// the object struct above.
class Neighbor_query {
const Objects& m_objects;
using Item = typename Objects::const_iterator;
using Region = std::vector<Item>;
public:
Neighbor_query(const Objects& objects) :
m_objects(objects)
{ }
void operator()(
const Item &query,
std::vector<Item>& neighbors) const {
for (auto it = m_objects.begin(); it != m_objects.end(); it++) {
if (it == query) {
neighbors = query->neighbors;
return;
}
}
}
};
// The Region_type class, where the function is_part_of_region() verifies
// a very specific condition that the first and second objects in the
// range are in fact neighbors; is_valid_region() function always
// returns true after the first call to the update() function.
// These are the only functions that have to be defined.
class Region_type {
bool m_is_valid = false;
const std::vector<Object>& m_input;
public:
Region_type(const std::vector<Object> &input) : m_input(input) { }
using Primitive = std::size_t;
using Item = std::vector<Object>::const_iterator;
using Region = std::vector<Item>;
struct hash_item {
std::size_t operator()(Item i) const {
using boost::hash_value;
return boost::hash_value(i.operator->());
}
};
using Region_unordered_map = std::unordered_map<Item, std::size_t, hash_item >;
using Region_index_map = boost::associative_property_map<Region_unordered_map>;
Region_index_map region_index_map() {
Region_index_map index_map(m_region_map);
return index_map;
}
bool is_part_of_region(
const Item query,
const Region& region) const {
if (region.size() == 0) return false;
auto it = m_input.begin();
if (query == it || query == (it + 1))
return true;
return false;
}
inline bool is_valid_region(const Region&) const {
return m_is_valid;
}
Primitive primitive() {
return Primitive();
}
bool update(const Region&) {
m_is_valid = true;
return m_is_valid;
}
private:
Region_unordered_map m_region_map;
};
} // namespace Custom
// Typedefs.
using Object = Custom::Object;
using Objects = Custom::Objects;
using Neighbor_query = Custom::Neighbor_query;
using Region_type = Custom::Region_type;
using Region_growing = CGAL::Shape_detection::Region_growing<Neighbor_query, Region_type>;
int main() {
// Define a range of objects, where the first two objects form
// the first region, while the third object forms the second region.
// Note that Objects is a random access container here, however the
// same algorithm/example can work with other containers, e.g. std::list.
Objects objects(3);
auto it = objects.begin();
// Region 1.
objects[0].neighbors.push_back(it+1);
objects[1].neighbors.push_back(it);
// The single third object constitutes the second region.
std::cout << "* number of input objects: " << objects.size() << std::endl;
assert(objects.size() == 3);
// Create instances of the classes Neighbor_query and Region_type.
Neighbor_query neighbor_query = Neighbor_query(objects);
Region_type region_type = Region_type(objects);
// Create an instance of the region growing class.
Region_growing region_growing(
objects, neighbor_query, region_type);
// Run the algorithm.
std::vector<typename Region_growing::Primitive_and_region> regions;
region_growing.detect(std::back_inserter(regions));
std::cout << "* number of found regions: " << regions.size() << std::endl;
assert(regions.size() == 2);
return EXIT_SUCCESS;
}
CGAL::Shape_detection::Region_growing
Main class/entry point for running the region growing algorithm.
Definition: Region_growing.h:72
boost::hash_value
std::size_t hash_value(const typename CGAL::internal::Point_set_3_index< Point, Vector > &i)

## Point Set

If one wants to detect lines (see 2D Example)

Figure 88.6 A 2D point set depicted with one color per detected line.

or planes (see 3D Example)

Figure 88.7 A 3D point set depicted with one color per detected plane.

in a 2D or 3D point set respectively, this CGAL component provides the corresponding models of the concepts `NeighborQuery` and `RegionType`. In particular, it provides two different ways to define neighbors of a point:

- Fuzzy sphere neighbors search via `CGAL::Shape_detection::Point_set::Sphere_neighbor_query`. This class creates a circle (in 2D case) or a sphere (in 3D case) centered at the query point with a user-specified sphere radius. All points, which belong to the sphere, will be treated as neighbors of the query point;
- Nearest neighbors search via `CGAL::Shape_detection::Point_set::K_neighbor_query`. This class finds K (specified by the user) nearest neighbors of the query point either 2D or 3D.

The component also provides

- `CGAL::Shape_detection::Point_set::Least_squares_line_fit_region` - least squares line fit type of region for 2D points
- `CGAL::Shape_detection::Point_set::Least_squares_circle_fit_region` - least squares circle fit type of region for 2D points
- `CGAL::Shape_detection::Point_set::Least_squares_plane_fit_region` - least squares plane fit type of region for 3D points
- `CGAL::Shape_detection::Point_set::Least_squares_sphere_fit_region` - least squares sphere fit type of region for 3D points
- `CGAL::Shape_detection::Point_set::Least_squares_cylinder_fit_region` - least squares cylinder fit type of region for 3D points

The program associates all points from a region to the best-fit object (2D line, 2D circle, 3D plane, 3D sphere, etc.) and controls the quality of this fit.

The quality of region growing in a point set (2D or 3D) can be improved by slightly increasing the running time. To achieve this, one can sort the input points with respect to some quality criteria. These quality criteria can be included through the seed range (see Framework for more details). We provide a quality sorting both for 2D and 3D points:

- `CGAL::Shape_detection::Point_set::Least_squares_line_fit_sorting` - the 2D input points are sorted with respect to the quality of the least squares line fit applied to the neighbors of each point;
- `CGAL::Shape_detection::Point_set::Least_squares_circle_fit_sorting` - the 2D input points are sorted with respect to the quality of the least squares circle fit applied to the neighbors of each point;
- `CGAL::Shape_detection::Point_set::Least_squares_plane_fit_sorting` - the 3D input points are sorted with respect to the quality of the least squares plane fit applied to the neighbors of each point.
- `CGAL::Shape_detection::Point_set::Least_squares_sphere_fit_sorting` - the 3D input points are sorted with respect to the quality of the least squares sphere fit applied to the neighbors of each point.
- `CGAL::Shape_detection::Point_set::Least_squares_cylinder_fit_sorting` - the 3D input points are sorted with respect to the quality of the least squares cylinder fit applied to the neighbors of each point.

### Parameters

The classes in Section Region Growing On Point Set depend on a few parameters that should be defined by the user. They also have default values, but these values do not necessarily guarantee to produce pleasant results. All parameter values can be specified using Named Parameters.

The `NeighborQuery` related classes depend on the parameters:

- `sphere_radius` in the class `CGAL::Shape_detection::Point_set::Sphere_neighbor_query`, defines the radius of the fuzzy search sphere centered at the query point;
- `k_neighbors` in the class `CGAL::Shape_detection::Point_set::K_neighbor_query`, defines the number K of nearest neighbors of the query point.

The right choice of the parameters above plays an important role in producing a good result. For example, if we consider the fuzzy sphere neighborhood, when `sphere_radius` is too large, we have fewer regions, and the details are not clearly separated. Meanwhile, if `sphere_radius` is too small, we produce more regions, and the point set may be over-segmented. Consider a 2D map of an intersection of streets in a city as in Figure 88.8. Each region is painted with a unique color. As `sphere_radius` increases, the details become less clear. When `sphere_radius` = 0.3 (c), the best visual result is produced.

Figure 88.8 (a) Input 2D point set; (b) 17 regions are found when `sphere_radius` = 0.1; (c) 8 regions are found when `sphere_radius` = 0.3; (d) 4 regions are found when `sphere_radius` = 1.2.

The `RegionType` related classes depend on the following parameters:

- `maximum_distance` - the maximum distance from a point to a line/plane;
- `maximum_angle` - the maximum angle in degrees between the normal associated with a point and the normal of a line/plane;
- `minimum_region_size` - the minimum number of points a region must have.

The first two parameters are used by the functions `RegionType::is_part_of_region()` and `RegionType::update()`, while the third parameter is used by the function `RegionType::is_valid_region()` from the `RegionType` concept.

The right choice of `maximum_distance` and `maximum_angle` parameters is also very important. For example, Figure 88.9 shows that the roof top of the house can be distinguished as two planes (painted in blue and dark red) when `maximum_angle` is strict enough (c), or it can be recognized as only one plane (painted in pale yellow) in the other case (b).

Figure 88.9 (a) Input 3D point cloud; (b) Result when `maximum_angle` = 60 degrees; (c) Result when `maximum_angle` = 25 degrees.

### Examples

#### Detecting 2D Lines

Typical usage of the Region Growing component consists of five steps:

1. Define an input range with points;
2. Create instances of the classes `NeighborQuery` and `RegionType` with the proper parameters;
3. Create an instance of the class `CGAL::Shape_detection::Region_growing`;
4. Detect;
5. Postprocess.

Given a 2D point set, we detect 2D lines using the fuzzy sphere neighborhood. We then color all points from the found regions and save them in a file (see Figure 88.6). A point set with normals is stored in `std::vector`.

File Shape_detection/region_growing_lines_on_point_set_2.cpp

#include <CGAL/Simple_cartesian.h>
#include <CGAL/Shape_detection/Region_growing/Region_growing.h>
#include <CGAL/Shape_detection/Region_growing/Point_set.h>
#include "include/utils.h"
#include <boost/range/irange.hpp>
// Typedefs.
using Kernel = CGAL::Simple_cartesian<double>;
using FT = typename Kernel::FT;
using Point_2 = typename Kernel::Point_2;
using Vector_2 = typename Kernel::Vector_2;
using Point_with_normal = std::pair<Point_2, Vector_2>;
using Point_set_2 = std::vector<Point_with_normal>;
// we use Compose_property_map as the property maps are expected to operate on the item type
// std::size_t passed as parameter to Sphere_neighbor_query and Least_squares_line_fit_region
using Point_map = CGAL::Compose_property_map<CGAL::Random_access_property_map<Point_set_2>,
CGAL::First_of_pair_property_map<Point_with_normal> >;
using Normal_map = CGAL::Compose_property_map<CGAL::Random_access_property_map<Point_set_2>,
CGAL::Second_of_pair_property_map<Point_with_normal> >;
using Neighbor_query = CGAL::Shape_detection::Point_set::Sphere_neighbor_query<Kernel, std::size_t, Point_map>;
using Region_type = CGAL::Shape_detection::Point_set::Least_squares_line_fit_region<Kernel, std::size_t, Point_map, Normal_map>;
using Region_growing = CGAL::Shape_detection::Region_growing<Neighbor_query, Region_type>;
int main(int argc, char *argv[]) {
// Load xyz data either from a local folder or a user-provided file.
const bool is_default_input = argc > 1 ? false : true;
std::ifstream in(is_default_input ? CGAL::data_file_path("points_3/buildings_outline.xyz") : argv[1]);
CGAL::IO::set_ascii_mode(in);
if (!in) {
std::cerr << "ERROR: cannot read the input file!" << std::endl;
return EXIT_FAILURE;
}
FT a, b, c, d, e, f;
Point_set_2 point_set_2;
while (in >> a >> b >> c >> d >> e >> f) {
point_set_2.push_back(
std::make_pair(Point_2(a, b), Vector_2(d, e)));
}
in.close();
std::cout << "* number of input points: " << point_set_2.size() << std::endl;
assert(!is_default_input || point_set_2.size() == 3634);
// Default parameter values for the data file buildings_outline.xyz.
const FT sphere_radius = FT(5);
const FT max_distance = FT(45) / FT(10);
const FT max_angle = FT(45);
const std::size_t min_region_size = 5;
Point_map point_map(CGAL::make_random_access_property_map(point_set_2));
Normal_map normal_map(CGAL::make_random_access_property_map(point_set_2));
// Create instances of the classes Neighbor_query and Region_type.
Neighbor_query neighbor_query(
boost::irange<std::size_t>(0,point_set_2.size()),
CGAL::parameters::sphere_radius(sphere_radius)
.point_map(point_map));
Region_type region_type(
CGAL::parameters::
maximum_distance(max_distance).
maximum_angle(max_angle).
minimum_region_size(min_region_size).
normal_map(normal_map).
point_map(point_map));
// Create an instance of the region growing class.
Region_growing region_growing(
boost::irange<std::size_t>(0,point_set_2.size()), neighbor_query, region_type);
// Run the algorithm.
std::vector<typename Region_growing::Primitive_and_region> regions;
region_growing.detect(std::back_inserter(regions));
std::cout << "* number of found lines: " << regions.size() << std::endl;
assert(!is_default_input || regions.size() == 72);
// Save regions to a file.
const std::string fullpath = (argc > 2 ? argv[2] : "lines_point_set_2.ply");
utils::save_point_regions_2<Kernel, std::vector<typename Region_growing::Primitive_and_region>, Point_map>(
regions, fullpath, point_map);
return EXIT_SUCCESS;
}
Point_set.h
A convenience header that includes all classes related to the region growing algorithm on a point set...
CGAL::Shape_detection::Point_set::Least_squares_line_fit_region
Region type based on the quality of the least squares line fit applied to 2D points.
Definition: Least_squares_line_fit_region.h:60
CGAL::Shape_detection::Point_set::Sphere_neighbor_query
Fuzzy sphere neighbors search in a set of 2D or 3D points.
Definition: Sphere_neighbor_query.h:57
Kernel::Point_2
Kernel::Vector_2
CGAL::IO::set_ascii_mode
Mode set_ascii_mode(std::ios &s)
CGAL::parameters
CGAL::Compose_property_map
CGAL::Simple_cartesian

#### Detecting 2D Circles

The following example shows a similar example, this time detecting circles instead of lines. In that case, we also preprocess points so that they are sorted according to their likelihood of belonging to a circle:

File Shape_detection/region_growing_circles_on_point_set_2.cpp

#include <CGAL/Point_set_3.h>
#include <CGAL/Point_set_3/IO.h>
#include <CGAL/Simple_cartesian.h>
#include <CGAL/Shape_detection/Region_growing/Region_growing.h>
#include <CGAL/Shape_detection/Region_growing/Point_set.h>
#include <boost/iterator/function_output_iterator.hpp>
#include "include/utils.h"
// Typedefs.
using Kernel = CGAL::Simple_cartesian<double>;
using FT = Kernel::FT;
using Point_2 = Kernel::Point_2;
using Point_3 = Kernel::Point_3;
using Vector_2 = Kernel::Vector_2;
using Vector_3 = Kernel::Vector_3;
using Point_set_2 = CGAL::Point_set_3<Point_2, Vector_2>;
using Point_set_3 = CGAL::Point_set_3<Point_3, Vector_3>;
using Point_map = Point_set_2::Point_map;
using Normal_map = Point_set_2::Vector_map;
using Neighbor_query = CGAL::Shape_detection::Point_set::K_neighbor_query_for_point_set<Point_set_2>;
using Region_type = CGAL::Shape_detection::Point_set::Least_squares_circle_fit_region_for_point_set<Point_set_2>;
using Sorting = CGAL::Shape_detection::Point_set::Least_squares_circle_fit_sorting_for_point_set<Point_set_2, Neighbor_query>;
using Region_growing = CGAL::Shape_detection::Region_growing<Neighbor_query, Region_type>;
int main(int argc, char** argv) {
// Load ply data either from a local folder or a user-provided file.
const bool is_default_input = argc > 1 ? false : true;
std::ifstream in(is_default_input ? CGAL::data_file_path("points_3/circles.ply") : argv[1]);
CGAL::IO::set_ascii_mode(in);
if (!in) {
std::cerr << "ERROR: cannot read the input file!" << std::endl;
return EXIT_FAILURE;
}
Point_set_3 point_set_3;
in >> point_set_3;
in.close();
std::cout << "* number of input points: " << point_set_3.size() << std::endl;
assert(!is_default_input || point_set_3.size() == 1101);
assert(point_set_3.has_normal_map()); // input should have normals
// Create a 2D point set.
Point_set_2 point_set_2;
point_set_2.add_normal_map();
for (const auto& idx : point_set_3) {
const Point_3& point = point_set_3.point(idx);
const Vector_3& normal = point_set_3.normal(idx);
point_set_2.insert(
Point_2(point.x(), point.y()),
Vector_2(normal.x(), normal.y()));
}
// Default parameter values for the data file circles.ply.
const std::size_t k = 16;
const FT max_distance = FT(1) / FT(100);
const FT max_angle = FT(10);
const std::size_t min_region_size = 20;
// Create instances of the classes Neighbor_query and Region_type.
Neighbor_query neighbor_query = CGAL::Shape_detection::Point_set::make_k_neighbor_query(
point_set_2, CGAL::parameters::k_neighbors(k));
Region_type region_type = CGAL::Shape_detection::Point_set::make_least_squares_circle_fit_region(
point_set_2,
CGAL::parameters::
maximum_distance(max_distance).
maximum_angle(max_angle).
minimum_region_size(min_region_size));
// Sort indices.
Sorting sorting = CGAL::Shape_detection::Point_set::make_least_squares_circle_fit_sorting(
point_set_2, neighbor_query, CGAL::parameters::k_neighbors(k));
sorting.sort();
// Create an instance of the region growing class.
Region_growing region_growing(
point_set_2, sorting.ordered(), neighbor_query, region_type);
// Add maps to get a colored output.
Point_set_2::Property_map<unsigned char>
red = point_set_2.add_property_map<unsigned char>("red" , 0).first,
green = point_set_2.add_property_map<unsigned char>("green", 0).first,
blue = point_set_2.add_property_map<unsigned char>("blue" , 0).first;
// Run the algorithm.
CGAL::Random random;
std::size_t num_circles = 0;
region_growing.detect(
boost::make_function_output_iterator(
[&](const std::pair< Region_type::Primitive, typename Region_growing::Region >& region) {
// Assign a random color to each region.
const unsigned char r = static_cast<unsigned char>(random.get_int(64, 192));
const unsigned char g = static_cast<unsigned char>(random.get_int(64, 192));
const unsigned char b = static_cast<unsigned char>(random.get_int(64, 192));
for (auto item : region.second) {
red[item] = r;
green[item] = g;
blue[item] = b;
}
++num_circles;
}
)
);
std::cout << "* number of found circles: " << num_circles << std::endl;
assert(!is_default_input || num_circles == 10);
// Save regions to a file.
std::ofstream out("circles_point_set_2.ply");
CGAL::IO::set_ascii_mode(out);
out << point_set_2;
return EXIT_SUCCESS;
}
CGAL::Point_set_3
Kernel::Point_3

#### Detecting 3D Planes

If we are given a 3D point set, then the example below shows how to detect 3D planes using the K nearest neighbors search. We color all points from the found regions and save them in a file (see Figure 88.7). The example also shows how to retrieve all points, which are not assigned to any region, and how to use a custom output iterator. A point set with normals is stored in `CGAL::Point_set_3`.

File Shape_detection/region_growing_planes_on_point_set_3.cpp

#include <CGAL/Point_set_3.h>
#include <CGAL/Point_set_3/IO.h>
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Shape_detection/Region_growing/Region_growing.h>
#include <CGAL/Shape_detection/Region_growing/Point_set.h>
#include <boost/iterator/function_output_iterator.hpp>
#include "include/utils.h"
// Typedefs.
using Kernel = CGAL::Exact_predicates_inexact_constructions_kernel;
using FT = typename Kernel::FT;
using Point_3 = typename Kernel::Point_3;
using Point_set = CGAL::Point_set_3<Point_3>;
using Output_range = CGAL::Point_set_3<Point_3>;
using Point_map = typename Point_set::Point_map;
using Normal_map = typename Point_set::Vector_map;
using Region_type = CGAL::Shape_detection::Point_set::Least_squares_plane_fit_region_for_point_set<Point_set>;
using Neighbor_query = CGAL::Shape_detection::Point_set::K_neighbor_query_for_point_set<Point_set>;
using Sorting = CGAL::Shape_detection::Point_set::Least_squares_plane_fit_sorting_for_point_set<Point_set, Neighbor_query>;
using Region_growing = CGAL::Shape_detection::Region_growing<Neighbor_query, Region_type>;
using Point_inserter = utils::Insert_point_colored_by_region_index<Point_set::Index, Output_range, Point_map, Kernel::Plane_3>;
int main(int argc, char *argv[]) {
// Load xyz data either from a local folder or a user-provided file.
const bool is_default_input = argc > 1 ? false : true;
std::ifstream in(is_default_input ? CGAL::data_file_path("points_3/building.xyz") : argv[1]);
CGAL::IO::set_ascii_mode(in);
if (!in) {
std::cerr << "ERROR: cannot read the input file!" << std::endl;
return EXIT_FAILURE;
}
const bool with_normal_map = true;
Point_set point_set(with_normal_map);
in >> point_set;
in.close();
std::cout << "* number of input points: " << point_set.size() << std::endl;
assert(!is_default_input || point_set.size() == 8075);
// Default parameter values for the data file building.xyz.
const std::size_t k = 12;
const FT max_distance = FT(2);
const FT max_angle = FT(20);
const std::size_t min_region_size = 50;
// Create instances of the classes Neighbor_query and Region_type.
Neighbor_query neighbor_query = CGAL::Shape_detection::Point_set::make_k_neighbor_query(
point_set, CGAL::parameters::k_neighbors(k));
Sorting sorting = CGAL::Shape_detection::Point_set::make_least_squares_plane_fit_sorting(point_set, neighbor_query);
sorting.sort();
Region_type region_type = CGAL::Shape_detection::Point_set::make_least_squares_plane_fit_region(
point_set,
CGAL::parameters::
maximum_distance(max_distance).
maximum_angle(max_angle).
minimum_region_size(min_region_size));
// Create an instance of the region growing class.
Region_growing region_growing(
point_set, sorting.ordered(), neighbor_query, region_type);
// Run the algorithm.
Output_range output_range;
std::size_t number_of_regions = 0;
Point_inserter inserter(
point_set.point_map(),
output_range, number_of_regions);
region_growing.detect(
boost::make_function_output_iterator(inserter));
std::cout << "* number of found planes: " << number_of_regions << std::endl;
assert(!is_default_input || number_of_regions == 7);
// Save regions to a file.
const std::string fullpath = (argc > 2 ? argv[2] : "planes_point_set_3.ply");
std::ofstream out(fullpath);
out << output_range;
out.close();
// Get all unassigned points.
std::vector<Region_type::Item> unassigned_items;
region_growing.unassigned_items(point_set, std::back_inserter(unassigned_items));
std::cout << "* number of unassigned points: " << unassigned_items.size() << std::endl;
assert(!is_default_input || unassigned_items.size() == 538);
// Store all unassigned points.
std::vector<Point_3> unassigned_points;
unassigned_points.reserve(unassigned_items.size());
for (const Region_type::Item &item : unassigned_items) {
const Point_3& point = get(point_set.point_map(), item);
unassigned_points.push_back(point);
}
return EXIT_SUCCESS;
}

#### Detecting 3D Spheres

The following example shows a similar example, this time detecting spheres instead of planes.

File Shape_detection/region_growing_spheres_on_point_set_3.cpp

#include <CGAL/Point_set_3.h>
#include <CGAL/Point_set_3/IO.h>
#include <CGAL/Simple_cartesian.h>
#include <CGAL/Shape_detection/Region_growing/Region_growing.h>
#include <CGAL/Shape_detection/Region_growing/Point_set.h>
#include <boost/iterator/function_output_iterator.hpp>
#include "include/utils.h"
// Typedefs.
using Kernel = CGAL::Simple_cartesian<double>;
using FT = Kernel::FT;
using Point_3 = Kernel::Point_3;
using Vector_3 = Kernel::Vector_3;
using Point_set = CGAL::Point_set_3<Point_3>;
using Point_map = typename Point_set::Point_map;
using Normal_map = typename Point_set::Vector_map;
using Neighbor_query = CGAL::Shape_detection::Point_set::K_neighbor_query_for_point_set<Point_set>;
using Region_type = CGAL::Shape_detection::Point_set::Least_squares_sphere_fit_region_for_point_set<Point_set>;
using Region_growing = CGAL::Shape_detection::Region_growing<Neighbor_query, Region_type>;
int main(int argc, char** argv) {
// Load ply data either from a local folder or a user-provided file.
const bool is_default_input = argc > 1 ? false : true;
std::ifstream in(is_default_input ? CGAL::data_file_path("points_3/spheres.ply") : argv[1]);
CGAL::IO::set_ascii_mode(in);
if (!in) {
std::cerr << "ERROR: cannot read the input file!" << std::endl;
return EXIT_FAILURE;
}
Point_set point_set;
in >> point_set;
in.close();
std::cout << "* number of input points: " << point_set.size() << std::endl;
assert(!is_default_input || point_set.size() == 5969);
assert(point_set.has_normal_map()); // input should have normals
// Default parameter values for the data file spheres.ply.
const std::size_t k = 12;
const FT max_distance = FT(1) / FT(100);
const FT max_angle = FT(10);
const std::size_t min_region_size = 50;
// Create instances of the classes Neighbor_query and Region_type.
Neighbor_query neighbor_query = CGAL::Shape_detection::Point_set::make_k_neighbor_query(
point_set, CGAL::parameters::k_neighbors(k));
Region_type region_type = CGAL::Shape_detection::Point_set::make_least_squares_sphere_fit_region(
point_set,
CGAL::parameters::
maximum_distance(max_distance).
maximum_angle(max_angle).
minimum_region_size(min_region_size));
// Create an instance of the region growing class.
Region_growing region_growing(
point_set, neighbor_query, region_type);
// Add maps to get a colored output.
Point_set::Property_map<unsigned char>
red = point_set.add_property_map<unsigned char>("red" , 0).first,
green = point_set.add_property_map<unsigned char>("green", 0).first,
blue = point_set.add_property_map<unsigned char>("blue" , 0).first;
// Run the algorithm.
std::size_t num_spheres = 0;
region_growing.detect(
boost::make_function_output_iterator(
[&](const std::pair< Region_type::Primitive, typename Region_growing::Region>& region) {
// Assign a random color to each region.
const unsigned char r = static_cast<unsigned char>(CGAL::get_default_random().get_int(64, 192));
const unsigned char g = static_cast<unsigned char>(CGAL::get_default_random().get_int(64, 192));
const unsigned char b = static_cast<unsigned char>(CGAL::get_default_random().get_int(64, 192));
for (auto item : region.second) {
red[item] = r;
green[item] = g;
blue[item] = b;
}
++num_spheres;
}
)
);
std::cout << "* number of found spheres: " << num_spheres << std::endl;
assert(!is_default_input || num_spheres == 10);
// Save regions to a file.
std::ofstream out("spheres_point_set_3.ply");
CGAL::IO::set_ascii_mode(out);
out << point_set;
return EXIT_SUCCESS;
}

#### Detecting 3D Cylinders

The following example shows another similar example, this time detecting (infinite) cylinders.

File Shape_detection/region_growing_cylinders_on_point_set_3.cpp

#include <CGAL/Point_set_3.h>
#include <CGAL/Point_set_3/IO.h>
#include <CGAL/Simple_cartesian.h>
#include <CGAL/Shape_detection/Region_growing/Region_growing.h>
#include <CGAL/Shape_detection/Region_growing/Point_set.h>
#include <boost/iterator/function_output_iterator.hpp>
#include "include/utils.h"
// Typedefs.
using Kernel = CGAL::Simple_cartesian<double>;
using FT = Kernel::FT;
using Point_3 = Kernel::Point_3;
using Vector_3 = Kernel::Vector_3;
using Point_set = CGAL::Point_set_3<Point_3>;
using Point_map = typename Point_set::Point_map;
using Normal_map = typename Point_set::Vector_map;
using Neighbor_query = CGAL::Shape_detection::Point_set::K_neighbor_query_for_point_set<Point_set>;
using Region_type = CGAL::Shape_detection::Point_set::Least_squares_cylinder_fit_region_for_point_set<Point_set>;
using Region_growing = CGAL::Shape_detection::Region_growing<Neighbor_query, Region_type>;
int main(int argc, char** argv) {
// Load ply data either from a local folder or a user-provided file.
const bool is_default_input = argc > 1 ? false : true;
std::ifstream in(is_default_input ? CGAL::data_file_path("points_3/cylinders.ply") : argv[1]);
CGAL::IO::set_ascii_mode(in);
if (!in) {
std::cerr << "ERROR: cannot read the input file!" << std::endl;
return EXIT_FAILURE;
}
Point_set point_set;
in >> point_set;
in.close();
std::cout << "* number of input points: " << point_set.size() << std::endl;
assert(!is_default_input || point_set.size() == 1813);
assert(point_set.has_normal_map()); // input should have normals
// Default parameter values for the data file cuble.pwn.
const std::size_t k = 20;
const FT max_distance = FT(1) / FT(10);
const FT max_angle = FT(25);
const std::size_t min_region_size = 20;
// Create instances of the classes Neighbor_query and Region_type.
Neighbor_query neighbor_query = CGAL::Shape_detection::Point_set::make_k_neighbor_query(point_set, CGAL::parameters::k_neighbors(k));
Region_type region_type = CGAL::Shape_detection::Point_set::make_least_squares_cylinder_fit_region(
point_set,
CGAL::parameters::
maximum_distance(max_distance).
maximum_angle(max_angle).
minimum_region_size(min_region_size));
// Create an instance of the region growing class.
Region_growing region_growing(
point_set, neighbor_query, region_type);
// Add maps to get a colored output.
Point_set::Property_map<unsigned char>
red = point_set.add_property_map<unsigned char>("red" , 0).first,
green = point_set.add_property_map<unsigned char>("green", 0).first,
blue = point_set.add_property_map<unsigned char>("blue" , 0).first;
// Run the algorithm.
CGAL::Random random;
std::size_t num_cylinders = 0;
region_growing.detect(
boost::make_function_output_iterator(
[&](const std::pair< Region_type::Primitive, std::vector<typename Point_set::Index> >& region) {
// Assign a random color to each region.
const unsigned char r = static_cast<unsigned char>(random.get_int(64, 192));
const unsigned char g = static_cast<unsigned char>(random.get_int(64, 192));
const unsigned char b = static_cast<unsigned char>(random.get_int(64, 192));
for (auto id : region.second) {
put(red, id, r);
put(green, id, g);
put(blue, id, b);
}
++num_cylinders;
}
)
);
std::cout << "* number of found cylinders: " << num_cylinders << std::endl;
assert(!is_default_input || num_cylinders == 2);
// Save regions to a file.
std::ofstream out("cylinders_point_set_3.ply");
CGAL::IO::set_ascii_mode(out);
out << point_set;
return EXIT_SUCCESS;
}

### Performance

The main parameter that affects the region growing algorithm on a point set is the neighborhood size at each retrieval (see `sphere_radius` or `k_neighbors`). Larger neighbor lists are often followed by a smaller number of regions, larger coverage (the ratio between the number of points assigned to regions and the total number of input points), and longer running time. For example, for a test of about 70k 2D points with the fuzzy sphere neighborhood, the following table is produced:

`sphere_radius` | Time (in seconds) | Number of regions | Number of assigned points
--- | --- | --- | ---
1.0 | 0.138831 | 794 | 4483
3.0 | 0.069098 | 3063 | 63038
6.0 | 0.077703 | 2508 | 64906
9.0 | 0.093415 | 2302 | 65334

If the neighborhood size is set too low, some points might be isolated, the region size would not reach a critical mass and so will be discarded. This does not only cause the latency in the program, but also reduces the coverage value, as can be seen when the `sphere_radius = 1.0`. A typical runtime measure for a 3D point set with the K nearest neighborhood and well-defined parameters can be found in the following table:

`Number of points` | Time (in seconds)
--- | ---
300k | 0.761617
600k | 1.68735
900k | 2.80346
1200k | 4.06246

## Polygon Mesh

If one wants to detect planes on a polygon mesh, this CGAL component provides the corresponding models of the concepts `NeighborQuery` and `RegionType`. In particular, it provides

- `CGAL::Shape_detection::Polygon_mesh::One_ring_neighbor_query` class that retrieves all edge-adjacent faces of a face, and
- `CGAL::Shape_detection::Polygon_mesh::Least_squares_plane_fit_region` class that fits a 3D plane to the vertices of all mesh faces, which have been added to the region so far, and controls the quality of this fit.

This component accepts any model of the concept `FaceListGraph` as a polygon mesh. Figure 88.10 gives an example of the region growing algorithm for detecting 3D planes on `CGAL::Surface_mesh`.

Figure 88.10 A surface mesh depicted with one color per detected plane.

The quality of region growing on a polygon mesh can be improved by slightly increasing the running time. To achieve this, one can sort the input faces with respect to some quality criteria. These quality criteria can be included in region growing through the seed range (see Framework for more details). We provide such a quality sorting:

- `CGAL::Shape_detection::Polygon_mesh::Least_squares_plane_fit_sorting` - the input faces are sorted with respect to the quality of the least squares plane fit applied to the neighbors of each face.

### Parameters

The `NeighborQuery` related class does not require any parameters, because edge-adjacent faces are found using the internal face graph connectivity, while the `RegionType` related class depends on three parameters:

- `maximum_distance` - the maximum distance from the furthest face vertex to a plane;
- `maximum_angle` - the maximum angle in degrees between the face normal and the normal of a plane;
- `minimum_region_size` - the minimum number of mesh faces a region must have.

The first two parameters are used by the functions `RegionType::is_part_of_region()` and `RegionType::update()`, while the third parameter is used by the function `RegionType::is_valid_region()`. The right choice of these parameters is as important as the one explained in Section Parameters For Region Growing On Point Set.

### Examples

In the example below, we show how to use region growing to detect planes on a polygon mesh that can be either stored as `CGAL::Surface_mesh` or `CGAL::Polyhedron_3`.

We can improve the quality of region growing by providing a different seeding order (analogously to Point Sets) that is why in this example we also sort the input faces using the `CGAL::Shape_detection::Polygon_mesh::Least_squares_plane_fit_sorting` and only then detect regions.

File Shape_detection/region_growing_planes_on_polygon_mesh.cpp

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#ifdef USE_POLYHEDRON
#include <CGAL/Polyhedron_3.h>
#else
#include <CGAL/Surface_mesh.h>
#endif
#include <CGAL/Shape_detection/Region_growing/Region_growing.h>
#include <CGAL/Shape_detection/Region_growing/Polygon_mesh.h>
#include <CGAL/IO/polygon_mesh_io.h>
#include "include/utils.h"
// Typedefs.
using Kernel = CGAL::Exact_predicates_inexact_constructions_kernel;
using FT = typename Kernel::FT;
using Point_3 = typename Kernel::Point_3;
#ifdef USE_POLYHEDRON
using Polygon_mesh = CGAL::Polyhedron_3<Kernel>;
#else
using Polygon_mesh = CGAL::Surface_mesh<Point_3>;
#endif
using Neighbor_query = CGAL::Shape_detection::Polygon_mesh::One_ring_neighbor_query<Polygon_mesh>;
using Region_type = CGAL::Shape_detection::Polygon_mesh::Least_squares_plane_fit_region<Kernel, Polygon_mesh>;
using Sorting = CGAL::Shape_detection::Polygon_mesh::Least_squares_plane_fit_sorting<Kernel, Polygon_mesh, Neighbor_query>;
using Region_growing = CGAL::Shape_detection::Region_growing<Neighbor_query, Region_type>;
int main(int argc, char *argv[]) {
// Load data either from a local folder or a user-provided file.
const bool is_default_input = argc > 1 ? false : true;
const std::string filename = is_default_input ? CGAL::data_file_path("meshes/building.off") : argv[1];
std::ifstream in(filename);
CGAL::IO::set_ascii_mode(in);
Polygon_mesh polygon_mesh;
if (!CGAL::IO::read_polygon_mesh(filename, polygon_mesh)) {
std::cerr << "ERROR: cannot read the input file!" << std::endl;
return EXIT_FAILURE;
}
const auto& face_range = faces(polygon_mesh);
std::cout << "* number of input faces: " << face_range.size() << std::endl;
assert(!is_default_input || face_range.size() == 32245);
// Default parameter values for the data file building.off.
const FT max_distance = FT(1);
const FT max_angle = FT(45);
const std::size_t min_region_size = 5;
// Create instances of the classes Neighbor_query and Region_type.
Neighbor_query neighbor_query(polygon_mesh);
Region_type region_type(
polygon_mesh,
CGAL::parameters::
maximum_distance(max_distance).
maximum_angle(max_angle).
minimum_region_size(min_region_size));
// Sort face indices.
Sorting sorting(
polygon_mesh, neighbor_query);
sorting.sort();
// Create an instance of the region growing class.
Region_growing region_growing(
face_range, sorting.ordered(), neighbor_query, region_type);
// Run the algorithm.
std::vector<typename Region_growing::Primitive_and_region> regions;
region_growing.detect(std::back_inserter(regions));
std::cout << "* number of found planes: " << regions.size() << std::endl;
assert(!is_default_input || regions.size() == 365);
const Region_growing::Region_map& map = region_growing.region_map();
for (std::size_t i = 0; i < regions.size(); i++)
for (auto& item : regions[i].second) {
if (i != get(map, item)) {
std::cout << "Region map incorrect" << std::endl;
}
}
std::vector<typename Region_growing::Item> unassigned;
region_growing.unassigned_items(face_range, std::back_inserter(unassigned));
for (auto& item : unassigned) {
if (std::size_t(-1) != get(map, item)) {
std::cout << "Region map for unassigned incorrect" << std::endl;
}
}
// Save regions to a file.
const std::string fullpath = (argc > 2 ? argv[2] : "planes_polygon_mesh.ply");
utils::save_polygon_mesh_regions(polygon_mesh, regions, fullpath);
return EXIT_SUCCESS;
}
Polygon_mesh.h
A convenience header that includes all classes related to the region growing algorithm on a polygon m...
CGAL::Polyhedron_3
CGAL::Shape_detection::Polygon_mesh::Least_squares_plane_fit_region
Region type based on the quality of the least squares plane fit applied to faces of a polygon mesh.
Definition: Least_squares_plane_fit_region.h:57
CGAL::Shape_detection::Polygon_mesh::Least_squares_plane_fit_sorting
Sorting of polygon mesh faces with respect to the local plane fit quality.
Definition: Least_squares_plane_fit_sorting.h:53
CGAL::Shape_detection::Polygon_mesh::One_ring_neighbor_query
Edge-adjacent faces connectivity in a polygon mesh.
Definition: One_ring_neighbor_query.h:45
CGAL::Surface_mesh
CGAL::IO::read_polygon_mesh
bool read_polygon_mesh(const std::string &fname, Graph &g, const NamedParameters &np=parameters::default_values())

### Performance

Since accessing neighbors of a face in a polygon mesh is fast, performance of the region growing on a polygon mesh mostly depends on the `RegionType` model's implementation, which is usually fast, too.

## Segment Set

If one wants to detect lines in a segment set, this CGAL component provides the corresponding models of the concept `RegionType`. In particular, it provides

- `CGAL::Shape_detection::Segment_set::Least_squares_line_fit_region` class that fits a 2D or 3D line to the chunks of 2D or 3D segments, respectively, and controls the quality of this fit.

The quality of region growing in a segment set (2D or 3D) can be improved by slightly increasing the running time. To achieve this, one can sort the input segments with respect to some quality criteria. These quality criteria can be included through the seed range (see Framework for more details). We provide such a quality sorting:

- `CGAL::Shape_detection::Segment_set::Least_squares_line_fit_sorting` - the input segments are sorted with respect to the quality of the least squares line fit applied to the neighbors of each segment.

If the user's input is a polygon mesh defined as a `FaceListGraph`, it is also possible to filter out all graph edges, which split this polygon mesh into a set of regions. To achieve that, one first needs to apply the region growing on a polygon mesh and extract a set of such regions and then provide these regions to the class `CGAL::Shape_detection::Polygon_mesh::Polyline_graph` that returns a range of the corresponding edges, which can later be split into a set of linear regions, as shown in the example below.

File Shape_detection/region_growing_lines_on_segment_set.cpp

#include <CGAL/IO/PLY.h>
#include <CGAL/Surface_mesh.h>
#include <CGAL/Shape_detection/Region_growing/Region_growing.h>
#include <CGAL/Shape_detection/Region_growing/Segment_set.h>
#include <CGAL/Shape_detection/Region_growing/Polygon_mesh.h>
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/IO/polygon_mesh_io.h>
#include "include/utils.h"
// Typedefs.
using Kernel = CGAL::Exact_predicates_inexact_constructions_kernel;
using Point_3 = typename Kernel::Point_3;
using Surface_mesh = CGAL::Surface_mesh<Point_3>;
using Face_range = typename Surface_mesh::Face_range;
using Edge_range = typename Surface_mesh::Edge_range;
using One_ring_query = CGAL::Shape_detection::Polygon_mesh::One_ring_neighbor_query<Surface_mesh>;
using Plane_region = CGAL::Shape_detection::Polygon_mesh::Least_squares_plane_fit_region<Kernel, Surface_mesh>;
using RG_planes = CGAL::Shape_detection::Region_growing<One_ring_query, Plane_region>;
using Polyline_graph = CGAL::Shape_detection::Polygon_mesh::Polyline_graph<Surface_mesh>;
using Segment_range = typename Polyline_graph::Segment_range;
using Segment_map = typename Polyline_graph::Segment_map;
using Line_region = CGAL::Shape_detection::Segment_set::Least_squares_line_fit_region<Kernel, Surface_mesh::Edge_index, Segment_map>;
using Line_sorting = CGAL::Shape_detection::Segment_set::Least_squares_line_fit_sorting<Kernel, Surface_mesh::Edge_index, Polyline_graph, Segment_map>;
using RG_lines = CGAL::Shape_detection::Region_growing<Polyline_graph, Line_region>;
int main(int argc, char *argv[]) {
// Load data either from a local folder or a user-provided file.
const bool is_default_input = argc > 1 ? false : true;
const std::string filename = is_default_input ? CGAL::data_file_path("meshes/am.off") : argv[1];
Surface_mesh surface_mesh;
if (!CGAL::IO::read_polygon_mesh(filename, surface_mesh)) {
std::cerr << "ERROR: cannot read the input file!" << std::endl;
return EXIT_FAILURE;
}
const Face_range face_range = faces(surface_mesh);
const Edge_range edge_range = edges(surface_mesh);
std::cout << "* number of input faces: " << face_range.size() << std::endl;
std::cout << "* number of input edges: " << edge_range.size() << std::endl;
assert(!is_default_input || face_range.size() == 7320);
assert(!is_default_input || edge_range.size() == 10980);
// Find planar regions.
One_ring_query one_ring_query(surface_mesh);
Plane_region plane_region(surface_mesh);
RG_planes rg_planes(face_range, one_ring_query, plane_region);
std::vector<typename RG_planes::Primitive_and_region> regions;
rg_planes.detect(std::back_inserter(regions));
std::cout << "* number of found planar regions: " << regions.size() << std::endl;
assert(!is_default_input || regions.size() == 9);
std::string fullpath = (argc > 2 ? argv[2] : "regions_sm.ply");
utils::save_polygon_mesh_regions(surface_mesh, regions, fullpath);
// Find linear regions.
Polyline_graph pgraph(surface_mesh, rg_planes.region_map());
const auto& segment_range = pgraph.segment_range();
std::cout << "* number of extracted segments: " << segment_range.size() << std::endl;
Line_region line_region(CGAL::parameters::segment_map(pgraph.segment_map()));
Line_sorting line_sorting(
segment_range, pgraph, CGAL::parameters::segment_map(pgraph.segment_map()));
line_sorting.sort();
RG_lines rg_lines(
segment_range, line_sorting.ordered(), pgraph, line_region);
std::vector<typename RG_lines::Primitive_and_region> subregions;
rg_lines.detect(std::back_inserter(subregions));
std::cout << "* number of found linear regions: " << subregions.size() << std::endl;
assert(!is_default_input || subregions.size() == 21);
fullpath = (argc > 2 ? argv[2] : "subregions_sm.ply");
utils::save_segment_regions_3<Kernel, std::vector<typename RG_lines::Primitive_and_region>, Segment_map>(
subregions, fullpath, pgraph.segment_map());
return EXIT_SUCCESS;
}
Segment_set.h
A convenience header that includes all classes related to the region growing algorithm on a segment s...
CGAL::Shape_detection::Polygon_mesh::Polyline_graph
Polygon mesh edges connected into a graph.
Definition: Polyline_graph.h:47
CGAL::Shape_detection::Segment_set::Least_squares_line_fit_region
Region type based on the quality of the least squares line fit applied to a segment set.
Definition: Least_squares_line_fit_region.h:56
CGAL::Shape_detection::Segment_set::Least_squares_line_fit_sorting
Sorting of segments with respect to the local line fit quality.
Definition: Least_squares_line_fit_sorting.h:53
CGAL::Surface_mesh::Face_range
unspecified_type Face_range
CGAL::Surface_mesh::Edge_range
unspecified_type Edge_range

### Parameters

The line fitting class above depends on three parameters:

- `maximum_distance` - the maximum distance from the furthest segment vertex to a line;
- `maximum_angle` - the maximum angle in degrees between the segment direction and the direction of a line;
- `minimum_region_size` - the minimum number of segments a region must have.

The right choice of these parameters is important for producing the good results.

# Comparison

The Efficient RANSAC algorithm is very quick, however, since it is not deterministic, some small shapes might be missed in the detection process.

Instead, the region growing algorithm usually takes longer to complete, but it may provide better quality output in the presence of large scenes with numerous small details. Since it iterates throughout all items of the scene, there are fewer chances to miss a shape. In addition, it is deterministic (for a given input and a given set of parameters, it always returns the same output, whereas the Efficient RANSAC algorithm is randomized and so the output may vary at each run), see Figure 88.11.

Figure 88.11 Comparison of the Efficient RANSAC and region growing algorithms. Top: the input point set. Bottom left: the output of the Efficient RANSAC, \(78\%\) of the shapes are correctly detected in 8 seconds. Bottom right: the output of the region growing, \(100\%\) of the shapes detected in 15 seconds. Unassigned points are in black in both output images.

# History

The Efficient RANSAC component was developed by Sven Oesau based on the prototype version created by Yannick Verdie, with the help of Clément Jamin and under the supervision of Pierre Alliez.

The region growing algorithm on a 3D point set was first implemented by Simon Giraudot based on the prototype version developed by Florent Lafarge and then generalized to arbitrary items including versions for a 2D point set, a 3D point set, and a polygon mesh by Thien Hoang during the Google Summer of Code 2018 under the supervision of Dmitry Anisimov. The version for segments was added later by Dmitry Anisimov.

# Acknowledgments

The authors wish to thank our reviewers Andreas Fabri, Sébastien Loriot, and Simon Giraudot for helpful comments and discussions.