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 complicated concave figures

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

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.

Go to:

Mail to: sky@interFOOBARgalact.com
© Copyright Intergalactic Reality, 1995.

Last updated: 5/22/95