[comp.graphics] Point inside of polygon revisited... A NEW APPROACH

klaiber@udel.EDU (Alexander Klaiber) (06/17/87)

Any other suggestions?  Well, there IS another way, which requires a little
preprocessing, but otherwise is quite efficient: (concave polygons only)
for n edges, we need 3*n multiplications and 4*n additions/subtractions/tests

Preprocessing is affordable, given it's only performed once whereas
intersection tests occur kinda frequently in raytracing. (It's rather
simple, anyway)

Now for those that aren't familiar with the normal-form of a plane: A plane
in 3-space can either be represented by one point on the plane and two
non-collinear vectors spanning the plane.
Example: r=r0+lambda*u+mu*v, 	where r0 is the point on the plane, lambda
and mu are parameters and u and v are the vectors spanning the plans.

Very often another representation, the normal form, is desirable. Here a
plane is defined by its normal vector and its distance from the origin of
the coordinate-system.
Example: r*n - d = 0, 	where n is the normal vector, '*' is the dot
product and d is the distance of the plane from the origin. Note that d may
be computed by simply picking a point on the plane, say r0, and computing
d=r0*n. Note that d will be positive if the origin lies "on the other side"
of the plane withrespect to the direction of the normal vector. I.e. the
direction of the normal vector defines an "up" and a "down" side of the
plane. This will be important.
To see where a point given by vector r lies with respect to the plane, all
you do is plug it in above formula, i.e. compute r*n-d. There are tree
cases for the result:
	>0	the point is on the side of the plane the normal points to
	=0	the point is on the plane
	<0	the point is on the "other" side of the plane

	       the plane normal
		    ^
		    |   (>0 here)
		    |
                    |
        ----------------------------- the plane (=0 on the plane)

			(<0 here)


NOW: assume the Polygon is specified by its vertices P1...Pn. Assume these all
lie in the same pane. Also, assume we have a point PI which is known to be
inside the polygon (use 1/n*(P1+P2+...+Pn), for example).

Now we'll represent the polygon by 
(1) the plane it lies in 
(2) planes, 1 per edge, that "stand up" perpendicular to the polygon-plane 
    and that are erected over "their" respective edges.

Think of the polygon lying flat on the table, then the (projection onto the
table of the) planes in (2) (which stand up vertically from the table)
represent the "cuts" you'd have to make to cut out the polygon from a sheet
of paper.

Now, representing each of these "walls" in the normal-form (normal ni and
distance di), all we have to do is select the direction of the ni's such
that our sample point PI is on the side the normal ni points to, that is,
	rPI * ni - di >0, 	where rPI is the space vector for PI
this is simple: Just plug PI into the formula and if the result is<0, jus
reverse the signs of ni and di.

FINALLY: Given a point P (vector r) that lies on the polygon-plane. To
detect whether it is also inside the polygon, do the following:
      For every edge i (now given by ni,di):
	 compute temp := r * ni - di
	 if temp<0, then exit with "point outside polygon"
      end
      "if we reach end of FOR: point is inside"

ANALYSIS: 
* time: For two 3-vectors x=(x1 x2 x3) and y=(y1 y2 y3), the dot product is
	x1*y1 + x2*y2 + x3*y3. This requires 3 multiplications and 2
	additions. Moreover, we have one subtraction (-di) and one test
	(<0) per edge. For n edges, this gives a worst case of
		3n mult's and 4n add's
	the worst case occurs, if the point indeed is inside the polygon,
	which is of course not always the case.
* space: After preprocessing, we can of course throw away the points P1..Pn
        and the inner point, PI. What we keep is:
	(1) the polygon plane in normal form: normal n and dist d
	(2) for each edge: ni and di
	{(1) is needed to find the intersection of the ray with the poly
	 plane in the first place}
	Assuming 3-vectors, a polygon with n edges requires
		4(n+1)
	real numbers.

The beauty of this method is that it's fast enough, works entirely in
world-coordinates and does not have any of the problems of the other
proposal.


---------------------------------------------------------------------------
Standard disclaimer etc etc etc.

dleigh@hplabsz.HPL.HP.COM (Darren Leigh) (06/19/87)

In article <284@louie.udel.EDU>, klaiber@udel.EDU (Alexander Klaiber) writes:
> 
> Any other suggestions?  Well, there IS another way, which requires a little
> preprocessing, but otherwise is quite efficient: (concave polygons only)
                                                   ^^^^^^^^^^^^^^^^^^^^^^^

I really like your solution.  It's neat and tidy.  I guess you mean
convex polygons, not concave.  This algorithm won't work for concave
polygons because one part of the polygon might be on one side of a
particular bounding plane while another part would be on the other side.

The other solutions suggested before also work well and efficiently
for convex polygons.  Using the original method of counting the number
of edge intersections between a known inside point and a known outside
point can be simplified, too.  We no longer have to count the number of
intersections:  if there is at least one, then the point is inside 
the polygon.

Darren Leigh
dlleigh@media-lab.mit.edu
or
dleigh@hplabs.hp.com

klaiber@udel.UUCP (06/20/87)

In article <538@hplabsz.HPL.HP.COM> dleigh@hplabsz.HPL.HP.COM (Darren Leigh) writes:
>[...] I guess you mean
>convex polygons, not concave.  This algorithm won't work for concave
>polygons because one part of the polygon might be on one side of a
>particular bounding plane while another part would be on the other side.

ooooops! Yup, I goofed on that one. Of course I meant to say "convex", but
it must have been late or maybe I had other things on my mind...
Alright everybody: "convex" is it -- or, to put it explicitly, polygons
with the property:
	for every two points P1,P2 inside the polygon, the line P1-P2 is
	also entirely inside the polygon.
(always had trouble with convex vs. concave in highschool.....)

Thanks for pointing out,


      Alexander Klaiber
      klaiber@dewey.udel.edu

ph@pixar.UUCP (Paul Heckbert) (06/20/87)

In article <284@louie.udel.EDU>, klaiber@udel.EDU (Alexander Klaiber) writes:
> ... there IS another way, which requires a little preprocessing,
> but otherwise is quite efficient: (concave polygons only) for n edges,
> we need 3*n multiplications and 4*n additions/subtractions/tests ...

Actually the method you describe works for convex polygons only, not concave.

There's a big difference between point-in-polygon testing for convex polyons
and for concave polygons.  It's easy to write robust software for the former
but the latter is often prone to numerical problems and special cases
(hence all the discussion here).

Here's a quick summary of some of the 2-d point-in-polygon methods I know:

POLYGON TYPE	TIME	METHOD

convex		O(n)	dot product with line equation for each edge
			(represent polygon as an intersection of half-planes)

star		O(logn)	binary search on angle from "center" of polygon, test
			inclusion in a single triangle
			(decompose polygon into triangles radiating from center)
			this method requires preprocessing.

concave		O(n)	Jordan Curve Theorem: intersect a ray out of the
			polygon with the edges, count intersections;
			odd means inside, even means outside.

concave		O(n)	winding number method
			Trace a ray out of the polygon and compute the total
			number of intersections which cross the ray going from
			right to left minus the number that cross from left to
			right.  The winding number=0 when outside the polygon.

The two concave methods work for non-simple polygons (polygons which intersect
themselves) too, although they give different results, as described in the
"PostScript Language Referance Manual".  These methods are described more fully
in Preparata & Shamos's book "Computational Geometry".

For 3d point-in-polygon testing it's cheaper to project the point and polygon
to 2d and test inclusion in the projected polygon than to test against planes,
as you did.  Project onto either the x-y, y-z, or z-x plane, use whichever one
has the largest projection.

Paul Heckbert
Pixar				415-499-3600
P.O. Box 13719			UUCP: {sun,ucbvax}!pixar!ph
San Rafael, CA 94913		ARPA: ph%pixar.uucp@ucbv Al Al cter unn

airey@unc.cs.unc.edu (John Airey) (06/26/87)

Paul Heckbert says use two-d solutions for three-d polygons by
choosing the largest projection onto the xy,yz,xz planes.
The largest projection may be determined by finding the component of the
normal to the poly with largest abs val and projecting onto
the "other" plane. i.e. if z normal component largest, project to xy.
sorry, I had to put in my 2 cents.