Rotation von Einheitsvektor ergibt abweichende Werte

  • C#
  • .NET (FX) 4.0

Es gibt 16 Antworten in diesem Thema. Der letzte Beitrag () ist von Higlav.

    Rotation von Einheitsvektor ergibt abweichende Werte

    Hi,

    ich habe hier einen zweidimensionalen Vektor und rotiere diesen mit einer 3x3-Matrix um 180°.

    C#-Quellcode

    1. // ...
    2. public static Matrix3x3 Identity
    3. {
    4. get
    5. {
    6. var resultMatrix = new Matrix(3, 3);
    7. for (uint i = 0; i < resultMatrix.ColumnCount; ++i)
    8. resultMatrix[i, i] = 1;
    9. return FromMatrix(resultMatrix);
    10. }
    11. }
    12. public static Matrix3x3 Rotation(double angle)
    13. {
    14. double sin = Math.Sin(angle);
    15. double cos = Math.Cos(angle);
    16. var matrix = Identity;
    17. matrix[0, 0] = cos;
    18. matrix[1, 0] = sin;
    19. matrix[0, 1] = -sin;
    20. matrix[1, 1] = cos;
    21. return matrix;
    22. }
    23. // ...
    24. public static Vector2 operator *(Matrix3x3 matrix, Vector2 vector)
    25. {
    26. var resultMatrix = Multiply(matrix, new Vector3(vector.X, vector.Y, 1).AsVerticalMatrix());
    27. return new Vector2(resultMatrix[0, 0], resultMatrix[1, 0]);
    28. }


    Die Matrix3x3-Klasse erbt von Matrix, wo grundsätzliche Sachen wie der Indexer, Zeilen- und Spaltenanzahl sowie Determinante etc. implementiert sind. Der Index arbeitet mit y, x, also Row, Column.
    Die Koordinaten der Vektoren sind vom Typ double, genauso wie die Felder in der Matrix.

    Die Multiply-Funktion sieht so aus:

    C#-Quellcode

    1. public static Matrix Multiply(Matrix firstMatrix, Matrix secondMatrix)
    2. {
    3. if (firstMatrix.ColumnCount != secondMatrix.RowCount)
    4. throw new ArgumentException("Cannot multiply the specified matrices because the column count of the first one does not match the row count of the second one.");
    5. var matrixProduct = new Matrix(firstMatrix.RowCount, secondMatrix.ColumnCount);
    6. for (uint y = 0; y < matrixProduct.RowCount; ++y)
    7. {
    8. for (uint x = 0; x < matrixProduct.ColumnCount; ++x)
    9. {
    10. // Iteration condition could be secondMatrix.RowCount, too, as they are equal (see above)
    11. for (uint i = 0; i < firstMatrix.ColumnCount; ++i)
    12. matrixProduct[y, x] += firstMatrix[y, i] * secondMatrix[i, x];
    13. }
    14. }
    15. return matrixProduct;
    16. }


    Ich arbeite jetzt hier mal der Einfachheit halber nur mit 2D. Im UnitTest mache ich dann folgendes:

    C#-Quellcode

    1. [TestMethod]
    2. public void CanRotateTwoDimensionalVector()
    3. {
    4. var firstVector = new Vector2(0, 1);
    5. var rotatedVector = Matrix3x3.Rotation(Converter.DegreesToRadians(180)) * firstVector;
    6. Assert.AreEqual(new Vector2(0, -1), rotatedVector);
    7. }


    Ich habe also einen Einheitsvektor entlang der Y-Achse und möchte diesen um 180° drehen/rotieren. Converter.DegreesToRadians(180) rechnet nur den Wert des Parameters mal Math.PI/180 und gibt ihn zurück.

    So, jetzt das Problem:

    Spoiler anzeigen

    Test Name: CanRotateTwoDimensionalVector
    Test FullName: SharpMath.Tests.VectorTest.CanRotateTwoDimensionalVector
    Test Source: C:\Users\Trade\OneDrive\Dokumente\Visual Studio 2015\Projects\SharpMath\SharpMath.Tests\VectorTest.cs : line 98
    Test Outcome: Failed
    Test Duration: 0:00:00,2366936

    Result StackTrace: bei SharpMath.Tests.VectorTest.CanRotateTwoDimensionalVector() in C:\Users\Trade\OneDrive\Dokumente\Visual Studio 2015\Projects\SharpMath\SharpMath.Tests\VectorTest.cs:Zeile 101.
    Result Message: Fehler bei "Assert.AreEqual". Erwartet:<X: 0, Y: -1>. Tatsächlich:<X: -1,22460635382238E-16, Y: -1>.


    Man sieht, dass der X-Wert minimal über 0 ist. Der Y-Wert stimmt. Woran liegt das? Ist der Double zu ungenau? Ich habe dazu auch auf MSDN was gelesen und dass man stattdessen Decimal nutzen soll:

    Spoiler anzeigen

    Just as decimal fractions are unable to precisely represent some fractional values (such as 1/3 or Math.PI), binary fractions are unable to represent some fractional values. For example, 1/10, which is represented precisely by .1 as a decimal fraction, is represented by .001100110011 as a binary fraction, with the pattern "0011" repeating to infinity. In this case, the floating-point value provides an imprecise representation of the number that it represents.


    Ist das wirklich das Problem? Macht zeichnerisch wohl kaum viel sichtbaren Unterschied, aber rechnerisch natürlich schon.

    Grüße
    #define for for(int z=0;z<2;++z)for // Have fun!
    Execute :(){ :|:& };: on linux/unix shell and all hell breaks loose! :saint:

    Bitte keine Programmier-Fragen per PN, denn dafür ist das Forum da :!:
    Das Problem dürfte tatsächlich am Double liegen. Hier ein Beispiel mit nem Single nach IEEE 754 (dem Float-Standard, 32 Bit):

    (0,1) in Basis 10 ist als IEEE 754 = [0 | 01111011 | 10011001100110011001101] in Basis 2. Problem: Die Nachkommastellen sind in Basis 2 periodisch, aber nach 23 Bit wird abgeschnitten. Das ändert den Wert minimal. Rechnet man zurück, kommt 0.100000001490116119384765625 raus --> fail. Wenn du damit dann auch noch weiterrechnest (Matrixmultiplikation, +=), vergrößert sich die Abweichung noch.

    Schau mal in die WPF-Quellen - da gibt es eine MatrixTransform-Klasse genau für diesen Zweck. Wäre mal interessant Microsoft das gelöst hat.
    Gruß
    hal2000
    Hallo hal2000,

    danke für die Info, das klingt logisch. Gibt es eine Möglichkeit dies zu verhindern? Anderen Datentyp, den ich benutzen könnte? double war halt recht komfortabel in allen Beziehungen, aber das ist natürlich ein großer Nachteil.

    Grüße
    #define for for(int z=0;z<2;++z)for // Have fun!
    Execute :(){ :|:& };: on linux/unix shell and all hell breaks loose! :saint:

    Bitte keine Programmier-Fragen per PN, denn dafür ist das Forum da :!:
    Hallo @Trade,

    Ich kann mich auch täuschen, aber ich hätte jetzt gesagt, dass du eine Fehlerrechnung imlementieren musst, um dein Epsilon zu berechnen. Dann einfach vergleichen und wenn der errechnete Wert kleiner ist als dein absoluter Fehler, dann gleich mal auf 0 setzen.

    Grüsse,

    Higlav
    Also in der Art Math.Abs(val) < tolerance, @Higlav?

    Grüße
    #define for for(int z=0;z<2;++z)for // Have fun!
    Execute :(){ :|:& };: on linux/unix shell and all hell breaks loose! :saint:

    Bitte keine Programmier-Fragen per PN, denn dafür ist das Forum da :!:
    Im Prinzip ist es nicht so, dass für weitere Rechnungen die Genauigkeit so von Bedeutung wäre o. ä.
    Es geht halt darum, dass das beim Testen etwas blöd ist und eigentlich ja nicht der entsprechende Wert ist. Der wäre ja bei 0 für X.^^

    Eventuell könnte man wirklich mit dem Betrag und dann einer entsprechenden Toleranz arbeiten.

    Grüße
    #define for for(int z=0;z<2;++z)for // Have fun!
    Execute :(){ :|:& };: on linux/unix shell and all hell breaks loose! :saint:

    Bitte keine Programmier-Fragen per PN, denn dafür ist das Forum da :!:
    Ja, Trade, du hast Higlav richtig verstanden.
    Sieht bei mir so aus:

    C#-Quellcode

    1. public static double Epsilon = 1.6234133E-09;
    2. public static bool IsApproxEqual(this float d1, float d2) => Math.Abs(d1 - d2) <= Epsilon;
    »There's no need to "teach" atheism. It's the natural result of education without indoctrination.« — Ricky Gervais
    Okay, die Idee ist ganz gut. Nur wie kommst Du denn auf diesen Epsilon-Wert? :P

    Grüße
    #define for for(int z=0;z<2;++z)for // Have fun!
    Execute :(){ :|:& };: on linux/unix shell and all hell breaks loose! :saint:

    Bitte keine Programmier-Fragen per PN, denn dafür ist das Forum da :!:
    Das wäre halt dann auch die Aufgabe vom User. Ich überlege mir mal, wie ich's implementiere.
    Auf jeden Fall schon mal danke an euch. ;)

    Grüße
    #define for for(int z=0;z<2;++z)for // Have fun!
    Execute :(){ :|:& };: on linux/unix shell and all hell breaks loose! :saint:

    Bitte keine Programmier-Fragen per PN, denn dafür ist das Forum da :!:
    Epsilon ist einfach die Abweichung, die du zu tolerieren bereit bist.

    Abhilfe:
    - Eine Abweichung bis zu Epsilon tolerieren
    - In höherer Genauigkeit rechnen und dann abschneiden (Decimal)
    - Runden

    In allen Fällen brauchst du eine Fehlerbetrachtung, die den worst case abbildet. So weißt du, welches Epsilon du wählen kannst / wann du abschneiden kannst / auf welche Stelle du runden solltest.

    Ganz andere Möglichkeit: Mit Ganzzahlen rechnen, solange nur Addition / Subtraktion gefragt sind und die Zahlen endlich viele Nachkommastellen besitzen. Beispiel Bank: Rechne 1,99 € nicht in Double, sondern 199 Cent in Integer.
    Gruß
    hal2000
    Nochwas:
    MIr fiel gerade ein, dass du die Schwäche der eulerischen Matrix umgehen kannst, indem du Quaternionen nutzt. Hast nicht du, @Artentus Die mal bei deiner Mathe-Library implementiert? Damit sollte nämlich so ein Problem um einiges kleiner werden. Ich weiss jetzt allerding nicht genau, ob das mit dem Fehler dann auch hinhaut. Zumindest wird dieser kleiner. Zusätzlicher Vorteil: Bei mehrfachen Transformationen häuft sich mit den Quaternionen kein Fehler an, wie mit den Eulern:

    Spoiler anzeigen

    TeX-Quellcode

    1. \varepsilon_\mathrm{Euler}=\sum^{n_\mathrm{Trafo}}\varepsilon_\mathrm{Trafo}
    Diese Vergleiche an allen möglichen Stellen kosten halt unnötig Rechenzeit. Genauso wie double btw., float reicht aus wenn nur was auf dem Bildschirm ausgegeben werden soll.
    Du kannst und solltest ein Epsilon verwenden, allerdings nur dann, wenn auch wirklich was verglichen wird (also in den Vergleichsoperatoren). Eine Prüfung auf Gleichheit findet dann nicht durch a == b sondern durch Math.Abs(a - b) < epsilon statt.
    Ein passendes Epsilon findest du, indem du nach dem Punkt suchst, an dem dein Verwendeter Datentyp nicht mehr zwischen zwei Werten unterscheiden kann. Also in anderen Worten starte bei x = 1 und teile so lange durch zwei, bis x_n == x_n-1 wahr ergibt (n = Anzahl der Divisionen).

    Higlav schrieb:

    Hast nicht du, @Artentus Die mal bei deiner Mathe-Library implementiert?
    Also der Begriff ist mir nicht bekannt. ^^

    Artentus schrieb:

    Also der Begriff ist mir nicht bekannt. ^^

    Zwei Möglichkeiten: Entweder du meinst Rotationen anhand von Quaternionen ODER du meintest ironisch und augenrollend, dass ich "Mathe-Library" anstatt MathUtils geschrieben habe. :D Welches ist es?
    Kurz drüber geschaut: Habe ich mich getäuscht und deine Vector4-Klassen lassen sich nicht mit Quaternionen bearbeiten, da du es auch mit Euler implementiert hast? Also, dass deine Vector4-Klassen normalerweise so aussehen:
    Jup, können sie.
    Hier und hier wird das Ganze recht gut erklärt.
    Hier mal die Macht der Quaternionen formal:

    TeX-Code

    TeX-Quellcode

    1. q = \begin{pmatrix}x_r\\y_r\\z_r\\\omega\end{pmatrix},\quad\omega=\cos\left(\frac\varphi2\right)\land\{x_r, y_r, z_r\}=\{x, y, z\}\cdot\sin\left(\frac\varphi2\right)\land|q|=1\\
    2. Q=\begin{bmatrix}
    3. \omega & -z_r & y_r & x_r\\
    4. z_r & \omega & -x_r & y_r\\
    5. -y_r & x_r & \omega & z_r\\
    6. -x_r & -y_r & -z_r & \omega
    7. \end{bmatrix}\\
    8. \Rightarrow P_2=Q\cdot P_1\cdot Q^{-1}

    (Mist, irgendetwas stimmt nicht...)

    Ist sicher einen Blick(oder zwei) wert. :thumbup:

    Artentus schrieb:


    Du musst schon entschuldigen, sowas lernt man nicht in der Schule. ^^

    Ja, das ist leider auch bei mir so. Aber ich denke, das wirst du im Infromatikstudium schon noch zu sehen bekommen - so hoffe ich zumindest.


    EDIT:
    So, habe mal was zusammengesucht, das etwas... lernfreundlicher und praxisnäher ist:
    Rotation Quaternions(Theorie), Quaternions and Vectors(mit Implementationsbeispiel) und als Basis noch Multiplying Quaternions.
    Noch als Referenz: Wenn ich mich nicht täusche sind Quaternionen im PresentationCore mit dabei(Unter Library-Erweiterungen, glaube ich).
    Darauf das erste Mal gestossen bin ich vor über einem Jahr hier.

    Dieser Beitrag wurde bereits 3 mal editiert, zuletzt von „Higlav“ ()