From: bsmith@mickey.NoSubdomain.NoDomain (Brian Smith)
Newsgroups: comp.programming
Subject: Re: How to decide if a point is inside a polygon?
Date: 29 Aug 1992 05:15:35 GMT
Organization: University of California, Berkeley
Distribution: world
NNTP-Posting-Host: mickey.cs.berkeley.edu

I got this off the net a while ago.  The original poster was Michael J.
Zehr, from MIT at the time.  It uses a clever version of the "scan in any
direction any see how many times the point crosses the polygon" test that, as
far as I know, handles the special cases correctly.

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

(in pseudo C:)

  edge[] is an array of edges
  an edge has two fields, end1 and end2, plus four precomputed values
to speed up the actual check: slope, intercept, minx and maxx
    slope = (end2.y - end1.y)/(end2.x - end1.x)
    intercept = end2.y - slope*end2.x
    minx = min(end1.x, end2.x)
    maxx = max(end1.x, end2.x)
  an "end" has two fields, x and y

{this isn't the same datastructure i'm using, so i'm paraphrasing
for simplicity and readability}

flag = FALSE;
for (k=0; k<num_edges; k++)
    if ((x > edge[k].end1.x) ^ (x > edge[k].end2.x))
        {
        if (y < (temp = edge[k].intercept + edge[k].slope*x)
            flag = !flag;
        else if (y==temp) 
            /* point is on an edge so do whatever you want with it */;
        }
    else if (edge[k].end1.x == edge[k].end2.x)
        if ((y > edge[k].end1.y) ^ (y > edge[k].end2.y))
            /* point is on an edge so do whatever you want with it */;

after it's checked each edge, flag == TRUE implies it's inside,
flag == FALSE implies it's outside the polygon.

the multiply only occurs when a vertical line extended from the point
being tested intersects with the edge. it then computes the intersection
point, and checks if the point is below the edge.  if so, it counts as
an intersection.  the multiply is only executed when there is a possible
intersection, which is zero, one or two times for a convex polygon. 
this works as-is for concave, but doesn't guarantee only two multiplies.
depending on your compiler and specific machine, it may be marginally 
faster to keep track of the lowest and highest end of each edge, and
put another if statement before the one that computes the actual
intersection, but i haven't tried this to see
then you'd have:

if (y < edge[k].min_y) flag^= TRUE;
else if (y <= edge[k].max_y) 
        if (y < (temp= ....

Also, if your polygon are usually short and fat then you'd want to 
change that to a horizontal line test instead of the current vertical
line test.  Also, if your polygons are guaranteed convex, you could count
the number of times the x tested true, and after the second one you
wouldn't have to look at any more edges.

BTW, there's a faster way to code the tests, at least for my
hardware/compiler combination:  precompute a max and min x and y for
each edge (uses up a lot of space, but does really help for a lot of
tests).  Then instead of:

((x > edge[k].end1.x) ^ (x > edge[k].end2.x))

use:

((x > edge[k].minx) && (x <= edge[k].maxx))

you avoid always having to do both compares, which is the where most of
the speed comes from.  in addition, some compilers don't produce good
code for the first code fragment.
 
-----
Brian C. Smith				arpa:  bsmith@cs.Berkeley.EDU
University of California, Berkeley	uucp: uunet!ucbvax!postgres!bsmith
Computer Sciences Department		phone: (510)642-9585