diff --git a/Packages/Spatial_searching/examples/Spatial_searching/Example01.C b/Packages/Spatial_searching/examples/Spatial_searching/Example01.C new file mode 100755 index 00000000000..c21fe6c7ae5 --- /dev/null +++ b/Packages/Spatial_searching/examples/Spatial_searching/Example01.C @@ -0,0 +1,72 @@ +// Example illustrating for each separate splitting rule +// building a kd-tree + +#include +#include +#include +#include + +#include +#include + +#include +#include +#include +#include +#include +#include + +// use CGAL Point type +#include +#include + +typedef CGAL::Cartesian R; +typedef R::Point_3 Point; + +typedef CGAL::Plane_separator Separator; +typedef CGAL::Kd_tree_traits_point Traits; +typedef CGAL::Creator_uniform_3 Creator; + +int generate_kd_tree(CGAL::Split_rule_enumeration::Split_rule s) { + + int bucket_size=10; + + const int data_point_number=1000; + + typedef std::list point_list; + point_list data_points; + + CGAL::Random_points_in_cube_3 g( 1.0); + CGAL::copy_n( g, data_point_number, std::back_inserter(data_points)); + + Traits tr(bucket_size, s, 3.0, true); + typedef CGAL::Kd_tree Tree; + Tree d(data_points.begin(), data_points.end(), tr); + + std::cout << "created kd tree using splitting rule " << s << " containing " + << data_point_number << " points. " << std::endl; + d.statistics(); + + return 0; +}; + +int main() { + + generate_kd_tree(CGAL::Split_rule_enumeration::MEDIAN_OF_MAX_SPREAD); + generate_kd_tree(CGAL::Split_rule_enumeration::MEDIAN_OF_BOX); + generate_kd_tree(CGAL::Split_rule_enumeration::MIDPOINT_OF_MAX_SPREAD); + generate_kd_tree(CGAL::Split_rule_enumeration::MIDPOINT_OF_BOX); + generate_kd_tree(CGAL::Split_rule_enumeration::FAIR); + generate_kd_tree(CGAL::Split_rule_enumeration::SLIDING_MIDPOINT); + generate_kd_tree(CGAL::Split_rule_enumeration::SLIDING_FAIR); + + /* + double dummy; + std::cout << "Enter input to stop: \n" ; + std::cin >> dummy; + */ + + return 0; +}; + + diff --git a/Packages/Spatial_searching/include/CGAL/Orthogonal_distance_search.h b/Packages/Spatial_searching/include/CGAL/Orthogonal_distance_search.h new file mode 100755 index 00000000000..08eb239613a --- /dev/null +++ b/Packages/Spatial_searching/include/CGAL/Orthogonal_distance_search.h @@ -0,0 +1,439 @@ +// ====================================================================== +// +// Copyright (c) 2002 The CGAL Consortium +// +// This software and related documentation is part of an INTERNAL release +// of the Computational Geometry Algorithms Library (CGAL). It is not +// intended for general use. +// +// ---------------------------------------------------------------------- +// +// release : +// release_date : +// +// file : include/CGAL/Orthogonal_distance_search.h +// package : ASPAS +// revision : 1.4 +// revision_date : 2002/16/08 +// authors : Hans Tangelder () +// maintainer : Hans Tangelder () +// coordinator : Utrecht University +// +// ====================================================================== + +#ifndef ORTHOGONAL_DISTANCE_SEARCH +#define ORTHOGONAL_DISTANCE_SEARCH +#include +#include +#include +#include +#include + +namespace CGAL { + +template +class Orthogonal_distance_search { + +public: + +typedef typename Tree_traits::Item Item; +typedef typename Tree_traits::NT NT; +typedef Item** Item_iterator; +typedef Kd_tree_node Node; +typedef Kd_tree Tree; + +typedef typename Tree_traits::Item_with_distance Item_with_distance; +typedef std::pair Node_with_distance; + +private: + +// this forward declaration may cause problems for g++ +class iterator; + + + + typedef std::vector Node_with_distance_vector; + + typedef std::vector Item_with_distance_vector; + + typedef std::vector Distance_vector; + + iterator *start; + iterator *past_the_end; + + public: + + // constructor + Orthogonal_distance_search(Tree& tree, + Orthogonal_Distance& tr, Item& q, bool Standard_search, NT Eps = 0.0, + bool Search_nearest=true) + { + if (!Standard_search) { + start = new iterator(tree,q,tr,Eps,Search_nearest); + past_the_end = new iterator(); + } + search_nearest_neighbour=Search_nearest; + }; + + // destructor + Orthogonal_distance_search() { + delete start; + delete past_the_end; + }; + + iterator begin() { + return *start; + } + + iterator end() { + return *past_the_end; + } + + class iterator { + + public: + + typedef std::input_iterator_tag iterator_category; + typedef Item_with_distance value_type; + typedef int distance_type; + + class Iterator_implementation; + Iterator_implementation *Ptr_implementation; + + + public: + + // default constructor + iterator() {Ptr_implementation=0;} + + int the_number_of_items_visited() { + return Ptr_implementation->number_of_items_visited; + } + + // constructor + iterator(Tree& tree, Item& q, Orthogonal_Distance& tr, NT eps=0.0, bool search_nearest=true){ + Ptr_implementation = + new Iterator_implementation(tree, q, tr, eps, search_nearest); + } + + // copy constructor + iterator(iterator& Iter) { + Ptr_implementation = Iter.Ptr_implementation; + if (Ptr_implementation != 0) Ptr_implementation->reference_count++; + } + + Item_with_distance& operator* () { + return *(*Ptr_implementation); + } + + // prefix operator + iterator& operator++() { + ++(*Ptr_implementation); + return *this; + } + + // postfix operator + std::auto_ptr operator++(int) { + std::auto_ptr result = (*Ptr_implementation)++; + return result; + } + + bool operator==(const iterator& It) const { + + if ( + ((Ptr_implementation == 0) || + Ptr_implementation->Item_PriorityQueue.empty()) && + ((It.Ptr_implementation == 0) || + It.Ptr_implementation->Item_PriorityQueue.empty()) + ) + return true; + // else + return (Ptr_implementation == It.Ptr_implementation); + } + + bool operator!=(const iterator& It) const { + return !(*this == It); + } + + ~iterator() { + if (Ptr_implementation != 0) { + Ptr_implementation->reference_count--; + if (Ptr_implementation->reference_count==0) { + delete Ptr_implementation; + Ptr_implementation = 0; + } + } + } + + + class Iterator_implementation { + + public: + + int number_of_neighbours_computed; + int number_of_internal_nodes_visited; + int number_of_leaf_nodes_visited; + int number_of_items_visited; + + private: + + + NT multiplication_factor; + + Item* query_point; + + int total_item_number; + + NT distance_to_root; + + bool search_nearest_neighbour; + + NT rd; + + class Priority_higher + { + public: + + bool search_nearest; + + Priority_higher(bool search_the_nearest_neighbour) { + search_nearest = search_the_nearest_neighbour; + } + + //highest priority is smallest distance + bool operator() (Node_with_distance* n1, Node_with_distance* n2) const + { + if (search_nearest) { return (n1->second > n2->second);} + else {return (n2->second > n1->second);} + } + }; + + class Distance_smaller + { + + public: + + bool search_nearest; + + Distance_smaller(bool search_the_nearest_neighbour) { + // Search_traits s; + search_nearest = search_the_nearest_neighbour; + } + + //highest priority is smallest distance + bool operator() (Item_with_distance* p1, Item_with_distance* p2) const + { + if (search_nearest) {return (p1->second > p2->second);} + else {return (p2->second > p1->second);} + } + }; + + std::priority_queue* PriorityQueue; + + std::priority_queue* Item_PriorityQueue; + + Orthogonal_Distance* Orthogonal_Distance_instance; + + public: + + int reference_count; + + + + // constructor + Iterator_implementation(Tree& tree, Item& q, Orthogonal_Distance& tr, + NT Eps=0.0, bool search_nearest=true) + { + PriorityQueue= new std::priority_queue(Priority_higher(search_nearest)); + + Item_PriorityQueue= new std::priority_queue(Distance_smaller(search_nearest)); + + search_nearest_neighbour=search_nearest; + reference_count=1; + Orthogonal_Distance_instance=&tr; + multiplication_factor= + Orthogonal_Distance_instance->transformed_distance(1.0+Eps); + + if (search_nearest) distance_to_root= + Orthogonal_Distance_instance->lower_bound_distance_to_box(q, + *(tree.bounding_box())); + else distance_to_root= + Orthogonal_Distance_instance->upper_bound_distance_to_box(q, + *(tree.bounding_box())); + + + query_point = &q; + + total_item_number=tree.item_number(); + + number_of_leaf_nodes_visited=0; + number_of_internal_nodes_visited=0; + number_of_items_visited=0; + number_of_neighbours_computed=0; + + + + Node_with_distance *The_Root = new Node_with_distance(tree.root(), + distance_to_root); + PriorityQueue->push(The_Root); + + // rd is the distance of the top of the priority queue to q + rd=The_Root->second; + Compute_the_next_nearest_neighbour(); + } + + // * operator + Item_with_distance& operator* () { + return *(Item_PriorityQueue->top()); + } + + // prefix operator + Iterator_implementation& operator++() { + // std::cout << "called ++" << std::endl; + Delete_the_current_item_top(); + Compute_the_next_nearest_neighbour(); + return *this; + } + + // postfix operator + std::auto_ptr operator++(int) { + Item_with_distance Value = *(Item_PriorityQueue->top()); + std::auto_ptr + result(new Item_with_distance(Value)); + ++*this; + return result; + } + + + //destructor + ~Iterator_implementation() { + + // std::cout << "called destructor" << std::endl; + while (PriorityQueue->size()>0) { + Node_with_distance* The_top=PriorityQueue->top(); + PriorityQueue->pop(); + delete The_top; + }; + while (Item_PriorityQueue->size()>0) { + Item_with_distance* The_top=Item_PriorityQueue->top(); + Item_PriorityQueue->pop(); + delete The_top; + }; + delete PriorityQueue; + delete Item_PriorityQueue; + } + + private: + + void Delete_the_current_item_top() { + Item_with_distance* The_item_top=Item_PriorityQueue->top(); + Item_PriorityQueue->pop(); + delete The_item_top; + } + + void Compute_the_next_nearest_neighbour() { + + // compute the next item + bool next_neighbour_found=false; + if (!(Item_PriorityQueue->empty())) { + if (search_nearest_neighbour) + next_neighbour_found= + (multiplication_factor*rd > Item_PriorityQueue->top()->second); + else + next_neighbour_found= + (multiplication_factor*rd < Item_PriorityQueue->top()->second); + } + // otherwise browse the tree further + while ((!next_neighbour_found) && (!PriorityQueue->empty())) { + Node_with_distance* The_node_top=PriorityQueue->top(); + Node* N= The_node_top->first; + PriorityQueue->pop(); + delete The_node_top; + + while (!(N->is_leaf())) { // compute new distance + number_of_internal_nodes_visited++; + int new_cut_dim=N->separator()->cutting_dimension(); + NT old_off, new_rd; + NT new_off = + (*query_point)[new_cut_dim] - + N->separator()->cutting_value(); + if ( ((new_off < 0.0) && (search_nearest_neighbour)) || + (( new_off >= 0.0) && (!search_nearest_neighbour)) ) { + old_off= + (*query_point)[new_cut_dim]-N->low_value(); + if (old_off>0.0) old_off=0.0; + new_rd= + Orthogonal_Distance_instance-> + new_distance(rd,old_off,new_off,new_cut_dim); + Node_with_distance *Upper_Child = + new Node_with_distance(N->upper(),new_rd); + PriorityQueue->push(Upper_Child); + N=N->lower(); + + } + else { // compute new distance + old_off= N->high_value() - + (*query_point)[new_cut_dim]; + if (old_off>0.0) old_off=0.0; + new_rd=Orthogonal_Distance_instance-> + new_distance(rd,old_off,new_off,new_cut_dim); + Node_with_distance *Lower_Child = + new Node_with_distance(N->lower(),new_rd); + PriorityQueue->push(Lower_Child); + N=N->upper(); + } + } + // n is a leaf + number_of_leaf_nodes_visited++; + for (Item_iterator it=N->begin(); it != N->end(); it++) { + number_of_items_visited++; + NT distance_to_query_point= + Orthogonal_Distance_instance-> + distance(*query_point,**it); + Item_with_distance *NN_Candidate= + new Item_with_distance(*it,distance_to_query_point); + Item_PriorityQueue->push(NN_Candidate); + } + // old top of PriorityQueue has been processed, + // hence update rd + if (!PriorityQueue->empty()) { + rd = PriorityQueue->top()->second; + if (search_nearest_neighbour) + next_neighbour_found = + (multiplication_factor*rd > Item_PriorityQueue->top()->second); + else + next_neighbour_found = + (multiplication_factor*rd < Item_PriorityQueue->top()->second); + } + else // priority queue empty => last neighbour found + { + next_neighbour_found=true; + }; + number_of_neighbours_computed++; + } // next_neighbour_found or priority queue is empty + // in the latter case also the item priority quee is empty + } +}; // class Iterator_implementaion +}; // class iterator +}; // class Nearest neighbour_L2 + +template +void swap (typename Orthogonal_distance_search::iterator& x, +typename Orthogonal_distance_search::iterator& y) { +typename Orthogonal_distance_search::iterator::Iterator_implementation + *tmp = x.Ptr_implementation; + x.Ptr_implementation = y.Ptr_implementation; + y.Ptr_implementation = tmp; +} + +} // namespace CGAL + + +#endif // ORTHOGONAL_DISTANCE_SEARCH_H diff --git a/Packages/Spatial_searching/include/CGAL/Point_container.h b/Packages/Spatial_searching/include/CGAL/Point_container.h index 80514fe4d6d..bb835003049 100644 --- a/Packages/Spatial_searching/include/CGAL/Point_container.h +++ b/Packages/Spatial_searching/include/CGAL/Point_container.h @@ -34,13 +34,14 @@ namespace CGAL { template class Point_container; + /* problem with platform sparc_SunOS-5.6_CC-5.30 template struct comp_coord_val { int coord; comp_coord_val(const int i) : coord(i) {} bool operator() (const Item *a, const Item *b) { return (*a)[coord] < (*b)[coord]; } - }; + }; */ template std::ostream& operator<< (std::ostream& s, Point_container& c) { @@ -56,7 +57,7 @@ namespace CGAL { private: Point_list p_list; // list of pointers to points int built_coord; // a coordinate for which the pointer list is built - Kd_tree_rectangle bbox; // bounding box, i.e. cell of node + Kd_tree_rectangle bbox; // bounding box, i.e. rectangle of node Kd_tree_rectangle tbox; // tight bounding box, i.e. minimal enclosing bounding // box of points public: @@ -300,11 +301,18 @@ namespace CGAL { // assert(c.is_valid()); } - NT median(const int split_coord) { - /* if (p_list.empty()) { - // copy p_list[built_coord] to p_list[split_coord] - p_list[split_coord]=p_list[built_coord]; */ - // sort p_list[split_coord] +template + struct comp_coord_val { + int coord; + comp_coord_val(const int i) : coord(i) {} + bool operator() (const Item_ *a, const Item_ *b) { + return (*a)[coord] < (*b)[coord]; + } + }; + + + NT median(const int split_coord) { + p_list.sort(comp_coord_val(split_coord)); // } typename Point_list::iterator