diff --git a/Spatial_searching/doc/Spatial_searching/CGAL/Weighted_Minkowski_distance.h b/Spatial_searching/doc/Spatial_searching/CGAL/Weighted_Minkowski_distance.h index e1575fb05f3..80dd8677699 100644 --- a/Spatial_searching/doc/Spatial_searching/CGAL/Weighted_Minkowski_distance.h +++ b/Spatial_searching/doc/Spatial_searching/CGAL/Weighted_Minkowski_distance.h @@ -10,6 +10,9 @@ defined by \f$ l_{\infty}(w)(r,q)=max \{w_i |r_i-q_i| \mid 1 \leq i \leq d\}\f$. For the purpose of the distance computations it is more efficient to compute the transformed distance \f$ {\sigma_{i=1}^{i=d} \, w_i(r_i-q_i)^p}\f$ instead of the actual distance. +\note As this distance involves the computation of a power it is not +done exact but with floating point arithmetic. + \tparam Traits must be a model of the concept `SearchTraits`, for example `Search_traits_2`. diff --git a/Spatial_searching/include/CGAL/Weighted_Minkowski_distance.h b/Spatial_searching/include/CGAL/Weighted_Minkowski_distance.h index f38f5b8cf37..fcc19d99841 100644 --- a/Spatial_searching/include/CGAL/Weighted_Minkowski_distance.h +++ b/Spatial_searching/include/CGAL/Weighted_Minkowski_distance.h @@ -18,15 +18,16 @@ #include +#include +#include +#include +#include +#include +#include #include #include -#include -#include -#include -#include - namespace CGAL { namespace internal { template @@ -56,10 +57,22 @@ namespace CGAL { private: typedef typename SearchTraits::Cartesian_const_iterator_d Coord_iterator; - FT power; + double power; + Weight_vector the_weights; + + FT pow(FT ft) const + { + return FT(std::pow(to_double(ft), power)); + } + + FT pow(FT ft, FT p) const + { + return FT(std::pow(to_double(ft), to_double(p))); + } + public: @@ -82,9 +95,9 @@ namespace CGAL { InputIterator begin, InputIterator CGAL_assertion_code(end), const SearchTraits& traits_=SearchTraits()) - : traits(traits_),power(pow) + : traits(traits_),power(to_double(pow)) { - CGAL_assertion(power >= FT(0)); + CGAL_assertion(is_positive(power)); Weight_vector_traits::resize(the_weights, dim); for (int i = 0; i < dim; ++i){ the_weights[i] = *begin; @@ -110,7 +123,7 @@ namespace CGAL { Coord_iterator qit = construct_it(q), qe = construct_it(q,1), pit = construct_it(p); - if (power == FT(0)) { + if (is_zero(power)) { for (unsigned int i = 0; qit != qe; ++qit, ++i) if (the_weights[i] * CGAL::abs((*qit) - (*pit)) > distance) distance = the_weights[i] * CGAL::abs((*qit)-(*pit)); @@ -118,7 +131,7 @@ namespace CGAL { else for (unsigned int i = 0; qit != qe; ++qit, ++i) distance += - the_weights[i] * std::pow(CGAL::abs((*qit)-(*pit)),power); + the_weights[i] * pow(CGAL::abs((*qit)-(*pit))); return distance; } @@ -133,7 +146,7 @@ namespace CGAL { Coord_iterator qit = construct_it(q), qe = construct_it(q,1), pit = construct_it(p); - if (power == FT(0)) { + if (is_zero(power)) { for (unsigned int i = 0; qit != qe; ++qit, ++i) if (the_weights[i] * CGAL::abs((*qit) - (*pit)) > distance) distance = the_weights[i] * CGAL::abs((*qit)-(*pit)); @@ -141,7 +154,7 @@ namespace CGAL { else for (unsigned int i = 0; qit != qe; ++qit, ++i) distance += - the_weights[i] * std::pow(CGAL::abs((*qit)-(*pit)),power); + the_weights[i] * pow(CGAL::abs((*qit)-(*pit))); return distance; } @@ -154,7 +167,7 @@ namespace CGAL { traits.construct_cartesian_const_iterator_d_object(); Coord_iterator qit = construct_it(q), pit = construct_it(p); - if (power == FT(0)) { + if (is_zero(power)) { if (the_weights[0] * CGAL::abs((*qit) - (*pit)) > distance) distance = the_weights[0] * CGAL::abs((*qit)-(*pit)); qit++;pit++; @@ -163,10 +176,10 @@ namespace CGAL { } else{ distance += - the_weights[0] * std::pow(CGAL::abs((*qit)-(*pit)),power); + the_weights[0] * pow(CGAL::abs((*qit)-(*pit))); qit++;pit++; distance += - the_weights[1] * std::pow(CGAL::abs((*qit)-(*pit)),power); + the_weights[1] * pow(CGAL::abs((*qit)-(*pit))); } return distance; } @@ -180,7 +193,7 @@ namespace CGAL { traits.construct_cartesian_const_iterator_d_object(); Coord_iterator qit = construct_it(q), pit = construct_it(p); - if (power == FT(0)) { + if (is_zero(power)) { if (the_weights[0] * CGAL::abs((*qit) - (*pit)) > distance) distance = the_weights[0] * CGAL::abs((*qit)-(*pit)); qit++;pit++; @@ -192,13 +205,13 @@ namespace CGAL { } else{ distance += - the_weights[0] * std::pow(CGAL::abs((*qit)-(*pit)),power); + the_weights[0] * pow(CGAL::abs((*qit)-(*pit))); qit++;pit++; distance += - the_weights[1] * std::pow(CGAL::abs((*qit)-(*pit)),power); + the_weights[1] * pow(CGAL::abs((*qit)-(*pit))); qit++;pit++; distance += - the_weights[2] * std::pow(CGAL::abs((*qit)-(*pit)),power); + the_weights[2] * pow(CGAL::abs((*qit)-(*pit))); } return distance; } @@ -212,7 +225,7 @@ namespace CGAL { typename SearchTraits::Construct_cartesian_const_iterator_d construct_it= traits.construct_cartesian_const_iterator_d_object(); Coord_iterator qit = construct_it(q), qe = construct_it(q,1); - if (power == FT(0)) + if (is_zero(power)) { for (unsigned int i = 0; qit != qe; ++qit, ++i) { if (the_weights[i]*(r.min_coord(i) - @@ -230,10 +243,10 @@ namespace CGAL { for (unsigned int i = 0; qit != qe; ++qit, ++i) { if ((*qit) < r.min_coord(i)) distance += the_weights[i] * - std::pow(r.min_coord(i)-(*qit),power); + pow(r.min_coord(i)-(*qit)); if ((*qit) > r.max_coord(i)) distance += the_weights[i] * - std::pow((*qit)-r.max_coord(i),power); + pow((*qit)-r.max_coord(i)); } }; return distance; @@ -247,7 +260,7 @@ namespace CGAL { typename SearchTraits::Construct_cartesian_const_iterator_d construct_it= traits.construct_cartesian_const_iterator_d_object(); Coord_iterator qit = construct_it(q), qe = construct_it(q,1); - if (power == FT(0)) + if (is_zero(power)) { for (unsigned int i = 0; qit != qe; ++qit, ++i) { if (the_weights[i]*(r.min_coord(i) - @@ -270,12 +283,12 @@ namespace CGAL { if ((*qit) < r.min_coord(i)){ dists[i] = r.min_coord(i)-(*qit); distance += the_weights[i] * - std::pow(dists[i],power); + pow(dists[i]); } if ((*qit) > r.max_coord(i)){ dists[i] = (*qit)-r.max_coord(i); distance += the_weights[i] * - std::pow(dists[i],power); + pow(dists[i]); } } }; @@ -290,7 +303,7 @@ namespace CGAL { typename SearchTraits::Construct_cartesian_const_iterator_d construct_it= traits.construct_cartesian_const_iterator_d_object(); Coord_iterator qit = construct_it(q), qe = construct_it(q,1); - if (power == FT(0)) + if (is_zero(power)) { for (unsigned int i = 0; qit != qe; ++qit, ++i) { if ((*qit) >= (r.min_coord(i) + @@ -311,9 +324,9 @@ namespace CGAL { { for (unsigned int i = 0; qit != qe; ++qit, ++i) { if ((*qit) <= (r.min_coord(i)+r.max_coord(i))/FT(2.0)) - distance += the_weights[i] * std::pow(r.max_coord(i)-(*qit),power); + distance += the_weights[i] * pow(r.max_coord(i)-(*qit)); else - distance += the_weights[i] * std::pow((*qit)-r.min_coord(i),power); + distance += the_weights[i] * pow((*qit)-r.min_coord(i)); } }; return distance; @@ -327,7 +340,7 @@ namespace CGAL { typename SearchTraits::Construct_cartesian_const_iterator_d construct_it= traits.construct_cartesian_const_iterator_d_object(); Coord_iterator qit = construct_it(q), qe = construct_it(q,1); - if (power == FT(0)) + if (is_zero(power)) { for (unsigned int i = 0; qit != qe; ++qit, ++i) { if ((*qit) >= (r.min_coord(i) + @@ -353,11 +366,11 @@ namespace CGAL { for (unsigned int i = 0; qit != qe; ++qit, ++i) { if ((*qit) <= (r.min_coord(i)+r.max_coord(i))/FT(2.0)){ dists[i] = r.max_coord(i)-(*qit); - distance += the_weights[i] * std::pow(dists[i],power); + distance += the_weights[i] * pow(dists[i]); } else{ dists[i] = (*qit)-r.min_coord(i); - distance += the_weights[i] * std::pow(dists[i],power); + distance += the_weights[i] * pow(dists[i]); } } }; @@ -370,7 +383,7 @@ namespace CGAL { int cutting_dimension) const { FT new_dist; - if (power == FT(0)) + if (is_zero(power)) { if (the_weights[cutting_dimension]*CGAL::abs(new_off) > dist) @@ -381,7 +394,7 @@ namespace CGAL { else { new_dist = dist + the_weights[cutting_dimension] * - (std::pow(CGAL::abs(new_off),power)-std::pow(CGAL::abs(old_off),power)); + (pow(CGAL::abs(new_off))-pow(CGAL::abs(old_off))); } return new_dist; } @@ -390,8 +403,8 @@ namespace CGAL { FT transformed_distance(FT d) const { - if (power <= FT(0)) return d; - else return std::pow(d,power); + if (! is_positive(power)) return d; + else return pow(d); } @@ -399,8 +412,8 @@ namespace CGAL { FT inverse_of_transformed_distance(FT d) const { - if (power <= FT(0)) return d; - else return std::pow(d,1/power); + if (! is_positive(power)) return d; + else return pow(d,1/power); }