From af343d926efbf10ae7d931c48ae21f799a48c6c2 Mon Sep 17 00:00:00 2001 From: "m.overtheil" Date: Tue, 2 Dec 2014 10:36:36 +0100 Subject: [PATCH] Added new traversal rule and dists vector to Orthogonal search Functions min_dist_to_rectangle and max... have an overload with std::vector& dists now. That is required for ortho searches, so we have no backwards compatibility. The values for extended internal nodes have also been changed. --- .../benchmark/Spatial_searching/nn3cgal.cpp | 7 +- .../include/CGAL/Euclidean_distance.h | 39 ++++++++ .../CGAL/Euclidean_distance_sphere_point.h | 45 ++++++++++ Spatial_searching/include/CGAL/Kd_tree.h | 11 ++- .../include/CGAL/Kd_tree_rectangle.h | 6 +- .../CGAL/Manhattan_distance_iso_box_point.h | 46 ++++++++++ .../CGAL/Orthogonal_k_neighbor_search.h | 82 +++++++++-------- .../include/CGAL/Search_traits_adapter.h | 11 +++ .../CGAL/Weighted_Minkowski_distance.h | 88 +++++++++++++++++++ .../test/Spatial_searching/Distance.h | 50 +++++++++++ 10 files changed, 340 insertions(+), 45 deletions(-) diff --git a/Spatial_searching/benchmark/Spatial_searching/nn3cgal.cpp b/Spatial_searching/benchmark/Spatial_searching/nn3cgal.cpp index 07c47cd27e1..26c780a07a7 100644 --- a/Spatial_searching/benchmark/Spatial_searching/nn3cgal.cpp +++ b/Spatial_searching/benchmark/Spatial_searching/nn3cgal.cpp @@ -87,10 +87,11 @@ int main(int argc,char *argv[]) for (Neighbor_search_3::iterator it = search.begin(); it != search.end();it++, i++) { result[i] = it->first; if(dump){ - std::cerr << result[i].x()<<" "<& r,std::vector& dists) { + FT distance = FT(0); + typename SearchTraits::Construct_cartesian_const_iterator_d construct_it=traits.construct_cartesian_const_iterator_d_object(); + typename SearchTraits::Cartesian_const_iterator_d qit = construct_it(q), + qe = construct_it(q,1); + for(unsigned int i = 0;qit != qe; i++, qit++){ + if((*qit) < r.min_coord(i)){ + dists[i] = (r.min_coord(i)-(*qit)); + distance += dists[i] * dists[i]; + } + else if ((*qit) > r.max_coord(i)){ + dists[i] = ((*qit)-r.max_coord(i)); + distance += dists[i] * dists[i]; + } + + } + return distance; + } + inline FT max_distance_to_rectangle(const Query_item& q, const Kd_tree_rectangle& r) const { FT distance=FT(0); @@ -123,6 +143,25 @@ namespace CGAL { return distance; } + inline FT max_distance_to_rectangle(const Query_item& q, + const Kd_tree_rectangle& r,std::vector& dists ) { + FT distance=FT(0); + typename SearchTraits::Construct_cartesian_const_iterator_d construct_it=traits.construct_cartesian_const_iterator_d_object(); + typename SearchTraits::Cartesian_const_iterator_d qit = construct_it(q), + qe = construct_it(q,1); + for(unsigned int i = 0;qit != qe; i++, qit++){ + if ((*qit) <= (r.min_coord(i)+r.max_coord(i))/FT(2.0)){ + dists[i] = (r.max_coord(i)-(*qit)); + distance += dists[i] * dists[i]; + } + else{ + dists[i] = ((*qit)-r.min_coord(i)); + distance += dists[i] * dists[i]; + } + }; + return distance; + } + inline FT new_distance(FT dist, FT old_off, FT new_off, int /* cutting_dimension */) const { diff --git a/Spatial_searching/include/CGAL/Euclidean_distance_sphere_point.h b/Spatial_searching/include/CGAL/Euclidean_distance_sphere_point.h index 6945b5fb8ec..1b602bf07e9 100644 --- a/Spatial_searching/include/CGAL/Euclidean_distance_sphere_point.h +++ b/Spatial_searching/include/CGAL/Euclidean_distance_sphere_point.h @@ -107,6 +107,28 @@ namespace CGAL { return distance; } + inline FT min_distance_to_rectangle(const Sphere_d& q, + const Kd_tree_rectangle& r,std::vector& dists) { + Point_d c= Construct_center_d()(q); + FT distance = FT(0); + Construct_cartesian_const_iterator_d construct_it=traits.construct_cartesian_const_iterator_d_object(); + Cartesian_const_iterator_d cit = construct_it(c), + ce = construct_it(c,1); + for (unsigned int i = 0; cit != ce; ++i, ++cit) { + if ((*cit) < r.min_coord(i)){ + dists[i] =(r.min_coord(i)-(*cit)); + distance += dists[i] * dists[i]; + } + else if ((*cit) > r.max_coord(i)){ + dists[i] = ((*cit)-r.max_coord(i)); + distance += dists[i] * dists[i]; + } + }; + distance += - Compute_squared_radius_d()(q); + if (distance<0) distance=FT(0); + return distance; + } + inline FT max_distance_to_rectangle(const Sphere_d& q, const Kd_tree_rectangle& r) const { Construct_center_d construct_center_d; @@ -126,6 +148,29 @@ namespace CGAL { return distance; } + inline FT max_distance_to_rectangle(const Sphere_d& q, + const Kd_tree_rectangle& r,std::vector& dists) { + Construct_center_d construct_center_d; + Point_d c = construct_center_d(q); + FT distance=FT(0); + Construct_cartesian_const_iterator_d construct_it=traits.construct_cartesian_const_iterator_d_object(); + Cartesian_const_iterator_d cit = construct_it(c), + ce = construct_it(c,1); + for (unsigned int i = 0; cit != ce; ++i, ++cit) { + if ((*cit) <= (r.min_coord(i)+r.max_coord(i))/FT(2.0)){ + dists[i] = (r.max_coord(i)-(*cit)); + distance += dists[i] * dists[i]; + } + else{ + dists[i] = ((*cit)-r.min_coord(i)); + distance += dists[i] * dists[i]; + } + }; + distance += - Compute_squared_radius_d()(q); + if (distance<0) distance=FT(0); + return distance; + } + inline FT transformed_distance(FT d) const { diff --git a/Spatial_searching/include/CGAL/Kd_tree.h b/Spatial_searching/include/CGAL/Kd_tree.h index 7d328e2987c..f08c1aaa052 100644 --- a/Spatial_searching/include/CGAL/Kd_tree.h +++ b/Spatial_searching/include/CGAL/Kd_tree.h @@ -175,9 +175,14 @@ private: split(node.separator(), c, c_low); int cd = node.separator().cutting_dimension(); - - node.low_val = c_low.bounding_box().min_coord(cd); - node.high_val = c.bounding_box().max_coord(cd); + if(!c_low.empty()) + node.low_val = c_low.tight_bounding_box().max_coord(cd); + else + node.low_val = c_low.bounding_box().min_coord(cd); + if(!c.empty()) + node.high_val = c.tight_bounding_box().min_coord(cd); + else + node.high_val = c.bounding_box().max_coord(cd); CGAL_assertion(node.separator().cutting_value() >= node.low_val); CGAL_assertion(node.separator().cutting_value() <= node.high_val); diff --git a/Spatial_searching/include/CGAL/Kd_tree_rectangle.h b/Spatial_searching/include/CGAL/Kd_tree_rectangle.h index 8d8a6f5e4d1..8b108621ef1 100644 --- a/Spatial_searching/include/CGAL/Kd_tree_rectangle.h +++ b/Spatial_searching/include/CGAL/Kd_tree_rectangle.h @@ -175,14 +175,14 @@ namespace CGAL { std::ostream& print(std::ostream& s) const { - s << "Rectangle dimension = " << D; + s << "Rectangle dimension = " << D::value; s << "\n lower: "; - for (int i=0; i < D; ++i) + for (int i=0; i < D::value; ++i) s << lower_[i] << " "; // std::copy(lower_, lower_ + D, // std::ostream_iterator(s," ")); s << "\n upper: "; - for (int j=0; j < D; ++j) + for (int j=0; j < D::value; ++j) s << upper_[j] << " "; // std::copy(upper_, upper_ + D, // std::ostream_iterator(s," ")); diff --git a/Spatial_searching/include/CGAL/Manhattan_distance_iso_box_point.h b/Spatial_searching/include/CGAL/Manhattan_distance_iso_box_point.h index 07f6ed292db..263cb338378 100644 --- a/Spatial_searching/include/CGAL/Manhattan_distance_iso_box_point.h +++ b/Spatial_searching/include/CGAL/Manhattan_distance_iso_box_point.h @@ -101,6 +101,28 @@ namespace CGAL { return distance; } + inline FT min_distance_to_rectangle(const Query_item& q, + const Kd_tree_rectangle& r,std::vector& dists) { + FT distance = FT(0); + typename SearchTraits::Construct_cartesian_const_iterator_d construct_it= + traits.construct_cartesian_const_iterator_d_object(); + typename SearchTraits::Construct_min_vertex_d construct_min_vertex; + typename SearchTraits::Construct_max_vertex_d construct_max_vertex; + typename SearchTraits::Cartesian_const_iterator_d qmaxit = construct_it(construct_max_vertex(q)), + qe = construct_it(construct_max_vertex(q),1), qminit = construct_it(construct_min_vertex(q)); + for (unsigned int i = 0; qmaxit != qe; ++ qmaxit, ++qminit, ++i) { + if (r.min_coord(i)>(*qmaxit)) { + dists[i]=(r.min_coord(i)-(*qmaxit)); + distance +=dists[i]; + } + if (r.max_coord(i)<(*qminit)) { + dists[i]=((*qminit)-r.max_coord(i)); + distance += dists[i]; + } + } + return distance; + } + inline FT max_distance_to_rectangle(const Query_item& q, @@ -120,6 +142,30 @@ namespace CGAL { } return distance; } + + inline + FT + max_distance_to_rectangle(const Query_item& q, + const Kd_tree_rectangle& r,std::vector& dists) { + FT distance=FT(0); + typename SearchTraits::Construct_cartesian_const_iterator_d construct_it= + traits.construct_cartesian_const_iterator_d_object(); + typename SearchTraits::Construct_min_vertex_d construct_min_vertex; + typename SearchTraits::Construct_max_vertex_d construct_max_vertex; + typename SearchTraits::Cartesian_const_iterator_d qmaxit = construct_it(construct_max_vertex(q)), + qe = construct_it(construct_max_vertex(q),1), qminit = construct_it(construct_min_vertex(q)); + for (unsigned int i = 0; qmaxit != qe; ++ qmaxit, ++qminit, ++i) { + if ( r.max_coord(i)-(*qminit) >(*qmaxit)-r.min_coord(i) ) { + dists[i]=(r.max_coord(i)-(*qminit)); + distance += dists[i]; + } + else { + dists[i]=((*qmaxit)-r.min_coord(i)); + distance += dists[i]; + } + } + return distance; + } inline FT diff --git a/Spatial_searching/include/CGAL/Orthogonal_k_neighbor_search.h b/Spatial_searching/include/CGAL/Orthogonal_k_neighbor_search.h index 94ad2d72a03..21e4bc4bd6d 100644 --- a/Spatial_searching/include/CGAL/Orthogonal_k_neighbor_search.h +++ b/Spatial_searching/include/CGAL/Orthogonal_k_neighbor_search.h @@ -43,23 +43,30 @@ public: { if (tree.empty()) return; - FT distance_to_root; - if (this->search_nearest) - distance_to_root = this->distance_instance.min_distance_to_rectangle(q, tree.bounding_box()); - else - distance_to_root = this->distance_instance.max_distance_to_rectangle(q, tree.bounding_box()); - - typename SearchTraits::Construct_cartesian_const_iterator_d construct_it=tree.traits().construct_cartesian_const_iterator_d_object(); query_object_it = construct_it(this->query_object); + + int dim = static_cast(std::distance(query_object_it, construct_it(this->query_object,0))); + std::vector dists(dim); + for(int i=0;isearch_nearest) + distance_to_root = this->distance_instance.min_distance_to_rectangle(q, tree.bounding_box(),dists); + else + distance_to_root = this->distance_instance.max_distance_to_rectangle(q, tree.bounding_box(),dists); + + + - compute_neighbors_orthogonally(tree.root(), distance_to_root); + compute_neighbors_orthogonally(tree.root(), distance_to_root,dists); if (sorted) this->queue.sort(); } private: - void compute_neighbors_orthogonally(typename Base::Node_const_handle N, FT rd) + void compute_neighbors_orthogonally(typename Base::Node_const_handle N, FT rd,std::vector& dists) { if (!(N->is_leaf())) { @@ -67,43 +74,46 @@ private: static_cast(N); this->number_of_internal_nodes_visited++; int new_cut_dim=node->cutting_dimension(); - FT old_off, new_rd; - FT new_off = *(query_object_it + new_cut_dim) - node->cutting_value(); - if ( ((new_off < FT(0.0)) && (this->search_nearest)) || - ((new_off >= FT(0.0)) && (!this->search_nearest)) ) + Base::Node_const_handle bestChild, otherChild; + FT new_off; + FT val = *(query_object_it + new_cut_dim); + FT diff1 = val - node->high_value(); + FT diff2 = val - node->low_value(); + if ( ((diff1 + diff2 < FT(0.0)) && (this->search_nearest)) || + ((diff1 + diff2 >= FT(0.0)) && (!this->search_nearest)) ) { - compute_neighbors_orthogonally(node->lower(),rd); if (this->search_nearest) { - old_off= *(query_object_it + new_cut_dim) - node->low_value(); - if (old_off>FT(0.0)) - old_off=FT(0.0); + new_off = diff1; + bestChild = node->lower(); + otherChild = node->upper(); } - else { - old_off= *(query_object_it + new_cut_dim) - node->high_value(); - if (old_offdistance_instance.new_distance(rd,old_off,new_off,new_cut_dim); - if (this->branch(new_rd)) - compute_neighbors_orthogonally(node->upper(), new_rd); + else { + new_off= diff2; + bestChild = node->upper(); + otherChild = node->lower(); + } } else // compute new distance { - compute_neighbors_orthogonally(node->upper(),rd); if (this->search_nearest) { - old_off= node->high_value() - *(query_object_it + new_cut_dim); - if (old_off>FT(0.0)) - old_off=FT(0.0); + new_off= diff2; + bestChild = node->upper(); + otherChild = node->lower(); } - else { - old_off= node->low_value() - *(query_object_it + new_cut_dim); - if (old_offlower(); + otherChild = node->upper(); } - new_rd = this->distance_instance. new_distance(rd,old_off,new_off,new_cut_dim); - if (this->branch(new_rd)) - compute_neighbors_orthogonally(node->lower(), new_rd); + } + compute_neighbors_orthogonally(bestChild, rd,dists); + FT dst=dists[new_cut_dim]; + FT new_rd = this->distance_instance.new_distance(rd,dst,new_off,new_cut_dim); + dists[new_cut_dim]=new_off; + if (this->branch(new_rd)) + compute_neighbors_orthogonally(otherChild, new_rd,dists); + dists[new_cut_dim]=dst; } else { diff --git a/Spatial_searching/include/CGAL/Search_traits_adapter.h b/Spatial_searching/include/CGAL/Search_traits_adapter.h index 315a093d7d2..3a0448888ad 100644 --- a/Spatial_searching/include/CGAL/Search_traits_adapter.h +++ b/Spatial_searching/include/CGAL/Search_traits_adapter.h @@ -157,11 +157,22 @@ public: return static_cast(this)->min_distance_to_rectangle(p,b); } + template + FT min_distance_to_rectangle(const Query_item& p, const CGAL::Kd_tree_rectangle& b,std::vector& dists) + { + return static_cast(this)->min_distance_to_rectangle(p,b,dists); + } + template FT max_distance_to_rectangle(const Query_item& p,const CGAL::Kd_tree_rectangle& b) const { return static_cast(this)->max_distance_to_rectangle(p,b); } + template + FT max_distance_to_rectangle(const Query_item& p,const CGAL::Kd_tree_rectangle& b,std::vector& dists) + { + return static_cast(this)->max_distance_to_rectangle(p,b,dists); + } }; }//namespace CGAL diff --git a/Spatial_searching/include/CGAL/Weighted_Minkowski_distance.h b/Spatial_searching/include/CGAL/Weighted_Minkowski_distance.h index 567c8fc8e6d..d53ed807ad8 100644 --- a/Spatial_searching/include/CGAL/Weighted_Minkowski_distance.h +++ b/Spatial_searching/include/CGAL/Weighted_Minkowski_distance.h @@ -178,6 +178,49 @@ namespace CGAL { return distance; } + inline + FT + min_distance_to_rectangle(const Query_item& q, + const Kd_tree_rectangle& r,std::vector& dists) { + FT distance = FT(0); + 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)) + { + for (unsigned int i = 0; qit != qe; ++qit, ++i) { + if (the_weights[i]*(r.min_coord(i) - + (*qit)) > distance){ + dists[i] = (r.min_coord(i)- + (*qit)); + distance = the_weights[i] * dists[i]; + } + if (the_weights[i] * ((*qit) - r.max_coord(i)) > + distance){ + dists[i] = + ((*qit)-r.max_coord(i)); + distance = the_weights[i] * dists[i]; + } + } + } + else + { + for (unsigned int i = 0; qit != qe; ++qit, ++i) { + if ((*qit) < r.min_coord(i)){ + dists[i] = r.min_coord(i)-(*qit); + distance += the_weights[i] * + std::pow(dists[i],power); + } + if ((*qit) > r.max_coord(i)){ + dists[i] = (*qit)-r.max_coord(i); + distance += the_weights[i] * + std::pow(dists[i],power); + } + } + }; + return distance; + } + inline FT max_distance_to_rectangle(const Query_item& q, @@ -214,6 +257,51 @@ namespace CGAL { }; return distance; } + + inline + FT + max_distance_to_rectangle(const Query_item& q, + const Kd_tree_rectangle& r,std::vector& dists) { + FT distance=FT(0); + 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)) + { + for (unsigned int i = 0; qit != qe; ++qit, ++i) { + if ((*qit) >= (r.min_coord(i) + + r.max_coord(i))/FT(2.0)) { + if (the_weights[i] * ((*qit) - + r.min_coord(i)) > distance){ + dists[i] = (*qit)-r.min_coord(i); + distance = the_weights[i] * + (dists[i]); + } + else + if (the_weights[i] * + (r.max_coord(i) - (*qit)) > distance){ + dists[i] = r.max_coord(i)-(*qit); + distance = the_weights[i] * + (dists[i]); + } + } + } + } + else + { + 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); + } + else{ + dists[i] = (*qit)-r.min_coord(i); + distance += the_weights[i] * std::pow(dists[i],power); + } + } + }; + return distance; + } inline FT diff --git a/Spatial_searching/test/Spatial_searching/Distance.h b/Spatial_searching/test/Spatial_searching/Distance.h index 04f270f06a3..772621b864b 100644 --- a/Spatial_searching/test/Spatial_searching/Distance.h +++ b/Spatial_searching/test/Spatial_searching/Distance.h @@ -23,6 +23,39 @@ struct Distance { if (h > b.max_coord(2)) distance += (h-b.max_coord(2))*(h-b.max_coord(2)); return distance; } + + template + double min_distance_to_rectangle(const Point& p, + const CGAL::Kd_tree_rectangle& b,std::vector& dists){ + double distance(0.0), h = p.x(); + if (h < b.min_coord(0)){ + dists[0] = (b.min_coord(0)-h); + distance += dists[0]*dists[0]; + } + if (h > b.max_coord(0)){ + dists[0] = (h-b.max_coord(0)); + distance += dists[0]*dists[0]; + } + h=p.y(); + if (h < b.min_coord(1)){ + dists[1] = (b.min_coord(1)-h); + distance += dists[1]*dists[1]; + } + if (h > b.max_coord(1)){ + dists[1] = (h-b.max_coord(1)); + distance += dists[1]*dists[1]; + } + h=p.z(); + if (h < b.min_coord(2)){ + dists[2] = (b.min_coord(2)-h); + distance += dists[2]*dists[2]; + } + if (h > b.max_coord(2)){ + dists[2] = (h-b.max_coord(2)); + distance += dists[2]*dists[2]; + } + return distance; + } template double max_distance_to_rectangle(const Point& p, @@ -40,6 +73,23 @@ struct Distance { (h-b.min_coord(2))*(h-b.min_coord(2)) : (b.max_coord(2)-h)*(b.max_coord(2)-h); return d0 + d1 + d2; } + + template + double max_distance_to_rectangle(const Point& p, + const CGAL::Kd_tree_rectangle& b,std::vector& dists){ + double h = p.x(); + + dists[0] = (h >= (b.min_coord(0)+b.max_coord(0))/2.0) ? + (h-b.min_coord(0)) : (b.max_coord(0)-h); + + h=p.y(); + dists[1] = (h >= (b.min_coord(1)+b.max_coord(1))/2.0) ? + (h-b.min_coord(1)) : (b.max_coord(1)-h); + h=p.z(); + dists[2] = (h >= (b.min_coord(2)+b.max_coord(2))/2.0) ? + (h-b.min_coord(2)) : (b.max_coord(2)-h); + return dists[0] * dists[0] + dists[1] * dists[1] + dists[2] * dists[2]; + } double new_distance(double& dist, double old_off, double new_off, int /*cutting_dimension*/) const {