Quaternion: Unterschied zwischen den Versionen

Version vom 10. März 2009, 20:45 Uhr

Hamiltonsche Quaternionen

Übersicht

Quaternionen bilden ein 4D-Zahlensystem ähnlich dem 2D-Zahlensystem der Komplexen Zahlen, jedoch sind sie bei der Multiplikation nicht kommutativ ( d.h. für Quaternionen q1, q2 gilt nicht immer: q1*q2 = q2*q1 ). Sie werden häufig zur Darstellung und einfachen Berechnung von Isometrien (Drehungen) im 3D-Raum verwendet, wobei sie hier deutlich anschaulicher sind als unitäre Matrizen (Rotationsmatrizen).

Definition

Ein Quaternion q hat folgende, eindeutige Gestalt:

q = a + b*i + c*j + d*k,

wobei a, b, c, d reele Zahlen darstellen und i, j, k Imaginäre Zahlen mit den Eigenschaften:

i*i=j*j=k*k=-1
i*j = k = -j*i
j*k = i = -k*j
k*i = j = -i*k

Die Quaternionen H bilden also einen 4-Dimensionalen Raum in den Komponenten a, b, c, d. Man trennt dabei häufig den 3D-Raum der reinen Quaternionen V der Form v = b*i + c*j + d*k heraus.

Grundlegende Arithmetik

Addition

Die Addition zweier Quaternionen geschieht komponentenweise, also:

q = a + bi+ cj + dk
r = e + fi+ gj + hk
q + r = (a + e) + (b + f)i + (c + g)j + (d + h)k

Multiplikation

Die Multiplikation zweier Quaternionen geschieht wie auf den reelen Zahlen gewohnt mit gedachten Variablen i, j, k, wobei die oben definierten Eigenschaften zum Einsatz kommen:

q*r = (a + bi+ cj + dk)*(e + fi+ gj + hk) =
    = ae + afi + agj + ahk  
      + bie + bifi + bigj + bihk
      + cje + cjfi + cjgj + cjhk
      + dke + dkfi + dkgj + dkhk
    = (a * e - b * f  - c * g - d * h)
      + (a * f + b * e + c * h - d * g)i
      + (a * g - b * h + c * e + d * f)j
      + (a * h + b * g - c * f + d * e)k 

Man beachte im Allgemeinen: q*r <> r*q !

Konjugation

Die Konjugation ist definiert durch:

konj(q) = a - bi - cj - dk,

also die Inversion des Vorzeichens der Imaginären Anteile i, j, k. Es lässt sich leicht zeigen, dass gilt:

konj(konj(q)) = q
konj(q*p) = konj(p)* konj(q)
konj(q + p) = konj(q) + konj(p)

Norm / Betrag

Die Norm ||q|| eines Quaternions Q ist definiert als:

 ||q|| = Sqrt(a*a + b*b + c*c + d*d) = Sqrt(q*konj(q))

Inversion

Mit den obigen Funktionen lässt sich sehr leicht ein Quaternion h zu q berechnen, so dass gilt: h*q=q*h=1, denn:

||q||^2 = q * konj(q)
<=> 
1 = (q * konj(q)) / ||q||2
h := konj(q) / ||q||2
=> 1 = q * h

Und weil offensichtlich ||q|| = ||konj(q)|| gilt, folgt

||q||2 = ||konj(q)||2 = konj(q)*konj(konj(q)) = konj(q)*q

Und damit:

1 = h * q.

h ist also das multiplikativ invers zu q: h = q^-1

Injektion der Reelen und Komplexen Zahlen

Die reellen Zahlen lassen sich injektiv in die Quaternionen abbilden:

x -> x + 0*i + 0*j + 0*k

Entsprechend für die komplexen Zahlen:

x + y*i -> x+ y*i + 0*j + 0*k

Die zugehörigen Rechenregeln bleiben dabei erhalten.

Isometrien auf Quaternionen

Isometrien auf H

Die Isometrien der Quaternionen auf sich selbst (also die Isometrien des 4D-Raumes) lassen sich in Quaternionen beschreiben durch zwei Quaternionen P,Q, mit ||P||=||Q|| = 1 und den Abbildungen

φ(P,Q,.): H -> H,
x -> φ(P,Q,x) := P*x*konj(Q) = P*x*Q-1

ψ(P,Q,.): H -> H
x -> ψ(P,Q,x) := P*konj(x)*konj(Q) = P*konj(x)*Q-1

Dabei gilt für Quaternionen P',Q' mit ||P'||=||Q'||=1:

φ(P,Q,φ(P',Q',.)) = φ(P*P',Q*Q',.)
φ(P,Q,ψ(P',Q',.)) = ψ(P*P',Q*Q',.)

ψ(P,Q,ψ(P',Q',.)) = φ(P*Q',Q*P',.)
ψ(P,Q,φ(P',Q',.)) = ψ(P*Q',Q*P',.),

was man ohne weiteres ausrechnen kann.

Isometrien auf V

Der Fall der Isometrien auf den reinen Quaternionen V (siehe Definition) ist der für die 3D Grafik interessante Fall. Sie stellen die Rotationen auf V dar und das in einer besonders leicht lesbaren Form. Sei also Q eine Quaternion der Norm 1, dann gehört dazu eine Isometrie auf V:

φ(Q,.): V -> V
x -> φ(Q,x) := φ(Q,Q,x) = Q*x*konj(Q) = Q*x*Q-1,

wobei offensichtlich φ(-Q,.) = φ(Q,.), Q und -Q beschreiben also die selbe Drehung. Und entsprechnd der Rechenregeln für φ auf H (siehe oben) gilt mit einem Quaternion Q' mit ||Q'||=1:

φ(Q, φ(Q', x)) = Q*Q'*x*konj(Q')*konj(Q) = (Q*Q')*x*konj(Q*Q') = φ(Q*Q',x)

Man kann Q sehr leicht zerlegen in:

Q = cos(α) + sin(α)*P, 0<=α<=π, P aus V, ||P|| = 1,

wobei dann φ(Q,x) eine Drehung von x (man beachte, dass dabei x aus V ist) um den Winkel 2*α um die Achse P darstellt.

Einige Implementationsdetails

Am besten speichert man ein Quaternion duch speichern der 4 Komponenten in einem Array (Reihenfolge: realteil, i, j, k - Komponente, bzw. a,b,c,d):

    public sealed class Quaternion : ILinearTransformation
    {
        private double[] coords = new double[4];

Ein besonders spannender Konstruktor für Quaternion ist der, der "rotations-Quaternionen" erstellt. Übergeben wird also ein Vektor, der die Drehachse beschreibt, sowie ein passender Winkel:

        public Quaternion(double alpha, Vertex3 rotAxis)
        {
            alpha = (2.0*Math.PI/(360.0*2.0))* alpha;

            coords[0] = 1.0;
            if (rotAxis.Magnitude < double.Epsilon)
                return;
            
            rotAxis.Normalize();
            double c = Math.Cos(alpha), s = Math.Sin(alpha);
            coords[0] = c; coords[1] = rotAxis[0] * s;
            coords[2] = rotAxis[1] * s; coords[3] = rotAxis[2] * s;
        }

Diese Funktion ergibt sich ganz einfach aus der obigen Formel:

Q = cos(α) + sin(α)*P, 0<=α<=π, P aus V, ||P|| = 1,

denn d.h. der Realteil der Drehung muss auf cos(α) gesetzt werden und die Drehachse wird mit sin(α) multipliziert und in den reinen Imaginärteil gesteckt (die Komponenten i,j,k werden gerade zu x,y,z im Koordinatensystem).

Ausserdem sollte man die wesentlichen Operationen implementeiren, etwa:

        public Quaternion Conjugate()
        {
            return new Quaternion(coords[0], -coords[1], -coords[2], -coords[3]);
        }

Die Konjugation dreht also wie oben definiert gerade die Vorzeichen der Imaginärteile um.

Genauso implementiert man die anderen Operationen nach:

 ...
        public static Quaternion Multiply(Quaternion q1, Quaternion q2)
        {
            if (q1 == null) throw new ArgumentNullException("q1");
            if (q2 == null) throw new ArgumentNullException("q2");
            Quaternion result = new Quaternion();
            result.coords[0] = q1.coords[0] * q2.coords[0]
                                - q1.coords[1] * q2.coords[1]
                                - q1.coords[2] * q2.coords[2]
                                - q1.coords[3] * q2.coords[3];
            result.coords[1] = q1.coords[0] * q2.coords[1]
                                + q1.coords[1] * q2.coords[0]
                                + q1.coords[2] * q2.coords[3]
                                - q1.coords[3] * q2.coords[2];
            result.coords[2] = q1.coords[0] * q2.coords[2]
                                - q1.coords[1] * q2.coords[3]
                                + q1.coords[2] * q2.coords[0]
                                + q1.coords[3] * q2.coords[1];
            result.coords[3] = q1.coords[0] * q2.coords[3]
                                + q1.coords[1] * q2.coords[2]
                                - q1.coords[2] * q2.coords[1]
                                + q1.coords[3] * q2.coords[0];
            return result;
        }
 ...
        public Quaternion Inverse()
        {
            double mag = Magnitude;
            mag = mag * mag;
            if (mag == 0.0) return null;
            return Conjugate().Scale(1.0 / mag);
        }

Besonders spannend ist sicherlich auch das Anwenden eines Rotations-Quaternions auf ein Vertex. Diese Funktion entsteht, wenn man die oben genannte Funktion φ explizit aufschreibt und ausrechnet. Es ergibt sich dabei eine wenig spannende, aber sehr ausladende Rechnung, um am Ende alles richtig ausmultipliziert und nach komponente sortiert zu haben (nicht vergessen: man darf dabei nur reele Zahlen miteinander Tauschen, sowie reele Zahlen mit i,j,k. i,j,k aber letztere nie miteinander, sonst stimmen die Vorzeichen nicht! Um ein Produkt aus i,j,k loszuwerden verwendet man die Tabelle, welche unter der Definition der Quaternionen gegeben ist -- die Multiplikation war nicht kommutativ:

        public Vertex3 Apply(Vertex3 v)
        {
            Vertex3 result = new Vertex3();
            double v0 = v[0];
            double v1 = v[1];
            double v2 = v[2];
            double a00 = coords[0] * coords[0];
            double a01 = coords[0] * coords[1];
            double a02 = coords[0] * coords[2];
            double a03 = coords[0] * coords[3];
            double a11 = coords[1] * coords[1];
            double a12 = coords[1] * coords[2];
            double a13 = coords[1] * coords[3];
            double a22 = coords[2] * coords[2];
            double a23 = coords[2] * coords[3];
            double a33 = coords[3] * coords[3];
            result[0] = v0 * (+a00 + a11 - a22 - a33) 
                + 2 * (a12 * v1 + a13 * v2 + a02 * v2 - a03 * v1);
            result[1] = v1 * (+a00 - a11 + a22 - a33) 
                + 2 * (a12 * v0 + a23 * v2 + a03 * v0 - a01 * v2);
            result[2] = v2 * (+a00 - a11 - a22 + a33) 
                + 2 * (a13 * v0 + a23 * v1 - a02 * v0 + a01 * v1);
            return result;
        }