1 // Copyright (c) 2017 GeometryFactory Sarl (France). 2 // All rights reserved. 3 // 4 // This file is part of CGAL (www.cgal.org). 5 // 6 // $URL: https://github.com/CGAL/cgal/blob/v5.3/Classification/include/CGAL/Classification/Feature/Vertical_dispersion.h $ 7 // $Id: Vertical_dispersion.h 0e934b1 2020-08-04T13:16:13+02:00 Simon Giraudot 8 // SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial 9 // 10 // Author(s) : Simon Giraudot 11 12 #ifndef CGAL_CLASSIFICATION_FEATURE_VERTICAL_DISPERSION_H 13 #define CGAL_CLASSIFICATION_FEATURE_VERTICAL_DISPERSION_H 14 15 #include <CGAL/license/Classification.h> 16 17 #include <vector> 18 19 #include <CGAL/number_utils.h> 20 #include <CGAL/Classification/Image.h> 21 #include <CGAL/Classification/compressed_float.h> 22 #include <CGAL/Classification/Planimetric_grid.h> 23 #include <boost/algorithm/minmax_element.hpp> 24 #include <CGAL/Classification/Feature_base.h> 25 #include <CGAL/int.h> 26 #include <CGAL/float.h> 27 #include <boost/tuple/tuple.hpp> 28 29 namespace CGAL { 30 31 namespace Classification { 32 33 namespace Feature { 34 35 /*! 36 \ingroup PkgClassificationFeatures 37 38 %Feature based on local vertical dispersion of points. Urban 39 scenes can often be decomposed as a set of 2D regions with 40 different heights. While these heights are usually piecewise 41 constant or piecewise linear, on some specific parts of the scene 42 such as vegetation, they can become extremely unstable. This 43 feature quantifies the vertical dispersion of the points on a 44 local Z-cylinder around the points. 45 46 Its default name is "vertical_dispersion". 47 48 \tparam GeomTraits model of \cgal Kernel. 49 \tparam PointRange model of `ConstRange`. Its iterator type 50 is `RandomAccessIterator` and its value type is the key type of 51 `PointMap`. 52 \tparam PointMap model of `ReadablePropertyMap` whose key 53 type is the value type of the iterator of `PointRange` and value type 54 is `GeomTraits::Point_3`. 55 */ 56 template <typename GeomTraits, typename PointRange, typename PointMap> 57 class Vertical_dispersion : public Feature_base 58 { 59 using Image_cfloat = Classification::Image<compressed_float>; 60 using Grid = Classification::Planimetric_grid<GeomTraits, PointRange, PointMap>; 61 62 const Grid& grid; 63 Image_cfloat Dispersion; 64 std::vector<compressed_float> values; 65 66 public: 67 /*! 68 \brief constructs the feature. 69 70 \param input point range. 71 \param point_map property map to access the input points. 72 \param grid precomputed `Planimetric_grid`. 73 \param radius_neighbors radius of local neighborhoods. 74 */ 75 Vertical_dispersion (const PointRange& input, 76 PointMap point_map, 77 const Grid& grid, 78 float radius_neighbors = -1.) grid(grid)79 : grid (grid) 80 { 81 this->set_name ("vertical_dispersion"); 82 if (radius_neighbors < 0.) 83 radius_neighbors = 5.f * grid.resolution(); 84 85 if (grid.width() * grid.height() > input.size()) 86 values.resize(input.size(), compressed_float(0)); 87 else 88 { 89 Dispersion = Image_cfloat(grid.width(), grid.height()); 90 for (std::size_t j = 0; j < grid.height(); j++) 91 for (std::size_t i = 0; i < grid.width(); i++) 92 if (grid.has_points(i,j)) 93 Dispersion(i,j) = compressed_float(0); 94 } 95 96 std::size_t square = (std::size_t)(0.5 * radius_neighbors / grid.resolution()) + 1; 97 typename GeomTraits::Vector_3 verti (0., 0., 1.); 98 99 std::vector<float> hori; 100 101 for (std::size_t j = 0; j < grid.height(); j++){ 102 for (std::size_t i = 0; i < grid.width(); i++){ 103 104 if(!(grid.has_points(i,j))) 105 continue; 106 hori.clear(); 107 108 std::size_t squareXmin = (i < square ? 0 : i - square); 109 std::size_t squareXmax = (std::min) (grid.width()-1, i + square); 110 std::size_t squareYmin = (j < square ? 0 : j - square); 111 std::size_t squareYmax = (std::min) (grid.height()-1, j + square); 112 113 float bound = (float)0.5*radius_neighbors/grid.resolution(); 114 bound = CGAL::square(bound); 115 for(std::size_t k = squareXmin; k <= squareXmax; k++) 116 for(std::size_t l = squareYmin; l <= squareYmax; l++) 117 { 118 if(CGAL::square((float)(k-i))+ CGAL::square((float)(l-j)) 119 <= bound) 120 { 121 for (typename Grid::iterator it = grid.indices_begin(k,l); it != grid.indices_end(k,l); ++ it) 122 hori.push_back (float(get(point_map, *(input.begin()+(*it))).z())); 123 } 124 } 125 126 if (hori.empty()) 127 continue; 128 129 std::vector<float>::iterator min_it, max_it; 130 boost::tie(min_it, max_it) 131 = boost::minmax_element (hori.begin(), hori.end()); 132 133 std::vector<bool> occupy (1 + (std::size_t)((*max_it - *min_it) / grid.resolution()), false); 134 135 for (std::size_t k = 0; k < hori.size(); ++ k) 136 { 137 std::size_t index = (std::size_t)((hori[k] - *min_it) / grid.resolution()); 138 occupy[index] = true; 139 } 140 141 std::size_t nb_occ = 0; 142 for (std::size_t k = 0; k < occupy.size(); ++ k) 143 if (occupy[k]) 144 ++ nb_occ; 145 146 compressed_float v = compress_float (1.f - (nb_occ / (float)(occupy.size()))); 147 if (values.empty()) 148 Dispersion(i,j) = v; 149 else 150 { 151 typename Grid::iterator end = grid.indices_end(i,j); 152 for (typename Grid::iterator it = grid.indices_begin(i,j); it != end; ++ it) 153 values[*it] = v; 154 } 155 } 156 157 } 158 } 159 /// \cond SKIP_IN_MANUAL value(std::size_t pt_index)160 virtual float value (std::size_t pt_index) 161 { 162 if (values.empty()) 163 { 164 std::size_t I = grid.x(pt_index); 165 std::size_t J = grid.y(pt_index); 166 return decompress_float (Dispersion(I,J)); 167 } 168 return decompress_float (values[pt_index]); 169 } 170 /// \endcond 171 }; 172 173 } 174 175 } 176 177 } 178 179 #endif // CGAL_CLASSIFICATION_FEATURE_VERTICAL_DISPERSION_H 180