Date: Mon, 31 Jul 95 10:30:56 CDT
From: rusin (Dave Rusin)
To: jrothman@carbon.cudenver.edu
Subject: Re: Need Algorithm - Is point (x,y) inside the polygon?

>In article <3vccok$5qr@watson.math.niu.edu> you wrote:
>: Now you just need to know that area can be computed from the vertices's
>: coordinates  Pi = ( x[i], y[i] ):
>: 	Area = | Sum (x[i] y[i+1] - x[i+1] y[i]) |/2
>
>I assume that it is a cyclic sum so for 3 points the 3th term is
>x3y1-x1y3. I haven't seen this formula before. For 3 points it is
>equivalent to the determinant formula I am familar with. 

Yes, it is supposed to be cyclic, and yes, it's similar to the determinantal
formula. Indeed, You can imagine computing area of the polygon by drawing
a bunch of triangles joining consecutive vertices to the origin (I'm picturing
sort of an ellipse floating in the first quadrant). As you pass along the far
edge of the polygon, the triangles you draw will have a union which will
cover the interior of the polygon; you'll also cover the wedge which lies 
outside the polygon, between it and the origin, but as you keep computing
you'll start subtracting away this extra area as you traverse the near
side of the polygon. (I could draw a very helpful picture, but this is ASCII)

>But for 5 points it doesn't always work:
>
>(2,0)  (2,2)  (0,2)  (1,1)  (0,0)
>
>has area 3 but the formula gives area = (4+4+2+0+0)/2 = 5
                                             ^
should be "-" (the absolute value is taken at the end, not for each summand)

>When does it work and can you give a reference for it?

One way to demonstrate its validity is by recognizing it as a line integral.
I'm sending you a file from my Web site -- a post I made the last time
such an issue arose. It's a handy tool for shapes whose edge is given by
vertices or, more generally, as a parametric curve. Let me know if there's
anything else you need.

Incidentally, a similar technique applies to three-dimensional shapes
(e.g. surface area of polyhedra).

dave
http://www.math.niu.edu/~rusin
==============================================================================

From: rusin@washington.math.niu.edu (Dave Rusin)
Newsgroups: sci.math,comp.lang.c++,comp.lang.c
Subject: Re: Algorithm to find center of mass.
Date: 5 Dec 1994 19:52:59 GMT

Keith Sibson (ks2@st-and.ac.uk) wrote:
> When generating the asteroid I have values for the polygon vertices
> coordinates relative to each other. What I need is to find the coordinate
> of the center of mass so that I can rotate the vertices around this point.
> Can anyone detail an algorithm that could calculate this point?

Did I miss any posts of the line integral method? Although the posted
methods are OK they're based on surface integrals e.g. mass=integral(1 dxdy)
But the likely method of presenting an asteroid in this setting makes
it useful to use Green's theorem:
         _                   _
       _/  P dx + Q dy  =  _/   (dQ/dx - dP/dy) dxdy

so we can get the surface integrals of  1,  x,  and  y  respectively by
integrating  x,  x^2/2,  and  xy  w.r.t. y along the curve describing the
exterior of the asteroid.

You'll need to know that the curve joining two points (x0 y0) and (x1, y1)
has a parameterization   (x,y) = ( (1-t) x0 + t x1,  (1-t) y0 + t y1 )
so that  dy = (y1 - y0) dt. Thus, the integral of   x dy along this
curve is  integral_[0,1]  ( (1-t) x0 + t x1 ) . (y1-y0) dt   which
comes out to  (x0+x1)/2 . (y1-y0). Compute this product over all consecutive
pairs of vertices (don't forget to rejoin the last to the first) and
you'll have the area of the polygon.

Likewise you can compute the integrals of  x^2/2  and  xy  on this curve;
If I haven't messed up, the first is  (x0^2 + x0 x1 + x1^2)/6 . (y1-y0)  
and the latter is  ( (y0 x1 + x0 y1) + 2(x1 y1 + x0 y0) )/6 . (y1 - y0).

So here's a procedure to find the enclosed area:
Clear counters  Area, Xcoord, Ycoord.
For each  i= index of the points x[i], y[i]:
	Let  x0=x[i], y0=y[i]
	Let  x1=x[i+1], y1=y[i+1], where  i+1=first index if  i=last
	Compute factors  a1=(x0+x1)/2,
		a2=(x0^2+x0x1+x1^2)/6, and
		a3=(x0y1+y0x1+2(x1y1+x0y0))/6, as well as
		b=y1-y0
	Adjust counters Area=Area+a1. b,
		Xcoord=Xcoord+a2. b,  and
		Ycoord=Ycoord+a3. b.
Then the center of mass is at  (Xcoord/Area,  Ycoord/Area).

If you need to do this computation very often (i.e. you need it fast each
time) it makes sense to keep various partial products handy; I'll
leave this to professional programmers.

dave

==============================================================================

From: Richard Mercer <richard@seuss.math.wright.edu>
Newsgroups: comp.soft-sys.math.mathematica
Subject: Re: finding irregular areas
Date: 11 May 1995 05:41:27 GMT

[deletia - djr]

Here is a formula for the area of a polygon with vertices 

{(xk,yk): k = 1,...,n}:

Area = 1/2 [(x1*y2 - x2*y1) + (x2*y3 - x3*y2) + ... + (xn*y1 - x1*yn)].

This formula appears in an Article by Gil Strang of MIT on p. 253 of the March  
1993 issue of The American Mathematical Monthly, with the note that it is  
"known, but not well known". There is also a very brief discussion of proofs  
and other references, including an article by Bart Braden of Northern Kentucky  
U., a known Mathematica enthusiast.

[deletia - djr]

Richard Mercer