I worked with a friend on an university project to do with mobile-phone application development, for which we wanted an algorithm to determine the country that had given GPS coordinates. Here we describe our algorithm one step at a time and explain how it works.
1. We need a non-self-intersecting polygon with vertices
, indexed to go around the polygon in either direction.
2. Regroup the vertices as
3. For each pair of vertices
, define a linear function
that passes through both vertices:
4. For each pair of vertices
, compute the midpoint
The gradient vectors
in this case are constant vectors, either all pointing to the interior of the polygon or all pointing to the exterior, depending on whether the polygon vertices were indexed counterclockwise or clockwise. This makes sense because we can treat a polygon as a discretized, closed contour curve of a function
, and a closed contour implies the presence of a local maximum or minimum. From experience, the vectors
all point into the interior if the vertices are indexed counterclockwise.
5. Find the bounding box of the polygon, a rectangular region identified by the highest, lowest, left-most, and right-most coordinates. Given coordinates
to test, if
is outside the bounding box then
is outside the polygon. If
is inside the bounding box, go to Step 6. (This step can alternatively be done right at the beginning.)
6. Compute the discriminants
for each pair of vertices.
7. Find the edge closest to
; if the edge closest to
is the one defined by
, then the distance between
will be the minimum. Find
8a. Now suppose that all
are pointing to the interior (counterclockwise-defined polygon): then
is inside if
; otherwise it is outside.
8b. If on the other hand all
are pointing to the exterior (for the clockwise orientation), then
is inside if
; otherwise it is outside.
Because the polygon can be nonconvex, a ray drawn from a point may cross more than one edge, but only the edge encountered first is relevant. Finding the edge closest to a point ensures that we are dealing with such an edge.
We can tell from the part of the boundary closest to the test point whether it is inside or outside because we know the
are either all pointing to the interior of the polygon or all pointing to the exterior.
8c. Finally, it may happen that
. In that case, perturb
by a very small amount in a random direction and go back to Step 7.
If the point is well inside or well outside the polygon, the result will be consistent no matter how we perturb
. The result will be inconsistent if
is close to the edge. Since we only perturb
by a very small amount, we can infer that
is "on the edge", for which we may still decide that
is in the interior.
From Step 4 and onwards, instead of calculating
and computing the closest edge, we draw a ray from the testing coordinate and count the number of edges intersected by the ray. If the ray intersects the polygon an even number of times (including zero) then the testing coordinate is outside the polygon; otherwise it is inside. This method is also known as the odd-even rule. In both cases the problem to be solved here is how to deal with the nonconvex part of the polygon.
Our goal is an algorithm that is able to handle the scenario where a country is inside a bigger country, for example Lesotho and South Africa. Suppose that the GPS location is in Lesotho: we want an algorithm to make an unambiguous decision. In this case one solution would be to find the distances between the current point and the closest edge for each polygon, and then take the polygon (country) that gives the minimum. This is why we initially came up with the alternative algorithm that involves finding the closest polygon edge and its geometric properties. However, the odd-even rule is much easier to implement and less computationally intensive, although in practice we still need to deal with the special cases in which the testing coordinate coincides with one of the vertices or the edge separately.