2D point in polygon algorithm

Game Mechanics
#1

I was searching for a good point-in-polygon algorithm. The most popular one for arbitrary concave and convex polygons without holes seems to be the “ray casting algorithm.”
It is excellently described on http://alienryderflex.com/polygon/.
Now this algorithm has a problem that it always linearly tests all polygon edges against the query point (which is really unnecessary) and when using polygons with many vertices (like even more than 256) it goes down on its knees.
I was searching for a better way to do point in polygon queries but couldn’t really find any. So I took that ray casting algorithm and simply extended it with an interval tree to dramatically reduce the number of polygon edges that need to be tested. Actually, only the edges that intersect the imaginary ray with the y co-ordinate of the query point are tested now.
The general idea was that the two y co-ordinates of a polygon edge would form an interval that the interval tree allows to query for very efficiently for a given point.
With that solution, a convex polygon with around 256,000 vertices can be tested against around 4,000,000 query points in about a second.
The runtime only really depends on the number of polygon edges that intersect a given horizontal line with a given y co-ordinate.

Now I’d like to ask for improvements to have a final version ready for JOML’s geometry framework.

The current source is here: https://raw.githubusercontent.com/JOML-CI/JOML/geometry/src/org/joml/PolygonPointIntersection.java

Especially bad is the allocations of the Interval and IntervalTree objects. Maybe someone has a nifty idea how to do that better and with good memory locality.

#2

This is pretty much the way I have done it (except for the interval tree, I only had small polygon counts to deal with), what I would suggest though is to use the winding number to determine in or out as the odd / even method breaks down in certain circumstances. Is your interval tree like a quad tree?

#3

Your optimisation is clever.
Just one small thing to keep in mind. Since your ray is a horizontal line, and many games use shapes with horizontal lines such as axis aligned bounding boxes aabb’s, the ray-polygon line test may often be a collinear overlapping line test which is a special case that will cause problems with your line line intersection method.

One method to overcome the problem but still keep your clever optimisation is to offset the testing point slightly by adding a small amount like 0.0001 to its y coordinate so then the casting ray is less likely to be collinear with the polygon line.

EDIT: an interesting trove of 2D geometry is JTS, java topology suite. It prioritises stability over speed, however the algorithms are interesting.

#4

Thank you @ziozio and @CommanderKeith for your suggestions and the reference to JTS.
Yes, the even/odd algorithm really does have a weakness with collinear edges. Either one can implement it that a point is inside the polygon when it lies on the bottom/left/minimum side of a polygon’s edge or on the top/right/maximum edge.
CommanderKeith, I’ll add your suggestion with the small delta/epsilon.

#5

I generalized the algorithm a bit to handle a set of polygons instead of just one polygon: https://raw.githubusercontent.com/JOML-CI/JOML/geometry-2d/src/org/joml/PolygonsPointIntersection.java
The algorithm handles all polygons with the same interval tree and the query method returns a BitSet where the i-th bit is set iff the query point lies inside the i-th polygon.

To handle non-overlapping polygons a simple bounding volume or binary space partitioning acceleration structure would probably be better to first cull away polygons by their bounding boxes, rather than to insert all edges of all polygons in the same interval tree.
The goal is to have an efficient algorithm for collision detection with a level containing many (hand-drawn) arbitrarily shaped polygons. A level designer would be able to simply “draw” such levels by drawing the polygons.
With every mouse move in the paint program (i.e. level editor) that editor would add a new vertex to the currently drawn polygon.
Therefore I wanted to be able to handle thousands of vertices per polygon to have really huge and detailed levels.

Next optimization I am going to do is to use different ray directions in the raycasting query. Currently with the original algorithm the ray is cast along the same y co-ordinate. This is suboptimal if the polygon(s) contain many (nearly) vertical edges with approximately the same y position, such as when drawing the side-view of a bumpy terrain with hills that go up and down in waves.
In this scenario it would be best to actually cast a ray along the same x co-ordinate.
But that can even be generalized to using an arbitrary ray direction and projecting the polygon edges onto that direction.
Always with the goal of minimizing the number of polygon edges that the interval tree may output and that then must be tested with the even/odd ray casting algorithm.

Let’s see what makes sense and what works best.

#6

Don’t forget to add usual bounding box fast path early out test.

#7

The joml-lwjgl3-demos repository now has a simple PolygonDrawer demo.
You can hold the mouse button down and drag the mouse around to draw any kind of polygon (convex, concave, self-intersecting, counter-clockwise or clockwise). Afterwards move the mouse in and out of the polygon to color the outline red when the mouse cursor is inside and black when outside.

#8

In case you cannot/wantnot :slight_smile: build, here is a jar of that: https://www.dropbox.com/s/lqt8xdkmg9euvq5/polydraw.jar?dl=0
Uses LWJGL 3 with OpenGL 1.1 and GLFW and works under Windows, Linux, OS X.

#9

the jar is missing the linux and os x lwjgl3 natives :slight_smile:

#10

Dang it… I should not claim something which I did not test myself. :smiley:

For Linux and OS X people, here is a video :smiley: -> http://i.imgur.com/hPTRIW7.gifv

EDIT: Updated with Linux and OS X libs.

#11

Thats very cool stuff, a polygon of that complexity is not something i’d have ever considered testing using its edges. An approach such as using bit masks (like Worms) to test such complex shapes might have been easier to implement (also would handle holes). Nonetheless very impressive results.

btw on OS X i get the following error message when running the jar

java -jar polydraw.jar
Exception in thread "main" java.lang.UnsatisfiedLinkError: org.lwjgl.system.MemoryAccess.getPointerSize()I
	at org.lwjgl.system.MemoryAccess.getPointerSize(Native Method)
	at org.lwjgl.system.Pointer.<clinit>(Pointer.java:24)
	at org.lwjgl.glfw.GLFW.<clinit>(GLFW.java:594)
	at org.joml.lwjgl.PolygonDrawer.run(PolygonDrawer.java:51)
	at org.joml.lwjgl.PolygonDrawer.main(PolygonDrawer.java:176)
#12

Thanks! :slight_smile:
It just keeps amazing me that apparently no one else (also on the various polygon/point intersection algorithms on GitHub) thought of actually accelerating the query for the edges that could possibly intersect the ray in that raycast algorithm…
It is simple after all. I implemented the interval tree after reading once about it on the Wikipedia site. Nothing fancy in there, actually. :slight_smile:

Handling holes is also easy with this approach:
Just check whether the point is inside the polygon of a hole. This can also be very efficiently implemented by first testing whether the point intersects the bounding sphere/rectangle of the hole, like with doing a normal polygon/point intersection.
If it is in the hole and also in the outer polygon, then -> no intersection. :slight_smile:
And because the algorithm is damn fast, we can throw hundreds and even thousands of holes and polygons at it without it even batting an eyelash. Having drawn a 6781 vertices polygon and measuring with System.nanoTime() the query time, it was 0.004 milliseconds…, so just 4 microseconds.

But that OS X issue is strange. I took the latest oss.sonatype.org deployments with Maven.

EDIT: Improved the performance of the tree traversal again by 10%: Before descending into the left child of the interval tree, check if it actually contains an interval whose maximum value is at least the queried ‘y’ (that maximum value is cached during tree construction).
Likewise for the right child, check if it contains an interval whose minimum value is at most ‘y’.
A 1,048,576 vertices polygon can now be queried around 4,664,358 times per second when all query points lie within the polygon.

#13

Kai

If you draw a 5 pointed star (https://en.wikipedia.org/wiki/List_of_regular_polytopes_and_compounds#/media/File:Star_polygon_5-2.svg) you get overlapping (but its a continuous drawn polygon), the centre part of the star is currently seen as “out” when it is actually “in” but the 5 points of the star show correctly “in”.

This is an example where the odd / even test fails.

The winding number algo is more accurate and potentially faster (http://geomalgorithms.com/a03-_inclusion.html). Most packages seem to offer both methods so it might be good to be able to select either algorithm depending on whether you have basic geometry or complex.

#14

That’s true. Such polygons don’t work with the even/odd/raycast/crosscutting algorithm. It’s questionable whether actually all faces of that polygon should indeed be seen as “inside” of the polygon.
Other authors argue quite differently: http://www.cs.berkeley.edu/~ddgarcia/cs184/r3/#02.1
It can be argued that modeling holes is possible using this.
It would indeed be good to have the option to choose from both.
But it will be quite some work to get the winding order algorithm at even the same order of magnitude of performance of what I currently have with the modified even/odd algorithm. :slight_smile: Some k-nearest neighbor search with an appropriate data structure seems to be called for with that algorithm.

#15

That’s interesting. But the exclusion of self-intersecting polygons is a reasonable restriction.

#16

Yepp, handling self-intersecting polygons would really be nice to have!
I think we can modify the crosscutting/raycast/even-odd algorithm to support self-intersecting polygons.
All that needs is to first find the two points of intersection that gets a side of the polygon from being treated as inside to being treated as outside and change the order of the vertices (reverse the edge direction) of that edge strip between the two intersections.
In order to detect the intersections, we can reuse the interval tree during its construction to query it for all edges with at least the same y interval as the edge we are trying to add. Or another structure that can be queried for a 2D interval/area, such as a quadtree, might be better suited for this.
This all must be done in the order the polygon edges are drawn, so it makes the currently easy interval tree construction a bit harder, because the tree now needs to be dynamic and balance itself accordingly as we add new edges to it, because we cannot simply sort the edges anymore and simply distribute it to the tree.
But I am somewhat determined to stay with the raycast algorithm, because querying is so freaking fast now. :slight_smile:

#17

I am proud to present that the algorithm now supports holes and multipolygons. :smiley:
The API is as follows:


PolygonPointIntersection(float[] verticesXY, int[] polygons, int count)

Invoke this constructor with the list of all vertices. These vertices can make up many polygons, including holes. The second parameter “polygons” is an int[] array which defines at which index into the verticesXY a new polygon starts.
The first polygon (which is always an “inner” polygon) starts always at index 0.
If there is only one polygon, the “polygons” argument is simply an empty int[] array.
The “count” parameter is simply how many of the vertices inside verticesXY should be used in total, since you may want to allocate one large float[] array first and then populate it with some vertices.
This actually works great!

EDIT:
Also the query method got simplified. Now it only looks like this:


boolean pointInPolygon(float x, float y)

is fully multithreaded/thread-safe and does not have that annoying “working” array parameter anymore, because now the actual even/odd/crosscutting/raycast algorithm has been merged into the interval tree traversal, so it does also not produce a list of found edges anymore but immediately performs the even/odd algorithm on the edges as they are found.
That does not give a noticeable performance increase, though.

EDIT2:
Here is a video of the polydraw tool showing drawing multipolygons and holes :slight_smile:

JAR Download: https://www.dropbox.com/s/lqt8xdkmg9euvq5/polydraw.jar?dl=0

#18

I’d disagree here, you are close to doing it already, I’ll write you some untested code


public boolean pointInPolygons(float x, float y) {
        // check bounding sphere first
        float dx = (x - centerX);
        float dy = (y - centerY);
        if (dx * dx + dy * dy > radiusSquared)
            return null;
        // check bounding box next
        if (maxX < x || maxY < y || minX > x || minY > y)
            return null;
        // ask interval tree for all polygon edges intersecting 'y'
        int c = tree.traverse(y, intervals, 0);
        // check the polygon edges
        int windingNumber = 0;
        for (int r = 0; r < c; r++) {
            Interval ival = intervals[r];
            float[] verticesXY = polygons[ival.polygonIndex];
            // order i and j in the order they are drawn in (assume verticesXY is ordered correctly)
            int i = Math.min(ival.i,ival.j);
            int j = Math.max(ival.i,ival.j);
            float y1 = verticesXY[2 * i + 1];
            float y2 = verticesXY[2 * j + 1];
            float x1 = verticesXY[2 * i + 0];
            float x2 = verticesXY[2 * j + 0];
            //no need to test intervals which are definitely on the wrong side of the point to be tested
            if (x1 && x2 < x)
            {
                continue;
            }
            // test that the interval is an upward interval
            // test also that y is to the left of this interval in the direction of travel (last logic condition)
            if ((y1 < y) && (y2 >= y) && (( (x2 - x1) * (y - y1)- (x - x1) * (y2 - y1) )>0))
            {
                windingNumber++;
            }
            // test that the interval is an downward interval
            // test also that y is to the right of this interval in the direction of travel (last logic condition)
            else if ((y2 < y) && (y1 >= y) && (( (x2 - x1) * (y - y1)- (x - x1) * (y2 - y1) )<0))
            {
                windingNumber--;
            }
        }
        return windingNumber !=0;
    }

If your code already filters out intervals that don’t need to be counted (wrong side of the point to be tested) then ignore the last logic statement in the 2 if’s. If you don’t then this could be an optimisation.

#19

Thanks! I’ll try that out.

In the meantime I added another feature I wanted to have for the physics in the soon-to-come 2D polygon game, which was already possible with a little extension to the crosscutting/even/odd algorithm.

-> if using multiple polygons, detect in which polygon the point is inside of

Here is a gifv: http://i.imgur.com/9azqeUZ.gifv

Was really easy to implement and can actually support millions of polygons very efficiently. If one had 1 million polygons drawn it would only take 122KB of RAM.

#20

what about making an object that holds the polygon as triangle pieces, rather than a polygon? This even gives me an idea for a data structure…holy cow. I’ll work on it later.