Simple questions on geometric primitives, such as intersection and distance computations, can themselves become quite expensive if the primitives are not so simple anymore, for example, three-dimensional triangles and facets of polyhedral surfaces. Thus algorithms operating on these primitives tend to be slow in practice. A common (heuristic) optimization approximates the geometric primitives with their axis-aligned bounding boxes, runs a suitable modification of the algorithm on the boxes, and whenever a pair of boxes has an interesting interaction, only then the exact answer is computed on the complicated geometric primitives contained in the boxes.
We provide an efficient algorithm [ZE02] for finding all
intersecting pairs for large numbers of iso-oriented boxes, i.e.,
typically these will be such bounding boxes of more complicated geometries.
One immediate application of this algorithm is the detection of all
intersections (and self-intersections) for polyhedral surfaces, i.e.,
applying the algorithm on a large set of triangles in space, we give
an example program later in this chapter. Not so obvious applications
are proximity queries and distance computations among such surfaces,
see Section 
for an example
and [ZE02] for more details.
A -dimensional iso-oriented box is defined as the Cartesian product of intervals. We call the box half-open if the intervals are half-open intervals, and we call the box closed if the intervals are closed intervals. Note that closed boxes support zero-width boxes and they can intersect at their boundaries, while non-empty half-open boxes always have a positive volume and they only intersect iff their interiors overlap. The distinction between closed and half-open boxes does not require a different representation of boxes, just a different interpretation when comparing boxes, which is selected with the two possible values for the topology parameter:
The number type of the interval boundaries must be one of the builtin types int, unsigned int, double or float.
In addition, a box has an unique id-number. It is used to order boxes consistently in each dimension even if boxes have identical coordinates. In consequence, the algorithm guarantees that a pair of intersecting boxes is reported only once. Note that boxes with equal id-number are not reported since they obviously intersect trivially.
The box intersection algorithm comes in two flavors: One algorithm works on a single sequence of boxes and computes all pairwise intersections, which is called the complete case, and used, for example, in the self-intersection test. The other algorithm works on two sequences of boxes and computes the pairwise intersections between boxes from the first sequence with boxes from the second sequence, which is called the bipartite case. For each pairwise intersection found a callback function is called with two arguments; the first argument is a box from the first sequence and the second argument a box from the second sequence. In the complete case, the second argument is a box from an internal copy of the first sequence.
The box intersection algorithm is implemented as a family of generic functions; the functions for the complete case accept one iterator range, and the functions for the bipartite case accept two iterator ranges. The callback function for reporting the intersecting pairs is provided as a template parameter of the BinaryFunction concept. The two principle function calls utilizing all default arguments look as follows:
#include <CGAL/box_intersection_d.h>
|
| ||||
|
|
| |||
|
| ||||
|
|
| |||
Additional parameters to the functions calls are a cutoff value to adjust performance tradeoffs, and a topology parameter selecting between topologically closed boxes (the default) and topologically half-open boxes.
The algorithm reorders the boxes in the course of the algorithm. Now, depending on the size of a box it can be faster to copy the boxes, or to work with pointers to boxes and copy only pointers. We offer automatic support for both options. To simplify the description, let us call the value_type of the iterator ranges box handle. The box handle can either be our box type itself or a pointer (or const pointer) to the box type; these choices represent both options from above.
In general, the algorithms treat the box type as opaque type and just assume that they are models of the Assignable concept, so that the algorithms can modify the input sequences and reorder the boxes. The access to the box dimension and box coordinates is mediated with a traits class of the BoxIntersectionTraits_d concept. A default traits class is provided that assumes that the box type is a model of the BoxIntersectionBox_d concept and that the box handle, i.e., the iterators value type, is identical to the box type or a pointer to the box type (see the previous paragraph for the value versus pointer nature of the box handle).
Two implementations of iso-oriented boxes are provided; CGAL::Box_intersection_d::Box_d as a plain box, and CGAL::Box_intersection_d::Box_with_handle_d as a box plus a handle that can be used to point to the full geometry that is approximated by the box. Both implementations have template parameters for the number type used for the interval bounds, for the fixed dimension of the box, and for a policy class [Ale01] selecting among several solutions for providing the id-number.
The function signatures for the bipartite case look as follows. The signatures for the complete case with the box_self_intersection_d function look the same except for the single iterator range.
#include <CGAL/box_intersection_d.h>
|
| ||||
|
|
| |||
|
| ||||
|
|
| |||
The box implementation provided with CGAL::Box_intersection_d::Box_d<double,2> has a dedicated constructor for the CGAL bounding box type CGAL::Bbox_2 (similar for dimension 3). We use this in our minimal example to create easily nine two-dimensional boxes in a grid layout of boxes. Additionally we pick the center box and the box in the upper-right corner as our second box sequence query.
The default policy of the box type implements the id-number with
an explicit counter in the boxes, which is the default choice since it
always works, but it costs space that could potentially be avoided,
see the example in the next section. We use the id-number in our
callback function to report the result of the intersection algorithm.
The result will be that the first query box intersects all nine
boxes and the second query box intersects the four boxes
in the upper-right quadrant. See Section
for the change of the topology parameter and its effect.
// file: examples/Box_intersection_d/minimal.C
#include <CGAL/box_intersection_d.h>
#include <CGAL/Bbox_2.h>
#include <iostream>
typedef CGAL::Box_intersection_d::Box_d<double,2> Box;
typedef CGAL::Bbox_2 Bbox;
// 9 boxes of a grid
Box boxes[9] = { Bbox( 0,0,1,1), Bbox( 1,0,2,1), Bbox( 2,0,3,1), // low
Bbox( 0,1,1,2), Bbox( 1,1,2,2), Bbox( 2,1,3,2), // middle
Bbox( 0,2,1,3), Bbox( 1,2,2,3), Bbox( 2,2,3,3)};// upper
// 2 selected boxes as query; center and upper right
Box query[2] = { Bbox( 1,1,2,2), Bbox( 2,2,3,3)};
void callback( const Box& a, const Box& b ) {
std::cout << "box " << a.id() << " intersects box " << b.id() << std::endl;
}
int main() {
CGAL::box_intersection_d( boxes, boxes+9, query, query+2, callback);
return 0;
}
The conventional application of the axis-aligned box intersection algorithm will start from complex geometry, here 3D triangles, approximate them with their bounding box, compute the intersecting pairs of boxes, and check only for those if the original triangles intersect as well.
We start in the main function and create ten triangles with endpoints chosen randomly in a cube . We store the triangles in a vector called triangles.
Next we create a vector for the bounding boxes of the triangles called boxes. For the boxes we choose the type Box_with_handle_d<double,3,Iterator> that works nicely together with the CGAL bounding boxes of type CGAL::Bbox_3. In addition, each box stores the iterator to the corresponding triangle.
The default policy of this box type uses for the id-number the address of the value of the iterator, i.e., the address of the triangle. This is a good choice that works correctly iff the boxes have unique iterators, i.e., there is a one-to-one mapping between boxes and approximated geometry, which is the case here. It saves us the extra space that was needed for the explicit id-number in the previous example.
We run the self intersection algorithm with the report_inters function as callback. This callback reports the intersecting boxes. It uses the handle and the global triangles vector to calculate the triangle numbers. Then it checks the triangles themselves for intersection and reports if not only the boxes but also the triangles intersect. We take some precautions before the intersection test in order to avoid problems, although unlikely, with degenerate triangles that we might have created with the random process.
This example can be easily extended to test polyhedral surfaces of the Polyhedron_3 class for (self-) intersections. The main difference are the numerous cases of incidences between triangles in the polyhedral surface that should not be reported as intersections, see the examples/Polyhedron/polyhedron_self_intersection.C example program in the CGAL distribution.
// file: examples/Box_intersection_d/triangle_self_intersect.C
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/intersections.h>
#include <CGAL/point_generators_3.h>
#include <CGAL/Bbox_3.h>
#include <CGAL/box_intersection_d.h>
#include <CGAL/function_objects.h>
#include <CGAL/Join_input_iterator.h>
#include <CGAL/copy_n.h>
#include <vector>
typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
typedef Kernel::Point_3 Point_3;
typedef Kernel::Triangle_3 Triangle_3;
typedef std::vector<Triangle_3> Triangles;
typedef Triangles::iterator Iterator;
typedef CGAL::Box_intersection_d::Box_with_handle_d<double,3,Iterator> Box;
Triangles triangles; // global vector of all triangles
// callback function that reports all truly intersecting triangles
void report_inters( const Box& a, const Box& b) {
std::cout << "Box " << (a.handle() - triangles.begin()) << " and "
<< (b.handle() - triangles.begin()) << " intersect";
if ( ! a.handle()->is_degenerate() && ! b.handle()->is_degenerate()
&& CGAL::do_intersect( *(a.handle()), *(b.handle()))) {
std::cout << ", and the triangles intersect also";
}
std::cout << '.' << std::endl;
}
int main() {
// Create 10 random triangles
typedef CGAL::Random_points_in_cube_3<Point_3> Pts;
typedef CGAL::Creator_uniform_3< Point_3, Triangle_3> Creator;
typedef CGAL::Join_input_iterator_3<Pts,Pts,Pts,Creator> Triangle_gen;
Pts points( 1); // in centered cube [-1,1)^3
Triangle_gen triangle_gen( points, points, points);
CGAL::copy_n( triangle_gen, 10, std::back_inserter(triangles));
// Create the corresponding vector of bounding boxes
std::vector<Box> boxes;
for ( Iterator i = triangles.begin(); i != triangles.end(); ++i)
boxes.push_back( Box( i->bbox(), i));
// Run the self intersection algorithm with all defaults
CGAL::box_self_intersection_d( boxes.begin(), boxes.end(), report_inters);
return 0;
}
We modify the previous example, finding intersecting 3D triangles, and add an additional vector ptr that stores pointers to the bounding boxes, so that the intersection algorithm will work on a sequence of pointers and not on a sequence of boxes. The change just affects the preparation of the additional vector and the call of the box intersection function. The box intersection function (actually its default traits class) detects automatically that the value type of the iterators is a pointer type and not a class type.
// Create the corresponding vector of pointers to bounding boxes
std::vector<Box *> ptr;
for ( std::vector<Box>::iterator i = boxes.begin(); i != boxes.end(); ++i)
ptr.push_back( &*i);
// Run the self intersection algorithm with all defaults on the
// indirect pointers to bounding boxes. Avoids copying the boxes.
CGAL::box_self_intersection_d( ptr.begin(), ptr.end(), report_inters);
In addition, the callback function report_inters needs to be changed to work with pointers to boxes. See the following file for the full example program.
examples/Box_intersection_d/triangle_self_intersect_pointers.C
A note on performance: The algorithm sorts and partitions the input sequences. It is clearly costly to copy a large box compared to a simple pointer. However, the algorithm benefits from memory locality in the later stages when it copies the boxes, while the pointers would refer to boxes that become wildly scattered in memory. These two effects, copying costs and memory locality, counteract each other. For small box sizes, i.e., small dimension, memory locality wins and one should work with boxes, while for larger box sizes one should work with pointers. The exact threshold depends on the memory hierarchy (caching) of the hardware platform and the size of the boxes, most notably the type used to represent the box coordinates. A concrete example; on a laptop with an Intel Mobile Pentium4 running at 1.80GHz with 512KB cache and 254MB main memory under Linux this version with pointers was 20% faster than the version above that copies the boxes for 10000 boxes, but the picture reversed for 100000 boxes, where the version above that copies the boxes becomes 300% faster.
Note that switching to the builtin type float is supported by the box intersection algorithm, but the interfacing with the CGAL bounding box CGAL::Bbox_3 would not be that easy. In particular, just converting from the double to the float representation incurs rounding that needs to be controlled properly, otherwise the box might shrink and one might miss intersections.
Boxes can be interpreted by the box intersection algorithm as closed
or as half-open boxes, see also Section . Closed
boxes support zero-width boxes and they can intersect at their
boundaries, while half-open boxes always have a positive volume and
they only intersect iff their interiors overlap. The choice between
closed or half-open boxes is selected with the topology
parameter and its two values:
The example program uses a two-dimensional box with int
coordinates and id-numbers that are by default explicitly
stored. We create the same boxes as in the minimal example in
Section . We create a grid
of boxes, and two boxes for the query sequence, namely the
box at the center and the box from the upper-right corner of the grid.
We write a more involved callback function object Report that stores an output iterator and writes the id-number of the box in the first argument to the output iterator. We also provide a small helper function report that simplifies the use of the function object.
We call the box intersection algorithm twice; once for the default topology, which is the closed box topology, and once for the half-open box topology. We sort the resulting output for better readability and verify its correctness with the check1 and check2 data. For the closed box topology, the center box in query intersects all boxes, and the upper-right box in query intersects the four boxes of the upper-right quadrant in boxes. Almost all intersections are with the box boundaries, thus, for the half-open topology only one intersection remains per query box, namely its corresponding box in boxes. So, the output of the algorithm will be:
0 1 2 3 4 4 5 5 6 7 7 8 8
4 8
For the second box intersection function call we have to specify the
cutoff parameter explicitly. See the
Section below for a detailed
discussion.
// file: examples/Box_intersection_d/box_grid.C
#include <CGAL/box_intersection_d.h>
#include <vector>
#include <algorithm>
#include <iterator>
#include <assert.h>
typedef CGAL::Box_intersection_d::Box_d<int,2> Box;
// coordinates for 9 boxes of a grid
int p[9*4] = { 0,0,1,1, 1,0,2,1, 2,0,3,1, // lower
0,1,1,2, 1,1,2,2, 2,1,3,2, // middle
0,2,1,3, 1,2,2,3, 2,2,3,3};// upper
// 9 boxes
Box boxes[9] = { Box( p, p+ 2), Box( p+ 4, p+ 6), Box( p+ 8, p+10),
Box( p+12, p+14), Box( p+16, p+18), Box( p+20, p+22),
Box( p+24, p+26), Box( p+28, p+30), Box( p+32, p+34)};
// 2 selected boxes as query; center and upper right
Box query[2] = { Box( p+16, p+18), Box( p+32, p+34)};
// callback function object writing results to an output iterator
template <class OutputIterator>
struct Report {
OutputIterator it;
Report( OutputIterator i) : it(i) {} // store iterator in object
// We write the id-number of box a to the output iterator assuming
// that box b (the query box) is not interesting in the result.
void operator()( const Box& a, const Box&) { *it++ = a.id(); }
};
template <class Iter> // helper function to create the function object
Report<Iter> report( Iter it) { return Report<Iter>(it); }
int main() {
// run the intersection algorithm and store results in a vector
std::vector<std::size_t> result;
CGAL::box_intersection_d( boxes, boxes+9, query, query+2,
report( std::back_inserter( result)));
// sort, check, and show result
std::sort( result.begin(), result.end());
std::size_t check1[13] = {0,1,2,3,4,4,5,5,6,7,7,8,8};
assert(result.size() == 13 && std::equal(check1,check1+13,result.begin()));
std::copy( result.begin(), result.end(),
std::ostream_iterator<std::size_t>( std::cout, " "));
std::cout << std::endl;
// run it again but for different cutoff value and half-open boxes
result.clear();
CGAL::box_intersection_d( boxes, boxes+9, query, query+2,
report( std::back_inserter( result)),
std::ptrdiff_t(1),
CGAL::Box_intersection_d::HALF_OPEN);
// sort, check, and show result
std::sort( result.begin(), result.end());
std::size_t check2[2] = {4,8};
assert(result.size() == 2 && std::equal(check2, check2+2, result.begin()));
std::copy( result.begin(), result.end(),
std::ostream_iterator<std::size_t>( std::cout, " "));
std::cout << std::endl;
return 0;
}
The implemented algorithm is described in [ZE02] as version two. Its performance depends on a cutoff parameter. When the size of both iterator ranges drops below the cutoff parameter the function switches from the streamed segment-tree algorithm to the two-way-scan algorithm, see [ZE02] for the details.
The streamed segment-tree algorithm needs log worst-case running time and space, where is the number of boxes in both input sequences, the (constant) dimension of the boxes, and the output complexity, i.e., the number of pairwise intersections of the boxes. The two-way-scan algorithm needs log worst-case running time and space, where is the number of pairwise overlapping intervals in one dimensions (the dimension where the algorithm is used instead of the segment tree). Note that is not necessarily related to and using the two-way-scan algorithm is a heuristic.
Unfortunately, we have no general method to automatically determine an optimal cutoff parameter, since it depends on the used hardware, the runtime ratio between callback runtime and segment-tree runtime, and of course the number of boxes to be checked and their distribution. In cases where the callback runtime is dominant, it may be best to make the threshold parameter small. Otherwise a cutoff can lead to acceptable results. For well distributed boxes the original paper [ZE02] gives optimal cutoffs in the thousands. Anyway, for optimal runtime some experiments to compare different cutoff parameters are recommended.
To demonstrate that box intersection can be done quite fast, different box sequences are intersected in the range between 4 and 800000 boxes total. We use three-dimensional default boxes of closed topology with float coordinates and without additional data fields. The algorithm works directly on the boxes, not on pointer to boxes. Each box intersection is reported to an empty dummy callback.
For each box set, a near-optimal cutoff parameter is determined using
an adaptive approximation. The runtime required for streaming is
compared against usual scanning. Results on a Xeon 2.4GHz with 4GB
main memory can be seen in Figure . For a small
number of boxes, pure scanning is still faster than streaming with
optimal cutoff, which would just delegate the box sets to the scanning
algorithm. As there are more and more boxes, the overhead becomes less
important.
The example in the previous Section uses
an array to provide the coordinates and then creates another array for
the boxes. In the following example we write our own box class
Box that we can initialize directly with the four coordinates and
create the array of boxes directly. We also omit the explicitly stored
id-number and use the address of the box itself as
id-number. This works only if the boxes do not change their
position, i.e., we work with pointers to the boxes in the intersection
algorithm.
We follow with our own box class Box the BoxIntersectionBox_d concept, which allows us to reuse the default traits implementation, i.e., we can use the same default function call to compute all intersections. See the example in the next section for a self-written traits class. So, in principle, the remainder of the example stays the same and we omit the part from the previous example for brevity that illustrates the half-open box topology.
The requirements for the box implementation are best studied on
page in the Reference Manual. In a
nutshell, we have to define the type NT for the box coordinates
and the type ID for the id-number. Member functions
give access to the coordinates and the id-number. A static
member function returns the dimension.
// file: examples/Box_intersection_d/custom_box_grid.C
#include <CGAL/box_intersection_d.h>
#include <vector>
#include <algorithm>
#include <iterator>
#include <assert.h>
struct Box {
typedef int NT;
typedef std::ptrdiff_t ID;
int x[2], y[2];
Box( int x0, int x1, int y0, int y1) { x[0]=x0; x[1]=x1; y[0]=y0; y[1]=y1;}
static int dimension() { return 2; }
int min_coord(int dim) const { return x[dim]; }
int max_coord(int dim) const { return y[dim]; }
// id-function using address of current box,
// requires to work with pointers to boxes later
std::ptrdiff_t id() const { return (std::ptrdiff_t)(this); }
};
// 9 boxes of a grid
Box boxes[9] = { Box( 0,0,1,1), Box( 1,0,2,1), Box( 2,0,3,1), // low
Box( 0,1,1,2), Box( 1,1,2,2), Box( 2,1,3,2), // middle
Box( 0,2,1,3), Box( 1,2,2,3), Box( 2,2,3,3)};// upper
// 2 selected boxes as query; center and upper right
Box query[2] = { Box( 1,1,2,2), Box( 2,2,3,3)};
// With the special id-function we need to work on box pointers
Box* b_ptr[9] = { boxes, boxes+1, boxes+2, boxes+3, boxes+4, boxes+5,
boxes+6, boxes+7, boxes+8};
Box* q_ptr[2] = { query, query+1};
// callback function object writing results to an output iterator
template <class OutputIterator>
struct Report {
OutputIterator it;
Report( OutputIterator i) : it(i) {} // store iterator in object
// We write the position with respect to 'boxes' to the output iterator
// assuming that box b (the query box) is not interesting in the result.
void operator()( const Box* a, const Box*) {
*it++ = ( reinterpret_cast<Box*>(a->id()) - boxes);
}
};
template <class Iter> // helper function to create the function object
Report<Iter> report( Iter it) { return Report<Iter>(it); }
int main() {
// run the intersection algorithm and store results in a vector
std::vector<std::size_t> result;
CGAL::box_intersection_d( b_ptr, b_ptr+9, q_ptr, q_ptr+2,
report( std::back_inserter( result)),
std::ptrdiff_t(0));
// sort and check result
std::sort( result.begin(), result.end());
std::size_t chk[13] = {0,1,2,3,4,4,5,5,6,7,7,8,8};
assert( result.size()==13 && std::equal(chk,chk+13,result.begin()));
return 0;
}
Given a set of 3D points, we want to find all pairs of points that are less than a certain distance apart. We use the box intersection algorithm to find good candidates, namely those that are less than this specified distance apart in the norm, which is a good approximation of the Euclidean norm.
We use an unusual representation for the box, namely pointers to the 3D points themselves. We implement a special box traits class that interprets the point as a box of the dimensions epseps centered at this point. The value for eps is half the specified distance from above, i.e., points are reported if their distance is smaller than 2*eps.
The requirements for the box traits class are best studied on
page in the Reference Manual. In a
nutshell, we have to define the type NT for the box coordinates,
the type ID for the id-number, and the type Box_parameter
similar to the box handle, here Point_3* since we work with the pointers.
All member functions in the traits class are static. Two functions give
access to the max and min coordinates that we compute from the point
coordinates plus or minus the eps value, respectively. For the
id-number function the address of the point itself is
sufficient, since the points stay stable. Another function
returns the dimension.
The report callback function computes than the Euclidean distance and prints a message for points that are close enough.
Note that we need to reserve sufficient space in the points vector to avoid reallocations while we create the points vector and the boxes vector in parallel, since otherwise the points vector might reallocate and invalidate all pointers stored in the boxes so far.
// file: examples/Box_intersection_d/proximity_custom_box_traits.C
#include <CGAL/Simple_cartesian.h>
#include <CGAL/box_intersection_d.h>
#include <CGAL/point_generators_3.h>
#include <CGAL/copy_n.h>
#include <vector>
#include <algorithm>
#include <iterator>
typedef CGAL::Simple_cartesian<float> Kernel;
typedef Kernel::Point_3 Point_3;
typedef CGAL::Random_points_on_sphere_3<Point_3> Points_on_sphere;
std::vector<Point_3> points;
std::vector<Point_3*> boxes; // boxes are just pointers to points
const float eps = 0.1; // finds point pairs of distance < 2*eps
// Boxes are just pointers to 3d points. The traits class adds the
// +- eps size to each interval around the point, effectively building
// on the fly a box of size 2*eps centered at the point.
struct Traits {
typedef float NT;
typedef Point_3* Box_parameter;
typedef std::ptrdiff_t ID;
static int dimension() { return 3; }
static float coord( Box_parameter b, int d) {
return (d == 0) ? b->x() : ((d == 1) ? b->y() : b->z());
}
static float min_coord( Box_parameter b, int d) { return coord(b,d)-eps;}
static float max_coord( Box_parameter b, int d) { return coord(b,d)+eps;}
// id-function using address of current box,
// requires to work with pointers to boxes later
static std::ptrdiff_t id(Box_parameter b) { return (std::ptrdiff_t)(b); }
};
// callback function reports pairs in close proximity
void report( const Point_3* a, const Point_3* b) {
float dist = sqrt( CGAL::squared_distance( *a, *b));
if ( dist < 2*eps) {
std::cout << "Point " << (a - &(points.front())) << " and Point "
<< (b - &(points.front())) << " have distance " << dist
<< "." << std::endl;
}
}
int main() {
// create some random points on the sphere of radius 1.0
Points_on_sphere generator( 1.0);
points.reserve( 50);
for ( int i = 0; i != 50; ++i) {
points.push_back( *generator++);
boxes.push_back( & points.back());
}
// run the intersection algorithm and report proximity pairs
CGAL::box_self_intersection_d( boxes.begin(), boxes.end(),
report, Traits());
return 0;
}
Lutz Kettner and Andreas Meyer implemented the algorithms starting from the publication [ZE02]. We had access to the original C implementation of Afra Zomorodian, which helped clarifying some questions, and we are grateful to the help of Afra Zomorodian in answering our questions during his visit. We thank Steve Robbins for an excellent review for this package. Steve Robbins provided an independent and earlier implementation of this algorithm, however, we learned too late about this implementation.