diff --git a/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/ComputeTriangleGeomShapes.cpp b/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/ComputeTriangleGeomShapes.cpp index 9afe29b2ed..0e28f6ad3e 100644 --- a/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/ComputeTriangleGeomShapes.cpp +++ b/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/ComputeTriangleGeomShapes.cpp @@ -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 eigenSolver(Cinertia); // pass in HQR Matrix for implicit calculation + Eigen::EigenSolver eigenSolver(Cinertia); // The primary axis is the largest eigenvalue Eigen::EigenSolver::EigenvalueType eigenvalues = eigenSolver.eigenvalues(); @@ -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(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(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 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(a / k_ScaleFactor); - axisLengths[3 * featureId + 1] = static_cast(b / k_ScaleFactor); - axisLengths[3 * featureId + 2] = static_cast(c / k_ScaleFactor); - auto bOverA = static_cast(b / a); - auto cOverA = static_cast(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 */ @@ -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 @@ -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(a / k_ScaleFactor); + axisLengths[3 * featureId + 1] = static_cast(b / k_ScaleFactor); + axisLengths[3 * featureId + 2] = static_cast(c / k_ScaleFactor); + auto bOverA = static_cast(b / a); + auto cOverA = static_cast(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 {};