Point in Polygon aka hit test

Question

What's your best, most elegant 2D "point inside polygon" or Polygon.contains(p:Point) algorithm?

Edit: There may be different answers for floats vs integers. Primary goal is speed.

Answer 1

The trivial solution would be to divide the polygon to triangles and hit test the triangles as explained here

If your polygon is CONVEX there might be a better approach though. Look at the polygon as a collection of infinite lines. Each line dividing space into two. for every point it's easy to say if its on the one side or the other side of the line. If a point is on the same side of all lines then it is inside the polygon.

Answer 2

I did some work on this back when I was a researcher under Michael Stonebraker - you know, the Professor who came up with Ingres, Postgres, etc.

We realized that the fastest way was to first do a bounding box because it's SUPER fast. If it's outside the bounding box, it's outside. Otherwise, you do the harder work...

If you want a great algorithm, look to the open source project Postgres source code for the geo work...

I want to point out, we never got any insight into right vs left handedness (also expressible as an "inside" vs "outside" problem...

---

UPDATE

BKB's link provided a good number of reasonable algorithms. I was working on Earth Science problems and therefore needed a solution that works in latitude/longitude, and it has the peculiar problem of handedness - is the area inside the smaller area or the bigger area? The answer is that the "direction" of the verticies matters - it's either left-handed or right handed and in this way you can indicate either area as "inside" any given polygon. As such, my work used solution three enumerated on that page.

In addition, my work used separate functions for "on the line" tests.

...Since someone asked: we figured out that bounding box tests were best when the number of verticies went beyond some number - do a very quick test before doing the longer test if necessary... A bounding box is created by simply taking the largest x, smallest x, largest y and smallest y and putting them together to make four points of a box...

Another tip for those that follow: we did all our more sophisticated and "light-dimming" computing in a grid space all in positive points on a plane and then re-projected back into "real" longitude/latitude, thus avoiding possible errors of wrapping around when one crossed line 180 of longitude and when handling polar regions. Worked great!

Answer 3

Compute the oriented sum of angles between the point p and each of the polygon apices. If the total oriented angle is 360 degrees, the point is inside. If the total is 0, the point is outside.

I like this method better because it is more robust and less dependent on numerical precision.

Methods that compute evenness of number of intersections are limited because you can 'hit' an apex during the computation of the number of intersections.

EDIT: By The Way, this method works with concave and convex polygons.

Answer 4

David Segond's answer is pretty much the standard general answer, and Richard T's is the most common optimization, though therre are some others. Other strong optimizations are based on less general solutions. For example if you are going to check the same polygon with lots of points, triangulating the polygon can speed things up hugely as there are a number of very fast TIN searching algorithms. Another is if the polygon and points are on a limited plane at low resolution, say a screen display, you can paint the polygon onto a memory mapped display buffer in a given colour, and check the color of a given pixel to see if it lies in the polygons.

Like many optimizations, these are based on specific rather than general cases, and yield beneifits based on amortized time rather than single usage.

Working in this field, i found Joeseph O'Rourkes 'Computation Geometry in C' ISBN 0-521-44034-3 to be a great help.

Answer 5

I too thought 360 summing only worked for convex polygons but this isn't true.

This site has a nice diagram showing exactly this, and a good explanation on hit testing: Gamasutra - Crashing into the New Year: Collision Detection

Answer 6

Eric Haines's survey of approaches

Answer 7

The Eric Haines article cited by bobobobo is really excellent. Particularly interesting are the tables comparing performance of the algorithms; the angle summation method is really bad compared to the others. Also interesting is that optimisations like using a lookup grid to further subdivide the polygon into "in" and "out" sectors can make the test incredibly fast even on polygons with > 1000 sides.

Anyway, it's early days but my vote goes to the "crossings" method, which is pretty much what Mecki describes I think. However I found it most succintly described and codified by David Bourke. I love that there is no real trigonometry required, and it works for convex and concave, and it performs reasonably well as the number of sides increases.

By the way, here's one of the performance tables from the Eric Haines' article for interest, testing on random polygons.

                       number of edges per polygon
                         3       4      10      100    1000
MacMartin               2.9     3.2     5.9     50.6    485
Crossings               3.1     3.4     6.8     60.0    624
Triangle Fan+edge sort  1.1     1.8     6.5     77.6    787
Triangle Fan            1.2     2.1     7.3     85.4    865
Barycentric             2.1     3.8    13.8    160.7   1665
Angle Summation        56.2    70.4   153.6   1403.8  14693

Grid (100x100)          1.5     1.5     1.6      2.1      9.8
Grid (20x20)            1.7     1.7     1.9      5.7     42.2
Bins (100)              1.8     1.9     2.7     15.1    117
Bins (20)               2.1     2.2     3.7     26.3    278

Answer 8

A 90 sided polygon is called a enneacontagon and a 10,000 sided polygon is called a myriagon.

Answer 9

I think the following piece of code is the best solution (taken from here):

int pnpoly(int nvert, float *vertx, float *verty, float testx, float testy)
{
  int i, j, c = 0;
  for (i = 0, j = nvert-1; i < nvert; j = i++) {
    if ( ((verty[i]>testy) != (verty[j]>testy)) &&
     (testx < (vertx[j]-vertx[i]) * (testy-verty[i]) / (verty[j]-verty[i]) + vertx[i]) )
       c = !c;
  }
  return c;
}

It's both short and efficient and works both for convex and concave polygons. As suggested before, you should check the bounding rectangle first and treat polygon holes separately.

The idea behind this is pretty simple. It is based on the observation that a test point is within a polygon if when projected on the y-axis it's x value is below odd number of polygon edges.