Pick’s theorem

How can the area of a polygon be found by counting points? For a polygon whose corners lie on a square grid of points, Pick’s theorem gives the area as \(A = I + B/2 - 1\), from the number \(I\) of grid points inside the polygon and the number \(B\) of grid points on its boundary. No lengths, angles or coordinates enter the formula. For a polygon with \(h\) holes, the area is \(A = I + B/2 + h - 1\).

The example draws a polygon on the grid, counts the grid points with Contains and compares the counts with the polygon’s exact area, first as it is and then with two holes cut out.

import math

import numpy as np
import planex as px

Counting the grid points

Contains tests an array of points against a face at once. With boundary=False, it leaves out the points on the boundary, so the difference between the two counts is \(B\):

def grid_points(face):
    """the grid points inside the face and those on its boundary"""
    (x0, y0), (x1, y1) = face.bounding_box
    grid = np.array([(x, y) for x in range(math.floor(x0), math.ceil(x1) + 1)
                     for y in range(math.floor(y0), math.ceil(y1) + 1)], float)
    inside = face.Contains(grid, boundary=False)
    return grid[inside], grid[face.Contains(grid) & ~inside]

A polygon on the grid

The example is a polygon with eight corners on the grid. Its area is compared with \(I + B/2 - 1\) (Figure 1):

polygon = px.FromPolyline([(0, 0), (8, 1), (11, 6), (7, 5), (6, 10), (2, 7), (-2, 9), (1, 4)], closed=True).Face()
inside, boundary = grid_points(polygon)
I, B = len(inside), len(boundary)
print("area", polygon.mass, "| I", I, "B", B, "| I + B/2 - 1", I + B / 2 - 1)

area 63.5 | I 60 B 9 | I + B/2 - 1 63.5

Figure 1: The grid points inside the polygon (dots) and on its boundary (crosses)

Holes

Each hole adds 1 to the formula. With two triangles cut out, both with their corners on the grid, the polygon has two holes, and the grid points on their edges count toward \(B\). The area is compared with \(I + B/2 + h - 1\) for \(h = 2\) (Figure 2):

holes = [px.FromPolyline([(3, 2), (6, 2), (5, 4)], closed=True).Face(),
         px.FromPolyline([(2, 5), (4, 6), (3, 7)], closed=True).Face()]
holed = polygon - holes[0] - holes[1]
inside, boundary = grid_points(holed)
I, B = len(inside), len(boundary)
print("area", holed.mass, "| I", I, "B", B, "| I + B/2 + h - 1", I + B / 2 + len(holes) - 1)

area 59.0 | I 49 B 18 | I + B/2 + h - 1 59.0

Figure 2: The polygon with two holes and its grid points

See Inside or outside for Contains.