Replaced llRot2Euler function.

The original function suffered from unexpected results due to rounding
errors. An error of only 1 or 2 ulps would cause the code to not detect
a singularity at Y rotation +/- PI/2 and take the non-singularity code
path. The replacement code does not suffer from wildly inaccurate
results at the +/- PI/2 singularity. The check in the code for the
singularity isn't strictly needed, but gives more consistent results
At the singularity, the X and Z rotations add. The if check simply
forces the X rotation to be zero so the entirety of the X+Z rotation is
carried by Z.

Additionally, the test code has been updated to include test cases that
caused the old code to fail. The test algorithm is also updated to
perform a more meaningful test. The original code checked if the values
against expected values. This could fail at +/- PI rotations since a
rotation around an axis by PI causes the identical effect as a rotation
by -PI. The new test code checks that the returned angles can be used
to recreate a quaternion that causes the same rotation.

2 files changed, 92 insertions, 36 deletions

diff --git a/OpenSim/Region/ScriptEngine/Shared/Api/Implementation/LSL_Api.cs b/OpenSim/Region/ScriptEngine/Shared/Api/Implementation/LSL_Api.cs
index 443e7a5..d6316b2 100644
--- a/OpenSim/Region/ScriptEngine/Shared/Api/Implementation/LSL_Api.cs
+++ b/OpenSim/Region/ScriptEngine/Shared/Api/Implementation/LSL_Api.cs
@@ -468,26 +468,21 @@ namespace OpenSim.Region.ScriptEngine.Shared.Api
 
         //Now we start getting into quaternions which means sin/cos, matrices and vectors. ckrinke
 
-        // Old implementation of llRot2Euler. Normalization not required as Atan2 function will
-        // only return values >= -PI (-180 degrees) and <= PI (180 degrees).
-
+        // Using algorithm based off http://www.euclideanspace.com/maths/geometry/rotations/conversions/quaternionToEuler/quat_2_euler_paper_ver2-1.pdf
+        // to avoid issues with singularity and rounding with Y rotation of +/- PI/2
         public LSL_Vector llRot2Euler(LSL_Rotation r)
         {
-            m_host.AddScriptLPS(1);
-            //This implementation is from http://lslwiki.net/lslwiki/wakka.php?wakka=LibraryRotationFunctions. ckrinke
-            LSL_Rotation t = new LSL_Rotation(r.x * r.x, r.y * r.y, r.z * r.z, r.s * r.s);
-            double m = (t.x + t.y + t.z + t.s);
-            if (m == 0) return new LSL_Vector();
-            double n = 2 * (r.y * r.s + r.x * r.z);
-            double p = m * m - n * n;
-            if (p > 0)
-                return new LSL_Vector(Math.Atan2(2.0 * (r.x * r.s - r.y * r.z), (-t.x - t.y + t.z + t.s)),
-                                             Math.Atan2(n, Math.Sqrt(p)),
-                                             Math.Atan2(2.0 * (r.z * r.s - r.x * r.y), (t.x - t.y - t.z + t.s)));
-            else if (n > 0)
-                return new LSL_Vector(0.0, Math.PI * 0.5, Math.Atan2((r.z * r.s + r.x * r.y), 0.5 - t.x - t.z));
-            else
-                return new LSL_Vector(0.0, -Math.PI * 0.5, Math.Atan2((r.z * r.s + r.x * r.y), 0.5 - t.x - t.z));
+            LSL_Vector v = new LSL_Vector(0.0, 0.0, 1.0) * r;   // Z axis unit vector unaffected by Z rotation component of r.
+            double m = LSL_Vector.Mag(v);                       // Just in case v isn't normalized, need magnitude for Asin() operation later.
+            if (m == 0.0) return new LSL_Vector();
+            double x = Math.Atan2(-v.y, v.z);
+            double sin = v.x / m;
+            if (sin < -0.999999 || sin > 0.999999) x = 0.0;     // Force X rotation to 0 at the singularities.
+            double y = Math.Asin(sin);
+            // Rotate X axis unit vector by r and unwind the X and Y rotations leaving only the Z rotation
+            v = new LSL_Vector(1.0, 0.0, 0.0) * ((r * new LSL_Rotation(Math.Sin(-x / 2.0), 0.0, 0.0, Math.Cos(-x / 2.0))) * new LSL_Rotation(0.0, Math.Sin(-y / 2.0), 0.0, Math.Cos(-y / 2.0)));
+            double z = Math.Atan2(v.y, v.x);
+            return new LSL_Vector(x, y, z);
         }
 
         /* From wiki:
diff --git a/OpenSim/Region/ScriptEngine/Shared/Tests/LSL_ApiTest.cs b/OpenSim/Region/ScriptEngine/Shared/Tests/LSL_ApiTest.cs
index 0cbad41..7594691 100644
--- a/OpenSim/Region/ScriptEngine/Shared/Tests/LSL_ApiTest.cs
+++ b/OpenSim/Region/ScriptEngine/Shared/Tests/LSL_ApiTest.cs
@@ -142,30 +142,91 @@ namespace OpenSim.Region.ScriptEngine.Shared.Tests
         public void TestllRot2Euler()
         {
             // 180, 90 and zero degree rotations.
-            CheckllRot2Euler(new LSL_Types.Quaternion(1.0f, 0.0f, 0.0f, 0.0f), new LSL_Types.Vector3(Math.PI, 0.0f, 0.0f));
-            CheckllRot2Euler(new LSL_Types.Quaternion(0.0f, 1.0f, 0.0f, 0.0f), new LSL_Types.Vector3(Math.PI, 0.0f, Math.PI));
-            CheckllRot2Euler(new LSL_Types.Quaternion(0.0f, 0.0f, 1.0f, 0.0f), new LSL_Types.Vector3(0.0f, 0.0f, Math.PI));
-            CheckllRot2Euler(new LSL_Types.Quaternion(0.0f, 0.0f, 0.0f, 1.0f), new LSL_Types.Vector3(0.0f, 0.0f, 0.0f));
-            CheckllRot2Euler(new LSL_Types.Quaternion(-0.5f, -0.5f, 0.5f, 0.5f), new LSL_Types.Vector3(0, -Math.PI / 2.0f, Math.PI / 2.0f));
-            CheckllRot2Euler(new LSL_Types.Quaternion(-0.707107f, 0.0f, 0.0f, -0.707107f), new LSL_Types.Vector3(Math.PI / 2.0f, 0.0f, 0.0f));
+            CheckllRot2Euler(new LSL_Types.Quaternion(0.0f, 0.0f, 0.0f, 1.0f));
+            CheckllRot2Euler(new LSL_Types.Quaternion(0.0f, 0.0f, 0.707107f, 0.707107f));
+            CheckllRot2Euler(new LSL_Types.Quaternion(0.0f, 0.0f, 1.0f, 0.0f));
+            CheckllRot2Euler(new LSL_Types.Quaternion(0.0f, 0.0f, 0.707107f, -0.707107f));
+            CheckllRot2Euler(new LSL_Types.Quaternion(0.707107f, 0.0f, 0.0f, 0.707107f));
+            CheckllRot2Euler(new LSL_Types.Quaternion(0.5f, -0.5f, 0.5f, 0.5f));
+            CheckllRot2Euler(new LSL_Types.Quaternion(0.0f, -0.707107f, 0.707107f, 0.0f));
+            CheckllRot2Euler(new LSL_Types.Quaternion(-0.5f, -0.5f, 0.5f, -0.5f));
+            CheckllRot2Euler(new LSL_Types.Quaternion(1.0f, 0.0f, 0.0f, 0.0f));
+            CheckllRot2Euler(new LSL_Types.Quaternion(0.707107f, -0.707107f, 0.0f, 0.0f));
+            CheckllRot2Euler(new LSL_Types.Quaternion(0.0f, -1.0f, 0.0f, 0.0f));
+            CheckllRot2Euler(new LSL_Types.Quaternion(-0.707107f, -0.707107f, 0.0f, 0.0f));
+            CheckllRot2Euler(new LSL_Types.Quaternion(0.707107f, 0.0f, 0.0f, -0.707107f));
+            CheckllRot2Euler(new LSL_Types.Quaternion(0.5f, -0.5f, -0.5f, -0.5f));
+            CheckllRot2Euler(new LSL_Types.Quaternion(0.0f, -0.707107f, -0.707107f, 0.0f));
+            CheckllRot2Euler(new LSL_Types.Quaternion(-0.5f, -0.5f, -0.5f, 0.5f));
+            CheckllRot2Euler(new LSL_Types.Quaternion(0.0f, -0.707107f, 0.0f, 0.707107f));
+            CheckllRot2Euler(new LSL_Types.Quaternion(-0.5f, -0.5f, 0.5f, 0.5f));
+            CheckllRot2Euler(new LSL_Types.Quaternion(-0.707107f, 0.0f, 0.707107f, 0.0f));
+            CheckllRot2Euler(new LSL_Types.Quaternion(-0.5f, 0.5f, 0.5f, -0.5f));
+            CheckllRot2Euler(new LSL_Types.Quaternion(0.0f, -0.707107f, 0.0f, -0.707107f));
+            CheckllRot2Euler(new LSL_Types.Quaternion(-0.5f, -0.5f, -0.5f, -0.5f));
+            CheckllRot2Euler(new LSL_Types.Quaternion(-0.707107f, 0.0f, -0.707107f, 0.0f));
+            CheckllRot2Euler(new LSL_Types.Quaternion(-0.5f, 0.5f, -0.5f, 0.5f));
+
             // A couple of messy rotations.
-            CheckllRot2Euler(new LSL_Types.Quaternion(1.0f, 5.651f, -3.1f, 67.023f), new LSL_Types.Vector3(0.037818f, 0.166447f, -0.095595f));
-            CheckllRot2Euler(new LSL_Types.Quaternion(0.719188f, -0.408934f, -0.363998f, -0.427841f), new LSL_Types.Vector3(-1.954769f, -0.174533f, 1.151917f));
+            CheckllRot2Euler(new LSL_Types.Quaternion(1.0f, 5.651f, -3.1f, 67.023f));
+            CheckllRot2Euler(new LSL_Types.Quaternion(0.719188f, -0.408934f, -0.363998f, -0.427841f));
+
+            // Some deliberately malicious rotations (intended on provoking singularity errors)
+            // The "f" suffexes are deliberately omitted.
+            CheckllRot2Euler(new LSL_Types.Quaternion(0.50001f, 0.50001f, 0.50001f, 0.50001f));
+            // More malice. The "f" suffixes are deliberately omitted.
+            CheckllRot2Euler(new LSL_Types.Quaternion(-0.701055, 0.092296, 0.701055, -0.092296));
+            CheckllRot2Euler(new LSL_Types.Quaternion(-0.183005, -0.683010, 0.183005, 0.683010));
+            CheckllRot2Euler(new LSL_Types.Quaternion(-0.430460, -0.560982, 0.430460, 0.560982));
+            CheckllRot2Euler(new LSL_Types.Quaternion(-0.701066, 0.092301, -0.701066, 0.092301));
+            CheckllRot2Euler(new LSL_Types.Quaternion(-0.183013, -0.683010, 0.183013, 0.683010));
+            CheckllRot2Euler(new LSL_Types.Quaternion(-0.183005, -0.683014, -0.183005, -0.683014));
+            CheckllRot2Euler(new LSL_Types.Quaternion(-0.353556, 0.612375, 0.353556, -0.612375));
+            CheckllRot2Euler(new LSL_Types.Quaternion(0.353554, -0.612385, -0.353554, 0.612385));
+            CheckllRot2Euler(new LSL_Types.Quaternion(-0.560989, 0.430450, 0.560989, -0.430450));
+            CheckllRot2Euler(new LSL_Types.Quaternion(-0.183013, 0.683009, -0.183013, 0.683009));
+            CheckllRot2Euler(new LSL_Types.Quaternion(0.430457, -0.560985, -0.430457, 0.560985));
+            CheckllRot2Euler(new LSL_Types.Quaternion(0.353552, 0.612360, -0.353552, -0.612360));
+            CheckllRot2Euler(new LSL_Types.Quaternion(-0.499991, 0.500003, 0.499991, -0.500003));
+            CheckllRot2Euler(new LSL_Types.Quaternion(-0.353555, -0.612385, -0.353555, -0.612385));
+            CheckllRot2Euler(new LSL_Types.Quaternion(0.701066, -0.092301, -0.701066, 0.092301));
+            CheckllRot2Euler(new LSL_Types.Quaternion(-0.499991, 0.500007, 0.499991, -0.500007));
+            CheckllRot2Euler(new LSL_Types.Quaternion(-0.683002, 0.183016, -0.683002, 0.183016));
+            CheckllRot2Euler(new LSL_Types.Quaternion(0.430458, 0.560982, 0.430458, 0.560982));
+            CheckllRot2Euler(new LSL_Types.Quaternion(0.499991, -0.500003, -0.499991, 0.500003));
+            CheckllRot2Euler(new LSL_Types.Quaternion(-0.183009, 0.683011, -0.183009, 0.683011));
+            CheckllRot2Euler(new LSL_Types.Quaternion(0.560975, -0.430457, 0.560975, -0.430457));
+            CheckllRot2Euler(new LSL_Types.Quaternion(0.701055, 0.092300, 0.701055, 0.092300));
+            CheckllRot2Euler(new LSL_Types.Quaternion(-0.560990, 0.430459, -0.560990, 0.430459));
+            CheckllRot2Euler(new LSL_Types.Quaternion(-0.092302, -0.701059, -0.092302, -0.701059));
         }
 
-        private void CheckllRot2Euler(LSL_Types.Quaternion rot, LSL_Types.Vector3 eulerCheck)
+        // Testing Rot2Euler this way instead of comparing against expected angles because
+        // 1. There are several ways to get to the original Quaternion. For example a rotation
+        //    of PI and -PI will give the same result. But PI and -PI aren't equal.
+        // 2. This method checks to see if the calculated angles from a quaternion can be used
+        //    to create a new quaternion to produce the same rotation.
+        // However, can't compare the newly calculated quaternion against the original because
+        // once again, there are multiple quaternions that give the same result. For instance
+        //  <X, Y, Z, S> == <-X, -Y, -Z, -S>.  Additionally, the magnitude of S can be changed
+        // and will still result in the same rotation if the values for X, Y, Z are also changed
+        // to compensate.
+        // However, if two quaternions represent the same rotation, then multiplying the first
+        // quaternion by the conjugate of the second, will give a third quaternion representing
+        // a zero rotation. This can be tested for by looking at the X, Y, Z values which should
+        // be zero.
+        private void CheckllRot2Euler(LSL_Types.Quaternion rot)
         {
             // Call LSL function to convert quaternion rotaion to euler radians.
             LSL_Types.Vector3 eulerCalc = m_lslApi.llRot2Euler(rot);
-            // Check upper and lower bounds of x, y and z.
-            // This type of check is performed as opposed to comparing for equal numbers, in order to allow slight
-            // differences in accuracy.
-            Assert.Greater(eulerCalc.x, eulerCheck.x - ANGLE_ACCURACY_IN_RADIANS, "TestllRot2Euler X lower bounds check fail");
-            Assert.Less(eulerCalc.x, eulerCheck.x + ANGLE_ACCURACY_IN_RADIANS, "TestllRot2Euler X upper bounds check fail");
-            Assert.Greater(eulerCalc.y, eulerCheck.y - ANGLE_ACCURACY_IN_RADIANS, "TestllRot2Euler Y lower bounds check fail");
-            Assert.Less(eulerCalc.y, eulerCheck.y + ANGLE_ACCURACY_IN_RADIANS, "TestllRot2Euler Y upper bounds check fail");
-            Assert.Greater(eulerCalc.z, eulerCheck.z - ANGLE_ACCURACY_IN_RADIANS, "TestllRot2Euler Z lower bounds check fail");
-            Assert.Less(eulerCalc.z, eulerCheck.z + ANGLE_ACCURACY_IN_RADIANS, "TestllRot2Euler Z upper bounds check fail");
+            // Now use the euler radians to recalculate a new quaternion rotation
+            LSL_Types.Quaternion newRot = m_lslApi.llEuler2Rot(eulerCalc);
+            // Multiple original quaternion by conjugate of quaternion calculated with angles.
+            LSL_Types.Quaternion check = rot * new LSL_Types.Quaternion(-newRot.x, -newRot.y, -newRot.z, newRot.s);
+            
+            Assert.AreEqual(0.0, check.x, VECTOR_COMPONENT_ACCURACY, "TestllRot2Euler X bounds check fail");
+            Assert.AreEqual(0.0, check.y, VECTOR_COMPONENT_ACCURACY, "TestllRot2Euler Y bounds check fail");
+            Assert.AreEqual(0.0, check.z, VECTOR_COMPONENT_ACCURACY, "TestllRot2Euler Z bounds check fail");
         }
 
         [Test]