Python中的仿射3D变换
我正在用Python在Autodesk Maya里编写一个函数(使用PyMel)。
我有三个三维点:p0、p1和p2。
然后它们进行了一次刚性变换,所以在变换之后,我得到了它们的新位置:q0、q1和q2。
我还有一个第四个点,变换之前是p3。我想计算它在同样变换后的位置,记作q4。
所以我需要计算一个变换矩阵,然后把它应用到p4上。但我不知道该怎么做。List = 一组对象的数组。
import pymel.core as pm
import pymel.core.datatypes as dt
p0 = dt.Vector(pm.getAttr(list[0]+".tx"), pm.getAttr(list[0]+".ty"), pm.getAttr(list[0]+".tz"))
p1 = dt.Vector(pm.getAttr(list[1]+".tx"), pm.getAttr(list[1]+".ty"), pm.getAttr(list[1]+".tz"))
p2 = dt.Vector(pm.getAttr(list[2]+".tx"), pm.getAttr(list[2]+".ty"), pm.getAttr(list[2]+".tz")
p3 = dt.Vector(pm.getAttr(list[3]+".tx"), pm.getAttr(list[3]+".ty"), pm.getAttr(list[3]+".tz"))
这些三维点是从Maya场景中的动画对象读取的。所以在另一个帧上,我运行这段代码来获取
q0 = dt.Vector(pm.getAttr(list[0]+".tx"), pm.getAttr(list[0]+".ty"), pm.getAttr(list[0]+".tz"))
q1 = dt.Vector(pm.getAttr(list[1]+".tx"), pm.getAttr(list[1]+".ty"), pm.getAttr(list[1]+".tz"))
q2 = dt.Vector(pm.getAttr(list[2]+".tx"), pm.getAttr(list[2]+".ty"), pm.getAttr(list[2]+".tz"))
#q3 = TransformationMatrix between (p0,p1,p2) and (q0,q1,q2), applied to p3
我尝试用向量计算,但因为出现了除以零的错误,最后搞得一团糟……所以我觉得用变换矩阵应该能顺利解决这个问题。
我的截止日期快到了,我真的需要帮助来解决这个问题!请帮帮我!
编辑:
如何用Python进行坐标的仿射变换?我需要这个函数“solve_affine”,但它应该只从每组中取3个点,而不是4个。而且我不能使用numpy……
3 个回答
我已经搞定了
p0p1 = p1-p0
p0p2 = p2-p0
p0p3 = p3-p0
q0q1 = q1-q0
q0q2 = q2-q0
q0q3 = q3-q0
before = dt.Matrix([p0.x, p0.y, p0.z, 0],[p1.x, p1.y, p1.z, 0],[p2.x, p2.y, p2.z, 0], [0,0,0,1]);
after = dt.Matrix([q0.x, q0.y, q0.z, 0],[q1.x, q1.y, q1.z, 0],[q2.x, q2.y, q2.z, 0], [0,0,0,1]);
normal = p0p1.cross(p0p2).normal()
dist = p0p3.dot(normal)
q3 = p3 - dist*normal
transformMatrix = before.inverse()*after
solve = dt.Matrix(q3.x, q3.y, q3.z, 1)*transformMatrix
q3 = dt.Vector(solve[0][0], solve[0][1], solve[0][2])
newNormal = q0q1.cross(q0q2).normal()
q3 = q3 + newNormal*dist
pm.move(list[3], q3, r=False)
这个变换矩阵只对在平面p0p1p2内的点有效。所以我通过先变换p3的投影点,然后再把它按相同的距离移出平面来解决这个问题。
如果你有只用矩阵就能解决的办法,欢迎分享,可能对我还有帮助! :)
你的描述和代码都有点让人困惑。描述有点模糊,而代码示例缺少一些重要的部分。所以我理解这个问题是这样的:
已知两个空间中的三个点,如何构建从空间A到空间B的变换?
图片1: 如何在两个空间之间形成变换。
答案取决于这两个空间的变换类型。你会发现,三个点总是形成一个平面。这意味着你可以知道新空间的旋转、变换和均匀缩放。你还可以知道平面上的剪切和非均匀缩放。不过,你无法知道在法线方向上的剪切或非均匀缩放。
因此,问题变成了如何旋转和移动这两个空间以使它们匹配?这其实很简单,移动的部分直接可以用:
trans = q0 - p0
这就剩下旋转的部分了,关于旋转的解释在几个帖子中都有:
在这之后,你还可以计算一个缩放因子。
这里有一个使用numpy和scipy的解决方案。scipy主要用来生成随机旋转,除了scipy.linalg.norm这个功能,自己写也很简单。numpy主要用到的功能是叉乘和矩阵乘法,这些也都不难自己实现。
基本的思路是这样的:给定三个不在同一条直线上的点x1、x2、x3,可以找到一组三个互相垂直的向量(坐标轴)v1、v2、v3。v1的方向是x2-x1,v2在由(x2-x1)和(x3-x1)构成的平面内,而v3则用来补全这组三个向量。
点y1、y2、y3是相对于x1、x2、x3进行旋转和移动的。由y1、y2、y3生成的坐标轴w1、w2、w3是从v1、v2、v3旋转过来的(也就是说没有移动)。这两组向量都是互相垂直的,所以很容易从中找到旋转关系:R = W * transpose(V)
一旦我们找到了旋转,求移动就简单了:y1 = R*x + t
,所以t = y1 - R*x
。可能更好的方法是使用最小二乘法求解,把三个点结合起来估算t。
import numpy
import scipy.linalg
def rand_rot():
"""Return a random rotation
Return a random orthogonal matrix with determinant 1"""
q, _ = scipy.linalg.qr(numpy.random.randn(3, 3))
if scipy.linalg.det(q) < 0:
# does this ever happen?
print "got a negative det"
q[:, 0] = -q[:, 0]
return q
def rand_noncollinear():
"""Return 3 random non-collinear vectors"""
while True:
b = numpy.random.randn(3, 3)
sigma = scipy.linalg.svdvals(b)
if sigma[2]/sigma[0] > 0.1:
# "very" non-collinear
break
# "nearly" collinear; try again
return b[:, 0], b[:, 1], b[:, 2]
def normalize(a):
"""Return argument normalized"""
return a/scipy.linalg.norm(a)
def colstack(a1, a2, a3):
"""Stack three vectors as columns"""
return numpy.hstack((a1[:, numpy.newaxis],
a2[:, numpy.newaxis],
a3[:, numpy.newaxis]))
def get_axes(a1, a2, a3):
"""Generate orthogonal axes from three non-collinear points"""
# I tried to do this with QR, but something didn't work
b1 = normalize(a2-a1)
b2 = normalize(a3-a1)
b3 = normalize(numpy.cross(b1, b2))
b4 = normalize(numpy.cross(b3, b1))
return b1, b4, b3
# random rotation and translation
r = rand_rot()
t = numpy.random.randn(3)
# three non-collinear points
x1, x2, x3 = rand_noncollinear()
# some other point
x4 = numpy.random.randn(3)
# the images of the above in the transformation.
# y4 is for checking only -- won't be used to estimate r or t
y1, y2, y3, y4 = [numpy.dot(r, x) + t
for x in x1, x2, x3, x4]
v1, v2, v3 = get_axes(x1, x2, x3)
w1, w2, w3 = get_axes(y1, y2, y3)
V = colstack(v1, v2, v3)
W = colstack(w1, w2, w3)
# W = R V, so R = W * inverse(V); but V orthogonal, so inverse(V) is
# transpose(V):
rfound = numpy.dot(W, V.T)
# y1 = R x1 + t, so...
tfound = y1-numpy.dot(r, x1)
# get error on images of x2 and x3, just in case
y2err = scipy.linalg.norm(numpy.dot(rfound, x2) + tfound - y2)
y3err = scipy.linalg.norm(numpy.dot(rfound, x3) + tfound - y3)
# and check error image of x4 -- getting an estimate of y4 is the
# point of all of this
y4err = scipy.linalg.norm(numpy.dot(rfound, x4) + tfound - y4)
print "y2 error: ", y2err
print "y3 error: ", y3err
print "y4 error: ", y4err