Signed Area of Polygon, in depth: orientation and winding, edge functions and barycentrics in rasterizers, backface culling, curved outlines and moments

What the sign means: orientation and winding

Green's theorem turns area into a boundary integral: A = 1/2 * closed integral of (x dy - y dx). On a straight edge from p to q the integral is cross(p, q) / 2, which is the signed area of the triangle formed by the origin, p and q. The polygon's signed area is the sum of those fan triangles; the parts outside the polygon are counted once with each sign and cancel. Nothing in that argument needs the polygon to be simple, which is why the formula still says something for self-intersecting and repeated rings.

For any closed ring the sum equals the integral of the winding number over the plane: each region counts its area times the number of times the ring circles it anticlockwise. That explains results that look like bugs.

def signed_area2(pts):
    # Twice the signed area; > 0 for counter-clockwise in a y-up frame.
    s = 0
    for i in range(len(pts)):
        x0, y0 = pts[i]
        x1, y1 = pts[(i + 1) % len(pts)]
        s += x0 * y1 - x1 * y0
    return s

signed_area2([(0, 0), (4, 0), (4, 3), (0, 3)])        #  24: 4 x 3 rectangle, CCW
signed_area2([(0, 3), (4, 3), (4, 0), (0, 0)])        # -24: same rectangle, CW
signed_area2([(0, 0), (2, 2), (2, 0), (0, 2)])        #   0: bow tie, lobes cancel
signed_area2([(0, 0), (1, 0), (1, 1), (0, 1)] * 2)    #   4: unit square traced twice

The bow tie's lobes have winding numbers +1 and -1, so its signed area is zero even though it covers two square units. A ring traced twice has winding number 2 and reports double area. A repeated closing vertex is harmless, because cross(p, p) = 0. So the signed area is an orientation test only for simple rings; for anything else it is a winding-weighted measure, and a validity check must come first, as in Polygon Algorithms, in depth.

Edge functions: the rasterizer's inner loop

Take a triangle a, b, c and a pixel centre p. The edge function E_ab(p) = (b.x - a.x)(p.y - a.y) - (b.y - a.y)(p.x - a.x) is twice the signed area of triangle a, b, p. It is positive when p is to the left of the directed edge. For a counter-clockwise triangle, p is inside exactly when all three edge functions are non-negative. Pineda described this formulation for parallel rasterization in 1988, and it suits hardware: each function is linear in p, so stepping one pixel right adds the constant -(b.y - a.y) and stepping one row adds b.x - a.x, and many pixels can be tested at once.

Ties need a rule. A pixel centre exactly on an edge shared by two triangles must be drawn by exactly one of them, or a mesh shows either cracks or double-blended pixels. The common convention, the top-left rule used by Direct3D and common elsewhere, keeps on-edge pixels only for top edges (horizontal, above the interior) and left edges. With integer coordinates it is just a bias of -1 on the other edges.

def edge(a, b, p):
    return (b[0] - a[0]) * (p[1] - a[1]) - (b[1] - a[1]) * (p[0] - a[0])

def is_top_left(a, b):          # CCW triangle, y up: top edges run right-to-left
    dx, dy = b[0] - a[0], b[1] - a[1]
    return (dy == 0 and dx < 0) or dy < 0

def raster(tri, width, height):
    a, b, c = tri
    if edge(a, b, c) == 0:
        return set()                                  # zero area: draws nothing
    if edge(a, b, c) < 0:
        b, c = c, b                                   # make it counter-clockwise
    edges = ((b, c), (c, a), (a, b))
    bias = [0 if is_top_left(p, q) else -1 for p, q in edges]
    covered = set()
    for y in range(height):
        for x in range(width):
            p = (2 * x + 1, 2 * y + 1)                # pixel centre, doubled coordinates
            if all(edge(e0, e1, p) + k >= 0 for (e0, e1), k in zip(edges, bias)):
                covered.add((x, y))
    return covered

Vertices here are in doubled integer coordinates so pixel centres are integers too, the same trick a GPU plays by snapping vertices to a fixed-point grid with a few subpixel bits. In tests, 200 random quadrilaterals split along a diagonal never had a pixel claimed by both halves, and an eight-triangle fan around an interior point covered each of the 576 pixel centres inside its outline exactly once, with no gaps.

Barycentric coordinates are area ratios

Divide each edge function by the triangle's own doubled area and you get barycentric coordinates: l0 = E_bc(p) / 2A, l1 = E_ca(p) / 2A, l2 = E_ab(p) / 2A. They sum to one, they are all non-negative inside, and p = l0 * a + l1 * b + l2 * c. For a = (0, 0), b = (8, 0), c = (0, 6) and p = (2, 3) they come out as 0.25, 0.25 and 0.5, and the weighted sum gives back (2, 3). A rasterizer computes these weights anyway, so interpolating colours, texture coordinates and depth across the triangle is almost free. Under perspective, attributes are interpolated as attribute / w together with 1 / w and divided per pixel, but the area ratios are still the starting point.

The same signed areas drive the orientation tests used to walk through a mesh to the triangle containing a point, as in Delaunay triangulation.

Backface culling and the y-axis trap

A closed mesh viewed from outside shows each triangle's front only when its projected vertices run one way, so the sign of the projected signed area says whether a triangle faces the camera. OpenGL computes the area from window coordinates; glFrontFace chooses whether counter-clockwise (GL_CCW, the default) or clockwise triangles are front-facing, and glCullFace with GL_CULL_FACE enabled decides which side is discarded. The spec also negates the sign when glClipControl sets the clip origin to GL_UPPER_LEFT.

That last detail is the general lesson. Flipping the y axis flips every sign. A mesh exported with one winding convention, drawn through a projection or framebuffer that flips y, culls all of its front faces and the screen goes blank. When that happens, check the winding convention and the axis direction before anything else. Zero-area triangles are discarded by the same test, which is why degenerate slivers simply vanish.

Curved boundaries: Green's theorem per segment

Green's theorem works segment by segment, so the area of a shape bounded by curves is the sum of each segment's 1/2 * integral of (x dy - y dx). For a quadratic Bezier with control points P0, P1, P2 the integral has a closed form, (2 cross(P0, P1) + 2 cross(P1, P2) + cross(P0, P2)) / 6. For a cubic, the integrand is a polynomial of degree 5 in the parameter, and three-point Gauss-Legendre quadrature is exact for degree 5, so three evaluations give the exact term.

GL3 = [(-(3 / 5) ** 0.5, 5 / 9), (0.0, 8 / 9), ((3 / 5) ** 0.5, 5 / 9)]

def curve_area_term(ctrl):
    # Exact for Bezier segments up to degree 3; ctrl is a list of (x, y).
    s = 0.0
    for xi, wt in GL3:
        t = (xi + 1) / 2
        x, y = bezier(ctrl, t)
        dx, dy = bezier_derivative(ctrl, t)
        s += wt * (x * dy - y * dx) / 4    # 1/2 for the area, 1/2 for the [0, 1] map
    return s

Checked against 200,000-segment polylines: the arch (0,0), (2,4), (4,0) closed by its chord gives -5.3333, which is minus two thirds of base times height (clockwise, hence negative), and a test cubic gives -6.7000 by both methods. A circle of radius one drawn with four cubic arcs using the usual control distance 4(sqrt 2 - 1)/3 has area 3.14247, which is 0.028 percent above pi. That error belongs to the curve approximation, not to the integration. Font outlines are made of such segments, and formats disagree on which direction outer contours run, so check your format's specification instead of assuming a sign.

Moments from the same sum

The same edge-by-edge decomposition gives higher moments, each weighted by the edge's cross product k = x0*y1 - x1*y0. The centroid is Cx = sum((x0 + x1) k) / (6A) and likewise for y. The second moment about the x axis is Ixx = sum((y0^2 + y0*y1 + y1^2) k) / 12, with Iyy the same in x. Because every term is signed, reversing the ring negates the area and the moments while the centroid stays put, so compute with signed values and take one absolute value at the end.

Worked example: the L-shape (0,0), (6,0), (6,2), (2,2), (2,5), (0,5). The formulas give area 18, centroid (7/3, 11/6), Ixx = 94 and Iyy = 152 about the origin. A 600 by 600 midpoint grid over the shape agreed to within 0.0002 on every quantity. Reversed, the area reads -18 and the centroid is unchanged. Section properties for beams, the inertia tensor of a 2D rigid body and the orientation of a blob in image analysis all come from these sums.

Failure modes

• Cancellation far from the origin. The triangle (1e8, 1e8), (1e8 + 1, 1e8), (1e8, 1e8 + 1) has area 0.5, but the shoelace sum in doubles returns 0.0, because each product is about 10¹⁶, where doubles are spaced 2 apart. Subtracting the first vertex from every point first gives 0.5.
• Overflow in fixed point. Edge functions multiply two coordinate differences, so 16-bit coordinates need a product of over 32 bits. Size the accumulator for the product, not the inputs.
• Axis convention. Screen coordinates usually run y down, which turns counter-clockwise into clockwise. Every orientation-dependent decision should name its frame.
• Sign near zero. Nearly collinear points can produce the wrong sign in floating point; use an exact or adaptive predicate when a wrong sign changes topology.
• Invalid rings treated as valid. A bow tie reports zero and a doubled ring reports double; neither raises an error. Validate before trusting a sign.
• Unclosed curved paths. Green's theorem needs a closed boundary; an open path's sum depends on where the origin is.

Trade-offs

For the related hull and ordering tasks where the same cross product is the only primitive, see Andrew's Monotone Chain, in depth.

What to do next

1. Write signed_area2 and run it on a counter-clockwise square, its reverse, a bow tie and a doubled ring; explain each result.
2. Pick and document the winding convention and axis direction for your codebase.
3. Implement the edge-function rasterizer with the top-left rule and test a split quad for double coverage.
4. Compute barycentric weights for a few points and reconstruct the points from them.
5. Add the curved-segment term and compare against a dense polyline on your own outlines.
6. Compute centroid and second moments for an L-shape and check them against a grid.
7. Shift coordinates to a local origin before summing, and switch to a robust predicate wherever a sign changes topology.