1 // Copyright (c) 2002,2011 Utrecht University (The Netherlands). 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/Spatial_searching/include/CGAL/Weighted_Minkowski_distance.h $ 7 // $Id: Weighted_Minkowski_distance.h 0779373 2020-03-26T13:31:46+01:00 Sébastien Loriot 8 // SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial 9 // 10 // 11 // Author(s) : Hans Tangelder (<hanst@cs.uu.nl>) 12 13 // Note: Use p=0 to denote the weighted Linf-distance 14 // For 0<p<1 Lp is not a metric 15 16 #ifndef CGAL_WEIGHTED_MINKOWSKI_DISTANCE_H 17 #define CGAL_WEIGHTED_MINKOWSKI_DISTANCE_H 18 19 #include <CGAL/license/Spatial_searching.h> 20 21 22 #include <cmath> 23 #include <vector> 24 25 #include <CGAL/array.h> 26 #include <CGAL/number_utils.h> 27 #include <CGAL/Kd_tree_rectangle.h> 28 #include <CGAL/internal/Get_dimension_tag.h> 29 30 namespace CGAL { 31 namespace internal { 32 template<class T, class Dim> 33 struct Array_or_vector_selector { 34 typedef std::vector<T> type; resizeArray_or_vector_selector35 static void resize(type&v, std::size_t d) { v.resize(d); } 36 }; 37 template<class T, int D> 38 struct Array_or_vector_selector<T, Dimension_tag<D> > { 39 typedef std::array<T,D> type; 40 static void resize(type&, std::size_t CGAL_assertion_code(d)) { CGAL_assertion(d==D); } 41 }; 42 } 43 44 template <class SearchTraits> 45 class Weighted_Minkowski_distance { 46 SearchTraits traits; 47 public: 48 49 typedef typename SearchTraits::Point_d Point_d; 50 typedef Point_d Query_item; 51 typedef typename SearchTraits::FT FT; 52 typedef typename internal::Get_dimension_tag<SearchTraits>::Dimension Dimension; 53 typedef internal::Array_or_vector_selector<FT,Dimension> Weight_vector_traits; 54 typedef typename Weight_vector_traits::type Weight_vector; 55 56 private: 57 58 typedef typename SearchTraits::Cartesian_const_iterator_d Coord_iterator; 59 FT power; 60 61 Weight_vector the_weights; 62 63 public: 64 65 66 // default constructor 67 Weighted_Minkowski_distance(const SearchTraits& traits_=SearchTraits()) 68 : traits(traits_),power(2) 69 {} 70 71 Weighted_Minkowski_distance(const int d,const SearchTraits& traits_=SearchTraits()) 72 : traits(traits_),power(2), the_weights(d) 73 { 74 for (int i = 0; i < d; ++i) the_weights[i]=FT(1); 75 } 76 77 //default copy constructor and destructor 78 79 80 template <class InputIterator> 81 Weighted_Minkowski_distance (FT pow, int dim, 82 InputIterator begin, 83 InputIterator CGAL_assertion_code(end), 84 const SearchTraits& traits_=SearchTraits()) 85 : traits(traits_),power(pow) 86 { 87 CGAL_assertion(power >= FT(0)); 88 Weight_vector_traits::resize(the_weights, dim); 89 for (int i = 0; i < dim; ++i){ 90 the_weights[i] = *begin; 91 ++begin; 92 CGAL_assertion(the_weights[i]>=FT(0)); 93 } 94 CGAL_assertion(begin == end); 95 } 96 97 98 inline FT transformed_distance(const Query_item& q, const Point_d& p) const { 99 return transformed_distance(q,p, Dimension()); 100 } 101 102 //Dynamic version for runtime dimension 103 inline 104 FT 105 transformed_distance(const Query_item& q, const Point_d& p, Dynamic_dimension_tag) const 106 { 107 FT distance = FT(0); 108 typename SearchTraits::Construct_cartesian_const_iterator_d construct_it= 109 traits.construct_cartesian_const_iterator_d_object(); 110 Coord_iterator qit = construct_it(q), 111 qe = construct_it(q,1), 112 pit = construct_it(p); 113 if (power == FT(0)) { 114 for (unsigned int i = 0; qit != qe; ++qit, ++i) 115 if (the_weights[i] * CGAL::abs((*qit) - (*pit)) > distance) 116 distance = the_weights[i] * CGAL::abs((*qit)-(*pit)); 117 } 118 else 119 for (unsigned int i = 0; qit != qe; ++qit, ++i) 120 distance += 121 the_weights[i] * std::pow(CGAL::abs((*qit)-(*pit)),power); 122 return distance; 123 } 124 125 //Generic version for DIM > 3 126 template <int DIM> 127 inline FT 128 transformed_distance(const Query_item& q, const Point_d& p, Dimension_tag<DIM>) const 129 { 130 FT distance = FT(0); 131 typename SearchTraits::Construct_cartesian_const_iterator_d construct_it= 132 traits.construct_cartesian_const_iterator_d_object(); 133 Coord_iterator qit = construct_it(q), 134 qe = construct_it(q,1), 135 pit = construct_it(p); 136 if (power == FT(0)) { 137 for (unsigned int i = 0; qit != qe; ++qit, ++i) 138 if (the_weights[i] * CGAL::abs((*qit) - (*pit)) > distance) 139 distance = the_weights[i] * CGAL::abs((*qit)-(*pit)); 140 } 141 else 142 for (unsigned int i = 0; qit != qe; ++qit, ++i) 143 distance += 144 the_weights[i] * std::pow(CGAL::abs((*qit)-(*pit)),power); 145 return distance; 146 } 147 148 //DIM = 2 loop unrolled 149 inline FT 150 transformed_distance(const Query_item& q, const Point_d& p, Dimension_tag<2>) const 151 { 152 FT distance = FT(0); 153 typename SearchTraits::Construct_cartesian_const_iterator_d construct_it= 154 traits.construct_cartesian_const_iterator_d_object(); 155 Coord_iterator qit = construct_it(q), 156 pit = construct_it(p); 157 if (power == FT(0)) { 158 if (the_weights[0] * CGAL::abs((*qit) - (*pit)) > distance) 159 distance = the_weights[0] * CGAL::abs((*qit)-(*pit)); 160 qit++;pit++; 161 if (the_weights[1] * CGAL::abs((*qit) - (*pit)) > distance) 162 distance = the_weights[1] * CGAL::abs((*qit)-(*pit)); 163 } 164 else{ 165 distance += 166 the_weights[0] * std::pow(CGAL::abs((*qit)-(*pit)),power); 167 qit++;pit++; 168 distance += 169 the_weights[1] * std::pow(CGAL::abs((*qit)-(*pit)),power); 170 } 171 return distance; 172 } 173 174 //DIM = 3 loop unrolled 175 inline FT 176 transformed_distance(const Query_item& q, const Point_d& p, Dimension_tag<3>) const 177 { 178 FT distance = FT(0); 179 typename SearchTraits::Construct_cartesian_const_iterator_d construct_it= 180 traits.construct_cartesian_const_iterator_d_object(); 181 Coord_iterator qit = construct_it(q), 182 pit = construct_it(p); 183 if (power == FT(0)) { 184 if (the_weights[0] * CGAL::abs((*qit) - (*pit)) > distance) 185 distance = the_weights[0] * CGAL::abs((*qit)-(*pit)); 186 qit++;pit++; 187 if (the_weights[1] * CGAL::abs((*qit) - (*pit)) > distance) 188 distance = the_weights[1] * CGAL::abs((*qit)-(*pit)); 189 qit++;pit++; 190 if (the_weights[2] * CGAL::abs((*qit) - (*pit)) > distance) 191 distance = the_weights[2] * CGAL::abs((*qit)-(*pit)); 192 } 193 else{ 194 distance += 195 the_weights[0] * std::pow(CGAL::abs((*qit)-(*pit)),power); 196 qit++;pit++; 197 distance += 198 the_weights[1] * std::pow(CGAL::abs((*qit)-(*pit)),power); 199 qit++;pit++; 200 distance += 201 the_weights[2] * std::pow(CGAL::abs((*qit)-(*pit)),power); 202 } 203 return distance; 204 } 205 206 inline 207 FT 208 min_distance_to_rectangle(const Query_item& q, 209 const Kd_tree_rectangle<FT,Dimension>& r) const 210 { 211 FT distance = FT(0); 212 typename SearchTraits::Construct_cartesian_const_iterator_d construct_it= 213 traits.construct_cartesian_const_iterator_d_object(); 214 Coord_iterator qit = construct_it(q), qe = construct_it(q,1); 215 if (power == FT(0)) 216 { 217 for (unsigned int i = 0; qit != qe; ++qit, ++i) { 218 if (the_weights[i]*(r.min_coord(i) - 219 (*qit)) > distance) 220 distance = the_weights[i] * (r.min_coord(i)- 221 (*qit)); 222 if (the_weights[i] * ((*qit) - r.max_coord(i)) > 223 distance) 224 distance = the_weights[i] * 225 ((*qit)-r.max_coord(i)); 226 } 227 } 228 else 229 { 230 for (unsigned int i = 0; qit != qe; ++qit, ++i) { 231 if ((*qit) < r.min_coord(i)) 232 distance += the_weights[i] * 233 std::pow(r.min_coord(i)-(*qit),power); 234 if ((*qit) > r.max_coord(i)) 235 distance += the_weights[i] * 236 std::pow((*qit)-r.max_coord(i),power); 237 } 238 }; 239 return distance; 240 } 241 242 inline 243 FT 244 min_distance_to_rectangle(const Query_item& q, 245 const Kd_tree_rectangle<FT,Dimension>& r,std::vector<FT>& dists) const { 246 FT distance = FT(0); 247 typename SearchTraits::Construct_cartesian_const_iterator_d construct_it= 248 traits.construct_cartesian_const_iterator_d_object(); 249 Coord_iterator qit = construct_it(q), qe = construct_it(q,1); 250 if (power == FT(0)) 251 { 252 for (unsigned int i = 0; qit != qe; ++qit, ++i) { 253 if (the_weights[i]*(r.min_coord(i) - 254 (*qit)) > distance){ 255 dists[i] = (r.min_coord(i)- 256 (*qit)); 257 distance = the_weights[i] * dists[i]; 258 } 259 if (the_weights[i] * ((*qit) - r.max_coord(i)) > 260 distance){ 261 dists[i] = 262 ((*qit)-r.max_coord(i)); 263 distance = the_weights[i] * dists[i]; 264 } 265 } 266 } 267 else 268 { 269 for (unsigned int i = 0; qit != qe; ++qit, ++i) { 270 if ((*qit) < r.min_coord(i)){ 271 dists[i] = r.min_coord(i)-(*qit); 272 distance += the_weights[i] * 273 std::pow(dists[i],power); 274 } 275 if ((*qit) > r.max_coord(i)){ 276 dists[i] = (*qit)-r.max_coord(i); 277 distance += the_weights[i] * 278 std::pow(dists[i],power); 279 } 280 } 281 }; 282 return distance; 283 } 284 285 inline 286 FT 287 max_distance_to_rectangle(const Query_item& q, 288 const Kd_tree_rectangle<FT,Dimension>& r) const { 289 FT distance=FT(0); 290 typename SearchTraits::Construct_cartesian_const_iterator_d construct_it= 291 traits.construct_cartesian_const_iterator_d_object(); 292 Coord_iterator qit = construct_it(q), qe = construct_it(q,1); 293 if (power == FT(0)) 294 { 295 for (unsigned int i = 0; qit != qe; ++qit, ++i) { 296 if ((*qit) >= (r.min_coord(i) + 297 r.max_coord(i))/FT(2.0)) { 298 if (the_weights[i] * ((*qit) - 299 r.min_coord(i)) > distance) 300 distance = the_weights[i] * 301 ((*qit)-r.min_coord(i)); 302 else 303 if (the_weights[i] * 304 (r.max_coord(i) - (*qit)) > distance) 305 distance = the_weights[i] * 306 ( r.max_coord(i)-(*qit)); 307 } 308 } 309 } 310 else 311 { 312 for (unsigned int i = 0; qit != qe; ++qit, ++i) { 313 if ((*qit) <= (r.min_coord(i)+r.max_coord(i))/FT(2.0)) 314 distance += the_weights[i] * std::pow(r.max_coord(i)-(*qit),power); 315 else 316 distance += the_weights[i] * std::pow((*qit)-r.min_coord(i),power); 317 } 318 }; 319 return distance; 320 } 321 322 inline 323 FT 324 max_distance_to_rectangle(const Query_item& q, 325 const Kd_tree_rectangle<FT,Dimension>& r,std::vector<FT>& dists) const { 326 FT distance=FT(0); 327 typename SearchTraits::Construct_cartesian_const_iterator_d construct_it= 328 traits.construct_cartesian_const_iterator_d_object(); 329 Coord_iterator qit = construct_it(q), qe = construct_it(q,1); 330 if (power == FT(0)) 331 { 332 for (unsigned int i = 0; qit != qe; ++qit, ++i) { 333 if ((*qit) >= (r.min_coord(i) + 334 r.max_coord(i))/FT(2.0)) { 335 if (the_weights[i] * ((*qit) - 336 r.min_coord(i)) > distance){ 337 dists[i] = (*qit)-r.min_coord(i); 338 distance = the_weights[i] * 339 (dists[i]); 340 } 341 else 342 if (the_weights[i] * 343 (r.max_coord(i) - (*qit)) > distance){ 344 dists[i] = r.max_coord(i)-(*qit); 345 distance = the_weights[i] * 346 (dists[i]); 347 } 348 } 349 } 350 } 351 else 352 { 353 for (unsigned int i = 0; qit != qe; ++qit, ++i) { 354 if ((*qit) <= (r.min_coord(i)+r.max_coord(i))/FT(2.0)){ 355 dists[i] = r.max_coord(i)-(*qit); 356 distance += the_weights[i] * std::pow(dists[i],power); 357 } 358 else{ 359 dists[i] = (*qit)-r.min_coord(i); 360 distance += the_weights[i] * std::pow(dists[i],power); 361 } 362 } 363 }; 364 return distance; 365 } 366 367 inline 368 FT 369 new_distance(FT dist, FT old_off, FT new_off, 370 int cutting_dimension) const 371 { 372 FT new_dist; 373 if (power == FT(0)) 374 { 375 if (the_weights[cutting_dimension]*CGAL::abs(new_off) 376 > dist) 377 new_dist= 378 the_weights[cutting_dimension]*CGAL::abs(new_off); 379 else new_dist=dist; 380 } 381 else 382 { 383 new_dist = dist + the_weights[cutting_dimension] * 384 (std::pow(CGAL::abs(new_off),power)-std::pow(CGAL::abs(old_off),power)); 385 } 386 return new_dist; 387 } 388 389 inline 390 FT 391 transformed_distance(FT d) const 392 { 393 if (power <= FT(0)) return d; 394 else return std::pow(d,power); 395 396 } 397 398 inline 399 FT 400 inverse_of_transformed_distance(FT d) const 401 { 402 if (power <= FT(0)) return d; 403 else return std::pow(d,1/power); 404 405 } 406 407 }; // class Weighted_Minkowski_distance 408 409 } // namespace CGAL 410 411 #endif // CGAL_WEIGHTED_MINKOWSKI_DISTANCE_H 412