PNPOLY  Point Inclusion in Polygon Test
W. Randolph Franklin (WRF)
This is an expansion of the answer in the
comp.graphics.algorithms FAQ
question 2.03, How do I find if a point lies within a polygon?
Table of Contents
 The C Code
 The Method
 Originality
 The Inequality Tests are Tricky
 C Semantics
 Point on a (Boundary) Edge
 Concave Components, Multiple Components, and Holes
 Testing Which One of Many Polygons Contains the Point
 Explanation of "for (i = 0, j = nvert1; i < nvert; j = i++)"
 Fortran Code for the Point in Polygon Test
 Converting the Code to All Integers
 License to Use
 Inclusion Testing by Subtended Angles
 Convex Polygons, such as Rectangles
 Almost Convex Polygons
 3D Polygons
 Testing Point Inclusion in a Polyhedron
 If You Have a Suggestion or Find an Error
 My Other Work
 Acknowledgements
The C Code
Here is the code, for reference. Excluding lines with only braces, there are only 7 lines of code.
int pnpoly(int nvert, float *vertx, float *verty, float testx, float testy) { int i, j, c = 0; for (i = 0, j = nvert1; i < nvert; j = i++) { if ( ((verty[i]>testy) != (verty[j]>testy)) && (testx < (vertx[j]vertx[i]) * (testyverty[i]) / (verty[j]verty[i]) + vertx[i]) ) c = !c; } return c; }
Argument  Meaning 

nvert  Number of vertices in the polygon. Whether to repeat the first vertex at the end is discussed below. 
vertx, verty  Arrays containing the x and ycoordinates of the polygon's vertices. 
testx, testy  X and ycoordinate of the test point. 
The Method
I run a semiinfinite 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 case of the ray going thru a vertex is handled correctly via a careful selection of inequalities. Don't mess with this code unless you're familiar with the idea of Simulation of Simplicity. This pretends to shift the ray infinitesimally down so that it either clearly intersects, or clearly doesn't touch. Since this is merely a conceptual, infinitesimal, shift, it never creates an intersection that didn't exist before, and never destroys an intersection that clearly existed before.
The ray is tested against each edge thus:
 Is the point in the halfplane to the left of the extended edge? and
 Is the point's Y coordinate within the edge's Yrange?
Originality
I make no claim to having invented the idea. However in 1970, I did produce the Fortran code given below on my own, and include it in a package of cartographic SW publiclydistributed by David Douglas, Dept of Geography, Simon Fraser U and U of Ottawa.
The C code in the FAQ is my code, lightly edited.
Earlier implementations of pointinpolygon testing presumably exist, tho the code might never have been released. Pointers to prior art, especially publicly available code, are welcome. One early publication, which doesn't handle the point on an edge, and has a typo, is this:
 M Shimrat, "Algorithm 112, Position of Point Relative to Polygon", Comm. ACM 5(8), Aug 1962, p 434.
A wellwritten recent summary is this:
The Inequality Tests are Tricky
If translating the program to another language, be sure to get the inequalities in the conditional correct. They were carefully chosen to make the program work correctly when the point is vertically below a vertex.
Several people have thought that my program was wrong, when really they had gotten the inequalities wrong.
C Semantics
My code uses the fact that, in the C language, when executing the
code a&&b
, if a
is false, then
b
must not be evaluated. If your compiler doesn't
do this, then it's not implementing C, and you will get a
dividebyzero, i.a., when the test point is vertically in line with a
vertical edge. When translating this code to another language with
different semantics, then you must implement this test explicitly.
Point on a (Boundary) Edge
PNPOLY partitions the plane into points inside the polygon and points outside the polygon. Points that are on the boundary are classified as either inside or outside.
Any particular point is always classified consistently the same way. In the following figure, consider what PNPOLY would say when the red point, P, is tested against the two triangles, T_{L} and T_{R}. Depending on internal roundoff errors, PNPOLY may say that P is in T_{L} or in T_{R}. However it will always give the same answer when P is tested against those triangles. That is, if PNPOLY finds that P is in T_{L}, then it will find that P is not T_{R}. If PNPOLY finds that P is not in T_{L}, then it will find that P is in T_{R}.
If you want to know when a point is exactly on the boundary, you need another program. This is only one of many functions that PNPOLY lacks; it also doesn't predict tomorrow's weather. You are free to extend PNPOLY's source code.
The first reason for this is the numerical analysis position that you should not be testing exact equality unless your input is exact. Even then, computational roundoff error would often make the result wrong.
The second reason is that, if you partition a region of the plane into polygons, i.e., form a planar graph, then PNPOLY will locate each point into exactly one polygon. In other words, PNPOLY considers each polygon to be topologically a semiopen set. This makes things simpler, i.e., causes fewer special cases, if you use PNPOLY as part of a larger system. Examples of this include locating a point in a planar graph, and intersecting two planar graphs.
Concave Components, Multiple Components, and Holes
The polygon may be concave. However, if a vertex is very close to an edge (that the vertex is not an end of) then beware of roundoff errors.
The direction that you list the vertices (clockwise or counterclockwise) does not matter.
The polygon may contain multiple separate components, and/or holes, which may be concave, provided that you separate the components and holes with a (0,0) vertex, as follows.
First, include a (0,0) vertex.
Then include the first component' vertices, repeating its first vertex after the last vertex.
Include another (0,0) vertex.
Include another component or hole, repeating its first vertex after the last vertex.
Repeat the above two steps for each component and hole.
Include a final (0,0) vertex.
For example, let three components' vertices be A1, A2, A3, B1, B2, B3, and C1, C2, C3. Let two holes be H1, H2, H3, and I1, I2, I3. Let O be the point (0,0). List the vertices thus:
O, A1, A2, A3, A1, O, B1, B2, B3, B1, O, C1, C2, C3, C1, O, H1, H2, H3, H1, O, I1, I2, I3, I1, O.
Each component or hole's vertices may be listed either clockwise or counterclockwise.

If there is only one connected component, then it is optional to repeat the first vertex at the end. It's also optional to surround the component with zero vertices.
Testing Which One of Many Polygons Contains the Point
This is called the point location problem in computational geometry, see [Preparata]. There are many methods with different requirements for preprocessing time and storage, versus query time. If the polygons form a planar graph with each edge labelled with its two polygon neighbors, then one simple way is to run a semiinfinite ray up from the point until it hits its first edge. Then the relevant containing polygon is one of that edge's neighbors.
Explanation of "for (i = 0, j = nvert1; i < nvert; j = i++)"
The intention is to execute the loop for each i from 0 to nvert1. For each iteration, j is i1. However that wraps, so if i=0 then j=nvert1. Therefore the current edge runs between verts j and i, and the loop is done once per edge. In detail:
 Start by setting i and j:
i = 0
j = nvert1  If i<nvert is false then exit the loop.
 Do the loop body.
 Set j=i and then
add 1 to i and then  Go back to step 2.
Fortran Code for the Point in Polygon Test
Here it is, in Fortran; I wrote it 47 years ago.
Converting the Code to All Integers
If you want to convert the code from floats to integers, consider these issues.
 On many current processors floats are at least as fast as ints.
 If you move the denominator over to the other side of the inequality, remember that, when the denominator is negative, the inequality will flip.
 If coordinates are large enough, the multiplication will silently overflow.
License to Use
Copyright (c) 19702003, Wm. Randolph Franklin
Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:
 Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimers.
 Redistributions in binary form must reproduce the above copyright notice in the documentation and/or other materials provided with the distribution.
 The name of W. Randolph Franklin may not be used to endorse or promote products derived from this Software without specific prior written permission.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
Inclusion Testing by Subtended Angles
This is another method of testing whether a point is inside a polygon.
 Join the test point to each edge in turn
 Calculate the angle that the point subtends at the edge.
 Note that the edge must be directed (that is, you know the starting point from the ending point), and this angle must be signed.
 Add the angles.
 The result must be either 0 or 2pi.
 The test point is inside iff the result is 2pi.
 The simple implementation requires an arctan evaluation for each edge, which is slow.
 The arctan can be avoided with a careful use of the arctan sum formula, but this is also tricky.
 Properly handling the case where the edge crosses the positive Xaxis is tricky.
 Finally, if you handle the above problems correctly and optimize the algorithm, then this method reduces to the Jordan curve formula.
Convex Polygons, such as Rectangles
If the polygon is convex, then an alternate method is this.
Find the equation of the infinite line that contains each edge.
Express each equation as an expression: d=ax+by+c. Points on the line will give zero.
Standardize each equation so that if a point inside the polygon is substituted in, the result is positive, or equivalently an outside point will be negative.
Now, test your test point against every line. It's inside the polygon iff it's on the inside side of every line.
Note that the equation is particularly simple for horizontal and vertical edges.
Almost Convex Polygons
What if you like testing your point against the line equations, but, unfortunately, your polygon is not convex?
It is possible to express the interior of any polygon as a boolean expression in the halfplanes defined by the edges. For the convex polygon below,
let A be a Boolean expression that is true if the test point is on the inside side of the infinite line containing this edge of the polygon. Ditto B, C, D.
Then, the Boolean expression for the polygon's interior is A.B.C.D, that is, the four line expressions anded together.
Now consider the following concave polygon,
As before, A, B, C, D, E, F are Boolean expressions for the test point being on the inside side of each edge in turn.
Then, the Boolean expression for the polygon's interior is this
A.B.C.(D+E+F) A similar formula exists for any 2D polygon.
However, not all 3D polyhedra have such formulae. This is yet another difference between 2D and 3D. For more differences, see here.
3D Polygons
By this, I mean a flat polygon embedded in a plane in 3D.
Simply project the polygon's vertices, and the test point, into 2D.
Do the projection by deleting one coordinate, e.g., project (x,y,z) to (x,z).
For numerical stability, delete the coordinate with the smallest range.
Testing Point Inclusion in a Polyhedron
Here is one way to test whether a 3D point in inside a polyhedron.
 Run a semiinfinite ray up from the point, P, and
 Count how many faces it crosses.
 P lies below the plane of F, and
 when P is projected (in the direction of the ray) onto the plane of F, then it is inside F. This is a 2D pointcontainment problem.
If You Have a Suggestion or Find an Error
I'm always interested in error reports, altho it might take awhile to respond. So far, there have been some good suggestions. However, all but one of the error reports have themselves been erroneous. I've incorporated their errors in this file. It's truly amazing how much trouble an 8line program can cause!
My Other Work
Since 1970 my students and I have devised many other algorithms to efficiently process the largest geometric and GIS datasets on parallel computers. They include point location in planar graphs and in 3D tesslations, and overlaying (intersecting) planar graphs and 3D tesslations and computing areas of the intersections, and using rational numbers to avoid roundoff errors. They parallelize well and have been tested on very large datasets. Details are on my web site.Acknowledgements
Thanks to
 Loren Albertazzi,
for finding an error in the treatment of multiple component polygons.  Mark Sullivan for suggesting the documentation change from "vertical" to "horizontal".
 Christopher Sargent for suggesting an optimization based on Joe O'Rourke's Computational Geometry in C.