python xy坐标_Python：从xyz坐标中查找多边形区域

Here is the derivation of a formula for calculating the area of a 3D planar polygon

这是实现它的Python代码：

#determinant of matrix a

def det(a):

return a[0][0]*a[1][1]*a[2][2] + a[0][1]*a[1][2]*a[2][0] + a[0][2]*a[1][0]*a[2][1] - a[0][2]*a[1][1]*a[2][0] - a[0][1]*a[1][0]*a[2][2] - a[0][0]*a[1][2]*a[2][1]

#unit normal vector of plane defined by points a, b, and c

def unit_normal(a, b, c):

x = det([[1,a[1],a[2]],

[1,b[1],b[2]],

[1,c[1],c[2]]])

y = det([[a[0],1,a[2]],

[b[0],1,b[2]],

[c[0],1,c[2]]])

z = det([[a[0],a[1],1],

[b[0],b[1],1],

[c[0],c[1],1]])

magnitude = (x**2 + y**2 + z**2)**.5

return (x/magnitude, y/magnitude, z/magnitude)

#dot product of vectors a and b

def dot(a, b):

return a[0]*b[0] + a[1]*b[1] + a[2]*b[2]

#cross product of vectors a and b

def cross(a, b):

x = a[1] * b[2] - a[2] * b[1]

y = a[2] * b[0] - a[0] * b[2]

z = a[0] * b[1] - a[1] * b[0]

return (x, y, z)

#area of polygon poly

def area(poly):

if len(poly) < 3: # not a plane - no area

return 0

total = [0, 0, 0]

for i in range(len(poly)):

vi1 = poly[i]

if i is len(poly)-1:

vi2 = poly[0]

else:

vi2 = poly[i+1]

prod = cross(vi1, vi2)

total[0] += prod[0]

total[1] += prod[1]

total[2] += prod[2]

result = dot(total, unit_normal(poly[0], poly[1], poly[2]))

return abs(result/2)

为了测试它,这里是一个10×5平方,倾斜：

>>> poly = [[0, 0, 0], [10, 0, 0], [10, 3, 4], [0, 3, 4]]

>>> poly_translated = [[0+5, 0+5, 0+5], [10+5, 0+5, 0+5], [10+5, 3+5, 4+5], [0+5, 3+5, 4+5]]

>>> area(poly)

50.0

>>> area(poly_translated)

50.0

>>> area([[0,0,0],[1,1,1]])

0

最初的问题是我过度简化了.它需要计算垂直于平面的单位矢量.该面积是该点积的一半和所有交叉积的总和,而不是交叉积的所有量的总和的一半.

这可以稍微清理一下(矩阵和矢量类可以使它更好,如果你有它们,或者决定性/交叉积/点积的标准实现),但它应该在概念上合理.