logaritmo built python math

python - built - Sistema de coordenadas rotativo a través de un cuaternión.



python math gdc (3)

Tenemos una gran cantidad de coordenadas espaciales (x, y, z) que representan los átomos en el espacio 3d, y estoy construyendo una función que traducirá estos puntos a un nuevo sistema de coordenadas. Cambiar las coordenadas a un origen arbitrario es simple, pero no puedo envolver mi cabeza en el siguiente paso: cálculos de rotación de puntos 3D. En otras palabras, estoy tratando de traducir los puntos de (x, y, z) a (x '', y'', z ''), donde x'', y ''y z'' están en términos de i '', j'' y k '', los nuevos vectores de eje que estoy haciendo con la ayuda del módulo python de euclides .

Creo que todo lo que necesito es un cuaternión de euclides para hacer esto, es decir,

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

pero para hacer ESO creo que necesito un vector de eje de rotación y un ángulo de rotación. Pero no tengo idea de cómo calcularlos a partir de i '', j'' y k ''. Esto parece un procedimiento simple para codificar desde cero, pero sospecho que algo como esto requiere álgebra lineal para resolverlo por mi cuenta. Muchas gracias por un empujón en la dirección correcta.

Esta pregunta y la respuesta dada por @senderle realmente me ayudaron con uno de mis proyectos. La respuesta es mínima y cubre el núcleo de la mayoría de los cálculos de cuaternión que uno podría necesitar realizar.

Para mi propio proyecto, me pareció tedioso tener funciones separadas para todas las operaciones e importarlas una por una cada vez que necesito una, así que implementé una versión orientada a objetos.

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)

En esta versión, solo se pueden usar los operadores sobrecargados para la multiplicación de cuaterniones-cuaterniones y vectores de cuaterniones

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)

No tenía la intención de implementar un módulo de cuaternión de pleno derecho, así que esto es de nuevo con fines de instrucción, como en la gran respuesta de @ senderle. Espero que esto ayude a aquellos que quieren entender y probar cosas nuevas con cuaterniones.

Tenga en cuenta que la inversión de la matriz no es tan trivial en absoluto! En primer lugar, todos los n (donde n es la dimensión de su espacio) los puntos deben estar en la posición general (es decir, ningún punto individual puede expresarse como una combinación lineal del resto de los puntos [advertencia: esto puede parecer ser un requisito simple de hecho, pero en el ámbito del álgebra lineal numérica, no es trivial; la decisión final, ya sea que exista o no dicha configuración, eventualmente se basará en el conocimiento específico del "dominio real").

Además, la "correspondencia" de los puntos nuevos y antiguos puede no ser exacta (y entonces debe utilizar el mejor aproximador posible de la "verdadera correspondencia", es decir, :). La pseudo inversa (en lugar de tratar de utilizar el inverso plano) se recomienda siempre cuando su lib lo proporciona.

El pseudo inverso tiene la ventaja de que podrá usar más puntos para su transformación, por lo que aumenta la probabilidad de que al menos n puntos estén en la posición general.

Aquí hay un ejemplo, rotación de unidad cuadrada 90 grados. ccw en 2D (pero obviamente esta determinación funciona en cualquier numpy ), con 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.]])

PS numpy también es rápido. Transformación de 1 millón de puntos en mi modesta computadora:

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

Usar cuaterniones para representar la rotación no es difícil desde un punto de vista algebraico. Personalmente, me resulta difícil razonar visualmente sobre los cuaterniones, pero las fórmulas involucradas en su uso para rotaciones son bastante simples. Proporcionaré un conjunto básico de funciones de referencia aquí. 1 (Vea también esta hermosa respuesta de hosolmaz , en la que los agrupa para crear una clase de Quaternion práctica).

Puede pensar en los cuaterniones (para nuestros propósitos) como un escalar más un vector tridimensional: de manera abstracta, w + xi + yj + zk , aquí representado por una tupla simple (w, x, y, z) . El espacio de las rotaciones tridimensionales está representado en su totalidad por un subespacio de los cuaterniones, el espacio de los cuaterniones unitarios , por lo que debe asegurarse de que sus cuaterniones estén normalizados. Puede hacerlo de la misma manera que normalizaría cualquier vector de 4 (es decir, la magnitud debería estar cerca de 1; si no lo está, reducir los valores por la magnitud):

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

Tenga en cuenta que, para simplificar, las siguientes funciones asumen que los valores de cuaternión ya están normalizados . En la práctica, tendrá que volver a normalizarlos de vez en cuando, pero la mejor manera de lidiar con eso dependerá del dominio del problema. Estas funciones proporcionan solo lo más básico, solo para fines de referencia.

Cada rotación está representada por un quaternion unitario, y las concatenaciones de rotaciones corresponden a multiplicaciones de cuaterniones unitarios. La fórmula 2 para esto es la siguiente:

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

Para rotar un vector por un cuaternión, también se necesita el conjugado del cuaternión. Eso es fácil:

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

Ahora la multiplicación de vector de cuaternión es tan simple como convertir un vector en un cuaternión (estableciendo w = 0 y dejando igual a x , y y z ) y luego multiplicando 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:]

Finalmente, necesita saber cómo convertir de rotaciones de ángulo de eje a cuaterniones. ¡También fácil! Tiene sentido "sanear" la entrada y la salida aquí llamando a 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

Y de vuelta

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

Aquí hay un ejemplo de uso rápido. Una secuencia de rotaciones de 90 grados sobre los ejes x, y y z devolverá un vector en el eje y a su posición original. Este código realiza esas rotaciones:

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)

Tenga en cuenta que esta secuencia de rotaciones no devolverá todos los vectores a la misma posición; por ejemplo, para un vector en el eje x, corresponderá a una rotación de 90 grados alrededor del eje y. (Tenga en cuenta la regla de la mano derecha aquí; una rotación positiva sobre el eje y empuja un vector en el eje x en la región z negativa ).

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)

Como siempre, avíseme si encuentra algún problema aquí.

1. Estos son adaptados de un tutorial de OpenGL archivado aquí .

2. La fórmula de multiplicación de cuaternión parece un nido de rata loca, pero la derivación es simple (aunque tediosa). Solo note primero que ii = jj = kk = -1 ; entonces que ij = k , jk = i , ki = j ; y finalmente que ji = -k , kj = -i , ik = -j . Luego multiplique los dos cuaterniones, distribuyendo los términos y reorganizándolos en función de los resultados de cada una de las 16 multiplicaciones. Esto también ayuda a ilustrar por qué puede usar los cuaterniones para representar la rotación; las últimas seis identidades siguen la regla de la mano derecha, creando biyecciones entre las rotaciones de i a j y las rotaciones alrededor de k , y así sucesivamente.