【问题标题】:Are there any numerically stable versions of the centroid finding algorithm for polygons?是否有任何数值稳定的多边形质心查找算法版本?
【发布时间】:2016-12-23 05:52:53
【问题描述】:

假设我有一个几乎退化的二维多边形,例如:

[[40.802,9.289],[40.875,9.394],[40.910000000000004,9.445],[40.911,9.446],[40.802,9.289]]

供参考,如下所示:

如果我使用on Wikipedia所示的标准质心算法,例如这个python代码:

pts = [[40.802,9.289],[40.875,9.394],[40.910000000000004,9.445], [40.911,9.446],[40.802,9.289]]
a = 0.0
c = [0.0, 0.0]
for i in range(0,4):
    k = pts[i][0] * pts[i + 1][1] - pts[i + 1][0] * pts[i][1]
    a += k
    c = [c[0] + k * (pts[i][0] + pts[i + 1][0]), c[1] + k * (pts[i][1] + pts[i + 1][1])]
c = [c[0] / (3 * a), c[1] / (3 * a)]

我收到c = [-10133071.666666666, -14636692.583333334]。在a == 0.0 的其他情况下,我也可能会被零除。

我最理想的情况是,在最坏的情况下,质心等于顶点之一或多边形内的某个位置,并且不应使用任意公差来避免这种情况。是否有一些巧妙的方法可以重写方程以使其在数值上更加稳定?

【问题讨论】:

  • 您可以尝试的一件简单的事情是将整个多边形移动到原点(即减去角的平均值)。这至少会给你一些更高的浮点精度。
  • 你的程序错了,c不能乘以k。 (但修复后,数值不稳定仍然存在。)
  • @NicoSchertler:我认为这没有帮助,问题是面积接近于零。
  • 另一种常见的方法是避免分割并使用齐次坐标。在这样的设置中,您会得到一个接近空向量的结果,它不代表一个点。原点的齐次坐标为(0,0,1)或其任意倍数。如果您可以在没有除法的情况下继续,这会有所帮助,但它只是为您提供了“未定义此质心”的不同表示,而不是您预期的多边形内的点,所以我只是将其作为评论发布。跨度>
  • @YvesDaoust 同意这不是一个有保证的解决方案(这就是为什么它只是一个评论)。但是,原点附近增加的精度可能足以表示接近零的区域。

标签: algorithm geometry polygon numeric


【解决方案1】:

我会说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 允许您使用精确的多精度数字类型而不是 doublefloat 来获得精确的计算,这将花费执行时间开销,这个想法在 The Exact Computation Paradigm 中有所描述。

一个commercial implementation声称使用格林定理,不知是否使用了多精度数类型:

面积和质心是通过应用格林定理计算得出的 只有轮廓或多边形上的点

我认为它指的是 Wikipedia 算法,因为 Wikipedia 中的公式是格林定理的应用,正如 here 所解释的那样。

【讨论】:

    【解决方案2】:

    当面积为零(或非常接近零,如果您无法进行精确算术)时,最好的选择可能是获取点集的周界质心。

    周长质心由多边形每条边的中点的加权和(权重是相应边的长度)与多边形周长的比率给出。

    在这种情况下,可以使用精确的算术计算质心。 红点是周界质心,绿点是真正的质心

    我使用 sage 精确计算质心https://cloud.sagemath.com/projects/f3149cab-2b4b-494a-b795-06d62ae133dd/files/2016-08-17-102024.sagews

    人们一直在寻找一种将这些观点相互关联的方法——https://math.stackexchange.com/questions/1173903/centroids-of-a-polygon

    【讨论】:

      【解决方案3】:

      我认为这个公式对于几乎退化的 2D 多边形来说不容易变得更稳定。问题在于面积 (A) 的计算依赖于减去梯形形状(请参阅Paul Bourke)。对于非常小的区域,您不可避免地会遇到数值精度。

      我看到了两种可能的解决方案:

      1.) 您可以检查该区域,如果它低于阈值,则假设多边形已退化,只需取最小和最大 x 和 y 值的平均值(线的中间)

      2.). 使用精度更高的浮点运算,比如mpmath

      顺便说一句。你的代码有错误。应该是:

      c = [c[0] + k * (pts[i][0] + pts[i + 1][0]), c[1] + k * (pts[i][1] + pts[i + 1][1])]
      

      但这并没有什么不同。

      【讨论】:

        猜你喜欢
        • 2013-11-22
        • 1970-01-01
        • 1970-01-01
        • 1970-01-01
        • 1970-01-01
        • 1970-01-01
        • 1970-01-01
        • 2011-05-24
        • 2011-01-22
        相关资源
        最近更新 更多