Point in Polygon Algorithm

Question

I saw the below algorithm works to check if a point is in a given polygon from this link:

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

I tried this algorithm and it actually works just perfect. But sadly I cannot understand it well after spending some time trying to get the idea of it.

So if someone is able to understand this algorithm, please explain it to me a little.

Thank you.

The answer you link to contains a link to its source, and the source has an extensive explanation of the algorithm: http://www.ecse.rpi.edu/Homepages/wrf/Research/Short_Notes/pnpoly.html Did you read that? — sleske, Mar 18 '13 at 22:54
I've found http://erich.realtimerendering.com/ptinpoly/ to be an excellent resource for understanding "point in polygon" strategies. — Jon Kirkman, Apr 25 '13 at 01:59

Chowlett · Accepted Answer · 2013-08-10T21:17:16.923

52

The algorithm is ray-casting to the right. Each iteration of the loop, the test point is checked against one of the polygon's edges. The first line of the if-test succeeds if the point's y-coord is within the edge's scope. The second line checks whether the test point is to the left of the line (I think - I haven't got any scrap paper to hand to check). If that is true the line drawn rightwards from the test point crosses that edge.

By repeatedly inverting the value of c, the algorithm counts how many times the rightward line crosses the polygon. If it crosses an odd number of times, then the point is inside; if an even number, the point is outside.

I would have concerns with a) the accuracy of floating-point arithmetic, and b) the effects of having a horizontal edge, or a test point with the same y-coord as a vertex, though.

edited Aug 10 '13 at 21:17

answered Jul 30 '12 at 06:31

Chowlett

44,716
18
112
146

1

Would the ordering of the vertices matter for this test? – paulm Jan 27 '14 at 23:33
@paulm - yes; in the code above, the polygon is formed by joining vertex 0 to vertex 1 to vertex 2... to vertex 0. – Chowlett Jan 28 '14 at 08:55
Sorry - I mean would the winding ordering matter? – paulm Jan 28 '14 at 09:28
@paulm - I'm afraid I don't follow; what do you mean by the "winding ordering"? – Chowlett Jan 28 '14 at 09:56
If the polygon verts are ordered clockwise or anti clockwise v0 -> v1 -> v3 or v3 -> v1 -> v0? – paulm Jan 28 '14 at 12:32
1

@paulm - Ah, I see. It... shouldn't do. I haven't checked the algorithm in depth to make sure, but each step only cares about one edge, at which point clockwise/anti-clockwise has no meaning. – Chowlett Jan 28 '14 at 13:18

Josh · Answer 2 · 2022-01-30T21:22:47.300

24

Edit 1/30/2022: I wrote this answer 9 years ago when I was in college. People in the chat conversation are indicating it's not accurate. You should probably look elsewhere. ‍♂️

Chowlett is correct in every way, shape, and form. The algorithm assumes that if your point is on the line of the polygon, then that is outside - for some cases, this is false. Changing the two '>' operators to '>=' and changing '<' to '<=' will fix that.

bool PointInPolygon(Point point, Polygon polygon) {
  vector<Point> points = polygon.getPoints();
  int i, j, nvert = points.size();
  bool c = false;
  
  for(i = 0, j = nvert - 1; i < nvert; j = i++) {
    if( ( (points[i].y >= point.y ) != (points[j].y >= point.y) ) &&
        (point.x <= (points[j].x - points[i].x) * (point.y - points[i].y) / (points[j].y - points[i].y) + points[i].x)
      )
      c = !c;
  }
  
  return c;
}

edited Jan 30 '22 at 21:22

answered Mar 24 '13 at 14:10

Josh

1,848
13
19

Is there a more efficient way of doing the Point in Polygon? MySQL spatial extension takes about 50ms to query a table of 100 polygons.. this function in PHP takes about 100ms. – kouton Mar 29 '13 at 08:22
1

@Reddox MySQL probably uses native (compiled) code, while PHP is interpreted, and not known for its speed... – PhiLho Jul 08 '13 at 12:15
Actually (without actually researching it much), I'd say that a simple optimization is first checking for the bounding box (if it's out of the bounding box no need to do the check that's more refined). – Fabio Zadrozny Jan 06 '14 at 16:12
@Josh, thanks for the fix. There is just a minor brace issue in the if condition, `(points[i].y) >= point.y` should be `(points[i].y >= point.y)` – Michal Fapso Mar 04 '14 at 09:20
2

@Josh, You answer is confusing. There is only one '' operators to '>=' as well as the **single** ' 200,100 -> 200,200 -> 100,200 fails. – Ash Apr 02 '14 at 03:40
7

This answer is not correct despite the up-voting. By having point.x <= ... means that you consider a point on the edge as having a ray that crosses that edge to the right. This will count points that fall on rightmost edges as "in" and points that fall on edges to the left as "out". The obvious correction would be to explicitly test for the equality and exit the loop. – Ecuador Aug 19 '14 at 16:49
does the order of the points in the points vector matter to the outcome of this test? – BenKoshy Nov 11 '15 at 07:04
Now if you are coincident with a vertex you double count that vertex and are on the outside, while in the polygon. And fixing that will often screw up the testcase where, an internal vertex comes down and touches the line coincident with what you are checking. Without actually being on the edge, just in the way of raycasting but not putting you outside the polygon. V_____ – Tatarize Apr 26 '16 at 00:18
@Ecuador is very much correct. It is best to check to see if the point lies on the edge of the polygon in the loop and simply return true in that case. – Borgboy Sep 15 '16 at 12:21
2

Like @Ecuador says, this is not a correct answer. The original **pnpoly** algorithm, and the point on an edge case, is described [here](https://wrf.ecse.rpi.edu//Research/Short_Notes/pnpoly.html) by its author. TL;DR, if the point is on the edge between two polygons, the **pnpoly** algorithm will _consistently_ classify the point as inside one polygon, and outside the other polygon. – alondono May 04 '17 at 01:59

apil.tamang · Answer 3 · 2014-03-27T14:09:57.610

This might be as detailed as it might get for explaining the ray-tracing algorithm in actual code. It might not be optimized but that must always come after a complete grasp of the system.

    //method to check if a Coordinate is located in a polygon
public boolean checkIsInPolygon(ArrayList<Coordinate> poly){
    //this method uses the ray tracing algorithm to determine if the point is in the polygon
    int nPoints=poly.size();
    int j=-999;
    int i=-999;
    boolean locatedInPolygon=false;
    for(i=0;i<(nPoints);i++){
        //repeat loop for all sets of points
        if(i==(nPoints-1)){
            //if i is the last vertex, let j be the first vertex
            j= 0;
        }else{
            //for all-else, let j=(i+1)th vertex
            j=i+1;
        }

        float vertY_i= (float)poly.get(i).getY();
        float vertX_i= (float)poly.get(i).getX();
        float vertY_j= (float)poly.get(j).getY();
        float vertX_j= (float)poly.get(j).getX();
        float testX  = (float)this.getX();
        float testY  = (float)this.getY();

        // following statement checks if testPoint.Y is below Y-coord of i-th vertex
        boolean belowLowY=vertY_i>testY;
        // following statement checks if testPoint.Y is below Y-coord of i+1-th vertex
        boolean belowHighY=vertY_j>testY;

        /* following statement is true if testPoint.Y satisfies either (only one is possible) 
        -->(i).Y < testPoint.Y < (i+1).Y        OR  
        -->(i).Y > testPoint.Y > (i+1).Y

        (Note)
        Both of the conditions indicate that a point is located within the edges of the Y-th coordinate
        of the (i)-th and the (i+1)- th vertices of the polygon. If neither of the above
        conditions is satisfied, then it is assured that a semi-infinite horizontal line draw 
        to the right from the testpoint will NOT cross the line that connects vertices i and i+1 
        of the polygon
        */
        boolean withinYsEdges= belowLowY != belowHighY;

        if( withinYsEdges){
            // this is the slope of the line that connects vertices i and i+1 of the polygon
            float slopeOfLine   = ( vertX_j-vertX_i )/ (vertY_j-vertY_i) ;

            // this looks up the x-coord of a point lying on the above line, given its y-coord
            float pointOnLine   = ( slopeOfLine* (testY - vertY_i) )+vertX_i;

            //checks to see if x-coord of testPoint is smaller than the point on the line with the same y-coord
            boolean isLeftToLine= testX < pointOnLine;

            if(isLeftToLine){
                //this statement changes true to false (and vice-versa)
                locatedInPolygon= !locatedInPolygon;
            }//end if (isLeftToLine)
        }//end if (withinYsEdges
    }

    return locatedInPolygon;
}

Just one word about optimization: It isn't true that the shortest (and/or the tersest) code is the fastest implemented. It is a much faster process to read and store an element from an array and use it (possibly) many times within the execution of the block of code than to access the array each time it is required. This is especially significant if the array is extremely large. In my opinion, by storing each term of an array in a well-named variable, it is also easier to assess its purpose and thus form a much more readable code. Just my two cents...

score 7 · Answer 4 · answered May 10 '17 at 15:32

I changed the original code to make it a little more readable (also this uses Eigen). The algorithm is identical.

// This uses the ray-casting algorithm to decide whether the point is inside
// the given polygon. See https://en.wikipedia.org/wiki/Point_in_polygon#Ray_casting_algorithm
bool pnpoly(const Eigen::MatrixX2d &poly, float x, float y)
{
    // If we never cross any lines we're inside.
    bool inside = false;

    // Loop through all the edges.
    for (int i = 0; i < poly.rows(); ++i)
    {
        // i is the index of the first vertex, j is the next one.
        // The original code uses a too-clever trick for this.
        int j = (i + 1) % poly.rows();

        // The vertices of the edge we are checking.
        double xp0 = poly(i, 0);
        double yp0 = poly(i, 1);
        double xp1 = poly(j, 0);
        double yp1 = poly(j, 1);

        // Check whether the edge intersects a line from (-inf,y) to (x,y).

        // First check if the line crosses the horizontal line at y in either direction.
        if ((yp0 <= y) && (yp1 > y) || (yp1 <= y) && (yp0 > y))
        {
            // If so, get the point where it crosses that line. This is a simple solution
            // to a linear equation. Note that we can't get a division by zero here -
            // if yp1 == yp0 then the above if be false.
            double cross = (xp1 - xp0) * (y - yp0) / (yp1 - yp0) + xp0;

            // Finally check if it crosses to the left of our test point. You could equally
            // do right and it should give the same result.
            if (cross < x)
                inside = !inside;
        }
    }
    return inside;
}

you seem to be the only one mentioning the zero division case ;) — Baedsch, Dec 13 '18 at 16:02

score 3 · Answer 5 · answered Dec 18 '13 at 11:35

The algorithm is stripped down to the most necessary elements. After it was developed and tested all unnecessary stuff has been removed. As result you can't undertand it easily but it does the job and also in very good performance.

I took the liberty to translate it to ActionScript-3:

// not optimized yet (nvert could be left out)
public static function pnpoly(nvert: int, vertx: Array, verty: Array, x: Number, y: Number): Boolean
{
    var i: int, j: int;
    var c: Boolean = false;
    for (i = 0, j = nvert - 1; i < nvert; j = i++)
    {
        if (((verty[i] > y) != (verty[j] > y)) && (x < (vertx[j] - vertx[i]) * (y - verty[i]) / (verty[j] - verty[i]) + vertx[i]))
            c = !c;
    }
    return c;
}

score 3 · Answer 6 · edited Mar 21 '14 at 22:20

This algorithm works in any closed polygon as long as the polygon's sides don't cross. Triangle, pentagon, square, even a very curvy piecewise-linear rubber band that doesn't cross itself.

1) Define your polygon as a directed group of vectors. By this it is meant that every side of the polygon is described by a vector that goes from vertex an to vertex an+1. The vectors are so directed so that the head of one touches the tail of the next until the last vector touches the tail of the first.

2) Select the point to test inside or outside of the polygon.

3) For each vector Vn along the perimeter of the polygon find vector Dn that starts on the test point and ends at the tail of Vn. Calculate the vector Cn defined as DnXVn/DN*VN (X indicates cross product; * indicates dot product). Call the magnitude of Cn by the name Mn.

4) Add all Mn and call this quantity K.

5) If K is zero, the point is outside the polygon.

6) If K is not zero, the point is inside the polygon.

Theoretically, a point lying ON the edge of the polygon will produce an undefined result.

The geometrical meaning of K is the total angle that the flea sitting on our test point "saw" the ant walking at the edge of the polygon walk to the left minus the angle walked to the right. In a closed circuit, the ant ends where it started. Outside of the polygon, regardless of location, the answer is zero.
Inside of the polygon, regardless of location, the answer is "one time around the point".

score 2 · Answer 7 · answered Jul 30 '12 at 07:17

This method check whether the ray from the point (testx, testy) to O (0,0) cut the sides of the polygon or not .

There's a well-known conclusion here: if a ray from 1 point and cut the sides of a polygon for a odd time, that point will belong to the polygon, otherwise that point will be outside the polygon.

score 2 · Answer 8 · edited May 23 '17 at 12:10

To expand on @chowlette's answer where the second line checks if the point is to the left of the line, No derivation is given but this is what I worked out: First it helps to imagine 2 basic cases:

the point is left of the line . / or
the point is right of the line / .

If our point were to shoot a ray out horizontally where would it strike the line segment. Is our point to the left or right of it? Inside or out? We know its y coordinate because it's by definition the same as the point. What would the x coordinate be?

Take your traditional line formula y = mx + b. m is the rise over the run. Here, instead we are trying to find the x coordinate of the point on that line segment that has the same height (y) as our point.

So we solve for x: x = (y - b)/m. m is rise over run, so this becomes run over rise or (yj - yi)/(xj - xi) becomes (xj - xi)/(yj - yi). b is the offset from origin. If we assume yi as the base for our coordinate system, b becomes yi. Our point testy is our input, subtracting yi turns the whole formula into an offset from yi.

We now have (xj - xi)/(yj - yi) or 1/m times y or (testy - yi): (xj - xi)(testy - yi)/(yj - yi) but testx isn't based to yi so we add it back in order to compare the two ( or zero testx as well )

score 1 · Answer 9 · answered Jul 30 '12 at 06:30

I think the basic idea is to calculate vectors from the point, one per edge of the polygon. If vector crosses one edge, then the point is within the polygon. By concave polygons if it crosses an odd number of edges it is inside as well (disclaimer: although not sure if it works for all concave polygons).

score 1 · Answer 10 · answered Jan 01 '14 at 13:05

This is the algorithm I use, but I added a bit of preprocessing trickery to speed it up. My polygons have ~1000 edges and they don't change, but I need to look up whether the cursor is inside one on every mouse move.

I basically split the height of the bounding rectangle to equal length intervals and for each of these intervals I compile the list of edges that lie within/intersect with it.

When I need to look up a point, I can calculate - in O(1) time - which interval it is in and then I only need to test those edges that are in the interval's list.

I used 256 intervals and this reduced the number of edges I need to test to 2-10 instead of ~1000.

score 1 · Answer 11 · edited Feb 22 '16 at 08:09

Here's a php implementation of this:

<?php
class Point2D {

    public $x;
    public $y;

    function __construct($x, $y) {
        $this->x = $x;
        $this->y = $y;
    }

    function x() {
        return $this->x;
    }

    function y() {
        return $this->y;
    }

}

class Point {

    protected $vertices;

    function __construct($vertices) {

        $this->vertices = $vertices;
    }

    //Determines if the specified point is within the polygon. 
    function pointInPolygon($point) {
        /* @var $point Point2D */
    $poly_vertices = $this->vertices;
    $num_of_vertices = count($poly_vertices);

    $edge_error = 1.192092896e-07;
    $r = false;

    for ($i = 0, $j = $num_of_vertices - 1; $i < $num_of_vertices; $j = $i++) {
        /* @var $current_vertex_i Point2D */
        /* @var $current_vertex_j Point2D */
        $current_vertex_i = $poly_vertices[$i];
        $current_vertex_j = $poly_vertices[$j];

        if (abs($current_vertex_i->y - $current_vertex_j->y) <= $edge_error && abs($current_vertex_j->y - $point->y) <= $edge_error && ($current_vertex_i->x >= $point->x) != ($current_vertex_j->x >= $point->x)) {
            return true;
        }

        if ($current_vertex_i->y > $point->y != $current_vertex_j->y > $point->y) {
            $c = ($current_vertex_j->x - $current_vertex_i->x) * ($point->y - $current_vertex_i->y) / ($current_vertex_j->y - $current_vertex_i->y) + $current_vertex_i->x;

            if (abs($point->x - $c) <= $edge_error) {
                return true;
            }

            if ($point->x < $c) {
                $r = !$r;
            }
        }
    }

    return $r;
}

Test Run:

        <?php
        $vertices = array();

        array_push($vertices, new Point2D(120, 40));
        array_push($vertices, new Point2D(260, 40));
        array_push($vertices, new Point2D(45, 170));
        array_push($vertices, new Point2D(335, 170));
        array_push($vertices, new Point2D(120, 300));
        array_push($vertices, new Point2D(260, 300));


        $Point = new Point($vertices);
        $point_to_find = new Point2D(190, 170);
        $isPointInPolygon = $Point->pointInPolygon($point_to_find);
        echo $isPointInPolygon;
        var_dump($isPointInPolygon);

score 0 · Answer 12 · answered Jan 24 '15 at 21:30

I modified the code to check whether the point is in a polygon, including the point is on an edge.

bool point_in_polygon_check_edge(const vec<double, 2>& v, vec<double, 2> polygon[], int point_count, double edge_error = 1.192092896e-07f)
{
    const static int x = 0;
    const static int y = 1;
    int i, j;
    bool r = false;
    for (i = 0, j = point_count - 1; i < point_count; j = i++)
    {
        const vec<double, 2>& pi = polygon[i);
        const vec<double, 2>& pj = polygon[j];
        if (fabs(pi[y] - pj[y]) <= edge_error && fabs(pj[y] - v[y]) <= edge_error && (pi[x] >= v[x]) != (pj[x] >= v[x]))
        {
            return true;
        }

        if ((pi[y] > v[y]) != (pj[y] > v[y]))
        {
            double c = (pj[x] - pi[x]) * (v[y] - pi[y]) / (pj[y] - pi[y]) + pi[x];
            if (fabs(v[x] - c) <= edge_error)
            {
                return true;
            }
            if (v[x] < c)
            {
                r = !r;
            }
        }
    }
    return r;
}

Changes are simple, If the point is on the edge, Check it as the point is in the polygon. If pi[y] == pj[y] == v[y] and v[x] is in between pi[x] and pj[x], point is on the edge of horizontal line of the polygon. if c == v[x] and v[y] is in between pi[y] and pj[y], the point is on the edge of vertical line of the polygon. — rawdev, Jan 26 '15 at 02:30

Point in Polygon Algorithm

12 Answers12

Linked

Related