Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Optimize barycentric coordinate computation #164

Open
wants to merge 3 commits into
base: master
Choose a base branch
from
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

why delete this?

Vector2 barycentricDisplacementToCartesian(const std::array<Vector2, 3>& vertCoords, Vector3 baryVec);

// Normalize to sum to 1 (does nothing else)
Expand Down
20 changes: 0 additions & 20 deletions include/geometrycentral/surface/barycentric_coordinate_helpers.ipp
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Expand Down
49 changes: 21 additions & 28 deletions src/surface/trace_geodesic.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down