Вращение 3D вектора?

Question

У меня есть два вектора в виде списков Python и угол.Например:

v = [3,5,0]
axis = [4,4,1]
theta = 1.2 #radian

Каков наилучший / самый простой способ получить результирующий вектор при вращении вектора v вокруг оси?

Вращение должно быть направлено против часовой стрелкина кого указывает вектор оси.Это называется Правило правой руки

Answer 1

Используя формулу Эйлера-Родригеса :

import numpy as np
import math

def rotation_matrix(axis, theta):
    """
    Return the rotation matrix associated with counterclockwise rotation about
    the given axis by theta radians.
    """
    axis = np.asarray(axis)
    axis = axis / math.sqrt(np.dot(axis, axis))
    a = math.cos(theta / 2.0)
    b, c, d = -axis * math.sin(theta / 2.0)
    aa, bb, cc, dd = a * a, b * b, c * c, d * d
    bc, ad, ac, ab, bd, cd = b * c, a * d, a * c, a * b, b * d, c * d
    return np.array([[aa + bb - cc - dd, 2 * (bc + ad), 2 * (bd - ac)],
                     [2 * (bc - ad), aa + cc - bb - dd, 2 * (cd + ab)],
                     [2 * (bd + ac), 2 * (cd - ab), aa + dd - bb - cc]])

v = [3, 5, 0]
axis = [4, 4, 1]
theta = 1.2 

print(np.dot(rotation_matrix(axis, theta), v)) 
# [ 2.74911638  4.77180932  1.91629719]

Answer 2

Однострочник с функциями numpy / scipy.

Мы используем следующее:

пусть a будет единичным вектором вдоль ось , т.е. a = ось / норма (ось)
и A = I × a - кососимметричная матрица, связанная с a то есть перекрестное произведение единичной матрицы на a

, тогда M = exp (θ A) - матрица вращения.

from numpy import cross, eye, dot
from scipy.linalg import expm, norm

def M(axis, theta):
    return expm(cross(eye(3), axis/norm(axis)*theta))

v, axis, theta = [3,5,0], [4,4,1], 1.2
M0 = M(axis, theta)

print(dot(M0,v))
# [ 2.74911638  4.77180932  1.91629719]

expm (код здесь) вычисляет ряд Тейлора экспоненты:
\sum_{k=0}^{20} \frac{1}{k!} (θ A)^k, так что это дорого, но читабельно и безопасно.Это может быть хорошим способом, если у вас мало поворотов, но много векторов.

Answer 3

Я просто хотел упомянуть, что если требуется скорость, завершение кода unutbu в weave.inline scipy и передача уже существующей матрицы в качестве параметра приводит к 20-кратному уменьшению времени выполнения.

Код(в revolution_matrix_test.py):

import numpy as np
import timeit

from math import cos, sin, sqrt
import numpy.random as nr

from scipy import weave

def rotation_matrix_weave(axis, theta, mat = None):
    if mat == None:
        mat = np.eye(3,3)

    support = "#include <math.h>"
    code = """
        double x = sqrt(axis[0] * axis[0] + axis[1] * axis[1] + axis[2] * axis[2]);
        double a = cos(theta / 2.0);
        double b = -(axis[0] / x) * sin(theta / 2.0);
        double c = -(axis[1] / x) * sin(theta / 2.0);
        double d = -(axis[2] / x) * sin(theta / 2.0);

        mat[0] = a*a + b*b - c*c - d*d;
        mat[1] = 2 * (b*c - a*d);
        mat[2] = 2 * (b*d + a*c);

        mat[3*1 + 0] = 2*(b*c+a*d);
        mat[3*1 + 1] = a*a+c*c-b*b-d*d;
        mat[3*1 + 2] = 2*(c*d-a*b);

        mat[3*2 + 0] = 2*(b*d-a*c);
        mat[3*2 + 1] = 2*(c*d+a*b);
        mat[3*2 + 2] = a*a+d*d-b*b-c*c;
    """

    weave.inline(code, ['axis', 'theta', 'mat'], support_code = support, libraries = ['m'])

    return mat

def rotation_matrix_numpy(axis, theta):
    mat = np.eye(3,3)
    axis = axis/sqrt(np.dot(axis, axis))
    a = cos(theta/2.)
    b, c, d = -axis*sin(theta/2.)

    return np.array([[a*a+b*b-c*c-d*d, 2*(b*c-a*d), 2*(b*d+a*c)],
                  [2*(b*c+a*d), a*a+c*c-b*b-d*d, 2*(c*d-a*b)],
                  [2*(b*d-a*c), 2*(c*d+a*b), a*a+d*d-b*b-c*c]])

Сроки:

>>> import timeit
>>> 
>>> setup = """
... import numpy as np
... import numpy.random as nr
... 
... from rotation_matrix_test import rotation_matrix_weave
... from rotation_matrix_test import rotation_matrix_numpy
... 
... mat1 = np.eye(3,3)
... theta = nr.random()
... axis = nr.random(3)
... """
>>> 
>>> timeit.repeat("rotation_matrix_weave(axis, theta, mat1)", setup=setup, number=100000)
[0.36641597747802734, 0.34883809089660645, 0.3459300994873047]
>>> timeit.repeat("rotation_matrix_numpy(axis, theta)", setup=setup, number=100000)
[7.180983066558838, 7.172032117843628, 7.180462837219238]

Answer 4

Вот элегантный метод, использующий кватернионы, которые являются невероятно быстрыми; Я могу рассчитать 10 миллионов оборотов в секунду с соответствующим векторным массивом. Он основан на расширении кватерниона для найденной numy здесь .

Теория кватернионов: Кватернион - это число с одним действительным и 3 мнимыми измерениями, обычно записываемое как q = w + xi + yj + zk, где «i», «j», «k» - мнимые измерения. Подобно тому, как единичное комплексное число «c» может представлять все 2d вращения на c=exp(i * theta), так и единичный кватернион «q» может представлять все 3d вращения на q=exp(p), где «p» представляет собой чисто воображаемый кватернион, заданный вашей осью и углом .

Мы начнем с преобразования вашей оси и угла в кватернион, мнимые размеры которого определяются вашей осью вращения, а величина равна половине угла поворота в радианах. Векторы из 4 элементов (w, x, y, z) строятся следующим образом:

import numpy as np
import quaternion as quat

v = [3,5,0]
axis = [4,4,1]
theta = 1.2 #radian

vector = np.array([0.] + v)
rot_axis = np.array([0.] + axis)
axis_angle = (theta*0.5) * rot_axis/np.linalg.norm(rot_axis)

Сначала строится массив из 4 элементов с вещественной составляющей w = 0 как для вектора, который должен вращаться vector, так и для оси вращения rot_axis. Представление угла оси затем строится путем нормализации и умножения на половину желаемого угла theta. См. здесь , почему требуется половина угла.

Теперь создайте кватернионы v и qlog, используя библиотеку, и получите кватернион вращения единицы q, взяв экспоненту.

vec = quat.quaternion(*v)
qlog = quat.quaternion(*axis_angle)
q = np.exp(qlog)

Наконец, вращение вектора вычисляется с помощью следующей операции.

v_prime = q * vec * np.conjugate(q)

print(v_prime) # quaternion(0.0, 2.7491163, 4.7718093, 1.9162971)

Теперь просто отбросьте реальный элемент, и у вас есть повернутый вектор!

v_prime_vec = v_prime.imag # [2.74911638 4.77180932 1.91629719] as a numpy array

Обратите внимание, что этот метод особенно эффективен, если вам нужно повернуть вектор на множество последовательных вращений, поскольку произведение кватернионов можно просто вычислить как q = q1 * q2 * q3 * q4 * ... * qn, а затем вектор вращается только на 'q' в самом конце, используя v '= q * v * con (q).

Этот метод обеспечивает плавное преобразование между оператором трехмерного поворота оси <---> просто с помощью функций exp и log (да log(q) просто возвращает представление угла оси!). Подробнее о том, как работает умножение кватернионов и т. Д., См. здесь

Answer 5

Взгляните на http://vpython.org/contents/docs/visual/VisualIntro.html.

Он предоставляет класс vector, который имеет метод A.rotate(theta,B).Он также предоставляет вспомогательную функцию rotate(A,theta,B), если вы не хотите вызывать метод для A.

http://vpython.org/contents/docs/visual/vector.html

Answer 6

Я сделал довольно полную библиотеку трехмерной математики для Python {2,3}. Он все еще не использует Cython, но сильно зависит от эффективности numpy. Вы можете найти его здесь с помощью пункта:

python[3] -m pip install math3d

Или посмотрите на мой gitweb http://git.automatics.dyndns.dk/?p=pymath3d.git, а теперь и на github: https://github.com/mortlind/pymath3d.

После установки в python вы можете создать объект ориентации, который может вращать векторы или быть частью объектов преобразования. Например. следующий фрагмент кода создает ориентацию, которая представляет вращение на 1 рад вокруг оси [1,2,3], применяет его к вектору [4,5,6] и печатает результат:

import math3d as m3d
r = m3d.Orientation.new_axis_angle([1,2,3], 1)
v = m3d.Vector(4,5,6)
print(r * v)

Вывод будет

<Vector: (2.53727, 6.15234, 5.71935)>

Насколько я могу судить, это примерно в четыре раза эффективнее, чем на однослойном с использованием scipy, опубликованном Б. М. выше. Однако, это требует установки моего пакета math3d.

Answer 7

Отказ от ответственности: я являюсь автором этого пакета

Хотя специальные классы для вращений могут быть удобными, в некоторых случаях нужны матрицы вращения (например, для работы с другими библиотеками, такими как функции affine_transform в scipy).Чтобы каждый не реализовывал свои собственные маленькие функции генерации матриц, существует крошечный чистый пакет python, который делает только удобные функции генерации матриц вращения.Пакет находится на github ( mgen ) и может быть установлен через pip:

pip install mgen

Пример использования, скопированный из readme:

import numpy as np
np.set_printoptions(suppress=True)

from mgen import rotation_around_axis
from mgen import rotation_from_angles
from mgen import rotation_around_x

matrix = rotation_from_angles([np.pi/2, 0, 0], 'XYX')
matrix.dot([0, 1, 0])
# array([0., 0., 1.])

matrix = rotation_around_axis([1, 0, 0], np.pi/2)
matrix.dot([0, 1, 0])
# array([0., 0., 1.])

matrix = rotation_around_x(np.pi/2)
matrix.dot([0, 1, 0])
# array([0., 0., 1.])

Обратите внимание, что матрицыявляются просто обычными массивами numpy, поэтому при использовании этого пакета новые структуры данных не вводятся.

Answer 8

Это также можно решить с помощью теории кватернионов:

def angle_axis_quat(theta, axis):
    """
    Given an angle and an axis, it returns a quaternion.
    """
    axis = np.array(axis) / np.linalg.norm(axis)
    return np.append([np.cos(theta/2)],np.sin(theta/2) * axis)

def mult_quat(q1, q2):
    """
    Quaternion multiplication.
    """
    q3 = np.copy(q1)
    q3[0] = q1[0]*q2[0] - q1[1]*q2[1] - q1[2]*q2[2] - q1[3]*q2[3]
    q3[1] = q1[0]*q2[1] + q1[1]*q2[0] + q1[2]*q2[3] - q1[3]*q2[2]
    q3[2] = q1[0]*q2[2] - q1[1]*q2[3] + q1[2]*q2[0] + q1[3]*q2[1]
    q3[3] = q1[0]*q2[3] + q1[1]*q2[2] - q1[2]*q2[1] + q1[3]*q2[0]
    return q3

def rotate_quat(quat, vect):
    """
    Rotate a vector with the rotation defined by a quaternion.
    """
    # Transfrom vect into an quaternion 
    vect = np.append([0],vect)
    # Normalize it
    norm_vect = np.linalg.norm(vect)
    vect = vect/norm_vect
    # Computes the conjugate of quat
    quat_ = np.append(quat[0],-quat[1:])
    # The result is given by: quat * vect * quat_
    res = mult_quat(quat, mult_quat(vect,quat_)) * norm_vect
    return res[1:]

v = [3, 5, 0]
axis = [4, 4, 1]
theta = 1.2 

print(rotate_quat(angle_axis_quat(theta, axis), v))
# [2.74911638 4.77180932 1.91629719]

Answer 9

Использование pyquaternion чрезвычайно просто;чтобы установить его (пока в Python), запустите в консоли:

import pip;
pip.main(['install','pyquaternion'])

После установки:

  from pyquaternion import Quaternion
  v = [3,5,0]
  axis = [4,4,1]
  theta = 1.2 #radian
  rotated_v = Quaternion(axis=axis,angle=theta).rotate(v)

Answer 10

Используйте scipy's Rotation.from_rotvec(). Аргументом является вектор вращения (единичный вектор), умноженный на угол поворота в радах.

from scipy.spatial.transform import Rotation
from numpy.linalg import norm


v = [3, 5, 0]
axis = [4, 4, 1]
theta = 1.2

axis = axis / norm(axis)  # normalize the rotation vector first
rot = Rotation.from_rotvec(theta * axis)

new_v = rot.apply(v)  
print(new_v)    # results in [2.74911638 4.77180932 1.91629719]

Существует еще несколько способов использования Rotation в зависимости от имеющихся у вас данных о ротации:

• from_quat Инициализируется из кватернионов.

• from_dcm Инициализируется из направляющих косинусных матриц.

• from_euler Инициализировано с углов Эйлера.

---

_{Примечание не по теме: Однострочный код не обязательно лучший код, как подразумевают некоторые пользователи.}