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 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