我使用了voronoi_finite_polygons_2d(vor, radius=None)函数，这个函数是我在StackOverflow上找到的。
我想修改它以显示每个voronoi单元格的质心。调试时发现一些质心明显错误（见绿色箭头指出的偏离中心），首先我确定的错误是：某些计算没有按照逆时针或顺时针的顺序处理顶点。
不确定为什么有些点没有正确排序，但在我调查之前，我发现另一个异常。
如果我按顺时针或逆时针方向走，应该得到相同面积（带有相反的符号）。在简单的例子中，我确实得到了相同的结果。但在我制作的随机多边形中，我得到了略微不同的结果。
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))
现在看一下这个几乎是三角形的四边形
#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）非常大，我怀疑这不是浮点舍入问题。但除此之外，我已经没有假设为什么这不起作用了。
===============================
我发现错误了...我的列表解析/索引技巧，用于实现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))