Announcement

Collapse
No announcement yet.

Inside or outside a polygonal region?

Collapse
X
 
  • Filter
  • Time
  • Show
Clear All
new posts

  • Inside or outside a polygonal region?

    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

  • #2
    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,
    Last edited by Manuel Valdes; 3 Sep 2009, 08:53 AM.

    Comment


    • #3
      CreatePolygonRgn () ==> Polygon
      PtInRegion() true/false
      Michael Mattias
      Tal Systems Inc. (retired)
      Racine WI USA
      [email protected]
      http://www.talsystems.com

      Comment


      • #4
        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:

        Code:
         
        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.

        Comment


        • #5
          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,
          Last edited by Manuel Valdes; 3 Sep 2009, 10:18 AM.

          Comment


          • #6
            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.

            Comment


            • #7
              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

              Comment


              • #8
                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...

                Comment


                • #9
                  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
                  Regards,
                  Marc

                  Comment


                  • #10
                    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.

                    Comment


                    • #11
                      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.
                      Last edited by Joseph Cote; 3 Sep 2009, 03:17 PM. Reason: spelling; clarification
                      The boy just ain't right.

                      Comment


                      • #12
                        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.

                        Code:
                        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
                        The boy just ain't right.

                        Comment


                        • #13
                          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:
                          Code:
                          '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.
                          Last edited by Paul Dixon; 4 Sep 2009, 07:19 AM.

                          Comment


                          • #14
                            Mike,
                            I don't have time to code an example now but the alternative which gives accurate results goes as follows.
                            Code:
                            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/pb...ad.php?t=23695

                            Paul.

                            Comment


                            • #15
                              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.

                              Comment


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

                                Code:
                                #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

                                Comment


                                • #17
                                  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

                                  Code:
                                     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

                                  Comment


                                  • #18
                                    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.

                                    Comment


                                    • #19
                                      Thanks, Arie... that is an important addition.

                                      Comment


                                      • #20
                                        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...

                                        Code:
                                        #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.
                                           ' http://astronomy.swin.edu.au/~pbourk...ry/insidepoly/
                                           '--------------------------------------------------------------------------------------
                                         
                                           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

                                        Comment

                                        Working...
                                        X