求多边形形心坐标算法

[xb4270293]

已知多边形的各顶点坐标(X,Y)，如何计算形心坐标？

[opentuxedo]

Sigma(X)/n,Sigma(Y)/n

[wqs6]

http://course.bitsde.com/tt/normal/cllx/kjzt/clkc/fulu/tm02/tm02.htm

[shaolunyuan]

/* 求任意简单多边形polygon的重心
需要调用下面几个函数：
void AddPosPart(); 增加右边区域的面积
void AddNegPart(); 增加左边区域的面积
void AddRegion(); 增加区域面积
在使用该程序时，如果把xtr,ytr,wtr,xtl,ytl,wtl设成全局变量就可以使这些函数的形式得到化简,但要注意函数的声明和调用要做相应变化
*/
void AddPosPart(double x, double y, double w, double &xtr, double &ytr, double &wtr)
{
if (abs(wtr + w)<1e-10 ) return; // detect zero regions
xtr = ( wtr*xtr + w*x ) / ( wtr + w );
ytr = ( wtr*ytr + w*y ) / ( wtr + w );
wtr = w + wtr;
return;
}
void AddNegPart(double x, ouble y, double w, double &xtl, double &ytl, double &wtl)
{
if ( abs(wtl + w)<1e-10 )
return; // detect zero regions

xtl = ( wtl*xtl + w*x ) / ( wtl + w );
ytl = ( wtl*ytl + w*y ) / ( wtl + w );
wtl = w + wtl;
return;
}

void AddRegion ( double x1, double y1, double x2, double y2, double &xtr, double &ytr,
double &wtr, double &xtl, double &ytl, double &wtl )
{
if ( abs (x1 - x2)< 1e-10 )
return;

if ( x2 > x1 )
{
AddPosPart ((x2+x1)/2, y1/2, (x2-x1) * y1,xtr,ytr,wtr); /* rectangle 全局
变量变化处 */
AddPosPart ((x1+x2+x2)/3, (y1+y1+y2)/3, (x2-x1)*(y2-y1)/2,xtr,ytr,wtr);
// triangle 全局变量变化处
}
else
{
AddNegPart ((x2+x1)/2, y1/2, (x2-x1) * y1,xtl,ytl,wtl);
// rectangle 全局变量变化处
AddNegPart ((x1+x2+x2)/3, (y1+y1+y2)/3, (x2-x1)*(y2-y1)/2,xtl,ytl,wtl);
// triangle 全局变量变化处
}
}

POINT cg_simple(int vcount,POINT polygon[])
{
double xtr,ytr,wtr,xtl,ytl,wtl;
//注意： 如果把xtr,ytr,wtr,xtl,ytl,wtl改成全局变量后这里要删去
POINT p1,p2,tp;
xtr = ytr = wtr = 0.0;
xtl = ytl = wtl = 0.0;
for(int i=0;i<vcount;i++)
{
p1=polygon[i];
p2=polygon[(i+1)%vcount];
AddRegion(p1.x,p1.y,p2.x,p2.y,xtr,ytr,wtr,xtl,ytl,wtl); //全局变量变化处
}
tp.x = (wtr*xtr + wtl*xtl) / (wtr + wtl);
tp.y = (wtr*ytr + wtl*ytl) / (wtr + wtl);
return tp;
}

[shaolunyuan]

// 求凸多边形的重心,要求输入多边形按逆时针排序
POINT gravitycenter(int vcount,POINT polygon[])
{
POINT tp;
double x,y,s,x0,y0,cs,k;
x=0;y=0;s=0;
for(int i=1;i<vcount-1;i++)
{
x0=(polygon[0].x+polygon[i].x+polygon[i+1].x)/3;
y0=(polygon[0].y+polygon[i].y+polygon[i+1].y)/3; //求当前三角形的重心
cs=multiply(polygon[i],polygon[i+1],polygon[0])/2;
//三角形面积可以直接利用该公式求解
if(abs(s)<1e-20)
{
x=x0;y=y0;s+=cs;continue;
}
k=cs/s; //求面积比例
x=(x+k*x0)/(1+k);
y=(y+k*y0)/(1+k);
s += cs;
}
tp.x=x;
tp.y=y;
return tp;
}