Point in Polygon Problem
Dick Fulton
I'm looking for a method/algorithm to determine if a given point lies
inside or outside of a predefined polygon. The point and the polygon
vertices are latitude/longitude points on the earth's surface. The
methods that come to mind seem difficult to implement as a computer
program algorithm, which is what I need to do.
Anyone been down this road or able to point me somewhere?
Thanks...
--
Dick Fulton
Palo Alto, CA
email: dfu...@dinc.com
Chris Monico
On Mon, 19 Jan 1998 22:40:30 -0800, dfu...@dinc.com (Dick Fulton)
wrote:
Good books on computer graphics usually deal with this problem as
Polygon clipping (i.e., clipping a point against a polygon). The
method usually boils down to computing a normal vector at some point
on each edge (bisecting point, for simplicity), computing vectors from
the given point to each of the edge points, and taking some dot
products (innner products) of the vectors. Then recall the formula
A*B = cos(angle between A & B)
---------------------------------------
|A| |B|
where * denotes the inner product.
Rainer Thonnes
In article <dfulton-1901...@dynamic27.pm02.mv.best.com>,
dfu...@dinc.com (Dick Fulton) writes:
> I'm looking for a method/algorithm to determine if a given point lies
> inside or outside of a predefined polygon. The point and the polygon
> vertices are latitude/longitude points on the earth's surface. The
> methods that come to mind seem difficult to implement as a computer
> program algorithm, which is what I need to do.
Go on, risk embarrassment by revealing the methods that came to mind.
I'll risk embarrassment by saying what the obvious method that comes
to *my* mind is. It exploits properties of vector products.
In particular, the 3-dimensional vector product AxB of two vectors A and B
in the xy plane, whose z components are therefore zero, is a vector whose
x and y components are zero. What matters is its z component, which for
A=(a,b), B=(c,d), is ad-bc. If the sign of this quantity is positive/zero/
negative then the rotation needed to get from A to B is anticlockwise/zero/
clockwise, i.e. in the zero case A and B are parallel.
This means that if the feather ends of the arrows represented by the two
vectors are tied together, and you're standing at their common point,
facing the tip of arrow A, then if ad-bc is +/0/- then the tip of arrow
B is (on your left)/(straight ahead or behind you or where you are)/
(on your right).
So all you need to do, given the coordinates of the point in question (call
this vector P), and those of the polygon's n vertices (call these vectors
Vi, numbered 0 to n-1, and for convenience define Vn=V0), is this:
Compute the n edge vectors Ei=V(i+1)-Vi, which join vertex i to i+1,
and the n radius vectors Ri=P-Vi, which join vertex i to the point. Then,
for each i, compute the z component of the vector product EixRi. If they
are all positive, then P is "to the left of" each of the edge vectors, and
therefore inside the polygon. You may like to allow the results to be zero
which would allow the point to be on the boundary.
This assumes, of course, that we're working our way anti-clockwise around
the polygon starting at vertex 0 and that the vertices are numbered in
proper sequence and that the polygon is convex, i.e. that EixE(i+1) also has
positive z components throughout.
If the polygon is not convex, this doesn't invalidate the method but
merely complicates it. You need to go into dentist mode, i.e. find and
fill all the cavities, also known as finding the convex hull of the original
polygon. Start at edge 0 and note whether edge 1 involves a left turn, i.e.
whether EixE(i+1) is positive. If it's zero, that's OK too, it simply
means V(i+1) is a funny vertex because edges i and i+1 form a straight line.
Anyway, if you find yourself having to hang a right at any vertex, then it
is in a cavity and should be marked "bad" and disregarded. Keep going round
until you find no more bad vertices during any whole round trip. The convex
hull consists of all the original vertices except those marked bad.
The next stage is to form a list of polygons which represent the cavities
you've filled. Do this by going around all the original vertixes until you
find one marked bad. The cavity is bounded by the last good vertex, the
bad one, any further bad ones adjacent to it, and the next good one.
The point is in the original polygon iff it is in its convex hull and not
in any of its cavities. So if the point is not in the polygon's convex
hull, the answer is no, but if it is then you need to check, recursively,
whether it is in any of its cavities. If it is in at least one of them,
the answer is still no, if it is in none of them, then the answer is yes
(i.e. the point is in the original polygon).
Since the cavities themselves may or may not be convex, your recursion may
need to go more than one level deep.
So there you are. It can be done this way, though it might be easier to
arrange for the vertex list to be generated in the first place as a simple
collection of convex portions, then you wouldn't need to faff around with
all this fun stuff, and could just boringly use the z-component-of-vector-
product algorithm on each of the portions. The answer would be yes iff it
is yes for any of them.
--
I don't know whether there really is a ukol.com.
In any case I'm not there, I'm at dcs.ed.ac.uk.
Lynn Killingbeck
Dick Fulton wrote:
>
> I'm looking for a method/algorithm to determine if a given point lies
> inside or outside of a predefined polygon. The point and the polygon
> vertices are latitude/longitude points on the earth's surface. The
> methods that come to mind seem difficult to implement as a computer
> program algorithm, which is what I need to do.
>
>
> Thanks...
>
> --
> Dick Fulton
> Palo Alto, CA
> email: dfu...@dinc.com
This is for a plane, but might give you some ideas.
"Computer Graphics Handbook" by Michael E. Mortenson, has a section on
'Point Containment Test'. If the polygon is convex, the idea is to use
a test-point which is known to be inside the polygon, such as the
midpoint of the X-Y extremes. Compute the signed distance to each edge,
for both your point and the test-point. If the signs are all the same,
then your point is also inside the polygon.
For a non-convex polygon, the idea is to determine the intersection of
any line (say, a horizontal line) through your point with all of the
edges of the polygon; then, pair the ordered intersections. If your
point is within any pair, then it's also within the polygon. I believe
there needs to be additional logic if you're exactly at a vertex, and
for any edge(s) parallel to your reference line.
Anyway, give a look here, or possibly in other standard references for
computer graphics.
By the way, on a sphere, whether or not the point is 'inside' depends,
really, on your viewpoint. For example, the polygon might be regarded
as the 'smaller area' - or, as the 'larger area' - both of which are
bounded by the identical lines.
Lynn Killingbeck
Brian Skinner
dfu...@dinc.com (Dick Fulton) wrote:
> I'm looking for a method/algorithm to determine if a given point lies
> inside or outside of a predefined polygon. The point and the polygon
> vertices are latitude/longitude points on the earth's surface. The
> methods that come to mind seem difficult to implement as a computer
> program algorithm, which is what I need to do.
I had this problem in a computer program over 30 years ago. My memory
is very poor but the basis of my method was something like:
P is the point and A[1], A[2] to A[n] are the vertices of the polygon
(in order). Eliminate the case where P is a vertex by explicit
checking.
Find the (signed) angle between the point and each successive side of
the polygon (using high precision routines) i.e. the angles
A[1]-P-A[2] to A[n]-P-A[1].
Add the angles up. If the answer is zero (with suitable tolerance)
then the point is outside. If the answer is +-2*pi (with suitable
tolerance) then the point is inside.
--
Brian Skinner (br...@brisk.demon.co.uk, http://www.brisk.demon.co.uk)