système de coordonnées rotatif via un quaternion

Notez que l'inversion de la matrice n'est pas si triviale du tout! Premièrement, tous les points n (où n est la dimension de votre espace) doivent être en position générale (c'est-à-dire qu'aucun point individuel ne peut être exprimé comme une combinaison linéaire du reste des points [mise en garde: cela peut sembler être une exigence simple, mais dans le domaine de l'algèbre linéaire numérique, c'est non trivial; la décision finale connaissance]).

Aussi la 'correspondance' des nouveaux et anciens points peut ne pas être exacte (et alors vous devriez utiliser le meilleur approximateur possible de la 'vraie correspondance', c'est-à-dire:). Pseudo inverse (au lieu d'essayer d'utiliser l'inverse simple) est recommandé toujours quand votre lib le fournit.

Le pseudo inverse a l'avantage de pouvoir utiliser plus de points pour votre transformation, augmentant ainsi la probabilité qu'au moins n points soient en général position.

Voici un exemple, la rotation de l'unité carré 90 deg. ccw en 2D (mais évidemment cette détermination fonctionne dans n'importe quel dim), avec 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 est aussi rapide. Transformation de 1 million de points dans mon ordinateur modeste:

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

Cette question et la réponse donnée par @senderle m'ont vraiment aidé dans l'un de mes projets. La réponse est minime et couvre le cœur de la plupart des calculs quaternion que l'on pourrait avoir besoin d'effectuer.

Pour mon propre projet, j'ai trouvé fastidieux d'avoir des fonctions séparées pour toutes les opérations et de les importer une par une chaque fois que j'en ai besoin, donc j'ai implémenté une version orientée objet.

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)

Dans cette version, on peut simplement utiliser le surchargé opérateurs pour la multiplication quaternion-quaternion et quaternion-vecteur

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)

Je n'avais pas l'intention de mettre en œuvre un module quaternion à part entière, donc c'est à nouveau à des fins pédagogiques, comme dans la grande réponse de @senderle. J'espère que cela aide à ceux qui veulent comprendre et essayer de nouvelles choses avec les quaternions.