Gdal库计算形心方法。

在Gdal库中计算形心的方法如下：

int OGRGeometry::Centroid( OGRPoint *poPoint ) const

其函数实现中，是调用的Geos库中的GEOSGetCentroid()方法，最终在Geos的bool

Geometry::getCentroid(Coordinate& ret) const函数中创建了CentroidArea对象，并将几何对象添加进去，获取到其形心。代码如下：

CentroidArea cent;

cent.add(this);

if ( ! cent.getCentroid(c) )

return false;

可见，关键就是这个类CentroidArea。

返回的形心坐标结果：

ret = Coordinate(cg3.x/3.0/areasum2, cg3.y/3.0/areasum2);

计算流程：

此类的方法add添加几何对象，最终会执行如下步骤（略去内环考虑）：

void

CentroidArea::addShell(const CoordinateSequence *pts)

{

bool isPositiveArea=!CGAlgorithms::isCCW(pts);

std::size_t const n=pts->getSize()-1;

for(std::size_t i=0; i<n; ++i)

{

addTriangle(basePt, pts->getAt(i), pts->getAt(i+1), isPositiveArea);

}

addLinearSegments(*pts);

}

函数中的basePt被初始化为环的第一个点。

如代码所述，依次将每一个点加进类中，调用私有方法addTriangle ，其代码如下：

CentroidArea::addTriangle(const Coordinate &p0, const Coordinate &p1,

const Coordinate &p2, bool isPositiveArea)

{

double sign=(isPositiveArea)?1.0:-1.0;

centroid3(p0,p1,p2,triangleCent3);

double area2res=area2(p0,p1,p2);

cg3.x+=sign*area2res*triangleCent3.x;

cg3.y+=sign*area2res*triangleCent3.y;

areasum2+=sign*area2res;

}

成员属性cg3和areasum2得到了值。

其中area2res和triangleCent3计算如下：

double CentroidArea::area2(const Coordinate &p1, const Coordinate &p2, const Coordinate &p3)

{

return (p2.x-p1.x)*(p3.y-p1.y)-(p3.x-p1.x)*(p2.y-p1.y);

}

计算的是面积的2倍。

void CentroidArea::centroid3(const Coordinate &p1, const Coordinate &p2,

const Coordinate &p3, Coordinate &c)

{

c.x=p1.x+p2.x+p3.x;

c.y=p1.y+p2.y+p3.y;

}

c的xy除以3的话就是p1,p2,p3的均值了。

还有其他分支，似是一些备选，这里略过了。

所以是相当遍历所有点，将环的0点和当前点及当前点的下一个点构成一个三角形，进行累计计算。对于每一步：

wiki上的多边形形心计算公式：

一个由N个顶点（x_i , y_i）确定的不自交闭多边形的中心能如下计算： ^[4]

记号（ x_N , y_N）与顶点（ x₀ , y₀）相同。多边形的面积为：

多边形的中心由下式给出：

以面积中心来求：

面积中心和质量中心非常类似，面积中心只取决于图形的几何形状。如果物体是均匀的，质量中心将位于面积中心。

对于两部分组成的图形，将有如下等式：

是特定部分的面积中心到所选参考系的距离。是特定部分的面积。

当一个复杂几何图形可以分成一些已知的简单几何图形时，先计算各部分的面积中心，然后通过下面一般的公式计算整个图形的面积中心：

这里从y-轴到中心的距离是，从x-轴到中心的距离是，中心的坐标是。

http://zh.wikipedia.org/wiki/%E5%87%A0%E4%BD%95%E4%B8%AD%E5%BF%83