jmn@cbnewsh.att.com (john.b.medamana) (07/28/90)
I am posting this article for a friend who works in the field of Urban Planning. I am looking for software which determines whether an arbitrary point lies within or outside a polygon. Please e-mail or post replys. Thanks in advance. John Medamana AT&T email: att!arch3!johnm
dga@cs.brown.edu (Daniel Gerardo Aliaga) (07/28/90)
In article <1990Jul27.201939.17816@cbnewsh.att.com> jmn@cbnewsh.att.com (john.b.medamana) writes: >I am posting this article for a friend who works in the field of >Urban Planning. > >I am looking for software which determines whether an arbitrary >point lies within or outside a polygon. Please e-mail or post >replys. > >Thanks in advance. > >John Medamana >AT&T >email: att!arch3!johnm One very simple algorithm (from some paper in SigGraph '89 I believe) is take the dot product of the normal of each side of the polygon and the vector formed by the point in question and a vertex of the current side. If all of these dot products are negative, the point is inside the object. This of course can be generalized to the 3-D case. The disadvantage is that this only works for convex polygons. graphically: | n | <---| /p | / |/v If all the dot products <n,(p-v)> < 0, then p is inside convex polygon. Alternatively, you could use an even odd rule if you have non convex polygons. Imagine the 2-D case, take a polygon and draw a line through the point. If that line intersects an even number of sides, it is inside else outside. This may produce some hairy cases when the point is right on a side, or when the line you draw coincides exactly with one of the sides. You will have to check for these cases. graphically: _____ / \ --- \_____ | \ | --- | xxxxxxxxx p xxxxxxxxxxxxxxxxx \_______/ \___| Notice that the line crosses the sides. It also coincides with a side. These should work, good luck ! P.S. I have software for the first method, though it uses a in house polyhedra representation package, I could give you a copy ... XXXXXXXXXXXXXXX XXXXXXXXXXXXXXXXXX Daniel G. Aliaga '91 XXX XXX XXXXX Brown University XXXXXXXXXXXXXXXXXX XXXXXXXXXXXXXXXXX "Real men write self modifying code"
BOYDJ@QUCDN.QueensU.CA (Jeff Boyd) (07/29/90)
What information have you on the polygon ? Eg. how is it defined ? Do you have a set of ordered co-ordinates for the vertices ? Is the polygon always convex ?
gordon@cs.tamu.edu (Dan Gordon) (07/30/90)
In article <46093@brunix.UUCP] dga@cs.brown.edu (Daniel Gerardo Aliaga) writes: ]In article <1990Jul27.201939.17816@cbnewsh.att.com> jmn@cbnewsh.att.com (john.b.medamana) writes: ]>I am posting this article for a friend who works in the field of ]>Urban Planning. ]> ]>I am looking for software which determines whether an arbitrary ]>point lies within or outside a polygon. Please e-mail or post ]>replys. ]> ]>John Medamana ]>AT&T ]>email: att!arch3!johnm ] ]One very simple algorithm (from some paper in SigGraph '89 I believe) is ]take the dot product of the normal of each side of the polygon and the ]vector formed by the point in question and a vertex of the current side. ]If all of these dot products are negative, the point is inside the object. ]This of course can be generalized to the 3-D case. The disadvantage is that ]this only works for convex polygons. In 2-D, this algorithm takes time O(N), where N is the no. of points. If you assume that the coordinates of the polygon are in arrays, then there is a method of doing it in O(logN) time. It is only advantageous if your polygons have "lots" of vertices, where "lots" depends on your computing environment and on how often you need to solve this problem.
gordon@cs.tamu.edu (Dan Gordon) (07/30/90)
I forgot to mention in my previous post that the logarithmic-time solution works only for convex polygons.
phil@hpsmdca.HP.COM (Philip Walden) (07/30/90)
For 2D, I sum the absolute angles from THE POINT to all the vertices on the polygon. It the sum will be either 2pi: point is inside 0: point is outside
gordon@cs.tamu.edu (Dan Gordon) (07/31/90)
Subject: Re:: Algorithm needed to determine if a point is inside a polygon >> I forgot to mention... >> works only for convex polygons. >You also forgot to mention the name of the algorithm and where to obtain >information on it. I, for one, would be interested in finding this out. >Thanks, >Lem Davis >lem@nosc.mil Oh, dear; me and my big mouth. To begin with, turn to "Computational Geometry" by Preparata & Shamos, Section 2.2.1. Theorem 2.2 shows how to do it in O(logN) time using O(N) preprocessing time. The idea in this proof is to divide the polygon into wedges and do a circular binary search on the wedges. The division into wedges takes O(N) time, but it is a one-time piece of work. The actual binary search, which can be expected to be repeated many times, takes O(logN) time. My variation on this scheme works on arrays. We assume that the coordinates are in arrays X[1..N], Y[1..N], with the vertices of the polygon in circular order (clockwise or anticlockwise). Preprocessing step: find index of vertex with maximal Y, and index of vertex with minimal Y. This can be easily done in O(N) time, and this is probably good enough for most applications. However, it can also be done in O(logN) time. Query step: Given a point (a,b), is it in the polygon? Geometrically, this is done by the parity method: consider the horizontal line y=b. If the line is above the maximal Y or below the minimal Y, the answer is no. Otherwise, find the edges of the polygon which the line intersects. Count the no. of intersections to the left of (a,b) - if it is one, then the point is inside; otherwise (0 or 2), it is outside. The logarithmic time is obtained by finding the edges intersecting y=b by doing a circular binary on the array Y. Assume for simplicity that the minimum Y is at Y[1] and the maximum is Y[m], where 1<m<=N. One intersection is found by doing a binary search of b among Y[1..m], and the other by doing a binary search on Y[m],Y[m+1],...,Y[N],Y[1]. Y[m] * / \ / \ * * | \ --------------------------*----------- y=b | \ (a,b) | \ Y[2]* *Y[N] \ / \ / \ / \ / \* Y[1] I hope this helps.
jcs@crash.cts.com (John Schultz) (07/31/90)
/*
Here's a line/convex polygon intersection routine. I haven't fully tested
it, so if you spot a bug, let me know.
This code is set up for 32 bit integer fixed point, but can easily be
modified for floats by removing the shifts and declaring floats.
Another trick is to rotate and translate the *ray* (two points) and not
the polygons/normals. The ray transform matrix is the inverse of the
polygon (local) transform. See the Siggraph article for more info.
John
*/
/* Collision.c */
/* Created by John C. Schultz 18-July-90 */
/* Idea from Siggraph Computer Graphics Vol 21, Number 4, July '87,
page 126 */
#include <math.h>
#include <stdio.h>
#define sgn(x) ((x)>0 ? 1: ((x)<0 ? -1: (0)))
#define FBITS 14
typedef struct point3d {long defx,defy,defz,x,y,z;} point3d;
/* Prototypes */
short linepolyint(point3d * ip,point3d * b,point3d * d,point3d * n,
point3d * poly,short count);
/* Functions */
short linepolyint(point3d * ip,point3d * b,point3d * d,point3d * n,
point3d * poly,short count) {
long k,t; /* k = Plane constant, t = parametric var */
long sn; /* sign max normal */
long absx,absy,absz; /* abs x,y,z of normal */
long nd; /* N dot D */
long beta; /* Result of cross product */
short i,i1;
point3d p,q; /* temps */
/* Compute:
t = parametric variable
k = plane constant.
N = normal (plane equation).
B,D = line base, direction.
k - (N dot B)
t = -----------
(N dot D)
*/
/* Compute plane contant k */
k = poly[0].x*n->x + poly[0].y*n->y + poly[0].z*n->z; /* Keep upshifted */
nd = (n->x*d->x + n->y*d->y + n->z*d->z) >> FBITS;
if (nd) {
/* Compute parametric variable */
t = (k - (n->x*b->x + n->y*b->y + n->z*b->z)) / nd;
} else {
return 0; /* Line is parallel to plane */
} /* if nd */
/* Now see if point is within polygon */
absx = abs(n->x);
absy = abs(n->y);
absz = abs(n->z);
/* Find max of normal{x,y,z}, do specific case */
if (absx > absy) {
if (absx > absz) { /************ X CASE ************/
sn = sgn(n->x);
/* Compute only intersection points being checked */
ip->y = b->y + ((t*d->y) >> FBITS);
ip->z = b->z + ((t*d->z) >> FBITS);
for (i=0; i<count; i++) {
if ( (i1 = i+1) >= count) i1 = 0; /* Modulo count */
p.y = poly[i1].y - poly[i].y;
p.z = poly[i1].z - poly[i].z;
q.y = ip->y - poly[i].y;
q.z = ip->z - poly[i].z;
beta = (p.y*q.z - p.z*q.y); /* cross product: makes an x */
/* beta == 0 when on polygon boundary */
if ((beta) && (sgn(beta) != sn)) return 0;
/* The above line may need to be modified to implement:
0 <= beta <= 1, if n->x is positive and
-1 <= beta <= 0, if n->x is negative.
I simplified it for speed.
The above code occurs twice more below.
*/
}/* for i */
/* Success, now compute last point */
ip->x = b->x + ((t*d->x) >> FBITS);
} else {
goto zmax; /* For speed and non-redunancy */
} /* if x */
} else {
if (absy > absz) { /************ Y CASE ************/
sn = sgn(n->y);
/* Compute only intersection points being checked */
ip->x = b->x + ((t*d->x) >> FBITS);
ip->z = b->z + ((t*d->z) >> FBITS);
for (i=0; i<count; i++) {
if ( (i1 = i+1) >= count) i1 = 0;
p.x = poly[i1].x - poly[i].x;
p.z = poly[i1].z - poly[i].z;
q.x = ip->x - poly[i].x;
q.z = ip->z - poly[i].z;
beta = (p.z*q.x - p.x*q.z); /* cross product: makes a y */
/* beta == 0 when on polygon boundary */
if ((beta) && (sgn(beta) != sn)) return 0;
}/* for i */
/* Success, now compute last point */
ip->y = b->y + ((t*d->y) >> FBITS);
} else {
zmax: /************ Z CASE ************/
sn = sgn(n->z);
/* Compute only intersection points being checked */
ip->x = b->x + ((t*d->x) >> FBITS);
ip->y = b->y + ((t*d->y) >> FBITS);
for (i=0; i<count; i++) {
if ( (i1 = i+1) >= count) i1 = 0;
p.x = poly[i1].x - poly[i].x;
p.y = poly[i1].y - poly[i].y;
q.x = ip->x - poly[i].x;
q.y = ip->y - poly[i].y;
beta = (p.x*q.y - p.y*q.x); /* cross product: makes a z */
/* beta == 0 when on polygon boundary */
if ((beta) && (sgn(beta) != sn)) return 0;
}/* for i */
/* Success, now compute last point */
ip->z = b->z + ((t*d->z) >> FBITS);
} /* if y/z */
} /* if x/y */
return 1; /* Line intersected polygon */
} /* end raypolyint */
void main(long argc,char ** argv) {
point3d b,d,n,ip,p1,p2;
point3d poly[3];
short count = 3;
short hit;
/* Initialize intersection point */
ip.x = 0;
ip.y = 0;
ip.z = 0;
/* Start of line */
p1.x = 0;
p1.y = -1000;
p1.z =-9992;
/* End of line */
p2.x = 0;
p2.y = 1000;
p2.z = -10002;
/* Normal */
n.x = 0;
n.y = -((1<<FBITS)-1); /* -1.0 */
n.z = 0;
/* Line base is p1 */
b = p1;
/* Line direction is p2 - p1 */
d.x = p2.x - p1.x;
d.y = p2.y - p1.y;
d.z = p2.z - p1.z;
/* Define polygon */
poly[0].x = 0;
poly[0].y = 500;
poly[0].z = 10000;
poly[2].x = 10000;
poly[2].y = 500;
poly[2].z = -10000;
poly[1].x = -10000;
poly[1].y = 500;
poly[1].z = -10000;
/* Test the code */
hit = linepolyint(&ip,&b,&d,&n,poly,count);
printf("hit %d, x %d y %d z %d\n",hit,ip.x,ip.y,ip.z);
} /* end main */
/* end collision.c */