在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个顶点(xi , yi)确定的不自交闭多边形的中心能如下计算: [4]
记号( xN , yN)与顶点( x0 , y0)相同。多边形的面积为:
多边形的中心由下式给出:
以面积中心来求:
面积中心和质量中心非常类似,面积中心只取决于图形的几何形状。如果物体是均匀的,质量中心将位于面积中心。
对于两部分组成的图形,将有如下等式:
是特定部分的面积中心到所选参考系的距离。是特定部分的面积。
当一个复杂几何图形可以分成一些已知的简单几何图形时,先计算各部分的面积中心,然后通过下面一般的公式计算整个图形的面积中心:
这里从y-轴到中心的距离是,从x-轴到中心的距离是,中心的坐标是。
http://zh.wikipedia.org/wiki/%E5%87%A0%E4%BD%95%E4%B8%AD%E5%BF%83