Fix QuaternionToAngleAxis to ensure rotations are between -pi and pi.

Thanks to Guoxuan Zhang for reporting this.

Change-Id: I2831ca3a04d5dc6467849c290461adbe23faaea3

include/ceres/rotation.h

diff --git a/include/ceres/rotation.h b/include/ceres/rotation.h
index 7f05187..0d8a390 100644
--- a/include/ceres/rotation.h
+++ b/include/ceres/rotation.h
@@ -148,9 +148,9 @@
 
 template<typename T>
 inline void AngleAxisToQuaternion(const T* angle_axis, T* quaternion) {
-  const T &a0 = angle_axis[0];
-  const T &a1 = angle_axis[1];
-  const T &a2 = angle_axis[2];
+  const T& a0 = angle_axis[0];
+  const T& a1 = angle_axis[1];
+  const T& a2 = angle_axis[2];
   const T theta_squared = a0 * a0 + a1 * a1 + a2 * a2;
 
   // For points not at the origin, the full conversion is numerically stable.
@@ -177,16 +177,35 @@
 
 template<typename T>
 inline void QuaternionToAngleAxis(const T* quaternion, T* angle_axis) {
-  const T &q1 = quaternion[1];
-  const T &q2 = quaternion[2];
-  const T &q3 = quaternion[3];
-  const T sin_squared = q1 * q1 + q2 * q2 + q3 * q3;
+  const T& q1 = quaternion[1];
+  const T& q2 = quaternion[2];
+  const T& q3 = quaternion[3];
+  const T sin_squared_theta = q1 * q1 + q2 * q2 + q3 * q3;
 
   // For quaternions representing non-zero rotation, the conversion
   // is numerically stable.
-  if (sin_squared > T(0.0)) {
-    const T sin_theta = sqrt(sin_squared);
-    const T k = T(2.0) * atan2(sin_theta, quaternion[0]) / sin_theta;
+  if (sin_squared_theta > T(0.0)) {
+    const T sin_theta = sqrt(sin_squared_theta);
+    const T& cos_theta = quaternion[0];
+
+    // If cos_theta is negative, theta is greater than pi/2, which
+    // means that angle for the angle_axis vector which is 2 * theta
+    // would be greater than pi.
+    //
+    // While this will result in the correct rotation, it does not
+    // result in a normalized angle-axis vector.
+    //
+    // In that case we observe that 2 * theta ~ 2 * theta - 2 * pi,
+    // which is equivalent saying
+    //
+    //   theta - pi = atan(sin(theta - pi), cos(theta - pi))
+    //              = atan(-sin(theta), -cos(theta))
+    //
+    const T two_theta =
+        T(2.0) * ((cos_theta < 0.0)
+                  ? atan2(-sin_theta, -cos_theta)
+                  : atan2(sin_theta, cos_theta));
+    const T k = two_theta / sin_theta;
     angle_axis[0] = q1 * k;
     angle_axis[1] = q2 * k;
     angle_axis[2] = q3 * k;

internal/ceres/rotation_test.cc

diff --git a/internal/ceres/rotation_test.cc b/internal/ceres/rotation_test.cc
index 65421fa..8e40507 100644
--- a/internal/ceres/rotation_test.cc
+++ b/internal/ceres/rotation_test.cc
@@ -87,14 +87,43 @@
     return false;
   }
 
+  // Quaternions are equivalent upto a sign change. So we will compare
+  // both signs before declaring failure.
+  bool near = true;
   for (int i = 0; i < 4; i++) {
     if (fabs(arg[i] - expected[i]) > kTolerance) {
-      *result_listener << "component " << i << " should be " << expected[i];
-      return false;
+      near = false;
+      break;
     }
   }
 
-  return true;
+  if (near) {
+    return true;
+  }
+
+  near = true;
+  for (int i = 0; i < 4; i++) {
+    if (fabs(arg[i] + expected[i]) > kTolerance) {
+      near = false;
+      break;
+    }
+  }
+
+  if (near) {
+    return true;
+  }
+
+  *result_listener << "expected : "
+                   << expected[0] << " "
+                   << expected[1] << " "
+                   << expected[2] << " "
+                   << expected[3] << " "
+                   << "actual : "
+                   << arg[0] << " "
+                   << arg[1] << " "
+                   << arg[2] << " "
+                   << arg[3];
+  return false;
 }
 
 // Use as:
@@ -259,6 +288,23 @@
   EXPECT_THAT(axis_angle, IsNearAngleAxis(expected));
 }
 
+TEST(Rotation, QuaternionToAngleAxisAngleIsLessThanPi) {
+  double quaternion[4];
+  double angle_axis[3];
+
+  const double half_theta = 0.75 * kPi;
+
+  quaternion[0] = cos(half_theta);
+  quaternion[1] = 1.0 * sin(half_theta);
+  quaternion[2] = 0.0;
+  quaternion[3] = 0.0;
+  QuaternionToAngleAxis(quaternion, angle_axis);
+  const double angle = sqrt(angle_axis[0] * angle_axis[0] +
+                            angle_axis[1] * angle_axis[1] +
+                            angle_axis[2] * angle_axis[2]);
+  EXPECT_LE(angle, kPi);
+}
+
 static const int kNumTrials = 10000;
 
 // Takes a bunch of random axis/angle values, converts them to quaternions,