Point-In-Polygon Testing Part 1

One of the most basic geometric tests is if a point, p, is inside a polygon poly.

The most common method for this test is casting a ray from the point p either horizontally or vertically and checking how many times it intersects the polygon. If the number of intersections is odd then p is inside poly. If the number of intersections is even then p is outside poly. In the below illustration the horizontal ray cast from black point intersects the polygon 3 times. Thus the black point is inside the polygon.

There are a few corner cases that have to be carefully handled.

Case 1. The point p is on the boundary of the polygon. In this case the algorithm may return true or false because of floating point issues.

Case 2. The casted ray intersects a vertex of the polygon. In this case the naive implementation of the algorithm would return two intersections because of the top line segment and the bottom line segment. The way to fix this is to only count an intersection if (1) it’s not an end-point intersection or (2) it’s an end-point intersection and the end-point of the line segment lies below the casted ray.

Case 3. The casted ray overlaps a line segment, l. In this case the casted ray intersects an end-point of the overlapping line segment, l and the line segment does not lie below the casted ray. So do not count the line segment l intersection(s).

By implementing the special cases this point-in-polygon method should work extremely well. It’s a fast method since the only operations used are compares. The big-O is O(n) where n is the number of line segments for the polygon.

More on Floating Point

In my last post I outlined one of the issues with floating point. My particular geometric problem is testing if two polygons are subsets of each other. In general I want to test if one polygon is a subset of another polygon. Also the polygons can be rotated, translated, or reflected in particular any rigid motion can be applied to the polygons.

To again illustrate the problem with floating point I’ll give an example:

Let polygon, T1 = \{(0, 0), (1, 0), (1, 1),\} (a triangle), and polygon T2 = \{ (0, 0), (0.1, 0.1), (0, 1)\} (another triangle). See the below illustration.
Two triangles
By rotating triangle T1 by exactly -\pi/4 radians, I can make triangle T1 a subset of triangle T2. But the problem with floating point is that -\pi/4 is not a floating point number so I can’t rotate by exactly -\pi/4. Which means that using the ordinary algorithms for subsetness there is no easy way to make T1 a subset of T2.

That why I wrote the previous post on using some sort of fuzzy \epsilon in my geometry algorithms. But I believe that there is another way to get the result I want. What I can do is rotate T1 by the closest floating point number to -\pi/4 and then do some boolean operations to T1 and T2. In particular take the boolean intersection of T1 and T2, call the resulting shape I (the intersection). And then take the boolean difference of T1 and I, call the resulting shape R (the remainder). If R is small enough (say the area is within \epsilon = 10^-5) then I can say that T1 is a subset of T2.

In this way I can avoid using \epsilon in the algorithm for testing if a point is contained in a polygon (which is what I was supposed to write about in this post).

Rigorous Geometry

I’ve mentioned before that floating numbers are not real numbers. In particular they are not associative or distributive link
and also most real numbers cannot be exactly represented. Thus a naive implementation for common geometry operations such as line intersection or testing if a point is contained (by contained I mean the point is either on the boundary or in the interior of the polygon) inside a polygon can have problems. I’ll give a short example of testing if a point is inside a square.

Let p = (0.5, 0) and let the square, s, have vertices at (0, 0), (1, 0), (1, 1), (0, 1) . Obviously p is contained in s and any reasonable algorithm for testing so would yield true. I now want to try rotating s by \pi/4 radians. This is the point where the pesky problems with floating point come in. The number \pi/4 is not a floating point number so it is approximated by 0.785398164 (I haven’t written all the digits of the floating point number but that doesn’t change the issue). It turns out that 0.785398164 is greater than \pi/4. That is the square was rotated by a little bit more than \pi/4 radians. Thus the point (0.5, 0.5) is not contained in the rotated square.

This is a bad thing because intuitively the point (0.5, 0.5) should be contained in the rotated square. After all we meant to rotate the square by \pi/4 radians. I don’t know of any way to fix the problem where we can’t rotate by the number of desired radians. What is possible is to modify the algorithms for testing if a point is contained inside a polygon so that they return true if the point is within \epsilon of the boundary of the polygon or is actually contained in the polygon. The \epsilon is a fixed number that is large enough to mask the floating point rounding problems.

My next post will give one of the algorithms for testing point containedness with the \epsilon modification.