1 #if defined (_MSC_VER) && !defined (_WIN64)
2 #pragma warning(disable:4244) // boost::number_distance::distance()
3                               // converts 64 to 32 bits integers
4 #endif
5 
6 #include "Cluster_classification.h"
7 #include "Color_ramp.h"
8 
9 #include <CGAL/Delaunay_triangulation_3.h>
10 #include <CGAL/Delaunay_triangulation_cell_base_3.h>
11 #include <CGAL/Triangulation_vertex_base_with_info_3.h>
12 #include <CGAL/Timer.h>
13 #include <CGAL/Memory_sizer.h>
14 
15 #include <QLineEdit>
16 
17 #include <CGAL/Three/Viewer_interface.h>
18 
19 #include <set>
20 #include <stack>
21 #include <algorithm>
22 #include <boost/array.hpp>
23 
Cluster_classification(Scene_points_with_normal_item * points)24 Cluster_classification::Cluster_classification(Scene_points_with_normal_item* points)
25   : m_points (points)
26   , m_input_is_las (false)
27 {
28   m_index_color = 1;
29 
30   reset_indices();
31 
32   backup_existing_colors_and_add_new();
33 
34   bool cluster_found = false;
35   boost::tie (m_cluster_id, cluster_found) = m_points->point_set()->property_map<int>("shape");
36   if (!cluster_found)
37   {
38     std::cerr << "Error! Cluster not found!" << std::endl;
39     abort();
40   }
41 
42   CGAL::Classification::create_clusters_from_indices (*(m_points->point_set()),
43                                                       m_points->point_set()->point_map(),
44                                                       m_cluster_id,
45                                                       m_clusters);
46 
47   std::cerr << m_clusters.size() << " cluster(s) found" << std::endl;
48 
49   bool training_found = false;
50   boost::tie (m_training, training_found) = m_points->point_set()->add_property_map<int>("training", -1);
51   bool classif_found = false;
52   boost::tie (m_classif, classif_found) = m_points->point_set()->add_property_map<int>("label", -1);
53 
54   training_found = !training_found; // add_property_map returns false if
55   classif_found = !classif_found;   // property was already there
56 
57   bool las_found = false;
58 
59   if (!classif_found)
60   {
61     Point_set::Property_map<unsigned char> las_classif;
62     boost::tie (las_classif, las_found) = m_points->point_set()->property_map<unsigned char>("classification");
63     if (las_found)
64     {
65       m_input_is_las = true;
66       for (Point_set::const_iterator it = m_points->point_set()->begin();
67            it != m_points->point_set()->first_selected(); ++ it)
68       {
69         unsigned char uc = las_classif[*it];
70         m_classif[*it] = int(uc);
71         if (!training_found)
72           m_training[*it] = int(uc);
73       }
74       m_points->point_set()->remove_property_map (las_classif);
75       classif_found = true;
76       training_found = true;
77     }
78   }
79 
80   if (training_found || classif_found)
81   {
82     std::vector<int> used_indices;
83 
84     for (Point_set::const_iterator it = m_points->point_set()->begin();
85          it != m_points->point_set()->first_selected(); ++ it)
86     {
87       if (training_found)
88       {
89         int l = m_training[*it];
90         if (l >= 0)
91         {
92           if (std::size_t(l) >= used_indices.size())
93             used_indices.resize(std::size_t(l + 1), -1);
94           used_indices[std::size_t(l)] = 0;
95         }
96       }
97       if (classif_found)
98       {
99         int l = m_classif[*it];
100         if (l >= 0)
101         {
102           if (std::size_t(l) >= used_indices.size())
103             used_indices.resize(std::size_t(l + 1), -1);
104           used_indices[std::size_t(l)] = 0;
105         }
106       }
107     }
108 
109     // map indices to filtered indices
110     int current_idx = 0;
111     for (std::size_t i = 0; i < used_indices.size(); ++ i)
112     {
113       if (las_found && (i < 2))
114       {
115         used_indices[i] = -1;
116         continue;
117       }
118       if (used_indices[i] == -1)
119         continue;
120 
121       used_indices[i] = current_idx;
122       ++ current_idx;
123     }
124 
125 
126     for (Point_set::const_iterator it = m_points->point_set()->begin();
127          it != m_points->point_set()->first_selected(); ++ it)
128     {
129       int c = m_cluster_id[*it];
130 
131       if (training_found)
132       {
133         if (std::size_t(current_idx) != used_indices.size()) // Empty indices -> reorder indices in point set
134         {
135           if (las_found && (m_training[*it] == 0 || m_training[*it] == 1)) // Unclassified class in LAS
136             m_training[*it] = -1;
137           else if (m_training[*it] != -1)
138             m_training[*it] = used_indices[std::size_t(m_training[*it])];
139         }
140         if (c != -1 && m_training[*it] != -1)
141           m_clusters[c].training() = m_training[*it];
142       }
143       if (classif_found)
144       {
145         if (std::size_t(current_idx) != used_indices.size()) // Empty indices -> reorder indices in point set
146         {
147           if (las_found && (m_classif[*it] == 0 || m_classif[*it] == 1)) // Unclassified class in LAS
148             m_classif[*it] = -1;
149           else if (m_classif[*it] != -1)
150             m_classif[*it] = used_indices[std::size_t(m_classif[*it])];
151         }
152         if (c != -1 && m_classif[*it] != -1)
153           m_clusters[c].label() = m_classif[*it];
154       }
155     }
156 
157     std::map<int, std::string> label_names;
158     if (las_found) // Use LAS standard
159     {
160       // label_names.insert (std::make_pair (0, std::string("never_clfied")));
161       // label_names.insert (std::make_pair (1, std::string("unclassified")));
162       label_names.insert (std::make_pair (2, std::string("ground")));
163       label_names.insert (std::make_pair (3, std::string("low_veget")));
164       label_names.insert (std::make_pair (4, std::string("med_veget")));
165       label_names.insert (std::make_pair (5, std::string("high_veget")));
166       label_names.insert (std::make_pair (6, std::string("building")));
167       label_names.insert (std::make_pair (7, std::string("noise")));
168       label_names.insert (std::make_pair (8, std::string("reserved")));
169       label_names.insert (std::make_pair (9, std::string("water")));
170       label_names.insert (std::make_pair (10, std::string("rail")));
171       label_names.insert (std::make_pair (11, std::string("road_surface")));
172       label_names.insert (std::make_pair (12, std::string("reserved_2")));
173       label_names.insert (std::make_pair (13, std::string("wire_guard")));
174       label_names.insert (std::make_pair (14, std::string("wire_conduct")));
175       label_names.insert (std::make_pair (15, std::string("trans_tower")));
176       label_names.insert (std::make_pair (16, std::string("wire_connect")));
177       label_names.insert (std::make_pair (17, std::string("bridge_deck")));
178       label_names.insert (std::make_pair (18, std::string("high_noise")));
179     }
180     else // Try to recover label names from PLY comments
181     {
182 
183       const std::string& comments = m_points->comments();
184       std::istringstream stream (comments);
185       std::string line;
186       while (getline(stream, line))
187       {
188         std::stringstream iss (line);
189         std::string tag;
190         if (iss >> tag && tag == "label")
191         {
192           int idx;
193           std::string name;
194           if (iss >> idx >> name)
195             label_names.insert (std::make_pair (idx, name));
196         }
197       }
198     }
199 
200     for (std::size_t i = 0; i < used_indices.size(); ++ i)
201     {
202       if (used_indices[i] == -1)
203         continue;
204 
205       Label_handle new_label;
206       std::map<int, std::string>::iterator found
207         = label_names.find (int(i));
208       if (found != label_names.end())
209         new_label = m_labels.add(found->second.c_str());
210       else
211       {
212         std::ostringstream oss;
213         oss << "label_" << i;
214         new_label = m_labels.add(oss.str().c_str());
215       }
216     }
217   }
218   else
219   {
220     m_labels.add("ground");
221     m_labels.add("vegetation");
222     m_labels.add("roof");
223     m_labels.add("facade");
224   }
225 
226   update_comments_of_point_set_item();
227 
228   m_sowf = new Sum_of_weighted_features (m_labels, m_features);
229   m_ethz = nullptr;
230 #ifdef CGAL_LINKED_WITH_OPENCV
231   m_random_forest = nullptr;
232 #endif
233 
234   // Compute neighborhood
235 #if 0
236   typedef CGAL::Triangulation_vertex_base_with_info_3<int, Kernel>         Vb;
237   typedef CGAL::Delaunay_triangulation_cell_base_3<Kernel>                 Cb;
238   typedef CGAL::Triangulation_data_structure_3<Vb, Cb>                     Tds;
239   typedef CGAL::Delaunay_triangulation_3<Kernel, Tds>                      Delaunay;
240 
241   Delaunay dt (boost::make_transform_iterator
242                (m_points->point_set()->begin(),
243                 Point_set_with_cluster_info (m_points->point_set(),
244                                              m_cluster_id)),
245                 boost::make_transform_iterator
246                (m_points->point_set()->end(),
247                 Point_set_with_cluster_info (m_points->point_set(),
248                                              m_cluster_id)));
249 
250   std::set<std::pair<int, int> > adjacencies;
251 
252   for (Delaunay::Finite_edges_iterator it = dt.finite_edges_begin();
253        it != dt.finite_edges_end(); ++ it)
254   {
255     Delaunay::Vertex_handle v0 = it->first->vertex (it->second);
256     Delaunay::Vertex_handle v1 = it->first->vertex (it->third);
257     int a = v0->info();
258     int b = v1->info();
259 
260     if (a == -1 || b == -1)
261       continue;
262 
263     if (a > b)
264       std::swap (a, b);
265 
266     adjacencies.insert (std::make_pair (a, b));
267   }
268 
269   for (std::set<std::pair<int, int> >::iterator it = adjacencies.begin();
270        it != adjacencies.end(); ++ it)
271   {
272     m_clusters[std::size_t(it->first)].neighbors->push_back (std::size_t(it->second));
273     m_clusters[std::size_t(it->second)].neighbors->push_back (std::size_t(it->first));
274   }
275 #endif
276 }
277 
278 
~Cluster_classification()279 Cluster_classification::~Cluster_classification()
280 {
281   if (m_sowf != nullptr)
282     delete m_sowf;
283   if (m_ethz != nullptr)
284     delete m_ethz;
285 #ifdef CGAL_LINKED_WITH_OPENCV
286   if (m_random_forest != nullptr)
287     delete m_random_forest;
288 #endif
289 #ifdef CGAL_LINKED_WITH_TENSORFLOW
290   if (m_neural_network != nullptr)
291     delete m_neural_network;
292 #endif
293   if (m_points != nullptr)
294   {
295     for (Point_set::const_iterator it = m_points->point_set()->begin();
296          it != m_points->point_set()->first_selected(); ++ it)
297     {
298       int c = m_cluster_id[*it];
299       if (c != -1)
300       {
301         m_training[*it] = m_clusters[c].training();
302         m_classif[*it] = m_clusters[c].label();
303       }
304       else
305       {
306         m_training[*it] = -1;
307         m_classif[*it] = -1;
308       }
309     }
310 
311     // For LAS saving, convert classification info in the LAS standard
312 //    if (m_input_is_las)
313     {
314       Point_set::Property_map<unsigned char> las_classif
315         = m_points->point_set()->add_property_map<unsigned char>("classification", 0).first;
316 
317       std::vector<unsigned char> label_indices;
318 
319       unsigned char custom = 19;
320       for (std::size_t i = 0; i < m_labels.size(); ++ i)
321       {
322         if (m_labels[i]->name() == "ground")
323           label_indices.push_back (2);
324         else if (m_labels[i]->name() == "low_veget")
325           label_indices.push_back (3);
326         else if (m_labels[i]->name() == "med_veget" || m_labels[i]->name() == "vegetation")
327           label_indices.push_back (4);
328         else if (m_labels[i]->name() == "high_veget")
329           label_indices.push_back (5);
330         else if (m_labels[i]->name() == "building" || m_labels[i]->name() == "roof")
331           label_indices.push_back (6);
332         else if (m_labels[i]->name() == "noise")
333           label_indices.push_back (7);
334         else if (m_labels[i]->name() == "reserved" || m_labels[i]->name() == "facade")
335           label_indices.push_back (8);
336         else if (m_labels[i]->name() == "water")
337           label_indices.push_back (9);
338         else if (m_labels[i]->name() == "rail")
339           label_indices.push_back (10);
340         else if (m_labels[i]->name() == "road_surface")
341           label_indices.push_back (11);
342         else if (m_labels[i]->name() == "reserved_2")
343           label_indices.push_back (12);
344         else if (m_labels[i]->name() == "wire_guard")
345           label_indices.push_back (13);
346         else if (m_labels[i]->name() == "wire_conduct")
347           label_indices.push_back (14);
348         else if (m_labels[i]->name() == "trans_tower")
349           label_indices.push_back (15);
350         else if (m_labels[i]->name() == "wire_connect")
351           label_indices.push_back (16);
352         else if (m_labels[i]->name() == "bridge_deck")
353           label_indices.push_back (17);
354         else if (m_labels[i]->name() == "high_noise")
355           label_indices.push_back (18);
356         else
357           label_indices.push_back (custom ++);
358       }
359 
360       for (Point_set::const_iterator it = m_points->point_set()->begin();
361            it != m_points->point_set()->end(); ++ it)
362       {
363         int c = m_classif[*it];
364         unsigned char lc = 1; // unclassified in LAS standard
365         if (c != -1)
366           lc = label_indices[std::size_t(c)];
367 
368         las_classif[*it] = lc;
369 
370         int t = m_training[*it];
371         unsigned char lt = 1; // unclassified in LAS standard
372         if (t != -1)
373           lt = label_indices[std::size_t(t)];
374 
375         m_training[*it] = int(lt);
376       }
377 
378       m_points->point_set()->remove_property_map (m_classif);
379     }
380 
381 
382     reset_colors();
383     erase_item();
384   }
385 }
386 
387 
backup_existing_colors_and_add_new()388 void Cluster_classification::backup_existing_colors_and_add_new()
389 {
390   if (m_points->point_set()->has_colors())
391   {
392     m_color = m_points->point_set()->add_property_map<CGAL::IO::Color>("real_color").first;
393     for (Point_set::const_iterator it = m_points->point_set()->begin();
394          it != m_points->point_set()->first_selected(); ++ it)
395       m_color[*it] = CGAL::IO::Color ((unsigned char)(255 * m_points->point_set()->red(*it)),
396                                   (unsigned char)(255 * m_points->point_set()->green(*it)),
397                                   (unsigned char)(255 * m_points->point_set()->blue(*it)));
398 
399     m_points->point_set()->remove_colors();
400   }
401 
402   m_points->point_set()->add_colors();
403 }
404 
reset_colors()405 void Cluster_classification::reset_colors()
406 {
407   if (m_color == Point_set::Property_map<CGAL::IO::Color>())
408     m_points->point_set()->remove_colors();
409   else
410   {
411     for (Point_set::const_iterator it = m_points->point_set()->begin();
412          it != m_points->point_set()->first_selected(); ++ it)
413       m_points->point_set()->set_color(*it, m_color[*it]);
414 
415     m_points->point_set()->remove_property_map(m_color);
416   }
417 }
418 
change_color(int index,float * vmin,float * vmax)419 void Cluster_classification::change_color (int index, float* vmin, float* vmax)
420 {
421   m_index_color = index;
422 
423   int index_color = real_index_color();
424 
425   // Colors
426   static Color_ramp ramp;
427   ramp.build_rainbow();
428   reset_indices();
429 
430   if (index_color == -1) // item color
431     m_points->point_set()->remove_colors();
432   else
433   {
434     if (!m_points->point_set()->has_colors())
435       m_points->point_set()->add_colors();
436 
437     if (index_color == 0) // real colors
438     {
439 
440       for (Point_set::const_iterator it = m_points->point_set()->begin();
441            it != m_points->point_set()->first_selected(); ++ it)
442         m_points->point_set()->set_color(*it, m_color[*it]);
443     }
444     else if (index_color == 1) // classif
445     {
446       for (Point_set::const_iterator it = m_points->point_set()->begin();
447            it != m_points->point_set()->first_selected(); ++ it)
448       {
449         QColor color (0, 0, 0);
450         int cid = m_cluster_id[*it];
451         if (cid != -1)
452         {
453           std::size_t c = m_clusters[cid].label();
454 
455           if (c != std::size_t(-1))
456             color = label_qcolor (m_labels[std::size_t(c)]);
457         }
458 
459         m_points->point_set()->set_color(*it, color);
460       }
461     }
462     else if (index_color == 2) // training
463     {
464       for (Point_set::const_iterator it = m_points->point_set()->begin();
465            it != m_points->point_set()->first_selected(); ++ it)
466       {
467         QColor color (0, 0, 0);
468         int cid = m_cluster_id[*it];
469         float div = 1;
470 
471         if (cid != -1)
472         {
473           int c = m_clusters[cid].training();
474           int c2 = m_clusters[cid].label();
475 
476           if (c != -1)
477             color = label_qcolor (m_labels[std::size_t(c)]);
478 
479           if (c != c2)
480             div = 2;
481         }
482         m_points->point_set()->set_color(*it, color.red() / div, color.green() / div, color.blue() / div);
483       }
484     }
485     else if (index_color == 3) // clusters
486     {
487       for (Point_set::const_iterator it = m_points->point_set()->begin();
488            it != m_points->point_set()->first_selected(); ++ it)
489       {
490         int cid = m_cluster_id[*it];
491 
492         if (cid != -1)
493         {
494           srand(cid);
495           m_points->point_set()->set_color(*it, 64 + rand() % 192, 64 + rand() % 192, 64 + rand() % 192);
496         }
497         else
498         {
499           m_points->point_set()->set_color(*it);
500         }
501       }
502     }
503     else
504     {
505       std::size_t corrected_index = index_color - 4;
506       if (corrected_index < m_labels.size()) // Display label probabilities
507       {
508         if (m_label_probabilities.size() <= corrected_index ||
509             m_label_probabilities[corrected_index].size() != m_clusters.size())
510         {
511           for (Point_set::const_iterator it = m_points->point_set()->begin();
512                it != m_points->point_set()->first_selected(); ++ it)
513             m_points->point_set()->set_color(*it);
514         }
515         else
516         {
517           for (Point_set::const_iterator it = m_points->point_set()->begin();
518                it != m_points->point_set()->first_selected(); ++ it)
519           {
520             int cid = m_cluster_id[*it];
521             if (cid != -1)
522             {
523               float v = (std::max) (0.f, (std::min)(1.f, m_label_probabilities[corrected_index][cid]));
524               m_points->point_set()->set_color(*it, ramp.r(v) * 255, ramp.g(v) * 255, ramp.b(v) * 255);
525             }
526             else
527               m_points->point_set()->set_color(*it);
528           }
529         }
530       }
531       else
532       {
533         corrected_index -= m_labels.size();
534         if (corrected_index >= m_features.size())
535         {
536           std::cerr << "Error: trying to access feature " << corrected_index << " out of " << m_features.size() << std::endl;
537           return;
538         }
539 
540         Feature_handle feature = m_features[corrected_index];
541 
542         float min = (std::numeric_limits<float>::max)();
543         float max = -(std::numeric_limits<float>::max)();
544 
545         if (vmin != nullptr && vmax != nullptr
546             && *vmin != std::numeric_limits<float>::infinity()
547             && *vmax != std::numeric_limits<float>::infinity())
548         {
549           min = *vmin;
550           max = *vmax;
551         }
552         else
553         {
554           for (Point_set::const_iterator it = m_points->point_set()->begin();
555                it != m_points->point_set()->first_selected(); ++ it)
556           {
557             int cid = m_cluster_id[*it];
558             if (cid != -1)
559             {
560               if (feature->value(cid) > max)
561                 max = feature->value(cid);
562               if (feature->value(cid) < min)
563                 min = feature->value(cid);
564             }
565           }
566         }
567 
568         for (Point_set::const_iterator it = m_points->point_set()->begin();
569              it != m_points->point_set()->first_selected(); ++ it)
570         {
571           int cid = m_cluster_id[*it];
572           if (cid != -1)
573           {
574             float v = (feature->value(cid) - min) / (max - min);
575             if (v < 0.f) v = 0.f;
576             if (v > 1.f) v = 1.f;
577 
578             m_points->point_set()->set_color(*it, ramp.r(v) * 255, ramp.g(v) * 255, ramp.b(v) * 255);
579           }
580           else
581             m_points->point_set()->set_color(*it);
582         }
583 
584         if (vmin != nullptr && vmax != nullptr)
585         {
586           *vmin = min;
587           *vmax = max;
588         }
589       }
590     }
591   }
592 
593   for (Point_set::const_iterator it = m_points->point_set()->first_selected();
594        it != m_points->point_set()->end(); ++ it)
595     m_points->point_set()->set_color(*it, 255, 0, 0);
596 }
597 
real_index_color() const598 int Cluster_classification::real_index_color() const
599 {
600   int out = m_index_color;
601 
602   if (out == 0 && m_color == Point_set::Property_map<CGAL::IO::Color>())
603     out = -1;
604   return out;
605 }
606 
reset_indices()607 void Cluster_classification::reset_indices ()
608 {
609   Point_set::Property_map<Point_set::Index> indices
610     = m_points->point_set()->property_map<Point_set::Index>("index").first;
611 
612   m_points->point_set()->unselect_all();
613   Point_set::Index idx;
614   ++ idx;
615   for (std::size_t i = 0; i < m_points->point_set()->size(); ++ i)
616     *(indices.begin() + i) = idx ++;
617 }
618 
compute_features(std::size_t nb_scales,float voxel_size)619 void Cluster_classification::compute_features (std::size_t nb_scales, float voxel_size)
620 {
621   CGAL_assertion (!(m_points->point_set()->empty()));
622 
623   reset_indices();
624 
625   std::cerr << "Computing pointwise features with " << nb_scales << " scale(s) and ";
626   if (voxel_size == -1)
627     std::cerr << "automatic voxel size" << std::endl;
628   else
629     std::cerr << "voxel size = " << voxel_size << std::endl;
630 
631   m_features.clear();
632 
633   Point_set::Vector_map normal_map;
634   bool normals = m_points->point_set()->has_normal_map();
635   if (normals)
636     normal_map = m_points->point_set()->normal_map();
637 
638   bool colors = (m_color != Point_set::Property_map<CGAL::IO::Color>());
639 
640   Point_set::Property_map<boost::uint8_t> echo_map;
641   bool echo;
642   boost::tie (echo_map, echo) = m_points->point_set()->template property_map<boost::uint8_t>("echo");
643   if (!echo)
644     boost::tie (echo_map, echo) = m_points->point_set()->template property_map<boost::uint8_t>("number_of_returns");
645 
646   Feature_set pointwise_features;
647 
648   Generator generator (*(m_points->point_set()), m_points->point_set()->point_map(), nb_scales, voxel_size);
649 
650   CGAL::Real_timer t;
651   t.start();
652 
653 #ifdef CGAL_LINKED_WITH_TBB
654   pointwise_features.begin_parallel_additions();
655 #endif
656 
657   generator.generate_point_based_features(pointwise_features);
658   if (normals)
659     generator.generate_normal_based_features (pointwise_features, normal_map);
660   if (colors)
661     generator.generate_color_based_features (pointwise_features, m_color);
662   if (echo)
663     generator.generate_echo_based_features (pointwise_features, echo_map);
664 
665 #ifdef CGAL_LINKED_WITH_TBB
666   pointwise_features.end_parallel_additions();
667 #endif
668 
669   add_remaining_point_set_properties_as_features(pointwise_features);
670 
671   t.stop();
672   std::cerr << pointwise_features.size() << " feature(s) computed in " << t.time() << " second(s)" << std::endl;
673   t.reset();
674 
675   std::cerr << "Computing cluster features" << std::endl;
676   t.start();
677 
678 #ifdef CGAL_LINKED_WITH_TBB
679   m_features.begin_parallel_additions();
680 #endif
681 
682   for (std::size_t i = 0; i < pointwise_features.size(); ++ i)
683   {
684     m_features.template add<Mean_of_feature> (m_clusters,
685                                               pointwise_features[i]);
686   }
687 
688 #ifdef CGAL_LINKED_WITH_TBB
689   m_features.end_parallel_additions();
690   m_features.begin_parallel_additions();
691 #endif
692 
693   for (std::size_t i = 0; i < pointwise_features.size(); ++ i)
694   {
695     m_features.template add<Variance_of_feature> (m_clusters,
696                                                   pointwise_features[i],
697                                                   m_features[i]);
698   }
699 
700   add_cluster_features();
701 
702 #ifdef CGAL_LINKED_WITH_TBB
703   m_features.end_parallel_additions();
704 #endif
705 
706   t.stop();
707   std::cerr << m_features.size() << " feature(s) computed in " << t.time() << " second(s)" << std::endl;
708 
709 
710   delete m_sowf;
711   m_sowf = new Sum_of_weighted_features (m_labels, m_features);
712   if (m_ethz != nullptr)
713   {
714     delete m_ethz;
715     m_ethz = nullptr;
716   }
717 #ifdef CGAL_LINKED_WITH_OPENCV
718   if (m_random_forest != nullptr)
719   {
720     delete m_random_forest;
721     m_random_forest = nullptr;
722   }
723 #endif
724 
725   std::cerr << "Features = " << m_features.size() << std::endl;
726 }
727 
select_random_region()728 void Cluster_classification::select_random_region()
729 {
730   std::size_t c = rand() % m_clusters.size();
731   std::set<Point_set::Index> in_cluster;
732   for (std::size_t i = 0; i < m_clusters[c].size(); ++ i)
733     in_cluster.insert (m_clusters[c].index(i));
734 
735   std::vector<std::size_t> selected;
736   std::vector<std::size_t> unselected;
737 
738   for (Point_set::const_iterator it = m_points->point_set()->begin();
739        it != m_points->point_set()->end(); ++ it)
740     if (in_cluster.find (*it) != in_cluster.end())
741       selected.push_back (*it);
742     else
743       unselected.push_back (*it);
744 
745   for (std::size_t i = 0; i < unselected.size(); ++ i)
746     *(m_points->point_set()->begin() + i) = unselected[i];
747   for (std::size_t i = 0; i < selected.size(); ++ i)
748     *(m_points->point_set()->begin() + (unselected.size() + i)) = selected[i];
749 
750   m_points->point_set()->set_first_selected
751     (m_points->point_set()->begin() + unselected.size());
752 
753 }
754 
add_remaining_point_set_properties_as_features(Feature_set & feature_set)755 void Cluster_classification::add_remaining_point_set_properties_as_features(Feature_set& feature_set)
756 {
757   const std::vector<std::string>& prop = m_points->point_set()->base().properties();
758 
759   for (std::size_t i = 0; i < prop.size(); ++ i)
760   {
761     if (prop[i] == "index" ||
762         prop[i] == "point" ||
763         prop[i] == "normal" ||
764         prop[i] == "echo" ||
765         prop[i] == "number_of_returns" ||
766         prop[i] == "training" ||
767         prop[i] == "label" ||
768         prop[i] == "classification" ||
769         prop[i] == "scan_direction_flag" ||
770         prop[i] == "real_color" ||
771         prop[i] == "shape" ||
772         prop[i] == "red" || prop[i] == "green" || prop[i] == "blue" ||
773         prop[i] == "r" || prop[i] == "g" || prop[i] == "b")
774       continue;
775 
776     if (try_adding_simple_feature<boost::int8_t>(feature_set, prop[i]))
777       continue;
778     if (try_adding_simple_feature<boost::uint8_t>(feature_set, prop[i]))
779       continue;
780     if (try_adding_simple_feature<boost::int16_t>(feature_set, prop[i]))
781       continue;
782     if (try_adding_simple_feature<boost::uint16_t>(feature_set, prop[i]))
783       continue;
784     if (try_adding_simple_feature<boost::int32_t>(feature_set, prop[i]))
785       continue;
786     if (try_adding_simple_feature<boost::uint32_t>(feature_set, prop[i]))
787       continue;
788     if (try_adding_simple_feature<float>(feature_set, prop[i]))
789       continue;
790     if (try_adding_simple_feature<double>(feature_set, prop[i]))
791       continue;
792   }
793 }
794 
train(int classifier,const QMultipleInputDialog & dialog)795 void Cluster_classification::train(int classifier, const QMultipleInputDialog& dialog)
796 {
797   if (m_features.size() == 0)
798   {
799     std::cerr << "Error: features not computed" << std::endl;
800     return;
801   }
802   reset_indices();
803 
804   m_label_probabilities.clear();
805   m_label_probabilities.resize (m_labels.size());
806   for (std::size_t i = 0; i < m_label_probabilities.size(); ++ i)
807     m_label_probabilities[i].resize (m_clusters.size(), -1);
808 
809   std::vector<std::size_t> nb_label (m_labels.size(), 0);
810   std::size_t nb_total = 0;
811 
812   std::vector<int> training;
813   training.reserve (m_clusters.size());
814   for (std::size_t i = 0; i < m_clusters.size(); ++ i)
815   {
816     training.push_back (m_clusters[i].training());
817     if (training.back() != -1)
818     {
819       nb_label[std::size_t(training.back())] ++;
820       ++ nb_total;
821     }
822   }
823 
824   std::cerr << nb_total << " cluster(s) used for training ("
825             << 100. * (nb_total / double(m_clusters.size())) << "% of the total):" << std::endl;
826   for (std::size_t i = 0; i < m_labels.size(); ++ i)
827     std::cerr << " * " << m_labels[i]->name() << ": " << nb_label[i] << " clusters(s)" << std::endl;
828 
829   std::vector<int> indices (m_clusters.size(), -1);
830 
831   if (classifier == CGAL_CLASSIFICATION_SOWF_NUMBER)
832   {
833     m_sowf->train<Concurrency_tag>(training, dialog.get<QSpinBox>("trials")->value());
834     CGAL::Classification::classify<Concurrency_tag> (m_clusters,
835                                                      m_labels, *m_sowf,
836                                                      indices, m_label_probabilities);
837   }
838   else if (classifier == CGAL_CLASSIFICATION_ETHZ_NUMBER)
839   {
840     if (m_ethz != nullptr)
841       delete m_ethz;
842     m_ethz = new ETHZ_random_forest (m_labels, m_features);
843     m_ethz->train<Concurrency_tag>(training, true,
844                                    dialog.get<QSpinBox>("num_trees")->value(),
845                                    dialog.get<QSpinBox>("max_depth")->value());
846     CGAL::Classification::classify<Concurrency_tag> (m_clusters,
847                                                      m_labels, *m_ethz,
848                                                      indices, m_label_probabilities);
849   }
850   else if (classifier == CGAL_CLASSIFICATION_OPENCV_NUMBER)
851   {
852 #ifdef CGAL_LINKED_WITH_OPENCV
853     if (m_random_forest != nullptr)
854       delete m_random_forest;
855     m_random_forest = new Random_forest (m_labels, m_features,
856                                          dialog.get<QSpinBox>("max_depth")->value(), 5, 15,
857                                          dialog.get<QSpinBox>("num_trees")->value());
858     m_random_forest->train (training);
859     CGAL::Classification::classify<Concurrency_tag> (m_clusters,
860                                                      m_labels, *m_random_forest,
861                                                      indices, m_label_probabilities);
862 #endif
863   }
864 
865   for (std::size_t i = 0; i < m_clusters.size(); ++ i)
866     m_clusters[i].label() = indices[i];
867 
868   if (m_index_color == 1 || m_index_color == 2)
869     change_color (m_index_color);
870 }
871 
run(int method,int classifier,std::size_t subdivisions,double smoothing)872 bool Cluster_classification::run (int method, int classifier,
873                                   std::size_t subdivisions,
874                                   double smoothing)
875 {
876   if (m_features.size() == 0)
877   {
878     std::cerr << "Error: features not computed" << std::endl;
879     return false;
880   }
881   reset_indices();
882 
883   if (classifier == CGAL_CLASSIFICATION_SOWF_NUMBER)
884     run (method, *m_sowf, subdivisions, smoothing);
885   else if (classifier == CGAL_CLASSIFICATION_ETHZ_NUMBER)
886   {
887     if (m_ethz == nullptr)
888     {
889       std::cerr << "Error: ETHZ Random Forest must be trained or have a configuration loaded first" << std::endl;
890       return false;
891     }
892     run (method, *m_ethz, subdivisions, smoothing);
893   }
894   else if (classifier == CGAL_CLASSIFICATION_OPENCV_NUMBER)
895   {
896 #ifdef CGAL_LINKED_WITH_OPENCV
897     if (m_random_forest == nullptr)
898     {
899       std::cerr << "Error: OpenCV Random Forest must be trained or have a configuration loaded first" << std::endl;
900       return false;
901     }
902     run (method, *m_random_forest, subdivisions, smoothing);
903 #endif
904   }
905 
906   return true;
907 }
908