Inside or outside a polygonal region?

[Michael Bell]

Please can anyone help me with a reasonably straightforward algorithm for determining whether a point is inside or outside a region defined by a polygon? I have searched previous posts, but what I've found (e.g. on determining whether a mouse click is within a region) doesn't seem to offer the answer I need for my purposes.

The context is a spatial simulation model for a fishery. Individual fish move within a spatial arena, parts of which may be closed to fishing. If the closed areas are simple shapes, such as rectangles, circles or triangles, it is easy to determine whether a fish at a given coordinate (x,y) is inside or outside a closed area. Convex polygons are slightly less easy, but I could still see how to determine a simple rule for inside/outside. Where I'm stuck is with defining rules for more complicated closed areas, the edges of which may contain both convexities and concavities.

One possible approach to these more complex shapes is to approximate them on a fine-scale grid. Inside/outside determination is then just a matter of determining the identity of the current grid square and looking up the status of this square in an array. This might cause some edge effects, since (except for exact N-S or E-W directions) there will no longer be straight lines between the vertices of the polygon. Possibly this doesn't matter provided that the grid squares are small with respect to the scale of fish movements within any time step of the model, but I can see that this might give rise to some very large arrays.

So far I have just been working with areas defined as blocks on a grid, but I would be interested to know if there is a more elegant and exact solution that I could use for modelling more complex shapes. I suspect that this type of problem may be fairly trivial for someone with experience of GIS and mapping applications ...? For my purposes it needs to be fast to execute, since the program involves simulation of the movements of millions of individuals over, say, daily time steps for a hundred or more years (required to equilibrate for a given fishery management scenario). I won't be wanting to display the movements on the screen (except for some basic checking work to make sure I'm doing what I think I am).

Thanks in advance for any enlightenment. For information: I program with PB/CC5, but I'm by no means a programmer.

Cheers

Mike

[Manuel Valdes]

May be you can use colors to differentiate inside and outside. Colors have a very precise definition so, if you can get the object sorroundig color or colors, you can determine whether it is inside or outside, even being the naked eye unable to differentiate.

Regards,

[Michael Mattias]

CreatePolygonRgn () ==> Polygon
PtInRegion() true/false

[Michael Bell]

Thanks Manuel, I can see that that would help if I was plotting the simulation results on a graphics screen - then I could use the pixels as my grid. However, in this case I don't want to use graphics at all - I just want the simulations to run invisibly, and to collect various outcomes (such as the fish being caught or dying of natural causes) along the way.

What I need is a simple function for inside or outside. For example, for a rectangular closed area it would be trivial:

 
FUNCTION InsideRectangle(BYVAL x as EXT, BYVAL y as EXT, _ 
                         BYREF Closure() as EXT) as LONG
  '--------------------------------------------------------
  ' x and y are coordinates to be tested
  ' Closure(1 to 2, 1 to 2) gives min and max
  ' x and y coordinates for the closed rectangle
  '--------------------------------------------------------
 
  IF x>=Closure(1,1) AND x<=Closure(1,2) AND _
     y>=Closure(2,1) AND y<=Closure(2,2) THEN
    FUNCTION=-1 ' true
  ELSE
    FUNCTION=0  ' false
  END IF
 
END FUNCTION

So in my code I could ask ISTRUE(InsideRectangle(x,y,Closure())). What I really want to ask is ISTRUE(InsidePolygon(x,y,Closure(),npoints)), where Closure() gives the npoints coordinates of the vertices of any polygon defining a closed area. Possibly this is trivial as well, but I really can't see an answer at present.

[Manuel Valdes]

May be you can achieve the same purpose invisibly using bitmaps instead. It is possible to simulate the same on a bitmap, with the bonus that you will not be constrained by the screen limits. Elaborating a bit, once the poligonal is defined, you can fill it with a colour different from the background, if you get the pixel in a point and its colour is the background's then the point is outside else inside.

Regards,

[Michael Bell]

Thank you Michael. I was afraid that WinAPI might come into it. Does that, like Manuel's suggestion, depend on me working with pixels on a screen? Or else with virtual pixels that are never displayed? I might well go along with this as a good fine-scale grid approximation, but if possible I'd like to find some kind of geometry-based logic solution.

[Barry Marks]

You can do that in a bitmap in memory without putting it on a screen. Draw your polygon in a memory bitmap in a different color than the background, then check the color of the point.

I suspect there's a better mathmatical answer than this but it should be very workable and it should be fairly fast.

You might want to google that, also. I googled "is a point inside a polygon" and found a lot of papers describing algorithms.

http://en.wikipedia.org/wiki/Point_in_polygon

and http://alienryderflex.com/polygon/

were a couple of likely looking ones.

Barry

[Michael Bell]

Manuel, Barry - Thanks for the pointer towards bitmaps. Memory bitmaps should definitely do the job. Thanks also for the links - at least two algorithms there, and implementations in C, Fortran and other languages. I should be able to make something of these...

[Marc Giao]

Hi Michael,

I believe what you're trying to do is often referd to as "Hit Testing". A Google search for that phrase should keep you busy for awhile

[Michael Bell]

Thank you Marc - yes, plenty out there on hit testing!

Seems to boil down to the 'ray casting' algorithm as the best candidate: how many times does a ray starting at the test point pass through the boundary of the polygon? If the number of intersections is odd then the point is inside, if it is zero or even then the point is outside. Of course points at vertices or on the edges cause complications. Microsoft suggest this as a better option than PtInRegion() for hit testing in a polygon (http://support.microsoft.com/kb/121960), though I can't pretend to follow the code given here.

Thanks folks for much help. I'll go ahead with the ray casting, and possibly test against bitmaps for speed.

[Joseph Cote]

If you can break your polygons down into adjacent triangular regions then there are a number of methods to determine if a given point is inside or outside.

The most elegant, and the one I can't seem to put my finger on at the moment, would be the Barycentric method. One formula, if result is positive, point is inside; negative means outside; and zero would mean the point was on the edge. I'll look some more when I get home tonight. If your math is strong, have a look for yourself.

The brute force solution, but more understandable in a visual way, is the sum-of-subtriangles method. Given triangle ABC and test point D. Calculate the areas of triangles ABD, BCD and ACD, if they sum to the area of ABC then the point is inside or on the edge else outside. (Ob warning re calling real values "equal.") This works because if the test point is inside, then you're cutting the original triangle into component pieces; if it's outside, you're defining new pieces that will be bigger than the original.

Either way, maintain your polygons as a list of triangles, run through the list until you find a hit inside or reach the end of the list. Timesaver: pre-calculate the areas of your original triangles.

[Joseph Cote]

Here's another way to see if a point is inside a triangle. I was seeing if the origin was enclosed, but it's pretty obvious where to change it.

FUNCTION containsOrigin(corners() AS LONG) AS LONG
    LOCAL x1,y1, x2,y2, x3,y3, x0, y0 AS EXT
    LOCAL b0, b1, b2, b3 AS EXT
    x0 = 0 : y0 = 0  'point to see if inside triangle
    x1 = corners(1)
    y1 = corners(2)
    x2 = corners(3)
    y2 = corners(4)
    x3 = corners(5)
    y3 = corners(6)
    b0 =  (x2-x1)*(y3-y1) - (x3-x1)*(y2-y1)
    b1 = ((x2-x0)*(y3-y0) - (x3-x0)*(y2-y0)) / b0
    b2 = ((x3-x0)*(y1-y0) - (x1-x0)*(y3-y0)) / b0
    b3 = ((x1-x0)*(y2-y0) - (x2-x0)*(y1-y0)) / b0
    IF b1 > 0 AND b2 > 0 AND b3 > 0 THEN
        FUNCTION = -1 'TRUE
    ELSE
        FUNCTION = 0  'FALSE
    END IF
END FUNCTION

[Paul Dixon]

Michael,
if there are lots of checks to make on the same polygon(s) then the best way to do this is to draw the polygon(s) in a graphic window or bitmap that matches the shape of the real polygon.
It takes very little code and very little time to do and, in effect, all of the calculation is done just once to draw the polygons and the problem becomes one of simply looking up the result for each of your million fish.

You fill the polygon with a colour defining the region leaving the remainder of the graphic black and your checks can be done like this:

'PBCC5.01 program
TYPE PolyPoint
  x AS SINGLE
  y AS SINGLE
END TYPE

TYPE PolyArray
  COUNT AS LONG
  xy(1 TO 1024) AS PolyPoint
END TYPE

%TotalFish = 100
%xSize=500
%ySize = 500

FUNCTION PBMAIN () AS LONG
    
'define the location of all the fish
LOCAL Fishes() AS PolyPoint
DIM Fishes(%TotalFish) AS PolyPoint
RANDOMIZE TIMER
FOR WhichFish& = 1 TO %TotalFish
    Fishes(WhichFish&).x=RND(1,%xSize)
    Fishes(WhichFish&).y=RND(1,%ySize)
NEXT


CONSOLE SET LOC 500,0
           
GRAPHIC WINDOW "Graphics", 0, 0, 500, 500  TO hWin???
GRAPHIC ATTACH hWin???, 0, REDRAW
GRAPHIC CLEAR  0 'set background to black

'define the coordinates of the fisrt polygon
DIM MyPolygon AS PolyArray
MyPolygon.count=9   'specify 9 points in this polygon
MyPolygon.xy(1).x=20    :MyPolygon.xy(1).y=90
MyPolygon.xy(2).x=30    :MyPolygon.xy(2).y=130
MyPolygon.xy(3).x=55    :MyPolygon.xy(3).y=100
MyPolygon.xy(4).x=40    :MyPolygon.xy(4).y=420
MyPolygon.xy(5).x=100   :MyPolygon.xy(5).y=400
MyPolygon.xy(6).x=160   :MyPolygon.xy(6).y=300
MyPolygon.xy(7).x=380   :MyPolygon.xy(7).y=231
MyPolygon.xy(8).x=290   :MyPolygon.xy(8).y=59
MyPolygon.xy(9).x=20    :MyPolygon.xy(9).y=50

GRAPHIC POLYGON MyPolygon  ,%RED, %RED ' [, [fillstyle&] [, fillmode&]]]]

'define the coordinates of the second polygon
MyPolygon.count=5   'specify 5 points in this polygon
MyPolygon.xy(1).x=200    :MyPolygon.xy(1).y=390
MyPolygon.xy(2).x=230    :MyPolygon.xy(2).y=430
MyPolygon.xy(3).x=355    :MyPolygon.xy(3).y=400
MyPolygon.xy(4).x=400    :MyPolygon.xy(4).y=420
MyPolygon.xy(5).x=250    :MyPolygon.xy(5).y=300
      
GRAPHIC POLYGON MyPolygon  ,%BLUE, %BLUE ' [, [fillstyle&] [, fillmode&]]]]

FOR WhichFish& = 1 TO %TotalFish
    GRAPHIC GET PIXEL (Fishes(WhichFish&).x,Fishes(WhichFish&).y) TO WhichRegion&
    SELECT CASE WhichRegion&
        CASE %RED
            'the fish is in the red zone
            INCR RedCnt&
            
        CASE %BLUE
            'the fish is in the blue zone
            INCR BlueCnt&
            
        CASE ELSE
            'the fish is not in any zone
            INCR NoneCnt&
            
    END SELECT
    
    'note, don't do this next line in real life as 2 fish in the same location will count as 1 fish in the region plus a white one
    GRAPHIC SET PIXEL (Fishes(WhichFish&).x,Fishes(WhichFish&).y ), %WHITE

NEXT

GRAPHIC REDRAW

PRINT  "Fish in red zone :";RedCnt&
PRINT  "Fish in blue zone:";BlueCnt&
PRINT  "Fish in no zone  :";NoneCnt&
PRINT
PRINT "Press a key to exit"

CONSOLE SET FOCUS

WAITKEY$
GRAPHIC WINDOW END

END FUNCTION

You don't need to draw the window on the screen, you could just as well use a bitmap which is never drawn but I chose to draw it so you can see what the program does.
There are more precise ways of doing what you want but they involve a lot of calculation for each fish which may be verry slow for a million fish.

Paul.

[Paul Dixon]

Mike,
I don't have time to code an example now but the alternative which gives accurate results goes as follows.

For each polygon area
    For each fish
        consider a line from outside the polygon to the fish
        For each line making up the polygon, check if the line to the fish crosses the line of the polygon 
            If it crosses then increment COUNT
        next line in the polygon
        
        If COUNT is ODD then the fish is inside the polygon, if it's even then it's outside.

    Next Fish
Next Polygon area

For a fast line crossing calculation see here:

http://www.powerbasic.com/support/pbforums/showthread.php?t=23695

Paul.

[Michael Bell]

Joseph: thanks for the suggestion on the triangle method. I can certainly see how that would work, though coding the break-up of a polygon into triangles might take a bit of thought.

Paul: thanks for the graphical and line intersection suggestions. I've tried out your code for the first one - very nice and straightforward. The line intersection method is the one I have been thinking about most. Your line-crossing calculation code should make this simple and fast to implement.

In summary, thanks to suggestions on this forum, I now have four possible methods of determining whether a point is inside a polygonal region:

1. Graphical method, based on a graphic screen or bit map, where location is determined by pixel colour. This is a smart way of implementing a fine-scale grid approximation.

2. Line intersection/ray-casting method, where a line cast from the point will cross a polygon boundary an odd number of times if the point is inside and an even number of times if it is outside the polygon.

3. Triangle method, based on breaking the polygon up into triangles and determining whether the point is inside one of these triangles.

4. Winding number algorithm, based on summation of the angles subtended by each side of the polygon.

My method of choice is likely to be No. 2, but I could see myself using No. 1 as a fast and easy method to get started with. No. 4 is the least favoured candidate. Once I've found a bit of time to work on this (next week) I'll post some code.

Thanks for all the help.

[Charles Dietz]

Here is a routine that I have used in my section properties program. I have also included a simple example of its use.

#COMPILE EXE
#DIM ALL

TYPE polygonPts
   x AS SINGLE
   y AS SINGLE
END TYPE

FUNCTION PBMAIN
   LOCAL s AS STRING
   LOCAL polygons() AS polygonPts
   DIM polygons(1 TO 6)
   polygons(1).x = 0: polygons(1).y = 0
   polygons(2).x = 2: polygons(1).y = 0
   polygons(3).x = 2: polygons(1).y = 2
   polygons(4).x = 4: polygons(1).y = 2
   polygons(5).x = 4: polygons(1).y = 4
   polygons(6).x = 0: polygons(1).y = 4
   IF pointInsidePolygon(polygons(), 1, 1) THEN s = "inside" ELSE s = "outside"
   MSGBOX "This point is " + s
   IF pointInsidePolygon(polygons(), 3, 1) THEN s = "inside" ELSE s = "outside"
   MSGBOX "This point is " + s
END FUNCTION

FUNCTION pointInsidePolygon(polygons() AS polygonPts, x AS DOUBLE, y AS DOUBLE) AS LONG
   'http://astronomy.swin.edu.au/~pbourke/geometry/insidepoly/
   LOCAL PI, twoPI AS DOUBLE
   LOCAL i, j, first, last AS LONG
   LOCAL x1, y1, x2, y2, b, a1, a2, a, aSum AS DOUBLE
   PI = 3.141592653589793##: twoPI = 2*PI
   first = LBOUND(polygons): last = UBOUND(polygons)
   FOR i = first TO last
      j = i + 1: IF j > last THEN j = first
      x1 = polygons(i).x - x: y1 = polygons(i).y - y
      x2 = polygons(j).x - x: y2 = polygons(j).y - y 
      a1 = ATN(y1/x1): IF x1 < 0 THEN a1 = PI + a1
      a2 = ATN(y2/x2): IF x2 < 0 THEN a2 = PI + a2
      a = a2 - a1
      WHILE (a > PI): a = a - twoPI: WEND
      WHILE (a < -PI): a = a + twoPI: WEND
      aSum = aSum + a
   NEXT i
   IF ABS(aSum) > PI THEN FUNCTION = 1
END FUNCTION

[Arie Verheul]

Charles, as i am interested in something similar i tested your algorithm,
and it works great, except when the testpoint coincides with a vertex.
The algorithm hangs in that case. I added two lines to intercept this.
Coincidence of testpoint and vertex is then reported as outside.

Arie Verheul

   For i = first To last

      j = i + 1: If j > last Then j = first

      x1 = polygons(i).x - x: y1 = polygons(i).y - y
      x2 = polygons(j).x - x: y2 = polygons(j).y - y

      If x1 = 0 And y1 = 0 Then Function = 0: Exit Function	' Added these two lines to intercept
      If x2 = 0 And y2 = 0 Then Function = 0: Exit Function	' coincidence of testpoint with a vertex

      a1 = Atn(y1/x1): If x1 < 0 Then a1 = PI + a1
      a2 = Atn(y2/x2): If x2 < 0 Then a2 = PI + a2

      a = a2 - a1

      While (a > PI) : a = a - twoPI: Wend
      While (a < -PI): a = a + twoPI: Wend

      aSum = aSum + a

   Next i

[Michael Bell]

Charles, Arie: tried the routine with Arie's modification - works just fine - thanks! It also works pretty fast. I'm also working on a ray-casting routine to do the same job, which should work even faster. I'll post that when I've finished debugging it.

[Charles Dietz]

Thanks, Arie... that is an important addition.

[Michael Bell]

Here's a PB/CC test program for two inside polygon algorithms. Random locations for 1 million fish are found within an arena of 100 x 100 units. The program counts the number of fish within a polygon inside that arena (or square, or diamond, or triangle - commented out).

The two algorithms are a ray-casting method (e.g. http://en.wikipedia.org/wiki/Point_in_polygon) and a winding number method as implemented by Charles Dietz and modified by Arie Verheul. Ray-casting performs a fair bit faster than winding number, but this depends on the polygon shape and number of edges. For a 10 edge polygon the advantage is about an order of magnitude, but this shrinks a lot if a range test is added to the winding number method (only perform the calculations if the point is within the min to max ranges of x- and y-values for the polygon vertices).

No doubt it would be possible to tidy up and optimize the ray-casting code... Also faster if I didn't insist on using EXT precision...

#COMPILE EXE
#DIM ALL
#OPTIMIZE SPEED
 
TYPE Cartesian
  x AS EXT
  y AS EXT
END TYPE
 
TYPE Area
  COUNT AS LONG
  minval AS Cartesian
  maxval AS Cartesian
  xy(1 TO 100) AS Cartesian
END TYPE
 
'-----------------------------------------------------------------------------------------
FUNCTION PBMAIN () AS LONG
 
  '---------------------------------------------------------------------------------------
  ' TESTING ROUTINES FOR DETERMINING WHETHER A FISH LOCATION IS INSIDE OR OUTSIDE A
  ' POLYGONAL CLOSED AREA
  '---------------------------------------------------------------------------------------
 
  LOCAL ClosedArea AS Area
  LOCAL fishloc AS Cartesian
  LOCAL fish, Ninside1, Ninside2 AS LONG
  LOCAL Cycle1, Cycles1, Cycle2, Cycles2 AS QUAD
 
  RANDOMIZE TIMER
 
  ' polygon defining closed area
  ClosedArea.Count=10
  ClosedArea.minval.x=20 : ClosedArea.minval.y=10
  ClosedArea.maxval.x=83 : ClosedArea.maxval.y=85
  ClosedArea.xy(1).x=20 :  ClosedArea.xy(1).y=20
  ClosedArea.xy(2).x=40 :  ClosedArea.xy(2).y=50
  ClosedArea.xy(3).x=30 :  ClosedArea.xy(3).y=55
  ClosedArea.xy(4).x=50 :  ClosedArea.xy(4).y=85
  ClosedArea.xy(5).x=52 :  ClosedArea.xy(5).y=60
  ClosedArea.xy(6).x=80 :  ClosedArea.xy(6).y=60
  ClosedArea.xy(7).x=60 :  ClosedArea.xy(7).y=50
  ClosedArea.xy(8).x=60 :  ClosedArea.xy(8).y=40
  ClosedArea.xy(9).x=83 :  ClosedArea.xy(9).y=38
  ClosedArea.xy(10).x=62 : ClosedArea.xy(10).y=10
 
  ' square defining closed area
'  ClosedArea.Count=4
'  ClosedArea.minval.x=25 : ClosedArea.minval.y=25
'  ClosedArea.maxval.x=75 : ClosedArea.maxval.y=75
'  ClosedArea.xy(1).x=25 :  ClosedArea.xy(1).y=25
'  ClosedArea.xy(2).x=25 :  ClosedArea.xy(2).y=75
'  ClosedArea.xy(3).x=75 :  ClosedArea.xy(3).y=75
'  ClosedArea.xy(4).x=75 :  ClosedArea.xy(4).y=25
 
  ' diamond defining closed area
'  ClosedArea.Count=4
'  ClosedArea.minval.x=0 : ClosedArea.minval.y=0
'  ClosedArea.maxval.x=100 : ClosedArea.maxval.y=100
'  ClosedArea.xy(1).x=0 :  ClosedArea.xy(1).y=50
'  ClosedArea.xy(2).x=50 :  ClosedArea.xy(2).y=100
'  ClosedArea.xy(3).x=100 :  ClosedArea.xy(3).y=50
'  ClosedArea.xy(4).x=50 :  ClosedArea.xy(4).y=0
 
  ' triangle defining closed area
'  ClosedArea.Count=3
' ClosedArea.minval.x=25 : ClosedArea.minval.y=25
' ClosedArea.maxval.x=75 : ClosedArea.maxval.y=75
' ClosedArea.xy(1).x=25 :  ClosedArea.xy(1).y=25
' ClosedArea.xy(2).x=50 :  ClosedArea.xy(2).y=75
' ClosedArea.xy(3).x=75 :  ClosedArea.xy(3).y=25
 
  ' 1000000 fish with uniform random distribution within 100 x 100 arena
  Ninside1=0 : Ninside2=0
  FOR fish=1 TO 1000000
 
    fishloc.x=RND*100
    fishloc.y=RND*100
 
    TIX Cycle1
    IF ISTRUE(InsidePoly(fishloc,ClosedArea)) THEN INCR Ninside1
    TIX END Cycle1
    Cycles1=Cycles1+Cycle1
 
    TIX Cycle2
    IF ISTRUE(pointInsidePolygon(fishloc,ClosedArea)) THEN INCR Ninside2
    TIX END Cycle2
    Cycles2=Cycles2+Cycle2
 
    IF fish MOD 10000 = 0 THEN
      LOCATE 2,3
      PRINT USING$("########## fish",fish)
      PRINT
      PRINT USING$("  ########## inside closed area by ray-casting algorithm",Ninside1)
      PRINT USING$("  ########## cycles",Cycles1)
      PRINT
      PRINT USING$("  ########## inside closed area by winding number algorithm",Ninside2)
      PRINT USING$("  ########## cycles",Cycles2)
    END IF
 
  NEXT fish
 
  PRINT
  PRINT "  Press any key to exit...";
  WAITKEY$
 
END FUNCTION
 
'-----------------------------------------------------------------------------------------
 
FUNCTION InsidePoly(BYVAL fish AS Cartesian, BYREF Closure AS Area) AS LONG
 
  '---------------------------------------------------------------------------------------
  ' Determines whether the location of the fish is within the polygonal region
  ' defined by Closure.  Uses ray-casting along a line of equal x-values.
  '---------------------------------------------------------------------------------------
 
  LOCAL i, j, ncrosses AS LONG
  LOCAL dx, dy, currenty, eps, a, b AS EXT
 
  ' first determine whether the fish location is within the range of the polygon
  IF fish.x<Closure.minval.x OR fish.x>Closure.maxval.x OR _
     fish.y<Closure.minval.y OR fish.y>Closure.maxval.y THEN
    FUNCTION=0 ' FALSE - point is outside the polygon's range, no further tests needed
  ELSE
 
    ' fish is within range, so cast a line in a westerly direction
    ncrosses=0&
    eps=1e-16##
    currenty=fish.y
 
    ' look at each edge, one at a time
    FOR i=1 TO Closure.count
 
      IF i=Closure.count THEN j=1 ELSE j=i+1
 
      ' test for coincidence with vertices, jigger away if necessary
      IF currenty=Closure.xy(i).y OR currenty=Closure.xy(j).y THEN currenty=currenty+eps
 
      ' for crossing, point must lie within y-range of the edge
      IF currenty<MAX(Closure.xy(i).y,Closure.xy(j).y) AND _
         currenty>MIN(Closure.xy(i).y,Closure.xy(j).y) THEN
 
        ' we can ignore horizontal (E-W) edges (dy=0)
        dy=Closure.xy(j).y-Closure.xy(i).y
        IF dy<>0## THEN
 
          dx=Closure.xy(j).x-Closure.xy(i).x
          IF dx=0## THEN
            ' vertical (N-S) line - crossing occurs if point is to the east
            IF fish.x<=Closure.xy(i).x THEN INCR ncrosses
          ELSE
            ' line is not vertical, so calculate equation y=a+bx
            b=dy/dx
            a=Closure.xy(i).y-b*Closure.xy(i).x
            ' find crossing point of edge and cast ray
            ' check that crossing point is to the west of the test point
            IF ((currenty-a)/b)>=fish.x THEN INCR ncrosses
          END IF ' dx test
 
        END IF ' horizontal line test
      END IF ' y-range test
 
    NEXT i
 
    'test for odd or even
    IF ncrosses MOD 2& > 0& THEN
      FUNCTION=-1 ' TRUE - odd number of line crosses
    ELSE
      FUNCTION=0  ' FALSE - even number of line crosses
    END IF
 
  END IF
 
END FUNCTION
 
'-----------------------------------------------------------------------------------------
 
FUNCTION pointInsidePolygon(BYVAL fish AS Cartesian, BYREF Closure AS Area) AS LONG
 
   '--------------------------------------------------------------------------------------
   ' Charles Dietz's algorithm for testing whether a point is within a polygon,
   ' using Arie Verheul's modification to deal with points falling on a vertex.
   ' Added range test for points outside max and min x- and y-values of polygon.
   ' [URL]http://astronomy.swin.edu.au/~pbourke/geometry/insidepoly/[/URL]
   '--------------------------------------------------------------------------------------
 
   LOCAL PI, twoPI AS EXT
   LOCAL i, j, first, last AS LONG
   LOCAL x1, y1, x2, y2, b, a1, a2, a, aSum AS EXT
 
   ' range test
   IF fish.x<Closure.minval.x OR fish.x>Closure.maxval.x OR _
      fish.y<Closure.minval.y OR fish.y>Closure.maxval.y THEN
     FUNCTION=0 ' FALSE - point is outside the polygon's range, no further tests needed
     EXIT FUNCTION
   END IF
 
   PI = 3.141592653589793##: twoPI = 2*PI
   first = 1: last = Closure.count
   FOR i = first TO last
      j = i + 1: IF j > last THEN j = first
      x1 = Closure.xy(i).x - fish.x: y1 = Closure.xy(i).y - fish.y
      x2 = Closure.xy(j).x - fish.x: y2 = Closure.xy(j).y - fish.y
 
      ' Arie Verheul's modification to test coincidence of test point with a vertex
      IF x1 = 0 AND y1 = 0 THEN FUNCTION = 0: EXIT FUNCTION
      IF x2 = 0 AND y2 = 0 THEN FUNCTION = 0: EXIT FUNCTION
 
      a1 = ATN(y1/x1): IF x1 < 0 THEN a1 = PI + a1
      a2 = ATN(y2/x2): IF x2 < 0 THEN a2 = PI + a2
      a = a2 - a1
      WHILE (a > PI): a = a - twoPI: WEND
      WHILE (a < -PI): a = a + twoPI: WEND
      aSum = aSum + a
   NEXT i
 
   IF ABS(aSum) > PI THEN
     FUNCTION = -1 ' TRUE
   ELSE
     FUNCTION = 0  ' FALSE
   END IF
END FUNCTION