Intersections, Distances, and Point-in-Polygon Queries

Everything Here Is Built on One Predicate

Part 1 introduced orientation(a, b, c) — the sign of a cross product — as the primitive almost every algorithm in this series reduces to. This part is the proof of that claim: segment intersection, point-to-segment distance, polygon area, and point-in-polygon tests are all short compositions of that one function, plus a handful of boundary-case policies.

Segment Intersection

Two segments $p_1p_2$ and $p_3p_4$ properly intersect (cross at a single interior point of both) exactly when $p_3$ and $p_4$ are on opposite sides of line $p_1p_2$, and $p_1$ and $p_2$ are on opposite sides of line $p_3p_4$ — four orientation calls, no coordinate solving required:

The return type matters. A Boolean answers only “do these sets share anything?” but discards whether the result is empty, one point, or an overlapping segment. Geometry libraries often return NONE, POINT(p), or OVERLAP(a, b). That richer contract becomes essential in subdivision overlay, where an overlapping edge cannot be represented honestly as a single point.

def orientation(a, b, c):
    return (b[0] - a[0]) * (c[1] - a[1]) - (b[1] - a[1]) * (c[0] - a[0])


def on_segment(a, b, p):
    """True if p is collinear with a,b AND within the segment's 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, p3, p4):
    d1 = orientation(p3, p4, p1)
    d2 = orientation(p3, p4, p2)
    d3 = orientation(p1, p2, p3)
    d4 = orientation(p1, p2, p4)

    if ((d1 > 0) != (d2 > 0)) and ((d3 > 0) != (d4 > 0)) and d1 != 0 and d2 != 0 and d3 != 0 and d4 != 0:
        return True  # proper intersection

    # Degenerate cases: an endpoint lies exactly on the other segment.
    if d1 == 0 and on_segment(p3, p4, p1):
        return True
    if d2 == 0 and on_segment(p3, p4, p2):
        return True
    if d3 == 0 and on_segment(p1, p2, p3):
        return True
    if d4 == 0 and on_segment(p1, p2, p4):
        return True
    return False


print(segments_intersect((0, 0), (4, 4), (0, 4), (4, 0)))  # True (proper cross)
print(segments_intersect((0, 0), (1, 1), (2, 2), (3, 3)))  # False (collinear, disjoint)

#include <algorithm>
#include <iostream>

struct Point { double x, y; };

double orientation(const Point& a, const Point& b, const Point& c) {
    return (b.x - a.x) * (c.y - a.y) - (b.y - a.y) * (c.x - a.x);
}

bool onSegment(const Point& a, const Point& b, const Point& p) {
    return p.x >= std::min(a.x, b.x) && p.x <= std::max(a.x, b.x) &&
           p.y >= std::min(a.y, b.y) && p.y <= std::max(a.y, b.y);
}

bool segmentsIntersect(Point p1, Point p2, Point p3, Point p4) {
    double d1 = orientation(p3, p4, p1);
    double d2 = orientation(p3, p4, p2);
    double d3 = orientation(p1, p2, p3);
    double d4 = orientation(p1, p2, p4);

    if (((d1 > 0) != (d2 > 0)) && ((d3 > 0) != (d4 > 0)) &&
        d1 != 0 && d2 != 0 && d3 != 0 && d4 != 0) {
        return true;
    }
    if (d1 == 0 && onSegment(p3, p4, p1)) return true;
    if (d2 == 0 && onSegment(p3, p4, p2)) return true;
    if (d3 == 0 && onSegment(p1, p2, p3)) return true;
    if (d4 == 0 && onSegment(p1, p2, p4)) return true;
    return false;
}

Point-to-Segment Distance

Project the query point $p$ onto the line through $a, b$, clamp the projection parameter $t$ to $[0, 1]$ so the closest point stays on the segment (not the infinite line), then measure the distance to that clamped point:

$$ t = \operatorname{clamp}\!\left(\frac{(p - a) \cdot (b - a)}{(b - a)\cdot(b - a)},\; 0,\; 1\right) \qquad \text{closest point} = a + t (b - a) $$

The clamp is the entire boundary-case story here: without it, a point whose perpendicular foot falls outside the segment would report a distance to a point that isn't actually on the segment.

There are three intuitive regions: before $a$ gives $t<0$, beside the segment gives $0\le t\le1$, and beyond $b$ gives $t>1$. Clamping collapses the first and third regions onto the appropriate endpoint.

from math import hypot


def point_segment_distance(p, a, b):
    ab = (b[0] - a[0], b[1] - a[1])
    ap = (p[0] - a[0], p[1] - a[1])
    length_squared = ab[0] * ab[0] + ab[1] * ab[1]
    if length_squared == 0:
        return hypot(ap[0], ap[1])
    t = (ap[0] * ab[0] + ap[1] * ab[1]) / length_squared
    t = max(0.0, min(1.0, t))
    closest = (a[0] + t * ab[0], a[1] + t * ab[1])
    return hypot(p[0] - closest[0], p[1] - closest[1])

The zero-length check is not cosmetic: duplicate endpoints make the projection denominator zero. Treating the “segment” as a point gives a meaningful answer.

Polygon Signed Area and Point-in-Polygon

The shoelace formula gives a simple polygon's signed area directly from its vertices $p_0, \dots, p_{n-1}$:

$$ A = \frac{1}{2} \sum_{i=0}^{n-1} \left( x_i y_{i+1} - x_{i+1} y_i \right) $$

The sign of $A$ tells you the polygon's orientation (positive for counterclockwise, negative for clockwise) — the same orientation test from Part 1, just summed over every edge instead of one triangle. The full derivation, plus the winding number and the closely related ray-crossing point-in-polygon test, is the subject of the "Polygon Area and Winding Number" algorithm spotlight later in this series; here it's enough to state the ray-crossing rule this part's code uses: cast a ray from the query point to infinity and count how many polygon edges it crosses — odd means inside, even means outside.

Parity works because every ordinary boundary crossing toggles the state from outside to inside or back again. A ray through a vertex is awkward: counting both incident edges can toggle twice at one location. The standard half-open rule counts an edge only when one endpoint is strictly above the query's horizontal line and the other is at or below it. The unequal > tests in the code encode that convention.

def point_in_polygon(point, polygon):
    """Ray-crossing test. polygon is a list of (x, y) vertices, no duplicate closing point."""
    x, y = point
    inside = False
    n = len(polygon)
    j = n - 1
    for i in range(n):
        xi, yi = polygon[i]
        xj, yj = polygon[j]
        if ((yi > y) != (yj > y)) and (x < (xj - xi) * (y - yi) / (yj - yi) + xi):
            inside = not inside
        j = i
    return inside


square = [(0, 0), (4, 0), (4, 4), (0, 4)]
print(point_in_polygon((2, 2), square))   # True
print(point_in_polygon((5, 2), square))   # False

#include <vector>

struct Point { double x, y; };

bool pointInPolygon(const Point& p, const std::vector<Point>& poly) {
    bool inside = false;
    int n = poly.size();
    for (int i = 0, j = n - 1; i < n; j = i++) {
        bool crosses = (poly[i].y > p.y) != (poly[j].y > p.y);
        double xIntersect = (poly[j].x - poly[i].x) * (p.y - poly[i].y) /
                             (poly[j].y - poly[i].y) + poly[i].x;
        if (crosses && p.x < xIntersect) inside = !inside;
    }
    return inside;
}

Circle-Line and Circle-Circle Intersection

Circle-line intersection reduces to substituting the line's parametric form into the circle equation and solving a quadratic in the parameter — two real roots (two crossing points), one repeated root (tangent), or no real roots (no intersection). Circle-circle intersection is the 2,000-year-old Problem of Apollonius: given two circles with centers $c_1, c_2$ and radii $r_1, r_2$, find their intersection points by first computing $d = |c_2 - c_1|$, then checking the triangle inequality on $r_1, r_2, d$ (no intersection if $d > r_1 + r_2$ or $d < |r_1 - r_2|$) before solving for the intersection chord.

The center distance $d$ classifies the relationship before any square root is taken. Besides the pictured cases, $d<|r_1-r_2|$ means one circle lies strictly inside the other, while $d=0$ with equal radii means the circles are coincident. Classify first, construct coordinates second: this avoids taking the square root of a tiny negative value created by rounding near tangency.

Apollonius of Perga posed and solved the general tangent-circles problem (construct a circle tangent to three given circles, points, or lines) around 200 BCE, using only compass and straightedge; the special case here — where the "circles" degenerate to the specific configuration of finding intersection points rather than a tangent circle — is the version every 2D physics engine and CAD kernel implements today, usually without realizing its lineage.

Boundary Cases This Part Commits To

• Tangent circles ($d$ exactly equal to $r_1 + r_2$ or $|r_1 - r_2|$): exactly one intersection point, not zero and not two — worth an explicit check rather than relying on the quadratic's discriminant landing on exactly zero in floating point.
• Coincident circles (same center and radius): infinitely many intersection points; this part's implementations treat this as a distinct degenerate case rather than "two intersection points."
• Point-on-boundary for point-in-polygon: as shown above, left as an explicit implementation choice rather than a guaranteed answer.

A dedicated "Polygon Area and Winding Number" algorithm spotlight, later in the series, covers the winding-number test and its relationship to ray-crossing in full depth. The next part builds the first real algorithm on top of these primitives: computing the convex hull of a point set.