python - Area of polygon calculation inconsistent -
i using voronoi_finite_polygons_2d(vor, radius=none) function found elsewhere on stackoverflow.
want modify show centroid of each voronoi cell. debugging why centroids appear dramatically wrong (see green arror pointing out centroid way off in weeds). first error identified: of calculations weren't processing vertices in proper all-clockwise or all-counterclockwise order.
not sure why points don't sorted correctly, before investigate that, found anomaly.
i should same area (with opposite sign) if go clockwise or counterclockwise. in simple examples, do. in random polygon made, /slightly/ different results.
import numpy np import matplotlib.pyplot plt 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]) 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]) 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]) in range(0, len(xs)-1))/(6.0*area) return (area, (centroid_x, centroid_y)) the first example work expect -- same area , 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]]) 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 @ 4-sided polygon (that 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]]) 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 different. , if areas different, centroids different. discrepancies of area (1280 versus 1278) large doubt it's floating point rounding thing. other that, i've run out of hypotheses why isn't working.
===============================
i found error.... list-comprehension/indexing hack enable y-1 , y+1 notation broken (in sinister way half-worked). correct routine follows:
def measure_polygon(vertices): xs = vertices[:,0] ys = vertices[:,1] #the first , last elements +1 -1 work @ 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 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]) 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]) 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]) in range(1, len(xs)-1))/(6.0*area) return (area, (centroid_x, centroid_y)) so 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]]) 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))
polygons have self closing first , last point equal. pretty standard. can use shoelace formula ( https://en.m.wikipedia.org/wiki/shoelace_formula ) regular coordinates, if dataset missing replicated last point, add it.. makes calculations easier. consider polygon without holes defined following coordinates (from reference). note first , last point same...if aren't alignment errors multipart polygons (polygons holes instance)
x = np.array([3,5,12,9,5,3]) # wikipedia y= np.array([4,11,8,5,6,4]) = 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)) yielding
roll area 30.0 slice area30.0 now polygon treated same, adding first point in last point give closure polygon
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 the area result differs yours confirmed using third method using einsum. portion of script below
def ein_area(a, b=none): """area calculation, using einsum. :requires: :-------- : - either 2d+ array of coordinates or array of x values : b - if < 2d, y values need supplied : outer rings ordered clockwise, inner holes 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]]) """ = np.array(a) if b 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 but can see largely based on slicing/rolling approach. check see if can confirm results including last point missing polygon lists assumed.
Comments
Post a Comment