How can I determine whether a 2D Point is within a Polygon?

I’m trying to create a fast 2D point inside polygon algorithm, for use in hit-testing (e.g. Polygon.contains(p:Point) ). Suggestions for effective techniques would be appreciated.




Answer

For graphics, I’d never go with Integers. Many OSes use Integers for painting UI (pixels are ints after all), but Mac OS X for example uses float for everything (it talks of points. A point can translate to one pixel, but depending on monitor resolution, it might translate to anything else but a pixel and if resolution is very high of a monitor, in theory a pixel can translate to half a point, so 1.51.5 can be a different pixel than 11 and 22). However, I never noticed that Mac UIs are significantly slower than other UIs. After all 3D (OpenGL or Direct3D) also works with floats all of the time and modern graphics libraries might very often take advantage of high power GPUs of a system without you even noticing it (as a GPU can also speed up 2D painting, not just 3D; 2D is only a subset of 3D, consider 2D to be 3D where the Z coordinates are always 0 for example).

Now you said speed is your main concern, okay, let’s go for speed. Before you run any sophisticated algorithm, first do a simple test. Create an axis aligned bounding box around your polygon. This is very easy, fast and can already safe you tons of CPU time. How does that work? Iterate over all points of the polygon (every point being X/Y) and find the min/max values of X and Y. E.g. you have the points (91), (43), (27), (82), (36). Xmin is 2, Xmax is 9, Ymin is 1 and Ymax is 7. Now you know that no point within your polygon can ever have a x value smaller than 2 and greater than 9, no point can have a y value smaller than 1 and greater than 7. This way you can quickly exclude many points not being within your polygon:

// p is your point, p.x is the x coord, p.y is the y coord
if (p.x < Xmin || p.x > Xmax || p.y < Ymin || p.y > Ymax) {
    // Definitely not within the polygon!
}

This is the first test to run on any point. If this test already excludes the point from the polygon (even though it’s a very coarse test!), then there is no use of running any other tests. As you can see, this test is ultra fast.

If this test won’t exclude the point, it is within the axis aligned bounding box, but that does not mean it is within the polygon; it only means it might be within the polygon. What we need next is a more sophisticated test to check if the point is really within the polygon or just within the bounding box. There are a couple of ways how this can be calculated. It also makes a huge difference is the polygon can have holes or whether its solid. Here are examples of solid ones (one convex, one concave):

Polygon without hole

And here’s one with a hole:

Polygon with hole

The green one has a hole in the middle!

The easiest way is to use ray casting , since it can handle all the polygons shown above correctly, no special handling is necessary and it still provides good speed. The idea of the algorithm is pretty simple: Draw a virtual ray from anywhere outside the polygon to your point and count how often it hits any side of the polygon. If the number of hits is even, it’s outside of the polygon, if it’s odd, it’s inside.

Demonstrating how the ray cuts through a polygon

The winding number algorithm is more accurate for points being very, very, very close to a polygon line; ray casting may fail here because of float precision and rounding issues, however winding number is much slower and if a point is accidentally detected to be outside of the polygon if it’s so close to a polygon line, that your eye can’t even tell if it’s inside or outside, is it really a problem? I don’t think so, so let’s keep things simple.

You still have the bounding box of above, remember? Okay, now let’s say your point is p again (p.x/p.y). Your ray might go from

  • (Xmin - e/p.y) to (p.x/p.y) or
  • (p.x/p.y) to (Xmax + e/p.y) or
  • (p.x/Ymin - e) to (p.x/p.y) or
  • (p.x/p.y) to (p.x/Ymax + e)

It won’t matter. Choose whatever sounds best to you. It’s only important that the ray starts definitely outside of the polygon and stops at the point.

So what is e anyway? Well, e (actually epsilon) gives the bounding box some padding . As I said, ray tracing fails if we start too close to a polygon line. Since the bounding box might equal the polygon (if the polygon is actually an axis aligned rectangle, the bounding box is equal to the polygon itself! And the ray would start directly on the polygon side). How big should you choose e? Not too big. It depends on the coordinate system scale you use for drawing. You could select e to be 1.0, however, if your polygons have coordinates much smaller than 1.0, selecting e to be 0.001 might be large enough. You could select e to be always 1% of the polygon size, e.g. when having a ray along the x axis, you could calculate e like this:

e = ((Xmax - Xmin) / 100)

Now that we have the ray with its start and end coordinates, the problem shifts from “is the point within the polygon” to “how often does the ray intersect the polygon side”. Therefore we can’t just work with the polygon points as before (for the bounding box), now we need the actual sides. A side is always defined by two points.

side 1: (X1/Y1)-(X2/Y2)
side 2: (X2/Y2)-(X3/Y3)
side 3: (X3/Y3)-(X4/Y4)
:

You need to test the ray against all sides. Consider the ray to be a vector and every side to be a vector. The ray has to hit each side exactly once or never at all. It can’t hit the same side twice (two lines in 2D space will always intersect exactly once, unless they are parallel, in which case they never intersect. However since vectors have a limited length, two vectors might not be parallel and still never intersect).

// Test the ray against all sides
int intersections = 0;
for (side = 0; side < numberOfSides; side++) {
    // Test if current side intersects with ray.
    // If yes, intersections++;
}
if ((intersections & 1) == 1) {
    // Inside of polygon
} else {
    // Outside of polygon
}

So far so well, but how do you test if two vectors intersect? Here’s some C code (not tested), that should do the trick:

#define NO 0
#define YES 1
#define COLLINEAR 2

int areIntersecting(
    float v1x1, float v1y1, float v1x2, float v1y2,
    float v2x1, float v2y1, float v2x2, float v2y2
) {
    float d1, d2;
    float a1, a2, b1, b2, c1, c2;

    // Convert vector 1 to a line (line 1) of infinite length.
    // We want the line in linear equation standard form: A*x + B*y + C = 0
    // See: http://en.wikipedia.org/wiki/Linear_equation
    a1 = v1y2 - v1y1;
    b1 = v1x1 - v1x2;
    c1 = (v1x2 * v1y1) - (v1x1 * v1y2);

    // Every point (x,y), that solves the equation above, is on the line,
    // every point that does not solve it, is not. The equation will have a
    // positive result if it is on one side of the line and a negative one 
    // if is on the other side of it. We insert (x1,y1) and (x2,y2) of vector
    // 2 into the equation above.
    d1 = (a1 * v2x1) + (b1 * v2y1) + c1;
    d2 = (a1 * v2x2) + (b1 * v2y2) + c1;

    // If d1 and d2 both have the same sign, they are both on the same side
    // of our line 1 and in that case no intersection is possible. Careful, 
    // 0 is a special case, that's why we don't test ">=" and "<=", 
    // but "<" and ">".
    if (d1 > 0 && d2 > 0) return NO;
    if (d1 < 0 && d2 < 0) return NO;

    // The fact that vector 2 intersected the infinite line 1 above doesn't 
    // mean it also intersects the vector 1. Vector 1 is only a subset of that
    // infinite line 1, so it may have intersected that line before the vector
    // started or after it ended. To know for sure, we have to repeat the
    // the same test the other way round. We start by calculating the 
    // infinite line 2 in linear equation standard form.
    a2 = v2y2 - v2y1;
    b2 = v2x1 - v2x2;
    c2 = (v2x2 * v2y1) - (v2x1 * v2y2);

    // Calculate d1 and d2 again, this time using points of vector 1.
    d1 = (a2 * v1x1) + (b2 * v1y1) + c2;
    d2 = (a2 * v1x2) + (b2 * v1y2) + c2;

    // Again, if both have the same sign (and neither one is 0),
    // no intersection is possible.
    if (d1 > 0 && d2 > 0) return NO;
    if (d1 < 0 && d2 < 0) return NO;

    // If we get here, only two possibilities are left. Either the two
    // vectors intersect in exactly one point or they are collinear, which
    // means they intersect in any number of points from zero to infinite.
    if ((a1 * b2) - (a2 * b1) == 0.0f) return COLLINEAR;

    // If they are not collinear, they must intersect in exactly one point.
    return YES;
}

The input values are the two endpoints of vector 1 (x and y) and vector 2 (x and y). So you have 2 vectors, 4 points, 8 coordinates. YES and NO are clear. YES increases intersections, NO does nothing. What about COLLINEAR? It means both vectors lie on the same infinite line, depending on position and length, they don’t intersect at all or they intersect in an endless number of points. I’m not absolutely sure how to handle this case, I would not count it as intersection either way. Well, this case is rather rare in practice anyway because of floating point rounding errors; better code would probably not test for == 0.0f but instead for something like < epsilon , where epsilon is a rather small number (so smaller rounding mistakes won’t lead to wrong results).

If you need to test a larger number of points, you can maybe speed up the whole thing a bit by keeping the linear equation standard forms of the polygon sides in memory, so you don’t have to recalculate these every time. This will save you two floating point multiplications and three floating point subtractions on every test in exchange for storing three floating point values per polygon side in memory; it’s the typical calculations vs memory trade off. On a very fast CPU with very slow memory and too little cache, this may in fact not be faster at all, so be sure to benchmark if that really pays off.

Last but not least: If you may use 3D hardware to solve the problem, forget about anything above. There is a much faster and much easier way. Just let the GPU do all the work for you. Create a painting surface that is off screen (so you can paint into it, without it appearing anywhere on the screen). Fill it completely with the color black. Now let OpenGL or Direct3D paint your polygon (or even all of your polygons if you just want to test if the point is within any of them, but you don’t care for which one) into this drawing surface and fill the polygon(s) with a different color, e.g. white. To check if a point is within the polygon, get the color of this point from the drawing surface. This is just a O(1) memory fetch. If it’s white, it’s inside, if it’s black, it’s outside. Easy, isn’t it? This method will pay off if you have very little polygons (e.g. 50-100), but a damn lot of points to test (> 1000), in which case this method is much speedier than anything else. It will only be a problem if your drawing surface must be huge because your polygons are. If your drawing surface needs to be 100 MB or more to make the polygons fit, this method might become very slow (despite the fact that it wastes tons of memory).