2D Geometry Algorithms, in depth: orientation, segment and line intersection, signed area, point in polygon, angle sorting and robust arithmetic

Vectors, dot and cross products

Represent a point as a pair (x, y) and a vector as the difference of two points. Two products do almost all the work. The dot product u.x*v.x + u.y*v.y is positive for an acute angle, zero for perpendicular vectors and negative for an obtuse one. The cross product u.x*v.y - u.y*v.x is the signed area of the parallelogram the two vectors span. Its sign says whether v lies counter-clockwise (positive) or clockwise (negative) from u, and it is zero exactly when the vectors are parallel.

Neither product divides or takes a square root, so with integer inputs both are exact, and so is any decision built from their signs. Division and square roots belong at the end of a computation, never inside a branch. The distance article, point-to-line and point-to-segment distance, derives the projection formula from the dot product and shows the same cross-product area trick from the distance side.

The orientation predicate

The orientation predicate takes three points and returns the sign of cross(b - a, c - a): +1 if the path a, b, c turns left, -1 if it turns right, 0 if the three are collinear. Hulls, sweeps and every test below reduce to it. The convex hull article shows how much of the monotone chain algorithm is this one sign.

For the worked polygon below, orient((0,0), (6,0), (6,4)) has cross product 24, a left turn, which is what a counter-clockwise boundary does at a convex corner. orient((6,4), (3,2), (0,4)) gives -12, a right turn, marking the reflex vertex at (3,2) where the polygon dents inward. orient((0,0), (2,2), (4,4)) gives exactly 0.

def sub(a, b): return (a[0] - b[0], a[1] - b[1])
def cross(u, v): return u[0] * v[1] - u[1] * v[0]
def dot(u, v): return u[0] * v[0] + u[1] * v[1]

def orient(a, b, c):
    """+1 if a->b->c turns left, -1 if right, 0 if collinear."""
    v = cross(sub(b, a), sub(c, a))
    return (v > 0) - (v < 0)

Segment intersection, including the degenerate cases

Do segments p1-p2 and q1-q2 intersect? In the general case they cross exactly when each segment's endpoints lie on opposite sides of the other segment's line: orient(q1, q2, p1) and orient(q1, q2, p2) have opposite signs, and orient(p1, p2, q1) and orient(p1, p2, q2) have opposite signs. That covers a proper X-shaped crossing. Everything else is a degenerate case, and real data is full of them: shared endpoints, junctions, long collinear boundaries.

The fix is to handle each zero separately. If an orientation is zero, that endpoint is collinear with the other segment, and the segments touch exactly when the endpoint also lies within the other segment's bounding box. The point is already on the line, so the box test suffices. This rule covers T-junctions, shared endpoints and collinear overlap.

def on_segment(p, a, b):
    """p is known collinear with a-b; is it within the bounding box?"""
    return (min(a[0], b[0]) <= p[0] <= max(a[0], b[0]) and
            min(a[1], b[1]) <= p[1] <= max(a[1], b[1]))

def segments_intersect(p1, p2, q1, q2):
    d1, d2 = orient(q1, q2, p1), orient(q1, q2, p2)
    d3, d4 = orient(p1, p2, q1), orient(p1, p2, q2)
    if d1 * d2 < 0 and d3 * d4 < 0:
        return True                       # proper crossing
    if d1 == 0 and on_segment(p1, q1, q2): return True
    if d2 == 0 and on_segment(p2, q1, q2): return True
    if d3 == 0 and on_segment(q1, p1, p2): return True
    if d4 == 0 and on_segment(q2, p1, p2): return True
    return False

Checked by running the code: (0,0)-(4,4) against (0,4)-(4,0) returns True, a proper crossing. The collinear pair (0,0)-(4,0) and (2,0)-(6,0) returns True through the overlap rule, while (0,0)-(4,0) and (5,0)-(6,0), collinear but separated, returns False. (0,0)-(2,2) against (3,3)-(5,1) returns False even though the infinite lines meet at (3,3): the point is on the second segment but past the end of the first. Testing only d1*d2 and d3*d4 misses the overlap case.

Where two lines meet

Knowing that two segments meet is a predicate; knowing where is a construction. Write the first segment as p1 + t*r with r = p2 - p1 and the second as q1 + u*s with s = q2 - q1. Crossing both sides with s eliminates u, giving t = cross(q1 - p1, s) / cross(r, s). When cross(r, s) is zero the lines are parallel, or collinear if cross(q1 - p1, r) is also zero, and there is no single intersection point to return.

from fractions import Fraction

def line_intersection(p1, p2, q1, q2):
    """Intersection of the infinite lines as exact fractions; None if parallel."""
    r, s = sub(p2, p1), sub(q2, q1)
    den = cross(r, s)
    if den == 0:
        return None                       # parallel or collinear
    t = Fraction(cross(sub(q1, p1), s), den)
    return (p1[0] + t * r[0], p1[1] + t * r[1])

For (0,0)-(4,4) and (0,4)-(4,0) the result is (2,2). For (0,0)-(4,2) and (0,3)-(3,0), r = (4,2), s = (3,-3), cross(r, s) = -18 and cross(q1 - p1, s) = -9, so t = 1/2 and the point is (2,1). The result is rational even for integer inputs, hence fractions. If t and the matching u both lie in [0, 1], the point is on both segments. Keep the exact rational when later decisions use the point, and convert to floating point only for display.

Signed area with the shoelace formula

The shoelace formula gives twice the signed area of a simple polygon as the sum, over consecutive vertices, of cross(v[i], v[i+1]), wrapping from the last vertex to the first. Each term is the signed area of the triangle from the origin to one edge; the parts outside the polygon cancel. The sign of the total tells you the winding direction: positive for counter-clockwise, negative for clockwise.

def area2(poly):
    """Twice the signed area (shoelace). Positive = counter-clockwise."""
    n = len(poly)
    return sum(cross(poly[i], poly[(i + 1) % n]) for i in range(n))

For the worked polygon (0,0), (6,0), (6,4), (3,2), (0,4) the terms are 0, 24, 0, 12 and 0. The sum is 36, so the area is 18, and the vertex list is counter-clockwise. Reversing the list gives -36. Keep the doubled value as an integer and halve only for output, so integer inputs stay exact. Many algorithms assume counter-clockwise input, so reverse the list when area2 is negative. The formula is wrong for self-intersecting polygons, so validate user input first.

Point in polygon with the winding number

Is a point inside a polygon? Ray crossing counts edge crossings of a ray from the point, odd meaning inside, but the textbook version mishandles rays through vertices and never reports boundary points. The winding number counts how many times the boundary winds around the point. Walking the edges, an upward edge that has the point strictly on its left adds one, and a downward edge that has the point strictly on its right subtracts one. Non-zero means inside.

Two details make it exact. Treat each edge as half-open in y (it includes its lower endpoint and excludes its upper endpoint), so a horizontal ray through a vertex is counted exactly once. And test for the boundary first, using orientation zero plus the bounding-box check from the segment test.

def point_in_polygon(p, poly):
    """'boundary', 'inside' or 'outside' using the winding number."""
    wn, n = 0, len(poly)
    for i in range(n):
        a, b = poly[i], poly[(i + 1) % n]
        if orient(a, b, p) == 0 and on_segment(p, a, b):
            return "boundary"
        if a[1] <= p[1]:
            if b[1] > p[1] and orient(a, b, p) > 0:
                wn += 1                   # upward edge, p strictly left
        elif b[1] <= p[1] and orient(a, b, p) < 0:
            wn -= 1                       # downward edge, p strictly right
    return "inside" if wn != 0 else "outside"

On the worked polygon the code reports (3,1) inside, (3,3) outside because it sits in the notch above the reflex vertex, (6,2) and the reflex vertex (3,2) on the boundary, (1,3) inside and (7,1) outside. Parity and winding differ only for self-overlapping shapes, which is the even-odd versus non-zero fill rule in SVG. For many polygons, index bounding boxes first; a structure like the one in the KD-tree article prunes candidates before the exact test runs.

Sorting directions by angle without trigonometry

Many algorithms need directions sorted by angle: Graham scan around a pivot, visibility sweeps. The obvious key, atan2(y, x), is a floating-point approximation, so two collinear directions can compare unequal and nearly equal ones can swap order. The exact alternative splits the plane into two half-planes and compares with the cross product inside each.

from functools import cmp_to_key

def half(v):
    """0 for angles in [0, pi), 1 for [pi, 2*pi)."""
    return 0 if (v[1] > 0 or (v[1] == 0 and v[0] > 0)) else 1

def polar_cmp(u, v):
    if half(u) != half(v):
        return half(u) - half(v)
    c = cross(u, v)
    return -1 if c > 0 else (1 if c < 0 else 0)

def polar_sort(vectors):
    return sorted(vectors, key=cmp_to_key(polar_cmp))

Within a half-plane, directions are under 180 degrees apart, so the cross-product sign is consistent; across halves the index decides. Sorting (1,0), (0,1), (-1,0), (0,-1), (1,1), (-1,-1), (1,-1), (-1,1) and (2,0) produces (1,0), (2,0), (1,1), (0,1), (-1,1), (-1,0), (-1,-1), (0,-1), (1,-1). (1,0) and (2,0) compare equal; break ties by length if needed, and remove the zero vector first.

Integers, overflow and floating-point robustness

Every predicate above is exact on integers, and that is the main reason to keep coordinates as integers whenever the problem allows: grid cells, pixels, fixed-point map coordinates in units of centimetres or 1e-7 degrees. With coordinates bounded by M in absolute value, differences are at most 2M, each product at most 4M squared, and a cross product at most 8M squared. For signed 64-bit integers that means M below 2 to the 30th power (about 1.07e9) is safe. Coordinates up to 1e9 give a worst case of 8e18, under the 9.22e18 limit. Outside Python, check the bound or use 128-bit intermediates.

Floating point breaks predicates in a specific way: orientation can return the wrong sign for nearly collinear points, and different evaluations can disagree with each other. Checked by running it: the doubles (0.1, 0.3), (0.7, 2.1) and (0.3, 0.9), taken at their exact stored binary values, are exactly collinear, yet the floating-point cross product evaluates to 5.55e-17, a confident left turn. Code that gets two answers to one question can loop, spike a hull or crash. An epsilon hides this but makes equality intransitive. The principled options are snapping input to an integer grid, exact rational arithmetic, or adaptive-precision predicates of the Shewchuk kind, which compute fast in floating point and fall back to exact arithmetic only when the result is too close to call.

Failure modes

• Testing only the general-position crossing condition in segment intersection, which silently misses T-junctions and shared endpoints.
• Using a ray-crossing test without a half-open edge rule, so rays through vertices are counted twice or not at all.
• Assuming a vertex order without checking the sign of the shoelace sum; clipping and offsetting algorithms flip their output.
• Sorting by atan2 and then running an algorithm that assumes collinear points sort together.
• Computing an intersection point in floating point and then asking orientation questions about it; the point is not exactly on either line.
• Overflow in a 32-bit cross product: by the same 8M squared bound, int32 is only safe for coordinates below about 16,000 (2 to the 14th).

Trade-offs and related reading

The primitives here are the leaves; the interesting algorithms compose them. Closest pair of points uses squared distances and sorting to reach O(n log n). The smallest enclosing circle uses circumcircle constructions inside a randomised incremental loop. Convex hulls, line sweeps and polygon clipping all reduce to orientation tests. Integers are fast but bound your range, rationals allocate, and adaptive predicates are best borrowed from a library. For production GIS or CAD work, a mature kernel such as CGAL, GEOS or JTS has already paid for the degenerate cases; write your own for small integer problems.

What to do next

1. Decide on your coordinate type first: integers with a documented bound, rationals, or floats plus an exact-predicate library.
2. Implement orient, on_segment and area2, and unit-test the collinear, touching and overlapping cases from this page.
3. Normalise every input polygon to counter-clockwise with area2 and reject self-intersecting input if your algorithms assume simplicity.
4. Use the winding-number test with an explicit boundary result and decide per application what boundary points mean.
5. Replace atan2-based sorting with the half-plane and cross-product comparator wherever an algorithm branches on the order.
6. Write down the 64-bit overflow bound for your data and add an assertion at the input boundary.
7. Fuzz your predicates with nearly collinear random points and check that results stay consistent across repeated and permuted calls.