diff --git a/glm/gtx/matrix_interpolation.inl b/glm/gtx/matrix_interpolation.inl
index 8e99ac4d..3ba49df0 100644
--- a/glm/gtx/matrix_interpolation.inl
+++ b/glm/gtx/matrix_interpolation.inl
@@ -1,6 +1,8 @@
 /// @ref gtx_matrix_interpolation
 /// @file glm/gtx/matrix_interpolation.hpp
 
+#include "../gtc/constants.hpp"
+
 namespace glm
 {
 	template<typename T, precision P>
@@ -72,7 +74,11 @@ namespace glm
 		T s = sqrt((mat[2][1] - mat[1][2]) * (mat[2][1] - mat[1][2]) + (mat[2][0] - mat[0][2]) * (mat[2][0] - mat[0][2]) + (mat[1][0] - mat[0][1]) * (mat[1][0] - mat[0][1]));
 		if (glm::abs(s) < T(0.001))
 			s = (T)1.0;
-		angle = acos((mat[0][0] + mat[1][1] + mat[2][2] - (T)1.0) * (T)0.5);
+		T const angleCos = (mat[0][0] + mat[1][1] + mat[2][2] - (T)1.0) * (T)0.5;
+		if (angleCos - static_cast<T>(1)) < epsilon)
+			angle = pi<T>() * static_cast<T>(0.25);
+		else
+			angle = acos(angleCos);
 		axis.x = (mat[1][2] - mat[2][1]) / s;
 		axis.y = (mat[2][0] - mat[0][2]) / s;
 		axis.z = (mat[0][1] - mat[1][0]) / s;
diff --git a/readme.md b/readme.md
index ce9c4be9..f50970e5 100644
--- a/readme.md
+++ b/readme.md
@@ -82,6 +82,7 @@ glm::mat4 camera(float Translate, glm::vec2 const & Rotate)
 - Fixed usused variable warning in GTX_spline #618
 - Fixed references to GLM_FORCE_RADIANS which was removed #642
 - Fixed glm::fastInverseSqrt to use fast inverse square #640
+- Fixed axisAngle NaN #638
 
 #### Deprecation:
 - Requires Visual Studio 2013, GCC 4.7, Clang 3.4, Cuda 7, ICC 2013 or a C++11 compiler
diff --git a/test/gtx/gtx_matrix_interpolation.cpp b/test/gtx/gtx_matrix_interpolation.cpp
index 6498339b..e35fa100 100644
--- a/test/gtx/gtx_matrix_interpolation.cpp
+++ b/test/gtx/gtx_matrix_interpolation.cpp
@@ -1,9 +1,45 @@
 #define GLM_ENABLE_EXPERIMENTAL
+#include <glm/gtc/quaternion.hpp>
 #include <glm/gtx/matrix_interpolation.hpp>
 
+#include <iostream>
+
+int test_axisAngle()
+{
+	int Error = 0;
+
+	float p = 0.171654f;
+	glm::mat4 m1(-0.9946f, 0.0f, -0.104531f, 0.0f,
+		0.0f, 1.0f, 0.0f, 0.0f,
+		0.104531f, 0.0f, -0.9946f, 0.0f,
+		0.0f, 0.0f, 0.0f, 1.0f);
+	glm::mat4 m2(-0.992624f, 0.0f, -0.121874f, 0.0f,
+		0.0f, 1.0f, 0.0f, 0.0f,
+		0.121874f, 0.0f, -0.992624f, 0.0f,
+		0.0f, 0.0f, 0.0f, 1.0f);
+
+	glm::mat4 const m1rot = glm::extractMatrixRotation(m1);
+	glm::mat4 const dltRotation = m2 * glm::transpose(m1rot);
+
+	glm::vec3 dltAxis(0.0f);
+	float dltAngle = 0.0f;
+	glm::axisAngle(dltRotation, dltAxis, dltAngle);
+
+	std::cout << "dltAngle: (" << dltAxis.x << ", " << dltAxis.y << ", " << dltAxis.z << "), dltAngle: " << dltAngle << std::endl;
+
+	glm::fquat q = glm::quat_cast(dltRotation);
+	std::cout << "q: (" << q.x << ", " << q.y << ", " << q.z << ", " << q.w << ")" << std::endl;
+	float yaw = glm::yaw(q);
+	std::cout << "Yaw: " << yaw << std::endl;
+
+	return Error;
+}
+
 int main()
 {
-	int Error(0);
+	int Error = 0;
+
+	Error += test_axisAngle();
 
 	return Error;
 }