From 77a590d33878c5fd8e4c8050a8e59cd2d9d20441 Mon Sep 17 00:00:00 2001 From: Simon Giraudot Date: Thu, 16 Mar 2017 14:16:44 +0100 Subject: [PATCH] Global change of API --- .../doc/Classification/PackageDescription.txt | 17 +- .../doc/Classification/examples.txt | 4 +- .../examples/Classification/CMakeLists.txt | 4 +- .../Classification/example_classification.cpp | 202 +++ .../Classification/example_classifier.cpp | 194 --- .../Classification/example_feature.cpp | 66 +- .../example_generation_and_training.cpp | 150 +++ .../example_point_set_classifier.cpp | 184 --- Classification/include/CGAL/Classification.h | 45 + .../include/CGAL/Classification/Evaluation.h | 181 +++ .../Feature/Distance_to_plane.h | 10 +- .../Classification/Feature/Echo_scatter.h | 5 +- .../CGAL/Classification/Feature/Effect.h | 42 - .../CGAL/Classification/Feature/Eigen.h | 20 +- .../CGAL/Classification/Feature/Elevation.h | 6 +- .../include/CGAL/Classification/Feature/Hsv.h | 7 +- .../Feature/Vertical_dispersion.h | 5 +- .../CGAL/Classification/Feature/Verticality.h | 11 +- .../CGAL/Classification/Feature_base.h | 53 +- .../include/CGAL/Classification/Feature_set.h | 182 +++ .../include/CGAL/Classification/Label.h | 45 +- .../include/CGAL/Classification/Label_set.h | 113 ++ .../Point_set_feature_generator.h | 752 +++++++++++ .../Sum_of_weighted_features_predicate.h | 731 +++++++++++ .../include/CGAL/Classification/classify.h | 438 +++++++ Classification/include/CGAL/Classifier.h | 1087 --------------- .../include/CGAL/Point_set_classifier.h | 1166 ----------------- .../Classification/Classification_plugin.cpp | 43 +- .../Classification/Item_classification_base.h | 107 +- .../Point_set_item_classification.cpp | 190 +-- .../Point_set_item_classification.h | 168 +-- 31 files changed, 3122 insertions(+), 3106 deletions(-) create mode 100644 Classification/examples/Classification/example_classification.cpp delete mode 100644 Classification/examples/Classification/example_classifier.cpp create mode 100644 Classification/examples/Classification/example_generation_and_training.cpp delete mode 100644 Classification/examples/Classification/example_point_set_classifier.cpp create mode 100644 Classification/include/CGAL/Classification.h create mode 100644 Classification/include/CGAL/Classification/Evaluation.h delete mode 100644 Classification/include/CGAL/Classification/Feature/Effect.h create mode 100644 Classification/include/CGAL/Classification/Feature_set.h create mode 100644 Classification/include/CGAL/Classification/Label_set.h create mode 100644 Classification/include/CGAL/Classification/Point_set_feature_generator.h create mode 100644 Classification/include/CGAL/Classification/Sum_of_weighted_features_predicate.h create mode 100644 Classification/include/CGAL/Classification/classify.h delete mode 100644 Classification/include/CGAL/Classifier.h delete mode 100644 Classification/include/CGAL/Point_set_classifier.h diff --git a/Classification/doc/Classification/PackageDescription.txt b/Classification/doc/Classification/PackageDescription.txt index 0e8c0f425ad..a47eb9771d5 100644 --- a/Classification/doc/Classification/PackageDescription.txt +++ b/Classification/doc/Classification/PackageDescription.txt @@ -35,29 +35,38 @@ ## Concepts ## +- `CGAL::ClassificationPredicate` - `CGAL::NeighborQuery` -## Classification ## +## Main Functions ## -- `CGAL::Classifier` -- `CGAL::Point_set_classifier` +- `CGAL::Classification::classify()` +- `CGAL::Classification::classify_with_local_smoothing()` +- `CGAL::Classification::classify_with_graphcut()` + +## Classification Predicate ## + +- `CGAL::Classification::Sum_of_weighted_features_predicate` ## Data Structures ## +- `CGAL::Classification::Point_set_feature_generator` - `CGAL::Classification::Point_set_neighborhood` - `CGAL::Classification::Local_eigen_analysis` - `CGAL::Classification::Planimetric_grid` -- `CGAL::Classification::Trainer` + ## Label ## - `CGAL::Classification::Label` - `CGAL::Classification::Label_handle` +- `CGAL::Classification::Label_set` ## Feature ## - `CGAL::Classification::Feature` - `CGAL::Classification::Feature_handle` +- `CGAL::Classification::Feature_set` ## Predefined Features ## diff --git a/Classification/doc/Classification/examples.txt b/Classification/doc/Classification/examples.txt index 1307af0c1f1..0b39c016e77 100644 --- a/Classification/doc/Classification/examples.txt +++ b/Classification/doc/Classification/examples.txt @@ -1,5 +1,5 @@ /*! -\example Classification/example_classifier.cpp +\example Classification/example_classification.cpp \example Classification/example_feature.cpp -\example Classification/example_point_set_classifier.cpp +\example Classification/example_generation_and_training.cpp */ diff --git a/Classification/examples/Classification/CMakeLists.txt b/Classification/examples/Classification/CMakeLists.txt index c20e05a4ca0..ee1b7185eb5 100644 --- a/Classification/examples/Classification/CMakeLists.txt +++ b/Classification/examples/Classification/CMakeLists.txt @@ -50,6 +50,6 @@ include_directories( BEFORE ../../include ) 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_classification.cpp" ) +create_single_source_cgal_program( "example_generation_and_training.cpp" ) create_single_source_cgal_program( "example_feature.cpp" ) diff --git a/Classification/examples/Classification/example_classification.cpp b/Classification/examples/Classification/example_classification.cpp new file mode 100644 index 00000000000..80d13429b84 --- /dev/null +++ b/Classification/examples/Classification/example_classification.cpp @@ -0,0 +1,202 @@ +#include +#include +#include +#include + +#include +#include +#include + +#include + +#include + +typedef CGAL::Simple_cartesian Kernel; +typedef Kernel::Point_3 Point; +typedef Kernel::Iso_cuboid_3 Iso_cuboid_3; +typedef std::vector Point_range; +typedef CGAL::Identity_property_map Pmap; + +namespace Classif = CGAL::Classification; + +typedef Classif::Sum_of_weighted_features_predicate Classification_predicate; + +typedef Classif::Planimetric_grid Planimetric_grid; +typedef Classif::Point_set_neighborhood Neighborhood; +typedef Classif::Local_eigen_analysis Local_eigen_analysis; + +typedef Classif::Label_handle Label_handle; +typedef Classif::Feature_handle Feature_handle; +typedef Classif::Label_set Label_set; +typedef Classif::Feature_set Feature_set; + +typedef Classif::Feature::Distance_to_plane Distance_to_plane; +typedef Classif::Feature::Linearity Linearity; +typedef Classif::Feature::Omnivariance Omnivariance; +typedef Classif::Feature::Planarity Planarity; +typedef Classif::Feature::Surface_variation Surface_variation; +typedef Classif::Feature::Elevation Elevation; +typedef Classif::Feature::Vertical_dispersion 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 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)); + + //! [Analysis] + /////////////////////////////////////////////////////////////////// + + /////////////////////////////////////////////////////////////////// + //! [Features] + + std::cerr << "Computing features" << std::endl; + Feature_set features; + Feature_handle d2p = features.add (pts, Pmap(), eigen); + Feature_handle lin = features.add (pts, eigen); + Feature_handle omni = features.add (pts, eigen); + Feature_handle plan = features.add (pts, eigen); + Feature_handle surf = features.add (pts, eigen); + Feature_handle disp = features.add (pts, Pmap(), grid, + grid_resolution, + radius_neighbors); + Feature_handle elev = features.add (pts, Pmap(), grid, + grid_resolution, + radius_dtm); + + //! [Features] + /////////////////////////////////////////////////////////////////// + + /////////////////////////////////////////////////////////////////// + //! [Labels] + + Label_set labels; + Label_handle ground = labels.add ("ground"); + Label_handle vege = labels.add ("vegetation"); + Label_handle roof = labels.add ("roof"); + + std::cerr << "Setting weights" << std::endl; + Classification_predicate predicate (labels, features); + predicate.set_weight (d2p, 6.75e-2); + predicate.set_weight (lin, 1.19); + predicate.set_weight (omni, 1.34e-1); + predicate.set_weight (plan, 7.32e-1); + predicate.set_weight (surf, 1.36e-1); + predicate.set_weight (disp, 5.45e-1); + predicate.set_weight (elev, 1.47e1); + + std::cerr << "Setting effects" << std::endl; + predicate.set_effect (ground, d2p, Classification_predicate::NEUTRAL); + predicate.set_effect (ground, lin, Classification_predicate::PENALIZING); + predicate.set_effect (ground, omni, Classification_predicate::NEUTRAL); + predicate.set_effect (ground, plan, Classification_predicate::FAVORING); + predicate.set_effect (ground, surf, Classification_predicate::PENALIZING); + predicate.set_effect (ground, disp, Classification_predicate::NEUTRAL); + predicate.set_effect (ground, elev, Classification_predicate::PENALIZING); + + predicate.set_effect (vege, d2p, Classification_predicate::FAVORING); + predicate.set_effect (vege, lin, Classification_predicate::NEUTRAL); + predicate.set_effect (vege, omni, Classification_predicate::FAVORING); + predicate.set_effect (vege, plan, Classification_predicate::NEUTRAL); + predicate.set_effect (vege, surf, Classification_predicate::NEUTRAL); + predicate.set_effect (vege, disp, Classification_predicate::FAVORING); + predicate.set_effect (vege, elev, Classification_predicate::NEUTRAL); + + predicate.set_effect (roof, d2p, Classification_predicate::NEUTRAL); + predicate.set_effect (roof, lin, Classification_predicate::PENALIZING); + predicate.set_effect (roof, omni, Classification_predicate::FAVORING); + predicate.set_effect (roof, plan, Classification_predicate::FAVORING); + predicate.set_effect (roof, surf, Classification_predicate::PENALIZING); + predicate.set_effect (roof, disp, Classification_predicate::NEUTRAL); + predicate.set_effect (roof, elev, Classification_predicate::FAVORING); + + //! [Labels] + /////////////////////////////////////////////////////////////////// + + // Run classification + std::cerr << "Classifying" << std::endl; + std::vector label_indices; + + CGAL::Real_timer t; + t.start(); + Classif::classify (pts, labels, predicate, label_indices); + t.stop(); + std::cerr << "Raw classification performed in " << t.time() << " second(s)" << std::endl; + t.reset(); + + t.start(); + Classif::classify_with_local_smoothing + (pts, Pmap(), labels, predicate, + neighborhood.range_neighbor_query(radius_neighbors), + label_indices); + t.stop(); + std::cerr << "Classification with local smoothing performed in " << t.time() << " second(s)" << std::endl; + t.reset(); + + t.start(); + Classif::classify_with_graphcut + (pts, Pmap(), Pmap(), labels, predicate, + neighborhood.k_neighbor_query(12), + 0.2, 1, label_indices); + 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 = labels[label_indices[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; +} diff --git a/Classification/examples/Classification/example_classifier.cpp b/Classification/examples/Classification/example_classifier.cpp deleted file mode 100644 index bc7f780c4c0..00000000000 --- a/Classification/examples/Classification/example_classifier.cpp +++ /dev/null @@ -1,194 +0,0 @@ -#include -#include -#include -#include - -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#include - -#include - -typedef CGAL::Simple_cartesian Kernel; -typedef Kernel::Point_3 Point; -typedef Kernel::Iso_cuboid_3 Iso_cuboid_3; -typedef std::vector Point_range; -typedef CGAL::Identity_property_map Pmap; - -typedef CGAL::Classifier Classifier; - -typedef CGAL::Classification::Planimetric_grid Planimetric_grid; -typedef CGAL::Classification::Point_set_neighborhood Neighborhood; -typedef CGAL::Classification::Local_eigen_analysis Local_eigen_analysis; - -typedef CGAL::Classification::Label_handle Label_handle; -typedef CGAL::Classification::Feature_handle Feature_handle; - -typedef CGAL::Classification::Feature::Distance_to_plane Distance_to_plane; -typedef CGAL::Classification::Feature::Linearity Linearity; -typedef CGAL::Classification::Feature::Omnivariance Omnivariance; -typedef CGAL::Classification::Feature::Planarity Planarity; -typedef CGAL::Classification::Feature::Surface_variation Surface_variation; -typedef CGAL::Classification::Feature::Elevation Elevation; -typedef CGAL::Classification::Feature::Vertical_dispersion 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 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 (Pmap(), eigen); - Feature_handle lin = classifier.add_feature (eigen); - Feature_handle omni = classifier.add_feature (eigen); - Feature_handle plan = classifier.add_feature (eigen); - Feature_handle surf = classifier.add_feature (eigen); - Feature_handle disp = classifier.add_feature (Pmap(), grid, - grid_resolution, - radius_neighbors); - Feature_handle elev = classifier.add_feature (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; -} diff --git a/Classification/examples/Classification/example_feature.cpp b/Classification/examples/Classification/example_feature.cpp index c105684ace1..a45724fc2b5 100644 --- a/Classification/examples/Classification/example_feature.cpp +++ b/Classification/examples/Classification/example_feature.cpp @@ -4,12 +4,7 @@ #include #include -#include -#include -#include -#include -#include - +#include #include typedef CGAL::Simple_cartesian Kernel; @@ -18,15 +13,19 @@ typedef Kernel::Iso_cuboid_3 Iso_cuboid_3; typedef std::vector Point_range; typedef CGAL::Identity_property_map Pmap; -typedef CGAL::Classifier Classifier; +namespace Classif = CGAL::Classification; -typedef CGAL::Classification::Point_set_neighborhood Neighborhood; -typedef CGAL::Classification::Local_eigen_analysis Local_eigen_analysis; +typedef Classif::Sum_of_weighted_features_predicate Classification_predicate; -typedef CGAL::Classification::Label_handle Label_handle; -typedef CGAL::Classification::Feature_handle Feature_handle; +typedef Classif::Point_set_neighborhood Neighborhood; +typedef Classif::Local_eigen_analysis Local_eigen_analysis; -typedef CGAL::Classification::Feature::Sphericity Sphericity; +typedef Classif::Label_handle Label_handle; +typedef Classif::Feature_handle Feature_handle; +typedef Classif::Label_set Label_set; +typedef Classif::Feature_set Feature_set; + +typedef Classif::Feature::Sphericity Sphericity; // User-defined feature that identifies a specific area of the 3D @@ -37,10 +36,12 @@ 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 + My_feature (const Point_range& range, double xmin, double xmax, double ymin, double ymax) : range (range), xmin(xmin), xmax(xmax), ymin(ymin), ymax(ymax) - { } + { + this->set_name ("my_feature"); + } double value (std::size_t pt_index) { @@ -70,29 +71,36 @@ int main (int argc, char** argv) Neighborhood neighborhood (pts, Pmap()); Local_eigen_analysis eigen (pts, Pmap(), neighborhood.k_neighbor_query(6)); - Classifier classifier (pts, Pmap()); - + Label_set labels; + Label_handle a = labels.add ("label_A"); + Label_handle b = labels.add ("label_B"); + std::cerr << "Computing features" << std::endl; - Feature_handle sphericity = classifier.add_feature (eigen); + Feature_set features; + + Feature_handle sphericity = features.add (pts, 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 (-20., 20., -15., 15.); - + Feature_handle my_feature = features.add (pts, -20., 20., -15., 15.); + + Classification_predicate predicate (labels, features); + std::cerr << "Setting weights" << std::endl; - sphericity->set_weight(0.5); - my_feature->set_weight(0.25); + predicate.set_weight(sphericity, 0.5); + predicate.set_weight(my_feature, 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); + predicate.set_effect (a, sphericity, Classification_predicate::FAVORING); + predicate.set_effect (a, my_feature, Classification_predicate::FAVORING); + predicate.set_effect (b, sphericity, Classification_predicate::PENALIZING); + predicate.set_effect (b, my_feature, Classification_predicate::PENALIZING); - 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::vector label_indices; + Classif::classify_with_graphcut + (pts, Pmap(), Pmap(), labels, predicate, + neighborhood.k_neighbor_query(12), + 0.5, 1, label_indices); std::cerr << "All done" << std::endl; return EXIT_SUCCESS; diff --git a/Classification/examples/Classification/example_generation_and_training.cpp b/Classification/examples/Classification/example_generation_and_training.cpp new file mode 100644 index 00000000000..8824a8e0d64 --- /dev/null +++ b/Classification/examples/Classification/example_generation_and_training.cpp @@ -0,0 +1,150 @@ +#include +#include +#include +#include + +//#define CGAL_CLASSIFICATION_VERBOSE + +#include +#include +#include + +#include + +typedef CGAL::Simple_cartesian Kernel; +typedef Kernel::Point_3 Point; +typedef Kernel::Iso_cuboid_3 Iso_cuboid_3; +typedef std::vector Point_range; +typedef CGAL::Identity_property_map Pmap; + +namespace Classif = CGAL::Classification; + +typedef Classif::Label_handle Label_handle; +typedef Classif::Feature_handle Feature_handle; +typedef Classif::Label_set Label_set; +typedef Classif::Feature_set Feature_set; + +typedef Classif::Sum_of_weighted_features_predicate Classification_predicate; + +typedef Classif::Point_set_feature_generator Feature_generator; + +/* + This interpreter is used to read a PLY input that contains training + attributes (with the PLY "label" property). +*/ +class My_ply_interpreter +{ + std::vector& points; + std::vector& labels; + +public: + My_ply_interpreter (std::vector& points, + std::vector& 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 ("x") + && reader.does_tag_exist ("y") + && reader.does_tag_exist ("z") + && reader.does_tag_exist ("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(std::size_t(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 pts; + std::vector ground_truth; + + std::cerr << "Reading input" << std::endl; + My_ply_interpreter interpreter (pts, ground_truth); + if (!in + || !(CGAL::read_ply_custom_points (in, interpreter, Kernel()))) + { + std::cerr << "Error: cannot read " << filename << std::endl; + return EXIT_FAILURE; + } + + Feature_set features; + + std::cerr << "Generating features" << std::endl; + CGAL::Real_timer t; + t.start(); + Feature_generator generator (features, 5, // using 5 scales + pts, Pmap()); + t.stop(); + std::cerr << "Done in " << t.time() << " second(s)" << std::endl; + + // Add types + Label_set labels; + Label_handle ground = labels.add ("ground"); + Label_handle vege = labels.add ("vegetation"); + Label_handle roof = labels.add ("roof"); + Label_handle facade = labels.add ("facade"); + + Classification_predicate predicate (labels, features); + + std::cerr << "Training" << std::endl; + t.reset(); + t.start(); + predicate.train (ground_truth, 400); + t.stop(); + std::cerr << "Done in " << t.time() << " second(s)" << std::endl; + + t.reset(); + t.start(); + std::vector label_indices; + Classif::classify_with_graphcut + (pts, Pmap(), Pmap(), labels, predicate, + generator.neighborhood().k_neighbor_query(12), + 0.2, 10, label_indices); + t.stop(); + std::cerr << "Classification with graphcut done in " << t.time() << " second(s)" << std::endl; + + std::cerr << "Precision, recall, F1 scores and IoU:" << std::endl; + Classif::Evaluation eval (labels, ground_truth, label_indices); + + for (std::size_t i = 0; i < labels.size(); ++ i) + { + std::cerr << " * " << labels[i]->name() << ": " + << eval.precision(labels[i]) << " ; " + << eval.recall(labels[i]) << " ; " + << eval.f1_score(labels[i]) << " ; " + << eval.intersection_over_union(labels[i]) << std::endl; + } + + std::cerr << "Accuracy = " << eval.accuracy() << std::endl + << "Mean F1 score = " << eval.mean_f1_score() << std::endl + << "Mean IoU = " << eval.mean_intersection_over_union() << std::endl; + + + /// Save the configuration to be able to reload it later + std::ofstream fconfig ("config.xml"); + predicate.save_configuration (fconfig); + fconfig.close(); + + std::cerr << "All done" << std::endl; + + return EXIT_SUCCESS; +} diff --git a/Classification/examples/Classification/example_point_set_classifier.cpp b/Classification/examples/Classification/example_point_set_classifier.cpp deleted file mode 100644 index 391c45a62ac..00000000000 --- a/Classification/examples/Classification/example_point_set_classifier.cpp +++ /dev/null @@ -1,184 +0,0 @@ -#include -#include -#include -#include - -//#define CGAL_CLASSIFICATION_VERBOSE - -#include -#include -#include -#include - -#include - -typedef CGAL::Simple_cartesian Kernel; -typedef Kernel::Point_3 Point; -typedef Kernel::Iso_cuboid_3 Iso_cuboid_3; -typedef std::vector Point_range; -typedef CGAL::Identity_property_map Pmap; - -typedef CGAL::Point_set_classifier Point_set_classifier; -typedef CGAL::Classification::Trainer 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& points; - std::vector& labels; - -public: - My_ply_interpreter (std::vector& points, - std::vector& 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 ("x") - && reader.does_tag_exist ("y") - && reader.does_tag_exist ("z") - && reader.does_tag_exist ("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 pts; - std::vector 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; -} diff --git a/Classification/include/CGAL/Classification.h b/Classification/include/CGAL/Classification.h new file mode 100644 index 00000000000..a95e3147e4d --- /dev/null +++ b/Classification/include/CGAL/Classification.h @@ -0,0 +1,45 @@ +// 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_H +#define CGAL_CLASSIFICATION_H + +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include + +#endif // CGAL_CLASSIFICATION_H diff --git a/Classification/include/CGAL/Classification/Evaluation.h b/Classification/include/CGAL/Classification/Evaluation.h new file mode 100644 index 00000000000..7056ce92f37 --- /dev/null +++ b/Classification/include/CGAL/Classification/Evaluation.h @@ -0,0 +1,181 @@ +// 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_EVALUATION_H +#define CGAL_CLASSIFICATION_EVALUATION_H_ + +#include +#include + +namespace CGAL { + +namespace Classification { + +class Evaluation +{ + mutable std::map m_map_labels; + + std::vector m_precision; + std::vector m_recall; + std::vector m_iou; // intersection over union + double m_accuracy; + double m_mean_iou; + double m_mean_f1; + +public: + + Evaluation (const Label_set& labels, + const std::vector& ground_truth, + const std::vector& result) + : m_precision (labels.size()), + m_recall (labels.size()), + m_iou (labels.size()) + { + for (std::size_t i = 0; i < labels.size(); ++ i) + m_map_labels[labels[i]] = i; + + std::vector true_positives (labels.size()); + std::vector false_positives (labels.size()); + std::vector false_negatives (labels.size()); + + std::size_t sum_true_positives = 0; + std::size_t total = 0; + + for (std::size_t j = 0; j < ground_truth.size(); ++ j) + { + std::size_t gt = ground_truth[j]; + std::size_t res = result[j]; + if (gt == std::size_t(-1) || res == std::size_t(-1)) + continue; + ++ total; + if (gt == res) + { + ++ true_positives[gt]; + ++ sum_true_positives; + continue; + } + ++ false_positives[res]; + ++ false_negatives[gt]; + } + + m_mean_iou = 0.; + m_mean_f1 = 0.; + + for (std::size_t j = 0; j < labels.size(); ++ j) + { + m_precision[j] = true_positives[j] / double(true_positives[j] + false_positives[j]); + m_recall[j] = true_positives[j] / double(true_positives[j] + false_negatives[j]); + m_iou[j] = true_positives[j] / double(true_positives[j] + false_positives[j] + false_negatives[j]); + + m_mean_iou += m_iou[j]; + m_mean_f1 += 2. * (m_precision[j] * m_recall[j]) + / (m_precision[j] + m_recall[j]); + } + + m_mean_iou /= labels.size(); + m_mean_f1 /= labels.size(); + m_accuracy = sum_true_positives / double(total); + } + + + /*! + + \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 + { + return m_precision[m_map_labels[label]]; + } + + /*! + + \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 + { + return m_recall[m_map_labels[label]]; + } + + /*! + + \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 = m_map_labels[label]; + 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 + { + return m_iou[m_map_labels[label]]; + } + + /*! + \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; } + +}; + + +} // namespace Classification + +} // namespace CGAL + +#endif // CGAL_CLASSIFICATION_EVALUATION_H_ diff --git a/Classification/include/CGAL/Classification/Feature/Distance_to_plane.h b/Classification/include/CGAL/Classification/Feature/Distance_to_plane.h index 52e53e720ad..a5bd18aea00 100644 --- a/Classification/include/CGAL/Classification/Feature/Distance_to_plane.h +++ b/Classification/include/CGAL/Classification/Feature/Distance_to_plane.h @@ -23,8 +23,6 @@ #include -#include - namespace CGAL { namespace Classification { @@ -78,7 +76,8 @@ public: : input(input), point_map(point_map), eigen(eigen) #endif { - this->set_weight(1.); + this->set_name ("distance_to_plane"); + #ifndef CGAL_CLASSIFICATION_PRECOMPUTE_FEATURES std::vector distance_to_plane_feature; #endif @@ -86,9 +85,6 @@ public: 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 @@ -101,8 +97,6 @@ public: (get(point_map, *(input.begin()+pt_index)), eigen.plane(pt_index))); #endif } - - virtual std::string name() { return "distance_to_plane"; } /// \endcond }; diff --git a/Classification/include/CGAL/Classification/Feature/Echo_scatter.h b/Classification/include/CGAL/Classification/Feature/Echo_scatter.h index 116f828fc33..748946c00ee 100644 --- a/Classification/include/CGAL/Classification/Feature/Echo_scatter.h +++ b/Classification/include/CGAL/Classification/Feature/Echo_scatter.h @@ -73,7 +73,7 @@ public: const double grid_resolution, double radius_neighbors = 1.) { - this->set_weight(1.); + this->set_name ("echo_scatter"); 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++) @@ -130,7 +130,6 @@ public: 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 @@ -138,8 +137,6 @@ public: { return echo_scatter[pt_index]; } - - virtual std::string name() { return "echo_scatter"; } /// \endcond }; diff --git a/Classification/include/CGAL/Classification/Feature/Effect.h b/Classification/include/CGAL/Classification/Feature/Effect.h deleted file mode 100644 index 3a158656f0f..00000000000 --- a/Classification/include/CGAL/Classification/Feature/Effect.h +++ /dev/null @@ -1,42 +0,0 @@ -// 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 diff --git a/Classification/include/CGAL/Classification/Feature/Eigen.h b/Classification/include/CGAL/Classification/Feature/Eigen.h index d61391e2688..f1a36ec6a3f 100644 --- a/Classification/include/CGAL/Classification/Feature/Eigen.h +++ b/Classification/include/CGAL/Classification/Feature/Eigen.h @@ -55,15 +55,12 @@ public: virtual void init (const PointRange& input, const Local_eigen_analysis& eigen) { - this->set_weight(1.); #ifndef CGAL_CLASSIFICATION_PRECOMPUTE_FEATURES std::vector 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; @@ -75,7 +72,6 @@ public: return get_value(eigen, pt_index); #endif } - virtual std::string name() { return "eigen_feature"; } }; @@ -115,6 +111,7 @@ public: Linearity (const PointRange& input, const Local_eigen_analysis& eigen) : Base (input, eigen) { + this->set_name("linearity"); this->init(input, eigen); } @@ -127,7 +124,6 @@ public: else return ((ev[2] - ev[1]) / ev[2]); } - virtual std::string name() { return "linearity"; } /// \endcond }; @@ -168,6 +164,7 @@ public: const Local_eigen_analysis& eigen) : Base(input, eigen) { + this->set_name("planarity"); this->init(input, eigen); } /// \cond SKIP_IN_MANUAL @@ -179,7 +176,6 @@ public: else return ((ev[1] - ev[0]) / ev[2]); } - virtual std::string name() { return "planarity"; } /// \endcond }; @@ -221,6 +217,7 @@ public: const Local_eigen_analysis& eigen) : Base(input, eigen) { + this->set_name("sphericity"); this->init(input, eigen); } /// \cond SKIP_IN_MANUAL @@ -232,7 +229,6 @@ public: else return (ev[0] / ev[2]); } - virtual std::string name() { return "sphericity"; } /// \endcond }; @@ -273,6 +269,7 @@ public: const Local_eigen_analysis& eigen) : Base(input, eigen) { + this->set_name("omnivariance"); this->init(input, eigen); } /// \cond SKIP_IN_MANUAL @@ -281,7 +278,6 @@ public: 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 }; @@ -322,6 +318,7 @@ public: const Local_eigen_analysis& eigen) : Base(input, eigen) { + this->set_name("anisotropy"); this->init(input, eigen); } /// \cond SKIP_IN_MANUAL @@ -333,7 +330,6 @@ public: else return ((ev[2] - ev[0]) / ev[2]); } - virtual std::string name() { return "anisotropy"; } /// \endcond }; @@ -374,6 +370,7 @@ public: const Local_eigen_analysis& eigen) : Base(input, eigen) { + this->set_name("eigentropy"); this->init(input, eigen); } /// \cond SKIP_IN_MANUAL @@ -389,7 +386,6 @@ public: - ev[1] * std::log(ev[1]) - ev[2] * std::log(ev[2])); } - virtual std::string name() { return "eigentropy"; } /// \endcond }; @@ -431,6 +427,7 @@ public: const Local_eigen_analysis& eigen) : Base(input, eigen) { + this->set_name("sum_eigen"); this->init(input, eigen); } /// \cond SKIP_IN_MANUAL @@ -438,7 +435,6 @@ public: { return eigen.sum_of_eigenvalues(i); } - virtual std::string name() { return "sum_eigen"; } /// \endcond }; @@ -480,6 +476,7 @@ public: const Local_eigen_analysis& eigen) : Base(input, eigen) { + this->set_name("surface_variation"); this->init(input, eigen); } /// \cond SKIP_IN_MANUAL @@ -491,7 +488,6 @@ public: else return (ev[0] / (ev[0] + ev[1] + ev[2])); } - virtual std::string name() { return "surface_variation"; } /// \endcond }; diff --git a/Classification/include/CGAL/Classification/Feature/Elevation.h b/Classification/include/CGAL/Classification/Feature/Elevation.h index ec3616f9db0..0645a558a2b 100644 --- a/Classification/include/CGAL/Classification/Feature/Elevation.h +++ b/Classification/include/CGAL/Classification/Feature/Elevation.h @@ -76,7 +76,7 @@ public: const double grid_resolution, double radius_dtm = -1.) { - this->set_weight(1.); + this->set_name ("elevation"); if (radius_dtm < 0.) radius_dtm = 100. * grid_resolution; @@ -147,7 +147,6 @@ public: 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 @@ -155,8 +154,7 @@ public: { return elevation_feature[pt_index]; } - - virtual std::string name() { return "elevation"; } + /// \endcond }; diff --git a/Classification/include/CGAL/Classification/Feature/Hsv.h b/Classification/include/CGAL/Classification/Feature/Hsv.h index 7f23751aab0..68eb821d931 100644 --- a/Classification/include/CGAL/Classification/Feature/Hsv.h +++ b/Classification/include/CGAL/Classification/Feature/Hsv.h @@ -81,7 +81,6 @@ class Hsv : public Feature_base double m_mean; double m_sd; #endif - std::string m_name; public: @@ -104,8 +103,6 @@ public: : input(input), color_map(color_map), m_channel(channel), m_mean(mean), m_sd(sd) #endif { - this->set_weight(1.); - #ifndef CGAL_CLASSIFICATION_PRECOMPUTE_FEATURES std::vector color_feature; #endif @@ -114,14 +111,13 @@ public: 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(); + this->set_name (oss.str()); } /// \cond SKIP_IN_MANUAL @@ -135,7 +131,6 @@ public: #endif } - virtual std::string name() { return m_name; } /// \endcond }; diff --git a/Classification/include/CGAL/Classification/Feature/Vertical_dispersion.h b/Classification/include/CGAL/Classification/Feature/Vertical_dispersion.h index 62642edd3b2..7f4a151f28c 100644 --- a/Classification/include/CGAL/Classification/Feature/Vertical_dispersion.h +++ b/Classification/include/CGAL/Classification/Feature/Vertical_dispersion.h @@ -72,7 +72,7 @@ public: const double grid_resolution, double radius_neighbors = -1.) { - this->set_weight(1.); + this->set_name ("vertical_dispersion"); if (radius_neighbors < 0.) radius_neighbors = 5. * grid_resolution; @@ -142,15 +142,12 @@ public: 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 }; diff --git a/Classification/include/CGAL/Classification/Feature/Verticality.h b/Classification/include/CGAL/Classification/Feature/Verticality.h index 69474ad2a76..d8636d50835 100644 --- a/Classification/include/CGAL/Classification/Feature/Verticality.h +++ b/Classification/include/CGAL/Classification/Feature/Verticality.h @@ -67,7 +67,7 @@ public: Verticality (const PointRange& input, const Local_eigen_analysis& eigen) { - this->set_weight(1.); + this->set_name ("verticality"); typename Geom_traits::Vector_3 vertical (0., 0., 1.); for (std::size_t i = 0; i < input.size(); i++) @@ -77,8 +77,6 @@ public: verticality_feature.push_back (1. - CGAL::abs(normal * vertical)); } - this->compute_mean_max (verticality_feature, this->mean, this->max); - // max *= 2; } /*! @@ -94,7 +92,7 @@ public: Verticality (const PointRange& input, VectorMap normal_map) { - this->set_weight(1.); + this->set_name ("verticality"); typename Geom_traits::Vector_3 vertical (0., 0., 1.); for (std::size_t i = 0; i < input.size(); i++) @@ -103,9 +101,6 @@ public: 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; } @@ -114,8 +109,6 @@ public: { return verticality_feature[pt_index]; } - - virtual std::string name() { return "verticality"; } /// \endcond }; diff --git a/Classification/include/CGAL/Classification/Feature_base.h b/Classification/include/CGAL/Classification/Feature_base.h index 93b61dab7b1..65c402e4802 100644 --- a/Classification/include/CGAL/Classification/Feature_base.h +++ b/Classification/include/CGAL/Classification/Feature_base.h @@ -45,32 +45,26 @@ namespace Classification { class Feature_base { - double m_weight; + std::string m_name; + public: /// \cond SKIP_IN_MANUAL - double mean; - double max; - + Feature_base() : m_name ("abstract_feature") { } virtual ~Feature_base() { } /// \endcond /*! - \brief Returns the weight of the feature. + \brief Returns the name of the feature (initialized to + `abstract_feature` for `Feature_base`). */ - double weight() const { return m_weight; } + const std::string& name() const { return m_name; } /*! - \brief Sets the weight of the feature (`weight` must be positive). + \brief Changes the name of the feature. */ - 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"; } - + void set_name (const std::string& name) { m_name = name; } + /*! \brief Returns the value taken by the feature for at the item at position `index`. This method must be implemented by inherited @@ -78,35 +72,6 @@ public: */ 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& vect, double& mean, double& max) - { - mean = 0.; - max = -std::numeric_limits::max(); - double min = std::numeric_limits::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 }; diff --git a/Classification/include/CGAL/Classification/Feature_set.h b/Classification/include/CGAL/Classification/Feature_set.h new file mode 100644 index 00000000000..99d8d35b08e --- /dev/null +++ b/Classification/include/CGAL/Classification/Feature_set.h @@ -0,0 +1,182 @@ +// 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_SET_H +#define CGAL_CLASSIFICATION_FEATURE_SET_H + +#include + +#ifdef CGAL_LINKED_WITH_TBB +#include +#endif // CGAL_LINKED_WITH_TBB + +#include + +namespace CGAL { + +namespace Classification { + +class Feature_set +{ + typedef std::vector Base; + Base m_features; + +#ifdef CGAL_LINKED_WITH_TBB + tbb::mutex m_mutex; + void mutex_lock() { m_mutex.lock(); } + void mutex_unlock() { m_mutex.unlock(); } +#else // CGAL_LINKED_WITH_TBB + void mutex_lock() { } + void mutex_unlock() { } +#endif // CGAL_LINKED_WITH_TBB + +public: + + Feature_set() { } + virtual ~Feature_set() { } + + /*! + \brief Adds a feature. + + \tparam Feature type of the feature, inherited from + `Classification::Feature_base`. + + \tparam T types of the parameters of the feature's constructor + (with the exception of the first parameter that is always of type + `ItemRange&` and that is automatically passed by the classifier to + the feature's constructor). + + \param t parameters of the feature's constructor (with the + exception of the first parameter that is always the input item + range and that is automatically passed by the classifier to the + feature's constructor). + + \return a handle to the newly added feature. + */ +#if (!defined(CGAL_CFG_NO_CPP0X_VARIADIC_TEMPLATES) && !defined(CGAL_CFG_NO_CPP0X_RVALUE_REFERENCE)) || DOXYGEN_RUNNING + template + Feature_handle add (T&& ... t) + { + Feature_handle fh (new Feature(std::forward(t)...)); + mutex_lock(); + m_features.push_back (fh); + mutex_unlock(); + return fh; + } +#else + template + Feature_handle add (T1& t1) + { + Feature_handle fh (new Feature(t1)); + mutex_lock(); + m_features.push_back (fh); + mutex_unlock(); + return fh; + } + template + Feature_handle add (T1& t1, T2& t2) + { + Feature_handle fh (new Feature(t1, t2)); + mutex_lock(); + m_features.push_back (fh); + mutex_unlock(); + return fh; + } + template + Feature_handle add (T1& t1, T2& t2, T3& t3) + { + Feature_handle fh (new Feature(t1, t2, t3)); + mutex_lock(); + m_features.push_back (fh); + mutex_unlock(); + return fh; + } + template + Feature_handle add (T1& t1, T2& t2, T3& t3, T4& t4) + { + Feature_handle fh (new Feature(t1, t2, t3, t4)); + mutex_lock(); + m_features.push_back (fh); + mutex_unlock(); + return fh; + } + template + Feature_handle add (T1& t1, T2& t2, T3& t3, T4& t4, T5& t5) + { + Feature_handle fh (new Feature(t1, t2, t3, t4, t5)); + mutex_lock(); + m_features.push_back (fh); + mutex_unlock(); + return fh; + } +#endif + + /*! + \brief Removes a feature. + + \param feature the handle to feature type that must be removed. + + \return `true` if the feature was correctly removed, `false` if + its handle was not found. + */ + bool remove (Feature_handle feature) + { + for (std::size_t i = 0; i < m_features.size(); ++ i) + if (m_features[i] == feature) + { + m_features.erase (m_features.begin() + i); + return true; + } + return false; + } + + /*! + \brief Returns how many features are defined. + */ + std::size_t size() const + { + return m_features.size(); + } + + + /*! + \brief Returns the \f$i^{th}\f$ feature. + */ + Feature_handle operator[](std::size_t i) const + { + return m_features[i]; + } + + /*! + \brief Removes all features. + */ + void clear () + { + m_features.clear(); + } + +}; + + + +} // namespace Classification + +} // namespace CGAL + +#endif // CGAL_CLASSIFICATION_FEATURE_SET_H diff --git a/Classification/include/CGAL/Classification/Label.h b/Classification/include/CGAL/Classification/Label.h index e81e24b2c2a..c8bb503351d 100644 --- a/Classification/include/CGAL/Classification/Label.h +++ b/Classification/include/CGAL/Classification/Label.h @@ -20,8 +20,6 @@ #ifndef CGAL_CLASSIFICATION_LABEL_H #define CGAL_CLASSIFICATION_LABEL_H -#include - #include namespace CGAL { @@ -37,13 +35,9 @@ defined as a set of relationships with classification features. */ class Label { -public: - private: - /// \cond SKIP_IN_MANUAL + std::string m_name; - std::map m_feature_effects; - /// \endcond public: @@ -52,44 +46,7 @@ public: */ 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::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::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 diff --git a/Classification/include/CGAL/Classification/Label_set.h b/Classification/include/CGAL/Classification/Label_set.h new file mode 100644 index 00000000000..8b3485439cc --- /dev/null +++ b/Classification/include/CGAL/Classification/Label_set.h @@ -0,0 +1,113 @@ +// 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_SET_H +#define CGAL_CLASSIFICATION_LABEL_SET_H + +#include + +#include + +namespace CGAL { + +namespace Classification { + +class Label_set +{ + typedef std::vector Base; + Base m_labels; + +public: + + Label_set() { } + virtual ~Label_set() { } + + /*! + \brief Adds a label. + + \param name name of the label. + + \return a handle to the newly added label. + */ + Label_handle add (const char* name) + { + Label_handle out (new Classification::Label (name)); + m_labels.push_back (out); + return out; + } + + /*! + \brief Removes a label. + + \param label the handle to the label that must be removed. + + \return `true` if the label was correctly removed, + `false` if its handle was not found. + */ + bool remove (Label_handle label) + { + std::size_t idx = (std::size_t)(-1); + for (std::size_t i = 0; i < m_labels.size(); ++ i) + if (m_labels[i] == label) + { + m_labels.erase (m_labels.begin() + i); + idx = i; + break; + } + if (idx == (std::size_t)(-1)) + return false; + + return true; + } + + /*! + \brief Returns how many labels are defined. + */ + std::size_t size () const + { + return m_labels.size(); + } + + /*! + \brief Returns the \f$i^{th}\f$ label. + */ + Label_handle operator[] (std::size_t i) const + { + return m_labels[i]; + } + + + /*! + \brief Removes all labels. + */ + void clear () + { + m_labels.clear(); + } + + +}; + + + +} // namespace Classification + +} // namespace CGAL + +#endif // CGAL_CLASSIFICATION_LABEL_SET_H diff --git a/Classification/include/CGAL/Classification/Point_set_feature_generator.h b/Classification/include/CGAL/Classification/Point_set_feature_generator.h new file mode 100644 index 00000000000..3c016722747 --- /dev/null +++ b/Classification/include/CGAL/Classification/Point_set_feature_generator.h @@ -0,0 +1,752 @@ +// 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_FEATURE_GENERATOR_H +#define CGAL_CLASSIFICATION_POINT_SET_FEATURE_GENERATOR_H + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#include +#include +#include +#include +#include +#include + +#include +#include + +#ifdef CGAL_LINKED_WITH_TBB +#include +#include +#endif // CGAL_LINKED_WITH_TBB + +namespace CGAL { + +namespace Classification { + +/*! + \ingroup PkgClassification + + \brief Classifies a point set based on a set of features and a set + of labels. + + This class specializes `Classifier` to point sets. It takes care of + generating necessary data structures and automatically generate a + set of generic features. Features can be generated at multiple + scales to increase the reliability of the classification. + + \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 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). + \tparam DiagonalizeTraits model of `DiagonalizeTraits` used + for matrix diagonalization. + +*/ +template > +class Point_set_feature_generator +{ + +public: + typedef typename Geom_traits::Iso_cuboid_3 Iso_cuboid_3; + + /// \cond SKIP_IN_MANUAL + typedef typename PointRange::const_iterator Iterator; + typedef typename PointMap::value_type Point; + /// \endcond + + typedef Classification::Planimetric_grid + Planimetric_grid; + typedef Classification::Point_set_neighborhood + Neighborhood; + typedef Classification::Local_eigen_analysis + Local_eigen_analysis; + + /// \cond SKIP_IN_MANUAL + typedef Classification::Feature_handle Feature_handle; + typedef Classification::Label Label; + typedef Classification::Label_handle Label_handle; + + typedef Classification::Feature::Anisotropy + Anisotropy; + typedef Classification::Feature::Distance_to_plane + Distance_to_plane; + typedef Classification::Feature::Eigentropy + Eigentropy; + typedef Classification::Feature::Elevation + Elevation; + typedef Classification::Feature::Linearity + Linearity; + typedef Classification::Feature::Omnivariance + Omnivariance; + typedef Classification::Feature::Planarity + Planarity; + typedef Classification::Feature::Sphericity + Sphericity; + typedef Classification::Feature::Sum_eigenvalues + Sum_eigen; + typedef Classification::Feature::Surface_variation + Surface_variation; + typedef Classification::Feature::Vertical_dispersion + Dispersion; + typedef Classification::Feature::Verticality + Verticality; + typedef typename Classification::RGB_Color RGB_Color; + /// \endcond + +private: + + struct Scale + { + Neighborhood* neighborhood; + Planimetric_grid* grid; + Local_eigen_analysis* eigen; + double voxel_size; + + Scale (const PointRange& input, PointMap point_map, + const Iso_cuboid_3& bbox, double voxel_size) + : voxel_size (voxel_size) + { + CGAL::Real_timer t; + t.start(); + if (voxel_size < 0.) + neighborhood = new Neighborhood (input, point_map); + else + neighborhood = new Neighborhood (input, point_map, voxel_size); + t.stop(); + + if (voxel_size < 0.) + CGAL_CLASSIFICATION_CERR << "Neighborhood computed in " << t.time() << " second(s)" << std::endl; + else + CGAL_CLASSIFICATION_CERR << "Neighborhood with voxel size " << voxel_size + << " computed in " << t.time() << " second(s)" << std::endl; + t.reset(); + t.start(); + + eigen = new Local_eigen_analysis (input, point_map, neighborhood->k_neighbor_query(6)); + double range = eigen->mean_range(); + if (this->voxel_size < 0) + this->voxel_size = range; + t.stop(); + CGAL_CLASSIFICATION_CERR << "Eigen values computed in " << t.time() << " second(s)" << std::endl; + CGAL_CLASSIFICATION_CERR << "Range = " << range << std::endl; + t.reset(); + t.start(); + + grid = new Planimetric_grid (input, point_map, bbox, this->voxel_size); + t.stop(); + CGAL_CLASSIFICATION_CERR << "Planimetric grid computed in " << t.time() << " second(s)" << std::endl; + t.reset(); + } + ~Scale() + { + delete neighborhood; + delete grid; + delete eigen; + } + + double grid_resolution() const { return voxel_size; } + double radius_neighbors() const { return voxel_size * 5; } + double radius_dtm() const { return voxel_size * 100; } + + }; + + Iso_cuboid_3 m_bbox; + std::vector m_scales; + +#ifdef CGAL_LINKED_WITH_TBB + tbb::task_group* m_tasks; +#endif + + struct Feature_adder + { + mutable Point_set_feature_generator* generator; + std::size_t scale; + + Feature_adder (Point_set_feature_generator* generator, std::size_t scale) + : generator (generator), scale (scale) { } + + virtual ~Feature_adder() { } + + virtual void operator()() const = 0; + }; + friend Feature_adder; + + const PointRange& m_input; + PointMap m_point_map; + std::vector m_adders; + Feature_set* m_features; + +public: + + + /// \name Constructor + /// @{ + + /*! + \brief Initializes a classification object. + + \param input input range. + + \param point_map property map to access the input points. + */ + /*! + \brief Generate all possible features from an input range. + + The size of the smallest scale is automatically estimated and the + data structures needed (`Neighborhood`, `Planimetric_grid` and + `Local_eigen_analysis`) are computed at `nb_scales` recursively + larger scales. At each scale, the following features are + generated: + + - `CGAL::Classification::Feature::Anisotropy` + - `CGAL::Classification::Feature::Distance_to_plane` + - `CGAL::Classification::Feature::Eigentropy` + - `CGAL::Classification::Feature::Elevation` + - `CGAL::Classification::Feature::Linearity` + - `CGAL::Classification::Feature::Omnivariance` + - `CGAL::Classification::Feature::Planarity` + - `CGAL::Classification::Feature::Sphericity` + - `CGAL::Classification::Feature::Sum_eigenvalues` + - `CGAL::Classification::Feature::Surface_variation` + - `CGAL::Classification::Feature::Vertical_dispersion` based on eigenvalues + + If normal vectors are provided (if `VectorMap` is different from + `CGAL::Default`), the following feature is generated at each + scale: + + - `CGAL::Classification::Feature::Vertical_dispersion` based on normal vectors + + If colors are provided (if `ColorMap` is different from + `CGAL::Default`), the following features are generated at each + scale: + + - 9 features `CGAL::Classification::Feature::Hsv` on + channel 0 (hue) with mean ranging from 0° to 360° and standard + deviation of 22.5. + + - 5 features `CGAL::Classification::Feature::Hsv` on + channel 1 (saturation) with mean ranging from 0 to 100 and standard + deviation of 12.5. + + - 5 features `CGAL::Classification::Feature::Hsv` on channel 2 + (value) with mean ranging from 0 to 100 and standard deviation + of 12.5. + + If echo numbers are provided (if `EchoMap` is different from + `CGAL::Default`), the following feature is computed at each + scale: + + - `CGAL::Classification::Feature::Echo_scatter` + + \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`. + \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`. + \tparam EchoMap model of `ReadablePropertyMap` whose key type is + the value type of the iterator of `PointRange` and value type is + `std::size_t`. + \param nb_scales number of scales to compute. + \param normal_map property map to access the normal vectors of the input points (if any). + \param color_map property map to access the colors of the input points (if any). + \param echo_map property map to access the echo values of the input points (if any). + */ + template + Point_set_feature_generator(Feature_set& features, + std::size_t nb_scales, + const PointRange& input, + PointMap point_map, + VectorMap normal_map = VectorMap(), + ColorMap color_map = ColorMap(), + EchoMap echo_map = EchoMap()) + : m_input (input), m_point_map (point_map), m_features (&features) + { + m_bbox = CGAL::bounding_box + (boost::make_transform_iterator (m_input.begin(), CGAL::Property_map_to_unary_function(m_point_map)), + boost::make_transform_iterator (m_input.end(), CGAL::Property_map_to_unary_function(m_point_map))); + + typedef typename Default::Get::type + Vmap; + typedef typename Default::Get::type + Cmap; + typedef typename Default::Get::type + Emap; + + generate_features_impl (nb_scales, + get_parameter(normal_map), + get_parameter(color_map), + get_parameter(echo_map)); + } + + /// @} + + /// \cond SKIP_IN_MANUAL + virtual ~Point_set_feature_generator() + { + clear(); + } + /// \endcond + + /// \name Data Structures and Parameters + /// @{ + + /*! + \brief Returns the bounding box of the input point set. + */ + const Iso_cuboid_3& bbox() const { return m_bbox; } + /*! + \brief Returns the neighborhood structure at scale `scale`. + + \note `generate_features()` must have been called before calling + this method. + */ + const Neighborhood& neighborhood(std::size_t scale = 0) const { return (*m_scales[scale]->neighborhood); } + /*! + \brief Returns the planimetric grid structure at scale `scale`. + + \note `generate_features()` must have been called before calling + this method. + */ + const Planimetric_grid& grid(std::size_t scale = 0) const { return *(m_scales[scale]->grid); } + /*! + \brief Returns the local eigen analysis structure at scale `scale`. + + \note `generate_features()` must have been called before calling + this method. + */ + const Local_eigen_analysis& eigen(std::size_t scale = 0) const { return *(m_scales[scale]->eigen); } + /*! + \brief Returns the number of scales that were computed. + */ + std::size_t number_of_scales() const { return m_scales.size(); } + + /*! + \brief Returns the grid resolution at scale `scale`. This + resolution is the length and width of a cell of the + `Planimetric_grid` defined at this scale. + + \note `generate_features()` must have been called before calling + this method. + */ + double grid_resolution(std::size_t scale = 0) const { return m_scales[scale]->grid_resolution(); } + /*! + + \brief Returns the radius used for neighborhood queries at scale + `scale`. This radius is the smallest radius that is relevant from + a geometric point of view at this scale (that is to say that + encloses a few cells of `Planimetric_grid`). + + \note `generate_features()` must have been called before calling + this method. + */ + double radius_neighbors(std::size_t scale = 0) const { return m_scales[scale]->radius_neighbors(); } + /*! + \brief Returns the radius used for digital terrain modeling at + scale `scale`. This radius represents the minimum size of a + building at this scale. + + \note `generate_features()` must have been called before calling + this method. + */ + double radius_dtm(std::size_t scale = 0) const { return m_scales[scale]->radius_dtm(); } + + + +private: + + void clear() + { + for (std::size_t i = 0; i < m_scales.size(); ++ i) + delete m_scales[i]; + m_scales.clear(); + } + + void generate_point_based_features () + { + + generate_multiscale_feature_variant_0 (); + generate_multiscale_feature_variant_0 (); + generate_multiscale_feature_variant_0 (); + generate_multiscale_feature_variant_0 (); + generate_multiscale_feature_variant_0 (); + generate_multiscale_feature_variant_0 (); + generate_multiscale_feature_variant_0 (); + generate_multiscale_feature_variant_0 (); + + generate_multiscale_feature_variant_1 (); + generate_multiscale_feature_variant_2 (); + generate_multiscale_feature_variant_3 (); + } + + template + void launch_feature_computation (FeatureAdder* adder) + { + m_adders.push_back (adder); + +#ifndef CGAL_LINKED_WITH_TBB + CGAL_static_assertion_msg (!(boost::is_convertible::value), + "Parallel_tag is enabled but TBB is unavailable."); +#else + if (boost::is_convertible::value) + { + m_tasks->run (*adder); + } + else +#endif + { + (*adder)(); + } + } + + template + struct Feature_adder_verticality : public Feature_adder + { + using Feature_adder::generator; + using Feature_adder::scale; + VectorMap normal_map; + + // TODO! + Feature_adder_verticality (Point_set_feature_generator* generator, VectorMap normal_map, std::size_t scale) + : Feature_adder (generator, scale), normal_map (normal_map) { } + + void operator() () const + { + Feature_handle fh = generator->m_features->template add (generator->m_input, normal_map); + std::ostringstream oss; + oss << fh->name() << "_" << scale; + fh->set_name (oss.str()); + } + }; + + template + void generate_normal_based_features(VectorMap normal_map) + { + launch_feature_computation (new Feature_adder_verticality (this, normal_map, 0)); + } + + void generate_normal_based_features(const CGAL::Default_property_map&) + { + generate_multiscale_feature_variant_0 (); + } + + template + struct Feature_adder_color : public Feature_adder + { + typedef Classification::Feature::Hsv Hsv; + + using Feature_adder::generator; + using Feature_adder::scale; + ColorMap color_map; + std::size_t channel; + double mean; + double sd; + + // TODO! + Feature_adder_color (Point_set_feature_generator* generator, ColorMap color_map, std::size_t scale, + std::size_t channel, double mean, double sd) + : Feature_adder (generator, scale), color_map (color_map), + channel (channel), mean (mean), sd (sd) { } + + void operator() () const + { + Feature_handle fh = generator->m_features->template add (generator->m_input, color_map, channel, mean, sd); + std::ostringstream oss; + oss << fh->name() << "_" << scale; + fh->set_name (oss.str()); + } + }; + + template + void generate_color_based_features(ColorMap color_map) + { + for (std::size_t i = 0; i <= 8; ++ i) + launch_feature_computation (new Feature_adder_color (this, color_map, 0, + 0, 45 * i, 22.5)); + + for (std::size_t i = 0; i <= 4; ++ i) + launch_feature_computation (new Feature_adder_color (this, color_map, 0, + 1, 25 * i, 12.5)); + + for (std::size_t i = 0; i <= 4; ++ i) + launch_feature_computation (new Feature_adder_color (this, color_map, 0, + 2, 25 * i, 12.5)); + } + + void generate_color_based_features(const CGAL::Default_property_map&) + { + } + + + template + struct Feature_adder_echo : public Feature_adder + { + typedef Classification::Feature::Echo_scatter Echo_scatter; + + using Feature_adder::generator; + using Feature_adder::scale; + EchoMap echo_map; + + // TODO! + Feature_adder_echo (Point_set_feature_generator* generator, EchoMap echo_map, std::size_t scale) + : Feature_adder (generator, scale), echo_map (echo_map) { } + + void operator() () const + { + Feature_handle fh = generator->m_features->template add (generator->m_input, + echo_map, + generator->grid(scale), + generator->grid_resolution(scale), + generator->radius_neighbors(scale)); + std::ostringstream oss; + oss << fh->name() << "_" << scale; + fh->set_name (oss.str()); + } + }; + + + template + void generate_echo_based_features(EchoMap echo_map) + { + for (std::size_t i = 0; i < m_scales.size(); ++ i) + launch_feature_computation (new Feature_adder_echo (this, echo_map, i)); + } + + void generate_echo_based_features(const CGAL::Default_property_map&) + { + } + + + template + const T& get_parameter (const T& t) + { + return t; + } + + template + Default_property_map + get_parameter (const Default&) + { + return Default_property_map(); + } + + template + void generate_features_impl (std::size_t nb_scales, + VectorMap normal_map, + ColorMap color_map, + EchoMap echo_map) + { + CGAL::Real_timer t; t.start(); + + m_scales.reserve (nb_scales); + double voxel_size = - 1.; + + m_scales.push_back (new Scale (m_input, m_point_map, m_bbox, voxel_size)); + voxel_size = m_scales[0]->grid_resolution(); + + for (std::size_t i = 1; i < nb_scales; ++ i) + { + voxel_size *= 2; + m_scales.push_back (new Scale (m_input, m_point_map, m_bbox, voxel_size)); + } + t.stop(); + CGAL_CLASSIFICATION_CERR << "Scales computed in " << t.time() << " second(s)" << std::endl; + t.reset(); + + t.start(); + +#ifdef CGAL_LINKED_WITH_TBB + m_tasks = new tbb::task_group; +#endif + + generate_point_based_features (); + generate_normal_based_features (normal_map); + generate_color_based_features (color_map); + generate_echo_based_features (echo_map); + +#ifdef CGAL_LINKED_WITH_TBB + m_tasks->wait(); + delete m_tasks; +#endif + + t.stop(); + CGAL_CLASSIFICATION_CERR << "Features computed in " << t.time() << " second(s)" << std::endl; + for (std::size_t i = 0; i < m_adders.size(); ++ i) + delete m_adders[i]; + } + + template + struct Feature_adder_variant_0 : public Feature_adder + { + using Feature_adder::generator; + using Feature_adder::scale; + + Feature_adder_variant_0 (Point_set_feature_generator* generator, std::size_t scale) + : Feature_adder (generator, scale) { } + + void operator() () const + { + Feature_handle fh = generator->m_features->template add (generator->m_input, generator->eigen(scale)); + std::ostringstream oss; + oss << fh->name() << "_" << scale; + fh->set_name (oss.str()); + } + }; + + template + void generate_multiscale_feature_variant_0 () + { + for (std::size_t i = 0; i < m_scales.size(); ++ i) + launch_feature_computation (new Feature_adder_variant_0 (this, i)); + } + + template + struct Feature_adder_variant_1 : public Feature_adder + { + using Feature_adder::generator; + using Feature_adder::scale; + PointMap point_map; + + // TODO! + Feature_adder_variant_1 (Point_set_feature_generator* generator, PointMap point_map, std::size_t scale) + : Feature_adder (generator, scale), point_map (point_map) { } + + void operator() () const + { + Feature_handle fh = generator->m_features->template add (generator->m_input, point_map, + generator->eigen(scale)); + std::ostringstream oss; + oss << fh->name() << "_" << scale; + fh->set_name (oss.str()); + } + }; + + template + void generate_multiscale_feature_variant_1 () + { + for (std::size_t i = 0; i < m_scales.size(); ++ i) + launch_feature_computation (new Feature_adder_variant_1 (this, m_point_map, i)); + } + + template + struct Feature_adder_variant_2 : public Feature_adder + { + using Feature_adder::generator; + using Feature_adder::scale; + PointMap point_map; + + // TODO! + Feature_adder_variant_2 (Point_set_feature_generator* generator, PointMap point_map, std::size_t scale) + : Feature_adder (generator, scale), point_map (point_map) { } + + void operator() () const + { + Feature_handle fh = generator->m_features->template add + (generator->m_input, point_map, + generator->grid(scale), + generator->grid_resolution(scale), + generator->radius_neighbors(scale)); + std::ostringstream oss; + oss << fh->name() << "_" << scale; + fh->set_name (oss.str()); + } + }; + + template + void generate_multiscale_feature_variant_2 () + { + for (std::size_t i = 0; i < m_scales.size(); ++ i) + launch_feature_computation (new Feature_adder_variant_2 (this, m_point_map, i)); + } + + template + struct Feature_adder_variant_3 : public Feature_adder + { + using Feature_adder::generator; + using Feature_adder::scale; + PointMap point_map; + + // TODO! + Feature_adder_variant_3 (Point_set_feature_generator* generator, PointMap point_map, std::size_t scale) + : Feature_adder (generator, scale), point_map (point_map) { } + + void operator() () const + { + Feature_handle fh = generator->m_features->template add (generator->m_input, + point_map, + generator->grid(scale), + generator->grid_resolution(scale), + generator->radius_dtm(scale)); + std::ostringstream oss; + oss << fh->name() << "_" << scale; + fh->set_name (oss.str()); + } + }; + + template + void generate_multiscale_feature_variant_3 () + { + for (std::size_t i = 0; i < m_scales.size(); ++ i) + launch_feature_computation (new Feature_adder_variant_3 (this, m_point_map, i)); + } + + +}; + + +} // namespace Classification + +} // namespace CGAL + + + +#endif // CGAL_CLASSIFICATION_POINT_SET_FEATURE_GENERATOR_H diff --git a/Classification/include/CGAL/Classification/Sum_of_weighted_features_predicate.h b/Classification/include/CGAL/Classification/Sum_of_weighted_features_predicate.h new file mode 100644 index 00000000000..2d2b700d703 --- /dev/null +++ b/Classification/include/CGAL/Classification/Sum_of_weighted_features_predicate.h @@ -0,0 +1,731 @@ +// 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 CLASSIFICATION_SUM_OF_WEIGHTED_FEATURES_PREDICATE_H +#define CLASSIFICATION_SUM_OF_WEIGHTED_FEATURES_PREDICATE_H + +#include +#include + +#include +#include +#include +#include +#include +#include + +#ifdef CGAL_LINKED_WITH_TBB +#include +#include +#include +#include +#endif // CGAL_LINKED_WITH_TBB + +//#define CGAL_CLASSIFICATION_VERBOSE +#if defined(CGAL_CLASSIFICATION_VERBOSE) +#define CGAL_CLASSIFICATION_CERR std::cerr +#else +#define CGAL_CLASSIFICATION_CERR std::ostream(0) +#endif + +//#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 { + +class Sum_of_weighted_features_predicate +{ +public: + + 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 + }; + +private: + +#ifdef CGAL_LINKED_WITH_TBB + class Compute_worst_score_and_confidence + { + std::vector& m_training_set; + const Sum_of_weighted_features_predicate& m_predicate; + 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& training_set, + const Sum_of_weighted_features_predicate& predicate, + std::size_t label, + double& confidence, + std::size_t& nb_okay, + tbb::mutex& mutex) + : m_training_set (training_set) + , m_predicate (predicate) + , m_label (label) + , m_confidence (confidence) + , m_nb_okay (nb_okay) + , m_mutex (mutex) + { } + + void operator()(const tbb::blocked_range& r) const + { + for (std::size_t k = r.begin(); k != r.end(); ++ k) + { + std::vector > values; + + std::vector v; + m_predicate.probabilities (m_training_set[k], v); + + for(std::size_t l = 0; l < v.size(); ++ l) + values.push_back (std::make_pair (v[l], 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 + + + Label_set& m_labels; + Feature_set& m_features; + std::vector m_weights; + std::vector > m_effect_table; + mutable std::map m_map_labels; + mutable std::map m_map_features; + +public: + Sum_of_weighted_features_predicate (Label_set& labels, + Feature_set& features) + : m_labels (labels), m_features (features), + m_weights (features.size(), 1.), + m_effect_table (labels.size(), std::vector + (features.size(), + NEUTRAL)) + { + for (std::size_t i = 0; i < labels.size(); ++ i) + m_map_labels[labels[i]] = i; + for (std::size_t i = 0; i < features.size(); ++ i) + m_map_features[features[i]] = i; + } + + /*! + \brief Sets the weight of the feature (`weight` must be positive). + */ + void set_weight (Feature_handle feature, double weight) + { + m_weights[m_map_features[feature]] = weight; + } + /// \cond SKIP_IN_MANUAL + void set_weight (std::size_t feature, double weight) + { + m_weights[feature] = weight; + } + /// \endcond + + /*! + \brief Returns the weight of the feature. + */ + double weight (Feature_handle feature) const + { + return m_weights[m_map_features[feature]]; + } + /// \cond SKIP_IN_MANUAL + double weight (std::size_t feature) const + { + return m_weights[feature]; + } + /// \endcond + + /*! + \brief Sets the `effect` of `feature` on `label`. + */ + void set_effect (Label_handle label, Feature_handle feature, + Effect effect) + { + m_effect_table[m_map_labels[label]][m_map_features[feature]] = effect; + } + /// \cond SKIP_IN_MANUAL + void set_effect (std::size_t label, std::size_t feature, + Effect effect) + { + m_effect_table[label][feature] = effect; + } + /// \endcond + + /*! + \brief Returns the `effect` of `feature` on `label`. + */ + Effect effect (Label_handle label, Feature_handle feature) const + { + return m_effect_table[m_map_labels[label]][m_map_features[feature]]; + } + /// \cond SKIP_IN_MANUAL + Effect effect (std::size_t label, std::size_t feature) const + { + return m_effect_table[label][feature]; + } + /// \endcond + + void probabilities (std::size_t item_index, + std::vector& out) const + { + out.resize (m_labels.size()); + for (std::size_t l = 0; l < m_labels.size(); ++ l) + { + out[l] = 0.; + for (std::size_t f = 0; f < m_features.size(); ++ f) + if (weight(f) != 0.) + out[l] += value (l, f, item_index); + } + } + + + /*! + \brief Saves the current configuration in the stream `output`. + + This allows to easily save and recover a specific classification + configuration, that is to say: + + - The size of the smallest scale + - The features and their respective weights + - The labels and the effects of the features on them + + The output file is written in an XML format that is readable by + the `load_configuration()` method. + */ + void save_configuration (std::ostream& output) + { + boost::property_tree::ptree tree; + + for (std::size_t i = 0; i < m_features.size(); ++ i) + { + if (weight(m_features[i]) == 0) + continue; + boost::property_tree::ptree ptr; + + ptr.put("name", m_features[i]->name()); + ptr.put("weight", weight(m_features[i])); + tree.add_child("classification.features.feature", ptr); + } + + + for (std::size_t i = 0; i < m_labels.size(); ++ i) + { + boost::property_tree::ptree ptr; + ptr.put("name", m_labels[i]->name()); + for (std::size_t j = 0; j < m_features.size(); ++ j) + { + if (weight(j) == 0) + continue; + boost::property_tree::ptree ptr2; + ptr2.put("name", m_features[j]->name()); + Effect e = effect(i, j); + if (e == PENALIZING) + ptr2.put("effect", "penalized"); + else if (e == NEUTRAL) + ptr2.put("effect", "neutral"); + else if (e == FAVORING) + ptr2.put("effect", "favored"); + ptr.add_child("feature", ptr2); + } + tree.add_child("classification.labels.label", ptr); + } + + // Write property tree to XML file + boost::property_tree::xml_writer_settings settings(' ', 3); + boost::property_tree::write_xml(output, tree, settings); + } + + /*! + \brief Loads a configuration from the stream `input`. + + All data structures, features and labels specified in the input + stream `input` are instantiated if possible (in particular, + property maps needed should be provided), similarly to what is + done in `generate_features()`. + + The input file should be in the XML format written by the + `save_configuration()` method. + + \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`. + \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`. + \tparam EchoMap model of `ReadablePropertyMap` whose key type is + the value type of the iterator of `PointRange` and value type is + `std::size_t`. + \param input input stream. + \param normal_map property map to access the normal vectors of the input points (if any). + \param color_map property map to access the colors of the input points (if any). + \param echo_map property map to access the echo values of the input points (if any). + */ + bool load_configuration (std::istream& input, bool verbose = false) + { + bool out = true; + std::map map_n2l; + std::map map_n2f; + for (std::size_t i = 0; i < m_labels.size(); ++ i) + map_n2l.insert (std::make_pair (m_labels[i]->name(), i)); + for (std::size_t i = 0; i < m_features.size(); ++ i) + map_n2f.insert (std::make_pair (m_features[i]->name(), i)); + + boost::property_tree::ptree tree; + boost::property_tree::read_xml(input, tree); + + BOOST_FOREACH(boost::property_tree::ptree::value_type &v, tree.get_child("classification.features")) + { + std::string name = v.second.get("name"); + typename std::map::iterator + found = map_n2f.find (name); + if (found != map_n2f.end()) + m_weights[found->second] = v.second.get("weight"); + else + { + if (verbose) + std::cerr << "Warning: feature \"" << name << "\" in configuration file not found" << std::endl; + out = false; + } + } + + BOOST_FOREACH(boost::property_tree::ptree::value_type &v, tree.get_child("classification.labels")) + { + std::string label_name = v.second.get("name"); + typename std::map::iterator + found = map_n2l.find (label_name); + std::size_t l = 0; + if (found != map_n2l.end()) + l = found->second; + else + { + if (verbose) + std::cerr << "Warning: label \"" << label_name << "\" in configuration file not found" << std::endl; + out = false; + continue; + } + + BOOST_FOREACH(boost::property_tree::ptree::value_type &v2, v.second) + { + if (v2.first == "name") + continue; + + std::string feature_name = v2.second.get("name"); + + typename std::map::iterator + found2 = map_n2f.find (feature_name); + std::size_t f = 0; + if (found2 != map_n2f.end()) + f = found2->second; + else if (verbose) + { + if (verbose) + std::cerr << "Warning: feature \"" << feature_name << "\" in configuration file not found" << std::endl; + out = false; + continue; + } + std::string e = v2.second.get("effect"); + if (e == "penalized") + set_effect (l, f, PENALIZING); + else if (e == "neutral") + set_effect (l, f, NEUTRAL); + else + set_effect (l, f, FAVORING); + } + } + return out; + } + /// \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. + */ + template + double train (const std::vector& ground_truth, + std::size_t nb_tests = 300) + { + std::vector > training_sets (m_labels.size()); + for (std::size_t i = 0; i < ground_truth.size(); ++ i) + if (ground_truth[i] != std::size_t(-1)) + training_sets[ground_truth[i]].push_back (i); + + for (std::size_t i = 0; i < m_labels.size(); ++ i) + if (training_sets.size() <= i || training_sets[i].empty()) + std::cerr << "WARNING: \"" << m_labels[i]->name() << "\" doesn't have a training set." << std::endl; + + std::vector best_weights (m_features.size(), 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_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_features.size(); ++ j) + { + Feature_handle feature = m_features[j]; + best_weights[j] = weight(j); + + std::size_t nb_useful = 0; + double min = (std::numeric_limits::max)(); + double max = -(std::numeric_limits::max)(); + + set_weight(j, wmin); + for (std::size_t i = 0; i < 100; ++ i) + { + estimate_feature_effect(j, training_sets); + if (feature_useful(j)) + { + CGAL_CLASSTRAINING_CERR << "#"; + nb_useful ++; + min = (std::min) (min, weight(j)); + max = (std::max) (max, weight(j)); + } + else + CGAL_CLASSTRAINING_CERR << "-"; + set_weight(j, factor * weight(j)); + } + 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) + { + set_weight(j, 0.); + best_weights[j] = weight(j); + 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.) + { + set_weight(j, 0.5 * (feature_train.back().wmin + feature_train.back().wmax)); + best_weights[j] = weight(j); + } + else + set_weight(j, best_weights[j]); + estimate_feature_effect(j, training_sets); + } + + 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., training_sets); + + CGAL_CLASSIFICATION_CERR << "TRAINING GLOBALLY: Best score evolution: " << std::endl; + + CGAL_CLASSIFICATION_CERR << 100. * best_score << "% (found at initialization)" << std::endl; + + 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_features[current_feature_changed]; + + std::size_t nb_used = 0; + for (std::size_t j = 0; j < m_features.size(); ++ j) + { + if (j == current_feature_changed) + continue; + + set_weight(j, best_weights[j]); + estimate_feature_effect(j, training_sets); + if (feature_useful(j)) + nb_used ++; + else + set_weight(j, 0.); + } + + set_weight(current_feature_changed, tr.wmin); + for (std::size_t j = 0; j < nb_trials_per_feature; ++ j) + { + estimate_feature_effect(current_feature_changed, training_sets); + + double worst_confidence = 0., worst_score = 0.; + boost::tie (worst_confidence, worst_score) + = compute_worst_confidence_and_score (best_confidence, best_score, training_sets); + + 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_changed) ? 1 : 0) + << "/" << m_features.size() << " feature(s) used)" << std::endl; + for (std::size_t k = 0; k < m_features.size(); ++ k) + best_weights[k] = weight(k); + } + set_weight(current_feature_changed, weight(current_feature_changed) * tr.factor); + } + } + + for (std::size_t i = 0; i < best_weights.size(); ++ i) + set_weight(i, best_weights[i]); + + estimate_features_effects(training_sets); + + 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_features[i]; + CGAL_CLASSTRAINING_CERR << "FEATURE " << feature->name() << ": " << best_weights[i] << std::endl; + set_weight(i, best_weights[i]); + + Effect side = effect(0, i); + bool to_remove = true; + for (std::size_t j = 0; j < m_labels.size(); ++ j) + { + Label_handle clabel = m_labels[j]; + if (effect(j,i) == FAVORING) + CGAL_CLASSTRAINING_CERR << " * Favored for "; + else if (effect(j,i) == PENALIZING) + CGAL_CLASSTRAINING_CERR << " * Penalized for "; + else + CGAL_CLASSTRAINING_CERR << " * Neutral for "; + if (effect(j,i) != 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_features.size() << " are useless" << std::endl; + + return best_score; + } + + /// @} + +private: + + double value (std::size_t label, std::size_t feature, std::size_t index) const + { + if (m_effect_table[label][feature] == FAVORING) + return favored (feature, index); + else if (m_effect_table[label][feature] == PENALIZING) + return penalized (feature, index); + else + return ignored (feature, index); + } + + double normalized (std::size_t feature, std::size_t index) const + { + return (std::max) (0., (std::min) (1., m_features[feature]->value(index) / m_weights[feature])); + } + double favored (std::size_t feature, std::size_t index) const + { + return (1. - normalized (feature, index)); + } + double penalized (std::size_t feature, std::size_t index) const + { + return normalized (feature, index); + } + double ignored (std::size_t, std::size_t) const + { + return 0.5; + } + + void estimate_features_effects(std::vector >& training_sets) + { + for (std::size_t i = 0; i < m_features.size(); ++ i) + estimate_feature_effect (i, training_sets); + } + + void estimate_feature_effect (std::size_t feature, + std::vector >& training_sets) + { + std::vector mean (m_labels.size(), 0.); + + for (std::size_t j = 0; j < m_labels.size(); ++ j) + { + for (std::size_t k = 0; k < training_sets[j].size(); ++ k) + { + double val = normalized(feature, training_sets[j][k]); + mean[j] += val; + } + mean[j] /= training_sets[j].size(); + } + + std::vector sd (m_labels.size(), 0.); + + for (std::size_t j = 0; j < m_labels.size(); ++ j) + { + Label_handle clabel = m_labels[j]; + + for (std::size_t k = 0; k < training_sets[j].size(); ++ k) + { + double val = normalized(feature, training_sets[j][k]); + sd[j] += (val - mean[j]) * (val - mean[j]); + } + sd[j] = std::sqrt (sd[j] / training_sets[j].size()); + if (mean[j] - sd[j] > 0.75) + set_effect (j, feature, FAVORING); + else if (mean[j] + sd[j] < 0.25) + set_effect (j, feature, PENALIZING); + else + set_effect (j, feature, NEUTRAL); + } + } + + template + std::pair compute_worst_confidence_and_score (double lower_conf, double lower_score, + std::vector >& training_sets) + { + double worst_confidence = (std::numeric_limits::max)(); + double worst_score = (std::numeric_limits::max)(); + + for (std::size_t j = 0; j < m_labels.size(); ++ j) + { + double confidence = 0.; + std::size_t nb_okay = 0; + +#ifndef CGAL_LINKED_WITH_TBB + CGAL_static_assertion_msg (!(boost::is_convertible::value), + "Parallel_tag is enabled but TBB is unavailable."); +#else + if (boost::is_convertible::value) + { + tbb::mutex mutex; + Compute_worst_score_and_confidence f(training_sets[j], *this, j, confidence, nb_okay, mutex); + tbb::parallel_for(tbb::blocked_range(0, training_sets[j].size ()), f); + } + else +#endif + { + for (std::size_t k = 0; k < training_sets[j].size(); ++ k) + { + std::vector > values; + + std::vector v; + probabilities (training_sets[j][k], v); + + for(std::size_t l = 0; l < m_labels.size(); ++ l) + values.push_back (std::make_pair (v[l], 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)(training_sets[j].size()); + confidence /= (double)(training_sets[j].size() * m_features.size()); + + 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 (std::size_t feature) + { + Effect side = effect(0, feature); + for (std::size_t k = 1; k < m_labels.size(); ++ k) + if (effect(k, feature) != side) + return true; + return false; + } + +}; + +} + +} + +#endif // CLASSIFICATION_SUM_OF_WEIGHTED_FEATURES_PREDICATE_H diff --git a/Classification/include/CGAL/Classification/classify.h b/Classification/include/CGAL/Classification/classify.h new file mode 100644 index 00000000000..05f8f796762 --- /dev/null +++ b/Classification/include/CGAL/Classification/classify.h @@ -0,0 +1,438 @@ +// 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_CLASSIFY_H +#define CGAL_CLASSIFICATION_CLASSIFY_H + +#include + +#include + +#ifdef CGAL_LINKED_WITH_TBB +#include +#include +#include +#include +#endif // CGAL_LINKED_WITH_TBB + +namespace CGAL { + +namespace Classification { + +namespace internal { + + template + class Classifier + { + const Label_set& m_labels; + const ClassificationPredicate& m_predicate; + std::vector& m_out; + + public: + + Classifier (const Label_set& labels, + const ClassificationPredicate& predicate, + std::vector& out) + : m_labels (labels), m_predicate (predicate), m_out (out) + { } + +#ifdef CGAL_LINKED_WITH_TBB + void operator()(const tbb::blocked_range& r) const + { + for (std::size_t s = r.begin(); s != r.end(); ++ s) + apply(s); + } +#endif // CGAL_LINKED_WITH_TBB + + inline void apply (std::size_t s) const + { + std::size_t nb_class_best=0; + std::vector values; + m_predicate.probabilities (s, values); + + double val_class_best = (std::numeric_limits::max)(); + for(std::size_t k = 0; k < m_labels.size(); ++ k) + { + if(val_class_best > values[k]) + { + val_class_best = values[k]; + nb_class_best = k; + } + } + + m_out[s] = nb_class_best; + } + + }; + + template + class Classifier_local_smoothing_preprocessing + { + const Label_set& m_labels; + const ClassificationPredicate& m_predicate; + std::vector >& m_values; + + public: + + Classifier_local_smoothing_preprocessing + (const Label_set& labels, + const ClassificationPredicate& predicate, + std::vector >& values) + : m_labels (labels), m_predicate (predicate), m_values (values) + { } + +#ifdef CGAL_LINKED_WITH_TBB + void operator()(const tbb::blocked_range& r) const + { + for (std::size_t s = r.begin(); s != r.end(); ++ s) + apply (s); + } +#endif + + inline void apply (std::size_t s) const + { + std::vector values; + m_predicate.probabilities(s, values); + for(std::size_t k = 0; k < m_labels.size(); ++ k) + m_values[k][s] = values[k]; + } + }; + + template + class Classifier_local_smoothing + { + const ItemRange& m_input; + const ItemMap m_item_map; + const Label_set& m_labels; + const std::vector >& m_values; + const NeighborQuery& m_neighbor_query; + std::vector& m_out; + + public: + + Classifier_local_smoothing (const ItemRange& input, + ItemMap item_map, + const Label_set& labels, + const std::vector >& values, + const NeighborQuery& neighbor_query, + std::vector& out) + : m_input (input), m_item_map (item_map), m_labels (labels), + m_values(values), + m_neighbor_query (neighbor_query), + m_out (out) + { } + +#ifdef CGAL_LINKED_WITH_TBB + void operator()(const tbb::blocked_range& r) const + { + for (std::size_t s = r.begin(); s != r.end(); ++ s) + apply (s); + } +#endif + + inline void apply (std::size_t s) const + { + std::vector neighbors; + m_neighbor_query (get (m_item_map, *(m_input.begin()+s)), std::back_inserter (neighbors)); + + std::vector mean (m_values.size(), 0.); + for (std::size_t n = 0; n < neighbors.size(); ++ n) + for (std::size_t j = 0; j < m_values.size(); ++ j) + mean[j] += m_values[j][neighbors[n]]; + + std::size_t nb_class_best=0; + double val_class_best = (std::numeric_limits::max)(); + for(std::size_t k = 0; k < mean.size(); ++ k) + { + mean[k] /= neighbors.size(); + if(val_class_best > mean[k]) + { + val_class_best = mean[k]; + nb_class_best = k; + } + } + + m_out[s] = nb_class_best; + } + + + }; + + template + class Classifier_graphcut + { + const ItemRange& m_input; + ItemMap m_item_map; + const Label_set& m_labels; + const ClassificationPredicate& m_predicate; + const NeighborQuery& m_neighbor_query; + double m_weight; + const std::vector >& m_indices; + const std::vector >& m_input_to_indices; + std::vector& m_out; + +#ifdef CGAL_DO_NOT_USE_BOYKOV_KOLMOGOROV_MAXFLOW_SOFTWARE + typedef CGAL::internal::Alpha_expansion_graph_cut_boost Alpha_expansion; +#else + typedef CGAL::internal::Alpha_expansion_graph_cut_boykov_kolmogorov Alpha_expansion; +#endif + + public: + + Classifier_graphcut (const ItemRange& input, + ItemMap item_map, + const Label_set& labels, + const ClassificationPredicate& predicate, + const NeighborQuery& neighbor_query, + double weight, + const std::vector >& indices, + const std::vector >& input_to_indices, + std::vector& out) + : m_input (input), m_item_map (item_map), m_labels (labels), + m_predicate (predicate), m_neighbor_query (neighbor_query), + m_weight (weight), m_indices (indices), m_input_to_indices (input_to_indices), m_out (out) + { } + +#ifdef CGAL_LINKED_WITH_TBB + void operator()(const tbb::blocked_range& r) const + { + for (std::size_t s = r.begin(); s != r.end(); ++ s) + apply(s); + } +#endif // CGAL_LINKED_WITH_TBB + + + inline void apply (std::size_t sub) const + { + if (m_indices[sub].empty()) + return; + + std::vector > edges; + std::vector edge_weights; + std::vector > probability_matrix + (m_labels.size(), std::vector(m_indices[sub].size(), 0.)); + std::vector assigned_label (m_indices[sub].size()); + + for (std::size_t j = 0; j < m_indices[sub].size(); ++ j) + { + std::size_t s = m_indices[sub][j]; + + std::vector neighbors; + + m_neighbor_query (get(m_item_map, *(m_input.begin()+s)), std::back_inserter (neighbors)); + + for (std::size_t i = 0; i < neighbors.size(); ++ i) + if (sub == m_input_to_indices[neighbors[i]].first + && j != m_input_to_indices[neighbors[i]].second) + { + edges.push_back (std::make_pair (j, m_input_to_indices[neighbors[i]].second)); + edge_weights.push_back (m_weight); + } + + std::vector values; + m_predicate.probabilities(s, values); + std::size_t nb_class_best = 0; + double val_class_best = (std::numeric_limits::max)(); + for(std::size_t k = 0; k < m_labels.size(); ++ k) + { + double value = values[k]; + probability_matrix[k][j] = value; + + if(val_class_best > value) + { + val_class_best = value; + nb_class_best = k; + } + } + assigned_label[j] = nb_class_best; + } + + Alpha_expansion graphcut; + graphcut(edges, edge_weights, probability_matrix, assigned_label); + + for (std::size_t i = 0; i < assigned_label.size(); ++ i) + m_out[m_indices[sub][i]] = assigned_label[i]; + } + + }; + +} // namespace internal + + + + template + void classify (const ItemRange& input, + const Label_set& labels, + const ClassificationPredicate& predicate, + std::vector& output) + { + output.resize(input.size()); + + internal::Classifier + f (labels, predicate, output); + +#ifndef CGAL_LINKED_WITH_TBB + CGAL_static_assertion_msg (!(boost::is_convertible::value), + "Parallel_tag is enabled but TBB is unavailable."); +#else + if (boost::is_convertible::value) + { + tbb::parallel_for(tbb::blocked_range(0, input.size ()), f); + } + else +#endif + { + for (std::size_t i = 0; i < input.size(); ++ i) + f.apply(i); + } + } + + template + void classify_with_local_smoothing (const ItemRange& input, + const ItemMap item_map, + const Label_set& labels, + const ClassificationPredicate& predicate, + const NeighborQuery& neighbor_query, + std::vector& output) + { + output.resize(input.size()); + + std::vector > values + (labels.size(), std::vector (input.size(), -1.)); + internal::Classifier_local_smoothing_preprocessing + f1 (labels, predicate, values); + internal::Classifier_local_smoothing + f2 (input, item_map, labels, values, neighbor_query, output); + +#ifndef CGAL_LINKED_WITH_TBB + CGAL_static_assertion_msg (!(boost::is_convertible::value), + "Parallel_tag is enabled but TBB is unavailable."); +#else + if (boost::is_convertible::value) + { + tbb::parallel_for(tbb::blocked_range(0, input.size ()), f1); + tbb::parallel_for(tbb::blocked_range(0, input.size ()), f2); + } + else +#endif + { + for (std::size_t i = 0; i < input.size(); ++ i) + f1.apply(i); + for (std::size_t i = 0; i < input.size(); ++ i) + f2.apply(i); + } + } + + template + void classify_with_graphcut (const ItemRange& input, + const ItemMap item_map, + const ItemWithBboxMap bbox_map, + const Label_set& labels, + const ClassificationPredicate& predicate, + const NeighborQuery& neighbor_query, + const double weight, + const std::size_t min_number_of_subdivisions, + std::vector& output) + { + CGAL::Bbox_3 bbox = CGAL::bbox_3 + (boost::make_transform_iterator (input.begin(), CGAL::Property_map_to_unary_function(bbox_map)), + boost::make_transform_iterator (input.end(), CGAL::Property_map_to_unary_function(bbox_map))); + + double Dx = bbox.xmax() - bbox.xmin(); + double Dy = bbox.ymax() - bbox.ymin(); + double A = Dx * Dy; + double a = A / min_number_of_subdivisions; + double l = std::sqrt(a); + std::size_t nb_x = std::size_t(Dx / l) + 1; + std::size_t nb_y = std::size_t((A / nb_x) / a) + 1; + std::size_t nb = nb_x * nb_y; + + std::vector bboxes; + bboxes.reserve(nb); + for (std::size_t x = 0; x < nb_x; ++ x) + for (std::size_t y = 0; y < nb_y; ++ y) + { + bboxes.push_back + (CGAL::Bbox_3 (bbox.xmin() + Dx * (x / double(nb_x)), + bbox.ymin() + Dy * (y / double(nb_y)), + bbox.zmin(), + bbox.xmin() + Dx * ((x+1) / double(nb_x)), + bbox.ymin() + Dy * ((y+1) / double(nb_y)), + bbox.zmax())); + } + + std::cerr << "Number of divisions = " << nb_x * nb_y << std::endl; + std::cerr << " -> Size of division: " << Dx / nb_x << " " << Dy / nb_y << std::endl; + + std::vector > indices (nb); + std::vector > input_to_indices(input.size()); + + for (std::size_t s = 0; s < input.size(); ++ s) + { + CGAL::Bbox_3 b = get(bbox_map, *(input.begin() + s)).bbox(); + + for (std::size_t i = 0; i < bboxes.size(); ++ i) + if (CGAL::do_overlap (b, bboxes[i])) + { + input_to_indices[s] = std::make_pair (i, indices[i].size()); + indices[i].push_back (s); + break; + } + } + + output.resize (input.size()); + + internal::Classifier_graphcut + f (input, item_map, labels, predicate, neighbor_query, weight, indices, input_to_indices, output); + +#ifndef CGAL_LINKED_WITH_TBB + CGAL_static_assertion_msg (!(boost::is_convertible::value), + "Parallel_tag is enabled but TBB is unavailable."); +#else + if (boost::is_convertible::value) + { + tbb::parallel_for(tbb::blocked_range(0, indices.size ()), f); + } + else +#endif + { + for (std::size_t sub = 0; sub < indices.size(); ++ sub) + f.apply (sub); + } + } + + +} + +} + +#endif // CGAL_CLASSIFICATION_CLASSIFY_H diff --git a/Classification/include/CGAL/Classifier.h b/Classification/include/CGAL/Classifier.h deleted file mode 100644 index d7fa6af9cc7..00000000000 --- a/Classification/include/CGAL/Classifier.h +++ /dev/null @@ -1,1087 +0,0 @@ -// 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_CLASSIFIER_H -#define CGAL_CLASSIFIER_H - -#include -#include -#include -#include -#include -#include -#include - -#include -#include -#include -#include - -#include -#include -#include - -#include -#include - -#include - -#ifdef CGAL_LINKED_WITH_TBB -#include -#include -#include -#include -#endif // CGAL_LINKED_WITH_TBB - -//#define CGAL_CLASSIFICATION_VERBOSE -#if defined(CGAL_CLASSIFICATION_VERBOSE) -#define CGAL_CLASSIFICATION_CERR std::cerr -#else -#define CGAL_CLASSIFICATION_CERR std::ostream(0) -#endif - -namespace CGAL { - -/*! -\ingroup PkgClassification - -\brief Classifies a data set based on a set of features and a set of -labels. - -This class implements the core of the classification algorithm -\cgalCite{cgal:lm-clscm-12} (section 2). It uses a data set as input -and assigns each input item to a label among a set of user defined -labels. To achieve this classification, a set of local geometric -features are used, such as planarity, elevation or vertical -dispersion. In addition, the user must define a set of labels such as -building, ground or vegetation. - -Each pair of feature and label must be assigned an -[Feature::Effect](@ref CGAL::Classification::Feature::Effect) (for -example, vegetation has a low planarity and a high vertical -dispersion) and each feature must be assigned a weight. These -parameters can be set up by hand or by automatic training, provided a -small user-defined set of inlier is given for each classification -label. - -\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 -#elif defined(CGAL_LINKED_WITH_TBB) - typename ConcurrencyTag = CGAL::Parallel_tag> -#else - typename ConcurrencyTag = CGAL::Sequential_tag> -#endif - -class Classifier -{ - - -public: - typedef typename Classification::Label_handle Label_handle; - typedef typename Classification::Feature_handle Feature_handle; - typedef Classification::Feature::Effect Feature_effect; - - /// \cond SKIP_IN_MANUAL - typedef typename ItemMap::value_type Item; - -#ifdef CGAL_DO_NOT_USE_BOYKOV_KOLMOGOROV_MAXFLOW_SOFTWARE - typedef internal::Alpha_expansion_graph_cut_boost Alpha_expansion; -#else - typedef internal::Alpha_expansion_graph_cut_boykov_kolmogorov Alpha_expansion; -#endif - -protected: - -#ifdef CGAL_LINKED_WITH_TBB - class Run - { - Classifier& m_classifier; - const std::vector >& m_effect_table; - std::vector& m_assigned_label; - std::vector& m_confidence; - - - public: - - Run (Classifier& classifier, - const std::vector >& effect_table, - std::vector& assigned_label, - std::vector& confidence) - : m_classifier (classifier), m_effect_table (effect_table), - m_assigned_label (assigned_label), m_confidence (confidence) - { } - - void operator()(const tbb::blocked_range& r) const - { - for (std::size_t s = r.begin(); s != r.end(); ++ s) - { - std::size_t nb_class_best=0; - - double val_class_best = (std::numeric_limits::max)(); - std::vector values; - - for(std::size_t k = 0; k < m_effect_table.size(); ++ k) - { - double value = m_classifier.classification_value (k, s); - values.push_back (value); - - if(val_class_best > value) - { - val_class_best = value; - nb_class_best=k; - } - } - - m_assigned_label[s] = nb_class_best; - - std::sort (values.begin(), values.end()); - m_confidence[s] = values[1] - values[0]; - } - } - - }; - - class Run_with_local_smoothing_preprocessing - { - Classifier& m_classifier; - std::vector >& m_values; - const std::vector >& m_effect_table; - - public: - - Run_with_local_smoothing_preprocessing (Classifier& classifier, - std::vector >& values, - const std::vector >& effect_table) - : m_classifier (classifier), m_values(values), m_effect_table (effect_table) - { } - - void operator()(const tbb::blocked_range& r) const - { - for (std::size_t s = r.begin(); s != r.end(); ++ s) - for(std::size_t k = 0; k < m_effect_table.size(); ++ k) - m_values[k][s] = m_classifier.classification_value (k, s); - } - }; - - template - class Run_with_local_smoothing - { - Classifier& m_classifier; - const ItemRange& m_input; - ItemMap m_item_map; - const std::vector >& m_values; - const NeighborQuery& m_neighbor_query; - std::vector& m_assigned_label; - std::vector& m_confidence; - - public: - - Run_with_local_smoothing (Classifier& classifier, - const ItemRange& input, - ItemMap item_map, - const std::vector >& values, - const NeighborQuery& neighbor_query, - std::vector& assigned_label, - std::vector& confidence) - : m_classifier (classifier), m_input (input), m_item_map (item_map), - m_values(values), - m_neighbor_query (neighbor_query), - m_assigned_label (assigned_label), m_confidence (confidence) - { } - - void operator()(const tbb::blocked_range& r) const - { - for (std::size_t s = r.begin(); s != r.end(); ++ s) - { - std::vector neighbors; - m_neighbor_query (get (m_item_map, *(m_input.begin()+s)), std::back_inserter (neighbors)); - - std::vector mean (m_values.size(), 0.); - for (std::size_t n = 0; n < neighbors.size(); ++ n) - for (std::size_t j = 0; j < m_values.size(); ++ j) - mean[j] += m_values[j][neighbors[n]]; - - std::size_t nb_class_best=0; - double val_class_best = (std::numeric_limits::max)(); - for(std::size_t k = 0; k < mean.size(); ++ k) - { - mean[k] /= neighbors.size(); - if(val_class_best > mean[k]) - { - val_class_best = mean[k]; - nb_class_best = k; - } - } - - m_assigned_label[s] = nb_class_best; - - std::sort (mean.begin(), mean.end()); - m_confidence[s] = mean[1] - mean[0]; - } - } - - }; - - template - class Run_with_graphcut - { - Classifier& m_classifier; - const ItemRange& m_input; - ItemMap m_item_map; - const NeighborQuery& m_neighbor_query; - double m_weight; - const std::vector >& m_effect_table; - const std::vector >& m_indices; - const std::vector >& m_input_to_indices; - std::vector& m_assigned_label; - - public: - - Run_with_graphcut (Classifier& classifier, - const ItemRange& input, - ItemMap item_map, - const NeighborQuery& neighbor_query, - double weight, - const std::vector >& effect_table, - const std::vector >& indices, - const std::vector >& input_to_indices, - std::vector& assigned_label) - : m_classifier (classifier), m_input (input), m_item_map (item_map), - m_neighbor_query (neighbor_query), m_weight (weight), m_effect_table (effect_table), - m_indices (indices), m_input_to_indices (input_to_indices), - m_assigned_label (assigned_label) - { } - - void operator()(const tbb::blocked_range& r) const - { - for (std::size_t sub = r.begin(); sub != r.end(); ++ sub) - { - if (m_indices[sub].empty()) - continue; - - std::vector > edges; - std::vector edge_weights; - std::vector > probability_matrix - (m_effect_table.size(), std::vector(m_indices[sub].size(), 0.)); - std::vector assigned_label (m_indices[sub].size()); - - for (std::size_t j = 0; j < m_indices[sub].size(); ++ j) - { - std::size_t s = m_indices[sub][j]; - - std::vector neighbors; - - m_neighbor_query (get(m_item_map, *(m_input.begin()+s)), std::back_inserter (neighbors)); - - for (std::size_t i = 0; i < neighbors.size(); ++ i) - if (sub == m_input_to_indices[neighbors[i]].first - && j != m_input_to_indices[neighbors[i]].second) - { - edges.push_back (std::make_pair (j, m_input_to_indices[neighbors[i]].second)); - edge_weights.push_back (m_weight); - } - - std::size_t nb_class_best = 0; - double val_class_best = (std::numeric_limits::max)(); - for(std::size_t k = 0; k < m_effect_table.size(); ++ k) - { - double value = m_classifier.classification_value (k, s); - probability_matrix[k][j] = value; - - if(val_class_best > value) - { - val_class_best = value; - nb_class_best = k; - } - } - assigned_label[j] = nb_class_best; - } - - Alpha_expansion graphcut; - graphcut(edges, edge_weights, probability_matrix, assigned_label); - - for (std::size_t i = 0; i < assigned_label.size(); ++ i) - m_assigned_label[m_indices[sub][i]] = assigned_label[i]; - - } - } - - }; - - tbb::mutex m_mutex; - void mutex_lock() { m_mutex.lock(); } - void mutex_unlock() { m_mutex.unlock(); } -#else // CGAL_LINKED_WITH_TBB - void mutex_lock() { } - void mutex_unlock() { } -#endif // CGAL_LINKED_WITH_TBB - - const ItemRange& m_input; - ItemMap m_item_map; - - std::vector m_assigned_label; - std::vector m_confidence; - - std::vector m_labels; - std::vector m_features; - - std::vector > m_effect_table; - /// \endcond - -public: - - - /// \name Constructor - /// @{ - - /*! - \brief Initializes a classification object. - - \param input input range. - - \param item_map property map to access the input items. - */ - Classifier (const ItemRange& input, - ItemMap item_map) - : m_input (input), m_item_map (item_map) - { - } - - /// @} - - /// \cond SKIP_IN_MANUAL - virtual ~Classifier() { } - /// \endcond - - /// \name Features - /// @{ - - /*! - \brief Adds a feature. - - \tparam Feature type of the feature, inherited from - `Classification::Feature_base`. - - \tparam T types of the parameters of the feature's constructor - (with the exception of the first parameter that is always of type - `ItemRange&` and that is automatically passed by the classifier to - the feature's constructor). - - \param t parameters of the feature's constructor (with the - exception of the first parameter that is always the input item - range and that is automatically passed by the classifier to the - feature's constructor). - - \return a handle to the newly added feature. - */ -#if (!defined(CGAL_CFG_NO_CPP0X_VARIADIC_TEMPLATES) && !defined(CGAL_CFG_NO_CPP0X_RVALUE_REFERENCE)) || DOXYGEN_RUNNING - template - Feature_handle add_feature (T&& ... t) - { - Feature_handle fh (new Feature(m_input, std::forward(t)...)); - mutex_lock(); - m_features.push_back (fh); - mutex_unlock(); - return fh; - } -#else - template - Feature_handle add_feature () - { - Feature_handle fh (new Feature(m_input)); - mutex_lock(); - m_features.push_back (fh); - mutex_unlock(); - return fh; - } - template - Feature_handle add_feature (T1& t1) - { - Feature_handle fh (new Feature(m_input, t1)); - mutex_lock(); - m_features.push_back (fh); - mutex_unlock(); - return fh; - } - template - Feature_handle add_feature (T1& t1, T2& t2) - { - Feature_handle fh (new Feature(m_input, t1, t2)); - mutex_lock(); - m_features.push_back (fh); - mutex_unlock(); - return fh; - } - template - Feature_handle add_feature (T1& t1, T2& t2, T3& t3) - { - Feature_handle fh (new Feature(m_input, t1, t2, t3)); - mutex_lock(); - m_features.push_back (fh); - mutex_unlock(); - return fh; - } - template - Feature_handle add_feature (T1& t1, T2& t2, T3& t3, T4& t4) - { - Feature_handle fh (new Feature(m_input, t1, t2, t3, t4)); - mutex_lock(); - m_features.push_back (fh); - mutex_unlock(); - return fh; - } - template - Feature_handle add_feature (T1& t1, T2& t2, T3& t3, T4& t4, T5& t5) - { - Feature_handle fh (new Feature(m_input, t1, t2, t3, t4, t5)); - mutex_lock(); - m_features.push_back (fh); - mutex_unlock(); - return fh; - } -#endif - - /*! - \brief Removes a feature. - - \param feature the handle to feature type that must be removed. - - \return `true` if the feature was correctly removed, `false` if - its handle was not found. - */ - bool remove_feature (Feature_handle feature) - { - for (std::size_t i = 0; i < m_features.size(); ++ i) - if (m_features[i] == feature) - { - m_features.erase (m_features.begin() + i); - return true; - } - return false; - } - - /*! - \brief Returns how many features are defined. - */ - std::size_t number_of_features() const - { - return m_features.size(); - } - - - /*! - \brief Returns the \f$i^{th}\f$ feature. - */ - Feature_handle feature(std::size_t i) - { - return m_features[i]; - } - - /*! - \brief Removes all features. - */ - void clear_features () - { - m_features.clear(); - } - - /// \endcond - - /// @} - - - /// \name Labels - /// @{ - - /*! - \brief Adds a label. - - \param name name of the label. - - \return a handle to the newly added label. - */ - Label_handle add_label (const char* name) - { - Label_handle out (new Classification::Label (name)); - m_labels.push_back (out); - return out; - } - - /*! - \brief Removes a label. - - \param label the handle to the label that must be removed. - - \return `true` if the label was correctly removed, - `false` if its handle was not found. - */ - bool remove_label (Label_handle label) - { - std::size_t idx = (std::size_t)(-1); - for (std::size_t i = 0; i < m_labels.size(); ++ i) - if (m_labels[i] == label) - { - m_labels.erase (m_labels.begin() + i); - idx = i; - break; - } - if (idx == (std::size_t)(-1)) - return false; - std::cerr << idx << std::endl; - - for (std::size_t i = 0; i < m_assigned_label.size(); ++ i) - if (m_assigned_label[i] == (std::size_t)(-1)) - continue; - else if (m_assigned_label[i] > idx) - m_assigned_label[i] --; - else if (m_assigned_label[i] == idx) - m_assigned_label[i] = (std::size_t)(-1); - - return true; - } - - /*! - \brief Returns how many labels are defined. - */ - std::size_t number_of_labels () const - { - return m_labels.size(); - } - - /*! - \brief Returns the \f$i^{th}\f$ label. - */ - Label_handle label (std::size_t i) const - { - return m_labels[i]; - } - - - /*! - \brief Removes all labels. - */ - void clear_labels () - { - m_labels.clear(); - } - - /// @} - - /// \name Classification - /// @{ - - - /*! - \brief Runs the classification algorithm without any regularization. - - There is no relationship between items, the classification energy - is only minimized itemwise. This method is quick but produce - suboptimal results. - */ - void run() - { - prepare_classification (); - -#ifndef CGAL_LINKED_WITH_TBB - CGAL_static_assertion_msg (!(boost::is_convertible::value), - "Parallel_tag is enabled but TBB is unavailable."); -#else - if (boost::is_convertible::value) - { - Run f (*this, m_effect_table, m_assigned_label, m_confidence); - tbb::parallel_for(tbb::blocked_range(0, m_input.size ()), f); - } - else -#endif - { - - for (std::size_t s = 0; s < m_input.size(); s++) - { - - std::size_t nb_class_best=0; - - double val_class_best = (std::numeric_limits::max)(); - std::vector values; - - for(std::size_t k = 0; k < m_effect_table.size(); ++ k) - { - double value = classification_value (k, s); - values.push_back (value); - - if(val_class_best > value) - { - val_class_best = value; - nb_class_best=k; - } - } - - m_assigned_label[s] = nb_class_best; - - std::sort (values.begin(), values.end()); - m_confidence[s] = values[1] - values[0]; - } - } - } - - - /*! - \brief Runs the classification algorithm with a local smoothing. - - The computed classification energy is smoothed on a user defined - local neighborhood of items. This method is a compromise between - efficiency and reliability. - - \tparam NeighborQuery model of `NeighborQuery`. - \param neighbor_query used to access neighborhoods of items. - */ - template - void run_with_local_smoothing (const NeighborQuery& neighbor_query) - { - prepare_classification (); - - // data term initialisation - CGAL_CLASSIFICATION_CERR << "Labeling... "; - - std::vector > values - (m_labels.size(), - std::vector (m_input.size(), -1.)); - -#ifndef CGAL_LINKED_WITH_TBB - CGAL_static_assertion_msg (!(boost::is_convertible::value), - "Parallel_tag is enabled but TBB is unavailable."); -#else - if (boost::is_convertible::value) - { - Run_with_local_smoothing_preprocessing f1 - (*this, values, m_effect_table); - tbb::parallel_for(tbb::blocked_range(0, m_input.size ()), f1); - - Run_with_local_smoothing f2 - (*this, m_input, m_item_map, values, neighbor_query,m_assigned_label, m_confidence); - tbb::parallel_for(tbb::blocked_range(0, m_input.size ()), f2); - } - else -#endif - { - for (std::size_t s=0; s < m_input.size(); ++ s) - { - std::vector neighbors; - neighbor_query (get (m_item_map, *(m_input.begin()+s)), std::back_inserter (neighbors)); - - std::vector mean (values.size(), 0.); - for (std::size_t n = 0; n < neighbors.size(); ++ n) - { - if (values[0][neighbors[n]] < 0.) - for(std::size_t k = 0; k < m_effect_table.size(); ++ k) - { - values[k][neighbors[n]] = classification_value (k, neighbors[n]); - mean[k] += values[k][neighbors[n]]; - } - else - for (std::size_t j = 0; j < values.size(); ++ j) - mean[j] += values[j][neighbors[n]]; - } - - std::size_t nb_class_best=0; - double val_class_best = (std::numeric_limits::max)(); - for(std::size_t k = 0; k < mean.size(); ++ k) - { - mean[k] /= neighbors.size(); - if(val_class_best > mean[k]) - { - val_class_best = mean[k]; - nb_class_best = k; - } - } - - m_assigned_label[s] = nb_class_best; - - std::sort (mean.begin(), mean.end()); - m_confidence[s] = mean[1] - mean[0]; - } - } - } - - - /// \cond SKIP_IN_MANUAL - // (see below for the documented method) - template - void run_with_one_graphcut (const NeighborQuery& neighbor_query, - const double weight) - { - prepare_classification (); - - // data term initialisation -#ifdef CGAL_DO_NOT_USE_BOYKOV_KOLMOGOROV_MAXFLOW_SOFTWARE - CGAL_CLASSIFICATION_CERR << "Labeling using Boost with regularization weight " << weight << "... "; -#else - CGAL_CLASSIFICATION_CERR << "Labeling using Boyvok Kolmogorov with regularization weight " << weight << "... "; -#endif - - std::vector > edges; - std::vector edge_weights; - std::vector > probability_matrix - (m_effect_table.size(), std::vector(m_input.size(), 0.)); - std::vector(m_input.size()).swap(m_assigned_label); - - std::cerr << "Size of probability matrix = " << m_effect_table.size() * m_input.size() << std::endl; - std::cerr << "Size of assigned labels = " << m_assigned_label.size() << std::endl; - - for (std::size_t s = 0; s < m_input.size(); ++ s) - { - std::vector neighbors; - - neighbor_query (get(m_item_map, *(m_input.begin()+s)), std::back_inserter (neighbors)); - - for (std::size_t i = 0; i < neighbors.size(); ++ i) - if (s != neighbors[i]) - { - edges.push_back (std::make_pair (s, neighbors[i])); - edge_weights.push_back (weight); - } - - std::size_t nb_class_best = 0; - double val_class_best = (std::numeric_limits::max)(); - for(std::size_t k = 0; k < m_effect_table.size(); ++ k) - { - double value = classification_value (k, s); - probability_matrix[k][s] = value; - - if(val_class_best > value) - { - val_class_best = value; - nb_class_best = k; - } - } - m_assigned_label[s] = nb_class_best; - } - - Alpha_expansion graphcut; - graphcut(edges, edge_weights, probability_matrix, m_assigned_label); - std::cerr << ((double)(CGAL::Memory_sizer().virtual_size()) / 1073741824.) << " GB allocated" << std::endl; - } - /// \endcond - - /*! - \brief Runs the classification algorithm with a global - regularization based on a graphcut. - - The computed classification energy is globally regularized through - an alpha-expansion algorithm. This method is slow but provides - the user with good quality results. - - \tparam NeighborQuery model of `NeighborQuery`. - \param neighbor_query used to access neighborhoods of items. - \param weight weight of the regularization with respect to the - classification energy. Higher values produce more regularized - output but may result in a loss of details. - - \param min_number_of_subdivisions used to make computation - faster. The graphcut algorithm is applied separately to a certain - number of subdivisions of the input set. If `ConcurrencyTag` is - `CGAL::Parallel_tag`, the graphcuts are computed in parallel. The - default value (1) implies a unique and therefore sequential - graphcut algorithm. - - */ - template - void run_with_graphcut (const NeighborQuery& neighbor_query, - const double weight, - std::size_t min_number_of_subdivisions = 1) - { - if (min_number_of_subdivisions <= 1) - return run_with_one_graphcut (neighbor_query, weight); - - - prepare_classification (); - - // data term initialisation -#ifdef CGAL_DO_NOT_USE_BOYKOV_KOLMOGOROV_MAXFLOW_SOFTWARE - CGAL_CLASSIFICATION_CERR << "Labeling using Boost with regularization weight " << weight << "... "; -#else - CGAL_CLASSIFICATION_CERR << "Labeling using Boyvok Kolmogorov with regularization weight " << weight << "... "; -#endif - - CGAL::Bbox_3 bbox = CGAL::bbox_3 - (boost::make_transform_iterator (m_input.begin(), CGAL::Property_map_to_unary_function(m_item_map)), - boost::make_transform_iterator (m_input.end(), CGAL::Property_map_to_unary_function(m_item_map))); - - double Dx = bbox.xmax() - bbox.xmin(); - double Dy = bbox.ymax() - bbox.ymin(); - double A = Dx * Dy; - double a = A / min_number_of_subdivisions; - double l = std::sqrt(a); - std::size_t nb_x = std::size_t(Dx / l) + 1; - std::size_t nb_y = std::size_t((A / nb_x) / a) + 1; - std::size_t nb = nb_x * nb_y; - - std::vector bboxes; - bboxes.reserve(nb); - for (std::size_t x = 0; x < nb_x; ++ x) - for (std::size_t y = 0; y < nb_y; ++ y) - { - bboxes.push_back - (CGAL::Bbox_3 (bbox.xmin() + Dx * (x / double(nb_x)), - bbox.ymin() + Dy * (y / double(nb_y)), - bbox.zmin(), - bbox.xmin() + Dx * ((x+1) / double(nb_x)), - bbox.ymin() + Dy * ((y+1) / double(nb_y)), - bbox.zmax())); - } - - - std::cerr << "Number of divisions = " << nb_x * nb_y << std::endl; - std::cerr << " -> Size of division: " << Dx / nb_x << " " << Dy / nb_y << std::endl; - - std::vector > indices (nb); - std::vector > input_to_indices(m_input.size()); - - for (std::size_t s = 0; s < m_input.size(); ++ s) - { - CGAL::Bbox_3 b = get(m_item_map, *(m_input.begin() + s)).bbox(); - - for (std::size_t i = 0; i < bboxes.size(); ++ i) - if (CGAL::do_overlap (b, bboxes[i])) - { - input_to_indices[s] = std::make_pair (i, indices[i].size()); - indices[i].push_back (s); - break; - } - } - - std::vector(m_input.size()).swap(m_assigned_label); - -#ifndef CGAL_LINKED_WITH_TBB - CGAL_static_assertion_msg (!(boost::is_convertible::value), - "Parallel_tag is enabled but TBB is unavailable."); -#else - if (boost::is_convertible::value) - { - Run_with_graphcut f - (*this, m_input, m_item_map, neighbor_query, weight, m_effect_table, indices, input_to_indices, - m_assigned_label); - tbb::parallel_for(tbb::blocked_range(0, indices.size ()), f); - } - else -#endif - { - for (std::size_t sub = 0; sub < indices.size(); ++ sub) - { - if (indices[sub].empty()) - continue; - - CGAL_CLASSIFICATION_CERR << "Subset #" << sub << ": " - << indices[sub].size() << " points" << std::endl; - - std::vector > edges; - std::vector edge_weights; - std::vector > probability_matrix - (m_effect_table.size(), std::vector(indices[sub].size(), 0.)); - std::vector assigned_label (indices[sub].size()); - - for (std::size_t j = 0; j < indices[sub].size(); ++ j) - { - std::size_t s = indices[sub][j]; - - std::vector neighbors; - - neighbor_query (get(m_item_map, *(m_input.begin()+s)), std::back_inserter (neighbors)); - - for (std::size_t i = 0; i < neighbors.size(); ++ i) - if (sub == input_to_indices[neighbors[i]].first - && j != input_to_indices[neighbors[i]].second) - { - edges.push_back (std::make_pair (j, input_to_indices[neighbors[i]].second)); - edge_weights.push_back (weight); - } - - std::size_t nb_class_best = 0; - double val_class_best = (std::numeric_limits::max)(); - for(std::size_t k = 0; k < m_effect_table.size(); ++ k) - { - double value = classification_value (k, s); - probability_matrix[k][j] = value; - - if(val_class_best > value) - { - val_class_best = value; - nb_class_best = k; - } - } - assigned_label[j] = nb_class_best; - } - - Alpha_expansion graphcut; - graphcut(edges, edge_weights, probability_matrix, assigned_label); - - for (std::size_t i = 0; i < assigned_label.size(); ++ i) - m_assigned_label[indices[sub][i]] = assigned_label[i]; - } - } - - } - /// @} - - - - /// \name Output - /// @{ - - /*! - \brief Returns the value of the energy of `label` at the item at - position `index`. - */ - double energy_of (Label_handle label, std::size_t index) const - { - double out = 0.; - for (std::size_t i = 0; i < m_features.size(); ++ i) - { - if (m_features[i]->weight() == 0.) - continue; - - Feature_effect eff = label->feature_effect (m_features[i]); - - if (eff == Classification::Feature::FAVORING) - out += m_features[i]->favored (index); - else if (eff == Classification::Feature::PENALIZING) - out += m_features[i]->penalized (index); - else if (eff == Classification::Feature::NEUTRAL) - out += m_features[i]->ignored (index); - } - return out; - } - - /*! - \brief Returns the label of the item at position - `index`. - - \note If classification was not performed (using `run()`, - `run_with_local_smoothing()` or `run_with_graphcut()`), this - function always returns the default `Label_handle`. - */ - Label_handle label_of (std::size_t index) const - { - if (m_assigned_label.size() <= index - || m_assigned_label[index] == (std::size_t)(-1)) - { - return Label_handle(); - } - return m_labels[m_assigned_label[index]]; - } - - /// \cond SKIP_IN_MANUAL - void set_label_of (std::size_t index, Label_handle label) - { - if (index >= m_assigned_label.size()) - m_assigned_label.resize (index + 1, (std::size_t)(-1)); - for (std::size_t i = 0; i < m_labels.size(); ++ i) - if (m_labels[i] == label) - { - m_assigned_label[index] = i; - return; - } - m_assigned_label[index] = (std::size_t)(-1); - } - /// \endcond - - /*! - - \brief Returns the confidence of the label of the - item at position `index`. - - \note If classification was not performed (using `run()`, - `run_with_local_smoothing()` or `run_with_graphcut()`), this - function always returns 0. - - \return confidence ranging from 0 (not confident at all) to 1 - (very confident). - */ - double confidence_of (std::size_t index) const - { - if (m_confidence.size() <= index) - return 0.; - return m_confidence[index]; - } - - /// @} - - - - -protected: - - /// \cond SKIP_IN_MANUAL - void prepare_classification () - { - // Reset data structure - std::vector(m_input.size(), (std::size_t)(-1)).swap (m_assigned_label); - std::vector(m_input.size()).swap (m_confidence); - - m_effect_table = std::vector > - (m_labels.size(), std::vector (m_features.size(), - Classification::Feature::NEUTRAL)); - - for (std::size_t i = 0; i < m_effect_table.size (); ++ i) - for (std::size_t j = 0; j < m_effect_table[i].size (); ++ j) - m_effect_table[i][j] = m_labels[i]->feature_effect (m_features[j]); - - } - - double classification_value (const std::size_t& label, - const std::size_t& pt_index) const - { - double out = 0.; - for (std::size_t i = 0; i < m_effect_table[label].size(); ++ i) - { - if (m_features[i]->weight() == 0.) - continue; - if (m_effect_table[label][i] == Classification::Feature::FAVORING) - out += m_features[i]->favored (pt_index); - else if (m_effect_table[label][i] == Classification::Feature::PENALIZING) - out += m_features[i]->penalized (pt_index); - else if (m_effect_table[label][i] == Classification::Feature::NEUTRAL) - out += m_features[i]->ignored (pt_index); - } - return out; - } - /// \endcond -}; - - - - - -} // namespace CGAL - -#endif // CGAL_CLASSIFIER_H - diff --git a/Classification/include/CGAL/Point_set_classifier.h b/Classification/include/CGAL/Point_set_classifier.h deleted file mode 100644 index 5daee100ffd..00000000000 --- a/Classification/include/CGAL/Point_set_classifier.h +++ /dev/null @@ -1,1166 +0,0 @@ -// 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_POINT_SET_CLASSIFIER_H -#define CGAL_POINT_SET_CLASSIFIER_H - -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#include -#include -#include -#include -#include -#include - -#include -#include - -#ifdef CGAL_LINKED_WITH_TBB -#include -#include -#endif // CGAL_LINKED_WITH_TBB - - -namespace CGAL { - -/*! - \ingroup PkgClassification - - \brief Classifies a point set based on a set of features and a set - of labels. - - This class specializes `Classifier` to point sets. It takes care of - generating necessary data structures and automatically generate a - set of generic features. Features can be generated at multiple - scales to increase the reliability of the classification. - - \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 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). - \tparam DiagonalizeTraits model of `DiagonalizeTraits` used - for matrix diagonalization. - -*/ -template > -class Point_set_classifier : public Classifier -{ - -public: - typedef typename Geom_traits::Iso_cuboid_3 Iso_cuboid_3; - - /// \cond SKIP_IN_MANUAL - typedef Classifier Base; - typedef typename PointRange::const_iterator Iterator; - using Base::m_input; - using Base::m_item_map; - /// \endcond - - typedef Classification::Planimetric_grid - Planimetric_grid; - typedef Classification::Point_set_neighborhood - Neighborhood; - typedef Classification::Local_eigen_analysis - Local_eigen_analysis; - - /// \cond SKIP_IN_MANUAL - typedef Classification::Feature_handle Feature_handle; - typedef Classification::Label Label; - typedef Classification::Label_handle Label_handle; - - typedef typename Geom_traits::Point_3 Point; - - typedef Classification::Feature::Anisotropy - Anisotropy; - typedef Classification::Feature::Distance_to_plane - Distance_to_plane; - typedef Classification::Feature::Eigentropy - Eigentropy; - typedef Classification::Feature::Elevation - Elevation; - typedef Classification::Feature::Linearity - Linearity; - typedef Classification::Feature::Omnivariance - Omnivariance; - typedef Classification::Feature::Planarity - Planarity; - typedef Classification::Feature::Sphericity - Sphericity; - typedef Classification::Feature::Sum_eigenvalues - Sum_eigen; - typedef Classification::Feature::Surface_variation - Surface_variation; - typedef Classification::Feature::Vertical_dispersion - Dispersion; - typedef Classification::Feature::Verticality - Verticality; - typedef typename Classification::RGB_Color RGB_Color; - /// \endcond - -private: - - struct Scale - { - Neighborhood* neighborhood; - Planimetric_grid* grid; - Local_eigen_analysis* eigen; - double voxel_size; - std::vector features; - - Scale (const PointRange& input, PointMap point_map, - const Iso_cuboid_3& bbox, double voxel_size) - : voxel_size (voxel_size) - { - CGAL::Real_timer t; - t.start(); - if (voxel_size < 0.) - neighborhood = new Neighborhood (input, point_map); - else - neighborhood = new Neighborhood (input, point_map, voxel_size); - t.stop(); - - if (voxel_size < 0.) - CGAL_CLASSIFICATION_CERR << "Neighborhood computed in " << t.time() << " second(s)" << std::endl; - else - CGAL_CLASSIFICATION_CERR << "Neighborhood with voxel size " << voxel_size - << " computed in " << t.time() << " second(s)" << std::endl; - t.reset(); - t.start(); - - eigen = new Local_eigen_analysis (input, point_map, neighborhood->k_neighbor_query(6)); - double range = eigen->mean_range(); - if (this->voxel_size < 0) - this->voxel_size = range; - t.stop(); - CGAL_CLASSIFICATION_CERR << "Eigen values computed in " << t.time() << " second(s)" << std::endl; - CGAL_CLASSIFICATION_CERR << "Range = " << range << std::endl; - t.reset(); - t.start(); - - grid = new Planimetric_grid (input, point_map, bbox, this->voxel_size); - t.stop(); - CGAL_CLASSIFICATION_CERR << "Planimetric grid computed in " << t.time() << " second(s)" << std::endl; - t.reset(); - } - ~Scale() - { - delete neighborhood; - delete grid; - delete eigen; - } - - double grid_resolution() const { return voxel_size; } - double radius_neighbors() const { return voxel_size * 5; } - double radius_dtm() const { return voxel_size * 100; } - - }; - - Iso_cuboid_3 m_bbox; - std::vector m_scales; - -#ifdef CGAL_LINKED_WITH_TBB - tbb::mutex m_scale_mutex; - tbb::task_group* m_tasks; -#endif - - struct Feature_adder - { - mutable Point_set_classifier* classifier; - mutable Scale* scale; - - Feature_adder (Point_set_classifier* classifier, Scale* scale) - : classifier (classifier), scale (scale) { } - - virtual ~Feature_adder() { } - -#ifdef CGAL_LINKED_WITH_TBB - void mutex_lock() const { classifier->m_scale_mutex.lock(); } - void mutex_unlock() const { classifier->m_scale_mutex.unlock(); } -#else - void mutex_lock() const { } - void mutex_unlock() const { } -#endif - - virtual void operator()() const = 0; - }; - friend Feature_adder; - - - std::vector m_adders; - - -public: - - - /// \name Constructor - /// @{ - - /*! - \brief Initializes a classification object. - - \param input input range. - - \param point_map property map to access the input points. - */ - Point_set_classifier(const PointRange& input, PointMap point_map) : Base (input, point_map) - { - m_bbox = CGAL::bounding_box - (boost::make_transform_iterator (m_input.begin(), CGAL::Property_map_to_unary_function(m_item_map)), - boost::make_transform_iterator (m_input.end(), CGAL::Property_map_to_unary_function(m_item_map))); - } - - /// @} - - /// \cond SKIP_IN_MANUAL - virtual ~Point_set_classifier() - { - clear(); - } - /// \endcond - - /// \name Features - /// @{ - - - /*! - \brief Generate all possible features from an input range. - - The size of the smallest scale is automatically estimated and the - data structures needed (`Neighborhood`, `Planimetric_grid` and - `Local_eigen_analysis`) are computed at `nb_scales` recursively - larger scales. At each scale, the following features are - generated: - - - `CGAL::Classification::Feature::Anisotropy` - - `CGAL::Classification::Feature::Distance_to_plane` - - `CGAL::Classification::Feature::Eigentropy` - - `CGAL::Classification::Feature::Elevation` - - `CGAL::Classification::Feature::Linearity` - - `CGAL::Classification::Feature::Omnivariance` - - `CGAL::Classification::Feature::Planarity` - - `CGAL::Classification::Feature::Sphericity` - - `CGAL::Classification::Feature::Sum_eigenvalues` - - `CGAL::Classification::Feature::Surface_variation` - - `CGAL::Classification::Feature::Vertical_dispersion` based on eigenvalues - - If normal vectors are provided (if `VectorMap` is different from - `CGAL::Default`), the following feature is generated at each - scale: - - - `CGAL::Classification::Feature::Vertical_dispersion` based on normal vectors - - If colors are provided (if `ColorMap` is different from - `CGAL::Default`), the following features are generated at each - scale: - - - 9 features `CGAL::Classification::Feature::Hsv` on - channel 0 (hue) with mean ranging from 0° to 360° and standard - deviation of 22.5. - - - 5 features `CGAL::Classification::Feature::Hsv` on - channel 1 (saturation) with mean ranging from 0 to 100 and standard - deviation of 12.5. - - - 5 features `CGAL::Classification::Feature::Hsv` on channel 2 - (value) with mean ranging from 0 to 100 and standard deviation - of 12.5. - - If echo numbers are provided (if `EchoMap` is different from - `CGAL::Default`), the following feature is computed at each - scale: - - - `CGAL::Classification::Feature::Echo_scatter` - - \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`. - \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`. - \tparam EchoMap model of `ReadablePropertyMap` whose key type is - the value type of the iterator of `PointRange` and value type is - `std::size_t`. - \param nb_scales number of scales to compute. - \param normal_map property map to access the normal vectors of the input points (if any). - \param color_map property map to access the colors of the input points (if any). - \param echo_map property map to access the echo values of the input points (if any). - */ - template - void generate_features (std::size_t nb_scales, - VectorMap normal_map = VectorMap(), - ColorMap color_map = ColorMap(), - EchoMap echo_map = EchoMap()) - { - typedef typename Default::Get::type - Vmap; - typedef typename Default::Get::type - Cmap; - typedef typename Default::Get::type - Emap; - - generate_features_impl (nb_scales, - get_parameter(normal_map), - get_parameter(color_map), - get_parameter(echo_map)); - } - - /// @} - - /// \name Data Structures and Parameters - /// @{ - - /*! - \brief Returns the bounding box of the input point set. - */ - const Iso_cuboid_3& bbox() const { return m_bbox; } - /*! - \brief Returns the neighborhood structure at scale `scale`. - - \note `generate_features()` must have been called before calling - this method. - */ - const Neighborhood& neighborhood(std::size_t scale = 0) const { return (*m_scales[scale]->neighborhood); } - /*! - \brief Returns the planimetric grid structure at scale `scale`. - - \note `generate_features()` must have been called before calling - this method. - */ - const Planimetric_grid& grid(std::size_t scale = 0) const { return *(m_scales[scale]->grid); } - /*! - \brief Returns the local eigen analysis structure at scale `scale`. - - \note `generate_features()` must have been called before calling - this method. - */ - const Local_eigen_analysis& eigen(std::size_t scale = 0) const { return *(m_scales[scale]->eigen); } - /*! - \brief Returns the number of scales that were computed. - */ - std::size_t number_of_scales() const { return m_scales.size(); } - - /*! - \brief Returns the grid resolution at scale `scale`. This - resolution is the length and width of a cell of the - `Planimetric_grid` defined at this scale. - - \note `generate_features()` must have been called before calling - this method. - */ - double grid_resolution(std::size_t scale = 0) const { return m_scales[scale]->grid_resolution(); } - /*! - - \brief Returns the radius used for neighborhood queries at scale - `scale`. This radius is the smallest radius that is relevant from - a geometric point of view at this scale (that is to say that - encloses a few cells of `Planimetric_grid`). - - \note `generate_features()` must have been called before calling - this method. - */ - double radius_neighbors(std::size_t scale = 0) const { return m_scales[scale]->radius_neighbors(); } - /*! - \brief Returns the radius used for digital terrain modeling at - scale `scale`. This radius represents the minimum size of a - building at this scale. - - \note `generate_features()` must have been called before calling - this method. - */ - double radius_dtm(std::size_t scale = 0) const { return m_scales[scale]->radius_dtm(); } - - /*! - \brief Clears all computed data structures. - */ - void clear() - { - for (std::size_t i = 0; i < m_scales.size(); ++ i) - delete m_scales[i]; - m_scales.clear(); - - this->clear_labels(); - this->clear_features(); - } - - /// @} - - /// @} - - /// \cond SKIP_IN_MANUAL - void info() const - { - CGAL_CLASSIFICATION_CERR << m_scales.size() << " scale(s) used:" << std::endl; - for (std::size_t i = 0; i < m_scales.size(); ++ i) - { - std::size_t nb_useful = 0; - for (std::size_t j = 0; j < m_scales[i]->features.size(); ++ j) - if (m_scales[i]->features[j]->weight() != 0.) - nb_useful ++; - CGAL_CLASSIFICATION_CERR << " * scale " << i << " with size " << m_scales[i]->voxel_size - << ", " << nb_useful << " useful feature(s)"; - if (nb_useful != 0) CGAL_CLASSIFICATION_CERR << ":" << std::endl; - else CGAL_CLASSIFICATION_CERR << std::endl; - for (std::size_t j = 0; j < m_scales[i]->features.size(); ++ j) - if (m_scales[i]->features[j]->weight() != 0.) - CGAL_CLASSIFICATION_CERR << " - " << m_scales[i]->features[j]->name() - << " (weight = " << m_scales[i]->features[j]->weight() << ")" << std::endl; - } - } - - void generate_point_based_features () - { - - generate_multiscale_feature_variant_0 (); - generate_multiscale_feature_variant_0 (); - generate_multiscale_feature_variant_0 (); - generate_multiscale_feature_variant_0 (); - generate_multiscale_feature_variant_0 (); - generate_multiscale_feature_variant_0 (); - generate_multiscale_feature_variant_0 (); - generate_multiscale_feature_variant_0 (); - - generate_multiscale_feature_variant_1 (); - generate_multiscale_feature_variant_2 (); - generate_multiscale_feature_variant_3 (); - } - - template - void launch_feature_computation (FeatureAdder* adder) - { - m_adders.push_back (adder); - -#ifndef CGAL_LINKED_WITH_TBB - CGAL_static_assertion_msg (!(boost::is_convertible::value), - "Parallel_tag is enabled but TBB is unavailable."); -#else - if (boost::is_convertible::value) - { - m_tasks->run (*adder); - } - else -#endif - { - (*adder)(); - } - } - - template - struct Feature_adder_verticality : public Feature_adder - { - using Feature_adder::classifier; - using Feature_adder::scale; - VectorMap normal_map; - - // TODO! - Feature_adder_verticality (Point_set_classifier* classifier, VectorMap normal_map, Scale* scale) - : Feature_adder (classifier, scale), normal_map (normal_map) { } - - void operator() () const - { - Feature_handle result = classifier->template add_feature (normal_map); - this->mutex_lock(); - scale->features.push_back (result); - this->mutex_unlock(); - } - }; - - template - void generate_normal_based_features(VectorMap normal_map) - { - launch_feature_computation (new Feature_adder_verticality (this, normal_map, m_scales[0])); - } - - void generate_normal_based_features(const CGAL::Default_property_map&) - { - generate_multiscale_feature_variant_0 (); - } - - template - struct Feature_adder_color : public Feature_adder - { - typedef Classification::Feature::Hsv Hsv; - - using Feature_adder::classifier; - using Feature_adder::scale; - ColorMap color_map; - std::size_t channel; - double mean; - double sd; - - // TODO! - Feature_adder_color (Point_set_classifier* classifier, ColorMap color_map, Scale* scale, - std::size_t channel, double mean, double sd) - : Feature_adder (classifier, scale), color_map (color_map), - channel (channel), mean (mean), sd (sd) { } - - void operator() () const - { - Feature_handle result = classifier->template add_feature (color_map, channel, mean, sd); - this->mutex_lock(); - scale->features.push_back (result); - this->mutex_unlock(); - } - }; - - template - void generate_color_based_features(ColorMap color_map) - { - for (std::size_t i = 0; i <= 8; ++ i) - launch_feature_computation (new Feature_adder_color (this, color_map, m_scales[0], - 0, 45 * i, 22.5)); - - for (std::size_t i = 0; i <= 4; ++ i) - launch_feature_computation (new Feature_adder_color (this, color_map, m_scales[0], - 1, 25 * i, 12.5)); - - for (std::size_t i = 0; i <= 4; ++ i) - launch_feature_computation (new Feature_adder_color (this, color_map, m_scales[0], - 2, 25 * i, 12.5)); - } - - void generate_color_based_features(const CGAL::Default_property_map&) - { - } - - - template - struct Feature_adder_echo : public Feature_adder - { - typedef Classification::Feature::Echo_scatter Echo_scatter; - - using Feature_adder::classifier; - using Feature_adder::scale; - EchoMap echo_map; - - // TODO! - Feature_adder_echo (Point_set_classifier* classifier, EchoMap echo_map, Scale* scale) - : Feature_adder (classifier, scale), echo_map (echo_map) { } - - void operator() () const - { - Feature_handle result = classifier->template add_feature (echo_map, - *(scale->grid), - scale->grid_resolution(), - scale->radius_neighbors()); - this->mutex_lock(); - scale->features.push_back (result); - this->mutex_unlock(); - } - }; - - - template - void generate_echo_based_features(EchoMap echo_map) - { - for (std::size_t i = 0; i < m_scales.size(); ++ i) - launch_feature_computation (new Feature_adder_echo (this, echo_map, m_scales[i])); - } - - void generate_echo_based_features(const CGAL::Default_property_map&) - { - } - - void get_map_scale (std::map& map_scale) - { - for (std::size_t i = 0; i < m_scales.size(); ++ i) - for (std::size_t j = 0; j < m_scales[i]->features.size(); ++ j) - map_scale[m_scales[i]->features[j]] = i; - } - - std::size_t scale_of_feature (Feature_handle att) - { - for (std::size_t i = 0; i < m_scales.size(); ++ i) - for (std::size_t j = 0; j < m_scales[i]->features.size(); ++ j) - if (m_scales[i]->features[j] == att) - return i; - return (std::size_t)(-1); - } - - std::string name_att (Feature_handle att, - std::map& map_scale) - { - std::ostringstream oss; - oss << att->name() << "_" << map_scale[att]; - return oss.str(); - } - /// \endcond - - /// \name Input/Output - /// @{ - - /*! - \brief Saves the current configuration in the stream `output`. - - This allows to easily save and recover a specific classification - configuration, that is to say: - - - The size of the smallest scale - - The features and their respective weights - - The labels and the effects of the features on them - - The output file is written in an XML format that is readable by - the `load_configuration()` method. - */ - void save_configuration (std::ostream& output) - { - boost::property_tree::ptree tree; - - // tree.put("classification.parameters.grid_resolution", m_grid_resolution); - // tree.put("classification.parameters.radius_neighbors", m_radius_neighbors); - tree.put("classification.parameters.voxel_size", m_scales[0]->voxel_size); - - std::map map_scale; - get_map_scale (map_scale); - - for (std::size_t i = 0; i < this->number_of_features(); ++ i) - { - Feature_handle att = this->feature(i); - if (att->weight() == 0) - continue; - boost::property_tree::ptree ptr; - - ptr.put("id", name_att (att, map_scale)); - ptr.put("weight", att->weight()); - tree.add_child("classification.features.feature", ptr); - } - - - for (std::size_t i = 0; i < this->number_of_labels(); ++ i) - { - Label_handle label = this->label(i); - boost::property_tree::ptree ptr; - ptr.put("id", label->name()); - for (std::size_t j = 0; j < this->number_of_features(); ++ j) - { - Feature_handle att = this->feature(j); - if (att->weight() == 0) - continue; - boost::property_tree::ptree ptr2; - ptr2.put("id", name_att (att, map_scale)); - Classification::Feature::Effect effect = label->feature_effect(att); - if (effect == Classification::Feature::PENALIZING) - ptr2.put("effect", "penalized"); - else if (effect == Classification::Feature::NEUTRAL) - ptr2.put("effect", "neutral"); - else if (effect == Classification::Feature::FAVORING) - ptr2.put("effect", "favored"); - ptr.add_child("feature", ptr2); - } - tree.add_child("classification.labels.label", ptr); - } - - // Write property tree to XML file - boost::property_tree::xml_writer_settings settings(' ', 3); - boost::property_tree::write_xml(output, tree, settings); - } - - - /*! - \brief Loads a configuration from the stream `input`. - - All data structures, features and labels specified in the input - stream `input` are instantiated if possible (in particular, - property maps needed should be provided), similarly to what is - done in `generate_features()`. - - The input file should be in the XML format written by the - `save_configuration()` method. - - \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`. - \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`. - \tparam EchoMap model of `ReadablePropertyMap` whose key type is - the value type of the iterator of `PointRange` and value type is - `std::size_t`. - \param input input stream. - \param normal_map property map to access the normal vectors of the input points (if any). - \param color_map property map to access the colors of the input points (if any). - \param echo_map property map to access the echo values of the input points (if any). - */ - template - bool load_configuration (std::istream& input, - VectorMap normal_map = VectorMap(), - ColorMap color_map = ColorMap(), - EchoMap echo_map = EchoMap()) - { - typedef typename Default::Get::type - Vmap; - typedef typename Default::Get::type - Cmap; - typedef typename Default::Get::type - Emap; - - return load_configuration_impl (input, - get_parameter(normal_map), - get_parameter(color_map), - get_parameter(echo_map)); - } - - - - /*! - - \brief Writes the classified point set in a colored and labeled - PLY format in the stream `output`. - - The input points are written in a PLY format with the addition of - the following PLY properties: - - - a property `label` that indicates which label is - assigned to the point. The labels are indexed from 0 to N (the - correspondancy is given as comments in the PLY header). - - - 3 properties `red`, `green` and `blue` to associate each label - to a color (this is useful to visualize the classification in a - viewer that supports PLY colors). Colors are picked randomly. - */ - void write_classification_to_ply (std::ostream& output) - { - output << "ply" << std::endl - << "format ascii 1.0" << std::endl - << "comment Generated by the CGAL library www.cgal.org" << std::endl - << "element vertex " << m_input.size() << std::endl - << "property double x" << std::endl - << "property double y" << std::endl - << "property double z" << std::endl - << "property uchar red" << std::endl - << "property uchar green" << std::endl - << "property uchar blue" << std::endl - << "property int label" << std::endl; - - std::vector colors; - - std::map map_labels; - output << "comment label -1 is (unclassified)" << std::endl; - - for (std::size_t i = 0; i < this->number_of_labels(); ++ i) - { - map_labels.insert (std::make_pair (this->label(i), i)); - output << "comment label " << i << " is " << this->label(i)->name() << std::endl; - - RGB_Color c = {{ (unsigned char)(64 + rand() % 128), - (unsigned char)(64 + rand() % 128), - (unsigned char)(64 + rand() % 128) }}; - - if (this->label(i)->name() == "vegetation") - c = {{ 0, 255, 27 }}; - else if (this->label(i)->name() == "ground") - c = {{ 245, 180, 0 }}; - else if (this->label(i)->name() == "roof") - c = {{ 255, 0, 170 }}; - else if (this->label(i)->name() == "facade") - c = {{ 100, 0, 255 }}; - colors.push_back (c); - } - map_labels.insert (std::make_pair (Label_handle(), this->number_of_labels())); - - output << "end_header" << std::endl; - - std::size_t i = 0; - for (Iterator it = m_input.begin(); it != m_input.end(); ++ it) - { - Label_handle t = this->label_of(i); - std::size_t idx = map_labels[t]; - - if (idx == this->number_of_labels()) - output << get(m_item_map, *it) << " 0 0 0 -1" << std::endl; - else - output << get(m_item_map, *it) << " " - << (int)(colors[idx][0]) << " " - << (int)(colors[idx][1]) << " " - << (int)(colors[idx][2]) << " " - << idx << std::endl; - ++ i; - } - } - - /// @} - -private: - template - const T& get_parameter (const T& t) - { - return t; - } - - template - Default_property_map - get_parameter (const Default&) - { - return Default_property_map(); - } - - template - void generate_features_impl (std::size_t nb_scales, - VectorMap normal_map, - ColorMap color_map, - EchoMap echo_map) - { - CGAL::Real_timer t; t.start(); - - m_scales.reserve (nb_scales); - double voxel_size = - 1.; - - m_scales.push_back (new Scale (m_input, m_item_map, m_bbox, voxel_size)); - voxel_size = m_scales[0]->grid_resolution(); - - for (std::size_t i = 1; i < nb_scales; ++ i) - { - voxel_size *= 2; - m_scales.push_back (new Scale (m_input, m_item_map, m_bbox, voxel_size)); - } - t.stop(); - CGAL_CLASSIFICATION_CERR << "Scales computed in " << t.time() << " second(s)" << std::endl; - t.reset(); - - t.start(); - -#ifdef CGAL_LINKED_WITH_TBB - m_tasks = new tbb::task_group; -#endif - - generate_point_based_features (); - generate_normal_based_features (normal_map); - generate_color_based_features (color_map); - generate_echo_based_features (echo_map); - -#ifdef CGAL_LINKED_WITH_TBB - m_tasks->wait(); - delete m_tasks; -#endif - - t.stop(); - CGAL_CLASSIFICATION_CERR << "Features computed in " << t.time() << " second(s)" << std::endl; - for (std::size_t i = 0; i < m_adders.size(); ++ i) - delete m_adders[i]; - } - - template - struct Feature_adder_variant_0 : public Feature_adder - { - using Feature_adder::classifier; - using Feature_adder::scale; - - Feature_adder_variant_0 (Point_set_classifier* classifier, Scale* scale) - : Feature_adder (classifier, scale) { } - - void operator() () const - { - Feature_handle result = classifier->template add_feature (*(scale->eigen)); - this->mutex_lock(); - scale->features.push_back (result); - this->mutex_unlock(); - } - }; - - template - void generate_multiscale_feature_variant_0 () - { - for (std::size_t i = 0; i < m_scales.size(); ++ i) - launch_feature_computation (new Feature_adder_variant_0 (this, m_scales[i])); - } - - template - struct Feature_adder_variant_1 : public Feature_adder - { - using Feature_adder::classifier; - using Feature_adder::scale; - PointMap point_map; - - // TODO! - Feature_adder_variant_1 (Point_set_classifier* classifier, PointMap point_map, Scale* scale) - : Feature_adder (classifier, scale), point_map (point_map) { } - - void operator() () const - { - Feature_handle result = classifier->template add_feature (point_map, - *(scale->eigen)); - this->mutex_lock(); - scale->features.push_back (result); - this->mutex_unlock(); - } - }; - - template - void generate_multiscale_feature_variant_1 () - { - for (std::size_t i = 0; i < m_scales.size(); ++ i) - launch_feature_computation (new Feature_adder_variant_1 (this, m_item_map, m_scales[i])); - } - - template - struct Feature_adder_variant_2 : public Feature_adder - { - using Feature_adder::classifier; - using Feature_adder::scale; - PointMap point_map; - - // TODO! - Feature_adder_variant_2 (Point_set_classifier* classifier, PointMap point_map, Scale* scale) - : Feature_adder (classifier, scale), point_map (point_map) { } - - void operator() () const - { - Feature_handle result = classifier->template add_feature - (point_map, - *(scale->grid), - scale->grid_resolution(), - scale->radius_neighbors()); - this->mutex_lock(); - scale->features.push_back (result); - this->mutex_unlock(); - } - }; - - template - void generate_multiscale_feature_variant_2 () - { - for (std::size_t i = 0; i < m_scales.size(); ++ i) - launch_feature_computation (new Feature_adder_variant_2 (this, m_item_map, m_scales[i])); - } - - template - struct Feature_adder_variant_3 : public Feature_adder - { - using Feature_adder::classifier; - using Feature_adder::scale; - PointMap point_map; - - // TODO! - Feature_adder_variant_3 (Point_set_classifier* classifier, PointMap point_map, Scale* scale) - : Feature_adder (classifier, scale), point_map (point_map) { } - - void operator() () const - { - Feature_handle result = classifier->template add_feature (point_map, - *(scale->grid), - scale->grid_resolution(), - scale->radius_dtm()); - this->mutex_lock(); - scale->features.push_back (result); - this->mutex_unlock(); - } - }; - - template - void generate_multiscale_feature_variant_3 () - { - for (std::size_t i = 0; i < m_scales.size(); ++ i) - launch_feature_computation (new Feature_adder_variant_3 (this, m_item_map, m_scales[i])); - } - - template - bool load_configuration_impl (std::istream& input, - VectorMap normal_map, - ColorMap color_map, - EchoMap echo_map) - { - typedef Classification::Feature::Echo_scatter Echo_scatter; - typedef Classification::Feature::Hsv Hsv; - - clear(); - - boost::property_tree::ptree tree; - boost::property_tree::read_xml(input, tree); - - double voxel_size = tree.get("classification.parameters.voxel_size"); - - m_scales.push_back (new Scale (m_input, m_item_map, m_bbox, voxel_size)); - - CGAL::Real_timer t; - std::map att_map; - BOOST_FOREACH(boost::property_tree::ptree::value_type &v, tree.get_child("classification.features")) - { - std::string full_id = v.second.get("id"); - - std::vector splitted_id; - boost::split(splitted_id, full_id, boost::is_any_of("_")); - std::string id = splitted_id[0]; - for (std::size_t i = 1; i < splitted_id.size() - 1; ++ i) - id = id + "_" + splitted_id[i]; - std::size_t scale = std::atoi (splitted_id.back().c_str()); - - while (m_scales.size() <= scale) - { - voxel_size *= 2; - m_scales.push_back (new Scale (m_input, m_item_map, m_bbox, voxel_size)); - } - - double weight = v.second.get("weight"); - - // Generate the right feature if possible - if (id == "anisotropy") - this->template add_feature(*(m_scales[scale]->eigen)); - else if (id == "distance_to_plane") - this->template add_feature(m_item_map, *(m_scales[scale]->eigen)); - else if (id == "eigentropy") - this->template add_feature(*(m_scales[scale]->eigen)); - else if (id == "elevation") - { - t.start(); - this->template add_feature(m_item_map, - *(m_scales[scale]->grid), - m_scales[scale]->grid_resolution(), - m_scales[scale]->radius_dtm()); - t.stop(); - } - else if (id == "linearity") - this->template add_feature(*(m_scales[scale]->eigen)); - else if (id == "omnivariance") - this->template add_feature(*(m_scales[scale]->eigen)); - else if (id == "planarity") - this->template add_feature(*(m_scales[scale]->eigen)); - else if (id == "sphericity") - this->template add_feature(*(m_scales[scale]->eigen)); - else if (id == "sum_eigen") - this->template add_feature(*(m_scales[scale]->eigen)); - else if (id == "surface_variation") - this->template add_feature(*(m_scales[scale]->eigen)); - else if (id == "vertical_dispersion") - this->template add_feature(m_item_map, - *(m_scales[scale]->grid), - m_scales[scale]->grid_resolution(), - m_scales[scale]->radius_neighbors()); - else if (id == "verticality") - { - if (boost::is_convertible >::value) - this->template add_feature(*(m_scales[scale]->eigen)); - else - this->template add_feature(normal_map); - } - else if (id == "echo_scatter") - { - if (boost::is_convertible >::value) - { - CGAL_CLASSIFICATION_CERR << "Warning: echo_scatter required but no echo map given." << std::endl; - continue; - } - this->template add_feature(echo_map, *(m_scales[scale]->grid), - m_scales[scale]->grid_resolution(), - m_scales[scale]->radius_neighbors()); - } - else if (boost::starts_with(id.c_str(), "hue") - || boost::starts_with(id.c_str(), "saturation") - || boost::starts_with(id.c_str(), "value")) - { - if (boost::is_convertible >::value) - { - CGAL_CLASSIFICATION_CERR << "Warning: color feature required but no color map given." << std::endl; - continue; - } - if (boost::starts_with(id.c_str(), "hue")) - { - double value = boost::lexical_cast(id.c_str() + 4); - this->template add_feature(color_map, 0, value, 22.5); - } - else if (boost::starts_with(id.c_str(), "saturation")) - { - double value = boost::lexical_cast(id.c_str() + 11); - this->template add_feature(color_map, 1, value, 12.5); - } - else if (boost::starts_with(id.c_str(), "value")) - { - double value = boost::lexical_cast(id.c_str() + 6); - this->template add_feature(color_map, 2, value, 12.5); - } - } - else - { - CGAL_CLASSIFICATION_CERR << "Warning: unknown feature \"" << id << "\"" << std::endl; - continue; - } - - Feature_handle att = this->feature (this->number_of_features() - 1); - m_scales[scale]->features.push_back (att); - att->set_weight(weight); - att_map[full_id] = att; - } - - BOOST_FOREACH(boost::property_tree::ptree::value_type &v, tree.get_child("classification.labels")) - { - std::string label_id = v.second.get("id"); - - Label_handle new_label = this->add_label (label_id.c_str()); - - BOOST_FOREACH(boost::property_tree::ptree::value_type &v2, v.second) - { - if (v2.first == "id") - continue; - std::string att_id = v2.second.get("id"); - std::map::iterator it = att_map.find(att_id); - if (it == att_map.end()) - continue; - Feature_handle att = it->second; - std::string effect = v2.second.get("effect"); - if (effect == "penalized") - new_label->set_feature_effect (att, Classification::Feature::PENALIZING); - else if (effect == "neutral") - new_label->set_feature_effect (att, Classification::Feature::NEUTRAL); - else - new_label->set_feature_effect (att, Classification::Feature::FAVORING); - } - } - - return true; - } -}; - -} // namespace CGAL - - -#endif // CGAL_POINT_SET_CLASSIFIER_H diff --git a/Polyhedron/demo/Polyhedron/Plugins/Classification/Classification_plugin.cpp b/Polyhedron/demo/Polyhedron/Plugins/Classification/Classification_plugin.cpp index bf1e9998b3b..5896fa973f3 100644 --- a/Polyhedron/demo/Polyhedron/Plugins/Classification/Classification_plugin.cpp +++ b/Polyhedron/demo/Polyhedron/Plugins/Classification/Classification_plugin.cpp @@ -281,9 +281,9 @@ public Q_SLOTS: class_rows.clear(); // Add labels - for (std::size_t i = 0; i < classif->labels().size(); ++ i) - add_new_label (ClassRow (dock_widget, classif->labels()[i].first->name().c_str(), - classif->labels()[i].second)); + for (std::size_t i = 0; i < classif->number_of_labels(); ++ i) + add_new_label (ClassRow (dock_widget, classif->label(i)->name().c_str(), + classif->label_color(i))); // Enabled classif if features computed if (!(classif->features_computed())) @@ -339,6 +339,8 @@ public Q_SLOTS: else if (Scene_points_with_normal_item* points_item = qobject_cast(scene->item(scene->mainSelectionIndex()))) return create_from_item(points_item); + + return NULL; } @@ -798,17 +800,15 @@ public Q_SLOTS: if (att == Item_classification_base::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)); + ui_widget.feature_weight->setValue ((int)(1001. * 2. * std::atan(classif->weight(att)) / CGAL_PI)); - for (std::size_t i = 0; i < classif->labels().size(); ++ i) + for (std::size_t i = 0; i < classif->number_of_labels(); ++ i) { - CGAL::Classification::Feature::Effect - eff = classif->labels()[i].first->feature_effect(att); - if (eff == CGAL::Classification::Feature::PENALIZING) + CGAL::Classification::Sum_of_weighted_features_predicate::Effect + eff = classif->effect (classif->label(i), att); + if (eff == CGAL::Classification::Sum_of_weighted_features_predicate::PENALIZING) class_rows[i].effect->setCurrentIndex(0); - else if (eff == CGAL::Classification::Feature::NEUTRAL) + else if (eff == CGAL::Classification::Sum_of_weighted_features_predicate::NEUTRAL) class_rows[i].effect->setCurrentIndex(1); else class_rows[i].effect->setCurrentIndex(2); @@ -831,11 +831,10 @@ public Q_SLOTS: if (att == Item_classification_base::Feature_handle()) return; - att->set_weight(std::tan ((CGAL_PI/2.) * v / 1001.)); - // std::cerr << att->weight() << std::endl; + classif->set_weight(att, std::tan ((CGAL_PI/2.) * v / 1001.)); for (std::size_t i = 0; i < class_rows.size(); ++ i) - class_rows[i].effect->setEnabled(att->weight() != 0.); + class_rows[i].effect->setEnabled(classif->weight(att) != 0.); } void on_effect_changed (int v) @@ -860,23 +859,19 @@ public Q_SLOTS: // std::cerr << att->id() << " is "; if (v == 0) { - classif->labels()[i].first->set_feature_effect - (att, CGAL::Classification::Feature::PENALIZING); - // std::cerr << " penalized for "; + classif->set_effect(classif->label(i), + att, CGAL::Classification::Sum_of_weighted_features_predicate::PENALIZING); } else if (v == 1) { - classif->labels()[i].first->set_feature_effect - (att, CGAL::Classification::Feature::NEUTRAL); - // std::cerr << " neutral for "; + classif->set_effect(classif->label(i), + att, CGAL::Classification::Sum_of_weighted_features_predicate::NEUTRAL); } else { - classif->labels()[i].first->set_feature_effect - (att, CGAL::Classification::Feature::FAVORING); - // std::cerr << " favored for "; + classif->set_effect(classif->label(i), + att, CGAL::Classification::Sum_of_weighted_features_predicate::FAVORING); } - // std::cerr << classif->labels()[i].first->id() << std::endl; break; } } diff --git a/Polyhedron/demo/Polyhedron/Plugins/Classification/Item_classification_base.h b/Polyhedron/demo/Polyhedron/Plugins/Classification/Item_classification_base.h index 0d27833da7e..6a324e9d61d 100644 --- a/Polyhedron/demo/Polyhedron/Plugins/Classification/Item_classification_base.h +++ b/Polyhedron/demo/Polyhedron/Plugins/Classification/Item_classification_base.h @@ -5,15 +5,19 @@ #include -#include -#include +#include +#include +#include class Item_classification_base { public: typedef CGAL::Classification::Label_handle Label_handle; typedef CGAL::Classification::Feature_handle Feature_handle; - + typedef CGAL::Classification::Label_set Label_set; + typedef CGAL::Classification::Feature_set Feature_set; + typedef CGAL::Classification::Sum_of_weighted_features_predicate Predicate; + public: Item_classification_base() { } @@ -23,47 +27,110 @@ public: virtual void erase_item() = 0; virtual void compute_features () = 0; - virtual bool features_computed() const = 0; - virtual std::size_t number_of_features() const = 0; - virtual Feature_handle feature(std::size_t i) = 0; - virtual void add_new_label (const char* name, const QColor& color) = 0; - virtual void remove_label (const char* name) = 0; - virtual void add_selection_to_training_set (const char* name) = 0; virtual void reset_training_sets() = 0; + virtual void validate_selection () = 0; virtual void train() = 0; virtual bool run (int method) = 0; virtual void change_color (int index) = 0; - virtual void fill_display_combo_box (QComboBox* cb, QComboBox* cb1) const = 0; virtual void generate_one_item_per_label(std::vector& items, const char* name) const = 0; virtual bool write_output(std::ostream& out) = 0; - virtual void save_config(const char* filename) = 0; - virtual void load_config(const char* filename) = 0; + + bool features_computed() const { return (m_features.size() != 0); } + std::size_t number_of_features() const { return m_features.size(); } + Feature_handle feature(std::size_t i) { return m_features[i]; } + double weight (Feature_handle f) const { return m_predicate->weight(f); } + void set_weight (Feature_handle f, double w) const { m_predicate->set_weight(f,w); } + Predicate::Effect effect (Label_handle l, Feature_handle f) const { return m_predicate->effect(l,f); } + void set_effect (Label_handle l, Feature_handle f, Predicate::Effect e) + { m_predicate->set_effect (l, f, e); } + + void add_new_label (const char* name, const QColor& color) + { + m_labels.add(name); + m_label_colors.push_back (color); + delete m_predicate; + m_predicate = new Predicate (m_labels, m_features); + } + void remove_label (const char* name) + { + for (std::size_t i = 0; i < m_labels.size(); ++ i) + if (m_labels[i]->name() == name) + { + m_labels.remove(m_labels[i]); + m_label_colors.erase (m_label_colors.begin() + i); + break; + } + delete m_predicate; + m_predicate = new Predicate (m_labels, m_features); + } + std::size_t number_of_labels() const { return m_labels.size(); } + Label_handle label(std::size_t i) { return m_labels[i]; } + + void fill_display_combo_box (QComboBox* cb, QComboBox* cb1) const + { + for (std::size_t i = 0; i < m_features.size(); ++ i) + { + std::ostringstream oss; + oss << "Feature " << m_features[i]->name(); + cb->addItem (oss.str().c_str()); + cb1->addItem (oss.str().c_str()); + } + } + + void save_config(const char* filename) + { + if (m_features.size() == 0) + { + std::cerr << "Error: features not computed" << std::endl; + return; + } + + std::ofstream f (filename); + m_predicate->save_configuration (f); + } + void load_config(const char* filename) + { + if (m_features.size() == 0) + { + std::cerr << "Error: features not computed" << std::endl; + return; + } + + std::ifstream f (filename); + m_predicate->load_configuration (f, true); + } std::size_t& nb_scales() { return m_nb_scales; } std::size_t& number_of_trials() { return m_nb_trials; } double& smoothing() { return m_smoothing; } std::size_t& subdivisions() { return m_subdivisions; } - std::vector >& labels() { return m_labels; } - void change_label_color (const char* name, const QColor& color) + const QColor& label_color(std::size_t i) const { return m_label_colors[i]; } + std::size_t get_label (const char* name) { for (std::size_t i = 0; i < m_labels.size(); ++ i) - if (m_labels[i].first->name() == name) - { - m_labels[i].second = color; - break; - } + if (m_labels[i]->name() == name) + return i; + return std::size_t(-1); + } + void change_label_color (const char* name, const QColor& color) + { + m_label_colors[get_label(name)] = color; } protected: + Label_set m_labels; + Feature_set m_features; + std::vector m_label_colors; + Predicate* m_predicate; + std::size_t m_nb_scales; - std::vector > m_labels; std::size_t m_nb_trials; double m_smoothing; std::size_t m_subdivisions; diff --git a/Polyhedron/demo/Polyhedron/Plugins/Classification/Point_set_item_classification.cpp b/Polyhedron/demo/Polyhedron/Plugins/Classification/Point_set_item_classification.cpp index 1d09be9a6db..84a5888583d 100644 --- a/Polyhedron/demo/Polyhedron/Plugins/Classification/Point_set_item_classification.cpp +++ b/Polyhedron/demo/Polyhedron/Plugins/Classification/Point_set_item_classification.cpp @@ -11,10 +11,15 @@ #include #include +#ifdef CGAL_LINKED_WITH_TBB +typedef CGAL::Parallel_tag Concurrency_tag; +#else +typedef CGAL::Sequential_tag Concurrency_tag; +#endif + Point_set_item_classification::Point_set_item_classification(Scene_points_with_normal_item* points) : m_points (points), - m_psc (NULL), - m_trainer (NULL) + m_generator (NULL) { m_nb_scales = 5; m_index_color = 1; @@ -24,32 +29,35 @@ Point_set_item_classification::Point_set_item_classification(Scene_points_with_n reset_indices(); - m_psc = new PSC(*(m_points->point_set()), m_points->point_set()->point_map()); - backup_existing_colors_and_add_new(); + m_training = m_points->point_set()->add_property_map("training", std::size_t(-1)).first; + m_classif = m_points->point_set()->add_property_map("label", std::size_t(-1)).first; + + m_labels.add("ground"); + m_labels.add("vegetation"); + m_labels.add("roof"); + m_labels.add("facade"); - 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))); + m_label_colors.push_back (QColor(245, 180, 0)); + m_label_colors.push_back (QColor(0, 255, 27)); + m_label_colors.push_back (QColor(255, 0, 170)); + m_label_colors.push_back (QColor(100, 0, 255)); + m_predicate = new Predicate (m_labels, m_features); } Point_set_item_classification::~Point_set_item_classification() { - if (m_psc != NULL) - delete m_psc; - if (m_trainer != NULL) - delete m_trainer; + if (m_predicate != NULL) + delete m_predicate; + if (m_generator != NULL) + delete m_generator; if (m_points != NULL) - reset_colors(); + { + reset_colors(); + m_points->point_set()->remove_property_map(m_training); + } } @@ -105,23 +113,23 @@ void Point_set_item_classification::reset_colors() // Write point set to .PLY file bool Point_set_item_classification::write_output(std::ostream& stream) { - if (m_psc->number_of_features() == 0) + if (m_features.size() == 0) return false; reset_indices(); stream.precision (std::numeric_limits::digits10 + 2); - std::vector 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); - } + // std::vector 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); + // m_psc->write_classification_to_ply (stream); return true; } @@ -159,18 +167,14 @@ void Point_set_item_classification::change_color (int index) } else if (index_color == 1) // classif { - std::map 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); + std::size_t c = m_classif[*it]; - if (c != Label_handle()) - color = map_colors[c]; + if (c != std::size_t(-1)) + color = m_label_colors[c]; m_red[*it] = color.red(); m_green[*it] = color.green(); @@ -179,19 +183,16 @@ void Point_set_item_classification::change_color (int index) } else if (index_color == 2) // training { - std::map 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); + std::size_t c = m_training[*it]; + std::size_t c2 = m_classif[*it]; + + if (c != std::size_t(-1)) + color = m_label_colors[c]; - if (c != Label_handle()) - color = map_colors[c]; double div = 1; if (c != c2) div = 2; @@ -203,17 +204,22 @@ void Point_set_item_classification::change_color (int index) } else { - Feature_handle att = m_psc->feature(index_color - 3); - double weight = att->weight(); - att->set_weight(att->max); + Feature_handle feature = m_features[index_color - 3]; + + double max = 0.; + for (Point_set::const_iterator it = m_points->point_set()->begin(); + it != m_points->point_set()->first_selected(); ++ it) + if (feature->value(*it) > max) + max = feature->value(*it); + for (Point_set::const_iterator it = m_points->point_set()->begin(); it != m_points->point_set()->first_selected(); ++ it) { - m_red[*it] = (unsigned char)(ramp.r(att->normalized(*it)) * 255); - m_green[*it] = (unsigned char)(ramp.g(att->normalized(*it)) * 255); - m_blue[*it] = (unsigned char)(ramp.b(att->normalized(*it)) * 255); + double v = std::max (0., feature->value(*it) / max); + m_red[*it] = (unsigned char)(ramp.r(v) * 255); + m_green[*it] = (unsigned char)(ramp.g(v) * 255); + m_blue[*it] = (unsigned char)(ramp.b(v) * 255); } - att->set_weight(weight); } for (Point_set::const_iterator it = m_points->point_set()->first_selected(); @@ -250,13 +256,14 @@ void Point_set_item_classification::reset_indices () void Point_set_item_classification::compute_features () { CGAL_assertion (!(m_points->point_set()->empty())); - CGAL_assertion (m_psc != NULL); - m_psc->clear_features(); + + if (m_generator != NULL) + delete m_generator; + reset_indices(); std::cerr << "Computing features with " << m_nb_scales << " scale(s)" << std::endl; - if (m_psc->number_of_features() != 0) - m_psc->clear(); + m_features.clear(); bool normals = m_points->point_set()->has_normal_map(); bool colors = (m_color != Point_set::Property_map()); @@ -265,60 +272,93 @@ void Point_set_item_classification::compute_features () boost::tie (echo_map, echo) = m_points->point_set()->template property_map("echo"); if (!normals && !colors && !echo) - m_psc->generate_features (m_nb_scales); + m_generator = new Generator (m_features, m_nb_scales, *(m_points->point_set()), m_points->point_set()->point_map()); else if (!normals && !colors && echo) - m_psc->generate_features (m_nb_scales, CGAL::Default(), CGAL::Default(), echo_map); + m_generator = new Generator (m_features, m_nb_scales, *(m_points->point_set()), m_points->point_set()->point_map(), + CGAL::Default(), CGAL::Default(), echo_map); else if (!normals && colors && !echo) - m_psc->generate_features (m_nb_scales, CGAL::Default(), m_color); + m_generator = new Generator (m_features, m_nb_scales, *(m_points->point_set()), m_points->point_set()->point_map(), + CGAL::Default(), m_color); else if (!normals && colors && echo) - m_psc->generate_features (m_nb_scales, CGAL::Default(), m_color, echo_map); + m_generator = new Generator (m_features, m_nb_scales, *(m_points->point_set()), m_points->point_set()->point_map(), + CGAL::Default(), m_color, echo_map); else if (normals && !colors && !echo) - m_psc->generate_features (m_nb_scales, m_points->point_set()->normal_map()); + m_generator = new Generator (m_features, m_nb_scales, *(m_points->point_set()), m_points->point_set()->point_map(), + 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); + m_generator = new Generator (m_features, m_nb_scales, *(m_points->point_set()), m_points->point_set()->point_map(), + 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(), m_color); + m_generator = new Generator (m_features, m_nb_scales, *(m_points->point_set()), m_points->point_set()->point_map(), + m_points->point_set()->normal_map(), m_color); else - m_psc->generate_features (m_nb_scales, m_points->point_set()->normal_map(), m_color, echo_map); + m_generator = new Generator (m_features, m_nb_scales, *(m_points->point_set()), m_points->point_set()->point_map(), + m_points->point_set()->normal_map(), m_color, echo_map); + delete m_predicate; + m_predicate = new Predicate (m_labels, m_features); + std::cerr << "Features = " << m_features.size() << std::endl; } void Point_set_item_classification::train() { - if (m_psc->number_of_features() == 0) + if (m_features.size() == 0) { std::cerr << "Error: features not computed" << std::endl; return; } reset_indices(); + + std::vector indices (m_points->point_set()->size(), std::size_t(-1)); + + for (Point_set::const_iterator it = m_points->point_set()->begin(); + it != m_points->point_set()->first_selected(); ++ it) + indices[*it] = m_training[*it]; + + m_predicate->train(indices, m_nb_trials); + CGAL::Classification::classify (*(m_points->point_set()), + m_labels, *m_predicate, + indices); + for (Point_set::const_iterator it = m_points->point_set()->begin(); + it != m_points->point_set()->first_selected(); ++ it) + m_classif[*it] = indices[*it]; - m_trainer->train(m_nb_trials); - m_psc->run(); - m_psc->info(); if (m_index_color == 1 || m_index_color == 2) - change_color (m_index_color); + change_color (m_index_color); } bool Point_set_item_classification::run (int method) { - if (m_psc->number_of_features() == 0) + if (m_features.size() == 0) { std::cerr << "Error: features not computed" << std::endl; return false; } reset_indices(); + std::vector indices; + if (method == 0) - { - m_psc->run(); - m_psc->info(); - } + CGAL::Classification::classify (*(m_points->point_set()), + m_labels, *m_predicate, + indices); else if (method == 1) - m_psc->run_with_local_smoothing (m_psc->neighborhood().range_neighbor_query(m_psc->radius_neighbors())); + CGAL::Classification::classify_with_local_smoothing + (*(m_points->point_set()), m_points->point_set()->point_map(), m_labels, *m_predicate, + m_generator->neighborhood().range_neighbor_query(m_generator->radius_neighbors()), + indices); else if (method == 2) - m_psc->run_with_graphcut (m_psc->neighborhood().k_neighbor_query(12), m_smoothing, m_subdivisions); + CGAL::Classification::classify_with_graphcut + (*(m_points->point_set()), m_points->point_set()->point_map(), + m_points->point_set()->point_map(), m_labels, *m_predicate, + m_generator->neighborhood().k_neighbor_query(12), + m_smoothing, m_subdivisions, indices); + + for (Point_set::const_iterator it = m_points->point_set()->begin(); + it != m_points->point_set()->first_selected(); ++ it) + m_classif[*it] = indices[*it]; if (m_index_color == 1 || m_index_color == 2) change_color (m_index_color); diff --git a/Polyhedron/demo/Polyhedron/Plugins/Classification/Point_set_item_classification.h b/Polyhedron/demo/Polyhedron/Plugins/Classification/Point_set_item_classification.h index 0c2d2ad331f..6de2b599b1f 100644 --- a/Polyhedron/demo/Polyhedron/Plugins/Classification/Point_set_item_classification.h +++ b/Polyhedron/demo/Polyhedron/Plugins/Classification/Point_set_item_classification.h @@ -5,16 +5,6 @@ #define CGAL_CLASSIFICATION_VERBOSE #include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include #include "Scene_points_with_normal_item.h" #include "Item_classification_base.h" @@ -22,6 +12,8 @@ #include "Kernel_type.h" #include "Point_set_3.h" +#include + #include @@ -36,16 +28,8 @@ class Point_set_item_classification : public Item_classification_base typedef Point_set::Point_map Point_map; typedef Point_set::Vector_map Vector_map; - typedef CGAL::Point_set_classifier PSC; - typedef CGAL::Classification::Trainer Trainer; - typedef CGAL::Classification::Label_handle Label_handle; - typedef CGAL::Classification::Feature_handle Feature_handle; - typedef CGAL::Classification::Feature::Vertical_dispersion Dispersion; - typedef CGAL::Classification::Feature::Elevation Elevation; - typedef CGAL::Classification::Feature::Verticality Verticality; - typedef CGAL::Classification::Feature::Distance_to_plane Distance_to_plane; - - + typedef CGAL::Classification::Point_set_feature_generator Generator; + public: Point_set_item_classification(Scene_points_with_normal_item* points); @@ -55,37 +39,17 @@ class Point_set_item_classification : public Item_classification_base void erase_item() { m_points = NULL; } void compute_features (); - bool features_computed() const { return (m_psc->number_of_features() != 0); } - 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; - } - } - void add_selection_to_training_set (const char* name) { - Label_handle label = get_label (name); + std::size_t 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_training[*it] = label; + m_classif[*it] = label; + } m_points->resetSelection(); if (m_index_color == 1 || m_index_color == 2) @@ -93,16 +57,15 @@ class Point_set_item_classification : public Item_classification_base } void reset_training_sets() { - m_trainer->reset_inlier_sets(); + for (Point_set::const_iterator it = m_points->point_set()->begin(); + it != m_points->point_set()->end(); ++ it) + m_training[*it] = std::size_t(-1); } 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_training[*it] = m_classif[*it]; m_points->resetSelection(); if (m_index_color == 1 || m_index_color == 2) @@ -112,114 +75,34 @@ class Point_set_item_classification : public Item_classification_base bool run (int method); void change_color (int index); - void fill_display_combo_box (QComboBox* cb, QComboBox* cb1) const - { - 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 generate_one_item_per_label(std::vector& items, const char* name) const { - std::map map_labels; + std::vector points_item + (m_labels.size(), NULL); for (std::size_t i = 0; i < m_labels.size(); ++ i) { - Scene_points_with_normal_item* new_item = new Scene_points_with_normal_item; - new_item->setName (QString("%1 (%2)").arg(name).arg(m_labels[i].first->name().c_str())); - new_item->setColor (m_labels[i].second); - map_labels[m_labels[i].first] = new_item; - items.push_back (new_item); + points_item[i] = new Scene_points_with_normal_item; + points_item[i]->setName (QString("%1 (%2)").arg(name).arg(m_labels[i]->name().c_str())); + points_item[i]->setColor (m_label_colors[i]); + items.push_back (points_item[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()) - map_labels[c]->point_set()->insert (m_points->point_set()->point(*it)); + std::size_t c = m_classif[*it]; + if (c != std::size_t(-1)) + points_item[c]->point_set()->insert (m_points->point_set()->point(*it)); } } bool write_output(std::ostream& out); - 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_color != Point_set::Property_map()); - Point_set::Property_map echo_map; - bool echo; - boost::tie (echo_map, echo) = m_points->point_set()->template property_map("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(), m_color); - else if (!normals && colors && echo) - m_psc->load_configuration (f, CGAL::Default(), m_color, 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(), m_color); - else - m_psc->load_configuration (f, m_points->point_set()->normal_map(), m_color, echo_map); - - std::vector > 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); - } - int real_index_color() const; void reset_indices(); void backup_existing_colors_and_add_new(); void reset_colors(); - 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(); - } private: @@ -229,9 +112,10 @@ class Point_set_item_classification : public Item_classification_base Point_set::Property_map m_green; Point_set::Property_map m_blue; Point_set::Property_map m_color; - - PSC* m_psc; - Trainer* m_trainer; + Point_set::Property_map m_training; + Point_set::Property_map m_classif; + + Generator* m_generator; int m_index_color;