This page is in response to a message that was posted on the sci.math.num-analysis news group. The post asked if there was an algorithm that would calculate the center of mass of a polygonal figure consisting of n points (x, y).
I thought this was an interesting problem, so I spent an evening and morning coming up with an algorithm, and writing a small Mac program that would allow one to draw an arbitray, irregular, non self-intersecting, polygon, and then calculate its center of mass. The results of that program are displayed below for various polygons.
The solution to this problem is fairly easy. It involves breaking up the polygon into a set of triangles, and then calculating the area and center of mass for each triangle. The first point (p(0)) of the polygon was chosen as a vertex of every triangle, with the points p(i) and p(i + 1) chosen as the other two vertices, for i = 1 to n - 2.
The area of each triangle is simply 1/2 the vector cross product of the sides {p(0), p(i)} and {p(0), p(i + 1)}. Note that this is a signed quantity, as it must be to take into account possible concavity of the polygon. The center of each triangle is the intersection of its medians, where each median is a line between a vertex of the triangle and the midpoint of the opposing side.
To compute the total area of the polygon, the areas of each triangle (including those negative ones (opposite cross product sense)) are summed. To compute the center of mass of the entire figure, the center of mass of each triagle, weighted by its area (also including negatives), is summed for both the x and y coordinates.
Note that this approach will only work for non self-intersecting polygons: i.e. those polygons whose perimeter does not cross itself. For a self-intersecting perimeter, a much more complicated algorithm which computes intersections and changes cross product sense is necessary. This algorithm is left as an exercise for the reader ;-).
This looks pretty much like you would expect (which is good!). The center of mass is at the intersection of the red lines. Note that the positive y direction is down from the top of the window:
Again, this looks pretty much like you would expect (including the concave parts):
This calculation works fine for figures with holes ("lacunae"), as long as those holes are connected to the perimeter by narrow channels. This is similar to the path of integration which would be used to integrate a complex function which had singularities ("poles") within the holes:
Multiple figures can be computed as long as they are linked by narrow bridges, so that again a continuous path of integration is formed:
Here is the C++ code for performing this calculation. This code was compiled and run on a 68040-Mac with the Symantec C++ 7.0 compiler and project manager.
typedef struct point_type
{
long double x; // x coordinate of point
long double y; // y coordinate of point
} point_type;
class mainwin_type: public win_type
{
public:
// ...other fields
// calculation methods
virtual void recalc(void);
virtual long double get_area(point_type p1, point_type p2, point_type p3);
virtual point_type get_center(point_type p1, point_type p2, point_type p3);
int points; // number of points
point_type **point; // handle to point array
point_type center; // center of mass
long double area; // area of figure
int state; // drawing state
};
void mainwin_type::recalc(void)
{
area = 0; // initialize area to zero
center.x = 0; // initialize center to zero
center.y = 0; // initialize center to zero
// loop through all traingles in figure
for (int i = 2; i < points; i ++)
{
point_type p1 = (*point)[0]; // first vertex of triangle
point_type p2 = (*point)[i - 1]; // second vertex of triangle
point_type p3 = (*point)[i]; // third vertex of triangle
// compute area of triangle
long double a = get_area(p1, p2, p3);
area += a; // sum area
// only consider non-degenerate triangles
if (a != 0)
{
// compute center of triangle
point_type cntr = get_center(p1, p2, p3);
// weight center by area, and add to center coordinates
center.x += a * cntr.x;
center.y += a * cntr.y;
}
}
if (area != 0)
{
center.x /= area; // renormalize coords by total area
center.y /= area; // renormalize coords by total area
}
// area is actually a positive value
area = fabs(area);
}
long double mainwin_type::get_area(point_type p1, point_type p2, point_type p3)
{
long double v1x, v1y, v2x, v2y;
v1x = p2.x - p1.x; // compute side vectors
v1y = p2.y - p1.y; // compute side vectors
v2x = p3.x - p1.x; // compute side vectors
v2y = p3.y - p1.y; // compute side vectors
// return cross product of side vectors / 2
return 0.5 * (v1x * v2y - v1y * v2x);
}
point_type mainwin_type::get_center(point_type p1, point_type p2, point_type p3)
{
point_type cntr;
cntr.x = 0;
cntr.y = 0;
long double x1, y1, x2, y2, x3, y3, x4, y4;
x1 = p1.x; // get endpoints of first median
y1 = p1.y; // get endpoints of first median
x2 = (p2.x + p3.x) / 2.0; // get endpoints of first median
y2 = (p2.y + p3.y) / 2.0; // get endpoints of first median
x3 = p2.x; // get endpoints of second median (only need two)
y3 = p2.y; // get endpoints of second median
x4 = (p3.x + p1.x) / 2.0; // get endpoints of second median
y4 = (p3.y + p1.y) / 2.0; // get endpoints of second median
// see if either median is vertical (slope == infinity)
if (x1 == x2) // if so...
{
x1 = p3.x; // use third median (can't be two vertical medians)
y1 = p3.y; // use third median
x2 = (p1.x + p2.x) / 2.0; // use third median
y2 = (p1.y + p2.y) / 2.0; // use third median
}
else if (x3 == x4)
{
x3 = p3.x;
y3 = p3.y;
x4 = (p1.x + p2.x) / 2.0;
y4 = (p1.y + p2.y) / 2.0;
}
long double a1, a2, b1, b2;
a1 = (y2 - y1) / (x2 - x1); // compute slope of first median
b1 = y1 - a1 * x1; // compute intercept of first median
a2 = (y4 - y3) / (x4 - x3); // compute slope of second median
b2 = y3 - a2 * x3; // compute intercept of second median
// solve a1 * x + b1 = a2 * x + b2
cntr.x = (b2 - b1) / (a1 - a2); // solve for x coordinate of intersection
cntr.y = a1 * cntr.x + b1; // solve for y coordinate of intersection
return cntr; // return center as a point
}
I make no claim that this is the best strategy to use for this problem, or that this particular implementation is an efficient one. However, the algorithm seems to work for all cases, and is very fast for a large (~100) number of points.
Download the file: cntr.sea.hqx, it is about 88K.
This is a binhex of a self-extracting archive for the Macintosh. You will need to unbinhex it once you download it (or your web client may do this for you). Just double-click on the resulting application to unstuff the program. There is a "READ-ME" file in the folder; please read it.