back to
contents
"Light Makes Right"
September 2, 1992
Volume 5, Number 3
Compiled by
All contents are copyright (c) 1992, all rights reserved by the individual authors
Archive locations: anonymous FTP at
ftp://ftp-graphics.stanford.edu/pub/Graphics/RTNews/,
wuarchive.wustl.edu:/graphics/graphics/RTNews, and many others.
You may also want to check out the Ray Tracing News issue guide and the ray tracing FAQ.
back to
contents
[To begin this issue, here is one of the better pieces of obfuscatory writing in the field of computer graphics. Reprinted from RTNews8 - EAH]
From: deb@charisma.graphics.cornell.edu
Newsgroups: comp.graphics
Keywords: P, NP, Jordan curve separation, Ursyhon Metrization Theorem
Organization: Program of Computer Graphics
In article [...] ncsmith@ndsuvax.UUCP (Timothy Lyle Smith) writes:
>
> I need to find a formula/algorithm to determine if a line intersects
> a polygon. I would prefer a method that would do this in as little
> time as possible. I need this for use in a forward raytracing
> program.
I think that this is a very difficult problem. To start with, lines and polygons are semi-algebraic sets which both contain uncountable number of points. Here are a few off-the-cuff ideas.
First, we need to check if the line and the polygon are separated. Now, the Jordan curve separation theorem says that the polygon divides the plane into exactly two open (and thus non-compact) regions. Thus, the line lies completely inside the polygon, the line lies completely outside the polygon, or possibly (but this will rarely happen) the line intersects the polyon.
Now, the phrasing of this question says "if a line intersects a polygon", so this is a decision problem. One possibility (the decision model approach) is to reduce the question to some other (well known) problem Q, and then try to solve Q. An answer to Q gives an answer to the original decision problem.
In recent years, many geometric problems have been successfully modeled in a new language called PostScript. (See "PostScript Language", by Adobe Systems Incorporated, ISBN # 0-201-10179-3, co. 1985).
So, given a line L and a polygon P, we can write a PostScript program that draws the line L and the polygon P, and then "outputs" the answer. By "output", we mean the program executes a command called "showpage", which actually prints a page of paper containing the line and the polygon. A quick examination of the paper provides an answer to the reduced problem Q, and thus the original problem.
There are two small problems with this approach.
(1) There is an infinite number of ways to encode L and P into the
reduced problem Q. So, we will be forced to invoke the Axiom of
Choice (or equivalently, Zorn's Lemma). But the use of the Axiom of
Choice is not regarded in a very serious light these days.
(2) More importantly, the question arises as to whether or not the
PostScript program Q will actually output a piece of paper; or in
other words, will it halt?
Now, PostScript is expressive enough to encode everything that a
Turing Machine might do; thus the halting problem (for PostScript) is
undecidable. It is quite possible that the original problem will turn
out to be undecidable.
I won't even begin to go into other difficulties, such as aliasing, finite precision and running out of ink, paper or both.
A couple of references might be:
1. Principia Mathematica. Newton, I. Cambridge University Press, Cambridge, England. (Sorry, I don't have an ISBN# for this).
2. An Introduction to Automata Theory, Languages, and Computation. Hopcroft, J and Ulman, J.
3. The C Programming Language. Kernighan, B and Ritchie, D.
4. A Tale of Two Cities. Dickens, C.
________
From: td@alice.UUCP (Tom Duff) Summary: Overkill. Organization: AT&T Bell Laboratories, Murray Hill NJ
The situation is not nearly as bleak as Baraff suggests (he should know better, he's hung around The Labs for long enough). By the well known Dobbin-Dullman reduction (see J. Dullman & D. Dobbin, J. Comp. Obfusc. 37,ii: pp. 33-947, lemma 17(a)) line-polygon intersection can be reduced to Hamiltonian Circuit, without(!) the use of Grobner bases, so LPI (to coin an acronym) is probably only NP-complete. Besides, Turing-completeness will no longer be a problem once our Cray-3 is delivered, since it will be able to complete an infinite loop in 4 milliseconds (with scatter-gather.)
________
From: deb@svax.cs.cornell.edu (David Baraff)
Well, sure it's no worse than NP-complete, but that's ONLY if you restrict yourself to the case where the line satisfies a Lipschitz condition on its second derivative. (I think there's an '89 SIGGRAPH paper from Caltech that deals with this).
back to
contents
What is the MOST EFFICIENT way to find if a point lies within a SIMPLE polygon given a set of vertices (in no particular order) and the point coordinate ? Are there anonymous ftp sites with such canned functions instead of re-inventing the wheel. Excuse me if that is a trivial question.
____
Philip Walden (pwalden@hpcc01.corp.hp.com) replies:
Another alternative to using a test ray is to sum the arc angles from the vertices to the test point. If the sum is 2pi, the point is inside, if 0 it's outside. i.e.
give a polygon with a set of vertices v[1:n] and test point p
if (SUM i=1 to n (arcangle(v[i]->p->v[i+1])) > pi
then p is inside
else it's outside
____
Eric Haines replies:
This method turns out to be exceedingly slow, about an order of magnitude slower than the "shoot a ray and count crossings" method. The problem is that the "arcangle" routine takes awhile: you need to normalize your vectors, you need to compute the dot product (for the angle) and cross product (for the sign of the angle), and you need the arccosine of the dot product for the angle itself. This is done for every vertex of the polygon. All of this adds up!
In fact, I wrote a program to test the two algorithms head-to-head when I was deciding which to use: code attached (for a switch!). The test program makes polygons from 3 to 7 (though you can change this) and random scatters the vertices from [0-1). It then tests each of these polygons with a series of random points from [0-1).
I won't claim that either implementation of the angle summation or segment crossing routines below are optimal (though I do try to share intermediate results between vertices). But I don't care how much you tune the angle summation algorithm: it's still going to be slower than crossing counting. Please do tell me of any hacks you make to get either algorithm faster, anyway.
Code's at the end of this issue; you should mess with *_TEST_TIMES however you see fit. I found you could set them pretty low, but I felt better when they were relatively high (i.e. the test times were greater than the granularity of the "time" command). The crossing counting inside/outside test is the one I wrote up in _An Introduction to Ray Tracing_, edited by Andrew Glassner, Academic Press, plus a trick or two not in the book.
I should also mention that the tests below assuming that an even winding number means that the point is outside, as this is how the Jordan Curve test is used. Think of the winding number as the number of loops around a point. As an example, think of a star shape: the interior pentagon is considered "outside" since points in it cross two lines. You can modify the tests as shown in the comments if you want to test for just a non-zero winding number (i.e. the star is fully filled).
Timings are shown after the next comment.
____
Tom Dickens notes in a related article:
Additional logic [in the barycentric test] can be added to avoid checking triangles which are located in a region which has already been determined to be beyond of a side of a previous triangle which returned a negative check.
[I don't know if such an enhancement is possible (one triangle's outside is another's inside) if the polygon is not convex, but mention it here for interest's sake; it is certainly not incorporated in my code. - EAH]
____
comments from Ron Capelli (capelli@vnet.ibm.com):
>the "arcangle" routine takes awhile: you need to normalize your vectors, you
>need to compute the dot product (for the angle) and cross product (for the
>sign of the angle), and you need the arccosine of the dot product for the
>angle itself. This is done for every vertex of the polygon.
See "Fast Approximation to the Arctangent" in Graphics Gems II (p.389). The maximum error of this approximation allows polygons with up to about 40 vertexes to be handled safely.
No need to normalize, take dot products, or cross products either. From the vector components (dx, dy) from a test point to a polygon vertex, get an angle using arctan (or the fast approximation, in which case do not scale to radians or degrees; use the angle value in the range 0 to 8). Calculate the signed angle between two successive vertexes by subtracting the angles for each vertex. Add up all the signed angle differences, and compare to the mid-range value.
The case of a test point coincident with a polygon vertex must be explicitly checked, as the arctangent for a zero-length vector is undefined. Another fine point is that angle differences must be wrapped in the range -180 to +180 degrees.
As fond as I am :-) of this method, I still won't argue that the "count ray crossings" method is not better. For a polygon with N edges, an optimized "count ray crossings" algorithm requires at most one multiply, one divide, and about a half-dozen add/subtract/compare operations per edge. The "add angle differences" algorithm using the fast approximation to the arctangent requires only one divide but more than twice as many add/subtract/compare operations per edge...
____
Sundar Narasimhan (sundar@fiber-one.ai.mit.edu) comments on Ron's post:
I think it is somewhat misleading to count worst case operations / edge and then multiply it blindly by N esp. since the "count ray crossings" algorithm can be optimized to consider only the subset of N edges that actually have a chance of intersecting the fired ray, whereas I haven't seen any such pre-filtering done to the add angle difference algorithm.
____
Christopher Jam (phillips@swanee.ee.uwa.oz.au) writes:
For the "count ray crossings" method, if your processor has no multiply or divide instructions, the multiply and divide may be replaced by a binary search for the point on the edge with the same y coordinate as the point being tested. Alternatively, if you are going to draw the polygon at some stage, generate a list of points on the line for plotting, and keep these in a table to completely remove any requirements for multiplies, divides or searches. Just look up a point on the same line ;)
____
Stuart MacMartin (sjm@bcrki65.bnr.ca) comments:
I know the purpose of your code was for a simple comparison of the two algorithms, but for production you might want to write optimal code.
The loop overhead is actually quite expensive, as I found out when I optimized someone else's implementation of this algorithm. Most of the time is spent looking for the next y that is on/below the ray or the next y that is above the ray. You should be able to get at lease a factor of 2 improvement by optimizing those two loops. Call a subroutine if necessary for the deleted stuff. (Note that subroutines can be MORE efficient because the compiler may make a better use of registers for the time-consuming part of the loop.) Use pointers to the y rather than to the x to avoid the extra adding, and don't keep track of so many variables in the loop.
Finally, is it necessary to use doubles? With appropriate units, will long suffice? Just be careful if you make this switch: if you have to have doubles in some places in your code, watch out for conversions from long to double and back. On some platforms (definitely Apollo; I don't know about others) such casts take so long that they totally swamp anything else that is inefficient in a loop.
My code for the loop:
stop = vertices + (numVertices << 1);
for(y = vertices[1], p = vertices+3; p < stop; y = *p,p+=2) {
/* Skip all lines above */
if (y > ay) while ((p < stop) && (*p > ay)) p += 2;
/* Skip all lines below */
else if (y < ay) while ((p < stop) && (*p < ay)) p += 2;
/* NOTE! Don't skip horizontal edges at ay! */
if (p >= stop) break;
The routine requires that the first point is duplicated as the last.
____
Final (for now) comments, by Eric Haines:
Other people noted the barycentric test as being a way to do inside/outside testing. The barycentric test is included below in Peter Shirley's article, and is also written up on page 390-393 of Graphics Gems, Didier Badouel's article. You can do better than figuring out all the intermediate differences (as shown in Badouel's article), instead computing them as needed. Also, my code worries about zero area triangles just to be safe, so it could be a bit faster without these sanity checks.
Badouel points out that triangle testing works for convex polygons: you simply test all of the triangles, and any intersection ends the test. Most people assume that concave triangles must be intersected using a different test. Not so! As Berlin pointed out (Vol 11, SIGGRAPH '85 course notes - sorry for the obscure reference!), you can test a concave polygon by checking all its triangles generated from one vertex. If the number of triangles hit is odd, the point is inside the (potentially concave) polygon - try it out, it actually works! You do have to test all triangles, however; if you knew the polygon was convex you can stop on the first intersection. Also, the barycentric coordinates become somewhat meaningless, since more than one triangle can be found to overlap your test point.
So, included in the new code (which is at the end of this issue) are two flavors of the crossings (segment) test, one which is fairly fast and readable, the other is "macmartinized" by me and runs a bit faster still. The barycentric polygon intersector is also included, and it works for all polygons (concave, self-intersecting, etc).
Random polygon testing:
number (or range) of vertices
3 4 5 3-7 20 100
angle 57.21 69.81 86.44 82.98 303.33 1435.78
barycentric 1.80 3.44 5.29 4.90 30.78 164.38
crossings 2.17 2.76 3.32 3.24 12.51 60.13
macmartin 1.86 2.35 2.87 2.64 11.04 52.27
inside % 6.00 12.00 13.78 9.33 26.89 41.22
Times are in microseconds per intersection test, on an HP 720 workstation. The "for" loop for repeating each test is counted into these test times (basically, I repeat each polygon/point test a number of times so that there's some reasonable amount of time to count) - not counting it just means the "crossings" tests is all that much faster. Basically, the macmartinized crossing test wins on my machine; your mileage may vary.
Note that the barycentric intersector could be faster if the polygons are known to be convex, since then a quick out could be taken (i.e. the first hit ends testing). I won't swear at all for the completeness of the barycentric code: I suspect you can get into roundoff error problems for intersection points which lay exactly on some internal, invisible triangle testing edge (rare, but possible).
There are also a few hacks I haven't done on the barycentric intersector: you could make the set of if statements into two big if statements, the variable "v0" is not all that useful, etc. Also, if the polygon is a triangle the test could be made a special case and the "for" loop and other baggage could be eliminated (this I tried, and it resulted in about a 20% reduction in run time). However, the code presented is at least readable, and the general behavior is the same: this intersector gets more expensive than the crossing methods as the number of polygon vertices rises.
Looking over the branch flow analysis, I find that my test case polygons (randomly generated points) cause a lot of algebraic tests to occur when using the crossing algorithms (i.e. the exact x intersection point has to be computed when a polygon edge goes from one quadrant to its diagonal opposite). For real polygons, this case is much less prevalent since vertices tend to be close to each other, so the polygons I've generated are more pathological than normally encountered. For example, if the polygons generated are regular polygons inscribed in a circle of radius 1.0 and given a random rotation about the origin, the timings are
Regular polygon testing:
number (or range) of vertices
3 4 5 3-7 20 100
angle 54.20 70.80 85.14 87.83 314.53 1522.00
barycentric 2.11 3.98 5.48 5.64 27.63 144.27
crossings 2.50 3.09 3.30 3.42 9.46 42.71
macmartin 2.32 2.81 2.91 2.92 6.29 23.91
inside % 33.22 51.00 60.22 55.22 77.44 80.22
I'm not sure why there's any noticeable difference in the "angle summation" tests, since this test should be the same for any test point (I suspect it's just the "times()" routine's granularity). The barycentric times change a little for better and worse, probably due to times() and to different polygon vs. point tests. The two crossings test noticeably improve for large polygons (29% and 54% savings, respectively) - this is due to less "difficult" edges, and for the macmartin test there are now longer series of edges fully above or fully below the test point. These tests show off the real improvement of the macmartin algorithm when testing many sided polygons.
The near-final word is that the crossings tests seems to be the most efficient overall, the angle test should generally be avoided like the plague, and the barycentric test is good for the special case of triangles (and you also get the vertex interpolant weights). Just to add to the confusion, I talk about intersecting complex polygons in an article below; if you're using a lot of these then another intersection algorithm might be a better route to go.
back to
contents
p = a*p0 + b*p2 + g*p3
and the constraint a + b + g = 1 will also hold. If you want to interpolate colors or normals across the triangle, you can use a similar formula. A little messing around will convince you that for points inside the triangle, a, b, and g are all positive, and outside the triangle, at least one is negative. We can take this info, and get rid of a:
p = p0 + b*(p1 - p0) + g*(p2-p0)
The point p is inside if b and g are positive, and (b+g) < 1.
A point on our ray can be represented as p' = o + t*v, where o is the point of origin on the ray, and v is the direction of propagation. Plugging this into the triangle formula, we can find whether any point on the ray is also a point on the plane:
o + t*v = p0 + b*(p1 - p0) + g*(p2-p0)
Note that this is really 3 linear equations (one for x, y and z). You can rewrite this as a 3 by 3 linear system and solve for (t, b, g) and there is a hit if and only if:
t > 0
b > 0
g > 0
b+g < 1
Alternatively, you can first find the equation of the plane containing the triangle and find t. This will give you an explicit point h on the plane. You now have:
h = p0 + b*(p1 - p0) + g*(p2-p0)
Now you have 3 equations and 2 unknowns. Choosing any two of the equations is equivalent to projecting the problem onto one of the cartesian planes. Watch out if the triangle is in one of the planes, because the projection might be a line segment. You can check the surface normal to see if it is one of those cases (I just use the dimensions associated with the smallest magnitude components of the normal).
Note-- none of this is new. GG II and the ray tracing book have other ways to do it. Also, if your scenes are large, you are probably better off optimizing your spatial partitioning code (the part which takes longer for bigger scenes).
[For more info, see Didier Badouel's article in "Graphics Gems", p. 390-393]
back to
contents
>For a project i was working on i needed a very fast 2-D point-in-polygon
>routine for a monte-carlo simulation, and came up with a variation of the
>standard jordan curve theorem algorithm which ran approximately on the
>order of O(n/20) when n > ~100 (for "mostly" round-ish polygons), instead
>of O(n). What I can't find, and no one here seems to know, is for large
>polygons what is the best 2-D algorithm. Everything they've heard of runs
>in O(n) (I know that n and n/20 are the same "order" but the you know what
>I mean). I'd like to know if there is something faster.
You might check Preparata & Shamos's _Computational Geometry_ book, I believe they talk about faster algorithms. Unfortunately, much of the stuff in CG tends to be things like O(K * log n), where K is something pretty large.
A simple way to test polygons with a large number of edges in O(n/K) time is to sort the edges into buckets, then find out which bucket the test point falls into. Then you do the usual ray test, but only against those edges in the bucket.
For example, you could take a hundred edge polygon and find the Y extent in the polygon's plane. Split this extent into, say, 20 buckets (or whatever number you like, more ==> more efficient, but also more memory and preprocessing). Now take each edge and note which buckets it's in. As an added speed-up, sort the edges by their leftmost X values in that bucket.
Now testing can begin. When you test a point, use the Y coordinate to find which bucket to check. Now test each edge by the normal Jordan test against the point. You only need to test edges to the left of the point (i.e. you are using a test ray in the -X direction), so as soon as you reach an edge whose leftmost X is greater than the point's X, you can stop testing.
The bucket sort immediately gets rid of a lot of edges, and the X test means often not having to test all edges in the bucket (it's optional icing).
If you're really into using memory and like preprocessing, you could use a grid structure to place buckets on the polygon. Each grid cell is inside, outside, or indeterminate. There are some tricks, but I leave these as an exercise for the reader... Anyway, this method can be very fast, since much of the polygon's area is classified as in or out and much of the time no edges have to be tested at all. The only limit is the grid cell resolution. As the resolution rises the actual testing approaches constant time ( O(1)!). There's more preprocessing involved, but if you're making a lot of tests then this is worth doing.
____
Benjamin Zhu (zhu@graphack.asd.sgi.com) comments on a similar question:
If your polygon is convex, (more generally, your polygon is star-shaped,) you can pre-process the polygon into triangle strips by connecting any point in the kernel of the polygon with each vertex of the polygon. This takes you O(n) pre-processing time. After that, for each point, you do a binary search in polar coordinates to locate the potential triangle where the point might lie. Then you can determine if the point lies within this triangle, which is trivial. See Preparata and Shamos' "Computational Geometry" for more details. Obviously, this algorithm will give you O(log(n)) time per point. So you might consider this one if you need to solve the point-inclusion problem many times.
back to
contents
/* Point in polygon inside/outside tester. Angle summation, barycentric * coordinates, and ray along x-axis (crossings testing) compared. * * copyright 1992 by Eric Haines, 3D/Eye Inc, erich@acm.org */
#include <math.h>
#include <sys/types.h>
#include <sys/param.h>
#include <sys/times.h>
#define X 0
#define Y 1
double drand48() ;
double AngleTimeTotal ;
double BarycentricTimeTotal ;
double CrossingsTimeTotal ;
double MacmartinTimeTotal ;
/* minimum & maximum number of polygon vertices to generate */
#define MIN_VERTS 3
#define MAX_VERTS 7
/* number of different polygons to try */
#define TEST_POLYGONS 30
/* number of different intersection points to try */
#define TEST_POINTS 30
/* number of times to try a single point vs. a polygon */
/* this should be > 1 / ( HZ * approx. single test time in seconds ) */
#define ANGLE_TEST_TIMES 1000
#define BARYCENTRIC_TEST_TIMES 10000
#define CROSSINGS_TEST_TIMES 10000
#define MACMARTIN_TEST_TIMES 10000
main(argc,argv)
int argc; char *argv[];
{
int i, j, k, numverts, inside_tot ;
int angle_flag, barycentric_flag, crossings_flag, macmartin_flag ;
double pgon[MAX_VERTS][2] ;
double point[2] ;
srand( 12345 ) ;
AngleTimeTotal = 0.0 ;
BarycentricTimeTotal = 0.0 ;
CrossingsTimeTotal = 0.0 ;
MacmartinTimeTotal = 0.0 ;
inside_tot = 0 ;
for ( i = 0 ; i < TEST_POLYGONS ; i++ ) {
#ifdef CENTERED_SQUARE
/* for debugging purposes, test against a square */
numverts = 4 ;
pgon[0][X] = 0.2 ;
pgon[0][Y] = 0.2 ;
pgon[1][X] = 0.8 ;
pgon[1][Y] = 0.2 ;
pgon[2][X] = 0.8 ;
pgon[2][Y] = 0.8 ;
pgon[3][X] = 0.2 ;
pgon[3][Y] = 0.8 ;
#else
/* make an arbitrary polygon fitting 0-1 range in x and y */
numverts = MIN_VERTS +
(int)(drand48() * (double)(MAX_VERTS-MIN_VERTS+1)) ;
for ( j = 0 ; j < numverts ; j++ ) {
pgon[j][X] = drand48() ;
pgon[j][Y] = drand48() ;
}
#endif
/* now try # of points against it */
for ( j = 0 ; j < TEST_POINTS ; j++ ) {
point[X] = drand48() ;
point[Y] = drand48() ;
angle_flag = angletest( pgon, numverts, point ) ;
barycentric_flag = barycentrictest( pgon, numverts, point ) ;
crossings_flag = crossingstest( pgon, numverts, point ) ;
macmartin_flag = macmartintest( pgon, numverts, point ) ;
/* reality check */
if ( angle_flag != crossings_flag ) {
printf( "angle test says %s, crossings test says %s\n",
angle_flag ? "INSIDE" : "OUTSIDE",
crossings_flag ? "INSIDE" : "OUTSIDE" ) ;
printf( "point %g %g\n", (float)point[X], (float)point[Y] ) ;
printf( "polygon:\n" ) ;
for ( k = 0 ; k < numverts ; k++ ) {
printf( " %g %g\n", (float)pgon[k][X], (float)pgon[k][Y]);
}
}
if ( barycentric_flag != macmartin_flag ) {
printf( "barycentric test says %s, macmartin test says %s\n",
barycentric_flag ? "INSIDE" : "OUTSIDE",
macmartin_flag ? "INSIDE" : "OUTSIDE" ) ;
printf( "point %g %g\n", (float)point[X], (float)point[Y] ) ;
printf( "polygon:\n" ) ;
for ( k = 0 ; k < numverts ; k++ ) {
printf( " %g %g\n", (float)pgon[k][X], (float)pgon[k][Y]);
}
}
inside_tot += crossings_flag ;
}
}
printf( "angle test time: %g microseconds per test\n",
(float)( AngleTimeTotal/(double)(TEST_POINTS*TEST_POLYGONS) ) ) ;
printf( "barycentric test time: %g microseconds per test\n",
(float)( BarycentricTimeTotal/(double)(TEST_POINTS*TEST_POLYGONS) ) ) ;
printf( "crossings test time: %g microseconds per test\n",
(float)( CrossingsTimeTotal/(double)(TEST_POINTS*TEST_POLYGONS) ) ) ;
printf( "macmartin crossings test time: %g microseconds per test\n",
(float)( MacmartinTimeTotal/(double)(TEST_POINTS*TEST_POLYGONS) ) ) ;
printf( "%g %% of all points were inside polygons\n",
(float)inside_tot * 100.0 / (float)(TEST_POINTS*TEST_POLYGONS) ) ;
}
/* sum angles of vtxN - point - vtxN+1, check if in PI to 3*PI range */
int
angletest( pgon, numverts, point )
double pgon[MAX_VERTS][2] ;
int numverts ;
double point[2] ;
{
int i, j, inside_flag ;
struct tms timebuf ;
long timestart ;
long timestop ;
double *vtx0, *vtx1, angle, len, vec0[2], vec1[2] ;
/* do the test a bunch of times to get a useful time reading */
timestart = times( &timebuf ) ;
for ( i = 0 ; i < ANGLE_TEST_TIMES ; i++ ) {
/* sum the angles and see if answer mod 2*PI > PI */
vtx0 = pgon[numverts-1] ;
vec0[X] = vtx0[X] - point[X] ;
vec0[Y] = vtx0[Y] - point[Y] ;
if ( (len = hypot( vec0[X], vec0[Y] )) <= 0.0 ) {
/* point and vertex coincide */
return( 1 ) ;
}
vec0[X] /= len ;
vec0[Y] /= len ;
angle = 0.0 ;
for ( j = 0 ; j < numverts ; j++ ) {
vtx1 = pgon[j] ;
vec1[X] = vtx1[X] - point[X] ;
vec1[Y] = vtx1[Y] - point[Y] ;
if ( (len = hypot( vec1[X], vec1[Y] )) <= 0.0 ) {
/* point and vertex coincide */
return( 1 ) ;
}
vec1[X] /= len ;
vec1[Y] /= len ;
/* check if vec1 is to "left" or "right" of vec0 */
if ( vec0[X] * vec1[Y] - vec1[X] * vec0[Y] >= 0.0 ) {
/* add angle due to dot product of vectors */
angle += acos( vec0[X] * vec1[X] + vec0[Y] * vec1[Y] ) ;
} else {
/* subtract angle due to dot product of vectors */
angle -= acos( vec0[X] * vec1[X] + vec0[Y] * vec1[Y] ) ;
}
/* get to next point */
vtx0 = vtx1 ;
vec0[X] = vec1[X] ;
vec0[Y] = vec1[Y] ;
}
/* test if between PI and 3*PI, 5*PI and 7*PI, etc */
/* if we care about is winding number > 0, then just:
inside_flag = fabs(angle) > M_PI ;
*/
inside_flag = fmod( fabs(angle) + M_PI, 4.0*M_PI ) > 2.0*M_PI ;
}
timestop = times( &timebuf ) ;
/* time in microseconds */
AngleTimeTotal += 1000000.0 * (double)(timestop - timestart) /
(double)(HZ * ANGLE_TEST_TIMES) ;
return (inside_flag) ;
}
/* barycentric, a la Gems I, with a little efficiency tuning */
int
barycentrictest( pgon, numverts, point )
double pgon[MAX_VERTS][2] ;
int numverts ;
double point[2] ;
{
int i, tris_hit, inside_flag, p1, p2 ;
struct tms timebuf ;
long timestart ;
long timestop ;
double tx, ty, u0, u1, u2, v0, v1, alpha, beta, denom ;
/* do the test a bunch of times to get a useful time reading */
timestart = times( &timebuf ) ;
for ( i = 0 ; i < BARYCENTRIC_TEST_TIMES ; i++ ) {
tx = point[X] ;
ty = point[Y] ;
tris_hit = 0 ;
for ( p1 = 1, p2 = 2 ; p2 < numverts ; p1++, p2++ ) {
if ( ( u1 = pgon[0][X] - pgon[p2][X] ) == 0.0 ) {
/* zero area test - optional */
if ( ( u2 = pgon[p1][X] - pgon[p2][X] ) == 0.0 ) {
goto NextTri;
}
/* Compute intersection point */
if ( ( ( beta = ( tx - pgon[p2][X] ) / u2 ) < 0.0 ) ||
( beta > 1.0 ) ) {
goto NextTri;
}
if ( ( v1 = pgon[0][Y] - pgon[p2][Y] ) == 0.0 ) {
goto NextTri;
}
if ( ( alpha = ( ty - pgon[p2][Y] - beta *
( pgon[p1][Y] - pgon[p2][Y] ) ) / v1 ) < 0.0 ) {
goto NextTri;
}
} else {
if ( ( denom = ( pgon[p1][Y] - pgon[p2][Y] ) * u1 -
( u2 = pgon[p1][X] - pgon[p2][X] ) *
( v1 = pgon[0][Y] - pgon[p2][Y] ) ) == 0.0 ) {
goto NextTri;
}
/* Compute intersection point & subtract 0 vertex */
u0 = tx - pgon[p2][X] ;
v0 = ty - pgon[p2][Y] ;
if ( ( ( ( beta = ( v0 * u1 - u0 * v1 ) / denom ) ) < 0.0 ) ||
( beta > 1.0 ) ) {
goto NextTri;
}
if ( ( alpha = ( u0 - beta * u2 ) / u1 ) < 0.0 ) {
goto NextTri;
}
}
/* check gamma */
if ( alpha + beta <= 1.0 ) {
/* survived */
tris_hit++ ;
}
NextTri: ;
}
inside_flag = tris_hit & 0x1 ;
}
timestop = times( &timebuf ) ;
/* time in microseconds */
BarycentricTimeTotal += 1000000.0 * (double)(timestop - timestart) /
(double)(HZ * BARYCENTRIC_TEST_TIMES) ;
return (inside_flag) ;
}
/* shoot a test ray along +X axis - slower version, less messy */
int
crossingstest( pgon, numverts, point )
double pgon[MAX_VERTS][2] ;
int numverts ;
double point[2] ;
{
int i, j, inside_flag, xflag0 ;
struct tms timebuf ;
long timestart ;
long timestop ;
double *vtx0, *vtx1, dv0 ;
int crossings, yflag0, yflag1 ;
/* do the test a bunch of times to get a useful time reading */
timestart = times( &timebuf ) ;
for ( i = 0 ; i < CROSSINGS_TEST_TIMES ; i++ ) {
vtx0 = pgon[numverts-1] ;
/* get test bit for above/below Y axis */
yflag0 = ( dv0 = vtx0[Y] - point[Y] ) >= 0.0 ;
crossings = 0 ;
for ( j = 0 ; j < numverts ; j++ ) {
/* cleverness: bobble between filling endpoints of edges, so
* that the previous edge's shared endpoint is maintained.
*/
if ( j & 0x1 ) {
vtx0 = pgon[j] ;
yflag0 = ( dv0 = vtx0[Y] - point[Y] ) >= 0.0 ;
} else {
vtx1 = pgon[j] ;
yflag1 = ( vtx1[Y] >= point[Y] ) ;
}
/* check if points not both above/below X axis - can't hit ray */
if ( yflag0 != yflag1 ) {
/* check if points on same side of Y axis */
if ( ( xflag0 = ( vtx0[X] >= point[X] ) ) ==
( vtx1[X] >= point[X] ) ) {
if ( xflag0 ) crossings++ ;
} else {
/* compute intersection of pgon segment with X ray, note
* if > point's X.
*/
crossings += (vtx0[X] -
dv0*( vtx1[X]-vtx0[X])/(vtx1[Y]-vtx0[Y])) >= point[X] ;
}
}
}
/* test if crossings is odd */
/* if we care about is winding number > 0, then just:
inside_flag = crossings > 0 ;
*/
inside_flag = crossings & 0x01 ;
}
timestop = times( &timebuf ) ;
/* time in microseconds */
CrossingsTimeTotal += 1000000.0 * (double)(timestop - timestart) /
(double)(HZ * CROSSINGS_TEST_TIMES) ;
return (inside_flag) ;
}
/* shoot a test ray along +X axis - macmartinized by me, a bit messier */
int
macmartintest( pgon, numverts, point )
double pgon[MAX_VERTS][2] ;
int numverts ;
double point[2] ;
{
int i, inside_flag, xflag0 ;
struct tms timebuf ;
long timestart ;
long timestop ;
double *p, *stop ;
int crossings ;
double tx, ty, y ;
/* do the test a bunch of times to get a useful time reading */
timestart = times( &timebuf ) ;
for ( i = 0 ; i < MACMARTIN_TEST_TIMES ; i++ ) {
crossings = 0 ;
tx = point[X] ;
ty = point[Y] ;
y = pgon[numverts-1][Y] ;
p = (double *)pgon + 1 ;
if ( ( y >= ty ) != ( *p >= ty ) ) {
/* x crossing */
if ( ( xflag0 = ( pgon[numverts-1][X] >= tx ) ) ==
( *(double *)pgon >= tx ) ) {
if ( xflag0 ) crossings++ ;
} else {
/* compute intersection of pgon segment with X ray, note
* if > point's X.
*/
crossings += ( pgon[numverts-1][X] -
(y-ty)*( *(double *)pgon - pgon[numverts-1][X])/(*p-y)) >= tx ;
}
}
stop = pgon[numverts] ;
for ( y=*p,p += 2 ; p < stop ; y=*p,p+=2) {
if ( y >= ty ) {
while ( (p < stop) && (*p >= ty) ) p+=2 ;
if ( p >= stop ) goto Exit ;
/* y crosses */
if ( ( xflag0 = ( *(p-3) >= tx ) ) ==
( *(p-1) >= tx ) ) {
if ( xflag0 ) crossings++ ;
} else {
/* compute intersection of pgon segment with X ray, note
* if > point's X.
*/
crossings += ( *(p-3) -
(*(p-2)-ty)*( *(p-1)-*(p-3))/(*p-*(p-2))) >= tx ;
}
} else {
while ( (p < stop) && (*p < ty)) p+=2 ;
if ( p >= stop ) goto Exit ;
/* y crosses */
if ( ( xflag0 = ( *(p-3) >= tx ) ) ==
( *(p-1) >= tx ) ) {
if ( xflag0 ) crossings++ ;
} else {
/* compute intersection of pgon segment with X ray, note
* if > point's X.
*/
crossings += ( *(p-3) -
(*(p-2)-ty)*( *(p-1)-*(p-3))/(*p-*(p-2))) >= tx ;
}
}
}
Exit:
/* test if crossings is odd */
/* if we care about is winding number > 0, then just:
inside_flag = crossings > 0 ;
*/
inside_flag = crossings & 0x01 ;
}
timestop = times( &timebuf ) ;
/* time in microseconds */
MacmartinTimeTotal += 1000000.0 * (double)(timestop - timestart) /
(double)(HZ * MACMARTIN_TEST_TIMES) ;
Eric Haines / erich@acm.org