A CLOSED SET OF ALGORITHMS FOR PERFORMING SET OPERATIONS ON POLYGONAL REGIONS IN THE PLANE...

Preface to this interpretation

The excellent paper by Michael V. Leonov and Alexey G. Nikitin at the A. P. Ershov Institute of Information Systems, Siberian Branch of Russian back from the 1990s clearly and precisely presents a surprisingly simple algorithm for the problem also called polygon clipping. Besides describing the algorithm, the paper also gives accurate definitions and an estimation on runtime.

While the paper, in its original form, is already rather useful, I've found that as a coder, I have a hard time processing some parts of it, because of the following factors:

• mathematically precise, but often complex wording for expressing otherwise simple concepts
• condensed mathematical expressions
• file format: the original pdf seems to be broken (I had two different pdf reader implementation choke on it)
• while figures and example cases are chosen very wisely to present interesting cases, more labeling and coloring could make it easier to understand them

This document is an English-to-English translation of the original paper, with the intended target audience being coders who want to understand and implement the algorithm but do not wish to prove or verify the maths behind it. In fact, this document is more of my personal interpretation than a translation.

This document closely follows the structure of the original document, so that the two can be read in parallel.

I heavily simplified or even skipped parts that are very important for mathematical correctness but not required for implementing the algorithm in code. Instead of the mathematical notations, beyond the most basic operators I use a named operator function; instead of single letter variables, I use multiple-letter ASCII-only variable names (no greek letters). Instead of condensing the formulas, I rather make them more verbose, explaining what each variable is for. These modifications may degrade pedantic mathematical correctness, but I believe it makes it easier to digest the information on first pass.

I've also redrawn all illustrations so that they have more labels, more colors; more information explicitly encoded in them that otherwise could have been implied or figured from context.

When an algorithm is specified in pseudo-code or as expressions in tables, the original paper uses a Pascal-like syntax mxed with non-ASCII math notations, this paper uses an ASCII-only C-like syntax.

Finally, this document is plain, hand written html, which is much more portable than pdf.

Abstract

See the original paper.

Introduction

See the original paper.

The basic concept of the algorithm is simple: insert a new vertex (corner) in the contour of polygons ('A' and 'B') where they intersect (marked as p1 and p2 on Figure 1). Then split up the polyline of 'A' into contour sections between intersection points, name these contour sections CA1, CA2. Do the same on 'B', resulting contour sections are called CB1 and CB2.

Figure 1. Intersection polygons.

With this the problem is converted from geometry to topology. A few simple properties can be defined, such as whether points left of a contour section are within 'A' or 'B'. A contour section is "boring", since it does not have any intersection thus properties can not change. All the algorithm needs to do is figuring which contour sections to include in the result (and whether to reverse their direction so the result is a continous loop of polyline). For example an union operator will keep CA1 and CB1, an intersection will keep CA2 and CB2 and an "A-B" will keep CA1 and CB2 (with reversed direction).

Basic definitions

Our coordinate system is a 2d Euclidian plane with the usual Cartesian system. A point is specified by its x;y coordinates, but points are not stored, only used as a reference for this document. A vertex is also specified by its x;y coordinates and is stored in polylines which are stored in [polygon] regions in turn.

(Definition 1 and 2)

A [polygon] region consists of zero or more islands. (The original paper calls this polygon; this document will call it an island to reduce confusion with similar terms.) An island is an outer polyline and zero or more inner polylines (which act as cutouts). A polyline is an ordererd, cyclic list of vertices, containing at least 3 vertices.

A contour is a section of a polyline that contains 1 or more edges. (Unlike in the original paper, this document never uses contour as a synonym a whole polyline but only to reference a section of the polyline between two intersections. In the corner case when the polyline has no intersection, the contour is the whole polyline.)

An edge is a straight line segment¹ that starts from such a vertex and ends at the next vertex; an edge is referenced by its starting vertex. This also means edges are directed from start vertex to end vertex. There are no zero-legth edges.

Given such a directed edge, there is a left neighborhood and a right neighborhood (of imaginary points) next to the edge. (This concept is used only to decide if something is left or right to the edge.)

The straight ordering of vertices or edges within a polyline is counter-clockwise (CCW), the reverse ordering is clockwise (CW).

Figure 2. Neighborhoods. a) is an edge E with its left and right neighborhood; b) shows a vertex on the right between two edges of an outline, dividing the plane into an inner neighborhood and an outer neighborhood (of the vertex and of the edges).

A high degree vertex is shared by the endpoints of more than two edges.

(Definition 3)

Given a polyline, there's an inner neighborhood (all points inside the polyline) and an outer neighborhood (all points outside). A valid polyline can be self touching (a vertex between two edges falling on a 3rd edge of the polyline) but can not be self intersecting.

The outline of an island is always straight (CCW) while cutouts are always reversed (CW). Each island has exactly one outline, and that is the outmost polyline of the island. Polylines of an island do not touch or intersect other polylines within the island.

Figure 3. The left region has one outer polyline and one self-touching inner (cutout) polyline. The middle region has two islands of a single outer polylines each. The third region has a rectangular island with a trangular cutout and a second triangular island placed within that cutout.

Naming convention: this document will name vertices as V or Vsomething, edges as E or Esomething, and index them starting from 0 from an arbitrarily chosen vertex on the polyline, so that (for example) Efoo[5] is the edge that starts at Vfoo[0] and ends at Vfoo[1]. The last edge of a polyline of n vertices, E[n-1] is between V[n-1] and V[0]. Polylines are named Psomething; when indexed, Pfoo[0] is the outline, indices starting from 1 refer to cutouts.

A sector is an imaginary area, a neighborhood of points between two adjacent edges of a polyline. It is roughly the points within the triangle formed by the two edges and a third, virtual line connecting their non-incident ends. (Sectors are not calculated but are used in the explanation and in conceptual relations between edges and polylines.)

The algorithm

Input:

• A and B: two regions
• op: an operation (one of: union, intersection, difference, xor²)

The algorithm is split into four steps:

• Step 1. Processing of the edge intersections.
• Step 2. Edge and contour labeling.
• Step 3. Collecting the resulting contours.
• Step 4. Placing the resulting contours in R.

Figure 4. The two input regions and the protoype for result R, with all edges included. The vertex circled on the three regions are in the same point on the plane.

Step 1: processing edge intersections

Find all intersections between edge objects of A and B (ignoring the shared endpoint of adjacent edges within a polyline but not ignoring other intersections happening between edges of the same polyline):

• for an edge-edge mid intersection, insert a new vertex in both regions, at the intersection point, splitting the edge of both objects
• if a vertex of an island is on the edge of another, split that edge inserting a new vertex in the intersection point
• if two edges overlap, split them ensuring there's a vertex on the start and end of the overlap
• (if two vertices fall on the same point, don't do anything)

Figuer 5 shows the result. Vertices inserted for intersections are marked with bold. Input regions A and B are arranged so that the circled vertex (a5, B2, B9), marked with y in the result, is a common point. Note: a5 is on an inner (cutout) polyline in A, while B consists of a single outer polyline in which B2 and B9 happens to be on the same point in space (self touching polyline).

Figure 5. The two input regions with inserted cross-vertices marked with bold and the protoype for result R, with all edges included.

By now any intersection has an explicit vertex on all related polylines. If there are multiple vertices falling on the same point, that means there is an intersection/touching there and all vertices on that point are called cross-vertices.

Comment: for each such point of intersection there are at least two cross-vertices for such a point.

Comment: this process does change the structure of input regions, when it is creating new vertices, which means if the original regions are to be preserved, a copy needs to be made. However, this process (in theory^{3, 4}) does not change the shape/meaning/topology of the input regions.)

Take each cross-vertex V[i], where i is the index of the given vertex within the polyline. There will be a previous edge Ep that ends on that vertex and a next edge En that starts at that vertex. Allocate two cross-vertex descriptors (CVD), Dp and Dn, and assign them to Ep and En respectively (cross-linking the corresponding E and D with pointers in both directions). A CVD contains angle⁵ information and flags.

For every intersection point, collecd all CVDs into a cyclic list called the connectivity list [of that point]. The connectivity list is ordered by the angle of the CVDs (fig. 6). Order of equal angles is unspecified.

Figure 6. Cross linking CVDs in the point circled on fig. 5; top cyclic list is the outer polyline of region B, bottom cyclic list is the inner polyline of region A; middle is the cyclic list of DVDs for this point, sorted by angle; green arrows represent cross-linking-pointers between the cross vertices in that point (B2, B9 and a5) and the corresponding CVDs.

Comment: for a geometric point, like the one curcled in fig. 5, there is only one connectivity list and each cross vertices linked to the CVDs of that connectivity list are in that point. The same cross vertices can not be linked to any other connectivity list, because that other connectivity list would be for a different point.

Step 2. Edge and contour labeling

Let Pc be the outer contour of either input region A or input region B. Let Pm be the outer contour of the other input region.

(Predicate 1.) After the first step:

• Any pair of edges of A and B intersect at their endpoints (and may or may not overlap, but if they overlap, they overlap in full length)
• Every edge of Pc:
  • either overlaps with an edge in Pm (and is called a shared edge)
  • or lies fully inside Pm (except for an endpoint maybe)
  • or lies fully outside Pm (except for an endpoint maybe)
• if Pc does not have any cross-vertex, it is either fully inside Pm or is outside of Pm

The label on an edge E (stored as E.flags) of Pc is one of:

• INSIDE: when E is inside Pm
• OUTSIDE: when E is outside Pm
• SHARED1: when E overlaps with an edge from Pm that, and the two edges have the same direction
• SHARED2: when E overlaps with an edge from Pm that, and the two edges have different direction

The label of Pc (stored as Pc.flags) is one of:

• ISECTED: if it contains at least one cross-vertex
• INSIDE: if it has no cross-vertex and it is fully inside Pm
• OUTSIDE: if it has no cross-vertex and is outside of Pm

Algorithm for edge (and polyline) labeling

• 1. If Pc does not contain any cross-vertice: label Pc INSIDE or OUTSIDE⁶. Do not label edges of Pc.
• 2. If Pc contains at least one cross-vertex: label Pc ISECTED, and label its each edge E[i] (i is between 0 and n, n is the number of vertices in Pc):
  • 2.1. if E does not contain cross-vertices at its endpoints:
    • 2.1.2. if i == 0: if starting vertex lies inside Pm label E[i] INSIDE, else label E[i] OUTSIDE
    • 2.1.1. if i != 0: copy label value from E[i-1]
  • 2.2. if E contains 1 or 2 cross-vertices at its endpoints:
    • 2.2.1. there is an overlapping edge in Pm that has the same starting point as E and the same ending point as E: label E as SHARED1
    • 2.2.2. there is an overlapping edge in Pm that has the same starting point as the endpoint of E and the same ending point as the start point of E: label E as SHARED2
    • 2.2.3. Pc's cross-vertices connectivity lists do not contain any vertex from Pm: C touches itself but doesn't intersect M; jump to 2.1
    • 2.2.4. Pc's cross-vertices connectivity lists contains vertices (called Vm[j]) from Pm:
      • take the sector defined by Vm[j-1], Vm[j] and Vm[j+1]
      • check, using the angle of E: if E is inside that sector, label E as INSIDE, else label E as OUTSIDE (see Fig. 7).

Figure 7. Deciding if an edge being labeled is INSIDE or OUTSIDE of region M

Step 3. Collecting the resulting contours

An essential idea is that the polylines of R are collected using edge and polyline labels instead of using coordinates⁷.

rules for edge inclusion in result

Consider an edge E belonging to a pline in either input region A or B. Let NL(E) be the left neighborhood of E and NR(E) be the right neighborhood of E, where neighborhood is a theoretical construction and as a simplifiation can be considered to consists of arbitrary points very close to the edge. Such a neighborhood is either fully inside or fully outside of A (or B).

The rules for inclusion of edge E is described in table 1 and actions described in Table 1b.

	Edge labels		action per operator
			union	inters.	diff.	xor
1.		E is in input region A, E's label is OUTSIDE, NL(E) is inside input region A, NR(E) is not inside any region.	A+S	O	A+S	A+S
2.		E is in input region B, E's label is OUTSIDE, NL(E) is inside input region B, NR(E) is not inside any region.	A+S	O	O	A+S
3.		E is in input region A, E's label is INSIDE, NL(E) is inside both input regions A and B, NR(E) is inside input region B.	O	A+S	O	A+R
4.		E is in input region B, E's label is INSIDE, NL(E) is inside both input regions A and B, NR(E) is inside input region A.	O	A+S	A+R	A+R
5.		E from region A is overlapping edge F from region B and is labeled SHARED2; NL(E) is inside input region A, NR(F) is inside input region B. Pick either E or F for inclusion, use E's original direction for straight.8	O	A+S	O	O
6.		E from region A is overlapping edge F from region B and is labeled SHARED1; NL(E)==NL(F) is inside input both regions A and B, NR(E)==NR(F) is not inside either input region. Pick either E or F for inclusion.8	A+S	O	A+S	O

Comment: In other words, according to the above table, E is added in the resulting R:

• 1. when NL(E) is inside R and NR(E) is outside of R, the E is added with its original direction (fig. 1. example: union, case 1: edge is part of CA1)
• 2. when NL(E) is outside R and NR(E) is inside of R, the E is added with reverse direction (fig. 1. example: xor, case 5: edge is part of CB2)
• 3. when both NL(E) and NR(E) are outside R, E is not added (fig. 1. example: difference, case 2: edge is part of CB1)
• 4. when both NL(E) and NR(E) are inside R, E is not added (fig. 1. example: union, case 5: edge is part of CB2)

Comment: For shared edges only one of the edges is added.

jumping rules at intersections

If the collect step reaches a cross-vertex, that means we have reached an intersection where more than 2 vertices (4 or more edges) meet. The number of CVDs (descriptors from step 1) is twice the number of vertices at that point. Using the connectivity list of step 1, that also means we have a set of edges ordered by angle. Let each edge of this set be F(j), where j is between 0 and length of the point's connectivity list minus 1.

Points in the sector between F[j] and F[j+1] are either all inside or outside of region A and/or B. If the points of two neigbouring sectors are in the result R, by definition 3 their outer polyline contour will be the same.

Using this rule, we loop over cross-linked edges F[] with j until we find an edge that satisfies the above rule. In other words, we jump to the closest (in angle) edge that would make up a polyline that includes all the same-polarity neighboring sectors, using the rules detailed below.

The collecting algorithm

The same algorithm is used for all operations, except for some differences in inclusion rules.

Let Pc be the polyline currently being considered. The label of Pc is stored in Pc.flags and the label of an edge E is stored in E.flags. Direction FORWARD means the original direction (of E or Pc), BACKWARD means reversed direction.

1. The simpler case is when Pc.flags != ISECTED, which means the polyline has no intersection. In this case whether Pc is included in R and whether it is reversed depends on conditions described in Table 2.

op	is Pc included?	direction
union	Pc.flags == INSIDE	FORWARD
inters.	Pc.flags == OUTSIDE	FORWARD
diff.	(Pc is in region A) && (Pc.flags == OUTSIDE)	FORWARD
diff.	(Pc is in region B) && (Pc.flags == INSIDE)	BACKWARD
xor	Pc.flags == OUTSIDE	FORWARD
xor	Pc.flags == INSIDE	BACKWARD

2. If Pc.flags == ISECTED, the polyline intersects with another polyline. We need to find all edges in Pc which resulting contour section can start from. n is the number of vertices in Pc, E[] are edges of Pc; each E has a field called Mark that inidicates if E has been included in a resulting contour already (it is set to false initially on all E[]); v[] is each vertex in the polyline. The algorithm for starting edges is:

for(i=0; i < n; i++) {
	if (EdgeRule(E[i], &dir) && !E[i].Mark) {
		contour r; /* a contour section is a list of edges */
		if (dir == FORWARD)
			r = Collect(v[i], dir);
		else
			r = Collect(v[i+1], dir);
		include_into_a_set_of_resulting_contours(r);
	}
}

contour Collect(vertex v, enum dir)
{
	edge E, Esh, Enext;
	contour r = alloc_empty_contour();

	do {
		add_vertex_to_contour(r, v);
		if (dir == FORWARD)
			E = edge_next_to(v);
		else
			E = edge_prev_to(v);
		E.mark = true;
		if ((E.flags == SHARED1) || (E.flags == SHARED2)) {
			Esh = get_edge_shared_with(E);
			Esh.mark = true;
		}
		if (dir == FORWARD)
			v = vertex_next_to(v);
		else
			v = vertex_prev_to(v);
		if (is_cross_vertex(v))
			Jump(&v, &dir);

		/* look at the next edge to terminate; this part was not well specified
		   in the psudo-code of the original paper */
		if (dir == FORWARD)
			Enext = edge_next_to(v);
		else
			Enext = edge_prev_to(v);
	} while(!Enext.mark)
	return r;
}

void Jump(vertex *v, enum *dir)
{
	if (*dir == FORWARD)
		d = prev_D_in_CW(Dp[v]);   /* DP[v] is the previous on the connectivity list of v */
	else
		d = prev_D_in_CW(Dn[v]);   /* DP[v] is the previous on the connectivity list of v */

	for(;;) {
		edge e = edge_of_descriptor(d);
		if (!e.mark && EdgeRule(e, &newdir)) { /* EdgeRule(): see Table 3 */
			*v = vertex_of_descriptor(d);
			if (edge_is_next_to_vertex(e, *v) && (newdir == FORWARD)) ||
			   (edge_is_prev_to_vertex(e, *v) && (newdir == BACKWARD)) {
				*dir = newdir;
				return;
			}
		}
		d = prev_D_in_CW(d);
	}
}

/* prev_D_in_CW(Din) is the neatest D to Din in clockwise direction */