从 xyz 坐标中查找多边形的面积
Posted
技术标签:
【中文标题】从 xyz 坐标中查找多边形的面积【英文标题】：Find area of polygon from xyz coordinates 【发布时间】：2012-09-20 11:18:33 【问题描述】：我正在尝试使用shapely.geometry.Polygon 模块来查找多边形的面积，但它在xy 平面上执行所有计算。这对我的一些多边形来说很好，但其他多边形也有z 维度，所以它不是我想要的。
是否有一个包可以为我提供来自xyz 坐标的平面多边形的面积，或者一个包或算法将多边形旋转到xy 平面以便我可以使用shapely.geometry.Polygon().area？
多边形以[(x1,y1,z1),(x2,y2,z3),...(xn,yn,zn)] 形式的元组列表表示。
【问题讨论】：
多边形是严格的二维图形。你到底想计算什么？ 我正在尝试从顶点的“xyz”坐标中找到建筑物屋顶和墙壁的表面积。 我还没有找到任何模块可以做到这一点，但你可以简单地将每个面向下投射到一个 xy 平面上，然后用你一直在使用的模块进行计算 你所说的“投降”是什么意思？ 只需旋转形状，直到它在 z 平面上平坦。 【参考方案1】：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)
为了测试它，这里有一个 10x5 的正方形，可以倾斜：
>>> 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
最初的问题是我过于简单化了。它需要计算垂直于平面的单位向量。面积是其点积和所有叉积总和的一半，而不是叉积所有量值总和的一半。
这可以稍微清理一下（如果你有矩阵和向量类，或者行列式/叉积/点积的标准实现，它会更好），但它应该在概念上是合理的。
【讨论】：
谢谢，汤姆。我找到了该页面以及一些将斯托克定理应用于 2D 多边形的示例代码，但无法使其适用于 3D。你的实现对我来说看起来不错。我只是在调整它以适应我的数据结构方式，即 [(x1,y1,z1),(x2,y2,z2),...]。area 函数应该是一样的。 cross_product_magnitude 将更改为 x = a[1] * b[2] - a[2] * b[1] 等。
是的，我知道了 - 但它会抛出太大的结果。我是否需要移动形状以使一个顶点位于原点？
你不应该这样做。我想我在某个地方搞砸了，我会调查一下。
为什么单位法线是通过行列式计算的？你不能只做多边形前两条边的叉积+归一化吗？【参考方案2】：
这是我使用的最终代码。它没有使用shapely，而是实现了Stoke's theorem直接计算面积。它建立在@Tom Smilack 的回答之上，该回答展示了如何在没有 numpy 的情况下做到这一点。
import numpy as np
#unit normal vector of plane defined by points a, b, and c
def unit_normal(a, b, c):
x = np.linalg.det([[1,a[1],a[2]],
[1,b[1],b[2]],
[1,c[1],c[2]]])
y = np.linalg.det([[a[0],1,a[2]],
[b[0],1,b[2]],
[c[0],1,c[2]]])
z = np.linalg.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)
#area of polygon poly
def poly_area(poly):
if len(poly) < 3: # not a plane - no area
return 0
total = [0, 0, 0]
N = len(poly)
for i in range(N):
vi1 = poly[i]
vi2 = poly[(i+1) % N]
prod = np.cross(vi1, vi2)
total[0] += prod[0]
total[1] += prod[1]
total[2] += prod[2]
result = np.dot(total, unit_normal(poly[0], poly[1], poly[2]))
return abs(result/2)
【讨论】：
我正在寻求实现这个解决方案，但不清楚的是为什么 unit_normal 函数实现了多边形的前 3 个点。 poly 是一个 3d 点列表，即原始问题中发布的元组列表。还是响应仅适用于 3 点多边形？谢谢 据我记忆，多边形上任意三个（非共线）点的单位法向量都是相同的，我们可以取前三个点并从中计算【参考方案3】：2D 多边形的面积可以使用 Numpy 作为单线计算...
poly_Area(vertices) = np.sum( [0.5, -0.5] * vertices * np.roll( np.roll(vertices, 1, axis=0), 1, axis=1) )
【讨论】：
这不适用于 3D 空间中的 2D 多边形，例如全部共面，但在 xyz 坐标中引用。【参考方案4】：仅供参考，这是 Mathematica 中的相同算法，带有一个婴儿单元测试
ClearAll[vertexPairs, testPoly, area3D, planeUnitNormal, pairwise];
pairwise[list_, fn_] := MapThread[fn, Drop[list, -1], Drop[list, 1]];
vertexPairs[Polygon[points___]] := Append[points, First[points]];
testPoly = Polygon[20, -30, 0, 40, -30, 0, 40, -30, 20, 20, -30, 20];
planeUnitNormal[Polygon[points___]] :=
With[ps = Take[points, 3],
With[p0 = First[ps],
With[qs = (# - p0) & /@ Rest[ps],
Normalize[Cross @@ qs]]]];
area3D[p : Polygon[polys___]] :=
With[n = planeUnitNormal[p], vs = vertexPairs[p],
With[areas = (Dot[n, #]) & /@ pairwise[vs, Cross],
Plus @@ areas/2]];
area3D[testPoly]
【讨论】：
planeUnitNormal 计算在前三个点共线的情况下不可靠。更智能的算法会选择三个不共线的点（由pairwise[...,Cross]=!=0 测试，如果找不到三个点则抛出。
@reb-cabin 为什么要扔？如果每三个点都是共线的，那么答案是零。【参考方案5】：
与@Tom Smilack 的回答相同，但在 javascript 中
//determinant of matrix a
function 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
function unit_normal(a, b, c)
let x = math.det([
[1, a[1], a[2]],
[1, b[1], b[2]],
[1, c[1], c[2]]
]);
let y = math.det([
[a[0], 1, a[2]],
[b[0], 1, b[2]],
[c[0], 1, c[2]]
]);
let z = math.det([
[a[0], a[1], 1],
[b[0], b[1], 1],
[c[0], c[1], 1]
]);
let magnitude = Math.pow(Math.pow(x, 2) + Math.pow(y, 2) + Math.pow(z, 2), 0.5);
return [x / magnitude, y / magnitude, z / magnitude];
// dot product of vectors a and b
function dot(a, b)
return a[0] * b[0] + a[1] * b[1] + a[2] * b[2];
// cross product of vectors a and b
function cross(a, b)
let x = (a[1] * b[2]) - (a[2] * b[1]);
let y = (a[2] * b[0]) - (a[0] * b[2]);
let z = (a[0] * b[1]) - (a[1] * b[0]);
return [x, y, z];
// area of polygon poly
function area(poly)
if (poly.length < 3)
console.log("not a plane - no area");
return 0;
else
let total = [0, 0, 0]
for (let i = 0; i < poly.length; i++)
var vi1 = poly[i];
if (i === poly.length - 1)
var vi2 = poly[0];
else
var vi2 = poly[i + 1];
let prod = cross(vi1, vi2);
total[0] = total[0] + prod[0];
total[1] = total[1] + prod[1];
total[2] = total[2] + prod[2];
let result = dot(total, unit_normal(poly[0], poly[1], poly[2]));
return Math.abs(result/2);
【讨论】：
"math.det" 应该只是 "det"【参考方案6】：#pythonn 3D 多边形区域代码（优化版）
def polygon_area(poly):
#shape (N, 3)
if isinstance(poly, list):
poly = np.array(poly)
#all edges
edges = poly[1:] - poly[0:1]
# row wise cross product
cross_product = np.cross(edges[:-1],edges[1:], axis=1)
#area of all triangles
area = np.linalg.norm(cross_product, axis=1)/2
return sum(area)
if __name__ == "__main__":
poly = [[0+5, 0+5, 0+5], [10+5, 0+5, 0+5], [10+5, 3+5, 4+5], [0+5, 3+5, 4+5]]
print(polygon_area(poly))
【讨论】：
【参考方案7】：感谢您提供详细的答案，但我有点惊讶没有简单的答案来获取该区域。
所以，我只是发布了一种使用 pyny3d 使用多边形或曲面的 3d 坐标计算面积的简化方法。
#Install pyny3d as:
pip install pyny3d
#Calculate area
import numpy as np
import pyny3d.geoms as pyny
coords_3d = np.array([[0, 0, 0],
[7, 0, 0],
[7, 10, 2],
[0, 10, 2]])
polygon = pyny.Polygon(coords_3d)
print(f'Area is : polygon.get_area()')
【讨论】：
以上是关于从 xyz 坐标中查找多边形的面积的主要内容，如果未能解决你的问题，请参考以下文章