The point location problem

Computational geometry is a field that studies efficient solution of geometric problems, which are critical in mapping, manufacturing, and particularly graphics applications. If you find data structures and algorithms interesting, it's likely you'd also find computational geometry interesting, because it often combines and adapts well-known data structures and algorithms in novel ways to solve natural problems that are easy to explain. As an example of this, I'm going to discuss the point location problem.

Imagine you have a political map of Africa with a red dot somewhere on it. Someone asks you, what country is that dot in? This is a simple question for even a child to answer, but how do we make a computer answer the same question? We can represent the map as a planar decomposition, meaning a division of the plane into polygonal regions. We store the vertices making up the boundary of each region. Stop and think for a little bit about how you might solve this.

Perhaps the simplest approach is to exploit the cross product. Suppose the vertices of the region are listed in clockwise (or counterclockwise) order. Compute the cross product of the vector from vertex k to vertex k − 1 with the vector from vertex k to the point being located. The point lies in the region if and only if all of these cross products (for all k) point in the same direction. Since we visit each vertex twice and the cross product takes constant time, this requires O(n) time, where n is the total number of edges.

Another simple solution that generalizes better to higher dimensions is to start at the point and draw a line emanating from it. The first time you hit an edge, you know that edge is on the boundary of the correct region. The edge may be on the boundary of two adjacent regions, but comparing the slope of the line to the slope of the edge will tell you which one is correct. To actually compute this, we define the line parametrically with an equation like this, where (x₀, y₀) is the point being located:

(x, y) = (x₀, y₀) + t(1, 1)

Then for each edge, we compute the t value at which the line intersects that edge, if any, and choose the edge with the smallest positive t value. This algorithm is also O(n).

These algorithms are fine if you're only locating one point, but if we wish to locate many points we can do better by creating a data structure which facilitates quick point location. I'll begin by solving a much simpler problem: say your map only has rectangles in it, and these rectangles are arranged in vertical "strips" as shown to the right. To solve this problem, we can construct a balanced binary search tree. This tree will first compare the x coordinate of the point against the x coordinates of the vertical lines to determine what strip of rectangles the point falls in. It then compares the y coordinate of the point against the y coordinates of the horizontal lines to determine what rectangle it falls in. Each comparison eliminates about half of the remaining rectangles from consideration, so this takes only O(log n) time in all for n rectangles. The search tree takes O(nlog n) time to construct and takes O(n) space, since there's at most one node for each edge and there are 4n edges.

We generalize this by allowing each "strip" to contain not just rectangles but trapezoids , as shown below (both rectangles and triangles are considered to be trapezoids here). This is called a trapezoidal map. This makes the second part of the search only a little more complicated: instead of comparing the y coordinate to a fixed value, we choose an edge and use a cross product to check which side of the edge the point falls on.

Finally, take any map, and draw a vertical line through each vertex; in the image to the left the gray vertical lines were added in this manner. You'll end up with a trapezoidal map. Since we can determine what region of this map contains the point in O(log n) time, and each of these regions is part of only one region in the original map, we've solved the problem for general maps in O(log n) time.

Unfortunately, adding all these vertical lines also creates a lot of new regions, potentially as much as squaring the number of regions. Search time is still good, since log(n²) = 2 log n, but the storage required can jump from O(n) up to O(n²). This is worst-case, but even in typical practical cases it can require a lot of space, and construction time for the data structure also increases. We want to retain a quick search capability but without having quite so many regions.

To achieve this, notice that we don't actually need to draw a complete vertical line through every point to divide the map into trapezoids. All we really have to do is draw a line up and down from each vertex until we reach the edges above and below. It can be shown that the number of trapezoids is now no more than about three times the number of original regions, so the worst-case space has been cut back down to O(n). The search procedure is similar: at each node we either compare the x coordinate of the point to that of a vertical line segment, or we check on which side of an edge the point lies. Either way we eliminate about half the remaining regions. The search "tree" is no longer a tree but a directed acyclic graph, because a single trapezoid may fall on both sides of a vertical line segment. See a demonstration.

Of course I've left a lot of details out here, like how to ensure the search tree remains well-balanced during construction, what to do if some points have the same x coordinate, and so on, but I hope this gives you a taste of how clever algorithms can efficiently solve geometric problems.