Merge branch 'Classification-GF-old' into Classification-GF

Classification/doc/Classification/Classification.txt

@@ -0,0 +1,292 @@
+namespace CGAL {
+/*!
+
+\mainpage User Manual
+\anchor Chapter_Classification
+\cgalAutoToc
+\author Simon Giraudot, Florent Lafarge
+
+This component implements the algorithm described in \cgalCite{cgal:lm-clscm-12} (section 2), generalized to handle multiple features and multiple labels. It classifies a data set into a user-defined set of labels, such as ground, vegetation and buildings. A flexible API is provided so that the user can classify any type of data, compute its own local features on the input data set and define its own labels based on these features.
+
+\section Classification_Organization Package Organization
+
+%Classification of data sets is achieved as follows (see Figure \cgalFigureRef{Classification_organization_fig}):
+
+- some analysis is performed on the input data set
+- features are computed based on this analysis
+- a set of labels (for example: ground, building, vegetation) is defined by the user
+- features are given weights and each pair of feature and label is assigned a relationship. This can either be done by hand or by automatic training if a set of inliers is given for each label (see \ref Classification_trainer)
+- classification is computed itemwise by minimizing an energy defined as the sum of the values taken by features on input items (which depend on the feature-label relationship)
+- additional regularization can be used by smoothing either locally or globally through an Alpha Expansion approach.
+
+\cgalFigureBegin{Classification_organization_fig,organization.png}
+Organization of the package.
+\cgalFigureEnd
+
+This package is designed to be easily extended by users: more specifically, features and labels can be defined by users to handle any data they need to classify (although \cgal provides a predefined framework for common urban scenes).
+
+Currently, \cgal provides data structures to handle classification of point sets.
+
+\section Classification_structures Data Structures
+
+\subsection Classification_analysis Analysis
+
+%Classification is based on the computation of local features. These features can take advantage of shared data structures that are precomputed and stored separately.
+
+\cgal provides the following structures:
+
+- [Point_set_neighborhood](@ref CGAL::Classification::Point_set_neighborhood) stores spatial searching structures and provides adapted queries for points
+- [Local_eigen_analysis](@ref CGAL::Classification::Local_eigen_analysis) precomputes covariance matrices on local neighborhoods of points and stores the associated eigenvectors and eigenvalues
+- [Planimetric_grid](@ref CGAL::Classification::Planimetric_grid) is a 2D grid used for digital terrain modeling.
+
+The following code snippet shows how to instantiate such data structures from an input PLY point set (the full example is given at the end of the manual). The class `CGAL::Classifier` that handles classification is also instantiated (see \ref Classification_classifier).
+
+\snippet Classification/example_classifier.cpp Analysis
+
+\subsection Classification_features Features
+
+Features are defined as scalar fields that associate each input item with a specific value. Users may want to define their own features, especially if the input data set comes with additional properties that were not anticipated by \cgal. A user-defined feature must inherit from `CGAL::Classification::Feature_base` and provide a method [value()](@ref CGAL::Classification::Feature_base::value) that associate a scalar value to each input item. Each feature has a weight that measure its strength with respect to the other features (see \ref Classification_labels).
+
+Features are accessed through `Handle` objects, `CGAL::Classification::Feature_handle`.
+
+\cgal provides some predefined features that are relevant for classification of urban point sets:
+
+- [Distance_to_plane](@ref CGAL::Classification::Feature::Distance_to_plane) measures how far away a point is from a locally estimated plane
+- [Elevation](@ref CGAL::Classification::Feature::Elevation) computes the local distance to an estimation of the ground
+- [Vertical_dispersion](@ref CGAL::Classification::Feature::Vertical_dispersion) computes how noisy the point set is on a local Z-cylinder
+- [Verticality](@ref CGAL::Classification::Feature::Verticality) compares the local normal vector to the vertical vector.
+
+For more details about how these different features can help to identify one label or the other, please refer to their associated reference manual pages. In addition, \cgal also provides features solely based on the eigenvalues \cgalCite{cgal:mbrsh-raofw-11} of the local covariance matrix:
+
+- [Anisotropy](@ref CGAL::Classification::Feature::Anisotropy)
+- [Eigentropy](@ref CGAL::Classification::Feature::Eigentropy)
+- [Linearity](@ref CGAL::Classification::Feature::Linearity)
+- [Omnivariance](@ref CGAL::Classification::Feature::Omnivariance)
+- [Planarity](@ref CGAL::Classification::Feature::Planarity)
+- [Sphericity](@ref CGAL::Classification::Feature::Sphericity)
+- [Sum_eigenvalues](@ref CGAL::Classification::Feature::Sum_eigenvalues)
+- [Surface_variation](@ref CGAL::Classification::Feature::Surface_variation)
+
+Finally, if the input data set has additional properties, these can also be used as features. For example, \cgal provides the following features:
+
+- [Echo_scatter](@ref CGAL::Classification::Feature::Echo_scatter) uses the number of returns (echo) provided by most LIDAR scanners if available
+- [Hsv](@ref CGAL::Classification::Feature::Hsv) uses input color information if available.
+
+In the following code snippet, a subset of these features are instantiated and their respective weights are set. Note that these weights can also be automatically set up by training (see \ref Classification_trainer).
+
+\snippet Classification/example_classifier.cpp Features
+
+Users can define their own feature classes. Such classes must fulfill the following requirements:
+
+- they must inherit `CGAL::Classification::Feature_base`
+- their constructor(s) must take as first argument the item range used by the classifier object
+- they must provide a method [value()](@ref CGAL::Classification::Feature_base::value) that associates each input point to a scalar value
+
+The following example shows how to define a feature that discriminates
+points that lie inside a 2D box from the others:
+
+\cgalExample{Classification/example_feature.cpp}
+
+
+\subsection Classification_labels Labels
+
+A label represents how an item should be classified, for example: vegetation, building, road, etc. It is defined by the values the features are expected to take for a specific label. For example, vegetation is expected to have a high distance to plane and have a color close to green (if colors are available); facades have a low distance to plane and a low verticality; etc.
+
+\cgal provides a class `CGAL::Classification::Label` to define such a set of feature effects, along with the associated `Handle` object: `CGAL::Classification::Label_handle`. Each label may define how a specific feature affects it:
+
+- [FAVORING](@ref CGAL::Classification::Feature::FAVORING): the label is favored by high values of the feature
+- [NEUTRAL](@ref CGAL::Classification::Feature::NEUTRAL): the label is not affected by the feature
+- [PENALIZING](@ref CGAL::Classification::Feature::PENALIZING): the label is favored by low values of the feature
+
+Let \f$x=(x_i)_{i=1..N_c}\f$ be a potential classification result with \f$N_c\f$ the number of input items and \f$x_i\f$ the class of the \f$i^{th}\f$ item (for example: vegetation, ground, etc.). Let \f$a_j(i)\f$ be the raw value of the \f$j^{th}\f$ feature at the \f$i^{th}\f$ item and \f$w_j\f$ be the weight of this feature. We define the normalized value \f$A_j(x_i) \in [0:1]\f$ of the \f$j^{th}\f$ feature at the \f$i^{th}\f$ item as follows:
+
+\f{eqnarray*}{
+    A_j(x_i) = & (1 - \min(\max(0,\frac{a_j(i)}{w_j}), 1)) & \mbox{if } a_j \mbox{ favors } x_i \\
+    & 0.5 & \mbox{if } a_j \mbox{ is neutral for } x_i \\
+    & \min(\max(0,\frac{a_j(i)}{w_j}), 1) & \mbox{if } a_j \mbox{ penalizes } x_i
+\f}
+
+
+The following code snippet shows how to add labels to the classification object and how to set the effects of the features on these labels. Note that these effects can also be automatically set up by training (see \ref Classification_trainer).
+
+\snippet Classification/example_classifier.cpp Labels
+
+\section Classification_classifier Classifier
+
+%Classification is performed by minizing an energy over the input data set that may include regularization. \cgal provides 3 different methods for classification, ranging from high speed / low quality to low speed / high quality:
+
+- `CGAL::Classifier::run()` 
+- `CGAL::Classifier::run_with_local_smoothing()`
+- `CGAL::Classifier::run_with_graphcut()`
+
+On a point set of 3 millions points, the first method takes about 4 seconds, the second about 40 seconds and the third about 2 minutes.
+
+\cgalFigureBegin{Classification_image,classif.png}
+Top-Left: input point set. Top-Right: raw output classification represented by a set of colors (ground is orange, facades are blue, roofs are pink and vegetation is green). Bottom-Left: output classification using local smoothing. Bottom-Right: output classification using graphcut.
+\cgalFigureEnd
+
+
+Mathematical details are provided hereafter.
+
+\subsection Classification_regularized_no No Regularization
+
+- `CGAL::Classifier::run()`: this is the fastest method
+that provides acceptable but usually noisy results (see Figure
+\cgalFigureRef{Classification_image}, top-right).
+
+Let \f$x=(x_i)_{i=1..N_c}\f$ be a potential classification result with \f$N_c\f$ the number of input items and \f$x_i\f$ the class of the \f$i^{th}\f$ item (for example: vegetation, ground, etc.). The classification is performed by minimizing the following energy:
+
+\f[
+  E(x) = \sum_{i = 1..N_c} E_{di}(x_i)
+\f]
+
+This energy is a sum of itemwise energies and involves no regularization. Let \f$A=(A_j)_{j=1..N_a}\f$ be the set of normalized features
+(see \ref Classification_labels). The itemwise energy measures the coherence of the class \f$x_i\f$ at the \f$i^{th}\f$ item and is defined as:
+
+\f[
+  E_{di}(x_i) = \sum_{j = 1..N_a} A_j(x_i)
+\f]
+
+
+
+\subsection Classification_regularized_local Local Smoothing
+
+- `CGAL::Classifier::run_with_local_smoothing()`: this
+method is a tradeoff between quality and efficiency (see Figure
+\cgalFigureRef{Classification_image}, bottom-left). The minimized
+energy is defined as follows:
+
+
+\f[
+  E(x) = \sum_{i = 1..N_c} E_{si}(x_i)
+\f]
+
+The itemwise energy \f$E_{di}(x_i)\f$ is replaced by
+\f$E_{si}(x_i)\f$ defined on a small local neighborhood \f$N(i)\f$ of
+the \f$i^{th}\f$ item:
+
+\f[
+  E_{si}(x_i) = \sum_{j = 1..N_a} \frac{\sum_{k \in N(i)} A_j(x_k)}{\left| N(i) \right|}
+\f]
+
+This allows to eliminate local noisy variations of assigned
+labels. Increasing the size of the neighborhood
+increases the noise reduction at the cost of higher computation times.
+
+
+\subsection Classification_regularized_graphcut Global Regularization (Graph Cut)
+
+- `CGAL::Classifier::run_with_graphcut()`: this method
+offers the best quality but requires longer computation time (see
+Figure \cgalFigureRef{Classification_image}, bottom-right). The
+total energy that is minimized is the sum of the partial data term
+\f$E_{di}(x_i)\f$ and of a pairwise interaction energy defined by the
+standard Potts model \cgalCite{cgal:l-mrfmi-09} :
+
+\f[
+  E(x) = \sum_{i = 1..N_c} E_{di}(x_i) + \gamma \sum_{i \sim j} \mathbf{1}_{x_i \neq x_j}
+\f]
+
+where \f$\gamma>0\f$ is the parameter of the Potts model that
+quantifies the strengh of the regularization, \f$i \sim j\f$
+represents the pairs of neighboring items and
+\f$\mathbf{1}_{\{.\}}\f$ the characteristic function.
+
+A Graph Cut based algorithm (Alpha Expansion) is used to quickly reach
+an approximate solution close to the global optimum of this energy.
+
+This method allows to consistently segment the input data set in
+piecewise constant parts and to correct large wrongly classified
+clusters. Increasing \f$\gamma\f$ produces more regular result with a
+constant computation time.
+
+
+\section Classification_point_set_classifier Point Set Classifier
+
+The classification algorithm is designed to be as flexible as possible: users may classify any type of item they need and define their own features and labels. Nevertheless, \cgal provides a predefined specialization of `CGAL::Classifier` to handle common urban point sets: the class `CGAL::Point_set_classifier`. In addition to all the methods inherited from `CGAL::Classifier`, it provides additional features:
+
+- it takes care of generating all needed analysis structures
+- it generates all possible features (among all the \cgal predefined ones) based on which property maps are available (it uses colors if available, etc.)
+- multiple scales can be used to increase the quality of the results \cgalCite{cgal:hws-fsso3-16}
+- input/ouput methods are provided to save and recover a specific configuration (with all features, labels and relationships between them)
+- classification can be saved as a PLY format with colors and labels.
+
+Note that using this class in order to classify point sets is not mandatory, as `CGAL::Classifier` can handle points as well. It is mainly provided to simplify the whole process of defining data structures and features in the specific case of urban point sets. Users can still add their own data structures and features within `CGAL::Point_set_classifier` similarly to what can be done with `CGAL::Classifier`.
+
+The example \ref Classification_example_training shows how to generate features and save the configuration and classification.
+
+\section Classification_trainer Trainer
+
+%Classification is based on relationships between features and labels. Each feature has a specific weight and each pair of feature-label has a specific effect. This means that the number of parameters to set up can quickly explode: if 6 features are used to classify between 4 labels, 30 parameters have to be set up (6 weights + 6x4 feature-label relationships).
+
+Though it is possible to set them up one by one, \cgal also provides a trainer class `CGAL::Classification::Trainer` that requires a small set of ground truth items provided by users. More specifically, users must provide, for each label they want to classify, a set of known inliers among the input data set (for example, selecting one roof, one tree and one section of the ground). The training algorithm works as follows:
+
+- for each feature, a range of weights is tested: the effect each feature have on each label is estimated. For a given weight, if a feature has the same effect on each label, it is non-relevant for classification. The range of weights such that the feature is relevant is estimated
+
+- for each feature, uniformly picked weight values are tested and their effects estimated
+
+- each inlier provided by the user is classified using this set of weights and effects
+
+- for each label, the ratio of correctly classified inliers is computed. The minimum of these ratios is used as a score for this set of weights and effects: a ratio of 0.8 means that for each label, at least 80\% of the provided inliers were correctly classified
+
+- the same mechanism is repeated until all features' ranges have been tested. Weights are only changed one by one, the other ones kept to the previous value that gave the best score.
+
+This usually converges to a satisfying solution (see Figure \cgalFigureRef{Classification_trainer_fig}). The number of trials is user defined, set to 300 by default. Using at least 10 times the number of features is advised (for example, at least 300 iterations if 30 attributes are used). If the solution is not satisfying, more inliers can be selected, for example, in a region that the user identifies as misclassified with the current configuration. The training algorithm keeps the best weights found as initialization and carries on trying new weights by taking new inliers into account.
+
+\cgalFigureBegin{Classification_trainer_fig,classif_training.png}
+Example of evolution of training score. The purple curve is the score computed at the current iteration, green curve is the best score found so far.
+\cgalFigureEnd
+
+Tools to evaluate the reliability of the training are provided, namely:
+
+- [precision()](@ref CGAL::Classification::Trainer::precision) computes, for one label, the ratio of true positives over the total number of detected positives
+- [recall()](@ref CGAL::Classification::Trainer::recall) computes, for one label, the ratio of true positives over the total number of provided inliers of this label
+- [f1_score()](@ref CGAL::Classification::Trainer::f1_score) is the harmonic mean of precision and recall
+- [intersection_over_union()](@ref CGAL::Classification::Trainer::intersection_over_union) computes the ratio of true positives over the union of the detected positives and of the provided inliers
+- [accuracy()](@ref CGAL::Classification::Trainer::accuracy) computes the ratio of all true positives over the total number of provided inliers
+- [mean_f1_score()](@ref CGAL::Classification::Trainer::mean_f1_score)
+- [mean_intersection_over_union()](@ref CGAL::Classification::Trainer::mean_intersection_over_union)
+
+
+Note that the [Trainer](@ref CGAL::Classification::Trainer) class only uses documented public methods of the [Classifier](@ref CGAL::Classifier), [Label](@ref CGAL::Classification::Label) and [Feature](@ref CGAL::Classification::Feature_base) classes (for example, [set_weight()](@ref CGAL::Classification::Feature_base::set_weight), [set_feature_effect()](@ref CGAL::Classification::Label::set_feature_effect) or [energy_of()](@ref CGAL::Classifier::energy_of)). Users can therefore define their own classes for training using these methods to set up weights and effects and to estimate what are the best parameters.
+
+
+
+The example \ref Classification_example_training shows how to set inliers and run the training algorithm.
+
+\section Classification_examples Examples
+
+\subsection Classification_example_general General classification
+
+The following example:
+
+- reads an input file (LIDAR point set in PLY format)
+- computes useful structures from this input
+- computes segmentation features from the input and the precomputed structures
+- defines 3 labels (vegetation, ground and roof) along with the effects of features on them
+- classifies the point set using `CGAL::Classifier`
+- saves the result in a colored PLY format.
+
+\cgalExample{Classification/example_classifier.cpp}
+
+\subsection Classification_example_training Point set classification with training
+
+The following example:
+
+- reads a point set with a training set (embedded as a PLY feature _label_)
+- uses `CGAL::Point_set_classifier` to generate features on 5 scales
+- trains the algorithm using 800 trials
+- runs the algorithm using the graphcut regularization
+- saves the output to PLY.
+
+\cgalExample{Classification/example_point_set_classifier.cpp}
+
+
+\section Classification_history History
+
+This package is based on a research code by Florent Lafarge that was generalized, extended and packaged by Simon Giraudot.
+
+
+
+*/
+} /* namespace CGAL */

Classification/doc/Classification/Concepts/NeighborQuery.h

@@ -0,0 +1,35 @@
+namespace CGAL
+{
+
+/*!
+\ingroup PkgClassificationConcepts
+\cgalConcept
+
+Concept describing a neighbor query used for classification.
+
+\cgalHasModel `CGAL::Classification::Point_set_neighborhood::K_neighbor_query`
+\cgalHasModel `CGAL::Classification::Point_set_neighborhood::Range_neighbor_query`
+
+*/
+class NeighborQuery
+{
+public:
+
+  /*!
+    \brief Type of the data that is classified.
+  */
+  typedef unspecified_type value_type;
+
+  /*!
+
+    \brief Puts in `output` the indices of the neighbors of `query`.
+
+    \tparam OutputIterator An output iterator accepting `std::size_t`
+    values.
+   */
+  template <typename OutputIterator>
+  OutputIterator operator() (const value_type& query, OutputIterator output) const;
+
+};
+
+} // namespace CGAL

Classification/doc/Classification/Doxyfile.in

@@ -0,0 +1,7 @@
+@INCLUDE = ${CGAL_DOC_PACKAGE_DEFAULTS}
+PROJECT_NAME = "CGAL ${CGAL_DOC_VERSION} - Classification"
+INPUT        = ${CMAKE_SOURCE_DIR}/Classification/doc/Classification/ \
+               ${CMAKE_SOURCE_DIR}/Classification/include/CGAL/Classifier.h \
+               ${CMAKE_SOURCE_DIR}/Classification/include/CGAL/Point_set_classifier.h \
+               ${CMAKE_SOURCE_DIR}/Classification/include/CGAL/Classification/Feature/ \
+               ${CMAKE_SOURCE_DIR}/Classification/include/CGAL/Classification/

Classification/doc/Classification/PackageDescription.txt

@@ -0,0 +1,80 @@
+/*!
+\defgroup PkgClassification Classification Reference
+
+\defgroup PkgClassificationConcepts Concepts
+\ingroup PkgClassification
+
+\defgroup PkgClassificationDataStructures Data Structures
+\ingroup PkgClassification
+
+\defgroup PkgClassificationAttributes Predefined Attributes
+\ingroup PkgClassification
+
+
+\addtogroup PkgClassification
+
+\cgalPkgDescriptionBegin{Classification, PkgClassificationSummary}
+\cgalPkgPicture{data_classif.png}
+
+\cgalPkgSummaryBegin
+\cgalPkgAuthors{Simon Giraudot, Florent Lafarge}
+\cgalPkgDesc{This component implements an algorithm that classifies a data set into a user-defined set of labels (such as ground, vegetation, buildings, etc.). A flexible API is provided so that the user can classify any type of data, compute its own local features on the data set and define its own labels based on these features.}
+\cgalPkgManuals{Chapter_Classification, PkgClassification}
+\cgalPkgSummaryEnd
+
+\cgalPkgShortInfoBegin
+\cgalPkgSince{4.9}
+\cgalPkgBib{cgal:lm-clscm-12}
+\cgalPkgLicense{\ref licensesGPL "GPL"}
+\cgalPkgDemo{Operations on Polyhedra,polyhedron_3.zip}
+\cgalPkgShortInfoEnd
+
+\cgalPkgDescriptionEnd
+
+\cgalClassifedRefPages
+
+## Concepts ##
+
+- `CGAL::NeighborQuery`
+
+## Classification ##
+
+- `CGAL::Classifier<ItemRange, ItemMap, ConcurrencyTag>`
+- `CGAL::Point_set_classifier<Geom_traits, PointRange, PointMap, ConcurrencyTag, DiagonalizeTraits>`
+
+## Data Structures ##
+
+- `CGAL::Classification::Point_set_neighborhood<Geom_traits, PointRange, PointMap>`
+- `CGAL::Classification::Local_eigen_analysis<Geom_traits, PointRange, PointMap, DiagonalizeTraits>`
+- `CGAL::Classification::Planimetric_grid<Geom_traits, PointRange, PointMap>`
+- `CGAL::Classification::Trainer<PointRange, PointMap, ConcurrencyTag>`
+
+## Label ##
+
+- `CGAL::Classification::Label`
+- `CGAL::Classification::Label_handle`
+
+## Feature ##
+
+- `CGAL::Classification::Feature`
+- `CGAL::Classification::Feature_handle`
+
+## Predefined Features ##
+
+- `CGAL::Classification::Feature::Anisotropy<Geom_traits, PointRange, PointMap, DiagonalizeTraits>`
+- `CGAL::Classification::Feature::Distance_to_plane<Geom_traits, PointRange, PointMap, DiagonalizeTraits>`
+- `CGAL::Classification::Feature::Echo_scatter<Geom_traits, PointRange, PointMap, EchoMap>`
+- `CGAL::Classification::Feature::Eigentropy<Geom_traits, PointRange, PointMap, DiagonalizeTraits>`
+- `CGAL::Classification::Feature::Elevation<Geom_traits, PointRange, PointMap>`
+- `CGAL::Classification::Feature::Hsv<Geom_traits, PointRange, ColorMap>`
+- `CGAL::Classification::Feature::Linearity<Geom_traits, PointRange, PointMap, DiagonalizeTraits>`
+- `CGAL::Classification::Feature::Omnivariance<Geom_traits, PointRange, PointMap, DiagonalizeTraits>`
+- `CGAL::Classification::Feature::Planarity<Geom_traits, PointRange, PointMap, DiagonalizeTraits>`
+- `CGAL::Classification::Feature::Sphericity<Geom_traits, PointRange, PointMap, DiagonalizeTraits>`
+- `CGAL::Classification::Feature::Sum_eigenvalues<Geom_traits, PointRange, PointMap, DiagonalizeTraits>`
+- `CGAL::Classification::Feature::Surface_variation<Geom_traits, PointRange, PointMap, DiagonalizeTraits>`
+- `CGAL::Classification::Feature::Vertical_dispersion<Geom_traits, PointRange, PointMap>`
+- `CGAL::Classification::Feature::Verticality<Geom_traits, PointRange, PointMap, DiagonalizeTraits>`
+
+
+*/

Classification/doc/Classification/dependencies

@@ -0,0 +1,7 @@
+Manual
+Kernel_23
+STL_Extension
+Algebraic_foundations
+Circulator
+Stream_support
+Solver_interface

Classification/doc/Classification/examples.txt

@@ -0,0 +1,5 @@
+/*!
+\example Classification/example_classifier.cpp
+\example Classification/example_feature.cpp
+\example Classification/example_point_set_classifier.cpp
+*/

Classification/examples/Classification/CMakeLists.txt

@@ -0,0 +1,55 @@
+# Created by the script cgal_create_CMakeLists
+# This is the CMake script for compiling a set of CGAL applications.
+
+project( Classification )
+
+
+cmake_minimum_required(VERSION 2.8.11)
+
+# CGAL and its components
+find_package( CGAL QUIET COMPONENTS  )
+
+if ( NOT CGAL_FOUND )
+
+  message(STATUS "This project requires the CGAL library, and will not be compiled.")
+  return()  
+
+endif()
+
+# include helper file
+include( ${CGAL_USE_FILE} )
+
+
+# Boost and its components
+find_package( Boost REQUIRED )
+
+if ( NOT Boost_FOUND )
+
+  message(STATUS "This project requires the Boost library, and will not be compiled.")
+
+  return()  
+
+endif()
+
+find_package( TBB REQUIRED )
+if( TBB_FOUND )
+  include(${TBB_USE_FILE})
+  list(APPEND CGAL_3RD_PARTY_LIBRARIES ${TBB_LIBRARIES})
+endif()
+
+
+# include for local directory
+include_directories( BEFORE include )
+
+# include for local package
+include_directories( BEFORE ../../include )
+
+
+# Creating entries for all C++ files with "main" routine
+# ##########################################################
+
+include( CGAL_CreateSingleSourceCGALProgram )
+
+create_single_source_cgal_program( "example_classifier.cpp" )
+create_single_source_cgal_program( "example_point_set_classifier.cpp" )
+create_single_source_cgal_program( "example_feature.cpp" )

Classification/examples/Classification/example_classifier.cpp

@@ -0,0 +1,194 @@
+#include <cstdlib>
+#include <fstream>
+#include <iostream>
+#include <string>
+
+#include <CGAL/Simple_cartesian.h>
+#include <CGAL/Classifier.h>
+#include <CGAL/Classification/Point_set_neighborhood.h>
+#include <CGAL/Classification/Planimetric_grid.h>
+#include <CGAL/Classification/Feature_base.h>
+#include <CGAL/Classification/Feature/Eigen.h>
+#include <CGAL/Classification/Feature/Distance_to_plane.h>
+#include <CGAL/Classification/Feature/Vertical_dispersion.h>
+#include <CGAL/Classification/Feature/Elevation.h>
+
+#include <CGAL/IO/read_ply_points.h>
+
+#include <CGAL/Real_timer.h>
+
+typedef CGAL::Simple_cartesian<double> Kernel;
+typedef Kernel::Point_3 Point;
+typedef Kernel::Iso_cuboid_3 Iso_cuboid_3;
+typedef std::vector<Point> Point_range;
+typedef CGAL::Identity_property_map<Point> Pmap;
+
+typedef CGAL::Classifier<Point_range, Pmap> Classifier;
+
+typedef CGAL::Classification::Planimetric_grid<Kernel, Point_range, Pmap>             Planimetric_grid;
+typedef CGAL::Classification::Point_set_neighborhood<Kernel, Point_range, Pmap>       Neighborhood;
+typedef CGAL::Classification::Local_eigen_analysis<Kernel, Point_range, Pmap>         Local_eigen_analysis;
+
+typedef CGAL::Classification::Label_handle                                            Label_handle;
+typedef CGAL::Classification::Feature_handle                                          Feature_handle;
+
+typedef CGAL::Classification::Feature::Distance_to_plane<Kernel, Point_range, Pmap>   Distance_to_plane;
+typedef CGAL::Classification::Feature::Linearity<Kernel, Point_range, Pmap>           Linearity;
+typedef CGAL::Classification::Feature::Omnivariance<Kernel, Point_range, Pmap>        Omnivariance;
+typedef CGAL::Classification::Feature::Planarity<Kernel, Point_range, Pmap>           Planarity;
+typedef CGAL::Classification::Feature::Surface_variation<Kernel, Point_range, Pmap>   Surface_variation;
+typedef CGAL::Classification::Feature::Elevation<Kernel, Point_range, Pmap>           Elevation;
+typedef CGAL::Classification::Feature::Vertical_dispersion<Kernel, Point_range, Pmap> Dispersion;
+
+
+///////////////////////////////////////////////////////////////////
+//! [Analysis]
+
+int main (int argc, char** argv)
+{
+  std::string filename (argc > 1 ? argv[1] : "data/b9.ply");
+  std::ifstream in (filename.c_str());
+  std::vector<Point> pts;
+
+  std::cerr << "Reading input" << std::endl;
+  if (!in
+      || !(CGAL::read_ply_points (in, std::back_inserter (pts))))
+    {
+      std::cerr << "Error: cannot read " << filename << std::endl;
+      return EXIT_FAILURE;
+    }
+
+  double grid_resolution = 0.34;
+  double radius_neighbors = 1.7;
+  double radius_dtm = 15.0;
+
+  std::cerr << "Computing useful structures" << std::endl;
+
+  Iso_cuboid_3 bbox = CGAL::bounding_box (pts.begin(), pts.end());
+
+  Planimetric_grid grid (pts, Pmap(), bbox, grid_resolution);
+  Neighborhood neighborhood (pts, Pmap());
+  Local_eigen_analysis eigen (pts, Pmap(), neighborhood.k_neighbor_query(6));
+
+  Classifier classifier (pts, Pmap());
+
+  //! [Analysis]
+  ///////////////////////////////////////////////////////////////////
+  
+  ///////////////////////////////////////////////////////////////////
+  //! [Features]
+
+  std::cerr << "Computing features" << std::endl;
+  Feature_handle d2p = classifier.add_feature<Distance_to_plane> (Pmap(), eigen);
+  Feature_handle lin = classifier.add_feature<Linearity> (eigen);
+  Feature_handle omni = classifier.add_feature<Omnivariance> (eigen);
+  Feature_handle plan = classifier.add_feature<Planarity> (eigen);
+  Feature_handle surf = classifier.add_feature<Surface_variation> (eigen);
+  Feature_handle disp = classifier.add_feature<Dispersion> (Pmap(), grid,
+                                                                grid_resolution,
+                                                                radius_neighbors);
+  Feature_handle elev = classifier.add_feature<Elevation> (Pmap(), grid,
+                                                               grid_resolution,
+                                                               radius_dtm);
+  
+  std::cerr << "Setting weights" << std::endl;
+  d2p->set_weight(6.75e-2);
+  lin->set_weight(1.19);
+  omni->set_weight(1.34e-1);
+  plan->set_weight(7.32e-1);
+  surf->set_weight(1.36e-1);
+  disp->set_weight(5.45e-1);
+  elev->set_weight(1.47e1);
+
+  
+  //! [Features]
+  ///////////////////////////////////////////////////////////////////
+
+  ///////////////////////////////////////////////////////////////////
+  //! [Labels]
+  
+  std::cerr << "Setting up labels" << std::endl;
+  
+  // Create label and define how features affect them
+  Label_handle ground = classifier.add_label ("ground");
+  ground->set_feature_effect (d2p,  CGAL::Classification::Feature::NEUTRAL);
+  ground->set_feature_effect (lin,  CGAL::Classification::Feature::PENALIZING);
+  ground->set_feature_effect (omni, CGAL::Classification::Feature::NEUTRAL);
+  ground->set_feature_effect (plan, CGAL::Classification::Feature::FAVORING);
+  ground->set_feature_effect (surf, CGAL::Classification::Feature::PENALIZING);
+  ground->set_feature_effect (disp, CGAL::Classification::Feature::NEUTRAL);
+  ground->set_feature_effect (elev, CGAL::Classification::Feature::PENALIZING);
+
+  Label_handle vege = classifier.add_label ("vegetation");
+  vege->set_feature_effect (d2p,  CGAL::Classification::Feature::FAVORING);
+  vege->set_feature_effect (lin,  CGAL::Classification::Feature::NEUTRAL);
+  vege->set_feature_effect (omni, CGAL::Classification::Feature::FAVORING);
+  vege->set_feature_effect (plan, CGAL::Classification::Feature::NEUTRAL);
+  vege->set_feature_effect (surf, CGAL::Classification::Feature::NEUTRAL);
+  vege->set_feature_effect (disp, CGAL::Classification::Feature::FAVORING);
+  vege->set_feature_effect (elev, CGAL::Classification::Feature::NEUTRAL);
+  
+  Label_handle roof = classifier.add_label ("roof");
+  roof->set_feature_effect (d2p,  CGAL::Classification::Feature::NEUTRAL);
+  roof->set_feature_effect (lin,  CGAL::Classification::Feature::PENALIZING);
+  roof->set_feature_effect (omni, CGAL::Classification::Feature::FAVORING);
+  roof->set_feature_effect (plan, CGAL::Classification::Feature::FAVORING);
+  roof->set_feature_effect (surf, CGAL::Classification::Feature::PENALIZING);
+  roof->set_feature_effect (disp, CGAL::Classification::Feature::NEUTRAL);
+  roof->set_feature_effect (elev, CGAL::Classification::Feature::FAVORING);
+
+  //! [Labels]
+  ///////////////////////////////////////////////////////////////////
+
+  // Run classification
+  CGAL::Real_timer t;
+  t.start();
+  classifier.run ();
+  t.stop();
+  std::cerr << "Raw classification performed in " << t.time() << " second(s)" << std::endl;
+  t.reset();
+  t.start();
+  classifier.run_with_local_smoothing (neighborhood.range_neighbor_query(radius_neighbors));
+  t.stop();
+  std::cerr << "Classification with local smoothing performed in " << t.time() << " second(s)" << std::endl;
+  t.reset();
+  t.start();
+  classifier.run_with_graphcut (neighborhood.k_neighbor_query(12), 0.2);
+  t.stop();
+  std::cerr << "Classification with graphcut performed in " << t.time() << " second(s)" << std::endl;
+  
+  // Save the output in a colored PLY format
+
+  std::ofstream f ("classification.ply");
+  f << "ply" << std::endl
+    << "format ascii 1.0" << std::endl
+    << "element vertex " << pts.size() << std::endl
+    << "property float x" << std::endl
+    << "property float y" << std::endl
+    << "property float z" << std::endl
+    << "property uchar red" << std::endl
+    << "property uchar green" << std::endl
+    << "property uchar blue" << std::endl
+    << "end_header" << std::endl;
+  
+  for (std::size_t i = 0; i < pts.size(); ++ i)
+    {
+      f << pts[i] << " ";
+      
+      Label_handle label = classifier.label_of (i);
+      if (label == ground)
+        f << "245 180 0" << std::endl;
+      else if (label == vege)
+        f << "0 255 27" << std::endl;
+      else if (label == roof)
+        f << "255 0 170" << std::endl;
+      else
+        {
+          f << "0 0 0" << std::endl;
+          std::cerr << "Error: unknown classification label" << std::endl;
+        }
+    }
+  
+  std::cerr << "All done" << std::endl;
+  return EXIT_SUCCESS;
+}

Classification/examples/Classification/example_feature.cpp

@@ -0,0 +1,99 @@
+#include <cstdlib>
+#include <fstream>
+#include <iostream>
+#include <string>
+
+#include <CGAL/Simple_cartesian.h>
+#include <CGAL/Classifier.h>
+#include <CGAL/Classification/Point_set_neighborhood.h>
+#include <CGAL/Classification/Planimetric_grid.h>
+#include <CGAL/Classification/Feature_base.h>
+#include <CGAL/Classification/Feature/Eigen.h>
+
+#include <CGAL/IO/read_ply_points.h>
+
+typedef CGAL::Simple_cartesian<double> Kernel;
+typedef Kernel::Point_3 Point;
+typedef Kernel::Iso_cuboid_3 Iso_cuboid_3;
+typedef std::vector<Point> Point_range;
+typedef CGAL::Identity_property_map<Point> Pmap;
+
+typedef CGAL::Classifier<Point_range, Pmap> Classifier;
+
+typedef CGAL::Classification::Point_set_neighborhood<Kernel, Point_range, Pmap>       Neighborhood;
+typedef CGAL::Classification::Local_eigen_analysis<Kernel, Point_range, Pmap>         Local_eigen_analysis;
+
+typedef CGAL::Classification::Label_handle                                            Label_handle;
+typedef CGAL::Classification::Feature_handle                                          Feature_handle;
+
+typedef CGAL::Classification::Feature::Sphericity<Kernel, Point_range, Pmap>          Sphericity;
+
+
+// User-defined feature that identifies a specific area of the 3D
+// space. This feature takes value 1 for points that lie inside the
+// area and 0 for the others.
+class My_feature : public CGAL::Classification::Feature_base
+{
+  const Point_range& range;
+  double xmin, xmax, ymin, ymax;
+public:
+  My_feature (const Point_range& range, // constructor should start with item range
+              double xmin, double xmax, double ymin, double ymax)
+    : range (range), xmin(xmin), xmax(xmax), ymin(ymin), ymax(ymax)
+  { }
+
+  double value (std::size_t pt_index)
+  {
+    if (xmin < range[pt_index].x() && range[pt_index].x() < xmax &&
+        ymin < range[pt_index].y() && range[pt_index].y() < ymax)
+      return 1.;
+    else
+      return 0.;
+  }
+    
+};
+
+int main (int argc, char** argv)
+{
+  std::string filename (argc > 1 ? argv[1] : "data/b9.ply");
+  std::ifstream in (filename.c_str());
+  std::vector<Point> pts;
+
+  std::cerr << "Reading input" << std::endl;
+  if (!in
+      || !(CGAL::read_ply_points (in, std::back_inserter (pts))))
+    {
+      std::cerr << "Error: cannot read " << filename << std::endl;
+      return EXIT_FAILURE;
+    }
+
+  Neighborhood neighborhood (pts, Pmap());
+  Local_eigen_analysis eigen (pts, Pmap(), neighborhood.k_neighbor_query(6));
+
+  Classifier classifier (pts, Pmap());
+
+  std::cerr << "Computing features" << std::endl;
+  Feature_handle sphericity = classifier.add_feature<Sphericity> (eigen);
+
+  // Feature that identifies points whose x coordinate is between -20
+  // and 20 and whose y coordinate is between -15 and 15
+  Feature_handle my_feature = classifier.add_feature<My_feature> (-20., 20., -15., 15.);
+                                                     
+  std::cerr << "Setting weights" << std::endl;
+  sphericity->set_weight(0.5);
+  my_feature->set_weight(0.25);
+
+  std::cerr << "Setting up labels" << std::endl;
+  Label_handle a = classifier.add_label ("label_A");
+  a->set_feature_effect (sphericity,  CGAL::Classification::Feature::FAVORING);
+  a->set_feature_effect (my_feature,  CGAL::Classification::Feature::FAVORING);
+
+  Label_handle b = classifier.add_label ("label_B");
+  b->set_feature_effect (sphericity,  CGAL::Classification::Feature::PENALIZING);
+  b->set_feature_effect (my_feature,  CGAL::Classification::Feature::PENALIZING);
+
+  classifier.run_with_graphcut (neighborhood.k_neighbor_query(12), 0.2);
+
+  std::cerr << "All done" << std::endl;
+  return EXIT_SUCCESS;
+}

Classification/examples/Classification/example_point_set_classifier.cpp

@@ -0,0 +1,184 @@
+#include <cstdlib>
+#include <fstream>
+#include <iostream>
+#include <string>
+
+//#define CGAL_CLASSIFICATION_VERBOSE
+
+#include <CGAL/Simple_cartesian.h>
+#include <CGAL/Point_set_classifier.h>
+#include <CGAL/Classification/Trainer.h>
+#include <CGAL/IO/read_ply_points.h>
+
+#include <CGAL/Real_timer.h>
+
+typedef CGAL::Simple_cartesian<double> Kernel;
+typedef Kernel::Point_3 Point;
+typedef Kernel::Iso_cuboid_3 Iso_cuboid_3;
+typedef std::vector<Point> Point_range;
+typedef CGAL::Identity_property_map<Point> Pmap;
+
+typedef CGAL::Point_set_classifier<Kernel, Point_range, Pmap> Point_set_classifier;
+typedef CGAL::Classification::Trainer<Point_range, Pmap> Trainer;
+
+/*
+ This interpreter is used to read a PLY input that contains training
+ attributes (with the PLY "label" property).
+*/
+class My_ply_interpreter
+{
+  std::vector<Point>& points;
+  std::vector<int>& labels;
+    
+public:
+  My_ply_interpreter (std::vector<Point>& points,
+                      std::vector<int>& labels)
+    : points (points), labels (labels)
+  { }
+
+  // Init and test if input file contains the right properties
+  bool is_applicable (CGAL::Ply_reader& reader)
+  {
+    return reader.does_tag_exist<double> ("x")
+      && reader.does_tag_exist<double> ("y")
+      && reader.does_tag_exist<double> ("z")
+      && reader.does_tag_exist<int> ("label");
+  }
+
+  // Describes how to process one line (= one point object)
+  void process_line (CGAL::Ply_reader& reader)
+  {
+    double x = 0., y = 0., z = 0.;
+    int l = 0;
+
+    reader.assign (x, "x");
+    reader.assign (y, "y");
+    reader.assign (z, "z");
+    reader.assign (l, "label");
+
+    points.push_back (Point (x, y, z));
+    labels.push_back(l);
+  }
+
+};
+
+
+int main (int argc, char** argv)
+{
+  std::string filename (argc > 1 ? argv[1] : "data/b9_training.ply");
+  std::ifstream in (filename.c_str());
+  std::vector<Point> pts;
+  std::vector<int> labels;
+
+  std::cerr << "Reading input" << std::endl;
+  My_ply_interpreter interpreter (pts, labels);
+  if (!in
+      || !(CGAL::read_ply_custom_points (in, interpreter, Kernel())))
+    {
+      std::cerr << "Error: cannot read " << filename << std::endl;
+      return EXIT_FAILURE;
+    }
+
+  Point_set_classifier psc (pts, Pmap());
+  
+  std::cerr << "Generating features" << std::endl;
+  CGAL::Real_timer t;
+  t.start();
+  psc.generate_features (5); // Using 5 scales
+  t.stop();
+  std::cerr << "Done in " << t.time() << " second(s)" << std::endl;
+  
+  // Add types to PSC
+  CGAL::Classification::Label_handle ground
+    = psc.add_label ("ground");
+  CGAL::Classification::Label_handle vege
+    = psc.add_label ("vegetation");
+  CGAL::Classification::Label_handle roof
+    = psc.add_label ("roof");
+  CGAL::Classification::Label_handle facade
+    = psc.add_label ("facade");
+  
+  Trainer trainer (psc);
+
+  // Set training sets
+  std::size_t nb_inliers = 0;
+  for (std::size_t i = 0; i < labels.size(); ++ i)
+    {
+      switch (labels[i])
+        {
+        case 0:
+          trainer.set_inlier(vege, i);
+          ++ nb_inliers;
+          break;
+        case 1:
+          trainer.set_inlier(ground, i);
+          ++ nb_inliers;
+          break;
+        case 2:
+          trainer.set_inlier(roof, i);
+          ++ nb_inliers;
+          break;
+        case 3:
+          trainer.set_inlier(facade, i);
+          ++ nb_inliers;
+          break;
+        default:
+          break;
+        }
+    }
+
+  std::cerr << "Training using " << nb_inliers << " inliers" << std::endl;
+  t.reset();
+  t.start();
+  trainer.train (400); // 800 trials
+  t.stop();
+  std::cerr << "Done in " << t.time() << " second(s)" << std::endl;
+
+  std::cerr << "Precision, recall, F1 scores and IoU:" << std::endl;
+  for (std::size_t i = 0; i < psc.number_of_labels(); ++ i)
+    {
+      std::cerr << " * " << psc.label(i)->name() << ": "
+                << trainer.precision(psc.label(i)) << " ; "
+                << trainer.recall(psc.label(i)) << " ; "
+                << trainer.f1_score(psc.label(i)) << " ; "
+                << trainer.intersection_over_union(psc.label(i)) << std::endl;
+    }
+
+  std::cerr << "Accuracy = " << trainer.accuracy() << std::endl
+            << "Mean F1 score = " << trainer.mean_f1_score() << std::endl
+            << "Mean IoU = " << trainer.mean_intersection_over_union() << std::endl;
+  
+  t.reset();
+  t.start();
+  psc.run_with_graphcut (psc.neighborhood().k_neighbor_query(12), 0.5);
+  t.stop();
+  std::cerr << "One graphcut done in " << t.time() << " second(s)" << std::endl;
+  
+  // Save the output in a colored PLY format
+  {
+    std::ofstream f ("classification_one.ply");
+    f.precision(18);
+    psc.write_classification_to_ply (f);
+  }
+
+  t.reset();
+  t.start();
+  psc.run_with_graphcut (psc.neighborhood().k_neighbor_query(12), 0.5, 30);
+  t.stop();
+  std::cerr << std::size_t(pts.size() / 25000) << " graphcuts done in " << t.time() << " second(s)" << std::endl;
+  
+  // Save the output in a colored PLY format
+  {
+    std::ofstream f ("classification_several.ply");
+    f.precision(18);
+    psc.write_classification_to_ply (f);
+  }
+  
+  /// Save the configuration to be able to reload it later
+  std::ofstream fconfig ("config.xml");
+  psc.save_configuration (fconfig);
+  
+  std::cerr << "All done" << std::endl;
+  
+  return EXIT_SUCCESS;
+}

Classification/include/CGAL/Classification/Color.h

@@ -0,0 +1,136 @@
+// Copyright (c) 2017 GeometryFactory Sarl (France).
+// All rights reserved.
+//
+// This file is part of CGAL (www.cgal.org).
+// You can redistribute it and/or modify it under the terms of the GNU
+// General Public License as published by the Free Software Foundation,
+// either version 3 of the License, or (at your option) any later version.
+//
+// Licensees holding a valid commercial license may use this file in
+// accordance with the commercial license agreement provided with the software.
+//
+// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
+// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
+//
+// $URL$
+// $Id$
+//
+// Author(s)     : Simon Giraudot
+
+#ifndef CGAL_CLASSIFICATION_COLOR_H
+#define CGAL_CLASSIFICATION_COLOR_H
+
+#include <CGAL/array.h>
+
+namespace CGAL {
+namespace Classification {
+
+  /*!
+    \ingroup PkgClassificationColor
+
+    Color described in red/green/blue space. Each component is stored
+    as an unsigned char ranging from 0 (no color) to 255 (
+   */  
+typedef CGAL::cpp11::array<unsigned char, 3> RGB_Color;
+  /*!
+    \ingroup PkgClassificationColor
+
+    Color described in hue/saturation/value space. Each component is stored
+    as a double:
+
+    - `hue` ranges from 0° to 360° (corresponding to the color tint)
+    - `saturation` ranges from 0.0 (gray) to 100.0 (full saturation)
+    - `value` ranges from 0.0 (black) to 100.0 (white)
+   */  
+typedef CGAL::cpp11::array<double, 3> HSV_Color;
+
+
+  /// \cond SKIP_IN_MANUAL
+inline HSV_Color rgb_to_hsv (const RGB_Color& c)
+{
+  double r = (double)(c[0]) / 255.;
+  double g = (double)(c[1]) / 255.;
+  double b = (double)(c[2]) / 255.;
+  double Cmax = (std::max) (r, (std::max) (g, b));
+  double Cmin = (std::min) (r, (std::min) (g, b));
+  double delta = Cmax - Cmin;
+  double H = 0.;
+  
+  if (delta != 0.)
+    {
+      if (Cmax == r)
+        H = 60. * ((g - b) / delta);
+      else if (Cmax == g)
+        H = 60. * (((b - r) / delta) + 2.);
+      else
+        H = 60. * (((r - g) / delta) + 4.);
+    }
+  if (H < 0.) H += 360.;
+  double S = (Cmax == 0. ? 0. : 100. * (delta / Cmax));
+  double V = 100. * Cmax;
+  HSV_Color out = {{ H, S, V }};
+  return out;
+}
+
+inline RGB_Color hsv_to_rgb (const HSV_Color& c)
+{
+  double h = c[0];
+  double s = c[1];
+  double v = c[2];
+  
+  s /= 100.;
+  v /= 100.;
+  double C = v*s;
+  int hh = (int)(h/60.);
+  double X = C * (1-CGAL::abs (hh % 2 - 1));
+  double r = 0, g = 0, b = 0;
+  
+  if( hh>=0 && hh<1 )
+    {
+      r = C;
+      g = X;
+    }
+  else if( hh>=1 && hh<2 )
+    {
+      r = X;
+      g = C;
+    }
+  else if( hh>=2 && hh<3 )
+    {
+      g = C;
+      b = X;
+    }
+  else if( hh>=3 && hh<4 )
+    {
+      g = X;
+      b = C;
+    }
+  else if( hh>=4 && hh<5 )
+    {
+      r = X;
+      b = C;
+    }
+  else
+    {
+      r = C;
+      b = X;
+    }
+  double m = v-C;
+  r += m;
+  g += m;
+  b += m;
+  r *= 255.0;
+  g *= 255.0;
+  b *= 255.0;
+
+  RGB_Color out = {{ (unsigned char)r, (unsigned char)g, (unsigned char)b }};
+  return out;
+}
+  /// \endcond
+
+} // namespace Classification
+} // namespace CGAL
+
+
+
+#endif // CGAL_CLASSIFICATION_COLOR_H

Classification/include/CGAL/Classification/Feature/Distance_to_plane.h

@@ -0,0 +1,115 @@
+// Copyright (c) 2012 INRIA Sophia-Antipolis (France).
+// Copyright (c) 2017 GeometryFactory Sarl (France).
+// All rights reserved.
+//
+// This file is part of CGAL (www.cgal.org).
+// You can redistribute it and/or modify it under the terms of the GNU
+// General Public License as published by the Free Software Foundation,
+// either version 3 of the License, or (at your option) any later version.
+//
+// Licensees holding a valid commercial license may use this file in
+// accordance with the commercial license agreement provided with the software.
+//
+// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
+// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
+//
+// $URL$
+// $Id$
+//
+// Author(s)     : Simon Giraudot, Florent Lafarge
+
+#ifndef CGAL_CLASSIFICATION_FEATURE_DISTANCE_TO_PLANE_H
+#define CGAL_CLASSIFICATION_FEATURE_DISTANCE_TO_PLANE_H
+
+#include <vector>
+
+#include <CGAL/Classifier.h>
+
+namespace CGAL {
+
+namespace Classification {
+
+namespace Feature {
+
+  /*!
+    \ingroup PkgClassificationFeatures
+
+    %Feature based on local distance to a fitted
+    plane. Characterizing a level of non-planarity can help to
+    identify noisy parts of the input such as vegetation. This
+    feature computes the distance of a point to a locally fitted
+    plane.
+    
+    \tparam Geom_traits model of \cgal Kernel.
+    \tparam PointRange model of `ConstRange`. Its iterator type
+    is `RandomAccessIterator`.
+    \tparam PointMap model of `ReadablePropertyMap` whose key
+    type is the value type of the iterator of `PointRange` and value type
+    is `Geom_traits::Point_3`.
+    \tparam DiagonalizeTraits model of `DiagonalizeTraits` used
+    for matrix diagonalization.
+  */
+template <typename Geom_traits, typename PointRange, typename PointMap,
+          typename DiagonalizeTraits = CGAL::Default_diagonalize_traits<double,3> >
+class Distance_to_plane : public Feature_base
+{
+  typedef Classification::Local_eigen_analysis<Geom_traits, PointRange,
+                                               PointMap, DiagonalizeTraits> Local_eigen_analysis;
+#ifdef CGAL_CLASSIFICATION_PRECOMPUTE_FEATURES
+  std::vector<double> distance_to_plane_feature;
+#else
+  const PointRange& input;
+  PointMap point_map;
+  const Local_eigen_analysis& eigen;
+#endif
+  
+public:
+  /*!
+    \brief Constructs the feature.
+
+    \param input input range.
+    \param point_map property map to access the input points.
+    \param eigen class with precomputed eigenvectors and eigenvalues.
+  */
+  Distance_to_plane (const PointRange& input,
+                     PointMap point_map,
+                     const Local_eigen_analysis& eigen)
+#ifndef CGAL_CLASSIFICATION_PRECOMPUTE_FEATURES
+    : input(input), point_map(point_map), eigen(eigen)
+#endif
+  {
+    this->set_weight(1.);
+#ifndef CGAL_CLASSIFICATION_PRECOMPUTE_FEATURES
+    std::vector<double> distance_to_plane_feature;
+#endif
+    
+    for(std::size_t i = 0; i < input.size(); i++)
+      distance_to_plane_feature.push_back
+        (CGAL::sqrt (CGAL::squared_distance (get(point_map, *(input.begin()+i)), eigen.plane(i))));
+    
+    this->compute_mean_max (distance_to_plane_feature, this->mean, this->max);
+    //    max *= 2;
+  }
+
+  /// \cond SKIP_IN_MANUAL
+  virtual double value (std::size_t pt_index)
+  {
+#ifdef CGAL_CLASSIFICATION_PRECOMPUTE_FEATURES
+    return distance_to_plane_feature[pt_index];
+#else
+    return CGAL::sqrt (CGAL::squared_distance
+                       (get(point_map, *(input.begin()+pt_index)), eigen.plane(pt_index)));
+#endif
+  }
+
+  virtual std::string name() { return "distance_to_plane"; }
+  /// \endcond
+};
+
+} // namespace Feature
+
+} // namespace Classification
+  
+} // namespace CGAL
+
+#endif // CGAL_CLASSIFICATION_FEATURE_DISTANCE_TO_PLANE_H

Classification/include/CGAL/Classification/Feature/Echo_scatter.h

@@ -0,0 +1,152 @@
+// Copyright (c) 2012 INRIA Sophia-Antipolis (France).
+// Copyright (c) 2017 GeometryFactory Sarl (France).
+// All rights reserved.
+//
+// This file is part of CGAL (www.cgal.org).
+// You can redistribute it and/or modify it under the terms of the GNU
+// General Public License as published by the Free Software Foundation,
+// either version 3 of the License, or (at your option) any later version.
+//
+// Licensees holding a valid commercial license may use this file in
+// accordance with the commercial license agreement provided with the software.
+//
+// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
+// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
+//
+// $URL$
+// $Id$
+//
+// Author(s)     : Simon Giraudot, Florent Lafarge
+
+#ifndef CGAL_CLASSIFICATION_FEATURE_ECHO_SCATTER_H
+#define CGAL_CLASSIFICATION_FEATURE_ECHO_SCATTER_H
+
+#include <vector>
+
+
+namespace CGAL {
+
+namespace Classification {
+
+namespace Feature {
+  
+  /*!
+    \ingroup PkgClassificationFeatures
+
+    %Feature based on echo scatter. The number of returns (echo
+    number) is a useful information provided by most LIDAR sensors. It
+    can help to identify trees.
+
+    \tparam Geom_traits model of \cgal Kernel.
+    \tparam PointRange model of `ConstRange`. Its iterator type
+    is `RandomAccessIterator`.
+    \tparam PointMap model of `ReadablePropertyMap` whose key
+    type is the value type of the iterator of `PointRange` and value type
+    is `Geom_traits::Point_3`.
+    \tparam EchoMap model of `ReadablePropertyMap` whose key
+    type is the value type of the iterator of `PointRange` and value type
+    is `std::size_t`.
+  */
+template <typename Geom_traits, typename PointRange, typename PointMap, typename EchoMap>
+class Echo_scatter : public Feature_base
+{
+public:
+  typedef Classification::Planimetric_grid<Geom_traits, PointRange, PointMap> Grid;
+private:  
+  typedef Classification::Image<float> Image_float;
+
+  std::vector<double> echo_scatter;
+  
+public:
+  /*!
+    \brief Constructs the feature.
+
+    \param input input range.
+    \param echo_map property map to access the echo values of the input points.
+    \param grid precomputed `Planimetric_grid`.
+    \param grid_resolution resolution of the planimetric grid.
+    \param radius_neighbors radius of local neighborhoods.
+  */
+  Echo_scatter (const PointRange& input,
+                EchoMap echo_map,
+                const Grid& grid,
+                const double grid_resolution,
+                double radius_neighbors = 1.)
+  {
+    this->set_weight(1.);
+    Image_float Scatter(grid.width(), grid.height());
+    for (std::size_t j = 0; j < grid.height(); j++)
+      for (std::size_t i = 0; i < grid.width(); i++)
+        Scatter(i,j)=0;
+
+    std::size_t square = (std::size_t)(0.5 * radius_neighbors / grid_resolution) + 1;
+
+    for (std::size_t j = 0; j < grid.height(); j++){	
+      for (std::size_t i = 0; i < grid.width(); i++){
+						
+        if(grid.mask(i,j)){
+
+          std::size_t squareXmin = (i < square ? 0 : i - square);
+          std::size_t squareXmax = (std::min) (grid.width()-1, i + square);
+          std::size_t squareYmin = (j < square ? 0 : j - square);
+          std::size_t squareYmax = (std::min) (grid.height()-1, j + square);
+			
+          std::size_t NB_echo_sup=0;
+          std::size_t NB_echo_total=0;
+
+          for(std::size_t k = squareXmin; k <= squareXmax; k++){
+            for(std::size_t l = squareYmin; l <= squareYmax; l++){
+									
+              if(CGAL::sqrt(pow((double)k-i,2)+pow((double)l-j,2))<=(double)0.5*radius_neighbors/grid_resolution){
+										
+                if(grid.indices(k,l).size()>0){
+									
+                  for(std::size_t t=0; t<grid.indices(k,l).size();t++){
+												
+                    std::size_t ip = grid.indices(k,l)[t]; 
+                    if(get(echo_map, *(input.begin()+ip)) > 1)
+                      NB_echo_sup++;
+                  }
+									
+                  NB_echo_total=NB_echo_total+grid.indices(k,l).size();
+									
+                }
+							
+              }
+						
+            }
+					
+          }
+					
+          Scatter(i,j)=(float)NB_echo_sup/NB_echo_total;
+				
+        }
+			
+      }
+		
+    }
+    for(std::size_t i = 0; i < input.size(); i++){
+      std::size_t I= grid.x(i);
+      std::size_t J= grid.y(i);
+      echo_scatter.push_back((double)Scatter(I,J));
+    }
+    this->compute_mean_max (echo_scatter, this->mean, this->max);
+  }
+
+  /// \cond SKIP_IN_MANUAL
+  virtual double value (std::size_t pt_index)
+  {
+    return echo_scatter[pt_index];
+  }
+
+  virtual std::string name() { return "echo_scatter"; }
+  /// \endcond
+};
+
+} // namespace Feature
+
+} // namespace Classification
+  
+} // namespace CGAL
+
+#endif // CGAL_CLASSIFICATION_FEATURE_ECHO_SCATTER_H

Classification/include/CGAL/Classification/Feature/Effect.h

@@ -0,0 +1,42 @@
+// Copyright (c) 2017 GeometryFactory Sarl (France).
+// All rights reserved.
+//
+// This file is part of CGAL (www.cgal.org).
+// You can redistribute it and/or modify it under the terms of the GNU
+// General Public License as published by the Free Software Foundation,
+// either version 3 of the License, or (at your option) any later version.
+//
+// Licensees holding a valid commercial license may use this file in
+// accordance with the commercial license agreement provided with the software.
+//
+// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
+// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
+//
+// $URL$
+// $Id$
+//
+// Author(s)     : Simon Giraudot
+
+#ifndef CGAL_CLASSIFICATION_FEATURE_EFFECT_H
+#define CGAL_CLASSIFICATION_FEATURE_EFFECT_H
+
+namespace CGAL {
+
+namespace Classification {
+
+namespace Feature {
+
+enum Effect /// Defines the effect of an feature on a type.
+    {
+      FAVORING = 0, ///< High values of the feature favor this type
+      NEUTRAL = 1, ///< The feature has no effect on this type
+      PENALIZING = 2 ///< Low values of the feature favor this type
+    };
+  
+} // Feature
+
+} // Classification
+
+} // CGAL
+
+#endif

Classification/include/CGAL/Classification/Feature/Eigen.h

@@ -0,0 +1,504 @@
+// Copyright (c) 2017 GeometryFactory Sarl (France).
+// All rights reserved.
+//
+// This file is part of CGAL (www.cgal.org).
+// You can redistribute it and/or modify it under the terms of the GNU
+// General Public License as published by the Free Software Foundation,
+// either version 3 of the License, or (at your option) any later version.
+//
+// Licensees holding a valid commercial license may use this file in
+// accordance with the commercial license agreement provided with the software.
+//
+// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
+// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
+//
+// $URL$
+// $Id$
+//
+// Author(s)     : Simon Giraudot
+
+#ifndef CGAL_CLASSIFICATION_FEATURES_EIGEN_H
+#define CGAL_CLASSIFICATION_FEATURES_EIGEN_H
+
+#include <vector>
+
+#include <CGAL/Classification/Local_eigen_analysis.h>
+
+namespace CGAL {
+
+namespace Classification {
+
+namespace Feature {
+
+template <typename Geom_traits, typename PointRange, typename PointMap,
+          typename DiagonalizeTraits = CGAL::Default_diagonalize_traits<double,3> >
+class Eigen_feature : public Feature_base
+{
+protected:
+  typedef Classification::Local_eigen_analysis<Geom_traits, PointRange,
+                                               PointMap, DiagonalizeTraits> Local_eigen_analysis;
+
+#ifdef CGAL_CLASSIFICATION_PRECOMPUTE_FEATURES
+  std::vector<double> attrib;
+#else
+  const Local_eigen_analysis& eigen;
+#endif
+  
+public:
+  Eigen_feature (const PointRange&,
+                   const Local_eigen_analysis& eigen)
+#ifndef CGAL_CLASSIFICATION_PRECOMPUTE_FEATURES
+    : eigen (eigen)
+#endif
+  {
+  }
+
+  virtual void init (const PointRange& input, const Local_eigen_analysis& eigen)
+  {
+    this->set_weight(1.);
+#ifndef CGAL_CLASSIFICATION_PRECOMPUTE_FEATURES
+    std::vector<double> attrib;
+#endif
+    attrib.reserve (input.size());
+    for (std::size_t i = 0; i < input.size(); ++ i)
+      attrib.push_back (get_value (eigen, i));
+
+    this->compute_mean_max (attrib, mean, this->max);
+  }
+  
+  virtual double get_value (const Local_eigen_analysis& eigen, std::size_t i) = 0;
+  virtual double value (std::size_t pt_index)
+  {
+#ifdef CGAL_CLASSIFICATION_PRECOMPUTE_FEATURES
+    return attrib[pt_index];
+#else
+    return get_value(eigen, pt_index);
+#endif
+  }
+  virtual std::string name() { return "eigen_feature"; }
+
+};
+  
+  /*!
+    \ingroup PkgClassificationFeatures
+
+    %Feature based on the eigenvalues of the covariance matrix of a
+    local neighborhood. Linearity is defined, for the 3 eigenvalues
+    \f$\lambda_1 \ge \lambda_2 \ge \lambda_3 \ge 0\f$, as:
+
+    \f[
+    \frac{\lambda_1 - \lambda_2}{\lambda_1}
+    \f]
+    \tparam Geom_traits model of \cgal Kernel.
+    \tparam PointRange model of `ConstRange`. Its iterator type
+    is `RandomAccessIterator`.
+    \tparam PointMap model of `ReadablePropertyMap` whose key
+    type is the value type of the iterator of `PointRange` and value type
+    is `Geom_traits::Point_3`.
+    \tparam DiagonalizeTraits model of `DiagonalizeTraits` used
+    for matrix diagonalization.
+  */
+template <typename Geom_traits, typename PointRange, typename PointMap,
+          typename DiagonalizeTraits = CGAL::Default_diagonalize_traits<double,3> >
+class Linearity : public Eigen_feature<Geom_traits, PointRange, PointMap, DiagonalizeTraits>
+{
+  typedef Classification::Local_eigen_analysis<Geom_traits, PointRange,
+                                               PointMap, DiagonalizeTraits> Local_eigen_analysis;
+  typedef Eigen_feature<Geom_traits, PointRange, PointMap, DiagonalizeTraits> Base;
+public:
+  /*!
+    Constructs the feature.
+
+    \param input input range.
+    \param eigen class with precomputed eigenvectors and eigenvalues.
+  */
+  Linearity (const PointRange& input,
+             const Local_eigen_analysis& eigen) : Base (input, eigen)
+  {
+    this->init(input, eigen);
+  }
+
+  /// \cond SKIP_IN_MANUAL
+  virtual double get_value (const Local_eigen_analysis& eigen, std::size_t i)
+  {
+    const typename Local_eigen_analysis::Eigenvalues& ev = eigen.eigenvalue(i);
+    if (ev[2] < 1e-15)
+      return 0.;
+    else
+      return ((ev[2] - ev[1]) / ev[2]);
+  }
+  virtual std::string name() { return "linearity"; }
+  /// \endcond
+};
+
+  /*!
+    \ingroup PkgClassificationFeatures
+
+    %Feature based on the eigenvalues of the covariance matrix of a
+    local neighborhood. Planarity is defined, for the 3 eigenvalues
+    \f$\lambda_1 \ge \lambda_2 \ge \lambda_3 \ge 0\f$, as:
+
+    \f[
+    \frac{\lambda_2 - \lambda_3}{\lambda_1}
+    \f]
+    \tparam Geom_traits model of \cgal Kernel.
+    \tparam PointRange model of `ConstRange`. Its iterator type
+    is `RandomAccessIterator`.
+    \tparam PointMap model of `ReadablePropertyMap` whose key
+    type is the value type of the iterator of `PointRange` and value type
+    is `Geom_traits::Point_3`.
+    \tparam DiagonalizeTraits model of `DiagonalizeTraits` used
+    for matrix diagonalization.
+  */
+template <typename Geom_traits, typename PointRange, typename PointMap,
+          typename DiagonalizeTraits = CGAL::Default_diagonalize_traits<double,3> >
+class Planarity : public Eigen_feature<Geom_traits, PointRange, PointMap, DiagonalizeTraits>
+{
+  typedef Classification::Local_eigen_analysis<Geom_traits, PointRange,
+                                               PointMap, DiagonalizeTraits> Local_eigen_analysis;
+  typedef Eigen_feature<Geom_traits, PointRange, PointMap, DiagonalizeTraits> Base;
+public:
+  /*!
+    Constructs the feature.
+
+    \param input input range.
+    \param eigen class with precomputed eigenvectors and eigenvalues.
+  */
+  Planarity (const PointRange& input,
+             const Local_eigen_analysis& eigen)
+    : Base(input, eigen)
+  {
+    this->init(input, eigen);
+  }
+  /// \cond SKIP_IN_MANUAL
+  virtual double get_value (const Local_eigen_analysis& eigen, std::size_t i)
+  {
+    const typename Local_eigen_analysis::Eigenvalues& ev = eigen.eigenvalue(i);
+    if (ev[2] < 1e-15)
+      return 0.;
+    else
+      return ((ev[1] - ev[0]) / ev[2]);
+  }
+  virtual std::string name() { return "planarity"; }
+  /// \endcond
+ 
+};
+
+  /*!
+    \ingroup PkgClassificationFeatures
+
+    %Feature based on the eigenvalues of the covariance matrix of a
+    local neighborhood. Sphericity is defined, for the 3 eigenvalues
+    \f$\lambda_1 \ge \lambda_2 \ge \lambda_3 \ge 0\f$, as:
+
+    \f[
+    \frac{\lambda_3}{\lambda_1}
+    \f]
+    \tparam Geom_traits model of \cgal Kernel.
+    \tparam PointRange model of `ConstRange`. Its iterator type
+    is `RandomAccessIterator`.
+    \tparam PointMap model of `ReadablePropertyMap` whose key
+    type is the value type of the iterator of `PointRange` and value type
+    is `Geom_traits::Point_3`.
+    \tparam DiagonalizeTraits model of `DiagonalizeTraits` used
+    for matrix diagonalization.
+  */
+template <typename Geom_traits, typename PointRange, typename PointMap,
+          typename DiagonalizeTraits = CGAL::Default_diagonalize_traits<double,3> >
+class Sphericity : public Eigen_feature<Geom_traits, PointRange, PointMap, DiagonalizeTraits>
+{
+  typedef Classification::Local_eigen_analysis<Geom_traits, PointRange,
+                                               PointMap, DiagonalizeTraits> Local_eigen_analysis;
+  typedef Eigen_feature<Geom_traits, PointRange, PointMap, DiagonalizeTraits> Base;
+public:
+  /*!
+    Constructs the feature.
+
+    \param input input range.
+    \param eigen class with precomputed eigenvectors and eigenvalues.
+  */
+  Sphericity (const PointRange& input,
+              const Local_eigen_analysis& eigen)
+    : Base(input, eigen)
+  {
+    this->init(input, eigen);
+  }
+  /// \cond SKIP_IN_MANUAL
+  virtual double get_value (const Local_eigen_analysis& eigen, std::size_t i)
+  {
+    const typename Local_eigen_analysis::Eigenvalues& ev = eigen.eigenvalue(i);
+    if (ev[2] < 1e-15)
+      return 0.;
+    else
+      return (ev[0] / ev[2]);
+  }
+  virtual std::string name() { return "sphericity"; }
+  /// \endcond
+};
+
+  /*!
+    \ingroup PkgClassificationFeatures
+
+    %Feature based on the eigenvalues of the covariance matrix of a
+    local neighborhood. Omnivariance is defined, for the 3 eigenvalues
+    \f$\lambda_1 \ge \lambda_2 \ge \lambda_3 \ge 0\f$, as:
+
+    \f[
+    (\lambda_1 \times \lambda_2 \times \lambda_3)^{\frac{1}{3}}
+    \f]
+    \tparam Geom_traits model of \cgal Kernel.
+    \tparam PointRange model of `ConstRange`. Its iterator type
+    is `RandomAccessIterator`.
+    \tparam PointMap model of `ReadablePropertyMap` whose key
+    type is the value type of the iterator of `PointRange` and value type
+    is `Geom_traits::Point_3`.
+    \tparam DiagonalizeTraits model of `DiagonalizeTraits` used
+    for matrix diagonalization.
+  */
+template <typename Geom_traits, typename PointRange, typename PointMap,
+          typename DiagonalizeTraits = CGAL::Default_diagonalize_traits<double,3> >
+class Omnivariance : public Eigen_feature<Geom_traits, PointRange, PointMap, DiagonalizeTraits>
+{
+  typedef Classification::Local_eigen_analysis<Geom_traits, PointRange,
+                                               PointMap, DiagonalizeTraits> Local_eigen_analysis;
+  typedef Eigen_feature<Geom_traits, PointRange, PointMap, DiagonalizeTraits> Base;
+public:
+  /*!
+    Constructs the feature.
+
+    \param input input range.
+    \param eigen class with precomputed eigenvectors and eigenvalues.
+  */
+  Omnivariance (const PointRange& input,
+                const Local_eigen_analysis& eigen)
+    : Base(input, eigen)
+  {
+    this->init(input, eigen);
+  }
+  /// \cond SKIP_IN_MANUAL
+  virtual double get_value (const Local_eigen_analysis& eigen, std::size_t i)
+  {
+    const typename Local_eigen_analysis::Eigenvalues& ev = eigen.eigenvalue(i);
+    return (std::pow (std::fabs(ev[0] * ev[1] * ev[2]), 0.333333333));
+  }
+  virtual std::string name() { return "omnivariance"; }
+  /// \endcond
+};
+
+  /*!
+    \ingroup PkgClassificationFeatures
+
+    %Feature based on the eigenvalues of the covariance matrix of a
+    local neighborhood. Anisotropy is defined, for the 3 eigenvalues
+    \f$\lambda_1 \ge \lambda_2 \ge \lambda_3 \ge 0\f$, as:
+
+    \f[
+    \frac{\lambda_1 - \lambda_3}{\lambda_1}
+    \f]
+    \tparam Geom_traits model of \cgal Kernel.
+    \tparam PointRange model of `ConstRange`. Its iterator type
+    is `RandomAccessIterator`.
+    \tparam PointMap model of `ReadablePropertyMap` whose key
+    type is the value type of the iterator of `PointRange` and value type
+    is `Geom_traits::Point_3`.
+    \tparam DiagonalizeTraits model of `DiagonalizeTraits` used
+    for matrix diagonalization.
+  */
+template <typename Geom_traits, typename PointRange, typename PointMap,
+          typename DiagonalizeTraits = CGAL::Default_diagonalize_traits<double,3> >
+class Anisotropy : public Eigen_feature<Geom_traits, PointRange, PointMap, DiagonalizeTraits>
+{
+  typedef Classification::Local_eigen_analysis<Geom_traits, PointRange,
+                                               PointMap, DiagonalizeTraits> Local_eigen_analysis;
+  typedef Eigen_feature<Geom_traits, PointRange, PointMap, DiagonalizeTraits> Base;
+public:
+  /*!
+    Constructs the feature.
+
+    \param input input range.
+    \param eigen class with precomputed eigenvectors and eigenvalues.
+  */
+  Anisotropy (const PointRange& input,
+              const Local_eigen_analysis& eigen)
+    : Base(input, eigen)
+  {
+    this->init(input, eigen);
+  }
+  /// \cond SKIP_IN_MANUAL
+  virtual double get_value (const Local_eigen_analysis& eigen, std::size_t i)
+  {
+    const typename Local_eigen_analysis::Eigenvalues& ev = eigen.eigenvalue(i);
+    if (ev[2] < 1e-15)
+      return 0.;
+    else
+      return ((ev[2] - ev[0]) / ev[2]);
+  }
+  virtual std::string name() { return "anisotropy"; }
+  /// \endcond
+};
+
+  /*!
+    \ingroup PkgClassificationFeatures
+
+    %Feature based on the eigenvalues of the covariance matrix of a
+    local neighborhood. Eigentropy is defined, for the 3 eigenvalues
+    \f$\lambda_1 \ge \lambda_2 \ge \lambda_3 \ge 0\f$, as:
+
+    \f[
+    - \sum_{i=1}^3 \lambda_i \times \log{\lambda_i}
+    \f]
+    \tparam Geom_traits model of \cgal Kernel.
+    \tparam PointRange model of `ConstRange`. Its iterator type
+    is `RandomAccessIterator`.
+    \tparam PointMap model of `ReadablePropertyMap` whose key
+    type is the value type of the iterator of `PointRange` and value type
+    is `Geom_traits::Point_3`.
+    \tparam DiagonalizeTraits model of `DiagonalizeTraits` used
+    for matrix diagonalization.
+  */
+template <typename Geom_traits, typename PointRange, typename PointMap,
+          typename DiagonalizeTraits = CGAL::Default_diagonalize_traits<double,3> >
+class Eigentropy : public Eigen_feature<Geom_traits, PointRange, PointMap, DiagonalizeTraits>
+{
+  typedef Classification::Local_eigen_analysis<Geom_traits, PointRange,
+                                                    PointMap, DiagonalizeTraits> Local_eigen_analysis;
+  typedef Eigen_feature<Geom_traits, PointRange, PointMap, DiagonalizeTraits> Base;
+public:
+  /*!
+    Constructs the feature.
+
+    \param input input range.
+    \param eigen class with precomputed eigenvectors and eigenvalues.
+  */
+  Eigentropy (const PointRange& input,
+              const Local_eigen_analysis& eigen)
+    : Base(input, eigen)
+  {
+    this->init(input, eigen);
+  }
+  /// \cond SKIP_IN_MANUAL
+  virtual double get_value (const Local_eigen_analysis& eigen, std::size_t i)
+  {
+    const typename Local_eigen_analysis::Eigenvalues& ev = eigen.eigenvalue(i);
+    if (ev[0] < 1e-15
+        || ev[1] < 1e-15
+        || ev[2] < 1e-15)
+      return 0.;
+    else
+      return (- ev[0] * std::log(ev[0])
+              - ev[1] * std::log(ev[1])
+              - ev[2] * std::log(ev[2]));
+  }
+  virtual std::string name() { return "eigentropy"; }
+  /// \endcond
+};
+
+  /*!
+    \ingroup PkgClassificationFeatures
+
+    %Feature based on the eigenvalues of the covariance matrix of a
+    local neighborhood. The sum of the eigenvalues is defined, for the
+    3 eigenvalues \f$\lambda_1 \ge \lambda_2 \ge \lambda_3 \ge 0\f$,
+    as:
+
+    \f[
+    \lambda_1 + \lambda_2 + \lambda_3
+    \f]
+    \tparam Geom_traits model of \cgal Kernel.
+    \tparam PointRange model of `ConstRange`. Its iterator type
+    is `RandomAccessIterator`.
+    \tparam PointMap model of `ReadablePropertyMap` whose key
+    type is the value type of the iterator of `PointRange` and value type
+    is `Geom_traits::Point_3`.
+    \tparam DiagonalizeTraits model of `DiagonalizeTraits` used
+    for matrix diagonalization.
+  */
+template <typename Geom_traits, typename PointRange, typename PointMap,
+          typename DiagonalizeTraits = CGAL::Default_diagonalize_traits<double,3> >
+class Sum_eigenvalues : public Eigen_feature<Geom_traits, PointRange, PointMap, DiagonalizeTraits>
+{
+  typedef Classification::Local_eigen_analysis<Geom_traits, PointRange,
+                                               PointMap, DiagonalizeTraits> Local_eigen_analysis;
+  typedef Eigen_feature<Geom_traits, PointRange, PointMap, DiagonalizeTraits> Base;
+public:
+  /*!
+    Constructs the feature.
+
+    \param input input range.
+    \param eigen class with precomputed eigenvectors and eigenvalues.
+  */
+  Sum_eigenvalues (const PointRange& input,
+                   const Local_eigen_analysis& eigen)
+    : Base(input, eigen)
+  {
+    this->init(input, eigen);
+  }
+  /// \cond SKIP_IN_MANUAL
+  virtual double get_value (const Local_eigen_analysis& eigen, std::size_t i)
+  {
+    return eigen.sum_of_eigenvalues(i);
+  }
+  virtual std::string name() { return "sum_eigen"; }
+  /// \endcond
+};
+
+  /*!
+    \ingroup PkgClassificationFeatures
+
+    %Feature based on the eigenvalues of the covariance
+    matrix of a local neighborhood. Surface variation is defined, for
+    the 3 eigenvalues \f$\lambda_1 \ge \lambda_2 \ge \lambda_3 \ge
+    0\f$, as:
+
+    \f[
+    \frac{\lambda_3}{\lambda_1 + \lambda_2 + \lambda_3}
+    \f]
+    \tparam Geom_traits model of \cgal Kernel.
+    \tparam PointRange model of `ConstRange`. Its iterator type
+    is `RandomAccessIterator`.
+    \tparam PointMap model of `ReadablePropertyMap` whose key
+    type is the value type of the iterator of `PointRange` and value type
+    is `Geom_traits::Point_3`.
+    \tparam DiagonalizeTraits model of `DiagonalizeTraits` used
+    for matrix diagonalization.
+  */
+template <typename Geom_traits, typename PointRange, typename PointMap,
+          typename DiagonalizeTraits = CGAL::Default_diagonalize_traits<double,3> >
+class Surface_variation : public Eigen_feature<Geom_traits, PointRange, PointMap, DiagonalizeTraits>
+{
+  typedef Classification::Local_eigen_analysis<Geom_traits, PointRange,
+                                               PointMap, DiagonalizeTraits> Local_eigen_analysis;
+  typedef Eigen_feature<Geom_traits, PointRange, PointMap, DiagonalizeTraits> Base;
+public:
+  /*!
+    Constructs the feature.
+
+    \param input input range.
+    \param eigen class with precomputed eigenvectors and eigenvalues.
+  */
+  Surface_variation (const PointRange& input,
+                     const Local_eigen_analysis& eigen)
+    : Base(input, eigen)
+  {
+    this->init(input, eigen);
+  }
+  /// \cond SKIP_IN_MANUAL
+  virtual double get_value (const Local_eigen_analysis& eigen, std::size_t i)
+  {
+    const typename Local_eigen_analysis::Eigenvalues& ev = eigen.eigenvalue(i);
+    if (ev[0] + ev[1] + ev[2] < 1e-15)
+      return 0.;
+    else
+      return (ev[0] / (ev[0] + ev[1] + ev[2]));
+  }
+  virtual std::string name() { return "surface_variation"; }
+  /// \endcond
+};
+
+} // namespace Feature
+
+} // namespace Classification
+
+} // namespace CGAL
+
+#endif // CGAL_CLASSIFICATION_FEATURES_EIGEN_H

Classification/include/CGAL/Classification/Feature/Elevation.h

@@ -0,0 +1,170 @@
+// Copyright (c) 2012 INRIA Sophia-Antipolis (France).
+// Copyright (c) 2017 GeometryFactory Sarl (France).
+// All rights reserved.
+//
+// This file is part of CGAL (www.cgal.org).
+// You can redistribute it and/or modify it under the terms of the GNU
+// General Public License as published by the Free Software Foundation,
+// either version 3 of the License, or (at your option) any later version.
+//
+// Licensees holding a valid commercial license may use this file in
+// accordance with the commercial license agreement provided with the software.
+//
+// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
+// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
+//
+// $URL$
+// $Id$
+//
+// Author(s)     : Florent Lafarge, Simon Giraudot
+
+#ifndef CGAL_CLASSIFICATION_FEATURE_ELEVATION_H
+#define CGAL_CLASSIFICATION_FEATURE_ELEVATION_H
+
+#include <vector>
+
+#include <CGAL/Classification/Feature_base.h>
+#include <CGAL/Classification/Image.h>
+#include <CGAL/Classification/Planimetric_grid.h>
+
+namespace CGAL {
+
+namespace Classification {
+
+namespace Feature {
+
+  /*!
+    \ingroup PkgClassificationFeatures
+
+    %Feature based on local elevation. The local position of the
+    ground can be computed for urban scenes. This feature computes
+    the distance to the local estimation of the ground. It is useful
+    to discriminate the ground from horizontal roofs.
+
+    \tparam Geom_traits model of \cgal Kernel.
+    \tparam PointRange model of `ConstRange`. Its iterator type
+    is `RandomAccessIterator`.
+    \tparam PointMap model of `ReadablePropertyMap` whose key
+    type is the value type of the iterator of `PointRange` and value type
+    is `Geom_traits::Point_3`.
+
+  */
+template <typename Geom_traits, typename PointRange, typename PointMap>
+class Elevation : public Feature_base
+{
+  typedef typename Geom_traits::Iso_cuboid_3 Iso_cuboid_3;
+
+  typedef Image<float> Image_float;
+  typedef Planimetric_grid<Geom_traits, PointRange, PointMap> Grid;
+   
+  std::vector<double> elevation_feature;
+  
+public:
+  /*!
+    \brief Constructs the feature.
+
+    \param input input range.
+    \param point_map property map to access the input points.
+    \param grid precomputed `Planimetric_grid`.
+    \param grid_resolution resolution of the planimetric grid.
+    \param radius_dtm radius for digital terrain modeling (should be
+    larger than the width and length of the largest building).
+  */
+  Elevation (const PointRange& input,
+             PointMap point_map,
+             const Grid& grid,
+             const double grid_resolution,
+             double radius_dtm = -1.)
+  {
+    this->set_weight(1.);
+    if (radius_dtm < 0.)
+      radius_dtm = 100. * grid_resolution;
+
+    //DEM
+    Image_float dem(grid.width(),grid.height());
+
+    for (std::size_t j = 0; j < grid.height(); ++ j)
+      for (std::size_t i = 0; i < grid.width(); ++ i)
+        {
+          if (grid.indices(i,j).empty())
+            continue;
+          double mean = 0.;
+          for (std::size_t k = 0; k < grid.indices(i,j).size(); ++ k)
+            mean += get(point_map, *(input.begin()+grid.indices(i,j)[k])).z();
+          mean /= grid.indices(i,j).size();
+          dem(i,j) = mean;
+        }
+
+    std::size_t square = (std::size_t)(0.5 * radius_dtm / grid_resolution) + 1;
+    
+    Image_float dtm_x(grid.width(),grid.height());
+    
+    for (std::size_t j = 0; j < grid.height(); ++ j)
+      {
+        for (std::size_t i = 0; i < grid.width(); ++ i)
+          {
+            std::size_t squareXmin = (i < square ? 0 : i - square);
+            std::size_t squareXmax = (std::min)(grid.width() - 1, i + square);
+
+            std::vector<double> z;
+            z.reserve(squareXmax - squareXmin +1 );
+            for(std::size_t k = squareXmin; k <= squareXmax; k++)
+              if (dem(k,j) != 0.)
+                z.push_back (dem(k,j));
+            if (z.empty())
+              continue;
+            std::nth_element (z.begin(), z.begin() + (z.size() / 10), z.end());
+            dtm_x(i,j) = z[z.size() / 10];
+          }
+      }
+    dem.free();
+
+    Image_float dtm(grid.width(),grid.height());
+    
+    for (std::size_t i = 0; i < grid.width(); ++ i)
+      {
+        for (std::size_t j = 0; j < grid.height(); ++ j)
+          {
+            std::size_t squareYmin = (j < square ? 0 : j - square);
+            std::size_t squareYmax = (std::min)(grid.height() - 1, j + square);
+            std::vector<double> z;
+            z.reserve(squareYmax - squareYmin +1 );
+            for(std::size_t l = squareYmin; l <= squareYmax; l++)
+              if (dtm_x(i,l) != 0.)
+                z.push_back (dtm_x(i,l));
+            if (z.empty())
+              continue;
+            std::nth_element (z.begin(), z.begin() + (z.size() / 10), z.end());
+            dtm(i,j) = z[z.size() / 10];
+          }
+      }
+    dtm_x.free();
+
+    elevation_feature.reserve(input.size());
+    for (std::size_t i = 0; i < input.size(); i++){
+      std::size_t I = grid.x(i);
+      std::size_t J = grid.y(i);
+      elevation_feature.push_back ((double)(get(point_map, *(input.begin()+i)).z()-dtm(I,J)));
+    }
+
+    this->compute_mean_max (elevation_feature, this->mean, this->max);
+  }
+
+  /// \cond SKIP_IN_MANUAL
+  virtual double value (std::size_t pt_index)
+  {
+    return elevation_feature[pt_index];
+  }
+  
+  virtual std::string name() { return "elevation"; }
+  /// \endcond
+};
+
+} // namespace Feature
+
+} // namespace Classification
+
+
+} // namespace CGAL
+
+#endif // CGAL_CLASSIFICATION_FEATURE_ELEVATION_H

Classification/include/CGAL/Classification/Feature/Hsv.h

@@ -0,0 +1,148 @@
+// Copyright (c) 2017 GeometryFactory Sarl (France).
+// All rights reserved.
+//
+// This file is part of CGAL (www.cgal.org).
+// You can redistribute it and/or modify it under the terms of the GNU
+// General Public License as published by the Free Software Foundation,
+// either version 3 of the License, or (at your option) any later version.
+//
+// Licensees holding a valid commercial license may use this file in
+// accordance with the commercial license agreement provided with the software.
+//
+// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
+// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
+//
+// $URL$
+// $Id$
+//
+// Author(s)     : Simon Giraudot
+
+#ifndef CGAL_CLASSIFICATION_FEATURE_HSV_H
+#define CGAL_CLASSIFICATION_FEATURE_HSV_H
+
+#include <vector>
+
+#include <CGAL/Classification/Color.h>
+
+namespace CGAL {
+
+namespace Classification {
+
+namespace Feature {
+  
+  /*!
+    \ingroup PkgClassificationFeatures
+
+    %Feature based on HSV colorimetric information. If the input
+    point cloud has colorimetric information, it can be used for
+    classification purposes. This feature is based on a Gaussian
+    probabilistic model on one of the three HSV channels (hue,
+    saturation or value). It computes the probability of the color of
+    the input point to match this specific color channel defined by a
+    mean and a standard deviation.
+
+    The HSV channels are defined this way:
+
+    - Hue ranges from 0 to 360 and measures the general "tint" of the
+      color (green, blue, pink, etc.)
+
+    - Saturation ranges from 0 to 100 and measures the "strength" of the
+      color (0 is gray and 100 is the fully saturated color)
+
+    - Value ranges from 0 to 100 and measures the "brightness" of the
+      color (0 is black and 100 is the fully bright color)
+
+    For example, such an feature using the channel 0 (hue) with a
+    mean of 90 (which corresponds to a green hue) can help to identify
+    trees.
+    
+    \note The user only needs to provide a map to standard (and more common)
+    RGB colors, the conversion to HSV is done internally.
+
+    \tparam Geom_traits model of \cgal Kernel.
+    \tparam PointRange model of `ConstRange`. Its iterator type
+    is `RandomAccessIterator`.
+    \tparam ColorMap model of `ReadablePropertyMap` whose key
+    type is the value type of the iterator of `PointRange` and value type
+    is `CGAL::Classification::RGB_Color`.
+  */
+template <typename Geom_traits, typename PointRange, typename ColorMap>
+class Hsv : public Feature_base
+{
+  typedef typename Classification::RGB_Color RGB_Color;
+  typedef typename Classification::HSV_Color HSV_Color;
+
+#ifdef CGAL_CLASSIFICATION_PRECOMPUTE_FEATURES
+  std::vector<double> color_feature;
+#else
+  const PointRange& input;
+  ColorMap color_map;
+  std::size_t channel;
+  double mean;
+  double sd;
+#endif  
+  std::string m_name;
+  
+public:
+  
+  /*!
+
+    \brief Constructs an feature based on the given color channel,
+    mean and standard deviation.
+
+    \param input input range.
+    \param color_map property map to access the colors of the input points.
+    \param channel chosen HSV channel (0 for hue, 1 for saturation, 2 for value).
+    \param mean mean value of the specified channel.
+    \param sd standard deviation of the specified channel.
+  */
+  Hsv (const PointRange& input,
+       ColorMap color_map,
+       std::size_t channel,
+       double mean, double sd)
+#ifndef CGAL_CLASSIFICATION_PRECOMPUTE_FEATURES
+    : input(input), color_map(color_map), channel(channel), mean(mean), sd(sd)
+#endif
+  {
+    this->set_weight(1.);
+
+#ifndef CGAL_CLASSIFICATION_PRECOMPUTE_FEATURES
+    std::vector<double> color_feature;
+#endif
+    for(std::size_t i = 0; i < input.size();i++)
+      {
+        HSV_Color c = Classification::rgb_to_hsv (get(color_map, *(input.begin()+i)));
+        color_feature.push_back (std::exp (-(c[channel] - mean) * (c[channel] - mean) / (2. * sd * sd)));
+      }
+    this->compute_mean_max (color_feature, this->mean, this->max);
+
+    std::ostringstream oss;
+    if (channel == 0) oss << "hue";
+    else if (channel == 1) oss << "saturation";
+    else if (channel == 2) oss << "value";
+    oss << "_" << mean;
+    m_name = oss.str();
+  }
+
+  /// \cond SKIP_IN_MANUAL
+  virtual double value (std::size_t pt_index)
+  {
+#ifdef CGAL_CLASSIFICATION_PRECOMPUTE_FEATURES
+    return color_feature[pt_index];
+#else
+    HSV_Color c = Classification::rgb_to_hsv (get(color_map, *(input.begin()+pt_index)));
+    return std::exp (-(c[channel] - mean) * (c[channel] - mean) / (2. * sd * sd));
+#endif
+  }
+
+  virtual std::string name() { return m_name; }
+  /// \endcond
+};
+
+} // namespace Feature
+
+} // namespace Classification
+
+} // namespace CGAL
+
+#endif // CGAL_CLASSIFICATION_FEATURE_HSV_H

Classification/include/CGAL/Classification/Feature/Vertical_dispersion.h

@@ -0,0 +1,163 @@
+// Copyright (c) 2017 GeometryFactory Sarl (France).
+// All rights reserved.
+//
+// This file is part of CGAL (www.cgal.org).
+// You can redistribute it and/or modify it under the terms of the GNU
+// General Public License as published by the Free Software Foundation,
+// either version 3 of the License, or (at your option) any later version.
+//
+// Licensees holding a valid commercial license may use this file in
+// accordance with the commercial license agreement provided with the software.
+//
+// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
+// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
+//
+// $URL$
+// $Id$
+//
+// Author(s)     : Simon Giraudot
+
+#ifndef CGAL_CLASSIFICATION_FEATURE_VERTICAL_DISPERSION_H
+#define CGAL_CLASSIFICATION_FEATURE_VERTICAL_DISPERSION_H
+
+#include <vector>
+
+#include <CGAL/Classification/Image.h>
+#include <CGAL/Classification/Planimetric_grid.h>
+
+namespace CGAL {
+
+namespace Classification {
+
+namespace Feature {
+  
+  /*!
+    \ingroup PkgClassificationFeatures
+
+    %Feature based on local vertical dispersion of points. Urban
+    scenes can often be decomposed as a set of 2D regions with
+    different heights. While these heights are usually piecewise
+    constant or piecewise linear, on some specific parts of the scene
+    such as vegetation, they can become extremely unstable. This
+    feature quantifies the vertical dispersion of the points on a
+    local Z-cylinder around the points.
+
+    \tparam Geom_traits model of \cgal Kernel.
+    \tparam PointRange model of `ConstRange`. Its iterator type
+    is `RandomAccessIterator`.
+    \tparam PointMap model of `ReadablePropertyMap` whose key
+    type is the value type of the iterator of `PointRange` and value type
+    is `Geom_traits::Point_3`.
+  */
+template <typename Geom_traits, typename PointRange, typename PointMap>
+class Vertical_dispersion : public Feature_base
+{
+  typedef Classification::Image<float> Image_float;
+  typedef Classification::Planimetric_grid<Geom_traits, PointRange, PointMap> Grid;
+  std::vector<double> vertical_dispersion;
+  
+public:
+  /*!
+    \brief Constructs the feature.
+
+    \param input input range.
+    \param point_map property map to access the input points.
+    \param grid precomputed `Planimetric_grid`.
+    \param grid_resolution resolution of the planimetric grid.
+    \param radius_neighbors radius of local neighborhoods.
+  */
+  Vertical_dispersion (const PointRange& input,
+                       PointMap point_map,
+                       const Grid& grid,
+                       const double grid_resolution,
+                       double radius_neighbors = -1.)
+  {
+    this->set_weight(1.);
+    if (radius_neighbors < 0.)
+      radius_neighbors = 5. * grid_resolution;
+    
+    Image_float Dispersion(grid.width(), grid.height());
+    for (std::size_t j = 0; j < grid.height(); j++)	
+      for (std::size_t i = 0; i < grid.width(); i++)
+        Dispersion(i,j)=0;
+    
+    std::size_t square = (std::size_t)(0.5 * radius_neighbors / grid_resolution) + 1;
+    typename Geom_traits::Vector_3 verti (0., 0., 1.);
+    
+    for (std::size_t j = 0; j < grid.height(); j++){	
+      for (std::size_t i = 0; i < grid.width(); i++){
+						
+        if(!(grid.mask(i,j)))
+          continue;
+        std::vector<double> hori;
+            
+        std::size_t squareXmin = (i < square ? 0 : i - square);
+        std::size_t squareXmax = (std::min) (grid.width()-1, i + square);
+        std::size_t squareYmin = (j < square ? 0 : j - square);
+        std::size_t squareYmax = (std::min) (grid.height()-1, j + square);
+
+        for(std::size_t k = squareXmin; k <= squareXmax; k++)
+          for(std::size_t l = squareYmin; l <= squareYmax; l++)
+            if(CGAL::sqrt(pow((double)k-i,2)+pow((double)l-j,2))
+               <=(double)0.5*radius_neighbors/grid_resolution
+               && (grid.indices(k,l).size()>0))
+              for(std::size_t t=0; t<grid.indices(k,l).size();t++)
+                {
+                  std::size_t ip = grid.indices(k,l)[t];
+                  hori.push_back (get(point_map, *(input.begin()+ip)).z());
+                }
+        if (hori.empty())
+          continue;
+              
+        std::sort (hori.begin(), hori.end());
+
+        std::size_t nb_layers = 1;
+
+        std::vector<bool> occupy (1 + (std::size_t)((hori.back() - hori.front()) / grid_resolution), false);
+              
+        std::size_t last_index = 0;
+        for (std::size_t k = 0; k < hori.size(); ++ k)
+          {
+            std::size_t index = (std::size_t)((hori[k] - hori.front()) / grid_resolution);
+            occupy[index] = true;
+            if (index > last_index + 1)
+              ++ nb_layers;
+            last_index = index;
+          }
+
+        std::size_t nb_occ = 0;
+        for (std::size_t k = 0; k < occupy.size(); ++ k)
+          if (occupy[k])
+            ++ nb_occ;
+					
+        Dispersion(i,j)= 1.f - (nb_occ / (float)(occupy.size()));
+			
+      }
+		
+    }
+    for (std::size_t i = 0; i < input.size(); i++)
+      {
+        std::size_t I= grid.x(i);
+        std::size_t J= grid.y(i);
+        vertical_dispersion.push_back((double)Dispersion(I,J));
+      }
+
+    this->compute_mean_max (vertical_dispersion, this->mean, this->max);
+  }
+  /// \cond SKIP_IN_MANUAL
+  virtual double value (std::size_t pt_index)
+  {
+    return vertical_dispersion[pt_index];
+  }
+
+  virtual std::string name() { return "vertical_dispersion"; }
+  /// \endcond
+};
+
+}
+  
+}
+
+}
+
+#endif // CGAL_CLASSIFICATION_FEATURE_VERTICAL_DISPERSION_H

Classification/include/CGAL/Classification/Feature/Verticality.h

@@ -0,0 +1,128 @@
+// Copyright (c) 2017 GeometryFactory Sarl (France).
+// All rights reserved.
+//
+// This file is part of CGAL (www.cgal.org).
+// You can redistribute it and/or modify it under the terms of the GNU
+// General Public License as published by the Free Software Foundation,
+// either version 3 of the License, or (at your option) any later version.
+//
+// Licensees holding a valid commercial license may use this file in
+// accordance with the commercial license agreement provided with the software.
+//
+// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
+// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
+//
+// $URL$
+// $Id$
+//
+// Author(s)     : Simon Giraudot
+
+#ifndef CGAL_CLASSIFICATION_FEATURE_VERTICALITY_H
+#define CGAL_CLASSIFICATION_FEATURE_VERTICALITY_H
+
+#include <vector>
+
+#include <CGAL/Classification/Local_eigen_analysis.h>
+
+namespace CGAL {
+
+namespace Classification {
+
+namespace Feature {
+
+  /*!
+    \ingroup PkgClassificationFeatures
+
+    %Feature based on local verticality. The orientation of the
+    local tangent plane of the considered point can be useful to
+    discriminate facades from the ground. This feature can use
+    normal vectors if available or eigen analysis that estimates
+    tangent planes from local neighborhoods.
+
+    \tparam Geom_traits model of \cgal Kernel.
+    \tparam PointRange model of `ConstRange`. Its iterator type
+    is `RandomAccessIterator`.
+    \tparam PointMap model of `ReadablePropertyMap` whose key
+    type is the value type of the iterator of `PointRange` and value type
+    is `Geom_traits::Point_3`.
+    \tparam DiagonalizeTraits model of `DiagonalizeTraits` used
+    for matrix diagonalization.
+  */
+template <typename Geom_traits, typename PointRange, typename PointMap,
+          typename DiagonalizeTraits = CGAL::Default_diagonalize_traits<double,3> >
+class Verticality : public Feature_base
+{
+  typedef Classification::Local_eigen_analysis<Geom_traits, PointRange,
+                                               PointMap, DiagonalizeTraits> Local_eigen_analysis;
+
+  std::vector<double> verticality_feature;
+  
+public:
+  /*!
+    \brief Constructs the feature using local eigen analysis.
+
+    \param input input range.
+    \param eigen class with precomputed eigenvectors and eigenvalues.
+  */
+  Verticality (const PointRange& input,
+               const Local_eigen_analysis& eigen)
+  {
+    this->set_weight(1.);
+    typename Geom_traits::Vector_3 vertical (0., 0., 1.);
+
+    for (std::size_t i = 0; i < input.size(); i++)
+      {
+        typename Geom_traits::Vector_3 normal = eigen.normal_vector(i);
+        normal = normal / CGAL::sqrt (normal * normal);
+        verticality_feature.push_back (1. - CGAL::abs(normal * vertical));
+      }
+    
+    this->compute_mean_max (verticality_feature, this->mean, this->max);
+    //    max *= 2;
+  }
+
+  /*!
+    \brief Constructs the feature using provided normals of points.
+
+    \tparam VectorMap model of `ReadablePropertyMap` whose key
+    type is the value type of the iterator of `PointRange` and value type
+    is `Geom_traits::Vector_3`.
+    \param input input range.
+    \param normal_map property map to access the normal vectors of the input points.
+  */
+  template <typename VectorMap>
+  Verticality (const PointRange& input,
+               VectorMap normal_map)
+  {
+    this->set_weight(1.);
+    typename Geom_traits::Vector_3 vertical (0., 0., 1.);
+
+    for (std::size_t i = 0; i < input.size(); i++)
+      {
+        typename Geom_traits::Vector_3 normal = get(normal_map, *(input.begin()+i));
+        normal = normal / CGAL::sqrt (normal * normal);
+        verticality_feature.push_back (1. - std::fabs(normal * vertical));
+      }
+    
+    this->compute_mean_max (verticality_feature, this->mean, this->max);
+    //    max *= 2;
+  }
+
+
+  /// \cond SKIP_IN_MANUAL
+  virtual double value (std::size_t pt_index)
+  {
+    return verticality_feature[pt_index];
+  }
+
+  virtual std::string name() { return "verticality"; }
+  /// \endcond
+};
+
+} // namespace Feature
+
+} // namespace Classification
+
+} // namespace CGAL
+
+#endif // CGAL_CLASSIFICATION_FEATURE_VERTICALITY_H

Classification/include/CGAL/Classification/Feature_base.h

@@ -0,0 +1,131 @@
+// Copyright (c) 2017 GeometryFactory Sarl (France).
+// All rights reserved.
+//
+// This file is part of CGAL (www.cgal.org).
+// You can redistribute it and/or modify it under the terms of the GNU
+// General Public License as published by the Free Software Foundation,
+// either version 3 of the License, or (at your option) any later version.
+//
+// Licensees holding a valid commercial license may use this file in
+// accordance with the commercial license agreement provided with the software.
+//
+// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
+// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
+//
+// $URL$
+// $Id$
+//
+// Author(s)     : Simon Giraudot
+
+#ifndef CGAL_CLASSIFICATION_FEATURE_BASE_H
+#define CGAL_CLASSIFICATION_FEATURE_BASE_H
+
+#include <boost/shared_ptr.hpp>
+
+#include <vector>
+
+namespace CGAL {
+
+namespace Classification {
+  
+/*!
+  \ingroup PkgClassification
+
+  \brief Abstract class describing a classification feature that
+  associates a scalar value to each item of the classification input.
+
+  User-defined features must inherit this class. Feature objects are
+  generated by `CGAL::Classifier` objects by the method
+  `CGAL::Classifier::add_feature()`. For this method to work, the
+  first argument of the feature's constructors must be of type
+  `ItemRange&` (where `ItemRange` is the first template argument of
+  `CGAL::Classifier`).
+*/
+
+
+class Feature_base
+{
+  double m_weight; 
+public:
+
+  /// \cond SKIP_IN_MANUAL
+  double mean;
+  double max;
+
+  virtual ~Feature_base() { }
+  /// \endcond
+
+  /*!
+    \brief Returns the weight of the feature.
+  */
+  double weight() const { return m_weight; }
+
+  /*!
+    \brief Sets the weight of the feature (`weight` must be positive).
+  */
+  void set_weight (double weight) { m_weight = weight; }
+
+  /*!
+    \brief Returns `abstract_feature` and should be overloaded
+    by an inherited class with a specific name.
+  */
+  virtual std::string name() { return "abstract_feature"; }
+
+  /*!
+    \brief Returns the value taken by the feature for at the item at
+    position `index`. This method must be implemented by inherited
+    classes.
+  */
+  virtual double value (std::size_t index) = 0;
+
+  /// \cond SKIP_IN_MANUAL
+  virtual double normalized (std::size_t index)
+  {
+    return (std::max) (0., (std::min) (1., value(index) / m_weight));
+  }
+  virtual double favored (std::size_t index) { return (1. - normalized (index)); }
+  virtual double penalized (std::size_t index) { return normalized (index); }
+  //  virtual double ignored (std::size_t index) { return (std::min) (favored(index), penalized(index)); }
+  virtual double ignored (std::size_t) { return 0.5; }
+
+  void compute_mean_max (std::vector<double>& vect, double& mean, double& max)
+  {
+    mean = 0.;
+    max = -std::numeric_limits<double>::max();
+    double min = std::numeric_limits<double>::max();
+    
+    for (std::size_t i = 0; i < vect.size(); ++ i)
+      {
+        mean += vect[i];
+        if (vect[i] > max)
+          max = vect[i];
+        if (vect[i] < min)
+          min = vect[i];
+      }
+    //    std::cerr << name() << " Min/max = " << min << " / " << max << std::endl;
+    mean /= vect.size();
+
+  }
+  /// \endcond
+};
+
+
+#ifdef DOXYGEN_RUNNING
+/*!
+  \ingroup PkgClassification
+
+  \brief %Handle to an `Feature_base`.
+
+  \cgalModels Handle
+*/
+class Feature_handle { };
+#else
+typedef boost::shared_ptr<Feature_base> Feature_handle;
+#endif
+
+
+} // namespace Classification
+
+} // namespace CGAL
+
+#endif // CGAL_CLASSIFICATION_FEATURE_BASE_H

Classification/include/CGAL/Classification/Image.h

@@ -0,0 +1,118 @@
+// Copyright (c) 2017 GeometryFactory Sarl (France).
+// All rights reserved.
+//
+// This file is part of CGAL (www.cgal.org).
+// You can redistribute it and/or modify it under the terms of the GNU
+// General Public License as published by the Free Software Foundation,
+// either version 3 of the License, or (at your option) any later version.
+//
+// Licensees holding a valid commercial license may use this file in
+// accordance with the commercial license agreement provided with the software.
+//
+// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
+// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
+//
+// $URL$
+// $Id$
+//
+// Author(s)     : Simon Giraudot
+
+#ifndef CGAL_CLASSIFICATION_IMAGE_H
+#define CGAL_CLASSIFICATION_IMAGE_H
+
+namespace CGAL {
+namespace Classification {
+
+  /// \cond SKIP_IN_MANUAL
+  
+template <typename Type>
+class Image
+{
+  std::size_t m_width;
+  std::size_t m_height;
+  Type* m_raw;
+
+
+public:
+
+  Image () : m_width(0), m_height(0), m_raw (NULL)
+  {
+  }
+  
+  Image (std::size_t width, std::size_t height)
+    : m_width (width),
+      m_height (height)
+  {
+    if (m_width * m_height > 0)
+      m_raw = new Type[width * height]();
+    else
+      m_raw = NULL;
+  }
+  
+  ~Image ()
+  {
+    free();
+  }
+
+  void free()
+  {
+    if (m_raw != NULL)
+      delete[] m_raw;
+    m_raw = NULL;
+  }
+
+  Image (const Image& other)
+    : m_width (other.width()),
+      m_height (other.height())
+
+  {
+    if (m_width * m_height > 0)
+      {
+        m_raw = new Type[m_width * m_height];
+        std::copy (other.m_raw, other.m_raw + (m_width * m_height), this->m_raw);
+      }
+    else
+      m_raw = NULL;
+  }
+  Image& operator= (const Image& other)
+  {
+    if (m_raw != NULL)
+      delete[] m_raw;
+
+    m_raw = NULL;
+    m_width = other.width();
+    m_height = other.height();
+    if (m_width * m_height > 0)
+      {
+        m_raw = new Type[m_width * m_height];
+        std::copy (other.m_raw, other.m_raw + (m_width * m_height), this->m_raw);
+      }
+    
+    return *this;
+  }
+  
+  std::size_t width() const { return m_width; }
+  std::size_t height() const { return m_height; }
+
+  Type& operator() (const std::size_t& x, const std::size_t& y)
+  {
+    //    return m_raw[y * m_width + x];
+    return m_raw[x * m_height + y];
+  }
+  const Type& operator() (const std::size_t& x, const std::size_t& y) const
+  {
+    //    return m_raw[y * m_width + x];
+    return m_raw[x * m_height + y];
+  }
+  
+
+};
+
+  /// \endcond
+
+} // namespace Classification
+} // namespace CGAL
+
+
+
+#endif // CGAL_CLASSIFICATION_IMAGE_H

Classification/include/CGAL/Classification/Label.h

@@ -0,0 +1,112 @@
+// Copyright (c) 2017 GeometryFactory Sarl (France).
+// All rights reserved.
+//
+// This file is part of CGAL (www.cgal.org).
+// You can redistribute it and/or modify it under the terms of the GNU
+// General Public License as published by the Free Software Foundation,
+// either version 3 of the License, or (at your option) any later version.
+//
+// Licensees holding a valid commercial license may use this file in
+// accordance with the commercial license agreement provided with the software.
+//
+// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
+// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
+//
+// $URL$
+// $Id$
+//
+// Author(s)     : Simon Giraudot
+
+#ifndef CGAL_CLASSIFICATION_LABEL_H
+#define CGAL_CLASSIFICATION_LABEL_H
+
+#include <CGAL/Classification/Feature/Effect.h>
+
+#include <boost/shared_ptr.hpp>
+
+namespace CGAL {
+
+namespace Classification {
+
+/*!
+\ingroup PkgClassification
+
+\brief %Classification label (for example: vegetation, ground, etc.)
+defined as a set of relationships with classification features.
+
+*/
+class Label
+{
+public:
+  
+private:
+  /// \cond SKIP_IN_MANUAL
+  std::string m_name;
+  std::map<Feature_handle, Feature::Effect> m_feature_effects;
+  /// \endcond
+
+public:
+
+  /*! 
+    \param name name of the classification label (e.g. vegetation).
+  */ 
+  Label (std::string name) : m_name (name) { }
+
+  /*! 
+    \brief Sets the effect of feature `att` on the classification label.
+  */ 
+  void set_feature_effect (Feature_handle att, Feature::Effect effect)
+  {
+    m_feature_effects[att] = effect;
+  }
+
+  /*!
+    \brief Returns the effect of feature `att` on the classification label.
+   */
+  Feature::Effect feature_effect (Feature_handle att) 
+  {
+    std::map<Feature_handle, Feature::Effect>::iterator
+      search = m_feature_effects.find (att);
+    return (search == m_feature_effects.end () ? Feature::NEUTRAL : search->second);
+  }
+
+  const std::string& name() const { return m_name; }
+  
+  /// \cond SKIP_IN_MANUAL
+  void info()
+  {
+    std::cerr << "Feature " << m_name << ": ";
+    for (std::map<Feature_handle, Feature::Effect>::iterator it = m_feature_effects.begin();
+         it != m_feature_effects.end(); ++ it)
+      {
+        if (it->second == Feature::NEUTRAL)
+          continue;
+        
+        std::cerr << it->first;
+        if (it->second == Feature::FAVORING) std::cerr << " (favored), ";
+        else if (it->second == Feature::PENALIZING) std::cerr << " (penalized), ";
+      }
+    std::cerr << std::endl;
+  }
+  /// \endcond
+
+};
+
+#ifdef DOXYGEN_RUNNING
+/*!
+  \ingroup PkgClassification
+
+  \brief %Handle to a classification `Label`.
+
+  \cgalModels Handle
+*/
+class Label_handle { };
+#else
+typedef boost::shared_ptr<Label> Label_handle;
+#endif
+
+} // namespace Classification
+
+} // namespace CGAL
+
+#endif // CGAL_CLASSIFICATION_LABEL_H

Classification/include/CGAL/Classification/Local_eigen_analysis.h

@@ -0,0 +1,282 @@
+// Copyright (c) 2017 GeometryFactory Sarl (France).
+// All rights reserved.
+//
+// This file is part of CGAL (www.cgal.org).
+// You can redistribute it and/or modify it under the terms of the GNU
+// General Public License as published by the Free Software Foundation,
+// either version 3 of the License, or (at your option) any later version.
+//
+// Licensees holding a valid commercial license may use this file in
+// accordance with the commercial license agreement provided with the software.
+//
+// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
+// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
+//
+// $URL$
+// $Id$
+//
+// Author(s)     : Simon Giraudot
+
+#ifndef CGAL_CLASSIFICATION_LOCAL_EIGEN_ANALYSIS_H
+#define CGAL_CLASSIFICATION_LOCAL_EIGEN_ANALYSIS_H
+
+#include <vector>
+
+#include <CGAL/Search_traits_3.h>
+#include <CGAL/Fuzzy_sphere.h>
+#include <CGAL/Orthogonal_k_neighbor_search.h>
+#include <CGAL/Default_diagonalize_traits.h>
+#include <CGAL/centroid.h>
+
+#ifdef CGAL_LINKED_WITH_TBB
+#include <tbb/parallel_for.h>
+#include <tbb/blocked_range.h>
+#include <tbb/scalable_allocator.h>
+#include <tbb/mutex.h>
+#endif // CGAL_LINKED_WITH_TBB
+
+namespace CGAL {
+
+namespace Classification {
+
+  /*!
+    \ingroup PkgClassificationDataStructures
+
+    \brief Class that precomputes and stored the eigenvectors and
+    eigenvalues of the covariance matrices of all points of a point
+    set using a local neighborhood.
+
+    \tparam Geom_traits model of \cgal Kernel.
+    \tparam PointRange model of `ConstRange`. Its iterator type is
+    `RandomAccessIterator`.
+    \tparam PointMap model of `ReadablePropertyMap` whose key
+    type is the value type of the iterator of `PointRange` and value type
+    is `Geom_traits::Point_3`.
+    \tparam DiagonalizeTraits model of `DiagonalizeTraits` used
+    for matrix diagonalization.
+  */
+template <typename Geom_traits, typename PointRange, typename PointMap,
+          typename DiagonalizeTraits = CGAL::Default_diagonalize_traits<double,3> >
+class Local_eigen_analysis
+{
+public:
+  typedef typename Geom_traits::Point_3 Point;
+  typedef typename Geom_traits::Vector_3 Vector;
+  typedef typename Geom_traits::Plane_3 Plane;
+  
+  typedef CGAL::cpp11::array<double, 3> Eigenvalues; ///< Eigenvalues (sorted in ascending order)
+
+
+private:
+
+#ifdef CGAL_LINKED_WITH_TBB
+  template <typename NeighborQuery>
+  class Compute_eigen_values
+  {
+    Local_eigen_analysis& m_eigen;
+    const PointRange& m_input;
+    PointMap m_point_map;
+    const NeighborQuery& m_neighbor_query;
+    double& m_mean_range;
+    tbb::mutex& m_mutex;
+    
+  public:
+    
+    Compute_eigen_values (Local_eigen_analysis& eigen,
+                          const PointRange& input,
+                          PointMap point_map,
+                          const NeighborQuery& neighbor_query,
+                          double& mean_range,
+                          tbb::mutex& mutex)
+      : m_eigen (eigen), m_input (input), m_point_map (point_map),
+        m_neighbor_query (neighbor_query), m_mean_range (mean_range), m_mutex (mutex)
+    { }
+    
+    void operator()(const tbb::blocked_range<std::size_t>& r) const
+    {
+      for (std::size_t i = r.begin(); i != r.end(); ++ i)
+        {
+          std::vector<std::size_t> neighbors;
+          m_neighbor_query (get(m_point_map, *(m_input.begin()+i)), std::back_inserter (neighbors));
+
+          std::vector<Point> neighbor_points;
+          for (std::size_t j = 0; j < neighbors.size(); ++ j)
+            neighbor_points.push_back (get(m_point_map, *(m_input.begin()+neighbors[j])));
+
+          m_mutex.lock();
+          m_mean_range += CGAL::sqrt (CGAL::squared_distance (get(m_point_map, *(m_input.begin() + i)),
+                                                              get(m_point_map, *(m_input.begin() + neighbors.back()))));
+          m_mutex.unlock();
+          
+          m_eigen.compute (i, get(m_point_map, *(m_input.begin()+i)), neighbor_points);
+        }
+    }
+
+  };
+#endif
+  
+  std::vector<Eigenvalues> m_eigenvalues;
+  std::vector<double> m_sum_eigenvalues;
+  std::vector<Point> m_centroids;
+  std::vector<Vector> m_smallest_eigenvectors;
+  std::vector<Vector> m_middle_eigenvectors;
+  std::vector<Vector> m_largest_eigenvectors;
+  double m_mean_range;
+
+  
+public:
+
+  /// \cond SKIP_IN_MANUAL
+  Local_eigen_analysis () { }
+  /// \endcond
+
+  /*! 
+
+    \brief Computes the local eigen analysis of an input range based
+    on a local neighborhood.
+
+    \tparam NeighborQuery model of `NeighborQuery`
+    \tparam ConcurrencyTag enables sequential versus parallel
+    algorithm. Possible values are `Parallel_tag` (default value is %CGAL
+    is linked with TBB) or `Sequential_tag` (default value otherwise).
+    \param input input range.
+    \param point_map property map to access the input points
+    \param neighbor_query object used to access neighborhoods of points.
+  */
+  template <typename NeighborQuery,
+#if defined(DOXYGEN_RUNNING)
+            typename ConcurrencyTag>
+#elif defined(CGAL_LINKED_WITH_TBB)
+            typename ConcurrencyTag = CGAL::Parallel_tag>
+#else
+            typename ConcurrencyTag = CGAL::Sequential_tag>
+#endif
+  Local_eigen_analysis (const PointRange& input,
+                        PointMap point_map,
+                        const NeighborQuery& neighbor_query)
+  {
+    m_eigenvalues.resize (input.size());
+    m_sum_eigenvalues.resize (input.size());
+    m_centroids.resize (input.size());
+    m_smallest_eigenvectors.resize (input.size());
+    m_middle_eigenvectors.resize (input.size());
+    m_largest_eigenvectors.resize (input.size());
+    
+    m_mean_range = 0.;
+      
+#ifndef CGAL_LINKED_WITH_TBB
+    CGAL_static_assertion_msg (!(boost::is_convertible<ConcurrencyTag, Parallel_tag>::value),
+                               "Parallel_tag is enabled but TBB is unavailable.");
+#else
+    if (boost::is_convertible<ConcurrencyTag,Parallel_tag>::value)
+      {
+        tbb::mutex mutex;
+        Compute_eigen_values<NeighborQuery> f(*this, input, point_map, neighbor_query, m_mean_range, mutex);
+        tbb::parallel_for(tbb::blocked_range<size_t>(0, input.size ()), f);
+      }
+    else
+#endif
+      {
+        for (std::size_t i = 0; i < input.size(); i++)
+          {
+            std::vector<std::size_t> neighbors;
+            neighbor_query (get(point_map, *(input.begin()+i)), std::back_inserter (neighbors));
+
+            std::vector<Point> neighbor_points;
+            for (std::size_t j = 0; j < neighbors.size(); ++ j)
+              neighbor_points.push_back (get(point_map, *(input.begin()+neighbors[j])));
+
+            m_mean_range += CGAL::sqrt (CGAL::squared_distance (get(point_map, *(input.begin() + i)),
+                                                                get(point_map, *(input.begin() + neighbors.back()))));
+        
+            compute (i, get(point_map, *(input.begin()+i)), neighbor_points);
+          }
+      }
+    m_mean_range /= input.size();
+  }
+
+  /*!
+    \brief Returns the estimated unoriented normal vector of the point at position `index`.
+  */
+  const Vector& normal_vector (std::size_t index) const { return m_smallest_eigenvectors[index]; }
+
+  /*!
+    \brief Returns the estimated local tangent plane of the point at position `index`.
+  */
+  Plane plane (std::size_t index) const { return Plane (m_centroids[index], m_smallest_eigenvectors[index]); }
+
+  /*!
+    \brief Returns the normalized eigenvalues of the point at position `index`.
+  */
+  const Eigenvalues& eigenvalue (std::size_t index) const { return m_eigenvalues[index]; }
+
+  /*!
+    \brief Returns the sum of eigenvalues of the point at position `index`.
+  */
+  const double& sum_of_eigenvalues (std::size_t index) const { return m_sum_eigenvalues[index]; }
+
+  /// \cond SKIP_IN_MANUAL
+  double mean_range() const { return m_mean_range; }
+  /// \endcond
+
+private:
+
+  void compute (std::size_t index, const Point& query, std::vector<Point>& neighbor_points)
+  {
+    if (neighbor_points.size() == 0)
+      {
+        Eigenvalues v = {{ 0., 0., 0. }};
+        m_eigenvalues[index] = v;
+        m_centroids[index] = query;
+        m_smallest_eigenvectors[index] = Vector (0., 0., 1.);
+        m_middle_eigenvectors[index] = Vector (0., 1., 0.);
+        m_largest_eigenvectors[index] = Vector (1., 0., 0.);
+        return;
+      }
+
+    Point centroid = CGAL::centroid (neighbor_points.begin(), neighbor_points.end());
+    m_centroids[index] = centroid;
+    
+    CGAL::cpp11::array<double, 6> covariance = {{ 0., 0., 0., 0., 0., 0. }};
+      
+    for (std::size_t i = 0; i < neighbor_points.size(); ++ i)
+      {
+        Vector d = neighbor_points[i] - centroid;
+        covariance[0] += d.x () * d.x ();
+        covariance[1] += d.x () * d.y ();
+        covariance[2] += d.x () * d.z ();
+        covariance[3] += d.y () * d.y ();
+        covariance[4] += d.y () * d.z ();
+        covariance[5] += d.z () * d.z ();
+      }
+
+    Eigenvalues evalues = {{ 0., 0., 0. }};
+    CGAL::cpp11::array<double, 9> evectors = {{ 0., 0., 0.,
+                                                0., 0., 0.,
+                                                0., 0., 0. }};
+
+    DiagonalizeTraits::diagonalize_selfadjoint_covariance_matrix
+      (covariance, evalues, evectors);
+
+    // Normalize
+    double sum = evalues[0] + evalues[1] + evalues[2];
+    if (sum > 0.)
+      for (std::size_t i = 0; i < 3; ++ i)
+        evalues[i] = evalues[i] / sum;
+    m_sum_eigenvalues[index] = sum;
+    m_eigenvalues[index] = evalues;
+    m_smallest_eigenvectors[index] = Vector (evectors[0], evectors[1], evectors[2]);
+    m_middle_eigenvectors[index] = Vector (evectors[3], evectors[4], evectors[5]);
+    m_largest_eigenvectors[index] = Vector (evectors[6], evectors[7], evectors[8]);
+
+  }
+
+};
+  
+
+}
+  
+}
+
+
+#endif // CGAL_CLASSIFICATION_LOCAL_EIGEN_ANALYSIS_H

Classification/include/CGAL/Classification/Planimetric_grid.h

@@ -0,0 +1,133 @@
+// Copyright (c) 2012 INRIA Sophia-Antipolis (France).
+// Copyright (c) 2017 GeometryFactory Sarl (France).
+// All rights reserved.
+//
+// This file is part of CGAL (www.cgal.org).
+// You can redistribute it and/or modify it under the terms of the GNU
+// General Public License as published by the Free Software Foundation,
+// either version 3 of the License, or (at your option) any later version.
+//
+// Licensees holding a valid commercial license may use this file in
+// accordance with the commercial license agreement provided with the software.
+//
+// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
+// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
+//
+// $URL$
+// $Id$
+//
+// Author(s)     : Simon Giraudot, Florent Lafarge
+
+#ifndef CGAL_CLASSIFICATION_PLANIMETRIC_GRID_H
+#define CGAL_CLASSIFICATION_PLANIMETRIC_GRID_H
+
+#include <vector>
+
+#include <CGAL/Classification/Image.h>
+
+namespace CGAL {
+
+namespace Classification {
+
+  /*!
+    \ingroup PkgClassificationDataStructures
+
+    \brief Class that precomputes a 2D planimetric grid.
+
+    The grid is composed of squared cells with a user-defined size,
+    each cell containing the list of indices of the points whose
+    projection along the Z-axis lies within this cell. The mapping
+    from each point to the cell it lies in is also stored.
+
+    \tparam Geom_traits model of \cgal Kernel.
+    \tparam PointRange model of `ConstRange`. Its iterator type is
+    `RandomAccessIterator`.
+    \tparam PointMap model of `ReadablePropertyMap` whose key
+    type is the value type of the iterator of `PointRange` and value type
+    is `Geom_traits::Point_3`.
+  */
+
+template <typename Geom_traits, typename PointRange, typename PointMap>
+class Planimetric_grid
+{
+public:
+  typedef typename Geom_traits::Point_3 Point_3;
+  typedef typename Geom_traits::Iso_cuboid_3 Iso_cuboid_3;
+
+private:
+  typedef Image<std::vector<std::size_t> > Image_indices;
+  typedef Image<bool> Image_bool;
+
+  Image_indices m_grid;
+  std::vector<std::size_t> m_x;
+  std::vector<std::size_t> m_y;
+  
+public:
+
+  /// \cond SKIP_IN_MANUAL
+  Planimetric_grid () { }
+  /// \endcond
+
+  /*
+    \brief Constructs a planimetric grid based on the input range.
+
+    \param input input range.
+    \param point_map property map to access the input points.
+    \param bbox bounding box of the input range.
+    \param grid_resolution resolution of the planimetric grid.
+  */
+  Planimetric_grid (const PointRange& input,
+                    PointMap point_map,
+                    const Iso_cuboid_3& bbox,
+                    const double grid_resolution)
+  {
+    std::size_t width = (std::size_t)((bbox.xmax() - bbox.xmin()) / grid_resolution) + 1;
+    std::size_t height = (std::size_t)((bbox.ymax() - bbox.ymin()) / grid_resolution) + 1;
+    m_grid = Image_indices (width, height);
+
+    for (std::size_t i = 0; i < input.size(); ++ i)
+      {
+        const Point_3& p = get(point_map, *(input.begin()+i));
+        m_x.push_back ((std::size_t)((p.x() - bbox.xmin()) / grid_resolution));
+        m_y.push_back ((std::size_t)((p.y() - bbox.ymin()) / grid_resolution));
+        m_grid(m_x.back(), m_y.back()).push_back (i);
+      }
+
+  }
+
+  /*!
+    \brief Returns the number of cells along the X-axis.
+  */
+  std::size_t width() const { return m_grid.width(); }
+  /*!
+    \brief Returns the number of cells along the Y-axis.
+  */
+  std::size_t height() const { return m_grid.height(); }
+
+  /*!
+    \brief Returns the indices of the points lying in the cell at position `(x,y)`.
+  */
+  const std::vector<std::size_t>& indices(std::size_t x, std::size_t y) const { return m_grid(x,y); }
+  
+  /*!
+    \brief Returns `false` if the cell at position `(x,y)` is empty, `true` otherwise.
+  */
+  bool mask(std::size_t x, std::size_t y) const { return (!(m_grid(x,y).empty())); }
+
+  /*!
+    \brief Returns the `x` grid coordinate of the point at position `index`.
+  */
+  std::size_t x(std::size_t index) const { return m_x[index]; }
+  /*!
+    \brief Returns the `y` grid coordinate of the point at position `index`.
+  */
+  std::size_t y(std::size_t index) const { return m_y[index]; }
+};
+  
+
+}
+  
+}
+
+
+#endif // CGAL_CLASSIFICATION_PLANIMETRIC_GRID_H

Classification/include/CGAL/Classification/Point_set_neighborhood.h

@@ -0,0 +1,316 @@
+// Copyright (c) 2017 GeometryFactory Sarl (France).
+// All rights reserved.
+//
+// This file is part of CGAL (www.cgal.org).
+// You can redistribute it and/or modify it under the terms of the GNU
+// General Public License as published by the Free Software Foundation,
+// either version 3 of the License, or (at your option) any later version.
+//
+// Licensees holding a valid commercial license may use this file in
+// accordance with the commercial license agreement provided with the software.
+//
+// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
+// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
+//
+// $URL$
+// $Id$
+//
+// Author(s)     : Simon Giraudot
+
+#ifndef CGAL_CLASSIFICATION_POINT_SET_NEIGHBORHOOD_H
+#define CGAL_CLASSIFICATION_POINT_SET_NEIGHBORHOOD_H
+
+#include <vector>
+
+#include <boost/iterator/counting_iterator.hpp>
+
+#include <CGAL/Search_traits_3.h>
+#include <CGAL/Fuzzy_sphere.h>
+#include <CGAL/Orthogonal_k_neighbor_search.h>
+#include <CGAL/Default_diagonalize_traits.h>
+#include <CGAL/centroid.h>
+#include <CGAL/grid_simplify_point_set.h>
+
+#include <CGAL/Classification/Image.h>
+
+namespace CGAL {
+
+namespace Classification {
+
+  /*!
+    \ingroup PkgClassificationDataStructures
+
+    \brief Class that precomputes spatial searching structures for an
+    input point set and gives access to the local neighborhood of a
+    point as a set of indices.
+
+    It allows the user to generate models of `NeighborQuery` based on
+    a fixed range neighborhood or on a fixed K number of neighbors. In
+    addition, the spatial searching structures can be computed on a
+    simplified version of the point set to allow for neighbor queries
+    at a higher scale.
+
+    \tparam Geom_traits is a model of \cgal Kernel.
+    \tparam PointRange model of `ConstRange`. Its iterator type is
+    `RandomAccessIterator`.
+    \tparam PointMap model of `ReadablePropertyMap` whose key
+    type is the value type of the iterator of `PointRange` and value type
+    is `Geom_traits::Point_3`.
+  */
+template <typename Geom_traits, typename PointRange, typename PointMap>
+class Point_set_neighborhood
+{
+  
+  typedef typename Geom_traits::FT FT;
+  typedef typename Geom_traits::Point_3 Point;
+  
+  class My_point_property_map{
+    const PointRange* input;
+    PointMap point_map;
+    
+  public:
+    typedef Point value_type;
+    typedef const value_type& reference;
+    typedef std::size_t key_type;
+    typedef boost::lvalue_property_map_tag category;
+    My_point_property_map () { }
+    My_point_property_map (const PointRange *input, PointMap point_map)
+      : input (input), point_map (point_map) { }
+    reference operator[] (key_type k) const { return get(point_map, *(input->begin()+k)); }
+    friend inline reference get (const My_point_property_map& ppmap, key_type i) 
+    { return ppmap[i]; }
+  };
+
+  typedef Search_traits_3<Geom_traits> SearchTraits_3;
+  typedef Search_traits_adapter <std::size_t, My_point_property_map, SearchTraits_3> Search_traits;
+  typedef Sliding_midpoint<Search_traits> Splitter;
+  typedef Distance_adapter<std::size_t, My_point_property_map, Euclidean_distance<SearchTraits_3> > Distance;
+  typedef Kd_tree<Search_traits, Splitter, Tag_true> Tree;
+  typedef Fuzzy_sphere<Search_traits> Sphere;
+  typedef Orthogonal_k_neighbor_search<Search_traits, Distance, Splitter, Tree> Knn;
+
+
+  Tree* m_tree;
+  Distance m_distance;
+
+  std::vector<std::vector<std::size_t> > m_precomputed_neighbors;
+  
+public:
+
+  /*!
+    Functor that computes the neighborhood of an input point with a
+    fixed number of neighbors.
+
+    \cgalModels NeighborQuery
+  */
+  class K_neighbor_query
+  {
+  public:
+    typedef Point_set_neighborhood::Point value_type; ///<
+  private:
+    const Point_set_neighborhood& neighborhood;
+    std::size_t k;
+  public:
+    /*!
+      \brief Constructs a K neighbor query object.
+      \param neighborhood point set neighborhood object.
+      \param k number of neighbors per query.
+     */
+    K_neighbor_query (const Point_set_neighborhood& neighborhood, std::size_t k)
+      : neighborhood (neighborhood), k(k) { }
+
+    /// \cond SKIP_IN_MANUAL
+    template <typename OutputIterator>
+    OutputIterator operator() (const value_type& query, OutputIterator output) const
+    {
+      neighborhood.k_neighbors (query, k, output);
+      return output;
+    }
+    /// \endcond
+  };
+
+  /*!
+    Functor that computes the neighborhood of an input point defined
+    as the points lying in a sphere of fixed radius centered at the
+    input point.
+
+    \cgalModels NeighborQuery
+  */
+  class Range_neighbor_query
+  {
+  public:
+    typedef Point_set_neighborhood::Point value_type; ///<
+  private:
+    const Point_set_neighborhood& neighborhood;
+    double radius;
+  public:
+    /*!
+      \brief Constructs a range neighbor query object.
+      \param neighborhood point set neighborhood object.
+      \param radius radius of the neighbor query sphere.
+    */
+    Range_neighbor_query (const Point_set_neighborhood& neighborhood, double radius)
+      : neighborhood (neighborhood), radius(radius) { }
+
+    /// \cond SKIP_IN_MANUAL
+    template <typename OutputIterator>
+    OutputIterator operator() (const value_type& query, OutputIterator output) const
+    {
+      neighborhood.range_neighbors (query, radius, output);
+      return output;
+    }
+    /// \endcond
+  };
+
+  friend class K_neighbor_query;
+  friend class Range_neighbor_query;
+
+    /// \cond SKIP_IN_MANUAL
+  Point_set_neighborhood () : m_tree (NULL) { }
+  /// \endcond
+
+  /*!
+    \brief Constructs a neighborhood object based on the input range.
+
+    \param input input range.
+    \param point_map property map to access the input points.
+  */
+  Point_set_neighborhood (const PointRange& input,
+                          PointMap point_map)
+    : m_tree (NULL)
+  {
+    My_point_property_map pmap (&input, point_map);
+    m_tree = new Tree (boost::counting_iterator<std::size_t> (0),
+                       boost::counting_iterator<std::size_t> (input.size()),
+                       Splitter(),
+                       Search_traits (pmap));
+    m_distance = Distance (pmap);
+    m_tree->build();
+  }
+
+  /*!
+    \brief Constructs a simplified neighborhood object based on the input range.
+
+    This method first computes a simplified version of the input point
+    set by voxelization: a 3D grid is defined and for each subset
+    present in one cell, only the point closest to the centroid of
+    this subset is used.
+
+    \param input input range.
+    \param point_map property map to access the input points.
+    \param voxel_size size of the cells of the 3D grid used for simplification.
+  */
+  Point_set_neighborhood (const PointRange& input,
+                          PointMap point_map,
+                          double voxel_size)
+    : m_tree (NULL)
+  {
+    // First, simplify
+    std::vector<std::size_t> indices (input.size());
+    for (std::size_t i = 0; i < indices.size(); ++ i)
+      indices[i] = i;
+    My_point_property_map pmap (&input, point_map);
+
+    voxelize_point_set(indices, pmap, voxel_size);
+    
+    m_tree = new Tree (indices.begin(), indices.end(),
+                       Splitter(),
+                       Search_traits (pmap));
+    m_distance = Distance (pmap);
+    m_tree->build();
+  }
+
+  /// \cond SKIP_IN_MANUAL
+  ~Point_set_neighborhood ()
+  {
+    if (m_tree != NULL)
+      delete m_tree;
+  }
+  /// \endcond
+
+  /*!
+    \brief Returns a neighbor query object with fixed number of neighbors `k`.
+  */
+  K_neighbor_query k_neighbor_query (const std::size_t k) const
+  {
+    return K_neighbor_query (*this, k);
+  }
+
+  /*!
+    \brief Returns a neighbor query object with fixed radius `radius`.
+  */
+  Range_neighbor_query range_neighbor_query (const double radius) const
+  {
+    return Range_neighbor_query (*this, radius);
+  }
+
+private:
+
+  template <typename OutputIterator>
+  void range_neighbors (const Point& query, const FT radius_neighbors, OutputIterator output) const
+  {
+    CGAL_assertion (m_tree != NULL);
+    Sphere fs (query, radius_neighbors, 0, m_tree->traits());
+    m_tree->search (output, fs);
+  }
+
+  template <typename OutputIterator>
+  void k_neighbors (const Point& query, const std::size_t k, OutputIterator output) const
+  {
+    CGAL_assertion (m_tree != NULL);
+    Knn search (*m_tree, query, (unsigned int)k, 0, true, m_distance);
+    for (typename Knn::iterator it = search.begin(); it != search.end(); ++ it)
+      *(output ++) = it->first;
+  }
+
+  template <typename Map>
+  void voxelize_point_set (std::vector<std::size_t>& indices, Map point_map,
+                           double voxel_size)
+  {
+    std::map<Point, std::vector<std::size_t> > grid;
+
+    for (std::size_t i = 0; i < indices.size(); ++ i)
+      {
+        const Point& p = get(point_map, indices[i]);
+        Point ref (std::floor(p.x() / voxel_size),
+                   std::floor(p.y() / voxel_size),
+                   std::floor(p.z() / voxel_size));
+        typename std::map<Point, std::vector<std::size_t> >::iterator it;
+        boost::tie (it, boost::tuples::ignore)
+          = grid.insert (std::make_pair (ref, std::vector<std::size_t>()));
+        it->second.push_back (indices[i]);
+      }
+    indices.clear();
+    for (typename std::map<Point, std::vector<std::size_t> >::iterator
+           it = grid.begin(); it != grid.end(); ++ it)
+      {
+        const std::vector<std::size_t>& pts = it->second;
+        Point centroid = CGAL::centroid (boost::make_transform_iterator
+                                         (pts.begin(),
+                                          CGAL::Property_map_to_unary_function<Map>(point_map)),
+                                         boost::make_transform_iterator
+                                         (pts.end(),
+                                          CGAL::Property_map_to_unary_function<Map>(point_map)));
+        std::size_t chosen = 0;
+        double min_dist = (std::numeric_limits<double>::max)();
+        for (std::size_t i = 0; i < pts.size(); ++ i)
+          {
+            double dist = CGAL::squared_distance(get(point_map, pts[i]), centroid);
+            if (dist < min_dist)
+              {
+                min_dist = dist;
+                chosen = pts[i];
+              }
+          }
+        indices.push_back (chosen);
+      }
+  }
+};
+  
+
+}
+  
+}
+
+
+#endif // CGAL_CLASSIFICATION_POINT_SET_POINT_SET_NEIGHBORHOOD_H

Classification/include/CGAL/Classification/Trainer.h

@@ -0,0 +1,790 @@
+// Copyright (c) 2017 GeometryFactory Sarl (France).
+// All rights reserved.
+//
+// This file is part of CGAL (www.cgal.org).
+// You can redistribute it and/or modify it under the terms of the GNU
+// General Public License as published by the Free Software Foundation,
+// either version 3 of the License, or (at your option) any later version.
+//
+// Licensees holding a valid commercial license may use this file in
+// accordance with the commercial license agreement provided with the software.
+//
+// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
+// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
+//
+// $URL$
+// $Id$
+//
+// Author(s)     : Simon Giraudot
+
+#ifndef CGAL_CLASSIFICATION_TRAINER_H
+#define CGAL_CLASSIFICATION_TRAINER_H
+
+#include <CGAL/Classification/Feature_base.h>
+#include <CGAL/Classification/Label.h>
+#include <CGAL/Classifier.h>
+
+#ifdef CGAL_LINKED_WITH_TBB
+#include <tbb/parallel_for.h>
+#include <tbb/blocked_range.h>
+#include <tbb/scalable_allocator.h>
+#include <tbb/mutex.h>
+#endif // CGAL_LINKED_WITH_TBB
+
+//#define CGAL_CLASSTRAINING_VERBOSE
+#if defined(CGAL_CLASSTRAINING_VERBOSE)
+#define CGAL_CLASSTRAINING_CERR std::cerr
+#else
+#define CGAL_CLASSTRAINING_CERR std::ostream(0)
+#endif
+
+namespace CGAL {
+
+namespace Classification {
+
+/*!
+\ingroup PkgClassification
+
+\brief Training algorithm to set up weights and effects of the
+features used for classification.
+
+Each label must have ben given a small set of user-defined inliers to
+provide the training algorithm with a ground truth (see `set_inlier()`
+and `set_inliers()`).
+
+This methods estimates the set of feature weights and of
+[effects](@ref Classification::Feature::Effect) that make the
+classifier succeed in correctly classifying the sets of inliers given
+by the user. These parameters are directly modified within the
+`Classification::Feature_base` and `Classification::Label` objects.
+
+\tparam ItemRange model of `ConstRange`. Its iterator type is
+`RandomAccessIterator`.
+
+\tparam ItemMap model of `ReadablePropertyMap` whose key
+type is the value type of the iterator of `ItemRange` and value type is
+the type of the items that are classified.
+
+\tparam ConcurrencyTag enables sequential versus parallel
+algorithm. Possible values are `Parallel_tag` (default value is %CGAL
+is linked with TBB) or `Sequential_tag` (default value otherwise).
+*/
+template <typename ItemRange, typename ItemMap,
+#if defined(DOXYGEN_RUNNING)
+          typename ConcurrencyTag>
+#elif defined(CGAL_LINKED_WITH_TBB)
+          typename ConcurrencyTag = CGAL::Parallel_tag>
+#else
+          typename ConcurrencyTag = CGAL::Sequential_tag>
+#endif
+class Trainer
+{
+
+public:
+  typedef CGAL::Classifier<ItemRange, ItemMap, ConcurrencyTag>  Classifier;
+  typedef typename Classification::Label_handle                 Label_handle;
+  typedef typename Classification::Feature_handle               Feature_handle;
+
+private:
+
+#ifdef CGAL_LINKED_WITH_TBB
+  class Compute_worst_score_and_confidence
+  {
+    std::vector<std::size_t>& m_training_set;
+    Classifier* m_classifier;
+    std::size_t m_label;
+    double& m_confidence;
+    std::size_t& m_nb_okay;
+    tbb::mutex& m_mutex;
+    
+  public:
+
+    Compute_worst_score_and_confidence (std::vector<std::size_t>& training_set,
+                                        Classifier* classifier,
+                                        std::size_t label,
+                                        double& confidence,
+                                        std::size_t& nb_okay,
+                                        tbb::mutex& mutex)
+      : m_training_set (training_set)
+      , m_classifier (classifier)
+      , m_label (label)
+      , m_confidence (confidence)
+      , m_nb_okay (nb_okay)
+      , m_mutex (mutex)
+    { }
+          
+    void operator()(const tbb::blocked_range<std::size_t>& r) const
+    {
+      for (std::size_t k = r.begin(); k != r.end(); ++ k)
+        {
+          std::vector<std::pair<double, std::size_t> > values;
+      
+          for(std::size_t l = 0; l < m_classifier->number_of_labels(); ++ l)
+            values.push_back (std::make_pair (m_classifier->energy_of (m_classifier->label(l),
+                                                                       m_training_set[k]),
+                                              l));
+          std::sort (values.begin(), values.end());
+
+          if (values[0].second == m_label)
+            {
+              m_mutex.lock();
+              m_confidence += values[1].first - values[0].first;
+              ++ m_nb_okay;
+              m_mutex.unlock();
+            }
+        }
+    }
+
+  };
+#endif // CGAL_LINKED_WITH_TBB
+
+  
+  Classifier* m_classifier;
+  std::vector<std::vector<std::size_t> > m_training_sets;
+  std::vector<double> m_precision;
+  std::vector<double> m_recall;
+  std::vector<double> m_iou; // intersection over union
+  double m_accuracy;
+  double m_mean_iou;
+  double m_mean_f1;
+
+public:
+
+  /// \name Constructor
+  /// @{
+  
+  Trainer (Classifier& classifier)
+    : m_classifier (&classifier)
+    , m_training_sets (classifier.number_of_labels())
+  {
+  }
+
+  /// @}
+
+
+  /// \name Inliers
+  /// @{
+  
+  /*!
+    \brief Adds the item at position `index` as an inlier of
+    `label` for the training algorithm.
+
+    \note This inlier is only used for training. There is no guarantee
+    that the item at position `index` will be classified as `label`
+    after calling `run()`, `run_with_local_smoothing()` or
+    `run_with_graphcut()`.
+
+    \return `true` if the inlier was correctly added, `false`
+    otherwise (if `label` was not found).
+  */
+  bool set_inlier (Label_handle label, std::size_t index)
+  {
+    std::size_t label_idx = (std::size_t)(-1);
+    for (std::size_t i = 0; i < m_classifier->number_of_labels(); ++ i)
+      if (m_classifier->label(i) == label)
+        {
+          label_idx = i;
+          break;
+        }
+    if (label_idx == (std::size_t)(-1))
+      return false;
+
+    if (label_idx >= m_training_sets.size())
+      m_training_sets.resize (label_idx + 1);
+
+    m_training_sets[label_idx].push_back (index);
+
+    return true;
+  }
+
+  /*!
+
+    \brief Adds the items at positions `indices` as inliers of
+    `label` for the training algorithm.
+
+    \note These inliers are only used for training. There is no
+    guarantee that the items at positions `indices` will be classified
+    as `label` after calling `run()`,
+    `run_with_local_smoothing()` or `run_with_graphcut()`.
+
+    \tparam IndexRange range of `std::size_t`, model of `ConstRange`.
+  */
+  template <class IndexRange>
+  bool set_inliers (Label_handle label,
+                    IndexRange indices)
+  {
+    std::size_t label_idx = (std::size_t)(-1);
+    for (std::size_t i = 0; i < m_classifier->number_of_labels(); ++ i)
+      if (m_classifier->label(i) == label)
+        {
+          label_idx = i;
+          break;
+        }
+    if (label_idx == (std::size_t)(-1))
+      return false;
+
+    if (label_idx >= m_training_sets.size())
+      m_training_sets.resize (label_idx + 1);
+
+    std::copy (indices.begin(), indices.end(), std::back_inserter (m_training_sets[label_idx]));
+
+    return true;
+  }
+  
+  /*!
+    \brief Resets inlier sets used for training.
+  */
+  void reset_inlier_sets()
+  {
+    std::vector<std::vector<std::size_t > >().swap (m_training_sets);
+  }
+
+  
+  /// @}
+
+
+  /// \name Training
+  /// @{
+
+  /*!
+    \brief Runs the training algorithm.
+
+    All the `Classification::Label` and `Classification::Feature`
+    necessary for classification should have been added before running
+    this function. After training, the user can call `run()`,
+    `run_with_local_smoothing()` or `run_with_graphcut()` to compute
+    the classification using the estimated parameters.
+
+    \param nb_tests number of tests to perform. Higher values may
+    provide the user with better results at the cost of a higher
+    computation time. Using a value of at least 10 times the number of
+    features is advised.
+
+    \return minimum ratio (over all labels) of provided
+    ground truth items correctly classified using the best
+    configuration found.
+  */
+  
+  double train (std::size_t nb_tests = 300)
+  {
+    if (m_training_sets.empty())
+      return 0.;
+
+    for (std::size_t i = 0; i < m_classifier->number_of_labels(); ++ i)
+      if (m_training_sets.size() <= i || m_training_sets[i].empty())
+        std::cerr << "WARNING: \"" << m_classifier->label(i)->name() << "\" doesn't have a training set." << std::endl;
+
+    std::vector<double> best_weights (m_classifier->number_of_features(), 1.);
+
+    struct Feature_training
+    {
+      std::size_t i;
+      double wmin;
+      double wmax;
+      double factor;
+
+      bool operator<(const Feature_training& other) const
+      {
+        return (wmin / wmax) < (other.wmin / other.wmax);
+      }
+    };
+    std::vector<Feature_training> feature_train;
+    std::size_t nb_trials = 100;
+    double wmin = 1e-5, wmax = 1e5;
+    double factor = std::pow (wmax/wmin, 1. / (double)nb_trials);
+    
+    for (std::size_t j = 0; j < m_classifier->number_of_features(); ++ j)
+      {
+        Feature_handle feature = m_classifier->feature(j);
+        best_weights[j] = feature->weight();
+
+        std::size_t nb_useful = 0;
+        double min = (std::numeric_limits<double>::max)();
+        double max = -(std::numeric_limits<double>::max)();
+
+        feature->set_weight(wmin);
+        for (std::size_t i = 0; i < 100; ++ i)
+          {
+            estimate_feature_effect(feature);
+            if (feature_useful(feature))
+              {
+                CGAL_CLASSTRAINING_CERR << "#";
+                nb_useful ++;
+                min = (std::min) (min, feature->weight());
+                max = (std::max) (max, feature->weight());
+              }
+            else
+              CGAL_CLASSTRAINING_CERR << "-";
+            feature->set_weight(factor * feature->weight());
+          }
+        CGAL_CLASSTRAINING_CERR << std::endl;
+        CGAL_CLASSTRAINING_CERR << feature->name() << " useful in "
+                                << nb_useful << "% of the cases, in interval [ "
+                                << min << " ; " << max << " ]" << std::endl;
+        if (nb_useful < 2)
+          {
+            feature->set_weight(0.);
+            best_weights[j] = feature->weight();
+            continue;
+          }
+
+        feature_train.push_back (Feature_training());
+        feature_train.back().i = j;
+        feature_train.back().wmin = min / factor;
+        feature_train.back().wmax = max * factor;
+
+        if (best_weights[j] == 1.)
+          {
+            feature->set_weight(0.5 * (feature_train.back().wmin + feature_train.back().wmax));
+            best_weights[j] = feature->weight();
+          }
+        else
+          feature->set_weight(best_weights[j]);
+        estimate_feature_effect(feature);
+      }
+
+    std::size_t nb_trials_per_feature = 1 + (std::size_t)(nb_tests / (double)(feature_train.size()));
+    CGAL_CLASSIFICATION_CERR << "Trials = " << nb_tests << ", features = " << feature_train.size()
+              << ", trials per feature = " << nb_trials_per_feature << std::endl;
+    for (std::size_t i = 0; i < feature_train.size(); ++ i)
+      feature_train[i].factor
+        = std::pow (feature_train[i].wmax / feature_train[i].wmin,
+                    1. / (double)nb_trials_per_feature);
+    
+    
+    double best_score = 0.;
+    double best_confidence = 0.;
+    boost::tie (best_confidence, best_score)
+      = compute_worst_confidence_and_score (0., 0.);
+    
+    CGAL_CLASSIFICATION_CERR << "TRAINING GLOBALLY: Best score evolution: " << std::endl;
+
+    CGAL_CLASSIFICATION_CERR << 100. * best_score << "% (found at initialization)" << std::endl;
+
+    CGAL::Timer teff, tconf, tscore;
+
+    std::sort (feature_train.begin(), feature_train.end());
+    for (std::size_t i = 0; i < feature_train.size(); ++ i)
+      {
+        const Feature_training& tr = feature_train[i];
+        std::size_t current_feature_changed = tr.i;
+        Feature_handle current_feature = m_classifier->feature(current_feature_changed);
+        
+        std::size_t nb_used = 0;
+        for (std::size_t j = 0; j < m_classifier->number_of_features(); ++ j)
+          {
+            if (j == current_feature_changed)
+              continue;
+            
+            m_classifier->feature(j)->set_weight(best_weights[j]);
+            teff.start();
+            estimate_feature_effect(m_classifier->feature(j));
+            teff.stop();
+            if (feature_useful(m_classifier->feature(j)))
+              nb_used ++;
+            else
+              m_classifier->feature(j)->set_weight(0.);
+          }
+        
+        current_feature->set_weight(tr.wmin);
+        for (std::size_t j = 0; j < nb_trials_per_feature; ++ j)
+          {
+            teff.start();
+            estimate_feature_effect(current_feature);
+            teff.stop();
+
+            tconf.start();
+            double worst_confidence = 0., worst_score = 0.;
+            boost::tie (worst_confidence, worst_score)
+              = compute_worst_confidence_and_score (best_confidence, best_score);
+            tconf.stop();
+
+            if (worst_score > best_score
+                && worst_confidence > best_confidence)
+              {
+                best_score = worst_score;
+                best_confidence = worst_confidence;
+                CGAL_CLASSIFICATION_CERR << 100. * best_score << "% (found at iteration "
+                          << (i * nb_trials_per_feature) + j << "/" << nb_tests << ", "
+                          << nb_used + (feature_useful(current_feature) ? 1 : 0)
+                          << "/" << m_classifier->number_of_features() << " feature(s) used)" << std::endl;
+                for (std::size_t k = 0; k < m_classifier->number_of_features(); ++ k)
+                  {
+                    Feature_handle feature = m_classifier->feature(k);
+                    best_weights[k] = feature->weight();
+                  }
+              }
+            current_feature->set_weight(current_feature->weight() * tr.factor);
+          }
+      }
+
+    std::cerr << "Estimation of effects = " << teff.time() << std::endl
+              << "Confidence computation = " << tconf.time() << std::endl
+              << "Score computation = " << tscore.time() << std::endl;
+    
+    for (std::size_t i = 0; i < best_weights.size(); ++ i)
+      {
+        Feature_handle feature = m_classifier->feature(i);
+        feature->set_weight(best_weights[i]);
+      }
+
+    estimate_features_effects();
+    
+    CGAL_CLASSIFICATION_CERR << std::endl << "Best score found is at least " << 100. * best_score
+              << "% of correct classification" << std::endl;
+
+    std::size_t nb_removed = 0;
+    for (std::size_t i = 0; i < best_weights.size(); ++ i)
+      {
+        Feature_handle feature = m_classifier->feature(i);
+        CGAL_CLASSTRAINING_CERR << "FEATURE " << feature->name() << ": " << best_weights[i] << std::endl;
+        feature->set_weight(best_weights[i]);
+
+        Classification::Feature::Effect side = m_classifier->label(0)->feature_effect(feature);
+        bool to_remove = true;
+        for (std::size_t j = 0; j < m_classifier->number_of_labels(); ++ j)
+          {
+            Label_handle clabel = m_classifier->label(j);
+            if (clabel->feature_effect(feature) == Classification::Feature::FAVORING)
+              CGAL_CLASSTRAINING_CERR << " * Favored for ";
+            else if (clabel->feature_effect(feature) == Classification::Feature::PENALIZING)
+              CGAL_CLASSTRAINING_CERR << " * Penalized for ";
+            else
+              CGAL_CLASSTRAINING_CERR << " * Neutral for ";
+            if (clabel->feature_effect(feature) != side)
+              to_remove = false;
+            CGAL_CLASSTRAINING_CERR << clabel->name() << std::endl;
+          }
+        if (to_remove)
+          {
+            CGAL_CLASSTRAINING_CERR << "   -> Useless! Should be removed" << std::endl;
+            ++ nb_removed;
+          }
+      }
+    CGAL_CLASSIFICATION_CERR << nb_removed
+              << " feature(s) out of " << m_classifier->number_of_features() << " are useless" << std::endl;
+
+    compute_precision_recall();
+    return best_score;
+  }
+
+  /// @}
+
+
+  /// \name Evaluation
+  /// @{
+
+  /*!
+
+    \brief Returns the precision of the training for the given label.
+
+    Precision is the number of true positives divided by the sum of
+    the true positives and the false positives.
+
+  */
+  double precision (Label_handle label) const
+  {
+    std::size_t label_idx = (std::size_t)(-1);
+    for (std::size_t i = 0; i < m_classifier->number_of_labels(); ++ i)
+      if (m_classifier->label(i) == label)
+        {
+          label_idx = i;
+          break;
+        }
+    if (label_idx == (std::size_t)(-1))
+      return 0.;
+
+    return m_precision[label_idx];
+  }
+
+  /*!
+
+    \brief Returns the recall of the training for the given label.
+
+    Recall is the number of true positives divided by the sum of
+    the true positives and the false negatives.
+
+  */
+  double recall (Label_handle label) const
+  {
+    std::size_t label_idx = (std::size_t)(-1);
+    for (std::size_t i = 0; i < m_classifier->number_of_labels(); ++ i)
+      if (m_classifier->label(i) == label)
+        {
+          label_idx = i;
+          break;
+        }
+    if (label_idx == (std::size_t)(-1))
+      return 0.;
+
+    return m_recall[label_idx];
+  }
+
+  /*!
+
+    \brief Returns the \f$F_1\f$ score of the training for the given label.
+
+    \f$F_1\f$ score is the harmonic mean of `precision()` and `recall()`:
+
+    \f[
+    F_1 = 2 \times \frac{precision \times recall}{precision + recall}
+    \f]
+
+  */
+  double f1_score (Label_handle label) const
+  {
+    std::size_t label_idx = (std::size_t)(-1);
+    for (std::size_t i = 0; i < m_classifier->number_of_labels(); ++ i)
+      if (m_classifier->label(i) == label)
+        {
+          label_idx = i;
+          break;
+        }
+    if (label_idx == (std::size_t)(-1))
+      return 0.;
+    
+    return 2. * (m_precision[label_idx] * m_recall[label_idx])
+      / (m_precision[label_idx] + m_recall[label_idx]);
+  }
+
+  /*!
+    \brief Returns the intersection over union of the training for the
+    given label.
+
+    Intersection over union is the number of true positives divided by
+    the sum of the true positives, of the false positives and of the
+    false negatives.
+  */
+  double intersection_over_union (Label_handle label) const
+  {
+    std::size_t label_idx = (std::size_t)(-1);
+    for (std::size_t i = 0; i < m_classifier->number_of_labels(); ++ i)
+      if (m_classifier->label(i) == label)
+        {
+          label_idx = i;
+          break;
+        }
+    if (label_idx == (std::size_t)(-1))
+      return 0.;
+    
+    return m_iou[label_idx];
+  }
+
+  /*!
+    \brief Returns the accuracy of the training.
+
+    Accuracy is the total number of true positives divided by the
+    total number of provided inliers.
+  */
+  double accuracy() const { return m_accuracy; }
+  
+  /*!
+    \brief Returns the mean \f$F_1\f$ score of the training over all
+    labels (see `f1_score()`).
+  */
+  double mean_f1_score() const { return m_mean_f1; }
+  
+  /*!
+    \brief Returns the mean intersection over union of the training
+    over all labels (see `intersection_over_union()`).
+  */
+  double mean_intersection_over_union() const { return m_mean_iou; }
+  
+  /// @}
+
+  /// \cond SKIP_IN_MANUAL
+  Label_handle training_label_of (std::size_t index) const
+  {
+    for (std::size_t i = 0; i < m_training_sets.size(); ++ i)
+      if (std::find (m_training_sets[i].begin(), m_training_sets[i].end(), index) != m_training_sets[i].end())
+        return m_classifier->label(i);
+    return Label_handle();
+  }
+
+  /// \endcond
+
+private:
+
+  void estimate_features_effects()
+  {
+    for (std::size_t i = 0; i < m_classifier->number_of_features(); ++ i)
+      estimate_feature_effect (m_classifier->feature(i));
+  }
+
+  void estimate_feature_effect (Feature_handle feature)
+  {
+    std::vector<double> mean (m_classifier->number_of_labels(), 0.);
+                                  
+    for (std::size_t j = 0; j < m_classifier->number_of_labels(); ++ j)
+      {
+        for (std::size_t k = 0; k < m_training_sets[j].size(); ++ k)
+          {
+            double val = feature->normalized(m_training_sets[j][k]);
+            mean[j] += val;
+          }
+        mean[j] /= m_training_sets[j].size();
+      }
+
+    std::vector<double> sd (m_classifier->number_of_labels(), 0.);
+        
+    for (std::size_t j = 0; j < m_classifier->number_of_labels(); ++ j)
+      {
+        Label_handle clabel = m_classifier->label(j);
+            
+        for (std::size_t k = 0; k < m_training_sets[j].size(); ++ k)
+          {
+            double val = feature->normalized(m_training_sets[j][k]);
+            sd[j] += (val - mean[j]) * (val - mean[j]);
+          }
+        sd[j] = std::sqrt (sd[j] / m_training_sets[j].size());
+        if (mean[j] - sd[j] > 0.75)
+          clabel->set_feature_effect (feature, Classification::Feature::FAVORING);
+        else if (mean[j] + sd[j] < 0.25)
+          clabel->set_feature_effect (feature, Classification::Feature::PENALIZING);
+        else
+          clabel->set_feature_effect (feature, Classification::Feature::NEUTRAL);
+      }
+  }
+
+  std::pair<double, double> compute_worst_confidence_and_score (double lower_conf, double lower_score)
+  {
+    double worst_confidence = (std::numeric_limits<double>::max)();
+    double worst_score = (std::numeric_limits<double>::max)();
+    
+    for (std::size_t j = 0; j < m_classifier->number_of_labels(); ++ j)
+      {
+        double confidence = 0.;
+        std::size_t nb_okay = 0;
+
+#ifndef CGAL_LINKED_WITH_TBB
+        CGAL_static_assertion_msg (!(boost::is_convertible<ConcurrencyTag, Parallel_tag>::value),
+                                   "Parallel_tag is enabled but TBB is unavailable.");
+#else
+        if (boost::is_convertible<ConcurrencyTag,Parallel_tag>::value)
+          {
+            tbb::mutex mutex;
+            Compute_worst_score_and_confidence f(m_training_sets[j], m_classifier, j, confidence, nb_okay, mutex);
+            tbb::parallel_for(tbb::blocked_range<size_t>(0, m_training_sets[j].size ()), f);
+          }
+        else
+#endif
+          {
+            for (std::size_t k = 0; k < m_training_sets[j].size(); ++ k)
+              {
+                std::vector<std::pair<double, std::size_t> > values;
+      
+                for(std::size_t l = 0; l < m_classifier->number_of_labels(); ++ l)
+                  values.push_back (std::make_pair (m_classifier->energy_of (m_classifier->label(l),
+                                                                             m_training_sets[j][k]),
+                                                    l));
+                std::sort (values.begin(), values.end());
+
+                if (values[0].second == j)
+                  {
+                    confidence += values[1].first - values[0].first;
+                    ++ nb_okay;
+                  }
+              }
+          }
+        
+        double score = nb_okay / (double)(m_training_sets[j].size());
+        confidence /= (double)(m_training_sets[j].size() * m_classifier->number_of_features());
+
+        if (confidence < worst_confidence)
+          worst_confidence = confidence;
+        if (score < worst_score)
+          worst_score = score;
+        
+        if (worst_confidence < lower_conf || worst_score < lower_score)
+          return std::make_pair (worst_confidence, worst_score);
+      }
+    return std::make_pair (worst_confidence, worst_score);
+  }
+
+  bool feature_useful (Feature_handle feature)
+  {
+    Classification::Feature::Effect side = m_classifier->label(0)->feature_effect(feature);
+    for (std::size_t k = 1; k < m_classifier->number_of_labels(); ++ k)
+      if (m_classifier->label(k)->feature_effect(feature) != side)
+        return true;
+    return false;
+  }
+
+  
+  void compute_precision_recall ()
+  {
+    std::vector<std::size_t> true_positives (m_classifier->number_of_labels());
+    std::vector<std::size_t> false_positives (m_classifier->number_of_labels());
+    std::vector<std::size_t> false_negatives (m_classifier->number_of_labels());
+
+    std::size_t sum_true_positives = 0;
+    std::size_t total = 0;
+    
+    for (std::size_t j = 0; j < m_classifier->number_of_labels(); ++ j)
+      {
+        for (std::size_t k = 0; k < m_training_sets[j].size(); ++ k)
+          {
+            std::size_t nb_class_best=0; 
+            double val_class_best = (std::numeric_limits<double>::max)();
+      
+            for(std::size_t l = 0; l < m_classifier->number_of_labels(); ++ l)
+              {
+                double value = m_classifier->energy_of (m_classifier->label(l),
+                                                        m_training_sets[j][k]);
+          
+                if(val_class_best > value)
+                  {
+                    val_class_best = value;
+                    nb_class_best = l;
+                  }
+              }
+            ++ total;
+            if (nb_class_best == j)
+              {
+                ++ true_positives[j];
+                ++ sum_true_positives;
+              }
+            else
+              {
+                ++ false_negatives[j];
+                for(std::size_t l = 0; l < m_classifier->number_of_labels(); ++ l)
+                  if (nb_class_best == l)
+                    ++ false_positives[l];
+              }
+          }
+      }
+
+
+    m_precision.clear();
+    m_recall.clear();
+
+    m_mean_iou = 0.;
+    m_mean_f1 = 0.;
+    
+    for (std::size_t j = 0; j < m_classifier->number_of_labels(); ++ j)
+      {
+        m_precision.push_back (true_positives[j] / double(true_positives[j] + false_positives[j]));
+        m_recall.push_back (true_positives[j] / double(true_positives[j] + false_negatives[j]));
+        m_iou.push_back (true_positives[j] / double(true_positives[j] + false_positives[j] + false_negatives[j]));
+
+        m_mean_iou += m_iou.back();
+        m_mean_f1 += 2. * (m_precision.back() * m_recall.back())
+          / (m_precision.back() + m_recall.back());
+      }
+
+    m_mean_iou /= m_classifier->number_of_labels();
+    m_mean_f1 /= m_classifier->number_of_labels();
+    m_accuracy = sum_true_positives / double(total);
+
+  }
+
+  
+
+};
+
+
+} // namespace Classification
+  
+} // namespace CGAL
+
+#endif // CGAL_CLASSIFICATION_TRAINER_H

Classification/test/Classification/CMakeLists.txt

@@ -0,0 +1,57 @@
+# Created by the script cgal_create_CMakeLists
+# This is the CMake script for compiling a set of CGAL applications.
+
+project( Classification )
+
+
+cmake_minimum_required(VERSION 2.8.11)
+
+# CGAL and its components
+find_package( CGAL QUIET COMPONENTS  )
+
+if ( NOT CGAL_FOUND )
+
+  message(STATUS "This project requires the CGAL library, and will not be compiled.")
+  return()  
+
+endif()
+
+# include helper file
+include( ${CGAL_USE_FILE} )
+
+
+# Boost and its components
+find_package( Boost REQUIRED )
+
+# Use Eigen
+find_package(Eigen3 3.1.0) #(requires 3.1.0 or greater)
+if (EIGEN3_FOUND)
+  include( ${EIGEN3_USE_FILE} )
+endif()
+
+if ( NOT Boost_FOUND )
+
+  message(STATUS "This project requires the Boost library, and will not be compiled.")
+
+  return()  
+
+endif()
+
+# include for local directory
+include_directories( BEFORE include )
+
+# include for local package
+include_directories( BEFORE ../../include )
+
+
+# Creating entries for all C++ files with "main" routine
+# ##########################################################
+
+include( CGAL_CreateSingleSourceCGALProgram )
+
+#add_definitions("-DCGAL_DO_NOT_USE_BOYKOV_KOLMOGOROV_MAXFLOW_SOFTWARE")
+
+
+create_single_source_cgal_program( "test_classifier.cpp" )
+create_single_source_cgal_program( "test_point_set_classifier.cpp" )
+

Classification/test/Classification/test_classifier.cpp

@@ -0,0 +1,176 @@
+#include <cstdlib>
+#include <fstream>
+#include <iostream>
+#include <string>
+
+#include <CGAL/Simple_cartesian.h>
+#include <CGAL/Random.h>
+#include <CGAL/tuple.h>
+#include <CGAL/Classifier.h>
+#include <CGAL/Classification/Local_eigen_analysis.h>
+#include <CGAL/Classification/Point_set_neighborhood.h>
+#include <CGAL/Classification/Planimetric_grid.h>
+#include <CGAL/Classification/Attribute_base.h>
+#include <CGAL/Classification/Attribute/Distance_to_plane.h>
+#include <CGAL/Classification/Attribute/Echo_scatter.h>
+#include <CGAL/Classification/Attribute/Eigen.h>
+#include <CGAL/Classification/Attribute/Elevation.h>
+#include <CGAL/Classification/Attribute/Hsv.h>
+#include <CGAL/Classification/Attribute/Vertical_dispersion.h>
+#include <CGAL/Classification/Attribute/Verticality.h>
+
+typedef CGAL::Simple_cartesian<double> Kernel;
+typedef Kernel::Point_3 Point;
+typedef Kernel::Vector_3 Vector;
+typedef CGAL::cpp11::array<unsigned char, 3> Color;
+typedef Kernel::Iso_cuboid_3 Iso_cuboid_3;
+
+typedef boost::tuple<Point, Vector, Color, std::size_t> Point_with_info;
+typedef std::vector<Point_with_info> Point_range;
+
+typedef CGAL::Nth_of_tuple_property_map<0, Point_with_info> Pmap;
+typedef CGAL::Nth_of_tuple_property_map<1, Point_with_info> Nmap;
+typedef CGAL::Nth_of_tuple_property_map<2, Point_with_info> Cmap;
+typedef CGAL::Nth_of_tuple_property_map<3, Point_with_info> Emap;
+
+typedef CGAL::Classifier<Point_range, Pmap> Classifier;
+
+typedef CGAL::Classification::Planimetric_grid<Kernel, Point_range, Pmap>       Planimetric_grid;
+typedef CGAL::Classification::Point_set_neighborhood<Kernel, Point_range, Pmap> Neighborhood;
+typedef CGAL::Classification::Local_eigen_analysis<Kernel, Point_range, Pmap>   Local_eigen_analysis;
+
+typedef CGAL::Classification::Type_handle                                           Type_handle;
+typedef CGAL::Classification::Attribute_handle                                      Attribute_handle;
+
+typedef CGAL::Classification::Attribute::Anisotropy
+<Kernel, Point_range, Pmap> Anisotropy;
+typedef CGAL::Classification::Attribute::Distance_to_plane
+<Kernel, Point_range, Pmap> Distance_to_plane;
+typedef CGAL::Classification::Attribute::Eigentropy
+<Kernel, Point_range, Pmap> Eigentropy;
+typedef CGAL::Classification::Attribute::Elevation
+<Kernel, Point_range, Pmap>                    Elevation;
+typedef CGAL::Classification::Attribute::Linearity
+<Kernel, Point_range, Pmap> Linearity;
+typedef CGAL::Classification::Attribute::Omnivariance
+<Kernel, Point_range, Pmap> Omnivariance;
+typedef CGAL::Classification::Attribute::Planarity
+<Kernel, Point_range, Pmap> Planarity;
+typedef CGAL::Classification::Attribute::Sphericity
+<Kernel, Point_range, Pmap> Sphericity;
+typedef CGAL::Classification::Attribute::Sum_eigenvalues
+<Kernel, Point_range, Pmap> Sum_eigen;
+typedef CGAL::Classification::Attribute::Surface_variation
+<Kernel, Point_range, Pmap> Surface_variation;
+typedef CGAL::Classification::Attribute::Vertical_dispersion
+<Kernel, Point_range, Pmap>                    Dispersion;
+typedef CGAL::Classification::Attribute::Verticality
+<Kernel, Point_range, Pmap> Verticality;
+typedef CGAL::Classification::RGB_Color RGB_Color;
+typedef CGAL::Classification::Attribute::Hsv
+<Kernel, Point_range, Cmap> Hsv;
+typedef CGAL::Classification::Attribute::Hsv
+<Kernel, Point_range, Cmap> Hsv;
+typedef CGAL::Classification::Attribute::Echo_scatter
+<Kernel, Point_range, Pmap, Emap> Echo_scatter;
+
+typedef CGAL::Classification::Attribute::Effect Attribute_effect;
+
+int main (int, char**)
+{
+  Point_range range;
+  range.reserve(10000);
+  for (std::size_t i = 0; i < 10000; ++ i)
+    range.push_back (boost::make_tuple
+                     (Point (CGAL::get_default_random().get_double(),
+                             CGAL::get_default_random().get_double(),
+                             CGAL::get_default_random().get_double()),
+                      Vector (CGAL::get_default_random().get_double(),
+                              CGAL::get_default_random().get_double(),
+                              CGAL::get_default_random().get_double()),
+                      CGAL::make_array ((unsigned char)(CGAL::get_default_random().get_int(0, 255)),
+                                        (unsigned char)(CGAL::get_default_random().get_int(0, 255)),
+                                        (unsigned char)(CGAL::get_default_random().get_int(0, 255))),
+                      std::size_t(CGAL::get_default_random().get_int(0, 10))));
+
+  double grid_resolution = 0.34;
+  double radius_neighbors = 1.7;
+  double radius_dtm = 15.0;
+  
+  Classifier classifier (range, Pmap());
+  Iso_cuboid_3 bbox = CGAL::bounding_box (boost::make_transform_iterator
+                                          (range.begin(),
+                                           CGAL::Property_map_to_unary_function<Pmap>()),
+                                          boost::make_transform_iterator
+                                          (range.end(),
+                                           CGAL::Property_map_to_unary_function<Pmap>()));
+
+  Planimetric_grid grid (range, Pmap(), bbox, grid_resolution);
+  Neighborhood neighborhood (range, Pmap());
+  Local_eigen_analysis eigen (range, Pmap(), neighborhood.k_neighbor_query(6));
+
+  std::vector<Attribute_handle> att;
+  
+  att.push_back (classifier.add_attribute<Anisotropy> (eigen));
+  att.push_back (classifier.add_attribute<Dispersion> (Pmap(), grid,
+                                                       grid_resolution,
+                                                       radius_neighbors));
+  att.push_back (classifier.add_attribute<Distance_to_plane> (Pmap(), eigen));
+  att.push_back (classifier.add_attribute<Echo_scatter> (Emap(), grid, grid_resolution, radius_neighbors));
+  att.push_back (classifier.add_attribute<Eigentropy> (eigen));
+  att.push_back (classifier.add_attribute<Elevation> (Pmap(), grid,
+                                                      grid_resolution,
+                                                      radius_dtm));
+  att.push_back (classifier.add_attribute<Hsv> (Cmap(), 0, 180, 90));
+  att.push_back (classifier.add_attribute<Hsv> (Cmap(), 1, 50, 25));
+  att.push_back (classifier.add_attribute<Hsv> (Cmap(), 2, 50, 25));
+  att.push_back (classifier.add_attribute<Linearity> (eigen));
+  att.push_back (classifier.add_attribute<Omnivariance> (eigen));
+  att.push_back (classifier.add_attribute<Planarity> (eigen));
+  att.push_back (classifier.add_attribute<Sphericity> (eigen));
+  att.push_back (classifier.add_attribute<Sum_eigen> (eigen));
+  att.push_back (classifier.add_attribute<Surface_variation> (eigen));
+  att.push_back (classifier.add_attribute<Verticality> (eigen));
+  att.push_back (classifier.add_attribute<Verticality> (Nmap()));
+  att.push_back (classifier.add_attribute<Omnivariance> (eigen));
+  assert (classifier.remove_attribute (att.back()));
+  att.pop_back();
+  assert (att.size() == classifier.number_of_attributes());
+  
+  std::vector<Type_handle> types;
+  types.push_back (classifier.add_classification_type ("type1"));
+  types.push_back (classifier.add_classification_type ("type2"));
+  types.push_back (classifier.add_classification_type ("type3"));
+  types.push_back (classifier.add_classification_type ("type4"));
+  types.push_back (classifier.add_classification_type ("type5"));
+  types.push_back (classifier.add_classification_type ("type6"));
+  assert (classifier.remove_classification_type (types.back()));
+  types.pop_back();
+  assert (types.size() == classifier.number_of_classification_types());
+
+  for (std::size_t i = 0; i < types.size(); ++ i)
+    for (std::size_t j = 0; j < att.size(); ++ j)
+      {
+        att[j]->set_weight (CGAL::get_default_random().get_double());
+        types[i]->set_attribute_effect
+          (att[j], (Attribute_effect)(CGAL::get_default_random().get_int(0,3)));
+      }
+
+  classifier.run();
+  classifier.run_with_local_smoothing(neighborhood.k_neighbor_query(6));
+  classifier.run_with_graphcut(neighborhood.k_neighbor_query(12), 0.5);
+
+  std::vector<Type_handle> output;
+  output.reserve (10000);
+  for (std::size_t i = 0; i < range.size(); ++ i)
+    {
+      output.push_back (classifier.classification_type_of(i));
+      assert (0 <= classifier.confidence_of(i));
+    }
+  
+
+  classifier.clear_classification_types();
+  classifier.clear_attributes();
+  
+  return EXIT_SUCCESS;
+}

Classification/test/Classification/test_point_set_classifier.cpp

@@ -0,0 +1,81 @@
+#include <cstdlib>
+#include <fstream>
+#include <iostream>
+#include <string>
+
+#include <CGAL/Simple_cartesian.h>
+#include <CGAL/Random.h>
+#include <CGAL/tuple.h>
+#include <CGAL/Point_set_classifier.h>
+
+typedef CGAL::Simple_cartesian<double> Kernel;
+typedef Kernel::Point_3 Point;
+typedef Kernel::Vector_3 Vector;
+typedef CGAL::cpp11::array<unsigned char, 3> Color;
+typedef Kernel::Iso_cuboid_3 Iso_cuboid_3;
+
+typedef boost::tuple<Point, Vector, Color, std::size_t> Point_with_info;
+typedef std::vector<Point_with_info> Point_range;
+
+typedef CGAL::Nth_of_tuple_property_map<0, Point_with_info> Pmap;
+typedef CGAL::Nth_of_tuple_property_map<1, Point_with_info> Nmap;
+typedef CGAL::Nth_of_tuple_property_map<2, Point_with_info> Cmap;
+typedef CGAL::Nth_of_tuple_property_map<3, Point_with_info> Emap;
+
+typedef CGAL::Point_set_classifier<Kernel, Point_range, Pmap> PSC;
+
+typedef CGAL::Classification::Type_handle                                           Type_handle;
+typedef CGAL::Classification::Attribute_handle                                      Attribute_handle;
+
+typedef CGAL::Classification::Attribute::Effect Attribute_effect;
+
+int main (int, char**)
+{
+  Point_range range;
+  range.reserve(10000);
+  for (std::size_t i = 0; i < 10000; ++ i)
+    range.push_back (boost::make_tuple
+                     (Point (CGAL::get_default_random().get_double(),
+                             CGAL::get_default_random().get_double(),
+                             CGAL::get_default_random().get_double()),
+                      Vector (CGAL::get_default_random().get_double(),
+                              CGAL::get_default_random().get_double(),
+                              CGAL::get_default_random().get_double()),
+                      CGAL::make_array ((unsigned char)(CGAL::get_default_random().get_int(0, 255)),
+                                        (unsigned char)(CGAL::get_default_random().get_int(0, 255)),
+                                        (unsigned char)(CGAL::get_default_random().get_int(0, 255))),
+                      std::size_t(CGAL::get_default_random().get_int(0, 10))));
+
+  PSC classifier (range, Pmap());
+  classifier.generate_attributes (5, Nmap(), Cmap(), Emap());
+  assert (classifier.number_of_scales() == 5);
+
+  std::vector<Type_handle> types;
+  types.push_back (classifier.add_classification_type ("type1"));
+  types.push_back (classifier.add_classification_type ("type2"));
+  types.push_back (classifier.add_classification_type ("type3"));
+  types.push_back (classifier.add_classification_type ("type4"));
+  types.push_back (classifier.add_classification_type ("type5"));
+
+  for (std::size_t i = 0; i < 100; ++ i)
+    classifier.set_inlier (types[CGAL::get_default_random().get_int(0, types.size())], i);
+
+  classifier.train(500);
+
+
+  std::ofstream fconfig ("config.xml");
+  classifier.save_configuration (fconfig);
+  fconfig.close();
+  
+  std::ofstream fclassif ("classif.ply");
+  classifier.write_classification_to_ply(fclassif);
+  fclassif.close();
+  
+  classifier.clear();
+
+  std::ifstream fconfig2 ("config.xml");
+  assert (classifier.load_configuration (fconfig2, Nmap(), Cmap(), Emap()));
+  
+  
+  return EXIT_SUCCESS;
+}

Documentation/doc/Documentation/dependencies

@@ -97,6 +97,6 @@ Barycentric_coordinates_2
 Surface_mesh
 Surface_mesh_shortest_path
 Polygon_mesh_processing
-
+Classification
 
 

Documentation/doc/Documentation/packages.txt

@@ -120,6 +120,7 @@ h1 {
 \package_listing{Point_set_processing_3}
 \package_listing{Point_set_shape_detection_3}
 \package_listing{Stream_lines_2}
+\package_listing{Classification}
 
 
 \section PartSearchStructures Spatial Searching and Sorting

Documentation/doc/biblio/cgal_manual.bib

@@ -998,6 +998,15 @@ note         = {\url{ttp://hal.inria.fr/inria-00090522}}
   , pages       = "209--225"
 }
 
+@article{cgal:hws-fsso3-16,
+  title={Fast semantic segmentation of 3D point clouds with strongly varying density},
+  author={Hackel, Timo and Wegner, Jan D and Schindler, Konrad},
+  journal={ISPRS Annals of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Prague, Czech Republic},
+  volume={3},
+  pages={177--184},
+  year={2016}
+}
+
 @book{ cgal:hw-vrml2h-96
   ,author       = {Jed Hartman and Josie Wernecke}
   ,title        = {The {VRML} 2.0 Handbook: Building Moving Worlds on the
@@ -1059,6 +1068,15 @@ note         = {\url{ttp://hal.inria.fr/inria-00090522}}
    ,pages        = "307--320"
 }
 
+@book{ cgal:l-mrfmi-09,
+ author = {Li, Stan Z.},
+ title = {Markov Random %Field Modeling in Image Analysis},
+ year = {2009},
+ isbn = {9781848002784},
+ edition = {3rd},
+ publisher = {Springer Publishing Company, Incorporated}
+} 
+
 @inproceedings { cgal:l-nmdgp-05,
   AUTHOR = {Bruno Levy},
   TITLE  = {Numerical Methods for Digital Geometry Processing},
@@ -1089,6 +1107,15 @@ note         = {\url{ttp://hal.inria.fr/inria-00090522}}
     address = {Girona, Spain}
 }
 
+@article{cgal:lm-clscm-12,
+    author = {Lafarge, Florent and Mallet, Clement},
+    title = {{Creating large-scale city models from 3D-point clouds: a robust approach with hybrid representation}},
+    journal = {International Journal of Computer Vision},
+    volume = {99},
+    number = {1},
+    pages = {69-85},
+    year = {2012},
+}
 
 @inproceedings{ cgal:lt-fmeps-98,
     author = "Peter Lindstrom and Greg Turk",
@@ -1325,6 +1352,20 @@ ABSTRACT = {We present the first complete, exact and efficient C++ implementatio
   update        = "09.11 penarand"
 }
 
+@article{cgal:mbrsh-raofw-11,
+    title = "Relevance assessment of full-waveform lidar data for urban area classification ",
+    journal = "\{ISPRS\} Journal of Photogrammetry and Remote Sensing ",
+    volume = "66",
+    number = "6, Supplement",
+    pages = "S71 - S84",
+    year = "2011",
+    note = "Advances in \{LIDAR\} Data Processing and Applications ",
+    issn = "0924-2716",
+    doi = "http://dx.doi.org/10.1016/j.isprsjprs.2011.09.008",
+    url = "http://www.sciencedirect.com/science/article/pii/S0924271611001055",
+    author = "Clément Mallet and Frédéric Bretar and Michel Roux and Uwe Soergel and Christian Heipke"
+}
+
 @article{cgal:ml-cfsg-00
 , author =	"G. Medioni and M. Lee and C. Tang"
 , title =	"A Computational Framework for Segmentation and Grouping"

Documentation/doc/resources/1.8.4/BaseDoxyfile.in

@@ -1395,7 +1395,7 @@ EXTRA_SEARCH_MAPPINGS  =
 # If the GENERATE_LATEX tag is set to YES (the default) Doxygen will
 # generate Latex output.
 
-GENERATE_LATEX         = NO
+GENERATE_LATEX         = YES
 
 # The LATEX_OUTPUT tag is used to specify where the LaTeX docs will be put.
 # If a relative path is entered the value of OUTPUT_DIRECTORY will be

Installation/changes.html

@@ -625,6 +625,13 @@ and <code>src/</code> directories).
       <code>CGAL::write_ply_points()</code>
       and <code>CGAL::write_ply_points_and_normals()</code>.</li>
   </ul>
+  <h3>Point Set Shape Detection</h3>
+  <ul>
+    <li>New post-processing
+      algorithm: <code>CGAL::regularize_planes()</code>. This allows the user
+      to favor parallelism, orthogonality, coplanarity and/or axial
+      symmetry between detected planes.</li>
+  </ul>
   <h3>Surface Mesh Parameterization</h3>
   <ul>
     <li><code>LSCM_parameterizer_3</code> now uses by default Eigen

Point_set_processing_3/include/CGAL/IO/read_ply_points.h

@@ -27,6 +27,7 @@
 #include <CGAL/value_type_traits.h>
 #include <CGAL/point_set_processing_assertions.h>
 #include <CGAL/Kernel_traits.h>
+#include <CGAL/IO/io.h>
 
 #include <boost/version.hpp>
 #include <boost/cstdint.hpp>
@@ -86,6 +87,16 @@ namespace internal {
       c = static_cast<unsigned char>(s);
     }
     
+    void read_ascii (std::istream& stream, float& t) const
+    {
+      stream >> iformat(t);
+    }
+
+    void read_ascii (std::istream& stream, double& t) const
+    {
+      stream >> iformat(t);
+    }
+    
     // Default template when Type is not a char type
     template <typename Type>
     void read_ascii (std::istream& stream, Type& t) const
@@ -93,6 +104,7 @@ namespace internal {
       stream >> t;
     }
     
+    
     template <typename Type>
     Type read (std::istream& stream) const
     {

Polyhedron/demo/Polyhedron/CMakeLists.txt

@@ -249,6 +249,9 @@ if(CGAL_Qt5_FOUND AND Qt5_FOUND AND OPENGL_FOUND AND QGLVIEWER_FOUND)
   add_item(scene_polyhedron_selection_item Scene_polyhedron_selection_item.cpp)
   target_link_libraries(scene_polyhedron_selection_item scene_polyhedron_item_decorator scene_polyhedron_item_k_ring_selection)
 
+  add_item(scene_point_set_classification_item Scene_point_set_classification_item.cpp)
+  target_link_libraries(scene_point_set_classification_item scene_points_with_normal_item scene_color_ramp ${TBB_LIBRARIES})
+
   add_item(scene_polyhedron_shortest_path_item Plugins/Surface_mesh/Scene_polyhedron_shortest_path_item.cpp)
 
   add_item(scene_surface_mesh_item Scene_surface_mesh_item.cpp)
@@ -277,9 +280,6 @@ if(CGAL_Qt5_FOUND AND Qt5_FOUND AND OPENGL_FOUND AND QGLVIEWER_FOUND)
 
   target_link_libraries( demo_framework gl_splat)
 
-  
-
-
   foreach( lib 
       demo_framework
       scene_basic_objects

Polyhedron/demo/Polyhedron/Plugins/IO/CMakeLists.txt

@@ -60,3 +60,13 @@ target_link_libraries(xyz_plugin scene_points_with_normal_item)
 
 polyhedron_demo_plugin(ply_to_xyz_plugin PLY_to_xyz_io_plugin)
 target_link_libraries(ply_to_xyz_plugin scene_points_with_normal_item)
+
+find_package(libLAS)
+if (libLAS_FOUND)
+  include_directories(${libLAS_INCLUDE_DIRS})
+  
+  polyhedron_demo_plugin(las_to_xyz_plugin LAS_to_xyz_io_plugin)
+  target_link_libraries(las_to_xyz_plugin scene_points_with_normal_item ${libLAS_LIBRARIES})
+else()
+  message(STATUS "NOTICE : the las IO plugin needs LAS libraries and will not be compiled.")
+endif()

Polyhedron/demo/Polyhedron/Plugins/IO/LAS_to_xyz_io_plugin.cpp

@@ -0,0 +1,63 @@
+#include "Scene_points_with_normal_item.h"
+#include <CGAL/Three/Polyhedron_demo_io_plugin_interface.h>
+
+#include "read_las_point_set.h"
+
+#include <fstream>
+
+class Polyhedron_demo_las_to_xyz_plugin :
+  public QObject,
+  public CGAL::Three::Polyhedron_demo_io_plugin_interface
+{
+  Q_OBJECT
+  Q_INTERFACES(CGAL::Three::Polyhedron_demo_io_plugin_interface)
+  Q_PLUGIN_METADATA(IID "com.geometryfactory.PolyhedronDemo.IOPluginInterface/1.0")
+
+public:
+  QString name() const { return "las_to_xyz_plugin"; }
+  QString nameFilters() const { return "LAS files as Point set (*.las);;Compressed LAZ files as Point set (*.laz)"; }
+  bool canLoad() const;
+  CGAL::Three::Scene_item* load(QFileInfo fileinfo);
+
+  bool canSave(const CGAL::Three::Scene_item*);
+  bool save(const CGAL::Three::Scene_item*, QFileInfo fileinfo);
+};
+
+bool Polyhedron_demo_las_to_xyz_plugin::canLoad() const {
+  return true;
+}
+
+CGAL::Three::Scene_item*
+Polyhedron_demo_las_to_xyz_plugin::load(QFileInfo fileinfo) {
+  std::ifstream in(fileinfo.filePath().toUtf8(), std::ios_base::binary);
+
+  if(!in)
+    std::cerr << "Error!\n";
+
+  Scene_points_with_normal_item* item;
+  item = new Scene_points_with_normal_item();
+  if(!read_las_point_set(in, *(item->point_set())) || item->isEmpty())
+    {
+      delete item;
+      return 0;
+    }
+  item->point_set()->check_colors();
+  std::cerr << item->point_set()->info();
+  
+  item->setName(fileinfo.completeBaseName());
+  return item;
+}
+
+bool Polyhedron_demo_las_to_xyz_plugin::canSave(const CGAL::Three::Scene_item*)
+{
+  // No save function yet
+  return false;
+}
+
+bool Polyhedron_demo_las_to_xyz_plugin::save(const CGAL::Three::Scene_item*, QFileInfo)
+{
+  return false;
+}
+
+
+#include "LAS_to_xyz_io_plugin.moc"

Polyhedron/demo/Polyhedron/Plugins/PCA/CMakeLists.txt

@@ -11,3 +11,6 @@ target_link_libraries(edit_box_plugin  scene_edit_box_item  scene_polyhedron_ite
 
 polyhedron_demo_plugin(create_bbox_mesh_plugin Create_bbox_mesh_plugin)
 target_link_libraries(create_bbox_mesh_plugin scene_polyhedron_item)
+
+polyhedron_demo_plugin(recenter_plugin Recenter_plugin)
+target_link_libraries(recenter_plugin scene_points_with_normal_item)

Polyhedron/demo/Polyhedron/Plugins/PCA/Recenter_plugin.cpp

@@ -0,0 +1,81 @@
+#include "config.h"
+#include "Scene_points_with_normal_item.h"
+#include <CGAL/Three/Polyhedron_demo_plugin_helper.h>
+#include <CGAL/Three/Polyhedron_demo_plugin_interface.h>
+
+#include <CGAL/centroid.h>
+
+
+#include <QObject>
+#include <QAction>
+#include <QMainWindow>
+#include <QApplication>
+#include <QtPlugin>
+#include <QMessageBox>
+
+using namespace CGAL::Three;
+class Polyhedron_demo_recenter_plugin :
+  public QObject,
+  public Polyhedron_demo_plugin_helper
+{
+  Q_OBJECT
+  Q_INTERFACES(CGAL::Three::Polyhedron_demo_plugin_interface)
+  Q_PLUGIN_METADATA(IID "com.geometryfactory.PolyhedronDemo.PluginInterface/1.0")
+
+private:
+  QAction* actionRecenter;
+  
+public:
+  void init(QMainWindow* mainWindow, CGAL::Three::Scene_interface* scene_interface, Messages_interface*) {
+    scene = scene_interface;
+    mw = mainWindow;
+    actionRecenter = new QAction(tr("Recenter"), mainWindow);
+    actionRecenter->setObjectName("actionRecenter");
+    autoConnectActions();
+  }
+
+  QList<QAction*> actions() const {
+    return QList<QAction*>() << actionRecenter;
+  }
+  
+  //! Applicable if the currently selected item is a
+  //! points_with_normal_item (can be extended to other items).
+  bool applicable(QAction*) const {
+    return qobject_cast<Scene_points_with_normal_item*>(scene->item(scene->mainSelectionIndex()));
+  }
+
+public Q_SLOTS:
+  void on_actionRecenter_triggered();
+}; // end Polyhedron_demo_recenter_plugin
+
+void Polyhedron_demo_recenter_plugin::on_actionRecenter_triggered()
+{
+  const CGAL::Three::Scene_interface::Item_id index = scene->mainSelectionIndex();
+
+  Scene_points_with_normal_item* item =
+    qobject_cast<Scene_points_with_normal_item*>(scene->item(index));
+
+  if(item)
+  {
+    Point_set* points = item->point_set();
+    if(points == NULL)
+        return;
+
+    QApplication::setOverrideCursor(Qt::WaitCursor);
+
+    Kernel::Point_3 centroid = CGAL::centroid (points->points().begin(),
+                                               points->points().end());
+    Kernel::Vector_3 recenter (centroid, CGAL::ORIGIN);
+    std::cerr << "Recenter by " << recenter << std::endl;
+    for (Point_set::iterator it = points->begin(); it != points->end(); ++ it)
+      points->point(*it) = points->point(*it) + recenter;
+
+    QApplication::restoreOverrideCursor();
+    points->invalidate_bounds();
+    item->invalidateOpenGLBuffers();
+    scene->itemChanged(index);
+  }
+}
+
+
+#include "Recenter_plugin.moc"

Polyhedron/demo/Polyhedron/Plugins/Point_set/CMakeLists.txt

@@ -39,6 +39,10 @@ endif()
   polyhedron_demo_plugin(point_set_selection_plugin Point_set_selection_plugin ${point_set_selectionUI_FILES})
   target_link_libraries(point_set_selection_plugin scene_points_with_normal_item scene_polylines_item scene_edit_box_item)
 
+  qt5_wrap_ui( point_set_classificationUI_FILES Point_set_classification_widget.ui)
+  polyhedron_demo_plugin(point_set_classification_plugin Point_set_classification_plugin ${point_set_classificationUI_FILES})
+  target_link_libraries(point_set_classification_plugin scene_point_set_classification_item scene_points_with_normal_item scene_polylines_item scene_polygon_soup_item)
+
   qt5_wrap_ui(point_set_shape_detectionUI_FILES Point_set_shape_detection_plugin.ui)
   polyhedron_demo_plugin(point_set_shape_detection_plugin Point_set_shape_detection_plugin ${point_set_shape_detectionUI_FILES})
   target_link_libraries(point_set_shape_detection_plugin scene_polyhedron_item scene_points_with_normal_item scene_polygon_soup_item)

Polyhedron/demo/Polyhedron/Plugins/Point_set/Point_set_classification_plugin.cpp

@@ -0,0 +1,903 @@
+#include <QtCore/qglobal.h>
+#include <QFileDialog>
+#include <QColorDialog> 
+#include <fstream>
+#include "opengl_tools.h"
+
+#include "Messages_interface.h"
+#include "Scene_points_with_normal_item.h"
+#include "Scene_point_set_classification_item.h"
+#include "Scene_polylines_item.h"
+#include "Scene_polygon_soup_item.h"
+
+#include <CGAL/Three/Scene_interface.h>
+#include <CGAL/Three/Polyhedron_demo_plugin_helper.h>
+
+#include <CGAL/Random.h>
+
+#include "ui_Point_set_classification_widget.h"
+
+#include <QAction>
+#include <QMainWindow>
+#include <QApplication>
+#include <QCheckBox>
+#include <QInputDialog>
+#include <QSlider>
+
+#include <map>
+
+#include <boost/graph/adjacency_list.hpp>
+#include <CGAL/boost/graph/split_graph_into_polylines.h>
+
+using namespace CGAL::Three;
+
+class Polyhedron_demo_point_set_classification_plugin :
+  public QObject,
+  public Polyhedron_demo_plugin_helper
+{
+  Q_OBJECT
+    Q_INTERFACES(CGAL::Three::Polyhedron_demo_plugin_interface)
+    Q_PLUGIN_METADATA(IID "com.geometryfactory.PolyhedronDemo.PluginInterface/1.0")
+    
+  struct ClassRow
+  {
+    QLabel* label;
+
+    QPushButton* color_button;
+    QPushButton* train;
+    QPushButton* remove;
+    QColor color;
+
+    QLabel* label2;
+    QComboBox* effect;
+
+    ClassRow (QWidget* parent, const char* name, const QColor& color)
+      : color (color)
+    {
+      label = new QLabel (name, parent);
+      color_button = new QPushButton ("", parent);
+      
+      QString s("background: #"
+                + QString(color.red() < 16? "0" : "") + QString::number(color.red(),16)
+                + QString(color.green() < 16? "0" : "") + QString::number(color.green(),16)
+                + QString(color.blue() < 16? "0" : "") + QString::number(color.blue(),16) + ";");
+      color_button->setStyleSheet(s);
+      
+      train = new QPushButton ("Add selection");
+      remove = new QPushButton ("Remove");
+      label2 = new QLabel (name, parent);
+      effect = new QComboBox;
+      effect->addItem("Penalized");
+      effect->addItem("Neutral");
+      effect->addItem("Favored");
+
+    }
+    ~ClassRow ()
+    {
+    }
+    void change_color (const QColor& color)
+    {
+      this->color = color;
+      QString s("background: #"
+                + QString(color.red() < 16? "0" : "") + QString::number(color.red(),16)
+                + QString(color.green() < 16? "0" : "") + QString::number(color.green(),16)
+                + QString(color.blue() < 16? "0" : "") + QString::number(color.blue(),16) + ";");
+      color_button->setStyleSheet(s);
+      color_button->update();
+    }
+  };
+
+
+public:
+  bool applicable(QAction*) const { 
+      return
+        qobject_cast<Scene_points_with_normal_item*>(scene->item(scene->mainSelectionIndex()));
+  }
+  void print_message(QString message) { messages->information(message); }
+  QList<QAction*> actions() const { return QList<QAction*>() << actionPointSetClassification; }
+  
+  using Polyhedron_demo_plugin_helper::init;
+  void init(QMainWindow* mainWindow, CGAL::Three::Scene_interface* scene_interface, Messages_interface* m) {
+    mw = mainWindow;
+    scene = scene_interface;
+    messages = m;
+    actionPointSetClassification = new QAction(tr("Point Set Classification"), mw);
+    connect(actionPointSetClassification, SIGNAL(triggered()), this, SLOT(point_set_classification_action()));
+
+    dock_widget = new QDockWidget("Point Set Classification", mw);
+    dock_widget->setVisible(false);
+
+    ui_widget.setupUi(dock_widget);
+    addDockWidget(dock_widget);
+
+    color_att = QColor (75, 75, 77);
+    
+    connect(ui_widget.create_from_item,  SIGNAL(clicked()), this,
+            SLOT(on_create_from_item_button_clicked()));
+    connect(ui_widget.compute_features,  SIGNAL(clicked()), this,
+            SLOT(on_compute_features_button_clicked()));
+    connect(ui_widget.display,  SIGNAL(currentIndexChanged(int)), this,
+            SLOT(on_display_button_clicked(int)));
+    connect(ui_widget.run,  SIGNAL(clicked()), this,
+            SLOT(on_run_button_clicked()));
+    connect(ui_widget.save_config,  SIGNAL(clicked()), this,
+            SLOT(on_save_config_button_clicked()));
+    connect(ui_widget.load_config,  SIGNAL(clicked()), this,
+            SLOT(on_load_config_button_clicked()));
+    connect(ui_widget.run_smoothed,  SIGNAL(clicked()), this,
+            SLOT(on_run_smoothed_button_clicked()));
+    connect(ui_widget.run_graphcut,  SIGNAL(clicked()), this,
+            SLOT(on_run_graphcut_button_clicked()));
+    connect(ui_widget.smoothingDoubleSpinBox,  SIGNAL(valueChanged(double)), this,
+            SLOT(on_smoothing_value_changed(double)));
+    connect(ui_widget.subdivisionsSpinBox,  SIGNAL(valueChanged(int)), this,
+            SLOT(on_subdivisions_value_changed(int)));
+    connect(ui_widget.save,  SIGNAL(clicked()), this,
+            SLOT(on_save_button_clicked()));
+    connect(ui_widget.generate_point_set_items,  SIGNAL(clicked()), this,
+            SLOT(on_generate_point_set_items_button_clicked()));
+
+    connect(ui_widget.numberOfScalesSpinBox,  SIGNAL(valueChanged(int)), this,
+            SLOT(on_update_nb_scales()));
+    connect(ui_widget.number_of_trials,  SIGNAL(valueChanged(int)), this,
+            SLOT(on_update_number_of_trials()));
+
+    connect(ui_widget.add_new_label,  SIGNAL(clicked()), this,
+            SLOT(on_add_new_label_clicked()));
+    connect(ui_widget.reset_training_sets,  SIGNAL(clicked()), this,
+            SLOT(on_reset_training_sets_clicked()));
+    connect(ui_widget.validate_selection,  SIGNAL(clicked()), this,
+            SLOT(on_validate_selection_clicked()));
+    connect(ui_widget.train,  SIGNAL(clicked()), this,
+            SLOT(on_train_clicked()));
+
+    connect(ui_widget.selected_feature,  SIGNAL(currentIndexChanged(int)), this,
+            SLOT(on_selected_feature_changed(int)));
+    connect(ui_widget.feature_weight,  SIGNAL(valueChanged(int)), this,
+            SLOT(on_feature_weight_changed(int)));
+
+    QObject* scene_obj = dynamic_cast<QObject*>(scene_interface);
+    if(scene_obj)
+      {
+        connect(scene_obj, SIGNAL(itemAboutToBeDestroyed(CGAL::Three::Scene_item*)), this,
+                SLOT(item_about_to_be_destroyed(CGAL::Three::Scene_item*)));
+        
+        connect(scene_obj, SIGNAL(itemIndexSelected(int)), this,
+                SLOT(update_plugin(int)));
+      }
+  }
+  virtual void closure()
+  {
+    dock_widget->hide();
+  }
+
+
+public Q_SLOTS:
+  
+  void item_about_to_be_destroyed(CGAL::Three::Scene_item* scene_item) {
+    // if points item
+    Scene_points_with_normal_item* points_item = qobject_cast<Scene_points_with_normal_item*>(scene_item);
+    if(points_item)
+      {
+        Item_map::iterator it = item_map.find(points_item);
+        if (it != item_map.end())
+          {
+            Scene_point_set_classification_item* classif_item = it->second;
+            item_map.erase(it); // first erase from map, because scene->erase will cause a call to this function
+            scene->erase( scene->item_id(classif_item) );
+          }
+      }
+
+    // if classification item
+    Scene_point_set_classification_item* classif_item = qobject_cast<Scene_point_set_classification_item*>(scene_item);
+    if(classif_item)
+      item_map.erase (classif_item->points_item());
+  }
+
+  void point_set_classification_action()
+  { 
+    dock_widget->show();
+    dock_widget->raise();
+    Scene_points_with_normal_item* points_item = getSelectedItem<Scene_points_with_normal_item>();
+
+    if (item_map.find(points_item) == item_map.end())
+      on_create_from_item_button_clicked();
+  }
+
+  void update_plugin(int)
+  {
+    update_plugin_from_item(get_classification_item());
+  }
+
+  void disable_everything ()
+  {
+    ui_widget.load_config->setEnabled(false);
+    ui_widget.save_config->setEnabled(false);
+    ui_widget.compute_features->setEnabled(false);
+    ui_widget.numberOfScalesSpinBox->setEnabled(false);
+    ui_widget.display->setEnabled(false);
+    ui_widget.tabWidget->setEnabled(false);
+    ui_widget.run->setEnabled(false);
+    ui_widget.run_smoothed->setEnabled(false);
+    ui_widget.frame->setEnabled(false);
+    ui_widget.save->setEnabled(false);
+    ui_widget.generate_point_set_items->setEnabled(false);
+  }
+
+  void enable_computation()
+  {
+    ui_widget.load_config->setEnabled(true);
+    ui_widget.compute_features->setEnabled(true);
+    ui_widget.numberOfScalesSpinBox->setEnabled(true);
+    ui_widget.display->setEnabled(true);
+  }
+
+  void enable_classif()
+  {
+    ui_widget.save_config->setEnabled(true);
+    ui_widget.tabWidget->setEnabled(true);
+    ui_widget.run->setEnabled(true);
+    ui_widget.run_smoothed->setEnabled(true);
+    ui_widget.frame->setEnabled(true);
+    ui_widget.save->setEnabled(true);
+    ui_widget.generate_point_set_items->setEnabled(true);
+  }
+
+
+  void update_plugin_from_item(Scene_point_set_classification_item* item)
+  {
+    if (item == NULL) // Deactivate plugin
+      {
+        disable_everything();
+        ui_widget.tabWidget->setCurrentIndex(0);
+      }
+    else
+      {
+        disable_everything();
+        enable_computation();
+        
+        ui_widget.numberOfScalesSpinBox->setValue((int)(item->nb_scales()));
+        ui_widget.number_of_trials->setValue((int)(item->number_of_trials()));
+        ui_widget.smoothingDoubleSpinBox->setValue((int)(item->smoothing()));
+
+        // Clear class labels
+        for (std::size_t i = 0; i < class_rows.size(); ++ i)
+          {
+
+            ui_widget.gridLayout_3->removeWidget (class_rows[i].label);
+            delete class_rows[i].label;
+            ui_widget.gridLayout_3->removeWidget (class_rows[i].color_button);
+            delete class_rows[i].color_button;
+            ui_widget.gridLayout_3->removeWidget (class_rows[i].train);
+            delete class_rows[i].train;
+            ui_widget.gridLayout_3->removeWidget (class_rows[i].remove);
+            delete class_rows[i].remove;
+            ui_widget.gridLayout->removeWidget (class_rows[i].label2);
+            delete class_rows[i].label2;
+            ui_widget.gridLayout->removeWidget (class_rows[i].effect);
+            delete class_rows[i].effect;
+          }
+        class_rows.clear();
+
+        // Add labels
+        for (std::size_t i = 0; i < item->labels().size(); ++ i)
+          add_new_label (ClassRow (dock_widget, item->labels()[i].first->name().c_str(),
+                                   item->labels()[i].second));
+
+        // Enabled classif if features computed
+        if (!(item->features_computed()))
+          ui_widget.tabWidget->setCurrentIndex(0);
+        else
+          enable_classif();
+
+        int index = ui_widget.display->currentIndex();
+        ui_widget.display->clear();
+        ui_widget.display->addItem("Real colors");
+        ui_widget.display->addItem("Classification");
+        ui_widget.display->addItem("Training sets");
+        ui_widget.selected_feature->clear();
+        item->fill_display_combo_box(ui_widget.display, ui_widget.selected_feature);
+        if (index >= ui_widget.display->count())
+          ui_widget.display->setCurrentIndex(1);
+        else
+          ui_widget.display->setCurrentIndex(index);
+        ui_widget.selected_feature->setCurrentIndex(0);
+      }
+  }
+
+  void on_update_nb_scales()
+  {
+    Scene_point_set_classification_item* classification_item
+      = get_classification_item();
+    if(!classification_item)
+      return; 
+    classification_item->nb_scales() = ui_widget.numberOfScalesSpinBox->value();
+  }
+  void on_update_number_of_trials()
+  {
+    Scene_point_set_classification_item* classification_item
+      = get_classification_item();
+    if(!classification_item)
+      return; 
+    classification_item->number_of_trials() = ui_widget.number_of_trials->value();
+  }
+                           
+  
+  Scene_point_set_classification_item* get_classification_item()
+  {
+    Scene_point_set_classification_item* out =
+      qobject_cast<Scene_point_set_classification_item*>(scene->item(scene->mainSelectionIndex()));
+    if (out)
+      return out;
+
+    Scene_points_with_normal_item* points_item =
+      qobject_cast<Scene_points_with_normal_item*>(scene->item(scene->mainSelectionIndex()));
+    if (!points_item)
+      return NULL;
+    
+    Item_map::iterator it = item_map.find(points_item);
+    if (it != item_map.end())
+      return it->second;
+
+    return NULL;
+  }
+  
+
+  void on_create_from_item_button_clicked()
+  {
+    Scene_points_with_normal_item* points_item = getSelectedItem<Scene_points_with_normal_item>();
+    if(!points_item)
+      {
+        print_message("Error: there is no selected point set item!");
+        return; 
+      }
+    if (item_map.find(points_item) != item_map.end())
+      return;
+
+    QApplication::setOverrideCursor(Qt::WaitCursor);
+    Scene_point_set_classification_item* new_item
+      = new Scene_point_set_classification_item (points_item);
+    scene->addItem(new_item);
+    new_item->setName(QString("%1 (classification)").arg(points_item->name()));
+    points_item->setVisible (false);
+    item_map.insert (std::make_pair (points_item, new_item));
+    QApplication::restoreOverrideCursor();
+    connect(points_item, &Scene_points_with_normal_item::itemChanged,
+            new_item, &Scene_point_set_classification_item::invalidateOpenGLBuffers);
+
+    update_plugin_from_item(new_item);
+  }
+
+  void run (Scene_point_set_classification_item* classification_item, int method)
+  {
+    classification_item->run (method);
+  }
+
+  void on_compute_features_button_clicked()
+  {
+    Scene_point_set_classification_item* classification_item
+      = get_classification_item();
+    if(!classification_item)
+      {
+        print_message("Error: there is no point set classification item!");
+        return; 
+      }
+    
+    QApplication::setOverrideCursor(Qt::WaitCursor);
+
+    classification_item->compute_features ();
+
+    update_plugin_from_item(classification_item);
+    QApplication::restoreOverrideCursor();
+    scene->itemChanged(classification_item);
+  }
+
+  void on_save_config_button_clicked()
+  {
+    Scene_point_set_classification_item* classification_item
+      = get_classification_item();
+    if(!classification_item)
+      {
+        print_message("Error: there is no point set classification item!");
+        return; 
+      }
+
+    QString filename = QFileDialog::getSaveFileName(mw,
+                                                    tr("Save classification configuration"),
+                                                    QString("config.xml"),
+                                                    "Config file (*.xml);;");
+    if (filename == QString())
+      return;
+
+    
+    QApplication::setOverrideCursor(Qt::WaitCursor);
+
+    classification_item->save_config (filename.toStdString().c_str());
+    
+    QApplication::restoreOverrideCursor();
+
+  }
+
+  void on_load_config_button_clicked()
+  {
+    Scene_point_set_classification_item* classification_item
+      = get_classification_item();
+    if(!classification_item)
+      {
+        print_message("Error: there is no point set classification item!");
+        return; 
+      }
+    QString filename = QFileDialog::getOpenFileName(mw,
+                                                    tr("Open classification configuration"),
+                                                    ".",
+                                                    "Config file (*.xml);;All Files (*)");
+
+    if (filename == QString())
+      return;
+
+    QApplication::setOverrideCursor(Qt::WaitCursor);
+
+    classification_item->load_config (filename.toStdString().c_str());
+    update_plugin_from_item(classification_item);
+    run (classification_item, 0);
+    
+    QApplication::restoreOverrideCursor();
+    scene->itemChanged(classification_item);
+  }
+
+
+  void on_display_button_clicked(int index)
+  {
+    Scene_point_set_classification_item* classification_item
+      = get_classification_item();
+    if(!classification_item)
+      return; 
+
+    classification_item->change_color (index);
+    scene->itemChanged(classification_item);
+  }
+
+  void on_run_button_clicked()
+  {
+    Scene_point_set_classification_item* classification_item
+      = get_classification_item();
+    if(!classification_item)
+      {
+        print_message("Error: there is no point set classification item!");
+        return; 
+      }
+    QApplication::setOverrideCursor(Qt::WaitCursor);
+    run (classification_item, 0);
+    QApplication::restoreOverrideCursor();
+    scene->itemChanged(classification_item);
+  }
+
+  void on_run_smoothed_button_clicked()
+  {
+    Scene_point_set_classification_item* classification_item
+      = get_classification_item();
+    if(!classification_item)
+      {
+        print_message("Error: there is no point set classification item!");
+        return; 
+      }
+
+    QApplication::setOverrideCursor(Qt::WaitCursor);
+    QTime time;
+    time.start();
+    run (classification_item, 1);
+    std::cerr << "Smoothed classification computed in " << time.elapsed() / 1000 << " second(s)" << std::endl;
+    QApplication::restoreOverrideCursor();
+    scene->itemChanged(classification_item);
+  }
+  
+  void on_run_graphcut_button_clicked()
+  {
+    Scene_point_set_classification_item* classification_item
+      = get_classification_item();
+    if(!classification_item)
+      {
+        print_message("Error: there is no point set classification item!");
+        return; 
+      }
+
+    QApplication::setOverrideCursor(Qt::WaitCursor);
+    QTime time;
+    time.start();
+    run (classification_item, 2);
+    std::cerr << "Graphcut classification computed in " << time.elapsed() / 1000 << " second(s)" << std::endl;
+    QApplication::restoreOverrideCursor();
+    scene->itemChanged(classification_item);
+  }
+
+  void on_smoothing_value_changed(double v)
+  {
+    Scene_point_set_classification_item* classification_item
+      = get_classification_item();
+    if(!classification_item)
+      return; 
+    classification_item->smoothing() = v;
+  }
+
+  void on_subdivisions_value_changed(int v)
+  {
+    Scene_point_set_classification_item* classification_item
+      = get_classification_item();
+    if(!classification_item)
+      return; 
+    classification_item->subdivisions() = v;
+  }
+
+ void on_save_button_clicked()
+  {
+    Scene_point_set_classification_item* classification_item
+      = get_classification_item();
+    if(!classification_item)
+      {
+        print_message("Error: there is no point set classification item!");
+        return; 
+      }
+
+    QString filename = QFileDialog::getSaveFileName(mw,
+                                                    tr("Save PLY classified point set"),
+                                                    QString("%1.ply").arg(classification_item->name()),
+                                                    "PLY point set (*.ply);;");
+
+    if (filename == QString())
+      return;
+
+    std::ofstream out(filename.toUtf8());
+    
+    QApplication::setOverrideCursor(Qt::WaitCursor);
+    classification_item->write_ply_point_set (out);
+    QApplication::restoreOverrideCursor();
+
+    out.close();
+  }
+
+  void on_generate_point_set_items_button_clicked()
+  {
+    Scene_point_set_classification_item* classification_item
+      = get_classification_item();
+    if(!classification_item)
+      {
+        print_message("Error: there is no point set classification item!");
+        return; 
+      }
+
+    
+    QApplication::setOverrideCursor(Qt::WaitCursor);
+
+    std::vector<Scene_points_with_normal_item*> new_items;
+    classification_item->generate_point_set_items<Scene_points_with_normal_item>
+      (new_items, classification_item->points_item()->name().toStdString().c_str());
+
+    for (std::size_t i = 0; i < new_items.size(); ++ i)
+      {
+        if (new_items[i]->point_set()->empty())
+          delete new_items[i];
+        else
+          scene->addItem (new_items[i]);
+      }
+    
+    QApplication::restoreOverrideCursor();
+  }
+
+  void on_add_new_label_clicked()
+  {
+    Scene_point_set_classification_item* classification_item
+      = get_classification_item();
+    if(!classification_item)
+      {
+        print_message("Error: there is no point set classification item!");
+        return; 
+      }
+
+    bool ok;
+    QString name =
+      QInputDialog::getText((QWidget*)mw,
+                            tr("Add new label"), // dialog title
+                            tr("Name:"), // field label
+                            QLineEdit::Normal,
+                            tr("label%1").arg(class_rows.size() + 1),
+                            &ok);
+    if (!ok)
+      return;
+
+    add_new_label (ClassRow (dock_widget, name.toStdString().c_str(),
+                                      QColor (192 + rand() % 60,
+                                              192 + rand() % 60,
+                                              192 + rand() % 60)));
+    classification_item->add_new_label (class_rows.back().label->text().toStdString().c_str(),
+                                        class_rows.back().color);
+    
+
+  }
+
+  void on_reset_training_sets_clicked()
+  {
+    Scene_point_set_classification_item* classification_item
+      = get_classification_item();
+    if(!classification_item)
+      {
+        print_message("Error: there is no point set classification item!");
+        return; 
+      }
+
+    classification_item->reset_training_sets();
+  }
+
+  void on_validate_selection_clicked()
+  {
+    Scene_point_set_classification_item* classification_item
+      = get_classification_item();
+    if(!classification_item)
+      {
+        print_message("Error: there is no point set classification item!");
+        return; 
+      }
+
+    classification_item->validate_selection();
+  }
+
+  void on_train_clicked()
+  {
+    Scene_point_set_classification_item* classification_item
+      = get_classification_item();
+    if(!classification_item)
+      {
+        print_message("Error: there is no point set classification item!");
+        return; 
+      }
+
+    std::vector<std::string> classes;
+    std::vector<QColor> colors;
+    
+    for (std::size_t i = 0; i < class_rows.size(); ++ i)
+      {
+        classes.push_back (class_rows[i].label->text().toStdString());
+        colors.push_back (class_rows[i].color);
+      }
+    QApplication::setOverrideCursor(Qt::WaitCursor);
+    classification_item->train();
+    QApplication::restoreOverrideCursor();
+    update_plugin_from_item(classification_item);
+  }
+
+  void add_new_label (const ClassRow& class_row)
+  {
+    class_rows.push_back (class_row);
+    int position = static_cast<int>(class_rows.size());
+
+    ui_widget.gridLayout_3->addWidget (class_rows.back().label, position, 0);
+    ui_widget.gridLayout_3->addWidget (class_rows.back().color_button, position, 1);
+    ui_widget.gridLayout_3->addWidget (class_rows.back().train, position, 3);
+    ui_widget.gridLayout_3->addWidget (class_rows.back().remove, position, 5);
+
+    connect(class_rows.back().remove,  SIGNAL(clicked()), this,
+            SLOT(on_remove_class_clicked()));
+    connect(class_rows.back().color_button,  SIGNAL(clicked()), this,
+            SLOT(on_color_changed_clicked()));
+    connect(class_rows.back().train,  SIGNAL(clicked()), this,
+            SLOT(on_add_selection_to_training_set_clicked()));
+
+    ui_widget.gridLayout->addWidget (class_rows.back().label2, position, 0);
+    ui_widget.gridLayout->addWidget (class_rows.back().effect, position, 2);
+
+    connect(class_rows.back().effect,  SIGNAL(currentIndexChanged(int)), this,
+            SLOT(on_effect_changed(int)));
+
+  }
+
+  void on_remove_class_clicked()
+  {
+    Scene_point_set_classification_item* classification_item
+      = get_classification_item();
+    if(!classification_item)
+      {
+        print_message("Error: there is no point set classification item!");
+        return; 
+      }
+
+    int index = ui_widget.gridLayout_3->indexOf(qobject_cast<QWidget*>(QObject::sender()));
+    
+    int row_index, column_index, row_span, column_span;
+    ui_widget.gridLayout_3->getItemPosition(index, &row_index, &column_index, &row_span, &column_span);
+    --row_index;
+
+    classification_item->remove_label (class_rows[row_index].label->text().toStdString().c_str());
+
+    ui_widget.gridLayout_3->removeWidget (class_rows[row_index].label);
+    delete class_rows[row_index].label;
+    ui_widget.gridLayout_3->removeWidget (class_rows[row_index].color_button);
+    delete class_rows[row_index].color_button;
+    ui_widget.gridLayout_3->removeWidget (class_rows[row_index].train);
+    delete class_rows[row_index].train;
+    ui_widget.gridLayout_3->removeWidget (class_rows[row_index].remove);
+    delete class_rows[row_index].remove;
+
+    ui_widget.gridLayout->removeWidget (class_rows[row_index].label2);
+    delete class_rows[row_index].label2;
+    ui_widget.gridLayout->removeWidget (class_rows[row_index].effect);
+    delete class_rows[row_index].effect;
+
+    if (class_rows.size() > 1)
+      for (std::size_t i = row_index + 1; i < class_rows.size(); ++ i)
+        {
+          ui_widget.gridLayout_3->addWidget (class_rows[i].label, (int)i, 0);
+          ui_widget.gridLayout_3->addWidget (class_rows[i].color_button, (int)i, 1);
+          ui_widget.gridLayout_3->addWidget (class_rows[i].train, (int)i, 3);
+          ui_widget.gridLayout_3->addWidget (class_rows[i].remove, (int)i, 5);
+          ui_widget.gridLayout->addWidget (class_rows[i].label2, (int)i, 0);
+          ui_widget.gridLayout->addWidget (class_rows[i].effect, (int)i, 2);
+        }
+
+    class_rows.erase (class_rows.begin() + row_index);
+    
+    scene->itemChanged(classification_item);
+  }
+
+  void on_color_changed_clicked()
+  {
+    Scene_point_set_classification_item* classification_item
+      = get_classification_item();
+    if(!classification_item)
+      {
+        print_message("Error: there is no point set classification item!");
+        return; 
+      }
+
+    QPushButton* color_button = qobject_cast<QPushButton*>(QObject::sender());
+    int index = ui_widget.gridLayout_3->indexOf(color_button);
+    int row_index, column_index, row_span, column_span;
+    ui_widget.gridLayout_3->getItemPosition(index, &row_index, &column_index, &row_span, &column_span);
+    -- row_index;
+    
+    QColor color = class_rows[row_index].color;
+    color = QColorDialog::getColor(color, (QWidget*)mw, "Change of color of label");
+    class_rows[row_index].change_color (color);
+    classification_item->change_label_color (class_rows[row_index].label->text().toStdString().c_str(),
+                                            color);
+
+    scene->itemChanged(classification_item);
+  }
+
+  void on_add_selection_to_training_set_clicked()
+  {
+    Scene_point_set_classification_item* classification_item
+      = get_classification_item();
+    if(!classification_item)
+      {
+        print_message("Error: there is no point set classification item!");
+        return; 
+      }
+
+    int index = ui_widget.gridLayout_3->indexOf(qobject_cast<QWidget*>(QObject::sender()));
+    int row_index, column_index, row_span, column_span;
+    ui_widget.gridLayout_3->getItemPosition(index, &row_index, &column_index, &row_span, &column_span);
+    --row_index;
+
+    classification_item->add_selection_to_training_set
+      (class_rows[row_index].label->text().toStdString().c_str());
+  }
+
+  void on_selected_feature_changed(int v)
+  {
+    Scene_point_set_classification_item* classification_item
+      = get_classification_item();
+    if(!classification_item)
+      {
+        print_message("Error: there is no point set classification item!");
+        return; 
+      }
+
+    if (classification_item->number_of_features() <= (std::size_t)v)
+      return;
+    
+    Scene_point_set_classification_item::Feature_handle
+      att = classification_item->feature(v);
+
+    if (att == Scene_point_set_classification_item::Feature_handle())
+      return;
+
+    // std::cerr << att->weight()
+    //           << " " << (int)(1001. * 2. * std::atan(att->weight()) / CGAL_PI) << std::endl;
+    ui_widget.feature_weight->setValue ((int)(1001. * 2. * std::atan(att->weight()) / CGAL_PI));
+
+    for (std::size_t i = 0; i < classification_item->labels().size(); ++ i)
+      {
+        CGAL::Classification::Feature::Effect
+          eff = classification_item->labels()[i].first->feature_effect(att);
+        if (eff == CGAL::Classification::Feature::PENALIZING)
+          class_rows[i].effect->setCurrentIndex(0);
+        else if (eff == CGAL::Classification::Feature::NEUTRAL)
+          class_rows[i].effect->setCurrentIndex(1);
+        else
+          class_rows[i].effect->setCurrentIndex(2);
+      }
+  }
+
+
+  void on_feature_weight_changed(int v)
+  {
+    Scene_point_set_classification_item* classification_item
+      = get_classification_item();
+    if(!classification_item)
+      {
+        print_message("Error: there is no point set classification item!");
+        return; 
+      }
+    Scene_point_set_classification_item::Feature_handle
+      att = classification_item->feature(ui_widget.selected_feature->currentIndex());
+
+    if (att == Scene_point_set_classification_item::Feature_handle())
+      return;
+
+    att->set_weight(std::tan ((CGAL_PI/2.) * v / 1001.));
+    // std::cerr << att->weight() << std::endl;
+
+    for (std::size_t i = 0; i < class_rows.size(); ++ i)
+      class_rows[i].effect->setEnabled(att->weight() != 0.);
+  }
+
+  void on_effect_changed (int v)
+  {
+    Scene_point_set_classification_item* classification_item
+      = get_classification_item();
+    if(!classification_item)
+      {
+        print_message("Error: there is no point set classification item!");
+        return; 
+      }
+    Scene_point_set_classification_item::Feature_handle
+      att = classification_item->feature(ui_widget.selected_feature->currentIndex());
+
+    if (att == Scene_point_set_classification_item::Feature_handle())
+      return;
+
+    QComboBox* combo = qobject_cast<QComboBox*>(QObject::sender());
+    for (std::size_t i = 0;i < class_rows.size(); ++ i)
+      if (class_rows[i].effect == combo)
+        {
+          //          std::cerr << att->id() << " is ";
+          if (v == 0)
+            {
+              classification_item->labels()[i].first->set_feature_effect
+                (att, CGAL::Classification::Feature::PENALIZING);
+              //              std::cerr << " penalized for ";
+            }
+          else if (v == 1)
+            {
+              classification_item->labels()[i].first->set_feature_effect
+                (att, CGAL::Classification::Feature::NEUTRAL);
+              //              std::cerr << " neutral for ";
+            }
+          else
+            {
+              classification_item->labels()[i].first->set_feature_effect
+                (att, CGAL::Classification::Feature::FAVORING);
+              //              std::cerr << " favored for ";
+            }
+          //          std::cerr << classification_item->labels()[i].first->id() << std::endl;
+          break;
+        }
+  }
+
+private:
+  Messages_interface* messages;
+  QAction* actionPointSetClassification;
+
+  QDockWidget* dock_widget;
+
+  std::vector<ClassRow> class_rows;
+  
+  Ui::PointSetClassification ui_widget;
+
+  QColor color_att;
+
+  typedef std::map<Scene_points_with_normal_item*, Scene_point_set_classification_item*> Item_map;
+  Item_map item_map;
+
+}; // end Polyhedron_demo_point_set_classification_plugin
+
+#include "Point_set_classification_plugin.moc"

Polyhedron/demo/Polyhedron/Plugins/Point_set/Point_set_classification_widget.ui

@@ -0,0 +1,757 @@
+<?xml version="1.0" encoding="UTF-8"?>
+<ui version="4.0">
+ <class>PointSetClassification</class>
+ <widget class="QDockWidget" name="PointSetClassification">
+  <property name="geometry">
+   <rect>
+    <x>0</x>
+    <y>0</y>
+    <width>414</width>
+    <height>670</height>
+   </rect>
+  </property>
+  <property name="windowTitle">
+   <string>Point set classification</string>
+  </property>
+  <widget class="QWidget" name="dockWidgetContents">
+   <layout class="QVBoxLayout" name="verticalLayout_4">
+    <item>
+     <spacer name="verticalSpacer_12">
+      <property name="orientation">
+       <enum>Qt::Vertical</enum>
+      </property>
+      <property name="sizeHint" stdset="0">
+       <size>
+        <width>20</width>
+        <height>40</height>
+       </size>
+      </property>
+     </spacer>
+    </item>
+    <item>
+     <layout class="QHBoxLayout" name="horizontalLayout_8">
+      <item>
+       <widget class="QPushButton" name="create_from_item">
+        <property name="font">
+         <font>
+          <weight>75</weight>
+          <bold>true</bold>
+         </font>
+        </property>
+        <property name="text">
+         <string>Create classification item</string>
+        </property>
+       </widget>
+      </item>
+      <item>
+       <widget class="QPushButton" name="load_config">
+        <property name="text">
+         <string>Load config</string>
+        </property>
+       </widget>
+      </item>
+      <item>
+       <widget class="QPushButton" name="save_config">
+        <property name="text">
+         <string>Save config</string>
+        </property>
+       </widget>
+      </item>
+     </layout>
+    </item>
+    <item>
+     <spacer name="verticalSpacer">
+      <property name="orientation">
+       <enum>Qt::Vertical</enum>
+      </property>
+      <property name="sizeHint" stdset="0">
+       <size>
+        <width>20</width>
+        <height>40</height>
+       </size>
+      </property>
+     </spacer>
+    </item>
+    <item>
+     <layout class="QHBoxLayout" name="horizontalLayout_7">
+      <item>
+       <widget class="QPushButton" name="compute_features">
+        <property name="font">
+         <font>
+          <weight>75</weight>
+          <bold>true</bold>
+         </font>
+        </property>
+        <property name="text">
+         <string>Compute features</string>
+        </property>
+       </widget>
+      </item>
+      <item>
+       <widget class="QSpinBox" name="numberOfScalesSpinBox">
+        <property name="suffix">
+         <string> scale(s)</string>
+        </property>
+        <property name="prefix">
+         <string>using </string>
+        </property>
+        <property name="minimum">
+         <number>1</number>
+        </property>
+        <property name="value">
+         <number>5</number>
+        </property>
+       </widget>
+      </item>
+     </layout>
+    </item>
+    <item>
+     <spacer name="verticalSpacer_3">
+      <property name="orientation">
+       <enum>Qt::Vertical</enum>
+      </property>
+      <property name="sizeHint" stdset="0">
+       <size>
+        <width>20</width>
+        <height>40</height>
+       </size>
+      </property>
+     </spacer>
+    </item>
+    <item>
+     <layout class="QHBoxLayout" name="horizontalLayout">
+      <item>
+       <widget class="QLabel" name="label">
+        <property name="sizePolicy">
+         <sizepolicy hsizetype="Maximum" vsizetype="Preferred">
+          <horstretch>0</horstretch>
+          <verstretch>0</verstretch>
+         </sizepolicy>
+        </property>
+        <property name="text">
+         <string>View:</string>
+        </property>
+       </widget>
+      </item>
+      <item>
+       <widget class="QComboBox" name="display">
+        <property name="currentIndex">
+         <number>0</number>
+        </property>
+        <item>
+         <property name="text">
+          <string>Real colors</string>
+         </property>
+        </item>
+        <item>
+         <property name="text">
+          <string>Classification</string>
+         </property>
+        </item>
+        <item>
+         <property name="text">
+          <string>Training sets</string>
+         </property>
+        </item>
+       </widget>
+      </item>
+     </layout>
+    </item>
+    <item>
+     <widget class="QTabWidget" name="tabWidget">
+      <property name="contextMenuPolicy">
+       <enum>Qt::NoContextMenu</enum>
+      </property>
+      <property name="currentIndex">
+       <number>0</number>
+      </property>
+      <widget class="QWidget" name="tab_classif">
+       <attribute name="title">
+        <string>Classification</string>
+       </attribute>
+       <layout class="QVBoxLayout" name="verticalLayout">
+        <item>
+         <spacer name="verticalSpacer_4">
+          <property name="orientation">
+           <enum>Qt::Vertical</enum>
+          </property>
+          <property name="sizeHint" stdset="0">
+           <size>
+            <width>20</width>
+            <height>40</height>
+           </size>
+          </property>
+         </spacer>
+        </item>
+        <item>
+         <widget class="QLabel" name="label_8">
+          <property name="text">
+           <string>Labels</string>
+          </property>
+         </widget>
+        </item>
+        <item>
+         <widget class="QFrame" name="frame_2">
+          <property name="frameShape">
+           <enum>QFrame::StyledPanel</enum>
+          </property>
+          <property name="frameShadow">
+           <enum>QFrame::Raised</enum>
+          </property>
+          <layout class="QHBoxLayout" name="horizontalLayout_12">
+           <item>
+            <layout class="QGridLayout" name="gridLayout_3">
+             <property name="sizeConstraint">
+              <enum>QLayout::SetDefaultConstraint</enum>
+             </property>
+             <item row="0" column="1">
+              <widget class="QLabel" name="label_10">
+               <property name="font">
+                <font>
+                 <weight>75</weight>
+                 <bold>true</bold>
+                </font>
+               </property>
+               <property name="text">
+                <string>Color</string>
+               </property>
+               <property name="alignment">
+                <set>Qt::AlignCenter</set>
+               </property>
+              </widget>
+             </item>
+             <item row="0" column="3">
+              <widget class="QLabel" name="label_12">
+               <property name="font">
+                <font>
+                 <weight>75</weight>
+                 <bold>true</bold>
+                </font>
+               </property>
+               <property name="text">
+                <string>Training set</string>
+               </property>
+               <property name="alignment">
+                <set>Qt::AlignCenter</set>
+               </property>
+              </widget>
+             </item>
+             <item row="0" column="0">
+              <widget class="QLabel" name="label_9">
+               <property name="font">
+                <font>
+                 <weight>75</weight>
+                 <bold>true</bold>
+                </font>
+               </property>
+               <property name="text">
+                <string>Name</string>
+               </property>
+               <property name="alignment">
+                <set>Qt::AlignLeading|Qt::AlignLeft|Qt::AlignVCenter</set>
+               </property>
+              </widget>
+             </item>
+             <item row="0" column="2">
+              <spacer name="horizontalSpacer_2">
+               <property name="orientation">
+                <enum>Qt::Horizontal</enum>
+               </property>
+               <property name="sizeHint" stdset="0">
+                <size>
+                 <width>40</width>
+                 <height>20</height>
+                </size>
+               </property>
+              </spacer>
+             </item>
+             <item row="0" column="5">
+              <widget class="QPushButton" name="add_new_label">
+               <property name="sizePolicy">
+                <sizepolicy hsizetype="Fixed" vsizetype="Fixed">
+                 <horstretch>0</horstretch>
+                 <verstretch>0</verstretch>
+                </sizepolicy>
+               </property>
+               <property name="font">
+                <font>
+                 <weight>75</weight>
+                 <bold>true</bold>
+                </font>
+               </property>
+               <property name="text">
+                <string>Add new</string>
+               </property>
+              </widget>
+             </item>
+             <item row="0" column="4">
+              <spacer name="horizontalSpacer_3">
+               <property name="orientation">
+                <enum>Qt::Horizontal</enum>
+               </property>
+               <property name="sizeHint" stdset="0">
+                <size>
+                 <width>40</width>
+                 <height>20</height>
+                </size>
+               </property>
+              </spacer>
+             </item>
+            </layout>
+           </item>
+          </layout>
+         </widget>
+        </item>
+        <item>
+         <spacer name="verticalSpacer_5">
+          <property name="orientation">
+           <enum>Qt::Vertical</enum>
+          </property>
+          <property name="sizeHint" stdset="0">
+           <size>
+            <width>20</width>
+            <height>40</height>
+           </size>
+          </property>
+         </spacer>
+        </item>
+        <item>
+         <widget class="QLabel" name="label_11">
+          <property name="text">
+           <string>Training</string>
+          </property>
+         </widget>
+        </item>
+        <item>
+         <widget class="QFrame" name="frame_3">
+          <property name="frameShape">
+           <enum>QFrame::StyledPanel</enum>
+          </property>
+          <property name="frameShadow">
+           <enum>QFrame::Raised</enum>
+          </property>
+          <layout class="QHBoxLayout" name="horizontalLayout_6">
+           <item>
+            <layout class="QGridLayout" name="gridLayout_2">
+             <item row="1" column="0">
+              <widget class="QPushButton" name="train">
+               <property name="font">
+                <font>
+                 <weight>75</weight>
+                 <bold>true</bold>
+                </font>
+               </property>
+               <property name="text">
+                <string>Train</string>
+               </property>
+              </widget>
+             </item>
+             <item row="1" column="1">
+              <widget class="QSpinBox" name="number_of_trials">
+               <property name="suffix">
+                <string> trials</string>
+               </property>
+               <property name="minimum">
+                <number>10</number>
+               </property>
+               <property name="maximum">
+                <number>1000000000</number>
+               </property>
+               <property name="singleStep">
+                <number>10</number>
+               </property>
+               <property name="value">
+                <number>300</number>
+               </property>
+              </widget>
+             </item>
+             <item row="0" column="0">
+              <widget class="QPushButton" name="reset_training_sets">
+               <property name="text">
+                <string>Reset training sets</string>
+               </property>
+              </widget>
+             </item>
+             <item row="0" column="1">
+              <widget class="QPushButton" name="validate_selection">
+               <property name="text">
+                <string>Validate selection for training</string>
+               </property>
+              </widget>
+             </item>
+            </layout>
+           </item>
+          </layout>
+         </widget>
+        </item>
+        <item>
+         <spacer name="verticalSpacer_11">
+          <property name="orientation">
+           <enum>Qt::Vertical</enum>
+          </property>
+          <property name="sizeHint" stdset="0">
+           <size>
+            <width>20</width>
+            <height>40</height>
+           </size>
+          </property>
+         </spacer>
+        </item>
+       </layout>
+      </widget>
+      <widget class="QWidget" name="tab_advanced">
+       <property name="sizePolicy">
+        <sizepolicy hsizetype="Preferred" vsizetype="Minimum">
+         <horstretch>0</horstretch>
+         <verstretch>0</verstretch>
+        </sizepolicy>
+       </property>
+       <attribute name="title">
+        <string>Advanced</string>
+       </attribute>
+       <layout class="QVBoxLayout" name="verticalLayout_3">
+        <item>
+         <spacer name="verticalSpacer_2">
+          <property name="orientation">
+           <enum>Qt::Vertical</enum>
+          </property>
+          <property name="sizeHint" stdset="0">
+           <size>
+            <width>20</width>
+            <height>40</height>
+           </size>
+          </property>
+         </spacer>
+        </item>
+        <item>
+         <layout class="QHBoxLayout" name="horizontalLayout_4" stretch="0,2">
+          <item>
+           <widget class="QLabel" name="label_2">
+            <property name="sizePolicy">
+             <sizepolicy hsizetype="MinimumExpanding" vsizetype="Preferred">
+              <horstretch>0</horstretch>
+              <verstretch>0</verstretch>
+             </sizepolicy>
+            </property>
+            <property name="text">
+             <string>Features:</string>
+            </property>
+           </widget>
+          </item>
+          <item>
+           <widget class="QComboBox" name="selected_feature">
+            <property name="sizePolicy">
+             <sizepolicy hsizetype="MinimumExpanding" vsizetype="Fixed">
+              <horstretch>0</horstretch>
+              <verstretch>0</verstretch>
+             </sizepolicy>
+            </property>
+           </widget>
+          </item>
+         </layout>
+        </item>
+        <item>
+         <spacer name="verticalSpacer_10">
+          <property name="orientation">
+           <enum>Qt::Vertical</enum>
+          </property>
+          <property name="sizeHint" stdset="0">
+           <size>
+            <width>20</width>
+            <height>40</height>
+           </size>
+          </property>
+         </spacer>
+        </item>
+        <item>
+         <widget class="QSlider" name="feature_weight">
+          <property name="sizePolicy">
+           <sizepolicy hsizetype="MinimumExpanding" vsizetype="Fixed">
+            <horstretch>0</horstretch>
+            <verstretch>0</verstretch>
+           </sizepolicy>
+          </property>
+          <property name="maximum">
+           <number>1000</number>
+          </property>
+          <property name="pageStep">
+           <number>50</number>
+          </property>
+          <property name="value">
+           <number>500</number>
+          </property>
+          <property name="orientation">
+           <enum>Qt::Horizontal</enum>
+          </property>
+          <property name="invertedAppearance">
+           <bool>false</bool>
+          </property>
+          <property name="invertedControls">
+           <bool>false</bool>
+          </property>
+          <property name="tickPosition">
+           <enum>QSlider::TicksBothSides</enum>
+          </property>
+         </widget>
+        </item>
+        <item>
+         <spacer name="verticalSpacer_7">
+          <property name="orientation">
+           <enum>Qt::Vertical</enum>
+          </property>
+          <property name="sizeHint" stdset="0">
+           <size>
+            <width>20</width>
+            <height>40</height>
+           </size>
+          </property>
+         </spacer>
+        </item>
+        <item>
+         <layout class="QGridLayout" name="gridLayout">
+          <item row="0" column="1">
+           <spacer name="horizontalSpacer">
+            <property name="orientation">
+             <enum>Qt::Horizontal</enum>
+            </property>
+            <property name="sizeHint" stdset="0">
+             <size>
+              <width>40</width>
+              <height>20</height>
+             </size>
+            </property>
+           </spacer>
+          </item>
+          <item row="0" column="0">
+           <widget class="QLabel" name="label_3">
+            <property name="font">
+             <font>
+              <weight>75</weight>
+              <bold>true</bold>
+             </font>
+            </property>
+            <property name="text">
+             <string>Label</string>
+            </property>
+           </widget>
+          </item>
+          <item row="0" column="2">
+           <widget class="QLabel" name="label_4">
+            <property name="font">
+             <font>
+              <weight>75</weight>
+              <bold>true</bold>
+             </font>
+            </property>
+            <property name="text">
+             <string>Effect</string>
+            </property>
+           </widget>
+          </item>
+          <item row="0" column="3">
+           <spacer name="horizontalSpacer_4">
+            <property name="orientation">
+             <enum>Qt::Horizontal</enum>
+            </property>
+            <property name="sizeHint" stdset="0">
+             <size>
+              <width>40</width>
+              <height>20</height>
+             </size>
+            </property>
+           </spacer>
+          </item>
+         </layout>
+        </item>
+        <item>
+         <spacer name="verticalSpacer_9">
+          <property name="orientation">
+           <enum>Qt::Vertical</enum>
+          </property>
+          <property name="sizeHint" stdset="0">
+           <size>
+            <width>20</width>
+            <height>40</height>
+           </size>
+          </property>
+         </spacer>
+        </item>
+       </layout>
+      </widget>
+     </widget>
+    </item>
+    <item>
+     <spacer name="verticalSpacer_8">
+      <property name="orientation">
+       <enum>Qt::Vertical</enum>
+      </property>
+      <property name="sizeHint" stdset="0">
+       <size>
+        <width>20</width>
+        <height>41</height>
+       </size>
+      </property>
+     </spacer>
+    </item>
+    <item>
+     <layout class="QHBoxLayout" name="horizontalLayout_9">
+      <item>
+       <widget class="QPushButton" name="run_smoothed">
+        <property name="text">
+         <string>Run (smoothed)</string>
+        </property>
+       </widget>
+      </item>
+      <item>
+       <widget class="QPushButton" name="run">
+        <property name="font">
+         <font>
+          <weight>50</weight>
+          <bold>false</bold>
+         </font>
+        </property>
+        <property name="text">
+         <string>Run (quick)</string>
+        </property>
+       </widget>
+      </item>
+      <item>
+       <widget class="QFrame" name="frame">
+        <property name="frameShape">
+         <enum>QFrame::StyledPanel</enum>
+        </property>
+        <property name="frameShadow">
+         <enum>QFrame::Raised</enum>
+        </property>
+        <layout class="QVBoxLayout" name="verticalLayout_2">
+         <item>
+          <widget class="QPushButton" name="run_graphcut">
+           <property name="font">
+            <font>
+             <weight>50</weight>
+             <bold>false</bold>
+            </font>
+           </property>
+           <property name="text">
+            <string>Run (graphcut)</string>
+           </property>
+          </widget>
+         </item>
+         <item>
+          <layout class="QHBoxLayout" name="horizontalLayout_2">
+           <item>
+            <widget class="QLabel" name="label_5">
+             <property name="text">
+              <string>Weight =</string>
+             </property>
+            </widget>
+           </item>
+           <item>
+            <widget class="QDoubleSpinBox" name="smoothingDoubleSpinBox">
+             <property name="sizePolicy">
+              <sizepolicy hsizetype="Minimum" vsizetype="Fixed">
+               <horstretch>0</horstretch>
+               <verstretch>0</verstretch>
+              </sizepolicy>
+             </property>
+             <property name="toolTip">
+              <string>Smoothing strenght</string>
+             </property>
+             <property name="maximum">
+              <double>10.000000000000000</double>
+             </property>
+             <property name="singleStep">
+              <double>0.100000000000000</double>
+             </property>
+             <property name="value">
+              <double>0.500000000000000</double>
+             </property>
+            </widget>
+           </item>
+          </layout>
+         </item>
+         <item>
+          <layout class="QHBoxLayout" name="horizontalLayout_3">
+           <item>
+            <widget class="QLabel" name="label_6">
+             <property name="text">
+              <string>Subdivisions =</string>
+             </property>
+            </widget>
+           </item>
+           <item>
+            <widget class="QSpinBox" name="subdivisionsSpinBox">
+             <property name="minimum">
+              <number>1</number>
+             </property>
+             <property name="maximum">
+              <number>10000</number>
+             </property>
+             <property name="value">
+              <number>16</number>
+             </property>
+            </widget>
+           </item>
+          </layout>
+         </item>
+        </layout>
+       </widget>
+      </item>
+     </layout>
+    </item>
+    <item>
+     <spacer name="verticalSpacer_6">
+      <property name="orientation">
+       <enum>Qt::Vertical</enum>
+      </property>
+      <property name="sizeHint" stdset="0">
+       <size>
+        <width>20</width>
+        <height>40</height>
+       </size>
+      </property>
+     </spacer>
+    </item>
+    <item>
+     <layout class="QHBoxLayout" name="horizontalLayout_5">
+      <item>
+       <widget class="QPushButton" name="save">
+        <property name="font">
+         <font>
+          <weight>75</weight>
+          <bold>true</bold>
+         </font>
+        </property>
+        <property name="text">
+         <string>Save classified point set</string>
+        </property>
+       </widget>
+      </item>
+      <item>
+       <widget class="QPushButton" name="generate_point_set_items">
+        <property name="text">
+         <string>Generate point set items</string>
+        </property>
+       </widget>
+      </item>
+     </layout>
+    </item>
+    <item>
+     <spacer name="verticalSpacer_14">
+      <property name="orientation">
+       <enum>Qt::Vertical</enum>
+      </property>
+      <property name="sizeHint" stdset="0">
+       <size>
+        <width>20</width>
+        <height>40</height>
+       </size>
+      </property>
+     </spacer>
+    </item>
+   </layout>
+  </widget>
+ </widget>
+ <resources/>
+ <connections/>
+</ui>

Polyhedron/demo/Polyhedron/Scene_point_set_classification_item.cpp

@@ -0,0 +1,535 @@
+#include "Scene_point_set_classification_item.h"
+#include "Color_ramp.h"
+
+#include <CGAL/Timer.h>
+#include <CGAL/Memory_sizer.h>
+
+#include <CGAL/Three/Viewer_interface.h>
+
+#include <QObject>
+#include <QMenu>
+#include <QGLViewer/manipulatedCameraFrame.h>
+ 
+#include <set>
+#include <stack>
+#include <algorithm>
+#include <boost/array.hpp>
+
+Scene_point_set_classification_item::Scene_point_set_classification_item(PSC* psc)
+: Scene_item(NbOfVbos,NbOfVaos),
+  m_points (NULL),
+  m_psc (psc),
+  m_trainer (new Trainer(*m_psc))
+{
+  setRenderingMode(PointsPlusNormals);
+  m_nb_scales = 5;
+  m_nb_trials = 300;
+  m_smoothing = 0.5;
+  m_subdivisions = 16;
+  m_index_color = 1;
+  is_selected = true;
+  nb_points = 0;
+  nb_lines = 0;
+}
+
+Scene_point_set_classification_item::Scene_point_set_classification_item(Scene_points_with_normal_item* points)
+  : Scene_item(NbOfVbos,NbOfVaos),
+    m_points (points),
+    m_psc (NULL),
+    m_trainer (NULL)
+{
+  setRenderingMode(PointsPlusNormals);
+  m_nb_scales = 5;
+  m_index_color = 1;
+  m_nb_trials = 300;
+  m_smoothing = 0.5;
+  m_subdivisions = 16;
+  
+  reset_indices();
+  
+  m_psc = new PSC(*(m_points->point_set()), m_points->point_set()->point_map());
+  m_trainer = new Trainer(*m_psc);
+
+  Label_handle ground = m_psc->add_label("ground");
+  Label_handle vegetation = m_psc->add_label("vegetation");
+  Label_handle roof = m_psc->add_label("roof");
+  Label_handle facade = m_psc->add_label("facade");
+  m_labels.push_back (std::make_pair(ground, QColor(245, 180, 0)));
+  m_labels.push_back (std::make_pair(vegetation, QColor(0, 255, 27)));
+  m_labels.push_back (std::make_pair(roof, QColor(255, 0, 170)));
+  m_labels.push_back (std::make_pair(facade, QColor(100, 0, 255)));
+  
+  is_selected = true;
+  nb_points = 0;
+  nb_lines = 0;
+}
+
+
+// Copy constructor
+Scene_point_set_classification_item::Scene_point_set_classification_item(const Scene_point_set_classification_item&)
+  :Scene_item(NbOfVbos,NbOfVaos), // do not call superclass' copy constructor
+   m_points (NULL),
+   m_psc (NULL),
+   m_trainer (NULL)
+{
+  setRenderingMode(PointsPlusNormals);
+  m_nb_scales = 5;
+  m_index_color = 1;
+  m_nb_trials = 300;
+  m_smoothing = 0.5;
+  m_subdivisions = 16;
+    
+  nb_points = 0;
+  invalidateOpenGLBuffers();
+}
+
+Scene_point_set_classification_item::~Scene_point_set_classification_item()
+{
+  if (m_psc != NULL)
+    delete m_psc;
+  if (m_trainer != NULL)
+    delete m_trainer;
+}
+
+void Scene_point_set_classification_item::initializeBuffers(CGAL::Three::Viewer_interface *viewer) const
+{
+  compute_normals_and_vertices();
+  //vao for the edges
+  {
+    program = getShaderProgram(PROGRAM_NO_SELECTION, viewer);
+    program->bind();
+
+    vaos[Edges]->bind();
+    buffers[Edges_vertices].bind();
+    buffers[Edges_vertices].allocate(positions_lines.data(),
+                                     static_cast<int>(positions_lines.size()*sizeof(double)));
+    program->enableAttributeArray("vertex");
+    program->setAttributeBuffer("vertex",GL_DOUBLE,0,3);
+    buffers[Edges_vertices].release();
+
+    vaos[Edges]->release();
+    nb_lines = positions_lines.size();
+    positions_lines.resize(0);
+    std::vector<double>(positions_lines).swap(positions_lines);
+    program->release();
+  }
+  //vao for the points
+  {
+    program = getShaderProgram(PROGRAM_NO_SELECTION, viewer);
+    program->bind();
+
+    vaos[ThePoints]->bind();
+    buffers[Points_vertices].bind();
+    buffers[Points_vertices].allocate(positions_points.data(),
+                                      static_cast<int>(positions_points.size()*sizeof(double)));
+    program->enableAttributeArray("vertex");
+    program->setAttributeBuffer("vertex",GL_DOUBLE,0,3);
+    buffers[Points_vertices].release();
+    if (!(colors_points.empty()))
+      {
+        buffers[Points_colors].bind();
+        buffers[Points_colors].allocate (colors_points.data(),
+                                         static_cast<int>(colors_points.size()*sizeof(double)));
+        program->enableAttributeArray("colors");
+        program->setAttributeBuffer("colors",GL_DOUBLE,0,3);
+        buffers[Points_colors].release();
+        std::vector<double>(colors_points).swap(colors_points);
+      }
+
+    vaos[ThePoints]->release();
+    nb_points = positions_points.size();
+    positions_points.resize(0);
+    std::vector<double>(positions_points).swap(positions_points);
+    program->release();
+  }
+  are_buffers_filled = true;
+}
+
+void Scene_point_set_classification_item::compute_normals_and_vertices() const
+{
+  int index_color = real_index_color();
+    
+  positions_lines.resize(0);
+  positions_points.resize(0);
+  colors_points.resize(0);
+  if (index_color == 9) // Show clusters centroids
+    {
+      // positions_points.reserve(m_psc->clusters().size() * 3);
+      // colors_points.reserve(m_psc->clusters().size() * 3);
+
+      // for (std::size_t i = 0; i < m_psc->clusters().size(); ++ i)
+      //   {
+      //     const Kernel::Point_3& p = m_psc->clusters()[i].centroid;
+      //     positions_points.push_back(p.x());
+      //     positions_points.push_back(p.y());
+      //     positions_points.push_back(p.z());
+
+      //     for (std::set<std::size_t>::iterator it = m_psc->clusters()[i].neighbors.begin ();
+      //          it != m_psc->clusters()[i].neighbors.end (); ++ it)
+      //       {
+      //         const Kernel::Point_3& q = m_psc->clusters()[*it].centroid;
+      //         positions_lines.push_back(p.x());
+      //         positions_lines.push_back(p.y());
+      //         positions_lines.push_back(p.z());
+
+      //         positions_lines.push_back(q.x());
+      //         positions_lines.push_back(q.y());
+      //         positions_lines.push_back(q.z());
+      //       }
+      //   }
+    }
+  else // Show points
+    {
+      positions_points.reserve(m_points->point_set()->size() * 3);
+      colors_points.reserve(m_points->point_set()->size() * 3);
+
+      for (Point_set::const_iterator it = m_points->point_set()->begin();
+           it != m_points->point_set()->end(); ++ it)
+        {
+          const Kernel::Point_3& p = m_points->point_set()->point(*it);
+          positions_points.push_back(p.x());
+          positions_points.push_back(p.y());
+          positions_points.push_back(p.z());
+        }
+    }
+
+  // Colors
+  static Color_ramp ramp;
+  ramp.build_red();
+
+  if (index_color == -1) // item color
+    {
+
+      for (Point_set::const_iterator it = m_points->point_set()->begin();
+           it != m_points->point_set()->first_selected(); ++ it)
+        {
+          colors_points.push_back (color().redF());
+          colors_points.push_back (color().greenF());
+          colors_points.push_back (color().blueF());
+        }
+    }
+  else if (index_color == 0) // real colors
+    {
+
+      for (Point_set::const_iterator it = m_points->point_set()->begin();
+           it != m_points->point_set()->first_selected(); ++ it)
+        {
+          colors_points.push_back (m_points->point_set()->red(*it));
+          colors_points.push_back (m_points->point_set()->green(*it));
+          colors_points.push_back (m_points->point_set()->blue(*it));
+        }
+    }
+  else if (index_color == 1) // classif
+    {
+      std::map<Label_handle, QColor> map_colors;
+      for (std::size_t i = 0; i < m_labels.size(); ++ i)
+        map_colors.insert (m_labels[i]);
+          
+      for (Point_set::const_iterator it = m_points->point_set()->begin();
+           it != m_points->point_set()->first_selected(); ++ it)
+        {
+          QColor color (0, 0, 0);
+          Label_handle c = m_psc->label_of(*it);
+          
+          if (c != Label_handle())
+            color = map_colors[c];
+
+          colors_points.push_back ((double)(color.red()) / 255.);
+          colors_points.push_back ((double)(color.green()) / 255.);
+          colors_points.push_back ((double)(color.blue()) / 255.);
+        }
+    }
+  else if (index_color == 2) // training
+    {
+      std::map<Label_handle, QColor> map_colors;
+      for (std::size_t i = 0; i < m_labels.size(); ++ i)
+        map_colors.insert (m_labels[i]);
+          
+      for (Point_set::const_iterator it = m_points->point_set()->begin();
+           it != m_points->point_set()->first_selected(); ++ it)
+        {
+          QColor color (0, 0, 0);
+          Label_handle c = m_trainer->training_label_of(*it);
+          Label_handle c2 = m_psc->label_of(*it);
+          
+          if (c != Label_handle())
+            color = map_colors[c];
+          double div = 255.;
+          if (c != c2)
+            div = 400.;
+          
+          colors_points.push_back ((double)(color.red()) / div);
+          colors_points.push_back ((double)(color.green()) / div);
+          colors_points.push_back ((double)(color.blue()) / div);
+        }
+    }
+  else
+    {
+      Feature_handle att = m_psc->feature(index_color - 3);
+      double weight = att->weight();
+      att->set_weight(att->max);
+      for (Point_set::const_iterator it = m_points->point_set()->begin();
+           it != m_points->point_set()->first_selected(); ++ it)
+        {
+          colors_points.push_back (ramp.r(att->normalized(*it)));
+          colors_points.push_back (ramp.g(att->normalized(*it)));
+          colors_points.push_back (ramp.b(att->normalized(*it)));
+        }
+      att->set_weight(weight);
+    }
+
+  for (Point_set::const_iterator it = m_points->point_set()->first_selected();
+       it != m_points->point_set()->end(); ++ it)
+    {
+      colors_points.push_back (1.);
+      colors_points.push_back (0.);
+      colors_points.push_back (0.);
+    }
+}
+
+
+// Duplicates scene item
+Scene_point_set_classification_item*
+Scene_point_set_classification_item::clone() const
+{
+  return new Scene_point_set_classification_item(*this);
+}
+
+// Write point set to .PLY file
+bool Scene_point_set_classification_item::write_ply_point_set(std::ostream& stream)
+{
+  if (m_psc->number_of_features() == 0)
+    return false;
+  
+  stream.precision (std::numeric_limits<double>::digits10 + 2);
+
+  reset_indices();
+
+  std::vector<Color> colors;
+  for (std::size_t i = 0; i < m_labels.size(); ++ i)
+    {
+      Color c = {{ (unsigned char)(m_labels[i].second.red()),
+                   (unsigned char)(m_labels[i].second.green()),
+                   (unsigned char)(m_labels[i].second.blue()) }};
+      colors.push_back (c);
+    }
+  
+  m_psc->write_classification_to_ply (stream);
+  return true;
+}
+
+QString
+Scene_point_set_classification_item::toolTip() const
+{
+  return QObject::tr("<p><b>%1</b><br />"
+                     "<i>Point set classification</i></p>"
+                     "<p>Number of points: %2</p>")
+    .arg(name())
+    .arg(m_points->point_set()->size());
+}
+
+bool Scene_point_set_classification_item::supportsRenderingMode(RenderingMode m) const 
+{
+  return m==PointsPlusNormals;
+
+}
+
+void Scene_point_set_classification_item::drawSplats(CGAL::Three::Viewer_interface*) const
+{
+}
+
+void Scene_point_set_classification_item::drawEdges(CGAL::Three::Viewer_interface* viewer) const
+{
+  if(!are_buffers_filled)
+    initializeBuffers(viewer);
+  vaos[Edges]->bind();
+  program=getShaderProgram(PROGRAM_NO_SELECTION);
+  attribBuffers(viewer,PROGRAM_NO_SELECTION);
+  program->bind();
+  program->setAttributeValue("colors", this->color());
+  viewer->glDrawArrays(GL_LINES, 0,
+                       static_cast<GLsizei>((nb_lines)/3));
+  vaos[Edges]->release();
+  program->release();
+}
+
+void Scene_point_set_classification_item::drawPoints(CGAL::Three::Viewer_interface* viewer) const
+{
+  if(!are_buffers_filled)
+    initializeBuffers(viewer);
+  GLfloat point_size;
+  viewer->glGetFloatv(GL_POINT_SIZE, &point_size);
+  if (m_index_color == 9)
+    viewer->glPointSize(15.f);
+  else
+    viewer->glPointSize(3.f);
+
+  vaos[ThePoints]->bind();
+  program=getShaderProgram(PROGRAM_NO_SELECTION);
+  attribBuffers(viewer,PROGRAM_NO_SELECTION);
+  program->bind();
+  if (colors_points.empty())
+    program->setAttributeValue("colors", this->color());
+  viewer->glDrawArrays(GL_POINTS, 0,
+                       static_cast<GLsizei>(((std::size_t)(nb_points)/3)));
+  vaos[ThePoints]->release();
+  program->release();
+  viewer->glPointSize(point_size);
+}
+
+bool
+Scene_point_set_classification_item::isEmpty() const
+{
+  return false;
+}
+
+void
+Scene_point_set_classification_item::compute_bbox() const
+{
+  if (m_points->point_set()->empty())
+    return;
+  
+  double xmin = std::numeric_limits<double>::max();
+  double ymin = std::numeric_limits<double>::max();
+  double zmin = std::numeric_limits<double>::max();
+  double xmax = -std::numeric_limits<double>::max();
+  double ymax = -std::numeric_limits<double>::max();
+  double zmax = -std::numeric_limits<double>::max();
+
+  for (Point_set::const_iterator it = m_points->point_set()->begin();
+       it != m_points->point_set()->end(); ++ it)
+    {
+      const Kernel::Point_3& p = m_points->point_set()->point(*it);
+      xmin = (std::min)(xmin, p.x());
+      ymin = (std::min)(ymin, p.y());
+      zmin = (std::min)(zmin, p.z());
+      xmax = (std::max)(xmax, p.x());
+      ymax = (std::max)(ymax, p.y());
+      zmax = (std::max)(zmax, p.z());
+    }
+  
+  _bbox = Bbox(xmin,ymin,zmin,
+               xmax,ymax,zmax);
+}
+
+QMenu* Scene_point_set_classification_item::contextMenu()
+{
+  QMenu* menu = Scene_item::contextMenu();
+
+  return menu;
+}
+
+void Scene_point_set_classification_item::setRenderingMode(RenderingMode m)
+{
+  Scene_item::setRenderingMode(m);
+}
+
+void Scene_point_set_classification_item::invalidateOpenGLBuffers()
+{
+  are_buffers_filled = false;
+  compute_bbox();
+}
+
+void Scene_point_set_classification_item::change_color (int index)
+{
+  m_index_color = index;
+  invalidateOpenGLBuffers();
+}
+
+int Scene_point_set_classification_item::real_index_color() const
+{
+  int out = m_index_color;
+  
+  if (out == 0 && !(m_points->point_set()->has_colors()))
+    out = -1;
+  return out;
+}
+
+void Scene_point_set_classification_item::reset_indices ()
+{
+  Point_set::Property_map<Point_set::Index> indices
+    = m_points->point_set()->property_map<Point_set::Index>("index").first;
+
+  m_points->point_set()->unselect_all();
+  Point_set::Index idx;
+  ++ idx;
+  for (std::size_t i = 0; i < m_points->point_set()->size(); ++ i)
+    *(indices.begin() + i) = idx ++;
+}
+
+void Scene_point_set_classification_item::compute_features ()
+{
+  Q_ASSERT (!(m_points->point_set()->empty()));
+  Q_ASSERT (m_psc != NULL);
+  m_psc->clear_features();
+  reset_indices();
+  
+  std::cerr << "Computing features with " << m_nb_scales << " scale(s)" << std::endl;
+  if (m_psc->number_of_features() != 0)
+    m_psc->clear();
+  compute_bbox();
+
+  bool normals = m_points->point_set()->has_normal_map();
+  bool colors = m_points->point_set()->has_colors();
+  Point_set::Property_map<boost::uint8_t> echo_map;
+  bool echo;
+  boost::tie (echo_map, echo) = m_points->point_set()->template property_map<boost::uint8_t>("echo");
+
+  if (!normals && !colors && !echo)
+    m_psc->generate_features (m_nb_scales);
+  else if (!normals && !colors && echo)
+    m_psc->generate_features (m_nb_scales, CGAL::Default(), CGAL::Default(), echo_map);
+  else if (!normals && colors && !echo)
+    m_psc->generate_features (m_nb_scales, CGAL::Default(), Color_map(m_points->point_set()));
+  else if (!normals && colors && echo)
+    m_psc->generate_features (m_nb_scales, CGAL::Default(), Color_map(m_points->point_set()), echo_map);
+  else if (normals && !colors && !echo)
+    m_psc->generate_features (m_nb_scales, m_points->point_set()->normal_map());
+  else if (normals && !colors && echo)
+    m_psc->generate_features (m_nb_scales, m_points->point_set()->normal_map(), CGAL::Default(), echo_map);
+  else if (normals && colors && !echo)
+    m_psc->generate_features (m_nb_scales, m_points->point_set()->normal_map(), Color_map(m_points->point_set()));
+  else
+    m_psc->generate_features (m_nb_scales, m_points->point_set()->normal_map(), Color_map(m_points->point_set()), echo_map);
+
+}
+
+
+
+void Scene_point_set_classification_item::train()
+{
+  if (m_psc->number_of_features() == 0)
+    {
+      std::cerr << "Error: features not computed" << std::endl;
+      return;
+    }
+  reset_indices();
+  
+  m_trainer->train(m_nb_trials);
+  m_psc->run();
+  m_psc->info();
+}
+
+bool Scene_point_set_classification_item::run (int method)
+{
+  if (m_psc->number_of_features() == 0)
+    {
+      std::cerr << "Error: features not computed" << std::endl;
+      return false;
+    }
+  reset_indices();
+
+  if (method == 0)
+    {
+      m_psc->run();
+      m_psc->info();
+    }
+  else if (method == 1)
+    m_psc->run_with_local_smoothing (m_psc->neighborhood().range_neighbor_query(m_psc->radius_neighbors()));
+  else if (method == 2)
+    m_psc->run_with_graphcut (m_psc->neighborhood().k_neighbor_query(12), m_smoothing, m_subdivisions);
+  invalidateOpenGLBuffers();
+  Q_EMIT itemChanged();
+  
+  return true;
+}

Polyhedron/demo/Polyhedron/Scene_point_set_classification_item.h

@@ -0,0 +1,363 @@
+#ifndef POINT_SET_CLASSIFICATION_ITEM_H
+#define POINT_SET_CLASSIFICATION_ITEM_H
+
+//#define CGAL_DO_NOT_USE_BOYKOV_KOLMOGOROV_MAXFLOW_SOFTWARE
+#define CGAL_CLASSIFICATION_VERBOSE
+
+#include <CGAL/Three/Scene_item.h>
+#include <CGAL/Point_set_classifier.h>
+#include <CGAL/Classification/Trainer.h>
+#include <CGAL/Classification/Feature_base.h>
+#include <CGAL/Classification/Feature/Vertical_dispersion.h>
+#include <CGAL/Classification/Feature/Elevation.h>
+#include <CGAL/Classification/Feature/Verticality.h>
+#include <CGAL/Classification/Feature/Distance_to_plane.h>
+#include <CGAL/Classification/Feature/Hsv.h>
+#include <CGAL/Classification/Feature/Echo_scatter.h>
+#include <CGAL/Classification/Feature/Eigen.h>
+
+
+#include "Scene_point_set_classification_item_config.h"
+#include "Scene_points_with_normal_item.h"
+#include "Polyhedron_type_fwd.h"
+#include "Kernel_type.h"
+#include "Point_set_3.h"
+
+#include <iostream>
+
+
+class QMenu;
+class QAction;
+
+// This class represents a point set in the OpenGL scene
+class SCENE_POINT_SET_CLASSIFICATION_ITEM_EXPORT Scene_point_set_classification_item
+  : public CGAL::Three::Scene_item
+{
+  Q_OBJECT
+
+  
+  template <typename PointSet>
+    class Point_set_color_map
+  {
+  public:
+    typedef CGAL::Classification::RGB_Color Color;
+    typedef typename PointSet::Index key_type;
+    typedef Color value_type;
+    typedef value_type& reference;
+    typedef boost::lvalue_property_map_tag category;
+    
+    PointSet* ps;
+
+    Point_set_color_map(PointSet* ps = NULL)
+      : ps(ps) { }
+
+    friend value_type get (const Point_set_color_map& pm, const key_type& i)
+    {
+      Color out = {{ (unsigned char)(255 * pm.ps->red(i)),
+                     (unsigned char)(255 * pm.ps->green(i)),
+                     (unsigned char)(255 * pm.ps->blue(i)) }};
+                     
+      return out;
+    }
+  };
+
+ public:
+  typedef Kernel::Point_3 Point_3;
+  typedef Kernel::Vector_3 Vector_3;
+  typedef CGAL::Classification::RGB_Color Color;
+  
+  typedef Point_set::Point_map Point_map;
+  typedef Point_set::Vector_map Vector_map;
+  typedef Point_set_color_map<Point_set> Color_map;
+
+  typedef CGAL::Point_set_classifier<Kernel, Point_set, Point_map>                     PSC;
+  typedef CGAL::Classification::Trainer<Point_set, Point_map>        Trainer;
+  typedef CGAL::Classification::Label_handle                                               Label_handle;
+  typedef CGAL::Classification::Feature_handle                                             Feature_handle;
+  typedef CGAL::Classification::Feature::Vertical_dispersion<Kernel, Point_set, Point_map> Dispersion;
+  typedef CGAL::Classification::Feature::Elevation<Kernel, Point_set, Point_map>           Elevation;
+  typedef CGAL::Classification::Feature::Verticality<Kernel, Point_set, Point_map>         Verticality;
+  typedef CGAL::Classification::Feature::Distance_to_plane<Kernel, Point_set, Point_map>   Distance_to_plane;
+
+
+ public:
+  
+  Scene_point_set_classification_item(PSC* psc = NULL);
+  Scene_point_set_classification_item(Scene_points_with_normal_item* points);
+  Scene_point_set_classification_item(const Scene_point_set_classification_item& toCopy);
+  ~Scene_point_set_classification_item();
+  
+  Scene_point_set_classification_item* clone() const;
+  
+  // Function to override the context menu
+  QMenu* contextMenu();
+
+  // IO
+  bool write_ply_point_set(std::ostream& out);
+  // Function for displaying meta-data of the item
+  virtual QString toolTip() const;
+
+  virtual void invalidateOpenGLBuffers();
+
+  // Indicate if rendering mode is supported
+  virtual bool supportsRenderingMode(RenderingMode m) const;
+
+  virtual void drawEdges(CGAL::Three::Viewer_interface* viewer) const;
+  virtual void drawPoints(CGAL::Three::Viewer_interface*) const;
+  virtual void drawSplats(CGAL::Three::Viewer_interface*) const;
+  
+
+  // Gets dimensions
+  virtual bool isFinite() const { return true; }
+  virtual bool isEmpty() const;
+  virtual void compute_bbox() const;
+
+  virtual void setRenderingMode(RenderingMode m);
+
+  void change_color (int index);
+  int real_index_color() const;
+
+  void reset_indices();
+  void compute_features ();
+  bool run (int method);
+
+  template <typename Item>
+  void generate_point_set_items(std::vector<Item*>& items,
+                                const char* name)
+  {
+    std::map<Label_handle, std::size_t> map_labels;
+    for (std::size_t i = 0; i < m_labels.size(); ++ i)
+      {
+        items.push_back (new Item);
+        items.back()->setName (QString("%1 (%2)").arg(name).arg(m_labels[i].first->name().c_str()));
+        items.back()->setColor (m_labels[i].second);
+        map_labels[m_labels[i].first] = i;
+      }
+
+    for (Point_set::const_iterator it = m_points->point_set()->begin();
+         it != m_points->point_set()->end(); ++ it)
+      {
+        Label_handle c = m_psc->label_of(*it);
+        if (c != Label_handle())
+          items[map_labels[c]]->point_set()->insert (m_points->point_set()->point(*it));
+      }
+  }
+
+  Scene_points_with_normal_item* points_item() { return m_points; }
+
+  void reset_training_sets()
+  {
+    m_trainer->reset_inlier_sets();
+    invalidateOpenGLBuffers();
+    Q_EMIT itemChanged();
+  }
+
+  void train();
+
+  void callback()
+  {
+    invalidateOpenGLBuffers();
+    Q_EMIT itemChanged();
+  }
+  
+
+  Label_handle get_label (const char* name)
+  {
+    for (std::size_t i = 0; i < m_labels.size(); ++ i)
+      if (m_labels[i].first->name() == name)
+        return m_labels[i].first;
+    return Label_handle();
+  }
+  
+  void add_selection_to_training_set (const char* name)
+  {
+    Label_handle label = get_label (name);
+
+    for (Point_set::const_iterator it = m_points->point_set()->first_selected();
+         it != m_points->point_set()->end(); ++ it)
+      m_psc->set_label_of(*it, label);
+
+    m_trainer->set_inliers(label,
+                           boost::make_iterator_range(m_points->point_set()->first_selected(),
+                                                      m_points->point_set()->end()));
+
+    m_points->resetSelection();
+  }
+
+  void validate_selection ()
+  {
+    for (Point_set::const_iterator it = m_points->point_set()->first_selected();
+         it != m_points->point_set()->end(); ++ it)
+      {
+        Label_handle t = m_psc->label_of(*it);
+        m_trainer->set_inlier (t, *it);
+      }
+
+    m_points->resetSelection();
+  }
+
+  double& smoothing() { return m_smoothing; }
+  std::size_t& subdivisions() { return m_subdivisions; }
+  std::size_t& nb_scales() { return m_nb_scales; }
+  std::size_t& number_of_trials() { return m_nb_trials; }
+  bool features_computed() const { return (m_psc->number_of_features() != 0); }
+
+  std::vector<std::pair<Label_handle, QColor> >& labels() { return m_labels; }
+  std::size_t number_of_features() const { return m_psc->number_of_features(); }
+  Feature_handle feature(std::size_t i) { return m_psc->feature(i); }
+  void add_new_label (const char* name, const QColor& color)
+  {
+    m_labels.push_back (std::make_pair (m_psc->add_label(name),
+                                       color));
+  }
+
+  void remove_label (const char* name)
+  {
+    for (std::size_t i = 0; i < m_labels.size(); ++ i)
+      if (m_labels[i].first->name() == name)
+        {
+          m_psc->remove_label (m_labels[i].first);
+          m_labels.erase (m_labels.begin() + i);
+          break;
+        }
+    invalidateOpenGLBuffers();
+    Q_EMIT itemChanged();
+  }
+
+  void change_label_color (const char* name, const QColor& color)
+  {
+    for (std::size_t i = 0; i < m_labels.size(); ++ i)
+      if (m_labels[i].first->name() == name)
+        {
+          m_labels[i].second = color;
+          break;
+        }
+  }
+
+  template <typename ComboBox>
+ void fill_display_combo_box (ComboBox* cb, ComboBox* cb1)
+  {
+    for (std::size_t i = 0; i < m_psc->number_of_features(); ++ i)
+      {
+        std::size_t scale = m_psc->scale_of_feature(m_psc->feature(i));
+        std::ostringstream oss;
+        oss << "Feature " << m_psc->feature(i)->name() << "_" << scale;
+        cb->addItem (oss.str().c_str());
+        cb1->addItem (oss.str().c_str());
+      }
+  }
+
+  void save_config(const char* filename)
+  {
+    if (m_psc->number_of_features() == 0)
+      {
+        std::cerr << "Error: features not computed" << std::endl;
+        return;
+      }
+
+    std::ofstream f (filename);
+    m_psc->save_configuration (f);
+  }
+
+  void load_config(const char* filename)
+  {
+    if (m_psc->number_of_features() != 0)
+      m_psc->clear();
+
+    reset_indices();
+    
+    bool normals = m_points->point_set()->has_normal_map();
+    bool colors = m_points->point_set()->has_colors();
+    Point_set::Property_map<boost::uint8_t> echo_map;
+    bool echo;
+    boost::tie (echo_map, echo) = m_points->point_set()->template property_map<boost::uint8_t>("echo");
+
+    std::ifstream f (filename);
+    if (!normals && !colors && !echo)
+      m_psc->load_configuration (f);
+    else if (!normals && !colors && echo)
+      m_psc->load_configuration (f, CGAL::Default(), CGAL::Default(), echo_map);
+    else if (!normals && colors && !echo)
+      m_psc->load_configuration (f, CGAL::Default(), Color_map(m_points->point_set()));
+    else if (!normals && colors && echo)
+      m_psc->load_configuration (f, CGAL::Default(), Color_map(m_points->point_set()), echo_map);
+    else if (normals && !colors && !echo)
+      m_psc->load_configuration (f, m_points->point_set()->normal_map());
+    else if (normals && !colors && echo)
+      m_psc->load_configuration (f, m_points->point_set()->normal_map(), CGAL::Default(), echo_map);
+    else if (normals && colors && !echo)
+      m_psc->load_configuration (f, m_points->point_set()->normal_map(), Color_map(m_points->point_set()));
+    else
+      m_psc->load_configuration (f, m_points->point_set()->normal_map(), Color_map(m_points->point_set()), echo_map);
+    
+    std::vector<std::pair<Label_handle, QColor> > new_labels;
+    for (std::size_t i = 0; i < m_psc->number_of_labels(); ++ i)
+      {
+        Label_handle t = m_psc->label(i);
+        QColor color (192 + rand() % 60,
+                      192 + rand() % 60,
+                      192 + rand() % 60);
+
+        for (std::size_t j = 0; j < m_labels.size(); ++ j)
+          if (t->name() == m_labels[j].first->name())
+            {
+              color = m_labels[j].second;
+              break;
+            }
+
+        new_labels.push_back (std::make_pair (t, color));
+      }
+    m_labels.swap (new_labels);
+  }
+
+ public Q_SLOTS:
+
+  // Data
+ private:
+  
+  Scene_points_with_normal_item* m_points;
+  PSC* m_psc;
+  Trainer* m_trainer;
+
+  std::size_t m_nb_scales;
+
+  std::vector<std::pair<Label_handle, QColor> > m_labels;
+  std::size_t m_nb_trials;
+  double m_smoothing;
+  std::size_t m_subdivisions;
+  
+  int m_index_color;
+  
+  enum VAOs {
+    Edges=0,
+    ThePoints,
+    NbOfVaos = ThePoints+1
+  };
+  enum VBOs {
+    Edges_vertices = 0,
+    Points_vertices,
+    Points_colors,
+    NbOfVbos = Points_colors+1
+  };
+
+  mutable std::vector<double> positions_lines;
+  mutable std::vector<double> positions_points;
+  mutable std::vector<double> colors_points;
+  mutable std::size_t nb_points;
+  mutable std::size_t nb_lines;
+   
+
+  mutable QOpenGLShaderProgram *program;
+
+  using CGAL::Three::Scene_item::initializeBuffers;
+  void initializeBuffers(CGAL::Three::Viewer_interface *viewer) const;
+
+  void compute_normals_and_vertices() const;
+
+
+}; // end class Scene_point_set_classification_item
+
+
+
+
+#endif // POINT_SET_CLASSIFICATION_ITEM_H

Polyhedron/demo/Polyhedron/Scene_point_set_classification_item_config.h

@@ -0,0 +1,10 @@
+#ifndef SCENE_POINT_SET_CLASSIFICATION_ITEM_CONFIG_H
+#define SCENE_POINT_SET_CLASSIFICATION_ITEM_CONFIG_H
+
+#ifdef scene_point_set_classification_item_EXPORTS
+#  define SCENE_POINT_SET_CLASSIFICATION_ITEM_EXPORT Q_DECL_EXPORT
+#else
+#  define SCENE_POINT_SET_CLASSIFICATION_ITEM_EXPORT Q_DECL_IMPORT
+#endif
+
+#endif // SCENE_POINT_SET_CLASSIFICATION_ITEM_CONFIG_H

Polyhedron/demo/Polyhedron/include/read_las_point_set.h

@@ -0,0 +1,87 @@
+#ifndef READ_LAS_POINT_SET_H
+#define READ_LAS_POINT_SET_H
+
+#include <liblas/liblas.hpp>
+
+template <typename Point, typename Vector>
+bool read_las_point_set(
+  std::istream& stream, ///< input stream.
+  CGAL::Point_set_3<Point, Vector>& point_set) ///< point set
+{
+  typename CGAL::Point_set_3<Point, Vector> Point_set;
+  
+  liblas::ReaderFactory f;
+  liblas::Reader reader = f.CreateWithStream(stream);
+
+  Point_set::Property_map<unsigned char> echo_map;
+  boost::tie(echo_map, boost::tuples::ignore)
+    = point_set.template add_property_map<unsigned char>("echo", 0);
+
+  Point_set::Property_map<unsigned char> red_map;
+  boost::tie(red_map, boost::tuples::ignore)
+    = point_set.template add_property_map<unsigned char>("red", 0);
+
+  Point_set::Property_map<unsigned char> green_map;
+  boost::tie(green_map, boost::tuples::ignore)
+    = point_set.template add_property_map<unsigned char>("green", 0);
+
+  Point_set::Property_map<unsigned char> blue_map;
+  boost::tie(blue_map, boost::tuples::ignore)
+    = point_set.template add_property_map<unsigned char>("blue", 0);
+
+  unsigned int bit_short_to_char = 0;
+
+  while (reader.ReadNextPoint())
+    {
+      const liblas::Point& p = reader.GetPoint();
+      Point_set::iterator it = point_set.insert (Kernel::Point_3 (p.GetX(), p.GetY(), p.GetZ()));
+      put(echo_map, *it, p.GetReturnNumber());
+      const liblas::Color& c = p.GetColor();
+      if (bit_short_to_char == 0
+          && (c[0] > 255 || c[1] > 255 || c[2] > 255))
+        {
+          // if a value is above 255, shorts are used (that should be
+          // default but not always the case, hence this test)
+          bit_short_to_char = 8;
+          for (Point_set::iterator it2 = point_set.begin(); it2 != it; ++ it2)
+            {
+              // Correct previously added values
+              put(red_map, *it2, ((get(red_map, *it2)) >> 8));
+              put(green_map, *it2, ((get(green_map, *it2)) >> 8));
+              put(blue_map, *it2, ((get(blue_map, *it2)) >> 8));
+            }
+        }
+      put(red_map, *it, (c.GetRed() >> bit_short_to_char));
+      put(green_map, *it, (c.GetGreen() >> bit_short_to_char));
+      put(blue_map, *it, (c.GetBlue() >> bit_short_to_char));
+    }
+
+  bool remove_echo = true;
+  bool remove_colors = true;
+  
+  for (Point_set::iterator it = point_set.begin(); it != point_set.end(); ++ it)
+    {
+      if (get(echo_map, *it) != 0)
+        remove_echo = false;
+      if (get(red_map, *it) != 0 ||
+          get(green_map, *it) != 0 ||
+          get(blue_map, *it) != 0)
+        remove_colors = false;
+      if (!remove_echo && !remove_colors)
+        break;
+    }
+
+  if (remove_echo)
+    point_set.remove_property_map(echo_map);
+  if (remove_colors)
+    {
+      point_set.remove_property_map(red_map);
+      point_set.remove_property_map(green_map);
+      point_set.remove_property_map(blue_map);
+    }
+  
+  return true;
+}
+
+
+#endif // READ_LAS_POINT_SET_H

Property_map/include/CGAL/property_map.h

@@ -408,6 +408,28 @@ make_property_map(const std::vector<T>& v)
   return make_property_map(&v[0]);
 }
 
+/// \ingroup PkgProperty_map
+/// Property map that only returns the default value type
+/// \cgalModels `ReadablePropertyMap`
+template<class InputIterator, class ValueType>
+struct Default_property_map{
+  const ValueType default_value;
+  
+  typedef typename InputIterator::value_type key_type;
+  typedef boost::readable_property_map_tag category;
+
+  Default_property_map(const ValueType& default_value = ValueType()) : default_value (default_value) { }
+  
+  /// Free function to use a get the value from an iterator using Input_iterator_property_map.
+  inline friend ValueType
+  get (const Default_property_map&, const key_type&){ return ValueType(); }
+};
+
+
+/// \endcond
+
 } // namespace CGAL
 
+
+
 #endif // CGAL_POINT_SET_PROPERTY_MAP_H

Stream_support/include/CGAL/IO/io.h

@@ -174,6 +174,59 @@ public:
     return is; 
   }
 };
+
+template <>
+class Input_rep<float> {
+    float& t;
+public:
+  //! initialize with a reference to \a t.
+  Input_rep( float& tt) : t(tt) {}
+
+  std::istream& operator()( std::istream& is) const 
+  {
+    typedef std::istream istream;
+    typedef istream::char_type char_type;
+    typedef istream::int_type int_type;
+    typedef istream::traits_type traits_type;
+
+    std::string buffer;
+    buffer.reserve(32);
+
+    char_type c;
+    do {
+      const int_type i = is.get();
+      if(i == traits_type::eof()) {
+	return is;
+      }
+      c = static_cast<char_type>(i);
+    }while (std::isspace(c));
+    if(c == '-'){
+      buffer += '-';
+    } else if(c != '+'){
+      is.unget();
+    }
+    do {
+      const int_type i = is.get();
+      if(i == traits_type::eof()) {
+	is.clear(is.rdstate() & ~std::ios_base::failbit);
+	break;
+      }
+      c = static_cast<char_type>(i);
+      if(std::isdigit(c) || (c =='.') || (c =='E') || (c =='e') || (c =='+') || (c =='-')){
+        buffer += c;
+      }else{
+	is.unget();
+	break;
+      }
+    }while(true);
+
+    if(sscanf(buffer.c_str(), "%f", &t) != 1) {
+      // if a 'buffer' does not contain a double, set the fail bit.
+      is.setstate(std::ios_base::failbit);
+    }
+    return is; 
+  }
+};
 #endif
 
 /*! \relates Input_rep

Surface_mesh_segmentation/include/CGAL/internal/auxiliary/graph.h

@@ -1751,4 +1751,6 @@ inline Graph::termtype Graph::what_segment(node_id i)
   return SINK;
 }
 
+#undef last_node
+
 #endif