Stosowanie rotacji kwaternionów do wektorów szeregów czasowych

Nov 24 2020

Mam serię czasową wektorów 3D w tablicy numpy w języku Python podobną do następującej:

array([[-0.062, -0.024,  1.   ],
       [-0.071, -0.03 ,  0.98 ],
       [-0.08 , -0.035,  0.991],
       [-0.083, -0.035,  0.98 ],
       [-0.083, -0.035,  0.977],
       [-0.082, -0.035,  0.993],
       [-0.08 , -0.034,  1.006],
       [-0.081, -0.032,  1.008],
       .......

Chcę obrócić każdy wektor wokół określonej osi pod określonym kątem theta. Używałem kwaternionów, aby to osiągnąć dla jednego wektora, jak znaleźć tutaj w odpowiedzi Henneray.

v1 = np.array ([1, -2, 0])
axis = np.array([-4, -2,  3])
theta = 1.5

rot_axis = np.insert(axis, 0, 0, axis=0)
axis_angle = (theta*0.5) * rot_axis/np.linalg.norm(rot_axis)
vec = quat.quaternion(*v1)
qlog = quat.quaternion(*axis_angle)
q = np.exp(qlog)
v_prime = q * vec * np.conjugate(q)
v_prime_vec = v_prime.imag

Moje pytanie brzmi: jaki jest najszybszy sposób na zastosowanie tego samego obrotu do każdego wektora w v1?

Nie można utworzyć quaternion z v1if v1zawiera tablicę wektorów 2D, więc mógłbym użyć pętli do obracania każdego elementu tablicy po kolei; Jednak w odpowiedzi Henneray w powyższym odsyłaczu wspomina się, że kwaternionie można zastosować do „odpowiednio wektoryzowanych macierzy numpy”. Czy ktoś ma jakieś sugestie, jak można to wdrożyć?

(Pytanie poboczne: jeśli my thetai axiszmienne byłyby tablicami o długości równej v1, czy można by użyć tej samej metody do obrócenia każdego wektora w v1 poprzez odpowiedni obrót?)

Odpowiedzi

1 henneray Nov 26 2020 at 11:52

Konieczne jest najpierw przekonwertowanie [x, y, z] wektorów kartezjańskich na 4 wektory o pierwszej składowej równej zero [0, x, y, z]. Następnie możesz rzutować to na tablicę kwaternionów, aby wykonać obliczenia wektoryzowane.

Poniższa funkcja pobiera tablicę wektorów kartezjańskich i obraca je wokół jednej osi obrotu. Będziesz musiał upewnić się, że norma tej osi jest równa kątowi obrotu theta.

def rotate_vectors(vecs, axis):
    """
    Rotate a list of 3D [x,y,z] vectors about corresponding 3D axis
    [x,y,z] with norm equal to the rotation angle in radians

    Parameters
    ----------
    vectors : numpy.ndarray with shape [n,3]
        list of [x,y,z] cartesian vector coordinates
    axis : numpy.ndarray with shape [3]
        [x,y,z] axis to rotate corresponding vectors about
    """
    # Make an 4 x n array of zeros
    vecs4 = np.zeros([vecs.shape[0],vecs.shape[1]+1])
    # Fill the imaginary i, j, k components with x, y, z values, leaving the real part w=0
    vecs4[:,1:] = vecs
    # Convert to quaternion array
    vecsq = quat.as_quat_array(vecs4)

    # Make a rotation quaternion
    qrot = quat.from_rotation_vector(axis)
    # Rotate vectors
    vecsq_rotated = qrot * vecsq * qrot.conjugate()
    # Cast quaternion array to float and return only imaginary components (ignore real part)
    return quat.as_float_array(vecsq_rotated)[:,1:]

Jako bonus, ta funkcja przyjmuje tablicę osi obrotu, aby obrócić każdy wektor o odpowiednie osie.

def rotate_vectors_each(vecs, axes):
    """
    Rotate a list of 3D [x,y,z] vectors about corresponding 3D axes
    [x,y,z] with norm equal to the rotation angle in radians

    Parameters
    ----------
    vectors : numpy.ndarray with shape [n,3]
        list of [x,y,z] cartesian vector coordinates
    axes : numpy.ndarray with shape [n,3]
        axes to rotate corresponding vectors about
        n = pulse shape time domain
        3 = [x,y,z]
    """
    # Make an 4 x n array of zeros
    vecs4 = np.zeros([vecs.shape[0],vecs.shape[1]+1])
    # Fill the imaginary i, j, k components with x, y, z values, leaving the real part w=0
    vecs4[:,1:] = vecs
    # Convert to quaternion array
    vecsq = quat.as_quat_array(vecs4)

    # Make an 4 x n array of zeros
    rots4 = np.zeros([rots.shape[0],rots.shape[1]+1])
    # Fill the imaginary i, j, k components with x, y, z values, leaving the real part w=0
    rots4[:,1:] = rots
    # Convert to quaternion array and take exponential
    qrots = np.exp(quat.as_quat_array(0.5 * rots4))

    # Rotate vectors
    vecsq_rotated = qrots * vecsq * qrots.conjugate()

    return quat.as_float_array(vecsq_rotated)[:,1:]

Zauważ, że przy tak wielu konwersjach między reprezentacją kąta osi a reprezentacją kwaternionów, da to niewielką poprawę wydajności w stosunku do algebry macierzy obrotu. Quaternions naprawdę przynoszą korzyści tylko wtedy, gdy obracasz wektor przez wiele sekwencyjnych obrotów, dzięki czemu możesz kumulować mnożenie kwaternionów.

1 JamesTursa Nov 25 2020 at 01:21

Jednym „szybkim” sposobem wykonania samego obliczenia rotacji byłoby przekształcenie kwaternionu w macierz cosinusów kierunków 3x3, umieszczenie wektorów w jednej ciągłej macierzy 3xN, a następnie wywołanie procedury z biblioteki BLAS (np. Dgemm) w celu wykonania standardu mnożenie macierzy. Dobra biblioteka BLAS z dużym N wykonałaby to wielowątkowe obliczenia.