Menelaos Karavelas
This chapter describes the two-dimensional segment Delaunay graph package of CGAL. We start with a few definitions in Section 26.1. The software design of the 2D segment Delaunay graph package is described in Section 26.2. In Section 26.3 we discuss the geometric traits of the 2D segment Delaunay graph package and in Section 26.4 the segment Delaunay graph hierarchy, a data structure suitable for fast nearest neighbor queries, is briefly described.
The 2D segment Delaunay graph package of CGAL is designed to compute the Delaunay graph of a set of possibly intersecting segments on the plane. Although we compute the Delaunay graph, we will often refer to its dual, the segment Voronoi diagram, since it is easier to explain and understand. The algorithm that has been implemented is incremental. The corresponding CGAL class is called Segment_Delaunay_graph_2<SegmentDelaunayGraphTraits_2,SegmentDelaunayGraphStructure_2> and will be discussed in more detail in the sequel. The interested reader may want to refer to the paper by Karavelas [Kar04] for the general idea as well as the details of the algorithm implemented.
Definitions. Before describing the details of the implementation we make a brief introduction to the theory of segment Delaunay graphs and segment Voronoi diagrams. The segment Voronoi diagram is defined over a set of non-intersecting sites, which can either be points or linear segments, which we assume that are given through their endpoints. The segment Voronoi diagram a subdivision of the plane into connected regions, called cells, associated with the sites. The dual graph of the segment Voronoi diagram is called the segment Delaunay graph. The cell of a site is the locus of points on the plane that are closer to than any other site , . The distance of a point in the plane to a site is defined as the minimum of the Euclidean distances of from the points in . Hence, if is a point , then
whereas if , is a segment, then
min
where denotes the Euclidean norm. It can easily be seen that it is a generalization of the Voronoi diagram for points.
In many applications the restriction that sites are non-intersecting is too strict. Often we want to allow segments that touch at their endpoints, or even segments that overlap or intersect properly at their interior (for example, see Fig. 26.1). Allowing such configurations poses certain problems. More specifically, when we allow segments to touch at their endpoints we may end up with pairs of segments whose bisector is two-dimensional. If we allow pairs of segments that intersect properly at their interior, the interiors of their Voronoi cells are no longer simply connected. In both cases above the resulting Voronoi diagrams are no longer instances of abstract Voronoi diagrams (cf. [Kle89]), which has a direct consequence on the efficient computation of the corresponding Voronoi diagram. The remedy to these problems is to consider linear segments not as one object, but rather as three, namely the two endpoints and the interior. This choice guarantees that all bisectors in the Voronoi diagram are one-dimensional and that all Voronoi cells are simply connected. Moreover, we further distinguish between two cases, according to the type of intersecting pair that our input data set contains. A pair of sites is called weakly intersecting if they a single common point and this common point does not lie in the interior of any of the two sites. A pair of sites is called strongly intersecting if they intersect and they either have more than one common point or their common point lies in the interior of at least one of the two sites. As it will be seen later the two cases have different representation (and thus storage) requirements, as well as they require a somehow different treatment on how the predicates are evaluated. Having made the distinction between weakly and strongly intersecting sites, and having said that segment sites are treated as three objects, we are now ready to precisely define the Delaunay graph we compute. Given a set of input sites, let be the set of points and (open) segments in the arrangement of . The 2D segment Delaunay graph package of CGAL computes the (triangulated) Delaunay graph that is dual to the Euclidean Voronoi diagram of the sites in the set .
The segment Delaunay graph is uniquely defined once we have the segment Voronoi diagram. If the all sites are in general position, then Delaunay graph is a graph with triangular faces away from the convex hull of the set of sites. To unify our approach and handling of the Delaunay graph we add to the set of (finite) sites a fictitious site at infinity, which we call the site at infinity. We can then connect all vertices of the outer face of the Delaunay graph to the site at infinity which gives us a graph with the property that all of its faces are now triangular. However, the Delaunay graph is not a triangulation for two main reasons: we cannot always embed it on the plane with straight line segments that yield a triangulation and, moreover, we may have two faces of the graph that have two edges in common, which is not allowed in a triangulation.
We would like to finish our brief introduction to the theory of segment Delaunay graphs and segment Voronoi diagrams by discussing the concept of general position. We say that a set of sites is in general position if no two triplets of sites have the same tritangent Voronoi circle. This statement is rather technical and it is best understood in the context of points. The equivalent statement for points is that we have no two triplets of points that define the same circumcircle, or equivalently that no four points are co-circular. The statement about general position made above is a direct generalization of the (much simpler to understand) statement about points. On the contrary, when we have sites in degenerate position, the Delaunay graph has faces with more than three edges on their boundary. We can get a triangulated version of the Delaunay graph by simply triangulating the corresponding faces in an arbitrary way. In fact the algorithm that has been implemented in CGAL has the property that it always returns a valid triangulated version of the segment Delaunay graph. By valid we mean that it contains the actual (non-triangulated) Delaunay graph, and whenever there are faces with more than three faces then they are triangulated. The way that they are triangulated depends on the order of insertion of the sites in the diagram.
One final remark has to be made with respect to the difference between the set of input sites and the set of output sites. The set of input sites consists of the closed sites that the user inserts in the diagram. Since segment sites are treated as three objects, internally our algorithm sees only points and open segments. As a result, from the point of view of the algorithm, the input sites have no real meaning. What has real meaning is the set of sites that correspond to cells of the Voronoi diagram and this is the set of output sites.
Degenerate Dimensions. The dimension of the segment Delaunay graph is in general 2. The exceptions to this rule are as follows:
The 2D segment Delaunay graph class Segment_Delaunay_graph_2<SegmentDelaunayGraphTraits_2,SegmentDelaunayGraphDataStructure_2> follows the design of the triangulation package of CGAL. It is parametrized by two arguments:
Strongly Intersecting Sites and their Representation. As we have mentioned above, the segment Delaunay graph package of CGAL is designed to support the computation of the segment Voronoi diagram even when the input segment sites are intersecting. This choice poses certain issues for the design of the software package. The major concern is the representation of the subsegments that appear in the arrangement of the these sites, because the sites in the arrangement are the ones over which the diagram is actually defined. A direct consequence of the choice of representation is the algebraic degree of the predicates involved in the computation of the segment Delaunay graph, as well as the storage requirements for the subsegments and points on intersection in the arrangement.
The case of weakly intersecting sites does not require any special treatment. We can simply represent points by their coordinates and segments by their endpoints. In the case of strongly intersecting sites, the obvious choice to use the afore-mentioned representation has severe disadvantages. Consider two strongly intersecting segments and , whose endpoints have homogeneous coordinates of size . Their intersection point will have homogeneous coordinates of bit size . This effect can be cascaded, which implies that after inserting (input) segments we can arrive at having points of intersection whose bit sizes are exponential with respect to , i.e., their homogeneous coordinates will have bit size . Not only the points of intersection, but also the adjacent subsegments will be represented by quantities of arbitrarily high bit size, and as a result we would not be able to give a bound on the bit sizes of the coordinates of the points of intersection. As a result, we would not be able to give a bound on the memory needed to store these coordinates. An equally important consequence is that we would also not be able to give a bound on the algebraic degree of the algebraic expressions involved in the evaluation of the predicates.
Such a behavior is obviously undesirable. For robustness, efficiency, and scalability purposes, it is critical that the bit size of the algebraic expressions in the predicates does not depend on the input size. For this reason, as well as for others to be discussed below, we decided to represent sites in a implicit manner, which somehow encodes the history of their construction. In particular, we exploit the fact that points of intersection always lie on two input segments, and that segments that are not part of the input are always supported by input segments.
For example, let us consider the configuration in Fig. 26.2.0.0.1. We assume that the segments , , are inserted in that order. Upon the insertion of , our algorithm will split the segment into the subsegments and , then add , and finally insert the subsegments and . How do we represent the five new sites? will be represented by its two defining segments and . The segment will be represented by two segments, a point, and a boolean. The first segment is , which is always the segment with the same support as the newly created segment. The second segment is and the point is . The boolean indicates whether the first endpoint of is an input point; in this case the boolean is equal to true. The segment will also be represented by two segments, a point, and a boolean, namely, (the supporting segment of ), and false (it is the second endpoint of that is an input point). Subsegments and are represented analogously. Consider now what happens when we insert . The point will again be represented by two segments, but not and . In fact, it will be represented by (the supporting segment of ) and . will be represented by two segments, a point, and a boolean (, and false), and similarly for and . On the other hand, both endpoints of are non-input points. In such a case we represent the segment by three input segments. More precisely, is represented by the segments (the supporting segment of ), (it defines along with ) and (it defines along with ).
The five different presentations, two for points (coordinates; two input segments) and three for segments (two input points; two input segments, an input point and a boolean; three input segments), form a closed set of representations and thus represent any point of intersection or subsegment regardless of the number of input segments. Moreover, every point (input or intersection) has homogeneous coordinates of bit size at most . The supporting lines of the segments (they are needed in some of the predicates) have coefficients which are always of bit size . As a result, the bit size of the expressions involved in our predicates will always be , independently of the size of the input. The SegmentDelaunayGraphSite_2 concept encapsulates the ideas presented above. A site is represented in this concept by up to four points and a boolean, or up to six points, depending on its type. The class Segment_Delaunay_graph_site_2<K> implements this concept.
Even this representation, however, has some degree of redundancy. The endpoint of a segment appears in both the representation of the (open) segment site as well as the representation of the point site itself. The situation becomes even worse in the presence of strongly intersecting sites: a point may appear in the representation of multiple subsegments and/or points of intersection. To avoid this redundancy, input points are stored in a container, and the various types of sites (input points and segments, points of intersection, subsegments with one or two points of intersection as endpoints) only store handles to the points in the container. This is achieved by the Segment_Delaunay_graph_storage_site_2<Gt> class which is a model of the corresponding concept: SegmentDelaunayGraphStorageSite_2. This concept enforces a site to be represented by up to 6 handles (which are very lightweight objects) instead of 6 points, which are, compared to handles of course, very heavy objects.
Optimizing Memory Allocation. There are applications where we know beforehand that the input consists of only weakly intersecting sites. In these cases the site representation described above poses a significant overhead in the memory requirements of our implementation: instead of representing sites with up to two points (or ultimately with to two handles), we require sites to store six points (respectively, six handles). To avoid this overhead we have introduced two series of traits classes:
The predicates required for the computation of the segment Voronoi diagram are rather complicated. It is not the purpose of this document to discuss them in detail. The interested reader may refer to Burnikel's thesis [Bur96], where it is shown that in the case of weakly intersecting sites represented in homogeneous coordinates of bit size , the maximum bit size of the algebraic expressions involved in the predicates is . Given our site representation given above we can guarantee that even in the case of strongly intersecting sites, the algebraic degree of the predicates remains , independently of the size of the input. What we want to focus in the remainder of this section are the different kinds of filtering techniques that we have employed in our implementation.
Geometric Filtering. Our representation of sites is coupled very naturally, with what we call geometric filtering. The technique amounts to performing simple geometric tests exploiting the representation of our data, as well as the geometric structure inherent in our problem, in order to evaluate predicates in seemingly degenerate configurations. Geometric filtering can be seen as a preprocessing step before performing arithmetic filtering. Roughly speaking, by arithmetic filtering we mean that we first try to evaluate the predicates using a fixed-precision floating-point number type (such as double), and at the same time keep error bounds on the numerical errors of the computations we perform. If the numerical errors are too big and do not permit us to evaluate the predicate, we switch to an exact number type, and repeat the evaluation of the predicate. Geometric filtering can help by eliminating situations in which the arithmetic filter will fail, thus decreasing the number of times we need to evaluate a predicate using exact arithmetic.
To illustrate the application and effectiveness of this approach, let us consider a very simple example usage. Suppose we want to determine if two non-input points are identical (we assume here that the input sites are represented by doubles). In order to do that we need to compute their coordinates and compare them. If the two points are identical, the answer to our question using double arithmetic may be wrong (due to numerical errors), in which case we will have to reside to the more expensive exact computation. Instead, before testing the coordinates for equality, we can use the representation of the points to potentially answer the question. More specifically, and this is the geometric filtering part of the computation, we can first test if the defining segments of the two points are the same. If they are not, then we proceed to comparing their coordinates as usual. Testing the defining segments for equality does not involve any arithmetic operations on the input, but rather only comparisons on doubles. By performing this very simple test we avoid a numerically difficult computation, which could be performed thousands of times during the computation of a Delaunay graph.
Geometric filtering has been implemented in all our models of the SegmentDelaunayGraphTraits_2 concept. These models are the classes: Segment_Delaunay_graph_traits_2<K,MTag>, Segment_Delaunay_graph_traits_without_intersections_2<K,MTag>, Segment_Delaunay_graph_filtered_traits_2<CK,CM,EK,EM,FK,FM> and Segment_Delaunay_graph_filtered_traits_without_intersections_2<CK,CM,EK,EM,FK,FM>.
Arithmetic Filtering. As mentioned above, performing computations with exact arithmetic can be very costly. For this reason we have devoted considerable effort in implementing different kinds of arithmetic filtering mechanisms. Presently, there two ways of performing arithmetic filtering for the predicates involved in the computation of segment Delaunay graphs:
Let's consider once more the classes Segment_Delaunay_graph_2<K,MTag> and Segment_Delaunay_graph_with_intersections_2<K,MTag>. The template parameter MTag provides another degree of freedom to the user, who can indicate the type of arithmetic operations to be used in the evaluation of the predicates. More specifically, in both classes, MTag can be CGAL::Sqrt_field_tag, in which case the predicates will be evaluated using all four basic arithmetic operations plus square roots; this requires, of course, that the number type used in the kernel K supports these operations exactly. The second choices are CGAL::Field_tag for the Segment_Delaunay_graph_2<K,MTag> class, and CGAL::Ring_tag for the Segment_Delaunay_graph_with_intersections_2<K,MTag> class. In the first case we indicate that we want the predicates to be computed using only the four basic arithmetic operations, whereas in the second case we evaluate the predicates using only ring operations. Again, for the predicates to be evaluated correctly, the number type used in the kernel K must support the corresponding operations exactly.
The semantics for the template parameters CM, FM and EM in the Segment_Delaunay_graph_filtered_traits_2<CK,CM,EK,EM,FK,FM> and Segment_Delaunay_graph_filtered_traits_without_intersections_2<CK,CM,EK,EM,FK,FM> classes are analogous. With each of these template parameters we can control the type of arithmetic operations that are going to be used in calculations involving each of the corresponding kernels CK, FK and EK. When the Segment_Delaunay_graph_filtered_traits_2<CK,CM,EK,EM,FK,FM> is used the possible values for CM, FM and EM are CGAL::Sqrt_field_tag and CGAL::Field_tag, whereas if the Segment_Delaunay_graph_filtered_traits_without_intersections_2<CK,CM,EK,EM,FK,FM> class is used, the possible values are CGAL::Sqrt_field_tag and CGAL::Ring_tag. The semantics are the same as in the case of the Segment_Delaunay_graph_2<K,MTag> and Segment_Delaunay_graph_with_intersections_2<K,MTag> classes.
The Segment_Delaunay_graph_hierarchy_2<SegmentDelaunayGraphTraits_2, SSTag, SegmentDelaunayGraphDataStructure_2> class is the analogue of the Triangulation_hierarchy_2 or the Apollonius_graph_hierarchy_2 classes, applied to the segment Delaunay graph. It consists of a hierarchy of segment Delaunay graphs constructed in a manner analogous to the Delaunay hierarchy by Devillers [Dev02]. Unlike the triangulation hierarchy or the Apollonius graph hierarchy, the situation here is more complicated because of two factors: firstly, segments are treated as three objects instead of one (the two endpoints and the interior of the segments), and secondly, the presence of strongly intersecting sites complicates significantly the way the hierarchy is constructed. The interested reader may refer to the paper by Karavelas [Kar04] for the details of the construction of the hierarchy. Another alternative is to have a hybrid hierarchy that consists of the segment Delaunay graph at the bottom-most level and point Voronoi diagrams at all other levels. This choice seems to work very well in practice , primarily because it avoids the overhead of maintaining a Delaunay graph for segments at the upper levels of the hierarchy. However, it seems much less likely to be possible to give any theoretical guarantees for its performance, in contrast to the hierarchy with segment Delaunay graphs at all levels (cf. [Kar04]). The user can choose between the two types of hierarchies by means of the template parameter SSTag. If SSTag is set to false (which is also the default value), the upper levels of the hierarchy consist of point Delaunay graphs. If SSTag is set to true, we have segment Delaunay graphs at all levels of the hierarchy.
The class Segment_Delaunay_graph_hierarchy_2<SegmentDelaunayGraphTraits_2, SSTag, SegmentDelaunayGraphDataStructure_2> has exactly the same interface and functionality as the Segment_Delaunay_graph_2<SegmentDelaunayGraphTraits_2,SegmentDelaunayGraphDataStructure_2> class. Using the segment Delaunay graph hierarchy involves an additional cost in space and time for maintaining the hierarchy. Our experiments have shown that it usually pays off to use the hierarchy for inputs consisting of more than about 1,000 sites.
The following example shows to use the segment Delaunay graph traits in conjunction with the Filtered_exact<CT,ET> mechanism. In addition it shows how to use a few of the iterators provided by the Segment_Delaunay_graph_2 class in order to count a few site-related quantities.
// file: sdg-count-sites.C
#include <CGAL/basic.h>
// standard includes
#include <iostream>
#include <fstream>
#include <cassert>
// define the exact number type
# include <CGAL/Quotient.h>
# include <CGAL/MP_Float.h>
typedef CGAL::Quotient<CGAL::MP_Float> ENT;
// define the kernels
#include <CGAL/Simple_cartesian.h>
typedef CGAL::Simple_cartesian<double> CK;
typedef CGAL::Simple_cartesian<ENT> EK;
// typedefs for the traits and the algorithm
#include <CGAL/Segment_Delaunay_graph_filtered_traits_2.h>
#include <CGAL/Segment_Delaunay_graph_2.h>
typedef CGAL::Segment_Delaunay_graph_filtered_traits_2<CK,
/* The construction kernel allows for / and sqrt */ CGAL::Sqrt_field_tag,
EK,
/* The exact kernel supports field ops exactly */ CGAL::Field_tag> Gt;
typedef CGAL::Segment_Delaunay_graph_2<Gt> SDG2;
using namespace std;
int main() {
ifstream ifs("data/sitesx.cin");
assert( ifs );
SDG2 sdg;
SDG2::Site_2 site;
while ( ifs >> site ) { sdg.insert( site ); }
ifs.close();
assert( sdg.is_valid(true, 1) );
cout << endl << endl;
// print the number of input and output sites
cout << "# of input sites : " << sdg.number_of_input_sites() << endl;
cout << "# of output sites: " << sdg.number_of_output_sites() << endl;
unsigned int n_ipt(0), n_iseg(0), n_opt(0), n_oseg(0), n_ptx(0);
// count the number of input points and input segments
SDG2::Input_sites_iterator iit;
for (iit = sdg.input_sites_begin(); iit != sdg.input_sites_end(); ++iit)
{
if ( iit->is_point() ) { n_ipt++; } else { n_iseg++; }
}
// count the number of output points and output segments, as well
// as the number of points that are points of intersection of pairs
// of strongly intersecting sites
SDG2::Output_sites_iterator oit;
for (oit = sdg.output_sites_begin(); oit != sdg.output_sites_end(); ++oit)
{
if ( oit->is_segment() ) { n_oseg++; } else {
n_opt++;
if ( !oit->is_input() ) { n_ptx++; }
}
}
cout << endl << "# of input segments: " << n_iseg << endl;
cout << "# of input points: " << n_ipt << endl << endl;
cout << "# of output segments: " << n_oseg << endl;
cout << "# of output points: " << n_opt << endl << endl;
cout << "# of intersection points: " << n_ptx << endl;
return 0;
}
The following example shows how to use the segment Delaunay graph hierarchy along with the filtered traits class that supports intersecting sites.
// file: sdg-filtered-traits.C
#include <CGAL/basic.h>
// standard includes
#include <iostream>
#include <fstream>
#include <cassert>
// example that uses the filtered traits and
// the segment Delaunay graph hierarchy
// choose the kernel
#include <CGAL/Simple_cartesian.h>
struct Rep : public CGAL::Simple_cartesian<double> {};
// typedefs for the traits and the algorithm
#include <CGAL/Segment_Delaunay_graph_hierarchy_2.h>
#include <CGAL/Segment_Delaunay_graph_filtered_traits_2.h>
struct Gt
: public CGAL::Segment_Delaunay_graph_filtered_traits_2<Rep> {};
typedef CGAL::Segment_Delaunay_graph_hierarchy_2<Gt> SDG2;
int main()
{
std::ifstream ifs("data/sites.cin");
assert( ifs );
SDG2 sdg;
SDG2::Site_2 site;
// read the sites and insert them in the segment Delaunay graph
while ( ifs >> site ) {
sdg.insert(site);
}
// validate the segment Delaunay graph
assert( sdg.is_valid(true, 1) );
return 0;
}
The following example demonstrates how to recover the defining sites for the edges of the Voronoi diagram (which are the duals of the edges of the segment Delaunay graph computed).
// file: sdg-voronoi-vertices.C
#include <CGAL/basic.h>
// standard includes
#include <iostream>
#include <fstream>
#include <cassert>
#include <string>
// define the kernel
#include <CGAL/Simple_cartesian.h>
#include <CGAL/Filtered_kernel.h>
typedef CGAL::Simple_cartesian<double> CK;
typedef CGAL::Filtered_kernel<CK> Kernel;
// typedefs for the traits and the algorithm
#include <CGAL/Segment_Delaunay_graph_traits_2.h>
#include <CGAL/Segment_Delaunay_graph_2.h>
typedef CGAL::Segment_Delaunay_graph_traits_2<Kernel> Gt;
typedef CGAL::Segment_Delaunay_graph_2<Gt> SDG2;
using namespace std;
int main()
{
ifstream ifs("data/sites2.cin");
assert( ifs );
SDG2 sdg;
SDG2::Site_2 site;
// read the sites from the stream and insert them in the diagram
while ( ifs >> site ) { sdg.insert( site ); }
ifs.close();
// validate the diagram
assert( sdg.is_valid(true, 1) );
cout << endl << endl;
/*
// now walk through the non-infinite edges of the segment Delaunay
// graphs (which are dual to the edges in the Voronoi diagram) and
// print the sites defining each Voronoi edge.
//
// Each oriented Voronoi edge (horizontal segment in the figure
// below) is defined by four sites A, B, C and D.
//
// \ /
// \ B /
// \ /
// C ----------------- D
// / \
// / A \
// / \
//
// The sites A and B define the (oriented) bisector on which the
// edge lies whereas the sites C and D, along with A and B define
// the two endpoints of the edge. These endpoints are the Voronoi
// vertices of the triples A, B, C and B, A, D.
// If one of these vertices is the vertex at infinity the string
// "infinite vertex" is printed; the corresponding Voronoi edge is
// actually a stright-line or parabolic ray.
// The sites below are printed in the order A, B, C, D.
*/
string inf_vertex("infinite vertex");
char vid[] = {'A', 'B', 'C', 'D'};
SDG2::Finite_edges_iterator eit = sdg.finite_edges_begin();
for (int k = 1; eit != sdg.finite_edges_end(); ++eit, ++k) {
SDG2::Edge e = *eit;
// get the vertices defining the Voronoi edge
SDG2::Vertex_handle v[] = { e.first->vertex( sdg.ccw(e.second) ),
e.first->vertex( sdg.cw(e.second) ),
e.first->vertex( e.second ),
sdg.tds().mirror_vertex(e.first, e.second) };
cout << "--- Edge " << k << " ---" << endl;
for (int i = 0; i < 4; i++) {
// check if the vertex is the vertex at infinity; if yes, print
// the corresponding string, otherwise print the site
if ( sdg.is_infinite(v[i]) ) {
cout << vid[i] << ": " << inf_vertex << endl;
} else {
cout << vid[i] << ": " << v[i]->site() << endl;
}
}
cout << endl;
}
return 0;
}