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: edp@mv.com (Eric Postpischil)
Subject: Re: Algorithm to find center of mass.
Newsgroups: sci.math,comp.lang.c++,comp.lang.c
Date: Thu, 8 Dec 1994 00:38:26 GMT
[post reformatted -- djr]
Here's a C program to compute center of mass. Somebody previously
showed how to integrate 1, x, and y over the area of the polygon by
evaluating line integrals around the perimeter. Kudos to them, who I
am unable to name at the moment; this news server appears to expire
articles quickly.
I was proceeding along the same lines, thinking of the problem in
terms of adding in the contribution of each trapezoid formed by
adjacent points on the polygon and their projections onto the x axis.
Using trapezoids yields the same formulae the earlier poster gave.
Once I had those, I changed it to triangles with vertices at the
origin and two adjacent vertices of the polygon. This yields the
slightly simpler formulae used below.
/* Find center of mass of polygon. Written by Eric Postpischil. */
#include
#include
#include
/* Define a point to be an ordered pair. */
typedef struct
{
double x, y;
} Point;
/* This is a safe allocation routine for points. Input is a pointer to an
array of points, a pointer to the current number of allocated points,
and the index (not number) of the last desired point.
If sufficiently many points are allocated, nothing is done. Otherwise,
enough memory is allocated to provide for the desired point, and then
some.
Initially, when no points are allocated, p must be NULL and current
must be zero. */
Point *spalloc(Point **p, int *current, int desired)
{
if (desired >= *current)
{
*current = desired < 8 ? 8 : 2 * desired;
*p = realloc(*p, *current * sizeof **p);
if (*p == NULL)
{
fprintf(stderr, "Out of memory\n");
exit(1);
}
}
}
/* Read points describing the vertices of a polygon and compute its center
of mass. */
void main(void)
{
Point *p = NULL;
double area, t, x, y;
int allocated = 0, i, n;
/* Read each pair of points. */
for (n = 0; 2 == scanf("%lg%lg", &x, &y); n++)
{
spalloc(&p, &allocated, n);
p[n].x = x;
p[n].y = y;
}
if (n == 0)
printf("Undefined\n");
else
{
/* Copy the 0-th point to the n-th point, to close the loop. */
spalloc(&p, &allocated, n);
p[n] = p[0];
/* Compute the center of mass. See comments below. */
area = x = y = 0;
for (i = 0; i < n; i++)
{
/* Define x0, x1, y0, and y1 to refer to the current and next x and y
values. These just make the following formulae a little less cumbersome. */
#define x0 (p[i].x)
#define x1 (p[i+1].x)
#define y0 (p[i].y)
#define y1 (p[i+1].y)
/* Consider a triangle with vertices at (0,0), (x0,y0), (x1,y1). The union of
each such triangle, as the points (x0,y0) and (x1,y1) are iterated through the
vertices of the polygon, are the entire polygon. This is clear if the origin
is inside the polygon, but, even if the polygon does not contain the origin,
some of these triangles will have negative area, and the sum will be correct
for the polygon.
Given this, we can compute the integral of a function over area of the polygon
by computing the integrals for each triangle and summing them. */
/* Sum integral of 1 over that triangle. */
area += (y1*x0-y0*x1)/2;
/* Sum integral of x over that triangle. */
x += (x0+x1)*(y1*x0-y0*x1)/6;
/* Sum integral of y over that triangle. */
y += (y1+y0)*(y1*x0-y0*x1)/6;
}
/* Now divide to get int x dA / int 1 dA and
int y dA / int 1 dA. */
x /= area;
y /= area;
printf("area = %lg, center = (%lg, %lg)\n", area, x, y);
}
}
-- edp (Eric Postpischil)
"Always mount a scratch monkey."
edp@mv.com
PGP key fingerprint: 8E AD 63 61 BA 0C 26 86 32 0A 7D 28 DB E7 6F 75