Home

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

  1. The C Code
  2. The Method
  3. Originality
  4. The Inequality Tests are Tricky
  5. C Semantics
  6. Point on a (Boundary) Edge
  7. Concave Components, Multiple Components, and Holes
  8. Testing Which One of Many Polygons Contains the Point
  9. Explanation of "for (i = 0, j = nvert-1; i < nvert; j = i++)"
  10. Fortran Code for the Point in Polygon Test
  11. Converting the Code to All Integers
  12. License to Use
  13. Inclusion Testing by Subtended Angles
  14. Convex Polygons, such as Rectangles
  15. Almost Convex Polygons
  16. 3D Polygons
  17. Testing Point Inclusion in a Polyhedron
  18. If You Have a Suggestion or Find an Error
  19. My Other Work
  20. 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 = 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;
}
ArgumentMeaning
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 y-coordinates of the polygon's vertices.
testx, testyX- and y-coordinate of the test point.

The Method

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 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:

  1. Is the point in the half-plane to the left of the extended edge? and
  2. Is the point's Y coordinate within the edge's Y-range?
Handling endpoints here is tricky.

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 publicly-distributed 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 point-in-polygon 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:

A well-written 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 divide-by-zero, 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.

  1. 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, TL and TR. Depending on internal roundoff errors, PNPOLY may say that P is in TL or in TR. However it will always give the same answer when P is tested against those triangles. That is, if PNPOLY finds that P is in TL, then it will find that P is not TR. If PNPOLY finds that P is not in TL, then it will find that P is in TR.

  2. 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.

  3. 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.

  4. 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 semi-open 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

  1. 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.

  2. The direction that you list the vertices (clockwise or counterclockwise) does not matter.

  3. 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.

    1. First, include a (0,0) vertex.

    2. Then include the first component' vertices, repeating its first vertex after the last vertex.

    3. Include another (0,0) vertex.

    4. Include another component or hole, repeating its first vertex after the last vertex.

    5. Repeat the above two steps for each component and hole.

    6. Include a final (0,0) vertex.

  4. 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.

  5. Each component or hole's vertices may be listed either clockwise or counter-clockwise.

  6. 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 semi-infinite 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 = nvert-1; i < nvert; j = i++)"

The intention is to execute the loop for each i from 0 to nvert-1. For each iteration, j is i-1. However that wraps, so if i=0 then j=nvert-1. Therefore the current edge runs between verts j and i, and the loop is done once per edge. In detail:

  1. Start by setting i and j:
    i = 0
    j = nvert-1
  2. If i<nvert is false then exit the loop.
  3. Do the loop body.
  4. Set j=i and then
    add 1 to i and then
  5. Go back to step 2.

Fortran Code for the Point in Polygon Test

Here it is, in Fortran; I wrote it 47 years ago.

C>>>PNP1 C C .................................................................. C C SUBROUTINE PNPOLY C C PURPOSE C TO DETERMINE WHETHER A POINT IS INSIDE A POLYGON C C USAGE C CALL PNPOLY (PX, PY, XX, YY, N, INOUT ) C C DESCRIPTION OF THE PARAMETERS C PX - X-COORDINATE OF POINT IN QUESTION. C PY - Y-COORDINATE OF POINT IN QUESTION. C XX - N LONG VECTOR CONTAINING X-COORDINATES OF C VERTICES OF POLYGON. C YY - N LONG VECTOR CONTAING Y-COORDINATES OF C VERTICES OF POLYGON. C N - NUMBER OF VERTICES IN THE POLYGON. C INOUT - THE SIGNAL RETURNED: C -1 IF THE POINT IS OUTSIDE OF THE POLYGON, C 0 IF THE POINT IS ON AN EDGE OR AT A VERTEX, C 1 IF THE POINT IS INSIDE OF THE POLYGON. C C REMARKS C THE VERTICES MAY BE LISTED CLOCKWISE OR ANTICLOCKWISE. C THE FIRST MAY OPTIONALLY BE REPEATED, IF SO N MAY C OPTIONALLY BE INCREASED BY 1. C THE INPUT POLYGON MAY BE A COMPOUND POLYGON CONSISTING C OF SEVERAL SEPARATE SUBPOLYGONS. IF SO, THE FIRST VERTEX C OF EACH SUBPOLYGON MUST BE REPEATED, AND WHEN CALCULATING C N, THESE FIRST VERTICES MUST BE COUNTED TWICE. C INOUT IS THE ONLY PARAMETER WHOSE VALUE IS CHANGED. C THE SIZE OF THE ARRAYS MUST BE INCREASED IF N > MAXDIM C WRITTEN BY RANDOLPH FRANKLIN, UNIVERSITY OF OTTAWA, 7/70. C C SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED C NONE C C METHOD C A VERTICAL LINE IS DRAWN THRU THE POINT IN QUESTION. IF IT C CROSSES THE POLYGON AN ODD NUMBER OF TIMES, THEN THE C POINT IS INSIDE OF THE POLYGON. C C .................................................................. C SUBROUTINE PNPOLY(PX,PY,XX,YY,N,INOUT) REAL X(200),Y(200),XX(N),YY(N) LOGICAL MX,MY,NX,NY INTEGER O C OUTPUT UNIT FOR PRINTED MESSAGES DATA O/6/ MAXDIM=200 IF(N.LE.MAXDIM)GO TO 6 WRITE(O,7) 7 FORMAT('0WARNING:',I5,' TOO GREAT FOR THIS VERSION OF PNPOLY. 1RESULTS INVALID') RETURN 6 DO 1 I=1,N X(I)=XX(I)-PX 1 Y(I)=YY(I)-PY INOUT=-1 DO 2 I=1,N J=1+MOD(I,N) MX=X(I).GE.0.0 NX=X(J).GE.0.0 MY=Y(I).GE.0.0 NY=Y(J).GE.0.0 IF(.NOT.((MY.OR.NY).AND.(MX.OR.NX)).OR.(MX.AND.NX)) GO TO 2 IF(.NOT.(MY.AND.NY.AND.(MX.OR.NX).AND..NOT.(MX.AND.NX))) GO TO 3 INOUT=-INOUT GO TO 2 3 IF((Y(I)*X(J)-X(I)*Y(J))/(X(J)-X(I))) 2,4,5 4 INOUT=0 RETURN 5 INOUT=-INOUT 2 CONTINUE RETURN END

Converting the Code to All Integers

If you want to convert the code from floats to integers, consider these issues.

  1. On many current processors floats are at least as fast as ints.
  2. If you move the denominator over to the other side of the inequality, remember that, when the denominator is negative, the inequality will flip.
  3. If coordinates are large enough, the multiplication will silently overflow.

License to Use

Copyright (c) 1970-2003, 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:

  1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimers.
  2. Redistributions in binary form must reproduce the above copyright notice in the documentation and/or other materials provided with the distribution.
  3. 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.

  1. Join the test point to each edge in turn
  2. Calculate the angle that the point subtends at the edge.
  3. Note that the edge must be directed (that is, you know the starting point from the ending point), and this angle must be signed.
  4. Add the angles.
  5. The result must be either 0 or 2pi.
  6. The test point is inside iff the result is 2pi.
This method is obvious but not good.
  1. The simple implementation requires an arctan evaluation for each edge, which is slow.
  2. The arctan can be avoided with a careful use of the arctan sum formula, but this is also tricky.
  3. Properly handling the case where the edge crosses the positive X-axis is tricky.
  4. 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.

  1. Find the equation of the infinite line that contains each edge.

  2. Express each equation as an expression: d=ax+by+c. Points on the line will give zero.

  3. 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.

  4. Now, test your test point against every line. It's inside the polygon iff it's on the inside side of every line.

  5. Note that the equation is particularly simple for horizontal and vertical edges.

Almost Convex Polygons

  1. What if you like testing your point against the line equations, but, unfortunately, your polygon is not convex?

  2. It is possible to express the interior of any polygon as a boolean expression in the half-planes 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.

  3. Then, the Boolean expression for the polygon's interior is A.B.C.D, that is, the four line expressions anded together.

  4. 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.

  5. Then, the Boolean expression for the polygon's interior is this

    A.B.C.(D+E+F)

  6. A similar formula exists for any 2D polygon.

  7. However, not all 3D polyhedra have such formulae. This is yet another difference between 2D and 3D. For more differences, see here.

3D Polygons

  1. By this, I mean a flat polygon embedded in a plane in 3D.

  2. Simply project the polygon's vertices, and the test point, into 2D.

  3. Do the projection by deleting one coordinate, e.g., project (x,y,z) to (x,z).

  4. 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.

  1. Run a semi-infinite ray up from the point, P, and
  2. Count how many faces it crosses.
The ray crosses face F iff
  1. P lies below the plane of F, and
  2. when P is projected (in the direction of the ray) onto the plane of F, then it is inside F. This is a 2D point-containment problem.
If the ray intersects a vertex or an edge, then you're in trouble. Consider perturbing P slightly. Simulation of Simplicity provides a better, more complicated, solution.

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 8-line 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