# Is Point Inside Polygon?

**QUESTION:** How can I determine if a given point is inside or outside
a polygon?

** ANSWER:** There are certainly several different ways to do this in IDL.
This solution is provided by Baard Krane (bard.krane@fys.uio.no) at the University of Oslo. He writes:

The code is based on an algorithm which takes the dot and cross product of vectors from the point itself to neighbouring points in the polygon. If we write the expressions for the dot and cross product, we see that by dividing the two we get the tangent of the angle between two vectors extending from the point to two neighbouring points in the polygon. A sum over all angles should then give a very small number if the point is outside the polygon. Rather than choosing an arbitrarily small number to check, I just choose to compare the total of the angles to Pi, which works adequately for my purposes.

The code (a program named **Inside **) is shown
below.

FUNCTION Inside, x, y, px, py ; x - The x coordinate of the point. ; y - The y coordinate of the point. ; px - The x coordinates of the polygon. ; py - The y coordinates of the polygon. ; ; The return value of the function is 1 if the point is inside the ; polygon and 0 if it is outside the polygon. sx = Size(px) sy = Size(py) IF (sx[0] EQ 1) THEN NX=sx[1] ELSE RETURN, -1 ; Error if px not a vector IF (sy[0] EQ 1) THEN NY=sy[1] ELSE RETURN, -1 ; Error if py not a vector IF (NX EQ NY) THEN N = NX ELSE RETURN, -1 ; Incompatible dimensions tmp_px = [px, px[0]] ; Close Polygon in x tmp_py = [py, py[0]] ; Close Polygon in y i = indgen(N) ; Counter (0:NX-1) ip = indgen(N)+1 ; Counter (1:nx) X1 = tmp_px(i) - x Y1 = tmp_py(i) - y X2 = tmp_px(ip) - x Y2 = tmp_py(ip) - y dp = X1*X2 + Y1*Y2 ; Dot-product cp = X1*Y2 - Y1*X2 ; Cross-product theta = Atan(cp,dp) IF (Abs(Total(theta)) GT !PI) THEN RETURN, 1 ELSE RETURN, 0 END

You can do a quick check of the program by typing the following commands:

IDL> x = 0.5 IDL> y = 0.5 IDL> px = [0.0, 1.0, 1.0, 0.0] IDL> py = [0.0, 0.0, 1.0, 1.0] IDL> Print, Inside(x,y,px,py) 1

If you are using a later version of IDL (later than
IDL 5.2, I think, or maybe 5.3), you can use the **IDLanROI**
object to find out if a point is inside a polygon. Use the
**ContainsPoints** method, like this:

IDL> object = Obj_New('IDLanROI', px, py) IDL> Print, object->ContainsPoints(x, y) 1

This method is more robust than the **Inside** method,
Since you can check a vector or points, and you can find out if
each point is inside, is outside, is on the line, or is a vertex point all
at once. See the **IDLanROI** documentation for details.

Be sure to destroy your objects when you are finished with them!

IDL> Obj_Destroy, object

# Updates

In early November 2003, William Connelly submitted
this vectorized version
of the **Inside** program discussed above.

FUNCTION Inside, x, y, px, py, Index=index ; Purpose: see if point is inside polygon ; Category: maths ; Input: x, y - [vector of] points ; px,py - points defining polygon (will be closed automatically) ; Output: vector of 1's and 0's ; OR ; indicies of points inside (if /Index is set) ; Author: "Bård Krane"; Mods: wmc - make it work with x, y as vectors ; More-info: posted to comp.lang.idl-pvwave on Wed, 01 Apr 1998 12:26:38 +0200 ; See-also: http://www.ecse.rpi.edu/Homepages/wrf/geom/pnpoly.html for another method, possibly better ; This is better than my routine "is_inside" ; Note: reduce test from 1e-8 to 1e-4, since we are usually in single precision. In fact, "0.1" ; would do as well. On_Error, 1 sx = Size(px) sy = Size(py) IF (sx[0] EQ 1) THEN NX = sx[1] ELSE Message,'Variable px is not a vector' IF (sy[0] EQ 1) THEN NY = sy[1] ELSE Message,'Varialbe py is not a vector' IF (NX EQ NY) THEN N = NX ELSE Message,'Incompatible vector dimensions' tmp_px = [px, px[0]] ; Close Polygon in x tmp_py = [py, py[0]] ; Close Polygon in y i = Indgen(N,/Long) ; indices 0...N-1 ip = Indgen(N,/Long) + 1 ; indices 1...N nn = N_Elements(x) X1 = tmp_px(i) # Replicate(1,nn) - Replicate(1,n) # Reform([x],nn) Y1 = tmp_py(i) # Replicate(1,nn) - Replicate(1,n) # Reform([y],nn) X2 = tmp_px(ip) # Replicate(1,nn) - Replicate(1,n) # Reform([x],nn) Y2 = tmp_py(ip) # Replicate(1,nn) - Replicate(1,n) # Reform([y],nn) dp = X2*X1 + Y1*Y2 ; Dot-product cp = X1*Y2 - Y1*X2 ; Cross-product theta = Atan(cp,dp) ret = Replicate(0L, N_Elements(x)) i = Where(Abs(Total(theta,1)) GT 0.01,count) IF (count GT 0) THEN ret(i)=1 IF (N_Elements(ret) EQ 1) THEN ret=ret[0] IF (Keyword_Set(index)) THEN ret=(Indgen(/Long, N_Elements(x)))(Where(ret eq 1)) RETURN, ret END

Copyright © 1998-2001 David W. Fanning

Last Updated 23 October 2001