Points, vectors and the cross product

Computational geometry solves problems about points, lines and shapes: maps, games, robots, chip layouts. Almost all of it rests on one small formula, the cross product, which tells you whether a turn goes left or right and how much area two steps enclose. With whole-number coordinates it is exact, and that matters more here than anywhere else in the course.

Points and vectors

A point is a place, given by two coordinates (x, y): x goes right and y goes up. A vector is a step: how far to move in x and in y. The step from point A to point B is B − A = (B.x − A.x, B.y − A.y). In Python a pair (x, y) works for both.

• Adding a vector to a point moves the point: A + (B − A) = B.
• Two vectors add coordinate by coordinate; multiplying by a number stretches a vector.
• The length of v = (x, y) is √(x² + y²), written math.hypot(x, y). To compare lengths, compare x² + y² instead: no square root, so the answer stays a whole number (DVCQ-4 did this for distances).

The dot product

The dot product of a = (a.x, a.y) and b = (b.x, b.y) is a·b = a.x·b.x + a.y·b.y. It equals |a| · |b| · cos θ, where |a| is the length of a, θ (theta) is the angle between the two vectors, and cos θ is the cosine of that angle. Only its sign is needed most of the time:

• positive: the angle is less than 90°, so the vectors point roughly the same way;
• zero: they are at right angles;
• negative: the angle is more than 90°.

For example, the point of a line piece AB that is closest to a point P: it lies beyond B if (P − B)·(B − A) > 0, before A if (P − A)·(B − A) < 0, and between them otherwise.

The cross product

The cross product of a and b, for vectors in the plane, is the single number

• Its size is the area of the parallelogram with sides a and b (|a| · |b| · sin θ, with sin θ the sine of the angle).
• Its sign is the direction of the turn from a to b. Positive: counterclockwise, so b points to the left of a. Negative: clockwise, to the right. Zero: a and b lie along one line (or one of them is zero).

Why the area? Start with a = (a.x, 0) along the x-axis. The parallelogram has base a.x and height b.y, so its area is a.x · b.y, which is the formula with a.y = 0. Turning both vectors by the same angle changes neither the area nor the formula’s value, so the formula holds for every a. The sign follows the same way: with a along the x-axis, b points left of a exactly when b.y > 0.

Swapping the two vectors flips the sign: cross(b, a) = −cross(a, b). For three points O, A, B, write cross(O, A, B) for cross(A − O, B − O). Step through the figure.

def cross(o, a, b):
    """> 0: o -> a -> b turns left;  < 0: right;  0: one line."""
    return (a[0] - o[0]) * (b[1] - o[1]) - (a[1] - o[1]) * (b[0] - o[0])

Whole numbers, not floats

Geometry programs often decide something from a sign: left or right, inside or outside. With floats a value that should be 0 can come out as 1e-17, or a tiny positive value can round to 0, and then the program takes the wrong branch. In 0.1 + 0.2 == 0.3 Python answers False.

• Keep coordinates as whole numbers when the input allows it. The cross product of whole numbers is a whole number, and Python’s whole numbers are exact at any size: with coordinates up to 10⁹ the products reach 10¹⁸, with no rounding at all.
• When division is needed (an intersection point), use fractions.Fraction, which keeps an exact numerator and denominator (lesson 2).
• When floats cannot be avoided (angles, square roots), compare with a small tolerance, such as abs(v) < 1e-9, and know that the tolerance is a guess.

A float has about 16 significant digits. For coordinates near 10¹⁷, 10¹⁷ · 10¹⁷ and (10¹⁷ + 1)(10¹⁷ − 1) are both stored as 1e34, although they differ by 1. Lesson 2’s practice problem has such a test.

The area of a polygon

A polygon is given by its corners P₀, P₁, …, P_n−1 in order around the edge; it is simple if its edges meet only at shared corners. Split it into the triangles P₀, P_k, P_k+1 and add their cross products. Triangles that cover a spot outside the polygon come in with a minus sign and cancel out (the fourth variant of the figure shows one). Using the point (0, 0) instead of P₀ gives the same total and a shorter formula, the shoelace formula. The bars | … | mean the value without its sign:

2 · area = | (x₀y₁ − x₁y₀) + (x₁y₂ − x₂y₁) + … + (x_n−1y₀ − x₀y_n−1) |

def twice_area(poly):
    s = 0
    for k in range(len(poly)):
        x1, y1 = poly[k]
        x2, y2 = poly[(k + 1) % len(poly)]   # the last edge closes the loop
        s += x1 * y2 - x2 * y1
    return s          # > 0 counterclockwise, < 0 clockwise

The sign tells the direction: positive when the corners go counterclockwise. Many algorithms need corners in a known direction, and this is how to check. It costs O(n). Keep the doubled area: it is always a whole number, while the area may end in .5.

Points inside: Pick’s theorem

For a simple polygon whose corners have whole-number coordinates, let A be its area, B the number of whole-number points on its edges, and I the number strictly inside. Georg Pick (1899) proved:

B is easy to count. An edge from (x₁, y₁) to (x₂, y₂) passes through gcd(|x₂ − x₁|, |y₂ − y₁|) + 1 whole-number points, counting both ends (gcd is the greatest common divisor, MATH-1). The points on it are evenly spaced, every 1/g of the way, where g is that gcd. Counting each corner once, B is the sum of the gcds over all edges.

Check it on the square with corners (0, 0) and (4, 4): 2A = 32, B = 16, so I = (32 − 16 + 2)/2 = 9, the points (1..3, 1..3).

Now try GEOM-P1, Points inside a field.