Computational Geometry

Geometry problems show up in every ICPC contest, and in a lot of online contests. Knowing how to solve geometry problems is a staple of any good competitive programmer's arsenal. They are often considered to be the hardest category of problems, as they are filled with corner cases, and rounding errors can be deadly. Today we will cover the fundamental techniques needed to solve basic geometry problems.

Notebook code

Due to the ubiquity of corner cases and rounding errors in geometry problems, it is highly recommended to have prewritten, pretested code to use to solve standard geometry problems. As a beginner, as you solve problems, you should save your reliable code to reuse for the future.

Pitfalls -- Overflow and rounding errors

If you are using integer values, watch out for overflow! If you are using floating-point values, watch out for rounding errors! Here are some useful tips for both of these. When comparing integer-valued equations for geometry problems, it is often necessary to compare two fractions, say \( a/b \) and \( c/d \). Do not simply convert these to doubles, or you may get rounding errors! Instead, note that the following works: $$ \frac{a}{b} < \frac{c}{d} \iff ad < bc $$ You can compare \( ad \) and \( bc \) without converting to doubles. Watch our for overflow, though, since you need to multiply.

Another useful trick to is realize that some computations that look like they require floating-point computations do not. A common example is comparing the magnitude of two point vectors. Recall that the magnitude (length) of a vector \( (x,y) \) is \( \sqrt{x^2 + y^2} \). Taking the square root will require using floating-point numbers. However, notice the following useful fact: $$ \sqrt{{x_1}^2 + {y_1}^2} < \sqrt{{x_2}^2 + {y_2}^2} \iff {x_1}^2 + {y_1}^2 < {x_2}^2 + {y_2}^2 $$ In other words, if you have to compare two square roots, you can simply avoid taking the square root and compare them anyway. This means that the entire computation can stay in integers.

If you need to compare floating-point values, you usually want to use some small tolerance since rounding errors can make two values that should be equal look not equal. So instead of checking whether \( a = b \), you usually want to check whether \( |a-b| < \varepsilon \), for some small \( \varepsilon \), usually something like \( \varepsilon = 10^{-8} \), depending on the precicison of the problem.

Primitives

The simplest objects that you will deal with in geometry problems are points and lines. Most commonly, problems will be in 2D, so points will have two coordinates, but occassionally you might be asked to solve 3D problems. For today, we will just cover 2D. Many of the ideas carry over the 3D, but some do not and become much harder.

Points

A point will either be represented by a pair of integers, or a pair of doubles, depending on whether you need floating-point values or whole numbers. There are many options to do this. One is to just write a custom Point class, and add methods for performing your favourite geometric operations on them. You can see an example here along with implementations of many of the operations we will talk about.

Another option if you use C++ is to make use of C++'s complex number class. Operations on complex numbers often mirror common geometry operations; for example, lengths, angles, dot products, and cross products can be computed in terms of complex numbers. If you're interested in this approach, read this.

Lines

A line can be defined by a pair of points. There are three kinds of lines that you might need to consider. Line segments are finite-length lines with two endpoints. A ray is an infinite line that starts at a point and goes off forever in some direction. A line usually refers to an infinite line that goes off forever in both directions. Segments are usually defined by their endpoints. Rays are defined by their endpoint and direction, and lines are defined by any pair of points that sit on them.

Polygons

A polygon is a sequence of at least three points defining some closed shape. A polygon is convex if for all pairs of points inside it, the line segment between them does not leave the polygon. Many problems are much easier to solve on convex polygons than on general polygons.

Common operations

Length (magnitude) of a vector

The length (or magnitude) of a vector \( v = (x,y) \) is defined to be $$ |v| = \sqrt{x^2 + y^2} $$

Dot product and cross product

Dot products and cross products are useful for computing many things, especially angles and areas. Given two vectors \( v_1 = (x_1, y_1) \) and \( v_2 = (x_2, y_2) \), their dot product is defined as $$ v_1 \cdot v_2 = x_1 x_2 + y_1 y_2 $$ A useful fact about dot products is the following relationship: $$ v_1 \cdot v_2 = |v_1| |v_2| \cos(\theta), $$ where \( \theta \) is the angle between the two vectors. Since we know how to compute the lengths of vectors, this means that we can compute the angle between two vectors with the following formula: $$ \theta = \cos^{-1}\left(\frac{v_1 \cdot v_2}{|v_1| |v_2|}\right) $$

The cross product of two vectors \( v_1, v_2 \) is defined as $$ v_1 \times v_2 = x_1 y_2 - x_2 y_1 $$ A useful fact is that the cross product gives you (double) the area of the triangle formed by the two vectors. $$ \text{Area} = \frac{1}{2} v_1 \times v_2 $$ Note that this assumes that the angle between \( v_1 \) and \( v_2 \) is positive. Otherwise you will get a negative area.

Orientation

Given three points \( a, b, c \), the orientation test tells us whether the path from \( a \to b \to c \) turns clockwise or anticlockwise.

Orientation can be computed using the sign of the cross product of \( (b-a) \times (c-b) \). Specifically, if this is positive, then the orientiation is anticlockwise. If it is negative, it is clockwise, and otherwise if it is zero, then the lines are parallel.

Useful algorithms

Polygon area

The area of a polygon can be computed by using the fact that cross products tell us the area of triangles, and then breaking the polygon into triangles. Consider a convex polygon centered at the origin. We can see that we can break it into triangles, and then use the cross product formula to compute the area as $$ \sum_{i=0}^{n-1} \frac{1}{2} p_i \times p_{i+1} $$ where \( p_i \) denotes the \( i^\text{th} \) point in the polygon in counterclockwise order (this is important. If we reverse the order we will get a negative area!)

Now, here's the magic! This formula actually works for all polygons, even if they are nonconvex and even if they are not centered at the origin! Try to figure out why this works. Hint: Remember that you get a negative area if the order of the sides is reversed.

Convex hull

The convex hull is probably the most well-known computational geometry algorithm, and a very useful one for programming competitions. Given a set of points, their convex hull is the smallest convex polygon that contains all of them. Note that the convex hull is a polygon whose points are from the input points -- it will never contain a vertex not from the input.

Understanding how the convex hull algorithm works is nice, but as a beginner competitive programmer, it is more important to understand how to use it. So we will only talk about how the algorithm works briefly. There are many algorithms for convex hulls; we will look at the Graham scan algorithm.

We can first sort all of the points by their \( x \) coordinate. We then simultaneously construct the top half and bottom half of the convex hull (called the upper hull and lower hull). The key observation is that for a convex polygon, the orientation of every pair of adjacent sides is the same, so we can test every new point to see whethe the orientation changes, and if it does, it should not be part of the hull.

// Compute the convex hull of the points P
// Assumes there are no duplicate points in P
function convex_hull(P) {
  sort(P)  // Sort by x value
  Upper = copy(P)  // Upper hull
  Lower = copy(P)  // Lower hull
  for (each point p in P) {
    while (Lower is not empty and cross(p - Lower[-1], Lower[-2] - p) < 0)
      Lower.pop()
    while (Upper is not empty and cross(p - Upper[-1], Upper[-2] - p) > 0)
      Upper.pop()
    Lower.push(p)
    Upper.push(p)
  }
  return Upper + Lower  // Make sure to remove duplicates
}

Here, Lower[-1] means the last element of Lower, and Lower[-2] means the second last, etc. Although the pseudocode makes the algorithm seem quite simple, I strongly recommend that as a beginner, you use someone else's prewritten and thoroughly tested onvex hull code until you are comfortable writing your own. A real implementation has to deal with rounding errors, which are actually pretty nontrivial, so it is harder to implement than the pseudocode suggests.

You can find an example implementation here.


Danny Sleator
Last modified: Mon Apr 6 19:53:26 2020