我认为以下是计算简单多边形重心的权威C语言实现，它是由《C语言计算几何》一书作者Joseph O'Rourke编写的。
#include <stdio.h>
#define DIM 2
typedef int tPointi[DIM];
typedef double tPointd[DIM];
#define PMAX 1000
typedef tPointi tPolygoni[PMAX];
int Area2( tPointi a, tPointi b, tPointi c );
void FindCG( int n, tPolygoni P, tPointd CG );
int ReadPoints( tPolygoni P );
void Centroid3( tPointi p1, tPointi p2, tPointi p3, tPointi c );
void PrintPoint( tPointd p );
int main()
{
int n;
tPolygoni P;
tPointd CG;
n = ReadPoints( P );
FindCG( n, P ,CG);
printf("The cg is ");
PrintPoint( CG );
}
int Area2( tPointi a, tPointi b, tPointi c )
{
return
(b[0] - a[0]) * (c[1] - a[1]) -
(c[0] - a[0]) * (b[1] - a[1]);
}
void FindCG( int n, tPolygoni P, tPointd CG)
{
int i;
double A2, Areasum2 = 0;
tPointi Cent3;
CG[0] = 0;
CG[1] = 0;
for (i = 1; i < n-1; i++) {
Centroid3( P[0], P[i], P[i+1], Cent3 );
A2 = Area2( P[0], P[i], P[i+1]);
CG[0] += A2 * Cent3[0];
CG[1] += A2 * Cent3[1];
Areasum2 += A2;
}
CG[0] /= 3 * Areasum2;
CG[1] /= 3 * Areasum2;
return;
}
void Centroid3( tPointi p1, tPointi p2, tPointi p3, tPointi c )
{
c[0] = p1[0] + p2[0] + p3[0];
c[1] = p1[1] + p2[1] + p3[1];
return;
}
void PrintPoint( tPointd p )
{
int i;
putchar('(');
for ( i=0; i<DIM; i++) {
printf("%f",p[i]);
if (i != DIM - 1) putchar(',');
}
putchar(')');
putchar('\n');
}
int ReadPoints( tPolygoni P )
{
int n = 0;
printf("Polygon:\n");
printf(" i x y\n");
while ( (n < PMAX) &&
(scanf("%d %d",&P[n][0],&P[n][1]) != EOF) ) {
printf("%3d%4d%4d\n", n, P[n][0], P[n][1]);
++n;
}
if (n < PMAX)
printf("n = %3d vertices read\n",n);
else printf("Error in ReadPoints:\too many points; max is %d\n",
PMAX);
putchar('\n');
return n;
}
该代码解决了书的第一版第47页上的练习12，这里有一个简短的解释
在这里：
Subject 2.02: How can the centroid of a polygon be computed?
The centroid (a.k.a. the center of mass, or center of gravity)
of a polygon can be computed as the weighted sum of the centroids
of a partition of the polygon into triangles. The centroid of a
triangle is simply the average of its three vertices, i.e., it
has coordinates (x1 + x2 + x3)/3 and (y1 + y2 + y3)/3. This
suggests first triangulating the polygon, then forming a sum
of the centroids of each triangle, weighted by the area of
each triangle, the whole sum normalized by the total polygon area.
This indeed works, but there is a simpler method: the triangulation
need not be a partition, but rather can use positively and
negatively oriented triangles (with positive and negative areas),
as is used when computing the area of a polygon. This leads to
a very simple algorithm for computing the centroid, based on a
sum of triangle centroids weighted with their signed area.
The triangles can be taken to be those formed by any fixed point,
e.g., the vertex v0 of the polygon, and the two endpoints of
consecutive edges of the polygon: (v1,v2), (v2,v3), etc. The area
of a triangle with vertices a, b, c is half of this expression:
(b[X] - a[X]) * (c[Y] - a[Y]) -
(c[X] - a[X]) * (b[Y] - a[Y]);
Code available at ftp://cs.smith.edu/pub/code/centroid.c (3K).
Reference: [Gems IV] pp.3-6; also includes code.
我没有学习过这个算法，也没有测试过它，但是乍一看，它似乎与维基百科上的算法略有不同。
来自书籍《图形宝石 IV》的代码在这里：
int polyCentroid(double x[], double y[], int n,
double *xCentroid, double *yCentroid, double *area)
{
register int i, j;
double ai, atmp = 0, xtmp = 0, ytmp = 0;
if (n < 3) return 1;
for (i = n-1, j = 0; j < n; i = j, j++)
{
ai = x[i] * y[j] - x[j] * y[i];
atmp += ai;
xtmp += (x[j] + x[i]) * ai;
ytmp += (y[j] + y[i]) * ai;
}
*area = atmp / 2;
if (atmp != 0)
{
*xCentroid = xtmp / (3 * atmp);
*yCentroid = ytmp / (3 * atmp);
return 0;
}
return 2;
}
CGAL允许您使用精确的多精度数字类型，而不是double或float，以获得精确的计算结果。这将会增加执行时间开销，其思想在精确计算范例中有所描述。
一种商业实现声称使用Green定理，但我不知道它是否使用了多精度数字类型：
通过仅使用轮廓或多边形上的点应用Green定理来计算面积和质心
我认为它指的是维基百科算法，因为维基百科中的公式是Green定理的一个应用，如此处所解释。
c不应该乘以k。（但修复后，数值不稳定性仍然存在。） - user1196549