Polygon Algorithms, in depth: validation, winding, convexity, centroids, ear-clipping triangulation and clipping

What a polygon is to a program

A polygon in code is an ordered ring of vertices, and the edges are implied: vertex i connects to vertex i+1, and the last connects back to the first. Formats disagree on whether that closing edge is stored. GeoJSON and WKT repeat the first vertex at the end; most algorithm code does not. Mixing the two silently adds a zero-length edge, which breaks any test that divides by edge length or looks at the turn at a vertex.

Orientation is the second convention. The sign of the shoelace sum gives the winding: positive for counter-clockwise (CCW) with y pointing up. Screen coordinates flip y, which reverses it. GeoJSON (RFC 7946) asks for CCW exteriors but tells parsers not to reject the other order, so input cannot be trusted to follow it. Pick one internal convention (here, CCW exteriors) and enforce it at the boundary.

Algorithms also promise different things for different kinds of polygon:

• Convex: every interior angle is below 180 degrees; clipping is easy and any vertex fans into a triangulation. See Convex Hull, in depth.
• Simple: edges meet only where neighbours share an endpoint. Most algorithms here require this.
• With holes: one exterior ring plus interior rings, like most real footprints.
• Self-intersecting: a bowtie. Area is ambiguous and many algorithms return nonsense, not an error.

Normalise and validate before anything else

Normalising and validating first is the cheapest bug fix in geometry. Normalisation removes the noise; validation rejects input the later algorithms cannot handle. The orientation primitive cross(o, a, b) returns twice the signed area of triangle o-a-b: positive for a left turn, negative for a right turn, zero when the three points are collinear.

def cross(o, a, b):
    return (a[0]-o[0])*(b[1]-o[1]) - (a[1]-o[1])*(b[0]-o[0])

def signed_area2(P):  # twice the signed area; > 0 means CCW
    n = len(P)
    return sum(P[i][0]*P[(i+1)%n][1] - P[(i+1)%n][0]*P[i][1] for i in range(n))

def normalise(P):
    """Drop the repeated closing vertex and collinear vertices; return CCW order."""
    if len(P) > 1 and P[0] == P[-1]:
        P = P[:-1]
    out = [P[i] for i in range(len(P))
           if cross(P[i-1], P[i], P[(i+1) % len(P)]) != 0]
    if signed_area2(out) < 0:
        out.reverse()
    return out

def is_simple(P):
    """O(n^2) check: no two non-adjacent edges touch."""
    n = len(P)
    for i in range(n):
        for j in range(i+1, n):
            if j == i+1 or (i == 0 and j == n-1):
                continue          # neighbours share a vertex by design
            if seg_intersect(P[i], P[(i+1)%n], P[j], P[(j+1)%n]):
                return False
    return True

Here seg_intersect is the standard four-orientation segment test, treating touching as intersecting, so a vertex lying on a non-adjacent edge is rejected. The pairwise loop is O(n squared); for large rings, a Shamos-Hoey sweep line answers the same question in O(n log n).

Run it on the worked example, an L shape given in GeoJSON style with its closing vertex: [(0,0),(4,0),(4,1),(1,1),(1,3),(0,3),(0,0)]. Normalisation drops the duplicate and returns six vertices v0 to v5, already CCW; the area is 6.0 and is_simple returns True. Feed it the bowtie [(0,0),(2,2),(2,0),(0,2)] and is_simple returns False, while the shoelace formula reports an area of 0.0, because the two lobes cancel: a self-intersecting zone looks like an empty one. Remove collinear vertices for computation only; stored boundaries shared with neighbouring parcels may need them.

What a polygon is to a program

A polygon in code is an ordered ring of vertices, and the edges are implied: vertex i connects to vertex i+1, and the last connects back to the first. Formats disagree on whether that closing edge is stored. GeoJSON and WKT repeat the first vertex at the end; most algorithm code does not. Mixing the two silently adds a zero-length edge, which breaks any test that divides by edge length or looks at the turn at a vertex.

Orientation is the second convention. The sign of the shoelace sum gives the winding: positive for counter-clockwise (CCW) with y pointing up. Screen coordinates flip y, which reverses it. GeoJSON (RFC 7946) asks for CCW exteriors but tells parsers not to reject the other order, so input cannot be trusted to follow it. Pick one internal convention (here, CCW exteriors) and enforce it at the boundary.

Algorithms also promise different things for different kinds of polygon:

• Convex: every interior angle is below 180 degrees; clipping is easy and any vertex fans into a triangulation. See Convex Hull, in depth.
• Simple: edges meet only where neighbours share an endpoint. Most algorithms here require this.
• With holes: one exterior ring plus interior rings, like most real footprints.
• Self-intersecting: a bowtie. Area is ambiguous and many algorithms return nonsense, not an error.

Normalise and validate before anything else

Normalising and validating first is the cheapest bug fix in geometry. Normalisation removes the noise; validation rejects input the later algorithms cannot handle. The orientation primitive cross(o, a, b) returns twice the signed area of triangle o-a-b: positive for a left turn, negative for a right turn, zero when the three points are collinear.

def cross(o, a, b):
    return (a[0]-o[0])*(b[1]-o[1]) - (a[1]-o[1])*(b[0]-o[0])

def signed_area2(P):  # twice the signed area; > 0 means CCW
    n = len(P)
    return sum(P[i][0]*P[(i+1)%n][1] - P[(i+1)%n][0]*P[i][1] for i in range(n))

def normalise(P):
    """Drop the repeated closing vertex and collinear vertices; return CCW order."""
    if len(P) > 1 and P[0] == P[-1]:
        P = P[:-1]
    out = [P[i] for i in range(len(P))
           if cross(P[i-1], P[i], P[(i+1) % len(P)]) != 0]
    if signed_area2(out) < 0:
        out.reverse()
    return out

def is_simple(P):
    """O(n^2) check: no two non-adjacent edges touch."""
    n = len(P)
    for i in range(n):
        for j in range(i+1, n):
            if j == i+1 or (i == 0 and j == n-1):
                continue          # neighbours share a vertex by design
            if seg_intersect(P[i], P[(i+1)%n], P[j], P[(j+1)%n]):
                return False
    return True

Here seg_intersect is the standard four-orientation segment test, treating touching as intersecting, so a vertex lying on a non-adjacent edge is rejected. The pairwise loop is O(n squared); for large rings, a Shamos-Hoey sweep line answers the same question in O(n log n).

Run it on the worked example, an L shape given in GeoJSON style with its closing vertex: [(0,0),(4,0),(4,1),(1,1),(1,3),(0,3),(0,0)]. Normalisation drops the duplicate and returns six vertices v0 to v5, already CCW; the area is 6.0 and is_simple returns True. Feed it the bowtie [(0,0),(2,2),(2,0),(0,2)] and is_simple returns False, while the shoelace formula reports an area of 0.0, because the two lobes cancel: a self-intersecting zone looks like an empty one. Remove collinear vertices for computation only; stored boundaries shared with neighbouring parcels may need them.

Convexity and the true centroid

A polygon is convex when every turn has the same sign. Skipping zero turns makes the test tolerate collinear points, and requiring simplicity beforehand matters: a star-shaped self-intersecting pentagram turns the same way at every vertex and would pass a turn-sign test alone.

def is_convex(P):
    sign = 0
    for i in range(len(P)):
        z = cross(P[i], P[(i+1) % len(P)], P[(i+2) % len(P)])
        if z != 0:
            if sign == 0:
                sign = 1 if z > 0 else -1
            elif (z > 0) != (sign > 0):
                return False
    return True

def centroid(P):
    """Area centroid of a simple polygon (not the vertex average)."""
    a2 = signed_area2(P)
    cx = cy = 0.0
    for i in range(len(P)):
        (x0, y0), (x1, y1) = P[i], P[(i+1) % len(P)]
        f = x0*y1 - x1*y0
        cx += (x0 + x1) * f
        cy += (y0 + y1) * f
    return (cx / (3*a2), cy / (3*a2))

On the L shape, is_convex returns False, because the turn at v3 = (1,1) is a right turn: v3 is a reflex vertex. The centroid comes out at (1.5, 1.0). You can check it by decomposition: the 4 by 1 bar has area 4 and centre (2, 0.5); the 1 by 2 upright has area 2 and centre (0.5, 2); the area-weighted mean is (9/6, 6/6) = (1.5, 1.0).

Two traps: the vertex average is not the centroid (densely digitised curves drag it), and the centroid of a concave polygon such as a C shape can lie outside it. For label anchors use a guaranteed interior point; Shapely calls it representative_point().

Ear-clipping triangulation

Triangulation feeds GPUs, meshes and anything easy on triangles, such as uniform random sampling. A simple polygon with n vertices always splits into n-2 triangles.

Ear clipping rests on the Two Ears Theorem (Meisters, 1975): every simple polygon with more than three vertices has at least two non-overlapping ears. An ear is a convex vertex b whose neighbours a and c can be joined by a diagonal lying inside the polygon, which holds exactly when no other vertex lies inside triangle a-b-c. Cut the ear off, and what remains is still simple, so it has ears too. Repeat until one triangle is left.

def in_triangle(p, a, b, c):          # inclusive: points on an edge count as inside
    return cross(a, b, p) >= 0 and cross(b, c, p) >= 0 and cross(c, a, p) >= 0

def ear_clip(P):
    """P: simple, CCW, no collinear vertices. Returns index triples."""
    idx, tris = list(range(len(P))), []
    while len(idx) > 3:
        n = len(idx)
        for k in range(n):
            i, j, l = idx[k-1], idx[k], idx[(k+1) % n]
            a, b, c = P[i], P[j], P[l]
            if cross(a, b, c) <= 0:
                continue                   # reflex or flat: not an ear
            if any(in_triangle(P[m], a, b, c) for m in idx if m not in (i, j, l)):
                continue                   # another vertex blocks the diagonal
            tris.append((i, j, l))
            del idx[k]
            break
        else:
            raise ValueError("no ear found: input not simple or not CCW")
    tris.append(tuple(idx))
    return tris

On the L shape this returns [(0,1,2), (0,2,3), (5,0,3), (3,4,5)]: four triangles, as n-2 predicts, with areas 2, 1.5, 1.5 and 1 summing to the polygon's 6. The first candidate, v0, is rejected because triangle v5-v0-v1 contains v3; v1 passes and is clipped. The else branch matters: on subtly invalid input, a version without it loops forever.

This version is O(n cubed) in the worst case. Keeping a linked list of current ears and re-testing only the neighbours of each clipped vertex gives O(n squared), and only reflex vertices need the containment check. Ear clipping produces slivers; when triangle quality matters, use constrained Delaunay triangulation instead.

Polygons with holes: bridge edges

Ear clipping needs one ring. A polygon with holes is converted by cutting a bridge: pick the hole's vertex with the largest x, cast a ray toward +x, find the outer-boundary edge it hits first, and choose a vertex on that side that is mutually visible with the hole vertex. Splice the hole (wound clockwise) into the outer ring by walking out along the bridge, around the hole, and back along the same bridge. The result is a single ring with two coincident edges; the containment test must tolerate the duplicated vertices. Process holes in decreasing order of maximum x so bridges cannot cross. David Eberly documents this method, and Mapbox's earcut library is built on the same idea.

Clipping against a convex window

Clipping keeps the part of a subject polygon inside a clip region, such as a parcel inside a flood zone. For a convex region, Sutherland-Hodgman (1974) clips against one window edge at a time as a half-plane, walking the subject's edges with four cases: inside to inside emits the endpoint, inside to outside emits the crossing, outside to inside emits the crossing and the endpoint, outside to outside emits nothing.

def clip(subject, window):
    """Sutherland-Hodgman. window must be convex and CCW."""
    out = subject
    for i in range(len(window)):
        a, b = window[i], window[(i+1) % len(window)]
        inp, out = out, []
        for j in range(len(inp)):
            cur, prev = inp[j], inp[j-1]
            cur_in, prev_in = cross(a, b, cur) >= 0, cross(a, b, prev) >= 0
            if cur_in:
                if not prev_in:
                    out.append(crossing(prev, cur, a, b))
                out.append(cur)
            elif prev_in:
                out.append(crossing(prev, cur, a, b))
        if not out:
            break
    return out

def crossing(p, q, a, b):
    d1, d2 = cross(a, b, p), cross(a, b, q)
    t = d1 / (d1 - d2)
    return (p[0] + t*(q[0]-p[0]), p[1] + t*(q[1]-p[1]))

Clip the L shape against the rectangle with corners (0.5, 0.5) and (3, 2) and you get [(0.5,2.0), (0.5,0.5), (3.0,0.5), (3.0,1.0), (1,1), (1.0,2.0)]: a smaller L with area 1.25 + 0.5 = 1.75, which you can confirm from the diagram.

Limits: a concave window gives the wrong region, and a concave subject cut into two pieces comes back as one ring joined by zero-width edges along the window boundary, harmless for rendering and wrong for topology.

General boolean operations and libraries

Intersection, union and difference of arbitrary polygons need every edge crossing computed and the output traced. Weiler-Atherton (1977) and Greiner-Hormann (1998) walk between rings at crossings, the latter assuming no degenerate touches in its original form. Vatti's sweep (1992) underlies the Clipper library, and GEOS, behind Shapely and PostGIS, uses a noding overlay. Degeneracies are most of the work, so use a library:

Robustness: where the sign of cross goes wrong

Every algorithm above decides from the sign of cross. With integers that sign is exact (Python integers never overflow). With doubles, nearly collinear points give rounding noise around zero, and two code paths can disagree about which side a point is on. That inconsistency is what breaks geometry: a vertex is inside an ear for one test and outside for the next. Defences, in order of preference:

1. Snap to a grid and use integers; map data at 1 cm precision fits easily in 64 bits. Clipper2 works this way.
2. For floats, use Shewchuk's adaptive-precision orientation predicate, which returns the exact sign.
3. Avoid ad hoc epsilons: they make one test lenient and its neighbours strict. Snap once, then compute exactly.
4. Re-validate constructed output; rounded intersection points can make a clip result self-intersect.

Failure modes and operational checks

Symptoms you will meet, and their usual causes:

Log invalid, repaired and rejected counts per source at ingest instead of repairing silently; a jump usually means an upstream editor or projection changed. For many polygons, index bounding boxes in an R-tree or a k-d tree so exact tests run only on candidates.

What to do next

Build the checks before the algorithms:

1. Run normalise and is_simple on a sample of real data and count closed, clockwise and self-intersecting rings.
2. Fix one internal convention (open rings, CCW exteriors, units) and enforce it at ingest with tests.
3. Re-run the L-shape example (area 6, centroid (1.5, 1.0), four triangles, clip area 1.75), then confirm the bowtie is rejected.
4. Property-test triangulation: for random simple polygons, check that there are n-2 triangles, all CCW, with areas that sum to the polygon area.
5. Decide where exactness comes from: an integer grid, robust predicates, or a library such as GEOS, Clipper2 or CGAL. Use it for boolean operations instead of extending Sutherland-Hodgman.
6. Keep learning: orientation and point-in-polygon in 2D geometry, in depth, exact hull construction in the Graham scan article.