// Copyright (c) 1997 Tel-Aviv University (Israel).
// All rights reserved.
//
// This file is part of CGAL (www.cgal.org); you may redistribute it under
// the terms of the Q Public License version 1.0.
// See the file LICENSE.QPL distributed with CGAL.
//
// Licensees holding a valid commercial license may use this file in
// accordance with the commercial license agreement provided with the software.
//
// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
//
// $Source: /CVSROOT/CGAL/Packages/kdtree/include/CGAL/kdtree_d.h,v $
// $Revision: 1.6 $ $Date: 2004/06/24 17:08:25 $
// $Name: $
//
// Author(s) : Sariel Har-Peled (sariel@math.tau.ac.il)
// Eyal Flato (flato@math.tau.ac.il)
#ifndef CGAL_KDTREE_D_H
#define CGAL_KDTREE_D_H
#include <cstring>
#include <list>
using std::list; // to avoid compiler crash on MSVC++
CGAL_BEGIN_NAMESPACE
/*=======================================================================
* Kdtree_interface -
* This is the default interface of point. It assume that PT (the
* point type have the following properties:
* default constructor
* int dimension()
* const coord_type & operator[]( int ) const
* ... operator=( const Pt & ) - copy operator
\*=======================================================================*/
template <class PT>
class Kdtree_interface
{
public:
typedef PT Point;
static int dimension( const PT & pnt )
{
// return pnt.dimensions();
return pnt.dimension();
}
static int compare( int d, const PT & a, const PT & b )
{
if ( a[ d ] < b[ d ] )
return -1;
if ( a[ d ] > b[ d ] )
return 1;
return 0;
}
static void copy_coord( int d, PT & a, const PT & b )
{
a[ d ] = b[ d ];
}
};
template <class PT>
class Kdtree_interface_2d
{
public:
typedef PT Point;
static int dimension( const PT & pnt )
{
// return pnt.dimensions();
return pnt.dimension();
}
static int compare( int d, const PT & a, const PT & b )
{
if ( a[ d ] < b[ d ] )
return -1;
if ( a[ d ] > b[ d ] )
return 1;
return 0;
}
static void copy_coord( int d, PT & a, const PT & b )
{
if ( d == 0 )
a = PT( b[ 0 ], a[1] );
else
if ( d == 1 )
a = PT( a[ 0 ], b[1] );
else {
assert( 0 );
}
}
};
template <class PT>
class Kdtree_interface_3d
{
public:
typedef PT Point;
static int dimension( const PT & pnt )
{
// return pnt.dimensions();
return pnt.dimension();
}
static int compare( int d, const PT & a, const PT & b )
{
if ( a[ d ] < b[ d ] )
return -1;
if ( a[ d ] > b[ d ] )
return 1;
return 0;
}
static void copy_coord( int d, PT & a, const PT & b )
{
if ( d == 0 )
a = PT( b[ 0 ], a[1], a[ 2 ] );
else
if ( d == 1 )
a = PT( a[ 0 ], b[1], a[ 2 ] );
else
if ( d == 2 )
a = PT( a[ 0 ], a[1], b[ 2 ] );
else {
assert( 0 );
}
}
};
/*=========================================================================
* kdtree_d -
* A the kdtree class.
*
* Remark: The kd-trees allocates all the memory it needs in advance,
* This results in a rather efficient memory management, and a fast
* iplementation.
\*=========================================================================*/
template <class Traits>
class Kdtree_d
{
class Plane;
public:
typedef typename Traits::Point Point;
typedef list<Point> List_points;
//-------------------------------------------------------------------------
// Extended_Point_d -
// A class for representing an extended d-dimal point: with
// ability to support +/- infinity. Templated by the kd-treee interface
// type.
//-------------------------------------------------------------------------
class ExtPoint
{
private:
class coordinate_type
{
public:
const Point * p_pnt;
int type;
//signed char type;
};
coordinate_type * p_arr;
int dim;
Point def_pnt;
void init( int _dim )
{
assert( _dim > 0 );
dim = _dim;
p_arr = (coordinate_type *)malloc( sizeof( coordinate_type )
* dim );
//printf( "p_arr(new): %p\n", (void *)p_arr );
assert( p_arr != NULL );
CGAL_CLIB_STD::memset( p_arr, 0, sizeof( coordinate_type ) * dim );
}
public:
enum { MINUS_INFINITY = -1, FINITE = 0, PLUS_INFINITY = 1 };
ExtPoint( int _type, int _dim )
{
assert( _type == MINUS_INFINITY || _type == PLUS_INFINITY );
init( _dim );
for ( int ind = 0; ind < dim; ind++ )
p_arr[ ind ].type = _type;
}
ExtPoint()
{
p_arr = NULL;
dim = -1;
}
ExtPoint( const ExtPoint & p )
{
init( p.dim );
def_pnt = p.def_pnt;
for ( int ind = 0; ind < dim; ind++ ) {
p_arr[ ind ] = p.p_arr[ ind ];
if ( p.p_arr[ ind ].p_pnt == &p.def_pnt )
p_arr[ ind ].p_pnt = &def_pnt;
}
}
ExtPoint & operator=( const ExtPoint & p )
{
term();
init( p.dim );
def_pnt = p.def_pnt;
for ( int ind = 0; ind < dim; ind++ ) {
p_arr[ ind ] = p.p_arr[ ind ];
if ( p.p_arr[ ind ].p_pnt == &p.def_pnt )
p_arr[ ind ].p_pnt = &def_pnt;
}
return *this;
}
ExtPoint( const Point & point, int dim )
{
init( dim );
def_pnt = point;
for ( int ind = 0; ind < dim; ind++ ) {
p_arr[ ind ].p_pnt = &def_pnt;
p_arr[ ind ].type = FINITE;
}
}
void term()
{
if ( p_arr != NULL ) {
//printf( "a: %p\n", p_arr );
free( p_arr );
//printf( "_a\n" );
p_arr = NULL;
}
dim = 0;
}
~ExtPoint()
{
term();
}
void set_coord( int k, Point & point )
{
assert( 0 <= k && k < dim );
p_arr[ k ].type = FINITE;
//p_arr[ k ].p_pnt = &point;
Traits::copy_coord( k, def_pnt, point );
p_arr[ k ].p_pnt = &def_pnt;
}
void set_coord( int k, const ExtPoint & point )
{
assert( 0 <= k && k < dim );
assert( 0 <= k && k < point.dim );
p_arr[ k ] = point.p_arr[ k ];
if ( p_arr[ k ].type == FINITE ) {
Traits::copy_coord( k, def_pnt, *(p_arr[ k ].p_pnt) );
p_arr[ k ].p_pnt = &def_pnt;
}
}
int compare( int k, const ExtPoint & point ) const
{
assert( 0 <= k && k < dim );
// the following does not compile on msvc++...
// coordinate_type & a( p_arr[ k ] ),
// & b( point.p_arr[ k ] );
coordinate_type & a = p_arr[ k ];
coordinate_type & b = point.p_arr[ k ];
if ( a.type != FINITE ) {
if ( b.type != FINITE ) {
return a.type - b.type;
} else {
return a.type;
}
} else {
if ( b.type != FINITE )
return -b.type;
else
return Traits::compare( k, *a.p_pnt, *b.p_pnt );
}
}
int compare_vector( const ExtPoint & point ) const
{
int ind, res;
for ( ind = 0; ind < dim; ind++ ) {
res = compare( ind, point );
if ( res != 0 )
return res;
}
return 0;
}
int compare( int k, const Point & point ) const
{
assert( 0 <= k && k < dim );
// coordinate_type & a( p_arr[ k ] );
coordinate_type & a = p_arr[ k ];
if ( a.type != FINITE )
return a.type;
else
return Traits::compare( k, *a.p_pnt, point );
}
int dimension() const
{
return dim;
}
int get_coord_status( int d ) const
{
assert( 0 <= d && d < dim );
return p_arr[ d ].type;
}
const Point * get_coord_point( int d ) const
{
assert( 0 <= d && d < dim );
return p_arr[ d ].p_pnt;
}
};
// Box - represents an axis parallel box.
class Box
{
public:
int dim;
ExtPoint left, right;
private:
friend class Plane;
ExtPoint & get_vertex( bool f_left )
{
return f_left? left : right;
}
public:
// constructors
//CHECK
Box()
{
}
Box( const Box & box )
{
dim = box.dim;
left = box.left;
right = box.right;
}
Box( const Point &l, const Point &r, int _dim )
{
dim = _dim;
left = ExtPoint( l, dim );
right = ExtPoint( r, dim );
}
Box( int _dim ) : left( ExtPoint::MINUS_INFINITY, _dim ),
right( ExtPoint::PLUS_INFINITY, _dim )
{
dim = _dim;
}
// data access
void set_left( Point &l)
{
left = ExtPoint( l, dim );
};
void set_right( Point &r)
{
right = ExtPoint( r, dim );
}
const ExtPoint &get_left() const
{
return left;
}
const ExtPoint &get_right() const
{
return right;
}
void set_coord_left( int k, Point & p )
{
left.set_coord( k, p );
}
void set_coord_right( int k, Point & p )
{
right.set_coord( k, p );
}
// operations
bool is_in( const Box &o) const
// checks if o is completely inside <this>
{
int dim = left.dimension();
for (int i = 0; i < dim; i++)
{
if ( (left.compare(i, o.get_left()) > 0 )
|| (right.compare(i, o.get_right()) < 0 ) )
return false;
}
return true;
}
bool is_in( const Point & o ) const
// checks if o is completely inside <this>
{
int dim = left.dimension();
for (int i = 0; i < dim; i++)
{
if ( (left.compare( i, o ) > 0 ) ||
(right.compare( i, o ) <= 0 ) )
return false;
}
return true;
}
bool is_coord_in_range( int k, const Point & o ) const
// checks if o is completely inside <this>
{
return ( ! ( (left.compare( k, o ) > 0 )
|| (right.compare( k, o ) <= 0 ) ) );
}
bool is_intersect( const Box &o) const
// checks if there is an intersection between o and this
{
int dim = left.dimension();
for (int i = 0; i < dim; i++)
{
if ( (left.compare(i, o.get_right()) >= 0) ||
(right.compare(i, o.get_left()) <= 0) )
return false;
}
return true;
}
// checks if there is an intersection between o and this box
// only in a specific coordinate...
bool is_intersect_in_dim( int d, const Box & o ) const
{
return (! ( (left.compare( d, o.get_right() ) >= 0 )
|| (right.compare( d, o.get_left() ) <= 0) ));
}
// checks if there is an intersection between o and this box
// only in a specific coordinate...
bool is_intersect_in_dim_closed( int d, const Box & o ) const
{
return (! ( (left.compare( d, o.get_right() ) > 0 )
|| (right.compare( d, o.get_left() ) < 0) ));
}
bool intersect(Box &o)
// intersects this with o. the intersection will be in this
// returns false if intersection is empty
{
int dim = left.dimension();
for (int i = 0; i < dim; i++)
{
// left is the maximal of the lefts
if (left.compare(i, o.get_left()) == -1)
left.set_coord(i, o.get_left());
// right is the minimal of the rights
if (right.compare(i, o.get_right()) == 1)
right.set_coord(i, o.get_right());
}
return !(is_empty());
}
bool is_empty() const
// return true if this is not an interval (left[k] > right[k])
{
int dim = left.dimension();
for (int i = 0; i < dim; i++)
{
if (left.compare(i, right) == 1)
return true;
}
return false;
}
bool is_empty_open() const
// return true if this is not an interval (left[k] > right[k])
{
int dim = left.dimension();
for (int i = 0; i < dim; i++)
{
if ( left.compare(i, right) >= 0 )
return true;
}
return false;
}
int comp( const Box & o ) const
{
int res;
res = left.compare_vector( o.left );
if ( res != 0 )
return res;
return right.compare_vector( o.right );
}
// destructor - needed in ...recursive ... DVP
~Box()
{
left.term();
right.term();
dim = 0;
}
};
private:
class Plane
{
private:
int coord;
Point * normal;
bool f_plus; // orientation of half space
// is (0, 0, ... , +inifinity, 0, ..., 0) inside plane
public:
Plane()
{
normal = NULL;
coord = 0;
}
Plane( int k, Point & p )
{
normal = &p;
coord = k;
}
Plane( const Plane & p )
{
coord = p.coord;
normal = p.normal;
}
void dump( void )
{
std::cout << "(" << coord << ": " << *normal << ")";
}
bool is_in( const Point & p ) const
{
int cmp;
cmp = Traits::compare( coord, p, *normal );
if ( ! f_plus )
cmp = -cmp;
return cmp >= 0;
}
void set_plane(int k, Point &p)
{
coord = k;
normal = &p;
//normal->copy( coord, p );
}
void split( Box & region, bool f_neg )
{
ExtPoint * p_p = &(region.get_vertex( ! f_neg ));
if ( f_neg ) {
if ( p_p->compare( coord, *normal ) > 0 )
p_p->set_coord( coord, *normal );
} else
if ( p_p->compare( coord, *normal ) < 0 )
p_p->set_coord( coord, *normal );
}
void orient_half_space( bool f_neg_side )
{
f_plus = ! f_neg_side;
}
int get_coord() const
{
return coord;
}
};
private:
class Node
{
public:
Plane plane;
Point * pnt;
Node * left, * right;
enum { LEFT, RIGHT };
const Plane & get_hs( int side ) const
{
((Plane *)&plane)->orient_half_space( side == LEFT );
return plane;
}
bool is_points_in_hs( const Plane & pl ) const
{
if ( is_point() )
return pl.is_in( *pnt );
if ( left != NULL && ( ! left->is_points_in_hs( pl ) ) )
return false;
if ( right != NULL && ( ! right->is_points_in_hs( pl ) ) )
return false;
return true;
}
bool is_valid() const
{
if ( is_point() )
return true;
if ( left != NULL )
if ( ! left->is_points_in_hs( get_hs( LEFT ) ) )
return false;
if ( right != NULL )
if ( ! right->is_points_in_hs( get_hs( RIGHT ) ) )
return false;
return true;
}
void dump( int depth )
{
int ind;
for ( ind = 0; ind < depth; ind++ )
std::cout << " ";
if ( is_point() ) {
std::cout << *pnt << "\n";
return;
}
plane.dump();
std::cout << "\n";
left->dump( depth + 1 );
for ( ind = 0; ind < depth; ind++ )
std::cout << " ";
std::cout << "!!!!!!!!!!!!\n";
right->dump( depth + 1 );
}
bool is_point() const
{
return ((left == NULL) && (right == NULL));
}
typedef std::back_insert_iterator<List_points> back_iter;
Node() : plane()
{
left = right = NULL;
}
void copy_subtree_points( back_iter & result,
const Box & rect )
{
if ( is_point() ) {
if ( rect.is_in( *pnt ) )
(*result++) = *pnt;
return;
}
if ( left != NULL )
left->copy_subtree_points( result, rect );
if ( right != NULL )
right->copy_subtree_points( result, rect );
}
static void search_recursive( back_iter & result,
Node * node,
const Box & rect,
Box & _region,
Plane & plane,
bool f_split_plus )
{
//printf( "search_recusrive\n" );
Box * p_r = new Box( _region );
//printf( "z" );
//fflush( stdout );
plane.split( *p_r, f_split_plus );
//printf( "c" );
//fflush( stdout );
assert( node != NULL );
//printf( "b" );
//fflush( stdout );
if ( rect.is_in( *p_r ) )
{
//printf( "5" );
//fflush( stdout );
node->copy_subtree_points( result, rect );
//printf( "\tsearch_recursive done...\n" );
//printf( "6" );
//fflush( stdout );
delete p_r;
return;
}
//printf( "v" );
//fflush( stdout );
if ( rect.is_intersect_in_dim_closed( plane.get_coord(), *p_r ) )
node->search( result, rect, *p_r );
//printf( "x" );
//fflush( stdout );
delete p_r;
}
void search( std::back_insert_iterator<List_points> result,
const Box &rect, Box ®ion )
{
if (is_point()) {
if ( rect.is_in( *pnt ) )
(*result++) = *pnt;
return;
}
//this is not a point so it is a hypeplane
if ( left != NULL )
search_recursive( result, left, rect,
region, plane, true );
if ( right != NULL )
search_recursive( result, right, rect,
region, plane, false );
}
};
// SunPro requires this :
friend class Box;
friend class Node;
typedef Point * Point_ptr;
int size;
Point * p_arr_pt;
Node *root;
int dim;
Node * p_node_arr;
int node_count;
Node * malloc_node( void )
{
Node * p_ret;
p_ret = &(p_node_arr[ node_count ]);
node_count++;
assert( node_count <= ( 2 * size ));
*p_ret = Node();
return p_ret;
}
static int comp( const Point & a, const Point & b, int dim )
{
return Traits::compare( dim, a, b );
}
static int partition( Point_ptr * arr, int left, int right,
Point * p_pivot, int dim )
{
int i, j;
Point_ptr tmp;
if ( left >= right )
return left;
i = left;
j = right;
while ( i < j ) {
if ( comp( *(arr[ i ]), *(arr[ j ]), dim ) > 0 ) {
tmp = arr[ i ];
arr[ i ] = arr[ j ];
arr[ j ] = tmp;
}
if ( comp( *(arr[ i ]), *p_pivot, dim ) < 0 ) {
i++;
} else
if ( comp( *p_pivot, *(arr[ j ]), dim ) <= 0 )
j--;
}
return (i > left)? i - 1 : left;
}
/* split the array into two sub-arrays, such that all the elements
* from left to pos_mid are smaller than the elements from pos+1 to
* right.
*/
static void split_arr( Point_ptr * arr,
int left,
int right,
int pos_mid,
int dim )
{
int pos;
if ( left >= right )
return;
pos = partition( arr, left, right, arr[ (left + right ) / 2 ],
dim );
if ( pos == pos_mid )
return;
if ( pos < pos_mid )
split_arr( arr, pos+1, right, pos_mid, dim );
else
split_arr( arr, left, pos, pos_mid, dim );
}
static Point_ptr get_max_element( Point_ptr * arr,
int left, int right, int d )
{
int max_pos = left;
Point mx = *(arr[ max_pos ]);
for ( int ind = left + 1; ind <= right; ind++ )
if ( comp( mx, *(arr[ ind ]), d ) < 0 ) {
mx = *(arr[ ind ]);
max_pos = ind;
}
return arr[ max_pos ];
}
Node *build_r( Point_ptr * arr, int left, int right,
int d )
{
int num, pos, next_d;
Node * n;
num = right - left + 1;
if ( num < 1)
return NULL;
// if the list contains only one point,
// construct a leaf for this node
if ( num == 1) {
//n = new node;
n = malloc_node();
n->pnt = arr[ left ];
return n;
}
// else divide space into two regions in
// dim-dim and cotinue recursively
pos = (left + right) / 2;
split_arr( arr, left, right, pos, d );
Point * p_median = get_max_element( arr, left, pos, d );
// create division plane;
Plane plane( d, *p_median );
//n = new node;
// assert( n != NULL );
n = malloc_node();
n->plane = plane;
next_d = d + 1;
if ( next_d >= dim )
next_d = 0;
// build left sub-tree
n->left = build_r( arr, left, pos, next_d );
// build right sub-tree
n->right = build_r( arr, pos + 1, right, next_d );
return n;
}
public:
typedef list<Point> list_points;
Kdtree_d(int k = 2)
{
dim = k;
root = NULL;
p_arr_pt = NULL;
p_node_arr = NULL;
}
~Kdtree_d()
{
delete_all();
}
bool is_valid( bool verbose = false, int level = 0 ) const
{
(void)verbose;
(void)level;
if ( root == NULL )
return true;
return root->is_valid();
}
void dump()
{
root->dump( 0);
}
void delete_all()
{
root = NULL;
if ( p_arr_pt != NULL )
delete[] p_arr_pt;
p_arr_pt = NULL;
if ( p_node_arr != NULL )
delete[] p_node_arr;
p_node_arr = NULL;
}
void search( std::back_insert_iterator<list_points> result,
Box & rect )
{
if (root == NULL)
return; // it is an empty tree - nothing to search in
Box region = Box( dim );
root->search( result, rect, region );
}
void build(list<Point> &l)
{
int i;
Point_ptr * p_arr;
size = l.size();
p_arr_pt = new Point[ size ];
assert( p_arr_pt != NULL );
p_arr = new Point_ptr[ size ];
assert( p_arr != NULL );
p_node_arr = new Node[ 2 * size ];
assert( p_node_arr != NULL );
node_count = 0;
/* fill the array */
i = 0;
for ( typename std::list<Point>::iterator j = l.begin();
j != l.end();
j++ )
{
p_arr_pt[ i ] = (*j);
p_arr[ i ] = p_arr_pt + i;
i++;
}
// recursively build the tree from the sorted list
// starting to divide it in coordinate 0
root = build_r( p_arr, 0, size - 1, 0 );
//printf( "b\n" );
delete[] p_arr;
}
};
CGAL_END_NAMESPACE
#else /* CGAL_KDTREE_D_H */
#error Header file kdtree_d.h included twice
#endif /* CGAL_KDTREE_D_H */
/*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*
*
* kdtree_d.h - End of File
\*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*/
syntax highlighted by Code2HTML, v. 0.9.1