Point in Polygon-Algorithm

Hi folks,
I'm looking for an algorithm determining if a point is within a given
polygon or not. I found some for plane geometry, but I need one in
shperical geometry - i.e. working with GPS points.
Any links or pointers appreciated!

tia
Roland

What are the longest distances between your points?

Uffe Kousgaard wrote:


I know there is not much failure if I do the calculation planar. But
then it's not "clean".

Regards,
Roland

Hi Roland,

Would you mind posting what you come up with? (If you'd rather, you
could just email me. :-)

If you were content to use a planar algorithm, you could use the
algorithm for ccw(p0,p1,p2) from Robert Sedgewick's book _Algorithms in
C_. This algorithm returns +1 if point p2 is to the left of or in front
of the line segment from p0 to p1, 0 if it is on the line segment, and
-1 if it is to the right of or behind the line segment. If the polygon
is described in a counterclockwise direction, then p2 is inside it if p2
returns a +1 when p0-p1 describes each of the line segments in the
polygon. If the polygon is described in a clockwise direction, you need
-1's. If the polygon is a collection of unordered sides, either find
another algorithm or preprocess all your polygons.

I do not guarantee this is the fastest way to answer the question,
but the ccw function is very simple and fast. If the algorithms you have
found are better, I'd be interested in hearing about them.

You could probably revise the ccw function for spherical trig, but I
suspect it would lose its simplicity. On the other hand, when the
problem is stated this way, it becomes an issue of topology, so the
spherical version might be simple too. A good source for spherical trig

For what it's worth, here is the ccw(p0,p1,p2) function, as described
in my notes. I don't have a copy of Sedgewick's book an hand at the
moment, so this version might be a little different than his. The x and
y components of point pn (where n = 0, 1, 2) are pxn and pyn,
respectively:

ccw( p0, p1, p2 )

dx1 = px1 - px0
dy1 = py1 - py0
dx2 = px2 - px0
dy2 = py2 - py0

elseif dx1 * dy2 < dx1 * dx2 return with -1
elseif dx1 * dx2 < 0 or dy1 * dy2 < 0 return with -1
elseif dx1^2 + dy1^2 < dx2^2 + dy2^2 return with 0
endif

end

Disclaimer: Test the algorithm before you use it! I may have made a
typo, something may have gotten garbled in transmission, you may have
misunderstood my pseudocode, or a gremlin may have made _you_ make a
typo.


Bob

Bob (Robert G.) Chamberlain | I have yet to see any problem, however
rgc@jpl.nasa.gov | complicated, which, when looked at in
Opinions & quips are mine | the right way, did not become still more
- or public - not JPL's | complicated. - Poul Anderson

Ooooooops: I did make a typo in the pseudocode for the ccw function. The
first if statement (at least) is wrong. Since it is wrong in my notes
and I do not have a copy of Sedgewick available, please ignore the
pseudocode I gave. * blush *


Bob (Robert G.) Chamberlain | If you need any _more_ help,
rgc@jpl.nasa.gov | don't hesitate to call.
Opinions & quips are mine | My number is 911.
- or public - not JPL's | ;-)


In article