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;