这是用于计算3D平面多边形面积的公式的推导
这是实现它的Python代码:
#determinant of matrix adef 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 cdef 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 bdef dot(a, b): return a[0]*b[0] + a[1]*b[1] + a[2]*b[2]#cross product of vectors a and bdef 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 polydef 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
最初的问题是我过于简化了。它需要计算垂直于平面的单位矢量。面积是该乘积和所有叉积的总和的一半,而不是叉积所有量值之和的一半。
可以稍微清理一下(如果有矩阵和向量类,或者行列式/叉积/点积的标准实现,矩阵和向量类会更好),但是从概念上讲应该是合理的。



