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.

Translate

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.

所有的回答

For graphics, I'd rather not prefer integers. Many systems use integers for UI painting (pixels are ints after all), but macOS for example uses float for everything. macOS only knows points and a point can translate to one pixel, but depending on monitor resolution, it might translate to something else. On retina screens half a point (0.5/0.5) is pixel. Still, I never noticed that macOS UIs are significantly slower than other UIs. After all 3D APIs (OpenGL or Direct3D) also works with floats and modern graphics libraries very often take advantage of GPU acceleration.

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 boxaround your polygon. This is very easy, fast and can already safe you a lot of calculations. How does that work? Iterate over all points of the polygon and find the min/max values of X and Y.E.g. you have the points

`(9/1), (4/3), (2/7), (8/2), (3/6)`

. This means Xmin is 2, Xmax is 9, Ymin is 1 and Ymax is 7. A point outside of the rectangle with the two edges (2/1) and (9/7) cannot be within the polygon.This is the first test to run for any point. As you can see, this test is ultra fast but it's also very coarse. To handle points that are within the bounding rectangle, we need a more sophisticated algorithm. There are a couple of ways how this can be calculated. Which method works also depends on the fact if the polygon can have holes or will always be solid. Here are examples of solid ones (one convex, one concave):

And here's one with a hole:

The green one has a hole in the middle!

The easiest algorithm, that can handle all three cases above and is still pretty fast is named

ray casting. 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 a side of the polygon. If the number of hits is even, it's outside of the polygon, if it's odd, it's inside.The

winding number algorithmwould be an alternative, it is more accurate for points being very close to a polygon line but it's also much slower. Ray casting may fail for points too close to a polygon side because of limited floating point precision and rounding issues, but in reality that is hardly a problem, as if a point lies that close to a side, it's often visually not even possible for a viewer to recognize if it is already inside or still outside.You still have the bounding box of above, remember? Just pick a point outside the bounding box and use it as starting point for your ray. E.g. the point

`(Xmin - e/p.y)`

is outside the polygon for sure.But what is

`e`

? Well,`e`

(actually epsilon) gives the bounding box somepadding. 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 an axis aligned rectangle, the bounding box is equal to the polygon itself!), we need some padding to make this safe, that's all. How big should you choose`e`

? Not too big. It depends on the coordinate system scale you use for drawing. If your pixel step width is 1.0, then just choose 1.0 (yet 0.1 would have worked as well)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 intersects a polygon side". Therefore we can't just work with the polygon points as before, now we need the actual sides. A side is always defined by two points.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 because they are too short to ever meet each other.

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:

The input values are the

two endpointsof vector 1 (`v1x1/v1y1`

and`v1x2/v1y2`

) and vector 2 (`v2x1/v2y1`

and`v2x2/v2y2`

). 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.If you need to test a larger number of points, you can certainly 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 a typical memory vs computation time trade off.

Last but not least: If you may use 3D hardware to solve the problem, there is an interesting alternative. Just let the GPU do all the work for you. Create a painting surface that is off 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) 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.

Of course this method is only usable if your drawing surface doesn't have to be huge. If it cannot fit into the GPU memory, this method is slower than doing it on the CPU. If it would have to be huge and your GPU supports modern shaders, you can still use the GPU by implementing the ray casting shown above as a GPU shader, which absolutely is possible. For a larger number of polygons or a large number of points to test, this will pay off, consider some GPUs will be able to test 64 to 256 points in parallel. Note however that transferring data from CPU to GPU and back is always expensive, so for just testing a couple of points against a couple of simple polygons, where either the points or the polygons are dynamic and will change frequently, a GPU approach will rarely pay off.

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

## Arguments

nvert: Number of vertices in the polygon. Whether to repeat the first vertex at the end has been discussed in the article referred above.vertx, verty: Arrays containing the x- and y-coordinates of the polygon's vertices.testx, testy: X- and y-coordinate of the test point.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. The author describes it as follows:

The variable c is switching from 0 to 1 and 1 to 0 each time the horizontal ray crosses any edge. So basically it's keeping track of whether the number of edges crossed are even or odd. 0 means even and 1 means odd.

Here is a C# version of the answer given by nirg, which comes from this RPI professor. Note that use of the code from that RPI source requires attribution.

A bounding box check has been added at the top. However, as James Brown points out, the main code is almost as fast as the bounding box check itself, so the bounding box check can actually slow the overall operation, in the case that most of the points you are checking are inside the bounding box. So you could leave the bounding box check out, or an alternative would be to precompute the bounding boxes of your polygons if they don't change shape too often.

Here is a JavaScript variant of the answer by M. Katz based on Nirg's approach:

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.

EDIT: I recently found a whole Wikipedia article on the topic.

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.

This question is so interesting. I have another workable idea different from other answers of this post.

The idea is to use the sum of angles to decide whether the target is inside or outside. If the target is inside an area, the sum of angle form by the target and every two border points will be 360. If the target is outside, the sum will not be 360. The angle has direction. If the angle going backward, the angle is a negative one. This is just like calculating the winding number.

Refer to this image to get a basic understanding of the idea:

My algorithm assumes the clockwise is the positive direction. Here is a potential input:

The following is the python code that implements the idea:

Swift version of the answer by nirg:

I did some work on this back when I was a researcher under Michael Stonebraker - you know, the professor who came up with Ingres, PostgreSQL, 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 PostgreSQL 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!

Really like the solution posted by Nirg and edited by bobobobo. I just made it javascript friendly and a little more legible for my use:

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.

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

If your polygon is

CONVEXthere 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.I realize this is old, but here is a ray casting algorithm implemented in Cocoa, in case anyone is interested. Not sure it is the most efficient way to do things, but it may help someone out.

Obj-C version of nirg's answer with sample method for testing points. Nirg's answer worked well for me.

C# version of nirg's answer is here: I'll just share the code. It may save someone some time.

There is nothing more beutiful than an inductive definition of a problem. For the sake of completeness here you have a version in prolog which might also clarify the thoughs behind

ray casting:Based on the simulation of simplicity algorithm in http://www.ecse.rpi.edu/Homepages/wrf/Research/Short_Notes/pnpoly.html

Some helper predicates:

The equation of a line given 2 points A and B (Line(A,B)) is:

It is important that the direction of rotation for the line is setted to clock-wise for boundaries and anti-clock-wise for holes. We are going to check whether the point (X,Y), i.e the tested point is at the left half-plane of our line (it is a matter of taste, it could also be the right side, but also the direction of boundaries lines has to be changed in that case), this is to project the ray from the point to the right (or left) and acknowledge the intersection with the line. We have chosen to project the ray in the horizontal direction (again it is a matter of taste, it could also be done in vertical with similar restrictions), so we have:

Now we need to know if the point is at the left (or right) side of the line segment only, not the entire plane, so we need to restrict the search only to this segment, but this is easy since to be inside the segment only one point in the line can be higher than Y in the vertical axis. As this is a stronger restriction it needs to be the first to check, so we take first only those lines meeting this requirement and then check its possition. By the Jordan Curve theorem any ray projected to a polygon must intersect at an even number of lines. So we are done, we will throw the ray to the right and then everytime it intersects a line, toggle its state. However in our implementation we are goint to check the lenght of the bag of solutions meeting the given restrictions and decide the innership upon it. for each line in the polygon this have to be done.

Java Version:

.Net port:

VBA VERSION:Note: Remember that if your polygon is an area within a map that Latitude/Longitude are Y/X values as opposed to X/Y (Latitude = Y, Longitude = X) due to from what I understand are historical implications from way back when Longitude was not a measurement.

CLASS MODULE: CPoint

MODULE:

I've made a Python implementation of nirg's c++ code:

Inputs

bounding_points:nodes that make up the polygon.bounding_box_positions:candidate points to filter. (In my implementation created from the bounding box.(The inputs are lists of tuples in the format:

`[(xcord, ycord), ...]`

)Returns

Again, the idea is taken from here

Heres a point in polygon test in C that isn't using ray-casting. And it can work for overlapping areas (self intersections), see the

`use_holes`

argument.Note: this is one of the less optimal methods since it includes a lot of calls to

`atan2f`

, but it may be of interest to developers reading this thread (in my tests its ~23x slower then using the line intersection method).For Detecting hit on Polygon we need to test two things:

To deal with the following special cases in Ray casting algorithm:

Check Determining Whether A Point Is Inside A Complex Polygon. The article provides an easy way to resolve them so there will be no special treatment required for the above cases.

You can do this by checking if the area formed by connecting the desired point to the vertices of your polygon matches the area of the polygon itself.

Or you could check if the sum of the inner angles from your point to each pair of two consecutive polygon vertices to your check point sums to 360, but I have the feeling that the first option is quicker because it doesn't involve divisions nor calculations of inverse of trigonometric functions.

I don't know what happens if your polygon has a hole inside it but it seems to me that the main idea can be adapted to this situation

You can as well post the question in a math community. I bet they have one million ways of doing that

If you are looking for a java-script library there's a javascript google maps v3 extension for the Polygon class to detect whether or not a point resides within it.

Google Extention Github

When using qt (Qt 4.3+), one can use QPolygon's function containsPoint

The answer depends on if you have the simple or complex polygons. Simple polygons must not have any line segment intersections. So they can have the holes but lines can't cross each other. Complex regions can have the line intersections - so they can have the overlapping regions, or regions that touch each other just by a single point.

For simple polygons the best algorithm is Ray casting (Crossing number) algorithm. For complex polygons, this algorithm doesn't detect points that are inside the overlapping regions. So for complex polygons you have to use Winding number algorithm.

Here is an excellent article with C implementation of both algorithms. I tried them and they work well.

http://geomalgorithms.com/a03-_inclusion.html

Scala version of solution by nirg (assumes bounding rectangle pre-check is done separately):

Here is golang version of @nirg answer (inspired by C# code by @@m-katz)

Surprised nobody brought this up earlier, but for the pragmatists requiring a database: MongoDB has excellent support for Geo queries including this one.

What you are looking for is:

`Neighborhoods`

is the collection that stores one or more polygons in standard GeoJson format. If the query returns null it is not intersected otherwise it is.Very well documented here: https://docs.mongodb.com/manual/tutorial/geospatial-tutorial/

The performance for more than 6,000 points classified in a 330 irregular polygon grid was less than one minute with no optimization at all and including the time to update documents with their respective polygon.