我会说following 是一个用于计算简单多边形质心的权威 C 实现,它由 Computational 一书的作者Joseph O'Rourke 编写C 中的几何图形。
/*
Written by Joseph O'Rourke
orourke@cs.smith.edu
October 27, 1995
Computes the centroid (center of gravity) of an arbitrary
simple polygon via a weighted sum of signed triangle areas,
weighted by the centroid of each triangle.
Reads x,y coordinates from stdin.
NB: Assumes points are entered in ccw order!
E.g., input for square:
0 0
10 0
10 10
0 10
This solves Exercise 12, p.47, of my text,
Computational Geometry in C. See the book for an explanation
of why this works. Follow links from
http://cs.smith.edu/~orourke/
*/
#include <stdio.h>
#define DIM 2 /* Dimension of points */
typedef int tPointi[DIM]; /* type integer point */
typedef double tPointd[DIM]; /* type double point */
#define PMAX 1000 /* Max # of pts in polygon */
typedef tPointi tPolygoni[PMAX];/* type integer polygon */
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 );
}
/*
Returns twice the signed area of the triangle determined by a,b,c,
positive if a,b,c are oriented ccw, and negative if cw.
*/
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]);
}
/*
Returns the cg in CG. Computes the weighted sum of
each triangle's area times its centroid. Twice area
and three times centroid is used to avoid division
until the last moment.
*/
void FindCG( int n, tPolygoni P, tPointd CG)
{
int i;
double A2, Areasum2 = 0; /* Partial area sum */
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;
}
/*
Returns three times the centroid. The factor of 3 is
left in to permit division to be avoided until later.
*/
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');
}
/*
Reads in the coordinates of the vertices of a polygon from stdin,
puts them into P, and returns n, the number of vertices.
The input is assumed to be pairs of whitespace-separated coordinates,
one pair per line. The number of points is not part of the input.
*/
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,简要说明为here:
课题2.02:如何计算多边形的质心?
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.
我没有研究过这个算法,也没有测试过,但乍一看,它与维基百科的算法略有不同。
Graphics Gems IV一书中的代码是here:
/*
* ANSI C code from the article
* "Centroid of a Polygon"
* by Gerard Bashein and Paul R. Detmer,
(gb@locke.hs.washington.edu, pdetmer@u.washington.edu)
* in "Graphics Gems IV", Academic Press, 1994
*/
/*********************************************************************
polyCentroid: Calculates the centroid (xCentroid, yCentroid) and area
of a polygon, given its vertices (x[0], y[0]) ... (x[n-1], y[n-1]). It
is assumed that the contour is closed, i.e., that the vertex following
(x[n-1], y[n-1]) is (x[0], y[0]). The algebraic sign of the area is
positive for counterclockwise ordering of vertices in x-y plane;
otherwise negative.
Returned values: 0 for normal execution; 1 if the polygon is
degenerate (number of vertices < 3); and 2 if area = 0 (and the
centroid is undefined).
**********************************************************************/
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 来获得精确的计算,这将花费执行时间开销,这个想法在 The Exact Computation Paradigm 中有所描述。
一个commercial implementation声称使用格林定理,不知是否使用了多精度数类型:
面积和质心是通过应用格林定理计算得出的
只有轮廓或多边形上的点
我认为它指的是 Wikipedia 算法,因为 Wikipedia 中的公式是格林定理的应用,正如 here 所解释的那样。