Подтвердить что ты не робот

Вращающаяся система координат через кватернион

У нас есть gazillion пространственные координаты (x, y и z), представляющие атомы в 3d-пространстве, и я строю функцию, которая переведет эти точки в новую систему координат. Смещение координат к произвольному происхождению прост, но я не могу обернуть голову вокруг следующего шага: вычисления трехмерного вращения точки. Другими словами, я пытаюсь перевести точки из (x, y, z) в (x ', y', z '), где x', y 'и z' находятся в терминах я ', j' и k ', новые векторы оси, которые я создаю с помощью euplid python module.

Я думаю, что все, что мне нужно, это евклидский кватернион, чтобы сделать это, т.е.

>>> q * Vector3(x, y, z)
Vector3(x', y', z')

но для этого я считаю, что мне нужен вектор оси вращения и угол поворота. Но я не знаю, как вычислить их из я ', j' и k '. Это похоже на простую процедуру кодирования с нуля, но я подозреваю, что что-то вроде этого требует, чтобы линейная алгебра выяснялась сама по себе. Большое спасибо за толкание в правильном направлении.

4b9b3361

ОТВЕТЫ

Ответ 1

Использование кватернионов для представления вращения не является трудным с алгебраической точки зрения. Лично мне нелегко рассуждать о кватернионах, но формулы, связанные с их использованием для вращения, довольно просты. Я приведу здесь базовый набор эталонных функций. 1 (см. Также этот прекрасный ответ hosolmaz, в котором он объединяет их вместе создайте удобный класс кватернионов.)

Вы можете думать о кватернионах (для наших целей) как скаляр плюс трехмерный вектор - абстрактно, w + xi + yj + zk, здесь представлен простой набор (w, x, y, z). Пространство трехмерных вращений представлено в полном объеме подпространством кватернионов - пространством единичных кватернионов, поэтому вы хотите убедиться, что ваши кватернионы нормализованы. Вы можете сделать это так, как если бы вы нормализовали любой 4-векторный (то есть величина должна быть близка к 1, если это не так, уменьшите значения по величине):

def normalize(v, tolerance=0.00001):
    mag2 = sum(n * n for n in v)
    if abs(mag2 - 1.0) > tolerance:
        mag = sqrt(mag2)
        v = tuple(n / mag for n in v)
    return v

Обратите внимание, что для простоты следующие функции предполагают, что значения кватернионов уже нормализованы. На практике вам нужно время от времени перенормировать их, но лучший способ справиться с этим будет зависеть от проблемной области. Эти функции предоставляют только самые основы, только для справочных целей.

Каждое вращение представлено единичным кватернионом, а конкатенации поворотов соответствуют умножениям единичных кватернионов. Формула 2 для этого выглядит так:

def q_mult(q1, q2):
    w1, x1, y1, z1 = q1
    w2, x2, y2, z2 = q2
    w = w1 * w2 - x1 * x2 - y1 * y2 - z1 * z2
    x = w1 * x2 + x1 * w2 + y1 * z2 - z1 * y2
    y = w1 * y2 + y1 * w2 + z1 * x2 - x1 * z2
    z = w1 * z2 + z1 * w2 + x1 * y2 - y1 * x2
    return w, x, y, z

Чтобы повернуть вектор по кватерниону, вам понадобится также кватернионное сопряжение. Это легко:

def q_conjugate(q):
    w, x, y, z = q
    return (w, -x, -y, -z)

Теперь умножение кватернионного вектора так же просто, как преобразование вектора в кватернион (путем установки w = 0 и оставления x, y и z того же самого), а затем умножения q * v * q_conjugate(q):

def qv_mult(q1, v1):
    q2 = (0.0,) + v1
    return q_mult(q_mult(q1, q2), q_conjugate(q1))[1:]

Наконец, вам нужно знать, как преобразовать из осевых поворотов в кватернионы. Также легко! Имеет смысл "санировать" вход и выход здесь, вызывая normalize.

def axisangle_to_q(v, theta):
    v = normalize(v)
    x, y, z = v
    theta /= 2
    w = cos(theta)
    x = x * sin(theta)
    y = y * sin(theta)
    z = z * sin(theta)
    return w, x, y, z

И обратно:

def q_to_axisangle(q):
    w, v = q[0], q[1:]
    theta = acos(w) * 2.0
    return normalize(v), theta

Вот пример быстрого использования. Последовательность поворотов на 90 градусов вокруг осей x, y и z вернет вектор по оси y в исходное положение. Этот код выполняет эти вращения:

x_axis_unit = (1, 0, 0)
y_axis_unit = (0, 1, 0)
z_axis_unit = (0, 0, 1)
r1 = axisangle_to_q(x_axis_unit, numpy.pi / 2)
r2 = axisangle_to_q(y_axis_unit, numpy.pi / 2)
r3 = axisangle_to_q(z_axis_unit, numpy.pi / 2)

v = qv_mult(r1, y_axis_unit)
v = qv_mult(r2, v)
v = qv_mult(r3, v)

print v
# output: (0.0, 1.0, 2.220446049250313e-16)

Имейте в виду, что эта последовательность вращений не вернет все векторы в одно и то же положение; например, для вектора на оси х, он будет соответствовать повороту на 90 градусов вокруг оси y. (Обратите внимание на правильное правило: положительное вращение вокруг оси y толкает вектор по оси x в отрицательную область z.)

v = qv_mult(r1, x_axis_unit)
v = qv_mult(r2, v)
v = qv_mult(r3, v)

print v
# output: (4.930380657631324e-32, 2.220446049250313e-16, -1.0)

Как всегда, пожалуйста, дайте мне знать, если вы найдете здесь какие-либо проблемы.

1. Они адаптированы из учебника OpenGL архивированного здесь.

2. Формула умножения кватернионов выглядит как сумасшедшее гнездо крысы, но вывод прост (если утомительный). Заметим сначала, что ii = jj = kk = -1; то ij = k, jk = i, ki = j; и, наконец, ji = -k, kj = -i, ik = -j. Затем умножьте два кватерниона, распределяя их и перестраивая их по результатам каждого из 16 умножений. Это также помогает проиллюстрировать, почему вы можете использовать кватернионы для представления вращения; последние шесть тождеств следуют правилу правой руки, создавая биекции между вращениями от i до j и вращения вокруг k и т.д.

Ответ 2

Заметим, что инверсия матрицы не является вообще тривиальной! Во-первых, все n (где n - размерность вашего пространства), точки должны быть в общем положении (т.е. никакая отдельная точка не может быть выражена как линейная комбинация остальной части точек [оговорка: это может показаться простым требованием, но в области числовой линейной алгебры это нетривиальный, окончательный децизон, который такой конфигурации действительно существует или нет, в конечном итоге будет основываться на специфическом значении "фактической области" ]).

Также "соответствие" новой и старой точек может быть неточным (и тогда вы должны использовать наилучший возможный аппроксиматор "истинного соответствия", т.е.:). Pseudo inverse (вместо того, чтобы пытаться использовать простой обратный), рекомендуется всегда, когда ваш lib предоставляет его.

Псевдо-обратное имеет то преимущество, что вы сможете использовать больше очков для своей трансформации, следовательно, увеличивая вероятность того, что по крайней мере n точек будут в общем положении.

Вот пример, поворот единицы квадрата 90 град. ccw в 2D (но, очевидно, это определение работает в любом dim), с numpy:

In []: P=  matrix([[0, 0, 1, 1],
                   [0, 1, 1, 0]])
In []: Pn= matrix([[0, -1, -1, 0],
                   [0,  0,  1, 1]])
In []: T= Pn* pinv(P)
In []: (T* P).round()
Out[]:
matrix([[ 0., -1., -1.,  0.],
        [ 0.,  0.,  1.,  1.]])

P.S. numpy также быстро. Преобразование 1 млн. Точек на моем скромном компьютере:

In []: P= matrix(rand(2, 1e6))
In []: %timeit T* P
10 loops, best of 3: 37.7 ms per loop

Ответ 3

Этот вопрос и ответ, данный @senderle, действительно помогли мне с одним из моих проектов. Ответ минимален и охватывает ядро ​​большинства вычислений кватернионов, которые, возможно, потребуется выполнить.

Для моего собственного проекта мне показалось утомительным, что у меня есть отдельные функции для всех операций и импортировать их один за другим каждый раз, когда он мне нужен, поэтому я реализовал объектно-ориентированную версию.

quaternion.py:

import numpy as np
from math import sin, cos, acos, sqrt

def normalize(v, tolerance=0.00001):
    mag2 = sum(n * n for n in v)
    if abs(mag2 - 1.0) > tolerance:
        mag = sqrt(mag2)
        v = tuple(n / mag for n in v)
    return np.array(v)

class Quaternion:

    def from_axisangle(theta, v):
        theta = theta
        v = normalize(v)

        new_quaternion = Quaternion()
        new_quaternion._axisangle_to_q(theta, v)
        return new_quaternion

    def from_value(value):
        new_quaternion = Quaternion()
        new_quaternion._val = value
        return new_quaternion

    def _axisangle_to_q(self, theta, v):
        x = v[0]
        y = v[1]
        z = v[2]

        w = cos(theta/2.)
        x = x * sin(theta/2.)
        y = y * sin(theta/2.)
        z = z * sin(theta/2.)

        self._val = np.array([w, x, y, z])

    def __mul__(self, b):

        if isinstance(b, Quaternion):
            return self._multiply_with_quaternion(b)
        elif isinstance(b, (list, tuple, np.ndarray)):
            if len(b) != 3:
                raise Exception(f"Input vector has invalid length {len(b)}")
            return self._multiply_with_vector(b)
        else:
            raise Exception(f"Multiplication with unknown type {type(b)}")

    def _multiply_with_quaternion(self, q2):
        w1, x1, y1, z1 = self._val
        w2, x2, y2, z2 = q2._val
        w = w1 * w2 - x1 * x2 - y1 * y2 - z1 * z2
        x = w1 * x2 + x1 * w2 + y1 * z2 - z1 * y2
        y = w1 * y2 + y1 * w2 + z1 * x2 - x1 * z2
        z = w1 * z2 + z1 * w2 + x1 * y2 - y1 * x2

        result = Quaternion.from_value(np.array((w, x, y, z)))
        return result

    def _multiply_with_vector(self, v):
        q2 = Quaternion.from_value(np.append((0.0), v))
        return (self * q2 * self.get_conjugate())._val[1:]

    def get_conjugate(self):
        w, x, y, z = self._val
        result = Quaternion.from_value(np.array((w, -x, -y, -z)))
        return result

    def __repr__(self):
        theta, v = self.get_axisangle()
        return f"((%.6f; %.6f, %.6f, %.6f))"%(theta, v[0], v[1], v[2])

    def get_axisangle(self):
        w, v = self._val[0], self._val[1:]
        theta = acos(w) * 2.0

        return theta, normalize(v)

    def tolist(self):
        return self._val.tolist()

    def vector_norm(self):
        w, v = self.get_axisangle()
        return np.linalg.norm(v)

В этой версии можно просто использовать перегруженные операторы для умножения кватерниона-кватерниона и кватерниона-вектора

from quaternion import Quaternion
import numpy as np

x_axis_unit = (1, 0, 0)
y_axis_unit = (0, 1, 0)
z_axis_unit = (0, 0, 1)

r1 = Quaternion.from_axisangle(np.pi / 2, x_axis_unit)
r2 = Quaternion.from_axisangle(np.pi / 2, y_axis_unit)
r3 = Quaternion.from_axisangle(np.pi / 2, z_axis_unit)

# Quaternion - vector multiplication
v = r1 * y_axis_unit
v = r2 * v
v = r3 * v

print(v)

# Quaternion - quaternion multiplication
r_total = r3 * r2 * r1
v = r_total * y_axis_unit

print(v)

Я не собирался реализовывать полноценный модуль кватернионов, поэтому это снова для учебных целей, как в отличном ответе @senderle. Надеюсь, это поможет тем, кто хочет понять и попробовать новые вещи с кватернионами.

Ответ 4

Кватернионы - это старая школа (очень старая школа, она ушла из поля зрения столетие назад), чтобы сделать такую ​​операцию. Стандартный способ сделать это в наши дни - использовать матрицы.

Самый простой способ генерации правильной матрицы - это следующее. Вы выяснили, что такое я ', j' и k '. Предположим, что вы пишете их как столбцы (координата я сверху, j посередине, k внизу). Тогда матрица A для преобразования из (i ', j', k ') координат в (i, j, k) является только матрицей [i' j 'k']. (Это, собственно, квадратная матрица 3x3.) Тогда обратная к этой матрице матрица, которая отображает из координат (i, j, k) координаты (i ', j', k ').