python - 多边形面积计算不一致

Question

我正在使用 voronoi_finite_polygons_2d(vor, radius=None) 函数 found elsewhere on StackOverflow .

我想修改它以显示 centroid of each voronoi cell .调试为什么某些质心出现严重错误(请参阅绿色箭头，指出质心在杂草中偏离)。我发现的第一个错误:一些计算没有按照正确的全顺时针或全逆时针顺序处理顶点。

不确定为什么有些点没有正确排序，但在我调查之前，我发现了另一个异常情况。

如果我顺时针或逆时针走，我应该得到相同的区域(符号相反)。在简单的例子中，我做到了。但是在我制作的随机多边形中，我得到/稍微/不同的结果。

import numpy as np
import matplotlib.pyplot as plt
from scipy.spatial import Voronoi
import random
import math

def measure_polygon(vertices):
    xs = vertices[:,0]
    ys = vertices[:,1]
    xs = np.append(xs,xs[0])
    ys = np.append(ys,ys[0])

    #https://en.wikipedia.org/wiki/Centroid#Centroid_of_polygon
    area = sum(xs[i]*(ys[i+1]-ys[i-1]) for i in range(0, len(xs)-1))/2.0
    centroid_x =  sum((xs[i]+xs[i+1])*(xs[i]*ys[i+1] - xs[i+1]*ys[i]) for i in range(0, len(xs)-1))/(6.0*area)
    centroid_y =  sum((ys[i]+ys[i+1])*(xs[i]*ys[i+1] - xs[i+1]*ys[i]) for i in range(0, len(xs)-1))/(6.0*area)

    return (area, (centroid_x, centroid_y))

第一个示例按预期工作 - 相同的区域和质心，无论处理顺序如何(顺时针或逆时针)。

d = [[0.0 ,  0.0], [1.0,3.0],[  5.0,3.0],[   4.0 ,  0.0] ] 
print len(d)

defects = [] 
defects.append([d[0], d[1], d[2], d[3]]) 
defects.append([d[3], d[2], d[1], d[0]])

for v in defects:
    print measure_polygon(np.array(v))

简单的平行四边形输出:

4 
(-12.0, (2.5, 1.5)) 
(12.0, (2.5, 1.5))

但现在看看这个 4 边形(几乎是三角形)

#original list of vertices
d = [[-148.35290745 ,  -1.95467472], [-124.93580616 ,  -2.09420039],[  -0.58281373,    1.32530292],[   8.77020932 ,  22.79390931] ]
print len(d)

defects = []
#cw
defects.append([d[0], d[2], d[3], d[1]])
#ccw
defects.append([d[1], d[3], d[2], d[0]])

for v in defects:
    print measure_polygon(np.array(v))

给我奇怪的输出:

4
(1280.4882517358433, (-36.609159411740798, 7.5961622623413145))
(-1278.8546083623708, (-36.655924939495335, 7.6058658049196115))

地区不同。如果区域不同，则质心也会不同。面积差异(1280 与 1278)如此之大，我怀疑这是一个 float 舍入问题。但除此之外，我已经想不出为什么这不起作用。

===============================

我发现了这个错误....我的列表理解/索引黑客来启用 y-1 和 y+1 符号被破坏了(以一种半途而废的险恶方式)。正确的套路如下:

def measure_polygon(vertices):
    xs = vertices[:,0]
    ys = vertices[:,1]

    #the first and last elements are for +1 -1 to work at end of range
    xs = vertices[-1:,0]
    xs = np.append(xs,vertices[:,0])
    xs = np.append(xs,vertices[:1,0])

    ys = vertices[-1:,1]
    ys = np.append(ys,vertices[:,1])
    ys = np.append(ys,vertices[:1,1])

    #for i in range(1, len(xs)-1):
    #    print ("digesting x, y+1, y-1 points: {0}/{1}/{2}".format(xs[i], ys[i+1], ys[i-1]))

    #https://en.wikipedia.org/wiki/Centroid#Centroid_of_polygon
    area = sum(xs[i]*(ys[i+1]-ys[i-1]) for i in range(1, len(xs)-1))/2.0
    centroid_x =  sum((xs[i]+xs[i+1])*(xs[i]*ys[i+1] - xs[i+1]*ys[i]) for i in range(1, len(xs)-1))/(6.0*area)
    centroid_y =  sum((ys[i]+ys[i+1])*(xs[i]*ys[i+1] - xs[i+1]*ys[i]) for i in range(1, len(xs)-1))/(6.0*area)

    return (area, (centroid_x, centroid_y))

现在 NaN 的例子可以正常工作了:

#NaN Example
d = [[3.0 ,  4], [5.0,11],[  12.0,8],[   9.0 ,  5],[5,6] ]
print "number of vertices: {0}".format(len(d))

defects = []
defects.append([d[0], d[1], d[2], d[3], d[4] ])
defects.append([ d[4], d[3], d[2], d[1], d[0]])

for v in defects:
    print measure_polygon(np.array(v))

结果:

number of vertices: 5
(-30.0, (7.166666666666667, 7.6111111111111107))
(30.0, (7.166666666666667, 7.6111111111111107))

Answer 1

多边形必须是自闭合的，所以第一个点和最后一个点是相等的。这是非常标准的。您可以使用鞋带公式 ( https://en.m.wikipedia.org/wiki/Shoelace_formula) 与常规坐标，但如果我得到一个数据集缺少复制的最后一个点，我只是添加它.. 这使得计算更容易。因此，考虑一个没有由以下坐标定义的孔的多边形(来自引用)。请注意第一个点和最后一个点是相同的...如果不是，您将得到多部分多边形(例如带孔的多边形)的对齐错误

x = np.array([3,5,12,9,5,3]) # wikipedia
y= np.array([4,11,8,5,6,4])
a = np.array(list(zip(x,y)))
area1 = 0.5*np.abs(np.dot(x, np.roll(y, 1))-np.dot(y, np.roll(x, 1)))
area2 =0.5*np.abs(np.dot(x[1:], y[:-1])-np.dot(y[1:], x[:-1]))
print("\nroll area {}\nslice area{}".format(area1, area2))

屈服

roll area 30.0
slice area30.0

现在您的多边形得到了相同的处理，将第一个点添加回作为最后一个点以关闭多边形

x = np.array([-148.35290745, -124.93580616, -0.58281373,  8.77029032, -148.35290745])
y = np.array([-1.95467472, -2.09420039,  1.32530292,  22.79390931, -1.95467472])
roll area  1619.5826480482792
slice area 1619.5826480482792

区域结果与您的不同，但我使用使用 einsum 的第三种方法确认了它。下面是脚本的一部分

def ein_area(a, b=None):
    """Area calculation, using einsum.
    :Requires:
    :--------
    :  a - either a 2D+ array of coordinates or an array of x values
    :  b - if a < 2D, then the y values need to be supplied
    :  Outer rings are ordered clockwise, inner holes are counter-clockwise
    :Notes:
    :  x => array([ 0.000,  0.000,  10.000,  10.000,  0.000])  .... OR ....
    :  t = x.reshape((1,) + x.shape)
    :      array([[ 0.000,  0.000,  10.000,  10.000,  0.000]]) .... OR .... 
    :  u = np.atleast_2d(x)
    :      array([[ 0.000,  0.000,  10.000,  10.000,  0.000]]) .... OR ....
    :  v = x[None, :]
    :      array([[ 0.000,  0.000,  10.000,  10.000,  0.000]])
    """
    a = np.array(a)  
    if b is None:
        xs = a[..., 0]
        ys = a[..., 1]
    else:
        xs, ys = a, b
    x0 = np.atleast_2d(xs[..., 1:])
    y0 = np.atleast_2d(ys[..., :-1])
    x1 = np.atleast_2d(xs[..., :-1])
    y1 = np.atleast_2d(ys[..., 1:])
    e0 = np.einsum('...ij,...ij->...i', x0, y0)
    e1 = np.einsum('...ij,...ij->...i', x1, y1)
    area = abs(np.sum((e0 - e1)*0.5))
    return area

但您可以看到，这主要基于切片/滚动方法。我会检查您是否可以通过包括多边形列表中通常缺少但假设的最后一个点来确认结果。