I was using voronoi_finite_polygons_2d(vor, radius=None) function that I found elsewhere on StackOverflow.
I want to modify it to show the centroid of each voronoi cell. Debugging why some centroids appear dramatically wrong (see green arror pointing out centroid way off in the weeds). First error I identified: some of the calculations weren't processing the vertices in proper all-clockwise or all-counterclockwise order.
Not sure why some points don't get sorted correctly, but before I investigate that, I found another anomaly.
I should get the same area (with opposite sign) if I go clockwise or counterclockwise. In simple examples, I do. But in a random polygon I made, I get /slightly/ different results.
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))
The first example work as expect -- same area and centroid regardless of processing order (cw or ccw).
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))
simple parallelogram output:
4
(-12.0, (2.5, 1.5))
(12.0, (2.5, 1.5))
But now look at this 4-sided polygon (that is almost a triangle)
#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))
Gives me weird output:
4
(1280.4882517358433, (-36.609159411740798, 7.5961622623413145))
(-1278.8546083623708, (-36.655924939495335, 7.6058658049196115))
The areas are different. And if areas are different, then the centroids will be different. The discrepancies of area (1280 versus 1278) are so large that I doubt it's a floating point rounding thing. But other than that, I've run out of hypotheses why this isn't working.
===============================
I found the error.... my list-comprehension/indexing hack to enable y-1 and y+1 notation was broken (in a sinister way that half-worked). The correct routine is as follows:
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))
So now NaN's example works right:
#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))
results:
number of vertices: 5
(-30.0, (7.166666666666667, 7.6111111111111107))
(30.0, (7.166666666666667, 7.6111111111111107))
