Spatial Searching

- both nearest and furthest neighbor searching
- both exact and approximate searching
- (approximate) range searching
- (approximate) $$
*k*-nearest and $$*k*-furthest neighbor searching - (approximate) incremental nearest and incremental furthest neighbor searching
- query items representing points and spatial objects.

*Midpoint_of_rectangle*- This splitting rule cuts a rectangle through its midpoint orthogonal to the longest side.
*Midpoint_of_max_spread*-
This splitting rule cuts a rectangle through $$
*(Mind+Maxd)/2*orthogonal to the dimension with the maximum point spread $$*[Mind,Maxd]*. *Sliding_midpoint*- This is a modification of the midpoint of rectangle splitting rule. It first attempts to perform a midpoint of rectangle split as described above. If data points lie on both sides of the separating plane the sliding midpoint rule computes the same separator as the midpoint of rectangle rule. If the data points lie only on one side it avoids this by sliding the separator, computed by the midpoint of rectangle rule, to the nearest datapoint.
*Median_of_rectangle*- The splitting dimension is the dimension of the longest side of the rectangle. The splitting value is defined by the median of the coordinates of the data points along this dimension.
*Median_of_max_spread*- The splitting dimension is the dimension of the longest side of the rectangle. The splitting value is defined by the median of the coordinates of the data points along this dimension.
*Fair*- This splitting rule is a compromise between the median of rectangle splitting rule and the midpoint of rectangle splitting rule. This splitting rule maintains an upper bound on the maximal allowed ratio of the longest and shortest side of a rectangle (the value of this upper bound is set in the constructor of the fair splitting rule). Among the splits that satisfy this bound, it selects the one in which the points have the largest spread. It then splits the points in the most even manner possible, subject to maintaining the bound on the ratio of the resulting rectangles.
*Sliding_fair*- This splitting rule is a compromise between the fair splitting rule and the sliding midpoint rule. Sliding fair-split is based on the theory that there are two types of splits that are good: balanced splits that produce fat rectangles, and unbalanced splits provided the rectangle with fewer points is fat. Also, this splitting rule maintains an upper bound on the maximal allowed ratio of the longest and shortest side of a rectangle (the value of this upper bound is set in the constructor of the fair splitting rule). Among the splits that satisfy this bound, it selects the one one in which the points have the largest spread. It then considers the most extreme cuts that would be allowed by the aspect ratio bound. This is done by dividing the longest side of the rectangle by the aspect ratio bound. If the median cut lies between these extreme cuts, then we use the median cut. If not, then consider the extreme cut that is closer to the median. If all the points lie to one side of this cut, then we slide the cut until it hits the first point. This may violate the aspect ratio bound, but will never generate empty cells.

// file: examples/Spatial_searching/Nearest_neighbor_searching.C #include <CGAL/Simple_cartesian.h> #include <CGAL/point_generators_2.h> #include <CGAL/Orthogonal_k_neighbor_search.h> #include <CGAL/Search_traits_2.h> #include <list> typedef CGAL::Simple_cartesian<double> K; typedef K::Point_2 Point_d; typedef CGAL::Search_traits_2<K> TreeTraits; typedef CGAL::Orthogonal_k_neighbor_search<TreeTraits> Neighbor_search; typedef Neighbor_search::Tree Tree; int main() { const int N = 1; std::list<Point_d> points; points.push_back(Point_d(0,0)); Tree tree(points.begin(), points.end()); Point_d query(0,0); // Initialize the search structure, and search all N points Neighbor_search search(tree, query, N); // report the N nearest neighbors and their distance // This should sort all N points by increasing distance from origin for(Neighbor_search::iterator it = search.begin(); it != search.end(); ++it){ std::cout << it->first << " "<< sqrt(it->second) << std::endl; } return 0; }

// file: examples/Spatial_searching/Distance_browsing.C #include <CGAL/Simple_cartesian.h> #include <CGAL/Orthogonal_incremental_neighbor_search.h> #include <CGAL/Search_traits_2.h> typedef CGAL::Simple_cartesian<double> K; typedef K::Point_2 Point_d; typedef CGAL::Search_traits_2<K> TreeTraits; typedef CGAL::Orthogonal_incremental_neighbor_search<TreeTraits> NN_incremental_search; typedef NN_incremental_search::iterator NN_iterator; typedef NN_incremental_search::Tree Tree; // A functor that returns true, iff the x-coordinate of a dD point is not positive struct X_not_positive { bool operator()(const NN_iterator& it) { return ((*it).first)[0]<0; } }; // An iterator that only enumerates dD points with positive x-coordinate typedef CGAL::Filter_iterator<NN_iterator, X_not_positive> NN_positive_x_iterator; int main() { Tree tree; tree.insert(Point_d(0,0)); tree.insert(Point_d(1,1)); tree.insert(Point_d(0,1)); tree.insert(Point_d(10,110)); tree.insert(Point_d(45,0)); tree.insert(Point_d(0,2340)); tree.insert(Point_d(0,30)); Point_d query(0,0); NN_incremental_search NN(tree, query); NN_positive_x_iterator it(NN.end(), X_not_positive(), NN.begin()), end(NN.end(), X_not_positive()); std::cout << "The first 5 nearest neighbours with positive x-coord are: " << std::endl; for (int j=0; (j < 5)&&(it!=end); ++j,++it) std::cout << (*it).first << " at squared distance = " << (*it).second << std::endl; return 0; }

// file: examples/Spatial_searching/General_neighbor_searching.C #include <CGAL/Cartesian_d.h> #include <CGAL/point_generators_2.h> #include <CGAL/Manhattan_distance_iso_box_point.h> #include <CGAL/K_neighbor_search.h> #include <CGAL/Search_traits_2.h> typedef CGAL::Cartesian_d<double> K; typedef K::Point_d Point_d; typedef CGAL::Random_points_in_square_2<Point_d> Random_points_iterator; typedef K::Iso_box_d Iso_box_d; typedef K TreeTraits; typedef CGAL::Manhattan_distance_iso_box_point<TreeTraits> Distance; typedef CGAL::K_neighbor_search<TreeTraits, Distance> Neighbor_search; typedef Neighbor_search::Tree Tree; int main() { const int N = 1000; const int K = 10; Tree tree; Random_points_iterator rpg; for(int i = 0; i < N; i++){ tree.insert(*rpg++); } Point_d pp(0.1,0.1); Point_d qq(0.2,0.2); Iso_box_d query(pp,qq); Distance tr_dist; Neighbor_search N1(tree, query, K, 0.0, false); // eps=10.0, nearest=false std::cout << "For query rectange = [0.1,0.2]^2 " << std::endl << "The " << K << " approximate furthest neighbors are: " << std::endl; for (Neighbor_search::iterator it = N1.begin();it != N1.end();it++) { std::cout << " Point " << it->first << " at distance = " << tr_dist.inverse_of_transformed_distance(it->second) << std::endl; } return 0; }

// file: examples/Spatial_searching/Fuzzy_range_query.C #include <CGAL/Cartesian_d.h> #include <CGAL/point_generators_d.h> #include <CGAL/Kd_tree.h> #include <CGAL/Fuzzy_sphere.h> #include <CGAL/Fuzzy_iso_box.h> #include <CGAL/Search_traits_d.h> typedef CGAL::Cartesian_d<double> K; typedef K::Point_d Point_d; typedef CGAL::Search_traits_d<K> Traits; typedef CGAL::Random_points_in_iso_box_d<Point_d> Random_points_iterator; typedef CGAL::Counting_iterator<Random_points_iterator> N_Random_points_iterator; typedef CGAL::Kd_tree<Traits> Tree; typedef CGAL::Fuzzy_sphere<Traits> Fuzzy_sphere; typedef CGAL::Fuzzy_iso_box<Traits> Fuzzy_iso_box; int main() { const int D = 4; const int N = 1000; // generator for random data points in the square ( (-1000,-1000), (1000,1000) ) Random_points_iterator rpit(4, 1000.0); // Insert N points in the tree Tree tree(N_Random_points_iterator(rpit,0), N_Random_points_iterator(N)); // define range query objects double pcoord[D] = { 300, 300, 300, 300 }; double qcoord[D] = { 900.0, 900.0, 900.0, 900.0 }; Point_d p(D, pcoord, pcoord+D); Point_d q(D, qcoord, qcoord+D); Fuzzy_sphere fs(p, 700.0, 100.0); Fuzzy_iso_box fib(p, q, 100.0); std::cout << "points approximately in fuzzy range query" << std::endl; std::cout << "with center (300.0, 300.0, 300.0, 300.0)" << std::endl; std::cout << "and fuzzy radius <200.0,400.0> are:" << std::endl; tree.search(std::ostream_iterator<Point_d>(std::cout, "\n"), fs); std::cout << "points approximately in fuzzy range query "; std::cout << "[<200,4000>,<800,1000>]]^4 are:" << std::endl; tree.search(std::ostream_iterator<Point_d>(std::cout, "\n"), fib); return 0; }

struct Point { double vec[3]; Point() { vec[0]= vec[1] = vec[2] = 0; } Point (double x, double y, double z) { vec[0]=x; vec[1]=y; vec[2]=z; } double x() const { return vec[ 0 ]; } double y() const { return vec[ 1 ]; } double z() const { return vec[ 2 ]; } double& x() { return vec[ 0 ]; } double& y() { return vec[ 1 ]; } double& z() { return vec[ 2 ]; } bool operator==(const Point& p) const { return (x() == p.x()) && (y() == p.y()) && (z() == p.z()) ; } bool operator!=(const Point& p) const { return ! (*this == p); } }; //end of class namespace CGAL { template <> struct Kernel_traits<Point> { struct Kernel { typedef double FT; typedef double RT; }; }; } struct Construct_coord_iterator { const double* operator()(const Point& p) const { return static_cast<const double*>(p.vec); } const double* operator()(const Point& p, int) const { return static_cast<const double*>(p.vec+3); } };We have put the glue layer in this file as well, that is a class that allows to iterate over the Cartesian coordinates of the point, and a class to construct such an iterator for a point. We next need a distance class

struct Distance { typedef Point Query_item; double transformed_distance(const Point& p1, const Point& p2) const { double distx= p1.x()-p2.x(); double disty= p1.y()-p2.y(); double distz= p1.z()-p2.z(); return distx*distx+disty*disty+distz*distz; } template <class TreeTraits> double min_distance_to_rectangle(const Point& p, const CGAL::Kd_tree_rectangle<TreeTraits>& b) const { double distance(0.0), h = p.x(); if (h < b.min_coord(0)) distance += (b.min_coord(0)-h)*(b.min_coord(0)-h); if (h > b.max_coord(0)) distance += (h-b.max_coord(0))*(h-b.max_coord(0)); h=p.y(); if (h < b.min_coord(1)) distance += (b.min_coord(1)-h)*(b.min_coord(1)-h); if (h > b.max_coord(1)) distance += (h-b.max_coord(1))*(h-b.min_coord(1)); h=p.z(); if (h < b.min_coord(2)) distance += (b.min_coord(2)-h)*(b.min_coord(2)-h); if (h > b.max_coord(2)) distance += (h-b.max_coord(2))*(h-b.max_coord(2)); return distance; } template <class TreeTraits> double max_distance_to_rectangle(const Point& p, const CGAL::Kd_tree_rectangle<TreeTraits>& b) const { double h = p.x(); double d0 = (h >= (b.min_coord(0)+b.max_coord(0))/2.0) ? (h-b.min_coord(0))*(h-b.min_coord(0)) : (b.max_coord(0)-h)*(b.max_coord(0)-h); h=p.y(); double d1 = (h >= (b.min_coord(1)+b.max_coord(1))/2.0) ? (h-b.min_coord(1))*(h-b.min_coord(1)) : (b.max_coord(1)-h)*(b.max_coord(1)-h); h=p.z(); double d2 = (h >= (b.min_coord(2)+b.max_coord(2))/2.0) ? (h-b.min_coord(2))*(h-b.min_coord(2)) : (b.max_coord(2)-h)*(b.max_coord(2)-h); return d0 + d1 + d2; } double new_distance(double& dist, double old_off, double new_off, int cutting_dimension) const { return dist + new_off*new_off - old_off*old_off; } double transformed_distance(double d) const { return d*d; } double inverse_of_transformed_distance(double d) { return sqrt(d); } }; // end of struct Distance

We are ready to put the pices together.
The class *Search_traits<..>* which you see in the next file is then a mere
wrapper for all these types. The searching itself works exactly as for CGAL kernels.

//file: examples/Spatial_searching/User_defined_point_and_distance.C #include <CGAL/basic.h> #include <CGAL/Search_traits.h> #include <CGAL/point_generators_3.h> #include <CGAL/Orthogonal_k_neighbor_search.h> #include "Point.h" // defines types Point, Construct_coord_iterator #include "Distance.h" typedef CGAL::Random_points_in_cube_3<Point> Random_points_iterator; typedef CGAL::Counting_iterator<Random_points_iterator> N_Random_points_iterator; typedef CGAL::Search_traits<double, Point, const double*, Construct_coord_iterator> Traits; typedef CGAL::Orthogonal_k_neighbor_search<Traits, Distance> K_neighbor_search; typedef K_neighbor_search::Tree Tree; int main() { const int N = 1000; const int K = 5; // generator for random data points in the cube ( (-1,-1,-1), (1,1,1) ) Random_points_iterator rpit( 1.0); // Insert number_of_data_points in the tree Tree tree(N_Random_points_iterator(rpit,0), N_Random_points_iterator(N)); Point query(0.0, 0.0, 0.0); Distance tr_dist; // search K nearest neighbours K_neighbor_search search(tree, query, K); for(K_neighbor_search::iterator it = search.begin(); it != search.end(); it++){ std::cout << " d(q, nearest neighbor)= " << tr_dist.inverse_of_transformed_distance(it->second) << std::endl; } // search K furthest neighbour searching, with eps=0, search_nearest=false K_neighbor_search search2(tree, query, K, 0.0, false); for(K_neighbor_search::iterator it = search2.begin(); it != search2.end(); it++){ std::cout << " d(q, furthest neighbor)= " << tr_dist.inverse_of_transformed_distance(it->second) << std::endl; } return 0; }

// file: examples/Spatial_searching/Using_fair_splitting_rule.C #include <CGAL/Simple_cartesian.h> #include <CGAL/point_generators_2.h> #include <CGAL/Search_traits_2.h> #include <CGAL/Orthogonal_k_neighbor_search.h> typedef CGAL::Simple_cartesian<double> R; typedef R::Point_2 Point_d; typedef CGAL::Random_points_in_square_2<Point_d> Random_points_iterator; typedef CGAL::Counting_iterator<Random_points_iterator> N_Random_points_iterator; typedef CGAL::Search_traits_2<R> Traits; typedef CGAL::Euclidean_distance<Traits> Distance; typedef CGAL::Fair<Traits> Fair; typedef CGAL::Orthogonal_k_neighbor_search<Traits,Distance,Fair> Neighbor_search; typedef Neighbor_search::Tree Tree; int main() { const int N = 1000; // generator for random data points in the square ( (-1,-1), (1,1) ) Random_points_iterator rpit( 1.0); Fair fair(5); // bucket size=5 // Insert number_of_data_points in the tree Tree tree(N_Random_points_iterator(rpit,0), N_Random_points_iterator(N), fair); Point_d query(0,0); // Initialize the search structure, and search all N points Neighbor_search search(tree, query, N); // report the N nearest neighbors and their distance // This should sort all N points by increasing distance from origin for(Neighbor_search::iterator it = search.begin(); it != search.end(); ++it){ std::cout << it->first << " "<< sqrt(it->second) << std::endl; } return 0; }

Next chapter: Spatial Searching

The CGAL Project .
Tue, December 21, 2004 .