home *** CD-ROM | disk | FTP | other *** search
/ NetNews Usenet Archive 1992 #19 / NN_1992_19.iso / spool / comp / graphics / 9429 < prev    next >
Encoding:
Text File  |  1992-09-03  |  37.9 KB  |  1,012 lines

  1. Path: sparky!uunet!cis.ohio-state.edu!rutgers!psinntp!psinntp!eye!erich
  2. From: erich@eye.com (Eric Haines)
  3. Newsgroups: comp.graphics
  4. Subject: Ray Tracing News, Volume 5, Number 3
  5. Summary: point in polygon revisited
  6. Message-ID: <1992Sep2.101609.9105@eye.com>
  7. Date: 2 Sep 92 14:16:09 GMT
  8. Sender: erich@eye.com (Eric Haines)
  9. Organization: 3D/EYE, Inc.  Ithaca, NY
  10. Lines: 1000
  11.  
  12.  
  13.  _ __                 ______                         _ __
  14. ' )  )                  /                           ' )  )
  15.  /--' __.  __  ,     --/ __  __.  _. o ____  _,      /  / _  , , , _
  16. /  \_(_/|_/ (_/_    (_/ / (_(_/|_(__<_/ / <_(_)_    /  (_</_(_(_/_/_)_
  17.              /                               /|
  18.             '                               |/
  19.  
  20.             "Light Makes Right"
  21.  
  22.              September 2, 1992
  23.             Volume 5, Number 3
  24.  
  25. Compiled by Eric Haines, 3D/Eye Inc, 2359 Triphammer Rd, Ithaca, NY 14850
  26.     erich@eye.com
  27. All contents are copyright (c) 1992 by the individual authors
  28. Archive locations: anonymous FTP at princeton.edu (128.112.128.1)
  29.     /pub/Graphics/RTNews, wuarchive.wustl.edu:/graphics/graphics/RTNews,
  30.     and many others.
  31. UUCP archive access: write Kory Hamzeh (quad.com!avatar!kory) for info.
  32.  
  33. Contents:
  34.     Introduction
  35.     Intersection Between a Line and a Polygon (UNDECIDABLE??), by Dave Baraff,
  36.     Tom Duff
  37.     Fastest Point in Polygon Test, by Aladdin Nassar, Philip Walden, Eric
  38.     Haines, Tom Dickens, Ron Capelli, Sundar Narasimhan, Christopher Jam,
  39.     and (last but not least) Stuart MacMartin
  40.     Polygon Intersection via Barycentric Coordinates, by Peter Shirley
  41.     Many-Sided Polygon Intersection, by Eric Haines, Benjamin Zhu
  42.     Code for Point in Polygon Intersectors, by Eric Haines
  43.  
  44. -------------------------------------------------------------------------------
  45.  
  46. Introduction
  47.  
  48. This special issue is dedicated to everyone's (well, at least my) favorite
  49. computer graphics FAQ:  how do you find if a given point is inside a polygon?
  50. John Griegg's FAQ posting points at _An Introduction to Ray Tracing_, edited
  51. by Andrew Glassner.  This issue speeds up that algorithm (count the crossings
  52. of a test ray), and hopefully lays to rest the idea of using the "sum the
  53. angles" algorithm.  Plus there's a little discussion of what to do for special
  54. cases such as triangles and complex polygons, and code for four different ways
  55. to do this test at the end of the issue.
  56.  
  57. -------------------------------------------------------------------------------
  58.  
  59. [To begin this issue, here is one of the better pieces of obfuscatory writing
  60. in the field of computer graphics.  Reprinted from RTNews8 - EAH]
  61.  
  62.  
  63. Intersection Between a Line and a Polygon (UNDECIDABLE??),
  64.     by Dave Baraff, Tom Duff
  65.  
  66.     From: deb@charisma.graphics.cornell.edu
  67.     Newsgroups: comp.graphics
  68.     Keywords: P, NP, Jordan curve separation, Ursyhon Metrization Theorem
  69.     Organization: Program of Computer Graphics
  70.  
  71. In article [...] ncsmith@ndsuvax.UUCP (Timothy Lyle Smith) writes:
  72. >
  73. >  I need to find a formula/algorithm to determine if a line intersects
  74. >  a polygon.  I would prefer a method that would do this in as little
  75. >  time as possible.  I need this for use in a forward raytracing
  76. >  program.
  77.  
  78. I think that this is a very difficult problem.  To start with, lines and
  79. polygons are semi-algebraic sets which both contain uncountable number of
  80. points.  Here are a few off-the-cuff ideas.
  81.  
  82. First, we need to check if the line and the polygon are separated.  Now, the
  83. Jordan curve separation theorem says that the polygon divides the plane into
  84. exactly two open (and thus non-compact) regions.  Thus, the line lies
  85. completely inside the polygon, the line lies completely outside the polygon,
  86. or possibly (but this will rarely happen) the line intersects the polyon.
  87.  
  88. Now, the phrasing of this question says "if a line intersects a polygon", so
  89. this is a decision problem.  One possibility (the decision model approach) is
  90. to reduce the question to some other (well known) problem Q, and then try to
  91. solve Q.  An answer to Q gives an answer to the original decision problem.
  92.  
  93. In recent years, many geometric problems have been successfully modeled in a
  94. new language called PostScript.  (See "PostScript Language", by Adobe Systems
  95. Incorporated, ISBN # 0-201-10179-3, co. 1985).
  96.  
  97. So, given a line L and a polygon P, we can write a PostScript program that
  98. draws the line L and the polygon P, and then "outputs" the answer.  By
  99. "output", we mean the program executes a command called "showpage", which
  100. actually prints a page of paper containing the line and the polygon.  A quick
  101. examination of the paper provides an answer to the reduced problem Q, and thus
  102. the original problem.
  103.  
  104. There are two small problems with this approach.
  105.  
  106.     (1) There is an infinite number of ways to encode L and P into the
  107.     reduced problem Q.  So, we will be forced to invoke the Axiom of
  108.     Choice (or equivalently, Zorn's Lemma).  But the use of the Axiom of
  109.     Choice is not regarded in a very serious light these days.
  110.  
  111.     (2) More importantly, the question arises as to whether or not the
  112.     PostScript program Q will actually output a piece of paper; or in
  113.     other words, will it halt?
  114.  
  115.     Now, PostScript is expressive enough to encode everything that a
  116.     Turing Machine might do; thus the halting problem (for PostScript) is
  117.     undecidable.  It is quite possible that the original problem will turn
  118.     out to be undecidable.
  119.  
  120.  
  121. I won't even begin to go into other difficulties, such as aliasing, finite
  122. precision and running out of ink, paper or both.
  123.  
  124. A couple of references might be:
  125.  
  126. 1. Principia Mathematica.  Newton, I.  Cambridge University Press, Cambridge,
  127.    England.  (Sorry, I don't have an ISBN# for this).
  128.  
  129. 2. An Introduction to Automata Theory, Languages, and Computation.  Hopcroft, J
  130.   and Ulman, J.
  131.  
  132. 3. The C Programming Language. Kernighan, B and Ritchie, D.
  133.  
  134. 4. A Tale of Two Cities. Dickens, C.
  135.  
  136. --------
  137.  
  138. From: td@alice.UUCP (Tom Duff)
  139. Summary: Overkill.
  140. Organization: AT&T Bell Laboratories, Murray Hill NJ
  141.  
  142. The situation is not nearly as bleak as Baraff suggests (he should know
  143. better, he's hung around The Labs for long enough).  By the well known
  144. Dobbin-Dullman reduction (see J. Dullman & D. Dobbin, J. Comp. Obfusc.
  145. 37,ii:  pp. 33-947, lemma 17(a)) line-polygon intersection can be reduced to
  146. Hamiltonian Circuit, without(!) the use of Grobner bases, so LPI (to coin an
  147. acronym) is probably only NP-complete.  Besides, Turing-completeness will no
  148. longer be a problem once our Cray-3 is delivered, since it will be able to
  149. complete an infinite loop in 4 milliseconds (with scatter-gather.)
  150.  
  151. --------
  152.  
  153. From: deb@svax.cs.cornell.edu (David Baraff)
  154.  
  155. Well, sure it's no worse than NP-complete, but that's ONLY if you restrict
  156. yourself to the case where the line satisfies a Lipschitz condition on its
  157. second derivative.  (I think there's an '89 SIGGRAPH paper from Caltech that
  158. deals with this).
  159.  
  160. -------------------------------------------------------------------------------
  161.  
  162. Fastest Point in Polygon Test, by Aladdin Nassar, Philip Walden, Eric Haines,
  163.     Tom Dickens, Ron Capelli, Sundar Narasimhan, Christopher Jam, and
  164.     (last but not least) Stuart MacMartin
  165.  
  166.  
  167. Aladdin Nassar asked a variant on the classic question:
  168.  
  169. What is the MOST EFFICIENT way to find if a point lies within a SIMPLE
  170. polygon given a set of vertices (in no particular order) and the point
  171. coordinate ?   Are there anonymous ftp sites with such canned functions
  172. instead of re-inventing the wheel.  Excuse me if that is a trivial question.
  173.  
  174. ----
  175.  
  176. Philip Walden (pwalden@hpcc01.corp.hp.com) replies:
  177.  
  178. Another alternative to using a test ray is to sum the arc angles from
  179. the vertices to the test point. If the sum is 2pi, the point is inside,
  180. if 0 it's outside. i.e.
  181.  
  182.     give a polygon with a set of vertices v[1:n] an test point p
  183.  
  184.     if (SUM i=1 to n (arcangle(v[i]->p->v[i+1])) > pi
  185.     then p is inside
  186.     else it's outside
  187.  
  188. ----
  189.  
  190. Eric Haines replies:
  191.  
  192. This method turns out to be exceedingly slow, about an order of magnitude
  193. slower than the "shoot a ray and count crossings" method.  The problem is that
  194. the "arcangle" routine takes awhile:  you need to normalize your vectors, you
  195. need to compute the dot product (for the angle) and cross product (for the
  196. sign of the angle), and you need the arccosine of the dot product for the
  197. angle itself.  This is done for every vertex of the polygon.  All of this adds
  198. up!
  199.  
  200. In fact, I wrote a program to test the two algorithms head-to-head when I was
  201. deciding which to use:  code attached (for a switch!).  The test program makes
  202. polygons from 3 to 7 (though you can change this) and random scatters the
  203. vertices from [0-1).  It then tests each of these polygons with a series of
  204. random points from [0-1).
  205.  
  206. I won't claim that either implementation of the angle summation or segment
  207. crossing routines below are optimal (though I do try to share intermediate
  208. results between vertices).  But I don't care how much you tune the angle
  209. summation algorithm:  it's still going to be slower than crossing counting.
  210. Please do tell me of any hacks you make to get either algorithm faster,
  211. anyway.
  212.  
  213. Code's at the end of this issue; you should mess with *_TEST_TIMES however you
  214. see fit.  I found you could set them pretty low, but I felt better when they
  215. were relatively high (i.e. the test times were greater than the granularity
  216. of the "time" command).  The crossing counting inside/outside test is the one
  217. I wrote up in _An Introduction to Ray Tracing_, edited by Andrew Glassner,
  218. Academic Press, plus a trick or two not in the book.
  219.  
  220. I should also mention that the tests below assuming that an even winding
  221. number means that the point is outside, as this is how the Jordan Curve test
  222. is used.  Think of the winding number as the number of loops around a point.
  223. As an example, think of a star shape:  the interior pentagon is considered
  224. "outside" since points in it cross two lines.  You can modify the tests as
  225. shown in the comments if you want to test for just a non-zero winding number
  226. (i.e. the star is fully filled).
  227.  
  228. Timings are shown after the next comment.
  229.  
  230. ----
  231.  
  232. Tom Dickens notes in a related article:
  233.  
  234. Additional logic [in the barycentric test] can be added to avoid checking
  235. triangles which are located in a region which has already been determined to
  236. be beyond of a side of a previous triangle which returned a negative check.
  237.  
  238. [I don't know if such an enhancement is possible (one triangle's outside is
  239. another's inside) if the polygon is not convex, but mention it here for
  240. interest's sake; it is certainly not incorporated in my code.  - EAH]
  241.  
  242. ----
  243.  
  244. comments from Ron Capelli (capelli@vnet.ibm.com):
  245.  
  246. >the "arcangle" routine takes awhile:  you need to normalize your vectors, you
  247. >need to compute the dot product (for the angle) and cross product (for the
  248. >sign of the angle), and you need the arccosine of the dot product for the
  249. >angle itself.  This is done for every vertex of the polygon.
  250.  
  251. See "Fast Approximation to the Arctangent" in Graphics Gems II (p.389).  The
  252. maximum error of this approximation allows polygons with up to about 40
  253. vertexes to be handled safely.
  254.  
  255. No need to normalize, take dot products, or cross products either.  From the
  256. vector components (dx, dy) from a test point to a polygon vertex, get an angle
  257. using arctan (or the fast approximation, in which case do not scale to radians
  258. or degrees; use the angle value in the range 0 to 8).  Calculate the signed
  259. angle between two successive vertexes by subtracting the angles for each
  260. vertex.  Add up all the signed angle differences, and compare to the mid-range
  261. value.
  262.  
  263. The case of a test point coincident with a polygon vertex must be explicitly
  264. checked, as the arctangent for a zero-length vector is undefined.  Another
  265. fine point is that angle differences must be wrapped in the range -180 to +180
  266. degrees.
  267.  
  268. As fond as I am :-) of this method, I still won't argue that the "count ray
  269. crossings" method is not better.  For a polygon with N edges, an optimized
  270. "count ray crossings" algorithm requires at most one multiply, one divide, and
  271. about a half-dozen add/subtract/compare operations per edge.  The "add angle
  272. differences" algorithm using the fast approximation to the arctangent requires
  273. only one divide but more than twice as many add/subtract/compare operations
  274. per edge...
  275.  
  276. ----
  277.  
  278. Sundar Narasimhan (sundar@fiber-one.ai.mit.edu) comments on Ron's post:
  279.  
  280. I think it is somewhat misleading to count worst case operations / edge and
  281. then multiply it blindly by N esp.  since the "count ray crossings" algorithm
  282. can be optimized to consider only the subset of N edges that actually have a
  283. chance of intersecting the fired ray, whereas I haven't seen any such
  284. pre-filtering done to the add angle difference algorithm.
  285.  
  286. ----
  287.  
  288. Christopher Jam (phillips@swanee.ee.uwa.oz.au) writes:
  289.  
  290. For the "count ray crossings" method, if your processor has no multiply or
  291. divide instructions, the multiply and divide may be replaced by a binary
  292. search for the point on the edge with the same y coordinate as the point being
  293. tested.  Alternatively, if you are going to draw the polygon at some stage,
  294. generate a list of points on the line for plotting, and keep these in a table
  295. to completely remove any requirements for multiplies, divides or searches.
  296. Just look up a point on the same line ;)
  297.  
  298. ----
  299.  
  300. Stuart MacMartin (sjm@bcrki65.bnr.ca) comments:
  301.  
  302. I know the purpose of your code was for a simple comparison of the two
  303. algorithms, but for production you might want to write optimal code.
  304.  
  305. The loop overhead is actually quite expensive, as I found out when I optimized
  306. someone else's implementation of this algorithm.  Most of the time is spent
  307. looking for the next y that is on/below the ray or the next y that is above
  308. the ray.  You should be able to get at lease a factor of 2 improvement by
  309. optimizing those two loops.  Call a subroutine if necessary for the deleted
  310. stuff.  (Note that subroutines can be MORE efficient because the compiler may
  311. make a better use of registers for the time-consuming part of the loop.)  Use
  312. pointers to the y rather than to the x to avoid the extra adding, and don't
  313. keep track of so many variables in the loop.
  314.  
  315. Finally, is it necessary to use doubles?  With appropriate units, will long
  316. suffice?  Just be careful if you make this switch:  if you have to have
  317. doubles in some places in your code, watch out for conversions from long to
  318. double and back.  On some platforms (definitely Apollo; I don't know about
  319. others) such casts take so long that they totally swamp anything else that is
  320. inefficient in a loop.
  321.  
  322. My code for the loop:
  323.  
  324.     stop = vertices + (numVertices << 1);
  325.  
  326.     for(y = vertices[1], p = vertices+3;  p < stop;  y = *p,p+=2) {
  327.     /* Skip all lines above */
  328.     if      (y > ay) while ((p < stop) && (*p > ay)) p += 2;
  329.     /* Skip all lines below */
  330.     else if (y < ay) while ((p < stop) && (*p < ay)) p += 2;
  331.         /* NOTE! Don't skip horizontal edges at ay! */
  332.     if (p >= stop) break;
  333.  
  334. The routine requires that the first point is duplicated as the last.
  335.  
  336. ----
  337.  
  338. Final (for now) comments, by Eric Haines:
  339.  
  340. Other people noted the barycentric test as being a way to do inside/outside
  341. testing.  The barycentric test is included below in Peter Shirley's article,
  342. and is also written up on page 390-393 of Graphics Gems, Didier Badouel's
  343. article.  You can do better than figuring out all the intermediate differences
  344. (as shown in Badouel's article), instead computing them as needed.  Also, my
  345. code worries about zero area triangles just to be safe, so it could be a bit
  346. faster without these sanity checks.
  347.  
  348. Badouel points out that triangle testing works for convex polygons:  you
  349. simply test all of the triangles, and any intersection ends the test.  Most
  350. people assume that concave triangles must be intersected using a different
  351. test.  Not so!  As Berlin pointed out (Vol 11, SIGGRAPH '85 course notes -
  352. sorry for the obscure reference!), you can test a concave polygon by checking
  353. all its triangles generated from one vertex.  If the number of triangles hit
  354. is odd, the point is inside the (potentially concave) polygon - try it out, it
  355. actually works!  You do have to test all triangles, however; if you knew the
  356. polygon was convex you can stop on the first intersection.  Also, the
  357. barycentric coordinates become somewhat meaningless, since more than one
  358. triangle can be found to overlap your test point.
  359.  
  360.  
  361. So, included in the new code (which is at the end of this issue) are two
  362. flavors of the crossings (segment) test, one which is fairly fast and
  363. readable, the other is "macmartinized" by me and runs a bit faster still.
  364. The barycentric polygon intersector is also included, and it works for all
  365. polygons (concave, self-intersecting, etc).
  366.  
  367. Random polygon testing:
  368.  
  369.              number (or range) of vertices
  370.             3        4        5        3-7        20        100
  371.  
  372. angle           57.21   69.81   86.44   82.98  303.33 1435.78
  373. barycentric        1.80    3.44    5.29    4.90   30.78  164.38
  374. crossings        2.17    2.76    3.32    3.24   12.51   60.13
  375. macmartin        1.86    2.35    2.87    2.64   11.04   52.27
  376.  
  377. inside %        6.00   12.00   13.78    9.33   26.89   41.22
  378.  
  379. Times are in microseconds per intersection test, on an HP 720 workstation.
  380. The "for" loop for repeating each test is counted into these test times
  381. (basically, I repeat each polygon/point test a number of times so that there's
  382. some reasonable amount of time to count) - not counting it just means the
  383. "crossings" tests is all that much faster.  Basically, the macmartinized
  384. crossing test wins on my machine; your mileage may vary.
  385.  
  386. Note that the barycentric intersector could be faster if the polygons are
  387. known to be convex, since then a quick out could be taken (i.e. the first hit
  388. ends testing).  I won't swear at all for the completeness of the barycentric
  389. code:  I suspect you can get into roundoff error problems for intersection
  390. points which lay exactly on some internal, invisible triangle testing edge
  391. (rare, but possible).
  392.  
  393. There are also a few hacks I haven't done on the barycentric intersector:  you
  394. could make the set of if statements into two big if statements, the variable
  395. "v0" is not all that useful, etc.  Also, if the polygon is a triangle the test
  396. could be made a special case and the "for" loop and other baggage could be
  397. eliminated (this I tried, and it resulted in about a 20% reduction in run
  398. time).  However, the code presented is at least readable, and the general
  399. behavior is the same:  this intersector gets more expensive than the crossing
  400. methods as the number of polygon vertices rises.
  401.  
  402. Looking over the branch flow analysis, I find that my test case polygons
  403. (randomly generated points) cause a lot of algebraic tests to occur when using
  404. the crossing algorithms (i.e.  the exact x intersection point has to be
  405. computed when a polygon edge goes from one quadrant to its diagonal opposite).
  406. For real polygons, this case is much less prevalent since vertices tend to be
  407. close to each other, so the polygons I've generated are more pathological than
  408. normally encountered.  For example, if the polygons generated are regular
  409. polygons inscribed in a circle of radius 1.0 and given a random rotation about
  410. the origin, the timings are
  411.  
  412. Regular polygon testing:
  413.  
  414.              number (or range) of vertices
  415.             3        4        5        3-7        20        100
  416.  
  417. angle           54.20   70.80   85.14   87.83  314.53 1522.00
  418. barycentric        2.11    3.98    5.48    5.64   27.63  144.27
  419. crossings        2.50    3.09    3.30    3.42    9.46   42.71
  420. macmartin        2.32    2.81    2.91    2.92    6.29   23.91
  421.  
  422. inside %       33.22   51.00   60.22   55.22   77.44   80.22
  423.  
  424. I'm not sure why there's any noticeable difference in the "angle summation"
  425. tests, since this test should be the same for any test point (I suspect it's
  426. just the "times()" routine's granularity).  The barycentric times change a
  427. little for better and worse, probably due to times() and to different polygon
  428. vs. point tests.  The two crossings test noticeably improve for large polygons
  429. (29% and 54% savings, respectively) - this is due to less "difficult" edges,
  430. and for the macmartin test there are now longer series of edges fully above or
  431. fully below the test point.  These tests show off the real improvement of the
  432. macmartin algorithm when testing many sided polygons.
  433.  
  434. The near-final word is that the crossings tests seems to be the most efficient
  435. overall, the angle test should generally be avoided like the plague, and the
  436. barycentric test is good for the special case of triangles (and you also
  437. get the vertex interpolant weights).  Just to add to the confusion, I talk
  438. about intersecting complex polygons in an article below; if you're using a lot
  439. of these then another intersection algorithm might be a better route to go.
  440.  
  441. -------------------------------------------------------------------------------
  442.  
  443. Polygon Intersection via Barycentric Coordinates, by Peter Shirley
  444.     (shirley@iuvax.cs.indiana.edu)
  445.  
  446. One good way to do ray triangle intersection is to use barycentric
  447. coordinates.  If the triangle has vertices p0, p1, p2, then any point in the
  448. plane containing the triangle can be represented as
  449.  
  450.     p = a*p0 + b*p2 + g*p3
  451.  
  452. and the constraint a + b + g = 1 will also hold.  If you want to interpolate
  453. colors or normals across the triangle, you can use a similar formula.  A
  454. little messing around will convince you that for points inside the triangle,
  455. a, b, and g are all positive, and outside the triangle, at least one is
  456. negative.  We can take this info, and get rid of a:
  457.  
  458.     p = p0 + b*(p1 - p0) + g*(p2-p0)
  459.  
  460. The point p is inside if b and g are positive, and (b+g) < 1.
  461.  
  462. A point on our ray can be represented as p' = o + t*v, where o is the point of
  463. origin on the ray, and v is the direction of propagation.  Plugging this into
  464. the triangle formula, we can find whether any point on the ray is also a point
  465. on the plane:
  466.  
  467.     o + t*v = p0 + b*(p1 - p0) + g*(p2-p0)
  468.  
  469. Note that this is really 3 linear equations (one for x, y and z).  You can
  470. rewrite this as a 3 by 3 linear system and solve for (t, b, g) and there is a
  471. hit if and only if:
  472.  
  473.     t > 0
  474.     b > 0
  475.     g > 0
  476.     b+g < 1
  477.  
  478. Alternatively, you can first find the equation of the plane containing the
  479. triangle and find t.  This will give you an explicit point h on the plane.
  480. You now have:
  481.  
  482.     h = p0 + b*(p1 - p0) + g*(p2-p0)
  483.  
  484. Now you have 3 equations and 2 unknowns.  Choosing any two of the equations is
  485. equivalent to projecting the problem onto one of the cartesian planes.  Watch
  486. out if the triangle is in one of the planes, because the projection might be a
  487. line segment.  You can check the surface normal to see if it is one of those
  488. cases (I just use the dimensions associated with the smallest magnitude
  489. components of the normal).
  490.  
  491. Note-- none of this is new.  GG II and the ray tracing book have other ways to
  492. do it.  Also, if your scenes are large, you are probably better off optimizing
  493. your spatial partitioning code (the part which takes longer for bigger
  494. scenes).
  495.  
  496. [For more info, see Didier Badouel's article in "Graphics Gems", p. 390-393]
  497.  
  498. -------------------------------------------------------------------------------
  499.  
  500. Many-Sided Polygon Intersection, by Eric Haines, Benjamin Zhu
  501.  
  502.  
  503. Jon Bennett (jb7m+@andrew.cmu.edu) writes:
  504.  
  505. >For a project i was working on i needed a very fast 2-D point-in-polygon
  506. >routine for a monte-carlo simulation, and came up with a variation of the
  507. >standard jordan curve theorem algorithm which ran approximately on the
  508. >order of O(n/20) when n > ~100 (for "mostly" round-ish polygons), instead
  509. >of O(n).  What I can't find, and no one here seems to know, is for large
  510. >polygons what is the best 2-D algorithm.  Everything they've heard of runs
  511. >in O(n) (I know that n and n/20 are the same "order" but the you know what
  512. >I mean). I'd like to know if there is something faster.
  513.  
  514.  
  515. You might check Preparata & Shamos's _Computational Geometry_ book, I believe
  516. they talk about faster algorithms.  Unfortunately, much of the stuff in CG
  517. tends to be things like O(K * log n), where K is something pretty large.
  518.  
  519. A simple way to test polygons with a large number of edges in O(n/K) time is
  520. to sort the edges into buckets, then find out which bucket the test point
  521. falls into.  Then you do the usual ray test, but only against those edges in
  522. the bucket.
  523.  
  524. For example, you could take a hundred edge polygon and find the Y extent in the
  525. polygon's plane.  Split this extent into, say, 20 buckets (or whatever number
  526. you like, more ==> more efficient, but also more memory and preprocessing).
  527. Now take each edge and note which buckets it's in.  As an added speed-up, sort
  528. the edges by their leftmost X values in that bucket.
  529.  
  530. Now testing can begin.  When you test a point, use the Y coordinate to find
  531. which bucket to check.  Now test each edge by the normal Jordan test against
  532. the point.  You only need to test edges to the left of the point (i.e. you are
  533. using a test ray in the -X direction), so as soon as you reach an edge whose
  534. leftmost X is greater than the point's X, you can stop testing.
  535.  
  536. The bucket sort immediately gets rid of a lot of edges, and the X test means
  537. often not having to test all edges in the bucket (it's optional icing).
  538.  
  539. If you're really into using memory and like preprocessing, you could use a
  540. grid structure to place buckets on the polygon.  Each grid cell is inside,
  541. outside, or indeterminate.  There are some tricks, but I leave these as an
  542. exercise for the reader...  Anyway, this method can be very fast, since much
  543. of the polygon's area is classified as in or out and much of the time no edges
  544. have to be tested at all.  The only limit is the grid cell resolution.  As the
  545. resolution rises the actual testing approaches constant time ( O(1)!).
  546. There's more preprocessing involved, but if you're making a lot of tests then
  547. this is worth doing.
  548.  
  549. ----
  550.  
  551. Benjamin Zhu (zhu@graphack.asd.sgi.com) comments on a similar question:
  552.  
  553. If your polygon is convex, (more generally, your polygon is star-shaped,) you
  554. can pre-process the polygon into triangle strips by connecting any point in
  555. the kernel of the polygon with each vertex of the polygon.  This takes you
  556. O(n) pre-processing time.  After that, for each point, you do a binary search
  557. in polar coordinates to locate the potential triangle where the point might
  558. lie.  Then you can determine if the point lies within this triangle, which is
  559. trivial.  See Preparata and Shamos' "Computational Geometry" for more details.
  560. Obviously, this algorithm will give you O(log(n)) time per point.  So you
  561. might consider this one if you need to solve the point-inclusion problem many
  562. times.
  563.  
  564. -------------------------------------------------------------------------------
  565.  
  566. Code for Point in Polygon Intersectors, by Eric Haines
  567.  
  568. Here it is, code for the angle test, the barycentric test, and the two
  569. crossings tests, with a main program to test their speeds.  You may need to
  570. use something other than srand()/drand48() for your random number generator,
  571. and the times() command for timing.  You'll also want to change the
  572. *_TEST_TIMES constants if you're using something slower than an HP 720
  573. workstation.  Feel free to complain about the slowness of any of the code -
  574. I'm always interested in new hacks.
  575.  
  576.  
  577. /* Point in polygon inside/outside tester.  Angle summation, barycentric
  578.  * coordinates, and ray along x-axis (crossings testing) compared.
  579.  *
  580.  * copyright 1992 by Eric Haines, 3D/Eye Inc, erich@eye.com
  581.  */
  582.  
  583. #include <math.h>
  584. #include <sys/types.h>
  585. #include <sys/param.h>
  586. #include <sys/times.h>
  587.  
  588. #define    X    0
  589. #define    Y    1
  590.  
  591. double    drand48() ;
  592.  
  593. double    AngleTimeTotal ;
  594. double    BarycentricTimeTotal ;
  595. double    CrossingsTimeTotal ;
  596. double    MacmartinTimeTotal ;
  597.  
  598. /* minimum & maximum number of polygon vertices to generate */
  599. #define    MIN_VERTS    3
  600. #define    MAX_VERTS    7
  601.  
  602. /* number of different polygons to try */
  603. #define    TEST_POLYGONS    30
  604.  
  605. /* number of different intersection points to try */
  606. #define    TEST_POINTS    30
  607.  
  608. /* number of times to try a single point vs. a polygon */
  609. /* this should be > 1 / ( HZ * approx. single test time in seconds ) */
  610. #define    ANGLE_TEST_TIMES    1000
  611. #define    BARYCENTRIC_TEST_TIMES    10000
  612. #define    CROSSINGS_TEST_TIMES    10000
  613. #define    MACMARTIN_TEST_TIMES    10000
  614.  
  615. main(argc,argv)
  616. int argc;  char *argv[];
  617. {
  618. int    i, j, numverts, inside_tot ;
  619. int    angle_flag, barycentric_flag, crossings_flag, macmartin_flag ;
  620. double    pgon[MAX_VERTS][2] ;
  621. double    point[2] ;
  622.  
  623.     srand( 12345 ) ;
  624.  
  625.     AngleTimeTotal = 0.0 ;
  626.     BarycentricTimeTotal = 0.0 ;
  627.     CrossingsTimeTotal = 0.0 ;
  628.     MacmartinTimeTotal = 0.0 ;
  629.     inside_tot = 0 ;
  630.  
  631.     for ( i = 0 ; i < TEST_POLYGONS ; i++ ) {
  632.  
  633. #ifdef CENTERED_SQUARE
  634.     /* for debugging purposes, test against a square */
  635.     numverts = 4 ;
  636.     pgon[0][X] = 0.2 ;
  637.     pgon[0][Y] = 0.2 ;
  638.     pgon[1][X] = 0.8 ;
  639.     pgon[1][Y] = 0.2 ;
  640.     pgon[2][X] = 0.8 ;
  641.     pgon[2][Y] = 0.8 ;
  642.     pgon[3][X] = 0.2 ;
  643.     pgon[3][Y] = 0.8 ;
  644. #else
  645.     /* make an arbitrary polygon fitting 0-1 range in x and y */
  646.     numverts = MIN_VERTS +
  647.         (int)(drand48() * (double)(MAX_VERTS-MIN_VERTS+1)) ;
  648.     for ( j = 0 ; j < numverts ; j++ ) {
  649.         pgon[j][X] = drand48() ;
  650.         pgon[j][Y] = drand48() ;
  651.     }
  652. #endif
  653.  
  654.     /* now try # of points against it */
  655.     for ( j = 0 ; j < TEST_POINTS ; j++ ) {
  656.         point[X] = drand48() ;
  657.         point[Y] = drand48() ;
  658.         angle_flag = angletest( pgon, numverts, point ) ;
  659.         barycentric_flag = barycentrictest( pgon, numverts, point ) ;
  660.         crossings_flag = crossingstest( pgon, numverts, point ) ;
  661.         macmartin_flag = macmartintest( pgon, numverts, point ) ;
  662.  
  663.         /* reality check */
  664.         if ( angle_flag != crossings_flag ) {
  665.         printf( "angle test says %s, crossings test says %s\n",
  666.             angle_flag ? "INSIDE" : "OUTSIDE",
  667.             crossings_flag ? "INSIDE" : "OUTSIDE" ) ;
  668.         printf( "point %g %g\n", (float)point[X], (float)point[Y] ) ;
  669.         printf( "polygon:\n" ) ;
  670.         for ( j = 0 ; j < numverts ; j++ ) {
  671.             printf( "  %g %g\n", (float)pgon[j][X], (float)pgon[j][Y]);
  672.         }
  673.         }
  674.         if ( barycentric_flag != macmartin_flag ) {
  675.         printf( "barycentric test says %s, macmartin test says %s\n",
  676.             barycentric_flag ? "INSIDE" : "OUTSIDE",
  677.             macmartin_flag ? "INSIDE" : "OUTSIDE" ) ;
  678.         printf( "point %g %g\n", (float)point[X], (float)point[Y] ) ;
  679.         printf( "polygon:\n" ) ;
  680.         for ( j = 0 ; j < numverts ; j++ ) {
  681.             printf( "  %g %g\n", (float)pgon[j][X], (float)pgon[j][Y]);
  682.         }
  683.         }
  684.  
  685.         inside_tot += crossings_flag ;
  686.     }
  687.     }
  688.     printf( "angle test time: %g microseconds per test\n",
  689.     (float)( AngleTimeTotal/(double)(TEST_POINTS*TEST_POLYGONS) ) ) ;
  690.     printf( "barycentric test time: %g microseconds per test\n",
  691.     (float)( BarycentricTimeTotal/(double)(TEST_POINTS*TEST_POLYGONS) ) ) ;
  692.     printf( "crossings test time: %g microseconds per test\n",
  693.     (float)( CrossingsTimeTotal/(double)(TEST_POINTS*TEST_POLYGONS) ) ) ;
  694.     printf( "macmartin crossings test time: %g microseconds per test\n",
  695.     (float)( MacmartinTimeTotal/(double)(TEST_POINTS*TEST_POLYGONS) ) ) ;
  696.  
  697.     printf( "%g %% of all points were inside polygons\n",
  698.     (float)inside_tot * 100.0 / (float)(TEST_POINTS*TEST_POLYGONS) ) ;
  699. }
  700.  
  701. /* sum angles of vtxN - point - vtxN+1, check if in PI to 3*PI range */
  702. int
  703. angletest( pgon, numverts, point )
  704. double    pgon[MAX_VERTS][2] ;
  705. int    numverts ;
  706. double    point[2] ;
  707. {
  708. int    i, j, inside_flag ;
  709. struct  tms     timebuf ;
  710. long    timestart ;
  711. long    timestop ;
  712. double    *vtx0, *vtx1, angle, len, vec0[2], vec1[2] ;
  713.  
  714.     /* do the test a bunch of times to get a useful time reading */
  715.     timestart = times( &timebuf ) ;
  716.     for ( i = 0 ; i < ANGLE_TEST_TIMES ; i++ ) {
  717.     /* sum the angles and see if answer mod 2*PI > PI */
  718.     vtx0 = pgon[numverts-1] ;
  719.     vec0[X] = vtx0[X] - point[X] ;
  720.     vec0[Y] = vtx0[Y] - point[Y] ;
  721.     if ( (len = hypot( vec0[X], vec0[Y] )) <= 0.0 ) {
  722.         /* point and vertex coincide */
  723.         return( 1 ) ;
  724.     }
  725.     vec0[X] /= len ;
  726.     vec0[Y] /= len ;
  727.  
  728.     angle = 0.0 ;
  729.     for ( j = 0 ; j < numverts ; j++ ) {
  730.         vtx1 = pgon[j] ;
  731.         vec1[X] = vtx1[X] - point[X] ;
  732.         vec1[Y] = vtx1[Y] - point[Y] ;
  733.         if ( (len = hypot( vec1[X], vec1[Y] )) <= 0.0 ) {
  734.         /* point and vertex coincide */
  735.         return( 1 ) ;
  736.         }
  737.         vec1[X] /= len ;
  738.         vec1[Y] /= len ;
  739.  
  740.         /* check if vec1 is to "left" or "right" of vec0 */
  741.         if ( vec0[X] * vec1[Y] - vec1[X] * vec0[Y] >= 0.0 ) {
  742.         /* add angle due to dot product of vectors */
  743.         angle += acos( vec0[X] * vec1[X] + vec0[Y] * vec1[Y] ) ;
  744.         } else {
  745.         /* subtract angle due to dot product of vectors */
  746.         angle -= acos( vec0[X] * vec1[X] + vec0[Y] * vec1[Y] ) ;
  747.         }
  748.  
  749.         /* get to next point */
  750.         vtx0 = vtx1 ;
  751.         vec0[X] = vec1[X] ;
  752.         vec0[Y] = vec1[Y] ;
  753.     }
  754.     /* test if between PI and 3*PI, 5*PI and 7*PI, etc */
  755.     /* if we care about is winding number > 0, then just:
  756.            inside_flag = fabs(angle) > M_PI ;
  757.      */
  758.     inside_flag = fmod( fabs(angle) + M_PI, 4.0*M_PI ) > 2.0*M_PI ;
  759.     }
  760.     timestop = times( &timebuf ) ;
  761.     /* time in microseconds */
  762.     AngleTimeTotal += 1000000.0 * (double)(timestop - timestart) /
  763.     (double)(HZ * ANGLE_TEST_TIMES) ;
  764.  
  765.     return (inside_flag) ;
  766. }
  767.  
  768. /* barycentric, a la Gems I, with a little efficiency tuning */
  769. int
  770. barycentrictest( pgon, numverts, point )
  771. double    pgon[MAX_VERTS][2] ;
  772. int    numverts ;
  773. double    point[2] ;
  774. {
  775. int    i, tris_hit, inside_flag, p1, p2 ;
  776. struct  tms     timebuf ;
  777. long    timestart ;
  778. long    timestop ;
  779. double    tx, ty, u0, u1, u2, v0, v1, alpha, beta, denom ;
  780.  
  781.     /* do the test a bunch of times to get a useful time reading */
  782.     timestart = times( &timebuf ) ;
  783.     for ( i = 0 ; i < BARYCENTRIC_TEST_TIMES ; i++ ) {
  784.  
  785.     tx = point[X] ;
  786.     ty = point[Y] ;
  787.  
  788.     tris_hit = 0 ;
  789.  
  790.     for ( p1 = 1, p2 = 2 ; p2 < numverts ; p1++, p2++ ) {
  791.  
  792.         if ( ( u1 = pgon[0][X] - pgon[p2][X] ) == 0.0 ) {
  793.  
  794.         /* zero area test - optional */
  795.         if ( ( u2 = pgon[p1][X] - pgon[p2][X] ) == 0.0 ) {
  796.             goto NextTri;
  797.         }
  798.  
  799.         /* Compute intersection point */
  800.         if ( ( ( beta = ( tx - pgon[p2][X] ) / u2 ) < 0.0 ) ||
  801.                ( beta > 1.0 ) ) {
  802.  
  803.             goto NextTri;
  804.         }
  805.  
  806.         if ( ( v1 = pgon[0][Y] - pgon[p2][Y] ) == 0.0 ) {
  807.             goto NextTri;
  808.         }
  809.  
  810.         if ( ( alpha = ( ty - pgon[p2][Y] - beta *
  811.             ( pgon[p1][Y] - pgon[p2][Y] ) ) / v1 ) < 0.0 ) {
  812.             goto NextTri;
  813.         }
  814.  
  815.         } else {
  816.         if ( ( denom = ( pgon[p1][Y] - pgon[p2][Y] ) * u1 -
  817.             ( u2 = pgon[p1][X] - pgon[p2][X] ) *
  818.             ( v1 = pgon[0][Y] - pgon[p2][Y] ) ) == 0.0 ) {
  819.  
  820.             goto NextTri;
  821.         }
  822.  
  823.         /* Compute intersection point & subtract 0 vertex */
  824.         u0 = tx - pgon[p2][X] ;
  825.         v0 = ty - pgon[p2][Y] ;
  826.  
  827.         if ( ( ( ( beta = ( v0 * u1 - u0 * v1 ) / denom ) ) < 0.0 ) ||
  828.                ( beta > 1.0 ) ) {
  829.  
  830.             goto NextTri;
  831.         }
  832.         if ( ( alpha = ( u0 - beta * u2 ) / u1 ) < 0.0 ) {
  833.             goto NextTri;
  834.         }
  835.         }
  836.  
  837.         /* check gamma */
  838.         if ( alpha + beta <= 1.0 ) {
  839.         /* survived */
  840.         tris_hit++ ;
  841.         }
  842.  
  843.         NextTri: ;
  844.     }
  845.     inside_flag = tris_hit & 0x1 ;
  846.     }
  847.     timestop = times( &timebuf ) ;
  848.     /* time in microseconds */
  849.     BarycentricTimeTotal += 1000000.0 * (double)(timestop - timestart) /
  850.     (double)(HZ * BARYCENTRIC_TEST_TIMES) ;
  851.  
  852.     return (inside_flag) ;
  853. }
  854.  
  855. /* shoot a test ray along +X axis - slower version, less messy */
  856. int
  857. crossingstest( pgon, numverts, point )
  858. double    pgon[MAX_VERTS][2] ;
  859. int    numverts ;
  860. double    point[2] ;
  861. {
  862. int    i, j, inside_flag, xflag0 ;
  863. struct  tms     timebuf ;
  864. long    timestart ;
  865. long    timestop ;
  866. double    *vtx0, *vtx1, dv0 ;
  867. int    crossings, yflag0, yflag1 ;
  868.  
  869.     /* do the test a bunch of times to get a useful time reading */
  870.     timestart = times( &timebuf ) ;
  871.     for ( i = 0 ; i < CROSSINGS_TEST_TIMES ; i++ ) {
  872.  
  873.     vtx0 = pgon[numverts-1] ;
  874.     /* get test bit for above/below Y axis */
  875.     yflag0 = ( dv0 = vtx0[Y] - point[Y] ) >= 0.0 ;
  876.  
  877.     crossings = 0 ;
  878.     for ( j = 0 ; j < numverts ; j++ ) {
  879.         /* cleverness:  bobble between filling endpoints of edges, so
  880.          * that the previous edge's shared endpoint is maintained.
  881.          */
  882.         if ( j & 0x1 ) {
  883.         vtx0 = pgon[j] ;
  884.         yflag0 = ( dv0 = vtx0[Y] - point[Y] ) >= 0.0 ;
  885.         } else {
  886.         vtx1 = pgon[j] ;
  887.         yflag1 = ( vtx1[Y] >= point[Y] ) ;
  888.         }
  889.  
  890.         /* check if points not both above/below X axis - can't hit ray */
  891.         if ( yflag0 != yflag1 ) {
  892.         /* check if points on same side of Y axis */
  893.         if ( ( xflag0 = ( vtx0[X] >= point[X] ) ) ==
  894.              ( vtx1[X] >= point[X] ) ) {
  895.  
  896.             if ( xflag0 ) crossings++ ;
  897.         } else {
  898.             /* compute intersection of pgon segment with X ray, note
  899.              * if > point's X.
  900.              */
  901.             crossings += (vtx0[X] -
  902.             dv0*( vtx1[X]-vtx0[X])/(vtx1[Y]-vtx0[Y])) >= point[X] ;
  903.         }
  904.         }
  905.     }
  906.     /* test if crossings is odd */
  907.     /* if we care about is winding number > 0, then just:
  908.            inside_flag = crossings > 0 ;
  909.      */
  910.     inside_flag = crossings & 0x01 ;
  911.     }
  912.     timestop = times( &timebuf ) ;
  913.     /* time in microseconds */
  914.     CrossingsTimeTotal += 1000000.0 * (double)(timestop - timestart) /
  915.     (double)(HZ * CROSSINGS_TEST_TIMES) ;
  916.  
  917.     return (inside_flag) ;
  918. }
  919.  
  920. /* shoot a test ray along +X axis - macmartinized by me, a bit messier */
  921. int
  922. macmartintest( pgon, numverts, point )
  923. double    pgon[MAX_VERTS][2] ;
  924. int    numverts ;
  925. double    point[2] ;
  926. {
  927. int    i, inside_flag, xflag0 ;
  928. struct  tms     timebuf ;
  929. long    timestart ;
  930. long    timestop ;
  931. double    *p, *stop ;
  932. int    crossings ;
  933. double    tx, ty, y ;
  934.  
  935.     /* do the test a bunch of times to get a useful time reading */
  936.     timestart = times( &timebuf ) ;
  937.     for ( i = 0 ; i < MACMARTIN_TEST_TIMES ; i++ ) {
  938.  
  939.     crossings = 0 ;
  940.  
  941.     tx = point[X] ;
  942.     ty = point[Y] ;
  943.     y = pgon[numverts-1][Y] ;
  944.  
  945.     p = (double *)pgon + 1 ;
  946.     if ( ( y >= ty ) != ( *p >= ty ) ) {
  947.         /* x crossing */
  948.         if ( ( xflag0 = ( pgon[numverts-1][X] >= tx ) ) ==
  949.          ( *(double *)pgon >= tx ) ) {
  950.  
  951.         if ( xflag0 ) crossings++ ;
  952.         } else {
  953.         /* compute intersection of pgon segment with X ray, note
  954.          * if > point's X.
  955.          */
  956.         crossings += ( pgon[numverts-1][X] -
  957.         (y-ty)*( *(double *)pgon - pgon[numverts-1][X])/(*p-y)) >= tx ;
  958.         }
  959.     }
  960.  
  961.     stop = pgon[numverts] ;
  962.  
  963.     for ( y=*p,p += 2 ; p < stop ; y=*p,p+=2) {
  964.  
  965.         if ( y >= ty ) {
  966.         while ( (p < stop) && (*p >= ty) ) p+=2 ;
  967.         if ( p >= stop ) goto Exit ;
  968.         /* y crosses */
  969.         if ( ( xflag0 = ( *(p-3) >= tx ) ) ==
  970.              ( *(p-1) >= tx ) ) {
  971.  
  972.             if ( xflag0 ) crossings++ ;
  973.         } else {
  974.             /* compute intersection of pgon segment with X ray, note
  975.              * if > point's X.
  976.              */
  977.             crossings += ( *(p-3) -
  978.             (*(p-2)-ty)*( *(p-1)-*(p-3))/(*p-*(p-2))) >= tx ;
  979.         }
  980.         } else {
  981.         while ( (p < stop) && (*p < ty)) p+=2 ;
  982.         if ( p >= stop ) goto Exit ;
  983.         /* y crosses */
  984.         if ( ( xflag0 = ( *(p-3) >= tx ) ) ==
  985.              ( *(p-1) >= tx ) ) {
  986.  
  987.             if ( xflag0 ) crossings++ ;
  988.         } else {
  989.             /* compute intersection of pgon segment with X ray, note
  990.              * if > point's X.
  991.              */
  992.             crossings += ( *(p-3) -
  993.             (*(p-2)-ty)*( *(p-1)-*(p-3))/(*p-*(p-2))) >= tx ;
  994.         }
  995.         }
  996.     }
  997.  
  998.     Exit:
  999.     /* test if crossings is odd */
  1000.     /* if we care about is winding number > 0, then just:
  1001.            inside_flag = crossings > 0 ;
  1002.      */
  1003.     inside_flag = crossings & 0x01 ;
  1004.     }
  1005.     timestop = times( &timebuf ) ;
  1006.     /* time in microseconds */
  1007.     MacmartinTimeTotal += 1000000.0 * (double)(timestop - timestart) /
  1008.     (double)(HZ * MACMARTIN_TEST_TIMES) ;
  1009.  
  1010.     return (inside_flag) ;
  1011. }
  1012.