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