diff --git a/docs/source/version_history.rst b/docs/source/version_history.rst
index 7aedf5f..cad4a03 100644
--- a/docs/source/version_history.rst
+++ b/docs/source/version_history.rst
@@ -77,6 +77,8 @@
 
 Bug Fixes
 ---------
+#. Fix RotationMatrixToAngleAxis when the angle of rotation is near
+   PI. (Tobias Strauss)
 #. Sometimes gradient norm based convergence would miss a step with a
    substantial solution quality improvement. (Rodney Hoskinson)
 #. Ignore warnings from within Eigen/SparseQR (3.2.2).
diff --git a/include/ceres/rotation.h b/include/ceres/rotation.h
index e3dbfe8..9c08c9d 100644
--- a/include/ceres/rotation.h
+++ b/include/ceres/rotation.h
@@ -47,6 +47,7 @@
 
 #include <algorithm>
 #include <cmath>
+#include "Eigen/Geometry"
 #include "glog/logging.h"
 
 namespace ceres {
@@ -361,27 +362,22 @@
   }
 
   // Case 3: theta ~ pi, this is the hard case. Since theta is large,
-  // and sin(theta) is small. Dividing by theta by sin(theta) will
-  // either give an overflow or worse still numerically meaningless
-  // results. Thus we use an alternate more complicated formula
-  // here.
-
-  // Since cos(theta) is negative, division by (1-cos(theta)) cannot
-  // overflow.
-  const T inv_one_minus_costheta = kOne / (kOne - costheta);
-
-  // We now compute the absolute value of coordinates of the axis
-  // vector using the diagonal entries of R. To resolve the sign of
-  // these entries, we compare the sign of angle_axis[i]*sin(theta)
-  // with the sign of sin(theta). If they are the same, then
-  // angle_axis[i] should be positive, otherwise negative.
+  // and sin(theta) is small. Dividing theta by sin(theta) will either
+  // give an overflow or worse still numerically meaningless
+  // results. Thus we use an alternate more complicated and expensive
+  // formula implemented by Eigen.
+  Eigen::Matrix<T, 3, 3> rotation_matrix;
   for (int i = 0; i < 3; ++i) {
-    angle_axis[i] = theta * sqrt((R(i, i) - costheta) * inv_one_minus_costheta);
-    if (((sintheta < 0.0) && (angle_axis[i] > 0.0)) ||
-        ((sintheta > 0.0) && (angle_axis[i] < 0.0))) {
-      angle_axis[i] = -angle_axis[i];
+    for (int j = 0; j < 3; ++j) {
+      rotation_matrix(i,j) = R(i,j);
     }
   }
+
+  Eigen::AngleAxis<T> aa;
+  aa.fromRotationMatrix(rotation_matrix);
+  angle_axis[0] = aa.angle() * aa.axis()[0];
+  angle_axis[1] = aa.angle() * aa.axis()[1];
+  angle_axis[2] = aa.angle() * aa.axis()[2];
 }
 
 template <typename T>
diff --git a/internal/ceres/rotation_test.cc b/internal/ceres/rotation_test.cc
index fab0a7a..84e07f7 100644
--- a/internal/ceres/rotation_test.cc
+++ b/internal/ceres/rotation_test.cc
@@ -127,8 +127,8 @@
 }
 
 // Use as:
-// double expected_axis_angle[4];
-// double actual_axis_angle[4];
+// double expected_axis_angle[3];
+// double actual_axis_angle[3];
 // EXPECT_THAT(actual_axis_angle, IsNearAngleAxis(expected_axis_angle));
 MATCHER_P(IsNearAngleAxis, expected, "") {
   if (arg == NULL) {
@@ -136,14 +136,41 @@
     return false;
   }
 
-  for (int i = 0; i < 3; i++) {
-    if (fabs(arg[i] - expected[i]) > kTolerance) {
-      *result_listener << "component " << i << " should be " << expected[i];
-      return false;
-    }
+  // The following reads a bit complicated, here is why:
+  //
+  // 1. We are computing the relative difference between arg and
+  // expected vectors and checking if it is smaller than kTolerance.
+  //
+  // |arg - expected| < kTolerance * |expected|
+  //
+  // 2. We are doing this comparison modulo Pi, i.e., we make sure
+  // that both angle axis vectors have norms less than Pi.
+  //
+  // Evaluating these two things without performing division leads to
+  // the following mildly unreadable expressions.
+
+  Eigen::Vector3d a(arg[0], arg[1], arg[2]);
+  Eigen::Vector3d e(expected[0], expected[1], expected[2]);
+  const double a_norm = a.norm();
+  const double e_norm = e.norm();
+  const double normalized_a_norm = fmod(a_norm, kPi);
+  const double normalized_e_norm = fmod(e_norm, kPi);
+  const double delta_norm =
+      (a * normalized_a_norm * e_norm  - e * normalized_e_norm * a_norm).norm();
+
+  if (delta_norm <= kLooseTolerance * normalized_e_norm * a_norm * e_norm) {
+    return true;
   }
 
-  return true;
+  *result_listener << " arg:"
+                   << " " << arg[0]
+                   << " " << arg[1]
+                   << " " << arg[2]
+                   << " was expected to be:"
+                   << " " << expected[0]
+                   << " " << expected[1]
+                   << " " << expected[2];
+  return false;
 }
 
 // Use as:
@@ -440,10 +467,7 @@
     }
     AngleAxisToRotationMatrix(in_axis_angle, matrix);
     RotationMatrixToAngleAxis(matrix, out_axis_angle);
-
-    for (int i = 0; i < 3; ++i) {
-      EXPECT_NEAR(out_axis_angle[i], in_axis_angle[i], kLooseTolerance);
-    }
+    EXPECT_THAT(in_axis_angle, IsNearAngleAxis(out_axis_angle));
   }
 }
 
@@ -480,7 +504,7 @@
     }
   }
   LOG(INFO) << "Rotation:";
-  LOG(INFO) << "EXPECTED      |        ACTUAL";
+  LOG(INFO) << "EXPECTED        |        ACTUAL";
   for (int i = 0; i < 3; ++i) {
     string line;
     for (int j = 0; j < 3; ++j) {
@@ -1056,5 +1080,38 @@
   }
 }
 
+TEST(RotationMatrixToAngleAxis, NearPiExampleOneFromTobiasStrauss) {
+  // Example from Tobias Strauss
+  const double rotation_matrix[] = {
+    -0.999807135425239,    -0.0128154391194470,   -0.0148814136745799,
+    -0.0128154391194470,   -0.148441438622958,     0.988838158557669,
+    -0.0148814136745799,    0.988838158557669,     0.148248574048196
+  };
+
+  double angle_axis[3];
+  RotationMatrixToAngleAxis(RowMajorAdapter3x3(rotation_matrix), angle_axis);
+  double round_trip[9];
+  AngleAxisToRotationMatrix(angle_axis, RowMajorAdapter3x3(round_trip));
+  EXPECT_THAT(rotation_matrix, IsNear3x3Matrix(round_trip));
+}
+
+TEST(RotationMatrixToAngleAxis, AtPi) {
+  const double theta = 0.0;
+  const double phi = M_PI / 180.0;
+  const double angle = M_PI;
+
+  double angle_axis[3];
+  angle_axis[0] = angle * sin(phi) * cos(theta);
+  angle_axis[1] = angle * sin(phi) * sin(theta);
+  angle_axis[2] = angle * cos(phi);
+
+  double rotation_matrix[9];
+  AngleAxisToRotationMatrix(angle_axis, rotation_matrix);
+
+  double angle_axis_round_trip[3];
+  RotationMatrixToAngleAxis(rotation_matrix, angle_axis_round_trip);
+  EXPECT_THAT(angle_axis_round_trip, IsNearAngleAxis(angle_axis));
+}
+
 }  // namespace internal
 }  // namespace ceres