Inside or outside a polygonal region?

[Michael Bell]

Just tried Arie's latest version of the PointInPolygon algorithm. It works really well - about 25-30% faster than my previous ray-casting version (depending on whether or not you test for polygon range) and about 85% faster than the first winding number algorithm given in this thread. Good work! That settles my choice of algorithm for future use. Thanks everyone for input on this thread.

[Arie Verheul]

Point_In_Polgygon algorithm - final update

I did some final tidying up, and for the convenience of possible future users
of this thread it might be useful to post a well commented version of the code.

As the winding algorithm is actually a disguised ray casting algorithm,
it is more clear to speak here about a ray casting algorithm.
Therefore i changed some names and comments accordingly and removed
references to the winding algorithm from which this algortithm originated.

I also removed the GoSub calls as this gives some speed gain.
Functionally the code below is the same that of the Winding algorithm-3rd revision.
Note that the function is direction sensitive, and a point inside the polygon may
therefore return +/- 1, depending on the way the polygon was defined.

For the moment i don't see any possibilities for further improvement, but if someone does ..

Thanks to all contributors to this thread.

---------------------------------------------------------------------------

Arie Verheul

In the last comparison i erroneously used : If y1 * (x1 - x2) / (y1 - y2) > x1 Then ..
This should of course be : If y1 * (x1 - x2) / (y1 - y2) < x1 Then ..
I corrected this. My apologies.

Arie Verheul

Type Cartesian
   x As Ext
   y As Ext
End Type

Dim Vertex(Points) as Cartesian

'
'
'

Function PointInPolygon(Vertex() As Cartesian, x As Ext, y As Ext) As Long

    Local   x1,y1,x2,y2             As Ext
    Local   i, IntersectionBalance  As Long
    Local   Coord                   As Ext Pointer

    Coord = VarPtr(Vertex(0).x)

    '-------------------------

    ' Function to test if a point (x,y) is inside or outside a polygonal area
    ' defined by the array Vertex().	

    ' Return values :
    '  0  if point is outside polygon
    '  1  if point is inside polygon, and polygon was defined clockwise
    ' -1  if point is inside polygon, and polygon was defined counter clockwise

    ' Polygon definition array must start at Vertex(0)
    ' Vertex(0) must be repeated as last point in Vertex()

    '------------------------- Read first point

    x1 = @Coord - x : Incr Coord
    y1 = @Coord - y : Incr Coord

    If x1 = 0 And y1 = 0 Then Function = 0 : Exit Function

    '------------------------- Investigate polygon contour for intersection with positive x-axis

    For i =  1 To UBound(Vertex)

        x2 = x1
        y2 = y1

        x1 = @Coord - x : Incr Coord                            ' Maps vertices to a local coord system
        y1 = @Coord - y : Incr Coord                            ' with test point as origin

        If x1 = 0 And y1 = 0 Then Function = 0 : Exit Function  ' Function fails if test point coincides
                                                                ' with a vertex

        If       y1 >= 0 And y2 < 0 Or _                        ' Test if both vertices are
                 y2 >= 0 And y1 < 0 Then                        ' on opposite sides of the x-axis

           If y1 * (x1 - x2) / (y1 - y2) < x1 Then              ' If so, test if connecting edge
                                                                ' intersects with +x half axis

               IntersectionBalance = IntersectionBalance + Sgn(y1 - y2)
           End If
        End If

    Next

    Function = IntersectionBalance

End Function

[Charles Dietz]

Arie, your last algorithm is very nice, in fact, its the one I'm going to use in the future. Like you said, I think it is really a ray casting algorithm, and since it is twice as fast as our best efforts with the winding method, I think we have to conclude that the ray casting algorithm wins. And I really do favor the way you redefined the coordinate system, allowing you to ignore many of the polygon sides.

Following this thread and reviewing the ideas of its contributors has been very rewarding. Thank you all.

[Arie Verheul]

Winding algorithm - third revision

In an earlier contribution it was argued that there is no need to
actually add up all individual angles, but that it is sufficient
to keep account of the 2*Pi adjustments applied.
Such adjustments are necessary if an individual angle encompasses
the arbitrarily placed gap in the angular scale.

Taking this concept one step further, it may be realised that there is
no need to actually calculate the individual angles either, but that
it is sufficicient to test if an individual angle encompasses the gap
in the angular scale. There are more efficient ways to do this.
This means that the time consuming ATN calculations may be eliminated
from the algorithm.

Assuming that the angular scale runs counterclockwise from the positive
x-axis, this test involves two simple steps :

Test if both vertices (x1,y1) and (x2,y2) that define the
individual angle are on either side of the x-axis,

And if so, test if the edge connecting them intersects the positive x-axis.
This can be done with a simple expression.

Further, as Charles proposed, an important bit of overhead may be avoided by
including the first vertex of the polygon definition at the end of the list again,
and by agreeing that the polygon definition array starts at 0.

The result is an elegant and pretty fast algorithm.
It is clear that this algorithm is almost a ray casting algorithm, however with
one difference: a ray casting algorithm counts the number of crossings,
while this algorithm keeps account of the sense of the crossings.
Because of this it can more easily handle special cases where a vertex,
or even an edge is on the x-axis.

Arie Verheul

Type Cartesian
   x As Ext
   y As Ext
End Type

Dim Vertex(Points) as Cartesian


Function PointInsidePolygon_Economy(Vertex() As Cartesian, x As Ext, y As Ext) As Long

    Local   x1,y1,x2,y2        As Ext
    Local   i, AdjustmentCount As Long
    Local   Coord              As Ext Pointer

    Coord = VarPtr(Vertex(0).x)
    
    '------------------------- Read first point

    GoSub MapPoint

    '------------------------- Investigate angles

    For i =  1 To UBound(Vertex)

        GoSub MapPoint
        GoSub CheckAngle
    Next

    '--------------------------

    Function = AdjustmentCount

Exit Function

'--------------------------

MapPoint:

    ' Maps vertices onto local coord system with test point as origin
    ' (x1,y1) is leading vertex, (x2,y2) follows

    x2 = x1
    y2 = y1

    x1 = @Coord - x : Incr Coord
    y1 = @Coord - y : Incr Coord

    If x1 = 0 And y1 = 0 Then Function = 0 : Exit Function

Return

CheckAngle:

    ' Tests if individual angle encompasses positive x-axis

    If y1 >= 0 And y2 < 0 Or _
       y2 >= 0 And y1 < 0 Then
       
       If y1 * (x2 - x1) / (y2 -y1) > x1 Then

           AdjustmentCount = AdjustmentCount + Sgn(y2 - y1)
       
       End If
    End If
Return

End Function

[Arie Verheul]

Thanks Charles, now it is all clear to me.

Arie Verheul

[Charles Dietz]

Arie, I should have been more clear... the polygon for this last routine must be redefined with its dimension increased by one to allow for the first node to be repeated as the ending node. Using my earlier case, the polygon would be:

DIM polygons(1 TO 7) AS polygonPts 'increased by one

polygons(1).x = 0: polygons(1).y = 0 'first node
polygons(2).x = 2: polygons(2).y = 0
polygons(3).x = 2: polygons(3).y = 2
polygons(4).x = 4: polygons(4).y = 2
polygons(5).x = 4: polygons(5).y = 4
polygons(6).x = 0: polygons(6).y = 4
polygons(7).x = 0: polygons(7).y = 0 'repeating 1st node

The old angle for node 6 to node 1 is replaced now by the angle for node 6 to node 7. As you can see, my new code is just the old code streamlined, avoiding the double calculations. I'm still getting the aSum answer from 0 to 2* PI.

[Michael Bell]

Thanks John. Yes, that does clarify things a lot. I wasn't trying to set the S and T values from the calling routine, but I did wonder whether I needed to do something with the return values. Knowing what they are means that, as you suggest, I can remove them from the argument list.

[John R. Heathcote]

Michael,

Oops, one more thing, the "S" and "T" parameters are calculated by the intersect function and the values are returned to the calling routine (sometimes I have further need to access this information). You do not need to set them from your calling routine. In fact, if you don't have a need for these two values then you can remove them as parameters from the argument list.

[John R. Heathcote]

Michael,

The parameters "S" and "T" are percentage indicators where the intersection point occurs on both line segments. "S" applies to the line segment from point 1 to 2, "T" are points 3 to 4. If "S" and "T" are both > 0 and < 1 the intersection occurs between the end points for both line segments.

If "S" = 0 and "T" > 0 and < 1 then the begin point of the first line segment falls on the second line segment. If "S" = 1 then the end point of the first line segment falls on the second line segment. Reverse this logic for the second line segment. This describes a possible ambiguity or degenerate condition with your line segment geometry, however, in your case I would assume you would want to check for the ambiguity condition all the time.

The lType parameter determines the type of line segment crossings you wish to check: 0 = actual line segment crossings, 1 = actual line segment crossings and end points falling on another line segment. I use this function for a variety of crossing conditions in my apps hence the general nature of the source code.

There are actually 3 return values: 0 = No intersection, 1 = both line segments parallel to each other (or on each other), and -1 = an intersection occurred. So the function doesn't return a specific %TRUE/%FALSE value.

Hope this clarifies things a bit.

[Michael Bell]

Arie, Charles: big improvements on the winding method. In my test (one million test points, ten edge polygon), for every 100 CPU cycles taken by the unmodified winding method, Arie's modification takes 57 and Charles's modification takes 56. This compares with 22 using the ray-casting method (but see reservation below). All algorithms are implemented with range testing for the test point being within the max and min x- and y-ranges for the polygon.

John: I take your point about near vertical and horizontal lines with extended precision variables. I've implemented your routine, but am not sure whether I'm using it right - can you tell me what lType and the S and T variables mean? Setting lType=1, and using your function to give me true or false I don't get the same answers as the other routines, but this is probably because I am not understanding (or using) the S and T values. My call to your function is

 
IF ISTRUE(isect2d_DoLinesCross_XY(fish.x,fish.y, _
          Closure.maxval.x+1,fish.y, _
          Closure.xy(i).x,Closure.xy(i).y, _
          Closure.xy(j).x,Closure.xy(j).y, _
          1,S,T)) THEN INCR ncrosses

I might also try Paul Dixon's line crossing routine http://www.powerbasic.com/support/pb...ad.php as suggested in post #14

Thanks everybody! I'm consistently impressed by the time and thought people give to questions posed on these forums.

[Arie Verheul]

Charles, your rewritten function is a miracle.
It is considerably faster than the original, and even faster than my approach.
Partly this is because you skipped the angle between vertices 1 and Last.
Surprisingly enough this does not seem to matter.
I tested the function, and so far it always produced the correct result.

However, i do not really understand the math behind this.
In the original algorithm the sum of all individual angles is either 0 or 2*Pi,
but in your function it can have all sorts of values.
I usually worry a bit about things that i do not understand, so if you have
any explicit thoughts about this, could you please elaborate a bit on it.

Arie Verheul

[Charles Dietz]

Following Arie's lead, and repeating the first polygon node with dimensions increased by one, I was able to save about 43% with this routine.

FUNCTION pointInsidePolygon(polygons() AS polygonPts, x AS DOUBLE, y AS DOUBLE) AS LONG
   'http://astronomy.swin.edu.au/~pbourke/geometry/insidepoly/
   LOCAL i, last AS LONG
   LOCAL xx, yy, aa, a0, a, aSum, PI, twoPI AS DOUBLE
   PI = 3.141592653589793##: twoPI = 2*PI
   last = UBOUND(polygons)
   xx = polygons(1).x - x: yy = polygons(1).y - y
   IF xx = 0 AND yy = 0 THEN EXIT FUNCTION
   a0 = ATN(yy/xx): IF xx < 0 THEN a0 = PI + a0
   FOR i = 2 TO last
      xx = polygons(i).x - x: yy = polygons(i).y - y
      IF xx = 0 AND yy = 0 THEN EXIT FUNCTION
      aa = ATN(yy/xx): IF xx < 0 THEN aa = PI + aa
      a = aa - a0
      WHILE (a > PI): a = a - twoPI: WEND
      WHILE (a < -PI): a = a + twoPI: WEND
      aSum = aSum + a: a0 = aa
   NEXT i
   IF ABS(aSum) > PI THEN FUNCTION = 1
END FUNCTION

[John R. Heathcote]

Michael,

In your original example you determine if a polygon edge crosses the ray casting segment by computing the slope of the closure edge. You may find that the closer such an edge approaches vertical or horizontal unexpected and unwanted side effects can occur in your slope calculation, especially using extended precision variables.

For computing the possibility of two line segments intersecting I use:

FUNCTION isect2d_DoLinesCross_XY ALIAS "isect2d_DoLinesCross_XY" _
                             (BYVAL X1 AS DOUBLE, _
                              BYVAL Y1 AS DOUBLE, _
                              BYVAL X2 AS DOUBLE, _
                              BYVAL Y2 AS DOUBLE, _
                              BYVAL X3 AS DOUBLE, _
                              BYVAL Y3 AS DOUBLE, _
                              BYVAL X4 AS DOUBLE, _
                              BYVAL Y4 AS DOUBLE, _
                              BYVAL lType AS LONG, _
                              S AS DOUBLE, _
                              T AS DOUBLE) EXPORT AS LONG

  LOCAL dNum AS DOUBLE
  LOCAL dDenom AS DOUBLE
  '-----------------------------------------------------------------------

  FUNCTION = 0      'Set default status, Line segments do NOT cross.

  dDenom = (((X2 - X1) * (Y4 - Y3)) - ((Y2 - Y1) * (X4 - X3)))

    IF dDenom = 0 THEN                'Line Segments are parallel?
      FUNCTION = 1

      EXIT FUNCTION
    END IF

  S = (((Y1 - Y3) * (X4 - X3)) - ((X1 - X3) * (Y4 - Y3))) / dDenom
  T = (((Y1 - Y3) * (X2 - X1)) - ((X1 - X3) * (Y2 - Y1))) / dDenom

  SELECT CASE lType    'Check segment crossing type?
    CASE 0              'Check if the line segments actually do cross (default).
        IF S > 0 AND S < 1 AND T > 0 AND T < 1 THEN
          FUNCTION = -1
        END IF

    CASE 1              'Check degenerate (possible ambiguity) condition.
        IF S >= 0 AND S <= 1 AND T >= 0 AND T <= 1 THEN
          FUNCTION = -1
        END IF

  END SELECT

END FUNCTION

This function uses doubles as paramters, but should be easily expanded to extended precision.

Since you are only interested in the closure segments that do cross your ray casting line rather than computing the actual intersection point, this FUNCTION should do nicely. The advantage of this function is the fact that it can determine the possibility of intersection regardless of the original orientation of the two line segments.

The parameters "S" and "T" will tell you what type of intersection has occured. In your case you will probably want to set lType to 1 to return the possibility of one closure end point falling on your ray cast line.

Just as a suggestion:
You can also set "ncrosses" as a single, then when you encounter the possibility of an end point falling on your ray cast line increment "ncrosses" by 0.5. The other closure line segment that shares this end point will also fall on the ray cast line and increment "ncrosses" by another 0.5. Otherwise increment "ncrosses" by 1 when a crossing occurs.

Set your ray cast line to use a cardinal axis direction for your test, for example the X axis. You can then trivially reject the closure segments where both ends of the closure line segment Y coordinates are either above or below your ray cast line because they will never cross.

Hope this helps.

[Arie Verheul]

Winding number algorithm - Economised version

I realised that in Charles' function much of the work is done twice,
and wondered how much could be earned by rewriting it in such a way
that this double work is avoided. Well, not very much.

My first attempt was to calculate the azimuths of all vertices first
and write them to an aray to be analysed later. Completely wrong.
The array slowed the proces down by an order of magnitude.

So i tried it the other way round, avoided the array, and dimensioned
certain variables as Register. This worked.

Second i tried the effect of floating point variable type. No effect.
PB does all floating point math with extended precision, and dimensoning
smaller variable types did not speed things up.

Third, as arrays had shown to slow things down, i accessed Polygons()
with a pointer. This resulted in a ~ 10% speed gain.

Fourth i rewrote the function so that double work was avoided.
This had a clear but limited effect.

Fifth i avoided conditional jumps in the code where possible, and wrote things out.
This gave some speed gain at the expense of a few lines of extra code.

Sixth i realised that there is no need to actually add up all individual angles as

[Phi(2) - Phi(1)] + [Phi(3) - Phi(2)] + ... + [ Phi(n) - Phi(n-1)] + [Phi(1) - Phi(n)] = 0

The result may differ from zero if one or more angles include the gap in the angular scale,
and an adjustment by 2*Pi has to be made. There is no need to actually calculate this,
and it is sufficient to keep count of the number and sign of 2*Pi adjustments.

So i actually tried to eliminate any inefficiency i could find. With limited result.
With floating point variables dimensioned as extended, all attempts to increase
the efficiency produced a speed gain of at most a factor of two, and often less.
I found in a separate test that about half of the time is used to do the ATN
calculations, so this sets a clear limit to the speed of the algorithm.

Arie Verheul

Type PolygonPts
   x As Extended
   y As Extended
End Type

'----------------------------------------------------------------------------------------

Function PointInsidePolygon (Polygons() As PolygonPts, x As Ext, y As Ext) As Long

    Register    i As Long
    Register    Pi As Ext, dx As Ext, dy As Ext, Subtended_Angle As Ext
    Local       Azimuth_First, Azimuth_1, Azimuth_2 As Ext
    Local       First, Last, AdjustmentCount As Long
    Local       Coord As Extended Pointer

    Pi    = 4*Atn(1##)

    First = LBound(Polygons)
    Last  = UBound(Polygons)

    Coord = VarPtr(Polygons(First).x)

    '-------------------------- Investigate first edge

    GoSub Calculate_Azimuth

    Azimuth_First = Azimuth_2       ' Store azimuth of first vertex to be reused later

    GoSub Calculate_Azimuth

    Subtended_Angle = Azimuth_2 - Azimuth_First

    ' Keep count of number and sign of 2*Pi adjustments

    If Abs(Subtended_Angle) > Pi Then AdjustmentCount = AdjustmentCount - Sgn(Subtended_Angle)

    Azimuth_1 = Azimuth_2           ' Use value to investigate next edge


    '-------------------------- Investigate further edges

    For i = First + 2 To Last

        GoSub Calculate_Azimuth

        Subtended_Angle = Azimuth_2 - Azimuth_1

        If Abs(Subtended_Angle) > Pi Then AdjustmentCount = AdjustmentCount - Sgn(Subtended_Angle)

        Azimuth_1 = Azimuth_2
   Next


    '-------------------------- Investigate last edge

    Subtended_Angle = Azimuth_First - Azimuth_1

    If Abs(Subtended_Angle) > Pi Then AdjustmentCount = AdjustmentCount - Sgn(Subtended_Angle)

    Function = AdjustmentCount

Exit Function

'--------------------------

Calculate_Azimuth:

    dx = @Coord - x

    Incr Coord

    dy = @Coord - y

    Incr Coord

    If dx = 0 And dy = 0 Then Function = 0: Exit Function

    ' In second and third quadrant add Pi, so Azimuth ranges from -.5*Pi to 1.5*Pi

    Azimuth_2 = IIf(dx < 0, Atn(dy/dx) + Pi, Atn(dy/dx))

Return

End Function

[Michael Bell]

Further to my last post...

I also tried using a memory bitmap, reading pixel colours, but ran into size problems. I tried a bitmap of 10000 by 10000, which for my purposes is the minimum precision I would be happy with (i.e. two decimal places for locations on a 100 x 100 grid) - Error 5. A little search of the forums turned up the reason why: http://www.powerbasic.com/support/pb...ad.php?t=39542

1000 by 1000 worked OK, but this is a very coarse grid approximation, too coarse for my purposes (I am working with movements at quite a small spatial scale for each time step of the simulation model). Reading pixel colour also seems to take longer than running the other two routines, though I don't see why.

Conclusion: memory bitmaps aren't the answer to my purpose, but could be OK for applications that did not require a high degree of spatial precision and in which speed of calculation was not an issue.

[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

[Charles Dietz]

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

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

[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

[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