


(A) Test if Simple. The ShamosHoey algorithm can be used to test if a polygon is simple or not. We give a C++ implementation simple_Polygon() for this algorithm below. Note that the shared endpoint between sequential edges does not count as a nonsimple intersection point, and the intersection test routine must check for that.
Also, we do not provide the code for a balanced binary tree, which is needed for the sweepline data structure. We did this to concentrate only on the geometric aspect of the algorithm, and not recommend a specific type of balanced tree, such as AVL, or redblack, or 23. These can be found in standard libraries. Nevertheless, there have often been requests to include a complete standalone algorithm.
Recently, a complete simple_Polygon() algorithm using an AVL tree has been developed by [Glenn Burkhardt, 2014], and provided to us for publication here. He used AVL tree code previously developed by [Brad Appleton, 1997], who gave his permission to republish his code here (including his License agreement). In addition to integrating these two software packages, Glenn also made modifications for greater efficiency, and to make the code compatible with the current C++ standards.
(B) Decompose into Simple Pieces. The BentleyOttmann algorithm can be used to decompose a nonsimple polygon into simple pieces. To do this, all intersection points are needed. One approach for a simple decomposition algorithm is to perform “surgery” at each intersection point. This proceedure can be incorporated into the sweepline algorithm as it discovers new intersection points, resulting in an decomposition algorithm. To do this, one must also output the edges that persist in a linked list, whose connected sublists will be the new simple polygons.
Let the original polygon be given by the set of segment endpoints , with . Its directed edges are given by the vectors .
Note that the resulting simple polygons may not be disjoint since one could be contained inside another. In fact, the decomposition inclusion hierarchy is based on the inclusion winding number of each simple polygon in the original nonsimple one (see Algorithm 3 about Winding Number Inclusion). For example:
The BentleyOttmann algorithm can be used to speed up computing the intersection, union, or difference of two general nonconvex simple polygons. Of course, before using any complicated algorithm to perform these operations, one should first test the bounding boxes or spheres of the polygons for overlap (see Algorithm 8 on Bounding Containers). If the bounding containers are disjoint, then so are the two polygons, and the set operations become trivial.
However, when two polygons overlap, the sweep line strategy of the BentleyOttmann algorithm can be adapted to perform a set operation on any two simple polygons. For further details see [O'Rourke, 1998, 266269]. If the two polygons are known to be simple, then one just needs intersections for segments from different polygons, which is a redblue intersection problem.
The BentleyOttmann algorithm can be used to efficiently compute the overlay of two planar subdivisions. For details, see [de Berg et al, 2000, 3339]. A planar subdivision is a planar graph with straight line segments for edges, and it divides the plane into a finite number of regions. For example, boundary lines divide a country into states. When two such planar graphs are overlaid (or superimposed), then their combined graph defines a subdivision refinement of each one. To compute this refinement, one needs to calculate all intersections between the line segments in both graphs. For a segment in one graph, we only need the intersections with segments in the other graph, and so this is another redblue intersection problem.
Here are some sample "C++" implementations of these algorithms.
// Copyright 2001 softSurfer, 2012 Dan Sunday
// This code may be freely used, distributed and modified for any purpose
// providing that this copyright notice is included with it.
// SoftSurfer makes no warranty for this code, and cannot be held
// liable for any real or imagined damage resulting from its use.
// Users of this code must verify correctness for their application.
// Assume that classes are already given for the objects:
// Point with 2D coordinates {float x, y;}
// Polygon with n vertices {int n; Point *V;} with V[n]=V[0]
// Tnode is a node element structure for a BBT
// BBT is a class for a Balanced Binary Tree
// such as an AVL, a 23, or a redblack tree
// with methods given by the placeholder code:
typedef struct _BBTnode Tnode;
struct _BBTnode {
void* val;
// plus node mgmt info ...
};
class BBT {
Tnode *root;
public:
BBT() {root = (Tnode*)0;} // constructor
~BBT() {freetree();} // destructor
Tnode* insert( void* ){}; // insert data into the tree
Tnode* find( void* ){}; // find data from the tree
Tnode* next( Tnode* ){}; // get next tree node
Tnode* prev( Tnode* ){}; // get previous tree node
void remove( Tnode* ){}; // remove node from the tree
void freetree(){}; // free all tree data structs
};
// NOTE:
// Code for these methods must be provided for the algorithm to work.
// We have not provided it since binary tree algorithms are wellknown
// and code is widely available. Further, we want to reduce the clutter
// accompanying the essential sweep line algorithm.
/**
Recently, a complete simple_Polygon() algorithm using an AVL tree has been developed by
[Glenn Burkhardt, 2014], and provided to us for publication here. He integrated our code
with AVL tree code previously developed by [Brad Appleton, 1997], and modified both to
improve the algorithm’s efficiency.
**/
//===================================================================
#define FALSE 0
#define TRUE 1
#define LEFT 0
#define RIGHT 1
extern void
qsort(void*, unsigned, unsigned, int(*)(const void*,const void*));
// xyorder(): determines the xy lexicographical order of two points
// returns: (+1) if p1 > p2; (1) if p1 < p2; and 0 if equal
int xyorder( Point* p1, Point* p2 )
{
// test the xcoord first
if (p1>x > p2>x) return 1;
if (p1>x < p2>x) return (1);
// and test the ycoord second
if (p1>y > p2>y) return 1;
if (p1>y < p2>y) return (1);
// when you exclude all other possibilities, what remains is...
return 0; // they are the same point
}
// isLeft(): tests if point P2 is LeftOnRight of the line P0 to P1.
// returns: >0 for left, 0 for on, and <0 for right of the line.
// (see Algorithm 1 on Area of Triangles)
inline float
isLeft( Point P0, Point P1, Point P2 )
{
return (P1.x  P0.x)*(P2.y  P0.y)  (P2.x  P0.x)*(P1.y  P0.y);
}
//===================================================================
// EventQueue Class
// Event element data struct
typedef struct _event Event;
struct _event {
int edge; // polygon edge i is V[i] to V[i+1]
int type; // event type: LEFT or RIGHT vertex
Point* eV; // event vertex
};
int E_compare( const void* v1, const void* v2 ) // qsort compare two events
{
Event** pe1 = (Event**)v1;
Event** pe2 = (Event**)v2;
return xyorder( (*pe1)>eV, (*pe2)>eV );
}
// the EventQueue is a presorted array (no insertions needed)
class EventQueue {
int ne; // total number of events in array
int ix; // index of next event on queue
Event* Edata; // array of all events
Event** Eq; // sorted list of event pointers
public:
EventQueue(Polygon P); // constructor
~EventQueue(void) // destructor
{ delete[] Eq; delete[] Edata;}
Event* next(); // next event on queue
};
// EventQueue Routines
EventQueue::EventQueue( Polygon P )
{
ix = 0;
ne = 2 * P.n; // 2 vertex events for each edge
Edata = (Event*)new Event[ne];
Eq = (Event**)new (Event*)[ne];
for (int i=0; i < ne; i++) // init Eq array pointers
Eq[i] = &Edata[i];
// Initialize event queue with edge segment endpoints
for (int i=0; i < P.n; i++) { // init data for edge i
Eq[2*i]>edge = i;
Eq[2*i+1]>edge = i;
Eq[2*i]>eV = &(P.V[i]);
Eq[2*i+1]>eV = &(P.V[i+1]);
if (xyorder( &P.V[i], &P.V[i+1]) < 0) { // determine type
Eq[2*i]>type = LEFT;
Eq[2*i+1]>type = RIGHT;
}
else {
Eq[2*i]>type = RIGHT;
Eq[2*i+1]>type = LEFT;
}
}
// Sort Eq[] by increasing x and y
qsort( Eq, ne, sizeof(Event*), E_compare );
}
Event* EventQueue::next()
{
if (ix >= ne)
return (Event*)0;
else
return Eq[ix++];
}
//===================================================================
// SweepLine Class
// SweepLine segment data struct
typedef struct _SL_segment SLseg;
struct _SL_segment {
int edge; // polygon edge i is V[i] to V[i+1]
Point lP; // leftmost vertex point
Point rP; // rightmost vertex point
SLseg* above; // segment above this one
SLseg* below; // segment below this one
};
// the Sweep Line itself
class SweepLine {
int nv; // number of vertices in polygon
Polygon* Pn; // initial Polygon
BBT Tree; // balanced binary tree
public:
SweepLine(Polygon P) // constructor
{ nv = P.n; Pn = &P; }
~SweepLine(void) // destructor
{ Tree.freetree();}
SLseg* add( Event* );
SLseg* find( Event* );
int intersect( SLseg*, SLseg* );
void remove( SLseg* );
};
SLseg* SweepLine::add( Event* E )
{
// fill in SLseg element data
SLseg* s = new SLseg;
s>edge = E>edge;
// if it is being added, then it must be a LEFT edge event
// but need to determine which endpoint is the left one
Point* v1 = &(Pn>V[s>edge]);
Point* v2 = &(Pn>V[s>edge+1]);
if (xyorder( v1, v2) < 0) { // determine which is leftmost
s>lP = *v1;
s>rP = *v2;
}
else {
s>rP = *v1;
s>lP = *v2;
}
s>above = (SLseg*)0;
s>below = (SLseg*)0;
// add a node to the balanced binary tree
Tnode* nd = Tree.insert(s);
Tnode* nx = Tree.next(nd);
Tnode* np = Tree.prev(nd);
if (nx != (Tnode*)0) {
s>above = (SLseg*)nx>val;
s>above>below = s;
}
if (np != (Tnode*)0) {
s>below = (SLseg*)np>val;
s>below>above = s;
}
return s;
}
SLseg* SweepLine::find( Event* E )
{
// need a segment to find it in the tree
SLseg* s = new SLseg;
s>edge = E>edge;
s>above = (SLseg*)0;
s>below = (SLseg*)0;
Tnode* nd = Tree.find(s);
delete s;
if (nd == (Tnode*)0)
return (SLseg*)0;
return (SLseg*)nd>val;
}
void SweepLine::remove( SLseg* s )
{
// remove the node from the balanced binary tree
Tnode* nd = Tree.find(s);
if (nd == (Tnode*)0)
return; // not there
// get the above and below segments pointing to each other
Tnode* nx = Tree.next(nd);
if (nx != (Tnode*)0) {
SLseg* sx = (SLseg*)(nx>val);
sx>below = s>below;
}
Tnode* np = Tree.prev(nd);
if (np != (Tnode*)0) {
SLseg* sp = (SLseg*)(np>val);
sp>above = s>above;
}
Tree.remove(nd); // now can safely remove it
delete s;
}
// test intersect of 2 segments and return: 0=none, 1=intersect
int SweepLine::intersect( SLseg* s1, SLseg* s2)
{
if (s1 == (SLseg*)0  s2 == (SLseg*)0)
return FALSE; // no intersect if either segment doesn't exist
// check for consecutive edges in polygon
int e1 = s1>edge;
int e2 = s2>edge;
if (((e1+1)%nv == e2)  (e1 == (e2+1)%nv))
return FALSE; // no nonsimple intersect since consecutive
// test for existence of an intersect point
float lsign, rsign;
lsign = isLeft(s1>lP, s1>rP, s2>lP); // s2 left point sign
rsign = isLeft(s1>lP, s1>rP, s2>rP); // s2 right point sign
if (lsign * rsign > 0) // s2 endpoints have same sign relative to s1
return FALSE; // => on same side => no intersect is possible
lsign = isLeft(s2>lP, s2>rP, s1>lP); // s1 left point sign
rsign = isLeft(s2>lP, s2>rP, s1>rP); // s1 right point sign
if (lsign * rsign > 0) // s1 endpoints have same sign relative to s2
return FALSE; // => on same side => no intersect is possible
// the segments s1 and s2 straddle each other
return TRUE; // => an intersect exists
}
//===================================================================
// simple_Polygon(): test if a Polygon is simple or not
// Input: Pn = a polygon with n vertices V[]
// Return: FALSE(0) = is NOT simple
// TRUE(1) = IS simple
int
simple_Polygon( Polygon Pn )
{
EventQueue Eq(Pn);
SweepLine SL(Pn);
Event* e; // the current event
SLseg* s; // the current SL segment
// This loop processes all events in the sorted queue
// Events are only left or right vertices since
// No new events will be added (an intersect => Done)
while (e = Eq.next()) { // while there are events
if (e>type == LEFT) { // process a left vertex
s = SL.add(e); // add it to the sweep line
if (SL.intersect( s, s>above))
return FALSE; // Pn is NOT simple
if (SL.intersect( s, s>below))
return FALSE; // Pn is NOT simple
}
else { // processs a right vertex
s = SL.find(e);
if (SL.intersect( s>above, s>below))
return FALSE; // Pn is NOT simple
SL.remove(s); // remove it from the sweep line
}
}
return TRUE; // Pn IS simple
}
//===================================================================
Brad Appleton, C++ code for an AVL balanced Tree, see: www.bradapp.com, oopweb.com/Algorithms/Documents/AvlTrees/VolumeFrames.html, www.bradapp.com/ftp/src/libs/C++/AvlTrees.html (1997)
I.J. Balaban, "An Optimal Algorithm for Finding Segment Intersections", Proc. 11th Ann. ACM Sympos. Comp. Geom., 211219 (1995)
Ulrike Bartuschka, Kurt Mehlhorn & Stefan Naher, "A Robust and Efficient Implementation of a Sweep Line Algorithm for the Straight Line Segment Intersection Problem", Proc. Workshop on Algor. Engineering, Venice, Italy, 124135 (1997)
Jon Bentley & Thomas Ottmann, "Algorithms for Reporting and Counting Geometric Intersections", IEEE Trans. Computers C28, 643647 (1979)
Mark de Berg et al, Computational Geometry : Algorithms and Applications (2nd Ed), Chapter 2 "Line Segment Intersection" (2010)
Glenn Burkhardt, personal communication (2014).
Tim Chan, "A Simple Trapezoid Sweep Algorithm for Reporting Red/Blue Segment Intersections", Proc. 6th Can. Conf. Comp. Geom., Saskatoon, Saskatchewan, Canada, 263268 (1994)
Bernard Chazelle & Herbert Edelsbrunner, "An Optimal Algorithm for Intersecting Line Segments in the Plane", Proc. 29th Ann. IEEE Sympos. Found. Comp. Sci., 590600 (1988)
Bernard Chazelle & Herbert Edelsbrunner, "An Optimal Algorithm for Intersecting Line Segments in the Plane", J. ACM 39, 154 (1992)
K.L. Clarkson & P.W. Shor, "Applications of Random Sampling in Computational Geometry, II", Discrete Comp. Geom. 4, 387421 (1989)
John Hobby, "Practical Segment Intersection with Finite Precision Output", Comp. Geom. : Theory & Applics. 13(4), (1999) [Note: the original Bell Labs paper appeared in 1993]
H.G. Mairson & J. Stolfi, "Reporting and Counting Intersections between Two Sets of Line Segments", in: Theoretic Found. of Comp. Graphics and CAD, NATO ASI Series Vol. F40, 307326 (1988)
E. Myers, "An O(E log E + I) Expected Time Algorithm for the Planar Segment Intersection Problem", SIAM J. Comput., 625636 (1985)
Joseph O'Rourke, Computational Geometry in C (2nd Edition), Section 7.7 "Intersection of Segments" (1998)
Franco Preparata & Michael Shamos, Computational Geometry: An Introduction, Chapter 7 "Intersections" (1985)
Michael Shamos & Dan Hoey, "Geometric Intersection Problems", Proc. 17th Ann. Conf. Found. Comp. Sci., 208215 (1976)
© Copyright 2012 Dan Sunday, 2001 softSurfer


