From 5603d1a572fb99b3ae47f5074cae2c51574a8057 Mon Sep 17 00:00:00 2001 From: Casey Primozic Date: Wed, 11 Oct 2023 18:03:20 -0700 Subject: [PATCH 1/3] Optimize `cartesianVectorToBarycentric` * Makes geodesic path tracing ~100x faster --- src/surface/trace_geodesic.cpp | 28 +++++++++++++++++++++++++--- 1 file changed, 25 insertions(+), 3 deletions(-) diff --git a/src/surface/trace_geodesic.cpp b/src/surface/trace_geodesic.cpp index 7a9f55ad..2a5d0295 100644 --- a/src/surface/trace_geodesic.cpp +++ b/src/surface/trace_geodesic.cpp @@ -52,9 +52,7 @@ inline Vector2 baryCoordsToFaceCoords(const std::array& vertCoords, return vertCoords[0] * baryCoord.x + vertCoords[1] * baryCoord.y + vertCoords[2] * baryCoord.z; } -inline Vector3 cartesianVectorToBarycentric(const std::array& vertCoords, Vector2 faceVec) { - - +inline Vector3 cartesianVectorToBarycentricOLD(const std::array& vertCoords, Vector2 faceVec) { // Build matrix for linear transform problem // (last constraint comes from chosing the displacement vector with sum = 0) Eigen::Matrix3d A; @@ -83,6 +81,30 @@ inline Vector3 cartesianVectorToBarycentric(const std::array& vertCo return resultBary; } +// Version of cartesianVectorToBarycentric which is immensely faster thanks to avoiding the eigen +// overhead or whatever else is going on with their systems of equations solver. +// +// From: https://gamedev.stackexchange.com/a/23745 +// +// Which is from: http://realtimecollisiondetection.net/ by Christer Ericson +inline Vector3 cartesianVectorToBarycentric(const std::array& vertCoords, Vector2 faceVec) { + 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 inline double convertTToEdge(Halfedge he, double tCross) { if (he == he.edge().halfedge()) return tCross; From 49c5e22679906cb4ce619a38c3937c21cd3f67a5 Mon Sep 17 00:00:00 2001 From: Casey Primozic Date: Tue, 17 Oct 2023 15:35:30 -0700 Subject: [PATCH 2/3] Clean up barycentric conversion code --- .../barycentric_coordinate_helpers.ipp | 20 ---------- src/surface/trace_geodesic.cpp | 37 ++----------------- 2 files changed, 4 insertions(+), 53 deletions(-) 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& return vertCoords[0] * baryVec.x + vertCoords[1] * baryVec.y + vertCoords[2] * baryVec.z; } -inline Vector3 cartesianVectorToBarycentric(const std::array& 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& 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 2a5d0295..754fe6d7 100644 --- a/src/surface/trace_geodesic.cpp +++ b/src/surface/trace_geodesic.cpp @@ -52,41 +52,12 @@ inline Vector2 baryCoordsToFaceCoords(const std::array& vertCoords, return vertCoords[0] * baryCoord.x + vertCoords[1] * baryCoord.y + vertCoords[2] * baryCoord.z; } -inline Vector3 cartesianVectorToBarycentricOLD(const std::array& 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& 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; -} - -// Version of cartesianVectorToBarycentric which is immensely faster thanks to avoiding the eigen -// overhead or whatever else is going on with their systems of equations solver. +// Specialized barycentric coordinate conversion by solving the system of equations directly +// via Cramer's rule. Significantly faster than using a generic solver from Eigen. // -// From: https://gamedev.stackexchange.com/a/23745 +// Adapted from: https://gamedev.stackexchange.com/a/23745 // -// Which is from: http://realtimecollisiondetection.net/ by Christer Ericson +// Which itself is from: http://realtimecollisiondetection.net/ by Christer Ericson inline Vector3 cartesianVectorToBarycentric(const std::array& vertCoords, Vector2 faceVec) { Vector2 v0 = vertCoords[1] - vertCoords[0]; Vector2 v1 = vertCoords[2] - vertCoords[0]; From 3dcfaad022589fc4b5104b5a4702ca6f4d3a3673 Mon Sep 17 00:00:00 2001 From: Casey Primozic Date: Tue, 17 Oct 2023 15:38:05 -0700 Subject: [PATCH 3/3] Remove unused declaration --- include/geometrycentral/surface/barycentric_coordinate_helpers.h | 1 - 1 file changed, 1 deletion(-) 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& vertCoords, Vector2 faceVec); Vector2 barycentricDisplacementToCartesian(const std::array& vertCoords, Vector3 baryVec); // Normalize to sum to 1 (does nothing else)