namespace CGAL {

/*!

\mainpage User Manual 
\anchor Chapter_Point_Set_Shape_Detection
\cgalAutoToc

\authors Sven Oesau, Yannick Verdie, Clément Jamin, Pierre Alliez, Florent Lafarge, Simon Giraudot

\section Point_set_shape_detection_3Introduction Introduction

This \cgal component implements two algorithms for shape detection:

- the efficient RANSAC method for shape detection, contributed by Schnabel et al. \cgalCite{schnabel2007efficient};
- the region growing method, contributed by Lafarge and Mallet \cgalCite{cgal:lm-clscm-12}.

From an unstructured point set with unoriented normals, these algorithms detect a set of shapes (see \cgalFigureRef{Point_set_shape_detection_3_overview}). Five types of primitive shapes are provided by this package: plane, sphere, cylinder, cone and torus. Other types of primitive shapes can also be added by the user.

Note that at the moment, region growing only detects planes.

\cgalFigureBegin{Point_set_shape_detection_3_overview,overview2.png}
(a) Input point set. (b) Point set depicted with one color per detected shape. 
\cgalFigureEnd

\section Point_set_shape_detection_3Method Methods

Both methods take as input a point set with unoriented normals and provide as output a set of detected shapes with associated input points. The output of the shape detection algorithms 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.

\subsection Point_set_shape_detection_3Method_RANSAC Efficient RANSAC

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 \a 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 \a 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 \ref Point_set_shape_detection_3Parameters for more details.

\subsection Point_set_shape_detection_3Method_region_growing Region Growing

Planes are detected by growing planar regions from estimated planar seeds. First, points are sorted by a local planarity measure (the fitting quality returned by `linear_least_squares_fitting_3()`), such that points whose local neighborhood is the most planar are treated first.

The following steps are repeated: (1) Pick the next available point in the sorted list. (2) Initialize a plane tangent to this point and perpendicular to the normal vector of this point. (3) Compute the local neighborhood of the point with a user-specified radius. (4) Find inliers among these points, i.e. points that are within a user-specified error tolerance to the plane (and with a normal deviation within a user-specified tolerance to the normal of the plane). Steps 3-4 are repeated until no other inlier can be added. The plane is periodically reestimated by principal component analysis using all the inliers. If, after the growing, the plane has less inliers than a user-specified minimum number of points, it is discarded and the points are made available for other shapes to take once again.

\subsection Point_set_shape_detection_3Method_comparison Comparisons

The efficient RANSAC is a very quick method. Because it is based on randomness, it is not deterministic and some small shapes might be missed in the detection process.

Region growing usually takes longer than the efficient RANSAC but may provide better quality output in the presence of large scenes with numerous small details: as it iterates throughout all points of the scene, there are less chances of missing a shape. In addition, it is deterministic (for a given input and a given set of parameters, it always returns the same output, whereas as RANSAC uses randomness to detect shapes, the output of a RANSAC algorithm may vary). See figure \cgalFigureRef{Point_set_shape_detection_3Method_comparison_figure}.

\cgalFigureBegin{Point_set_shape_detection_3Method_comparison_figure,comparison.png}
Comparison of efficient RANSAC and region growing. Top: input point set. Bottom left: output of efficient RANSAC, \f$78\%\f$ of the shapes correctly detected in 8 seconds. Bottom right: output of region growing, \f$100\%\f$ of the shapes detected in 15 seconds. Unassigned points are in black in both output images.
\cgalFigureEnd



\section Point_set_shape_detection_3Parameters Parameters

The algorithms have four parameters that are common. They have the same effects on the output for both RANSAC and region growing.

- `epsilon` and `normal_threshold`:
The error between a point-with-normal \f$p\f$ and a shape \f$S\f$ is defined by its Euclidean distance and normal deviation to \f$S\f$. The normal deviation is computed between the normal at \f$p\f$ and the normal of \f$S\f$ at the closest projection of \f$p\f$ onto \f$S\f$. 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 \cgalFigureRef{Point_set_shape_detection_3_parameter_epsilon_variation}. Its impact on performance is evaluated in Section \ref Point_set_shape_detection_3Performance.
\cgalFigureBegin{Point_set_shape_detection_3_parameter_epsilon_variation,epsilon_variation2.png}
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.
\cgalFigureEnd

- `cluster_epsilon`:
The region growing algorithm uses `cluster_epsilon` to compute the neighborhood around points when growing the regions while 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. Parameter `cluster_epsilon` defines the spacing between two cells of the regular grid, so that two points separated by a distance of at most \f$2\sqrt{2}\f$, `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 \cgalFigureRef{Point_set_shape_detection_3_parameter_connectivity}.

\cgalFigureBegin{Point_set_shape_detection_3_parameter_connectivity,varying_connectivity.png}
For efficient RANSAC, 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.
\cgalFigureEnd

- `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 for efficient RANSAC, this parameter is not strict: depending on the chosen probability, shapes may be extracted with a number of points lower than the specified parameter.

In addition, the efficient RANSAC has a fifth 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 times due to higher search endurance. 


\section Point_set_shape_detection_3Usage Examples
The main classes `Shape_detection_3::Efficient_RANSAC` and `Shape_detection_3::Region_growing` take a template parameter `Shape_detection_3::Shape_detection_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::Shape_detection_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::Shape_detection_traits` class. The concept behind property maps is detailed in \ref chapterProperty_map "CGAL and Propery Maps".

Typical usage consists of five steps:

-# Provide input data via a range iterator
-# Register shape factories
-# Choose parameters 
-# Detect
-# Retrieve detected shapes

\subsection Point_set_shape_detection_3Usage_minimal 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.

\cgalExample{Point_set_shape_detection_3/shape_detection_basic.cpp}

\subsection Point_set_shape_detection_3Usage_callback Planar Shape Detection with Callback

Both Region Growing and Efficient RANSAC provide a callback mechanism that enables the user to track the progress of the algorithms. 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:

\cgalExample{Point_set_shape_detection_3/shape_detection_with_callback.cpp}

\subsection Point_set_shape_detection_3Usage_parameters Setting Parameters and Using Different Types of Shape

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

Basic information of the detected shapes is written to the standard output: if the plane 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 shape.


\cgalExample{Point_set_shape_detection_3/efficient_RANSAC_parameters.cpp}

\subsection Point_set_shape_detection_3Usage_point_access 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.

\cgalExample{Point_set_shape_detection_3/efficient_RANSAC_point_access.cpp}

\section Point_set_shape_detection_3Arbitrary_shapes Adding More Shapes to RANSAC

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 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 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 to 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  \ref Point_set_shape_detection_3/efficient_RANSAC_custom_shape.cpp.

\cgalExample{Point_set_shape_detection_3/efficient_RANSAC_custom_shape.h}

\section Point_set_shape_detection_3Plane_regularization Plane Regularization

Shape detection applies to man-made scenes or objects such as urban scenes or scans of mechanical parts. Such scenes often contain a wide range of geometric regularities such as parallelism, orthogonality, or symmetry. This package offers a function to reinforce four types of regularities for planar shapes: `CGAL::regularize_planes()`:

- Planes that are near __parallel__ are made parallel: normal vectors of planes that form angles smaller than a user-defined threshold are made equal.

- Parallel planes that are near __coplanar__ are made coplanar.

- Planes that are near __orthogonal__ are made exactly orthogonal.

- Planes that are near __symmetrical__ with respect to a user-defined axis are made symmetrical.

The user can choose to only regularize one or several of these 4 properties (see \link CGAL::regularize_planes() reference manual\endlink)). The process is greedy and based on a hierarchical decomposition (coplanar clusters are subgroups of parallel clusters which are subgroups of axis-symmetric and orthogonal clusters) as described by Verdie et al. \cgalCite{cgal:vla-lod-15}

\cgalExample{Point_set_shape_detection_3/plane_regularization.cpp}

\section Point_set_shape_detection_3Performance Performances

\subsection Point_set_shape_detection_3PerformanceRANSAC Efficient RANSAC
 
The running time and detection performance of the efficient RANSAC 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 \cgalFigureRef{Point_set_shape_detection_3_performace_epsilon}. 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 \cgalFigureRef{Point_set_shape_detection_3_performace_probability}.

\cgalFigureBegin{Point_set_shape_detection_3_performace_epsilon,epsilon_graph.png}
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.
\cgalFigureEnd

\cgalFigureBegin{Point_set_shape_detection_3_performace_probability,prob_graph.png}
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.
\cgalFigureEnd

\subsection Point_set_shape_detection_3PerformanceRegionGrowing Region Growing

Region growing iterates through every single point in the input, which makes it usually slower (although more robust). The main parameter that has an effect on running times is the epsilon used for clustering: internally, range queries are perform using spheres with radius cluster_epsilon. Using larger values means using larger spheres which are more computationally demanding.

We plot again the detection performance against the epsilon parameter for detecting planes, this time in a scene of 500k points, see figure \cgalFigureRef{Point_set_shape_detection_3_performace_epsilon_rg}.

\cgalFigureBegin{Point_set_shape_detection_3_performace_epsilon_rg,epsilon_graph_2.png}
The graph depicts the number of detected shapes (purple) and the coverage (green), i.e., the ratio assignedPoints / totalPoints, against the epsilon parameter (for simplicity, we use the same value for epsilon and cluster_epsilon). A higher value for cluster_epsilon, i.e., larger neighborhood spheres, leads to longer computation and to larger shapes.
\cgalFigureEnd

\section Point_set_shape_detection_3History History

The efficient RANSAC implementation was developed by Sven Oesau based on a prototype version by Yannick Verdie, with the help of Clément Jamin and under the supervision of Pierre Alliez. Plane regularization and region growing were added by Simon Giraudot based on prototype versions developed by Florent Lafarge.


*/ 
} /* namespace CGAL */

