Space-efficient geometric algorithms and data structures By Ilya Katz and Hervé Brönnimann

static_kd_tree.hpp

Go to the documentation of this file.

#ifndef STATIC_KD_TREE_HPP
#define STATIC_KD_TREE_HPP

#include <iostream>
#include <algorithm>

#include <boost/array.hpp>
#include "cgl.hpp"

namespace inplaceds
{
        
        typedef cgl<int>::point_type point_type;
        
        #define STATIC_KD_TREE_CUTOFF  16
        
        template <class RandomAccessIterator, class AxisAlignedGeometry>
        void
        build_static_kdtree(RandomAccessIterator first, RandomAccessIterator last, AxisAlignedGeometry geom)
        {
          //stop recursing if minimum size is reached
          if(last - first < STATIC_KD_TREE_CUTOFF)
          {
            return;
          }
        
          RandomAccessIterator mid = first + (last-first+1)/2 ;
          RandomAccessIterator Lmid = first + (mid-first+1)/2;
          RandomAccessIterator Rmid = mid + (last-mid+1)/2;
        
        
          //split by x coordinate
          std::nth_element (first, mid, last, geom.less_x_object() );
        
        
          //split left half by y-coordniate
          std::nth_element (first, Lmid ,mid, geom.less_y_object() );
        
          //split right half by y-coordinate
          std::nth_element (mid+1, Rmid, last, geom.less_y_object() );
        
          //----recursive calls----
        
          //recursively partition Lbot
          build_static_kdtree(first, Lmid, geom);
        
          //recursively partition Ltop
          build_static_kdtree(Lmid+1, mid, geom);
        
          //recursively partition Rbot
          build_static_kdtree(mid+1, Rmid, geom);
        
          //recursively partition Rtop
          build_static_kdtree(Rmid+1, last, geom);
        
        }
         
        template <class Box, class RandomAccessIterator, class OutputIterator, class AxisAlignedGeometry>
        OutputIterator
        window_query_static_kdtree(RandomAccessIterator first, RandomAccessIterator last, AxisAlignedGeometry geom, Box Q, OutputIterator result)
         {
          //if reached the lower limit, check each point individually
          if(last-first<STATIC_KD_TREE_CUTOFF)
          {
                for(RandomAccessIterator it = first; it!=last; it++)
                {
                        if(geom.point_in_box(*it,Q))
                        {
                                result = std::copy(it, it+1, result);
                        }
                }
            return result;
          }
                
          RandomAccessIterator mid = first + (last-first+1)/2 ;
          RandomAccessIterator Lmid = first + (mid-first+1)/2;
          RandomAccessIterator Rmid = mid + (last-mid+1)/2;

          //in all of the following cases we need to use >= or <= because
          //it is possible that points that have the same coordinates
          //appear in different regions (eg. Lbot and Rbot) since the division
          //into region is done by sizes of each region
          //simple example, is when all points in input are the same
        
        
          //check if need to search quadrant Ltop
          //if the left x coordinate of Q is to the left of intersection point
          //and the y coordinate is above the intersection point
          if(geom.get_x_coordinate(*mid)>=geom.get_left(Q) && geom.get_y_coordinate(*Lmid)<=geom.get_top(Q))
          {
            result = window_query_static_kdtree(Lmid+1, mid, geom, Q, result);
          }
        
          //check if need to search quadrant Rtop
          //if the right x coordniate of Q is to the right of intersection point
          //and the top y coordinate is above the intersection point
          if(geom.get_x_coordinate(*mid)<= geom.upper_right(Q)[0] && geom.get_y_coordinate(*Rmid)<=geom.upper_right(Q)[1])
          {
            result = window_query_static_kdtree(Rmid+1, last, geom, Q, result);
          }
        
          //check if need to search quadrant Lbottom
          //if the left x coordniate of Q is to the left of intersection point
          //and the bottom y coordinate is below the interscation point
          if(geom.get_x_coordinate(*mid)>=geom.get_left(Q) && geom.get_y_coordinate(*Lmid)>=geom.get_bottom(Q))
          {
            result = window_query_static_kdtree(first, Lmid, geom,Q, result);
          }
        
          //check if need to search quadrant Rbottom
          //if the right of Q is to the right of intersection point
          //and the bottom of Q is below the intersection point
          if(geom.get_x_coordinate(*mid)<=geom.get_right(Q) && geom.get_y_coordinate(*Rmid)>=geom.get_bottom(Q))
          {
            result = window_query_static_kdtree(mid+1, Rmid, geom,Q, result);
          }
          
          //mid, Lmid, and Rmid are considered independently
          //because these points are the immovable roots of the corresponding regions
          if(geom.point_in_box(*mid, Q))
          {
                result = std::copy(mid, mid+1, result);
          }
          if(geom.point_in_box(*Rmid, Q))
          {
                result = std::copy(Rmid, Rmid+1, result);
          }
          if(geom.point_in_box(*Lmid, Q))
          {
                result = std::copy(Lmid, Lmid+1, result);
          }
        
          return result;
        }
        
        template <class Halfplane, class RandomAccessIterator, class OutputIterator, class AxisAlignedGeometry>
        OutputIterator
        halfplane_query_static_kdtree(RandomAccessIterator first, RandomAccessIterator last,
                                                                 AxisAlignedGeometry geom, Halfplane H, OutputIterator result)
        {
         
          #define H_p1 geom.get_halfplane_points(H)[0]
          #define H_p2 geom.get_halfplane_points(H)[1]
         
        
         //if reached the lower limit, check each point individually
         
         
          if(last-first<STATIC_KD_TREE_CUTOFF)
          {
                for(RandomAccessIterator it = first; it!=last; it++)
                {
                                                //cout<<geom.get_x_coordinate(*it)<<","<<       geom.get_y_coordinate(*it)<<") - "<<it-first;
                        
                        if(geom.point_in_halfplane(*it,H))
                        {
                                result = std::copy(it, it+1, result);
                                //cout<<"#";
                        }
                        //cout<<endl;
                }
                //cout<<endl;
            return result;
          }
        
          RandomAccessIterator mid = first + (last-first+1)/2 ;
          RandomAccessIterator Lmid = first + (mid-first+1)/2;
          RandomAccessIterator Rmid = mid + (last-mid+1)/2;
        
          #define left_intersection_point  geom.make_point(geom.get_x_coordinate(*mid), geom.get_y_coordinate(*Lmid))
          #define right_intersection_point geom.make_point(geom.get_x_coordinate(*mid), geom.get_y_coordinate(*Rmid))
                                 
          //check if need to search quadrant Ltop
          //if the point is not in the plane and the plane is poiting down
          //do not needs to check Ltop if the first point is lower then the
          //second points
          if( !(
                        !geom.point_in_halfplane(left_intersection_point, H) &&
                         geom.is_BOTTOM(H) &&
                         geom.less_y_object()(H_p1,H_p2)
                         )
                )
          {
            result = halfplane_query_static_kdtree(Lmid+1, mid, geom, H, result);
          }
          
          //check if need to search quadrant Rtop
          if(!(
                        !geom.point_in_halfplane(right_intersection_point, H) &&
                         geom.is_BOTTOM(H) &&
                         geom.less_y_object()(H_p2,H_p1)
                        )
                 )
          {
            result = halfplane_query_static_kdtree(Rmid+1, last, geom, H, result);
          }
        
          //check if need to search quadrant Lbottom
          if( !(
                        !geom.point_in_halfplane(left_intersection_point, H) &&
                         geom.is_TOP(H) &&
                         geom.less_y_object()(H_p2,H_p1)
                         )
                )
          {
            result = halfplane_query_static_kdtree(first+1, Lmid, geom,H, result);
          }
        
          //check if need to search quadrant Rbottom
          if(!(
                        !geom.point_in_halfplane(right_intersection_point, H) &&
                         geom.is_TOP(H) &&
                         geom.less_y_object()(H_p1,H_p2)
                  )
                 )
          {
            result = halfplane_query_static_kdtree(mid+1, Rmid, geom, H, result);
          }     
          //critical points are considered independently
          //because these points are the immovable roots of the corresponding regions
          if(geom.point_in_halfplane(*mid, H))
          {
                result = std::copy(mid, mid+1, result);
          }
          if(geom.point_in_halfplane(*first, H))
          {
                result = std::copy(first, first+1, result);
          }
          if(geom.point_in_halfplane(*Rmid, H))
          {
                result = std::copy(Rmid, Rmid+1, result);
          }
          if(geom.point_in_halfplane(*Lmid, H))
          {
                result = std::copy(Lmid, Lmid+1, result);
          }
        
          #undef H_p1
          #undef H_p2
          #undef left_intersection_point
          #undef right_intersection_point
                
          return result;        
        
        }
        
        template <class Box, class RandomAccessIterator, class OutputIterator, class AxisAlignedGeometry>
        OutputIterator
        closed_triangle_query_static_kdtree(RandomAccessIterator first, RandomAccessIterator last, AxisAlignedGeometry geom, Box Q, OutputIterator result)
         {
          //if reached the lower limit, check each point individually
          if(last-first<STATIC_KD_TREE_CUTOFF)
          {
                for(RandomAccessIterator it = first; it!=last; it++)
                {
                        if(geom.point_in_closed_triangle(*it,T))
                        {
                                result = std::copy(it, it+1, result);
                        }
                }
            return result;
          }
                
          RandomAccessIterator mid = first + (last-first+1)/2 ;
          RandomAccessIterator Lmid = first + (mid-first+1)/2;
          RandomAccessIterator Rmid = mid + (last-mid+1)/2;

          //check if need to search quadrant Ltop
          if()
          {
            result = closed_triangle_query_static_kdtree(Lmid, mid, geom, Q, result);
          }
        
          //check if need to search quadrant Rtop
          if()
          {
            result = window_query_static_kdtree(Rmid, last, geom, Q, result);
          }
        
          //check if need to search quadrant Lbottom
          if()
          {
            result = window_query_static_kdtree(first, Lmid, geom,Q, result);
          }
        
          //check if need to search quadrant Rbottom
          if()
          {
            result = window_query_static_kdtree(mid, Rmid, geom,Q, result);
          }
        
          return result;
        }
        
}//namespace
#endif

---

Code Documentation generated Using Doxygen

Copyright © Ilya Katz and Hervé Brönnimann, 2005, 2006.