Polygon Center of Mass

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 ;-).

### Center of mass for a quadrilateral

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:

### Center of mass for an irregular figure

Again, this looks pretty much like you would expect (including the concave parts):

### Center of mass for irregular figure with lacunae

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:

### Center of mass for multiple linked irregular figures

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.

### First, a struct for a point type is defined:

```typedef struct point_type
{
long double  x;	// x coordinate of point
long double  y;	// y coordinate of point
} point_type;
```

### Here is part of the window class which contains the point set and the final area and center of mass:

```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
};
```

### Here is the method that does the actual calculations:

```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);
}
```

### Here is the method that computes the area of a triangle, given three points:

```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);
}

```

### Here is the method that computes the center of a triangle, given three points:

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