Geometry Algorithms ( geo_alg )

All functions listed in this section work for geometric objects based on both floating-point and exact (rational) arithmetic. In particular, point can be replace by rat_point, segment by rat_segment, and circle by rat_circle.

The floating point versions are faster but unreliable. They may produce incorrect results, abort, or run forever. Only the rational versions will produce correct results for all inputs.

The include-file for the rational version is rat_geo_alg.h, the include-file for the floating point version is float_geo_alg.h, and geo_alg.h includes both versions. Including both versions increases compile time. An alternative name for geo_alg.h is plane_alg.h.

$\bullet$ Convex Hulls

list<point>	CONVEX_HULL(list<point> L)
		CONVEX_HULL takes as argument a list of points and returns the polygon representing the convex hull of L. The cyclic order of the vertices in the result list corresponds to counter-clockwise order of the vertices on the hull. The algorithm calls our current favorite of the algorithms below.
polygon	CONVEX_HULL_POLY(list<point> L)
		as above, but returns the convex hull of L as a polygon.
list<point>	UPPER_CONVEX_HULL(list<point> L)
		returns the upper convex hull of L.
list<point>	LOWER_CONVEX_HULL(list<point> L)
		returns the lower convex hull of L.
list<point>	CONVEX_HULL_S(list<point> L)
		as above, but the algorithm is based on the sweep paradigm. Running time is O(nlog n) in the worst and in the best case.
list<point>	CONVEX_HULL_IC(list<point> L)
		as above, but the algorithm is based on incremental construction. The running time is O(n2) worst case and is O(nlog n) expected case. The expectation is computed as the average over all permutations of L. The running time is linear in the best case.
list<point>	CONVEX_HULL_RIC(list<point> L)
		as above. The algorithm permutes L randomly and then calls the preceding function.
double	WIDTH(list<point> L, line& l1, line& l2)
		returns the square of the minimum width of a stripe covering all points in Land the two boundaries of the stripe. Precondition L is non-empty

$\bullet$ Halfplane intersections

void	HALFPLANE_INTERSECTION(list<line> L, list<line>& Lout)
		For every line $\ell$ $\in$ L let h$\ell$ be the closed halfplane lying on the positive side of $\ell$, i.e., h$\ell$ = { p $\in$ $\hbox$R2 | orientation($\ell$, p) > = 0 }, and let H = $\cap_{\ell\in L}^{}$ h$\ell$. Then HALFPLANE_INTERSECTION computes the list of lines Lout defining the boundary of H in counter-clockwise ordering.

$\bullet$ Triangulations

edge	TRIANGULATE_POINTS(list<point> L, GRAPH<point,int>& T)
		computes a triangulation (planar map) T of the points in L and returns an edge of the outer face (convex hull).
void	DELAUNAY_TRIANG(list<point> L, GRAPH<point,int>& DT)
		computes the delaunay triangulation DT of the points in L.
void	DELAUNAY_DIAGRAM(list<point> L, GRAPH<point,int>& DD)
		computes the delaunay diagram DD of the points in L.
void	F_DELAUNAY_TRIANG(list<point> L, GRAPH<point,int>& FDT)
		computes the furthest point delaunay triangulation FDT of the points in L.
void	F_DELAUNAY_DIAGRAM(list<point> L, GRAPH<point,int>& FDD)
		computes the furthest point delaunay diagram FDD of the points in L.

$\bullet$ Constraint Triangulations

edge	TRIANGULATE_SEGMENTS(list<segment> L, GRAPH<point,segment>& G)
		computes a constrained triangulation (planar map) T of the segments in L (trivial segments representing points are allowed). The function returns an edge of the outer face (convex hull).
edge	DELAUNAY_TRIANG(list<segment> L, GRAPH<point,segment>& G)
		computes a constrained Delaunay triangulation T of the segments in L. The function returns an edge of the outer face (convex hull).
edge	TRIANGULATE_PLANE_MAP(GRAPH<point,segment>& G)
		computes a constrained triangulation T of the plane map (counter-clockwise straight-line embedded Graph) G. The function returns an edge of the outer face (convex hull). Precondition G is simple.
edge	DELAUNAY_TRIANG(GRAPH<point,segment>& G)
		computes a constrained Delaunay triangulation T of the plane map G. The function returns an edge of the outer face (convex hull). Precondition G is simple.
edge	TRIANGULATE_POLYGON(polygon P, GRAPH<point,segment>& G, list<edge>& inner_edges, list<edge>& outer_edges, list<edge>& boundary_edges)
		triangulates the interior and exterior of the simple polygon P and stores all edges of the inner (outer) triangulation in inner_edges (outer_edges) and the edges of the polygon boundary in boundary_edges. The function returns an edge of the convex hull of P if P is simple and nil otherwise.
edge	TRIANGULATE_POLYGON(gen_polygon GP, GRAPH<point,segment>& G, list<edge>& inner_edges, list<edge>& outer_edges, list<edge>& boundary_edges, list<edge>& hole_edges)
		triangulates the interior and exterior of the generalized polygon GP and stores all edges of the inner (outer) triangulation in inner_edges (outer_edges). The function returns nil if GP is trivial, and an edge of the convex hull otherwise. boundary_edges contains the edges of every counter-clockwise oriented boundary cycle of GP, and hole_edges contains the edges on every clockwise oriented boundary cycle of GP. Note that the reversals of boundary and hole edges will be returned in inner_edges. Precondition GP is simple.
edge	CONVEX_COMPONENTS(polygon P, GRAPH<point,segment>& G, list<edge>& inner_edges, list<edge>& boundary)
		if P is a bounded and non-trivial simple polygon its interior is decomposed into convex parts. All inner edges of the constructed decomposition are returned in inner_edges. boundary_edges contains the edges of the polygon boundary Note that the reversals of boundary edges will be stored in inner_edges. The function returns an edge of the convex hull if P is simple and non-trivial and nil otherwise.
edge	CONVEX_COMPONENTS(gen_polygon GP, GRAPH<point,segment>& G, list<edge>& inner_edges, list<edge>& boundary_edges, list<edge>& hole_edges)
		if GP is a bounded and non-trivial generalized polygon, its interior is decomposed into convex parts. All inner edges of the constructed decomposition are returned in inner_edges. boundary_edges contains the edges of every counter-clockwise oriented boundary cycle of GP, and hole_edges contains the edges of every clockwise oriented boundary cycle of GP. Note that the reversals of boundary and hole edges will be stored in inner_edges. The function returns an edge of the convex hull if GP is a bounded and non-trivial and nil otherwise. Precondition GP must be simple.

$\bullet$ Minkowski Sums

Please note that the Minkowski sums only work reliable for the rational kernel.

gen_polygon	MINKOWSKI_SUM(polygon P, polygon R)
		computes the Minkowski sum of P and R.
gen_polygon	MINKOWSKI_DIFF(polygon P, polygon R)
		computes the Minkowski difference of P and R, i.e. the Minkowski sum of P and R.reflect(point(0,0)).
gen_polygon	MINKOWSKI_SUM(gen_polygon P, polygon R)
		computes the Minkowski sum of P and R.
gen_polygon	MINKOWSKI_DIFF(gen_polygon P, polygon R)
		computes the Minkowski difference of P and R, i.e. the Minkowski sum of P and R.reflect(point(0,0)).

The following variants of the MINKOWSKI functions take two additional call-back function arguments conv_partition and conv_unite which are used by the algorithm to partition the input polygons into convex parts and for computing the union of a list of convex polygons, respectively (instead of using the default methods).

gen_polygon	MINKOWSKI_SUM(polygon P, polygon R, void (*conv_partition)(gen_polygon p, polygon r, list<polygon>& lp, list<polygon>& lr), gen_polygon (*conv_unite)(list<gen_polygon> ))
gen_polygon	MINKOWSKI_DIFF(polygon P, polygon R, void (*conv_partition)(gen_polygon p, polygon r, list<polygon>& lp, list<polygon>& lr), gen_polygon (*conv_unite)(list<gen_polygon> ))
gen_polygon	MINKOWSKI_SUM(gen_polygon P, polygon R, void (*conv_partition)(gen_polygon p, polygon r, list<polygon>& lp, list<polygon>& lr), gen_polygon (*conv_unite)(list<gen_polygon> ))
gen_polygon	MINKOWSKI_DIFF(gen_polygon P, polygon R, void (*conv_partition)(gen_polygon p, polygon r, list<polygon>& lp, list<polygon>& lr), gen_polygon (*conv_unite)(list<gen_polygon> ))

$\bullet$ Euclidean Spanning Trees

void	MIN_SPANNING_TREE(list<point> L, GRAPH<point,int>& T)
		computes the Euclidian minimum spanning tree T of the points in L.

$\bullet$ Triangulation Checker

bool	Is_Convex_Subdivision(GRAPH<point,int> G)
		returns true if G is a convex planar subdivision.
bool	Is_Triangulation(GRAPH<point,int> G)
		returns true if G is convex planar subdivision in which every bounded face is a triangle or if all nodes of G lie on a common line.
bool	Is_Delaunay_Triangulation(GRAPH<point,int> G, delaunay_voronoi_kind kind)
		checks whether G is a nearest (kind = NEAREST) or furthest (kind = FURTHEST) site Delaunay triangulation of its vertex set. G is a Delaunay triangulation iff it is a triangulation and all triangles have the Delaunay property. A triangle has the Delaunay property if no vertex of an adjacent triangle is contained in the interior (kind = NEAREST) or exterior (kind = FURTHEST) of the triangle.
bool	Is_Delaunay_Diagram(GRAPH<point,int> G, delaunay_voronoi_kind kind)
		checks whether G is a nearest (kind = NEAREST) or furthest (kind = FURTHEST) site Delaunay diagram of its vertex set. G is a Delaunay diagram if it is a convex subdivision, if the vertices of any bounded face are co-circular, and if every triangulation of G is a Delaunay triangulation.

$\bullet$ Voronoi Diagrams

void	VORONOI(list<point> L, GRAPH<circle,point>& VD)
		VORONOI takes as input a list of points (sites) L. It computes a directed graph VD representing the planar subdivision defined by the Voronoi diagram of L. For each node v of VD G[v] is the corresponding Voronoi vertex (point) and for each edge e G[e] is the site (point) whose Voronoi region is bounded by e. The algorithm has running time O(n2) in the worst case and O(nlog n) with high probability, where n is the number of sites.
void	F_VORONOI(list<point> L, GRAPH<circle,point>& FVD)
		computes the farthest point Voronoi Diagram FVD of the points in L.
circle	LARGEST_EMPTY_CIRCLE(list<point> L)
		computes a largest circle whose center lies inside the convex hull of L that contains no point of L in its interior. Returns the trivial circle if L is empty.
circle	SMALLEST_ENCLOSING_CIRCLE(list<point> L)
		computes a smallest circle containing all points of L in its interior.
void	ALL_EMPTY_CIRCLES(list<point> L, list<circle>& CL)
		computes the list CL of all empty circles passing through three or more points of L.
void	ALL_ENCLOSING_CIRCLES(list<point> L, list<circle>& CL)
		computes the list CL of all enclosing circles passing through three or more points of L.

An annulus is either the region between two concentric circles or the region between two parallel lines.

bool	MIN_AREA_ANNULUS(list<point> L, point& center, point& ipoint, point& opoint, line& l1)
		computes the minimum area annulus containing the points of L. The function returns false if all points in L are collinear and returns true otherwise. In the former case a line passing through the points in L is returned in l1, and in the latter case the annulus is returned by its center and a point on the inner and the outer circle, respectively.
bool	MIN_WIDTH_ANNULUS(list<point> L, point& center, point& ipoint, point& opoint, line& l1, line& l2)
		computes the minimum width annulus containing the points of L. The function returns false if the minimum width annulus is a stripe and returns true otherwise. In the former case the boundaries of the stripes are returned in l1 and l2 and in the latter case the annulus is returned by its center and a point on the inner and the outer circle, respectively.
void	CRUST(list<point> L0, GRAPH<point,int>& G)
		takes a list L0 of points and traces to guess the curve(s) from which L0 are sampled. The algorithm is due to Amenta, Bern, and Eppstein []. The algorithm is guaranteed to succeed if L0 is a sufficiently dense sample from a smooth closed curve.
bool	Is_Voronoi_Diagram(GRAPH<circle,point> G, delaunay_voronoi_kind kind)
		checks whether G represents a nearest (kind = NEAREST) or furthest (kind = FURTHEST) site Voronoi diagram. Voronoi diagrams of point sites are represented as planar maps as follows: There is a vertex for each vertex of the Voronoi diagram and, in addition, a vertex ``at infinity'' for each ray of the Voronoi diagram. Vertices at infinity have degree one. The edges of the graph correspond to the edges of the Voronoi diagram. Figure Voronoi diagram shows some examples and the chapter on Voronoi diagrams of the LEDA-book [58] contains more details. Each edge is labeled with the site (class POINT) owning the region to its left and each vertex is labeled with a triple of points (= the three defining points of a CIRCLE). For a ``finite'' vertex the three points are any three sites associated with regions incident to the vertex (and hence the center of the circle is the position of the vertex in the plane) and for a vertex at infinity the three points are collinear and the first point and the third point of the triple are the sites whose regions are incident to the vertex at infinity. Let a and c be the first and third point of the triple respectively; a and c encode the geometric position of the vertex at infinity as follows: the vertex lies on the perpendicular bisector of a and c and to the left of the segment ac.

$\bullet$ Line Segment Intersection

void	SEGMENT_INTERSECTION(list<segment> S, GRAPH<point,segment>& G, bool embed=false)
		takes a list of segments S as input and computes the planar graph G induced by the set of straight line segments in S. The nodes of G are all endpoints and all proper intersection points of segments in S. The edges of G are the maximal relatively open subsegments of segments in S that contain no node of G. The edges are directed as the corresponding segments. If the flag embed is true, the corresponding planar map is computed. Note that for each edge eG[e] is the input segment that contains e (see the LEDA book for details)
void	SWEEP_SEGMENTS(list<segment> S, GRAPH<point,segment>& G, bool embed=false, bool use_optimization = true)
		as above. The algorithm ([10]) runs in time O((n + s)log n) + m), where n is the number of segments, s is the number of vertices of the graph G, and m is the number of edges of G. If S contains no overlapping segments then m = O(n + s). If embed is true the running time increases by O(mlog m). If use_optimization is true an optimization described in the LEDA book is used.
void	MULMULEY_SEGMENTS(list<segment> S, GRAPH<point,segment>& G, bool embed=false)
		as above. There is one additional output convention. If G is an undirected graph, the undirected planar map corresponding to G(s) is computed. The computation follows the incremental algorithm of Mulmuley ([62]) whose expected running time is O(M + s + nlog n), where n is the number of segments, s is the number of vertices of the graph G, and m is the number of edges.
void	SEGMENT_INTERSECTION(list<segment> S, void (*report)(segment, segment ))
		takes a list of segments S as input and executes for every pair (s1, s2) of intersecting segments report(s1, s2). The algorithm ([6]) has running time O(nlog2n + k), where n is the number of segments and k is the number intersecting pairs of segments.
void	SEGMENT_INTERSECTION(list<segment> S, list<point>& P)
		takes a list of segments S as input, computes the set of (proper) intersection points between all segments in S and stores this set in P. The algorithm ([10]) has running time O((| P| + | S|)log| S|).

$\bullet$ Red-Blue Line Segment Intersection

void	SEGMENT_INTERSECTION(list<segment> S1, list<segment> S2, GRAPH<point,segment>& G, bool embed=false)
		takes two lists of segments S1 and S2 as input and computes the planar graph G induced by the set of straight line segments in S1 $\cup$ S2 (as defined above). Precondition Any pair of segments in S1 or S2, respectively, does not intersect in a point different from one of the endpoints of the segments, i.e. segments of S1 or S2 are either pairwise disjoint or have a common endpoint.

$\bullet$ Closest Pairs

double	CLOSEST_PAIR(list<point>& L, point& r1, point& r2)
		CLOSEST_PAIR takes as input a list of points L. It computes a pair of points r1, r2 $\in$ L with minimal Euclidean distance and returns the squared distance between r1 and r2. The algorithm ([70]) has running time O(nlog n) where n is the number of input points.

$\bullet$ Miscellaneous Functions

void	Bounding_Box(list<point> L, point& pl, point& pb, point& pr, point& pt)
		computes four points pl, pb, pr, pt from L such that (xleft, ybot, xright, ytop) with xleft = pl.xcoord(), ybot = pb.ycoord(), xright = pr.xcoord() and ytop = pt.ycoord() is the smallest iso-oriented rectangle containing all points of L. Precondition L is not empty.
bool	Is_Simple_Polygon(list<point> L)
		takes as input a list of points L and returns true if Lis the vertex sequence of a simple polygon and false otherwise. The algorithms has running time O(nlog n), where n is the number of points in L.

$\bullet$ Properties of Geometric Graphs

We give procedures to check properties of geometric graphs. We give procedures to verify properties of geometric graphs. A geometric graph is a straight-line embedded map. Every node is mapped to a point in the plane and every dart is mapped to the line segment connecting its endpoints.

We use geo_graph as a template parameter for geometric graphs. Any instantiation of geo_graph must provide a function

VECTOR edge_vector(const geo_graph& G, const edge& e)

that returns a vector from the source to the target of e. In order to use any of these template functions the file /LEDA/templates/geo_check.t must be included.

template <class geo_graph>
bool	Is_CCW_Ordered(geo_graph G)
		returns true if for all nodes v the neighbors of v are in increasing counter-clockwise order around v.
template <class geo_graph>
bool	Is_CCW_Weakly_Ordered(geo_graph G)
		returns true if for all nodes v the neighbors of v are in non-decreasing counter-clockwise order around v.
template <class geo_graph>
bool	Is_CCW_Ordered_Plane_Map(geo_graph G)
		Equivalent to Is_Plane_Map(G) and Is_CCW_Ordered(G).
template <class geo_graph>
bool	Is_CCW_Weakly_Ordered_Plane_Map(geo_graph G)
		Equivalent to Is_Plane_Map(G) and Is_CCW_Weakly_Ordered(G).
template <class geo_graph>
void	SORT_EDGES(geo_graph& G)	Reorders the edges of G such that for every node v the edges in A(v) are in non-decreasing order by angle.
template <class geo_graph>
bool	Is_CCW_Convex_Face_Cycle(geo_graph G, edge e)
		returns true if the face cycle of G containing edefines a counter-clockwise convex polygon, i.e, if the face cycle forms a cyclically increasing sequence of edges according to the compare-by-angles ordering.
template <class geo_graph>
bool	Is_CCW_Weakly_Convex_Face_Cycle(geo_graph G, edge e)
		returns true if the face cycle of G containing edefines a counter-clockwise weakly convex polygon, i.e, if the face cycle forms a cyclically non-decreasing sequence of edges according to the compare-by-angles ordering.
template <class geo_graph>
bool	Is_CW_Convex_Face_Cycle(geo_graph G, edge e)
		returns true if the face cycle of G containing edefines a clockwise convex polygon, i.e, if the face cycle forms a cyclically decreasing sequence of edges according to the compare-by-angles ordering.
template <class geo_graph>
bool	Is_CW_Weakly_Convex_Face_Cycle(geo_graph G, edge e)
		returns true if the face cycle of G containing edefines a clockwise weakly convex polygon, i.e, if the face cycle forms a cyclically non-increasing sequence of edges according to the compare-by-angles ordering.