Source code for inpoly.pro:

function inpoly,Px,Py,xl,yl

; Aug 2004 - Written by Alexie Leithauld

;----------------------------------------------------------------
;   inpoly.pro 
;
;   Determines if a point P(px, py) in inside or outside a polygon.
;   The method used is the ray intersection method based on the
;   Jordan curve theorem. Basically, we trace a "ray" (a semi line)
;   from the point P and we calculate the number of edges that it 
;   intersects. If this number is odd, P is inside.
;    
;   (Px,Py)    Coordinates of point
;   xv         List of x coordinates of vertices of polygon
;   yv         List of y coordinates of vertices of polygon
;
;   The x and y lists may be given clockwise or anticlockwise.
;   The algorithm works for both convex and concave polygons.
;----------------------------------------------------------------

;---------------------- Parameters  -----------------------------

N=n_elements(xl)
xv=dblarr(N)
yv=dblarr(N)

for i=0.,N-1 do begin
 xv[i]=xl[i]
 yv[i]=yl[i]
endfor

nc=0                     ;Number of edge crossings
N=n_elements(xv);        ;Number of vertices

; Test input
if (N lt 3) then print,'A polygon must have at least three vertices'
if (n_elements(xv) ne n_elements(yv)) then print, 'Must have same number of X and Y coordinates'

;---------------------- Change coordinate system -----------------

; Place P at the center of the coordinate system.

for i=0,N-1 do begin 
    xv[i]=xv[i]-Px
    yv[i]=yv[i]-Py  
endfor

;---------------------- Calculate crossings ----------------------

; The positive half of the x axis is chosen as the ray
; We must determine how many edges cross the x axis with x>0

for i=0,N-1 do begin

Ax=xv[i]      ; First vertice of edge
Ay=yv[i]

if (i eq N-1) then begin
  Bx=xv[0]
  By=yv[0]
endif else begin
  Bx=xv[i+1]    ; Second vertice of edge
  By=yv[i+1]
endelse

; We define two regions in the plan: R1/ y<0 and R2/ y>=0. Note that
; the case y=0 (vertice on the ray) is included in R2.

if (Ay lt 0) then signA = -1 else signA=+1
if (By lt 0) then signB = -1 else signB=+1

; The edge crosses the ray only if A and B are in different regions.
; If a vertice is only the ray, the number of crossings will still be
; correct.

if( (signA*signB) lt 0) then begin

    ; if Ax>0 and Bx>o then the edge crosses the positive x axis
    if((Ax gt 0) and (Bx gt 0)) then begin
    nc=nc+1

    ; Otherwise (end points are in diagonally opposite quadrants)
    ; we must calculate the intersection
    endif else begin
	x=Ax-(Ay*(Bx-Ax))/(By-Ay)
        if (x gt 0) then nc=nc+1
    endelse
endif
endfor  ;i

;if inside then uneven
;if outside then even

nc=nc mod 2
return,nc


end

Return to the shapelets web page or the code help menu.