Point in Polygon Algorithm

Question

I saw that the algorithm below 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 really works just fine. But, unfortunately, I cannot understand it well, having spent some time trying to figure it out.

So, if someone can understand this algorithm, please explain this to me a little.

Thank.

Answer 1

The algorithm is beam-casting to the right. Each iteration of the loop, the control point is checked at one of the edges of the polygon. The first line of the if-test succeeds if the point-coordinate of the point is within the border area. The second line checks if the control point is to the left of the line (I think I don’t have scrap paper to check it). If so, a line to the right of the control point crosses this edge.

By inverting the value of c again, the algorithm calculates how many times the line to the right intersects the polygon. If it crosses an odd number of times, then the point is inside; if an even number, the point is out.

I would have problems with: a) the precision of floating point arithmetic, and b) the consequences of having a horizontal edge or control point with the same y-coordinate as the vertex.

Answer 2

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

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

Answer 3

This can be as detailed as possible to explain the ray tracing algorithm in real code. This may not be optimized, but it should always happen after a complete understanding 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 is not true that the shortest (and / or tersest) code is the fastest. This is a much faster process for reading and storing an element from an array and using it (possibly) many times in the process of executing a block of code than for accessing an array every time it is required. This is especially important if the array is extremely large. In my opinion, storing each member of the array in a well-defined variable, it is also easier to evaluate its purpose and, therefore, generate much more readable code. Only my two cents ...

Answer 4

The algorithm is divided into the most necessary elements. After it was developed and tested, all unnecessary things were deleted. As a result, you cannot easily handle this, but it does the job, as well as a very good job.

---

I took the liberty of translating 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; } 

Answer 5

This algorithm works in any closed polygon if the sides of the polygon do not intersect. Triangle, pentagon, square, even a very curved piecewise-linear rubber band that does not cross itself.

1) Define your polygon as a directed group of vectors. This implies that each side of the polygon is described by a vector that goes from the vertex a to the vertex a + 1. The vectors are so directed that the head of one touches the tail of the next until the last vector touches the tail of the first.

2) Select a point to check inside or outside the polygon.

3) For each vector Vn along the perimeter of the search vector of the polygon Dn, starting at the test point and ending at the tail of Vn. Calculate the vector Cn defined as DnXVn / DN * VN (X denotes the transverse product, * indicates the point product). Name the value Cn by the name Mn.

4) Add all Mn and call this value K.

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

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

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

The geometric meaning of K is the full angle at which the flea sitting on our test point “saw” ant, going along the edge of the polygon, going to the left minus the corner going to the right. In a closed loop, ant ends where it begins. Outside the polygon, regardless of location, the answer is zero.
Inside the polygon, regardless of location, the answer is "once around a point."

---

Answer 6

This method checks if the ray cut from the point (testx, testy) to O (0,0) of the side of the polygon or not.

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

Answer 7

I think the main idea is to compute vectors from a point, one at a time to the edge of the polygon. If the vector intersects one edge, then the point is inside the polygon. According to concave polygons, if it intersects an odd number of edges, it is also inside (disclaimer: although not sure if it works for all concave polygons).

Answer 8

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

Testing:

  <?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); 

Answer 9

Expand on @chowlette answer , where the second line checks to see if the point is to the left of the line. No output is given, but this is what I developed: First, this allows us to present two main cases:

• the dot remains from the line . / . / or
• dot to the right of the line / .

If our point was to shoot the rays horizontally, wherever he would hit a segment of the line. Is our point heading left or right? Inside or outside? We know its y coordinate, because by definition it is the same as a point. What would be the x coordinate?

Take the traditional line formula y = mx + b . m - rise above the span. Here, instead, we try to find the x coordinate of a point on this segment of a line that has the same height (y) as our point .

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

Now we have (xj - xi)/(yj - yi) or 1/m times y or (testy - yi) : (xj - xi)(testy - yi)/(yj - yi) , but testx is not based on yi , so we add it back to compare two (or zero testx as well)

Answer 10

This is the algorithm I use, but I added a bit of preprocessing workaround to speed it up. My polygons have ~ 1000 edges, and they do not change, but I need to see if the cursor is inside one at each mouse movement.

I basically split the height of the bounding box into equal length intervals, and for each of these intervals I compile a list of edges that lie inside / intersect with it.

When I need to find a point, I can calculate - in O (1) time - what interval it is in, and then I only need to check those edges that are in the list of intervals.

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

Answer 11

I changed the code to check if the point is in the polygon, including the point on the 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; } 

Answer 12

I modified the original code to make it more readable (it also 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; }