Skip to content

Commit

Permalink
Correct CAD shape finder, Axis lengths are actually un-normalized ratios
Browse files Browse the repository at this point in the history
  • Loading branch information
@nyoungbq
nyoungbq committed Jan 2, 2025
1 parent 9c48092 commit 7479296
Showing 1 changed file with 55 additions and 60 deletions.
115 changes: 55 additions & 60 deletions ...entationAnalysis/src/OrientationAnalysis/Filters/Algorithms/ComputeTriangleGeomShapes.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -241,12 +241,7 @@ Result<> ComputeTriangleGeomShapes::operator()()
* The main goal is to derive the eigenvalues from the moment of inertia tensor therein finding the eigenvectors,
* which are the angular velocity vectors.
*/

// !!! DO NOT REMOVE ZEROING !!! It is integral to numerical stability in the following calculations
// Zero out small numbers in Cinertia: https://stackoverflow.com/a/54505281
Cinertia = (0.000000001 < Cinertia.array().abs()).select(Cinertia, 0.0);

Eigen::EigenSolver<Matrix3x3> eigenSolver(Cinertia); // pass in HQR Matrix for implicit calculation
Eigen::EigenSolver<Matrix3x3> eigenSolver(Cinertia);

// The primary axis is the largest eigenvalue
Eigen::EigenSolver<Matrix3x3>::EigenvalueType eigenvalues = eigenSolver.eigenvalues();
Expand All @@ -257,63 +252,26 @@ Result<> ComputeTriangleGeomShapes::operator()()
/**
* Following section for debugging
*/
// std::cout << "Eigenvalues:\n" << eigenvalues << std::endl;
// std::cout << "\n Eigenvectors:\n" << eigenvectors << std::endl;
// std::cout << "Eigenvalues:\n" << eigenvalues << std::endl;
// std::cout << "\n Eigenvectors:\n" << eigenvectors << std::endl;
//
// constexpr char k_BaselineAxisLabel = 'x'; // x
// char axisLabel = 'x';
// double primaryAxis = eigenvalues[0].real();
// for(usize i = 1; i < eigenvalues.size(); i++)
// {
// if(primaryAxis < eigenvalues[i].real())
// {
// axisLabel = k_BaselineAxisLabel + static_cast<char>(i);
// primaryAxis = eigenvalues[i].real();
// }
// }
// std::cout << "\nPrimary Axis: " << axisLabel << " | Associated Eigenvalue: " << primaryAxis << std::endl;
// constexpr char k_BaselineAxisLabel = 'x'; // x
// char axisLabel = 'x';
// double primaryAxis = eigenvalues[0].real();
// for(usize i = 1; i < eigenvalues.size(); i++)
// {
// if(primaryAxis < eigenvalues[i].real())
// {
// axisLabel = k_BaselineAxisLabel + static_cast<char>(i);
// primaryAxis = eigenvalues[i].real();
// }
// }
// std::cout << "\nPrimary Axis: " << axisLabel << " | Associated Eigenvalue: " << primaryAxis << std::endl;

// Presort eigen ordering for following sections
// Returns the argument order sorted high to low
std::array<size_t, 3> idxs = ::TripletSort(eigenvalues[0].real(), eigenvalues[1].real(), eigenvalues[2].real(), false);

/**
* The following section finds axes
*/
{
if(m_ShouldCancel)
{
return {};
}
// Formula: I = (15.0 * eigenvalue) / (4.0 * Pi);
// in the below implementation the original divisor has been put under one to avoid repeated division during execution
double I1 = (15.0 * eigenvalues[idxs[0]].real()) * k_Multiplier;
double I2 = (15.0 * eigenvalues[idxs[1]].real()) * k_Multiplier;
double I3 = (15.0 * eigenvalues[idxs[2]].real()) * k_Multiplier;

// Adjust to ABC of ellipsoid volume
double aRatio = (I1 + I2 - I3) * 0.5;
double bRatio = (I1 + I3 - I2) * 0.5;
double cRatio = (I2 + I3 - I1) * 0.5;
double a = (aRatio * aRatio * aRatio * aRatio) / (bRatio * cRatio);
a = std::pow(a, 0.1);
double b = bRatio / aRatio;
b = std::sqrt(b) * a;
double c = aRatio / (a * a * a * b);

axisLengths[3 * featureId] = static_cast<float32>(a / k_ScaleFactor);
axisLengths[3 * featureId + 1] = static_cast<float32>(b / k_ScaleFactor);
axisLengths[3 * featureId + 2] = static_cast<float32>(c / k_ScaleFactor);
auto bOverA = static_cast<float32>(b / a);
auto cOverA = static_cast<float32>(c / a);
if(aRatio == 0 || bRatio == 0 || cRatio == 0)
{
bOverA = 0.0f, cOverA = 0.0f;
}
aspectRatios[2 * featureId] = bOverA;
aspectRatios[2 * featureId + 1] = cOverA;
}

/**
* The following section calculates the axis eulers
*/
Expand All @@ -335,9 +293,9 @@ Result<> ComputeTriangleGeomShapes::operator()()
// insert principal unit vectors into rotation matrix representing Feature reference frame within the sample reference frame
//(Note that the 3 direction is actually the long axis and the 1 direction is actually the short axis)
// clang-format off
double g[3][3] = {{col1(0).real(), col1(1).real(), col1(2).real()},
{col2(0).real(), col2(1).real(), col2(2).real()},
{col3(0).real(), col3(1).real(), col3(2).real()}};
double g[3][3] = {{col1(0).real(), col1(1).real(), col1(2).real()},
{col2(0).real(), col2(1).real(), col2(2).real()},
{col3(0).real(), col3(1).real(), col3(2).real()}};
// clang-format on

// check for right-handedness
Expand All @@ -355,6 +313,43 @@ Result<> ComputeTriangleGeomShapes::operator()()
axisEulerAngles[3 * featureId + 1] = euler[1];
axisEulerAngles[3 * featureId + 2] = euler[2];
}

/**
* The following section finds axes
*/
{
if(m_ShouldCancel)
{
return {};
}
// Formula: I = (15.0 * eigenvalue) / (4.0 * Pi);
// in the below implementation the original divisor has been put under one to avoid repeated division during execution
double I1 = (15.0 * eigenvalues[idxs[0]].real()) * k_Multiplier;
double I2 = (15.0 * eigenvalues[idxs[1]].real()) * k_Multiplier;
double I3 = (15.0 * eigenvalues[idxs[2]].real()) * k_Multiplier;

// Adjust to ABC of ellipsoid volume
double aRatio = (I1 + I2 - I3) * 0.5;
double bRatio = (I1 + I3 - I2) * 0.5;
double cRatio = (I2 + I3 - I1) * 0.5;
double a = (aRatio * aRatio * aRatio * aRatio) / (bRatio * cRatio);
a = std::pow(a, 0.1);
double b = bRatio / aRatio;
b = std::sqrt(b) * a;
double c = aRatio / (a * a * a * b);

axisLengths[3 * featureId] = static_cast<float32>(a / k_ScaleFactor);
axisLengths[3 * featureId + 1] = static_cast<float32>(b / k_ScaleFactor);
axisLengths[3 * featureId + 2] = static_cast<float32>(c / k_ScaleFactor);
auto bOverA = static_cast<float32>(b / a);
auto cOverA = static_cast<float32>(c / a);
if(aRatio == 0 || bRatio == 0 || cRatio == 0)
{
bOverA = 0.0f, cOverA = 0.0f;
}
aspectRatios[2 * featureId] = bOverA;
aspectRatios[2 * featureId + 1] = cOverA;
}
}

return {};
Expand Down

0 comments on commit 7479296

Please sign in to comment.