ポリゴン計算の範囲が一貫していません

Question

のvoronoi_finite_polygons_2d(vor, radius=None)関数を使用していました。

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)) 

最初の例は、処理順序（cwまたは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)) 

簡単な平行四辺形出力：今

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）は非常に大きいので、浮動小数点丸めの対象ではないと考えられます。しかし、それ以外に、私はこれがなぜ機能していないのかという仮説を使い果たしました。

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

私は私のリスト - ....エラーを検出しました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の例では、右の作品：最初と最後の点があるので、自己閉鎖する必要が

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

Answer 1

ポリゴン：

#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)) 

結果を等しい。これはかなり標準的です。通常の座標で靴ひも式（https://en.m.wikipedia.org/wiki/Shoelace_formula ）を使用することはできますが、複製された最後の点が欠落しているデータセットを取得したら、追加するだけで計算が簡単になります。したがって、以下の座標（参照からの）で定義された穴のないポリゴンを考えてみましょう。最初と最後のポイントは同じであることに注意してください...そうでない場合は、今すぐあなたのポリゴンが裏で最初のポイントを追加することで、同じように処理した

roll area 30.0 
slice area30.0 

をもたらすマルチパートポリゴン（例えば穴のあるポリゴン）

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)) 

のためのアライメントエラーを取得しますポリゴン

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を用いた第3の方法を用いて確認しました。スクリプトの一部は次のとおりです

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 

しかし、これは主にスライス/ローリングのアプローチに基づいています。私はあなたが通常ポリゴンリストから欠落していると思われる最後のポイントを含めることによって結果を確認できるかどうかを確認します。

Answer 2

理由は最後の点がありません。 最初のものとほぼ同じですが、ポリゴンは回路でなければなりません。

データベースでは、通常、最初のデータベースと同じにすることが標準であるため、データベースは省略されます。