多角形の重心 - trinoの走り書き

2次元平面状の多角形の頂点座標 r_t = (x_t,y_t) がある。この内部の重心はどうすれば求まるか？仕事で必要になった。
最初は適当な画像を用意してそこに多角形を描いて中を塗ってそのピクセルと位置を掛け算して足し合わせて最後にピクセル数で割れば一応出てくるなと思って実装始めたけど、なんだか色々遠回りしている気がしてもっと直接求める方法はないか考えました。

まずは簡単のために原点の周りを左回りに後戻りなく一周回る場合を考えよう。この場合、各頂点から原点に向かって線を引けば三角形にばらされる。それぞれの三角形に対して計算し、後で合成すればよい。
重心の合成は、それぞれの重心×質量(今回の場合は面積)を足し合わせ、最後に全質量(全面積)で割れば求まる。
$G = \frac{ \sum_t G_t S_t}{ \sum_t S_t }$
三角形の重心は3つの頂点の平均位置なので簡単に求まる。
$G_t = \frac{1}{3}( r_t + r_{t-1} )$
2点と原点を囲む三角形の面積は、外積が平行四辺形の面積になるのでその半分である。
$S_t = \frac{1}{2} r_t \times r_{t-1}$
但し、これは原点に対して左回りの場合である。右回りだと符号が逆になる。とりあえず左回りを前提に話を進める。
これにより、重心は
$G = \frac{1}{3} \frac{\sum_t ( r_t + r_{t-1} )( r_t \times r_{t-1} )}{\sum_t r_t \times r_{t-1}}$
ちなみに、右回りの場合は分母と分子で符号がキャンセルされるので結果として正しい値が出る。

さて、これで原点に対して一周している場合に解くことが出来た。あとは一般に拡張するわけだが、結論を先に書いてしまえば上記の式は実は一般に成り立っている。証明は少々面倒なので割愛するが、直感的には逆向きになったら外積の値が逆になり適切に引いてくからです。

以下、簡単なプログラム

bool CenterOfMass( const int* aiBoundaryX, const int* aiBoundaryY, int iBoundaryCount, double& gx, double& gy )
{
	int sum(0), sumX(0), sumY(0);
	int xm = aiBoundaryX[iBoundaryCount-1];
	int ym = aiBoundaryY[iBoundaryCount-1];
	for( int iBoundary = 0; iBoundary < iBoundaryCount; ++iBoundary ){
		const int xp = aiBoundaryX[iBoundary];
		const int yp = aiBoundaryY[iBoundary];

		const int s = xp * ym - yp * xm;
		const int x = (xm + xp) * s;
		const int y = (ym + yp) * s;

		sum += s;  //面積×２
		sumX += x; //Ｘ重心×６
		sumY += y; //Ｙ重心×６

		xm = xp;
		ym = yp;
	}
	gx = sumX / ( 3.0 * sum );
	gy = sumY / ( 3.0 * sum );

	return 0 < sum;
}

ちなみに戻り値は右回りだとtrue、左回りだとfalseが返ります。