Numpy - 坐标系之间的转换

T81*_*T81 6 python 3d numpy linear-algebra coordinate-systems

使用Numpy我想在坐标系之间转换位置矢量.

为了帮助可视化问题:http://tube.geogebra.org/student/m1097765

我在3D空间有两架飞机.每个平面由其中心定义:

C[0] = (X0, Y0, Z0)

C[1] = (X1, Y1, Z1)
Run Code Online (Sandbox Code Playgroud)

(X,Y,Z称为全局坐标系)

C = np.array([[0,0,0],[-4,2,1]])
Run Code Online (Sandbox Code Playgroud)

和它的法线向量:

H[0] = (cos(alpha[0])*sin(A[0]), cos(alpha[0])*cos(A[0]), sin(A[0])

H[1] = (cos(alpha[1])*sin(A[1]), cos(alpha[1])*cos(A[1]), sin(A[1])
Run Code Online (Sandbox Code Playgroud)

alpha =仰角

A =方位角

H = np.array([[-0.23, -0.45, 0.86], [-0.12, -0.24, 0.86]])
Run Code Online (Sandbox Code Playgroud)

我有一个点p(xp, yp, 0)趴在平面0(xp,yp被称为局部坐标系与中心C[0]和它的xyz轴与全球一致XYZ的轴时alpha = A = 0)

我使用以下函数从平面0的局部坐标系转换为全局坐标系:

import numpy as np

def rotateAxisX(alpha):
    '''
    Rotation about x axis
    :param alpha: plane altitude angle in degrees
    :return: x-axis rotation matrix
    '''
    rotX = np.array([[1, 0, 0], [0, np.cos(np.deg2rad(alpha)), np.sin(np.deg2rad(alpha))], [0, -np.sin(np.deg2rad(alpha)), np.cos(np.deg2rad(alpha))]])
    return rotX

def rotateAxisZ(A):
    '''
    Rotation about z axis
    :param A: plane azimuth angle in degrees
    :return: z-axis rotation matrix
    '''
    rotZ = np.array([[np.cos(np.deg2rad(A)), np.sin(np.deg2rad(A)), 0], [-np.sin(np.deg2rad(A)), np.cos(np.deg2rad(A)), 0], [0, 0, 1]])
    return rotZ

def local2Global(positionVector, planeNormalVector, positionVectorLocal):
    '''
    Convert point from plane's local coordinate system to global coordinate system
    :param positionVector: plane center in global coordinates
    :param planeNormalVector: the normal vector of the plane
    :param positionVectorLocal: a point on plane (xp,yp,0) with respect to the local coordinate system of the plane
    :return: the position vector of the point in global coordinates 
    >>> C = np.array([-10,20,1200]) 
    >>> H = np.array([-0.23, -0.45, 0.86])
    >>> p = np.array([-150, -1.5, 0])
    >>> P = local2Global(C, H, p)
    >>> np.linalg.norm(P-C) == np.linalg.norm(p)
    True
    '''
    alpha = np.rad2deg(np.arcsin(planeNormalVector[2]))
    A = np.where(planeNormalVector[1] > 0, np.rad2deg(np.arccos(planeNormalVector[1] / np.cos(np.deg2rad(alpha)))), 360 - np.rad2deg(np.arccos(planeNormalVector[1] / np.cos(np.deg2rad(alpha)))))
    positionVectorGlobal = positionVector + np.dot(np.dot(rotateAxisZ(A), rotateAxisX(90 - alpha)), positionVectorLocal)
    return positionVectorGlobal
Run Code Online (Sandbox Code Playgroud)

以上似乎按预期工作.

然后我计算从平面0上的点经过的线的交点,p(xp,yp,0)并且具有方向向量S = (0.56, -0.77, 0.3)

>>> C = np.array([[0,0,0],[-4,2,1]]) # plane centers
>>> H = np.array([[-0.23, -0.45, 0.86], [-0.12, -0.24, 0.86]]) # plane normal vectors
>>> S = np.array([0.56, -0.77, 0.3]) # a direction vector
>>> p = np.array([-1.5, -1.5, 0]) # a point on a plane

>>> intersectingPlaneIndex = 0 # choose intersecting plane, this plane has the point p on it
>>> intersectedPlaneIndex = 1 # this plane intersects with the line passing from p with direction vector s

>>> P = local2Global(C[intersectingPlaneIndex], H[intersectingPlaneIndex], p)   # point p in global coordinates
>>> np.isclose(np.linalg.norm(p), np.linalg.norm(P - C[intersectingPlaneIndex]), 10e-8)
True
Run Code Online (Sandbox Code Playgroud)

所以第一次转型是成功的.

现在让我们在全局坐标中找到交叉点E.

>>> t = np.dot(H[intersectedPlaneIndex], C[intersectedPlaneIndex, :] - P) / np.dot(H[intersectedPlaneIndex], S)
>>> E = P + S * t
>>> np.around(E, 2)
array([ 2.73, -0.67,  1.19])
Run Code Online (Sandbox Code Playgroud)

到目前为止,我找到了E位于平面1上的点(全局坐标).

问题:

如何将点E从全局坐标转换为平面1的坐标系并获得e(xe, ye, 0)

我试过了:

def global2Local(positionVector, planeNormalVector, positionVectorGlobal):
    '''
    Convert point from global coordinate system to plane's local coordinate system
    :param positionVector: plane center in global coordinates
    :param planeNormalVector: the normal vector of the plane
    :param positionVectorGlobal: a point in global coordinates
    :note: This function translates the given position vector by the positionVector and rotates the basis axis in order to obtain the positionVectorCoordinates in plane's coordinate system
    :warning: it does not function as it should
    '''
    alpha = np.rad2deg(np.arcsin(planeNormalVector[2]))
    A = np.where(planeNormalVector[1] > 0, np.rad2deg(np.arccos(planeNormalVector[1] / np.cos(np.deg2rad(alpha)))), 360 - np.rad2deg(np.arccos(planeNormalVector[1] / np.cos(np.deg2rad(alpha)))))
    positionVectorLocal = np.dot(np.dot(np.linalg.inv(rotateAxisZ(A)), np.linalg.inv(rotateAxisX(90 - alpha))), positionVectorGlobal - positionVector) + positionVectorGlobal
    return positionVectorLocal
Run Code Online (Sandbox Code Playgroud)

和:

>>> e = global2Local(C[intersectedPlaneIndex], H[intersectedPlaneIndex], E)
>>> e
array([ -2.54839059e+00,  -5.48380179e+00,  -1.42292121e-03])
Run Code Online (Sandbox Code Playgroud)

首先看,这似乎没问题,只要e [2]接近零,但是,

>>> np.linalg.norm(E-C[intersectedPlaneIndex])
7.2440723159783182
>>> np.linalg.norm(e)
6.0470140356703537
Run Code Online (Sandbox Code Playgroud)

所以转型是错误的.有任何想法吗?

小智 1

我建议阅读这个这个。对于第一个,看看齐次坐标的概念,对于不同原点的空间变换,这是有点需要的。对于第二个,看看相机“注视”变换是如何执行的。只要你有正交基向量(很容易从角度获得),你就可以使用第二个方程中的方程来进行变换。评论中链接的帖子似乎涵盖了类似的材料。