如何校正 3D 平面多边形?
How do you rectify a 3D planar polygon?
我有一个 3D 平面(所有顶点都位于某个平面内)多边形,其顶点为:[(x1, y1, z1) ... (x1, y1, z1)]。
我想转换这个多边形,以便我可以正交地查看它(就像我在直视它一样)。
如何在 Python 中完成?
我假设你除了顶点坐标之外没有其他信息。
取三个非共线(可能是后续)顶点C, A, B
。计算归一化向量(除以向量长度)
b = (B - A) / |B - A|
然后法向量(使用vector/cross乘法)
N = b.cross.(A-C) and normalize it
un = N / |N|
和多边形平面中的另一个单位向量
v = b.cross.n
现在我们要找到这样的仿射变换矩阵,将顶点 A
变换为点 (0,0,0)
,边 AB
将与 OX
轴共线,法线将与 OZ
轴共线,矢量 q
将与 OY
轴共线。这一切都意味着旋转的多边形将位于 OXY
平面内。
数学上:点 A, u=A+b, v=A+q, n=A+un
应转换为四元组 (0,0,0), (1,0,0), (0,1,0), (0,0,1)
。矩阵形式
[Ax ux vx nx] [0 1 0 0]
M * [Ay uy vy ny] = [0 0 1 0]
[Az uz vz nz] [0 0 0 1]
[1 1 1 1 ] [1 1 1 1]
或
M * S = D
使用矩阵求逆
M * S * Sinv = D * Sinv
最后
M = D * Sinv
因此计算矩阵 M
并将其与每个顶点坐标相乘。新坐标应具有零 Z 分量(或由于数值错误而非常小)。
您可以使用 numpy
库执行所有描述的操作,只需一点代码
简单的快速实现Python供参考
import math
def calcMatrix(ax, bx, cx, ay, by, cy, az, bz, cz):
ux, uy, uz = bx - ax, by - ay, bz - az
mag = math.sqrt(ux*ux+uy*uy +uz*uz)
ux, uy, uz = ux / mag, uy / mag, uz / mag
Cx, Cy, Cz = ax - cx, ay - cy, az - cz
nx, ny, nz = uy * Cz - uz * Cy, uz * Cx - ux * Cz, ux * Cy - uy * Cx
mag = math.sqrt(nx*nx+ny*ny+nz*nz)
nx, ny, nz = nx / mag, ny / mag, nz / mag
vx, vy, vz = uy * nz - uz * ny, uz * nx - ux * nz, ux * ny - uy * nx
denom = 1.0 / (ux*ux+uy*uy + uz*uz)
M = [[0.0]*4 for _ in range(4)]
M[3][3] = 1.0
M[0][0] = denom*(ux)
M[0][1] = denom*(uy)
M[0][2] = denom*(uz)
M[0][3] = denom*(-ax*ux-ay*uy+az*uz)
M[1][0] = denom*(vx)
M[1][1] = denom*(vy)
M[1][2] = denom*(vz)
M[1][3] = -denom*(ax*vx-ay*vy+az*vz)
M[2][0] = denom*(nx)
M[2][1] = denom*(ny)
M[2][2] = denom*(nz)
M[2][3] = denom*(-ax*nx-ay*ny+az*nz)
return M
def mult(M, vec):
res = [0]*4
for k in range(4):
for i in range(4):
res[k] += M[k][i] * vec[i]
return res
#test corners and middle point
M = calcMatrix(1, 0, 0, 0, 1, 0, 0, 0, 1)
#print(M)
p = [1, 0, 0, 1]
print(mult(M, p))
p = [0, 1, 0, 1]
print(mult(M, p))
p = [0, 0, 1, 1]
print(mult(M, p))
p = [1/3, 1/3, 1/3, 1]
print(mult(M, p))
测试结果:
[0.0, 0.0, 0.0, 1.0]
[1.4142135623730951, 0.0, 0.0, 1.0]
[0.7071067811865476, 1.2247448713915892, 0.0, 1.0]
[0.7071067811865476, 0.4082482904638631, 1.1102230246251565e-16, 1.0]
通过两个非平行边之间的叉积,找到多边形的法线 n
。将 n
与垂直向量进行叉积,得到水平向量 u
。然后取 n
和 u
的叉积得到 v
,并对向量进行归一化。 u
和 v
平行于多边形的平面并且彼此正交。
最后,对于每个顶点 p
计算 2D 坐标 (p.u, p.v)
,它会显示平面中的多边形。
numpy
提供 cross
和 dot
向量函数。还有 linalg.norm
(或 sqrt(dot(v, v))
)。
这是使用 NumPy 的可靠方法(project()
;其余为测试代码)。
import numpy
import scipy.spatial
def project(x):
# Center the plane on the origin
x = x - numpy.mean(x, axis=0)
# Compute the Singular Value Decomposition
u, s, v = numpy.linalg.svd(x)
# Return the top two principal components
return u[:, :2] @ numpy.diag(s[:2])
def test():
n = 10
x = (numpy.random.rand(n, 2) @ numpy.random.rand(2, 3)) + numpy.random.rand(3)
y = project(x)
print(x.shape, y.shape)
print(
numpy.max(
numpy.abs(
scipy.spatial.distance_matrix(x, x)
- scipy.spatial.distance_matrix(y, y)
)
)
)
if __name__ == "__main__":
test()
示例输出:
(10, 3) (10, 2)
5.551115123125783e-16
我有一个 3D 平面(所有顶点都位于某个平面内)多边形,其顶点为:[(x1, y1, z1) ... (x1, y1, z1)]。
我想转换这个多边形,以便我可以正交地查看它(就像我在直视它一样)。
如何在 Python 中完成?
我假设你除了顶点坐标之外没有其他信息。
取三个非共线(可能是后续)顶点C, A, B
。计算归一化向量(除以向量长度)
b = (B - A) / |B - A|
然后法向量(使用vector/cross乘法)
N = b.cross.(A-C) and normalize it
un = N / |N|
和多边形平面中的另一个单位向量
v = b.cross.n
现在我们要找到这样的仿射变换矩阵,将顶点 A
变换为点 (0,0,0)
,边 AB
将与 OX
轴共线,法线将与 OZ
轴共线,矢量 q
将与 OY
轴共线。这一切都意味着旋转的多边形将位于 OXY
平面内。
数学上:点 A, u=A+b, v=A+q, n=A+un
应转换为四元组 (0,0,0), (1,0,0), (0,1,0), (0,0,1)
。矩阵形式
[Ax ux vx nx] [0 1 0 0]
M * [Ay uy vy ny] = [0 0 1 0]
[Az uz vz nz] [0 0 0 1]
[1 1 1 1 ] [1 1 1 1]
或
M * S = D
使用矩阵求逆
M * S * Sinv = D * Sinv
最后
M = D * Sinv
因此计算矩阵 M
并将其与每个顶点坐标相乘。新坐标应具有零 Z 分量(或由于数值错误而非常小)。
您可以使用 numpy
库执行所有描述的操作,只需一点代码
简单的快速实现Python供参考
import math
def calcMatrix(ax, bx, cx, ay, by, cy, az, bz, cz):
ux, uy, uz = bx - ax, by - ay, bz - az
mag = math.sqrt(ux*ux+uy*uy +uz*uz)
ux, uy, uz = ux / mag, uy / mag, uz / mag
Cx, Cy, Cz = ax - cx, ay - cy, az - cz
nx, ny, nz = uy * Cz - uz * Cy, uz * Cx - ux * Cz, ux * Cy - uy * Cx
mag = math.sqrt(nx*nx+ny*ny+nz*nz)
nx, ny, nz = nx / mag, ny / mag, nz / mag
vx, vy, vz = uy * nz - uz * ny, uz * nx - ux * nz, ux * ny - uy * nx
denom = 1.0 / (ux*ux+uy*uy + uz*uz)
M = [[0.0]*4 for _ in range(4)]
M[3][3] = 1.0
M[0][0] = denom*(ux)
M[0][1] = denom*(uy)
M[0][2] = denom*(uz)
M[0][3] = denom*(-ax*ux-ay*uy+az*uz)
M[1][0] = denom*(vx)
M[1][1] = denom*(vy)
M[1][2] = denom*(vz)
M[1][3] = -denom*(ax*vx-ay*vy+az*vz)
M[2][0] = denom*(nx)
M[2][1] = denom*(ny)
M[2][2] = denom*(nz)
M[2][3] = denom*(-ax*nx-ay*ny+az*nz)
return M
def mult(M, vec):
res = [0]*4
for k in range(4):
for i in range(4):
res[k] += M[k][i] * vec[i]
return res
#test corners and middle point
M = calcMatrix(1, 0, 0, 0, 1, 0, 0, 0, 1)
#print(M)
p = [1, 0, 0, 1]
print(mult(M, p))
p = [0, 1, 0, 1]
print(mult(M, p))
p = [0, 0, 1, 1]
print(mult(M, p))
p = [1/3, 1/3, 1/3, 1]
print(mult(M, p))
测试结果:
[0.0, 0.0, 0.0, 1.0]
[1.4142135623730951, 0.0, 0.0, 1.0]
[0.7071067811865476, 1.2247448713915892, 0.0, 1.0]
[0.7071067811865476, 0.4082482904638631, 1.1102230246251565e-16, 1.0]
通过两个非平行边之间的叉积,找到多边形的法线 n
。将 n
与垂直向量进行叉积,得到水平向量 u
。然后取 n
和 u
的叉积得到 v
,并对向量进行归一化。 u
和 v
平行于多边形的平面并且彼此正交。
最后,对于每个顶点 p
计算 2D 坐标 (p.u, p.v)
,它会显示平面中的多边形。
numpy
提供 cross
和 dot
向量函数。还有 linalg.norm
(或 sqrt(dot(v, v))
)。
这是使用 NumPy 的可靠方法(project()
;其余为测试代码)。
import numpy
import scipy.spatial
def project(x):
# Center the plane on the origin
x = x - numpy.mean(x, axis=0)
# Compute the Singular Value Decomposition
u, s, v = numpy.linalg.svd(x)
# Return the top two principal components
return u[:, :2] @ numpy.diag(s[:2])
def test():
n = 10
x = (numpy.random.rand(n, 2) @ numpy.random.rand(2, 3)) + numpy.random.rand(3)
y = project(x)
print(x.shape, y.shape)
print(
numpy.max(
numpy.abs(
scipy.spatial.distance_matrix(x, x)
- scipy.spatial.distance_matrix(y, y)
)
)
)
if __name__ == "__main__":
test()
示例输出:
(10, 3) (10, 2)
5.551115123125783e-16