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

354,246

Solution 1

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 work 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 box around your polygon. This is very easy, fast and can already save 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.

// 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 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 whether the polygon can have holes or will always be 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 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.

Demonstrating how the ray cuts through a polygon

The winding number algorithm would 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 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 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.

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 because they are too short to ever meet each other.

// 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 (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.

Solution 2

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;
}

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:

I run a semi-infinite ray horizontally (increasing x, fixed y) out from the test point, and count how many edges it crosses. At each crossing, the ray switches between inside and outside. This is called the Jordan curve theorem.

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.

Solution 3

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.

public bool IsPointInPolygon( Point p, Point[] polygon )
{
    double minX = polygon[ 0 ].X;
    double maxX = polygon[ 0 ].X;
    double minY = polygon[ 0 ].Y;
    double maxY = polygon[ 0 ].Y;
    for ( int i = 1 ; i < polygon.Length ; i++ )
    {
        Point q = polygon[ i ];
        minX = Math.Min( q.X, minX );
        maxX = Math.Max( q.X, maxX );
        minY = Math.Min( q.Y, minY );
        maxY = Math.Max( q.Y, maxY );
    }

    if ( p.X < minX || p.X > maxX || p.Y < minY || p.Y > maxY )
    {
        return false;
    }

    // https://wrf.ecse.rpi.edu/Research/Short_Notes/pnpoly.html
    bool inside = false;
    for ( int i = 0, j = polygon.Length - 1 ; i < polygon.Length ; j = i++ )
    {
        if ( ( polygon[ i ].Y > p.Y ) != ( polygon[ j ].Y > p.Y ) &&
             p.X < ( polygon[ j ].X - polygon[ i ].X ) * ( p.Y - polygon[ i ].Y ) / ( polygon[ j ].Y - polygon[ i ].Y ) + polygon[ i ].X )
        {
            inside = !inside;
        }
    }

    return inside;
}

Solution 4

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

function pointIsInPoly(p, polygon) {
    var isInside = false;
    var minX = polygon[0].x, maxX = polygon[0].x;
    var minY = polygon[0].y, maxY = polygon[0].y;
    for (var n = 1; n < polygon.length; n++) {
        var q = polygon[n];
        minX = Math.min(q.x, minX);
        maxX = Math.max(q.x, maxX);
        minY = Math.min(q.y, minY);
        maxY = Math.max(q.y, maxY);
    }

    if (p.x < minX || p.x > maxX || p.y < minY || p.y > maxY) {
        return false;
    }

    var i = 0, j = polygon.length - 1;
    for (i, j; i < polygon.length; j = i++) {
        if ( (polygon[i].y > p.y) != (polygon[j].y > p.y) &&
                p.x < (polygon[j].x - polygon[i].x) * (p.y - polygon[i].y) / (polygon[j].y - polygon[i].y) + polygon[i].x ) {
            isInside = !isInside;
        }
    }

    return isInside;
}

Solution 5

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.

Share:
354,246
Scott Evernden
Author by

Scott Evernden

40+ years programming everything. Virtual and shared online environments is my current focus. User-interface design and human-facing applications has been my main specialty for well over 30 years contributing core IP to several startups. I have specific experience in 2D, 3D, VR and AR applications over a wide array of projects for several young startups dating back to days at MIT's Media Lab. In the past, I've been involved with unique work utilizing Papervision, Alternativa, Torque and Unity. Most recently I've made many exploratory 3d environments with WebGL, and shared environment experiments built utilizing nodejs, a-frame, and reactjs also. Today I am creating Infinite Stores - a spatial CMS platform and editing tools for making dynamic content-rich 3D and virtual reality shopping experiences. Infinite Stores provides individuals as well as groups of artisans a new way to sell products online via virtual store-spaces within an endless shopping universe that they can maintain and visually merchandise - just like real-world stores. The layout and mgmt apps are super easy to use and run in any web browser and do not require any special equipment. Engaging store experiences can be enjoyed by anyone who loves to shop - via client apps that run on the web, on mobile, and immersively via head-mounted VR systems.

Updated on July 08, 2022

Comments

  • Scott Evernden
    Scott Evernden almost 2 years

    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.

  • Scott Evernden
    Scott Evernden over 15 years
    What if I don't happen to have the bounding box? :)
  • Richard T
    Richard T over 15 years
    You can easily create a bounding box - it's just the four points which use the greatest and least x and greatest and least y. More soon.
  • Faisal Al-Harbi
    Faisal Al-Harbi over 15 years
    mathematically elegant, but takes a bit of trig, making it painfully slow.
  • Faisal Al-Harbi
    Faisal Al-Harbi over 15 years
    very fast, and can be applied to more general shapes. back around 1990, we had "curvigons" where some sides were circular arcs. By analyzing those sides into circular wedges and a pair of triangles joining the wedge to the origin (polygon centroid), hit testing was fast and easy.
  • shoosh
    shoosh over 15 years
    got any references on these curvigons?
  • Jasper Bekkers
    Jasper Bekkers almost 15 years
    @DarenW: Only one acos per vertex; on the other hand, this algorithm should be the easiest to implement in SIMD because it has absolutely no branching.
  • Emilio M Bumachar
    Emilio M Bumachar over 14 years
    This doesn't work at all. Suppose the polygon is a triangle and the point is far below it (vertical distance much larger than both horizontal distance and the size of the triangle). Then the angle between p and all tree apices is aprox 90 deg, and the sum is aprox 270 deg. This is not a pathological example, it's just easy to visualize. If the point is inside, the sum is indeed 360 deg, but if it is outside the sum might coincidentally be 360 deg
  • David Segonds
    David Segonds over 14 years
    @emilio, if the point is far away from the triangle, I don't see how the angle formed by the point and two apices of the triangle will be 90 degrees.
  • iecanfly
    iecanfly over 14 years
    I would love a ref for the curvigons too.
  • Gavin
    Gavin over 14 years
    +1 Fantastic answer. On using the hardware to do it, I've read in other places that it can be slow because you have to get data back from the graphics card. But I do like the principle of taking load off the CPU a lot. Does anyone have any good references for how this might be done in OpenGL?
  • RMorrisey
    RMorrisey over 14 years
    If I'm not mistaken, the ray casting algorithm you describe only works for convex polygons, unless you consider a concave section to be a "hole" in the polygon
  • Jared Updike
    Jared Updike over 14 years
    +1 because this is so simple! The main problem is if your polygon and test point line up on a grid (not uncommon) then you have to deal with "duplicitous" intersections, for example, straight through a polygon point! (yielding a parity of two instead of one). Gets into this weird area: stackoverflow.com/questions/2255842/… . Computer Graphics is full of these special cases: simple in theory, hairy in practice.
  • JOM
    JOM about 14 years
    First use bounding box check to solve "point is far" cases. For trig, you could use pregenerated tables.
  • Mecki
    Mecki about 13 years
    @RMorrisey: Why do you think so? I fail to see how it would fail for a concave polygon. The ray may leave and re-enter the polygon multiple times when the polygon is concave, but in the end, the hit counter will be odd if the point is within and even if it is outside, also for concave polygons.
  • S P
    S P over 12 years
    The 'Fast Winding Number Algorithm', described at softsurfer.com/Archive/algorithm_0103/algorithm_0103.htm works pretty fast...
  • ThomasW
    ThomasW over 12 years
    Note that if you are really doing it in Cocoa, then you can use the [NSBezierPath containsPoint:] method.
  • Mick
    Mick over 12 years
    Don't you have to watch out for divide by zero if verty[j]==verty[i]?
  • aljabear
    aljabear about 12 years
    the post links to a page which discusses further details, including the divide-by-zero issue.
  • Peter Wood
    Peter Wood about 12 years
    @Mick It checks that verty[i] and verty[j] are either side of testy, so they are never equal.
  • gongzhitaao
    gongzhitaao over 11 years
    There seems to be a bug when the polygon is concave.
  • Nicolas Miari
    Nicolas Miari about 11 years
    I have only ONE polygon, which is guaranteed to be smaller than the screen of the iPhone/iPad. Perhaps I should use the GPU approach?
  • Mecki
    Mecki about 11 years
    @NicolasMiari For a single polygon, I guess the GPU approach will be slower, since for raw pixel access to a drawing buffer, the whole buffer must first be copied in memory, which is a pretty slow operation on iOS devices; unless you want to test thousands of points to be inside or outside your polygon in which case the faster test after the copy may pay off in the end.
  • Nicolas Miari
    Nicolas Miari about 11 years
    I ended up using vectors. A bit tough to debug, bit pays off.
  • bobobobo
    bobobobo about 11 years
    This actually works. I created a github gist to test it (GLUT).
  • hguser
    hguser about 11 years
    Any body can re-factor this method to accept parameters like this check([testpointx,testpointy],[p1x,p1y,p2x,p2y....])?
  • Philipp Lenssen
    Philipp Lenssen almost 11 years
    Works great, thanks, I converted to JavaScript. stackoverflow.com/questions/217578/…
  • Lukas Hanacek
    Lukas Hanacek almost 11 years
    this works in most of the cases but it is wrong and doesnt work properly always ! use the solution from M Katz which is correct
  • Dhara
    Dhara almost 11 years
    How can we find Number of vertices in the polygon?
  • Ash
    Ash almost 11 years
    Your usage of / to separate x and y coordinates is confusing, it reads as x divided by y. It's much more clear to use x, y (i.e. x comma y) Overall, a useful answer.
  • Jack Franzen
    Jack Franzen over 10 years
    does this use the MacMartin test as described in the artice bobobo cited? I can hardly read it it's so condensed and good looking. It was the fastest by far (70% faster than non-macmartin ray casting) when dealing with large polygons so I'm trying to figure out if this covers that
  • Ed Fryed
    Ed Fryed over 10 years
    Just used the "hardware" method to implement this using js on a canvas. Worked a charm. So mush easier to understand than all this maths stuff ;) Plus no gpu used as im not using hardware acceleration. And on top of that its really easy to debug as i can just unhide the second black and white canvas. Genius.
  • Cody Piersall
    Cody Piersall over 10 years
    Relevant Wikipedia link: en.wikipedia.org/wiki/Point_in_polygon. One of the pictures is from this article.
  • kaka
    kaka over 10 years
    @DarenW: According to wikipedia it doesn't have to be slow. You can simply check the quadrants instead of calculating the angles.
  • paulm
    paulm over 10 years
    Would 3 tests of a point VS an OBB using SAT be faster than using this method by combing all 3 OOB's as concave polygon?
  • Mecki
    Mecki over 10 years
    @paulm: I don't know. Maybe. I'd implement both and benchmark them. This will not only tell you which one is faster but also if one of them is significant faster or not.
  • Mikola
    Mikola about 10 years
    This code is not robust, and I would not recommend using it. Here is a link giving some detailed analysis: www-ma2.upc.es/geoc/Schirra-pointPolygon.pdf
  • Zev Eisenberg
    Zev Eisenberg about 10 years
    Of course, in Objective-C, CGPathContainsPoint() is your friend.
  • Jon
    Jon about 10 years
    @ZevEisenberg but that would be too easy! Thanks for the note. I'll dig up that project at some point to see why I used a custom solution. I likely didn't know about CGPathContainsPoint()
  • Maury Markowitz
    Maury Markowitz almost 10 years
    I realize this is old, I'm not sure if anyone will see it, but... David, what is the "oriented sum of angles"? Is this simply the angle between me and the point in question, 0..360? If so, don't you need to consider the number of points in the poly? Isn't it 360 only for four-point polys?
  • David Segonds
    David Segonds over 9 years
    @MauryMarkowitz, oriented in this context means that a counter-clockwise angle will be positive while a clockwise angle will be negative. This algorithm works regardless of the number of points in the polygon.
  • Alexander Pacha
    Alexander Pacha over 9 years
    This approach actually DOES have limitations (edge-cases): Checking the Point (15,20) in the Polygon [(10,10),(10,20),(20,20),(20,10)] will return false instead of true. Same with (10,20) or (20,15). For all other cases, the algorithm works perfectly fine and the false-negatives in edge-cases are ok for my application.
  • wardw
    wardw about 9 years
    @Alexander, this is in fact by design: by handling left & bottom boundaries in the opposite sense to top & right boundaries, should two distinct polygons share an edge, any point along this edge will be located in one and only one polygon. ..a useful property.
  • James Brown
    James Brown over 8 years
    This is >1000x faster than using GraphicsPath.IsVisible!! The bounding box check makes the function about 70% slower.
  • user3344977
    user3344977 over 8 years
    @o0'. Not sure what you mean? Under the "License to Use" section, it basically gives you full rights to the code.
  • Adam
    Adam over 8 years
  • t2k32316
    t2k32316 over 8 years
    stackoverflow.com/a/23223947/229311 provides a very nice explanation of this approach.
  • M Arfan
    M Arfan over 8 years
    This is very nice solution if coordinates draw in polygon. What about if coordinates line passed from polygon. It's not worked for that. Any have have idea it's also alert us if coordinate line pass in other drawn polygon. Thanks
  • Dominic Cerisano
    Dominic Cerisano almost 8 years
    This is the optimal solution, since it is O(n), and requires minimal computation. Works for all polygons. I researched this solution 30 years ago at my first IBM job. They signed off on it, and still use it today in their GIS technologies.
  • Sphinxxx
    Sphinxxx almost 8 years
    @Sathvik - It looks like that's the improved Winding Number algorithm (mentioned on Wikipedia, which uses signed crossings instead of trigonometric functions, and is just as fast as Ray Casting.
  • Rufus
    Rufus over 7 years
    It should be noted that there are some caveats when handling polygons with edges parallel to the direction of ray casting
  • Mecki
    Mecki over 7 years
    @Woofas I did, that's the case named "COLLINEAR" and I talked about it in the 3rd paragraph from bottom.
  • Nicolas Castro
    Nicolas Castro about 7 years
    This is the best answer. Tested in PHP and worked very fast!
  • Yves Daoust
    Yves Daoust over 6 years
    You missed an important solution for the case of convex polygons: by comparing the point against a diagonal, you can halve the number of vertices. And repeating this process, you reduce to a single triangle in Log(N) operations rather than N.
  • Ivan P.
    Ivan P. almost 6 years
    Thanks for insight. I have implemented the approach in yellow highlighted box, but I have a problem when a sweeping line comes straight through the angle created by intersection of edges.
  • Ivan P.
    Ivan P. almost 6 years
    (in addition to prev. comment) Example: polygon is formed by set of points: P => [(4;1), (2;4), (4;7), (6;4)] . If we test point A(4;4), we find that ray moves through point (2;4) (or through (6;4) ) intersecting two edges at once at polygon vertex, which counts as two intesections
  • amitp
    amitp over 5 years
    In the comments here, Mikola says the code is not robust and provides a URL which is now a broken link. The paper is Schirra S. (2008) How Reliable Are Practical Point-in-Polygon Strategies? DOI doi.org/10.1007/978-3-540-87744-8_62 and I found a mirror at cs.tau.ac.il/~danha/courses/CG2008/pointPolygon.pdf . Mikola also has a Javascript library robust-point-in-polygon.
  • spurra
    spurra over 5 years
    The link to the source of the code is broken. The page can be found under the following link: wrf.ecse.rpi.edu//Research/Short_Notes/pnpoly.html
  • tiritea
    tiritea over 5 years
    "...avoiding possible errors of wrapping around when one crossed line 180 of longitude and when handling polar regions." can you perhaps describe this in more detail? I think I can figure out how to move everything 'up' to avoid my polygon crossing 0 longitude, but I'm not clear on how to handle polygons that contain either of the poles...
  • tetrahydrocannon
    tetrahydrocannon over 4 years
    Not only GraphicsPath.IsVisible() is way slower but also it doesn't work well with very small polygons with side in the 0.01f range
  • anorskdev
    anorskdev over 4 years
    This has a potential divide by zero problem in calculating b. Need to only calculate "b" and the next line with "c" if the calculation for "a" shows that there is a possibility of intersecting. No possibility if both points are above, or both points below - which is described by the calculation for "a".
  • Celdor
    Celdor about 4 years
    I’d like to use this method in an non-profit open source code. Can I do that without breaking law?
  • Liam Clink
    Liam Clink about 4 years
    I implemented this, and have been doing pathological case testing, and I've determined that this algorithm needs to take vertex intersection into account. If you intersect a vertex, it counts as intersecting two edges.
  • Mecki
    Mecki about 4 years
    @LiamClink Vertex intersection of different polygons play no role as you need to test each polygon separately with this algorithm. And the vertices or a regular polygon cannot intersect with each other. If you refer to a complex polygon, as the one shown to the very right at cutt.ly/sysGTuT please note that there are in fact three regular polygons that could be tested individually. Personally I would not allow complex polygons in software that simulates real world effects as they are unnatural (a shape like this cannot exist in real life).
  • psksvp
    psksvp about 4 years
    I upvote this. It works for the cases that I need. The algorithm is easy to understand and simple. psksvp
  • shieldgenerator7
    shieldgenerator7 about 4 years
    If you always cast the ray from the left or right, then the line is always in the form of y=b, where b is the y coordinate of your point. This simplifies a lot of the calculations
  • GDavoli
    GDavoli almost 4 years
    What's the point of the first for and if? The last for works just fine for all cases.
  • M Katz
    M Katz almost 4 years
    @GDavoli It's an efficiency thing. If the point is not inside the bounding box of the polygon. it can't be in the polygon. This is based on the assumption that the loop to find the bounding box of the polygon is significantly faster than the second for loop. That might not be true on modern hardware. But in a real system it may make sense to cache the bounding box of each polygon, in which case the bounding box check definitely makes sense.
  • hypehuman
    hypehuman almost 4 years
    for the COLLINEAR case... I think I have a solution. If the point is on that edge (simple to test), then we're on the edge of the polygon, return a special value that causes the main method to return true immediately. Otherwise, our ray skips right over the segment, which counts as entering and leaving (or leaving and entering), so it's like two intersections. But you don't actually need to add those two intersections to the total, since it wouldn't flip the parity and therefore doesn't affect the final answer. So you could just return NO.
  • Techiemanu
    Techiemanu over 3 years
    Do not give me correct output while checking the point exists or not in polygon. Not all cases are covered
  • R. Navega
    R. Navega over 3 years
    Another name for the ray casting algorithm is "Even-Odd", I think the latter is more popular.
  • Liam Clink
    Liam Clink over 3 years
    @Mecki I mean when testing within an individual polygon. Every vertex in a polygon has two edges, so if your ray intersects a vertex, it intersects two edges. And then there’s two cases, either the ray crossed inside and it should be counted as one intersection, or it just touched the vertex from the outside, where it should be counted as 0 or 2 intersections.
  • Mecki
    Mecki over 3 years
    @LiamClink There is no "touching". Only lines may touch (tangent) circles. A ray, which is a line, can only cross a vertex, never touch a vertex at all or be colinear to it. The case you are talking about is the colinear case (the entire vertex lies on the ray or the ray stops while being on the vertex). The code detects that case and leaves it up to you how to handle it. If the point you are looking for lies on the polygon side, is it inside the polygon or not? Does the side belong to the inside or to the outside of the polygon? You decide that yourself.
  • Liam Clink
    Liam Clink over 3 years
    @Mecki How can you say that a ray cannot pass through a point? That is what I mean by touching. The issue is whether the ray passes through the vertex while crossing the perimeter or not crossing, and I'm not sure how to handle discriminating these two cases
  • Mecki
    Mecki over 3 years
    @LiamClink The algorithm I presented doesn't test if a the ray hits points/vertices, it tests if a ray crosses connections ("lines") between vertices. And here the ray either crosses it, doesn't cross it or is collinear, there is no fourth case. You will not be able to feed any data into my test code that will not report one of these three cases. And two cases are obvious to handle and for the third one, I already explained in my post that you are free to handle this case as you like.
  • Liam Clink
    Liam Clink over 3 years
    Right and what I am saying is that if a ray goes through a vertex, that is equivalent to crossing all edges connected to it as best I can tell.
  • Ahmad
    Ahmad over 3 years
    This is the best answer. all the other answers ignore corner cases.
  • aferriss
    aferriss over 2 years
    I was also pretty confused by the / and - notation used for the sides. I think it would be more clearly expressed as: side1: [(x1, y1), (x2, y2)] side2: [(x2, y2), (x3, y3)]
  • MaVCArt
    MaVCArt over 2 years
    This worked for me better than any other algorithm I found
  • AndrejH
    AndrejH over 2 years
    Fastest and handles the edge cases!
  • GMSL
    GMSL over 2 years
    For the colinear case, as it's very rare you could probably just choose another nearby point and try again with that, e.g. change point tracing ray from to ((x - 0.01), (y - 0.01)), check this new point is still outside the bounding box, and run the function from there. As you've changed x and y it won'tbe colinear this time, and the added computing time is probably minimal.
  • BIS Tech
    BIS Tech over 2 years
    @MKatz Why didn't you break the 2nd loop?
  • M Katz
    M Katz over 2 years
    Break after inside = !inside;? No can do. The algorithm changes that multiple times as it hits the various polygon segments.