Announcement
Collapse
No announcement yet.
Inside or outside a polygonal region?
Collapse
X

Just tried Arie's latest version of the PointInPolygon algorithm. It works really well  about 2530% faster than my previous raycasting 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.


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 algorithm3rd 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
Code: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 xaxis 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 xaxis 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
Leave a comment:


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.
Leave a comment:


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
xaxis, 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 xaxis,
And if so, test if the edge connecting them intersects the positive xaxis.
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 xaxis.
Arie Verheul
Code: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 xaxis 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
Leave a comment:


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.
Leave a comment:


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.
Leave a comment:


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.
Leave a comment:


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.
Leave a comment:


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 raycasting method (but see reservation below). All algorithms are implemented with range testing for the test point being within the max and min x and yranges 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
Code: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
Thanks everybody! I'm consistently impressed by the time and thought people give to questions posed on these forums.
Leave a comment:


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
Leave a comment:


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.
Code: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
Leave a comment:


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:
Code: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
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.
Leave a comment:


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(n1)] + [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
Code: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
Leave a comment:


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.
Leave a comment:


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 raycasting 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. Raycasting 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 yvalues for the polygon vertices).
No doubt it would be possible to tidy up and optimize the raycasting 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 raycasting 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 raycasting along a line of equal xvalues. ' 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=1e16## 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 yrange 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 (EW) edges (dy=0) dy=Closure.xy(j).yClosure.xy(i).y IF dy<>0## THEN dx=Closure.xy(j).xClosure.xy(i).x IF dx=0## THEN ' vertical (NS) 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).yb*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 ((currentya)/b)>=fish.x THEN INCR ncrosses END IF ' dx test END IF ' horizontal line test END IF ' yrange 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 yvalues 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
Leave a comment:


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


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
Leave a comment:


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
Leave a comment:

Leave a comment: