;+
; :Description:
;    Generates a planar representation of a convex polygon.  These planes can then be
;    used to determine if points are inside or outside the polygon
;    Multiple polygons can be passed in to this routine.  See testpolygons.pro for 
;    an example
; :Params:
;    x - x vertices
;    y - y vertices
;    polygons - connectivity of the vertices
;
; :Author: Ronn Kling and Jerry Lefever
;-
function pointInsideApolygon, x,y,polygons, data, count=count

sizeOfPolygon =  size(polygons,/dimen)
numEdges = sizeOfPolygon[0]
numPolygons = sizeOfPolygon[1]
;there is a plane for every edge in every polygon
planarSpace = fltarr(numEdges,3,numPolygons)

xVertArray = fltarr(numEdges, numPolygons)
yVertArray = fltarr(numEdges, numPolygons)
;get the vertices for each edge
for i=0,numEdges-1 do begin
  xVertArray[i,*] = x[polygons[i,*]]
  yVertArray[i,*] = y[polygons[i,*]]
endfor
  
aArray = fltarr(numEdges, numPolygons)
bArray = fltarr(numEdges, numPolygons)
cArray = fltarr(numEdges, numPolygons)

;get the constants for each plane defined by an edge
for i=0L, numEdges-1 do begin
  uInd = (i+1) mod numEdges
  edgeLength = sqrt( (xVertArray[uInd,*] - xVertArray[i,*])^2  + $
                     (yVertArray[uInd,*] - yVertArray[i,*])^2)
  aArray[i,*] = (yVertArray[uInd,*] - yVertArray[i,*])/edgeLength
  bArray[i,*] = -(xVertArray[uInd,*] - xVertArray[i,*])/edgeLength
  cArray[i,*] = -aArray[i,*]*xVertArray[i,*] - bArray[i,*]*yVertArray[i,*]  
endfor

;get the center of the polygon
xc = total(xVertArray,1)/numEdges
yc = total(yVertArray,1)/numEdges

;check the normal and switch signs if necessary
for i=0L,numEdges-1 do begin
  ind = where((xc*aArray[i,*] + yc*bArray[i,*] + cArray[i,*]) lt 0, nCount)
  if nCount gt 0 then begin
    aArray[i,ind] = -aArray[i,ind]
    bArray[i,ind] = -bArray[i,ind]
    cArray[i,ind] = -cArray[i,ind]
  endif
  planarSpace[i,*,*] = [aArray[i,*], bArray[i,*], cArray[i,*]]
endfor 

x = reform(data[0,*])
y = reform(data[1,*])

origIndices = lonarr(n_elements(x))

for i=0,n_elements(x)-1 do begin
  for j=0,numPolygons-1 do begin
    plane = planarSpace[*,*,j]
    Q = min( plane # [x[i],y[i],1])
    if Q ge 0 then begin ;means that the point is in the polygon
      origIndices[i] = 1
      break
    endif
  endfor
endfor  

goodIndices = where( origIndices eq 1, count)

;for i=0L,numEdges-1 do begin
;  posInd = where((x*aArray[i,*] + y*bArray[i,*] + cArray[i,*]) gt 0, nCount)
;  if nCount eq 0 then continue
;  aArray[i,*] = aArray[i,posInd]
;  x = x[posInd]
;  y = y[posInd]
;  z = z[posInd]
;  origIndices = origIndices[posInd]
;endfor
  
;  for j=0,numPolygons-1 do begin
;    plane = planarSpace[*,*,j]
;    Q = min( plane # [x,y,1])
;    if Q ge 0 then begin ;means that the point is in the polygon
;      plots, x,y, psym=3,color=4
;      break
;    endif
;  endfor

return,goodIndices

end