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, threedimensional 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 axisaligned 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 isooriented 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 selfintersections) 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 $$ddimensional isooriented box is defined as the Cartesian product of $$d intervals. We call the box halfopen if the $$d intervals $${ [lo_{i},hi_{i})  0 i < d} are halfopen intervals, and we call the box closed if the $$d intervals $${ [lo_{i},hi_{i}]  0 i < d} are closed intervals. Note that closed boxes support zerowidth boxes and they can intersect at their boundaries, while nonempty halfopen boxes always have a positive volume and they only intersect iff their interiors overlap. The distinction between closed and halfopen 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 idnumber. 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 idnumber 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 selfintersection 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 halfopen 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 isooriented 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 idnumber.
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 twodimensional boxes in a grid layout of $$3 × 3 boxes. Additionally we pick the center box and the box in the upperright corner as our second box sequence query.
The default policy of the box type implements the idnumber 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 idnumber 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 upperright 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 axisaligned 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 $$[1,+1)^{3}. 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 idnumber 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 onetoone mapping between boxes and approximated geometry, which is the case here. It saves us the extra space that was needed for the explicit idnumber 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 halfopen boxes, see also Section . Closed boxes support zerowidth boxes and they can intersect at their boundaries, while halfopen boxes always have a positive volume and they only intersect iff their interiors overlap. The choice between closed or halfopen boxes is selected with the topology parameter and its two values:
The example program uses a twodimensional box with int coordinates and idnumbers that are by default explicitly stored. We create the same boxes as in the minimal example in Section . We create a $$3 × 3 grid of boxes, and two boxes for the query sequence, namely the box at the center and the box from the upperright corner of the grid.
We write a more involved callback function object Report that stores an output iterator and writes the idnumber 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 halfopen 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 upperright box in query intersects the four boxes of the upperright quadrant in boxes. Almost all intersections are with the box boundaries, thus, for the halfopen 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 idnumber 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 halfopen 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 segmenttree algorithm to the twowayscan algorithm, see [ZE02] for the details.
The streamed segmenttree algorithm needs $$O(n log$$^{d} (n) + k) worstcase running time and $$O(n) space, where $$n is the number of boxes in both input sequences, $$d the (constant) dimension of the boxes, and $$k the output complexity, i.e., the number of pairwise intersections of the boxes. The twowayscan algorithm needs $$O(n log$$(n) + l) worstcase running time and $$O(n) space, where $$l is the number of pairwise overlapping intervals in one dimensions (the dimension where the algorithm is used instead of the segment tree). Note that $$l is not necessarily related to $$k and using the twowayscan 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 segmenttree 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$$=sqrt(n) 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 threedimensional 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 nearoptimal 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 idnumber and use the address of the box itself as idnumber. 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 selfwritten 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 halfopen 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 idnumber. Member functions give access to the coordinates and the idnumber. 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]; } // idfunction 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 idfunction 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 $$L_{ } 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 $$[eps$$,+eps$$]^{3} 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 idnumber, 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 idnumber 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;} // idfunction 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.