diff --git a/include/geometrycentral/surface/barycentric_coordinate_helpers.h b/include/geometrycentral/surface/barycentric_coordinate_helpers.h
index bac14dd5..35435581 100644
--- a/include/geometrycentral/surface/barycentric_coordinate_helpers.h
+++ b/include/geometrycentral/surface/barycentric_coordinate_helpers.h
@@ -23,7 +23,6 @@ double displacementLength2(Vector3 displacement, Vector3 triangleLengths);
 double displacementLength(Vector3 displacement, Vector3 triangleLengths);
 
 // Convert between cartesian and barycentric coordinates
-Vector3 cartesianVectorToBarycentric(const std::array<Vector2, 3>& vertCoords, Vector2 faceVec);
 Vector2 barycentricDisplacementToCartesian(const std::array<Vector2, 3>& vertCoords, Vector3 baryVec);
 
 // Normalize to sum to 1 (does nothing else)
diff --git a/include/geometrycentral/surface/barycentric_coordinate_helpers.ipp b/include/geometrycentral/surface/barycentric_coordinate_helpers.ipp
index 5e502356..44f645fb 100644
--- a/include/geometrycentral/surface/barycentric_coordinate_helpers.ipp
+++ b/include/geometrycentral/surface/barycentric_coordinate_helpers.ipp
@@ -25,26 +25,6 @@ inline Vector2 barycentricDisplacementToCartesian(const std::array<Vector2, 3>&
   return vertCoords[0] * baryVec.x + vertCoords[1] * baryVec.y + vertCoords[2] * baryVec.z;
 }
 
-inline Vector3 cartesianVectorToBarycentric(const std::array<Vector2, 3>& vertCoords, Vector2 faceVec, bool verbose) {
-
-
-  // Build matrix for linear transform problem
-  // (last constraint comes from chosing the displacement vector with sum = 0)
-  Eigen::Matrix3d A;
-  Eigen::Vector3d rhs;
-  const std::array<Vector2, 3>& c = vertCoords; // short name
-  A << c[0].x, c[1].x, c[2].x, c[0].y, c[1].y, c[2].y, 1., 1., 1.;
-  rhs << faceVec.x, faceVec.y, 0.;
-
-  // Solve
-  Eigen::Vector3d result = A.colPivHouseholderQr().solve(rhs);
-  Vector3 resultBary{result(0), result(1), result(2)};
-
-  resultBary = normalizeBarycentricDisplacement(resultBary);
-
-  return resultBary;
-}
-
 // Allows not-normalized input
 inline Vector3 normalizeBarycentric(Vector3 baryCoords) {
   double s = sum(baryCoords);
diff --git a/src/surface/trace_geodesic.cpp b/src/surface/trace_geodesic.cpp
index 7a9f55ad..754fe6d7 100644
--- a/src/surface/trace_geodesic.cpp
+++ b/src/surface/trace_geodesic.cpp
@@ -52,35 +52,28 @@ inline Vector2 baryCoordsToFaceCoords(const std::array<Vector2, 3>& vertCoords,
   return vertCoords[0] * baryCoord.x + vertCoords[1] * baryCoord.y + vertCoords[2] * baryCoord.z;
 }
 
+// Specialized barycentric coordinate conversion by solving the system of equations directly
+// via Cramer's rule.  Significantly faster than using a generic solver from Eigen.
+//
+// Adapted from: https://gamedev.stackexchange.com/a/23745
+//
+// Which itself is from: http://realtimecollisiondetection.net/ by Christer Ericson
 inline Vector3 cartesianVectorToBarycentric(const std::array<Vector2, 3>& vertCoords, Vector2 faceVec) {
-
-
-  // Build matrix for linear transform problem
-  // (last constraint comes from chosing the displacement vector with sum = 0)
-  Eigen::Matrix3d A;
-  Eigen::Vector3d rhs;
-  const std::array<Vector2, 3>& c = vertCoords; // short name
-  A << c[0].x, c[1].x, c[2].x, c[0].y, c[1].y, c[2].y, 1., 1., 1.;
-  rhs << faceVec.x, faceVec.y, 0.;
-
-  // Solve
-  Eigen::Vector3d result = A.colPivHouseholderQr().solve(rhs);
-  Vector3 resultBary{result(0), result(1), result(2)};
-
-  resultBary = normalizeBarycentricDisplacement(resultBary);
-
-  if (TRACE_PRINT) {
-    cout << "       cartesianVectorToBarycentric() " << endl;
-    cout << "         input = " << faceVec << endl;
-    cout << "         positions = " << vertCoords[0] << " " << vertCoords[1] << " " << vertCoords[2] << endl;
-    cout << "         transform result = " << resultBary << endl;
-    cout << "         transform back = " << barycentricDisplacementToCartesian(vertCoords, resultBary) << endl;
-    cout << " A = " << endl << A << endl;
-    cout << " rhs = " << endl << rhs << endl;
-    cout << " Ax = " << endl << A * result << endl;
-  }
-
-  return resultBary;
+  Vector2 v0 = vertCoords[1] - vertCoords[0];
+  Vector2 v1 = vertCoords[2] - vertCoords[0];
+  Vector2 v2 = faceVec - vertCoords[0];
+
+  double d00 = dot(v0, v0);
+  double d01 = dot(v0, v1);
+  double d11 = dot(v1, v1);
+  double d20 = dot(v2, v0);
+  double d21 = dot(v2, v1);
+  double denom = d00 * d11 - d01 * d01;
+
+  double v = (d11 * d20 - d01 * d21) / denom;
+  double w = (d00 * d21 - d01 * d20) / denom;
+  double u = 0. - v - w;
+  return Vector3{u, v, w};
 }
 
 // Converts tCross from halfedge to edge coordinates, handling sign conventions